Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
\documentclass[../main/NEMO_manual]{subfiles}
\begin{document}
\chapter{Ocean Dynamics (DYN)}
\label{chap:DYN}
\chaptertoc
\paragraph{Changes record} ~\\
{\footnotesize
\begin{tabularx}{\textwidth}{l||X|X}
Release & Author(s) & Modifications \\
\hline
{\em 4.0} & {\em ...} & {\em ...} \\
{\em 3.6} & {\em ...} & {\em ...} \\
{\em 3.4} & {\em ...} & {\em ...} \\
{\em <=3.4} & {\em ...} & {\em ...}
\end{tabularx}
}
\clearpage
Using the representation described in \autoref{chap:DOM},
several semi-discrete space forms of the dynamical equations are available depending on
the vertical coordinate used and on the conservation properties of the vorticity term.
In all the equations presented here, the masking has been omitted for simplicity.
One must be aware that all the quantities are masked fields and
that each time an average or difference operator is used, the resulting field is multiplied by a mask.
The prognostic ocean dynamics equation can be summarized as follows:
\[
\text{NXT} = \dbinom {\text{VOR} + \text{KEG} + \text {ZAD} }
{\text{COR} + \text{ADV} }
+ \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
\]
NXT stands for next, referring to the time-stepping.
The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that
are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) and
a vertical advection part (ZAD) in the vector invariant formulation,
or a Coriolis and advection part (COR+ADV) in the flux formulation.
The terms following these are the pressure gradient contributions
(HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient);
and contributions from lateral diffusion (LDF) and vertical diffusion (ZDF),
which are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules.
The vertical diffusion term includes the surface and bottom stresses.
The external forcings and parameterisations require complex inputs
(surface wind stress calculation using bulk formulae, estimation of mixing coefficients)
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 of the velocities (\emph{divhor} module) and the vertical velocity (\emph{sshwzv} 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"?}
%% =================================================================================================
%\section{Continuity equation and diagnostic variables ($\chi$, $w$)}
Andrew Coward
committed
%\label{sec:DYN_divcur_wzv}
%% =================================================================================================
\section[Continuity equation (\textit{sshwzv.F90}, $w$ )]{Continuity equation (\protect\mdl{sshwzv})}
\label{sec:DYN_wzv}
The evolution of sea surface height ($\eta$) and vertical velocity ($w$) in ocean modeling
is fundamentally governed by the continuity equation, which ensures the conservation of volume.
Through this equation, we derive expressions for $\eta$ and $w$ that are essential for maintaining
the balance between horizontal and vertical transport in the ocean's volume-conserving framework.
Both are deduced from the horizontal divergence $\chi$ and require the horizontal divergence calculation.
Because $\eta$ evolution is related to the external mode it is described in \autoref{chap:D2D}
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
while $w$ evolution is described here after.
%% =================================================================================================
\subsection[Horizontal divergence (\textit{divhor.F90}, $\chi$ )]{Horizontal divergence (\protect\mdl{divhor})}
\label{subsec:DYN_divhor}
The horizontal divergence is defined at a $T$-point.
It is given by:
\[
% \label{eq:DYN_divcur_div}
\chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
\left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right]
+\delta_j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right)
\]
\noindent Besides the velocity three other sources may be added to the horizontal divergence :
(1) When \np[=.true.]{ln_rnf}{ln\_rnf} (see \autoref{sec:SBC_rnf}),
the divergence caused by river runoff is included.
(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.
(3) 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 contribution from IAU weighted ssh increments is included.
\vskip 0.5cm
\noindent In a leapfrog time-stepping scheme, the divergence at \textit{now} time step is used
to calculate both nonlinear advection and vertical velocity.
In an RK3 time-stepping scheme, the divergence at \textit{before} time step is applied
during the first stage, while the \textit{now} time step divergence is used
for calculating nonlinear advection and vertical velocity in the following stages.
%% =================================================================================================
\subsection[Vertical velocity (\textit{sshwzv.F90})]{Vertical velocity (\protect\mdl{sshwzv})}
\label{subsec:DYN_wzv}
The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom,
taking into account the change of the thickness of the levels:
\begin{equation}
\label{eq:DYN_wzv}
\left\{
\begin{aligned}
&\left. w \right|_{k_b-1/2} \quad= 0 \qquad \text{where } k_b \text{ is the level just above the sea floor } \\
&\left. w \right|_{k+1/2} = \left. w \right|_{k-1/2} + \left. e_{3t} \right|_{k}\; \left. \chi \right|_k
- \frac{1} {2 \rdt} \left( { \left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k} }\right)
\end{aligned}
\right.
\end{equation}
In the case of a linear free surface(\key{linssh}), the time derivative in \autoref{eq:DYN_wzv} disappears.
The upper boundary condition applies at a fixed level $z=0$.
The top vertical velocity is thus equal to the divergence of the barotropic transport
(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}).
Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates,
its physical meaning is not the same:
in the second case, $w$ is the velocity normal to the $s$-surfaces.
Note also that the $k$-axis is re-orientated downwards in the \fortran\ code compared to
the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv}
(see \autoref{subsec:DOM_Num_Index_vertical}).
When \np[=.true.]{ln_zad_Aimp}{ln\_zad\_Aimp},
a proportion of the vertical advection can be treated implicitly (see \autoref{sec:DYN_zdf}) depending on the Courant number.
This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
%% =================================================================================================
\section{Coriolis and advection: vector invariant form}
\label{sec:DYN_adv_cor_vect}
\begin{listing}
\nlst{namdyn_adv}
\caption{\forcode{&namdyn_adv}}
\label{lst:namdyn_adv}
\end{listing}
Andrew Coward
committed
The vector invariant form of the momentum equation is most commonly used in
Andrew Coward
committed
coarse-resolution (1$^\textrm{o}$) applications of the \NEMO\ ocean model. For higher resolutions it
Andrew Coward
committed
requires the activation of Hollingsworth correction (\np[=1]{nn_dynkeg}{nn\_dynkeg})
following Arakawa (2001) The flux form option (see next section) has been present since
version $2$. By structuring the equations in vector invariant form, the dynamics are
expressed in terms of intrinsic geometric properties like gradients, curls, and
divergences. This ensures that the physics remain consistent and interpretable regardless
of the underlying curvilinear grid or coordinate system. It highlights key physical terms
like the kinetic energy advection and the relative vorticity.
\vskip 0.5cm
\noindent Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and
momentum advection terms are evaluated either using a leapfrog scheme or a RK3 scheme.
In the leapfrog case the velocity appearing in these expressions is centred in time (\textit{now} velocity).
In the RK3 case the velocity appearing in these expressions is forward in time (\textit{before} velocity) at stage 1,
it is is centred in time (\textit{now} velocity) at stage 2 and 3.
At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following
\autoref{chap:LBC}.
%% =================================================================================================
\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})}
\label{subsec:DYN_vor}
\begin{listing}
\nlst{namdyn_vor}
\caption{\forcode{&namdyn_vor}}
\label{lst:namdyn_vor}
\end{listing}
Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables.
Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available:
conserving potential enstrophy of horizontally non-divergent flow (ENS scheme);
conserving horizontal kinetic energy (ENE scheme);
conserving potential enstrophy for the relative vorticity term and
horizontal kinetic energy for the planetary vorticity term (MIX scheme);
or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy
(EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}).
In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of
vorticity term with analytical equations (\np[=.true.]{ln_dynvor_msk}{ln\_dynvor\_msk}).
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module.
% enstrophy conserving scheme
%% =================================================================================================
\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})}
\label{subsec:DYN_vor_ens}
In the enstrophy conserving case (ENS scheme),
the discrete formulation of the vorticity term provides a global conservation of the enstrophy
($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie\ $\chi$=$0$),
but does not conserve the total kinetic energy.
It is given by:
\begin{equation}
\label{eq:DYN_vor_ens}
\left\{
\begin{aligned}
{+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i}
& {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2} \\
{- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}
& {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j}
\end{aligned}
\right.
\end{equation}
% energy conserving scheme
%% =================================================================================================
\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})}
\label{subsec:DYN_vor_ene}
The kinetic energy conserving scheme (ENE scheme) conserves the global kinetic energy but not the global enstrophy.
It is given by:
\begin{equation}
\label{eq:DYN_vor_ene}
\left\{
\begin{aligned}
{+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
\; \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\
{- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
\; \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} }
\end{aligned}
\right.
\end{equation}
% mix energy/enstrophy conserving scheme
%% =================================================================================================
\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})}
\label{subsec:DYN_vor_mix}
For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used.
It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term,
and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term.
\[
% \label{eq:DYN_vor_mix}
\left\{ {
\begin{aligned}
{+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i}
\; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} }
\; {\overline {\left( {\frac{f}{e_{3f} }} \right)
\;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\
{-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j
\; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} }
\; {\overline {\left( {\frac{f}{e_{3f} }} \right)
\;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill
\end{aligned}
} \right.
\]
% energy and enstrophy conserving scheme
%% =================================================================================================
\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})}
\label{subsec:DYN_vor_een}
In both the ENS and ENE schemes,
it is apparent that the combination of $i$ and $j$ averages of the velocity allows for
the presence of grid point oscillation structures that will be invisible to the operator.
These structures are \textit{computational modes} that will be at least partly damped by
the momentum diffusion operator (\ie\ the subgrid-scale advection), but not by the resolved advection term.
The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field.
Such noise would result in more noise in the vertical velocity field, an undesirable feature.
This is a well-known characteristic of $C$-grid discretization where
$u$ and $v$ are located at different grid points,
a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid.
\cmtgm{ To circumvent this, Adcroft (ADD REF HERE)
Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....}
A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}.
The idea is to get rid of the double averaging by considering triad combinations of vorticity.
It is noteworthy that this solution is conceptually quite similar to the one proposed by
\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}).
The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified
for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme.
First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
\[
% \label{eq:DYN_pot_vor}
q = \frac{\zeta +f} {e_{3f} }
\]
where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}),
the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:
\begin{equation}
\label{eq:DYN_een_e3f}
e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
\end{equation}
\begin{figure}[!ht]
\centering
\includegraphics[width=0.66\textwidth]{DYN_een_triad}
\caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{
Triads used in the energy and enstrophy conserving scheme (EEN) for
$u$-component (upper panel) and $v$-component (lower panel).}
\label{fig:DYN_een_triad}
\end{figure}
A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks
(\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=0]{nn_een_e3f}{nn\_een\_e3f}).
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and
extends by continuity the value of $e_{3f}$ into the land areas.
This case introduces a sub-grid-scale topography at f-points
(with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry)
that tends to reinforce the topostrophy of the flow
(\ie\ the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}.
Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as
the following triad combinations of the neighbouring potential vorticities defined at f-points
(\autoref{fig:DYN_een_triad}):
\begin{equation}
\label{eq:DYN_Q_triads}
_i^j \mathbb{Q}^{i_p}_{j_p}
= \frac{1}{12} \ \left( q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p} \right)
\end{equation}
where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$.
Finally, the vorticity terms are represented as:
\begin{equation}
\label{eq:DYN_vor_een}
\left\{ {
\begin{aligned}
+q\,e_3 \, v &\equiv +\frac{1}{e_{1u} } \sum_{\substack{i_p,\,k_p}}
{^{i+1/2-i_p}_j} \mathbb{Q}^{i_p}_{j_p} \left( e_{1v}\,e_{3v} \;v \right)^{i+1/2-i_p}_{j+j_p} \\
- q\,e_3 \, u &\equiv -\frac{1}{e_{2v} } \sum_{\substack{i_p,\,k_p}}
{^i_{j+1/2-j_p}} \mathbb{Q}^{i_p}_{j_p} \left( e_{2u}\,e_{3u} \;u \right)^{i+i_p}_{j+1/2-j_p} \\
\end{aligned}
} \right.
\end{equation}
This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.
It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow
(\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}).
Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}.
Furthermore, used in combination with a partial steps representation of bottom topography,
it improves the interaction between current and topography,
leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}.
%% =================================================================================================
\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})}
\label{subsec:DYN_keg}
As demonstrated in \autoref{apdx:INVARIANTS},
there is a single discrete formulation of the kinetic energy gradient term that,
together with the formulation chosen for the vertical advection (see below),
conserves the total kinetic energy:
\[
% \label{eq:DYN_keg}
\left\{
\begin{aligned}
-\frac{1}{2 \; e_{1u} } & \ \delta_{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] \\
-\frac{1}{2 \; e_{2v} } & \ \delta_{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]
\end{aligned}
\right.
\]
%% =================================================================================================
\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})}
\label{subsec:DYN_zad}
The discrete formulation of the vertical advection, t
ogether with the formulation chosen for the gradient of kinetic energy (KE) term,
conserves the total kinetic energy.
Indeed, the change of KE due to the vertical advection is exactly balanced by
the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}).
\[
% \label{eq:DYN_zad}
\left\{
\begin{aligned}
-\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2} \;\delta_{k+1/2} \left[ u \right]\ }^{\,k} \\
-\frac{1} {e_{1v}\,e_{2v}\,e_{3v}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,j+1/2} \;\delta_{k+1/2} \left[ u \right]\ }^{\,k}
\end{aligned}
\right.
\]
%% =================================================================================================
\section{Coriolis and advection: flux form}
\label{sec:DYN_adv_cor_flux}
Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables.
In the flux form (as in the vector invariant form),
the Coriolis and momentum advection terms are evaluated using either a leapfrog scheme or a RK3 scheme.
In the leapfrog case the velocity appearing in these expressions is centred in time (\textit{now} velocity).
In the RK3 case the velocity appearing in these expressions is forward in time (\textit{before} velocity) at stage 1,
it is is centred in time (\textit{now} velocity) at stage 2 and 3.
At the lateral boundaries either free slip,
no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}.
%% =================================================================================================
\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
\label{subsec:DYN_cor_flux}
In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term.
This altered Coriolis parameter is thus discretised at $f$-points.
It is given by:
Andrew Coward
committed
\[
\begin{aligned}
% \label{eq:DYN_cor_metric}
f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right)
\equiv f + \frac{1}{e_{1f} e_{2f} } \left( \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right]
- \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right] \right)
Andrew Coward
committed
\end{aligned}
\]
% energy conserving scheme at T-point
%% =================================================================================================
\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_enT})]{Energy conserving scheme (\protect\np{ln_dynvor_enT}{ln\_dynvor\_enT})}
\label{subsec:DYN_vor_enT}
The kinetic energy conserving scheme at T-point (ENT scheme) conserves the global kinetic energy but not the global enstrophy.
In this case the altered Coriolis parameter is discretised at $t$-points.
ENT scheme is given by:
\begin{equation}
\label{eq:DYN_vor_enT}
\left\{
\begin{aligned}
&+\frac{1}{e_{1u}\,e_{2u}\,e_{3u}} \, \overline{ \left( f^T \right)
e_{1t}\,e_{2t}\,e_{3t} \, \overline{v}^{\,j}}^{\,i+1/2} \\
&-\frac{1}{e_{1v}\,e_{2v}\,e_{3v}} \, \overline{ \left( f^T \right)
e_{1t}\,e_{2t}\,e_{3t} \, \overline{u}^{\,i}}^{\,j+1/2}
\end{aligned}
\right.
\end{equation}
\noindent Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}), (\autoref{eq:DYN_vor_enT}) and (\autoref{eq:DYN_vor_een})
schemes can be used to
compute the product of the Coriolis parameter and the vorticity.
However, the energy-conserving schemes (\autoref{eq:DYN_vor_een} and \autoref{eq:DYN_vor_enT})
have exclusively been used to date.
\vskip 0.5cm
\noindent This term is evaluated using either a leapfrog scheme or a RK3 scheme.
In the leapfrog case it is centred in time (\textit{now} velocity).
In the RK3 case it is forward in time (\textit{before} velocity) at stage 1,
it is is centred in time (\textit{now} velocity) at stage 2 and 3.
%% =================================================================================================
\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})}
\label{subsec:DYN_adv_flux}
The discrete expression of the advection term is given by:
Andrew Coward
committed
\begin{equation}
% \label{eq:DYN_adv}
\left\{
\begin{aligned}
\frac{1}{e_{1u}\,e_{2u}\,e_{3u}}
Andrew Coward
committed
\left( \delta_{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i} \ u_t \right] \right.
+ & \delta_{j} \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2} \ u_f \right] \\
+ & \left. \delta_{k} \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2} \ u_{uw} \right] \right) \\[10pt]
Andrew Coward
committed
\left( \delta_{i} \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f \right] \right.
+ & \delta_{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i} \ v_t \right] \\
+ & \left. \delta_{k} \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw} \right] \right) \\
Andrew Coward
committed
\end{equation}
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
Two advection schemes are available:
a $2^{nd}$ order centered finite difference scheme, CEN2,
or a $3^{rd}$ order upstream biased scheme, UP3.
The latter is described in \citet{shchepetkin.mcwilliams_OM05}.
The schemes are selected using the namelist logicals \np{ln_dynadv_cen2}{ln\_dynadv\_cen2} and \np{ln_dynadv_up3}{ln\_dynadv\_up3}.
In flux form, the schemes differ by the choice of a space and time interpolation to define the value of
$u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-,
and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$.
% 2nd order centred scheme
%% =================================================================================================
\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})}
\label{subsec:DYN_adv_cen2}
In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points:
\begin{equation}
\label{eq:DYN_adv_cen2}
\left\{
\begin{aligned}
u_T^{cen2} &=\overline u^{i } \quad & u_F^{cen2} &=\overline u^{j+1/2} \quad & u_{uw}^{cen2} &=\overline u^{k+1/2} \\
v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j \quad & v_{vw}^{cen2} &=\overline v ^{k+1/2} \\
\end{aligned}
\right.
\end{equation}
The scheme is non diffusive (\ie\ conserves the kinetic energy) but dispersive (\ie\ it may create false extrema).
It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to
produce a sensible solution.
% UP3 scheme
%% =================================================================================================
\subsubsection[UP3: Upstream Biased Scheme (\forcode{ln_dynadv_up3})]{UP3: Upstream Biased Scheme (\protect\np{ln_dynadv_up3}{ln\_dynadv\_up3})}
\label{subsec:DYN_adv_up3}
The UP3 advection scheme is an upstream biased third order scheme based on
an upstream-biased parabolic interpolation.
For example, the evaluation of $u_T^{up3} $ is done as follows:
\begin{equation}
\label{eq:DYN_adv_up3}
u_T^{up3} =\overline u ^i-\;\frac{1}{6}
\begin{cases}
u"_{i-1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i \geqslant 0$ } \\
u"_{i+1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i < 0$ }
\end{cases}
\end{equation}
where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$.
This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error
\citep{shchepetkin.mcwilliams_OM05}.
The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}.
It is a relatively good compromise between accuracy and smoothness.
It is not a \emph{positive} scheme, meaning that false extrema are permitted.
But the amplitudes of the false extrema are significantly reduced over those in the centred second order method.
As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum
(\ie\ \np[=.true.]{ln_dynldf_OFF}{ln\_dynldf\_OFF}),
and it is recommended to do so.
The UP3 scheme is used in all directions.
UP3 is diffusive and is associated with vertical mixing of momentum. \cmtgm{ gm pursue the
sentence:Since vertical mixing of momentum is a source term of the TKE equation... }
In a leapfrog environment, for stability reasons, the first term in (\autoref{eq:DYN_adv_up3}),
which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time),
while the second term, which is the diffusion part of the scheme,
is evaluated using the \textit{before} velocity (forward in time).
In an RK3 environment, the first term in (\autoref{eq:DYN_adv_up3}),
which corresponds to a second order centred scheme, is evaluated using the \textit{before} velocity at stage 1
and using the \textit{before} velocity (centred in time) at stage 2 and 3,
while the second term, which is the diffusion part of the scheme,
is evaluated using the \textit{before} velocity (forward in time).
This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
Note that the UP3 and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
one coefficient.
Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_up3}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_up3} module and obtain a QUICK scheme.
%% =================================================================================================
\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
\label{sec:DYN_hpg}
\begin{listing}
\nlst{namdyn_hpg}
\caption{\forcode{&namdyn_hpg}}
\label{lst:namdyn_hpg}
\end{listing}
%Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables.
%The key distinction between the different algorithms used for
%the hydrostatic pressure gradient is the vertical coordinate used,
%since HPG is a \emph{horizontal} pressure gradient, \ie\ computed along geopotential surfaces.
%As a result, any tilt of the surface of the computational levels will require a specific treatment to
%compute the hydrostatic pressure gradient.
%In v5.0 partial cells are considered as penalized cells. A penalized cell has a solid fraction and
%a liquid fraction, both are distributed homogeneouly across the cell. Therefore adjacent cells are not
%expected to live at distinct depth because of their content differences.
NEMO offers a selection of different algorithms to compute the hydrostatic pressure gradient (HPG) term in the momentum equation.
Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables.
Since HPGs are computed along geopotential surfaces, a key distinction between the various algorithms is the type of vertical
coordinate they target. In particular, NEMO offers a number of options to compute HPGs with generalised $s$-coordinates that
may be not aligned with geopotentials.
The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme,
\ie\ the density appearing in its expression is centred in time (\emph{now} $\rho$),
or a RK3 scheme \ie\ the density appearing in its expression is forward in time (\emph{before} $\rho$),
it is centred in time (\emph{now} $\rho$) at stage 2 and 3.
At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied.
%% =================================================================================================
\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})}
\label{subsec:DYN_hpg_zco}
When using standard geopotential coordinates (\np[=.true.]{ln_zco}{ln\_zco}), the hydrostatic pressure can be directly obtained
by vertically integrating the hydrostatic equation from the surface to the bottom. However, pressure is large at great depths while
its horizontal gradient is several orders of magnitude smaller. This may lead to large truncation errors in the pressure gradient terms.
Thus, the two horizontal components of the hydrostatic pressure gradient are computed directly as follows:
for $k=km$ (surface layer, $jk=1$ in the code)
\begin{equation}
\label{eq:DYN_hpg_zco_surf}
\left\{
\begin{aligned}
\left. \delta_{i+1/2} \left[ p^h \right] \right|_{k=km}
&= \frac{1}{2} g \ \left. \delta_{i+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\
\left. \delta_{j+1/2} \left[ p^h \right] \right|_{k=km}
&= \frac{1}{2} g \ \left. \delta_{j+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\
\end{aligned}
\right.
\end{equation}
for $1<k<km$ (interior layer)
\begin{equation}
\label{eq:DYN_hpg_zco}
\left\{
\begin{aligned}
\left. \delta_{i+1/2} \left[ p^h \right] \right|_{k}
&= \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k-1}
+ g\; \left. \delta_{i+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\
\left. \delta_{j+1/2} \left[ p^h \right] \right|_{k}
&= \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k-1}
+ g\; \left. \delta_{j+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\
\end{aligned}
\right.
\end{equation}
Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
the vertical derivative of the scale factor at the surface level ($z=0$).
%% =================================================================================================
\subsection{Generalised $S$-coordinates}
\label{subsec:DYN_hpg_sco}
Pressure gradient formulations with a generalised $s(x,y,z,t)$ coordinate have been the subject of a vast number of papers
(\eg, \citet{song_MWR98, shchepetkin.mcwilliams_jgro03}). A number of different pressure gradient options are available in NEMO:
\begin{itemize}
\item
\textbf{Traditional coding} (\np[=.true.]{ln_hpg_sco}{ln\_hpg\_sco}, e.g. \citet{madec.delecluse.ea_JPO96}):
%\end{itemize}
\begin{equation}
\label{eq:DYN_hpg_sco}
\left\{
\begin{aligned}
- \frac{1} {\rho_o \, e_{1u}} \; \delta_{i+1/2} \left[ p^h \right]
+ \frac{g\; \overline {\rho}^{i+1/2}} {\rho_o \, e_{1u}} \; \delta_{i+1/2} \left[ z_t \right], \\
- \frac{1} {\rho_o \, e_{2v}} \; \delta_{j+1/2} \left[ p^h \right]
+ \frac{g\; \overline {\rho}^{j+1/2}} {\rho_o \, e_{2v}} \; \delta_{j+1/2} \left[ z_t \right], \\
\end{aligned}
\right.
\end{equation}
where the first term is the pressure gradient along coordinates (computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco})
and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $W$-point ($e_{3w}$).
Note that this scheme is not recommended when using steeply inclined computational levels (e.g., terrain-following or hybrid generalised
vertical coordinates, i.e., \np[=.true.]{ln_sco}{ln\_sco}) - see e.g. \citet{shchepetkin.mcwilliams_jgro03}.
However, it should be the standard choice when using $z$-coordinates (\np[=.true.]{ln_zco}{ln\_zps} or \np[=.true.]{ln_zps}{ln\_zps})
with the non-linear free surface (\np[=.false.]{ln_linssh}{ln\_linssh} and \key{qco}), since in this case model levels will follow the barotropic motion of the ocean \citep{levier.treguier.ea_trpt07}.
\textbf{Traditional coding with adaptation for ice shelf cavities} (\np[=.true.]{ln_hpg_isf}{ln\_hpg\_isf}):
In the presence of ice shelves, the traditional coding has been adapted to accommodate the load provided by the ice shelves.
This scheme must be used when ice shelf cavities are activated (\np[=.true.]{ln_isfcav}{ln\_isfcav} and the inclusion of \key{isf}.
All the details on the modification are provided in \autoref{subsec:ISF_hpg}.
\textbf{Pressure Jacobian scheme} (\np[=.true.]{ln_hpg_prj}{ln\_hpg\_prj}):
this scheme uses a constrained cubic spline to reconstruct the vertical density profile within a water column.
This method maintains the monotonicity between the density nodes. The pressure is calculated by analytical integration of the density
profile. For the force in the $i$-direction, it calculates the difference of the pressures on the $i+\tfrac{1}{2}$ and $i-\tfrac{1}{2}$
faces of the cell using pressures calculated at the same height. In grid cells just above the bathymetry, this height is higher than the
cells’ centre. This scheme works well for moderately steep computational levels but produces large velocities in the SEAMOUNT test case
when model levels are steeply inclined.
\textbf{Density Jacobian with cubic polynomial scheme} (\np[=.true.]{ln_hpg_djc}{ln\_hpg\_djc}, \cite{shchepetkin.mcwilliams_jgro03}):
the ROMS-like, density Jacobian with cubic polynomial method has been debugged and from vn4.2 is available as an option.
This scheme is based on section 5 of \cite{shchepetkin.mcwilliams_jgro03} For the force in the $i$-direction, it uses constrained cubic
splines to re-construct the density along lines of constant $s$ and constant $i$ in the $(i,s)$ plane. It calculates a line integral of
$\rho$ and then integrates vertically to obtain the horizontal pressure gradient. The constrained cubic splines require boundary conditions
to be specified at the upper and lower boundaries and at points where model levels encrop the model bathymetry (i.e., with geopotential
or hybrid vertical coordinates). The user can choose between von Neumann and linear extrapolation boundary conditions via the
\texttt{ln\_hpg\_djc\_vnh} and \texttt{ln\_hpg\_djc\_vnv} namelist switches, respectively. This scheme can be used with any type of
generalised $s$-coordinates - i.e., $z$ or $z^*$, terrain-following or hybrids of these two (e.g., via the vanishing quasi-sigma or
multi-envelope methods, see e.g. \cite{shapiro.luneva.ea_OS13,bruciaferri.shapiro.wobus_OD18,wise.harle.ea_OM21}) - but at the moment
can not be used with ice shelf cavities.
\end{itemize}
Starting from version 4.2, the density field used by dyn\_hpg is the density anomaly field rhd rather than $1+\mathrm{rhd}$.
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
The calculation of the source term for the free surface has been adjusted to take this into account.
The true in situ density $\rho= \rho_0 (1 + r_0(z) + rhd )$ where $r_0(z)$ accounts for the variation of density
with depth for water with a potential temperature of $4^{\circ}$C and salinity of $35.16504$g/kg
(see (13) and (14) of \citet{roquet.madec.ea_OM15}).
%% =================================================================================================
%% \subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})}
%% \label{subsec:DYN_hpg_imp}
%%
%% The default time differencing scheme used for the horizontal pressure gradient is a leapfrog scheme and
%% therefore the density used in all discrete expressions given above is the \textit{now} density,
%% computed from the \textit{now} temperature and salinity.
%% In some specific cases
%% (usually high resolution simulations over an ocean domain which includes weakly stratified regions)
%% the physical phenomenon that controls the time-step is internal gravity waves (IGWs).
%% A semi-implicit scheme for doubling the stability limit associated with IGWs can be used
%% \citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}.
%% It involves the evaluation of the hydrostatic pressure gradient as
%% an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$
%% (\ie\ \textit{before}, \textit{now} and \textit{after} time-steps),
%% rather than at the central time level $t$ only, as in the standard leapfrog scheme.
%%
%% $\bullet$ leapfrog scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
%%
%% \begin{equation}
%% \label{eq:DYN_hpg_lf}
%% \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
%% -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right]
%% \end{equation}
%%
%% $\bullet$ semi-implicit scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
%% \begin{equation}
%% \label{eq:DYN_hpg_imp}
%% \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
%% -\frac{1}{4\,\rho_o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt} \right]
%% \end{equation}
%%
%% The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without
%% significant additional computation since the density can be updated to time level $t+\rdt$ before
%% computing the horizontal hydrostatic pressure gradient.
%% It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using
%% \autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}.
%% Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to
%% eliminate high frequency IGWs.
%% Obviously, when using \autoref{eq:DYN_hpg_imp},
%% the doubling of the time-step is achievable only if no other factors control the time-step,
%% such as the stability limits associated with advection or diffusion.
%%
%% In practice, the semi-implicit scheme is used when \np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}.
%% In this case, we choose to apply the time filter to temperature and salinity used in the equation of state,
%% instead of applying it to the hydrostatic pressure or to the density,
%% so that no additional storage array has to be defined.
%% The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows:
%% \[
%% % \label{eq:DYN_rho_flt}
%% \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
%% \quad \text{with} \quad
%% \widetilde{X} = 1 / 4 \left( X^{t+\rdt} +2 \,X^t + X^{t-\rdt} \right)
%% \]
%%
%% Note that in the semi-implicit case, it is necessary to save the filtered density,
%% an extra three-dimensional field, in the restart file to restart the model with exact reproducibility.
%% This option is controlled by \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter.
%%
Andrew Coward
committed
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
\section[Surface pressure gradient (\textit{dynspg.F90}/\textit{stp2d.F90})]
{Surface pressure gradient (\protect\mdl{dynspg}/\protect\mdl{stp2d})}
\label{sec:DYN_spg}
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
(\autoref{subsec:MB_free_surface}) the vertical scale factors $e_{3}$ are fixed in time,
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} (MLF time-stepping only)). With explicit
time-stepping, the surface pressure gradient is evaluated using the leap-frog scheme
(\ie centred in time) and is thus simply given by:
\begin{equation}
\label{eq:DYN_spg_exp}
\left\{
\begin{aligned}
- \frac{1}{e_{1u}\,\rho_o} \; \delta_{i+1/2} \left[ \,\rho \,\eta\, \right] \\
- \frac{1}{e_{2v}\,\rho_o} \; \delta_{j+1/2} \left[ \,\rho \,\eta\, \right]
\end{aligned}
\right.
\end{equation}
where values are evaluated at the now timestep (\textit{dynspg\_exp.F90}).
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,
a quasi-linear form of 2d barotropic equations is substepped with a small time
increment. Details of this were provided in \autoref{chap:D2D}. Options are defined
through the \nam{dyn_spg}{dyn\_spg} namelist variables.
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
%% =================================================================================================
\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
\label{sec:DYN_ldf}
\begin{listing}
\nlst{namdyn_ldf}
\caption{\forcode{&namdyn_ldf}}
\label{lst:namdyn_ldf}
\end{listing}
Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables.
The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators.
The coefficients may be constant or spatially variable;
the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}).
The lateral diffusion of momentum is evaluated using a forward scheme,
\ie\ the velocity appearing in its expression is the \textit{before} velocity in time,
except for the pure vertical component that appears when a tensor of rotation is used.
This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}).
At the lateral boundaries either free slip,
no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}).
\cmtgm{
Hyperviscous operators are frequently used in the simulation of turbulent flows to
control the dissipation of unresolved small scale features.
Their primary role is to provide strong dissipation at the smallest scale supported by
the grid while minimizing the impact on the larger scale features.
Hyperviscous operators are thus designed to be more scale selective than the traditional,
physically motivated Laplace operator.
In finite difference methods,
the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since
its damping time (\ie\ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$
(so that short waves damped more rapidelly than long ones),
whereas the Laplace operator damping time scales only like $\lambda^{-2}$.
}
%% =================================================================================================
\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}
\label{subsec:DYN_ldf_lap}
For lateral iso-level diffusion (\np[=0]{nn_dynldf_typ}{nn\_dynldf\_typ}), the discrete operator is:
\begin{equation}
\label{eq:DYN_ldf_lap}
\left\{
\begin{aligned}
D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
\;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[
{A_f^{lm} \;e_{3f} \zeta } \right] \\ \\
D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
\;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[
{A_f^{lm} \;e_{3f} \zeta } \right]
\end{aligned}
\right.
\end{equation}
As explained in \autoref{sec:DIFFOPERS_3}, this formulation (as the gradient of a divergence and curl of the vorticity) preserves
symmetry and ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.
In v5.0 a symetrical lateral iso-level operator (\np[=1]{nn_dynldf_typ}{nn\_dynldf\_typ}) has been introduced :
\begin{equation}
\label{eq:DYN_ldf_lap_sym}
\left\{
\begin{aligned}
D_u^{l{\mathrm {\mathbf U}}} &= \frac{1}{e_{1u}\,e_{2u}\,e_{3u} } \left(
\frac{1}{e_{2u}} \delta_{i+1/2} \left[ e_{2t}\,e_{2t}\,e_{3t}\, A_T^{lm} \,\epsilon_T \right]
- \frac{1}{e_{1u}} \delta_{j+1/2} \left[ e_{1f}\,e_{1f}\,e_{3f}\,A_F^{lm} \epsilon_F \right] \right) \\ \\
D_v^{l{\mathrm {\mathbf U}}} &= \frac{1}{e_{1v}\,e_{2v}\,e_{3v} } \left(
\frac{1}{e_{2v}} \delta_{j+1/2} \left[ e_{2f}\,e_{2f}\,e_{3f}\, A_F^{lm} \,\epsilon_F \right]
- \frac{1}{e_{1v}} \delta_{i+1/2} \left[ e_{1t}\,e_{1t}\,e_{3t}\,A_T^{lm} \epsilon_T \right] \right)
Where $\epsilon_F$ and $\epsilon_T$ are respectively the shearing stress component (F-point) and the
tension stress component (T-point) defined as :
\begin{equation}
\label{eq:DYN_ldf_lap_sheten}
\left\{
\begin{aligned}
\epsilon_F &= \frac{e_{1f}}{e_{2f}}\; \delta_{j+1/2} \left[ \frac{u}{e_{1u}} \right] + \frac{e_{2f}}{e_{1f}}\;
\delta_{i+1/2} \left[ \frac{v}{e_{2v}} \right] \\ \\
\epsilon_T &= \frac{e_{2t}}{e_{1t}}\; \delta_{i} \left[ \frac{u}{e_{2u}} \right] - \frac{e_{1t}}{e_{2t}}\;
\delta_j \left[ \frac{v}{e_{1v}} \right]
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
%% =================================================================================================
\subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})}
\label{subsec:DYN_ldf_iso}
A rotation of the lateral momentum diffusion operator is needed in several cases:
for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and
for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or
geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate.
In the partial step case, coordinates are horizontal except at the deepest level and
no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on
each momentum component.
It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry.
The resulting discrete representation is:
\begin{equation}
\label{eq:DYN_ldf_iso}
\begin{split}
D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
& \left\{\quad {\delta_{i+1/2} \left[ {A_T^{lm} \left(
{\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u]
-e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}}
\right)} \right]} \right. \\
& \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f}
}\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f}
\,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}}
\right)} \right] \\
&\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline
{\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} }
\right.} \right. \\
& \ \qquad \qquad \qquad \quad\
- e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\
& \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
+\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2}
\right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\
D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\
& \left\{\quad {\delta_{i+1/2} \left[ {A_f^{lm} \left(
{\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v]
-e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}}
\right)} \right]} \right. \\
& \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t}
}\,\delta_{j} [v] - e_{1t}\, r_{2t}
\,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}}
\right)} \right] \\
& \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline
{\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\
& \ \qquad \qquad \qquad \quad\
- e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\
& \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
+\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2}
\right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\}
\end{split}
\end{equation}
where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and
the surface of computation ($z$- or $s$-surfaces).
The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).
%% =================================================================================================
\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})}
\label{subsec:DYN_ldf_bilap}
The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice.
It requires an additional assumption on boundary conditions:
the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,
while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}).
\cmtgm{add a remark on the the change in the position of the coefficient}
%% =================================================================================================
\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}
\label{sec:DYN_zdf}
Options are defined through the \nam{zdf}{zdf} namelist variables.
The large vertical diffusion coefficient found in the surface mixed layer together with high vertical resolution implies that in the case of explicit time stepping there would be too restrictive a constraint on the time step.
In v5.0 only a backward (or implicit) time differencing scheme can be used for the vertical diffusion term.
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
(see \autoref{chap:TD}).
The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.
The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form:
\[
% \label{eq:DYN_zdf}
\left\{
\begin{aligned}
D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
\ \delta_{k+1/2} [\,u\,] \right] \\
\\
D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
\ \delta_{k+1/2} [\,v\,] \right]
\end{aligned}
\right.
\]
where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients.
The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}).
The surface boundary condition on momentum is the stress exerted by the wind.
At the surface, the momentum fluxes are prescribed as the boundary condition on
the vertical turbulent momentum fluxes,
\begin{equation}
\label{eq:DYN_zdf_sbc}
\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
= \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }
\end{equation}
where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in
the (\textbf{i},\textbf{j}) coordinate system.
The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in
the vertical over the mixed layer depth.
If the vertical mixing coefficient is small (when no mixed layer scheme is used)
the surface stress enters only the top model level, as a body force.
The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}).
The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
(see \autoref{sec:ZDF_drg})
When activated (\np[=.true.]{ln_zad_Aimp}{ln\_zad\_Aimp}) vertical advection of momentum is
done partly implicitly in areas where there is potential to breach the vertical CFL condition
(see \autoref{subsec:ZDF_aimp}). The following expressions are included in the implicit solvers in
\textit{dynzdf.F90} and \textit{trazdf.F90}: