source: trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex @ 6997

Last change on this file since 6997 was 6997, checked in by nicolasmartin, 5 years ago

Duplication of changes in DOC directory for the trunk

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