From 2a9941171a6b405235a4126c9b09fd80151d8fa9 Mon Sep 17 00:00:00 2001
From: Andrew Coward <acc@noc.ac.uk>
Date: Mon, 11 Nov 2024 13:34:36 +0000
Subject: [PATCH] Additions and corrections to the new chap_D2D.tex file
---
latex/NEMO/subfiles/chap_D2D.tex | 278 ++++++++++++-------------------
1 file changed, 106 insertions(+), 172 deletions(-)
diff --git a/latex/NEMO/subfiles/chap_D2D.tex b/latex/NEMO/subfiles/chap_D2D.tex
index b9f056f..ca4e727 100644
--- a/latex/NEMO/subfiles/chap_D2D.tex
+++ b/latex/NEMO/subfiles/chap_D2D.tex
@@ -50,35 +50,15 @@ The external forcings and parameterisations require complex inputs
that are carried out in modules SBC, LDF and ZDF and are described in
\autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.
-In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence,
-curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module).
-
-The different options available to the user are managed by namelist variables.
-For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx},
-where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.
-%If a CPP key is used for this term its name is \key{ttt}.
-The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory,
-and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine.
-
-The user has the option of extracting and outputting each tendency term from the 3D momentum equations
-(\texttt{trddyn?} defined), as described in \autoref{chap:MISC}.
-Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined)
-can be derived from the 3D terms.
-\cmtgm{STEVEN: not quite sure I've got the sense of the last sentence.
- Does MISC correspond to "extracting tendency terms" or "vorticity balance"?}
+In the present chapter we describe the equations used to compute the barotropic velocities and
+the sea surface height that provides the surface pressure gradient.
%% =================================================================================================
-\section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})}
+\section[Surface pressure gradient (\textit{dynspg.F90}/\textit{stp2d.F90})]
+{Surface pressure gradient (\protect\mdl{dynspg}/\protect\mdl{stp2d})}
\label{sec:D2D_spg}
-\begin{listing}
- \nlst{namdyn_spg}
- \caption{\forcode{&namdyn_spg}}
- \label{lst:namdyn_spg}
-\end{listing}
-
-Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables. The surface
-pressure gradient term is related to the representation of the free surface
+The surface pressure gradient term is related to the representation of the free surface
(\autoref{sec:MB_hor_pg}). The main distinction is between the fixed volume case (linear
free surface, \ie\ with the inclusion of \key{linssh}) and the variable volume case
(nonlinear free surface, \key{qco}). In the linear free surface case
@@ -86,20 +66,28 @@ free surface, \ie\ with the inclusion of \key{linssh}) and the variable volume c
while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}).
With both linear and nonlinear free surface, external gravity waves are allowed in the
equations, which imposes a very small time step when an explicit time stepping is used
-(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). To allow a longer time step for the
-three-dimensional equations, one can use a split-explicit free surface
-(\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}). In that case, only a quasi-linear
+(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp} (MLF time-stepping only)). To allow a
+longer time step for the three-dimensional equations, one can use a split-explicit free
+surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}). In that case, only a quasi-linear
form of 2d barotropic equations is substepped with a small time increment.
+Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables.
+
+\begin{listing}
+ \nlst{namdyn_spg}
+ \caption{\forcode{&namdyn_spg}}
+ \label{lst:namdyn_spg}
+\end{listing}
%% =================================================================================================
-\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
+\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]
+{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
\label{subsec:D2D_spg_exp}
-In the explicit free surface formulation (\np{ln_dynspg_exp}{ln\_dynspg\_exp} set to true),
-the model time step is chosen to be small enough to resolve the external gravity waves
-(typically a few tens of seconds).
-The surface pressure gradient, evaluated using a leap-frog scheme (\ie\ centered in time),
-is thus simply given by :
+With MLF time-stepping there is an explicit free surface formulation
+(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}), for which the model time step must be
+chosen to be small enough to resolve the external gravity waves (typically a few tens of
+seconds). The surface pressure gradient, evaluated using a leap-frog scheme (\ie\
+centered in time), is thus simply given by :
\begin{equation}
\label{eq:DYN_spg_exp}
\left\{
@@ -110,7 +98,7 @@ is thus simply given by :
\right.
\end{equation}
-Note that this option is not yet compatible with RK3 time-stepping.
+Note that this option is not yet available with RK3 time-stepping.
%% =================================================================================================
\subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})}
@@ -149,11 +137,44 @@ where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containi
coupling term between modes, surface atmospheric forcing as well as slowly varying
barotropic terms not explicitly computed to gain efficiency. The third term on the right
hand side of \autoref{eq:D2D_BT_dyn} represents the bottom stress (see section
-\autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration. Temporal
-discretization of the system above follows a three-time step Generalized Forward Backward
-algorithm detailed in \citet{shchepetkin.mcwilliams_OM05}. AB3-AM4 coefficients used in
-\NEMO\ follow the second-order accurate, "multi-purpose" stability compromise as defined
-in \citet{shchepetkin.mcwilliams_ibk09} (see their figure 12, lower left).
+\autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.
+
+For 'single 1$^{\mathrm st}$', RK3 time-stepping, barotropic velocities and sea surface
+height are calculated for the {\textit after} timestep from the {\textit before} values
+prior to the first stage of the main RK3 timestepping. This occurs in the new module
+\forcode{stp2d.F90} which computes the RHS forcing terms from {\textit before} fields and
+then calls the \forcode{dyn_spg_ts} routine to perform the time-stepping as detailed
+below. In the MLF time-stepping case, RHS forcing terms are computed inside
+\forcode{dyn_spg_ts} itself and the time index used will depend on options chosen in
+\forcode{namspgts}.
+
+For clarity, in RK3 the RHS forcing terms on momentum are vertically averaged 3D trends of:
+\[
+ \text{HPG} + \text{LDF} + \text {(COR + RVO)} + \text{KEG} + \text{ZAD}
+\]
+for the Vector invariant formulation (for which the 3D fields are also the RHS terms for
+the 1$^{\mathrm st}$ stage RK3 time stepping) and:
+\[
+ \text{HPG} + \text{LDF} + \text {(COR + MET)} + \text{ADV}
+\]
+for the flux form. In the latter case ADV must be recomputed at the 1$^{\mathrm st}$ stage.
+The new terms in these summaries are: RVO - Relative Vorticity and MET - Metric term.
+External forcing is also applied in the form of baroclinic bottom drag and surface winds.
+Other external forcings may be optionally applied as detailed in
+\autoref{subsec:D2D_forcing}.
+
+For ssh, the external forcing is the net column-average freshwater flux. This is
+evaporation minus prepcitation plus optional contributions from other sources as listed in
+\autoref{subsec:D2D_forcing}.
+
+%% =================================================================================================
+\subsection{Split-explicit time-stepping (\protect\mdl{dynspg\_ts})}
+\label{subsec:D2D_spg_ts_ts}
+
+Temporal discretization of the system above follows a three-time step Generalized Forward
+Backward algorithm detailed in \citet{shchepetkin.mcwilliams_OM05}. AB3-AM4 coefficients
+used in \NEMO\ follow the second-order accurate, "multi-purpose" stability compromise as
+defined in \citet{shchepetkin.mcwilliams_ibk09} (see their figure 12, lower left).
\begin{figure}[!t]
\centering
@@ -194,8 +215,8 @@ barotropic equations starting from \textit{before} time step
\np{nn_e}{nn\_e} additional iterations are indeed necessary), the baroclinic to barotropic
forcing term given at \textit{now} time step become centred in the middle of the
integration window. It can easily be shown that this property removes part of splitting
-errors between modes, which increases the overall numerical robustness.
-%references to Patrick Marsaleix' work here. Also work done by SHOM group.
+errors between modes, which increases the overall numerical robustness. %references to
+Patrick Marsaleix' work here. Also work done by SHOM group.
As far as tracer conservation is concerned, barotropic velocities used to advect tracers
@@ -209,151 +230,64 @@ is here the key to obtain exact conservation.
One can eventually choose to feedback instantaneous values by not using any time filter
-(\protect\np[=3]{nn_bt_flt}{nn\_bt\_flt}).
-In that case, external mode equations are continuous in time,
-\ie\ they are not re-initialized when starting a new sub-stepping sequence.
-This is the method used in the POM model for example, the stability being maintained by
-refreshing at (almost) each barotropic time step advection and horizontal diffusion terms.
-Since the latter terms have not been added in \NEMO\ for computational efficiency,
-removing time filtering would be inevitably unstable. One can however add some dissipation, but in the time domain, by slightly modifying the barotropic time stepping coefficients (\citet{demange_JCP19}). This is implemented here through an additional parameter (\np{rn_bt_alpha}{rn\_bt\_alpha}), which controls the amount of temporal diffusion.
+(\protect\np[=3]{nn_bt_flt}{nn\_bt\_flt}). In that case, external mode equations are
+continuous in time, \ie\ they are not re-initialized when starting a new sub-stepping
+sequence. This is the method used in the POM model for example, the stability being
+maintained by refreshing at (almost) each barotropic time step advection and horizontal
+diffusion terms. Since the latter terms have not been added in \NEMO\ for computational
+efficiency, removing time filtering would be inevitably unstable. One can however add some
+dissipation, but in the time domain, by slightly modifying the barotropic time stepping
+coefficients (\citet{demange_JCP19}). This is implemented here through an additional
+parameter (\np{rn_bt_alpha}{rn\_bt\_alpha}), which controls the amount of temporal
+diffusion.
+%% =================================================================================================
+\subsection{External forcings}
+\label{subsec:D2D_forcing}
-\cmtgm{ %%% copy from griffies Book
-
-\textbf{title: Time stepping the barotropic system }
-
-Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$.
-Hence, we can update the surface height and vertically integrated velocity with a leap-frog scheme using
-the small barotropic time step $\rdt$.
-We have
+Besides the surface and bottom baroclinic stresses four other forcings may enter the barotropic
+momentum equations.
-\[
- % \label{eq:DYN_spg_ts_eta}
- \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1})
- = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right]
-\]
-\begin{multline*}
- % \label{eq:DYN_spg_ts_u}
- \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}) \\
- = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n})
- - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right]
-\end{multline*}
-\
-
-In these equations, araised (b) denotes values of surface height and vertically integrated velocity updated with
-the barotropic time steps.
-The $\tau$ time label on $\eta^{(b)}$ and $U^{(b)}$ denotes the baroclinic time at which
-the vertically integrated forcing $\textbf{M}(\tau)$
-(note that this forcing includes the surface freshwater forcing),
-the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$,
-and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over
-a single cycle.
-This is also the time that sets the barotropic time steps via
-\[
- % \label{eq:DYN_spg_ts_t}
- t_n=\tau+n\rdt
-\]
-with $n$ an integer.
-The density scaled surface pressure is evaluated via
-\[
- % \label{eq:DYN_spg_ts_ps}
- p_s^{(b)}(\tau,t_{n}) =
- \begin{cases}
- g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o & \text{non-linear case} \\
- g \;\eta_s^{(b)}(\tau,t_{n}) & \text{linear case}
- \end{cases}
-\]
-To get started, we assume the following initial conditions
-\[
- % \label{eq:DYN_spg_ts_eta}
- \begin{split}
- \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} \\
- \eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}
- \end{split}
-\]
-with
-\[
- % \label{eq:DYN_spg_ts_etaF}
- \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n})
-\]
-the time averaged surface height taken from the previous barotropic cycle.
-Likewise,
-\[
- % \label{eq:DYN_spg_ts_u}
- \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ \\
- \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}
-\]
-with
-\[
- % \label{eq:DYN_spg_ts_u}
- \overline{\textbf{U}^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n})
-\]
-the time averaged vertically integrated transport.
-Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration.
-
-Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ ,
-the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at
-baroclinic time $\tau + \rdt \tau$
-\[
- % \label{eq:DYN_spg_ts_u}
- \textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n})
-\]
-The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using
-the following form
+(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
+the atmospheric pressure is taken into account when computing the surface pressure gradient.
-\begin{equation}
- \label{eq:DYN_spg_ts_ssh}
- \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right]
-\end{equation}
+(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see
+\autoref{sec:SBC_TDE}), the tidal potential is taken into account when computing the
+surface pressure gradient.
-The use of this "big-leap-frog" scheme for the surface height ensures compatibility between
-the mass/volume budgets and the tracer budgets.
-More discussion of this point is provided in Chapter 10 (see in particular Section 10.2).
+(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and SI3 is used (\ie\ when the sea-ice is
+embedded in the ocean), the snow-ice mass is taken into account when computing the surface
+pressure gradient.
-In general, some form of time filter is needed to maintain integrity of the surface height field due to
-the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
-We have tried various forms of such filtering,
-with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to
-its stability and reasonably good maintenance of tracer conservation properties (see ??).
+(4) When \np[=.true.]{ln_wave .AND. ln_bern_srfc}{ln\_wave .AND. ln\_bern\_srfc} (see
+\autoref{subsec:SBC_wave_bhd}) , a depth-uniform wave-induced kinematic pressure term (the
+Bernoulli head) is added to the mean pressure.
-\begin{equation}
- \label{eq:DYN_spg_ts_sshf}
- \eta^{F}(\tau-\Delta) = \overline{\eta^{(b)}(\tau)}
-\end{equation}
-Another approach tried was
+\vskip 0.5cm
-\[
- % \label{eq:DYN_spg_ts_sshf2}
- \eta^{F}(\tau-\Delta) = \eta(\tau)
- + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt)
- + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right]
-\]
+\noindent The net column-average freshwater flux applied to the ssh equation can also be modified by
+the following options:
-which is useful since it isolates all the time filtering aspects into the term multiplied by $\alpha$.
-This isolation allows for an easy check that tracer conservation is exact when
-eliminating tracer and surface height time filtering (see ?? for more complete discussion).
-However, in the general case with a non-zero $\alpha$,
-the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.
+(1) When \np[=.true.]{ln_rnf}{ln\_rnf} (see \autoref{sec:SBC_rnf}),
+river runoff is taken into account when computing the net freshwater flux.
-} %%end gm comment (copy of griffies book)
+(2) When \np[=.true.]{ln_isf}{ln\_isf} (see \autoref{sec:SBC_isf}), explicit or
+parameterised contributions from ice-shelf cavities are taken into account when computing
+the net freshwater flux.
-%% =================================================================================================
-\section{External forcings}
-\label{sec:D2D_forcing}
+(3) When \np[=.true.]{ln_sdw}{ln\_sdw} (see \autoref{sec:SBC_wave_sdw}), the contribution from
+divergence due to Stoke's drift is taken into account when computing the net freshwater
+flux.
-Besides the surface and bottom stresses (see the above section)
-which are introduced as boundary conditions on the vertical mixing,
-three other forcings may enter the dynamical equations by affecting the surface pressure gradient.
+(4) When \np[=.true.]{ln_isf .AND. ln_isfcpl}{ln\_isf .AND. ln\_isfcpl} (see
+\autoref{sec:SBC_isf}), the conribution from coupled ice-sheets is taken into account when
+computing the net freshwater flux.
-(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
-the atmospheric pressure is taken into account when computing the surface pressure gradient.
+(5) When \np[=.true.]{lk_asminc .AND. ln_sshinc .AND. ln_asmiau}{lk\_asminc .AND.
+ln\_sshinc .AND. ln\_asmiau} (see \autoref{sec:ASM_IAU}), the conribution from IAU
+weighted ssh increments is taken into account when computing the net freshwater flux.
-(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_TDE}),
-the tidal potential is taken into account when computing the surface pressure gradient.
-(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and SI3 is used
-(\ie\ when the sea-ice is embedded in the ocean),
-the snow-ice mass is taken into account when computing the surface pressure gradient.
\cmtgm{ missing : the lateral boundary condition !!! another external forcing
}
--
GitLab