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