Add steps for ISOMIP+ init and forcing

polaris/ocean/ice_shelf/__init__.py

@@ -1,6 +1,12 @@
 from typing import Dict
 
-from polaris import Step, Task
+from polaris import Task
+from polaris.ocean.ice_shelf.freeze import compute_freezing_temperature
+from polaris.ocean.ice_shelf.pressure import (
+    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
 

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

polaris/ocean/ice_shelf/freeze.py

@@ -0,0 +1,33 @@
+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

polaris/ocean/ice_shelf/pressure.py

@@ -0,0 +1,90 @@
+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.)
+    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.)
+    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

polaris/ocean/tasks/isomip_plus/__init__.py

@@ -45,9 +45,25 @@ def add_isomip_plus_tasks(component, mesh_type):
         shared_steps = _get_shared_steps(
             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 = ['inception', 'wetting', 'drying']
+        vert_coords = ['z-star', 'sigma', 'single-layer']
+
+        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.:
+            for experiment in extra_expts:
+                vertical_coordinate = 'z-star'
                 task = IsomipPlusTest(
                     component=component,
                     resdir=resdir,
@@ -58,6 +74,19 @@ def add_isomip_plus_tasks(component, mesh_type):
                     shared_steps=shared_steps[experiment])
                 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):

polaris/ocean/tasks/isomip_plus/init/__init__.py

@@ -0,0 +1,2 @@
+from polaris.ocean.tasks.isomip_plus.init.forcing import Forcing
+from polaris.ocean.tasks.isomip_plus.init.init import Init

polaris/ocean/tasks/isomip_plus/init/forcing.py

@@ -0,0 +1,191 @@
+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. - ds_init.zMid) / (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.)
+        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.
+        # C*m/s to W/m^2
+        hflux_factor = 1. / (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')