[707] | 1 | % ================================================================ |
---|
| 2 | % Chapter Ñ Ocean Dynamics (DYN) |
---|
| 3 | % ================================================================ |
---|
| 4 | \chapter{Ocean Dynamics (DYN)} |
---|
| 5 | \label{DYN} |
---|
| 6 | \minitoc |
---|
| 7 | |
---|
| 8 | % add a figure for dynvor ens, ene latices |
---|
| 9 | |
---|
[1224] | 10 | %\vspace{2.cm} |
---|
[707] | 11 | $\ $\newline %force an empty line |
---|
| 12 | |
---|
[2282] | 13 | Using the representation described in Chapter \ref{DOM}, several semi-discrete |
---|
[817] | 14 | space forms of the dynamical equations are available depending on the vertical |
---|
| 15 | coordinate used and on the conservation properties of the vorticity term. In all |
---|
| 16 | the equations presented here, the masking has been omitted for simplicity. |
---|
[2282] | 17 | One must be aware that all the quantities are masked fields and that each time an |
---|
[817] | 18 | average or difference operator is used, the resulting field is multiplied by a mask. |
---|
[707] | 19 | |
---|
| 20 | The prognostic ocean dynamics equation can be summarized as follows: |
---|
| 21 | \begin{equation*} |
---|
| 22 | \text{NXT} = \dbinom {\text{VOR} + \text{KEG} + \text {ZAD} } |
---|
| 23 | {\text{COR} + \text{ADV} } |
---|
| 24 | + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF} |
---|
| 25 | \end{equation*} |
---|
[817] | 26 | NXT stands for next, referring to the time-stepping. The first group of terms on |
---|
[2285] | 27 | the rhs of this equation corresponds to the Coriolis and advection |
---|
| 28 | terms that are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) |
---|
| 29 | and a vertical advection part (ZAD) in the vector invariant formulation, or a Coriolis |
---|
[2282] | 30 | and advection part (COR+ADV) in the flux formulation. The terms following these |
---|
[817] | 31 | are the pressure gradient contributions (HPG, Hydrostatic Pressure Gradient, |
---|
| 32 | and SPG, Surface Pressure Gradient); and contributions from lateral diffusion |
---|
| 33 | (LDF) and vertical diffusion (ZDF), which are added to the rhs in the \mdl{dynldf} |
---|
| 34 | and \mdl{dynzdf} modules. The vertical diffusion term includes the surface and |
---|
| 35 | bottom stresses. The external forcings and parameterisations require complex |
---|
| 36 | inputs (surface wind stress calculation using bulk formulae, estimation of mixing |
---|
| 37 | coefficients) that are carried out in modules SBC, LDF and ZDF and are described |
---|
| 38 | in Chapters \ref{SBC}, \ref{LDF} and \ref{ZDF}, respectively. |
---|
[707] | 39 | |
---|
[817] | 40 | In the present chapter we also describe the diagnostic equations used to compute |
---|
[2282] | 41 | the horizontal divergence, curl of the velocities (\emph{divcur} module) and |
---|
| 42 | the vertical velocity (\emph{wzvmod} module). |
---|
[707] | 43 | |
---|
[817] | 44 | The different options available to the user are managed by namelist variables. |
---|
[2282] | 45 | For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, |
---|
[817] | 46 | where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. |
---|
| 47 | If a CPP key is used for this term its name is \textbf{key\_ttt}. The corresponding |
---|
| 48 | code can be found in the \textit{dynttt\_xxx} module in the DYN directory, and it is |
---|
| 49 | usually computed in the \textit{dyn\_ttt\_xxx} subroutine. |
---|
[707] | 50 | |
---|
[2282] | 51 | The user has the option of extracting and outputting each tendency term from the |
---|
| 52 | 3D momentum equations (\key{trddyn} defined), as described in |
---|
| 53 | Chap.\ref{MISC}. Furthermore, the tendency terms associated with the 2D |
---|
| 54 | barotropic vorticity balance (when \key{trdvor} is defined) can be derived from the |
---|
[817] | 55 | 3D terms. |
---|
| 56 | %%% |
---|
[996] | 57 | \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does |
---|
| 58 | MISC correspond to "extracting tendency terms" or "vorticity balance"?} |
---|
[707] | 59 | |
---|
[1224] | 60 | $\ $\newline % force a new ligne |
---|
[2282] | 61 | |
---|
[707] | 62 | % ================================================================ |
---|
[2282] | 63 | % Sea Surface Height evolution & Diagnostics variables |
---|
| 64 | % ================================================================ |
---|
| 65 | \section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} |
---|
| 66 | \label{DYN_divcur_wzv} |
---|
| 67 | |
---|
| 68 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 69 | % Horizontal divergence and relative vorticity |
---|
| 70 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 71 | \subsection [Horizontal divergence and relative vorticity (\textit{divcur})] |
---|
| 72 | {Horizontal divergence and relative vorticity (\mdl{divcur})} |
---|
| 73 | \label{DYN_divcur} |
---|
| 74 | |
---|
| 75 | The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: |
---|
| 76 | \begin{equation} \label{Eq_divcur_cur} |
---|
| 77 | \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right] |
---|
| 78 | -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) |
---|
| 79 | \end{equation} |
---|
| 80 | |
---|
| 81 | The horizontal divergence is defined at a $T$-point. It is given by: |
---|
| 82 | \begin{equation} \label{Eq_divcur_div} |
---|
| 83 | \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } |
---|
| 84 | \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right] |
---|
| 85 | +\delta _j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right) |
---|
| 86 | \end{equation} |
---|
| 87 | |
---|
[2285] | 88 | Note that although the vorticity has the same discrete expression in $z$- |
---|
[2282] | 89 | and $s$-coordinates, its physical meaning is not identical. $\zeta$ is a pseudo |
---|
| 90 | vorticity along $s$-surfaces (only pseudo because $(u,v)$ are still defined along |
---|
| 91 | geopotential surfaces, but are not necessarily defined at the same depth). |
---|
| 92 | |
---|
| 93 | The vorticity and divergence at the \textit{before} step are used in the computation |
---|
| 94 | of the horizontal diffusion of momentum. Note that because they have been |
---|
| 95 | calculated prior to the Asselin filtering of the \textit{before} velocities, the |
---|
| 96 | \textit{before} vorticity and divergence arrays must be included in the restart file |
---|
| 97 | to ensure perfect restartability. The vorticity and divergence at the \textit{now} |
---|
| 98 | time step are used for the computation of the nonlinear advection and of the |
---|
| 99 | vertical velocity respectively. |
---|
| 100 | |
---|
| 101 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 102 | % Sea Surface Height evolution |
---|
| 103 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 104 | \subsection [Sea surface height evolution and vertical velocity (\textit{sshwzv})] |
---|
| 105 | {Horizontal divergence and relative vorticity (\mdl{sshwzv})} |
---|
| 106 | \label{DYN_sshwzv} |
---|
| 107 | |
---|
| 108 | The sea surface height is given by : |
---|
| 109 | \begin{equation} \label{Eq_dynspg_ssh} |
---|
| 110 | \begin{aligned} |
---|
| 111 | \frac{\partial \eta }{\partial t} |
---|
[2285] | 112 | &\equiv \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left\{ \delta _i \left[ {e_{2u}\,e_{3u}\;u} \right] |
---|
| 113 | +\delta _j \left[ {e_{1v}\,e_{3v}\;v} \right] \right\} } |
---|
[2282] | 114 | - \frac{\textit{emp}}{\rho _w } \\ |
---|
| 115 | &\equiv \sum\limits_k {\chi \ e_{3t}} - \frac{\textit{emp}}{\rho _w } |
---|
| 116 | \end{aligned} |
---|
| 117 | \end{equation} |
---|
| 118 | where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), |
---|
[2285] | 119 | expressed in Kg/m$^2$/s (which is equal to mm/s), and $\rho _w$=1,035~Kg/m$^3$ |
---|
| 120 | is the reference density of sea water (Boussinesq approximation). If river runoff is |
---|
| 121 | expressed as a surface freshwater flux (see \S\ref{SBC}) then \textit{emp} can be |
---|
| 122 | written as the evaporation minus precipitation, minus the river runoff. |
---|
| 123 | The sea-surface height is evaluated using exactly the same time stepping scheme |
---|
| 124 | as the tracer equation \eqref{Eq_tra_nxt}: |
---|
[2282] | 125 | a leapfrog scheme in combination with an Asselin time filter, $i.e.$ the velocity appearing |
---|
| 126 | in \eqref{Eq_dynspg_ssh} is centred in time (\textit{now} velocity). |
---|
| 127 | This is of paramount importance. Replacing $T$ by the number $1$ in the tracer equation and summing |
---|
| 128 | over the water column must lead to the sea surface height equation otherwise tracer content |
---|
| 129 | will not be conserved \ref{Griffies_al_MWR01, LeclairMadec2009}. |
---|
| 130 | |
---|
| 131 | The vertical velocity is computed by an upward integration of the horizontal |
---|
| 132 | divergence starting at the bottom, taking into account the change of the thickness of the levels : |
---|
| 133 | \begin{equation} \label{Eq_wzv} |
---|
| 134 | \left\{ \begin{aligned} |
---|
[2285] | 135 | &\left. w \right|_{k_b-1/2} \quad= 0 \qquad \text{where } k_b \text{ is the level just above the sea floor } \\ |
---|
| 136 | &\left. w \right|_{k+1/2} = \left. w \right|_{k-1/2} + \left. e_{3t} \right|_{k}\; \left. \chi \right|_k |
---|
| 137 | - \frac{1} {2 \rdt} \left( \left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k}\right) |
---|
[2282] | 138 | \end{aligned} \right. |
---|
| 139 | \end{equation} |
---|
| 140 | |
---|
| 141 | In the case of a non-linear free surface (\key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$, |
---|
| 142 | as changes in the divergence of the barotropic transport are absorbed into the change |
---|
| 143 | of the level thicknesses, re-orientated downward. |
---|
[2285] | 144 | \gmcomment{not sure of this... to be modified with the change in emp setting} |
---|
[2282] | 145 | In the case of a linear free surface, the time derivative in \eqref{Eq_wzv} disappears. |
---|
| 146 | The upper boundary condition applies at a fixed level $z=0$. The top vertical velocity |
---|
| 147 | is thus equal to the divergence of the barotropic transport ($i.e.$ the first term in the |
---|
| 148 | right-hand-side of \eqref{Eq_dynspg_ssh}). |
---|
| 149 | |
---|
| 150 | Note also that whereas the vertical velocity has the same discrete |
---|
| 151 | expression in $z$- and $s$-coordinates, its physical meaning is not the same: |
---|
| 152 | in the second case, $w$ is the velocity normal to the $s$-surfaces. |
---|
| 153 | Note also that the $k$-axis is re-orientated downwards in the \textsc{fortran} code compared |
---|
| 154 | to the indexing used in the semi-discrete equations such as \eqref{Eq_wzv} |
---|
| 155 | (see \S\ref{DOM_Num_Index_vertical}). |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | % ================================================================ |
---|
[707] | 159 | % Coriolis and Advection terms: vector invariant form |
---|
| 160 | % ================================================================ |
---|
| 161 | \section{Coriolis and Advection: vector invariant form} |
---|
| 162 | \label{DYN_adv_cor_vect} |
---|
| 163 | %-----------------------------------------nam_dynadv---------------------------------------------------- |
---|
[2282] | 164 | \namdisplay{namdyn_adv} |
---|
[707] | 165 | %------------------------------------------------------------------------------------------------------------- |
---|
| 166 | |
---|
[817] | 167 | The vector invariant form of the momentum equations is the one most |
---|
[2282] | 168 | often used in applications of the \NEMO ocean model. The flux form option |
---|
| 169 | (see next section) has been present since version $2$. |
---|
[1224] | 170 | Coriolis and momentum advection terms are evaluated using a leapfrog |
---|
| 171 | scheme, $i.e.$ the velocity appearing in these expressions is centred in |
---|
| 172 | time (\textit{now} velocity). |
---|
[817] | 173 | At the lateral boundaries either free slip, no slip or partial slip boundary |
---|
| 174 | conditions are applied following Chap.\ref{LBC}. |
---|
[707] | 175 | |
---|
| 176 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 177 | % Vorticity term |
---|
| 178 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 179 | \subsection [Vorticity term (\textit{dynvor}) ] |
---|
| 180 | {Vorticity term (\mdl{dynvor})} |
---|
[707] | 181 | \label{DYN_vor} |
---|
| 182 | %------------------------------------------nam_dynvor---------------------------------------------------- |
---|
[2282] | 183 | \namdisplay{namdyn_vor} |
---|
[707] | 184 | %------------------------------------------------------------------------------------------------------------- |
---|
| 185 | |
---|
[2282] | 186 | Four discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=true) are available: |
---|
| 187 | conserving potential enstrophy of horizontally non-divergent flow (ENS scheme) ; |
---|
| 188 | conserving horizontal kinetic energy (ENE scheme) ; conserving potential enstrophy for |
---|
| 189 | the relative vorticity term and horizontal kinetic energy for the planetary vorticity |
---|
| 190 | term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent |
---|
[3294] | 191 | flow and horizontal kinetic energy (EEN scheme) (see Appendix~\ref{Apdx_C_vor_zad}). In the |
---|
| 192 | case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the |
---|
| 193 | consistency of vorticity term with analytical equations (\textit{ln\_dynvor\_con}=true). |
---|
[2285] | 194 | The vorticity terms are all computed in dedicated routines that can be found in |
---|
[2282] | 195 | the \mdl{dynvor} module. |
---|
[707] | 196 | |
---|
| 197 | %------------------------------------------------------------- |
---|
| 198 | % enstrophy conserving scheme |
---|
| 199 | %------------------------------------------------------------- |
---|
[2282] | 200 | \subsubsection{Enstrophy conserving scheme (\np{ln\_dynvor\_ens}=true)} |
---|
[707] | 201 | \label{DYN_vor_ens} |
---|
| 202 | |
---|
[817] | 203 | In the enstrophy conserving case (ENS scheme), the discrete formulation of the |
---|
| 204 | vorticity term provides a global conservation of the enstrophy |
---|
| 205 | ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent |
---|
[2282] | 206 | flow ($i.e.$ $\chi$=$0$), but does not conserve the total kinetic energy. It is given by: |
---|
[707] | 207 | \begin{equation} \label{Eq_dynvor_ens} |
---|
| 208 | \left\{ |
---|
| 209 | \begin{aligned} |
---|
[1224] | 210 | {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} |
---|
[2282] | 211 | & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2} \\ |
---|
[1224] | 212 | {- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} |
---|
[2282] | 213 | & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j} |
---|
[707] | 214 | \end{aligned} |
---|
| 215 | \right. |
---|
| 216 | \end{equation} |
---|
| 217 | |
---|
| 218 | %------------------------------------------------------------- |
---|
| 219 | % energy conserving scheme |
---|
| 220 | %------------------------------------------------------------- |
---|
[2282] | 221 | \subsubsection{Energy conserving scheme (\np{ln\_dynvor\_ene}=true)} |
---|
[707] | 222 | \label{DYN_vor_ene} |
---|
| 223 | |
---|
[817] | 224 | The kinetic energy conserving scheme (ENE scheme) conserves the global |
---|
| 225 | kinetic energy but not the global enstrophy. It is given by: |
---|
[707] | 226 | \begin{equation} \label{Eq_dynvor_ene} |
---|
[1224] | 227 | \left\{ \begin{aligned} |
---|
| 228 | {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
[2282] | 229 | \; \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
[1224] | 230 | {- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
[2282] | 231 | \; \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } |
---|
[1224] | 232 | \end{aligned} \right. |
---|
[707] | 233 | \end{equation} |
---|
| 234 | |
---|
| 235 | %------------------------------------------------------------- |
---|
| 236 | % mix energy/enstrophy conserving scheme |
---|
| 237 | %------------------------------------------------------------- |
---|
[2282] | 238 | \subsubsection{Mixed energy/enstrophy conserving scheme (\np{ln\_dynvor\_mix}=true) } |
---|
[707] | 239 | \label{DYN_vor_mix} |
---|
| 240 | |
---|
[2282] | 241 | For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the |
---|
[817] | 242 | two previous schemes is used. It consists of the ENS scheme (\ref{Eq_dynvor_ens}) |
---|
[2282] | 243 | for the relative vorticity term, and of the ENE scheme (\ref{Eq_dynvor_ene}) applied |
---|
[817] | 244 | to the planetary vorticity term. |
---|
[707] | 245 | \begin{equation} \label{Eq_dynvor_mix} |
---|
[2282] | 246 | \left\{ { \begin{aligned} |
---|
[994] | 247 | {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} |
---|
[2282] | 248 | \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } |
---|
[707] | 249 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
[2282] | 250 | \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
[994] | 251 | {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j |
---|
[2282] | 252 | \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } |
---|
[707] | 253 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
[2282] | 254 | \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill |
---|
| 255 | \end{aligned} } \right. |
---|
[707] | 256 | \end{equation} |
---|
| 257 | |
---|
| 258 | %------------------------------------------------------------- |
---|
| 259 | % energy and enstrophy conserving scheme |
---|
| 260 | %------------------------------------------------------------- |
---|
[2282] | 261 | \subsubsection{Energy and enstrophy conserving scheme (\np{ln\_dynvor\_een}=true) } |
---|
[707] | 262 | \label{DYN_vor_een} |
---|
| 263 | |
---|
[2282] | 264 | In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$ |
---|
| 265 | averages of the velocity allows for the presence of grid point oscillation structures |
---|
| 266 | that will be invisible to the operator. These structures are \textit{computational modes} |
---|
| 267 | that will be at least partly damped by the momentum diffusion operator ($i.e.$ the |
---|
| 268 | subgrid-scale advection), but not by the resolved advection term. The ENS and ENE schemes |
---|
[2285] | 269 | therefore do not contribute to dump any grid point noise in the horizontal velocity field. |
---|
| 270 | Such noise would result in more noise in the vertical velocity field, an undesirable feature. |
---|
| 271 | This is a well-known characteristic of $C$-grid discretization where $u$ and $v$ are located |
---|
| 272 | at different grid points, a price worth paying to avoid a double averaging in the pressure |
---|
| 273 | gradient term as in the $B$-grid. |
---|
[2282] | 274 | \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) |
---|
| 275 | Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} |
---|
| 276 | |
---|
[2285] | 277 | A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}. |
---|
| 278 | The idea is to get rid of the double averaging by considering triad combinations of vorticity. |
---|
[2282] | 279 | It is noteworthy that this solution is conceptually quite similar to the one proposed by |
---|
[2285] | 280 | \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see App.\ref{Apdx_C}). |
---|
[2282] | 281 | |
---|
| 282 | The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified |
---|
| 283 | for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme. |
---|
| 284 | First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: |
---|
[707] | 285 | \begin{equation} \label{Eq_pot_vor} |
---|
[2282] | 286 | q = \frac{\zeta +f} {e_{3f} } |
---|
[707] | 287 | \end{equation} |
---|
[817] | 288 | where the relative vorticity is defined by (\ref{Eq_divcur_cur}), the Coriolis parameter |
---|
| 289 | is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: |
---|
[707] | 290 | \begin{equation} \label{Eq_een_e3f} |
---|
| 291 | e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} |
---|
| 292 | \end{equation} |
---|
| 293 | |
---|
| 294 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
[2376] | 295 | \begin{figure}[!ht] \begin{center} |
---|
[998] | 296 | \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_DYN_een_triad.pdf} |
---|
[2376] | 297 | \caption{ \label{Fig_DYN_een_triad} |
---|
| 298 | Triads used in the energy and enstrophy conserving scheme (een) for |
---|
[817] | 299 | $u$-component (upper panel) and $v$-component (lower panel).} |
---|
[2376] | 300 | \end{center} \end{figure} |
---|
[707] | 301 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 302 | |
---|
[2282] | 303 | Note that a key point in \eqref{Eq_een_e3f} is that the averaging in the \textbf{i}- and |
---|
[817] | 304 | \textbf{j}- directions uses the masked vertical scale factor but is always divided by |
---|
[2282] | 305 | $4$, not by the sum of the masks at the four $T$-points. This preserves the continuity of |
---|
| 306 | $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and |
---|
| 307 | extends by continuity the value of $e_{3f}$ into the land areas. This feature is essential for |
---|
| 308 | the $z$-coordinate with partial steps. |
---|
[707] | 309 | |
---|
[2282] | 310 | Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as |
---|
| 311 | the following triad combinations of the neighbouring potential vorticities defined at f-points |
---|
| 312 | (Fig.~\ref{Fig_DYN_een_triad}): |
---|
| 313 | \begin{equation} \label{Q_triads} |
---|
| 314 | _i^j \mathbb{Q}^{i_p}_{j_p} |
---|
| 315 | = \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) |
---|
| 316 | \end{equation} |
---|
| 317 | 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$. |
---|
| 318 | |
---|
| 319 | Finally, the vorticity terms are represented as: |
---|
[707] | 320 | \begin{equation} \label{Eq_dynvor_een} |
---|
| 321 | \left\{ { |
---|
| 322 | \begin{aligned} |
---|
[2282] | 323 | +q\,e_3 \, v &\equiv +\frac{1}{e_{1u} } \sum_{\substack{i_p,\,k_p}} |
---|
| 324 | {^{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} \\ |
---|
| 325 | - q\,e_3 \, u &\equiv -\frac{1}{e_{2v} } \sum_{\substack{i_p,\,k_p}} |
---|
| 326 | {^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} \\ |
---|
[707] | 327 | \end{aligned} |
---|
| 328 | } \right. |
---|
| 329 | \end{equation} |
---|
| 330 | |
---|
[2282] | 331 | This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. |
---|
| 332 | It conserves both total energy and potential enstrophy in the limit of horizontally |
---|
| 333 | nondivergent flow ($i.e.$ $\chi$=$0$) (see Appendix~\ref{Apdx_C_vor_zad}). |
---|
[2285] | 334 | Applied to a realistic ocean configuration, it has been shown that it leads to a significant |
---|
| 335 | reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}. |
---|
[2282] | 336 | Furthermore, used in combination with a partial steps representation of bottom topography, |
---|
| 337 | it improves the interaction between current and topography, leading to a larger |
---|
| 338 | topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}. |
---|
| 339 | |
---|
[707] | 340 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 341 | % Kinetic Energy Gradient term |
---|
| 342 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 343 | \subsection [Kinetic Energy Gradient term (\textit{dynkeg})] |
---|
| 344 | {Kinetic Energy Gradient term (\mdl{dynkeg})} |
---|
[707] | 345 | \label{DYN_keg} |
---|
| 346 | |
---|
[2282] | 347 | As demonstrated in Appendix~\ref{Apdx_C}, there is a single discrete formulation |
---|
[817] | 348 | of the kinetic energy gradient term that, together with the formulation chosen for |
---|
| 349 | the vertical advection (see below), conserves the total kinetic energy: |
---|
[707] | 350 | \begin{equation} \label{Eq_dynkeg} |
---|
| 351 | \left\{ \begin{aligned} |
---|
[2282] | 352 | -\frac{1}{2 \; e_{1u} } & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] \\ |
---|
| 353 | -\frac{1}{2 \; e_{2v} } & \ \delta _{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] |
---|
[707] | 354 | \end{aligned} \right. |
---|
| 355 | \end{equation} |
---|
| 356 | |
---|
| 357 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 358 | % Vertical advection term |
---|
| 359 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 360 | \subsection [Vertical advection term (\textit{dynzad}) ] |
---|
| 361 | {Vertical advection term (\mdl{dynzad}) } |
---|
[707] | 362 | \label{DYN_zad} |
---|
| 363 | |
---|
[817] | 364 | The discrete formulation of the vertical advection, together with the formulation |
---|
| 365 | chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic |
---|
| 366 | energy. Indeed, the change of KE due to the vertical advection is exactly |
---|
| 367 | balanced by the change of KE due to the gradient of KE (see Appendix~\ref{Apdx_C}). |
---|
[707] | 368 | \begin{equation} \label{Eq_dynzad} |
---|
| 369 | \left\{ \begin{aligned} |
---|
[2282] | 370 | -\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2} \;\delta _{k+1/2} \left[ u \right]\ }^{\,k} \\ |
---|
| 371 | -\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} |
---|
| 372 | \end{aligned} \right. |
---|
[707] | 373 | \end{equation} |
---|
| 374 | |
---|
| 375 | % ================================================================ |
---|
| 376 | % Coriolis and Advection : flux form |
---|
| 377 | % ================================================================ |
---|
| 378 | \section{Coriolis and Advection: flux form} |
---|
| 379 | \label{DYN_adv_cor_flux} |
---|
| 380 | %------------------------------------------nam_dynadv---------------------------------------------------- |
---|
[2282] | 381 | \namdisplay{namdyn_adv} |
---|
[707] | 382 | %------------------------------------------------------------------------------------------------------------- |
---|
| 383 | |
---|
[817] | 384 | In the flux form (as in the vector invariant form), the Coriolis and momentum |
---|
| 385 | advection terms are evaluated using a leapfrog scheme, $i.e.$ the velocity |
---|
| 386 | appearing in their expressions is centred in time (\textit{now} velocity). At the |
---|
| 387 | lateral boundaries either free slip, no slip or partial slip boundary conditions |
---|
| 388 | are applied following Chap.\ref{LBC}. |
---|
[707] | 389 | |
---|
| 390 | |
---|
| 391 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 392 | % Coriolis plus curvature metric terms |
---|
| 393 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 394 | \subsection [Coriolis plus curvature metric terms (\textit{dynvor}) ] |
---|
| 395 | {Coriolis plus curvature metric terms (\mdl{dynvor}) } |
---|
[707] | 396 | \label{DYN_cor_flux} |
---|
| 397 | |
---|
[817] | 398 | In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis |
---|
| 399 | parameter has been modified to account for the "metric" term. This altered |
---|
| 400 | Coriolis parameter is thus discretised at $f$-points. It is given by: |
---|
[707] | 401 | \begin{multline} \label{Eq_dyncor_metric} |
---|
| 402 | f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right) \\ |
---|
[2282] | 403 | \equiv f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right] |
---|
| 404 | - \overline u ^{j+1/2}\delta _{j+1/2} \left[ {e_{1u} } \right] } \ \right) |
---|
[707] | 405 | \end{multline} |
---|
| 406 | |
---|
[817] | 407 | Any of the (\ref{Eq_dynvor_ens}), (\ref{Eq_dynvor_ene}) and (\ref{Eq_dynvor_een}) |
---|
| 408 | schemes can be used to compute the product of the Coriolis parameter and the |
---|
| 409 | vorticity. However, the energy-conserving scheme (\ref{Eq_dynvor_een}) has |
---|
| 410 | exclusively been used to date. This term is evaluated using a leapfrog scheme, |
---|
| 411 | $i.e.$ the velocity is centred in time (\textit{now} velocity). |
---|
[707] | 412 | |
---|
| 413 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 414 | % Flux form Advection term |
---|
| 415 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 416 | \subsection [Flux form Advection term (\textit{dynadv}) ] |
---|
| 417 | {Flux form Advection term (\mdl{dynadv}) } |
---|
[707] | 418 | \label{DYN_adv_flux} |
---|
| 419 | |
---|
| 420 | The discrete expression of the advection term is given by : |
---|
| 421 | \begin{equation} \label{Eq_dynadv} |
---|
| 422 | \left\{ |
---|
| 423 | \begin{aligned} |
---|
| 424 | \frac{1}{e_{1u}\,e_{2u}\,e_{3u}} |
---|
[2282] | 425 | \left( \delta _{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i } \ u_t \right] |
---|
| 426 | + \delta _{j } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2} \ u_f \right] \right. \ \; \\ |
---|
| 427 | \left. + \delta _{k } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2} \ u_{uw} \right] \right) \\ |
---|
[707] | 428 | \\ |
---|
| 429 | \frac{1}{e_{1v}\,e_{2v}\,e_{3v}} |
---|
[2282] | 430 | \left( \delta _{i } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f \right] |
---|
| 431 | + \delta _{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i } \ v_t \right] \right. \ \, \, \\ |
---|
| 432 | \left. + \delta _{k } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw} \right] \right) \\ |
---|
[707] | 433 | \end{aligned} |
---|
| 434 | \right. |
---|
| 435 | \end{equation} |
---|
| 436 | |
---|
[817] | 437 | Two advection schemes are available: a $2^{nd}$ order centered finite |
---|
| 438 | difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS. |
---|
[2282] | 439 | The latter is described in \citet{Shchepetkin_McWilliams_OM05}. The schemes are |
---|
| 440 | selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. |
---|
| 441 | In flux form, the schemes differ by the choice of a space and time interpolation to |
---|
[817] | 442 | define the value of $u$ and $v$ at the centre of each face of $u$- and $v$-cells, |
---|
| 443 | $i.e.$ at the $T$-, $f$-, and $uw$-points for $u$ and at the $f$-, $T$- and |
---|
| 444 | $vw$-points for $v$. |
---|
[707] | 445 | |
---|
| 446 | %------------------------------------------------------------- |
---|
| 447 | % 2nd order centred scheme |
---|
| 448 | %------------------------------------------------------------- |
---|
[2282] | 449 | \subsubsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_dynadv\_cen2}=true)} |
---|
[707] | 450 | \label{DYN_adv_cen2} |
---|
| 451 | |
---|
| 452 | In the centered $2^{nd}$ order formulation, the velocity is evaluated as the |
---|
| 453 | mean of the two neighbouring points : |
---|
| 454 | \begin{equation} \label{Eq_dynadv_cen2} |
---|
| 455 | \left\{ \begin{aligned} |
---|
[2282] | 456 | u_T^{cen2} &=\overline u^{i } \quad & u_F^{cen2} &=\overline u^{j+1/2} \quad & u_{uw}^{cen2} &=\overline u^{k+1/2} \\ |
---|
| 457 | v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j \quad & v_{vw}^{cen2} &=\overline v ^{k+1/2} \\ |
---|
| 458 | \end{aligned} \right. |
---|
[707] | 459 | \end{equation} |
---|
| 460 | |
---|
[817] | 461 | The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive |
---|
[1224] | 462 | ($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be |
---|
| 463 | used in conjunction with an explicit diffusion operator to produce a sensible solution. |
---|
| 464 | The associated time-stepping is performed using a leapfrog scheme in conjunction |
---|
| 465 | with an Asselin time-filter, so $u$ and $v$ are the \emph{now} velocities. |
---|
[707] | 466 | |
---|
| 467 | %------------------------------------------------------------- |
---|
| 468 | % UBS scheme |
---|
| 469 | %------------------------------------------------------------- |
---|
[2282] | 470 | \subsubsection{Upstream Biased Scheme (UBS) (\np{ln\_dynadv\_ubs}=true)} |
---|
[707] | 471 | \label{DYN_adv_ubs} |
---|
| 472 | |
---|
| 473 | The UBS advection scheme is an upstream biased third order scheme based on |
---|
| 474 | an upstream-biased parabolic interpolation. For example, the evaluation of |
---|
| 475 | $u_T^{ubs} $ is done as follows: |
---|
| 476 | \begin{equation} \label{Eq_dynadv_ubs} |
---|
| 477 | u_T^{ubs} =\overline u ^i-\;\frac{1}{6} \begin{cases} |
---|
| 478 | u"_{i-1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i \geqslant 0$ } \\ |
---|
| 479 | u"_{i+1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i < 0$ } |
---|
| 480 | \end{cases} |
---|
| 481 | \end{equation} |
---|
[817] | 482 | where $u"_{i+1/2} =\delta _{i+1/2} \left[ {\delta _i \left[ u \right]} \right]$. This results |
---|
[2282] | 483 | in a dissipatively dominant ($i.e.$ hyper-diffusive) truncation error \citep{Shchepetkin_McWilliams_OM05}. |
---|
[817] | 484 | The overall performance of the advection scheme is similar to that reported in |
---|
| 485 | \citet{Farrow1995}. It is a relatively good compromise between accuracy and |
---|
| 486 | smoothness. It is not a \emph{positive} scheme, meaning that false extrema are |
---|
| 487 | permitted. But the amplitudes of the false extrema are significantly reduced over |
---|
[2282] | 488 | those in the centred second order method. As the scheme already includes |
---|
| 489 | a diffusion component, it can be used without explicit lateral diffusion on momentum |
---|
| 490 | ($i.e.$ \np{ln\_dynldf\_lap}=\np{ln\_dynldf\_bilap}=false), and it is recommended to do so. |
---|
[707] | 491 | |
---|
[817] | 492 | The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ |
---|
| 493 | order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and |
---|
| 494 | $u_{vw}^{ubs}$ in \eqref{Eq_dynadv_cen2} are used. UBS is diffusive and is |
---|
| 495 | associated with vertical mixing of momentum. \gmcomment{ gm pursue the |
---|
| 496 | sentence:Since vertical mixing of momentum is a source term of the TKE equation... } |
---|
[707] | 497 | |
---|
[817] | 498 | For stability reasons, the first term in (\ref{Eq_dynadv_ubs}), which corresponds |
---|
| 499 | to a second order centred scheme, is evaluated using the \textit{now} velocity |
---|
[2282] | 500 | (centred in time), while the second term, which is the diffusion part of the scheme, |
---|
[817] | 501 | is evaluated using the \textit{before} velocity (forward in time). This is discussed |
---|
[2282] | 502 | by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. |
---|
[707] | 503 | |
---|
[2282] | 504 | Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) |
---|
| 505 | schemes only differ by one coefficient. Replacing $1/6$ by $1/8$ in |
---|
| 506 | (\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. |
---|
[817] | 507 | This option is not available through a namelist parameter, since the $1/6$ coefficient |
---|
[2282] | 508 | is hard coded. Nevertheless it is quite easy to make the substitution in the |
---|
[817] | 509 | \mdl{dynadv\_ubs} module and obtain a QUICK scheme. |
---|
[707] | 510 | |
---|
[817] | 511 | Note also that in the current version of \mdl{dynadv\_ubs}, there is also the |
---|
| 512 | possibility of using a $4^{th}$ order evaluation of the advective velocity as in |
---|
| 513 | ROMS. This is an error and should be suppressed soon. |
---|
| 514 | %%% |
---|
| 515 | \gmcomment{action : this have to be done} |
---|
| 516 | %%% |
---|
[707] | 517 | |
---|
| 518 | % ================================================================ |
---|
| 519 | % Hydrostatic pressure gradient term |
---|
| 520 | % ================================================================ |
---|
[817] | 521 | \section [Hydrostatic pressure gradient (\textit{dynhpg})] |
---|
| 522 | {Hydrostatic pressure gradient (\mdl{dynhpg})} |
---|
[707] | 523 | \label{DYN_hpg} |
---|
| 524 | %------------------------------------------nam_dynhpg--------------------------------------------------- |
---|
[2282] | 525 | \namdisplay{namdyn_hpg} |
---|
[707] | 526 | %------------------------------------------------------------------------------------------------------------- |
---|
| 527 | |
---|
[817] | 528 | The key distinction between the different algorithms used for the hydrostatic |
---|
| 529 | pressure gradient is the vertical coordinate used, since HPG is a \emph{horizontal} |
---|
| 530 | pressure gradient, $i.e.$ computed along geopotential surfaces. As a result, any |
---|
| 531 | tilt of the surface of the computational levels will require a specific treatment to |
---|
| 532 | compute the hydrostatic pressure gradient. |
---|
[707] | 533 | |
---|
[817] | 534 | The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, |
---|
[2282] | 535 | $i.e.$ the density appearing in its expression is centred in time (\emph{now} $\rho$), or |
---|
[817] | 536 | a semi-implcit scheme. At the lateral boundaries either free slip, no slip or partial slip |
---|
| 537 | boundary conditions are applied. |
---|
[707] | 538 | |
---|
| 539 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 540 | % z-coordinate with full step |
---|
| 541 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 542 | \subsection [$z$-coordinate with full step (\np{ln\_dynhpg\_zco}) ] |
---|
[2282] | 543 | {$z$-coordinate with full step (\np{ln\_dynhpg\_zco}=true)} |
---|
[707] | 544 | \label{DYN_hpg_zco} |
---|
| 545 | |
---|
[817] | 546 | The hydrostatic pressure can be obtained by integrating the hydrostatic equation |
---|
| 547 | vertically from the surface. However, the pressure is large at great depth while its |
---|
| 548 | horizontal gradient is several orders of magnitude smaller. This may lead to large |
---|
| 549 | truncation errors in the pressure gradient terms. Thus, the two horizontal components |
---|
| 550 | of the hydrostatic pressure gradient are computed directly as follows: |
---|
[707] | 551 | |
---|
| 552 | for $k=km$ (surface layer, $jk=1$ in the code) |
---|
| 553 | \begin{equation} \label{Eq_dynhpg_zco_surf} |
---|
| 554 | \left\{ \begin{aligned} |
---|
| 555 | \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k=km} |
---|
| 556 | &= \frac{1}{2} g \ \left. \delta _{i+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
| 557 | \left. \delta _{j+1/2} \left[ p^h \right] \right|_{k=km} |
---|
| 558 | &= \frac{1}{2} g \ \left. \delta _{j+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
| 559 | \end{aligned} \right. |
---|
| 560 | \end{equation} |
---|
| 561 | |
---|
| 562 | for $1<k<km$ (interior layer) |
---|
| 563 | \begin{equation} \label{Eq_dynhpg_zco} |
---|
| 564 | \left\{ \begin{aligned} |
---|
| 565 | \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k} |
---|
| 566 | &= \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k-1} |
---|
| 567 | + \frac{1}{2}\;g\; \left. \delta _{i+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
| 568 | \left. \delta _{j+1/2} \left[ p^h \right] \right|_{k} |
---|
| 569 | &= \left. \delta _{j+1/2} \left[ p^h \right] \right|_{k-1} |
---|
| 570 | + \frac{1}{2}\;g\; \left. \delta _{j+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
| 571 | \end{aligned} \right. |
---|
| 572 | \end{equation} |
---|
| 573 | |
---|
[817] | 574 | Note that the $1/2$ factor in (\ref{Eq_dynhpg_zco_surf}) is adequate because of |
---|
| 575 | the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface |
---|
[2282] | 576 | level ($z=0$). Note also that in case of variable volume level (\key{vvl} defined), the |
---|
| 577 | surface pressure gradient is included in \eqref{Eq_dynhpg_zco_surf} and \eqref{Eq_dynhpg_zco} |
---|
| 578 | through the space and time variations of the vertical scale factor $e_{3w}$. |
---|
[707] | 579 | |
---|
| 580 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 581 | % z-coordinate with partial step |
---|
| 582 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 583 | \subsection [$z$-coordinate with partial step (\np{ln\_dynhpg\_zps})] |
---|
[2282] | 584 | {$z$-coordinate with partial step (\np{ln\_dynhpg\_zps}=true)} |
---|
[707] | 585 | \label{DYN_hpg_zps} |
---|
| 586 | |
---|
[817] | 587 | With partial bottom cells, tracers in horizontally adjacent cells generally live at |
---|
| 588 | different depths. Before taking horizontal gradients between these tracer points, |
---|
| 589 | a linear interpolation is used to approximate the deeper tracer as if it actually lived |
---|
| 590 | at the depth of the shallower tracer point. |
---|
[707] | 591 | |
---|
[817] | 592 | Apart from this modification, the horizontal hydrostatic pressure gradient evaluated |
---|
| 593 | in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case. |
---|
| 594 | As explained in detail in section \S\ref{TRA_zpshde}, the nonlinearity of pressure |
---|
| 595 | effects in the equation of state is such that it is better to interpolate temperature and |
---|
| 596 | salinity vertically before computing the density. Horizontal gradients of temperature |
---|
| 597 | and salinity are needed for the TRA modules, which is the reason why the horizontal |
---|
| 598 | gradients of density at the deepest model level are computed in module \mdl{zpsdhe} |
---|
| 599 | located in the TRA directory and described in \S\ref{TRA_zpshde}. |
---|
[707] | 600 | |
---|
| 601 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 602 | % s- and s-z-coordinates |
---|
| 603 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 604 | \subsection{$s$- and $z$-$s$-coordinates} |
---|
[707] | 605 | \label{DYN_hpg_sco} |
---|
| 606 | |
---|
[2282] | 607 | Pressure gradient formulations in an $s$-coordinate have been the subject of a vast |
---|
| 608 | number of papers ($e.g.$, \citet{Song1998, Shchepetkin_McWilliams_OM05}). |
---|
[3294] | 609 | A number of different pressure gradient options are coded but the ROMS-like, density Jacobian with |
---|
| 610 | cubic polynomial method is currently disabled whilst known bugs are under investigation. |
---|
[707] | 611 | |
---|
[3294] | 612 | $\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true) |
---|
[707] | 613 | \begin{equation} \label{Eq_dynhpg_sco} |
---|
| 614 | \left\{ \begin{aligned} |
---|
[817] | 615 | - \frac{1} {\rho_o \, e_{1u}} \; \delta _{i+1/2} \left[ p^h \right] |
---|
[2282] | 616 | + \frac{g\; \overline {\rho}^{i+1/2}} {\rho_o \, e_{1u}} \; \delta _{i+1/2} \left[ z_t \right] \\ |
---|
[817] | 617 | - \frac{1} {\rho_o \, e_{2v}} \; \delta _{j+1/2} \left[ p^h \right] |
---|
[2282] | 618 | + \frac{g\; \overline {\rho}^{j+1/2}} {\rho_o \, e_{2v}} \; \delta _{j+1/2} \left[ z_t \right] \\ |
---|
[707] | 619 | \end{aligned} \right. |
---|
| 620 | \end{equation} |
---|
| 621 | |
---|
[817] | 622 | Where the first term is the pressure gradient along coordinates, computed as in |
---|
| 623 | \eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of |
---|
| 624 | the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point |
---|
[3294] | 625 | ($e_{3w}$). |
---|
[707] | 626 | |
---|
[3294] | 627 | $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true) |
---|
[707] | 628 | |
---|
[2282] | 629 | $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} |
---|
[3294] | 630 | (\np{ln\_dynhpg\_djc}=true) (currently disabled; under development) |
---|
[707] | 631 | |
---|
[3294] | 632 | Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume formulation is |
---|
| 633 | activated (\key{vvl}) because in that case, even with a flat bottom, the coordinate surfaces are not |
---|
| 634 | horizontal but follow the free surface \citep{Levier2007}. The pressure jacobian scheme |
---|
| 635 | (\np{ln\_dynhpg\_prj}=true) is available as an improved option to \np{ln\_dynhpg\_sco}=true when |
---|
| 636 | \key{vvl} is active. The pressure Jacobian scheme uses a constrained cubic spline to reconstruct |
---|
| 637 | the density profile across the water column. This method maintains the monotonicity between the |
---|
| 638 | density nodes The pressure can be calculated by analytical integration of the density profile and a |
---|
| 639 | pressure Jacobian method is used to solve the horizontal pressure gradient. This method can provide |
---|
| 640 | a more accurate calculation of the horizontal pressure gradient than the standard scheme. |
---|
[707] | 641 | |
---|
| 642 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 643 | % Time-scheme |
---|
| 644 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 645 | \subsection [Time-scheme (\np{ln\_dynhpg\_imp}) ] |
---|
[2282] | 646 | {Time-scheme (\np{ln\_dynhpg\_imp}= true/false)} |
---|
[707] | 647 | \label{DYN_hpg_imp} |
---|
| 648 | |
---|
[817] | 649 | The default time differencing scheme used for the horizontal pressure gradient is |
---|
| 650 | a leapfrog scheme and therefore the density used in all discrete expressions given |
---|
| 651 | above is the \textit{now} density, computed from the \textit{now} temperature and |
---|
| 652 | salinity. In some specific cases (usually high resolution simulations over an ocean |
---|
[2282] | 653 | domain which includes weakly stratified regions) the physical phenomenon that |
---|
[817] | 654 | controls the time-step is internal gravity waves (IGWs). A semi-implicit scheme for |
---|
[2282] | 655 | doubling the stability limit associated with IGWs can be used \citep{Brown_Campana_MWR78, |
---|
[817] | 656 | Maltrud1998}. It involves the evaluation of the hydrostatic pressure gradient as an |
---|
[2282] | 657 | average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ ($i.e.$ |
---|
| 658 | \textit{before}, \textit{now} and \textit{after} time-steps), rather than at the central |
---|
[817] | 659 | time level $t$ only, as in the standard leapfrog scheme. |
---|
[707] | 660 | |
---|
[2282] | 661 | $\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}=true): |
---|
[707] | 662 | |
---|
| 663 | \begin{equation} \label{Eq_dynhpg_lf} |
---|
[2282] | 664 | \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; |
---|
| 665 | -\frac{1}{\rho _o \,e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right] |
---|
[707] | 666 | \end{equation} |
---|
| 667 | |
---|
[2282] | 668 | $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}=true): |
---|
[707] | 669 | \begin{equation} \label{Eq_dynhpg_imp} |
---|
[2282] | 670 | \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \; |
---|
| 671 | -\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] | 672 | \end{equation} |
---|
| 673 | |
---|
[817] | 674 | The semi-implicit time scheme \eqref{Eq_dynhpg_imp} is made possible without |
---|
| 675 | significant additional computation since the density can be updated to time level |
---|
[2282] | 676 | $t+\rdt$ before computing the horizontal hydrostatic pressure gradient. It can |
---|
[817] | 677 | be easily shown that the stability limit associated with the hydrostatic pressure |
---|
| 678 | gradient doubles using \eqref{Eq_dynhpg_imp} compared to that using the |
---|
| 679 | standard leapfrog scheme \eqref{Eq_dynhpg_lf}. Note that \eqref{Eq_dynhpg_imp} |
---|
| 680 | is equivalent to applying a time filter to the pressure gradient to eliminate high |
---|
| 681 | frequency IGWs. Obviously, when using \eqref{Eq_dynhpg_imp}, the doubling of |
---|
| 682 | the time-step is achievable only if no other factors control the time-step, such as |
---|
| 683 | the stability limits associated with advection or diffusion. |
---|
[707] | 684 | |
---|
[2282] | 685 | In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}=true. |
---|
[817] | 686 | In this case, we choose to apply the time filter to temperature and salinity used in |
---|
| 687 | the equation of state, instead of applying it to the hydrostatic pressure or to the |
---|
| 688 | density, so that no additional storage array has to be defined. The density used to |
---|
| 689 | compute the hydrostatic pressure gradient (whatever the formulation) is evaluated |
---|
| 690 | as follows: |
---|
[707] | 691 | \begin{equation} \label{Eq_rho_flt} |
---|
[2282] | 692 | \rho^t = \rho( \widetilde{T},\widetilde {S},z_t) |
---|
[707] | 693 | \quad \text{with} \quad |
---|
[2282] | 694 | \widetilde{X} = 1 / 4 \left( X^{t+\rdt} +2 \,X^t + X^{t-\rdt} \right) |
---|
[707] | 695 | \end{equation} |
---|
| 696 | |
---|
[817] | 697 | Note that in the semi-implicit case, it is necessary to save the filtered density, an |
---|
| 698 | extra three-dimensional field, in the restart file to restart the model with exact |
---|
[2282] | 699 | reproducibility. This option is controlled by \np{nn\_dynhpg\_rst}, a namelist parameter. |
---|
[707] | 700 | |
---|
| 701 | % ================================================================ |
---|
| 702 | % Surface Pressure Gradient |
---|
| 703 | % ================================================================ |
---|
[817] | 704 | \section [Surface pressure gradient (\textit{dynspg}) ] |
---|
| 705 | {Surface pressure gradient (\mdl{dynspg})} |
---|
[2282] | 706 | \label{DYN_spg} |
---|
[707] | 707 | %-----------------------------------------nam_dynspg---------------------------------------------------- |
---|
[2282] | 708 | \namdisplay{namdyn_spg} |
---|
[707] | 709 | %------------------------------------------------------------------------------------------------------------ |
---|
| 710 | |
---|
[2282] | 711 | $\ $\newline %force an empty line |
---|
[817] | 712 | |
---|
[2282] | 713 | %%% |
---|
| 714 | The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}). The main distinction is between the fixed volume case (linear free surface) and the variable volume case (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface}) the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case (\S\ref{PE_free_surface}). With both linear and nonlinear free surface, external gravity waves are allowed in the equations, which imposes a very small time step when an explicit time stepping is used. Two methods are proposed to allow a longer time step for the three-dimensional equations: the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}), and the split-explicit free surface described below. The extra term introduced in the filtered method is calculated implicitly, so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
| 715 | |
---|
| 716 | %%% |
---|
| 717 | |
---|
| 718 | |
---|
| 719 | The form of the surface pressure gradient term depends on how the user wants to handle |
---|
| 720 | the fast external gravity waves that are a solution of the analytical equation (\S\ref{PE_hor_pg}). |
---|
| 721 | Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx): |
---|
| 722 | an explicit formulation which requires a small time step ; |
---|
| 723 | a filtered free surface formulation which allows a larger time step by adding a filtering |
---|
| 724 | term into the momentum equation ; |
---|
| 725 | and a split-explicit free surface formulation, described below, which also allows a larger time step. |
---|
| 726 | |
---|
| 727 | The extra term introduced in the filtered method is calculated |
---|
| 728 | implicitly, so that a solver is used to compute it. As a consequence the update of the $next$ |
---|
| 729 | velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
| 730 | |
---|
| 731 | |
---|
| 732 | |
---|
[707] | 733 | %-------------------------------------------------------------------------------------------------------------- |
---|
[2282] | 734 | % Explicit free surface formulation |
---|
[707] | 735 | %-------------------------------------------------------------------------------------------------------------- |
---|
[2282] | 736 | \subsection{Explicit free surface (\key{dynspg\_exp})} |
---|
[707] | 737 | \label{DYN_spg_exp} |
---|
| 738 | |
---|
[2282] | 739 | In the explicit free surface formulation (\key{dynspg\_exp} defined), the model time step |
---|
| 740 | is chosen to be small enough to resolve the external gravity waves (typically a few tens of seconds). |
---|
| 741 | The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time), |
---|
| 742 | is thus simply given by : |
---|
[707] | 743 | \begin{equation} \label{Eq_dynspg_exp} |
---|
| 744 | \left\{ \begin{aligned} |
---|
[2282] | 745 | - \frac{1}{e_{1u}\,\rho_o} \; \delta _{i+1/2} \left[ \,\rho \,\eta\, \right] \\ |
---|
| 746 | - \frac{1}{e_{2v}\,\rho_o} \; \delta _{j+1/2} \left[ \,\rho \,\eta\, \right] |
---|
[707] | 747 | \end{aligned} \right. |
---|
| 748 | \end{equation} |
---|
| 749 | |
---|
[2282] | 750 | Note that in the non-linear free surface case ($i.e.$ \key{vvl} defined), the surface pressure |
---|
| 751 | gradient is already included in the momentum tendency through the level thickness variation |
---|
| 752 | allowed in the computation of the hydrostatic pressure gradient. Thus, nothing is done in the \mdl{dynspg\_exp} module. |
---|
[707] | 753 | |
---|
[2282] | 754 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 755 | % Split-explict free surface formulation |
---|
| 756 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 757 | \subsection{Split-Explicit free surface (\key{dynspg\_ts})} |
---|
[707] | 758 | \label{DYN_spg_ts} |
---|
| 759 | |
---|
[2282] | 760 | The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined), |
---|
| 761 | also called the time-splitting formulation, follows the one |
---|
| 762 | proposed by \citet{Griffies_Bk04}. The general idea is to solve the free surface |
---|
| 763 | equation and the associated barotropic velocity equations with a smaller time |
---|
| 764 | step than $\rdt$, the time step used for the three dimensional prognostic |
---|
[2349] | 765 | variables (Fig.~\ref{Fig_DYN_dynspg_ts}). |
---|
| 766 | The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) |
---|
[2282] | 767 | is provided through the \np{nn\_baro} namelist parameter as: |
---|
[2349] | 768 | $\rdt_e = \rdt / nn\_baro$. |
---|
[2282] | 769 | |
---|
[817] | 770 | |
---|
[707] | 771 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
[2376] | 772 | \begin{figure}[!t] \begin{center} |
---|
[998] | 773 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf} |
---|
[2376] | 774 | \caption{ \label{Fig_DYN_dynspg_ts} |
---|
| 775 | Schematic of the split-explicit time stepping scheme for the external |
---|
[994] | 776 | and internal modes. Time increases to the right. |
---|
| 777 | Internal mode time steps (which are also the model time steps) are denoted |
---|
[2282] | 778 | by $t-\rdt$, $t, t+\rdt$, and $t+2\rdt$. |
---|
[994] | 779 | The curved line represents a leap-frog time step, and the smaller time |
---|
[2282] | 780 | steps $N \rdt_e=\frac{3}{2}\rdt$ are denoted by the zig-zag line. |
---|
[1224] | 781 | The vertically integrated forcing \textbf{M}(t) computed at the model time step $t$ |
---|
[994] | 782 | represents the interaction between the external and internal motions. |
---|
[1224] | 783 | While keeping \textbf{M} and freshwater forcing field fixed, a leap-frog |
---|
| 784 | integration carries the external mode variables (surface height and vertically |
---|
[2282] | 785 | integrated velocity) from $t$ to $t+\frac{3}{2} \rdt$ using N external time |
---|
| 786 | steps of length $\rdt_e$. Time averaging the external fields over the |
---|
[1224] | 787 | $\frac{2}{3}N+1$ time steps (endpoints included) centers the vertically integrated |
---|
[2282] | 788 | velocity and the sea surface height at the model timestep $t+\rdt$. |
---|
[1224] | 789 | These averaged values are used to update \textbf{M}(t) with both the surface |
---|
[2282] | 790 | pressure gradient and the Coriolis force, therefore providing the $t+\rdt$ |
---|
[1224] | 791 | velocity. The model time stepping scheme can then be achieved by a baroclinic |
---|
[2282] | 792 | leap-frog time step that carries the surface height from $t-\rdt$ to $t+\rdt$. } |
---|
[2376] | 793 | \end{center} \end{figure} |
---|
[707] | 794 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
| 795 | |
---|
[817] | 796 | The split-explicit formulation has a damping effect on external gravity waves, |
---|
[2282] | 797 | which is weaker damping than that for the filtered free surface but still significant, as |
---|
[817] | 798 | shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave. |
---|
[707] | 799 | |
---|
[2282] | 800 | %>>>>>=============== |
---|
| 801 | \gmcomment{ %%% copy from griffies Book |
---|
[707] | 802 | |
---|
[2282] | 803 | \textbf{title: Time stepping the barotropic system } |
---|
[707] | 804 | |
---|
[2282] | 805 | Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$. Hence, |
---|
| 806 | we can update the surface height and vertically integrated velocity with a leap-frog |
---|
| 807 | scheme using the small barotropic time step $\rdt$. We have |
---|
[707] | 808 | |
---|
[2282] | 809 | \begin{equation} \label{DYN_spg_ts_eta} |
---|
| 810 | \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1}) |
---|
| 811 | = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right] |
---|
| 812 | \end{equation} |
---|
| 813 | \begin{multline} \label{DYN_spg_ts_u} |
---|
| 814 | \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}) \\ |
---|
| 815 | = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n}) |
---|
| 816 | - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right] |
---|
| 817 | \end{multline} |
---|
| 818 | \ |
---|
[707] | 819 | |
---|
[2282] | 820 | In these equations, araised (b) denotes values of surface height and vertically integrated velocity updated with the barotropic time steps. The $\tau$ time label on $\eta^{(b)}$ |
---|
| 821 | and $U^{(b)}$ denotes the baroclinic time at which the vertically integrated forcing $\textbf{M}(\tau)$ (note that this forcing includes the surface freshwater forcing), the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$, and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over a single cycle. This is also the time |
---|
| 822 | that sets the barotropic time steps via |
---|
| 823 | \begin{equation} \label{DYN_spg_ts_t} |
---|
| 824 | t_n=\tau+n\rdt |
---|
| 825 | \end{equation} |
---|
| 826 | with $n$ an integer. The density scaled surface pressure is evaluated via |
---|
| 827 | \begin{equation} \label{DYN_spg_ts_ps} |
---|
| 828 | p_s^{(b)}(\tau,t_{n}) = \begin{cases} |
---|
| 829 | g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o & \text{non-linear case} \\ |
---|
| 830 | g \;\eta_s^{(b)}(\tau,t_{n}) & \text{linear case} |
---|
| 831 | \end{cases} |
---|
| 832 | \end{equation} |
---|
| 833 | To get started, we assume the following initial conditions |
---|
| 834 | \begin{equation} \label{DYN_spg_ts_eta} |
---|
| 835 | \begin{split} |
---|
| 836 | \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)} |
---|
| 837 | \\ |
---|
| 838 | \eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} |
---|
| 839 | \end{split} |
---|
| 840 | \end{equation} |
---|
| 841 | with |
---|
| 842 | \begin{equation} \label{DYN_spg_ts_etaF} |
---|
| 843 | \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n}) |
---|
| 844 | \end{equation} |
---|
| 845 | the time averaged surface height taken from the previous barotropic cycle. Likewise, |
---|
| 846 | \begin{equation} \label{DYN_spg_ts_u} |
---|
| 847 | \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ |
---|
| 848 | \\ |
---|
| 849 | \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} |
---|
| 850 | \end{equation} |
---|
| 851 | with |
---|
| 852 | \begin{equation} \label{DYN_spg_ts_u} |
---|
| 853 | \overline{\textbf{U}^{(b)}(\tau)} |
---|
| 854 | = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n}) |
---|
| 855 | \end{equation} |
---|
| 856 | the time averaged vertically integrated transport. Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration. |
---|
[707] | 857 | |
---|
[2349] | 858 | Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ , the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at baroclinic time $\tau + \rdt \tau$ |
---|
[2282] | 859 | \begin{equation} \label{DYN_spg_ts_u} |
---|
| 860 | \textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} |
---|
| 861 | = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n}) |
---|
| 862 | \end{equation} |
---|
| 863 | The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using the following form |
---|
| 864 | |
---|
| 865 | \begin{equation} \label{DYN_spg_ts_ssh} |
---|
| 866 | \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] |
---|
| 867 | \end{equation} |
---|
| 868 | |
---|
| 869 | The use of this "big-leap-frog" scheme for the surface height ensures compatibility between the mass/volume budgets and the tracer budgets. More discussion of this point is provided in Chapter 10 (see in particular Section 10.2). |
---|
| 870 | |
---|
| 871 | In general, some form of time filter is needed to maintain integrity of the surface |
---|
| 872 | height field due to the leap-frog splitting mode in equation \ref{DYN_spg_ts_ssh}. We |
---|
| 873 | have tried various forms of such filtering, with the following method discussed in |
---|
| 874 | \cite{Griffies_al_MWR01} chosen due to its stability and reasonably good maintenance of |
---|
| 875 | tracer conservation properties (see Section ??) |
---|
| 876 | |
---|
| 877 | \begin{equation} \label{DYN_spg_ts_sshf} |
---|
| 878 | \eta^{F}(\tau-\Delta) = \overline{\eta^{(b)}(\tau)} |
---|
| 879 | \end{equation} |
---|
| 880 | Another approach tried was |
---|
| 881 | |
---|
| 882 | \begin{equation} \label{DYN_spg_ts_sshf2} |
---|
| 883 | \eta^{F}(\tau-\Delta) = \eta(\tau) |
---|
| 884 | + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt) |
---|
| 885 | + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right] |
---|
| 886 | \end{equation} |
---|
| 887 | |
---|
| 888 | which is useful since it isolates all the time filtering aspects into the term multiplied |
---|
| 889 | by $\alpha$. This isolation allows for an easy check that tracer conservation is exact when |
---|
| 890 | eliminating tracer and surface height time filtering (see Section ?? for more complete discussion). However, in the general case with a non-zero $\alpha$, the filter \ref{DYN_spg_ts_sshf} was found to be more conservative, and so is recommended. |
---|
| 891 | |
---|
| 892 | } %%end gm comment (copy of griffies book) |
---|
| 893 | |
---|
| 894 | %>>>>>=============== |
---|
| 895 | |
---|
| 896 | |
---|
[707] | 897 | %-------------------------------------------------------------------------------------------------------------- |
---|
[2282] | 898 | % Filtered free surface formulation |
---|
[707] | 899 | %-------------------------------------------------------------------------------------------------------------- |
---|
[2282] | 900 | \subsection{Filtered free surface (\key{dynspg\_flt})} |
---|
| 901 | \label{DYN_spg_fltp} |
---|
[707] | 902 | |
---|
[2282] | 903 | The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation. |
---|
[2541] | 904 | The extra term introduced in the equations (see \S\ref{PE_free_surface}) is solved implicitly. |
---|
[2282] | 905 | The elliptic solvers available in the code are documented in \S\ref{MISC}. |
---|
[707] | 906 | |
---|
[2282] | 907 | %% gm %%======>>>> given here the discrete eqs provided to the solver |
---|
[2541] | 908 | \gmcomment{ %%% copy from chap-model basics |
---|
| 909 | \begin{equation} \label{Eq_spg_flt} |
---|
| 910 | \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}} |
---|
| 911 | - g \nabla \left( \tilde{\rho} \ \eta \right) |
---|
| 912 | - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) |
---|
| 913 | \end{equation} |
---|
| 914 | where $T_c$, is a parameter with dimensions of time which characterizes the force, |
---|
| 915 | $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ |
---|
| 916 | represents the collected contributions of the Coriolis, hydrostatic pressure gradient, |
---|
| 917 | non-linear and viscous terms in \eqref{Eq_PE_dyn}. |
---|
| 918 | } %end gmcomment |
---|
[707] | 919 | |
---|
[2282] | 920 | Note that in the linear free surface formulation (\key{vvl} not defined), the ocean depth |
---|
| 921 | is time-independent and so is the matrix to be inverted. It is computed once and for all and applies to all ocean time steps. |
---|
[707] | 922 | |
---|
| 923 | % ================================================================ |
---|
| 924 | % Lateral diffusion term |
---|
| 925 | % ================================================================ |
---|
[817] | 926 | \section [Lateral diffusion term (\textit{dynldf})] |
---|
| 927 | {Lateral diffusion term (\mdl{dynldf})} |
---|
[707] | 928 | \label{DYN_ldf} |
---|
| 929 | %------------------------------------------nam_dynldf---------------------------------------------------- |
---|
[2282] | 930 | \namdisplay{namdyn_ldf} |
---|
[707] | 931 | %------------------------------------------------------------------------------------------------------------- |
---|
| 932 | |
---|
[2282] | 933 | The options available for lateral diffusion are to use either laplacian |
---|
[817] | 934 | (rotated or not) or biharmonic operators. The coefficients may be constant |
---|
| 935 | or spatially variable; the description of the coefficients is found in the chapter |
---|
[2282] | 936 | on lateral physics (Chap.\ref{LDF}). The lateral diffusion of momentum is |
---|
| 937 | evaluated using a forward scheme, $i.e.$ the velocity appearing in its expression |
---|
[817] | 938 | is the \textit{before} velocity in time, except for the pure vertical component |
---|
| 939 | that appears when a tensor of rotation is used. This latter term is solved |
---|
| 940 | implicitly together with the vertical diffusion term (see \S\ref{DOM_nxt}) |
---|
| 941 | |
---|
[707] | 942 | At the lateral boundaries either free slip, no slip or partial slip boundary |
---|
[817] | 943 | conditions are applied according to the user's choice (see Chap.\ref{LBC}). |
---|
[707] | 944 | |
---|
| 945 | % ================================================================ |
---|
[817] | 946 | \subsection [Iso-level laplacian operator (\np{ln\_dynldf\_lap}) ] |
---|
[2282] | 947 | {Iso-level laplacian operator (\np{ln\_dynldf\_lap}=true)} |
---|
[707] | 948 | \label{DYN_ldf_lap} |
---|
| 949 | |
---|
| 950 | For lateral iso-level diffusion, the discrete operator is: |
---|
| 951 | \begin{equation} \label{Eq_dynldf_lap} |
---|
| 952 | \left\{ \begin{aligned} |
---|
| 953 | D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta _{i+1/2} \left[ {A_T^{lm} |
---|
| 954 | \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta _j \left[ |
---|
| 955 | {A_f^{lm} \;e_{3f} \zeta } \right] \\ |
---|
| 956 | \\ |
---|
| 957 | D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta _{j+1/2} \left[ {A_T^{lm} |
---|
| 958 | \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta _i \left[ |
---|
| 959 | {A_f^{lm} \;e_{3f} \zeta } \right] \\ |
---|
| 960 | \end{aligned} \right. |
---|
| 961 | \end{equation} |
---|
| 962 | |
---|
| 963 | As explained in \S\ref{PE_ldf}, this formulation (as the gradient of a divergence |
---|
| 964 | and curl of the vorticity) preserves symmetry and ensures a complete |
---|
[2282] | 965 | separation between the vorticity and divergence parts of the momentum diffusion. |
---|
[707] | 966 | |
---|
| 967 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 968 | % Rotated laplacian operator |
---|
| 969 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 970 | \subsection [Rotated laplacian operator (\np{ln\_dynldf\_iso}) ] |
---|
[2282] | 971 | {Rotated laplacian operator (\np{ln\_dynldf\_iso}=true)} |
---|
[707] | 972 | \label{DYN_ldf_iso} |
---|
| 973 | |
---|
[2282] | 974 | A rotation of the lateral momentum diffusion operator is needed in several cases: |
---|
| 975 | for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}=true) and for |
---|
| 976 | either iso-neutral (\np{ln\_dynldf\_iso}=true) or geopotential |
---|
| 977 | (\np{ln\_dynldf\_hor}=true) diffusion in the $s$-coordinate. In the partial step |
---|
| 978 | case, coordinates are horizontal except at the deepest level and no |
---|
| 979 | rotation is performed when \np{ln\_dynldf\_hor}=true. The diffusion operator |
---|
[707] | 980 | is defined simply as the divergence of down gradient momentum fluxes on each |
---|
[817] | 981 | momentum component. It must be emphasized that this formulation ignores |
---|
[707] | 982 | constraints on the stress tensor such as symmetry. The resulting discrete |
---|
| 983 | representation is: |
---|
| 984 | \begin{equation} \label{Eq_dyn_ldf_iso} |
---|
| 985 | \begin{split} |
---|
| 986 | D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ |
---|
| 987 | & \left\{\quad {\delta _{i+1/2} \left[ {A_T^{lm} \left( |
---|
[2282] | 988 | {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta _{i}[u] |
---|
| 989 | -e_{2t} \; r_{1t} \,\overline{\overline {\delta _{k+1/2}[u]}}^{\,i,\,k}} |
---|
[707] | 990 | \right)} \right]} \right. |
---|
| 991 | \\ |
---|
| 992 | & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} |
---|
| 993 | }\,\delta _{j+1/2} [u] - e_{1f}\, r_{2f} |
---|
| 994 | \,\overline{\overline {\delta _{k+1/2} [u]}} ^{\,j+1/2,\,k}} |
---|
| 995 | \right)} \right] |
---|
| 996 | \\ |
---|
| 997 | &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline |
---|
| 998 | {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } |
---|
| 999 | \right.} \right. |
---|
| 1000 | \\ |
---|
| 1001 | & \ \qquad \qquad \qquad \quad\ |
---|
| 1002 | - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} |
---|
| 1003 | \\ |
---|
| 1004 | & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ |
---|
| 1005 | +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} |
---|
| 1006 | \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} |
---|
| 1007 | \\ |
---|
| 1008 | \\ |
---|
| 1009 | D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\ |
---|
| 1010 | & \left\{\quad {\delta _{i+1/2} \left[ {A_f^{lm} \left( |
---|
| 1011 | {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta _{i+1/2}[v] |
---|
| 1012 | -e_{2f} \; r_{1f} \,\overline{\overline {\delta _{k+1/2}[v]}}^{\,i+1/2,\,k}} |
---|
| 1013 | \right)} \right]} \right. |
---|
| 1014 | \\ |
---|
[2282] | 1015 | & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} |
---|
| 1016 | }\,\delta _{j} [v] - e_{1t}\, r_{2t} |
---|
[707] | 1017 | \,\overline{\overline {\delta _{k+1/2} [v]}} ^{\,j,\,k}} |
---|
| 1018 | \right)} \right] |
---|
| 1019 | \\ |
---|
| 1020 | & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline |
---|
| 1021 | {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. |
---|
| 1022 | \\ |
---|
| 1023 | & \ \qquad \qquad \qquad \quad\ |
---|
| 1024 | - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} |
---|
| 1025 | \\ |
---|
| 1026 | & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ |
---|
| 1027 | +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} |
---|
| 1028 | \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} |
---|
| 1029 | \end{split} |
---|
| 1030 | \end{equation} |
---|
[817] | 1031 | where $r_1$ and $r_2$ are the slopes between the surface along which the |
---|
[2282] | 1032 | diffusion operator acts and the surface of computation ($z$- or $s$-surfaces). |
---|
[817] | 1033 | The way these slopes are evaluated is given in the lateral physics chapter |
---|
| 1034 | (Chap.\ref{LDF}). |
---|
[707] | 1035 | |
---|
| 1036 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1037 | % Iso-level bilaplacian operator |
---|
| 1038 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 1039 | \subsection [Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap})] |
---|
[2282] | 1040 | {Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap}=true)} |
---|
[707] | 1041 | \label{DYN_ldf_bilap} |
---|
| 1042 | |
---|
| 1043 | The lateral fourth order operator formulation on momentum is obtained by |
---|
[817] | 1044 | applying \eqref{Eq_dynldf_lap} twice. It requires an additional assumption on |
---|
| 1045 | boundary conditions: the first derivative term normal to the coast depends on |
---|
| 1046 | the free or no-slip lateral boundary conditions chosen, while the third |
---|
| 1047 | derivative terms normal to the coast are set to zero (see Chap.\ref{LBC}). |
---|
| 1048 | %%% |
---|
| 1049 | \gmcomment{add a remark on the the change in the position of the coefficient} |
---|
| 1050 | %%% |
---|
[707] | 1051 | |
---|
| 1052 | % ================================================================ |
---|
| 1053 | % Vertical diffusion term |
---|
| 1054 | % ================================================================ |
---|
[817] | 1055 | \section [Vertical diffusion term (\mdl{dynzdf})] |
---|
| 1056 | {Vertical diffusion term (\mdl{dynzdf})} |
---|
[707] | 1057 | \label{DYN_zdf} |
---|
| 1058 | %----------------------------------------------namzdf------------------------------------------------------ |
---|
| 1059 | \namdisplay{namzdf} |
---|
| 1060 | %------------------------------------------------------------------------------------------------------------- |
---|
| 1061 | |
---|
[817] | 1062 | The large vertical diffusion coefficient found in the surface mixed layer together |
---|
| 1063 | with high vertical resolution implies that in the case of explicit time stepping there |
---|
| 1064 | would be too restrictive a constraint on the time step. Two time stepping schemes |
---|
| 1065 | can be used for the vertical diffusion term : $(a)$ a forward time differencing |
---|
[2282] | 1066 | scheme (\np{ln\_zdfexp}=true) using a time splitting technique |
---|
| 1067 | (\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme |
---|
| 1068 | (\np{ln\_zdfexp}=false) (see \S\ref{DOM_nxt}). Note that namelist variables |
---|
| 1069 | \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. |
---|
[707] | 1070 | |
---|
| 1071 | The formulation of the vertical subgrid scale physics is the same whatever |
---|
[817] | 1072 | the vertical coordinate is. The vertical diffusion operators given by |
---|
[707] | 1073 | \eqref{Eq_PE_zdf} take the following semi-discrete space form: |
---|
| 1074 | \begin{equation} \label{Eq_dynzdf} |
---|
| 1075 | \left\{ \begin{aligned} |
---|
| 1076 | D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta _k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } |
---|
| 1077 | \ \delta _{k+1/2} [\,u\,] \right] \\ |
---|
| 1078 | \\ |
---|
| 1079 | D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta _k \left[ \frac{A_{vw}^{vm} }{e_{3vw} } |
---|
| 1080 | \ \delta _{k+1/2} [\,v\,] \right] |
---|
| 1081 | \end{aligned} \right. |
---|
| 1082 | \end{equation} |
---|
[817] | 1083 | where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and |
---|
| 1084 | diffusivity coefficients. The way these coefficients are evaluated |
---|
[707] | 1085 | depends on the vertical physics used (see \S\ref{ZDF}). |
---|
| 1086 | |
---|
[2282] | 1087 | The surface boundary condition on momentum is the stress exerted by |
---|
[707] | 1088 | the wind. At the surface, the momentum fluxes are prescribed as the boundary |
---|
| 1089 | condition on the vertical turbulent momentum fluxes, |
---|
| 1090 | \begin{equation} \label{Eq_dynzdf_sbc} |
---|
| 1091 | \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} |
---|
| 1092 | = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v } |
---|
| 1093 | \end{equation} |
---|
[817] | 1094 | where $\left( \tau _u ,\tau _v \right)$ are the two components of the wind stress |
---|
| 1095 | vector in the (\textbf{i},\textbf{j}) coordinate system. The high mixing coefficients |
---|
| 1096 | in the surface mixed layer ensure that the surface wind stress is distributed in |
---|
[707] | 1097 | the vertical over the mixed layer depth. If the vertical mixing coefficient |
---|
| 1098 | is small (when no mixed layer scheme is used) the surface stress enters only |
---|
| 1099 | the top model level, as a body force. The surface wind stress is calculated |
---|
[817] | 1100 | in the surface module routines (SBC, see Chap.\ref{SBC}) |
---|
[707] | 1101 | |
---|
[817] | 1102 | The turbulent flux of momentum at the bottom of the ocean is specified through |
---|
[1224] | 1103 | a bottom friction parameterisation (see \S\ref{ZDF_bfr}) |
---|
[707] | 1104 | |
---|
| 1105 | % ================================================================ |
---|
| 1106 | % External Forcing |
---|
| 1107 | % ================================================================ |
---|
| 1108 | \section{External Forcings} |
---|
| 1109 | \label{DYN_forcing} |
---|
| 1110 | |
---|
| 1111 | Besides the surface and bottom stresses (see the above section) which are |
---|
| 1112 | introduced as boundary conditions on the vertical mixing, two other forcings |
---|
| 1113 | enter the dynamical equations. |
---|
| 1114 | |
---|
[2282] | 1115 | One is the effect of atmospheric pressure on the ocean dynamics. |
---|
| 1116 | Another forcing term is the tidal potential. |
---|
| 1117 | Both of which will be introduced into the reference version soon. |
---|
[707] | 1118 | |
---|
[2285] | 1119 | \gmcomment{atmospheric pressure is there!!!! include its description } |
---|
| 1120 | |
---|
[707] | 1121 | % ================================================================ |
---|
| 1122 | % Time evolution term |
---|
| 1123 | % ================================================================ |
---|
[817] | 1124 | \section [Time evolution term (\textit{dynnxt})] |
---|
| 1125 | {Time evolution term (\mdl{dynnxt})} |
---|
[707] | 1126 | \label{DYN_nxt} |
---|
| 1127 | |
---|
| 1128 | %----------------------------------------------namdom---------------------------------------------------- |
---|
| 1129 | \namdisplay{namdom} |
---|
| 1130 | %------------------------------------------------------------------------------------------------------------- |
---|
| 1131 | |
---|
[817] | 1132 | The general framework for dynamics time stepping is a leap-frog scheme, |
---|
| 1133 | $i.e.$ a three level centred time scheme associated with an Asselin time filter |
---|
[2282] | 1134 | (cf. Chap.\ref{STP}). The scheme is applied to the velocity, except when using |
---|
| 1135 | the flux form of momentum advection (cf. \S\ref{DYN_adv_cor_flux}) in the variable |
---|
| 1136 | volume case (\key{vvl} defined), where it has to be applied to the thickness |
---|
| 1137 | weighted velocity (see \S\ref{Apdx_A_momentum}) |
---|
[707] | 1138 | |
---|
[2282] | 1139 | $\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}=true ; \key{vvl} not defined): |
---|
| 1140 | \begin{equation} \label{Eq_dynnxt_vec} |
---|
| 1141 | \left\{ \begin{aligned} |
---|
| 1142 | &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt \ \text{RHS}_u^t \\ |
---|
| 1143 | &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right] |
---|
| 1144 | \end{aligned} \right. |
---|
[707] | 1145 | \end{equation} |
---|
| 1146 | |
---|
[2282] | 1147 | $\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}=false ; \key{vvl} defined): |
---|
| 1148 | \begin{equation} \label{Eq_dynnxt_flux} |
---|
[707] | 1149 | \left\{ \begin{aligned} |
---|
[2282] | 1150 | &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t \\ |
---|
| 1151 | &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t |
---|
| 1152 | +\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] |
---|
[707] | 1153 | \end{aligned} \right. |
---|
| 1154 | \end{equation} |
---|
[2282] | 1155 | where RHS is the right hand side of the momentum equation, the subscript $f$ |
---|
| 1156 | denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is |
---|
| 1157 | initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = $10^{-3}$. |
---|
| 1158 | In both cases, the modified Asselin filter is not applied since perfect conservation |
---|
| 1159 | is not an issue for the momentum equations. |
---|
[707] | 1160 | |
---|
[2282] | 1161 | Note that with the filtered free surface, the update of the \textit{after} velocities |
---|
| 1162 | is done in the \mdl{dynsp\_flt} module, and only array swapping |
---|
| 1163 | and Asselin filtering is done in \mdl{dynnxt}. |
---|
[707] | 1164 | |
---|
[1224] | 1165 | % ================================================================ |
---|
[3294] | 1166 | % Neptune effect |
---|
| 1167 | % ================================================================ |
---|
| 1168 | \section [Neptune effect (\textit{dynnept})] |
---|
| 1169 | {Neptune effect (\mdl{dynnept})} |
---|
| 1170 | \label{DYN_nept} |
---|
| 1171 | |
---|
| 1172 | The "Neptune effect" (thus named in \citep{HollowayOM86}) is a |
---|
| 1173 | parameterisation of the potentially large effect of topographic form stress |
---|
| 1174 | (caused by eddies) in driving the ocean circulation. Originally developed for |
---|
| 1175 | low-resolution models, in which it was applied via a Laplacian (second-order) |
---|
| 1176 | diffusion-like term in the momentum equation, it can also be applied in eddy |
---|
| 1177 | permitting or resolving models, in which a more scale-selective bilaplacian |
---|
| 1178 | (fourth-order) implementation is preferred. This mechanism has a |
---|
| 1179 | significant effect on boundary currents (including undercurrents), and the |
---|
| 1180 | upwelling of deep water near continental shelves. |
---|
| 1181 | |
---|
| 1182 | The theoretical basis for the method can be found in |
---|
| 1183 | \citep{HollowayJPO92}, including the explanation of why form stress is not |
---|
| 1184 | necessarily a drag force, but may actually drive the flow. |
---|
| 1185 | \citep{HollowayJPO94} demonstrate the effects of the parameterisation in |
---|
| 1186 | the GFDL-MOM model, at a horizontal resolution of about 1.8 degrees. |
---|
| 1187 | \citep{HollowayOM08} demonstrate the biharmonic version of the |
---|
| 1188 | parameterisation in a global run of the POP model, with an average horizontal |
---|
| 1189 | grid spacing of about 32km. |
---|
| 1190 | |
---|
| 1191 | The NEMO implementation is a simplified form of that supplied by |
---|
| 1192 | Greg Holloway, the testing of which was described in \citep{HollowayJGR09}. |
---|
| 1193 | The major simplification is that a time invariant Neptune velocity |
---|
| 1194 | field is assumed. This is computed only once, during start-up, and |
---|
| 1195 | made available to the rest of the code via a module. Vertical |
---|
| 1196 | diffusive terms are also ignored, and the model topography itself |
---|
| 1197 | is used, rather than a separate topographic dataset as in |
---|
| 1198 | \citep{HollowayOM08}. This implementation is only in the iso-level |
---|
| 1199 | formulation, as is the case anyway for the bilaplacian operator. |
---|
| 1200 | |
---|
| 1201 | The velocity field is derived from a transport stream function given by: |
---|
| 1202 | |
---|
| 1203 | \begin{equation} \label{Eq_dynnept_sf} |
---|
| 1204 | \psi = -fL^2H |
---|
| 1205 | \end{equation} |
---|
| 1206 | |
---|
| 1207 | where $L$ is a latitude-dependant length scale given by: |
---|
| 1208 | |
---|
| 1209 | \begin{equation} \label{Eq_dynnept_ls} |
---|
| 1210 | L = l_1 + (l_2 -l_1)\left ( {1 + \cos 2\phi \over 2 } \right ) |
---|
| 1211 | \end{equation} |
---|
| 1212 | |
---|
| 1213 | where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively. |
---|
| 1214 | Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as: |
---|
| 1215 | |
---|
| 1216 | \begin{equation} \label{Eq_dynnept_vel} |
---|
| 1217 | u^* = -{1\over H} {\partial \psi \over \partial y}\ \ \ ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x} |
---|
| 1218 | \end{equation} |
---|
| 1219 | |
---|
| 1220 | \smallskip |
---|
| 1221 | %----------------------------------------------namdom---------------------------------------------------- |
---|
| 1222 | \namdisplay{namdyn_nept} |
---|
| 1223 | %-------------------------------------------------------------------------------------------------------- |
---|
| 1224 | \smallskip |
---|
| 1225 | |
---|
| 1226 | The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false). |
---|
| 1227 | \np{ln\_smooth\_neptvel} controls whether a scale-selective smoothing is applied |
---|
| 1228 | to the Neptune effect flow field (default=false) (this smoothing method is as |
---|
| 1229 | used by Holloway). \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and |
---|
| 1230 | polar values respectively of the length-scale parameter $L$ used in determining |
---|
| 1231 | the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}. |
---|
| 1232 | Values at intermediate latitudes are given by a cosine fit, mimicking the |
---|
| 1233 | variation of the deformation radius with latitude. The default values of 12km |
---|
| 1234 | and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse |
---|
| 1235 | resolution model. The finer resolution study of \citep{HollowayOM08} increased |
---|
| 1236 | the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the |
---|
| 1237 | stream function for a given topography. |
---|
| 1238 | |
---|
| 1239 | The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities |
---|
| 1240 | in shallow water, and \citep{HollowayOM08} add an offset to the depth in the |
---|
| 1241 | denominator to control this problem. In this implementation we offer instead (at |
---|
| 1242 | the suggestion of G. Madec) the option of ramping down the Neptune flow field to |
---|
| 1243 | zero over a finite depth range. The switch \np{ln\_neptramp} activates this |
---|
| 1244 | option (default=false), in which case velocities at depths greater than |
---|
| 1245 | \np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a |
---|
| 1246 | depth of \np{rn\_htrmin} (and shallower). |
---|
| 1247 | |
---|
| 1248 | % ================================================================ |
---|