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

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

Reinstate svn keywords

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