Changeset 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/doc/latex/NEMO/subfiles/chap_model_basics.tex
- Timestamp:
- 2019-09-19T11:18:03+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/doc/latex/NEMO/subfiles/chap_model_basics.tex
r11151 r11573 1 1 2 \documentclass[../main/NEMO_manual]{subfiles} 2 3 … … 7 8 % ================================================================ 8 9 \chapter{Model Basics} 9 \label{chap:PE} 10 \minitoc 10 \label{chap:MB} 11 12 \chaptertoc 11 13 12 14 \newpage … … 16 18 % ================================================================ 17 19 \section{Primitive equations} 18 \label{sec: PE_PE}19 20 % ------------------------------------------------------------------------------------------------------------- 21 % Vector Invariant Formulation 20 \label{sec:MB_PE} 21 22 % ------------------------------------------------------------------------------------------------------------- 23 % Vector Invariant Formulation 22 24 % ------------------------------------------------------------------------------------------------------------- 23 25 24 26 \subsection{Vector invariant formulation} 25 \label{subsec: PE_Vector}27 \label{subsec:MB_PE_vector} 26 28 27 29 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 which30 \ie\ the Navier-Stokes equations along with a nonlinear equation of state which 29 31 couples the two active tracers (temperature and salinity) to the fluid velocity, 30 32 plus the following additional assumptions made from scale considerations: … … 32 34 \begin{enumerate} 33 35 \item 34 \textit{spherical earth approximation}: the geopotential surfaces are assumed to be spheres so that 35 gravity (local vertical) is parallel to the earth's radius 36 \textit{spherical Earth approximation}: the geopotential surfaces are assumed to be oblate spheriods 37 that follow the Earth's bulge; these spheroids are approximated by spheres with 38 gravity locally vertical (parallel to the Earth's radius) and independent of latitude 39 \citep[][section 2]{white.hoskins.ea_QJRMS05}. 36 40 \item 37 41 \textit{thin-shell approximation}: the ocean depth is neglected compared to the earth's radius … … 44 48 the buoyancy force 45 49 \begin{equation} 46 \label{eq: PE_eos}50 \label{eq:MB_PE_eos} 47 51 \rho = \rho \ (T,S,p) 48 52 \end{equation} … … 53 57 convective processes must be parameterized instead) 54 58 \begin{equation} 55 \label{eq: PE_hydrostatic}59 \label{eq:MB_PE_hydrostatic} 56 60 \pd[p]{z} = - \rho \ g 57 61 \end{equation} … … 60 64 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} 69 \item 70 \textit{Neglect of additional Coriolis terms}: the Coriolis terms that vary with the cosine of latitude are neglected. 71 These terms may be non-negligible where the Brunt-Vaisala frequency $N$ is small, either in the deep ocean or 72 in the sub-mesoscale motions of the mixed layer, or near the equator \citep[][section 1]{white.hoskins.ea_QJRMS05}. 73 They can be consistently included as part of the ocean dynamics \citep[][section 3(d)]{white.hoskins.ea_QJRMS05} and are 74 retained in the MIT ocean model. 65 75 \end{enumerate} 66 76 67 77 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 that78 it is useful to choose an orthogonal set of unit vectors $(i,j,k)$ linked to the Earth such that 69 79 $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 70 \ie tangent to the geopotential surfaces.80 \ie\ tangent to the geopotential surfaces. 71 81 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),82 (the subscript $h$ denotes the local horizontal vector, \ie\ over the $(i,j)$ plane), 73 83 $T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density. 74 84 The vector invariant form of the primitive equations in the $(i,j,k)$ vector system provides 75 85 the following equations: 76 86 \begin{subequations} 77 \label{eq: PE}87 \label{eq:MB_PE} 78 88 \begin{gather} 79 89 \intertext{$-$ the momentum balance} 80 \label{eq: PE_dyn}90 \label{eq:MB_PE_dyn} 81 91 \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 92 - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p 83 93 + \vect D^{\vect U} + \vect F^{\vect U} \\ 84 94 \intertext{$-$ the heat and salt conservation equations} 85 \label{eq: PE_tra_T}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 105 (where $\vect \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 96 106 $\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, 97 107 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}.108 Their nature and formulation are discussed in \autoref{sec:MB_zdf_ldf} and \autoref{subsec:MB_boundary_condition}. 99 109 100 110 % ------------------------------------------------------------------------------------------------------------- … … 102 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}). 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 (\ie\ a mean sea surface height) on which $z = 0$. 122 (\autoref{fig:MB_ocean_bc}). 112 123 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with 113 124 the solid earth, the continental margins, the sea ice and the atmosphere. … … 119 130 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 120 131 \begin{figure}[!ht] 121 \begin{center} 122 \includegraphics[width=\textwidth]{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} 132 \centering 133 \includegraphics[width=0.66\textwidth]{Fig_I_ocean_bc} 134 \caption[Ocean boundary conditions]{ 135 The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, 136 where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. 137 Both $H$ and $\eta$ are referenced to $z = 0$.} 138 \label{fig:MB_ocean_bc} 130 139 \end{figure} 131 140 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 144 153 \footnote{ 145 154 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 ocean155 (\ie\ the geothermal heating) is not negligible for the thermohaline circulation of the world ocean 147 156 (see \autoref{subsec:TRA_bbc}). 148 157 }. 149 158 The boundary condition is thus set to no flux of heat and salt across solid boundaries. 150 159 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,160 \ie\ the velocity normal to the ocean bottom and coastlines is zero (in other words, 152 161 the bottom velocity is parallel to solid boundaries). This kinematic boundary condition 153 162 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} … … 160 169 It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. 161 170 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.171 $\vect D^{\vect U}$ in \autoref{eq:MB_PE_dyn}. 172 It is discussed in \autoref{eq:MB_zdf}.% and Chap. III.6 to 9. 164 173 \item[Atmosphere - ocean interface:] 165 174 the kinematic surface condition plus the mass flux of fresh water PE (the precipitation minus evaporation budget) 166 175 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 \] … … 183 192 % ================================================================ 184 193 \section{Horizontal pressure gradient} 185 \label{sec: PE_hor_pg}194 \label{sec:MB_hor_pg} 186 195 187 196 % ------------------------------------------------------------------------------------------------------------- … … 189 198 % ------------------------------------------------------------------------------------------------------------- 190 199 \subsection{Pressure formulation} 191 \label{subsec: PE_p_formulation}200 \label{subsec:MB_p_formulation} 192 201 193 202 The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at 194 203 a reference geopotential surface ($z = 0$) and a hydrostatic pressure $p_h$ such that: 195 204 $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}).205 The latter is computed by integrating (\autoref{eq:MB_PE_hydrostatic}), 206 assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:MB_PE_eos}). 198 207 The hydrostatic pressure is then given by: 199 208 \[ 200 % \label{eq: PE_pressure}209 % \label{eq:MB_pressure} 201 210 p_h (i,j,z,t) = \int_{\varsigma = z}^{\varsigma = 0} g \; \rho (T,S,\varsigma) \; d \varsigma 202 211 \] … … 210 219 The flow is barotropic and the surface moves up and down with gravity as the restoring force. 211 220 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.221 the time step has to be very short when they are present in the model. 213 222 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$.223 \ie\ the sea surface is the surface $z = 0$. 215 224 This well known approximation increases the surface wave speed to infinity and 216 modifies certain other longwave dynamics (\eg barotropic Rossby or planetary waves).225 modifies certain other longwave dynamics (\eg\ barotropic Rossby or planetary waves). 217 226 The rigid-lid hypothesis is an obsolescent feature in modern OGCMs. 218 227 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 th e this document (see the next sub-section).228 Only the free surface formulation is now described in this document (see the next sub-section). 220 229 221 230 % ------------------------------------------------------------------------------------------------------------- … … 223 232 % ------------------------------------------------------------------------------------------------------------- 224 233 \subsection{Free surface formulation} 225 \label{subsec: PE_free_surface}234 \label{subsec:MB_free_surface} 226 235 227 236 In the free surface formulation, a variable $\eta$, the sea-surface height, 228 237 is introduced which describes the shape of the air-sea interface. 229 238 This variable is solution of a prognostic equation which is established by forming the vertical average of 230 the kinematic surface condition (\autoref{eq: PE_w_bbc}):239 the kinematic surface condition (\autoref{eq:MB_w_bbc}): 231 240 \begin{equation} 232 \label{eq: PE_ssh}241 \label{eq:MB_ssh} 233 242 \pd[\eta]{t} = - D + P - E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt] 234 243 \end{equation} 235 and using (\autoref{eq: PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$.244 and using (\autoref{eq:MB_PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. 236 245 237 246 Allowing the air-sea interface to move introduces the external gravity waves (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 … … 246 255 the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 247 256 then a non linear free surface is the most appropriate. 248 This means that no approximation is made in \autoref{eq: PE_ssh} and that257 This means that no approximation is made in \autoref{eq:MB_ssh} and that 249 258 the variation of the ocean volume is fully taken into account. 250 259 Note that in order to study the fast time scales associated with EGWs it is necessary to … … 257 266 not altering the slow barotropic Rossby waves. 258 267 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},268 then it is sufficient to solve a linearized version of \autoref{eq:MB_ssh}, 260 269 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{roullet.madec_JGR00}. 261 270 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. … … 266 275 the implicit scheme \citep{dukowicz.smith_JGR94} or 267 276 the addition of a filtering force in the momentum equation \citep{roullet.madec_JGR00}. 268 With the present release, \NEMO offers the choice between277 With the present release, \NEMO\ offers the choice between 269 278 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or 270 279 a split-explicit scheme strongly inspired the one proposed by \citet{shchepetkin.mcwilliams_OM05} … … 275 284 % ================================================================ 276 285 \section{Curvilinear \textit{z-}coordinate system} 277 \label{sec: PE_zco}286 \label{sec:MB_zco} 278 287 279 288 % ------------------------------------------------------------------------------------------------------------- … … 281 290 % ------------------------------------------------------------------------------------------------------------- 282 291 \subsection{Tensorial formalism} 283 \label{subsec: PE_tensorial}292 \label{subsec:MB_tensorial} 284 293 285 294 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).295 (\ie\ surface layers, western boundary currents, equatorial currents, or ocean fronts). 287 296 The representation of such dynamical processes can be improved by 288 297 specifically increasing the model resolution in these regions. … … 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}323 \label{eq:MB_scale_factors} 315 324 \begin{aligned} 316 325 e_1 &= (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \\ … … 322 331 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 323 332 \begin{figure}[!tb] 324 \begin{center} 325 \includegraphics[width=\textwidth]{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} 333 \centering 334 \includegraphics[width=0.66\textwidth]{Fig_I_earth_referential} 335 \caption[Geographical and curvilinear coordinate systems]{ 336 the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 337 coordinate system $(i,j,k)$.} 338 \label{fig:MB_referential} 332 339 \end{figure} 333 340 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 334 341 335 342 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).343 (\autoref{eq:MB_scale_factors}) (thin-shell approximation). 337 344 The resulting horizontal scale factors $e_1$, $e_2$ are independent of $k$ while 338 345 the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$. 339 346 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,347 (\autoref{eq:MB_PE_dyn} to \autoref{eq:MB_PE_eos}) can then be written in the tensorial form, 341 348 invariant in any orthogonal horizontal curvilinear coordinate system transformation: 342 349 \begin{subequations} 343 % \label{eq: PE_discrete_operators}350 % \label{eq:MB_discrete_operators} 344 351 \begin{gather} 345 \label{eq: PE_grad}352 \label{eq:MB_grad} 346 353 \nabla q = \frac{1}{e_1} \pd[q]{i} \; \vect i 347 354 + \frac{1}{e_2} \pd[q]{j} \; \vect j 348 355 + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ 349 \label{eq: PE_div}356 \label{eq:MB_div} 350 357 \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 358 + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] 352 359 \end{gather} 353 360 \begin{multline} 354 \label{eq: PE_curl}361 \label{eq:MB_curl} 355 362 \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 363 + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i} \rt] \vect j \\ … … 358 365 \end{multline} 359 366 \begin{gather} 360 \label{eq: PE_lap}367 \label{eq:MB_lap} 361 368 \Delta q = \nabla \cdot (\nabla q) \\ 362 \label{eq: PE_lap_vector}369 \label{eq:MB_lap_vector} 363 370 \Delta \vect A = \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A) 364 371 \end{gather} … … 370 377 % ------------------------------------------------------------------------------------------------------------- 371 378 \subsection{Continuous model equations} 372 \label{subsec: PE_zco_Eq}379 \label{subsec:MB_zco_Eq} 373 380 374 381 In order to express the Primitive Equations in tensorial formalism, 375 382 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}.383 \autoref{eq:MB_grad}) to \autoref{eq:MB_lap_vector}. 377 384 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 385 define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 379 386 \begin{gather} 380 \label{eq: PE_curl_Uh}387 \label{eq:MB_curl_Uh} 381 388 \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}389 \label{eq:MB_div_Uh} 383 390 \chi = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, u)]{i} + \pd[(e_1 \, v)]{j} \rt] 384 391 \end{gather} 385 392 386 Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that393 Using again the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that 387 394 $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:395 $NLT$ the nonlinear term of \autoref{eq:MB_PE_dyn} can be transformed as follows: 389 396 \begin{alignat*}{2} 390 397 &NLT &= &\lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ … … 427 434 \end{alignat*} 428 435 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:436 \autoref{eq:MB_PE_dyn} is written in the $(i,j,k)$ coordinate system: 430 437 \begin{equation} 431 \label{eq: PE_vector_form}438 \label{eq:MB_vector_form} 432 439 NLT = \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) 433 440 + \frac{1}{e_3} w \pd[\vect U_h]{k} … … 436 443 This is the so-called \textit{vector invariant form} of the momentum advection term. 437 444 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:445 \ie\ to write it as the divergence of fluxes. 446 For example, the first component of \autoref{eq:MB_vector_form} (the $i$-component) is transformed as follows: 440 447 \begin{alignat*}{2} 441 448 &NLT_i &= &- \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\ … … 456 463 & &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u 457 464 + \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:}465 \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:} 459 466 & &= &\, \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 467 \end{alignat*} … … 462 469 The flux form of the momentum advection term is therefore given by: 463 470 \begin{equation} 464 \label{eq: PE_flux_form}471 \label{eq:MB_flux_form} 465 472 NLT = \nabla \cdot \lt( 466 473 \begin{array}{*{20}c} … … 475 482 the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) 476 483 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: 484 The latter is called the \textit{metric} term and can be viewed as a modification of the Coriolis parameter: 478 485 \[ 479 % \label{eq: PE_cor+metric}486 % \label{eq:MB_cor+metric} 480 487 f \to f + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) 481 488 \] 482 489 483 490 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)$,491 \ie\ when $(i,j) \to (\lambda,\varphi)$ and $(e_1,e_2) \to (a \, \cos \varphi,a)$, 485 492 we recover the commonly used modification of the Coriolis parameter $f \to f + (u / a) \tan \varphi$. 486 493 … … 492 499 \textbf{Vector invariant form of the momentum equations}: 493 500 \begin{equation} 494 \label{eq: PE_dyn_vect}501 \label{eq:MB_dyn_vect} 495 502 \begin{split} 496 % \label{eq: PE_dyn_vect_u}503 % \label{eq:MB_dyn_vect_u} 497 504 \pd[u]{t} = &+ (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) 498 505 - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 499 506 &+ D_u^{\vect U} + F_u^{\vect U} \\ 500 507 \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) \\ 508 - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 502 509 &+ D_v^{\vect U} + F_v^{\vect U} 503 510 \end{split} … … 505 512 \item 506 513 \textbf{flux form of the momentum equations}: 507 % \label{eq: PE_dyn_flux}514 % \label{eq:MB_dyn_flux} 508 515 \begin{multline*} 509 % \label{eq: PE_dyn_flux_u}516 % \label{eq:MB_dyn_flux_u} 510 517 \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 518 - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\ … … 514 521 \end{multline*} 515 522 \begin{multline*} 516 % \label{eq: PE_dyn_flux_v}523 % \label{eq:MB_dyn_flux_v} 517 524 \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) \\525 - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\ 519 526 - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 520 527 + D_v^{\vect U} + F_v^{\vect U} 521 528 \end{multline*} 522 where $\zeta$, the relative vorticity, is given by \autoref{eq: PE_curl_Uh} and $p_s$, the surface pressure,529 where $\zeta$, the relative vorticity, is given by \autoref{eq:MB_curl_Uh} and $p_s$, the surface pressure, 523 530 is given by: 524 531 \[ 525 % \label{eq: PE_spg}532 % \label{eq:MB_spg} 526 533 p_s = \rho \,g \, \eta 527 534 \] 528 with $\eta$ is solution of \autoref{eq:PE_ssh}.535 and $\eta$ is the solution of \autoref{eq:MB_ssh}. 529 536 530 537 The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: 531 538 \[ 532 % \label{eq: w_diag}533 \pd[w]{k} = - \chi \; e_3 \qquad 534 % \label{eq: hp_diag}539 % \label{eq:MB_w_diag} 540 \pd[w]{k} = - \chi \; e_3 \qquad 541 % \label{eq:MB_hp_diag} 535 542 \pd[p_h]{k} = - \rho \; g \; e_3 536 543 \] 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] 544 where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:MB_div_Uh}. 545 546 \item 547 \textbf{tracer equations}: 548 \begin{equation} 549 \begin{split} 550 \pd[T]{t} = & - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] 542 551 - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 543 %\label{eq:T} 544 \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] 545 - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S 546 %\label{eq:rho} 547 \rho = \rho \big( T,S,z(k) \big) 548 \] 552 \pd[S]{t} = & - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] 553 - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\ 554 \rho = & \rho \big( T,S,z(k) \big) 555 \end{split} 556 \end{equation} 549 557 \end{itemize} 550 558 551 559 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}.560 It will be defined in \autoref{eq:MB_zdf}. 553 561 The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms, 554 562 are discussed in \autoref{chap:SBC}. … … 560 568 % ================================================================ 561 569 \section{Curvilinear generalised vertical coordinate system} 562 \label{sec: PE_gco}570 \label{sec:MB_gco} 563 571 564 572 The ocean domain presents a huge diversity of situation in the vertical. … … 569 577 Therefore, in order to represent the ocean with respect to 570 578 the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height 571 \eg an \zstar-coordinate;579 \eg\ an \zstar-coordinate; 572 580 for the second point, a space variation to fit the change of bottom topography 573 \eg a terrain-following or $\sigma$-coordinate;581 \eg\ a terrain-following or $\sigma$-coordinate; 574 582 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 constrain s one can even be tempted to mixed these coordinate systems, as in583 follows the isopycnal surfaces, \eg\ an isopycnic coordinate. 584 585 In order to satisfy two or more constraints one can even be tempted to mixed these coordinate systems, as in 578 586 HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at 579 587 the ocean bottom) \citep{chassignet.smith.ea_JPO03} or … … 584 592 introducing an arbitrary vertical coordinate : 585 593 \begin{equation} 586 \label{eq: PE_s}594 \label{eq:MB_s} 587 595 s = s(i,j,k,t) 588 596 \end{equation} 589 597 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, 590 598 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 into599 \autoref{eq:MB_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into 592 600 the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through 593 \autoref{eq: PE_s}.601 \autoref{eq:MB_s}. 594 602 This so-called \textit{generalised vertical coordinate} \citep{kasahara_MWR74} is in fact 595 603 an Arbitrary Lagrangian--Eulerian (ALE) coordinate. 596 Indeed, choosing an expression for $s$ is an arbitrary choice thatdetermines604 Indeed, one has a great deal of freedom in the choice of expression for $s$. The choice determines 597 605 which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and 598 606 which part will be used to move them (Lagrangian part). … … 601 609 Its most often used implementation is via an ALE algorithm, 602 610 in which a pure lagrangian step is followed by regridding and remapping steps, 603 the lat er step implicitly embedding the vertical advection611 the latter step implicitly embedding the vertical advection 604 612 \citep{hirt.amsden.ea_JCP74, chassignet.smith.ea_JPO03, white.adcroft.ea_JCP09}. 605 613 Here we follow the \citep{kasahara_MWR74} strategy: 606 a regridding step (an update of the vertical coordinate) followed by an eulerian step with614 a regridding step (an update of the vertical coordinate) followed by an Eulerian step with 607 615 an explicit computation of vertical advection relative to the moving s-surfaces. 608 616 609 617 %\gmcomment{ 610 618 %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,619 The generalized vertical coordinates used in ocean modelling are not orthogonal, 612 620 which contrasts with many other applications in mathematical physics. 613 621 Hence, it is useful to keep in mind the following properties that may seem odd on initial encounter. … … 615 623 The horizontal velocity in ocean models measures motions in the horizontal plane, 616 624 perpendicular to the local gravitational field. 617 That is, horizontal velocity is mathematically the same regardless the vertical coordinate, be it geopotential,625 That is, horizontal velocity is mathematically the same regardless of the vertical coordinate, be it geopotential, 618 626 isopycnal, pressure, or terrain following. 619 627 The key motivation for maintaining the same horizontal velocity component is that … … 644 652 \subsection{\textit{S}-coordinate formulation} 645 653 646 Starting from the set of equations established in \autoref{sec: PE_zco} for the special case $k = z$ and654 Starting from the set of equations established in \autoref{sec:MB_zco} for the special case $k = z$ and 647 655 thus $e_3 = 1$, we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$, 648 656 which includes $z$-, \zstar- and $\sigma$-coordinates as special cases 649 657 ($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}.658 A formal derivation of the transformed equations is given in \autoref{apdx:SCOORD}. 651 659 Let us define the vertical scale factor by $e_3 = \partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ), 652 660 and the slopes in the $(i,j)$ directions between $s$- and $z$-surfaces by: 653 661 \begin{equation} 654 \label{eq: PE_sco_slope}662 \label{eq:MB_sco_slope} 655 663 \sigma_1 = \frac{1}{e_1} \; \lt. \pd[z]{i} \rt|_s \quad \text{and} \quad 656 664 \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt|_s 657 665 \end{equation} 658 We also introduce $\omega$, a dia-surface velocity component, defined as the velocity 666 We also introduce $\omega$, a dia-surface velocity component, defined as the velocity 659 667 relative to the moving $s$-surfaces and normal to them: 660 668 \[ 661 % \label{eq: PE_sco_w}662 \omega = w - e_3 \, \pd[s]{t}- \sigma_1 \, u - \sigma_2 \, v669 % \label{eq:MB_sco_w} 670 \omega = w - \, \lt. \pd[z]{t} \rt|_s - \sigma_1 \, u - \sigma_2 \, v 663 671 \] 664 672 665 The equations solved by the ocean model \autoref{eq: PE} in $s$-coordinate can be written as follows666 (see \autoref{sec: A_momentum}):673 The equations solved by the ocean model \autoref{eq:MB_PE} in $s$-coordinate can be written as follows 674 (see \autoref{sec:SCOORD_momentum}): 667 675 668 676 \begin{itemize} 669 677 \item \textbf{Vector invariant form of the momentum equation}: 670 678 \begin{multline*} 671 % \label{eq: PE_sco_u_vector}679 % \label{eq:MB_sco_u_vector} 672 680 \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_1681 - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 674 682 + D_u^{\vect U} + F_u^{\vect U} 675 683 \end{multline*} 676 684 \begin{multline*} 677 % \label{eq: PE_sco_v_vector}685 % \label{eq:MB_sco_v_vector} 678 686 \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_2687 - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_2 680 688 + D_v^{\vect U} + F_v^{\vect U} 681 689 \end{multline*} 682 690 \item \textbf{Flux form of the momentum equation}: 683 691 \begin{multline*} 684 % \label{eq: PE_sco_u_flux}692 % \label{eq:MB_sco_u_flux} 685 693 \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 694 - \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 695 - \frac{1}{e_3} \pd[(\omega \, u)]{k} 688 696 - \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}697 - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 690 698 \end{multline*} 691 699 \begin{multline*} 692 % \label{eq: PE_sco_v_flux}700 % \label{eq:MB_sco_v_flux} 693 701 \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 702 - \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 703 - \frac{1}{e_3} \pd[(\omega \, v)]{k} 696 704 - \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}705 - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 698 706 \end{multline*} 699 707 where the relative vorticity, $\zeta$, the surface pressure gradient, 700 708 and the hydrostatic pressure have the same expressions as in $z$-coordinates although 701 709 they do not represent exactly the same quantities. 702 $\omega$ is provided by the continuity equation (see \autoref{apdx: A}):710 $\omega$ is provided by the continuity equation (see \autoref{apdx:SCOORD}): 703 711 \[ 704 % \label{eq: PE_sco_continuity}712 % \label{eq:MB_sco_continuity} 705 713 \pd[e_3]{t} + e_3 \; \chi + \pd[\omega]{s} = 0 \quad \text{with} \quad 706 714 \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) … … 708 716 \item \textit{tracer equations}: 709 717 \begin{multline*} 710 % \label{eq: PE_sco_t}718 % \label{eq:MB_sco_t} 711 719 \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 720 + \pd[(e_1 e_3 \, v \, T)]{j} \rt) \\ … … 714 722 \end{multline*} 715 723 \begin{multline} 716 % \label{eq: PE_sco_s}724 % \label{eq:MB_sco_s} 717 725 \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 726 + \pd[(e_1 e_3 \, v \, S)]{j} \rt) \\ … … 733 741 % ------------------------------------------------------------------------------------------------------------- 734 742 \subsection{Curvilinear \zstar-coordinate system} 735 \label{subsec: PE_zco_star}743 \label{subsec:MB_zco_star} 736 744 737 745 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 738 746 \begin{figure}[!b] 739 \begin{center} 740 \includegraphics[width=\textwidth]{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 coordinate 746 (become popular as the \zstar-coordinate \citep{adcroft.campin_OM04}). 747 } 748 \end{center} 747 \centering 748 \includegraphics[width=0.66\textwidth]{Fig_z_zstar} 749 \caption[Curvilinear z-coordinate systems (\{non-\}linear free-surface cases and re-scaled \zstar)]{ 750 (a) $z$-coordinate in linear free-surface case ; 751 (b) $z$-coordinate in non-linear free surface case ; 752 (c) re-scaled height coordinate 753 (become popular as the \zstar-coordinate \citep{adcroft.campin_OM04}).} 754 \label{fig:MB_z_zstar} 749 755 \end{figure} 750 756 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 751 757 752 In th atcase, 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{levier.treguier.ea_rpt07} available on the \NEMO web site.758 In this case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. 759 These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on the \NEMO\ web site. 754 760 755 761 The \zstar coordinate approach is an unapproximated, non-linear free surface implementation which allows one to … … 759 765 as in the $z$-coordinate formulation, but is equally distributed over the full water column. 760 766 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.767 as illustrated by \autoref{fig:MB_z_zstar}. 768 Note that with a flat bottom, such as in \autoref{fig:MB_z_zstar}, the bottom-following $z$ coordinate and \zstar are equivalent. 763 769 The definition and modified oceanic equations for the rescaled vertical coordinate \zstar, 764 770 including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004). … … 766 772 The position (\zstar) and vertical discretization (\zstar) are expressed as: 767 773 \[ 768 % \label{eq: z-star}774 % \label{eq:MB_z-star} 769 775 H + \zstar = (H + z) / r \quad \text{and} \quad \delta \zstar 770 776 = \delta z / r \quad \text{with} \quad r 771 = \frac{H + \eta}{H} 777 = \frac{H + \eta}{H} . 778 \] 779 Simple re-organisation of the above expressions gives 780 \[ 781 % \label{eq:MB_zstar_2} 782 \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) . 772 783 \] 773 784 Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar, … … 776 787 Also the divergence of the flow field is no longer zero as shown by the continuity equation: 777 788 \[ 778 \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) (r \; w *) = 0789 \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 . 779 790 \] 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 791 This \zstar coordinate is closely related to the "eta" coordinate used in many atmospheric models 789 792 (see Black (1994) for a review of eta coordinate atmospheric models). 790 793 It was originally used in ocean models by Stacey et al. (1995) for studies of tides next to shelves, … … 797 800 it is clear that surfaces constant \zstar are very similar to the depth surfaces. 798 801 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 \zstarwhen $\eta = 0$,802 terrain following sigma models discussed in \autoref{subsec:MB_sco}. 803 Additionally, since $\zstar = z$ when $\eta = 0$, 801 804 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography. 802 805 This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of … … 804 807 depending on the sophistication of the pressure gradient solver. 805 808 The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of 806 neutral physics parameterizations in \zstar models using the same techniques as in $z$-models809 neutral physics parameterizations in \zstar models using the same techniques as in $z$-models 807 810 (see Chapters 13-16 of \cite{griffies_bk04}) for a discussion of neutral physics in $z$-models, 808 811 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). 809 812 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$.813 The range over which \zstar varies is time independent $-H \leq \zstar \leq 0$. 814 Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > -H$. 812 815 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 813 816 814 Because \zstar has a time independent range, all grid cells have static increments ds,815 and the sum of the ver 817 Because \zstar has a time independent range, all grid cells have static increments ds, 818 and the sum of the vertical increments yields the time independent ocean depth. %k ds = H. 816 819 The \zstar coordinate is therefore invisible to undulations of the free surface, 817 820 since it moves along with the free surface. 818 This proper ty means that no spurious vertical transport is induced across surfaces of constant \zstarby821 This property means that no spurious vertical transport is induced across surfaces of constant \zstar by 819 822 the motion of external gravity waves. 820 Such spurious transpor 821 Quite generally, the time independent range for the \zstar coordinate is a very convenient property that822 allows for a nearly arbitrary ver 823 Such spurious transport can be a problem in z-models, especially those with tidal forcing. 824 Quite generally, the time independent range for the \zstar coordinate is a very convenient property that 825 allows for a nearly arbitrary vertical resolution even in the presence of large amplitude fluctuations of 823 826 the surface height, again so long as $\eta > -H$. 824 827 %end MOM doc %%% 825 828 826 \newpage 829 \newpage 827 830 828 831 % ------------------------------------------------------------------------------------------------------------- … … 830 833 % ------------------------------------------------------------------------------------------------------------- 831 834 \subsection{Curvilinear terrain-following \textit{s}--coordinate} 832 \label{subsec: PE_sco}835 \label{subsec:MB_sco} 833 836 834 837 % ------------------------------------------------------------------------------------------------------------- … … 842 845 For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. 843 846 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}),847 In the $z$-coordinate system presented in the previous section (\autoref{sec:MB_zco}), 845 848 $z$-surfaces are geopotential surfaces. 846 849 The bottom topography is discretised by steps. … … 866 869 The main two problems come from the truncation error in the horizontal pressure gradient and 867 870 a possibly increased diapycnal diffusion. 868 The horizontal pressure force in $s$-coordinate consists of two terms (see \autoref{apdx: A}),871 The horizontal pressure force in $s$-coordinate consists of two terms (see \autoref{apdx:SCOORD}), 869 872 870 873 \begin{equation} 871 \label{eq: PE_p_sco}872 \nabla p |_z = \nabla p |_s - \ pd[p]{s} \nabla z |_s874 \label{eq:MB_p_sco} 875 \nabla p |_z = \nabla p |_s - \frac{1}{e_3} \pd[p]{s} \nabla z |_s 873 876 \end{equation} 874 877 875 The second term in \autoref{eq: PE_p_sco} depends on the tilt of the coordinate surface and876 introducesa truncation error that is not present in a $z$-model.878 The second term in \autoref{eq:MB_p_sco} depends on the tilt of the coordinate surface and 879 leads to a truncation error that is not present in a $z$-model. 877 880 In the special case of a $\sigma$-coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 878 881 \citet{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of the magnitude of this truncation error. … … 887 890 However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for 888 891 a realistic bottom topography: 889 a envelope topography is defined in $s$-coordinate on which a full or892 an envelope topography is defined in $s$-coordinate on which a full or 890 893 partial step bottom topography is then applied in order to adjust the model depth to the observed one 891 (see \autoref{s ec:DOM_zgr}.894 (see \autoref{subsec:DOM_zgr}. 892 895 893 896 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of … … 904 907 In contrast, the ocean will stay at rest in a $z$-model. 905 908 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.delecluse.ea_JPO96}.909 the strongly stratified portion of the water column (\ie\ the main thermocline) \citep{madec.delecluse.ea_JPO96}. 907 910 An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces 908 (see \autoref{subsec: PE_ldf}).911 (see \autoref{subsec:MB_ldf}). 909 912 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 910 913 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). … … 919 922 % ------------------------------------------------------------------------------------------------------------- 920 923 \subsection{\texorpdfstring{Curvilinear \ztilde-coordinate}{}} 921 \label{subsec: PE_zco_tilde}924 \label{subsec:MB_zco_tilde} 922 925 923 926 The \ztilde -coordinate has been developed by \citet{leclair.madec_OM11}. 924 It is available in \NEMO since the version 3.4.927 It is available in \NEMO\ since the version 3.4 and is more robust in version 4.0 than previously. 925 928 Nevertheless, it is currently not robust enough to be used in all possible configurations. 926 929 Its use is therefore not recommended. 927 930 928 \newpage 931 \newpage 929 932 930 933 % ================================================================ … … 932 935 % ================================================================ 933 936 \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 than937 \label{sec:MB_zdf_ldf} 938 939 The hydrostatic primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than 937 940 a few kilometres in the horizontal, a few meters in the vertical and a few minutes. 938 941 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. … … 940 943 must be represented entirely in terms of large-scale patterns to close the equations. 941 944 These effects appear in the equations as the divergence of turbulent fluxes 942 (\ie fluxes associated with the mean correlation of small scale perturbations).945 (\ie\ fluxes associated with the mean correlation of small scale perturbations). 943 946 Assuming a turbulent closure hypothesis is equivalent to choose a formulation for these fluxes. 944 947 It is usually called the subgrid scale physics. … … 949 952 The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions. 950 953 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 into954 \autoref{eq:MB_PE_dyn}, \autoref{eq:MB_PE_tra_T} and \autoref{eq:MB_PE_tra_S} are divided into 952 955 a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and 953 956 a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$. … … 958 961 % ------------------------------------------------------------------------------------------------------------- 959 962 \subsection{Vertical subgrid scale physics} 960 \label{subsec: PE_zdf}963 \label{subsec:MB_zdf} 961 964 962 965 The model resolution is always larger than the scale at which the major sources of vertical turbulence occur … … 972 975 The resulting vertical momentum and tracer diffusive operators are of second order: 973 976 \begin{equation} 974 \label{eq: PE_zdf}977 \label{eq:MB_zdf} 975 978 \begin{gathered} 976 979 \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \ , \\ … … 984 987 All the vertical physics is embedded in the specification of the eddy coefficients. 985 988 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...),989 (\eg\ Richardson number, Brunt-Vais\"{a}l\"{a} frequency, distance from the boundary ...), 987 990 or computed from a turbulent closure model. 988 The choices available in \NEMO are discussed in \autoref{chap:ZDF}).991 The choices available in \NEMO\ are discussed in \autoref{chap:ZDF}). 989 992 990 993 % ------------------------------------------------------------------------------------------------------------- … … 992 995 % ------------------------------------------------------------------------------------------------------------- 993 996 \subsection{Formulation of the lateral diffusive and viscous operators} 994 \label{subsec: PE_ldf}997 \label{subsec:MB_ldf} 995 998 996 999 Lateral turbulence can be roughly divided into a mesoscale turbulence associated with eddies … … 999 1002 and a sub mesoscale turbulence which is never explicitly solved even partially, but always parameterized. 1000 1003 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).1004 (\ie\ the model is eddy-resolving or not). 1002 1005 1003 1006 In non-eddy-resolving configurations, the closure is similar to that used for the vertical physics. … … 1007 1010 (or more precisely neutral surfaces \cite{mcdougall_JPO87}) rather than across them. 1008 1011 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.1012 the `lateral' direction is the horizontal, \ie\ the lateral mixing is performed along geopotential surfaces. 1010 1013 This leads to a geopotential second order operator for lateral subgrid scale physics. 1011 1014 This assumption can be relaxed: the eddy-induced turbulent fluxes can be better approached by assuming that … … 1014 1017 it has components in the three space directions. 1015 1018 However, 1016 both horizontal and isoneutral operators have no effect on mean (\ie large scale) potential energy whereas1019 both horizontal and isoneutral operators have no effect on mean (\ie\ large scale) potential energy whereas 1017 1020 potential energy is a main source of turbulence (through baroclinic instabilities). 1018 \citet{gent.mcwilliams_JPO90} haveproposed a parameterisation of mesoscale eddy-induced turbulence which1021 \citet{gent.mcwilliams_JPO90} proposed a parameterisation of mesoscale eddy-induced turbulence which 1019 1022 associates an eddy-induced velocity to the isoneutral diffusion. 1020 1023 Its mean effect is to reduce the mean potential energy of the ocean. … … 1033 1036 Another approach is becoming more and more popular: 1034 1037 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.1038 one uses an advective scheme which is diffusive enough to maintain the model stability. 1036 1039 It must be emphasised that then, all the sub-grid scale physics is included in the formulation of 1037 1040 the advection scheme. 1038 1041 1039 1042 All these parameterisations of subgrid scale physics have advantages and drawbacks. 1040 The reare not all available in \NEMO. For active tracers (temperature and salinity) the main ones are:1043 They are not all available in \NEMO. For active tracers (temperature and salinity) the main ones are: 1041 1044 Laplacian and bilaplacian operators acting along geopotential or iso-neutral surfaces, 1042 1045 \citet{gent.mcwilliams_JPO90} parameterisation, and various slightly diffusive advection schemes. … … 1046 1049 \subsubsection{Lateral laplacian tracer diffusive operator} 1047 1050 1048 The lateral Laplacian tracer diffusive operator is defined by (see \autoref{apdx: B}):1051 The lateral Laplacian tracer diffusive operator is defined by (see \autoref{apdx:DIFFOPERS}): 1049 1052 \begin{equation} 1050 \label{eq: PE_iso_tensor}1053 \label{eq:MB_iso_tensor} 1051 1054 D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad 1052 1055 \Re = … … 1058 1061 \end{equation} 1059 1062 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 for1063 the model level (\eg\ $z$- or $s$-surfaces). 1064 Note that the formulation \autoref{eq:MB_iso_tensor} is exact for 1062 1065 the rotation between geopotential and $s$-surfaces, 1063 1066 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{cox_OM87}.1067 Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:MB_iso_tensor} \citep{cox_OM87}. 1065 1068 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and 1066 1069 dia-neutral diffusive coefficients is known to be several orders of magnitude smaller than unity. 1067 1070 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}).1071 the slopes are generally less than $10^{-2}$ in the ocean (see \autoref{apdx:DIFFOPERS}). 1069 1072 1070 1073 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. … … 1073 1076 For \textit{geopotential} diffusion, 1074 1077 $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}).1078 they are equal to $\sigma_1$ and $\sigma_2$, respectively (see \autoref{eq:MB_sco_slope}). 1076 1079 1077 1080 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral and computational surfaces. … … 1079 1082 In $z$-coordinates: 1080 1083 \begin{equation} 1081 \label{eq: PE_iso_slopes}1084 \label{eq:MB_iso_slopes} 1082 1085 r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 1083 1086 r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} … … 1090 1093 an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: 1091 1094 \[ 1092 % \label{eq: PE_iso+eiv}1095 % \label{eq:MB_iso+eiv} 1093 1096 D^{lT} = \nabla \cdot \lt( A^{lT} \; \Re \; \nabla T \rt) + \nabla \cdot \lt( \vect U^\ast \, T \rt) 1094 1097 \] … … 1096 1099 eddy-induced transport velocity. This velocity field is defined by: 1097 1100 \begin{gather} 1098 % \label{eq: PE_eiv}1101 % \label{eq:MB_eiv} 1099 1102 u^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_1 \rt) \\ 1100 1103 v^\ast = \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_2 \rt) \\ … … 1105 1108 (or equivalently the isoneutral thickness diffusivity coefficient), 1106 1109 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: 1110 Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate: 1108 1111 \begin{align} 1109 \label{eq: PE_slopes_eiv}1112 \label{eq:MB_slopes_eiv} 1110 1113 \tilde{r}_n = 1111 1114 \begin{cases} … … 1119 1122 This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of 1120 1123 the boundaries. 1121 The latter strategy is used in \NEMO (cf. \autoref{chap:LDF}).1124 The latter strategy is used in \NEMO\ (cf. \autoref{chap:LDF}). 1122 1125 1123 1126 \subsubsection{Lateral bilaplacian tracer diffusive operator} … … 1125 1128 The lateral bilaplacian tracer diffusive operator is defined by: 1126 1129 \[ 1127 % \label{eq: PE_bilapT}1130 % \label{eq:MB_bilapT} 1128 1131 D^{lT}= - \Delta \; (\Delta T) \quad \text{where} \quad 1129 1132 \Delta \bullet = \nabla \lt( \sqrt{B^{lT}} \; \Re \; \nabla \bullet \rt) 1130 1133 \] 1131 It is the Laplacian operator given by \autoref{eq: PE_iso_tensor} applied twice with1134 It is the Laplacian operator given by \autoref{eq:MB_iso_tensor} applied twice with 1132 1135 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1133 1136 … … 1135 1138 1136 1139 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}):1140 applying \autoref{eq:MB_lap_vector} to the horizontal velocity vector (see \autoref{apdx:DIFFOPERS}): 1138 1141 \begin{align*} 1139 % \label{eq: PE_lapU}1142 % \label{eq:MB_lapU} 1140 1143 \vect D^{l \vect U} &= \nabla_h \big( A^{lm} \chi \big) 1141 1144 - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\ 1142 1145 &= \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} 1146 - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j} , 1144 1147 \frac{1}{e_2} \pd[ \lt( A^{lm} \chi \rt) ]{j} 1145 1148 \lt. + \frac{1}{e_1 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{i} \rt) … … 1147 1150 1148 1151 Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields 1149 (see \autoref{apdx: C}).1152 (see \autoref{apdx:INVARIANTS}). 1150 1153 Unfortunately, it is only available in \textit{iso-level} direction. 1151 1154 When a rotation is required 1152 (\ie geopotential diffusion in $s$-coordinates or isoneutral diffusion in both $z$- and $s$-coordinates),1155 (\ie\ geopotential diffusion in $s$-coordinates or isoneutral diffusion in both $z$- and $s$-coordinates), 1153 1156 the $u$ and $v$-fields are considered as independent scalar fields, so that the diffusive operator is given by: 1154 1157 \begin{gather*} 1155 % \label{eq: PE_lapU_iso}1158 % \label{eq:MB_lapU_iso} 1156 1159 D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \\ 1157 1160 D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt) 1158 1161 \end{gather*} 1159 where $\Re$ is given by \autoref{eq: PE_iso_tensor}.1162 where $\Re$ is given by \autoref{eq:MB_iso_tensor}. 1160 1163 It is the same expression as those used for diffusive operator on tracers. 1161 1164 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.1165 \ie\ on a $f$- or $\beta$-plane, not on the sphere. 1163 1166 It is also a very good approximation in vicinity of the Equator in 1164 1167 a geographical coordinate system \citep{lengaigne.madec.ea_JGR03}. 1165 1168 1166 \subsubsection{ lateral bilaplacian momentum diffusive operator}1169 \subsubsection{Lateral bilaplacian momentum diffusive operator} 1167 1170 1168 1171 As for tracers, the bilaplacian order momentum diffusive operator is a re-entering Laplacian operator with
Note: See TracChangeset
for help on using the changeset viewer.
