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

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

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

Last change on this file since 11435 was 11435, checked in by nicolasmartin, 5 years ago

Various corrections on chapters

Cleaning the indexes by fixing/removing wrong entries (or appending a ? to unknown items) and
improve the classification with new index definitions for CPP keys and namelist blocks:

  • from \key{...} cmd, key_ prefix no longer precedes the index entry
  • namelist block declaration moves from \ngn{nam...} to \nam{...} (i.e. \ngn{namtra\_ldf} -> \nam{tra\_ldf}) The expected prefix nam is added to the printed word but not the index entry.

Now we have indexes with a better sorting instead of all CPP keys under 'K' and namelists blocks under 'N'.

Fix missing space issues with alias commands by adding a trailing backslash (\NEMO\, \ie\, \eg\, ...).
There is no perfect solution for this, and I prefer not using a particular package to solve it.

Review the initial LaTeX code snippet for the historic changes in chapters

Finally, for readability and future diff visualisations, please avoid writing paragraphs with continuous lines.
Break the lines around 80 to 100 characters long

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