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