Add inactive top cell global ocean test case

compass/ocean/inactive_top_cells.py

@@ -0,0 +1,45 @@
+import os
+
+import xarray
+from mpas_tools.io import write_netcdf
+
+
+def remove_inactive_top_cells_output(in_filename, out_filename=None,
+                                     mesh_filename=None):
+    """
+    Remove inactive top cells from the output netCDF file
+
+    Parameters
+    ----------
+    in_filename : str
+        Filename for the netCDF file to be cropped
+
+    out_filename : str, optional
+        Filename for the netCDF file after cropping.  Tbe default name is the
+        original file name with ``_crop`` appended before the extension
+
+    mesh_filename : str, optional
+        Filename for an MPAS mesh if not included in the file to be cropped
+
+    """
+    if not os.path.exists(in_filename):
+        raise OSError(f'File {in_filename} does not exist.')
+
+    if out_filename is None:
+        basename, ext = os.path.splitext(in_filename)
+        out_filename = f'{basename}_crop{ext}'
+
+    with xarray.open_dataset(in_filename) as ds_in:
+        if mesh_filename is None:
+            ds_mesh = ds_in
+        else:
+            ds_mesh = xarray.open_dataset(mesh_filename)
+        minLevelCell = ds_mesh.minLevelCell
+        minval = minLevelCell.min().values
+        maxval = minLevelCell.max().values
+        if minval != maxval:
+            raise ValueError('Expected minLevelCell to have a constant '
+                             'value for inactive top cell tests')
+        ds_out = ds_in.isel(nVertLevels=slice(minval - 1, None))
+
+    write_netcdf(ds_out, out_filename)

compass/ocean/suites/inactive_top.txt

@@ -0,0 +1,7 @@
+ocean/global_ocean/QU240/mesh
+ocean/global_ocean/QU240/WOA23/init
+ocean/global_ocean/QU240/WOA23/init_inactive_top
+ocean/global_ocean/QU240/WOA23/performance_test
+ocean/global_ocean/QU240/WOA23/performance_test_inactive_top
+ocean/global_ocean/QU240/WOA23/RK4/performance_test
+ocean/global_ocean/QU240/WOA23/RK4/performance_test_inactive_top

compass/ocean/tests/global_ocean/__init__.py

@@ -56,7 +56,8 @@ def __init__(self, mpas_core):
                         include_rk4=True,
                         include_regression=True,
                         include_phc=True,
-                        include_en4_1900=True)
+                        include_en4_1900=True,
+                        with_inactive_top_cells=True)
 
         # for other meshes, we do fewer tests
         self._add_tests(mesh_names=['QU', 'Icos', 'QUwISC', 'IcoswISC'],
@@ -88,7 +89,7 @@ def __init__(self, mpas_core):
     def _add_tests(self, mesh_names, DynamicAdjustment,
                    high_res_topography=True, include_rk4=False,
                    include_regression=False, include_phc=False,
-                   include_en4_1900=False):
+                   include_en4_1900=False, with_inactive_top_cells=False):
         """ Add test cases for the given mesh(es) """
 
         default_ic = 'WOA23'
@@ -100,7 +101,8 @@ def _add_tests(self, mesh_names, DynamicAdjustment,
             self.add_test_case(mesh_test)
 
             init_test = Init(test_group=self, mesh=mesh_test,
-                             initial_condition=default_ic)
+                             initial_condition=default_ic,
+                             with_inactive_top_cells=False)
             self.add_test_case(init_test)
 
             time_integrator = default_time_int
@@ -163,6 +165,21 @@ def _add_tests(self, mesh_names, DynamicAdjustment,
                         test_group=self, mesh=mesh_test, init=init_test,
                         time_integrator=time_integrator))
 
+            if with_inactive_top_cells:
+                init_test = Init(test_group=self, mesh=mesh_test,
+                                 initial_condition=default_ic,
+                                 with_inactive_top_cells=True)
+                self.add_test_case(init_test)
+                self.add_test_case(
+                    PerformanceTest(
+                        test_group=self, mesh=mesh_test, init=init_test,
+                        time_integrator=default_time_int))
+                if include_rk4:
+                    self.add_test_case(
+                        PerformanceTest(
+                            test_group=self, mesh=mesh_test, init=init_test,
+                            time_integrator='RK4'))
+
             initial_conditions = []
             if include_phc:
                 initial_conditions.append('PHC')

compass/ocean/tests/global_ocean/forward.py

@@ -3,6 +3,7 @@
 from importlib.resources import contents
 
 from compass.model import run_model
+from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output
 from compass.ocean.tests.global_ocean.metadata import (
     add_mesh_and_init_metadata,
 )
@@ -155,6 +156,9 @@ def __init__(self, test_case, mesh, time_integrator, init=None,
                 filename='graph.info',
                 work_dir_target=f'{init.path}/initial_state/graph.info')
 
+            if init.with_inactive_top_cells:
+                self.add_output_file(filename='output_crop.nc')
+
         self.add_model_as_input()
 
     def setup(self):
@@ -214,6 +218,12 @@ def run(self):
         update_pio = self.config.getboolean('global_ocean',
                                             'forward_update_pio')
         run_model(self, update_pio=update_pio)
+
+        if self.init.with_inactive_top_cells:
+            remove_inactive_top_cells_output(in_filename='output.nc',
+                                             out_filename='output_crop.nc',
+                                             mesh_filename='init.nc')
+
         add_mesh_and_init_metadata(self.outputs, self.config,
                                    init_filename='init.nc')
 

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

@@ -24,8 +24,14 @@ class Init(TestCase):
     init_subdir : str
         The subdirectory within the test group for all test cases with this
         initial condition
+
+    inactive_top_comp_subdir : str
+        If ``with_inactive_top_cells == True``, the subdirectory equivalent to
+        ``init_subdir`` for test cases without inactive top cells for use for
+        validation between tests with and without inactive top cells
     """
-    def __init__(self, test_group, mesh, initial_condition):
+    def __init__(self, test_group, mesh, initial_condition,
+                 with_inactive_top_cells=False):
         """
         Create the test case
 
@@ -43,17 +49,25 @@ def __init__(self, test_group, mesh, initial_condition):
         name = 'init'
         mesh_name = mesh.mesh_name
         ic_dir = initial_condition
-        self.init_subdir = os.path.join(mesh_name, ic_dir)
+        init_subdir = os.path.join(mesh_name, ic_dir)
+        if with_inactive_top_cells:
+            # Add the name of the subdir without inactive top cells whether or
+            # not is has or will be run
+            self.inactive_top_comp_subdir = init_subdir
+            name = f'{name}_inactive_top'
+        self.init_subdir = init_subdir
         subdir = os.path.join(self.init_subdir, name)
         super().__init__(test_group=test_group, name=name, subdir=subdir)
 
         self.mesh = mesh
         self.initial_condition = initial_condition
+        self.with_inactive_top_cells = with_inactive_top_cells
 
         self.add_step(
             InitialState(
                 test_case=self, mesh=mesh,
-                initial_condition=initial_condition))
+                initial_condition=initial_condition,
+                with_inactive_top_cells=with_inactive_top_cells))
 
         if mesh.with_ice_shelf_cavities:
             self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh))
@@ -84,6 +98,15 @@ def configure(self, config=None):
         config.set('global_ocean', 'init_description',
                    descriptions[initial_condition])
 
+        if self.with_inactive_top_cells:
+            # Since we start at minLevelCell = 2, we need to increase the
+            # number of vertical levels in the cfg file to end up with the
+            # intended number in the initial state
+            vert_levels = config.getint('vertical_grid', 'vert_levels')
+            config.set('vertical_grid', 'vert_levels', f'{vert_levels + 1}',
+                       comment='active vertical levels + inactive_top_cells')
+            config.set('vertical_grid', 'inactive_top_cells', '1')
+
     def validate(self):
         """
         Test cases can override this method to perform validation of variables
@@ -93,6 +116,27 @@ def validate(self):
         compare_variables(test_case=self, variables=variables,
                           filename1='initial_state/initial_state.nc')
 
+        if self.with_inactive_top_cells:
+            # construct the work directory for the other test
+            filename2 = os.path.join(self.base_work_dir, self.mpas_core.name,
+                                     self.test_group.name,
+                                     self.inactive_top_comp_subdir,
+                                     'init/initial_state/initial_state.nc')
+            if os.path.exists(filename2):
+                variables = ['temperature', 'salinity', 'layerThickness',
+                             'normalVelocity']
+                compare_variables(
+                    test_case=self, variables=variables,
+                    filename1='initial_state/initial_state_crop.nc',
+                    filename2=filename2,
+                    quiet=False, check_outputs=False,
+                    skip_if_step_not_run=False)
+
+            else:
+                self.logger.warn('The version of "init" without inactive top '
+                                 'cells was not run.  Skipping\n'
+                                 'validation.')
+
         if self.mesh.with_ice_shelf_cavities:
             variables = ['ssh', 'landIcePressure']
             compare_variables(test_case=self, variables=variables,

compass/ocean/tests/global_ocean/init/init.cfg

@@ -0,0 +1,5 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# no inactive top cells by default
+inactive_top_cells = 0

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

@@ -1,7 +1,11 @@
 import os
 from importlib.resources import contents
 
+import xarray
+from mpas_tools.io import write_netcdf
+
 from compass.model import run_model
+from compass.ocean.inactive_top_cells import remove_inactive_top_cells_output
 from compass.ocean.plot import plot_initial_state, plot_vertical_grid
 from compass.ocean.tests.global_ocean.metadata import (
     add_mesh_and_init_metadata,
@@ -23,7 +27,8 @@ class InitialState(Step):
     initial_condition : {'WOA23', 'PHC', 'EN4_1900'}
         The initial condition dataset to use
     """
-    def __init__(self, test_case, mesh, initial_condition):
+    def __init__(self, test_case, mesh, initial_condition,
+                 with_inactive_top_cells):
         """
         Create the step
 
@@ -44,6 +49,7 @@ def __init__(self, test_case, mesh, initial_condition):
         super().__init__(test_case=test_case, name='initial_state')
         self.mesh = mesh
         self.initial_condition = initial_condition
+        self.with_inactive_top_cells = with_inactive_top_cells
 
         package = 'compass.ocean.tests.global_ocean.init'
 
@@ -134,6 +140,9 @@ def __init__(self, test_case, mesh, initial_condition):
                      'graph.info']:
             self.add_output_file(filename=file)
 
+        if with_inactive_top_cells:
+            self.add_output_file(filename='initial_state_crop.nc')
+
     def setup(self):
         """
         Get resources at setup from config options
@@ -162,6 +171,7 @@ def run(self):
         Run this step of the testcase
         """
         config = self.config
+        logger = self.logger
         interfaces = generate_1d_grid(config=config)
 
         write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc')
@@ -171,6 +181,39 @@ def run(self):
         update_pio = config.getboolean('global_ocean', 'init_update_pio')
         run_model(self, update_pio=update_pio)
 
+        if self.with_inactive_top_cells:
+
+            logger.info("   * Updating minLevelCell for inactive top cells")
+
+            in_filename = 'initial_state.nc'
+            out_filename = in_filename
+
+            with xarray.open_dataset(in_filename) as ds:
+                ds.load()
+
+                # keep the data set with Time for output
+                ds_out = ds
+
+                ds = ds.isel(Time=0)
+
+                offset = config.getint('vertical_grid',
+                                       'inactive_top_cells')
+
+                if 'minLevelCell' in ds:
+                    minLevelCell = ds.minLevelCell + offset
+                    ds_out['minLevelCell'] = minLevelCell
+                else:
+                    logger.info("   - Variable minLevelCell, needed for "
+                                "inactive top cells, is missing from the "
+                                "initial condition")
+
+            write_netcdf(ds_out, out_filename)
+
+            remove_inactive_top_cells_output(
+                in_filename=in_filename, out_filename='initial_state_crop.nc')
+
+            logger.info("   - Complete")
+
         add_mesh_and_init_metadata(self.outputs, config,
                                    init_filename='initial_state.nc')
 

compass/ocean/tests/global_ocean/init/streams.init

@@ -23,6 +23,7 @@
     <var name="bottomDepth"/>
     <var name="bottomDepthObserved"/>
     <var name="oceanFracObserved"/>
+    <var name="minLevelCell"/>
     <var name="maxLevelCell"/>
     <var name="vertCoordMovementWeights"/>
     <var name="edgeMask"/>

compass/ocean/tests/global_ocean/mesh/qu/icos.cfg

@@ -9,3 +9,9 @@ prefix = Icos
 mesh_description = MPAS subdivided icosahedral mesh for E3SM version
                    ${e3sm_version} at ${min_res}-km global resolution with
                    <<<levels>>> vertical level
+
+# Options related to the vertical grid
+[vertical_grid]
+
+# Number of inactive top cell layers
+inactive_top_cells = 0

compass/ocean/tests/global_ocean/mesh/qu/qu.cfg

@@ -20,6 +20,8 @@ max_layer_thickness = 250.0
 # the min and max occurs
 transition_levels = 28
 
+# Number of inactive top cell layers
+inactive_top_cells = 0
 
 # options for global ocean testcases
 [global_ocean]

compass/ocean/tests/global_ocean/mesh/qu240/qu240.cfg

@@ -16,6 +16,9 @@ min_layer_thickness = 3.0
 # The maximum layer thickness
 max_layer_thickness = 500.0
 
+# Number of inactive top cell layers
+inactive_top_cells = 0
+
 
 # options for spherical meshes
 [spherical_mesh]