# source:trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex@2541

Last change on this file since 2541 was 2541, checked in by gm, 10 years ago

v3.3beta: #658 hopefully the final update, including MPP sum, SIGN, IOM, fldread documentation

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