Merge pull request #723 from cbegeman/ocn-fix-grounded-land-ice-pressure

Fixup ISOMIP+ land ice pressure:

This PR changes the calculation of land ice pressure for ISOMIP+ thin-film test cases so that it is a function of the land ice thickness rather than land ice draft. The land ice pressure is not perfectly identical between Ocean0 and thin_film_Ocean0 test cases, mostly because in the former case we smooth the draft and in the latter we smooth ice thickness. This PR also fixes a bug wherein the thin film region was culled from the domain.

Author: cbegeman

compass/ocean/iceshelf.py

@@ -11,34 +11,94 @@
 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 * -land_ice_pressure / (ref_density * gravity)
+    land_ice_draft = land_ice_draft_array.transpose(*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)))

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

@@ -1,12 +1,5 @@
-import os
-import shutil
-
-import numpy as np
-import xarray
-
 from compass.model import run_model
 from compass.ocean.tests.isomip_plus.evap import update_evaporation_flux
-from compass.ocean.tests.isomip_plus.viz.plot import MoviePlotter
 from compass.ocean.time import get_time_interval_string
 from compass.step import Step
 
@@ -157,78 +150,6 @@ def run(self):
 
         run_model(self)
 
-        if self.name == 'performance':
-            # plot a few fields
-            plot_folder = f'{self.work_dir}/plots'
-            if os.path.exists(plot_folder):
-                shutil.rmtree(plot_folder)
-
-            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
-            plotter = MoviePlotter(inFolder=self.work_dir,
-                                   streamfunctionFolder=self.work_dir,
-                                   outFolder=plot_folder, sectionY=section_y,
-                                   dsMesh=dsMesh, ds=ds, expt=self.experiment,
-                                   showProgress=show_progress)
-
-            plotter.plot_horiz_series(ds.ssh, 'ssh', 'ssh',
-                                      True, vmin=-700, vmax=0)
-            plotter.plot_3d_field_top_bot_section(
-                ds.temperature, nameInTitle='temperature', prefix='temp',
-                units='C', vmin=-2., vmax=1., cmap='cmo.thermal')
-
-            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,
-                             'land_ice_fluxes.nc'))
-            if 'topDragMagnitude' in dsIce.keys():
-                plotter.plot_horiz_series(
-                    dsIce.topDragMagnitude,
-                    'topDragMagnitude', 'topDragMagnitude', True,
-                    vmin=0 + tol, vmax=np.max(dsIce.topDragMagnitude.values),
-                    cmap_set_under='k')
-            if 'landIceHeatFlux' in dsIce.keys():
-                plotter.plot_horiz_series(
-                    dsIce.landIceHeatFlux,
-                    'landIceHeatFlux', 'landIceHeatFlux', True,
-                    vmin=np.min(dsIce.landIceHeatFlux.values),
-                    vmax=np.max(dsIce.landIceHeatFlux.values))
-            if 'landIceInterfaceTemperature' in dsIce.keys():
-                plotter.plot_horiz_series(
-                    dsIce.landIceInterfaceTemperature,
-                    'landIceInterfaceTemperature',
-                    'landIceInterfaceTemperature',
-                    True,
-                    vmin=np.min(dsIce.landIceInterfaceTemperature.values),
-                    vmax=np.max(dsIce.landIceInterfaceTemperature.values))
-            if 'landIceFreshwaterFlux' in dsIce.keys():
-                plotter.plot_horiz_series(
-                    dsIce.landIceFreshwaterFlux,
-                    'landIceFreshwaterFlux', 'landIceFreshwaterFlux', True,
-                    vmin=0 + tol, vmax=1e-4,
-                    cmap_set_under='k', cmap_scale='log')
-            if 'landIceFraction' in dsIce.keys():
-                plotter.plot_horiz_series(
-                    dsIce.landIceFraction,
-                    'landIceFraction', 'landIceFraction', True,
-                    vmin=0 + tol, vmax=1 - tol,
-                    cmap='cmo.balance',
-                    cmap_set_under='k', cmap_set_over='r')
-            if 'landIceFloatingFraction' in dsIce.keys():
-                plotter.plot_horiz_series(
-                    dsIce.landIceFloatingFraction,
-                    'landIceFloatingFraction', 'landIceFloatingFraction',
-                    True, vmin=0 + tol, vmax=1 - tol,
-                    cmap='cmo.balance', cmap_set_under='k', cmap_set_over='r')
-
         if self.name == 'simulation':
             update_evaporation_flux(
                 in_forcing_file='forcing_data_init.nc',

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