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 @ 11578

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

New LaTeX commands \nam and \np to mention namelist content (step 2)
So far so good, the manual compiles successfully 4 times and everything seems in its right place

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