source: branches/2017/dev_merge_2017/DOC/TexFiles/Chapters/Chap_TRA.tex @ 9350

Last change on this file since 9350 was 9350, checked in by nicolasmartin, 2 years ago

Add extra message filters and fix too long section titles for ToC

File size: 87.2 KB
Line 
1\documentclass[NEMO_book]{subfiles}
2\begin{document}
3% ================================================================
4% Chapter 1 ——— Ocean Tracers (TRA)
5% ================================================================
6\chapter{Ocean Tracers (TRA)}
7\label{TRA}
8\minitoc
9
10% missing/update
11% traqsr: need to coordinate with SBC module
12
13%STEVEN :  is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? I added a comment to this effect on some instances of this below
14
15%\newpage
16\vspace{2.cm}
17%$\ $\newline    % force a new ligne
18
19Using the representation described in Chap.~\ref{DOM}, several semi-discrete
20space forms of the tracer equations are available depending on the vertical
21coordinate used and on the physics used. In all the equations presented
22here, the masking has been omitted for simplicity. One must be aware that
23all the quantities are masked fields and that each time a mean or difference
24operator is used, the resulting field is multiplied by a mask.
25
26The two active tracers are potential temperature and salinity. Their prognostic
27equations can be summarized as follows:
28\begin{equation*}
29\text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC}
30                   \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP})
31\end{equation*}
32
33NXT stands for next, referring to the time-stepping. From left to right, the terms
34on the rhs of the tracer equations are the advection (ADV), the lateral diffusion
35(LDF), the vertical diffusion (ZDF), the contributions from the external forcings
36(SBC: Surface Boundary Condition, QSR: penetrative Solar Radiation, and BBC:
37Bottom Boundary Condition), the contribution from the bottom boundary Layer
38(BBL) parametrisation, and an internal damping (DMP) term. The terms QSR,
39BBC, BBL and DMP are optional. The external forcings and parameterisations
40require complex inputs and complex calculations ($e.g.$ bulk formulae, estimation
41of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and
42described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively.
43Note that \mdl{tranpc}, the non-penetrative convection module, although
44located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields,
45is described with the model vertical physics (ZDF) together with other available
46parameterization of convection.
47
48In the present chapter we also describe the diagnostic equations used to compute
49the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and
50freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}).
51
52The different options available to the user are managed by namelist logicals or CPP keys.
53For each equation term  \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx},
54where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.
55The CPP key (when it exists) is \textbf{key\_traTTT}. The equivalent code can be
56found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory.
57
58The user has the option of extracting each tendency term on the RHS of the tracer
59equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{DIA}.
60
61$\ $\newline    % force a new ligne
62% ================================================================
63% Tracer Advection
64% ================================================================
65\section  [Tracer Advection (\textit{traadv})]
66      {Tracer Advection (\mdl{traadv})}
67\label{TRA_adv}
68%------------------------------------------namtra_adv-----------------------------------------------------
69\namdisplay{namtra_adv}
70%-------------------------------------------------------------------------------------------------------------
71
72When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \textit{true}),
73the advection tendency of a tracer is expressed in flux form,
74$i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by :
75\begin{equation} \label{Eq_tra_adv}
76ADV_\tau =-\frac{1}{b_t} \left(
77\;\delta _i \left[ e_{2u}\,e_{3u} \;  u\; \tau _u  \right]
78+\delta _j \left[ e_{1v}\,e_{3v}  \;  v\; \tau _v  \right] \; \right)
79-\frac{1}{e_{3t}} \;\delta _k \left[ w\; \tau _w \right]
80\end{equation}
81where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells.
82The flux form in \eqref{Eq_tra_adv} 
83implicitly requires the use of the continuity equation. Indeed, it is obtained
84by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
85which results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 
86(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}=true).
87Therefore it is of paramount importance to design the discrete analogue of the
88advection tendency so that it is consistent with the continuity equation in order to
89enforce the conservation properties of the continuous equations. In other words,
90by setting $\tau = 1$ in (\ref{Eq_tra_adv}) we recover the discrete form of
91the continuity equation which is used to calculate the vertical velocity.
92%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
93\begin{figure}[!t]    \begin{center}
94\includegraphics[width=0.9\textwidth]{Fig_adv_scheme}
95\caption{   \label{Fig_adv_scheme} 
96Schematic representation of some ways used to evaluate the tracer value
97at $u$-point and the amount of tracer exchanged between two neighbouring grid
98points. Upsteam biased scheme (ups): the upstream value is used and the black
99area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation
100is used and the black and dark grey areas are exchanged. Monotonic upstream
101scheme for conservative laws (muscl):  a parabolic interpolation is used and black,
102dark grey and grey areas are exchanged. Second order scheme (cen2): the mean
103value is used and black, dark grey, grey and light grey areas are exchanged. Note
104that this illustration does not include the flux limiter used in ppm and muscl schemes.}
105\end{center}   \end{figure}
106%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
107
108The key difference between the advection schemes available in \NEMO is the choice
109made in space and time interpolation to define the value of the tracer at the
110velocity points (Fig.~\ref{Fig_adv_scheme}).
111
112Along solid lateral and bottom boundaries a zero tracer flux is automatically
113specified, since the normal velocity is zero there. At the sea surface the
114boundary condition depends on the type of sea surface chosen:
115\begin{description}
116\item [linear free surface:] (\np{ln\_linssh}=true) the first level thickness is constant in time:
117the vertical boundary condition is applied at the fixed surface $z=0$ 
118rather than on the moving surface $z=\eta$. There is a non-zero advective
119flux which is set for all advection schemes as
120$\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 
121the product of surface velocity (at $z=0$) by the first level tracer value.
122\item [non-linear free surface:] (\np{ln\_linssh}=false)
123convergence/divergence in the first ocean level moves the free surface
124up/down. There is no tracer advection through it so that the advective
125fluxes through the surface are also zero
126\end{description}
127In all cases, this boundary condition retains local conservation of tracer.
128Global conservation is obtained in non-linear free surface case,
129but \textit{not} in the linear free surface case. Nevertheless, in the latter case,
130it is achieved to a good approximation since the non-conservative
131term is the product of the time derivative of the tracer and the free surface
132height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}.
133
134The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco})
135is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity
136(see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv})
137and/or the mixed layer eddy induced velocity (\textit{eiv})
138when those parameterisations are used (see Chap.~\ref{LDF}).
139
140Several tracer advection scheme are proposed, namely
141a $2^{nd}$ or $4^{th}$ order centred schemes (CEN),
142a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT),
143a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL),
144a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), and
145a Quadratic Upstream Interpolation for Convective Kinematics with
146Estimated Streaming Terms scheme (QUICKEST).
147The choice is made in the \textit{\ngn{namtra\_adv}} namelist, by
148setting to \textit{true} one of the logicals \textit{ln\_traadv\_xxx}.
149The corresponding code can be found in the \textit{traadv\_xxx.F90} module,
150where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme.
151By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals
152are set to \textit{false}. If the user does not select an advection scheme
153in the configuration namelist (\ngn{namelist\_cfg}), the tracers will \textit{not} be advected !
154
155Details of the advection schemes are given below. The choosing an advection scheme
156is a complex matter which depends on the model physics, model resolution,
157type of tracer, as well as the issue of numerical cost. In particular, we note that
158(1) CEN and FCT schemes require an explicit diffusion operator
159while the other schemes are diffusive enough so that they do not necessarily need additional diffusion ;
160(2) CEN and UBS are not \textit{positive} schemes
161\footnote{negative values can appear in an initially strictly positive tracer field
162which is advected}
163, implying that false extrema are permitted. Their use is not recommended on passive tracers ;
164(3) It is recommended that the same advection-diffusion scheme is
165used on both active and passive tracers. Indeed, if a source or sink of a
166passive tracer depends on an active one, the difference of treatment of
167active and passive tracers can create very nice-looking frontal structures
168that are pure numerical artefacts. Nevertheless, most of our users set a different
169treatment on passive and active tracers, that's the reason why this possibility
170is offered. We strongly suggest them to perform a sensitivity experiment
171using a same treatment to assess the robustness of their results.
172
173% -------------------------------------------------------------------------------------------------------------
174%        2nd and 4th order centred schemes
175% -------------------------------------------------------------------------------------------------------------
176\subsection [Centred schemes (CEN) (\np{ln\_traadv\_cen})]
177            {Centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
178\label{TRA_adv_cen}
179
180%        2nd order centred scheme 
181
182The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~=~\textit{true}.
183Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level)
184and vertical direction by setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$.
185CEN implementation can be found in the \mdl{traadv\_cen} module.
186
187In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points
188is evaluated as the mean of the two neighbouring $T$-point values.
189For example, in the $i$-direction :
190\begin{equation} \label{Eq_tra_adv_cen2}
191\tau _u^{cen2} =\overline T ^{i+1/2}
192\end{equation}
193
194CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 
195but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously
196noisy and must be used in conjunction with an explicit diffusion operator to
197produce a sensible solution. The associated time-stepping is performed using
198a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in
199(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
200
201Note that using the CEN2, the overall tracer advection is of second
202order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2})
203have this order of accuracy.
204
205%        4nd order centred scheme 
206
207In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as
208a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points.
209For example, in the $i$-direction:
210\begin{equation} \label{Eq_tra_adv_cen4}
211\tau _u^{cen4} 
212=\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
213\end{equation}
214In the vertical direction (\np{nn\_cen\_v}=$4$), a $4^{th}$ COMPACT interpolation
215has been prefered \citep{Demange_PhD2014}.
216In the COMPACT scheme, both the field and its derivative are interpolated,
217which leads, after a matrix inversion, spectral characteristics
218similar to schemes of higher order \citep{Lele_JCP1992}.
219 
220
221Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme
222but a $4^{th}$ order evaluation of advective fluxes, since the divergence of
223advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order.
224The expression \textit{$4^{th}$ order scheme} used in oceanographic literature
225is usually associated with the scheme presented here.
226Introducing a \textit{true} $4^{th}$ order advection scheme is feasible but,
227for consistency reasons, it requires changes in the discretisation of the tracer
228advection together with changes in the continuity equation,
229and the momentum advection and pressure terms. 
230
231A direct consequence of the pseudo-fourth order nature of the scheme is that
232it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using CEN4.
233Furthermore, it must be used in conjunction with an explicit diffusion operator
234to produce a sensible solution.
235As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction
236with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
237
238At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface),
239an additional hypothesis must be made to evaluate $\tau _u^{cen4}$.
240This hypothesis usually reduces the order of the scheme.
241Here we choose to set the gradient of $T$ across the boundary to zero.
242Alternative conditions can be specified, such as a reduction to a second order scheme
243for these near boundary grid points.
244
245% -------------------------------------------------------------------------------------------------------------
246%        FCT scheme 
247% -------------------------------------------------------------------------------------------------------------
248\subsection   [Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})]
249         {Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct}=true)}
250\label{TRA_adv_tvd}
251
252The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~=~\textit{true}.
253Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level)
254and vertical direction by setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.
255FCT implementation can be found in the \mdl{traadv\_fct} module.
256
257In FCT formulation, the tracer at velocity points is evaluated using a combination of
258an upstream and a centred scheme. For example, in the $i$-direction :
259\begin{equation} \label{Eq_tra_adv_fct}
260\begin{split}
261\tau _u^{ups}&= \begin{cases}
262               T_{i+1}  & \text{if $\ u_{i+1/2} <     0$} \hfill \\
263               T_i         & \text{if $\ u_{i+1/2} \geq 0$} \hfill \\
264              \end{cases}     \\
265\\
266\tau _u^{fct}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen} -\tau _u^{ups} } \right)
267\end{split}
268\end{equation}
269where $c_u$ is a flux limiter function taking values between 0 and 1.
270The FCT order is the one of the centred scheme used ($i.e.$ it depends on the setting of
271\np{nn\_fct\_h} and \np{nn\_fct\_v}.
272There exist many ways to define $c_u$, each corresponding to a different
273FCT scheme. The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.
274$c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field.
275The resulting scheme is quite expensive but \emph{positive}.
276It can be used on both active and passive tracers.
277A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}.
278
279An additional option has been added controlled by \np{nn\_fct\_zts}. By setting this integer to
280a value larger than zero, a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction,
281but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals
282to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited
283by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit
284time stepping should be used on vertical advection of momentum to insure a better stability
285(see \S\ref{DYN_zad}).
286
287For stability reasons (see \S\ref{STP}), $\tau _u^{cen}$ is evaluated in (\ref{Eq_tra_adv_fct})
288using the \textit{now} tracer while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words,
289the advective part of the scheme is time stepped with a leap-frog scheme
290while a forward scheme is used for the diffusive part.
291
292% -------------------------------------------------------------------------------------------------------------
293%        MUSCL scheme 
294% -------------------------------------------------------------------------------------------------------------
295\subsection[MUSCL scheme  (\np{ln\_traadv\_mus})]
296   {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus}=T)}
297\label{TRA_adv_mus}
298
299The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~=~\textit{true}.
300MUSCL implementation can be found in the \mdl{traadv\_mus} module.
301
302MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points
303is evaluated assuming a linear tracer variation between two $T$-points
304(Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
305\begin{equation} \label{Eq_tra_adv_mus}
306   \tau _u^{mus} = \left\{      \begin{aligned}
307         &\tau _&+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
308         &\ \widetilde{\partial _i \tau}  & \quad \text{if }\;u_{i+1/2} \geqslant 0      \\
309         &\tau _{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\rdt}{e_{1u} } \right)
310         &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0
311   \end{aligned}    \right.
312\end{equation}
313where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation
314is imposed to ensure the \textit{positive} character of the scheme.
315
316The time stepping is performed using a forward scheme, that is the \textit{before} 
317tracer field is used to evaluate $\tau _u^{mus}$.
318
319For an ocean grid point adjacent to land and where the ocean velocity is
320directed toward land, an upstream flux is used. This choice ensure
321the \textit{positive} character of the scheme.
322In addition, fluxes round a grid-point where a runoff is applied can optionally be
323computed using upstream fluxes (\np{ln\_mus\_ups}~=~\textit{true}).
324
325% -------------------------------------------------------------------------------------------------------------
326%        UBS scheme 
327% -------------------------------------------------------------------------------------------------------------
328\subsection   [Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs})]
329         {Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=true)}
330\label{TRA_adv_ubs}
331
332The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~=~\textit{true}.
333UBS implementation can be found in the \mdl{traadv\_mus} module.
334
335The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme
336(Quadratic Upstream Interpolation for Convective Kinematics). It is an upstream-biased
337third order scheme based on an upstream-biased parabolic interpolation. 
338For example, in the $i$-direction :
339\begin{equation} \label{Eq_tra_adv_ubs}
340   \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{     
341   \begin{aligned}
342         &\tau"_i          & \quad \text{if }\ u_{i+1/2} \geqslant 0      \\
343         &\tau"_{i+1}   & \quad \text{if }\ u_{i+1/2}       <       0
344   \end{aligned}    \right.
345\end{equation}
346where $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$.
347
348This results in a dissipatively dominant (i.e. hyper-diffusive) truncation
349error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of
350 the advection scheme is similar to that reported in \cite{Farrow1995}.
351It is a relatively good compromise between accuracy and smoothness.
352Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted,
353but the amplitude of such are significantly reduced over the centred second
354or fourth order method. therefore it is not recommended that it should be
355applied to a passive tracer that requires positivity.
356
357The intrinsic diffusion of UBS makes its use risky in the vertical direction
358where the control of artificial diapycnal fluxes is of paramount importance \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}.
359Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme
360or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}=2 or 4).
361
362For stability reasons  (see \S\ref{STP}),
363the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order
364centred scheme) is evaluated using the \textit{now} tracer (centred in time)
365while the second term (which is the diffusive part of the scheme), is
366evaluated using the \textit{before} tracer (forward in time).
367This choice is discussed by \citet{Webb_al_JAOT98} in the context of the
368QUICK advection scheme. UBS and QUICK schemes only differ
369by one coefficient. Replacing 1/6 with 1/8 in \eqref{Eq_tra_adv_ubs} 
370leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
371This option is not available through a namelist parameter, since the
3721/6 coefficient is hard coded. Nevertheless it is quite easy to make the
373substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
374
375Note that it is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
376\begin{equation} \label{Eq_traadv_ubs2}
377\tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{ 
378   \begin{aligned}
379   & + \tau"_i       & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
380   &  - \tau"_{i+1}     & \quad \text{if }\ u_{i+1/2}       <       0
381   \end{aligned}    \right.
382\end{equation}
383or equivalently
384\begin{equation} \label{Eq_traadv_ubs2b}
385u_{i+1/2} \ \tau _u^{ubs} 
386=u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta _i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}
387- \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
388\end{equation}
389
390\eqref{Eq_traadv_ubs2} has several advantages. Firstly, it clearly reveals
391that the UBS scheme is based on the fourth order scheme to which an
392upstream-biased diffusion term is added. Secondly, this emphasises that the
393$4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has
394to be evaluated at the \emph{now} time step using \eqref{Eq_tra_adv_ubs}.
395Thirdly, the diffusion term is in fact a biharmonic operator with an eddy
396coefficient which is simply proportional to the velocity:
397 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO uses
398the computationally more efficient formulation \eqref{Eq_tra_adv_ubs}.
399
400% -------------------------------------------------------------------------------------------------------------
401%        QCK scheme 
402% -------------------------------------------------------------------------------------------------------------
403\subsection   [QUICKEST scheme (QCK) (\np{ln\_traadv\_qck})]
404         {QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=true)}
405\label{TRA_adv_qck}
406
407The Quadratic Upstream Interpolation for Convective Kinematics with
408Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 
409is used when \np{ln\_traadv\_qck}~=~\textit{true}.
410QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
411
412QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST
413limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray
414(MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
415The resulting scheme is quite expensive but \emph{positive}.
416It can be used on both active and passive tracers.
417However, the intrinsic diffusion of QCK makes its use risky in the vertical
418direction where the control of artificial diapycnal fluxes is of paramount importance.
419Therefore the vertical flux is evaluated using the CEN2 scheme.
420This no longer guarantees the positivity of the scheme.
421The use of FCT in the vertical direction (as for the UBS case) should be implemented
422to restore this property.
423
424%%%gmcomment   :  Cross term are missing in the current implementation....
425
426
427% ================================================================
428% Tracer Lateral Diffusion
429% ================================================================
430\section  [Tracer Lateral Diffusion (\textit{traldf})]
431      {Tracer Lateral Diffusion (\mdl{traldf})}
432\label{TRA_ldf}
433%-----------------------------------------nam_traldf------------------------------------------------------
434\namdisplay{namtra_ldf}
435%-------------------------------------------------------------------------------------------------------------
436 
437Options are defined through the \ngn{namtra\_ldf} namelist variables.
438They are regrouped in four items, allowing to specify
439$(i)$   the type of operator used (none, laplacian, bilaplacian),
440$(ii)$  the direction along which the operator acts (iso-level, horizontal, iso-neutral),
441$(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and
442$(iv)$  the specification of eddy diffusivity coefficient (either constant or variable in space and time).
443Item $(iv)$ will be described in Chap.\ref{LDF} .
444The direction along which the operators act is defined through the slope between this direction and the iso-level surfaces.
445The slope is computed in the \mdl{ldfslp} module and will also be described in Chap.~\ref{LDF}.
446
447The lateral diffusion of tracers is evaluated using a forward scheme,
448$i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time,
449except for the pure vertical component that appears when a rotation tensor is used.
450This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}).
451When \np{ln\_traldf\_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which
452the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
453
454% -------------------------------------------------------------------------------------------------------------
455%        Type of operator
456% -------------------------------------------------------------------------------------------------------------
457\subsection   [Type of operator (\np{ln\_traldf\{\_NONE, \_lap, \_blp\}})]
458              {Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, or \np{ln\_traldf\_blp} = true) } 
459\label{TRA_ldf_op}
460
461Three operator options are proposed and, one and only one of them must be selected:
462\begin{description}
463\item [\np{ln\_traldf\_NONE}] = true : no operator selected, the lateral diffusive tendency will not be
464applied to the tracer equation. This option can be used when the selected advection scheme
465is diffusive enough (MUSCL scheme for example).
466\item [ \np{ln\_traldf\_lap}] = true : a laplacian operator is selected. This harmonic operator
467takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $,
468where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}),
469and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}).
470\item [\np{ln\_traldf\_blp}] = true : a bilaplacian operator is selected. This biharmonic operator
471takes the following expression: 
472$\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 
473where the gradient operats along the selected direction,
474and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$  (see Chap.~\ref{LDF}).
475In the code, the bilaplacian operator is obtained by calling the laplacian twice.
476\end{description}
477
478Both laplacian and bilaplacian operators ensure the total tracer variance decrease.
479Their primary role is to provide strong dissipation at the smallest scale supported
480by the grid while minimizing the impact on the larger scale features.
481The main difference between the two operators is the scale selectiveness.
482The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{-4}$ 
483for disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones),
484whereas the laplacian damping time scales only like $\lambda^{-2}$.
485
486
487% -------------------------------------------------------------------------------------------------------------
488%        Direction of action
489% -------------------------------------------------------------------------------------------------------------
490\subsection   [Direction of action (\np{ln\_traldf\{\_lev, \_hor, \_iso, \_triad\}})]
491              {Direction of action (\np{ln\_traldf\_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) } 
492\label{TRA_ldf_dir}
493
494The choice of a direction of action determines the form of operator used.
495The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane
496when iso-level option is used (\np{ln\_traldf\_lev}~=~\textit{true})
497or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate
498(\np{ln\_traldf\_hor} and \np{ln\_zco} equal \textit{true}).
499The associated code can be found in the \mdl{traldf\_lap\_blp} module.
500The operator is a rotated (re-entrant) laplacian when the direction along which it acts
501does not coincide with the iso-level surfaces,
502that is when standard or triad iso-neutral option is used (\np{ln\_traldf\_iso} or
503 \np{ln\_traldf\_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.),
504or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate
505(\np{ln\_traldf\_hor} and \np{ln\_sco} equal \textit{true})
506\footnote{In this case, the standard iso-neutral operator will be automatically selected}.
507In that case, a rotation is applied to the gradient(s) that appears in the operator
508so that diffusive fluxes acts on the three spatial direction.
509
510The resulting discret form of the three operators (one iso-level and two rotated one)
511is given in the next two sub-sections.
512
513
514% -------------------------------------------------------------------------------------------------------------
515%       iso-level operator
516% -------------------------------------------------------------------------------------------------------------
517\subsection   [Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso})]
518         {Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso}) }
519\label{TRA_ldf_lev}
520
521The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by:
522\begin{equation} \label{Eq_tra_ldf_lap}
523D_t^{lT} =\frac{1}{b_t} \left( \;
524   \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right]
525+ \delta _{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right\;\right)
526\end{equation}
527where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells
528and where zero diffusive fluxes is assumed across solid boundaries,
529first (and third in bilaplacian case) horizontal tracer derivative are masked.
530It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module.
531The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} 
532in order to compute the iso-level bilaplacian operator.
533
534It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate
535with or without partial steps, but is simply an iso-level operator in the $s$-coordinate.
536It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~=~\textit{true},
537we have \np{ln\_traldf\_lev}~=~\textit{true} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~=~\textit{true}.
538In both cases, it significantly contributes to diapycnal mixing.
539It is therefore never recommended, even when using it in the bilaplacian case.
540
541Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally
542adjacent cells are located at different depths in the vicinity of the bottom.
543In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level
544require a specific treatment. They are calculated in the \mdl{zpshde} module,
545described in \S\ref{TRA_zpshde}.
546
547
548% -------------------------------------------------------------------------------------------------------------
549%         Rotated laplacian operator
550% -------------------------------------------------------------------------------------------------------------
551\subsection   [Standard and triad rotated (bi-)laplacian operator]
552               {Standard and triad (bi-)laplacian operator}
553\label{TRA_ldf_iso_triad}
554
555%&&    Standard rotated (bi-)laplacian operator
556%&& ----------------------------------------------
557\subsubsection   [Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})]
558                 {Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})}
559\label{TRA_ldf_iso}
560The general form of the second order lateral tracer subgrid scale physics
561(\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates:
562\begin{equation} \label{Eq_tra_ldf_iso}
563\begin{split}
564 D_T^{lT} = \frac{1}{b_t}   & \left\{   \,\;\delta_i \left[   A_u^{lT}  \left(
565     \frac{e_{2u}\,e_{3u}}{e_{1u}} \,\delta_{i+1/2}[T]
566   - e_{2u}\;r_{1u} \,\overline{\overline{ \delta_{k+1/2}[T] }}^{\,i+1/2,k}
567                                                     \right)   \right]   \right.    \\ 
568&             +\delta_j \left[ A_v^{lT} \left(
569          \frac{e_{1v}\,e_{3v}}{e_{2v}}  \,\delta_{j+1/2} [T]
570        - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k} 
571                                                    \right)   \right]                 \\ 
572& +\delta_k \left[ A_w^{lT} \left(
573       -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2}
574                                                    \right.   \right.                 \\ 
575& \qquad \qquad \quad 
576        - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2}     \\
577& \left. {\left. {   \qquad \qquad \ \ \ \left. {
578        +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right)
579           \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 
580 \end{split}
581 \end{equation}
582where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells,
583$r_1$ and $r_2$ are the slopes between the surface of computation
584($z$- or $s$-surfaces) and the surface along which the diffusion operator
585acts ($i.e.$ horizontal or iso-neutral surfaces).  It is thus used when,
586in addition to \np{ln\_traldf\_lap}= true, we have \np{ln\_traldf\_iso}=true,
587or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true. The way these
588slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom
589and lateral boundaries, the turbulent fluxes of heat and salt are set to zero
590using the mask technique (see \S\ref{LBC_coast}).
591
592The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical
593derivatives. For numerical stability, the vertical second derivative must
594be solved using the same implicit time scheme as that used in the vertical
595physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term
596is not computed in the \mdl{traldf\_iso} module, but in the \mdl{trazdf} module
597where, if iso-neutral mixing is used, the vertical mixing coefficient is simply
598increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$.
599
600This formulation conserves the tracer but does not ensure the decrease
601of the tracer variance. Nevertheless the treatment performed on the slopes
602(see \S\ref{LDF}) allows the model to run safely without any additional
603background horizontal diffusion \citep{Guilyardi_al_CD01}.
604
605Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal derivatives
606at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific treatment.
607They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
608
609%&&     Triad rotated (bi-)laplacian operator
610%&&  -------------------------------------------
611\subsubsection   [Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})]
612                 {Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})}
613\label{TRA_ldf_triad}
614
615If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}=true ; see App.\ref{sec:triad})
616
617An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
618is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of
619the algorithm is given in App.\ref{sec:triad}.
620
621The lateral fourth order bilaplacian operator on tracers is obtained by
622applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption
623on boundary conditions: both first and third derivative terms normal to the
624coast are set to zero.
625
626The lateral fourth order operator formulation on tracers is obtained by
627applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption
628on boundary conditions: first and third derivative terms normal to the
629coast, normal to the bottom and normal to the surface are set to zero.
630
631%&&    Option for the rotated operators
632%&& ----------------------------------------------
633\subsubsection   [Option for the rotated operators]
634                 {Option for the rotated operators}
635\label{TRA_ldf_options}
636
637\np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
638
639\np{rn\_slpmax} = slope limit (both operators)
640
641\np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
642
643\np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only)
644
645\np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
646
647% ================================================================
648% Tracer Vertical Diffusion
649% ================================================================
650\section  [Tracer Vertical Diffusion (\textit{trazdf})]
651      {Tracer Vertical Diffusion (\mdl{trazdf})}
652\label{TRA_zdf}
653%--------------------------------------------namzdf---------------------------------------------------------
654\namdisplay{namzdf}
655%--------------------------------------------------------------------------------------------------------------
656
657Options are defined through the \ngn{namzdf} namelist variables.
658The formulation of the vertical subgrid scale tracer physics is the same
659for all the vertical coordinates, and is based on a laplacian operator.
660The vertical diffusion operator given by (\ref{Eq_PE_zdf}) takes the
661following semi-discrete space form:
662\begin{equation} \label{Eq_tra_zdf}
663\begin{split}
664D^{vT}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T] \;\right]
665\\
666D^{vS}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S] \;\right]
667\end{split}
668\end{equation}
669where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity
670coefficients on temperature and salinity, respectively. Generally,
671$A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is
672parameterised ($i.e.$ \key{zdfddm} is defined). The way these coefficients
673are evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when
674iso-neutral mixing is used, both mixing coefficients are increased
675by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ 
676to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}.
677
678At the surface and bottom boundaries, the turbulent fluxes of
679heat and salt must be specified. At the surface they are prescribed
680from the surface forcing and added in a dedicated routine (see \S\ref{TRA_sbc}),
681whilst at the bottom they are set to zero for heat and salt unless
682a geothermal flux forcing is prescribed as a bottom boundary
683condition (see \S\ref{TRA_bbc}).
684
685The large eddy coefficient found in the mixed layer together with high
686vertical resolution implies that in the case of explicit time stepping
687(\np{ln\_zdfexp}=true) there would be too restrictive a constraint on
688the time step. Therefore, the default implicit time stepping is preferred
689for the vertical diffusion since it overcomes the stability constraint.
690A forward time differencing scheme (\np{ln\_zdfexp}=true) using a time
691splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative.
692Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both
693tracers and dynamics.
694
695% ================================================================
696% External Forcing
697% ================================================================
698\section{External Forcing}
699\label{TRA_sbc_qsr_bbc}
700
701% -------------------------------------------------------------------------------------------------------------
702%        surface boundary condition
703% -------------------------------------------------------------------------------------------------------------
704\subsection   [Surface boundary condition (\textit{trasbc})]
705         {Surface boundary condition (\mdl{trasbc})}
706\label{TRA_sbc}
707
708The surface boundary condition for tracers is implemented in a separate
709module (\mdl{trasbc}) instead of entering as a boundary condition on the vertical
710diffusion operator (as in the case of momentum). This has been found to
711enhance readability of the code. The two formulations are completely
712equivalent; the forcing terms in trasbc are the surface fluxes divided by
713the thickness of the top model layer.
714
715Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components
716($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer
717of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
718and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$,
719the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details).
720By doing this, the forcing formulation is the same for any tracer (including temperature and salinity).
721
722The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following
723forcing fields (used on tracers):
724
725$\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface
726(i.e. the difference between the total surface heat flux and the fraction of the short wave flux that
727penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with
728of the mass exchange with the atmosphere and lands.
729
730$\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...)
731
732$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)
733 and possibly with the sea-ice and ice-shelves.
734
735$\bullet$ \textit{rnf}, the mass flux associated with runoff
736(see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
737
738$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt,
739(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied).
740
741The surface boundary condition on temperature and salinity is applied as follows:
742\begin{equation} \label{Eq_tra_sbc}
743\begin{aligned}
744 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^& \\ 
745& F^S =\frac{ 1 }{\rho _\,      \left. e_{3t} \right|_{k=1} }  &\overline{ \textit{sfx} }^t   & \\   
746 \end{aligned}
747\end{equation} 
748where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps
749($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the
750divergence of odd and even time step (see \S\ref{STP}).
751
752In the linear free surface case (\np{ln\_linssh}~=~\textit{true}),
753an additional term has to be added on both temperature and salinity.
754On temperature, this term remove the heat content associated with mass exchange
755that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that
756would have resulted from a change in the volume of the first level.
757The resulting surface boundary condition is applied as follows:
758\begin{equation} \label{Eq_tra_sbc_lin}
759\begin{aligned}
760 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }   
761           &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^& \\ 
762%
763& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} } 
764           &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1}  \right) }^t   & \\   
765 \end{aligned}
766\end{equation} 
767Note that an exact conservation of heat and salt content is only achieved with non-linear free surface.
768In the linear free surface case, there is a small imbalance. The imbalance is larger
769than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.
770This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}).
771
772% -------------------------------------------------------------------------------------------------------------
773%        Solar Radiation Penetration
774% -------------------------------------------------------------------------------------------------------------
775\subsection   [Solar Radiation Penetration (\textit{traqsr})]
776         {Solar Radiation Penetration (\mdl{traqsr})}
777\label{TRA_qsr}
778%--------------------------------------------namqsr--------------------------------------------------------
779\namdisplay{namtra_qsr}
780%--------------------------------------------------------------------------------------------------------------
781
782Options are defined through the  \ngn{namtra\_qsr} namelist variables.
783When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true),
784the solar radiation penetrates the top few tens of meters of the ocean. If it is not used
785(\np{ln\_flxqsr}=false) all the heat flux is absorbed in the first ocean level.
786Thus, in the former case a term is added to the time evolution equation of
787temperature \eqref{Eq_PE_tra_T} and the surface boundary condition is
788modified to take into account only the non-penetrative part of the surface
789heat flux:
790\begin{equation} \label{Eq_PE_qsr}
791\begin{split}
792\frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}   \\
793Q_{ns} &= Q_\text{Total} - Q_{sr}
794\end{split}
795\end{equation}
796where $Q_{sr}$ is the penetrative part of the surface heat flux ($i.e.$ the shortwave radiation)
797and $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$).
798The additional term in \eqref{Eq_PE_qsr} is discretized as follows:
799\begin{equation} \label{Eq_tra_qsr}
800\frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k} \equiv \frac{1}{\rho_o\, C_p\, e_{3t}} \delta_k \left[ I_w \right]
801\end{equation}
802
803The shortwave radiation,  $Q_{sr}$, consists of energy distributed across a wide spectral range.
804The ocean is strongly absorbing for wavelengths longer than 700~nm and these
805wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$ 
806that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified
807through namelist parameter \np{rn\_abs}).  It is assumed to penetrate the ocean
808with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$,
809of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namelist).
810For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy
811propagates to larger depths where it contributes to
812local heating.
813The way this second part of the solar energy penetrates into the ocean depends on
814which formulation is chosen. In the simple 2-waveband light penetration scheme  (\np{ln\_qsr\_2bd}=true)
815a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,
816leading to the following expression  \citep{Paulson1977}:
817\begin{equation} \label{Eq_traqsr_iradiance}
818I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right]
819\end{equation}
820where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 
821It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter.
822The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in
823Jerlov's (1968) classification (oligotrophic waters).
824
825Such assumptions have been shown to provide a very crude and simplistic
826representation of observed light penetration profiles (\cite{Morel_JGR88}, see also
827Fig.\ref{Fig_traqsr_irradiance}). Light absorption in the ocean depends on
828particle concentration and is spectrally selective. \cite{Morel_JGR88} has shown
829that an accurate representation of light penetration can be provided by a 61 waveband
830formulation. Unfortunately, such a model is very computationally expensive.
831Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this
832formulation in which visible light is split into three wavebands: blue (400-500 nm),
833green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependent
834attenuation coefficient is fitted to the coefficients computed from the full spectral model
835of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), assuming
836the same power-law relationship. As shown in Fig.\ref{Fig_traqsr_irradiance},
837this formulation, called RGB (Red-Green-Blue), reproduces quite closely
838the light penetration profiles predicted by the full spectal model, but with much greater
839computational efficiency. The 2-bands formulation does not reproduce the full model very well.
840
841The RGB formulation is used when \np{ln\_qsr\_rgb}=true. The RGB attenuation coefficients
842($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform
843chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 
844in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation:
845\begin{description} 
846\item[\np{nn\_chdta}=0] 
847a constant 0.05 g.Chl/L value everywhere ;
848\item[\np{nn\_chdta}=1] 
849an observed time varying chlorophyll deduced from satellite surface ocean color measurement
850spread uniformly in the vertical direction ;
851\item[\np{nn\_chdta}=2] 
852same as previous case except that a vertical profile of chlorophyl is used.
853Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ;
854\item[\np{ln\_qsr\_bio}=true] 
855simulated time varying chlorophyll by TOP biogeochemical model.
856In this case, the RGB formulation is used to calculate both the phytoplankton
857light limitation in PISCES or LOBSTER and the oceanic heating rate.
858\end{description} 
859The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation
860is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}.
861
862When the $z$-coordinate is preferred to the $s$-coordinate, the depth of $w-$levels does
863not significantly vary with location. The level at which the light has been totally
864absorbed ($i.e.$ it is less than the computer precision) is computed once,
865and the trend associated with the penetration of the solar radiation is only added down to that level.
866Finally, note that when the ocean is shallow ($<$ 200~m), part of the
867solar radiation can reach the ocean floor. In this case, we have
868chosen that all remaining radiation is absorbed in the last ocean
869level ($i.e.$ $I$ is masked).
870
871%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
872\begin{figure}[!t]     \begin{center}
873\includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance}
874\caption{    \label{Fig_traqsr_irradiance}
875Penetration profile of the downward solar irradiance calculated by four models.
876Two waveband chlorophyll-independent formulation (blue), a chlorophyll-dependent
877monochromatic formulation (green), 4 waveband RGB formulation (red),
87861 waveband Morel (1988) formulation (black) for a chlorophyll concentration of
879(a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. From \citet{Lengaigne_al_CD07}.}
880\end{center}   \end{figure}
881%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
882
883% -------------------------------------------------------------------------------------------------------------
884%        Bottom Boundary Condition
885% -------------------------------------------------------------------------------------------------------------
886\subsection   [Bottom Boundary Condition (\textit{trabbc})]
887         {Bottom Boundary Condition (\mdl{trabbc})}
888\label{TRA_bbc}
889%--------------------------------------------nambbc--------------------------------------------------------
890\namdisplay{nambbc}
891%--------------------------------------------------------------------------------------------------------------
892%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
893\begin{figure}[!t]     \begin{center}
894\includegraphics[width=1.0\textwidth]{Fig_TRA_geoth}
895\caption{   \label{Fig_geothermal}
896Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
897It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.}
898\end{center}   \end{figure}
899%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
900
901Usually it is assumed that there is no exchange of heat or salt through
902the ocean bottom, $i.e.$ a no flux boundary condition is applied on active
903tracers at the bottom. This is the default option in \NEMO, and it is
904implemented using the masking technique. However, there is a
905non-zero heat flux across the seafloor that is associated with solid
906earth cooling. This flux is weak compared to surface fluxes (a mean
907global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), but it warms
908systematically the ocean and acts on the densest water masses.
909Taking this flux into account in a global ocean model increases
910the deepest overturning cell ($i.e.$ the one associated with the Antarctic
911Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OS09}.
912
913Options are defined through the  \ngn{namtra\_bbc} namelist variables.
914The presence of geothermal heating is controlled by setting the namelist
915parameter  \np{ln\_trabbc} to true. Then, when \np{nn\_geoflx} is set to 1,
916a constant geothermal heating is introduced whose value is given by the
917\np{nn\_geoflx\_cst}, which is also a namelist parameter.
918When  \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is
919introduced which is provided in the \ifile{geothermal\_heating} NetCDF file
920(Fig.\ref{Fig_geothermal}) \citep{Emile-Geay_Madec_OS09}.
921
922% ================================================================
923% Bottom Boundary Layer
924% ================================================================
925\section  [Bottom Boundary Layer (\mdl{trabbl} - \key{trabbl})]
926      {Bottom Boundary Layer (\mdl{trabbl} - \key{trabbl})}
927\label{TRA_bbl}
928%--------------------------------------------nambbl---------------------------------------------------------
929\namdisplay{nambbl}
930%--------------------------------------------------------------------------------------------------------------
931
932Options are defined through the  \ngn{nambbl} namelist variables.
933In a $z$-coordinate configuration, the bottom topography is represented by a
934series of discrete steps. This is not adequate to represent gravity driven
935downslope flows. Such flows arise either downstream of sills such as the Strait of
936Gibraltar or Denmark Strait, where dense water formed in marginal seas flows
937into a basin filled with less dense water, or along the continental slope when dense
938water masses are formed on a continental shelf. The amount of entrainment
939that occurs in these gravity plumes is critical in determining the density
940and volume flux of the densest waters of the ocean, such as Antarctic Bottom Water,
941or North Atlantic Deep Water. $z$-coordinate models tend to overestimate the
942entrainment, because the gravity flow is mixed vertically by convection
943as it goes ''downstairs'' following the step topography, sometimes over a thickness
944much larger than the thickness of the observed gravity plume. A similar problem
945occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly
946downstream of a sill \citep{Willebrand_al_PO01}, and the thickness
947of the plume is not resolved.
948
949The idea of the bottom boundary layer (BBL) parameterisation, first introduced by
950\citet{Beckmann_Doscher1997}, is to allow a direct communication between
951two adjacent bottom cells at different levels, whenever the densest water is
952located above the less dense water. The communication can be by a diffusive flux
953(diffusive BBL), an advective flux (advective BBL), or both. In the current
954implementation of the BBL, only the tracers are modified, not the velocities.
955Furthermore, it only connects ocean bottom cells, and therefore does not include
956all the improvements introduced by \citet{Campin_Goosse_Tel99}.
957
958% -------------------------------------------------------------------------------------------------------------
959%        Diffusive BBL
960% -------------------------------------------------------------------------------------------------------------
961\subsection{Diffusive Bottom Boundary layer (\np{nn\_bbl\_ldf}=1)}
962\label{TRA_bbl_diff}
963
964When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1),
965the diffusive flux between two adjacent cells at the ocean floor is given by
966\begin{equation} \label{Eq_tra_bbl_diff}
967{\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T
968\end{equation} 
969with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells,
970and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997},
971the latter is prescribed with a spatial dependence, $i.e.$ in the conditional form
972\begin{equation} \label{Eq_tra_bbl_coef}
973A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
974 A_{bbl}  \quad \quad   \mbox{if}  \quad   \nabla_\sigma \rho  \cdot  \nabla H<0 \\ 
975 \\
976 0\quad \quad \;\,\mbox{otherwise} \\ 
977 \end{array}} \right.
978\end{equation} 
979where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist
980parameter \np{rn\_ahtbbl} and usually set to a value much larger
981than the one used for lateral mixing in the open ocean. The constraint in \eqref{Eq_tra_bbl_coef} 
982implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of
983the slope, is larger than in the deeper ocean (see green arrow in Fig.\ref{Fig_bbl}).
984In practice, this constraint is applied separately in the two horizontal directions,
985and the density gradient in \eqref{Eq_tra_bbl_coef} is evaluated with the log gradient formulation:
986\begin{equation} \label{Eq_tra_bbl_Drho}
987   \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S
988\end{equation} 
989where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$,
990$\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature,
991salinity and depth, respectively.
992
993% -------------------------------------------------------------------------------------------------------------
994%        Advective BBL
995% -------------------------------------------------------------------------------------------------------------
996\subsection   {Advective Bottom Boundary Layer  (\np{nn\_bbl\_adv}= 1 or 2)}
997\label{TRA_bbl_adv}
998
999\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following
1000if this is not what is meant then "downwards sloping flow" is also a possibility"}
1001
1002%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1003\begin{figure}[!t]   \begin{center}
1004\includegraphics[width=0.7\textwidth]{Fig_BBL_adv}
1005\caption{   \label{Fig_bbl} 
1006Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is
1007activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$.
1008Red arrows indicate the additional overturning circulation due to the advective BBL.
1009The transport of the downslope flow is defined either as the transport of the bottom
1010ocean cell (black arrow), or as a function of the along slope density gradient.
1011The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$
1012ocean bottom cells.
1013connection}
1014\end{center}   \end{figure}
1015%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1016
1017
1018%!!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
1019%!!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
1020%!!        i.e. transport proportional to the along-slope density gradient
1021
1022%%%gmcomment   :  this section has to be really written
1023
1024When applying an advective BBL (\np{nn\_bbl\_adv} = 1 or 2), an overturning
1025circulation is added which connects two adjacent bottom grid-points only if dense
1026water overlies less dense water on the slope. The density difference causes dense
1027water to move down the slope.
1028
1029\np{nn\_bbl\_adv} = 1 : the downslope velocity is chosen to be the Eulerian
1030ocean velocity just above the topographic step (see black arrow in Fig.\ref{Fig_bbl})
1031\citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection
1032is allowed only if dense water overlies less dense water on the slope ($i.e.$ 
1033$\nabla_\sigma \rho  \cdot  \nabla H<0$) and if the velocity is directed towards
1034greater depth ($i.e.$ $\vect{U}  \cdot  \nabla H>0$).
1035
1036\np{nn\_bbl\_adv} = 2 : the downslope velocity is chosen to be proportional to $\Delta \rho$,
1037the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
1038The advection is allowed only  if dense water overlies less dense water on the slope ($i.e.$ 
1039$\nabla_\sigma \rho  \cdot  \nabla H<0$). For example, the resulting transport of the
1040downslope flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by the
1041following expression:
1042\begin{equation} \label{Eq_bbl_Utr}
1043 u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right)
1044\end{equation}
1045where $\gamma$, expressed in seconds, is the coefficient of proportionality
1046provided as \np{rn\_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn} 
1047are the vertical index of the higher and lower cells, respectively.
1048The parameter $\gamma$ should take a different value for each bathymetric
1049step, but for simplicity, and because no direct estimation of this parameter is
1050available, a uniform value has been assumed. The possible values for $\gamma$ 
1051range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}
1052
1053Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ 
1054using the upwind scheme. Such a diffusive advective scheme has been chosen
1055to mimic the entrainment between the downslope plume and the surrounding
1056water at intermediate depths. The entrainment is replaced by the vertical mixing
1057implicit in the advection scheme. Let us consider as an example the
1058case displayed in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is
1059larger than the one at level $(i,kdwn)$. The advective BBL scheme
1060modifies the tracer time tendency of the ocean cells near the
1061topographic step by the downslope flow \eqref{Eq_bbl_dw},
1062the horizontal \eqref{Eq_bbl_hor}  and the upward \eqref{Eq_bbl_up} 
1063return flows as follows:
1064\begin{align} 
1065\partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw}
1066                                     +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right\label{Eq_bbl_dw} \\
1067%
1068\partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 
1069               + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{Eq_bbl_hor} \\
1070%
1071\intertext{and for $k =kdw-1,\;..., \; kup$ :} 
1072%
1073\partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k}
1074               + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{Eq_bbl_up}
1075\end{align}
1076where $b_t$ is the $T$-cell volume.
1077
1078Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in
1079the model outputs. It has to be used to compute the effective velocity
1080as well as the effective overturning circulation.
1081
1082% ================================================================
1083% Tracer damping
1084% ================================================================
1085\section  [Tracer damping (\textit{tradmp})]
1086      {Tracer damping (\mdl{tradmp})}
1087\label{TRA_dmp}
1088%--------------------------------------------namtra_dmp-------------------------------------------------
1089\namdisplay{namtra_dmp}
1090%--------------------------------------------------------------------------------------------------------------
1091
1092In some applications it can be useful to add a Newtonian damping term
1093into the temperature and salinity equations:
1094\begin{equation} \label{Eq_tra_dmp}
1095\begin{split}
1096 \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right\\
1097 \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right)
1098 \end{split}
1099 \end{equation} 
1100where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ 
1101are given temperature and salinity fields (usually a climatology).
1102Options are defined through the  \ngn{namtra\_dmp} namelist variables.
1103The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true.
1104It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true
1105in \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are
1106correctly set  ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read
1107using \mdl{fldread}, see \S\ref{SBC_fldread}).
1108The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. The file name is specified by the namelist variable \np{cn\_resto}. The DMP\_TOOLS tool is provided to allow users to generate the netcdf file.
1109
1110The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} 
1111the specification of the boundary conditions along artificial walls of a
1112limited domain basin and \textit{(b)} the computation of the velocity
1113field associated with a given $T$-$S$ field (for example to build the
1114initial state of a prognostic simulation, or to use the resulting velocity
1115field for a passive tracer study). The first case applies to regional
1116models that have artificial walls instead of open boundaries.
1117In the vicinity of these walls, $\gamma$ takes large values (equivalent to
1118a time scale of a few days) whereas it is zero in the interior of the
1119model domain. The second case corresponds to the use of the robust
1120diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity
1121field consistent with the model dynamics whilst having a $T$, $S$ field
1122close to a given climatological field ($T_o$, $S_o$).
1123
1124The robust diagnostic method is very efficient in preventing temperature
1125drift in intermediate waters but it produces artificial sources of heat and salt
1126within the ocean. It also has undesirable effects on the ocean convection.
1127It tends to prevent deep convection and subsequent deep-water formation,
1128by stabilising the water column too much.
1129
1130The namelist parameter \np{nn\_zdmp} sets whether the damping should be applied in the whole water column or only below the mixed layer (defined either on a density or $S_o$ criterion). It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here \citep{Madec_al_JPO96}.
1131
1132\subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS}
1133
1134DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$.
1135Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled
1136and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input.
1137This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1.
1138The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work.
1139The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
1140
1141%--------------------------------------------nam_dmp_create-------------------------------------------------
1142\namtools{namelist_dmp}
1143%-------------------------------------------------------------------------------------------------------
1144
1145\np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list.
1146
1147The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations.
1148\np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain.
1149\np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea
1150for the ORCA4, ORCA2 and ORCA05 configurations.
1151If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as
1152a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference
1153configurations with previous model versions.
1154\np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines.
1155This option only has an effect if \np{ln\_full\_field} is true.
1156\np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer.
1157Finally \np{ln\_custom} specifies that the custom module will be called.
1158This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
1159
1160The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}.
1161Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to
1162the full values of a 10\deg latitud band.
1163This is often used because of the short adjustment time scale in the equatorial region
1164\citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a
1165hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}
1166
1167% ================================================================
1168% Tracer time evolution
1169% ================================================================
1170\section  [Tracer time evolution (\textit{tranxt})]
1171      {Tracer time evolution (\mdl{tranxt})}
1172\label{TRA_nxt}
1173%--------------------------------------------namdom-----------------------------------------------------
1174\namdisplay{namdom}
1175%--------------------------------------------------------------------------------------------------------------
1176
1177Options are defined through the  \ngn{namdom} namelist variables.
1178The general framework for tracer time stepping is a modified leap-frog scheme
1179\citep{Leclair_Madec_OM09}, $i.e.$ a three level centred time scheme associated
1180with a Asselin time filter (cf. \S\ref{STP_mLF}):
1181\begin{equation} \label{Eq_tra_nxt}
1182\begin{aligned}
1183(e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt  \,e_{3t}^t\ \text{RHS}^t & \\
1184\\
1185(e_{3t}T)_f^\;\ \quad &= (e_{3t}T)^t \;\quad 
1186                                    &+\gamma \,\left[ {(e_{3t}T)_f^{t-\rdt} -2(e_{3t}T)^t+(e_{3t}T)^{t+\rdt}} \right] &  \\
1187                                 & &- \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]  &                     
1188\end{aligned}
1189\end{equation} 
1190where RHS is the right hand side of the temperature equation,
1191the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient,
1192and $S$ is the total forcing applied on $T$ ($i.e.$ fluxes plus content in mass exchanges).
1193$\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter).
1194Its default value is \np{rn\_atfp}=$10^{-3}$. Note that the forcing correction term in the filter
1195is not applied in linear free surface (\jp{lk\_vvl}=false) (see \S\ref{TRA_sbc}.
1196Not also that in constant volume case, the time stepping is performed on $T$,
1197not on its content, $e_{3t}T$.
1198
1199When the vertical mixing is solved implicitly, the update of the \textit{next} tracer
1200fields is done in module \mdl{trazdf}. In this case only the swapping of arrays
1201and the Asselin filtering is done in the \mdl{tranxt} module.
1202
1203In order to prepare for the computation of the \textit{next} time step,
1204a swap of tracer arrays is performed: $T^{t-\rdt} = T^t$ and $T^t = T_f$.
1205
1206% ================================================================
1207% Equation of State (eosbn2)
1208% ================================================================
1209\section  [Equation of State (\textit{eosbn2}) ]
1210      {Equation of State (\mdl{eosbn2}) }
1211\label{TRA_eosbn2}
1212%--------------------------------------------nameos-----------------------------------------------------
1213\namdisplay{nameos}
1214%--------------------------------------------------------------------------------------------------------------
1215
1216% -------------------------------------------------------------------------------------------------------------
1217%        Equation of State
1218% -------------------------------------------------------------------------------------------------------------
1219\subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)}
1220\label{TRA_eos}
1221
1222The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship
1223linking seawater density, $\rho$, to a number of state variables,
1224most typically temperature, salinity and pressure.
1225Because density gradients control the pressure gradient force through the hydrostatic balance,
1226the equation of state provides a fundamental bridge between the distribution of active tracers
1227and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular
1228influencing the circulation through determination of the static stability below the mixed layer,
1229thus controlling rates of exchange between the atmosphere  and the ocean interior \citep{Roquet_JPO2015}.
1230Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983})
1231or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real
1232ocean circulation is attempted \citep{Roquet_JPO2015}.
1233The use of TEOS-10 is highly recommended because
1234\textit{(i)} it is the new official EOS,
1235\textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and
1236\textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature
1237and practical salinity for EOS-980, both variables being more suitable for use as model variables
1238\citep{TEOS10, Graham_McDougall_JPO13}.
1239EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility.
1240For process studies, it is often convenient to use an approximation of the EOS. To that purposed,
1241a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available.
1242
1243In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$,
1244is computed, with $\rho_o$ a reference density. Called \textit{rau0} 
1245in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$.
1246This is a sensible choice for the reference density used in a Boussinesq ocean
1247climate model, as, with the exception of only a small percentage of the ocean,
1248density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
1249
1250Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 
1251which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS).
1252\begin{description}
1253
1254\item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 
1255The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
1256but it is optimized for a boussinesq fluid and the polynomial expressions have simpler
1257and more computationally efficient expressions for their derived quantities
1258which make them more adapted for use in ocean models.
1259Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10
1260rational function approximation for hydrographic data analysis  \citep{TEOS10}.
1261A key point is that conservative state variables are used:
1262Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$).
1263The pressure in decibars is approximated by the depth in meters.
1264With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to
1265$C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}.
1266
1267Choosing polyTEOS10-bsq implies that the state variables used by the model are
1268$\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as
1269\textit{Conservative} Temperature and \textit{Absolute} Salinity.
1270In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST
1271prior to either computing the air-sea and ice-sea fluxes (forced mode)
1272or sending the SST field to the atmosphere (coupled mode).
1273
1274\item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used.
1275It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized
1276to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80
1277and the ocean model are:
1278the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $^{\circ}C$, notation: $\theta$).
1279The pressure in decibars is approximated by the depth in meters. 
1280With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature,
1281salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to
1282have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant
1283value, the TEOS10 value.
1284 
1285\item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,
1286the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.)
1287(see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both
1288cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS
1289in theoretical studies \citep{Roquet_JPO2015}.
1290With such an equation of state there is no longer a distinction between
1291\textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute} 
1292and \textit{practical} salinity.
1293S-EOS takes the following expression:
1294\begin{equation} \label{Eq_tra_S-EOS}
1295\begin{split}
1296  d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_\\
1297                                & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_\\
1298                                & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\
1299  with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3
1300\end{split}
1301\end{equation} 
1302where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}.
1303In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing
1304the associated coefficients.
1305Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS.
1306setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS.
1307Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S.
1308
1309\end{description}
1310
1311
1312%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1313\begin{table}[!tb]
1314\begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|}
1315\hline
1316coeff.   & computer name   & S-EOS     &  description                      \\ \hline
1317$a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff.    \\ \hline
1318$b_0$       & \np{rn\_b0}     & 7.6554 $10^{-1}$ &  linear haline  expansion coeff.    \\ \hline
1319$\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$          \\ \hline
1320$\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$       \\ \hline
1321$\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$       \\ \hline
1322$\mu_1$     & \np{rn\_mu1}    & 1.4970 $10^{-4}$ &  thermobaric coeff. in T         \\ \hline
1323$\mu_2$     & \np{rn\_mu2}    & 1.1090 $10^{-5}$ &  thermobaric coeff. in S            \\ \hline
1324\end{tabular}
1325\caption{ \label{Tab_SEOS}
1326Standard value of S-EOS coefficients. }
1327\end{center}
1328\end{table}
1329%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1330
1331
1332% -------------------------------------------------------------------------------------------------------------
1333%        Brunt-V\"{a}is\"{a}l\"{a} Frequency
1334% -------------------------------------------------------------------------------------------------------------
1335\subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
1336\label{TRA_bn2}
1337
1338An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a}
1339 frequency) is of paramount importance as determine the ocean stratification and
1340 is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent
1341 vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing
1342 parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure
1343 (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 
1344 is given by:
1345\begin{equation} \label{Eq_tra_bn2}
1346N^2 = \frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
1347\end{equation} 
1348where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS,
1349and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.
1350The coefficients are a polynomial function of temperature, salinity and depth which expression
1351depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran} 
1352function that can be found in \mdl{eosbn2}.
1353
1354% -------------------------------------------------------------------------------------------------------------
1355%        Freezing Point of Seawater
1356% -------------------------------------------------------------------------------------------------------------
1357\subsection   [Freezing Point of Seawater]
1358         {Freezing Point of Seawater}
1359\label{TRA_fzp}
1360
1361The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}:
1362\begin{equation} \label{Eq_tra_eos_fzp}
1363   \begin{split}
1364T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} 
1365                       -  2.154996 \;10^{-4} \,\right) \ S    \\
1366               - 7.53\,10^{-3} \ \ p
1367   \end{split}
1368\end{equation}
1369
1370\eqref{Eq_tra_eos_fzp} is only used to compute the potential freezing point of
1371sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent
1372terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing
1373point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found
1374in \mdl{eosbn2}
1375
1376
1377% -------------------------------------------------------------------------------------------------------------
1378%        Potential Energy     
1379% -------------------------------------------------------------------------------------------------------------
1380%\subsection{Potential Energy anomalies}
1381%\label{TRA_bn2}
1382
1383%    =====>>>>> TO BE written
1384%
1385
1386
1387% ================================================================
1388% Horizontal Derivative in zps-coordinate
1389% ================================================================
1390\section  [Horizontal Derivative in \textit{zps}-coordinate (\textit{zpshde})]
1391      {Horizontal Derivative in \textit{zps}-coordinate (\mdl{zpshde})}
1392\label{TRA_zpshde}
1393
1394\gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators,
1395                   I've changed "derivative" to "difference" and "mean" to "average"}
1396
1397With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general,
1398tracers in horizontally adjacent cells live at different depths.
1399Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module)
1400and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module).
1401The partial cell properties at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom.
1402So, only the bottom interpolation is explained below.
1403
1404Before taking horizontal gradients between the tracers next to the bottom, a linear
1405interpolation in the vertical is used to approximate the deeper tracer as if it actually
1406lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}).
1407For example, for temperature in the $i$-direction the needed interpolated
1408temperature, $\widetilde{T}$, is:
1409
1410%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1411\begin{figure}[!p]    \begin{center}
1412\includegraphics[width=0.9\textwidth]{Partial_step_scheme}
1413\caption{   \label{Fig_Partial_step_scheme} 
1414Discretisation of the horizontal difference and average of tracers in the $z$-partial
1415step coordinate (\np{ln\_zps}=true) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$.
1416A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value
1417at the depth of the shallower tracer point of the two adjacent bottom $T$-points.
1418The horizontal difference is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ 
1419and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$}
1420\end{center}   \end{figure}
1421%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1422\begin{equation*}
1423\widetilde{T}= \left\{  \begin{aligned} 
1424&T^{\,i+1}      -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta _k T^{i+1} 
1425                        && \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$   }  \\
1426                              \\
1427&T^{\,i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i       }\;\delta _k T^{i+1}
1428                        && \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
1429            \end{aligned}   \right.
1430\end{equation*}
1431and the resulting forms for the horizontal difference and the horizontal average
1432value of $T$ at a $U$-point are:
1433\begin{equation} \label{Eq_zps_hde}
1434\begin{aligned}
1435 \delta _{i+1/2} T=  \begin{cases}
1436\ \ \ \widetilde {T}\quad\ -T^i     & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\
1437                              \\
1438\ \ \ T^{\,i+1}-\widetilde{T}    & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
1439                  \end{cases}     \\
1440\\
1441\overline {T}^{\,i+1/2} \ =   \begin{cases}
1442( \widetilde {T}\ \ \;\,-T^{\,i})    / 2  & \;\ \ \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\
1443                              \\
1444( T^{\,i+1}-\widetilde{T} ) / 2     & \;\ \ \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
1445            \end{cases}
1446\end{aligned}
1447\end{equation}
1448
1449The computation of horizontal derivative of tracers as well as of density is
1450performed once for all at each time step in \mdl{zpshde} module and stored
1451in shared arrays to be used when needed. It has to be emphasized that the
1452procedure used to compute the interpolated density, $\widetilde{\rho}$, is not
1453the same as that used for $T$ and $S$. Instead of forming a linear approximation
1454of density, we compute $\widetilde{\rho }$ from the interpolated values of $T$ 
1455and $S$, and the pressure at a $u$-point (in the equation of state pressure is
1456approximated by depth, see \S\ref{TRA_eos} ) :
1457\begin{equation} \label{Eq_zps_hde_rho}
1458\widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u })
1459\quad \text{where }\  z_u = \min \left( {z_T^{i+1} ,z_T^i } \right)
1460\end{equation} 
1461
1462This is a much better approximation as the variation of $\rho$ with depth (and
1463thus pressure) is highly non-linear with a true equation of state and thus is badly
1464approximated with a linear interpolation. This approximation is used to compute
1465both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral
1466surfaces (\S\ref{LDF_slp})
1467
1468Note that in almost all the advection schemes presented in this Chapter, both
1469averaging and differencing operators appear. Yet \eqref{Eq_zps_hde} has not
1470been used in these schemes: in contrast to diffusion and pressure gradient
1471computations, no correction for partial steps is applied for advection. The main
1472motivation is to preserve the domain averaged mean variance of the advected
1473field when using the $2^{nd}$ order centred scheme. Sensitivity of the advection
1474schemes to the way horizontal averages are performed in the vicinity of partial
1475cells should be further investigated in the near future.
1476%%%
1477\gmcomment{gm :   this last remark has to be done}
1478%%%
1479\end{document}
Note: See TracBrowser for help on using the repository browser.