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

Last change on this file since 11335 was 11335, checked in by mikebell, 2 years ago

review of chap_model_basics, annex_A and annex_B

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