Ignore:
Timestamp:
2012-01-18T11:28:11+01:00 (9 years ago)
Author:
agn
Message:

dev_NEMO_MERGE_2011:(DOC) Finish griffies description, remove need for pstricks

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_ISO.tex

    r3218 r3267  
    22% Iso-neutral diffusion : 
    33% ================================================================ 
    4 \chapter{Griffies's iso-neutral diffusion} 
     4\chapter[Iso-neutral diffusion and eddy advection using 
     5triads]{Iso-neutral diffusion and eddy advection using triads} 
    56\label{sec:triad} 
    67\minitoc 
    7  
    8 \section{Griffies's formulation of iso-neutral diffusion} 
     8\pagebreak 
     9\section{Choice of namelist parameters} 
     10%-----------------------------------------nam_traldf------------------------------------------------------ 
     11\namdisplay{namtra_ldf} 
     12%--------------------------------------------------------------------------------------------------------- 
     13If the namelist variable \np{ln\_traldf\_grif} is set true (and 
     14\key{ldfslp} is set), \NEMO updates both active and passive tracers 
     15using the Griffies triad representation of iso-neutral diffusion and 
     16the eddy-induced advective skew (GM) fluxes. Otherwise (by default) the 
     17filtered version of Cox's original scheme is employed 
     18(\S\ref{LDF_slp}). In the present implementation of the Griffies 
     19scheme, the advective skew fluxes are implemented even if 
     20\key{traldf\_eiv} is not set. 
     21 
     22Values of iso-neutral diffusivity and GM coefficient are set as 
     23described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd}, 
     24N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and 
     25GM 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 
     28scale 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} 
     33is 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 
     35instead 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 
     39unless it is zero. 
     40 
     41The 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} 
    965\label{sec:triad:iso} 
    10  
    11 We define a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO 
     66We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO 
    1267framework, using scale factors rather than grid-sizes. 
    1368 
     
    613668or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is 
    614669masked. 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}), 
     670masked if \np{ln\_botmix\_grif}=false, but left unmasked, 
     671giving bottom mixing, if \np{ln\_botmix\_grif}=true. 
     672 
     673The default option \np{ln\_botmix\_grif}=false is suitable when the 
     674bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1), 
    620675or  for simple idealized  problems. For setups with topography without 
    621 bbl mixing, \nlv{ln\_botmix\_grif=.true.} may be necessary. 
     676bbl mixing, \np{ln\_botmix\_grif}=true may be necessary. 
    622677% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    623678\begin{figure}[h] \begin{center} 
     
    636691      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point 
    637692      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.}} 
     693      line) are masked if \np{botmix\_grif}=.false., but left 
     694      unmasked, giving bottom mixing, if \np{botmix\_grif}=.true.} 
    640695 \end{center} \end{figure} 
    641696% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    665720iso-neutral density flux that drives dianeutral mixing.  In particular this iso-neutral density flux 
    666721is always downwards, and so acts to reduce gravitational potential energy. 
    667 \subsection{Tapering within the surface mixed layer} 
     722\subsection{Tapering within the surface mixed layer}\label{sec:triad:taper} 
     723 
    668724Additional tapering of the iso-neutral fluxes is necessary within the 
    669725surface mixed layer. When the Griffies triads are used, we offer two 
     
    671727\subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper} 
    672728This is the option activated by the default choice 
    673 \nlv{ln\_triad\_iso=.false.}. Slopes $\tilde{r}_i$ relative to 
     729\np{ln\_triad\_iso}=false. Slopes $\tilde{r}_i$ relative to 
    674730geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the 
    675731surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values 
     
    798854% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    799855 
    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 
     856\subsubsection{Additional truncation of skew iso-neutral flux 
     857  components} 
     858\label{sec:triad:Gerdes-taper} 
     859The alternative option is activated by setting \np{ln\_triad\_iso} = 
     860  true. This retains the same tapered slope $\rML$  described above for the 
    803861calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the 
    804862vertical tracer flux driven by vertical tracer gradients), but 
     
    839897% Skew flux formulation for Eddy Induced Velocity : 
    840898% ================================================================ 
    841 \section{Eddy induced advection and its formulation as a skew flux} 
     899\section{Eddy induced advection formulated as a skew flux}\label{sec:triad:skew-flux} 
    842900 
    843901\subsection{The continuous skew flux formulation}\label{sec:triad:continuous-skew-flux} 
    844902 
    845  When Gent and McWilliams's [1990] diffusion is used (\key{traldf\_eiv} defined), 
     903 When Gent and McWilliams's [1990] diffusion is used, 
    846904an additional advection term is added. The associated velocity is the so called 
    847905eddy induced velocity, the formulation of which depends on the slopes of iso- 
     
    852910 
    853911The eddy induced velocity is given by: 
    854 \begin{equation} \label{eq:triad:eiv_v} 
     912\begin{subequations} \label{eq:triad:eiv} 
     913\begin{equation}\label{eq:triad:eiv_v} 
    855914\begin{split} 
    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\} 
     915 u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\ 
     916 v^* & = - \frac{1}{e_{3}}\;          \partial_j\psi_2,    \\ 
     917w^* & =    \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\}, 
    860919\end{split} 
    861920\end{equation} 
    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 
     921where 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} 
     929with $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 
     931The traditional way to implement this additional advection is to add 
     932it to the Eulerian velocity prior to computing the tracer 
     933advection. This is implemented if \key{traldf\_eiv} is set in the 
     934default implementation, where \np{ln\_traldf\_grif} is set 
     935false. This allows us to take advantage of all the advection schemes 
     936offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$ 
     937order advection scheme. This is particularly useful for passive 
     938tracers where \emph{positivity} of the advection scheme is of 
     939paramount importance. 
     940 
     941However, when \np{ln\_traldf\_grif} is set true, \NEMO instead 
     942implements eddy induced advection according to the so-called skew form 
     943\citep{Griffies_JPO98}. It is based on a transformation of the advective fluxes 
    872944using the non-divergent nature of the eddy induced velocity. 
    873945For example in the (\textbf{i},\textbf{k}) plane, the tracer advective 
     
    883955&= 
    884956\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 \;}  \\ 
     957           { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;}     \\ 
     958      {+ \partial_i  \left( e_{2} \, \psi_1 \right) \; T \;}    \\ 
    887959\end{pmatrix}        \\ 
    888960&= 
    889961\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) \;}    \\ 
     962           { - \partial_k \left( e_{2} \, \psi_1  \; T \right) \;}  \\ 
     963      {+ \partial_i  \left( e_{2} \,\psi_1 \; T \right) \;}  \\ 
    892964\end{pmatrix} 
    893965 + 
    894966\begin{pmatrix} 
    895            {+ e_{2} \, A_{e} \; \tilde{r}_1  \; \partial_k T}  \\ 
    896       { - e_{2} \, A_{e} \; \tilde{r}_1  \; \partial_i  T}   \\ 
     967           {+ e_{2} \, \psi_1  \; \partial_k T}  \\ 
     968      { - e_{2} \, \psi_1  \; \partial_i  T}  \\ 
    897969\end{pmatrix} 
    898970\end{split} 
     
    902974\begin{equation} \label{eq:triad:eiv_skew_ijk} 
    903975\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}   \\ 
     976           {+ e_{2} \, \psi_1  \; \partial_k T}   \\ 
     977      { - e_{2} \, \psi_1  \; \partial_i  T}  \\ 
    906978                                 \end{pmatrix} 
    907979\end{equation} 
     
    909981\begin{equation}\label{eq:triad:eiv_skew_physical} 
    910982\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\}. \\ 
     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\}. \\ 
    915987\end{split} 
    916988\end{equation} 
    917989Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the 
    918990vertical coordinate, though of course the slopes 
    919 $\tilde{r}_i$ are relative to geopotentials. 
     991$\tilde{r}_i$ which define the $\psi_i$ in \eqref{eq:triad:eiv_psi} are relative to geopotentials. 
    920992The tendency associated with eddy induced velocity is then simply the convergence 
    921993of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so 
    922994\begin{equation} \label{eq:triad:skew_eiv_conv} 
    923995\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] 
     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] 
    9291001\end{equation} 
    9301002 It naturally conserves the tracer content, as it is expressed in flux 
     
    9761048The discretization conserves tracer variance, $i.e.$ it does not 
    9771049include 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 
     1050be seen 
     1051%either from Appendix \ref{Apdx_eiv_skew} or 
     1052by considering the 
    9791053fluxes associated with a given triad slope 
    9801054$_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following 
     
    10641138and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the 
    10651139$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 
     1140$u$-point is masked. The namelist parameter \np{ln\_botmix\_grif} has 
    10671141no effect on the eddy-induced skew-fluxes. 
    10681142 
     
    10791153option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the 
    10801154slopes used to calculate the eddy-induced fluxes is 
    1081 unaffected by the value of \nlv{ln\_triad\_iso}. 
     1155unaffected by the value of \np{ln\_triad\_iso}. 
    10821156 
    10831157The justification for this linear slope tapering is that, for $A_e$ 
     
    10941168 
    10951169\subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag} 
    1096 Where the namelist parameter \nlv{ln\_botmix\_grif=.true.}, diagnosed 
     1170Where the namelist parameter \np{ln\_botmix\_grif}=true, diagnosed 
    10971171mean eddy-induced velocities are output. Each time step, 
    10981172streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at 
     
    11041178\begin{equation} 
    11051179  \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} 
    1109  
    1110 \newpage      %force an empty line 
    1111 % ================================================================ 
    1112 % Discrete Invariants of the skew flux formulation 
    1113 % ================================================================ 
    1114 \subsection{Discrete Invariants of the skew flux formulation} 
    1115 \label{Apdx_eiv_skew} 
    1116  
    1117  
    1118 Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. 
    1119  
    1120 This have to be moved in an Appendix. 
    1121  
    1122 The continuous property to be demonstrated is : 
    1123 \begin{align*} 
    1124 \int_D \nabla \cdot \textbf{F}_\mathrm{eiv}(T) \; T \;dv  \equiv 0 
    1125 \end{align*} 
    1126 The discrete form of its left hand side is obtained using \eqref{eq:triad:allskewflux} 
    1127 \begin{align*} 
    1128  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\; 
    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}] 
    1132    \right] \; T_i^k      \\ 
    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\} 
    1137 \end{align*} 
    1138 apply the adjoint of delta operator, it becomes 
    1139 \begin{align*} 
    1140  \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\; 
    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}] 
    1144    \right) \; \delta_{i+1/2}[T^{k}]      \\ 
    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\} 
    1149 \end{align*} 
    1150 Expending the summation on $i_p$ and $k_p$, it becomes: 
    1151 \begin{align*} 
    1152  \begin{matrix} 
    1153 &\sum\limits_{i,k}   \Bigl\{ 
    1154   &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k} 
    1155   &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\ 
    1156 &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:} 
    1157   &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}]  &\delta_{i+1/2}[T^{k}] &\\ 
    1158 &&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k} 
    1159   &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\ 
    1160 &&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:} 
    1161     &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\ 
    1162 % 
    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\ \ \ \:} 
    1166   &{\ \ \;_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}] 
    1171 &\Bigr\}  \\ 
    1172 \end{matrix} 
    1173 \end{align*} 
    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 
    1182 tracer is preserved by the discretisation of the skew fluxes. 
    1183  
    1184 %%% Local Variables: 
    1185 %%% TeX-master: "../../NEMO_book-luatex.tex" 
    1186 %%% End: 
     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} 
     1183The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 
     1184The eddy-induced velocities are then calculated from the 
     1185straightforward discretisation of \eqref{eq:triad:eiv_v}: 
     1186\begin{equation}\label{eq:triad:eiv_v} 
     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.