diff --git a/polaris/ocean/ice_shelf/__init__.py b/polaris/ocean/ice_shelf/__init__.py index fc88b4106a..bf8854f600 100644 --- a/polaris/ocean/ice_shelf/__init__.py +++ b/polaris/ocean/ice_shelf/__init__.py @@ -1,14 +1,13 @@ -from typing import Dict as Dict - -from polaris import ( - Step as Step, -) -from polaris import ( - Task as Task, +from polaris import Task +from polaris.ocean.ice_shelf.freeze import ( # noqa: F401 + compute_freezing_temperature, ) -from polaris.ocean.ice_shelf.ssh_adjustment import ( - SshAdjustment as SshAdjustment, +from polaris.ocean.ice_shelf.pressure import ( # noqa: F401 + compute_land_ice_draft_from_pressure, + compute_land_ice_pressure_from_draft, + compute_land_ice_pressure_from_thickness, ) +from polaris.ocean.ice_shelf.ssh_adjustment import SshAdjustment from polaris.ocean.ice_shelf.ssh_forward import SshForward as SshForward @@ -100,7 +99,7 @@ def setup_ssh_adjustment_steps( yaml_filename : str, optional the yaml filename used for ssh_forward steps - yaml_replacements : Dict, optional + yaml_replacements : dict, optional key, string combinations for templated replacements in the yaml file diff --git a/polaris/ocean/ice_shelf/freeze.cfg b/polaris/ocean/ice_shelf/freeze.cfg new file mode 100644 index 0000000000..f614827be7 --- /dev/null +++ b/polaris/ocean/ice_shelf/freeze.cfg @@ -0,0 +1,14 @@ +# Options related to freezing temperature below ice shelves +[ice_shelf_freeze] + +# The freezing temperature (C) at zero pressure and salinity +coeff_0 = 6.22e-2 + +# The coefficient (1e3 C) for the term proportional to salinity +coeff_S = -5.63e-2 + +# The coefficient (C Pa^-1) for the term proportional to the pressure +coeff_p = -7.43e-8 + +# The coefficient (1e3 C Pa^-1) for the term proportional to salinity times pressure +coeff_pS = -1.74e-10 diff --git a/polaris/ocean/ice_shelf/freeze.py b/polaris/ocean/ice_shelf/freeze.py new file mode 100644 index 0000000000..7e97b8bc17 --- /dev/null +++ b/polaris/ocean/ice_shelf/freeze.py @@ -0,0 +1,35 @@ +def compute_freezing_temperature(config, salinity, pressure): + """ + Get the freezing temperature in an ice-shelf cavity using the same equation + of state as in the ocean model + + Parameters + ---------- + config : polaris.config.PolarisConfigParser + Config options + + salinity : xarray.DataArray + The salinity field + + pressure + The pressure field + + Returns + ------- + freezing_temp : xarray.DataArray + The freezing temperature + """ + section = config['ice_shelf_freeze'] + coeff_0 = section.getfloat('coeff_0') + coeff_S = section.getfloat('coeff_S') + coeff_p = section.getfloat('coeff_p') + coeff_pS = section.getfloat('coeff_pS') + + freezing_temp = ( + coeff_0 + + coeff_S * salinity + + coeff_p * pressure + + coeff_pS * pressure * salinity + ) + + return freezing_temp diff --git a/polaris/ocean/ice_shelf/pressure.py b/polaris/ocean/ice_shelf/pressure.py new file mode 100644 index 0000000000..65305d2f0f --- /dev/null +++ b/polaris/ocean/ice_shelf/pressure.py @@ -0,0 +1,96 @@ +import numpy as np +from mpas_tools.cime.constants import constants + + +def compute_land_ice_pressure_from_thickness( + land_ice_thickness, modify_mask, land_ice_density=None +): + """ + Compute the pressure from an overlying ice shelf from ice thickness + + Parameters + ---------- + land_ice_thickness: xarray.DataArray + The ice thickness + + modify_mask : xarray.DataArray + A mask that is 1 where ``landIcePressure`` can be deviate from 0 + + land_ice_density : float, optional + A reference density for land ice + + Returns + ------- + land_ice_pressure : xarray.DataArray + The pressure from the overlying land ice on the ocean + """ + gravity = constants['SHR_CONST_G'] + if land_ice_density is None: + land_ice_density = constants['SHR_CONST_RHOICE'] + land_ice_pressure = modify_mask * np.maximum( + land_ice_density * gravity * land_ice_thickness, 0.0 + ) + return land_ice_pressure + + +def compute_land_ice_pressure_from_draft( + land_ice_draft, modify_mask, ref_density=None +): + """ + Compute the pressure from an overlying ice shelf from ice draft + + Parameters + ---------- + land_ice_draft : xarray.DataArray + The ice draft (sea surface height) + + modify_mask : xarray.DataArray + A mask that is 1 where ``landIcePressure`` can be deviate from 0 + + ref_density : float, optional + A reference density for seawater displaced by the ice shelf + + Returns + ------- + land_ice_pressure : xarray.DataArray + The pressure from the overlying land ice on the ocean + """ + gravity = constants['SHR_CONST_G'] + if ref_density is None: + ref_density = constants['SHR_CONST_RHOSW'] + land_ice_pressure = modify_mask * np.maximum( + -ref_density * gravity * land_ice_draft, 0.0 + ) + return land_ice_pressure + + +def compute_land_ice_draft_from_pressure( + land_ice_pressure, modify_mask, ref_density=None +): + """ + Compute the ice-shelf draft associated with the pressure from an overlying + ice shelf + + Parameters + ---------- + land_ice_pressure : xarray.DataArray + The pressure from the overlying land ice on the ocean + + modify_mask : xarray.DataArray + A mask that is 1 where ``landIcePressure`` can be deviate from 0 + + ref_density : float, optional + A reference density for seawater displaced by the ice shelf + + Returns + ------- + land_ice_draft : xarray.DataArray + The ice draft + """ + gravity = constants['SHR_CONST_G'] + if ref_density is None: + ref_density = constants['SHR_CONST_RHOSW'] + land_ice_draft = -( + modify_mask * land_ice_pressure / (ref_density * gravity) + ) + return land_ice_draft diff --git a/polaris/ocean/tasks/isomip_plus/__init__.py b/polaris/ocean/tasks/isomip_plus/__init__.py index eb3d8bbe3d..a015a67f88 100644 --- a/polaris/ocean/tasks/isomip_plus/__init__.py +++ b/polaris/ocean/tasks/isomip_plus/__init__.py @@ -72,17 +72,26 @@ def add_isomip_plus_tasks(component, mesh_type): mesh_type, resolution, mesh_name, resdir, component, config ) - for experiment in [ - 'ocean0', - 'ocean1', - 'ocean2', - 'ocean3', - 'ocean4', - 'inception', - 'wetting', - 'drying', - ]: - for vertical_coordinate in ['z-star']: + basic_expts = ['ocean0', 'ocean1', 'ocean2', 'ocean3', 'ocean4'] + extra_expts = ['wetting', 'drying'] + vert_coords = ['z-star', 'sigma'] + + for experiment in basic_expts: + vertical_coordinate = 'z-star' + task = IsomipPlusTest( + component=component, + resdir=resdir, + resolution=resolution, + experiment=experiment, + vertical_coordinate=vertical_coordinate, + planar=planar, + shared_steps=shared_steps[experiment], + ) + component.add_task(task) + + if planar and resolution >= 2.0: + for experiment in extra_expts: + vertical_coordinate = 'z-star' task = IsomipPlusTest( component=component, resdir=resdir, @@ -94,6 +103,20 @@ def add_isomip_plus_tasks(component, mesh_type): ) component.add_task(task) + for vertical_coordinate in vert_coords: + for experiment in ['ocean0'] + extra_expts: + task = IsomipPlusTest( + component=component, + resdir=resdir, + resolution=resolution, + experiment=experiment, + vertical_coordinate=vertical_coordinate, + planar=planar, + shared_steps=shared_steps[experiment], + thin_film=True, + ) + component.add_task(task) + def _get_shared_steps( mesh_type, resolution, mesh_name, resdir, component, config @@ -193,7 +216,7 @@ def _get_shared_steps( # ocean0 and ocean1 use the same topography shared_steps['ocean0'] = shared_steps['ocean1'] - for experiment in ['inception', 'wetting', 'drying']: + for experiment in ['wetting', 'drying']: shared_steps[experiment] = dict(shared_steps['ocean1']) subdir = f'{resdir}/topo/scale/{experiment}' topo_scale = TopoScale( diff --git a/polaris/ocean/tasks/isomip_plus/init/__init__.py b/polaris/ocean/tasks/isomip_plus/init/__init__.py new file mode 100644 index 0000000000..ccfd12247e --- /dev/null +++ b/polaris/ocean/tasks/isomip_plus/init/__init__.py @@ -0,0 +1,2 @@ +from polaris.ocean.tasks.isomip_plus.init.forcing import Forcing as Forcing +from polaris.ocean.tasks.isomip_plus.init.init import Init as Init diff --git a/polaris/ocean/tasks/isomip_plus/init/forcing.py b/polaris/ocean/tasks/isomip_plus/init/forcing.py new file mode 100644 index 0000000000..8483d5b475 --- /dev/null +++ b/polaris/ocean/tasks/isomip_plus/init/forcing.py @@ -0,0 +1,222 @@ +import os + +import numpy as np +import xarray as xr +from mpas_tools.cime.constants import constants +from mpas_tools.io import write_netcdf + +from polaris.step import Step + + +class Forcing(Step): + """ + A step for creating forcing files for ISOMIP+ experiments + + Attributes + ---------- + resolution : float + The horizontal resolution (km) of the test case + + experiment : {'ocean0', 'ocean1', 'ocean2', 'ocean3', 'ocean4'} + The ISOMIP+ experiment + + vertical_coordinate : str + The type of vertical coordinate (``z-star``, ``z-level``, etc.) + + thin_film: bool + Whether the run includes a thin film below grounded ice + """ + + def __init__( + self, + component, + indir, + culled_mesh, + topo, + resolution, + experiment, + vertical_coordinate, + thin_film, + ): + """ + Create the step + + Parameters + ---------- + component : polaris.Component + The component this step belongs to + + indir : str, optional + the directory the step is in, to which ``name`` will be appended + + culled_mesh : polaris.Step + The step that culled the MPAS mesh + + topo : polaris.Step + The step with topography data on the culled mesh + + resolution : float + The horizontal resolution (km) of the test case + + experiment : {'ocean0', 'ocean1', 'ocean2', 'ocean3', 'ocean4'} + The ISOMIP+ experiment + + vertical_coordinate : str + The type of vertical coordinate (``z-star``, ``z-level``, etc.) + + thin_film: bool + Whether the run includes a thin film below grounded ice + """ + super().__init__(component=component, name='forcing', indir=indir) + self.resolution = resolution + self.experiment = experiment + self.vertical_coordinate = vertical_coordinate + self.thin_film = thin_film + + self.add_input_file( + filename='mesh.nc', + work_dir_target=os.path.join(culled_mesh.path, 'culled_mesh.nc'), + ) + + self.add_input_file( + filename='topo.nc', + work_dir_target=os.path.join(topo.path, 'topography_remapped.nc'), + ) + + self.add_input_file(filename='init.nc', target='../init/init.nc') + + self.add_output_file('init.nc') + + def run(self): + """ + Create restoring and other forcing files + """ + self._compute_restoring() + self._write_time_varying_forcing() + + def _compute_restoring(self): + config = self.config + experiment = self.experiment + + ref_density = constants['SHR_CONST_RHOSW'] + + ds_mesh = xr.open_dataset('mesh.nc') + ds_init = xr.open_dataset('init.nc') + + ds_forcing = xr.Dataset() + + section = config['isomip_plus'] + if experiment in ['ocean0', 'ocean1', 'ocean3']: + top_temp = section.getfloat('warm_top_temp') + bot_temp = section.getfloat('warm_bot_temp') + top_sal = section.getfloat('warm_top_sal') + bot_sal = section.getfloat('warm_bot_sal') + else: + top_temp = section.getfloat('cold_top_temp') + bot_temp = section.getfloat('cold_bot_temp') + top_sal = section.getfloat('cold_top_sal') + bot_sal = section.getfloat('cold_bot_sal') + + section = config['isomip_plus_forcing'] + restore_rate = section.getfloat('restore_rate') + restore_xmin = section.getfloat('restore_xmin') + restore_xmax = section.getfloat('restore_xmax') + restore_evap_rate = section.getfloat('restore_evap_rate') + + max_bottom_depth = -config.getfloat('vertical_grid', 'bottom_depth') + z_frac = (0.0 - ds_init.zMid) / (0.0 - max_bottom_depth) + + ds_forcing['temperatureInteriorRestoringValue'] = ( + 1.0 - z_frac + ) * top_temp + z_frac * bot_temp + ds_forcing['salinityInteriorRestoringValue'] = ( + 1.0 - z_frac + ) * top_sal + z_frac * bot_sal + + x_frac = np.maximum( + ( + (ds_mesh.xIsomipCell - restore_xmin) + / (restore_xmax - restore_xmin) + ), + 0.0, + ) + x_frac = x_frac.broadcast_like( + ds_forcing.temperatureInteriorRestoringValue + ) + + # convert from 1/days to 1/s + ds_forcing['temperatureInteriorRestoringRate'] = ( + x_frac * restore_rate / constants['SHR_CONST_CDAY'] + ) + ds_forcing['salinityInteriorRestoringRate'] = ( + ds_forcing.temperatureInteriorRestoringRate + ) + + # compute "evaporation" + mask = np.logical_and( + ds_mesh.xIsomipCell >= restore_xmin, + ds_mesh.xIsomipCell <= restore_xmax, + ) + mask = mask.expand_dims(dim='Time', axis=0) + # convert to m/s, negative for evaporation rather than precipitation + evap_rate = -restore_evap_rate / (constants['SHR_CONST_CDAY'] * 365) + # PSU*m/s to kg/m^2/s + sflux_factor = 1.0 + # C*m/s to W/m^2 + hflux_factor = 1.0 / (ref_density * constants['SHR_CONST_CPSW']) + ds_forcing['evaporationFlux'] = mask * ref_density * evap_rate + ds_forcing['seaIceSalinityFlux'] = ( + mask * evap_rate * top_sal / sflux_factor + ) + ds_forcing['seaIceHeatFlux'] = ( + mask * evap_rate * top_temp / hflux_factor + ) + + if self.vertical_coordinate == 'single-layer': + x_max = np.max(ds_mesh.xIsomipCell.values) + ds_forcing['tidalInputMask'] = xr.where( + ds_mesh.xIsomipCell > (x_max - 0.6 * self.resolution * 1e3), + 1.0, + 0.0, + ) + else: + ds_forcing['tidalInputMask'] = xr.zeros_like(ds_mesh.xIsomipCell) + + write_netcdf(ds_forcing, 'init_mode_forcing_data.nc') + + @staticmethod + def _write_time_varying_forcing(): + """ + Write time-varying land-ice forcing and update the initial condition + """ + + ds_topo = xr.open_dataset('topo.nc') + + if 'Time' not in ds_topo.dims: + # no time-varying forcing needed + return + + ds_out = xr.Dataset() + ds_out['xtime'] = ds_topo.xtime + + ds_out['landIceDraftForcing'] = ds_topo.landIceDraft + ds_out.landIceDraftForcing.attrs['units'] = 'm' + ds_out.landIceDraftForcing.attrs['long_name'] = ( + 'The approximate elevation of the land ice-ocean interface' + ) + ds_out['landIcePressureForcing'] = ds_topo.landIcePressure + ds_out.landIcePressureForcing.attrs['units'] = 'm' + ds_out.landIcePressureForcing.attrs['long_name'] = ( + 'Pressure from the weight of land ice at the ice-ocean interface' + ) + ds_out['landIceFractionForcing'] = ds_topo.landIceFraction + ds_out.landIceFractionForcing.attrs['long_name'] = ( + 'The fraction of each cell covered by land ice' + ) + ds_out['landIceFloatingFractionForcing'] = ( + ds_topo.landIceFloatingFraction + ) + ds_out.landIceFloatingFractionForcing.attrs['long_name'] = ( + 'The fraction of each cell covered by floating land ice' + ) + write_netcdf(ds_out, 'land_ice_forcing.nc') diff --git a/polaris/ocean/tasks/isomip_plus/init/init.py b/polaris/ocean/tasks/isomip_plus/init/init.py new file mode 100644 index 0000000000..9fddb7ebcd --- /dev/null +++ b/polaris/ocean/tasks/isomip_plus/init/init.py @@ -0,0 +1,579 @@ +import os +import shutil + +import cmocean # noqa: F401 +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf + +from polaris.ocean.ice_shelf import ( + compute_freezing_temperature, + compute_land_ice_draft_from_pressure, +) +from polaris.ocean.vertical import init_vertical_coord +from polaris.ocean.viz import compute_transect, plot_transect +from polaris.step import Step +from polaris.viz import plot_horiz_field + + +class Init(Step): + """ + A step for creating an initial condition for ISOMIP+ experiments + + Attributes + ---------- + experiment : {'ocean0', 'ocean1', 'ocean2', 'ocean3', 'ocean4'} + The ISOMIP+ experiment + + vertical_coordinate : str + The type of vertical coordinate (``z-star``, ``z-level``, etc.) + + thin_film: bool + Whether the run includes a thin film below grounded ice + """ + + def __init__( + self, + component, + indir, + culled_mesh, + topo, + experiment, + vertical_coordinate, + thin_film, + ): + """ + Create the step + + Parameters + ---------- + component : polaris.Component + The component this step belongs to + + indir : str, optional + the directory the step is in, to which ``name`` will be appended + + culled_mesh : polaris.Step + The step that culled the MPAS mesh + + topo : polaris.Step + The step with topography data on the culled mesh + + experiment : {'ocean0', 'ocean1', 'ocean2', 'ocean3', 'ocean4'} + The ISOMIP+ experiment + + vertical_coordinate : str + The type of vertical coordinate (``z-star``, ``z-level``, etc.) + + thin_film: bool + Whether the run includes a thin film below grounded ice + """ + super().__init__(component=component, name='init', indir=indir) + self.experiment = experiment + self.vertical_coordinate = vertical_coordinate + self.thin_film = thin_film + + self.add_input_file( + filename='mesh.nc', + work_dir_target=os.path.join(culled_mesh.path, 'culled_mesh.nc'), + ) + + self.add_input_file( + filename='topo.nc', + work_dir_target=os.path.join(topo.path, 'topography_remapped.nc'), + ) + + self.add_output_file('init.nc') + + def run(self): + """ + Create and plot the initial condition + """ + config = self.config + self._compute_init_topo() + + with xr.open_dataset('init_topo.nc') as ds_init: + init_vertical_coord(config, ds_init) + write_netcdf(ds_init, 'init_vert_coord.nc') + + self._compute_t_s_vel_coriolis() + self._plot() + + def _compute_init_topo(self): + """Compute fractions, masks, draft, land-ice pressure, bot. depth""" + config = self.config + logger = self.logger + thin_film = self.thin_film + + section = config['isomip_plus_init'] + min_levels = section.getint('minimum_levels') + ocean_density = config.getfloat('isomip_plus', 'ocean_density') + min_layer_thickness = section.getfloat('min_layer_thickness') + if thin_film: + min_column_thickness = section.getfloat( + 'min_column_thickness_thin_film' + ) + else: + min_column_thickness = section.getfloat( + 'min_column_thickness_no_thin_film' + ) + min_land_ice_fraction = section.getfloat('min_land_ice_fraction') + + ds_topo = xr.open_dataset('topo.nc') + if 'Time' in ds_topo.dims: + ds_topo = ds_topo.isel(Time=0) + ds_mesh = xr.open_dataset('mesh.nc') + + ds_init = xr.Dataset(ds_mesh) + ds_init.attrs = ds_mesh.attrs + + ds_init['landIceFraction'] = ds_topo['landIceFraction'] + ds_init['landIceFloatingFraction'] = ds_topo['landIceFloatingFraction'] + + # This inequality needs to be > rather than >= to ensure correctness + # when min_land_ice_fraction = 0 + land_ice_mask = ds_init.landIceFraction > min_land_ice_fraction + + floating_mask = np.logical_and( + land_ice_mask, ds_init.landIceFloatingFraction > 0 + ) + + ds_init['landIceMask'] = land_ice_mask.astype(int) + ds_init['landIceFloatingMask'] = floating_mask.astype(int) + + ds_init['landIceFraction'] = xr.where( + land_ice_mask, ds_init.landIceFraction, 0.0 + ) + ds_init['landIceFloatingFraction'] = xr.where( + land_ice_mask, ds_init.landIceFloatingFraction, 0.0 + ) + + ds_init['landIceGroundedFraction'] = ds_topo['landIceGroundedFraction'] + + # start from landIcePressure and compute landIceDraft + ds_init['landIcePressure'] = ds_topo['landIcePressure'] + + land_ice_draft = compute_land_ice_draft_from_pressure( + land_ice_pressure=ds_topo.landIcePressure, + modify_mask=ds_topo.landIcePressure > 0.0, + ref_density=ocean_density, + ) + + land_ice_draft = np.maximum(ds_topo.bedrockTopography, land_ice_draft) + + ds_init['landIceDraft'] = land_ice_draft + + ds_init['ssh'] = ds_init.landIceDraft + + ds_init['bottomDepth'] = -ds_topo.bedrockTopography + + thin_film_mask = np.logical_and( + land_ice_mask, np.logical_not(floating_mask) + ) + + if thin_film: + active_ocean = xr.ones_like(ds_init.bottomDepth) + else: + active_ocean = np.logical_not(thin_film_mask) + + logger.info( + f'Not using a thin film so {np.sum(thin_film_mask).values} ' + f'cells are being deactivated.' + ) + + # need to move bottomDepth and ssh up to the sea surface where we + # are grounded so we will get maxLevelCell == 0 + ds_init['bottomDepth'] = xr.where( + thin_film_mask, 0.0, ds_init.bottomDepth + ) + ds_init['ssh'] = xr.where(thin_film_mask, 0.0, ds_init.ssh) + + # Deepen the bottom depth to maintain the minimum water-column + # thickness + min_column_thickness = max( + min_column_thickness, min_levels * min_layer_thickness + ) + min_depth = -ds_init.ssh + min_column_thickness + too_thin = np.logical_and( + active_ocean, ds_init.bottomDepth < min_depth + ) + ds_init['bottomDepth'] = xr.where( + too_thin, min_depth, ds_init.bottomDepth + ) + too_thin_count = np.sum(too_thin).values + if too_thin_count > 0: + logger.info( + f'Adjusted bottomDepth for {too_thin_count} cells to achieve ' + f'minimum column thickness of {min_column_thickness}' + ) + + ds_init['activeOcean'] = active_ocean.astype(int) + ds_init['thinFilmMask'] = thin_film_mask.astype(int) + + write_netcdf(ds_init, 'init_topo.nc') + + def _compute_t_s_vel_coriolis(self): + config = self.config + experiment = self.experiment + thin_film = self.thin_film + + section = config['isomip_plus_init'] + + if experiment in ['ocean0', 'ocean3']: + top_temp = config.getfloat('isomip_plus', 'warm_top_temp') + bot_temp = config.getfloat('isomip_plus', 'warm_bot_temp') + top_sal = config.getfloat('isomip_plus', 'warm_top_sal') + bot_sal = config.getfloat('isomip_plus', 'warm_bot_sal') + else: + top_temp = config.getfloat('isomip_plus', 'cold_top_temp') + bot_temp = config.getfloat('isomip_plus', 'cold_bot_temp') + top_sal = config.getfloat('isomip_plus', 'cold_top_sal') + bot_sal = config.getfloat('isomip_plus', 'cold_bot_sal') + + ds_mesh = xr.open_dataset('mesh.nc') + ds_init = xr.open_dataset('init_vert_coord.nc') + + active_ocean = ds_init.activeOcean == 1 + thin_film_mask = ds_init.thinFilmMask == 1 + + if not thin_film: + # need to set maxLevelCell == 0 where ocean isn't active + ds_init['maxLevelCell'] = ds_init.maxLevelCell.where( + active_ocean, 0 + ) + + max_bottom_depth = -config.getfloat('vertical_grid', 'bottom_depth') + frac = (0.0 - ds_init.zMid) / (0.0 - max_bottom_depth) + + # compute T, S + ds_init['temperature'] = (1.0 - frac) * top_temp + frac * bot_temp + ds_init['salinity'] = (1.0 - frac) * top_sal + frac * bot_sal + + # Note that using the land ice pressure rather than the pressure at + # floatation will mean that there is a small amount of cooling from + # grounding line retreat. However, the thin film should be thin enough + # that this effect isn't significant. + + freezing_temp = compute_freezing_temperature( + config=config, + salinity=ds_init.salinity, + pressure=ds_init.landIcePressure, + ) + ds_init['temperature'] = ds_init.temperature.where( + np.logical_not(thin_film_mask), freezing_temp + ) + + # compute coriolis + coriolis_parameter = section.getfloat('coriolis_parameter') + + ds_init['fCell'] = coriolis_parameter * xr.ones_like(ds_mesh.xCell) + ds_init['fEdge'] = coriolis_parameter * xr.ones_like(ds_mesh.xEdge) + ds_init['fVertex'] = coriolis_parameter * xr.ones_like(ds_mesh.xVertex) + + normalVelocity = xr.zeros_like(ds_mesh.xEdge) + normalVelocity = normalVelocity.broadcast_like(ds_init.refBottomDepth) + normalVelocity = normalVelocity.transpose('nEdges', 'nVertLevels') + ds_init['normalVelocity'] = normalVelocity.expand_dims( + dim='Time', axis=0 + ) + + write_netcdf(ds_init, 'init.nc') + + def _plot(self): + """ + Plot several fields from the initial condition + """ + section = self.config['isomip_plus_init'] + thin_film = self.thin_film + if thin_film: + min_column_thickness = section.getfloat( + 'min_column_thickness_thin_film' + ) + else: + min_column_thickness = section.getfloat( + 'min_column_thickness_no_thin_film' + ) + min_layer_thickness = section.getfloat('min_layer_thickness') + + tol = 1e-10 + + plot_folder = 'plots' + try: + shutil.rmtree(plot_folder) + except FileNotFoundError: + pass + + ds = xr.open_dataset('init.nc').isel(Time=0) + ds_mesh = xr.open_dataset('mesh.nc') + + for ds_fix in [ds, ds_mesh]: + # use the planar projection coordinates for horizontal plots + ds_fix['xCell'] = ds_fix['xIsomipCell'] + ds_fix['yCell'] = ds_fix['yIsomipCell'] + ds_fix['zCell'] = xr.zeros_like(ds_fix.xCell) + ds_fix['xVertex'] = ds_fix['xIsomipVertex'] + ds_fix['yVertex'] = ds_fix['yIsomipVertex'] + ds_fix['zVertex'] = xr.zeros_like(ds_fix.xVertex) + + x = np.linspace(320e3, 800e3, 2) + y = 40.0e3 * np.ones_like(x) + + x = xr.DataArray(data=x, dims=('nPoints',)) + y = xr.DataArray(data=y, dims=('nPoints',)) + + ds_transect = compute_transect( + x=x, + y=y, + ds_horiz_mesh=ds_mesh, + layer_thickness=ds.layerThickness, + bottom_depth=ds.bottomDepth, + min_level_cell=ds.minLevelCell - 1, + max_level_cell=ds.maxLevelCell - 1, + spherical=False, + ) + + ds['totalColThickness'] = ds['layerThickness'].sum(dim='nVertLevels') + + ds['H'] = ds.ssh + ds.bottomDepth + + figsize = (8, 4.5) + + cell_mask = ds.maxLevelCell >= 1 + descriptor = plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['maxLevelCell'], + out_file_name='plots/maxLevelCell.png', + figsize=figsize, + field_mask=cell_mask, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['landIceMask'], + out_file_name='plots/landIceMask.png', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['landIceFloatingMask'], + out_file_name='plots/landIceFloatingMask.png', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['landIcePressure'], + out_file_name='plots/landIcePressure.png', + vmin=1e4, + vmax=1e7, + cmap_scale='log', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['ssh'], + out_file_name='plots/ssh.png', + vmin=-720, + vmax=0, + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['bottomDepth'], + out_file_name='plots/bottomDepth.png', + vmin=0, + vmax=720, + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['totalColThickness'], + out_file_name='plots/totalColThickness.png', + vmin=min_column_thickness + tol, + vmax=720, + cmap_set_under='r', + cmap_scale='log', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['H'], + out_file_name='plots/H.png', + vmin=min_column_thickness + tol, + vmax=720, + cmap_set_under='r', + cmap_scale='log', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['landIceFraction'], + out_file_name='plots/landIceFraction.png', + vmin=0 + tol, + vmax=1 - tol, + cmap='cmo.balance', + cmap_set_under='k', + cmap_set_over='r', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['landIceFloatingFraction'], + out_file_name='plots/landIceFloatingFraction.png', + vmin=0 + tol, + vmax=1 - tol, + cmap='cmo.balance', + cmap_set_under='k', + cmap_set_over='r', + figsize=figsize, + descriptor=descriptor, + ) + + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds['landIceGroundedFraction'], + out_file_name='plots/landIceGroundedFraction.png', + vmin=0 + tol, + vmax=1 - tol, + cmap='cmo.balance', + cmap_set_under='k', + cmap_set_over='r', + figsize=figsize, + descriptor=descriptor, + ) + + plot_transect( + ds_transect=ds_transect, + title='layer interfaces at y=40 km', + out_filename='plots/layerInterfacesSection.png', + figsize=figsize, + interface_color='black', + ssh_color='blue', + seafloor_color='red', + ) + + _plot_top_bot_slice( + ds=ds, + ds_mesh=ds_mesh, + ds_transect=ds_transect, + descriptor=descriptor, + field_name='layerThickness', + figsize=figsize, + vmin=min_layer_thickness + tol, + vmax=50, + cmap='cmo.deep_r', + units='m', + transect_x=x, + transect_y=y, + under='r', + over='r', + ) + + _plot_top_bot_slice( + ds=ds, + ds_mesh=ds_mesh, + ds_transect=ds_transect, + descriptor=descriptor, + field_name='zMid', + figsize=figsize, + vmin=-720.0, + vmax=0.0, + cmap='cmo.deep_r', + units='m', + transect_x=x, + transect_y=y, + under='r', + over='r', + ) + + _plot_top_bot_slice( + ds=ds, + ds_mesh=ds_mesh, + ds_transect=ds_transect, + descriptor=descriptor, + field_name='temperature', + figsize=figsize, + vmin=-2.0, + vmax=1.0, + cmap='cmo.thermal', + units=r'$^\circ$C', + transect_x=x, + transect_y=y, + ) + + _plot_top_bot_slice( + ds=ds, + ds_mesh=ds_mesh, + ds_transect=ds_transect, + descriptor=descriptor, + field_name='salinity', + figsize=figsize, + vmin=33.8, + vmax=34.7, + cmap='cmo.haline', + units='PSU', + transect_x=x, + transect_y=y, + ) + + +def _plot_top_bot_slice( + ds, + ds_mesh, + ds_transect, + descriptor, + field_name, + figsize, + vmin, + vmax, + cmap, + units, + transect_x, + transect_y, + under=None, + over=None, +): + for suffix, z_index in [['Top', 0], ['Bot', ds.maxLevelCell - 1]]: + plot_horiz_field( + ds_mesh=ds_mesh, + field=ds[field_name], + title=f'{field_name} {suffix}', + z_index=z_index, + out_file_name=f'plots/{field_name}{suffix}.png', + vmin=vmin, + vmax=vmax, + cmap=cmap, + cmap_set_under=under, + cmap_set_over=over, + figsize=figsize, + descriptor=descriptor, + transect_x=transect_x, + transect_y=transect_y, + ) + + plot_transect( + ds_transect=ds_transect, + mpas_field=ds[field_name], + title=f'{field_name} at y=40 km', + out_filename=f'plots/{field_name}Section.png', + vmin=vmin, + vmax=vmax, + cmap=cmap, + figsize=figsize, + colorbar_label=units, + ) diff --git a/polaris/ocean/tasks/isomip_plus/isomip_plus.cfg b/polaris/ocean/tasks/isomip_plus/isomip_plus.cfg index aa38f46ca0..71d8168e05 100644 --- a/polaris/ocean/tasks/isomip_plus/isomip_plus.cfg +++ b/polaris/ocean/tasks/isomip_plus/isomip_plus.cfg @@ -3,3 +3,26 @@ # the density of ice prescribed in ISOMIP+ ice_density = 918 + +# the density of ocean water prescribed in ISOMIP+ +ocean_density = 1028 + +# the WARM profile +# the initial temperature at the sea surface +warm_top_temp = -1.9 +# the initial temperature at the sea floor +warm_bot_temp = 1.0 +# the initial salinity at the sea surface +warm_top_sal = 33.8 +# the initial salinity at the sea floor +warm_bot_sal = 34.7 + +# the COLD profile +# the initial temperature at the sea surface +cold_top_temp = -1.9 +# the initial temperature at the sea floor +cold_bot_temp = -1.9 +# the initial salinity at the sea surface +cold_top_sal = 33.8 +# the initial salinity at the sea floor +cold_bot_sal = 34.55 diff --git a/polaris/ocean/tasks/isomip_plus/isomip_plus_init.cfg b/polaris/ocean/tasks/isomip_plus/isomip_plus_init.cfg new file mode 100644 index 0000000000..418fc3defe --- /dev/null +++ b/polaris/ocean/tasks/isomip_plus/isomip_plus_init.cfg @@ -0,0 +1,57 @@ +# Options related to the vertical grid +[vertical_grid] + +# the type of vertical grid +grid_type = uniform + +# Number of vertical levels (altered for sigma tests in code) +vert_levels = 36 + +# Depth of the bottom of the ocean +bottom_depth = 720.0 + +# The type of vertical coordinate (e.g. z-level, z-star) +coord_type = z-star + +# Whether to use "partial" or "full", or "None" to not alter the topography +partial_cell_type = None + +# The minimum fraction of a layer for partial cells +min_pc_fraction = 0.1 + + +# config options for ISOMIP+ init +[isomip_plus_init] + +# Minimum number of vertical levels in a column +minimum_levels = 3 + +# The minimum allowable layer thickness +min_layer_thickness = 0.0 + +# Minimum thickness of the initial ocean column (to prevent 'drying') without +# and with a thin film region +min_column_thickness_no_thin_film = 1.1e-2 +min_column_thickness_thin_film = 1e-3 + +# Minimum fraction of a cell that contains land ice in order for it to be +# considered a land-ice cell by MPAS-Ocean (landIceMask == 1). +min_land_ice_fraction = 0.0 + +# Coriolis parameter (1/s) for entire domain +coriolis_parameter = -1.409e-4 + +# config options for ISOMIP+ forcing +[isomip_plus_forcing] + +# restoring rate (1/days) at the open-ocean boundary +restore_rate = 10.0 + +# the "evaporation" rate (m/yr) near the open-ocean boundary used to keep sea +# level from rising +restore_evap_rate = 200.0 + +# southern boundary of restoring region (m) +restore_xmin = 790e3 +# northern boundary of restoring region (m) +restore_xmax = 800e3 diff --git a/polaris/ocean/tasks/isomip_plus/isomip_plus_test.py b/polaris/ocean/tasks/isomip_plus/isomip_plus_test.py index d32f9198de..1923e839a4 100644 --- a/polaris/ocean/tasks/isomip_plus/isomip_plus_test.py +++ b/polaris/ocean/tasks/isomip_plus/isomip_plus_test.py @@ -1,4 +1,5 @@ from polaris import Task +from polaris.ocean.tasks.isomip_plus.init import Forcing, Init class IsomipPlusTest(Task): @@ -92,11 +93,51 @@ def __init__( continue self.add_step(step, symlink=symlink) + self.add_step( + Init( + component=component, + indir=subdir, + culled_mesh=shared_steps['topo/cull_mesh'], + topo=shared_steps['topo_final'], + experiment=experiment, + vertical_coordinate=vertical_coordinate, + thin_film=thin_film, + ) + ) + + self.add_step( + Forcing( + component=component, + indir=subdir, + culled_mesh=shared_steps['topo/cull_mesh'], + topo=shared_steps['topo_final'], + resolution=resolution, + experiment=experiment, + vertical_coordinate=vertical_coordinate, + thin_film=thin_film, + ) + ) + def configure(self): """ Modify the configuration options for this test case. """ config = self.config + config.add_from_package('polaris.ocean.ice_shelf', 'freeze.cfg') config.add_from_package( 'polaris.ocean.tasks.isomip_plus', 'isomip_plus.cfg' ) + config.add_from_package( + 'polaris.ocean.tasks.isomip_plus', 'isomip_plus_init.cfg' + ) + vertical_coordinate = self.vertical_coordinate + + # for most coordinates, use the config options + levels = None + + if vertical_coordinate == 'sigma': + levels = 10 + + config.set('vertical_grid', 'coord_type', vertical_coordinate) + if levels is not None: + config.set('vertical_grid', 'vert_levels', f'{levels}') diff --git a/polaris/ocean/tasks/isomip_plus/isomip_plus_topo.cfg b/polaris/ocean/tasks/isomip_plus/isomip_plus_topo.cfg index 8e83ccc49c..9d81f350fa 100644 --- a/polaris/ocean/tasks/isomip_plus/isomip_plus_topo.cfg +++ b/polaris/ocean/tasks/isomip_plus/isomip_plus_topo.cfg @@ -27,17 +27,12 @@ expand_factor = 2.0 # simple thickening and thinning experiments that involve scaling the Ocean1 # landIcePressure and landIceDraft over time # -# "inception" reference dates -inception_dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00 -# scaling at each date -inception_scales = 0.0, 1.0, 1.0 - # "drying" reference dates drying_dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00 # scaling at each date -drying_scales = 1.0, 2.0, 2.0 +drying_scales = 1.0, 1.05, 1.10 # "wetting" reference dates wetting_dates = 0001-01-01_00:00:00, 0002-01-01_00:00:00, 0003-01-01_00:00:00 # scaling at each date -wetting_scales = 1.0, 0.0, 0.0 +wetting_scales = 1.0, 0.95, 0.9 diff --git a/polaris/ocean/tasks/isomip_plus/topo/scale.py b/polaris/ocean/tasks/isomip_plus/topo/scale.py index 19cfb399e1..59345c7473 100644 --- a/polaris/ocean/tasks/isomip_plus/topo/scale.py +++ b/polaris/ocean/tasks/isomip_plus/topo/scale.py @@ -14,7 +14,7 @@ class TopoScale(Step): Attributes ---------- - experiment : {'inception', 'drying', 'wetting'} + experiment : {'drying', 'wetting'} The scaling experiment """ @@ -36,10 +36,10 @@ def __init__(self, component, subdir, config, topo_remap, experiment): topo_remap : polaris.ocean.tasks.isomip_plus.topo_map.TopoRemap The step that produced the remapped topography to be scaled - experiment : {'inception', 'drying', 'wetting'} + experiment : {'drying', 'wetting'} The scaling experiment """ - if experiment not in ['inception', 'drying', 'wetting']: + if experiment not in ['drying', 'wetting']: raise ValueError(f'Unexpected experiment {experiment}') super().__init__(component=component, name='topo_scale', subdir=subdir)