gen_nc_grid.ipynb gen_cesm_files.ipynb

gen_cesm_files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating CESM Grid, Domain, and Mapping files for MOM6\n",
    "\n",
    "This notebook describes how to generate SCRIP file, binary grid files (for CICE), and domain and mapping files (for CIME) to run MOM6 within CESM framework. Overview of the workflow for generating these files is as follows:\n",
    "\n",
    "1. Begin with an existing netcdf grid file.\n",
    "2. Generate the SCRIP grid file. (to be used when creating mapping files)\n",
    "3. Generate the binary grid files. (to be used by CICE)\n",
    "4. Generate the mapping and domain files. (to be used by CIME)\n",
    "5. Configure CIME\n",
    "6. Create a new case with the new resolution\n",
    "7. Prepare MOM6 input files\n",
    "8. Copy SourceMods and user_nl files\n",
    "\n",
    "Steps 1 through 5 needs to be done only once for a new grid, while steps 6, 7, and 8 need to be repeated for each new case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Begin with an existing netcdf grid file.\n",
    "\n",
    "First, a netcdf grid file containing the grid variables, e.g., coordinates, areas, masks, etc. is required. (This file is NOT the SCRIP grid file or MOM6 supergrid file). If you don't have the netcdf grid file ready, follow the steps described in the notebook \"gen_nc_grid.ipynb\" to generate one.The header of an example netcdf grid file is as follows:\n",
    "\n",
    "See [Appendix-A](#Appendix---A:-Example-netcdf-grid-file-header) for an example netcdf grid file header\n",
    "\n",
    "\n",
    "## 2. Generate the SCRIP grid file\n",
    "\n",
    "### 2.1 Global Ocean SCRIP file\n",
    "\n",
    "Given the NetCDF grid file and a rectilinear togopraphy file, e.g., ETOPO1, generate the SCRIP file. This is done by running the following NCL script, which is composed of three main parts:\n",
    "\n",
    "1. Convert the original source topography grid to a SCRIP convention file.\n",
    "2. Convert the netcdf grid file to a SCRIP convention file.\n",
    "3. Generate the weights that map from the original topography grid to the ocean grid to be generated.\n",
    "\n",
    "See NCL script: [gen_scrip.ncl](#B.1:-gen_scrip.ncl)\n",
    "\n",
    "### 2.2 Coastal SCRIP file\n",
    "\n",
    "The coastal SCRIP file is used as an input when generating the rof->ocn map. To generate this file, customize and run the following NCL script.\n",
    "\n",
    "See NCL script: [gen_coastal_scrip.ncl](#B.2:-gen_coastal_scrip.ncl)\n",
    "\n",
    "\n",
    "## 3. Generate the binary grid files\n",
    "\n",
    "CICE requires two binary grid files, horizontal grid file and topography file, which should be in the same format as the corresponding POP input files. Generate these files by running the following ncl script:\n",
    "\n",
    "**Note: Update the file paths following the ```begin``` keyword in the script accordingly.**\n",
    "\n",
    "See NCL script: [gen_cice_gridfiles.ncl](#B.3:-gen_cice_gridfiles.ncl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Generate the mapping and domain files. \n",
    "\n",
    "CESM comes with a suite of auxiliary tools including scripts to generate domain and mapping files for its components. Usage of these tools are explained in this section. \n",
    "\n",
    "Note: ```$SRCROOT``` in the following commands denotes the full pathname of CESM source root directory.\n",
    "\n",
    "### 4.1 Generate atm<->ocn mapping files\n",
    "\n",
    "- First, navigate to check_maps directory to build the tool that checks the quality of the maps after they are generated.\n",
    "```\n",
    "cd $SRCROOT/cime/tools/mapping/check_maps/\n",
    "```\n",
    "\n",
    "- Follow the instructions in ```README``` file to build the check_maps tool.\n",
    "\n",
    "- Navigate to ```gen_mapping_files``` directory and run ```gen_cesm_maps.sh``` script to generate atm<->ocn mapping files as follows.\n",
    "```\n",
    "cd $SRCROOT/cime/tools/mapping/gen_mapping_files/\n",
    "./gen_cesm_maps.sh -fatm ATM_SCRIP_FILE -natm ATM_GRID_NAME -focn OCN_SCRIP_FILE -nocn OCN_GRID_NAME\n",
    "```\n",
    "where ```ATM_SCRIP_FILE``` is the path to the atmosphere SCRIP grid file (e.g., T62 grid on cheyenne: ```/glade/p/cesmdata/cseg/mapping/grids/T62_040121.nc```), ```ATM_GRID_NAME``` is the name of the atmosphere grid, e.g., T62, ```OCN_SCRIP_FILE``` is the ocean SCRIP grid file path that you created in Step 2, and ```OCN_GRID_NAME``` is the name of the ocean grid, e.g., tx0.66v1.\n",
    "\n",
    "The tool generates three atm->ocn and two ocn->atm maps. These will be described in later sections.\n",
    "\n",
    "### 4.2 Generate rof->ocn mapping files\n",
    "\n",
    "- Navigate to runoff_to_ocn directory:\n",
    "```\n",
    "cd $SRCROOT/cime/tools/mapping/gen_mapping_files/runoff_to_ocn/\n",
    "```\n",
    "\n",
    "- Follow the instructions in ```INSTALL``` file to install the tool.\n",
    "- Create a directory where you want to generate the rof2ocn mapping file. In the following commands ```$ROF2OCN``` denotes this directory.\n",
    "- Create a namelist file inside ```$ROF2OCN``` directory. The namelist entries to be provided in this file are described in ```$SRCROOT/cime/tools/mapping/gen_mapping_files/runoff_to_ocn/README``` file. An example file (for rx1 runoff grid to tx0.66v1 ocn grid) is as follows:\n",
    "```\n",
    "&input_nml\n",
    "   gridtype     = 'obs'\n",
    "   file_roff    = '/glade/p/cesm/cseg/inputdata/lnd/dlnd7/RX1/runoff.daitren.annual.090225.nc'\n",
    "   file_ocn     = '/glade/p/work/altuntas/cesm.input/grid/scrip/tx0.66v1_SCRIP_171107.nc'\n",
    "   file_ocn_coastal_mask = '/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_171220/tx0.66v1_SCRIP_coastal_171220.nc'\n",
    "   file_nn      = 'map_rx1_to_tx0.66v1_nn_171226.nc'\n",
    "   file_smooth  = 'map_rx1_to_tx0.66v1_sm_e1000r300_171226.nc'\n",
    "   file_new     = 'map_rx1_to_tx0.66v1_nnsm_e1000r300_171226.nc'\n",
    "   title        = 'runoff map: rx1 -> rx0.66v1, nearest neighbor and smoothed '\n",
    "   eFold        = 1000000.0\n",
    "   rMax         =  300000.0\n",
    "   step1 = .true.\n",
    "   step2 = .true.\n",
    "   step3 = .true.\n",
    "/\n",
    "```\n",
    "\n",
    "- Run ```runoff_map``` tool as follows:\n",
    "```\n",
    "$SRCROOT/cime/tools/mapping/gen_mapping_files/runoff_to_ocn/runoff_map < NML_FILE_PATH\n",
    "```\n",
    "where ```NML_FILE_PATH``` is the path to the namelist file.\n",
    "The tool generates three mapping files in three steps, as described in the README file. In MOM6, we only use the map that is generated in the final step, i.e., whose path is specified in ```file_new``` entry in the namelist file. \n",
    "\n",
    "### 4.2 Generate the domain files\n",
    "\n",
    "Domain files, read by CIME, contain mask information about the grids to be mapped. These files are similarly generated using the CIME tool, called gen_domain_files. \n",
    "\n",
    "\n",
    "- First, navigate to gen_domain_files directory to build the tool:\n",
    "```\n",
    "cd $SRCROOT/cime/tools/mapping/gen_domain_files/\n",
    "```\n",
    "\n",
    "- Follow the instructions in ```INSTALL``` file to build the gen_domain_files tool.\n",
    "\n",
    "- To generate the domain files, follow the instructions in ```README``` file. The following is an example usage where the domain files (three in total) for T62_t061 resolution are generated.\n",
    "\n",
    "$SRCROOT/cime/tools/mapping/gen_domain_files/gen_domain -m /glade/p/work/altuntas/cesm.input/map/T62_tx0.66v1/map_tx0.66v1_TO_T62_aave.171017.nc -o tx0.66v1 -l T62\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Configure CIME\n",
    "\n",
    "The new grid needs to be added to CIME. To do this, the resolution alias (e.g., T62_t061), and the paths to domain and mapping files (generated in step 4) should be added to config_grids.xml file in cime/config/cesm/ directory as follows.\n",
    " \n",
    "#### Note: See the following github commit for an example of how a grid gets added to CIME: https://github.com/ESMCI/cime/commit/1879acb7b07cd13f64d237f9573a6ca8b902f997\n",
    " \n",
    "- Add the new grid alias, e.g.,\n",
    " ```\n",
    "    <model_grid alias=\"T62_t061\" not_compset=\"_CAM\">\n",
    "      <grid name=\"atm\">T62</grid>\n",
    "      <grid name=\"lnd\">T62</grid>\n",
    "      <grid name=\"ocnice\">tx0.66v1</grid>\n",
    "    </model_grid>\n",
    "```\n",
    "\n",
    "- Add the atm|lnd maks file directory, e.g.:\n",
    "```\n",
    "<file grid=\"atm|lnd\" mask=\"tx0.66v1\">/glade/p/work/altuntas/cesm.input/domain/T62_tx0.66v1/domain.lnd.T62_tx0.66v1.171017.nc</file>\n",
    "```\n",
    "\n",
    "- Add the ocn mask file directory, e.g.:\n",
    "```\n",
    "<file grid=\"ocnice\" mask=\"tx0.66v1\">/glade/p/work/altuntas/cesm.input/domain/T62_tx0.66v1/domain.ocn.T62_tx0.66v1.171017.nc</file>\n",
    "```\n",
    "\n",
    "- Define the new domain, e.g.,\n",
    "```\n",
    "<domain name=\"tx0.66v1\">\n",
    "  <nx>540</nx> <ny>458</ny>\n",
    "  <file grid=\"ocnice\">/glade/p/work/altuntas/cesm.input/domain/T62_tx0.66v1/domain.ocn.tx0.66v1.171017.nc</file>\n",
    "  <desc>tx0.66v1 is tripole v1 0.66-deg MOM6 grid:</desc>\n",
    "  <support>Experimental for MOM6 experiments</support>\n",
    "</domain>\n",
    "```\n",
    "\n",
    "- Add atm<->ocn mapping file directories:\n",
    "```\n",
    "<gridmap atm_grid=\"T62\" ocn_grid=\"tx0.66v1\">\n",
    " <map name=\"ATM2OCN_FMAPNAME\">/glade/p/work/altuntas/cesm.input/map/T62_tx0.66v1/map_T62_TO_tx0.66v1_aave.171017.nc</map>\n",
    " <map name=\"ATM2OCN_SMAPNAME\">/glade/p/work/altuntas/cesm.input/map/T62_tx0.66v1/map_T62_TO_tx0.66v1_blin.171017.nc</map>\n",
    " <map name=\"ATM2OCN_VMAPNAME\">/glade/p/work/altuntas/cesm.input/map/T62_tx0.66v1/map_T62_TO_tx0.66v1_blin.171017.nc</map>\n",
    " <map name=\"OCN2ATM_FMAPNAME\">/glade/p/work/altuntas/cesm.input/map/T62_tx0.66v1/map_tx0.66v1_TO_T62_aave.171017.nc</map>\n",
    " <map name=\"OCN2ATM_SMAPNAME\">/glade/p/work/altuntas/cesm.input/map/T62_tx0.66v1/map_tx0.66v1_TO_T62_aave.171017.nc</map>\n",
    "</gridmap\n",
    "```\n",
    "\n",
    "- Add rof->ocn mapping file directories:\n",
    "```\n",
    "<gridmap rof_grid=\"rx1\" ocn_grid=\"tx0.66v1\" >\n",
    "   <map name=\"ROF2OCN_LIQ_RMAPNAME\">/glade/p/work/altuntas/cesm.input/map/rx1_tx0.66v1/map_rx1_to_tx0.66v1_nnsm_e1000r300_171010.nc</map>\n",
    "   <map name=\"ROF2OCN_ICE_RMAPNAME\">/glade/p/work/altuntas/cesm.input/map/rx1_tx0.66v1/map_rx1_to_tx0.66v1_nnsm_e1000r300_171010.nc</map>\n",
    "</gridmap>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Note: the steps above need to be done only once for a new grid. The steps below need to be repeated for each new case.\n",
    "\n",
    "## 6. Create a new case with the new resolution\n",
    "\n",
    "- Create the new case, e.g.:\n",
    "```\n",
    "create_newcase --case gmom.e20.GMOM.T62_t061.001 --compset 2000_DATM%NYF_SLND_CICE_MOM6_DROF%NYF_SGLC_SWAV --res T62_t061 --user-compset --run-unsupported\n",
    "```\n",
    "- Run ./case.setup\n",
    "\n",
    "\n",
    "## 7. Prepare MOM6 input files\n",
    "\n",
    "Collect the following set of input files and copy them into RUNDIR directory of the case:\n",
    "- diag_table\n",
    "- input.nml\n",
    "- MOM_input\n",
    "- MOM_override\n",
    "- INPUT/ocean_hgrid.nc\n",
    "- INPUT/ocean_topog.nc\n",
    "- INPUT/ocean_vgrid.nc\n",
    "- WOA05_pottemp_salt.nc\n",
    "\n",
    "Make sure that ocean_hgrid.nc, ocean_topog.nc, and ocean_vgrid.nc files are compatible with the MOM6 supergrid file (and so all other files generated by following this notebook.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Copy SourceMods:\n",
    "\n",
    "### 8.1 CIME SourceMod:\n",
    "\n",
    "At the initialization of all case runs, CIME checks for inconsistencies between grids, mapping files, domain files, etc. When comparing ice and ocean grids, for instance, CIME checks whether the difference between the areas of overlapping ice and ocean cells are less then a predefined tolerance. For an experimental grid, this tolerance may need to be relaxed until the grid is finalized. The following SourceMod is an example of how this may be done:\n",
    "```\n",
    "/glade/p/work/altuntas/cesm.sourcemods/cesm2_0_beta06.mom6/src.drv/seq_domain_mct.F90\n",
    "```\n",
    "To run a case with T62_t061 grid, copy the file above to SourceMods/src.drv/ directory within your CASEROOT before building your case. \n",
    "\n",
    "### 8.2 CICE SourceMods:\n",
    "\n",
    "#### 8.2.1 Add the new grid:\n",
    "\n",
    "Add the new grid and the paths to its files (topography and horizontal grid) to components/cice/cime_config/namelist_definition_cice.xml file. An example SourceMod can be found in:\n",
    "```\n",
    "/glade/p/work/altuntas/cesm.sourcemods/cesm2_0_beta06.mom6/src.cice/namelist_definition_cice.xml\n",
    "```\n",
    "To run a case with T62_t061 grid, copy the file above to SourceMods/src.cice/ directory within your CASEROOT before building your case. \n",
    "\n",
    "#### 8.2.2 Fix the cell coordinates interval:\n",
    "\n",
    "CICE assumes that the coordinates of cells read from the horizontal grid file are confined to be between -180 and 180 degrees (because of the way the intrinsic function atan2 works). If this is not the case, (like tx0.66v1), copy the following SourceMod to fix this issue:\n",
    "```\n",
    "/glade/p/work/altuntas/cesm.sourcemods/cesm2_0_beta06.mom6/src.cice/ice_grid.F90\n",
    "```\n",
    "To run a case with T62_t061 grid, copy the file above to SourceMods/src.cice/ directory within your CASEROOT before building your case. \n",
    "\n",
    "### 8.2 CICE user_nl_cice changes:\n",
    "\n",
    "If there is not an initial conditions file availabe for CICE yet, paste the following line to the end of user_nl_cice file to experiment with the new grid using the default initial conditions for CICE:\n",
    "```\n",
    "ice_ic = 'UNSET'\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "---------\n",
    "## Final Steps:\n",
    "\n",
    "After copying the SourceMods and user_nl files as described, build and submit the case.\n",
    "\n",
    "----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix - A: Example NetCDF grid file header\n",
    "\n",
    "```\n",
    "netcdf tx0.66v1_grid {\n",
    "dimensions:\n",
    "\tnyp = 459 ;\n",
    "\tnxp = 541 ;\n",
    "\tny = 458 ;\n",
    "\tnx = 540 ;\n",
    "variables:\n",
    "\tdouble tlon(ny, nx) ;\n",
    "\t\ttlon:long_name = \"array of t-grid longitudes\" ;\n",
    "\t\ttlon:units = \"degrees_east\" ;\n",
    "\tdouble tlat(ny, nx) ;\n",
    "\t\ttlat:long_name = \"array of t-grid latitudes\" ;\n",
    "\t\ttlat:units = \"degrees_north\" ;\n",
    "\tdouble ulon(ny, nxp) ;\n",
    "\t\tulon:long_name = \"array of u-grid longitudes\" ;\n",
    "\t\tulon:units = \"degrees_east\" ;\n",
    "\tdouble ulat(ny, nxp) ;\n",
    "\t\tulat:long_name = \"array of u-grid latitudes\" ;\n",
    "\t\tulat:units = \"degrees_north\" ;\n",
    "\tdouble vlon(nyp, nx) ;\n",
    "\t\tvlon:long_name = \"array of v-grid longitudes\" ;\n",
    "\t\tvlon:units = \"degrees_east\" ;\n",
    "\tdouble vlat(nyp, nx) ;\n",
    "\t\tvlat:long_name = \"array of v-grid latitudes\" ;\n",
    "\t\tvlat:units = \"degrees_north\" ;\n",
    "\tdouble qlon(nyp, nxp) ;\n",
    "\t\tqlon:long_name = \"array of q-grid longitudes\" ;\n",
    "\t\tqlon:units = \"degrees_east\" ;\n",
    "\tdouble qlat(nyp, nxp) ;\n",
    "\t\tqlat:long_name = \"array of q-grid latitudes\" ;\n",
    "\t\tqlat:units = \"degrees_north\" ;\n",
    "\tdouble dxt(ny, nx) ;\n",
    "\t\tdxt:long_name = \"x-distance between u-points, centered at t\" ;\n",
    "\t\tdxt:units = \"meters\" ;\n",
    "\tdouble dyt(ny, nx) ;\n",
    "\t\tdyt:long_name = \"y-distance between v-points, centered at t\" ;\n",
    "\t\tdyt:units = \"meters\" ;\n",
    "\tdouble dxCv(ny, nx) ;\n",
    "\t\tdxCv:long_name = \"x-distance between  q-points, centered at v\" ;\n",
    "\t\tdxCv:units = \"meters\" ;\n",
    "\tdouble dyCu(ny, nx) ;\n",
    "\t\tdyCu:long_name = \"y-distance between  q-points, centered at u\" ;\n",
    "\t\tdyCu:units = \"meters\" ;\n",
    "\tdouble dxCu(ny, nx) ;\n",
    "\t\tdxCu:long_name = \"x-distance between  t-points, centered at u\" ;\n",
    "\t\tdxCu:units = \"meters\" ;\n",
    "\tdouble dyCv(ny, nx) ;\n",
    "\t\tdyCv:long_name = \"y-distance between  t-points, centered at v\" ;\n",
    "\t\tdyCv:units = \"meters\" ;\n",
    "\tdouble tarea(ny, nx) ;\n",
    "\t\ttarea:long_name = \"area of t-cells\" ;\n",
    "\t\ttarea:units = \"meters^2\" ;\n",
    "\tdouble tmask(ny, nx) ;\n",
    "\t\ttmask:long_name = \"ocean fraction at t-cell centers\" ;\n",
    "\t\ttmask:units = \"none\" ;\n",
    "\tdouble angle(ny, nx) ;\n",
    "\t\tangle:long_name = \"angle grid makes with latitude line\" ;\n",
    "\t\tangle:units = \"degrees\" ;\n",
    "\tdouble depth(ny, nx) ;\n",
    "\t\tdepth:long_name = \"depth at h points\" ;\n",
    "\t\tdepth:units = \"meters\" ;\n",
    "\tdouble ar(ny, nx) ;\n",
    "\t\tar:long_name = \"grid aspect ratio (dyt/dxt)\" ;\n",
    "\t\tar:units = \"none\" ;\n",
    "\tdouble egs(ny, nx) ;\n",
    "\t\tegs:long_name = \"grid effective grid spacing\" ;\n",
    "\t\tegs:units = \"degrees\" ;\n",
    "\n",
    "// global attributes:\n",
    "\t\t:Description = \"MOM6 2/3 degree tripolar grid (ORCA type)\" ;\n",
    "\t\t:Author = \"Fred Castruccio (fredc@ucar.edu)\" ;\n",
    "\t\t:Created = \"2017-10-26T09:28:05.530173\" ;\n",
    "\t\t:type = \"tx0.66v1 file\" ;\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Appendix - B: NCL Scripts\n",
    "\n",
    "### B.1: gen_scrip.ncl\n",
    "\n",
    "**Note: Update the file paths following the ```begin``` keyword in the script accordingly.**\n",
    "\n",
    "\n",
    "```\n",
    "; gen_scrip.ncl : \n",
    "\n",
    "load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl\"\n",
    "load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl\"\n",
    "load \"$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl\"\n",
    "load \"$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl\"\n",
    "setvalues NhlGetWorkspaceObjectId()\n",
    "        \"wsMaximumSize\" : 100000000000\n",
    "end setvalues\n",
    "\n",
    "begin\n",
    ";---Input files\n",
    "    topoFilePath  = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_180604/ETOPO1_Bed_c_gmt4.grd\"\n",
    "    ncGrdFilePath = \"/glade/p/work/altuntas/scripts/notebooks/gen_tx0.66v1_180604/tx0.66v1_grid.nc\"\n",
    "    maskFilePath  = \"/glade/p/work/gmarques/cesm/grid/MOM6/tx0.66v1/tx0.66v1_ocean_topo_edited_180405.nc\"\n",
    "\n",
    ";---Output (and input) files\n",
    "    srcGridPath = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_180604/ETOPO1_SCRIP_180604.nc\"\n",
    "    dstGridPath = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_180604/tx0.66v1_SCRIP_180604.nc\"\n",
    "    wgtFilePath = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_180604/ETOPO1_2_tx0.66v1.nc\"\n",
    "\n",
    ";----------------------------------------------------------------------\n",
    "; Convert original source ETOPO1 grid to a SCRIP convention file.\n",
    ";----------------------------------------------------------------------\n",
    "\n",
    "    src_file = addfile(topoFilePath,\"r\")\n",
    "    srclon      = src_file->x(:)\n",
    "    srclat      = src_file->y(:)\n",
    "    srctopo     = src_file->z(:,:)\n",
    "\n",
    "    Opt                = True\n",
    "    Opt@ForceOverwrite = True\n",
    "    Opt@PrintTimings   = True\n",
    "    Opt@NetCDFType     = \"netcdf4\"\n",
    "    Opt@Title          = \"ETOPO1 Grid\"\n",
    "\n",
    "    rectilinear_to_SCRIP(srcGridPath,srclat,srclon,Opt)\n",
    "\n",
    ";---Clean up\n",
    "    delete(Opt)\n",
    "\n",
    ";----------------------------------------------------------------------\n",
    "; Convert destination grid to a SCRIP convention file.\n",
    ";----------------------------------------------------------------------\n",
    "    mom6_grd_file = addfile(ncGrdFilePath,\"r\")\n",
    "    tlon     = mom6_grd_file->tlon(:,:)\n",
    "    tlat     = mom6_grd_file->tlat(:,:)\n",
    "\n",
    "    Opt                = True\n",
    "    Opt@ForceOverwrite = True\n",
    "    Opt@PrintTimings   = True\n",
    "    Opt@NetCDFType     = \"netcdf4\"\n",
    "    Opt@Title          = \"tx0.66v1 Grid\"\n",
    "\n",
    "\n",
    "    ; Corners' lat lon --------------------------------------------------\n",
    "\n",
    "    grid_dimsizes = getfiledimsizes(mom6_grd_file)\n",
    "    nyp = grid_dimsizes(0)\n",
    "    nxp = grid_dimsizes(1)\n",
    "    ny = grid_dimsizes(2)\n",
    "    nx = grid_dimsizes(3)\n",
    "    n = nx*ny\n",
    "\n",
    "    Opt@GridCornerLat = new(4*n,double)\n",
    "    Opt@GridCornerLon = new(4*n,double)\n",
    "\n",
    "    do j=0,ny-1\n",
    "      do i=0,nx-1\n",
    "        bi = (j*nx+i)*4 ; index of 1st corner in cell \n",
    "\n",
    "        Opt@GridCornerLat(bi)   = mom6_grd_file->qlat(j,i)\n",
    "        Opt@GridCornerLat(bi+1) = mom6_grd_file->qlat(j,i+1)\n",
    "        Opt@GridCornerLat(bi+2) = mom6_grd_file->qlat(j+1,i)\n",
    "        Opt@GridCornerLat(bi+3) = mom6_grd_file->qlat(j+1,i+1)\n",
    "        Opt@GridCornerLon(bi)   = mom6_grd_file->qlon(j,i)\n",
    "        Opt@GridCornerLon(bi+1) = mom6_grd_file->qlon(j,i+1)\n",
    "        Opt@GridCornerLon(bi+2) = mom6_grd_file->qlon(j+1,i)\n",
    "        Opt@GridCornerLon(bi+3) = mom6_grd_file->qlon(j+1,i+1)\n",
    "      end do\n",
    "    end do\n",
    "\n",
    "    ; Mask --------------------------------------------------\n",
    "\n",
    "    maskFile = addfile(maskFilePath,\"r\")\n",
    "    Opt@GridMask = toint(maskFile->mask)\n",
    "\n",
    "    ; Generate the script file ------------------------------\n",
    "\n",
    "    curvilinear_to_SCRIP(dstGridPath,tlat,tlon,Opt)\n",
    "\n",
    ";---Clean up\n",
    "    delete(Opt)\n",
    "\n",
    "    ; Append Area to SCRIP file------------------------------\n",
    "\n",
    "    scripFile = addfile(dstGridPath,\"w\")\n",
    "\n",
    "    grid_size = dimsizes(scripFile->grid_center_lat)\n",
    "    grid_area = new(grid_size,double)\n",
    "    grid_area!0 = \"grid_size\"\n",
    "  \n",
    "    do i=0,grid_size-1\n",
    "      temp_tlat = (/ scripFile->grid_corner_lat(i,3), \\\n",
    "                scripFile->grid_corner_lat(i,1), \\\n",
    "                scripFile->grid_corner_lat(i,0), \\\n",
    "                scripFile->grid_corner_lat(i,2)    /)\n",
    "      temp_tlon = (/ scripFile->grid_corner_lon(i,3), \\\n",
    "                scripFile->grid_corner_lon(i,1), \\\n",
    "                scripFile->grid_corner_lon(i,0), \\\n",
    "                scripFile->grid_corner_lon(i,2)    /)\n",
    "  \n",
    "      grid_area(i) = area_poly_sphere(temp_tlat, temp_tlon, 1)\n",
    "  \n",
    "    end do\n",
    "  \n",
    "    scripFile->grid_area = grid_area\n",
    "\n",
    "\n",
    ";----------------------------------------------------------------------\n",
    "; Generate the weights that take you from the ETOPO1 grid to a\n",
    "; tx0.66v1 degree grid.\n",
    ";----------------------------------------------------------------------\n",
    "    Opt                      = True\n",
    "    Opt@InterpMethod         = \"bilinear\"     ; default\n",
    "    Opt@ForceOverwrite       = True\n",
    "    Opt@PrintTimings         = True\n",
    "    Opt@NetCDFType           = \"netcdf4\"\n",
    "\n",
    "    ESMF_regrid_gen_weights(srcGridPath,dstGridPath,wgtFilePath,Opt)\n",
    "\n",
    "    delete(Opt)\n",
    "\n",
    "end\n",
    "```\n",
    "\n",
    "### B.2: gen_coastal_scrip.ncl\n",
    "\n",
    "**Note: Update the file paths and grid dimensions following the ```begin``` keyword in the script accordingly.**\n",
    "\n",
    "\n",
    "```\n",
    "; gen_coastal_scrip.ncl : \n",
    "\n",
    "begin\n",
    "\n",
    "  ; old grid file\n",
    "  scripFile = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_171220/tx0.66v1_SCRIP_171220.nc\"\n",
    "  coastalScripFile = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_171220/tx0.66v1_SCRIP_coastal_171220.nc\"\n",
    "\n",
    "  nx = 540\n",
    "  ny = 458\n",
    "\n",
    "  ; -----------------------------------------------------------------------\n",
    "\n",
    "  inputFile = addfile(scripFile,\"r\")\n",
    "  system(\"/bin/rm -f \" + coastalScripFile)\n",
    "  outputFile = addfile(coastalScripFile,\"c\")\n",
    "\n",
    "  ; copy global file attributes:\n",
    "  varatts = getvaratts(inputFile)\n",
    "  if (.not. ismissing(varatts(0))) then\n",
    "    fileattdef(outputFile,inputFile)\n",
    "  end if\n",
    "\n",
    "  ; copy all variables:\n",
    "  varNames_in = getfilevarnames(inputFile);\n",
    "  nNames = dimsizes(varNames_in)\n",
    "  do i=0,nNames-1\n",
    "    outputFile->$varNames_in(i)$ = inputFile->$varNames_in(i)$\n",
    "  end do\n",
    "\n",
    "  grid_size = dimsizes(inputFile->grid_center_lat)\n",
    "  outputFile->grid_imask = inputFile->grid_imask\n",
    "\n",
    "  do i=0,nx-1\n",
    "    do j=0,ny-1\n",
    "\n",
    "      c = inputFile->grid_imask(i+nx*j)                     ; center\n",
    "\n",
    "      if (c.eq.1) then\n",
    "        n = inputFile->grid_imask(i+nx*( mod(j+1, ny) ) )     ; north\n",
    "        s = inputFile->grid_imask(i+nx*( mod(j-1+ny, ny) ) )  ; south\n",
    "        e = inputFile->grid_imask( (mod(i-1+nx,nx) +nx*j) )   ; east\n",
    "        w = inputFile->grid_imask( (mod(i+1,nx) +nx*j) )      ; west\n",
    "\n",
    "        if (j.eq.(ny-1)) then\n",
    "           n=1\n",
    "        end if\n",
    "\n",
    "        outputFile->grid_imask(i+nx*j) = 1-(n*s*e*w)\n",
    "\n",
    "      else\n",
    "        outputFile->grid_imask(i+nx*j) = c\n",
    "      end if\n",
    "\n",
    "    end do\n",
    "  end do\n",
    "\n",
    "end\n",
    "\n",
    "```\n",
    "\n",
    "### B.3: gen_cice_gridfiles.ncl\n",
    "\n",
    "**Note: Update the file paths following the ```begin``` keyword in the script accordingly.**\n",
    "\n",
    "```\n",
    "; gen_cice_gridfiles.ncl\n",
    "\n",
    "; Input/Output File Paths:\n",
    "; ===================================================================\n",
    "\n",
    "; Input Vertical Grid File:\n",
    "in_vgridPath      = \"/glade/p/work/altuntas/mom.input/tx0.66v1/case_input_171220/INPUT/ocean_vgrid.nc\"\n",
    "\n",
    "; Input Topography File:\n",
    "in_topoPath       = \"/glade/p/work/altuntas/mom.input/tx0.66v1/case_input_171220/INPUT/ocean_topog.nc\"\n",
    "\n",
    "; Input Horizontal Grid File:\n",
    "in_superGridPath  = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_171220/tx0.66v1_grid.nc\"\n",
    "\n",
    "; Input Mask File:\n",
    "in_maskPath       = \"/glade/p/work/altuntas/mom.input/tx0.66v1/gen_grid_171220/ocean_mask_edited.nc\"\n",
    "\n",
    "; Output Horizontal Grid File (binary):\n",
    "out_hgridPath     = \"/glade/p/work/altuntas/mom.input/tx0.66v1/cice_input_171220/horiz_grid_20171220.ieeer8\"\n",
    "\n",
    "; Output Topography file (binary):\n",
    "out_topoPath      = \"/glade/p/work/altuntas/mom.input/tx0.66v1/cice_input_171220/topography_20171220.ieeei4\"\n",
    "\n",
    "; ===================================================================\n",
    "\n",
    "in_vgridFile = addfile(in_vgridPath,\"r\")\n",
    "in_topoFile = addfile(in_topoPath,\"r\")\n",
    "in_superGridFile = addfile(in_superGridPath,\"r\")\n",
    "\n",
    "sizes = dimsizes(in_topoFile->depth)\n",
    "ny = sizes(0)\n",
    "nx = sizes(1)\n",
    "\n",
    "vgridSizes = dimsizes(in_vgridFile->dz)\n",
    "ndepth = vgridSizes(0)\n",
    "depths = new(ndepth, double)\n",
    "\n",
    "depths(0) = in_vgridFile->dz(0)\n",
    "do k=1,ndepth-1\n",
    "  depths(k) = in_vgridFile->dz(k) + depths(k-1)\n",
    "end do \n",
    "\n",
    "; fill in kmt:\n",
    "kmt = new((/ny,nx/), integer)\n",
    "kmt = ndepth\n",
    "kmt!0 = \"ny\"\n",
    "kmt!1 = \"nx\"\n",
    "\n",
    "in_maskFile = addfile(in_maskPath,\"r\")\n",
    "msk = toint(in_maskFile->mask)\n",
    "\n",
    "do j=0,ny-1\n",
    "  do i=0,nx-1\n",
    "    do k=0,ndepth-1\n",
    "      if (in_topoFile->depth(j,i) .lt. depths(k))\n",
    "        kmt(j,i) = k\n",
    "        if ( msk(j,i).lt.1 .and. kmt(j,i).gt.0 ) then\n",
    "          print(\" \"+j+\" \"+i+\" \"+in_topoFile->depth(j,i)+\" \"+depths(k)+\" \"+(k-1) )\n",
    "        end if\n",
    "        break    \n",
    "      end if\n",
    "    end do    \n",
    "  end do\n",
    "end do\n",
    "\n",
    "dtr = get_d2r(\"double\")\n",
    "PI  = get_pi(\"double\")\n",
    "\n",
    "ulat = dtr * in_superGridFile->qlat(1:ny,1:nx)\n",
    "ulon = dtr * in_superGridFile->qlon(1:ny,1:nx) \n",
    "htn = (in_superGridFile->dxCv) * 100.0 ; convert to cm\n",
    "hte = (in_superGridFile->dyCu) * 100.0 ; convert to cm\n",
    "hus = (in_superGridFile->dxCu) * 100.0 ; convert to cm\n",
    "huw = (in_superGridFile->dyCv) * 100.0 ; convert to cm\n",
    "angle = dtr * (in_superGridFile->angle)\n",
    "\n",
    "print(\" \"+max(htn)+\" \"+min(htn))\n",
    "print(\" \"+max(hte)+\" \"+min(hte))\n",
    "\n",
    "; write the topography file\n",
    "setfileoption(\"bin\",\"WriteByteOrder\",\"BigEndian\")\n",
    "system(\"/bin/rm -f \" + out_topoPath)\n",
    "fbindirwrite(out_topoPath, (kmt))\n",
    "\n",
    "; write the horizontal grid file\n",
    "system(\"/bin/rm -f \" + out_hgridPath)\n",
    "fbindirwrite(out_hgridPath, (ulat))\n",
    "fbindirwrite(out_hgridPath, (ulon))\n",
    "fbindirwrite(out_hgridPath, (htn))\n",
    "fbindirwrite(out_hgridPath, (hte))\n",
    "fbindirwrite(out_hgridPath, (hus))\n",
    "fbindirwrite(out_hgridPath, (huw))\n",
    "fbindirwrite(out_hgridPath, (angle))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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
}

gen_nc_grid.ipynb