diff --git a/compass/ocean/tests/global_ocean/__init__.py b/compass/ocean/tests/global_ocean/__init__.py index 69f3e9666..da065fa12 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 358b987f7..cf210b68c 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 0ed632447..e8ff89682 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,15 @@ 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) + self.add_step(remap_step) self.add_step(CullMeshStep( @@ -142,6 +178,13 @@ def configure(self, config=None): config.add_from_package('compass.ocean.mesh', 'remap_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 not self.high_res_topography: filename = 'BedMachineAntarctica_v2_and_GEBCO_2022_0.05_degree_20220729.nc' # noqa: E501 description = 'Bathymetry is from GEBCO 2022, combined with ' \ 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 000000000..226c2f24f --- /dev/null +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py @@ -0,0 +1,252 @@ +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 + """ + + def __init__(self, test_case, base_mesh_step, mesh_name, mali_ais_topo): + """ + 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 + """ + super().__init__(test_case=test_case, base_mesh_step=base_mesh_step, + mesh_name=mesh_name) + self.mali_ais_topo = mali_ais_topo + + 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.nc', + target=mali_filename, + database='mali_topo') + + def run(self): + """ + Run this step of the test case + """ + super().run() + self._remap_mali_topo() + self._combine_topo() + + def _remap_mali_topo(self): + config = self.config + logger = self.logger + + in_mesh_name = self.mali_ais_topo + in_descriptor = MpasCellMeshDescriptor(fileName='mali_topography.nc', + meshName=in_mesh_name) + in_descriptor.format = 'NETCDF3_64BIT' + in_descriptor.to_scrip('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') + + args = ['mbtempest', + '--type', '5', + '--load', 'mali_scrip.nc', + '--load', 'mpaso_scrip.nc', + '--intx', 'moab_intx_mali_mpaso.h5m', + '--rrmgrids'] + + run_command(args=args, cpus_per_task=self.cpus_per_task, + ntasks=self.ntasks, openmp_threads=self.openmp_threads, + config=config, logger=logger) + + ds_mali = xr.open_dataset('mali_topography.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') + + mali_ice_density = ds_mali.attrs['config_ice_density'] + mali_ocean_density = ds_mali.attrs['config_ocean_density'] + sea_level = ds_mali.attrs['config_sea_level'] + + g = constants['SHR_CONST_G'] + ocean_density = constants['SHR_CONST_RHOSW'] + + if ice_density != mali_ice_density: + raise ValueError('Ice density from the config option in ' + '[remap_topography] does not match the value ' + 'from MALI config_ice_density') + if ocean_density != mali_ocean_density: + logger.warn('\nWARNING: Ocean density from SHR_CONST_RHOSW does ' + 'not match the value from MALI config_ocean_density\n') + + draft = - (ice_density / ocean_density) * thickness + + ice_mask = ds_mali.thickness > 0 + floating_mask = np.logical_and( + ice_mask, + ice_density / mali_ocean_density * thickness <= sea_level - bed) + grounded_mask = np.logical_and(ice_mask, np.logical_not(floating_mask)) + ocean_mask = np.logical_and(np.logical_not(ice_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') + + ds_in = {} + # we will remap topography bilinearly + ds_in['bilinear'] = xr.Dataset() + ds_in['bilinear']['bed_elevation'] = bed + ds_in['bilinear']['landIceThkObserved'] = thickness + ds_in['bilinear']['landIcePressureObserved'] = lithop + ds_in['bilinear']['landIceDraftObserved'] = draft + ds_in['bilinear']['landIceFracBilinear'] = ice_frac + ds_in['bilinear']['landIceGroundedFracBilinear'] = grounded_frac + ds_in['bilinear']['oceanFracBilinear'] = ocean_frac + ds_in['bilinear']['maliFracBilinear'] = mali_frac + # we will remap masks conservatively + ds_in['conserve'] = xr.Dataset() + ds_in['conserve']['landIceFracConserve'] = ice_frac + ds_in['conserve']['landIceGroundedFracConserve'] = grounded_frac + ds_in['conserve']['oceanFracConserve'] = ocean_frac + ds_in['conserve']['maliFracConserve'] = mali_frac + + mbtempest_args = {'conserve': ['--order', '1', + '--order', '1', + '--fvmethod', 'none', + '--rrmgrids'], + 'bilinear': ['--order', '1', + '--order', '1', + '--fvmethod', 'bilin']} + + suffix = {'conserve': 'mbtraave', + 'bilinear': 'mbtrbilin'} + + ds_out = xr.Dataset() + for method in ['conserve', 'bilinear']: + mapping_file_name = \ + f'map_{in_mesh_name}_to_{out_mesh_name}_{suffix[method]}.nc' + + # split the parallel executable into constituents in case it + # includes flags + args = ['mbtempest', + '--type', '5', + '--load', 'mali_scrip.nc', + '--load', 'mpaso_scrip.nc', + '--intx', 'moab_intx_mali_mpaso.h5m', + '--weights', + '--method', 'fv', + '--method', 'fv', + '--file', mapping_file_name] + mbtempest_args[method] + + if method == 'bilinear': + # unhappy in parallel for now + check_call(args, logger) + else: + + run_command(args=args, cpus_per_task=self.cpus_per_task, + ntasks=self.ntasks, + openmp_threads=self.openmp_threads, + config=config, logger=logger) + + in_filename = f'mali_topography_{method}.nc' + out_filename = f'mali_topography_ncremap_{method}.nc' + write_netcdf(ds_in[method], in_filename) + + # remapping with the -P mpas flag leads to undesired + # renormalization + args = ['ncremap', + '-m', mapping_file_name, + '--vrb=1', + in_filename, + out_filename] + check_call(args, logger) + + ds_remapped = xr.open_dataset(out_filename) + 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') + + 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] + + # for now, blend topography at calving fronts, but we may want a + # smoother blend in the future + alpha = ds_mali.landIceFracBilinear + ds_out['bed_elevation'] = (alpha * ds_mali.bed_elevation + + (1.0 - alpha) * ds_bedmachine.bed_elevation) + + for var in ['maliFracBilinear', 'maliFracConserve']: + mali_field = ds_mali[var] + mali_field = xr.where(mali_field.notnull(), mali_field, 0.) + ds_out[var] = mali_field + + 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_4to20km.cfg b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_4to20km.cfg new file mode 100644 index 000000000..e85b4a0ea --- /dev/null +++ b/compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_4to20km.cfg @@ -0,0 +1,5 @@ +# 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_20230105_relaxed_10yrs.nc