Remove landIceDraft as a forcing variable from ISOMIP+

compass/ocean/iceshelf.py

@@ -11,34 +11,96 @@
 from compass.model import partition, run_model
 
 
-def compute_land_ice_pressure_and_draft(ssh, modify_mask, ref_density):
+def compute_land_ice_pressure_from_thickness(land_ice_thickness, modify_mask,
+                                             land_ice_density=None):
     """
-    Compute the pressure from and overlying ice shelf and the ice-shelf draft
+    Compute the pressure from an overlying ice shelf from ice thickness
 
     Parameters
     ----------
-    ssh : xarray.DataArray
-        The sea surface height (the ice draft)
+    land_ice_thickness: xarray.DataArray
+        The ice thickness
 
     modify_mask : xarray.DataArray
         A mask that is 1 where ``landIcePressure`` can be deviate from 0
 
-    ref_density : float
+    land_ice_density : float, optional
+        A reference density for land ice
+
+    Returns
+    -------
+    land_ice_pressure : xarray.DataArray
+        The pressure from the overlying land ice on the ocean
+    """
+    gravity = constants['SHR_CONST_G']
+    if land_ice_density is None:
+        land_ice_density = constants['SHR_CONST_RHOICE']
+    land_ice_pressure = modify_mask * \
+        numpy.maximum(land_ice_density * gravity * land_ice_thickness, 0.)
+    return land_ice_pressure
+
+
+def compute_land_ice_pressure_from_draft(land_ice_draft, modify_mask,
+                                         ref_density=None):
+    """
+    Compute the pressure from an overlying ice shelf from ice draft
+
+    Parameters
+    ----------
+    land_ice_draft : xarray.DataArray
+        The ice draft (sea surface height)
+
+    modify_mask : xarray.DataArray
+        A mask that is 1 where ``landIcePressure`` can be deviate from 0
+
+    ref_density : float, optional
         A reference density for seawater displaced by the ice shelf
 
     Returns
     -------
-    landIcePressure : xarray.DataArray
+    land_ice_pressure : xarray.DataArray
         The pressure from the overlying land ice on the ocean
+    """
+    gravity = constants['SHR_CONST_G']
+    if ref_density is None:
+        ref_density = constants['SHR_CONST_RHOSW']
+    land_ice_pressure = \
+        modify_mask * numpy.maximum(-ref_density * gravity * land_ice_draft,
+                                    0.)
+    return land_ice_pressure
+
+
+def compute_land_ice_draft_from_pressure(land_ice_pressure, modify_mask,
+                                         ref_density=None):
+    """
+    Compute the ice-shelf draft associated with the pressure from an overlying
+    ice shelf
+
+    Parameters
+    ----------
+    land_ice_pressure : xarray.DataArray
+        The pressure from the overlying land ice on the ocean
+
+    modify_mask : xarray.DataArray
+        A mask that is 1 where ``landIcePressure`` can be deviate from 0
 
-    landIceDraft : xarray.DataArray
-        The ice draft, equal to the initial ``ssh``
+    ref_density : float, optional
+        A reference density for seawater displaced by the ice shelf
+
+    Returns
+    -------
+    land_ice_draft : xarray.DataArray
+        The ice draft
     """
     gravity = constants['SHR_CONST_G']
-    landIcePressure = \
-        modify_mask * numpy.maximum(-ref_density * gravity * ssh, 0.)
-    landIceDraft = ssh
-    return landIcePressure, landIceDraft
+    if ref_density is None:
+        ref_density = constants['SHR_CONST_RHOSW']
+    land_ice_draft_array = \
+        - (modify_mask.values *
+           land_ice_pressure.values / (ref_density * gravity))
+    land_ice_draft = xarray.DataArray(data=land_ice_draft_array,
+                                      dims=(land_ice_pressure.dims))
+    return land_ice_draft
 
 
 def adjust_ssh(variable, iteration_count, step, update_pio=True,

compass/ocean/tests/ice_shelf_2d/initial_state.py

@@ -5,7 +5,7 @@
 from mpas_tools.mesh.conversion import convert, cull
 from mpas_tools.planar_hex import make_planar_hex_mesh
 
-from compass.ocean.iceshelf import compute_land_ice_pressure_and_draft
+from compass.ocean.iceshelf import compute_land_ice_pressure_from_draft
 from compass.ocean.vertical import init_vertical_coord
 from compass.step import Step
 
@@ -103,8 +103,10 @@ def run(self):
         landIceFloatingMask = landIceMask.copy()
 
         ref_density = constants['SHR_CONST_RHOSW']
-        landIcePressure, landIceDraft = compute_land_ice_pressure_and_draft(
-            ssh=ds.ssh, modify_mask=modify_mask, ref_density=ref_density)
+        landIceDraft = ds.ssh
+        landIcePressure = compute_land_ice_pressure_from_draft(
+            land_ice_draft=landIceDraft, modify_mask=modify_mask,
+            ref_density=ref_density)
 
         salinity = surface_salinity + ((bottom_salinity - surface_salinity) *
                                        (ds.zMid / (-bottom_depth)))
@@ -128,7 +130,6 @@ def run(self):
         ds['landIceMask'] = landIceMask
         ds['landIceFloatingMask'] = landIceFloatingMask
         ds['landIcePressure'] = landIcePressure
-        ds['landIceDraft'] = landIceDraft
 
         write_netcdf(ds, 'initial_state.nc')
 

compass/ocean/tests/isomip_plus/__init__.py

@@ -36,7 +36,7 @@ def __init__(self, mpas_core):
                             experiment=experiment,
                             vertical_coordinate=vertical_coordinate,
                             time_varying_forcing=True))
-                for vertical_coordinate in ['z-star', 'sigma', 'single_layer']:
+                for vertical_coordinate in ['z-star', 'sigma']:
                     self.add_test_case(
                         IsomipPlusTest(
                             test_group=self, resolution=resolution,
@@ -50,7 +50,7 @@ def __init__(self, mpas_core):
                             vertical_coordinate=vertical_coordinate,
                             time_varying_forcing=True,
                             thin_film_present=True))
-                for vertical_coordinate in ['sigma', 'single_layer']:
+                for vertical_coordinate in ['z-star', 'sigma']:
                     self.add_test_case(
                         IsomipPlusTest(
                             test_group=self, resolution=resolution,
@@ -67,22 +67,7 @@ def __init__(self, mpas_core):
                             time_varying_forcing=True,
                             time_varying_load='decreasing',
                             thin_film_present=True))
-                for vertical_coordinate in ['sigma']:
-                    self.add_test_case(
-                        IsomipPlusTest(
-                            test_group=self, resolution=resolution,
-                            experiment=experiment,
-                            vertical_coordinate=vertical_coordinate,
-                            tidal_forcing=True,
-                            thin_film_present=True))
-                for vertical_coordinate in ['single_layer']:
-                    self.add_test_case(
-                        IsomipPlusTest(
-                            test_group=self, resolution=resolution,
-                            experiment=experiment,
-                            vertical_coordinate=vertical_coordinate,
-                            tidal_forcing=True,
-                            thin_film_present=False))
+                for vertical_coordinate in ['sigma', 'z-star']:
                     self.add_test_case(
                         IsomipPlusTest(
                             test_group=self, resolution=resolution,

compass/ocean/tests/isomip_plus/forward.py

@@ -162,11 +162,12 @@ def run(self):
             plot_folder = f'{self.work_dir}/plots'
             if os.path.exists(plot_folder):
                 shutil.rmtree(plot_folder)
+            min_column_thickness = self.config.getfloat('isomip_plus',
+                                                        'min_column_thickness')
 
             dsMesh = xarray.open_dataset(os.path.join(self.work_dir,
                                                       'init.nc'))
             ds = xarray.open_dataset(os.path.join(self.work_dir, 'output.nc'))
-
             section_y = self.config.getfloat('isomip_plus_viz', 'section_y')
             # show progress only if we're not writing to a log file
             show_progress = self.log_filename is None
@@ -176,6 +177,11 @@ def run(self):
                                    dsMesh=dsMesh, ds=ds, expt=self.experiment,
                                    showProgress=show_progress)
 
+            bottomDepth = ds.bottomDepth.expand_dims(dim='Time', axis=0)
+            plotter.plot_horiz_series(ds.ssh + bottomDepth,
+                                      'H', 'H', True,
+                                      vmin=min_column_thickness, vmax=700,
+                                      cmap_set_under='r', cmap_scale='log')
             plotter.plot_horiz_series(ds.ssh, 'ssh', 'ssh',
                                       True, vmin=-700, vmax=0)
             plotter.plot_3d_field_top_bot_section(
@@ -185,6 +191,7 @@ def run(self):
             plotter.plot_3d_field_top_bot_section(
                 ds.salinity, nameInTitle='salinity', prefix='salin',
                 units='PSU', vmin=33.8, vmax=34.7, cmap='cmo.haline')
+
             tol = 1e-10
             dsIce = xarray.open_dataset(
                 os.path.join(self.work_dir,

compass/ocean/tests/isomip_plus/geom.py

@@ -111,11 +111,12 @@ def interpolate_geom(ds_mesh, ds_geom, min_ocean_fraction, thin_film_present):
                                   ds_geom.landIceGroundedFraction)
 
     # mask the topography to the ocean region before interpolation
-    for var in ['Z_bed', 'Z_ice_draft', 'landIceFraction',
+    for var in ['Z_bed', 'Z_ice_surface', 'Z_ice_draft', 'landIceFraction',
                 'landIceFloatingFraction', 'smoothedDraftMask']:
         ds_geom[var] = ds_geom[var] * ds_geom['oceanFraction']
 
     fields = {'bottomDepthObserved': 'Z_bed',
+              'landIceThickness': 'iceThickness',
               'ssh': 'Z_ice_draft',
               'oceanFracObserved': 'oceanFraction',
               'landIceFraction': 'landIceFraction',
@@ -156,7 +157,7 @@ def _get_geom_fields(ds_geom, ds_mesh, thin_film_present):
     y_cell = ds_mesh.yIsomipCell.values
 
     if thin_film_present:
-        ocean_fraction = - ds_geom['landFraction'] + 1.0
+        ocean_fraction = xarray.where(ds_geom['Z_bed'] > 0., 0., 1.)
     else:
         ocean_fraction = (ds_geom['landIceFloatingFraction'] +
                           ds_geom['openOceanFraction'])