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

Last change on this file since 11693 was 11693, checked in by nicolasmartin, 13 months ago

Macros renaming

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