source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex @ 10442

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

Front page edition, cleaning in custom LaTeX commands and add index for single subfile compilation

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