Switch from restart file to mesh file in all tasks

docs/tutorials/dev_add_task.rst

@@ -997,7 +997,7 @@ script into the ``customize_masked_climatology()`` function:
         """
         logger = self.logger
 
-        ds_mesh = xr.open_dataset(self.restartFileName)
+        ds_mesh = xr.open_dataset(self.meshFilename)
         ds_mesh = ds_mesh[['cellsOnEdge', 'cellsOnVertex', 'nEdgesOnCell',
                            'edgesOnCell', 'verticesOnCell', 'verticesOnEdge',
                            'dcEdge', 'dvEdge']]
@@ -1313,7 +1313,7 @@ described in this tutorial:
             """
             logger = self.logger
 
-            ds_mesh = xr.open_dataset(self.restartFileName)
+            ds_mesh = xr.open_dataset(self.meshFilename)
             ds_mesh = ds_mesh[['cellsOnEdge', 'cellsOnVertex', 'nEdgesOnCell',
                                'edgesOnCell', 'verticesOnCell', 'verticesOnEdge',
                                'dcEdge', 'dvEdge']]

docs/tutorials/dev_understand_a_task.rst

@@ -779,8 +779,8 @@ at that before we continue with ``customize_masked_climatology()``.
             Compute the OHC from the temperature and layer thicknesses in a given
             climatology data sets.
             """
-            ds_restart = xr.open_dataset(self.restartFileName)
-            ds_restart = ds_restart.isel(Time=0)
+            ds_mesh = xr.open_dataset(self.meshFilename)
+            ds_mesh = ds_mesh.isel(Time=0)
 
             # specific heat [J/(kg*degC)]
             cp = self.namelist.getfloat('config_specific_heat_sea_water')
@@ -789,18 +789,18 @@ at that before we continue with ``customize_masked_climatology()``.
 
             units_scale_factor = 1e-9
 
-            n_vert_levels = ds_restart.sizes['nVertLevels']
+            n_vert_levels = ds_mesh.sizes['nVertLevels']
 
-            z_mid = compute_zmid(ds_restart.bottomDepth, ds_restart.maxLevelCell-1,
-                                 ds_restart.layerThickness)
+            z_mid = compute_zmid(ds_mesh.bottomDepth, ds_mesh.maxLevelCell-1,
+                                 ds_mesh.layerThickness)
 
             vert_index = xr.DataArray.from_dict(
                 {'dims': ('nVertLevels',), 'data': np.arange(n_vert_levels)})
 
             temperature = climatology['timeMonthly_avg_activeTracers_temperature']
             layer_thickness = climatology['timeMonthly_avg_layerThickness']
 
-            masks = [vert_index < ds_restart.maxLevelCell,
+            masks = [vert_index < ds_mesh.maxLevelCell,
                      z_mid <= self.min_depth,
                      z_mid >= self.max_depth]
             for mask in masks:
@@ -812,7 +812,7 @@ at that before we continue with ``customize_masked_climatology()``.
             return ohc
 
 This function uses a combination of mesh information taken from an MPAS
-restart file (available from the ``self.restartFileName`` attribute inherited
+mesh file (available from the ``self.meshFilename`` attribute inherited
 from :py:class:`~mpas_analysis.shared.climatology.RemapMpasClimatologySubtask`),
 namelist options available from the ``self.namelist`` reader (inherited from
 :py:class:`~mpas_analysis.shared.AnalysisTask`), and ``temperature`` and
@@ -1160,8 +1160,8 @@ here is the full analysis task as described in this tutorial:
             Compute the OHC from the temperature and layer thicknesses in a given
             climatology data sets.
             """
-            ds_restart = xr.open_dataset(self.restartFileName)
-            ds_restart = ds_restart.isel(Time=0)
+            ds_mesh = xr.open_dataset(self.meshFilename)
+            ds_mesh = ds_mesh.isel(Time=0)
 
             # specific heat [J/(kg*degC)]
             cp = self.namelist.getfloat('config_specific_heat_sea_water')
@@ -1170,18 +1170,18 @@ here is the full analysis task as described in this tutorial:
 
             units_scale_factor = 1e-9
 
-            n_vert_levels = ds_restart.sizes['nVertLevels']
+            n_vert_levels = ds_mesh.sizes['nVertLevels']
 
-            z_mid = compute_zmid(ds_restart.bottomDepth, ds_restart.maxLevelCell-1,
-                                 ds_restart.layerThickness)
+            z_mid = compute_zmid(ds_mesh.bottomDepth, ds_mesh.maxLevelCell-1,
+                                 ds_mesh.layerThickness)
 
             vert_index = xr.DataArray.from_dict(
                 {'dims': ('nVertLevels',), 'data': np.arange(n_vert_levels)})
 
             temperature = climatology['timeMonthly_avg_activeTracers_temperature']
             layer_thickness = climatology['timeMonthly_avg_layerThickness']
 
-            masks = [vert_index < ds_restart.maxLevelCell,
+            masks = [vert_index < ds_mesh.maxLevelCell,
                      z_mid <= self.min_depth,
                      z_mid >= self.max_depth]
             for mask in masks:

mpas_analysis/default.cfg

@@ -527,6 +527,13 @@ lonLines = np.arange(-180., 181., 30.)
 generate = True
 
 
+[ocean]
+## options related to ocean analysis
+
+# the name of the stream that points to the MPAS mesh file.
+meshStream = mesh
+
+
 [oceanObservations]
 ## options related to ocean observations with which the results will be
 ## compared
@@ -585,6 +592,14 @@ remappedClimSubdirectory = clim/obs/remapped
 baseDirectory = /dir/to/ocean/reference
 
 
+[seaIce]
+## options related to sea-ice analysis
+
+# the name of the stream that points to the MPAS mesh file.  The "mesh" stream
+# in MPAS-seaice points to a dummy file, so it is not useful for this purpose.
+meshStream = landIceMasks
+
+
 [seaIceObservations]
 ## options related to sea ice observations with which the results will be
 ## compared

mpas_analysis/ocean/climatology_map_antarctic_melt.py

@@ -333,8 +333,8 @@ def run_task(self):
         # -------
         # Xylar Asay-Davis
 
-        # first, load the land-ice mask from the restart file
-        dsLandIceMask = xr.open_dataset(self.restartFileName)
+        # first, load the land-ice mask from the mesh file
+        dsLandIceMask = xr.open_dataset(self.meshFilename)
         dsLandIceMask = dsLandIceMask[['landIceMask']]
         dsLandIceMask = dsLandIceMask.isel(Time=0)
         self.landIceMask = dsLandIceMask.landIceMask > 0.
@@ -555,12 +555,11 @@ def run_task(self):
                 cellMasks = \
                     dsRegionMask.regionCellMasks.chunk({'nRegions': 10})
 
-                restartFileName = \
-                    self.runStreams.readpath('restart')[0]
+                meshFilename = self.get_mesh_filename()
 
-                dsRestart = xr.open_dataset(restartFileName)
-                landIceFraction = dsRestart.landIceFraction.isel(Time=0)
-                areaCell = dsRestart.areaCell
+                dsMesh = xr.open_dataset(meshFilename)
+                landIceFraction = dsMesh.landIceFraction.isel(Time=0)
+                areaCell = dsMesh.areaCell
 
                 # convert from kg/s to kg/yr
                 totalMeltFlux = constants.sec_per_year * \

mpas_analysis/ocean/climatology_map_bsf.py

@@ -315,7 +315,7 @@ def customize_masked_climatology(self, climatology, season):
         logger = self.logger
         config = self.config
 
-        ds_mesh = xr.open_dataset(self.restartFileName)
+        ds_mesh = xr.open_dataset(self.meshFilename)
         var_list = [
             'cellsOnEdge',
             'cellsOnVertex',
@@ -332,6 +332,10 @@ def customize_masked_climatology(self, climatology, season):
             'latVertex',
             'areaTriangle',
         ]
+        if 'minLevelCell' not in ds_mesh:
+            # some older meshes don't have this one
+            ds_mesh['minLevelCell'] = xr.ones_like(ds_mesh.maxLevelCell)
+
         ds_mesh = ds_mesh[var_list].as_numpy()
 
         masked_filename = self.get_masked_file_name(season)

mpas_analysis/ocean/climatology_map_custom.py

@@ -316,10 +316,10 @@ def _add_thermal_forcing(self, climatology, derivedVars):
         press = dp.cumsum(dim='nVertLevels') - 0.5*dp
 
         # add land ice pressure if available
-        ds_restart = xr.open_dataset(self.restartFileName)
-        ds_restart = ds_restart.isel(Time=0)
-        if 'landIcePressure' in ds_restart:
-            press += ds_restart.landIcePressure
+        ds_mesh = xr.open_dataset(self.meshFilename)
+        ds_mesh = ds_mesh.isel(Time=0)
+        if 'landIcePressure' in ds_mesh:
+            press += ds_mesh.landIcePressure
 
         tempFreeze = c0 + cs*salin + cp*press + cps*press*salin
 

mpas_analysis/ocean/climatology_map_ohc_anomaly.py

@@ -295,8 +295,8 @@ def _compute_ohc(self, climatology):
         Compute the OHC from the temperature and layer thicknesses in a given
         climatology data sets.
         """
-        ds_restart = xr.open_dataset(self.restartFileName)
-        ds_restart = ds_restart.isel(Time=0)
+        ds_mesh = xr.open_dataset(self.meshFilename)
+        ds_mesh = ds_mesh.isel(Time=0)
 
         # specific heat [J/(kg*degC)]
         cp = self.namelist.getfloat('config_specific_heat_sea_water')
@@ -305,18 +305,18 @@ def _compute_ohc(self, climatology):
 
         units_scale_factor = 1e-9
 
-        n_vert_levels = ds_restart.sizes['nVertLevels']
+        n_vert_levels = ds_mesh.sizes['nVertLevels']
 
-        z_mid = compute_zmid(ds_restart.bottomDepth, ds_restart.maxLevelCell-1,
-                             ds_restart.layerThickness)
+        z_mid = compute_zmid(ds_mesh.bottomDepth, ds_mesh.maxLevelCell-1,
+                             ds_mesh.layerThickness)
 
         vert_index = xr.DataArray.from_dict(
             {'dims': ('nVertLevels',), 'data': np.arange(n_vert_levels)})
 
         temperature = climatology['timeMonthly_avg_activeTracers_temperature']
         layer_thickness = climatology['timeMonthly_avg_layerThickness']
 
-        masks = [vert_index < ds_restart.maxLevelCell,
+        masks = [vert_index < ds_mesh.maxLevelCell,
                  z_mid <= self.min_depth,
                  z_mid >= self.max_depth]
         for mask in masks:

mpas_analysis/ocean/climatology_map_wind_stress_curl.py

@@ -174,7 +174,7 @@ def customize_masked_climatology(self, climatology, season):
         """
         logger = self.logger
 
-        ds_mesh = xr.open_dataset(self.restartFileName)
+        ds_mesh = xr.open_dataset(self.meshFilename)
         var_list = [
             'verticesOnEdge',
             'cellsOnVertex',

mpas_analysis/ocean/compute_transects_subtask.py

@@ -247,7 +247,7 @@ def run_task(self):
 
         # first, get maxLevelCell and zMid, needed for masking
 
-        dsMesh = xr.open_dataset(self.restartFileName)
+        dsMesh = xr.open_dataset(self.meshFilename)
         dsMesh = dsMesh.isel(Time=0)
 
         self.maxLevelCell = dsMesh.maxLevelCell - 1

mpas_analysis/ocean/histogram.py

@@ -247,10 +247,10 @@ def run_task(self):
         ds_mask = ds_region_mask.isel(nRegions=region_index)
         cell_mask = ds_mask.regionCellMasks == 1
 
-        # Open the restart file, which contains unmasked weight variables
-        restart_filename = self.runStreams.readpath('restart')[0]
-        ds_restart = xarray.open_dataset(restart_filename)
-        ds_restart = ds_restart.isel(Time=0)
+        # Open the mesh file, which contains unmasked weight variables
+        mesh_filename = self.get_mesh_filename()
+        ds_mesh = xarray.open_dataset(mesh_filename)
+        ds_mesh = ds_mesh.isel(Time=0)
 
         # Save the cell mask only for the region in its own file, which may be
         # referenced by future analysis (i.e., as a control run)
@@ -263,13 +263,15 @@ def run_task(self):
             # Fetch the weight variables and mask them for each region
             for index, var in enumerate(self.variableList):
                 weight_var_name = self.weightList[index]
-                if weight_var_name in ds_restart.keys():
+                if weight_var_name in ds_mesh.keys():
                     var_name = f'timeMonthly_avg_{var}'
                     ds_weights[f'{var_name}_weight'] = \
-                        ds_restart[weight_var_name].where(cell_mask, drop=True)
+                        ds_mesh[weight_var_name].where(cell_mask, drop=True)
                 else:
-                    self.logger.warn(f'Weight variable {weight_var_name} is '
-                                     f'not in the restart file, skipping')
+                    self.logger.warning(
+                        f'Weight variable {weight_var_name} is '
+                        f'not in the mesh file, skipping'
+                    )
 
         weights_filename = \
             f'{base_directory}/{self.filePrefix}_{self.regionName}_weights.nc'

mpas_analysis/ocean/meridional_heat_transport.py

@@ -181,14 +181,10 @@ def run_task(self):
 
             # Read in depth and MHT latitude points
             # Latitude is from binBoundaryMerHeatTrans
-            try:
-                restartFileName = self.runStreams.readpath('restart')[0]
-            except ValueError:
-                raise IOError('No MPAS-O restart file found: need at least '
-                              'one for MHT calcuation')
-
-            with xr.open_dataset(restartFileName) as dsRestart:
-                refBottomDepth = dsRestart.refBottomDepth
+            meshFilename = self.get_mesh_filename()
+
+            with xr.open_dataset(meshFilename) as dsMesh:
+                refBottomDepth = dsMesh.refBottomDepth
 
             nVertLevels = refBottomDepth.sizes['nVertLevels']
             refLayerThickness = np.zeros(nVertLevels)

mpas_analysis/ocean/ocean_regional_profiles.py

@@ -344,22 +344,21 @@ def run_task(self):
             return
 
         # get areaCell
-        restartFileName = \
-            self.runStreams.readpath('restart')[0]
+        meshFilename = self.get_mesh_filename()
 
-        dsRestart = xr.open_dataset(restartFileName)
-        dsRestart = dsRestart.isel(Time=0)
-        areaCell = dsRestart.areaCell
+        dsMesh = xr.open_dataset(meshFilename)
+        dsMesh = dsMesh.isel(Time=0)
+        areaCell = dsMesh.areaCell
 
-        nVertLevels = dsRestart.sizes['nVertLevels']
+        nVertLevels = dsMesh.sizes['nVertLevels']
 
         vertIndex = \
             xr.DataArray.from_dict({'dims': ('nVertLevels',),
                                     'data': np.arange(nVertLevels)})
 
-        vertMask = vertIndex < dsRestart.maxLevelCell
+        vertMask = vertIndex < dsMesh.maxLevelCell
         if self.max_bottom_depth is not None:
-            depthMask = dsRestart.bottomDepth < self.max_bottom_depth
+            depthMask = dsMesh.bottomDepth < self.max_bottom_depth
             vertDepthMask = np.logical_and(vertMask, depthMask)
         else:
             vertDepthMask = vertMask
@@ -436,14 +435,10 @@ def run_task(self):
 
         # Note: restart file, not a mesh file because we need refBottomDepth,
         # not in a mesh file
-        try:
-            restartFile = self.runStreams.readpath('restart')[0]
-        except ValueError:
-            raise IOError('No MPAS-O restart file found: need at least one '
-                          'restart file for plotting time series vs. depth')
+        meshFilename = self.get_mesh_filename()
 
-        with xr.open_dataset(restartFile) as dsRestart:
-            depths = dsRestart.refBottomDepth.values
+        with xr.open_dataset(meshFilename) as dsMesh:
+            depths = dsMesh.refBottomDepth.values
             z = np.zeros(depths.shape)
             z[0] = -0.5 * depths[0]
             z[1:] = -0.5 * (depths[0:-1] + depths[1:])

mpas_analysis/ocean/plot_hovmoller_subtask.py

@@ -259,19 +259,15 @@ def run_task(self):
         ds = ds.set_xindex('regionNames')
         ds = ds.sel(regionNames=self.regionName)
 
-        # Note: restart file, not a mesh file because we need refBottomDepth,
+        # Note: mesh file, not a mesh file because we need refBottomDepth,
         # not in a mesh file
-        try:
-            restartFile = self.runStreams.readpath('restart')[0]
-        except ValueError:
-            raise IOError('No MPAS-O restart file found: need at least one '
-                          'restart file for plotting time series vs. depth')
+        meshFilename = self.get_mesh_filename()
 
         # Define/read in general variables
         self.logger.info('  Read in depth...')
-        with xr.open_dataset(restartFile) as dsRestart:
+        with xr.open_dataset(meshFilename) as dsMesh:
             # reference depth [m]
-            depths = dsRestart.refBottomDepth.values
+            depths = dsMesh.refBottomDepth.values
             z = np.zeros(depths.shape)
             z[0] = -0.5 * depths[0]
             z[1:] = -0.5 * (depths[0:-1] + depths[1:])