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 10368 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_E.tex – NEMO

Ignore:
Timestamp:
2018-12-03T12:45:01+01:00 (5 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10365, see #2133

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/annex_E.tex

    r9407 r10368  
    1010\newpage 
    1111$\ $\newline    % force a new ligne 
    12   
    13  This appendix some on going consideration on algorithms used or planned to be used 
    14 in \NEMO.  
     12 
     13This appendix some on going consideration on algorithms used or planned to be used in \NEMO.  
    1514 
    1615$\ $\newline    % force a new ligne 
    … …  
    2221\label{sec:TRA_adv_ubs} 
    2322 
    24 The UBS advection scheme is an upstream biased third order scheme based on  
    25 an upstream-biased parabolic interpolation. It is also known as Cell Averaged  
    26 QUICK scheme (Quadratic Upstream Interpolation for Convective  
    27 Kinematics). For example, in the $i$-direction : 
     23The UBS advection scheme is an upstream biased third order scheme based on 
     24an upstream-biased parabolic interpolation. 
     25It is also known as Cell Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective Kinematics). 
     26For example, in the $i$-direction: 
    2827\begin{equation} \label{eq:tra_adv_ubs2} 
    2928\tau _u^{ubs} = \left\{  \begin{aligned} 
    … …  
    3837- \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] 
    3938\end{equation} 
    40 where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and  
    41 $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$.  
    42 By choosing this expression for $\tau "$ we consider a fourth order approximation  
    43 of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$).  
     39where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and 
     40$\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. 
     41By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with 
     42a constant i-grid spacing ($\Delta i=1$). 
    4443 
    4544Alternative choice: introduce the scale factors:   
    … …  
    4746 
    4847 
    49 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation  
    50 error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the  
    51 advection scheme is similar to that reported in \cite{Farrow1995}.  
    52 It is a relatively good compromise between accuracy and smoothness. It is  
    53 not a \emph{positive} scheme meaning false extrema are permitted but the  
    54 amplitude of such are significantly reduced over the centred second order  
    55 method. Nevertheless it is not recommended to apply it to a passive tracer  
    56 that requires positivity.  
    57  
    58 The intrinsic diffusion of UBS makes its use risky in the vertical direction  
    59 where the control of artificial diapycnal fluxes is of paramount importance.  
    60 It has therefore been preferred to evaluate the vertical flux using the TVD  
    61 scheme when \np{ln\_traadv\_ubs}\forcode{ = .true.}. 
    62  
    63 For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds  
    64 to a second order centred scheme is evaluated using the \textit{now} velocity  
    65 (centred in time) while the second term which is the diffusive part of the scheme,  
    66 is evaluated using the \textit{before} velocity (forward in time. This is discussed  
    67 by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. UBS and QUICK  
    68 schemes only differ by one coefficient. Substituting 1/6 with 1/8 in  
    69 (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.  
    70 This option is not available through a namelist parameter, since the 1/6  
    71 coefficient is hard coded. Nevertheless it is quite easy to make the  
    72 substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme 
    73  
    74 NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can  
    75 be controlled by vertical advection (not vertical diffusion which is usually  
    76 solved using an implicit scheme). Computer time can be saved by using a  
    77 time-splitting technique on vertical advection. This possibility have been  
    78 implemented and validated in ORCA05-L301. It is not currently offered in the  
    79 current reference version.  
    80  
    81 NB 2 : In a forthcoming release four options will be proposed for the  
    82 vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be  
    83 evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme ,  
    84 or  \textit{(b)} a TVD scheme, or  \textit{(c)} an interpolation based on conservative  
    85 parabolic splines following \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS,  
    86 or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an  
    87 eight-order accurate conventional scheme. 
    88  
    89 NB 3 : It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: 
     48This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error 
     49\citep{Shchepetkin_McWilliams_OM05}. 
     50The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. 
     51It is a relatively good compromise between accuracy and smoothness. 
     52It is not a \emph{positive} scheme meaning false extrema are permitted but 
     53the amplitude of such are significantly reduced over the centred second order method. 
     54Nevertheless it is not recommended to apply it to a passive tracer that requires positivity.  
     55 
     56The intrinsic diffusion of UBS makes its use risky in the vertical direction where 
     57the control of artificial diapycnal fluxes is of paramount importance. 
     58It has therefore been preferred to evaluate the vertical flux using the TVD scheme when 
     59\np{ln\_traadv\_ubs}\forcode{ = .true.}. 
     60 
     61For stability reasons, in \autoref{eq:tra_adv_ubs}, the first term which corresponds to 
     62a second order centred scheme is evaluated using the \textit{now} velocity (centred in time) while 
     63the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity 
     64(forward in time). 
     65This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. 
     66UBS and QUICK schemes only differ by one coefficient. 
     67Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
     68This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 
     69Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme. 
     70 
     71NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can be controlled by vertical advection 
     72(not vertical diffusion which is usually solved using an implicit scheme). 
     73Computer time can be saved by using a time-splitting technique on vertical advection. 
     74This possibility have been implemented and validated in ORCA05-L301. 
     75It is not currently offered in the current reference version.  
     76 
     77NB 2: In a forthcoming release four options will be proposed for the vertical component used in the UBS scheme. 
     78$\tau _w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme, 
     79or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following 
     80\citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS. 
     81The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme. 
     82 
     83NB 3: It is straight forward to rewrite \autoref{eq:tra_adv_ubs} as follows: 
    9084\begin{equation} \label{eq:tra_adv_ubs2} 
    9185\tau _u^{ubs} = \left\{  \begin{aligned} 
    … …  
    10296\end{split} 
    10397\end{equation} 
    104 \autoref{eq:tra_adv_ubs2} has several advantages. First it clearly evidence that  
    105 the UBS scheme is based on the fourth order scheme to which is added an  
    106 upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order  
    107 part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order  
    108 part as stated above using \autoref{eq:tra_adv_ubs}. Third, the diffusive term is  
    109 in fact a biharmonic operator with a eddy coefficient with is simply proportional  
    110 to the velocity. 
     98\autoref{eq:tra_adv_ubs2} has several advantages. 
     99First it clearly evidences that the UBS scheme is based on the fourth order scheme to which 
     100is added an upstream biased diffusive term. 
     101Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, 
     102not only the $2^{th}$ order part as stated above using \autoref{eq:tra_adv_ubs}. 
     103Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient which 
     104is simply proportional to the velocity. 
    111105 
    112106laplacian diffusion: 
    … …  
    135129with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$,  
    136130$i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$ 
    137 it comes : 
     131it comes: 
    138132\begin{equation} \label{eq:tra_ldf_lap} 
    139133\begin{split} 
    … …  
    163157\end{split} 
    164158\end{equation} 
    165 if the velocity is uniform ($i.e.$ $|u|=cst$) and choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right]$ 
     159if the velocity is uniform ($i.e.$ $|u|=cst$) and 
     160choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right]$ 
    166161 
    167162sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$): 
    … …  
    191186\label{sec:LF} 
    192187 
    193 We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are: 
     188We adopt the following semi-discrete notation for time derivative. 
     189Given the values of a variable $q$ at successive time step, 
     190the time derivation and averaging operators at the mid time step are: 
    194191\begin{subequations} \label{eq:dt_mt} 
    195192\begin{align} 
    … …  
    198195\end{align} 
    199196\end{subequations} 
    200 As for space operator, the adjoint of the derivation and averaging time operators are  
    201 $\delta_t^*=\delta_{t+\rdt/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ 
    202 , respectively.  
     197As for space operator, 
     198the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\rdt/2}$ and 
     199$\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$, respectively. 
    203200 
    204201The Leap-frog time stepping given by \autoref{eq:DOM_nxt} can be defined as: 
    … …  
    208205      =         \frac{q^{t+\rdt}-q^{t-\rdt}}{2\rdt} 
    209206\end{equation}  
    210 Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, not $2\rdt$  
    211 as it can be found sometime in literature.  
    212 The leap-Frog time stepping is a second order centered scheme. As such it respects  
    213 the quadratic invariant in integral forms, $i.e.$ the following continuous property, 
     207Note that \autoref{chap:LF} shows that the leapfrog time step is $\rdt$, 
     208not $2\rdt$ as it can be found sometimes in literature. 
     209The leap-Frog time stepping is a second order centered scheme. 
     210As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property, 
    214211\begin{equation} \label{eq:Energy} 
    215212\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}  
    … …  
    217214   =  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) , 
    218215\end{equation} 
    219 is satisfied in discrete form. Indeed,  
     216is satisfied in discrete form. 
     217Indeed,  
    220218\begin{equation} \begin{split} 
    221219\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt}  
    … …  
    228226      \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) 
    229227\end{split} \end{equation} 
    230 NB here pb of boundary condition when applying the adjoin! In space, setting to 0  
    231 the quantity in land area is sufficient to get rid of the boundary condition  
    232 (equivalently of the boundary value of the integration by part). In time this boundary  
    233 condition is not physical and \textbf{add something here!!!} 
     228NB here pb of boundary condition when applying the adjoint! 
     229In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition  
     230(equivalently of the boundary value of the integration by part). 
     231In time this boundary condition is not physical and \textbf{add something here!!!} 
    234232 
    235233 
    … …  
    249247\subsection{Griffies iso-neutral diffusion operator} 
    250248 
    251 Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} 
    252 scheme, but is formulated within the \NEMO framework ($i.e.$ using scale  
    253 factors rather than grid-size and having a position of $T$-points that is not  
    254 necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). 
    255  
    256 In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO,  
    257 the off-diagonal terms of the small angle diffusion tensor contain several double  
    258 spatial averages of a gradient, for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$.  
    259 It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer  
    260 allows for the presence of grid point oscillation structures that will be invisible  
    261 to the operator. These structures are \textit{computational modes}. They 
    262 will not be damped by the iso-neutral operator, and even possibly amplified by it.  
    263 In other word, the operator applied to a tracer does not warranties the decrease of  
    264 its global average variance. To circumvent this, we have introduced a smoothing of  
    265 the slopes of the iso-neutral surfaces (see \autoref{chap:LDF}). Nevertheless, this technique  
    266 works fine for $T$ and $S$ as they are active tracers ($i.e.$ they enter the computation  
    267 of density), but it does not work for a passive tracer.   \citep{Griffies_al_JPO98} introduce  
    268 a different way to discretise the off-diagonal terms that nicely solve the problem.  
    269 The idea is to get rid of combinations of an averaged in one direction combined  
    270 with a derivative in the same direction by considering triads. For example in the  
    271 (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows: 
     249Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme, 
     250but is formulated within the \NEMO framework 
     251($i.e.$ using scale factors rather than grid-size and having a position of $T$-points that 
     252is not necessary in the middle of vertical velocity points, see \autoref{fig:zgr_e3}). 
     253 
     254In the formulation \autoref{eq:tra_ldf_iso} introduced in 1995 in OPA, the ancestor of \NEMO, 
     255the off-diagonal terms of the small angle diffusion tensor contain several double spatial averages of a gradient, 
     256for example $\overline{\overline{\delta_k \cdot}}^{\,i,k}$. 
     257It is apparent that the combination of a $k$ average and a $k$ derivative of the tracer allows for 
     258the presence of grid point oscillation structures that will be invisible to the operator. 
     259These structures are \textit{computational modes}. 
     260They will not be damped by the iso-neutral operator, and even possibly amplified by it. 
     261In other word, the operator applied to a tracer does not warranties the decrease of its global average variance. 
     262To circumvent this, we have introduced a smoothing of the slopes of the iso-neutral surfaces 
     263(see \autoref{chap:LDF}). 
     264Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers 
     265($i.e.$ they enter the computation of density), but it does not work for a passive tracer. 
     266\citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that 
     267nicely solve the problem. 
     268The idea is to get rid of combinations of an averaged in one direction combined with 
     269a derivative in the same direction by considering triads. 
     270For example in the (\textbf{i},\textbf{k}) plane, the four triads are defined at the $(i,k)$ $T$-point as follows: 
    272271\begin{equation} \label{eq:Gf_triads} 
    273272_i^k \mathbb{T}_{i_p}^{k_p} (T) 
    … …  
    277276                             \right) 
    278277\end{equation} 
    279 where the indices $i_p$ and $k_p$ define the four triads and take the following value:  
    280 $i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$,  
    281 $b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells,  
     278where the indices $i_p$ and $k_p$ define the four triads and take the following value: 
     279$i_p = -1/2$ or $1/2$ and $k_p = -1/2$ or $1/2$, 
     280$b_u= e_{1u}\,e_{2u}\,e_{3u}$ is the volume of $u$-cells, 
    282281$A_i^k$ is the lateral eddy diffusivity coefficient defined at $T$-point, 
    283 and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad : 
     282and $_i^k \mathbb{R}_{i_p}^{k_p}$ is the slope associated with each triad: 
    284283\begin{equation} \label{eq:Gf_slopes} 
    285284_i^k \mathbb{R}_{i_p}^{k_p}  
    … …  
    288287{\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] } 
    289288\end{equation} 
    290 Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of  
    291 multiplying the temperature derivative by $\alpha$ and the salinity derivative  
    292 by $\beta$. This is more efficient as the ratio $\alpha / \beta$ can to be  
    293 evaluated directly. 
    294  
    295 Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of  
    296 ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. This choice has been motivated by the decrease  
    297 of tracer variance and the presence of partial cell at the ocean bottom  
    298 (see \autoref{apdx:Gf_operator}). 
     289Note that in \autoref{eq:Gf_slopes} we use the ratio $\alpha / \beta$ instead of 
     290multiplying the temperature derivative by $\alpha$ and the salinity derivative by $\beta$. 
     291This is more efficient as the ratio $\alpha / \beta$ can to be evaluated directly. 
     292 
     293Note that in \autoref{eq:Gf_triads}, we chose to use ${b_u}_{\,i+i_p}^{\,k}$ instead of ${b_{uw}}_{\,i+i_p}^{\,k+k_p}$. 
     294This choice has been motivated by the decrease of tracer variance and 
     295the presence of partial cell at the ocean bottom (see \autoref{apdx:Gf_operator}). 
    299296 
    300297%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    301298\begin{figure}[!ht] \begin{center} 
    302299\includegraphics[width=0.70\textwidth]{Fig_ISO_triad} 
    303 \caption{  \protect\label{fig:ISO_triad}    
    304 Triads used in the Griffies's like iso-neutral diffision scheme for  
    305 $u$-component (upper panel) and $w$-component (lower panel).} 
     300\caption{  \protect\label{fig:ISO_triad} 
     301  Triads used in the Griffies's like iso-neutral diffision scheme for 
     302  $u$-component (upper panel) and $w$-component (lower panel).} 
    306303\end{center} 
    307304\end{figure} 
    … …  
    309306 
    310307The four iso-neutral fluxes associated with the triads are defined at $T$-point.  
    311 They take the following expression : 
     308They take the following expression: 
    312309\begin{flalign} \label{eq:Gf_fluxes} 
    313310\begin{split} 
    … …  
    320317\end{flalign} 
    321318 
    322 The resulting iso-neutral fluxes at $u$- and $w$-points are then given by the  
    323 sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): 
     319The resulting iso-neutral fluxes at $u$- and $w$-points are then given by 
     320the sum of the fluxes that cross the $u$- and $w$-face (\autoref{fig:ISO_triad}): 
    324321\begin{flalign} \label{eq:iso_flux}  
    325322\textbf{F}_{iso}(T)  
    … …  
    350347%    \end{pmatrix}       
    351348\end{flalign} 
    352 resulting in a iso-neutral diffusion tendency on temperature given by the divergence  
    353 of the sum of all the four triad fluxes : 
     349resulting in a iso-neutral diffusion tendency on temperature given by 
     350the divergence of the sum of all the four triad fluxes: 
    354351\begin{equation} \label{eq:Gf_operator} 
    355352D_l^T = \frac{1}{b_T}  \sum_{\substack{i_p,\,k_p}} \left\{   
    … …  
    359356where $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells.  
    360357 
    361 This expression of the iso-neutral diffusion has been chosen in order to satisfy  
    362 the following six properties: 
     358This expression of the iso-neutral diffusion has been chosen in order to satisfy the following six properties: 
    363359\begin{description} 
    364 \item[$\bullet$ horizontal diffusion] The discretization of the diffusion operator  
    365 recovers the traditional five-point Laplacian in the limit of flat iso-neutral direction : 
     360\item[$\bullet$ horizontal diffusion] 
     361  The discretization of the diffusion operator recovers the traditional five-point Laplacian in 
     362  the limit of flat iso-neutral direction: 
    366363\begin{equation} \label{eq:Gf_property1a} 
    367364D_l^T = \frac{1}{b_T}  \ \delta_{i}  
    … …  
    371368\end{equation} 
    372369 
    373 \item[$\bullet$ implicit treatment in the vertical]  In the diagonal term associated  
    374 with the vertical divergence of the iso-neutral fluxes (i.e. the term associated  
    375 with a second order vertical derivative) appears only tracer values associated  
    376 with a single water column. This is of paramount importance since it means 
    377 that the implicit in time algorithm for solving the vertical diffusion equation can  
    378 be used to evaluate this term. It is a necessity since the vertical eddy diffusivity  
    379 associated with this term,   
     370\item[$\bullet$ implicit treatment in the vertical] 
     371  In the diagonal term associated with the vertical divergence of the iso-neutral fluxes 
     372  (i.e. the term associated with a second order vertical derivative) 
     373  appears only tracer values associated with a single water column. 
     374  This is of paramount importance since it means that 
     375  the implicit in time algorithm for solving the vertical diffusion equation can be used to evaluate this term. 
     376  It is a necessity since the vertical eddy diffusivity associated with this term,   
    380377\begin{equation} 
    381378    \sum_{\substack{i_p, \,k_p}} \left\{   
    … …  
    385382can be quite large. 
    386383 
    387 \item[$\bullet$ pure iso-neutral operator]  The iso-neutral flux of locally referenced  
    388 potential density is zero, $i.e.$ 
     384\item[$\bullet$ pure iso-neutral operator] 
     385  The iso-neutral flux of locally referenced potential density is zero, $i.e.$ 
    389386\begin{align} \label{eq:Gf_property2} 
    390387\begin{matrix} 
    … …  
    397394\end{matrix} 
    398395\end{align} 
    399 This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$  
    400 and the definition of the triads' slopes \autoref{eq:Gf_slopes}. 
    401  
    402 \item[$\bullet$ conservation of tracer] The iso-neutral diffusion term conserve the  
    403 total tracer content, $i.e.$ 
     396This result is trivially obtained using the \autoref{eq:Gf_triads} applied to $T$ and $S$ and 
     397the definition of the triads' slopes \autoref{eq:Gf_slopes}. 
     398 
     399\item[$\bullet$ conservation of tracer] 
     400  The iso-neutral diffusion term conserve the total tracer content, $i.e.$ 
    404401\begin{equation} \label{eq:Gf_property1} 
    405402\sum_{i,j,k} \left\{ D_l^T \ b_T \right\} = 0 
    406403\end{equation} 
    407 This property is trivially satisfied since the iso-neutral diffusive operator  
    408 is written in flux form. 
    409  
    410 \item[$\bullet$ decrease of tracer variance] The iso-neutral diffusion term does  
    411 not increase the total tracer variance, $i.e.$ 
     404This property is trivially satisfied since the iso-neutral diffusive operator is written in flux form. 
     405 
     406\item[$\bullet$ decrease of tracer variance] 
     407  The iso-neutral diffusion term does not increase the total tracer variance, $i.e.$ 
    412408\begin{equation} \label{eq:Gf_property1} 
    413409\sum_{i,j,k} \left\{ T \ D_l^T \ b_T \right\} \leq 0 
    414410\end{equation} 
    415 The property is demonstrated in the \autoref{apdx:Gf_operator}. It is a  
    416 key property for a diffusion term. It means that the operator is also a dissipation  
    417 term, $i.e.$ it is a sink term for the square of the quantity on which it is applied.  
    418 It therfore ensure that, when the diffusivity coefficient is large enough, the field  
    419 on which it is applied become free of grid-point noise. 
    420  
    421 \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion operator is self-adjoint,  
    422 $i.e.$ 
     411The property is demonstrated in the \autoref{apdx:Gf_operator}. 
     412It is a key property for a diffusion term. 
     413It means that the operator is also a dissipation term, 
     414$i.e.$ it is a sink term for the square of the quantity on which it is applied. 
     415It therfore ensures that, when the diffusivity coefficient is large enough, 
     416the field on which it is applied become free of grid-point noise. 
     417 
     418\item[$\bullet$ self-adjoint operator] 
     419  The iso-neutral diffusion operator is self-adjoint, $i.e.$ 
    423420\begin{equation} \label{eq:Gf_property1} 
    424421\sum_{i,j,k} \left\{ S \ D_l^T \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}  
    425422\end{equation} 
    426 In other word, there is no needs to develop a specific routine from the adjoint of this  
    427 operator. We just have to apply the same routine. This properties can be demonstrated  
    428 quite easily in a similar way the "non increase of tracer variance" property has been  
    429 proved (see \autoref{apdx:Gf_operator}). 
     423In other word, there is no needs to develop a specific routine from the adjoint of this operator. 
     424We just have to apply the same routine. 
     425This properties can be demonstrated quite easily in a similar way the "non increase of tracer variance" property 
     426has been proved (see \autoref{apdx:Gf_operator}). 
    430427\end{description} 
    431428 
    … …  
    437434\subsection{Eddy induced velocity and skew flux formulation} 
    438435 
    439 When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined),  
    440 an additional advection term is added. The associated velocity is the so called  
    441 eddy induced velocity, the formulation of which depends on the slopes of iso- 
    442 neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used  
    443 here are referenced to the geopotential surfaces, $i.e.$ \autoref{eq:ldfslp_geo}  
    444 is used in $z$-coordinate, and the sum \autoref{eq:ldfslp_geo} 
    445 + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.  
     436When Gent and McWilliams [1990] diffusion is used (\key{traldf\_eiv} defined), 
     437an additional advection term is added. 
     438The associated velocity is the so called eddy induced velocity, 
     439the formulation of which depends on the slopes of iso-neutral surfaces. 
     440Contrary to the case of iso-neutral mixing, the slopes used here are referenced to the geopotential surfaces, 
     441$i.e.$ \autoref{eq:ldfslp_geo} is used in $z$-coordinate, 
     442and the sum \autoref{eq:ldfslp_geo} + \autoref{eq:ldfslp_iso} in $z^*$ or $s$-coordinates.  
    446443 
    447444The eddy induced velocity is given by:  
    … …  
    456453\end{split} 
    457454\end{equation} 
    458 where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the  
    459 slopes between the iso-neutral and the geopotential surfaces. 
     455where $A_{e}$ is the eddy induced velocity coefficient, 
     456and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces. 
    460457%%gm wrong: to be modified with 2 2D streamfunctions 
    461  In other words, 
    462 the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, which 
    463 is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$ 
     458In other words, the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, 
     459which is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$. 
    464460%%end gm 
    465461 
    466 A traditional way to implement this additional advection is to add it to the eulerian  
    467 velocity prior to compute the tracer advection. This allows us to take advantage of  
    468 all the advection schemes offered for the tracers (see \autoref{sec:TRA_adv}) and not just  
    469 a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers  
    470 where \emph{positivity} of the advection scheme is of paramount importance.  
     462A traditional way to implement this additional advection is to add it to the eulerian velocity prior to 
     463compute the tracer advection. 
     464This allows us to take advantage of all the advection schemes offered for the tracers 
     465(see \autoref{sec:TRA_adv}) and not just a $2^{nd}$ order advection scheme. 
     466This is particularly useful for passive tracers where 
     467\emph{positivity} of the advection scheme is of paramount importance.  
    471468% give here the expression using the triads. It is different from the one given in \autoref{eq:ldfeiv} 
    472469% see just below a copy of this equation: 
    … …  
    490487\end{equation} 
    491488 
    492 \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection,  
    493 the so-called skew form. It is based on a transformation of the advective fluxes  
    494 using the non-divergent nature of the eddy induced velocity.  
    495 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be  
    496 transformed as follows: 
     489\citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form. 
     490It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity. 
     491For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows: 
    497492\begin{flalign*} 
    498493\begin{split} 
    … …  
    519514\end{split} 
    520515\end{flalign*} 
    521 and since the eddy induces velocity field is no-divergent, we end up with the skew  
    522 form of the eddy induced advective fluxes: 
     516and since the eddy induces velocity field is no-divergent, 
     517we end up with the skew form of the eddy induced advective fluxes: 
    523518\begin{equation} \label{eq:eiv_skew_continuous} 
    524519\textbf{F}_{eiv}^T = \begin{pmatrix}  
    … …  
    527522                                 \end{pmatrix} 
    528523\end{equation} 
    529 The tendency associated with eddy induced velocity is then simply the divergence  
    530 of the \autoref{eq:eiv_skew_continuous} fluxes. It naturally conserves the tracer  
    531 content, as it is expressed in flux form and, as the advective form, it preserve the  
    532 tracer variance. Another interesting property of \autoref{eq:eiv_skew_continuous}  
    533 form is that when $A=A_e$, a simplification occurs in the sum of the iso-neutral  
    534 diffusion and eddy induced velocity terms: 
     524The tendency associated with eddy induced velocity is then simply the divergence of 
     525the \autoref{eq:eiv_skew_continuous} fluxes. 
     526It naturally conserves the tracer content, as it is expressed in flux form and, 
     527as the advective form, it preserves the tracer variance. 
     528Another interesting property of \autoref{eq:eiv_skew_continuous} form is that when $A=A_e$, 
     529a simplification occurs in the sum of the iso-neutral diffusion and eddy induced velocity terms: 
    535530\begin{flalign} \label{eq:eiv_skew+eiv_continuous} 
    536531\textbf{F}_{iso}^T + \textbf{F}_{eiv}^T &=  
    … …  
    549544\end{pmatrix} 
    550545\end{flalign} 
    551 The horizontal component reduces to the one use for an horizontal laplacian  
    552 operator and the vertical one keep the same complexity, but not more. This property 
    553 has been used to reduce the computational time \citep{Griffies_JPO98}, but it is  
    554 not of practical use as usually $A \neq A_e$. Nevertheless this property can be used to  
    555 choose a discret form of  \autoref{eq:eiv_skew_continuous} which is consistent with the  
    556 iso-neutral operator \autoref{eq:Gf_operator}. Using the slopes \autoref{eq:Gf_slopes}  
    557 and defining $A_e$ at $T$-point($i.e.$ as $A$, the eddy diffusivity coefficient), 
    558 the resulting discret form is given by: 
     546The horizontal component reduces to the one use for an horizontal laplacian operator and 
     547the vertical one keeps the same complexity, but not more. 
     548This property has been used to reduce the computational time \citep{Griffies_JPO98}, 
     549but it is not of practical use as usually $A \neq A_e$. 
     550Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which 
     551is consistent with the iso-neutral operator \autoref{eq:Gf_operator}. 
     552Using the slopes \autoref{eq:Gf_slopes} and defining $A_e$ at $T$-point($i.e.$ as $A$, 
     553the eddy diffusivity coefficient), the resulting discret form is given by: 
    559554\begin{equation} \label{eq:eiv_skew}   
    560555\textbf{F}_{eiv}^T   \equiv   \frac{1}{4} \left( \begin{aligned}                                 
    … …  
    569564\end{equation} 
    570565Note that \autoref{eq:eiv_skew} is valid in $z$-coordinate with or without partial cells.  
    571 In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces  
    572 must be added to $\mathbb{R}$ for the discret form to be exact.  
    573  
    574 Such a choice of discretisation is consistent with the iso-neutral operator as it uses the  
    575 same definition for the slopes. It also ensures the conservation of the tracer variance  
    576 (see Appendix \autoref{apdx:eiv_skew}), $i.e.$ it does not include a diffusive component  
    577 but is a "pure" advection term. 
     566In $z^*$ or $s$-coordinate, the slope between the level and the geopotential surfaces must be added to 
     567$\mathbb{R}$ for the discret form to be exact.  
     568 
     569Such a choice of discretisation is consistent with the iso-neutral operator as 
     570it uses the same definition for the slopes. 
     571It also ensures the conservation of the tracer variance (see Appendix \autoref{apdx:eiv_skew}), 
     572$i.e.$ it does not include a diffusive component but is a "pure" advection term. 
    578573 
    579574 
    … …  
    591586This part will be moved in an Appendix. 
    592587 
    593 The continuous property to be demonstrated is : 
     588The continuous property to be demonstrated is: 
    594589\begin{align*} 
    595590\int_D  D_l^T \; T \;dv   \leq 0 
    … …  
    642637% 
    643638\allowdisplaybreaks 
    644 \intertext{The summation is done over all $i$ and $k$ indices, it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to regroup all the terms of the summation by triad at a ($i$,$k$) point. In other words, we regroup all the terms in the neighbourhood  that contain a triad at the same ($i$,$k$) indices. It becomes: } 
     639  \intertext{The summation is done over all $i$ and $k$ indices, 
     640  it is therefore possible to introduce a shift of $-1$ either in $i$ or $k$ direction in order to 
     641  regroup all the terms of the summation by triad at a ($i$,$k$) point. 
     642  In other words, we regroup all the terms in the neighbourhood that contain a triad at the same ($i$,$k$) indices. 
     643  It becomes: } 
    645644% 
    646645&\equiv -\sum_{i,k} 
    … …  
    672671% 
    673672\allowdisplaybreaks 
    674 \intertext{Then outing in factor the triad in each of the four terms of the summation and substituting the triads by their expression given in \autoref{eq:Gf_triads}. It becomes: } 
     673  \intertext{Then outing in factor the triad in each of the four terms of the summation and 
     674  substituting the triads by their expression given in \autoref{eq:Gf_triads}. 
     675  It becomes: } 
    675676% 
    676677&\equiv -\sum_{i,k} 
    … …  
    710711The last inequality is obviously obtained as we succeed in obtaining a negative summation of square quantities. 
    711712 
    712 Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$, then the previous demonstration would have let to: 
     713Note that, if instead of multiplying $D_l^T$ by $T$, we were using another tracer field, let say $S$, 
     714then the previous demonstration would have let to: 
    713715\begin{align*} 
    714716\int_D  S \; D_l^T  \;dv &\equiv  \sum_{i,k} \left\{ S \ D_l^T \ b_T \right\}    \\ 
    … …  
    729731&\equiv  \sum_{i,k} \left\{ D_l^S \ T \ b_T \right\}    
    730732\end{align*}  
    731 This means that the iso-neutral operator is self-adjoint. There is no need to develop a specific to obtain it. 
     733This means that the iso-neutral operator is self-adjoint. 
     734There is no need to develop a specific to obtain it. 
    732735 
    733736 
    … …  
    745748This have to be moved in an Appendix. 
    746749 
    747 The continuous property to be demonstrated is : 
     750The continuous property to be demonstrated is: 
    748751\begin{align*} 
    749752\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0 
    … …  
    797800\end{matrix}    
    798801\end{align*} 
    799 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the  
    800 same but of opposite signs, they cancel out.  
    801 Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.  
    802 The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the  
    803 same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$  
    804 they cancel out with the neighbouring grid points.  
    805 Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the  
    806 $i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the  
    807 tracer is preserved by the discretisation of the skew fluxes. 
     802The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the same but of opposite signs, 
     803they cancel out.  
     804Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$. 
     805The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the same but both of opposite signs and 
     806shifted by 1 in $k$ direction. 
     807When summing over $k$ they cancel out with the neighbouring grid points. 
     808Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the $i$ direction. 
     809Therefore the sum over the domain is zero, 
     810$i.e.$ the variance of the tracer is preserved by the discretisation of the skew fluxes. 
    808811 
    809812\end{document} 
Note: See TracChangeset for help on using the changeset viewer.