For initialization with ice-shelf cavities, adjust SSH not land ice pressure

compass/ocean/iceshelf.py

@@ -87,7 +87,7 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
         on_a_sphere = ds.attrs['on_a_sphere'].lower() == 'yes'
 
         modify_mask = numpy.logical_and(ds.maxLevelCell > 0,
-                                        ds.modifyLandIcePressureMask == 1)
+                                        ds.sshAdjustmentMask == 1)
 
     for iter_index in range(iteration_count):
         logger.info(f" * Iteration {iter_index + 1}/{iteration_count}")

compass/ocean/mesh/cull.py

@@ -501,9 +501,7 @@ def _land_mask_from_topo(with_cavities, topo_filename, mask_filename):
         # we want the mask to be 1 where there's not ocean
         cull_mask = xr.where(ocean_frac < 0.5, 1, 0)
     else:
-        land_ice_frac = ds_topo.landIceFracObserved
-        grounded_ice_frac = ds_topo.landIceGroundedFracObserved
-        floating_ice_frac = land_ice_frac - grounded_ice_frac
+        floating_ice_frac = ds_topo.landIceFloatingFracObserved
         no_cavities_ocean_frac = ocean_frac - floating_ice_frac
 
         # we want the mask to be 1 where there's not open ocean
@@ -558,7 +556,7 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh, logger):
     land_ice_draft_mask = ds_mask.cellSeedMask
     ds_topo['landIceDraftMask'] = land_ice_draft_mask
     ds_topo['landIceDraftObserved'] = (
-        land_ice_mask * ds_topo.landIceDraftObserved)
+        land_ice_draft_mask * ds_topo.landIceDraftObserved)
 
 
 def _land_mask_from_geojson(with_cavities, process_count, logger,

compass/ocean/mesh/remap_topography.cfg

@@ -2,17 +2,17 @@
 [remap_topography]
 
 # the name of the topography file in the bathymetry database
-topo_filename = BedMachineAntarctica_v3_and_GEBCO_2023_0.0125_degree_20230831.nc
+topo_filename = BedMachineAntarctica_v3_and_GEBCO_2023_0.0125_degree_20240611.nc
 
 # variable names in topo_filename
 lon_var = lon
 lat_var = lat
 bathymetry_var = bathymetry
-ice_draft_var = ice_draft
 ice_thickness_var = thickness
 ice_frac_var = ice_mask
 grounded_ice_frac_var = grounded_mask
 ocean_frac_var = ocean_mask
+bathy_frac_var = bathymetry_mask
 
 # the description to include in metadata
 description = Bathymetry is from GEBCO 2023, combined with BedMachine
@@ -28,3 +28,6 @@ method = conserve
 # threshold of what fraction of an MPAS cell must contain ocean in order to
 # perform renormalization of elevation variables
 renorm_threshold = 0.01
+
+# the density of land ice from MALI (kg/m^3)
+ice_density = 910.0

compass/ocean/mesh/remap_topography.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import xarray as xr
+from mpas_tools.cime.constants import constants
 from mpas_tools.io import write_netcdf
 from pyremap import LatLonGridDescriptor, MpasCellMeshDescriptor, Remapper
 
@@ -57,7 +58,8 @@ def setup(self):
         dependencies.
         """
         super().setup()
-        topo_filename = self.config.get('remap_topography', 'topo_filename')
+        config = self.config
+        topo_filename = config.get('remap_topography', 'topo_filename')
         self.add_input_file(
             filename='topography.nc',
             target=topo_filename,
@@ -69,7 +71,6 @@ def setup(self):
         target = os.path.join(base_path, base_filename)
         self.add_input_file(filename='base_mesh.nc', work_dir_target=target)
 
-        config = self.config
         self.ntasks = config.getint('remap_topography', 'ntasks')
         self.min_tasks = config.getint('remap_topography', 'min_tasks')
 
@@ -101,6 +102,9 @@ def run(self):
         method = config.get('remap_topography', 'method')
         renorm_threshold = config.getfloat('remap_topography',
                                            'renorm_threshold')
+        ice_density = config.getfloat('remap_topography', 'ice_density')
+        ocean_density = constants['SHR_CONST_RHOSW']
+        g = constants['SHR_CONST_G']
 
         in_descriptor = LatLonGridDescriptor.read(fileName='topography.nc',
                                                   lonVarName=lon_var,
@@ -128,27 +132,44 @@ def run(self):
         ds_in = ds_in.rename({'ncol': 'nCells'})
         ds_out = xr.Dataset()
         rename = {'bathymetry_var': 'bed_elevation',
-                  'ice_draft_var': 'landIceDraftObserved',
                   'ice_thickness_var': 'landIceThkObserved',
                   'ice_frac_var': 'landIceFracObserved',
                   'grounded_ice_frac_var': 'landIceGroundedFracObserved',
-                  'ocean_frac_var': 'oceanFracObserved'}
+                  'ocean_frac_var': 'oceanFracObserved',
+                  'bathy_frac_var': 'bathyFracObserved'}
 
-        for option in rename:
+        for option, out_var in rename.items():
             in_var = config.get('remap_topography', option)
-            out_var = rename[option]
             ds_out[out_var] = ds_in[in_var]
 
+        ds_out['landIceFloatingFracObserved'] = \
+            ds_out.landIceFracObserved - ds_out.landIceGroundedFracObserved
+
         # make sure fractions don't exceed 1
         for var in ['landIceFracObserved', 'landIceGroundedFracObserved',
-                    'oceanFracObserved']:
+                    'landIceFloatingFracObserved', 'oceanFracObserved',
+                    'bathyFracObserved']:
             ds_out[var] = np.minimum(ds_out[var], 1.)
 
         # renormalize elevation variables
-        norm = ds_out.oceanFracObserved
+        norm = ds_out.bathyFracObserved
         valid = norm > renorm_threshold
-        for var in ['bed_elevation', 'landIceDraftObserved',
-                    'landIceThkObserved']:
+        for var in ['bed_elevation', 'landIceThkObserved']:
             ds_out[var] = xr.where(valid, ds_out[var] / norm, 0.)
 
+        thickness = ds_out.landIceThkObserved
+        ds_out['landIcePressureObserved'] = ice_density * g * thickness
+
+        # compute the ice draft to be consistent with the land ice pressure
+        # and using E3SM's density of seawater
+        draft = - (ice_density / ocean_density) * thickness
+        bed = ds_out.bed_elevation
+
+        # can't be deeper than the bed
+        draft = xr.where(draft >= bed, draft, bed)
+
+        ds_out['landIceDraftObserved'] = draft
+
+        ds_out['ssh'] = draft
+
         write_netcdf(ds_out, 'topography_remapped.nc')

compass/ocean/tests/global_ocean/forward.py

@@ -139,12 +139,8 @@ def __init__(self, test_case, mesh, time_integrator, init=None,
             work_dir_target=f'{mesh_path}/culled_mesh.nc')
 
         if init is not None:
-            if mesh.with_ice_shelf_cavities:
-                initial_state_target = \
-                    f'{init.path}/ssh_adjustment/adjusted_init.nc'
-            else:
-                initial_state_target = \
-                    f'{init.path}/initial_state/initial_state.nc'
+            initial_state_target = \
+                f'{init.path}/initial_state/initial_state.nc'
             self.add_input_file(filename='init.nc',
                                 work_dir_target=initial_state_target)
             self.add_input_file(

compass/ocean/tests/global_ocean/global_ocean.cfg

@@ -21,7 +21,7 @@ sweep_count = 20
 [ssh_adjustment]
 
 # the number of iterations of ssh adjustment to perform
-iterations = 10
+iterations = 4
 
 # Whether to convert adjusted initial condition files to CDF5 format during
 # ssh adjustment under ice shelves
@@ -54,8 +54,10 @@ btr_dt_per_km = 1.5
 min_levels = 3
 cavity_min_levels = ${min_levels}
 
-# minimum thickness of layers in ice-shelf cavities
-cavity_min_layer_thickness = 1.0
+# minimum thickness of layers in ice-shelf cavities at the beginning and end
+# of iterative ssh init
+cavity_min_layer_thickness_initial = 10.0
+cavity_min_layer_thickness_final = 3.0
 
 # Maximum allowed Haney number for configurations with ice-shelf cavities
 rx1_max = 20.0

compass/ocean/tests/global_ocean/init/__init__.py

@@ -5,6 +5,9 @@
     RemapIceShelfMelt,
 )
 from compass.ocean.tests.global_ocean.init.ssh_adjustment import SshAdjustment
+from compass.ocean.tests.global_ocean.init.ssh_from_surface_density import (
+    SshFromSurfaceDensity,
+)
 from compass.testcase import TestCase
 from compass.validate import compare_variables
 
@@ -50,29 +53,21 @@ def __init__(self, test_group, mesh, initial_condition):
         self.mesh = mesh
         self.initial_condition = initial_condition
 
-        self.add_step(
-            InitialState(
-                test_case=self, mesh=mesh,
-                initial_condition=initial_condition))
-
-        if mesh.with_ice_shelf_cavities:
-            self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh))
-
-            self.add_step(
-                SshAdjustment(test_case=self))
-
     def configure(self, config=None):
         """
         Modify the configuration options for this test case
 
         config : compass.config.CompassConfigParser, optional
             Configuration options to update if not those for this test case
         """
+        add_steps = config is None
         if config is None:
             config = self.config
 
+        mesh = self.mesh
+
         # set mesh-relate config options
-        self.mesh.configure(config=config)
+        mesh.configure(config=config)
 
         initial_condition = self.initial_condition
         descriptions = {'WOA23': 'World Ocean Atlas 2023 climatology '
@@ -84,6 +79,72 @@ def configure(self, config=None):
         config.set('global_ocean', 'init_description',
                    descriptions[initial_condition])
 
+        if add_steps:
+            # add the steps for ssh adjustment
+            if mesh.with_ice_shelf_cavities:
+                step_index = 1
+                name = \
+                    f'{step_index:02d}_init_with_draft_from_constant_density'
+                subdir = f'adjust_ssh/{name}'
+                init_const_rho = InitialState(
+                    test_case=self, mesh=mesh,
+                    initial_condition=initial_condition,
+                    name=name, subdir=subdir,
+                    adjustment_fraction=0.)
+                self.add_step(init_const_rho)
+
+                # Recompute ssh using surface density
+                step_index += 1
+                name = f'{step_index:02d}_ssh_from_surface_density'
+                subdir = f'adjust_ssh/{name}'
+                ssh_from_surf_rho = SshFromSurfaceDensity(
+                    test_case=self, init_path=init_const_rho.path,
+                    name=name, subdir=subdir)
+                self.add_step(ssh_from_surf_rho)
+
+                culled_topo_path = ssh_from_surf_rho.path
+
+                iteration_count = config.getint('ssh_adjustment', 'iterations')
+                for iter_index in range(iteration_count):
+                    fraction = iter_index / iteration_count
+
+                    step_index += 1
+                    name = f'{step_index:02d}_init'
+                    subdir = f'adjust_ssh/{name}'
+                    init_step = InitialState(
+                        test_case=self, mesh=mesh,
+                        initial_condition=initial_condition,
+                        culled_topo_path=culled_topo_path,
+                        name=name, subdir=subdir,
+                        adjustment_fraction=fraction)
+                    self.add_step(init_step)
+
+                    step_index += 1
+                    name = f'{step_index:02d}_adjust_ssh'
+                    subdir = f'adjust_ssh/{name}'
+                    adjust_ssh = SshAdjustment(
+                        test_case=self, init_path=init_step.path,
+                        name=name, subdir=subdir)
+                    self.add_step(adjust_ssh)
+                    culled_topo_path = adjust_ssh.path
+
+                name = 'initial_state'
+                subdir = 'initial_state'
+                init_step = InitialState(
+                    test_case=self, mesh=mesh,
+                    initial_condition=initial_condition,
+                    culled_topo_path=culled_topo_path,
+                    name=name, subdir=subdir,
+                    adjustment_fraction=1.0)
+                self.add_step(init_step)
+
+                self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh))
+            else:
+                self.add_step(
+                    InitialState(
+                        test_case=self, mesh=mesh,
+                        initial_condition=initial_condition))
+
     def validate(self):
         """
         Test cases can override this method to perform validation of variables
@@ -92,8 +153,3 @@ def validate(self):
         variables = ['temperature', 'salinity', 'layerThickness']
         compare_variables(test_case=self, variables=variables,
                           filename1='initial_state/initial_state.nc')
-
-        if self.mesh.with_ice_shelf_cavities:
-            variables = ['ssh', 'landIcePressure']
-            compare_variables(test_case=self, variables=variables,
-                              filename1='ssh_adjustment/adjusted_init.nc')

compass/ocean/tests/global_ocean/init/initial_state.py

@@ -28,8 +28,14 @@ class InitialState(Step):
 
     initial_condition : {'WOA23', 'PHC', 'EN4_1900'}
         The initial condition dataset to use
+
+    adjustment_fraction : float
+        The fraction of the way through iterative ssh adjustment for this
+        step
     """
-    def __init__(self, test_case, mesh, initial_condition):
+    def __init__(self, test_case, mesh, initial_condition,
+                 culled_topo_path=None, name='initial_state', subdir=None,
+                 adjustment_fraction=None):
         """
         Create the step
 
@@ -43,13 +49,30 @@ def __init__(self, test_case, mesh, initial_condition):
 
         initial_condition : {'WOA23', 'PHC', 'EN4_1900'}
             The initial condition dataset to use
+
+        culled_topo_path : str, optional
+            The path to a step where ``topography_culled.nc`` is provided
+
+        name : str, optional
+            The name of the step
+
+        subdir : str, optional
+            The subdirectory for the step
+
+        adjustment_fraction : float, optional
+            The fraction of the way through iterative ssh adjustment for this
+            step
         """
         if initial_condition not in ['WOA23', 'PHC', 'EN4_1900']:
             raise ValueError(f'Unknown initial_condition {initial_condition}')
 
-        super().__init__(test_case=test_case, name='initial_state')
+        super().__init__(test_case=test_case, name=name, subdir=subdir)
         self.mesh = mesh
         self.initial_condition = initial_condition
+        if mesh.with_ice_shelf_cavities and adjustment_fraction is None:
+            raise ValueError('Must provide adjustment_fraction for '
+                             'initializing meshes with ice-shelf cavities')
+        self.adjustment_fraction = adjustment_fraction
 
         package = 'compass.ocean.tests.global_ocean.init'
 
@@ -83,8 +106,10 @@ def __init__(self, test_case, mesh, initial_condition):
         self.add_namelist_options(options, mode='init')
         self.add_streams_file(package, 'streams.topo', mode='init')
 
-        cull_step = self.mesh.steps['cull_mesh']
-        target = os.path.join(cull_step.path, 'topography_culled.nc')
+        if culled_topo_path is None:
+            cull_step = self.mesh.steps['cull_mesh']
+            culled_topo_path = cull_step.path
+        target = os.path.join(culled_topo_path, 'topography_culled.nc')
         self.add_input_file(filename='topography_culled.nc',
                             work_dir_target=target)
 
@@ -168,6 +193,7 @@ def run(self):
         Run this step of the testcase
         """
         config = self.config
+        section = config['global_ocean']
         self._smooth_topography()
 
         interfaces = generate_1d_grid(config=config)
@@ -185,10 +211,14 @@ def run(self):
         namelist = {'config_global_ocean_minimum_depth': f'{min_depth}'}
 
         if self.mesh.with_ice_shelf_cavities:
-            cavity_min_levels = \
-                config.getint('global_ocean', 'cavity_min_levels')
-            cavity_min_layer_thickness = \
-                config.getfloat('global_ocean', 'cavity_min_layer_thickness')
+            frac = self.adjustment_fraction
+            cavity_min_levels = section.getint('cavity_min_levels')
+            min_thick_init = section.getfloat(
+                'cavity_min_layer_thickness_initial')
+            min_thick_final = section.getfloat(
+                'cavity_min_layer_thickness_final')
+            cavity_min_layer_thickness = (
+                (1.0 - frac) * min_thick_init + frac * min_thick_final)
             namelist['config_rx1_min_levels'] = f'{cavity_min_levels}'
             namelist['config_rx1_min_layer_thickness'] = \
                 f'{cavity_min_layer_thickness}'