Use land-ice pressure to determine ice draft

This was previously only the case for thin-film cases.

To accommodate this, we make the minimum ocean column thickness
10 m for non-thin-film cases.  To make the different values for
thin-film and non-thin-film cases more obvious, we use separate
config options for the two.

Author: xylar

polaris/ocean/tasks/isomip_plus/init/init.py

@@ -9,7 +9,6 @@
 from polaris.ocean.ice_shelf import (
     compute_freezing_temperature,
     compute_land_ice_draft_from_pressure,
-    compute_land_ice_pressure_from_draft,
 )
 from polaris.ocean.vertical import init_vertical_coord
 from polaris.ocean.viz import compute_transect, plot_transect
@@ -109,7 +108,14 @@ def _compute_init_topo(self):
         section = config['isomip_plus_init']
         min_levels = section.getint('minimum_levels')
         min_layer_thickness = section.getfloat('min_layer_thickness')
-        min_column_thickness = section.getfloat('min_column_thickness')
+        if thin_film:
+            min_column_thickness = section.getfloat(
+                'min_column_thickness_thin_film'
+            )
+        else:
+            min_column_thickness = section.getfloat(
+                'min_column_thickness_no_thin_film'
+            )
         min_land_ice_fraction = section.getfloat('min_land_ice_fraction')
 
         ds_topo = xr.open_dataset('topo.nc')
@@ -143,31 +149,17 @@ def _compute_init_topo(self):
 
         ds_init['landIceGroundedFraction'] = ds_topo['landIceGroundedFraction']
 
-        if thin_film:
-            # start from landIcePressure and compute landIceDraft
-            ds_init['landIcePressure'] = ds_topo['landIcePressure']
-
-            land_ice_draft = compute_land_ice_draft_from_pressure(
-                land_ice_pressure=ds_topo.landIcePressure,
-                modify_mask=ds_topo.landIcePressure > 0.0,
-            )
-
-            land_ice_draft = np.maximum(
-                ds_topo.bedrockTopography, land_ice_draft
-            )
-
-            ds_init['landIceDraft'] = land_ice_draft
+        # start from landIcePressure and compute landIceDraft
+        ds_init['landIcePressure'] = ds_topo['landIcePressure']
 
-        else:
-            # start form landIceDraft and compute landIcePressure
-            ds_init['landIceDraft'] = ds_topo['landIceDraft']
+        land_ice_draft = compute_land_ice_draft_from_pressure(
+            land_ice_pressure=ds_topo.landIcePressure,
+            modify_mask=ds_topo.landIcePressure > 0.0,
+        )
 
-            land_ice_pressure = compute_land_ice_pressure_from_draft(
-                land_ice_draft=ds_topo.landIceDraft,
-                modify_mask=ds_topo.landIceDraft < 0.0,
-            )
+        land_ice_draft = np.maximum(ds_topo.bedrockTopography, land_ice_draft)
 
-            ds_init['landIcePressure'] = land_ice_pressure
+        ds_init['landIceDraft'] = land_ice_draft
 
         ds_init['ssh'] = ds_init.landIceDraft
 
@@ -303,7 +295,15 @@ def _plot(self):
         Plot several fields from the initial condition
         """
         section = self.config['isomip_plus_init']
-        min_column_thickness = section.getfloat('min_column_thickness')
+        thin_film = self.thin_film
+        if thin_film:
+            min_column_thickness = section.getfloat(
+                'min_column_thickness_thin_film'
+            )
+        else:
+            min_column_thickness = section.getfloat(
+                'min_column_thickness_no_thin_film'
+            )
         min_layer_thickness = section.getfloat('min_layer_thickness')
 
         tol = 1e-10

polaris/ocean/tasks/isomip_plus/isomip_plus_init.cfg

@@ -52,8 +52,10 @@ minimum_levels = 3
 # The minimum allowable layer thickness
 min_layer_thickness = 0.0
 
-# Minimum thickness of the initial ocean column (to prevent 'drying')
-min_column_thickness = 1e-3
+# Minimum thickness of the initial ocean column (to prevent 'drying') without
+# and with a thin film region
+min_column_thickness_no_thin_film = 10.0
+min_column_thickness_thin_film = 1e-3
 
 # Minimum fraction of a cell that contains land ice in order for it to be
 # considered a land-ice cell by MPAS-Ocean (landIceMask == 1).

polaris/ocean/tasks/isomip_plus/isomip_plus_test.py

@@ -131,12 +131,8 @@ def configure(self):
             'polaris.ocean.tasks.isomip_plus', 'isomip_plus_init.cfg'
         )
         vertical_coordinate = self.vertical_coordinate
-        thin_film = self.thin_film
         experiment = self.experiment
 
-        if thin_film:
-            config.set('isomip_plus', 'min_column_thickness', '1e-3')
-
         # for most coordinates, use the config options, which is 36 layers
         levels = None
         if vertical_coordinate == 'single-layer':