New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Chap_ZDF.tex in branches/DEV_r1826_DOC/DOC/TexFiles/Chapters – NEMO

source: branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_ZDF.tex @ 1831

Last change on this file since 1831 was 1831, checked in by gm, 14 years ago

cover, namelist, rigid-lid, e3t, appendices, see ticket: #658

  • Property svn:executable set to *
File size: 40.4 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{np\_zdfexp}=true) or a backward time stepping
43scheme (\np{np\_zdfexp}=false) depending on the magnitude of the mixing
44coefficients, and thus of the formulation used (see \S\ref{DOM_nxt}).
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}  \\
61\\
62A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
63\end{align*}
64
65These values are set through the \np{avm0} and \np{avt0} namelist parameters.
66In all cases, do not use values smaller that those associated with the molecular
67viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum,
68$\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity.
69
70
71% -------------------------------------------------------------------------------------------------------------
72%        Richardson Number Dependent
73% -------------------------------------------------------------------------------------------------------------
74\subsection{Richardson Number Dependent (\key{zdfric})}
75\label{ZDF_ric}
76
77%--------------------------------------------namric---------------------------------------------------------
78\namdisplay{namzdf_ric}
79%--------------------------------------------------------------------------------------------------------------
80
81When \key{zdfric} is defined, a local Richardson number dependent formulation
82for the vertical momentum and tracer eddy coefficients is set. The vertical mixing
83coefficients are diagnosed from the large scale variables computed by the model.
84\textit{In situ} measurements have been used to link vertical turbulent activity to
85large scale ocean structures. The hypothesis of a mixing mainly maintained by the
86growth of Kelvin-Helmholtz like instabilities leads to a dependency between the
87vertical eddy coefficients and the local Richardson number ($i.e.$ the
88ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following
89formulation has been implemented:
90\begin{equation} \label{Eq_zdfric}
91   \left\{      \begin{aligned}
92         A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
93         \\
94         A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
95   \end{aligned}    \right.
96\end{equation}
97where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson
98number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),
99$A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the
100constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ 
101is the maximum value that can be reached by the coefficient when $Ri\leq 0$,
102$a=5$ and $n=2$. The last three values can be modified by setting the
103\np{avmri}, \np{alp} and \np{nric} namelist parameters, respectively.
104
105% -------------------------------------------------------------------------------------------------------------
106%        TKE Turbulent Closure Scheme
107% -------------------------------------------------------------------------------------------------------------
108\subsection{TKE Turbulent Closure Scheme (\key{zdftke})}
109\label{ZDF_tke}
110
111%--------------------------------------------namzdf_tke--------------------------------------------------
112\namdisplay{namzdf_tke}
113%--------------------------------------------------------------------------------------------------------------
114
115The vertical eddy viscosity and diffusivity coefficients are computed from a TKE
116turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent
117kinetic energy, and a closure assumption for the turbulent length scales. This
118turbulent closure model has been developed by \citet{Bougeault1989} in the
119atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and
120embedded in OPA, the ancestor of NEMO, by \citet{Blanke1993} for equatorial Atlantic
121simulations. Since then, significant modifications have been introduced by
122\citet{Madec1998} in both the implementation and the formulation of the mixing
123length scale. The time evolution of $\bar{e}$ is the result of the production of
124$\bar{e}$ through vertical shear, its destruction through stratification, its vertical
125diffusion, and its dissipation of \citet{Kolmogorov1942} type:
126\begin{equation} \label{Eq_zdftke_e}
127\frac{\partial \bar{e}}{\partial t} =
128\frac{A^{vm}}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
129                     +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
130-A^{vT}\,N^2
131+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
132            \;\frac{\partial \bar{e}}{\partial k}} \right]
133- c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
134\end{equation}
135\begin{equation} \label{Eq_zdftke_kz}
136   \begin{split}
137         A^{vm} &= C_k\  l_k\  \sqrt {\bar{e}\; }     \\
138         A^{vT} &= A^{vm} / P_{rt}
139   \end{split}
140\end{equation}
141where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),
142$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales,
143$P_{rt}$ is the Prandtl number. The constants $C_k =  0.1$ and
144$C_\epsilon = \sqrt {2} /2$  $\approx 0.7$ are designed to deal with vertical mixing
145at any depth \citep{Gaspar1990}. They are set through namelist parameters
146\np{nn\_ediff} and \np{nn\_ediss}. $P_{rt}$ can be set to unity or, following
147\citet{Blanke1993}, be a function of the local Richardson number, $R_i$:
148\begin{align*} \label{Eq_prt}
149P_{rt} = \begin{cases}
150                    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}  \\
151                    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}  \\
152                    \ \ 10 &      \text{if $\ 2 \leq R_i$} 
153            \end{cases}
154\end{align*}
155The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist parameter.
156
157For computational efficiency, the original formulation of the turbulent length
158scales proposed by \citet{Gaspar1990} has been simplified. Four formulations
159are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist
160parameter. The first two are based on the following first order approximation
161\citep{Blanke1993}:
162\begin{equation} \label{Eq_tke_mxl0_1}
163l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
164\end{equation}
165which is valid in a stable stratified region with constant values of the Brunt-
166Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance
167to the surface or to the bottom (\np{nn\_mxl}=0) or by the local vertical scale factor
168(\np{nn\_mxl}=1). \citet{Blanke1993} notice that this simplification has two major
169drawbacks: it makes no sense for locally unstable stratification and the
170computation no longer uses all the information contained in the vertical density
171profile. To overcome these drawbacks, \citet{Madec1998} introduces the
172\np{nn\_mxl}=2 or 3 cases, which add an extra assumption concerning the vertical
173gradient of the computed length scale. So, the length scales are first evaluated
174as in \eqref{Eq_tke_mxl0_1} and then bounded such that:
175\begin{equation} \label{Eq_tke_mxl_constraint}
176\frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
177\qquad \text{with }\  l =  l_k = l_\epsilon
178\end{equation}
179\eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length
180scale cannot be larger than the variations of depth. It provides a better
181approximation of the \citet{Gaspar1990} formulation while being much less
182time consuming. In particular, it allows the length scale to be limited not only
183by the distance to the surface or to the ocean bottom but also by the distance
184to a strongly stratified portion of the water column such as the thermocline
185(Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} 
186constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$,
187the upward and downward length scales, and evaluate the dissipation and
188mixing length scales as (and note that here we use numerical indexing):
189%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
190\begin{figure}[!t] \label{Fig_mixing_length}  \begin{center}
191\includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_mixing_length.pdf}
192\caption {Illustration of the mixing length computation. }
193\end{center} 
194\end{figure}
195%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
196\begin{equation} \label{Eq_tke_mxl2}
197\begin{aligned}
198 l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
199    \quad &\text{ from $k=1$ to $jpk$ }\ \\
200 l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)}  \right)   
201    \quad &\text{ from $k=jpk$ to $1$ }\ \\
202\end{aligned}
203\end{equation}
204where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1},
205$i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
206
207In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same
208value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the
209\np{nmxl}=2 case, the dissipation and mixing turbulent length scales are give
210as in \citet{Gaspar1990}:
211\begin{equation} \label{Eq_tke_mxl_gaspar}
212\begin{aligned}
213& l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }    \\
214& l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)
215\end{aligned}
216\end{equation}
217
218At the sea surface the value of $\bar{e}$ is prescribed from the wind
219stress field: $\bar{e}=rn\_ebb\;\left| \tau \right|$ (\np{rn\_ebb}=60 by default)
220with a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist
221parameters). Its value at the bottom of the ocean is assumed to be
222equal to the value of the level just above. The time integration of the
223$\bar{e}$ equation may formally lead to negative values because the
224numerical scheme does not ensure its positivity. To overcome this
225problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 
226namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set
227to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations
228to match that of \citet{Gargett1984} for the diffusion in the thermocline and
229deep ocean :  $(A^{vT} = 10^{-3} / N)$.
230In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical
231instabilities associated with too weak vertical diffusion. They must be
232specified at least larger than the molecular values, and are set through
233\np{avm0} and \np{avt0} (namelist parameters).
234
235% -------------------------------------------------------------------------------------------------------------
236%        TKE Turbulent Closure Scheme : new organization to energetic considerations
237% -------------------------------------------------------------------------------------------------------------
238\subsection{TKE Turbulent Closure Scheme - time integration (\key{zdftke} and (\key{zdftke2})}
239\label{ZDF_tke2}
240
241%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
242\begin{figure}[!t] \label{Fig_TKE_time_scheme}  \begin{center}
243\includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_ZDF_TKE_time_scheme.pdf}
244\caption {Illustration of the TKE time integration and its links to the momentum and tracer time integration. }
245\end{center} 
246\end{figure}
247%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
248
249The production of turbulence by vertical shear (the first term of the right hand side
250of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with
251the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}).
252The total loss of kinetic energy (in 1D for the demonstration)
253due to this term is obtained by multiply this quantity by $u^n$ and verticaly integrating:   
254
255\begin{equation} \label{Eq_energ1}
256\int_{k_b}^{k_s} {u^t \frac{1}{e_3}
257                                 \frac{\partial   } 
258                                        {\partial k} \left( \frac{A^{vm}}{e_3} 
259                                                                 \frac{\partial{u^{t+1}}}
260                                                                        {\partial k           }     \right) \; e_3 \; dk }
261= \left[  \frac{u^t}{e_3}   A^{vm} \frac{\partial{u^{t+1}}}{\partial k} \right]_{k_b}^{k_s}
262 - \int_{k_b}^{k_s}{\frac{A^{vm}}{{e_3}}\frac{\partial{u^t}}{\partial k}\frac{\partial{u^{t+1}}}{\partial k}} \ dk
263\end{equation}
264
265The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic
266energy transfer at the surface (atmospheric forcing) and at the bottom (friction effect).
267The second term is always negative and have to balance the term of \eqref{Eq_zdftke_e} 
268previously identified.
269
270The sink term (possibly a source term in statically unstable situations) of turbulence
271by buoyancy (second term of the right hand side of \eqref{Eq_zdftke_e}) must balance
272the source of potential energy associated with the vertical diffusion
273in the density equation (second line in \eqref{Eq_PE_zdf}). The source of potential
274energy (in 1D for the demonstration) due to this term is obtained by multiply this quantity
275by $gz{\rho_r}^{-1}$ and verticaly integrating:
276
277\begin{equation} \label{Eq_energ2}
278\begin{aligned}
279\int_{k_b}^{k_s}{\frac{g\;z}{e_3} \frac{\partial }{\partial k}
280                        \left( \frac{A^{vT}}{e_3} 
281                        \frac{\partial{\rho^{t+1}}}{\partial k}   \right)} \; e_3 \; dk
282=\left[  g\;z \frac{A^{vT}}{e_3}
283                 \frac{\partial{\rho^{t+1}}}{\partial k} \right]_{k_b}^{k_s} 
284- \int_{k_b}^{k_s}{  \frac{A^{vT}}{e_3} g \frac{\partial{\rho^{t+1}}}{\partial k}} \; dk\\
285\\
286= - \left[  z\,A^{vT} {N_{t+1}}^2 \right]_{k_b}^{k_s}
287+ \int_{k_b}^{k_s}{  \rho^{t+1} \; A^{vT}{N_{t+1}}^2 \; e_3 \; dk  }\\
288\end{aligned}
289\end{equation}
290where $N^2_{t+1}$ is  the Brunt-Vaissala frequency at $t+1$ 
291and noting that $\frac{\partial z}{\partial k} = e_3$.
292
293The first term is always zero because the Brunt Vaissala frequency is set to zero at the
294surface and the bottom. The second term is of opposite sign than the buoyancy term
295identified previously.
296
297Under these energetic considerations, \eqref{Eq_zdftke_e} have to be written like this
298to be consistant:
299
300\begin{equation} \label{Eq_zdftke_ene}
301\frac{\partial \bar{e}}{\partial t} =
302\frac{A^{vm}}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u^a}{\partial k}} \right)\left( {\frac{\partial u^n}{\partial k}} \right)
303                                                        +\left( {\frac{\partial v^a}{\partial k}} \right)\left( {\frac{\partial v^n}{\partial k}} \right)
304} \right]-A^{vT}\,{N_n}^2+...
305\end{equation}
306
307Note that during a time step, the equation \eqref{Eq_zdftke_e} is resolved before those
308of momentum and density. So, the indice "a" (after) become "n" (now) and the indice "n"
309(now) become "b" (before).
310
311Moreover, the vertical shear have to be multiply by the appropriate viscosity for numerical
312stability. Thus, the vertical shear at U-point have to be multiply by the viscosity avmu and
313the vertical shear at V-point have to be multiply by the viscosity avmv. Next, these two
314quantities are averaged to obtain a production term by vertical shear at W-point :
315
316\begin{equation} \label{Eq_zdftke_ene2}
317\frac{\partial \bar{e}}{\partial t} =
318\frac{1}{{e_3}^2 }\;\left[ { 
319    A^{vmu}\left({\frac{\partial u^a}{\partial k}} \right)\left( {\frac{\partial u^n}{\partial k}} \right)
320   +A^{vmv}\left({\frac{\partial v^a}{\partial k}} \right)\left( {\frac{\partial v^n}{\partial k}} \right)
321                           } \right]-A^{vT}\,{N_n}^2+...
322\end{equation} 
323
324The TKE equation is resolved before the mixing length, the viscosity and diffusivity. Two tabs
325are then declared : dissl (dissipation length) (Remark : it's not only the dissipation lenght,
326it's the root of the TKE divided by the dissipation lenght) and avmt (viscosity at the points T)
327used for the vertical diffusion of the TKE.
328
329This new organization needs also a reorganization of the routine step.F90 (controled by
330the key \key{ene\_cons}). The bigger change is the estimation of the Brunt-Vaissala
331frequency at "n" instead of "b". Moreover for energetic considerations, the call of tranxt.F90
332is done after the density update. On the contrary, the density is updated with scalars fields
333filtered by the Asselin filter.
334
335This new organisation of the routine zdftke force to save five three dimensionnal tabs in
336the restart : avmu, avmv, avt, avmt and dissl are needed to calculate $e_n$. At the end
337of the run (time step = nitend), the alternative is to save only $e_n$ estimated at the
338following time step (nitend+1). The next run using this restart file, the mixing length
339and turbulents coefficients are directly calculated using $e_n$. It is the same thing
340for the intermediate restart.
341
342
343%%GM for figure of the time scheme:
344\begin{equation} 
345  \rho^{t+1} \; A^{vT}{N_{t+1}}^2 \; dk  \\
346\end{equation}
347
348%%end GM
349
350
351% -------------------------------------------------------------------------------------------------------------
352%        K Profile Parametrisation (KPP)
353% -------------------------------------------------------------------------------------------------------------
354\subsection{K Profile Parametrisation (KPP) (\key{zdfkpp}) }
355\label{ZDF_kpp}
356
357%--------------------------------------------namkpp--------------------------------------------------------
358\namdisplay{namzdf_kpp}
359%--------------------------------------------------------------------------------------------------------------
360
361The KKP scheme has been implemented by J. Chanut ...
362
363\colorbox{yellow}{Add a description of KPP here.}
364
365
366% ================================================================
367% Convection
368% ================================================================
369\section{Convection}
370\label{ZDF_conv}
371
372%--------------------------------------------namzdf--------------------------------------------------------
373\namdisplay{namzdf}
374%--------------------------------------------------------------------------------------------------------------
375
376Static instabilities (i.e. light potential densities under heavy ones) may
377occur at particular ocean grid points. In nature, convective processes
378quickly re-establish the static stability of the water column. These
379processes have been removed from the model via the hydrostatic
380assumption so they must be parameterized. Three parameterisations
381are available to deal with convective processes: a non-penetrative
382convective adjustment or an enhanced vertical diffusion, or/and the
383use of a turbulent closure scheme.
384
385% -------------------------------------------------------------------------------------------------------------
386%       Non-Penetrative Convective Adjustment
387% -------------------------------------------------------------------------------------------------------------
388\subsection   [Non-Penetrative Convective Adjustment (\np{ln\_tranpc}) ]
389         {Non-Penetrative Convective Adjustment (\np{ln\_tranpc}=.true.) }
390\label{ZDF_npc}
391
392%--------------------------------------------namzdf--------------------------------------------------------
393\namdisplay{namzdf}
394%--------------------------------------------------------------------------------------------------------------
395
396
397%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
398\begin{figure}[!htb] \label{Fig_npc}   \begin{center}
399\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_npc.pdf}
400\caption {Example of an unstable density profile treated by the non penetrative
401convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from
402the surface to the bottom. It is found to be unstable between levels 3 and 4.
403They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5
404are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are
405mixed. The $1^{st}$ step ends since the density profile is then stable below
406the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same
407procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile
408is checked. It is found stable: end of algorithm.}
409\end{center}   \end{figure}
410%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
411
412The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true.
413It is applied at each \np{nnpc1} time step and mixes downwards instantaneously
414the statically unstable portion of the water column, but only until the density
415structure becomes neutrally stable ($i.e.$ until the mixed portion of the water
416column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}.
417The associated algorithm is an iterative process used in the following way
418(Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is
419found. Assume in the following that the instability is located between levels
420$k$ and $k+1$. The potential temperature and salinity in the two levels are
421vertically mixed, conserving the heat and salt contents of the water column.
422The new density is then computed by a linear approximation. If the new
423density profile is still unstable between levels $k+1$ and $k+2$, levels $k$,
424$k+1$ and $k+2$ are then mixed. This process is repeated until stability is
425established below the level $k$ (the mixing process can go down to the
426ocean bottom). The algorithm is repeated to check if the density profile
427between level $k-1$ and $k$ is unstable and/or if there is no deeper instability.
428
429This algorithm is significantly different from mixing statically unstable levels
430two by two. The latter procedure cannot converge with a finite number
431of iterations for some vertical profiles while the algorithm used in \NEMO 
432converges for any profile in a number of iterations which is less than the
433number of vertical levels. This property is of paramount importance as
434pointed out by \citet{Killworth1989}: it avoids the existence of permanent
435and unrealistic static instabilities at the sea surface. This non-penetrative
436convective algorithm has been proved successful in studies of the deep
437water formation in the north-western Mediterranean Sea
438\citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
439
440Note that in the current implementation of this algorithm presents several
441limitations. First, potential density referenced to the sea surface is used to
442check whether the density profile is stable or not. This is a strong
443simplification which leads to large errors for realistic ocean simulations.
444Indeed, many water masses of the world ocean, especially Antarctic Bottom
445Water, are unstable when represented in surface-referenced potential density.
446The scheme will erroneously mix them up. Second, the mixing of potential
447density is assumed to be linear. This assures the convergence of the algorithm
448even when the equation of state is non-linear. Small static instabilities can thus
449persist due to cabbeling: they will be treated at the next time step.
450Third, temperature and salinity, and thus density, are mixed, but the
451corresponding velocity fields remain unchanged. When using a Richardson
452Number dependent eddy viscosity, the mixing of momentum is done through
453the vertical diffusion: after a static adjustment, the Richardson Number is zero
454and thus the eddy viscosity coefficient is at a maximum. When this convective
455adjustment algorithm is used with constant vertical eddy viscosity, spurious
456solutions can occur since the vertical momentum diffusion remains small even
457after a static adjustment. In that case, we recommend the addition of momentum
458mixing in a manner that mimics the mixing in temperature and salinity
459\citep{Speich_PhD92, Speich_al_JPO96}.
460
461% -------------------------------------------------------------------------------------------------------------
462%       Enhanced Vertical Diffusion
463% -------------------------------------------------------------------------------------------------------------
464\subsection   [Enhanced Vertical Diffusion (\np{ln\_zdfevd})]
465         {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=.true.)}
466\label{ZDF_evd}
467
468%--------------------------------------------namzdf--------------------------------------------------------
469\namdisplay{namzdf}
470%--------------------------------------------------------------------------------------------------------------
471
472The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true.
473In this case, the vertical eddy mixing coefficients are assigned very large values
474(a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable
475($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar_PhD97, Lazar_al_JPO99}.
476This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and
477tracers (\np{n\_evdm}=1).
478
479In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and
480if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
481values also, are set equal to the namelist parameter \np{avevd}. A typical value
482for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of
483convective processes is less time consuming than the convective adjustment
484algorithm presented above when mixing both tracers and momentum in the
485case of static instabilities. It requires the use of an implicit time stepping on
486vertical diffusion terms (i.e. \np{ln\_zdfexp}=false).
487
488% -------------------------------------------------------------------------------------------------------------
489%       Turbulent Closure Scheme
490% -------------------------------------------------------------------------------------------------------------
491\subsection{Turbulent Closure Scheme (\key{zdftke})}
492\label{ZDF_tcs}
493
494The TKE turbulent closure scheme presented in \S\ref{ZDF_tke} and used
495when the \key{zdftke} is defined, in theory solves the problem of statically
496unstable density profiles. In such a case, the term corresponding to the
497destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e} 
498becomes a source term, since $N^2$ is negative. It results in large values of
499$A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring
500$A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values
501restore the static stability of the water column in a way similar to that of the
502enhanced vertical diffusion parameterisation (\S\ref{ZDF_evd}). However,
503in the vicinity of the sea surface (first ocean layer), the eddy coefficients
504computed by the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$,
505because the mixing length scale is bounded by the distance to the sea surface
506(see \S\ref{ZDF_tke}). It can thus be useful to combine the enhanced vertical
507diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} 
508namelist parameter to true and defining the \key{zdftke} CPP key all together.
509
510The KPP turbulent closure scheme already includes enhanced vertical diffusion
511in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ 
512found in \mdl{zdfkpp}, therefore \np{np\_zdfevd} should not be used with the KPP
513scheme. %gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
514
515% ================================================================
516% Double Diffusion Mixing
517% ================================================================
518\section  [Double Diffusion Mixing (\textit{zdfddm} - \key{zdfddm})]
519      {Double Diffusion Mixing (\mdl{zdfddm} module - \key{zdfddm})}
520\label{ZDF_ddm}
521
522%-------------------------------------------namddm--------------------------------------------------------
523\namdisplay{namzdf_ddm}
524%--------------------------------------------------------------------------------------------------------------
525
526Double diffusion occurs when relatively warm, salty water overlies cooler, fresher
527water, or vice versa. The former condition leads to salt fingering and the latter
528to diffusive convection. Double-diffusive phenomena contribute to diapycnal
529mixing in extensive regions of the ocean.  \citet{Merryfield1999} include a
530parameterisation of such phenomena in a global ocean model and show that
531it leads to relatively minor changes in circulation but exerts significant regional
532influences on temperature and salinity.
533
534Diapycnal mixing of S and T are described by diapycnal diffusion coefficients
535\begin{align*} % \label{Eq_zdfddm_Kz}
536    &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT}  \\
537    &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
538\end{align*}
539where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,
540and $o$ by processes other than double diffusion. The rates of double-diffusive mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$,
541where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline
542contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt
543fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981):
544\begin{align} \label{Eq_zdfddm_f}
545A_f^{vS} &=    \begin{cases}
546   \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
547   0                              &\text{otherwise} 
548            \end{cases}   
549\\           \label{Eq_zdfddm_f_T}
550A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
551\end{align}
552
553%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
554\begin{figure}[!t] \label{Fig_zdfddm}  \begin{center}
555\includegraphics[width=0.99\textwidth]{./TexFiles/Figures/Fig_zdfddm.pdf}
556\caption {From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ 
557and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy
558curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves
559$A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and
560$A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy
561curves denote the Federov parameterisation and thin curves the Kelley
562parameterisation. The latter is not implemented in \NEMO. }
563\end{center}    \end{figure}
564%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
565
566The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio
567$\alpha F_T /\beta F_S \approx  0.7$ of buoyancy flux of heat to buoyancy
568flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following  \citet{Merryfield1999},
569we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
570
571To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by Federov (1988) is used:
572\begin{align}  \label{Eq_zdfddm_d}
573A_d^{vT} &=    \begin{cases}
574   1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
575                           &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
576   0                       &\text{otherwise} 
577            \end{cases}   
578\\          \label{Eq_zdfddm_d_S}
579A_d^{vS} &=    \begin{cases}
580   A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right)
581                           &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
582   A_d^{vT} \ 0.15 \ R_\rho
583                           &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
584   0                       &\text{otherwise} 
585            \end{cases}   
586\end{align}
587
588The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ 
589are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing
590$R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the
591same time as $N^2$ is computed. This avoids duplication in the computation of
592$\alpha$ and $\beta$ (which is usually quite expensive).
593
594% ================================================================
595% Bottom Friction
596% ================================================================
597\section  [Bottom Friction (\textit{zdfbfr})]
598      {Bottom Friction (\mdl{zdfbfr} module)}
599\label{ZDF_bfr}
600
601%--------------------------------------------nambfr--------------------------------------------------------
602\namdisplay{namzdf_bfr}
603%--------------------------------------------------------------------------------------------------------------
604
605Both the surface momentum flux (wind stress) and the bottom momentum
606flux (bottom friction) enter the equations as a condition on the vertical
607diffusive flux. For the bottom boundary layer, one has:
608\begin{equation} \label{Eq_zdfbfr_flux}
609A^{vm} \left( \partial \textbf{U}_h / \partial z \right) = \textbf{F}_h
610\end{equation}
611where $\textbf{F}_h$ is supposed to represent the horizontal momentum flux
612outside the logarithmic turbulent boundary layer (thickness of the order of
6131~m in the ocean). How $\textbf{F}_h$ influences the interior depends on the
614vertical resolution of the model near the bottom relative to the Ekman layer
615depth. For example, in order to obtain an Ekman layer depth
616$d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient
617$A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis frequency
618$f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient
619$A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m.
620When the vertical mixing coefficient is this small, using a flux condition is
621equivalent to entering the viscous forces (either wind stress or bottom friction)
622as a body force over the depth of the top or bottom model layer. To illustrate
623this, consider the equation for $u$ at $k$, the last ocean level:
624\begin{equation} \label{Eq_zdfbfr_flux2}
625\frac{\partial u \; (k)}{\partial t} = \frac{1}{e_{3u}} \left[ A^{vm} \; (k) \frac{U \; (k-1) - U \; (k)}{e_{3uw} \; (k-1)} - F_u \right] \approx - \frac{F_u}{e_{3u}}
626\end{equation}
627For example, if the bottom layer thickness is 200~m, the Ekman transport will
628be distributed over that depth. On the other hand, if the vertical resolution
629is high (1~m or less) and a turbulent closure model is used, the turbulent
630Ekman layer will be represented explicitly by the model. However, the
631logarithmic layer is never represented in current primitive equation model
632applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $.
633Two choices are available in \NEMO: a linear and a quadratic bottom friction.
634Note that in both cases, the rotation between the interior velocity and the
635bottom friction is neglected in the present release of \NEMO.
636
637% -------------------------------------------------------------------------------------------------------------
638%       Linear Bottom Friction
639% -------------------------------------------------------------------------------------------------------------
640\subsection{Linear Bottom Friction (\np{nbotfr} = 1) }
641\label{ZDF_bfr_linear}
642
643The linear bottom friction parameterisation assumes that the bottom friction
644is proportional to the interior velocity (i.e. the velocity of the last model level):
645\begin{equation} \label{Eq_zdfbfr_linear}
646\textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
647\end{equation}
648where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean
649layer and $r$ is a friction coefficient expressed in m.s$^{-1}$. This coefficient
650is generally estimated by setting a typical decay time $\tau$ in the deep ocean,
651and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted
652values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}.
653A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used
654in quasi-geostrophic models. One may consider the linear friction as an
655approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982},
656Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed
657of tidal currents of $U_{av} =0.1$~m.s$^{-1}$, and assuming an ocean depth
658$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$.
659This is the default value used in \NEMO. It corresponds to a decay time scale
660of 115~days. It can be changed by specifying \np{bfric1} (namelist parameter).
661
662In the code, the bottom friction is imposed by updating the value of the
663vertical eddy coefficient at the bottom level. Indeed, the discrete formulation
664of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that
665$\textbf {U}_h =0$ below the ocean floor, leads to
666\begin{equation} \label{Eq_zdfbfr_linKz}
667\begin{split}
668A_u^{vm} &= r\;e_{3uw}\\
669A_v^{vm} &= r\;e_{3vw}\\
670\end{split}
671\end{equation}
672
673This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$ 
674used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and
675leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets
676$r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background
677vertical eddy coefficient, and a no-slip boundary condition is imposed.
678Note that this latter choice generally leads to an underestimation of the
679bottom friction: for example with a deepest level thickness of $200~m$ 
680and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient
681is only $r=10^{-6}$m.s$^{-1}$.
682
683% -------------------------------------------------------------------------------------------------------------
684%       Non-Linear Bottom Friction
685% -------------------------------------------------------------------------------------------------------------
686\subsection{Non-Linear Bottom Friction (\np{nbotfr} = 2)}
687\label{ZDF_bfr_nonlinear}
688
689The non-linear bottom friction parameterisation assumes that the bottom
690friction is quadratic:
691\begin{equation} \label{Eq_zdfbfr_nonlinear}
692\textbf {F}_h = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h
693}{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b
694\end{equation}
695
696with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$ 
697the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient,
698and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves
699breaking and other short time scale currents. A typical value of the drag
700coefficient is $C_D = 10^{-3} $. As an example, the CME experiment
701\citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$,
702while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 
703and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been
704set as default values (\np{bfric2} and \np{bfeb2} namelist parameters).
705
706As for the linear case, the bottom friction is imposed in the code by
707updating the value of the vertical eddy coefficient at the bottom level:
708\begin{equation} \label{Eq_zdfbfr_nonlinKz}
709\begin{split}
710A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^
711{1/2}\\
712A_v^{vm} &=C_D\; e_{3uw} \left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
713\end{split}
714\end{equation}
715
716This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the
717non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2},
718and $e_b$ =\np{bfeb2}.
719
720% ================================================================
Note: See TracBrowser for help on using the repository browser.