Merge pull request #807 from xylar/create-sowisc12to30e3r3

New Mesh: SOwISC12to30E3r3

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,7 @@ 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=['WC14', 'WCwISC14'])
 

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

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

@@ -16,7 +16,7 @@
     QUMeshFromConfigStep,
 )
 from compass.ocean.tests.global_ocean.mesh.rrs6to18 import RRS6to18BaseMesh
-from compass.ocean.tests.global_ocean.mesh.so12to60 import SO12to60BaseMesh
+from compass.ocean.tests.global_ocean.mesh.so12to30 import SO12to30BaseMesh
 from compass.ocean.tests.global_ocean.mesh.wc14 import WC14BaseMesh
 from compass.ocean.tests.global_ocean.metadata import (
     get_author_and_email_from_git,
@@ -98,8 +98,8 @@ def __init__(self, test_group, mesh_name, high_res_topography):
             base_mesh_step = ARRM10to60BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['RRS6to18', 'RRSwISC6to18']:
             base_mesh_step = RRS6to18BaseMesh(self, name=name, subdir=subdir)
-        elif mesh_name in ['SO12to60', 'SOwISC12to60']:
-            base_mesh_step = SO12to60BaseMesh(self, name=name, subdir=subdir)
+        elif mesh_name in ['SO12to30', 'SOwISC12to30']:
+            base_mesh_step = SO12to30BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['FRIS01to60', 'FRISwISC01to60']:
             base_mesh_step = FRIS01to60BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['FRIS02to60', 'FRISwISC02to60']:

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

@@ -0,0 +1,70 @@
+import numpy as np
+from geometric_features import read_feature_collection
+from mpas_tools.cime.constants import constants
+from mpas_tools.mesh.creation.signed_distance import (
+    signed_distance_from_geojson,
+)
+
+from compass.mesh import QuasiUniformSphericalMeshStep
+
+
+class SO12to30BaseMesh(QuasiUniformSphericalMeshStep):
+    """
+    A step for creating SO12to30 meshes
+    """
+    def setup(self):
+        """
+        Add some input files
+        """
+
+        self.add_input_file(filename='high_res_region.geojson',
+                            package=self.__module__)
+
+        super().setup()
+
+    def build_cell_width_lat_lon(self):
+        """
+        Create cell width array for this mesh on a regular latitude-longitude
+        grid
+
+        Returns
+        -------
+        cellWidth : numpy.array
+            m x n array of cell width in km
+
+        lon : numpy.array
+            longitude in degrees (length n and between -180 and 180)
+
+        lat : numpy.array
+            longitude in degrees (length m and between -90 and 90)
+        """
+
+        dlon = 0.1
+        dlat = dlon
+        earth_radius = constants['SHR_CONST_REARTH']
+        nlon = int(360. / dlon) + 1
+        nlat = int(180. / dlat) + 1
+        lon = np.linspace(-180., 180., nlon)
+        lat = np.linspace(-90., 90., nlat)
+
+        # start with a uniform 30 km background resolution
+        dx_max = 30.
+        cell_width = dx_max * np.ones((nlat, nlon))
+
+        fc = read_feature_collection('high_res_region.geojson')
+
+        so_signed_distance = signed_distance_from_geojson(fc, lon, lat,
+                                                          earth_radius,
+                                                          max_length=0.25)
+
+        # Equivalent to 20 degrees latitude
+        trans_width = 1600e3
+        trans_start = 500e3
+        dx_min = 12.
+
+        weights = 0.5 * (1 + np.tanh((so_signed_distance - trans_start) /
+                                     trans_width))
+
+        cell_width = dx_min * (1 - weights) + cell_width * weights
+
+        return cell_width, lon, lat

compass/ocean/tests/global_ocean/mesh/so12to60/high_res_region.geojson → compass/ocean/tests/global_ocean/mesh/so12to30/high_res_region.geojson

@@ -4,7 +4,7 @@
     {
       "type": "Feature",
       "properties": {
-        "name": "SO12to60 high res region",
+        "name": "SO12to30 high res region",
         "component": "ocean",
         "object": "region",
         "author": "Xylar Asay-Davis"

compass/ocean/tests/global_ocean/mesh/so12to30/namelist.init

@@ -0,0 +1,8 @@
+config_rx1_inner_iter_count = 20
+config_rx1_horiz_smooth_weight = 1.0
+config_rx1_vert_smooth_weight = 1.0
+config_rx1_slope_weight = 1e-2
+config_rx1_zstar_weight = 1.0
+config_rx1_min_levels = 5
+config_rx1_min_layer_thickness = 2.0
+config_global_ocean_minimum_depth = 10.0