source: branches/2015/nemo_v3_6_STABLE/DOC/TexFiles/Chapters/Chap_TRA.tex @ 6039

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

#1629: DOC of v3.6_stable : update EOS presentation

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