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

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

source: branches/2017/dev_merge_2017/DOC/tex_sub/chap_TRA.tex @ 9407

Last change on this file since 9407 was 9407, checked in by nicolasmartin, 6 years ago

Complete refactoring of cross-referencing

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