Use new land ice pressure function in isomip_plus

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',

compass/ocean/tests/isomip_plus/initial_state.py

@@ -7,7 +7,11 @@
 from mpas_tools.cime.constants import constants
 from mpas_tools.io import write_netcdf
 
-from compass.ocean.iceshelf import compute_land_ice_pressure_from_draft
+from compass.ocean.iceshelf import (
+    compute_land_ice_density_from_draft,
+    compute_land_ice_pressure_from_draft,
+    compute_land_ice_pressure_from_thickness,
+)
 from compass.ocean.tests.isomip_plus.geom import interpolate_geom
 from compass.ocean.tests.isomip_plus.viz.plot import MoviePlotter
 from compass.ocean.vertical import init_vertical_coord
@@ -125,14 +129,20 @@ def _compute_initial_condition(self):
 
         ds['landIceFraction'] = xr.where(mask, ds.landIceFraction, 0.)
 
-        ref_density = constants['SHR_CONST_RHOSW']
-        landIceDraft = ds.ssh
-        landIcePressure = compute_land_ice_pressure_from_draft(
-            land_ice_draft=landIceDraft, modify_mask=ds.ssh < 0.,
-            ref_density=ref_density)
+        land_ice_draft = ds.ssh
+        if thin_film_present:
+            land_ice_thickness = ds.landIceThickness
+            land_ice_density = compute_land_ice_density_from_draft(
+                land_ice_draft, land_ice_thickness,
+                floating_mask=ds.landIceFloatingMask)
+            land_ice_pressure = compute_land_ice_pressure_from_thickness(
+                land_ice_thickness=land_ice_thickness, modify_mask=ds.ssh < 0.,
+                land_ice_density=land_ice_density)
+        else:
+            land_ice_pressure = compute_land_ice_pressure_from_draft(
+                land_ice_draft=land_ice_draft, modify_mask=land_ice_draft < 0.)
 
-        ds['landIcePressure'] = landIcePressure
-        ds['landIceDraft'] = landIceDraft
+        ds['landIcePressure'] = land_ice_pressure
 
         if self.time_varying_forcing:
             self._write_time_varying_forcing(ds_init=ds)
@@ -173,8 +183,8 @@ def _compute_initial_condition(self):
         # that this effect isn't signicant.
         freezing_temp = (6.22e-2 +
                          -5.63e-2 * init_bot_sal +
-                         -7.43e-8 * landIcePressure +
-                         -1.74e-10 * landIcePressure * init_bot_sal)
+                         -7.43e-8 * land_ice_pressure +
+                         -1.74e-10 * land_ice_pressure * init_bot_sal)
         _, thin_film_temp = np.meshgrid(ds.refZMid, freezing_temp)
         _, thin_film_mask = np.meshgrid(ds.refZMid, thin_film_mask)
         thin_film_temp = np.expand_dims(thin_film_temp, axis=0)
@@ -233,6 +243,8 @@ def _plot(self, ds):
             ds['oceanFracObserved'].expand_dims(dim='Time', axis=0)
         ds['landIcePressure'] = \
             ds['landIcePressure'].expand_dims(dim='Time', axis=0)
+        ds['landIceThickness'] = \
+            ds['landIceThickness'].expand_dims(dim='Time', axis=0)
         ds['landIceGroundedFraction'] = \
             ds['landIceGroundedFraction'].expand_dims(dim='Time', axis=0)
         ds['bottomDepth'] = ds['bottomDepth'].expand_dims(dim='Time', axis=0)
@@ -248,7 +260,10 @@ def _plot(self, ds):
                                   True)
         plotter.plot_horiz_series(ds.landIcePressure,
                                   'landIcePressure', 'landIcePressure',
-                                  True, vmin=1, vmax=1e6, cmap_scale='log')
+                                  True, vmin=1e5, vmax=1e7, cmap_scale='log')
+        plotter.plot_horiz_series(ds.landIceThickness,
+                                  'landIceThickness', 'landIceThickness',
+                                  True, vmin=0, vmax=1e3)
         plotter.plot_horiz_series(ds.ssh,
                                   'ssh', 'ssh',
                                   True, vmin=-700, vmax=0)

compass/ocean/tests/isomip_plus/process_geom.py

@@ -171,10 +171,18 @@ def _smooth_geometry(land_fraction, ds, filter_sigma, threshold=0.01):
                                       mode='constant', cval=0.)
         bed[mask] /= smoothed_mask[mask]
 
+        # We only use the ice surface where the ocean is present to compute
+        # land ice pressure
+        ice_thickness = filters.gaussian_filter(
+            ds.iceThickness * ocean_fraction, filter_sigma, mode='constant',
+            cval=0.)
+        ice_thickness[mask] /= smoothed_mask[mask]
+
         smoothed_draft_mask = filters.gaussian_filter(
             ds.landIceFloatingFraction, filter_sigma, mode='constant', cval=0.)
         smoothed_draft_mask[mask] /= smoothed_mask[mask]
 
         ds['Z_ice_draft'] = (('y', 'x'), draft)
         ds['Z_bed'] = (('y', 'x'), bed)
+        ds['iceThickness'] = (('y', 'x'), ice_thickness)
         ds['smoothedDraftMask'] = (('y', 'x'), smoothed_draft_mask)