New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10414 for NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex – NEMO

Ignore:
Timestamp:
2018-12-19T00:02:00+01:00 (5 years ago)
Author:
nicolasmartin
Message:
  • Comment \label commands on maths environments for unreferenced equations and adapt the unnumbered math container accordingly (mainly switch to shortanded LateX syntax with \[ ... \])
  • Add a code trick to build subfile with its own bibliography
  • Fix right path for main LaTeX document in first line of subfiles (\documentclass[...]{subfiles})
  • Rename abstract_foreword.tex to foreword.tex
  • Fix some non-ASCII codes inserted here or there in LaTeX (\[0-9]*)
  • Made a first iteration on the indentation and alignement within math, figure and table environments to improve source code readability
File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/doc/latex/NEMO/subfiles/annex_iso.tex

    r10406 r10414  
    1 \documentclass[../tex_main/NEMO_manual]{subfiles} 
     1\documentclass[../main/NEMO_manual]{subfiles} 
     2 
    23\begin{document} 
    34% ================================================================ 
    … …  
    78         {\texorpdfstring{Iso-Neutral Diffusion and\\ Eddy Advection using Triads}{Iso-Neutral Diffusion and Eddy Advection using Triads}} 
    89\label{apdx:triad} 
     10 
    911\minitoc 
    10 \pagebreak 
     12 
     13\newpage 
     14 
    1115\section{Choice of \protect\ngn{namtra\_ldf} namelist parameters} 
    1216%-----------------------------------------nam_traldf------------------------------------------------------ 
    … …  
    5963\section{Triad formulation of iso-neutral diffusion} 
    6064\label{sec:iso} 
     65 
    6166We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, 
    6267but formulated within the \NEMO framework, using scale factors rather than grid-sizes. 
    6368 
    6469\subsection{Iso-neutral diffusion operator} 
     70 
    6571The iso-neutral second order tracer diffusive operator for small angles between 
    6672iso-neutral surfaces and geopotentials is given by \autoref{eq:iso_tensor_1}: 
    67 \begin{subequations} \label{eq:iso_tensor_1} 
     73\begin{subequations} 
     74  \label{eq:iso_tensor_1} 
    6875  \begin{equation} 
    6976    D^{lT}=-\Div\vect{f}^{lT}\equiv 
    … …  
    7986    \mbox{with}\quad \;\;\Re = 
    8087    \begin{pmatrix} 
    81        1   &  0   & -r_1           \mystrut \\ 
    82        0   &  1   & -r_2           \mystrut \\ 
     88      1   &  0   & -r_1           \mystrut \\ 
     89      0   &  1   & -r_2           \mystrut \\ 
    8390      -r_1 & -r_2 &  r_1 ^2+r_2 ^2 \mystrut 
    8491    \end{pmatrix} 
    … …  
    8895      \frac{1}{e_2} \pd[T]{j} \mystrut \\ 
    8996      \frac{1}{e_3} \pd[T]{k} \mystrut 
    90     \end{pmatrix}. 
     97    \end{pmatrix} 
     98    . 
    9199  \end{equation} 
    92100\end{subequations} 
    … …  
    99107\begin{align*} 
    100108  r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i} 
    101   \right) 
    102   \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \\ 
    103   &=-\frac{e_3 }{e_1 } \left( -\alpha\frac{\partial T }{\partial i} + 
    104     \beta\frac{\partial S }{\partial i} \right) \left( 
    105     -\alpha\frac{\partial T }{\partial k} + \beta\frac{\partial S 
    106     }{\partial k} \right)^{-1} 
     109        \right) 
     110        \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \\ 
     111      &=-\frac{e_3 }{e_1 } \left( -\alpha\frac{\partial T }{\partial i} + 
     112        \beta\frac{\partial S }{\partial i} \right) \left( 
     113        -\alpha\frac{\partial T }{\partial k} + \beta\frac{\partial S 
     114        }{\partial k} \right)^{-1} 
    107115\end{align*} 
    108116is the $i$-component of the slope of the iso-neutral surface relative to the computational surface, 
    … …  
    110118 
    111119We will find it useful to consider the fluxes per unit area in $i,j,k$ space; we write 
    112 \begin{equation} 
    113   \label{eq:Fijk} 
     120\[ 
     121  % \label{eq:Fijk} 
    114122  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right). 
    115 \end{equation} 
     123\] 
    116124Additionally, we will sometimes write the contributions towards the fluxes $\vect{f}$ and 
    117125$\vect{F}_{\mathrm{iso}}$ from the component $R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, 
    … …  
    124132  \label{eq:i13c} 
    125133  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}\\ 
    126 \intertext{and in the k-direction resulting from the lateral tracer gradients} 
     134  \intertext{and in the k-direction resulting from the lateral tracer gradients} 
    127135  \label{eq:i31c} 
    128  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} 
     136  f_{31}+f_{32}=& \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i} 
    129137\end{align} 
    130138 
    … …  
    147155 
    148156\subsection{Standard discretization} 
     157 
    149158The straightforward approach to discretize the lateral skew flux 
    150159\autoref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995 into OPA, 
    … …  
    163172\[ 
    164173  \overline{\overline 
    165    r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k} 
     174    r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k} 
    166175  \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}}, 
    167176\] 
    … …  
    177186 
    178187\subsection{Expression of the skew-flux in terms of triad slopes} 
     188 
    179189\citep{Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that 
    180190nicely solves the problem. 
    … …  
    182192% the mean vertical gradient at the $u$-point, 
    183193% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    184 \begin{figure}[tb] \begin{center} 
     194\begin{figure}[tb] 
     195  \begin{center} 
    185196    \includegraphics[width=1.05\textwidth]{Fig_GRIFF_triad_fluxes} 
    186     \caption{ \protect\label{fig:ISO_triad} 
     197    \caption{ 
     198      \protect\label{fig:ISO_triad} 
    187199      (a) Arrangement of triads $S_i$ and tracer gradients to 
    188            give lateral tracer flux from box $i,k$ to $i+1,k$ 
     200      give lateral tracer flux from box $i,k$ to $i+1,k$ 
    189201      (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from 
    190             box $i,k$ to $i,k+1$.} 
    191  \end{center} \end{figure} 
     202      box $i,k$ to $i,k+1$. 
     203    } 
     204  \end{center} 
     205\end{figure} 
    192206% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    193207They get the skew flux from the products of the vertical gradients at each $w$-point surrounding the $u$-point with 
    … …  
    204218  _{k+\frac{1}{2}} \left[ T^i 
    205219  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\ 
    206    +\Alts _{i+1}^k a_3 s_3 \delta_{k-\frac{1}{2}} \left[ T^{i+1} 
     220  +\Alts _{i+1}^k a_3 s_3 \delta_{k-\frac{1}{2}} \left[ T^{i+1} 
    207221  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  +\Alts _i^k a_4 s_4 \delta 
    208222  _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}, 
    … …  
    219233  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \Alts_i^{k+1} a_{1}' 
    220234  s_{1}' \delta_{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1} 
    221    +\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}\\ 
     235  +\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} \\ 
    222236  + \Alts_i^k a_{3}' s_{3}' \delta_{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k 
    223237  +\Alts_i^k a_{4}' s_{4}' \delta_{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k. 
    … …  
    242256 
    243257% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    244 \begin{figure}[tb] \begin{center} 
     258\begin{figure}[tb] 
     259  \begin{center} 
    245260    \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells} 
    246     \caption{   \protect\label{fig:qcells} 
     261    \caption{ 
     262      \protect\label{fig:qcells} 
    247263      Triad notation for quarter cells. $T$-cells are inside boxes, 
    248264      while the  $i+\half,k$ $u$-cell is shaded in green and 
    249       the $i,k+\half$ $w$-cell is shaded in pink.} 
    250   \end{center} \end{figure} 
     265      the $i,k+\half$ $w$-cell is shaded in pink. 
     266    } 
     267  \end{center} 
     268\end{figure} 
    251269% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    252270 
    … …  
    266284 
    267285\subsection{Full triad fluxes} 
     286 
    268287A key property of iso-neutral diffusion is that it should not affect the (locally referenced) density. 
    269288In particular there should be no lateral or vertical density flux. 
    … …  
    306325  \label{eq:i33} 
    307326  \left( F_w^{33} \right) _i^{k+\frac{1}{2}} = 
    308     - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2 
     327  - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2 
    309328    + \Alts_i^{k+1} a_{2}' s_{2}'^2 
    310329    + \Alts_i^k a_{3}' s_{3}'^2 
    … …  
    318337  _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 
    319338  &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p} 
    320   \left( 
     339    \left( 
    321340    {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    322341    -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 
    323342    \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    324   \right) \\ 
     343    \right) \\ 
    325344  &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right) 
    326    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} 
     345    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2} 
    327346\end{align} 
    328347may be associated with each triad. 
    … …  
    338357the iso-neutral fluxes at $u$- and $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces: 
    339358%(\autoref{fig:ISO_triad}): 
    340 \begin{flalign} \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
     359\begin{flalign} 
     360  \label{eq:iso_flux} \vect{F}_{\mathrm{iso}}(T) &\equiv 
    341361  \sum_{\substack{i_p,\,k_p}} 
    342362  \begin{pmatrix} 
    343     {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\ 
    344     \\ 
    345     {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)      \\ 
     363    {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ \\ 
     364    {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) \\ 
    346365  \end{pmatrix}. 
    347366\end{flalign} 
    … …  
    369388  \quad {b_T}_{i+i_p+1/2}^k\left(\frac{\partial T}{\partial 
    370389      t}T\right)_{i+i_p+1/2}^k \\ 
    371  \begin{aligned} 
    372   &= -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 
    373   {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ 
    374   &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:dvar_iso_i} 
    375  \end{aligned} 
     390  \begin{aligned} 
     391    &= -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 
     392    {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\ 
     393    &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:dvar_iso_i} 
     394  \end{aligned} 
    376395\end{multline} 
    377396while the vertical flux similarly drives a net rate of change of variance summed over 
    378397the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 
    379398\begin{equation} 
    380 \label{eq:dvar_iso_k} 
     399  \label{eq:dvar_iso_k} 
    381400  _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    382401\end{equation} 
    … …  
    385404\autoref{eq:latflux-triad} and \autoref{eq:vertflux-triad}, it is 
    386405\begin{multline*} 
    387 -\Alts_i^k\left \{ 
    388 { } _i^k{\mathbb{A}_u}_{i_p}^{k_p} 
    389   \left( 
    390     \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    391     - {_i^k\mathbb{R}_{i_p}^{k_p}} \ 
    392     \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }\right)\,\delta_{i+ i_p}[T^k] \right.\\ 
    393 - \left. { } _i^k{\mathbb{A}_w}_{i_p}^{k_p} 
    394   \left( 
    395     \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    396     -{\:}_i^k\mathbb{R}_{i_p}^{k_p} 
    397     \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    398   \right) {\,}_i^k\mathbb{R}_{i_p}^{k_p}\delta_{k+ k_p}[T^i] 
    399 \right \}. 
     406  -\Alts_i^k\left \{ 
     407    { } _i^k{\mathbb{A}_u}_{i_p}^{k_p} 
     408    \left( 
     409      \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
     410      - {_i^k\mathbb{R}_{i_p}^{k_p}} \ 
     411      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }\right)\,\delta_{i+ i_p}[T^k] \right.\\ 
     412  - \left. { } _i^k{\mathbb{A}_w}_{i_p}^{k_p} 
     413    \left( 
     414      \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
     415      -{\:}_i^k\mathbb{R}_{i_p}^{k_p} 
     416      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
     417    \right) {\,}_i^k\mathbb{R}_{i_p}^{k_p}\delta_{k+ k_p}[T^i] 
     418  \right \}. 
    400419\end{multline*} 
    401420The key point is then that if we require $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$ and $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$ to 
    … …  
    431450where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the lateral and vertical fluxes/unit area 
    432451\[ 
    433 \mathbf{F}=\left( 
    434 \left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p}, 
    435 \left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p} 
    436  \right) 
     452  \mathbf{F}=\left( 
     453    \left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p}, 
     454    \left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p} 
     455  \right) 
    437456\] 
    438457and the gradient 
    439  \[\nabla T = \left( 
    440 \left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k}, 
    441 \left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p} 
    442 \right) 
     458\[ 
     459  \nabla T = \left( 
     460    \left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k}, 
     461    \left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p} 
     462  \right) 
    443463\] 
    444464 
    445465\subsection{Triad volumes in Griffes's scheme and in \NEMO} 
     466 
    446467To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. 
    447468\citet{Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells, 
    … …  
    460481\begin{equation} 
    461482  \label{eq:lat-normal} 
    462 -\overline\Alts_{\,i+1/2}^k\; 
    463 \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    464 \;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}} 
    465  = -\overline\Alts_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{\,i + 1/2}^{\,k}}. 
     483  -\overline\Alts_{\,i+1/2}^k\; 
     484  \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     485  \;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}} 
     486  = -\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}}. 
    466487\end{equation} 
    467488In fact if the diffusive coefficient is defined at $u$-points, 
    … …  
    471492 
    472493\subsection{Summary of the scheme} 
     494 
    473495The iso-neutral fluxes at $u$- and $w$-points are the sums of the triad fluxes that 
    474496cross the $u$- and $w$-faces \autoref{eq:iso_flux}: 
    475 \begin{subequations}\label{eq:alltriadflux} 
    476   \begin{flalign}\label{eq:vect_isoflux} 
     497\begin{subequations} 
     498  % \label{eq:alltriadflux} 
     499  \begin{flalign*} 
     500    % \label{eq:vect_isoflux} 
    477501    \vect{F}_{\mathrm{iso}}(T) &\equiv 
    478502    \sum_{\substack{i_p,\,k_p}} 
    479503    \begin{pmatrix} 
    480       {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\ 
    481       \\ 
    482       {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)  
     504      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T) \\ \\ 
     505      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T) 
    483506    \end{pmatrix}, 
    484   \end{flalign} 
     507  \end{flalign*} 
    485508  where \autoref{eq:latflux-triad}: 
    486509  \begin{align} 
    487510    \label{eq:triadfluxu} 
    488511    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{ 
    489       \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    490     \left( 
    491       \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    492       -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \ 
    493       \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    494     \right),\\ 
     512                                          \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     513                                          \left( 
     514                                          \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
     515                                          -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \ 
     516                                          \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
     517                                          \right),\\ 
    495518    \intertext{and} 
    496519    _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T) 
    497     &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}} 
    498     \left( 
    499       {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    500       -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 
    501       \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    502     \right),\label{eq:triadfluxw} 
     520                                        &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}} 
     521                                          \left( 
     522                                          {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
     523                                          -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \ 
     524                                          \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
     525                                          \right),\label{eq:triadfluxw} 
    503526  \end{align} 
    504527  with \autoref{eq:V-NEMO} 
    505   \begin{equation} 
    506     \label{eq:V-NEMO2} 
     528  \[ 
     529    % \label{eq:V-NEMO2} 
    507530    _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k. 
    508   \end{equation} 
     531  \] 
    509532\end{subequations} 
    510533 
    511534The divergence of the expression \autoref{eq:iso_flux} for the fluxes gives the iso-neutral diffusion tendency at 
    512535each tracer point: 
    513 \begin{equation} \label{eq:iso_operator} D_l^T = \frac{1}{b_T} 
     536\[ 
     537  % \label{eq:iso_operator} 
     538  D_l^T = \frac{1}{b_T} 
    514539  \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k 
    515540        {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[ 
    516541      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right] \right\} 
    517 \end{equation} 
     542\] 
    518543where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. 
    519544The diffusion scheme satisfies the following six properties: 
    … …  
    522547  The discretization of the diffusion operator recovers the traditional five-point Laplacian 
    523548  \autoref{eq:lat-normal} in the limit of flat iso-neutral direction: 
    524   \begin{equation} \label{eq:iso_property0} D_l^T = \frac{1}{b_T} \ 
     549  \[ 
     550    % \label{eq:iso_property0} 
     551    D_l^T = \frac{1}{b_T} \ 
    525552    \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \; 
    526553      \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad 
    527554    \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0 
    528   \end{equation} 
     555  \] 
    529556 
    530557\item[$\bullet$ implicit treatment in the vertical] 
    … …  
    534561  solve the vertical diffusion equation. 
    535562  This is necessary since the vertical eddy diffusivity associated with this term, 
    536   \begin{equation} 
     563  \[ 
    537564    \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{ 
    538565      {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 
    … …  
    541568      {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2 
    542569    \right\}, 
    543   \end{equation} 
     570  \] 
    544571  (where $b_w= e_{1w}\,e_{2w}\,e_{3w}$ is the volume of $w$-cells) can be quite large. 
    545572 
    … …  
    550577\item[$\bullet$ conservation of tracer] 
    551578  The iso-neutral diffusion conserves tracer content, $i.e.$ 
    552   \begin{equation} \label{eq:iso_property1} \sum_{i,j,k} \left\{ D_l^T \ 
    553       b_T \right\} = 0 
    554   \end{equation} 
     579  \[ 
     580    % \label{eq:iso_property1} 
     581    \sum_{i,j,k} \left\{ D_l^T \      b_T \right\} = 0 
     582  \] 
    555583  This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form. 
    556584 
    557585\item[$\bullet$ no increase of tracer variance] 
    558586  The iso-neutral diffusion does not increase the tracer variance, $i.e.$ 
    559   \begin{equation} \label{eq:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T 
    560       \ b_T \right\} \leq 0 
    561   \end{equation} 
     587  \[ 
     588    % \label{eq:iso_property2} 
     589    \sum_{i,j,k} \left\{ T \ D_l^T      \ b_T \right\} \leq 0 
     590  \] 
    562591  The property is demonstrated in \autoref{subsec:variance} above. 
    563592  It is a key property for a diffusion term. 
    … …  
    569598\item[$\bullet$ self-adjoint operator] 
    570599  The iso-neutral diffusion operator is self-adjoint, $i.e.$ 
    571   \begin{equation} \label{eq:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T 
    572       \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
     600  \begin{equation} 
     601    \label{eq:iso_property3} 
     602    \sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\} 
    573603  \end{equation} 
    574604  In other word, there is no need to develop a specific routine from the adjoint of this operator. 
    … …  
    578608  can be found by replacing $\delta[T]$ by $\delta[S]$ in \autoref{eq:dvar_iso_i} and \autoref{eq:dvar_iso_k}. 
    579609  This results in a term similar to \autoref{eq:perfect-square}, 
    580 \begin{equation} 
    581   \label{eq:TScovar} 
    582   - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
    583   \left( 
    584     \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    585     -{\:}_i^k\mathbb{R}_{i_p}^{k_p} 
    586     \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    587   \right) 
    588   \left( 
    589     \frac{ \delta_{i+ i_p}[S^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
    590     -{\:}_i^k\mathbb{R}_{i_p}^{k_p} 
    591     \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
    592   \right). 
    593 \end{equation} 
     610  \[ 
     611    % \label{eq:TScovar} 
     612    - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p} 
     613    \left( 
     614      \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
     615      -{\:}_i^k\mathbb{R}_{i_p}^{k_p} 
     616      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
     617    \right) 
     618    \left( 
     619      \frac{ \delta_{i+ i_p}[S^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} } 
     620      -{\:}_i^k\mathbb{R}_{i_p}^{k_p} 
     621      \frac{ \delta_{k+k_p} [S^i] }{{e_{3w}}_{\,i}^{\,k+k_p} } 
     622    \right). 
     623  \] 
    594624This is symmetrical in $T $ and $S$, so exactly the same term arises from 
    595625the discretization of this triad's contribution towards the RHS of \autoref{eq:iso_property3}. 
    596626\end{description} 
    597627 
    598 \subsection{Treatment of the triads at the boundaries}\label{sec:iso_bdry} 
     628\subsection{Treatment of the triads at the boundaries} 
     629\label{sec:iso_bdry} 
     630 
    599631The triad slope can only be defined where both the grid boxes centred at the end of the arms exist. 
    600632Triads that would poke up through the upper ocean surface into the atmosphere, 
    … …  
    618650For setups with topography without bbl mixing, \np{ln\_botmix\_triad}\forcode{ = .true.} may be necessary. 
    619651% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    620 \begin{figure}[h] \begin{center} 
     652\begin{figure}[h] 
     653  \begin{center} 
    621654    \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads} 
    622     \caption{  \protect\label{fig:bdry_triads} 
     655    \caption{ 
     656      \protect\label{fig:bdry_triads} 
    623657      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots), 
    624658      and $i+1/2,1$ $u$-point (blue square). 
    … …  
    627661      However, the lateral $_{11}$ contributions towards $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and 
    628662      $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$ (yellow line) are still applied, 
    629       giving diapycnal diffusive fluxes.\newline 
     663      giving diapycnal diffusive fluxes. 
     664      \newline 
    630665      (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and 
    631666      $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$ or $i+1,k+1$ tracer points is masked, 
    … …  
    633668      The associated lateral fluxes (grey-black dashed line) are masked if 
    634669      \protect\np{botmix\_triad}\forcode{ = .false.}, but left unmasked, 
    635       giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.}} 
    636  \end{center} \end{figure} 
     670      giving bottom mixing, if \protect\np{botmix\_triad}\forcode{ = .true.} 
     671    } 
     672  \end{center} 
     673\end{figure} 
    637674% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    638675 
    639 \subsection{ Limiting of the slopes within the interior}\label{sec:limit} 
     676\subsection{ Limiting of the slopes within the interior} 
     677\label{sec:limit} 
     678 
    640679As discussed in \autoref{subsec:LDF_slp_iso}, 
    641680iso-neutral slopes relative to geopotentials must be bounded everywhere, 
    … …  
    663702and so acts to reduce gravitational potential energy. 
    664703 
    665 \subsection{Tapering within the surface mixed layer}\label{sec:taper} 
     704\subsection{Tapering within the surface mixed layer} 
     705\label{sec:taper} 
     706 
    666707Additional tapering of the iso-neutral fluxes is necessary within the surface mixed layer. 
    667708When the Griffies triads are used, we offer two options for this. 
    668709 
    669 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:lintaper} 
     710\subsubsection{Linear slope tapering within the surface mixed layer} 
     711\label{sec:lintaper} 
     712 
    670713This is the option activated by the default choice \np{ln\_triad\_iso}\forcode{ = .false.}. 
    671714Slopes $\tilde{r}_i$ relative to geopotentials are tapered linearly from their value immediately below 
    672715the mixed layer to zero at the surface, as described in option (c) of \autoref{fig:eiv_slp}, to values 
    673 \begin{subequations} 
    674   \begin{equation} 
    675     \label{eq:rmtilde} 
    676     \rMLt = 
    677     -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h, 
    678   \end{equation} 
    679   and then the $r_i$ relative to vertical coordinate surfaces are appropriately adjusted to 
    680   \begin{equation} 
    681     \label{eq:rm} 
    682     \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h. 
    683   \end{equation} 
    684 \end{subequations} 
     716\begin{equation} 
     717  \label{eq:rmtilde} 
     718  \rMLt = -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h, 
     719\end{equation} 
     720and then the $r_i$ relative to vertical coordinate surfaces are appropriately adjusted to 
     721\[ 
     722  % \label{eq:rm} 
     723  \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h. 
     724\] 
    685725Thus the diffusion operator within the mixed layer is given by: 
    686 \begin{equation} \label{eq:iso_tensor_ML} 
    687 D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    688 \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 
    689  1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\ 
    690  0 \hfill & 1 \hfill & {-\rML[2]} \hfill \\ 
    691  {-\rML[1]}\hfill &   {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill 
    692 \end{array} }} \right) 
    693 \end{equation} 
     726\[ 
     727  % \label{eq:iso_tensor_ML} 
     728  D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
     729  \mbox{with}\quad \;\;\Re =\left( {{ 
     730        \begin{array}{*{20}c} 
     731          1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\ 
     732          0 \hfill & 1 \hfill & {-\rML[2]} \hfill \\ 
     733          {-\rML[1]}\hfill &   {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill 
     734        \end{array} 
     735      }} \right) 
     736\] 
    694737 
    695738This slope tapering gives a natural connection between tracer in the mixed-layer and 
    … …  
    726769  these basal triad slopes ${\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}$ are representative of the thermocline. 
    727770  The four basal triads defined in the bottom part of \autoref{fig:MLB_triad} are then 
    728 \begin{align} 
    729   {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &= 
    730  {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, \label{eq:Rbase} 
    731 \\ 
    732 \intertext{with e.g.\ the green triad} 
    733 {\:}_i{\mathbb{R}_{\mathrm{base}}}_{1/2}^{-1/2}&= 
    734 {\:}^{k_{\mathrm{ML}}}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}. \notag 
    735 \end{align} 
     771  \begin{align*} 
     772    {\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p} &= 
     773                                                       {\:}^{k_{\mathrm{ML}}-k_p-1/2}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}, 
     774                                                       % \label{eq:Rbase} 
     775    \\ 
     776    \intertext{with e.g.\ the green triad} 
     777    {\:}_i{\mathbb{R}_{\mathrm{base}}}_{1/2}^{-1/2}&= 
     778                                                     {\:}^{k_{\mathrm{ML}}}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2}. 
     779  \end{align*} 
    736780The vertical flux associated with each of these triads passes through 
    737781the $w$-point $i,k_{\mathrm{ML}}-1/2$ lying \emph{below} the $i,k_{\mathrm{ML}}$ tracer point, so it is this depth 
    738 \begin{equation} 
    739   \label{eq:zbase} 
     782\[ 
     783  % \label{eq:zbase} 
    740784  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2} 
    741 \end{equation} 
     785\] 
    742786one gridbox deeper than the diagnosed ML depth $z_{\mathrm{ML}})$ that sets the $h$ used to taper the slopes in 
    743787\autoref{eq:rmtilde}. 
    … …  
    747791  the ratio of the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_{\mathrm{base}}}_{\,i}$. 
    748792  For instance the green triad centred on $i,k$ 
    749 \begin{align} 
    750   {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,1/2}^{-1/2} &= 
    751 \frac{{z_w}_{k-1/2}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2} 
    752 \notag \\ 
    753 \intertext{and more generally} 
    754  {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &= 
    755 \frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}.\label{eq:RML} 
    756 \end{align} 
     793  \begin{align*} 
     794    {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,1/2}^{-1/2} &= 
     795                                                        \frac{{z_w}_{k-1/2}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,1/2}^{-1/2} \\ 
     796    \intertext{and more generally} 
     797    {\:}_i^k{\mathbb{R}_{\mathrm{ML}}}_{\,i_p}^{k_p} &= 
     798                                                       \frac{{z_w}_{k+k_p}}{{z_{\mathrm{base}}}_{\,i}}{\:}_i{\mathbb{R}_{\mathrm{base}}}_{\,i_p}^{k_p}. 
     799                                                       % \label{eq:RML} 
     800  \end{align*} 
    757801\end{enumerate} 
    758802 
    … …  
    760804\begin{figure}[h] 
    761805%  \fcapside { 
    762   \caption{\protect\label{fig:MLB_triad} 
     806  \caption{ 
     807    \protect\label{fig:MLB_triad} 
    763808    Definition of mixed-layer depth and calculation of linearly tapered triads. 
    764809    The figure shows a water column at a given $i,j$ (simplified to $i$), with the ocean surface at the top. 
    765810    Tracer points are denoted by bullets, and black lines the edges of the tracer cells; 
    766     $k$ increases upwards. \newline 
     811    $k$ increases upwards. 
     812    \newline 
    767813    \hspace{5 em} 
    768814    We define the mixed-layer by setting the vertical index of the tracer point immediately below the mixed layer, 
    … …  
    776822    Triads with different $i_p,k_p$, denoted by different colours, 
    777823    (e.g. the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.} 
    778 %} 
    779   {\includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads}} 
     824  % } 
     825  \includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads} 
    780826\end{figure} 
    781827% >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    783829\subsubsection{Additional truncation of skew iso-neutral flux components} 
    784830\label{subsec:Gerdes-taper} 
     831 
    785832The alternative option is activated by setting \np{ln\_triad\_iso} = true. 
    786833This retains the same tapered slope $\rML$  described above for the calculation of the $_{33}$ term of 
    … …  
    792839\end{equation} 
    793840giving a ML diffusive operator 
    794 \begin{equation} \label{eq:iso_tensor_ML2} 
    795 D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
    796 \mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c} 
    797  1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\ 
    798  0 \hfill & 1 \hfill & {-\rML[2]^*} \hfill \\ 
    799  {-\rML[1]^*}\hfill &   {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\ 
    800 \end{array} }} \right). 
    801 \end{equation} 
     841\[ 
     842  % \label{eq:iso_tensor_ML2} 
     843  D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad 
     844  \mbox{with}\quad \;\;\Re =\left( {{ 
     845        \begin{array}{*{20}c} 
     846          1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\ 
     847          0 \hfill & 1 \hfill & {-\rML[2]^*} \hfill \\ 
     848          {-\rML[1]^*}\hfill &   {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\ 
     849        \end{array} 
     850      }} \right). 
     851\] 
    802852This operator 
    803853\footnote{ 
    804854  To ensure good behaviour where horizontal density gradients are weak, 
    805855  we in fact follow \citet{Gerdes1991} and 
    806   set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.} 
     856  set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$. 
     857} 
    807858then has the property it gives no vertical density flux, and so does not change the potential energy. 
    808859This approach is similar to multiplying the iso-neutral diffusion coefficient by 
    … …  
    821872% Skew flux formulation for Eddy Induced Velocity : 
    822873% ================================================================ 
    823 \section{Eddy induced advection formulated as a skew flux}\label{sec:skew-flux} 
    824  
    825 \subsection{Continuous skew flux formulation}\label{sec:continuous-skew-flux} 
     874\section{Eddy induced advection formulated as a skew flux} 
     875\label{sec:skew-flux} 
     876 
     877\subsection{Continuous skew flux formulation} 
     878\label{sec:continuous-skew-flux} 
    826879 
    827880When Gent and McWilliams's [1990] diffusion is used, an additional advection term is added. 
    … …  
    833886 
    834887The eddy induced velocity is given by: 
    835 \begin{subequations} \label{eq:eiv} 
    836 \begin{equation}\label{eq:eiv_v} 
    837 \begin{split} 
    838  u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\ 
    839  v^* & = - \frac{1}{e_{3}}\;          \partial_j\psi_2,    \\ 
    840 w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i  \left( e_{2} \, \psi_1\right) 
    841                          + \partial_j  \left( e_{1} \, \psi_2\right) \right\}, 
    842 \end{split} 
    843 \end{equation} 
    844 where the streamfunctions $\psi_i$ are given by 
    845 \begin{equation} \label{eq:eiv_psi} 
    846 \begin{split} 
    847 \psi_1 & = A_{e} \; \tilde{r}_1,   \\ 
    848 \psi_2 & = A_{e} \; \tilde{r}_2, 
    849 \end{split} 
    850 \end{equation} 
     888\begin{subequations} 
     889  % \label{eq:eiv} 
     890  \begin{equation} 
     891    \label{eq:eiv_v} 
     892    \begin{split} 
     893      u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\ 
     894      v^* & = - \frac{1}{e_{3}}\;          \partial_j\psi_2,    \\ 
     895      w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i  \left( e_{2} \, \psi_1\right) 
     896        + \partial_j  \left( e_{1} \, \psi_2\right) \right\}, 
     897    \end{split} 
     898  \end{equation} 
     899  where the streamfunctions $\psi_i$ are given by 
     900  \begin{equation} 
     901    \label{eq:eiv_psi} 
     902    \begin{split} 
     903      \psi_1 & = A_{e} \; \tilde{r}_1,   \\ 
     904      \psi_2 & = A_{e} \; \tilde{r}_2, 
     905    \end{split} 
     906  \end{equation} 
    851907\end{subequations} 
    852908with $A_{e}$ the eddy induced velocity coefficient, 
    … …  
    868924the tracer advective fluxes per unit area in $ijk$ space can be transformed as follows: 
    869925\begin{flalign*} 
    870 \begin{split} 
    871 \textbf{F}_{\mathrm{eiv}}^T = 
    872 \begin{pmatrix} 
    873            {e_{2}\,e_{3}\;  u^*}       \\ 
    874       {e_{1}\,e_{2}\; w^*}  \\ 
    875 \end{pmatrix}   \;   T 
    876 &= 
    877 \begin{pmatrix} 
    878            { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;}     \\ 
    879       {+ \partial_i  \left( e_{2} \, \psi_1 \right) \; T \;}    \\ 
    880 \end{pmatrix}        \\ 
    881 &= 
    882 \begin{pmatrix} 
    883            { - \partial_k \left( e_{2} \, \psi_1  \; T \right) \;} \\ 
    884       {+ \partial_i  \left( e_{2} \,\psi_1 \; T \right) \;}  \\ 
    885 \end{pmatrix} 
    886  + 
    887 \begin{pmatrix} 
    888            {+ e_{2} \, \psi_1  \; \partial_k T} \\ 
    889       { - e_{2} \, \psi_1  \; \partial_i  T}  \\ 
    890 \end{pmatrix} 
    891 \end{split} 
     926  \begin{split} 
     927    \textbf{F}_{\mathrm{eiv}}^T = 
     928    \begin{pmatrix} 
     929      {e_{2}\,e_{3}\;  u^*} \\ 
     930      {e_{1}\,e_{2}\; w^*} 
     931    \end{pmatrix}   \;   T 
     932    &= 
     933    \begin{pmatrix} 
     934      { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;} \\ 
     935      {+ \partial_i  \left( e_{2} \, \psi_1 \right) \; T \;} 
     936    \end{pmatrix}          \\ 
     937    &= 
     938    \begin{pmatrix} 
     939      { - \partial_k \left( e_{2} \, \psi_1  \; T \right) \;} \\ 
     940      {+ \partial_i  \left( e_{2} \,\psi_1 \; T \right) \;} 
     941    \end{pmatrix} 
     942    + 
     943    \begin{pmatrix} 
     944      {+ e_{2} \, \psi_1  \; \partial_k T} \\ 
     945      { - e_{2} \, \psi_1  \; \partial_i  T} 
     946    \end{pmatrix} 
     947  \end{split} 
    892948\end{flalign*} 
    893949and since the eddy induced velocity field is non-divergent, 
    894950we end up with the skew form of the eddy induced advective fluxes per unit area in $ijk$ space: 
    895 \begin{equation} \label{eq:eiv_skew_ijk} 
    896 \textbf{F}_\mathrm{eiv}^T = \begin{pmatrix} 
    897            {+ e_{2} \, \psi_1  \; \partial_k T}   \\ 
    898       { - e_{2} \, \psi_1  \; \partial_i  T}  \\ 
    899                                  \end{pmatrix} 
     951\begin{equation} 
     952  \label{eq:eiv_skew_ijk} 
     953  \textbf{F}_\mathrm{eiv}^T = 
     954  \begin{pmatrix} 
     955    {+ e_{2} \, \psi_1  \; \partial_k T}   \\ 
     956    { - e_{2} \, \psi_1  \; \partial_i  T} 
     957  \end{pmatrix} 
    900958\end{equation} 
    901959The total fluxes per unit physical area are then 
    902 \begin{equation}\label{eq:eiv_skew_physical} 
    903 \begin{split} 
    904  f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\ 
    905  f^*_2 & = \frac{1}{e_{3}}\; \psi_2 \partial_k T   \\ 
    906  f^*_3 & =  -\frac{1}{e_{1}e_{2}}\; \left\{ e_{2} \psi_1 \partial_i T 
    907    + e_{1} \psi_2 \partial_j T \right\}. \\ 
     960\begin{equation} 
     961  \label{eq:eiv_skew_physical} 
     962  \begin{split} 
     963    f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\ 
     964    f^*_2 & = \frac{1}{e_{3}}\; \psi_2 \partial_k T   \\ 
     965    f^*_3 & =  -\frac{1}{e_{1}e_{2}}\; \left\{ e_{2} \psi_1 \partial_i T + e_{1} \psi_2 \partial_j T \right\}. 
    908966\end{split} 
    909967\end{equation} 
    … …  
    913971The tendency associated with eddy induced velocity is then simply the convergence of the fluxes 
    914972(\autoref{eq:eiv_skew_ijk}, \autoref{eq:eiv_skew_physical}), so 
    915 \begin{equation} \label{eq:skew_eiv_conv} 
    916 \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[ 
    917   \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right) 
    918   + \frac{\partial}{\partial j} \left( e_1  \; 
    919     \psi_2 \partial_k T\right) 
    920  -  \frac{\partial}{\partial k} \left( e_{2} \psi_1 \partial_i T 
    921    + e_{1} \psi_2 \partial_j T \right)  \right] 
    922 \end{equation} 
     973\[ 
     974  % \label{eq:skew_eiv_conv} 
     975  \frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[ 
     976    \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right) 
     977    + \frac{\partial}{\partial j} \left( e_1  \; 
     978      \psi_2 \partial_k T\right) 
     979    -  \frac{\partial}{\partial k} \left( e_{2} \psi_1 \partial_i T 
     980      + e_{1} \psi_2 \partial_j T \right)  \right] 
     981\] 
    923982It naturally conserves the tracer content, as it is expressed in flux form. 
    924983Since it has the same divergence as the advective form it also preserves the tracer variance. 
    925984 
    926985\subsection{Discrete skew flux formulation} 
     986 
    927987The skew fluxes in (\autoref{eq:eiv_skew_physical}, \autoref{eq:eiv_skew_ijk}), 
    928988like the off-diagonal terms (\autoref{eq:i13c}, \autoref{eq:i31c}) of the small angle diffusion tensor, 
    … …  
    934994defining $A_e$ at $T$-points is then given by: 
    935995 
    936  
    937 \begin{subequations}\label{eq:allskewflux} 
    938   \begin{flalign}\label{eq:vect_skew_flux} 
    939     \vect{F}_{\mathrm{eiv}}(T) &\equiv 
    940     \sum_{\substack{i_p,\,k_p}} 
     996\begin{subequations} 
     997  % \label{eq:allskewflux} 
     998  \begin{flalign*} 
     999    % \label{eq:vect_skew_flux} 
     1000    \vect{F}_{\mathrm{eiv}}(T) &\equiv    \sum_{\substack{i_p,\,k_p}} 
    9411001    \begin{pmatrix} 
    942       {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T)      \\ 
    943       \\ 
     1002      {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T)      \\      \\ 
    9441003      {_i^{k+1/2-k_p} {\mathbb{S}_w}_{i_p}^{k_p} } (T)      \\ 
    9451004    \end{pmatrix}, 
    946   \end{flalign} 
     1005  \end{flalign*} 
    9471006  where the skew flux in the $i$-direction associated with a given triad is (\autoref{eq:latflux-triad}, 
    9481007  \autoref{eq:triadfluxu}): 
    … …  
    9501009    \label{eq:skewfluxu} 
    9511010    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{ 
    952       \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    953      \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \ 
    954       \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, 
    955    \\ 
    956     \intertext{and \autoref{eq:triadfluxw} in the $k$-direction, changing the sign 
    957       to be consistent with \autoref{eq:eiv_skew_ijk}:} 
     1011                                          \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     1012                                          \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \ 
     1013                                          \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }, \\ 
     1014    \intertext{ 
     1015    and \autoref{eq:triadfluxw} in the $k$-direction, changing the sign 
     1016    to be consistent with \autoref{eq:eiv_skew_ijk}: 
     1017    } 
    9581018    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T) 
    959     &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} 
    960      {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:skewfluxw} 
     1019                                        &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}} 
     1020                                          {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:skewfluxw} 
    9611021  \end{align} 
    9621022\end{subequations} 
    … …  
    9661026 
    9671027\subsubsection{No change in tracer variance} 
     1028 
    9681029The discretization conserves tracer variance, $i.e.$ it does not include a diffusive component but is a `pure' advection term. 
    9691030This can be seen %either from Appendix \autoref{apdx:eiv_skew} or 
    … …  
    9731034summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of 
    9741035\begin{equation} 
    975 \label{eq:dvar_eiv_i} 
     1036  \label{eq:dvar_eiv_i} 
    9761037  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], 
    9771038\end{equation} 
    … …  
    9791040the $T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of 
    9801041\begin{equation} 
    981 \label{eq:dvar_eiv_k} 
     1042  \label{eq:dvar_eiv_k} 
    9821043  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i]. 
    9831044\end{equation} 
    … …  
    9871048 
    9881049\subsubsection{Reduction in gravitational PE} 
     1050 
    9891051The vertical density flux associated with the vertical skew-flux always has the same sign as 
    9901052the vertical density gradient; 
    … …  
    9991061    {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k 
    10001062    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\ 
    1001 \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from 
    1002   \autoref{eq:skewfluxw}, gives} 
    1003 % and separating out 
    1004 % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, 
    1005 % gives two terms. The 
    1006 % first $\rtriad{R}$ term (the only term for $z$-coordinates) is: 
    1007  &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} 
    1008 \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 \\ 
    1009  &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 
    1010      \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}} 
    1011 \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, 
     1063  \intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from \autoref{eq:skewfluxw}, gives} 
     1064  % and separating out 
     1065  % $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$, 
     1066  % gives two terms. The 
     1067  % first $\rtriad{R}$ term (the only term for $z$-coordinates) is: 
     1068  &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} 
     1069    \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 \\ 
     1070  &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 
     1071    \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}} 
     1072    \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}}, 
    10121073\end{align} 
    10131074using the definition of the triad slope $\rtriad{R}$, \autoref{eq:R} to 
    … …  
    10181079\begin{multline} 
    10191080  \label{eq:lat_densityPE} 
    1020  g \delta_{i+i_p}[z_T^k] 
    1021 \left[ 
    1022 -\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) 
    1023 \right] \\ 
    1024 = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 
    1025      \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} 
    1026 \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) 
    1027 \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}, 
     1081  g \delta_{i+i_p}[z_T^k] 
     1082  \left[ 
     1083    -\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) 
     1084  \right] \\ 
     1085  = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 
     1086  \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}} 
     1087  \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) 
     1088  \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}}, 
    10281089\end{multline} 
    10291090(using \autoref{eq:skewfluxu}) and so the total PE change \autoref{eq:vert_densityPE} + 
    10301091\autoref{eq:lat_densityPE} associated with the triad fluxes is 
    1031 \begin{multline} 
    1032   \label{eq:tot_densityPE} 
     1092\begin{multline*} 
     1093  % \label{eq:tot_densityPE} 
    10331094  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) + 
    1034 g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ 
    1035 = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 
    1036      \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 
    1037 \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}}. 
    1038 \end{multline} 
     1095  g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\ 
     1096  = +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k 
     1097  \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)^2 
     1098  \frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}. 
     1099\end{multline*} 
    10391100Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+ 
    10401101\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative. 
    10411102 
    1042 \subsection{Treatment of the triads at the boundaries}\label{sec:skew_bdry} 
     1103\subsection{Treatment of the triads at the boundaries} 
     1104\label{sec:skew_bdry} 
     1105 
    10431106Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes are masked at the boundaries  
    10441107in exactly the same way as are the triad slopes \rtriad{R} used for the iso-neutral diffusive fluxes,  
    … …  
    10491112The namelist parameter \np{ln\_botmix\_triad} has no effect on the eddy-induced skew-fluxes. 
    10501113 
    1051 \subsection{Limiting of the slopes within the interior}\label{sec:limitskew} 
     1114\subsection{Limiting of the slopes within the interior} 
     1115\label{sec:limitskew} 
     1116 
    10521117Presently, the iso-neutral slopes $\tilde{r}_i$ relative to geopotentials are limited to be less than $1/100$,  
    10531118exactly as in calculating the iso-neutral diffusion, \S \autoref{sec:limit}.  
    10541119Each individual triad \rtriadt{R} is so limited. 
    10551120 
    1056 \subsection{Tapering within the surface mixed layer}\label{sec:taperskew} 
     1121\subsection{Tapering within the surface mixed layer} 
     1122\label{sec:taperskew} 
     1123 
    10571124The slopes $\tilde{r}_i$ relative to geopotentials (and thus the individual triads \rtriadt{R})  
    10581125are always tapered linearly from their value immediately below the mixed layer to zero at the surface  
    … …  
    10721139(the horizontal flux convergence is relatively insignificant within the mixed-layer). 
    10731140 
    1074 \subsection{Streamfunction diagnostics}\label{sec:sfdiag} 
     1141\subsection{Streamfunction diagnostics} 
     1142\label{sec:sfdiag} 
     1143 
    10751144Where the namelist parameter \np{ln\_traldf\_gdia}\forcode{ = .true.}, 
    10761145diagnosed mean eddy-induced velocities are output. 
    … …  
    10801149We follow \citep{Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from 
    10811150the surrounding four triads according to: 
    1082 \begin{equation} 
    1083   \label{eq:sfdiagi} 
     1151\[ 
     1152  % \label{eq:sfdiagi} 
    10841153  {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}} 
    10851154  {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}. 
    1086 \end{equation} 
     1155\] 
    10871156The streamfunction $\psi_1$ is calculated similarly at $vw$ points. 
    10881157The eddy-induced velocities are then calculated from the straightforward discretisation of \autoref{eq:eiv_v}: 
    1089 \begin{equation}\label{eq:eiv_v_discrete} 
    1090 \begin{split} 
    1091  {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),   \\ 
    1092  {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),   \\ 
    1093  {w^*}_{i,j}^{k+1/2} & =    \frac{1}{e_{1t}e_{2t}}\; \left\{ 
    1094  {e_{2u}}_{i+1/2}^{k+1/2} \,{\psi_1}_{i+1/2}^{k+1/2} - 
    1095  {e_{2u}}_{i-1/2}^{k+1/2} \,{\psi_1}_{i-1/2}^{k+1/2} \right. + \\ 
    1096 \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\}, 
    1097 \end{split} 
    1098 \end{equation} 
     1158\[ 
     1159  % \label{eq:eiv_v_discrete} 
     1160  \begin{split} 
     1161    {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),   \\ 
     1162    {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),   \\ 
     1163    {w^*}_{i,j}^{k+1/2} & =    \frac{1}{e_{1t}e_{2t}}\; \left\{ 
     1164      {e_{2u}}_{i+1/2}^{k+1/2} \,{\psi_1}_{i+1/2}^{k+1/2} - 
     1165      {e_{2u}}_{i-1/2}^{k+1/2} \,{\psi_1}_{i-1/2}^{k+1/2} \right. + \\ 
     1166    \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\}, 
     1167  \end{split} 
     1168\] 
     1169 
     1170\biblio 
     1171 
    10991172\end{document} 
Note: See TracChangeset for help on using the changeset viewer.