[10419] | 1 | \documentclass[../main/NEMO_manual]{subfiles} |
---|
| 2 | |
---|
[6997] | 3 | \begin{document} |
---|
[707] | 4 | % ================================================================ |
---|
[6140] | 5 | % Chapter ——— Ocean Dynamics (DYN) |
---|
[707] | 6 | % ================================================================ |
---|
| 7 | \chapter{Ocean Dynamics (DYN)} |
---|
[9407] | 8 | \label{chap:DYN} |
---|
[10419] | 9 | |
---|
[707] | 10 | \minitoc |
---|
| 11 | |
---|
[10368] | 12 | Using the representation described in \autoref{chap:DOM}, |
---|
| 13 | several semi-discrete space forms of the dynamical equations are available depending on |
---|
| 14 | the vertical coordinate used and on the conservation properties of the vorticity term. |
---|
| 15 | In all the equations presented here, the masking has been omitted for simplicity. |
---|
| 16 | One must be aware that all the quantities are masked fields and |
---|
| 17 | that each time an average or difference operator is used, the resulting field is multiplied by a mask. |
---|
[707] | 18 | |
---|
| 19 | The prognostic ocean dynamics equation can be summarized as follows: |
---|
[10419] | 20 | \[ |
---|
| 21 | \text{NXT} = \dbinom {\text{VOR} + \text{KEG} + \text {ZAD} } |
---|
| 22 | {\text{COR} + \text{ADV} } |
---|
| 23 | + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF} |
---|
| 24 | \] |
---|
[10368] | 25 | NXT stands for next, referring to the time-stepping. |
---|
| 26 | The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that |
---|
| 27 | are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) and |
---|
| 28 | a vertical advection part (ZAD) in the vector invariant formulation, |
---|
| 29 | or a Coriolis and advection part (COR+ADV) in the flux formulation. |
---|
| 30 | The terms following these are the pressure gradient contributions |
---|
| 31 | (HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient); |
---|
| 32 | and contributions from lateral diffusion (LDF) and vertical diffusion (ZDF), |
---|
| 33 | which are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules. |
---|
| 34 | The vertical diffusion term includes the surface and bottom stresses. |
---|
| 35 | The external forcings and parameterisations require complex inputs |
---|
| 36 | (surface wind stress calculation using bulk formulae, estimation of mixing coefficients) |
---|
| 37 | that are carried out in modules SBC, LDF and ZDF and are described in |
---|
| 38 | \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. |
---|
[707] | 39 | |
---|
[10368] | 40 | In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence, |
---|
| 41 | curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module). |
---|
[707] | 42 | |
---|
[817] | 43 | The different options available to the user are managed by namelist variables. |
---|
[2282] | 44 | For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, |
---|
[10368] | 45 | where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. |
---|
| 46 | If a CPP key is used for this term its name is \key{ttt}. |
---|
| 47 | The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory, |
---|
| 48 | and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine. |
---|
[707] | 49 | |
---|
[10368] | 50 | The user has the option of extracting and outputting each tendency term from the 3D momentum equations |
---|
| 51 | (\key{trddyn} defined), as described in \autoref{chap:MISC}. |
---|
| 52 | Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \key{trdvor} is defined) |
---|
| 53 | can be derived from the 3D terms. |
---|
[817] | 54 | %%% |
---|
[996] | 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"?} |
---|
[707] | 57 | |
---|
| 58 | % ================================================================ |
---|
[2282] | 59 | % Sea Surface Height evolution & Diagnostics variables |
---|
| 60 | % ================================================================ |
---|
| 61 | \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} |
---|
[9407] | 62 | \label{sec:DYN_divcur_wzv} |
---|
[2282] | 63 | |
---|
| 64 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 65 | % Horizontal divergence and relative vorticity |
---|
| 66 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 67 | \subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} |
---|
[9407] | 68 | \label{subsec:DYN_divcur} |
---|
[2282] | 69 | |
---|
| 70 | The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: |
---|
[10419] | 71 | \begin{equation} |
---|
| 72 | \label{eq:divcur_cur} |
---|
| 73 | \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right] |
---|
| 74 | -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) |
---|
[2282] | 75 | \end{equation} |
---|
| 76 | |
---|
[10368] | 77 | The horizontal divergence is defined at a $T$-point. |
---|
| 78 | It is given by: |
---|
[10419] | 79 | \[ |
---|
| 80 | % \label{eq:divcur_div} |
---|
| 81 | \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } |
---|
| 82 | \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right] |
---|
| 83 | +\delta_j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right) |
---|
| 84 | \] |
---|
[2282] | 85 | |
---|
[10368] | 86 | Note that although the vorticity has the same discrete expression in $z$- and $s$-coordinates, |
---|
| 87 | its physical meaning is not identical. |
---|
| 88 | $\zeta$ is a pseudo vorticity along $s$-surfaces |
---|
| 89 | (only pseudo because $(u,v)$ are still defined along geopotential surfaces, |
---|
| 90 | but are not necessarily defined at the same depth). |
---|
[2282] | 91 | |
---|
[10368] | 92 | The vorticity and divergence at the \textit{before} step are used in the computation of |
---|
| 93 | the horizontal diffusion of momentum. |
---|
| 94 | Note that because they have been calculated prior to the Asselin filtering of the \textit{before} velocities, |
---|
| 95 | the \textit{before} vorticity and divergence arrays must be included in the restart file to |
---|
| 96 | ensure perfect restartability. |
---|
| 97 | The vorticity and divergence at the \textit{now} time step are used for the computation of |
---|
| 98 | the nonlinear advection and of the vertical velocity respectively. |
---|
[2282] | 99 | |
---|
| 100 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 101 | % Sea Surface Height evolution |
---|
| 102 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 103 | \subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} |
---|
[9407] | 104 | \label{subsec:DYN_sshwzv} |
---|
[2282] | 105 | |
---|
[10368] | 106 | The sea surface height is given by: |
---|
[10419] | 107 | \begin{equation} |
---|
| 108 | \label{eq:dynspg_ssh} |
---|
| 109 | \begin{aligned} |
---|
| 110 | \frac{\partial \eta }{\partial t} |
---|
| 111 | &\equiv \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left\{ \delta_i \left[ {e_{2u}\,e_{3u}\;u} \right] |
---|
| 112 | +\delta_j \left[ {e_{1v}\,e_{3v}\;v} \right] \right\} } |
---|
| 113 | - \frac{\textit{emp}}{\rho_w } \\ |
---|
| 114 | &\equiv \sum\limits_k {\chi \ e_{3t}} - \frac{\textit{emp}}{\rho_w } |
---|
| 115 | \end{aligned} |
---|
[2282] | 116 | \end{equation} |
---|
| 117 | where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), |
---|
[10368] | 118 | expressed in Kg/m$^2$/s (which is equal to mm/s), |
---|
[10419] | 119 | and $\rho_w$=1,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation). |
---|
[10368] | 120 | If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then |
---|
| 121 | \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. |
---|
| 122 | The sea-surface height is evaluated using exactly the same time stepping scheme as |
---|
| 123 | the tracer equation \autoref{eq:tra_nxt}: |
---|
| 124 | a leapfrog scheme in combination with an Asselin time filter, |
---|
| 125 | $i.e.$ the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity). |
---|
| 126 | This is of paramount importance. |
---|
| 127 | Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to |
---|
| 128 | the sea surface height equation otherwise tracer content will not be conserved |
---|
| 129 | \citep{Griffies_al_MWR01, Leclair_Madec_OM09}. |
---|
[2282] | 130 | |
---|
[10368] | 131 | The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, |
---|
| 132 | taking into account the change of the thickness of the levels: |
---|
[10419] | 133 | \begin{equation} |
---|
| 134 | \label{eq:wzv} |
---|
| 135 | \left\{ |
---|
| 136 | \begin{aligned} |
---|
| 137 | &\left. w \right|_{k_b-1/2} \quad= 0 \qquad \text{where } k_b \text{ is the level just above the sea floor } \\ |
---|
| 138 | &\left. w \right|_{k+1/2} = \left. w \right|_{k-1/2} + \left. e_{3t} \right|_{k}\; \left. \chi \right|_k |
---|
| 139 | - \frac{1} {2 \rdt} \left( \left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k}\right) |
---|
| 140 | \end{aligned} |
---|
| 141 | \right. |
---|
[2282] | 142 | \end{equation} |
---|
| 143 | |
---|
| 144 | In the case of a non-linear free surface (\key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$, |
---|
[10368] | 145 | as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses, |
---|
| 146 | re-orientated downward. |
---|
[2285] | 147 | \gmcomment{not sure of this... to be modified with the change in emp setting} |
---|
[9407] | 148 | In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears. |
---|
[10368] | 149 | The upper boundary condition applies at a fixed level $z=0$. |
---|
| 150 | The top vertical velocity is thus equal to the divergence of the barotropic transport |
---|
| 151 | ($i.e.$ the first term in the right-hand-side of \autoref{eq:dynspg_ssh}). |
---|
[2282] | 152 | |
---|
[10368] | 153 | Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates, |
---|
| 154 | its physical meaning is not the same: |
---|
| 155 | in the second case, $w$ is the velocity normal to the $s$-surfaces. |
---|
| 156 | Note also that the $k$-axis is re-orientated downwards in the \textsc{fortran} code compared to |
---|
| 157 | the indexing used in the semi-discrete equations such as \autoref{eq:wzv} |
---|
[9407] | 158 | (see \autoref{subsec:DOM_Num_Index_vertical}). |
---|
[2282] | 159 | |
---|
| 160 | |
---|
| 161 | % ================================================================ |
---|
[707] | 162 | % Coriolis and Advection terms: vector invariant form |
---|
| 163 | % ================================================================ |
---|
[9393] | 164 | \section{Coriolis and advection: vector invariant form} |
---|
[9407] | 165 | \label{sec:DYN_adv_cor_vect} |
---|
[707] | 166 | %-----------------------------------------nam_dynadv---------------------------------------------------- |
---|
[10146] | 167 | |
---|
| 168 | \nlst{namdyn_adv} |
---|
[707] | 169 | %------------------------------------------------------------------------------------------------------------- |
---|
| 170 | |
---|
[10368] | 171 | The vector invariant form of the momentum equations is the one most often used in |
---|
| 172 | applications of the \NEMO ocean model. |
---|
| 173 | The flux form option (see next section) has been present since version $2$. |
---|
| 174 | Options are defined through the \ngn{namdyn\_adv} namelist variables Coriolis and |
---|
| 175 | momentum advection terms are evaluated using a leapfrog scheme, |
---|
| 176 | $i.e.$ the velocity appearing in these expressions is centred in time (\textit{now} velocity). |
---|
| 177 | At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following |
---|
| 178 | \autoref{chap:LBC}. |
---|
[707] | 179 | |
---|
| 180 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 181 | % Vorticity term |
---|
| 182 | % ------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 183 | \subsection{Vorticity term (\protect\mdl{dynvor})} |
---|
[9407] | 184 | \label{subsec:DYN_vor} |
---|
[707] | 185 | %------------------------------------------nam_dynvor---------------------------------------------------- |
---|
[10146] | 186 | |
---|
| 187 | \nlst{namdyn_vor} |
---|
[707] | 188 | %------------------------------------------------------------------------------------------------------------- |
---|
| 189 | |
---|
[4147] | 190 | Options are defined through the \ngn{namdyn\_vor} namelist variables. |
---|
[10368] | 191 | Four discretisations of the vorticity term (\np{ln\_dynvor\_xxx}\forcode{ = .true.}) are available: |
---|
| 192 | conserving potential enstrophy of horizontally non-divergent flow (ENS scheme); |
---|
| 193 | conserving horizontal kinetic energy (ENE scheme); |
---|
| 194 | conserving potential enstrophy for the relative vorticity term and |
---|
| 195 | horizontal kinetic energy for the planetary vorticity term (MIX scheme); |
---|
| 196 | or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy |
---|
| 197 | (EEN scheme) (see \autoref{subsec:C_vorEEN}). |
---|
| 198 | 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.}). |
---|
| 200 | The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. |
---|
[707] | 201 | |
---|
| 202 | %------------------------------------------------------------- |
---|
| 203 | % enstrophy conserving scheme |
---|
| 204 | %------------------------------------------------------------- |
---|
[9393] | 205 | \subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} |
---|
[9407] | 206 | \label{subsec:DYN_vor_ens} |
---|
[707] | 207 | |
---|
[10368] | 208 | In the enstrophy conserving case (ENS scheme), |
---|
| 209 | 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 ($i.e.$ $\chi$=$0$), |
---|
| 211 | but does not conserve the total kinetic energy. |
---|
| 212 | It is given by: |
---|
[10419] | 213 | \begin{equation} |
---|
| 214 | \label{eq:dynvor_ens} |
---|
| 215 | \left\{ |
---|
| 216 | \begin{aligned} |
---|
| 217 | {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} |
---|
| 218 | & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2} \\ |
---|
| 219 | {- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} |
---|
| 220 | & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j} |
---|
| 221 | \end{aligned} |
---|
| 222 | \right. |
---|
[707] | 223 | \end{equation} |
---|
| 224 | |
---|
| 225 | %------------------------------------------------------------- |
---|
| 226 | % energy conserving scheme |
---|
| 227 | %------------------------------------------------------------- |
---|
[9393] | 228 | \subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} |
---|
[9407] | 229 | \label{subsec:DYN_vor_ene} |
---|
[707] | 230 | |
---|
[10368] | 231 | The kinetic energy conserving scheme (ENE scheme) conserves the global kinetic energy but not the global enstrophy. |
---|
| 232 | It is given by: |
---|
[10419] | 233 | \begin{equation} |
---|
| 234 | \label{eq:dynvor_ene} |
---|
| 235 | \left\{ |
---|
| 236 | \begin{aligned} |
---|
| 237 | {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
| 238 | \; \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
| 239 | {- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
| 240 | \; \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } |
---|
| 241 | \end{aligned} |
---|
| 242 | \right. |
---|
[707] | 243 | \end{equation} |
---|
| 244 | |
---|
| 245 | %------------------------------------------------------------- |
---|
| 246 | % mix energy/enstrophy conserving scheme |
---|
| 247 | %------------------------------------------------------------- |
---|
[9393] | 248 | \subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } |
---|
[9407] | 249 | \label{subsec:DYN_vor_mix} |
---|
[707] | 250 | |
---|
[10368] | 251 | 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. |
---|
[10419] | 254 | \[ |
---|
| 255 | % \label{eq:dynvor_mix} |
---|
| 256 | \left\{ { |
---|
| 257 | \begin{aligned} |
---|
| 258 | {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} |
---|
| 259 | \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } |
---|
| 260 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
| 261 | \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
| 262 | {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j |
---|
| 263 | \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } |
---|
| 264 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
| 265 | \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill |
---|
| 266 | \end{aligned} |
---|
| 267 | } \right. |
---|
| 268 | \] |
---|
[707] | 269 | |
---|
| 270 | %------------------------------------------------------------- |
---|
| 271 | % energy and enstrophy conserving scheme |
---|
| 272 | %------------------------------------------------------------- |
---|
[9393] | 273 | \subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } |
---|
[9407] | 274 | \label{subsec:DYN_vor_een} |
---|
[707] | 275 | |
---|
[10368] | 276 | In both the ENS and ENE schemes, |
---|
| 277 | it is apparent that the combination of $i$ and $j$ averages of the velocity allows for |
---|
| 278 | the presence of grid point oscillation structures that will be invisible to the operator. |
---|
| 279 | These structures are \textit{computational modes} that will be at least partly damped by |
---|
| 280 | the momentum diffusion operator ($i.e.$ the subgrid-scale advection), but not by the resolved advection term. |
---|
| 281 | The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field. |
---|
| 282 | Such noise would result in more noise in the vertical velocity field, an undesirable feature. |
---|
| 283 | This is a well-known characteristic of $C$-grid discretization where |
---|
| 284 | $u$ and $v$ are located at different grid points, |
---|
| 285 | a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid. |
---|
[2282] | 286 | \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) |
---|
| 287 | Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} |
---|
| 288 | |
---|
[10368] | 289 | A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}. |
---|
| 290 | The idea is to get rid of the double averaging by considering triad combinations of vorticity. |
---|
[2282] | 291 | It is noteworthy that this solution is conceptually quite similar to the one proposed by |
---|
[9407] | 292 | \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). |
---|
[2282] | 293 | |
---|
| 294 | The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified |
---|
| 295 | 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: |
---|
[10419] | 297 | \[ |
---|
| 298 | % \label{eq:pot_vor} |
---|
| 299 | q = \frac{\zeta +f} {e_{3f} } |
---|
| 300 | \] |
---|
[10368] | 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: |
---|
[10419] | 303 | \begin{equation} |
---|
| 304 | \label{eq:een_e3f} |
---|
| 305 | e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} |
---|
[707] | 306 | \end{equation} |
---|
| 307 | |
---|
| 308 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
[10419] | 309 | \begin{figure}[!ht] |
---|
| 310 | \begin{center} |
---|
| 311 | \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} |
---|
| 312 | \caption{ |
---|
| 313 | \protect\label{fig:DYN_een_triad} |
---|
| 314 | Triads used in the energy and enstrophy conserving scheme (een) for |
---|
| 315 | $u$-component (upper panel) and $v$-component (lower panel). |
---|
| 316 | } |
---|
| 317 | \end{center} |
---|
| 318 | \end{figure} |
---|
| 319 | % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
[707] | 320 | |
---|
[9407] | 321 | A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. |
---|
[10368] | 322 | 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.}). |
---|
| 324 | The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and |
---|
| 325 | extends by continuity the value of $e_{3f}$ into the land areas. |
---|
| 326 | This case introduces a sub-grid-scale topography at f-points |
---|
| 327 | (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) |
---|
| 328 | that tends to reinforce the topostrophy of the flow |
---|
[6289] | 329 | ($i.e.$ the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}. |
---|
[707] | 330 | |
---|
[10368] | 331 | Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as |
---|
| 332 | the following triad combinations of the neighbouring potential vorticities defined at f-points |
---|
[9407] | 333 | (\autoref{fig:DYN_een_triad}): |
---|
[10419] | 334 | \begin{equation} |
---|
| 335 | \label{eq:Q_triads} |
---|
| 336 | _i^j \mathbb{Q}^{i_p}_{j_p} |
---|
| 337 | = \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) |
---|
[2282] | 338 | \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: |
---|
[10419] | 342 | \begin{equation} |
---|
| 343 | \label{eq:dynvor_een} |
---|
| 344 | \left\{ { |
---|
| 345 | \begin{aligned} |
---|
| 346 | +q\,e_3 \, v &\equiv +\frac{1}{e_{1u} } \sum_{\substack{i_p,\,k_p}} |
---|
| 347 | {^{i+1/2-i_p}_j} \mathbb{Q}^{i_p}_{j_p} \left( e_{1v}\,e_{3v} \;v \right)^{i+1/2-i_p}_{j+j_p} \\ |
---|
| 348 | - q\,e_3 \, u &\equiv -\frac{1}{e_{2v} } \sum_{\substack{i_p,\,k_p}} |
---|
| 349 | {^i_{j+1/2-j_p}} \mathbb{Q}^{i_p}_{j_p} \left( e_{2u}\,e_{3u} \;u \right)^{i+i_p}_{j+1/2-j_p} \\ |
---|
| 350 | \end{aligned} |
---|
| 351 | } \right. |
---|
[707] | 352 | \end{equation} |
---|
| 353 | |
---|
[10368] | 354 | This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. |
---|
| 355 | It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow |
---|
| 356 | ($i.e.$ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}). |
---|
| 357 | 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}. |
---|
[2282] | 359 | Furthermore, used in combination with a partial steps representation of bottom topography, |
---|
[10368] | 360 | it improves the interaction between current and topography, |
---|
| 361 | leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}. |
---|
[2282] | 362 | |
---|
[707] | 363 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 364 | % Kinetic Energy Gradient term |
---|
| 365 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 366 | \subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} |
---|
[9407] | 367 | \label{subsec:DYN_keg} |
---|
[707] | 368 | |
---|
[10368] | 369 | As demonstrated in \autoref{apdx:C}, |
---|
| 370 | there is a single discrete formulation of the kinetic energy gradient term that, |
---|
| 371 | together with the formulation chosen for the vertical advection (see below), |
---|
| 372 | conserves the total kinetic energy: |
---|
[10419] | 373 | \[ |
---|
| 374 | % \label{eq:dynkeg} |
---|
| 375 | \left\{ |
---|
| 376 | \begin{aligned} |
---|
| 377 | -\frac{1}{2 \; e_{1u} } & \ \delta_{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] \\ |
---|
| 378 | -\frac{1}{2 \; e_{2v} } & \ \delta_{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] |
---|
| 379 | \end{aligned} |
---|
| 380 | \right. |
---|
| 381 | \] |
---|
[707] | 382 | |
---|
| 383 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 384 | % Vertical advection term |
---|
| 385 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 386 | \subsection{Vertical advection term (\protect\mdl{dynzad}) } |
---|
[9407] | 387 | \label{subsec:DYN_zad} |
---|
[707] | 388 | |
---|
[10368] | 389 | The discrete formulation of the vertical advection, t |
---|
| 390 | ogether with the formulation chosen for the gradient of kinetic energy (KE) term, |
---|
| 391 | conserves the total kinetic energy. |
---|
| 392 | 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}). |
---|
[10419] | 394 | \[ |
---|
| 395 | % \label{eq:dynzad} |
---|
| 396 | \left\{ |
---|
| 397 | \begin{aligned} |
---|
| 398 | -\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2} \;\delta_{k+1/2} \left[ u \right]\ }^{\,k} \\ |
---|
| 399 | -\frac{1} {e_{1v}\,e_{2v}\,e_{3v}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,j+1/2} \;\delta_{k+1/2} \left[ u \right]\ }^{\,k} |
---|
| 400 | \end{aligned} |
---|
| 401 | \right. |
---|
| 402 | \] |
---|
[10368] | 403 | When \np{ln\_dynzad\_zts}\forcode{ = .true.}, |
---|
| 404 | a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. |
---|
[6289] | 405 | This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}. |
---|
[10368] | 406 | Note that in this case, |
---|
| 407 | a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, |
---|
[9407] | 408 | an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). |
---|
[707] | 409 | |
---|
[6289] | 410 | |
---|
[707] | 411 | % ================================================================ |
---|
| 412 | % Coriolis and Advection : flux form |
---|
| 413 | % ================================================================ |
---|
[9393] | 414 | \section{Coriolis and advection: flux form} |
---|
[9407] | 415 | \label{sec:DYN_adv_cor_flux} |
---|
[707] | 416 | %------------------------------------------nam_dynadv---------------------------------------------------- |
---|
[10146] | 417 | |
---|
| 418 | \nlst{namdyn_adv} |
---|
[707] | 419 | %------------------------------------------------------------------------------------------------------------- |
---|
| 420 | |
---|
[4147] | 421 | Options are defined through the \ngn{namdyn\_adv} namelist variables. |
---|
[10368] | 422 | In the flux form (as in the vector invariant form), |
---|
| 423 | the Coriolis and momentum advection terms are evaluated using a leapfrog scheme, |
---|
| 424 | $i.e.$ the velocity appearing in their expressions is centred in time (\textit{now} velocity). |
---|
| 425 | At the lateral boundaries either free slip, |
---|
| 426 | no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. |
---|
[707] | 427 | |
---|
| 428 | |
---|
| 429 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 430 | % Coriolis plus curvature metric terms |
---|
| 431 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 432 | \subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } |
---|
[9407] | 433 | \label{subsec:DYN_cor_flux} |
---|
[707] | 434 | |
---|
[10368] | 435 | 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 | This altered Coriolis parameter is thus discretised at $f$-points. |
---|
| 437 | It is given by: |
---|
[10419] | 438 | \begin{multline*} |
---|
| 439 | % \label{eq:dyncor_metric} |
---|
| 440 | f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right) \\ |
---|
| 441 | \equiv f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right] |
---|
| 442 | - \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right] } \ \right) |
---|
| 443 | \end{multline*} |
---|
[707] | 444 | |
---|
[10368] | 445 | Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to |
---|
| 446 | 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, $i.e.$ the velocity is centred in time (\textit{now} velocity). |
---|
[707] | 449 | |
---|
| 450 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 451 | % Flux form Advection term |
---|
| 452 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 453 | \subsection{Flux form advection term (\protect\mdl{dynadv}) } |
---|
[9407] | 454 | \label{subsec:DYN_adv_flux} |
---|
[707] | 455 | |
---|
[10368] | 456 | The discrete expression of the advection term is given by: |
---|
[10419] | 457 | \[ |
---|
| 458 | % \label{eq:dynadv} |
---|
| 459 | \left\{ |
---|
| 460 | \begin{aligned} |
---|
| 461 | \frac{1}{e_{1u}\,e_{2u}\,e_{3u}} |
---|
| 462 | \left( \delta_{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i } \ u_t \right] |
---|
| 463 | + \delta_{j } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2} \ u_f \right] \right. \ \; \\ |
---|
| 464 | \left. + \delta_{k } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2} \ u_{uw} \right] \right) \\ |
---|
| 465 | \\ |
---|
| 466 | \frac{1}{e_{1v}\,e_{2v}\,e_{3v}} |
---|
| 467 | \left( \delta_{i } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f \right] |
---|
| 468 | + \delta_{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i } \ v_t \right] \right. \ \, \, \\ |
---|
| 469 | \left. + \delta_{k } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw} \right] \right) \\ |
---|
| 470 | \end{aligned} |
---|
| 471 | \right. |
---|
| 472 | \] |
---|
[707] | 473 | |
---|
[10368] | 474 | Two advection schemes are available: |
---|
| 475 | a $2^{nd}$ order centered finite difference scheme, CEN2, |
---|
| 476 | 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}. |
---|
| 479 | 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, $i.e.$ at the $T$-, $f$-, |
---|
| 481 | and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. |
---|
[707] | 482 | |
---|
| 483 | %------------------------------------------------------------- |
---|
| 484 | % 2nd order centred scheme |
---|
| 485 | %------------------------------------------------------------- |
---|
[9393] | 486 | \subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} |
---|
[9407] | 487 | \label{subsec:DYN_adv_cen2} |
---|
[707] | 488 | |
---|
[10368] | 489 | In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points: |
---|
[10419] | 490 | \begin{equation} |
---|
| 491 | \label{eq:dynadv_cen2} |
---|
| 492 | \left\{ |
---|
| 493 | \begin{aligned} |
---|
| 494 | u_T^{cen2} &=\overline u^{i } \quad & u_F^{cen2} &=\overline u^{j+1/2} \quad & u_{uw}^{cen2} &=\overline u^{k+1/2} \\ |
---|
| 495 | v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j \quad & v_{vw}^{cen2} &=\overline v ^{k+1/2} \\ |
---|
| 496 | \end{aligned} |
---|
| 497 | \right. |
---|
[707] | 498 | \end{equation} |
---|
| 499 | |
---|
[10368] | 500 | The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive ($i.e.$ it may create false extrema). |
---|
| 501 | It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to |
---|
| 502 | produce a sensible solution. |
---|
| 503 | The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, |
---|
| 504 | so $u$ and $v$ are the \emph{now} velocities. |
---|
[707] | 505 | |
---|
| 506 | %------------------------------------------------------------- |
---|
| 507 | % UBS scheme |
---|
| 508 | %------------------------------------------------------------- |
---|
[9393] | 509 | \subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} |
---|
[9407] | 510 | \label{subsec:DYN_adv_ubs} |
---|
[707] | 511 | |
---|
[10368] | 512 | The UBS advection scheme is an upstream biased third order scheme based on |
---|
| 513 | an upstream-biased parabolic interpolation. |
---|
| 514 | For example, the evaluation of $u_T^{ubs} $ is done as follows: |
---|
[10419] | 515 | \begin{equation} |
---|
| 516 | \label{eq:dynadv_ubs} |
---|
| 517 | u_T^{ubs} =\overline u ^i-\;\frac{1}{6} |
---|
| 518 | \begin{cases} |
---|
| 519 | u"_{i-1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i \geqslant 0$ } \\ |
---|
| 520 | u"_{i+1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i < 0$ } |
---|
| 521 | \end{cases} |
---|
[707] | 522 | \end{equation} |
---|
[10419] | 523 | where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. |
---|
[10368] | 524 | This results in a dissipatively dominant ($i.e.$ hyper-diffusive) truncation error |
---|
| 525 | \citep{Shchepetkin_McWilliams_OM05}. |
---|
| 526 | The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}. |
---|
| 527 | It is a relatively good compromise between accuracy and smoothness. |
---|
| 528 | It is not a \emph{positive} scheme, meaning that false extrema are permitted. |
---|
| 529 | 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 | ($i.e.$ \np{ln\_dynldf\_lap}\forcode{ = }\np{ln\_dynldf\_bilap}\forcode{ = .false.}), |
---|
| 532 | and it is recommended to do so. |
---|
[707] | 533 | |
---|
[10368] | 534 | The UBS scheme is not used in all directions. |
---|
| 535 | In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and |
---|
| 536 | $u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used. |
---|
| 537 | UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm pursue the |
---|
[817] | 538 | sentence:Since vertical mixing of momentum is a source term of the TKE equation... } |
---|
[707] | 539 | |
---|
[10368] | 540 | For stability reasons, the first term in (\autoref{eq:dynadv_ubs}), |
---|
| 541 | which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time), |
---|
| 542 | while the second term, which is the diffusion part of the scheme, |
---|
| 543 | 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. |
---|
[707] | 545 | |
---|
[10368] | 546 | Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by |
---|
| 547 | one coefficient. |
---|
| 548 | Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. |
---|
| 549 | This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. |
---|
| 550 | Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. |
---|
[707] | 551 | |
---|
[10368] | 552 | Note also that in the current version of \mdl{dynadv\_ubs}, |
---|
| 553 | there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS. |
---|
| 554 | This is an error and should be suppressed soon. |
---|
[817] | 555 | %%% |
---|
| 556 | \gmcomment{action : this have to be done} |
---|
| 557 | %%% |
---|
[707] | 558 | |
---|
| 559 | % ================================================================ |
---|
| 560 | % Hydrostatic pressure gradient term |
---|
| 561 | % ================================================================ |
---|
[9393] | 562 | \section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} |
---|
[9407] | 563 | \label{sec:DYN_hpg} |
---|
[707] | 564 | %------------------------------------------nam_dynhpg--------------------------------------------------- |
---|
[10146] | 565 | |
---|
| 566 | \nlst{namdyn_hpg} |
---|
[707] | 567 | %------------------------------------------------------------------------------------------------------------- |
---|
| 568 | |
---|
[4147] | 569 | Options are defined through the \ngn{namdyn\_hpg} namelist variables. |
---|
[10368] | 570 | The key distinction between the different algorithms used for |
---|
| 571 | the hydrostatic pressure gradient is the vertical coordinate used, |
---|
| 572 | since HPG is a \emph{horizontal} pressure gradient, $i.e.$ computed along geopotential surfaces. |
---|
| 573 | As a result, any tilt of the surface of the computational levels will require a specific treatment to |
---|
[817] | 574 | compute the hydrostatic pressure gradient. |
---|
[707] | 575 | |
---|
[10368] | 576 | The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, |
---|
| 577 | $i.e.$ the density appearing in its expression is centred in time (\emph{now} $\rho$), |
---|
| 578 | or a semi-implcit scheme. |
---|
| 579 | At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied. |
---|
[707] | 580 | |
---|
| 581 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 582 | % z-coordinate with full step |
---|
| 583 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 584 | \subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})} |
---|
[9407] | 585 | \label{subsec:DYN_hpg_zco} |
---|
[707] | 586 | |
---|
[10368] | 587 | The hydrostatic pressure can be obtained by integrating the hydrostatic equation vertically from the surface. |
---|
| 588 | However, the pressure is large at great depth while its horizontal gradient is several orders of magnitude smaller. |
---|
| 589 | This may lead to large truncation errors in the pressure gradient terms. |
---|
| 590 | Thus, the two horizontal components of the hydrostatic pressure gradient are computed directly as follows: |
---|
[707] | 591 | |
---|
| 592 | for $k=km$ (surface layer, $jk=1$ in the code) |
---|
[10419] | 593 | \begin{equation} |
---|
| 594 | \label{eq:dynhpg_zco_surf} |
---|
| 595 | \left\{ |
---|
| 596 | \begin{aligned} |
---|
| 597 | \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k=km} |
---|
| 598 | &= \frac{1}{2} g \ \left. \delta_{i+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
| 599 | \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k=km} |
---|
| 600 | &= \frac{1}{2} g \ \left. \delta_{j+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
| 601 | \end{aligned} |
---|
| 602 | \right. |
---|
[707] | 603 | \end{equation} |
---|
| 604 | |
---|
| 605 | for $1<k<km$ (interior layer) |
---|
[10419] | 606 | \begin{equation} |
---|
| 607 | \label{eq:dynhpg_zco} |
---|
| 608 | \left\{ |
---|
| 609 | \begin{aligned} |
---|
| 610 | \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k} |
---|
| 611 | &= \left. \delta_{i+1/2} \left[ p^h \right] \right|_{k-1} |
---|
| 612 | + \frac{1}{2}\;g\; \left. \delta_{i+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
| 613 | \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k} |
---|
| 614 | &= \left. \delta_{j+1/2} \left[ p^h \right] \right|_{k-1} |
---|
| 615 | + \frac{1}{2}\;g\; \left. \delta_{j+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
| 616 | \end{aligned} |
---|
| 617 | \right. |
---|
[707] | 618 | \end{equation} |
---|
| 619 | |
---|
[10368] | 620 | Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as |
---|
| 621 | 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}$. |
---|
[707] | 625 | |
---|
| 626 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 627 | % z-coordinate with partial step |
---|
| 628 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 629 | \subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})} |
---|
[9407] | 630 | \label{subsec:DYN_hpg_zps} |
---|
[707] | 631 | |
---|
[10368] | 632 | With partial bottom cells, tracers in horizontally adjacent cells generally live at different depths. |
---|
| 633 | Before taking horizontal gradients between these tracer points, |
---|
| 634 | a linear interpolation is used to approximate the deeper tracer as if |
---|
| 635 | it actually lived at the depth of the shallower tracer point. |
---|
[707] | 636 | |
---|
[10368] | 637 | Apart from this modification, |
---|
| 638 | the horizontal hydrostatic pressure gradient evaluated in the $z$-coordinate with partial step is exactly as in |
---|
| 639 | the pure $z$-coordinate case. |
---|
| 640 | As explained in detail in section \autoref{sec:TRA_zpshde}, |
---|
| 641 | the nonlinearity of pressure effects in the equation of state is such that |
---|
| 642 | it is better to interpolate temperature and salinity vertically before computing the density. |
---|
| 643 | Horizontal gradients of temperature and salinity are needed for the TRA modules, |
---|
| 644 | which is the reason why the horizontal gradients of density at the deepest model level are computed in |
---|
| 645 | module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}. |
---|
[707] | 646 | |
---|
| 647 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 648 | % s- and s-z-coordinates |
---|
| 649 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 650 | \subsection{$S$- and $Z$-$S$-coordinates} |
---|
[9407] | 651 | \label{subsec:DYN_hpg_sco} |
---|
[707] | 652 | |
---|
[10368] | 653 | Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers |
---|
| 654 | ($e.g.$, \citet{Song1998, Shchepetkin_McWilliams_OM05}). |
---|
| 655 | A number of different pressure gradient options are coded but the ROMS-like, |
---|
| 656 | density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation. |
---|
[707] | 657 | |
---|
[9393] | 658 | $\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.}) |
---|
[10419] | 659 | \begin{equation} |
---|
| 660 | \label{eq:dynhpg_sco} |
---|
| 661 | \left\{ |
---|
| 662 | \begin{aligned} |
---|
| 663 | - \frac{1} {\rho_o \, e_{1u}} \; \delta_{i+1/2} \left[ p^h \right] |
---|
| 664 | + \frac{g\; \overline {\rho}^{i+1/2}} {\rho_o \, e_{1u}} \; \delta_{i+1/2} \left[ z_t \right] \\ |
---|
| 665 | - \frac{1} {\rho_o \, e_{2v}} \; \delta_{j+1/2} \left[ p^h \right] |
---|
| 666 | + \frac{g\; \overline {\rho}^{j+1/2}} {\rho_o \, e_{2v}} \; \delta_{j+1/2} \left[ z_t \right] \\ |
---|
| 667 | \end{aligned} |
---|
| 668 | \right. |
---|
[707] | 669 | \end{equation} |
---|
| 670 | |
---|
[10368] | 671 | 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 |
---|
[5120] | 674 | ($e_{3w}$). |
---|
| 675 | |
---|
[9393] | 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.}). |
---|
[707] | 678 | |
---|
[9393] | 679 | $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) |
---|
[707] | 680 | |
---|
[2282] | 681 | $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} |
---|
[9393] | 682 | (\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development) |
---|
[707] | 683 | |
---|
[10368] | 684 | Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated |
---|
| 685 | (\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 as |
---|
| 688 | an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active. |
---|
| 689 | The pressure Jacobian scheme uses a constrained cubic spline to |
---|
| 690 | reconstruct the density profile across the water column. |
---|
| 691 | This method maintains the monotonicity between the density nodes. |
---|
| 692 | The pressure can be calculated by analytical integration of the density profile and |
---|
| 693 | a pressure Jacobian method is used to solve the horizontal pressure gradient. |
---|
| 694 | This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme. |
---|
[707] | 695 | |
---|
[6320] | 696 | \subsection{Ice shelf cavity} |
---|
[9407] | 697 | \label{subsec:DYN_hpg_isf} |
---|
[6320] | 698 | Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and |
---|
[10368] | 699 | the pressure gradient due to the ocean load. |
---|
| 700 | If cavity opened (\np{ln\_isfcav}\forcode{ = .true.}) these 2 terms can be calculated by |
---|
| 701 | setting \np{ln\_dynhpg\_isf}\forcode{ = .true.}. |
---|
| 702 | No other scheme are working with the ice shelf.\\ |
---|
[6320] | 703 | |
---|
| 704 | $\bullet$ The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium. |
---|
[10368] | 705 | The top pressure is computed integrating from surface to the base of the ice shelf a reference density profile |
---|
| 706 | (prescribed as density of a water at 34.4 PSU and -1.9\degC) and |
---|
| 707 | corresponds to the water replaced by the ice shelf. |
---|
| 708 | This top pressure is constant over time. |
---|
| 709 | A detailed description of this method is described in \citet{Losch2008}.\\ |
---|
[6320] | 710 | |
---|
[10368] | 711 | $\bullet$ The ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in |
---|
| 712 | \autoref{subsec:DYN_hpg_sco}. |
---|
[6320] | 713 | |
---|
[707] | 714 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 715 | % Time-scheme |
---|
| 716 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 717 | \subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})} |
---|
[9407] | 718 | \label{subsec:DYN_hpg_imp} |
---|
[707] | 719 | |
---|
[10368] | 720 | The default time differencing scheme used for the horizontal pressure gradient is a leapfrog scheme and |
---|
| 721 | therefore the density used in all discrete expressions given above is the \textit{now} density, |
---|
| 722 | computed from the \textit{now} temperature and salinity. |
---|
| 723 | In some specific cases |
---|
| 724 | (usually high resolution simulations over an ocean domain which includes weakly stratified regions) |
---|
| 725 | the physical phenomenon that controls the time-step is internal gravity waves (IGWs). |
---|
| 726 | A semi-implicit scheme for doubling the stability limit associated with IGWs can be used |
---|
| 727 | \citep{Brown_Campana_MWR78, Maltrud1998}. |
---|
| 728 | It involves the evaluation of the hydrostatic pressure gradient as |
---|
| 729 | an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ |
---|
| 730 | ($i.e.$ \textit{before}, \textit{now} and \textit{after} time-steps), |
---|
| 731 | rather than at the central time level $t$ only, as in the standard leapfrog scheme. |
---|
[707] | 732 | |
---|
[9393] | 733 | $\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}): |
---|
[707] | 734 | |
---|
[10419] | 735 | \begin{equation} |
---|
| 736 | \label{eq:dynhpg_lf} |
---|
| 737 | \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; |
---|
| 738 | -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right] |
---|
[707] | 739 | \end{equation} |
---|
| 740 | |
---|
[9393] | 741 | $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}): |
---|
[10419] | 742 | \begin{equation} |
---|
| 743 | \label{eq:dynhpg_imp} |
---|
| 744 | \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; |
---|
| 745 | -\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] |
---|
[707] | 746 | \end{equation} |
---|
| 747 | |
---|
[10368] | 748 | The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without |
---|
| 749 | significant additional computation since the density can be updated to time level $t+\rdt$ before |
---|
| 750 | computing the horizontal hydrostatic pressure gradient. |
---|
| 751 | It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using |
---|
| 752 | \autoref{eq:dynhpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:dynhpg_lf}. |
---|
| 753 | Note that \autoref{eq:dynhpg_imp} is equivalent to applying a time filter to the pressure gradient to |
---|
| 754 | eliminate high frequency IGWs. |
---|
| 755 | Obviously, when using \autoref{eq:dynhpg_imp}, |
---|
| 756 | the doubling of the time-step is achievable only if no other factors control the time-step, |
---|
| 757 | such as the stability limits associated with advection or diffusion. |
---|
[707] | 758 | |
---|
[10368] | 759 | In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}\forcode{ = .true.}. |
---|
| 760 | In this case, we choose to apply the time filter to temperature and salinity used in the equation of state, |
---|
| 761 | instead of applying it to the hydrostatic pressure or to the density, |
---|
| 762 | so that no additional storage array has to be defined. |
---|
| 763 | The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows: |
---|
[10419] | 764 | \[ |
---|
| 765 | % \label{eq:rho_flt} |
---|
| 766 | \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) |
---|
| 767 | \quad \text{with} \quad |
---|
| 768 | \widetilde{X} = 1 / 4 \left( X^{t+\rdt} +2 \,X^t + X^{t-\rdt} \right) |
---|
| 769 | \] |
---|
[707] | 770 | |
---|
[10368] | 771 | Note that in the semi-implicit case, it is necessary to save the filtered density, |
---|
| 772 | an extra three-dimensional field, in the restart file to restart the model with exact reproducibility. |
---|
| 773 | This option is controlled by \np{nn\_dynhpg\_rst}, a namelist parameter. |
---|
[707] | 774 | |
---|
| 775 | % ================================================================ |
---|
| 776 | % Surface Pressure Gradient |
---|
| 777 | % ================================================================ |
---|
[9393] | 778 | \section{Surface pressure gradient (\protect\mdl{dynspg})} |
---|
[9407] | 779 | \label{sec:DYN_spg} |
---|
[707] | 780 | %-----------------------------------------nam_dynspg---------------------------------------------------- |
---|
[10146] | 781 | |
---|
| 782 | \nlst{namdyn_spg} |
---|
[707] | 783 | %------------------------------------------------------------------------------------------------------------ |
---|
| 784 | |
---|
[4147] | 785 | Options are defined through the \ngn{namdyn\_spg} namelist variables. |
---|
[10368] | 786 | The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}). |
---|
| 787 | The main distinction is between the fixed volume case (linear free surface) and |
---|
| 788 | the variable volume case (nonlinear free surface, \key{vvl} is defined). |
---|
| 789 | In the linear free surface case (\autoref{subsec:PE_free_surface}) |
---|
| 790 | the vertical scale factors $e_{3}$ are fixed in time, |
---|
| 791 | while they are time-dependent in the nonlinear case (\autoref{subsec:PE_free_surface}). |
---|
[6289] | 792 | With both linear and nonlinear free surface, external gravity waves are allowed in the equations, |
---|
[10368] | 793 | which imposes a very small time step when an explicit time stepping is used. |
---|
[6289] | 794 | Two methods are proposed to allow a longer time step for the three-dimensional equations: |
---|
[9407] | 795 | the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}), |
---|
[10368] | 796 | and the split-explicit free surface described below. |
---|
[6289] | 797 | The extra term introduced in the filtered method is calculated implicitly, |
---|
| 798 | so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
[2282] | 799 | |
---|
| 800 | |
---|
[10368] | 801 | The form of the surface pressure gradient term depends on how the user wants to |
---|
| 802 | handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}). |
---|
[2282] | 803 | Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): |
---|
[10368] | 804 | an explicit formulation which requires a small time step; |
---|
| 805 | a filtered free surface formulation which allows a larger time step by |
---|
| 806 | adding a filtering term into the momentum equation; |
---|
[2282] | 807 | and a split-explicit free surface formulation, described below, which also allows a larger time step. |
---|
| 808 | |
---|
[10368] | 809 | The extra term introduced in the filtered method is calculated implicitly, so that a solver is used to compute it. |
---|
| 810 | As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
[2282] | 811 | |
---|
| 812 | |
---|
[707] | 813 | %-------------------------------------------------------------------------------------------------------------- |
---|
[2282] | 814 | % Explicit free surface formulation |
---|
[707] | 815 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9363] | 816 | \subsection{Explicit free surface (\protect\key{dynspg\_exp})} |
---|
[9407] | 817 | \label{subsec:DYN_spg_exp} |
---|
[707] | 818 | |
---|
[10368] | 819 | In the explicit free surface formulation (\key{dynspg\_exp} defined), |
---|
| 820 | the model time step is chosen to be small enough to resolve the external gravity waves |
---|
| 821 | (typically a few tens of seconds). |
---|
[2282] | 822 | The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time), |
---|
| 823 | is thus simply given by : |
---|
[10419] | 824 | \begin{equation} |
---|
| 825 | \label{eq:dynspg_exp} |
---|
| 826 | \left\{ |
---|
| 827 | \begin{aligned} |
---|
| 828 | - \frac{1}{e_{1u}\,\rho_o} \; \delta_{i+1/2} \left[ \,\rho \,\eta\, \right] \\ |
---|
| 829 | - \frac{1}{e_{2v}\,\rho_o} \; \delta_{j+1/2} \left[ \,\rho \,\eta\, \right] |
---|
| 830 | \end{aligned} |
---|
| 831 | \right. |
---|
[707] | 832 | \end{equation} |
---|
| 833 | |
---|
[10368] | 834 | Note that in the non-linear free surface case ($i.e.$ \key{vvl} defined), |
---|
| 835 | the surface pressure gradient is already included in the momentum tendency through |
---|
| 836 | the level thickness variation allowed in the computation of the hydrostatic pressure gradient. |
---|
| 837 | Thus, nothing is done in the \mdl{dynspg\_exp} module. |
---|
[707] | 838 | |
---|
[2282] | 839 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 840 | % Split-explict free surface formulation |
---|
| 841 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 842 | \subsection{Split-explicit free surface (\protect\key{dynspg\_ts})} |
---|
[9407] | 843 | \label{subsec:DYN_spg_ts} |
---|
[4560] | 844 | %------------------------------------------namsplit----------------------------------------------------------- |
---|
[10146] | 845 | % |
---|
| 846 | %\nlst{namsplit} |
---|
[4560] | 847 | %------------------------------------------------------------------------------------------------------------- |
---|
[707] | 848 | |
---|
[2282] | 849 | The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined), |
---|
[10368] | 850 | also called the time-splitting formulation, follows the one proposed by \citet{Shchepetkin_McWilliams_OM05}. |
---|
| 851 | The general idea is to solve the free surface equation and the associated barotropic velocity equations with |
---|
| 852 | a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables |
---|
| 853 | (\autoref{fig:DYN_dynspg_ts}). |
---|
| 854 | The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through |
---|
| 855 | the \np{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$. |
---|
| 856 | This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}\forcode{ = .true.}) considering that |
---|
| 857 | the stability of the barotropic system is essentially controled by external waves propagation. |
---|
[6289] | 858 | Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry. |
---|
[9393] | 859 | Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}. |
---|
[817] | 860 | |
---|
[4560] | 861 | %%% |
---|
| 862 | The barotropic mode solves the following equations: |
---|
[10419] | 863 | % \begin{subequations} |
---|
| 864 | % \label{eq:BT} |
---|
| 865 | \begin{equation} |
---|
| 866 | \label{eq:BT_dyn} |
---|
| 867 | \frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}= |
---|
| 868 | -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h} |
---|
| 869 | -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}} |
---|
| 870 | \end{equation} |
---|
| 871 | \[ |
---|
| 872 | % \label{eq:BT_ssh} |
---|
| 873 | \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E |
---|
| 874 | \] |
---|
| 875 | % \end{subequations} |
---|
[10368] | 876 | where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes, |
---|
| 877 | surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency. |
---|
| 878 | The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress |
---|
| 879 | (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration. |
---|
| 880 | Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm |
---|
| 881 | detailed in \citet{Shchepetkin_McWilliams_OM05}. |
---|
| 882 | AB3-AM4 coefficients used in \NEMO follow the second-order accurate, |
---|
| 883 | "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08} |
---|
| 884 | (see their figure 12, lower left). |
---|
[4560] | 885 | |
---|
[707] | 886 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
[10419] | 887 | \begin{figure}[!t] |
---|
| 888 | \begin{center} |
---|
| 889 | \includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts} |
---|
| 890 | \caption{ |
---|
| 891 | \protect\label{fig:DYN_dynspg_ts} |
---|
| 892 | Schematic of the split-explicit time stepping scheme for the external and internal modes. |
---|
| 893 | Time increases to the right. In this particular exemple, |
---|
| 894 | a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$. |
---|
| 895 | Internal mode time steps (which are also the model time steps) are denoted by $t-\rdt$, $t$ and $t+\rdt$. |
---|
| 896 | Variables with $k$ superscript refer to instantaneous barotropic variables, |
---|
| 897 | $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary (red vertical bars) and |
---|
| 898 | secondary weights (blue vertical bars). |
---|
| 899 | The former are used to obtain time filtered quantities at $t+\rdt$ while |
---|
| 900 | the latter are used to obtain time averaged transports to advect tracers. |
---|
| 901 | a) Forward time integration: \protect\np{ln\_bt\_fw}\forcode{ = .true.}, |
---|
| 902 | \protect\np{ln\_bt\_av}\forcode{ = .true.}. |
---|
| 903 | b) Centred time integration: \protect\np{ln\_bt\_fw}\forcode{ = .false.}, |
---|
| 904 | \protect\np{ln\_bt\_av}\forcode{ = .true.}. |
---|
| 905 | c) Forward time integration with no time filtering (POM-like scheme): |
---|
| 906 | \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .false.}. |
---|
| 907 | } |
---|
| 908 | \end{center} |
---|
| 909 | \end{figure} |
---|
[707] | 910 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
| 911 | |
---|
[10368] | 912 | In the default case (\np{ln\_bt\_fw}\forcode{ = .true.}), |
---|
| 913 | the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps |
---|
| 914 | (\autoref{fig:DYN_dynspg_ts}a). |
---|
| 915 | To avoid aliasing of fast barotropic motions into three dimensional equations, |
---|
| 916 | time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{ = .true.}). |
---|
| 917 | In that case, the integration is extended slightly beyond \textit{after} time step to |
---|
| 918 | provide time filtered quantities. |
---|
| 919 | These are used for the subsequent initialization of the barotropic mode in the following baroclinic step. |
---|
[4560] | 920 | Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme, |
---|
[10368] | 921 | asselin filtering is not applied to barotropic quantities.\\ |
---|
| 922 | Alternatively, one can choose to integrate barotropic equations starting from \textit{before} time step |
---|
| 923 | (\np{ln\_bt\_fw}\forcode{ = .false.}). |
---|
| 924 | Although more computationaly expensive ( \np{nn\_baro} additional iterations are indeed necessary), |
---|
| 925 | the baroclinic to barotropic forcing term given at \textit{now} time step become centred in |
---|
| 926 | the middle of the integration window. |
---|
| 927 | It can easily be shown that this property removes part of splitting errors between modes, |
---|
| 928 | which increases the overall numerical robustness. |
---|
[4560] | 929 | %references to Patrick Marsaleix' work here. Also work done by SHOM group. |
---|
[707] | 930 | |
---|
[4560] | 931 | %%% |
---|
| 932 | |
---|
[10368] | 933 | As far as tracer conservation is concerned, |
---|
| 934 | barotropic velocities used to advect tracers must also be updated at \textit{now} time step. |
---|
| 935 | This implies to change the traditional order of computations in \NEMO: |
---|
| 936 | most of momentum trends (including the barotropic mode calculation) updated first, tracers' after. |
---|
| 937 | This \textit{de facto} makes semi-implicit hydrostatic pressure gradient |
---|
| 938 | (see section \autoref{subsec:DYN_hpg_imp}) |
---|
| 939 | and time splitting not compatible. |
---|
| 940 | Advective barotropic velocities are obtained by using a secondary set of filtering weights, |
---|
| 941 | uniquely defined from the filter coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}). |
---|
| 942 | Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to |
---|
| 943 | obtain exact conservation. |
---|
[4560] | 944 | |
---|
| 945 | %%% |
---|
| 946 | |
---|
[10368] | 947 | One can eventually choose to feedback instantaneous values by not using any time filter |
---|
| 948 | (\np{ln\_bt\_av}\forcode{ = .false.}). |
---|
| 949 | In that case, external mode equations are continuous in time, |
---|
| 950 | $i.e.$ they are not re-initialized when starting a new sub-stepping sequence. |
---|
| 951 | This is the method used so far in the POM model, the stability being maintained by |
---|
| 952 | refreshing at (almost) each barotropic time step advection and horizontal diffusion terms. |
---|
| 953 | Since the latter terms have not been added in \NEMO for computational efficiency, |
---|
| 954 | removing time filtering is not recommended except for debugging purposes. |
---|
| 955 | This may be used for instance to appreciate the damping effect of the standard formulation on |
---|
| 956 | external gravity waves in idealized or weakly non-linear cases. |
---|
| 957 | Although the damping is lower than for the filtered free surface, |
---|
| 958 | it is still significant as shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave. |
---|
[4560] | 959 | |
---|
[2282] | 960 | %>>>>>=============== |
---|
| 961 | \gmcomment{ %%% copy from griffies Book |
---|
[707] | 962 | |
---|
[2282] | 963 | \textbf{title: Time stepping the barotropic system } |
---|
[707] | 964 | |
---|
[10368] | 965 | Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$. |
---|
| 966 | Hence, we can update the surface height and vertically integrated velocity with a leap-frog scheme using |
---|
| 967 | the small barotropic time step $\rdt$. |
---|
| 968 | We have |
---|
[707] | 969 | |
---|
[10419] | 970 | \[ |
---|
| 971 | % \label{eq:DYN_spg_ts_eta} |
---|
| 972 | \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1}) |
---|
| 973 | = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right] |
---|
| 974 | \] |
---|
| 975 | \begin{multline*} |
---|
| 976 | % \label{eq:DYN_spg_ts_u} |
---|
| 977 | \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}) \\ |
---|
| 978 | = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n}) |
---|
| 979 | - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right] |
---|
| 980 | \end{multline*} |
---|
[2282] | 981 | \ |
---|
[707] | 982 | |
---|
[10368] | 983 | In these equations, araised (b) denotes values of surface height and vertically integrated velocity updated with |
---|
| 984 | the barotropic time steps. |
---|
| 985 | The $\tau$ time label on $\eta^{(b)}$ and $U^{(b)}$ denotes the baroclinic time at which |
---|
| 986 | the vertically integrated forcing $\textbf{M}(\tau)$ |
---|
| 987 | (note that this forcing includes the surface freshwater forcing), |
---|
| 988 | the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$, |
---|
| 989 | and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over |
---|
| 990 | a single cycle. |
---|
| 991 | This is also the time that sets the barotropic time steps via |
---|
[10419] | 992 | \[ |
---|
| 993 | % \label{eq:DYN_spg_ts_t} |
---|
| 994 | t_n=\tau+n\rdt |
---|
| 995 | \] |
---|
[10368] | 996 | with $n$ an integer. |
---|
| 997 | The density scaled surface pressure is evaluated via |
---|
[10419] | 998 | \[ |
---|
| 999 | % \label{eq:DYN_spg_ts_ps} |
---|
| 1000 | p_s^{(b)}(\tau,t_{n}) = |
---|
| 1001 | \begin{cases} |
---|
| 1002 | g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o & \text{non-linear case} \\ |
---|
| 1003 | g \;\eta_s^{(b)}(\tau,t_{n}) & \text{linear case} |
---|
| 1004 | \end{cases} |
---|
| 1005 | \] |
---|
[2282] | 1006 | To get started, we assume the following initial conditions |
---|
[10419] | 1007 | \[ |
---|
| 1008 | % \label{eq:DYN_spg_ts_eta} |
---|
| 1009 | \begin{split} |
---|
| 1010 | \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} \\ |
---|
| 1011 | \eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} |
---|
| 1012 | \end{split} |
---|
| 1013 | \] |
---|
[2282] | 1014 | with |
---|
[10419] | 1015 | \[ |
---|
| 1016 | % \label{eq:DYN_spg_ts_etaF} |
---|
| 1017 | \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n}) |
---|
| 1018 | \] |
---|
[10368] | 1019 | the time averaged surface height taken from the previous barotropic cycle. |
---|
| 1020 | Likewise, |
---|
[10419] | 1021 | \[ |
---|
| 1022 | % \label{eq:DYN_spg_ts_u} |
---|
| 1023 | \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ \\ |
---|
| 1024 | \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} |
---|
| 1025 | \] |
---|
[2282] | 1026 | with |
---|
[10419] | 1027 | \[ |
---|
| 1028 | % \label{eq:DYN_spg_ts_u} |
---|
| 1029 | \overline{\textbf{U}^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n}) |
---|
| 1030 | \] |
---|
[10368] | 1031 | the time averaged vertically integrated transport. |
---|
| 1032 | Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration. |
---|
[707] | 1033 | |
---|
[10368] | 1034 | Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ , |
---|
| 1035 | the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at |
---|
| 1036 | baroclinic time $\tau + \rdt \tau$ |
---|
[10419] | 1037 | \[ |
---|
| 1038 | % \label{eq:DYN_spg_ts_u} |
---|
| 1039 | \textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n}) |
---|
| 1040 | \] |
---|
[10368] | 1041 | The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using |
---|
| 1042 | the following form |
---|
[2282] | 1043 | |
---|
[10419] | 1044 | \begin{equation} |
---|
| 1045 | \label{eq:DYN_spg_ts_ssh} |
---|
| 1046 | \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] |
---|
[2282] | 1047 | \end{equation} |
---|
| 1048 | |
---|
[10368] | 1049 | The use of this "big-leap-frog" scheme for the surface height ensures compatibility between |
---|
| 1050 | the mass/volume budgets and the tracer budgets. |
---|
| 1051 | More discussion of this point is provided in Chapter 10 (see in particular Section 10.2). |
---|
[2282] | 1052 | |
---|
[10368] | 1053 | In general, some form of time filter is needed to maintain integrity of the surface height field due to |
---|
| 1054 | the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}. |
---|
| 1055 | We have tried various forms of such filtering, |
---|
| 1056 | with the following method discussed in \cite{Griffies_al_MWR01} chosen due to |
---|
| 1057 | its stability and reasonably good maintenance of tracer conservation properties (see ??). |
---|
[2282] | 1058 | |
---|
[10419] | 1059 | \begin{equation} |
---|
| 1060 | \label{eq:DYN_spg_ts_sshf} |
---|
| 1061 | \eta^{F}(\tau-\Delta) = \overline{\eta^{(b)}(\tau)} |
---|
[2282] | 1062 | \end{equation} |
---|
| 1063 | Another approach tried was |
---|
| 1064 | |
---|
[10419] | 1065 | \[ |
---|
| 1066 | % \label{eq:DYN_spg_ts_sshf2} |
---|
| 1067 | \eta^{F}(\tau-\Delta) = \eta(\tau) |
---|
| 1068 | + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt) |
---|
| 1069 | + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right] |
---|
| 1070 | \] |
---|
[2282] | 1071 | |
---|
[10368] | 1072 | which is useful since it isolates all the time filtering aspects into the term multiplied by $\alpha$. |
---|
| 1073 | This isolation allows for an easy check that tracer conservation is exact when |
---|
| 1074 | eliminating tracer and surface height time filtering (see ?? for more complete discussion). |
---|
| 1075 | However, in the general case with a non-zero $\alpha$, |
---|
| 1076 | the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended. |
---|
[2282] | 1077 | |
---|
| 1078 | } %%end gm comment (copy of griffies book) |
---|
| 1079 | |
---|
| 1080 | %>>>>>=============== |
---|
| 1081 | |
---|
| 1082 | |
---|
[707] | 1083 | %-------------------------------------------------------------------------------------------------------------- |
---|
[2282] | 1084 | % Filtered free surface formulation |
---|
[707] | 1085 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9363] | 1086 | \subsection{Filtered free surface (\protect\key{dynspg\_flt})} |
---|
[9407] | 1087 | \label{subsec:DYN_spg_fltp} |
---|
[707] | 1088 | |
---|
[2282] | 1089 | The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation. |
---|
[9407] | 1090 | The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly. |
---|
| 1091 | The elliptic solvers available in the code are documented in \autoref{chap:MISC}. |
---|
[707] | 1092 | |
---|
[2282] | 1093 | %% gm %%======>>>> given here the discrete eqs provided to the solver |
---|
[2541] | 1094 | \gmcomment{ %%% copy from chap-model basics |
---|
[10419] | 1095 | \[ |
---|
| 1096 | % \label{eq:spg_flt} |
---|
| 1097 | \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} |
---|
| 1098 | - g \nabla \left( \tilde{\rho} \ \eta \right) |
---|
| 1099 | - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) |
---|
| 1100 | \] |
---|
| 1101 | where $T_c$, is a parameter with dimensions of time which characterizes the force, |
---|
| 1102 | $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, |
---|
| 1103 | and $\rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, |
---|
| 1104 | non-linear and viscous terms in \autoref{eq:PE_dyn}. |
---|
[2541] | 1105 | } %end gmcomment |
---|
[707] | 1106 | |
---|
[10368] | 1107 | Note that in the linear free surface formulation (\key{vvl} not defined), |
---|
| 1108 | the ocean depth is time-independent and so is the matrix to be inverted. |
---|
| 1109 | It is computed once and for all and applies to all ocean time steps. |
---|
[707] | 1110 | |
---|
| 1111 | % ================================================================ |
---|
| 1112 | % Lateral diffusion term |
---|
| 1113 | % ================================================================ |
---|
[9393] | 1114 | \section{Lateral diffusion term and operators (\protect\mdl{dynldf})} |
---|
[9407] | 1115 | \label{sec:DYN_ldf} |
---|
[707] | 1116 | %------------------------------------------nam_dynldf---------------------------------------------------- |
---|
[10146] | 1117 | |
---|
| 1118 | \nlst{namdyn_ldf} |
---|
[707] | 1119 | %------------------------------------------------------------------------------------------------------------- |
---|
| 1120 | |
---|
[4147] | 1121 | Options are defined through the \ngn{namdyn\_ldf} namelist variables. |
---|
[10368] | 1122 | The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators. |
---|
| 1123 | The coefficients may be constant or spatially variable; |
---|
| 1124 | the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}). |
---|
| 1125 | The lateral diffusion of momentum is evaluated using a forward scheme, |
---|
| 1126 | $i.e.$ the velocity appearing in its expression is the \textit{before} velocity in time, |
---|
| 1127 | except for the pure vertical component that appears when a tensor of rotation is used. |
---|
| 1128 | This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). |
---|
[817] | 1129 | |
---|
[10368] | 1130 | At the lateral boundaries either free slip, |
---|
| 1131 | no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}). |
---|
[707] | 1132 | |
---|
[6289] | 1133 | \gmcomment{ |
---|
[10368] | 1134 | Hyperviscous operators are frequently used in the simulation of turbulent flows to |
---|
| 1135 | control the dissipation of unresolved small scale features. |
---|
| 1136 | Their primary role is to provide strong dissipation at the smallest scale supported by |
---|
| 1137 | the grid while minimizing the impact on the larger scale features. |
---|
| 1138 | Hyperviscous operators are thus designed to be more scale selective than the traditional, |
---|
| 1139 | physically motivated Laplace operator. |
---|
| 1140 | In finite difference methods, |
---|
| 1141 | the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since |
---|
| 1142 | its damping time ($i.e.$ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$ |
---|
| 1143 | (so that short waves damped more rapidelly than long ones), |
---|
| 1144 | whereas the Laplace operator damping time scales only like $\lambda^{-2}$. |
---|
[6289] | 1145 | } |
---|
| 1146 | |
---|
[707] | 1147 | % ================================================================ |
---|
[9393] | 1148 | \subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})] |
---|
| 1149 | {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})} |
---|
[9407] | 1150 | \label{subsec:DYN_ldf_lap} |
---|
[707] | 1151 | |
---|
| 1152 | For lateral iso-level diffusion, the discrete operator is: |
---|
[10419] | 1153 | \begin{equation} |
---|
| 1154 | \label{eq:dynldf_lap} |
---|
| 1155 | \left\{ |
---|
| 1156 | \begin{aligned} |
---|
| 1157 | D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} |
---|
| 1158 | \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[ |
---|
| 1159 | {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\ |
---|
| 1160 | D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} |
---|
| 1161 | \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[ |
---|
| 1162 | {A_f^{lm} \;e_{3f} \zeta } \right] |
---|
| 1163 | \end{aligned} |
---|
| 1164 | \right. |
---|
[707] | 1165 | \end{equation} |
---|
| 1166 | |
---|
[10368] | 1167 | As explained in \autoref{subsec:PE_ldf}, |
---|
| 1168 | this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and |
---|
| 1169 | ensures a complete separation between the vorticity and divergence parts of the momentum diffusion. |
---|
[707] | 1170 | |
---|
| 1171 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1172 | % Rotated laplacian operator |
---|
| 1173 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 1174 | \subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})] |
---|
| 1175 | {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})} |
---|
[9407] | 1176 | \label{subsec:DYN_ldf_iso} |
---|
[707] | 1177 | |
---|
[10368] | 1178 | A rotation of the lateral momentum diffusion operator is needed in several cases: |
---|
| 1179 | for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}\forcode{ = .true.}) and |
---|
| 1180 | for either iso-neutral (\np{ln\_dynldf\_iso}\forcode{ = .true.}) or |
---|
| 1181 | geopotential (\np{ln\_dynldf\_hor}\forcode{ = .true.}) diffusion in the $s$-coordinate. |
---|
| 1182 | In the partial step case, coordinates are horizontal except at the deepest level and |
---|
| 1183 | no rotation is performed when \np{ln\_dynldf\_hor}\forcode{ = .true.}. |
---|
| 1184 | The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on |
---|
| 1185 | each momentum component. |
---|
| 1186 | It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry. |
---|
| 1187 | The resulting discrete representation is: |
---|
[10419] | 1188 | \begin{equation} |
---|
| 1189 | \label{eq:dyn_ldf_iso} |
---|
| 1190 | \begin{split} |
---|
| 1191 | D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ |
---|
| 1192 | & \left\{\quad {\delta_{i+1/2} \left[ {A_T^{lm} \left( |
---|
| 1193 | {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u] |
---|
| 1194 | -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}} |
---|
| 1195 | \right)} \right]} \right. \\ |
---|
| 1196 | & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} |
---|
| 1197 | }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f} |
---|
| 1198 | \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}} |
---|
| 1199 | \right)} \right] \\ |
---|
| 1200 | &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline |
---|
| 1201 | {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } |
---|
| 1202 | \right.} \right. \\ |
---|
| 1203 | & \ \qquad \qquad \qquad \quad\ |
---|
| 1204 | - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\ |
---|
| 1205 | & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ |
---|
| 1206 | +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} |
---|
| 1207 | \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\ |
---|
| 1208 | D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\ |
---|
| 1209 | & \left\{\quad {\delta_{i+1/2} \left[ {A_f^{lm} \left( |
---|
| 1210 | {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v] |
---|
| 1211 | -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}} |
---|
| 1212 | \right)} \right]} \right. \\ |
---|
| 1213 | & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} |
---|
| 1214 | }\,\delta_{j} [v] - e_{1t}\, r_{2t} |
---|
| 1215 | \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}} |
---|
| 1216 | \right)} \right] \\ |
---|
| 1217 | & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline |
---|
| 1218 | {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\ |
---|
| 1219 | & \ \qquad \qquad \qquad \quad\ |
---|
| 1220 | - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\ |
---|
| 1221 | & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ |
---|
| 1222 | +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} |
---|
| 1223 | \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} |
---|
| 1224 | \end{split} |
---|
[707] | 1225 | \end{equation} |
---|
[10368] | 1226 | where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and |
---|
| 1227 | the surface of computation ($z$- or $s$-surfaces). |
---|
| 1228 | The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}). |
---|
[707] | 1229 | |
---|
| 1230 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1231 | % Iso-level bilaplacian operator |
---|
| 1232 | %-------------------------------------------------------------------------------------------------------------- |
---|
[9393] | 1233 | \subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})] |
---|
| 1234 | {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})} |
---|
[9407] | 1235 | \label{subsec:DYN_ldf_bilap} |
---|
[707] | 1236 | |
---|
[10368] | 1237 | The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:dynldf_lap} twice. |
---|
| 1238 | It requires an additional assumption on boundary conditions: |
---|
| 1239 | the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen, |
---|
| 1240 | while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}). |
---|
[817] | 1241 | %%% |
---|
| 1242 | \gmcomment{add a remark on the the change in the position of the coefficient} |
---|
| 1243 | %%% |
---|
[707] | 1244 | |
---|
| 1245 | % ================================================================ |
---|
| 1246 | % Vertical diffusion term |
---|
| 1247 | % ================================================================ |
---|
[9393] | 1248 | \section{Vertical diffusion term (\protect\mdl{dynzdf})} |
---|
[9407] | 1249 | \label{sec:DYN_zdf} |
---|
[707] | 1250 | %----------------------------------------------namzdf------------------------------------------------------ |
---|
[10146] | 1251 | |
---|
| 1252 | \nlst{namzdf} |
---|
[707] | 1253 | %------------------------------------------------------------------------------------------------------------- |
---|
| 1254 | |
---|
[4147] | 1255 | Options are defined through the \ngn{namzdf} namelist variables. |
---|
[10368] | 1256 | 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. |
---|
| 1257 | Two time stepping schemes can be used for the vertical diffusion term: |
---|
| 1258 | $(a)$ a forward time differencing scheme |
---|
| 1259 | (\np{ln\_zdfexp}\forcode{ = .true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or |
---|
| 1260 | $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{ = .false.}) |
---|
| 1261 | (see \autoref{chap:STP}). |
---|
| 1262 | Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. |
---|
[707] | 1263 | |
---|
[10368] | 1264 | The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is. |
---|
| 1265 | The vertical diffusion operators given by \autoref{eq:PE_zdf} take the following semi-discrete space form: |
---|
[10419] | 1266 | \[ |
---|
| 1267 | % \label{eq:dynzdf} |
---|
| 1268 | \left\{ |
---|
| 1269 | \begin{aligned} |
---|
| 1270 | D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } |
---|
| 1271 | \ \delta_{k+1/2} [\,u\,] \right] \\ |
---|
| 1272 | \\ |
---|
| 1273 | D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} } |
---|
| 1274 | \ \delta_{k+1/2} [\,v\,] \right] |
---|
| 1275 | \end{aligned} |
---|
| 1276 | \right. |
---|
| 1277 | \] |
---|
[10368] | 1278 | where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients. |
---|
| 1279 | The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}). |
---|
[707] | 1280 | |
---|
[10368] | 1281 | The surface boundary condition on momentum is the stress exerted by the wind. |
---|
| 1282 | At the surface, the momentum fluxes are prescribed as the boundary condition on |
---|
| 1283 | the vertical turbulent momentum fluxes, |
---|
[10419] | 1284 | \begin{equation} |
---|
| 1285 | \label{eq:dynzdf_sbc} |
---|
| 1286 | \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} |
---|
| 1287 | = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v } |
---|
[707] | 1288 | \end{equation} |
---|
[10419] | 1289 | where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in |
---|
[10368] | 1290 | the (\textbf{i},\textbf{j}) coordinate system. |
---|
| 1291 | The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in |
---|
| 1292 | the vertical over the mixed layer depth. |
---|
| 1293 | If the vertical mixing coefficient is small (when no mixed layer scheme is used) |
---|
| 1294 | the surface stress enters only the top model level, as a body force. |
---|
| 1295 | The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}). |
---|
[707] | 1296 | |
---|
[10368] | 1297 | The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation |
---|
| 1298 | (see \autoref{sec:ZDF_bfr}) |
---|
[707] | 1299 | |
---|
| 1300 | % ================================================================ |
---|
| 1301 | % External Forcing |
---|
| 1302 | % ================================================================ |
---|
[9393] | 1303 | \section{External forcings} |
---|
[9407] | 1304 | \label{sec:DYN_forcing} |
---|
[707] | 1305 | |
---|
[10368] | 1306 | Besides the surface and bottom stresses (see the above section) |
---|
| 1307 | which are introduced as boundary conditions on the vertical mixing, |
---|
| 1308 | three other forcings may enter the dynamical equations by affecting the surface pressure gradient. |
---|
[707] | 1309 | |
---|
[10368] | 1310 | (1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \autoref{sec:SBC_apr}), |
---|
| 1311 | the atmospheric pressure is taken into account when computing the surface pressure gradient. |
---|
[707] | 1312 | |
---|
[10368] | 1313 | (2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \autoref{sec:SBC_tide}), |
---|
[6289] | 1314 | the tidal potential is taken into account when computing the surface pressure gradient. |
---|
[2285] | 1315 | |
---|
[10368] | 1316 | (3) When \np{nn\_ice\_embd}\forcode{ = 2} and LIM or CICE is used |
---|
| 1317 | ($i.e.$ when the sea-ice is embedded in the ocean), |
---|
[6289] | 1318 | the snow-ice mass is taken into account when computing the surface pressure gradient. |
---|
| 1319 | |
---|
| 1320 | |
---|
| 1321 | \gmcomment{ missing : the lateral boundary condition !!! another external forcing |
---|
| 1322 | } |
---|
| 1323 | |
---|
[707] | 1324 | % ================================================================ |
---|
| 1325 | % Time evolution term |
---|
| 1326 | % ================================================================ |
---|
[9393] | 1327 | \section{Time evolution term (\protect\mdl{dynnxt})} |
---|
[9407] | 1328 | \label{sec:DYN_nxt} |
---|
[707] | 1329 | |
---|
| 1330 | %----------------------------------------------namdom---------------------------------------------------- |
---|
[10146] | 1331 | |
---|
| 1332 | \nlst{namdom} |
---|
[707] | 1333 | %------------------------------------------------------------------------------------------------------------- |
---|
| 1334 | |
---|
[4147] | 1335 | Options are defined through the \ngn{namdom} namelist variables. |
---|
[10368] | 1336 | The general framework for dynamics time stepping is a leap-frog scheme, |
---|
| 1337 | $i.e.$ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:STP}). |
---|
| 1338 | The scheme is applied to the velocity, except when |
---|
| 1339 | using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) |
---|
| 1340 | in the variable volume case (\key{vvl} defined), |
---|
| 1341 | where it has to be applied to the thickness weighted velocity (see \autoref{sec:A_momentum}) |
---|
[707] | 1342 | |
---|
[10368] | 1343 | $\bullet$ vector invariant form or linear free surface |
---|
| 1344 | (\np{ln\_dynhpg\_vec}\forcode{ = .true.} ; \key{vvl} not defined): |
---|
[10419] | 1345 | \[ |
---|
| 1346 | % \label{eq:dynnxt_vec} |
---|
| 1347 | \left\{ |
---|
| 1348 | \begin{aligned} |
---|
| 1349 | &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt \ \text{RHS}_u^t \\ |
---|
| 1350 | &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right] |
---|
| 1351 | \end{aligned} |
---|
| 1352 | \right. |
---|
| 1353 | \] |
---|
[707] | 1354 | |
---|
[10368] | 1355 | $\bullet$ flux form and nonlinear free surface |
---|
| 1356 | (\np{ln\_dynhpg\_vec}\forcode{ = .false.} ; \key{vvl} defined): |
---|
[10419] | 1357 | \[ |
---|
| 1358 | % \label{eq:dynnxt_flux} |
---|
| 1359 | \left\{ |
---|
| 1360 | \begin{aligned} |
---|
| 1361 | &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t \\ |
---|
| 1362 | &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t |
---|
| 1363 | +\gamma \,\left[ {\left(e_{3u}\,u\right)_f^{t-\rdt} -2\left(e_{3u}\,u\right)^t+\left(e_{3u}\,u\right)^{t+\rdt}} \right] |
---|
| 1364 | \end{aligned} |
---|
| 1365 | \right. |
---|
| 1366 | \] |
---|
[10368] | 1367 | where RHS is the right hand side of the momentum equation, |
---|
| 1368 | the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient. |
---|
| 1369 | $\gamma$ is initialized as \np{nn\_atfp} (namelist parameter). |
---|
| 1370 | Its default value is \np{nn\_atfp}\forcode{ = 10.e-3}. |
---|
| 1371 | In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for |
---|
| 1372 | the momentum equations. |
---|
[707] | 1373 | |
---|
[10368] | 1374 | Note that with the filtered free surface, |
---|
| 1375 | the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module, |
---|
| 1376 | and only array swapping and Asselin filtering is done in \mdl{dynnxt}. |
---|
[707] | 1377 | |
---|
[1224] | 1378 | % ================================================================ |
---|
[10419] | 1379 | \biblio |
---|
| 1380 | |
---|
[6997] | 1381 | \end{document} |
---|