- Timestamp:
- 2019-01-10T17:30:42+01:00 (5 years ago)
- Location:
- NEMO/trunk/doc/latex/NEMO/subfiles
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
r10442 r10501 2 2 3 3 \begin{document} 4 4 5 % ================================================================ 5 % Chapter 1 Model Basics6 % Chapter 1 Model Basics 6 7 % ================================================================ 7 8 8 \chapter{Model Basics} 9 9 \label{chap:PE} 10 11 10 \minitoc 12 11 … … 26 25 \label{subsec:PE_Vector} 27 26 28 29 27 The ocean is a fluid that can be described to a good approximation by the primitive equations, 30 28 \ie the Navier-Stokes equations along with a nonlinear equation of state which … … 32 30 plus the following additional assumptions made from scale considerations: 33 31 34 \textit{(1) spherical earth approximation:} the geopotential surfaces are assumed to be spheres so that 35 gravity (local vertical) is parallel to the earth's radius 36 37 \textit{(2) thin-shell approximation:} the ocean depth is neglected compared to the earth's radius 38 39 \textit{(3) turbulent closure hypothesis:} the turbulent fluxes 40 (which represent the effect of small scale processes on the large-scale) are expressed in terms of 41 large-scale features 42 43 \textit{(4) Boussinesq hypothesis:} density variations are neglected except in their contribution to 44 the buoyancy force 45 46 \textit{(5) Hydrostatic hypothesis:} the vertical momentum equation is reduced to a balance between 47 the vertical pressure gradient and the buoyancy force 48 (this removes convective processes from the initial Navier-Stokes equations and so 49 convective processes must be parameterized instead) 50 51 \textit{(6) Incompressibility hypothesis:} the three dimensional divergence of the velocity vector is assumed to 52 be zero. 32 \begin{enumerate} 33 \item 34 \textit{spherical earth approximation}: the geopotential surfaces are assumed to be spheres so that 35 gravity (local vertical) is parallel to the earth's radius 36 \item 37 \textit{thin-shell approximation}: the ocean depth is neglected compared to the earth's radius 38 \item 39 \textit{turbulent closure hypothesis}: the turbulent fluxes 40 (which represent the effect of small scale processes on the large-scale) 41 are expressed in terms of large-scale features 42 \item 43 \textit{Boussinesq hypothesis}: density variations are neglected except in their contribution to 44 the buoyancy force 45 \begin{equation} 46 \label{eq:PE_eos} 47 \rho = \rho \ (T,S,p) 48 \end{equation} 49 \item 50 \textit{Hydrostatic hypothesis}: the vertical momentum equation is reduced to a balance between 51 the vertical pressure gradient and the buoyancy force 52 (this removes convective processes from the initial Navier-Stokes equations and so 53 convective processes must be parameterized instead) 54 \begin{equation} 55 \label{eq:PE_hydrostatic} 56 \pd[p]{z} = - \rho \ g 57 \end{equation} 58 \item 59 \textit{Incompressibility hypothesis}: the three dimensional divergence of the velocity vector $\vect U$ 60 is assumed to be zero. 61 \begin{equation} 62 \label{eq:PE_continuity} 63 \nabla \cdot \vect U = 0 64 \end{equation} 65 \end{enumerate} 53 66 54 67 Because the gravitational force is so dominant in the equations of large-scale motions, 55 it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to 56 the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to 57 \textbf{k}, \ie tangent to the geopotential surfaces. 58 Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ 59 (the subscript $h$ denotes the local horizontal vector, \ie over the (\textbf{i},\textbf{j}) plane), 60 $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. 61 The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system 62 provides the following six equations 63 (namely the momentum balance, the hydrostatic equilibrium, the incompressibility equation, 64 the heat and salt conservation equations and an equation of state): 68 it is useful to choose an orthogonal set of unit vectors $(i,j,k)$ linked to the earth such that 69 $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 70 \ie tangent to the geopotential surfaces. 71 Let us define the following variables: $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ 72 (the subscript $h$ denotes the local horizontal vector, \ie over the $(i,j)$ plane), 73 $T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density. 74 The vector invariant form of the primitive equations in the $(i,j,k)$ vector system provides 75 the following equations: 65 76 \begin{subequations} 66 77 \label{eq:PE} 67 \begin{equation} 78 \begin{gather} 79 \intertext{$-$ the momentum balance} 68 80 \label{eq:PE_dyn} 69 \frac{\partial {\rm {\bf U}}_h }{\partial t}= 70 -\left[ {\left( {\nabla \times {\rm {\bf U}}} \right)\times {\rm {\bf U}} 71 +\frac{1}{2}\nabla \left( {{\rm {\bf U}}^2} \right)} \right]_h 72 -f\;{\rm {\bf k}}\times {\rm {\bf U}}_h 73 -\frac{1}{\rho_o }\nabla _h p + {\rm {\bf D}}^{\rm {\bf U}} + {\rm {\bf F}}^{\rm {\bf U}} 74 \end{equation} 75 \begin{equation} 76 \label{eq:PE_hydrostatic} 77 \frac{\partial p }{\partial z} = - \rho \ g 78 \end{equation} 79 \begin{equation} 80 \label{eq:PE_continuity} 81 \nabla \cdot {\bf U}= 0 82 \end{equation} 83 \begin{equation} 81 \pd[\vect U_h]{t} = - \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h 82 - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p 83 + \vect D^{\vect U} + \vect F^{\vect U} \\ 84 \intertext{$-$ the heat and salt conservation equations} 84 85 \label{eq:PE_tra_T} 85 \frac{\partial T}{\partial t} = - \nabla \cdot \left( T \ \rm{\bf U} \right) + D^T + F^T 86 \end{equation} 87 \begin{equation} 86 \pd[T]{t} = - \nabla \cdot (T \ \vect U) + D^T + F^T \\ 88 87 \label{eq:PE_tra_S} 89 \frac{\partial S}{\partial t} = - \nabla \cdot \left( S \ \rm{\bf U} \right) + D^S + F^S 90 \end{equation} 91 \begin{equation} 92 \label{eq:PE_eos} 93 \rho = \rho \left( T,S,p \right) 94 \end{equation} 88 \pd[S]{t} = - \nabla \cdot (S \ \vect U) + D^S + F^S 89 \end{gather} 95 90 \end{subequations} 96 where $\nabla$ is the generalised derivative vector operator in $( \bf i,\bf j, \bfk)$ directions, $t$ is the time,97 $z$ is the vertical coordinate, $\rho 91 where $\nabla$ is the generalised derivative vector operator in $(i,j,k)$ directions, $t$ is the time, 92 $z$ is the vertical coordinate, $\rho$ is the \textit{in situ} density given by the equation of state 98 93 (\autoref{eq:PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 99 $f =2 \bf \Omega \cdot \bfk$ is the Coriolis acceleration100 (where $\ bf\Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration.101 $ {\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum,102 temperature and salinity, and $ {\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$ surface forcing terms.94 $f = 2 \vect \Omega \cdot k$ is the Coriolis acceleration 95 (where $\vect \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 96 $\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, 97 temperature and salinity, and $\vect F^{\vect U}$, $F^T$ and $F^S$ surface forcing terms. 103 98 Their nature and formulation are discussed in \autoref{sec:PE_zdf_ldf} and \autoref{subsec:PE_boundary_condition}. 104 105 106 99 107 100 % ------------------------------------------------------------------------------------------------------------- … … 113 106 An ocean is bounded by complex coastlines, bottom topography at its base and 114 107 an air-sea or ice-sea interface at its top. 115 These boundaries can be defined by two surfaces, $z =-H(i,j)$ and $z=\eta(i,j,k,t)$,108 These boundaries can be defined by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,k,t)$, 116 109 where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface. 117 Both $H$ and $\eta$ are usually referenced to a given surface, $z =0$, chosen as a mean sea surface110 Both $H$ and $\eta$ are usually referenced to a given surface, $z = 0$, chosen as a mean sea surface 118 111 (\autoref{fig:ocean_bc}). 119 112 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with … … 127 120 \begin{figure}[!ht] 128 121 \begin{center} 129 \includegraphics[width=0.90\textwidth]{Fig_I_ocean_bc} 130 \caption{ \protect\label{fig:ocean_bc} 131 The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,t)$, 122 \includegraphics[]{Fig_I_ocean_bc} 123 \caption{ 124 \protect\label{fig:ocean_bc} 125 The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, 132 126 where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. 133 Both $H$ and $\eta$ are referenced to $z =0$.127 Both $H$ and $\eta$ are referenced to $z = 0$. 134 128 } 135 129 \end{center} 136 130 \end{figure} 137 131 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 138 139 132 140 133 \begin{description} … … 161 154 \begin{equation} 162 155 \label{eq:PE_w_bbc} 163 w = - {\rm {\bf U}}_h \cdot \nabla _h \left( H \right)156 w = - \vect U_h \cdot \nabla_h (H) 164 157 \end{equation} 165 158 In addition, the ocean exchanges momentum with the earth through frictional processes. … … 167 160 It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. 168 161 Its specification depends on the nature of the physical parameterisation used for 169 $ {\rm {\bf D}}^{\rm {\bf U}}$ in \autoref{eq:PE_dyn}.162 $\vect D^{\vect U}$ in \autoref{eq:PE_dyn}. 170 163 It is discussed in \autoref{eq:PE_zdf}.% and Chap. III.6 to 9. 171 164 \item[Atmosphere - ocean interface:] … … 174 167 \[ 175 168 % \label{eq:PE_w_sbc} 176 w = \frac{\partial \eta }{\partial t} 177 + \left. {{\rm {\bf U}}_h } \right|_{z=\eta } \cdot \nabla _h \left( \eta \right) 178 + P-E 169 w = \pd[\eta]{t} + \lt. \vect U_h \rt|_{z = \eta} \cdot \nabla_h (\eta) + P - E 179 170 \] 180 171 The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system) 181 leads to the continuity of pressure across the interface $z =\eta$.172 leads to the continuity of pressure across the interface $z = \eta$. 182 173 The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. 183 174 \item[Sea ice - ocean interface:] 184 175 the ocean and sea ice exchange heat, salt, fresh water and momentum. 185 176 The sea surface temperature is constrained to be at the freezing point at the interface. 186 Sea ice salinity is very low ($\sim4-6 \, psu$) compared to those of the ocean ($\sim34 \,psu$).177 Sea ice salinity is very low ($\sim4-6 \, psu$) compared to those of the ocean ($\sim34 \, psu$). 187 178 The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected. 188 179 \end{description} 189 190 191 %\newpage192 180 193 181 % ================================================================ 194 182 % The Horizontal Pressure Gradient 195 183 % ================================================================ 196 \section{Horizontal pressure gradient 184 \section{Horizontal pressure gradient} 197 185 \label{sec:PE_hor_pg} 198 186 … … 204 192 205 193 The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at 206 a reference geopotential surface ($z =0$) and a hydrostatic pressure $p_h$ such that:207 $p(i,j,k,t) =p_s(i,j,t)+p_h(i,j,k,t)$.194 a reference geopotential surface ($z = 0$) and a hydrostatic pressure $p_h$ such that: 195 $p(i,j,k,t) = p_s(i,j,t) + p_h(i,j,k,t)$. 208 196 The latter is computed by integrating (\autoref{eq:PE_hydrostatic}), 209 197 assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:PE_eos}). … … 211 199 \[ 212 200 % \label{eq:PE_pressure} 213 p_h \left( {i,j,z,t} \right) 214 = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma } 201 p_h (i,j,z,t) = \int_{\varsigma = z}^{\varsigma = 0} g \; \rho (T,S,\varsigma) \; d \varsigma 215 202 \] 216 203 Two strategies can be considered for the surface pressure term: … … 224 211 The phase speed of such waves is high (some hundreds of metres per second) so that 225 212 the time step would have to be very short if they were present in the model. 226 The latter strategy filters out these waves since the rigid lid approximation implies $\eta =0$,227 \ie the sea surface is the surface $z =0$.213 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$. 228 215 This well known approximation increases the surface wave speed to infinity and 229 216 modifies certain other longwave dynamics (\eg barotropic Rossby or planetary waves). 230 217 The rigid-lid hypothesis is an obsolescent feature in modern OGCMs. 231 It has been available until the release 3.1 of 218 It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings. 232 219 Only the free surface formulation is now described in the this document (see the next sub-section). 233 220 … … 244 231 \begin{equation} 245 232 \label{eq:PE_ssh} 246 \frac{\partial \eta }{\partial t}=-D+P-E 247 \quad \text{where} \ 248 D=\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right] 233 \pd[\eta]{t} = - D + P - E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt] 249 234 \end{equation} 250 235 and using (\autoref{eq:PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. … … 256 241 257 242 Two choices can be made regarding the implementation of the free surface in the model, 258 depending on the physical processes of interest. 243 depending on the physical processes of interest. 259 244 260 245 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction with 261 246 the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 262 247 then a non linear free surface is the most appropriate. 263 This means that no approximation is made in (\autoref{eq:PE_ssh})and that248 This means that no approximation is made in \autoref{eq:PE_ssh} and that 264 249 the variation of the ocean volume is fully taken into account. 265 250 Note that in order to study the fast time scales associated with EGWs it is necessary to … … 272 257 not altering the slow barotropic Rossby waves. 273 258 If further, an approximative conservation of heat and salt contents is sufficient for the problem solved, 274 then it is sufficient to solve a linearized version of (\autoref{eq:PE_ssh}),259 then it is sufficient to solve a linearized version of \autoref{eq:PE_ssh}, 275 260 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{Roullet_Madec_JGR00}. 276 261 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. … … 286 271 (see \autoref{subsec:DYN_spg_ts}). 287 272 288 %\newpage289 290 273 % ================================================================ 291 274 % Curvilinear z-coordinate System … … 293 276 \section{Curvilinear \textit{z-}coordinate system} 294 277 \label{sec:PE_zco} 295 296 278 297 279 % ------------------------------------------------------------------------------------------------------------- … … 318 300 The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics. 319 301 320 Let (\textit{i},\textit{j},\textit{k})be a set of orthogonal curvilinear coordinates on302 Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on 321 303 the sphere associated with the positively oriented orthogonal set of unit vectors 322 (\textbf{i},\textbf{j},\textbf{k})linked to the earth such that323 \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k},304 $(i,j,k)$ linked to the earth such that 305 $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 324 306 \ie along geopotential surfaces (\autoref{fig:referential}). 325 307 Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by 326 308 the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and 327 the distance from the centre of the earth $a +z(k)$ where $a$ is the earth's radius and309 the distance from the centre of the earth $a + z(k)$ where $a$ is the earth's radius and 328 310 $z$ the altitude above a reference sea level (\autoref{fig:referential}). 329 311 The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, … … 332 314 \label{eq:scale_factors} 333 315 \begin{aligned} 334 e_1 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda}{\partial i}\cos \varphi } \right)^2 335 +\left( {\frac{\partial \varphi }{\partial i}} \right)^2} \right]^{1/2} \\ 336 e_2 &=\left( {a+z} \right)\;\left[ {\left( {\frac{\partial \lambda }{\partial j}\cos \varphi } \right)^2+ 337 \left( {\frac{\partial \varphi }{\partial j}} \right)^2} \right]^{1/2} \\ 338 e_3 &=\left( {\frac{\partial z}{\partial k}} \right) \\ 316 e_1 &= (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \\ 317 e_2 &= (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \\ 318 e_3 &= \lt( \pd[z]{k} \rt) 339 319 \end{aligned} 340 320 \end{equation} … … 343 323 \begin{figure}[!tb] 344 324 \begin{center} 345 \includegraphics[width=0.60\textwidth]{Fig_I_earth_referential} 346 \caption{ \protect\label{fig:referential} 325 \includegraphics[]{Fig_I_earth_referential} 326 \caption{ 327 \protect\label{fig:referential} 347 328 the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 348 coordinate system (\textbf{i},\textbf{j},\textbf{k}). } 329 coordinate system $(i,j,k)$. 330 } 349 331 \end{center} 350 332 \end{figure} 351 333 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 352 334 353 Since the ocean depth is far smaller than the earth's radius, $a +z$, can be replaced by $a$ in335 Since the ocean depth is far smaller than the earth's radius, $a + z$, can be replaced by $a$ in 354 336 (\autoref{eq:scale_factors}) (thin-shell approximation). 355 337 The resulting horizontal scale factors $e_1$, $e_2$ are independent of $k$ while 356 the vertical scale factor is a single function of $k$ as \textbf{k} is parallel to \textbf{z}.338 the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$. 357 339 The scalar and vector operators that appear in the primitive equations 358 340 (\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can be written in the tensorial form, … … 360 342 \begin{subequations} 361 343 % \label{eq:PE_discrete_operators} 362 \begin{ equation}344 \begin{gather} 363 345 \label{eq:PE_grad} 364 \nabla q=\frac{1}{e_1 }\frac{\partial q}{\partial i}\;{\rm {\bf 365 i}}+\frac{1}{e_2 }\frac{\partial q}{\partial j}\;{\rm {\bf j}}+\frac{1}{e_3 366 }\frac{\partial q}{\partial k}\;{\rm {\bf k}} \\ 367 \end{equation} 368 \begin{equation} 346 \nabla q = \frac{1}{e_1} \pd[q]{i} \; \vect i 347 + \frac{1}{e_2} \pd[q]{j} \; \vect j 348 + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ 369 349 \label{eq:PE_div} 370 \nabla \cdot {\rm {\bf A}} 371 = \frac{1}{e_1 \; e_2} \left[ 372 \frac{\partial \left(e_2 \; a_1\right)}{\partial i } 373 +\frac{\partial \left(e_1 \; a_2\right)}{\partial j } \right] 374 + \frac{1}{e_3} \left[ \frac{\partial a_3}{\partial k } \right] 375 \end{equation} 376 \begin{equation} 350 \nabla \cdot \vect A = \frac{1}{e_1 \; e_2} \lt[ \pd[(e_2 \; a_1)]{\partial i} + \pd[(e_1 \; a_2)]{j} \rt] 351 + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] 352 \end{gather} 353 \begin{multline} 377 354 \label{eq:PE_curl} 378 \begin{split} 379 \nabla \times \vect{A} = 380 \left[ {\frac{1}{e_2 }\frac{\partial a_3}{\partial j} 381 -\frac{1}{e_3 }\frac{\partial a_2 }{\partial k}} \right] \; \vect{i} 382 &+\left[ {\frac{1}{e_3 }\frac{\partial a_1 }{\partial k} 383 -\frac{1}{e_1 }\frac{\partial a_3 }{\partial i}} \right] \; \vect{j} \\ 384 &+\frac{1}{e_1 e_2 } \left[ {\frac{\partial \left( {e_2 a_2 } \right)}{\partial i} 385 -\frac{\partial \left( {e_1 a_1 } \right)}{\partial j}} \right] \; \vect{k} 386 \end{split} 387 \end{equation} 388 \begin{equation} 355 \nabla \times \vect{A} = \lt[ \frac{1}{e_2} \pd[a_3]{j} - \frac{1}{e_3} \pd[a_2]{k} \rt] \vect i \\ 356 + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i} \rt] \vect j \\ 357 + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i} - \pd[(e_1 a_1)]{j} \rt] \vect k 358 \end{multline} 359 \begin{gather} 389 360 \label{eq:PE_lap} 390 \Delta q = \nabla \cdot \left( \nabla q \right) 391 \end{equation} 392 \begin{equation} 361 \Delta q = \nabla \cdot (\nabla q) \\ 393 362 \label{eq:PE_lap_vector} 394 \Delta {\rm {\bf A}} = 395 \nabla \left( \nabla \cdot {\rm {\bf A}} \right) 396 - \nabla \times \left( \nabla \times {\rm {\bf A}} \right) 397 \end{equation} 363 \Delta \vect A = \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A) 364 \end{gather} 398 365 \end{subequations} 399 where $q$ is a scalar quantity and $ {\rm {\bf A}}=(a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinatesystem.366 where $q$ is a scalar quantity and $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system. 400 367 401 368 % ------------------------------------------------------------------------------------------------------------- … … 408 375 it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using 409 376 \autoref{eq:PE_grad}) to \autoref{eq:PE_lap_vector}. 410 Let us set $\vect U =(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinatesystem and377 Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, the velocity in the $(i,j,k)$ coordinates system and 411 378 define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 412 \begin{ equation}379 \begin{gather} 413 380 \label{eq:PE_curl_Uh} 414 \zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v} 415 \right)}{\partial i}-\frac{\partial \left( {e_1 \,u} \right)}{\partial j}} 416 \right] 417 \end{equation} 418 \begin{equation} 381 \zeta = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, v)]{i} - \pd[(e_1 \, u)]{j} \rt] \\ 419 382 \label{eq:PE_div_Uh} 420 \chi =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,u} 421 \right)}{\partial i}+\frac{\partial \left( {e_1 \,v} \right)}{\partial j}} 422 \right] 423 \end{equation} 383 \chi = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, u)]{i} + \pd[(e_1 \, v)]{j} \rt] 384 \end{gather} 424 385 425 386 Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that 426 387 $e_3$ is a function of the single variable $k$, 427 the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows: 428 \begin{flalign*} 429 &\left[ {\left( { \nabla \times {\rm {\bf U}} } \right) \times {\rm {\bf U}} 430 +\frac{1}{2} \nabla \left( {{\rm {\bf U}}^2} \right)} \right]_h & 431 \end{flalign*} 432 \begin{flalign*} 433 &\qquad=\left( {{ 434 \begin{array}{*{20}c} 435 {\left[ { \frac{1}{e_3} \frac{\partial u }{\partial k} 436 -\frac{1}{e_1} \frac{\partial w }{\partial i} } \right] w - \zeta \; v } \\ 437 {\zeta \; u - \left[ { \frac{1}{e_2} \frac{\partial w}{\partial j} 438 -\frac{1}{e_3} \frac{\partial v}{\partial k} } \right] \ w} \\ 439 \end{array} 440 }} \right) 441 +\frac{1}{2} \left( {{ 442 \begin{array}{*{20}c} 443 { \frac{1}{e_1} \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial i}} \hfill \\ 444 { \frac{1}{e_2} \frac{\partial \left( u^2+v^2+w^2 \right)}{\partial j}} \hfill \\ 445 \end{array} 446 }} \right) & 447 \end{flalign*} 448 \begin{flalign*} 449 & \qquad =\left( {{ 450 \begin{array}{*{20}c} 451 {-\zeta \; v} \hfill \\ 452 { \zeta \; u} \hfill \\ 453 \end{array} 454 }} \right) 455 +\frac{1}{2}\left( {{ 456 \begin{array}{*{20}c} 457 {\frac{1}{e_1 }\frac{\partial \left( {u^2+v^2} \right)}{\partial i}} \hfill \\ 458 {\frac{1}{e_2 }\frac{\partial \left( {u^2+v^2} \right)}{\partial j}} \hfill \\ 459 \end{array} 460 }} \right) 461 +\frac{1}{e_3 }\left( {{ 462 \begin{array}{*{20}c} 463 { w \; \frac{\partial u}{\partial k}} \\ 464 { w \; \frac{\partial v}{\partial k}} \\ 465 \end{array} 466 }} \right) 467 -\left( {{ 468 \begin{array}{*{20}c} 469 {\frac{w}{e_1}\frac{\partial w}{\partial i} -\frac{1}{2e_1}\frac{\partial w^2}{\partial i}} \hfill \\ 470 {\frac{w}{e_2}\frac{\partial w}{\partial j} -\frac{1}{2e_2}\frac{\partial w^2}{\partial j}} \hfill \\ 471 \end{array} 472 }} \right) & 473 \end{flalign*} 474 388 $NLT$ the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows: 389 \begin{alignat*}{2} 390 &NLT &= &\lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ 391 & &= &\lt( 392 \begin{array}{*{20}c} 393 \lt[ \frac{1}{e_3} \pd[u]{k} - \frac{1}{e_1} \pd[w]{i} \rt] w - \zeta \; v \\ 394 \zeta \; u - \lt[ \frac{1}{e_2} \pd[w]{j} - \frac{1}{e_3} \pd[v]{k} \rt] \ w 395 \end{array} 396 \rt) 397 + \frac{1}{2} \lt( 398 \begin{array}{*{20}c} 399 \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\ 400 \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j} 401 \end{array} 402 \rt) \\ 403 & &= &\lt( 404 \begin{array}{*{20}c} 405 -\zeta \; v \\ 406 \zeta \; u 407 \end{array} 408 \rt) 409 + \frac{1}{2} \lt( 410 \begin{array}{*{20}c} 411 \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\ 412 \frac{1}{e_2} \pd[(u^2 + v^2)]{j} 413 \end{array} 414 \rt) \\ 415 & & &+ \frac{1}{e_3} \lt( 416 \begin{array}{*{20}c} 417 w \; \pd[u]{k} \\ 418 w \; \pd[v]{k} 419 \end{array} 420 \rt) 421 - \lt( 422 \begin{array}{*{20}c} 423 \frac{w}{e_1} \pd[w]{i} - \frac{1}{2 e_1} \pd[w^2]{i} \\ 424 \frac{w}{e_2} \pd[w]{j} - \frac{1}{2 e_2} \pd[w^2]{j} 425 \end{array} 426 \rt) 427 \end{alignat*} 475 428 The last term of the right hand side is obviously zero, and thus the nonlinear term of 476 429 \autoref{eq:PE_dyn} is written in the $(i,j,k)$ coordinate system: 477 430 \begin{equation} 478 431 \label{eq:PE_vector_form} 479 \left[ {\left( { \nabla \times {\rm {\bf U}} } \right) \times {\rm {\bf U}} 480 +\frac{1}{2} \nabla \left( {{\rm {\bf U}}^2} \right)} \right]_h 481 =\zeta 482 \;{\rm {\bf k}}\times {\rm {\bf U}}_h +\frac{1}{2}\nabla _h \left( {{\rm 483 {\bf U}}_h^2 } \right)+\frac{1}{e_3 }w\frac{\partial {\rm {\bf U}}_h 484 }{\partial k} 432 NLT = \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) 433 + \frac{1}{e_3} w \pd[\vect U_h]{k} 485 434 \end{equation} 486 435 … … 489 438 \ie to write it as the divergence of fluxes. 490 439 For example, the first component of \autoref{eq:PE_vector_form} (the $i$-component) is transformed as follows: 491 \begin{flalign*} 492 &{ 493 \begin{array}{*{20}l} 494 \left[ {\left( {\nabla \times \vect{U}} \right)\times \vect{U} 495 +\frac{1}{2}\nabla \left( {\vect{U}}^2 \right)} \right]_i % \\ 496 % \\ 497 = - \zeta \;v 498 + \frac{1}{2\;e_1 } \frac{\partial \left( {u^2+v^2} \right)}{\partial i} 499 + \frac{1}{e_3}w \ \frac{\partial u}{\partial k} \\ \\ 500 \qquad =\frac{1}{e_1 \; e_2} \left( -v\frac{\partial \left( {e_2 \,v} \right)}{\partial i} 501 +v\frac{\partial \left( {e_1 \,u} \right)}{\partial j} \right) 502 +\frac{1}{e_1 e_2 }\left( +e_2 \; u\frac{\partial u}{\partial i} 503 +e_2 \; v\frac{\partial v}{\partial i} \right) 504 +\frac{1}{e_3} \left( w\;\frac{\partial u}{\partial k} \right) \\ 440 \begin{alignat*}{2} 441 &NLT_i &= &- \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\ 442 & &= &\frac{1}{e_1 \; e_2} \lt( -v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt) 443 + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) \\ 444 & & &+ \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\ 445 & &= &\frac{1}{e_1 \; e_2} \lt[ - \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt) 446 + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j} - e_1 \, u \pd[v]{j} \rt) \rt. \\ 447 & & &\lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i} - u \pd[ \lt( e_2 u \rt)]{i} \rt) 448 + e_2 v \pd[v]{i} \rt] \\ 449 & & &+ \frac{1}{e_3} \lt( \pd[(w \, u)]{k} - u \pd[w]{k} \rt) \\ 450 & &= &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt) 451 + \frac{1}{e_3} \pd[(w \, u)]{k} \\ 452 & & &+ \frac{1}{e_1 e_2} \lt[ - u \lt( \pd[(e_1 v)]{j} - v \, \pd[e_1]{j} \rt) 453 - u \pd[(e_2 u)]{i} \rt] 454 - \frac{1}{e_3} \pd[w]{k} u \\ 455 & & &+ \frac{1}{e_1 e_2} \lt( - v^2 \pd[e_2]{i} \rt) \\ 456 & &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u 457 + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ 458 \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it comes:} 459 & &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i} - u \; \pd[e_1]{j} \rt) (-v) 460 \end{alignat*} 461 462 The flux form of the momentum advection term is therefore given by: 463 \begin{equation} 464 \label{eq:PE_flux_form} 465 NLT = \nabla \cdot \lt( 466 \begin{array}{*{20}c} 467 \vect U \, u \\ 468 \vect U \, v 505 469 \end{array} 506 } & 507 \end{flalign*} 508 \begin{flalign*} 509 &{ 510 \begin{array}{*{20}l} 511 \qquad =\frac{1}{e_1 \; e_2} \left\{ 512 -\left( v^2 \frac{\partial e_2 }{\partial i} 513 +e_2 \,v \frac{\partial v }{\partial i} \right) 514 +\left( \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j} 515 -e_1 \,u \frac{\partial v }{\partial j} \right) \right. \\ 516 \left. \qquad \qquad \quad 517 +\left( \frac{\partial \left( {e_2 u\,u} \right)}{\partial i} 518 -u \frac{\partial \left( {e_2 u} \right)}{\partial i} \right) 519 +e_2 v \frac{\partial v }{\partial i} 520 \right\} 521 +\frac{1}{e_3} \left( 522 \frac{\partial \left( {w\,u} \right) }{\partial k} 523 -u \frac{\partial w }{\partial k} \right) \\ 524 \end{array} 525 } & 526 \end{flalign*} 527 \begin{flalign*} 528 & 529 { 530 \begin{array}{*{20}l} 531 \qquad =\frac{1}{e_1 \; e_2} \left( 532 \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 533 + \frac{\partial \left( {e_1 \,u\,v} \right)}{\partial j} \right) 534 +\frac{1}{e_3 } \frac{\partial \left( {w\,u } \right)}{\partial k} \\ 535 \qquad \qquad \quad 536 +\frac{1}{e_1 e_2 } \left( 537 -u \left( \frac{\partial \left( {e_1 v } \right)}{\partial j} 538 -v\,\frac{\partial e_1 }{\partial j} \right) 539 -u \frac{\partial \left( {e_2 u } \right)}{\partial i} 540 \right) 541 -\frac{1}{e_3 } \frac{\partial w}{\partial k} u 542 +\frac{1}{e_1 e_2 }\left( -v^2\frac{\partial e_2 }{\partial i} \right) 543 \end{array} 544 } & 545 \end{flalign*} 546 \begin{flalign*} 547 &{ 548 \begin{array}{*{20}l} 549 \qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right) 550 - \left( \nabla \cdot {\rm {\bf U}} \right) \ u 551 +\frac{1}{e_1 e_2 }\left( 552 -v^2 \frac{\partial e_2 }{\partial i} 553 +uv \, \frac{\partial e_1 }{\partial j} \right) \\ 554 \end{array} 555 } & 556 \end{flalign*} 557 as $\nabla \cdot {\rm {\bf U}}\;=0$ (incompressibility) it comes: 558 \begin{flalign*} 559 &{ 560 \begin{array}{*{20}l} 561 \qquad = \nabla \cdot \left( {{\rm {\bf U}}\,u} \right) 562 + \frac{1}{e_1 e_2 } \left( v \; \frac{\partial e_2}{\partial i} 563 -u \; \frac{\partial e_1}{\partial j} \right) \left( -v \right) 564 \end{array} 565 } & 566 \end{flalign*} 567 568 The flux form of the momentum advection term is therefore given by: 569 \begin{multline} 570 \label{eq:PE_flux_form} 571 \left[ 572 \left( {\nabla \times {\rm {\bf U}}} \right) \times {\rm {\bf U}} 573 +\frac{1}{2} \nabla \left( {{\rm {\bf U}}^2} \right) 574 \right]_h \\ 575 = \nabla \cdot \left( {{ 576 \begin{array}{*{20}c} 577 {\rm {\bf U}} \, u \hfill \\ 578 {\rm {\bf U}} \, v \hfill \\ 579 \end{array} 580 }} 581 \right) 582 +\frac{1}{e_1 e_2 } \left( 583 v\frac{\partial e_2}{\partial i} 584 -u\frac{\partial e_1}{\partial j} 585 \right) {\rm {\bf k}} \times {\rm {\bf U}}_h 586 \end{multline} 470 \rt) 471 + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \vect k \times \vect U_h 472 \end{equation} 587 473 588 474 The flux form has two terms, 589 475 the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) 590 476 and the second one is due to the curvilinear nature of the coordinate system used. 591 The latter is called the \ emph{metric} term and can be viewed as a modification of the Coriolis parameter:477 The latter is called the \textit{metric} term and can be viewed as a modification of the Coriolis parameter: 592 478 \[ 593 479 % \label{eq:PE_cor+metric} 594 f \to f + \frac{1}{e_1\;e_2} \left( v \frac{\partial e_2}{\partial i} 595 -u \frac{\partial e_1}{\partial j} \right) 480 f \to f + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) 596 481 \] 597 482 598 483 Note that in the case of geographical coordinate, 599 \ie when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi,a)$,600 we recover the commonly used modification of the Coriolis parameter $f \to f +(u/a) \tan \varphi$.484 \ie when $(i,j) \to (\lambda,\varphi)$ and $(e_1,e_2) \to (a \, \cos \varphi,a)$, 485 we recover the commonly used modification of the Coriolis parameter $f \to f + (u / a) \tan \varphi$. 601 486 602 487 To sum up, the curvilinear $z$-coordinate equations solved by the ocean model can be written in 603 488 the following tensorial formalism: 604 489 605 \ vspace{+10pt}606 $\bullet$ \textbf{Vector invariant form of the momentum equations} : 607 608 \begin{subequations}609 \label{eq:PE_dyn_vect}610 \[490 \begin{itemize} 491 \item 492 \textbf{Vector invariant form of the momentum equations}: 493 \begin{equation} 494 \label{eq:PE_dyn_vect} 495 \begin{split} 611 496 % \label{eq:PE_dyn_vect_u} 612 \begin{split} 613 \frac{\partial u}{\partial t} 614 = + \left( {\zeta +f} \right)\,v 615 - \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left( u^2+v^2 \right) 616 - \frac{1}{e_3 } w \frac{\partial u}{\partial k} & \\ 617 - \frac{1}{e_1 } \frac{\partial}{\partial i} \left( \frac{p_s+p_h }{\rho_o} \right) 618 &+ D_u^{\vect{U}} + F_u^{\vect{U}} \\ \\ 619 \frac{\partial v}{\partial t} = 620 - \left( {\zeta +f} \right)\,u 621 - \frac{1}{2\,e_2 } \frac{\partial }{\partial j}\left( u^2+v^2 \right) 622 - \frac{1}{e_3 } w \frac{\partial v}{\partial k} & \\ 623 - \frac{1}{e_2 } \frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o} \right) 624 &+ D_v^{\vect{U}} + F_v^{\vect{U}} 497 \pd[u]{t} = &+ (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) 498 - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 499 &+ D_u^{\vect U} + F_u^{\vect U} \\ 500 \pd[v]{t} = &- (\zeta + f) \, u - \frac{1}{2 e_2} \pd[]{j} (u^2 + v^2) 501 - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 502 &+ D_v^{\vect U} + F_v^{\vect U} 625 503 \end{split} 626 \] 627 \end{subequations} 628 629 630 \vspace{+10pt} 631 $\bullet$ \textbf{flux form of the momentum equations} : 632 \begin{subequations} 504 \end{equation} 505 \item 506 \textbf{flux form of the momentum equations}: 633 507 % \label{eq:PE_dyn_flux} 634 508 \begin{multline*} 635 509 % \label{eq:PE_dyn_flux_u} 636 \frac{\partial u}{\partial t}= 637 + \left( { f + \frac{1}{e_1 \; e_2} 638 \left( v \frac{\partial e_2}{\partial i} 639 -u \frac{\partial e_1}{\partial j} \right)} \right) \, v \\ 640 - \frac{1}{e_1 \; e_2} \left( 641 \frac{\partial \left( {e_2 \,u\,u} \right)}{\partial i} 642 + \frac{\partial \left( {e_1 \,v\,u} \right)}{\partial j} \right) 643 - \frac{1}{e_3 }\frac{\partial \left( { w\,u} \right)}{\partial k} \\ 644 - \frac{1}{e_1 }\frac{\partial}{\partial i}\left( \frac{p_s+p_h }{\rho_o} \right) 645 + D_u^{\vect{U}} + F_u^{\vect{U}} 510 \pd[u]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v \\ 511 - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\ 512 - \frac{1}{e_3} \pd[(w \, u)]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 513 + D_u^{\vect U} + F_u^{\vect U} 646 514 \end{multline*} 647 515 \begin{multline*} 648 516 % \label{eq:PE_dyn_flux_v} 649 \frac{\partial v}{\partial t}= 650 - \left( { f + \frac{1}{e_1 \; e_2} 651 \left( v \frac{\partial e_2}{\partial i} 652 -u \frac{\partial e_1}{\partial j} \right)} \right) \, u \\ 653 \frac{1}{e_1 \; e_2} \left( 654 \frac{\partial \left( {e_2 \,u\,v} \right)}{\partial i} 655 + \frac{\partial \left( {e_1 \,v\,v} \right)}{\partial j} \right) 656 - \frac{1}{e_3 } \frac{\partial \left( { w\,v} \right)}{\partial k} \\ 657 - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o} \right) 658 + D_v^{\vect{U}} + F_v^{\vect{U}} 517 \pd[v]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\ 518 + \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\ 519 - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 520 + D_v^{\vect U} + F_v^{\vect U} 659 521 \end{multline*} 660 \end{subequations} 661 where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and 662 $p_s $, the surface pressure, is given by: 663 \[ 522 where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and $p_s$, the surface pressure, 523 is given by: 524 \[ 664 525 % \label{eq:PE_spg} 665 p_s = \rho \,g \,\eta666 \]667 with $\eta$ is solution of \autoref{eq:PE_ssh}.668 669 The vertical velocity and the hydrostatic pressure are diagnosed from the following equations:670 \[526 p_s = \rho \,g \, \eta 527 \] 528 with $\eta$ is solution of \autoref{eq:PE_ssh}. 529 530 The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: 531 \[ 671 532 % \label{eq:w_diag} 672 \frac{\partial w}{\partial k}=-\chi \;e_3 673 \] 674 \[ 533 \pd[w]{k} = - \chi \; e_3 \qquad 675 534 % \label{eq:hp_diag} 676 \frac{\partial p_h }{\partial k}=-\rho \;g\;e_3 677 \] 678 where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}. 679 680 \vspace{+10pt} 681 $\bullet$ \textit{tracer equations} : 682 \[ 683 % \label{eq:S} 684 \frac{\partial T}{\partial t} = 685 -\frac{1}{e_1 e_2 }\left[ { \frac{\partial \left( {e_2 T\,u} \right)}{\partial i} 686 +\frac{\partial \left( {e_1 T\,v} \right)}{\partial j}} \right] 687 -\frac{1}{e_3 }\frac{\partial \left( {T\,w} \right)}{\partial k} + D^T + F^T 688 \] 689 \[ 690 % \label{eq:T} 691 \frac{\partial S}{\partial t} = 692 -\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 S\,u} \right)}{\partial i} 693 +\frac{\partial \left( {e_1 S\,v} \right)}{\partial j}} \right] 694 -\frac{1}{e_3 }\frac{\partial \left( {S\,w} \right)}{\partial k} + D^S + F^S 695 \] 696 \[ 697 % \label{eq:rho} 698 \rho =\rho \left( {T,S,z(k)} \right) 699 \] 700 701 The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 535 \pd[p_h]{k} = - \rho \; g \; e_3 536 \] 537 where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}. 538 \item \textit{tracer equations}: 539 \[ 540 %\label{eq:S} 541 \pd[T]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] 542 - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 543 %\label{eq:T} 544 \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] 545 - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S 546 %\label{eq:rho} 547 \rho = \rho \big( T,S,z(k) \big) 548 \] 549 \end{itemize} 550 551 The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 702 552 It will be defined in \autoref{eq:PE_zdf}. 703 The nature and formulation of $ {\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms,553 The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms, 704 554 are discussed in \autoref{chap:SBC}. 705 555 706 707 \newpage 556 \newpage 708 557 709 558 % ================================================================ … … 717 566 Second the ocean floor depends on the geographical position, 718 567 varying from more than 6,000 meters in abyssal trenches to zero at the coast. 719 Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 568 Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 720 569 Therefore, in order to represent the ocean with respect to 721 570 the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height … … 736 585 \begin{equation} 737 586 \label{eq:PE_s} 738 s =s(i,j,k,t)587 s = s(i,j,k,t) 739 588 \end{equation} 740 589 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, … … 759 608 760 609 %\gmcomment{ 761 762 610 %A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient... 763 764 611 the generalized vertical coordinates used in ocean modelling are not orthogonal, 765 612 which contrasts with many other applications in mathematical physics. … … 772 619 The key motivation for maintaining the same horizontal velocity component is that 773 620 the hydrostatic and geostrophic balances are dominant in the large-scale ocean. 774 Use of an alternative quasi -horizontal velocity, for example one oriented parallel to the generalized surface,621 Use of an alternative quasi -horizontal velocity, for example one oriented parallel to the generalized surface, 775 622 would lead to unacceptable numerical errors. 776 Correspondingly, the vertical direction is anti -parallel to the gravitational force in623 Correspondingly, the vertical direction is anti -parallel to the gravitational force in 777 624 all of the coordinate systems. 778 We do not choose the alternative of a quasi -vertical direction oriented normal to779 the surface of a constant generalized vertical coordinate. 625 We do not choose the alternative of a quasi -vertical direction oriented normal to 626 the surface of a constant generalized vertical coordinate. 780 627 781 628 It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between … … 786 633 volume or mass conservation. 787 634 In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about 788 the physical processes producing a flux across the layer interfaces. 789 635 the physical processes producing a flux across the layer interfaces. 790 636 791 637 In this section we first establish the PE in the generalised vertical $s$-coordinate, 792 then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde. 638 then we discuss the particular cases available in \NEMO, namely $z$, \zstar, $s$, and \ztilde. 793 639 %} 794 640 … … 796 642 % The s-coordinate Formulation 797 643 % ------------------------------------------------------------------------------------------------------------- 798 \subsection{\textit{S -}coordinate formulation}799 800 Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k =z$ and thus $e_3=1$,801 we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$,802 which includes $z$-, \zstar- and $\sigma -$coordinates as special cases803 ($s =z$, $s=\zstar$, and $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.).644 \subsection{\textit{S}-coordinate formulation} 645 646 Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k = z$ and 647 thus $e_3 = 1$, we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$, 648 which includes $z$-, \zstar- and $\sigma$-coordinates as special cases 649 ($s = z$, $s = \zstar$, and $s = \sigma = z / H$ or $ = z / \lt( H + \eta \rt)$, resp.). 804 650 A formal derivation of the transformed equations is given in \autoref{apdx:A}. 805 Let us define the vertical scale factor by $e_3 =\partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ),806 and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by:651 Let us define the vertical scale factor by $e_3 = \partial_s z$ ($e_3$ is now a function of $(i,j,k,t)$ ), 652 and the slopes in the $(i,j)$ directions between $s$- and $z$-surfaces by: 807 653 \begin{equation} 808 654 \label{eq:PE_sco_slope} 809 \sigma_1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s 810 \quad \text{, and } \quad 811 \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s 655 \sigma_1 = \frac{1}{e_1} \; \lt. \pd[z]{i} \rt|_s \quad \text{and} \quad 656 \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt|_s 812 657 \end{equation} 813 We also introduce $\omega$, a dia-surface velocity component, defined as the velocity658 We also introduce $\omega$, a dia-surface velocity component, defined as the velocity 814 659 relative to the moving $s$-surfaces and normal to them: 815 660 \[ 816 661 % \label{eq:PE_sco_w} 817 \omega = w - e_3 \, \frac{\partial s}{\partial t} - \sigma_1 \,u - \sigma_2 \,v \\662 \omega = w - e_3 \, \pd[s]{t} - \sigma_1 \, u - \sigma_2 \, v 818 663 \] 819 664 820 The equations solved by the ocean model \autoref{eq:PE} in $s -$coordinate can be written as follows665 The equations solved by the ocean model \autoref{eq:PE} in $s$-coordinate can be written as follows 821 666 (see \autoref{sec:A_momentum}): 822 667 823 \vspace{0.5cm}824 $\bullet$ Vector invariant form of the momentum equation:825 \begin{multline*}668 \begin{itemize} 669 \item \textbf{Vector invariant form of the momentum equation}: 670 \begin{multline*} 826 671 % \label{eq:PE_sco_u_vector} 827 \frac{\partial u }{\partial t}= 828 + \left( {\zeta +f} \right)\,v 829 - \frac{1}{2\,e_1} \frac{\partial}{\partial i} \left( u^2+v^2 \right) 830 - \frac{1}{e_3} \omega \frac{\partial u}{\partial k} \\ 831 - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho_o} \right) 832 + g\frac{\rho }{\rho_o}\sigma_1 833 + D_u^{\vect{U}} + F_u^{\vect{U}} \quad 834 \end{multline*} 835 \begin{multline*} 672 \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 \, e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[u]{k} \\ 673 - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_1 674 + D_u^{\vect U} + F_u^{\vect U} 675 \end{multline*} 676 \begin{multline*} 836 677 % \label{eq:PE_sco_v_vector} 837 \frac{\partial v }{\partial t}= 838 - \left( {\zeta +f} \right)\,u 839 - \frac{1}{2\,e_2 }\frac{\partial }{\partial j}\left( u^2+v^2 \right) 840 - \frac{1}{e_3 } \omega \frac{\partial v}{\partial k} \\ 841 - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o} \right) 842 + g\frac{\rho }{\rho_o }\sigma_2 843 + D_v^{\vect{U}} + F_v^{\vect{U}} \quad 844 \end{multline*} 845 846 \vspace{0.5cm} 847 $\bullet$ Flux form of the momentum equation : 848 \begin{multline*} 678 \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 \, e_2} \pd[]{j}(u^2 + v^2) - \frac{1}{e_3} \omega \pd[v]{k} \\ 679 - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_2 680 + D_v^{\vect U} + F_v^{\vect U} 681 \end{multline*} 682 \item \textbf{Flux form of the momentum equation}: 683 \begin{multline*} 849 684 % \label{eq:PE_sco_u_flux} 850 \frac{1}{e_3} \frac{\partial \left( e_3\,u \right) }{\partial t}= 851 + \left( { f + \frac{1}{e_1 \; e_2 } 852 \left( v \frac{\partial e_2}{\partial i} 853 -u \frac{\partial e_1}{\partial j} \right)} \right) \, v \\ 854 - \frac{1}{e_1 \; e_2 \; e_3 } \left( 855 \frac{\partial \left( {e_2 \, e_3 \, u\,u} \right)}{\partial i} 856 + \frac{\partial \left( {e_1 \, e_3 \, v\,u} \right)}{\partial j} \right) 857 - \frac{1}{e_3 }\frac{\partial \left( { \omega\,u} \right)}{\partial k} \\ 858 - \frac{1}{e_1} \frac{\partial}{\partial i} \left( \frac{p_s + p_h}{\rho_o} \right) 859 + g\frac{\rho }{\rho_o}\sigma_1 860 + D_u^{\vect{U}} + F_u^{\vect{U}} \quad 861 \end{multline*} 862 \begin{multline*} 685 \frac{1}{e_3} \pd[(e_3 \, u)]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v \\ 686 - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) \\ 687 - \frac{1}{e_3} \pd[(\omega \, u)]{k} 688 - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 689 + g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 690 \end{multline*} 691 \begin{multline*} 863 692 % \label{eq:PE_sco_v_flux} 864 \frac{1}{e_3} \frac{\partial \left( e_3\,v \right) }{\partial t}= 865 - \left( { f + \frac{1}{e_1 \; e_2} 866 \left( v \frac{\partial e_2}{\partial i} 867 -u \frac{\partial e_1}{\partial j} \right)} \right) \, u \\ 868 - \frac{1}{e_1 \; e_2 \; e_3 } \left( 869 \frac{\partial \left( {e_2 \; e_3 \,u\,v} \right)}{\partial i} 870 + \frac{\partial \left( {e_1 \; e_3 \,v\,v} \right)}{\partial j} \right) 871 - \frac{1}{e_3 } \frac{\partial \left( { \omega\,v} \right)}{\partial k} \\ 872 - \frac{1}{e_2 }\frac{\partial }{\partial j}\left( \frac{p_s+p_h }{\rho_o} \right) 873 + g\frac{\rho }{\rho_o }\sigma_2 874 + D_v^{\vect{U}} + F_v^{\vect{U}} \quad 875 \end{multline*} 876 877 where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, 878 and the hydrostatic pressure have the same expressions as in $z$-coordinates although 879 they do not represent exactly the same quantities. 880 $\omega$ is provided by the continuity equation (see \autoref{apdx:A}): 881 \[ 693 \frac{1}{e_3} \pd[(e_3 \, v)]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\ 694 - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) \\ 695 - \frac{1}{e_3} \pd[(\omega \, v)]{k} 696 - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 697 + g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 698 \end{multline*} 699 where the relative vorticity, $\zeta$, the surface pressure gradient, 700 and the hydrostatic pressure have the same expressions as in $z$-coordinates although 701 they do not represent exactly the same quantities. 702 $\omega$ is provided by the continuity equation (see \autoref{apdx:A}): 703 \[ 882 704 % \label{eq:PE_sco_continuity} 883 \frac{\partial e_3}{\partial t} + e_3 \; \chi + \frac{\partial \omega }{\partial s} = 0 884 \qquad \text{with }\;\; 885 \chi =\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3 \,u} 886 \right)}{\partial i}+\frac{\partial \left( {e_1 e_3 \,v} \right)}{\partial 887 j}} \right] 888 \] 889 890 \vspace{0.5cm} 891 $\bullet$ tracer equations: 892 \begin{multline*} 705 \pd[e_3]{t} + e_3 \; \chi + \pd[\omega]{s} = 0 \quad \text{with} \quad 706 \chi = \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u)]{i} + \pd[(e_1 e_3 \, v)]{j} \rt) 707 \] 708 \item \textit{tracer equations}: 709 \begin{multline*} 893 710 % \label{eq:PE_sco_t} 894 \frac{1}{e_3} \frac{\partial \left( e_3\,T \right) }{\partial t}= 895 -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,T} \right)}{\partial i} 896 +\frac{\partial \left( {e_1 e_3\,v\,T} \right)}{\partial j}} \right] \\ 897 -\frac{1}{e_3 }\frac{\partial \left( {T\,\omega } \right)}{\partial k} + D^T + F^S \qquad 898 \end{multline*} 899 900 \begin{multline*} 711 \frac{1}{e_3} \pd[(e_3 \, T)]{t} = - \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u \, T)]{i} 712 + \pd[(e_1 e_3 \, v \, T)]{j} \rt) \\ 713 - \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S 714 \end{multline*} 715 \begin{multline} 901 716 % \label{eq:PE_sco_s} 902 \frac{1}{e_3} \frac{\partial \left( e_3\,S \right) }{\partial t}= 903 -\frac{1}{e_1 e_2 e_3 }\left[ {\frac{\partial \left( {e_2 e_3\,u\,S} \right)}{\partial i} 904 +\frac{\partial \left( {e_1 e_3\,v\,S} \right)}{\partial j}} \right] \\ 905 -\frac{1}{e_3 }\frac{\partial \left( {S\,\omega } \right)}{\partial k} + D^S + F^S \qquad 906 \end{multline*} 907 717 \frac{1}{e_3} \pd[(e_3 \, S)]{t} = - \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u \, S)]{i} 718 + \pd[(e_1 e_3 \, v \, S)]{j} \rt) \\ 719 - \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S 720 \end{multline} 721 \end{itemize} 908 722 The equation of state has the same expression as in $z$-coordinate, 909 723 and similar expressions are used for mixing and forcing terms. 910 724 911 725 \gmcomment{ 912 \colorbox{yellow}{ to be updated $= = >$} 913 Add a few works on z and zps and s and underlies the differences between all of them 914 \colorbox{yellow}{ $< = =$ end update} } 915 916 726 \colorbox{yellow}{ to be updated $= = >$} 727 Add a few works on z and zps and s and underlies the differences between all of them 728 \colorbox{yellow}{$< = =$ end update} 729 } 917 730 918 731 % ------------------------------------------------------------------------------------------------------------- 919 732 % Curvilinear \zstar-coordinate System 920 733 % ------------------------------------------------------------------------------------------------------------- 921 \subsection{Curvilinear \zstar- -coordinate system}734 \subsection{Curvilinear \zstar-coordinate system} 922 735 \label{subsec:PE_zco_star} 923 736 … … 925 738 \begin{figure}[!b] 926 739 \begin{center} 927 \includegraphics[width=1.0\textwidth]{Fig_z_zstar} 928 \caption{ \protect\label{fig:z_zstar} 740 \includegraphics[]{Fig_z_zstar} 741 \caption{ 742 \protect\label{fig:z_zstar} 929 743 (a) $z$-coordinate in linear free-surface case ; 930 (b) $z -$coordinate in non-linear free surface case ;931 (c) re-scaled height coordinate (become popular as the \zstar-coordinate932 \citep{Adcroft_Campin_OM04}).744 (b) $z$-coordinate in non-linear free surface case ; 745 (c) re-scaled height coordinate 746 (become popular as the \zstar-coordinate \citep{Adcroft_Campin_OM04}). 933 747 } 934 748 \end{center} … … 936 750 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 937 751 938 939 752 In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. 940 These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site. 941 942 %\gmcomment{ 753 These coordinates systems is presented in a report \citep{Levier2007} available on the \NEMO web site. 754 943 755 The \zstar coordinate approach is an unapproximated, non-linear free surface implementation which allows one to 944 756 deal with large amplitude free-surface variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. … … 947 759 as in the $z$-coordinate formulation, but is equally distributed over the full water column. 948 760 Thus vertical levels naturally follow sea-surface variations, with a linear attenuation with depth, 949 as illustrated by figure fig.1c.950 Note that with a flat bottom, such as in fig.1c, the bottom-following $z$ coordinate and \zstar are equivalent.951 The definition and modified oceanic equations for the rescaled vertical coordinate 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. 763 The definition and modified oceanic equations for the rescaled vertical coordinate \zstar, 952 764 including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004). 953 765 The major points are summarized here. 954 The position ( 766 The position (\zstar) and vertical discretization (\zstar) are expressed as: 955 767 \[ 956 768 % \label{eq:z-star} 957 H + \zstar = (H + z) / r \quad \text{and} \ \delta \zstar = \delta z / r \quad \text{with} \ r = \frac{H+\eta} {H} 958 \] 769 H + \zstar = (H + z) / r \quad \text{and} \quad \delta \zstar 770 = \delta z / r \quad \text{with} \quad r 771 = \frac{H + \eta}{H} 772 \] 959 773 Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar, 960 774 the upper and lower boundaries are at fixed \zstar position, 961 $\zstar = 0$ and 775 $\zstar = 0$ and $\zstar = -H$ respectively. 962 776 Also the divergence of the flow field is no longer zero as shown by the continuity equation: 963 \[ 964 \frac{\partial r}{\partial t} = \nabla_{\zstar} \cdot \left( r \; \rm{\bf U}_h \right) 965 \left( r \; w\textit{*} \right) = 0 966 \] 967 %} 968 777 \[ 778 \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) (r \; w *) = 0 779 \] 969 780 970 781 % from MOM4p1 documentation 971 972 782 To overcome problems with vanishing surface and/or bottom cells, we consider the zstar coordinate 973 783 \[ 974 784 % \label{eq:PE_} 975 z^\star = H \left( \frac{z-\eta}{H+\eta} \right)785 \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) 976 786 \] 977 787 … … 981 791 and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling. 982 792 983 The surfaces of constant $z^\star$ are quasi-horizontal.984 Indeed, the $z^\star$coordinate reduces to $z$ when $\eta$ is zero.793 The surfaces of constant \zstar are quasi -horizontal. 794 Indeed, the \zstar coordinate reduces to $z$ when $\eta$ is zero. 985 795 In general, when noting the large differences between 986 796 undulations of the bottom topography versus undulations in the surface height, 987 it is clear that surfaces constant $z^\star$are very similar to the depth surfaces.797 it is clear that surfaces constant \zstar are very similar to the depth surfaces. 988 798 These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to 989 799 terrain following sigma models discussed in \autoref{subsec:PE_sco}. 990 Additionally, since $z^\star$when $\eta = 0$,800 Additionally, since \zstar when $\eta = 0$, 991 801 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography. 992 802 This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of 993 803 nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, 994 804 depending on the sophistication of the pressure gradient solver. 995 The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of996 neutral physics parameterizations in $z^\star$models using the same techniques as in $z$-models805 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$-models 997 807 (see Chapters 13-16 of \cite{Griffies_Bk04}) for a discussion of neutral physics in $z$-models, 998 808 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). 999 809 1000 The range over which $z^\star$ varies is time independent $-H \leq z^\star \leq 0$.810 The range over which \zstar varies is time independent $-H \leq \zstar \leq 0$. 1001 811 Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$. 1002 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 1003 1004 Because $z^\star$has a time independent range, all grid cells have static increments ds,1005 and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H. 1006 The $z^\star$coordinate is therefore invisible to undulations of the free surface,812 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 813 814 Because \zstar has a time independent range, all grid cells have static increments ds, 815 and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H. 816 The \zstar coordinate is therefore invisible to undulations of the free surface, 1007 817 since it moves along with the free surface. 1008 This proper ty means that no spurious vertical transport is induced across surfaces of constant $z^\star$by818 This proper ty means that no spurious vertical transport is induced across surfaces of constant \zstar by 1009 819 the motion of external gravity waves. 1010 820 Such spurious transpor t can be a problem in z-models, especially those with tidal forcing. 1011 Quite generally, the time independent range for the $z^\star$coordinate is a very convenient property that821 Quite generally, the time independent range for the \zstar coordinate is a very convenient property that 1012 822 allows for a nearly arbitrary ver tical resolution even in the presence of large amplitude fluctuations of 1013 823 the surface height, again so long as $\eta > -H$. 1014 1015 824 %end MOM doc %%% 1016 825 1017 1018 1019 \newpage 826 \newpage 1020 827 1021 828 % ------------------------------------------------------------------------------------------------------------- … … 1035 842 For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. 1036 843 Topographic Rossby waves can be excited and can interact with the mean current. 1037 In the $z -$coordinate system presented in the previous section (\autoref{sec:PE_zco}),1038 $z -$surfaces are geopotential surfaces.844 In the $z$-coordinate system presented in the previous section (\autoref{sec:PE_zco}), 845 $z$-surfaces are geopotential surfaces. 1039 846 The bottom topography is discretised by steps. 1040 847 This often leads to a misrepresentation of a gradually sloping bottom and to … … 1043 850 One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of 1044 851 a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 1045 Another solution is to introduce a terrain-following coordinate system (hereafter $s -$coordinate).852 Another solution is to introduce a terrain-following coordinate system (hereafter $s$-coordinate). 1046 853 1047 854 The $s$-coordinate avoids the discretisation error in the depth field since the layers of … … 1050 857 which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths, 1051 858 can easily be represented (with relatively low vertical resolution). 1052 A terrain-following model (hereafter $s -$model) also facilitates the modelling of the boundary layer flows over859 A terrain-following model (hereafter $s$-model) also facilitates the modelling of the boundary layer flows over 1053 860 a large depth range, which in the framework of the $z$-model would require high vertical resolution over 1054 861 the whole depth range. … … 1063 870 \begin{equation} 1064 871 \label{eq:PE_p_sco} 1065 \left. {\nabla p} \right|_z =\left. {\nabla p} \right|_s -\frac{\partial 1066 p}{\partial s}\left. {\nabla z} \right|_s 872 \nabla p |_z = \nabla p |_s - \pd[p]{s} \nabla z |_s 1067 873 \end{equation} 1068 874 1069 875 The second term in \autoref{eq:PE_p_sco} depends on the tilt of the coordinate surface and 1070 876 introduces a truncation error that is not present in a $z$-model. 1071 In the special case of a $\sigma -$coordinate (\iedepth-normalised coordinate system $\sigma = z/H$),877 In the special case of a $\sigma$-coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 1072 878 \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error. 1073 879 It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, … … 1085 891 (see \autoref{sec:DOM_zgr}. 1086 892 1087 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of any finite difference model. 893 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of 894 any finite difference model. 1088 895 It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. 1089 896 This is the case for a $z$-model as well as for a $s$-model. 1090 However, density varies more strongly on $s -$surfaces than on horizontal surfaces in regions of897 However, density varies more strongly on $s$-surfaces than on horizontal surfaces in regions of 1091 898 large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model. 1092 899 Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus … … 1101 908 (see \autoref{subsec:PE_ldf}). 1102 909 Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 1103 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 1104 1105 The $s -$coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from910 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 911 912 The $s$-coordinates introduced here \citep{Lott_al_OM90,Madec_al_JPO96} differ mainly in two aspects from 1106 913 similar models: 1107 914 it allows a representation of bottom topography with mixed full or partial step-like/terrain following topography; 1108 915 It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. 1109 916 1110 1111 \newpage1112 1113 917 % ------------------------------------------------------------------------------------------------------------- 1114 918 % Curvilinear z-tilde coordinate System 1115 919 % ------------------------------------------------------------------------------------------------------------- 1116 \subsection{\texorpdfstring{Curvilinear \ztilde- -coordinate}{}}920 \subsection{\texorpdfstring{Curvilinear \ztilde-coordinate}{}} 1117 921 \label{subsec:PE_zco_tilde} 1118 922 1119 The \ztilde -coordinate has been developed by \citet{Leclair_Madec_OM11}.923 The \ztilde -coordinate has been developed by \citet{Leclair_Madec_OM11}. 1120 924 It is available in \NEMO since the version 3.4. 1121 925 Nevertheless, it is currently not robust enough to be used in all possible configurations. 1122 926 Its use is therefore not recommended. 1123 1124 927 1125 928 \newpage … … 1134 937 a few kilometres in the horizontal, a few meters in the vertical and a few minutes. 1135 938 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. 1136 The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) must be represented entirely in terms of large-scale patterns to close the equations. 939 The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) 940 must be represented entirely in terms of large-scale patterns to close the equations. 1137 941 These effects appear in the equations as the divergence of turbulent fluxes 1138 942 (\ie fluxes associated with the mean correlation of small scale perturbations). … … 1144 948 1145 949 The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions. 1146 Therefore subgrid-scale physics \textbf{D}$^{\vect {U}}$, $D^{S}$ and $D^{T}$ in950 Therefore subgrid-scale physics \textbf{D}$^{\vect U}$, $D^{S}$ and $D^{T}$ in 1147 951 \autoref{eq:PE_dyn}, \autoref{eq:PE_tra_T} and \autoref{eq:PE_tra_S} are divided into 1148 a lateral part \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and1149 a vertical part \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$.952 a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and 953 a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$. 1150 954 The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 1151 955 … … 1160 964 Turbulent motions are thus never explicitly solved, even partially, but always parameterized. 1161 965 The vertical turbulent fluxes are assumed to depend linearly on the gradients of large-scale quantities 1162 (for example, the turbulent heat flux is given by $\overline{T' w'}=-A^{vT} \partial_z \overline T$,1163 where $A^{v T}$ is an eddy coefficient).966 (for example, the turbulent heat flux is given by $\overline{T' w'} = -A^{v T} \partial_z \overline T$, 967 where $A^{v T}$ is an eddy coefficient). 1164 968 This formulation is analogous to that of molecular diffusion and dissipation. 1165 969 This is quite clearly a necessary compromise: considering only the molecular viscosity acting on … … 1169 973 \begin{equation} 1170 974 \label{eq:PE_zdf} 1171 \begin{split} 1172 {\vect{D}}^{v \vect{U}} &=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\ 1173 D^{vT} &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ , 1174 \quad 1175 D^{vS}=\frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial S}{\partial z}} \right) 1176 \end{split} 975 \begin{gathered} 976 \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \ , \\ 977 D^{vT} = \pd[]{z} \lt( A^{vT} \pd[T]{z} \rt) \quad \text{and} \quad 978 D^{vS} = \pd[]{z} \lt( A^{vT} \pd[S]{z} \rt) 979 \end{gathered} 1177 980 \end{equation} 1178 981 where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. … … 1232 1035 one uses a advective scheme which is diffusive enough to maintain the model stability. 1233 1036 It must be emphasised that then, all the sub-grid scale physics is included in the formulation of 1234 the advection scheme. 1037 the advection scheme. 1235 1038 1236 1039 All these parameterisations of subgrid scale physics have advantages and drawbacks. … … 1246 1049 \begin{equation} 1247 1050 \label{eq:PE_iso_tensor} 1248 D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 1249 \mbox{with}\quad \;\;\Re =\left( {{ 1250 \begin{array}{*{20}c} 1251 1 \hfill & 0 \hfill & {-r_1 } \hfill \\ 1252 0 \hfill & 1 \hfill & {-r_2 } \hfill \\ 1253 {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 1254 \end{array} 1255 }} \right) 1051 D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad 1052 \Re = 1053 \begin{pmatrix} 1054 1 & 0 & -r_1 \\ 1055 0 & 1 & -r_2 \\ 1056 -r_1 & -r_2 & r_1^2 + r_2^2 \\ 1057 \end{pmatrix} 1256 1058 \end{equation} 1257 where $r_1 \;\mbox{and}\;r_2$ are the slopes between the surface along which the diffusive operator acts and1258 the model level ( $e. g.$$z$- or $s$-surfaces).1059 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). 1259 1061 Note that the formulation \autoref{eq:PE_iso_tensor} is exact for 1260 1062 the rotation between geopotential and $s$-surfaces, … … 1267 1069 1268 1070 For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. 1269 $\Re $ reduces to the identity in the horizontal direction, no rotation is applied.1071 $\Re$ reduces to the identity in the horizontal direction, no rotation is applied. 1270 1072 1271 1073 For \textit{geopotential} diffusion, … … 1278 1080 \begin{equation} 1279 1081 \label{eq:PE_iso_slopes} 1280 r_1 = \frac{e_3 }{e_1 } \left( \pd[\rho]{i} \right) \left( \pd[\rho]{k} \right)^{-1} \,\quad1281 r_2 = \frac{e_3 }{e_2 } \left( \pd[\rho]{j} \right) \left( \pd[\rho]{k} \right)^{-1} \,1082 r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 1083 r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 1282 1084 \end{equation} 1283 1085 while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. 1284 1086 1285 1087 \subsubsection{Eddy induced velocity} 1088 1286 1089 When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used, 1287 1090 an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: 1288 1091 \[ 1289 1092 % \label{eq:PE_iso+eiv} 1290 D^{lT}=\nabla \cdot \left( {A^{lT}\;\Re \;\nabla T} \right) 1291 +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right) 1093 D^{lT} = \nabla \cdot \lt( A^{lT} \; \Re \; \nabla T \rt) + \nabla \cdot \lt( \vect U^\ast \, T \rt) 1292 1094 \] 1293 where $ {\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent,1095 where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a non-divergent, 1294 1096 eddy-induced transport velocity. This velocity field is defined by: 1295 1097 \[ 1296 1098 % \label{eq:PE_eiv} 1297 \begin{split} 1298 u^\ast &= +\frac{1}{e_3 }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_1 } \right] \\ 1299 v^\ast &= +\frac{1}{e_3 }\frac{\partial }{\partial k}\left[ {A^{eiv}\;\tilde{r}_2 } \right] \\ 1300 w^\ast &= -\frac{1}{e_1 e_2 }\left[ 1301 \frac{\partial }{\partial i}\left( {A^{eiv}\;e_2\,\tilde{r}_1 } \right) 1302 +\frac{\partial }{\partial j}\left( {A^{eiv}\;e_1\,\tilde{r}_2 } \right) \right] 1303 \end{split} 1099 u^\ast &= \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_1 \rt) \\ 1100 v^\ast &= \frac{1}{e_3} \pd[]{k} \lt( A^{eiv} \; \tilde{r}_2 \rt) \\ 1101 w^\ast &= - \frac{1}{e_1 e_2} \lt[ \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt) 1102 + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt] 1304 1103 \] 1305 1104 where $A^{eiv}$ is the eddy induced velocity coefficient 1306 1105 (or equivalently the isoneutral thickness diffusivity coefficient), 1307 and $\tilde {r}_1$ and $\tilde{r}_2$ are the slopes between isoneutral and \emph{geopotential} surfaces.1106 and $\tilde r_1$ and $\tilde r_2$ are the slopes between isoneutral and \textit{geopotential} surfaces. 1308 1107 Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate: 1309 1108 \begin{align} 1310 1109 \label{eq:PE_slopes_eiv} 1311 1110 \tilde{r}_n = 1312 \begin{cases}1313 r_n & \text{in $z$-coordinate}\\1314 r_n + \sigma_n & \text{in \zstarand $s$-coordinates}1315 \end{cases}1316 \quad \text{where } n=1,21111 \begin{cases} 1112 r_n & \text{in $z$-coordinate} \\ 1113 r_n + \sigma_n & \text{in \zstar- and $s$-coordinates} 1114 \end{cases} 1115 \quad \text{where~} n = 1, 2 1317 1116 \end{align} 1318 1117 1319 1118 The normal component of the eddy induced velocity is zero at all the boundaries. 1320 This can be achieved in a model by tapering either the eddy coefficient or 1321 the slopes to zero in the vicinity of theboundaries.1119 This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of 1120 the boundaries. 1322 1121 The latter strategy is used in \NEMO (cf. \autoref{chap:LDF}). 1323 1122 … … 1327 1126 \[ 1328 1127 % \label{eq:PE_bilapT} 1329 D^{lT}= - \Delta \ left( \;\Delta T \right)1330 \ qquad \text{where} \;\; \Delta \bullet = \nabla \left( {\sqrt{B^{lT}\,}\;\Re \;\nabla \bullet} \right)1128 D^{lT}= - \Delta \; (\Delta T) \quad \text{where} \quad 1129 \Delta \bullet = \nabla \lt( \sqrt{B^{lT}} \; \Re \; \nabla \bullet \rt) 1331 1130 \] 1332 1131 It is the Laplacian operator given by \autoref{eq:PE_iso_tensor} applied twice with 1333 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1334 1132 the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 1335 1133 1336 1134 \subsubsection{Lateral Laplacian momentum diffusive operator} … … 1338 1136 The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by 1339 1137 applying \autoref{eq:PE_lap_vector} to the horizontal velocity vector (see \autoref{apdx:B}): 1340 \ [1138 \begin{align*} 1341 1139 % \label{eq:PE_lapU} 1342 \begin{split} 1343 {\rm {\bf D}}^{l{\rm {\bf U}}} 1344 &= \quad \ \nabla _h \left( {A^{lm}\chi } \right) 1345 \ - \ \nabla _h \times \left( {A^{lm}\,\zeta \;{\rm {\bf k}}} \right) \\ 1346 &= \left( 1347 \begin{aligned} 1348 \frac{1}{e_1 } \frac{\partial \left( A^{lm} \chi \right)}{\partial i} 1349 &-\frac{1}{e_2 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial j} \\ 1350 \frac{1}{e_2 }\frac{\partial \left( {A^{lm} \chi } \right)}{\partial j} 1351 &+\frac{1}{e_1 e_3}\frac{\partial \left( {A^{lm} \;e_3 \zeta} \right)}{\partial i} 1352 \end{aligned} 1353 \right) 1354 \end{split} 1355 \] 1140 \vect D^{l \vect U} &= \nabla_h \big( A^{lm} \chi \big) 1141 - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\ 1142 &= \lt( \frac{1}{e_1} \pd[ \lt( A^{lm} \chi \rt) ]{i} \rt. 1143 - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j} 1144 \frac{1}{e_2} \pd[ \lt( A^{lm} \chi \rt) ]{j} 1145 \lt. + \frac{1}{e_1 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{i} \rt) 1146 \end{align*} 1356 1147 1357 1148 Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields … … 1359 1150 Unfortunately, it is only available in \textit{iso-level} direction. 1360 1151 When a rotation is required 1361 (\ie geopotential diffusion in $s -$coordinates or isoneutral diffusion in both $z$- and $s$-coordinates),1362 the $u$ and $v -$fields are considered as independent scalar fields, so that the diffusive operator is given by:1363 \ [1152 (\ie geopotential diffusion in $s$-coordinates or isoneutral diffusion in both $z$- and $s$-coordinates), 1153 the $u$ and $v$-fields are considered as independent scalar fields, so that the diffusive operator is given by: 1154 \begin{gather*} 1364 1155 % \label{eq:PE_lapU_iso} 1365 \begin{split} 1366 D_u^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla u} \right) \\ 1367 D_v^{l{\rm {\bf U}}} &= \nabla .\left( {A^{lm} \;\Re \;\nabla v} \right) 1368 \end{split} 1369 \] 1156 D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \\ 1157 D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt) 1158 \end{gather*} 1370 1159 where $\Re$ is given by \autoref{eq:PE_iso_tensor}. 1371 1160 It is the same expression as those used for diffusive operator on tracers. 1372 1161 It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, 1373 \ie on a $f -$ or $\beta-$plane, not on the sphere.1162 \ie on a $f$- or $\beta$-plane, not on the sphere. 1374 1163 It is also a very good approximation in vicinity of the Equator in 1375 1164 a geographical coordinate system \citep{Lengaigne_al_JGR03}. … … 1386 1175 1387 1176 \end{document} 1388 -
NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex
r10442 r10501 4 4 5 5 % ================================================================ 6 % Chapter 2 Time Domain (step.F90)6 % Chapter 2 ——— Time Domain (step.F90) 7 7 % ================================================================ 8 \chapter{Time Domain (STP) 8 \chapter{Time Domain (STP)} 9 9 \label{chap:STP} 10 11 10 \minitoc 12 11 … … 35 34 \begin{equation} 36 35 \label{eq:STP} 37 x^{t +\rdt} = x^{t-\rdt} + 2 \, \rdt \ \text{RHS}_x^{t-\rdt,\,t,\,t+\rdt}36 x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t - \rdt, \, t, \, t + \rdt} 38 37 \end{equation} 39 38 where $x$ stands for $u$, $v$, $T$ or $S$; … … 85 84 \begin{equation} 86 85 \label{eq:STP_asselin} 87 x_F^t = x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right]88 \end{equation} 86 x_F^t = x^t + \gamma \, \lt[ x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt] 87 \end{equation} 89 88 where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient. 90 89 $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter). 91 Its default value is \np{rn\_atfp} \forcode{= 10.e-3} (see \autoref{sec:STP_mLF}),90 Its default value is \np{rn\_atfp}~\forcode{= 10.e-3} (see \autoref{sec:STP_mLF}), 92 91 causing only a weak dissipation of high frequency motions (\citep{Farge1987}). 93 92 The addition of a time filter degrades the accuracy of the calculation from second to first order. … … 111 110 (when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : 112 111 \[ 113 % 114 x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \ {D_x}^{t-\rdt}115 \] 112 %\label{eq:STP_euler} 113 x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ D_x^{t - \rdt} 114 \] 116 115 117 116 This is diffusive in time and conditionally stable. … … 119 118 \begin{equation} 120 119 \label{eq:STP_euler_stability} 121 A^h < \left\{ 122 \begin{aligned} 123 &\frac{e^2}{ 8 \; \rdt } &&\quad \text{laplacian diffusion} \\ 124 &\frac{e^4}{64 \; \rdt } &&\quad \text{bilaplacian diffusion} 125 \end{aligned} 126 \right. 120 A^h < 121 \begin{cases} 122 \frac{e^2}{ 8 \, \rdt} & \text{laplacian diffusion} \\ 123 \frac{e^4}{64 \, \rdt} & \text{bilaplacian diffusion} 124 \end{cases} 127 125 \end{equation} 128 126 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. … … 134 132 but usually the numerical stability condition imposes a strong constraint on the time step. 135 133 Two solutions are available in \NEMO to overcome the stability constraint: 136 $(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp} \forcode{= .true.}) or137 $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp} \forcode{= .false.}).138 In $(a)$, the master time step $\Delta 134 $(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp}~\forcode{= .true.}) or 135 $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}~\forcode{= .false.}). 136 In $(a)$, the master time step $\Delta$t is cut into $N$ fractional time steps so that 139 137 the stability criterion is reduced by a factor of $N$. 140 138 The computation is performed as follows: 141 \ [139 \begin{alignat*}{2} 142 140 % \label{eq:STP_ts} 143 \begin{split} 144 & x_\ast ^{t-\rdt} = x^{t-\rdt} \\ 145 & x_\ast ^{t-\rdt+L\frac{2\rdt}{N}}=x_\ast ^{t-\rdt+\left( {L-1} 146 \right)\frac{2\rdt}{N}}+\frac{2\rdt}{N}\;\text{DF}^{t-\rdt+\left( {L-1} \right)\frac{2\rdt}{N}} 147 \quad \text{for $L=1$ to $N$} \\ 148 & x^{t+\rdt} = x_\ast^{t+\rdt} 149 \end{split} 150 \] 141 &x_\ast^{t - \rdt} &= &x^{t - \rdt} \\ 142 &x_\ast^{t - \rdt + L \frac{2 \rdt}{N}} &= &x_\ast ^{t - \rdt + (L - 1) \frac{2 \rdt}{N}} 143 + \frac{2 \rdt}{N} \; DF^{t - \rdt + (L - 1) \frac{2 \rdt}{N}} 144 \quad \text{for $L = 1$ to $N$} \\ 145 &x^{t + \rdt} &= &x_\ast^{t + \rdt} 146 \end{alignat*} 151 147 with DF a vertical diffusion term. 152 148 The number of fractional time steps, $N$, is given by setting \np{nn\_zdfexp}, (namelist parameter). … … 154 150 \begin{equation} 155 151 \label{eq:STP_imp} 156 x^{t +\rdt} = x^{t-\rdt} + 2 \, \rdt \ \text{RHS}_x^{t+\rdt}157 \end{equation} 152 x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t + \rdt} 153 \end{equation} 158 154 159 155 %%gm … … 166 162 \[ 167 163 % \label{eq:STP_imp_zdf} 168 \frac{T(k)^{t +1}-T(k)^{t-1}}{2\;\rdt}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta169 _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}170 \ right]164 \frac{T(k)^{t + 1} - T(k)^{t - 1}}{2 \; \rdt} 165 \equiv 166 \text{RHS} + \frac{1}{e_{3t}} \delta_k \lt[ \frac{A_w^{vT}}{e_{3w} } \delta_{k + 1/2} \lt[ T^{t + 1} \rt] \rt] 171 167 \] 172 168 where RHS is the right hand side of the equation except for the vertical diffusion term. … … 174 170 \begin{equation} 175 171 \label{eq:STP_imp_mat} 176 -c(k +1)\;T^{t+1}(k+1) + d(k)\;T^{t+1}(k) - \;c(k)\;T^{t+1}(k-1) \equiv b(k)172 -c(k + 1) \; T^{t + 1}(k + 1) + d(k) \; T^{t + 1}(k) - \; c(k) \; T^{t + 1}(k - 1) \equiv b(k) 177 173 \end{equation} 178 174 where 179 \begin{align*} 175 \begin{align*} 180 176 c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k) \\ 181 d(k) &= e_{3t} (k) \, / \, (2\rdt) + c_k + c_{k+1} \\182 b(k) &= e_{3t} (k) \; \left( T^{t-1}(k) \, / \, (2\rdt) + \text{RHS} \right)177 d(k) &= e_{3t} (k) \, / \, (2 \rdt) + c_k + c_{k + 1} \\ 178 b(k) &= e_{3t} (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt) 183 179 \end{align*} 184 180 … … 189 185 (see for example \citet{Richtmyer1967}). 190 186 191 192 193 187 % ------------------------------------------------------------------------------------------------------------- 194 188 % Surface Pressure gradient … … 203 197 \begin{figure}[!t] 204 198 \begin{center} 205 \includegraphics[ width=0.7\textwidth]{Fig_TimeStepping_flowchart}199 \includegraphics[]{Fig_TimeStepping_flowchart} 206 200 \caption{ 207 201 \protect\label{fig:TimeStep_flowchart} 208 202 Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}. 209 The use of a semi -implicit computation of the hydrostatic pressure gradient requires the tracer equation to203 The use of a semi -implicit computation of the hydrostatic pressure gradient requires the tracer equation to 210 204 be stepped forward prior to the momentum equation. 211 205 The need for knowledge of the vertical scale factor (here denoted as $h$) requires the sea surface height and … … 225 219 \label{sec:STP_mLF} 226 220 227 Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 221 Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to 222 ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 228 223 The modifications affect both the forcing and filtering treatments in the LF-RA scheme. 229 224 230 225 In a classical LF-RA environment, the forcing term is centred in time, 231 \ie it is time-stepped over a $2 \rdt$ period:232 $x^t = x^t + 2\rdt Q^t$ where $Q$ is the forcing applied to $x$,226 \ie it is time-stepped over a $2 \rdt$ period: 227 $x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, 233 228 and the time filter is given by \autoref{eq:STP_asselin} so that $Q$ is redistributed over several time step. 234 229 In the modified LF-RA environment, these two formulations have been replaced by: 235 \begin{align} 236 x^{t+\rdt} &= x^{t-\rdt} + \rdt \left( Q^{t-\rdt/2} + Q^{t+\rdt/2} \right) \label{eq:STP_forcing} \\ 237 % 238 x_F^t &= x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right] 239 - \gamma\,\rdt \, \left[ Q^{t+\rdt/2} - Q^{t-\rdt/2} \right] \label{eq:STP_RA} 240 \end{align} 230 \begin{gather} 231 \label{eq:STP_forcing} 232 x^{t + \rdt} = x^{t - \rdt} + \rdt \lt( Q^{t - \rdt / 2} + Q^{t + \rdt / 2} \rt) \\ 233 \label{eq:STP_RA} 234 x_F^t = x^t + \gamma \, \lt( x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt) 235 - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) 236 \end{gather} 241 237 The change in the forcing formulation given by \autoref{eq:STP_forcing} (see \autoref{fig:MLF_forcing}) 242 238 has a significant effect: 243 the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}. 239 the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}. 244 240 % forcing seen by the model.... 245 241 This property improves the LF-RA scheme in two respects. … … 250 246 Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, 251 247 the modified formulation becomes conservative \citep{Leclair_Madec_OM09}. 252 Second, the LF-RA becomes a truly quasi -second order scheme.248 Second, the LF-RA becomes a truly quasi -second order scheme. 253 249 Indeed, \autoref{eq:STP_forcing} used in combination with a careful treatment of static instability 254 250 (\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}), 255 251 the two other main sources of time step divergence, 256 allows a reduction by two orders of magnitude of the Asselin filter parameter. 252 allows a reduction by two orders of magnitude of the Asselin filter parameter. 257 253 258 254 Note that the forcing is now provided at the middle of a time step: 259 $Q^{t +\rdt/2}$ is the forcing applied over the $[t,t+\rdt]$ time interval.255 $Q^{t + \rdt / 2}$ is the forcing applied over the $[t,t + \rdt]$ time interval. 260 256 This and the change in the time filter, \autoref{eq:STP_RA}, 261 257 allows an exact evaluation of the contribution due to the forcing term between any two time steps, … … 265 261 \begin{figure}[!t] 266 262 \begin{center} 267 \includegraphics[ width=0.90\textwidth]{Fig_MLF_forcing}263 \includegraphics[]{Fig_MLF_forcing} 268 264 \caption{ 269 265 \protect\label{fig:MLF_forcing} … … 271 267 (top) ''Traditional'' formulation: 272 268 the forcing is defined at the same time as the variable to which it is applied 273 (integer value of the time step index) and it is applied over a $2 \rdt$ period.269 (integer value of the time step index) and it is applied over a $2 \rdt$ period. 274 270 (bottom) modified formulation: 275 271 the forcing is defined in the middle of the time (integer and a half value of the time step index) and 276 the mean of two successive forcing values ($n -1/2$, $n+1/2$) is applied over a $2\rdt$ period.272 the mean of two successive forcing values ($n - 1 / 2$, $n + 1 / 2$) is applied over a $2 \rdt$ period. 277 273 } 278 274 \end{center} … … 297 293 \] 298 294 This is done simply by keeping the leapfrog environment (\ie the \autoref{eq:STP} three level time stepping) but 299 setting all $x^0$ (\textit{before}) and $x^ {1}$ (\textit{now}) fields equal at the first time step and295 setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and 300 296 using half the value of $\rdt$. 301 297 … … 305 301 running the model for $2N$ time steps in one go, 306 302 or by performing two consecutive experiments of $N$ steps with a restart. 307 This requires saving two time levels and many auxiliary data in the restart files in machine precision. 308 309 Note that when a semi -implicit scheme is used to evaluate the hydrostatic pressure gradient303 This requires saving two time levels and many auxiliary data in the restart files in machine precision. 304 305 Note that when a semi -implicit scheme is used to evaluate the hydrostatic pressure gradient 310 306 (see \autoref{subsec:DYN_hpg_imp}), an extra three-dimensional field has to 311 307 be added to the restart file to ensure an exact restartability. 312 This is done optionally via the \np{nn\_dynhpg\_rst} namelist parameter,308 This is done optionally via the \np{nn\_dynhpg\_rst} namelist parameter, 313 309 so that the size of the restart file can be reduced when restartability is not a key issue 314 310 (operational oceanography or in ensemble simulations for seasonal forecasting). 315 311 316 312 Note the size of the time step used, $\rdt$, is also saved in the restart file. 317 When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed. 318 Options are defined through the \ngn{namrun} namelist variables.313 When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed. 314 Options are defined through the \ngn{namrun} namelist variables. 319 315 %%% 320 316 \gmcomment{ … … 360 356 361 357 \begin{flalign*} 362 &\frac{\l eft( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\rdt}363 \equiv \text{RHS}+ \delta_k \l eft[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}364 \r ight] \\365 &\l eft( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}366 \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \l eft[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}367 \r ight] \\368 &\l eft( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}358 &\frac{\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}}{2\rdt} 359 \equiv \text{RHS}+ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]} 360 \rt] \\ 361 &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1} 362 \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]} 363 \rt] \\ 364 &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1} 369 365 \equiv 2\rdt \ \text{RHS} 370 + 2\rdt \ \l eft\{ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} [ T_{k+1}^{t+1} - T_k ^{t+1} ]371 - \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k ^{t+1} - T_{k-1}^{t+1} ] \right\} \\366 + 2\rdt \ \lt\{ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} [ T_{k +1}^{t+1} - T_k ^{t+1} ] 367 - \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} [ T_k ^{t+1} - T_{k -1}^{t+1} ] \rt\} \\ 372 368 &\\ 373 &\l eft( e_{3t}\,T \right)_k^{t+1}374 - {2\rdt} \ \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} T_{k+1}^{t+1}375 + {2\rdt} \ \l eft\{ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}376 + \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \right\} T_{k }^{t+1}377 - {2\rdt} \ \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} T_{k-1}^{t+1} \\378 &\equiv \l eft( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS} \\379 % 369 &\lt( e_{3t}\,T \rt)_k^{t+1} 370 - {2\rdt} \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} T_{k +1}^{t+1} 371 + {2\rdt} \ \lt\{ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} 372 + \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \rt\} T_{k }^{t+1} 373 - {2\rdt} \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} T_{k -1}^{t+1} \\ 374 &\equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS} \\ 375 % 380 376 \end{flalign*} 381 382 377 \begin{flalign*} 383 378 \allowdisplaybreaks 384 379 \intertext{ Tracer case } 385 380 % 386 & \qquad \qquad \quad - {2\rdt} \ \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}387 \qquad \qquad \qquad \qquad T_{k +1}^{t+1} \\388 &+ {2\rdt} \ \biggl\{ (e_{3t})_{k }^{t+1} \bigg. + \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}389 + \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\} \ \ \ T_{k }^{t+1} &&\\390 & \qquad \qquad \qquad \qquad \qquad \quad \ \ - {2\rdt} \ \l eft[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \quad \ \ T_{k-1}^{t+1}391 \ \equiv \ \l eft( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS} \\381 & \qquad \qquad \quad - {2\rdt} \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} 382 \qquad \qquad \qquad \qquad T_{k +1}^{t+1} \\ 383 &+ {2\rdt} \ \biggl\{ (e_{3t})_{k }^{t+1} \bigg. + \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} 384 + \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \bigg. \biggr\} \ \ \ T_{k }^{t+1} &&\\ 385 & \qquad \qquad \qquad \qquad \qquad \quad \ \ - {2\rdt} \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \quad \ \ T_{k -1}^{t+1} 386 \ \equiv \ \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS} \\ 392 387 % 393 388 \end{flalign*} … … 396 391 \intertext{ Tracer content case } 397 392 % 398 & - {2\rdt} \ & \frac{(A_w^{vt})_{k +1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_{k+1}^{t+1}} && \ \left( e_{3t}\,T \right)_{k+1}^{t+1} &\\399 & + {2\rdt} \ \l eft[ 1 \right.+ & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_k^{t+1}}400 + & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_k^{t+1}} \left. \right] & \left( e_{3t}\,T \right)_{k }^{t+1} &\\401 & - {2\rdt} \ & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_{k-1}^{t+1}} &\ \left( e_{3t}\,T \right)_{k-1}^{t+1}402 \equiv \l eft( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS} &393 & - {2\rdt} \ & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_{k +1}^{t+1}} && \ \lt( e_{3t}\,T \rt)_{k +1}^{t+1} &\\ 394 & + {2\rdt} \ \lt[ 1 \rt.+ & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_k^{t+1}} 395 + & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_k^{t+1}} \lt. \rt] & \lt( e_{3t}\,T \rt)_{k }^{t+1} &\\ 396 & - {2\rdt} \ & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_{k -1}^{t+1}} &\ \lt( e_{3t}\,T \rt)_{k -1}^{t+1} 397 \equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS} & 403 398 \end{flalign*} 404 399 405 400 %% 406 401 } 407 %% 402 408 403 \biblio 409 404
Note: See TracChangeset
for help on using the changeset viewer.