gen_nc_grid.ipynb gen_cesm_files.ipynb

gen_cesm_files.ipynb

gen_nc_grid.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating the NetCDF Grid file for MOM6+CESM\n",
    "\n",
    "\n",
    "This notebook describes how to generate the **netcdf grid file** for MOM6. This grid file can then be used to generate SCRIP file (for ESMF mapping tools), binary grid files (for CICE), mapping files (for CESM). See ```gen_cesm_files.ipynb``` notebook file for generating those files. First, to generate the netcdf grid file, run the Python script presented in this notebook.\n",
    "\n",
    "## 0. Prerequisities\n",
    "\n",
    "First, load the following modules:\n",
    " - Python\n",
    " - numpy\n",
    " - netCDF4\n",
    " - matplotlib\n",
    " - pyngl\n",
    " \n",
    "On cheyenne, for example:\n",
    "\n",
    "```\n",
    "module load python/2.7.13\n",
    "module load numpy\n",
    "module load netcdf4-python/1.2.7\n",
    "module load matplotlib\n",
    "module load pyngl\n",
    "```\n",
    "\n",
    "## 1. Generate the netcdf grid file\n",
    "\n",
    "MOM6 itself reads in a supergrid, a refined version of the actual horizontal grid, to generate its own data structures. This gives the flexibility of aggregating the supergrid into the B, T, and U grids (?). CESM and CIME, however, requires the actual model (ocn) grid for coupling. \n",
    "\n",
    "We first generate a netcdf file containing grid variables, e.g., coordinates, areas, masks, etc. with the normal resolution (coarser than super grid). This file will be used to generate the SCRIP file. The following Python script generates this netcdf file. NOTE: Make sure to update the input file directories and names if necessary.\n",
    "\n",
    "Import all the necessary Python modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import netCDF4\n",
    "import numpy as np\n",
    "from mpl_toolkits.basemap import pyproj\n",
    "from datetime import datetime\n",
    "import Ngl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Read in the MOM6 supergrid, which is usually four times the resolution of the actual grid (4n)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# MOM6 supergrid:\n",
    "nc_sgrd = netCDF4.Dataset('/glade/p/work/gmarques/cesm/grid/MOM6/tx0.66v1/tx0.66v1_hgrid_2018-03-27.nc')\n",
    "x = nc_sgrd.variables['x'][:]\n",
    "y = nc_sgrd.variables['y'][:]\n",
    "dx = nc_sgrd.variables['dx'][:]\n",
    "dy = nc_sgrd.variables['dy'][:]\n",
    "area = nc_sgrd.variables['area'][:]\n",
    "angle_dx = nc_sgrd.variables['angle_dx'][:]\n",
    "nc_sgrd.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in the mask and the topography file which should have the normal resolution (n). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Mask:\n",
    "nc_msk = netCDF4.Dataset('/glade/p/work/gmarques/cesm/grid/MOM6/tx0.66v1/tx0.66v1_ocean_topo_edited_180405.nc')\n",
    "tmask = nc_msk.variables['mask'][:]\n",
    "nc_msk.close()\n",
    "\n",
    "# Topography\n",
    "nc_topo = netCDF4.Dataset('/glade/p/work/gmarques/cesm/grid/MOM6/tx0.66v1/tx0.66v1_ocean_topo_edited_180405.nc')\n",
    "depth = nc_topo.variables['depth'][:]\n",
    "nc_topo.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "By copying every 2nd element, reduce the resolution of the supergrid fields from (4n) to (n)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# T point locations\n",
    "tlon = x[1::2,1::2]\n",
    "tlat = y[1::2,1::2]\n",
    "\n",
    "# U point locations\n",
    "ulon = x[1::2,::2]\n",
    "ulat = y[1::2,::2]\n",
    "\n",
    "# V point locations\n",
    "vlon = x[::2,1::2]\n",
    "vlat = y[::2,1::2]\n",
    "\n",
    "# Corner point locations\n",
    "qlon = x[::2,::2]\n",
    "qlat = y[::2,::2]\n",
    "\n",
    "# T cell area (sum of four supergrid cells)\n",
    "tarea = area[::2,::2] + area[1::2,1::2] + area[::2,1::2] + area[::2,1::2]\n",
    "\n",
    "# t-point angle\n",
    "angle = angle_dx[1::2,1::2]\n",
    "\n",
    "# x-distance between u-points, centered at t\n",
    "dxt = dx[1::2,::2] + dx[1::2,1::2]\n",
    "\n",
    "# y-distance between v-points, centered at t\n",
    "dyt = dy[::2,1::2] + dy[1::2,1::2]\n",
    "\n",
    "# x-distance between  q-points, centered at v\n",
    "dxCv = dx[2::2,::2] + dx[2::2,1::2]\n",
    "\n",
    "# y-distance between  q-points, centered at u\n",
    "dyCu = dy[::2,2::2] + dy[1::2,2::2]\n",
    "\n",
    "# x-distance between t-points, centered at u\n",
    "dxCu = dx[1::2,1::2] + np.roll(dx[1::2,1::2], -1, axis=-1)\n",
    "\n",
    "# y-distance between t-points, centered at v\n",
    "dyCv = dy[1::2,1::2] + np.roll(dy[1::2,1::2], -1, axis=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to write variables to the netcdf file to be created:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def write_nc_var(var, name, dimensions, long_name=None, units=None, coordinates=None):\n",
    "    nc.createVariable(name, 'f8', dimensions)\n",
    "    if long_name is not None:\n",
    "        nc.variables[name].long_name = long_name\n",
    "    if units is not None:\n",
    "        nc.variables[name].units = units\n",
    "    if coordinates is not None:\n",
    "        nc.variables[name].coordinates = coordinates\n",
    "    nc.variables[name][:] = var\n",
    "    print ' ... wrote variable', name\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write the netcdf file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " ... wrote variable tlon\n",
      " ... wrote variable tlat\n",
      " ... wrote variable ulon\n",
      " ... wrote variable ulat\n",
      " ... wrote variable vlon\n",
      " ... wrote variable vlat\n",
      " ... wrote variable qlon\n",
      " ... wrote variable qlat\n",
      " ... wrote variable dxt\n",
      " ... wrote variable dyt\n",
      " ... wrote variable dxCv\n",
      " ... wrote variable dyCu\n",
      " ... wrote variable dxCu\n",
      " ... wrote variable dyCv\n",
      " ... wrote variable tarea\n",
      " ... wrote variable tmask\n",
      " ... wrote variable angle\n",
      " ... wrote variable depth\n",
      " ... wrote variable ar\n",
      " ... wrote variable egs\n"
     ]
    }
   ],
   "source": [
    "filename = 'tx0.66v1_grid.nc'\n",
    "nc = netCDF4.Dataset(filename, 'w', format='NETCDF3_64BIT')\n",
    "nc.Description = 'MOM6 2/3 degree tripolar grid (ORCA type) with tropical stretching'\n",
    "# write the author, e.g.,:\n",
    "# nc.Author = 'Fred Castruccio (fredc@ucar.edu)'\n",
    "nc.Created = datetime.now().isoformat()\n",
    "nc.type = 'tx0.66v1 file'\n",
    "\n",
    "# grid aspect ratio                    \n",
    "ar = dyt / dxt\n",
    "\n",
    "# grid effective grid spacing\n",
    "# A = 4*pi*r^2 , area of sphere of radius r\n",
    "# dA = (r*cos(theta)*dlambda)*(r*dtheta), differential area on sphere                                              \n",
    "#    = r^2*domega                                                                                                  \n",
    "# domega = dA/r^2, differential solid angle  (steradians, sr)                                                      \n",
    "# 1 sr = (180./pi)^2 square degrees                                                                                \n",
    "costheta = np.cos(tlat*np.pi/180.)\n",
    "rearth = 637122000 # Earth radius in centimeter                                                                    \n",
    "domega = tarea / rearth**2\n",
    "egs  = np.sqrt(domega * (180./np.pi)**2)\n",
    "\n",
    "# create netcdf dimension\n",
    "M, L = qlon.shape\n",
    "nc.createDimension('nyp', M)\n",
    "nc.createDimension('nxp', L)\n",
    "nc.createDimension('ny', M-1)\n",
    "nc.createDimension('nx', L-1)\n",
    "\n",
    "\n",
    "# write variable\n",
    "write_nc_var(tlon, 'tlon', ('ny', 'nx'), 'array of t-grid longitudes', 'degrees_east')\n",
    "write_nc_var(tlat, 'tlat', ('ny', 'nx'), 'array of t-grid latitudes', 'degrees_north')\n",
    "write_nc_var(ulon, 'ulon', ('ny', 'nxp'), 'array of u-grid longitudes', 'degrees_east')\n",
    "write_nc_var(ulat, 'ulat', ('ny', 'nxp'), 'array of u-grid latitudes', 'degrees_north')\n",
    "write_nc_var(vlon, 'vlon', ('nyp', 'nx'), 'array of v-grid longitudes', 'degrees_east')\n",
    "write_nc_var(vlat, 'vlat', ('nyp', 'nx'), 'array of v-grid latitudes', 'degrees_north')\n",
    "write_nc_var(qlon, 'qlon', ('nyp', 'nxp'), 'array of q-grid longitudes', 'degrees_east')\n",
    "write_nc_var(qlat, 'qlat', ('nyp', 'nxp'), 'array of q-grid latitudes', 'degrees_north')\n",
    "\n",
    "write_nc_var(dxt, 'dxt', ('ny', 'nx'), 'x-distance between u-points, centered at t', 'meters')\n",
    "write_nc_var(dyt, 'dyt', ('ny', 'nx'), 'y-distance between v-points, centered at t', 'meters')\n",
    "write_nc_var(dxCv, 'dxCv', ('ny', 'nx'), 'x-distance between  q-points, centered at v', 'meters')\n",
    "write_nc_var(dyCu, 'dyCu', ('ny', 'nx'), 'y-distance between  q-points, centered at u', 'meters')\n",
    "write_nc_var(dxCu, 'dxCu', ('ny', 'nx'), 'x-distance between  t-points, centered at u', 'meters')\n",
    "write_nc_var(dyCv, 'dyCv', ('ny', 'nx'), 'y-distance between  t-points, centered at v', 'meters')\n",
    "write_nc_var(tarea, 'tarea', ('ny', 'nx'), 'area of t-cells', 'meters^2')\n",
    "write_nc_var(tmask, 'tmask', ('ny', 'nx'), 'ocean fraction at t-cell centers', 'none')\n",
    "write_nc_var(angle, 'angle', ('ny', 'nx'), 'angle grid makes with latitude line', 'degrees')\n",
    "write_nc_var(depth, 'depth', ('ny', 'nx'), 'depth at h points', 'meters')\n",
    "write_nc_var(ar, 'ar', ('ny', 'nx'), 'grid aspect ratio (dyt/dxt)', 'none')\n",
    "write_nc_var(egs, 'egs', ('ny', 'nx'), 'grid effective grid spacing', 'degrees')\n",
    "\n",
    "# close files\n",
    "nc.close()\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}