From e0c3149e87cf2f7b4d812370086c763f3fb9a816 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 5 Aug 2024 14:26:52 -0500 Subject: [PATCH 1/3] Fix units for CF compliance --- compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py b/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py index 4980b81acf..fc010135d3 100644 --- a/compass/ocean/tests/global_ocean/init/remap_ice_shelf_melt.py +++ b/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 From dd199861da814998881a6cd91d8ae502bf4f1ce9 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 1 Jul 2024 08:59:27 -0500 Subject: [PATCH 2/3] Add a step for normalizing data iceberg and ice-shelf fluxes This is needed for an Antarctic balance approach in which all Antarctic solid runoff is converted into iceberg and ice-shelf melt with the patterns from the Merino et al. (2020) and Paolo et al. (2023) datasets. --- .../global_ocean/files_for_e3sm/__init__.py | 10 +- .../normalize_iceberg_ice_shelf_melt.py | 125 ++++++++++++++++++ .../files_for_e3sm/remap_ice_shelf_melt.py | 4 +- 3 files changed, 136 insertions(+), 3 deletions(-) create mode 100644 compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py b/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py index afa6ead955..fd220a1d16 100644 --- a/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py +++ b/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py @@ -10,6 +10,9 @@ 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, ) @@ -113,12 +116,15 @@ def __init__(self, test_group, mesh=None, init=None, 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( + self.add_step(NormalizeIcebergIceShelfMelt( test_case=self)) - self.add_step(RemapIcebergClimatology( + self.add_step(RemapSeaSurfaceSalinityRestoring( test_case=self)) self.add_step(RemapTidalMixing(test_case=self)) diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py b/compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py new file mode 100644 index 0000000000..654ff9848a --- /dev/null +++ b/compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py @@ -0,0 +1,125 @@ +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 NormalizeIcebergIceShelfMelt(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`` + """ + 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='normalize_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 input 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') + + def run(self): + """ + Run this step of the test case + """ + super().run() + + if not self.with_ice_shelf_cavities: + return + + 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 + + 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 var in ['bergFreshwaterFluxData']: + ds_dib[var] = ds_dib[var] / total_flux + + write_netcdf(ds_dib, 'dib_merino_2020_normalized.nc') + + for var in ['dataLandIceFreshwaterFlux', 'dataLandIceHeatFlux']: + ds_dismf[var] = ds_dismf[var] / total_flux + + write_netcdf(ds_dismf, 'dismf_paolo2023_normalized.nc') + + norm_total_dib_flux = (ds_dib.bergFreshwaterFluxData * weights * + area_cell).sum() + + norm_total_dismf_flux = (ds_dismf.dataLandIceFreshwaterFlux * + area_cell).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('') + + prefix = 'Iceberg_Climatology_Merino_normalized' + dest_filename = f'{prefix}.{suffix}.nc' + + symlink( + os.path.abspath('dib_merino_2020_normalized.nc'), + f'{self.ocean_inputdata_dir}/{dest_filename}') + + prefix = 'prescribed_ismf_paolo2023_normalized' + dest_filename = f'{prefix}.{suffix}.nc' + + symlink( + os.path.abspath('dismf_paolo2023_normalized.nc'), + 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 index c49cc941ab..fe5e1b09cd 100644 --- 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 @@ -40,7 +40,9 @@ def setup(self): setup input files based on config options """ super().setup() - if not self.with_ice_shelf_cavities: + if self.init is not None: + # we don't need any files, since we already did this remapping + # during init return filename = 'prescribed_ismf_paolo2023.nc' From 8b35df7f3969eeaab4b07e5fca0be906f83b0948 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 3 Jul 2024 04:46:03 -0500 Subject: [PATCH 3/3] 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. --- .../global_ocean/files_for_e3sm/__init__.py | 20 ++---- ...py => add_total_iceberg_ice_shelf_melt.py} | 67 +++++++++++-------- .../files_for_e3sm/remap_ice_shelf_melt.py | 51 +++++--------- .../remap_iceberg_climatology.py | 20 +++--- 4 files changed, 69 insertions(+), 89 deletions(-) rename compass/ocean/tests/global_ocean/files_for_e3sm/{normalize_iceberg_ice_shelf_melt.py => add_total_iceberg_ice_shelf_melt.py} (59%) diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py b/compass/ocean/tests/global_ocean/files_for_e3sm/__init__.py index fd220a1d16..6ffc8384ed 100644 --- 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/compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py b/compass/ocean/tests/global_ocean/files_for_e3sm/add_total_iceberg_ice_shelf_melt.py similarity index 59% rename from compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py rename to compass/ocean/tests/global_ocean/files_for_e3sm/add_total_iceberg_ice_shelf_melt.py index 654ff9848a..8d4a6da7f4 100644 --- a/compass/ocean/tests/global_ocean/files_for_e3sm/normalize_iceberg_ice_shelf_melt.py +++ b/compass/ocean/tests/global_ocean/files_for_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 index fe5e1b09cd..ee8c47a204 100644 --- 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,13 +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 self.init is not None: - # we don't need any files, since we already did this remapping - # during init - return filename = 'prescribed_ismf_paolo2023.nc' @@ -72,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) 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 index 24d4532b4d..952ab13569 100644 --- 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,