source: branches/UKMO/v3_6_extra_CMIP6_diagnostics/DOC/TexFiles/Chapters/Chap_TRA.tex @ 7492

Last change on this file since 7492 was 7492, checked in by timgraham, 5 years ago

Final merge of head of v3.6 STABLE branch into this branch

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