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