source: trunk/DOC/TexFiles/Chapters/Chap_TRA.tex @ 6497

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

#1720 - trunk: add Casimir tidal parameterization

File size: 87.4 KB
1% ================================================================
2% Chapter 1 ——— Ocean Tracers (TRA)
3% ================================================================
4\chapter{Ocean Tracers (TRA)}
8% missing/update
9% traqsr: need to coordinate with SBC module
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
15%$\ $\newline    % force a new ligne
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.
24The two active tracers are potential temperature and salinity. Their prognostic
25equations can be summarized as follows:
27\text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC}
28                   \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP})
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.
46In the present chapter we also describe the diagnostic equations used to compute
47the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and
48freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}).
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.
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{DIA}.
59$\ $\newline    % force a new ligne
60% ================================================================
61% Tracer Advection
62% ================================================================
63\section  [Tracer Advection (\textit{traadv})]
64      {Tracer Advection (\mdl{traadv})}
70When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \textit{true}),
71the advection tendency of a tracer is expressed in flux form,
72$i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by :
73\begin{equation} \label{Eq_tra_adv}
74ADV_\tau =-\frac{1}{b_t} \left(
75\;\delta _i \left[ e_{2u}\,e_{3u} \;  u\; \tau _u  \right]
76+\delta _j \left[ e_{1v}\,e_{3v}  \;  v\; \tau _v  \right] \; \right)
77-\frac{1}{e_{3t}} \;\delta _k \left[ w\; \tau _w \right]
79where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells.
80The flux form in \eqref{Eq_tra_adv} 
81implicitly requires the use of the continuity equation. Indeed, it is obtained
82by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
83which results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 
84(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}=true).
85Therefore it is of paramount importance to design the discrete analogue of the
86advection tendency so that it is consistent with the continuity equation in order to
87enforce the conservation properties of the continuous equations. In other words,
88by setting $\tau = 1$ in (\ref{Eq_tra_adv}) we recover the discrete form of
89the continuity equation which is used to calculate the vertical velocity.
91\begin{figure}[!t]    \begin{center}
93\caption{   \label{Fig_adv_scheme} 
94Schematic representation of some ways used to evaluate the tracer value
95at $u$-point and the amount of tracer exchanged between two neighbouring grid
96points. Upsteam biased scheme (ups): the upstream value is used and the black
97area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation
98is used and the black and dark grey areas are exchanged. Monotonic upstream
99scheme for conservative laws (muscl):  a parabolic interpolation is used and black,
100dark grey and grey areas are exchanged. Second order scheme (cen2): the mean
101value is used and black, dark grey, grey and light grey areas are exchanged. Note
102that this illustration does not include the flux limiter used in ppm and muscl schemes.}
103\end{center}   \end{figure}
106The key difference between the advection schemes available in \NEMO is the choice
107made in space and time interpolation to define the value of the tracer at the
108velocity points (Fig.~\ref{Fig_adv_scheme}).
110Along solid lateral and bottom boundaries a zero tracer flux is automatically
111specified, since the normal velocity is zero there. At the sea surface the
112boundary condition depends on the type of sea surface chosen:
114\item [linear free surface:] (\np{ln\_linssh}=true) the first level thickness is constant in time:
115the vertical boundary condition is applied at the fixed surface $z=0$ 
116rather than on the moving surface $z=\eta$. There is a non-zero advective
117flux which is set for all advection schemes as
118$\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 
119the product of surface velocity (at $z=0$) by the first level tracer value.
120\item [non-linear free surface:] (\np{ln\_linssh}=false)
121convergence/divergence in the first ocean level moves the free surface
122up/down. There is no tracer advection through it so that the advective
123fluxes through the surface are also zero
125In all cases, this boundary condition retains local conservation of tracer.
126Global conservation is obtained in non-linear free surface case,
127but \textit{not} in the linear free surface case. Nevertheless, in the latter case,
128it is achieved to a good approximation since the non-conservative
129term is the product of the time derivative of the tracer and the free surface
130height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}.
132The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco})
133is the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity
134(see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv})
135and/or the mixed layer eddy induced velocity (\textit{eiv})
136when those parameterisations are used (see Chap.~\ref{LDF}).
138Several tracer advection scheme are proposed, namely
139a $2^{nd}$ or $4^{th}$ order centred schemes (CEN),
140a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT),
141a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL),
142a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), and
143a Quadratic Upstream Interpolation for Convective Kinematics with
144Estimated Streaming Terms scheme (QUICKEST).
145The choice is made in the \textit{\ngn{namtra\_adv}} namelist, by
146setting to \textit{true} one of the logicals \textit{ln\_traadv\_xxx}.
147The corresponding code can be found in the \textit{traadv\_xxx.F90} module,
148where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme.
149By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals
150are set to \textit{false}. If the user does not select an advection scheme
151in the configuration namelist (\ngn{namelist\_cfg}), the tracers will \textit{not} be advected !
153Details of the advection schemes are given below. The choosing an advection scheme
154is a complex matter which depends on the model physics, model resolution,
155type of tracer, as well as the issue of numerical cost. In particular, we note that
156(1) CEN and FCT schemes require an explicit diffusion operator
157while the other schemes are diffusive enough so that they do not necessarily need additional diffusion ;
158(2) CEN and UBS are not \textit{positive} schemes
159\footnote{negative values can appear in an initially strictly positive tracer field
160which is advected}
161, implying that false extrema are permitted. Their use is not recommended on passive tracers ;
162(3) It is recommended that the same advection-diffusion scheme is
163used on both active and passive tracers. Indeed, if a source or sink of a
164passive tracer depends on an active one, the difference of treatment of
165active and passive tracers can create very nice-looking frontal structures
166that are pure numerical artefacts. Nevertheless, most of our users set a different
167treatment on passive and active tracers, that's the reason why this possibility
168is offered. We strongly suggest them to perform a sensitivity experiment
169using a same treatment to assess the robustness of their results.
171% -------------------------------------------------------------------------------------------------------------
172%        2nd and 4th order centred schemes
173% -------------------------------------------------------------------------------------------------------------
174\subsection [Centred schemes (CEN) (\np{ln\_traadv\_cen})]
175            {Centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
178%        2nd order centred scheme 
180The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~=~\textit{true}.
181Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level)
182and vertical direction by setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$.
183CEN implementation can be found in the \mdl{traadv\_cen} module.
185In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points
186is evaluated as the mean of the two neighbouring $T$-point values.
187For example, in the $i$-direction :
188\begin{equation} \label{Eq_tra_adv_cen2}
189\tau _u^{cen2} =\overline T ^{i+1/2}
192CEN2 is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ 
193but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously
194noisy and must be used in conjunction with an explicit diffusion operator to
195produce a sensible solution. The associated time-stepping is performed using
196a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in
197(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
199Note that using the CEN2, the overall tracer advection is of second
200order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2})
201have this order of accuracy.
203%        4nd order centred scheme 
205In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as
206a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points.
207For example, in the $i$-direction:
208\begin{equation} \label{Eq_tra_adv_cen4}
209\tau _u^{cen4} 
210=\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
212In the vertical direction (\np{nn\_cen\_v}=$4$), a $4^{th}$ COMPACT interpolation
213has been prefered \citep{Demange_PhD2014}.
214In the COMPACT scheme, both the field and its derivative are interpolated,
215which leads, after a matrix inversion, spectral characteristics
216similar to schemes of higher order \citep{Lele_JCP1992}.
219Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme
220but a $4^{th}$ order evaluation of advective fluxes, since the divergence of
221advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order.
222The expression \textit{$4^{th}$ order scheme} used in oceanographic literature
223is usually associated with the scheme presented here.
224Introducing a \textit{true} $4^{th}$ order advection scheme is feasible but,
225for consistency reasons, it requires changes in the discretisation of the tracer
226advection together with changes in the continuity equation,
227and the momentum advection and pressure terms. 
229A direct consequence of the pseudo-fourth order nature of the scheme is that
230it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using CEN4.
231Furthermore, it must be used in conjunction with an explicit diffusion operator
232to produce a sensible solution.
233As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction
234with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
236At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface),
237an additional hypothesis must be made to evaluate $\tau _u^{cen4}$.
238This hypothesis usually reduces the order of the scheme.
239Here we choose to set the gradient of $T$ across the boundary to zero.
240Alternative conditions can be specified, such as a reduction to a second order scheme
241for these near boundary grid points.
243% -------------------------------------------------------------------------------------------------------------
244%        FCT scheme 
245% -------------------------------------------------------------------------------------------------------------
246\subsection   [Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})]
247         {Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct}=true)}
250The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~=~\textit{true}.
251Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level)
252and vertical direction by setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.
253FCT implementation can be found in the \mdl{traadv\_fct} module.
255In FCT formulation, the tracer at velocity points is evaluated using a combination of
256an upstream and a centred scheme. For example, in the $i$-direction :
257\begin{equation} \label{Eq_tra_adv_fct}
259\tau _u^{ups}&= \begin{cases}
260               T_{i+1}  & \text{if $\ u_{i+1/2} <     0$} \hfill \\
261               T_i         & \text{if $\ u_{i+1/2} \geq 0$} \hfill \\
262              \end{cases}     \\
264\tau _u^{fct}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen} -\tau _u^{ups} } \right)
267where $c_u$ is a flux limiter function taking values between 0 and 1.
268The FCT order is the one of the centred scheme used ($i.e.$ it depends on the setting of
269\np{nn\_fct\_h} and \np{nn\_fct\_v}.
270There exist many ways to define $c_u$, each corresponding to a different
271FCT scheme. The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.
272$c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field.
273The resulting scheme is quite expensive but \emph{positive}.
274It can be used on both active and passive tracers.
275A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}.
277An additional option has been added controlled by \np{nn\_fct\_zts}. By setting this integer to
278a value larger than zero, a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction,
279but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals
280to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited
281by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit
282time stepping should be used on vertical advection of momentum to insure a better stability
283(see \S\ref{DYN_zad}).
285For stability reasons (see \S\ref{STP}), $\tau _u^{cen}$ is evaluated in (\ref{Eq_tra_adv_fct})
286using the \textit{now} tracer while $\tau _u^{ups}$ is evaluated using the \textit{before} tracer. In other words,
287the advective part of the scheme is time stepped with a leap-frog scheme
288while a forward scheme is used for the diffusive part.
290% -------------------------------------------------------------------------------------------------------------
291%        MUSCL scheme 
292% -------------------------------------------------------------------------------------------------------------
293\subsection[MUSCL scheme  (\np{ln\_traadv\_mus})]
294   {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus}=T)}
297The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~=~\textit{true}.
298MUSCL implementation can be found in the \mdl{traadv\_mus} module.
300MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points
301is evaluated assuming a linear tracer variation between two $T$-points
302(Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
303\begin{equation} \label{Eq_tra_adv_mus}
304   \tau _u^{mus} = \left\{      \begin{aligned}
305         &\tau _&+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
306         &\ \widetilde{\partial _i \tau}  & \quad \text{if }\;u_{i+1/2} \geqslant 0      \\
307         &\tau _{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\rdt}{e_{1u} } \right)
308         &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0
309   \end{aligned}    \right.
311where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation
312is imposed to ensure the \textit{positive} character of the scheme.
314The time stepping is performed using a forward scheme, that is the \textit{before} 
315tracer field is used to evaluate $\tau _u^{mus}$.
317For an ocean grid point adjacent to land and where the ocean velocity is
318directed toward land, an upstream flux is used. This choice ensure
319the \textit{positive} character of the scheme.
320In addition, fluxes round a grid-point where a runoff is applied can optionally be
321computed using upstream fluxes (\np{ln\_mus\_ups}~=~\textit{true}).
323% -------------------------------------------------------------------------------------------------------------
324%        UBS scheme 
325% -------------------------------------------------------------------------------------------------------------
326\subsection   [Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs})]
327         {Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=true)}
330The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~=~\textit{true}.
331UBS implementation can be found in the \mdl{traadv\_mus} module.
333The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme
334(Quadratic Upstream Interpolation for Convective Kinematics). It is an upstream-biased
335third order scheme based on an upstream-biased parabolic interpolation. 
336For example, in the $i$-direction :
337\begin{equation} \label{Eq_tra_adv_ubs}
338   \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{     
339   \begin{aligned}
340         &\tau"_i          & \quad \text{if }\ u_{i+1/2} \geqslant 0      \\
341         &\tau"_{i+1}   & \quad \text{if }\ u_{i+1/2}       <       0
342   \end{aligned}    \right.
344where $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$.
346This results in a dissipatively dominant (i.e. hyper-diffusive) truncation
347error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of
348 the advection scheme is similar to that reported in \cite{Farrow1995}.
349It is a relatively good compromise between accuracy and smoothness.
350Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted,
351but the amplitude of such are significantly reduced over the centred second
352or fourth order method. therefore it is not recommended that it should be
353applied to a passive tracer that requires positivity.
355The intrinsic diffusion of UBS makes its use risky in the vertical direction
356where the control of artificial diapycnal fluxes is of paramount importance \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}.
357Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme
358or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}=2 or 4).
360For stability reasons  (see \S\ref{STP}),
361the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order
362centred scheme) is evaluated using the \textit{now} tracer (centred in time)
363while the second term (which is the diffusive part of the scheme), is
364evaluated using the \textit{before} tracer (forward in time).
365This choice is discussed by \citet{Webb_al_JAOT98} in the context of the
366QUICK advection scheme. UBS and QUICK schemes only differ
367by one coefficient. Replacing 1/6 with 1/8 in \eqref{Eq_tra_adv_ubs} 
368leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
369This option is not available through a namelist parameter, since the
3701/6 coefficient is hard coded. Nevertheless it is quite easy to make the
371substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
373Note that it is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
374\begin{equation} \label{Eq_traadv_ubs2}
375\tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{ 
376   \begin{aligned}
377   & + \tau"_i       & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
378   &  - \tau"_{i+1}     & \quad \text{if }\ u_{i+1/2}       <       0
379   \end{aligned}    \right.
381or equivalently
382\begin{equation} \label{Eq_traadv_ubs2b}
383u_{i+1/2} \ \tau _u^{ubs} 
384=u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta _i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}
385- \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
388\eqref{Eq_traadv_ubs2} has several advantages. Firstly, it clearly reveals
389that the UBS scheme is based on the fourth order scheme to which an
390upstream-biased diffusion term is added. Secondly, this emphasises that the
391$4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has
392to be evaluated at the \emph{now} time step using \eqref{Eq_tra_adv_ubs}.
393Thirdly, the diffusion term is in fact a biharmonic operator with an eddy
394coefficient which is simply proportional to the velocity:
395 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO uses
396the computationally more efficient formulation \eqref{Eq_tra_adv_ubs}.
398% -------------------------------------------------------------------------------------------------------------
399%        QCK scheme 
400% -------------------------------------------------------------------------------------------------------------
401\subsection   [QUICKEST scheme (QCK) (\np{ln\_traadv\_qck})]
402         {QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=true)}
405The Quadratic Upstream Interpolation for Convective Kinematics with
406Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 
407is used when \np{ln\_traadv\_qck}~=~\textit{true}.
408QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
410QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST
411limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray
412(MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
413The resulting scheme is quite expensive but \emph{positive}.
414It can be used on both active and passive tracers.
415However, the intrinsic diffusion of QCK makes its use risky in the vertical
416direction where the control of artificial diapycnal fluxes is of paramount importance.
417Therefore the vertical flux is evaluated using the CEN2 scheme.
418This no longer guarantees the positivity of the scheme.
419The use of FCT in the vertical direction (as for the UBS case) should be implemented
420to restore this property.
422%%%gmcomment   :  Cross term are missing in the current implementation....
425% ================================================================
426% Tracer Lateral Diffusion
427% ================================================================
428\section  [Tracer Lateral Diffusion (\textit{traldf})]
429      {Tracer Lateral Diffusion (\mdl{traldf})}
435Options are defined through the \ngn{namtra\_ldf} namelist variables.
436They are regrouped in four items, allowing to specify
437$(i)$   the type of operator used (none, laplacian, bilaplacian),
438$(ii)$  the direction along which the operator acts (iso-level, horizontal, iso-neutral),
439$(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and
440$(iv)$  the specification of eddy diffusivity coefficient (either constant or variable in space and time).
441Item $(iv)$ will be described in Chap.\ref{LDF} .
442The direction along which the operators act is defined through the slope between this direction and the iso-level surfaces.
443The slope is computed in the \mdl{ldfslp} module and will also be described in Chap.~\ref{LDF}.
445The lateral diffusion of tracers is evaluated using a forward scheme,
446$i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time,
447except for the pure vertical component that appears when a rotation tensor is used.
448This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}).
449When \np{ln\_traldf\_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which
450the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
452% -------------------------------------------------------------------------------------------------------------
453%        Type of operator
454% -------------------------------------------------------------------------------------------------------------
455\subsection   [Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, \np{ln\_traldf\_blp})]
456              {Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, or \np{ln\_traldf\_blp} = true) } 
459Three operator options are proposed and, one and only one of them must be selected:
461\item [\np{ln\_traldf\_NONE}] = true : no operator selected, the lateral diffusive tendency will not be
462applied to the tracer equation. This option can be used when the selected advection scheme
463is diffusive enough (MUSCL scheme for example).
464\item [ \np{ln\_traldf\_lap}] = true : a laplacian operator is selected. This harmonic operator
465takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $,
466where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}),
467and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}).
468\item [\np{ln\_traldf\_blp}] = true : a bilaplacian operator is selected. This biharmonic operator
469takes the following expression: 
470$\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 
471where the gradient operats along the selected direction,
472and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$  (see Chap.~\ref{LDF}).
473In the code, the bilaplacian operator is obtained by calling the laplacian twice.
476Both laplacian and bilaplacian operators ensure the total tracer variance decrease.
477Their primary role is to provide strong dissipation at the smallest scale supported
478by the grid while minimizing the impact on the larger scale features.
479The main difference between the two operators is the scale selectiveness.
480The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{-4}$ 
481for disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones),
482whereas the laplacian damping time scales only like $\lambda^{-2}$.
485% -------------------------------------------------------------------------------------------------------------
486%        Direction of action
487% -------------------------------------------------------------------------------------------------------------
488\subsection   [Direction of action (\np{ln\_traldf\_lev}, \np{ln\_traldf\_hor}, \np{ln\_traldf\_iso}, \np{ln\_traldf\_triad})]
489              {Direction of action (\np{ln\_traldf\_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) } 
492The choice of a direction of action determines the form of operator used.
493The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane
494when iso-level option is used (\np{ln\_traldf\_lev}~=~\textit{true})
495or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate
496(\np{ln\_traldf\_hor} and \np{ln\_zco} equal \textit{true}).
497The associated code can be found in the \mdl{traldf\_lap\_blp} module.
498The operator is a rotated (re-entrant) laplacian when the direction along which it acts
499does not coincide with the iso-level surfaces,
500that is when standard or triad iso-neutral option is used (\np{ln\_traldf\_iso} or
501 \np{ln\_traldf\_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.),
502or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate
503(\np{ln\_traldf\_hor} and \np{ln\_sco} equal \textit{true})
504\footnote{In this case, the standard iso-neutral operator will be automatically selected}.
505In that case, a rotation is applied to the gradient(s) that appears in the operator
506so that diffusive fluxes acts on the three spatial direction.
508The resulting discret form of the three operators (one iso-level and two rotated one)
509is given in the next two sub-sections.
512% -------------------------------------------------------------------------------------------------------------
513%       iso-level operator
514% -------------------------------------------------------------------------------------------------------------
515\subsection   [Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso})]
516         {Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso}) }
519The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by:
520\begin{equation} \label{Eq_tra_ldf_lap}
521D_t^{lT} =\frac{1}{b_t} \left( \;
522   \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right]
523+ \delta _{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right\;\right)
525where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells
526and where zero diffusive fluxes is assumed across solid boundaries,
527first (and third in bilaplacian case) horizontal tracer derivative are masked.
528It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module.
529The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} 
530in order to compute the iso-level bilaplacian operator.
532It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate
533with or without partial steps, but is simply an iso-level operator in the $s$-coordinate.
534It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~=~\textit{true},
535we have \np{ln\_traldf\_lev}~=~\textit{true} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~=~\textit{true}.
536In both cases, it significantly contributes to diapycnal mixing.
537It is therefore never recommended, even when using it in the bilaplacian case.
539Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally
540adjacent cells are located at different depths in the vicinity of the bottom.
541In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level
542require a specific treatment. They are calculated in the \mdl{zpshde} module,
543described in \S\ref{TRA_zpshde}.
546% -------------------------------------------------------------------------------------------------------------
547%         Rotated laplacian operator
548% -------------------------------------------------------------------------------------------------------------
549\subsection   [Standard and triad rotated (bi-)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad})]
550               {Standard and triad (bi-)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad}))}
553%&&    Standard rotated (bi-)laplacian operator
554%&& ----------------------------------------------
555\subsubsection   [Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})]
556                 {Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})}
558The general form of the second order lateral tracer subgrid scale physics
559(\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates:
560\begin{equation} \label{Eq_tra_ldf_iso}
562 D_T^{lT} = \frac{1}{b_t}   & \left\{   \,\;\delta_i \left[   A_u^{lT}  \left(
563     \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}
565                                                     \right)   \right]   \right.    \\ 
566&             +\delta_j \left[ A_v^{lT} \left(
567          \frac{e_{1v}\,e_{3v}}{e_{2v}}  \,\delta_{j+1/2} [T]
568        - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k} 
569                                                    \right)   \right]                 \\ 
570& +\delta_k \left[ A_w^{lT} \left(
571       -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2}
572                                                    \right.   \right.                 \\ 
573& \qquad \qquad \quad 
574        - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2}     \\
575& \left. {\left. {   \qquad \qquad \ \ \ \left. {
576        +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right)
577           \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 
578 \end{split}
579 \end{equation}
580where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells,
581$r_1$ and $r_2$ are the slopes between the surface of computation
582($z$- or $s$-surfaces) and the surface along which the diffusion operator
583acts ($i.e.$ horizontal or iso-neutral surfaces).  It is thus used when,
584in addition to \np{ln\_traldf\_lap}= true, we have \np{ln\_traldf\_iso}=true,
585or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true. The way these
586slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom
587and lateral boundaries, the turbulent fluxes of heat and salt are set to zero
588using the mask technique (see \S\ref{LBC_coast}).
590The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical
591derivatives. For numerical stability, the vertical second derivative must
592be solved using the same implicit time scheme as that used in the vertical
593physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term
594is not computed in the \mdl{traldf\_iso} module, but in the \mdl{trazdf} module
595where, if iso-neutral mixing is used, the vertical mixing coefficient is simply
596increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$.
598This formulation conserves the tracer but does not ensure the decrease
599of the tracer variance. Nevertheless the treatment performed on the slopes
600(see \S\ref{LDF}) allows the model to run safely without any additional
601background horizontal diffusion \citep{Guilyardi_al_CD01}.
603Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal derivatives
604at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific treatment.
605They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
607%&&     Triad rotated (bi-)laplacian operator
608%&&  -------------------------------------------
609\subsubsection   [Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})]
610                 {Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})}
613If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}=true ; see App.\ref{sec:triad})
615An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
616is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of
617the algorithm is given in App.\ref{sec:triad}.
619The lateral fourth order bilaplacian operator on tracers is obtained by
620applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption
621on boundary conditions: both first and third derivative terms normal to the
622coast are set to zero.
624The lateral fourth order operator formulation on tracers is obtained by
625applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption
626on boundary conditions: first and third derivative terms normal to the
627coast, normal to the bottom and normal to the surface are set to zero.
629%&&    Option for the rotated operators
630%&& ----------------------------------------------
631\subsubsection   [Option for the rotated operators]
632                 {Option for the rotated operators}
635\np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
637\np{rn\_slpmax} = slope limit (both operators)
639\np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
641\np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only)
643\np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
645% ================================================================
646% Tracer Vertical Diffusion
647% ================================================================
648\section  [Tracer Vertical Diffusion (\textit{trazdf})]
649      {Tracer Vertical Diffusion (\mdl{trazdf})}
655Options are defined through the \ngn{namzdf} namelist variables.
656The formulation of the vertical subgrid scale tracer physics is the same
657for all the vertical coordinates, and is based on a laplacian operator.
658The vertical diffusion operator given by (\ref{Eq_PE_zdf}) takes the
659following semi-discrete space form:
660\begin{equation} \label{Eq_tra_zdf}
662D^{vT}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T] \;\right]
664D^{vS}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S] \;\right]
667where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity
668coefficients on temperature and salinity, respectively. Generally,
669$A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is
670parameterised ($i.e.$ \key{zdfddm} is defined). The way these coefficients
671are evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when
672iso-neutral mixing is used, both mixing coefficients are increased
673by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ 
674to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}.
676At the surface and bottom boundaries, the turbulent fluxes of
677heat and salt must be specified. At the surface they are prescribed
678from the surface forcing and added in a dedicated routine (see \S\ref{TRA_sbc}),
679whilst at the bottom they are set to zero for heat and salt unless
680a geothermal flux forcing is prescribed as a bottom boundary
681condition (see \S\ref{TRA_bbc}).
683The large eddy coefficient found in the mixed layer together with high
684vertical resolution implies that in the case of explicit time stepping
685(\np{ln\_zdfexp}=true) there would be too restrictive a constraint on
686the time step. Therefore, the default implicit time stepping is preferred
687for the vertical diffusion since it overcomes the stability constraint.
688A forward time differencing scheme (\np{ln\_zdfexp}=true) using a time
689splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative.
690Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both
691tracers and dynamics.
693% ================================================================
694% External Forcing
695% ================================================================
696\section{External Forcing}
699% -------------------------------------------------------------------------------------------------------------
700%        surface boundary condition
701% -------------------------------------------------------------------------------------------------------------
702\subsection   [Surface boundary condition (\textit{trasbc})]
703         {Surface boundary condition (\mdl{trasbc})}
706The surface boundary condition for tracers is implemented in a separate
707module (\mdl{trasbc}) instead of entering as a boundary condition on the vertical
708diffusion operator (as in the case of momentum). This has been found to
709enhance readability of the code. The two formulations are completely
710equivalent; the forcing terms in trasbc are the surface fluxes divided by
711the thickness of the top model layer.
713Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components
714($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer
715of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
716and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$,
717the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details).
718By doing this, the forcing formulation is the same for any tracer (including temperature and salinity).
720The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following
721forcing fields (used on tracers):
723$\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface
724(i.e. the difference between the total surface heat flux and the fraction of the short wave flux that
725penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with
726of the mass exchange with the atmosphere and lands.
728$\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...)
730$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)
731 and possibly with the sea-ice and ice-shelves.
733$\bullet$ \textit{rnf}, the mass flux associated with runoff
734(see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
736$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt,
737(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied).
739The surface boundary condition on temperature and salinity is applied as follows:
740\begin{equation} \label{Eq_tra_sbc}
742 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^& \\ 
743& F^S =\frac{ 1 }{\rho _\,      \left. e_{3t} \right|_{k=1} }  &\overline{ \textit{sfx} }^t   & \\   
744 \end{aligned}
746where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps
747($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the
748divergence of odd and even time step (see \S\ref{STP}).
750In the linear free surface case (\np{ln\_linssh}~=~\textit{true}),
751an additional term has to be added on both temperature and salinity.
752On temperature, this term remove the heat content associated with mass exchange
753that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that
754would have resulted from a change in the volume of the first level.
755The resulting surface boundary condition is applied as follows:
756\begin{equation} \label{Eq_tra_sbc_lin}
758 &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }   
759           &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^& \\ 
761& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} } 
762           &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1}  \right) }^t   & \\   
763 \end{aligned}
765Note that an exact conservation of heat and salt content is only achieved with non-linear free surface.
766In the linear free surface case, there is a small imbalance. The imbalance is larger
767than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.
768This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}).
770% -------------------------------------------------------------------------------------------------------------
771%        Solar Radiation Penetration
772% -------------------------------------------------------------------------------------------------------------
773\subsection   [Solar Radiation Penetration (\textit{traqsr})]
774         {Solar Radiation Penetration (\mdl{traqsr})}
780Options are defined through the  \ngn{namtra\_qsr} namelist variables.
781When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true),
782the solar radiation penetrates the top few tens of meters of the ocean. If it is not used
783(\np{ln\_flxqsr}=false) all the heat flux is absorbed in the first ocean level.
784Thus, in the former case a term is added to the time evolution equation of
785temperature \eqref{Eq_PE_tra_T} and the surface boundary condition is
786modified to take into account only the non-penetrative part of the surface
787heat flux:
788\begin{equation} \label{Eq_PE_qsr}
790\frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}   \\
791Q_{ns} &= Q_\text{Total} - Q_{sr}
794where $Q_{sr}$ is the penetrative part of the surface heat flux ($i.e.$ the shortwave radiation)
795and $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$).
796The additional term in \eqref{Eq_PE_qsr} is discretized as follows:
797\begin{equation} \label{Eq_tra_qsr}
798\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]
801The shortwave radiation,  $Q_{sr}$, consists of energy distributed across a wide spectral range.
802The ocean is strongly absorbing for wavelengths longer than 700~nm and these
803wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$ 
804that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified
805through namelist parameter \np{rn\_abs}).  It is assumed to penetrate the ocean
806with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$,
807of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namelist).
808For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy
809propagates to larger depths where it contributes to
810local heating.
811The way this second part of the solar energy penetrates into the ocean depends on
812which formulation is chosen. In the simple 2-waveband light penetration scheme  (\np{ln\_qsr\_2bd}=true)
813a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,
814leading to the following expression  \citep{Paulson1977}:
815\begin{equation} \label{Eq_traqsr_iradiance}
816I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right]
818where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 
819It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter.
820The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in
821Jerlov's (1968) classification (oligotrophic waters).
823Such assumptions have been shown to provide a very crude and simplistic
824representation of observed light penetration profiles (\cite{Morel_JGR88}, see also
825Fig.\ref{Fig_traqsr_irradiance}). Light absorption in the ocean depends on
826particle concentration and is spectrally selective. \cite{Morel_JGR88} has shown
827that an accurate representation of light penetration can be provided by a 61 waveband
828formulation. Unfortunately, such a model is very computationally expensive.
829Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this
830formulation in which visible light is split into three wavebands: blue (400-500 nm),
831green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependent
832attenuation coefficient is fitted to the coefficients computed from the full spectral model
833of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), assuming
834the same power-law relationship. As shown in Fig.\ref{Fig_traqsr_irradiance},
835this formulation, called RGB (Red-Green-Blue), reproduces quite closely
836the light penetration profiles predicted by the full spectal model, but with much greater
837computational efficiency. The 2-bands formulation does not reproduce the full model very well.
839The RGB formulation is used when \np{ln\_qsr\_rgb}=true. The RGB attenuation coefficients
840($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform
841chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 
842in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation:
845a constant 0.05 g.Chl/L value everywhere ;
847an observed time varying chlorophyll deduced from satellite surface ocean color measurement
848spread uniformly in the vertical direction ;
850same as previous case except that a vertical profile of chlorophyl is used.
851Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ;
853simulated time varying chlorophyll by TOP biogeochemical model.
854In this case, the RGB formulation is used to calculate both the phytoplankton
855light limitation in PISCES or LOBSTER and the oceanic heating rate.
857The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation
858is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}.
860When the $z$-coordinate is preferred to the $s$-coordinate, the depth of $w-$levels does
861not significantly vary with location. The level at which the light has been totally
862absorbed ($i.e.$ it is less than the computer precision) is computed once,
863and the trend associated with the penetration of the solar radiation is only added down to that level.
864Finally, note that when the ocean is shallow ($<$ 200~m), part of the
865solar radiation can reach the ocean floor. In this case, we have
866chosen that all remaining radiation is absorbed in the last ocean
867level ($i.e.$ $I$ is masked).
870\begin{figure}[!t]     \begin{center}
872\caption{    \label{Fig_traqsr_irradiance}
873Penetration profile of the downward solar irradiance calculated by four models.
874Two waveband chlorophyll-independent formulation (blue), a chlorophyll-dependent
875monochromatic formulation (green), 4 waveband RGB formulation (red),
87661 waveband Morel (1988) formulation (black) for a chlorophyll concentration of
877(a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. From \citet{Lengaigne_al_CD07}.}
878\end{center}   \end{figure}
881% -------------------------------------------------------------------------------------------------------------
882%        Bottom Boundary Condition
883% -------------------------------------------------------------------------------------------------------------
884\subsection   [Bottom Boundary Condition (\textit{trabbc})]
885         {Bottom Boundary Condition (\mdl{trabbc})}
891\begin{figure}[!t]     \begin{center}
893\caption{   \label{Fig_geothermal}
894Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
895It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.}
896\end{center}   \end{figure}
899Usually it is assumed that there is no exchange of heat or salt through
900the ocean bottom, $i.e.$ a no flux boundary condition is applied on active
901tracers at the bottom. This is the default option in \NEMO, and it is
902implemented using the masking technique. However, there is a
903non-zero heat flux across the seafloor that is associated with solid
904earth cooling. This flux is weak compared to surface fluxes (a mean
905global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), but it warms
906systematically the ocean and acts on the densest water masses.
907Taking this flux into account in a global ocean model increases
908the deepest overturning cell ($i.e.$ the one associated with the Antarctic
909Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OS09}.
911Options are defined through the  \ngn{namtra\_bbc} namelist variables.
912The presence of geothermal heating is controlled by setting the namelist
913parameter  \np{ln\_trabbc} to true. Then, when \np{nn\_geoflx} is set to 1,
914a constant geothermal heating is introduced whose value is given by the
915\np{nn\_geoflx\_cst}, which is also a namelist parameter.
916When  \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is
917introduced which is provided in the \ifile{geothermal\_heating} NetCDF file
918(Fig.\ref{Fig_geothermal}) \citep{Emile-Geay_Madec_OS09}.
920% ================================================================
921% Bottom Boundary Layer
922% ================================================================
923\section  [Bottom Boundary Layer (\mdl{trabbl} - \key{trabbl})]
924      {Bottom Boundary Layer (\mdl{trabbl} - \key{trabbl})}
930Options are defined through the  \ngn{nambbl} namelist variables.
931In a $z$-coordinate configuration, the bottom topography is represented by a
932series of discrete steps. This is not adequate to represent gravity driven
933downslope flows. Such flows arise either downstream of sills such as the Strait of
934Gibraltar or Denmark Strait, where dense water formed in marginal seas flows
935into a basin filled with less dense water, or along the continental slope when dense
936water masses are formed on a continental shelf. The amount of entrainment
937that occurs in these gravity plumes is critical in determining the density
938and volume flux of the densest waters of the ocean, such as Antarctic Bottom Water,
939or North Atlantic Deep Water. $z$-coordinate models tend to overestimate the
940entrainment, because the gravity flow is mixed vertically by convection
941as it goes ''downstairs'' following the step topography, sometimes over a thickness
942much larger than the thickness of the observed gravity plume. A similar problem
943occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly
944downstream of a sill \citep{Willebrand_al_PO01}, and the thickness
945of the plume is not resolved.
947The idea of the bottom boundary layer (BBL) parameterisation, first introduced by
948\citet{Beckmann_Doscher1997}, is to allow a direct communication between
949two adjacent bottom cells at different levels, whenever the densest water is
950located above the less dense water. The communication can be by a diffusive flux
951(diffusive BBL), an advective flux (advective BBL), or both. In the current
952implementation of the BBL, only the tracers are modified, not the velocities.
953Furthermore, it only connects ocean bottom cells, and therefore does not include
954all the improvements introduced by \citet{Campin_Goosse_Tel99}.
956% -------------------------------------------------------------------------------------------------------------
957%        Diffusive BBL
958% -------------------------------------------------------------------------------------------------------------
959\subsection{Diffusive Bottom Boundary layer (\np{nn\_bbl\_ldf}=1)}
962When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1),
963the diffusive flux between two adjacent cells at the ocean floor is given by
964\begin{equation} \label{Eq_tra_bbl_diff}
965{\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T
967with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells,
968and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997},
969the latter is prescribed with a spatial dependence, $i.e.$ in the conditional form
970\begin{equation} \label{Eq_tra_bbl_coef}
971A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
972 A_{bbl}  \quad \quad   \mbox{if}  \quad   \nabla_\sigma \rho  \cdot  \nabla H<0 \\ 
973 \\
974 0\quad \quad \;\,\mbox{otherwise} \\ 
975 \end{array}} \right.
977where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist
978parameter \np{rn\_ahtbbl} and usually set to a value much larger
979than the one used for lateral mixing in the open ocean. The constraint in \eqref{Eq_tra_bbl_coef} 
980implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of
981the slope, is larger than in the deeper ocean (see green arrow in Fig.\ref{Fig_bbl}).
982In practice, this constraint is applied separately in the two horizontal directions,
983and the density gradient in \eqref{Eq_tra_bbl_coef} is evaluated with the log gradient formulation:
984\begin{equation} \label{Eq_tra_bbl_Drho}
985   \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S
987where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$,
988$\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature,
989salinity and depth, respectively.
991% -------------------------------------------------------------------------------------------------------------
992%        Advective BBL
993% -------------------------------------------------------------------------------------------------------------
994\subsection   {Advective Bottom Boundary Layer  (\np{nn\_bbl\_adv}= 1 or 2)}
997\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following
998if this is not what is meant then "downwards sloping flow" is also a possibility"}
1001\begin{figure}[!t]   \begin{center}
1003\caption{   \label{Fig_bbl} 
1004Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is
1005activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$.
1006Red arrows indicate the additional overturning circulation due to the advective BBL.
1007The transport of the downslope flow is defined either as the transport of the bottom
1008ocean cell (black arrow), or as a function of the along slope density gradient.
1009The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$
1010ocean bottom cells.
1012\end{center}   \end{figure}
1016%!!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
1017%!!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
1018%!!        i.e. transport proportional to the along-slope density gradient
1020%%%gmcomment   :  this section has to be really written
1022When applying an advective BBL (\np{nn\_bbl\_adv} = 1 or 2), an overturning
1023circulation is added which connects two adjacent bottom grid-points only if dense
1024water overlies less dense water on the slope. The density difference causes dense
1025water to move down the slope.
1027\np{nn\_bbl\_adv} = 1 : the downslope velocity is chosen to be the Eulerian
1028ocean velocity just above the topographic step (see black arrow in Fig.\ref{Fig_bbl})
1029\citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection
1030is allowed only if dense water overlies less dense water on the slope ($i.e.$ 
1031$\nabla_\sigma \rho  \cdot  \nabla H<0$) and if the velocity is directed towards
1032greater depth ($i.e.$ $\vect{U}  \cdot  \nabla H>0$).
1034\np{nn\_bbl\_adv} = 2 : the downslope velocity is chosen to be proportional to $\Delta \rho$,
1035the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
1036The advection is allowed only  if dense water overlies less dense water on the slope ($i.e.$ 
1037$\nabla_\sigma \rho  \cdot  \nabla H<0$). For example, the resulting transport of the
1038downslope flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by the
1039following expression:
1040\begin{equation} \label{Eq_bbl_Utr}
1041 u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right)
1043where $\gamma$, expressed in seconds, is the coefficient of proportionality
1044provided as \np{rn\_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn} 
1045are the vertical index of the higher and lower cells, respectively.
1046The parameter $\gamma$ should take a different value for each bathymetric
1047step, but for simplicity, and because no direct estimation of this parameter is
1048available, a uniform value has been assumed. The possible values for $\gamma$ 
1049range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}
1051Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ 
1052using the upwind scheme. Such a diffusive advective scheme has been chosen
1053to mimic the entrainment between the downslope plume and the surrounding
1054water at intermediate depths. The entrainment is replaced by the vertical mixing
1055implicit in the advection scheme. Let us consider as an example the
1056case displayed in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is
1057larger than the one at level $(i,kdwn)$. The advective BBL scheme
1058modifies the tracer time tendency of the ocean cells near the
1059topographic step by the downslope flow \eqref{Eq_bbl_dw},
1060the horizontal \eqref{Eq_bbl_hor}  and the upward \eqref{Eq_bbl_up} 
1061return flows as follows:
1063\partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw}
1064                                     +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right\label{Eq_bbl_dw} \\
1066\partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 
1067               + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{Eq_bbl_hor} \\
1069\intertext{and for $k =kdw-1,\;..., \; kup$ :} 
1071\partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k}
1072               + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{Eq_bbl_up}
1074where $b_t$ is the $T$-cell volume.
1076Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in
1077the model outputs. It has to be used to compute the effective velocity
1078as well as the effective overturning circulation.
1080% ================================================================
1081% Tracer damping
1082% ================================================================
1083\section  [Tracer damping (\textit{tradmp})]
1084      {Tracer damping (\mdl{tradmp})}
1090In some applications it can be useful to add a Newtonian damping term
1091into the temperature and salinity equations:
1092\begin{equation} \label{Eq_tra_dmp}
1094 \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right\\
1095 \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right)
1096 \end{split}
1097 \end{equation} 
1098where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ 
1099are given temperature and salinity fields (usually a climatology).
1100Options are defined through the  \ngn{namtra\_dmp} namelist variables.
1101The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true.
1102It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true
1103in \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are
1104correctly set  ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read
1105using \mdl{fldread}, see \S\ref{SBC_fldread}).
1106The 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.
1108The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} 
1109the specification of the boundary conditions along artificial walls of a
1110limited domain basin and \textit{(b)} the computation of the velocity
1111field associated with a given $T$-$S$ field (for example to build the
1112initial state of a prognostic simulation, or to use the resulting velocity
1113field for a passive tracer study). The first case applies to regional
1114models that have artificial walls instead of open boundaries.
1115In the vicinity of these walls, $\gamma$ takes large values (equivalent to
1116a time scale of a few days) whereas it is zero in the interior of the
1117model domain. The second case corresponds to the use of the robust
1118diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity
1119field consistent with the model dynamics whilst having a $T$, $S$ field
1120close to a given climatological field ($T_o$, $S_o$).
1122The robust diagnostic method is very efficient in preventing temperature
1123drift in intermediate waters but it produces artificial sources of heat and salt
1124within the ocean. It also has undesirable effects on the ocean convection.
1125It tends to prevent deep convection and subsequent deep-water formation,
1126by stabilising the water column too much.
1128The 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}.
1130\subsection[DMP\_TOOLS]{Generating using DMP\_TOOLS}
1132DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$.
1133Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled
1134and run on the same machine as the NEMO model. A mesh\ file for the model configuration is required as an input.
1135This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1.
1136The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work.
1137The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
1143\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.
1145The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations.
1146\np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain.
1147\np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea
1148for the ORCA4, ORCA2 and ORCA05 configurations.
1149If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as
1150a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference
1151configurations with previous model versions.
1152\np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines.
1153This option only has an effect if \np{ln\_full\_field} is true.
1154\np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer.
1155Finally \np{ln\_custom} specifies that the custom module will be called.
1156This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
1158The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}.
1159Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to
1160the full values of a 10$^{\circ}$ latitud band.
1161This is often used because of the short adjustment time scale in the equatorial region
1162\citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a
1163hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}
1165% ================================================================
1166% Tracer time evolution
1167% ================================================================
1168\section  [Tracer time evolution (\textit{tranxt})]
1169      {Tracer time evolution (\mdl{tranxt})}
1175Options are defined through the  \ngn{namdom} namelist variables.
1176The general framework for tracer time stepping is a modified leap-frog scheme
1177\citep{Leclair_Madec_OM09}, $i.e.$ a three level centred time scheme associated
1178with a Asselin time filter (cf. \S\ref{STP_mLF}):
1179\begin{equation} \label{Eq_tra_nxt}
1181(e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt  \,e_{3t}^t\ \text{RHS}^t & \\
1183(e_{3t}T)_f^\;\ \quad &= (e_{3t}T)^t \;\quad 
1184                                    &+\gamma \,\left[ {(e_{3t}T)_f^{t-\rdt} -2(e_{3t}T)^t+(e_{3t}T)^{t+\rdt}} \right] &  \\
1185                                 & &- \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]  &                     
1188where RHS is the right hand side of the temperature equation,
1189the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient,
1190and $S$ is the total forcing applied on $T$ ($i.e.$ fluxes plus content in mass exchanges).
1191$\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter).
1192Its default value is \np{rn\_atfp}=$10^{-3}$. Note that the forcing correction term in the filter
1193is not applied in linear free surface (\jp{lk\_vvl}=false) (see \S\ref{TRA_sbc}.
1194Not also that in constant volume case, the time stepping is performed on $T$,
1195not on its content, $e_{3t}T$.
1197When the vertical mixing is solved implicitly, the update of the \textit{next} tracer
1198fields is done in module \mdl{trazdf}. In this case only the swapping of arrays
1199and the Asselin filtering is done in the \mdl{tranxt} module.
1201In order to prepare for the computation of the \textit{next} time step,
1202a swap of tracer arrays is performed: $T^{t-\rdt} = T^t$ and $T^t = T_f$.
1204% ================================================================
1205% Equation of State (eosbn2)
1206% ================================================================
1207\section  [Equation of State (\textit{eosbn2}) ]
1208      {Equation of State (\mdl{eosbn2}) }
1214% -------------------------------------------------------------------------------------------------------------
1215%        Equation of State
1216% -------------------------------------------------------------------------------------------------------------
1217\subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)}
1220The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship
1221linking seawater density, $\rho$, to a number of state variables,
1222most typically temperature, salinity and pressure.
1223Because density gradients control the pressure gradient force through the hydrostatic balance,
1224the equation of state provides a fundamental bridge between the distribution of active tracers
1225and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular
1226influencing the circulation through determination of the static stability below the mixed layer,
1227thus controlling rates of exchange between the atmosphere  and the ocean interior \citep{Roquet_JPO2015}.
1228Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983})
1229or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real
1230ocean circulation is attempted \citep{Roquet_JPO2015}.
1231The use of TEOS-10 is highly recommended because
1232\textit{(i)} it is the new official EOS,
1233\textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and
1234\textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature
1235and practical salinity for EOS-980, both variables being more suitable for use as model variables
1236\citep{TEOS10, Graham_McDougall_JPO13}.
1237EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility.
1238For process studies, it is often convenient to use an approximation of the EOS. To that purposed,
1239a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available.
1241In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$,
1242is computed, with $\rho_o$ a reference density. Called \textit{rau0} 
1243in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$.
1244This is a sensible choice for the reference density used in a Boussinesq ocean
1245climate model, as, with the exception of only a small percentage of the ocean,
1246density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
1248Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos} 
1249which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS).
1252\item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 
1253The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
1254but it is optimized for a boussinesq fluid and the polynomial expressions have simpler
1255and more computationally efficient expressions for their derived quantities
1256which make them more adapted for use in ocean models.
1257Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10
1258rational function approximation for hydrographic data analysis  \citep{TEOS10}.
1259A key point is that conservative state variables are used:
1260Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: $\degres C$, notation: $\Theta$).
1261The pressure in decibars is approximated by the depth in meters.
1262With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to
1263$C_p=3991.86795711963~J\,Kg^{-1}\,\degres K^{-1}$, according to \citet{TEOS10}.
1265Choosing polyTEOS10-bsq implies that the state variables used by the model are
1266$\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as
1267\textit{Conservative} Temperature and \textit{Absolute} Salinity.
1268In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST
1269prior to either computing the air-sea and ice-sea fluxes (forced mode)
1270or sending the SST field to the atmosphere (coupled mode).
1272\item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used.
1273It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized
1274to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80
1275and the ocean model are:
1276the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $\degres C$, notation: $\theta$).
1277The pressure in decibars is approximated by the depth in meters. 
1278With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature,
1279salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to
1280have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant
1281value, the TEOS10 value.
1283\item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,
1284the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.)
1285(see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both
1286cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS
1287in theoretical studies \citep{Roquet_JPO2015}.
1288With such an equation of state there is no longer a distinction between
1289\textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute} 
1290and \textit{practical} salinity.
1291S-EOS takes the following expression:
1292\begin{equation} \label{Eq_tra_S-EOS}
1294  d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_\\
1295                                & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_\\
1296                                & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\
1297  with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3
1300where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}.
1301In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing
1302the associated coefficients.
1303Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS.
1304setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS.
1305Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S.
1312\begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|}
1314coeff.   & computer name   & S-EOS     &  description                      \\ \hline
1315$a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff.    \\ \hline
1316$b_0$       & \np{rn\_b0}     & 7.6554 $10^{-1}$ &  linear haline  expansion coeff.    \\ \hline
1317$\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$          \\ \hline
1318$\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$       \\ \hline
1319$\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$       \\ \hline
1320$\mu_1$     & \np{rn\_mu1}    & 1.4970 $10^{-4}$ &  thermobaric coeff. in T         \\ \hline
1321$\mu_2$     & \np{rn\_mu2}    & 1.1090 $10^{-5}$ &  thermobaric coeff. in S            \\ \hline
1323\caption{ \label{Tab_SEOS}
1324Standard value of S-EOS coefficients. }
1330% -------------------------------------------------------------------------------------------------------------
1331%        Brunt-V\"{a}is\"{a}l\"{a} Frequency
1332% -------------------------------------------------------------------------------------------------------------
1333\subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
1336An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a}
1337 frequency) is of paramount importance as determine the ocean stratification and
1338 is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent
1339 vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing
1340 parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure
1341 (pressure in decibar being approximated by the depth in meters). The expression for $N^2$ 
1342 is given by:
1343\begin{equation} \label{Eq_tra_bn2}
1344N^2 = \frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
1346where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS,
1347and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.
1348The coefficients are a polynomial function of temperature, salinity and depth which expression
1349depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran} 
1350function that can be found in \mdl{eosbn2}.
1352% -------------------------------------------------------------------------------------------------------------
1353%        Freezing Point of Seawater
1354% -------------------------------------------------------------------------------------------------------------
1355\subsection   [Freezing Point of Seawater]
1356         {Freezing Point of Seawater}
1359The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}:
1360\begin{equation} \label{Eq_tra_eos_fzp}
1361   \begin{split}
1362T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} 
1363                       -  2.154996 \;10^{-4} \,\right) \ S    \\
1364               - 7.53\,10^{-3} \ \ p
1365   \end{split}
1368\eqref{Eq_tra_eos_fzp} is only used to compute the potential freezing point of
1369sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent
1370terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing
1371point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found
1372in \mdl{eosbn2}
1375% -------------------------------------------------------------------------------------------------------------
1376%        Potential Energy     
1377% -------------------------------------------------------------------------------------------------------------
1378%\subsection{Potential Energy anomalies}
1381%    =====>>>>> TO BE written
1385% ================================================================
1386% Horizontal Derivative in zps-coordinate
1387% ================================================================
1388\section  [Horizontal Derivative in \textit{zps}-coordinate (\textit{zpshde})]
1389      {Horizontal Derivative in \textit{zps}-coordinate (\mdl{zpshde})}
1392\gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators,
1393                   I've changed "derivative" to "difference" and "mean" to "average"}
1395With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general,
1396tracers in horizontally adjacent cells live at different depths.
1397Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module)
1398and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module).
1399The partial cell properties at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom.
1400So, only the bottom interpolation is explained below.
1402Before taking horizontal gradients between the tracers next to the bottom, a linear
1403interpolation in the vertical is used to approximate the deeper tracer as if it actually
1404lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}).
1405For example, for temperature in the $i$-direction the needed interpolated
1406temperature, $\widetilde{T}$, is:
1409\begin{figure}[!p]    \begin{center}
1411\caption{   \label{Fig_Partial_step_scheme} 
1412Discretisation of the horizontal difference and average of tracers in the $z$-partial
1413step coordinate (\np{ln\_zps}=true) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$.
1414A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value
1415at the depth of the shallower tracer point of the two adjacent bottom $T$-points.
1416The horizontal difference is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ 
1417and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$}
1418\end{center}   \end{figure}
1421\widetilde{T}= \left\{  \begin{aligned} 
1422&T^{\,i+1}      -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta _k T^{i+1} 
1423                        && \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$   }  \\
1424                              \\
1425&T^{\,i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i       }\;\delta _k T^{i+1}
1426                        && \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
1427            \end{aligned}   \right.
1429and the resulting forms for the horizontal difference and the horizontal average
1430value of $T$ at a $U$-point are:
1431\begin{equation} \label{Eq_zps_hde}
1433 \delta _{i+1/2} T=  \begin{cases}
1434\ \ \ \widetilde {T}\quad\ -T^i     & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\
1435                              \\
1436\ \ \ T^{\,i+1}-\widetilde{T}    & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
1437                  \end{cases}     \\
1439\overline {T}^{\,i+1/2} \ =   \begin{cases}
1440( \widetilde {T}\ \ \;\,-T^{\,i})    / 2  & \;\ \ \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\
1441                              \\
1442( T^{\,i+1}-\widetilde{T} ) / 2     & \;\ \ \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
1443            \end{cases}
1447The computation of horizontal derivative of tracers as well as of density is
1448performed once for all at each time step in \mdl{zpshde} module and stored
1449in shared arrays to be used when needed. It has to be emphasized that the
1450procedure used to compute the interpolated density, $\widetilde{\rho}$, is not
1451the same as that used for $T$ and $S$. Instead of forming a linear approximation
1452of density, we compute $\widetilde{\rho }$ from the interpolated values of $T$ 
1453and $S$, and the pressure at a $u$-point (in the equation of state pressure is
1454approximated by depth, see \S\ref{TRA_eos} ) :
1455\begin{equation} \label{Eq_zps_hde_rho}
1456\widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u })
1457\quad \text{where }\  z_u = \min \left( {z_T^{i+1} ,z_T^i } \right)
1460This is a much better approximation as the variation of $\rho$ with depth (and
1461thus pressure) is highly non-linear with a true equation of state and thus is badly
1462approximated with a linear interpolation. This approximation is used to compute
1463both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral
1464surfaces (\S\ref{LDF_slp})
1466Note that in almost all the advection schemes presented in this Chapter, both
1467averaging and differencing operators appear. Yet \eqref{Eq_zps_hde} has not
1468been used in these schemes: in contrast to diffusion and pressure gradient
1469computations, no correction for partial steps is applied for advection. The main
1470motivation is to preserve the domain averaged mean variance of the advected
1471field when using the $2^{nd}$ order centred scheme. Sensitivity of the advection
1472schemes to the way horizontal averages are performed in the vicinity of partial
1473cells should be further investigated in the near future.
1475\gmcomment{gm :   this last remark has to be done}
Note: See TracBrowser for help on using the repository browser.