Last change on this file since 14467 was 14467, checked in by nicolasmartin, 10 months ago

Fix documentation generation (manual and guide)

File size: 97.4 KB
Line
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5\chapter{Vertical Ocean Physics (ZDF)}
6\label{chap:ZDF}
7
8\thispagestyle{plain}
9
10\chaptertoc
11
12\paragraph{Changes record} ~\\
13
14{\footnotesize
15  \begin{tabularx}{\textwidth}{l||X|X}
16    Release & Author(s) & Modifications \\
17    \hline
18    {\em   4.0} & {\em ...} & {\em ...} \\
19    {\em   3.6} & {\em ...} & {\em ...} \\
20    {\em   3.4} & {\em ...} & {\em ...} \\
21    {\em <=3.4} & {\em ...} & {\em ...}
22  \end{tabularx}
23}
24
25\clearpage
26
27\cmtgm{ Add here a small introduction to ZDF and naming of the different physics
28(similar to what have been written for TRA and DYN).}
29
30%% =================================================================================================
31\section{Vertical mixing}
32\label{sec:ZDF}
33
34The discrete form of the ocean subgrid scale physics has been presented in
35\autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}.
36At the surface and bottom boundaries, the turbulent fluxes of momentum, heat and salt have to be defined.
37At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}),
38while at the bottom they are set to zero for heat and salt,
39unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie\ \np{ln_trabbc}{ln\_trabbc} defined,
40see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum
41(see \autoref{sec:ZDF_drg}).
42
43In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and
44diffusivity coefficients, $A_u^{vm}$ , $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- points,
45respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}).
46These coefficients can be assumed to be either constant, or a function of the local Richardson number,
47or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation).
48The computation of these coefficients is initialized in the \mdl{zdfphy} module and performed in
49the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules.
50The trends due to the vertical momentum and tracer diffusion, including the surface forcing,
51are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.
52%These trends can be computed using either a forward time stepping scheme
53%(namelist parameter \np[=.true.]{ln_zdfexp}{ln\_zdfexp}) or a backward time stepping scheme
54%(\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) depending on the magnitude of the mixing coefficients,
55%and thus of the formulation used (see \autoref{chap:TD}).
56
57\begin{listing}
58  \nlst{namzdf}
59  \caption{\forcode{&namzdf}}
60  \label{lst:namzdf}
61\end{listing}
62
63%% =================================================================================================
64\subsection[Constant (\forcode{ln_zdfcst})]{Constant (\protect\np{ln_zdfcst}{ln\_zdfcst})}
65\label{subsec:ZDF_cst}
66
67Options are defined through the \nam{zdf}{zdf} namelist variables.
68When \np{ln_zdfcst}{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to
69constant values over the whole ocean.
70This is the crudest way to define the vertical ocean physics.
71It is recommended to use this option only in process studies, not in basin scale simulations.
72Typical values used in this case are:
73\begin{align*}
74  A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}   \\
75  A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
76\end{align*}
77
78These values are set through the \np{rn_avm0}{rn\_avm0} and \np{rn_avt0}{rn\_avt0} namelist parameters.
79In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity,
80that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and
81$\sim10^{-9}~m^2.s^{-1}$ for salinity.
82
83%% =================================================================================================
84\subsection[Richardson number dependent (\forcode{ln_zdfric})]{Richardson number dependent (\protect\np{ln_zdfric}{ln\_zdfric})}
85\label{subsec:ZDF_ric}
86
87\begin{listing}
88  \nlst{namzdf_ric}
89  \caption{\forcode{&namzdf_ric}}
90  \label{lst:namzdf_ric}
91\end{listing}
92
93When \np[=.true.]{ln_zdfric}{ln\_zdfric}, a local Richardson number dependent formulation for the vertical momentum and
94tracer eddy coefficients is set through the \nam{zdf_ric}{zdf\_ric} namelist variables.
95The vertical mixing coefficients are diagnosed from the large scale variables computed by the model.
96\textit{In situ} measurements have been used to link vertical turbulent activity to large scale ocean structures.
97The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to
98a dependency between the vertical eddy coefficients and the local Richardson number
99(\ie\ the ratio of stratification to vertical shear).
100Following \citet{pacanowski.philander_JPO81}, the following formulation has been implemented:
101102 % \label{eq:ZDF_ric} 103 \left\{ 104 \begin{aligned} 105 A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT} \\ 106 A^{vm} &= \frac{A^{vT} }{\left( 1+ a \;Ri \right) } + A_b^{vm} 107 \end{aligned} 108 \right. 109
110where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number,
111$N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
112$A_b^{vT}$ and $A_b^{vm}$ are the constant background values set as in the constant case
113(see \autoref{subsec:ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ is the maximum value that
114can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$.
115The last three values can be modified by setting the \np{rn_avmri}{rn\_avmri}, \np{rn_alp}{rn\_alp} and
116\np{nn_ric}{nn\_ric} namelist parameters, respectively.
117
118A simple mixing-layer model to transfer and dissipate the atmospheric forcings
119(wind-stress and buoyancy fluxes) can be activated setting the \np[=.true.]{ln_mldw}{ln\_mldw} in the namelist.
120
121In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ is evaluated and
122the vertical eddy coefficients prescribed within this layer.
123
124This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
125$126 h_{e} = Ek \frac {u^{*}} {f_{0}} 127$
128where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis parameter.
129
130In this similarity height relationship, the turbulent friction velocity:
131$132 u^{*} = \sqrt \frac {|\tau|} {\rho_o} 133$
134is computed from the wind stress vector $|\tau|$ and the reference density $\rho_o$.
135The final $h_{e}$ is further constrained by the adjustable bounds \np{rn_mldmin}{rn\_mldmin} and \np{rn_mldmax}{rn\_mldmax}.
136Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to
137the empirical values \np{rn_wtmix}{rn\_wtmix} and \np{rn_wvmix}{rn\_wvmix} \citep{lermusiaux_JMS01}.
138
139%% =================================================================================================
140\subsection[TKE turbulent closure scheme (\forcode{ln_zdftke})]{TKE turbulent closure scheme (\protect\np{ln_zdftke}{ln\_zdftke})}
141\label{subsec:ZDF_tke}
142
143\begin{listing}
144  \nlst{namzdf_tke}
145  \caption{\forcode{&namzdf_tke}}
146  \label{lst:namzdf_tke}
147\end{listing}
148
149The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on
150a prognostic equation for $\bar{e}$, the turbulent kinetic energy,
151and a closure assumption for the turbulent length scales.
152This turbulent closure model has been developed by \citet{bougeault.lacarrere_MWR89} in the atmospheric case,
153adapted by \citet{gaspar.gregoris.ea_JGR90} for the oceanic case, and embedded in OPA, the ancestor of \NEMO,
154by \citet{blanke.delecluse_JPO93} for equatorial Atlantic simulations.
155Since then, significant modifications have been introduced by \citet{madec.delecluse.ea_NPM98} in both the implementation and
156the formulation of the mixing length scale.
157The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear,
158its destruction through stratification, its vertical diffusion, and its dissipation of \citet{kolmogorov_IANS42} type:
159\begin{equation}
160  \label{eq:ZDF_tke_e}
161  \frac{\partial \bar{e}}{\partial t} =
162  \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
163      +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
164  -K_\rho\,N^2
165  +\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
166      \;\frac{\partial \bar{e}}{\partial k}} \right]
167  - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
168\end{equation}
169$170 % \label{eq:ZDF_tke_kz} 171 \begin{split} 172 K_m &= C_k\ l_k\ \sqrt {\bar{e}\; } \\ 173 K_\rho &= A^{vm} / P_{rt} 174 \end{split} 175$
176where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
177$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales,
178$P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients.
179The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with
180vertical mixing at any depth \citep{gaspar.gregoris.ea_JGR90}.
181They are set through namelist parameters \np{nn_ediff}{nn\_ediff} and \np{nn_ediss}{nn\_ediss}.
182$P_{rt}$ can be set to unity or, following \citet{blanke.delecluse_JPO93}, be a function of the local Richardson number, $R_i$:
183\begin{align*}
184  % \label{eq:ZDF_prt}
185  P_{rt} =
186  \begin{cases}
187    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}   \\
188    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}   \\
189    \ \ 10 &      \text{if $\ 2 \leq R_i$}
190  \end{cases}
191\end{align*}
192The choice of $P_{rt}$ is controlled by the \np{nn_pdl}{nn\_pdl} namelist variable.
193
194At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as
195$\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter.
196The default value of $e_{bb}$ is 3.75. \citep{gaspar.gregoris.ea_JGR90}), however a much larger value can be used when
197taking into account the surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}).
198The bottom value of TKE is assumed to be equal to the value of the level just above.
199The time integration of the $\bar{e}$ equation may formally lead to negative values because
200the numerical scheme does not ensure its positivity.
201To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn_emin}{rn\_emin} namelist parameter).
202Following \citet{gaspar.gregoris.ea_JGR90}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.
203This allows the subsequent formulations to match that of \citet{gargett_JMR84} for the diffusion in
204the thermocline and deep ocean :  $K_\rho = 10^{-3} / N$.
205In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with
206too weak vertical diffusion.
207They must be specified at least larger than the molecular values, and are set through \np{rn_avm0}{rn\_avm0} and
208\np{rn_avt0}{rn\_avt0} (\nam{zdf}{zdf} namelist, see \autoref{subsec:ZDF_cst}).
209
210%% =================================================================================================
211\subsubsection{Turbulent length scale}
212
213For computational efficiency, the original formulation of the turbulent length scales proposed by
214\citet{gaspar.gregoris.ea_JGR90} has been simplified.
215Four formulations are proposed, the choice of which is controlled by the \np{nn_mxl}{nn\_mxl} namelist parameter.
216The first two are based on the following first order approximation \citep{blanke.delecluse_JPO93}:
217\begin{equation}
218  \label{eq:ZDF_tke_mxl0_1}
219  l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
220\end{equation}
221which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency.
222The resulting length scale is bounded by the distance to the surface or to the bottom
223(\np[=0]{nn_mxl}{nn\_mxl}) or by the local vertical scale factor (\np[=1]{nn_mxl}{nn\_mxl}).
224\citet{blanke.delecluse_JPO93} notice that this simplification has two major drawbacks:
225it makes no sense for locally unstable stratification and the computation no longer uses all
226the information contained in the vertical density profile.
227To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np[=2, 3]{nn_mxl}{nn\_mxl} cases,
228which add an extra assumption concerning the vertical gradient of the computed length scale.
229So, the length scales are first evaluated as in \autoref{eq:ZDF_tke_mxl0_1} and then bounded such that:
230\begin{equation}
231  \label{eq:ZDF_tke_mxl_constraint}
232  \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
233  \qquad \text{with }\  l =  l_k = l_\epsilon
234\end{equation}
235\autoref{eq:ZDF_tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than
236the variations of depth.
237It provides a better approximation of the \citet{gaspar.gregoris.ea_JGR90} formulation while being much less
238time consuming.
239In particular, it allows the length scale to be limited not only by the distance to the surface or
240to the ocean bottom but also by the distance to a strongly stratified portion of the water column such as
241the thermocline (\autoref{fig:ZDF_mixing_length}).
242In order to impose the \autoref{eq:ZDF_tke_mxl_constraint} constraint, we introduce two additional length scales:
243$l_{up}$ and $l_{dwn}$, the upward and downward length scales, and
244evaluate the dissipation and mixing length scales as
245(and note that here we use numerical indexing):
246\begin{figure}[!t]
247  \centering
248  \includegraphics[width=0.66\textwidth]{ZDF_mixing_length}
249  \caption[Mixing length computation]{Illustration of the mixing length computation}
250  \label{fig:ZDF_mixing_length}
251\end{figure}
252253 % \label{eq:ZDF_tke_mxl2} 254 \begin{aligned} 255 l_{up\ \ }^{(k)} &= \min \left( l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \; \right) 256 \quad &\text{ from k=1 to jpk }\ \\ 257 l_{dwn}^{(k)} &= \min \left( l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)\right) 258 \quad &\text{ from k=jpk to 1 }\ \\ 259 \end{aligned} 260
261where $l^{(k)}$ is computed using \autoref{eq:ZDF_tke_mxl0_1}, \ie\ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
262
263In the \np[=2]{nn_mxl}{nn\_mxl} case, the dissipation and mixing length scales take the same value:
264$l_k= l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right)$, while in the \np[=3]{nn_mxl}{nn\_mxl} case,
265the dissipation and mixing turbulent length scales are give as in \citet{gaspar.gregoris.ea_JGR90}:
266267 % \label{eq:ZDF_tke_mxl_gaspar} 268 \begin{aligned} 269 & l_k = \sqrt{\ l_{up} \ \ l_{dwn}\ } \\ 270 & l_\epsilon = \min \left(\ l_{up} \;,\; l_{dwn}\ \right) 271 \end{aligned} 272
273
274At the ocean surface, a non zero length scale is set through the  \np{rn_mxl0}{rn\_mxl0} namelist parameter.
275Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and
276$z_o$ the roughness parameter of the surface.
277Assuming $z_o=0.1$~m \citep{craig.banner_JPO94} leads to a 0.04~m, the default value of \np{rn_mxl0}{rn\_mxl0}.
278In the ocean interior a minimum length scale is set to recover the molecular viscosity when
279$\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
280
281%% =================================================================================================
282\subsubsection{Surface wave breaking parameterization}
283
284Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to
285include the effect of surface wave breaking energetics.
286This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow.
287The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and
288air-sea drag coefficient.
289The latter concerns the bulk formulae and is not discussed here.
290
291Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is :
292\begin{equation}
293  \label{eq:ZDF_Esbc}
294  \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
295\end{equation}
296where $\alpha_{CB}$ is the \citet{craig.banner_JPO94} constant of proportionality which depends on the ''wave age'',
297ranging from 57 for mature waves to 146 for younger waves \citep{mellor.blumberg_JPO04}.
298The boundary condition on the turbulent length scale follows the Charnock's relation:
299\begin{equation}
300  \label{eq:ZDF_Lsbc}
301  l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
302\end{equation}
303where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
304\citet{mellor.blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by
305\citet{stacey_JPO99} citing observation evidence, and
306$\alpha_{CB} = 100$ the Craig and Banner's value.
307As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$,
308with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter, setting \np[=67.83]{rn_ebb}{rn\_ebb} corresponds
309to $\alpha_{CB} = 100$.
310Further setting  \np[=.true.]{ln_mxl0}{ln\_mxl0},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,
311with $\beta$ hard coded to the Stacey's value.
312Note that a minimal threshold of \np{rn_emin0}{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the
313surface $\bar{e}$ value.
314
315%% =================================================================================================
316\subsubsection{Langmuir cells}
317
318Langmuir circulations (LC) can be described as ordered large-scale vertical motions in
319the surface layer of the oceans.
320Although LC have nothing to do with convection, the circulation pattern is rather similar to
321so-called convective rolls in the atmospheric boundary layer.
322The detailed physics behind LC is described in, for example, \citet{craik.leibovich_JFM76}.
323The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and
324wind drift currents.
325
326Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by
327\citep{axell_JGR02} for a $k-\epsilon$ turbulent closure.
328The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in
329an extra source term of TKE, $P_{LC}$.
330The presence of $P_{LC}$ in \autoref{eq:ZDF_tke_e}, the TKE equation, is controlled by setting \np{ln_lc}{ln\_lc} to
331\forcode{.true.} in the \nam{zdf_tke}{zdf\_tke} namelist.
332
333By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}),
334$P_{LC}$ is assumed to be :
335$336P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}} 337$
338where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth.
339With no information about the wave field, $w_{LC}$ is assumed to be proportional to
340the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module
341\footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as
342  $u_s = 0.016 \,|U_{10m}|$.
343  Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
344  $1.5~10^{-3}$ give the expression used of $u_s$ as a function of the module of surface stress
345}.
346For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at
347a finite depth $H_{LC}$ (which is often close to the mixed layer depth),
348and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures).
349The resulting expression for $w_{LC}$ is :
350$351 w_{LC} = 352 \begin{cases} 353 c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} ) & \text{if -z \leq H_{LC}} \\ 354 0 & \text{otherwise} 355 \end{cases} 356$
357where $c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data.
358The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second.
359The value of $c_{LC}$ is set through the \np{rn_lc}{rn\_lc} namelist parameter,
360having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}.
361
362The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations:
363$H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by
364converting its kinetic energy to potential energy, according to
365$366- \int_{-H_{LC}}^0 { N^2\;\;dz} = \frac{1}{2} u_s^2 367$
368
369%% =================================================================================================
370\subsubsection{Mixing just below the mixed layer}
371
372Vertical mixing parameterizations commonly used in ocean general circulation models tend to
373produce mixed-layer depths that are too shallow during summer months and windy conditions.
374This bias is particularly acute over the Southern Ocean.
375To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{rodgers.aumont.ea_B14}.
376The parameterization is an empirical one, \ie\ not derived from theoretical considerations,
377but rather is meant to account for observed processes that affect the density structure of
378the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme
379(\ie\ near-inertial oscillations and ocean swells and waves).
380
381When using this parameterization (\ie\ when \np[=1]{nn_etau}{nn\_etau}),
382the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations,
383swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition,
384plus a depth depend one given by:
385\begin{equation}
386  \label{eq:ZDF_Ehtau}
387  S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau}
388\end{equation}
389where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that
390penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of
391the penetration, and $f_i$ is the ice concentration
392(no penetration if $f_i=1$, \ie\ if the ocean is entirely covered by sea-ice).
393The value of $f_r$, usually a few percents, is specified through \np{rn_efr}{rn\_efr} namelist parameter.
394The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np[=0]{nn_etau}{nn\_etau}) or
395a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m at high latitudes
396(\np[=1]{nn_etau}{nn\_etau}).
397
398Note that two other option exist, \np[=2, 3]{nn_etau}{nn\_etau}.
399They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer,
400or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean.
401Those two options are obsolescent features introduced for test purposes.
402They will be removed in the next release.
403
404% This should be explain better below what this rn_eice parameter is meant for:
405In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn_eice}{rn\_eice} namelist parameter.
406This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration
407is greater than 25\%.
408
409% from Burchard et al OM 2008 :
410% the most critical process not reproduced by statistical turbulence models is the activity of
411% internal waves and their interaction with turbulence. After the Reynolds decomposition,
412% internal waves are in principle included in the RANS equations, but later partially
413% excluded by the hydrostatic assumption and the model resolution.
414% Thus far, the representation of internal wave mixing in ocean models has been relatively crude
415% (\eg\ Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
416
417%% =================================================================================================
418\subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls})]{GLS: Generic Length Scale (\protect\np{ln_zdfgls}{ln\_zdfgls})}
419\label{subsec:ZDF_gls}
420
421\begin{listing}
422  \nlst{namzdf_gls}
423  \caption{\forcode{&namzdf_gls}}
424  \label{lst:namzdf_gls}
425\end{listing}
426
427The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations:
428one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale,
429$\psi$ \citep{umlauf.burchard_JMR03, umlauf.burchard_CSR05}.
430This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$,
431where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:ZDF_GLS} allows to recover a number of
432well-known turbulent closures ($k$-$kl$ \citep{mellor.yamada_RG82}, $k$-$\epsilon$ \citep{rodi_JGR87},
433$k$-$\omega$ \citep{wilcox_AJ88} among others \citep{umlauf.burchard_JMR03,kantha.carniel_JMR03}).
434The GLS scheme is given by the following set of equations:
435\begin{equation}
436  \label{eq:ZDF_gls_e}
437  \frac{\partial \bar{e}}{\partial t} =
438  \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
439      +\left( \frac{\partial v}{\partial k} \right)^2} \right]
440  -K_\rho \,N^2
441  +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
442  - \epsilon
443\end{equation}
444
445$446 % \label{eq:ZDF_gls_psi} 447 \begin{split} 448 \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{ 449 \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2 450 +\left( \frac{\partial v}{\partial k} \right)^2} \right] 451 - C_3 \,K_\rho\,N^2 - C_2 \,\epsilon \,Fw \right\} \\ 452 &+\frac{1}{e_3\;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 } 453 \;\frac{\partial \psi}{\partial k}} \right]\; 454 \end{split} 455$
456
457$458 % \label{eq:ZDF_gls_kz} 459 \begin{split} 460 K_m &= C_{\mu} \ \sqrt {\bar{e}} \ l \\ 461 K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l 462 \end{split} 463$
464
465$466 % \label{eq:ZDF_gls_eps} 467 {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \; 468$
469where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and
470$\epsilon$ the dissipation rate.
471The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of
472the choice of the turbulence model.
473Four different turbulent models are pre-defined (\autoref{tab:ZDF_GLS}).
474They are made available through the \np{nn_clo}{nn\_clo} namelist parameter.
475
476\begin{table}[htbp]
477  \centering
478  % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
479  \begin{tabular}{ccccc}
480    &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\
481    % & \citep{mellor.yamada_RG82} &  \citep{rodi_JGR87}       & \citep{wilcox_AJ88} &                 \\
482    \hline
483    \hline
484    \np{nn_clo}{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\
485    \hline
486    $( p , n , m )$         &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
487    $\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
488    $\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
489    $C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
490    $C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
491    $C_3$              &      1.           &     1.              &      1.                &       1.           \\
492    $F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
493    \hline
494    \hline
495  \end{tabular}
496  \caption[Set of predefined GLS parameters or equivalently predefined turbulence models available]{
497    Set of predefined GLS parameters, or equivalently predefined turbulence models available with
498    \protect\np[=.true.]{ln_zdfgls}{ln\_zdfgls} and controlled by
499    the \protect\np{nn_clos}{nn\_clos} namelist variable in \protect\nam{zdf_gls}{zdf\_gls}.}
500  \label{tab:ZDF_GLS}
501\end{table}
502
503In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of
504the mixing length towards $\kappa z_b$ ($\kappa$ is the Von Karman constant and $z_b$ the rugosity length scale) value near physical boundaries
505(logarithmic boundary layer law).
506$C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88},
507or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01}
508(\np[=0, 3]{nn_stab_func}{nn\_stab\_func}, resp.).
509The value of $C_{0\mu}$ depends on the choice of the stability function.
510
511The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or
512Neumann condition through \np{nn_bc_surf}{nn\_bc\_surf} and \np{nn_bc_bot}{nn\_bc\_bot}, resp.
513As for TKE closure, the wave effect on the mixing is considered when
514\np[ > 0.]{rn_crban}{rn\_crban} \citep{craig.banner_JPO94, mellor.blumberg_JPO04}.
515The \np{rn_crban}{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and
516\np{rn_charn}{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}.
517
518The $\psi$ equation is known to fail in stably stratified flows, and for this reason
519almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy.
520With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$.
521A value of $c_{lim} = 0.53$ is often used \citep{galperin.kantha.ea_JAS88}.
522\cite{umlauf.burchard_CSR05} show that the value of the clipping factor is of crucial importance for
523the entrainment depth predicted in stably stratified situations,
524and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes.
525The clipping is only activated if \np[=.true.]{ln_length_lim}{ln\_length\_lim},
526and the $c_{lim}$ is set to the \np{rn_clim_galp}{rn\_clim\_galp} value.
527
528The time and space discretization of the GLS equations follows the same energetic consideration as for
529the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}.
530Evaluation of the 4 GLS turbulent closure schemes can be found in \citet{warner.sherwood.ea_OM05} in ROMS model and
531 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model.
532
533% -------------------------------------------------------------------------------------------------------------
534%        OSM OSMOSIS BL Scheme
535% -------------------------------------------------------------------------------------------------------------
536\subsection[OSM: OSMOSIS boundary layer scheme (\forcode{ln_zdfosm = .true.})]
537{OSM: OSMOSIS boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})}
538\label{subsec:ZDF_osm}
539
540\begin{listing}
541  \nlst{namzdf_osm}
542  \caption{\forcode{&namzdf_osm}}
543  \label{lst:namzdf_osm}
544\end{listing}
545
546%--------------------------------------------------------------------------------------------------------------
547\paragraph{Namelist choices}
548Most of the namelist options refer to how to specify the Stokes
549surface drift and penetration depth. There are three options:
550\begin{description}
551  \item \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave} Default value in \texttt{namelist\_ref}. In this case the Stokes drift is
552      assumed to be parallel to the surface wind stress, with
553      magnitude consistent with a constant turbulent Langmuir number
554    $\mathrm{La}_t=$ \protect\np{rn_m_la} {rn\_m\_la} i.e.\
555    $u_{s0}=\tau/(\mathrm{La}_t^2\rho_0)$.  Default value of
556    \protect\np{rn_m_la}{rn\_m\_la} is 0.3. The Stokes penetration
557      depth $\delta =$ \protect\np{rn_osm_dstokes}{rn\_osm\_dstokes}; this has default value
558      of 5~m.
559
560  \item \protect\np[=1]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is
561      assumed to be parallel to the surface wind stress, with
562      magnitude as in the classical Pierson-Moskowitz wind-sea
563      spectrum.  Significant wave height and
564      wave-mean period taken from this spectrum are used to calculate the Stokes penetration
565      depth, following the approach set out in  \citet{breivik.janssen.ea_JPO14}.
566
567    \item \protect\np[=2]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is
568      taken from  ECMWF wave model output, though only the component parallel
569      to the wind stress is retained. Significant wave height and
570      wave-mean period from ECMWF wave model output are used to calculate the Stokes penetration
571      depth, again following \citet{breivik.janssen.ea_JPO14}.
572
573    \end{description}
574
575    Others refer to the treatment of diffusion and viscosity beneath
576    the surface boundary layer:
577\begin{description}
578   \item \protect\np{ln_kpprimix} {ln\_kpprimix}  Default is \np{.true.}. Switches on KPP-style Ri \#-dependent
579     mixing below the surface boundary layer. If this is set
580     \texttt{.true.}  the following variable settings are honoured:
581    \item \protect\np{rn_riinfty}{rn\_riinfty} Critical value of local Ri \# below which
582      shear instability increases vertical mixing from background value.
583    \item \protect\np{rn_difri} {rn\_difri} Maximum value of Ri \#-dependent mixing at $\mathrm{Ri}=0$.
584    \item \protect\np{ln_convmix}{ln\_convmix} If \texttt{.true.} then, where water column is unstable, specify
585       diffusivity equal to \protect\np{rn_dif_conv}{rn\_dif\_conv} (default value is 1 m~s$^{-2}$).
586 \end{description}
587 Diagnostic output is controlled by:
588  \begin{description}
589    \item \protect\np{ln_dia_osm}{ln\_dia\_osm} Default is \np{.false.}; allows XIOS output of OSMOSIS internal fields.
590  \end{description}
591Obsolete namelist parameters include:
592\begin{description}
593   \item \protect\np{ln_use_osm_la}\np{ln\_use\_osm\_la} With \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave},
594      \protect\np{rn_osm_dstokes} {rn\_osm\_dstokes} is always used to specify the Stokes
595      penetration depth.
596   \item \protect\np{nn_ave} {nn\_ave} Choice of averaging method for KPP-style Ri \#
597      mixing. Not taken account of.
598   \item \protect\np{rn_osm_hbl0} {rn\_osm\_hbl0} Depth of initial boundary layer is now set
599     by a density criterion similar to that used in calculating \emph{hmlp} (output as \texttt{mldr10\_1}) in \mdl{zdfmxl}.
600\end{description}
601
602\subsubsection{Summary}
603Much of the time the turbulent motions in the ocean surface boundary
604layer (OSBL) are not given by
605classical shear turbulence. Instead they are in a regime known as
606Langmuir turbulence',  dominated by an
607interaction between the currents and the Stokes drift of the surface waves \citep[e.g.][]{mcwilliams.ea_JFM97}.
608This regime is characterised by strong vertical turbulent motion, and appears when the surface Stokes drift $u_{s0}$ is much greater than the friction velocity $u_{\ast}$. More specifically Langmuir turbulence is thought to be crucial where the turbulent Langmuir number $\mathrm{La}_{t}=(u_{\ast}/u_{s0}) > 0.4$.
609
610The OSMOSIS model is fundamentally based on results of Large Eddy
611Simulations (LES) of Langmuir turbulence and aims to fully describe
612this Langmuir regime. The description in this section is of necessity incomplete and further details are available in Grant. A (2019); in prep.
613
614The OSMOSIS turbulent closure scheme is a similarity-scale scheme in
615the same spirit as the K-profile
616parameterization (KPP) scheme of \citet{large.ea_RG97}.
617A specified shape of diffusivity, scaled by the (OSBL) depth
618$h_{\mathrm{BL}}$ and a turbulent velocity scale, is imposed throughout the
619boundary layer
620$-h_{\mathrm{BL}}<z<\eta$. The turbulent closure model
621also includes fluxes of tracers and momentum that are`non-local'' (independent of the local property gradient).
622
623Rather than the OSBL
624depth being diagnosed in terms of a bulk Richardson number criterion,
625as in KPP, it is set by a prognostic equation that is informed by
626energy budget considerations reminiscent of the classical mixed layer
627models of \citet{kraus.turner_tellus67}.
628The model also includes an explicit parametrization of the structure
629of the pycnocline (the stratified region at the bottom of the OSBL).
630
631Presently, mixing below the OSBL is handled by the Richardson
632number-dependent mixing scheme used in \citet{large.ea_RG97}.
633
634Convective parameterizations such as described in \ref{sec:ZDF_conv}
635below should not be used with the OSMOSIS-OBL model: instabilities
636within the OSBL are part of the model, while instabilities below the
637ML are handled by the Ri \# dependent scheme.
638
639\subsubsection{Depth and velocity scales}
640The model supposes a boundary layer of thickness $h_{\mathrm{bl}}$ enclosing a well-mixed layer of thickness $h_{\mathrm{ml}}$ and a relatively thin pycnocline at the base of thickness $\Delta h$; Fig.~\ref{fig: OSBL_structure} shows typical (a) buoyancy structure and (b) turbulent buoyancy flux profile for the unstable boundary layer (losing buoyancy at the surface; e.g.\ cooling).
641\begin{figure}[!t]
642  \begin{center}
643    %\includegraphics[width=0.7\textwidth]{ZDF_OSM_structure_of_OSBL}
644    \caption{
645      \protect\label{fig: OSBL_structure}
646     The structure of the entraining  boundary layer. (a) Mean buoyancy profile. (b) Profile of the buoyancy flux.
647    }
648  \end{center}
649\end{figure}
650The pycnocline in the OSMOSIS scheme is assumed to have a finite thickness, and may include a number of model levels. This means that the OSMOSIS scheme must parametrize both the thickness of the pycnocline, and the turbulent fluxes within the pycnocline.
651
652Consideration of the power input by wind acting on the Stokes drift suggests that the Langmuir turbulence has velocity scale:
653\begin{equation}\label{eq:w_La}
654w_{*L}= \left(u_*^2 u_{s\,0}\right)^{1/3};
655\end{equation}
656but at times the Stokes drift may be weak due to e.g.\ ice cover, short fetch, misalignment with the surface stress, etc.\ so  a composite velocity scale is assumed for the stable (warming) boundary layer:
657\begin{equation}\label{eq:composite-nu}
658  \nu_{\ast}= \left\{ u_*^3 \left[1-\exp(-.5 \mathrm{La}_t^2)\right]+w_{*L}^3\right\}^{1/3}.
659\end{equation}
660For the unstable boundary layer this is merged with the standard convective velocity scale $w_{*C}=\left(\overline{w^\prime b^\prime}_0 \,h_\mathrm{ml}\right)^{1/3}$, where $\overline{w^\prime b^\prime}_0$ is the upwards surface buoyancy flux, to give:
661\begin{equation}\label{eq:vel-scale-unstable}
662\omega_* = \left(\nu_*^3 + 0.5 w_{*C}^3\right)^{1/3}.
663\end{equation}
664
665\subsubsection{The flux gradient model}
666The flux-gradient relationships used in the OSMOSIS scheme take the form:
667%
669\overline{w^\prime\chi^\prime}=-K\frac{\partial\overline{\chi}}{\partial z} + N_{\chi,s} +N_{\chi,b} +N_{\chi,t},
670\end{equation}
671%
672where $\chi$ is a general variable and $N_{\chi,s}, N_{\chi,b} \mathrm{and} N_{\chi,t}$  are the non-gradient terms, and represent the effects of the different terms in the turbulent flux-budget on the transport of $\chi$. $N_{\chi,s}$ represents the effects that the Stokes shear has on the transport of $\chi$, $N_{\chi,b}$  the effect of buoyancy, and $N_{\chi,t}$ the effect of the turbulent transport.  The same general form for the flux-gradient relationship is used to parametrize the transports of momentum, heat and salinity.
673
674In terms of the non-dimensionalized depth variables
675%
676\begin{equation}\label{eq:sigma}
677\sigma_{\mathrm{ml}}= -z/h_{\mathrm{ml}}; \;\sigma_{\mathrm{bl}}= -z/h_{\mathrm{bl}},
678\end{equation}
679%
680in unstable conditions the eddy diffusivity ($K_d$) and eddy viscosity ($K_\nu$) profiles are parametrized as:
681%
682\begin{align}\label{eq:diff-unstable}
683K_d=&0.8\, \omega_*\, h_{\mathrm{ml}} \, \sigma_{\mathrm{ml}} \left(1-\beta_d \sigma_{\mathrm{ml}}\right)^{3/2}
684\\\label{eq:visc-unstable}
685K_\nu =& 0.3\, \omega_* \,h_{\mathrm{ml}}\, \sigma_{\mathrm{ml}} \left(1-\beta_\nu \sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right)
686\end{align}
687%
688where $\beta_d$ and $\beta_\nu$ are parameters that are determined by matching Eqs \ref{eq:diff-unstable} and \ref{eq:visc-unstable} to the eddy diffusivity and viscosity at the base of the well-mixed layer, given by
689%
690\begin{equation}\label{eq:diff-wml-base}
691K_{d,\mathrm{ml}}=K_{\nu,\mathrm{ml}}=\,0.16\,\omega_* \Delta h.
692\end{equation}
693%
694For stable conditions the eddy diffusivity/viscosity profiles are given by:
695%
696\begin{align}\label{diff-stable}
697K_d= & 0.75\,\, \nu_*\, h_{\mathrm{ml}}\,\,  \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right]\sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)^{3/2} \\\label{eq:visc-stable}
698K_\nu = & 0.375\,\,  \nu_*\, h_{\mathrm{ml}} \,\, \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right] \sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right).
699\end{align}
700%
701The shape of the eddy viscosity and diffusivity profiles is the same as the shape in the unstable OSBL. The eddy diffusivity/viscosity depends on the stability parameter $h_{\mathrm{bl}}/{L_L}$ where $L_L$ is analogous to the Obukhov length, but for Langmuir turbulence:
702\begin{equation}\label{eq:L_L}
703  L_L=-w_{*L}^3/\left<\overline{w^\prime b^\prime}\right>_L,
704\end{equation}
705with the mean turbulent buoyancy flux averaged over the boundary layer given in terms of its surface value $\overline{w^\prime b^\prime}_0$ and (downwards) )solar irradiance $I(z)$ by
706\begin{equation} \label{eq:stable-av-buoy-flux}
707\left<\overline{w^\prime b^\prime}\right>_L = \tfrac{1}{2} {\overline{w^\prime b^\prime}}_0-g\alpha_E\left[\tfrac{1}{2}(I(0)+I(-h))-\left<I\right>\right].
708\end{equation}
709%
710In unstable conditions the eddy diffusivity and viscosity depend on stability through the velocity scale $\omega_*$, which depends on the two velocity scales $\nu_*$ and $w_{*C}$.
711
712Details of the non-gradient terms in \eqref{eq:flux-grad-gen} and of the fluxes within the pycnocline $-h_{\mathrm{bl}}<z<h_{\mathrm{ml}}$ can be found in Grant (2019).
713
714\subsubsection{Evolution of the boundary layer depth}
715
716The prognostic equation for the depth of the neutral/unstable boundary layer is given by \citep{grant+etal18},
717
718\begin{equation} \label{eq:dhdt-unstable}
719%\frac{\partial h_\mathrm{bl}}{\partial t} + \mathbf{U}_b\cdot\nabla h_\mathrm{bl}= W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}}
720\frac{\partial h_\mathrm{bl}}{\partial t} = W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}}
721\end{equation}
722where $h_\mathrm{bl}$ is the horizontally-varying depth of the OSBL,
723$\mathbf{U}_b$ and $W_b$ are the mean horizontal and vertical
724velocities at the base of the OSBL, ${\overline{w^\prime 725 b^\prime}}_\mathrm{ent}$ is the buoyancy flux due to entrainment
726and $\Delta B_\mathrm{bl}$ is the difference between the buoyancy
727averaged over the depth of the OSBL (i.e.\ including the ML and
728pycnocline) and the buoyancy just below the base of the OSBL. This
729equation for the case when the pycnocline has a finite thickness,
730based on the potential energy budget of the OSBL, is the leading term
731\citep{grant+etal18} of a generalization of that used in mixed-layer
732models e.g.\ \citet{kraus.turner_tellus67}, in which the thickness of the pycnocline is taken to be zero.
733
734The entrainment flux for the combination of convective and Langmuir turbulence is given by
735\begin{equation} \label{eq:entrain-flux}
736  {\overline{w^\prime b^\prime}}_\mathrm{ent} = -\alpha_{\mathrm{B}} {\overline{w^\prime b^\prime}}_0 - \alpha_{\mathrm{S}} \frac{u_*^3}{h_{\mathrm{ml}}}
737  + G\left(\delta/h_{\mathrm{ml}} \right)\left[\alpha_{\mathrm{S}}e^{-1.5\, \mathrm{La}_t}-\alpha_{\mathrm{L}} \frac{w_{\mathrm{*L}}^3}{h_{\mathrm{ml}}}\right]
738\end{equation}
739where the factor $G\equiv 1 - \mathrm{e}^ {-25\delta/h_{\mathrm{bl}}}(1-4\delta/h_{\mathrm{bl}})$ models the lesser efficiency of Langmuir mixing when the boundary-layer depth is much greater than the Stokes depth, and $\alpha_{\mathrm{B}}$, $\alpha_{S}$  and $\alpha_{\mathrm{L}}$ depend on the ratio of the appropriate eddy turnover time to the inertial timescale $f^{-1}$. Results from the LES suggest $\alpha_{\mathrm{B}}=0.18 F(fh_{\mathrm{bl}}/w_{*C})$, $\alpha_{S}=0.15 F(fh_{\mathrm{bl}}/u_*$  and $\alpha_{\mathrm{L}}=0.035 F(fh_{\mathrm{bl}}/u_{*L})$, where $F(x)\equiv\tanh(x^{-1})^{0.69}$.
740
741For the stable boundary layer, the equation for the depth of the OSBL is:
742
743\begin{equation}\label{eq:dhdt-stable}
744\max\left(\Delta B_{bl},\frac{w_{*L}^2}{h_\mathrm{bl}}\right)\frac{\partial h_\mathrm{bl}}{\partial t} = \left(0.06 + 0.52\,\frac{ h_\mathrm{bl}}{L_L}\right) \frac{w_{*L}^3}{h_\mathrm{bl}} +\left<\overline{w^\prime b^\prime}\right>_L.
745\end{equation}
746
747Equation. \ref{eq:dhdt-unstable} always leads to the depth of the entraining OSBL increasing (ignoring the effect of the mean vertical motion), but the change in the thickness of the stable OSBL given by Eq. \ref{eq:dhdt-stable} can be positive or negative, depending on the magnitudes of $\left<\overline{w^\prime b^\prime}\right>_L$ and $h_\mathrm{bl}/L_L$. The rate at which the depth of the OSBL can decrease is limited by choosing an effective buoyancy $w_{*L}^2/h_\mathrm{bl}$, in place of $\Delta B_{bl}$ which will be $\approx 0$ for the collapsing OSBL.
748
749
750%% =================================================================================================
751\subsection[ Discrete energy conservation for TKE and GLS schemes]{Discrete energy conservation for TKE and GLS schemes}
752\label{subsec:ZDF_tke_ene}
753
754\begin{figure}[!t]
755  \centering
756  \includegraphics[width=0.66\textwidth]{ZDF_TKE_time_scheme}
757  \caption[Subgrid kinetic energy integration in GLS and TKE schemes]{
758    Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and
759    its links to the momentum and tracer time integration.}
760  \label{fig:ZDF_TKE_time_scheme}
761\end{figure}
762
763The production of turbulence by vertical shear (the first term of the right hand side of
764\autoref{eq:ZDF_tke_e}) and  \autoref{eq:ZDF_gls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion
765(first line in \autoref{eq:MB_zdf}).
766To do so a special care has to be taken for both the time and space discretization of
767the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}.
768
769Let us first address the time stepping issue. \autoref{fig:ZDF_TKE_time_scheme} shows how
770the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with
771the one-level forward time stepping of the equation for $\bar{e}$.
772With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to
773the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and
774summing the result vertically:
775\begin{equation}
776  \label{eq:ZDF_energ1}
777  \begin{split}
778    \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
779    &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}
780    - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
781  \end{split}
782\end{equation}
783Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$,
784known at time $t$ (\autoref{fig:ZDF_TKE_time_scheme}), as it is required when using the TKE scheme
785(see \autoref{sec:TD_forward_imp}).
786The first term of the right hand side of \autoref{eq:ZDF_energ1} represents the kinetic energy transfer at
787the surface (atmospheric forcing) and at the bottom (friction effect).
788The second term is always negative.
789It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$.
790\autoref{eq:ZDF_energ1} implies that, to be energetically consistent,
791the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as
792${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$
793(and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).
794
795A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification
796(second term of the right hand side of \autoref{eq:ZDF_tke_e} and \autoref{eq:ZDF_gls_e}).
797This term must balance the input of potential energy resulting from vertical mixing.
798The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by
799multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically:
800\begin{equation}
801  \label{eq:ZDF_energ2}
802  \begin{split}
803    \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
804    &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta}
805    - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
806    &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
807    + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
808  \end{split}
809\end{equation}
810where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.
811The first term of the right hand side of \autoref{eq:ZDF_energ2} is always zero because
812there is no diffusive flux through the ocean surface and bottom).
813The second term is minus the destruction rate of  $\bar{e}$ due to stratification.
814Therefore \autoref{eq:ZDF_energ1} implies that, to be energetically consistent,
815the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:ZDF_tke_e} and  \autoref{eq:ZDF_gls_e}.
816
817Let us now address the space discretization issue.
818The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in
819the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:DOM_cell}).
820A space averaging is thus required to obtain the shear TKE production term.
821By redoing the \autoref{eq:ZDF_energ1} in the 3D case, it can be shown that the product of eddy coefficient by
822the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
823Furthermore, the time variation of $e_3$ has be taken into account.
824
825The above energetic considerations leads to the following final discrete form for the TKE equation:
826\begin{equation}
827  \label{eq:ZDF_tke_ene}
828  \begin{split}
829    \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv
830    \Biggl\{ \Biggr.
831    &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} }
832        \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
833    +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} }
834        \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j}
835    \Biggr. \Biggr\}   \\
836    %
837    - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
838    %
839    +&\frac{1}{{e_3w}^{t+\rdt}}  \;\delta_{k+1/2} \left[   {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}}   \right]   \\
840    %
841    - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
842  \end{split}
843\end{equation}
844where the last two terms in \autoref{eq:ZDF_tke_ene} (vertical diffusion and Kolmogorov dissipation)
845are time stepped using a backward scheme (see\autoref{sec:TD_forward_imp}).
846Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible.
847%The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as
848%they all appear in the right hand side of \autoref{eq:ZDF_tke_ene}.
849%For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.
850
851%% =================================================================================================
852\section{Convection}
853\label{sec:ZDF_conv}
854
855Static instabilities (\ie\ light potential densities under heavy ones) may occur at particular ocean grid points.
856In nature, convective processes quickly re-establish the static stability of the water column.
857These processes have been removed from the model via the hydrostatic assumption so they must be parameterized.
858Three parameterisations are available to deal with convective processes:
859a non-penetrative convective adjustment or an enhanced vertical diffusion,
860or/and the use of a turbulent closure scheme.
861
862%% =================================================================================================
863\subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc})]{Non-penetrative convective adjustment (\protect\np{ln_tranpc}{ln\_tranpc})}
864\label{subsec:ZDF_npc}
865
866\begin{figure}[!htb]
867  \centering
868  \includegraphics[width=0.66\textwidth]{ZDF_npc}
869  \caption[Unstable density profile treated by the non penetrative convective adjustment algorithm]{
870    Example of an unstable density profile treated by
871    the non penetrative convective adjustment algorithm.
872    $1^{st}$ step: the initial profile is checked from the surface to the bottom.
873    It is found to be unstable between levels 3 and 4.
874    They are mixed.
875    The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed.
876    The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed.
877    The $1^{st}$ step ends since the density profile is then stable below the level 3.
878    $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step:
879    levels 2 to 5 are mixed.
880    The new density profile is checked.
881    It is found stable: end of algorithm.}
882  \label{fig:ZDF_npc}
883\end{figure}
884
885Options are defined through the \nam{zdf}{zdf} namelist variables.
886The non-penetrative convective adjustment is used when \np[=.true.]{ln_zdfnpc}{ln\_zdfnpc}.
887It is applied at each \np{nn_npc}{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of
888the water column, but only until the density structure becomes neutrally stable
889(\ie\ until the mixed portion of the water column has \textit{exactly} the density of the water just below)
891The associated algorithm is an iterative process used in the following way (\autoref{fig:ZDF_npc}):
892starting from the top of the ocean, the first instability is found.
893Assume in the following that the instability is located between levels $k$ and $k+1$.
894The temperature and salinity in the two levels are vertically mixed, conserving the heat and salt contents of
895the water column.
896The new density is then computed by a linear approximation.
897If the new density profile is still unstable between levels $k+1$ and $k+2$,
898levels $k$, $k+1$ and $k+2$ are then mixed.
899This process is repeated until stability is established below the level $k$
900(the mixing process can go down to the ocean bottom).
901The algorithm is repeated to check if the density profile between level $k-1$ and $k$ is unstable and/or
902if there is no deeper instability.
903
904This algorithm is significantly different from mixing statically unstable levels two by two.
905The latter procedure cannot converge with a finite number of iterations for some vertical profiles while
906the algorithm used in \NEMO\ converges for any profile in a number of iterations which is less than
907the number of vertical levels.
908This property is of paramount importance as pointed out by \citet{killworth_iprc89}:
909it avoids the existence of permanent and unrealistic static instabilities at the sea surface.
910This non-penetrative convective algorithm has been proved successful in studies of the deep water formation in
912
913The current implementation has been modified in order to deal with any non linear equation of seawater
914(L. Brodeau, personnal communication).
915Two main differences have been introduced compared to the original algorithm:
916$(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency
917(not the difference in potential density);
918$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in
919the same way their temperature and salinity has been mixed.
920These two modifications allow the algorithm to perform properly and accurately with TEOS10 or EOS-80 without
921having to recompute the expansion coefficients at each mixing iteration.
922
923%% =================================================================================================
924\subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd})]{Enhanced vertical diffusion (\protect\np{ln_zdfevd}{ln\_zdfevd})}
925\label{subsec:ZDF_evd}
926
927Options are defined through the  \nam{zdf}{zdf} namelist variables.
928The enhanced vertical diffusion parameterisation is used when \np[=.true.]{ln_zdfevd}{ln\_zdfevd}.
929In this case, the vertical eddy mixing coefficients are assigned very large values
930in regions where the stratification is unstable
931(\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}.
932This is done either on tracers only (\np[=0]{nn_evdm}{nn\_evdm}) or
933on both momentum and tracers (\np[=1]{nn_evdm}{nn\_evdm}).
934
935In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np[=1]{nn_evdm}{nn\_evdm},
936the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to
937the namelist parameter \np{rn_avevd}{rn\_avevd}.
938A typical value for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$.
939This parameterisation of convective processes is less time consuming than
940the convective adjustment algorithm presented above when mixing both tracers and
941momentum in the case of static instabilities.
942
943Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$.
944This removes a potential source of divergence of odd and even time step in
945a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:TD_mLF}).
946
947%% =================================================================================================
948\subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf_}\{\forcode{tke,gls,osm}\})]{Handling convection with turbulent closure schemes (\forcode{ln_zdf{tke,gls,osm}})}
949\label{subsec:ZDF_tcs}
950
951The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and
952\autoref{subsec:ZDF_osm} (\ie\ \np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} or \np{ln_zdfosm}{ln\_zdfosm} defined) deal, in theory,
953with statically unstable density profiles.
954In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in
955\autoref{eq:ZDF_tke_e} or \autoref{eq:ZDF_gls_e} becomes a source term, since $N^2$ is negative.
956It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also of the four neighboring values at
957velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$).
958These large values restore the static stability of the water column in a way similar to that of
959the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}).
960However, in the vicinity of the sea surface (first ocean layer), the eddy coefficients computed by
961the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$,
962because the mixing length scale is bounded by the distance to the sea surface.
963It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme,
964\ie\ setting the \np{ln_zdfnpc}{ln\_zdfnpc} namelist parameter to true and
965defining the turbulent closure (\np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} = \forcode{.true.}) all together.
966
967The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection,
968%as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},
969therefore \np[=.false.]{ln_zdfevd}{ln\_zdfevd} should be used with the OSMOSIS scheme.
970% gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
971
972%% =================================================================================================
973\section[Double diffusion mixing (\forcode{ln_zdfddm})]{Double diffusion mixing (\protect\np{ln_zdfddm}{ln\_zdfddm})}
974\label{subsec:ZDF_ddm}
975
976%\nlst{namzdf_ddm}
977
978This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter
979\np{ln_zdfddm}{ln\_zdfddm} in \nam{zdf}{zdf}.
980Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa.
981The former condition leads to salt fingering and the latter to diffusive convection.
982Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean.
983\citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that
984it leads to relatively minor changes in circulation but exerts significant regional influences on
985temperature and salinity.
986
987Diapycnal mixing of S and T are described by diapycnal diffusion coefficients
988\begin{align*}
989  % \label{eq:ZDF_ddm_Kz}
990  &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\
991  &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
992\end{align*}
993where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,
994and $o$ by processes other than double diffusion.
995The rates of double-diffusive mixing depend on the buoyancy ratio
996$R_\rho = \alpha \partial_z T / \beta \partial_z S$, where $\alpha$ and $\beta$ are coefficients of
997thermal expansion and saline contraction (see \autoref{subsec:TRA_eos}).
998To represent mixing of $S$ and $T$ by salt fingering, we adopt the diapycnal diffusivities suggested by Schmitt
999(1981):
1000\begin{align}
1001  \label{eq:ZDF_ddm_f}
1002  A_f^{vS} &=
1003             \begin{cases}
1004               \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
1005               0                              &\text{otherwise}
1006             \end{cases}
1007  \\         \label{eq:ZDF_ddm_f_T}
1008  A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho
1009\end{align}
1010
1011\begin{figure}[!t]
1012  \centering
1013  \includegraphics[width=0.66\textwidth]{ZDF_ddm}
1014  \caption[Diapycnal diffusivities for temperature and salt in regions of salt fingering and
1015  diffusive convection]{
1016    From \citet{merryfield.holloway.ea_JPO99}:
1017    (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in
1018    regions of salt fingering.
1019    Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and
1020    thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$;
1021    (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in
1022    regions of diffusive convection.
1023    Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation.
1024    The latter is not implemented in \NEMO.}
1025  \label{fig:ZDF_ddm}
1026\end{figure}
1027
1028The factor 0.7 in \autoref{eq:ZDF_ddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx 0.7$ of
1029buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}).
1030Following  \citet{merryfield.holloway.ea_JPO99}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
1031
1032To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by
1033Federov (1988) is used:
1034\begin{align}
1035  % \label{eq:ZDF_ddm_d}
1036  A_d^{vT} &=
1037             \begin{cases}
1038               1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
1039               &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
1040               0                       &\text{otherwise}
1041             \end{cases}
1042                                       \nonumber \\
1043  \label{eq:ZDF_ddm_d_S}
1044  A_d^{vS} &=
1045             \begin{cases}
1046               A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
1047               A_d^{vT} \ 0.15 \ R_\rho               &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
1048               0                       &\text{otherwise}
1049             \end{cases}
1050\end{align}
1051
1052The dependencies of \autoref{eq:ZDF_ddm_f} to \autoref{eq:ZDF_ddm_d_S} on $R_\rho$ are illustrated in
1053\autoref{fig:ZDF_ddm}.
1054Implementing this requires computing $R_\rho$ at each grid point on every time step.
1055This is done in \mdl{eosbn2} at the same time as $N^2$ is computed.
1056This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive).
1057
1058%% =================================================================================================
1059\section[Bottom and top friction (\textit{zdfdrg.F90})]{Bottom and top friction (\protect\mdl{zdfdrg})}
1060\label{sec:ZDF_drg}
1061
1062\begin{listing}
1063  \nlst{namdrg}
1064  \caption{\forcode{&namdrg}}
1065  \label{lst:namdrg}
1066\end{listing}
1067\begin{listing}
1068  \nlst{namdrg_top}
1069  \caption{\forcode{&namdrg_top}}
1070  \label{lst:namdrg_top}
1071\end{listing}
1072\begin{listing}
1073  \nlst{namdrg_bot}
1074  \caption{\forcode{&namdrg_bot}}
1075  \label{lst:namdrg_bot}
1076\end{listing}
1077
1078Options to define the top and bottom friction are defined through the \nam{drg}{drg} namelist variables.
1079The bottom friction represents the friction generated by the bathymetry.
1080The top friction represents the friction generated by the ice shelf/ocean interface.
1081As the friction processes at the top and the bottom are treated in and identical way,
1082the description below considers mostly the bottom friction case, if not stated otherwise.
1083
1084Both the surface momentum flux (wind stress) and the bottom momentum flux (bottom friction) enter the equations as
1085a condition on the vertical diffusive flux.
1086For the bottom boundary layer, one has:
1087 $1088 % \label{eq:ZDF_bfr_flux} 1089 A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U} 1090$
1091where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside
1092the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean).
1093How ${\cal F}_h^{\textbf U}$ influences the interior depends on the vertical resolution of the model near
1094the bottom relative to the Ekman layer depth.
1095For example, in order to obtain an Ekman layer depth $d = \sqrt{2\;A^{vm}} / f = 50$~m,
1096one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$
1097(for a Coriolis frequency $f = 10^{-4}$~m$^2$s$^{-1}$).
1098With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m.
1099When the vertical mixing coefficient is this small, using a flux condition is equivalent to
1100entering the viscous forces (either wind stress or bottom friction) as a body force over the depth of the top or
1101bottom model layer.
1102To illustrate this, consider the equation for $u$ at $k$, the last ocean level:
1103\begin{equation}
1104  \label{eq:ZDF_drg_flux2}
1105  \frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}}
1106\end{equation}
1107If the bottom layer thickness is 200~m, the Ekman transport will be distributed over that depth.
1108On the other hand, if the vertical resolution is high (1~m or less) and a turbulent closure model is used,
1109the turbulent Ekman layer will be represented explicitly by the model.
1110However, the logarithmic layer is never represented in current primitive equation model applications:
1111it is \emph{necessary} to parameterize the flux ${\cal F}^u_h$.
1112Two choices are available in \NEMO: a linear and a quadratic bottom friction.
1113Note that in both cases, the rotation between the interior velocity and the bottom friction is neglected in
1114the present release of \NEMO.
1115
1116In the code, the bottom friction is imposed by adding the trend due to the bottom friction to
1117 the general momentum trend in \mdl{dynzdf}.
1118For the time-split surface pressure gradient algorithm, the momentum trend due to
1119the barotropic component needs to be handled separately.
1120For this purpose it is convenient to compute and store coefficients which can be simply combined with
1121bottom velocities and geometric values to provide the momentum trend due to bottom friction.
1122 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where:
1123\begin{equation}
1124  \label{eq:ZDF_bfr_bdef}
1125  \frac{\partial {\textbf U_h}}{\partial t} =
1126  - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
1127\end{equation}
1128where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
1129Note than from \NEMO\ 4.0, drag coefficients are only computed at cell centers (\ie\ at T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points.
1130
1131%% =================================================================================================
1132\subsection[Linear top/bottom friction (\forcode{ln_lin})]{Linear top/bottom friction (\protect\np{ln_lin}{ln\_lin})}
1133\label{subsec:ZDF_drg_linear}
1134
1135The linear friction parameterisation (including the special case of a free-slip condition) assumes that
1136the friction is proportional to the interior velocity (\ie\ the velocity of the first/last model level):
1137$1138 % \label{eq:ZDF_bfr_linear} 1139 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 1140$
1141where $r$ is a friction coefficient expressed in $m s^{-1}$.
1142This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean,
1143and setting $r = H / \tau$, where $H$ is the ocean depth.
1144Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}.
1145A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used in quasi-geostrophic models.
1146One may consider the linear friction as an approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$
1147(\citet{gill_bk82}, Eq. 9.6.6).
1148For example, with a drag coefficient $C_D = 0.002$, a typical speed of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$,
1149and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.
1150This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days.
1151It can be changed by specifying \np{rn_Uc0}{rn\_Uc0} (namelist parameter).
1152
1153 For the linear friction case the drag coefficient used in the general expression \autoref{eq:ZDF_bfr_bdef} is:
1154$1155 % \label{eq:ZDF_bfr_linbfr_b} 1156 c_b^T = - r 1157$
1158When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}.
1159Setting \np[=.true.]{ln_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.
1160
1161These values are assigned in \mdl{zdfdrg}.
1162Note that there is support for local enhancement of these values via an externally defined 2D mask array
1163(\np[=.true.]{ln_boost}{ln\_boost}) given in the \ifile{bfr\_coef} input NetCDF file.
1164The mask values should vary from 0 to 1.
1165Locations with a non-zero mask value will have the friction coefficient increased by
1166$mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.
1167
1168%% =================================================================================================
1169\subsection[Non-linear top/bottom friction (\forcode{ln_non_lin})]{Non-linear top/bottom friction (\protect\np{ln_non_lin}{ln\_non\_lin})}
1170\label{subsec:ZDF_drg_nonlinear}
1171
1172The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic:
1173$1174 % \label{eq:ZDF_drg_nonlinear} 1175 {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h 1176 }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b 1177$
1178where $C_D$ is a drag coefficient, and $e_b$ a top/bottom turbulent kinetic energy due to tides,
1179internal waves breaking and other short time scale currents.
1180A typical value of the drag coefficient is $C_D = 10^{-3}$.
1181As an example, the CME experiment \citep{treguier_JGR92} uses $C_D = 10^{-3}$ and
1182$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{killworth_JPO92} uses $C_D = 1.4\;10^{-3}$ and
1183$e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.
1184The CME choices have been set as default values (\np{rn_Cd0}{rn\_Cd0} and \np{rn_ke0}{rn\_ke0} namelist parameters).
1185
1186As for the linear case, the friction is imposed in the code by adding the trend due to
1187the friction to the general momentum trend in \mdl{dynzdf}.
1188For the non-linear friction case the term computed in \mdl{zdfdrg} is:
1189$1190 % \label{eq:ZDF_drg_nonlinbfr} 1191 c_b^T = - \; C_D\;\left[ \left(\bar{u_b}^{i}\right)^2 + \left(\bar{v_b}^{j}\right)^2 + e_b \right]^{1/2} 1192$
1193
1194The coefficients that control the strength of the non-linear friction are initialised as namelist parameters:
1195$C_D$= \np{rn_Cd0}{rn\_Cd0}, and $e_b$ =\np{rn_bfeb2}{rn\_bfeb2}.
1196Note that for applications which consider tides explicitly, a low or even zero value of \np{rn_bfeb2}{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array
1197(\np[=.true.]{ln_boost}{ln\_boost}).
1198This works in the same way as for the linear friction case with non-zero masked locations increased by
1199$mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.
1200
1201%% =================================================================================================
1202\subsection[Log-layer top/bottom friction (\forcode{ln_loglayer})]{Log-layer top/bottom friction (\protect\np{ln_loglayer}{ln\_loglayer})}
1203\label{subsec:ZDF_drg_loglayer}
1204
1205In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using
1206a "law of the wall" scaling. This assumes that the model vertical resolution can capture the logarithmic layer which typically occur for layers thinner than 1 m or so.
1207If  \np[=.true.]{ln_loglayer}{ln\_loglayer}, $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness):
1208$1209 C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2 1210$
1211
1212\noindent where $\kappa$ is the von-Karman constant and \np{rn_z0}{rn\_z0} is a roughness length provided via the namelist.
1213
1214The drag coefficient is bounded such that it is kept greater or equal to
1215the base \np{rn_Cd0}{rn\_Cd0} value which occurs where layer thicknesses become large and presumably logarithmic layers are not resolved at all. For stability reason, it is also not allowed to exceed the value of an additional namelist parameter:
1216\np{rn_Cdmax}{rn\_Cdmax}, \ie
1217$1218 rn\_Cd0 \leq C_D \leq rn\_Cdmax 1219$
1220
1221\noindent The log-layer enhancement can also be applied to the top boundary friction if
1222under ice-shelf cavities are activated (\np[=.true.]{ln_isfcav}{ln\_isfcav}).
1223%In this case, the relevant namelist parameters are \np{rn_tfrz0}{rn\_tfrz0}, \np{rn_tfri2}{rn\_tfri2} and \np{rn_tfri2_max}{rn\_tfri2\_max}.
1224
1225%% =================================================================================================
1226\subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np[=.false.]{ln_drgimp}{ln\_drgimp})}
1227\label{subsec:ZDF_drg_stability}
1228
1229Setting \np[=.false.]{ln_drgimp}{ln\_drgimp} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads:
1230
1231At the top (below an ice shelf cavity):
1232$1233 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1234 = c_{t}^{\textbf{U}}\textbf{u}^{n-1}_{t} 1235$
1236
1237At the bottom (above the sea floor):
1238$1239 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1240 = c_{b}^{\textbf{U}}\textbf{u}^{n-1}_{b} 1241$
1242
1243Since this is conditionally stable, some care needs to exercised over the choice of parameters to ensure that the implementation of explicit top/bottom friction does not induce numerical instability.
1244For the purposes of stability analysis, an approximation to \autoref{eq:ZDF_drg_flux2} is:
1245\begin{equation}
1246  \label{eq:ZDF_Eqn_drgstab}
1247  \begin{split}
1248    \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
1249    &= -\frac{ru}{e_{3u}}\;2\rdt\\
1250  \end{split}
1251\end{equation}
1252\noindent where linear friction and a leapfrog timestep have been assumed.
1253To ensure that the friction cannot reverse the direction of flow it is necessary to have:
1254$1255 |\Delta u| < \;|u| 1256$
1257\noindent which, using \autoref{eq:ZDF_Eqn_drgstab}, gives:
1258$1259 r\frac{2\rdt}{e_{3u}} < 1 \qquad \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\ 1260$
1261This same inequality can also be derived in the non-linear bottom friction case if
1262a velocity of 1 m.s$^{-1}$ is assumed.
1263Alternatively, this criterion can be rearranged to suggest a minimum bottom box thickness to ensure stability:
1264$1265 e_{3u} > 2\;r\;\rdt 1266$
1267\noindent which it may be necessary to impose if partial steps are being used.
1268For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then $e_{3u}$ should be greater than 3.6 m.
1269For most applications, with physically sensible parameters these restrictions should not be of concern.
1270But caution may be necessary if attempts are made to locally enhance the bottom friction parameters.
1271To ensure stability limits are imposed on the top/bottom friction coefficients both
1272during initialisation and at each time step.
1273Checks at initialisation are made in \mdl{zdfdrg} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
1274The number of breaches of the stability criterion are reported as well as
1275the minimum and maximum values that have been set.
1276The criterion is also checked at each time step, using the actual velocity, in \mdl{dynzdf}.
1277Values of the friction coefficient are reduced as necessary to ensure stability;
1278these changes are not reported.
1279
1280Limits on the top/bottom friction coefficient are not imposed if the user has elected to
1281handle the friction implicitly (see \autoref{subsec:ZDF_drg_imp}).
1282The number of potential breaches of the explicit stability criterion are still reported for information purposes.
1283
1284%% =================================================================================================
1285\subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np[=.true.]{ln_drgimp}{ln\_drgimp})}
1286\label{subsec:ZDF_drg_imp}
1287
1288An optional implicit form of bottom friction has been implemented to improve model stability.
1289We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting.
1290This option can be invoked by setting \np{ln_drgimp}{ln\_drgimp} to \forcode{.true.} in the \nam{drg}{drg} namelist.
1291%This option requires \np{ln_zdfexp}{ln\_zdfexp} to be \forcode{.false.} in the \nam{zdf}{zdf} namelist.
1292
1293This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step:
1294
1295At the top (below an ice shelf cavity):
1296$1297 % \label{eq:ZDF_dynZDF__drg_top} 1298 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t} 1299 = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t} 1300$
1301
1302At the bottom (above the sea floor):
1303$1304 % \label{eq:ZDF_dynZDF__drg_bot} 1305 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b} 1306 = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b} 1307$
1308
1309where $t$ and $b$ refers to top and bottom layers respectively.
1310Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit.
1311
1312%% =================================================================================================
1313\subsection[Bottom friction with split-explicit free surface]{Bottom friction with split-explicit free surface}
1314\label{subsec:ZDF_drg_ts}
1315
1316With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np[=.false.]{ln_drgimp}{ln\_drgimp} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie\ \forcode{ln_bt_fw=.false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np[=.true.]{ln_drgimp}{ln\_drgimp},  stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions.
1317
1318The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows:
1319\begin{enumerate}
1320\item To extend the stability of the barotropic sub-stepping, bottom stresses are refreshed at each sub-iteration. The baroclinic part of the flow entering the stresses is frozen at the initial time of the barotropic iteration. In case of non-linear friction, the drag coefficient is also constant.
1321\item In case of an implicit drag, specific computations are performed in \mdl{dynzdf} which renders the overall scheme mixed explicit/implicit: the barotropic components of 3d velocities are removed before seeking for the implicit vertical diffusion result. Top/bottom stresses due to the barotropic components are explicitly accounted for thanks to the updated values of barotropic velocities. Then the implicit solution of 3d velocities is obtained. Lastly, the residual barotropic component is replaced by the time split estimate.
1322\end{enumerate}
1323
1324Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration.
1325
1326%% =================================================================================================
1327\section[Internal wave-driven mixing (\forcode{ln_zdfiwm})]{Internal wave-driven mixing (\protect\np{ln_zdfiwm}{ln\_zdfiwm})}
1328\label{subsec:ZDF_tmx_new}
1329
1330\begin{listing}
1331  \nlst{namzdf_iwm}
1332  \caption{\forcode{&namzdf_iwm}}
1333  \label{lst:namzdf_iwm}
1334\end{listing}
1335
1336The parameterization of mixing induced by breaking internal waves is a generalization of
1337the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}.
1338A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,
1339and the resulting diffusivity is obtained as
1340$1341 % \label{eq:ZDF_Kwave} 1342 A^{vT}_{wave} = R_f \,\frac{ \epsilon }{ \rho \, N^2 } 1343$
1344where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of
1345the energy available for mixing.
1346If the \np{ln_mevar}{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and
1347equal to 1/6 \citep{osborn_JPO80}.
1348In the opposite (recommended) case, $R_f$ is instead a function of
1349the turbulence intensity parameter $Re_b = \frac{ \epsilon}{\nu \, N^2}$,
1350with $\nu$ the molecular viscosity of seawater, following the model of \cite{bouffard.boegman_DAO13} and
1351the implementation of \cite{de-lavergne.madec.ea_JPO16}.
1352Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when
1353the mixing efficiency is constant.
1354
1355In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary
1356as a function of $Re_b$ by setting the \np{ln_tsdiff}{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice.
1357This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14},
1358is implemented as in \cite{de-lavergne.madec.ea_JPO16}.
1359
1360The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$,
1361is constructed from three static maps of column-integrated internal wave energy dissipation,
1362$E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures:
1363
1364\begin{align*}
1365  F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
1366  F_{pyc}(i,j,k) &\propto N^{n_p}\\
1367  F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
1368\end{align*}
1369In the above formula, $h_{ab}$ denotes the height above bottom,
1370$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
1371$1372 h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz' } \; , 1373$
1374The $n_p$ parameter (given by \np{nn_zpyc}{nn\_zpyc} in \nam{zdf_iwm}{zdf\_iwm} namelist)
1375controls the stratification-dependence of the pycnocline-intensified dissipation.
1376It can take values of $1$ (recommended) or $2$.
1377Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of
1378the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.
1379$h_{cri}$ is related to the large-scale topography of the ocean (etopo2) and
1380$h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of
1381the abyssal hill topography \citep{goff_JGR10} and the latitude.
1382% Jc: input files names ?
1383
1384%% =================================================================================================
1385\section[Surface wave-induced mixing (\forcode{ln_zdfswm})]{Surface wave-induced mixing (\protect\np{ln_zdfswm}{ln\_zdfswm})}
1386\label{subsec:ZDF_swm}
1387
1388Surface waves produce an enhanced mixing through wave-turbulence interaction.
1389In addition to breaking waves induced turbulence (\autoref{subsec:ZDF_tke}),
1390the influence of non-breaking waves can be accounted introducing
1391wave-induced viscosity and diffusivity as a function of the wave number spectrum.
1392Following \citet{qiao.yuan.ea_OD10}, a formulation of wave-induced mixing coefficient
1393is provided  as a function of wave amplitude, Stokes Drift and wave-number:
1394
1395\begin{equation}
1396  \label{eq:ZDF_Bv}
1397  B_{v} = \alpha {A} {U}_{st} {exp(3kz)}
1398\end{equation}
1399
1400Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude,
1401${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$
1402is a constant which should be determined by observations or
1403numerical experiments and is set to be 1.
1404
1405The coefficient $B_{v}$ is then directly added to the vertical viscosity
1406and diffusivity coefficients.
1407
1408In order to account for this contribution set: \forcode{ln_zdfswm=.true.},
1409then wave interaction has to be activated through \forcode{ln_wave=.true.},
1410the Stokes Drift can be evaluated by setting \forcode{ln_sdw=.true.}
1411(see \autoref{subsec:SBC_wave_sdw})
1412and the needed wave fields can be provided either in forcing or coupled mode
1413(for more information on wave parameters and settings see \autoref{sec:SBC_wave})
1414
1415%% =================================================================================================
1417\label{subsec:ZDF_aimp}
1418
1419The adaptive-implicit vertical advection option in NEMO is based on the work of
1420\citep{shchepetkin_OM15}.  In common with most ocean models, the timestep used with NEMO
1421needs to satisfy multiple criteria associated with different physical processes in order
1422to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical
1423CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the
1424constraints for a range of time and space discretizations and provide the CFL stability
1425criteria for a range of advection schemes. The values for the Leap-Frog with Robert
1426asselin filter time-stepping (as used in NEMO) are reproduced in
1427\autoref{tab:ZDF_zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these
1428restrictions but at the cost of large dispersive errors and, possibly, large numerical
1429viscosity. The adaptive-implicit vertical advection option provides a targetted use of the
1430implicit scheme only when and where potential breaches of the vertical CFL condition
1431occur. In many practical applications these events may occur remote from the main area of
1432interest or due to short-lived conditions such that the extra numerical diffusion or
1433viscosity does not greatly affect the overall solution. With such applications, setting:
1434\forcode{ln_zad_Aimp=.true.} should allow much longer model timesteps to be used whilst
1435retaining the accuracy of the high order explicit schemes over most of the domain.
1436
1437\begin{table}[htbp]
1438  \centering
1439  % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}}
1440  \begin{tabular}{r|ccc}
1441    \hline
1442    spatial discretization  & 2$^nd$ order centered & 3$^rd$ order upwind & 4$^th$ order compact \\
1443    advective CFL criterion &                 0.904 &              0.472  &                0.522 \\
1444    \hline
1445  \end{tabular}
1446  \caption[Advective CFL criteria for the leapfrog with Robert Asselin filter time-stepping]{
1447    The advective CFL criteria for a range of spatial discretizations for
1448    the leapfrog with Robert Asselin filter time-stepping
1449    ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}.}
1451\end{table}
1452
1453In particular, the advection scheme remains explicit everywhere except where and when
1454local vertical velocities exceed a threshold set just below the explicit stability limit.
1455Once the threshold is reached a tapered transition towards an implicit scheme is used by
1456partitioning the vertical velocity into a part that can be treated explicitly and any
1457excess that must be treated implicitly. The partitioning is achieved via a Courant-number
1458dependent weighting algorithm as described in \citep{shchepetkin_OM15}.
1459
1460The local cell Courant number ($Cu$) used for this partitioning is:
1461
1462\begin{equation}
1464  \begin{split}
1465    Cu &= {2 \rdt \over e^n_{3t_{ijk}}} \bigg (\big [ \texttt{Max}(w^n_{ijk},0.0) - \texttt{Min}(w^n_{ijk+1},0.0) \big ]    \\
1466       &\phantom{=} +\big [ \texttt{Max}(e_{{2_u}ij}e^n_{{3_{u}}ijk}u^n_{ijk},0.0) - \texttt{Min}(e_{{2_u}i-1j}e^n_{{3_{u}}i-1jk}u^n_{i-1jk},0.0) \big ]
1467                     \big / e_{{1_t}ij}e_{{2_t}ij}            \\
1468       &\phantom{=} +\big [ \texttt{Max}(e_{{1_v}ij}e^n_{{3_{v}}ijk}v^n_{ijk},0.0) - \texttt{Min}(e_{{1_v}ij-1}e^n_{{3_{v}}ij-1k}v^n_{ij-1k},0.0) \big ]
1469                     \big / e_{{1_t}ij}e_{{2_t}ij} \bigg )    \\
1470  \end{split}
1471\end{equation}
1472
1473\noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as:
1474
1475\begin{align}
1477Cu_{min} &= 0.15 \nonumber \\
1478Cu_{max} &= 0.3  \nonumber \\
1479Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\
1480Fcu    &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\
1481\cf &=
1482     \begin{cases}
1483        0.0                                                        &\text{if $Cu \leq Cu_{min}$} \\
1484        (Cu - Cu_{min})^2 / (Fcu +  (Cu - Cu_{min})^2)             &\text{else if $Cu < Cu_{cut}$} \\
1485        (Cu - Cu_{max}) / Cu                                       &\text{else}
1486     \end{cases}
1487\end{align}
1488
1489\begin{figure}[!t]
1490  \centering
1492  \caption[Partitioning coefficient used to partition vertical velocities into parts]{
1493    The value of the partitioning coefficient (\cf) used to partition vertical velocities into
1494    parts to be treated implicitly and explicitly for a range of typical Courant numbers
1497\end{figure}
1498
1499\noindent The partitioning coefficient is used to determine the part of the vertical
1500velocity that must be handled implicitly ($w_i$) and to subtract this from the total
1501vertical velocity ($w_n$) to leave that which can continue to be handled explicitly:
1502
1503\begin{align}
1505    w_{i_{ijk}} &= \cf_{ijk} w_{n_{ijk}}     \nonumber \\
1506    w_{n_{ijk}} &= (1-\cf_{ijk}) w_{n_{ijk}}
1507\end{align}
1508
1509\noindent Note that the coefficient is such that the treatment is never fully implicit;
1510the three cases from \autoref{eq:ZDF_Eqn_zad_Aimp_partition} can be considered as:
1511fully-explicit; mixed explicit/implicit and mostly-implicit.  With the settings shown the
1512coefficient (\cf) varies as shown in \autoref{fig:ZDF_zad_Aimp_coeff}. Note with these values
1513the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly
1514implicit' is 0.45 which is just below the stability limited given in
1515\autoref{tab:ZDF_zad_Aimp_CFLcrit}  for a 3rd order scheme.
1516
1517The $w_i$ component is added to the implicit solvers for the vertical mixing in
1518\mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}.  This is
1519sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further
1520intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}).
1521For these schemes the implicit upstream fluxes must be added to both the monotonic guess
1522and to the higher order solution when calculating the antidiffusive fluxes. The implicit
1523vertical fluxes are then removed since they are added by the implicit solver later on.
1524
1525The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be
1526used in a wide range of simulations. The following test simulation, however, does illustrate
1527the potential benefits and will hopefully encourage further testing and feedback from users:
1528
1529\begin{figure}[!t]
1530  \centering
1532  \caption[OVERFLOW: time-series of temperature vertical cross-sections]{
1533    A time-series of temperature vertical cross-sections for the OVERFLOW test case.
1534    These results are for the default settings with \forcode{nn_rdt=10.0} and
1537\end{figure}
1538
1539%% =================================================================================================
1540\subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case}
1541
1542The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case}
1543provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case
1544by only a few extra physics choices namely:
1545
1546\begin{verbatim}
1547     ln_dynldf_OFF = .false.
1548     ln_dynldf_lap = .true.
1549     ln_dynldf_hor = .true.
1550     ln_zdfnpc     = .true.
1551     ln_traadv_fct = .true.
1552        nn_fct_h   =  2
1553        nn_fct_v   =  2
1554\end{verbatim}
1555
1556\noindent which were chosen to provide a slightly more stable and less noisy solution. The
1557result when using the default value of \forcode{nn_rdt=10.} without adaptive-implicit
1558vertical velocity is illustrated in \autoref{fig:ZDF_zad_Aimp_overflow_frames}. The mass of
1559cold water, initially sitting on the shelf, moves down the slope and forms a
1560bottom-trapped, dense plume. Even with these extra physics choices the model is close to
1561stability limits and attempts with \forcode{nn_rdt=30.} will fail after about 5.5 hours
1562with excessively high horizontal velocities. This time-scale corresponds with the time the
1563plume reaches the steepest part of the topography and, although detected as a horizontal
1564CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good
1565candidate, therefore, for use of the adaptive-implicit vertical advection scheme.
1566
1567The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps
1568are shown in \autoref{fig:ZDF_zad_Aimp_overflow_all_rdt} (together with the equivalent
1569frames from the base run).  In this simple example the use of the adaptive-implicit
1570vertcal advection scheme has enabled a 12x increase in the model timestep without
1571significantly altering the solution (although at this extreme the plume is more diffuse
1572and has not travelled so far).  Notably, the solution with and without the scheme is
1573slightly different even with \forcode{nn_rdt=10.}; suggesting that the base run was
1574close enough to instability to trigger the scheme despite completing successfully.
1575To assist in diagnosing how active the scheme is, in both location and time, the 3D
1576implicit and explicit components of the vertical velocity are available via XIOS as
1577\texttt{wimp} and \texttt{wexp} respectively.  Likewise, the partitioning coefficient
1578(\cf) is also available as \texttt{wi\_cff}. For a quick oversight of
1579the schemes activity the global maximum values of the absolute implicit component
1580of the vertical velocity and the partitioning coefficient are written to the netCDF
1581version of the run statistics file (\texttt{run.stat.nc}) if this is active (see
1582\autoref{sec:MISC_opt} for activation details).
1583
1584\autoref{fig:ZDF_zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for
1585the various overflow tests.  Note that the adaptive-implicit vertical advection scheme is
1586active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the
1587test case is close to stability limits even with this value. At the larger timesteps, the
1588vertical velocity is treated mostly implicitly at some location throughout the run. The
1589oscillatory nature of this measure appears to be linked to the progress of the plume front
1590as each cusp is associated with the location of the maximum shifting to the adjacent cell.
1591This is illustrated in \autoref{fig:ZDF_zad_Aimp_maxCf_loc} where the i- and k- locations of the
1592maximum have been overlaid for the base run case.
1593
1594\medskip
1595\noindent Only limited tests have been performed in more realistic configurations. In the
1596ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes
1597restartability and reproducibility tests but it is unable to improve the model's stability
1598enough to allow an increase in the model time-step. A view of the time-series of maximum
1599partitioning coefficient (not shown here)  suggests that the default time-step of 5400s is
1600already pushing at stability limits, especially in the initial start-up phase. The
1601time-series does not, however, exhibit any of the 'cuspiness' found with the overflow
1602tests.
1603
1604\medskip
1605\noindent A short test with an eORCA1 configuration promises more since a test using a
1606time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the
1607time-step is limited to 2700s without.
1608
1609\begin{figure}[!t]
1610  \centering
1612  \caption[OVERFLOW: sample temperature vertical cross-sections from mid- and end-run]{
1613    Sample temperature vertical cross-sections from mid- and end-run using
1614    different values for \forcode{nn_rdt} and with or without adaptive implicit vertical advection.
1615    Without the adaptive implicit vertical advection
1616    only the run with the shortest timestep is able to run to completion.
1617    Note also that the colour-scale has been chosen to confirm that
1618    temperatures remain within the original range of 10$^o$ to 20$^o$.}
1620\end{figure}
1621
1622\begin{figure}[!t]
1623  \centering
1625  \caption[OVERFLOW: maximum partitioning coefficient during a series of test runs]{
1626    The maximum partitioning coefficient during a series of test runs with
1627    increasing model timestep length.
1628    At the larger timesteps,
1629    the vertical velocity is treated mostly implicitly at some location throughout the run.}
1631\end{figure}
1632
1633\begin{figure}[!t]
1634  \centering