The following notes document the process used to create
mksrf_GlacierRegion_10x10min_nomask_c170616.nc

The 10' resolution was chosen somewhat arbitrarily, but considerations were: (a)
wanted to use a resolution for which we already have at least one raw data file,
so that no new mapping files would be needed, and (b) wanted to use a high
enough resolution, without going to higher resolution than is really needed for
this purpose.

Short summary:

We need 3 masks:

- Glacier region = 1: Points in the Greenland grid that are outside CISM's
  initial ice mask.

- Glacier region = 2: Points in the Greenland grid that are within CISM's
  initial ice mask.

- Glacier region = 3: Points in Antarctica.

Glacier region 3 (Antarctica) is straightforward. We generated this mask by
assuming that any grid cell south of 60 deg S is Antarctica.

The distinction between regions 1 & 2 is subtle: For some purposes, we need a
mask that encompasses any area which might ever be inside CISM's ice mask (i.e.,
any point where ice could potentially form in CISM). This can be determined by
the full Greenland grid, and is given by the union of regions 1 and 2. However,
for the sake of determining where we should convert snow capping to liquid
(which we do for glaciers outside of the major ice sheet regions), the full
Greenland grid is too broad. In particular, the full Greenland grid overlaps
with the Canadian archipelago, which is the primary region where we DO want to
convert snow capping to liquid (to avoid runaway sea ice growth). So we break
apart the full Greenland grid region into two sub-regions to handle these
different needs. We delineate these two regions based on CISM's initial ice
mask, as a convenient way to separate the island of Greenland from other points
within the CISM grid.

This is complicated further by the fact that we do not want to change answers at
this point, relative to the somewhat kludgey implementation that has been used
for the CESM2 testing. So, to be consistent with these earlier runs, we continue
to use the CISM1 5-km grid to determine the union of regions 1 & 2, but use the
CISM2 4-km grid to delineate regions 1 vs. 2.

Note that, in regridding this field, mksurfdata_map takes the maximum value from
any source point. So the ordering of regions 1 and 2 is important: If a
destination cell has any overlapping source cells within CISM's ice mask, we
want it to be classified as inside the ice mask.

Complete set of steps:

(1) Created ESMF mapping file between the 5-km Greenland grid and the 10' grid:

From cime5.3.0-alpha.27: tools/mapping/gen_mapping_files/gen_ESMF_mapping_file

Used SCRIP grid file that had already been created for the 5-km Greenland grid
(by Jeremy Fyke, using a different method than the one documented in the CISM
tools directory).

Used the following batch script:

   #!/bin/bash
   #
   #
   # Batch script to submit to create ESMF mapping file
   #
   # Set up for yellowstone
   # 
   # yellowstone-specific batch commands:
   #BSUB -P P93300601        # project number
   #BSUB -n 8               # number of processors
   #BSUB -R "span[ptile=16]" 
   #BSUB -W 1:00             # wall-clock limit
   #BSUB -q caldera            # queue
   #BSUB -o regrid.%J.out    # ouput filename
   #BSUB -e regrid.%J.err    # error filename
   #BSUB -J create_ESMF_map  # job name
   #BSUB -N                  # send email upon job completion

   #----------------------------------------------------------------------

   #----------------------------------------------------------------------
   # Set user-defined parameters here
   #----------------------------------------------------------------------

   filesrc="$CESMDATAROOT/inputdata/lnd/clm2/mappingdata/grids/SCRIPgrid_10x10min_nomask_c110228.nc"
   filedst="/glade/p/cesmdata/cseg/inputdata/glc/cism/griddata/SCRIPgrid_gland_5km_c150511.nc"

   namesrc='10x10min'
   namedst='gland5km'

   typesrc='global'
   typedst='regional'
   maptype='aave'

   #----------------------------------------------------------------------
   # Done setting user-defined parameters
   #----------------------------------------------------------------------

   #----------------------------------------------------------------------
   # Stuff done in a machine-specific way
   #----------------------------------------------------------------------

   # Determine number of processors we're running on
   host_array=($LSB_HOSTS)
   REGRID_PROC=${#host_array[@]}

   #----------------------------------------------------------------------
   # Begin general script
   #----------------------------------------------------------------------

   cmdargs="--filesrc $filesrc --filedst $filedst --namesrc $namesrc --namedst $namedst --typesrc $typesrc --typedst $typedst --maptype $maptype --batch"
   env REGRID_PROC=$REGRID_PROC ./create_ESMF_map.sh $cmdargs

(2) Created a glcmask-like file

From https://svn-ccsm-models.cgd.ucar.edu/glc/trunk_tags/cism2_1_37 : tools

with these diffs:

   Index: scrip2CLMoverlap.ncl
   ===================================================================
   --- scrip2CLMoverlap.ncl	(revision 76273)
   +++ scrip2CLMoverlap.ncl	(working copy)
   @@ -11,7 +11,7 @@
    ; must contain regrid information from a CLM grid to a GLC grid
    ; -------------------------------------------------------------
    
   -infileS = addfile("map_fv0.9x1.25_TO_gland_aave.141105.nc","r")
   +infileS = addfile("map_10x10min_TO_gland5km_aave.170615.nc","r")
    
    ; ---------------------------------------
    ; read in corresponding CLM fracdata file
   @@ -19,7 +19,8 @@
    
    ; infileF    = addfile ("$CESMDATAROOT/inputdata/lnd/clm2/griddata/fracdata_48x96_gx3v7_c090915.nc","r")
    ; infileF    = addfile ("$CESMDATAROOT/inputdata/lnd/clm2/griddata/fracdata_1.9x2.5_gx1v6_c090206.nc","r")
   -infileF    = addfile ("$CESMDATAROOT/inputdata/lnd/clm2/griddata/fracdata_0.9x1.25_gx1v6_c090317.nc","r")
   +; infileF    = addfile ("$CESMDATAROOT/inputdata/lnd/clm2/griddata/fracdata_0.9x1.25_gx1v6_c090317.nc","r")
   +infileF    = addfile("/glade/p/cesmdata/cseg/inputdata/lnd/clm2/griddata/fracdata_10min_USGS_071205.nc","r")
    
    ; ------------------------------------------------
    ; read in necessary information from the two files


Renamed the resulting file to mksrf_GlacierRegion_10x10min_nomask_c170615.nc

(3) Renamed GLCMASK to GLACIER_REGION:

ncrename -v GLCMASK,GLACIER_REGION mksrf_GlacierRegion_10x10min_nomask_c170615.nc

(4) Created a field on the 4-km grid that gives the initial icemask:

In python (using netCDF4):

# This file was a copy of the original file in inputdata
dat = Dataset('greenland_4km_epsg3413_bheatflxPos_c170501.nc', 'a')
dat.createVariable('icemask', 'f4', ('y1','x1'))
topg = dat.variables['topg'][:][0,...]
thk = dat.variables['thk'][:][0,...]
dat.variables['icemask'][:] = ((thk > 0) | (topg > 0))
dat.close()

(I confirmed that this definition of icemask exactly agrees with the initial
g2x_Sg_icemask in an IG run with CISM2 at 4km.)

(5) Created ESMF mapping file from the 4-km Greenland grid to the 10' grid:

From cime5.3.0-alpha.27: tools/mapping/gen_mapping_files/gen_ESMF_mapping_file

Used the following batch script:

   #!/bin/bash
   #
   #
   # Batch script to submit to create ESMF mapping file
   #
   # Set up for yellowstone
   # 
   # yellowstone-specific batch commands:
   #BSUB -P P93300601        # project number
   #BSUB -n 8               # number of processors
   #BSUB -R "span[ptile=16]" 
   #BSUB -W 1:00             # wall-clock limit
   #BSUB -q caldera            # queue
   #BSUB -o regrid.%J.out    # ouput filename
   #BSUB -e regrid.%J.err    # error filename
   #BSUB -J create_ESMF_map  # job name
   #BSUB -N                  # send email upon job completion

   #----------------------------------------------------------------------

   #----------------------------------------------------------------------
   # Set user-defined parameters here
   #----------------------------------------------------------------------

   filesrc="/glade/p/cesmdata/cseg/inputdata/glc/cism/griddata/SCRIPgrid_greenland_4km_epsg3413_c170414.nc"
   filedst="$CESMDATAROOT/inputdata/lnd/clm2/mappingdata/grids/SCRIPgrid_10x10min_nomask_c110228.nc"

   namesrc='gland4km'
   namedst='10x10min'

   typesrc='regional'
   typedst='global'
   maptype='aave'

   #----------------------------------------------------------------------
   # Done setting user-defined parameters
   #----------------------------------------------------------------------

   #----------------------------------------------------------------------
   # Stuff done in a machine-specific way
   #----------------------------------------------------------------------

   # Determine number of processors we're running on
   host_array=($LSB_HOSTS)
   REGRID_PROC=${#host_array[@]}

   #----------------------------------------------------------------------
   # Begin general script
   #----------------------------------------------------------------------

   cmdargs="--filesrc $filesrc --filedst $filedst --namesrc $namesrc --namedst $namedst --typesrc $typesrc --typedst $typedst --maptype $maptype --batch"
   env REGRID_PROC=$REGRID_PROC ./create_ESMF_map.sh $cmdargs

(6) Remapped icemask onto the 10' grid using cime's map_field tool

First, built the map_field tool (in cime/tools/mapping/map_field), by following
the directions there.

Then ran:

./map_field -m /glade/p/work/sacks/cime/tools/mapping/gen_mapping_files/gen_ESMF_mapping_file/map_gland4km_TO_10x10min_aave.170615.nc -if /glade/scratch/sacks/cesm_code/cism2_1_37_withMods/tools/greenland_4km_epsg3413_bheatflxPos_c170501.nc -iv icemask -of icemask_10min.nc -ov icemask

(7) Modified the GLACIER_REGION field in the file created in (2) as follows:

In python (using netCDF4):

dat = Dataset('mksrf_GlacierRegion_10x10min_nomask_c170615.nc', 'a')
dat_icemask = Dataset('icemask_10min.nc')
icemask = dat_icemask.variables['icemask'][:]
glacier_region = dat.variables['GLACIER_REGION'][:]
region2 = np.where((glacier_region == 1) & (icemask <= 0))
dat.variables['GLACIER_REGION'][region2] = 2
antarctica = np.where(dat.variables['LATIXY'][:] <= -60)
dat.variables['GLACIER_REGION'][antarctica] = 3
dat.close()

(8) Revised GLACIER_REGION field: swapped regions 1 & 2

(Next time around, I should just revise step 7 to work correctly.)

First, copied mksrf_GlacierRegion_10x10min_nomask_c170615.nc to
mksrf_GlacierRegion_10x10min_nomask_c170616.nc

Then, in python:

dat = Dataset('mksrf_GlacierRegion_10x10min_nomask_c170616.nc', 'a')
glacier_region = dat.variables['GLACIER_REGION'][:]
dat.variables['GLACIER_REGION'][np.where(glacier_region == 1)] = 2
dat.variables['GLACIER_REGION'][np.where(glacier_region == 2)] = 1
dat.close()

(9) Added / modified metadata on netcdf file: see netcdf history
global attribute for details