Add area integrated, annual mean fluxes to DIB and DISMF

Instead of normalizing the fluxes, simply add the totals as a new field to each file so they can be normalized later in code.
MPAS-Dev · Aug 5, 2024 · 6464610 · 6464610
1 parent dc6343f
commit 6464610
Show file tree

Hide file tree

Showing 4 changed files with 69 additions and 85 deletions.
diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py b/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,
 )
@@ -10,9 +13,6 @@
 from compass.ocean.tests.global_ocean.files_for_e3sm.e3sm_to_cmip_maps import (
     E3smToCmipMaps,
 )
-from compass.ocean.tests.global_ocean.files_for_e3sm.normalize_iceberg_ice_shelf_melt import (  # noqa: E501
-    NormalizeIcebergIceShelfMelt,
-)
 from compass.ocean.tests.global_ocean.files_for_e3sm.ocean_graph_partition import (  # noqa: E501
     OceanGraphPartition,
 )
@@ -115,18 +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(RemapIcebergClimatology(test_case=self))
         self.add_step(RemapIceShelfMelt(test_case=self, init=init))
-
-        self.add_step(NormalizeIcebergIceShelfMelt(
-            test_case=self))
-
-        self.add_step(RemapSeaSurfaceSalinityRestoring(
-            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:

diff --git a/..._e3sm/normalize_iceberg_ice_shelf_melt.py → ..._e3sm/add_total_iceberg_ice_shelf_melt.py b/..._e3sm/normalize_iceberg_ice_shelf_melt.py → ..._e3sm/add_total_iceberg_ice_shelf_melt.py
@@ -10,10 +10,11 @@
 )
 
 
-class NormalizeIcebergIceShelfMelt(FilesForE3SMStep):
+class AddTotalIcebergIceShelfMelt(FilesForE3SMStep):
     """
-    A step for for normalizing data iceberg and ice-shelf melt rates on the
-    MPAS grid to a total flux of 1.0 and staging them in ``assembled_files``
+    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):
         """
@@ -24,7 +25,7 @@ def __init__(self, test_case):
         test_case : compass.TestCase
             The test case this step belongs to
         """
-        super().__init__(test_case, name='normalize_iceberg_ice_shelf_melt',
+        super().__init__(test_case, name='add_total_iceberg_ice_shelf_melt',
                          ntasks=1, min_tasks=1)
 
         filename = 'Iceberg_Climatology_Merino_MPAS.nc'
@@ -41,12 +42,14 @@ def __init__(self, test_case):
 
     def setup(self):
         """
-        setup input files based on config options
+        setup output files based on config options
         """
         super().setup()
         if self.with_ice_shelf_cavities:
-            self.add_output_file(filename='dib_merino_2020_normalized.nc')
-            self.add_output_file(filename='dismf_paolo2023_normalized.nc')
+            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):
         """
@@ -59,11 +62,8 @@ def run(self):
 
         logger = self.logger
 
-        suffix = f'{self.mesh_short_name}.{self.creation_date}'
-
         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
@@ -87,39 +87,48 @@ def run(self):
         logger.info(f'total_flux:         {total_flux:.1f}')
         logger.info('')
 
-        for var in ['bergFreshwaterFluxData']:
-            ds_dib[var] = ds_dib[var] / total_flux
+        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'
 
-        write_netcdf(ds_dib, 'dib_merino_2020_normalized.nc')
+        dib_filename = 'Iceberg_Climatology_Merino_MPAS_with_totals.nc'
+        write_netcdf(ds_dib, dib_filename)
 
-        for var in ['dataLandIceFreshwaterFlux', 'dataLandIceHeatFlux']:
-            ds_dismf[var] = ds_dismf[var] / total_flux
-
-        write_netcdf(ds_dismf, 'dismf_paolo2023_normalized.nc')
+        dismf_filename = 'prescribed_ismf_paolo2023_with_totals.nc'
+        write_netcdf(ds_dismf, dismf_filename)
 
         norm_total_dib_flux = (ds_dib.bergFreshwaterFluxData * weights *
-                               area_cell).sum()
+                               area_cell / total_flux).sum()
 
         norm_total_dismf_flux = (ds_dismf.dataLandIceFreshwaterFlux *
-                                 area_cell).sum()
+                                 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:.3f}')
-        logger.info(f'norm_total_dismf_flux:   {norm_total_dismf_flux:.3f}')
-        logger.info(f'norm_total_flux:         {norm_total_flux:.3f}')
+        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_normalized'
+        prefix = 'Iceberg_Climatology_Merino'
+        suffix = f'{self.mesh_short_name}.{self.creation_date}'
         dest_filename = f'{prefix}.{suffix}.nc'
-
         symlink(
-            os.path.abspath('dib_merino_2020_normalized.nc'),
-            f'{self.ocean_inputdata_dir}/{dest_filename}')
+            os.path.abspath(dib_filename),
+            f'{self.seaice_inputdata_dir}/{dest_filename}')
 
-        prefix = 'prescribed_ismf_paolo2023_normalized'
+        prefix = 'prescribed_ismf_paolo2023'
+        suffix = f'{self.mesh_short_name}.{self.creation_date}'
         dest_filename = f'{prefix}.{suffix}.nc'
-
         symlink(
-            os.path.abspath('dismf_paolo2023_normalized.nc'),
+            os.path.abspath(dismf_filename),
             f'{self.ocean_inputdata_dir}/{dest_filename}')
diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/remap_ice_shelf_melt.py b/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,7 +34,7 @@ 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 self.init is not None:
@@ -72,34 +69,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)
diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/remap_iceberg_climatology.py b/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,