# source:NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex@11571 Last change on this file since 11571 was 11571, checked in by nicolasmartin, 13 months ago

Cosmetic modifications on section titles

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