diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index 69f3e96669..da065fa122 100644 --- a/compass/ocean/tests/global_ocean/__init__.py +++ b/compass/ocean/tests/global_ocean/__init__.py @@ -45,6 +45,8 @@ def __init__(self, mpas_core): self._add_tests(mesh_names=['ARRM10to60', 'ARRMwISC10to60']) self._add_tests(mesh_names=['SO12to30', 'SOwISC12to30']) + self._add_tests(mesh_names=['SOwISC12to30'], + mali_ais_topo='AIS_4to20km') self._add_tests(mesh_names=['WC14', 'WCwISC14']) @@ -65,7 +67,7 @@ def __init__(self, mpas_core): def _add_tests(self, mesh_names, high_res_topography=True, include_rk4=False, include_regression=False, include_phc=False, - include_en4_1900=False): + include_en4_1900=False, mali_ais_topo=None): """ Add test cases for the given mesh(es) """ default_ic = 'WOA23' @@ -73,7 +75,8 @@ def _add_tests(self, mesh_names, high_res_topography=True, for mesh_name in mesh_names: mesh_test = Mesh(test_group=self, mesh_name=mesh_name, - high_res_topography=high_res_topography) + high_res_topography=high_res_topography, + mali_ais_topo=mali_ais_topo) self.add_test_case(mesh_test) init_test = Init(test_group=self, mesh=mesh_test, diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index 358b987f74..cf210b68c8 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -44,9 +44,8 @@ def __init__(self, test_group, mesh, initial_condition): The initial condition dataset to use """ name = 'init' - mesh_name = mesh.mesh_name ic_dir = initial_condition - self.init_subdir = os.path.join(mesh_name, ic_dir) + self.init_subdir = os.path.join(mesh.mesh_subdir, ic_dir) subdir = os.path.join(self.init_subdir, name) super().__init__(test_group=test_group, name=name, subdir=subdir) diff --git a/compass/ocean/tests/global_ocean/mesh/__init__.py b/compass/ocean/tests/global_ocean/mesh/__init__.py index 2d082c3b83..e5092e0134 100644 --- a/compass/ocean/tests/global_ocean/mesh/__init__.py +++ b/compass/ocean/tests/global_ocean/mesh/__init__.py @@ -1,3 +1,5 @@ +import os + from compass.mesh.spherical import ( IcosahedralMeshStep, QuasiUniformSphericalMeshStep, @@ -15,6 +17,9 @@ IcosMeshFromConfigStep, QUMeshFromConfigStep, ) +from compass.ocean.tests.global_ocean.mesh.remap_mali_topography import ( + RemapMaliTopography, +) from compass.ocean.tests.global_ocean.mesh.rrs6to18 import RRS6to18BaseMesh from compass.ocean.tests.global_ocean.mesh.so12to30 import SO12to30BaseMesh from compass.ocean.tests.global_ocean.mesh.wc14 import WC14BaseMesh @@ -31,6 +36,13 @@ class Mesh(TestCase): Attributes ---------- + mesh_name : str + The name of the mesh + + mesh_subdir : str + The subdirectory within the test group for all test cases with this + mesh and topography + package : str The python package for the mesh @@ -43,8 +55,14 @@ class Mesh(TestCase): high_res_topography : bool Whether to remap a high resolution topography data set. A lower res data set is used for low resolution meshes. + + mali_ais_topo : str + Short name for the MALI dataset to use for Antarctic Ice Sheet + topography """ - def __init__(self, test_group, mesh_name, high_res_topography): + + def __init__(self, test_group, mesh_name, # noqa: C901 + high_res_topography, mali_ais_topo=None): """ Create test case for creating a global MPAS-Ocean mesh @@ -59,9 +77,20 @@ def __init__(self, test_group, mesh_name, high_res_topography): high_res_topography : bool Whether to remap a high resolution topography data set. A lower res data set is used for low resolution meshes. + + mali_ais_topo : str, optional + Short name for the MALI dataset to use for Antarctic Ice Sheet + topography """ name = 'mesh' - subdir = f'{mesh_name}/{name}' + if mali_ais_topo is None: + self.mesh_subdir = mesh_name + else: + self.mesh_subdir = os.path.join(mesh_name, + f'MALI_topo_{mali_ais_topo}') + + subdir = os.path.join(self.mesh_subdir, name) + super().__init__(test_group=test_group, name=name, subdir=subdir) with_ice_shelf_cavities = 'wISC' in mesh_name @@ -75,6 +104,7 @@ def __init__(self, test_group, mesh_name, high_res_topography): self.mesh_config_filename = f'{mesh_lower}.cfg' self.mesh_name = mesh_name + self.mali_ais_topo = mali_ais_topo self.with_ice_shelf_cavities = with_ice_shelf_cavities self.high_res_topography = high_res_topography @@ -117,9 +147,16 @@ def __init__(self, test_group, mesh_name, high_res_topography): self.add_step(base_mesh_step) - remap_step = RemapTopography(test_case=self, - base_mesh_step=base_mesh_step, - mesh_name=mesh_name) + if mali_ais_topo is None: + remap_step = RemapTopography(test_case=self, + base_mesh_step=base_mesh_step, + mesh_name=mesh_name) + else: + remap_step = RemapMaliTopography( + test_case=self, base_mesh_step=base_mesh_step, + mesh_name=mesh_name, mali_ais_topo=mali_ais_topo, + ocean_includes_grounded=False) + self.add_step(remap_step) self.add_step(CullMeshStep( @@ -146,6 +183,13 @@ def configure(self, config=None): 'low_res_topography.cfg', exception=True) + if self.mali_ais_topo is not None: + package = 'compass.ocean.tests.global_ocean.mesh.' \ + 'remap_mali_topography' + config.add_from_package(package, + f'{self.mali_ais_topo.lower()}.cfg', + exception=True) + if self.mesh_name.startswith('Kuroshio'): # add the config options for all kuroshio meshes config.add_from_package( diff --git a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py new file mode 100644 index 0000000000..0293715bab --- /dev/null +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py @@ -0,0 +1,268 @@ +import os + +import numpy as np +import xarray as xr +from mpas_tools.cime.constants import constants +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call +from pyremap import MpasCellMeshDescriptor + +from compass.ocean.mesh.remap_topography import RemapTopography +from compass.parallel import run_command + + +class RemapMaliTopography(RemapTopography): + """ + A step for remapping bathymetry and ice-shelf topography from a MALI + datase to a global MPAS-Ocean mesh and combining it with a global + topography dataset + + Attributes + ---------- + mali_ais_topo : str + Short name for the MALI dataset to use for Antarctic Ice Sheet + topography + + ocean_includes_grounded : bool, optional + Whether to include grounded cells that are below sea level in the + ocean domain + """ + + def __init__(self, test_case, base_mesh_step, mesh_name, mali_ais_topo, + ocean_includes_grounded): + """ + Create a new step + + Parameters + ---------- + test_case : compass.ocean.tests.global_ocean.mesh.Mesh + The test case this step belongs to + + base_mesh_step : compass.mesh.spherical.SphericalBaseStep + The base mesh step containing input files to this step + + remap_topography : compass.ocean.mesh.remap_topography.RemapTopography + A step for remapping topography. If provided, the remapped + topography is used to determine the land mask + + mesh_name : str + The name of the MPAS mesh to include in the mapping file + + mali_ais_topo : str, optional + Short name for the MALI dataset to use for Antarctic Ice Sheet + topography + + ocean_includes_grounded : bool + Whether to include grounded cells that are below sea level in the + ocean domain + """ + super().__init__(test_case=test_case, base_mesh_step=base_mesh_step, + mesh_name=mesh_name) + self.mali_ais_topo = mali_ais_topo + self.ocean_includes_grounded = ocean_includes_grounded + + self.add_output_file(filename='mali_topography_remapped.nc') + + def setup(self): + """ + Set up the step in the work directory, including downloading any + dependencies. + """ + super().setup() + + mali_filename = self.config.get('remap_mali_topography', + 'mali_filename') + self.add_input_file( + filename='mali_topography_orig.nc', + target=mali_filename, + database='mali_topo') + + if not self.config.has_option('remap_topography', 'ocean_density'): + raise ValueError( + 'You must be using a mesh that defines [remap_topography]/' + 'ocean_density in the config file') + + def run(self): + """ + Run this step of the test case + """ + super().run() + self._remap_mali_topo() + self._combine_topo() + + def _remap_mali_topo(self): + in_mesh_name = self.mali_ais_topo + in_descriptor = MpasCellMeshDescriptor( + fileName='mali_topography_orig.nc', + meshName=in_mesh_name) + in_descriptor.format = 'NETCDF3_64BIT' + in_descriptor.to_scrip('mali.scrip.nc') + self._partition_scrip_file('mali.scrip.nc') + + out_mesh_name = self.mesh_name + out_descriptor = MpasCellMeshDescriptor( + fileName='base_mesh.nc', + meshName=out_mesh_name) + out_descriptor.format = 'NETCDF3_64BIT' + out_descriptor.to_scrip('mpaso.scrip.nc') + self._partition_scrip_file('mpaso.scrip.nc') + + map_filename = \ + f'map_{in_mesh_name}_to_{out_mesh_name}_mbtraave.nc' + + self._create_mali_to_mpaso_weights(map_filename) + + self._modify_mali_topo() + + self._remap_mali_to_mpaso(map_filename) + + def _modify_mali_topo(self): + """ + Modify MALI topography to have desired fields and names + """ + logger = self.logger + config = self.config + logger.info('Modifying MALI topography fields and names') + + ds_mali = xr.open_dataset('mali_topography_orig.nc') + if 'Time' in ds_mali.dims: + ds_mali = ds_mali.isel(Time=0) + bed = ds_mali.bedTopography + thickness = ds_mali.thickness + + ice_density = config.getfloat('remap_topography', 'ice_density') + ocean_density = config.getfloat('remap_topography', 'ocean_density') + sea_level = config.getfloat('remap_topography', 'sea_level') + + g = constants['SHR_CONST_G'] + + draft = - (ice_density / ocean_density) * thickness + + ice_mask = ds_mali.thickness > 0 + floating_mask = np.logical_and( + ice_mask, + ice_density / ocean_density * thickness <= sea_level - bed) + grounded_mask = np.logical_and(ice_mask, np.logical_not(floating_mask)) + + if self.ocean_includes_grounded: + ocean_mask = bed < sea_level + else: + ocean_mask = np.logical_and(np.logical_not(grounded_mask), + bed < sea_level) + + lithop = ice_density * g * thickness + ice_frac = xr.where(ice_mask, 1., 0.) + grounded_frac = xr.where(grounded_mask, 1., 0.) + ocean_frac = xr.where(ocean_mask, 1., 0.) + mali_frac = xr.DataArray(data=np.ones(ds_mali.sizes['nCells']), + dims='nCells') + + # we will remap conservatively + ds_in = xr.Dataset() + ds_in['bed_elevation'] = bed + ds_in['landIceThkObserved'] = thickness + ds_in['landIcePressureObserved'] = lithop + ds_in['landIceDraftObserved'] = draft + ds_in['landIceFrac'] = ice_frac + ds_in['landIceGroundedFrac'] = grounded_frac + ds_in['oceanFrac'] = ocean_frac + ds_in['maliFrac'] = mali_frac + + write_netcdf(ds_in, 'mali_topography_mod.nc') + + logger.info(' Done.') + + def _create_mali_to_mpaso_weights(self, map_filename): + """ + Create mapping weights file using TempestRemap + """ + logger = self.logger + logger.info('Create weights file') + + args = [ + 'mbtempest', '--type', '5', + '--load', f'mali.scrip.p{self.ntasks}.h5m', + '--load', f'mpaso.scrip.p{self.ntasks}.h5m', + '--file', map_filename, + '--weights', '--gnomonic', + '--boxeps', '1e-9', + ] + + run_command( + args, self.cpus_per_task, self.ntasks, + self.openmp_threads, self.config, self.logger, + ) + + logger.info(' Done.') + + def _remap_mali_to_mpaso(self, map_filename): + """ + Remap combined bathymetry onto MPAS target mesh + """ + logger = self.logger + logger.info('Remap to target') + + args = [ + 'ncremap', + '-m', map_filename, + '--vrb=1', + 'mali_topography_mod.nc', 'mali_topography_ncremap.nc', + ] + check_call(args, logger) + + ds_remapped = xr.open_dataset('mali_topography_ncremap.nc') + ds_out = xr.Dataset() + for var_name in ds_remapped: + var = ds_remapped[var_name] + if 'Frac' in var: + # copy the fractional variable, making sure it doesn't + # exceed 1 + var = np.minimum(var, 1.) + ds_out[var_name] = var + + write_netcdf(ds_out, 'mali_topography_remapped.nc') + + logger.info(' Done.') + + def _combine_topo(self): + os.rename('topography_remapped.nc', + 'bedmachine_topography_remapped.nc') + ds_mali = xr.open_dataset('mali_topography_remapped.nc') + ds_mali = ds_mali.reset_coords(['lat', 'lon'], drop=True) + ds_bedmachine = xr.open_dataset('bedmachine_topography_remapped.nc') + ds_bedmachine = ds_bedmachine.reset_coords(['lat', 'lon'], drop=True) + + ds_out = xr.Dataset(ds_bedmachine) + # for now, just use the land-ice pressure, draft and thickness from + # MALI + for var in ['landIcePressureObserved', 'landIceDraftObserved', + 'landIceThkObserved']: + ds_out[var] = ds_mali[var] + + ds_out['landIceFracObserved'] = ds_mali['landIceFrac'] + + ds_out['landIceFloatingFracObserved'] = ( + ds_mali['landIceFrac'] - + ds_mali['landIceGroundedFrac']) + + for var in ['maliFrac']: + mali_field = ds_mali[var] + mali_field = xr.where(mali_field.notnull(), mali_field, 0.) + ds_out[var] = mali_field + + # for now, blend topography at calving fronts, but we may want a + # smoother blend in the future + alpha = ds_mali.landIceFrac + ds_out['bed_elevation'] = ( + alpha * ds_mali.bed_elevation + + (1.0 - alpha) * ds_bedmachine.bed_elevation) + + alpha = ds_out.maliFrac + # NOTE: MALI's ocean fraction is already scaled by the MALI fraction + ds_out['oceanFracObserved'] = ( + ds_mali.oceanFrac + + (1.0 - alpha) * ds_bedmachine.oceanFracObserved) + + ds_out['ssh'] = ds_out.landIceDraftObserved + + write_netcdf(ds_out, 'topography_remapped.nc') diff --git a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_2to10km.cfg b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_2to10km.cfg new file mode 100644 index 0000000000..571248ab29 --- /dev/null +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_2to10km.cfg @@ -0,0 +1,18 @@ +# config options related to remapping topography to an MPAS-Ocean mesh +[remap_topography] + +# the density of land ice from MALI (kg/m^3) +ice_density = 910.0 + +# the density of ocean water (kg/m^3), equivalent to config_ocean_density in MALI +ocean_density = 1028.0 + +# the elevation of sea level, equivalent to config_sea_level in MALI +sea_level = 0.0 + + +# config options related to remapping MALI topography to an MPAS-Ocean mesh +[remap_mali_topography] + +# The name of the MALI topogrpahy file +mali_filename = ais_2to10km_20231222_modifiedBathymetry-drop20m-smooth6x_bulldozedTroughs-75-400m_BEDMAP2-ASE-mjh.nc diff --git a/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_4to20km.cfg b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_4to20km.cfg new file mode 100644 index 0000000000..c21671567f --- /dev/null +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_4to20km.cfg @@ -0,0 +1,18 @@ +# config options related to remapping topography to an MPAS-Ocean mesh +[remap_topography] + +# the density of land ice from MALI (kg/m^3) +ice_density = 910.0 + +# the density of ocean water (kg/m^3), equivalent to config_ocean_density in MALI +ocean_density = 1028.0 + +# the elevation of sea level, equivalent to config_sea_level in MALI +sea_level = 0.0 + + +# config options related to remapping MALI topography to an MPAS-Ocean mesh +[remap_mali_topography] + +# The name of the MALI topogrpahy file +mali_filename = AIS_4to20km_20241118_relaxation_0TGmelt_10yr.20241224.nc diff --git a/compass/ocean/tests/global_ocean/mesh/so12to30/so12to30.cfg b/compass/ocean/tests/global_ocean/mesh/so12to30/so12to30.cfg index 3f7551ae34..794bebef06 100644 --- a/compass/ocean/tests/global_ocean/mesh/so12to30/so12to30.cfg +++ b/compass/ocean/tests/global_ocean/mesh/so12to30/so12to30.cfg @@ -59,13 +59,13 @@ mesh_description = MPAS Southern Ocean regionally refined mesh for E3SM version e3sm_version = 3 # The revision number of the mesh, which should be incremented each time the # mesh is revised -mesh_revision = 3 +mesh_revision = 4 # the minimum (finest) resolution in the mesh min_res = 12 # the maximum (coarsest) resolution in the mesh, can be the same as min_res max_res = 30 # The URL of the pull request documenting the creation of the mesh -pull_request = https://github.com/MPAS-Dev/compass/pull/807 +pull_request = https://github.com/MPAS-Dev/compass/pull/829 # config options related to initial condition and diagnostics support files diff --git a/docs/conf.py b/docs/conf.py index 8fbbddf907..7679e5667a 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -53,10 +53,15 @@ # General information about the project. project = u'compass' -copyright = u'Copyright (c) 2013-2021, Los Alamos National Security, LLC (LANS) (Ocean: LA-CC-13-047;' \ - u'Land Ice: LA-CC-13-117) and the University Corporation for Atmospheric Research (UCAR).' -author = u'Xylar Asay-Davis, Matt Hoffman, Doug Jacobsen, Mark Petersen, ' \ - u'Philip Wolfram, Luke Van Roekel, Tong Zhang' +copyright = u'Copyright (c) 2013-2024, Triad National Security, LLC ' \ + u'(Ocean: LA-CC-13-047; Land Ice: LA-CC-13-117) and the ' \ + u'University Corporation for Atmospheric Research (UCAR).' +author = u'Xylar Asay-Davis, Alice Barthel, Carolyn Begeman, Pete Bosler, ' \ + u'Riley Brady, Steven Brus, Sara Calandrini, Zhendong Cao, ' \ + u'Giacomo Capodaglio, Max Carlson, Althea Denlinger, ' \ + u'Yaris Eidenbenz, Darren Engwirda, Holly Han, Jeremy Lilly, ' \ + u'Mauro Perego, Mark Petersen, Cameron Smith, Yohei Takano, ' \ + u'Irena Vankova, Luke Van Roekel, Philip Wolfram, Tong Zhang' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the diff --git a/docs/developers_guide/ocean/api.rst b/docs/developers_guide/ocean/api.rst index a2e1845d54..d0d2b588ca 100644 --- a/docs/developers_guide/ocean/api.rst +++ b/docs/developers_guide/ocean/api.rst @@ -284,6 +284,10 @@ test cases and steps mesh.Mesh.configure mesh.Mesh.run + mesh.remap_mali_topography.RemapMaliTopography + mesh.remap_mali_topography.RemapMaliTopography.setup + mesh.remap_mali_topography.RemapMaliTopography.run + mesh.ec30to60.EC30to60BaseMesh mesh.ec30to60.EC30to60BaseMesh.build_cell_width_lat_lon diff --git a/docs/developers_guide/ocean/test_groups/global_ocean.rst b/docs/developers_guide/ocean/test_groups/global_ocean.rst index 85fc30a5a2..df4fb7f98d 100644 --- a/docs/developers_guide/ocean/test_groups/global_ocean.rst +++ b/docs/developers_guide/ocean/test_groups/global_ocean.rst @@ -874,6 +874,41 @@ The vertical grid is an ``index_tanh_dz`` profile (see :ref:`dev_ocean_framework_vertical`) with 64 vertical levels ranging in thickness from 10 to 250 m. +.. _dev_ocean_global_ocean_sowisc12to30_mali_topo: + +SOwISC12to30 with MALI topography ++++++++++++++++++++++++++++++++++ + +Versions of the ``SOwISC12to30`` test cases exist with topography that comes +from the MALI model instead of from GEBCO and BedMachine Antarctica. + +For now, the only supported MALI initial condition is a 4-20 km resoution +mesh with the following config options: + +.. code-block:: cfg + + # config options related to remapping topography to an MPAS-Ocean mesh + [remap_topography] + + # the density of ocean water (kg/m^3), equivalent to config_ocean_density in MALI + ocean_density = 1028.0 + + # the elevation of sea level, equivalent to config_sea_level in MALI + sea_level = 0.0 + + + # config options related to remapping MALI topography to an MPAS-Ocean mesh + [remap_mali_topography] + + # The name of the MALI topogrpahy file + mali_filename = ais_4to20km.20240708.nc + + +The topography is remapped in a step defined in the +:py:class:`compass.ocean.tests.global_ocean.mesh.remap_mali_topography.RemapMaliTopography` +class, which descends from +:py:class:`compass.ocean.mesh.remap_topography.RemapTopography`. + .. _dev_ocean_global_ocean_wc14: WC14 and WCwISC14 diff --git a/docs/users_guide/ocean/test_groups/global_ocean.rst b/docs/users_guide/ocean/test_groups/global_ocean.rst index 94be78db92..3e15b613b1 100644 --- a/docs/users_guide/ocean/test_groups/global_ocean.rst +++ b/docs/users_guide/ocean/test_groups/global_ocean.rst @@ -476,6 +476,10 @@ The mesh has 12-km resolution around Antarctica and 30-km resolution elsewhere. The mesh includes the :ref:`global_ocean_ice_shelf_cavities` around Antarctica in the ocean domain. +A version of the SOwISC12to30 mesh also exist with topography that comes +from a 4-20 km MALI initial condition instead of from GEBCO and BedMachine +Antarctica. + .. image:: images/sowisc12to60.png :width: 500 px :align: center @@ -537,7 +541,7 @@ As a user, your main way of altering forward runs is by changing namelist options directly in ``namelist.ocean`` or modifying streams in ``streams.ocean``. However, there are a few parameters related to forward runs you can change in the config file for a test case. Since some test cases like -:ref:`global_ocean_restart_test` and :ref`global_ocean_dynamic_adjustment` have +:ref:`global_ocean_restart_test` and :ref:`global_ocean_dynamic_adjustment` have more than one forward run, it is convenient to change options like ``forward_ntasks`` once in the config file, knowing that this will change the target number of cores of all forward model runs in the test case. The same