Changeset 11608 for NEMO/trunk/doc
 Timestamp:
 20190927T12:24:49+02:00 (20 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
r11598 r11608 16 16 Release & Author(s) & Modifications \\ 17 17 \hline 18 {\em 4.0} & {\em ...} & {\em ...} \\ 19 {\em 3.6} & {\em ...} & {\em ...} \\ 20 {\em 3.4} & {\em ...} & {\em ...} \\ 21 {\em <=3.4} & {\em ...} & {\em ...} 18 {\em 4.0} & {\em Mike Bell } & {\em Update } \\ 19 {\em 3.6} & {\em Gurvan Madec } & {\em Minor changes} \\ 20 {\em <=3.4} & {\em Gurvan Madec and Steven Alderson} & {\em First version} \\ 22 21 \end{tabularx} 23 22 } … … 38 37 plus the following additional assumptions made from scale considerations: 39 38 40 \begin{enumerate} 41 \item \textit{spherical Earth approximation}: the geopotential surfaces are assumed to be oblate spheriods 42 that follow the Earth's bulge; these spheroids are approximated by spheres with 43 gravity locally vertical (parallel to the Earth's radius) and independent of latitude 39 \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 44 \citep[][section 2]{white.hoskins.ea_QJRMS05}. 45 \item \textit{thinshell approximation}: the ocean depth is neglected compared to the earth's radius46 \item \textit{turbulent closure hypothesis}: the turbulent fluxes45 \item [\textit{Thinshell approximation}] The ocean depth is neglected compared to the earth's radius 46 \item [\textit{Turbulent closure hypothesis}] The turbulent fluxes 47 47 (which represent the effect of small scale processes on the largescale) 48 48 are expressed in terms of largescale features 49 \item \textit{Boussinesq hypothesis}: density variations are neglected except in their contribution to50 the buoyancy force49 \item [\textit{Boussinesq hypothesis}] Density variations are neglected except in 50 their contribution to the buoyancy force 51 51 \begin{equation} 52 52 \label{eq:MB_PE_eos} 53 53 \rho = \rho \ (T,S,p) 54 54 \end{equation} 55 \item \textit{Hydrostatic hypothesis}: the vertical momentum equation is reduced to a balance between56 the vertical pressure gradient and the buoyancy force55 \item [\textit{Hydrostatic hypothesis}] The vertical momentum equation is reduced to 56 a balance between the vertical pressure gradient and the buoyancy force 57 57 (this removes convective processes from the initial NavierStokes equations and so 58 58 convective processes must be parameterized instead) … … 61 61 \pd[p]{z} =  \rho \ g 62 62 \end{equation} 63 \item \textit{Incompressibility hypothesis}: the three dimensional divergence of the velocity vector $\vect U$64 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. 65 65 \begin{equation} 66 66 \label{eq:MB_PE_continuity} 67 67 \nabla \cdot \vect U = 0 68 68 \end{equation} 69 \item \textit{Neglect of additional Coriolis terms}: the Coriolis terms that vary with the cosine of latitude are neglected. 70 These terms may be nonnegligible where the BruntVaisala frequency $N$ is small, either in the deep ocean or 71 in the submesoscale motions of the mixed layer, or near the equator \citep[][section 1]{white.hoskins.ea_QJRMS05}. 72 They can be consistently included as part of the ocean dynamics \citep[][section 3(d)]{white.hoskins.ea_QJRMS05} and are 73 retained in the MIT ocean model. 74 \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 nonnegligible where the BruntV\"{a}is\"{a}l\"{a} frequency $N$ is small, 72 either in the deep ocean or in the submesoscale 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} 75 77 76 78 Because the gravitational force is so dominant in the equations of largescale motions, … … 78 80 $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 79 81 \ie\ tangent to the geopotential surfaces. 80 Let us define the following variables: $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ 82 Let us define the following variables: 83 $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ 81 84 (the subscript $h$ denotes the local horizontal vector, \ie\ over the $(i,j)$ plane), 82 85 $T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density. … … 86 89 \label{eq:MB_PE} 87 90 \begin{gather} 88 \ intertext{$$ the momentum balance}91 \shortintertext{$$ the momentum balance} 89 92 \label{eq:MB_PE_dyn} 90 \pd[\vect U_h]{t} =  \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h 91  f \; k \times \vect U_h  \frac{1}{\rho_o} \nabla_h p 92 + \vect D^{\vect U} + \vect F^{\vect U} \\ 93 \intertext{$$ the heat and salt conservation equations} 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} 94 95 \label{eq:MB_PE_tra_T} 95 96 \pd[T]{t} =  \nabla \cdot (T \ \vect U) + D^T + F^T \\ … … 102 103 (\autoref{eq:MB_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 103 104 $f = 2 \vect \Omega \cdot k$ is the Coriolis acceleration 104 (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. 105 107 $\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of smallscale physics for momentum, 106 108 temperature and salinity, and $\vect F^{\vect U}$, $F^T$ and $F^S$ surface forcing terms. 107 Their nature and formulation are discussed in \autoref{sec:MB_zdf_ldf} and \autoref{subsec:MB_boundary_condition}. 109 Their nature and formulation are discussed in \autoref{sec:MB_zdf_ldf} and 110 \autoref{subsec:MB_boundary_condition}. 108 111 109 112 %% ================================================================================================= … … 116 119 where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface 117 120 (discretisation can introduce additional artificial ``sidewall'' boundaries). 118 Both $H$ and $\eta$ are referenced to a surface of constant geopotential (\ie\ a mean sea surface height) on which $z = 0$.119 (\ autoref{fig:MB_ocean_bc}).120 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with121 the solid earth, the continental margins, the sea ice and the atmosphere.122 However, some of these fluxes are so weak that even on climatic time scales of thousands of years123 they can be neglected.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. 124 127 In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and 125 128 the other components of the earth system. 126 129 127 \begin{figure} [!ht]130 \begin{figure} 128 131 \centering 129 132 \includegraphics[width=0.66\textwidth]{Fig_I_ocean_bc} … … 136 139 137 140 \begin{description} 138 \item [Land  ocean interface:] the major flux between continental margins and the ocean is a mass exchange of fresh 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. 139 143 Such an exchange modifies the sea surface salinity especially in the vicinity of major river mouths. 140 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 141 146 it influences the characteristics of water masses formed (especially at high latitudes). 142 147 It is required in order to close the water cycle of the climate system. 143 It is usually specified as a fresh water flux at the airsea interface in the vicinity of river mouths. 144 \item [Solid earth  ocean interface:] heat and salt fluxes through the sea floor are small, except in special areas of little extent. 145 They are usually neglected in the model 146 \footnote{ 148 It is usually specified as a fresh water flux at the airsea 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{ 147 153 In fact, it has been shown that the heat flux associated with the solid Earth cooling 148 (\ie\ the geothermal heating) is not negligible for the thermohaline circulation of the world ocean149 (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}). 150 156 }. 151 157 The boundary condition is thus set to no flux of heat and salt across solid boundaries. 152 For momentum, the situation is different. There is no flow across solid boundaries, 153 \ie\ the velocity normal to the ocean bottom and coastlines is zero (in other words, 154 the bottom velocity is parallel to solid boundaries). This kinematic boundary condition 155 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: 156 163 \begin{equation} 157 164 \label{eq:MB_w_bbc} … … 160 167 In addition, the ocean exchanges momentum with the earth through frictional processes. 161 168 Such momentum transfer occurs at small scales in a boundary layer. 162 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. 163 171 Its specification depends on the nature of the physical parameterisation used for 164 172 $\vect D^{\vect U}$ in \autoref{eq:MB_PE_dyn}. 165 It is discussed in \autoref{eq:MB_zdf}. % and Chap. III.6 to 9.166 \item [Atmosphere  ocean interface:] the kinematic surface condition plus the mass flux of fresh water PE (the precipitation minus evaporation budget)167 leads to: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: 168 176 \[ 169 177 % \label{eq:MB_w_sbc} 170 178 w = \pd[\eta]{t} + \lt. \vect U_h \rt_{z = \eta} \cdot \nabla_h (\eta) + P  E 171 179 \] 172 The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system) 173 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$. 174 183 The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. 175 \item [Sea ice  ocean interface:] 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 ($\sim46 \, 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 … … 198 208 \] 199 209 Two strategies can be considered for the surface pressure term: 200 $(a)$ introduce of a new variable $\eta$, the freesurface elevation, 210 \begin{enumerate*}[label={(\alph*)}] 211 \item introduce of a new variable $\eta$, the freesurface elevation, 201 212 for which a prognostic equation can be established and solved; 202 $(b)$assume that the ocean surface is a rigid lid,213 \item assume that the ocean surface is a rigid lid, 203 214 on which the pressure (or its horizontal gradient) can be diagnosed. 215 \end{enumerate*} 204 216 When the former strategy is used, one solution of the freesurface elevation consists of 205 217 the excitation of external gravity waves. … … 212 224 modifies certain other longwave dynamics (\eg\ barotropic Rossby or planetary waves). 213 225 The rigidlid hypothesis is an obsolescent feature in modern OGCMs. 214 It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings. 226 It has been available until the release 3.1 of \NEMO, 227 and it has been removed in release 3.2 and followings. 215 228 Only the free surface formulation is now described in this document (see the next subsection). 216 229 … … 221 234 In the free surface formulation, a variable $\eta$, the seasurface height, 222 235 is introduced which describes the shape of the airsea interface. 223 This variable is solution of a prognostic equation which is established by forming the vertical average of224 the kinematic surface condition (\autoref{eq:MB_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}): 225 238 \begin{equation} 226 239 \label{eq:MB_ssh} 227 240 \pd[\eta]{t} =  D + P  E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt] 228 241 \end{equation} 229 and using (\autoref{eq:MB_PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. 230 231 Allowing the airsea 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 airsea interface to move introduces 246 the \textbf{E}xternal \textbf{G}ravity \textbf{W}aves (EGWs) as 232 247 a class of solution of the primitive equations. 233 248 These waves are barotropic (\ie\ nearly independent of depth) and their phase speed is quite high. … … 236 251 Two choices can be made regarding the implementation of the free surface in the model, 237 252 depending on the physical processes of interest. 238 239 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction with 240 the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 241 then a non linear free surface is the most appropriate. 242 This means that no approximation is made in \autoref{eq:MB_ssh} and that 243 the variation of the ocean volume is fully taken into account. 244 Note that in order to study the fast time scales associated with EGWs it is necessary to 245 minimize time filtering effects 246 (use an explicit time scheme with very small time step, or a splitexplicit scheme with reasonably small time step, 247 see \autoref{subsec:DYN_spg_exp} or \autoref{subsec:DYN_spg_ts}). 248 249 $\bullet$ If one is not interested in EGW but rather sees them as high frequency noise, 250 it is possible to apply an explicit filter to slow down the fastest waves while 251 not altering the slow barotropic Rossby waves. 252 If further, an approximative conservation of heat and salt contents is sufficient for the problem solved, 253 then it is sufficient to solve a linearized version of \autoref{eq:MB_ssh}, 254 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{roullet.madec_JGR00}. 255 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. 256 257 The filtering of EGWs in models with a free surface is usually a matter of discretisation of 258 the temporal derivatives, 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 splitexplicit 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, 259 275 using a splitexplicit method \citep{killworth.webb.ea_JPO91, zhang.endoh_JGR92} or 260 the implicit scheme \citep{dukowicz.smith_JGR94} or 261 the addition of a filtering force in the momentum equation \citep{roullet.madec_JGR00}. 262 With the present release, \NEMO\ offers the choice between 263 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or 264 a splitexplicit scheme strongly inspired the one proposed by \citet{shchepetkin.mcwilliams_OM05} 265 (see \autoref{subsec:DYN_spg_ts}). 266 267 %% ================================================================================================= 268 \section{Curvilinear \textit{z}coordinate system} 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 splitexplicit 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} 269 284 \label{sec:MB_zco} 270 285 … … 283 298 A solution consists of introducing an appropriate coordinate transformation that 284 299 shifts the singular point onto land \citep{madec.imbard_CD96, murray_JCP96}. 285 As a consequence, it is important to solve the primitive equations in various curvilinear coordinate systems. 286 An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism. 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. 287 304 This formalism is suited to any multidimensional curvilinear coordinate system. 288 Ocean modellers mainly use threedimensional orthogonal grids on the sphere (spherical earth approximation), 289 with preservation of the local vertical. Here we give the simplified equations for this particular case. 290 The general case is detailed by \citet{eiseman.stone_SR80} in their survey of the conservation laws of fluid dynamics. 305 Ocean modellers mainly use threedimensional 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. 291 310 292 311 Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on … … 303 322 \begin{equation} 304 323 \label{eq:MB_scale_factors} 305 \begin{aligned} 306 e_1 &= (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \\ 307 e_2 &= (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \\ 308 e_3 &= \lt( \pd[z]{k} \rt) 309 \end{aligned} 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) 310 325 \end{equation} 311 326 312 \begin{figure} [!tb]327 \begin{figure} 313 328 \centering 314 \includegraphics[width=0. 66\textwidth]{Fig_I_earth_referential}329 \includegraphics[width=0.33\textwidth]{Fig_I_earth_referential} 315 330 \caption[Geographical and curvilinear coordinate systems]{ 316 331 the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear … … 328 343 \begin{subequations} 329 344 % \label{eq:MB_discrete_operators} 330 \begin{ gather}345 \begin{align} 331 346 \label{eq:MB_grad} 332 \nabla q = \frac{1}{e_1} \pd[q]{i} \; \vect i 333 + \frac{1}{e_2} \pd[q]{j} \; \vect j 334 + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ 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 \\ 335 348 \label{eq:MB_div} 336 \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] 337 + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] 338 \end{gather} 339 \begin{multline} 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] \\ 340 350 \label{eq:MB_curl} 341 \nabla \times \vect{A} = \lt[ \frac{1}{e_2} \pd[a_3]{j}  \frac{1}{e_3} \pd[a_2]{k} \rt] \vect i \\ 342 + \lt[ \frac{1}{e_3} \pd[a_1]{k}  \frac{1}{e_1} \pd[a_3]{i} \rt] \vect j \\ 343 + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i}  \pd[(e_1 a_1)]{j} \rt] \vect k 344 \end{multline} 345 \begin{gather} 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 \\ 346 352 \label{eq:MB_lap} 347 \Delta q = \nabla \cdot (\nabla q) \\353 \Delta q &= \nabla \cdot (\nabla q) \\ 348 354 \label{eq:MB_lap_vector} 349 \Delta \vect A = \nabla (\nabla \cdot \vect A)  \nabla \times (\nabla \times \vect A)350 \end{ gather}355 \Delta \vect A &= \nabla (\nabla \cdot \vect A)  \nabla \times (\nabla \times \vect A) 356 \end{align} 351 357 \end{subequations} 352 where $q$ is a scalar quantity and $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system. 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. 353 360 354 361 %% ================================================================================================= … … 357 364 358 365 In order to express the Primitive Equations in tensorial formalism, 359 it is necessary to compute the horizontal component of the nonlinear and viscous terms of the equation using 360 \autoref{eq:MB_grad}) to \autoref{eq:MB_lap_vector}. 361 Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, the velocity in the $(i,j,k)$ coordinates system and 362 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 nonlinear 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: 363 372 \begin{gather} 364 373 \label{eq:MB_curl_Uh} … … 371 380 $e_3$ is a function of the single variable $k$, 372 381 $NLT$ the nonlinear term of \autoref{eq:MB_PE_dyn} can be transformed as follows: 373 \begin{alignat*}{2} 374 &NLT &= &\lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ 375 & &= &\lt( 376 \begin{array}{*{20}c} 377 \lt[ \frac{1}{e_3} \pd[u]{k}  \frac{1}{e_1} \pd[w]{i} \rt] w  \zeta \; v \\ 378 \zeta \; u  \lt[ \frac{1}{e_2} \pd[w]{j}  \frac{1}{e_3} \pd[v]{k} \rt] \ w 379 \end{array} 380 \rt) 381 + \frac{1}{2} \lt( 382 \begin{array}{*{20}c} 383 \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\ 384 \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j} 385 \end{array} 386 \rt) \\ 387 & &= &\lt( 388 \begin{array}{*{20}c} 389 \zeta \; v \\ 390 \zeta \; u 391 \end{array} 392 \rt) 393 + \frac{1}{2} \lt( 394 \begin{array}{*{20}c} 395 \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\ 396 \frac{1}{e_2} \pd[(u^2 + v^2)]{j} 397 \end{array} 398 \rt) \\ 399 & & &+ \frac{1}{e_3} \lt( 400 \begin{array}{*{20}c} 401 w \; \pd[u]{k} \\ 402 w \; \pd[v]{k} 403 \end{array} 404 \rt) 405  \lt( 406 \begin{array}{*{20}c} 407 \frac{w}{e_1} \pd[w]{i}  \frac{1}{2 e_1} \pd[w^2]{i} \\ 408 \frac{w}{e_2} \pd[w]{j}  \frac{1}{2 e_2} \pd[w^2]{j} 409 \end{array} 410 \rt) 411 \end{alignat*} 412 The last term of the right hand side is obviously zero, and thus the nonlinear term of 413 \autoref{eq:MB_PE_dyn} is written in the $(i,j,k)$ coordinate system: 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: 414 424 \begin{equation} 415 425 \label{eq:MB_vector_form} 416 NLT = \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) 417 + \frac{1}{e_3} w \pd[\vect U_h]{k} 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} 418 427 \end{equation} 419 428 … … 421 430 For some purposes, it can be advantageous to write this term in the socalled flux form, 422 431 \ie\ to write it as the divergence of fluxes. 423 For example, the first component of \autoref{eq:MB_vector_form} (the $i$component) is transformed as follows: 424 \begin{alignat*}{2} 432 For example, 433 the first component of \autoref{eq:MB_vector_form} (the $i$component) is transformed as follows: 434 \begin{alignat*}{3} 425 435 &NLT_i &= & \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\ 426 & &= &\frac{1}{e_1 \; e_2} \lt( v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt) 427 + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) \\ 428 & & &+ \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\ 429 & &= &\frac{1}{e_1 \; e_2} \lt[  \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt) 430 + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j}  e_1 \, u \pd[v]{j} \rt) \rt. \\ 431 & & &\lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i}  u \pd[ \lt( e_2 u \rt)]{i} \rt) 432 + 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] \\ 433 438 & & &+ \frac{1}{e_3} \lt( \pd[(w \, u)]{k}  u \pd[w]{k} \rt) \\ 434 & &= &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt) 435 + \frac{1}{e_3} \pd[(w \, u)]{k} \\ 436 & & &+ \frac{1}{e_1 e_2} \lt[  u \lt( \pd[(e_1 v)]{j}  v \, \pd[e_1]{j} \rt) 437  u \pd[(e_2 u)]{i} \rt] 438  \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 \\ 439 440 & & &+ \frac{1}{e_1 e_2} \lt(  v^2 \pd[e_2]{i} \rt) \\ 440 & &= &\nabla \cdot (\vect U \, u)  (\nabla \cdot \vect U) \ u 441 + \frac{1}{e_1 e_2} \lt( v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ 442 \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:} 441 & &= &\nabla \cdot (\vect U \, u)  (\nabla \cdot \vect U) \ u + \frac{1}{e_1 e_2} \lt( v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ 442 \shortintertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:} 443 443 & &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i}  u \; \pd[e_1]{j} \rt) (v) 444 444 \end{alignat*} … … 447 447 \begin{equation} 448 448 \label{eq:MB_flux_form} 449 NLT = 450 451 452 453 454 455 449 NLT = \nabla \cdot \lt( 450 \begin{array}{*{20}c} 451 \vect U \, u \\ 452 \vect U \, v 453 \end{array} 454 \rt) 455 + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i}  u \pd[e_1]{j} \rt) \vect k \times \vect U_h 456 456 \end{equation} 457 457 458 458 The flux form has two terms, 459 the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) 460 and the second one is due to the curvilinear nature of the coordinate system used. 461 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: 462 464 \[ 463 465 % \label{eq:MB_cor+metric} … … 472 474 the following tensorial formalism: 473 475 474 \begin{ itemize}475 \item \textbf{Vector invariant form of the momentum equations}:476 \begin{description} 477 \item [Vector invariant form of the momentum equations] 476 478 \begin{equation} 477 479 \label{eq:MB_dyn_vect} 478 \begin{ split}480 \begin{gathered} 479 481 % \label{eq:MB_dyn_vect_u} 480 \pd[u]{t} = &+ (\zeta + f) \, v  \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) 481  \frac{1}{e_3} w \pd[u]{k}  \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 482 &+ 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) 484  \frac{1}{e_3} w \pd[v]{k}  \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 485 &+ D_v^{\vect U} + F_v^{\vect U} 486 \end{split} 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} 487 485 \end{equation} 488 \item \textbf{flux form of the momentum equations}:486 \item [Flux form of the momentum equations] 489 487 % \label{eq:MB_dyn_flux} 490 488 \begin{multline*} 491 489 % \label{eq:MB_dyn_flux_u} 492 \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 \\ 493  \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\ 494  \frac{1}{e_3} \pd[(w \, u)]{k}  \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 495 + D_u^{\vect U} + F_u^{\vect 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) \\ 491  \frac{1}{e_3} \pd[(w \, 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} 496 492 \end{multline*} 497 493 \begin{multline*} 498 494 % \label{eq:MB_dyn_flux_v} 499 \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 \\ 500  \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\ 501  \frac{1}{e_3} \pd[(w \, 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} 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) \\ 496  \frac{1}{e_3} \pd[(w \, 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} 503 497 \end{multline*} 504 where $\zeta$, the relative vorticity, is given by \autoref{eq:MB_curl_Uh} and $p_s$, the surface pressure,505 is given by:498 where $\zeta$, the relative vorticity, is given by \autoref{eq:MB_curl_Uh} and 499 $p_s$, the surface pressure, is given by: 506 500 \[ 507 501 % \label{eq:MB_spg} … … 518 512 \] 519 513 where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:MB_div_Uh}. 520 521 \item \textbf{tracer equations}: 522 \begin{equation} 523 \begin{split} 524 \pd[T]{t} = &  \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] 525  \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 526 \pd[S]{t} = &  \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] 527  \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\ 528 \rho = & \rho \big( T,S,z(k) \big) 529 \end{split} 530 \end{equation} 531 \end{itemize} 532 533 The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 514 \item [Tracer equations] 515 \begin{gather*} 516 \pd[T]{t} =  \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt]  \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 517 \pd[S]{t} =  \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt]  \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\ 518 \rho = \rho \big( T,S,z(k) \big) 519 \end{gather*} 520 \end{description} 521 522 The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on 523 the subgrid scale parameterisation used. 534 524 It will be defined in \autoref{eq:MB_zdf}. 535 525 The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms, … … 545 535 varying from more than 6,000 meters in abyssal trenches to zero at the coast. 546 536 Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 547 Therefore, in order to represent the ocean with respect to 548 the first pointa space and time dependent vertical coordinate that follows the variation of the sea surface height537 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 549 539 \eg\ an \zstarcoordinate; 550 for the second point, a space variation to fit the change of bottom topography 540 for the second point, 541 a space variation to fit the change of bottom topography 551 542 \eg\ a terrainfollowing or $\sigma$coordinate; 552 and for the third point, one will be tempted to use a space and time dependent coordinate that 553 follows the isopycnal surfaces, \eg\ an isopycnic coordinate. 554 555 In order to satisfy two or more constraints one can even be tempted to mixed these coordinate systems, as in 556 HYCOM (mixture of $z$coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at 557 the ocean bottom) \citep{chassignet.smith.ea_JPO03} or 558 OPA (mixture of $z$coordinate in vicinity the surface and steep topography areas and $\sigma$coordinate elsewhere) 559 \citep{madec.delecluse.ea_JPO96} among others. 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. 560 552 561 553 In fact one is totally free to choose any space and time vertical coordinate by … … 565 557 s = s(i,j,k,t) 566 558 \end{equation} 567 with the restriction that the above equation gives a singlevalued monotonic relationship between $s$ and $k$, 568 when $i$, $j$ and $t$ are held fixed. 569 \autoref{eq:MB_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 singlevalued 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 570 563 the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through 571 564 \autoref{eq:MB_s}. 572 565 This socalled \textit{generalised vertical coordinate} \citep{kasahara_MWR74} is in fact 573 an Arbitrary LagrangianEulerian (ALE) coordinate. 574 Indeed, one has a great deal of freedom in the choice of expression for $s$. The choice determines 575 which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and 576 which part will be used to move them (Lagrangian part). 577 The coordinate is also sometime referenced as an adaptive coordinate \citep{hofmeister.burchard.ea_OM10}, 578 since the coordinate system is adapted in the course of the simulation. 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. 579 573 Its most often used implementation is via an ALE algorithm, 580 574 in which a pure lagrangian step is followed by regridding and remapping steps, … … 593 587 The horizontal velocity in ocean models measures motions in the horizontal plane, 594 588 perpendicular to the local gravitational field. 595 That is, horizontal velocity is mathematically the same regardless of the vertical coordinate, be it geopotential,596 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. 597 591 The key motivation for maintaining the same horizontal velocity component is that 598 592 the hydrostatic and geostrophic balances are dominant in the largescale ocean. 599 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, 600 595 would lead to unacceptable numerical errors. 601 596 Correspondingly, the vertical direction is anti parallel to the gravitational force in … … 604 599 the surface of a constant generalized vertical coordinate. 605 600 606 It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between 607 the vertical coordinate choices. 608 That is, computation of the diasurface 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 diasurface velocity component represents the fundamental distinction between 609 605 the various coordinates. 610 In some models, such as geopotential, pressure, and terrain following, this transport is typically diagnosed from 611 volume or mass conservation. 612 In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about 613 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. 614 611 615 612 In this section we first establish the PE in the generalised vertical $s$coordinate, … … 620 617 \subsection{\textit{S}coordinate formulation} 621 618 622 Starting from the set of equations established in \autoref{sec:MB_zco} for the special case $k = z$ and 623 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)$, 624 622 which includes $z$, \zstar and $\sigma$coordinates as special cases 625 623 ($s = z$, $s = \zstar$, and $s = \sigma = z / H$ or $ = z / \lt( H + \eta \rt)$, resp.). 626 624 A formal derivation of the transformed equations is given in \autoref{apdx:SCOORD}. 627 Let us define the vertical scale factor by $e_3 = \partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ), 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)$ ), 628 627 and the slopes in the $(i,j)$ directions between $s$ and $z$surfaces by: 629 628 \begin{equation} … … 632 631 \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt_s 633 632 \end{equation} 634 We also introduce $\omega$, a diasurface velocity component, defined as the velocity635 relative to the moving $s$surfaces and normal to them:633 We also introduce $\omega$, a diasurface velocity component, 634 defined as the velocity relative to the moving $s$surfaces and normal to them: 636 635 \[ 637 636 % \label{eq:MB_sco_w} … … 642 641 (see \autoref{sec:SCOORD_momentum}): 643 642 644 \begin{itemize} 645 \item \textbf{Vector invariant form of the momentum equation}: 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] 646 652 \begin{multline*} 647 % \label{eq:MB_sco_u_vector} 648 \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} \\ 649  \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt)  g \frac{\rho}{\rho_o} \sigma_1 650 + D_u^{\vect U} + F_u^{\vect U} 651 \end{multline*} 652 \begin{multline*} 653 % \label{eq:MB_sco_v_vector} 654 \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} \\ 655  \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt)  g \frac{\rho}{\rho_o} \sigma_2 656 + D_v^{\vect U} + F_v^{\vect U} 657 \end{multline*} 658 \item \textbf{Flux form of the momentum equation}: 659 \begin{multline*} 660 % \label{eq:MB_sco_u_flux} 661 \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 \\ 662  \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) \\ 663  \frac{1}{e_3} \pd[(\omega \, u)]{k} 664  \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 665  g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 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) \\ 655  \frac{1}{e_3} \pd[(\omega \, 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} 666 656 \end{multline*} 667 657 \begin{multline*} 668 658 % \label{eq:MB_sco_v_flux} 669 \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 \\ 670  \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) \\ 671  \frac{1}{e_3} \pd[(\omega \, v)]{k} 672  \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 673  g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 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) \\ 660  \frac{1}{e_3} \pd[(\omega \, 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} 674 661 \end{multline*} 675 662 where the relative vorticity, $\zeta$, the surface pressure gradient, … … 682 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) 683 670 \] 684 \item \textit{tracer equations}: 685 \begin{multline*} 686 % \label{eq:MB_sco_t} 687 \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} 688 + \pd[(e_1 e_3 \, v \, T)]{j} \rt) \\ 689  \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S 690 \end{multline*} 691 \begin{multline} 692 % \label{eq:MB_sco_s} 693 \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} 694 + \pd[(e_1 e_3 \, v \, S)]{j} \rt) \\ 695  \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S 696 \end{multline} 697 \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} 698 679 The equation of state has the same expression as in $z$coordinate, 699 680 and similar expressions are used for mixing and forcing terms. … … 709 690 \label{subsec:MB_zco_star} 710 691 711 \begin{figure} [!b]692 \begin{figure} 712 693 \centering 713 694 \includegraphics[width=0.66\textwidth]{Fig_z_zstar} 714 695 \caption[Curvilinear zcoordinate systems (\{non\}linear freesurface cases and rescaled \zstar)]{ 715 (a) $z$coordinate in linear freesurface case ; 716 (b) $z$coordinate in nonlinear free surface case ; 717 (c) rescaled height coordinate 718 (become popular as the \zstarcoordinate \citep{adcroft.campin_OM04}).} 696 \begin{enumerate*}[label={(\alph*)}] 697 \item $z$coordinate in linear freesurface case; 698 \item $z$coordinate in nonlinear free surface case; 699 \item rescaled height coordinate 700 (become popular as the \zstarcoordinate \citep{adcroft.campin_OM04}). 701 \end{enumerate*} 702 } 719 703 \label{fig:MB_z_zstar} 720 704 \end{figure} 721 705 722 In this case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. 723 These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on the \NEMO\ web site. 724 725 The \zstar coordinate approach is an unapproximated, nonlinear free surface implementation which allows one to 726 deal with large amplitude freesurface variations relative to the vertical resolution \citep{adcroft.campin_OM04}. 727 In the \zstar formulation, 728 the variation of the column thickness due to seasurface undulations is not concentrated in the surface level, 729 as in the $z$coordinate formulation, but is equally distributed over the full water column. 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, nonlinear free surface implementation which 712 allows one to deal with large amplitude freesurface variations relative to the vertical resolution 713 \citep{adcroft.campin_OM04}. 714 In the \zstar\ formulation, the variation of the column thickness due to seasurface 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. 730 717 Thus vertical levels naturally follow seasurface variations, with a linear attenuation with depth, 731 718 as illustrated by \autoref{fig:MB_z_zstar}. 732 Note that with a flat bottom, such as in \autoref{fig:MB_z_zstar}, the bottomfollowing $z$ coordinate and \zstar are equivalent. 719 Note that with a flat bottom, such as in \autoref{fig:MB_z_zstar}, 720 the bottomfollowing $z$ coordinate and \zstar\ are equivalent. 733 721 The definition and modified oceanic equations for the rescaled vertical coordinate \zstar, 734 722 including the treatment of freshwater flux at the surface, are detailed in Adcroft and Campin (2004). 735 723 The major points are summarized here. 736 The position (\zstar) and vertical discretization ( \zstar) are expressed as:724 The position (\zstar) and vertical discretization ($\delta \zstar$) are expressed as: 737 725 \[ 738 726 % \label{eq:MB_zstar} 739 H + \zstar = (H + z) / r \quad \text{and} \quad \delta \zstar 740 = \delta z / r \quad \text{with} \quad r 741 = \frac{H + \eta}{H} . 727 H + \zstar = (H + z) / r \quad \text{and} \quad \delta \zstar = \delta z / r \quad \text{with} \quad r = \frac{H + \eta}{H} 742 728 \] 743 729 Simple reorganisation of the above expressions gives 744 730 \[ 745 731 % \label{eq:MB_zstar_2} 746 \zstar = H \lt( \frac{z  \eta}{H + \eta} \rt) .732 \zstar = H \lt( \frac{z  \eta}{H + \eta} \rt) 747 733 \] 748 734 Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar, 749 the upper and lower boundaries are at fixed \zstarposition,735 the upper and lower boundaries are at fixed \zstar\ position, 750 736 $\zstar = 0$ and $\zstar = H$ respectively. 751 737 Also the divergence of the flow field is no longer zero as shown by the continuity equation: 752 738 \[ 753 \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 .739 \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 754 740 \] 755 This \zstar coordinate is closely related to the "eta"coordinate used in many atmospheric models756 (see Black (1994) for a review of etacoordinate 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). 757 743 It was originally used in ocean models by Stacey et al. (1995) for studies of tides next to shelves, 758 744 and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling. 759 745 760 The surfaces of constant \zstar are quasihorizontal.761 Indeed, the \zstar coordinate reduces to $z$ when $\eta$ is zero.746 The surfaces of constant \zstar\ are quasihorizontal. 747 Indeed, the \zstar\ coordinate reduces to $z$ when $\eta$ is zero. 762 748 In general, when noting the large differences between 763 749 undulations of the bottom topography versus undulations in the surface height, 764 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. 765 751 These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to 766 752 terrain following sigma models discussed in \autoref{subsec:MB_sco}. 767 753 Additionally, since $\zstar = z$ when $\eta = 0$, 768 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography. 769 This behaviour is in contrast to the case with "s"models, where pressure gradient errors in the presence of 770 nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, 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, 771 759 depending on the sophistication of the pressure gradient solver. 772 The quasi 773 neutral physics parameterizations in \zstar 760 The quasihorizontal 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 774 762 (see Chapters 1316 of \cite{griffies_bk04}) for a discussion of neutral physics in $z$models, 775 763 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). 776 764 777 The range over which \zstar 765 The range over which \zstar\ varies is time independent $H \leq \zstar \leq 0$. 778 766 Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > H$. 779 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z  \eta$. 780 781 Because \zstar has a time independent range, all grid cells have static increments ds, 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, 782 771 and the sum of the vertical increments yields the time independent ocean depth. %k ds = H. 783 The \zstar coordinate is therefore invisible to undulations of the free surface,772 The \zstar\ coordinate is therefore invisible to undulations of the free surface, 784 773 since it moves along with the free surface. 785 This property means that no spurious vertical transport is induced across surfaces of constant \zstar by786 the motion of external gravity waves.774 This property means that no spurious vertical transport is induced across 775 surfaces of constant \zstar\ by the motion of external gravity waves. 787 776 Such spurious transport can be a problem in zmodels, especially those with tidal forcing. 788 Quite generally, the time independent range for the \zstar coordinate is a very convenient property that 789 allows for a nearly arbitrary vertical resolution even in the presence of large amplitude fluctuations of 790 the surface height, again so long as $\eta > H$. 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$. 791 781 %end MOM doc %%% 792 782 … … 795 785 \label{subsec:MB_sco} 796 786 797 %% =================================================================================================798 \subsubsection{Introduction}799 800 787 Several important aspects of the ocean circulation are influenced by bottom topography. 801 Of course, the most important is that bottom topography determines deep ocean subbasins, barriers, sills and 802 channels that strongly constrain the path of water masses, but more subtle effects exist. 803 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 subbasins, 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. 804 792 Topographic Rossby waves can be excited and can interact with the mean current. 805 793 In the $z$coordinate system presented in the previous section (\autoref{sec:MB_zco}), … … 809 797 large localized depth gradients associated with large localized vertical velocities. 810 798 The response to such a velocity field often leads to numerical dispersion effects. 811 One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of812 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}. 813 801 Another solution is to introduce a terrainfollowing coordinate system (hereafter $s$coordinate). 814 802 815 The $s$coordinate avoids the discretisation error in the depth field since the layers of 816 computation are gradually adjusted with depth to the ocean bottom. 817 Relatively small topographic features as well as gentle, largescale slopes of the sea floor in the deep ocean, 818 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, largescale 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, 819 809 can easily be represented (with relatively low vertical resolution). 820 A terrainfollowing model (hereafter $s$model) also facilitates the modelling of the boundary layer flows over 821 a large depth range, which in the framework of the $z$model would require high vertical resolution over 810 A terrainfollowing 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 822 813 the whole depth range. 823 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 824 816 the only boundaries of the domain (no more lateral boundary condition to specify). 825 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, 826 819 it has strong limitations as soon as stratification is introduced. 827 820 The main two problems come from the truncation error in the horizontal pressure gradient and 828 821 a possibly increased diapycnal diffusion. 829 822 The horizontal pressure force in $s$coordinate consists of two terms (see \autoref{apdx:SCOORD}), 830 831 823 \begin{equation} 832 824 \label{eq:MB_p_sco} … … 836 828 The second term in \autoref{eq:MB_p_sco} depends on the tilt of the coordinate surface and 837 829 leads to a truncation error that is not present in a $z$model. 838 In the special case of a $\sigma$coordinate (i.e. a depthnormalised coordinate system $\sigma = z/H$), 839 \citet{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of the magnitude of this truncation error. 840 It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, 841 and the finite difference scheme. 830 In the special case of a $\sigma$coordinate 831 (i.e. a depthnormalised 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. 842 836 This error limits the possible topographic slopes that a model can handle at 843 837 a given horizontal and vertical resolution. 844 838 This is a severe restriction for largescale applications using realistic bottom topography. 845 The largescale slopes require high horizontal resolution, and the computational cost becomes prohibitive. 839 The largescale slopes require high horizontal resolution, 840 and the computational cost becomes prohibitive. 846 841 This problem can be at least partially overcome by mixing $s$coordinate and 847 steplike representation of bottom topography \citep{gerdes_JGR93*a,gerdes_JGR93*b,madec.delecluse.ea_JPO96}. 842 steplike representation of bottom topography 843 \citep{gerdes_JGR93*a,gerdes_JGR93*b,madec.delecluse.ea_JPO96}. 848 844 However, the definition of the model domain vertical coordinate becomes then a nontrivial thing for 849 845 a realistic bottom topography: 850 an envelope topography is defined in $s$coordinate on which a full or851 partial step bottom topography is then applied in order to adjust the model depth to the observed one 852 (see \autoref{subsec:DOM_zgr}.853 854 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of855 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. 856 852 It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. 857 853 This is the case for a $z$model as well as for a $s$model. 858 However, density varies more strongly on $s$surfaces than on horizontal surfaces in regions of 859 large topographic slopes, implying larger diapycnal diffusion in a $s$model than in a $z$model. 860 Whereas such a diapycnal diffusion in a $z$model tends to weaken horizontal density (pressure) gradients and thus 861 the horizontal circulation, it usually reinforces these gradients in a $s$model, creating spurious circulation. 862 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. 863 862 Spurious diffusion along $s$surfaces will induce a bump of isoneutral surfaces over the topography, 864 863 and thus will generate there a baroclinic eddy. 865 864 In contrast, the ocean will stay at rest in a $z$model. 866 As for the truncation error, the problem can be reduced by introducing the terrainfollowing coordinate below 867 the strongly stratified portion of the water column (\ie\ the main thermocline) \citep{madec.delecluse.ea_JPO96}. 868 An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces 869 (see \autoref{subsec:MB_ldf}). 865 As for the truncation error, the problem can be reduced by 866 introducing the terrainfollowing 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}). 870 870 Unfortunately, the slope of isoneutral surfaces relative to the $s$surfaces can very large, 871 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 872 873 The $s$coordinates introduced here \citep{lott.madec.ea_OM90,madec.delecluse.ea_JPO96} differ mainly in two aspects from 874 similar models: 875 it allows a representation of bottom topography with mixed full or partial steplike/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 steplike/terrain following topography; 876 878 It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. 877 879 878 880 %% ================================================================================================= 879 \subsection{ \texorpdfstring{Curvilinear \ztildecoordinate}{}}881 \subsection{Curvilinear \ztildecoordinate} 880 882 \label{subsec:MB_zco_tilde} 881 883 882 The \ztilde 884 The \ztildecoordinate has been developed by \citet{leclair.madec_OM11}. 883 885 It is available in \NEMO\ since the version 3.4 and is more robust in version 4.0 than previously. 884 886 Nevertheless, it is currently not robust enough to be used in all possible configurations. … … 889 891 \label{sec:MB_zdf_ldf} 890 892 891 The hydrostatic primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than 892 a few kilometres in the horizontal, a few meters in the vertical and a few minutes. 893 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. 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. 894 898 The effects of smaller scale motions (coming from the advective terms in the NavierStokes equations) 895 899 must be represented entirely in terms of largescale patterns to close the equations. … … 902 906 small scale processes \textit{in fine} balance the surface input of kinetic energy and heat. 903 907 904 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. 905 910 Therefore subgridscale physics \textbf{D}$^{\vect U}$, $D^{S}$ and $D^{T}$ in 906 911 \autoref{eq:MB_PE_dyn}, \autoref{eq:MB_PE_tra_T} and \autoref{eq:MB_PE_tra_S} are divided into 907 a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and912 a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and 908 913 a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$. 909 The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 914 The formulation of these terms and their underlying physics are briefly discussed in 915 the next two subsections. 910 916 911 917 %% ================================================================================================= … … 913 919 \label{subsec:MB_zdf} 914 920 915 The model resolution is always larger than the scale at which the major sources of vertical turbulence occur916 (shear instability, internal wave breaking...).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...). 917 923 Turbulent motions are thus never explicitly solved, even partially, but always parameterized. 918 The vertical turbulent fluxes are assumed to depend linearly on the gradients of largescale quantities 919 (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 largescale quantities (for example, 926 the turbulent heat flux is given by $\overline{T' w'} = A^{v T} \partial_z \overline T$, 920 927 where $A^{v T}$ is an eddy coefficient). 921 928 This formulation is analogous to that of molecular diffusion and dissipation. … … 927 934 \label{eq:MB_zdf} 928 935 \begin{gathered} 929 \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \ , \\930 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} \ 931 938 D^{vS} = \pd[]{z} \lt( A^{vT} \pd[S]{z} \rt) 932 939 \end{gathered} 933 940 \end{equation} 934 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. 935 943 At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must be specified 936 944 (see \autoref{chap:SBC} and \autoref{chap:ZDF} and \autoref{sec:TRA_bbl}). 937 945 All the vertical physics is embedded in the specification of the eddy coefficients. 938 946 They can be assumed to be either constant, or function of the local fluid properties 939 (\eg\ Richardson number, BruntV ais\"{a}l\"{a} frequency, distance from the boundary ...),947 (\eg\ Richardson number, BruntV\"{a}is\"{a}l\"{a} frequency, distance from the boundary ...), 940 948 or computed from a turbulent closure model. 941 949 The choices available in \NEMO\ are discussed in \autoref{chap:ZDF}). … … 948 956 (which can be solved explicitly if the resolution is sufficient since 949 957 their underlying physics are included in the primitive equations), 950 and a sub mesoscale turbulence which is never explicitly solved even partially, but always parameterized. 951 The formulation of lateral eddy fluxes depends on whether the mesoscale is below or above the gridspacing 952 (\ie\ the model is eddyresolving 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 gridspacing (\ie\ the model is eddyresolving or not). 953 962 954 963 In noneddyresolving configurations, the closure is similar to that used for the vertical physics. 955 The lateral turbulent fluxes are assumed to depend linearly on the lateral gradients of largescale quantities. 964 The lateral turbulent fluxes are assumed to depend linearly on 965 the lateral gradients of largescale quantities. 956 966 The resulting lateral diffusive and dissipative operators are of second order. 957 Observations show that lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces 967 Observations show that 968 lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces 958 969 (or more precisely neutral surfaces \cite{mcdougall_JPO87}) rather than across them. 959 As the slope of neutral surfaces is small in the ocean, a common approximation is to assume that 960 the `lateral' direction is the horizontal, \ie\ the lateral mixing is performed along geopotential surfaces. 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. 961 973 This leads to a geopotential second order operator for lateral subgrid scale physics. 962 This assumption can be relaxed: the eddyinduced turbulent fluxes can be better approached by assuming that 974 This assumption can be relaxed: 975 the eddyinduced turbulent fluxes can be better approached by assuming that 963 976 they depend linearly on the gradients of largescale quantities computed along neutral surfaces. 964 977 In such a case, the diffusive operator is an isoneutral second order operator and 965 978 it has components in the three space directions. 966 However, 967 both horizontal and isoneutral operators have no effect onmean (\ie\ large scale) potential energy whereas979 However, both horizontal and isoneutral operators have no effect on 980 mean (\ie\ large scale) potential energy whereas 968 981 potential energy is a main source of turbulence (through baroclinic instabilities). 969 982 \citet{gent.mcwilliams_JPO90} proposed a parameterisation of mesoscale eddyinduced turbulence which 970 983 associates an eddyinduced velocity to the isoneutral diffusion. 971 984 Its mean effect is to reduce the mean potential energy of the ocean. 972 This leads to a formulation of lateral subgridscale physics made up of an isoneutral second order operator and973 an eddy induced advective part.985 This leads to a formulation of lateral subgridscale physics made up of 986 an isoneutral second order operator and an eddy induced advective part. 974 987 In all these lateral diffusive formulations, 975 988 the specification of the lateral eddy coefficients remains the problematic point as 976 there is no really satisfactory formulation of these coefficients as a function of largescale features. 989 there is no really satisfactory formulation of these coefficients as 990 a function of largescale features. 977 991 978 992 In eddyresolving configurations, a second order operator can be used, … … 983 997 not interfering with the resolved mesoscale activity. 984 998 Another approach is becoming more and more popular: 985 instead of specifying explicitly a subgrid scale term in the momentum and tracer time evolution equations, 999 instead of specifying explicitly a subgrid scale term in 1000 the momentum and tracer time evolution equations, 986 1001 one uses an advective scheme which is diffusive enough to maintain the model stability. 987 It must be emphasised that then, all the subgrid scale physics is included in the formulation of988 the advection scheme.1002 It must be emphasised that then, 1003 all the subgrid scale physics is included in the formulation of the advection scheme. 989 1004 990 1005 All these parameterisations of subgrid scale physics have advantages and drawbacks. 991 They 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: 992 1008 Laplacian and bilaplacian operators acting along geopotential or isoneutral surfaces, 993 1009 \citet{gent.mcwilliams_JPO90} parameterisation, and various slightly diffusive advection schemes. 994 For momentum, the main ones are: Laplacian and bilaplacian operators acting along geopotential surfaces, 1010 For momentum, the main ones are: 1011 Laplacian and bilaplacian operators acting along geopotential surfaces, 995 1012 and UBS advection schemes when flux form is chosen for the momentum advection. 996 1013 … … 1001 1018 \begin{equation} 1002 1019 \label{eq:MB_iso_tensor} 1003 D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad 1004 \Re = 1005 \begin{pmatrix} 1006 1 & 0 & r_1 \\ 1007 0 & 1 & r_2 \\ 1008 r_1 & r_2 & r_1^2 + r_2^2 \\ 1009 \end{pmatrix} 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} 1010 1026 \end{equation} 1011 1027 where $r_1$ and $r_2$ are the slopes between the surface along which the diffusive operator acts and … … 1014 1030 the rotation between geopotential and $s$surfaces, 1015 1031 while it is only an approximation for the rotation between isoneutral and $z$ or $s$surfaces. 1016 Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:MB_iso_tensor} \citep{cox_OM87}. 1017 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and 1018 dianeutral 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 dianeutral diffusive coefficients is known to be 1036 several orders of magnitude smaller than unity. 1019 1037 Second, the two isoneutral directions of diffusion are assumed to be independent since 1020 1038 the slopes are generally less than $10^{2}$ in the ocean (see \autoref{apdx:DIFFOPERS}). … … 1027 1045 they are equal to $\sigma_1$ and $\sigma_2$, respectively (see \autoref{eq:MB_sco_slope}). 1028 1046 1029 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral and computational surfaces. 1047 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between 1048 the isoneutral and computational surfaces. 1030 1049 Therefore, they are different quantities, but have similar expressions in $z$ and $s$coordinates. 1031 1050 In $z$coordinates: … … 1048 1067 where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a nondivergent, 1049 1068 eddyinduced transport velocity. This velocity field is defined by: 1050 \begin{gather }1069 \begin{gather*} 1051 1070 % \label{eq:MB_eiv} 1052 u^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_1 \rt) \ \1071 u^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_1 \rt) \quad 1053 1072 v^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_2 \rt) \\ 1054 1073 w^\ast =  \frac{1}{e_1 e_2} \lt[ \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt) 1055 1074 + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt] 1056 \end{gather }1075 \end{gather*} 1057 1076 where $A^{eiv}$ is the eddy induced velocity coefficient 1058 1077 (or equivalently the isoneutral thickness diffusivity coefficient), 1059 and $\tilde r_1$ and $\tilde r_2$ are the slopes between isoneutral and \textit{geopotential} surfaces. 1060 Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate: 1061 \begin{align} 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} 1062 1083 \label{eq:MB_slopes_eiv} 1063 1084 \tilde{r}_n = … … 1067 1088 \end{cases} 1068 1089 \quad \text{where~} n = 1, 2 1069 \end{ align}1090 \end{equation} 1070 1091 1071 1092 The normal component of the eddy induced velocity is zero at all the boundaries. 1072 This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of1073 the boundaries.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. 1074 1095 The latter strategy is used in \NEMO\ (cf. \autoref{chap:LDF}). 1075 1096 … … 1101 1122 \end{align*} 1102 1123 1103 Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields1104 (see \autoref{apdx:INVARIANTS}).1124 Such a formulation ensures a complete separation between 1125 the vorticity and horizontal divergence fields (see \autoref{apdx:INVARIANTS}). 1105 1126 Unfortunately, it is only available in \textit{isolevel} direction. 1106 1127 When a rotation is required 1107 (\ie\ geopotential diffusion in $s$coordinates or isoneutral diffusion in both $z$ and $s$coordinates), 1108 the $u$ and $v$fields are considered as independent scalar fields, so that the diffusive operator is given by: 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: 1109 1132 \begin{gather*} 1110 1133 % \label{eq:MB_lapU_iso} … … 1122 1145 \subsubsection{Lateral bilaplacian momentum diffusive operator} 1123 1146 1124 As for tracers, the bilaplacian order momentum diffusive operator is a reentering Laplacian operator with 1147 As for tracers, 1148 the bilaplacian order momentum diffusive operator is a reentering Laplacian operator with 1125 1149 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1126 1150 Nevertheless it is currently not available in the isoneutral case.
Note: See TracChangeset
for help on using the changeset viewer.