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