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

Last change on this file since 14257 was 14257, checked in by nicolasmartin, 4 months ago

Overall review of LaTeX sources (not tested completely as of now):

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