source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/doc/latex/NEMO/subfiles/chap_ZDF.tex @ 11263

Last change on this file since 11263 was 11263, checked in by smasson, 15 months ago

dev_r10984_HPC-13 : merge with trunk@11242, see #2285

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