adv_split_u - zonal advection of 2D quantities

INTERFACE:

    subroutine adv_split_u(dt,f,fi,Di,adv,U,DU,   &
 #if defined(SPHERICAL) || defined(CURVILINEAR)
                           dxu,dyu,arcd1,         &
 #endif
                           splitfac,scheme,AH,    &
                           mask_flux,mask_update)
    Note (KK): Keep in sync with interface in advection.F90
DESCRIPTION:

Executes an advection step in zonal direction for a 2D quantity. The 1D advection equation

$\displaystyle D^n_{i,j} c^n_{i,j} =
D^o_{i,j} c^o_{i,j}
\displaystyle
- \Delta ...
...i-1,j}\tilde c^u_{i-1,j}\Delta y^u_{i-1,j}
}{\Delta x^c_{i,j}\Delta y^c_{i,j}},$ (72)

is accompanied by an fractional step for the 1D continuity equation

$\displaystyle D^n_{i,j} =
D^o_{i,j}
\displaystyle
- \Delta t
\frac{
U_{i,j}\Delta y^u_{i,j}-
U_{i-1,j}\Delta y^u_{i-1,j}
}{\Delta x^c_{i,j}\Delta y^c_{i,j}}.$ (73)

Here, $n$ and $o$ denote values before and after this operation, respectively, $n$ denote intermediate values when other 1D advection steps come after this and $o$ denotes intermediate values when other 1D advection steps came before this. Furthermore, when this $u$-directional split step is repeated during the total time step (Strang splitting), the time step $\Delta t$ denotes a fraction of the full time step.

The interfacial fluxes $\tilde c^u_{i,j}$ are calculated according to the third-order polynomial scheme (so-called P$_2$ scheme), cast in Lax-Wendroff form by:

\begin{displaymath}\tilde c_{i,j}=
\left\{
\begin{array}{ll}
\left(c_{i,j}+\frac...
...
(c_{i,j}-c_{i+1,j})\right) & \mbox{ else, }
\end{array}\right.\end{displaymath} (74)

with the Courant number $C_{i,j}=u_{i,j}\Delta t/\Delta x$ and

$\displaystyle \tilde c_{i,j}^+=\alpha_{i,j}+\beta_{i,j}r^+_{i,j}, \quad
\tilde c_{i,j}^-=\alpha_{i,j}+\beta_{i,j}r^-_{i,j},$ (75)

where

$\displaystyle \alpha_{i,j}=\frac12 +\frac16(1-2\vert C_{i,j}\vert),\quad
\beta _{i,j}=\frac12 -\frac16(1-2\vert C_{i,j}\vert),$ (76)

and

$\displaystyle r^+_{i,j}=\frac{c_{i,j}-c_{i-1,j}}{c_{i+1,j}-c_{i,j}},
\quad
r^-_{i,j}=\frac{c_{i+2,j}-c_{i+1,j}}{c_{i+1,j}-c_{i,j}}.$ (77)

It should be noted, that for $\tilde c_{i,j}^+=\tilde c_{i,j}^-=1$ the original Lax-Wendroff scheme and for $\tilde c_{i,j}^+=\tilde c_{i,j}^-=0$ the first-oder upstream scheme can be recovered.

In order to obtain a monotonic and positive scheme, the factors $\tilde c_{i,j}^+$ are limited in the following way:

$\displaystyle \tilde c_{i,j}^+ \rightarrow \max \left[
0,\min\left(\tilde c_{i,...
...2}{1-\vert C_{i,j}\vert},
\frac{2r^+_{i,j}}{\vert C_{i,j}\vert}\right)
\right],$ (78)

and, equivalently, for $\tilde c_{i,j}^-$. This so-called PDM-limiter has been described in detail by Leonard (1991), who named the PDM-limited P$_2$ scheme also ULTIMATE QUICKEST (quadratic upstream interpolation for convective kinematics with estimated stream terms).

Some simpler limiters which do not exploit the third-order polynomial properties of the discretisation (74) have been listed by Zalezak (1987). Among those are the MUSCL scheme by van Leer (1979),

$\displaystyle \tilde c_{i,j}^+ \rightarrow \max \left[
0,\min\left(
2,2r^+_{i,j},\frac{1+r^+_{i,j}}{2}
\right)
\right],$ (79)

and the Superbee scheme by Roe (1985),

$\displaystyle \tilde c_{i,j}^+ \rightarrow \max \left[
0,\min(1,2r^+_{i,j}),\min(r^+_{i,j},2)
\right].$ (80)

The selector for the schemes is scheme:



scheme = UPSTREAM: first-order upstream (monotone)
scheme = P2: third-order polynomial (non-monotone)
scheme = SUPERBEE: second-order TVD (monotone)
scheme = MUSCL: second-order TVD (monotone)
scheme = P2_PDM: third-order ULTIMATE-QUICKEST (monotone)



Furthermore, the horizontal diffusion in zonal direction with the constant diffusion coefficient AH is carried out here by means of a central difference second-order scheme. USES:

    use domain, only: imin,imax,jmin,jmax
 #if !( defined(SPHERICAL) || defined(CURVILINEAR) )
    use domain, only: dx,dy,ard1
 #endif
    use advection, only: adv_interfacial_reconstruction
    use advection, only: UPSTREAM
  $ use omp_lib
    IMPLICIT NONE
INPUT PARAMETERS:
    Note (KK): in general dxu, dyu and mask_flux do only have valid data
               within (_IRANGE_HALO_-1,_JRANGE_HALO_). In some cases the
               original field extension may even be _IRANGE_HALO_. Then
               explicit declared array bounds _IRANGE_HALO_-1 require a
               provision of the corresponding subarray and will cause
               copying of the non-contiguously data into a temporarily
               array. Therefore they are declared as pointers here. This
               however requires, that the provided pointers already carry
               the correct bounds.
    REALTYPE,intent(in)                        :: dt,splitfac,AH
    REALTYPE,dimension(E2DFIELD),intent(in)    :: f,U,DU
 #if defined(SPHERICAL) || defined(CURVILINEAR)
    REALTYPE,dimension(:,:),pointer,contiguous,intent(in) :: dxu,dyu
    REALTYPE,dimension(E2DFIELD),intent(in)    :: arcd1
 #endif
    integer,intent(in)                         :: scheme
    logical,dimension(:,:),pointer,contiguous,intent(in)  :: mask_flux
    logical,dimension(E2DFIELD),intent(in)     :: mask_update
INPUT/OUTPUT PARAMETERS:
    REALTYPE,dimension(E2DFIELD),intent(inout) :: fi,Di,adv
LOCAL VARIABLES:
    REALTYPE,dimension(E2DFIELD) :: uflux
    logical            :: use_limiter,use_AH
    integer            :: i,j,isub
    REALTYPE           :: dti,Dio,advn,cfl,fuu,fu,fd
REVISION HISTORY:
    Original author(s): Hans Burchard & Karsten Bolding