Reorganize SSH adjustment

This merge reorganizes SSH adjustment substantially. It is now
a sequence of alternating init and forward runs, with decreasing
minimum water column thickness.

The sea surface height is adjusted instead of `landIcePressure`
for better consistency with planned coupling to MALI.

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)
+
+                # TODO: Recompute landIceDraft 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}'