Changeset 10354 for NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex
- Timestamp:
- 2018-11-21T17:59:55+01:00 (5 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex
r10146 r10354 15 15 %--------------------------------------------------------------------------------------------------------- 16 16 17 Two scheme are available to perform the iso-neutral diffusion. 18 If the namelist logical \np{ln\_traldf\_triad} is set true, 19 \NEMO updates both active and passive tracers using the Griffies triad representation 20 of iso-neutral diffusion and the eddy-induced advective skew (GM) fluxes. 21 If the namelist logical \np{ln\_traldf\_iso} is set true, 22 the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}). 23 In the present implementation of the Griffies scheme, 17 Two scheme are available to perform the iso-neutral diffusion. 18 If the namelist logical \np{ln\_traldf\_triad} is set true, 19 \NEMO updates both active and passive tracers using the Griffies triad representation of iso-neutral diffusion and 20 the eddy-induced advective skew (GM) fluxes. 21 If the namelist logical \np{ln\_traldf\_iso} is set true, 22 the filtered version of Cox's original scheme (the Standard scheme) is employed (\autoref{sec:LDF_slp}). 23 In the present implementation of the Griffies scheme, 24 24 the advective skew fluxes are implemented even if \np{ln\_traldf\_eiv} is false. 25 25 26 Values of iso-neutral diffusivity and GM coefficient are set as 27 described in \autoref{sec:LDF_coef}. Note that when GM fluxes are used, 28 the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS, 26 Values of iso-neutral diffusivity and GM coefficient are set as described in \autoref{sec:LDF_coef}. 27 Note that when GM fluxes are used, the eddy-advective (GM) velocities are output for diagnostic purposes using xIOS, 29 28 even though the eddy advection is accomplished by means of the skew fluxes. 30 29 … … 32 31 The options specific to the Griffies scheme include: 33 32 \begin{description}[font=\normalfont] 34 \item[\np{ln\_triad\_iso}] See \autoref{sec:taper}. If this is set false (the default), then 35 `iso-neutral' mixing is accomplished within the surface mixed-layer 36 along slopes linearly decreasing with depth from the value immediately below 37 the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}). 38 This is the same treatment as used in the default implementation \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}. 39 Where \np{ln\_triad\_iso} is set true, the vertical skew flux is further reduced 40 to ensure no vertical buoyancy flux, giving an almost pure 41 horizontal diffusive tracer flux within the mixed layer. This is similar to 42 the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper} 43 \item[\np{ln\_botmix\_triad}] See \autoref{sec:iso_bdry}. 33 \item[\np{ln\_triad\_iso}] 34 See \autoref{sec:taper}. 35 If this is set false (the default), 36 then `iso-neutral' mixing is accomplished within the surface mixed-layer along slopes linearly decreasing with 37 depth from the value immediately below the mixed-layer to zero (flat) at the surface (\autoref{sec:lintaper}). 38 This is the same treatment as used in the default implementation 39 \autoref{subsec:LDF_slp_iso}; \autoref{fig:eiv_slp}. 40 Where \np{ln\_triad\_iso} is set true, 41 the vertical skew flux is further reduced to ensure no vertical buoyancy flux, 42 giving an almost pure horizontal diffusive tracer flux within the mixed layer. 43 This is similar to the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper} 44 \item[\np{ln\_botmix\_triad}] 45 See \autoref{sec:iso_bdry}. 44 46 If this is set false (the default) then the lateral diffusive fluxes 45 47 associated with triads partly masked by topography are neglected. 46 48 If it is set true, however, then these lateral diffusive fluxes are applied, 47 49 giving smoother bottom tracer fields at the cost of introducing diapycnal mixing. 48 \item[\np{rn\_sw\_triad}] blah blah to be added.... 50 \item[\np{rn\_sw\_triad}] 51 blah blah to be added.... 49 52 \end{description} 50 53 The options shared with the Standard scheme include: … … 56 59 \section{Triad formulation of iso-neutral diffusion} 57 60 \label{sec:iso} 58 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, 61 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, 59 62 but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 60 63 61 64 \subsection{Iso-neutral diffusion operator} 62 The iso-neutral second order tracer diffusive operator for small 63 angles between iso-neutral surfaces and geopotentials is given by 64 \autoref{eq:iso_tensor_1}: 65 The iso-neutral second order tracer diffusive operator for small angles between 66 iso-neutral surfaces and geopotentials is given by \autoref{eq:iso_tensor_1}: 65 67 \begin{subequations} \label{eq:iso_tensor_1} 66 68 \begin{equation} … … 94 96 % {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 95 97 % \end{array} }} \right) 96 98 Here \autoref{eq:PE_iso_slopes} 97 99 \begin{align*} 98 100 r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} … … 104 106 }{\partial k} \right)^{-1} 105 107 \end{align*} 106 is the $i$-component of the slope of the iso-neutral surface relative to the computational 107 surface, and $r_2$ is the $j$-component. 108 109 We will find it useful to consider the fluxes per unit area in $i,j,k$ 110 space; we write 108 is the $i$-component of the slope of the iso-neutral surface relative to the computational surface, 109 and $r_2$ is the $j$-component. 110 111 We will find it useful to consider the fluxes per unit area in $i,j,k$ space; we write 111 112 \begin{equation} 112 113 \label{eq:Fijk} 113 114 \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). 114 115 \end{equation} 115 Additionally, we will sometimes write the contributions towards the 116 fluxes $\vect{f}$ and $\vect{F}_{\mathrm{iso}}$ from the component 117 $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, with 118 $f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc. 116 Additionally, we will sometimes write the contributions towards the fluxes $\vect{f}$ and 117 $\vect{F}_{\mathrm{iso}}$ from the component $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, 118 with $f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc. 119 119 120 120 The off-diagonal terms of the small angle diffusion tensor 121 \autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along the122 $i$- and $j$-directions resulting from the vertical tracer gradient:121 \autoref{eq:iso_tensor_1}, \autoref{eq:iso_tensor_2} produce skew-fluxes along 122 the $i$- and $j$-directions resulting from the vertical tracer gradient: 123 123 \begin{align} 124 124 \label{eq:i13c} … … 129 129 \end{align} 130 130 131 The vertical diffusive flux associated with the $_{33}$ 132 component of the small angle diffusion tensor is 131 The vertical diffusive flux associated with the $_{33}$ component of the small angle diffusion tensor is 133 132 \begin{equation} 134 133 \label{eq:i33c} … … 136 135 \end{equation} 137 136 138 Since there are no cross terms involving $r_1$ and $r_2$ in the above, we can 139 consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$ 140 planes, just adding together the vertical components from each 141 plane. The following description will describe the fluxes on the $i$-$k$ 142 plane. 143 144 There is no natural discretization for the $i$-component of the 145 skew-flux, \autoref{eq:i13c}, as 146 although it must be evaluated at $u$-points, it involves vertical 147 gradients (both for the tracer and the slope $r_1$), defined at 148 $w$-points. Similarly, the vertical skew flux, \autoref{eq:i31c}, is evaluated at 149 $w$-points but involves horizontal gradients defined at $u$-points. 137 Since there are no cross terms involving $r_1$ and $r_2$ in the above, 138 we can consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$ planes, 139 just adding together the vertical components from each plane. 140 The following description will describe the fluxes on the $i$-$k$ plane. 141 142 There is no natural discretization for the $i$-component of the skew-flux, \autoref{eq:i13c}, 143 as although it must be evaluated at $u$-points, 144 it involves vertical gradients (both for the tracer and the slope $r_1$), defined at $w$-points. 145 Similarly, the vertical skew flux, \autoref{eq:i31c}, 146 is evaluated at $w$-points but involves horizontal gradients defined at $u$-points. 150 147 151 148 \subsection{Standard discretization} 152 149 The straightforward approach to discretize the lateral skew flux 153 \autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 154 into OPA, \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical 155 gradient at the $u$-point from the average of the four surrounding 156 vertical tracer gradients, and multiply this by a mean slope at the 157 $u$-point, calculated from the averaged surrounding vertical density 158 gradients. The total area-integrated skew-flux (flux per unit area in 159 $ijk$ space) from tracer cell $i,k$ 160 to $i+1,k$, noting that the $e_{{3}_{i+1/2}^k}$ in the area 161 $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with 162 the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer 163 gradient, is then \autoref{eq:tra_ldf_iso} 150 \autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA, 151 \autoref{eq:tra_ldf_iso}, is to calculate a mean vertical gradient at the $u$-point from 152 the average of the four surrounding vertical tracer gradients, and multiply this by a mean slope at the $u$-point, 153 calculated from the averaged surrounding vertical density gradients. 154 The total area-integrated skew-flux (flux per unit area in $ijk$ space) from tracer cell $i,k$ to $i+1,k$, 155 noting that the $e_{{3}_{i+1/2}^k}$ in the area $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with 156 the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer gradient, is then \autoref{eq:tra_ldf_iso} 164 157 \begin{equation*} 165 158 \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k … … 173 166 \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}}, 174 167 \end{equation*} 175 and here and in the following we drop the $^{lT}$ superscript from 176 $\Alt$ for simplicity. 177 Unfortunately the resulting combination $\overline{\overline{\delta_k 178 \bullet}}^{\,i,k}$ of a $k$ average and a $k$ difference %of the tracer 179 reduces to $\bullet_{k+1}-\bullet_{k-1}$, so two-grid-point oscillations are 180 invisible to this discretization of the iso-neutral operator. These 181 \emph{computational modes} will not be damped by this operator, and 182 may even possibly be amplified by it. Consequently, applying this 183 operator to a tracer does not guarantee the decrease of its 184 global-average variance. To correct this, we introduced a smoothing of 185 the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). This 186 technique works for $T$ and $S$ in so far as they are active tracers 187 ($i.e.$ they enter the computation of density), but it does not work 188 for a passive tracer. 168 and here and in the following we drop the $^{lT}$ superscript from $\Alt$ for simplicity. 169 Unfortunately the resulting combination $\overline{\overline{\delta_k\bullet}}^{\,i,k}$ of a $k$ average and 170 a $k$ difference of the tracer reduces to $\bullet_{k+1}-\bullet_{k-1}$, 171 so two-grid-point oscillations are invisible to this discretization of the iso-neutral operator. 172 These \emph{computational modes} will not be damped by this operator, and may even possibly be amplified by it. 173 Consequently, applying this operator to a tracer does not guarantee the decrease of its global-average variance. 174 To correct this, we introduced a smoothing of the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). 175 This technique works for $T$ and $S$ in so far as they are active tracers 176 ($i.e.$ they enter the computation of density), but it does not work for a passive tracer. 189 177 190 178 \subsection{Expression of the skew-flux in terms of triad slopes} 191 \citep{Griffies_al_JPO98} introduce a different discretization of the 192 off-diagonal terms thatnicely solves the problem.179 \citep{Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that 180 nicely solves the problem. 193 181 % Instead of multiplying the mean slope calculated at the $u$-point by 194 182 % the mean vertical gradient at the $u$-point, … … 203 191 \end{center} \end{figure} 204 192 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 205 They get the skew flux from the products of the vertical gradients at 206 each $w$-point surrounding the $u$-point with the corresponding `triad' 207 slope calculated from the lateral density gradient across the $u$-point 208 divided by the vertical density gradient at the same $w$-point as the 209 tracer gradient. See \autoref{fig:ISO_triad}a, where the thick lines 210 denote the tracer gradients, and the thin lines the corresponding 211 triads, with slopes $s_1, \dotsc s_4$. The total area-integrated 212 skew-flux from tracer cell $i,k$ to $i+1,k$ 193 They get the skew flux from the products of the vertical gradients at each $w$-point surrounding the $u$-point with 194 the corresponding `triad' slope calculated from the lateral density gradient across the $u$-point divided by 195 the vertical density gradient at the same $w$-point as the tracer gradient. 196 See \autoref{fig:ISO_triad}a, where the thick lines denote the tracer gradients, 197 and the thin lines the corresponding triads, with slopes $s_1, \dotsc s_4$. 198 The total area-integrated skew-flux from tracer cell $i,k$ to $i+1,k$ 213 199 \begin{multline} 214 200 \label{eq:i13} … … 222 208 _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}, 223 209 \end{multline} 224 where the contributions of the triad fluxes are weighted by areas 225 $a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points 226 rather than the $u$-points. This discretization gives a much closer 227 stencil, and disallows the two-point computational modes. 228 229 The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at the 230 $w$-point $i,k+\hhalf$ is constructed similarly (\autoref{fig:ISO_triad}b) 231 by multiplying lateral tracer gradients from each of the four 232 surrounding $u$-points by the appropriate triad slope: 210 where the contributions of the triad fluxes are weighted by areas $a_1, \dotsc a_4$, 211 and $\Alts$ is now defined at the tracer points rather than the $u$-points. 212 This discretization gives a much closer stencil, and disallows the two-point computational modes. 213 214 The vertical skew flux \autoref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at 215 the $w$-point $i,k+\hhalf$ is constructed similarly (\autoref{fig:ISO_triad}b) by 216 multiplying lateral tracer gradients from each of the four surrounding $u$-points by the appropriate triad slope: 233 217 \begin{multline} 234 218 \label{eq:i31} … … 241 225 242 226 We notate the triad slopes $s_i$ and $s'_i$ in terms of the `anchor point' $i,k$ 243 (appearing in both the vertical and lateral gradient), and the $u$- and244 $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the 245 triad as follows(see also \autoref{fig:ISO_triad}):227 (appearing in both the vertical and lateral gradient), 228 and the $u$- and $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the triad as follows 229 (see also \autoref{fig:ISO_triad}): 246 230 \begin{equation} 247 231 \label{eq:R} … … 253 237 { \alpha_i^k \ \delta_{k+k_p}[T^i] - \beta_i^k \ \delta_{k+k_p}[S^i] }. 254 238 \end{equation} 255 In calculating the slopes of the local neutral surfaces, 256 the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad, 239 In calculating the slopes of the local neutral surfaces, 240 the expansion coefficients $\alpha$ and $\beta$ are evaluated at the anchor points of the triad, 257 241 while the metrics are calculated at the $u$- and $w$-points on the arms. 258 242 … … 261 245 \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells} 262 246 \caption{ \protect\label{fig:qcells} 263 Triad notation for quarter cells. $T$-cells are inside264 boxes, while the $i+\half,k$ $u$-cell is shaded in green and the265 $i,k+\half$ $w$-cell is shaded in pink.}247 Triad notation for quarter cells. $T$-cells are inside boxes, 248 while the $i+\half,k$ $u$-cell is shaded in green and 249 the $i,k+\half$ $w$-cell is shaded in pink.} 266 250 \end{center} \end{figure} 267 251 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 268 252 269 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter 270 cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell. 271 Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation, 272 we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. 273 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) 274 to calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$, 275 and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$. 276 Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used 277 to calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms 278 of a unique triad, and we notate these areas, similarly to the triad slopes, 279 as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, 280 where $e.g.$ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, 253 Each triad $\{_i^{k}\:_{i_p}^{k_p}\}$ is associated (\autoref{fig:qcells}) with the quarter cell that is 254 the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ $u$-cell and the $i,k+k_p$ $w$-cell. 255 Expressing the slopes $s_i$ and $s'_i$ in \autoref{eq:i13} and \autoref{eq:i31} in this notation, 256 we have $e.g.$ \ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. 257 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to 258 calculate the lateral flux along its $u$-arm, at $(i+i_p,k)$, 259 and then again as an $s'$ to calculate the vertical flux along its $w$-arm at $(i,k+k_p)$. 260 Each vertical area $a_i$ used to calculate the lateral flux and horizontal area $a'_i$ used to 261 calculate the vertical flux can also be identified as the area across the $u$- and $w$-arms of a unique triad, 262 and we notate these areas, similarly to the triad slopes, 263 as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, 264 where $e.g.$ in \autoref{eq:i13} $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, 281 265 and in \autoref{eq:i31} $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 282 266 283 267 \subsection{Full triad fluxes} 284 A key property of iso-neutral diffusion is that it should not affect 285 the (locally referenced) density. In particular there should be no 286 lateral or vertical density flux. The lateral density flux disappears so long as the 287 area-integrated lateral diffusive flux from tracer cell $i,k$ to 288 $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the 289 form 268 A key property of iso-neutral diffusion is that it should not affect the (locally referenced) density. 269 In particular there should be no lateral or vertical density flux. 270 The lateral density flux disappears so long as the area-integrated lateral diffusive flux from 271 tracer cell $i,k$ to $i+1,k$ coming from the $_{11}$ term of the diffusion tensor takes the form 290 272 \begin{equation} 291 273 \label{eq:i11} … … 295 277 \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, 296 278 \end{equation} 297 where the areas $a_i$ are as in \autoref{eq:i13}. In this case, 298 separating the total lateral flux, the sum of \autoref{eq:i13} and 299 \autoref{eq:i11}, into triad components, a lateral tracer 300 flux 279 where the areas $a_i$ are as in \autoref{eq:i13}. 280 In this case, separating the total lateral flux, the sum of \autoref{eq:i13} and \autoref{eq:i11}, 281 into triad components, a lateral tracer flux 301 282 \begin{equation} 302 283 \label{eq:latflux-triad} … … 308 289 \right) 309 290 \end{equation} 310 can be identified with each triad. Then, because the 311 same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and 312 ${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both the density gradients 313 in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, the lateral 314 density flux associated with each triad separately disappears. 291 can be identified with each triad. 292 Then, because the same metric factors ${e_{3w}}_{\,i}^{\,k+k_p}$ and ${e_{1u}}_{\,i+i_p}^{\,k}$ are employed for both 293 the density gradients in $ _i^k \mathbb{R}_{i_p}^{k_p}$ and the tracer gradients, 294 the lateral density flux associated with each triad separately disappears. 315 295 \begin{equation} 316 296 \label{eq:latflux-rho} 317 297 {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0 318 298 \end{equation} 319 Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} + 320 \left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from tracer cell $i,k$ 321 to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. 322 323 The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the 324 $_{33}$ component is also expressed in terms of area-weighted 325 squared triad slopes, so the area-integrated vertical flux from tracer 326 cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is 299 Thus the total flux $\left( F_u^{31} \right) ^i _{i,k+\frac{1}{2}} + \left( F_u^{11} \right) ^i _{i,k+\frac{1}{2}}$ from 300 tracer cell $i,k$ to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. 301 302 The squared slope $r_1^2$ in the expression \autoref{eq:i33c} for the $_{33}$ component is also expressed in 303 terms of area-weighted squared triad slopes, 304 so the area-integrated vertical flux from tracer cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is 327 305 \begin{equation} 328 306 \label{eq:i33} … … 333 311 + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right], 334 312 \end{equation} 335 where the areas $a'$ and slopes $s'$ are the same as in 336 \autoref{eq:i31}. 337 Then, separating the total vertical flux, the sum of \autoref{eq:i31} and 338 \autoref{eq:i33}, into triad components, a vertical flux 313 where the areas $a'$ and slopes $s'$ are the same as in \autoref{eq:i31}. 314 Then, separating the total vertical flux, the sum of \autoref{eq:i31} and \autoref{eq:i33}, 315 into triad components, a vertical flux 339 316 \begin{align} 340 317 \label{eq:vertflux-triad} … … 349 326 {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} 350 327 \end{align} 351 may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ 352 associated with a triad then separately disappears (because the 353 lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$ 354 disappears). Consequently the total vertical density flux $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} + 355 \left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from tracer cell $i,k$ 356 to $i,k+1$ must also vanish since it is a sum of four such triad 357 fluxes. 358 359 We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of 360 the $u$-fluxes and $w$-fluxes in 361 \autoref{eq:i31}, \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and 362 \autoref{fig:ISO_triad} to write out the iso-neutral fluxes at $u$- and 363 $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 328 may be associated with each triad. 329 Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ associated with a triad then 330 separately disappears (because the lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (\rho)$ disappears). 331 Consequently the total vertical density flux 332 $\left( F_w^{31} \right)_i ^{k+\frac{1}{2}} + \left( F_w^{33} \right)_i^{k+\frac{1}{2}}$ from 333 tracer cell $i,k$ to $i,k+1$ must also vanish since it is a sum of four such triad fluxes. 334 335 We can explicitly identify (\autoref{fig:qcells}) the triads associated with the $s_i$, $a_i$, 336 and $s'_i$, $a'_i$ used in the definition of the $u$-fluxes and $w$-fluxes in \autoref{eq:i31}, 337 \autoref{eq:i13}, \autoref{eq:i11} \autoref{eq:i33} and \autoref{fig:ISO_triad} to write out 338 the iso-neutral fluxes at $u$- and $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 364 339 %(\autoref{fig:ISO_triad}): 365 340 \begin{flalign} \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv … … 375 350 \label{subsec:variance} 376 351 377 We now require that this operator should not increase the 378 globally-integrated tracer variance. 352 We now require that this operator should not increase the globally-integrated tracer variance. 379 353 %This changes according to 380 354 % \begin{align*} … … 387 361 % + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \ \delta_{k+1/2} [T] \right\} \\ 388 362 % \end{align*} 389 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux 390 $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and 391 a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the 392 $w$-point $i,k+k_p$. The lateral flux drives a net rate of change of 393 variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 363 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across 364 the $u$-point $i+i_p,k$ and a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the $w$-point $i,k+k_p$. 365 The lateral flux drives a net rate of change of variance, 366 summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 394 367 \begin{multline} 395 368 {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+ … … 402 375 \end{aligned} 403 376 \end{multline} 404 while the vertical flux similarly drives a net rate of change of 405 variance summed over the $T$-points $i,k+k_p-\half$ (above) and 406 $i,k+k_p+\half$ (below) of 377 while the vertical flux similarly drives a net rate of change of variance summed over 378 the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 407 379 \begin{equation} 408 380 \label{eq:dvar_iso_k} 409 381 _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 410 382 \end{equation} 411 The total variance tendency driven by the triad is the sum of these 412 two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and 413 $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \autoref{eq:latflux-triad} and 414 \autoref{eq:vertflux-triad}, it is 383 The total variance tendency driven by the triad is the sum of these two. 384 Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with 385 \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, it is 415 386 \begin{multline*} 416 387 -\Alts_i^k\left \{ … … 428 399 \right \}. 429 400 \end{multline*} 430 The key point is then that if we require 431 $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$ 432 to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by 401 The key point is then that if we require $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$ to 402 be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by 433 403 \begin{equation} 434 404 \label{eq:V-A} … … 447 417 \right)^2\leq 0. 448 418 \end{equation} 449 Thus, the constraint \autoref{eq:V-A} ensures that the fluxes (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated 450 with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase 451 the net variance. Since the total fluxes are sums of such fluxes from 452 the various triads, this constraint, applied to all triads, is 453 sufficient to ensure that the globally integrated variance does not 454 increase. 455 456 The expression \autoref{eq:V-A} can be interpreted as a discretization 457 of the global integral 419 Thus, the constraint \autoref{eq:V-A} ensures that the fluxes 420 (\autoref{eq:latflux-triad}, \autoref{eq:vertflux-triad}) associated with 421 a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase the net variance. 422 Since the total fluxes are sums of such fluxes from the various triads, this constraint, applied to all triads, 423 is sufficient to ensure that the globally integrated variance does not increase. 424 425 The expression \autoref{eq:V-A} can be interpreted as a discretization of the global integral 458 426 \begin{equation} 459 427 \label{eq:cts-var} … … 461 429 \int\!\mathbf{F}\cdot\nabla T\, dV, 462 430 \end{equation} 463 where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the 464 lateral and vertical fluxes/unit area 431 where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the lateral and vertical fluxes/unit area 465 432 \[ 466 433 \mathbf{F}=\left( … … 477 444 478 445 \subsection{Triad volumes in Griffes's scheme and in \NEMO} 479 To complete the discretization we now need only specify the triad 480 volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify 481 these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter 482 cells, defined in terms of the distances between $T$, $u$,$f$ and 483 $w$-points. This is the natural discretization of 484 \autoref{eq:cts-var}. The \NEMO model, however, operates with scale 485 factors instead of grid sizes, and scale factors for the quarter 486 cells are not defined. Instead, therefore we simply choose 446 To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. 447 \citet{Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, 448 defined in terms of the distances between $T$, $u$,$f$ and $w$-points. 449 This is the natural discretization of \autoref{eq:cts-var}. 450 The \NEMO model, however, operates with scale factors instead of grid sizes, 451 and scale factors for the quarter cells are not defined. 452 Instead, therefore we simply choose 487 453 \begin{equation} 488 454 \label{eq:V-NEMO} 489 455 _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k, 490 456 \end{equation} 491 as a quarter of the volume of the $u$-cell inside which the triad 492 quarter-cell lies. This has the nice property that when the slopes 493 $\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to 494 $i+1,k$ reduces to the classical form 457 as a quarter of the volume of the $u$-cell inside which the triad quarter-cell lies. 458 This has the nice property that when the slopes $\mathbb{R}$ vanish, 459 the lateral flux from tracer cell $i,k$ to $i+1,k$ reduces to the classical form 495 460 \begin{equation} 496 461 \label{eq:lat-normal} … … 500 465 = -\overline\Alts_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{\,i + 1/2}^{\,k}}. 501 466 \end{equation} 502 In fact if the diffusive coefficient is defined at $u$-points, so that503 we employ $\Alts_{i+i_p}^k$ instead of $\Alts_i^k$ in the definitions of the 504 triad fluxes\autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad},467 In fact if the diffusive coefficient is defined at $u$-points, 468 so that we employ $\Alts_{i+i_p}^k$ instead of $\Alts_i^k$ in the definitions of the triad fluxes 469 \autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, 505 470 we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. 506 471 507 472 \subsection{Summary of the scheme} 508 The iso-neutral fluxes at $u$- and 509 $w$-points are the sums of the triad fluxes that cross the $u$- and 510 $w$-faces \autoref{eq:iso_flux}: 473 The iso-neutral fluxes at $u$- and $w$-points are the sums of the triad fluxes that 474 cross the $u$- and $w$-faces \autoref{eq:iso_flux}: 511 475 \begin{subequations}\label{eq:alltriadflux} 512 476 \begin{flalign}\label{eq:vect_isoflux} … … 545 509 \end{subequations} 546 510 547 511 The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 548 512 each tracer point: 549 513 \begin{equation} \label{eq:iso_operator} D_l^T = \frac{1}{b_T} … … 555 519 The diffusion scheme satisfies the following six properties: 556 520 \begin{description} 557 \item[$\bullet$ horizontal diffusion] The discretization of the558 diffusion operator recovers \autoref{eq:lat-normal} the traditional five-point Laplacian in559 the limit of flat iso-neutral direction:521 \item[$\bullet$ horizontal diffusion] 522 The discretization of the diffusion operator recovers the traditional five-point Laplacian 523 \autoref{eq:lat-normal} in the limit of flat iso-neutral direction: 560 524 \begin{equation} \label{eq:iso_property0} D_l^T = \frac{1}{b_T} \ 561 525 \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; … … 564 528 \end{equation} 565 529 566 \item[$\bullet$ implicit treatment in the vertical] Only tracer values 567 associated with a single water column appear in the expression 568 \autoref{eq:i33} for the $_{33}$ fluxes, vertical fluxes driven by 569 vertical gradients. This is of paramount importance since it means 570 that a time-implicit algorithm can be used to solve the vertical 571 diffusion equation. This is necessary 572 since the vertical eddy 573 diffusivity associated with this term, 530 \item[$\bullet$ implicit treatment in the vertical] 531 Only tracer values associated with a single water column appear in the expression \autoref{eq:i33} for 532 the $_{33}$ fluxes, vertical fluxes driven by vertical gradients. 533 This is of paramount importance since it means that a time-implicit algorithm can be used to 534 solve the vertical diffusion equation. 535 This is necessary since the vertical eddy diffusivity associated with this term, 574 536 \begin{equation} 575 537 \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{ … … 579 541 {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 580 542 \right\}, 581 \end{equation}543 \end{equation} 582 544 (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large. 583 545 584 \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of585 locally referenced potential density is zero. See586 \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}.587 588 \item[$\bullet$ conservation of tracer] The iso-neutral diffusion589 conserves tracer content, $i.e.$546 \item[$\bullet$ pure iso-neutral operator] 547 The iso-neutral flux of locally referenced potential density is zero. 548 See \autoref{eq:latflux-rho} and \autoref{eq:vertflux-triad2}. 549 550 \item[$\bullet$ conservation of tracer] 551 The iso-neutral diffusion conserves tracer content, $i.e.$ 590 552 \begin{equation} \label{eq:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 591 553 b_T \right\} = 0 592 554 \end{equation} 593 This property is trivially satisfied since the iso-neutral diffusive 594 operator is written in flux form. 595 596 \item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion 597 does not increase the tracer variance, $i.e.$ 555 This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form. 556 557 \item[$\bullet$ no increase of tracer variance] 558 The iso-neutral diffusion does not increase the tracer variance, $i.e.$ 598 559 \begin{equation} \label{eq:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 599 560 \ b_T \right\} \leq 0 600 561 \end{equation} 601 The property is demonstrated in 602 \autoref{subsec:variance} above. It is a key property for a diffusion 603 term. It means that it is also a dissipation term, $i.e.$ it 604 dissipates the square of the quantity on which it is applied. It 605 therefore ensures that, when the diffusivity coefficient is large 606 enough, the field on which it is applied becomes free of grid-point 607 noise. 608 609 \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion 610 operator is self-adjoint, $i.e.$ 562 The property is demonstrated in \autoref{subsec:variance} above. 563 It is a key property for a diffusion term. 564 It means that it is also a dissipation term, 565 $i.e.$ it dissipates the square of the quantity on which it is applied. 566 It therefore ensures that, when the diffusivity coefficient is large enough, 567 the field on which it is applied becomes free of grid-point noise. 568 569 \item[$\bullet$ self-adjoint operator] 570 The iso-neutral diffusion operator is self-adjoint, $i.e.$ 611 571 \begin{equation} \label{eq:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 612 572 \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 613 573 \end{equation} 614 In other word, there is no need to develop a specific routine from 615 the adjoint of this operator. We just have to apply the same 616 routine. This property can be demonstrated similarly to the proof of 617 the `no increase of tracer variance' property. The contribution by a 618 single triad towards the left hand side of \autoref{eq:iso_property3}, can 619 be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} 620 and \autoref{eq:dvar_iso_k}. This results in a term similar to 621 \autoref{eq:perfect-square}, 574 In other word, there is no need to develop a specific routine from the adjoint of this operator. 575 We just have to apply the same routine. 576 This property can be demonstrated similarly to the proof of the `no increase of tracer variance' property. 577 The contribution by a single triad towards the left hand side of \autoref{eq:iso_property3}, 578 can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} and \autoref{eq:dvar_iso_k}. 579 This results in a term similar to \autoref{eq:perfect-square}, 622 580 \begin{equation} 623 581 \label{eq:TScovar} … … 634 592 \right). 635 593 \end{equation} 636 This is symmetrical in $T $ and $S$, so exactly the same term arises 637 from the discretization of this triad's contribution towards the 638 RHS of \autoref{eq:iso_property3}. 594 This is symmetrical in $T $ and $S$, so exactly the same term arises from 595 the discretization of this triad's contribution towards the RHS of \autoref{eq:iso_property3}. 639 596 \end{description} 640 597 641 598 \subsection{Treatment of the triads at the boundaries}\label{sec:iso_bdry} 642 The triad slope can only be defined where both the grid boxes centred at 643 the end of the arms exist. Triads that would poke up 644 through the upper ocean surface into the atmosphere, or down into the 645 ocean floor, must be masked out. See \autoref{fig:bdry_triads}. Surface layer triads 646 $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and 647 $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be 648 specified above the ocean surface are masked (\autoref{fig:bdry_triads}a): this ensures that lateral 649 tracer gradients produce no flux through the ocean surface. However, 650 to prevent surface noise, it is customary to retain the $_{11}$ contributions towards 651 the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and 652 $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer 653 fluxes. Similar comments apply to triads that would intersect the 654 ocean floor (\autoref{fig:bdry_triads}b). Note that both near bottom 655 triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 656 $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ 657 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is 658 masked. The associated lateral fluxes (grey-black dashed line) are 659 masked if \np{ln\_botmix\_triad}\forcode{ = .false.}, but left unmasked, 660 giving bottom mixing, if \np{ln\_botmix\_triad}\forcode{ = .true.}. 661 662 The default option \np{ln\_botmix\_triad}\forcode{ = .false.} is suitable when the 663 bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}\forcode{ = 1}), 664 or for simple idealized problems. For setups with topography without 665 bbl mixing, \np{ln\_botmix\_triad}\forcode{ = .true.} may be necessary. 599 The triad slope can only be defined where both the grid boxes centred at the end of the arms exist. 600 Triads that would poke up through the upper ocean surface into the atmosphere, 601 or down into the ocean floor, must be masked out. 602 See \autoref{fig:bdry_triads}. 603 Surface layer triads $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that 604 require density to be specified above the ocean surface are masked (\autoref{fig:bdry_triads}a): 605 this ensures that lateral tracer gradients produce no flux through the ocean surface. 606 However, to prevent surface noise, it is customary to retain the $_{11}$ contributions towards 607 the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; 608 this drives diapycnal tracer fluxes. 609 Similar comments apply to triads that would intersect the ocean floor (\autoref{fig:bdry_triads}b). 610 Note that both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 611 $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, 612 i.e.\ the $i,k+1$ $u$-point is masked. 613 The associated lateral fluxes (grey-black dashed line) are masked if \np{ln\_botmix\_triad}\forcode{ = .false.}, 614 but left unmasked, giving bottom mixing, if \np{ln\_botmix\_triad}\forcode{ = .true.}. 615 616 The default option \np{ln\_botmix\_triad}\forcode{ = .false.} is suitable when the bbl mixing option is enabled 617 (\key{trabbl}, with \np{nn\_bbl\_ldf}\forcode{ = 1}), or for simple idealized problems. 618 For setups with topography without bbl mixing, \np{ln\_botmix\_triad}\forcode{ = .true.} may be necessary. 666 619 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 667 620 \begin{figure}[h] \begin{center} 668 621 \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads} 669 622 \caption{ \protect\label{fig:bdry_triads} 670 (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer 671 points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad 672 slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ 673 (blue) poking through the ocean surface are masked (faded in 674 figure). However, the lateral $_{11}$ contributions towards 675 $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$ 676 (yellow line) are still applied, giving diapycnal diffusive 677 fluxes.\newline 623 (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots), 624 and $i+1/2,1$ $u$-point (blue square). 625 Triad slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) poking through 626 the ocean surface are masked (faded in figure). 627 However, the lateral $_{11}$ contributions towards $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and 628 $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$ (yellow line) are still applied, 629 giving diapycnal diffusive fluxes.\newline 678 630 (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 679 $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ 680 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point681 is masked. The associated lateral fluxes (grey-black dashed682 line) are masked if \protect\np{botmix\_triad}\forcode{ = .false.}, but left683 unmasked,giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}}631 $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, 632 i.e.\ the $i,k+1$ $u$-point is masked. 633 The associated lateral fluxes (grey-black dashed line) are masked if 634 \protect\np{botmix\_triad}\forcode{ = .false.}, but left unmasked, 635 giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}} 684 636 \end{center} \end{figure} 685 637 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 686 638 687 639 \subsection{ Limiting of the slopes within the interior}\label{sec:limit} 688 As discussed in \autoref{subsec:LDF_slp_iso}, iso-neutral slopes relative to689 geopotentials must be bounded everywhere, both for consistency with the small-slope 690 approximation and for numerical stability \citep{Cox1987, 691 Griffies_Bk04}. The bound chosen in \NEMO is applied to each 692 component of the slope separately andhas a value of $1/100$ in the ocean interior.640 As discussed in \autoref{subsec:LDF_slp_iso}, 641 iso-neutral slopes relative to geopotentials must be bounded everywhere, 642 both for consistency with the small-slope approximation and for numerical stability \citep{Cox1987, Griffies_Bk04}. 643 The bound chosen in \NEMO is applied to each component of the slope separately and 644 has a value of $1/100$ in the ocean interior. 693 645 %, ramping linearly down above 70~m depth to zero at the surface 694 It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to 695 geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to 696 geopotentials) \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate 697 surfaces, so we require 646 It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to geopotentials 647 (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to geopotentials) 648 \autoref{eq:PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate surfaces, so we require 698 649 \begin{equation*} 699 650 |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01. … … 701 652 and then recalculate the slopes $r_i$ relative to coordinates. 702 653 Each individual triad slope 703 \begin{equation} 704 \label{eq:Rtilde} 705 _i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p} + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} 706 \end{equation} 707 is limited like this and then the corresponding 708 $_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and combined to form the fluxes. 709 Note that where the slopes have been limited, there is now a non-zero 710 iso-neutral density flux that drives dianeutral mixing. In particular this iso-neutral density flux 711 is always downwards, and so acts to reduce gravitational potential energy. 654 \begin{equation} 655 \label{eq:Rtilde} 656 _i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p} + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} 657 \end{equation} 658 is limited like this and then the corresponding $_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and 659 combined to form the fluxes. 660 Note that where the slopes have been limited, there is now a non-zero iso-neutral density flux that 661 drives dianeutral mixing. 662 In particular this iso-neutral density flux is always downwards, 663 and so acts to reduce gravitational potential energy. 712 664 713 665 \subsection{Tapering within the surface mixed layer}\label{sec:taper} 714 Additional tapering of the iso-neutral fluxes is necessary within the 715 surface mixed layer. When the Griffies triads are used, we offer two 716 options for this. 666 Additional tapering of the iso-neutral fluxes is necessary within the surface mixed layer. 667 When the Griffies triads are used, we offer two options for this. 717 668 718 669 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:lintaper} 719 This is the option activated by the default choice 720 \np{ln\_triad\_iso}\forcode{ = .false.}. Slopes $\tilde{r}_i$ relative to 721 geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the 722 surface, as described in option (c) of \autoref{fig:eiv_slp}, to values 670 This is the option activated by the default choice \np{ln\_triad\_iso}\forcode{ = .false.}. 671 Slopes $\tilde{r}_i$ relative to geopotentials are tapered linearly from their value immediately below 672 the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:eiv_slp}, to values 723 673 \begin{subequations} 724 674 \begin{equation} 725 \label{eq:rmtilde}726 727 -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for } z>-h,675 \label{eq:rmtilde} 676 \rMLt = 677 -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for } z>-h, 728 678 \end{equation} 729 and then the $r_i$ relative to vertical coordinate surfaces are appropriately 730 adjusted to 679 and then the $r_i$ relative to vertical coordinate surfaces are appropriately adjusted to 731 680 \begin{equation} 732 \label{eq:rm}733 \rML =\rMLt -\sigma_i \quad \text{ for } z>-h.681 \label{eq:rm} 682 \rML =\rMLt -\sigma_i \quad \text{ for } z>-h. 734 683 \end{equation} 735 684 \end{subequations} … … 744 693 \end{equation} 745 694 746 This slope tapering gives a natural connection between tracer in the 747 mixed-layer and in isopycnal layers immediately below, in the 748 thermocline. It is consistent with the way the $\tilde{r}_i$ are 749 tapered within the mixed layer (see \autoref{sec:taperskew} below) 750 so as to ensure a uniform GM eddy-induced velocity throughout the 751 mixed layer. However, it gives a downwards density flux and so acts so 752 as to reduce potential energy in the same way as does the slope 753 limiting discussed above in \autoref{sec:limit}. 695 This slope tapering gives a natural connection between tracer in the mixed-layer and 696 in isopycnal layers immediately below, in the thermocline. 697 It is consistent with the way the $\tilde{r}_i$ are tapered within the mixed layer 698 (see \autoref{sec:taperskew} below) so as to ensure a uniform GM eddy-induced velocity throughout the mixed layer. 699 However, it gives a downwards density flux and so acts so as to reduce potential energy in the same way as 700 does the slope limiting discussed above in \autoref{sec:limit}. 754 701 755 As in \autoref{sec:limit} above, the tapering 756 \autoref{eq:rmtilde} is applied separately to each triad 757 $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the 758 $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume 759 $z$-coordinates in the following; the conversion from 760 $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described 761 above by \autoref{eq:Rtilde}. 702 As in \autoref{sec:limit} above, the tapering \autoref{eq:rmtilde} is applied separately to 703 each triad $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. 704 For clarity, we assume $z$-coordinates in the following; 705 the conversion from $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as 706 described above by \autoref{eq:Rtilde}. 762 707 \begin{enumerate} 763 \item Mixed-layer depth is defined so as to avoid including regions of weak 764 vertical stratification in the slope definition. 765 At each $i,j$ (simplified to $i$ in 766 \autoref{fig:MLB_triad}), we define the mixed-layer by setting 767 the vertical index of the tracer point immediately below the mixed 768 layer, $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) 769 such that the potential density 770 ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is 771 the tracer gridbox within which the depth reaches 10~m. See the left 772 side of \autoref{fig:MLB_triad}. We use the $k_{10}$-gridbox 773 instead of the surface gridbox to avoid problems e.g.\ with thin 774 daytime mixed-layers. Currently we use the same 775 $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is 776 used to output the diagnosed mixed-layer depth 777 $h_{\mathrm{ML}}=|z_{W}|_{k_{\mathrm{ML}}+1/2}$, the depth of the $w$-point 778 above the $i,k_{\mathrm{ML}}$ tracer point. 779 780 \item We define `basal' triad slopes 781 ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ as the slopes 782 of those triads whose vertical `arms' go down from the 783 $i,k_{\mathrm{ML}}$ tracer point to the $i,k_{\mathrm{ML}}-1$ tracer point 784 below. This is to ensure that the vertical density gradients 785 associated with these basal triad slopes 786 ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are 787 representative of the thermocline. The four basal triads defined in the bottom part 788 of \autoref{fig:MLB_triad} are then 708 \item 709 Mixed-layer depth is defined so as to avoid including regions of weak vertical stratification in 710 the slope definition. 711 At each $i,j$ (simplified to $i$ in \autoref{fig:MLB_triad}), 712 we define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, 713 $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) such that 714 the potential density ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, 715 where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10~m. 716 See the left side of \autoref{fig:MLB_triad}. 717 We use the $k_{10}$-gridbox instead of the surface gridbox to avoid problems e.g.\ with thin daytime mixed-layers. 718 Currently we use the same $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is used to 719 output the diagnosed mixed-layer depth $h_{\mathrm{ML}}=|z_{W}|_{k_{\mathrm{ML}}+1/2}$, 720 the depth of the $w$-point above the $i,k_{\mathrm{ML}}$ tracer point. 721 \item 722 We define `basal' triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ as 723 the slopes of those triads whose vertical `arms' go down from the $i,k_{\mathrm{ML}}$ tracer point to 724 the $i,k_{\mathrm{ML}}-1$ tracer point below. 725 This is to ensure that the vertical density gradients associated with 726 these basal triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are representative of the thermocline. 727 The four basal triads defined in the bottom part of \autoref{fig:MLB_triad} are then 789 728 \begin{align} 790 729 {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &= … … 795 734 {\:}^{k_{\mathrm{ML}}}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}. \notag 796 735 \end{align} 797 The vertical flux associated with each of these triads passes through the $w$-point 798 $i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point, 799 so it is this depth 736 The vertical flux associated with each of these triads passes through 737 the $w$-point $i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point, so it is this depth 800 738 \begin{equation} 801 739 \label{eq:zbase} 802 740 {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} 803 741 \end{equation} 804 (one gridbox deeper than the 805 diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper 806 the slopes in \autoref{eq:rmtilde}. 807 \item Finally, we calculate the adjusted triads 808 ${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within the mixed 809 layer, by multiplying the appropriate 810 ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ by the ratio of 811 the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$. For 812 instance the green triad centred on $i,k$ 742 one gridbox deeper than the diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper the slopes in 743 \autoref{eq:rmtilde}. 744 \item 745 Finally, we calculate the adjusted triads ${\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p}$ within 746 the mixed layer, by multiplying the appropriate ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ by 747 the ratio of the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$. 748 For instance the green triad centred on $i,k$ 813 749 \begin{align} 814 750 {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,1/2}^{-1/2} &= … … 824 760 \begin{figure}[h] 825 761 % \fcapside { 826 \caption{\protect\label{fig:MLB_triad} Definition of 827 mixed-layer depth and calculation of linearly tapered 828 triads. The figure shows a water column at a given $i,j$ 829 (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by 830 bullets, and black lines the edges of the tracer cells; $k$ 831 increases upwards. \newline 832 \hspace{5 em}We define the mixed-layer by setting the vertical index 833 of the tracer point immediately below the mixed layer, 834 $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) 835 such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, 836 where $i,k_{10}$ is the tracer gridbox within which the depth 837 reaches 10~m. We calculate the triad slopes within the mixed 838 layer by linearly tapering them from zero (at the surface) to 839 the `basal' slopes, the slopes of the four triads passing through the 840 $w$-point $i,k_{\mathrm{ML}}-1/2$ (blue square), 841 ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$. Triads with 842 different $i_p,k_p$, denoted by different colours, (e.g. the green 843 triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.} 762 \caption{\protect\label{fig:MLB_triad} 763 Definition of mixed-layer depth and calculation of linearly tapered triads. 764 The figure shows a water column at a given $i,j$ (simplified to $i$), with the ocean surface at the top. 765 Tracer points are denoted by bullets, and black lines the edges of the tracer cells; 766 $k$ increases upwards. \newline 767 \hspace{5 em} 768 We define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, 769 $k_{\mathrm{ML}}$, as the maximum $k$ (shallowest tracer point) such that 770 ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, 771 where $i,k_{10}$ is the tracer gridbox within which the depth reaches 10~m. 772 We calculate the triad slopes within the mixed layer by linearly tapering them from zero 773 (at the surface) to the `basal' slopes, 774 the slopes of the four triads passing through the $w$-point $i,k_{\mathrm{ML}}-1/2$ (blue square), 775 ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$. 776 Triads with different $i_p,k_p$, denoted by different colours, 777 (e.g. the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.} 844 778 %} 845 779 {\includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads}} … … 849 783 \subsubsection{Additional truncation of skew iso-neutral flux components} 850 784 \label{subsec:Gerdes-taper} 851 The alternative option is activated by setting \np{ln\_triad\_iso} = 852 true. This retains the same tapered slope $\rML$ described above for the 853 calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the 854 vertical tracer flux driven by vertical tracer gradients), but 855 replaces the $\rML$ in the skew term by 785 The alternative option is activated by setting \np{ln\_triad\_iso} = true. 786 This retains the same tapered slope $\rML$ described above for the calculation of the $_{33}$ term of 787 the iso-neutral diffusion tensor (the vertical tracer flux driven by vertical tracer gradients), 788 but replaces the $\rML$ in the skew term by 856 789 \begin{equation} 857 790 \label{eq:rm*} … … 868 801 \end{equation} 869 802 This operator 870 \footnote{ To ensure good behaviour where horizontal density871 gradients are weak, we in fact follow \citet{Gerdes1991} and set872 $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.} 873 then has the property it gives no vertical density flux, and so does 874 not change the potential energy.875 This approach is similar to multiplying the iso-neutral diffusion876 coefficient by $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep 877 slopes,as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}).803 \footnote{ 804 To ensure good behaviour where horizontal density gradients are weak, 805 we in fact follow \citet{Gerdes1991} and 806 set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.} 807 then has the property it gives no vertical density flux, and so does not change the potential energy. 808 This approach is similar to multiplying the iso-neutral diffusion coefficient by 809 $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep slopes, 810 as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}). 878 811 Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$ 879 812 880 In practice, this approach gives weak vertical tracer fluxes through 881 the mixed-layer, as well as vanishing density fluxes. While it is 882 theoretically advantageous that it does not change the potential 883 energy, it may give a discontinuity between the 884 fluxes within the mixed-layer (purely horizontal) and just below (along 885 iso-neutral surfaces). 813 In practice, this approach gives weak vertical tracer fluxes through the mixed-layer, 814 as well as vanishing density fluxes. 815 While it is theoretically advantageous that it does not change the potential energy, 816 it may give a discontinuity between the fluxes within the mixed-layer (purely horizontal) and 817 just below (along iso-neutral surfaces). 886 818 % This may give strange looking results, 887 819 % particularly where the mixed-layer depth varies strongly laterally. … … 893 825 \subsection{Continuous skew flux formulation}\label{sec:continuous-skew-flux} 894 826 895 When Gent and McWilliams's [1990] diffusion is used, 896 an additional advection term is added. The associated velocity is the so called 897 eddy induced velocity, the formulation of which depends on the slopes of iso- 898 neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used 899 here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo} 900 is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} 901 + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 827 When Gent and McWilliams's [1990] diffusion is used, an additional advection term is added. 828 The associated velocity is the so called eddy induced velocity, 829 the formulation of which depends on the slopes of iso-neutral surfaces. 830 Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, 831 $i.e.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinate, 832 and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates. 902 833 903 834 The eddy induced velocity is given by: … … 919 850 \end{equation} 920 851 \end{subequations} 921 with $A_{e}$ the eddy induced velocity coefficient, and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces. 922 923 The traditional way to implement this additional advection is to add 924 it to the Eulerian velocity prior to computing the tracer 925 advection. This is implemented if \key{traldf\_eiv} is set in the 926 default implementation, where \np{ln\_traldf\_triad} is set 927 false. This allows us to take advantage of all the advection schemes 928 offered for the tracers (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ 929 order advection scheme. This is particularly useful for passive 930 tracers where \emph{positivity} of the advection scheme is of 931 paramount importance. 932 933 However, when \np{ln\_traldf\_triad} is set true, \NEMO instead 934 implements eddy induced advection according to the so-called skew form 935 \citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes 936 using the non-divergent nature of the eddy induced velocity. 937 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective 938 fluxes per unit area in $ijk$ space can be 939 transformed as follows: 852 with $A_{e}$ the eddy induced velocity coefficient, 853 and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces. 854 855 The traditional way to implement this additional advection is to add it to the Eulerian velocity prior to 856 computing the tracer advection. 857 This is implemented if \key{traldf\_eiv} is set in the default implementation, 858 where \np{ln\_traldf\_triad} is set false. 859 This allows us to take advantage of all the advection schemes offered for the tracers 860 (see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme. 861 This is particularly useful for passive tracers where 862 \emph{positivity} of the advection scheme is of paramount importance. 863 864 However, when \np{ln\_traldf\_triad} is set true, 865 \NEMO instead implements eddy induced advection according to the so-called skew form \citep{Griffies_JPO98}. 866 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. 867 For example in the (\textbf{i},\textbf{k}) plane, 868 the tracer advective fluxes per unit area in $ijk$ space can be transformed as follows: 940 869 \begin{flalign*} 941 870 \begin{split} … … 962 891 \end{split} 963 892 \end{flalign*} 964 and since the eddy induced velocity field is non-divergent, we end up with the skew965 form of the eddy induced advective fluxes per unit area in $ijk$ space:893 and since the eddy induced velocity field is non-divergent, 894 we end up with the skew form of the eddy induced advective fluxes per unit area in $ijk$ space: 966 895 \begin{equation} \label{eq:eiv_skew_ijk} 967 896 \textbf{F}_\mathrm{eiv}^T = \begin{pmatrix} … … 979 908 \end{split} 980 909 \end{equation} 981 Note that \autoref{eq:eiv_skew_physical} takes the same form whatever the982 vertical coordinate, though of course the slopes 983 $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative togeopotentials.984 The tendency associated with eddy induced velocity is then simply the convergence 985 of the fluxes(\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so910 Note that \autoref{eq:eiv_skew_physical} takes the same form whatever the vertical coordinate, 911 though of course the slopes $\tilde{r}_i$ which define the $\psi_i$ in \autoref{eq:eiv_psi} are relative to 912 geopotentials. 913 The tendency associated with eddy induced velocity is then simply the convergence of the fluxes 914 (\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so 986 915 \begin{equation} \label{eq:skew_eiv_conv} 987 916 \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 } \left[ … … 992 921 + e_{1} \psi_2 \partial_j T \right) \right] 993 922 \end{equation} 994 It naturally conserves the tracer content, as it is expressed in flux 995 form. Since it has the same divergence as the advective form it also 996 preserves the tracer variance. 923 It naturally conserves the tracer content, as it is expressed in flux form. 924 Since it has the same divergence as the advective form it also preserves the tracer variance. 997 925 998 926 \subsection{Discrete skew flux formulation} 999 The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), like the off-diagonal terms1000 (\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, are best 1001 expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} 1002 and (\autoref{eq:i13}, \autoref{eq:i31}); but now in terms of the triad slopes 1003 $\tilde{\mathbb{R}}$ relative to geopotentials instead of the 1004 $\mathbb{R}$ relative to coordinate surfaces. The discrete form of 1005 \autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and927 The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), 928 like the off-diagonal terms (\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, 929 are best expressed in terms of the triad slopes, as in \autoref{fig:ISO_triad} and 930 (\autoref{eq:i13}, \autoref{eq:i31}); 931 but now in terms of the triad slopes $\tilde{\mathbb{R}}$ relative to geopotentials instead of 932 the $\mathbb{R}$ relative to coordinate surfaces. 933 The discrete form of \autoref{eq:eiv_skew_ijk} using the slopes \autoref{eq:R} and 1006 934 defining $A_e$ at $T$-points is then given by: 1007 935 … … 1017 945 \end{pmatrix}, 1018 946 \end{flalign} 1019 where the skew flux in the $i$-direction associated with a given 1020 triad is (\autoref{eq:latflux-triad},\autoref{eq:triadfluxu}):947 where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:latflux-triad}, 948 \autoref{eq:triadfluxu}): 1021 949 \begin{align} 1022 950 \label{eq:skewfluxu} … … 1034 962 \end{subequations} 1035 963 1036 Such a discretisation is consistent with the iso-neutral 1037 operator as it uses the same definition for the slopes. It also 1038 ensures the following two key properties. 964 Such a discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes. 965 It also ensures the following two key properties. 1039 966 1040 967 \subsubsection{No change in tracer variance} 1041 The discretization conserves tracer variance, $i.e.$ it does not 1042 include a diffusive component but is a `pure' advection term. This can 1043 be seen 1044 %either from Appendix \autoref{apdx:eiv_skew} or 1045 by considering the 1046 fluxes associated with a given triad slope 1047 $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following 1048 \autoref{subsec:variance} and \autoref{eq:dvar_iso_i}, the 1049 associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ 1050 drives a net rate of change of variance, summed over the two 1051 $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 968 The discretization conserves tracer variance, $i.e.$ it does not include a diffusive component but is a `pure' advection term. 969 This can be seen %either from Appendix \autoref{apdx:eiv_skew} or 970 by considering the fluxes associated with a given triad slope $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. 971 For, following \autoref{subsec:variance} and \autoref{eq:dvar_iso_i}, 972 the associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ drives a net rate of change of variance, 973 summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 1052 974 \begin{equation} 1053 975 \label{eq:dvar_eiv_i} 1054 976 _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], 1055 977 \end{equation} 1056 while the associated vertical skew-flux gives a variance change summed over the1057 $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of978 while the associated vertical skew-flux gives a variance change summed over 979 the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 1058 980 \begin{equation} 1059 981 \label{eq:dvar_eiv_k} 1060 982 _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 1061 983 \end{equation} 1062 Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) 1063 shows that these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) 1064 sum to zero. Hence the two fluxes associated with each triad make no 1065 net contribution to the variance budget. 984 Inspection of the definitions (\autoref{eq:skewfluxu}, \autoref{eq:skewfluxw}) shows that 985 these two variance changes (\autoref{eq:dvar_eiv_i}, \autoref{eq:dvar_eiv_k}) sum to zero. 986 Hence the two fluxes associated with each triad make no net contribution to the variance budget. 1066 987 1067 988 \subsubsection{Reduction in gravitational PE} 1068 The vertical density flux associated with the vertical skew-flux 1069 always has the same sign as the vertical density gradient; thus, so 1070 long as the fluid is stable (the vertical density gradient is 1071 negative) the vertical density flux is negative (downward) and hence 1072 reduces the gravitational PE. 989 The vertical density flux associated with the vertical skew-flux always has the same sign as 990 the vertical density gradient; 991 thus, so long as the fluid is stable (the vertical density gradient is negative) 992 the vertical density flux is negative (downward) and hence reduces the gravitational PE. 1073 993 1074 994 For the change in gravitational PE driven by the $k$-flux is … … 1091 1011 \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, 1092 1012 \end{align} 1093 using the definition of the triad slope $\rtriad{R}$, 1094 \autoref{eq:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 1095 \beta_i^k\delta_{i+ i_p}[S^k]$ in terms of $-\alpha_i^k \delta_{k+ 1096 k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. 1013 using the definition of the triad slope $\rtriad{R}$, \autoref{eq:R} to 1014 express $-\alpha _i^k\delta_{i+ i_p}[T^k]+\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of 1015 $-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. 1097 1016 1098 1017 Where the coordinates slope, the $i$-flux gives a PE change … … 1108 1027 \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, 1109 1028 \end{multline} 1110 (using \autoref{eq:skewfluxu}) and so the total PE change 1111 \autoref{eq: vert_densityPE} + \autoref{eq:lat_densityPE} associated with the triad fluxes is1029 (using \autoref{eq:skewfluxu}) and so the total PE change \autoref{eq:vert_densityPE} + 1030 \autoref{eq:lat_densityPE} associated with the triad fluxes is 1112 1031 \begin{multline} 1113 1032 \label{eq:tot_densityPE} … … 1122 1041 1123 1042 \subsection{Treatment of the triads at the boundaries}\label{sec:skew_bdry} 1124 Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes 1125 are masked at the boundaries in exactly the same way as are the triad 1126 slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as 1127 described in \autoref{sec:iso_bdry} and 1128 \autoref{fig:bdry_triads}. Thus surface layer triads 1129 $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are 1130 masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ 1131 and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the 1132 $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ 1133 $u$-point is masked. The namelist parameter \np{ln\_botmix\_triad} has 1134 no effect on the eddy-induced skew-fluxes. 1043 Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries 1044 in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes, 1045 as described in \autoref{sec:iso_bdry} and \autoref{fig:bdry_triads}. 1046 Thus surface layer triads $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are masked, 1047 and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when 1048 either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is masked. 1049 The namelist parameter \np{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes. 1135 1050 1136 1051 \subsection{Limiting of the slopes within the interior}\label{sec:limitskew} 1137 Presently, the iso-neutral slopes $\tilde{r}_i$ relative 1138 to geopotentials are limited to be less than $1/100$, exactly as in 1139 calculating the iso-neutral diffusion, \S \autoref{sec:limit}. Each 1140 individual triad \rtriadt{R} is so limited. 1052 Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$, 1053 exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:limit}. 1054 Each individual triad \rtriadt{R} is so limited. 1141 1055 1142 1056 \subsection{Tapering within the surface mixed layer}\label{sec:taperskew} 1143 The slopes $\tilde{r}_i$ relative to 1144 geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the 1145 surface \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}. This is 1146 option (c) of \autoref{fig:eiv_slp}. This linear tapering for the 1147 slopes used to calculate the eddy-induced fluxes is 1148 unaffected by the value of \np{ln\_triad\_iso}. 1149 1150 The justification for this linear slope tapering is that, for $A_e$ 1151 that is constant or varies only in the horizontal (the most commonly 1152 used options in \NEMO: see \autoref{sec:LDF_coef}), it is 1153 equivalent to a horizontal eiv (eddy-induced velocity) that is uniform 1154 within the mixed layer \autoref{eq:eiv_v}. This ensures that the 1155 eiv velocities do not restratify the mixed layer \citep{Treguier1997, 1156 Danabasoglu_al_2008}. Equivantly, in terms 1157 of the skew-flux formulation we use here, the 1158 linear slope tapering within the mixed-layer gives a linearly varying 1159 vertical flux, and so a tracer convergence uniform in depth (the 1160 horizontal flux convergence is relatively insignificant within the mixed-layer). 1057 The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R}) 1058 are always tapered linearly from their value immediately below the mixed layer to zero at the surface 1059 \autoref{eq:rmtilde}, as described in \autoref{sec:lintaper}. 1060 This is option (c) of \autoref{fig:eiv_slp}. 1061 This linear tapering for the slopes used to calculate the eddy-induced fluxes is unaffected by 1062 the value of \np{ln\_triad\_iso}. 1063 1064 The justification for this linear slope tapering is that, for $A_e$ that is constant or varies only in 1065 the horizontal (the most commonly used options in \NEMO: see \autoref{sec:LDF_coef}), 1066 it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer 1067 \autoref{eq:eiv_v}. 1068 This ensures that the eiv velocities do not restratify the mixed layer \citep{Treguier1997,Danabasoglu_al_2008}. 1069 Equivantly, in terms of the skew-flux formulation we use here, 1070 the linear slope tapering within the mixed-layer gives a linearly varying vertical flux, 1071 and so a tracer convergence uniform in depth 1072 (the horizontal flux convergence is relatively insignificant within the mixed-layer). 1161 1073 1162 1074 \subsection{Streamfunction diagnostics}\label{sec:sfdiag} 1163 Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, diagnosed 1164 mean eddy-induced velocities are output. Each time step, 1165 streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at 1166 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ 1167 (integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table 1168 \autoref{tab:cell}) respectively. We follow \citep{Griffies_Bk04} and 1169 calculate the streamfunction at a given $uw$-point from the 1170 surrounding four triads according to: 1075 Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, 1076 diagnosed mean eddy-induced velocities are output. 1077 Each time step, streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at 1078 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$) 1079 points (see Table \autoref{tab:cell}) respectively. 1080 We follow \citep{Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from 1081 the surrounding four triads according to: 1171 1082 \begin{equation} 1172 1083 \label{eq:sfdiagi} … … 1175 1086 \end{equation} 1176 1087 The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 1177 The eddy-induced velocities are then calculated from the 1178 straightforward discretisation of \autoref{eq:eiv_v}: 1088 The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:eiv_v}: 1179 1089 \begin{equation}\label{eq:eiv_v_discrete} 1180 1090 \begin{split}
Note: See TracChangeset
for help on using the changeset viewer.
