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

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

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

Last change on this file since 11598 was 11598, checked in by nicolasmartin, 3 years ago

Add template of versioning record at the beginning of chapters

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