source: branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex @ 6040

Last change on this file since 6040 was 6040, checked in by gm, 5 years ago

#1613: vvl by default : start to update the DOC for change in vvl, LDF and solvers

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