Create a thin film below grounded ice for meshes with ice-shelf cavities

compass/ocean/mesh/cull.py

@@ -308,7 +308,7 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
 
     has_remapped_topo = os.path.exists('unsmoothed_topography.nc')
     if with_cavities and not has_remapped_topo:
-        raise ValueError('Mesh culling with caviites must be from '
+        raise ValueError('Mesh culling with cavities must be from '
                          'remapped topography.')
 
     if has_remapped_topo:
@@ -469,7 +469,9 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
         landIceMask = dsMask.regionCellMasks.isel(nRegions=0)
         dsLandIceMask = xr.Dataset()
         dsLandIceMask['landIceMask'] = landIceMask
-        dsLandIceMask['landIceFloatingMask'] = landIceMask
+        dsLandIceMask['landIceFloatingMask'] = xr.logical_and(
+            landIceMask, dsMask.landIceFloatingFraction > 0.
+        )
 
         write_netcdf(dsLandIceMask, 'land_ice_mask.nc')
 

compass/ocean/mesh/remap_topography.py

@@ -129,9 +129,10 @@ def run(self):
         Run this step of the test case
         """
         super().run()
-        self._symlink_unsmoothed()
+        if not self.smoothing:
+            self._symlink_unsmoothed()
         if self.symlinked_to_unsmoothed:
-            # we symlinked to the unsmoothed topography and we're done!
+            print('we symlinked to the unsmoothed topography and we\'re done!')
             return
 
         config = self.config
@@ -362,6 +363,8 @@ def _modify_remapped_bathymetry(self):
         for var in ['bed_elevation', 'landIceThkObserved']:
             ds_out[var] = xr.where(valid, ds_out[var] / norm, 0.)
 
+        # Note: the following variables are not used when MALI topography is
+        # used, they derive from the MALI topography file
         thickness = ds_out.landIceThkObserved
         bed = ds_out.bed_elevation
         flotation_thickness = - (ocean_density / ice_density) * bed

compass/ocean/plot.py

@@ -47,11 +47,14 @@ def plot_initial_state(input_file_name='initial_state.nc',
     varName = 'maxLevelCell'
     var = ds[varName]
     maxLevelCell = var.values - 1
+    if 'nVertLevels' in var.sizes:
+        for iCell in range(var.sizes['nCells']):
+            var[maxLevelCell[iCell]:, iCell] = np.nan
     xarray.plot.hist(var, bins=nVertLevels - 4)
     plt.ylabel('frequency')
     plt.xlabel(varName)
-    txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
-                                          var.max().values, varName)
+    txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, np.nanmin(var.values),
+                                          np.nanmax(var.values), varName)
 
     plt.subplot(4, 3, 3)
     varName = 'bottomDepth'

compass/ocean/tests/global_ocean/__init__.py

@@ -39,6 +39,8 @@ def __init__(self, mpas_core):
 
         # for other meshes, we do fewer tests
         self._add_tests(mesh_names=['QU', 'Icos', 'QUwISC', 'IcoswISC'])
+        self._add_tests(mesh_names=['IcoswISC'],
+                        mali_ais_topo='AIS_4to20km')
 
         self._add_tests(mesh_names=['EC30to60', 'ECwISC30to60'])
 

compass/ocean/tests/global_ocean/init/namelist.wisc

@@ -3,3 +3,4 @@ config_init_vertical_grid_type = 'haney-number'
 config_global_ocean_depress_by_land_ice = .true.
 config_global_ocean_use_constant_land_ice_cavity_temperature = .true.
 config_global_ocean_constant_land_ice_cavity_temperature = -1.8
+config_rx1_max = 20.0

compass/ocean/tests/global_ocean/mesh/__init__.py

@@ -175,7 +175,7 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
                 mesh_name=mesh_name,
                 smoothing=False,
                 mali_ais_topo=mali_ais_topo,
-                ocean_includes_grounded=False)
+                ocean_includes_grounded=True)
 
             self.add_step(unsmoothed_topo)
 
@@ -187,7 +187,7 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
                 smoothing=True,
                 unsmoothed_topo=unsmoothed_topo,
                 mali_ais_topo=mali_ais_topo,
-                ocean_includes_grounded=False)
+                ocean_includes_grounded=True)
 
             self.add_step(smoothed_topo)
 

compass/ocean/tests/global_ocean/mesh/remap_mali_topography/__init__.py

@@ -244,9 +244,9 @@ def _remap_mali_to_mpaso(self, map_filename):
 
         ds_remapped = xr.open_dataset('mali_topography_ncremap.nc')
         ds_out = xr.Dataset()
-        for var_name in ds_remapped:
+        for var_name in ds_remapped.keys():
             var = ds_remapped[var_name]
-            if 'Frac' in var:
+            if 'Frac' in var_name:
                 # copy the fractional variable, making sure it doesn't
                 # exceed 1
                 var = np.minimum(var, 1.)
@@ -315,7 +315,9 @@ def _combine_topo(self):
 
         alpha = ds_out.maliFrac
         # NOTE: MALI's ocean fraction is already scaled by the MALI fraction
-        ds_out['oceanFracObserved'] = (
+        ds_out['oceanFracObserved'] = xr.where(
+            ds_mali.bed_elevation > 0.,
+            0.,
             ds_mali.oceanFrac +
             (1.0 - alpha) * ds_bedmachine.oceanFracObserved)
 

compass/ocean/vertical/grid_1d/__init__.py

@@ -139,7 +139,9 @@ def add_1d_grid(config, ds):
     """
 
     interfaces = generate_1d_grid(config=config)
-
+    plot_1d_grid = config.get('vertical_grid', 'plot_1d_grid')
+    if plot_1d_grid:
+        write_1d_grid(interfaces=interfaces, out_filename='vertical_grid.nc')
     ds['refTopDepth'] = ('nVertLevels', interfaces[0:-1])
     ds['refZMid'] = ('nVertLevels', -0.5 * (interfaces[1:] + interfaces[0:-1]))
     ds['refBottomDepth'] = ('nVertLevels', interfaces[1:])