Changeset 3294 for trunk/DOC/TexFiles/Chapters/Annex_ISO.tex
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/DOC/TexFiles/Chapters/Annex_ISO.tex (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/DOC/TexFiles/Chapters/Annex_ISO.tex
r2376 r3294 1 1 % ================================================================ 2 % Iso-neutral diffusion : 2 % Iso-neutral diffusion : 3 3 % ================================================================ 4 \chapter{Griffies's iso-neutral diffusion} 5 \label{Apdx_C} 4 \chapter[Iso-neutral diffusion and eddy advection using 5 triads]{Iso-neutral diffusion and eddy advection using triads} 6 \label{sec:triad} 6 7 \minitoc 7 8 \section{Griffies's formulation of iso-neutral diffusion} 9 10 \subsection{Introduction} 11 We define a scheme that get its inspiration from the scheme of 12 \citet{Griffies_al_JPO98}, but is formulated within the \NEMO 8 \pagebreak 9 \section{Choice of namelist parameters} 10 %-----------------------------------------nam_traldf------------------------------------------------------ 11 \namdisplay{namtra_ldf} 12 %--------------------------------------------------------------------------------------------------------- 13 If the namelist variable \np{ln\_traldf\_grif} is set true (and 14 \key{ldfslp} is set), \NEMO updates both active and passive tracers 15 using the Griffies triad representation of iso-neutral diffusion and 16 the eddy-induced advective skew (GM) fluxes. Otherwise (by default) the 17 filtered version of Cox's original scheme is employed 18 (\S\ref{LDF_slp}). In the present implementation of the Griffies 19 scheme, the advective skew fluxes are implemented even if 20 \key{traldf\_eiv} is not set. 21 22 Values of iso-neutral diffusivity and GM coefficient are set as 23 described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd}, 24 N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and 25 GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and 26 \np{rn\_aeiv\_0}. If 2D-varying coefficients are set with 27 \key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal 28 scale factor according to \eqref{Eq_title} \footnote{Except in global 29 $0.5^{\circ}$ runs (\key{orca\_r05}) with \key{traldf\_eiv}, where 30 $A_l$ is set like $A_e$ but with a minimum vale of 31 $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$}. In idealised setups with 32 \key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv} 33 is set in the global configurations \key{orca\_r2}, \key{orca\_r1} or 34 \key{orca\_r05} with \key{traldf\_c2d}, a horizontally varying $A_e$ is 35 instead set from the Held-Larichev parameterisation\footnote{In this 36 case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further 37 reduced by a factor $|f/f_{20}|$, where $f_{20}$ is the value of $f$ 38 at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored 39 unless it is zero. 40 41 The options specific to the Griffies scheme include: 42 \begin{description}[font=\normalfont] 43 \item[\np{ln\_traldf\_gdia}] Default value is false. See \S\ref{sec:triad:sfdiag}. If this is set true, time-mean 44 eddy-advective (GM) velocities are output for diagnostic purposes, even 45 though the eddy advection is accomplished by means of the skew 46 fluxes. 47 \item[\np{ln\_traldf\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then 48 `iso-neutral' mixing is accomplished within the surface mixed-layer 49 along slopes linearly decreasing with depth from the value immediately below 50 the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}). This is the same 51 treatment as used in the default implementation 52 \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}. Where 53 \np{ln\_traldf\_iso} is set true, the vertical skew flux is further 54 reduced to ensure no vertical buoyancy flux, giving an almost pure 55 horizontal diffusive tracer flux within the mixed layer. This is similar to 56 the tapering suggested by \citet{Gerdes1991}. See \S\ref{sec:triad:Gerdes-taper} 57 \item[\np{ln\_traldf\_botmix}] See \S\ref{sec:triad:iso_bdry}. If this 58 is set false (the default) then the lateral diffusive fluxes 59 associated with triads partly masked by topography are neglected. If 60 it is set true, however, then these lateral diffusive fluxes are 61 applied, giving smoother bottom tracer fields at the cost of 62 introducing diapycnal mixing. 63 \end{description} 64 \section{Triad formulation of iso-neutral diffusion} 65 \label{sec:triad:iso} 66 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO 13 67 framework, using scale factors rather than grid-sizes. 14 68 15 The off-diagonal terms of the small angle diffusion tensor 16 \eqref{Eq_PE_iso_tensor} produce skew-fluxes along the 17 i- and j-directions resulting from the vertical tracer gradient: 18 \begin{align} 19 \label{eq:i13c} 20 &+\kappa r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad+\kappa r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\ 21 \intertext{and in the k-direction resulting from the lateral tracer gradients} 22 \label{eq:i31c} 23 & \kappa r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\kappa r_2\frac{1}{e_1}\frac{\partial T}{\partial i} 24 \end{align} 25 where \eqref{Eq_PE_iso_slopes} 69 \subsection{The iso-neutral diffusion operator} 70 The iso-neutral second order tracer diffusive operator for small 71 angles between iso-neutral surfaces and geopotentials is given by 72 \eqref{Eq_PE_iso_tensor}: 73 \begin{subequations} \label{eq:triad:PE_iso_tensor} 74 \begin{equation} 75 D^{lT}=-\Div\vect{f}^{lT}\equiv 76 -\frac{1}{e_1e_2e_3}\left[\pd{i}\left (f_1^{lT}e_2e_3\right) + 77 \pd{j}\left (f_2^{lT}e_2e_3\right) + \pd{k}\left (f_3^{lT}e_1e_2\right)\right], 78 \end{equation} 79 where the diffusive flux per unit area of physical space 80 \begin{equation} 81 \vect{f}^{lT}=-\Alt\Re\cdot\grad T, 82 \end{equation} 83 \begin{equation} 84 \label{eq:triad:PE_iso_tensor:c} 85 \mbox{with}\quad \;\;\Re = 86 \begin{pmatrix} 87 1&0&-r_1\mystrut \\ 88 0&1&-r_2\mystrut \\ 89 -r_1&-r_2&r_1 ^2+r_2 ^2\mystrut 90 \end{pmatrix} 91 \quad \text{and} \quad\grad T= 92 \begin{pmatrix} 93 \frac{1}{e_1}\pd[T]{i}\mystrut \\ 94 \frac{1}{e_2}\pd[T]{j}\mystrut \\ 95 \frac{1}{e_3}\pd[T]{k}\mystrut 96 \end{pmatrix}. 97 \end{equation} 98 \end{subequations} 99 % \left( {{\begin{array}{*{20}c} 100 % 1 \hfill & 0 \hfill & {-r_1 } \hfill \\ 101 % 0 \hfill & 1 \hfill & {-r_2 } \hfill \\ 102 % {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\ 103 % \end{array} }} \right) 104 Here \eqref{Eq_PE_iso_slopes} 26 105 \begin{align*} 27 106 r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} … … 33 112 }{\partial k} \right)^{-1} 34 113 \end{align*} 35 is the i-component of the slope of the isoneutral surface relative to the computational 36 surface, and $r_2$ is the j-component. 37 38 The extra vertical diffusive flux associated with the $_{33}$ 114 is the $i$-component of the slope of the iso-neutral surface relative to the computational 115 surface, and $r_2$ is the $j$-component. 116 117 We will find it useful to consider the fluxes per unit area in $i,j,k$ 118 space; we write 119 \begin{equation} 120 \label{eq:triad:Fijk} 121 \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). 122 \end{equation} 123 Additionally, we will sometimes write the contributions towards the 124 fluxes $\vect{f}$ and $\vect{F}_\mathrm{iso}$ from the component 125 $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, with 126 $f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc. 127 128 The off-diagonal terms of the small angle diffusion tensor 129 \eqref{Eq_PE_iso_tensor}, \eqref{eq:triad:PE_iso_tensor:c} produce skew-fluxes along the 130 $i$- and $j$-directions resulting from the vertical tracer gradient: 131 \begin{align} 132 \label{eq:triad:i13c} 133 f_{13}=&+\Alt r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+\Alt r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\ 134 \intertext{and in the k-direction resulting from the lateral tracer gradients} 135 \label{eq:triad:i31c} 136 f_{31}+f_{32}=& \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i} 137 \end{align} 138 139 The vertical diffusive flux associated with the $_{33}$ 39 140 component of the small angle diffusion tensor is 40 141 \begin{equation} 41 \label{eq: i33c}42 -\kappa(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.142 \label{eq:triad:i33c} 143 f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}. 43 144 \end{equation} 44 145 45 146 Since there are no cross terms involving $r_1$ and $r_2$ in the above, we can 46 consider the iso neutral diffusive fluxes separately in the i-k and j-k147 consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$ 47 148 planes, just adding together the vertical components from each 48 plane. The following description will describe the fluxes on the i-k149 plane. The following description will describe the fluxes on the $i$-$k$ 49 150 plane. 50 151 51 There is no natural discretization for the i-component of the52 skew-flux, \eqref{eq: i13c}, as53 although it must be evaluated at u-points, it involves vertical152 There is no natural discretization for the $i$-component of the 153 skew-flux, \eqref{eq:triad:i13c}, as 154 although it must be evaluated at $u$-points, it involves vertical 54 155 gradients (both for the tracer and the slope $r_1$), defined at 55 w-points. Similarly, the vertical skew flux, \eqref{eq:i31c}, is evaluated at56 w-points but involves horizontal gradients defined at u-points.156 $w$-points. Similarly, the vertical skew flux, \eqref{eq:triad:i31c}, is evaluated at 157 $w$-points but involves horizontal gradients defined at $u$-points. 57 158 58 159 \subsection{The standard discretization} 59 160 The straightforward approach to discretize the lateral skew flux 60 \eqref{eq: i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995161 \eqref{eq:triad:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 61 162 into OPA, \eqref{Eq_tra_ldf_iso}, is to calculate a mean vertical 62 gradient at the u-point from the average of the four surrounding163 gradient at the $u$-point from the average of the four surrounding 63 164 vertical tracer gradients, and multiply this by a mean slope at the 64 u-point, calculated from the averaged surrounding vertical density 65 gradients. The total area-integrated skew-flux from tracer cell $i,k$ 66 to $i+1,k$, noting that the $e_{{3u}_{i+1/2}^k}$ in the area 67 $e_{{3u}_{i+1/2}^k}e_{{2u}_{i+1/2}^k}$ at the u-point cancels out with 68 the $1/e_{{3u}_{i+1/2}^k}$ associated with the vertical tracer 165 $u$-point, calculated from the averaged surrounding vertical density 166 gradients. The total area-integrated skew-flux (flux per unit area in 167 $ijk$ space) from tracer cell $i,k$ 168 to $i+1,k$, noting that the $e_{{3}_{i+1/2}^k}$ in the area 169 $e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with 170 the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer 69 171 gradient, is then \eqref{Eq_tra_ldf_iso} 70 172 \begin{equation*} 71 \left(F_u^{ \mathrm{skew}} \right)_{i+\hhalf}^k = \kappa_{i+\hhalf}^k72 e_{{2u}_{i+1/2}^k}\overline{\overline173 \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k 174 {e_{2}}_{i+1/2}^k \overline{\overline 73 175 r_1} ^{\,i,k}\,\overline{\overline{\delta_k T}}^{\,i,k}, 74 176 \end{equation*} … … 76 178 \begin{equation*} 77 179 \overline{\overline 78 r_1} ^{\,i,k} = -\frac{e_{{3u}_{i+1/2}^k}}{e_{{1u}_{i+1/2}^k}}79 \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}} .180 r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k} 181 \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}}, 80 182 \end{equation*} 81 183 and here and in the following we drop the $^{lT}$ superscript from 184 $\Alt$ for simplicity. 82 185 Unfortunately the resulting combination $\overline{\overline{\delta_k 83 186 \bullet}}^{\,i,k}$ of a $k$ average and a $k$ difference %of the tracer … … 89 192 global-average variance. To correct this, we introduced a smoothing of 90 193 the slopes of the iso-neutral surfaces (see \S\ref{LDF}). This 91 technique works f ine for $T$ and $S$as they are active tracers194 technique works for $T$ and $S$ in so far as they are active tracers 92 195 ($i.e.$ they enter the computation of density), but it does not work 93 196 for a passive tracer. … … 95 198 \citep{Griffies_al_JPO98} introduce a different discretization of the 96 199 off-diagonal terms that nicely solves the problem. 97 % Instead of multiplying the mean slope calculated at the u-point by98 % the mean vertical gradient at the u-point,200 % Instead of multiplying the mean slope calculated at the $u$-point by 201 % the mean vertical gradient at the $u$-point, 99 202 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 100 203 \begin{figure}[h] \begin{center} 101 204 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_triad_fluxes} 102 \caption{ \label{Fig_ISO_triad}205 \caption{ \label{fig:triad:ISO_triad} 103 206 (a) Arrangement of triads $S_i$ and tracer gradients to 104 give lateral tracer flux from box $i,k$ to $i+1,k$ 207 give lateral tracer flux from box $i,k$ to $i+1,k$ 105 208 (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from 106 209 box $i,k$ to $i,k+1$.} 107 \label{Fig_triad} 108 \end{center} \end{figure} 210 \end{center} \end{figure} 109 211 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 110 212 They get the skew flux from the products of the vertical gradients at 111 each w-point surrounding the u-point with the corresponding `triad'112 slope calculated from the lateral density gradient across the u-point113 divided by the vertical density gradient at the same w-point as the114 tracer gradient. See Fig.~\ref{ Fig_triad}a, where the thick lines213 each $w$-point surrounding the $u$-point with the corresponding `triad' 214 slope calculated from the lateral density gradient across the $u$-point 215 divided by the vertical density gradient at the same $w$-point as the 216 tracer gradient. See Fig.~\ref{fig:triad:ISO_triad}a, where the thick lines 115 217 denote the tracer gradients, and the thin lines the corresponding 116 218 triads, with slopes $s_1, \dotsc s_4$. The total area-integrated 117 219 skew-flux from tracer cell $i,k$ to $i+1,k$ 118 220 \begin{multline} 119 \label{eq: i13}120 \left( F_u^{13} \right)_{i+\frac{1}{2}}^k = \ kappa_{i+1}^k a_1 s_1221 \label{eq:triad:i13} 222 \left( F_u^{13} \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1 121 223 \delta _{k+\frac{1}{2}} \left[ T^{i+1} 122 \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} + \ kappa_i^k a_2 s_2 \delta224 \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} + \Alts _i^k a_2 s_2 \delta 123 225 _{k+\frac{1}{2}} \left[ T^i 124 226 \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\ 125 +\ kappa_{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1}126 \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} +\ kappa_i^k a_4 s_4 \delta227 +\Alts _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1} 228 \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} +\Alts _i^k a_4 s_4 \delta 127 229 _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}, 128 230 \end{multline} 129 231 where the contributions of the triad fluxes are weighted by areas 130 $a_1, \dotsc a_4$, and $\ kappa$ is now defined at the tracer points131 rather than the u-points. This discretization gives a much closer232 $a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points 233 rather than the $u$-points. This discretization gives a much closer 132 234 stencil, and disallows the two-point computational modes. 133 235 134 The vertical skew flux \eqref{eq: i31c} from tracer cell $i,k$ to $i,k+1$ at the135 w-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{Fig_triad}b)236 The vertical skew flux \eqref{eq:triad:i31c} from tracer cell $i,k$ to $i,k+1$ at the 237 $w$-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{fig:triad:ISO_triad}b) 136 238 by multiplying lateral tracer gradients from each of the four 137 surrounding u-points by the appropriate triad slope:239 surrounding $u$-points by the appropriate triad slope: 138 240 \begin{multline} 139 \label{eq: i31}140 \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} = \ kappa_i^{k+1} a_{1}'241 \label{eq:triad:i31} 242 \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} = \Alts_i^{k+1} a_{1}' 141 243 s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} 142 +\ kappa_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\143 + \ kappa_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k144 +\ kappa_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k.244 +\Alts_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\ 245 + \Alts_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k 246 +\Alts_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k. 145 247 \end{multline} 146 147 We notate the triad slopes in terms of the `anchor point' $i,k$148 (appearing in both the vertical and lateral gradient), and the u- and149 w-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the150 triad as follows (see also Fig.~\ref{ Fig_triad}):151 \begin{equation} 152 \label{ Gf_slopes}153 _i^k \mathbb{R}_{i_p}^{k_p} 154 = \frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}}248 249 We notate the triad slopes $s_i$ and $s'_i$ in terms of the `anchor point' $i,k$ 250 (appearing in both the vertical and lateral gradient), and the $u$- and 251 $w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the 252 triad as follows (see also Fig.~\ref{fig:triad:ISO_triad}): 253 \begin{equation} 254 \label{eq:triad:R} 255 _i^k \mathbb{R}_{i_p}^{k_p} 256 =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}} 155 257 \ 156 \frac 258 \frac 157 259 {\left(\alpha / \beta \right)_i^k \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] } 158 260 {\left(\alpha / \beta \right)_i^k \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }. … … 160 262 In calculating the slopes of the local neutral 161 263 surfaces, the expansion coefficients $\alpha$ and $\beta$ are 162 evaluated at the anchor points of the triad \footnote{Note that in \eqref{ Gf_slopes} we use the ratio $\alpha / \beta$264 evaluated at the anchor points of the triad \footnote{Note that in \eqref{eq:triad:R} we use the ratio $\alpha / \beta$ 163 265 instead of multiplying the temperature derivative by $\alpha$ and the 164 266 salinity derivative by $\beta$. This is more efficient as the ratio 165 267 $\alpha / \beta$ can to be evaluated directly}, while the metrics are 166 calculated at the u- and w-points on the arms.268 calculated at the $u$- and $w$-points on the arms. 167 269 168 270 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 169 271 \begin{figure}[h] \begin{center} 170 272 \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_qcells} 171 \caption{ \label{Fig_ISO_triad_notation} 172 Triad notation for quarter cells.T-cells are inside 173 boxes, while the $i+\half,k$ u-cell is shaded in green and the 174 $i,k+\half$ w-cell is shaded in pink.} 175 \label{qcells} 273 \caption{ \label{fig:triad:qcells} 274 Triad notation for quarter cells.$T$-cells are inside 275 boxes, while the $i+\half,k$ $u$-cell is shaded in green and the 276 $i,k+\half$ $w$-cell is shaded in pink.} 176 277 \end{center} \end{figure} 177 278 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 178 279 179 Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{ qcells}) with the quarter180 cell that is the intersection of the $i,k$ T-cell, the $i+i_p,k$181 u-cell and the $i,k+k_p$ w-cell. Expressing the slopes $s_i$ and182 $s'_i$ in \eqref{eq: i13} and \eqref{eq:i31} in this notation, we have280 Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter 281 cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$ 282 $u$-cell and the $i,k+k_p$ $w$-cell. Expressing the slopes $s_i$ and 283 $s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation, we have 183 284 e.g.\ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. Each triad slope $_i^k 184 285 \mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to calculate the 185 lateral flux along its u-arm, at $(i+i_p,k)$, and then again as an186 $s'$ to calculate the vertical flux along its w-arm at286 lateral flux along its $u$-arm, at $(i+i_p,k)$, and then again as an 287 $s'$ to calculate the vertical flux along its $w$-arm at 187 288 $(i,k+k_p)$. Each vertical area $a_i$ used to calculate the lateral 188 289 flux and horizontal area $a'_i$ used to calculate the vertical flux 189 can also be identified as the area across the u- and w-arms of a190 unique triad, and we cannotate these areas, similarly to the triad290 can also be identified as the area across the $u$- and $w$-arms of a 291 unique triad, and we notate these areas, similarly to the triad 191 292 slopes, as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$, 192 $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq: i13}193 $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq: i31}293 $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:triad:i13} 294 $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:triad:i31} 194 295 $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$. 195 296 196 297 \subsection{The full triad fluxes} 197 A key property of iso neutral diffusion is that it should not affect298 A key property of iso-neutral diffusion is that it should not affect 198 299 the (locally referenced) density. In particular there should be no 199 300 lateral or vertical density flux. The lateral density flux disappears so long as the … … 202 303 form 203 304 \begin{equation} 204 \label{eq: i11}305 \label{eq:triad:i11} 205 306 \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} = 206 - \left( \ kappa_i^{k+1} a_{1} + \kappa_i^{k+1} a_{2} + \kappa_i^k207 a_{3} + \ kappa_i^k a_{4} \right)307 - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k 308 a_{3} + \Alts_i^k a_{4} \right) 208 309 \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}}, 209 310 \end{equation} 210 where the areas $a_i$ are as in \eqref{eq: i13}. In this case,211 separating the total lateral flux, the sum of \eqref{eq: i13} and212 \eqref{eq: i11}, into triad components, a lateral tracer311 where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case, 312 separating the total lateral flux, the sum of \eqref{eq:triad:i13} and 313 \eqref{eq:triad:i11}, into triad components, a lateral tracer 213 314 flux 214 315 \begin{equation} 215 \label{eq: latflux-triad}216 _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - A_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p}316 \label{eq:triad:latflux-triad} 317 _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p} 217 318 \left( 218 319 \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } … … 227 328 density flux associated with each triad separately disappears. 228 329 \begin{equation} 229 \label{eq: latflux-rho}330 \label{eq:triad:latflux-rho} 230 331 {\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 231 332 \end{equation} … … 234 335 to $i+1,k$ must also vanish since it is a sum of four such triad fluxes. 235 336 236 The squared slope $r_1^2$ in the expression \eqref{eq: i33c} for the337 The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the 237 338 $_{33}$ component is also expressed in terms of area-weighted 238 339 squared triad slopes, so the area-integrated vertical flux from tracer 239 340 cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is 240 341 \begin{equation} 241 \label{eq: i33}342 \label{eq:triad:i33} 242 343 \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = 243 - \left( \ kappa_i^{k+1} a_{1}' s_{1}'^2244 + \ kappa_i^{k+1} a_{2}' s_{2}'^2245 + \ kappa_i^k a_{3}' s_{3}'^2246 + \ kappa_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],344 - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2 345 + \Alts_i^{k+1} a_{2}' s_{2}'^2 346 + \Alts_i^k a_{3}' s_{3}'^2 347 + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right], 247 348 \end{equation} 248 349 where the areas $a'$ and slopes $s'$ are the same as in 249 \eqref{eq: i31}.250 Then, separating the total vertical flux, the sum of \eqref{eq: i31} and251 \eqref{eq: i33}, into triad components, a vertical flux350 \eqref{eq:triad:i31}. 351 Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and 352 \eqref{eq:triad:i33}, into triad components, a vertical flux 252 353 \begin{align} 253 \label{eq: vertflux-triad}354 \label{eq:triad:vertflux-triad} 254 355 _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 255 &= A_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p}356 &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} 256 357 \left( 257 358 {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } … … 260 361 \right) \\ 261 362 &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) 262 {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq: vertflux-triad2}363 {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:triad:vertflux-triad2} 263 364 \end{align} 264 365 may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$ … … 270 371 fluxes. 271 372 272 We can explicitly identify (Fig.~\ref{ qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of273 the u-fluxes and w-fluxes in274 \eqref{eq: i31}, \eqref{eq:i13}, \eqref{eq:i11} \eqref{eq:i33} and275 Fig.~\ref{ Fig_triad} to write out the iso-neutral fluxes at $u$- and373 We can explicitly identify (Fig.~\ref{fig:triad:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of 374 the $u$-fluxes and $w$-fluxes in 375 \eqref{eq:triad:i31}, \eqref{eq:triad:i13}, \eqref{eq:triad:i11} \eqref{eq:triad:i33} and 376 Fig.~\ref{fig:triad:ISO_triad} to write out the iso-neutral fluxes at $u$- and 276 377 $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 277 378 %(Fig.~\ref{Fig_ISO_triad}): 278 \begin{flalign} \label{Eq_iso_flux} \ textbf{F}_{iso}(T) &\equiv379 \begin{flalign} \label{Eq_iso_flux} \vect{F}_\mathrm{iso}(T) &\equiv 279 380 \sum_{\substack{i_p,\,k_p}} 280 381 \begin{pmatrix} … … 284 385 \end{pmatrix}. 285 386 \end{flalign} 286 \subsection{Ensuring the scheme cannot increase tracer variance}287 \label{sec: variance}288 289 We now require that this operator cannot increase the387 \subsection{Ensuring the scheme does not increase tracer variance} 388 \label{sec:triad:variance} 389 390 We now require that this operator should not increase the 290 391 globally-integrated tracer variance. 291 392 %This changes according to 292 393 % \begin{align*} 293 394 % &\int_D D_l^T \; T \;dv \equiv \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\} \\ 294 % &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 295 % \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] 395 % &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 396 % \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right] 296 397 % + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \ T \right\} \\ 297 % &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 398 % &\equiv - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{ 298 399 % {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T] 299 400 % + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \ \delta_{k+1/2} [T] \right\} \\ 300 401 % \end{align*} 301 402 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux 302 $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the u-point $i+i_p,k$ and403 $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and 303 404 a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the 304 w-point $i,k+k_p$. The lateral flux drives a net rate of change of305 variance at points $i+i_p-\half,k$ and $i+i_p+\half,k$of405 $w$-point $i,k+k_p$. The lateral flux drives a net rate of change of 406 variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 306 407 \begin{multline} 307 408 {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+ … … 311 412 &= -T_{i+i_p-1/2}^k{\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \quad + \quad T_{i+i_p+1/2}^k 312 413 {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ 313 &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{ locFdtdx}414 &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:triad:dvar_iso_i} 314 415 \end{split} 315 416 \end{multline} 316 while the vertical flux similarly drives a net rate of change of variance at points $i,k+k_p-\half$ and 317 $i,k+k_p+\half$ of 318 \begin{equation} 319 \label{locFdtdz} 417 while the vertical flux similarly drives a net rate of change of 418 variance summed over the $T$-points $i,k+k_p-\half$ (above) and 419 $i,k+k_p+\half$ (below) of 420 \begin{equation} 421 \label{eq:triad:dvar_iso_k} 320 422 _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 321 423 \end{equation} 322 424 The total variance tendency driven by the triad is the sum of these 323 425 two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and 324 $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq: latflux-triad} and325 \eqref{eq: vertflux-triad}, it is426 $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:triad:latflux-triad} and 427 \eqref{eq:triad:vertflux-triad}, it is 326 428 \begin{multline*} 327 - A_i^k\left \{429 -\Alts_i^k\left \{ 328 430 { } _i^k{\mathbb{A}_u}_{i_p}^{k_p} 329 431 \left( … … 343 445 to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by 344 446 \begin{equation} 345 \label{eq: V-A}447 \label{eq:triad:V-A} 346 448 _i^k\mathbb{V}_{i_p}^{k_p} 347 449 ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k} … … 350 452 the variance tendency reduces to the perfect square 351 453 \begin{equation} 352 \label{eq: perfect-square}353 - A_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}454 \label{eq:triad:perfect-square} 455 -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 354 456 \left( 355 457 \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } … … 358 460 \right)^2\leq 0. 359 461 \end{equation} 360 Thus, the constraint \eqref{eq: V-A} ensures that the fluxesassociated462 Thus, the constraint \eqref{eq:triad:V-A} ensures that the fluxes (\ref{eq:triad:latflux-triad}, \ref{eq:triad:vertflux-triad}) associated 361 463 with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase 362 464 the net variance. Since the total fluxes are sums of such fluxes from … … 365 467 increase. 366 468 367 The expression \eqref{eq: V-A} can be interpreted as a discretization469 The expression \eqref{eq:triad:V-A} can be interpreted as a discretization 368 470 of the global integral 369 471 \begin{equation} 370 \label{eq: cts-var}371 \frac{\partial}{\partial t}\int\ half T^2\, dV =372 \int\ mathbf{F}\cdot\nabla T\, dV,472 \label{eq:triad:cts-var} 473 \frac{\partial}{\partial t}\int\!\half T^2\, dV = 474 \int\!\mathbf{F}\cdot\nabla T\, dV, 373 475 \end{equation} 374 476 where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the … … 376 478 \[ 377 479 \mathbf{F}=\left( 378 _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)/_i^k{\mathbb{A}_u}_{i_p}^{k_p},379 {\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)/_i^k{\mathbb{A}_w}_{i_p}^{k_p}480 \left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p}, 481 \left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p} 380 482 \right) 381 483 \] 382 484 and the gradient 383 485 \[\nabla T = \left( 384 \ delta_{i+ i_p}[T^k]/ {e_{1u}}_{\,i + i_p}^{\,k},385 \ delta_{k+ k_p}[T^i]/ {e_{3w}}_{\,i}^{\,k + k_p}486 \left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k}, 487 \left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p} 386 488 \right) 387 489 \] … … 390 492 volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify 391 493 these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter 392 cells, defined in terms of the distances between T, u,fand393 w-points. This is the natural discretization of394 \eqref{eq: cts-var}. The \NEMO model, however, operates with scale494 cells, defined in terms of the distances between $T$, $u$,$f$ and 495 $w$-points. This is the natural discretization of 496 \eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale 395 497 factors instead of grid sizes, and scale factors for the quarter 396 498 cells are not defined. Instead, therefore we simply choose 397 499 \begin{equation} 398 \label{eq: V-NEMO}500 \label{eq:triad:V-NEMO} 399 501 _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k, 400 502 \end{equation} 401 as a quarter of the volume of the u-cell inside which the triad503 as a quarter of the volume of the $u$-cell inside which the triad 402 504 quarter-cell lies. This has the nice property that when the slopes 403 505 $\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to 404 506 $i+1,k$ reduces to the classical form 405 507 \begin{equation} 406 \label{eq: lat-normal}407 -\overline {A}_{\,i+1/2}^k\;508 \label{eq:triad:lat-normal} 509 -\overline\Alts_{\,i+1/2}^k\; 408 510 \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 409 511 \;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}} 410 = -\overline {A}_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}}{{e_{1u}}_{\,i + 1/2}^{\,k}}.411 \end{equation} 412 In fact if the diffusive coefficient is defined at u-points, so that413 we employ $ A_{i+i_p}^k$ instead of $A_i^k$ in the definitions of the414 triad fluxes \eqref{eq: latflux-triad} and \eqref{eq:vertflux-triad},512 = -\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}}. 513 \end{equation} 514 In fact if the diffusive coefficient is defined at $u$-points, so that 515 we employ $\Alts_{i+i_p}^k$ instead of $\Alts_i^k$ in the definitions of the 516 triad fluxes \eqref{eq:triad:latflux-triad} and \eqref{eq:triad:vertflux-triad}, 415 517 we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above. 416 518 417 519 \subsection{Summary of the scheme} 418 The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 520 The iso-neutral fluxes at $u$- and 521 $w$-points are the sums of the triad fluxes that cross the $u$- and 522 $w$-faces \eqref{Eq_iso_flux}: 523 \begin{subequations}\label{eq:triad:alltriadflux} 524 \begin{flalign}\label{eq:triad:vect_isoflux} 525 \vect{F}_\mathrm{iso}(T) &\equiv 526 \sum_{\substack{i_p,\,k_p}} 527 \begin{pmatrix} 528 {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ 529 \\ 530 {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) 531 \end{pmatrix}, 532 \end{flalign} 533 where \eqref{eq:triad:latflux-triad}: 534 \begin{align} 535 \label{eq:triad:triadfluxu} 536 _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{ 537 \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} 538 \left( 539 \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 540 -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \ 541 \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 542 \right),\\ 543 \intertext{and} 544 _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 545 &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}} 546 \left( 547 {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 548 -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 549 \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 550 \right),\label{eq:triad:triadfluxw} 551 \end{align} 552 with \eqref{eq:triad:V-NEMO} 553 \begin{equation} 554 \label{eq:triad:V-NEMO2} 555 _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k. 556 \end{equation} 557 \end{subequations} 558 559 The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 419 560 each tracer point: 420 \begin{equation} \label{ Gf_operator} D_l^T = \frac{1}{b_T}561 \begin{equation} \label{eq:triad:iso_operator} D_l^T = \frac{1}{b_T} 421 562 \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k 422 563 {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[ … … 427 568 \begin{description} 428 569 \item[$\bullet$ horizontal diffusion] The discretization of the 429 diffusion operator recovers \eqref{eq: lat-normal} the traditional five-point Laplacian in570 diffusion operator recovers \eqref{eq:triad:lat-normal} the traditional five-point Laplacian in 430 571 the limit of flat iso-neutral direction : 431 \begin{equation} \label{ Gf_property1a} D_l^T = \frac{1}{b_T} \572 \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \ 432 573 \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; 433 \overline {A}^{\,i} \; \delta_{i+1/2}[T] \right] \qquad574 \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad 434 575 \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0 435 576 \end{equation} … … 437 578 \item[$\bullet$ implicit treatment in the vertical] Only tracer values 438 579 associated with a single water column appear in the expression 439 \eqref{eq: i33} for the $_{33}$ fluxes, vertical fluxes driven by580 \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by 440 581 vertical gradients. This is of paramount importance since it means 441 that an implicit in time algorithm can be used to solve the vertical 442 diffusion equation. This is a necessity since the vertical eddy 582 that a time-implicit algorithm can be used to solve the vertical 583 diffusion equation. This is necessary 584 since the vertical eddy 443 585 diffusivity associated with this term, 444 586 \begin{equation} 445 \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{ 446 {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: A_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2447 \right\} = 448 \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{ 449 {b_u}_{i+i_p}^k\: A_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2587 \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{ 588 {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 589 \right\} = 590 \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{ 591 {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 450 592 \right\}, 451 593 \end{equation} … … 454 596 \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of 455 597 locally referenced potential density is zero. See 456 \eqref{eq: latflux-rho} and \eqref{eq:vertflux-triad2}.598 \eqref{eq:triad:latflux-rho} and \eqref{eq:triad:vertflux-triad2}. 457 599 458 600 \item[$\bullet$ conservation of tracer] The iso-neutral diffusion 459 601 conserves tracer content, $i.e.$ 460 \begin{equation} \label{ Gf_property1} \sum_{i,j,k} \left\{ D_l^T \602 \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 461 603 b_T \right\} = 0 462 604 \end{equation} … … 466 608 \item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion 467 609 does not increase the tracer variance, $i.e.$ 468 \begin{equation} \label{ Gf_property1} \sum_{i,j,k} \left\{ T \ D_l^T610 \begin{equation} \label{eq:triad:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 469 611 \ b_T \right\} \leq 0 470 612 \end{equation} 471 613 The property is demonstrated in 472 \ ref{sec:variance} above. It is a key property for a diffusion473 term. It means that it is also a dissipation term, $i.e.$ it is a474 di ffusion ofthe square of the quantity on which it is applied. It614 \S\ref{sec:triad:variance} above. It is a key property for a diffusion 615 term. It means that it is also a dissipation term, $i.e.$ it 616 dissipates the square of the quantity on which it is applied. It 475 617 therefore ensures that, when the diffusivity coefficient is large 476 enough, the field on which it is applied become free of grid-point618 enough, the field on which it is applied becomes free of grid-point 477 619 noise. 478 620 479 621 \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion 480 622 operator is self-adjoint, $i.e.$ 481 \begin{equation} \label{ Gf_property1} \sum_{i,j,k} \left\{ S \ D_l^T623 \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 482 624 \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 483 625 \end{equation} … … 486 628 routine. This property can be demonstrated similarly to the proof of 487 629 the `no increase of tracer variance' property. The contribution by a 488 single triad towards the left hand side of \eqref{ Gf_property1}, can489 be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{ locFdtdx}490 and \eqref{ locFdtdx}. This results in a term similar to491 \eqref{eq: perfect-square},492 \begin{equation} 493 \label{eq: TScovar}494 - A_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}630 single triad towards the left hand side of \eqref{eq:triad:iso_property3}, can 631 be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{eq:triad:dvar_iso_i} 632 and \eqref{eq:triad:dvar_iso_k}. This results in a term similar to 633 \eqref{eq:triad:perfect-square}, 634 \begin{equation} 635 \label{eq:triad:TScovar} 636 - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 495 637 \left( 496 638 \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } … … 506 648 This is symmetrical in $T $ and $S$, so exactly the same term arises 507 649 from the discretization of this triad's contribution towards the 508 RHS of \eqref{ Gf_property1}.650 RHS of \eqref{eq:triad:iso_property3}. 509 651 \end{description} 510 511 512 $\ $\newline %force an empty line 652 \subsection{Treatment of the triads at the boundaries}\label{sec:triad:iso_bdry} 653 The triad slope can only be defined where both the grid boxes centred at 654 the end of the arms exist. Triads that would poke up 655 through the upper ocean surface into the atmosphere, or down into the 656 ocean floor, must be masked out. See Fig.~\ref{fig:triad:bdry_triads}. Surface layer triads 657 $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and 658 $\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be 659 specified above the ocean surface are masked (Fig.~\ref{fig:triad:bdry_triads}a): this ensures that lateral 660 tracer gradients produce no flux through the ocean surface. However, 661 to prevent surface noise, it is customary to retain the $_{11}$ contributions towards 662 the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and 663 $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer 664 fluxes. Similar comments apply to triads that would intersect the 665 ocean floor (Fig.~\ref{fig:triad:bdry_triads}b). Note that both near bottom 666 triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 667 $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ 668 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is 669 masked. The associated lateral fluxes (grey-black dashed line) are 670 masked if \np{ln\_botmix\_grif}=false, but left unmasked, 671 giving bottom mixing, if \np{ln\_botmix\_grif}=true. 672 673 The default option \np{ln\_botmix\_grif}=false is suitable when the 674 bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1), 675 or for simple idealized problems. For setups with topography without 676 bbl mixing, \np{ln\_botmix\_grif}=true may be necessary. 677 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 678 \begin{figure}[h] \begin{center} 679 \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_bdry_triads} 680 \caption{ \label{fig:triad:bdry_triads} 681 (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer 682 points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad 683 slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$ 684 (blue) poking through the ocean surface are masked (faded in 685 figure). However, the lateral $_{11}$ contributions towards 686 $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$ 687 (yellow line) are still applied, giving diapycnal diffusive 688 fluxes.\\ 689 (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 690 $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ 691 or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point 692 is masked. The associated lateral fluxes (grey-black dashed 693 line) are masked if \np{botmix\_grif}=.false., but left 694 unmasked, giving bottom mixing, if \np{botmix\_grif}=.true.} 695 \end{center} \end{figure} 696 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 697 \subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit} 698 As discussed in \S\ref{LDF_slp_iso}, iso-neutral slopes relative to 699 geopotentials must be bounded everywhere, both for consistency with the small-slope 700 approximation and for numerical stability \citep{Cox1987, 701 Griffies_Bk04}. The bound chosen in \NEMO is applied to each 702 component of the slope separately and has a value of $1/100$ in the ocean interior. 703 %, ramping linearly down above 70~m depth to zero at the surface 704 It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to 705 geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to 706 geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate 707 surfaces, so we require 708 \begin{equation*} 709 |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01. 710 \end{equation*} 711 and then recalculate the slopes $r_i$ relative to coordinates. 712 Each individual triad slope 713 \begin{equation} 714 \label{eq:triad:Rtilde} 715 _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}} 716 \end{equation} 717 is limited like this and then the corresponding 718 $_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and combined to form the fluxes. 719 Note that where the slopes have been limited, there is now a non-zero 720 iso-neutral density flux that drives dianeutral mixing. In particular this iso-neutral density flux 721 is always downwards, and so acts to reduce gravitational potential energy. 722 \subsection{Tapering within the surface mixed layer}\label{sec:triad:taper} 723 724 Additional tapering of the iso-neutral fluxes is necessary within the 725 surface mixed layer. When the Griffies triads are used, we offer two 726 options for this. 727 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper} 728 This is the option activated by the default choice 729 \np{ln\_triad\_iso}=false. Slopes $\tilde{r}_i$ relative to 730 geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the 731 surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values 732 \begin{subequations} 733 \begin{equation} 734 \label{eq:triad:rmtilde} 735 \rMLt = 736 -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for } z>-h, 737 \end{equation} 738 and then the $r_i$ relative to vertical coordinate surfaces are appropriately 739 adjusted to 740 \begin{equation} 741 \label{eq:triad:rm} 742 \rML =\rMLt -\sigma_i \quad \text{ for } z>-h. 743 \end{equation} 744 \end{subequations} 745 Thus the diffusion operator within the mixed layer is given by: 746 \begin{equation} \label{eq:triad:iso_tensor_ML} 747 D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 748 \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 749 1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\ 750 0 \hfill & 1 \hfill & {-\rML[2]} \hfill \\ 751 {-\rML[1]}\hfill & {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill 752 \end{array} }} \right) 753 \end{equation} 754 755 This slope tapering gives a natural connection between tracer in the 756 mixed-layer and in isopycnal layers immediately below, in the 757 thermocline. It is consistent with the way the $\tilde{r}_i$ are 758 tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below) 759 so as to ensure a uniform GM eddy-induced velocity throughout the 760 mixed layer. However, it gives a downwards density flux and so acts so 761 as to reduce potential energy in the same way as does the slope 762 limiting discussed above in \S\ref{sec:triad:limit}. 763 764 As in \S\ref{sec:triad:limit} above, the tapering 765 \eqref{eq:triad:rmtilde} is applied separately to each triad 766 $_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the 767 $_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume 768 $z$-coordinates in the following; the conversion from 769 $\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described 770 above by \eqref{eq:triad:Rtilde}. 771 \begin{enumerate} 772 \item Mixed-layer depth is defined so as to avoid including regions of weak 773 vertical stratification in the slope definition. 774 At each $i,j$ (simplified to $i$ in 775 Fig.~\ref{fig:triad:MLB_triad}), we define the mixed-layer by setting 776 the vertical index of the tracer point immediately below the mixed 777 layer, $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point) 778 such that the potential density 779 ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is 780 the tracer gridbox within which the depth reaches 10~m. See the left 781 side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox 782 instead of the surface gridbox to avoid problems e.g.\ with thin 783 daytime mixed-layers. Currently we use the same 784 $\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is 785 used to output the diagnosed mixed-layer depth 786 $h_\mathrm{ML}=|z_{W}|_{k_\mathrm{ML}+1/2}$, the depth of the $w$-point 787 above the $i,k_\mathrm{ML}$ tracer point. 788 789 \item We define `basal' triad slopes 790 ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ as the slopes 791 of those triads whose vertical `arms' go down from the 792 $i,k_\mathrm{ML}$ tracer point to the $i,k_\mathrm{ML}-1$ tracer point 793 below. This is to ensure that the vertical density gradients 794 associated with these basal triad slopes 795 ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ are 796 representative of the thermocline. The four basal triads defined in the bottom part 797 of Fig.~\ref{fig:triad:MLB_triad} are then 798 \begin{align} 799 {\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p} &= 800 {\:}^{k_\mathrm{ML}-k_p-1/2}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}, \label{eq:triad:Rbase} 801 \\ 802 \intertext{with e.g.\ the green triad} 803 {\:}_i{\mathbb{R}_\mathrm{base}}_{1/2}^{-1/2}&= 804 {\:}^{k_\mathrm{ML}}_i{\mathbb{R}_\mathrm{base}}_{\,1/2}^{-1/2}. \notag 805 \end{align} 806 The vertical flux associated with each of these triads passes through the $w$-point 807 $i,k_\mathrm{ML}-1/2$ lying \emph{below} the $i,k_\mathrm{ML}$ tracer point, 808 so it is this depth 809 \begin{equation} 810 \label{eq:triad:zbase} 811 {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} 812 \end{equation} 813 (one gridbox deeper than the 814 diagnosed ML depth $z_\mathrm{ML})$ that sets the $h$ used to taper 815 the slopes in \eqref{eq:triad:rmtilde}. 816 \item Finally, we calculate the adjusted triads 817 ${\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,i_p}^{k_p}$ within the mixed 818 layer, by multiplying the appropriate 819 ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ by the ratio of 820 the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_\mathrm{base}}_{\,i}$. For 821 instance the green triad centred on $i,k$ 822 \begin{align} 823 {\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,1/2}^{-1/2} &= 824 \frac{{z_w}_{k-1/2}}{{z_\mathrm{base}}_{\,i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{\,1/2}^{-1/2} 825 \notag \\ 826 \intertext{and more generally} 827 {\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,i_p}^{k_p} &= 828 \frac{{z_w}_{k+k_p}}{{z_\mathrm{base}}_{\,i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}.\label{eq:triad:RML} 829 \end{align} 830 \end{enumerate} 831 832 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 833 \begin{figure}[h] 834 \fcapside {\caption{\label{fig:triad:MLB_triad} Definition of 835 mixed-layer depth and calculation of linearly tapered 836 triads. The figure shows a water column at a given $i,j$ 837 (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by 838 bullets, and black lines the edges of the tracer cells; $k$ 839 increases upwards. \\ 840 \hspace{5 em}We define the mixed-layer by setting the vertical index 841 of the tracer point immediately below the mixed layer, 842 $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point) 843 such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, 844 where $i,k_{10}$ is the tracer gridbox within which the depth 845 reaches 10~m. We calculate the triad slopes within the mixed 846 layer by linearly tapering them from zero (at the surface) to 847 the `basal' slopes, the slopes of the four triads passing through the 848 $w$-point $i,k_\mathrm{ML}-1/2$ (blue square), 849 ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$. Triads with 850 different $i_p,k_p$, denoted by different colours, (e.g. the green 851 triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}} 852 {\includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_triad_MLB}} 853 \end{figure} 854 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 855 856 \subsubsection{Additional truncation of skew iso-neutral flux 857 components} 858 \label{sec:triad:Gerdes-taper} 859 The alternative option is activated by setting \np{ln\_triad\_iso} = 860 true. This retains the same tapered slope $\rML$ described above for the 861 calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the 862 vertical tracer flux driven by vertical tracer gradients), but 863 replaces the $\rML$ in the skew term by 864 \begin{equation} 865 \label{eq:triad:rm*} 866 \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i, 867 \end{equation} 868 giving a ML diffusive operator 869 \begin{equation} \label{eq:triad:iso_tensor_ML2} 870 D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 871 \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 872 1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\ 873 0 \hfill & 1 \hfill & {-\rML[2]^*} \hfill \\ 874 {-\rML[1]^*}\hfill & {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\ 875 \end{array} }} \right). 876 \end{equation} 877 This operator 878 \footnote{To ensure good behaviour where horizontal density 879 gradients are weak, we in fact follow \citet{Gerdes1991} and set 880 $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.} 881 then has the property it gives no vertical density flux, and so does 882 not change the potential energy. 883 This approach is similar to multiplying the iso-neutral diffusion 884 coefficient by $\tilde{r}_\mathrm{max}^{-2}\tilde{r}_i^{-2}$ for steep 885 slopes, as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}). 886 Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$ 887 888 In practice, this approach gives weak vertical tracer fluxes through 889 the mixed-layer, as well as vanishing density fluxes. While it is 890 theoretically advantageous that it does not change the potential 891 energy, it may give a discontinuity between the 892 fluxes within the mixed-layer (purely horizontal) and just below (along 893 iso-neutral surfaces). 894 % This may give strange looking results, 895 % particularly where the mixed-layer depth varies strongly laterally. 513 896 % ================================================================ 514 897 % Skew flux formulation for Eddy Induced Velocity : 515 898 % ================================================================ 516 \section{ Eddy induced velocity and Skew flux formulation} 517 518 When Gent and McWilliams's [1990] diffusion is used (\key{traldf\_eiv} defined), 519 an additional advection term is added. The associated velocity is the so called 899 \section{Eddy induced advection formulated as a skew flux}\label{sec:triad:skew-flux} 900 901 \subsection{The continuous skew flux formulation}\label{sec:triad:continuous-skew-flux} 902 903 When Gent and McWilliams's [1990] diffusion is used, 904 an additional advection term is added. The associated velocity is the so called 520 905 eddy induced velocity, the formulation of which depends on the slopes of iso- 521 neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used 522 here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} 906 neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used 907 here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo} 523 908 is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo} 524 + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates. 525 526 The eddy induced velocity is given by: 527 \begin{equation} \label{Eq_eiv_v} 909 + \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates. 910 911 The eddy induced velocity is given by: 912 \begin{subequations} \label{eq:triad:eiv} 913 \begin{equation}\label{eq:triad:eiv_v} 528 914 \begin{split} 529 u^* & = - \frac{1}{e_{ 2}e_{3}}\; \partial_k \left( e_{2} \, A_{e} \; r_i \right)\\530 v^* & = - \frac{1}{e_{ 1}e_{3}}\; \partial_k \left( e_{1} \, A_{e} \; r_j \right)\\531 w^* & = \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i \left( e_{2} \, A_{e} \; r_i \right)532 + \partial_j \left( e_{1} \, A_{e} \;r_j \right) \right\} \\915 u^* & = - \frac{1}{e_{3}}\; \partial_i\psi_1, \\ 916 v^* & = - \frac{1}{e_{3}}\; \partial_j\psi_2, \\ 917 w^* & = \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i \left( e_{2} \, \psi_1\right) 918 + \partial_j \left( e_{1} \, \psi_2\right) \right\}, 533 919 \end{split} 534 920 \end{equation} 535 where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces. 536 537 The traditional way to implement this additional advection is to add it to the Eulerian 538 velocity prior to computing the tracer advection. This allows us to take advantage of 539 all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just 540 a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers 541 where \emph{positivity} of the advection scheme is of paramount importance. 542 543 \citet{Griffies_JPO98} introduces another way to implement the eddy induced advection, 544 the so-called skew form. It is based on a transformation of the advective fluxes 545 using the non-divergent nature of the eddy induced velocity. 546 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be 921 where the streamfunctions $\psi_i$ are given by 922 \begin{equation} \label{eq:triad:eiv_psi} 923 \begin{split} 924 \psi_1 & = A_{e} \; \tilde{r}_1, \\ 925 \psi_2 & = A_{e} \; \tilde{r}_2, 926 \end{split} 927 \end{equation} 928 \end{subequations} 929 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. 930 931 The traditional way to implement this additional advection is to add 932 it to the Eulerian velocity prior to computing the tracer 933 advection. This is implemented if \key{traldf\_eiv} is set in the 934 default implementation, where \np{ln\_traldf\_grif} is set 935 false. This allows us to take advantage of all the advection schemes 936 offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$ 937 order advection scheme. This is particularly useful for passive 938 tracers where \emph{positivity} of the advection scheme is of 939 paramount importance. 940 941 However, when \np{ln\_traldf\_grif} is set true, \NEMO instead 942 implements eddy induced advection according to the so-called skew form 943 \citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes 944 using the non-divergent nature of the eddy induced velocity. 945 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective 946 fluxes per unit area in $ijk$ space can be 547 947 transformed as follows: 548 948 \begin{flalign*} 549 949 \begin{split} 550 \textbf{F}_ {eiv}^T =551 \begin{pmatrix} 950 \textbf{F}_\mathrm{eiv}^T = 951 \begin{pmatrix} 552 952 {e_{2}\,e_{3}\; u^*} \\ 553 953 {e_{1}\,e_{2}\; w^*} \\ 554 954 \end{pmatrix} \; T 555 955 &= 556 \begin{pmatrix} 557 { - \partial_k \left( e_{2} \, A_{e} \; r_i\right) \; T \;} \\558 {+ \partial_i \left( e_{2} \, A_{e} \; r_i\right) \; T \;} \\956 \begin{pmatrix} 957 { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;} \\ 958 {+ \partial_i \left( e_{2} \, \psi_1 \right) \; T \;} \\ 559 959 \end{pmatrix} \\ 560 &= 561 \begin{pmatrix} 562 { - \partial_k \left( e_{2} \, A_{e} \; r_i\; T \right) \;} \\563 {+ \partial_i \left( e_{2} \, A_{e} \; r_i\; T \right) \;} \\564 \end{pmatrix} 565 + 566 \begin{pmatrix} 567 {+ e_{2} \, A_{e} \; r_i\; \partial_k T} \\568 { - e_{2} \, A_{e} \; r_i\; \partial_i T} \\569 \end{pmatrix} 960 &= 961 \begin{pmatrix} 962 { - \partial_k \left( e_{2} \, \psi_1 \; T \right) \;} \\ 963 {+ \partial_i \left( e_{2} \,\psi_1 \; T \right) \;} \\ 964 \end{pmatrix} 965 + 966 \begin{pmatrix} 967 {+ e_{2} \, \psi_1 \; \partial_k T} \\ 968 { - e_{2} \, \psi_1 \; \partial_i T} \\ 969 \end{pmatrix} 570 970 \end{split} 571 971 \end{flalign*} 572 and since the eddy induce s velocity field is no-divergent, we end up with the skew573 form of the eddy induced advective fluxes :574 \begin{equation} \label{ Eq_eiv_skew_continuous}575 \textbf{F}_ {eiv}^T = \begin{pmatrix}576 {+ e_{2} \, A_{e} \; r_i\; \partial_k T} \\577 { - e_{2} \, A_{e} \; r_i\; \partial_i T} \\972 and since the eddy induced velocity field is non-divergent, we end up with the skew 973 form of the eddy induced advective fluxes per unit area in $ijk$ space: 974 \begin{equation} \label{eq:triad:eiv_skew_ijk} 975 \textbf{F}_\mathrm{eiv}^T = \begin{pmatrix} 976 {+ e_{2} \, \psi_1 \; \partial_k T} \\ 977 { - e_{2} \, \psi_1 \; \partial_i T} \\ 578 978 \end{pmatrix} 579 979 \end{equation} 580 The tendency associated with eddy induced velocity is then simply the divergence 581 of the \eqref{Eq_eiv_skew_continuous} fluxes. It naturally conserves the tracer 582 content, as it is expressed in flux form. It also preserve the tracer variance. 583 584 The discrete form of \eqref{Eq_eiv_skew_continuous} using the slopes \eqref{Gf_slopes} and defining $A_e$ at $T$-point is then given by: 585 \begin{flalign} \label{Eq_eiv_skew} \begin{split} 586 \textbf{F}_{eiv}^T \equiv 587 \sum_{\substack{i_p,\,k_p}} \begin{pmatrix} 588 +{e_{2u}}_{i+1/2-i_p}^{k} \ {A_{e}}_{i+1/2-i_p}^{k} 589 \ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} } \ \ \delta_{k+k_p}[T_{i+1/2-i_p}] \\ 590 \\ 591 - {e_{2u}}_i^{k+1/2-k_p} \ {A_{e}}_i^{k+1/2-k_p} 592 \ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} } \ \delta_{i+i_p}[T^{k+1/2-k_p}] \\ 593 \end{pmatrix} 594 \end{split} \end{flalign} 595 Note that \eqref{Eq_eiv_skew} is valid in $z$-coordinate with or without partial cells. In $z^*$ or $s$-coordinate, the the slope between the level and the geopotential surfaces must be added to $\mathbb{R}$ for the discret form to be exact. 596 597 Such a choice of discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes. It also ensures the conservation of the tracer variance (see Appendix \ref{Apdx_eiv_skew}), $i.e.$ it does not include a diffusive component but is a `pure' advection term. 598 599 600 601 602 \newpage %force an empty line 603 % ================================================================ 604 % Discrete Invariants of the skew flux formulation 605 % ================================================================ 606 \subsection{Discrete Invariants of the skew flux formulation} 607 \label{Apdx_eiv_skew} 608 609 610 Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. 611 612 This have to be moved in an Appendix. 613 614 The continuous property to be demonstrated is : 615 \begin{align*} 616 \int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv \equiv 0 617 \end{align*} 618 The discrete form of its left hand side is obtained using \eqref{Eq_eiv_skew} 619 \begin{align*} 620 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; 621 \delta_i &\left[ 622 {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} 623 \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] 624 \right] \; T_i^k \\ 625 - \delta_k &\left[ 626 {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} 627 \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] 628 \right] \; T_i^k \ \Biggr\} 629 \end{align*} 630 apply the adjoint of delta operator, it becomes 631 \begin{align*} 632 \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}} \Biggl\{ \;\; 633 &\left( 634 {e_{2u}}_{i+i_p+1/2}^{k} \;\ \ {A_{e}}_{i+i_p+1/2}^{k} 635 \ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} } \quad \delta_{k+k_p}[T_{i+i_p+1/2}] 636 \right) \; \delta_{i+1/2}[T^{k}] \\ 637 - &\left( 638 {e_{2u}}_i^{k+k_p+1/2} \ \ {A_{e}}_i^{k+k_p+1/2} 639 \ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} } \ \ \delta_{i+i_p}[T^{k+k_p+1/2}] 640 \right) \; \delta_{k+1/2}[T_{i}] \ \Biggr\} 641 \end{align*} 642 Expending the summation on $i_p$ and $k_p$, it becomes: 643 \begin{align*} 644 \begin{matrix} 645 &\sum\limits_{i,k} \Bigl\{ 646 &+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} 647 &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ 648 &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} 649 &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}} &\delta_{k-1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ 650 &&+{e_{2u}}_{i+1}^{k} &{A_{e}}_{i+1 }^{k} 651 &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}] &\delta_{i+1/2}[T^{k}] &\\ 652 &&+{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_{i}^{k\ \ \ \:} 653 &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ 654 % 655 &&-{e_{2u}}_i^{k+1} &{A_{e}}_i^{k+1} 656 &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}} &\delta_{i-1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ 657 &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} 658 &{\ \ \;_i^k \mathbb{R}_{-1/2}^{+1/2}} &\delta_{i-1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] &\\ 659 &&-{e_{2u}}_i^{k+1 } &{A_{e}}_i^{k+1} 660 &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}} &\delta_{i+1/2}[T^{k+1}] &\delta_{k+1/2}[T_{i}] &\\ 661 &&-{e_{2u}}_i^{k\ \ \ \:} &{A_{e}}_i^{k\ \ \ \:} 662 &{\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{i+1/2}[T^{k\ \ \ \:}] &\delta_{k+1/2}[T_{i}] 663 &\Bigr\} \\ 664 \end{matrix} 665 \end{align*} 666 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the 667 same but of opposite signs, they cancel out. 668 Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. 669 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the 670 same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$ 671 they cancel out with the neighbouring grid points. 672 Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the 673 $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the 674 tracer is preserved by the discretisation of the skew fluxes. 675 676 %%% Local Variables: 677 %%% TeX-master: "../../NEMO_book.tex" 678 %%% End: 980 The total fluxes per unit physical area are then 981 \begin{equation}\label{eq:triad:eiv_skew_physical} 982 \begin{split} 983 f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T \\ 984 f^*_2 & = \frac{1}{e_{3}}\; \psi_2 \partial_k T \\ 985 f^*_3 & = -\frac{1}{e_{1}e_{2}}\; \left\{ e_{2} \psi_1 \partial_i T 986 + e_{1} \psi_2 \partial_j T \right\}. \\ 987 \end{split} 988 \end{equation} 989 Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the 990 vertical coordinate, though of course the slopes 991 $\tilde{r}_i$ which define the $\psi_i$ in \eqref{eq:triad:eiv_psi} are relative to geopotentials. 992 The tendency associated with eddy induced velocity is then simply the convergence 993 of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so 994 \begin{equation} \label{eq:triad:skew_eiv_conv} 995 \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 } \left[ 996 \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right) 997 + \frac{\partial}{\partial j} \left( e_1 \; 998 \psi_2 \partial_k T\right) 999 - \frac{\partial}{\partial k} \left( e_{2} \psi_1 \partial_i T 1000 + e_{1} \psi_2 \partial_j T \right) \right] 1001 \end{equation} 1002 It naturally conserves the tracer content, as it is expressed in flux 1003 form. Since it has the same divergence as the advective form it also 1004 preserves the tracer variance. 1005 1006 \subsection{The discrete skew flux formulation} 1007 The skew fluxes in (\ref{eq:triad:eiv_skew_physical}, \ref{eq:triad:eiv_skew_ijk}), like the off-diagonal terms 1008 (\ref{eq:triad:i13c}, \ref{eq:triad:i31c}) of the small angle diffusion tensor, are best 1009 expressed in terms of the triad slopes, as in Fig.~\ref{fig:triad:ISO_triad} 1010 and Eqs~(\ref{eq:triad:i13}, \ref{eq:triad:i31}); but now in terms of the triad slopes 1011 $\tilde{\mathbb{R}}$ relative to geopotentials instead of the 1012 $\mathbb{R}$ relative to coordinate surfaces. The discrete form of 1013 \eqref{eq:triad:eiv_skew_ijk} using the slopes \eqref{eq:triad:R} and 1014 defining $A_e$ at $T$-points is then given by: 1015 1016 1017 \begin{subequations}\label{eq:triad:allskewflux} 1018 \begin{flalign}\label{eq:triad:vect_skew_flux} 1019 \vect{F}_\mathrm{eiv}(T) &\equiv 1020 \sum_{\substack{i_p,\,k_p}} 1021 \begin{pmatrix} 1022 {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T) \\ 1023 \\ 1024 {_i^{k+1/2-k_p} {\mathbb{S}_w}_{i_p}^{k_p} } (T) \\ 1025 \end{pmatrix}, 1026 \end{flalign} 1027 where the skew flux in the $i$-direction associated with a given 1028 triad is (\ref{eq:triad:latflux-triad}, \ref{eq:triad:triadfluxu}): 1029 \begin{align} 1030 \label{eq:triad:skewfluxu} 1031 _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{ 1032 \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 1033 \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \ 1034 \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, 1035 \\ 1036 \intertext{and \eqref{eq:triad:triadfluxw} in the $k$-direction, changing the sign 1037 to be consistent with \eqref{eq:triad:eiv_skew_ijk}:} 1038 _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T) 1039 &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} 1040 {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:triad:skewfluxw} 1041 \end{align} 1042 \end{subequations} 1043 1044 Such a discretisation is consistent with the iso-neutral 1045 operator as it uses the same definition for the slopes. It also 1046 ensures the following two key properties. 1047 \subsubsection{No change in tracer variance} 1048 The discretization conserves tracer variance, $i.e.$ it does not 1049 include a diffusive component but is a `pure' advection term. This can 1050 be seen 1051 %either from Appendix \ref{Apdx_eiv_skew} or 1052 by considering the 1053 fluxes associated with a given triad slope 1054 $_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following 1055 \S\ref{sec:triad:variance} and \eqref{eq:triad:dvar_iso_i}, the 1056 associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$ 1057 drives a net rate of change of variance, summed over the two 1058 $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 1059 \begin{equation} 1060 \label{eq:triad:dvar_eiv_i} 1061 _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], 1062 \end{equation} 1063 while the associated vertical skew-flux gives a variance change summed over the 1064 $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 1065 \begin{equation} 1066 \label{eq:triad:dvar_eiv_k} 1067 _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 1068 \end{equation} 1069 Inspection of the definitions (\ref{eq:triad:skewfluxu}, \ref{eq:triad:skewfluxw}) 1070 shows that these two variance changes (\ref{eq:triad:dvar_eiv_i}, \ref{eq:triad:dvar_eiv_k}) 1071 sum to zero. Hence the two fluxes associated with each triad make no 1072 net contribution to the variance budget. 1073 1074 \subsubsection{Reduction in gravitational PE} 1075 The vertical density flux associated with the vertical skew-flux 1076 always has the same sign as the vertical density gradient; thus, so 1077 long as the fluid is stable (the vertical density gradient is 1078 negative) the vertical density flux is negative (downward) and hence 1079 reduces the gravitational PE. 1080 1081 For the change in gravitational PE driven by the $k$-flux is 1082 \begin{align} 1083 \label{eq:triad:vert_densityPE} 1084 g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) 1085 &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k 1086 {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k 1087 {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\ 1088 \intertext{Substituting ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from 1089 \eqref{eq:triad:skewfluxw}, gives} 1090 % and separating out 1091 % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, 1092 % gives two terms. The 1093 % first $\rtriad{R}$ term (the only term for $z$-coordinates) is: 1094 &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} 1095 \frac{ -\alpha _i^k\delta_{i+ i_p}[T^k]+ \beta_i^k\delta_{i+ i_p}[S^k]} { {e_{1u}}_{\,i + i_p}^{\,k} } \notag \\ 1096 &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 1097 \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) {_i^k\mathbb{R}_{i_p}^{k_p}} 1098 \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}}, 1099 \end{align} 1100 using the definition of the triad slope $\rtriad{R}$, 1101 \eqref{eq:triad:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+ 1102 \beta_i^k\delta_{i+ i_p}[S^k]$ in terms of $-\alpha_i^k \delta_{k+ 1103 k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$. 1104 1105 Where the coordinates slope, the $i$-flux gives a PE change 1106 \begin{multline} 1107 \label{eq:triad:lat_densityPE} 1108 g \delta_{i+i_p}[z_T^k] 1109 \left[ 1110 -\alpha _i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (S) 1111 \right] \\ 1112 = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 1113 \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} 1114 \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) 1115 \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}}, 1116 \end{multline} 1117 (using \eqref{eq:triad:skewfluxu}) and so the total PE change 1118 \eqref{eq:triad:vert_densityPE} + \eqref{eq:triad:lat_densityPE} associated with the triad fluxes is 1119 \begin{multline} 1120 \label{eq:triad:tot_densityPE} 1121 g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) + 1122 g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ 1123 = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 1124 \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)^2 1125 \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}}. 1126 \end{multline} 1127 Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+ 1128 \beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative. 1129 1130 \subsection{Treatment of the triads at the boundaries}\label{sec:triad:skew_bdry} 1131 Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes 1132 are masked at the boundaries in exactly the same way as are the triad 1133 slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as 1134 described in \S\ref{sec:triad:iso_bdry} and 1135 Fig.~\ref{fig:triad:bdry_triads}. Thus surface layer triads 1136 $\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are 1137 masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$ 1138 and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the 1139 $i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ 1140 $u$-point is masked. The namelist parameter \np{ln\_botmix\_grif} has 1141 no effect on the eddy-induced skew-fluxes. 1142 1143 \subsection{ Limiting of the slopes within the interior}\label{sec:triad:limitskew} 1144 Presently, the iso-neutral slopes $\tilde{r}_i$ relative 1145 to geopotentials are limited to be less than $1/100$, exactly as in 1146 calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each 1147 individual triad \rtriadt{R} is so limited. 1148 1149 \subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew} 1150 The slopes $\tilde{r}_i$ relative to 1151 geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the 1152 surface \eqref{eq:triad:rmtilde}, as described in \S\ref{sec:triad:lintaper}. This is 1153 option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the 1154 slopes used to calculate the eddy-induced fluxes is 1155 unaffected by the value of \np{ln\_triad\_iso}. 1156 1157 The justification for this linear slope tapering is that, for $A_e$ 1158 that is constant or varies only in the horizontal (the most commonly 1159 used options in \NEMO: see \S\ref{LDF_coef}), it is 1160 equivalent to a horizontal eiv (eddy-induced velocity) that is uniform 1161 within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the 1162 eiv velocities do not restratify the mixed layer \citep{Treguier1997, 1163 Danabasoglu_al_2008}. Equivantly, in terms 1164 of the skew-flux formulation we use here, the 1165 linear slope tapering within the mixed-layer gives a linearly varying 1166 vertical flux, and so a tracer convergence uniform in depth (the 1167 horizontal flux convergence is relatively insignificant within the mixed-layer). 1168 1169 \subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag} 1170 Where the namelist parameter \np{ln\_traldf\_gdia}=true, diagnosed 1171 mean eddy-induced velocities are output. Each time step, 1172 streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at 1173 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ 1174 (integer $i$, integer +1/2 $j$, integer +1/2 $k$) points (see Table 1175 \ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and 1176 calculate the streamfunction at a given $uw$-point from the 1177 surrounding four triads according to: 1178 \begin{equation} 1179 \label{eq:triad:sfdiagi} 1180 {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}} 1181 {A_e}_{i+1/2-i_p}^{k+1/2-k_p}\:\triadd{i+1/2-i_p}{k+1/2-k_p}{R}{i_p}{k_p}. 1182 \end{equation} 1183 The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 1184 The eddy-induced velocities are then calculated from the 1185 straightforward discretisation of \eqref{eq:triad:eiv_v}: 1186 \begin{equation}\label{eq:triad:eiv_v_discrete} 1187 \begin{split} 1188 {u^*}_{i+1/2}^{k} & = - \frac{1}{{e_{3u}}_{i}^{k}}\left({\psi_1}_{i+1/2}^{k+1/2}-{\psi_1}_{i+1/2}^{k+1/2}\right), \\ 1189 {v^*}_{j+1/2}^{k} & = - \frac{1}{{e_{3v}}_{j}^{k}}\left({\psi_2}_{j+1/2}^{k+1/2}-{\psi_2}_{j+1/2}^{k+1/2}\right), \\ 1190 {w^*}_{i,j}^{k+1/2} & = \frac{1}{e_{1t}e_{2t}}\; \left\{ 1191 {e_{2u}}_{i+1/2}^{k+1/2} \,{\psi_1}_{i+1/2}^{k+1/2} - 1192 {e_{2u}}_{i-1/2}^{k+1/2} \,{\psi_1}_{i-1/2}^{k+1/2} \right. + \\ 1193 \phantom{=} & \qquad\qquad\left. {e_{2v}}_{j+1/2}^{k+1/2} \,{\psi_2}_{j+1/2}^{k+1/2} - {e_{2v}}_{j-1/2}^{k+1/2} \,{\psi_2}_{j-1/2}^{k+1/2} \right\}, 1194 \end{split} 1195 \end{equation}
Note: See TracChangeset
for help on using the changeset viewer.
