Skip to content

Commit

Permalink
Merge branch 'master' into darwin
Browse files Browse the repository at this point in the history
  • Loading branch information
jahn committed Jul 27, 2024
2 parents 3a1d328 + 6519d60 commit a9e81b6
Show file tree
Hide file tree
Showing 119 changed files with 10,956 additions and 8,554 deletions.
101 changes: 77 additions & 24 deletions doc/phys_pkgs/obcs.rst
Original file line number Diff line number Diff line change
Expand Up @@ -548,10 +548,10 @@ OBCS\_APPLY\_*:
:filelink:`OBCS\_SPONGE <pkg/obcs/obcs_sponge.F>`:
##################################################

The sponge layer code (turned on with CPP option :varlink:`ALLOW_OBCS_SPONGE` and run-time parameter
:varlink:`useOBCSsponge`) adds a relaxation term to the right-hand-side of
the momentum and tracer equations. The variables are relaxed towards
the boundary values with a relaxation time scale that increases
The sponge layer code (turned on with CPP option :varlink:`ALLOW_OBCS_SPONGE`
and run-time parameter :varlink:`useOBCSsponge`) adds a relaxation term to the
right-hand-side of the momentum and tracer equations. The variables are relaxed
towards the boundary values with a relaxation time scale that increases
linearly with distance from the boundary

.. math::
Expand All @@ -563,27 +563,42 @@ linearly with distance from the boundary
{[(1-l)\tau_{b}+l\tau_{i}]}
where :math:`\chi` is the model variable (U/V/T/S) in the interior,
:math:`\chi_{BC}` the boundary value, :math:`L` the thickness of the
sponge layer (runtime parameter :varlink:`spongeThickness` in number
of grid points), :math:`\delta{L}\in[0,L]`
(:math:`\frac{\delta{L}}{L}=l\in[0,1]`) the distance from the boundary
(also in grid points), and :math:`\tau_{b}` (runtime parameters
:varlink:`Urelaxobcsbound` and :varlink:`Vrelaxobcsbound`) and
:math:`\chi_{BC}` the boundary value, :math:`L` the thickness of the sponge
layer (runtime parameter :varlink:`spongeThickness` in number of grid points),
:math:`\delta{L}\in[0,L]` (:math:`\frac{\delta{L}}{L}=l\in[0,1]`) the distance
from the boundary (also in grid points), and :math:`\tau_{b}` (runtime
parameters :varlink:`Urelaxobcsbound` and :varlink:`Vrelaxobcsbound`) and
:math:`\tau_{i}` (runtime parameters :varlink:`Urelaxobcsinner` and
:varlink:`Vrelaxobcsinner`) the relaxation time scales on the boundary
and at the interior termination of the sponge layer. The parameters
:varlink:`Urelaxobcsbound` and :varlink:`Urelaxobcsinner` set the relaxation time scales for
the Eastern and Western boundaries, :varlink:`Vrelaxobcsbound` and :varlink:`Vrelaxobcsinner`
for the Northern and Southern boundaries.
:varlink:`Vrelaxobcsinner`) the relaxation time scales on the boundary and at
the interior termination of the sponge layer. The parameters
:varlink:`Urelaxobcsbound` and :varlink:`Urelaxobcsinner` set the relaxation
time scales for the Eastern and Western boundaries, :varlink:`Vrelaxobcsbound`
and :varlink:`Vrelaxobcsinner` for the Northern and Southern boundaries.


OB's with nonlinear free surface
################################


OB's with sea ice
#################

Simple Dirichlet boundary conditions for sea ice parameters can be specified in
anology to the ocean variables via filenames ``OB[N/S/E/W][a/h/sl/sn/u/v]File``
(sea ice concentration, cell averaged sea ice thickness, salinity, cell
averaged snow thickness, ice drift components). With CPP-flag
:varlink:`ALLOW_OBCS_SEAICE_SPONGE` and runtime flags
:varlink:`useSeaiceSponge`, :varlink:`seaiceSpongeThickness`, and
``[A/H/SL/SN]relaxobcs[inner/bound]`` are available in analogy to the sponge
parameters for the ocean variables.

Neumann boundary conditions :math:`\frac{\partial\phi}{\partial{n}}=0` for all
sea ice variables can be applied with runtime flag
:varlink:`SEAICEuseNeumannBC`, which overrides the input files for the
Dirichlet values.

Defining CPP-flag :varlink:`OBCS_SEAICE_SMOOTH_EDGE` allows to smooth the
tracer sea-ice variables near the edges.


.. _ssub_phys_pkg_obcs_flowchart:

Expand All @@ -595,21 +610,59 @@ Flow chart


C !CALLING SEQUENCE:
c ...
C [...]
C | |-MAIN_DO_LOOP :: Open-AD case: Main timestepping loop routine
C | \ otherwise: just call FORWARD_STEP
C | |
C/\ | |-FORWARD_STEP :: Step forward a time-step ( AT LAST !!! )
C [...]
C/\ | | |-DO_OCEANIC_PHYS :: Oceanic (& seaice) physics computation
C/\ | | | |
C/\ | | | |-OBCS_CALC :: Open boundary. package (see pkg/obcs).
C/\ | | | |
C [...]
C/\ | | | |-SEAICE_MODEL :: pkg/seaice
C/\ | | | | |-SEAICE_DYNSOLVER :: pkg/seaice
C/\ | | | | | |-OBCS_APPLY_UVICE :: apply uIce/vIce boudnary conditions
C/\ | | | | |-OBCS_ADJUST_UVICE :: (Only for OBCS_UVICE_OLD)
C/\ | | | | |-SEAICE_GROWTH
C/\ | | | | |-SEAICE_APPLY_SEAICE :: add OBCS for scalar variables
C [...]
C/\ | | |-THERMODYNAMICS :: theta, salt + tracer equations driver.
C/\ | | | | (synchronous time-stepping case)
C [...]
C/\ | | | |-TEMP_INTEGRATE :: Step forward Prognostic Eq for Temperature.
C/\ | | | |
C/\ | | | |-SALT_INTEGRATE :: Step forward Prognostic Eq for Salinity.
C/\ | | | | same sequence of calls as in TEMP_INTEGRATE
C/\ | | | |
C/\ | | | |-PTRACERS_INTEGRATE :: Integrate other tracer(s) (see pkg/ptracers).
C/\ | | | | | same sequence of calls as in TEMP_INTEGRATE
C/\ | | | | |-OBCS_APPLY_PTRACER :: Open boundary package for pTracers
C/\ | | | |
C/\ | | | |-OBCS_APPLY_TS :: Open boundary package (see pkg/obcs ).
C/\ | | |
C [...]
C/\ | | |
C/\ | | |-DYNAMICS :: Momentum equations driver.
C/\ | | | |
C [...]
C/\ | | | |-OBCS_APPLY_UV :: Apply Open bndary Conditions to provisional U,V
C [...]
C/\ | | |-MOMENTUM_CORRECTION_STEP :: Finalise momentum stepping
C [...]
C/\ | | | |-OBCS_APPLY_UV :: Open boundary package (see pkg/obcs).



.. _ssub_phys_pkg_obcs_diagnostics:

OBCS diagnostics
++++++++++++++++

Diagnostics output is available via the diagnostics package (see :numref:`sub_outp_pkg_diagnostics`). Available output fields are summarized below:

::

------------------------------------------------------
<-Name->|Levs|grid|<-- Units -->|<- Tile (max=80c)
------------------------------------------------------
Diagnostics output is available via the diagnostics package (see
:numref:`sub_outp_pkg_diagnostics`). Currently there are no OBCS-specific
diagnostics available.


.. _ssub_phys_pkg_obcs_experiments:
Expand Down
102 changes: 60 additions & 42 deletions doc/phys_pkgs/seaice.rst
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ automatically undefines more recent features, see :filelink:`SEAICE_OPTIONS.h
:varlink:`SEAICE_ALLOW_FREEDRIFT`, #undef, enable solve approximate sea ice momentum equation and bypass solving for sea ice internal stress
:varlink:`SEAICE_EXTERNAL_FLUXES`, #define, use :filelink:`pkg/exf`-computed fluxes as starting point
:varlink:`SEAICE_ZETA_SMOOTHREG`, #define, use differentiable regularization for viscosities
:varlink:`SEAICE_DELTA_SMOOTHREG`, #undef, use differentiable regularization for :math:`1/\Delta`
:varlink:`SEAICE_DELTA_SMOOTHREG`, #undef, use differentiable regularization :math:`\Delta_{\mathrm{reg}}=\sqrt{\Delta^2+\Delta_{\min}}` instead of :math:`\max`-function for :math:`1/\Delta_{\mathrm{reg}}`
:varlink:`SEAICE_ALLOW_BOTTOMDRAG`, #undef, enable grounding parameterization for improved fastice in shallow seas
:varlink:`SEAICE_BGRID_DYNAMICS`, #undef, use sea ice dynamics code on legacy B-grid; most of the previous flags are not available with B-grid
:varlink:`SEAICE_BICE_STRESS`, #undef, B-grid only for backward compatiblity: turn on ice-stress on ocean; defined by default if :varlink:`SEAICE_BGRID_DYNAMICS` is defined
Expand Down Expand Up @@ -321,6 +321,8 @@ General flags and parameters
+------------------------------------+------------------------------+-------------------------------------------------------------------------+
| :varlink:`SEAICE_EPS` | 1.0E-10 | a "small number" used in various routines |
+------------------------------------+------------------------------+-------------------------------------------------------------------------+
| :varlink:`SEAICE_deltaMin` | :varlink:`SEAICE_EPS` | minimum to regularize :math:`\Delta` |
+------------------------------------+------------------------------+-------------------------------------------------------------------------+
| :varlink:`SEAICE_area_reg` | 1.0E-5 | minimum concentration to regularize ice thickness |
+------------------------------------+------------------------------+-------------------------------------------------------------------------+
| :varlink:`SEAICE_hice_reg` | 0.05 | minimum ice thickness (m) for regularization |
Expand Down Expand Up @@ -501,21 +503,27 @@ The momentum equation of the sea-ice model is
- m \nabla{\phi(0)} + \mathbf{F}
:label: eq_momseaice
where :math:`m=m_{i}+m_{s}` is the ice and snow mass per unit area;
:math:`\mathbf{u}=u\hat{\mathbf{i}}+v\hat{\mathbf{j}}` is the ice velocity vector;
:math:`\hat{\mathbf{i}}`, :math:`\hat{\mathbf{j}}`, and :math:`\hat{\mathbf{k}}` are unit vectors
in the :math:`x`, :math:`y`, and :math:`z` directions, respectively; :math:`f`
is the Coriolis parameter; :math:`\mathbf{\tau}_\mathrm{air}` and
:math:`\mathbf{\tau}_\mathrm{ocean}` are the wind-ice and ocean-ice stresses,
respectively; :math:`g` is the gravity accelation; :math:`\nabla\phi(0)` is the
gradient (or tilt) of the sea surface height; :math:`\phi(0) = g\eta +
p_{a}/\rho_{0} + mg/\rho_{0}` is the sea surface height potential in response
to ocean dynamics (:math:`g\eta`), to atmospheric pressure loading
(:math:`p_{a}/\rho_{0}`, where :math:`\rho_{0}` is a reference density) and a
term due to snow and ice loading ; and :math:`\mathbf{F}= \nabla \cdot\sigma` is
the divergence of the internal ice stress tensor :math:`\sigma_{ij}`.
Advection of sea-ice momentum is neglected. The wind and ice-ocean stress terms
are given by
where :math:`m=m_{i}+m_{s}` is the ice and snow mass per unit area. The ice
mass per grid cell is :math:`m_i=\rho_{\mathrm{ice}} h\,c` with the mean ice
density :math:`\rho_{\mathrm{ice}}` and the mean thickness :math:`h\,c = `
volume per grid cell area that is the product of the actual thickness :math:`h`
of the ice covered part of the cell and the fractional ice cover :math:`c =
[0,1]`, sloppily also called ice concentration. A similar relationship defines
the snow mass per grid cell :math:`m_s`.
:math:`\mathbf{u}=u\hat{\mathbf{i}}+v\hat{\mathbf{j}}` is the ice velocity
vector; :math:`\hat{\mathbf{i}}`, :math:`\hat{\mathbf{j}}`, and
:math:`\hat{\mathbf{k}}` are unit vectors in the :math:`x`, :math:`y`, and
:math:`z` directions, respectively; :math:`f` is the Coriolis parameter;
:math:`\mathbf{\tau}_\mathrm{air}` and :math:`\mathbf{\tau}_\mathrm{ocean}` are
the wind-ice and ocean-ice stresses, respectively; :math:`g` is the gravity
accelation; :math:`\nabla\phi(0)` is the gradient (or tilt) of the sea surface
height; :math:`\phi(0) = g\eta + p_{a}/\rho_{0} + mg/\rho_{0}` is the sea
surface height potential in response to ocean dynamics (:math:`g\eta`),
atmospheric pressure loading (:math:`p_{a}/\rho_{0}`, where :math:`\rho_{0}` is
a reference density), and a term due to snow and ice loading; and
:math:`\mathbf{F}= \nabla \cdot\sigma` is the divergence of the internal ice
stress tensor :math:`\sigma_{ij}`. Advection of sea-ice momentum is
neglected. The wind and ice-ocean stress terms are given by

.. math::
\begin{aligned}
Expand Down Expand Up @@ -568,21 +576,27 @@ compactness (concentration) :math:`c`:
with the constants :math:`P^{\ast}` (run-time parameter
:varlink:`SEAICE_strength`) and :math:`C^{\ast}` (run-time parameter
:varlink:`SEAICE_cStar`). By default, :math:`P` (variable :varlink:`PRESS` in
the code) is the replacement pressure
:varlink:`SEAICE_cStar`). Note that Hibler (1979) :cite:`hibler:79` defines
:math:`h` as the "mean thickness" or an "equivalent ice thickness" for mass,
which is :math:`c\,h` with our definitions. By default, :math:`P` (variable
:varlink:`PRESS` in the code) is the replacement pressure

.. math::
:label: eq_pressrepl
P = (1-k_t)\,P_{\max} \left( (1 - f_{r})
+ f_{r} \frac{\Delta}{\Delta_{\rm reg}} \right)
where :math:`f_{r}` is run-time parameter :varlink:`SEAICEpressReplFac`
where :math:`f_{r}` is a run-time parameter :varlink:`SEAICEpressReplFac`
(default = 1.0), and :math:`\Delta_{\rm reg}` is a regularized form of
:math:`\Delta = \left[ \left(\dot{\epsilon}_{11}+\dot{\epsilon}_{22}\right)^2 +
e^{-2}\left( \left(\dot{\epsilon}_{11}-\dot{\epsilon}_{22} \right)^2 +
\dot{\epsilon}_{12}^2 \right) \right]^{\frac{1}{2}}`, for example
:math:`\Delta_{\rm reg} = \max(\Delta,\Delta_{\min})`.
4\,\dot{\epsilon}_{12}^2 \right) \right]^{\frac{1}{2}}`. By default
:math:`\Delta_{\mathrm{reg}}=\max(\Delta,\Delta_{\min})`. If CPP-flag
:varlink:`SEAICE_DELTA_SMOOTHREG` is defined,
:math:`\Delta_{\mathrm{reg}}=\sqrt{\Delta^2+\Delta^2_{\min}}`. Run-time
parameter :varlink:`SEAICE_deltaMin` :math:`= \Delta_{\min} = 10^{-10}` by
default.

The tensile strength factor :math:`k_t` (run-time parameter
:varlink:`SEAICE_tensilFac`) determines the ice tensile strength :math:`T =
Expand Down Expand Up @@ -656,7 +670,7 @@ with the ratio of major to minor axis :math:`e = 2.0` (run-time parameter

.. math::
\begin{aligned}
\zeta =& \min\left(\frac{(1+k_t)P_{\max}}{2\max(\Delta,\Delta_{\min})},
\zeta =& \min\left(\frac{(1+k_t)P_{\max}}{2\Delta_\mathrm{reg}},
\zeta_{\max}\right) \\
\eta =& \frac{\zeta}{e^2}
\end{aligned}
Expand All @@ -669,16 +683,18 @@ with the abbreviation
\Delta = \left[
\left(\dot{\epsilon}_{11}+\dot{\epsilon}_{22}\right)^2
+ e^{-2}\left( \left(\dot{\epsilon}_{11}-\dot{\epsilon}_{22} \right)^2
+ \dot{\epsilon}_{12}^2 \right)
+ 4\,\dot{\epsilon}_{12}^2 \right)
\right]^{\frac{1}{2}}
The bulk viscosities are bounded above by imposing both a minimum
:math:`\Delta_{\min}` (for numerical reasons, run-time parameter
:math:`\Delta_{\min}` and replacing :math:`\Delta` by the regularized version
:math:`\Delta_\mathrm{reg}` (for historical reasons, run-time parameter
:varlink:`SEAICE_deltaMin` is set to a default value of
:math:`10^{-10}\,\text{s}^{-1}`, the value of :varlink:`SEAICE_EPS`) and a
maximum :math:`\zeta_{\max} = P_{\max}/(2\Delta^\ast)`, where
:math:`\Delta^\ast=(2\times10^4/5\times10^{12})\,\text{s}^{-1}` :math:`=
2\times10^{-9}\,\text{s}^{-1}`. Obviously, this corresponds to regularizing
:math:`\Delta^\ast=(2\times10^4/5\times10^{12})\,\text{s}^{-1} =
2\times10^{-9}\,\text{s}^{-1}` (:varlink:`SEAICE_zetaMaxFac`
:math:`=\frac{1}{2\Delta^\ast}`). Obviously, this corresponds to regularizing
:math:`\Delta` with the typical value of :varlink:`SEAICE_deltaMin` :math:`=
2\times10^{-9}`. Clearly, some of this regularization is redundant. (There is
also the option of bounding :math:`\zeta` from below by setting run-time
Expand All @@ -695,9 +711,9 @@ expression:
.. math::
\begin{split}
\zeta &= \zeta_{\max}\tanh\left(\frac{(1+k_t)P_{\max}}{2\,
\min(\Delta,\Delta_{\min}) \,\zeta_{\max}}\right)\\
\Delta_\mathrm{reg} \,\zeta_{\max}}\right)\\
&= \frac{(1+k_t)P_{\max}}{2\Delta^\ast}
\tanh\left(\frac{\Delta^\ast}{\min(\Delta,\Delta_{\min})}\right)
\tanh\left(\frac{\Delta^\ast}{\Delta_\mathrm{reg}}\right)
\end{split}
:label: eq_zetaregsmooth
Expand Down Expand Up @@ -737,7 +753,7 @@ with the abbreviation
.. math::
\Delta = \sqrt{(\dot{\epsilon}_{11}-\dot{\epsilon}_{22})^2
+\frac{e_F^2}{e_G^4}((\dot{\epsilon}_{11}
-\dot{\epsilon}_{22})^2+4\dot{\epsilon}_{12}^2)}.
-\dot{\epsilon}_{22})^2+4\,\dot{\epsilon}_{12}^2)}.
Note that if :math:`e_G=e_F=e`, these formulae reduce to the normal flow rule.

Expand Down Expand Up @@ -1468,7 +1484,7 @@ In the zero-layer model of Semtner (1976) :cite:`semtner:76`, the conductive
heat flux depends strongly on the ice thickness :math:`h`. However, the ice
thickness in the model represents a mean over a potentially very heterogeneous
thickness distribution. In order to parameterize a sub-grid scale distribution
for heat flux computations, the mean ice thickness :math:`h` is split into
for heat flux computations, the ice thickness :math:`h` is split into
:math:`N` thickness categories :math:`H_{n}` that are equally distributed
between :math:`2h` and a minimum imposed ice thickness of :math:`5\,\text{cm}`
by :math:`H_n= \frac{2n-1}{7}\,h` for :math:`n\in[1,N]`. The heat fluxes
Expand Down Expand Up @@ -1519,25 +1535,27 @@ parameter :varlink:`SEAICEuseFlooding` set to ``.TRUE.``.
Advection of thermodynamic variables
------------------------------------

Effective ice thickness (ice volume per unit area, :math:`c h`),
concentration :math:`c` and effective snow thickness (:math:`c h_s`)
are advected by ice velocities:
Mean ice thickness (ice volume per unit area, :math:`c h`, model variable
:varlink:`HEFF`, which implies the misleading name "effective thickness"),
concentration :math:`c` (model variable :varlink:`AREA`) and mean snow
thickness (:math:`c h_s`, model variable :varlink:`HSNOW`) are advected by ice
velocities:

.. math::
\frac{\partial{X}}{\partial{t}} =
- \nabla \cdot\left(\mathbf{u}\,X\right) + \Gamma_{X} + D_{X}
:label: eq_advection
where :math:`\Gamma_X` are the thermodynamic source terms and :math:`D_{X}` the
diffusive terms for quantities :math:`X= c h, c, c h_s`. From
the various advection schemes that are available in MITgcm, we recommend
flux-limited schemes to preserve sharp gradients and edges that are typical of
sea ice distributions and to rule out unphysical over- and undershoots
(negative thickness or concentration). These schemes conserve volume and
horizontal area and are unconditionally stable, so that we can set
:math:`D_{X}=0`. Run-time flags: :varlink:`SEAICEadvScheme` (default=77, is a
2nd-order flux limited scheme), :varlink:`DIFF1` = :math:`D_{X}/\Delta{x}`
(default=0).
diffusive terms for quantities :math:`X= c h, c, c h_s` or any other tracer,
such as sea ice salinity. From the various advection schemes that are available
in MITgcm, we recommend flux-limited schemes (runtime flag
:varlink:`SEAICEadvScheme`; default=77, a 2nd-order flux limited scheme) to
preserve sharp gradients and edges that are typical of sea ice distributions
and to rule out unphysical over- and undershoots (negative thickness or
concentration). These schemes conserve volume and horizontal area and are
unconditionally stable, so that we can set :math:`D_{X}=0` (runtime flag
:varlink:`DIFF1` = :math:`D_{X}/\Delta{x}`; default=0).

The MITgcm sea ice model provides the option to use the thermodynamics model of
Winton (2000) :cite:`winton:00`, which in turn is based on the 3-layer model of
Expand Down
24 changes: 24 additions & 0 deletions doc/tag-index
Original file line number Diff line number Diff line change
@@ -1,6 +1,30 @@
Notes on tags used in MITgcmUV
==============================

checkpoint68z (2024/07/27)
o pkg/diagnostics:
- to allow to use ADJ-diags in DIVA runs, change call sequence regarding
TURNOFF_MODEL_IO and DIAGSTATS_INI_IO and simplify (fewer tests for
"costfinal") the_model_main.F (+ shorten pkg/openad version) ;
- restore the use of "diag_pkgStatus" (switch & check) in AD backward sweep
with new small S/R for DIVA runs ;
- turn on FWD and ADJ-diags in secondary lab_sea AD test "noseaice".
o pkg/bling:
- clean-up "bling_light.F" (fix typo related to PHYTO_SELF_SHADING option,
move k-loop outside) ; add CHL to bling_ad_check_lev{2,3,4}_dir.h and
fix TAF storage directives ;
- change "global_oce_biogeo_bling" AD test exp. (both TAF & Tapenade) to
use PHYTO_SELF_SHADING (with #undef USE_QSW) and update ref. output.
- also clean-up MDSIO_BUFF_WH.h (avoid common block that mixes r4 & r8 vars).
o pkg/obcs:
- Neumann boundary conditions for sea ice variables, activate by setting
new runtime parameter useSeaiceNeumann=.TRUE. in data.obcs ;
- fewer (hidden) recomputations and no recomputation warnings when
both ALLOW_OBCS_BALANCE and ALLOW_OBCS_STEVENS are defined.
o pkg/cost:
- new runtime flag for pkg/cost "cost_mask_file", which can be used
for flexible cost function definition.

checkpoint68y (2024/06/05)
o pkg/bling:
- fix bigR setting when option CARBONCHEM_TOTALPHSCALE is defined.
Expand Down
Loading

0 comments on commit a9e81b6

Please sign in to comment.