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

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

Ignore:
Timestamp:
2008-02-09T15:13:48+01:00 (16 years ago)
Author:
gm
Message:

trunk - update including Steven correction of the first 5 chapters (until DYN) and activation of Appendix A & B

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/DOC/BETA/Chapters/Chap_Model_Basics.tex

    r707 r817  
    2222 
    2323 
    24 The ocean is a fluid that can be described to a good approximation by the primitive equations, i.e. the Navier-Stokes equations along with a nonlinear equation of state which couples the two active tracers (temperature and salinity) to the fluid velocity, plus the following additional assumptions made from scale considerations: 
     24The ocean is a fluid that can be described to a good approximation by the primitive equations, $i.e.$ the Navier-Stokes equations along with a nonlinear equation of state which couples the two active tracers (temperature and salinity) to the fluid velocity, plus the following additional assumptions made from scale considerations: 
    2525 
    2626\textit{(1) spherical earth approximation: }the geopotential surfaces are assumed to be spheres so that gravity (local vertical) is parallel to the earth's radius 
    … …  
    3232\textit{(4) Boussinesq hypothesis:} density variations are neglected except in their contribution to the buoyancy force 
    3333 
    34 \textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a balance between the vertical pressure gradient and buoyancy force (this removes convective processes from  
    35 the initial Navier-Stokes equations: they must be parameterized) 
     34\textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a balance between the vertical pressure gradient and the buoyancy force (this removes convective processes from  
     35the initial Navier-Stokes equations and so convective processes must be parameterized instead) 
    3636 
    3737\textit{(6) Incompressibility hypothesis: }the three dimensional divergence of the velocity vector is assumed to be zero. 
    3838 
    39 Because the gravitational force is so dominant in the equations of large-scale motions, it is quite useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, i.e. tangent to the geopotential surfaces. Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ (the subscript $h$ denotes the local horizontal vector, i.e. over the (\textbf{i},\textbf{j}) plan), $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system provides the following six equations (namely the momentum balance, the hydrostatic equilibrium, the incompressibility, the heat and salt conservation and an equation of state): 
     39Because the gravitational force is so dominant in the equations of large-scale motions, it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, $i.e.$ tangent to the geopotential surfaces. Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ (the subscript $h$ denotes the local horizontal vector, $i.e.$ over the (\textbf{i},\textbf{j}) plane), $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system provides the following six equations (namely the momentum balance, the hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation equations and an equation of state): 
    4040\begin{subequations} \label{Eq_PE} 
    4141  \begin{equation}     \label{Eq_PE_dyn} 
    … …  
    4444            +\frac{1}{2}\nabla \left( {{\rm {\bf U}}^2} \right)}    \right]_h 
    4545 -f\;{\rm {\bf k}}\times {\rm {\bf U}}_h  
    46 -\frac{1}{\rho _o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} 
     46-\frac{1}{\rho _o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} + {\rm {\bf F}}^{\rm {\bf U}} 
    4747  \end{equation} 
    4848  \begin{equation}     \label{Eq_PE_hydrostatic} 
    … …  
    5353  \end{equation} 
    5454\begin{equation} \label{Eq_PE_tra_T} 
    55 \frac{\partial T}{\partial t} = - \nabla \cdot  \left( T \ \rm{\bf U} \right) + D^T 
     55\frac{\partial T}{\partial t} = - \nabla \cdot  \left( T \ \rm{\bf U} \right) + D^T + F^T 
    5656  \end{equation} 
    5757  \begin{equation}     \label{Eq_PE_tra_S} 
    58 \frac{\partial S}{\partial t} = - \nabla \cdot  \left( S \ \rm{\bf U} \right) + D^S 
     58\frac{\partial S}{\partial t} = - \nabla \cdot  \left( S \ \rm{\bf U} \right) + D^S + F^S 
    5959  \end{equation} 
    6060  \begin{equation}     \label{Eq_PE_eos} 
    … …  
    6262  \end{equation} 
    6363\end{subequations} 
    64 where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, $t$ the time, $z$ the vertical coordinate, $\rho $ the \textit{in situ} density given by the equation of state (\ref{Eq_PE_eos}), $\rho_o$ a reference density, $p$ the pressure, $f=2 \bf \Omega \cdot \bf k$ the Coriolis acceleration (where $\bf \Omega$ is the Earth angular velocity vector), and $g$ the gravitational acceleration. ${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterizations of small-scale physics for momentum, temperature and salinity, including surface forcing terms. Their nature and formulation are discussed in \S\ref{PE_zdf_ldf}, page~\pageref{PE_zdf_ldf}. 
     64where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, $t$ is the time, $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by the equation of state (\ref{Eq_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, $f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration (where $\bf \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration.  
     65${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterizations of small-scale physics for momentum, temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$ surface forcing terms. Their nature and formulation are discussed in \S\ref{PE_zdf_ldf} and page \S\ref{PE_boundary_condition}. 
    6566 
    6667. 
    … …  
    7273\label{PE_boundary_condition} 
    7374 
    74 An ocean is bounded by complex coastlines and bottom topography at its base and by an air-sea or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ the height of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, chosen as a mean sea surface (Fig.~\ref{Fig_ocean_bc}). Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth, the continental surfaces, the sea ice and the atmosphere. However, some of these fluxes are so weak that even on climatic time scales of thousands of years they can be neglected. In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and the other components of the earth system. 
     75An ocean is bounded by complex coastlines, bottom topography at its base and an air-sea or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, chosen as a mean sea surface (Fig.~\ref{Fig_ocean_bc}). Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth, the continental margins, the sea ice and the atmosphere. However, some of these fluxes are so weak that even on climatic time scales of thousands of years they can be neglected. In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and the other components of the earth system. 
    7576 
    7677%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    8182%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    8283 
     84 
    8385\begin{description} 
    84 \item[Land - ocean interface:] the major flux between continental surfaces and the ocean is a mass exchange of fresh water through river runoff. Such an exchange modifies locally the sea surface salinity especially in the vicinity of major river mouths. It can be neglected for short range integrations but has to be taken into account for long term integrations as it influences the characteristics of water masses formed (especially at high latitudes). It is required to close the water cycle of the climatic system. It is usually specified as a fresh water flux at the air-sea interface in the vicinity of river mouths. 
    85 \item[Solid earth - ocean interface:] heat and salt fluxes across the sea floor are negligibly small, except in special areas of little extent. They are always neglected in the model. The boundary condition is thus set to no flux of heat and salt across solid boundaries. For momentum, the situation is different. There is no flow across solid boundaries, i.e. the velocity normal to the ocean bottom and coastlines is zero (in other words, the bottom velocity is parallel to solid boundaries). This kinematic boundary condition can be expressed as: 
     86\item[Land - ocean interface:] the major flux between continental margins and the ocean is a mass exchange of fresh water through river runoff. Such an exchange modifies the sea surface salinity especially in the vicinity of major river mouths. It can be neglected for short range integrations but has to be taken into account for long term integrations as it influences the characteristics of water masses formed (especially at high latitudes). It is required in order to close the water cycle of the climate system. It is usually specified as a fresh water flux at the air-sea interface in the vicinity of river mouths. 
     87\item[Solid earth - ocean interface:] heat and salt fluxes through the sea floor are small, except in special areas of little extent. They are usually neglected in the model \footnote{In fact, it has been shown that the heat flux associated with the solid Earth cooling ($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world ocean (see \ref{TRA_bbc}).}. The boundary condition is thus set to no flux of heat and salt across solid boundaries. For momentum, the situation is different. There is no flow across solid boundaries, $i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words, the bottom velocity is parallel to solid boundaries). This kinematic boundary condition can be expressed as: 
    8688\begin{equation} \label{Eq_PE_w_bbc} 
    8789w = -{\rm {\bf U}}_h \cdot  \nabla _h \left( H \right) 
    8890\end{equation} 
    89 In addition, the ocean exchanges momentum with the earth through friction processes. Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized in terms of turbulent fluxes through bottom and/or lateral boundary conditions. Its specification depends on the nature of the physical parameterization used for ${\rm {\bf D}}^{\rm {\bf U}}$ in \eqref{Eq_PE_dyn}. They are discussed in in \S\ref{PE_zdf}, page~\pageref{PE_zdf}.% and Chap. III.6 to 9. 
     91In addition, the ocean exchanges momentum with the earth through frictional processes. Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. Its specification depends on the nature of the physical parameterization used for ${\rm {\bf D}}^{\rm {\bf U}}$ in \eqref{Eq_PE_dyn}. It is discussed in \S\ref{PE_zdf}, page~\pageref{PE_zdf}.% and Chap. III.6 to 9. 
    9092\item[Atmosphere - ocean interface:] the kinematic surface condition plus the mass flux of fresh water PE  (the precipitation minus evaporation budget) leads to:  
    9193\begin{equation} \label{Eq_PE_w_sbc} 
    … …  
    9597\end{equation} 
    9698The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system) leads to the continuity of pressure across the interface $z=\eta$. The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. 
    97 \item[Sea ice - ocean interface:] the two media exchange heat, salt, fresh water and momentum. The sea surface temperature is constrained to be at the freezing point at the interface. Sea ice salinity is very low ($\sim5 \,psu$) compared to those of the ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected. 
     99\item[Sea ice - ocean interface:] the ocean and sea ice exchange heat, salt, fresh water and momentum. The sea surface temperature is constrained to be at the freezing point at the interface. Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected. 
    98100\end{description} 
    99101 
    … …  
    111113\label{PE_p_formulation} 
    112114 
    113 The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a reference geopotential surface ($z=0$) and an hydrostatic pressure $p_h$ such that: $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\ref{Eq_PE_hydrostatic}), assuming that pressure in decibars can be approximated by depth in meters in (\ref{Eq_PE_eos}). The hydrostatic pressure is then given by: 
     115The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\ref{Eq_PE_hydrostatic}), assuming that pressure in decibars can be approximated by depth in meters in (\ref{Eq_PE_eos}). The hydrostatic pressure is then given by: 
    114116\begin{equation} \label{Eq_PE_pressure} 
    115117p_h \left( {i,j,z,t} \right) 
    116118 = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,z} \right)\;d\varsigma }  
    117119\end{equation} 
    118 The surface pressure requires a more specific treatment. Two strategies can be considered: $(a)$ the introduction of a new variable $\eta$, the free-surface elevation, for which a prognostic equation can be established and solved; $(b)$ the assumption that the ocean surface is a rigid lid, on which the pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used, a solution of the free-surface elevation consists in the excitation of external gravity waves. The flow is barotropic and the surface moves up and down with gravity as the restoring force. The phase speed of such waves is high (some hundreds of metres per second) so that the time step would have to be very short if they were present in the model. The latter strategy filters these waves as the rigid lid approximation implies $\eta=0$, i.e. the sea surface is the surface $z=0$. This well known approximation increases the surface wave speed to infinity and modifies certain other longwave dynamics (e.g. barotropic Rossby or planetary waves). In the present release of OPA, both strategies are still available. They are further described in the next two sub-sections. 
     120 Two strategies can be considered for the surface pressure term: $(a)$ introduce of a new variable $\eta$, the free-surface elevation, for which a prognostic equation can be established and solved; $(b)$ assume that the ocean surface is a rigid lid, on which the pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used, one solution of the free-surface elevation consists of the excitation of external gravity waves. The flow is barotropic and the surface moves up and down with gravity as the restoring force. The phase speed of such waves is high (some hundreds of metres per second) so that the time step would have to be very short if they were present in the model. The latter strategy filters out these waves since the rigid lid approximation implies $\eta=0$, $i.e.$ the sea surface is the surface $z=0$. This well known approximation increases the surface wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic Rossby or planetary waves). In the present release of \NEMO, both strategies are still available. They are further described in the next two sub-sections. 
    119121 
    120122% ------------------------------------------------------------------------------------------------------------- 
    … …  
    124126\label{PE_free_surface} 
    125127 
    126 In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced which describes the shape of the air-sea interface. This variable is solution of a prognostic equation which is established by forming the vertical average of the kinematics surface condition (\ref{Eq_PE_w_bbc}): 
     128In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced which describes the shape of the air-sea interface. This variable is solution of a prognostic equation which is established by forming the vertical average of the kinematic surface condition (\ref{Eq_PE_w_bbc}): 
    127129\begin{equation} \label{Eq_PE_ssh} 
    128130\frac{\partial \eta }{\partial t}=-D+P-E 
    … …  
    138140$\bullet$ If one is interested in EGWs, in particular the tides and their interaction  
    139141with the baroclinic structure of the ocean (internal waves) possibly in  
    140 shallow seas, then a non linear free surface is the most adequate: this  
     142shallow seas, then a non linear free surface is the most appropriate. This  
    141143means that no approximation is made in (\ref{Eq_PE_ssh}) and that the variation of the  
    142144ocean volume is fully taken into account. Note that in order to study the  
    … …  
    147149 
    148150$\bullet$ If one is not interested in EGW but rather sees them as high frequency  
    149 noise, it is possible to apply a filter to slow down the fastest waves while  
    150 not altering the slow barotropic Rossby waves. In that case it is also  
    151 generally sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which allows  
     151noise, it is possible to apply an explicit filter to slow down the fastest waves while  
     152not altering the slow barotropic Rossby waves. If further, an approximative conservation  
     153of heat and salt contents is sufficient for the problem solved, then it is  
     154sufficient to solve a linearized version of (\ref{Eq_PE_ssh}), which still allows  
    152155to take into account freshwater fluxes applied at the ocean surface \citep{Roullet2000}. 
    153156 
    154157$\bullet$ For process studies not involving external waves nor surface freshwater  
    155158fluxes, it is possible to use the rigid lid approximation see (next  
    156 section). The ocean surface is considered as a fixed surface, so that all  
     159section). The ocean surface is then considered as a fixed surface, so that all  
    157160external waves are removed from the system.  
    158161 
    159 The filtering of EGWs in models with a free surface is usually a matter of discretisation of the temporal derivatives, using the time splitting method \citep{Killworth1991, Zhang1992} or the implicit scheme \citep{Dukowicz1994}. In OPA, we use a slightly different approach developed by \citet{Roullet2000}: the damping of EGWs is ensured by introducing an additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes:  
     162The filtering of EGWs in models with a free surface is usually a matter of discretisation of the temporal derivatives, using the time splitting method \citep{Killworth1991, Zhang1992} or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach developed by \citet{Roullet2000}: the damping of EGWs is ensured by introducing an additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes:  
    160163\begin{equation} \label{Eq_PE_flt} 
    161164\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} 
    … …  
    163166- g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)  
    164167\end{equation} 
    165 where $T_c$, is a parameter homogeneous to a time which characterizes the force, $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and  represents the collected contributions of the Coriolis, hydrostatic pressure gradient, non-linear and viscous terms in \eqref{Eq_PE_dyn}. 
    166  
    167 The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagates, $i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than $T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs can be damped by choosing $T_c > \Delta t$. \citet{Roullet2000} demonstrate that (\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which has to be computed implicitly. This is not surprising since the use of a large time step has a necessarily numerical cost. Two gains arise in comparison with the previous formulations. Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as soon as $T_c > \Delta t$. 
    168  
    169 When the variations of free surface elevation are small compared to the thickness of the model layers, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized by \citet{Roullet2000} the linearization of (\ref{Eq_PE_ssh}) has consequences on the conservation of salt in the model. With the nonlinear free surface equation, the time evolution of the total salt content is  
     168where $T_c$, is a parameter with dimensions of time which characterizes the force, $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, non-linear and viscous terms in \eqref{Eq_PE_dyn}. 
     169 
     170The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate, $i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than $T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs can be damped by choosing $T_c > \Delta t$. \citet{Roullet2000} demonstrate that (\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which has to be computed implicitly. This is not surprising since the use of a large time step has a necessarily numerical cost. Two gains arise in comparison with the previous formulations. Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as soon as $T_c > \Delta t$. 
     171 
     172When the variations of free surface elevation are small compared to the thickness of the first model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized by \citet{Roullet2000} the linearization of (\ref{Eq_PE_ssh}) has consequences on the conservation of salt in the model. With the nonlinear free surface equation, the time evolution of the total salt content is  
    170173\begin{equation} \label{Eq_PE_salt_content} 
    171174\frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} =\int\limits_S  
    172175{S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds} 
    173176\end{equation} 
    174 where $S$ is the salinity, and the total salt is integrated in the whole ocean volume $D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) is satisfied, so that the salt is perfectly conserved. When the free surface equation is linearized, \citet{Roullet2000} show that the total salt content integrated in the fixed volume $D$ (bounded by the surface $z=0$) is no longer conserved: 
     177where $S$ is the salinity, and the total salt is integrated over the whole ocean volume $D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) is satisfied, so that the salt is perfectly conserved. When the free surface equation is linearized, \citet{Roullet2000} show that the total salt content integrated in the fixed volume $D$ (bounded by the surface $z=0$) is no longer conserved: 
    175178\begin{equation} \label{Eq_PE_salt_content_linear} 
    176179\frac{\partial }{\partial t}\int\limits_D {S\;dv} =-\int\limits_S  
    … …  
    186189 
    187190Although the total salt content is not exactly conserved with the linearized free surface, its variations are driven by correlations of the time variation of surface salinity with the sea surface height, which is a negligible term. This situation contrasts with  
    188 the case of the rigid lid approximation (following section) in which case freshwater forcing is represented by a virtual salt flux, leading to spurious sources or sinks of salt \citep{Roullet2000}. 
     191the case of the rigid lid approximation (following section) in which case freshwater forcing is represented by a virtual salt flux, leading to a spurious source of salt at the ocean surface \citep{Roullet2000}. 
    189192 
    190193% ------------------------------------------------------------------------------------------------------------- 
    … …  
    194197\label{PE_rigid_lid} 
    195198 
    196 With the rigid lid approximation, we assume that the ocean surface ($z=0$) is a rigid lid on which a pressure $p_s$ is exerted. This implies that the vertical velocity at the surface is equal to zero. From the continuity equation (\ref{Eq_PE_continuity}) and the kinematic condition at the bottom (\ref{Eq_PE_w_bbc}) (no flux across the bottom), it can be shown that the vertically integrated flow $H{\rm {\bf \bar {U}}}_h$ is non-divergent (where the overbar indicates a vertical average over the whole water column, i.e. from $z=-H$ , the ocean bottom, to $z=0$ , the rigid-lid). Thus, can be derived from a volume transport streamfunction $\psi$: 
     199With the rigid lid approximation, we assume that the ocean surface ($z=0$) is a rigid lid on which a pressure $p_s$ is exerted. This implies that the vertical velocity at the surface is equal to zero. From the continuity equation \eqref{Eq_PE_continuity} and the kinematic condition at the bottom \eqref{Eq_PE_w_bbc} (no flux across the bottom), it can be shown that the vertically integrated flow $H{\rm {\bf \bar {U}}}_h$ is non-divergent (where the overbar indicates a vertical average over the whole water column, i.e. from $z=-H$, the ocean bottom, to $z=0$ , the rigid-lid). Thus, $\rm {\bf \bar {U}}_h$ can be derived from a volume transport streamfunction $\psi$: 
    197200\begin{equation} \label{Eq_PE_u_psi} 
    198201\overline{\vect{U}}_h =\frac{1}{H}\left(   \vect{k} \times \nabla \psi   \right) 
    … …  
    215218\end{equation} 
    216219 
    217 Using the proper boundary conditions, (\ref{Eq_PE_psi}) can be solved to find $\partial_t \psi$ and thus using (\ref{Eq_PE_u_barotrope}) the horizontal surface pressure gradient. It should be noted that $p_s$ can be computed by taking the divergence of (\ref{Eq_PE_u_barotrope}) and solving the resulting elliptic equation. Thus the surface pressure is a diagnostic quantity that can be recovered for analysis purposes. 
    218  
    219 A difficulty lies in the determination of the boundary condition on $\partial_t \psi$. The boundary condition on velocity is that there is no flow normal to a solid wall, i.e. the coastlines are streamlines. Therefore (I.2.7 is solved with the following Dirichlet boundary condition: $\partial_t \psi$ is constant along each coastline of the same continent or of the same island. When all the coastlines are connected (there are no islands), the constant value of $\partial_t \psi$ along the coast can be arbitrarily chosen to be zero. When islands are present in the domain, the value of the barotropic streamfunction will generally be different for each island and for the continent, and will vary with respect to time. So the boundary condition is: $\psi=0$ along the continent and $\psi=\mu_n$ along island $n$ ($1 \leq n \leq Q$), where $Q$ is the number of islands present in the domain and $\mu_n$ is a time dependent variable. A time evolution equation of the unknown $\mu_n$ can be found by evaluating the circulation of the time derivative of the vertical average (barotropic) velocity field along a closed contour around each island. Since the circulation of a  
    220 gradient field along a closed contour is zero, from (\ref{Eq_PE_u_barotrope}) we have: 
     220Using the proper boundary conditions, \eqref{Eq_PE_psi} can be solved to find $\partial_t \psi$ and thus using \eqref{Eq_PE_u_barotrope} the horizontal surface pressure gradient. It should be noted that $p_s$ can be computed by taking the divergence of \eqref{Eq_PE_u_barotrope} and solving the resulting elliptic equation. Thus the surface pressure is a diagnostic quantity that can be recovered for analysis purposes. 
     221 
     222A difficulty lies in the determination of the boundary condition on $\partial_t \psi$. The boundary condition on velocity is that there is no flow normal to a solid wall, $i.e.$ the coastlines are streamlines. Therefore \eqref{Eq_PE_psi} is solved with the following Dirichlet boundary condition: $\partial_t \psi$ is constant along each coastline of the same continent or of the same island. When all the coastlines are connected (there are no islands), the constant value of $\partial_t \psi$ along the coast can be arbitrarily chosen to be zero. When islands are present in the domain, the value of the barotropic streamfunction will generally be different for each island and for the continent, and will vary with respect to time. So the boundary condition is: $\psi=0$ along the continent and $\psi=\mu_n$ along island $n$ ($1 \leq n \leq Q$), where $Q$ is the number of islands present in the domain and $\mu_n$ is a time dependent variable. A time evolution equation of the unknown $\mu_n$ can be found by evaluating the circulation of the time derivative of the vertical average (barotropic) velocity field along a closed contour around each island. Since the circulation of a gradient field along a closed contour is zero, from \eqref{Eq_PE_u_barotrope} we have: 
    221223\begin{equation} \label{Eq_PE_isl_circulation} 
    222224\oint_n {\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \left(  
    … …  
    226228\end{equation} 
    227229 
    228 Since (\ref{Eq_PE_psi}) is linear, its solution \textit{$\psi $} can be decomposed as follows: 
     230Since (\ref{Eq_PE_psi}) is linear, its solution \textit{$\psi $} can be decomposed  
     231as follows: 
    229232\begin{equation} \label{Eq_PE_psi_isl} 
    230233\psi =\psi _o +\sum\limits_{n=1}^{n=Q} {\mu _n \psi _n }  
    231234\end{equation} 
    232 where $\psi _o$ is the solution of (\ref{Eq_PE_psi}) with $\psi _o=0$ long all the coastlines, and where $\psi _n$ is the solution of (\ref{Eq_PE_psi}) with the right-hand side equal to $0$, and with $\psi _n =1$ long the island $n$, $\psi _n =0$ along the other boundaries. The function $\psi _n$ is thus independent of time. Introducing (\ref{Eq_PE_psi_isl}) into (\ref{Eq_PE_isl_circulation}) yields: 
     235where $\psi _o$ is the solution of \eqref{Eq_PE_psi} with $\psi _o=0$ long all  
     236the coastlines, and where $\psi _n$ is the solution of \eqref{Eq_PE_psi} with  
     237the right-hand side equal to $0$, and with $\psi _n =1$ long the island $n$,  
     238$\psi _n =0$ along the other boundaries. The function $\psi _n$ is thus  
     239independent of time. Introducing \eqref{Eq_PE_psi_isl} into  
     240\eqref{Eq_PE_isl_circulation} yields: 
    233241\begin{multline} \label{Eq_PE_psi_isl_circulation} 
    234 \left[ {\oint_n {\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \psi _m }  
    235 \right]\cdot {\rm {\bf d}}\ell } } \right]_{1\leq m\leqslant Q \atop 1\leq n\leqslant Q  } 
     242\left[ {\oint_n {\frac{1}{H}  \left[ {{\rm {\bf k}}\times \nabla \psi _m } \right]\cdot {\rm {\bf d}}\ell } } \right]_{1\leq m\leqslant Q \atop 1\leq n\leqslant Q  } 
    236243 \left( {\frac{\partial \mu _n }{\partial t}}  
    237244\right)_{1\leqslant n\leqslant Q}        \\ 
    … …  
    246253\right)_{1\leqslant n\leqslant Q} ={\rm {\bf B}} 
    247254\end{equation} 
    248 where \textbf{A} is a $Q  \times Q$ matrix and \textbf{B} is a time dependent vector. As \textbf{A} is independent of time, it can be calculated and inverted once. The time derivative of the  streamfunction when islands are present is thus given by: 
     255where \textbf{A} is a $Q  \times Q$ matrix and \textbf{B} is a time dependent vector. As \textbf{A} is independent of time, it can be calculated and inverted once. The time derivative of the streamfunction when islands are present is thus given by: 
    249256\begin{equation} \label{Eq_PE_psi_isl_dt} 
    250257\frac{\partial \psi }{\partial t}=\frac{\partial \psi _o }{\partial  
    … …  
    261268 
    262269 
    263  
    264270% ------------------------------------------------------------------------------------------------------------- 
    265271% Tensorial Formalism 
    … …  
    268274\label{PE_tensorial} 
    269275 
    270 In many ocean circulation problems, the flow field has regions of enhanced dynamics (i.e. surface layers, western boundary currents, equatorial currents, or ocean fronts). The representation of such dynamical processes can be improved by specifically increasing the model resolution in these regions. As well, it may be convenient to use a lateral boundary-following coordinate system to better represent coastal dynamics. Moreover, the common geographical coordinate system has a singular point at the North Pole that cannot be easily treated in a global model without filtering. A solution consists in introducing an appropriate coordinate transformation that shifts the singular point on land \citep{MadecImb1996, Murray1996}. As a conclusion, it is important to solve the primitive equations in various curvilinear coordinate systems. An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism. This formalism is suited to any multidimensional curvilinear coordinate system. Ocean modellers mainly use three-dimensional orthogonal grids on the sphere, with conservation of the local vertical. Here we give the simplified equations for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics. 
    271  
    272 Let (\textbf{i},\textbf{j},\textbf{k}) be a set of orthogonal curvilinear coordinates on the sphere associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, i.e. along geopotential surfaces (\ref{Fig_referential}). Let $(\lambda,\varphi,z)$ be the geographical coordinates system in which a position is defined by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of the  
    273 earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea level (\ref{Fig_referential}). The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, the three scale factors: 
     276In many ocean circulation problems, the flow field has regions of enhanced dynamics ($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts). The representation of such dynamical processes can be improved by specifically increasing the model resolution in these regions. As well, it may be convenient to use a lateral boundary-following coordinate system to better represent coastal dynamics. Moreover, the common geographical coordinate system has a singular point at the North Pole that cannot be easily treated in a global model without filtering. A solution consists of introducing an appropriate coordinate transformation that shifts the singular point onto land \citep{MadecImb1996, Murray1996}. As a consequence, it is important to solve the primitive equations in various curvilinear coordinate systems. An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism. This formalism is suited to any multidimensional curvilinear coordinate system. Ocean modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth approximation), with preservation of the local vertical. Here we give the simplified equations for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics. 
     277 
     278Let (\textbf{i},\textbf{j},\textbf{k}) be a set of orthogonal curvilinear coordinates on the sphere associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, $i.e.$ along geopotential surfaces (Fig.\ref{Fig_referential}). Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of the earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea level (Fig.\ref{Fig_referential}). The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, the three scale factors: 
    274279\begin{equation} \label{Eq_scale_factors} 
    275280\begin{aligned} 
    … …  
    325330\end{equation} 
    326331\end{subequations} 
    327 where $q$ is a scalar quantity and ${\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in  
    328 the $(i,j,k)$ coordinate system. 
     332where $q$ is a scalar quantity and ${\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinate system. 
    329333 
    330334% ------------------------------------------------------------------------------------------------------------- 
    … …  
    334338\label{PE_zco_Eq} 
    335339 
    336 In order to express the primitive equations in tensorial formalism, it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using \eqref{Eq_PE_grad}) to \eqref{Eq_PE_lap_vector}. Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 
     340In order to express the Primitive Equations in tensorial formalism, it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using \eqref{Eq_PE_grad}) to \eqref{Eq_PE_lap_vector}. Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 
    337341\begin{equation} \label{Eq_PE_curl_Uh} 
    338342\zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v}  
    … …  
    346350\end{equation} 
    347351 
    348 Using the fact that horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that $e_3$  is a function of the single variable $k$, the nonlinear term of \eqref{Eq_PE_dyn} can be transformed as follows: 
     352Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that $e_3$  is a function of the single variable $k$, the nonlinear term of \eqref{Eq_PE_dyn} can be transformed as follows: 
    349353\begin{flalign*} 
    350354&\left[ {\left( { \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}} 
    … …  
    352356\end{flalign*} 
    353357\begin{flalign*} 
    354 &=\left( {{\begin{array}{*{20}c} 
     358&\qquad=\left( {{\begin{array}{*{20}c} 
    355359 {\left[    {   \frac{1}{e_3} \frac{\partial u  }{\partial k} 
    356360         -\frac{1}{e_1} \frac{\partial w  }{\partial i} } \right] w - \zeta \; v }     \\ 
    … …  
    364368\end{flalign*} 
    365369\begin{flalign*} 
    366 &=\left( {{ \begin{array}{*{20}c} 
     370& \qquad =\left( {{  \begin{array}{*{20}c} 
    367371 {-\zeta \; v} \hfill \\ 
    368372 { \zeta \; u} \hfill \\ 
    … …  
    394398\end{equation} 
    395399 
    396 This is the so-called \textit{vector invariant form} of the momentum advection. For some purposes, it can be advantageous to write this term in the so-called flux form, i.e. to write it as the divergence of fluxes. For example, the first component of (\ref{Eq_PE_vector_form}) (the $i$-component) is transformed as follows: 
     400This is the so-called \textit{vector invariant form} of the momentum advection term. For some purposes, it can be advantageous to write this term in the so-called flux form, $i.e.$ to write it as the divergence of fluxes. For example, the first component of \eqref{Eq_PE_vector_form} (the $i$-component) is transformed as follows: 
    397401\begin{flalign*} 
    398402&{ \begin{array}{*{20}l} 
    399403\left[ {\left( {\nabla \times \vect{U}} \right)\times \vect{U} 
    400           +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i    \\ 
    401 \\ 
     404          +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i   % \\ 
     405%\\ 
    402406     = - \zeta \;v  
    403407     + \frac{1}{2\;e_1 } \frac{\partial \left( {u^2+v^2} \right)}{\partial i} 
    404408     + \frac{1}{e_3}w \ \frac{\partial u}{\partial k}          \\ 
    405409\\ 
    406 =\frac{1}{e_1 \; e_2} \left(  -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i} 
     410\qquad =\frac{1}{e_1 \; e_2} \left(    -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i} 
    407411                     +v\frac{\partial \left( {e_1 \,u} \right)}{\partial j}    \right) 
    408412+\frac{1}{e_1 e_2 }\left(  +e_2 \; u\frac{\partial u}{\partial i} 
    … …  
    413417\begin{flalign*} 
    414418&{ \begin{array}{*{20}l} 
    415 =\frac{1}{e_1 \; e_2}   \left\{  
     419\qquad =\frac{1}{e_1 \; e_2}  \left\{  
    416420 -\left(        v^2  \frac{\partial e_2                                }{\partial i}  
    417421      +e_2 \,v    \frac{\partial v                                   }{\partial i}     \right) 
    … …  
    424428                  \right\}  
    425429+\frac{1}{e_3} \left( 
    426                \frac{\partial \left( {w\,v} \right)         }{\partial k} 
     430               \frac{\partial \left( {w\,u} \right)         }{\partial k} 
    427431       -u         \frac{\partial w                    }{\partial k}  \right) \\ 
    428432\end{array} }     & 
    … …  
    430434\begin{flalign*} 
    431435&{ \begin{array}{*{20}l} 
    432 =\frac{1}{e_1 \; e_2}   \left(  
     436\qquad =\frac{1}{e_1 \; e_2}  \left(  
    433437               \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 
    434438      +        \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j}  \right) 
    435 +\frac{1}{e_3 }      \frac{\partial \left( {w\,v       } \right)}{\partial k} 
     439+\frac{1}{e_3 }      \frac{\partial \left( {w\,u       } \right)}{\partial k} 
    436440\\  \qquad \qquad \quad 
    437441+\frac{1}{e_1 e_2 }     \left(  
    … …  
    446450\begin{flalign*} 
    447451&{ \begin{array}{*{20}l} 
    448 = \nabla \cdot \left( {{\rm {\bf U}}\,u}  \right) 
    449 -  \nabla \cdot {\rm {\bf U}} \ u 
     452\qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right) 
     453-   \left( \nabla \cdot {\rm {\bf U}} \right) \ u 
    450454+\frac{1}{e_1 e_2 }\left(  
    451455      -v^2     \frac{\partial e_2 }{\partial i} 
    … …  
    456460\begin{flalign*} 
    457461&{ \begin{array}{*{20}l} 
    458 = \nabla \cdot \left(   {{\rm {\bf U}}\,u}      \right) 
     462\qquad = \nabla \cdot \left(  {{\rm {\bf U}}\,u}      \right) 
    459463+  \frac{1}{e_1 e_2 }   \left( v \; \frac{\partial e_2}{\partial i} 
    460464                         -u \; \frac{\partial e_1}{\partial j}    \right)  \left( -v \right)  
    … …  
    462466\end{flalign*} 
    463467 
    464 The flux form of the momentum advection is therefore given by: 
     468The flux form of the momentum advection term is therefore given by: 
    465469\begin{multline} \label{Eq_PE_flux_form} 
    466470      \left[  
    … …  
    468472+\frac{1}{2}   \nabla \left(  {{\rm {\bf U}}^2}    \right) 
    469473      \right]_h  
    470 =  \nabla \cdot   \left( {{\begin{array}{*{20}c}   {\rm {\bf U}} \, u   \hfill \\ 
     474\\ 
     475= \nabla \cdot    \left( {{\begin{array}{*{20}c}   {\rm {\bf U}} \, u   \hfill \\ 
    471476                                    {\rm {\bf U}} \, v   \hfill \\ 
    472477                  \end{array} }}     
    473478            \right) 
    474 \\ 
    475479+\frac{1}{e_1 e_2 }     \left(  
    476480       v\frac{\partial e_2}{\partial i} 
    … …  
    479483\end{multline} 
    480484 
    481 The flux form has two terms, the first one is expressed as the divergence of momentum fluxes (so the flux form name given to this formulation) and the second one is due to the curvilinear nature of the coordinate system used. The latter is called the \emph{metric} term and can be viewed as a modification of the Coriolis parameter:  
     485The flux form has two terms, the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) and the second one is due to the curvilinear nature of the coordinate system used. The latter is called the \emph{metric} term and can be viewed as a modification of the Coriolis parameter:  
    482486\begin{equation} \label{Eq_PE_cor+metric} 
    483487f \to f + \frac{1}{e_1 \; e_2}   \left(    v \frac{\partial e_2}{\partial i} 
    … …  
    485489\end{equation} 
    486490 
    487 Note that in the case of geographical coordinate, i.e. when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of the Coriolis parameter $f \to f+(u/a) \tan \varphi$. 
    488  
    489 The equations solved by the ocean model can be written in the following tensorial formalism: 
    490  
    491 $\bullet$ vector invariant form of the momentum equations : 
     491Note that in the case of geographical coordinate, $i.e.$ when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of the Coriolis parameter $f \to f+(u/a) \tan \varphi$. 
     492 
     493To sum up, the equations solved by the ocean model can be written in the following tensorial formalism: 
     494 
     495\vspace{+10pt} 
     496$\bullet$ \textit{momentum equations} : 
     497 
     498vector invariant form : 
    492499\begin{subequations} \label{Eq_PE_dyn_vect} 
    493 \begin{equation} \label{Eq_PE_dyn_vect_u} 
     500\begin{multline} \label{Eq_PE_dyn_vect_u} 
    494501\frac{\partial u}{\partial t}= 
    495 +   \left( {\zeta +f} \right)\,v   -   \frac{1}{e_3 }w\frac{\partial u}{\partial k} 
    496 -   \frac{1}{e_1 }\frac{\partial}{\partial i}\left(  
    497       {\frac{1}{2}\left( {u^2+v^2} \right)+\frac{p_h+p_s }{\rho _o }}    \right) 
    498 +   D_u^{\vect{U}}  
    499 \end{equation} 
    500 \begin{equation} \label{Eq_PE_dyn_vect_v} 
     502   +   \left( {\zeta +f} \right)\,v                                     
     503   -   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right)  
     504   -   \frac{1}{e_3} w \frac{\partial u}{\partial k}       \\ 
     505   -   \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho _o}    \right)     
     506+   D_u^{\vect{U}}  +   F_u^{\vect{U}} 
     507\end{multline} 
     508\begin{multline} \label{Eq_PE_dyn_vect_v} 
    501509\frac{\partial v}{\partial t}= 
    502 -   \left( {\zeta +f} \right)\,u   -   \frac{1}{e_3 }w\frac{\partial v}{\partial k} 
    503 -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left(  
    504       {\frac{1}{2}\left( {u^2+v^2} \right)+\frac{p_h+p_s }{\rho _o }}    \right) 
    505 +  D_v^{\vect{U}}  
    506 \end{equation} 
     510   -   \left( {\zeta +f} \right)\,u    
     511   -   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^2  \right)        -   \frac{1}{e_3 }w\frac{\partial v}{\partial k}         \\ 
     512   -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}  \right)     
     513+  D_v^{\vect{U}}  +   F_v^{\vect{U}} 
     514\end{multline} 
    507515\end{subequations} 
    508516 
    509 where $\zeta$ is given by (\ref{Eq_PE_curl_Uh}) and the surface pressure gradient is given by: 
     517flux form: 
     518\begin{subequations} \label{Eq_PE_dyn_flux} 
     519\begin{multline} \label{Eq_PE_dyn_flux_u} 
     520\frac{\partial u}{\partial t}= 
     521+   \left( { f + \frac{1}{e_1 \; e_2} 
     522               \left(    v \frac{\partial e_2}{\partial i} 
     523                  -u \frac{\partial e_1}{\partial j}  \right)}    \right) \, v    \\ 
     524- \frac{1}{e_1 \; e_2}  \left(  
     525               \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 
     526      +        \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j}  \right) 
     527- \frac{1}{e_3 }     \frac{\partial \left( {w\,u       } \right)}{\partial k}    \\ 
     528-   \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho _o}   \right) 
     529+   D_u^{\vect{U}} +   F_u^{\vect{U}} 
     530\end{multline} 
     531\begin{multline} \label{Eq_PE_dyn_flux_v} 
     532\frac{\partial v}{\partial t}= 
     533-   \left( { f + \frac{1}{e_1 \; e_2} 
     534               \left(    v \frac{\partial e_2}{\partial i} 
     535                  -u \frac{\partial e_1}{\partial j}  \right)}    \right) \, u   \\ 
     536 \frac{1}{e_1 \; e_2}   \left(  
     537               \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i} 
     538      +        \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j}  \right) 
     539- \frac{1}{e_3 }     \frac{\partial \left( {w\,v       } \right)}{\partial k}    \\ 
     540-   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}    \right) 
     541+  D_v^{\vect{U}} +  F_v^{\vect{U}}  
     542\end{multline} 
     543\end{subequations} 
     544where $\zeta$ is given by \eqref{Eq_PE_curl_Uh} and the surface pressure gradient formulation depends on the one of the free surface: 
    510545 
    511546$*$ free surface formulation 
    … …  
    528563\end{equation} 
    529564where ${\vect{M}}= \left( M_u,M_v \right)$ represents the collected contributions of nonlinear,  
    530 viscous and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn_vect} and the overbar indicates a vertical average over the whole water column ($i.e.$ from $z=-H$, the ocean bottom, to $z=0$, the rigid-lid), and where the time derivative of $\psi$ is the solution of an elliptic equation: 
     565viscosity and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn_vect} and the overbar indicates a vertical average over the whole water column ($i.e.$ from $z=-H$, the ocean bottom, to $z=0$, the rigid-lid), and where the time derivative of $\psi$ is the solution of an elliptic equation: 
    531566\begin{multline} \label{Eq_psi_total} 
    532567  \frac{\partial }{\partial i}\left[ {\frac{e_2 }{H\,e_1}\frac{\partial}{\partial i} 
    … …  
    547582\end{equation} 
    548583 
    549 where the divergence of the horizontal velocity, $\chi$ is given by (I.3.8). 
    550  
    551 $\bullet$ tracer equations: 
     584where the divergence of the horizontal velocity, $\chi$ is given by \eqref{Eq_PE_div_Uh}. 
     585 
     586\vspace{+10pt} 
     587$\bullet$ \textit{tracer equations} : 
    552588\begin{equation} \label{Eq_S} 
    553 \frac{\partial T}{\partial t}=-\frac{1}{e_1 e_2 }\left[ {\frac{\partial  
    554 \left( {e_2 T\,u} \right)}{\partial i}+\frac{\partial \left( {e_1 T\,v}  
    555 \right)}{\partial j}} \right]-\frac{1}{e_3 }\frac{\partial \left( {T\,w}  
    556 \right)}{\partial k}+D^T 
     589\frac{\partial T}{\partial t} =  
     590-\frac{1}{e_1 e_2 }\left[ {      \frac{\partial \left( {e_2 T\,u} \right)}{\partial i} 
     591                  +\frac{\partial \left( {e_1 T\,v} \right)}{\partial j}} \right] 
     592-\frac{1}{e_3 }\frac{\partial \left( {T\,w} \right)}{\partial k} + D^T + F^T 
    557593\end{equation} 
    558594\begin{equation} \label{Eq_T} 
    559 \frac{\partial S}{\partial t}=-\frac{1}{e_1 e_2 }\left[ {\frac{\partial  
    560 \left( {e_2 S\,u} \right)}{\partial i}+\frac{\partial \left( {e_1 S\,v}  
    561 \right)}{\partial j}} \right]-\frac{1}{e_3 }\frac{\partial \left( {S\,w}  
    562 \right)}{\partial k}+D^S 
     595\frac{\partial S}{\partial t} =  
     596-\frac{1}{e_1 e_2 }\left[    {\frac{\partial \left( {e_2 S\,u} \right)}{\partial i} 
     597                  +\frac{\partial \left( {e_1 S\,v} \right)}{\partial j}} \right] 
     598-\frac{1}{e_3 }\frac{\partial \left( {S\,w} \right)}{\partial k} + D^S + F^S 
    563599\end{equation} 
    564600\begin{equation} \label{Eq_rho} 
    … …  
    566602\end{equation} 
    567603 
    568 The expression of \textbf{D}$^{U}$, $D^{S}$ and$ D^{T}$ depends on the subgrid  
    569 scale parameterization used. It will be defined in \S\ref{PE_zdf}. 
     604The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid  
     605scale parameterization used. It will be defined in \S\ref{PE_zdf}. The nature and formulation of ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, are discussed in Chapter~\ref{SBC}. 
     606 
     607\newpage  
     608% ================================================================ 
     609% Curvilinear z*-coordinate System 
     610% ================================================================ 
     611\section{Curvilinear \textit{z*}-coordinate System} 
     612 
     613%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     614\begin{figure}[!b] \label{Fig_z_zstar}  \begin{center} 
     615\includegraphics[width=1.0\textwidth]{./Figures/Fig_z_zstar.pdf} 
     616\caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate \citep{Adcroft_Campin_OM04} ).} 
     617\end{center}   \end{figure} 
     618%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     619 
     620 
     621In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. These coordinates systems is presented in a report  
     622\citep{Levier2007} available on the \NEMO web site.  
     623 
     624\gmcomment{ 
     625The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation which allows one to deal with large amplitude free-surface 
     626variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In 
     627the  \textit{z*} formulation, the variation of the column thickness due to sea-surface 
     628undulations is not concentrated in the surface level, as in the $z$-coordinate formulation, 
     629but is equally distributed over the full water column. Thus vertical 
     630levels naturally follow sea-surface variations, with a linear attenuation with 
     631depth, as illustrated by figure fig.1c . Note that with a flat bottom, such as in 
     632fig.1c, the bottom-following  $z$ coordinate and  \textit{z*} are equivalent. 
     633The definition and modified oceanic equations for the rescaled vertical coordinate 
     634 \textit{z*}, including the treatment of fresh-water flux at the surface, are 
     635detailed in Adcroft and Campin (2004). The major points are summarized 
     636here. The position ( \textit{z*}) and vertical discretization (\textit{z*}) are expressed as: 
     637 
     638$H +  \textit{z*} = (H + z) / r$ and  $\delta \textit{z*} = \delta z / r$ with $r = \frac{H+\eta} {H}$ 
     639 
     640Since the vertical displacement of the free surface is incorporated in the vertical 
     641coordinate  \textit{z*}, the upper and lower boundaries are at fixed  \textit{z*} position,  $\textit{z*} = 0$ and  $\textit{z*} = ?H$ respectively. Also the divergence of the flow field is no longer zero as shown by the continuity equation: 
     642 
     643$\frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right) 
     644      \left( r \; w\textit{*} \right) = 0 $ 
     645 
     646} 
     647 
    570648 
    571649\newpage  
    … …  
    581659\subsection{Introduction} 
    582660 
    583 Several important aspects of the ocean circulation are influenced by bottom topography. Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, sills and channels that strongly constrain the path of water masses, but more subtle effects exist. For  
    584 example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. Topographic Rossby waves can be excited and can interact with the mean current. In the $z-$coordinate system presented in the previous section (\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom and to large localized depth gradients associated with large localized vertical velocities. The response to such a velocity field often leads to numerical dispersion effects. 
    585  
    586 A terrain-following coordinate system (hereafter $s-$coordinates) avoids the discretisation error in the depth field since the layers of computation are gradually adjusted with depth to the ocean bottom. Relatively shallow topographic features in the deep ocean, which would be ignored in typical $z-$model applications with the largest grid spacing at greatest depths, can easily be represented (with relatively low vertical resolution) as can gentle, large-scale slopes of the sea floor. A terrain-following model (hereafter $s-$model) also facilitates the modelling of the boundary layer flows over a large depth range, which in the framework of the $z-$model would require high vertical resolution over the whole depth range. Moreover, with $s-$coordinates it is possible, at least in principle, to have the bottom and the sea surface as the only boundaries of the domain. Nevertheless, $s-$coordinates also have its drawbacks. Perfectly adapted to a homogeneous ocean, it has strong limitations as soon as stratification is introduced. The main two problems come from the truncation error in the horizontal pressure gradient and a possibly increased diapycnal diffusion. The horizontal pressure force in $s-$coordinates consists of two terms (see \colorbox{yellow}{Appendix A}), 
     661Several important aspects of the ocean circulation are influenced by bottom topography. Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, sills and channels that strongly constrain the path of water masses, but more subtle effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. Topographic Rossby waves can be excited and can interact with the mean current. In the $z-$coordinate system presented in the previous section (\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom and to large localized depth gradients associated with large localized vertical velocities. The response to such a velocity field often leads to numerical dispersion effects. One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate)  
     662 
     663The $s$-coordinate avoids the discretisation error in the depth field since the layers of computation are gradually adjusted with depth to the ocean bottom. Relatively small topographic features as well as  gentle, large-scale slopes of the sea floor in the deep ocean, which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths, can easily be represented (with relatively low vertical resolution). A terrain-following model (hereafter $s-$model) also facilitates the modelling of the boundary layer flows over a large depth range, which in the framework of the $z$-model would require high vertical resolution over the whole depth range. Moreover, with a $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface as the only boundaries of the domain (nomore lateral boundary condition to specify). Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a homogeneous ocean, it has strong limitations as soon as stratification is introduced. The main two problems come from the truncation error in the horizontal pressure gradient and a possibly increased diapycnal diffusion. The horizontal pressure force in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}), 
    587664 
    588665\begin{equation} \label{Eq_PE_p_sco} 
    … …  
    591668\end{equation} 
    592669 
    593 The second term in (\ref{Eq_PE_p_sco}) depends on the tilt of the coordinate surface and  
    594 introduces a truncation error that is not present in a $z-$model. In the special case of $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), \cite{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error. It depends on topographic slope, stratification, horizontal and vertical resolution, and the finite difference scheme. This error limits the possible topographic slopes that a model can handle at a given horizontal and vertical resolution. This is a severe restriction for large-scale applications using realistic bottom topography. The large-scale slopes require high horizontal resolution, and the computational cost becomes prohibitive. This problem can be, at least partially, overcome by mixing $s-$coordinates and step-like representation of bottom topography \citep{Madec1996}. However, another problem is then raised in the definition of the model domain.  
    595  
    596 \colorbox{yellow}{Aike Beckmann's solution} 
    597  
    598 A minimum of diffusion along the coordinate surfaces of any finite difference model is always required for numerical reasons. It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. This is the case for a $z-$model as well as for a $s-$model.  
    599 However, density varies more strongly on $s-$surfaces than on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal diffusion in a $s-$model than in a $z-$model. Whereas such a diapycnal diffusion in a $z-$model tends to weaken horizontal density (pressure) gradients and thus the horizontal circulation, it usually reinforces these gradients in a $s-$model, creating spurious circulation. For example, imagine an isolated bump of topography in an ocean at rest with a horizontally uniform stratification. Spurious diffusion along $s-$surfaces will induce a bump of isoneutral surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, the ocean will stay at rest in a $z-$model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the  
    600 water column (i.e. the main thermocline) \citep{Madec1996}. An alternate solution consists in rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_zdf} and \colorbox{yellow}{Appendix B}. 
    601  
    602 %AMT 
    603 \amtcomment{  
    604 The $s-$coordinates introduced here \citep{Lott1990,Madec1996}  
    605 \colorbox{yellow}{ to be update $==>$} 
    606 differ mainly in two aspects from  
    607 similar models. It combines the properties which make OPA suitable for  
    608 climate applications with a good representation of bottom topography  
    609 allowing mixed step-like/terrain following topography. It also offers a  
    610 completely general transformation, $s=s(i,j,z)$, for the vertical coordinate which goes beyond those of previous hybrid models except the GFDL version developed by \citep{Gerdes1993a,Gerdes1993b}  
    611 which has similar properties as the OPA release presented here. 
    612 \colorbox{yellow}{ <== end update} 
    613 } 
     670The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface and introduces a truncation error that is not present in a $z$-model. In the special case of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error. It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, and the finite difference scheme. This error limits the possible topographic slopes that a model can handle at a given horizontal and vertical resolution. This is a severe restriction for large-scale applications using realistic bottom topography. The large-scale slopes require high horizontal resolution, and the computational cost becomes prohibitive. This problem can be at least partially overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec1996}. However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for a realistic bottom topography: a envelope topography is defined in $s$-coordinate on which a full or partial step bottom topography is then applied in order to adjust the model depth to the observed one (see \S\ref{DOM_zgr}. 
     671 
     672For numerical reasons a minimum of diffusion is required along the coordinate surfaces of any finite difference model. It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as well as for a $s$-model.  
     673However, density varies more strongly on $s-$surfaces than on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. For example, imagine an isolated bump of topography in an ocean at rest with a horizontally uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column (i.e. the main thermocline) \citep{Madec1996}. An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}. Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}).  
     674 
     675The $s-$coordinates introduced here \citep{Lott1990,Madec1996} differ mainly in two aspects from similar models:  it allows  a representation of bottom topography with mixed full or partial step-like/terrain following topography ; It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. 
     676 
    614677% ------------------------------------------------------------------------------------------------------------- 
    615678% The s-coordinate Formulation 
    … …  
    617680\subsection{The \textit{s-}coordinate Formulation} 
    618681 
    619 Starting from the set of equations established in \S\ref{PE_zco} for the special case $k=z$ and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k)$, which includes $z-$ and $\sigma-$coordinates as special cases ($s=z$ and $s=\sigma=z/H$, resp.). A formal derivation of the transformed equations is given in \colorbox{yellow}{Appendix A}. Let us define the vertical scale factor by $e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k)$ ), and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by : 
     682Starting from the set of equations established in \S\ref{PE_zco} for the special case $k=z$ and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, which includes $z$-, \textit{z*}- and $\sigma-$coordinates as special cases ($s=z$, $s=\textit{z*}$, and $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). A formal derivation of the transformed equations is given in Appendix~\ref{Apdx_A}. Let us define the vertical scale factor by $e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by : 
    620683\begin{equation} \label{Eq_PE_sco_slope} 
    621684\sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s  
    … …  
    623686\sigma _2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s 
    624687\end{equation} 
    625  
    626 We also introduce a "vertical" velocity $\omega$ defined as the velocity normal to $s-$surfaces: 
     688We also introduce  $\omega $, a dia-surface velocity component, defined as the velocity relative to the moving $s$-surfaces and normal to them: 
    627689\begin{equation} \label{Eq_PE_sco_w} 
    628 \omega = w - \sigma_1 \; u - \sigma_2 \; v 
    629 \end{equation} 
    630 The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinates can be written as follows: 
    631  
     690\omega  = w - e_3 \, \frac{\partial s}{\partial t} - \sigma _1 \,u - \sigma _2 \,v    \\ 
     691\end{equation} 
     692 
     693The equations solved by the ocean model \eqref{Eq_PE} in $s-$coordinate can be written as follows: 
     694 
     695 \vspace{0.5cm} 
    632696* momentum equation: 
    633697\begin{multline} \label{Eq_PE_sco_u} 
    634 \frac{\partial u}{\partial t}=+\left( {\zeta +f} \right)\,v-\frac{1}{e_3  
    635 }\omega \frac{\partial u}{\partial k}-\frac{1}{e_1 }\frac{\partial  
    636 }{\partial i}\left( {\frac{1}{2}\left( {u^2+v^2} \right)+\frac{p_h }{\rho _o}} \right)   \\ 
    637 +g\frac{\rho }{\rho _o }\sigma _1 -\frac{1}{\rho _o e_1  
    638 }\frac{\partial p_s }{\partial i}+D_u^{\rm {\bf U}}  
     698\frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t}= 
     699   +   \left( {\zeta +f} \right)\,v                                     
     700   -   \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left(  u^2+v^2   \right)  
     701   -   \frac{1}{e_3} \omega \frac{\partial u}{\partial k}       \\ 
     702   -   \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho _o}    \right)     
     703   +  g\frac{\rho }{\rho _o}\sigma _1  
     704   +   D_u^{\vect{U}}  +   F_u^{\vect{U}} \quad 
    639705\end{multline} 
    640706\begin{multline} \label{Eq_PE_sco_v} 
    641 \frac{\partial v}{\partial t}=-\left( {\zeta +f} \right)\,u-\frac{1}{e_3  
    642 }\omega \frac{\partial v}{\partial k}-\frac{1}{e_2 }\frac{\partial  
    643 }{\partial j}\left( {\frac{1}{2}\left( {u^2+v^2} \right)+\frac{p_h }{\rho _o}} \right)   \\ 
    644 +g\frac{\rho }{\rho _o }\sigma _2 -\frac{1}{\rho _o e_2  
    645 }\frac{\partial p_s }{\partial j}+D_v^{\rm {\bf U}}  
    646 \end{multline} 
    647 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic pressure have the same expressions as in $z-$coordinates although they do not represent exactly the same quantities. $\omega $ is provided by the same equation as w, i.e. (\ref{Eq_w_diag}), with $\chi$, the divergence of the horizontal velocity field given by: 
    648  
    649 \begin{equation} \label{Eq_PE_sco_div} 
     707\frac{1}{e_3} \frac{\partial \left(  e_3\,v  \right) }{\partial t}= 
     708   -   \left( {\zeta +f} \right)\,u    
     709   -   \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left(  u^2+v^2  \right)         
     710   -   \frac{1}{e_3 } \omega \frac{\partial v}{\partial k}         \\ 
     711   -   \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho _o}  \right)  
     712    +  g\frac{\rho }{\rho _o }\sigma _2    
     713   +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad 
     714\end{multline} 
     715where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic pressure have the same expressions as in $z$-coordinates although they do not represent exactly the same quantities. $\omega$ is provided by the continuity equation (see Appendix~\ref{Apdx_A}): 
     716 
     717\begin{equation} \label{Eq_PE_sco_continuity} 
     718\frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0    
     719\qquad \text{with }\;\;   
    650720\chi =\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3 \,u}  
    651721\right)}{\partial i}+\frac{\partial \left( {e_1 e_3 \,v} \right)}{\partial  
    … …  
    653723\end{equation} 
    654724 
     725 \vspace{0.5cm} 
    655726* tracer equations: 
    656  
    657 \begin{equation} \label{Eq_PE_sco_t} 
    658 \frac{\partial T}{\partial t}=-\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial  
    659 \left( {e_2 e_3 T\,u} \right)}{\partial i}+\frac{\partial \left( {e_1 e_3  
    660 T\,v} \right)}{\partial j}} \right]-\frac{1}{e_3 }\frac{\partial \left(  
    661 {T\,\omega } \right)}{\partial k}+D^T 
    662 \end{equation} 
    663  
    664 \begin{equation} \label{Eq_PE_sco_s} 
    665 \frac{\partial S}{\partial t}=-\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial  
    666 \left( {e_2 e_3 S\,u} \right)}{\partial i}+\frac{\partial \left( {e_1 e_3  
    667 S\,v} \right)}{\partial j}} \right]-\frac{1}{e_3 }\frac{\partial \left(  
    668 {S\,\omega } \right)}{\partial k}+D^S 
    669 \end{equation} 
    670  
    671 The equation of state has the same expression as in $z-$coordinates. The expression of \textbf{D}$^{U}$, $D^{S}$ and$ D^{T}$ depends on the subgrid scale parameterization used. It will be defined in \S\ref{PE_zdf_ldf}.  
    672  
     727\begin{multline} \label{Eq_PE_sco_t} 
     728\frac{1}{e_3} \frac{\partial \left(  e_3\,T  \right) }{\partial t}= 
     729-\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,T} \right)}{\partial i} 
     730                                           +\frac{\partial \left( {e_1 e_3\,v\,T} \right)}{\partial j}} \right]   \\ 
     731-\frac{1}{e_3 }\frac{\partial \left( {T\,\omega } \right)}{\partial k}   + D^T + F^S   \qquad 
     732\end{multline} 
     733 
     734\begin{multline} \label{Eq_PE_sco_s} 
     735\frac{1}{e_3} \frac{\partial \left(  e_3\,S  \right) }{\partial t}= 
     736-\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,S} \right)}{\partial i} 
     737                                           +\frac{\partial \left( {e_1 e_3\,v\,S} \right)}{\partial j}} \right]    \\ 
     738-\frac{1}{e_3 }\frac{\partial \left( {S\,\omega } \right)}{\partial k}     + D^S + F^S   \qquad 
     739\end{multline} 
     740 
     741The equation of state has the same expression as in $z$-coordinate, and similar expressions are used for mixing and forcing terms. 
     742 
     743\gmcomment{ 
    673744\colorbox{yellow}{ to be updated $= = >$} 
    674 The whole set of the continuous equations solved by the model in the $s-$coordinate system is summarised in Table~I.2.  
    675  
    676745Add a few works on z and zps and s and underlies the differences between all of them 
    677 \colorbox{yellow}{ $< = =$ end update} 
    678  
    679 \newpage  
    680 % ================================================================ 
    681 % Curvilinear z*- s*-coordinate System 
    682 % ================================================================ 
    683 \section{Curvilinear \textit{z*}- or \textit{s*} coordinate System} 
    684  
    685 % ------------------------------------------------------------------------------------------------------------- 
    686 % ???? 
    687 % ------------------------------------------------------------------------------------------------------------- 
    688  
    689 \colorbox{yellow}{ to be updated $= = >$} 
    690  
    691 In that case, the free surface equation is nonlinear, and the variations of  
    692 volume are fully taken into account. These coordinates systems is presented in a report  
    693 \citep{Levier2007} available on the \NEMO web site.  
    694  
    695 \colorbox{yellow}{ $< = =$ end update} 
     746\colorbox{yellow}{ $< = =$ end update}  } 
    696747 
    697748\newpage  
    … …  
    705756space and time scales larger than a few kilometres in the horizontal, a few  
    706757meters in the vertical and a few minutes. They are usually solved at larger  
    707 scales, the specified grid spacing and time step of the numerical model. The  
     758scales: the specified grid spacing and time step of the numerical model. The  
    708759effects of smaller scale motions (coming from the advective terms in the  
    709760Navier-Stokes equations) must be represented entirely in terms of  
    710761large-scale patterns to close the equations. These effects appear in the  
    711 equations as the divergence of turbulent fluxes (i.e. fluxes associated with  
     762equations as the divergence of turbulent fluxes ($i.e.$ fluxes associated with  
    712763the mean correlation of small scale perturbations). Assuming a turbulent  
    713764closure hypothesis is equivalent to choose a formulation for these fluxes.  
    … …  
    718769 
    719770The control exerted by gravity on the flow induces a strong anisotropy  
    720 between the lateral and vertical motions. Therefore subgrid-scale physics  \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in \eqref{Eq_PE_dyn}, \eqref{Eq_PE_tra_T} and \eqref{Eq_PE_tra_S} are divided into a lateral part  \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part  \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$ . The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 
     771between the lateral and vertical motions. Therefore subgrid-scale physics  \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in \eqref{Eq_PE_dyn}, \eqref{Eq_PE_tra_T} and \eqref{Eq_PE_tra_S} are divided into a lateral part  \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part  \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 
    721772 
    722773% ------------------------------------------------------------------------------------------------------------- 
    … …  
    727778 
    728779The model resolution is always larger than the scale at which the major  
    729 sources of vertical turbulence occurs (shear instability, internal wave  
     780sources of vertical turbulence occur (shear instability, internal wave  
    730781breaking...). Turbulent motions are thus never explicitly solved, even  
    731782partially, but always parameterized. The vertical turbulent fluxes are  
    … …  
    746797   \end{split} 
    747798\end{equation} 
    748 where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must  
    749 be specified (see Chap.~\ref{SBC}). All the vertical physics is embedded in the specification of the eddy coefficients. They can be assumed to be either constant, or function of the local fluid properties (as Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a turbulent closure model. The choices available in OPA are discussed in \S\ref{ZDF}). 
     799where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must be specified (see Chap.~\ref{SBC} and \ref{ZDF} and \S\ref{TRA_bbl}). All the vertical physics is embedded in the specification of the eddy coefficients. They can be assumed to be either constant, or function of the local fluid properties ($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a turbulent closure model. The choices available in \NEMO are discussed in \S\ref{ZDF}). 
    750800 
    751801% ------------------------------------------------------------------------------------------------------------- 
    … …  
    756806 
    757807Lateral turbulence can be roughly divided into a mesoscale turbulence  
    758 associated to eddies which can be solved explicitly if the resolution is  
    759 sufficient as their underlying physics are included in the primitive  
    760 equations, and a sub mesoscale turbulence which is never explicitly solved  
     808associated with eddies (which can be solved explicitly if the resolution is  
     809sufficient since their underlying physics are included in the primitive  
     810equations), and a sub mesoscale turbulence which is never explicitly solved  
    761811even partially, but always parameterized. The formulation of lateral eddy  
    762812fluxes depends on whether the mesoscale is below or above the grid-spacing  
    763 (i.e. the model is eddy-resolving or not). 
    764  
    765 In non-eddy- resolving configurations, the closure is similar to that used  
     813($i.e.$ the model is eddy-resolving or not). 
     814 
     815In non-eddy-resolving configurations, the closure is similar to that used  
    766816for the vertical physics. The lateral turbulent fluxes are assumed to depend  
    767817linearly on the lateral gradients of large-scale quantities. The resulting  
    768818lateral diffusive and dissipative operators are of second order.  
    769819Observations show that lateral mixing induced by mesoscale turbulence tends  
    770 to be along isoneutral surfaces (or more precisely neutral surfaces, i.e.  
    771 isoneutral surfaces referenced at the local depth) rather than across them.  
    772 As the slope of isoneutral surfaces is small in the ocean, a common  
     820to be along isopycnal surfaces (or more precisely neutral surfaces \cite{McDougall1987}) rather than across them.  
     821As the slope of neutral surfaces is small in the ocean, a common  
    773822approximation is to assume that the `lateral' direction is the horizontal,  
    774 i.e. the lateral mixing is performed along geopotential surfaces. This leads  
     823$i.e.$ the lateral mixing is performed along geopotential surfaces. This leads  
    775824to a geopotential second order operator for lateral subgrid scale physics.  
    776825This assumption can be relaxed: the eddy-induced turbulent fluxes can be  
    777826better approached by assuming that they depend linearly on the gradients of  
    778 large-scale quantities computed along isoneutral surfaces. In such a case,  
     827large-scale quantities computed along neutral surfaces. In such a case,  
    779828the diffusive operator is an isoneutral second order operator and it has  
    780829components in the three space directions. However, both horizontal and  
    781 isoneutral operators have no effect on mean (i.e. large scale) potential  
     830isoneutral operators have no effect on mean ($i.e.$ large scale) potential  
    782831energy whereas potential energy is a main source of turbulence (through  
    783832baroclinic instabilities). \citet{Gent1990} have proposed a  
    … …  
    788837operator and an eddy induced advective part. In all these lateral diffusive  
    789838formulations, the specification of the lateral eddy coefficients remains the  
    790 problematic point as there is no satisfactory formulation of these  
     839problematic point as there is no really satisfactory formulation of these  
    791840coefficients as a function of large-scale features. 
    792841 
    … …  
    796845eddies. The role devoted to the subgrid-scale physics is to dissipate the  
    797846energy that cascades toward the grid scale and thus ensures the stability of  
    798 the model while not interfering with the solved mesoscale activity. 
     847the model while not interfering with the solved mesoscale activity. Another approach  
     848is becoming more and more popular: instead of specifying explicitly a sub-grid scale  
     849term in the momentum and tracer time evolution equations, one uses a advective scheme which is diffusive enough to maintain the model stability. It must be emphasised 
     850that then, all the sub-grid scale physics is in this case include in the formulation of the 
     851advection scheme.  
    799852 
    800853All these parameterizations of subgrid scale physics present advantages and  
    801 disadvantages. There are not all available in OPA. In the $z-$coordinate  
    802 formulation, four options are offered for active tracers (temperature and  
     854drawbacks. There are not all available in \NEMO. In the $z$-coordinate  
     855formulation, five options are offered for active tracers (temperature and  
    803856salinity): second order geopotential operator, second order isoneutral  
    804 operator, \citet{Gent1990} parameterization and fourth order  
    805 geopotential operator. The same options are available for momentum, except  
    806 \citet{Gent1990} parameterization which only involves tracers. In  
    807 $s-$coordinate formulation, an additional option is offered for tracers: second  
     857operator, \citet{Gent1990} parameterization, fourth order  
     858geopotential operator, and various slightly diffusive advection schemes. The same options are available for momentum, except  
     859\citet{Gent1990} parameterization which only involves tracers. In the 
     860$s$-coordinate formulation, additional options are offered for tracers: second  
    808861order operator acting along $s-$surfaces, and for momentum: fourth order  
    809862operator acting along $s-$surfaces (see \S\ref{LDF}). 
    … …  
    811864\subsubsection{lateral second order tracer diffusive operator} 
    812865 
    813 The lateral second order tracer diffusive operator is defined by (see \colorbox{yellow}{Appendix B}): 
     866The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}): 
    814867\begin{equation} \label{Eq_PE_iso_tensor} 
    815868D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad  
    … …  
    821874\end{equation} 
    822875where $r_1 \;\mbox{and}\;r_2 $ are the slopes between the surface along  
    823 which the diffusive operator acts and the surface of computation ($z-$ or  
    824 $s-$surfaces), and $r_1$ and $r_2$ is the differential operator defined  
    825 in \S\ref{PE_zco} or \S\ref{PE_sco} depending on the vertical coordinate used. Note that the formulation of $r_1$ and $r_2 $ is exact for the  
    826 slopes between geopotential and $s-$surfaces, while it is only an approximation  
    827 for the slopes between isoneutral and $z$ or $s-$surfaces. Indeed, in the latter  
    828 case, two assumptions are made to simplify $r_1$ and $r_2$ \citep{Cox1987}: the ratio between lateral and vertical diffusive coefficients is  
    829 known to be several orders of magnitude smaller than unity, and the slopes  
    830 are, generally less than $10^{2}$ in the ocean (see \colorbox{yellow}{Appendix B}). This leads to  
    831 the linear tensor (\ref{Eq_PE_iso_tensor}) where the two isoneutral directions of diffusion  
    832 are independent and where the diapycnal diffusivity contribution is solely along the vertical. 
     876which the diffusive operator acts and the model level ($e. g.$ $z$- or  
     877$s$-surfaces). Note that the formulation \eqref{Eq_PE_iso_tensor} is exact for the  
     878rotation between geopotential and $s$-surfaces, while it is only an approximation  
     879for the rotation between isoneutral and $z$- or $s$-surfaces. Indeed, in the latter  
     880case, two assumptions are made to simplify  \eqref{Eq_PE_iso_tensor} \citep{Cox1987}. First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and dia-neutral diffusive coefficients is known to be several orders of magnitude smaller than unity. Second, the two isoneutral directions of diffusion are assumed to be independent since the slopes are generally less than $10^{-2}$ in the ocean (see Appendix~\ref{Apdx_B}). 
    833881 
    834882For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the  
    835 geopotential and computational surfaces: in $z$-coordinates they are zero ($r_1$ and $r_2)$ while in $s-$coordinate they are equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ). 
    836  
    837 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral and computational surfaces. Therefore, they have a same expression in $z-$and $s-$coor-dinates: 
     883geopotential and computational surfaces: in $z$-coordinates they are zero ($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ). 
     884 
     885For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral and computational surfaces. Therefore, they have a same expression in $z$- and $s$-coordinates: 
    838886\begin{equation} \label{Eq_PE_iso_slopes} 
    839887r_1 =\frac{e_3 }{e_1 }  \left( {\frac{\partial \rho }{\partial i}} \right) 
    … …  
    861909\begin{align} \label{Eq_PE_slopes_eiv} 
    862910\tilde{r}_n = \begin{cases} 
    863    r_n                  &      \text{in $z-$coordinate}    \\ 
    864    r_n + \sigma_n &      \text{in $s-$coordinate}   
     911   r_n                  &      \text{in $z$-coordinate}    \\ 
     912   r_n + \sigma_n &      \text{in \textit{z*} and $s$-coordinates}   
    865913                   \end{cases} 
    866914\quad \text{where } n=1,2 
    867915\end{align} 
    868916 
    869 The normal component of the eddy induced velocity is zero at all the boundaries. this can be achieved by tapering either the eddy coefficient or the slopes to zero in the vicinity of the boundaries.  
     917The normal component of the eddy induced velocity is zero at all the boundaries. this can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}). 
    870918 
    871919\subsubsection{lateral fourth order tracer diffusive operator} 
    … …  
    877925 \end{equation} 
    878926 
    879 It is the second order operator given by (\ref{Eq_PE_iso_tensor}) applied twice with the eddy  
    880 diffusion coefficient correctly placed.  
     927It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with the eddy diffusion coefficient correctly placed.  
     928 
    881929 
    882930\subsubsection{lateral second order momentum diffusive operator} 
    883931 
    884 The second order momentum diffusive operator along $z-$ or $s-$surfaces is found by  
    885 applying \eqref{Eq_PE_curl} to the horizontal velocity vector (see Appendix B): 
     932The second order momentum diffusive operator along $z$- or $s$-surfaces is found by  
     933applying \eqref{Eq_PE_lap_vector} to the horizontal velocity vector (see Appendix~\ref{Apdx_B}): 
    886934\begin{equation} \label{Eq_PE_lapU} 
    887935\begin{split} 
    … …  
    899947 
    900948Such a formulation ensures a complete separation between the vorticity and  
    901 horizontal divergence fields ({\S}~II.4c). Unfortunately, it is not  
     949horizontal divergence fields (see Appendix~\ref{Apdx_C}). Unfortunately, it is not  
    902950available for geopotential diffusion in $s-$coordinates and for isoneutral  
    903 diffusion. In these two cases, the $u$ and $v-$fields are considered as independent  
    904 scalar fields, so that the diffusive operator is given by: 
     951diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required). In these two cases, the $u$ and $v-$fields are considered as independent scalar fields, so that the diffusive operator is given by: 
    905952\begin{equation} \label{Eq_PE_lapU_iso} 
    906953\begin{split} 
    … …  
    909956 \end{split} 
    910957 \end{equation} 
    911 where $\Re$ is given by (I.5.2). It is the same expression as those used  
    912 for diffusive operator on tracers. 
     958where $\Re$ is given by  \eqref{Eq_PE_iso_tensor}. It is the same expression as those used for diffusive operator on tracers. It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, $i.e.$ on a $f-$ or $\beta-$plane, not on the sphere. It is also a very good approximation in vicinity of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}. 
    913959 
    914960\subsubsection{lateral fourth order momentum diffusive operator} 
    915961 
    916 As for tracers, the fourth order momentum diffusive operator along $z$ or $s-$surfaces is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} with the eddy viscosity coefficient correctly placed: 
    917  
    918 geopotential diffusion in $z-$coordinates: 
     962As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} with the eddy viscosity coefficient correctly placed: 
     963 
     964geopotential diffusion in $z$-coordinate: 
    919965\begin{equation} \label{Eq_PE_bilapU} 
    920966\begin{split} 
    … …  
    928974\end{equation} 
    929975 
    930 geopotential diffusion in $s-$coordinates: 
     976\gmcomment{  change the position of the coefficient, both here and in the code} 
     977 
     978geopotential diffusion in $s$-coordinate: 
    931979\begin{equation} \label{Eq_bilapU_iso} 
    932980   \left\{   \begin{aligned} 
Note: See TracChangeset for help on using the changeset viewer.