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

$$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

$$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:

$$\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.$$ (74)

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

$$\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

$$\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

$$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:

$$\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),

$$\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),

$$\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