Add a step for remapping topography from MALI

compass/ocean/mesh/cull.py

@@ -488,7 +488,7 @@ def _cull_topo(with_cavities, process_count, logger, latitude_threshold,
     write_netcdf(ds_map_culled_to_base, 'map_culled_to_base.nc')
 
     if with_cavities:
-        _add_land_ice_mask_and_mask_draft(ds_topo, ds_base,
+        _flood_fill_and_add_land_ice_mask(ds_topo, ds_base,
                                           ds_map_culled_to_base, ds_preserve,
                                           logger, latitude_threshold,
                                           sweep_count)
@@ -525,7 +525,7 @@ def _land_mask_from_topo(with_cavities, topo_filename, mask_filename):
     write_netcdf(ds_mask, mask_filename)
 
 
-def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
+def _flood_fill_and_add_land_ice_mask(ds_topo, ds_base_mesh,
                                       ds_map_culled_to_base, ds_perserve,
                                       logger, latitude_threshold, sweep_count):
 
@@ -534,7 +534,8 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
         not_preserve = ds.transectCellMasks.sum(dim='nTransects') == 0
 
         for var in ['landIceFracObserved', 'landIcePressureObserved',
-                    'landIceDraftObserved', 'landIceGroundedFracObserved',
+                    'landIceDraftObserved', 'ssh',
+                    'landIceGroundedFracObserved',
                     'landIceFloatingFracObserved', 'landIceThkObserved']:
             ds_topo[var] = ds_topo[var].where(not_preserve, 0.0)
 
@@ -543,22 +544,30 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
 
     land_ice_frac = ds_topo.landIceFracObserved
 
-    # we want the mask to be 1 where there's at least half land-ice
-    land_ice_mask = xr.where(land_ice_frac > 0.5, 1, 0)
+    land_ice_present = xr.where(land_ice_frac > 0.0, 1, 0)
 
     gf = GeometricFeatures()
     fc_ocean_seed = gf.read(componentName='ocean', objectType='point',
                             tags=['seed_point'])
 
     fc_south_pole_seed = read_feature_collection('south_pole.geojson')
 
-    # flood fill the ice portion to remove isolated land ice
+    # flood fill anywhere land ice is present to remove isolated land ice
 
     ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh,
-                                           daGrow=land_ice_mask,
+                                           daGrow=land_ice_present,
                                            fcSeed=fc_south_pole_seed,
                                            logger=logger)
-    land_ice_mask = ds_mask.cellSeedMask
+    land_ice_present = ds_mask.cellSeedMask
+
+    # update land-ice variables and ocean fraction accordingly
+    for var in ['landIceFracObserved', 'landIcePressureObserved',
+                'landIceDraftObserved', 'ssh', 'landIceGroundedFracObserved',
+                'landIceFloatingFracObserved', 'landIceThkObserved']:
+        ds_topo[var] = ds_topo[var].where(land_ice_present, 0.0)
+
+    land_ice_frac = ds_topo.landIceFracObserved
+    land_ice_mask = xr.where(land_ice_frac > 0.5, 1, 0)
 
     # now, remove land-locked or land-ice-locked cells
 
@@ -600,19 +609,6 @@ def _add_land_ice_mask_and_mask_draft(ds_topo, ds_base_mesh,
     region_cell_mask = land_ice_mask.expand_dims(dim='nRegions', axis=1)
     ds_topo['regionCellMasks'] = region_cell_mask
 
-    # we also want to mask out the land-ice draft for cells detatched from the
-    # ice sheet, this time anywhere the land-ice fraction > 0
-    land_ice_draft_mask = xr.where(land_ice_frac > 0.0, 1, 0)
-    ds_mask = compute_mpas_flood_fill_mask(dsMesh=ds_base_mesh,
-                                           daGrow=land_ice_draft_mask,
-                                           fcSeed=fc_south_pole_seed,
-                                           logger=logger)
-
-    land_ice_draft_mask = ds_mask.cellSeedMask
-    ds_topo['landIceDraftMask'] = land_ice_draft_mask
-    ds_topo['landIceDraftObserved'] = (
-        land_ice_draft_mask * ds_topo.landIceDraftObserved)
-
 
 def _land_mask_from_geojson(with_cavities, process_count, logger,
                             mesh_filename, geojson_filename, mask_filename):

compass/ocean/plot.py

@@ -29,7 +29,7 @@ def plot_initial_state(input_file_name='initial_state.nc',
     nVertLevels = ds.sizes['nVertLevels']
 
     fig = plt.figure()
-    fig.set_size_inches(16.0, 12.0)
+    fig.set_size_inches(16.0, 16.0)
     plt.clf()
 
     print('plotting histograms of the initial condition')
@@ -43,7 +43,7 @@ def plot_initial_state(input_file_name='initial_state.nc',
         'number layers: {}\n\n'.format(nVertLevels) + \
         '  min val   max val  variable name\n'
 
-    plt.subplot(3, 3, 2)
+    plt.subplot(4, 3, 2)
     varName = 'maxLevelCell'
     var = ds[varName]
     maxLevelCell = var.values - 1
@@ -53,7 +53,7 @@ def plot_initial_state(input_file_name='initial_state.nc',
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
-    plt.subplot(3, 3, 3)
+    plt.subplot(4, 3, 3)
     varName = 'bottomDepth'
     var = ds[varName]
     xarray.plot.hist(var, bins=nVertLevels - 4)
@@ -75,7 +75,7 @@ def plot_initial_state(input_file_name='initial_state.nc',
     cellMask = xarray.DataArray(data=cellMask, dims=('nCells', 'nVertLevels'))
     edgeMask = xarray.DataArray(data=edgeMask, dims=('nEdges', 'nVertLevels'))
 
-    plt.subplot(3, 3, 4)
+    plt.subplot(4, 3, 4)
     varName = 'temperature'
     var = ds[varName].isel(Time=0).where(cellMask)
     xarray.plot.hist(var, bins=100, log=True)
@@ -84,23 +84,23 @@ def plot_initial_state(input_file_name='initial_state.nc',
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
-    plt.subplot(3, 3, 5)
+    plt.subplot(4, 3, 5)
     varName = 'salinity'
     var = ds[varName].isel(Time=0).where(cellMask)
     xarray.plot.hist(var, bins=100, log=True)
     plt.xlabel(varName)
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
-    plt.subplot(3, 3, 6)
+    plt.subplot(4, 3, 6)
     varName = 'layerThickness'
     var = ds[varName].isel(Time=0).where(cellMask)
     xarray.plot.hist(var, bins=100, log=True)
     plt.xlabel(varName)
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
-    plt.subplot(3, 3, 7)
+    plt.subplot(4, 3, 7)
     varName = 'rx1Edge'
     var = ds[varName].isel(Time=0).where(edgeMask)
     maxRx1Edge = var.max().values
@@ -110,7 +110,7 @@ def plot_initial_state(input_file_name='initial_state.nc',
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
-    plt.subplot(3, 3, 8)
+    plt.subplot(4, 3, 8)
     varName = 'areaCell'
     var = ds[varName]
     xarray.plot.hist(1e-6 * var, bins=100, log=True)
@@ -119,7 +119,7 @@ def plot_initial_state(input_file_name='initial_state.nc',
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
-    plt.subplot(3, 3, 9)
+    plt.subplot(4, 3, 9)
     varName = 'dcEdge'
     var = ds[varName]
     xarray.plot.hist(1e-3 * var, bins=100, log=True)
@@ -128,6 +128,26 @@ def plot_initial_state(input_file_name='initial_state.nc',
     txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
                                           var.max().values, varName)
 
+    plt.subplot(4, 3, 10)
+    var = ds.ssh - (-ds.bottomDepth)
+    if 'landIceMask' in ds:
+        mask = ds.landIceMask == 0
+        var = var.where(mask)
+    xarray.plot.hist(var, bins=100, log=True)
+    plt.ylabel('frequency')
+    plt.xlabel(r'open ocean water-column thickness (m)')
+    txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
+                                          var.max().values, 'open ocean wc')
+
+    if 'landIceMask' in ds:
+        plt.subplot(4, 3, 11)
+        var = (ds.ssh - (-ds.bottomDepth)).where(ds.landIceMask == 1)
+        xarray.plot.hist(var, bins=100, log=True)
+        plt.ylabel('frequency')
+        plt.xlabel(r'ice-shelf cavity water-column thickness (m)')
+        txt = '{}{:9.2e} {:9.2e} {}\n'.format(txt, var.min().values,
+                                              var.max().values, 'cavity wc')
+
     font = FontProperties()
     font.set_family('monospace')
     font.set_size(12)

compass/ocean/suites/so12to30.txt

@@ -0,0 +1,5 @@
+ocean/global_ocean/SO12to30/mesh
+ocean/global_ocean/SO12to30/WOA23/init
+ocean/global_ocean/SO12to30/WOA23/performance_test
+ocean/global_ocean/SO12to30/WOA23/dynamic_adjustment
+ocean/global_ocean/SO12to30/WOA23/files_for_e3sm

compass/ocean/suites/sowisc12to30.txt

@@ -0,0 +1,5 @@
+ocean/global_ocean/SOwISC12to30/mesh
+ocean/global_ocean/SOwISC12to30/WOA23/init
+ocean/global_ocean/SOwISC12to30/WOA23/performance_test
+ocean/global_ocean/SOwISC12to30/WOA23/dynamic_adjustment
+ocean/global_ocean/SOwISC12to30/WOA23/files_for_e3sm

compass/ocean/tests/global_ocean/__init__.py

@@ -44,7 +44,9 @@ def __init__(self, mpas_core):
 
         self._add_tests(mesh_names=['ARRM10to60', 'ARRMwISC10to60'])
 
-        self._add_tests(mesh_names=['SO12to60', 'SOwISC12to60'])
+        self._add_tests(mesh_names=['SO12to30', 'SOwISC12to30'])
+        self._add_tests(mesh_names=['SOwISC12to30'],
+                        mali_ais_topo='AIS_4to20km')
 
         self._add_tests(mesh_names=['WC14', 'WCwISC14'])
 
@@ -65,15 +67,16 @@ def __init__(self, mpas_core):
     def _add_tests(self, mesh_names, high_res_topography=True,
                    include_rk4=False,
                    include_regression=False, include_phc=False,
-                   include_en4_1900=False):
+                   include_en4_1900=False, mali_ais_topo=None):
         """ Add test cases for the given mesh(es) """
 
         default_ic = 'WOA23'
         default_time_int = 'split_explicit_ab2'
 
         for mesh_name in mesh_names:
             mesh_test = Mesh(test_group=self, mesh_name=mesh_name,
-                             high_res_topography=high_res_topography)
+                             high_res_topography=high_res_topography,
+                             mali_ais_topo=mali_ais_topo)
             self.add_test_case(mesh_test)
 
             init_test = Init(test_group=self, mesh=mesh_test,

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

@@ -44,9 +44,8 @@ def __init__(self, test_group, mesh, initial_condition):
             The initial condition dataset to use
         """
         name = 'init'
-        mesh_name = mesh.mesh_name
         ic_dir = initial_condition
-        self.init_subdir = os.path.join(mesh_name, ic_dir)
+        self.init_subdir = os.path.join(mesh.mesh_subdir, ic_dir)
         subdir = os.path.join(self.init_subdir, name)
         super().__init__(test_group=test_group, name=name, subdir=subdir)
 

compass/ocean/tests/global_ocean/init/initial_state.py

@@ -1,6 +1,7 @@
 import os
 from importlib.resources import contents, read_text
 
+import numpy as np
 import xarray as xr
 from jinja2 import Template
 from mpas_tools.io import write_netcdf
@@ -194,7 +195,7 @@ def run(self):
         """
         config = self.config
         section = config['global_ocean']
-        self._smooth_topography()
+        topo_filename = self._smooth_topography()
 
         interfaces = generate_1d_grid(config=config)
 
@@ -223,8 +224,16 @@ def run(self):
             namelist['config_rx1_min_layer_thickness'] = \
                 f'{cavity_min_layer_thickness}'
 
+            min_water_column_thickness = \
+                cavity_min_layer_thickness * cavity_min_levels
+
+            topo_filename = self._dig_cavity_bed_elevation(
+                topo_filename, min_water_column_thickness)
+
         self.update_namelist_at_runtime(namelist)
 
+        symlink(target=topo_filename, link_name='topography.nc')
+
         update_pio = config.getboolean('global_ocean', 'init_update_pio')
         run_model(self, update_pio=update_pio)
 
@@ -250,10 +259,8 @@ def _smooth_topography(self):
         section = config['global_ocean']
         num_passes = section.getint('topo_smooth_num_passes')
         if num_passes == 0:
-            # just symlink the culled topography to be the topography used for
-            # the initial condition
-            symlink(target='topography_culled.nc', link_name='topography.nc')
-            return
+            # just return the culled topography file name
+            return 'topography_culled.nc'
 
         distance_limit = section.getfloat('topo_smooth_distance_limit')
         std_deviation = section.getfloat('topo_smooth_std_deviation')
@@ -274,12 +281,32 @@ def _smooth_topography(self):
         check_call(args=['ocean_smooth_topo_before_init'],
                    logger=self.logger)
 
-        with (xr.open_dataset('topography_culled.nc') as ds_topo):
+        out_filename = 'topography_smoothed.nc'
+        with xr.open_dataset('topography_culled.nc') as ds_topo:
             with xr.open_dataset('topography_orig_and_smooth.nc') as ds_smooth:
                 for field in ['bed_elevation', 'landIceDraftObserved',
                               'landIceThkObserved']:
                     attrs = ds_topo[field].attrs
                     ds_topo[field] = ds_smooth[f'{field}New']
                     ds_topo[field].attrs = attrs
 
-            write_netcdf(ds_topo, 'topography.nc')
+            write_netcdf(ds_topo, out_filename)
+        return out_filename
+
+    def _dig_cavity_bed_elevation(self, in_filename,
+                                  min_water_column_thickness):
+        """ Dig bed elevation to preserve minimum water-column thickness """
+
+        out_filename = 'topography_dig_bed.nc'
+        with xr.open_dataset(in_filename) as ds_topo:
+            bed = ds_topo.bed_elevation
+            attrs = bed.attrs
+            draft = ds_topo.landIceDraftObserved
+            max_bed = draft - min_water_column_thickness
+            mask = np.logical_or(draft == 0., bed < max_bed)
+            bed = xr.where(mask, bed, max_bed)
+            ds_topo['bed_elevation'] = bed
+            ds_topo['bed_elevation'].attrs = attrs
+
+            write_netcdf(ds_topo, out_filename)
+        return out_filename