# source:trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex@3294

Last change on this file since 3294 was 3294, checked in by rblod, 9 years ago

Merge of 3.4beta into the trunk

File size: 75.1 KB
Line
1% ================================================================
2% Chapter 1 Ñ Model Basics
3% ================================================================
4
5\chapter{Model basics}
6\label{PE}
7\minitoc
8
9\newpage
10$\$\newline    % force a new ligne
11
12% ================================================================
13% Primitive Equations
14% ================================================================
15\section{Primitive Equations}
16\label{PE_PE}
17
18% -------------------------------------------------------------------------------------------------------------
19%        Vector Invariant Formulation
20% -------------------------------------------------------------------------------------------------------------
21
22\subsection{Vector Invariant Formulation}
23\label{PE_Vector}
24
25
26The ocean is a fluid that can be described to a good approximation by the primitive
27equations, $i.e.$ the Navier-Stokes equations along with a nonlinear equation of
28state which couples the two active tracers (temperature and salinity) to the fluid
30
31\textit{(1) spherical earth approximation: }the geopotential surfaces are assumed to
32be spheres so that gravity (local vertical) is parallel to the earth's radius
33
34\textit{(2) thin-shell approximation: }the ocean depth is neglected compared to the earth's radius
35
36\textit{(3) turbulent closure hypothesis: }the turbulent fluxes (which represent the effect
37of small scale processes on the large-scale) are expressed in terms of large-scale features
38
39\textit{(4) Boussinesq hypothesis:} density variations are neglected except in their
40contribution to the buoyancy force
41
42\textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a
43balance between the vertical pressure gradient and the buoyancy force (this removes
44convective processes from the initial Navier-Stokes equations and so convective processes
46
47\textit{(6) Incompressibility hypothesis: }the three dimensional divergence of the velocity
48vector is assumed to be zero.
49
50Because the gravitational force is so dominant in the equations of large-scale motions,
51it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked
52to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two
53vectors orthogonal to \textbf{k}, $i.e.$ tangent to the geopotential surfaces. Let us define
54the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$
55(the subscript $h$ denotes the local horizontal vector, $i.e.$ over the (\textbf{i},\textbf{j}) plane),
56$T$ the potential temperature, $S$ the salinity, \textit{$\rho$} the \textit{in situ} density.
57The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k})
58vector system provides the following six equations (namely the momentum balance, the
59hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation
60equations and an equation of state):
61\begin{subequations} \label{Eq_PE}
62  \begin{equation}     \label{Eq_PE_dyn}
63\frac{\partial {\rm {\bf U}}_h }{\partial t}=
64-\left[    {\left( {\nabla \times {\rm {\bf U}}} \right)\times {\rm {\bf U}}
65            +\frac{1}{2}\nabla \left( {{\rm {\bf U}}^2} \right)}    \right]_h
66 -f\;{\rm {\bf k}}\times {\rm {\bf U}}_h
67-\frac{1}{\rho _o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} + {\rm {\bf F}}^{\rm {\bf U}}
68  \end{equation}
69  \begin{equation}     \label{Eq_PE_hydrostatic}
70\frac{\partial p }{\partial z} = - \rho \ g
71  \end{equation}
72  \begin{equation}     \label{Eq_PE_continuity}
73\nabla \cdot {\bf U}=  0
74  \end{equation}
75\begin{equation} \label{Eq_PE_tra_T}
76\frac{\partial T}{\partial t} = - \nabla \cdot  \left( T \ \rm{\bf U} \right) + D^T + F^T
77  \end{equation}
78  \begin{equation}     \label{Eq_PE_tra_S}
79\frac{\partial S}{\partial t} = - \nabla \cdot  \left( S \ \rm{\bf U} \right) + D^S + F^S
80  \end{equation}
81  \begin{equation}     \label{Eq_PE_eos}
82\rho = \rho \left( T,S,p \right)
83  \end{equation}
84\end{subequations}
85where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions,
86$t$ is the time, $z$ is the vertical coordinate, $\rho$ is the \textit{in situ} density given by
87the equation of state (\ref{Eq_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure,
88$f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration (where $\bf \Omega$ is the Earth's
89angular velocity vector), and $g$ is the gravitational acceleration.
90${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale
91physics for momentum, temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$
92and $F^S$ surface forcing terms. Their nature and formulation are discussed in
93\S\ref{PE_zdf_ldf} and page \S\ref{PE_boundary_condition}.
94
95.
96
97% -------------------------------------------------------------------------------------------------------------
98% Boundary condition
99% -------------------------------------------------------------------------------------------------------------
100\subsection{Boundary Conditions}
101\label{PE_boundary_condition}
102
103An ocean is bounded by complex coastlines, bottom topography at its base and an air-sea
104or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$
105and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ is the height
106of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$,
107chosen as a mean sea surface (Fig.~\ref{Fig_ocean_bc}). Through these two boundaries,
108the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth,
109the continental margins, the sea ice and the atmosphere. However, some of these fluxes are
110so weak that even on climatic time scales of thousands of years they can be neglected.
111In the following, we briefly review the fluxes exchanged at the interfaces between the ocean
112and the other components of the earth system.
113
114%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
115\begin{figure}[!ht]   \begin{center}
116\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_I_ocean_bc.pdf}
117\caption{    \label{Fig_ocean_bc}
118The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$, where $H$
119is the depth of the sea floor and $\eta$ the height of the sea surface.
120Both $H$ and $\eta$ are referenced to $z=0$.}
121\end{center}   \end{figure}
122%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
123
124
125\begin{description}
126\item[Land - ocean interface:] the major flux between continental margins and the ocean is
127a mass exchange of fresh water through river runoff. Such an exchange modifies the sea
128surface salinity especially in the vicinity of major river mouths. It can be neglected for short
129range integrations but has to be taken into account for long term integrations as it influences
130the characteristics of water masses formed (especially at high latitudes). It is required in order
131to close the water cycle of the climate system. It is usually specified as a fresh water flux at
132the air-sea interface in the vicinity of river mouths.
133\item[Solid earth - ocean interface:] heat and salt fluxes through the sea floor are small,
134except in special areas of little extent. They are usually neglected in the model
135\footnote{In fact, it has been shown that the heat flux associated with the solid Earth cooling
136($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world
137ocean (see \ref{TRA_bbc}).}.
138The boundary condition is thus set to no flux of heat and salt across solid boundaries.
139For momentum, the situation is different. There is no flow across solid boundaries,
140$i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words,
141the bottom velocity is parallel to solid boundaries). This kinematic boundary condition
142can be expressed as:
143\begin{equation} \label{Eq_PE_w_bbc}
144w = -{\rm {\bf U}}_h \cdot  \nabla _h \left( H \right)
145\end{equation}
146In addition, the ocean exchanges momentum with the earth through frictional processes.
147Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized
148in terms of turbulent fluxes using bottom and/or lateral boundary conditions. Its specification
149depends on the nature of the physical parameterisation used for ${\rm {\bf D}}^{\rm {\bf U}}$
150in \eqref{Eq_PE_dyn}. It is discussed in \S\ref{PE_zdf}, page~\pageref{PE_zdf}.% and Chap. III.6 to 9.
151\item[Atmosphere - ocean interface:] the kinematic surface condition plus the mass flux
152of fresh water PE  (the precipitation minus evaporation budget) leads to:
153\begin{equation} \label{Eq_PE_w_sbc}
154w = \frac{\partial \eta }{\partial t}
155    + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot  \nabla _h \left( \eta \right)
156    + P-E
157\end{equation}
158The dynamic boundary condition, neglecting the surface tension (which removes capillary
159waves from the system) leads to the continuity of pressure across the interface $z=\eta$.
160The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat.
161\item[Sea ice - ocean interface:] the ocean and sea ice exchange heat, salt, fresh water
162and momentum. The sea surface temperature is constrained to be at the freezing point
163at the interface. Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the
164ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and
165salt fluxes that cannot be neglected.
166\end{description}
167
168
169%\newpage
170%$\$\newline    % force a new ligne
171
172% ================================================================
174% ================================================================
176\label{PE_hor_pg}
177
178% -------------------------------------------------------------------------------------------------------------
179% Pressure Formulation
180% -------------------------------------------------------------------------------------------------------------
181\subsection{Pressure Formulation}
182\label{PE_p_formulation}
183
184The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a
185reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that:
186$p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\ref{Eq_PE_hydrostatic}),
187assuming that pressure in decibars can be approximated by depth in meters in (\ref{Eq_PE_eos}).
188The hydrostatic pressure is then given by:
189\begin{equation} \label{Eq_PE_pressure}
190p_h \left( {i,j,z,t} \right)
191 = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma }
192\end{equation}
193 Two strategies can be considered for the surface pressure term: $(a)$ introduce of a
194 new variable $\eta$, the free-surface elevation, for which a prognostic equation can be
195 established and solved; $(b)$ assume that the ocean surface is a rigid lid, on which the
196 pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used,
197 one solution of the free-surface elevation consists of the excitation of external gravity waves.
198 The flow is barotropic and the surface moves up and down with gravity as the restoring force.
199 The phase speed of such waves is high (some hundreds of metres per second) so that
200 the time step would have to be very short if they were present in the model. The latter
201 strategy filters out these waves since the rigid lid approximation implies $\eta=0$, $i.e.$
202 the sea surface is the surface $z=0$. This well known approximation increases the surface
203 wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic
204 Rossby or planetary waves). The rigid-lid hypothesis is an obsolescent feature in modern
205 OGCMs. It has been available until the release 3.1 of  \NEMO, and it has been removed
206 in release 3.2 and followings. Only the free surface formulation is now described in the
207 this document (see the next sub-section).
208
209% -------------------------------------------------------------------------------------------------------------
210% Free Surface Formulation
211% -------------------------------------------------------------------------------------------------------------
212\subsection{Free Surface Formulation}
213\label{PE_free_surface}
214
215In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced
216which describes the shape of the air-sea interface. This variable is solution of a
217prognostic equation which is established by forming the vertical average of the kinematic
218surface condition (\ref{Eq_PE_w_bbc}):
219\begin{equation} \label{Eq_PE_ssh}
220\frac{\partial \eta }{\partial t}=-D+P-E
222D=\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]
223\end{equation}
224and using (\ref{Eq_PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$.
225
226Allowing the air-sea interface to move introduces the external gravity waves (EGWs)
227as a class of solution of the primitive equations. These waves are barotropic because
228of hydrostatic assumption, and their phase speed is quite high. Their time scale is
229short with respect to the other processes described by the primitive equations.
230
231Two choices can be made regarding the implementation of the free surface in the model,
232depending on the physical processes of interest.
233
234$\bullet$ If one is interested in EGWs, in particular the tides and their interaction
235with the baroclinic structure of the ocean (internal waves) possibly in shallow seas,
236then a non linear free surface is the most appropriate. This means that no
237approximation is made in (\ref{Eq_PE_ssh}) and that the variation of the ocean
238volume is fully taken into account. Note that in order to study the fast time scales
239associated with EGWs it is necessary to minimize time filtering effects (use an
240explicit time scheme with very small time step, or a split-explicit scheme with
241reasonably small time step, see \S\ref{DYN_spg_exp} or \S\ref{DYN_spg_ts}.
242
243$\bullet$ If one is not interested in EGW but rather sees them as high frequency
244noise, it is possible to apply an explicit filter to slow down the fastest waves while
245not altering the slow barotropic Rossby waves. If further, an approximative conservation
246of heat and salt contents is sufficient for the problem solved, then it is
247sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows
248to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}.
249
250The filtering of EGWs in models with a free surface is usually a matter of discretisation
251of the temporal derivatives, using the time splitting method \citep{Killworth_al_JPO91, Zhang_Endoh_JGR92}
252or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach
253developed by \citet{Roullet_Madec_JGR00}: the damping of EGWs is ensured by introducing an
254additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes:
255\begin{equation} \label{Eq_PE_flt}
256\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
257- g \nabla \left( \tilde{\rho} \ \eta \right)
258- g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)
259\end{equation}
260where $T_c$, is a parameter with dimensions of time which characterizes the force,
261$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$
262represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
263non-linear and viscous terms in \eqref{Eq_PE_dyn}.
264
265The new force can be interpreted as a diffusion of vertically integrated volume flux divergence.
266The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$
267and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime
268in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate,
269$i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than
270$T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs
271can be damped by choosing $T_c > \rdt$. \citet{Roullet_Madec_JGR00} demonstrate that
272(\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which
273has to be computed implicitly. This is not surprising since the use of a large time step has a
274necessarily numerical cost. Two gains arise in comparison with the previous formulations.
275Firstly, the damping of EGWs can be quantified through the magnitude of the additional term.
276Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as
277soon as $T_c > \rdt$.
278
279When the variations of free surface elevation are small compared to the thickness of the first
280model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized
281by \citet{Roullet_Madec_JGR00} the linearization of (\ref{Eq_PE_ssh}) has consequences on the
282conservation of salt in the model. With the nonlinear free surface equation, the time evolution
283of the total salt content is
284\begin{equation} \label{Eq_PE_salt_content}
285    \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv}
286                        =\int\limits_S {S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds}
287\end{equation}
288where $S$ is the salinity, and the total salt is integrated over the whole ocean volume
289$D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an
290integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh})
291is satisfied, so that the salt is perfectly conserved. When the free surface equation is
292linearized, \citet{Roullet_Madec_JGR00} show that the total salt content integrated in the fixed
293volume $D$ (bounded by the surface $z=0$) is no longer conserved:
294\begin{equation} \label{Eq_PE_salt_content_linear}
295         \frac{\partial }{\partial t}\int\limits_D {S\;dv}
296               = - \int\limits_S {S\;\frac{\partial \eta }{\partial t}ds}
297\end{equation}
298
299The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions
300\citep{Roullet_Madec_JGR00}. It can be significant when the freshwater forcing is not balanced and
301the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta$}
302results in a decrease of the salinity in the fixed volume $D$. Even in that case though,
303the total salt integrated in the variable volume $D_{\eta}$ varies much less, since
304(\ref{Eq_PE_salt_content_linear}) can be rewritten as
305\begin{equation} \label{Eq_PE_salt_content_corrected}
306\frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv}
307=\frac{\partial}{\partial t} \left[ \;{\int\limits_D {S\;dv} +\int\limits_S {S\eta \;ds} } \right]
308=\int\limits_S {\eta \;\frac{\partial S}{\partial t}ds}
309\end{equation}
310
311Although the total salt content is not exactly conserved with the linearized free surface,
312its variations are driven by correlations of the time variation of surface salinity with the
313sea surface height, which is a negligible term. This situation contrasts with the case of
314the rigid lid approximation in which case freshwater forcing is represented by a virtual
315salt flux, leading to a spurious source of salt at the ocean surface
317
318\newpage
319$\$\newline    % force a new ligne
320
321% ================================================================
322% Curvilinear z-coordinate System
323% ================================================================
324\section{Curvilinear \textit{z-}coordinate System}
325\label{PE_zco}
326
327
328% -------------------------------------------------------------------------------------------------------------
329% Tensorial Formalism
330% -------------------------------------------------------------------------------------------------------------
331\subsection{Tensorial Formalism}
332\label{PE_tensorial}
333
334In many ocean circulation problems, the flow field has regions of enhanced dynamics
335($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts).
336The representation of such dynamical processes can be improved by specifically increasing
337the model resolution in these regions. As well, it may be convenient to use a lateral
338boundary-following coordinate system to better represent coastal dynamics. Moreover,
339the common geographical coordinate system has a singular point at the North Pole that
340cannot be easily treated in a global model without filtering. A solution consists of introducing
341an appropriate coordinate transformation that shifts the singular point onto land
342\citep{Madec_Imbard_CD96, Murray_JCP96}. As a consequence, it is important to solve the primitive
343equations in various curvilinear coordinate systems. An efficient way of introducing an
344appropriate coordinate transform can be found when using a tensorial formalism.
345This formalism is suited to any multidimensional curvilinear coordinate system. Ocean
346modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth
347approximation), with preservation of the local vertical. Here we give the simplified equations
348for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey
349of the conservation laws of fluid dynamics.
350
351Let (\textit{i},\textit{j},\textit{k}) be a set of orthogonal curvilinear coordinates on the sphere
352associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k})
353linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are
354two vectors orthogonal to \textbf{k}, $i.e.$ along geopotential surfaces (Fig.\ref{Fig_referential}).
355Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined
356by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of
357the earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea
358level (Fig.\ref{Fig_referential}). The local deformation of the curvilinear coordinate system is
359given by $e_1$, $e_2$ and $e_3$, the three scale factors:
360\begin{equation} \label{Eq_scale_factors}
361\begin{aligned}
362 e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda
363}{\partial i}\cos \varphi } \right)^2+\left( {\frac{\partial \varphi
364}{\partial i}} \right)^2} \right]^{1/2} \\
365 e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda
366}{\partial j}\cos \varphi } \right)^2+\left( {\frac{\partial \varphi
367}{\partial j}} \right)^2} \right]^{1/2} \\
368 e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\
369 \end{aligned}
370 \end{equation}
371
372%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
373\begin{figure}[!tb]   \begin{center}
374\includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_I_earth_referential.pdf}
375\caption{   \label{Fig_referential}
376the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear
377coordinate system (\textbf{i},\textbf{j},\textbf{k}). }
378\end{center}   \end{figure}
379%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
380
381Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by
382$a$ in (\ref{Eq_scale_factors}) (thin-shell approximation). The resulting horizontal scale
383factors $e_1$, $e_2$  are independent of $k$ while the vertical scale factor is a single
384function of $k$ as \textbf{k} is parallel to \textbf{z}. The scalar and vector operators that
385appear in the primitive equations (Eqs. \eqref{Eq_PE_dyn} to \eqref{Eq_PE_eos}) can
386be written in the tensorial form, invariant in any orthogonal horizontal curvilinear coordinate
387system transformation:
388\begin{subequations} \label{Eq_PE_discrete_operators}
390\nabla q=\frac{1}{e_1 }\frac{\partial q}{\partial i}\;{\rm {\bf
391i}}+\frac{1}{e_2 }\frac{\partial q}{\partial j}\;{\rm {\bf j}}+\frac{1}{e_3
392}\frac{\partial q}{\partial k}\;{\rm {\bf k}}    \\
393\end{equation}
394\begin{equation} \label{Eq_PE_div}
395\nabla \cdot {\rm {\bf A}}
396= \frac{1}{e_1 \; e_2} \left[
397  \frac{\partial \left(e_2 \; a_1\right)}{\partial i }
398+\frac{\partial \left(e_1 \; a_2\right)}{\partial j }       \right]
399+ \frac{1}{e_3} \left[ \frac{\partial a_3}{\partial k }   \right]
400\end{equation}
401\begin{equation} \label{Eq_PE_curl}
402   \begin{split}
403\nabla \times \vect{A} =
404    \left[ {\frac{1}{e_2 }\frac{\partial a_3}{\partial j}
405            -\frac{1}{e_3 }\frac{\partial a_2 }{\partial k}} \right] \; \vect{i}
406&+\left[ {\frac{1}{e_3 }\frac{\partial a_1 }{\partial k}
407           -\frac{1}{e_1 }\frac{\partial a_3 }{\partial i}} \right] \; \vect{j}     \\
408&+\frac{1}{e_1 e_2 } \left[ {\frac{\partial \left( {e_2 a_2 } \right)}{\partial i}
409                                       -\frac{\partial \left( {e_1 a_1 } \right)}{\partial j}} \right] \; \vect{k}
410   \end{split}
411\end{equation}
412\begin{equation} \label{Eq_PE_lap}
413\Delta q = \nabla \cdot \left\nabla q \right)
414\end{equation}
415\begin{equation} \label{Eq_PE_lap_vector}
416\Delta {\rm {\bf A}} =
417  \nabla \left( \nabla \cdot {\rm {\bf A}} \right)
418- \nabla \times \left\nabla \times {\rm {\bf A}} \right)
419\end{equation}
420\end{subequations}
421where $q$ is a scalar quantity and ${\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinate system.
422
423% -------------------------------------------------------------------------------------------------------------
424% Continuous Model Equations
425% -------------------------------------------------------------------------------------------------------------
426\subsection{Continuous Model Equations}
427\label{PE_zco_Eq}
428
429In order to express the Primitive Equations in tensorial formalism, it is necessary to compute
430the horizontal component of the non-linear and viscous terms of the equation using
432Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate
433system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity
434field $\chi$, by:
435\begin{equation} \label{Eq_PE_curl_Uh}
436\zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v}
437\right)}{\partial i}-\frac{\partial \left( {e_1 \,u} \right)}{\partial j}}
438\right]
439\end{equation}
440\begin{equation} \label{Eq_PE_div_Uh}
441\chi =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,u}
442\right)}{\partial i}+\frac{\partial \left( {e_1 \,v} \right)}{\partial j}}
443\right]
444\end{equation}
445
446Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$
447and that $e_3$  is a function of the single variable $k$, the nonlinear term of
448\eqref{Eq_PE_dyn} can be transformed as follows:
449\begin{flalign*}
450&\left[ {\left( { \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}}
451+\frac{1}{2}   \nabla \left( {{\rm {\bf U}}^2} \right)}   \right]_h        &
452\end{flalign*}
453\begin{flalign*}
455 {\left[    {   \frac{1}{e_3} \frac{\partial u  }{\partial k}
456         -\frac{1}{e_1} \frac{\partial w  }{\partial i} } \right] w - \zeta \; v }     \\
457      {\zeta \; u - \left[ {   \frac{1}{e_2} \frac{\partial w}{\partial j}
458                     -\frac{1}{e_3} \frac{\partial v}{\partial k} } \right] \ w}  \\
459       \end{array} }} \right)
460+\frac{1}{2}   \left( {{\begin{array}{*{20}c}
461       { \frac{1}{e_1}  \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial i}}  \hfill    \\
462       { \frac{1}{e_2}  \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial j}}  \hfill    \\
463       \end{array} }} \right)       &
464\end{flalign*}
465\begin{flalign*}
467 {-\zeta \; v} \hfill \\
468 { \zeta \; u} \hfill \\
469         \end{array} }} \right)
470+\frac{1}{2}\left( {{   \begin{array}{*{20}c}
471 {\frac{1}{e_1 }\frac{\partial \left( {u^2+v^2} \right)}{\partial i}} \hfill  \\
472 {\frac{1}{e_2 }\frac{\partial \left( {u^2+v^2} \right)}{\partial j}} \hfill  \\
473                  \end{array} }} \right)
474+\frac{1}{e_3 }\left( {{      \begin{array}{*{20}c}
475 { w \; \frac{\partial u}{\partial k}}    \\
476 { w \; \frac{\partial v}{\partial k}}    \\
477                     \end{array} }} \right
478-\left( {{  \begin{array}{*{20}c}
479 {\frac{w}{e_1}\frac{\partial w}{\partial i}
480 -\frac{1}{2e_1}\frac{\partial w^2}{\partial i}} \hfill \\
481 {\frac{w}{e_2}\frac{\partial w}{\partial j}
482  -\frac{1}{2e_2}\frac{\partial w^2}{\partial j}} \hfill \\
483         \end{array} }} \right)        &
484\end{flalign*}
485
486The last term of the right hand side is obviously zero, and thus the nonlinear term of
487\eqref{Eq_PE_dyn} is written in the $(i,j,k)$ coordinate system:
488\begin{equation} \label{Eq_PE_vector_form}
489\left[ {\left( {  \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}}
490+\frac{1}{2}   \nabla \left( {{\rm {\bf U}}^2} \right)}   \right]_h
491=\zeta
492\;{\rm {\bf k}}\times {\rm {\bf U}}_h +\frac{1}{2}\nabla _h \left( {{\rm
493{\bf U}}_h^2 } \right)+\frac{1}{e_3 }w\frac{\partial {\rm {\bf U}}_h
494}{\partial k}
495\end{equation}
496
497This is the so-called \textit{vector invariant form} of the momentum advection term.
498For some purposes, it can be advantageous to write this term in the so-called flux form,
499$i.e.$ to write it as the divergence of fluxes. For example, the first component of
500\eqref{Eq_PE_vector_form} (the $i$-component) is transformed as follows:
501\begin{flalign*}
502&{ \begin{array}{*{20}l}
503\left[ {\left( {\nabla \times \vect{U}} \right)\times \vect{U}
504          +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i   % \\
505%\\
506     = - \zeta \;v
507     + \frac{1}{2\;e_1 } \frac{\partial \left( {u^2+v^2} \right)}{\partial i}
508     + \frac{1}{e_3}w \ \frac{\partial u}{\partial k}          \\
509\\
510\qquad =\frac{1}{e_1 \; e_2} \left(    -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i}
511                     +v\frac{\partial \left( {e_1 \,u} \right)}{\partial j}    \right)
512+\frac{1}{e_1 e_2 }\left(  +e_2 \; u\frac{\partial u}{\partial i}
513                     +e_2 \; v\frac{\partial v}{\partial i}              \right)
514+\frac{1}{e_3}       \left(   w\;\frac{\partial u}{\partial k}       \right)   \\
515\end{array} }        &
516\end{flalign*}
517\begin{flalign*}
518&{ \begin{array}{*{20}l}
520 -\left(        v^\frac{\partial e_2                                }{\partial i}
521      +e_2 \,v    \frac{\partial v                                   }{\partial i}     \right)
522+\left(           \frac{\partial \left( {e_1 \,u\,v}  \right)}{\partial j}
523      -e_1 \,u    \frac{\partial v                                   }{\partial j}  \right\right.
525+\left(           \frac{\partial \left( {e_2 u\,u}     \right)}{\partial i}
526      -u       \frac{\partial \left( {e_2 u}         \right)}{\partial i}  \right)
527+e_2 v            \frac{\partial v                                    }{\partial i}
528                  \right\}
529+\frac{1}{e_3} \left(
530               \frac{\partial \left( {w\,u} \right)         }{\partial k}
531       -u         \frac{\partial w                    }{\partial k}  \right) \\
532\end{array} }     &
533\end{flalign*}
534\begin{flalign*}
535&{ \begin{array}{*{20}l}
537               \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i}
538      +        \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j}  \right)
539+\frac{1}{e_3 }      \frac{\partial \left( {w\,u       } \right)}{\partial k}
541+\frac{1}{e_1 e_2 }     \left(
542      -u \left(   \frac{\partial \left( {e_1 v   } \right)}{\partial j}
543               -v\,\frac{\partial e_1 }{\partial j}             \right)
544      -u       \frac{\partial \left( {e_2 u   } \right)}{\partial i}
545                  \right)
546 -\frac{1}{e_3 }     \frac{\partial w}{\partial k} u
547 +\frac{1}{e_1 e_2 }\left(    -v^2\frac{\partial e_2   }{\partial i}     \right)
548\end{array} }     &
549\end{flalign*}
550\begin{flalign*}
551&{ \begin{array}{*{20}l}
552\qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right)
553-   \left( \nabla \cdot {\rm {\bf U}} \right) \ u
554+\frac{1}{e_1 e_2 }\left(
555      -v^2     \frac{\partial e_2 }{\partial i}
556      +uv   \,    \frac{\partial e_1 }{\partial j}    \right) \\
557\end{array} }     &
558\end{flalign*}
559as $\nabla \cdot {\rm {\bf U}}\;=0$ (incompressibility) it comes:
560\begin{flalign*}
561&{ \begin{array}{*{20}l}
562\qquad = \nabla \cdot \left{{\rm {\bf U}}\,u}      \right)
563\frac{1}{e_1 e_2 }   \left( v \; \frac{\partial e_2}{\partial i}
564                         -u \; \frac{\partial e_1}{\partial j}    \right\left( -v \right)
565\end{array} }     &
566\end{flalign*}
567
568The flux form of the momentum advection term is therefore given by:
569\begin{multline} \label{Eq_PE_flux_form}
570      \left[
571  \left(    {\nabla \times {\rm {\bf U}}}    \right) \times {\rm {\bf U}}
572+\frac{1}{2}   \nabla \left{{\rm {\bf U}}^2}    \right)
573      \right]_h
574\\
575= \nabla \cdot    \left( {{\begin{array}{*{20}c}   {\rm {\bf U}} \, u   \hfill \\
576                                    {\rm {\bf U}} \, v   \hfill \\
577                  \end{array} }}
578            \right)
579+\frac{1}{e_1 e_2 }     \left(
580       v\frac{\partial e_2}{\partial i}
581      -u\frac{\partial e_1}{\partial j}
582                  \right) {\rm {\bf k}} \times {\rm {\bf U}}_h
583\end{multline}
584
585The flux form has two terms, the first one is expressed as the divergence of momentum
586fluxes (hence the flux form name given to this formulation) and the second one is due to
587the curvilinear nature of the coordinate system used. The latter is called the \emph{metric}
588term and can be viewed as a modification of the Coriolis parameter:
589\begin{equation} \label{Eq_PE_cor+metric}
590f \to f + \frac{1}{e_1\;e_2}  \left(  v \frac{\partial e_2}{\partial i}
591                        -u \frac{\partial e_1}{\partial j}  \right)
592\end{equation}
593
594Note that in the case of geographical coordinate, $i.e.$ when $(i,j) \to (\lambda ,\varphi )$
595and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of
596the Coriolis parameter $f \to f+(u/a) \tan \varphi$.
597
598
599$\$\newline    % force a new ligne
600
601To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be
602written in the following tensorial formalism:
603
604\vspace{+10pt}
605$\bullet$ \textbf{Vector invariant form of the momentum equations} :
606
607\begin{subequations} \label{Eq_PE_dyn_vect}
608\begin{equation} \label{Eq_PE_dyn_vect_u} \begin{split}
609\frac{\partial u}{\partial t}
610= +   \left( {\zeta +f} \right)\,v
611   -   \frac{1}{2\,e_1}           \frac{\partial}{\partial i} \left(  u^2+v^2   \right)
612   -   \frac{1}{e_3    }  w     \frac{\partial u}{\partial k}      &      \\
613   -   \frac{1}{e_1    }            \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho _o}    \right)
614   &+   D_u^{\vect{U}}  +   F_u^{\vect{U}}      \\
615\\
616\frac{\partial v}{\partial t} =
617       -   \left( {\zeta +f} \right)\,u
618       -   \frac{1}{2\,e_2 }        \frac{\partial }{\partial j}\left(  u^2+v^\right)
619       -   \frac{1}{e_3     }   w  \frac{\partial v}{\partial k}     &      \\
620       -   \frac{1}{e_2     }        \frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}  \right)
621    &+  D_v^{\vect{U}}  +   F_v^{\vect{U}}
622\end{split} \end{equation}
623\end{subequations}
624
625
626\vspace{+10pt}
627$\bullet$ \textbf{flux form of the momentum equations} :
628\begin{subequations} \label{Eq_PE_dyn_flux}
629\begin{multline} \label{Eq_PE_dyn_flux_u}
630\frac{\partial u}{\partial t}=
631+   \left( { f + \frac{1}{e_1 \; e_2}
632               \left(    v \frac{\partial e_2}{\partial i}
633                  -u \frac{\partial e_1}{\partial j}  \right)}    \right) \, v    \\
634- \frac{1}{e_1 \; e_2}  \left(
635               \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i}
636      +        \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j}  \right)
637                 - \frac{1}{e_3 }\frac{\partial \left( {         w\,u} \right)}{\partial k}    \\
638-   \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho _o}   \right)
639+   D_u^{\vect{U}} +   F_u^{\vect{U}}
640\end{multline}
641\begin{multline} \label{Eq_PE_dyn_flux_v}
642\frac{\partial v}{\partial t}=
643-   \left( { f + \frac{1}{e_1 \; e_2}
644               \left(    v \frac{\partial e_2}{\partial i}
645                  -u \frac{\partial e_1}{\partial j}  \right)}    \right) \, u   \\
646 \frac{1}{e_1 \; e_2}   \left(
647               \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i}
648      +        \frac{\partial \left( {e_1 \,v\,v} \right)}{\partial j}  \right)
649                 - \frac{1}{e_3 } \frac{\partial \left( {        w\,v} \right)}{\partial k}    \\
650-   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}    \right)
651+  D_v^{\vect{U}} +  F_v^{\vect{U}}
652\end{multline}
653\end{subequations}
654where $\zeta$, the relative vorticity, is given by \eqref{Eq_PE_curl_Uh} and $p_s$,
655the surface pressure, is given by:
656\begin{equation} \label{Eq_PE_spg}
657p_s = \left\{ \begin{split}
659\rho \,g \,\eta &+ \rho_o \,\mu \,\frac{\partial \eta }{\partial t}      \qquad \text{ filtered     free surface}
660\end{split}
661\right.
662\end{equation}
663with $\eta$ is solution of \eqref{Eq_PE_ssh}
664
665The vertical velocity and the hydrostatic pressure are diagnosed from the following equations:
666\begin{equation} \label{Eq_w_diag}
667\frac{\partial w}{\partial k}=-\chi \;e_3
668\end{equation}
669\begin{equation} \label{Eq_hp_diag}
670\frac{\partial p_h }{\partial k}=-\rho \;g\;e_3
671\end{equation}
672where the divergence of the horizontal velocity, $\chi$ is given by \eqref{Eq_PE_div_Uh}.
673
674\vspace{+10pt}
675$\bullet$ \textit{tracer equations} :
676\begin{equation} \label{Eq_S}
677\frac{\partial T}{\partial t} =
678-\frac{1}{e_1 e_2 }\left[ {      \frac{\partial \left( {e_2 T\,u} \right)}{\partial i}
679                  +\frac{\partial \left( {e_1 T\,v} \right)}{\partial j}} \right]
680-\frac{1}{e_3 }\frac{\partial \left( {T\,w} \right)}{\partial k} + D^T + F^T
681\end{equation}
682\begin{equation} \label{Eq_T}
683\frac{\partial S}{\partial t} =
684-\frac{1}{e_1 e_2 }\left[    {\frac{\partial \left( {e_2 S\,u} \right)}{\partial i}
685                  +\frac{\partial \left( {e_1 S\,v} \right)}{\partial j}} \right]
686-\frac{1}{e_3 }\frac{\partial \left( {S\,w} \right)}{\partial k} + D^S + F^S
687\end{equation}
688\begin{equation} \label{Eq_rho}
689\rho =\rho \left( {T,S,z(k)} \right)
690\end{equation}
691
692The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale
693parameterisation used. It will be defined in \S\ref{PE_zdf}. The nature and formulation of
694${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, are discussed
695in Chapter~\ref{SBC}.
696
697
698\newpage
699$\$\newline    % force a new ligne
700% ================================================================
701% Curvilinear generalised vertical coordinate System
702% ================================================================
703\section{Curvilinear generalised vertical coordinate System}
704\label{PE_gco}
705
706The ocean domain presents a huge diversity of situation in the vertical. First the ocean surface is a time dependent surface (moving surface). Second the ocean floor depends on the geographical position, varying from more than 6,000 meters in abyssal trenches to zero at the coast. Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing.
707Therefore, in order to represent the ocean with respect to the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height $e.g.$ an $z$*-coordinate; for the second point, a space variation to fit the change of bottom topography $e.g.$ a terrain-following or $\sigma$-coordinate; and for the third point, one will be tempted to use a space and time dependent coordinate that follows the isopycnal surfaces, $e.g.$ an isopycnic coordinate.
708
709In order to satisfy two or more constrains one can even be tempted to mixed these coordinate systems, as in HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at the ocean bottom) \citep{Chassignet_al_JPO03}  or OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and $\sigma$-coordinate elsewhere) \citep{Madec_al_JPO96} among others.
710
711In fact one is totally free to choose any space and time vertical coordinate by introducing an arbitrary vertical coordinate :
712\begin{equation} \label{Eq_s}
713s=s(i,j,k,t)
714\end{equation}
715with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, when $i$, $j$ and $t$ are held fixed. \eqref{Eq_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through \eqref{Eq_s}.
716This so-called \textit{generalised vertical coordinate} \citep{Kasahara_MWR74} is in fact an Arbitrary Lagrangian--Eulerian (ALE) coordinate. Indeed, choosing an expression for $s$ is an arbitrary choice that determines which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and which part will be used to move them (Lagrangian part).
717The coordinate is also sometime referenced as an adaptive coordinate \citep{Hofmeister_al_OM09}, since the coordinate system is adapted in the course of the simulation. Its most often used implementation is via an ALE algorithm, in which a pure lagrangian step is followed by regridding and remapping steps, the later step implicitly embedding the vertical advection \citep{Hirt_al_JCP74, Chassignet_al_JPO03, White_al_JCP09}. Here we follow the \citep{Kasahara_MWR74} strategy : a regridding step (an update of the vertical coordinate) followed by an eulerian step with an explicit computation of vertical advection relative to the moving s-surfaces.
718
719%\gmcomment{
720
721%A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient...
722
723the generalized vertical coordinates used in ocean modelling are not orthogonal,
724which contrasts with many other applications in mathematical physics.
725Hence, it is useful to keep in mind the following properties that may seem
726odd on initial encounter.
727
728The horizontal velocity in ocean models measures motions in the horizontal plane,
729perpendicular to the local gravitational field. That is, horizontal velocity is mathematically
730the same regardless the vertical coordinate, be it geopotential, isopycnal, pressure,
731or terrain following. The key motivation for maintaining the same horizontal velocity
732component is that the hydrostatic and geostrophic balances are dominant in the large-scale ocean.
733Use of an alternative quasi-horizontal velocity, for example one oriented parallel
734to the generalized surface, would lead to unacceptable numerical errors.
735Correspondingly, the vertical direction is anti-parallel to the gravitational force in all
736of the coordinate systems. We do not choose the alternative of a quasi-vertical
737direction oriented normal to the surface of a constant generalized vertical coordinate.
738
739It is the method used to measure transport across the generalized vertical coordinate
740surfaces which differs between the vertical coordinate choices. That is, computation
741of the dia-surface velocity component represents the fundamental distinction between
742the various coordinates. In some models, such as geopotential, pressure, and
743terrain following, this transport is typically diagnosed from volume or mass conservation.
744In other models, such as isopycnal layered models, this transport is prescribed based
745on assumptions about the physical processes producing a flux across the layer interfaces.
746
747
748In this section we first establish the PE in the generalised vertical $s$-coordinate,
749then we discuss the particular cases available in \NEMO, namely $z$, $z$*, $s$, and $\tilde z$
750%}
751
752% -------------------------------------------------------------------------------------------------------------
753% The s-coordinate Formulation
754% -------------------------------------------------------------------------------------------------------------
755\subsection{The \textit{s-}coordinate Formulation}
756
757Starting from the set of equations established in \S\ref{PE_zco} for the special case $k=z$
758and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, which includes
759$z$-, \textit{z*}- and $\sigma-$coordinates as special cases ($s=z$, $s=\textit{z*}$, and
760$s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). A formal derivation of the transformed
761equations is given in Appendix~\ref{Apdx_A}. Let us define the vertical scale factor by
762$e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), and the slopes in the
763(\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by :
764\begin{equation} \label{Eq_PE_sco_slope}
765\sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s
767\sigma _2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s
768\end{equation}
769We also introduce  $\omega$, a dia-surface velocity component, defined as the velocity
770relative to the moving $s$-surfaces and normal to them:
771\begin{equation} \label{Eq_PE_sco_w}
772\omega  = w - e_3 \, \frac{\partial s}{\partial t} - \sigma _1 \,u - \sigma _2 \,v    \\
773\end{equation}
774
775The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows:
776
777 \vspace{0.5cm}
778* momentum equation:
779\begin{multline} \label{Eq_PE_sco_u}
780\frac{1}{e_3} \frac{\partial \left(  e_3\,\right) }{\partial t}=
781   +   \left( {\zeta +f} \right)\,v
782   -   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right)
783   -   \frac{1}{e_3} \omega \frac{\partial u}{\partial k}       \\
784   -   \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho _o}    \right)
785   +  g\frac{\rho }{\rho _o}\sigma _1
786   +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad
787\end{multline}
788\begin{multline} \label{Eq_PE_sco_v}
789\frac{1}{e_3} \frac{\partial \left(  e_3\,\right) }{\partial t}=
790   -   \left( {\zeta +f} \right)\,u
791   -   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^\right)
792   -   \frac{1}{e_3 } \omega \frac{\partial v}{\partial k}         \\
793   -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}  \right)
794    +  g\frac{\rho }{\rho _o }\sigma _2
795   +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad
796\end{multline}
797where the relative vorticity, \textit{$\zeta$}, the surface pressure gradient, and the hydrostatic
798pressure have the same expressions as in $z$-coordinates although they do not represent
799exactly the same quantities. $\omega$ is provided by the continuity equation
800(see Appendix~\ref{Apdx_A}):
801
802\begin{equation} \label{Eq_PE_sco_continuity}
803\frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0
805\chi =\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3 \,u}
806\right)}{\partial i}+\frac{\partial \left( {e_1 e_3 \,v} \right)}{\partial
807j}} \right]
808\end{equation}
809
810 \vspace{0.5cm}
811* tracer equations:
812\begin{multline} \label{Eq_PE_sco_t}
813\frac{1}{e_3} \frac{\partial \left(  e_3\,\right) }{\partial t}=
814-\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,T} \right)}{\partial i}
815                                           +\frac{\partial \left( {e_1 e_3\,v\,T} \right)}{\partial j}} \right]   \\
816-\frac{1}{e_3 }\frac{\partial \left( {T\,\omega } \right)}{\partial k}   + D^T + F^S   \qquad
817\end{multline}
818
819\begin{multline} \label{Eq_PE_sco_s}
820\frac{1}{e_3} \frac{\partial \left(  e_3\,\right) }{\partial t}=
821-\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,S} \right)}{\partial i}
822                                           +\frac{\partial \left( {e_1 e_3\,v\,S} \right)}{\partial j}} \right]    \\
823-\frac{1}{e_3 }\frac{\partial \left( {S\,\omega } \right)}{\partial k}     + D^S + F^S   \qquad
824\end{multline}
825
826The equation of state has the same expression as in $z$-coordinate, and similar expressions
827are used for mixing and forcing terms.
828
829\gmcomment{
830\colorbox{yellow}{ to be updated $= = >$}
831Add a few works on z and zps and s and underlies the differences between all of them
832\colorbox{yellow}{ $< = =$ end update}  }
833
834
835
836% -------------------------------------------------------------------------------------------------------------
837% Curvilinear z*-coordinate System
838% -------------------------------------------------------------------------------------------------------------
839\subsection{Curvilinear \textit{z*}--coordinate System}
840\label{PE_zco_star}
841
842%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
843\begin{figure}[!b]    \begin{center}
844\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zstar.pdf}
845\caption{   \label{Fig_z_zstar}
846(a) $z$-coordinate in linear free-surface case ;
847(b) $z-$coordinate in non-linear free surface case ;
848(c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate
850\end{center}   \end{figure}
851%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
852
853
854In that case, the free surface equation is nonlinear, and the variations of volume are fully
855taken into account. These coordinates systems is presented in a report \citep{Levier2007}
856available on the \NEMO web site.
857
858%\gmcomment{
859The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation
860which allows one to deal with large amplitude free-surface
861variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In
862the  \textit{z*} formulation, the variation of the column thickness due to sea-surface
863undulations is not concentrated in the surface level, as in the $z$-coordinate formulation,
864but is equally distributed over the full water column. Thus vertical
865levels naturally follow sea-surface variations, with a linear attenuation with
866depth, as illustrated by figure fig.1c . Note that with a flat bottom, such as in
867fig.1c, the bottom-following  $z$ coordinate and  \textit{z*} are equivalent.
868The definition and modified oceanic equations for the rescaled vertical coordinate
869 \textit{z*}, including the treatment of fresh-water flux at the surface, are
870detailed in Adcroft and Campin (2004). The major points are summarized
871here. The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as:
872\begin{equation} \label{Eq_z-star}
873H +  \textit{z*} = (H + z) / r \quad \text{and} \ \delta \textit{z*} = \delta z / r \quad \text{with} \ r = \frac{H+\eta} {H}
874\end{equation}
875Since the vertical displacement of the free surface is incorporated in the vertical
876coordinate  \textit{z*}, the upper and lower boundaries are at fixed  \textit{z*} position,
877$\textit{z*} = 0$ and  $\textit{z*} = -H$ respectively. Also the divergence of the flow field
878is no longer zero as shown by the continuity equation:
879\begin{equation*}
880\frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right)
881      \left( r \; w\textit{*} \right) = 0
882\end{equation*}
883%}
884
885
886% from MOM4p1 documentation
887
888To overcome problems with vanishing surface and/or bottom cells, we consider the
889zstar coordinate
890\begin{equation} \label{PE_}
891   z^\star = H \left( \frac{z-\eta}{H+\eta} \right)
892\end{equation}
893
894This coordinate is closely related to the "eta" coordinate used in many atmospheric
895models (see Black (1994) for a review of eta coordinate atmospheric models). It
896was originally used in ocean models by Stacey et al. (1995) for studies of tides
897next to shelves, and it has been recently promoted by Adcroft and Campin (2004)
898for global climate modelling.
899
900The surfaces of constant $z^\star$ are quasi-horizontal. Indeed, the $z^\star$ coordinate reduces to $z$ when $\eta$ is zero. In general, when noting the large differences between
901undulations of the bottom topography versus undulations in the surface height, it
902is clear that surfaces constant $z^\star$ are very similar to the depth surfaces. These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to terrain following sigma models discussed in \S\ref{PE_sco}.
903Additionally, since $z^\star$ when $\eta = 0$, no flow is spontaneously generated in an
904unforced ocean starting from rest, regardless the bottom topography. This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in
905the presence of nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, depending on the sophistication of the pressure
906gradient solver. The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of neutral physics parameterizations in $z^\star$ models using
907the same techniques as in $z$-models (see Chapters 13-16 of \cite{Griffies_Bk04}) for a
908discussion of neutral physics in $z$-models, as well as Section \S\ref{LDF_slp}
909in this document for treatment in \NEMO).
910
911The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$. Hence, all
912cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$. This
913is a minor constraint relative to that encountered on the surface height when using
914$s = z$ or $s = z - \eta$.
915
916Because $z^\star$ has a time independent range, all grid cells have static increments
917ds, and the sum of the ver tical increments yields the time independent ocean
918depth %·k ds = H.
919The $z^\star$ coordinate is therefore invisible to undulations of the
920free surface, since it moves along with the free surface. This proper ty means that
921no spurious ver tical transpor t is induced across surfaces of constant $z^\star$ by the
922motion of external gravity waves. Such spurious transpor t can be a problem in
923z-models, especially those with tidal forcing. Quite generally, the time independent
924range for the $z^\star$ coordinate is a very convenient proper ty that allows for a nearly
925arbitrary ver tical resolution even in the presence of large amplitude fluctuations of
926the surface height, again so long as $\eta > -H$.
927
928%end MOM doc %%%
929
930
931
932\newpage
933% -------------------------------------------------------------------------------------------------------------
934% Terrain following  coordinate System
935% -------------------------------------------------------------------------------------------------------------
936\subsection{Curvilinear Terrain-following \textit{s}--coordinate}
937\label{PE_sco}
938
939% -------------------------------------------------------------------------------------------------------------
940% Introduction
941% -------------------------------------------------------------------------------------------------------------
942\subsubsection{Introduction}
943
944Several important aspects of the ocean circulation are influenced by bottom topography.
945Of course, the most important is that bottom topography determines deep ocean sub-basins,
946barriers, sills and channels that strongly constrain the path of water masses, but more subtle
947effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary
948one along continental slopes. Topographic Rossby waves can be excited and can interact
949with the mean current. In the $z-$coordinate system presented in the previous section
950(\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is
951discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom
952and to large localized depth gradients associated with large localized vertical velocities.
953The response to such a velocity field often leads to numerical dispersion effects.
954One solution to strongly reduce this error is to use a partial step representation of bottom
956Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate)
957
958The $s$-coordinate avoids the discretisation error in the depth field since the layers of
959computation are gradually adjusted with depth to the ocean bottom. Relatively small
960topographic features as well as  gentle, large-scale slopes of the sea floor in the deep
961ocean, which would be ignored in typical $z$-model applications with the largest grid
962spacing at greatest depths, can easily be represented (with relatively low vertical resolution).
963A terrain-following model (hereafter $s-$model) also facilitates the modelling of the
964boundary layer flows over a large depth range, which in the framework of the $z$-model
965would require high vertical resolution over the whole depth range. Moreover, with a
966$s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface
967as the only boundaries of the domain (nomore lateral boundary condition to specify).
968Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a
969homogeneous ocean, it has strong limitations as soon as stratification is introduced.
970The main two problems come from the truncation error in the horizontal pressure
971gradient and a possibly increased diapycnal diffusion. The horizontal pressure force
972in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}),
973
974\begin{equation} \label{Eq_PE_p_sco}
975\left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial
976p}{\partial s}\left. {\nabla z} \right|_s
977\end{equation}
978
979The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface
980and introduces a truncation error that is not present in a $z$-model. In the special case
981of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$),
982\citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude
983of this truncation error. It depends on topographic slope, stratification, horizontal and
984vertical resolution, the equation of state, and the finite difference scheme. This error
985limits the possible topographic slopes that a model can handle at a given horizontal
986and vertical resolution. This is a severe restriction for large-scale applications using
987realistic bottom topography. The large-scale slopes require high horizontal resolution,
988and the computational cost becomes prohibitive. This problem can be at least partially
989overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec_al_JPO96}. However, the definition of the model
990domain vertical coordinate becomes then a non-trivial thing for a realistic bottom
991topography: a envelope topography is defined in $s$-coordinate on which a full or
992partial step bottom topography is then applied in order to adjust the model depth to
993the observed one (see \S\ref{DOM_zgr}.
994
995For numerical reasons a minimum of diffusion is required along the coordinate surfaces
996of any finite difference model. It causes spurious diapycnal mixing when coordinate
997surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as
998well as for a $s$-model. However, density varies more strongly on $s-$surfaces than
999on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal
1000diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a
1001$z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal
1002circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation.
1003For example, imagine an isolated bump of topography in an ocean at rest with a horizontally
1004uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral
1005surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast,
1006the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column
1007($i.e.$ the main thermocline) \citep{Madec_al_JPO96}. An alternate solution consists of rotating
1008the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}.
1009Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large,
1010strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}).
1011
1012The $s-$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two
1013aspects from similar models:  it allows  a representation of bottom topography with mixed
1014full or partial step-like/terrain following topography ; It also offers a completely general
1015transformation, $s=s(i,j,z)$ for the vertical coordinate.
1016
1017
1018\newpage
1019% -------------------------------------------------------------------------------------------------------------
1020% Curvilinear z-tilde coordinate System
1021% -------------------------------------------------------------------------------------------------------------
1022\subsection{Curvilinear $\tilde{z}$--coordinate}
1023\label{PE_zco_tilde}
1024
1025The $\tilde{z}$-coordinate has been developed by \citet{Leclair_Madec_OM10s}.
1026It is not available in the current version of \NEMO.
1027
1028\newpage
1029% ================================================================
1030% Subgrid Scale Physics
1031% ================================================================
1032\section{Subgrid Scale Physics}
1033\label{PE_zdf_ldf}
1034
1035The primitive equations describe the behaviour of a geophysical fluid at
1036space and time scales larger than a few kilometres in the horizontal, a few
1037meters in the vertical and a few minutes. They are usually solved at larger
1038scales: the specified grid spacing and time step of the numerical model. The
1039effects of smaller scale motions (coming from the advective terms in the
1040Navier-Stokes equations) must be represented entirely in terms of
1041large-scale patterns to close the equations. These effects appear in the
1042equations as the divergence of turbulent fluxes ($i.e.$ fluxes associated with
1043the mean correlation of small scale perturbations). Assuming a turbulent
1044closure hypothesis is equivalent to choose a formulation for these fluxes.
1045It is usually called the subgrid scale physics. It must be emphasized that
1046this is the weakest part of the primitive equations, but also one of the
1047most important for long-term simulations as small scale processes \textit{in fine}
1048balance the surface input of kinetic energy and heat.
1049
1050The control exerted by gravity on the flow induces a strong anisotropy
1051between the lateral and vertical motions. Therefore subgrid-scale physics
1052\textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in \eqref{Eq_PE_dyn},
1053\eqref{Eq_PE_tra_T} and \eqref{Eq_PE_tra_S} are divided into a lateral part
1054\textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part
1055\textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. The formulation of these terms
1056and their underlying physics are briefly discussed in the next two subsections.
1057
1058% -------------------------------------------------------------------------------------------------------------
1059% Vertical Subgrid Scale Physics
1060% -------------------------------------------------------------------------------------------------------------
1061\subsection{Vertical Subgrid Scale Physics}
1062\label{PE_zdf}
1063
1064The model resolution is always larger than the scale at which the major
1065sources of vertical turbulence occur (shear instability, internal wave
1066breaking...). Turbulent motions are thus never explicitly solved, even
1067partially, but always parameterized. The vertical turbulent fluxes are
1068assumed to depend linearly on the gradients of large-scale quantities (for
1069example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$,
1070where $A^{vT}$ is an eddy coefficient). This formulation is
1071analogous to that of molecular diffusion and dissipation. This is quite
1072clearly a necessary compromise: considering only the molecular viscosity
1073acting on large scale severely underestimates the role of turbulent
1074diffusion and dissipation, while an accurate consideration of the details of
1075turbulent motions is simply impractical. The resulting vertical momentum and
1076tracer diffusive operators are of second order:
1077\begin{equation} \label{Eq_PE_zdf}
1078   \begin{split}
1079{\vect{D}}^{v \vect{U}} &=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\
1080D^{vT}                         &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ ,
1082D^{vS}=\frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial S}{\partial z}} \right)
1083   \end{split}
1084\end{equation}
1085where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients,
1086respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat
1087and salt must be specified (see Chap.~\ref{SBC} and \ref{ZDF} and \S\ref{TRA_bbl}).
1088All the vertical physics is embedded in the specification of the eddy coefficients.
1089They can be assumed to be either constant, or function of the local fluid properties
1090($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a
1091turbulent closure model. The choices available in \NEMO are discussed in \S\ref{ZDF}).
1092
1093% -------------------------------------------------------------------------------------------------------------
1094% Lateral Diffusive and Viscous Operators Formulation
1095% -------------------------------------------------------------------------------------------------------------
1096\subsection{Formulation of the Lateral Diffusive and Viscous Operators}
1097\label{PE_ldf}
1098
1099Lateral turbulence can be roughly divided into a mesoscale turbulence
1100associated with eddies (which can be solved explicitly if the resolution is
1101sufficient since their underlying physics are included in the primitive
1102equations), and a sub mesoscale turbulence which is never explicitly solved
1103even partially, but always parameterized. The formulation of lateral eddy
1104fluxes depends on whether the mesoscale is below or above the grid-spacing
1105($i.e.$ the model is eddy-resolving or not).
1106
1107In non-eddy-resolving configurations, the closure is similar to that used
1108for the vertical physics. The lateral turbulent fluxes are assumed to depend
1109linearly on the lateral gradients of large-scale quantities. The resulting
1110lateral diffusive and dissipative operators are of second order.
1111Observations show that lateral mixing induced by mesoscale turbulence tends
1112to be along isopycnal surfaces (or more precisely neutral surfaces \cite{McDougall1987})
1113rather than across them.
1114As the slope of neutral surfaces is small in the ocean, a common
1115approximation is to assume that the `lateral' direction is the horizontal,
1116$i.e.$ the lateral mixing is performed along geopotential surfaces. This leads
1117to a geopotential second order operator for lateral subgrid scale physics.
1118This assumption can be relaxed: the eddy-induced turbulent fluxes can be
1119better approached by assuming that they depend linearly on the gradients of
1120large-scale quantities computed along neutral surfaces. In such a case,
1121the diffusive operator is an isoneutral second order operator and it has
1122components in the three space directions. However, both horizontal and
1123isoneutral operators have no effect on mean ($i.e.$ large scale) potential
1124energy whereas potential energy is a main source of turbulence (through
1125baroclinic instabilities). \citet{Gent1990} have proposed a
1126parameterisation of mesoscale eddy-induced turbulence which associates an
1127eddy-induced velocity to the isoneutral diffusion. Its mean effect is to
1128reduce the mean potential energy of the ocean. This leads to a formulation
1129of lateral subgrid-scale physics made up of an isoneutral second order
1130operator and an eddy induced advective part. In all these lateral diffusive
1131formulations, the specification of the lateral eddy coefficients remains the
1132problematic point as there is no really satisfactory formulation of these
1133coefficients as a function of large-scale features.
1134
1135In eddy-resolving configurations, a second order operator can be used, but
1136usually the more scale selective biharmonic operator is preferred as the
1137grid-spacing is usually not small enough compared to the scale of the
1138eddies. The role devoted to the subgrid-scale physics is to dissipate the
1139energy that cascades toward the grid scale and thus to ensure the stability of
1140the model while not interfering with the resolved mesoscale activity. Another approach
1141is becoming more and more popular: instead of specifying explicitly a sub-grid scale
1142term in the momentum and tracer time evolution equations, one uses a advective
1143scheme which is diffusive enough to maintain the model stability. It must be emphasised
1144that then, all the sub-grid scale physics is included in the formulation of the
1146
1147All these parameterisations of subgrid scale physics have advantages and
1148drawbacks. There are not all available in \NEMO. In the $z$-coordinate
1149formulation, five options are offered for active tracers (temperature and
1150salinity): second order geopotential operator, second order isoneutral
1151operator, \citet{Gent1990} parameterisation, fourth order
1152geopotential operator, and various slightly diffusive advection schemes.
1153The same options are available for momentum, except
1154\citet{Gent1990} parameterisation which only involves tracers. In the
1155$s$-coordinate formulation, additional options are offered for tracers: second
1156order operator acting along $s-$surfaces, and for momentum: fourth order
1157operator acting along $s-$surfaces (see \S\ref{LDF}).
1158
1159\subsubsection{Lateral second order tracer diffusive operator}
1160
1161The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}):
1162\begin{equation} \label{Eq_PE_iso_tensor}
1163D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
1165 1 \hfill & 0 \hfill & {-r_1 } \hfill \\
1166 0 \hfill & 1 \hfill & {-r_2 } \hfill \\
1167 {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
1168\end{array} }} \right)
1169\end{equation}
1170where $r_1 \;\mbox{and}\;r_2$ are the slopes between the surface along
1171which the diffusive operator acts and the model level ($e. g.$ $z$- or
1172$s$-surfaces). Note that the formulation \eqref{Eq_PE_iso_tensor} is exact for the
1173rotation between geopotential and $s$-surfaces, while it is only an approximation
1174for the rotation between isoneutral and $z$- or $s$-surfaces. Indeed, in the latter
1175case, two assumptions are made to simplify  \eqref{Eq_PE_iso_tensor} \citep{Cox1987}.
1176First, the horizontal contribution of the dianeutral mixing is neglected since the ratio
1177between iso and dia-neutral diffusive coefficients is known to be several orders of
1178magnitude smaller than unity. Second, the two isoneutral directions of diffusion are
1179assumed to be independent since the slopes are generally less than $10^{-2}$ in the
1180ocean (see Appendix~\ref{Apdx_B}).
1181
1182For \textit{geopotential} diffusion, $r_1$ and $r_2$ are the slopes between the
1183geopotential and computational surfaces: in $z$-coordinates they are zero
1184($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are
1185equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ).
1186
1187For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral
1188and computational surfaces. Therefore, they are different quantities,
1189but have similar expressions in $z$- and $s$-coordinates. In $z$-coordinates:
1190\begin{equation} \label{Eq_PE_iso_slopes}
1191r_1 =\frac{e_3 }{e_1 }  \left( {\frac{\partial \rho }{\partial i}} \right)
1192                  \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \ , \quad
1193r_1 =\frac{e_3 }{e_1 }  \left( {\frac{\partial \rho }{\partial i}} \right)
1194                  \left( {\frac{\partial \rho }{\partial k}} \right)^{-1},
1195\end{equation}
1196while in $s$-coordinates $\partial/\partial k$ is replaced by
1197$\partial/\partial s$.
1198
1199\subsubsection{Eddy induced velocity}
1200 When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used,
1201an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers:
1202\begin{equation} \label{Eq_PE_iso+eiv}
1203D^{lT}=\nabla \cdot \left( {A^{lT}\;\Re \;\nabla T} \right)
1204           +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right)
1205\end{equation}
1206where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent,
1207eddy-induced transport velocity. This velocity field is defined by:
1208\begin{equation} \label{Eq_PE_eiv}
1209   \begin{split}
1210 u^\ast  &= +\frac{1}{e_3       }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_1 } \right] \\
1211 v^\ast  &= +\frac{1}{e_3       }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_2 } \right] \\
1212 w^\ast &=  -\frac{1}{e_1 e_2 }\left[
1213                      \frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right)
1214                    +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right)      \right]
1215   \end{split}
1216\end{equation}
1217where $A^{eiv}$ is the eddy induced velocity coefficient (or equivalently the isoneutral
1218thickness diffusivity coefficient), and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes
1219between isoneutral and \emph{geopotential} surfaces. Their values are
1220thus independent of the vertical coordinate, but their expression depends on the coordinate:
1221\begin{align} \label{Eq_PE_slopes_eiv}
1222\tilde{r}_n = \begin{cases}
1223   r_n                  &      \text{in $z$-coordinate}    \\
1224   r_n + \sigma_n &      \text{in \textit{z*} and $s$-coordinates}
1225                   \end{cases}
1227\end{align}
1228
1229The normal component of the eddy induced velocity is zero at all the boundaries.
1230This can be achieved in a model by tapering either the eddy coefficient or the slopes
1231to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}).
1232
1233\subsubsection{Lateral fourth order tracer diffusive operator}
1234
1235The lateral fourth order tracer diffusive operator is defined by:
1236\begin{equation} \label{Eq_PE_bilapT}
1237D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right)
1238\qquad \text{where} \  D^{lT}=\Delta \left( {A^{lT}\;\Delta T} \right)
1239 \end{equation}
1240
1241It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with
1242the eddy diffusion coefficient correctly placed.
1243
1244
1245\subsubsection{Lateral second order momentum diffusive operator}
1246
1247The second order momentum diffusive operator along $z$- or $s$-surfaces is found by
1248applying \eqref{Eq_PE_lap_vector} to the horizontal velocity vector (see Appendix~\ref{Apdx_B}):
1249\begin{equation} \label{Eq_PE_lapU}
1250\begin{split}
1251{\rm {\bf D}}^{l{\rm {\bf U}}}
1252&= \quad \  \nabla _h \left( {A^{lm}\chi } \right)
1253   \ - \ \nabla _h \times \left( {A^{lm}\,\zeta \;{\rm {\bf k}}} \right)     \\
1254&=   \left(      \begin{aligned}
1255             \frac{1}{e_1      } \frac{\partial \left( A^{lm} \chi          \right)}{\partial i}
1256         &-\frac{1}{e_2 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial j}  \\
1257             \frac{1}{e_2      }\frac{\partial \left( {A^{lm} \chi         } \right)}{\partial j}
1258         &+\frac{1}{e_1 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial i}
1259        \end{aligned}    \right)
1260\end{split}
1261\end{equation}
1262
1263Such a formulation ensures a complete separation between the vorticity and
1264horizontal divergence fields (see Appendix~\ref{Apdx_C}). Unfortunately, it is not
1265available for geopotential diffusion in $s-$coordinates and for isoneutral
1266diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required).
1267In these two cases, the $u$ and $v-$fields are considered as independent scalar
1268fields, so that the diffusive operator is given by:
1269\begin{equation} \label{Eq_PE_lapU_iso}
1270\begin{split}
1271 D_u^{l{\rm {\bf U}}} &= \nabla .\left( {\Re \;\nabla u} \right) \\
1272 D_v^{l{\rm {\bf U}}} &= \nabla .\left( {\Re \;\nabla v} \right)
1273 \end{split}
1274 \end{equation}
1275where $\Re$ is given by  \eqref{Eq_PE_iso_tensor}. It is the same expression as
1276those used for diffusive operator on tracers. It must be emphasised that such a
1277formulation is only exact in a Cartesian coordinate system, $i.e.$ on a $f-$ or
1278$\beta-$plane, not on the sphere. It is also a very good approximation in vicinity
1279of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}.
1280
1281\subsubsection{lateral fourth order momentum diffusive operator}
1282
1283As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces
1284is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU}
1285with the eddy viscosity coefficient correctly placed:
1286
1287geopotential diffusion in $z$-coordinate:
1288\begin{equation} \label{Eq_PE_bilapU}
1289\begin{split}
1290{\rm {\bf D}}^{l{\rm {\bf U}}} &=\nabla _h \left\{ {\;\nabla _h {\rm {\bf
1291.}}\left[ {A^{lm}\,\nabla _h \left( \chi \right)} \right]\;}
1292\right\}\;   \\
1293&+\nabla _h \times \left\{ {\;{\rm {\bf k}}\cdot \nabla \times
1294\left[ {A^{lm}\,\nabla _h \times \left( {\zeta \;{\rm {\bf k}}} \right)}
1295\right]\;} \right\}
1296\end{split}
1297\end{equation}
1298
1299\gmcomment{  change the position of the coefficient, both here and in the code}
1300
1301geopotential diffusion in $s$-coordinate:
1302\begin{equation} \label{Eq_bilapU_iso}
1303   \left\{   \begin{aligned}
1304         D_u^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta u} \right) \\
1305         D_v^{l{\rm {\bf U}}} =\Delta \left( {A^{lm}\;\Delta v} \right)
1306   \end{aligned}    \right.