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

Last change on this file since 11690 was 11690, checked in by nicolasmartin, 13 months ago

Update chapters according to previous commit

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