Fixup ISOMIP+ land ice pressure

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/suites/isomip_rk4.txt

@@ -0,0 +1,11 @@
+# ocean/isomip_plus/planar/5km/z-star/Ocean0
+# ocean/isomip_plus/planar/5km/z-star/time_varying_Ocean0
+# ocean/isomip_plus/planar/5km/z-star/thin_film_Ocean0
+# ocean/isomip_plus/planar/5km/z-star/thin_film_time_varying_Ocean0
+
+#ocean/isomip_plus/planar/5km/sigma/thin_film_Ocean0
+# ocean/isomip_plus/planar/5km/sigma/Ocean0
+# ocean/isomip_plus/planar/5km/sigma/thin_film_time_varying_Ocean0
+ocean/isomip_plus/planar/5km/sigma/thin_film_drying_Ocean0
+ocean/isomip_plus/planar/5km/sigma/thin_film_wetting_Ocean0
+ocean/isomip_plus/planar/5km/sigma/thin_film_tidal_forcing_Ocean0

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)))

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'])