Back to home page

MITgcm

 
 

    


Warning, /doc/phys_pkgs/shelfice.rst is written in an unsupported language. File is not indexed.

view on githubraw file Latest commit 00f81e67 on 2024-05-24 21:00:12 UTC
b938a3c63b antn*0001 .. _sub_phys_pkg_shelfice:
                0002 
274298f12b Jeff*0003 SHELFICE Package
                0004 ----------------
                0005 
                0006 Authors: Martin Losch, Jean-Michel Campin
                0007 
0a8794f5ee Mart*0008 Introduction
274298f12b Jeff*0009 ~~~~~~~~~~~~
                0010 
                0011 :filelink:`pkg/shelfice` provides a thermodynamic model for basal melting
                0012 underneath floating ice shelves.
                0013 
                0014 CPP options enable or disable different aspects of the package
0a8794f5ee Mart*0015 (:numref:`shelfice_config`). Run-time options, flags, filenames and
274298f12b Jeff*0016 field-related dates/times are described in :numref:`shelfice_runtime`. A description of key subroutines is given
0a8794f5ee Mart*0017 in :numref:`shelfice_subroutines`. Available diagnostics output is listed in
274298f12b Jeff*0018 :numref:`shelfice_diagnostics`.
                0019 
                0020 
                0021 .. _shelfice_config:
                0022 
                0023 SHELFICE configuration
                0024 ~~~~~~~~~~~~~~~~~~~~~~
                0025  
                0026 As with all MITgcm packages, :filelink:`pkg/shelfice` can be turned on or off at compile
                0027 time:
                0028 
                0029 -  using the ``packages.conf`` file by adding ``shelfice`` to it,
                0030 
                0031 -  or using :filelink:`genmake2 <tools/genmake2>` adding ``-enable=shelfice`` or ``disable=shelfice`` switches
                0032 
                0033 :filelink:`pkg/shelfice` does not require any additional packages, but it will only
                0034 work with conventional vertical :math:`z`-coordinates (pressure
                0035 coordinates are not implemented). If you use it together with
                0036 vertical mixing schemes, be aware that non-local parameterizations
                0037 are turned off, e.g., such as :filelink:`pkg/kpp`.
                0038 
                0039 Parts of the :filelink:`pkg/shelfice` code can be enabled or disabled at compile time
                0040 via CPP preprocessor flags. These options are set in :filelink:`SHELFICE_OPTIONS.h <pkg/shelfice/SHELFICE_OPTIONS.h>`:
                0041 
7b8b86ab99 Timo*0042 .. tabularcolumns:: |\Y{.32}|\Y{.1}|\Y{.570}|
                0043 
                0044 .. table:: Compile-time parameters
                0045    :name: tab_phys_pkg_shelfice_compileparms
                0046 
                0047    +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
                0048    | CPP Flag Name                                 | Default | Description                                                                                                          |
                0049    +===============================================+=========+======================================================================================================================+
                0050    | :varlink:`ALLOW_SHELFICE_DEBUG`               | #undef  | include code for enhanced diagnostics and debug output                                                               |
                0051    +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
                0052    | :varlink:`ALLOW_ISOMIP_TD`                    | #define | include code for for simplified ISOMIP thermodynamics                                                                |
                0053    +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
                0054    | :varlink:`SHI_ALLOW_GAMMAFRICT`               | #define | allow friction velocity-dependent transfer coefficient following Holland and Jenkins (1999) :cite:`holland:99`       |
                0055    +-----------------------------------------------+---------+----------------------------------------------------------------------------------------------------------------------+
274298f12b Jeff*0056 
                0057 .. _shelfice_runtime:
                0058 
0a8794f5ee Mart*0059 SHELFICE run-time parameters
274298f12b Jeff*0060 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                0061 
7b8b86ab99 Timo*0062 :filelink:`pkg/shelfice` is switched on/off at run time by setting :varlink:`useSHELFICE` to ``.TRUE.`` in file ``data.pkg``.
274298f12b Jeff*0063 Run-time parameters are set in file ``data.shelfice`` (read in :filelink:`pkg/shelfice/shelfice_readparms.F`),as listed below.
                0064 
                0065 The data file specifying under-ice topography of ice shelves (:varlink:`SHELFICEtopoFile`) is in meters; upwards is positive,
                0066 and as for the bathymetry files, negative values are required for topography below the sea-level.
                0067 The data file for the pressure load anomaly at the bottom of the ice shelves :varlink:`SHELFICEloadAnomalyFile` is in pressure
                0068 units (Pa). This field is absolutely required to avoid large
0a8794f5ee Mart*0069 excursions of the free surface during initial adjustment processes,
274298f12b Jeff*0070 obtained by integrating an approximate density from the surface at
                0071 :math:`z=0` down to the bottom of the last fully dry cell within the
                0072 ice shelf, see :eq:`surfacepressure`. Note however the file :varlink:`SHELFICEloadAnomalyFile` must
                0073 not be :math:`p_{top}`, but
67dae37801 Eli *0074 :math:`p_{top}-g\sum_{k'=1}^{n-1}\rho_c \Delta{z}_{k'}`, with
                0075 :math:`\rho_c =` :varlink:`rhoConst`, so that in the absence of a :math:`\rho^{*}`
                0076 that is different from :math:`\rho_c`, the anomaly is zero.
274298f12b Jeff*0077 
7b8b86ab99 Timo*0078 .. tabularcolumns:: |\Y{.275}|\Y{.28}|\Y{.455}|
274298f12b Jeff*0079 
7b8b86ab99 Timo*0080 .. table:: Run-time parameters and default values; all parameters are in namelist group ``SHELFICE_PARM01``
0a8794f5ee Mart*0081    :name: tab_phys_pkg_shelfice_runtimeparms
                0082    :class: longtable
                0083 
7b8b86ab99 Timo*0084    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0085    | Parameter                              | Default                                    | Description                                                                                             |
                0086    +========================================+============================================+=========================================================================================================+
                0087    | :varlink:`useISOMIPTD`                 | FALSE                                      | use simplified ISOMIP thermodynamics on/off flag                                                        |
                0088    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0089    | :varlink:`SHELFICEconserve`            | FALSE                                      | use conservative form of temperature boundary conditions on/off flag                                    |
                0090    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0091    | :varlink:`SHELFICEboundaryLayer`       | FALSE                                      | use simple boundary layer mixing parameterization on/off flag                                           |
                0092    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0093    | :varlink:`SHI_withBL_realFWflux`       | FALSE                                      | with :varlink:`SHELFICEboundaryLayer`, allow to use real-FW flux                                        |
                0094    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0095    | :varlink:`SHI_withBL_uStarTopDz`       | FALSE                                      | with :varlink:`SHELFICEboundaryLayer`, compute uStar from uVel,vVel averaged over top Dz thickness      |
                0096    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0097    | :varlink:`SHELFICEloadAnomalyFile`     | :kbd:`' '`                                 | initial geopotential anomaly                                                                            |
                0098    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0099    | :varlink:`SHELFICEtopoFile`            | :kbd:`' '`                                 | filename for under-ice topography of ice shelves                                                        |
                0100    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0101    | :varlink:`SHELFICEmassFile`            | :kbd:`' '`                                 | filename for mass of ice shelves                                                                        |
                0102    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0103    | :varlink:`SHELFICEMassDynTendFile`     | :kbd:`' '`                                 | filename for mass tendency of ice shelves                                                               |
                0104    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0105    | :varlink:`SHELFICETransCoeffTFile`     | :kbd:`' '`                                 | filename for spatially varying transfer coefficients                                                    |
                0106    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0107    | :varlink:`SHELFICElatentHeat`          | 334.0E+03                                  | latent heat of fusion (J/kg)                                                                            |
                0108    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0109    | :varlink:`SHELFICEHeatCapacity_Cp`     | 2000.0E+00                                 | specific heat capacity of ice (J/kg/K)                                                                  |
                0110    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0111    | :varlink:`rhoShelfIce`                 | 917.0E+00                                  | (constant) mean density of ice shelf (kg/m\ :sup:`3`)                                                   |
                0112    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
00f81e6785 Ou W*0113    | :varlink:`SHELFICEsalinity`            | 0.0E+00                                    | (constant) salinity of ice shelf                                                                        |
                0114    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
7b8b86ab99 Timo*0115    | :varlink:`SHELFICEheatTransCoeff`      | 1.0E-04                                    | transfer coefficient (exchange velocity) for temperature (m/s)                                          |
                0116    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0117    | :varlink:`SHELFICEsaltTransCoeff`      | :varlink:`SHELFICEsaltToHeatRatio` *       | transfer coefficient (exchange velocity) for salinity (m/s)                                             |
                0118    |                                        | :varlink:`SHELFICEheatTransCoeff`          |                                                                                                         |
                0119    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0120    | :varlink:`SHELFICEsaltToHeatRatio`     | 5.05E-03                                   | ratio of salinity to temperature transfer coefficients (non-dim.)                                       |
                0121    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0122    | :varlink:`SHELFICEkappa`               | 1.54E-06                                   | temperature diffusion coefficient of the ice shelf (m\ :sup:`2`\ /s)                                    |
                0123    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0124    | :varlink:`SHELFICEthetaSurface`        | -20.0E+00                                  | (constant) surface temperature above the ice shelf (:sup:`o`\ C)                                        |
                0125    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0126    | :varlink:`no_slip_shelfice`            | :varlink:`no_slip_bottom`                  | slip along bottom of ice shelf on/off flag                                                              |
                0127    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0128    | :varlink:`SHELFICEDragLinear`          | :varlink:`bottomDragLinear`                | linear drag coefficient at bottom ice shelf (m/s)                                                       |
                0129    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0130    | :varlink:`SHELFICEDragQuadratic`       | :varlink:`bottomDragQuadratic`             | quadratic drag coefficient at bottom ice shelf (non-dim.)                                               |
                0131    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0132    | :varlink:`SHELFICEselectDragQuadr`     | -1                                         | select form of quadratic drag coefficient (non-dim.)                                                    |
                0133    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0134    | :varlink:`SHELFICEMassStepping`        | FALSE                                      | recalculate ice shelf mass at every time step                                                           |
                0135    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0136    | :varlink:`SHELFICEDynMassOnly`         | FALSE                                      | if :varlink:`SHELFICEmassStepping` = TRUE, exclude freshwater flux contribution                         |
                0137    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0138    | :varlink:`SHELFICEadvDiffHeatFlux`     | FALSE                                      | use advective-diffusive heat flux into ice shelf instead of default diffusive heat flux                 |
                0139    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0140    | :varlink:`SHELFICEuseGammaFrict`       | FALSE                                      | use velocity dependent exchange coefficients (Holland and Jenkins 1999 :cite:`holland:99`)              |
                0141    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0142    | :varlink:`SHELFICE_oldCalcUStar`       | FALSE                                      | use old uStar averaging expression                                                                      |
                0143    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0144    | :varlink:`SHELFICEwriteState`          | FALSE                                      | write ice shelf state to file on/off flag                                                               |
                0145    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0146    | :varlink:`SHELFICE_dumpFreq`           | :varlink:`dumpFreq`                        | dump frequency (s)                                                                                      |
                0147    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
                0148    | :varlink:`SHELFICE_dump_mnc`           | :varlink:`snapshot_mnc`                    | write snapshot using MNC  on/off flag                                                                   |
                0149    +----------------------------------------+--------------------------------------------+---------------------------------------------------------------------------------------------------------+
0a8794f5ee Mart*0150 
                0151 SHELFICE description
274298f12b Jeff*0152 ~~~~~~~~~~~~~~~~~~~~
                0153 
                0154 In the light of isomorphic equations for pressure and height
                0155 coordinates, the ice shelf topography on top of the water column has a
bd7056e9aa Jeff*0156 similar role as (and in the language of Marshall et al. (2004) :cite:`marshall:04`,
274298f12b Jeff*0157 is isomorphic to) the orography and the pressure boundary conditions at
                0158 the bottom of the fluid for atmospheric and oceanic models in pressure
0bad585a21 Navi*0159 coordinates. The total pressure :math:`p_{\rm tot}` in the ocean can be
274298f12b Jeff*0160 divided into the pressure at the top of the water column
0bad585a21 Navi*0161 :math:`p_{\rm top}`, the hydrostatic pressure and the non-hydrostatic
                0162 pressure contribution :math:`p_{\rm nh}`:
274298f12b Jeff*0163 
                0164 .. math::
0bad585a21 Navi*0165    p_{\rm tot} = p_{\rm top} + \int_z^{\eta-h} g\,\rho\,dz + p_{\rm nh}
274298f12b Jeff*0166    :label: pressureocean
                0167 
                0168 
                0169 with the gravitational acceleration :math:`g`, the density
                0170 :math:`\rho`, the vertical coordinate :math:`z` (positive upwards), and
                0171 the dynamic sea-surface height :math:`\eta`. For the open ocean,
0bad585a21 Navi*0172 :math:`p_{\rm top}=p_{a}` (atmospheric pressure) and :math:`h=0`. Underneath
274298f12b Jeff*0173 an ice-shelf that is assumed to be floating in isostatic equilibrium,
0bad585a21 Navi*0174 :math:`p_{\rm top}` at the top of the water column is the atmospheric
274298f12b Jeff*0175 pressure :math:`p_{a}` plus the weight of the ice-shelf. It is this
                0176 weight of the ice-shelf that has to be provided as a boundary condition
                0177 at the top of the water column (in run-time parameter :varlink:`SHELFICEloadAnomalyFile`). The weight is
                0178 conveniently computed by integrating a density profile :math:`\rho^*`,
                0179 that is constant in time and corresponds to the sea-water replaced by
                0180 ice, from :math:`z=0` to a “reference” ice-shelf draft at :math:`z=-h` (Beckmann et al. (1999)
                0181 :cite:`beckmann:99`), so that
                0182 
                0183 .. math::
0bad585a21 Navi*0184    p_{\rm top} = p_{a} + \int_{-h}^{0}g\,\rho^{*}\,dz
274298f12b Jeff*0185    :label: ptop
                0186 
                0187 Underneath the ice shelf, the “sea-surface height” :math:`\eta` is the
                0188 deviation from the “reference” ice-shelf draft :math:`h`. During a model
                0189 integration, :math:`\eta` adjusts so that the isostatic equilibrium is
                0190 maintained for sufficiently slow and large scale motion.
                0191 
0bad585a21 Navi*0192 In MITgcm, the total pressure anomaly :math:`p'_{\rm tot}` which is used
67dae37801 Eli *0193 for pressure gradient computations is defined by subtracting a purely
                0194 depth dependent contribution :math:`-g\rho_c z` using constant
0bad585a21 Navi*0195 reference density :math:`\rho_c` from :math:`p_{\rm tot}`.
274298f12b Jeff*0196 :eq:`pressureocean` becomes
                0197 
                0198 .. math::
0bad585a21 Navi*0199      p_{\rm tot} = p_{\rm top} - g \rho_c (z+h)  + g \rho_c \eta + \, \int_z^{\eta-h}{ g (\rho-\rho_c) \, dz} + \, p_{\rm nh}
274298f12b Jeff*0200      :label: pressure
                0201 
                0202 and after rearranging
                0203 
                0204 .. math::
0bad585a21 Navi*0205    p'_{\rm tot} = p'_{\rm top} + g \rho_c \eta + \, \int_z^{\eta-h}{g (\rho-\rho_c) \, dz} + \, p_{\rm nh}
274298f12b Jeff*0206 
0bad585a21 Navi*0207 with :math:`p'_{\rm tot} = p_{\rm tot} + g\,\rho_c\,z` and
                0208 :math:`p'_{\rm top} = p_{\rm top} - g\,\rho_c\,h`.
                0209 The non-hydrostatic pressure contribution :math:`p_{\rm nh}`
274298f12b Jeff*0210 is neglected in the following.
                0211 
0bad585a21 Navi*0212 In practice, the ice shelf contribution to :math:`p_{\rm top}` is computed
274298f12b Jeff*0213 by integrating :eq:`ptop` from :math:`z=0` to the bottom of the
                0214 last fully dry cell within the ice shelf:
                0215 
                0216 .. math::
0bad585a21 Navi*0217    p_{\rm top} = g\,\sum_{k'=1}^{n-1}\rho_{k'}^{*}\Delta{z_{k'}} + p_{a}
274298f12b Jeff*0218    :label: surfacepressure
                0219 
                0220 where :math:`n` is the vertical index of the first (at least partially)
                0221 “wet” cell and :math:`\Delta{z_{k'}}` is the thickness of the
                0222 :math:`k'`-th layer (counting downwards). The pressure anomaly for
                0223 evaluating the pressure gradient is computed in the center of the “wet”
                0224 cell :math:`k` as
                0225 
                0226 .. math::
0bad585a21 Navi*0227    p'_{k} = p'_{\rm top} + g\rho_{n}\eta +
67dae37801 Eli *0228    g\,\sum_{k'=n}^{k}\left((\rho_{k'}-\rho_c)\Delta{z_{k'}}
274298f12b Jeff*0229      \frac{1+H(k'-k)}{2}\right)
                0230    :label: discretizedpressure
                0231 
                0232 where :math:`H(k'-k)=1` for :math:`k'<k` and :math:`0` otherwise.
                0233 
bd7056e9aa Jeff*0234  .. figure:: figs/gridschematic.*
                0235     :width: 80%
                0236     :align: center
                0237     :alt: schematic of vertical section of grid
                0238     :name: shelfice_grid
                0239 
0a8794f5ee Mart*0240     Schematic of a vertical section of the grid at the base of an ice shelf. Grid lines are thin;
                0241     the thick line is the model’s representation of the ice shelf-water interface. Plus signs mark the position
                0242     of pressure points for pressure gradient computations. The letters A, B, and C mark specific grid cells for
                0243     reference. :math:`h_k` is the fractional cell thickness so that :math:`h_k \Delta z_k` is the actual cell thickness.
bd7056e9aa Jeff*0244 
                0245 
274298f12b Jeff*0246 Setting :varlink:`SHELFICEboundaryLayer` ``=.TRUE.`` introduces a simple boundary layer that reduces the potential
                0247 noise problem at the cost of increased vertical mixing. For this purpose
                0248 the water temperature at the :math:`k`-th layer abutting ice shelf
                0249 topography for use in the heat flux parameterizations is computed as a
                0250 mean temperature :math:`\overline{\theta}_{k}` over a boundary layer of
                0251 the same thickness as the layer thickness :math:`\Delta{z}_{k}`:
                0252 
                0253 .. math::
                0254    \overline{\theta}_{k} = \theta_{k} h_{k} + \theta_{k+1} (1-h_{k})
                0255    :label: thetabl
                0256 
                0257 where :math:`h_{k}\in[0,1]` is the fractional layer thickness of the
bd7056e9aa Jeff*0258 :math:`k`-th layer (see :numref:`shelfice_grid`). The original contributions due to ice shelf-ocean
274298f12b Jeff*0259 interaction :math:`g_{\theta}` to the total tendency terms
                0260 :math:`G_{\theta}` in the time-stepping equation
                0261 :math:`\theta^{n+1} = f(\theta^{n},\Delta{t},G_{\theta}^{n})` are
                0262 
                0263 .. math::
67dae37801 Eli *0264    g_{\theta,k}   = \frac{Q}{\rho_c c_{p} h_{k} \Delta{z}_{k}}
274298f12b Jeff*0265    \text{ and } g_{\theta,k+1} = 0
                0266    :label: orgtendency
                0267 
                0268 for layers :math:`k` and :math:`k+1` (:math:`c_{p}` is the heat
                0269 capacity). Averaging these terms over a layer thickness
                0270 :math:`\Delta{z_{k}}` (e.g., extending from the ice shelf base down to
                0271 the dashed line in cell C) and applying the averaged tendency to cell A
                0272 (in layer :math:`k`) and to the appropriate fraction of cells C (in
                0273 layer :math:`k+1`) yields
                0274 
                0275 .. math::
67dae37801 Eli *0276    g_{\theta,k}^*   = \frac{Q}{\rho_c c_{p} \Delta{z}_{k}}
274298f12b Jeff*0277    :label: tendencyk
                0278 
                0279 .. math::
                0280    g_{\theta,k+1}^*
67dae37801 Eli *0281    = \frac{Q}{\rho_c c_{p} \Delta{z}_{k}}
274298f12b Jeff*0282    \frac{ \Delta{z}_{k} ( 1- h_{k} )}{\Delta{z}_{k+1}}
                0283    :label: tendencykp1
                0284 
                0285 :eq:`tendencykp1` describes averaging over the part of the grid
                0286 cell :math:`k+1` that is part of the boundary layer with tendency
                0287 :math:`g_{\theta,k}^*` and the part with no tendency. Salinity is
                0288 treated in the same way. The momentum equations are not modified.
                0289 
0a8794f5ee Mart*0290 Three-equations thermodynamics
274298f12b Jeff*0291 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                0292 
67dae37801 Eli *0293 Freezing and melting form a boundary layer between the ice shelf and the
                0294 ocean that is represented in the model by an infinitesimal layer used to
                0295 calculate the exchanges between the ocean and the ice. Melting and
                0296 freezing at the boundary between saline water and ice imply a freshwater
                0297 mass flux :math:`q` (:math:`<0` for melting); the relevant heat fluxes
                0298 into and out of the boundary layer therefore include a diffusive flux
                0299 through the ice, the latent heat flux due to melting and freezing, and
                0300 the advective heat that is carried by the freshwater mass flux. There is
                0301 a salt flux carried by the mass flux if the ice has a non-zero salinity
274298f12b Jeff*0302 :math:`S_I`. Further, the position of the interface between ice and
                0303 ocean changes because of :math:`q`, so that, say, in the case of melting
67dae37801 Eli *0304 the volume of sea water increases. As a consequence of these fluxes,
                0305 salinity and temperature are affected.
                0306 
                0307 The turbulent tracer exchanges between the infinitesimal boundary layer and
                0308 the ocean are expressed as diffusive fluxes. Following Jenkins et
                0309 al. (2001) :cite:`jenkins:01`, the boundary conditions for a tracer take
                0310 into account that this boundary may not move with the ice ocean interface
                0311 (for example, in a linear free surface model). The implied upward freshwater
                0312 flux :math:`q` is therefore included in the boundary conditions for the
                0313 temperature and salinity equation as an advective flux.
                0314 
                0315 The boundary conditions for tracer :math:`X=S,T` (tracer :math:`X` stands
                0316 for either in-situ temperature :math:`T` or salinity :math:`S`, located at
                0317 the first interior ocean grid point) in the ocean are expressed as the sum
                0318 of advective and diffusive fluxes
274298f12b Jeff*0319 
                0320 .. math::
67dae37801 Eli *0321    F_X = (\rho_c \, \gamma_{X} -q ) ( X_{b} - X )
274298f12b Jeff*0322    :label: jenkinsbc
                0323 
67dae37801 Eli *0324 where the diffusive flux has been parameterized as a turbulent exchange
                0325 :math:`\rho_c \, \gamma_{X}( X_{b} - X )` following Holland and
                0326 Jenkins (1999) :cite:`holland:99` or Jenkins et al. (2001)
                0327 :cite:`jenkins:01`. :math:`X_b` indicates the tracer in the boundary layer,
                0328 :math:`\rho_c` the density of seawater (parameter :varlink:`rhoConst`), and
                0329 :math:`\gamma_X` is the turbulent exchange (or transfer) coefficient
                0330 (parameters :varlink:`SHELFICEheatTransCoeff` and
                0331 :varlink:`SHELFICEsaltTransCoeff`), in units of an exchange
                0332 velocity. In-situ temperature, computed locally from tracer potential
                0333 temperature, is required here to accurately compute the in-situ freezing
                0334 point of seawater in order to determine ice melt.
                0335 
                0336 The tracer budget for the infinitesimal boundary layer takes the general
                0337 form:
                0338 
                0339 .. math::
                0340    {\rho_I} K_{I,X} \frac{\partial{X_I}}{\partial{z}}\biggl|_{b}
                0341    = \rho_c \, \gamma_{X} ( X_{b} - X ) -  q ( X_{b} - X_{I} )
                0342    :label: jenkinsgenbudget
                0343 
                0344 where the LHS represents diffusive flux from the ice evaluated at the
                0345 interface between the infinitesimal boundary layer and the ice, and the RHS
                0346 represents the turbulent and advective exchanges between the infinitesimal
                0347 layer and the ocean and the advective exchange between the boundary layer
                0348 and the ice (:math:`qX_{I}`, this flux will be zero if the ice contains no
                0349 tracer: :math:`X_I=0`). The temperature in the boundary layer (:math:`T_b`)
                0350 is taken to be at the freezing point, which is a function of pressure and
                0351 salinity, :math:`\rho_I` is ice density (:varlink:`rhoShelfIce`), and
                0352 :math:`K_{I,X}` the appropriate ice diffusivity.
                0353 
                0354 For any material tracer such as salinity, the LHS in :eq:`jenkinsgenbudget`
                0355 vanishes (no material diffusion into the ice), while for temperature, the
                0356 term :math:`q\,( T_{b}-T_{I} )` vanishes, because both the boundary layer
                0357 and the ice are at the freezing point. Instead, the latent heat of freezing
                0358 is included as an additional term to take into account the conversion of
                0359 ice to water:
274298f12b Jeff*0360 
                0361 .. math::
67dae37801 Eli *0362    {\rho_I} \, c_{p,I} \, \kappa_{I,T}
                0363    \frac{\partial{T_I}}{\partial{z}}\biggl|_{b}
                0364    = c_{p} \, \rho_c \, \gamma_{T} ( T_{b} - T )+ L q.
                0365    :label: jenkinsheatbudget
                0366 
                0367 where :math:`\kappa_{I,T}` is the thermal ice diffusivity
                0368 (:varlink:`SHELFICEkappa`), :math:`L` is the latent heat of fusion
                0369 (:varlink:`SHELFICElatentHeat`), :math:`c_{p}` is the specific heat
                0370 capacity of water (:varlink:`HeatCapacity_Cp`), and :math:`c_{p,I}` the
                0371 heat capacity of the ice shelf (:varlink:`SHELFICEHeatCapacity_Cp`).  A
                0372 reasonable choice for :math:`\gamma_T` (:varlink:`SHELFICEheatTransCoeff`),
                0373 the turbulent exchange coefficient of temperature, is discussed in Holland
                0374 and Jenkins (1999) :cite:`holland:99` (see
                0375 :numref:`shelfice_exchange_coeff`).  The temperature at the interface
                0376 :math:`T_{b}` is assumed to be the in-situ freezing point temperature of
                0377 sea-water :math:`T_{f}`, which is computed from a linear equation of state:
274298f12b Jeff*0378 
                0379 .. math::
67dae37801 Eli *0380    T_{f} = (0.0901 - 0.0575\ S_{b})
                0381    - 7.61 \times 10^{-4} p_{b}.
0a8794f5ee Mart*0382    :label: hellmerfreeze
274298f12b Jeff*0383 
67dae37801 Eli *0384 where :math:`T_f` is given in :sup:`o`\ C and :math:`p_{b}` is in dBar. In
                0385 :eq:`jenkinsheatbudget`, the diffusive heat flux at the ice-ocean interface
                0386 can be appproximated by assuming a linear temperature profile in the ice
                0387 and approximating the vertical derivative of temperature in the ice as the
                0388 difference between the ice surface and ice bottom temperatures divided by
                0389 the ice thickness, so that the left-hand-side of :eq:`jenkinsheatbudget`
                0390 becomes
274298f12b Jeff*0391 
                0392 .. math::
67dae37801 Eli *0393    {\rho_I} \, c_{p,I} \, \kappa_{I,T}
                0394    \frac{\partial{T_I}}{\partial{z}}\biggl|_{b}
                0395    \approx \rho_{I} \, c_{p,I} \, \kappa_{I,T} \frac{(T_{S} - T_{b})}{h}
                0396    :label: dTdzdiffus
                0397 
                0398 where :math:`T_{S}` the (surface) temperature of the ice shelf
                0399 (:varlink:`SHELFICEthetaSurface`) and :math:`h` is the ice-shelf
                0400 draft. Alternatively, assuming that the ice is "advected" vertically as
                0401 implied by the meltflux :math:`q`, the diffusive flux can be approximated
                0402 as :math:`\min(q,0)\,c_{p,I} (T_{S} - T_{b})` (runtime flag
                0403 :varlink:`SHELFICEadvDiffHeatFlux`; see Holland and Jenkins, 1999
                0404 :cite:`holland:99` for details).
                0405 
                0406 From the salt budget, the salt flux across the shelf ice-ocean interface is
                0407 equal to the salt flux due to melting and freezing:
                0408 
                0409 .. math::
                0410    \rho_c \, \gamma_{S} (S - S_{b}) = - q\,(S_{b}-S_{I})
274298f12b Jeff*0411    :label: hellmersaltbalance
                0412 
67dae37801 Eli *0413 where :math:`\gamma_S =` :varlink:`SHELFICEsaltToHeatRatio` :math:`*
                0414 \gamma_T` is the turbulent salinity exchange coefficient.  Note, it is
                0415 assumed that :math:`\kappa_{I,S} =0`; moreover, the salinity of the ice
                0416 shelf is generally neglected (:math:`S_{I}=0`).
                0417 
c717173819 Matt*0418 The budget equations for temperature :eq:`jenkinsheatbudget` (with
                0419 :eq:`dTdzdiffus`) and salinity
67dae37801 Eli *0420 :eq:`hellmersaltbalance`, together with the freezing point temperature of
                0421 sea-water :eq:`hellmerfreeze`, form the so-called three-equation-model
                0422 (e.g., Hellmer and Olbers (1989) :cite:`hellmer:89`, Jenkins et al. (2001)
                0423 :cite:`jenkins:01`). These equations are solved to obtain :math:`S_b, T_b,
                0424 q`, which are then substituted into :eq:`jenkinsbc` to obtain the boundary
                0425 conditions for the temperature and salinity equations of the ocean
                0426 model. Note that with :math:`S_{I}=0` and :eq:`hellmersaltbalance`, the
                0427 boundary flux for salinity becomes :math:`F_S = q\,S`, which is the flux
                0428 that is necessary to account for the dilution of salinity in the case of
                0429 melting.
                0430 
                0431 The three-equation-model tends to yield smaller melt rates than the simpler
                0432 formulation of the ISOMIP protocol (:numref:`shelfice_isomip`) because the
                0433 freshwater flux (due to melting) decreases the salinity, which in turn
                0434 raises the freezing point temperature and thus leads to less melting at the
                0435 interface. For a simpler thermodynamics model where :math:`S_b` is not
                0436 computed explicitly, for example as in the ISOMIP protocol, :eq:`jenkinsbc`
                0437 cannot be applied directly. In this case :eq:`hellmersaltbalance` can be
                0438 used with :eq:`jenkinsbc` to obtain:
                0439 
                0440 .. math:: F_X = q\,(S-S_I)
274298f12b Jeff*0441 
                0442 This formulation can be used for all cases for which
67dae37801 Eli *0443 :eq:`hellmersaltbalance` is valid. Further, in this formulation it is
                0444 obvious that melting (:math:`q<0`) leads to a reduction of salinity.
274298f12b Jeff*0445 
67dae37801 Eli *0446 The default value of :varlink:`SHELFICEconserve` ``=.FALSE.`` removes the
                0447 contribution :math:`q\, ( X_{b}-X )` from :eq:`jenkinsbc`, making the
                0448 boundary conditions non-conservative.
274298f12b Jeff*0449 
0a8794f5ee Mart*0450 Solving the three-equations system
                0451 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                0452 
67dae37801 Eli *0453 There has been some confusion about the three-equations system, so we
                0454 document the solution in the code here: We use :eq:`hellmerfreeze`
                0455 :math:`T_{b} = a_{0} S_{b} + \epsilon_{4}` to eliminate :math:`T_{b}`
                0456 from :eq:`jenkinsheatbudget` with :eq:`dTdzdiffus` and find an
                0457 expression for the freshwater flux :math:`q`:
0a8794f5ee Mart*0458 
                0459 .. math::
                0460    \begin{aligned}
                0461    -Lq &= \epsilon_{1} (T - a_{0} S_{b} - \epsilon_{4})
                0462    + \epsilon_{3} (T_{S} - a_{0} S_{b} - \epsilon_{4}) \\
                0463    \Leftrightarrow Lq &=  a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}
                0464      + \epsilon_{q}
                0465    \end{aligned}
                0466    :label: solvedmeltrate
                0467 
                0468 to be substituted into :eq:`hellmersaltbalance`:
                0469 
                0470 .. math::
                0471    \begin{aligned}
                0472    \epsilon_{2}\,(S - S_{b}) &= - Lq\,(S_{b}-S_{I})
                0473    = - (a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}
                0474      + \epsilon_{q})\,(S_{b}-S_{I}) \\
                0475    \Leftrightarrow 0 &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}^{2}
                0476    + \{ \epsilon_{q}  - \epsilon_{2}
                0477      - a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{I} \}\,S_{b}
00f81e6785 Ou W*0478      + \epsilon_{2}\,S - \epsilon_{q}\,S_{I} \\
                0479    \Leftrightarrow 0 &= A\,S_{b}^{2} + B\,S_{b} + C \\
                0480    \Rightarrow S_{b} &= \frac{-B \pm \sqrt{ B^{2} - 4AC }}{2A}
                0481    \end{aligned}
                0482 
                0483 with the abbrevations
                0484 
                0485 .. math::
                0486    \begin{aligned}
                0487    \epsilon_{1} &= c_{p} \, \rho_c \, \gamma_{T}, \quad
                0488    \epsilon_{2} = \rho_c L \, \gamma_{S}, \quad
                0489    \epsilon_{3} = \frac{\rho_{I} \, c_{p,I} \, \kappa}{h}, \quad
                0490    \epsilon_{4} = b_{0}p + c_{0}, \\
                0491    \epsilon_{q} &= \epsilon_{1}\,(\epsilon_{4} - T)
                0492                   + \epsilon_{3}\,(\epsilon_{4} - T_{S}), \\
                0493    A &= a_{0}\,(\epsilon_{1} + \epsilon_{3}), \quad
                0494    B = \epsilon_{q} - \epsilon_{2} - a_{0}\,(\epsilon_{1} +
                0495    \epsilon_{3})\,S_{I}, \quad
                0496    C = \epsilon_{2}\,S -\epsilon_{q}\,S_{I}.
0a8794f5ee Mart*0497    \end{aligned}
                0498 
00f81e6785 Ou W*0499 The smaller non-negative root of the quadratic equation in :math:`S_{b}` is
                0500 used. By default, the ice shelf salinity :math:`S_{I}` is zero and the
                0501 quadratic equation simplifies to
0a8794f5ee Mart*0502 
                0503 .. math::
                0504    \begin{aligned}
                0505    0 &= a_{0}\,(\epsilon_{1} + \epsilon_{3})\,S_{b}^{2}
                0506    + (\epsilon_{q}  - \epsilon_{2}) \,S_{b} + \epsilon_{2}\,S \\
                0507      S_{b} &= \frac{\epsilon_{2} - \epsilon_{q}\mp
                0508      \sqrt{(\epsilon_{q}  - \epsilon_{2})^2
b8f9311cc6 Matt*0509      - 4\, a_{0}\,(\epsilon_{1} + \epsilon_{3})\,\epsilon_{2}\,S}}
0a8794f5ee Mart*0510      {2\,a_{0}\,(\epsilon_{1} + \epsilon_{3})}
                0511    \end{aligned}
                0512 
00f81e6785 Ou W*0513 With :math:`S_b`, the boundary layer temperature :math:`T_b` and the melt rate
                0514 :math:`q` are known through :eq:`hellmerfreeze` and :eq:`solvedmeltrate`.
0a8794f5ee Mart*0515 
67dae37801 Eli *0516 .. _shelfice_isomip:
                0517 
0a8794f5ee Mart*0518 ISOMIP thermodynamics
274298f12b Jeff*0519 ^^^^^^^^^^^^^^^^^^^^^
                0520 
                0521 A simpler formulation follows the ISOMIP protocol. The
                0522 freezing and melting in the boundary layer between ice shelf and ocean
                0523 is parameterized following Grosfeld et al. (1997) :cite:`grosfeld:97`. In this
67dae37801 Eli *0524 formulation :eq:`jenkinsheatbudget` reduces to
274298f12b Jeff*0525 
                0526 .. math::
67dae37801 Eli *0527    c_{p} \, \rho_c \, \gamma_T (T - T_{b})  = -L q
274298f12b Jeff*0528    :label: isomipheatbalance
                0529 
                0530 and the fresh water flux :math:`q` is computed from
                0531 
                0532 .. math::
67dae37801 Eli *0533    q = - \frac{c_{p} \, \rho_c \, \gamma_T (T - T_{b})}{L}
274298f12b Jeff*0534    :label: isomipfwflx
                0535 
67dae37801 Eli *0536 In order to use this formulation, set runtime parameter
                0537 :varlink:`useISOMIPTD` ``=.TRUE.`` in ``data.shelfice``.
                0538 
                0539 .. _shelfice_exchange_coeff:
274298f12b Jeff*0540 
0a8794f5ee Mart*0541 Exchange coefficients
                0542 ^^^^^^^^^^^^^^^^^^^^^
                0543 
                0544 The default exchange coefficents :math:`\gamma_{T/S}` are constant and
7b8b86ab99 Timo*0545 set by the run-time parameters :varlink:`SHELFICEheatTransCoeff` and
0a8794f5ee Mart*0546 :varlink:`SHELFICEsaltTransCoeff` (see
                0547 :numref:`tab_phys_pkg_shelfice_runtimeparms`). If
                0548 :varlink:`SHELFICEuseGammaFrict` ``=.TRUE.``, exchange coefficients
                0549 are computed from drag laws and friction velocities estimated from
                0550 ocean speeds following Holland and Jenkins (1999)
67dae37801 Eli *0551 :cite:`holland:99`. This computation can be modified using runtime
                0552 parameters and the user is referred to
0a8794f5ee Mart*0553 :filelink:`pkg/shelfice/shelfice_readparms.F` for details.
                0554 
274298f12b Jeff*0555 Remark
                0556 ^^^^^^
                0557 
                0558 The shelfice package and experiments demonstrating its strengths and
67dae37801 Eli *0559 weaknesses are also described in Losch (2008)
                0560 :cite:`losch:08`. Unfortunately, the description of the
                0561 thermodynamics in the appendix of Losch (2008) is wrong.
274298f12b Jeff*0562 
                0563 .. _shelfice_subroutines:
                0564 
                0565 Key subroutines
                0566 ~~~~~~~~~~~~~~~
                0567 
                0568 The main routine is :filelink:`shelfice_thermodynamics.F <pkg/shelfice/shelfice_thermodynamics.F>`
                0569 but note that :filelink:`/pkg/shelfice` routines are also called when solving the momentum equations.
                0570 
                0571 ::
                0572 
                0573     C     !CALLING SEQUENCE:
                0574     C ...
                0575     C |-FORWARD_STEP           :: Step forward a time-step ( AT LAST !!! )
                0576     C ...
                0577     C | |-DO_OCEANIC_PHY       :: Control oceanic physics and parameterization
                0578     C ...
                0579     C | | |-SHELFICE_THERMODYNAMICS :: main routine for thermodynamics
                0580     C                                  with diagnostics
                0581     C ...
                0582     C | |-THERMODYNAMICS       :: theta, salt + tracer equations driver.
                0583     C ...
                0584     C | | |-EXTERNAL_FORCING_T :: Problem specific forcing for temperature.
                0585     C | | |-SHELFICE_FORCING_T :: apply heat fluxes from ice shelf model
                0586     C ...
                0587     C | | |-EXTERNAL_FORCING_S :: Problem specific forcing for salinity.
                0588     C | | |-SHELFICE_FORCING_S :: apply fresh water fluxes from ice shelf model
                0589     C ...
                0590     C | |-DYNAMICS             :: Momentum equations driver.
                0591     C ...
                0592     C | | |-MOM_FLUXFORM       :: Flux form mom eqn. package ( see
                0593     C ...
                0594     C | | | |-SHELFICE_U_DRAG  :: apply drag along ice shelf to u-equation
                0595     C                             with diagnostics
                0596     C ...
                0597     C | | |-MOM_VECINV         :: Vector invariant form mom eqn. package ( see
                0598     C ...
                0599     C | | | |-SHELFICE_V_DRAG  :: apply drag along ice shelf to v-equation
                0600     C                             with diagnostics
                0601     C ...
                0602     C  o
                0603 
                0604 
                0605 
                0606 .. _shelfice_diagnostics:
                0607 
                0608 SHELFICE diagnostics
                0609 ~~~~~~~~~~~~~~~~~~~~
                0610 
                0611 Diagnostics output is available via the diagnostics package (see
                0612 :numref:`outp_pack`). Available output fields are summarized as follows:
                0613 
                0614 
                0615 ::
                0616 
                0617     ---------+----+----+----------------+-----------------
                0618      <-Name->|Levs|grid|<--  Units   -->|<- Tile (max=80c)
                0619     ---------+----+----+----------------+-----------------
                0620      SHIfwFlx|  1 |SM  |kg/m^2/s        |Ice shelf fresh water flux (positive upward)
                0621      SHIhtFlx|  1 |SM  |W/m^2           |Ice shelf heat flux  (positive upward)
                0622      SHIUDrag| 30 |UU  |m/s^2           |U momentum tendency from ice shelf drag
                0623      SHIVDrag| 30 |VV  |m/s^2           |V momentum tendency from ice shelf drag
                0624      SHIForcT|  1 |SM  |W/m^2           |Ice shelf forcing for theta, >0 increases theta
                0625      SHIForcS|  1 |SM  |g/m^2/s         |Ice shelf forcing for salt, >0 increases salt
                0626 
                0627 Experiments and tutorials that use shelfice
                0628 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                0629 
bd7056e9aa Jeff*0630 See the verification experiment :filelink:`isomip <verification/isomip>` for example usage of :filelink:`pkg/shelfice`.