Merge pull request #836 from xylar/add-normalized-dib-dismf

Add area integrated annual mean to data iceberg and ice-shelf flux files

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

@@ -1,6 +1,9 @@
 import os
 
 from compass.io import package_path, symlink
+from compass.ocean.tests.global_ocean.files_for_e3sm.add_total_iceberg_ice_shelf_melt import (  # noqa: E501
+    AddTotalIcebergIceShelfMelt,
+)
 from compass.ocean.tests.global_ocean.files_for_e3sm.diagnostic_maps import (
     DiagnosticMaps,
 )
@@ -112,15 +115,10 @@ def __init__(self, test_group, mesh=None, init=None,
         self.add_step(E3smToCmipMaps(test_case=self))
         self.add_step(DiagnosticMaps(test_case=self))
         self.add_step(DiagnosticMasks(test_case=self))
-
+        self.add_step(RemapIcebergClimatology(test_case=self))
         self.add_step(RemapIceShelfMelt(test_case=self, init=init))
-
-        self.add_step(RemapSeaSurfaceSalinityRestoring(
-            test_case=self))
-
-        self.add_step(RemapIcebergClimatology(
-            test_case=self))
-
+        self.add_step(AddTotalIcebergIceShelfMelt(test_case=self))
+        self.add_step(RemapSeaSurfaceSalinityRestoring(test_case=self))
         self.add_step(RemapTidalMixing(test_case=self))
 
         if mesh is not None and init is not None:

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

@@ -0,0 +1,134 @@
+import os
+
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+
+from compass.io import symlink
+from compass.ocean.tests.global_ocean.files_for_e3sm.files_for_e3sm_step import (  # noqa: E501
+    FilesForE3SMStep,
+)
+
+
+class AddTotalIcebergIceShelfMelt(FilesForE3SMStep):
+    """
+    A step for for adding the total data iceberg and ice-shelf melt rates to
+    to the data iceberg and ice-shelf melt files and staging them in
+    ``assembled_files``
+    """
+    def __init__(self, test_case):
+        """
+        Create a new step
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+        """
+        super().__init__(test_case, name='add_total_iceberg_ice_shelf_melt',
+                         ntasks=1, min_tasks=1)
+
+        filename = 'Iceberg_Climatology_Merino_MPAS.nc'
+        subdir = 'remap_iceberg_climatology'
+        self.add_input_file(
+            filename=filename,
+            target=f'../{subdir}/{filename}')
+
+        filename = 'prescribed_ismf_paolo2023.nc'
+        subdir = 'remap_ice_shelf_melt'
+        self.add_input_file(
+            filename=filename,
+            target=f'../{subdir}/{filename}')
+
+    def setup(self):
+        """
+        setup output files based on config options
+        """
+        super().setup()
+        if self.with_ice_shelf_cavities:
+            self.add_output_file(
+                filename='Iceberg_Climatology_Merino_MPAS_with_totals.nc')
+            self.add_output_file(
+                filename='prescribed_ismf_paolo2023_with_totals.nc')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        super().run()
+
+        if not self.with_ice_shelf_cavities:
+            return
+
+        logger = self.logger
+
+        ds_dib = xr.open_dataset('Iceberg_Climatology_Merino_MPAS.nc')
+        ds_dismf = xr.open_dataset('prescribed_ismf_paolo2023.nc')
+        ds_mesh = xr.open_dataset('restart.nc')
+
+        area_cell = ds_mesh.areaCell
+
+        days_in_month = np.array(
+            [31., 28., 31., 30., 31., 30., 31., 31., 30., 31., 30., 31.])
+
+        weights = xr.DataArray(data=days_in_month / 365.,
+                               dims=('Time',))
+
+        total_dib_flux = (ds_dib.bergFreshwaterFluxData * weights *
+                          area_cell).sum()
+
+        total_dismf_flux = (ds_dismf.dataLandIceFreshwaterFlux *
+                            area_cell).sum()
+
+        total_flux = total_dib_flux + total_dismf_flux
+
+        logger.info(f'total_dib_flux:     {total_dib_flux:.1f}')
+        logger.info(f'total_dismf_flux:   {total_dismf_flux:.1f}')
+        logger.info(f'total_flux:         {total_flux:.1f}')
+        logger.info('')
+
+        for ds in [ds_dib, ds_dismf]:
+            ntime = ds.sizes['Time']
+            field = 'areaIntegAnnMeanDataIcebergFreshwaterFlux'
+            ds[field] = (('Time',), np.ones(ntime) * total_dib_flux.values)
+            ds[field].attrs['units'] = 'kg s-1'
+            field = 'areaIntegAnnMeanDataIceShelfFreshwaterFlux'
+            ds[field] = (('Time',), np.ones(ntime) * total_dismf_flux.values)
+            ds[field].attrs['units'] = 'kg s-1'
+            field = 'areaIntegAnnMeanDataIcebergIceShelfFreshwaterFlux'
+            ds[field] = (('Time',), np.ones(ntime) * total_flux.values)
+            ds[field].attrs['units'] = 'kg s-1'
+
+        dib_filename = 'Iceberg_Climatology_Merino_MPAS_with_totals.nc'
+        write_netcdf(ds_dib, dib_filename)
+
+        dismf_filename = 'prescribed_ismf_paolo2023_with_totals.nc'
+        write_netcdf(ds_dismf, dismf_filename)
+
+        norm_total_dib_flux = (ds_dib.bergFreshwaterFluxData * weights *
+                               area_cell / total_flux).sum()
+
+        norm_total_dismf_flux = (ds_dismf.dataLandIceFreshwaterFlux *
+                                 area_cell / total_flux).sum()
+
+        norm_total_flux = norm_total_dib_flux + norm_total_dismf_flux
+
+        logger.info(f'norm_total_dib_flux:     {norm_total_dib_flux:.16f}')
+        logger.info(f'norm_total_dismf_flux:   {norm_total_dismf_flux:.16f}')
+        logger.info(f'norm_total_flux:         {norm_total_flux:.16f}')
+        logger.info(f'1 - norm_total_flux:     {1 - norm_total_flux:.16g}')
+        logger.info('')
+
+        prefix = 'Iceberg_Climatology_Merino'
+        suffix = f'{self.mesh_short_name}.{self.creation_date}'
+        dest_filename = f'{prefix}.{suffix}.nc'
+        symlink(
+            os.path.abspath(dib_filename),
+            f'{self.seaice_inputdata_dir}/{dest_filename}')
+
+        prefix = 'prescribed_ismf_paolo2023'
+        suffix = f'{self.mesh_short_name}.{self.creation_date}'
+        dest_filename = f'{prefix}.{suffix}.nc'
+        symlink(
+            os.path.abspath(dismf_filename),
+            f'{self.ocean_inputdata_dir}/{dest_filename}')

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

@@ -1,6 +1,3 @@
-import os
-
-from compass.io import symlink
 from compass.ocean.tests.global_ocean.files_for_e3sm.files_for_e3sm_step import (  # noqa: E501
     FilesForE3SMStep,
 )
@@ -37,11 +34,9 @@ def __init__(self, test_case, init):
 
     def setup(self):
         """
-        setup input files based on config options
+        setup input and output files based on config options
         """
         super().setup()
-        if not self.with_ice_shelf_cavities:
-            return
 
         filename = 'prescribed_ismf_paolo2023.nc'
 
@@ -70,34 +65,24 @@ def run(self):
         """
         super().run()
 
-        if not self.with_ice_shelf_cavities:
+        if not self.with_ice_shelf_cavities or self.init is not None:
             return
 
-        prefix = 'prescribed_ismf_paolo2023'
-        suffix = f'{self.mesh_short_name}.{self.creation_date}'
-
-        remapped_filename = f'{prefix}.nc'
-        dest_filename = f'{prefix}.{suffix}.nc'
-
-        if self.init is None:
-            logger = self.logger
-            config = self.config
-            ntasks = self.ntasks
-            in_filename = 'Paolo_2023_ANT_G1920V01_IceShelfMelt.nc'
-
-            parallel_executable = config.get('parallel', 'parallel_executable')
+        logger = self.logger
+        config = self.config
+        ntasks = self.ntasks
+        in_filename = 'Paolo_2023_ANT_G1920V01_IceShelfMelt.nc'
+        remapped_filename = 'prescribed_ismf_paolo2023.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'
+        parallel_executable = config.get('parallel', 'parallel_executable')
 
-            remap_paolo(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)
+        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'
 
-        symlink(
-            os.path.abspath(remapped_filename),
-            f'{self.ocean_inputdata_dir}/{dest_filename}')
+        remap_paolo(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/files_for_e3sm/remap_iceberg_climatology.py

@@ -5,7 +5,6 @@
 from mpas_tools.io import write_netcdf
 from pyremap import LatLonGridDescriptor, MpasCellMeshDescriptor, Remapper
 
-from compass.io import symlink
 from compass.ocean.tests.global_ocean.files_for_e3sm.files_for_e3sm_step import (  # noqa: E501
     FilesForE3SMStep,
 )
@@ -34,7 +33,13 @@ def __init__(self, test_case):
             target='Iceberg_Interannual_Merino.nc',
             database='initial_condition_database')
 
-        self.add_output_file(filename='Iceberg_Climatology_Merino_MPAS.nc')
+    def setup(self):
+        """
+        setup output files based on config options
+        """
+        super().setup()
+        if self.with_ice_shelf_cavities:
+            self.add_output_file(filename='Iceberg_Climatology_Merino_MPAS.nc')
 
     def run(self):
         """
@@ -48,12 +53,7 @@ def run(self):
         ntasks = self.ntasks
 
         in_filename = 'Iceberg_Interannual_Merino.nc'
-
-        prefix = 'Iceberg_Climatology_Merino'
-        suffix = f'{self.mesh_short_name}.{self.creation_date}'
-
-        remapped_filename = f'{prefix}_MPAS.nc'
-        dest_filename = f'{prefix}.{suffix}.nc'
+        remapped_filename = 'Iceberg_Climatology_Merino_MPAS.nc'
 
         parallel_executable = config.get('parallel', 'parallel_executable')
 
@@ -69,10 +69,6 @@ def run(self):
                             logger=logger, mpi_tasks=ntasks,
                             parallel_executable=parallel_executable)
 
-        symlink(
-            os.path.abspath(remapped_filename),
-            f'{self.seaice_inputdata_dir}/{dest_filename}')
-
 
 def remap_iceberg_climo(in_filename, mesh_filename, mesh_name,
                         land_ice_mask_filename, out_filename, logger,

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

@@ -261,10 +261,10 @@ def remap_paolo(in_filename, base_mesh_filename, culled_mesh_filename,
 
     field = 'dataLandIceFreshwaterFlux'
     ds_remap[field] = area_ratio * sphere_fwf
-    ds_remap[field].attrs['units'] = 'kg m^-2 s^-1'
+    ds_remap[field].attrs['units'] = 'kg m-2 s-1'
     field = 'dataLandIceHeatFlux'
     ds_remap[field] = area_ratio * ds_remap[field]
-    ds_remap[field].attrs['units'] = 'W m^-2'
+    ds_remap[field].attrs['units'] = 'W m-2'
 
     mpas_flux = (ds_remap.dataLandIceFreshwaterFlux *
                  mpas_area_cell).sum().values