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 afa6ead95..6ffc8384e 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, ) @@ -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: diff --git a/compass/ocean/tests/global_ocean/files_for_e3sm/add_total_iceberg_ice_shelf_melt.py b/compass/ocean/tests/global_ocean/files_for_e3sm/add_total_iceberg_ice_shelf_melt.py new file mode 100644 index 000000000..8d4a6da7f --- /dev/null +++ b/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}') 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 c49cc941a..ee8c47a20 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,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) 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 24d4532b4..952ab1356 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, 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 4980b81ac..fc010135d 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