diff --git a/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_e3sm_coupling.xml b/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_e3sm_coupling.xml index 138bef8220..d3467e19ae 100644 --- a/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_e3sm_coupling.xml +++ b/testing_and_setup/compass/ocean/global_ocean/config_files_ISC/config_e3sm_coupling.xml @@ -8,6 +8,7 @@ + diff --git a/testing_and_setup/compass/ocean/global_ocean/scripts/config_E3SM_coupling_files.ini b/testing_and_setup/compass/ocean/global_ocean/scripts/config_E3SM_coupling_files.ini index b8005d85b7..f9b5648501 100644 --- a/testing_and_setup/compass/ocean/global_ocean/scripts/config_E3SM_coupling_files.ini +++ b/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 diff --git a/testing_and_setup/compass/ocean/global_ocean/scripts/create_E3SM_coupling_files.py b/testing_and_setup/compass/ocean/global_ocean/scripts/create_E3SM_coupling_files.py index 6dffe193f8..f129385867 100755 --- a/testing_and_setup/compass/ocean/global_ocean/scripts/create_E3SM_coupling_files.py +++ b/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()