source: trunk/NEMO/DOC/BETA/Chapters/Annex_E.tex @ 707

% ================================================================
% Appendix E : UBS scheme
% ================================================================
\chapter{UBS scheme}
\label{AN_E}
\minitoc


% -------------------------------------------------------------------------------------------------------------
%        UBS scheme  
% -------------------------------------------------------------------------------------------------------------
\section{Upstream Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=T)}
\label{TRA_adv_ubs}

The UBS advection scheme is an upstream biased third order scheme based on 
an upstream-biased parabolic interpolation. It is also known as Cell 
Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective 
Kinematics). For example, in the $i$-direction :
\begin{equation} \label{Eq_tra_adv_ubs2}
\tau _u^{ubs} = \left\{  \begin{aligned}
  & \tau _u^{cen4} + \frac{1}{12} \,\tau"_i     & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
  & \tau _u^{cen4} - \frac{1}{12} \,\tau"_{i+1} & \quad \text{if }\ u_{i+1/2}       <       0
                   \end{aligned}    \right.
\end{equation}
or equivalently, the advective flux is
\begin{equation} \label{Eq_tra_adv_ubs2}
U_{i+1/2} \ \tau _u^{ubs} 
=U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2}
- \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
\end{equation}
where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. By choosing this expression for $\tau "$ we consider a fourth order approximation of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$). 

Alternative choice: introduce the scale factors:  $\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]$.


This results in a dissipatively dominant (i.e. hyper-diffusive) truncation 
error \citep{Sacha2005}. The overall performance of the 
advection scheme is similar to that reported in \cite{Farrow1995}. 
It is a relatively good compromise between accuracy and smoothness. It is 
not a \emph{positive} scheme meaning false extrema are permitted but the 
amplitude of such are significantly reduced over the centred second order 
method. Nevertheless it is not recommended to apply it to a passive tracer 
that requires positivity. 

The intrinsic diffusion of UBS makes its use risky in the vertical direction 
where the control of artificial diapycnal fluxes is of paramount importance. 
It has therefore been preferred to evaluate the vertical flux using the TVD 
scheme when \np{ln\_traadv\_ubs}=T.

For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds to a 
second order centred scheme is evaluated using the \textit{now} velocity (centred in 
time) while the second term which is the diffusive part of the scheme, is 
evaluated using the \textit{before} velocity (forward in time. This is discussed by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK 
schemes only differ by one coefficient. Substituting 1/6 with 1/8 in 
(\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}. This 
option is not available through a namelist parameter, since the 1/6 
coefficient is hard coded. Nevertheless it is quite easy to make the 
substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme

NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can 
be controlled by vertical advection (not vertical diffusion which is usually 
solved using an implicit scheme). Computer time can be saved by using a 
time-splitting technique on vertical advection. This possibility have been 
implemented and validated in ORCA05-L301. It is not currently offered in the 
current reference version. 

NB 2 : In a forthcoming release four options will be proposed for the 
vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be 
evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme , or  \textit{(b)} a TVD 
scheme, or  \textit{(c)} an interpolation based on conservative parabolic splines 
following \citet{Sacha2005} implementation of UBS in ROMS, 
or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an 
eight-order accurate conventional scheme.

NB 3 : It is straight forward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
\begin{equation} \label{Eq_tra_adv_ubs2}
\tau _u^{ubs} = \left\{  \begin{aligned}
   & \tau _u^{cen4} + \frac{1}{12} \tau"_i      & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
   & \tau _u^{cen4} - \frac{1}{12} \tau"_{i+1}  & \quad \text{if }\ u_{i+1/2}       <       0
                   \end{aligned}    \right.
\end{equation}
or equivalently 
\begin{equation} \label{Eq_tra_adv_ubs2}
\begin{split}
e_{2u} e_{3u}\,u_{i+1/2} \ \tau _u^{ubs} 
&= e_{2u} e_{3u}\,u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\tau"_i }^{\,i+1/2} \\
& - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
\end{split}
\end{equation}
\eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that the UBS scheme is based on the fourth order scheme to which is added an upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order part as stated above using \eqref{Eq_tra_adv_ubs}. Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient with is simply proportional to the velocity.


laplacian diffusion:
\begin{equation} \label{Eq_tra_ldf_lap}
\begin{split}
D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta _i 
\left[ {A_u^{lT} \frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} 
\left[ T \right]} \right]} \right.
\\
&\ \left. {+\; \delta _j \left[ 
{A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta _{j+1/2} \left[ T 
\right]} \right)} \right]\quad } \right]
\end{split}
\end{equation}

bilaplacian:
\begin{equation} \label{Eq_tra_ldf_lap}
\begin{split}
D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
& \delta _i \left[  \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta _{i+1/2} 
        \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}}
    \delta _i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta _{i+1/2} 
           [T] \right] \right] \right]
\end{split}
\end{equation}
with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, 
$i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$
it comes :
\begin{equation} \label{Eq_tra_ldf_lap}
\begin{split}
D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\;  e_{3T}} \\
& \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta _{i+1/2} 
       \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} 
    \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta _{i+1/2} 
         [T] \right] \right] \right]
\end{split}
\end{equation}
if the velocity is uniform ($i.e.$ $|u|=cst$) then the diffusive flux is
\begin{equation} \label{Eq_tra_ldf_lap}
\begin{split}
F_u^{lT} = - \frac{1}{12}
 e_{2u}\,e_{3u}\,|u| \;\sqrt{ e_{1u}}\,\delta _{i+1/2} 
       \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} 
    \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}}\:\delta _{i+1/2} 
         [T] \right] \right]
\end{split}
\end{equation}
beurk....  reverte the logic: starting from the diffusive part of the advective flux it comes:

\begin{equation} \label{Eq_tra_adv_ubs2}
\begin{split}
F_u^{lT}
&= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
\end{split}
\end{equation}
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]$

sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$):
\begin{equation} \label{Eq_tra_adv_ubs2}
\begin{split}
F_u^{lT}
&= - \frac{1}{12} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{e_{1T}^3\,|u|}{e_{1T}e_{2T}\,e_{3T}}\,\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right] \right]
\end{split}
\end{equation}
which leads to ${A_T^{lT}}^2 = \frac{1}{12} {e_{1T}}^3\ \overline{|u|}^{\,i+1/2}$

sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$
\begin{equation} \label{Eq_tra_adv_ubs2}
\begin{split}
F_u^{lT}
&= - \frac{1}{12} {e_{1u}}^1 \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{2T}\,e_{3T}}\,\delta _i \left[ \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right] \right] \\
&= - \frac{1}{12} e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{1T}\,e_{2T}\,e_{3T}}\,\delta _i \left[ e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u} }{e_{1u}} \delta _{i+1/2}[\tau] \right] \right]
\end{split}
\end{equation}
which leads to ${A_u^{lT}} = \frac{1}{12} {e_{1u}}^3\ |u|$


% -------------------------------------------------------------------------------------------------------------
%        Leap-Frog energetic  
% -------------------------------------------------------------------------------------------------------------
\section{Leap-Frog energetic }
\label{LF}

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:
\begin{subequations} \label{dt_mt}
\begin{align}
 \delta _{t+\Delta t/2} [q]     &=  \  \ \,   q^{t+\Delta t}  - q^{t}      \\
 \overline q^{\,t+\Delta t/2} &= \left\{ q^{t+\Delta t} + q^{t} \right\} \; / \; 2
\end{align}
\end{subequations}
As for space operator, the adjoint of the derivation and averaging time operators are $\delta_t^*=\delta_{t+\Delta t/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$, respectively. 

The Leap-frog time stepping given by \eqref{Eq_DOM_nxt} can be defined as:
\begin{equation} \label{LF}
   \frac{\partial q}{\partial t} 
         \equiv \frac{1}{\Delta t} \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} 
      =         \frac{q^{t+\Delta t}-q^{t-\Delta t}}{2\Delta t}
\end{equation} 
Note that \eqref{LF} shows that the leapfrog time step is $\Delta t$, not $2\Delta t$ as it can be found sometime in literature. 
The leap-Frog time stepping is a second order centered scheme. As such it respects the quadratic invariant in integral forms, $i.e.$ the following continuous property,
\begin{equation} \label{Energy}
\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} 
   =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt} 
   =  \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) ,
\end{equation}
is satisfied in discrete form. Indeed, 
\begin{equation} \begin{split}
\int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} 
   &\equiv \sum\limits_{0}^{N} 
         {\frac{1}{\Delta t} q^t \ \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} \ \Delta t} 
      \equiv \sum\limits_{0}^{N}  { q^t \ \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} } \\
   &\equiv \sum\limits_{0}^{N}  { \overline{q}^{\,t+\Delta/2}{ \delta _{t+\Delta t/2}[q]}}
      \equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta _{t+\Delta t/2}[q^2] }\\
   &\equiv \sum\limits_{0}^{N}  { \frac{1}{2} \delta _{t+\Delta t/2}[q^2] }
      \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right)
\end{split} \end{equation}
NB here pb of boundary condition when applying the adjoin! In space, setting to 0 the quantity in land area is sufficient to get rid of the boundary condition (equivalently of the boundary value of the integration by part). In time this boundary condition is not physical and \textbf{add something here!!!}