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 NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex @ 10406

Last change on this file since 10406 was 10406, checked in by nicolasmartin, 5 years ago

Edition of math environments

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