source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex @ 10414

Last change on this file since 10414 was 10414, checked in by nicolasmartin, 2 years ago
  • Comment \label commands on maths environments for unreferenced equations and adapt the unnumbered math container accordingly (mainly switch to shortanded LateX syntax with \[ ... \])
  • Add a code trick to build subfile with its own bibliography
  • Fix right path for main LaTeX document in first line of subfiles (\documentclass[...]{subfiles})
  • Rename abstract_foreword.tex to foreword.tex
  • Fix some non-ASCII codes inserted here or there in LaTeX (\[0-9]*)
  • Made a first iteration on the indentation and alignement within math, figure and table environments to improve source code readability
File size: 79.0 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 ($i.e.$ \key{trabbl} defined,
28see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum
29(see \autoref{sec:ZDF_bfr}).
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 formulation).
36The computation of these coefficients is initialized in the \mdl{zdfini} module and performed in
37the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} 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.
40These 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,
43and thus of the formulation used (see \autoref{chap:STP}).
44
45% -------------------------------------------------------------------------------------------------------------
46%        Constant
47% -------------------------------------------------------------------------------------------------------------
48\subsection{Constant (\protect\key{zdfcst})}
49\label{subsec:ZDF_cst}
50%--------------------------------------------namzdf---------------------------------------------------------
51
52\nlst{namzdf}
53%--------------------------------------------------------------------------------------------------------------
54
55Options are defined through the \ngn{namzdf} namelist variables.
56When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to
57constant values over the whole ocean.
58This is the crudest way to define the vertical ocean physics.
59It is recommended that this option is only used in process studies, not in basin scale simulations.
60Typical values used in this case are:
61\begin{align*}
62  A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}   \\
63  A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
64\end{align*}
65
66These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters.
67In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity,
68that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and
69$\sim10^{-9}~m^2.s^{-1}$ for salinity.
70
71% -------------------------------------------------------------------------------------------------------------
72%        Richardson Number Dependent
73% -------------------------------------------------------------------------------------------------------------
74\subsection{Richardson number dependent (\protect\key{zdfric})}
75\label{subsec:ZDF_ric}
76
77%--------------------------------------------namric---------------------------------------------------------
78
79\nlst{namzdf_ric}
80%--------------------------------------------------------------------------------------------------------------
81
82When \key{zdfric} is defined, a local Richardson number dependent formulation for the vertical momentum and
83tracer eddy coefficients is set through the \ngn{namzdf\_ric} namelist variables.
84The vertical mixing coefficients are diagnosed from the large scale variables computed by the model.
85\textit{In situ} measurements have been used to link vertical turbulent activity to large scale ocean structures.
86The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to
87a dependency between the vertical eddy coefficients and the local Richardson number
88($i.e.$ the ratio of stratification to vertical shear).
89Following \citet{Pacanowski_Philander_JPO81}, the following formulation has been implemented:
90\[
91  % \label{eq:zdfric}
92  \left\{
93    \begin{aligned}
94      A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
95      A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
96    \end{aligned}
97  \right.
98\]
99where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number,
100$N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
101$A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the constant case
102(see \autoref{subsec:ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ is the maximum value that
103can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$.
104The last three values can be modified by setting the \np{rn\_avmri}, \np{rn\_alp} and
105\np{nn\_ric} namelist parameters, respectively.
106
107A simple mixing-layer model to transfer and dissipate the atmospheric forcings
108(wind-stress and buoyancy fluxes) can be activated setting the \np{ln\_mldw}\forcode{ = .true.} in the namelist.
109
110In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ is evaluated and
111the vertical eddy coefficients prescribed within this layer.
112
113This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
114\[
115  h_{e} = Ek \frac {u^{*}} {f_{0}}
116\]
117where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis parameter.
118
119In this similarity height relationship, the turbulent friction velocity:
120\[
121  u^{*} = \sqrt \frac {|\tau|} {\rho_o}
122\]
123is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$.
124The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
125Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to
126the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}.
127
128% -------------------------------------------------------------------------------------------------------------
129%        TKE Turbulent Closure Scheme
130% -------------------------------------------------------------------------------------------------------------
131\subsection{TKE turbulent closure scheme (\protect\key{zdftke})}
132\label{subsec:ZDF_tke}
133
134%--------------------------------------------namzdf_tke--------------------------------------------------
135
136\nlst{namzdf_tke}
137%--------------------------------------------------------------------------------------------------------------
138
139The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on
140a prognostic equation for $\bar{e}$, the turbulent kinetic energy,
141and a closure assumption for the turbulent length scales.
142This turbulent closure model has been developed by \citet{Bougeault1989} in the atmospheric case,
143adapted by \citet{Gaspar1990} for the oceanic case, and embedded in OPA, the ancestor of NEMO,
144by \citet{Blanke1993} for equatorial Atlantic simulations.
145Since then, significant modifications have been introduced by \citet{Madec1998} in both the implementation and
146the formulation of the mixing length scale.
147The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear,
148its destruction through stratification, its vertical diffusion, and its dissipation of \citet{Kolmogorov1942} type:
149\begin{equation}
150  \label{eq:zdftke_e}
151  \frac{\partial \bar{e}}{\partial t} =
152  \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
153      +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
154  -K_\rho\,N^2
155  +\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
156      \;\frac{\partial \bar{e}}{\partial k}} \right]
157  - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
158\end{equation}
159\[
160  % \label{eq:zdftke_kz}
161  \begin{split}
162    K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }    \\
163    K_\rho &= A^{vm} / P_{rt}
164  \end{split}
165\]
166where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
167$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales,
168$P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients.
169The constants $C_k =  0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with
170vertical mixing at any depth \citep{Gaspar1990}.
171They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}.
172$P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function of the local Richardson number, $R_i$:
173\begin{align*}
174  % \label{eq:prt}
175  P_{rt} =
176  \begin{cases}
177    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}   \\
178    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}   \\
179    \ \ 10 &      \text{if $\ 2 \leq R_i$}
180  \end{cases}
181\end{align*}
182Options are defined through the  \ngn{namzdfy\_tke} namelist variables.
183The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable.
184
185At the sea surface, the value of $\bar{e}$ is prescribed from the wind stress field as
186$\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} namelist parameter.
187The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}), however a much larger value can be used when
188taking into account the surface wave breaking (see below Eq. \autoref{eq:ZDF_Esbc}).
189The bottom value of TKE is assumed to be equal to the value of the level just above.
190The time integration of the $\bar{e}$ equation may formally lead to negative values because
191the numerical scheme does not ensure its positivity.
192To overcome this problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} namelist parameter).
193Following \citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.
194This allows the subsequent formulations to match that of \citet{Gargett1984} for the diffusion in
195the thermocline and deep ocean :  $K_\rho = 10^{-3} / N$.
196In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with
197too weak vertical diffusion.
198They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and
199\np{rn\_avt0} (namzdf namelist, see \autoref{subsec:ZDF_cst}).
200
201\subsubsection{Turbulent length scale}
202
203For computational efficiency, the original formulation of the turbulent length scales proposed by
204\citet{Gaspar1990} has been simplified.
205Four formulations are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist parameter.
206The first two are based on the following first order approximation \citep{Blanke1993}:
207\begin{equation}
208  \label{eq:tke_mxl0_1}
209  l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
210\end{equation}
211which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency.
212The resulting length scale is bounded by the distance to the surface or to the bottom
213(\np{nn\_mxl}\forcode{ = 0}) or by the local vertical scale factor (\np{nn\_mxl}\forcode{ = 1}).
214\citet{Blanke1993} notice that this simplification has two major drawbacks:
215it makes no sense for locally unstable stratification and the computation no longer uses all
216the information contained in the vertical density profile.
217To overcome these drawbacks, \citet{Madec1998} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases,
218which add an extra assumption concerning the vertical gradient of the computed length scale.
219So, the length scales are first evaluated as in \autoref{eq:tke_mxl0_1} and then bounded such that:
220\begin{equation}
221  \label{eq:tke_mxl_constraint}
222  \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
223  \qquad \text{with }\  l =  l_k = l_\epsilon
224\end{equation}
225\autoref{eq:tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than
226the variations of depth.
227It provides a better approximation of the \citet{Gaspar1990} formulation while being much less
228time consuming.
229In particular, it allows the length scale to be limited not only by the distance to the surface or
230to the ocean bottom but also by the distance to a strongly stratified portion of the water column such as
231the thermocline (\autoref{fig:mixing_length}).
232In order to impose the \autoref{eq:tke_mxl_constraint} constraint, we introduce two additional length scales:
233$l_{up}$ and $l_{dwn}$, the upward and downward length scales, and
234evaluate the dissipation and mixing length scales as
235(and note that here we use numerical indexing):
236%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
237\begin{figure}[!t]
238  \begin{center}
239    \includegraphics[width=1.00\textwidth]{Fig_mixing_length}
240    \caption{
241      \protect\label{fig:mixing_length}
242      Illustration of the mixing length computation.
243    }
244  \end{center}
245\end{figure}
246%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
247\[
248  % \label{eq:tke_mxl2}
249  \begin{aligned}
250    l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
251    \quad &\text{ from $k=1$ to $jpk$ }\ \\
252    l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)\right)
253    \quad &\text{ from $k=jpk$ to $1$ }\ \\
254  \end{aligned}
255\]
256where $l^{(k)}$ is computed using \autoref{eq:tke_mxl0_1}, $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
257
258In the \np{nn\_mxl}\forcode{ = 2} case, the dissipation and mixing length scales take the same value:
259$ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the \np{nn\_mxl}\forcode{ = 3} case,
260the dissipation and mixing turbulent length scales are give as in \citet{Gaspar1990}:
261\[
262  % \label{eq:tke_mxl_gaspar}
263  \begin{aligned}
264    & l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }   \\
265    & l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)
266  \end{aligned}
267\]
268
269At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist parameter.
270Usually the surface scale is given by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant and
271$z_o$ the roughness parameter of the surface.
272Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} leads to a 0.04~m, the default value of \np{rn\_mxl0}.
273In the ocean interior a minimum length scale is set to recover the molecular viscosity when
274$\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
275
276\subsubsection{Surface wave breaking parameterization}
277%-----------------------------------------------------------------------%
278
279Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified to
280include the effect of surface wave breaking energetics.
281This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow.
282The \citet{Mellor_Blumberg_JPO04} modifications acts on surface length scale and TKE values and
283air-sea drag coefficient.
284The latter concerns the bulk formulea and is not discussed here.
285
286Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is :
287\begin{equation}
288  \label{eq:ZDF_Esbc}
289  \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
290\end{equation}
291where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality which depends on the ''wave age'',
292ranging from 57 for mature waves to 146 for younger waves \citep{Mellor_Blumberg_JPO04}.
293The boundary condition on the turbulent length scale follows the Charnock's relation:
294\begin{equation}
295  \label{eq:ZDF_Lsbc}
296  l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
297\end{equation}
298where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
299\citet{Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by
300\citet{Stacey_JPO99} citing observation evidence, and
301$\alpha_{CB} = 100$ the Craig and Banner's value.
302As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$,
303with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds
304to $\alpha_{CB} = 100$.
305Further setting  \np{ln\_mxl0} to true applies \autoref{eq:ZDF_Lsbc} as surface boundary condition on length scale,
306with $\beta$ hard coded to the Stacey's value.
307Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on
308surface $\bar{e}$ value.
309
310\subsubsection{Langmuir cells}
311%--------------------------------------%
312
313Langmuir circulations (LC) can be described as ordered large-scale vertical motions in
314the surface layer of the oceans.
315Although LC have nothing to do with convection, the circulation pattern is rather similar to
316so-called convective rolls in the atmospheric boundary layer.
317The detailed physics behind LC is described in, for example, \citet{Craik_Leibovich_JFM76}.
318The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and
319wind drift currents.
320
321Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by
322\citep{Axell_JGR02} for a $k-\epsilon$ turbulent closure.
323The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in
324an extra source terms of TKE, $P_{LC}$.
325The presence of $P_{LC}$ in \autoref{eq:zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to
326\forcode{.true.} in the namtke namelist.
327 
328By making an analogy with the characteristic convective velocity scale ($e.g.$, \citet{D'Alessio_al_JPO98}),
329$P_{LC}$ is assumed to be :
330\[
331P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}}
332\]
333where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth.
334With no information about the wave field, $w_{LC}$ is assumed to be proportional to
335the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module
336\footnote{Following \citet{Li_Garrett_JMR93}, the surface Stoke drift velocity may be expressed as
337  $u_s =  0.016 \,|U_{10m}|$.
338  Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
339  $1.5~10^{-3}$ give the expression used of $u_s$ as a function of the module of surface stress
340}.
341For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at
342a finite depth $H_{LC}$ (which is often close to the mixed layer depth),
343and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures).
344The resulting expression for $w_{LC}$ is :
345\[
346  w_{LC}  =
347  \begin{cases}
348    c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$}    \\
349    0                             &      \text{otherwise}
350  \end{cases}
351\]
352where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise to fit LES data.
353The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second.
354The value of $c_{LC}$ is set through the \np{rn\_lc} namelist parameter,
355having in mind that it should stay between 0.15 and 0.54 \citep{Axell_JGR02}.
356
357The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations:
358$H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by
359converting its kinetic energy to potential energy, according to
360\[
361- \int_{-H_{LC}}^0 { N^2\;\;dz} = \frac{1}{2} u_s^2
362\]
363
364\subsubsection{Mixing just below the mixed layer}
365%--------------------------------------------------------------%
366
367Vertical mixing parameterizations commonly used in ocean general circulation models tend to
368produce mixed-layer depths that are too shallow during summer months and windy conditions.
369This bias is particularly acute over the Southern Ocean.
370To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme \cite{Rodgers_2014}.
371The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations,
372but rather is meant to account for observed processes that affect the density structure of
373the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme
374($i.e.$ near-inertial oscillations and ocean swells and waves).
375
376When using this parameterization ($i.e.$ when \np{nn\_etau}\forcode{ = 1}),
377the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations,
378swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition,
379plus a depth depend one given by:
380\begin{equation}
381  \label{eq:ZDF_Ehtau}
382  S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau}
383\end{equation}
384where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that
385penetrate in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of
386the penetration, and $f_i$ is the ice concentration
387(no penetration if $f_i=1$, that is if the ocean is entirely covered by sea-ice).
388The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter.
389The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{ = 0}) or
390a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m at high latitudes
391(\np{nn\_etau}\forcode{ = 1}).
392
393Note that two other option existe, \np{nn\_etau}\forcode{ = 2..3}.
394They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer,
395or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrate the ocean.
396Those two options are obsolescent features introduced for test purposes.
397They will be removed in the next release.
398
399% from Burchard et al OM 2008 :
400% the most critical process not reproduced by statistical turbulence models is the activity of
401% internal waves and their interaction with turbulence. After the Reynolds decomposition,
402% internal waves are in principle included in the RANS equations, but later partially
403% excluded by the hydrostatic assumption and the model resolution.
404% Thus far, the representation of internal wave mixing in ocean models has been relatively crude
405% (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
406
407% -------------------------------------------------------------------------------------------------------------
408%        TKE discretization considerations
409% -------------------------------------------------------------------------------------------------------------
410\subsection{TKE discretization considerations (\protect\key{zdftke})}
411\label{subsec:ZDF_tke_ene}
412
413%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
414\begin{figure}[!t]
415  \begin{center}
416    \includegraphics[width=1.00\textwidth]{Fig_ZDF_TKE_time_scheme}
417    \caption{
418      \protect\label{fig:TKE_time_scheme}
419      Illustration of the TKE time integration and its links to the momentum and tracer time integration.
420    }
421  \end{center} 
422\end{figure}
423%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
424
425The production of turbulence by vertical shear (the first term of the right hand side of
426\autoref{eq:zdftke_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion
427(first line in \autoref{eq:PE_zdf}).
428To do so a special care have to be taken for both the time and space discretization of
429the TKE equation \citep{Burchard_OM02,Marsaleix_al_OM08}.
430
431Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how
432the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with
433the one-level forward time stepping of TKE equation.
434With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to
435the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and
436summing the result vertically:   
437\begin{equation}
438  \label{eq:energ1}
439  \begin{split}
440    \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
441    &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}
442    - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
443  \end{split}
444\end{equation}
445Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$,
446known at time $t$ (\autoref{fig:TKE_time_scheme}), as it is required when using the TKE scheme
447(see \autoref{sec:STP_forward_imp}).
448The first term of the right hand side of \autoref{eq:energ1} represents the kinetic energy transfer at
449the surface (atmospheric forcing) and at the bottom (friction effect).
450The second term is always negative.
451It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$.
452\autoref{eq:energ1} implies that, to be energetically consistent,
453the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as
454${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$
455(and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).
456
457A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification
458(second term of the right hand side of \autoref{eq:zdftke_e}).
459This term must balance the input of potential energy resulting from vertical mixing.
460The rate of change of potential energy (in 1D for the demonstration) due vertical mixing is obtained by
461multiplying vertical density diffusion tendency by $g\,z$ and and summing the result vertically:
462\begin{equation}
463  \label{eq:energ2}
464  \begin{split}
465    \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
466    &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta}
467    - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
468    &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
469    + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
470  \end{split}
471\end{equation}
472where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.
473The first term of the right hand side of \autoref{eq:energ2} is always zero because
474there is no diffusive flux through the ocean surface and bottom).
475The second term is minus the destruction rate of  $\bar{e}$ due to stratification.
476Therefore \autoref{eq:energ1} implies that, to be energetically consistent,
477the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e}, the TKE equation.
478
479Let us now address the space discretization issue.
480The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in
481the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:cell}).
482A space averaging is thus required to obtain the shear TKE production term.
483By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of eddy coefficient by
484the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
485Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into account.
486
487The above energetic considerations leads to the following final discrete form for the TKE equation:
488\begin{equation}
489  \label{eq:zdftke_ene}
490  \begin{split}
491    \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv
492    \Biggl\{ \Biggr.
493    &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} }
494        \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
495    +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} }
496        \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j}
497    \Biggr. \Biggr\}   \\
498    %
499    - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
500    %
501    +&\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]   \\
502    %
503    - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
504  \end{split}
505\end{equation}
506where the last two terms in \autoref{eq:zdftke_ene} (vertical diffusion and Kolmogorov dissipation)
507are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}).
508Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible.
509The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as
510they all appear in the right hand side of \autoref{eq:zdftke_ene}.
511For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.
512
513% -------------------------------------------------------------------------------------------------------------
514%        GLS Generic Length Scale Scheme
515% -------------------------------------------------------------------------------------------------------------
516\subsection{GLS: Generic Length Scale (\protect\key{zdfgls})}
517\label{subsec:ZDF_gls}
518
519%--------------------------------------------namzdf_gls---------------------------------------------------------
520
521\nlst{namzdf_gls}
522%--------------------------------------------------------------------------------------------------------------
523
524The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations:
525one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale,
526$\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}.
527This later variable is defined as: $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$,
528where the triplet $(p, m, n)$ value given in Tab.\autoref{tab:GLS} allows to recover a number of
529well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982}, $k$-$\epsilon$ \citep{Rodi_1987},
530$k$-$\omega$ \citep{Wilcox_1988} among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}).
531The GLS scheme is given by the following set of equations:
532\begin{equation}
533  \label{eq:zdfgls_e}
534  \frac{\partial \bar{e}}{\partial t} =
535  \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
536      +\left( \frac{\partial v}{\partial k} \right)^2} \right]
537  -K_\rho \,N^2
538  +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
539  - \epsilon
540\end{equation}
541
542\[
543  % \label{eq:zdfgls_psi}
544  \begin{split}
545    \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
546      \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
547          +\left( \frac{\partial v}{\partial k} \right)^2} \right]
548      - C_3 \,K_\rho\,N^2   - C_2 \,\epsilon \,Fw   \right\}             \\
549    &+\frac{1}{e_3\;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 }
550        \;\frac{\partial \psi}{\partial k}} \right]\;
551  \end{split}
552\]
553
554\[
555  % \label{eq:zdfgls_kz}
556  \begin{split}
557    K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
558    K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l
559  \end{split}
560\]
561
562\[
563  % \label{eq:zdfgls_eps}
564  {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
565\]
566where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and
567$\epsilon$ the dissipation rate.
568The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of
569the choice of the turbulence model.
570Four different turbulent models are pre-defined (Tab.\autoref{tab:GLS}).
571They are made available through the \np{nn\_clo} namelist parameter.
572
573%--------------------------------------------------TABLE--------------------------------------------------
574\begin{table}[htbp]
575  \begin{center}
576    % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
577    \begin{tabular}{ccccc}
578      &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\
579      % & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\
580      \hline
581      \hline
582      \np{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\
583      \hline
584      $( p , n , m )$          &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
585      $\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
586      $\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
587      $C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
588      $C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
589      $C_3$              &      1.           &     1.              &      1.                &       1.           \\
590      $F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
591      \hline
592      \hline
593    \end{tabular}
594    \caption{
595      \protect\label{tab:GLS}
596      Set of predefined GLS parameters, or equivalently predefined turbulence models available with
597      \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}.
598    }
599  \end{center}
600\end{table}
601%--------------------------------------------------------------------------------------------------------------
602
603In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of
604the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) value near physical boundaries
605(logarithmic boundary layer law).
606$C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{Galperin_al_JAS88},
607or by \citet{Kantha_Clayson_1994} or one of the two functions suggested by \citet{Canuto_2001}
608(\np{nn\_stab\_func}\forcode{ = 0..3}, resp.).
609The value of $C_{0\mu}$ depends of the choice of the stability function.
610
611The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or
612Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp.
613As for TKE closure, the wave effect on the mixing is considered when
614\np{ln\_crban}\forcode{ = .true.} \citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}.
615The \np{rn\_crban} namelist parameter is $\alpha_{CB}$ in \autoref{eq:ZDF_Esbc} and
616\np{rn\_charn} provides the value of $\beta$ in \autoref{eq:ZDF_Lsbc}.
617
618The $\psi$ equation is known to fail in stably stratified flows, and for this reason
619almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy.
620With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$.
621A value of $c_{lim} = 0.53$ is often used \citep{Galperin_al_JAS88}.
622\cite{Umlauf_Burchard_CSR05} show that the value of the clipping factor is of crucial importance for
623the entrainment depth predicted in stably stratified situations,
624and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes.
625The clipping is only activated if \np{ln\_length\_lim}\forcode{ = .true.},
626and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
627
628The time and space discretization of the GLS equations follows the same energetic consideration as for
629the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{Burchard_OM02}.
630Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}.
631
632% -------------------------------------------------------------------------------------------------------------
633%        OSM OSMOSIS BL Scheme
634% -------------------------------------------------------------------------------------------------------------
635\subsection{OSM: OSMOSIS boundary layer scheme (\protect\key{zdfosm})}
636\label{subsec:ZDF_osm}
637
638%--------------------------------------------namzdf_osm---------------------------------------------------------
639
640\nlst{namzdf_osm}
641%--------------------------------------------------------------------------------------------------------------
642
643The OSMOSIS turbulent closure scheme is based on......   TBC
644
645% ================================================================
646% Convection
647% ================================================================
648\section{Convection}
649\label{sec:ZDF_conv}
650
651%--------------------------------------------namzdf--------------------------------------------------------
652
653\nlst{namzdf}
654%--------------------------------------------------------------------------------------------------------------
655
656Static instabilities (i.e. light potential densities under heavy ones) may occur at particular ocean grid points.
657In nature, convective processes quickly re-establish the static stability of the water column.
658These processes have been removed from the model via the hydrostatic assumption so they must be parameterized.
659Three parameterisations are available to deal with convective processes:
660a non-penetrative convective adjustment or an enhanced vertical diffusion,
661or/and the use of a turbulent closure scheme.
662
663% -------------------------------------------------------------------------------------------------------------
664%       Non-Penetrative Convective Adjustment
665% -------------------------------------------------------------------------------------------------------------
666\subsection[Non-penetrative convective adjmt (\protect\np{ln\_tranpc}\forcode{ = .true.})]
667            {Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})}
668\label{subsec:ZDF_npc}
669
670%--------------------------------------------namzdf--------------------------------------------------------
671
672\nlst{namzdf}
673%--------------------------------------------------------------------------------------------------------------
674
675%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
676\begin{figure}[!htb]
677  \begin{center}
678    \includegraphics[width=0.90\textwidth]{Fig_npc}
679    \caption{
680      \protect\label{fig:npc}
681      Example of an unstable density profile treated by the non penetrative convective adjustment algorithm.
682      $1^{st}$ step: the initial profile is checked from the surface to the bottom.
683      It is found to be unstable between levels 3 and 4.
684      They are mixed.
685      The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed.
686      The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed.
687      The $1^{st}$ step ends since the density profile is then stable below the level 3.
688      $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step:
689      levels 2 to 5 are mixed.
690      The new density profile is checked.
691      It is found stable: end of algorithm.
692    }
693  \end{center}
694\end{figure}
695%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
696
697Options are defined through the \ngn{namzdf} namelist variables.
698The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}\forcode{ = .true.}.
699It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of
700the water column, but only until the density structure becomes neutrally stable
701($i.e.$ until the mixed portion of the water column has \textit{exactly} the density of the water just below)
702\citep{Madec_al_JPO91}.
703The associated algorithm is an iterative process used in the following way (\autoref{fig:npc}):
704starting from the top of the ocean, the first instability is found.
705Assume in the following that the instability is located between levels $k$ and $k+1$.
706The temperature and salinity in the two levels are vertically mixed, conserving the heat and salt contents of
707the water column.
708The new density is then computed by a linear approximation.
709If the new density profile is still unstable between levels $k+1$ and $k+2$,
710levels $k$, $k+1$ and $k+2$ are then mixed.
711This process is repeated until stability is established below the level $k$
712(the mixing process can go down to the ocean bottom).
713The algorithm is repeated to check if the density profile between level $k-1$ and $k$ is unstable and/or
714if there is no deeper instability.
715
716This algorithm is significantly different from mixing statically unstable levels two by two.
717The latter procedure cannot converge with a finite number of iterations for some vertical profiles while
718the algorithm used in \NEMO converges for any profile in a number of iterations which is less than
719the number of vertical levels.
720This property is of paramount importance as pointed out by \citet{Killworth1989}:
721it avoids the existence of permanent and unrealistic static instabilities at the sea surface.
722This non-penetrative convective algorithm has been proved successful in studies of the deep water formation in
723the north-western Mediterranean Sea \citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
724
725The current implementation has been modified in order to deal with any non linear equation of seawater
726(L. Brodeau, personnal communication).
727Two main differences have been introduced compared to the original algorithm:
728$(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency
729(not the the difference in potential density);
730$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in
731the same way their temperature and salinity has been mixed.
732These two modifications allow the algorithm to perform properly and accurately with TEOS10 or EOS-80 without
733having to recompute the expansion coefficients at each mixing iteration.
734
735% -------------------------------------------------------------------------------------------------------------
736%       Enhanced Vertical Diffusion
737% -------------------------------------------------------------------------------------------------------------
738\subsection{Enhanced vertical diffusion (\protect\np{ln\_zdfevd}\forcode{ = .true.})}
739\label{subsec:ZDF_evd}
740
741%--------------------------------------------namzdf--------------------------------------------------------
742
743\nlst{namzdf}
744%--------------------------------------------------------------------------------------------------------------
745
746Options are defined through the  \ngn{namzdf} namelist variables.
747The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{ = .true.}.
748In this case, the vertical eddy mixing coefficients are assigned very large values
749(a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable
750($i.e.$ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar_PhD97, Lazar_al_JPO99}.
751This is done either on tracers only (\np{nn\_evdm}\forcode{ = 0}) or
752on both momentum and tracers (\np{nn\_evdm}\forcode{ = 1}).
753
754In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np{nn\_evdm}\forcode{ = 1},
755the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to
756the namelist parameter \np{rn\_avevd}.
757A typical value for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$.
758This parameterisation of convective processes is less time consuming than
759the convective adjustment algorithm presented above when mixing both tracers and
760momentum in the case of static instabilities.
761It requires the use of an implicit time stepping on vertical diffusion terms
762(i.e. \np{ln\_zdfexp}\forcode{ = .false.}).
763
764Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$.
765This removes a potential source of divergence of odd and even time step in
766a leapfrog environment \citep{Leclair_PhD2010} (see \autoref{sec:STP_mLF}).
767
768% -------------------------------------------------------------------------------------------------------------
769%       Turbulent Closure Scheme
770% -------------------------------------------------------------------------------------------------------------
771\subsection[Turbulent closure scheme (\protect\key{zdf}\{tke,gls,osm\})]{Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})}
772\label{subsec:ZDF_tcs}
773
774The turbulent closure scheme presented in \autoref{subsec:ZDF_tke} and \autoref{subsec:ZDF_gls}
775(\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically unstable density profiles.
776In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in
777\autoref{eq:zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative.
778It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring $A_u^{vm} {and}\;A_v^{vm}$
779(up to $1\;m^2s^{-1}$).
780These large values restore the static stability of the water column in a way similar to that of
781the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}).
782However, in the vicinity of the sea surface (first ocean layer), the eddy coefficients computed by
783the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$,
784because the mixing length scale is bounded by the distance to the sea surface.
785It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme,
786$i.e.$ setting the \np{ln\_zdfnpc} namelist parameter to true and
787defining the turbulent closure CPP key all together.
788
789The KPP turbulent closure scheme already includes enhanced vertical diffusion in the case of convection,
790as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},
791therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the KPP scheme.
792% gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
793
794% ================================================================
795% Double Diffusion Mixing
796% ================================================================
797\section{Double diffusion mixing (\protect\key{zdfddm})}
798\label{sec:ZDF_ddm}
799
800%-------------------------------------------namzdf_ddm-------------------------------------------------
801%
802%\nlst{namzdf_ddm}
803%--------------------------------------------------------------------------------------------------------------
804
805Options are defined through the  \ngn{namzdf\_ddm} namelist variables.
806Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa.
807The former condition leads to salt fingering and the latter to diffusive convection.
808Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean.
809\citet{Merryfield1999} include a parameterisation of such phenomena in a global ocean model and show that
810it leads to relatively minor changes in circulation but exerts significant regional influences on
811temperature and salinity.
812This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the \key{zdfddm} CPP key.
813
814Diapycnal mixing of S and T are described by diapycnal diffusion coefficients
815\begin{align*}
816  % \label{eq:zdfddm_Kz}
817  &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\
818  &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
819\end{align*}
820where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,
821and $o$ by processes other than double diffusion.
822The rates of double-diffusive mixing depend on the buoyancy ratio
823$R_\rho = \alpha \partial_z T / \beta \partial_z S$, where $\alpha$ and $\beta$ are coefficients of
824thermal expansion and saline contraction (see \autoref{subsec:TRA_eos}).
825To represent mixing of $S$ and $T$ by salt fingering, we adopt the diapycnal diffusivities suggested by Schmitt
826(1981):
827\begin{align}
828  \label{eq:zdfddm_f}
829  A_f^{vS} &=
830             \begin{cases}
831               \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
832               0                              &\text{otherwise}
833             \end{cases}
834  \\         \label{eq:zdfddm_f_T}
835  A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
836\end{align}
837
838%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
839\begin{figure}[!t]
840  \begin{center}
841    \includegraphics[width=0.99\textwidth]{Fig_zdfddm}
842    \caption{
843      \protect\label{fig:zdfddm}
844      From \citet{Merryfield1999} :
845      (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering.
846      Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$;
847      (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in regions of
848      diffusive convection.
849      Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation.
850      The latter is not implemented in \NEMO.
851    }
852  \end{center}
853\end{figure}
854%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
855
856The factor 0.7 in \autoref{eq:zdfddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx  0.7$ of
857buoyancy flux of heat to buoyancy flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}).
858Following  \citet{Merryfield1999}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
859
860To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by
861Federov (1988) is used:
862\begin{align}
863  % \label{eq:zdfddm_d}
864  A_d^{vT} &=
865             \begin{cases}
866               1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
867               &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
868               0                       &\text{otherwise}
869             \end{cases}
870                                       \nonumber \\
871  \label{eq:zdfddm_d_S}
872  A_d^{vS} &=
873             \begin{cases}
874               A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
875               A_d^{vT} \ 0.15 \ R_\rho               &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
876               0                       &\text{otherwise}
877             \end{cases}
878\end{align}
879
880The dependencies of \autoref{eq:zdfddm_f} to \autoref{eq:zdfddm_d_S} on $R_\rho$ are illustrated in
881\autoref{fig:zdfddm}.
882Implementing this requires computing $R_\rho$ at each grid point on every time step.
883This is done in \mdl{eosbn2} at the same time as $N^2$ is computed.
884This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive).
885
886% ================================================================
887% Bottom Friction
888% ================================================================
889\section{Bottom and top friction (\protect\mdl{zdfbfr})}
890\label{sec:ZDF_bfr}
891
892%--------------------------------------------nambfr--------------------------------------------------------
893%
894%\nlst{nambfr}
895%--------------------------------------------------------------------------------------------------------------
896
897Options to define the top and bottom friction are defined through the \ngn{nambfr} namelist variables.
898The bottom friction represents the friction generated by the bathymetry.
899The top friction represents the friction generated by the ice shelf/ocean interface.
900As the friction processes at the top and bottom are treated in similar way,
901only the bottom friction is described in detail below.
902
903
904Both the surface momentum flux (wind stress) and the bottom momentum flux (bottom friction) enter the equations as
905a condition on the vertical diffusive flux.
906For the bottom boundary layer, one has:
907\[
908  % \label{eq:zdfbfr_flux}
909  A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
910\]
911where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside
912the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean).
913How ${\cal F}_h^{\textbf U}$ influences the interior depends on the vertical resolution of the model near
914the bottom relative to the Ekman layer depth.
915For example, in order to obtain an Ekman layer depth $d = \sqrt{2\;A^{vm}} / f = 50$~m,
916one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$
917(for a Coriolis frequency $f = 10^{-4}$~m$^2$s$^{-1}$).
918With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m.
919When the vertical mixing coefficient is this small, using a flux condition is equivalent to
920entering the viscous forces (either wind stress or bottom friction) as a body force over the depth of the top or
921bottom model layer.
922To illustrate this, consider the equation for $u$ at $k$, the last ocean level:
923\begin{equation}
924  \label{eq:zdfbfr_flux2}
925  \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}}
926\end{equation}
927If the bottom layer thickness is 200~m, the Ekman transport will be distributed over that depth.
928On the other hand, if the vertical resolution is high (1~m or less) and a turbulent closure model is used,
929the turbulent Ekman layer will be represented explicitly by the model.
930However, the logarithmic layer is never represented in current primitive equation model applications:
931it is \emph{necessary} to parameterize the flux ${\cal F}^u_h $.
932Two choices are available in \NEMO: a linear and a quadratic bottom friction.
933Note that in both cases, the rotation between the interior velocity and the bottom friction is neglected in
934the present release of \NEMO.
935
936In the code, the bottom friction is imposed by adding the trend due to the bottom friction to
937the general momentum trend in \mdl{dynbfr}.
938For the time-split surface pressure gradient algorithm, the momentum trend due to
939the barotropic component needs to be handled separately.
940For this purpose it is convenient to compute and store coefficients which can be simply combined with
941bottom velocities and geometric values to provide the momentum trend due to bottom friction.
942These coefficients are computed in \mdl{zdfbfr} and generally take the form $c_b^{\textbf U}$ where:
943\begin{equation}
944  \label{eq:zdfbfr_bdef}
945  \frac{\partial {\textbf U_h}}{\partial t} =
946  - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
947\end{equation}
948where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
949
950% -------------------------------------------------------------------------------------------------------------
951%       Linear Bottom Friction
952% -------------------------------------------------------------------------------------------------------------
953\subsection{Linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 0..1})}
954\label{subsec:ZDF_bfr_linear}
955
956The linear bottom friction parameterisation (including the special case of a free-slip condition) assumes that
957the bottom friction is proportional to the interior velocity (i.e. the velocity of the last model level):
958\[
959  % \label{eq:zdfbfr_linear}
960  {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
961\]
962where $r$ is a friction coefficient expressed in ms$^{-1}$.
963This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean,
964and setting $r = H / \tau$, where $H$ is the ocean depth.
965Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}.
966A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used in quasi-geostrophic models.
967One may consider the linear friction as an approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$
968(\citet{Gill1982}, Eq. 9.6.6).
969For example, with a drag coefficient $C_D = 0.002$, a typical speed of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$,
970and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.
971This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days.
972It can be changed by specifying \np{rn\_bfri1} (namelist parameter).
973
974For the linear friction case the coefficients defined in the general expression \autoref{eq:zdfbfr_bdef} are:
975\[
976  % \label{eq:zdfbfr_linbfr_b}
977  \begin{split}
978    c_b^u &= - r\\
979    c_b^v &= - r\\
980  \end{split}
981\]
982When \np{nn\_botfr}\forcode{ = 1}, the value of $r$ used is \np{rn\_bfri1}.
983Setting \np{nn\_botfr}\forcode{ = 0} is equivalent to setting $r=0$ and
984leads to a free-slip bottom boundary condition.
985These values are assigned in \mdl{zdfbfr}.
986From v3.2 onwards there is support for local enhancement of these values via an externally defined 2D mask array
987(\np{ln\_bfr2d}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file.
988The mask values should vary from 0 to 1.
989Locations with a non-zero mask value will have the friction coefficient increased by
990$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}.
991
992% -------------------------------------------------------------------------------------------------------------
993%       Non-Linear Bottom Friction
994% -------------------------------------------------------------------------------------------------------------
995\subsection{Non-linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 2})}
996\label{subsec:ZDF_bfr_nonlinear}
997
998The non-linear bottom friction parameterisation assumes that the bottom friction is quadratic:
999\[
1000  % \label{eq:zdfbfr_nonlinear}
1001  {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h
1002  }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b
1003\]
1004where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy due to tides,
1005internal waves breaking and other short time scale currents.
1006A typical value of the drag coefficient is $C_D = 10^{-3} $.
1007As an example, the CME experiment \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and
1008$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} uses $C_D = 1.4\;10^{-3}$ and
1009$e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.
1010The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} namelist parameters).
1011
1012As for the linear case, the bottom friction is imposed in the code by adding the trend due to
1013the bottom friction to the general momentum trend in \mdl{dynbfr}.
1014For the non-linear friction case the terms computed in \mdl{zdfbfr} are:
1015\[
1016  % \label{eq:zdfbfr_nonlinbfr}
1017  \begin{split}
1018    c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\
1019    c_b^v &= - \; C_D\;\left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
1020  \end{split}
1021\]
1022
1023The coefficients that control the strength of the non-linear bottom friction are initialised as namelist parameters:
1024$C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.
1025Note for applications which treat tides explicitly a low or even zero value of \np{rn\_bfeb2} is recommended.
1026From v3.2 onwards a local enhancement of $C_D$ is possible via an externally defined 2D mask array
1027(\np{ln\_bfr2d}\forcode{ = .true.}).
1028This works in the same way as for the linear bottom friction case with non-zero masked locations increased by
1029$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}.
1030
1031% -------------------------------------------------------------------------------------------------------------
1032%       Bottom Friction Log-layer
1033% -------------------------------------------------------------------------------------------------------------
1034\subsection[Log-layer btm frict enhncmnt (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})]
1035            {Log-layer bottom friction enhancement (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})}
1036\label{subsec:ZDF_bfr_loglayer}
1037
1038In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally enhanced using
1039a "law of the wall" scaling.
1040If  \np{ln\_loglayer} = .true., $C_D$ is no longer constant but is related to the thickness of
1041the last wet layer in each column by:
1042\[
1043  C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2
1044\]
1045
1046\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness length provided via
1047the namelist.
1048
1049For stability, the drag coefficient is bounded such that it is kept greater or equal to
1050the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional namelist parameter:
1051\np{rn\_bfri2\_max}, i.e.:
1052\[
1053  rn\_bfri2 \leq C_D \leq rn\_bfri2\_max
1054\]
1055
1056\noindent Note also that a log-layer enhancement can also be applied to the top boundary friction if
1057under ice-shelf cavities are in use (\np{ln\_isfcav}\forcode{ = .true.}).
1058In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}.
1059
1060% -------------------------------------------------------------------------------------------------------------
1061%       Bottom Friction stability
1062% -------------------------------------------------------------------------------------------------------------
1063\subsection{Bottom friction stability considerations}
1064\label{subsec:ZDF_bfr_stability}
1065
1066Some care needs to exercised over the choice of parameters to ensure that the implementation of
1067bottom friction does not induce numerical instability.
1068For the purposes of stability analysis, an approximation to \autoref{eq:zdfbfr_flux2} is:
1069\begin{equation}
1070  \label{eq:Eqn_bfrstab}
1071  \begin{split}
1072    \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
1073    &= -\frac{ru}{e_{3u}}\;2\rdt\\
1074  \end{split}
1075\end{equation}
1076\noindent where linear bottom friction and a leapfrog timestep have been assumed.
1077To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have:
1078\[
1079  |\Delta u| < \;|u|
1080\]
1081\noindent which, using \autoref{eq:Eqn_bfrstab}, gives:
1082\[
1083  r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
1084\]
1085This same inequality can also be derived in the non-linear bottom friction case if
1086a velocity of 1 m.s$^{-1}$ is assumed.
1087Alternatively, this criterion can be rearranged to suggest a minimum bottom box thickness to ensure stability:
1088\[
1089  e_{3u} > 2\;r\;\rdt
1090\]
1091\noindent which it may be necessary to impose if partial steps are being used.
1092For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then $e_{3u}$ should be greater than 3.6 m.
1093For most applications, with physically sensible parameters these restrictions should not be of concern.
1094But caution may be necessary if attempts are made to locally enhance the bottom friction parameters.
1095To ensure stability limits are imposed on the bottom friction coefficients both
1096during initialisation and at each time step.
1097Checks at initialisation are made in \mdl{zdfbfr} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
1098The number of breaches of the stability criterion are reported as well as
1099the minimum and maximum values that have been set.
1100The criterion is also checked at each time step, using the actual velocity, in \mdl{dynbfr}.
1101Values of the bottom friction coefficient are reduced as necessary to ensure stability;
1102these changes are not reported.
1103
1104Limits on the bottom friction coefficient are not imposed if the user has elected to
1105handle the bottom friction implicitly (see \autoref{subsec:ZDF_bfr_imp}).
1106The number of potential breaches of the explicit stability criterion are still reported for information purposes.
1107
1108% -------------------------------------------------------------------------------------------------------------
1109%       Implicit Bottom Friction
1110% -------------------------------------------------------------------------------------------------------------
1111\subsection{Implicit bottom friction (\protect\np{ln\_bfrimp}\forcode{ = .true.})}
1112\label{subsec:ZDF_bfr_imp}
1113
1114An optional implicit form of bottom friction has been implemented to improve model stability.
1115We recommend this option for shelf sea and coastal ocean applications, especially for split-explicit time splitting.
1116This option can be invoked by setting \np{ln\_bfrimp} to \forcode{.true.} in the \textit{nambfr} namelist.
1117This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist.
1118
1119This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp},
1120the bottom boundary condition is implemented implicitly.
1121
1122\[
1123  % \label{eq:dynzdf_bfr}
1124  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
1125  = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
1126\]
1127
1128where $mbk$ is the layer number of the bottom wet layer.
1129Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so, it is implicit.
1130
1131If split-explicit time splitting is used, care must be taken to avoid the double counting of the bottom friction in
1132the 2-D barotropic momentum equations.
1133As NEMO only updates the barotropic pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation,
1134we need to remove the bottom friction induced by these two terms which has been included in the 3-D momentum trend
1135and update it with the latest value.
1136On the other hand, the bottom friction contributed by the other terms
1137(e.g. the advection term, viscosity term) has been included in the 3-D momentum equations and
1138should not be added in the 2-D barotropic mode.
1139
1140The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the following:
1141
1142\[
1143  % \label{eq:dynspg_ts_bfr1}
1144  \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
1145  \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
1146\]
1147\[
1148  \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
1149  \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
1150  2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
1151\]
1152
1153where $\textbf{T}$ is the vertical integrated 3-D momentum trend.
1154We assume the leap-frog time-stepping is used here.
1155$\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step.
1156$c_{b}$ is the friction coefficient.
1157$\eta$ is the sea surface level calculated in the barotropic loops while $\eta^{'}$ is the sea surface level used in
1158the 3-D baroclinic mode.
1159$\textbf{u}_{b}$ is the bottom layer horizontal velocity.
1160
1161% -------------------------------------------------------------------------------------------------------------
1162%       Bottom Friction with split-explicit time splitting
1163% -------------------------------------------------------------------------------------------------------------
1164\subsection[Bottom friction w/ split-explicit time splitting (\protect\np{ln\_bfrimp})]
1165            {Bottom friction with split-explicit time splitting (\protect\np{ln\_bfrimp})}
1166\label{subsec:ZDF_bfr_ts}
1167
1168When calculating the momentum trend due to bottom friction in \mdl{dynbfr},
1169the bottom velocity at the before time step is used.
1170This velocity includes both the baroclinic and barotropic components which is appropriate when
1171using either the explicit or filtered surface pressure gradient algorithms
1172(\key{dynspg\_exp} or \key{dynspg\_flt}).
1173Extra attention is required, however, when using split-explicit time stepping (\key{dynspg\_ts}).
1174In this case the free surface equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro},
1175while the three dimensional prognostic variables are solved with the longer time step of \np{rn\_rdt} seconds.
1176The trend in the barotropic momentum due to bottom friction appropriate to this method is that given by
1177the selected parameterisation ($i.e.$ linear or non-linear bottom friction) computed with
1178the evolving velocities at each barotropic timestep.
1179
1180In the case of non-linear bottom friction, we have elected to partially linearise the problem by
1181keeping the coefficients fixed throughout the barotropic time-stepping to those computed in
1182\mdl{zdfbfr} using the now timestep.
1183This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to:
1184
1185\begin{enumerate}
1186\item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before barotropic velocity to
1187  the bottom friction component of the vertically integrated momentum trend.
1188  Note the same stability check that is carried out on the bottom friction coefficient in \mdl{dynbfr} has to
1189  be applied here to ensure that the trend removed matches that which was added in \mdl{dynbfr}.
1190\item At each barotropic step, compute the contribution of the current barotropic velocity to
1191  the trend due to bottom friction.
1192  Add this contribution to the vertically integrated momentum trend.
1193  This contribution is handled implicitly which eliminates the need to impose a stability criteria on
1194  the values of the bottom friction coefficient within the barotropic loop.
1195\end{enumerate}
1196
1197Note that the use of an implicit formulation within the barotropic loop for the bottom friction trend means that
1198any limiting of the bottom friction coefficient in \mdl{dynbfr} does not adversely affect the solution when
1199using split-explicit time splitting.
1200This is because the major contribution to bottom friction is likely to come from the barotropic component which
1201uses the unrestricted value of the coefficient.
1202However, if the limiting is thought to be having a major effect
1203(a more likely prospect in coastal and shelf seas applications) then
1204the fully implicit form of the bottom friction should be used (see \autoref{subsec:ZDF_bfr_imp})
1205which can be selected by setting \np{ln\_bfrimp} $=$ \forcode{.true.}.
1206
1207Otherwise, the implicit formulation takes the form:
1208\[
1209  % \label{eq:zdfbfr_implicitts}
1210  \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ]
1211\]
1212where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height),
1213$c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and
1214$RHS$ represents all the components to the vertically integrated momentum trend except for
1215that due to bottom friction.
1216
1217% ================================================================
1218% Tidal Mixing
1219% ================================================================
1220\section{Tidal mixing (\protect\key{zdftmx})}
1221\label{sec:ZDF_tmx}
1222
1223%--------------------------------------------namzdf_tmx--------------------------------------------------
1224%
1225%\nlst{namzdf_tmx}
1226%--------------------------------------------------------------------------------------------------------------
1227
1228
1229% -------------------------------------------------------------------------------------------------------------
1230%        Bottom intensified tidal mixing
1231% -------------------------------------------------------------------------------------------------------------
1232\subsection{Bottom intensified tidal mixing}
1233\label{subsec:ZDF_tmx_bottom}
1234
1235Options are defined through the  \ngn{namzdf\_tmx} namelist variables.
1236The parameterization of tidal mixing follows the general formulation for the vertical eddy diffusivity proposed by
1237\citet{St_Laurent_al_GRL02} and first introduced in an OGCM by \citep{Simmons_al_OM04}.
1238In this formulation an additional vertical diffusivity resulting from internal tide breaking,
1239$A^{vT}_{tides}$ is expressed as a function of $E(x,y)$,
1240the energy transfer from barotropic tides to baroclinic tides:
1241\begin{equation}
1242  \label{eq:Ktides}
1243  A^{vT}_{tides} =  q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 }
1244\end{equation}
1245where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
1246$\rho$ the density, $q$ the tidal dissipation efficiency, and $F(z)$ the vertical structure function.
1247
1248The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter) and
1249is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980).
1250The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter)
1251represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally,
1252with the remaining $1-q$ radiating away as low mode internal waves and
1253contributing to the background internal wave field.
1254A value of $q=1/3$ is typically used \citet{St_Laurent_al_GRL02}.
1255The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical.
1256It is implemented as a simple exponential decaying upward away from the bottom,
1257with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter,
1258with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04},
1259\[
1260  % \label{eq:Fz}
1261  F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
1262\]
1263and is normalized so that vertical integral over the water column is unity.
1264
1265The associated vertical viscosity is calculated from the vertical diffusivity assuming a Prandtl number of 1,
1266$i.e.$ $A^{vm}_{tides}=A^{vT}_{tides}$.
1267In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity is capped at $300\,cm^2/s$ and
1268impose a lower limit on $N^2$ of \np{rn\_n2min} usually set to $10^{-8} s^{-2}$.
1269These bounds are usually rarely encountered.
1270
1271The internal wave energy map, $E(x, y)$ in \autoref{eq:Ktides}, is derived from a barotropic model of
1272the tides utilizing a parameterization of the conversion of barotropic tidal energy into internal waves.
1273The essential goal of the parameterization is to represent the momentum exchange between the barotropic tides and
1274the unrepresented internal waves induced by the tidal flow over rough topography in a stratified ocean.
1275In the current version of \NEMO, the map is built from the output of
1276the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
1277This model provides the dissipation associated with internal wave energy for the M2 and K1 tides component
1278(\autoref{fig:ZDF_M2_K1_tmx}).
1279The S2 dissipation is simply approximated as being $1/4$ of the M2 one.
1280The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$.
1281Its global mean value is $1.1$ TW,
1282in agreement with independent estimates \citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}.
1283
1284%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1285\begin{figure}[!t]
1286  \begin{center}
1287    \includegraphics[width=0.90\textwidth]{Fig_ZDF_M2_K1_tmx}
1288    \caption{
1289      \protect\label{fig:ZDF_M2_K1_tmx}
1290      (a) M2 and (b) K1 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$).
1291    }
1292  \end{center}
1293\end{figure}
1294%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1295 
1296% -------------------------------------------------------------------------------------------------------------
1297%        Indonesian area specific treatment
1298% -------------------------------------------------------------------------------------------------------------
1299\subsection{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})}
1300\label{subsec:ZDF_tmx_itf}
1301
1302When the Indonesian Through Flow (ITF) area is included in the model domain,
1303a specific treatment of tidal induced mixing in this area can be used.
1304It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide an input NetCDF file,
1305\ifile{mask\_itf}, which contains a mask array defining the ITF area where the specific treatment is applied.
1306
1307When \np{ln\_tmx\_itf}\forcode{ = .true.}, the two key parameters $q$ and $F(z)$ are adjusted following
1308the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}:
1309
1310First, the Indonesian archipelago is a complex geographic region with a series of
1311large, deep, semi-enclosed basins connected via numerous narrow straits.
1312Once generated, internal tides remain confined within this semi-enclosed area and hardly radiate away.
1313Therefore all the internal tides energy is consumed within this area.
1314So it is assumed that $q = 1$, $i.e.$ all the energy generated is available for mixing.
1315Note that for test purposed, the ITF tidal dissipation efficiency is a namelist parameter (\np{rn\_tfe\_itf}).
1316A value of $1$ or close to is this recommended for this parameter.
1317
1318Second, the vertical structure function, $F(z)$, is no more associated with a bottom intensification of the mixing,
1319but with a maximum of energy available within the thermocline.
1320\citet{Koch-Larrouy_al_GRL07} have suggested that the vertical distribution of
1321the energy dissipation proportional to $N^2$ below the core of the thermocline and to $N$ above.
1322The resulting $F(z)$ is:
1323\[
1324  % \label{eq:Fz_itf}
1325  F(i,j,k) \sim     \left\{
1326    \begin{aligned}
1327      \frac{q\,\Gamma E(i,j) } {\rho N \, \int N     dz}    \qquad \text{when $\partial_z N < 0$} \\
1328      \frac{q\,\Gamma E(i,j) } {\rho     \, \int N^2 dz}    \qquad \text{when $\partial_z N > 0$}
1329    \end{aligned}
1330  \right.
1331\]
1332
1333Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$,
1334which agrees with the independent estimates inferred from observations.
1335Introduced in a regional OGCM, the parameterization improves the water mass characteristics in
1336the different Indonesian seas, suggesting that the horizontal and vertical distributions of
1337the mixing are adequately prescribed \citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}.
1338Note also that such a parameterisation has a significant impact on the behaviour of
1339global coupled GCMs \citep{Koch-Larrouy_al_CD10}.
1340
1341% ================================================================
1342% Internal wave-driven mixing
1343% ================================================================
1344\section{Internal wave-driven mixing (\protect\key{zdftmx\_new})}
1345\label{sec:ZDF_tmx_new}
1346
1347%--------------------------------------------namzdf_tmx_new------------------------------------------
1348%
1349%\nlst{namzdf_tmx_new}
1350%--------------------------------------------------------------------------------------------------------------
1351
1352The parameterization of mixing induced by breaking internal waves is a generalization of
1353the approach originally proposed by \citet{St_Laurent_al_GRL02}.
1354A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,
1355and the resulting diffusivity is obtained as
1356\[
1357  % \label{eq:Kwave}
1358  A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
1359\]
1360where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of
1361the energy available for mixing.
1362If the \np{ln\_mevar} namelist parameter is set to false, the mixing efficiency is taken as constant and
1363equal to 1/6 \citep{Osborn_JPO80}.
1364In the opposite (recommended) case, $R_f$ is instead a function of
1365the turbulence intensity parameter $Re_b = \frac{ \epsilon}{\nu \, N^2}$,
1366with $\nu$ the molecular viscosity of seawater, following the model of \cite{Bouffard_Boegman_DAO2013} and
1367the implementation of \cite{de_lavergne_JPO2016_efficiency}.
1368Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when
1369the mixing efficiency is constant.
1370
1371In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary
1372as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice.
1373This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014},
1374is implemented as in \cite{de_lavergne_JPO2016_efficiency}.
1375
1376The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$,
1377is constructed from three static maps of column-integrated internal wave energy dissipation,
1378$E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures
1379(de Lavergne et al., in prep):
1380\begin{align*}
1381  F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
1382  F_{pyc}(i,j,k) &\propto N^{n\_p}\\
1383  F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
1384\end{align*} 
1385In the above formula, $h_{ab}$ denotes the height above bottom,
1386$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
1387\[
1388  h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
1389\]
1390The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)
1391controls the stratification-dependence of the pycnocline-intensified dissipation.
1392It can take values of 1 (recommended) or 2.
1393Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of
1394the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.
1395$h_{cri}$ is related to the large-scale topography of the ocean (etopo2) and
1396$h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of
1397the abyssal hill topography \citep{Goff_JGR2010} and the latitude.
1398
1399% ================================================================
1400
1401\biblio
1402
1403\end{document}
Note: See TracBrowser for help on using the repository browser.