Merge pull request #777 from xylar/reroute-dismf

Reroute missing fluxes in data ice-shelf melt fluxes (DISMF)

Author: xylar

compass/ocean/mesh/cull.py

@@ -125,6 +125,7 @@ def setup(self):
             self.add_input_file(filename='topography.nc',
                                 work_dir_target=target)
             self.add_output_file('topography_culled.nc')
+            self.add_output_file('map_culled_to_base.nc')
 
         config = self.config
         self.cpus_per_task = config.getint('spherical_mesh',
@@ -267,6 +268,39 @@ def cull_mesh(with_cavities=False, with_critical_passages=False,
             convert_to_cdf5, latitude_threshold, sweep_count)
 
 
+def write_map_culled_to_base(base_mesh_filename, culled_mesh_filename,
+                             out_filename):
+    ds_base = xr.open_dataset(base_mesh_filename)
+    ncells_base = ds_base.sizes['nCells']
+    lon_base = ds_base.lonCell.values
+    lat_base = ds_base.latCell.values
+
+    ds_culled = xr.open_dataset(culled_mesh_filename)
+    ncells_culled = ds_culled.sizes['nCells']
+    lon_culled = ds_culled.lonCell.values
+    lat_culled = ds_culled.latCell.values
+
+    # create a map from lat-lon pairs to base-mesh cell indices
+    map_base = dict()
+    for cell_index in range(ncells_base):
+        lon = lon_base[cell_index]
+        lat = lat_base[cell_index]
+        map_base[(lon, lat)] = cell_index
+
+    # create a map from culled- to base-mesh cell indices
+    map_culled_to_base = list()
+    for cell_index in range(ncells_culled):
+        lon = lon_culled[cell_index]
+        lat = lat_culled[cell_index]
+        # each (lon, lat) for a culled cell *must* be in the base mesh
+        map_culled_to_base.append(map_base[(lon, lat)])
+
+    ds_map_culled_to_base = xr.Dataset()
+    ds_map_culled_to_base['mapCulledToBase'] = (('nCells',),
+                                                map_culled_to_base)
+    write_netcdf(ds_map_culled_to_base, out_filename)
+
+
 def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
                             custom_critical_passages, custom_land_blockages,
                             preserve_floodplain, use_progress_bar,
@@ -425,6 +459,9 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
         check_call(args, logger=logger)
 
     if has_remapped_topo:
+        write_map_culled_to_base(base_mesh_filename='base_mesh.nc',
+                                 culled_mesh_filename='culeed_mesh.nc',
+                                 out_filename='map_culled_to_base.nc')
         _cull_topo()
 
     if with_cavities:
@@ -471,30 +508,8 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
 
 
 def _cull_topo():
-    ds_base = xr.open_dataset('base_mesh.nc')
-    ncells_base = ds_base.sizes['nCells']
-    lon_base = ds_base.lonCell.values
-    lat_base = ds_base.latCell.values
-
-    ds_culled = xr.open_dataset('culled_mesh.nc')
-    ncells_culled = ds_culled.sizes['nCells']
-    lon_culled = ds_culled.lonCell.values
-    lat_culled = ds_culled.latCell.values
-
-    # create a map from lat-lon pairs to base-mesh cell indices
-    map_base = dict()
-    for cell_index in range(ncells_base):
-        lon = lon_base[cell_index]
-        lat = lat_base[cell_index]
-        map_base[(lon, lat)] = cell_index
-
-    # create a map from culled- to base-mesh cell indices
-    map_culled_to_base = list()
-    for cell_index in range(ncells_culled):
-        lon = lon_culled[cell_index]
-        lat = lat_culled[cell_index]
-        # each (lon, lat) for a culled cell *must* be in the base mesh
-        map_culled_to_base.append(map_base[(lon, lat)])
+    ds_map = xr.open_dataset('map_culled_to_base.nc')
+    map_culled_to_base = ds_map.mapCulledToBase.values
 
     ds_topo = xr.open_dataset('topography.nc')
     ds_topo = ds_topo.isel(nCells=map_culled_to_base)

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

@@ -154,6 +154,13 @@ def configure(self):
             graph_filename = os.path.abspath(graph_filename)
             config.set('files_for_e3sm', 'graph_filename', graph_filename)
 
+            base_mesh_filename = os.path.join(
+                self.base_work_dir, mesh.steps['base_mesh'].path,
+                'base_mesh.nc')
+            base_mesh_filename = os.path.abspath(base_mesh_filename)
+            config.set('files_for_e3sm', 'base_mesh_filename',
+                       base_mesh_filename)
+
         if init is not None:
             if mesh.with_ice_shelf_cavities:
                 initial_state_filename = \

compass/ocean/tests/global_ocean/files_for_e3sm/files_for_e3sm_step.py

@@ -109,25 +109,20 @@ def setup(self):
         """
         setup input files based on config options
         """
+        config = self.config
         self.add_input_file(filename='README', target='../README')
 
-        initial_state_filename = self.config.get(
-            'files_for_e3sm', 'ocean_initial_state_filename')
-        if initial_state_filename != 'autodetect':
-            initial_state_filename = os.path.normpath(os.path.join(
-                self.test_case.work_dir, initial_state_filename))
-            self.add_input_file(filename='initial_state.nc',
-                                target=initial_state_filename)
-
-        restart_filename = self.config.get('files_for_e3sm',
-                                           'ocean_restart_filename')
-        if restart_filename != 'autodetect':
-            restart_filename = os.path.normpath(os.path.join(
-                self.test_case.work_dir, restart_filename))
-            self.add_input_file(filename='restart.nc', target=restart_filename)
-
-        with_ice_shelf_cavities = self.config.get('files_for_e3sm',
-                                                  'with_ice_shelf_cavities')
+        for prefix in ['base_mesh', 'initial_state', 'restart']:
+            filename = config.get('files_for_e3sm',
+                                  f'ocean_{prefix}_filename')
+            if filename != 'autodetect':
+                filename = os.path.normpath(os.path.join(
+                    self.test_case.work_dir, filename))
+                self.add_input_file(filename='base_mesh.nc',
+                                    target=filename)
+
+        with_ice_shelf_cavities = config.get('files_for_e3sm',
+                                             'with_ice_shelf_cavities')
         if with_ice_shelf_cavities != 'autodetect':
             self.with_ice_shelf_cavities = \
                 (with_ice_shelf_cavities.lower() == 'true')
@@ -138,7 +133,7 @@ def run(self):  # noqa: C901
         Run this step of the testcase
         """
         config = self.config
-        for prefix in ['initial_state', 'restart']:
+        for prefix in ['base_mesh', 'initial_state', 'restart']:
             if not os.path.exists(f'{prefix}.nc'):
                 filename = config.get('files_for_e3sm',
                                       f'ocean_{prefix}_filename')

compass/ocean/tests/global_ocean/files_for_e3sm/ocean_mesh.py

@@ -120,3 +120,8 @@ def run(self):
 
         symlink(os.path.abspath(dest_filename),
                 f'{self.ocean_mesh_dir}/{dest_filename}')
+
+        dest_filename = f'{self.mesh_short_name}_base.{self.creation_date}.nc'
+
+        symlink(os.path.abspath('base_mesh.nc'),
+                f'{self.ocean_mesh_dir}/{dest_filename}')

compass/ocean/tests/global_ocean/files_for_e3sm/remap_ice_shelf_melt.py

@@ -81,11 +81,13 @@ def run(self):
 
             parallel_executable = config.get('parallel', 'parallel_executable')
 
-            mesh_filename = 'initial_state.nc'
+            base_mesh_filename = 'base_mesh.nc'
+            culled_mesh_filename = 'initial_state.nc'
             mesh_name = self.mesh_short_name
             land_ice_mask_filename = 'initial_state.nc'
 
-            remap_adusumilli(in_filename, mesh_filename, mesh_name,
+            remap_adusumilli(in_filename, base_mesh_filename,
+                             culled_mesh_filename, mesh_name,
                              land_ice_mask_filename, remapped_filename,
                              logger=logger, mpi_tasks=ntasks,
                              parallel_executable=parallel_executable)

compass/ocean/tests/global_ocean/global_ocean.cfg

@@ -181,6 +181,10 @@ mesh_short_name = autodetect
 # directory of an ocean restart file on the given mesh
 ocean_restart_filename = autodetect
 
+# the base mesh before culling for remapping and rerouting data ice-shelf melt
+# fluxes
+ocean_base_mesh_filename = autodetect
+
 # the initial state used to extract the ocean and sea-ice meshes
 ocean_initial_state_filename = ${ocean_restart_filename}