Changeset 10501


Ignore:
Timestamp:
2019-01-10T17:30:42+01:00 (15 months ago)
Author:
nicolasmartin
Message:

Global work on math environnements for equations (partial commits)

  • Enhance chosen environnement regarding alignement and numbering needs (alignat, gather and respective alignedat and gathered)
  • Simplify the syntax by using new commands:
    • \pd[…]{…} for \frac{\partial …}{\partial …}
    • \lt-\rt in place of \right-\left for automatic sizing of delimiters
    • \vect … for {\rm {\bf …}}} and similar ones
  • Align continuation lines to improve code readability
Location:
NEMO/trunk/doc/latex/NEMO/subfiles
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex

    r10442 r10501  
    22 
    33\begin{document} 
     4 
    45% ================================================================ 
    5 % Chapter 1 Model Basics 
     6% Chapter 1  Model Basics 
    67% ================================================================ 
    7  
    88\chapter{Model Basics} 
    99\label{chap:PE} 
    10  
    1110\minitoc 
    1211 
     
    2625\label{subsec:PE_Vector} 
    2726 
    28  
    2927The ocean is a fluid that can be described to a good approximation by the primitive equations, 
    3028\ie the Navier-Stokes equations along with a nonlinear equation of state which 
     
    3230plus the following additional assumptions made from scale considerations: 
    3331 
    34 \textit{(1) spherical earth approximation:} the geopotential surfaces are assumed to be spheres so that 
    35 gravity (local vertical) is parallel to the earth's radius 
    36  
    37 \textit{(2) thin-shell approximation:} the ocean depth is neglected compared to the earth's radius 
    38  
    39 \textit{(3) turbulent closure hypothesis:} the turbulent fluxes 
    40 (which represent the effect of small scale processes on the large-scale) are expressed in terms of 
    41 large-scale features 
    42  
    43 \textit{(4) Boussinesq hypothesis:} density variations are neglected except in their contribution to 
    44 the buoyancy force 
    45  
    46 \textit{(5) Hydrostatic hypothesis:} the vertical momentum equation is reduced to a balance between 
    47 the vertical pressure gradient and the buoyancy force 
    48 (this removes convective processes from the initial Navier-Stokes equations and so 
    49 convective processes must be parameterized instead) 
    50  
    51 \textit{(6) Incompressibility hypothesis:} the three dimensional divergence of the velocity vector is assumed to 
    52 be zero. 
     32\begin{enumerate} 
     33\item 
     34  \textit{spherical earth approximation}: the geopotential surfaces are assumed to be spheres so that 
     35  gravity (local vertical) is parallel to the earth's radius 
     36\item 
     37  \textit{thin-shell approximation}: the ocean depth is neglected compared to the earth's radius 
     38\item 
     39  \textit{turbulent closure hypothesis}: the turbulent fluxes 
     40  (which represent the effect of small scale processes on the large-scale) 
     41  are expressed in terms of large-scale features 
     42\item 
     43  \textit{Boussinesq hypothesis}: density variations are neglected except in their contribution to 
     44  the buoyancy force 
     45  \begin{equation} 
     46    \label{eq:PE_eos} 
     47    \rho = \rho \ (T,S,p) 
     48  \end{equation} 
     49\item 
     50  \textit{Hydrostatic hypothesis}: the vertical momentum equation is reduced to a balance between 
     51  the vertical pressure gradient and the buoyancy force 
     52  (this removes convective processes from the initial Navier-Stokes equations and so 
     53  convective processes must be parameterized instead) 
     54  \begin{equation} 
     55    \label{eq:PE_hydrostatic} 
     56    \pd[p]{z} = - \rho \ g 
     57  \end{equation} 
     58\item 
     59  \textit{Incompressibility hypothesis}: the three dimensional divergence of the velocity vector $\vect U$ 
     60  is assumed to be zero. 
     61  \begin{equation} 
     62    \label{eq:PE_continuity} 
     63    \nabla \cdot \vect U = 0 
     64  \end{equation} 
     65\end{enumerate} 
    5366 
    5467Because the gravitational force is so dominant in the equations of large-scale motions, 
    55 it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to 
    56 the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to 
    57 \textbf{k}, \ie tangent to the geopotential surfaces. 
    58 Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$  
    59 (the subscript $h$ denotes the local horizontal vector, \ie over the (\textbf{i},\textbf{j}) plane),  
    60 $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. 
    61 The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system 
    62 provides the following six equations 
    63 (namely the momentum balance, the hydrostatic equilibrium, the incompressibility equation, 
    64 the heat and salt conservation equations and an equation of state): 
     68it is useful to choose an orthogonal set of unit vectors $(i,j,k)$ linked to the earth such that 
     69$k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 
     70\ie tangent to the geopotential surfaces. 
     71Let us define the following variables: $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ 
     72(the subscript $h$ denotes the local horizontal vector, \ie over the $(i,j)$ plane),  
     73$T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density. 
     74The vector invariant form of the primitive equations in the $(i,j,k)$ vector system provides 
     75the following equations: 
    6576\begin{subequations} 
    6677  \label{eq:PE} 
    67   \begin{equation} 
     78  \begin{gather} 
     79    \intertext{$-$ the momentum balance} 
    6880    \label{eq:PE_dyn} 
    69     \frac{\partial {\rm {\bf U}}_h }{\partial t}= 
    70     -\left[    {\left( {\nabla \times {\rm {\bf U}}} \right)\times {\rm {\bf U}} 
    71         +\frac{1}{2}\nabla \left( {{\rm {\bf U}}^2} \right)}    \right]_h 
    72     -f\;{\rm {\bf k}}\times {\rm {\bf U}}_h 
    73     -\frac{1}{\rho_o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} + {\rm {\bf F}}^{\rm {\bf U}} 
    74   \end{equation} 
    75   \begin{equation} 
    76     \label{eq:PE_hydrostatic} 
    77     \frac{\partial p }{\partial z} = - \rho \ g 
    78   \end{equation} 
    79   \begin{equation} 
    80     \label{eq:PE_continuity} 
    81     \nabla \cdot {\bf U}=  0 
    82   \end{equation} 
    83   \begin{equation} 
     81    \pd[\vect U_h]{t} = - \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h 
     82                        - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p 
     83                        + \vect D^{\vect U} + \vect F^{\vect U} \\ 
     84    \intertext{$-$ the heat and salt conservation equations} 
    8485    \label{eq:PE_tra_T} 
    85     \frac{\partial T}{\partial t} = - \nabla \cdot  \left( T \ \rm{\bf U} \right) + D^T + F^T 
    86   \end{equation} 
    87   \begin{equation} 
     86    \pd[T]{t} = - \nabla \cdot (T \ \vect U) + D^T + F^T \\ 
    8887    \label{eq:PE_tra_S} 
    89     \frac{\partial S}{\partial t} = - \nabla \cdot  \left( S \ \rm{\bf U} \right) + D^S + F^S 
    90   \end{equation} 
    91   \begin{equation} 
    92     \label{eq:PE_eos} 
    93     \rho = \rho \left( T,S,p \right) 
    94   \end{equation} 
     88    \pd[S]{t} = - \nabla \cdot (S \ \vect U) + D^S + F^S 
     89  \end{gather} 
    9590\end{subequations} 
    96 where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, $t$ is the time, 
    97 $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by the equation of state 
     91where $\nabla$ is the generalised derivative vector operator in $(i,j,k)$ directions, $t$ is the time, 
     92$z$ is the vertical coordinate, $\rho$ is the \textit{in situ} density given by the equation of state 
    9893(\autoref{eq:PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 
    99 $f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration 
    100 (where $\bf \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 
    101 ${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, 
    102 temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$ surface forcing terms. 
     94$f = 2 \vect \Omega \cdot k$ is the Coriolis acceleration 
     95(where $\vect \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 
     96$\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, 
     97temperature and salinity, and $\vect F^{\vect U}$, $F^T$ and $F^S$ surface forcing terms. 
    10398Their nature and formulation are discussed in \autoref{sec:PE_zdf_ldf} and \autoref{subsec:PE_boundary_condition}. 
    104  
    105  
    10699 
    107100% ------------------------------------------------------------------------------------------------------------- 
     
    113106An ocean is bounded by complex coastlines, bottom topography at its base and 
    114107an air-sea or ice-sea interface at its top. 
    115 These boundaries can be defined by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, 
     108These boundaries can be defined by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,k,t)$, 
    116109where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface. 
    117 Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, chosen as a mean sea surface 
     110Both $H$ and $\eta$ are usually referenced to a given surface, $z = 0$, chosen as a mean sea surface 
    118111(\autoref{fig:ocean_bc}). 
    119112Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with 
     
    127120\begin{figure}[!ht] 
    128121  \begin{center} 
    129     \includegraphics[width=0.90\textwidth]{Fig_I_ocean_bc} 
    130     \caption{   \protect\label{fig:ocean_bc} 
    131       The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$, 
     122    \includegraphics[]{Fig_I_ocean_bc} 
     123    \caption{ 
     124      \protect\label{fig:ocean_bc} 
     125      The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, 
    132126      where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. 
    133       Both $H$ and $\eta$ are referenced to $z=0$. 
     127      Both $H$ and $\eta$ are referenced to $z = 0$. 
    134128    } 
    135129  \end{center} 
    136130\end{figure} 
    137131%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    138  
    139132 
    140133\begin{description} 
     
    161154  \begin{equation} 
    162155    \label{eq:PE_w_bbc} 
    163     w = -{\rm {\bf U}}_h \cdot  \nabla _h \left( H \right) 
     156    w = - \vect U_h \cdot \nabla_h (H) 
    164157  \end{equation} 
    165158  In addition, the ocean exchanges momentum with the earth through frictional processes. 
     
    167160  It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. 
    168161  Its specification depends on the nature of the physical parameterisation used for 
    169   ${\rm {\bf D}}^{\rm {\bf U}}$ in \autoref{eq:PE_dyn}. 
     162  $\vect D^{\vect U}$ in \autoref{eq:PE_dyn}. 
    170163  It is discussed in \autoref{eq:PE_zdf}.% and Chap. III.6 to 9. 
    171164\item[Atmosphere - ocean interface:] 
     
    174167  \[ 
    175168    % \label{eq:PE_w_sbc} 
    176     w = \frac{\partial \eta }{\partial t} 
    177     + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot  \nabla _h \left( \eta \right) 
    178     + P-E 
     169    w = \pd[\eta]{t} + \lt. \vect U_h \rt|_{z = \eta} \cdot \nabla_h (\eta) + P - E 
    179170  \] 
    180171  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$. 
     172  leads to the continuity of pressure across the interface $z = \eta$. 
    182173  The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. 
    183174\item[Sea ice - ocean interface:] 
    184175  the ocean and sea ice exchange heat, salt, fresh water and momentum. 
    185176  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$). 
     177  Sea ice salinity is very low ($\sim4-6 \, psu$) compared to those of the ocean ($\sim34 \, psu$). 
    187178  The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected. 
    188179\end{description} 
    189  
    190  
    191 %\newpage 
    192180 
    193181% ================================================================ 
    194182% The Horizontal Pressure Gradient 
    195183% ================================================================ 
    196 \section{Horizontal pressure gradient } 
     184\section{Horizontal pressure gradient} 
    197185\label{sec:PE_hor_pg} 
    198186 
     
    204192 
    205193The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at 
    206 a reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: 
    207 $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. 
     194a reference geopotential surface ($z = 0$) and a hydrostatic pressure $p_h$ such that: 
     195$p(i,j,k,t) = p_s(i,j,t) + p_h(i,j,k,t)$. 
    208196The latter is computed by integrating (\autoref{eq:PE_hydrostatic}), 
    209197assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:PE_eos}). 
     
    211199\[ 
    212200  % \label{eq:PE_pressure} 
    213   p_h \left( {i,j,z,t} \right) 
    214   = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma } 
     201  p_h (i,j,z,t) = \int_{\varsigma = z}^{\varsigma = 0} g \; \rho (T,S,\varsigma) \; d \varsigma 
    215202\] 
    216203Two strategies can be considered for the surface pressure term: 
     
    224211The phase speed of such waves is high (some hundreds of metres per second) so that 
    225212the time step would have to be very short if they were present in the model. 
    226 The latter strategy filters out these waves since the rigid lid approximation implies $\eta=0$, 
    227 \ie the sea surface is the surface $z=0$. 
     213The latter strategy filters out these waves since the rigid lid approximation implies $\eta = 0$, 
     214\ie the sea surface is the surface $z = 0$. 
    228215This well known approximation increases the surface wave speed to infinity and 
    229216modifies certain other longwave dynamics (\eg barotropic Rossby or planetary waves). 
    230217The rigid-lid hypothesis is an obsolescent feature in modern OGCMs. 
    231 It has been available until the release 3.1 of  \NEMO, and it has been removed in release 3.2 and followings. 
     218It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings. 
    232219Only the free surface formulation is now described in the this document (see the next sub-section). 
    233220 
     
    244231\begin{equation} 
    245232  \label{eq:PE_ssh} 
    246   \frac{\partial \eta }{\partial t}=-D+P-E 
    247   \quad \text{where} \ 
    248   D=\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right] 
     233  \pd[\eta]{t} = - D + P - E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt] 
    249234\end{equation} 
    250235and using (\autoref{eq:PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. 
     
    256241 
    257242Two choices can be made regarding the implementation of the free surface in the model, 
    258 depending on the physical processes of interest.  
     243depending on the physical processes of interest. 
    259244 
    260245$\bullet$ If one is interested in EGWs, in particular the tides and their interaction with 
    261246the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 
    262247then a non linear free surface is the most appropriate. 
    263 This means that no approximation is made in (\autoref{eq:PE_ssh}) and that 
     248This means that no approximation is made in \autoref{eq:PE_ssh} and that 
    264249the variation of the ocean volume is fully taken into account. 
    265250Note that in order to study the fast time scales associated with EGWs it is necessary to 
     
    272257not altering the slow barotropic Rossby waves. 
    273258If further, an approximative conservation of heat and salt contents is sufficient for the problem solved, 
    274 then it is sufficient to solve a linearized version of (\autoref{eq:PE_ssh}), 
     259then it is sufficient to solve a linearized version of \autoref{eq:PE_ssh}, 
    275260which still allows to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}. 
    276261Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. 
     
    286271(see \autoref{subsec:DYN_spg_ts}). 
    287272 
    288 %\newpage 
    289  
    290273% ================================================================ 
    291274% Curvilinear z-coordinate System 
     
    293276\section{Curvilinear \textit{z-}coordinate system} 
    294277\label{sec:PE_zco} 
    295  
    296278 
    297279% ------------------------------------------------------------------------------------------------------------- 
     
    318300The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics. 
    319301 
    320 Let (\textit{i},\textit{j},\textit{k}) be a set of orthogonal curvilinear coordinates on 
     302Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on 
    321303the sphere associated with the positively oriented orthogonal set of unit vectors 
    322 (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that 
    323 \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, 
     304$(i,j,k)$ linked to the earth such that 
     305$k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 
    324306\ie along geopotential surfaces (\autoref{fig:referential}). 
    325307Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by 
    326308the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and 
    327 the distance from the centre of the earth $a+z(k)$ where $a$ is the earth's radius and 
     309the distance from the centre of the earth $a + z(k)$ where $a$ is the earth's radius and 
    328310$z$ the altitude above a reference sea level (\autoref{fig:referential}). 
    329311The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, 
     
    332314  \label{eq:scale_factors} 
    333315  \begin{aligned} 
    334     e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda}{\partial i}\cos \varphi } \right)^2 
    335         +\left( {\frac{\partial \varphi }{\partial i}} \right)^2} \right]^{1/2} \\ 
    336     e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda }{\partial j}\cos \varphi } \right)^2+ 
    337         \left( {\frac{\partial \varphi }{\partial j}} \right)^2} \right]^{1/2} \\ 
    338     e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\ 
     316    e_1 &= (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \\ 
     317    e_2 &= (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \\ 
     318    e_3 &= \lt( \pd[z]{k} \rt) 
    339319  \end{aligned} 
    340320\end{equation} 
     
    343323\begin{figure}[!tb] 
    344324  \begin{center} 
    345     \includegraphics[width=0.60\textwidth]{Fig_I_earth_referential} 
    346     \caption{  \protect\label{fig:referential} 
     325    \includegraphics[]{Fig_I_earth_referential} 
     326    \caption{ 
     327      \protect\label{fig:referential} 
    347328      the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
    348       coordinate system (\textbf{i},\textbf{j},\textbf{k}). } 
     329      coordinate system $(i,j,k)$. 
     330    } 
    349331  \end{center} 
    350332\end{figure} 
    351333%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    352334 
    353 Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by $a$ in 
     335Since the ocean depth is far smaller than the earth's radius, $a + z$, can be replaced by $a$ in 
    354336(\autoref{eq:scale_factors}) (thin-shell approximation). 
    355337The resulting horizontal scale factors $e_1$, $e_2$  are independent of $k$ while 
    356 the vertical scale factor is a single function of $k$ as \textbf{k} is parallel to \textbf{z}. 
     338the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$. 
    357339The scalar and vector operators that appear in the primitive equations 
    358340(\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can be written in the tensorial form, 
     
    360342\begin{subequations} 
    361343  % \label{eq:PE_discrete_operators} 
    362   \begin{equation} 
     344  \begin{gather} 
    363345    \label{eq:PE_grad} 
    364     \nabla q=\frac{1}{e_1 }\frac{\partial q}{\partial i}\;{\rm {\bf 
    365         i}}+\frac{1}{e_2 }\frac{\partial q}{\partial j}\;{\rm {\bf j}}+\frac{1}{e_3 
    366     }\frac{\partial q}{\partial k}\;{\rm {\bf k}}    \\ 
    367   \end{equation} 
    368   \begin{equation} 
     346    \nabla q =   \frac{1}{e_1} \pd[q]{i} \; \vect i 
     347               + \frac{1}{e_2} \pd[q]{j} \; \vect j 
     348               + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ 
    369349    \label{eq:PE_div} 
    370     \nabla \cdot {\rm {\bf A}} 
    371     = \frac{1}{e_1 \; e_2} \left[ 
    372       \frac{\partial \left(e_2 \; a_1\right)}{\partial i } 
    373       +\frac{\partial \left(e_1 \; a_2\right)}{\partial j }       \right] 
    374     + \frac{1}{e_3} \left[ \frac{\partial a_3}{\partial k }   \right] 
    375   \end{equation} 
    376   \begin{equation} 
     350    \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] 
     351                           + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] 
     352  \end{gather} 
     353  \begin{multline} 
    377354    \label{eq:PE_curl} 
    378     \begin{split} 
    379       \nabla \times \vect{A} = 
    380       \left[ {\frac{1}{e_2 }\frac{\partial a_3}{\partial j} 
    381           -\frac{1}{e_3 }\frac{\partial a_2 }{\partial k}} \right] \; \vect{i} 
    382       &+\left[ {\frac{1}{e_3 }\frac{\partial a_1 }{\partial k} 
    383           -\frac{1}{e_1 }\frac{\partial a_3 }{\partial i}} \right] \; \vect{j}      \\ 
    384       &+\frac{1}{e_1 e_2 } \left[ {\frac{\partial \left( {e_2 a_2 } \right)}{\partial i} 
    385           -\frac{\partial \left( {e_1 a_1 } \right)}{\partial j}} \right] \; \vect{k} 
    386     \end{split} 
    387   \end{equation} 
    388   \begin{equation} 
     355      \nabla \times \vect{A} =   \lt[ \frac{1}{e_2} \pd[a_3]{j} - \frac{1}{e_3} \pd[a_2]{k}   \rt] \vect i \\ 
     356                               + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i}   \rt] \vect j \\ 
     357                               + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i} - \pd[(e_1 a_1)]{j} \rt] \vect k 
     358  \end{multline} 
     359  \begin{gather} 
    389360    \label{eq:PE_lap} 
    390     \Delta q = \nabla \cdot \left(  \nabla q \right) 
    391   \end{equation} 
    392   \begin{equation} 
     361    \Delta q = \nabla \cdot (\nabla q) \\ 
    393362    \label{eq:PE_lap_vector} 
    394     \Delta {\rm {\bf A}} = 
    395     \nabla \left( \nabla \cdot {\rm {\bf A}} \right) 
    396     - \nabla \times \left(  \nabla \times {\rm {\bf A}} \right) 
    397   \end{equation} 
     363    \Delta \vect A = \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A) 
     364  \end{gather} 
    398365\end{subequations} 
    399 where $q$ is a scalar quantity and ${\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinate system. 
     366where $q$ is a scalar quantity and $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system. 
    400367 
    401368% ------------------------------------------------------------------------------------------------------------- 
     
    408375it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using 
    409376\autoref{eq:PE_grad}) to \autoref{eq:PE_lap_vector}. 
    410 Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate system and 
     377Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, the velocity in the $(i,j,k)$ coordinates system and 
    411378define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 
    412 \begin{equation} 
     379\begin{gather} 
    413380  \label{eq:PE_curl_Uh} 
    414   \zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v} 
    415         \right)}{\partial i}-\frac{\partial \left( {e_1 \,u} \right)}{\partial j}} 
    416   \right] 
    417 \end{equation} 
    418 \begin{equation} 
     381  \zeta = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, v)]{i} - \pd[(e_1 \, u)]{j} \rt] \\ 
    419382  \label{eq:PE_div_Uh} 
    420   \chi =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,u} 
    421         \right)}{\partial i}+\frac{\partial \left( {e_1 \,v} \right)}{\partial j}} 
    422   \right] 
    423 \end{equation} 
     383  \chi  = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, u)]{i} + \pd[(e_1 \, v)]{j} \rt] 
     384\end{gather} 
    424385 
    425386Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that 
    426387$e_3$  is a function of the single variable $k$, 
    427 the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows: 
    428 \begin{flalign*} 
    429   &\left[ {\left( {  \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}} 
    430       +\frac{1}{2}   \nabla \left( {{\rm {\bf U}}^2} \right)}   \right]_h        & 
    431 \end{flalign*} 
    432 \begin{flalign*} 
    433   &\qquad=\left( {{ 
    434         \begin{array}{*{20}c} 
    435           {\left[    {   \frac{1}{e_3} \frac{\partial u  }{\partial k} 
    436           -\frac{1}{e_1}   \frac{\partial w  }{\partial i} } \right] w - \zeta \; v }     \\ 
    437           {\zeta \; u - \left[ {  \frac{1}{e_2} \frac{\partial w}{\partial j} 
    438           -\frac{1}{e_3} \frac{\partial v}{\partial k} } \right] \ w}    \\ 
    439         \end{array} 
    440       }} \right) 
    441   +\frac{1}{2} \left( {{ 
    442         \begin{array}{*{20}c} 
    443           { \frac{1}{e_1}  \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial i}}  \hfill    \\ 
    444           { \frac{1}{e_2}  \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial j}}  \hfill    \\ 
    445         \end{array} 
    446       }} \right)        & 
    447 \end{flalign*} 
    448 \begin{flalign*} 
    449   & \qquad =\left( {{ 
    450         \begin{array}{*{20}c} 
    451           {-\zeta \; v} \hfill \\ 
    452           { \zeta \; u} \hfill \\ 
    453         \end{array} 
    454       }} \right) 
    455   +\frac{1}{2}\left( {{ 
    456         \begin{array}{*{20}c} 
    457           {\frac{1}{e_1 }\frac{\partial \left( {u^2+v^2} \right)}{\partial i}} \hfill  \\ 
    458           {\frac{1}{e_2 }\frac{\partial \left( {u^2+v^2} \right)}{\partial j}} \hfill  \\ 
    459         \end{array} 
    460       }} \right) 
    461   +\frac{1}{e_3 }\left( {{ 
    462         \begin{array}{*{20}c} 
    463           { w \; \frac{\partial u}{\partial k}}    \\ 
    464           { w \; \frac{\partial v}{\partial k}}    \\ 
    465         \end{array} 
    466       }} \right) 
    467   -\left( {{ 
    468         \begin{array}{*{20}c} 
    469           {\frac{w}{e_1}\frac{\partial w}{\partial i} -\frac{1}{2e_1}\frac{\partial w^2}{\partial i}} \hfill \\ 
    470           {\frac{w}{e_2}\frac{\partial w}{\partial j} -\frac{1}{2e_2}\frac{\partial w^2}{\partial j}} \hfill \\ 
    471         \end{array} 
    472       }} \right)        & 
    473 \end{flalign*} 
    474  
     388$NLT$ the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows: 
     389\begin{alignat*}{2} 
     390  &NLT &=   &\lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ 
     391  &    &=   &\lt( 
     392    \begin{array}{*{20}c} 
     393                \lt[ \frac{1}{e_3} \pd[u]{k} - \frac{1}{e_1} \pd[w]{i} \rt] w - \zeta \; v   \\ 
     394                \zeta \; u - \lt[ \frac{1}{e_2} \pd[w]{j} - \frac{1}{e_3} \pd[v]{k} \rt] \ w 
     395    \end{array} 
     396                                                                                             \rt) 
     397          + \frac{1}{2} \lt( 
     398    \begin{array}{*{20}c} 
     399                             \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\ 
     400                             \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j} 
     401    \end{array} 
     402                                                                     \rt) \\ 
     403  &    &=   &\lt( 
     404    \begin{array}{*{20}c} 
     405                  -\zeta \; v \\ 
     406                   \zeta \; u 
     407    \end{array} 
     408                              \rt) 
     409          + \frac{1}{2} \lt( 
     410    \begin{array}{*{20}c} 
     411                             \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\ 
     412                             \frac{1}{e_2} \pd[(u^2 + v^2)]{j} 
     413    \end{array} 
     414                                                               \rt) \\ 
     415  &    &  &+ \frac{1}{e_3} \lt( 
     416    \begin{array}{*{20}c} 
     417                                w \; \pd[u]{k} \\ 
     418                                w \; \pd[v]{k} 
     419    \end{array} 
     420                                               \rt) 
     421           - \lt( 
     422    \begin{array}{*{20}c} 
     423                  \frac{w}{e_1} \pd[w]{i} - \frac{1}{2 e_1} \pd[w^2]{i} \\ 
     424                  \frac{w}{e_2} \pd[w]{j} - \frac{1}{2 e_2} \pd[w^2]{j} 
     425    \end{array} 
     426                                                                        \rt) 
     427\end{alignat*} 
    475428The last term of the right hand side is obviously zero, and thus the nonlinear term of 
    476429\autoref{eq:PE_dyn} is written in the $(i,j,k)$ coordinate system: 
    477430\begin{equation} 
    478431  \label{eq:PE_vector_form} 
    479   \left[ {\left( {   \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}} 
    480       +\frac{1}{2}   \nabla \left( {{\rm {\bf U}}^2} \right)}   \right]_h 
    481   =\zeta 
    482   \;{\rm {\bf k}}\times {\rm {\bf U}}_h +\frac{1}{2}\nabla _h \left( {{\rm 
    483         {\bf U}}_h^2 } \right)+\frac{1}{e_3 }w\frac{\partial {\rm {\bf U}}_h 
    484   }{\partial k} 
     432  NLT =   \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) 
     433        + \frac{1}{e_3} w \pd[\vect U_h]{k} 
    485434\end{equation} 
    486435 
     
    489438\ie to write it as the divergence of fluxes. 
    490439For example, the first component of \autoref{eq:PE_vector_form} (the $i$-component) is transformed as follows: 
    491 \begin{flalign*} 
    492   &{ 
    493     \begin{array}{*{20}l} 
    494       \left[ {\left( {\nabla \times \vect{U}} \right)\times \vect{U} 
    495       +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i   % \\ 
    496   % \\ 
    497       = - \zeta \;v 
    498       + \frac{1}{2\;e_1 } \frac{\partial \left( {u^2+v^2} \right)}{\partial i} 
    499       + \frac{1}{e_3}w \ \frac{\partial u}{\partial k}         \\ \\ 
    500       \qquad =\frac{1}{e_1 \; e_2} \left(    -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i} 
    501       +v\frac{\partial \left( {e_1 \,u} \right)}{\partial j}    \right) 
    502       +\frac{1}{e_1 e_2 }\left(  +e_2 \; u\frac{\partial u}{\partial i} 
    503       +e_2 \; v\frac{\partial v}{\partial i}              \right) 
    504       +\frac{1}{e_3}       \left(   w\;\frac{\partial u}{\partial k}       \right)   \\ 
     440\begin{alignat*}{2} 
     441  &NLT_i &= &- \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\ 
     442  &      &=  &\frac{1}{e_1 \; e_2} \lt( -v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt) 
     443            + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) \\ 
     444  &      & &+ \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\ 
     445  &      &=  &\frac{1}{e_1 \; e_2} \lt[ - \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt) 
     446                                     + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j} -         e_1 \, u \pd[v]{j} \rt) \rt. \\ 
     447  &      &                       &\lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i} - u \pd[ \lt( e_2 u \rt)]{i} \rt) 
     448                                     + e_2 v \pd[v]{i}                                                         \rt] \\ 
     449  &      & &+ \frac{1}{e_3} \lt( \pd[(w \, u)]{k} - u \pd[w]{k} \rt) \\ 
     450  &      &=  &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt) 
     451            + \frac{1}{e_3} \pd[(w \, u)]{k} \\ 
     452  &      & &+ \frac{1}{e_1 e_2} \lt[ - u \lt( \pd[(e_1 v)]{j} - v \, \pd[e_1]{j} \rt) 
     453                                  - u \pd[(e_2 u)]{i}                              \rt] 
     454            - \frac{1}{e_3} \pd[w]{k} u \\ 
     455  &      & &+ \frac{1}{e_1 e_2} \lt( - v^2 \pd[e_2]{i} \rt) \\ 
     456  &      &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u 
     457            + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ 
     458  \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it comes:} 
     459  &      &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i} - u \; \pd[e_1]{j} \rt) (-v) 
     460\end{alignat*} 
     461 
     462The flux form of the momentum advection term is therefore given by: 
     463\begin{equation} 
     464  \label{eq:PE_flux_form} 
     465  NLT =   \nabla \cdot \lt( 
     466    \begin{array}{*{20}c} 
     467                            \vect U \, u \\ 
     468                            \vect U \, v 
    505469    \end{array} 
    506   }         & 
    507 \end{flalign*} 
    508 \begin{flalign*} 
    509   &{ 
    510     \begin{array}{*{20}l} 
    511       \qquad =\frac{1}{e_1 \; e_2}  \left\{ 
    512       -\left(         v^2  \frac{\partial e_2                                }{\partial i} 
    513       +e_2 \,v    \frac{\partial v                                   }{\partial i}     \right) 
    514       +\left(           \frac{\partial \left( {e_1 \,u\,v}  \right)}{\partial j} 
    515       -e_1 \,u    \frac{\partial v                                   }{\partial j}  \right)  \right. \\ 
    516       \left.   \qquad \qquad \quad 
    517       +\left(           \frac{\partial \left( {e_2 u\,u}     \right)}{\partial i} 
    518       -u       \frac{\partial \left( {e_2 u}         \right)}{\partial i}  \right) 
    519       +e_2 v            \frac{\partial v                                    }{\partial i} 
    520       \right\} 
    521       +\frac{1}{e_3} \left( 
    522       \frac{\partial \left( {w\,u} \right)         }{\partial k} 
    523       -u       \frac{\partial w                    }{\partial k}  \right) \\ 
    524     \end{array} 
    525   }      & 
    526 \end{flalign*} 
    527 \begin{flalign*} 
    528   & 
    529   { 
    530     \begin{array}{*{20}l} 
    531       \qquad =\frac{1}{e_1 \; e_2}  \left( 
    532       \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 
    533       +        \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j}  \right) 
    534       +\frac{1}{e_3 }      \frac{\partial \left( {w\,u       } \right)}{\partial k} \\ 
    535       \qquad \qquad \quad 
    536       +\frac{1}{e_1 e_2 }     \left( 
    537       -u \left(   \frac{\partial \left( {e_1 v   } \right)}{\partial j} 
    538       -v\,\frac{\partial e_1 }{\partial j}             \right) 
    539       -u       \frac{\partial \left( {e_2 u   } \right)}{\partial i} 
    540       \right) 
    541       -\frac{1}{e_3 }      \frac{\partial w}{\partial k} u 
    542       +\frac{1}{e_1 e_2 }\left(  -v^2\frac{\partial e_2   }{\partial i}     \right) 
    543     \end{array} 
    544   }      & 
    545 \end{flalign*} 
    546 \begin{flalign*} 
    547   &{ 
    548     \begin{array}{*{20}l} 
    549       \qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right) 
    550       -   \left( \nabla \cdot {\rm {\bf U}} \right) \ u 
    551       +\frac{1}{e_1 e_2 }\left( 
    552       -v^2     \frac{\partial e_2 }{\partial i} 
    553       +uv   \,    \frac{\partial e_1 }{\partial j}    \right) \\ 
    554     \end{array} 
    555   }      & 
    556 \end{flalign*} 
    557 as $\nabla \cdot {\rm {\bf U}}\;=0$ (incompressibility) it comes: 
    558 \begin{flalign*} 
    559   &{ 
    560     \begin{array}{*{20}l} 
    561       \qquad = \nabla \cdot \left(  {{\rm {\bf U}}\,u}      \right) 
    562       +  \frac{1}{e_1 e_2 }   \left( v \; \frac{\partial e_2}{\partial i} 
    563       -u \; \frac{\partial e_1}{\partial j}  \right)  \left( -v \right) 
    564     \end{array} 
    565   }      & 
    566 \end{flalign*} 
    567  
    568 The flux form of the momentum advection term is therefore given by: 
    569 \begin{multline} 
    570   \label{eq:PE_flux_form} 
    571   \left[ 
    572     \left(  {\nabla \times {\rm {\bf U}}}    \right) \times {\rm {\bf U}} 
    573     +\frac{1}{2}  \nabla \left(  {{\rm {\bf U}}^2}    \right) 
    574   \right]_h \\ 
    575   = \nabla \cdot  \left( {{ 
    576         \begin{array}{*{20}c} 
    577           {\rm {\bf U}} \, u  \hfill \\ 
    578           {\rm {\bf U}} \, v  \hfill \\ 
    579         \end{array} 
    580       }} 
    581   \right) 
    582   +\frac{1}{e_1 e_2 }      \left( 
    583     v\frac{\partial e_2}{\partial i} 
    584     -u\frac{\partial e_1}{\partial j} 
    585   \right) {\rm {\bf k}} \times {\rm {\bf U}}_h 
    586 \end{multline} 
     470                                         \rt) 
     471        + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \vect k \times \vect U_h 
     472\end{equation} 
    587473 
    588474The flux form has two terms, 
    589475the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) 
    590476and the second one is due to the curvilinear nature of the coordinate system used. 
    591 The latter is called the \emph{metric} term and can be viewed as a modification of the Coriolis parameter:  
     477The latter is called the \textit{metric} term and can be viewed as a modification of the Coriolis parameter:  
    592478\[ 
    593479  % \label{eq:PE_cor+metric} 
    594   f \to f + \frac{1}{e_1\;e_2}  \left(  v \frac{\partial e_2}{\partial i} 
    595     -u \frac{\partial e_1}{\partial j}  \right) 
     480  f \to f + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) 
    596481\] 
    597482 
    598483Note that in the case of geographical coordinate, 
    599 \ie when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, 
    600 we recover the commonly used modification of the Coriolis parameter $f \to f+(u/a) \tan \varphi$. 
     484\ie when $(i,j) \to (\lambda,\varphi)$ and $(e_1,e_2) \to (a \, \cos \varphi,a)$, 
     485we recover the commonly used modification of the Coriolis parameter $f \to f + (u / a) \tan \varphi$. 
    601486 
    602487To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be written in 
    603488the following tensorial formalism: 
    604489 
    605 \vspace{+10pt} 
    606 $\bullet$ \textbf{Vector invariant form of the momentum equations} : 
    607  
    608 \begin{subequations} 
    609   \label{eq:PE_dyn_vect} 
    610   \[ 
     490\begin{itemize} 
     491\item 
     492  \textbf{Vector invariant form of the momentum equations}: 
     493  \begin{equation} 
     494    \label{eq:PE_dyn_vect} 
     495    \begin{split} 
    611496    % \label{eq:PE_dyn_vect_u} 
    612     \begin{split} 
    613       \frac{\partial u}{\partial t} 
    614       = +   \left( {\zeta +f} \right)\,v 
    615       -   \frac{1}{2\,e_1}           \frac{\partial}{\partial i} \left(  u^2+v^2   \right) 
    616       -   \frac{1}{e_3    }  w     \frac{\partial u}{\partial k}      &      \\ 
    617       -   \frac{1}{e_1    }            \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho_o}    \right) 
    618       &+   D_u^{\vect{U}}  +   F_u^{\vect{U}}      \\ \\ 
    619       \frac{\partial v}{\partial t} = 
    620       -   \left( {\zeta +f} \right)\,u 
    621       -   \frac{1}{2\,e_2 }        \frac{\partial }{\partial j}\left(  u^2+v^2  \right) 
    622       -   \frac{1}{e_3     }   w  \frac{\partial v}{\partial k}     &      \\ 
    623       -   \frac{1}{e_2     }        \frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}  \right) 
    624       &+  D_v^{\vect{U}}  +   F_v^{\vect{U}} 
     497      \pd[u]{t} = &+ (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) 
     498                   - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 
     499                  &+ D_u^{\vect U} + F_u^{\vect U} \\ 
     500      \pd[v]{t} = &- (\zeta + f) \, u - \frac{1}{2 e_2} \pd[]{j} (u^2 + v^2) 
     501                   - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\  
     502                  &+ D_v^{\vect U} + F_v^{\vect U} 
    625503    \end{split} 
    626   \] 
    627 \end{subequations} 
    628  
    629  
    630 \vspace{+10pt} 
    631 $\bullet$ \textbf{flux form of the momentum equations} : 
    632 \begin{subequations} 
     504  \end{equation} 
     505\item 
     506  \textbf{flux form of the momentum equations}: 
    633507  % \label{eq:PE_dyn_flux} 
    634508  \begin{multline*} 
    635509    % \label{eq:PE_dyn_flux_u} 
    636     \frac{\partial u}{\partial t}= 
    637     +   \left( { f + \frac{1}{e_1 \; e_2} 
    638         \left(     v \frac{\partial e_2}{\partial i} 
    639           -u \frac{\partial e_1}{\partial j}    \right)}    \right) \, v    \\ 
    640     - \frac{1}{e_1 \; e_2}    \left( 
    641       \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 
    642       +        \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j}  \right) 
    643     - \frac{1}{e_3 }\frac{\partial \left( {         w\,u} \right)}{\partial k}    \\ 
    644     -   \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho_o}   \right) 
    645     +   D_u^{\vect{U}} +   F_u^{\vect{U}} 
     510    \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 \\ 
     511                - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\ 
     512                - \frac{1}{e_3} \pd[(w \, u)]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
     513                + D_u^{\vect U} + F_u^{\vect U} 
    646514  \end{multline*} 
    647515  \begin{multline*} 
    648516    % \label{eq:PE_dyn_flux_v} 
    649     \frac{\partial v}{\partial t}= 
    650     -   \left( { f + \frac{1}{e_1 \; e_2} 
    651         \left(     v \frac{\partial e_2}{\partial i} 
    652           -u \frac{\partial e_1}{\partial j}    \right)}    \right) \, u   \\ 
    653     \frac{1}{e_1 \; e_2}   \left( 
    654       \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i} 
    655       +        \frac{\partial \left( {e_1 \,v\,v} \right)}{\partial j}  \right) 
    656     - \frac{1}{e_3 } \frac{\partial \left( {        w\,v} \right)}{\partial k}    \\ 
    657     -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}    \right) 
    658     +  D_v^{\vect{U}} +  F_v^{\vect{U}} 
     517    \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 \\ 
     518                + \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\ 
     519                - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
     520                + D_v^{\vect U} + F_v^{\vect U} 
    659521  \end{multline*} 
    660 \end{subequations} 
    661 where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and 
    662 $p_s $, the surface pressure, is given by: 
    663 \[ 
     522  where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and $p_s$, the surface pressure, 
     523  is given by: 
     524  \[ 
    664525  % \label{eq:PE_spg} 
    665   p_s =  \rho \,g \,\eta 
    666 \] 
    667 with $\eta$ is solution of \autoref{eq:PE_ssh}. 
    668  
    669 The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: 
    670 \[ 
     526    p_s = \rho \,g \, \eta 
     527  \] 
     528  with $\eta$ is solution of \autoref{eq:PE_ssh}. 
     529 
     530  The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: 
     531  \[ 
    671532  % \label{eq:w_diag} 
    672   \frac{\partial w}{\partial k}=-\chi \;e_3 
    673 \] 
    674 \[ 
     533    \pd[w]{k} = - \chi \; e_3 \qquad  
    675534  % \label{eq:hp_diag} 
    676   \frac{\partial p_h }{\partial k}=-\rho \;g\;e_3 
    677 \] 
    678 where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}. 
    679  
    680 \vspace{+10pt} 
    681 $\bullet$ \textit{tracer equations} : 
    682 \[ 
    683   % \label{eq:S} 
    684   \frac{\partial T}{\partial t} = 
    685   -\frac{1}{e_1 e_2 }\left[ {    \frac{\partial \left( {e_2 T\,u} \right)}{\partial i} 
    686       +\frac{\partial \left( {e_1 T\,v} \right)}{\partial j}} \right] 
    687   -\frac{1}{e_3 }\frac{\partial \left( {T\,w} \right)}{\partial k} + D^T + F^T 
    688 \] 
    689 \[ 
    690   % \label{eq:T} 
    691   \frac{\partial S}{\partial t} = 
    692   -\frac{1}{e_1 e_2 }\left[     {\frac{\partial \left( {e_2 S\,u} \right)}{\partial i} 
    693       +\frac{\partial \left( {e_1 S\,v} \right)}{\partial j}} \right] 
    694   -\frac{1}{e_3 }\frac{\partial \left( {S\,w} \right)}{\partial k} + D^S + F^S 
    695 \] 
    696 \[ 
    697   % \label{eq:rho} 
    698   \rho =\rho \left( {T,S,z(k)} \right) 
    699 \] 
    700  
    701 The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 
     535    \pd[p_h]{k} = - \rho \; g \; e_3 
     536  \] 
     537  where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}. 
     538\item \textit{tracer equations}: 
     539  \[ 
     540    %\label{eq:S} 
     541    \pd[T]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] 
     542                - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 
     543    %\label{eq:T} 
     544    \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] 
     545                - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S 
     546    %\label{eq:rho} 
     547    \rho = \rho \big( T,S,z(k) \big) 
     548  \] 
     549\end{itemize} 
     550 
     551The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 
    702552It will be defined in \autoref{eq:PE_zdf}. 
    703 The nature and formulation of ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, 
     553The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms, 
    704554are discussed in \autoref{chap:SBC}. 
    705555 
    706  
    707 \newpage  
     556\newpage 
    708557 
    709558% ================================================================ 
     
    717566Second the ocean floor depends on the geographical position, 
    718567varying from more than 6,000 meters in abyssal trenches to zero at the coast. 
    719 Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing.  
     568Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 
    720569Therefore, in order to represent the ocean with respect to 
    721570the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height 
     
    736585\begin{equation} 
    737586  \label{eq:PE_s} 
    738   s=s(i,j,k,t) 
     587  s = s(i,j,k,t) 
    739588\end{equation} 
    740589with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, 
     
    759608 
    760609%\gmcomment{ 
    761  
    762610%A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient... 
    763  
    764611the generalized vertical coordinates used in ocean modelling are not orthogonal, 
    765612which contrasts with many other applications in mathematical physics. 
     
    772619The key motivation for maintaining the same horizontal velocity component is that 
    773620the hydrostatic and geostrophic balances are dominant in the large-scale ocean. 
    774 Use of an alternative quasi-horizontal velocity, for example one oriented parallel to the generalized surface, 
     621Use of an alternative quasi -horizontal velocity, for example one oriented parallel to the generalized surface, 
    775622would lead to unacceptable numerical errors. 
    776 Correspondingly, the vertical direction is anti-parallel to the gravitational force in 
     623Correspondingly, the vertical direction is anti -parallel to the gravitational force in 
    777624all of the coordinate systems. 
    778 We do not choose the alternative of a quasi-vertical direction oriented normal to 
    779 the surface of a constant generalized vertical coordinate.  
     625We do not choose the alternative of a quasi -vertical direction oriented normal to 
     626the surface of a constant generalized vertical coordinate. 
    780627 
    781628It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between 
     
    786633volume or mass conservation. 
    787634In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about 
    788 the physical processes producing a flux across the layer interfaces.  
    789  
     635the physical processes producing a flux across the layer interfaces. 
    790636 
    791637In this section we first establish the PE in the generalised vertical $s$-coordinate, 
    792 then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde.   
     638then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde. 
    793639%} 
    794640 
     
    796642% The s-coordinate Formulation 
    797643% ------------------------------------------------------------------------------------------------------------- 
    798 \subsection{\textit{S-}coordinate formulation} 
    799  
    800 Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k=z$ and thus $e_3=1$, 
    801 we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, 
    802 which includes $z$-, \zstar- and $\sigma-$coordinates as special cases 
    803 ($s=z$, $s=\zstar$, and $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). 
     644\subsection{\textit{S}-coordinate formulation} 
     645 
     646Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k = z$ and 
     647thus $e_3 = 1$, we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$, 
     648which includes $z$-, \zstar- and $\sigma$-coordinates as special cases 
     649($s = z$, $s = \zstar$, and $s = \sigma = z / H$ or $ = z / \lt( H + \eta \rt)$, resp.). 
    804650A formal derivation of the transformed equations is given in \autoref{apdx:A}. 
    805 Let us define the vertical scale factor by $e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), 
    806 and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by: 
     651Let us define the vertical scale factor by $e_3 = \partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), 
     652and the slopes in the $(i,j)$ directions between $s$- and $z$-surfaces by: 
    807653\begin{equation} 
    808654  \label{eq:PE_sco_slope} 
    809   \sigma_1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s 
    810   \quad \text{, and } \quad 
    811   \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s 
     655  \sigma_1 = \frac{1}{e_1} \; \lt. \pd[z]{i} \rt|_s \quad \text{and} \quad 
     656  \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt|_s 
    812657\end{equation} 
    813 We also introduce  $\omega $, a dia-surface velocity component, defined as the velocity  
     658We also introduce $\omega$, a dia-surface velocity component, defined as the velocity  
    814659relative to the moving $s$-surfaces and normal to them: 
    815660\[ 
    816661  % \label{eq:PE_sco_w} 
    817   \omega  = w - e_3 \, \frac{\partial s}{\partial t} - \sigma_1 \,u - \sigma_2 \,v    \\ 
     662  \omega = w - e_3 \, \pd[s]{t} - \sigma_1 \, u - \sigma_2 \, v 
    818663\] 
    819664 
    820 The equations solved by the ocean model \autoref{eq:PE} in $s-$coordinate can be written as follows 
     665The equations solved by the ocean model \autoref{eq:PE} in $s$-coordinate can be written as follows 
    821666(see \autoref{sec:A_momentum}): 
    822667 
    823  \vspace{0.5cm} 
    824 $\bullet$ Vector invariant form of the momentum equation : 
    825 \begin{multline*} 
     668\begin{itemize} 
     669\item \textbf{Vector invariant form of the momentum equation}: 
     670  \begin{multline*} 
    826671  % \label{eq:PE_sco_u_vector} 
    827   \frac{\partial  u   }{\partial t}= 
    828   +   \left( {\zeta +f} \right)\,v 
    829   -   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right) 
    830   -   \frac{1}{e_3} \omega \frac{\partial u}{\partial k}       \\ 
    831   -   \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho_o}    \right) 
    832   +  g\frac{\rho }{\rho_o}\sigma_1 
    833   +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad 
    834 \end{multline*} 
    835 \begin{multline*} 
     672    \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} \\ 
     673                - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_1 
     674                + D_u^{\vect U} + F_u^{\vect U} 
     675  \end{multline*} 
     676  \begin{multline*} 
    836677  % \label{eq:PE_sco_v_vector} 
    837   \frac{\partial v }{\partial t}= 
    838   -   \left( {\zeta +f} \right)\,u 
    839   -   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^2  \right) 
    840   -   \frac{1}{e_3 } \omega \frac{\partial v}{\partial k}         \\ 
    841   -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}  \right) 
    842   +  g\frac{\rho }{\rho_o }\sigma_2 
    843   +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad 
    844 \end{multline*} 
    845  
    846  \vspace{0.5cm} 
    847 $\bullet$ Flux form of the momentum equation : 
    848 \begin{multline*} 
     678    \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} \\ 
     679                - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_2 
     680                + D_v^{\vect U} + F_v^{\vect U} 
     681  \end{multline*} 
     682\item \textbf{Flux form of the momentum equation}: 
     683  \begin{multline*} 
    849684  % \label{eq:PE_sco_u_flux} 
    850   \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}= 
    851   +   \left( { f + \frac{1}{e_1 \; e_2 } 
    852       \left(    v \frac{\partial e_2}{\partial i} 
    853         -u \frac{\partial e_1}{\partial j}   \right)}    \right) \, v    \\ 
    854   - \frac{1}{e_1 \; e_2 \; e_3 }    \left( 
    855     \frac{\partial \left( {e_2 \, e_3 \, u\,u} \right)}{\partial i} 
    856     +       \frac{\partial \left( {e_1 \, e_3 \, v\,u} \right)}{\partial j}   \right) 
    857   - \frac{1}{e_3 }\frac{\partial \left( { \omega\,u} \right)}{\partial k}    \\ 
    858   - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho_o}    \right) 
    859   +  g\frac{\rho }{\rho_o}\sigma_1 
    860   +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad 
    861 \end{multline*} 
    862 \begin{multline*} 
     685    \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 \\ 
     686                                       - \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) \\ 
     687                                       - \frac{1}{e_3} \pd[(\omega \, u)]{k} 
     688                                       - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
     689                                       + g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 
     690  \end{multline*} 
     691  \begin{multline*} 
    863692  % \label{eq:PE_sco_v_flux} 
    864   \frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
    865   -   \left( { f + \frac{1}{e_1 \; e_2} 
    866       \left(    v \frac{\partial e_2}{\partial i} 
    867         -u \frac{\partial e_1}{\partial j}   \right)}    \right) \, u   \\ 
    868   - \frac{1}{e_1 \; e_2 \; e_3 }    \left( 
    869     \frac{\partial \left( {e_2 \; e_3  \,u\,v} \right)}{\partial i} 
    870     +       \frac{\partial \left( {e_1 \; e_3  \,v\,v} \right)}{\partial j}   \right) 
    871   - \frac{1}{e_3 } \frac{\partial \left( { \omega\,v} \right)}{\partial k}    \\ 
    872   -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o}  \right) 
    873   +  g\frac{\rho }{\rho_o }\sigma_2 
    874   +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad 
    875 \end{multline*} 
    876  
    877 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, 
    878 and the hydrostatic pressure have the same expressions as in $z$-coordinates although 
    879 they do not represent exactly the same quantities. 
    880 $\omega$ is provided by the continuity equation (see \autoref{apdx:A}): 
    881 \[ 
     693    \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 \\ 
     694                                       - \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) \\ 
     695                                       - \frac{1}{e_3} \pd[(\omega \, v)]{k} 
     696                                       - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
     697                                       + g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 
     698  \end{multline*} 
     699  where the relative vorticity, $\zeta$, the surface pressure gradient, 
     700  and the hydrostatic pressure have the same expressions as in $z$-coordinates although 
     701  they do not represent exactly the same quantities. 
     702  $\omega$ is provided by the continuity equation (see \autoref{apdx:A}): 
     703  \[ 
    882704  % \label{eq:PE_sco_continuity} 
    883   \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0 
    884   \qquad \text{with }\;\; 
    885   \chi =\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3 \,u} 
    886         \right)}{\partial i}+\frac{\partial \left( {e_1 e_3 \,v} \right)}{\partial 
    887         j}} \right] 
    888 \] 
    889  
    890  \vspace{0.5cm} 
    891 $\bullet$ tracer equations: 
    892 \begin{multline*} 
     705    \pd[e_3]{t} + e_3 \; \chi + \pd[\omega]{s} = 0 \quad \text{with} \quad 
     706    \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) 
     707  \] 
     708\item \textit{tracer equations}: 
     709  \begin{multline*} 
    893710  % \label{eq:PE_sco_t} 
    894   \frac{1}{e_3} \frac{\partial \left(  e_3\,T  \right) }{\partial t}= 
    895   -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,T} \right)}{\partial i} 
    896       +\frac{\partial \left( {e_1 e_3\,v\,T} \right)}{\partial j}} \right]   \\ 
    897   -\frac{1}{e_3 }\frac{\partial \left( {T\,\omega } \right)}{\partial k}   + D^T + F^S   \qquad 
    898 \end{multline*} 
    899  
    900 \begin{multline*} 
     711    \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} 
     712                                                                    + \pd[(e_1 e_3 \, v \, T)]{j} \rt) \\ 
     713                                       - \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S 
     714  \end{multline*} 
     715  \begin{multline} 
    901716  % \label{eq:PE_sco_s} 
    902   \frac{1}{e_3} \frac{\partial \left(  e_3\,S  \right) }{\partial t}= 
    903   -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,S} \right)}{\partial i} 
    904       +\frac{\partial \left( {e_1 e_3\,v\,S} \right)}{\partial j}} \right]    \\ 
    905   -\frac{1}{e_3 }\frac{\partial \left( {S\,\omega } \right)}{\partial k}     + D^S + F^S   \qquad 
    906 \end{multline*} 
    907  
     717    \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} 
     718                                                                    + \pd[(e_1 e_3 \, v \, S)]{j} \rt) \\ 
     719                                       - \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S 
     720  \end{multline} 
     721\end{itemize} 
    908722The equation of state has the same expression as in $z$-coordinate, 
    909723and similar expressions are used for mixing and forcing terms. 
    910724 
    911725\gmcomment{ 
    912 \colorbox{yellow}{ to be updated $= = >$} 
    913 Add a few works on z and zps and s and underlies the differences between all of them 
    914 \colorbox{yellow}{ $< = =$ end update}  } 
    915  
    916  
     726  \colorbox{yellow}{ to be updated $= = >$} 
     727  Add a few works on z and zps and s and underlies the differences between all of them 
     728  \colorbox{yellow}{$< = =$ end update} 
     729} 
    917730 
    918731% ------------------------------------------------------------------------------------------------------------- 
    919732% Curvilinear \zstar-coordinate System 
    920733% ------------------------------------------------------------------------------------------------------------- 
    921 \subsection{Curvilinear \zstar--coordinate system} 
     734\subsection{Curvilinear \zstar-coordinate system} 
    922735\label{subsec:PE_zco_star} 
    923736 
     
    925738\begin{figure}[!b] 
    926739  \begin{center} 
    927     \includegraphics[width=1.0\textwidth]{Fig_z_zstar} 
    928     \caption{  \protect\label{fig:z_zstar} 
     740    \includegraphics[]{Fig_z_zstar} 
     741    \caption{ 
     742      \protect\label{fig:z_zstar} 
    929743      (a) $z$-coordinate in linear free-surface case ; 
    930       (b) $z-$coordinate in non-linear free surface case ; 
    931       (c) re-scaled height coordinate (become popular as the \zstar-coordinate 
    932       \citep{Adcroft_Campin_OM04} ). 
     744      (b) $z$-coordinate in non-linear free surface case ; 
     745      (c) re-scaled height coordinate 
     746      (become popular as the \zstar-coordinate \citep{Adcroft_Campin_OM04}). 
    933747    } 
    934748  \end{center} 
     
    936750%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    937751 
    938  
    939752In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. 
    940 These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site.  
    941  
    942 %\gmcomment{ 
     753These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site. 
     754 
    943755The \zstar coordinate approach is an unapproximated, non-linear free surface implementation which allows one to 
    944756deal with large amplitude free-surface variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. 
     
    947759as in the $z$-coordinate formulation, but is equally distributed over the full water column. 
    948760Thus vertical levels naturally follow sea-surface variations, with a linear attenuation with depth, 
    949 as illustrated by figure fig.1c. 
    950 Note that with a flat bottom, such as in fig.1c, the bottom-following $z$ coordinate and \zstar are equivalent. 
    951 The definition and modified oceanic equations for the rescaled vertical coordinate  \zstar, 
     761as illustrated by \autoref{fig:z_zstar}. 
     762Note that with a flat bottom, such as in \autoref{fig:z_zstar}, the bottom-following $z$ coordinate and \zstar are equivalent. 
     763The definition and modified oceanic equations for the rescaled vertical coordinate \zstar, 
    952764including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004). 
    953765The major points are summarized here. 
    954 The position ( \zstar) and vertical discretization (\zstar) are expressed as: 
     766The position (\zstar) and vertical discretization (\zstar) are expressed as: 
    955767\[ 
    956768  % \label{eq:z-star} 
    957   H +  \zstar = (H + z) / r \quad \text{and} \ \delta \zstar = \delta z / r \quad \text{with} \ r = \frac{H+\eta} {H} 
    958 \]  
     769  H + \zstar = (H + z)  / r \quad \text{and}  \quad \delta \zstar 
     770              = \delta z / r \quad \text{with} \quad r 
     771              = \frac{H + \eta}{H} 
     772\] 
    959773Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar, 
    960774the upper and lower boundaries are at fixed  \zstar position, 
    961 $\zstar = 0$ and  $\zstar = -H$ respectively. 
     775$\zstar = 0$ and $\zstar = -H$ respectively. 
    962776Also the divergence of the flow field is no longer zero as shown by the continuity equation: 
    963 \[  
    964   \frac{\partial r}{\partial t} = \nabla_{\zstar} \cdot \left( r \; \rm{\bf U}_h \right) 
    965   \left( r \; w\textit{*} \right) = 0  
    966 \]  
    967 %} 
    968  
     777\[ 
     778  \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) (r \; w *) = 0 
     779\] 
    969780 
    970781% from MOM4p1 documentation 
    971  
    972782To overcome problems with vanishing surface and/or bottom cells, we consider the zstar coordinate  
    973783\[ 
    974784  % \label{eq:PE_} 
    975   z^\star = H \left( \frac{z-\eta}{H+\eta} \right) 
     785  \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) 
    976786\] 
    977787 
     
    981791and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling. 
    982792 
    983 The surfaces of constant $z^\star$ are quasi-horizontal. 
    984 Indeed, the $z^\star$ coordinate reduces to $z$ when $\eta$ is zero. 
     793The surfaces of constant \zstar are quasi -horizontal. 
     794Indeed, the \zstar coordinate reduces to $z$ when $\eta$ is zero. 
    985795In general, when noting the large differences between 
    986796undulations of the bottom topography versus undulations in the surface height, 
    987 it is clear that surfaces constant $z^\star$ are very similar to the depth surfaces. 
     797it is clear that surfaces constant \zstar are very similar to the depth surfaces. 
    988798These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to 
    989799terrain following sigma models discussed in \autoref{subsec:PE_sco}. 
    990 Additionally, since $z^\star$ when $\eta = 0$, 
     800Additionally, since \zstar when $\eta = 0$, 
    991801no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography. 
    992802This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of 
    993803nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, 
    994804depending on the sophistication of the pressure gradient solver. 
    995 The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of 
    996 neutral physics parameterizations in $z^\star$ models using the same techniques as in $z$-models 
     805The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of 
     806neutral physics parameterizations in \zstar models using the same techniques as in $z$-models 
    997807(see Chapters 13-16 of \cite{Griffies_Bk04}) for a discussion of neutral physics in $z$-models, 
    998808as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). 
    999809 
    1000 The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$. 
     810The range over which \zstar varies is time independent $-H \leq \zstar \leq 0$. 
    1001811Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$. 
    1002 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$.  
    1003  
    1004 Because $z^\star$ has a time independent range, all grid cells have static increments ds, 
    1005 and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H.  
    1006 The $z^\star$ coordinate is therefore invisible to undulations of the free surface, 
     812This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 
     813 
     814Because \zstar has a time independent range, all grid cells have static increments ds, 
     815and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H. 
     816The \zstar coordinate is therefore invisible to undulations of the free surface, 
    1007817since it moves along with the free surface. 
    1008 This proper ty means that no spurious vertical transport is induced across surfaces of constant $z^\star$ by 
     818This proper ty means that no spurious vertical transport is induced across surfaces of constant \zstar by 
    1009819the motion of external gravity waves. 
    1010820Such spurious transpor t can be a problem in z-models, especially those with tidal forcing. 
    1011 Quite generally, the time independent range for the $z^\star$ coordinate is a very convenient property that 
     821Quite generally, the time independent range for the \zstar coordinate is a very convenient property that 
    1012822allows for a nearly arbitrary ver tical resolution even in the presence of large amplitude fluctuations of 
    1013823the surface height, again so long as $\eta > -H$. 
    1014  
    1015824%end MOM doc %%% 
    1016825 
    1017  
    1018  
    1019 \newpage 
     826\newpage  
    1020827 
    1021828% ------------------------------------------------------------------------------------------------------------- 
     
    1035842For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. 
    1036843Topographic Rossby waves can be excited and can interact with the mean current. 
    1037 In the $z-$coordinate system presented in the previous section (\autoref{sec:PE_zco}), 
    1038 $z-$surfaces are geopotential surfaces. 
     844In the $z$-coordinate system presented in the previous section (\autoref{sec:PE_zco}), 
     845$z$-surfaces are geopotential surfaces. 
    1039846The bottom topography is discretised by steps. 
    1040847This often leads to a misrepresentation of a gradually sloping bottom and to 
     
    1043850One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of 
    1044851a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 
    1045 Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate). 
     852Another solution is to introduce a terrain-following coordinate system (hereafter $s$-coordinate). 
    1046853 
    1047854The $s$-coordinate avoids the discretisation error in the depth field since the layers of 
     
    1050857which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths, 
    1051858can easily be represented (with relatively low vertical resolution). 
    1052 A terrain-following model (hereafter $s-$model) also facilitates the modelling of the boundary layer flows over 
     859A terrain-following model (hereafter $s$-model) also facilitates the modelling of the boundary layer flows over 
    1053860a large depth range, which in the framework of the $z$-model would require high vertical resolution over 
    1054861the whole depth range. 
     
    1063870\begin{equation} 
    1064871  \label{eq:PE_p_sco} 
    1065   \left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial 
    1066     p}{\partial s}\left. {\nabla z} \right|_s  
     872  \nabla p |_z = \nabla p |_s - \pd[p]{s} \nabla z |_s 
    1067873\end{equation} 
    1068874 
    1069875The second term in \autoref{eq:PE_p_sco} depends on the tilt of the coordinate surface and 
    1070876introduces a truncation error that is not present in a $z$-model. 
    1071 In the special case of a $\sigma-$coordinate (\ie depth-normalised coordinate system $\sigma = z/H$), 
     877In the special case of a $\sigma$-coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 
    1072878\citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error. 
    1073879It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, 
     
    1085891(see \autoref{sec:DOM_zgr}. 
    1086892 
    1087 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of any finite difference model. 
     893For numerical reasons a minimum of diffusion is required along the coordinate surfaces of 
     894any finite difference model. 
    1088895It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. 
    1089896This is the case for a $z$-model as well as for a $s$-model. 
    1090 However, density varies more strongly on $s-$surfaces than on horizontal surfaces in regions of 
     897However, density varies more strongly on $s$-surfaces than on horizontal surfaces in regions of 
    1091898large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model. 
    1092899Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus 
     
    1101908(see \autoref{subsec:PE_ldf}). 
    1102909Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 
    1103 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}).  
    1104  
    1105 The $s-$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from 
     910strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 
     911 
     912The $s$-coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from 
    1106913similar models: 
    1107914it allows a representation of bottom topography with mixed full or partial step-like/terrain following topography; 
    1108915It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. 
    1109916 
    1110  
    1111 \newpage 
    1112  
    1113917% ------------------------------------------------------------------------------------------------------------- 
    1114918% Curvilinear z-tilde coordinate System 
    1115919% ------------------------------------------------------------------------------------------------------------- 
    1116 \subsection{\texorpdfstring{Curvilinear \ztilde--coordinate}{}} 
     920\subsection{\texorpdfstring{Curvilinear \ztilde-coordinate}{}} 
    1117921\label{subsec:PE_zco_tilde} 
    1118922 
    1119 The \ztilde-coordinate has been developed by \citet{Leclair_Madec_OM11}. 
     923The \ztilde -coordinate has been developed by \citet{Leclair_Madec_OM11}. 
    1120924It is available in \NEMO since the version 3.4. 
    1121925Nevertheless, it is currently not robust enough to be used in all possible configurations. 
    1122926Its use is therefore not recommended. 
    1123  
    1124927 
    1125928\newpage  
     
    1134937a few kilometres in the horizontal, a few meters in the vertical and a few minutes. 
    1135938They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. 
    1136 The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) must be represented entirely in terms of large-scale patterns to close the equations. 
     939The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) 
     940must be represented entirely in terms of large-scale patterns to close the equations. 
    1137941These effects appear in the equations as the divergence of turbulent fluxes 
    1138942(\ie fluxes associated with the mean correlation of small scale perturbations). 
     
    1144948 
    1145949The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions. 
    1146 Therefore subgrid-scale physics \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in 
     950Therefore subgrid-scale physics \textbf{D}$^{\vect U}$, $D^{S}$ and $D^{T}$  in 
    1147951\autoref{eq:PE_dyn}, \autoref{eq:PE_tra_T} and \autoref{eq:PE_tra_S} are divided into 
    1148 a lateral part  \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and 
    1149 a vertical part  \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. 
     952a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and 
     953a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$. 
    1150954The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 
    1151955 
     
    1160964Turbulent motions are thus never explicitly solved, even partially, but always parameterized. 
    1161965The vertical turbulent fluxes are assumed to depend linearly on the gradients of large-scale quantities 
    1162 (for example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$, 
    1163 where $A^{vT}$ is an eddy coefficient). 
     966(for example, the turbulent heat flux is given by $\overline{T' w'} = -A^{v T} \partial_z \overline T$, 
     967where $A^{v T}$ is an eddy coefficient). 
    1164968This formulation is analogous to that of molecular diffusion and dissipation. 
    1165969This is quite clearly a necessary compromise: considering only the molecular viscosity acting on 
     
    1169973\begin{equation} 
    1170974  \label{eq:PE_zdf} 
    1171   \begin{split} 
    1172     {\vect{D}}^{v \vect{U}} &=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\ 
    1173     D^{vT}                         &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ , 
    1174     \quad 
    1175     D^{vS}=\frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial S}{\partial z}} \right) 
    1176   \end{split} 
     975  \begin{gathered} 
     976    \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \ , \\ 
     977          D^{vT}       = \pd[]{z} \lt( A^{vT} \pd[T]{z}         \rt) \quad \text{and} \quad 
     978          D^{vS}       = \pd[]{z} \lt( A^{vT} \pd[S]{z}         \rt) 
     979  \end{gathered} 
    1177980\end{equation} 
    1178981where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. 
     
    12321035one uses a advective scheme which is diffusive enough to maintain the model stability. 
    12331036It must be emphasised that then, all the sub-grid scale physics is included in the formulation of 
    1234 the advection scheme.  
     1037the advection scheme. 
    12351038 
    12361039All these parameterisations of subgrid scale physics have advantages and drawbacks. 
     
    12461049\begin{equation} 
    12471050  \label{eq:PE_iso_tensor} 
    1248   D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    1249   \mbox{with}\quad \;\;\Re =\left( {{ 
    1250         \begin{array}{*{20}c} 
    1251           1 \hfill & 0 \hfill & {-r_1 } \hfill \\ 
    1252           0 \hfill & 1 \hfill & {-r_2 } \hfill \\ 
    1253           {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 
    1254         \end{array} 
    1255       }} \right) 
     1051  D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad 
     1052  \Re = 
     1053    \begin{pmatrix} 
     1054      1    & 0    & -r_1          \\ 
     1055      0    & 1    & -r_2          \\ 
     1056      -r_1 & -r_2 & r_1^2 + r_2^2 \\ 
     1057    \end{pmatrix} 
    12561058\end{equation} 
    1257 where $r_1 \;\mbox{and}\;r_2 $ are the slopes between the surface along which the diffusive operator acts and 
    1258 the model level ($e. g.$ $z$- or $s$-surfaces). 
     1059where $r_1$ and $r_2$ are the slopes between the surface along which the diffusive operator acts and 
     1060the model level (\eg $z$- or $s$-surfaces). 
    12591061Note that the formulation \autoref{eq:PE_iso_tensor} is exact for 
    12601062the rotation between geopotential and $s$-surfaces, 
     
    12671069 
    12681070For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. 
    1269 $\Re $ reduces to the identity in the horizontal direction, no rotation is applied.  
     1071$\Re$ reduces to the identity in the horizontal direction, no rotation is applied. 
    12701072 
    12711073For \textit{geopotential} diffusion, 
     
    12781080\begin{equation} 
    12791081  \label{eq:PE_iso_slopes} 
    1280   r_1 =\frac{e_3 }{e_1 }   \left( \pd[\rho]{i} \right) \left( \pd[\rho]{k} \right)^{-1} \, \quad 
    1281   r_2 =\frac{e_3 }{e_2 }   \left( \pd[\rho]{j} \right) \left( \pd[\rho]{k} \right)^{-1} \, 
     1082  r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
     1083  r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
    12821084\end{equation} 
    12831085while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. 
    12841086 
    12851087\subsubsection{Eddy induced velocity} 
     1088 
    12861089When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used, 
    12871090an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: 
    12881091\[ 
    12891092  % \label{eq:PE_iso+eiv} 
    1290   D^{lT}=\nabla \cdot \left( {A^{lT}\;\Re \;\nabla T} \right) 
    1291   +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right) 
     1093  D^{lT} = \nabla \cdot \lt( A^{lT} \; \Re \; \nabla T \rt) + \nabla \cdot \lt( \vect U^\ast \, T \rt) 
    12921094\] 
    1293 where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent, 
     1095where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a non-divergent, 
    12941096eddy-induced transport velocity. This velocity field is defined by: 
    12951097\[ 
    12961098  % \label{eq:PE_eiv} 
    1297   \begin{split} 
    1298     u^\ast  &= +\frac{1}{e_3       }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_1 } \right] \\ 
    1299     v^\ast  &= +\frac{1}{e_3       }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_2 } \right] \\ 
    1300     w^\ast &=  -\frac{1}{e_1 e_2 }\left[ 
    1301       \frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right) 
    1302       +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right)      \right] 
    1303   \end{split} 
     1099  u^\ast &=   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_1 \rt) \\ 
     1100  v^\ast &=   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_2 \rt) \\ 
     1101  w^\ast &= - \frac{1}{e_1 e_2} \lt[   \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt) 
     1102                                     + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt] 
    13041103\] 
    13051104where $A^{eiv}$ is the eddy induced velocity coefficient 
    13061105(or equivalently the isoneutral thickness diffusivity coefficient), 
    1307 and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes between isoneutral and \emph{geopotential} surfaces. 
     1106and $\tilde r_1$ and $\tilde r_2$ are the slopes between isoneutral and \textit{geopotential} surfaces. 
    13081107Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate:  
    13091108\begin{align} 
    13101109  \label{eq:PE_slopes_eiv} 
    13111110  \tilde{r}_n = 
    1312   \begin{cases} 
    1313     r_n            &      \text{in $z$-coordinate}    \\ 
    1314     r_n + \sigma_n &      \text{in \zstar and $s$-coordinates} 
    1315   \end{cases} 
    1316                      \quad \text{where } n=1,2 
     1111    \begin{cases} 
     1112      r_n            & \text{in $z$-coordinate}                \\ 
     1113      r_n + \sigma_n & \text{in \zstar- and $s$-coordinates} 
     1114    \end{cases} 
     1115  \quad \text{where~} n = 1, 2 
    13171116\end{align} 
    13181117 
    13191118The normal component of the eddy induced velocity is zero at all the boundaries. 
    1320 This can be achieved in a model by tapering either the eddy coefficient or 
    1321 the slopes to zero in the vicinity of the boundaries. 
     1119This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of 
     1120the boundaries. 
    13221121The latter strategy is used in \NEMO (cf. \autoref{chap:LDF}). 
    13231122 
     
    13271126\[ 
    13281127  % \label{eq:PE_bilapT} 
    1329   D^{lT}= - \Delta \left( \;\Delta T \right) 
    1330   \qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right) 
     1128  D^{lT}= - \Delta \; (\Delta T) \quad \text{where} \quad 
     1129  \Delta \bullet = \nabla \lt( \sqrt{B^{lT}} \; \Re \; \nabla \bullet \rt) 
    13311130\] 
    13321131It is the Laplacian operator given by \autoref{eq:PE_iso_tensor} applied twice with 
    1333 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one.  
    1334  
     1132the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 
    13351133 
    13361134\subsubsection{Lateral Laplacian momentum diffusive operator} 
     
    13381136The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by 
    13391137applying \autoref{eq:PE_lap_vector} to the horizontal velocity vector (see \autoref{apdx:B}): 
    1340 \[ 
     1138\begin{align*} 
    13411139  % \label{eq:PE_lapU} 
    1342   \begin{split} 
    1343     {\rm {\bf D}}^{l{\rm {\bf U}}} 
    1344     &= \quad \  \nabla _h \left( {A^{lm}\chi } \right) 
    1345     \ - \ \nabla _h \times \left( {A^{lm}\,\zeta \;{\rm {\bf k}}} \right)     \\ 
    1346     &=   \left( 
    1347       \begin{aligned} 
    1348         \frac{1}{e_1      } \frac{\partial \left( A^{lm} \chi          \right)}{\partial i} 
    1349         &-\frac{1}{e_2 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial j}  \\ 
    1350         \frac{1}{e_2      }\frac{\partial \left( {A^{lm} \chi         } \right)}{\partial j} 
    1351         &+\frac{1}{e_1 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial i} 
    1352       \end{aligned} 
    1353     \right) 
    1354   \end{split} 
    1355 \] 
     1140  \vect D^{l \vect U} &=   \nabla_h        \big( A^{lm}    \chi             \big) 
     1141                         - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\ 
     1142                      &= \lt(   \frac{1}{e_1}     \pd[ \lt( A^{lm}    \chi      \rt) ]{i} \rt. 
     1143                              - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j} 
     1144                                \frac{1}{e_2}     \pd[ \lt( A^{lm}    \chi      \rt) ]{j} 
     1145                         \lt. + \frac{1}{e_1 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{i} \rt) 
     1146\end{align*} 
    13561147 
    13571148Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields 
     
    13591150Unfortunately, it is only available in \textit{iso-level} direction. 
    13601151When a rotation is required 
    1361 (\ie geopotential diffusion in $s-$coordinates or isoneutral diffusion in both $z$- and $s$-coordinates), 
    1362 the $u$ and $v-$fields are considered as independent scalar fields, so that the diffusive operator is given by: 
    1363 \[ 
     1152(\ie geopotential diffusion in $s$-coordinates or isoneutral diffusion in both $z$- and $s$-coordinates), 
     1153the $u$ and $v$-fields are considered as independent scalar fields, so that the diffusive operator is given by: 
     1154\begin{gather*} 
    13641155  % \label{eq:PE_lapU_iso} 
    1365   \begin{split} 
    1366     D_u^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla u} \right) \\ 
    1367     D_v^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla v} \right) 
    1368   \end{split} 
    1369 \] 
     1156    D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \\ 
     1157    D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt) 
     1158\end{gather*} 
    13701159where $\Re$ is given by \autoref{eq:PE_iso_tensor}. 
    13711160It is the same expression as those used for diffusive operator on tracers. 
    13721161It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, 
    1373 \ie on a $f-$ or $\beta-$plane, not on the sphere. 
     1162\ie on a $f$- or $\beta$-plane, not on the sphere. 
    13741163It is also a very good approximation in vicinity of the Equator in 
    13751164a geographical coordinate system \citep{Lengaigne_al_JGR03}. 
     
    13861175 
    13871176\end{document} 
    1388  
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex

    r10442 r10501  
    44 
    55% ================================================================ 
    6 % Chapter 2 Time Domain (step.F90) 
     6% Chapter 2 ——— Time Domain (step.F90) 
    77% ================================================================ 
    8 \chapter{Time Domain (STP) } 
     8\chapter{Time Domain (STP)} 
    99\label{chap:STP} 
    10  
    1110\minitoc 
    1211 
     
    3534\begin{equation} 
    3635  \label{eq:STP} 
    37   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t-\rdt,\,t,\,t+\rdt} 
     36  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t - \rdt, \, t, \, t + \rdt} 
    3837\end{equation}  
    3938where $x$ stands for $u$, $v$, $T$ or $S$; 
     
    8584\begin{equation} 
    8685  \label{eq:STP_asselin} 
    87   x_F^t  = x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right] 
    88 \end{equation}  
     86  x_F^t = x^t + \gamma \, \lt[ x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt] 
     87\end{equation} 
    8988where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient. 
    9089$\gamma$ is initialized as \np{rn\_atfp} (namelist parameter). 
    91 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec:STP_mLF}), 
     90Its default value is \np{rn\_atfp}~\forcode{= 10.e-3} (see \autoref{sec:STP_mLF}), 
    9291causing only a weak dissipation of high frequency motions (\citep{Farge1987}). 
    9392The addition of a time filter degrades the accuracy of the calculation from second to first order. 
     
    111110(when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : 
    112111\[ 
    113   % \label{eq:STP_euler} 
    114    x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \ {D_x}^{t-\rdt} 
    115 \]  
     112  %\label{eq:STP_euler} 
     113  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ D_x^{t - \rdt} 
     114\] 
    116115 
    117116This is diffusive in time and conditionally stable. 
     
    119118\begin{equation} 
    120119  \label{eq:STP_euler_stability} 
    121   A^h < \left\{ 
    122     \begin{aligned} 
    123       &\frac{e^2}{  8 \; \rdt }  &&\quad \text{laplacian diffusion}  \\ 
    124       &\frac{e^4}{64 \; \rdt }   &&\quad \text{bilaplacian diffusion} 
    125     \end{aligned} 
    126   \right. 
     120  A^h < 
     121  \begin{cases} 
     122    \frac{e^2}{ 8 \, \rdt} & \text{laplacian diffusion} \\ 
     123    \frac{e^4}{64 \, \rdt} & \text{bilaplacian diffusion} 
     124  \end{cases} 
    127125\end{equation} 
    128126where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. 
     
    134132but usually the numerical stability condition imposes a strong constraint on the time step. 
    135133Two solutions are available in \NEMO to overcome the stability constraint: 
    136 $(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp}\forcode{ = .true.}) or 
    137 $(b)$ a backward (or implicit) time differencing scheme                   (\np{ln\_zdfexp}\forcode{ = .false.}). 
    138 In $(a)$, the master time step $\Delta $t is cut into $N$ fractional time steps so that 
     134$(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp}~\forcode{= .true.}) or 
     135$(b)$ a backward (or implicit) time differencing scheme                   (\np{ln\_zdfexp}~\forcode{= .false.}). 
     136In $(a)$, the master time step $\Delta$t is cut into $N$ fractional time steps so that 
    139137the stability criterion is reduced by a factor of $N$. 
    140138The computation is performed as follows: 
    141 \[ 
     139\begin{alignat*}{2} 
    142140  % \label{eq:STP_ts} 
    143   \begin{split} 
    144     & x_\ast ^{t-\rdt} = x^{t-\rdt}   \\ 
    145     & x_\ast ^{t-\rdt+L\frac{2\rdt}{N}}=x_\ast ^{t-\rdt+\left( {L-1} 
    146       \right)\frac{2\rdt}{N}}+\frac{2\rdt}{N}\;\text{DF}^{t-\rdt+\left( {L-1} \right)\frac{2\rdt}{N}} 
    147     \quad \text{for $L=1$ to $N$}      \\ 
    148     & x^{t+\rdt} = x_\ast^{t+\rdt} 
    149   \end{split} 
    150 \] 
     141  &x_\ast^{t - \rdt}                      &= &x^{t - \rdt} \\ 
     142  &x_\ast^{t - \rdt + L \frac{2 \rdt}{N}} &=   &x_\ast ^{t - \rdt + (L - 1) \frac{2 \rdt}{N}} 
     143                                             + \frac{2 \rdt}{N} \; DF^{t - \rdt + (L - 1) \frac{2 \rdt}{N}} 
     144  \quad \text{for $L = 1$ to $N$} \\ 
     145  &x^{t + \rdt}                           &= &x_\ast^{t + \rdt} 
     146\end{alignat*} 
    151147with DF a vertical diffusion term. 
    152148The number of fractional time steps, $N$, is given by setting \np{nn\_zdfexp}, (namelist parameter). 
     
    154150\begin{equation} 
    155151  \label{eq:STP_imp} 
    156   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t+\rdt} 
    157 \end{equation}  
     152  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t + \rdt} 
     153\end{equation} 
    158154 
    159155%%gm 
     
    166162\[ 
    167163  % \label{eq:STP_imp_zdf} 
    168   \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\rdt}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta 
    169   _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta_{k+1/2} \left[ {T^{t+1}} \right]} 
    170   \right] 
     164  \frac{T(k)^{t + 1} - T(k)^{t - 1}}{2 \; \rdt} 
     165  \equiv 
     166  \text{RHS} + \frac{1}{e_{3t}} \delta_k \lt[ \frac{A_w^{vT}}{e_{3w} } \delta_{k + 1/2} \lt[ T^{t + 1} \rt] \rt] 
    171167\] 
    172168where RHS is the right hand side of the equation except for the vertical diffusion term. 
     
    174170\begin{equation} 
    175171  \label{eq:STP_imp_mat} 
    176   -c(k+1)\;T^{t+1}(k+1) + d(k)\;T^{t+1}(k) - \;c(k)\;T^{t+1}(k-1) \equiv b(k) 
     172  -c(k + 1) \; T^{t + 1}(k + 1) + d(k) \; T^{t + 1}(k) - \; c(k) \; T^{t + 1}(k - 1) \equiv b(k) 
    177173\end{equation} 
    178174where  
    179 \begin{align*} 
     175\begin{align*}  
    180176  c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k)     \\ 
    181   d(k) &= e_{3t} (k)       \, / \, (2\rdt) + c_k + c_{k+1}    \\ 
    182   b(k) &= e_{3t} (k) \; \left( T^{t-1}(k) \, / \, (2\rdt) + \text{RHS} \right)     
     177  d(k) &= e_{3t}   (k)       \, / \, (2 \rdt) + c_k + c_{k + 1}    \\ 
     178  b(k) &= e_{3t}   (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt) 
    183179\end{align*} 
    184180 
     
    189185(see for example \citet{Richtmyer1967}). 
    190186 
    191  
    192  
    193187% ------------------------------------------------------------------------------------------------------------- 
    194188%        Surface Pressure gradient 
     
    203197\begin{figure}[!t] 
    204198  \begin{center} 
    205     \includegraphics[width=0.7\textwidth]{Fig_TimeStepping_flowchart} 
     199    \includegraphics[]{Fig_TimeStepping_flowchart} 
    206200    \caption{ 
    207201      \protect\label{fig:TimeStep_flowchart} 
    208202      Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}. 
    209       The use of a semi-implicit computation of the hydrostatic pressure gradient requires the tracer equation to 
     203      The use of a semi -implicit computation of the hydrostatic pressure gradient requires the tracer equation to 
    210204      be stepped forward prior to the momentum equation. 
    211205      The need for knowledge of the vertical scale factor (here denoted as $h$) requires the sea surface height and 
     
    225219\label{sec:STP_mLF} 
    226220 
    227 Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 
     221Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to 
     222ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 
    228223The modifications affect both the forcing and filtering treatments in the LF-RA scheme. 
    229224 
    230225In a classical LF-RA environment, the forcing term is centred in time, 
    231 \ie it is time-stepped over a $2\rdt$ period: 
    232 $x^t  = x^t + 2\rdt Q^t $ where $Q$ is the forcing applied to $x$, 
     226\ie it is time-stepped over a $2 \rdt$ period: 
     227$x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, 
    233228and the time filter is given by \autoref{eq:STP_asselin} so that $Q$ is redistributed over several time step. 
    234229In the modified LF-RA environment, these two formulations have been replaced by: 
    235 \begin{align} 
    236   x^{t+\rdt}  &= x^{t-\rdt} + \rdt \left( Q^{t-\rdt/2} + Q^{t+\rdt/2} \right)                   \label{eq:STP_forcing} \\ 
    237   % 
    238   x_F^t  &= x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right] 
    239            - \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]                          \label{eq:STP_RA} 
    240 \end{align} 
     230\begin{gather} 
     231  \label{eq:STP_forcing} 
     232  x^{t + \rdt} = x^{t - \rdt} + \rdt \lt( Q^{t - \rdt / 2} + Q^{t + \rdt / 2} \rt)  \\ 
     233  \label{eq:STP_RA} 
     234  x_F^t       = x^t + \gamma \, \lt( x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt) 
     235                    - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) 
     236\end{gather} 
    241237The change in the forcing formulation given by \autoref{eq:STP_forcing} (see \autoref{fig:MLF_forcing}) 
    242238has a significant effect: 
    243 the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}.  
     239the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}. 
    244240% forcing seen by the model.... 
    245241This property improves the LF-RA scheme in two respects. 
     
    250246Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, 
    251247the modified formulation becomes conservative \citep{Leclair_Madec_OM09}. 
    252 Second, the LF-RA becomes a truly quasi-second order scheme. 
     248Second, the LF-RA becomes a truly quasi -second order scheme. 
    253249Indeed, \autoref{eq:STP_forcing} used in combination with a careful treatment of static instability 
    254250(\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}), 
    255251the two other main sources of time step divergence, 
    256 allows a reduction by two orders of magnitude of the Asselin filter parameter.  
     252allows a reduction by two orders of magnitude of the Asselin filter parameter. 
    257253 
    258254Note that the forcing is now provided at the middle of a time step: 
    259 $Q^{t+\rdt/2}$ is the forcing applied over the $[t,t+\rdt]$ time interval. 
     255$Q^{t + \rdt / 2}$ is the forcing applied over the $[t,t + \rdt]$ time interval. 
    260256This and the change in the time filter, \autoref{eq:STP_RA}, 
    261257allows an exact evaluation of the contribution due to the forcing term between any two time steps, 
     
    265261\begin{figure}[!t] 
    266262  \begin{center} 
    267     \includegraphics[width=0.90\textwidth]{Fig_MLF_forcing} 
     263    \includegraphics[]{Fig_MLF_forcing} 
    268264    \caption{ 
    269265      \protect\label{fig:MLF_forcing} 
     
    271267      (top) ''Traditional'' formulation: 
    272268      the forcing is defined at the same time as the variable to which it is applied 
    273       (integer value of the time step index) and it is applied over a $2\rdt$ period. 
     269      (integer value of the time step index) and it is applied over a $2 \rdt$ period. 
    274270      (bottom)  modified formulation: 
    275271      the forcing is defined in the middle of the time (integer and a half value of the time step index) and 
    276       the mean of two successive forcing values ($n-1/2$, $n+1/2$) is applied over a $2\rdt$ period. 
     272      the mean of two successive forcing values ($n - 1 / 2$, $n + 1 / 2$) is applied over a $2 \rdt$ period. 
    277273    } 
    278274  \end{center} 
     
    297293\] 
    298294This is done simply by keeping the leapfrog environment (\ie the \autoref{eq:STP} three level time stepping) but 
    299 setting all $x^0$ (\textit{before}) and $x^{1}$ (\textit{now}) fields equal at the first time step and 
     295setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and 
    300296using half the value of $\rdt$. 
    301297 
     
    305301running the model for $2N$ time steps in one go, 
    306302or by performing two consecutive experiments of $N$ steps with a restart. 
    307 This requires saving two time levels and many auxiliary data in the restart files in machine precision.  
    308  
    309 Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure gradient 
     303This requires saving two time levels and many auxiliary data in the restart files in machine precision. 
     304 
     305Note that when a semi -implicit scheme is used to evaluate the hydrostatic pressure gradient 
    310306(see \autoref{subsec:DYN_hpg_imp}), an extra three-dimensional field has to 
    311307be added to the restart file to ensure an exact restartability. 
    312 This is done optionally via the \np{nn\_dynhpg\_rst} namelist parameter, 
     308This is done optionally via the  \np{nn\_dynhpg\_rst} namelist parameter, 
    313309so that the size of the restart file can be reduced when restartability is not a key issue 
    314310(operational oceanography or in ensemble simulations for seasonal forecasting). 
    315311 
    316312Note the size of the time step used, $\rdt$, is also saved in the restart file. 
    317 When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed.  
    318 Options are defined through the \ngn{namrun} namelist variables. 
     313When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed. 
     314Options are defined through the  \ngn{namrun} namelist variables. 
    319315%%% 
    320316\gmcomment{ 
     
    360356 
    361357\begin{flalign*} 
    362   &\frac{\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\rdt} 
    363   \equiv \text{RHS}+ \delta_k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]} 
    364   \right]      \\ 
    365   &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1} 
    366   \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]} 
    367   \right]      \\ 
    368   &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1} 
     358  &\frac{\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}}{2\rdt} 
     359  \equiv \text{RHS}+ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]} 
     360  \rt]      \\ 
     361  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1} 
     362  \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]} 
     363  \rt]      \\ 
     364  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1} 
    369365  \equiv 2\rdt \ \text{RHS} 
    370   + 2\rdt \ \left\{ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} [ T_{k+1}^{t+1} - T_k      ^{t+1} ] 
    371     - \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k       ^{t+1} - T_{k-1}^{t+1} ]  \right\}     \\ 
     366  + 2\rdt \ \lt\{ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} [ T_{k +1}^{t+1} - T_k      ^{t+1} ] 
     367    - \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} [ T_k       ^{t+1} - T_{k -1}^{t+1} ]  \rt\}     \\ 
    372368  &\\ 
    373   &\left( e_{3t}\,T \right)_k^{t+1} 
    374   -  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}                  T_{k+1}^{t+1} 
    375   + {2\rdt} \ \left\{  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} 
    376     +  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}     \right\}   T_{k    }^{t+1} 
    377   -  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                  T_{k-1}^{t+1}      \\ 
    378   &\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}    \\ 
    379   %  
     369  &\lt( e_{3t}\,T \rt)_k^{t+1} 
     370  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}                  T_{k +1}^{t+1} 
     371  + {2\rdt} \ \lt\{  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} 
     372    +  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}     \rt\}   T_{k    }^{t+1} 
     373  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                  T_{k -1}^{t+1}      \\ 
     374  &\equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}    \\ 
     375  % 
    380376\end{flalign*} 
    381  
    382377\begin{flalign*} 
    383378  \allowdisplaybreaks 
    384379  \intertext{ Tracer case } 
    385380  % 
    386   &  \qquad \qquad  \quad   -  {2\rdt}                  \ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} 
    387   \qquad \qquad \qquad  \qquad  T_{k+1}^{t+1}   \\ 
    388   &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} 
    389   +   \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\ 
    390   & \qquad \qquad  \qquad \qquad \qquad \quad \ \ -  {2\rdt} \                          \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                          \quad \ \ T_{k-1}^{t+1} 
    391   \ \equiv \ \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  \\ 
     381  &  \qquad \qquad  \quad   -  {2\rdt}                  \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} 
     382  \qquad \qquad \qquad  \qquad  T_{k +1}^{t+1}   \\ 
     383  &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} 
     384  +   \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\ 
     385  & \qquad \qquad  \qquad \qquad \qquad \quad \ \ -  {2\rdt} \                          \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                          \quad \ \ T_{k -1}^{t+1} 
     386  \ \equiv \ \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  \\ 
    392387  % 
    393388\end{flalign*} 
     
    396391  \intertext{ Tracer content case } 
    397392  % 
    398   & -  {2\rdt} \              & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_{k+1}^{t+1}}  && \  \left( e_{3t}\,T \right)_{k+1}^{t+1}   &\\ 
    399   & + {2\rdt} \ \left[ 1  \right.+ & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_k^{t+1}} 
    400   + & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \left.  \right]  & \left( e_{3t}\,T \right)_{k   }^{t+1}  &\\ 
    401   & -  {2\rdt} \               & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_{k-1}^{t+1}}     &\  \left( e_{3t}\,T \right)_{k-1}^{t+1} 
    402   \equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  & 
     393  & -  {2\rdt} \              & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_{k +1}^{t+1}}  && \  \lt( e_{3t}\,T \rt)_{k +1}^{t+1}   &\\ 
     394  & + {2\rdt} \ \lt[ 1  \rt.+ & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_k^{t+1}} 
     395  + & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \lt.  \rt]  & \lt( e_{3t}\,T \rt)_{k   }^{t+1}  &\\ 
     396  & -  {2\rdt} \               & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_{k -1}^{t+1}}     &\  \lt( e_{3t}\,T \rt)_{k -1}^{t+1} 
     397  \equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  & 
    403398\end{flalign*} 
    404399 
    405400%% 
    406401} 
    407 %% 
     402 
    408403\biblio 
    409404 
Note: See TracChangeset for help on using the changeset viewer.