From deb03b282bfef8857a3eb51b2c677c86b3240dd1 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 10 Jun 2024 03:42:57 -0500 Subject: [PATCH] Add a step for remapping topography from MALI MALI topogrpahy gets combined with BedMachine Antarctica/GEBCO. For now, the masks for the grounding line, calving front, etc. are form BedMachine, not MALI, so that the same horizontal mesh and mapping files can be used with BedMachine and MALI topogrpahy. --- compass/ocean/tests/global_ocean/__init__.py | 7 +- .../ocean/tests/global_ocean/init/__init__.py | 3 +- .../ocean/tests/global_ocean/mesh/__init__.py | 53 +++- .../mesh/remap_mali_topography/__init__.py | 252 ++++++++++++++++++ .../remap_mali_topography/ais_4to20km.cfg | 5 + 5 files changed, 311 insertions(+), 9 deletions(-) create mode 100644 compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py create mode 100644 compass/ocean/tests/global_ocean/mesh/remap_mali_topography/ais_4to20km.cfg 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 5d092dc506..b907407412 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 0ed6324477..e8ff89682f 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 0000000000..226c2f24fe --- /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 0000000000..e85b4a0ea5 --- /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