Merge pull request #685 from cbegeman/add-ice-shelf-2d-tidal-cases

Add additional ice_shelf_2d test cases

compass/ocean/iceshelf.py

@@ -42,7 +42,7 @@ def compute_land_ice_pressure_and_draft(ssh, modify_mask, ref_density):
 
 
 def adjust_ssh(variable, iteration_count, step, update_pio=True,
-               convert_to_cdf5=False):
+               convert_to_cdf5=False, delta_ssh_threshold=None):
     """
     Adjust the sea surface height or land-ice pressure to be dynamically
     consistent with one another.  A series of short model runs are performed,
@@ -168,6 +168,7 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
                 iCell = numpy.abs(deltaSSH.where(mask)).argmax().values
 
                 ds_cell = ds.isel(nCells=iCell)
+                deltaSSHMax = deltaSSH.isel(nCells=iCell).values
 
                 if on_a_sphere:
                     coords = (f'lon/lat: '
@@ -177,7 +178,7 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
                     coords = (f'x/y: {1e-3 * ds_cell.xCell.values:f} '
                               f'{1e-3 * ds_cell.yCell.values:f}')
                 string = (f'deltaSSHMax: '
-                          f'{deltaSSH.isel(nCells=iCell).values:g}, {coords}')
+                          f'{deltaSSHMax:g}, {coords}')
                 logger.info(f'     {string}')
                 log_file.write(f'{string}\n')
                 string = (f'ssh: {finalSSH.isel(nCells=iCell).values:g}, '
@@ -186,6 +187,10 @@ def adjust_ssh(variable, iteration_count, step, update_pio=True,
                 logger.info(f'     {string}')
                 log_file.write(f'{string}\n')
 
+                if delta_ssh_threshold is not None:
+                    if abs(deltaSSHMax) < delta_ssh_threshold:
+                        break
+
         logger.info("   - Complete\n")
 
     if out_filename is not None:

compass/ocean/tests/ice_shelf_2d/__init__.py

@@ -14,11 +14,14 @@ def __init__(self, mpas_core):
         """
         super().__init__(mpas_core=mpas_core, name='ice_shelf_2d')
 
-        for resolution in ['5km']:
-            for coord_type in ['z-star', 'z-level']:
+        for resolution in [2e3, 5e3]:
+            for coord_type in ['z-star', 'z-level', 'single_layer']:
                 self.add_test_case(
                     Default(test_group=self, resolution=resolution,
                             coord_type=coord_type))
+                self.add_test_case(
+                    Default(test_group=self, resolution=resolution,
+                            coord_type=coord_type, tidal_forcing=True))
                 self.add_test_case(
                     RestartTest(test_group=self, resolution=resolution,
                                 coord_type=coord_type))
@@ -30,25 +33,31 @@ def configure(resolution, coord_type, config):
 
     Parameters
     ----------
-    resolution : str
-        The resolution of the test case
+    resolution : float
+        The resolution of the test case in meters
 
     coord_type : str
         The type of vertical coordinate (``z-star``, ``z-level``, etc.)
 
     config : compass.config.CompassConfigParser
         Configuration options for this test case
     """
-    res_params = {'5km': {'nx': 10, 'ny': 44, 'dc': 5e3}}
+    dx = 50e3  # width of domain in m
+    dy = 220e3  # length of domain in m
+    dc = resolution
+    nx = int(dx / resolution)
+    # ny needs to be even because it is nonperiodic
+    ny = 2 * int(dy / (2. * resolution))
 
-    if resolution not in res_params:
-        raise ValueError('Unsupported resolution {}. Supported values are: '
-                         '{}'.format(resolution, list(res_params)))
-    res_params = res_params[resolution]
-    for param in res_params:
-        config.set('ice_shelf_2d', param, '{}'.format(res_params[param]))
+    config.set('ice_shelf_2d', 'nx', f'{nx}')
+    config.set('ice_shelf_2d', 'ny', f'{ny}')
+    config.set('ice_shelf_2d', 'dc', f'{dc}')
 
     config.set('vertical_grid', 'coord_type', coord_type)
     if coord_type == 'z-level':
         # we need more vertical resolution
         config.set('vertical_grid', 'vert_levels', '100')
+    elif coord_type == 'single_layer':
+        config.set('vertical_grid', 'vert_levels', '1')
+        config.set('vertical_grid', 'coord_type', 'z-level')
+        config.set('vertical_grid', 'partial_cell_type', 'None')

compass/ocean/tests/ice_shelf_2d/default/__init__.py

@@ -15,14 +15,15 @@ class Default(TestCase):
 
     Attributes
     ----------
-    resolution : str
-        The horizontal resolution of the test case
+    resolution : float
+        The horizontal resolution of the test case in m
 
     coord_type : str
         The type of vertical coordinate (``z-star``, ``z-level``, etc.)
     """
 
-    def __init__(self, test_group, resolution, coord_type):
+    def __init__(self, test_group, resolution, coord_type,
+                 tidal_forcing=False):
         """
         Create the test case
 
@@ -38,19 +39,31 @@ def __init__(self, test_group, resolution, coord_type):
             The type of vertical coordinate (``z-star``, ``z-level``, etc.)
         """
         name = 'default'
+        if tidal_forcing:
+            name = 'tidal_forcing'
+            with_frazil = False
+        else:
+            with_frazil = True
         self.resolution = resolution
         self.coord_type = coord_type
-        subdir = '{}/{}/{}'.format(resolution, coord_type, name)
+        if resolution >= 1e3:
+            res_name = f'{resolution / 1e3:g}km'
+        else:
+            res_name = f'{resolution:g}m'
+        subdir = f'{res_name}/{coord_type}/{name}'
         super().__init__(test_group=test_group, name=name,
                          subdir=subdir)
 
         self.add_step(
             InitialState(test_case=self, resolution=resolution))
         self.add_step(
-            SshAdjustment(test_case=self, ntasks=4, openmp_threads=1))
+            SshAdjustment(test_case=self, coord_type=coord_type,
+                          resolution=resolution, ntasks=4, openmp_threads=1,
+                          tidal_forcing=tidal_forcing))
         self.add_step(
             Forward(test_case=self, ntasks=4, openmp_threads=1,
-                    resolution=resolution, with_frazil=True))
+                    coord_type=coord_type, resolution=resolution,
+                    with_frazil=with_frazil, tidal_forcing=tidal_forcing))
         self.add_step(Viz(test_case=self), run_by_default=False)
 
     def configure(self):

compass/ocean/tests/ice_shelf_2d/forward.py

@@ -1,3 +1,5 @@
+import time
+
 from compass.model import run_model
 from compass.step import Step
 
@@ -9,11 +11,16 @@ class Forward(Step):
 
     Attributes
     ----------
-    resolution : str
-        The resolution of the test case
+    resolution : float
+        The resolution of the test case in m
+
+    coord_type: str
+        The coordinate type (e.g., 'z-star', 'single_layer', etc.)
+
     """
-    def __init__(self, test_case, resolution, name='forward', subdir=None,
-                 ntasks=1, min_tasks=None, openmp_threads=1, with_frazil=True):
+    def __init__(self, test_case, resolution, coord_type, name='forward',
+                 subdir=None, ntasks=1, min_tasks=None, openmp_threads=1,
+                 with_frazil=True, tidal_forcing=False):
         """
         Create a new test case
 
@@ -22,8 +29,11 @@ def __init__(self, test_case, resolution, name='forward', subdir=None,
         test_case : compass.TestCase
             The test case this step belongs to
 
-        resolution : str
-            The resolution of the test case
+        coord_type: str
+            The coordinate type (e.g., 'z-star', 'single_layer', etc.)
+
+        resolution : float
+            The resolution of the test case in m
 
         name : str
             the name of the test case
@@ -58,6 +68,13 @@ def __init__(self, test_case, resolution, name='forward', subdir=None,
 
         self.add_namelist_file('compass.ocean.tests.ice_shelf_2d',
                                'namelist.forward')
+        if coord_type == 'single_layer':
+            self.add_namelist_file(
+                'compass.ocean.tests.ice_shelf_2d',
+                'namelist.single_layer.forward_and_ssh_adjust')
+        if tidal_forcing:
+            self.add_namelist_file('compass.ocean.tests.ice_shelf_2d',
+                                   'namelist.tidal_forcing.forward')
         if with_frazil:
             options = {'config_use_frazil_ice_formation': '.true.',
                        'config_frazil_maximum_depth': '2000.0'}
@@ -71,6 +88,9 @@ def __init__(self, test_case, resolution, name='forward', subdir=None,
         self.add_streams_file('compass.ocean.tests.ice_shelf_2d',
                               'streams.forward')
 
+        self.add_input_file(filename='forcing_data.nc',
+                            target=('../initial_state/'
+                                    'init_mode_forcing_data.nc'))
         self.add_input_file(filename='init.nc',
                             target='../ssh_adjustment/adjusted_init.nc')
         self.add_input_file(filename='graph.info',
@@ -87,4 +107,9 @@ def run(self):
         """
         Run this step of the test case
         """
+        config = self.config
+        dt_per_km = config.getfloat('ice_shelf_2d', 'dt_per_km')
+        dt = dt_per_km * self.resolution / 1.e3
+        dt_str = time.strftime('%H:%M:%S', time.gmtime(dt))
+        self.update_namelist_at_runtime({'config_dt': dt_str})
         run_model(self)

compass/ocean/tests/ice_shelf_2d/ice_shelf_2d.cfg

@@ -8,7 +8,7 @@ grid_type = uniform
 vert_levels = 20
 
 # Depth of the bottom of the ocean
-bottom_depth = 2000.0
+bottom_depth = 1050.0
 
 # The type of vertical coordinate (e.g. z-level, z-star)
 coord_type = z-star
@@ -23,20 +23,23 @@ min_pc_fraction = 0.1
 # config options for 2D ice-shelf testcases
 [ice_shelf_2d]
 
-# Vertical thickness of ocean sub-ice cavity
-cavity_thickness = 50.0
+# Time step in seconds as a function of resolution
+dt_per_km = 5
 
-# Vertical thickness of fixed slope
-slope_height = 500.0
+# Vertical thickness of ocean sub-ice cavity at GL
+y1_water_column_thickness = 50.0
 
-# Horizontal width of angled part of the ice
-edge_width = 15.0e3
+# Vertical thickness of water column thickness at y2
+y2_water_column_thickness = 1050.0
 
-# cavity edge in y
+# GL location in y
 y1 = 30.0e3
 
-# shelf edge in y
-y2 = 60.0e3
+# ice shelf inflection point in y
+y2 = 90.0e3
+
+# ice shelf front location in y
+y3 = 90.0e3
 
 # Temperature of the surface in the northern half of the domain
 temperature = 1.0

compass/ocean/tests/ice_shelf_2d/initial_state.py

@@ -1,3 +1,4 @@
+import numpy as np
 import xarray
 from mpas_tools.cime.constants import constants
 from mpas_tools.io import write_netcdf
@@ -71,9 +72,9 @@ def run(self):
         # d variables are total water-column thickness below ice shelf
         y1 = section.getfloat('y1')
         y2 = section.getfloat('y2')
-        y3 = y2 + section.getfloat('edge_width')
-        d1 = section.getfloat('cavity_thickness')
-        d2 = d1 + section.getfloat('slope_height')
+        y3 = section.getfloat('y3')
+        d1 = section.getfloat('y1_water_column_thickness')
+        d2 = section.getfloat('y2_water_column_thickness')
         d3 = bottom_depth
 
         ds = dsMesh.copy()
@@ -130,3 +131,10 @@ def run(self):
         ds['landIceDraft'] = landIceDraft
 
         write_netcdf(ds, 'initial_state.nc')
+
+        # Generate the tidal forcing dataset whether it is used or not
+        ds_forcing = xarray.Dataset()
+        y_max = np.max(ds.yCell.values)
+        ds_forcing['tidalInputMask'] = xarray.where(
+            ds.yCell > (y_max - 0.6 * 5.0e3), 1.0, 0.0)
+        write_netcdf(ds_forcing, 'init_mode_forcing_data.nc')

compass/ocean/tests/ice_shelf_2d/namelist.forward

@@ -1,4 +1,3 @@
-config_dt = '00:05:00'
 config_btr_dt = '00:00:15'
 config_run_duration = '0000_00:10:00'
 config_use_mom_del2 = .true.

compass/ocean/tests/ice_shelf_2d/namelist.single_layer.forward_and_ssh_adjust

@@ -0,0 +1,14 @@
+config_time_integrator='RK4'
+config_bottom_drag_mode = 'explicit'
+config_explicit_bottom_drag_coeff = 3.0e-3
+config_AM_globalStats_enable = .false.
+config_AM_globalStats_compute_on_startup = .false.
+config_AM_globalStats_write_on_startup = .false.
+config_compute_active_tracer_budgets = .false.
+config_disable_thick_vadv = .true.
+config_disable_thick_sflux = .true.
+config_disable_vel_hmix = .true.
+config_disable_vel_vmix = .true.
+config_disable_vel_vadv = .true.
+config_disable_tr_all_tend = .true.
+config_disable_redi_k33 = .true.

compass/ocean/tests/ice_shelf_2d/namelist.tidal_forcing.forward

@@ -0,0 +1,8 @@
+config_use_tidal_forcing = .true.
+config_use_tidal_forcing_tau = 10000
+config_tidal_forcing_model = 'monochromatic'
+config_tidal_forcing_type = 'thickness_source'
+config_tidal_forcing_monochromatic_amp = 1.0
+config_tidal_forcing_monochromatic_period = 10.0
+config_eos_type = 'linear'
+config_land_ice_flux_mode = 'pressure_only'

compass/ocean/tests/ice_shelf_2d/restart_test/__init__.py

@@ -40,20 +40,26 @@ def __init__(self, test_group, resolution, coord_type):
         name = 'restart_test'
         self.resolution = resolution
         self.coord_type = coord_type
-        subdir = '{}/{}/{}'.format(resolution, coord_type, name)
+        if resolution >= 1e3:
+            res_name = f'{resolution / 1e3:g}km'
+        else:
+            res_name = f'{resolution:g}m'
+        subdir = f'{res_name}/{coord_type}/{name}'
         super().__init__(test_group=test_group, name=name,
                          subdir=subdir)
 
         self.add_step(
             InitialState(test_case=self, resolution=resolution))
         self.add_step(
-            SshAdjustment(test_case=self, ntasks=4, openmp_threads=1))
+            SshAdjustment(test_case=self, resolution=resolution,
+                          coord_type=coord_type, ntasks=4,
+                          openmp_threads=1))
 
         for part in ['full', 'restart']:
             name = '{}_run'.format(part)
-            step = Forward(test_case=self, name=name, subdir=name, ntasks=4,
-                           openmp_threads=1, resolution=resolution,
-                           with_frazil=True)
+            step = Forward(test_case=self, name=name, subdir=name,
+                           coord_type=coord_type, ntasks=4, openmp_threads=1,
+                           resolution=resolution, with_frazil=True)
 
             step.add_namelist_file(
                 'compass.ocean.tests.ice_shelf_2d.restart_test',