[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 | |
---|
[817] | 13 | Using the representation described in Chap.\ref{DOM}, several semi-discrete |
---|
| 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. |
---|
| 17 | One must be aware that all the quantities are masked fields and that each time a |
---|
| 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*} |
---|
| 26 | |
---|
[817] | 27 | NXT stands for next, referring to the time-stepping. The first group of terms on |
---|
| 28 | the rhs of the momentum equations corresponds to the Coriolis and advection |
---|
| 29 | terms that are decomposed into a vorticity part (VOR), a kinetic energy part (KEG) |
---|
| 30 | and, a vertical advection part (ZAD) in the vector invariant formulation or a Coriolis |
---|
| 31 | and advection part(COR+ADV) in the flux formulation. The terms following these |
---|
| 32 | are the pressure gradient contributions (HPG, Hydrostatic Pressure Gradient, |
---|
| 33 | and SPG, Surface Pressure Gradient); and contributions from lateral diffusion |
---|
| 34 | (LDF) and vertical diffusion (ZDF), which are added to the rhs in the \mdl{dynldf} |
---|
| 35 | and \mdl{dynzdf} modules. The vertical diffusion term includes the surface and |
---|
| 36 | bottom stresses. The external forcings and parameterisations require complex |
---|
| 37 | inputs (surface wind stress calculation using bulk formulae, estimation of mixing |
---|
| 38 | coefficients) that are carried out in modules SBC, LDF and ZDF and are described |
---|
| 39 | in Chapters \ref{SBC}, \ref{LDF} and \ref{ZDF}, respectively. |
---|
[707] | 40 | |
---|
[817] | 41 | In the present chapter we also describe the diagnostic equations used to compute |
---|
| 42 | the horizontal divergence and curl of the velocities (\emph{divcur} module) as well |
---|
| 43 | as the vertical velocity (\emph{wzvmod} module). |
---|
[707] | 44 | |
---|
[817] | 45 | The different options available to the user are managed by namelist variables. |
---|
| 46 | For equation term \textit{ttt}, the logical namelist variables are \textit{ln\_dynttt\_xxx}, |
---|
| 47 | where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. |
---|
| 48 | If a CPP key is used for this term its name is \textbf{key\_ttt}. The corresponding |
---|
| 49 | code can be found in the \textit{dynttt\_xxx} module in the DYN directory, and it is |
---|
| 50 | usually computed in the \textit{dyn\_ttt\_xxx} subroutine. |
---|
[707] | 51 | |
---|
[817] | 52 | The user has the option of extracting each tendency term of both the rhs of the |
---|
| 53 | 3D momentum equation (\key{trddyn} defined) for output, as described in |
---|
| 54 | Chap.\ref{MISC}. Furthermore, the tendency terms associated to the 2D |
---|
| 55 | barotropic vorticity balance (\key{trdvor} defined) can be derived on-line from the |
---|
| 56 | 3D terms. |
---|
| 57 | %%% |
---|
[996] | 58 | \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does |
---|
| 59 | MISC correspond to "extracting tendency terms" or "vorticity balance"?} |
---|
[707] | 60 | |
---|
[1224] | 61 | $\ $\newline % force a new ligne |
---|
[707] | 62 | % ================================================================ |
---|
| 63 | % Coriolis and Advection terms: vector invariant form |
---|
| 64 | % ================================================================ |
---|
| 65 | \section{Coriolis and Advection: vector invariant form} |
---|
| 66 | \label{DYN_adv_cor_vect} |
---|
| 67 | %-----------------------------------------nam_dynadv---------------------------------------------------- |
---|
| 68 | \namdisplay{nam_dynadv} |
---|
| 69 | %------------------------------------------------------------------------------------------------------------- |
---|
| 70 | |
---|
[817] | 71 | The vector invariant form of the momentum equations is the one most |
---|
| 72 | often used in applications of \NEMO ocean model. The flux form option |
---|
[1224] | 73 | (see next section) has been introduced since version $2$. |
---|
| 74 | Coriolis and momentum advection terms are evaluated using a leapfrog |
---|
| 75 | scheme, $i.e.$ the velocity appearing in these expressions is centred in |
---|
| 76 | time (\textit{now} velocity). |
---|
[817] | 77 | At the lateral boundaries either free slip, no slip or partial slip boundary |
---|
| 78 | conditions are applied following Chap.\ref{LBC}. |
---|
[707] | 79 | |
---|
| 80 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 81 | % Vorticity term |
---|
| 82 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 83 | \subsection [Vorticity term (\textit{dynvor}) ] |
---|
| 84 | {Vorticity term (\mdl{dynvor})} |
---|
[707] | 85 | \label{DYN_vor} |
---|
| 86 | %------------------------------------------nam_dynvor---------------------------------------------------- |
---|
| 87 | \namdisplay{nam_dynvor} |
---|
| 88 | %------------------------------------------------------------------------------------------------------------- |
---|
| 89 | |
---|
[817] | 90 | Different discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=.true.) are |
---|
[1224] | 91 | available: conserving potential enstrophy of horizontally non-divergent flow ; |
---|
| 92 | conserving horizontal kinetic energy ; or conserving potential enstrophy for the |
---|
| 93 | relative vorticity term and horizontal kinetic energy for the planetary vorticity |
---|
| 94 | term (see Appendix~\ref{Apdx_C}). The vorticity terms are given below for the |
---|
| 95 | general case, but note that in the full step $z$-coordinate (\key{zco} is defined), |
---|
| 96 | $e_{3u} =e_{3v} =e_{3f}$ so that the vertical scale factors disappear. They are |
---|
| 97 | all computed in dedicated routines that can be found in the \mdl{dynvor} module. |
---|
[707] | 98 | |
---|
| 99 | %------------------------------------------------------------- |
---|
| 100 | % enstrophy conserving scheme |
---|
| 101 | %------------------------------------------------------------- |
---|
[817] | 102 | \subsubsection{Enstrophy conserving scheme (\np{ln\_dynvor\_ens}=.true.)} |
---|
[707] | 103 | \label{DYN_vor_ens} |
---|
| 104 | |
---|
[817] | 105 | In the enstrophy conserving case (ENS scheme), the discrete formulation of the |
---|
| 106 | vorticity term provides a global conservation of the enstrophy |
---|
| 107 | ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent |
---|
| 108 | flow ($i.e.$ $\chi=0$), but does not conserve the total kinetic energy. It is given by: |
---|
[707] | 109 | \begin{equation} \label{Eq_dynvor_ens} |
---|
| 110 | \left\{ |
---|
| 111 | \begin{aligned} |
---|
[1224] | 112 | {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} |
---|
| 113 | & {\overline{\overline {\left( {e_{1v} e_{3v} v} \right)}} }^{\,i, j+1/2} \\ |
---|
| 114 | {- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} |
---|
| 115 | & {\overline{\overline {\left( {e_{2u} e_{3u} u} \right)}} }^{\,i+1/2, j} |
---|
[707] | 116 | \end{aligned} |
---|
| 117 | \right. |
---|
| 118 | \end{equation} |
---|
| 119 | |
---|
| 120 | %------------------------------------------------------------- |
---|
| 121 | % energy conserving scheme |
---|
| 122 | %------------------------------------------------------------- |
---|
[817] | 123 | \subsubsection{Energy conserving scheme (\np{ln\_dynvor\_ene}=.true.)} |
---|
[707] | 124 | \label{DYN_vor_ene} |
---|
| 125 | |
---|
[817] | 126 | The kinetic energy conserving scheme (ENE scheme) conserves the global |
---|
| 127 | kinetic energy but not the global enstrophy. It is given by: |
---|
[707] | 128 | \begin{equation} \label{Eq_dynvor_ene} |
---|
[1224] | 129 | \left\{ \begin{aligned} |
---|
| 130 | {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
| 131 | \; \overline {\left( {e_{1v} e_{3v} v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
| 132 | {- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right) |
---|
| 133 | \; \overline {\left( {e_{2u} e_{3u} u} \right)} ^{\,j+1/2}} }^{\,i} } |
---|
| 134 | \end{aligned} \right. |
---|
[707] | 135 | \end{equation} |
---|
| 136 | |
---|
| 137 | %------------------------------------------------------------- |
---|
| 138 | % mix energy/enstrophy conserving scheme |
---|
| 139 | %------------------------------------------------------------- |
---|
[817] | 140 | \subsubsection{Mixed energy/enstrophy conserving scheme (\np{ln\_dynvor\_mix}=.true.) } |
---|
[707] | 141 | \label{DYN_vor_mix} |
---|
| 142 | |
---|
[817] | 143 | The mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the |
---|
| 144 | two previous schemes is used. It consists of the ENS scheme (\ref{Eq_dynvor_ens}) |
---|
| 145 | to the relative vorticity term, and of the ENE scheme (\ref{Eq_dynvor_ene}) applied |
---|
| 146 | to the planetary vorticity term. |
---|
[707] | 147 | \begin{equation} \label{Eq_dynvor_mix} |
---|
| 148 | \left\{ { |
---|
| 149 | \begin{aligned} |
---|
[994] | 150 | {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} |
---|
[707] | 151 | \; {\overline{\overline {\left( {e_{1v} \; e_{3v} \ v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} } |
---|
| 152 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
| 153 | \;\overline {\left( {e_{1v} \; e_{3v} \ v} \right)} ^{\,i+1/2}} }^{\,j} } \\ |
---|
[994] | 154 | {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j |
---|
[707] | 155 | \; {\overline{\overline {\left( {e_{2u} \; e_{3u} \ u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} } |
---|
| 156 | \; {\overline {\left( {\frac{f}{e_{3f} }} \right) |
---|
| 157 | \;\overline {\left( {e_{2u}\; e_{3u} \ u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill |
---|
| 158 | \end{aligned} |
---|
| 159 | } \right. |
---|
| 160 | \end{equation} |
---|
| 161 | |
---|
| 162 | %------------------------------------------------------------- |
---|
| 163 | % energy and enstrophy conserving scheme |
---|
| 164 | %------------------------------------------------------------- |
---|
[817] | 165 | \subsubsection{Energy and enstrophy conserving scheme (\np{ln\_dynvor\_een}=.true.) } |
---|
[707] | 166 | \label{DYN_vor_een} |
---|
| 167 | |
---|
[817] | 168 | In the energy and enstrophy conserving scheme (EEN scheme), the vorticity term |
---|
| 169 | is evaluated using the vorticity advection scheme of \citet{Arakawa1990}. |
---|
| 170 | This scheme conserves both total energy and potential enstrophy in the limit of |
---|
| 171 | horizontally nondivergent flow ($i.e. \ \chi=0$). While EEN is more complicated |
---|
| 172 | than ENS or ENE and does not conserve potential enstrophy and total energy in |
---|
| 173 | general flow, it tolerates arbitrarily thin layers. This feature is essential for |
---|
| 174 | $z$-coordinate with partial step. |
---|
| 175 | %%% |
---|
| 176 | \gmcomment{gm : it actually conserve kinetic energy ! show that in appendix C } |
---|
| 177 | %%% |
---|
[707] | 178 | |
---|
[817] | 179 | The \citet{Arakawa1990} vorticity advection scheme for a single layer is modified |
---|
| 180 | for spherical coordinates as described by \citet{Arakawa1981} to obtain the EEN |
---|
| 181 | scheme. The potential vorticity, defined at an $f$-point, is: |
---|
[707] | 182 | \begin{equation} \label{Eq_pot_vor} |
---|
| 183 | q_f = \frac{\zeta +f} {e_{3f} } |
---|
| 184 | \end{equation} |
---|
[817] | 185 | where the relative vorticity is defined by (\ref{Eq_divcur_cur}), the Coriolis parameter |
---|
| 186 | is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: |
---|
[707] | 187 | \begin{equation} \label{Eq_een_e3f} |
---|
| 188 | e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} |
---|
| 189 | \end{equation} |
---|
| 190 | |
---|
| 191 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 192 | \begin{figure}[!ht] \label{Fig_DYN_een_triad} |
---|
| 193 | \begin{center} |
---|
[998] | 194 | \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_DYN_een_triad.pdf} |
---|
[817] | 195 | \caption{Triads used in the energy and enstrophy conserving scheme (een) for |
---|
| 196 | $u$-component (upper panel) and $v$-component (lower panel).} |
---|
[707] | 197 | \end{center} |
---|
| 198 | \end{figure} |
---|
| 199 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 200 | |
---|
[817] | 201 | Note that a key point in \eqref{Eq_een_e3f} is that the averaging in \textbf{i}- and |
---|
| 202 | \textbf{j}- directions uses the masked vertical scale factor but is always divided by |
---|
| 203 | $4$, not by the sum of the mask at $T$-point. This preserves the continuity of |
---|
| 204 | $e_{3f}$ when one or more of the neighbouring $e_{3T}$ tends to zero and |
---|
| 205 | extends by continuity the value of $e_{3f}$ in the land areas. |
---|
| 206 | %%% |
---|
| 207 | \gmcomment{this has to be further investigate in case of several step topography} |
---|
| 208 | %%% |
---|
[707] | 209 | |
---|
| 210 | The vorticity terms are represented as: |
---|
| 211 | \begin{equation} \label{Eq_dynvor_een} |
---|
| 212 | \left\{ { |
---|
| 213 | \begin{aligned} |
---|
| 214 | +q\,e_3 \, v &\equiv +\frac{1}{e_{1u} } \left[ |
---|
| 215 | {{\begin{array}{*{20}c} |
---|
| 216 | {\,\ \ a_{j+1/2}^{i } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i+1/2} } |
---|
| 217 | {\,+\,b_{j+1/2}^{i } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i-1/2} } \\ |
---|
| 218 | \\ |
---|
| 219 | { +\,c_{j-1/2}^{i } \left( {e_{1v} e_{3v} \ v} \right)_{j }^{i+1/2} } |
---|
| 220 | {\,+\,d_{j+1/2}^{i } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i+1/2} } \\ |
---|
| 221 | \end{array} }} \right] \\ |
---|
| 222 | \\ |
---|
| 223 | -q\,e_3 \,u &\equiv -\frac{1}{e_{2v} } \left[ |
---|
| 224 | {{\begin{array}{*{20}c} |
---|
[994] | 225 | {\,\ \ a_{j-1/2}^{i } \left( {e_{2u} e_{3u} \ u} \right)_{j+1}^{i+1/2} } |
---|
| 226 | {\,+\,b_{j-1/2}^{i+1} \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i+1} } \\ |
---|
[707] | 227 | \\ |
---|
[994] | 228 | { +\,c_{j+1/2}^{i+1} \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i+1} } |
---|
| 229 | {\,+\,d_{j+1/2}^{i } \left( {e_{2u} e_{3u} \ u} \right)_{j+1/2}^{i } } \\ |
---|
[707] | 230 | \end{array} }} \right] |
---|
| 231 | \end{aligned} |
---|
| 232 | } \right. |
---|
| 233 | \end{equation} |
---|
[1224] | 234 | where $a$, $b$, $c$ and $d$ are the following triad combinations of the |
---|
| 235 | neighbouring potential vorticities (Fig.~\ref{Fig_DYN_een_triad}): |
---|
[707] | 236 | \begin{equation} \label{Eq_een_triads} |
---|
| 237 | \left\{ |
---|
| 238 | \begin{aligned} |
---|
| 239 | a_{\,j+1/2}^i & = \frac{1} {12} \left( q_{j+1/2}^{i+1} + q_{j+1 /2}^i + q_{j-1/2}^i \right) \\ |
---|
| 240 | \\ |
---|
| 241 | b_{\,j+1/2}^i & = \frac{1} {12} \left( q_{j+1/2}^{i-1} +q_{j+1/2}^i +q_{j-1/2}^i \right) \\ |
---|
| 242 | \\ |
---|
| 243 | c_{\,j+1/2}^i & = \frac{1} {12} \left( q_{j-1/2}^{i-1} +q_{j+1/2}^i +q_{j-1/2}^i \right) \\ |
---|
| 244 | \\ |
---|
| 245 | d_{\,j+1/2}^i & = \frac{1} {12} \left( q_{j-1/2}^{i+1} +q_{j+1/2}^i +q_{j-1/2}^i \right) \\ |
---|
| 246 | \end{aligned} |
---|
| 247 | \right. |
---|
| 248 | \end{equation} |
---|
| 249 | |
---|
| 250 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 251 | % Kinetic Energy Gradient term |
---|
| 252 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 253 | \subsection [Kinetic Energy Gradient term (\textit{dynkeg})] |
---|
| 254 | {Kinetic Energy Gradient term (\mdl{dynkeg})} |
---|
[707] | 255 | \label{DYN_keg} |
---|
| 256 | |
---|
[817] | 257 | As demonstarted in Appendix~\ref{Apdx_C}, there is a single discrete formulation |
---|
| 258 | of the kinetic energy gradient term that, together with the formulation chosen for |
---|
| 259 | the vertical advection (see below), conserves the total kinetic energy: |
---|
[707] | 260 | \begin{equation} \label{Eq_dynkeg} |
---|
| 261 | \left\{ \begin{aligned} |
---|
| 262 | -\frac{1}{2 \; e_{1u} } |
---|
| 263 | & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] \\ |
---|
| 264 | -\frac{1}{2 \; e_{2v} } |
---|
| 265 | & \ \delta _{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right] |
---|
| 266 | \end{aligned} \right. |
---|
| 267 | \end{equation} |
---|
| 268 | |
---|
| 269 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 270 | % Vertical advection term |
---|
| 271 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 272 | \subsection [Vertical advection term (\textit{dynzad}) ] |
---|
| 273 | {Vertical advection term (\mdl{dynzad}) } |
---|
[707] | 274 | \label{DYN_zad} |
---|
| 275 | |
---|
[817] | 276 | The discrete formulation of the vertical advection, together with the formulation |
---|
| 277 | chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic |
---|
| 278 | energy. Indeed, the change of KE due to the vertical advection is exactly |
---|
| 279 | balanced by the change of KE due to the gradient of KE (see Appendix~\ref{Apdx_C}). |
---|
[707] | 280 | \begin{equation} \label{Eq_dynzad} |
---|
| 281 | \left\{ \begin{aligned} |
---|
| 282 | -\frac{1} { e_{1u}\,e_{2u}\,e_{3u} } & |
---|
| 283 | \ {\overline {\overline{ e_{1T}\,e_{2T}\,w } ^{\,i+1/2} \;\delta _{k+1/2} \left[ u \right] }^{\,k} } \\ |
---|
| 284 | -\frac{1} { e_{1v}\,e_{2v}\,e_{3v} } & |
---|
| 285 | \ {\overline {\overline{ e_{1T}\,e_{2T}\,w } ^{\,j+1/2} \;\delta _{k+1/2} \left[ u \right] }^{\,k} } |
---|
| 286 | \end{aligned} \right. |
---|
| 287 | \end{equation} |
---|
| 288 | |
---|
| 289 | % ================================================================ |
---|
| 290 | % Coriolis and Advection : flux form |
---|
| 291 | % ================================================================ |
---|
| 292 | \section{Coriolis and Advection: flux form} |
---|
| 293 | \label{DYN_adv_cor_flux} |
---|
| 294 | %------------------------------------------nam_dynadv---------------------------------------------------- |
---|
| 295 | \namdisplay{nam_dynadv} |
---|
| 296 | %------------------------------------------------------------------------------------------------------------- |
---|
| 297 | |
---|
[817] | 298 | In the flux form (as in the vector invariant form), the Coriolis and momentum |
---|
| 299 | advection terms are evaluated using a leapfrog scheme, $i.e.$ the velocity |
---|
| 300 | appearing in their expressions is centred in time (\textit{now} velocity). At the |
---|
| 301 | lateral boundaries either free slip, no slip or partial slip boundary conditions |
---|
| 302 | are applied following Chap.\ref{LBC}. |
---|
[707] | 303 | |
---|
| 304 | |
---|
| 305 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 306 | % Coriolis plus curvature metric terms |
---|
| 307 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 308 | \subsection [Coriolis plus curvature metric terms (\textit{dynvor}) ] |
---|
| 309 | {Coriolis plus curvature metric terms (\mdl{dynvor}) } |
---|
[707] | 310 | \label{DYN_cor_flux} |
---|
| 311 | |
---|
[817] | 312 | In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis |
---|
| 313 | parameter has been modified to account for the "metric" term. This altered |
---|
| 314 | Coriolis parameter is thus discretised at $f$-points. It is given by: |
---|
[707] | 315 | \begin{multline} \label{Eq_dyncor_metric} |
---|
| 316 | f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i} - u\frac{\partial e_1 }{\partial j}} \right) \\ |
---|
| 317 | \equiv f + \frac{1}{e_{1f} e_{2f} } |
---|
| 318 | \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right] |
---|
| 319 | - \overline u ^{j+1/2}\delta _{j+1/2} \left[ {e_{1u} } \right] } \ \right) |
---|
| 320 | \end{multline} |
---|
| 321 | |
---|
[817] | 322 | Any of the (\ref{Eq_dynvor_ens}), (\ref{Eq_dynvor_ene}) and (\ref{Eq_dynvor_een}) |
---|
| 323 | schemes can be used to compute the product of the Coriolis parameter and the |
---|
| 324 | vorticity. However, the energy-conserving scheme (\ref{Eq_dynvor_een}) has |
---|
| 325 | exclusively been used to date. This term is evaluated using a leapfrog scheme, |
---|
| 326 | $i.e.$ the velocity is centred in time (\textit{now} velocity). |
---|
[707] | 327 | |
---|
| 328 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 329 | % Flux form Advection term |
---|
| 330 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 331 | \subsection [Flux form Advection term (\textit{dynadv}) ] |
---|
| 332 | {Flux form Advection term (\mdl{dynadv}) } |
---|
[707] | 333 | \label{DYN_adv_flux} |
---|
| 334 | |
---|
| 335 | The discrete expression of the advection term is given by : |
---|
| 336 | \begin{equation} \label{Eq_dynadv} |
---|
| 337 | \left\{ |
---|
| 338 | \begin{aligned} |
---|
| 339 | \frac{1}{e_{1u}\,e_{2u}\,e_{3u}} |
---|
| 340 | \left( \delta _{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\ u }^{i } \ u_T \right] |
---|
| 341 | + \delta _{j } \left[ \overline{e_{1u}\,e_{3u}\ v }^{i+1/2} \ u_F \right] \right. \ \; \\ |
---|
| 342 | \left. + \delta _{k } \left[ \overline{e_{1w}\,e_{2w} w}^{i+1/2} \ u_{uw} \right] \right) \\ |
---|
| 343 | \\ |
---|
| 344 | \frac{1}{e_{1v}\,e_{2v}\,e_{3v}} |
---|
| 345 | \left( \delta _{i } \left[ \overline{e_{2u}\,e_{3u } \ u }^{j+1/2} \ v_F \right] |
---|
| 346 | + \delta _{j+1/2} \left[ \overline{e_{1u}\,e_{3u } \ v }^{i } \ v_T \right] \right. \ \, \\ |
---|
| 347 | \left. + \delta _{k } \left[ \overline{e_{1w}\,e_{2w} \ w}^{j+1/2} \ v_{vw} \right] \right) \\ |
---|
| 348 | \end{aligned} |
---|
| 349 | \right. |
---|
| 350 | \end{equation} |
---|
| 351 | |
---|
[817] | 352 | Two advection schemes are available: a $2^{nd}$ order centered finite |
---|
| 353 | difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS. |
---|
| 354 | The latter is described in \citet{Sacha2005}. The schemes are selected using |
---|
| 355 | the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. In flux |
---|
| 356 | form, the schemes differ by the choice of a space and time interpolation to |
---|
| 357 | define the value of $u$ and $v$ at the centre of each face of $u$- and $v$-cells, |
---|
| 358 | $i.e.$ at the $T$-, $f$-, and $uw$-points for $u$ and at the $f$-, $T$- and |
---|
| 359 | $vw$-points for $v$. |
---|
[707] | 360 | |
---|
| 361 | %------------------------------------------------------------- |
---|
| 362 | % 2nd order centred scheme |
---|
| 363 | %------------------------------------------------------------- |
---|
[817] | 364 | \subsubsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_dynadv\_cen2}=.true.)} |
---|
[707] | 365 | \label{DYN_adv_cen2} |
---|
| 366 | |
---|
| 367 | In the centered $2^{nd}$ order formulation, the velocity is evaluated as the |
---|
| 368 | mean of the two neighbouring points : |
---|
| 369 | \begin{equation} \label{Eq_dynadv_cen2} |
---|
| 370 | \left\{ \begin{aligned} |
---|
| 371 | u_T^{cen2} &=\overline u^{i } \quad & |
---|
[817] | 372 | u_F^{cen2} &=\overline u^{j+1/2} \quad & |
---|
[707] | 373 | u_{uw}^{cen2} &=\overline u^{k+1/2} \\ |
---|
| 374 | v_F^{cen2} &=\overline v ^{i+1/2} \quad & |
---|
| 375 | v_F^{cen2} &=\overline v^j \quad & |
---|
| 376 | v_{vw}^{cen2} &=\overline v ^{k+1/2} \\ |
---|
| 377 | \end{aligned} \right. |
---|
| 378 | \end{equation} |
---|
| 379 | |
---|
[817] | 380 | The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive |
---|
[1224] | 381 | ($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be |
---|
| 382 | used in conjunction with an explicit diffusion operator to produce a sensible solution. |
---|
| 383 | The associated time-stepping is performed using a leapfrog scheme in conjunction |
---|
| 384 | with an Asselin time-filter, so $u$ and $v$ are the \emph{now} velocities. |
---|
[707] | 385 | |
---|
| 386 | %------------------------------------------------------------- |
---|
| 387 | % UBS scheme |
---|
| 388 | %------------------------------------------------------------- |
---|
[817] | 389 | \subsubsection{Upstream Biased Scheme (UBS) (\np{ln\_dynadv\_ubs}=.true.)} |
---|
[707] | 390 | \label{DYN_adv_ubs} |
---|
| 391 | |
---|
| 392 | The UBS advection scheme is an upstream biased third order scheme based on |
---|
| 393 | an upstream-biased parabolic interpolation. For example, the evaluation of |
---|
| 394 | $u_T^{ubs} $ is done as follows: |
---|
| 395 | \begin{equation} \label{Eq_dynadv_ubs} |
---|
| 396 | u_T^{ubs} =\overline u ^i-\;\frac{1}{6} \begin{cases} |
---|
| 397 | u"_{i-1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i \geqslant 0$ } \\ |
---|
| 398 | u"_{i+1/2}& \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i < 0$ } |
---|
| 399 | \end{cases} |
---|
| 400 | \end{equation} |
---|
[817] | 401 | where $u"_{i+1/2} =\delta _{i+1/2} \left[ {\delta _i \left[ u \right]} \right]$. This results |
---|
| 402 | in a dissipatively dominant ($i.e.$ hyper-diffusive) truncation error \citep{Sacha2005}. |
---|
| 403 | The overall performance of the advection scheme is similar to that reported in |
---|
| 404 | \citet{Farrow1995}. It is a relatively good compromise between accuracy and |
---|
| 405 | smoothness. It is not a \emph{positive} scheme, meaning that false extrema are |
---|
| 406 | permitted. But the amplitudes of the false extrema are significantly reduced over |
---|
| 407 | those in the centred second order method. |
---|
[707] | 408 | |
---|
[817] | 409 | The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ |
---|
| 410 | order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and |
---|
| 411 | $u_{vw}^{ubs}$ in \eqref{Eq_dynadv_cen2} are used. UBS is diffusive and is |
---|
| 412 | associated with vertical mixing of momentum. \gmcomment{ gm pursue the |
---|
| 413 | sentence:Since vertical mixing of momentum is a source term of the TKE equation... } |
---|
[707] | 414 | |
---|
[817] | 415 | For stability reasons, the first term in (\ref{Eq_dynadv_ubs}), which corresponds |
---|
| 416 | to a second order centred scheme, is evaluated using the \textit{now} velocity |
---|
| 417 | (centred in time), while the second term, which is the diffusive part of the scheme, |
---|
| 418 | is evaluated using the \textit{before} velocity (forward in time). This is discussed |
---|
| 419 | by \citet{Webb1998} in the context of the Quick advection scheme. |
---|
[707] | 420 | |
---|
[817] | 421 | Note that the UBS and Quadratic Upstream Interpolation for Convective Kinematics |
---|
| 422 | (QUICK) schemes only differ by one coefficient. Substituting $1/6$ with $1/8$ in |
---|
| 423 | (\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}. |
---|
| 424 | This option is not available through a namelist parameter, since the $1/6$ coefficient |
---|
| 425 | is hard coded. Nevertheless it is quite easy to make the substitution in |
---|
| 426 | \mdl{dynadv\_ubs} module and obtain a QUICK scheme. |
---|
[707] | 427 | |
---|
[817] | 428 | Note also that in the current version of \mdl{dynadv\_ubs}, there is also the |
---|
| 429 | possibility of using a $4^{th}$ order evaluation of the advective velocity as in |
---|
| 430 | ROMS. This is an error and should be suppressed soon. |
---|
| 431 | %%% |
---|
| 432 | \gmcomment{action : this have to be done} |
---|
| 433 | %%% |
---|
[707] | 434 | |
---|
| 435 | % ================================================================ |
---|
| 436 | % Hydrostatic pressure gradient term |
---|
| 437 | % ================================================================ |
---|
[817] | 438 | \section [Hydrostatic pressure gradient (\textit{dynhpg})] |
---|
| 439 | {Hydrostatic pressure gradient (\mdl{dynhpg})} |
---|
[707] | 440 | \label{DYN_hpg} |
---|
| 441 | %------------------------------------------nam_dynhpg--------------------------------------------------- |
---|
| 442 | \namdisplay{nam_dynhpg} |
---|
| 443 | \namdisplay{namflg} |
---|
| 444 | %------------------------------------------------------------------------------------------------------------- |
---|
[817] | 445 | %%% |
---|
| 446 | \gmcomment{Suppress the namflg namelist and incorporate it in the namhpg namelist} |
---|
| 447 | %%% |
---|
[707] | 448 | |
---|
[817] | 449 | The key distinction between the different algorithms used for the hydrostatic |
---|
| 450 | pressure gradient is the vertical coordinate used, since HPG is a \emph{horizontal} |
---|
| 451 | pressure gradient, $i.e.$ computed along geopotential surfaces. As a result, any |
---|
| 452 | tilt of the surface of the computational levels will require a specific treatment to |
---|
| 453 | compute the hydrostatic pressure gradient. |
---|
[707] | 454 | |
---|
[817] | 455 | The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, |
---|
| 456 | $i.e.$ the density appearing in its expression is centred in time (\emph{now} rho), or |
---|
| 457 | a semi-implcit scheme. At the lateral boundaries either free slip, no slip or partial slip |
---|
| 458 | boundary conditions are applied. |
---|
[707] | 459 | |
---|
| 460 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 461 | % z-coordinate with full step |
---|
| 462 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 463 | \subsection [$z$-coordinate with full step (\np{ln\_dynhpg\_zco}) ] |
---|
| 464 | {$z$-coordinate with full step (\np{ln\_dynhpg\_zco}=.true.)} |
---|
[707] | 465 | \label{DYN_hpg_zco} |
---|
| 466 | |
---|
[817] | 467 | The hydrostatic pressure can be obtained by integrating the hydrostatic equation |
---|
| 468 | vertically from the surface. However, the pressure is large at great depth while its |
---|
| 469 | horizontal gradient is several orders of magnitude smaller. This may lead to large |
---|
| 470 | truncation errors in the pressure gradient terms. Thus, the two horizontal components |
---|
| 471 | of the hydrostatic pressure gradient are computed directly as follows: |
---|
[707] | 472 | |
---|
| 473 | for $k=km$ (surface layer, $jk=1$ in the code) |
---|
| 474 | \begin{equation} \label{Eq_dynhpg_zco_surf} |
---|
| 475 | \left\{ \begin{aligned} |
---|
| 476 | \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k=km} |
---|
| 477 | &= \frac{1}{2} g \ \left. \delta _{i+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
| 478 | \left. \delta _{j+1/2} \left[ p^h \right] \right|_{k=km} |
---|
| 479 | &= \frac{1}{2} g \ \left. \delta _{j+1/2} \left[ e_{3w} \ \rho \right] \right|_{k=km} \\ |
---|
| 480 | \end{aligned} \right. |
---|
| 481 | \end{equation} |
---|
| 482 | |
---|
| 483 | for $1<k<km$ (interior layer) |
---|
| 484 | \begin{equation} \label{Eq_dynhpg_zco} |
---|
| 485 | \left\{ \begin{aligned} |
---|
| 486 | \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k} |
---|
| 487 | &= \left. \delta _{i+1/2} \left[ p^h \right] \right|_{k-1} |
---|
| 488 | + \frac{1}{2}\;g\; \left. \delta _{i+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
| 489 | \left. \delta _{j+1/2} \left[ p^h \right] \right|_{k} |
---|
| 490 | &= \left. \delta _{j+1/2} \left[ p^h \right] \right|_{k-1} |
---|
| 491 | + \frac{1}{2}\;g\; \left. \delta _{j+1/2} \left[ e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k} \\ |
---|
| 492 | \end{aligned} \right. |
---|
| 493 | \end{equation} |
---|
| 494 | |
---|
[817] | 495 | Note that the $1/2$ factor in (\ref{Eq_dynhpg_zco_surf}) is adequate because of |
---|
| 496 | the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface |
---|
| 497 | level ($z=0)$. |
---|
[707] | 498 | |
---|
| 499 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 500 | % z-coordinate with partial step |
---|
| 501 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 502 | \subsection [$z$-coordinate with partial step (\np{ln\_dynhpg\_zps})] |
---|
| 503 | {$z$-coordinate with partial step (\np{ln\_dynhpg\_zps}=.true.)} |
---|
[707] | 504 | \label{DYN_hpg_zps} |
---|
| 505 | |
---|
[817] | 506 | With partial bottom cells, tracers in horizontally adjacent cells generally live at |
---|
| 507 | different depths. Before taking horizontal gradients between these tracer points, |
---|
| 508 | a linear interpolation is used to approximate the deeper tracer as if it actually lived |
---|
| 509 | at the depth of the shallower tracer point. |
---|
[707] | 510 | |
---|
[817] | 511 | Apart from this modification, the horizontal hydrostatic pressure gradient evaluated |
---|
| 512 | in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case. |
---|
| 513 | As explained in detail in section \S\ref{TRA_zpshde}, the nonlinearity of pressure |
---|
| 514 | effects in the equation of state is such that it is better to interpolate temperature and |
---|
| 515 | salinity vertically before computing the density. Horizontal gradients of temperature |
---|
| 516 | and salinity are needed for the TRA modules, which is the reason why the horizontal |
---|
| 517 | gradients of density at the deepest model level are computed in module \mdl{zpsdhe} |
---|
| 518 | located in the TRA directory and described in \S\ref{TRA_zpshde}. |
---|
[707] | 519 | |
---|
| 520 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 521 | % s- and s-z-coordinates |
---|
| 522 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 523 | \subsection{$s$- and $z$-$s$-coordinates} |
---|
[707] | 524 | \label{DYN_hpg_sco} |
---|
| 525 | |
---|
[817] | 526 | Pressure gradient formulations in $s$-coordinate have been the subject of a vast |
---|
| 527 | literature ($e.g.$, \citet{Song1998, Sacha2003}). A number of different pressure |
---|
| 528 | gradient options are coded, but they are not yet fully documented or tested. |
---|
[707] | 529 | |
---|
[817] | 530 | $\bullet$ Traditional coding (see for example \citet{Madec1996}: (\np{ln\_dynhpg\_sco}=.true., |
---|
| 531 | \np{ln\_dynhpg\_hel}=.true.) |
---|
[707] | 532 | \begin{equation} \label{Eq_dynhpg_sco} |
---|
| 533 | \left\{ \begin{aligned} |
---|
[817] | 534 | - \frac{1} {\rho_o \, e_{1u}} \; \delta _{i+1/2} \left[ p^h \right] |
---|
[707] | 535 | + \frac{g\; \overline {\rho}^{i+1/2}} {\rho_o \, e_{1u}} \; \delta _{i+1/2} \left[ z_T \right] \\ |
---|
[817] | 536 | - \frac{1} {\rho_o \, e_{2v}} \; \delta _{j+1/2} \left[ p^h \right] |
---|
[707] | 537 | + \frac{g\; \overline {\rho}^{j+1/2}} {\rho_o \, e_{2v}} \; \delta _{j+1/2} \left[ z_T \right] \\ |
---|
| 538 | \end{aligned} \right. |
---|
| 539 | \end{equation} |
---|
| 540 | |
---|
[817] | 541 | Where the first term is the pressure gradient along coordinates, computed as in |
---|
| 542 | \eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of |
---|
| 543 | the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point |
---|
| 544 | ($e_{3w}$). The version \np{ln\_dynhpg\_hel}=.true. has been added by Aike |
---|
| 545 | Beckmann and involves a redefinition of the relative position of $T$-points relative |
---|
| 546 | to $w$-points. |
---|
[707] | 547 | |
---|
[817] | 548 | $\bullet$ Weighted density Jacobian (WDJ) \citep{Song1998} (\np{ln\_dynhpg\_wdj}=.true.) |
---|
[707] | 549 | |
---|
[817] | 550 | $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Sacha2003} |
---|
| 551 | (\np{ln\_dynhpg\_djc}=.true.) |
---|
[707] | 552 | |
---|
[817] | 553 | $\bullet$ Rotated axes scheme (rot) \citep{Thiem2006} (\np{ln\_dynhpg\_rot}=.true.) |
---|
[707] | 554 | |
---|
[817] | 555 | Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume |
---|
| 556 | formulation is activated (\key{vvl}) because in that case, even with a flat bottom, |
---|
| 557 | the coordinate surfaces are not horizontal but follow the free surface |
---|
| 558 | \citep{Levier2007}. The other pressure gradient options are not yet available. |
---|
[707] | 559 | |
---|
| 560 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 561 | % Time-scheme |
---|
| 562 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 563 | \subsection [Time-scheme (\np{ln\_dynhpg\_imp}) ] |
---|
| 564 | {Time-scheme (\np{ln\_dynhpg\_imp}=.true./.false.)} |
---|
[707] | 565 | \label{DYN_hpg_imp} |
---|
| 566 | |
---|
[817] | 567 | The default time differencing scheme used for the horizontal pressure gradient is |
---|
| 568 | a leapfrog scheme and therefore the density used in all discrete expressions given |
---|
| 569 | above is the \textit{now} density, computed from the \textit{now} temperature and |
---|
| 570 | salinity. In some specific cases (usually high resolution simulations over an ocean |
---|
| 571 | domain which includes weakly stratified regions) the physical phenomenum that |
---|
| 572 | controls the time-step is internal gravity waves (IGWs). A semi-implicit scheme for |
---|
| 573 | doubling the stability limit associated with IGWs can be used \citep{Brown1978, |
---|
| 574 | Maltrud1998}. It involves the evaluation of the hydrostatic pressure gradient as an |
---|
| 575 | average over the three time levels $t-\Delta t$, $t$, and $t+\Delta t$ ($i.e.$ |
---|
| 576 | \textit{before}, \textit{now} and \textit{after} time-steps), rather than at central |
---|
| 577 | time level $t$ only, as in the standard leapfrog scheme. |
---|
[707] | 578 | |
---|
[817] | 579 | $\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}=.true.): |
---|
[707] | 580 | |
---|
| 581 | \begin{equation} \label{Eq_dynhpg_lf} |
---|
| 582 | \frac{u^{t+\Delta t}-u^{t-\Delta t}}{2\Delta t} |
---|
| 583 | =\;\cdots \;-\frac{1}{\rho _o \,e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right] |
---|
| 584 | \end{equation} |
---|
| 585 | |
---|
[817] | 586 | $\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}=.true.): |
---|
[707] | 587 | \begin{equation} \label{Eq_dynhpg_imp} |
---|
| 588 | \frac{u^{t+\Delta t}-u^{t-\Delta t}}{2\Delta t} |
---|
| 589 | =\;\cdots \;-\frac{1}{\rho _o \,e_{1u} } \delta _{i+1/2} \left[ \frac{ p_h^{t+\Delta t} +2p_h^t |
---|
| 590 | +p_h^{t-\Delta t} } { 4 } \right] |
---|
| 591 | \end{equation} |
---|
| 592 | |
---|
[817] | 593 | The semi-implicit time scheme \eqref{Eq_dynhpg_imp} is made possible without |
---|
| 594 | significant additional computation since the density can be updated to time level |
---|
| 595 | $t+\Delta t$ before computing the horizontal hydrostatic pressure gradient. It can |
---|
| 596 | be easily shown that the stability limit associated with the hydrostatic pressure |
---|
| 597 | gradient doubles using \eqref{Eq_dynhpg_imp} compared to that using the |
---|
| 598 | standard leapfrog scheme \eqref{Eq_dynhpg_lf}. Note that \eqref{Eq_dynhpg_imp} |
---|
| 599 | is equivalent to applying a time filter to the pressure gradient to eliminate high |
---|
| 600 | frequency IGWs. Obviously, when using \eqref{Eq_dynhpg_imp}, the doubling of |
---|
| 601 | the time-step is achievable only if no other factors control the time-step, such as |
---|
| 602 | the stability limits associated with advection or diffusion. |
---|
[707] | 603 | |
---|
[817] | 604 | In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}=.true.. |
---|
| 605 | In this case, we choose to apply the time filter to temperature and salinity used in |
---|
| 606 | the equation of state, instead of applying it to the hydrostatic pressure or to the |
---|
| 607 | density, so that no additional storage array has to be defined. The density used to |
---|
| 608 | compute the hydrostatic pressure gradient (whatever the formulation) is evaluated |
---|
| 609 | as follows: |
---|
[707] | 610 | \begin{equation} \label{Eq_rho_flt} |
---|
| 611 | \rho^t = \rho( \widetilde{T},\widetilde {S},z_T) |
---|
| 612 | \quad \text{with} \quad |
---|
[817] | 613 | \widetilde{\,\cdot\,} = \frac{ \,\cdot\,^{t+\Delta t} +2 \,\,\cdot\,^t + \,\cdot\,^{t-\Delta t} } {4} |
---|
[707] | 614 | \end{equation} |
---|
[817] | 615 | \gmcomment{STEVEN: bullets look odd in this, could use X} |
---|
[707] | 616 | |
---|
[817] | 617 | Note that in the semi-implicit case, it is necessary to save the filtered density, an |
---|
| 618 | extra three-dimensional field, in the restart file to restart the model with exact |
---|
| 619 | reproducibility. This option is controlled by the namelist parameter |
---|
| 620 | \np{nn\_dynhpg\_rst}=.true.. |
---|
[707] | 621 | |
---|
| 622 | % ================================================================ |
---|
| 623 | % Surface Pressure Gradient |
---|
| 624 | % ================================================================ |
---|
[817] | 625 | \section [Surface pressure gradient (\textit{dynspg}) ] |
---|
| 626 | {Surface pressure gradient (\mdl{dynspg})} |
---|
[707] | 627 | \label{DYN_hpg_spg} |
---|
| 628 | %-----------------------------------------nam_dynspg---------------------------------------------------- |
---|
| 629 | \namdisplay{nam_dynspg} |
---|
| 630 | %------------------------------------------------------------------------------------------------------------ |
---|
| 631 | |
---|
[817] | 632 | The form of the surface pressure gradient term is dependent on the representation |
---|
| 633 | of the free surface (\S\ref{PE_hor_pg}). The main distinction is between the fixed |
---|
| 634 | volume case (linear free surface or rigid lid) and the variable volume case |
---|
| 635 | (nonlinear free surface, \key{vvl} is defined). In the linear free surface case |
---|
| 636 | (\S\ref{PE_free_surface}) and the rigid lid case (\S\ref{PE_rigid_lid}), the vertical |
---|
| 637 | scale factors $e_{3}$ are fixed in time, whilst in the nonlinear case |
---|
| 638 | (\S\ref{PE_free_surface}) they are time-dependent. With both linear and nonlinear |
---|
| 639 | free surface, external gravity waves are allowed in the equations, which imposes |
---|
| 640 | a very small time step when an explicit time stepping is used. Two methods are |
---|
| 641 | proposed to allow a longer time step for the three-dimensional equations: the |
---|
| 642 | filtered free surface method, which involves a modification of the continuous |
---|
| 643 | equations (see \eqref{Eq_PE_flt}), and the split-explicit free surface method |
---|
| 644 | described below. The extra term introduced in the filtered method is calculated |
---|
| 645 | implicitly, so that the update of the $next$ velocities is done in module |
---|
| 646 | \mdl{dynspg\_flt} and not in \mdl{dynnxt}. |
---|
| 647 | |
---|
[707] | 648 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 649 | % Linear free surface formulation |
---|
| 650 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 651 | \subsection{Linear free surface formulation (\key{exp} or \textbf{\_ts} or \textbf{\_flt})} |
---|
[707] | 652 | \label{DYN_spg_linear} |
---|
| 653 | |
---|
[817] | 654 | In the linear free surface formulation, the sea surface height is assumed to be |
---|
| 655 | small compared to the thickness of the ocean levels, so that $(a)$ the time |
---|
| 656 | evolution of the sea surface height becomes a linear equation, and $(b)$ the |
---|
| 657 | thickness of the ocean levels is assumed to be constant with time. |
---|
| 658 | As mentioned in (\S\ref{PE_free_surface}) the linearization affects the |
---|
| 659 | conservation of tracers. |
---|
[707] | 660 | |
---|
| 661 | %------------------------------------------------------------- |
---|
| 662 | % Explicit |
---|
| 663 | %------------------------------------------------------------- |
---|
| 664 | \subsubsection{Explicit (\key{dynspg\_exp})} |
---|
| 665 | \label{DYN_spg_exp} |
---|
| 666 | |
---|
[817] | 667 | In the explicit free surface formulation, the model time step is chosen to be |
---|
| 668 | small enough to describe the external gravity waves (typically a few tens of |
---|
| 669 | seconds). The sea surface height is given by : |
---|
[707] | 670 | \begin{equation} \label{Eq_dynspg_ssh} |
---|
| 671 | \frac{\partial \eta }{\partial t}\equiv \frac{\text{EMP}}{\rho _w }+\frac{1}{e_{1T} |
---|
| 672 | e_{2T} }\sum\limits_k {\left( {\delta _i \left[ {e_{2u} e_{3u} u} |
---|
| 673 | \right]+\delta _j \left[ {e_{1v} e_{3v} v} \right]} \right)} |
---|
| 674 | \end{equation} |
---|
[994] | 675 | where EMP is the surface freshwater budget, expressed in Kg/m$^2$/s |
---|
| 676 | (which is equal to mm/s), and $\rho _w$=1,000~Kg/m$^3$ is the volumic |
---|
| 677 | mass of pure water. If river runoff is expressed as a surface freshwater flux |
---|
| 678 | (see \S\ref{SBC}) then EMP can be written as the evaporation minus |
---|
| 679 | precipitation, minus the river runoff. The sea-surface height is evaluated |
---|
| 680 | using a leapfrog scheme in combination with an Asselin time filter, $i.e.$ |
---|
| 681 | the velocity appearing in \eqref{Eq_dynspg_ssh} is centred in time |
---|
| 682 | (\textit{now} velocity). |
---|
[707] | 683 | |
---|
[817] | 684 | The surface pressure gradient, also evaluated using a leap-frog scheme, is |
---|
| 685 | then simply given by : |
---|
[707] | 686 | \begin{equation} \label{Eq_dynspg_exp} |
---|
| 687 | \left\{ \begin{aligned} |
---|
[994] | 688 | - \frac{1}{e_{1u}} \; \delta _{i+1/2} \left[ \,\eta\, \right] \\ |
---|
| 689 | - \frac{1}{e_{2v}} \; \delta _{j+1/2} \left[ \,\eta\, \right] |
---|
[707] | 690 | \end{aligned} \right. |
---|
| 691 | \end{equation} |
---|
| 692 | |
---|
[817] | 693 | Consistent with the linearization, a factor of $\left. \rho \right|_{k=1} / \rho _o$ |
---|
| 694 | is omitted in \eqref{Eq_dynspg_exp}. |
---|
[707] | 695 | |
---|
| 696 | %------------------------------------------------------------- |
---|
| 697 | % Split-explicit time-stepping |
---|
| 698 | %------------------------------------------------------------- |
---|
| 699 | \subsubsection{Split-explicit time-stepping (\key{dynspg\_ts})} |
---|
| 700 | \label{DYN_spg_ts} |
---|
| 701 | %--------------------------------------------namdom---------------------------------------------------- |
---|
| 702 | \namdisplay{namdom} |
---|
| 703 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 704 | |
---|
[817] | 705 | The split-explicit free surface formulation used in \NEMO follows the one |
---|
| 706 | proposed by \citet{Griffies2004}. The general idea is to solve the free surface |
---|
| 707 | equation with a small time step \np{rdtbt}, while the three dimensional |
---|
[994] | 708 | prognostic variables are solved with a longer time step that is a multiple of |
---|
| 709 | \np{rdtbt} (Fig.\ref {Fig_DYN_dynspg_ts}). |
---|
[817] | 710 | |
---|
[707] | 711 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
| 712 | \begin{figure}[!t] \label{Fig_DYN_dynspg_ts} |
---|
| 713 | \begin{center} |
---|
[998] | 714 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf} |
---|
[994] | 715 | \caption{Schematic of the split-explicit time stepping scheme for the external |
---|
| 716 | and internal modes. Time increases to the right. |
---|
| 717 | Internal mode time steps (which are also the model time steps) are denoted |
---|
| 718 | by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$. |
---|
| 719 | The curved line represents a leap-frog time step, and the smaller time |
---|
[1224] | 720 | steps $N \Delta t_e=\frac{3}{2}\Delta t$ are denoted by the zig-zag line. |
---|
| 721 | The vertically integrated forcing \textbf{M}(t) computed at the model time step $t$ |
---|
[994] | 722 | represents the interaction between the external and internal motions. |
---|
[1224] | 723 | While keeping \textbf{M} and freshwater forcing field fixed, a leap-frog |
---|
| 724 | integration carries the external mode variables (surface height and vertically |
---|
| 725 | integrated velocity) from $t$ to $t+\frac{3}{2} \Delta t$ using N external time |
---|
| 726 | steps of length $\Delta t_e$. Time averaging the external fields over the |
---|
| 727 | $\frac{2}{3}N+1$ time steps (endpoints included) centers the vertically integrated |
---|
| 728 | velocity and the sea surface height at the model timestep $t+\Delta t$. |
---|
| 729 | These averaged values are used to update \textbf{M}(t) with both the surface |
---|
| 730 | pressure gradient and the Coriolis force, therefore providing the $t+\Delta t$ |
---|
| 731 | velocity. The model time stepping scheme can then be achieved by a baroclinic |
---|
| 732 | leap-frog time step that carries the surface height from $t-\Delta t$ to $t+\Delta t$. } |
---|
[707] | 733 | \end{center} |
---|
| 734 | \end{figure} |
---|
| 735 | %> > > > > > > > > > > > > > > > > > > > > > > > > > > > |
---|
| 736 | |
---|
[817] | 737 | The split-explicit formulation has a damping effect on external gravity waves, |
---|
| 738 | which is weaker damping than for the filtered free surface but still significant as |
---|
| 739 | shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave. |
---|
[707] | 740 | |
---|
| 741 | %------------------------------------------------------------- |
---|
| 742 | % Filtered formulation |
---|
| 743 | %------------------------------------------------------------- |
---|
| 744 | \subsubsection{Filtered formulation (\key{dynspg\_flt})} |
---|
| 745 | \label{DYN_spg_flt} |
---|
| 746 | |
---|
[817] | 747 | The filtered formulation follows the \citet{Roullet2000} implementation. The extra |
---|
| 748 | term introduced in the equations (see {\S}I.2.2) is solved implicitly. The elliptic |
---|
| 749 | solvers available in the code are documented in \S\ref{MISC}. The amplitude of |
---|
| 750 | the extra term is given by the namelist variable \np{rnu}. The default value is 1, |
---|
| 751 | as recommended by \citet{Roullet2000} |
---|
[707] | 752 | |
---|
[817] | 753 | \gmcomment{\np{rnu}=1 to be suppressed from namelist !} |
---|
[707] | 754 | |
---|
| 755 | %------------------------------------------------------------- |
---|
| 756 | % Non-linear free surface formulation |
---|
| 757 | %------------------------------------------------------------- |
---|
| 758 | \subsection{Non-linear free surface formulation (\key{vvl})} |
---|
| 759 | \label{DYN_spg_vvl} |
---|
| 760 | |
---|
[817] | 761 | In the non-linear free surface formulation, the variations of volume are fully |
---|
| 762 | taken into account. This option is presented in a report \citep{Levier2007} |
---|
| 763 | available on the \NEMO web site. The three time-stepping methods (explicit, |
---|
| 764 | split-explicit and filtered) are the same as in \S\ref{DYN_spg_linear} except |
---|
| 765 | that the ocean depth is now time-dependent. In particular, this means that |
---|
| 766 | in the filtered case, the matrix to be inverted has to be recomputed at each |
---|
| 767 | time-step. |
---|
[707] | 768 | |
---|
| 769 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 770 | % Rigid-lid formulation |
---|
| 771 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 772 | \subsection{Rigid-lid formulation (\key{dynspg\_rl})} |
---|
| 773 | \label{DYN_spg_rl} |
---|
| 774 | |
---|
| 775 | With the rigid lid formulation, an elliptic equation has to be solved for |
---|
| 776 | the barotropic streamfunction. For consistency this equation is obtained by |
---|
| 777 | taking the discrete curl of the discrete vertical sum of the discrete |
---|
| 778 | momentum equation: |
---|
| 779 | \begin{equation}\label{Eq_dynspg_rl} |
---|
| 780 | \frac{1}{\rho _o }\nabla _h p_s \equiv \left( {{\begin{array}{*{20}c} |
---|
| 781 | {\overline M_u +\frac{1}{H\;e_2 } \delta_ j \left[ \partial_t \psi \right]} \\ |
---|
| 782 | \\ |
---|
| 783 | {\overline M_v -\frac{1}{H\;e_1 } \delta_i \left[ \partial_t \psi \right]} \\ |
---|
| 784 | \end{array} }} \right) |
---|
| 785 | \end{equation} |
---|
| 786 | |
---|
[817] | 787 | Here ${\rm {\bf M}}= \left( M_u,M_v \right)$ represents the collected |
---|
| 788 | contributions of nonlinear, viscous and hydrostatic pressure gradient terms in |
---|
| 789 | \eqref{Eq_PE_dyn} and the overbar indicates a vertical average over the |
---|
| 790 | whole water column (i.e. from $z=-H$, the ocean bottom, to $z=0$, the rigid-lid). |
---|
| 791 | The time derivative of $\psi$ is the solution of an elliptic equation: |
---|
[707] | 792 | \begin{multline} \label{Eq_bsf} |
---|
| 793 | \delta_{i+1/2} \left[ \frac{e_{2v}}{H_v\;e_{1v}} \delta_{i} \left[ \partial_t \psi \right] \right] |
---|
| 794 | + \delta_{j+1/2} \left[ \frac{e_{1u}}{H_u\;e_{2u}} \delta_{j} \left[ \partial_t \psi \right] \right] |
---|
| 795 | \\ = |
---|
| 796 | \delta_{i+1/2} \left[ e_{2v} M_v \right] |
---|
| 797 | - \delta_{j+1/2} \left[ e_{1u} M_u \right] |
---|
| 798 | \end{multline} |
---|
| 799 | |
---|
[817] | 800 | The elliptic solvers available in the code are documented in \S\ref{MISC}). |
---|
| 801 | The boundary conditions must be given on all separate landmasses (islands), |
---|
| 802 | which is done by integrating the vorticity equation around each island. This |
---|
| 803 | requires identifying each island in the bathymetry file, a cumbersome task. |
---|
| 804 | This explains why the rigid lid option is not recommended with complex |
---|
| 805 | domains such as the global ocean. Parameters jpisl (number of islands) and |
---|
| 806 | jpnisl (maximum number of points per island) of the \hf{par\_oce} file are |
---|
| 807 | related to this option. |
---|
[707] | 808 | |
---|
| 809 | |
---|
| 810 | % ================================================================ |
---|
| 811 | % Lateral diffusion term |
---|
| 812 | % ================================================================ |
---|
[817] | 813 | \section [Lateral diffusion term (\textit{dynldf})] |
---|
| 814 | {Lateral diffusion term (\mdl{dynldf})} |
---|
[707] | 815 | \label{DYN_ldf} |
---|
| 816 | %------------------------------------------nam_dynldf---------------------------------------------------- |
---|
| 817 | \namdisplay{nam_dynldf} |
---|
| 818 | %------------------------------------------------------------------------------------------------------------- |
---|
| 819 | |
---|
[817] | 820 | The options available for lateral diffusion are for the choice of laplacian |
---|
| 821 | (rotated or not) or biharmonic operators. The coefficients may be constant |
---|
| 822 | or spatially variable; the description of the coefficients is found in the chapter |
---|
| 823 | on lateralphysics (Chap.\ref{LDF}). The lateral diffusion of momentum is |
---|
| 824 | evaluated using a forward scheme, i.e. the velocity appearing in its expression |
---|
| 825 | is the \textit{before} velocity in time, except for the pure vertical component |
---|
| 826 | that appears when a tensor of rotation is used. This latter term is solved |
---|
| 827 | implicitly together with the vertical diffusion term (see \S\ref{DOM_nxt}) |
---|
| 828 | |
---|
[707] | 829 | At the lateral boundaries either free slip, no slip or partial slip boundary |
---|
[817] | 830 | conditions are applied according to the user's choice (see Chap.\ref{LBC}). |
---|
[707] | 831 | |
---|
| 832 | % ================================================================ |
---|
[817] | 833 | \subsection [Iso-level laplacian operator (\np{ln\_dynldf\_lap}) ] |
---|
| 834 | {Iso-level laplacian operator (\np{ln\_dynldf\_lap}=.true.)} |
---|
[707] | 835 | \label{DYN_ldf_lap} |
---|
| 836 | |
---|
| 837 | For lateral iso-level diffusion, the discrete operator is: |
---|
| 838 | \begin{equation} \label{Eq_dynldf_lap} |
---|
| 839 | \left\{ \begin{aligned} |
---|
| 840 | D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta _{i+1/2} \left[ {A_T^{lm} |
---|
| 841 | \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta _j \left[ |
---|
| 842 | {A_f^{lm} \;e_{3f} \zeta } \right] \\ |
---|
| 843 | \\ |
---|
| 844 | D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta _{j+1/2} \left[ {A_T^{lm} |
---|
| 845 | \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta _i \left[ |
---|
| 846 | {A_f^{lm} \;e_{3f} \zeta } \right] \\ |
---|
| 847 | \end{aligned} \right. |
---|
| 848 | \end{equation} |
---|
| 849 | |
---|
| 850 | As explained in \S\ref{PE_ldf}, this formulation (as the gradient of a divergence |
---|
| 851 | and curl of the vorticity) preserves symmetry and ensures a complete |
---|
[817] | 852 | separation between the vorticity and divergence parts. Note that in the full step |
---|
| 853 | $z$-coordinate (\key{zco} is defined), $e_{3u} =e_{3v} =e_{3f}$ so that they |
---|
| 854 | cancel in the rotational part of \eqref{Eq_dynldf_lap}. |
---|
[707] | 855 | |
---|
| 856 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 857 | % Rotated laplacian operator |
---|
| 858 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 859 | \subsection [Rotated laplacian operator (\np{ln\_dynldf\_iso}) ] |
---|
| 860 | {Rotated laplacian operator (\np{ln\_dynldf\_iso}=.true.)} |
---|
[707] | 861 | \label{DYN_ldf_iso} |
---|
| 862 | |
---|
[817] | 863 | A rotation of the lateral momentum diffusive operator is needed in several cases: |
---|
| 864 | for iso-neutral diffusion in $z$-coordinate (\np{ln\_dynldf\_iso}=.true.) and for |
---|
| 865 | either iso-neutral (\np{ln\_dynldf\_iso}=.true.) or geopotential |
---|
| 866 | (\np{ln\_dynldf\_hor}=.true.) diffusion in $s$-coordinate. In the partial step |
---|
[707] | 867 | case, coordinates are horizontal excepted at the deepest level and no |
---|
[817] | 868 | rotation is performed when \np{ln\_dynldf\_hor}=.true.. The diffusive operator |
---|
[707] | 869 | is defined simply as the divergence of down gradient momentum fluxes on each |
---|
[817] | 870 | momentum component. It must be emphasized that this formulation ignores |
---|
[707] | 871 | constraints on the stress tensor such as symmetry. The resulting discrete |
---|
| 872 | representation is: |
---|
| 873 | \begin{equation} \label{Eq_dyn_ldf_iso} |
---|
| 874 | \begin{split} |
---|
| 875 | D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ |
---|
| 876 | & \left\{\quad {\delta _{i+1/2} \left[ {A_T^{lm} \left( |
---|
| 877 | {\frac{e_{2T} \; e_{3T} }{e_{1T} } \,\delta _{i}[u] |
---|
| 878 | -e_{2T} \; r_{1T} \,\overline{\overline {\delta _{k+1/2}[u]}}^{\,i,\,k}} |
---|
| 879 | \right)} \right]} \right. |
---|
| 880 | \\ |
---|
| 881 | & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} |
---|
| 882 | }\,\delta _{j+1/2} [u] - e_{1f}\, r_{2f} |
---|
| 883 | \,\overline{\overline {\delta _{k+1/2} [u]}} ^{\,j+1/2,\,k}} |
---|
| 884 | \right)} \right] |
---|
| 885 | \\ |
---|
| 886 | &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline |
---|
| 887 | {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } |
---|
| 888 | \right.} \right. |
---|
| 889 | \\ |
---|
| 890 | & \ \qquad \qquad \qquad \quad\ |
---|
| 891 | - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} |
---|
| 892 | \\ |
---|
| 893 | & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ |
---|
| 894 | +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} |
---|
| 895 | \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} |
---|
| 896 | \\ |
---|
| 897 | \\ |
---|
| 898 | D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\ |
---|
| 899 | & \left\{\quad {\delta _{i+1/2} \left[ {A_f^{lm} \left( |
---|
| 900 | {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta _{i+1/2}[v] |
---|
| 901 | -e_{2f} \; r_{1f} \,\overline{\overline {\delta _{k+1/2}[v]}}^{\,i+1/2,\,k}} |
---|
| 902 | \right)} \right]} \right. |
---|
| 903 | \\ |
---|
| 904 | & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1T}\,e_{3T} }{e_{2T} |
---|
| 905 | }\,\delta _{j} [v] - e_{1T}\, r_{2T} |
---|
| 906 | \,\overline{\overline {\delta _{k+1/2} [v]}} ^{\,j,\,k}} |
---|
| 907 | \right)} \right] |
---|
| 908 | \\ |
---|
| 909 | & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline |
---|
| 910 | {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. |
---|
| 911 | \\ |
---|
| 912 | & \ \qquad \qquad \qquad \quad\ |
---|
| 913 | - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} |
---|
| 914 | \\ |
---|
| 915 | & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ |
---|
| 916 | +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} |
---|
| 917 | \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} |
---|
| 918 | \end{split} |
---|
| 919 | \end{equation} |
---|
[817] | 920 | where $r_1$ and $r_2$ are the slopes between the surface along which the |
---|
| 921 | diffusive operator acts and the surface of computation ($z$- or $s$-surfaces). |
---|
| 922 | The way these slopes are evaluated is given in the lateral physics chapter |
---|
| 923 | (Chap.\ref{LDF}). |
---|
[707] | 924 | |
---|
| 925 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 926 | % Iso-level bilaplacian operator |
---|
| 927 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 928 | \subsection [Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap})] |
---|
| 929 | {Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap}=.true.)} |
---|
[707] | 930 | \label{DYN_ldf_bilap} |
---|
| 931 | |
---|
| 932 | The lateral fourth order operator formulation on momentum is obtained by |
---|
[817] | 933 | applying \eqref{Eq_dynldf_lap} twice. It requires an additional assumption on |
---|
| 934 | boundary conditions: the first derivative term normal to the coast depends on |
---|
| 935 | the free or no-slip lateral boundary conditions chosen, while the third |
---|
| 936 | derivative terms normal to the coast are set to zero (see Chap.\ref{LBC}). |
---|
| 937 | %%% |
---|
| 938 | \gmcomment{add a remark on the the change in the position of the coefficient} |
---|
| 939 | %%% |
---|
[707] | 940 | |
---|
| 941 | % ================================================================ |
---|
| 942 | % Vertical diffusion term |
---|
| 943 | % ================================================================ |
---|
[817] | 944 | \section [Vertical diffusion term (\mdl{dynzdf})] |
---|
| 945 | {Vertical diffusion term (\mdl{dynzdf})} |
---|
[707] | 946 | \label{DYN_zdf} |
---|
| 947 | %----------------------------------------------namzdf------------------------------------------------------ |
---|
| 948 | \namdisplay{namzdf} |
---|
| 949 | %------------------------------------------------------------------------------------------------------------- |
---|
| 950 | |
---|
[817] | 951 | The large vertical diffusion coefficient found in the surface mixed layer together |
---|
| 952 | with high vertical resolution implies that in the case of explicit time stepping there |
---|
| 953 | would be too restrictive a constraint on the time step. Two time stepping schemes |
---|
| 954 | can be used for the vertical diffusion term : $(a)$ a forward time differencing |
---|
| 955 | scheme (\np{ln\_zdfexp}=.true.) using a time splitting technique |
---|
| 956 | (\np{n\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme |
---|
| 957 | (\np{ln\_zdfexp}=.false.) (see \S\ref{DOM_nxt}). Note that namelist variables |
---|
| 958 | \np{ln\_zdfexp} and \np{n\_zdfexp} apply to both tracers and dynamics. |
---|
[707] | 959 | |
---|
| 960 | The formulation of the vertical subgrid scale physics is the same whatever |
---|
[817] | 961 | the vertical coordinate is. The vertical diffusion operators given by |
---|
[707] | 962 | \eqref{Eq_PE_zdf} take the following semi-discrete space form: |
---|
| 963 | \begin{equation} \label{Eq_dynzdf} |
---|
| 964 | \left\{ \begin{aligned} |
---|
| 965 | D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta _k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } |
---|
| 966 | \ \delta _{k+1/2} [\,u\,] \right] \\ |
---|
| 967 | \\ |
---|
| 968 | D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta _k \left[ \frac{A_{vw}^{vm} }{e_{3vw} } |
---|
| 969 | \ \delta _{k+1/2} [\,v\,] \right] |
---|
| 970 | \end{aligned} \right. |
---|
| 971 | \end{equation} |
---|
[817] | 972 | where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and |
---|
| 973 | diffusivity coefficients. The way these coefficients are evaluated |
---|
[707] | 974 | depends on the vertical physics used (see \S\ref{ZDF}). |
---|
| 975 | |
---|
| 976 | The surface boundary condition on momentum is given by the stress exerted by |
---|
| 977 | the wind. At the surface, the momentum fluxes are prescribed as the boundary |
---|
| 978 | condition on the vertical turbulent momentum fluxes, |
---|
| 979 | \begin{equation} \label{Eq_dynzdf_sbc} |
---|
| 980 | \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} |
---|
| 981 | = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v } |
---|
| 982 | \end{equation} |
---|
[817] | 983 | where $\left( \tau _u ,\tau _v \right)$ are the two components of the wind stress |
---|
| 984 | vector in the (\textbf{i},\textbf{j}) coordinate system. The high mixing coefficients |
---|
| 985 | in the surface mixed layer ensure that the surface wind stress is distributed in |
---|
[707] | 986 | the vertical over the mixed layer depth. If the vertical mixing coefficient |
---|
| 987 | is small (when no mixed layer scheme is used) the surface stress enters only |
---|
| 988 | the top model level, as a body force. The surface wind stress is calculated |
---|
[817] | 989 | in the surface module routines (SBC, see Chap.\ref{SBC}) |
---|
[707] | 990 | |
---|
[817] | 991 | The turbulent flux of momentum at the bottom of the ocean is specified through |
---|
[1224] | 992 | a bottom friction parameterisation (see \S\ref{ZDF_bfr}) |
---|
[707] | 993 | |
---|
| 994 | % ================================================================ |
---|
| 995 | % External Forcing |
---|
| 996 | % ================================================================ |
---|
| 997 | \section{External Forcings} |
---|
| 998 | \label{DYN_forcing} |
---|
| 999 | |
---|
| 1000 | Besides the surface and bottom stresses (see the above section) which are |
---|
| 1001 | introduced as boundary conditions on the vertical mixing, two other forcings |
---|
| 1002 | enter the dynamical equations. |
---|
| 1003 | |
---|
[817] | 1004 | One is the effect of atmospheric pressure on the ocean dynamics (to be |
---|
| 1005 | introduced later). |
---|
[707] | 1006 | |
---|
[817] | 1007 | Another forcing term is the tidal potential, which will be introduced in the |
---|
| 1008 | reference version soon. |
---|
[707] | 1009 | |
---|
| 1010 | % ================================================================ |
---|
| 1011 | % Time evolution term |
---|
| 1012 | % ================================================================ |
---|
[817] | 1013 | \section [Time evolution term (\textit{dynnxt})] |
---|
| 1014 | {Time evolution term (\mdl{dynnxt})} |
---|
[707] | 1015 | \label{DYN_nxt} |
---|
| 1016 | |
---|
| 1017 | %----------------------------------------------namdom---------------------------------------------------- |
---|
| 1018 | \namdisplay{namdom} |
---|
| 1019 | %------------------------------------------------------------------------------------------------------------- |
---|
| 1020 | |
---|
[817] | 1021 | The general framework for dynamics time stepping is a leap-frog scheme, |
---|
| 1022 | $i.e.$ a three level centred time scheme associated with an Asselin time filter |
---|
| 1023 | (cf. \S\ref{DOM_nxt}) |
---|
[707] | 1024 | \begin{equation} \label{Eq_dynnxt} |
---|
| 1025 | \begin{split} |
---|
| 1026 | &u^{t+\Delta t} = u^{t-\Delta t} + 2 \, \Delta t \ \text{RHS}_u^t \\ |
---|
| 1027 | \\ |
---|
| 1028 | &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\Delta t} -2u^t+u^{t+\Delta t}} \right] |
---|
| 1029 | \end{split} |
---|
| 1030 | \end{equation} |
---|
[817] | 1031 | where RHS is the right hand side of the momentum equation, the subscript $f$ |
---|
| 1032 | denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is |
---|
| 1033 | initialized as \np{atfp} (namelist parameter). Its default value is \np{atfp} = 0.1. |
---|
[707] | 1034 | |
---|
[817] | 1035 | Note that whith the filtered free surface, the update of the \textit{next} velocities |
---|
| 1036 | is done in the \mdl{dynsp\_flt} module, and only the swap of arrays |
---|
| 1037 | and Asselin filtering is done in \mdl{dynnxt.} |
---|
[707] | 1038 | |
---|
| 1039 | % ================================================================ |
---|
| 1040 | % Diagnostic variables |
---|
| 1041 | % ================================================================ |
---|
| 1042 | \section{Diagnostic variables ($\zeta$, $\chi$, $w$)} |
---|
| 1043 | \label{DYN_divcur_wzv} |
---|
| 1044 | |
---|
| 1045 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 1046 | % Horizontal divergence and relative vorticity |
---|
[707] | 1047 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 1048 | \subsection [Horizontal divergence and relative vorticity (\textit{divcur})] |
---|
| 1049 | {Horizontal divergence and relative vorticity (\mdl{divcur})} |
---|
[707] | 1050 | \label{DYN_divcur} |
---|
| 1051 | |
---|
[817] | 1052 | The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows: |
---|
[707] | 1053 | \begin{equation} \label{Eq_divcur_cur} |
---|
| 1054 | \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right] |
---|
[817] | 1055 | -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) |
---|
[707] | 1056 | \end{equation} |
---|
| 1057 | |
---|
[817] | 1058 | The horizontal divergence is defined at a $T$-point. It is given by: |
---|
[707] | 1059 | \begin{equation} \label{Eq_divcur_div} |
---|
| 1060 | \chi =\frac{1}{e_{1T}\,e_{2T}\,e_{3T} } |
---|
| 1061 | \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right] |
---|
| 1062 | +\delta _j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right) |
---|
| 1063 | \end{equation} |
---|
| 1064 | |
---|
[1224] | 1065 | Note that in the $z$-coordinate with full step (\key{zco} is defined), |
---|
| 1066 | $e_{3u} =e_{3v} =e_{3f}$ so that they cancel in \eqref{Eq_divcur_div}. |
---|
[707] | 1067 | |
---|
[817] | 1068 | Note also that whereas the vorticity have the same discrete expression in $z$- |
---|
| 1069 | and $s$-coordinate, its physical meaning is not identical. $\zeta$ is a pseudo |
---|
| 1070 | vorticity along $s$-surfaces (only pseudo because $(u,v)$ are still defined along |
---|
| 1071 | geopotential surfaces, but are no more necessary defined at the same depth). |
---|
[707] | 1072 | |
---|
[817] | 1073 | The vorticity and divergence at the \textit{before} step are used in the computation |
---|
| 1074 | of the horizontal diffusion of momentum. Note that because they have been |
---|
| 1075 | calculated prior to the Asselin filtering of the \textit{before} velocities, the |
---|
| 1076 | \textit{before} vorticity and divergence arrays must be included in the restart file |
---|
| 1077 | to ensure perfect restartability. The vorticity and divergence at the \textit{now} |
---|
| 1078 | time step are used for the computation of the nonlinear advection and of the |
---|
| 1079 | vertical velocity respectively. |
---|
[707] | 1080 | |
---|
| 1081 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1082 | % Vertical Velocity |
---|
| 1083 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 1084 | \subsection [Vertical velocity (\textit{wzvmod})] |
---|
| 1085 | {Vertical velocity (\mdl{wzvmod})} |
---|
[707] | 1086 | \label{DYN_wzv} |
---|
| 1087 | |
---|
| 1088 | The vertical velocity is computed by an upward integration of the horizontal |
---|
| 1089 | divergence from the bottom : |
---|
| 1090 | |
---|
| 1091 | \begin{equation} \label{Eq_wzv} |
---|
| 1092 | \left\{ \begin{aligned} |
---|
| 1093 | &\left. w \right|_{3/2} \quad= 0 \\ |
---|
| 1094 | \\ |
---|
| 1095 | &\left. w \right|_{k+1/2} = \left. w \right|_{k+1/2} + e_{3t}\; \left. \chi \right|_k |
---|
| 1096 | \end{aligned} \right. |
---|
| 1097 | \end{equation} |
---|
| 1098 | |
---|
[817] | 1099 | With a free surface, the top vertical velocity is non-zero, due to the |
---|
| 1100 | freshwater forcing and the variations of the free surface elevation. With a |
---|
| 1101 | linear free surface or with a rigid lid, the upper boundary condition |
---|
| 1102 | applies at a fixed level $z=0$. Note that in the rigid-lid case (\key{dynspg\_rl} |
---|
| 1103 | is defined), the surface boundary condition ($\left. w \right|_\text{surface}=0)$ is |
---|
| 1104 | automatically achieved at least at computer accuracy, due to the the way the |
---|
| 1105 | surface pressure gradient is expressed in discrete form (Appendix~\ref{Apdx_C}). |
---|
[707] | 1106 | |
---|
| 1107 | Note also that whereas the vertical velocity has the same discrete |
---|
[817] | 1108 | expression in $z$- and $s$-coordinate, its physical meaning is not the same: |
---|
| 1109 | in the second case, $w$ is the velocity normal to the $s$-surfaces. |
---|
[707] | 1110 | |
---|
| 1111 | With the variable volume option, the calculation of the vertical velocity is |
---|
[817] | 1112 | modified (see \citet{Levier2007}, report available on the \NEMO web site). |
---|
[707] | 1113 | |
---|
[1224] | 1114 | % ================================================================ |
---|