Changeset 11564 for NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex/NEMO/subfiles/chap_DYN.tex
- Timestamp:
- 2019-09-18T16:11:52+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex
- Property svn:ignore
-
old new 1 *.aux 2 *.bbl 3 *.blg 4 *.dvi 5 *.fdb* 6 *.fls 7 *.idx 8 *.ilg 9 *.ind 10 *.log 11 *.maf 12 *.mtc* 13 *.out 14 *.pdf 15 *.toc 16 _minted-* 1 figures
-
- Property svn:ignore
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex/NEMO
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/doc/latex/NEMO/subfiles/chap_DYN.tex
r10499 r11564 8 8 \label{chap:DYN} 9 9 10 \ minitoc10 \chaptertoc 11 11 12 12 Using the representation described in \autoref{chap:DOM}, … … 36 36 (surface wind stress calculation using bulk formulae, estimation of mixing coefficients) 37 37 that are carried out in modules SBC, LDF and ZDF and are described in 38 \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 38 \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 the horizontal divergence, 41 41 curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module). 42 42 43 The different options available to the user are managed by namelist variables. 44 For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, 43 The different options available to the user are managed by namelist variables. 44 For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, 45 45 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 46 If a CPP key is used for this term its name is \key{ttt}.46 %If a CPP key is used for this term its name is \key{ttt}. 47 47 The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory, 48 48 and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine. 49 49 50 50 The user has the option of extracting and outputting each tendency term from the 3D momentum equations 51 (\ key{trddyn} defined), as described in \autoref{chap:MISC}.52 Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \ key{trdvor} is defined)51 (\texttt{trddyn?} defined), as described in \autoref{chap:MISC}. 52 Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined) 53 53 can be derived from the 3D terms. 54 54 %%% 55 \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does 55 \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does 56 56 MISC correspond to "extracting tendency terms" or "vorticity balance"?} 57 57 … … 65 65 % Horizontal divergence and relative vorticity 66 66 %-------------------------------------------------------------------------------------------------------------- 67 \subsection {Horizontal divergence and relative vorticity (\protect\mdl{divcur})}67 \subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 68 68 \label{subsec:DYN_divcur} 69 69 70 The vorticity is defined at an $f$-point (\ie corner point) as follows:71 \begin{equation} 72 \label{eq: divcur_cur}70 The vorticity is defined at an $f$-point (\ie\ corner point) as follows: 71 \begin{equation} 72 \label{eq:DYN_divcur_cur} 73 73 \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right] 74 74 -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) 75 \end{equation} 75 \end{equation} 76 76 77 77 The horizontal divergence is defined at a $T$-point. 78 78 It is given by: 79 79 \[ 80 % \label{eq: divcur_div}80 % \label{eq:DYN_divcur_div} 81 81 \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 82 82 \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right] … … 96 96 ensure perfect restartability. 97 97 The vorticity and divergence at the \textit{now} time step are used for the computation of 98 the nonlinear advection and of the vertical velocity respectively. 98 the nonlinear advection and of the vertical velocity respectively. 99 99 100 100 %-------------------------------------------------------------------------------------------------------------- 101 101 % Sea Surface Height evolution 102 102 %-------------------------------------------------------------------------------------------------------------- 103 \subsection {Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}103 \subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 104 104 \label{subsec:DYN_sshwzv} 105 105 106 106 The sea surface height is given by: 107 107 \begin{equation} 108 \label{eq: dynspg_ssh}108 \label{eq:DYN_spg_ssh} 109 109 \begin{aligned} 110 110 \frac{\partial \eta }{\partial t} … … 115 115 \end{aligned} 116 116 \end{equation} 117 where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), 117 where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), 118 118 expressed in Kg/m$^2$/s (which is equal to mm/s), 119 119 and $\rho_w$=1,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation). 120 120 If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then 121 \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. 121 \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. 122 122 The sea-surface height is evaluated using exactly the same time stepping scheme as 123 the tracer equation \autoref{eq: tra_nxt}:123 the tracer equation \autoref{eq:TRA_nxt}: 124 124 a leapfrog scheme in combination with an Asselin time filter, 125 \ie the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity).125 \ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity). 126 126 This is of paramount importance. 127 127 Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to 128 128 the sea surface height equation otherwise tracer content will not be conserved 129 \citep{ Griffies_al_MWR01, Leclair_Madec_OM09}.129 \citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}. 130 130 131 131 The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, 132 132 taking into account the change of the thickness of the levels: 133 133 \begin{equation} 134 \label{eq: wzv}134 \label{eq:DYN_wzv} 135 135 \left\{ 136 136 \begin{aligned} … … 142 142 \end{equation} 143 143 144 In the case of a non-linear free surface (\ key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$,144 In the case of a non-linear free surface (\texttt{vvl?}), the top vertical velocity is $-\textit{emp}/\rho_w$, 145 145 as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses, 146 146 re-orientated downward. 147 147 \gmcomment{not sure of this... to be modified with the change in emp setting} 148 In the case of a linear free surface, the time derivative in \autoref{eq: wzv} disappears.148 In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears. 149 149 The upper boundary condition applies at a fixed level $z=0$. 150 150 The top vertical velocity is thus equal to the divergence of the barotropic transport 151 (\ie the first term in the right-hand-side of \autoref{eq:dynspg_ssh}).151 (\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}). 152 152 153 153 Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates, 154 154 its physical meaning is not the same: 155 155 in the second case, $w$ is the velocity normal to the $s$-surfaces. 156 Note also that the $k$-axis is re-orientated downwards in the \fortran code compared to157 the indexing used in the semi-discrete equations such as \autoref{eq: wzv}158 (see \autoref{subsec:DOM_Num_Index_vertical}). 156 Note also that the $k$-axis is re-orientated downwards in the \fortran\ code compared to 157 the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv} 158 (see \autoref{subsec:DOM_Num_Index_vertical}). 159 159 160 160 … … 166 166 %-----------------------------------------nam_dynadv---------------------------------------------------- 167 167 168 \nlst{namdyn_adv} 168 \begin{listing} 169 \nlst{namdyn_adv} 170 \caption{\texttt{namdyn\_adv}} 171 \label{lst:namdyn_adv} 172 \end{listing} 169 173 %------------------------------------------------------------------------------------------------------------- 170 174 171 175 The vector invariant form of the momentum equations is the one most often used in 172 applications of the \NEMO ocean model.176 applications of the \NEMO\ ocean model. 173 177 The flux form option (see next section) has been present since version $2$. 174 Options are defined through the \n gn{namdyn\_adv} namelist variables Coriolis and178 Options are defined through the \nam{dyn\_adv} namelist variables Coriolis and 175 179 momentum advection terms are evaluated using a leapfrog scheme, 176 \ie the velocity appearing in these expressions is centred in time (\textit{now} velocity).180 \ie\ the velocity appearing in these expressions is centred in time (\textit{now} velocity). 177 181 At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following 178 182 \autoref{chap:LBC}. 179 183 180 184 % ------------------------------------------------------------------------------------------------------------- 181 % Vorticity term 185 % Vorticity term 182 186 % ------------------------------------------------------------------------------------------------------------- 183 \subsection {Vorticity term (\protect\mdl{dynvor})}187 \subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} 184 188 \label{subsec:DYN_vor} 185 189 %------------------------------------------nam_dynvor---------------------------------------------------- 186 190 187 \nlst{namdyn_vor} 191 \begin{listing} 192 \nlst{namdyn_vor} 193 \caption{\texttt{namdyn\_vor}} 194 \label{lst:namdyn_vor} 195 \end{listing} 188 196 %------------------------------------------------------------------------------------------------------------- 189 197 190 Options are defined through the \n gn{namdyn\_vor} namelist variables.191 Four discretisations of the vorticity term (\ np{ln\_dynvor\_xxx}\forcode{ =.true.}) are available:198 Options are defined through the \nam{dyn\_vor} namelist variables. 199 Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available: 192 200 conserving potential enstrophy of horizontally non-divergent flow (ENS scheme); 193 201 conserving horizontal kinetic energy (ENE scheme); … … 195 203 horizontal kinetic energy for the planetary vorticity term (MIX scheme); 196 204 or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy 197 (EEN scheme) (see \autoref{subsec: C_vorEEN}).205 (EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}). 198 206 In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of 199 vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ =.true.}).207 vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{=.true.}). 200 208 The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. 201 209 … … 203 211 % enstrophy conserving scheme 204 212 %------------------------------------------------------------- 205 \subsubsection {Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}213 \subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})]{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 206 214 \label{subsec:DYN_vor_ens} 207 215 208 216 In the enstrophy conserving case (ENS scheme), 209 217 the discrete formulation of the vorticity term provides a global conservation of the enstrophy 210 ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie $\chi$=$0$),218 ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie\ $\chi$=$0$), 211 219 but does not conserve the total kinetic energy. 212 220 It is given by: 213 221 \begin{equation} 214 \label{eq: dynvor_ens}222 \label{eq:DYN_vor_ens} 215 223 \left\{ 216 224 \begin{aligned} … … 221 229 \end{aligned} 222 230 \right. 223 \end{equation} 231 \end{equation} 224 232 225 233 %------------------------------------------------------------- 226 234 % energy conserving scheme 227 235 %------------------------------------------------------------- 228 \subsubsection {Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}236 \subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})]{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 229 237 \label{subsec:DYN_vor_ene} 230 238 … … 232 240 It is given by: 233 241 \begin{equation} 234 \label{eq: dynvor_ene}242 \label{eq:DYN_vor_ene} 235 243 \left\{ 236 244 \begin{aligned} … … 241 249 \end{aligned} 242 250 \right. 243 \end{equation} 251 \end{equation} 244 252 245 253 %------------------------------------------------------------- 246 254 % mix energy/enstrophy conserving scheme 247 255 %------------------------------------------------------------- 248 \subsubsection {Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})}256 \subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})} 249 257 \label{subsec:DYN_vor_mix} 250 258 251 259 For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used. 252 It consists of the ENS scheme (\autoref{eq: dynvor_ens}) for the relative vorticity term,253 and of the ENE scheme (\autoref{eq: dynvor_ene}) applied to the planetary vorticity term.254 \[ 255 % \label{eq: dynvor_mix}260 It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term, 261 and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term. 262 \[ 263 % \label{eq:DYN_vor_mix} 256 264 \left\{ { 257 265 \begin{aligned} … … 271 279 % energy and enstrophy conserving scheme 272 280 %------------------------------------------------------------- 273 \subsubsection {Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}281 \subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})]{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})} 274 282 \label{subsec:DYN_vor_een} 275 283 … … 278 286 the presence of grid point oscillation structures that will be invisible to the operator. 279 287 These structures are \textit{computational modes} that will be at least partly damped by 280 the momentum diffusion operator (\ie the subgrid-scale advection), but not by the resolved advection term.288 the momentum diffusion operator (\ie\ the subgrid-scale advection), but not by the resolved advection term. 281 289 The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field. 282 290 Such noise would result in more noise in the vertical velocity field, an undesirable feature. … … 284 292 $u$ and $v$ are located at different grid points, 285 293 a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid. 286 \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) 294 \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) 287 295 Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} 288 296 289 A very nice solution to the problem of double averaging was proposed by \citet{ Arakawa_Hsu_MWR90}.297 A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}. 290 298 The idea is to get rid of the double averaging by considering triad combinations of vorticity. 291 299 It is noteworthy that this solution is conceptually quite similar to the one proposed by 292 \citep{ Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).293 294 The \citet{ Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified295 for spherical coordinates as described by \citet{ Arakawa_Lamb_MWR81} to obtain the EEN scheme.296 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 297 \[ 298 % \label{eq: pot_vor}300 \citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}). 301 302 The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified 303 for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme. 304 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 305 \[ 306 % \label{eq:DYN_pot_vor} 299 307 q = \frac{\zeta +f} {e_{3f} } 300 308 \] 301 where the relative vorticity is defined by (\autoref{eq: divcur_cur}),302 the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: 303 \begin{equation} 304 \label{eq: een_e3f}309 where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}), 310 the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: 311 \begin{equation} 312 \label{eq:DYN_een_e3f} 305 313 e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} 306 314 \end{equation} … … 308 316 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 309 317 \begin{figure}[!ht] 310 \begin{center} 311 \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} 312 \caption{ 313 \protect\label{fig:DYN_een_triad} 314 Triads used in the energy and enstrophy conserving scheme (een) for 315 $u$-component (upper panel) and $v$-component (lower panel). 316 } 317 \end{center} 318 \centering 319 \includegraphics[width=\textwidth]{Fig_DYN_een_triad} 320 \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{ 321 Triads used in the energy and enstrophy conserving scheme (EEN) for 322 $u$-component (upper panel) and $v$-component (lower panel).} 323 \label{fig:DYN_een_triad} 318 324 \end{figure} 319 325 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 320 326 321 A key point in \autoref{eq: een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.327 A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 322 328 It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks 323 (\np{nn\_een\_e3f}\forcode{ = 1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{ =.true.}).329 (\np{nn\_een\_e3f}\forcode{=1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{=.true.}). 324 330 The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and 325 331 extends by continuity the value of $e_{3f}$ into the land areas. … … 327 333 (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) 328 334 that tends to reinforce the topostrophy of the flow 329 (\ie the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.335 (\ie\ the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}. 330 336 331 337 Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as 332 338 the following triad combinations of the neighbouring potential vorticities defined at f-points 333 (\autoref{fig:DYN_een_triad}): 334 \begin{equation} 335 \label{eq: Q_triads}339 (\autoref{fig:DYN_een_triad}): 340 \begin{equation} 341 \label{eq:DYN_Q_triads} 336 342 _i^j \mathbb{Q}^{i_p}_{j_p} 337 343 = \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) 338 344 \end{equation} 339 where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. 340 341 Finally, the vorticity terms are represented as: 342 \begin{equation} 343 \label{eq: dynvor_een}345 where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. 346 347 Finally, the vorticity terms are represented as: 348 \begin{equation} 349 \label{eq:DYN_vor_een} 344 350 \left\{ { 345 351 \begin{aligned} … … 350 356 \end{aligned} 351 357 } \right. 352 \end{equation} 358 \end{equation} 353 359 354 360 This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. 355 361 It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow 356 (\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).362 (\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}). 357 363 Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of 358 the noise in the vertical velocity field \citep{ Le_Sommer_al_OM09}.364 the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. 359 365 Furthermore, used in combination with a partial steps representation of bottom topography, 360 366 it improves the interaction between current and topography, 361 leading to a larger topostrophy of the flow \citep{ Barnier_al_OD06, Penduff_al_OS07}.367 leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. 362 368 363 369 %-------------------------------------------------------------------------------------------------------------- 364 370 % Kinetic Energy Gradient term 365 371 %-------------------------------------------------------------------------------------------------------------- 366 \subsection {Kinetic energy gradient term (\protect\mdl{dynkeg})}372 \subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} 367 373 \label{subsec:DYN_keg} 368 374 369 As demonstrated in \autoref{apdx: C},375 As demonstrated in \autoref{apdx:INVARIANTS}, 370 376 there is a single discrete formulation of the kinetic energy gradient term that, 371 377 together with the formulation chosen for the vertical advection (see below), 372 378 conserves the total kinetic energy: 373 379 \[ 374 % \label{eq: dynkeg}380 % \label{eq:DYN_keg} 375 381 \left\{ 376 382 \begin{aligned} … … 384 390 % Vertical advection term 385 391 %-------------------------------------------------------------------------------------------------------------- 386 \subsection {Vertical advection term (\protect\mdl{dynzad})}392 \subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} 387 393 \label{subsec:DYN_zad} 388 394 … … 391 397 conserves the total kinetic energy. 392 398 Indeed, the change of KE due to the vertical advection is exactly balanced by 393 the change of KE due to the gradient of KE (see \autoref{apdx: C}).394 \[ 395 % \label{eq: dynzad}399 the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}). 400 \[ 401 % \label{eq:DYN_zad} 396 402 \left\{ 397 403 \begin{aligned} … … 401 407 \right. 402 408 \] 403 When \np{ln\_dynzad\_zts}\forcode{ =.true.},409 When \np{ln\_dynzad\_zts}\forcode{=.true.}, 404 410 a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. 405 This option can be useful when the value of the timestep is limited by vertical advection \citep{ Lemarie_OM2015}.411 This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}. 406 412 Note that in this case, 407 413 a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, … … 416 422 %------------------------------------------nam_dynadv---------------------------------------------------- 417 423 418 \nlst{namdyn_adv}419 424 %------------------------------------------------------------------------------------------------------------- 420 425 421 Options are defined through the \n gn{namdyn\_adv} namelist variables.426 Options are defined through the \nam{dyn\_adv} namelist variables. 422 427 In the flux form (as in the vector invariant form), 423 428 the Coriolis and momentum advection terms are evaluated using a leapfrog scheme, 424 \ie the velocity appearing in their expressions is centred in time (\textit{now} velocity).429 \ie\ the velocity appearing in their expressions is centred in time (\textit{now} velocity). 425 430 At the lateral boundaries either free slip, 426 431 no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. … … 430 435 % Coriolis plus curvature metric terms 431 436 %-------------------------------------------------------------------------------------------------------------- 432 \subsection {Coriolis plus curvature metric terms (\protect\mdl{dynvor})}437 \subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 433 438 \label{subsec:DYN_cor_flux} 434 439 435 440 In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term. 436 441 This altered Coriolis parameter is thus discretised at $f$-points. 437 It is given by: 442 It is given by: 438 443 \begin{multline*} 439 % \label{eq: dyncor_metric}444 % \label{eq:DYN_cor_metric} 440 445 f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right) \\ 441 446 \equiv f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right] 442 447 - \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right] } \ \right) 443 \end{multline*} 444 445 Any of the (\autoref{eq: dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to448 \end{multline*} 449 450 Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to 446 451 compute the product of the Coriolis parameter and the vorticity. 447 However, the energy-conserving scheme (\autoref{eq: dynvor_een}) has exclusively been used to date.448 This term is evaluated using a leapfrog scheme, \ie the velocity is centred in time (\textit{now} velocity).452 However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date. 453 This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). 449 454 450 455 %-------------------------------------------------------------------------------------------------------------- 451 456 % Flux form Advection term 452 457 %-------------------------------------------------------------------------------------------------------------- 453 \subsection {Flux form advection term (\protect\mdl{dynadv})}458 \subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} 454 459 \label{subsec:DYN_adv_flux} 455 460 456 461 The discrete expression of the advection term is given by: 457 462 \[ 458 % \label{eq: dynadv}463 % \label{eq:DYN_adv} 459 464 \left\{ 460 465 \begin{aligned} … … 475 480 a $2^{nd}$ order centered finite difference scheme, CEN2, 476 481 or a $3^{rd}$ order upstream biased scheme, UBS. 477 The latter is described in \citet{ Shchepetkin_McWilliams_OM05}.478 The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. 482 The latter is described in \citet{shchepetkin.mcwilliams_OM05}. 483 The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. 479 484 In flux form, the schemes differ by the choice of a space and time interpolation to define the value of 480 $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie at the $T$-, $f$-,481 and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. 485 $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-, 486 and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. 482 487 483 488 %------------------------------------------------------------- 484 489 % 2nd order centred scheme 485 490 %------------------------------------------------------------- 486 \subsubsection {CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})}491 \subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 487 492 \label{subsec:DYN_adv_cen2} 488 493 489 494 In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points: 490 495 \begin{equation} 491 \label{eq: dynadv_cen2}496 \label{eq:DYN_adv_cen2} 492 497 \left\{ 493 498 \begin{aligned} … … 496 501 \end{aligned} 497 502 \right. 498 \end{equation} 499 500 The scheme is non diffusive (\ie conserves the kinetic energy) but dispersive (\ieit may create false extrema).503 \end{equation} 504 505 The scheme is non diffusive (\ie\ conserves the kinetic energy) but dispersive (\ie\ it may create false extrema). 501 506 It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 502 507 produce a sensible solution. … … 507 512 % UBS scheme 508 513 %------------------------------------------------------------- 509 \subsubsection {UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})}514 \subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})]{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 510 515 \label{subsec:DYN_adv_ubs} 511 516 … … 514 519 For example, the evaluation of $u_T^{ubs} $ is done as follows: 515 520 \begin{equation} 516 \label{eq: dynadv_ubs}521 \label{eq:DYN_adv_ubs} 517 522 u_T^{ubs} =\overline u ^i-\;\frac{1}{6} 518 523 \begin{cases} … … 522 527 \end{equation} 523 528 where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. 524 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error525 \citep{ Shchepetkin_McWilliams_OM05}.526 The overall performance of the advection scheme is similar to that reported in \citet{ Farrow1995}.529 This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error 530 \citep{shchepetkin.mcwilliams_OM05}. 531 The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}. 527 532 It is a relatively good compromise between accuracy and smoothness. 528 533 It is not a \emph{positive} scheme, meaning that false extrema are permitted. 529 534 But the amplitudes of the false extrema are significantly reduced over those in the centred second order method. 530 As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum 531 (\ie \np{ln\_dynldf\_lap}\forcode{ = }\np{ln\_dynldf\_bilap}\forcode{ =.false.}),535 As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum 536 (\ie\ \np{ln\_dynldf\_lap}\forcode{=}\np{ln\_dynldf\_bilap}\forcode{=.false.}), 532 537 and it is recommended to do so. 533 538 534 539 The UBS scheme is not used in all directions. 535 In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie $u_{uw}^{ubs}$ and536 $u_{vw}^{ubs}$ in \autoref{eq: dynadv_cen2} are used.537 UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm pursue the 540 In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and 541 $u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used. 542 UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm pursue the 538 543 sentence:Since vertical mixing of momentum is a source term of the TKE equation... } 539 544 540 For stability reasons, the first term in (\autoref{eq: dynadv_ubs}),545 For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}), 541 546 which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time), 542 547 while the second term, which is the diffusion part of the scheme, 543 548 is evaluated using the \textit{before} velocity (forward in time). 544 This is discussed by \citet{ Webb_al_JAOT98} in the context of the Quick advection scheme.549 This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 545 550 546 551 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by 547 552 one coefficient. 548 Replacing $1/6$ by $1/8$ in (\autoref{eq: dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.553 Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 549 554 This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. 550 555 Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. … … 560 565 % Hydrostatic pressure gradient term 561 566 % ================================================================ 562 \section {Hydrostatic pressure gradient (\protect\mdl{dynhpg})}567 \section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 563 568 \label{sec:DYN_hpg} 564 569 %------------------------------------------nam_dynhpg--------------------------------------------------- 565 570 566 \nlst{namdyn_hpg} 571 \begin{listing} 572 \nlst{namdyn_hpg} 573 \caption{\texttt{namdyn\_hpg}} 574 \label{lst:namdyn_hpg} 575 \end{listing} 567 576 %------------------------------------------------------------------------------------------------------------- 568 577 569 Options are defined through the \n gn{namdyn\_hpg} namelist variables.578 Options are defined through the \nam{dyn\_hpg} namelist variables. 570 579 The key distinction between the different algorithms used for 571 580 the hydrostatic pressure gradient is the vertical coordinate used, 572 since HPG is a \emph{horizontal} pressure gradient, \ie computed along geopotential surfaces.581 since HPG is a \emph{horizontal} pressure gradient, \ie\ computed along geopotential surfaces. 573 582 As a result, any tilt of the surface of the computational levels will require a specific treatment to 574 583 compute the hydrostatic pressure gradient. 575 584 576 585 The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, 577 \ie the density appearing in its expression is centred in time (\emph{now} $\rho$),586 \ie\ the density appearing in its expression is centred in time (\emph{now} $\rho$), 578 587 or a semi-implcit scheme. 579 588 At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied. … … 582 591 % z-coordinate with full step 583 592 %-------------------------------------------------------------------------------------------------------------- 584 \subsection {Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})}593 \subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})]{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} 585 594 \label{subsec:DYN_hpg_zco} 586 595 … … 592 601 for $k=km$ (surface layer, $jk=1$ in the code) 593 602 \begin{equation} 594 \label{eq: dynhpg_zco_surf}603 \label{eq:DYN_hpg_zco_surf} 595 604 \left\{ 596 605 \begin{aligned} … … 601 610 \end{aligned} 602 611 \right. 603 \end{equation} 612 \end{equation} 604 613 605 614 for $1<k<km$ (interior layer) 606 615 \begin{equation} 607 \label{eq: dynhpg_zco}616 \label{eq:DYN_hpg_zco} 608 617 \left\{ 609 618 \begin{aligned} … … 616 625 \end{aligned} 617 626 \right. 618 \end{equation} 619 620 Note that the $1/2$ factor in (\autoref{eq: dynhpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as627 \end{equation} 628 629 Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as 621 630 the vertical derivative of the scale factor at the surface level ($z=0$). 622 Note also that in case of variable volume level (\ key{vvl} defined),623 the surface pressure gradient is included in \autoref{eq: dynhpg_zco_surf} and624 \autoref{eq: dynhpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.631 Note also that in case of variable volume level (\texttt{vvl?} defined), 632 the surface pressure gradient is included in \autoref{eq:DYN_hpg_zco_surf} and 633 \autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$. 625 634 626 635 %-------------------------------------------------------------------------------------------------------------- 627 636 % z-coordinate with partial step 628 637 %-------------------------------------------------------------------------------------------------------------- 629 \subsection {Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})}638 \subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})]{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} 630 639 \label{subsec:DYN_hpg_zps} 631 640 … … 633 642 Before taking horizontal gradients between these tracer points, 634 643 a linear interpolation is used to approximate the deeper tracer as if 635 it actually lived at the depth of the shallower tracer point. 644 it actually lived at the depth of the shallower tracer point. 636 645 637 646 Apart from this modification, … … 652 661 653 662 Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers 654 (\eg, \citet{ Song1998, Shchepetkin_McWilliams_OM05}).663 (\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}). 655 664 A number of different pressure gradient options are coded but the ROMS-like, 656 665 density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation. 657 666 658 $\bullet$ Traditional coding (see for example \citet{ Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ =.true.})659 \begin{equation} 660 \label{eq: dynhpg_sco}667 $\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{=.true.}) 668 \begin{equation} 669 \label{eq:DYN_hpg_sco} 661 670 \left\{ 662 671 \begin{aligned} … … 667 676 \end{aligned} 668 677 \right. 669 \end{equation} 678 \end{equation} 670 679 671 680 Where the first term is the pressure gradient along coordinates, 672 computed as in \autoref{eq: dynhpg_zco_surf} - \autoref{eq:dynhpg_zco},673 and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 681 computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco}, 682 and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 674 683 ($e_{3w}$). 675 676 $\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}\forcode{ =.true.}).677 This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}\forcode{ =.true.}).678 679 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ =.true.})680 681 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{ Shchepetkin_McWilliams_OM05}682 (\np{ln\_dynhpg\_djc}\forcode{ =.true.}) (currently disabled; under development)683 684 Note that expression \autoref{eq: dynhpg_sco} is commonly used when the variable volume formulation is activated685 (\ key{vvl}) because in that case, even with a flat bottom,686 the coordinate surfaces are not horizontal but follow the free surface \citep{ Levier2007}.687 The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ =.true.}) is available as688 an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active.684 685 $\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}\forcode{=.true.}). 686 This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}\forcode{=.true.}). 687 688 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{=.true.}) 689 690 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05} 691 (\np{ln\_dynhpg\_djc}\forcode{=.true.}) (currently disabled; under development) 692 693 Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated 694 (\texttt{vvl?}) because in that case, even with a flat bottom, 695 the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}. 696 The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{=.true.}) is available as 697 an improved option to \np{ln\_dynhpg\_sco}\forcode{=.true.} when \texttt{vvl?} is active. 689 698 The pressure Jacobian scheme uses a constrained cubic spline to 690 699 reconstruct the density profile across the water column. … … 696 705 \subsection{Ice shelf cavity} 697 706 \label{subsec:DYN_hpg_isf} 707 698 708 Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and 699 the pressure gradient due to the ocean load (\np{ln\_dynhpg\_isf}\forcode{ =.true.}).\\709 the pressure gradient due to the ocean load (\np{ln\_dynhpg\_isf}\forcode{=.true.}).\\ 700 710 701 711 The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium. … … 704 714 corresponds to the water replaced by the ice shelf. 705 715 This top pressure is constant over time. 706 A detailed description of this method is described in \citet{ Losch2008}.\\707 708 The pressure gradient due to ocean load is computed using the expression \autoref{eq: dynhpg_sco} described in709 \autoref{subsec:DYN_hpg_sco}. 716 A detailed description of this method is described in \citet{losch_JGR08}.\\ 717 718 The pressure gradient due to ocean load is computed using the expression \autoref{eq:DYN_hpg_sco} described in 719 \autoref{subsec:DYN_hpg_sco}. 710 720 711 721 %-------------------------------------------------------------------------------------------------------------- 712 722 % Time-scheme 713 723 %-------------------------------------------------------------------------------------------------------------- 714 \subsection {Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})}724 \subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})]{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)} 715 725 \label{subsec:DYN_hpg_imp} 716 726 … … 722 732 the physical phenomenon that controls the time-step is internal gravity waves (IGWs). 723 733 A semi-implicit scheme for doubling the stability limit associated with IGWs can be used 724 \citep{ Brown_Campana_MWR78, Maltrud1998}.734 \citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}. 725 735 It involves the evaluation of the hydrostatic pressure gradient as 726 736 an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ 727 (\ie \textit{before}, \textit{now} and \textit{after} time-steps),728 rather than at the central time level $t$ only, as in the standard leapfrog scheme. 729 730 $\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{ =.true.}):731 732 \begin{equation} 733 \label{eq: dynhpg_lf}737 (\ie\ \textit{before}, \textit{now} and \textit{after} time-steps), 738 rather than at the central time level $t$ only, as in the standard leapfrog scheme. 739 740 $\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{=.true.}): 741 742 \begin{equation} 743 \label{eq:DYN_hpg_lf} 734 744 \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 735 745 -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right] 736 746 \end{equation} 737 747 738 $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{ =.true.}):739 \begin{equation} 740 \label{eq: dynhpg_imp}748 $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{=.true.}): 749 \begin{equation} 750 \label{eq:DYN_hpg_imp} 741 751 \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; 742 752 -\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] 743 753 \end{equation} 744 754 745 The semi-implicit time scheme \autoref{eq: dynhpg_imp} is made possible without755 The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without 746 756 significant additional computation since the density can be updated to time level $t+\rdt$ before 747 757 computing the horizontal hydrostatic pressure gradient. 748 758 It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using 749 \autoref{eq: dynhpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:dynhpg_lf}.750 Note that \autoref{eq: dynhpg_imp} is equivalent to applying a time filter to the pressure gradient to759 \autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}. 760 Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to 751 761 eliminate high frequency IGWs. 752 Obviously, when using \autoref{eq: dynhpg_imp},762 Obviously, when using \autoref{eq:DYN_hpg_imp}, 753 763 the doubling of the time-step is achievable only if no other factors control the time-step, 754 764 such as the stability limits associated with advection or diffusion. 755 765 756 In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}\forcode{ =.true.}.766 In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}\forcode{=.true.}. 757 767 In this case, we choose to apply the time filter to temperature and salinity used in the equation of state, 758 768 instead of applying it to the hydrostatic pressure or to the density, … … 760 770 The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows: 761 771 \[ 762 % \label{eq: rho_flt}772 % \label{eq:DYN_rho_flt} 763 773 \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) 764 774 \quad \text{with} \quad … … 773 783 % Surface Pressure Gradient 774 784 % ================================================================ 775 \section {Surface pressure gradient (\protect\mdl{dynspg})}785 \section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})} 776 786 \label{sec:DYN_spg} 777 787 %-----------------------------------------nam_dynspg---------------------------------------------------- 778 788 779 \nlst{namdyn_spg} 789 \begin{listing} 790 \nlst{namdyn_spg} 791 \caption{\texttt{namdyn\_spg}} 792 \label{lst:namdyn_spg} 793 \end{listing} 780 794 %------------------------------------------------------------------------------------------------------------ 781 795 782 Options are defined through the \n gn{namdyn\_spg} namelist variables.783 The surface pressure gradient term is related to the representation of the free surface (\autoref{sec: PE_hor_pg}).796 Options are defined through the \nam{dyn\_spg} namelist variables. 797 The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}). 784 798 The main distinction is between the fixed volume case (linear free surface) and 785 the variable volume case (nonlinear free surface, \ key{vvl} is defined).786 In the linear free surface case (\autoref{subsec: PE_free_surface})799 the variable volume case (nonlinear free surface, \texttt{vvl?} is defined). 800 In the linear free surface case (\autoref{subsec:MB_free_surface}) 787 801 the vertical scale factors $e_{3}$ are fixed in time, 788 while they are time-dependent in the nonlinear case (\autoref{subsec: PE_free_surface}).789 With both linear and nonlinear free surface, external gravity waves are allowed in the equations, 802 while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}). 803 With both linear and nonlinear free surface, external gravity waves are allowed in the equations, 790 804 which imposes a very small time step when an explicit time stepping is used. 791 Two methods are proposed to allow a longer time step for the three-dimensional equations: 792 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq: PE_flt}),805 Two methods are proposed to allow a longer time step for the three-dimensional equations: 806 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:MB_flt?}), 793 807 and the split-explicit free surface described below. 794 The extra term introduced in the filtered method is calculated implicitly, 808 The extra term introduced in the filtered method is calculated implicitly, 795 809 so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 796 810 797 811 798 812 The form of the surface pressure gradient term depends on how the user wants to 799 handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec: PE_hor_pg}).813 handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}). 800 814 Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): 801 815 an explicit formulation which requires a small time step; 802 816 a filtered free surface formulation which allows a larger time step by 803 adding a filtering term into the momentum equation; 817 adding a filtering term into the momentum equation; 804 818 and a split-explicit free surface formulation, described below, which also allows a larger time step. 805 819 … … 811 825 % Explicit free surface formulation 812 826 %-------------------------------------------------------------------------------------------------------------- 813 \subsection {Explicit free surface (\protect\key{dynspg\_exp})}827 \subsection[Explicit free surface (\texttt{ln\_dynspg\_exp}\forcode{ = .true.})]{Explicit free surface (\protect\np{ln\_dynspg\_exp}\forcode{ = .true.})} 814 828 \label{subsec:DYN_spg_exp} 815 829 816 In the explicit free surface formulation (\ key{dynspg\_exp} defined),830 In the explicit free surface formulation (\np{ln\_dynspg\_exp} set to true), 817 831 the model time step is chosen to be small enough to resolve the external gravity waves 818 832 (typically a few tens of seconds). 819 The surface pressure gradient, evaluated using a leap-frog scheme (\ie centered in time),833 The surface pressure gradient, evaluated using a leap-frog scheme (\ie\ centered in time), 820 834 is thus simply given by : 821 835 \begin{equation} 822 \label{eq: dynspg_exp}836 \label{eq:DYN_spg_exp} 823 837 \left\{ 824 838 \begin{aligned} … … 827 841 \end{aligned} 828 842 \right. 829 \end{equation} 830 831 Note that in the non-linear free surface case (\ie \key{vvl} defined),843 \end{equation} 844 845 Note that in the non-linear free surface case (\ie\ \texttt{vvl?} defined), 832 846 the surface pressure gradient is already included in the momentum tendency through 833 847 the level thickness variation allowed in the computation of the hydrostatic pressure gradient. … … 837 851 % Split-explict free surface formulation 838 852 %-------------------------------------------------------------------------------------------------------------- 839 \subsection {Split-explicit free surface (\protect\key{dynspg\_ts})}853 \subsection[Split-explicit free surface (\texttt{ln\_dynspg\_ts}\forcode{ = .true.})]{Split-explicit free surface (\protect\np{ln\_dynspg\_ts}\forcode{ = .true.})} 840 854 \label{subsec:DYN_spg_ts} 841 855 %------------------------------------------namsplit----------------------------------------------------------- … … 844 858 %------------------------------------------------------------------------------------------------------------- 845 859 846 The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined),847 also called the time-splitting formulation, follows the one proposed by \citet{ Shchepetkin_McWilliams_OM05}.860 The split-explicit free surface formulation used in \NEMO\ (\np{ln\_dynspg\_ts} set to true), 861 also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}. 848 862 The general idea is to solve the free surface equation and the associated barotropic velocity equations with 849 863 a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables 850 (\autoref{fig:DYN_ dynspg_ts}).864 (\autoref{fig:DYN_spg_ts}). 851 865 The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through 852 866 the \np{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$. 853 This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}\forcode{ =.true.}) considering that867 This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}\forcode{=.true.}) considering that 854 868 the stability of the barotropic system is essentially controled by external waves propagation. 855 869 Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry. … … 859 873 The barotropic mode solves the following equations: 860 874 % \begin{subequations} 861 % \label{eq: BT}862 \begin{equation} 863 \label{eq: BT_dyn}864 \frac{\partial {\ rm \overline{{\bf U}}_h} }{\partial t}=865 -f\;{\ rm {\bf k}}\times {\rm \overline{{\bf U}}_h}866 -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \ rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}}867 \end{equation} 868 \[ 869 % \label{eq: BT_ssh}870 \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\ rm{\bf \overline{U}}}_h \,} \right]+P-E875 % \label{eq:DYN_BT} 876 \begin{equation} 877 \label{eq:DYN_BT_dyn} 878 \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}= 879 -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h} 880 -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}} 881 \end{equation} 882 \[ 883 % \label{eq:DYN_BT_ssh} 884 \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E 871 885 \] 872 886 % \end{subequations} 873 where $\ rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes,887 where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes, 874 888 surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. 875 The third term on the right hand side of \autoref{eq: BT_dyn} represents the bottom stress876 (see section \autoref{sec:ZDF_ bfr}), explicitly accounted for at each barotropic iteration.889 The third term on the right hand side of \autoref{eq:DYN_BT_dyn} represents the bottom stress 890 (see section \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration. 877 891 Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm 878 detailed in \citet{ Shchepetkin_McWilliams_OM05}.879 AB3-AM4 coefficients used in \NEMO follow the second-order accurate,880 "multi-purpose" stability compromise as defined in \citet{ Shchepetkin_McWilliams_Bk08}881 (see their figure 12, lower left). 892 detailed in \citet{shchepetkin.mcwilliams_OM05}. 893 AB3-AM4 coefficients used in \NEMO\ follow the second-order accurate, 894 "multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09} 895 (see their figure 12, lower left). 882 896 883 897 %> > > > > > > > > > > > > > > > > > > > > > > > > > > > 884 898 \begin{figure}[!t] 885 \begin{center} 886 \includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts} 887 \caption{ 888 \protect\label{fig:DYN_dynspg_ts} 889 Schematic of the split-explicit time stepping scheme for the external and internal modes. 890 Time increases to the right. In this particular exemple, 891 a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$. 892 Internal mode time steps (which are also the model time steps) are denoted by $t-\rdt$, $t$ and $t+\rdt$. 893 Variables with $k$ superscript refer to instantaneous barotropic variables, 894 $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary (red vertical bars) and 895 secondary weights (blue vertical bars). 896 The former are used to obtain time filtered quantities at $t+\rdt$ while 897 the latter are used to obtain time averaged transports to advect tracers. 898 a) Forward time integration: \protect\np{ln\_bt\_fw}\forcode{ = .true.}, 899 \protect\np{ln\_bt\_av}\forcode{ = .true.}. 900 b) Centred time integration: \protect\np{ln\_bt\_fw}\forcode{ = .false.}, 901 \protect\np{ln\_bt\_av}\forcode{ = .true.}. 902 c) Forward time integration with no time filtering (POM-like scheme): 903 \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .false.}. 904 } 905 \end{center} 899 \centering 900 \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts} 901 \caption[Split-explicit time stepping scheme for the external and internal modes]{ 902 Schematic of the split-explicit time stepping scheme for the external and internal modes. 903 Time increases to the right. 904 In this particular exemple, 905 a boxcar averaging window over \np{nn\_baro} barotropic time steps is used 906 (\np{nn\_bt\_flt}\forcode{=1}) and \np{nn\_baro}\forcode{=5}. 907 Internal mode time steps (which are also the model time steps) are denoted by 908 $t-\rdt$, $t$ and $t+\rdt$. 909 Variables with $k$ superscript refer to instantaneous barotropic variables, 910 $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary 911 (red vertical bars) and secondary weights (blue vertical bars). 912 The former are used to obtain time filtered quantities at $t+\rdt$ while 913 the latter are used to obtain time averaged transports to advect tracers. 914 a) Forward time integration: 915 \protect\np{ln\_bt\_fw}\forcode{=.true.}, \protect\np{ln\_bt\_av}\forcode{=.true.}. 916 b) Centred time integration: 917 \protect\np{ln\_bt\_fw}\forcode{=.false.}, \protect\np{ln\_bt\_av}\forcode{=.true.}. 918 c) Forward time integration with no time filtering (POM-like scheme): 919 \protect\np{ln\_bt\_fw}\forcode{=.true.}, \protect\np{ln\_bt\_av}\forcode{=.false.}.} 920 \label{fig:DYN_spg_ts} 906 921 \end{figure} 907 922 %> > > > > > > > > > > > > > > > > > > > > > > > > > > > 908 923 909 In the default case (\np{ln\_bt\_fw}\forcode{ =.true.}),924 In the default case (\np{ln\_bt\_fw}\forcode{=.true.}), 910 925 the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps 911 (\autoref{fig:DYN_ dynspg_ts}a).926 (\autoref{fig:DYN_spg_ts}a). 912 927 To avoid aliasing of fast barotropic motions into three dimensional equations, 913 time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{ =.true.}).928 time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{=.true.}). 914 929 In that case, the integration is extended slightly beyond \textit{after} time step to 915 930 provide time filtered quantities. 916 931 These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. 917 Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, 932 Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, 918 933 asselin filtering is not applied to barotropic quantities.\\ 919 934 Alternatively, one can choose to integrate barotropic equations starting from \textit{before} time step 920 (\np{ln\_bt\_fw}\forcode{ =.false.}).935 (\np{ln\_bt\_fw}\forcode{=.false.}). 921 936 Although more computationaly expensive ( \np{nn\_baro} additional iterations are indeed necessary), 922 937 the baroclinic to barotropic forcing term given at \textit{now} time step become centred in … … 936 951 and time splitting not compatible. 937 952 Advective barotropic velocities are obtained by using a secondary set of filtering weights, 938 uniquely defined from the filter coefficients used for the time averaging (\citet{ Shchepetkin_McWilliams_OM05}).953 uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}). 939 954 Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to 940 955 obtain exact conservation. … … 943 958 944 959 One can eventually choose to feedback instantaneous values by not using any time filter 945 (\np{ln\_bt\_av}\forcode{ = .false.}).960 (\np{ln\_bt\_av}\forcode{=.false.}). 946 961 In that case, external mode equations are continuous in time, 947 \ie they are not re-initialized when starting a new sub-stepping sequence.962 \ie\ they are not re-initialized when starting a new sub-stepping sequence. 948 963 This is the method used so far in the POM model, the stability being maintained by 949 964 refreshing at (almost) each barotropic time step advection and horizontal diffusion terms. 950 Since the latter terms have not been added in \NEMO for computational efficiency,965 Since the latter terms have not been added in \NEMO\ for computational efficiency, 951 966 removing time filtering is not recommended except for debugging purposes. 952 967 This may be used for instance to appreciate the damping effect of the standard formulation on 953 968 external gravity waves in idealized or weakly non-linear cases. 954 969 Although the damping is lower than for the filtered free surface, 955 it is still significant as shown by \citet{ Levier2007} in the case of an analytical barotropic Kelvin wave.970 it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave. 956 971 957 972 %>>>>>=============== 958 \gmcomment{ %%% copy from griffies Book 973 \gmcomment{ %%% copy from griffies Book 959 974 960 975 \textbf{title: Time stepping the barotropic system } … … 963 978 Hence, we can update the surface height and vertically integrated velocity with a leap-frog scheme using 964 979 the small barotropic time step $\rdt$. 965 We have 980 We have 966 981 967 982 \[ … … 986 1001 and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over 987 1002 a single cycle. 988 This is also the time that sets the barotropic time steps via 1003 This is also the time that sets the barotropic time steps via 989 1004 \[ 990 1005 % \label{eq:DYN_spg_ts_t} … … 992 1007 \] 993 1008 with $n$ an integer. 994 The density scaled surface pressure is evaluated via 1009 The density scaled surface pressure is evaluated via 995 1010 \[ 996 1011 % \label{eq:DYN_spg_ts_ps} … … 1001 1016 \end{cases} 1002 1017 \] 1003 To get started, we assume the following initial conditions 1018 To get started, we assume the following initial conditions 1004 1019 \[ 1005 1020 % \label{eq:DYN_spg_ts_eta} … … 1009 1024 \end{split} 1010 1025 \] 1011 with 1026 with 1012 1027 \[ 1013 1028 % \label{eq:DYN_spg_ts_etaF} … … 1015 1030 \] 1016 1031 the time averaged surface height taken from the previous barotropic cycle. 1017 Likewise, 1032 Likewise, 1018 1033 \[ 1019 1034 % \label{eq:DYN_spg_ts_u} … … 1021 1036 \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} 1022 1037 \] 1023 with 1038 with 1024 1039 \[ 1025 1040 % \label{eq:DYN_spg_ts_u} … … 1027 1042 \] 1028 1043 the time averaged vertically integrated transport. 1029 Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration. 1044 Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration. 1030 1045 1031 1046 Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ , 1032 1047 the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at 1033 baroclinic time $\tau + \rdt \tau$ 1048 baroclinic time $\tau + \rdt \tau$ 1034 1049 \[ 1035 1050 % \label{eq:DYN_spg_ts_u} … … 1037 1052 \] 1038 1053 The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using 1039 the following form 1054 the following form 1040 1055 1041 1056 \begin{equation} 1042 1057 \label{eq:DYN_spg_ts_ssh} 1043 \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] 1058 \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] 1044 1059 \end{equation} 1045 1060 1046 1061 The use of this "big-leap-frog" scheme for the surface height ensures compatibility between 1047 1062 the mass/volume budgets and the tracer budgets. 1048 More discussion of this point is provided in Chapter 10 (see in particular Section 10.2). 1049 1063 More discussion of this point is provided in Chapter 10 (see in particular Section 10.2). 1064 1050 1065 In general, some form of time filter is needed to maintain integrity of the surface height field due to 1051 1066 the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. 1052 1067 We have tried various forms of such filtering, 1053 with the following method discussed in \cite{ Griffies_al_MWR01} chosen due to1068 with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to 1054 1069 its stability and reasonably good maintenance of tracer conservation properties (see ??). 1055 1070 … … 1058 1073 \eta^{F}(\tau-\Delta) = \overline{\eta^{(b)}(\tau)} 1059 1074 \end{equation} 1060 Another approach tried was 1075 Another approach tried was 1061 1076 1062 1077 \[ … … 1071 1086 eliminating tracer and surface height time filtering (see ?? for more complete discussion). 1072 1087 However, in the general case with a non-zero $\alpha$, 1073 the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended. 1088 the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended. 1074 1089 1075 1090 } %%end gm comment (copy of griffies book) … … 1081 1096 % Filtered free surface formulation 1082 1097 %-------------------------------------------------------------------------------------------------------------- 1083 \subsection {Filtered free surface (\protect\key{dynspg\_flt})}1098 \subsection[Filtered free surface (\texttt{dynspg\_flt?})]{Filtered free surface (\protect\texttt{dynspg\_flt?})} 1084 1099 \label{subsec:DYN_spg_fltp} 1085 1100 1086 The filtered formulation follows the \citet{ Roullet_Madec_JGR00} implementation.1087 The extra term introduced in the equations (see \autoref{subsec: PE_free_surface}) is solved implicitly.1101 The filtered formulation follows the \citet{roullet.madec_JGR00} implementation. 1102 The extra term introduced in the equations (see \autoref{subsec:MB_free_surface}) is solved implicitly. 1088 1103 The elliptic solvers available in the code are documented in \autoref{chap:MISC}. 1089 1104 1090 1105 %% gm %%======>>>> given here the discrete eqs provided to the solver 1091 \gmcomment{ %%% copy from chap-model basics 1106 \gmcomment{ %%% copy from chap-model basics 1092 1107 \[ 1093 % \label{eq: spg_flt}1094 \frac{\partial {\ rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}1108 % \label{eq:DYN_spg_flt} 1109 \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}} 1095 1110 - g \nabla \left( \tilde{\rho} \ \eta \right) 1096 1111 - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) … … 1098 1113 where $T_c$, is a parameter with dimensions of time which characterizes the force, 1099 1114 $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, 1100 and $\ rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,1101 non-linear and viscous terms in \autoref{eq: PE_dyn}.1115 and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 1116 non-linear and viscous terms in \autoref{eq:MB_dyn}. 1102 1117 } %end gmcomment 1103 1118 1104 Note that in the linear free surface formulation (\ key{vvl} not defined),1119 Note that in the linear free surface formulation (\texttt{vvl?} not defined), 1105 1120 the ocean depth is time-independent and so is the matrix to be inverted. 1106 It is computed once and for all and applies to all ocean time steps. 1121 It is computed once and for all and applies to all ocean time steps. 1107 1122 1108 1123 % ================================================================ 1109 1124 % Lateral diffusion term 1110 1125 % ================================================================ 1111 \section {Lateral diffusion term and operators (\protect\mdl{dynldf})}1126 \section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})} 1112 1127 \label{sec:DYN_ldf} 1113 1128 %------------------------------------------nam_dynldf---------------------------------------------------- 1114 1129 1115 \nlst{namdyn_ldf} 1130 \begin{listing} 1131 \nlst{namdyn_ldf} 1132 \caption{\texttt{namdyn\_ldf}} 1133 \label{lst:namdyn_ldf} 1134 \end{listing} 1116 1135 %------------------------------------------------------------------------------------------------------------- 1117 1136 1118 Options are defined through the \n gn{namdyn\_ldf} namelist variables.1137 Options are defined through the \nam{dyn\_ldf} namelist variables. 1119 1138 The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators. 1120 1139 The coefficients may be constant or spatially variable; 1121 1140 the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}). 1122 1141 The lateral diffusion of momentum is evaluated using a forward scheme, 1123 \ie the velocity appearing in its expression is the \textit{before} velocity in time,1142 \ie\ the velocity appearing in its expression is the \textit{before} velocity in time, 1124 1143 except for the pure vertical component that appears when a tensor of rotation is used. 1125 This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap: STP}).1144 This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}). 1126 1145 1127 1146 At the lateral boundaries either free slip, … … 1137 1156 In finite difference methods, 1138 1157 the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since 1139 its damping time (\ie its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$1158 its damping time (\ie\ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$ 1140 1159 (so that short waves damped more rapidelly than long ones), 1141 1160 whereas the Laplace operator damping time scales only like $\lambda^{-2}$. … … 1143 1162 1144 1163 % ================================================================ 1145 \subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})] 1146 {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 1164 \subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})]{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} 1147 1165 \label{subsec:DYN_ldf_lap} 1148 1166 1149 For lateral iso-level diffusion, the discrete operator is: 1150 \begin{equation} 1151 \label{eq: dynldf_lap}1167 For lateral iso-level diffusion, the discrete operator is: 1168 \begin{equation} 1169 \label{eq:DYN_ldf_lap} 1152 1170 \left\{ 1153 1171 \begin{aligned} 1154 D_u^{l{\ rm {\bf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}1155 \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[ 1172 D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} 1173 \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[ 1156 1174 {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\ 1157 D_v^{l{\ rm {\bf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}1158 \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[ 1175 D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} 1176 \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[ 1159 1177 {A_f^{lm} \;e_{3f} \zeta } \right] 1160 1178 \end{aligned} 1161 1179 \right. 1162 \end{equation} 1163 1164 As explained in \autoref{subsec: PE_ldf},1180 \end{equation} 1181 1182 As explained in \autoref{subsec:MB_ldf}, 1165 1183 this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and 1166 ensures a complete separation between the vorticity and divergence parts of the momentum diffusion. 1184 ensures a complete separation between the vorticity and divergence parts of the momentum diffusion. 1167 1185 1168 1186 %-------------------------------------------------------------------------------------------------------------- 1169 1187 % Rotated laplacian operator 1170 1188 %-------------------------------------------------------------------------------------------------------------- 1171 \subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})] 1172 {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 1189 \subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})]{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} 1173 1190 \label{subsec:DYN_ldf_iso} 1174 1191 1175 1192 A rotation of the lateral momentum diffusion operator is needed in several cases: 1176 for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}\forcode{ =.true.}) and1177 for either iso-neutral (\np{ln\_dynldf\_iso}\forcode{ =.true.}) or1178 geopotential (\np{ln\_dynldf\_hor}\forcode{ =.true.}) diffusion in the $s$-coordinate.1193 for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}\forcode{=.true.}) and 1194 for either iso-neutral (\np{ln\_dynldf\_iso}\forcode{=.true.}) or 1195 geopotential (\np{ln\_dynldf\_hor}\forcode{=.true.}) diffusion in the $s$-coordinate. 1179 1196 In the partial step case, coordinates are horizontal except at the deepest level and 1180 no rotation is performed when \np{ln\_dynldf\_hor}\forcode{ =.true.}.1197 no rotation is performed when \np{ln\_dynldf\_hor}\forcode{=.true.}. 1181 1198 The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on 1182 1199 each momentum component. … … 1184 1201 The resulting discrete representation is: 1185 1202 \begin{equation} 1186 \label{eq: dyn_ldf_iso}1203 \label{eq:DYN_ldf_iso} 1187 1204 \begin{split} 1188 1205 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ … … 1208 1225 -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}} 1209 1226 \right)} \right]} \right. \\ 1210 & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 1227 & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 1211 1228 }\,\delta_{j} [v] - e_{1t}\, r_{2t} 1212 1229 \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}} 1213 1230 \right)} \right] \\ 1214 & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 1231 & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 1215 1232 {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\ 1216 1233 & \ \qquad \qquad \qquad \quad\ 1217 1234 - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\ 1218 & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 1235 & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 1219 1236 +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} 1220 \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 1237 \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 1221 1238 \end{split} 1222 1239 \end{equation} 1223 1240 where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and 1224 the surface of computation ($z$- or $s$-surfaces). 1241 the surface of computation ($z$- or $s$-surfaces). 1225 1242 The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}). 1226 1243 … … 1228 1245 % Iso-level bilaplacian operator 1229 1246 %-------------------------------------------------------------------------------------------------------------- 1230 \subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})] 1231 {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 1247 \subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})]{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} 1232 1248 \label{subsec:DYN_ldf_bilap} 1233 1249 1234 The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq: dynldf_lap} twice.1250 The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice. 1235 1251 It requires an additional assumption on boundary conditions: 1236 1252 the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen, … … 1243 1259 % Vertical diffusion term 1244 1260 % ================================================================ 1245 \section {Vertical diffusion term (\protect\mdl{dynzdf})}1261 \section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})} 1246 1262 \label{sec:DYN_zdf} 1247 1263 %----------------------------------------------namzdf------------------------------------------------------ 1248 1264 1249 \nlst{namzdf}1250 1265 %------------------------------------------------------------------------------------------------------------- 1251 1266 1252 Options are defined through the \n gn{namzdf} namelist variables.1267 Options are defined through the \nam{zdf} namelist variables. 1253 1268 The large vertical diffusion coefficient found in the surface mixed layer together with high vertical resolution implies that in the case of explicit time stepping there would be too restrictive a constraint on the time step. 1254 1269 Two time stepping schemes can be used for the vertical diffusion term: 1255 1270 $(a)$ a forward time differencing scheme 1256 (\np{ln\_zdfexp}\forcode{ =.true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or1257 $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{ =.false.})1258 (see \autoref{chap: STP}).1259 Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 1271 (\np{ln\_zdfexp}\forcode{=.true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or 1272 $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{=.false.}) 1273 (see \autoref{chap:TD}). 1274 Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 1260 1275 1261 1276 The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is. 1262 The vertical diffusion operators given by \autoref{eq: PE_zdf} take the following semi-discrete space form:1263 \[ 1264 % \label{eq: dynzdf}1277 The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form: 1278 \[ 1279 % \label{eq:DYN_zdf} 1265 1280 \left\{ 1266 1281 \begin{aligned} … … 1280 1295 the vertical turbulent momentum fluxes, 1281 1296 \begin{equation} 1282 \label{eq: dynzdf_sbc}1297 \label{eq:DYN_zdf_sbc} 1283 1298 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 1284 1299 = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v } … … 1286 1301 where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in 1287 1302 the (\textbf{i},\textbf{j}) coordinate system. 1288 The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in 1303 The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in 1289 1304 the vertical over the mixed layer depth. 1290 1305 If the vertical mixing coefficient is small (when no mixed layer scheme is used) … … 1293 1308 1294 1309 The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation 1295 (see \autoref{sec:ZDF_ bfr})1310 (see \autoref{sec:ZDF_drg}) 1296 1311 1297 1312 % ================================================================ … … 1303 1318 Besides the surface and bottom stresses (see the above section) 1304 1319 which are introduced as boundary conditions on the vertical mixing, 1305 three other forcings may enter the dynamical equations by affecting the surface pressure gradient. 1306 1307 (1) When \np{ln\_apr\_dyn}\forcode{ =.true.} (see \autoref{sec:SBC_apr}),1320 three other forcings may enter the dynamical equations by affecting the surface pressure gradient. 1321 1322 (1) When \np{ln\_apr\_dyn}\forcode{=.true.} (see \autoref{sec:SBC_apr}), 1308 1323 the atmospheric pressure is taken into account when computing the surface pressure gradient. 1309 1324 1310 (2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ =.true.} (see \autoref{sec:SBC_tide}),1325 (2) When \np{ln\_tide\_pot}\forcode{=.true.} and \np{ln\_tide}\forcode{=.true.} (see \autoref{sec:SBC_tide}), 1311 1326 the tidal potential is taken into account when computing the surface pressure gradient. 1312 1327 1313 (3) When \np{nn\_ice\_embd}\forcode{ =2} and LIM or CICE is used1314 (\ie when the sea-ice is embedded in the ocean),1328 (3) When \np{nn\_ice\_embd}\forcode{=2} and LIM or CICE is used 1329 (\ie\ when the sea-ice is embedded in the ocean), 1315 1330 the snow-ice mass is taken into account when computing the surface pressure gradient. 1316 1331 … … 1320 1335 1321 1336 % ================================================================ 1322 % Wetting and drying 1337 % Wetting and drying 1323 1338 % ================================================================ 1324 1339 \section{Wetting and drying } 1325 1340 \label{sec:DYN_wetdry} 1341 1326 1342 There are two main options for wetting and drying code (wd): 1327 1343 (a) an iterative limiter (il) and (b) a directional limiter (dl). 1328 The directional limiter is based on the scheme developed by \cite{ WarnerEtal13} for RO1344 The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO 1329 1345 MS 1330 which was in turn based on ideas developed for POM by \cite{ Oey06}. The iterative1346 which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative 1331 1347 limiter is a new scheme. The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 1332 1348 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated 1333 1349 by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$. 1334 1350 1335 \nlst{namwad} 1351 \begin{listing} 1352 \nlst{namwad} 1353 \caption{\texttt{namwad}} 1354 \label{lst:namwad} 1355 \end{listing} 1336 1356 1337 1357 The following terminology is used. The depth of the topography (positive downwards) 1338 at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the NEMOcode.1358 at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code. 1339 1359 The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign 1340 1360 conventions used, the water depth, $h$, is the height of the free surface plus the depth of the … … 1344 1364 covered by water. They require the topography specified with a model 1345 1365 configuration to have negative depths at points where the land is higher than the 1346 topography's reference sea-level. The vertical grid in NEMOis normally computed relative to an1366 topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an 1347 1367 initial state with zero sea surface height elevation. 1348 1368 The user can choose to compute the vertical grid and heights in the model relative to … … 1363 1383 All these configurations have used pure sigma coordinates. It is expected that 1364 1384 the wetting and drying code will work in domains with more general s-coordinates provided 1365 the coordinates are pure sigma in the region where wetting and drying actually occurs. 1385 the coordinates are pure sigma in the region where wetting and drying actually occurs. 1366 1386 1367 1387 The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter. … … 1372 1392 % Iterative limiters 1373 1393 %----------------------------------------------------------------------------------------- 1374 \subsection [Directional limiter (\textit{wet\_dry})] 1375 {Directional limiter (\mdl{wet\_dry})} 1394 \subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})} 1376 1395 \label{subsec:DYN_wd_directional_limiter} 1396 1377 1397 The principal idea of the directional limiter is that 1378 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than rn\_wdmin1).1398 water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn\_wdmin1}). 1379 1399 1380 1400 All the changes associated with this option are made to the barotropic solver for the non-linear … … 1386 1406 1387 1407 The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj). 1388 If the user sets ln\_wd\_dl\_ramp = .False.then zuwdmask is 1 when the1389 flux is from a cell with water depth greater than rn\_wdmin1and 0 otherwise. If the user sets1390 ln\_wd\_dl\_ramp = .True.the flux across the face is ramped down as the water depth decreases1391 from 2 * rn\_wdmin1 to rn\_wdmin1. The use of this ramp reduced grid-scale noise in idealised test cases.1408 If the user sets \np{ln\_wd\_dl\_ramp}\forcode{=.false.} then zuwdmask is 1 when the 1409 flux is from a cell with water depth greater than \np{rn\_wdmin1} and 0 otherwise. If the user sets 1410 \np{ln\_wd\_dl\_ramp}\forcode{=.true.} the flux across the face is ramped down as the water depth decreases 1411 from 2 * \np{rn\_wdmin1} to \np{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases. 1392 1412 1393 1413 At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen … … 1400 1420 1401 1421 1402 \cite{ WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic1422 \cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic 1403 1423 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 1404 1424 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer 1405 1425 fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because 1406 1426 the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts 1407 to equal their mean value during the barotropic steps. If the user sets ln\_wd\_dl\_bc = .True., the1408 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask. 1427 to equal their mean value during the barotropic steps. If the user sets \np{ln\_wd\_dl\_bc}\forcode{=.true.}, the 1428 baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask. 1409 1429 1410 1430 %----------------------------------------------------------------------------------------- … … 1412 1432 %----------------------------------------------------------------------------------------- 1413 1433 1414 \subsection [Iterative limiter (\textit{wet\_dry})] 1415 {Iterative limiter (\mdl{wet\_dry})} 1434 \subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})} 1416 1435 \label{subsec:DYN_wd_iterative_limiter} 1417 1436 1418 \subsubsection [Iterative flux limiter (\textit{wet\_dry})] 1419 {Iterative flux limiter (\mdl{wet\_dry})} 1420 \label{subsubsec:DYN_wd_il_spg_limiter} 1437 \subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})} 1438 \label{subsec:DYN_wd_il_spg_limiter} 1421 1439 1422 1440 The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry'' … … 1426 1444 1427 1445 The continuity equation for the total water depth in a column 1428 \begin{equation} \label{dyn_wd_continuity} 1429 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 1446 \begin{equation} 1447 \label{eq:DYN_wd_continuity} 1448 \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 1430 1449 \end{equation} 1431 1450 can be written in discrete form as 1432 1451 1433 \begin{align} \label{dyn_wd_continuity_2} 1434 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 1435 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 1436 &= \mathrm{zzflx}_{i,j} . 1452 \begin{align} 1453 \label{eq:DYN_wd_continuity_2} 1454 \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 1455 &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j} + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 1456 &= \mathrm{zzflx}_{i,j} . 1437 1457 \end{align} 1438 1458 … … 1447 1467 (zzflxp) and fluxes that are into the cell (zzflxn). Clearly 1448 1468 1449 \begin{equation} \label{dyn_wd_zzflx_p_n_1} 1450 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 1469 \begin{equation} 1470 \label{eq:DYN_wd_zzflx_p_n_1} 1471 \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 1451 1472 \end{equation} 1452 1473 … … 1459 1480 $\mathrm{zcoef}_{i,j}^{(m)}$ such that: 1460 1481 1461 \begin{equation} \label{dyn_wd_continuity_coef} 1462 \begin{split} 1463 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 1464 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 1465 \end{split} 1482 \begin{equation} 1483 \label{eq:DYN_wd_continuity_coef} 1484 \begin{split} 1485 \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 1486 \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 1487 \end{split} 1466 1488 \end{equation} 1467 1489 … … 1471 1493 The iteration is initialised by setting 1472 1494 1473 \begin{equation} \label{dyn_wd_zzflx_initial} 1474 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 1495 \begin{equation} 1496 \label{eq:DYN_wd_zzflx_initial} 1497 \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 1475 1498 \end{equation} 1476 1499 1477 1500 The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the 1478 1501 cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell 1479 times the timestep divided by the cell area. Using (\ ref{dyn_wd_continuity_2}) this1502 times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this 1480 1503 condition is 1481 1504 1482 \begin{equation} \label{dyn_wd_continuity_if} 1483 h_{i,j}(t_e) - \mathrm{rn\_wdmin1} < \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) . 1484 \end{equation} 1485 1486 Rearranging (\ref{dyn_wd_continuity_if}) we can obtain an expression for the maximum 1505 \begin{equation} 1506 \label{eq:DYN_wd_continuity_if} 1507 h_{i,j}(t_e) - \mathrm{rn\_wdmin1} < \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) . 1508 \end{equation} 1509 1510 Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum 1487 1511 outward flux that can be allowed and still maintain the minimum wet depth: 1488 1512 1489 \begin{equation} \label{dyn_wd_max_flux} 1490 \begin{split} 1491 \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 1492 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] 1493 \end{split} 1494 \end{equation} 1495 1496 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is 1497 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an 1513 \begin{equation} 1514 \label{eq:DYN_wd_max_flux} 1515 \begin{split} 1516 \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 1517 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] 1518 \end{split} 1519 \end{equation} 1520 1521 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is 1522 this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an 1498 1523 expression for the coefficient needed to multiply the outward flux at this cell in order 1499 1524 to avoid drying. 1500 1525 1501 \begin{equation} \label{dyn_wd_continuity_nxtcoef} 1502 \begin{split} 1503 \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 1504 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 1505 \end{split} 1526 \begin{equation} 1527 \label{eq:DYN_wd_continuity_nxtcoef} 1528 \begin{split} 1529 \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2}) \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 1530 \phantom{[} & - \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 1531 \end{split} 1506 1532 \end{equation} 1507 1533 … … 1522 1548 % Surface pressure gradients 1523 1549 %---------------------------------------------------------------------------------------- 1524 \subsubsection [Modification of surface pressure gradients (\textit{dynhpg})] 1525 {Modification of surface pressure gradients (\mdl{dynhpg})} 1526 \label{subsubsec:DYN_wd_il_spg} 1550 \subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})} 1551 \label{subsec:DYN_wd_il_spg} 1527 1552 1528 1553 At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the … … 1537 1562 neighbouring $(i+1,j)$ and $(i,j)$ tracer points. zcpx is calculated using two logicals 1538 1563 variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid 1539 column. The three possible combinations are illustrated in figure \ref{Fig_WAD_dynhpg}.1564 column. The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}. 1540 1565 1541 1566 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1542 \begin{figure}[!ht] \begin{center} 1543 \includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg} 1544 \caption{ \label{Fig_WAD_dynhpg} 1545 Illustrations of the three possible combinations of the logical variables controlling the 1546 limiting of the horizontal pressure gradient in wetting and drying regimes} 1547 \end{center}\end{figure} 1567 \begin{figure}[!ht] 1568 \centering 1569 \includegraphics[width=\textwidth]{Fig_WAD_dynhpg} 1570 \caption[Combinations controlling the limiting of the horizontal pressure gradient in 1571 wetting and drying regimes]{ 1572 Three possible combinations of the logical variables controlling the 1573 limiting of the horizontal pressure gradient in wetting and drying regimes} 1574 \label{fig:DYN_WAD_dynhpg} 1575 \end{figure} 1548 1576 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 1549 1577 … … 1553 1581 of the topography at the two points: 1554 1582 1555 \begin{equation} \label{dyn_ll_tmp1} 1556 \begin{split} 1557 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 1583 \begin{equation} 1584 \label{eq:DYN_ll_tmp1} 1585 \begin{split} 1586 \mathrm{ll\_tmp1} = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 1558 1587 & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\ .and.} \\ 1559 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\1560 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\1561 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }1562 \end{split}1588 & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\ 1589 & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\ 1590 & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 } 1591 \end{split} 1563 1592 \end{equation} 1564 1593 … … 1567 1596 at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ 1568 1597 1569 \begin{equation} \label{dyn_ll_tmp2} 1570 \begin{split} 1571 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 1572 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 1573 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 1574 \end{split} 1598 \begin{equation} 1599 \label{eq:DYN_ll_tmp2} 1600 \begin{split} 1601 \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 1602 & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 1603 & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 1604 \end{split} 1575 1605 \end{equation} 1576 1606 … … 1588 1618 conditions. 1589 1619 1590 \subsubsection [Additional considerations (\textit{usrdef\_zgr})] 1591 {Additional considerations (\mdl{usrdef\_zgr})} 1592 \label{subsubsec:WAD_additional} 1620 \subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})} 1621 \label{subsec:DYN_WAD_additional} 1593 1622 1594 1623 In the very shallow water where wetting and drying occurs the parametrisation of … … 1603 1632 % The WAD test cases 1604 1633 %---------------------------------------------------------------------------------------- 1605 \subsection [The WAD test cases (\textit{usrdef\_zgr})] 1606 {The WAD test cases (\mdl{usrdef\_zgr})} 1607 \label{WAD_test_cases} 1634 \subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})} 1635 \label{subsec:DYN_WAD_test_cases} 1608 1636 1609 1637 See the WAD tests MY\_DOC documention for details of the WAD test cases. … … 1612 1640 1613 1641 % ================================================================ 1614 % Time evolution term 1615 % ================================================================ 1616 \section {Time evolution term (\protect\mdl{dynnxt})}1642 % Time evolution term 1643 % ================================================================ 1644 \section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})} 1617 1645 \label{sec:DYN_nxt} 1618 1646 1619 1647 %----------------------------------------------namdom---------------------------------------------------- 1620 1648 1621 \nlst{namdom}1622 1649 %------------------------------------------------------------------------------------------------------------- 1623 1650 1624 Options are defined through the \n gn{namdom} namelist variables.1651 Options are defined through the \nam{dom} namelist variables. 1625 1652 The general framework for dynamics time stepping is a leap-frog scheme, 1626 \ie a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:STP}).1653 \ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}). 1627 1654 The scheme is applied to the velocity, except when 1628 1655 using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) 1629 in the variable volume case (\ key{vvl} defined),1630 where it has to be applied to the thickness weighted velocity (see \autoref{sec: A_momentum})1656 in the variable volume case (\texttt{vvl?} defined), 1657 where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum}) 1631 1658 1632 1659 $\bullet$ vector invariant form or linear free surface 1633 (\np{ln\_dynhpg\_vec}\forcode{ = .true.} ; \key{vvl} not defined):1634 \[ 1635 % \label{eq: dynnxt_vec}1660 (\np{ln\_dynhpg\_vec}\forcode{=.true.} ; \texttt{vvl?} not defined): 1661 \[ 1662 % \label{eq:DYN_nxt_vec} 1636 1663 \left\{ 1637 1664 \begin{aligned} … … 1643 1670 1644 1671 $\bullet$ flux form and nonlinear free surface 1645 (\np{ln\_dynhpg\_vec}\forcode{ = .false.} ; \key{vvl} defined):1646 \[ 1647 % \label{eq: dynnxt_flux}1672 (\np{ln\_dynhpg\_vec}\forcode{=.false.} ; \texttt{vvl?} defined): 1673 \[ 1674 % \label{eq:DYN_nxt_flux} 1648 1675 \left\{ 1649 1676 \begin{aligned} … … 1657 1684 the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient. 1658 1685 $\gamma$ is initialized as \np{nn\_atfp} (namelist parameter). 1659 Its default value is \np{nn\_atfp}\forcode{ =10.e-3}.1686 Its default value is \np{nn\_atfp}\forcode{=10.e-3}. 1660 1687 In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for 1661 1688 the momentum equations.
Note: See TracChangeset
for help on using the changeset viewer.