source: trunk/DOC/BETA/Chapters/Chap_TRA.tex @ 781

Last change on this file since 781 was 781, checked in by smasson, 13 years ago

doc update, see ticket:1

  • Property svn:executable set to *
File size: 64.4 KB
1% ================================================================
2% Chapter 1 Ñ Ocean Tracers (TRA)
3% ================================================================
4\chapter{Ocean Tracers (TRA)}
8% missing/update
9% traqsr: need to coordinate with SBC module
10% trabbl : advective case to be discussed
11%        diffusive case : add : only the bottom ocean cell is concerned
12%        ==> addfigure on bbl
15$\ $\newline    % force a new ligne
17Using the representation described in Chap.~\ref{DOM}, several semi-discrete space
18forms of the tracer equations are available depending on the vertical
19coordinate used and on the physics used. In all the equations presented
20here, the masking has been omitted for simplicity. One must be aware that
21all the quantities are masked fields and that each time a mean or difference
22operator is used, the resulting field is multiplied by a mask.
24The two active tracers are potential temperature and salinity. Their prognostic equation can be summarized as follows:
26\text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC}
27                   \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP})
30NXT stands for next, referring to the time-stepping. From left to right, the terms on the rhs of the tracer equations are the advection (ADV), the lateral diffusion (LDF), the vertical diffusion (ZDF), the contributions from the external forcings (SBC: Surface Boundary Condition, QSR: Solar Radiation penetration, and BBC: Bottom Boundary Condition), the contribution from the bottom boundary Layer (BBL) parametrisation, and an internal damping (DMP) term. The last four have been put inside brackets as they are optional. The external forcings and parameterizations require complex inputs and calculations (bulk formulae, estimation of mixing coefficients) that are carried out in modules of the SBC, LDF and ZDF categories and described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively. Note that \mdl{tranpc}, the non-penetrative convection module,  although (temporarily) located in the NEMO/OPA//TRA directory, is described with the model vertical physics (ZDF).
32In the present chapter we also describe the diagnostic equations used to compute the sea-water properties (density, Brunt-Vais\"{a}l\"{a} frequency, specific heat and freezing point) although the associated modules ($i.e.$ \mdl{eosbn2}, \mdl{ocfzpt} and \mdl{phycst}) are (temporarily) located in the NEMO/OPA directory.
34The different options available to the user are managed by namelist logical or CPP keys. For each equation term ttt, the namelist logicals are \textit{ln\_trattt\_xxx}, where \textit{xxx} is a 3 or 4 letter acronym accounting for each optional scheme. The CPP key (when it exists) is \textbf{key\_trattt}. The corresponding code can be found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.
36The user has the option of extracting each tendency term on the rhs of the tracer equation (\key {trdtra} defined), as described in Chap.~\ref{MISC}.
38% ================================================================
39% Tracer Advection
40% ================================================================
41\section{Tracer Advection (\mdl{traadv})}
47The advection tendency in flux form is the divergence of the advective
48fluxes. Its discrete expression is given by :
49\begin{equation} \label{Eq_tra_adv}
50ADV_\tau =-\frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}e_{3T} }\left(
51{\;\delta _i \left[ {e_{2u} {\kern 1pt}e_{3u} {\kern 1pt}\;u\;\tau _u } 
52\right]+\delta _j \left[ {e_{1v} {\kern 1pt}e_{3v} {\kern 1pt}v\;\tau _v } 
53\right]\;} \right)-\frac{1}{\mathop e\nolimits_{3T} }\delta _k \left[
54{w\;\tau _w } \right]
56which, in pure z-coordinate (\key{zco} defined), reduces to :
57\begin{equation} \label{Eq_tra_adv_zco}
58ADV_\tau =-\frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}}\left( {\;\delta _i
59\left[ {e_{2u} {\kern 1pt}{\kern 1pt}\;u\;\tau _u } \right]+\delta _j \left[
60{e_{1v} {\kern 1pt}v\;\tau _v } \right]\;} \right)-\frac{1}{\mathop 
61e\nolimits_{3T} }\delta _k \left[ {w\;\tau _w } \right]
63as the vertical scale factors are function of $k$ only, and thus $e_{3u}
64=e_{3v} =e_{3T} $.
66The flux form requires implicitly the use of the continuity equation:
67$\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
68using $\nabla \cdot \vect{U}=0)$ or $\partial _t e_3 +\nabla \cdot \vect{U}=0$ in variable volume case ($i.e.$ \key{vvl} defined). Therefore it is of
69paramount importance to design the discrete analogue of the advection
70tendency so that it is consistent with the continuity equation in order to
71enforce conservation properties of the continuous equations. In other words,
72by substituting $\tau$ by 1 in (\ref{Eq_tra_adv}) we recover discrete form of the
73continuity equation which is used to calculate the vertical velocity.
75\begin{figure}[!t] \label{Fig_adv_scheme}  \begin{center}
77\caption{Schematic representation of some ways used to evaluate the tracer value at $u$-point and the amount of tracer exchanged between two neighbouring grid points. Upsteam biased scheme (ups): the upstream value is used and the black area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation is used and black + dark grey areas is exchanged. Monotonic upstream scheme for conservative laws (muscl):  a parabolic interpolation is used and black + dark grey + grey areas are exchanged. Second order scheme (cen2): the mean value is used and black + dark grey + grey + light grey areas are exchanged. Note that this illustration does not include the flux limiter used in ppm and muscl schemes.}
78\end{center}   \end{figure}
80The advection schemes used in OPA differ by the choice made in space and
81time interpolation to define the value of the tracer at the velocity points (\ref{Fig_adv_scheme}).
82Along solid lateral and bottom boundaries a zero tracer flux is naturally
83specified, since the normal velocity is zero there. At the sea surface the
84boundary condition depends on the type of sea surface chosen: (1) in
85rigid-lid formulation, $w=0$ at the surface, so the advective fluxes through the
86surface is zero ; (2) in non-linear free surface (variable volume case,
87\key{vvl} defined), convergence/divergence in the first ocean level moves
88up/down the free surface: there is no tracer advection through it so that
89the advective fluxes through the surface is also zero ; (3) in the linear
90free surface, the first level thickness is constant in time. The vertical
91boundary condition is applied at the fixed surface $z=0$ rather than on the
92moving surface $z=\eta$. There is a non-zero advective flux which is set for all
93advection schemes as the product of surface velocity (at $z=0$) by the first level
94tracer value: $\left. {\tau _w } \right|_{k=1 \mathord{\left/ {\vphantom {1 
952}} \right. \kern-\nulldelimiterspace} 2} =T_{k=1} $. This boundary
96condition retains local conservation of tracer. Strict global conservation
97is not possible in linear free surface but is achieved to a good
98approximation since the non-conservative term is the product of the time
99derivative of the tracer and the free surface height, two quantities that are not correlated (see \S\ref{PE_free_surface}, and also \citet{Roullet2000,Griffies2001,Campin2004}).
101The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco}) is the centred (\textit{now}) \textit{eulerian} ocean velocity (see \S\ref{DYN}). Nevertheless, when advective bottom boundary layer (\textit{bbl}) and/or eddy induced velocity (\textit{eiv}) parameterisations are used it is the \textit{now} \textit{effective} velocity (i.e. the sum of the eulerian, the bbl and/or the eiv velocities) which is used.
103The choice of an advection scheme is made in the \np{nam\_traadv} namelist, by
104setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The
105corresponding code can be found in \textit{traadv\_xxx.F90} module, where
106\textit{xxx} is a 3 or 4 letter acronym accounting for each scheme. Details
107of the advection schemes are given below. The choice of an advection scheme
108is a complex matter which depends on the model physics, model resolution,
109type of tracer, as well as the issue of numerical cost.
111Note that (1) cen2, cen4 and TVD schemes require an explicit diffusion
112operator while the other schemes are diffusive enough so that they do not
113require additional diffusion ; (2) cen2, cen4, MUSCL2, and UBS are not
114\textit{positive} schemes, meaning false extrema are permitted. It is not recommended to use
115them on passive tracers ; (3) It is highly recommended to use the same
116advection-diffusion scheme on both active and passive tracers. In
117particular, if a source or sink of a passive tracer depends on a active one,
118the difference of treatment of active and passive tracers can create very
119nice-looking frontal structures that are pure numerical artefacts.
122% -------------------------------------------------------------------------------------------------------------
123%        2nd order centred scheme 
124% -------------------------------------------------------------------------------------------------------------
125\subsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=T)}
128In the centred second order formulation, the tracer at velocity points is
129evaluated as the mean of the two neighbouring $T$-points. For example, in the
130$i$-direction :
131\begin{equation} \label{Eq_tra_adv_cen2}
132\tau _u^{cen2} =\overline T ^{i+1/2}
135The scheme is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to produce a sensible solution. The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
137Note that using cen2 scheme, the overall tracer advection is of second order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) have this order of accuracy.
139% -------------------------------------------------------------------------------------------------------------
140%        4nd order centred scheme 
141% -------------------------------------------------------------------------------------------------------------
142\subsection{$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=T)}
145In the $4^{th}$ order formulation (to be implemented), tracer is evaluated
146at velocity points as the $4^{th}$ order interpolation of $T$, and thus use the
147four neighbouring $T$-points. For example, in the $i$-direction:
148\begin{equation} \label{Eq_tra_adv_cen4}
149\tau _u^{cen4} 
150=\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
153Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme
154but a $4^{th}$ order evaluation of advective fluxes since the divergence of
155advective fluxes, (\ref{Eq_tra_adv}), is kept at $2^{nd}$ order. The ``$4^{th}$ order
156scheme'' denomination used in oceanographic literature is usually associated
157with the scheme presented here. Introducing a \textit{true} $4^{th}$ order advection
158scheme is feasible but, for consistency reasons, it requires changes in the
159discretisation of the tracer advection together with changes in both the
160continuity equation and the momentum advection.
162A direct consequence of the pseudo-fourth order nature of the scheme is that
163it is not non-diffusive, i.e. the global variance of a tracer is not
164preserved through \textit{cen4}. Furthermore, it must be used in conjunction with an
165explicit diffusion operator to produce a sensible solution. The
166time-stepping is also performed using a leapfrog scheme in conjunction with
167an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
169At $T$-grid cell abutted to a boundary (coastline, bottom and surface), an
170additional hypothesis must be made to evaluate $\tau _u^{cen4}$. This
171hypothesis usually reduces the order of the scheme. Here we choose to set
172the gradient of $T$ across the boundary to zero. Alternative conditions can be
173specified such as the reduction to a second order scheme for near boundary
174grid point.
176% -------------------------------------------------------------------------------------------------------------
177%        TVD scheme 
178% -------------------------------------------------------------------------------------------------------------
179\subsection{Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=T)}
182In the Total Variance Dissipation (TVD) formulation, the tracer at velocity
183points is evaluated as a combination of upstream and centred scheme. For
184example, in the $i$-direction :
185\begin{equation} \label{Eq_tra_adv_tvd}
187\tau _u^{ups}&= \begin{cases}
188               T_{i+1}  & \text{if $\ u_{i+1/2} <     0$} \hfill \\
189               T_i         & \text{if $\ u_{i+1/2} \geq 0$} \hfill \\
190              \end{cases}     \\
192\tau _u^{tvd}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen2} -\tau _u^{ups} } \right)
195where $c_u$ is a flux limiter function taking values between 0 and 1. There
196exists many ways to define $c_u$., each correcponding to a different total
197variance decreasing scheme. The one chosen in OPA is described in \citet{Zalesak1979}. $c_u$ only departs from $1$ when the advective term produces a local
198extremum in the tracer field. The resulting scheme is quite expensive but
199\emph{positive}. It can be used on both active and passive tracers. This scheme is tested
200and compared with MUSCL and the MPDATA scheme in \citet{Levy2001}; note that
201in this paper it is referred to as "FCT" (Flux corrected transport)
202rather than TVD.
204For stability reasons in (\ref{Eq_tra_adv_tvd}) $\tau _u^{cen2}$ is evaluated using the
205\textit{now} velocity (leap-frog environment: centred in time) while $\tau _u^{ups}$ is
206evaluated using the \textit{before} velocity (diffusive part: forward in time).
208% -------------------------------------------------------------------------------------------------------------
209%        MUSCL scheme 
210% -------------------------------------------------------------------------------------------------------------
211\subsection[MUSCL scheme  (\np{ln\_traadv\_muscl}=T)]
212   {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_muscl}=T)}
215The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been
216implemented by \citet{Levy2001}. In its formulation, the tracer at
217velocity points is evaluated assuming a linear tracer variation between two
218$T$-points (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
219\begin{equation} \label{Eq_tra_adv_muscl}
220   \tau _u^{mus} = \left\{      \begin{aligned}
221         &\tau _&+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\Delta t}{e_{1u}} \right)
222         &\ \widetilde{\partial _i \tau}  & \quad \text{if }\;u_{i+1/2} \geqslant 0      \\
223         &\tau _{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\Delta t}{e_{1u} } \right)
224         &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0
225   \end{aligned}    \right.
227where $\widetilde{\partial _i \tau}$ is the slope of the tracer
228on which a limitation is imposed to ensure the \textit{positive} character of the scheme.
230The time stepping is performed using a forward scheme, that is the \textit{before} tracer
231field is used to evaluate $\tau _u^{mus}$.
233For an ocean grid point abutted to land and where the ocean velocity is
234toward land, two choices are available: use of an upstream flux
235(\np{ln\_traadv\_muscl}=T) or use of second order flux
236(\np{ln\_traadv\_muscl2}=T). Note that the latter choice does not insure the
237\textit{positive} character of the scheme. Only the former can be used on both active and
238passive tracers.
240% -------------------------------------------------------------------------------------------------------------
241%        UBS scheme 
242% -------------------------------------------------------------------------------------------------------------
243\subsection{Upstream Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=T)}
246The UBS advection scheme is an upstream biased third order scheme based on
247an upstream-biased parabolic interpolation. It is also known as Cell
248Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective
249Kinematics). For example, in the $i$-direction :
250\begin{equation} \label{Eq_tra_adv_ubs}
251   \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{     
252   \begin{aligned}
253         &\tau"_i          & \quad \text{if }\ u_{i+1/2} \geqslant 0      \\
254         &\tau"_{i+1}   & \quad \text{if }\ u_{i+1/2}       <       0
255   \end{aligned}    \right.
257where $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$.
259This results in a dissipatively dominant (i.e. hyper-diffusive) truncation
260error \citep{Sacha2005}. The overall performance of the
261advection scheme is similar to that reported in \cite{Farrow1995}.
262It is a relatively good compromise between accuracy and smoothness. It is
263not a \emph{positive} scheme meaning false extrema are permitted but the
264amplitude of such are significantly reduced over the centred second order
265method. Nevertheless it is not recommended to apply it to a passive tracer
266that requires positivity.
268The intrinsic diffusion of UBS makes its use risky in the vertical direction
269where the control of artificial diapycnal fluxes is of paramount importance.
270It has therefore been preferred to evaluate the vertical flux using the TVD
271scheme when \np{ln\_traadv\_ubs}=T.
273For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds to a
274second order centred scheme is evaluated using the \textit{now} velocity (centred in
275time) while the second term which is the diffusive part of the scheme, is
276evaluated using the \textit{before} velocity (forward in time. This is discussed by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK
277schemes only differ by one coefficient. Substituting 1/6 with 1/8 in
278(\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}. This
279option is not available through a namelist parameter, since the 1/6
280coefficient is hard coded. Nevertheless it is quite easy to make the
281substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme
283NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can
284be controlled by vertical advection (not vertical diffusion which is usually
285solved using an implicit scheme). Computer time can be saved by using a
286time-splitting technique on vertical advection. This possibility have been
287implemented and validated in ORCA05-L301. It is not currently offered in the
288current reference version.
290NB 2 : In a forthcoming release four options will be proposed for the
291vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be
292evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme , or  \textit{(b)} a TVD
293scheme, or  \textit{(c)} an interpolation based on conservative parabolic splines
294following \citet{Sacha2005} implementation of UBS in ROMS,
295or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an
296eight-order accurate conventional scheme.
298NB 3 : It is straight forward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
299\begin{equation} \label{Eq_tra_adv_ubs2}
300\tau _u^{ubs} = \left\{  \begin{aligned}
301   & \tau _u^{cen4} + \frac{1}{12} \tau"_i      & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
302   & \tau _u^{cen4} - \frac{1}{12} \tau"_{i+1}  & \quad \text{if }\ u_{i+1/2}       <       0
303                   \end{aligned}    \right.
305or equivalently
306\begin{equation} \label{Eq_tra_adv_ubs2b}
307u_{i+1/2} \ \tau _u^{ubs} 
308=u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta _i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}
309- \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
311\eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that the UBS scheme is based on the fourth order scheme to which is added an upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order part as stated above using \eqref{Eq_tra_adv_ubs} and also as it is coded in NEMO v2.3. Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient with is simply proportional to the velocity: $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that the current version of NEMO uses \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_tra_adv_ubs2}.
314% -------------------------------------------------------------------------------------------------------------
315%        QCK scheme 
316% -------------------------------------------------------------------------------------------------------------
317\subsection{QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=T)}
320The Quadratic Upstream Interpolation for Convective Kinematics with
321Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} is
322the third order Godunov scheme. It is associated with ULTIMATE QUICKEST
323limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray
326The resulting scheme is quite expensive but \emph{positive}. It can be used on both active
327and passive tracers.
329% -------------------------------------------------------------------------------------------------------------
330%        PPM scheme 
331% -------------------------------------------------------------------------------------------------------------
332\subsection{Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=T)}
335The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984)
336is based on a quadradic piecewise rebuilding. As QCK scheme, it is associated with
337ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the current reference
340% ================================================================
341% Tracer Lateral Diffusion
342% ================================================================
343\section{Tracer Lateral Diffusion (\mdl{traldf})}
349The options available for lateral diffusion are laplacian (rotated or not)
350or biharmonic operators, that latter being more scale-selective (more
351diffusive at small scales). The specification of eddy diffusivity
352coefficients (either constant, variable in space and time) as well as the
353computation of the slope along which the operators act are performed in
354\mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}. The lateral diffusion of tracers is evaluated using a forward scheme,
355i.e. the tracers appearing in its expression are the \textit{before} tracers in time,
356except for the pure vertical component that appears when a tensor of
357rotation is used. This latter term is solved implicitly together with the
358vertical diffusion term (see \S\ref{DOM_nxt})
360% -------------------------------------------------------------------------------------------------------------
361%        Iso-level laplacian operator
362% -------------------------------------------------------------------------------------------------------------
363\subsection{Iso-level laplacian operator (\mdl{traldf\_lap}, \np{ln\_traldf\_lap}) }
366A laplacian diffusive operator (i.e. a harmonic operator) acting along the model surfaces is given by:
367\begin{equation} \label{Eq_tra_ldf_lap}
369D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta _i
370\left[ {A_u^{lT} \left( {\frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} 
371\left[ T \right]} \right)} \right]} \right.
373&\ \left. {+\; \delta _j \left[
374{A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta _{j+1/2} \left[ T
375\right]} \right)} \right]\quad } \right]
379This lateral operator is a \emph{horizontal} one (i.e. acting along geopotential surfaces) in
380$z$-coordinate with or without partial step, but it is simply an iso-level operator in $s$-coordinate.
381It is thus used when, in addition to \np{ln\_traldf\_lap}=T, we have \np{ln\_traldf\_level}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=F. In both cases, it significantly contributes to diapycnal mixing. It is therefore not recommended to use it.
383\textit{Notes} : In pure z-coordinate (\key{zco} defined), $e_{3u}=e_{3v}=e_{3T}$, so that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}).
385\textit{Notes} : In partial step $z$-coordinate (\np{ln\_zps}=T), tracers in horizontally adjacent cells are located at different depths in vicinity of the bottom. In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level require a specific treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
387% -------------------------------------------------------------------------------------------------------------
388%        Rotated laplacian operator
389% -------------------------------------------------------------------------------------------------------------
390\subsection{Rotated laplacian operator (\mdl{traldf\_iso}, \np{ln\_traldf\_lap})}
393The general form of the second order lateral tracer subgrid scale physics (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates:
395\begin{equation} \label{Eq_tra_ldf_iso}
397 D_T^{lT} =& \frac{1}{e_{1T}\,e_{2T}\,e_{3T} }
398 \\
399& \left\{ {\delta _i \left[ {A_u^{lT}  \left(
400    {\frac{e_{2u} \; e_{3u} }{e_{1u} } \,\delta _{i+1/2}[T]
401   -e_{2u} \; r_{1u} \,\overline{\overline {\delta _{k+1/2}[T]}}^{\,i+1/2,k}}
402 \right)} \right]} \right.
404& +\delta 
405_j \left[ {A_v^{lT} \left( {\frac{e_{1v}\,e_{3v} }{e_{2v} 
406}\,\delta _{j+1/2} \left[ T \right]-e_{1v}\,r_{2v} 
407\,\overline{\overline {\delta _{k+1/2} \left[ T \right]}} ^{\,j+1/2,k}} 
408\right)} \right]
410& +\delta 
411_k \left[ {A_w^{lT} \left(
412-e_{2w}\,r_{1w} \,\overline{\overline {\delta _{i+1/2} \left[ T \right]}} ^{\,i,k+1/2}
413\right.} \right.
415& \qquad \qquad \quad 
416-e_{1w}\,r_{2w} \,\overline{\overline {\delta _{j+1/2} \left[ T \right]}} ^{\,j,k+1/2}
418& \left. {\left. { 
419 \quad \quad \quad \left. {{\kern 
4201pt}+\frac{e_{1w}\,e_{2w} }{e_{3w} }\,\left( {r_{1w} ^2+r_{2w} ^2} 
421\right)\,\delta _{k+1/2} \left[ T \right]} \right)} \right]\;\;\;} \right\} 
422 \end{split}
423 \end{equation}
424where $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and the surface along which the diffusive operator acts ($i.e.$ horizontal or iso-neutral surfaces).
425 It is thus used when, in addition to \np{ln\_traldf\_lap}=T, we have \np{ln\_traldf\_iso}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=T. The way these slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using the mask technique (see \S\ref{LBC_coast}).
427The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical derivatives. For numerical stability, the vertical second derivative must be solved using the same implicit time scheme as those used in the vertical physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term is not computed in \mdl{traldf} module, but in \mdl{trazdf} module where, if iso-neutral mixing is used, the vertical mixing coefficient is simply increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$.
429This formulation conserves the tracer but does not ensure the decrease of the tracer variance. Nevertheless the treatment performed on the slopes (see \S\ref{LDF}) allows to run safely without any additional background horizontal diffusion \citep{Guily2001}. An alternate scheme \citep{Griffies1998} which preserves both tracer and its variance is currently been tested in \NEMO.
431Note that in partial step $z$-coordinate (\np{ln\_zps}=T), the horizontal derivatives in \eqref{Eq_tra_ldf_iso} at the bottom level require a specific treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
433% -------------------------------------------------------------------------------------------------------------
434%        Iso-level bilaplacian operator
435% -------------------------------------------------------------------------------------------------------------
436\subsection{Iso-level bilaplacian operator (\mdl{traldf\_bilap}, \np{ln\_traldf\_bilap})}
439The lateral fourth order operator formulation on tracers is obtained by applying (\ref{Eq_tra_ldf_lap}) twice. It requires an additional assumption on boundary conditions: first and third derivative terms normal to the coast are set to zero.
441It is used when, in addition to \np{ln\_traldf\_bilap}=T, we have \np{ln\_traldf\_level}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=F. In both cases, it can contributes to diapycnal mixing even if it should be less than in the laplacian case. It is therefore not recommended to use it.
443\textit{Notes:} In the code, the bilaplacian routine does not call twice the laplacian
444routine but is rather a specific routine. This is due to the fact that we
445introduce the eddy diffusivity coefficient, A, in the operator as: $\nabla 
446\cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$ and instead of
447$-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$ where
448$a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations
449ensure the total variance decrease, but the former requires a larger number
450of code-lines. It will be corrected in a forthcoming release.
452% -------------------------------------------------------------------------------------------------------------
453%        Rotated bilaplacian operator
454% -------------------------------------------------------------------------------------------------------------
455\subsection{Rotated bilaplacian operator (\mdl{traldf\_bilapg}, \np{ln\_traldf\_bilap})}
458The lateral fourth order operator formulation on tracers is obtained by applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption on boundary conditions: first and third derivative terms normal to the coast, the bottom and the surface are set to zero.
460 It is used when, in addition to \np{ln\_traldf\_bilap}=T, we have \np{ln\_traldf\_iso}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=T. Nevertheless, this rotated bilaplacian operator has never been seriously tested. No warranties that it is neither free of bugs or correctly formulated.
461Moreover, the stability range of such an operator will be probably quite narrow, requiring a significantly smaller time-step than the one used on unrotated operator.
463% ================================================================
464% Tracer Vertical Diffusion
465% ================================================================
466\section{Tracer Vertical Diffusion (\mdl{trazdf})}
472The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, based on a laplacian operator. The vertical diffusive operator given by (\ref{Eq_PE_zdf}) takes the following semi-discrete space form:
473\begin{equation} \label{Eq_tra_zdf}
475D^{vT}_T &= \frac{1}{e_{3T}} \; \delta_k \left[
476\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T]   \right]
478D^{vS}_T &= \frac{1}{e_{3T}} \; \delta_k \left[
479\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S]   \right]
483where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity
484coefficients on Temperature and Salinity, respectively. Generally, $A_w^{vT}=A_w^{vS}$ ecept when double diffusion mixing is parameterised (\key{zdfddm} defined). The way these coefficients can be evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when iso-neutral mixing is used, the both mixing coefficient are increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}.
486At the surface and bottom boundaries, the turbulent fluxes of momentum, heat and salt must be specified. At the surface they are prescribed from the surface forcing (see \S\ref{TRA_sbc}), while at the bottom they are set to zero for heat and salt, unless a geothermal flux forcing is prescribed as a bottom boundary condition (\S\ref{TRA_bbc}).
488The large eddy coefficient found in the mixed layer together with high
489vertical resolution implies a too restrictive constraint on the time step in
490explicit time stepping case (\np{ln\_zdfexp}=True). Therefore, the default implicit time stepping is generally preferred for the vertical diffusion as it overcomes the stability
491constraint. A forward time differencing scheme (\np{ln\_zdfexp}=T) using a time splitting technique (\np{n\_zdfexp} $> 1$) is provided as an alternative. Namelist variables
492\np{ln\_zdfexp} and \np{n\_zdfexp} apply to both tracers and dynamics.
494% ================================================================
495% External Forcing
496% ================================================================
497\section{External Forcing}
500% -------------------------------------------------------------------------------------------------------------
501%        surface boundary condition
502% -------------------------------------------------------------------------------------------------------------
503\subsection{surface boundary condition (\mdl{trasbc})}
506The surface boundary condition for tracers is implemented in a separate
507module (\mdl{trasbc}) instead of entering as a boundary condition on the vertical
508diffusion operator (as in the case of momentum). This has been found to
509enhance readability of the code. The two formulations are completely
510equivalent; the forcing terms in trasbc are the surface fluxes divided by
511the thickness of the top model layer. Following \citet{Roullet2000} the
512forcing on an ocean tracer, $c$, can be split into two parts: $F_{ext}^C $, the
513flux of tracer crossing the sea surface and not linked with the water
514exchange d at the surface with the atmosphere, and $F_{wf}^C$ the forcing
515on the concentration associated with the water flux. The latter forcing has
516also two components: a direct effect of change in concentration associated
517with the tracer carried by the water flux, and an indirect concentration/dilution effect :
520 F^C &= F_{ext} + F_{wf}^d                                          +F_{wf}^i    \\
522        &= F_{ext} - \left( c_E \, E - c_p \,P - c_R \,R \right) +c\left( E-P-R \right)
526Two cases must be distinguished, the nonlinear free surface case (\key{vvl} defined) and the linear free surface case. The first case is simpler, because the indirect concentration/dilution effect is naturally taken into
527account by letting the vertical scale factors vary in time. The salinity of
528water exchanged at the surface is assumed to be zero, so there is no salt
529flux at the free surface, excepted in the presence of sea ice. The heat flux
530at the free surface is the sum of $F_{ext}$, the direct heating/cooling
531(by the total non-penetrative heat flux) and $F_{wf}^e$ the heat carried by
532the water exchanged through the surface (evaporation, precipitation,
533runoff). The temperature of precipitations is not well known. In the model
534we assume that this water has the same temperature as the sea surface
535temperature, The resulting forcing terms for temperature T and salinity S
537\begin{equation} \label{Eq_tra_forcing}
539 &F^T =\frac{Q_{ns} }{\rho _o \;C_p \,e_{3T} }-\frac{\text{EMP}\;\left. T
540\right|_{k=1} }{e_{3T} }  & \\ 
542& F^S =\frac{\text{EMP}_S\;\left. S \right|_{k=1} }{e_{3T} }   &
543 \end{aligned}
546where EMP is the freshwater budget (evaporation minus precipitation minus river runoff) which forces the ocean volume, $Q_{ns}$ is the non-penetrative part of the net surface heat flux (difference between the total surface heat flux and the fraction of the short wave flux that penetrates in the water column), the product EMP$_S\;.\left. S \right|_{k=1}$ is  the ice-ocean salt flux, and $\left. S\right|_{k=1}$ is the sea surface salinity (\textit{SSS}). The total salt content is conserved in this formulation (excepted for the effect of the Asselin filter).
547%AMT note: the ice-ocean flux had been forgotten in the first release of the key_vvl option,
548%AMT  has this been corrected in the code?
550In the second case (linear free surface), the vertical scale factors are fixed in time so that the concentration/dilution effect must be added in trasbc. Because of the hypothesis made for the temperature of precipitation and runoffs, for temperature $F_{wf}^e +F_{wf}^i =0$. The resulting forcing term for temperature is:
552\begin{equation} \label{Eq_tra_forcing_q}
553F^T=\frac{Q_{ns} }{\rho _o \;C_p \,e_{3T} }
556The salinity forcing is still given by \eqref{Eq_tra_forcing} but the definition of EMP$_S$ is
557different: it is the total surface freshwater budget (evaporation minus
558precipitation minus river runoff plus the rate of change of the sea ice
559thickness). The total salt content is not exactly conserved (\citet{Roullet2000}, see also
562In the case of the rigid lid approximation, the surface salinity forcing $F^s$ 
563is also expressed by \eqref{Eq_tra_forcing} but now the global integral of the product
564EMP*S is not compensated by the advection of fluid through the top level: in
565the rigid lid case (contrary to the linear free surface), because \textit{w(k=1) = 0}. As a
566result, even if the budget of \textit{EMP} is zero in average over the whole ocean
567domain, the associated salt flux is not, as sea-surface salinity and \textit{EMP} are
568intrinsically correlated (high \textit{SSS} are found where evaporation is strong while
569low \textit{SSS} is usually associated with high precipitation or river runoff input).
571The $Q_{ns} $ and \textit{EMP} fields are defined and updated in \mdl{sbcmod} module (see
574% -------------------------------------------------------------------------------------------------------------
575%        Solar Radiation Penetration
576% -------------------------------------------------------------------------------------------------------------
577\subsection{Solar Radiation Penetration (\mdl{traqsr})}
583When the penetrative solar radiation option is used (\np{ln\_flxqsr}=T, the solar radiation penetrates the top few meters of the ocean, otherwise all the heat flux is absorbed in the first ocean level (\np{ln\_flxqsr}=F). A term is thus added to the time evolution equation of temperature \eqref{Eq_PE_tra_T} while the surface boundary condition is modified to take into account only the non-penetrative part of the surface heat flux:
584\begin{equation} \label{Eq_PE_qsr}
586\frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}   \\
587Q_{ns} &= Q_\text{Total} - Q_{sr}
591where $I$ is the downward irradiance. The additional term in \eqref{Eq_PE_qsr} is discretized as follows:
592\begin{equation} \label{Eq_tra_qsr}
593\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]
596A formulation including extinction coefficients is assumed for the downward irradiance $I$ 
598\begin{equation} \label{Eq_traqsr_iradiance}
599I(z) = Q_{sr} \left[Re^{-z / \xi_1} + \left( 1-R\right) e^{-z / \xi_2} \right]
601where $Q_{sr}$ is the penetrative part of the surface heat flux,
602$\xi_1$ and $\xi_2$ are two extinction length scales and $R$ determines the relative
603contribution of the two terms. The default values used correspond to a Type
604I water in Jerlov's [1968] classification: $\xi_1 = 0.35m$, $\xi_2 = 0.23m$ and
605$R = 0.58$ ((corresponding to \np{xsi1}, \np{xsi2} and \np{rabs} namelist parameters, respectively). $I$ is masked (no flux through the ocean bottom), so all the solar radiation that reaches the last ocean level is absorbed in that level. The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation is added to the temperature trend and the surface heat flux modified in routine \mdl{traqsr}. Note that in $z$-coordinates, the depth of $T-$levels depends on the single variable $k$. A one dimensional array of the coefficients $gdsr(k) = Re^{-z_w (k)/\xi_1} + (1-R)e^{-z_w (k)/\xi_2}$ can then be computed once and saved in central memory. Moreover \textit{nksr}, the level at which $gdrs$ becomes negligible (less than the computer precision) is computed once and the trend associated with the penetration of the solar radiation is only added until that level. At last, note that when the ocean is shallow (< 200~m), the part of the solar radiation can reach the ocean floor. In this case, we have chosen that all the radiation is absorbed at the last ocean level ($i.e.$ $I_w$ is masked).
607When coupling with a biology model (PISCES or LOBSTER), it is possible to calculate the light attenuation using information from the biology model. At the time of this writing, reading the light attenuation from a file is not implemented yet in the reference version.
609\colorbox{yellow}{case 4 bands and bio-coupling to add !!!}
611% -------------------------------------------------------------------------------------------------------------
612%        Bottom Boundary Condition
613% -------------------------------------------------------------------------------------------------------------
614\subsection{Bottom Boundary Condition (\mdl{trabbc} + \key{bbc})}
620\begin{figure}[!t] \label{Fig_geothermal}  \begin{center}
622\caption{Geothermal Heat flux (in $mW.m^{-2}$) as inferred from the age of the sea floor and the formulae of \citet{Stein1992}.}
623\end{center}   \end{figure}
626Usually it is considered that there is no exchange of heat nor salt through
627the ocean bottom, i.e. a no flux boundary condition is applied on active
628tracers at the bottom. This is the default option in NEMO, and it is
629implemented using the masking technique. Nevertheless, there exists a
630non-zero heat flux across the seafloor that is associated with the solid
631earth cooling. This flux is weak compared with surface fluxes --- a mean
632global value of $\sim0.1\;W/m^2$ \citep{Stein1992} --- but it is
633systematically positive and it acts only on the densest water masses. Taking
634this flux into account in a global ocean model increases by a few Sverdrups
635the deepest overturning cell (i.e. the one associated with the Antarctic
636Bottom Water).
638The presence or not of a geothermal heating is controlled by the namelist
639parameter  \np{ngeo\_flux}. Set to 1, a constant geothermal heatingis
640introducted which value is given by the \np{ngeo\_flux\_const}, also a
641namelist parameter. Set to 2, a spatially varying geothermal heat flux is
642introducted which is provided in the geothermal\ NetCDF
643file (Fig.\ref{Fig_geothermal}).
645% ================================================================
646% Bottom Boundary Layer
647% ================================================================
648\section{Bottom Boundary Layer (\mdl{trabbl} + \key{bbl\_diff} or \key{bbl\_adv})}
654In z-coordinate configuration, the bottom topography is represented as a
655series of discrete steps. This is not adequate to represent gravity driven
656downslope flows. Such flows arise downstream of sills such as the Strait of
657Gibraltar, Bab El Mandeb, or Denmark Strait, where dense water formed in
658marginal seas flows into a basin filled with less dense water. The amount of
659entrainment that occurs in those gravity plumes is critical to determine the
660density and volume flux of the densest waters of the ocean, such as the
661Antarctic Bottom water, or the North Atlantic Deep Water. $z$-coordinate
662models tend to overestimate the entrainment because the gravity flow is
663mixed down vertically by convection as it goes ``downstairs'' following the
664step topography, sometimes over a thickness much larger than the thickness
665of the observed gravity plume. A similar problem occurs in $s$-coordinate when
666the thickness of the bottom level varies in large proportions downstream of
667a sill \citep{Willebrand2001}, and the thickness of the plume is not
670The idea of the bottom boundary layer parameterization first introduced by
671\citet{BeckDos1998} is to allow a direct communication between
672two adjacent bottom cells at varying level, whenever the densest water is
673located above the less dense water. The communication can be by diffusive
674fluxes (diffusive BBL), advective fluxes (advective BBL) or both. Only
675tracers are modified, not the velocities. Implementing a BBL
676parameterization for momentum is a more complex problem because of the
677pressure gradient errors.
679% -------------------------------------------------------------------------------------------------------------
680%        Diffusive BBL
681% -------------------------------------------------------------------------------------------------------------
682\subsection{Diffusive Bottom Boundary layer (\mdl{trabbl})}
685The lateral diffusivity $A_l^\sigma $ in the BBL can be prescribed with a
686spatial dependence, e.g., in the conditional form
687\begin{equation} \label{Eq_tra_bbl}
688A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
689 \mbox{large}\quad if\;\nabla \rho \cdot \nabla H<0 \\ 
690 \\
691 0\quad \quad \;\,\mbox{otherwise} \\ 
692 \end{array}} \right.
695The large value of the coefficient when the diffusive BBL is active is given
696by the namelist parameter \np{atrbbl.}
698% -------------------------------------------------------------------------------------------------------------
699%        Advective BBL
700% -------------------------------------------------------------------------------------------------------------
701\subsection{Advective Bottom Boundary Layer (\mdl{trabb\_adv})}
704Implemented in NEMO v2.
706\colorbox{yellow} {Documentation to be added here }
708% ================================================================
709% Tracer damping
710% ================================================================
711\section{Tracer damping (\mdl{tradmp})}
717In some applications it can be useful to add a Newtonian damping term in the
718temperature and salinity equations:
719\begin{equation} \label{Eq_tra_dmp}
721 \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right\\
723 \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right)
724 \end{split}
725 \end{equation} 
726where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields (usually a climatology). The restoring term is added when \key{tradmp} is defined. It also requires that both \key{temdta} and \key{saldta} are defined ($i.e.$ that $T_o$ and $S_o$ are read). The restoring coefficient $S_o$ is a three-dimensional array initialized by the user in
727\rou{dtacof} routine also located in  \mdl{tradmp}.
729The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} the specification of
730the boundary conditions along artificial walls of a limited domain basin and
731\textit{(b)} the computation of the velocity field associated with a given $T$-$S$ field
732(for example to build the initial state of a prognostic simulation, or to
733use the resulting velocity field for a passive tracer study). The first case
734applies to regional models that have artificial walls instead of open
735boundaries. In the vicinity of these walls, $S_o$ takes large values
736(equivalent to a few day time scale) whereas it is zero in the interior of
737the model domain. The second case corresponds to the use of the robust
738diagnostic method \citep{Sarmiento1982}. It allows to find the velocity
739field consistent with the model dynamics while having a $T$-$S$ field close to a
740given climatology field ($T_o -S_o$). The time scale associated with
741$S_o$ is generally not a constant but spatially varying in order to respect
742some considerations. For example, it is usually set to zero in the mixed
743layer (defined either on a density or $S_o$ criterion) \citep{Madec1996} and in the equatorial region \citep{Reverdin1991, Fujio1991, MartiTh1992} as those two regions have a small time scale of adjustment,
744while smaller $S_o$ are used in the deep ocean where the typical time scale
745is long \citep{Sarmiento1982}. In addition it is reduced (and even
746zero) along the western boundary to allow the model to reconstruct its own
747western boundary structure in equilibrium with its physics. The choice of a
748Newtonian damping acting in the mixed layer or not is controlled by \np{nmldmp} 
749(\textbf{namelist} \np{nmldmp}parameter).
751The robust diagnostic method is very efficient to prevent the temperature drift in intermediate waters but it produces artificial sources of heat and salt within the ocean. It has also undesirable effects on the ocean convection. It tends to prevent deep convection and subsequent deep-water formation by stabilising too much the water columns.
753An example of computation of $S_o$ for robust diagnostic experiments with the ORCA2 model is provided in the \mdl{tradmp} module (subroutines \rou{dtacof} and \rou{cofdis} which compute coefficient and the distance to the bathymetry, respectively). Those routines are provided as examples and can be customised by the user.
755% ================================================================
756% Tracer time evolution
757% ================================================================
758\section{Tracer time evolution (\mdl{tranxt})}
764The general framework of dynamics time stepping is a leap-frog scheme, $i.e.$ a three level centred time scheme associated with a Asselin time filter (cf. \S\ref{DOM_nxt}):
765\begin{equation} \label{Eq_tra_nxt}
767T^{t+\Delta t} &= T^{t-\Delta t} + 2 \, \Delta t  \ \text{RHS}_T^t   \\
769T_f^\;\ \quad &= T^t \;\quad +\gamma \,\left[ {T_f^{t-\Delta t} -2T^t+T^{t+\Delta t}} \right]
773where $\text{RHS}_T$ is the right hand side of the temperature equation, the subscript $f$ denotes
774filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is initialized as \np{atfp} (\textbf{namelist} parameter). Its default value is \np{atfp=0.1}.
776When the vertical mixing is solved implicitly, the update of the next tracer
777fields is done in module \mdl{trazdf}. In that case only the swap of arrays
778and the Asselin filtering is done in \mdl{tranxt} module.
780In order to prepare the computation of the next time step, a swap of tracer arrays is performed: $T^{t-\Delta t} = T^t$ and $T^t = T_f$.
782% ================================================================
783% Equation of State (eosbn2)
784% ================================================================
785\section{Equation of State (\mdl{eosbn2}) }
791% -------------------------------------------------------------------------------------------------------------
792%        Equation of State
793% -------------------------------------------------------------------------------------------------------------
794\subsection{Equation of State (\np{neos} = 0, 1 or 2)}
797It is necessary to know the equation of state for the ocean very accurately to determine stability properties (especially the Brunt-Vais\"{a}l\"{a} frequency), particularly in the deep ocean. The ocean density is a non linear empirical function of \textit{in situ }temperature, salinity and pressure. The reference is the equation of state defined by the Joint Panel on Oceanographic Tables
798and Standards \citep{UNESCO1983}. It was the standard equation of state used in early releases of OPA. Even though this computation is fully vectorised, it is quite time consuming ($15$ to $20${\%} of the total CPU time) as it requires the prior computation of the \textit{in situ} temperature from the model \textit{potential} temperature using the \citep{Bryden1973} polynomial for adiabatic lapse rate and a $4^th$ order Runge-Kutta integration scheme. Since OPA6, we have chosen the \citet{JackMcD1995} equation of state for seawater. It allows the computation of the \textit{in situ} ocean density directly as a function of \textit{potential} temperature relative to the sea surface (an OPA variable), the practical salinity (another OPA variable) and the pressure (assuming no pressure variation along geopotential surfaces, i.e. the pressure in decibars is
799approximated by the depth in meters). Both the \citet{UNESCO1983} and \citet{JackMcD1995} equations of state have the same expression except that the values of the various coefficients have been adjusted by \citet{JackMcD1995} in order to use directly the \textit{potential} temperature instead of the \textit{in situ} one. This reduces the CPU time of the in situ density computation to about $3${\%} of the total CPU time, while maintaining a quite accurate equation of state.
801In the computer code, a \textit{true} density, $d$, is computed, i.e. the ratio of seawater volumic mass over $\rho_o$, a reference volumic mass (\textit{rau0} defined in \mdl{phycst}, usually $rau0= 1,020~Kg/m^3$). The default option (\np{neos}=0) is the \citet{JackMcD1995} equation of state. It is highly recommended to use it. Nevertheless, for process studies, it is often convenient to use a linear approximation of the density$^{\ast}$\footnote{$^{\ast }$ With the linear equation of state there is no longer
802a distinction between \textit{in situ} and \textit{potential} density. Cabling and thermobaric effects are also removed.}. Two linear formulations are available: a function of $T$ only (\np{neos}=1) and a function of both $T$ and $S$ (\np{neos}=2):
803\begin{equation} \label{Eq_tra_eos_linear}
805 d(T)    &= {\rho (T)} / {\rho _0 } &&= 1.028 - \alpha \;T     \\ 
806 d(T,S) &= {\rho (T,S)}                &&= \ \ \ \beta \;S - \alpha \;T
809where $\alpha $ and $\beta $ are the thermal and haline expansion coefficients, and $\rho_o$, the reference volumic mass, $rau0$. $\alpha $ and $\beta$ can be modified through \np{ralpha} and \np{rbeta} namelist parameters). Note that when $d$ is a function of $T$ only (\np{neos}=1), the salinity is a passive tracer and can be used as such.
812% -------------------------------------------------------------------------------------------------------------
813%        Brunt-Vais\"{a}l\"{a} Frequency
814% -------------------------------------------------------------------------------------------------------------
815\subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{neos} = 0, 1 or 2)}
818An accurate computation of the ocean stability (i.e. of $N$, the brunt-Vais\"{a}l\"{a} frequency) is of paramount importance as it is used in several ocean parameterisations (namely TKE, KPP, Richardson number dependent vertical diffusion, enhanced vertical diffusion, non-penetrative convection, iso-neutral diffusion). In particular, one must be aware that $N^2$ has to be computed with an \textit{in situ} reference. The expression of $N^2$ depends on the type of equation of state used (\np{neos} namelist parameter).
820For \np{neos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987} 
821polynomial expression is used with the pressure in decibar approximated by
822the depth in meters:
823\begin{multline} \label{Eq_tra_bn2}
824N^2 = \frac{g}{e_{3w}} \; \beta ( \overline{T}^{\,k+1/2},\widetilde{S},z_w )   \\
825      \left\{  \alpha / \beta ( \overline{T}^{\,k+1/2},\widetilde{S},z_w )
826         \ \delta_{k+1/2}[T]     - \delta_{k+1/2}[S]
827       \right\}
829where $T$ is the \textit{potential} temperature, $\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$ a salinity anomaly, and $\alpha$ ($\beta\,$) the thermal (haline) expansion coefficient. Both $\alpha$ and $\beta$ depend on \textit{potential} temperature, salinity which are averaged at $w$-points prior to the computation.
831When a linear equation of state is used (\np{neos}=1 or 2, \eqref{Eq_tra_bn2} reduces to:
832\begin{equation} \label{Eq_tra_bn2_linear}
833N^2 = \frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
835where $\alpha$ and $\beta $ are the constant coefficients used to defined the linear
836equation of state \eqref{Eq_tra_eos_linear}.
838% -------------------------------------------------------------------------------------------------------------
839%        Specific Heat
840% -------------------------------------------------------------------------------------------------------------
841\subsection{Specific Heat (\textit{rcp}, \mdl{phycst})}
844The specific heat of sea water, $C_p$, is a function of temperature, salinity and pressure \citep{UNESCO1983}. It is only used in the model to convert surface heat fluxes into surface temperature increase, thus the pressure dependence is neglected. The dependence on $T$ and $S$ is weak. For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$ when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$. Its computer name is \textit{rcp} and its value is set in \mdl{phycst} module.
846% -------------------------------------------------------------------------------------------------------------
847%        Freezing Point of Seawater
848% -------------------------------------------------------------------------------------------------------------
849\subsection{Freezing Point of Seawater (\mdl{ocfzpt})}
852The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}:
853\begin{equation} \label{Eq_tra_eos_fzp}
854   \begin{split}
855T_f (S,p) &= \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} 
856                       -  2.154996 \;10^{-4} \,\right) \ S    \\
857               & - 7.53\,10^{-3}\,p
858   \end{split}
861\eqref{Eq_tra_eos_fzp} is only used to compute the potential freezing point of sea water
862($i.e.$ referenced to the surface $p=0$), thus the pressure dependent terms in \eqref{Eq_tra_eos_fzp} (last term) has been dropped. The \textit{before} and \textit{now} surface freezing point is introduced in the code as $fzptb$ and $fzptn$ 2D arrays together with a  \textit{now} mask (\textit{freezn}) which takes 0 or 1 whether the ocean temperature is above or at the freezing point. Caution: do not confuse \textit{freezn} with the fraction of lead (\textit{frld}) defined in LIM.
864% ================================================================
865% Horizontal Derivative in zps-coordinate
866% ================================================================
867\section{Horizontal Derivative in \textit{zps}-coordinate (\mdl{zpshde})}
870With partial bottom cells (\np{ln\_zps}=T), tracers in horizontally adjacent cells generally live at different depths. Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure gradient (\mdl{dynhpg} module). Before taking horizontal gradients between the tracers next to the bottom, a linear interpolation is used to approximate the deeper tracer as if it actually lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}). For example on temperature in the i-direction, the needed interpolated temperature, $\widetilde{T}$, is:
872\begin{figure}[!p] \label{Fig_Partial_step_scheme}  \begin{center}
874\caption{ Discretisation of horizontal derivative and mean of tracers in z-partial step coordinate (\np{ln\_zps}=T) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. The horizontal derivative is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and the mean by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$}
875\end{center}   \end{figure}
878\widetilde{T}= \left\{  \begin{aligned} 
879&T^{\,i+1}      -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta _k T^{i+1} 
880                        && \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$   }  \\
881                              \\
882&T^{\,i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i       }\;\delta _k T^{i+1}
883                        && \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
884            \end{aligned}   \right.
886and the resulting formulation of horizontal derivative and horizontal mean value of $T$ at $U$-point are:
887\begin{equation} \label{Eq_zps_hde}
889 \delta _{i+1/2} T=  \begin{cases}
890\ \ \ \widetilde {T}\quad\ -T^i     & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\
891                              \\
892\ \ \ T^{\,i+1}-\widetilde{T}    & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
893                  \end{cases}     \\
895\overline {T}^{\,i+1/2} \ =   \begin{cases}
896( \widetilde {T}\ \ \;\,-T^{\,i})    / 2  & \;\ \ \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\
897                              \\
898( T^{\,i+1}-\widetilde{T} ) / 2     & \;\ \ \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
899            \end{cases}
903The computation of horizontal derivative of tracers as well as of density is performed once for all at each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed. It has to be emphasized that the computation of the interpolated density, $\widetilde{\rho}$, is not identical to the one of $T$ and $S$. Instead of forming a linear approximation of density, we compute $\widetilde{\rho }$ from the interpolated value of $T$ and $S$, and the pressure of at $u$-point (in the equation of state pressure is approximated by depth, see \S\ref{TRA_eos} ) :
904\begin{equation} \label{Eq_zps_hde_rho}
905\widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u })
906\quad \text{where }\  z_u = \min \left( {z_T^{i+1} ,z_T^i } \right)
909This is a much better approximation as the variation of $\rho$ with depth (and thus pressure) is highly non-linear with a true equation of state and thus is badly approximated with a linear interpolation. This approximation is used to compute both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral surfaces (\S\ref{LDF_slp})
911\textit{Notes}: in almost all the advection schemes presented in this Chapter, both mean and derivative operators appear. Yet, it has been chosen not to use \eqref{Eq_zps_hde} in those schemes.: contrary to diffusion and pressure gradient computation, no correction for partial steps is applied for advection.The main motivation was to preserve the domain averaged mean variance of the field advected when using $2^{nd}$ order centred scheme. Sensitivity of the advection schemes to the way horizontal means are performed in the vicinity of partial cells should be further investigated in a near future.
Note: See TracBrowser for help on using the repository browser.