Changeset 9407 for branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex
- Timestamp:
- 2018-03-15T17:40:35+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/DOC/tex_sub/chap_DYN.tex
r9394 r9407 5 5 % ================================================================ 6 6 \chapter{Ocean Dynamics (DYN)} 7 \label{ DYN}7 \label{chap:DYN} 8 8 \minitoc 9 9 … … 11 11 $\ $\newline %force an empty line 12 12 13 Using the representation described in Chapter \ref{DOM}, several semi-discrete13 Using the representation described in \autoref{chap:DOM}, several semi-discrete 14 14 space forms of the dynamical equations are available depending on the vertical 15 15 coordinate used and on the conservation properties of the vorticity term. In all … … 36 36 inputs (surface wind stress calculation using bulk formulae, estimation of mixing 37 37 coefficients) that are carried out in modules SBC, LDF and ZDF and are described 38 in Chapters \ref{SBC}, \ref{LDF} and \ref{ZDF}, respectively.38 in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 39 39 40 40 In the present chapter we also describe the diagnostic equations used to compute … … 51 51 The user has the option of extracting and outputting each tendency term from the 52 52 3D momentum equations (\key{trddyn} defined), as described in 53 Chap.\ref{MISC}. Furthermore, the tendency terms associated with the 2D53 \autoref{chap:MISC}. Furthermore, the tendency terms associated with the 2D 54 54 barotropic vorticity balance (when \key{trdvor} is defined) can be derived from the 55 55 3D terms. … … 64 64 % ================================================================ 65 65 \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 66 \label{ DYN_divcur_wzv}66 \label{sec:DYN_divcur_wzv} 67 67 68 68 %-------------------------------------------------------------------------------------------------------------- … … 70 70 %-------------------------------------------------------------------------------------------------------------- 71 71 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 72 \label{ DYN_divcur}72 \label{subsec:DYN_divcur} 73 73 74 74 The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: 75 \begin{equation} \label{ Eq_divcur_cur}75 \begin{equation} \label{eq:divcur_cur} 76 76 \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right] 77 77 -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) … … 79 79 80 80 The horizontal divergence is defined at a $T$-point. It is given by: 81 \begin{equation} \label{ Eq_divcur_div}81 \begin{equation} \label{eq:divcur_div} 82 82 \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 83 83 \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right] … … 102 102 %-------------------------------------------------------------------------------------------------------------- 103 103 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 104 \label{ DYN_sshwzv}104 \label{subsec:DYN_sshwzv} 105 105 106 106 The sea surface height is given by : 107 \begin{equation} \label{ Eq_dynspg_ssh}107 \begin{equation} \label{eq:dynspg_ssh} 108 108 \begin{aligned} 109 109 \frac{\partial \eta }{\partial t} … … 117 117 expressed in Kg/m$^2$/s (which is equal to mm/s), and $\rho _w$=1,035~Kg/m$^3$ 118 118 is the reference density of sea water (Boussinesq approximation). If river runoff is 119 expressed as a surface freshwater flux (see \ S\ref{SBC}) then \textit{emp} can be119 expressed as a surface freshwater flux (see \autoref{chap:SBC}) then \textit{emp} can be 120 120 written as the evaporation minus precipitation, minus the river runoff. 121 121 The sea-surface height is evaluated using exactly the same time stepping scheme 122 as the tracer equation \ eqref{Eq_tra_nxt}:122 as the tracer equation \autoref{eq:tra_nxt}: 123 123 a leapfrog scheme in combination with an Asselin time filter, $i.e.$ the velocity appearing 124 in \ eqref{Eq_dynspg_ssh} is centred in time (\textit{now} velocity).124 in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity). 125 125 This is of paramount importance. Replacing $T$ by the number $1$ in the tracer equation and summing 126 126 over the water column must lead to the sea surface height equation otherwise tracer content … … 129 129 The vertical velocity is computed by an upward integration of the horizontal 130 130 divergence starting at the bottom, taking into account the change of the thickness of the levels : 131 \begin{equation} \label{ Eq_wzv}131 \begin{equation} \label{eq:wzv} 132 132 \left\{ \begin{aligned} 133 133 &\left. w \right|_{k_b-1/2} \quad= 0 \qquad \text{where } k_b \text{ is the level just above the sea floor } \\ … … 141 141 of the level thicknesses, re-orientated downward. 142 142 \gmcomment{not sure of this... to be modified with the change in emp setting} 143 In the case of a linear free surface, the time derivative in \ eqref{Eq_wzv} disappears.143 In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears. 144 144 The upper boundary condition applies at a fixed level $z=0$. The top vertical velocity 145 145 is thus equal to the divergence of the barotropic transport ($i.e.$ the first term in the 146 right-hand-side of \ eqref{Eq_dynspg_ssh}).146 right-hand-side of \autoref{eq:dynspg_ssh}). 147 147 148 148 Note also that whereas the vertical velocity has the same discrete … … 150 150 in the second case, $w$ is the velocity normal to the $s$-surfaces. 151 151 Note also that the $k$-axis is re-orientated downwards in the \textsc{fortran} code compared 152 to the indexing used in the semi-discrete equations such as \ eqref{Eq_wzv}153 (see \S\ref{DOM_Num_Index_vertical}).152 to the indexing used in the semi-discrete equations such as \autoref{eq:wzv} 153 (see \autoref{subsec:DOM_Num_Index_vertical}). 154 154 155 155 … … 158 158 % ================================================================ 159 159 \section{Coriolis and advection: vector invariant form} 160 \label{ DYN_adv_cor_vect}160 \label{sec:DYN_adv_cor_vect} 161 161 %-----------------------------------------nam_dynadv---------------------------------------------------- 162 162 \forfile{../namelists/namdyn_adv} … … 171 171 time (\textit{now} velocity). 172 172 At the lateral boundaries either free slip, no slip or partial slip boundary 173 conditions are applied following Chap.\ref{LBC}.173 conditions are applied following \autoref{chap:LBC}. 174 174 175 175 % ------------------------------------------------------------------------------------------------------------- … … 177 177 % ------------------------------------------------------------------------------------------------------------- 178 178 \subsection{Vorticity term (\protect\mdl{dynvor})} 179 \label{ DYN_vor}179 \label{subsec:DYN_vor} 180 180 %------------------------------------------nam_dynvor---------------------------------------------------- 181 181 \forfile{../namelists/namdyn_vor} … … 188 188 the relative vorticity term and horizontal kinetic energy for the planetary vorticity 189 189 term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent 190 flow and horizontal kinetic energy (EEN scheme) (see Appendix~\ref{Apdx_C_vorEEN}). In the190 flow and horizontal kinetic energy (EEN scheme) (see \autoref{subsec:C_vorEEN}). In the 191 191 case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the 192 192 consistency of vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ = .true.}). … … 198 198 %------------------------------------------------------------- 199 199 \subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 200 \label{ DYN_vor_ens}200 \label{subsec:DYN_vor_ens} 201 201 202 202 In the enstrophy conserving case (ENS scheme), the discrete formulation of the … … 204 204 ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent 205 205 flow ($i.e.$ $\chi$=$0$), but does not conserve the total kinetic energy. It is given by: 206 \begin{equation} \label{ Eq_dynvor_ens}206 \begin{equation} \label{eq:dynvor_ens} 207 207 \left\{ 208 208 \begin{aligned} … … 219 219 %------------------------------------------------------------- 220 220 \subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 221 \label{ DYN_vor_ene}221 \label{subsec:DYN_vor_ene} 222 222 223 223 The kinetic energy conserving scheme (ENE scheme) conserves the global 224 224 kinetic energy but not the global enstrophy. It is given by: 225 \begin{equation} \label{ Eq_dynvor_ene}225 \begin{equation} \label{eq:dynvor_ene} 226 226 \left\{ \begin{aligned} 227 227 {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) … … 236 236 %------------------------------------------------------------- 237 237 \subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } 238 \label{ DYN_vor_mix}238 \label{subsec:DYN_vor_mix} 239 239 240 240 For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the 241 two previous schemes is used. It consists of the ENS scheme (\ ref{Eq_dynvor_ens})242 for the relative vorticity term, and of the ENE scheme (\ ref{Eq_dynvor_ene}) applied241 two previous schemes is used. It consists of the ENS scheme (\autoref{eq:dynvor_ens}) 242 for the relative vorticity term, and of the ENE scheme (\autoref{eq:dynvor_ene}) applied 243 243 to the planetary vorticity term. 244 \begin{equation} \label{ Eq_dynvor_mix}244 \begin{equation} \label{eq:dynvor_mix} 245 245 \left\{ { \begin{aligned} 246 246 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} … … 259 259 %------------------------------------------------------------- 260 260 \subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } 261 \label{ DYN_vor_een}261 \label{subsec:DYN_vor_een} 262 262 263 263 In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$ … … 277 277 The idea is to get rid of the double averaging by considering triad combinations of vorticity. 278 278 It is noteworthy that this solution is conceptually quite similar to the one proposed by 279 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see App.\ref{Apdx_C}).279 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 280 280 281 281 The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified 282 282 for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme. 283 283 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 284 \begin{equation} \label{ Eq_pot_vor}284 \begin{equation} \label{eq:pot_vor} 285 285 q = \frac{\zeta +f} {e_{3f} } 286 286 \end{equation} 287 where the relative vorticity is defined by (\ ref{Eq_divcur_cur}), the Coriolis parameter287 where the relative vorticity is defined by (\autoref{eq:divcur_cur}), the Coriolis parameter 288 288 is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: 289 \begin{equation} \label{ Eq_een_e3f}289 \begin{equation} \label{eq:een_e3f} 290 290 e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} 291 291 \end{equation} … … 294 294 \begin{figure}[!ht] \begin{center} 295 295 \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} 296 \caption{ \protect\label{ Fig_DYN_een_triad}296 \caption{ \protect\label{fig:DYN_een_triad} 297 297 Triads used in the energy and enstrophy conserving scheme (een) for 298 298 $u$-component (upper panel) and $v$-component (lower panel).} … … 300 300 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 301 301 302 A key point in \ eqref{Eq_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.302 A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 303 303 It uses the sum of masked t-point vertical scale factor divided either 304 304 by the sum of the four t-point masks (\np{nn\_een\_e3f}\forcode{ = 1}), … … 312 312 Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as 313 313 the following triad combinations of the neighbouring potential vorticities defined at f-points 314 ( Fig.~\ref{Fig_DYN_een_triad}):315 \begin{equation} \label{ Q_triads}314 (\autoref{fig:DYN_een_triad}): 315 \begin{equation} \label{eq:Q_triads} 316 316 _i^j \mathbb{Q}^{i_p}_{j_p} 317 317 = \frac{1}{12} \ \left( q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p} \right) … … 320 320 321 321 Finally, the vorticity terms are represented as: 322 \begin{equation} \label{ Eq_dynvor_een}322 \begin{equation} \label{eq:dynvor_een} 323 323 \left\{ { 324 324 \begin{aligned} … … 333 333 This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. 334 334 It conserves both total energy and potential enstrophy in the limit of horizontally 335 nondivergent flow ($i.e.$ $\chi$=$0$) (see Appendix~\ref{Apdx_C_vorEEN}).335 nondivergent flow ($i.e.$ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}). 336 336 Applied to a realistic ocean configuration, it has been shown that it leads to a significant 337 337 reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}. … … 344 344 %-------------------------------------------------------------------------------------------------------------- 345 345 \subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} 346 \label{ DYN_keg}347 348 As demonstrated in Appendix~\ref{Apdx_C}, there is a single discrete formulation346 \label{subsec:DYN_keg} 347 348 As demonstrated in \autoref{apdx:C}, there is a single discrete formulation 349 349 of the kinetic energy gradient term that, together with the formulation chosen for 350 350 the vertical advection (see below), conserves the total kinetic energy: 351 \begin{equation} \label{ Eq_dynkeg}351 \begin{equation} \label{eq:dynkeg} 352 352 \left\{ \begin{aligned} 353 353 -\frac{1}{2 \; e_{1u} } & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] \\ … … 360 360 %-------------------------------------------------------------------------------------------------------------- 361 361 \subsection{Vertical advection term (\protect\mdl{dynzad}) } 362 \label{ DYN_zad}362 \label{subsec:DYN_zad} 363 363 364 364 The discrete formulation of the vertical advection, together with the formulation 365 365 chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic 366 366 energy. Indeed, the change of KE due to the vertical advection is exactly 367 balanced by the change of KE due to the gradient of KE (see Appendix~\ref{Apdx_C}).368 \begin{equation} \label{ Eq_dynzad}367 balanced by the change of KE due to the gradient of KE (see \autoref{apdx:C}). 368 \begin{equation} \label{eq:dynzad} 369 369 \left\{ \begin{aligned} 370 370 -\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2} \;\delta _{k+1/2} \left[ u \right]\ }^{\,k} \\ … … 377 377 Note that in this case, a similar split-explicit time stepping should be used on 378 378 vertical advection of tracer to ensure a better stability, 379 an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \ S\ref{TRA_adv_tvd}).379 an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). 380 380 381 381 … … 384 384 % ================================================================ 385 385 \section{Coriolis and advection: flux form} 386 \label{ DYN_adv_cor_flux}386 \label{sec:DYN_adv_cor_flux} 387 387 %------------------------------------------nam_dynadv---------------------------------------------------- 388 388 \forfile{../namelists/namdyn_adv} … … 394 394 appearing in their expressions is centred in time (\textit{now} velocity). At the 395 395 lateral boundaries either free slip, no slip or partial slip boundary conditions 396 are applied following Chap.\ref{LBC}.396 are applied following \autoref{chap:LBC}. 397 397 398 398 … … 401 401 %-------------------------------------------------------------------------------------------------------------- 402 402 \subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } 403 \label{ DYN_cor_flux}403 \label{subsec:DYN_cor_flux} 404 404 405 405 In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis 406 406 parameter has been modified to account for the "metric" term. This altered 407 407 Coriolis parameter is thus discretised at $f$-points. It is given by: 408 \begin{multline} \label{ Eq_dyncor_metric}408 \begin{multline} \label{eq:dyncor_metric} 409 409 f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right) \\ 410 410 \equiv f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right] … … 412 412 \end{multline} 413 413 414 Any of the (\ ref{Eq_dynvor_ens}), (\ref{Eq_dynvor_ene}) and (\ref{Eq_dynvor_een})414 Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) 415 415 schemes can be used to compute the product of the Coriolis parameter and the 416 vorticity. However, the energy-conserving scheme (\ ref{Eq_dynvor_een}) has416 vorticity. However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has 417 417 exclusively been used to date. This term is evaluated using a leapfrog scheme, 418 418 $i.e.$ the velocity is centred in time (\textit{now} velocity). … … 422 422 %-------------------------------------------------------------------------------------------------------------- 423 423 \subsection{Flux form advection term (\protect\mdl{dynadv}) } 424 \label{ DYN_adv_flux}424 \label{subsec:DYN_adv_flux} 425 425 426 426 The discrete expression of the advection term is given by : 427 \begin{equation} \label{ Eq_dynadv}427 \begin{equation} \label{eq:dynadv} 428 428 \left\{ 429 429 \begin{aligned} … … 454 454 %------------------------------------------------------------- 455 455 \subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 456 \label{ DYN_adv_cen2}456 \label{subsec:DYN_adv_cen2} 457 457 458 458 In the centered $2^{nd}$ order formulation, the velocity is evaluated as the 459 459 mean of the two neighbouring points : 460 \begin{equation} \label{ Eq_dynadv_cen2}460 \begin{equation} \label{eq:dynadv_cen2} 461 461 \left\{ \begin{aligned} 462 462 u_T^{cen2} &=\overline u^{i } \quad & u_F^{cen2} &=\overline u^{j+1/2} \quad & u_{uw}^{cen2} &=\overline u^{k+1/2} \\ … … 475 475 %------------------------------------------------------------- 476 476 \subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 477 \label{ DYN_adv_ubs}477 \label{subsec:DYN_adv_ubs} 478 478 479 479 The UBS advection scheme is an upstream biased third order scheme based on 480 480 an upstream-biased parabolic interpolation. For example, the evaluation of 481 481 $u_T^{ubs} $ is done as follows: 482 \begin{equation} \label{ Eq_dynadv_ubs}482 \begin{equation} \label{eq:dynadv_ubs} 483 483 u_T^{ubs} =\overline u ^i-\;\frac{1}{6} \begin{cases} 484 484 u"_{i-1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i \geqslant 0$ } \\ … … 498 498 The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ 499 499 order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and 500 $u_{vw}^{ubs}$ in \ eqref{Eq_dynadv_cen2} are used. UBS is diffusive and is500 $u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used. UBS is diffusive and is 501 501 associated with vertical mixing of momentum. \gmcomment{ gm pursue the 502 502 sentence:Since vertical mixing of momentum is a source term of the TKE equation... } 503 503 504 For stability reasons, the first term in (\ref{Eq_dynadv_ubs}), which corresponds504 For stability reasons, the first term in (\autoref{eq:dynadv_ubs}), which corresponds 505 505 to a second order centred scheme, is evaluated using the \textit{now} velocity 506 506 (centred in time), while the second term, which is the diffusion part of the scheme, … … 510 510 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) 511 511 schemes only differ by one coefficient. Replacing $1/6$ by $1/8$ in 512 (\ ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.512 (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 513 513 This option is not available through a namelist parameter, since the $1/6$ coefficient 514 514 is hard coded. Nevertheless it is quite easy to make the substitution in the … … 526 526 % ================================================================ 527 527 \section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 528 \label{ DYN_hpg}528 \label{sec:DYN_hpg} 529 529 %------------------------------------------nam_dynhpg--------------------------------------------------- 530 530 \forfile{../namelists/namdyn_hpg} … … 547 547 %-------------------------------------------------------------------------------------------------------------- 548 548 \subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 549 \label{ DYN_hpg_zco}549 \label{subsec:DYN_hpg_zco} 550 550 551 551 The hydrostatic pressure can be obtained by integrating the hydrostatic equation … … 556 556 557 557 for $k=km$ (surface layer, $jk=1$ in the code) 558 \begin{equation} \label{ Eq_dynhpg_zco_surf}558 \begin{equation} \label{eq:dynhpg_zco_surf} 559 559 \left\{ \begin{aligned} 560 560 \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k=km} … … 566 566 567 567 for $1<k<km$ (interior layer) 568 \begin{equation} \label{ Eq_dynhpg_zco}568 \begin{equation} \label{eq:dynhpg_zco} 569 569 \left\{ \begin{aligned} 570 570 \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k} … … 577 577 \end{equation} 578 578 579 Note that the $1/2$ factor in (\ ref{Eq_dynhpg_zco_surf}) is adequate because of579 Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of 580 580 the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface 581 581 level ($z=0$). Note also that in case of variable volume level (\key{vvl} defined), the 582 surface pressure gradient is included in \ eqref{Eq_dynhpg_zco_surf} and \eqref{Eq_dynhpg_zco}582 surface pressure gradient is included in \autoref{eq:dynhpg_zco_surf} and \autoref{eq:dynhpg_zco} 583 583 through the space and time variations of the vertical scale factor $e_{3w}$. 584 584 … … 587 587 %-------------------------------------------------------------------------------------------------------------- 588 588 \subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 589 \label{ DYN_hpg_zps}589 \label{subsec:DYN_hpg_zps} 590 590 591 591 With partial bottom cells, tracers in horizontally adjacent cells generally live at … … 596 596 Apart from this modification, the horizontal hydrostatic pressure gradient evaluated 597 597 in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case. 598 As explained in detail in section \ S\ref{TRA_zpshde}, the nonlinearity of pressure598 As explained in detail in section \autoref{sec:TRA_zpshde}, the nonlinearity of pressure 599 599 effects in the equation of state is such that it is better to interpolate temperature and 600 600 salinity vertically before computing the density. Horizontal gradients of temperature 601 601 and salinity are needed for the TRA modules, which is the reason why the horizontal 602 602 gradients of density at the deepest model level are computed in module \mdl{zpsdhe} 603 located in the TRA directory and described in \ S\ref{TRA_zpshde}.603 located in the TRA directory and described in \autoref{sec:TRA_zpshde}. 604 604 605 605 %-------------------------------------------------------------------------------------------------------------- … … 607 607 %-------------------------------------------------------------------------------------------------------------- 608 608 \subsection{$S$- and $Z$-$S$-coordinates} 609 \label{ DYN_hpg_sco}609 \label{subsec:DYN_hpg_sco} 610 610 611 611 Pressure gradient formulations in an $s$-coordinate have been the subject of a vast … … 615 615 616 616 $\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.}) 617 \begin{equation} \label{ Eq_dynhpg_sco}617 \begin{equation} \label{eq:dynhpg_sco} 618 618 \left\{ \begin{aligned} 619 619 - \frac{1} {\rho_o \, e_{1u}} \; \delta _{i+1/2} \left[ p^h \right] … … 625 625 626 626 Where the first term is the pressure gradient along coordinates, computed as in 627 \ eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of627 \autoref{eq:dynhpg_zco_surf} - \autoref{eq:dynhpg_zco}, and $z_T$ is the depth of 628 628 the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 629 629 ($e_{3w}$). … … 637 637 (\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development) 638 638 639 Note that expression \ eqref{Eq_dynhpg_sco} is commonly used when the variable volume formulation is639 Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is 640 640 activated (\key{vvl}) because in that case, even with a flat bottom, the coordinate surfaces are not 641 641 horizontal but follow the free surface \citep{Levier2007}. The pressure jacobian scheme … … 648 648 649 649 \subsection{Ice shelf cavity} 650 \label{ DYN_hpg_isf}650 \label{subsec:DYN_hpg_isf} 651 651 Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and 652 652 the pressure gradient due to the ocean load. If cavity opened (\np{ln\_isfcav}\forcode{ = .true.}) these 2 terms can be … … 658 658 This top pressure is constant over time. A detailed description of this method is described in \citet{Losch2008}.\\ 659 659 660 $\bullet$ The ocean load is computed using the expression \ eqref{Eq_dynhpg_sco} described in \ref{DYN_hpg_sco}.660 $\bullet$ The ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in \autoref{subsec:DYN_hpg_sco}. 661 661 662 662 %-------------------------------------------------------------------------------------------------------------- … … 664 664 %-------------------------------------------------------------------------------------------------------------- 665 665 \subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})} 666 \label{ DYN_hpg_imp}666 \label{subsec:DYN_hpg_imp} 667 667 668 668 The default time differencing scheme used for the horizontal pressure gradient is … … 680 680 $\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}): 681 681 682 \begin{equation} \label{ Eq_dynhpg_lf}682 \begin{equation} \label{eq:dynhpg_lf} 683 683 \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 684 684 -\frac{1}{\rho _o \,e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right] … … 686 686 687 687 $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}): 688 \begin{equation} \label{ Eq_dynhpg_imp}688 \begin{equation} \label{eq:dynhpg_imp} 689 689 \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 690 690 -\frac{1}{4\,\rho _o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt} \right] 691 691 \end{equation} 692 692 693 The semi-implicit time scheme \ eqref{Eq_dynhpg_imp} is made possible without693 The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without 694 694 significant additional computation since the density can be updated to time level 695 695 $t+\rdt$ before computing the horizontal hydrostatic pressure gradient. It can 696 696 be easily shown that the stability limit associated with the hydrostatic pressure 697 gradient doubles using \ eqref{Eq_dynhpg_imp} compared to that using the698 standard leapfrog scheme \ eqref{Eq_dynhpg_lf}. Note that \eqref{Eq_dynhpg_imp}697 gradient doubles using \autoref{eq:dynhpg_imp} compared to that using the 698 standard leapfrog scheme \autoref{eq:dynhpg_lf}. Note that \autoref{eq:dynhpg_imp} 699 699 is equivalent to applying a time filter to the pressure gradient to eliminate high 700 frequency IGWs. Obviously, when using \ eqref{Eq_dynhpg_imp}, the doubling of700 frequency IGWs. Obviously, when using \autoref{eq:dynhpg_imp}, the doubling of 701 701 the time-step is achievable only if no other factors control the time-step, such as 702 702 the stability limits associated with advection or diffusion. … … 708 708 compute the hydrostatic pressure gradient (whatever the formulation) is evaluated 709 709 as follows: 710 \begin{equation} \label{ Eq_rho_flt}710 \begin{equation} \label{eq:rho_flt} 711 711 \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) 712 712 \quad \text{with} \quad … … 722 722 % ================================================================ 723 723 \section{Surface pressure gradient (\protect\mdl{dynspg})} 724 \label{ DYN_spg}724 \label{sec:DYN_spg} 725 725 %-----------------------------------------nam_dynspg---------------------------------------------------- 726 726 \forfile{../namelists/namdyn_spg} … … 730 730 731 731 Options are defined through the \ngn{namdyn\_spg} namelist variables. 732 The surface pressure gradient term is related to the representation of the free surface (\ S\ref{PE_hor_pg}).732 The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}). 733 733 The main distinction is between the fixed volume case (linear free surface) and the variable volume case 734 (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\ S\ref{PE_free_surface})734 (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\autoref{subsec:PE_free_surface}) 735 735 the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case 736 (\ S\ref{PE_free_surface}).736 (\autoref{subsec:PE_free_surface}). 737 737 With both linear and nonlinear free surface, external gravity waves are allowed in the equations, 738 738 which imposes a very small time step when an explicit time stepping is used. 739 739 Two methods are proposed to allow a longer time step for the three-dimensional equations: 740 the filtered free surface, which is a modification of the continuous equations (see \ eqref{Eq_PE_flt}),740 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}), 741 741 and the split-explicit free surface described below. 742 742 The extra term introduced in the filtered method is calculated implicitly, … … 745 745 746 746 The form of the surface pressure gradient term depends on how the user wants to handle 747 the fast external gravity waves that are a solution of the analytical equation (\ S\ref{PE_hor_pg}).747 the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}). 748 748 Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): 749 749 an explicit formulation which requires a small time step ; … … 761 761 %-------------------------------------------------------------------------------------------------------------- 762 762 \subsection{Explicit free surface (\protect\key{dynspg\_exp})} 763 \label{ DYN_spg_exp}763 \label{subsec:DYN_spg_exp} 764 764 765 765 In the explicit free surface formulation (\key{dynspg\_exp} defined), the model time step … … 767 767 The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time), 768 768 is thus simply given by : 769 \begin{equation} \label{ Eq_dynspg_exp}769 \begin{equation} \label{eq:dynspg_exp} 770 770 \left\{ \begin{aligned} 771 771 - \frac{1}{e_{1u}\,\rho_o} \; \delta _{i+1/2} \left[ \,\rho \,\eta\, \right] \\ … … 782 782 %-------------------------------------------------------------------------------------------------------------- 783 783 \subsection{Split-explicit free surface (\protect\key{dynspg\_ts})} 784 \label{ DYN_spg_ts}784 \label{subsec:DYN_spg_ts} 785 785 %------------------------------------------namsplit----------------------------------------------------------- 786 786 %\forfile{../namelists/namsplit} … … 792 792 equation and the associated barotropic velocity equations with a smaller time 793 793 step than $\rdt$, the time step used for the three dimensional prognostic 794 variables ( Fig.~\ref{Fig_DYN_dynspg_ts}).794 variables (\autoref{fig:DYN_dynspg_ts}). 795 795 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) 796 796 is provided through the \np{nn\_baro} namelist parameter as: … … 802 802 %%% 803 803 The barotropic mode solves the following equations: 804 \begin{subequations} \label{ Eq_BT}805 \begin{equation} \label{ Eq_BT_dyn}804 \begin{subequations} \label{eq:BT} 805 \begin{equation} \label{eq:BT_dyn} 806 806 \frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}= 807 807 -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h} … … 809 809 \end{equation} 810 810 811 \begin{equation} \label{ Eq_BT_ssh}811 \begin{equation} \label{eq:BT_ssh} 812 812 \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E 813 813 \end{equation} 814 814 \end{subequations} 815 where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. The third term on the right hand side of \ eqref{Eq_BT_dyn} represents the bottom stress (see section \S\ref{ZDF_bfr}), explicitly accounted for at each barotropic iteration. Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm detailed in \citet{Shchepetkin_McWilliams_OM05}. AB3-AM4 coefficients used in \NEMO follow the second-order accurate, "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} (see their figure 12, lower left).815 where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration. Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm detailed in \citet{Shchepetkin_McWilliams_OM05}. AB3-AM4 coefficients used in \NEMO follow the second-order accurate, "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} (see their figure 12, lower left). 816 816 817 817 %> > > > > > > > > > > > > > > > > > > > > > > > > > > > 818 818 \begin{figure}[!t] \begin{center} 819 819 \includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts} 820 \caption{ \protect\label{ Fig_DYN_dynspg_ts}820 \caption{ \protect\label{fig:DYN_dynspg_ts} 821 821 Schematic of the split-explicit time stepping scheme for the external 822 822 and internal modes. Time increases to the right. In this particular exemple, … … 827 827 The former are used to obtain time filtered quantities at $t+\rdt$ while the latter are used to obtain time averaged 828 828 transports to advect tracers. 829 a) Forward time integration: \ np{ln\_bt\_fw}\forcode{ = .true.},\np{ln\_bt\_av}\forcode{ = .true.}.830 b) Centred time integration: \ np{ln\_bt\_fw}\forcode{ = .false.},\np{ln\_bt\_av}\forcode{ = .true.}.831 c) Forward time integration with no time filtering (POM-like scheme): \ np{ln\_bt\_fw}\forcode{ = .true.},\np{ln\_bt\_av}\forcode{ = .false.}. }829 a) Forward time integration: \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .true.}. 830 b) Centred time integration: \protect\np{ln\_bt\_fw}\forcode{ = .false.}, \protect\np{ln\_bt\_av}\forcode{ = .true.}. 831 c) Forward time integration with no time filtering (POM-like scheme): \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .false.}. } 832 832 \end{center} \end{figure} 833 833 %> > > > > > > > > > > > > > > > > > > > > > > > > > > > 834 834 835 835 In the default case (\np{ln\_bt\_fw}\forcode{ = .true.}), the external mode is integrated 836 between \textit{now} and \textit{after} baroclinic time-steps ( Fig.~\ref{Fig_DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic836 between \textit{now} and \textit{after} baroclinic time-steps (\autoref{fig:DYN_dynspg_ts}a). To avoid aliasing of fast barotropic motions into three dimensional equations, time filtering is eventually applied on barotropic 837 837 quantities (\np{ln\_bt\_av}\forcode{ = .true.}). In that case, the integration is extended slightly beyond \textit{after} time step to provide time filtered quantities. 838 838 These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. … … 850 850 at \textit{now} time step. This implies to change the traditional order of computations in \NEMO: most of momentum 851 851 trends (including the barotropic mode calculation) updated first, tracers' after. This \textit{de facto} makes semi-implicit hydrostatic 852 pressure gradient (see section \ S\ref{DYN_hpg_imp}) and time splitting not compatible.852 pressure gradient (see section \autoref{subsec:DYN_hpg_imp}) and time splitting not compatible. 853 853 Advective barotropic velocities are obtained by using a secondary set of filtering weights, uniquely defined from the filter 854 854 coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}). Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to obtain exact conservation. … … 872 872 scheme using the small barotropic time step $\rdt$. We have 873 873 874 \begin{equation} \label{ DYN_spg_ts_eta}874 \begin{equation} \label{eq:DYN_spg_ts_eta} 875 875 \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1}) 876 876 = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right] 877 877 \end{equation} 878 \begin{multline} \label{ DYN_spg_ts_u}878 \begin{multline} \label{eq:DYN_spg_ts_u} 879 879 \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}) \\ 880 880 = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n}) … … 886 886 and $U^{(b)}$ denotes the baroclinic time at which the vertically integrated forcing $\textbf{M}(\tau)$ (note that this forcing includes the surface freshwater forcing), the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$, and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over a single cycle. This is also the time 887 887 that sets the barotropic time steps via 888 \begin{equation} \label{ DYN_spg_ts_t}888 \begin{equation} \label{eq:DYN_spg_ts_t} 889 889 t_n=\tau+n\rdt 890 890 \end{equation} 891 891 with $n$ an integer. The density scaled surface pressure is evaluated via 892 \begin{equation} \label{ DYN_spg_ts_ps}892 \begin{equation} \label{eq:DYN_spg_ts_ps} 893 893 p_s^{(b)}(\tau,t_{n}) = \begin{cases} 894 894 g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o & \text{non-linear case} \\ … … 897 897 \end{equation} 898 898 To get started, we assume the following initial conditions 899 \begin{equation} \label{ DYN_spg_ts_eta}899 \begin{equation} \label{eq:DYN_spg_ts_eta} 900 900 \begin{split} 901 901 \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} … … 905 905 \end{equation} 906 906 with 907 \begin{equation} \label{ DYN_spg_ts_etaF}907 \begin{equation} \label{eq:DYN_spg_ts_etaF} 908 908 \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n}) 909 909 \end{equation} 910 910 the time averaged surface height taken from the previous barotropic cycle. Likewise, 911 \begin{equation} \label{ DYN_spg_ts_u}911 \begin{equation} \label{eq:DYN_spg_ts_u} 912 912 \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ 913 913 \\ … … 915 915 \end{equation} 916 916 with 917 \begin{equation} \label{ DYN_spg_ts_u}917 \begin{equation} \label{eq:DYN_spg_ts_u} 918 918 \overline{\textbf{U}^{(b)}(\tau)} 919 919 = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n}) … … 922 922 923 923 Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ , the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at baroclinic time $\tau + \rdt \tau$ 924 \begin{equation} \label{ DYN_spg_ts_u}924 \begin{equation} \label{eq:DYN_spg_ts_u} 925 925 \textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} 926 926 = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n}) … … 928 928 The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using the following form 929 929 930 \begin{equation} \label{ DYN_spg_ts_ssh}930 \begin{equation} \label{eq:DYN_spg_ts_ssh} 931 931 \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] 932 932 \end{equation} … … 935 935 936 936 In general, some form of time filter is needed to maintain integrity of the surface 937 height field due to the leap-frog splitting mode in equation \ ref{DYN_spg_ts_ssh}. We937 height field due to the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. We 938 938 have tried various forms of such filtering, with the following method discussed in 939 939 \cite{Griffies_al_MWR01} chosen due to its stability and reasonably good maintenance of 940 tracer conservation properties (see Section??)941 942 \begin{equation} \label{ DYN_spg_ts_sshf}940 tracer conservation properties (see ??) 941 942 \begin{equation} \label{eq:DYN_spg_ts_sshf} 943 943 \eta^{F}(\tau-\Delta) = \overline{\eta^{(b)}(\tau)} 944 944 \end{equation} 945 945 Another approach tried was 946 946 947 \begin{equation} \label{ DYN_spg_ts_sshf2}947 \begin{equation} \label{eq:DYN_spg_ts_sshf2} 948 948 \eta^{F}(\tau-\Delta) = \eta(\tau) 949 949 + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt) … … 953 953 which is useful since it isolates all the time filtering aspects into the term multiplied 954 954 by $\alpha$. This isolation allows for an easy check that tracer conservation is exact when 955 eliminating tracer and surface height time filtering (see Section ?? for more complete discussion). However, in the general case with a non-zero $\alpha$, the filter \ref{DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.955 eliminating tracer and surface height time filtering (see ?? for more complete discussion). However, in the general case with a non-zero $\alpha$, the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended. 956 956 957 957 } %%end gm comment (copy of griffies book) … … 964 964 %-------------------------------------------------------------------------------------------------------------- 965 965 \subsection{Filtered free surface (\protect\key{dynspg\_flt})} 966 \label{ DYN_spg_fltp}966 \label{subsec:DYN_spg_fltp} 967 967 968 968 The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation. 969 The extra term introduced in the equations (see \ S\ref{PE_free_surface}) is solved implicitly.970 The elliptic solvers available in the code are documented in \ S\ref{MISC}.969 The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly. 970 The elliptic solvers available in the code are documented in \autoref{chap:MISC}. 971 971 972 972 %% gm %%======>>>> given here the discrete eqs provided to the solver 973 973 \gmcomment{ %%% copy from chap-model basics 974 \begin{equation} \label{ Eq_spg_flt}974 \begin{equation} \label{eq:spg_flt} 975 975 \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} 976 976 - g \nabla \left( \tilde{\rho} \ \eta \right) … … 980 980 $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ 981 981 represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 982 non-linear and viscous terms in \ eqref{Eq_PE_dyn}.982 non-linear and viscous terms in \autoref{eq:PE_dyn}. 983 983 } %end gmcomment 984 984 … … 990 990 % ================================================================ 991 991 \section{Lateral diffusion term and operators (\protect\mdl{dynldf})} 992 \label{ DYN_ldf}992 \label{sec:DYN_ldf} 993 993 %------------------------------------------nam_dynldf---------------------------------------------------- 994 994 \forfile{../namelists/namdyn_ldf} … … 999 999 (rotated or not) or biharmonic operators. The coefficients may be constant 1000 1000 or spatially variable; the description of the coefficients is found in the chapter 1001 on lateral physics ( Chap.\ref{LDF}). The lateral diffusion of momentum is1001 on lateral physics (\autoref{chap:LDF}). The lateral diffusion of momentum is 1002 1002 evaluated using a forward scheme, $i.e.$ the velocity appearing in its expression 1003 1003 is the \textit{before} velocity in time, except for the pure vertical component 1004 1004 that appears when a tensor of rotation is used. This latter term is solved 1005 implicitly together with the vertical diffusion term (see \ S\ref{STP})1005 implicitly together with the vertical diffusion term (see \autoref{chap:STP}) 1006 1006 1007 1007 At the lateral boundaries either free slip, no slip or partial slip boundary 1008 conditions are applied according to the user's choice (see Chap.\ref{LBC}).1008 conditions are applied according to the user's choice (see \autoref{chap:LBC}). 1009 1009 1010 1010 \gmcomment{ … … 1025 1025 \subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})] 1026 1026 {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 1027 \label{ DYN_ldf_lap}1027 \label{subsec:DYN_ldf_lap} 1028 1028 1029 1029 For lateral iso-level diffusion, the discrete operator is: 1030 \begin{equation} \label{ Eq_dynldf_lap}1030 \begin{equation} \label{eq:dynldf_lap} 1031 1031 \left\{ \begin{aligned} 1032 1032 D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta _{i+1/2} \left[ {A_T^{lm} … … 1040 1040 \end{equation} 1041 1041 1042 As explained in \ S\ref{PE_ldf}, this formulation (as the gradient of a divergence1042 As explained in \autoref{subsec:PE_ldf}, this formulation (as the gradient of a divergence 1043 1043 and curl of the vorticity) preserves symmetry and ensures a complete 1044 1044 separation between the vorticity and divergence parts of the momentum diffusion. … … 1049 1049 \subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})] 1050 1050 {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 1051 \label{ DYN_ldf_iso}1051 \label{subsec:DYN_ldf_iso} 1052 1052 1053 1053 A rotation of the lateral momentum diffusion operator is needed in several cases: … … 1061 1061 constraints on the stress tensor such as symmetry. The resulting discrete 1062 1062 representation is: 1063 \begin{equation} \label{ Eq_dyn_ldf_iso}1063 \begin{equation} \label{eq:dyn_ldf_iso} 1064 1064 \begin{split} 1065 1065 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ … … 1111 1111 diffusion operator acts and the surface of computation ($z$- or $s$-surfaces). 1112 1112 The way these slopes are evaluated is given in the lateral physics chapter 1113 ( Chap.\ref{LDF}).1113 (\autoref{chap:LDF}). 1114 1114 1115 1115 %-------------------------------------------------------------------------------------------------------------- … … 1118 1118 \subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})] 1119 1119 {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 1120 \label{ DYN_ldf_bilap}1120 \label{subsec:DYN_ldf_bilap} 1121 1121 1122 1122 The lateral fourth order operator formulation on momentum is obtained by 1123 applying \ eqref{Eq_dynldf_lap} twice. It requires an additional assumption on1123 applying \autoref{eq:dynldf_lap} twice. It requires an additional assumption on 1124 1124 boundary conditions: the first derivative term normal to the coast depends on 1125 1125 the free or no-slip lateral boundary conditions chosen, while the third 1126 derivative terms normal to the coast are set to zero (see Chap.\ref{LBC}).1126 derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}). 1127 1127 %%% 1128 1128 \gmcomment{add a remark on the the change in the position of the coefficient} … … 1133 1133 % ================================================================ 1134 1134 \section{Vertical diffusion term (\protect\mdl{dynzdf})} 1135 \label{ DYN_zdf}1135 \label{sec:DYN_zdf} 1136 1136 %----------------------------------------------namzdf------------------------------------------------------ 1137 1137 \forfile{../namelists/namzdf} … … 1145 1145 scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using a time splitting technique 1146 1146 (\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme 1147 (\np{ln\_zdfexp}\forcode{ = .false.}) (see \ S\ref{STP}). Note that namelist variables1147 (\np{ln\_zdfexp}\forcode{ = .false.}) (see \autoref{chap:STP}). Note that namelist variables 1148 1148 \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 1149 1149 1150 1150 The formulation of the vertical subgrid scale physics is the same whatever 1151 1151 the vertical coordinate is. The vertical diffusion operators given by 1152 \ eqref{Eq_PE_zdf} take the following semi-discrete space form:1153 \begin{equation} \label{ Eq_dynzdf}1152 \autoref{eq:PE_zdf} take the following semi-discrete space form: 1153 \begin{equation} \label{eq:dynzdf} 1154 1154 \left\{ \begin{aligned} 1155 1155 D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta _k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } … … 1162 1162 where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and 1163 1163 diffusivity coefficients. The way these coefficients are evaluated 1164 depends on the vertical physics used (see \ S\ref{ZDF}).1164 depends on the vertical physics used (see \autoref{chap:ZDF}). 1165 1165 1166 1166 The surface boundary condition on momentum is the stress exerted by 1167 1167 the wind. At the surface, the momentum fluxes are prescribed as the boundary 1168 1168 condition on the vertical turbulent momentum fluxes, 1169 \begin{equation} \label{ Eq_dynzdf_sbc}1169 \begin{equation} \label{eq:dynzdf_sbc} 1170 1170 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 1171 1171 = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v } … … 1177 1177 is small (when no mixed layer scheme is used) the surface stress enters only 1178 1178 the top model level, as a body force. The surface wind stress is calculated 1179 in the surface module routines (SBC, see Chap.\ref{SBC})1179 in the surface module routines (SBC, see \autoref{chap:SBC}) 1180 1180 1181 1181 The turbulent flux of momentum at the bottom of the ocean is specified through 1182 a bottom friction parameterisation (see \ S\ref{ZDF_bfr})1182 a bottom friction parameterisation (see \autoref{sec:ZDF_bfr}) 1183 1183 1184 1184 % ================================================================ … … 1186 1186 % ================================================================ 1187 1187 \section{External forcings} 1188 \label{ DYN_forcing}1188 \label{sec:DYN_forcing} 1189 1189 1190 1190 Besides the surface and bottom stresses (see the above section) which are … … 1192 1192 may enter the dynamical equations by affecting the surface pressure gradient. 1193 1193 1194 (1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \ S\ref{SBC_apr}), the atmospheric pressure is taken1194 (1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \autoref{sec:SBC_apr}), the atmospheric pressure is taken 1195 1195 into account when computing the surface pressure gradient. 1196 1196 1197 (2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \ S\ref{SBC_tide}),1197 (2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \autoref{sec:SBC_tide}), 1198 1198 the tidal potential is taken into account when computing the surface pressure gradient. 1199 1199 … … 1209 1209 % ================================================================ 1210 1210 \section{Time evolution term (\protect\mdl{dynnxt})} 1211 \label{ DYN_nxt}1211 \label{sec:DYN_nxt} 1212 1212 1213 1213 %----------------------------------------------namdom---------------------------------------------------- … … 1218 1218 The general framework for dynamics time stepping is a leap-frog scheme, 1219 1219 $i.e.$ a three level centred time scheme associated with an Asselin time filter 1220 (cf. Chap.\ref{STP}). The scheme is applied to the velocity, except when using1221 the flux form of momentum advection (cf. \ S\ref{DYN_adv_cor_flux}) in the variable1220 (cf. \autoref{chap:STP}). The scheme is applied to the velocity, except when using 1221 the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) in the variable 1222 1222 volume case (\key{vvl} defined), where it has to be applied to the thickness 1223 weighted velocity (see \ S\ref{Apdx_A_momentum})1223 weighted velocity (see \autoref{sec:A_momentum}) 1224 1224 1225 1225 $\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}\forcode{ = .true.} ; \key{vvl} not defined): 1226 \begin{equation} \label{ Eq_dynnxt_vec}1226 \begin{equation} \label{eq:dynnxt_vec} 1227 1227 \left\{ \begin{aligned} 1228 1228 &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt \ \text{RHS}_u^t \\ … … 1232 1232 1233 1233 $\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}\forcode{ = .false.} ; \key{vvl} defined): 1234 \begin{equation} \label{ Eq_dynnxt_flux}1234 \begin{equation} \label{eq:dynnxt_flux} 1235 1235 \left\{ \begin{aligned} 1236 1236 &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t \\
Note: See TracChangeset
for help on using the changeset viewer.
