Add COMPASS step for remapping prescribed ISMF

testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_e3sm_coupling.xml

@@ -8,6 +8,7 @@
     <add_link source_path="script_configuration_dir" source="scripts/copy_cell_indices_ISC.py" dest="copy_cell_indices_ISC.py"/>
     <copy_file source_path="script_configuration_dir" source="scripts/config_E3SM_coupling_files.ini" dest="config_E3SM_coupling_files.ini"/>
     <add_link source_path="script_configuration_dir" source="scripts/readme_E3SM_coupling_files" dest="readme"/>
+	<add_link source_path="initial_condition_database" source="MeltRatesRignot2013.nc" dest="MeltRatesRignot2013.nc"/>
 
 	<run_script name="run.py">
 		<step executable="python">

testing_and_setup/compass/ocean/global_ocean/scripts/config_E3SM_coupling_files.ini

@@ -73,3 +73,10 @@ enable = true
 grid_Levitus_1x1_scrip_file = /usr/projects/climate/mpeterse/mapping_files/test_SSS_mapping_190821/EC60to30Rev4/genRemapFiles/grid_Levitus_1x1_scrip.nc
 # This file needs to be added to a standard repo. Local copy for now:
 salinity_restoring_input_file = /usr/projects/climate/mpeterse/mapping_files/test_SSS_mapping_190821/EC60to30Rev4/interpSSS/PHC2_salx.2004_08_03.filled_double_precision.nc
+
+[prescribed_ismf]
+# does nothing if ice-shelf cavities are not present
+enable = true
+# path to the Rignot2013 data set.
+# The default is a symlink to the initial conditions database
+rignot_2013_file = ../MeltRatesRignot2013.nc

testing_and_setup/compass/ocean/global_ocean/scripts/create_E3SM_coupling_files.py

@@ -21,8 +21,9 @@
 from mpas_tools.ocean.moc import make_moc_basins_and_transects
 from mpas_tools.io import write_netcdf
 import mpas_tools.conversion
-from pyremap import MpasMeshDescriptor, Remapper, get_lat_lon_descriptor, \
-    get_polar_descriptor
+from pyremap import MpasMeshDescriptor, ProjectionGridDescriptor, Remapper, \
+    get_lat_lon_descriptor, get_polar_descriptor
+import pyproj
 # }}}
 
 
@@ -47,9 +48,10 @@ def main():
                      domain_JRA_Gcase,
                      domain_ne30,
                      mapping_runoff,
-                     salinity_restoring]
+                     salinity_restoring,
+                     prescribed_ismf]
 
-    # clean: Delete all directories
+# clean: Delete all directories
     parser = argparse.ArgumentParser()
     parser.add_argument('--clean', action='store_true')
     parser.add_argument('--ice_shelf_cavities', action='store_true')
@@ -655,6 +657,30 @@ def salinity_restoring(config, mesh_name, date_string, ice_shelf_cavities):
 # }}}
 
 
+def prescribed_ismf(config, mesh_name, date_string, ice_shelf_cavities):  # {{{
+    if not ice_shelf_cavities:
+        return
+
+    in_filename = config.get('prescribed_ismf', 'rignot_2013_file')
+
+    out_filename = 'prescriped_ismf_rignot2013.{}.{}.nc'.format(
+        mesh_name, date_string)
+
+    mpiTasks = config.getint('main', 'nprocs')
+
+    remap_rignot(inFileName=in_filename, meshFileName='../init.nc',
+                 meshName=mesh_name, outFileName=out_filename,
+                 mappingDirectory='.', method='conserve',
+                 renormalizationThreshold=None, inVarName='melt_actual',
+                 mpiTasks=mpiTasks)
+
+    output_dir = '../assembled_files_for_upload/inputdata/ocn/mpas-o/{}'.format(
+        mesh_name)
+    make_link('../../../../../prescribed_ismf/{}'.format(out_filename),
+              '{}/{}'.format(output_dir, out_filename))
+    # }}}
+
+
 def make_dir(dirName):
 # {{{
     try:
@@ -1013,6 +1039,97 @@ def make_analysis_polar_map(config, mesh_name, projection):
     # }}}
 
 
+def remap_rignot(inFileName, meshFileName, meshName, outFileName,
+                 mappingDirectory='.', method='conserve',
+                 renormalizationThreshold=None, inVarName='melt_actual',
+                 mpiTasks=1):
+    # {{{
+
+    """
+    Remap the Rignot et al. (2013) melt rates at 1 km resolution to an MPAS
+    mesh
+
+    Parameters
+    ----------
+    inFileName : str
+        The original Rignot et al. (2013) melt rates
+
+    meshFileName : str
+        The MPAS mesh
+
+    meshName : str
+        The name of the mesh (e.g. oEC60to30wISC), used in the name of the
+        mapping file
+
+    outFileName : str
+        The melt rates interpolated to the MPAS mesh with ocean sensible heat
+        fluxes added on (assuming insulating ice)
+
+    mappingDirectory : str
+        The directory where the mapping file should be stored (if it is to be
+        computed) or where it already exists (if not)
+
+    method : {'bilinear', 'neareststod', 'conserve'}, optional
+        The method of interpolation used, see documentation for
+        `ESMF_RegridWeightGen` for details.
+
+    renormalizationThreshold : float, optional
+        The minimum weight of a denstination cell after remapping, below
+        which it is masked out, or ``None`` for no renormalization and
+        masking.
+
+    inVarName : {'melt_actual', 'melt_steadystate'}
+        Whether to use the melt rate for the time period covered in Rignot et
+        al. (2013) with observed thinning/thickening or the melt rates that
+        would be required if ice shelves were in steady state.
+
+    mpiTasks : int, optional
+        The number of MPI tasks to use to compute the mapping file
+    """
+
+    ds = xr.open_dataset(inFileName)
+    lx = np.abs(1e-3 * (ds.xaxis.values[-1] - ds.xaxis.values[0]))
+    ly = np.abs(1e-3 * (ds.yaxis.values[-1] - ds.yaxis.values[0]))
+
+    inGridName = '{}x{}km_1.0km_Antarctic_stereo'.format(lx, ly)
+
+    projection = pyproj.Proj('+proj=stere +lat_ts=-71.0 +lat_0=-90 +lon_0=0.0 '
+                             '+k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84')
+
+    inDescriptor = ProjectionGridDescriptor.read(
+        projection,  inFileName, xVarName='xaxis', yVarName='yaxis',
+        meshName=inGridName)
+
+    # convert to the units and variable names expected in MPAS-O
+    rho_fw = 1000.
+    s_per_yr = 365.*24.*60.*60.
+    latent_heat_of_fusion = 3.337e5
+    ds['prescribedLandIceFreshwaterFlux'] = ds[inVarName]*rho_fw/s_per_yr
+    ds['prescribedLandIceHeatFlux'] = (latent_heat_of_fusion *
+                                       ds['prescribedLandIceFreshwaterFlux'])
+    ds = ds.drop_vars(['melt_actual', 'melt_steadystate', 'lon', 'lat'])
+
+    outDescriptor = MpasMeshDescriptor(meshFileName, meshName)
+
+    mappingFileName = '{}/map_{}_to_{}.nc'.format(mappingDirectory, inGridName,
+                                                  meshName)
+
+    remapper = Remapper(inDescriptor, outDescriptor, mappingFileName)
+
+    remapper.build_mapping_file(method=method, mpiTasks=mpiTasks)
+
+    dsRemap = remapper.remap(
+        ds, renormalizationThreshold=renormalizationThreshold)
+
+    for field in ['prescribedLandIceFreshwaterFlux',
+                  'prescribedLandIceHeatFlux']:
+        # zero out the field where it's currently NaN
+        dsRemap[field] = dsRemap[field].where(dsRemap[field].nonnull(), 0.)
+
+    dsRemap.attrs['history'] = ' '.join(sys.argv)
+    write_netcdf(dsRemap, outFileName)  # }}}
+
+
 if __name__ == '__main__':
     # If called as a primary module, run main
     main()