From 0a856203ed411fa84bfbbccb1fb434ca57aeeb08 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 9 Apr 2024 08:12:39 -0500 Subject: [PATCH] Reorganize SSH adjustment This merge reorganizes SSH adjustment substantially. It is now a sequence of alternating init and forward runs, with decreasing minimum water column thickness. The sea surface height is adjusted instead of `landIcePressure` for better consistency with planned coupling to MALI. --- compass/ocean/tests/global_ocean/forward.py | 8 +- .../ocean/tests/global_ocean/global_ocean.cfg | 8 +- .../ocean/tests/global_ocean/init/__init__.py | 90 ++++++++++++--- .../tests/global_ocean/init/initial_state.py | 46 ++++++-- .../tests/global_ocean/init/ssh_adjustment.py | 109 ++++++++++++++++-- .../init/ssh_from_surface_density.py | 102 ++++++++++++++++ .../tests/global_ocean/init/streams.init | 1 + .../global_ocean/init/streams.ssh_adjust | 26 +++++ .../tests/global_ocean/init/streams.topo | 1 + 9 files changed, 346 insertions(+), 45 deletions(-) create mode 100644 compass/ocean/tests/global_ocean/init/ssh_from_surface_density.py diff --git a/compass/ocean/tests/global_ocean/forward.py b/compass/ocean/tests/global_ocean/forward.py index 43cc8d9ada..22c4ca3ca3 100644 --- a/compass/ocean/tests/global_ocean/forward.py +++ b/compass/ocean/tests/global_ocean/forward.py @@ -139,12 +139,8 @@ def __init__(self, test_case, mesh, time_integrator, init=None, work_dir_target=f'{mesh_path}/culled_mesh.nc') if init is not None: - if mesh.with_ice_shelf_cavities: - initial_state_target = \ - f'{init.path}/ssh_adjustment/adjusted_init.nc' - else: - initial_state_target = \ - f'{init.path}/initial_state/initial_state.nc' + initial_state_target = \ + f'{init.path}/initial_state/initial_state.nc' self.add_input_file(filename='init.nc', work_dir_target=initial_state_target) self.add_input_file( diff --git a/compass/ocean/tests/global_ocean/global_ocean.cfg b/compass/ocean/tests/global_ocean/global_ocean.cfg index 3c7a0e6de5..3a3bcda6a2 100644 --- a/compass/ocean/tests/global_ocean/global_ocean.cfg +++ b/compass/ocean/tests/global_ocean/global_ocean.cfg @@ -21,7 +21,7 @@ sweep_count = 20 [ssh_adjustment] # the number of iterations of ssh adjustment to perform -iterations = 10 +iterations = 4 # Whether to convert adjusted initial condition files to CDF5 format during # ssh adjustment under ice shelves @@ -54,8 +54,10 @@ btr_dt_per_km = 1.5 min_levels = 3 cavity_min_levels = ${min_levels} -# minimum thickness of layers in ice-shelf cavities -cavity_min_layer_thickness = 1.0 +# minimum thickness of layers in ice-shelf cavities at the beginning and end +# of iterative ssh init +cavity_min_layer_thickness_initial = 10.0 +cavity_min_layer_thickness_final = 3.0 # Maximum allowed Haney number for configurations with ice-shelf cavities rx1_max = 20.0 diff --git a/compass/ocean/tests/global_ocean/init/__init__.py b/compass/ocean/tests/global_ocean/init/__init__.py index c59841abeb..5d092dc506 100644 --- a/compass/ocean/tests/global_ocean/init/__init__.py +++ b/compass/ocean/tests/global_ocean/init/__init__.py @@ -5,6 +5,9 @@ RemapIceShelfMelt, ) from compass.ocean.tests.global_ocean.init.ssh_adjustment import SshAdjustment +from compass.ocean.tests.global_ocean.init.ssh_from_surface_density import ( + SshFromSurfaceDensity, +) from compass.testcase import TestCase from compass.validate import compare_variables @@ -50,17 +53,6 @@ def __init__(self, test_group, mesh, initial_condition): self.mesh = mesh self.initial_condition = initial_condition - self.add_step( - InitialState( - test_case=self, mesh=mesh, - initial_condition=initial_condition)) - - if mesh.with_ice_shelf_cavities: - self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh)) - - self.add_step( - SshAdjustment(test_case=self)) - def configure(self, config=None): """ Modify the configuration options for this test case @@ -68,11 +60,14 @@ def configure(self, config=None): config : compass.config.CompassConfigParser, optional Configuration options to update if not those for this test case """ + add_steps = config is None if config is None: config = self.config + mesh = self.mesh + # set mesh-relate config options - self.mesh.configure(config=config) + mesh.configure(config=config) initial_condition = self.initial_condition descriptions = {'WOA23': 'World Ocean Atlas 2023 climatology ' @@ -84,6 +79,72 @@ def configure(self, config=None): config.set('global_ocean', 'init_description', descriptions[initial_condition]) + if add_steps: + # add the steps for ssh adjustment + if mesh.with_ice_shelf_cavities: + step_index = 1 + name = \ + f'{step_index:02d}_init_with_draft_from_constant_density' + subdir = f'adjust_ssh/{name}' + init_const_rho = InitialState( + test_case=self, mesh=mesh, + initial_condition=initial_condition, + name=name, subdir=subdir, + adjustment_fraction=0.) + self.add_step(init_const_rho) + + # TODO: Recompute landIceDraft using surface density + step_index += 1 + name = f'{step_index:02d}_ssh_from_surface_density' + subdir = f'adjust_ssh/{name}' + ssh_from_surf_rho = SshFromSurfaceDensity( + test_case=self, init_path=init_const_rho.path, + name=name, subdir=subdir) + self.add_step(ssh_from_surf_rho) + + culled_topo_path = ssh_from_surf_rho.path + + iteration_count = config.getint('ssh_adjustment', 'iterations') + for iter_index in range(iteration_count): + fraction = iter_index / iteration_count + + step_index += 1 + name = f'{step_index:02d}_init' + subdir = f'adjust_ssh/{name}' + init_step = InitialState( + test_case=self, mesh=mesh, + initial_condition=initial_condition, + culled_topo_path=culled_topo_path, + name=name, subdir=subdir, + adjustment_fraction=fraction) + self.add_step(init_step) + + step_index += 1 + name = f'{step_index:02d}_adjust_ssh' + subdir = f'adjust_ssh/{name}' + adjust_ssh = SshAdjustment( + test_case=self, init_path=init_step.path, + name=name, subdir=subdir) + self.add_step(adjust_ssh) + culled_topo_path = adjust_ssh.path + + name = 'initial_state' + subdir = 'initial_state' + init_step = InitialState( + test_case=self, mesh=mesh, + initial_condition=initial_condition, + culled_topo_path=culled_topo_path, + name=name, subdir=subdir, + adjustment_fraction=1.0) + self.add_step(init_step) + + self.add_step(RemapIceShelfMelt(test_case=self, mesh=mesh)) + else: + self.add_step( + InitialState( + test_case=self, mesh=mesh, + initial_condition=initial_condition)) + def validate(self): """ Test cases can override this method to perform validation of variables @@ -92,8 +153,3 @@ def validate(self): variables = ['temperature', 'salinity', 'layerThickness'] compare_variables(test_case=self, variables=variables, filename1='initial_state/initial_state.nc') - - if self.mesh.with_ice_shelf_cavities: - variables = ['ssh', 'landIcePressure'] - compare_variables(test_case=self, variables=variables, - filename1='ssh_adjustment/adjusted_init.nc') diff --git a/compass/ocean/tests/global_ocean/init/initial_state.py b/compass/ocean/tests/global_ocean/init/initial_state.py index 17511846c1..71f47e94c6 100644 --- a/compass/ocean/tests/global_ocean/init/initial_state.py +++ b/compass/ocean/tests/global_ocean/init/initial_state.py @@ -28,8 +28,14 @@ class InitialState(Step): initial_condition : {'WOA23', 'PHC', 'EN4_1900'} The initial condition dataset to use + + adjustment_fraction : float + The fraction of the way through iterative ssh adjustment for this + step """ - def __init__(self, test_case, mesh, initial_condition): + def __init__(self, test_case, mesh, initial_condition, + culled_topo_path=None, name='initial_state', subdir=None, + adjustment_fraction=None): """ Create the step @@ -43,13 +49,30 @@ def __init__(self, test_case, mesh, initial_condition): initial_condition : {'WOA23', 'PHC', 'EN4_1900'} The initial condition dataset to use + + culled_topo_path : str, optional + The path to a step where ``topography_culled.nc`` is provided + + name : str, optional + The name of the step + + subdir : str, optional + The subdirectory for the step + + adjustment_fraction : float, optional + The fraction of the way through iterative ssh adjustment for this + step """ if initial_condition not in ['WOA23', 'PHC', 'EN4_1900']: raise ValueError(f'Unknown initial_condition {initial_condition}') - super().__init__(test_case=test_case, name='initial_state') + super().__init__(test_case=test_case, name=name, subdir=subdir) self.mesh = mesh self.initial_condition = initial_condition + if mesh.with_ice_shelf_cavities and adjustment_fraction is None: + raise ValueError('Must provide adjustment_fraction for ' + 'initializing meshes with ice-shelf cavities') + self.adjustment_fraction = adjustment_fraction package = 'compass.ocean.tests.global_ocean.init' @@ -83,8 +106,10 @@ def __init__(self, test_case, mesh, initial_condition): self.add_namelist_options(options, mode='init') self.add_streams_file(package, 'streams.topo', mode='init') - cull_step = self.mesh.steps['cull_mesh'] - target = os.path.join(cull_step.path, 'topography_culled.nc') + if culled_topo_path is None: + cull_step = self.mesh.steps['cull_mesh'] + culled_topo_path = cull_step.path + target = os.path.join(culled_topo_path, 'topography_culled.nc') self.add_input_file(filename='topography_culled.nc', work_dir_target=target) @@ -168,6 +193,7 @@ def run(self): Run this step of the testcase """ config = self.config + section = config['global_ocean'] self._smooth_topography() interfaces = generate_1d_grid(config=config) @@ -185,10 +211,14 @@ def run(self): namelist = {'config_global_ocean_minimum_depth': f'{min_depth}'} if self.mesh.with_ice_shelf_cavities: - cavity_min_levels = \ - config.getint('global_ocean', 'cavity_min_levels') - cavity_min_layer_thickness = \ - config.getfloat('global_ocean', 'cavity_min_layer_thickness') + frac = self.adjustment_fraction + cavity_min_levels = section.getint('cavity_min_levels') + min_thick_init = section.getfloat( + 'cavity_min_layer_thickness_initial') + min_thick_final = section.getfloat( + 'cavity_min_layer_thickness_final') + cavity_min_layer_thickness = ( + (1.0 - frac) * min_thick_init + frac * min_thick_final) namelist['config_rx1_min_levels'] = f'{cavity_min_levels}' namelist['config_rx1_min_layer_thickness'] = \ f'{cavity_min_layer_thickness}' diff --git a/compass/ocean/tests/global_ocean/init/ssh_adjustment.py b/compass/ocean/tests/global_ocean/init/ssh_adjustment.py index 2d988a887e..d6ce70a948 100644 --- a/compass/ocean/tests/global_ocean/init/ssh_adjustment.py +++ b/compass/ocean/tests/global_ocean/init/ssh_adjustment.py @@ -1,6 +1,11 @@ from importlib.resources import contents -from compass.ocean.iceshelf import adjust_ssh +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call + +from compass.model import partition, run_model from compass.ocean.tests.global_ocean.forward import ForwardStep @@ -9,7 +14,8 @@ class SshAdjustment(ForwardStep): A step for iteratively adjusting the pressure from the weight of the ice shelf to match the sea-surface height as part of ice-shelf 2D test cases """ - def __init__(self, test_case): + def __init__(self, test_case, init_path, name='ssh_adjustment', + subdir=None): """ Create the step @@ -17,16 +23,24 @@ def __init__(self, test_case): ---------- test_case : compass.ocean.tests.global_ocean.init.Init The test case this step belongs to + + init_path : str + The path to the initial state to use for forward runs + + name : str, optional + The name of the step + + subdir : str, optional + The subdirectory for the step """ super().__init__(test_case=test_case, mesh=test_case.mesh, time_integrator='split_explicit_ab2', - name='ssh_adjustment') + name=name, subdir=subdir) self.add_namelist_options({'config_AM_globalStats_enable': '.false.'}) self.add_namelist_file('compass.ocean.namelists', 'namelist.ssh_adjust') - self.add_streams_file('compass.ocean.streams', 'streams.ssh_adjust') self.add_streams_file('compass.ocean.tests.global_ocean.init', 'streams.ssh_adjust') @@ -41,10 +55,9 @@ def __init__(self, test_case): self.add_streams_file(mesh_package, mesh_streams) mesh_path = test_case.mesh.get_cull_mesh_path() - init_path = test_case.steps['initial_state'].path self.add_input_file( - filename='adjusting_init0.nc', + filename='init.nc', work_dir_target=f'{init_path}/initial_state.nc') self.add_input_file( filename='forcing_data.nc', @@ -53,17 +66,91 @@ def __init__(self, test_case): filename='graph.info', work_dir_target=f'{mesh_path}/culled_graph.info') - self.add_output_file(filename='adjusted_init.nc') + self.add_input_file( + filename='original_topograpy.nc', + work_dir_target=f'{mesh_path}/topography_culled.nc') + + self.add_output_file(filename='topography_culled.nc') def run(self): """ Run this step of the testcase """ config = self.config - iteration_count = config.getint('ssh_adjustment', 'iterations') update_pio = config.getboolean('global_ocean', 'forward_update_pio') convert_to_cdf5 = config.getboolean('ssh_adjustment', 'convert_to_cdf5') - adjust_ssh(variable='landIcePressure', iteration_count=iteration_count, - step=self, update_pio=update_pio, - convert_to_cdf5=convert_to_cdf5) + + self._adjust_ssh(update_pio=update_pio, + convert_to_cdf5=convert_to_cdf5) + + def _adjust_ssh(self, update_pio, convert_to_cdf5): + """ + Adjust the sea surface height to be dynamically consistent with + land-ice pressure. + """ + ntasks = self.ntasks + config = self.config + logger = self.logger + + if update_pio: + self.update_namelist_pio('namelist.ocean') + partition(ntasks, config, logger) + + with xr.open_dataset('init.nc') as ds: + ds = ds.isel(Time=0) + init_ssh = ds.ssh + modify_mask = np.logical_and(ds.maxLevelCell > 0, + ds.modifyLandIcePressureMask == 1) + land_ice_pressure = ds.landIcePressure + + logger.info(" * Running forward model") + run_model(self, update_pio=False, partition_graph=False) + logger.info(" - Complete") + + logger.info(" * Updating SSH") + + with xr.open_dataset('output_ssh.nc') as ds_ssh: + # get the last time entry + ds_ssh = ds_ssh.isel(Time=ds_ssh.sizes['Time'] - 1) + final_ssh = ds_ssh.ssh + + delta_ssh = modify_mask * (final_ssh - init_ssh) + + with xr.open_dataset('original_topograpy.nc') as ds_out: + + ds_out['ssh'] = modify_mask * final_ssh + + if convert_to_cdf5: + write_filename = 'topography_before_cdf5.nc' + else: + write_filename = 'topography_culled.nc' + write_netcdf(ds_out, write_filename) + if convert_to_cdf5: + args = ['ncks', '-O', '-5', 'topography_before_cdf5.nc', + 'topography_culled.nc'] + check_call(args, logger=logger) + + # Write the largest change in SSH and its lon/lat to a file + with open('maxDeltaSSH.log', 'w') as log_file: + + icell = np.abs(delta_ssh).argmax().values + + ds_cell = ds.isel(nCells=icell) + delta_ssh_max = delta_ssh.isel(nCells=icell).values + + coords = (f'lon/lat: ' + f'{np.rad2deg(ds_cell.lonCell.values):f} ' + f'{np.rad2deg(ds_cell.latCell.values):f}') + string = (f'delta_ssh_max: ' + f'{delta_ssh_max:g}, {coords}') + + logger.info(f' {string}') + log_file.write(f'{string}\n') + string = (f'ssh: {final_ssh.isel(nCells=icell).values:g}, ' + f'landIcePressure: ' + f'{land_ice_pressure.isel(nCells=icell).values:g}') + logger.info(f' {string}') + log_file.write(f'{string}\n') + + logger.info(" - Complete\n") diff --git a/compass/ocean/tests/global_ocean/init/ssh_from_surface_density.py b/compass/ocean/tests/global_ocean/init/ssh_from_surface_density.py new file mode 100644 index 0000000000..08f38c68b4 --- /dev/null +++ b/compass/ocean/tests/global_ocean/init/ssh_from_surface_density.py @@ -0,0 +1,102 @@ + +import numpy as np +import xarray as xr +from mpas_tools.cime.constants import constants +from mpas_tools.io import write_netcdf +from mpas_tools.logging import check_call + +from compass.step import Step + + +class SshFromSurfaceDensity(Step): + """ + Compute the sea surface height that is in approximate hydrostatic balance + with a given land-ice pressure using the density in the top layer of the + ocean + + Attributes + ---------- + mesh : compass.ocean.tests.global_ocean.mesh.mesh.MeshStep + The step for creating the mesh + """ + def __init__(self, test_case, init_path, name, subdir): + """ + Create the step + + Parameters + ---------- + test_case : compass.ocean.tests.global_ocean.init.Init + The test case this step belongs to + + init_path : str + The path to the initial state from which to get the land ice + pressure and surface density + + name : str + The name of the step + + subdir : str + The subdirectory for the step + """ + super().__init__(test_case=test_case, name=name, subdir=subdir) + self.mesh = test_case.mesh + self.add_input_file( + filename='init.nc', + work_dir_target=f'{init_path}/initial_state.nc') + + mesh_path = self.mesh.get_cull_mesh_path() + + self.add_input_file( + filename='mesh.nc', + work_dir_target=f'{mesh_path}/culled_mesh.nc') + + self.add_input_file( + filename='original_topograpy.nc', + work_dir_target=f'{mesh_path}/topography_culled.nc') + + self.add_output_file(filename='topography_culled.nc') + + def run(self): + """ + Run this step of the testcase + """ + config = self.config + logger = self.logger + + convert_to_cdf5 = config.getboolean('ssh_adjustment', + 'convert_to_cdf5') + + g = constants['SHR_CONST_G'] + + with xr.open_dataset('init.nc') as ds_init: + ds_init = ds_init.isel(Time=0) + modify_mask = np.logical_and( + ds_init.maxLevelCell > 0, + ds_init.modifyLandIcePressureMask == 1) + land_ice_pressure = ds_init.landIcePressure + + if 'minLevelCell' in ds_init: + min_level_cell = ds_init.minLevelCell - 1 + else: + min_level_cell = xr.zeros_like(ds_init.maxLevelCell) + + top_density = ds_init.density.isel(nVertLevels=min_level_cell) + + ssh = xr.where(modify_mask, + - land_ice_pressure / (top_density * g), + 0.) + + with xr.open_dataset('original_topograpy.nc') as ds_topo: + + ds_topo['ssh'] = ssh + + if convert_to_cdf5: + write_filename = 'topography_before_cdf5.nc' + else: + write_filename = 'topography_culled.nc' + write_netcdf(ds_topo, write_filename) + + if convert_to_cdf5: + args = ['ncks', '-O', '-5', 'topography_before_cdf5.nc', + 'topography_culled.nc'] + check_call(args, logger=logger) diff --git a/compass/ocean/tests/global_ocean/init/streams.init b/compass/ocean/tests/global_ocean/init/streams.init index 1dfe341489..8d2232c9c8 100644 --- a/compass/ocean/tests/global_ocean/init/streams.init +++ b/compass/ocean/tests/global_ocean/init/streams.init @@ -45,6 +45,7 @@ + diff --git a/compass/ocean/tests/global_ocean/init/streams.ssh_adjust b/compass/ocean/tests/global_ocean/init/streams.ssh_adjust index 6199d4cbb8..f752dd1c93 100644 --- a/compass/ocean/tests/global_ocean/init/streams.ssh_adjust +++ b/compass/ocean/tests/global_ocean/init/streams.ssh_adjust @@ -1,5 +1,31 @@ + + + + + + + + + + + + + + diff --git a/compass/ocean/tests/global_ocean/init/streams.topo b/compass/ocean/tests/global_ocean/init/streams.topo index e0b0de9563..13fd2e36a1 100644 --- a/compass/ocean/tests/global_ocean/init/streams.topo +++ b/compass/ocean/tests/global_ocean/init/streams.topo @@ -11,6 +11,7 @@ +