adv_split_u - zonal advection of 2D quantities


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

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:

$\displaystyle \tilde c_{i,j}= \left\{ \begin{array}{ll} \left(c_{i,j}+\frac12 \...
...t C_{i,j}\vert) (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

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


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


$\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
    use advection, only: adv_interfacial_reconstruction
    use advection, only: UPSTREAM
  $ use omp_lib
    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,intent(in) :: dxu,dyu
    REALTYPE,dimension(E2DFIELD),intent(in)    :: arcd1
    integer,intent(in)                         :: scheme
    logical,dimension(:,:),pointer,intent(in)  :: mask_flux
    logical,dimension(E2DFIELD),intent(in)     :: mask_update
    REALTYPE,dimension(E2DFIELD),intent(inout) :: fi,Di,adv
    REALTYPE,dimension(E2DFIELD) :: uflux
    logical            :: use_limiter,use_AH
    integer            :: i,j,isub
    REALTYPE           :: dti,Dio,advn,cfl,fuu,fu,fd
    Original author(s): Hans Burchard & Karsten Bolding

kklingbe 2017-10-02