# New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Chap_Model_Basics.tex in branches/2015/dev_r5003_MERCATOR6_CRS/DOC/TexFiles/Chapters – NEMO

# source:branches/2015/dev_r5003_MERCATOR6_CRS/DOC/TexFiles/Chapters/Chap_Model_Basics.tex@7260

Last change on this file since 7260 was 7260, checked in by cbricaud, 7 years ago

phaze DOC/ directory of the CRS branch with nemo_v3_6_STABLE branch at rev 7213 (09-09-2016) (merge -r 5519:7213 )

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