!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE ADVPTFCT                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE advptfct(nx,ny,nz,dtbig1,ptprt,u,v,w,wcont,                  & 1,12
           rhostr,rhostri,ptbar,mapfct,j3,j3inv,                        &
           ptadv,                                                       &
           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                     &
           ustr,vstr,wstr,tem1_0,tem2_0,tem3_0,mp_tem)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the advection of total potential temperature
!  using flux-correctd transport scheme.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue and Yifeng Tang
!  1/6/93.
!
!  10/12/1996 (M. Xue)
!  Code clean up. prprt*rhostr is now passed into advctsfct.
!
!  7/10/1997 (Fanyou Kong -- CMRP)
!  Include MPDCD simple positive definite advection scheme
!  (sadvopt = 5) through subroutine 'advmpdcd'
!
!  11/20/1997 (Fanyou Kong -- CMRP)
!  Added map factor into both FCT and MPDCD schemes
!
!  7/17/1998 (Ming Xue)
!  Significant modifications to call a new version of advctsfct.
!
!  10/5/1998 (Dan Weber)
!  Added rhostri, mapfct 7-8 array's.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity (m/s)
!    v        y component of velocity (m/s)
!    w        Vertical component of Cartesian velocity (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!    ptprt    Perturbation potential temperature at all time levels (K)
!
!    rhostr   Base state air density times j3 (kg/m**3)
!    rhostri  Inverse base state air density times j3 (kg/m**3)
!    ptbar    Base state potential temperature (K)
!    mapfct   Map factor for scalar point
!
!  OUTPUT:
!
!    ptadv    Advection term of potential temperature eqn (kg/m**3)*K/s
!
!  WORK ARRAYS:
!
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!    tem6     Temporary work array.
!    tem7     Temporary work array.
!    tem8     Temporary work array.
!    ustr     Work array
!    vstr     Work array
!    wstr     Work array
!
!    tem1_0   Temporary work array.
!    tem2_0   Temporary work array.
!    tem3_0   Temporary work array.
!
!    mp_tem   Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'timelvls.inc'

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  REAL :: dtbig1

  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: u     (nx,ny,nz,nt)  ! u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! w-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Contravariant w-velocity (m/s)
!
  REAL :: rhostr(nx,ny,nz)     ! Base state air density times j3 (kg/m**3)
  REAL :: rhostri(nx,ny,nz)    ! Inv. base state air density times j3 (kg/m**3)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: mapfct(nx,ny,8)      ! Map factor for scalar point
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ).
  REAL :: j3inv (nx,ny,nz)     ! 1/J3.
!
  REAL :: ptadv(nx,ny,nz)      ! Advection of potential temperature
!
!
  REAL :: tem1  (nx,ny,nz)     ! temporary work array
  REAL :: tem2  (nx,ny,nz)     ! temporary work array
  REAL :: tem3  (nx,ny,nz)     ! temporary work array
  REAL :: tem4  (nx,ny,nz)     ! temporary work array
  REAL :: tem5  (nx,ny,nz)     ! temporary work array
  REAL :: tem6  (nx,ny,nz)     ! temporary work array
  REAL :: tem7  (nx,ny,nz)     ! temporary work array
  REAL :: tem8  (nx,ny,nz)     ! temporary work array
  REAL :: ustr  (nx,ny,nz)     ! Work array holding u*rhostr/mapfct
  REAL :: vstr  (nx,ny,nz)     ! Work array holding v*rhostr/mapfct
  REAL :: wstr  (nx,ny,nz)     ! Work array holding wcont*rhostr

  REAL :: tem1_0(0:nx,0:ny,0:nz) ! Work array
  REAL :: tem2_0(0:nx,0:ny,0:nz) ! Work array
  REAL :: tem3_0(0:nx,0:ny,0:nz) ! Work array

  REAL :: mp_tem(MAX(nx+1,ny+1)*(nz+1))  ! Temporary message passing array.

!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  REAL :: ptmin, ptmin1, ptmax1, ptmin2, ptmax2, tema

  INTEGER :: advdiv     ! = 0, neglect anelastic correction in FCT advection
                        ! = 1, include anelastic correction in FCT advection
  INTEGER :: mptag1,mptag2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'       ! Global model control parameters
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'bndry.inc'          ! Boundary Condition Flags
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  advdiv    = 1  ! Include the anelastic correction term? (recommended!)

  IF( fctadvptprt == 1 .AND. sadvopt == 5 ) fctadvptprt = 2

  IF( ptsmlstp == 1 ) fctadvptprt = 1  ! Do not reset this one

  IF( fctadvptprt == 0 ) THEN

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1_0(i,j,k)=ptprt(i,j,k,1)+ptbar(i,j,k)
          tem2_0(i,j,k)=ptprt(i,j,k,2)+ptbar(i,j,k)
        END DO
      END DO
    END DO

  ELSE IF( fctadvptprt == 1 ) THEN

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1_0(i,j,k)=ptprt(i,j,k,1)
          tem2_0(i,j,k)=ptprt(i,j,k,2)
        END DO
      END DO
    END DO

  ELSE IF( fctadvptprt == 2 ) THEN

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1(i,j,k)=ptprt(i,j,k,1)+ptbar(i,j,k)
          tem2(i,j,k)=ptprt(i,j,k,2)+ptbar(i,j,k)
        END DO
      END DO
    END DO

    CALL a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,                          &
                 1,nz,1,nz-1,ptmax1,ptmin1)
    CALL a3dmax0(tem2,1,nx,1,nx-1,1,ny,1,ny-1,                          &
                 1,nz,1,nz-1,ptmax2,ptmin2)
    ptmin = AMIN1(ptmin1,ptmin2)

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1_0(i,j,k)=tem1(i,j,k)-ptmin
          tem2_0(i,j,k)=tem2(i,j,k)-ptmin
        END DO
      END DO
    END DO

  END IF

!  IF (mp_opt > 0) THEN
!    CALL acct_interrupt(mp_acct)
!    CALL mpsendextew(tem1_0,nx,ny,nz,ebc,wbc,mptag1,mp_tem)
!    CALL mpsendextew(tem2_0,nx,ny,nz,ebc,wbc,mptag2,mp_tem)
!    CALL mprecvextew(tem1_0,nx,ny,nz,ebc,wbc,mptag1,mp_tem)
!    CALL mprecvextew(tem2_0,nx,ny,nz,ebc,wbc,mptag2,mp_tem)
!    CALL mpsendextns(tem1_0,nx,ny,nz,nbc,sbc,mptag1,mp_tem)
!    CALL mpsendextns(tem2_0,nx,ny,nz,nbc,sbc,mptag2,mp_tem)
!    CALL mprecvextns(tem1_0,nx,ny,nz,nbc,sbc,mptag1,mp_tem)
!    CALL mprecvextns(tem2_0,nx,ny,nz,nbc,sbc,mptag2,mp_tem)
!  END IF
  CALL acct_interrupt(bc_acct)
  CALL extndsbc(tem1_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL extndsbc(tem2_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL acct_stop_inter
!test tmp: call test_dump (tem1_0,"XXXAtem1_0",nx+1,ny+1,nz+1,4,1)
!test tmp: call test_dump (tem2_0,"XXXAtem2_0",nx+1,ny+1,nz+1,4,1)
!
!-----------------------------------------------------------------------
!
!  Calculate ustr, vstr and wstr.
!
!-----------------------------------------------------------------------
!

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem3_0(i,j,k)=rhostr(i,j,k)
      END DO
    END DO
  END DO

!  IF (mp_opt > 0) THEN
!    CALL acct_interrupt(mp_acct)
!    CALL mpsendextew(tem3_0,nx,ny,nz,ebc,wbc,mptag1,mp_tem)
!    CALL mprecvextew(tem3_0,nx,ny,nz,ebc,wbc,mptag1,mp_tem)
!    CALL mpsendextns(tem3_0,nx,ny,nz,nbc,sbc,mptag1,mp_tem)
!    CALL mprecvextns(tem3_0,nx,ny,nz,nbc,sbc,mptag1,mp_tem)
!  END IF
  CALL acct_interrupt(bc_acct)
  CALL extndsbc(tem3_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL acct_stop_inter
!test tmp: call test_dump (tem3_0,"XXXAtem3_0",nx+1,ny+1,nz+1,4,1)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        ustr(i,j,k)=u(i,j,k,2)*(tem3_0(i-1,j,k)+tem3_0(i,j,k))          &
                   *mapfct(i,j,5)*0.5
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        vstr(i,j,k)=v(i,j,k,2)*(tem3_0(i,j-1,k)+tem3_0(i,j,k))          &
                   *mapfct(i,j,6)*0.5
      END DO
    END DO
  END DO

  DO k=1,nz
    DO j=1,ny-1
      DO i=1,nx-1
        wstr(i,j,k)=wcont(i,j,k)                                        &
                   *(tem3_0(i,j,k-1)+tem3_0(i,j,k))*0.5
      END DO
    END DO
  END DO

  IF(advdiv == 1) THEN
!
!-----------------------------------------------------------------------
!
!  Add the anelastic divergence correction term:
! -[m**2{d(rhostr*u/m)/dx+ d(rhostr*v/m)/dy}+ d(rhostr*wcont)/dz]*ptprt.
!
!-----------------------------------------------------------------------
!
    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          ptadv(i,j,k)= - tem2_0(i,j,k)*(                               &
               ((ustr(i+1,j,k)-ustr(i,j,k))*dxinv                       &
               +(vstr(i,j+1,k)-vstr(i,j,k))*dyinv)*mapfct(i,j,7)        &
               +(wstr(i,j,k+1)-wstr(i,j,k))*dzinv )
        END DO
      END DO
    END DO

  ELSE

    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          ptadv(i,j,k)= 0.0
        END DO
      END DO
    END DO

  END IF

  IF (sadvopt == 4) THEN

    IF(fctadvptprt == 1.AND.ptsmlstp /= 1 ) THEN  ! Advection of ptbar

      tema = 0.25/dz

      DO k=2,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem2(i,j,k)=(rhostr(i,j,k  )*j3inv(i,j,k  )                 &
                        +rhostr(i,j,k-1)*j3inv(i,j,k-1))*w(i,j,k,2)     &
                       *(ptbar(i,j,k)-ptbar(i,j,k-1))*tema
          END DO
        END DO
      END DO

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            ptadv(i,j,k)=ptadv(i,j,k)+tem2(i,j,k)+tem2(i,j,k+1)
          END DO
        END DO
      END DO

    END IF

    CALL advctsfct(nx,ny,nz,nt,dtbig1,ustr,vstr,wstr,    &
         mapfct,tem1_0,tem2_0,tem3_0,rhostr,rhostri,        &
         ptadv,                                                         &
         tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

  ELSE IF(sadvopt == 5) THEN

    IF(fctadvptprt /= 0 ) THEN
      CALL arpsstop(                                                    &
          'Subourinte ADVMPDCD should never be used to advect '//       &
          'non-positive definite variable PTPRT.'//                     &
          'Job Stopped in subroutine ADVPTFCT.',1)
    END IF

    CALL advmpdcd(nx,ny,nz,nt,dtbig1,                          &
         ustr,vstr,wstr,mapfct, tem1_0,tem2_0, rhostr,      &
         ptadv,                                                         &
         tem1,tem2,tem3,tem4, tem3_0, mp_tem)

  END IF

  RETURN
END SUBROUTINE advptfct

!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE ADVQFCT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE advqfct(nx,ny,nz,dtbig1,q,u,v,wcont, ustr,vstr,wstr,         & 3,6
           rhostr,rhostri,mapfct,j3,qadv,                               &
           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                     &
           tem1_0,tem2_0,tem3_0)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the advection of water substance
!  using flux-correctd transport scheme.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yifeng Tang and M. Xue
!  1/6/93.
!
!  10/12/1996 (M. Xue)
!  Code clean up. q*rhostr is now passed into advctsfct.
!
!  7/10/1997 (Fanyou Kong -- CMRP)
!  Include MPDCD simple positive definite advection scheme
!  (sadvopt = 5) through subroutine 'advmpdcd'
!
!  11/20/1997 (Fanyou Kong -- CMRP)
!  Added map factor into FCT and MPDCD schemes
!
!  7/17/1998 (Ming Xue)
!  Significant modifications to call a new version of advctsfct.
!
!  9/28/1998 (Dan Weber)
!  Moved calculation of u,v,wstr out of this subroutine and added
!  rhostri and mapfct 7-8.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity (m/s)
!    v        y component of velocity (m/s)
!    w        Vertical component of Cartesian velocity (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!    ustr     Appropriate ustr for advection
!    vstr     Appropriate ustr for advection
!    wstr     Appropriate ustr for advection
!    q        Water substance at all time levels (K)
!
!    rhostr   Base state air density times j3 (kg/m**3)
!    rhostri  Inverse base state air density times j3 (kg/m**3)
!    mapfct   Map factor for scalar point
!
!  OUTPUT:
!
!    qadv    Advection term of water substance eqn (kg/m**3)*(kg/kg)/s
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!    tem6     Temporary work array.
!    tem7     Temporary work array.
!    tem8     Temporary work array.
!    tem1_0   Temporary work array.
!    tem2_0   Temporary work array.
!    tem3_0   Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'timelvls.inc'
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  REAL :: dtbig1

  REAL :: q     (nx,ny,nz,nt)  ! Water substance  (kg/kg)
  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Total w-velocity (m/s)
  REAL :: ustr  (nx,ny,nz)     ! u*rhostr/mapfct
  REAL :: vstr  (nx,ny,nz)     ! v*rhostr/mapfct
  REAL :: wstr  (nx,ny,nz)     ! wcont*rhostr
!
  REAL :: rhostr(nx,ny,nz)     ! Base state air density times j3 (kg/m**3)
  REAL :: rhostri(nx,ny,nz)    ! Inv. base state air density times j3 (kg/m**3)
  REAL :: mapfct(nx,ny,8)      ! Map factor for scalar point
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ).
!
  REAL :: qadv  (nx,ny,nz)     ! Advection of water substance
!
  REAL :: tem1  (nx,ny,nz)     ! temporary work array
  REAL :: tem2  (nx,ny,nz)     ! temporary work array
  REAL :: tem3  (nx,ny,nz)     ! temporary work array
  REAL :: tem4  (nx,ny,nz)     ! temporary work array
  REAL :: tem5  (nx,ny,nz)     ! temporary work array
  REAL :: tem6  (nx,ny,nz)     ! temporary work array
  REAL :: tem7  (nx,ny,nz)     ! temporary work array
  REAL :: tem8  (nx,ny,nz)     ! temporary work array
  REAL :: tem1_0(0:nx,0:ny,0:nz) ! automatic work array
  REAL :: tem2_0(0:nx,0:ny,0:nz) ! automatic work array
  REAL :: tem3_0(0:nx,0:ny,0:nz) ! automatic work array

!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  INTEGER :: advdiv

  INTEGER :: mptag1,mptag2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'       ! Global model control parameters
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'bndry.inc'         ! Boundary Condition Flags
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  advdiv = 1  ! Include anelastic correction term in advection
!
!-----------------------------------------------------------------------
!
!  Calculate the q(tfurture) with FCT correction.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1_0(i,j,k)=q(i,j,k,1)
        tem2_0(i,j,k)=q(i,j,k,2)
      END DO
    END DO
  END DO

!  IF (mp_opt > 0) THEN
!    CALL acct_interrupt(mp_acct)
!    CALL mpsendextew(tem1_0,nx,ny,nz,ebc,wbc,mptag1,tem1)
!    CALL mpsendextew(tem2_0,nx,ny,nz,ebc,wbc,mptag2,tem1)
!    CALL mprecvextew(tem1_0,nx,ny,nz,ebc,wbc,mptag1,tem1)
!    CALL mprecvextew(tem2_0,nx,ny,nz,ebc,wbc,mptag2,tem1)
!    CALL mpsendextns(tem1_0,nx,ny,nz,nbc,sbc,mptag1,tem1)
!    CALL mpsendextns(tem2_0,nx,ny,nz,nbc,sbc,mptag2,tem1)
!    CALL mprecvextns(tem1_0,nx,ny,nz,nbc,sbc,mptag1,tem1)
!    CALL mprecvextns(tem2_0,nx,ny,nz,nbc,sbc,mptag2,tem1)
!  END IF
  CALL acct_interrupt(bc_acct)
  CALL extndsbc(tem1_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL extndsbc(tem2_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL acct_stop_inter
!test tmp: call test_dump (tem1_0,"XXXA1tem1_0",nx+1,ny+1,nz+1,4,1)
!test tmp: call test_dump (tem2_0,"XXXA1tem2_0",nx+1,ny+1,nz+1,4,1)
!
!-----------------------------------------------------------------------
!
!  Calculate ustr, vstr and wstr.
!
!-----------------------------------------------------------------------
!

  IF( advdiv == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  Add the anelastic divergence correction term:
!  -[d(rhostr*u)/dx+ d(rhostr*v)/dy+ d(rhostr*wcont)/dz]*q
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qadv(i,j,k)= - tem2_0(i,j,k)*(                                &
              ((ustr(i+1,j,k)-ustr(i,j,k))*dxinv                        &
              +(vstr(i,j+1,k)-vstr(i,j,k))*dyinv)*mapfct(i,j,7)         &
              +(wstr(i,j,k+1)-wstr(i,j,k))*dzinv)
        END DO
      END DO
    END DO
!test tmp: call test_dump (qadv,"XXXqadv",nx,ny,nz,0,0)
!test tmp: call test_dump (mapfct(1,1,7),"XXXqadv_mapfct",nx,ny,1,0,0)
!test tmp: call test_dump (tem2_0,"XXXqadv_tem2_0",nx+1,ny+1,nz+1,0,0)
!test tmp: call test_dump (ustr,"XXXqadv_ustr",nx,ny,nz,1,0)
!test tmp: call test_dump (vstr,"XXXqadv_vstr",nx,ny,nz,2,0)
!test tmp: call test_dump (wstr,"XXXqadv_wstr",nx,ny,nz,3,0)

  ELSE

    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          qadv(i,j,k)= 0.0
        END DO
      END DO
    END DO

  END IF

  IF(sadvopt == 4) THEN

    CALL advctsfct(nx,ny,nz,nt,dtbig1,ustr,vstr,wstr,    &
         mapfct(1,1,7),tem1_0,tem2_0,tem3_0,rhostr,rhostri,             &
         qadv,                                                          &
         tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

  ELSE IF(sadvopt == 5) THEN  ! MPDCD positive definite advection

    CALL advmpdcd(nx,ny,nz,nt,dtbig1,ustr,vstr,wstr,           &
         mapfct(1,1,7),tem1_0,tem2_0,rhostr,                &
         qadv,                                                          &
         tem1,tem2,tem3,tem4,tem3_0,tem5)

  END IF
!test tmp: call test_dump (qadv,"XXX2qadv",nx,ny,nz,0,0)

  RETURN
END SUBROUTINE advqfct
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE ADVCTSFCT                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE advctsfct(nx,ny,nz,nt,dtbig1,                 & 2,10
           ustr,vstr,wstr,mapfct2, s1,s2, s3,rhostr,rhostri,            &
           sadv,                                                        &
           fluxx,fluxy,fluxz,cx,cy,cz,tem1,tem2)
!
!----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Integrate scalar field forward on time step using Zalesak's
!  multi-dimensional version of FCT.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue. Adaptied for ARPS by Y.F. Tang.
!  11/29/1992.
!
!  3/10/1993 (M. Xue)
!  Code cleanup.
!
!  11/20/1997 (Fanyou Kong -- CMRP)
!  Added map factor into FCT schemes
!
!  7/17/1998 (Ming Xue)
!
!  The formulation of advective fluxes are now consistent with
!  the regular advection.  Array sadv is now accumulative.
!
!  10/5/1998 (Dan Weber)
!  Added rhostri array.
!
!  2000/09/15 (Gene Bassett)
!  Removed dx,dy,dz,dtbig & fctorderopt from argument list.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!    nt       number of the time level
!    dtbig1   Time step size
!
!    ustr     u*rhostr
!    vstr     v*rhostr
!    wstr     wcont*rhostr
!    mapfct2  map factor for the scalar ** 2
!    s1       scalar field at a given time n-1 level (scalar units)
!    s2       scalar field at a given time n level (scalar units)
!    rhostr   rhobar*j3
!    rhostri  inverse rhobar*j3
!    sadv     Part of the advection term for scalar s
!
!  OUTPUT:
!
!    sadv     Advection for scalar s
!    s1,s2 and s3 modified.
!
!  WORK ARRAYS:
!
!    fluxx     Temporary work array.
!    fluxy     Temporary work array.
!    fluxz     Temporary work array.
!    cx        Temporary work array.
!    cy        Temporary work array.
!    cz        Temporary work array.
!    tem1      Temporary work array.
!    tem2      Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
  INTEGER :: nt              ! The no. of t-levels of t-dependent arrays.
  REAL :: dtbig1

  REAL :: ustr  (nx,ny,nz)   ! u*rhostr
  REAL :: vstr  (nx,ny,nz)   ! v*rhostr
  REAL :: wstr  (nx,ny,nz)   ! wcont*rhostr
  REAL :: mapfct2(nx,ny)     ! map factor for scalar **2
  REAL :: s1    (0:nx,0:ny,0:nz)! Time lev. 1 of a scalar field to be advected.
  REAL :: s2    (0:nx,0:ny,0:nz)! Time lev. 2 of a scalar field to be advected.
  REAL :: s3    (0:nx,0:ny,0:nz)! Time lev. 3 of a scalar field to be advected.
  REAL :: rhostr(nx,ny,nz)   ! rhostr*j3
  REAL :: rhostri(nx,ny,nz)  ! Inverse rhostr*j3
  REAL :: sadv  (nx,ny,nz)   ! Advection of scalar s
!
  REAL :: fluxx (nx,ny,nz)   ! temporary work array
  REAL :: fluxy (nx,ny,nz)   ! temporary work array
  REAL :: fluxz (nx,ny,nz)   ! temporary work array
  REAL :: cx    (nx,ny,nz)   ! temporary work array
  REAL :: cy    (nx,ny,nz)   ! temporary work array
  REAL :: cz    (nx,ny,nz)   ! temporary work array
  REAL :: tem1  (nx,ny,nz)   ! temporary work array
  REAL :: tem2  (nx,ny,nz)   ! temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  REAL :: tem, rdxyz,dxyz, tema,rdt,tem43,tem16

  INTEGER :: mptag1,mptag2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'       ! Global model control parameters
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'bndry.inc'         ! Boundary Condition Flags
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  dxyz=dx*dy*dz
  rdxyz=1./(dxyz)

  tem43 = 4.0/3.0
  tem16 = 1.0/6.0
!
!-----------------------------------------------------------------------
!
!  Compute the second-order advective fluxes
!
!-----------------------------------------------------------------------
!
  CALL advflxs(nx,ny,nz,nt,dx,dy,dz,dtbig1,                             &
               ustr,vstr,wstr,s2,fctorderopt,                           &
               fluxx,fluxy,fluxz,                                       &
               s3,tem1)  ! Attention: s3 is used as a work array here!
!
!-----------------------------------------------------------------------
!
!  Compute the temperary scalar field at tfuture
!
!-----------------------------------------------------------------------
!
  tem = 2.0 * rdxyz
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        s3(i,j,k)=s1(i,j,k)-((fluxx(i+1,j,k)-fluxx(i,j,k)+              &
                              fluxy(i,j+1,k)-fluxy(i,j,k))*mapfct2(i,j)+ &
                              fluxz(i,j,k+1)-fluxz(i,j,k))*tem          &
                  *rhostri(i,j,k)
      END DO
    END DO
  END DO

!  IF (mp_opt > 0) THEN
!    CALL acct_interrupt(mp_acct)
!    CALL mpsendextew(s3,nx,ny,nz,ebc,wbc,mptag1,tem1)
!    CALL mprecvextew(s3,nx,ny,nz,ebc,wbc,mptag1,tem1)
!    CALL mpsendextns(s3,nx,ny,nz,nbc,sbc,mptag1,tem1)
!    CALL mprecvextns(s3,nx,ny,nz,nbc,sbc,mptag1,tem1)
!  END IF
  CALL acct_interrupt(bc_acct)
  CALL extndsbc(s3,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL acct_stop_inter
!test tmp: call test_dump (s3,"XXXA1s3",nx+1,ny+1,nz+1,4,1)

!-----------------------------------------------------------------------
!
!  Trapezoidal step:
!
!-----------------------------------------------------------------------

  DO k=0,nz
    DO j=0,ny
      DO i=0,nx
        s3(i,j,k)=(s3(i,j,k)+s2(i,j,k))*0.5
      END DO
    END DO
  END DO

  CALL advflxs(nx,ny,nz,nt,dx,dy,dz,dtbig,                              &
               ustr,vstr,wstr,s3,fctorderopt,                           &
               fluxx,fluxy,fluxz,                                       &
               s1,tem1)  ! Attention: s1 is used as a work array here!
!
!-----------------------------------------------------------------------
!
!  Low order scheme: Donnor Cell.
!  Calculate the donnor-cell fluxes
!
!-----------------------------------------------------------------------

  tem = dtbig*dy*dz*0.5

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        cx(i,j,k)=tem*                                                  &
            ((ustr(i,j,k)+ABS(ustr(i,j,k)))*s2(i-1,j,k)                 &
            +(ustr(i,j,k)-ABS(ustr(i,j,k)))*s2(i  ,j,k))
        fluxx(i,j,k)=fluxx(i,j,k)-cx(i,j,k)
      END DO
    END DO
  END DO
!
  tem = dtbig*dx*dz*0.5

  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        cy(i,j,k)=tem*                                                  &
            ((vstr(i,j,k)+ABS(vstr(i,j,k)))*s2(i,j-1,k)                 &
            +(vstr(i,j,k)-ABS(vstr(i,j,k)))*s2(i,j  ,k))
        fluxy(i,j,k)=fluxy(i,j,k)-cy(i,j,k)
      END DO
    END DO
  END DO
!
  tem = dtbig*dx*dy*0.5

  DO k=1,nz
    DO j=1,ny-1
      DO i=1,nx-1
        cz(i,j,k)=tem*                                                  &
            ((wstr(i,j,k)+ABS(wstr(i,j,k)))*s2(i,j,k-1)                 &
            +(wstr(i,j,k)-ABS(wstr(i,j,k)))*s2(i,j,k  ))
        fluxz(i,j,k)=fluxz(i,j,k)-cz(i,j,k)
      END DO
    END DO
  END DO
!
!  Compute the low order time advanced solution for scalar
!
  rdt=1.0/dtbig

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tema = ((cx(i+1,j,k)-cx(i,j,k)                                  &
                +cy(i,j+1,k)-cy(i,j,k))*mapfct2(i,j)                    &
                +cz(i,j,k+1)-cz(i,j,k))*rdxyz
        s3(i,j,k)=s2(i,j,k)-tema*rhostri(i,j,k)

        sadv(i,j,k) = sadv(i,j,k)+ tema * rdt

      END DO
    END DO
  END DO
!test tmp: call test_dump (sadv,"XXXfct_sadv",nx,ny,nz,0,0)

!-----------------------------------------------------------------------
!
!  Determine the flux limitor by calling fctlim.
!  fctlim returns the limited/corrected anti-diffusive fluxes
!  which are added to the lower order fluxes.
!
!c#######################################################################

!  IF (mp_opt > 0) THEN
!    CALL acct_interrupt(mp_acct)
!    CALL mpsendextew(s2,nx,ny,nz,ebc,wbc,mptag1,tem1)
!    CALL mpsendextew(s3,nx,ny,nz,ebc,wbc,mptag2,tem1)
!    CALL mprecvextew(s2,nx,ny,nz,ebc,wbc,mptag1,tem1)
!    CALL mprecvextew(s3,nx,ny,nz,ebc,wbc,mptag2,tem1)
!    CALL mpsendextns(s2,nx,ny,nz,nbc,sbc,mptag1,tem1)
!    CALL mpsendextns(s3,nx,ny,nz,nbc,sbc,mptag2,tem1)
!    CALL mprecvextns(s2,nx,ny,nz,nbc,sbc,mptag1,tem1)
!    CALL mprecvextns(s3,nx,ny,nz,nbc,sbc,mptag2,tem1)
!  END IF
  CALL acct_interrupt(bc_acct)
  CALL extndsbc(s2,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL extndsbc(s3,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL acct_stop_inter
!test tmp: call test_dump (s2,"XXXA2s2",nx+1,ny+1,nz+1,4,1)
!test tmp: call test_dump (s3,"XXXA2s3",nx+1,ny+1,nz+1,4,1)

  CALL fctlim(nx,ny,nz,dxyz,s3,s2,rhostr,mapfct2,fluxx,fluxy,fluxz,     &
              cx,cy,cz,tem1,tem2)

!-----------------------------------------------------------------------
!
!  Add the flux-limitted adtidifusive fluxes.
!  Note the flux nomal to the lateral boundaries are not applied.
!
!c#######################################################################

  tem = 1.0/(dtbig*dx*dy*dz)

  DO k=2,nz-2

    DO j=1,ny-1
      DO i=2,nx-2
        sadv(i,j,k)=sadv(i,j,k)+(fluxx(i+1,j,k)-fluxx(i,j,k))           &
                   *mapfct2(i,j)*tem
      END DO
    END DO

    DO j=2,ny-2
      DO i=1,nx-1
        sadv(i,j,k)=sadv(i,j,k)+(fluxy(i,j+1,k)-fluxy(i,j,k))           &
                   *mapfct2(i,j)*tem
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        sadv(i,j,k)=sadv(i,j,k)+(fluxz(i,j,k+1)-fluxz(i,j,k))*tem
      END DO
    END DO

  END DO
!test tmp: call test_dump (fluxx,"XXXfct_fluxx",nx,ny,nz,1,0)
!test tmp: call test_dump (fluxy,"XXXfct_fluxy",nx,ny,nz,2,0)
!test tmp: call test_dump (fluxz,"XXXfct_fluxz",nx,ny,nz,3,0)
!test tmp: call test_dump (sadv,"XXX1fct_sadv",nx,ny,nz,0,0)

  RETURN
END SUBROUTINE advctsfct
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE ADVFLXS                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE advflxs(nx,ny,nz,nt,dx,dy,dz,dt,                             & 3,9
           ustr,vstr,wstr,s,advodropt,                                  &
           fluxx,fluxy,fluxz,                                           &
           tem1_0,mp_tem)
!
!----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate 2nd- and 4th-order fluxes for the advection of a scalar.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue.
!  7/18/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!    nt       number of the time level
!    dx       Grid interval in x direction
!    dy       Grid interval in x direction
!    dz       Grid interval in x direction
!    dt       Time step size
!
!    ustr     u*rhostr
!    vstr     v*rhostr
!    wstr     wcont*rhostr
!    s        Scalar whose advective fluxes are to be found
!
!  OUTPUT:
!
!    fluxx    Advective flux in x direction.
!    fluxy    Advective flux in y direction.
!    fluxz    Advective flux in z direction.
!
!  Work array:
!    tem1_0   Work array
!    mp_tem   Work array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
  INTEGER :: nt              ! The no. of t-levels of t-dependent arrays.
  REAL :: dx,dy,dz,dt

  REAL :: ustr  (nx,ny,nz)   ! u*rhostr
  REAL :: vstr  (nx,ny,nz)   ! v*rhostr
  REAL :: wstr  (nx,ny,nz)   ! wcont*rhostr
  REAL :: s(0:nx,0:ny,0:nz)  ! Scalar whose advective fluxes are to be found
!
  REAL :: fluxx (nx,ny,nz)   ! Advective flux in x direction.
  REAL :: fluxy (nx,ny,nz)   ! Advective flux in y direction.
  REAL :: fluxz (nx,ny,nz)   ! Advective flux in z direction.

  REAL :: tem1_0(0:nx,0:ny,0:nz) ! Work array
  REAL :: mp_tem(MAX(nx+1,ny+1)*(nz+1))  ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  REAL :: tem, rdxyz,dxyz, tem43,tem16

  INTEGER :: advodropt

  INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'bndry.inc'     ! the boundary flags
  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  dxyz=dx*dy*dz
  rdxyz=1./(dxyz)

  tem43 = 4.0/3.0
  tem16 = 1.0/6.0
!
!-----------------------------------------------------------------------
!
!  Compute the second-order advective fluxes
!
!-----------------------------------------------------------------------
!
  tem = dt*dy*dz*0.5
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        fluxx(i,j,k)=tem*ustr(i,j,k)*(s(i,j,k)+s(i-1,j,k))
      END DO
    END DO
  END DO
!test tmp: call test_dump (fluxx,"XXXadvflxs_fluxx",nx,ny,nz,1,0)

  IF( advodropt == 2) THEN  ! For fourth order fluxes

    tem = dt*dy*dz*0.25
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1_0(i,j,k)=                                                &
              tem*(ustr(i,j,k)+ustr(i+1,j,k))*(s(i+1,j,k)+s(i-1,j,k))
        END DO
      END DO
    END DO

!    IF (mp_opt > 0) THEN
!      CALL acct_interrupt(mp_acct)
!      CALL mpsendextew(tem1_0,nx,ny,nz,ebc,wbc,mptag,mp_tem)
!      CALL mprecvextew(tem1_0,nx,ny,nz,ebc,wbc,mptag,mp_tem)
!    END IF
    CALL acct_interrupt(bc_acct)
    CALL extndsbc(tem1_0,nx,ny,nz,1,ebc,wbc,0,0,0,0)
    CALL acct_stop_inter
!test tmp: call test_dump2(tem1_0,"XXXA3tem1_0",nx+1,ny+1,nz+1,1,nx+1,2,ny,2,nz)

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx
          fluxx(i,j,k)=tem43*fluxx(i,j,k)                               &
                      -tem16*(tem1_0(i-1,j,k)+tem1_0(i,j,k))
        END DO
      END DO
    END DO
!test tmp: call test_dump (fluxx,"XXX1advflxs_fluxx",nx,ny,nz,1,0)

  END IF

  tem = dt*dx*dz*0.5
  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        fluxy(i,j,k)=tem*vstr(i,j,k)*(s(i,j,k)+s(i,j-1,k))
      END DO
    END DO
  END DO
!test tmp: call test_dump (fluxy,"XXXadvflxs_fluxy",nx,ny,nz,2,0)

  IF( advodropt == 2) THEN  ! For fourth order fluxes

    tem = dt*dx*dz*0.25
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1_0(i,j,k)=                                                &
              tem*(vstr(i,j,k)+vstr(i,j+1,k))*(s(i,j+1,k)+s(i,j-1,k))
        END DO
      END DO
    END DO

!    IF (mp_opt > 0) THEN
!      CALL acct_interrupt(mp_acct)
!      CALL mpsendextns(tem1_0,nx,ny,nz,nbc,sbc,mptag,mp_tem)
!      CALL mprecvextns(tem1_0,nx,ny,nz,nbc,sbc,mptag,mp_tem)
!    END IF
    CALL acct_interrupt(bc_acct)
    CALL extndsbc(tem1_0,nx,ny,nz,1,0,0,nbc,sbc,0,0)
    CALL acct_stop_inter
!test tmp: call test_dump (tem1_0,"XXXA4tem1_0",nx+1,ny+1,nz+1,4,1)

    DO k=1,nz-1
      DO j=1,ny
        DO i=1,nx-1
          fluxy(i,j,k)=tem43*fluxy(i,j,k)                               &
                      -tem16*(tem1_0(i,j-1,k)+tem1_0(i,j,k))
        END DO
      END DO
    END DO
!test tmp: call test_dump (fluxy,"XXX1advflxs_fluxy",nx,ny,nz,2,0)

  END IF

  tem = dt*dx*dy*0.5
  DO k=1,nz
    DO j=1,ny-1
      DO i=1,nx-1
        fluxz(i,j,k)=tem*wstr(i,j,k)*(s(i,j,k)+s(i,j,k-1))
      END DO
    END DO
  END DO
!test tmp: call test_dump (fluxz,"XXXadvflxs_fluxz",nx,ny,nz,3,0)

  IF( advodropt == 2) THEN  ! For fourth order fluxes

    tem = dt*dx*dy*0.25
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem1_0(i,j,k)=                                                &
              tem*(wstr(i,j,k)+wstr(i,j,k+1))*(s(i,j,k+1)+s(i,j,k-1))
        END DO
      END DO
    END DO

    CALL acct_interrupt(bc_acct)
    CALL extndsbc(tem1_0,nx,ny,nz,1,0,0,0,0,tbc,bbc)
    CALL acct_stop_inter
!test tmp: call test_dump (tem1_0,"XXXA5tem1_0",nx+1,ny+1,nz+1,4,1)

    DO k=1,nz
      DO j=1,ny-1
        DO i=1,nx-1
          fluxz(i,j,k)=tem43*fluxz(i,j,k)                               &
                      -tem16*(tem1_0(i,j,k-1)+tem1_0(i,j,k))
        END DO
      END DO
    END DO
!test tmp: call test_dump (fluxz,"XXX1advflxs_fluxz",nx,ny,nz,3,0)

  END IF

  RETURN
END SUBROUTINE advflxs
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE ADVMPDCD                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Coastal Meteorology Research Program (CMRP)      ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE advmpdcd(nx,ny,nz,nt,dtbig1,                                 & 2,4
           ustr,vstr,wstr,mapfct2,s1,s2,rhostr,                         &
           sadv,                                                        &
           fluxx,fluxy,fluxz,fout,bout,mp_tem)
!
!----------------------------------------------------------------------
!
!  PURPOSE:
!
!  MPDCD (Multidimensional Positive Definite Centered Difference
!  Scheme) positive definite advection for scalars
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Fanyou Kong (CMRP)
!  07/10/1997.
!
!  11/20/1997 (Fanyou Kong -- CMRP)
!  - Added map factor into MPDCD schemes
!  - Added Loop 50 to record possible negative scalar values in "s1"
!    and then force them positive to prevent possible floating
!    point error, but report the situation in output file
!
!  7/17/1998 (Ming Xue)
!    Modified the formulation of fluxes. sadv is now accumulative.
!
!  9/12/1999 (Pengfei Zhang & Ming Xue)
!    Added a constant "eps" in the denominator of an equation to
!    prevent overflow.
!
!  2000/09/15 (Gene Bassett)
!  Removed dx,dy,dz & fctorderopt from argument list.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!    nt       number of the time level
!    dx       Grid interval in x direction
!    dy       Grid interval in x direction
!    dz       Grid interval in x direction
!    dtbig1   Time step size
!
!    ustr     u*rhostr
!    vstr     v*rhostr
!    wstr     wcont*rhostr
!
!    mapfct2  map factor for scalar **2
!    s1       scalar field at a given time n-1 level (scalar units)
!    s2       scalar field at a given time n level (scalar units)
!    rhostr   rhobar*j3
!    fctorderopt Option parameter for 2nd or 4th order high-order fluxes
!
!  OUTPUT:
!
!    sadv     advection term              [(kg/m**3) scalar unit/s]
!
!  WORK ARRAYS:
!
!    fluxx     Temporary work array.
!    fluxy     Temporary work array.
!    fluxz     Temporary work array.
!    bout      Temporary work array.
!    fout      Temporary work array.
!    mp_tem    Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE
!
  INTEGER :: nx,ny,nz
  INTEGER :: nt              ! The no. of t-levels of t-dependent arrays.
  REAL :: dtbig1

  REAL :: ustr  (nx,ny,nz)   ! u*rhostr
  REAL :: vstr  (nx,ny,nz)   ! v*rhostr
  REAL :: wstr  (nx,ny,nz)   ! wcont*rhostr

  REAL :: mapfct2(nx,ny)     ! map factor for scalar **2
  REAL :: s1    (0:nx,0:ny,0:nz)! Time lev. 1 of a scalar field to be advected.
  REAL :: s2    (0:nx,0:ny,0:nz)! Time lev. 2 of a scalar field to be advected.
  REAL :: rhostr(nx,ny,nz)   ! rhobar*j3

  REAL :: sadv  (nx,ny,nz)   ! advection term

!
  REAL :: fluxx (nx,ny,nz)   ! temporary work array
  REAL :: fluxy (nx,ny,nz)   ! temporary work array
  REAL :: fluxz (nx,ny,nz)   ! temporary work array
  REAL :: fout  (nx,ny,nz)   ! temporary work array
  REAL :: bout  (0:nx,0:ny,0:nz)   ! temporary work array

  REAL :: mp_tem(MAX(nx+1,ny+1)*(nz+1))  ! Temporary message passing array.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, negnum
  REAL :: tem, eps
  REAL :: dtxy,dtxz,dtyz,rdxyz,dxyz

  INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'bndry.inc'     ! the boundary flags
  INCLUDE 'globcst.inc'       ! Global model control parameters
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  dtxy=dtbig1*dx*dy
  dtxz=dtbig1*dz*dx
  dtyz=dtbig1*dz*dy
  dxyz=dx*dy*dz
  rdxyz=1./(dxyz)
  eps=1.e-20
!
!-----------------------------------------------------------------------
!
!  Compute the second-order advective fluxes d(rhostr*u*s)/dx
!
!-----------------------------------------------------------------------
!
  CALL advflxs(nx,ny,nz,nt,dx,dy,dz,dtbig1,                             &
               ustr,vstr,wstr,s2,fctorderopt,                           &
               fluxx,fluxy,fluxz,                                       &
               bout,fout)
! Attention: bout & fout are used as a work arrays here!
!test tmp: call test_dump (fluxx,"XXX1fct5_fluxx",nx,ny,nz,1,0)
!test tmp: call test_dump (fluxy,"XXX1fct5_fluxy",nx,ny,nz,2,0)
!test tmp: call test_dump (fluxz,"XXX1fct5_fluxz",nx,ny,nz,3,0)

!  compute the total outgoing flux for each grid cell (FOUT)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        fout(i,j,k)=rdxyz*( mapfct2(i,j)*(                              &
                           (ABS(fluxx(i  ,j,k))-fluxx(i  ,j,k)) +       &
                           (ABS(fluxx(i+1,j,k))+fluxx(i+1,j,k)) +       &
                           (ABS(fluxy(i  ,j,k))-fluxy(i  ,j,k)) +       &
                           (ABS(fluxy(i,j+1,k))+fluxy(i,j+1,k)))+       &
                           (ABS(fluxz(i,j  ,k))-fluxz(i,j  ,k)) +       &
                           (ABS(fluxz(i,j,k+1))+fluxz(i,j,k+1)))
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1

        s1(i,j,k) = MAX(0.0, s1(i,j,k)) ! Set negative values to zero

        IF(fout(i,j,k) <= s1(i,j,k)*rhostr(i,j,k)) THEN
          bout(i,j,k) = 1.0
        ELSE
          bout(i,j,k) = (s1(i,j,k)*rhostr(i,j,k))/(fout(i,j,k)+eps)
        END IF

      END DO
    END DO
  END DO

!  IF (mp_opt > 0) THEN
!    CALL acct_interrupt(mp_acct)
!    CALL mpsendextew(bout,nx,ny,nz,ebc,wbc,mptag,mp_tem)
!    CALL mprecvextew(bout,nx,ny,nz,ebc,wbc,mptag,mp_tem)
!    CALL mpsendextns(bout,nx,ny,nz,nbc,sbc,mptag,mp_tem)
!    CALL mprecvextns(bout,nx,ny,nz,nbc,sbc,mptag,mp_tem)
!  END IF
  CALL acct_interrupt(bc_acct)
  CALL extndsbc(bout,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
  CALL acct_stop_inter
!test tmp: call test_dump (bout,"XXXAbout",nx+1,ny+1,nz+1,4,1)
!
!-----------------------------------------------------------------------
!
!  Application of the flux limiter to fluxes
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx
        IF(fluxx(i,j,k) > 0.0) THEN
          fluxx(i,j,k)=bout(i-1,j,k)*fluxx(i,j,k)
        ELSE
          fluxx(i,j,k)=bout(i,j,k)  *fluxx(i,j,k)
        END IF
      END DO
    END DO
  END DO
!
  DO k=1,nz-1
    DO j=1,ny
      DO i=1,nx-1
        IF(fluxy(i,j,k) > 0.0) THEN
          fluxy(i,j,k)=bout(i,j-1,k)*fluxy(i,j,k)
        ELSE
          fluxy(i,j,k)=bout(i,j,k)  *fluxy(i,j,k)
        END IF
      END DO
    END DO
  END DO
  DO k=1,nz
    DO j=1,ny-1
      DO i=1,nx-1
        IF(fluxz(i,j,k) > 0.0) THEN
          fluxz(i,j,k)=bout(i,j,k-1)*fluxz(i,j,k)
        ELSE
          fluxz(i,j,k)=bout(i,j,k)  *fluxz(i,j,k)
        END IF
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate advection term (sadv)
!
!-----------------------------------------------------------------------

!test tmp: call test_dump (sadv,"XXXfct5_sadv",nx,ny,nz,0,0)
  tem = 1.0/(dtbig1*dx*dy*dz)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        sadv(i,j,k)=sadv(i,j,k)+                                        &
                    ((fluxx(i+1,j,k)-fluxx(i,j,k) +                     &
                      fluxy(i,j+1,k)-fluxy(i,j,k))*mapfct2(i,j) +       &
                      fluxz(i,j,k+1)-fluxz(i,j,k))*tem
      END DO
    END DO
  END DO
!test tmp: call test_dump (fluxx,"XXXfct5_fluxx",nx,ny,nz,1,0)
!test tmp: call test_dump (fluxy,"XXXfct5_fluxy",nx,ny,nz,2,0)
!test tmp: call test_dump (fluxz,"XXXfct5_fluxz",nx,ny,nz,3,0)
!test tmp: call test_dump (sadv,"XXX1fct5_sadv",nx,ny,nz,0,0)

  RETURN
END SUBROUTINE advmpdcd
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE FCTLIM                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE fctlim(nx,ny,nz,dxyz,std,sa,rhostr,mapfct2,                  & 1
           fluxx,fluxy,fluxz,cx,cy,cz,rminus,rplus)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Perform 3-D flux limiting (Zalesak, 1979) on fluxx, fluxy, fluxz
!  which are fluxes in x, y, and z direction respectivly.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  1/1/1988
!
!  (Originally written for the sigma-coordinate model of
!  Xue and Thorpe, 1991).
!
!  MODIFICATION HISTORY
!
!  11/20/92  (Yifeng Tang)
!  Adapted for ARPS, taking into account of various boundary condition
!  options.
!
!  10/12/1996 (M. Xue)
!  Cleaned up and portions rewritten for ARPS.
!
!  11/20/1997 (Fanyou Kong - CMRP)
!  Added map factor
!
!  7/16/1998 (M. Xue and R. Richardson)
!  Significant changes made. The extrema in the advected
!  (instead of the previous rhostr multipled) field are now used
!  in determining the flux correction coefficients.
!  A one-dimensional prelimiter is included as an option,
!  although it's not recommended for common use.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    dxyz     dx*dy*dz
!
!    std      the scalar at n+1 computed from donnor cell fluxes
!    sa       the scalar at time n
!             std and sa are unchanged on return.
!    rhostr   rhobar * j3
!    mapfct2  Map factor for scalars **2
!    fluxx    the antidiffusive flux in x direction
!    fluxy    the antidiffusive flux in y direction
!    fluxz    the antidiffusive flux in z direction
!    cx       the donnor cell flux in x direction
!    cy       the donnor cell flux in y direction
!    cz       the donnor cell flux in z direction
!
!  OUTPUT:
!
!    fluxx    the limited antidiffusive flux in x direction
!    fluxy    the limited antidiffusive flux at y
!    fluxz    the limited antidiffusive flux at z
!
!  WORK ARRAYS:
!
!    rminus    the work array
!    rplus     the work array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: std(0:nx,0:ny,0:nz)
  REAL :: sa (0:nx,0:ny,0:nz)

  REAL :: rhostr(nx,ny,nz)
  REAL :: mapfct2(nx,ny)
  REAL :: fluxx (nx,ny,nz)
  REAL :: fluxy (nx,ny,nz)
  REAL :: fluxz (nx,ny,nz)
  REAL :: cx    (nx,ny,nz)
  REAL :: cy    (nx,ny,nz)
  REAL :: cz    (nx,ny,nz)
  REAL :: rminus(nx,ny,nz)
  REAL :: rplus (nx,ny,nz)
  REAL :: dxyz
  REAL :: qplus,qminus,pplus,pminus

  INTEGER :: fct1dlim  ! Switch for 1D prelimiter
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'bndry.inc'           ! the boundary flags
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  fct1dlim = 0

  IF( fct1dlim == 1) THEN
!
!-----------------------------------------------------------------------
!
!  Pre-limiting in x-direciton.
!
!-----------------------------------------------------------------------
!

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pplus=(MAX(0.,fluxx(i,j,k))-MIN(0.,fluxx(i+1,j,k)))           &
                 *mapfct2(i,j)
          qplus=( MAX(std(i-1,j,k),std(i+1,j,k),std(i,j,k),             &
                      sa (i-1,j,k),sa (i+1,j,k),sa (i,j,k))             &
                      -std(i,j,k) )*dxyz * rhostr(i,j,k)
          IF( pplus < 1.0E-20) THEN
            rplus(i,j,k)=0.0
          ELSE
            rplus(i,j,k)=MIN(1.,qplus/pplus )
          END IF

        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pminus=(MAX(0.,fluxx(i+1,j,k))-MIN(0.,fluxx(i,j,k)))          &
                  *mapfct2(i,j)
          qminus=-( MIN(std(i-1,j,k),std(i+1,j,k),std(i,j,k),           &
                        sa (i-1,j,k),sa (i+1,j,k),sa (i,j,k))           &
                       -std(i,j,k) )*dxyz * rhostr(i,j,k)

          IF( pminus < 1.e-20 ) THEN
            rminus(i,j,k)=0.0
          ELSE
            rminus(i,j,k)= MIN(1.,qminus/pminus)
          END IF

        END DO
      END DO
    END DO
!
!  Determine the correction coefficient C
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          IF( fluxx(i,j,k) >= 0.0) THEN
            cx(i,j,k)=MIN(rplus (i,j,k),rminus(i-1,j,k))
          ELSE
            cx(i,j,k)=MIN(rminus(i,j,k),rplus (i-1,j,k))
          END IF
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Pre-limiting in y-direciton.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pplus=(MAX(0.,fluxy(i,j,k))-MIN(0.,fluxy(i,j+1,k)))           &
                *mapfct2(i,j)
          qplus=( MAX(std(i,j,k),std(i,j-1,k),std(i,j+1,k),             &
                      sa (i,j,k),sa (i,j-1,k),sa (i,j+1,k))             &
                     -std(i,j,k) )*dxyz * rhostr(i,j,k)
          IF( pplus < 1.0E-20) THEN
            rplus(i,j,k)=0.0
          ELSE
            rplus(i,j,k)=MIN(1.,qplus/pplus )
          END IF
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pminus=(MAX(0.,fluxy(i,j+1,k))-MIN(0.,fluxy(i,j,k)))          &
                  *mapfct2(i,j)
          qminus=-( MIN(std(i,j,k),std(i,j-1,k),std(i,j+1,k),           &
                      sa (i,j,k),sa (i,j-1,k),sa (i,j+1,k))             &
                      -std(i,j,k) )*dxyz * rhostr(i,j,k)
          IF( pminus < 1.e-20 ) THEN
            rminus(i,j,k)=0.0
          ELSE
            rminus(i,j,k)= MIN(1.,qminus/pminus)
          END IF
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          IF( fluxy(i,j,k) >= 0.0) THEN
            cy(i,j,k)=MIN(rplus (i,j,k),rminus(i,j-1,k))
          ELSE
            cy(i,j,k)=MIN(rminus(i,j,k),rplus (i,j-1,k))
          END IF
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Pre-limiting in z-direciton.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pplus=(MAX(0.,fluxz(i,j,k))-MIN(0.,fluxz(i,j,k+1)))
          qplus=( MAX(std(i,j,k),std(i,j,k-1),std(i,j,k+1),             &
                      sa (i,j,k),sa (i,j,k-1),sa (i,j,k+1))             &
                     -std(i,j,k) )*dxyz * rhostr(i,j,k)
          IF( pplus < 1.0E-20) THEN
            rplus(i,j,k)=0.0
          ELSE
            rplus(i,j,k)=MIN(1.,qplus/pplus )
          END IF
        END DO
      END DO
    END DO

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pminus=(MAX(0.,fluxz(i,j,k+1))-MIN(0.,fluxz(i,j,k)))
          qminus=-( MIN(std(i,j,k),std(i,j,k-1),std(i,j,k+1),           &
                      sa (i,j,k),sa (i,j,k-1),sa (i,j,k+1))             &
                      -std(i,j,k) )*dxyz * rhostr(i,j,k)

          IF( pminus < 1.e-20 ) THEN
            rminus(i,j,k)=0.0
          ELSE
            rminus(i,j,k)= MIN(1.,qminus/pminus)
          END IF

        END DO
      END DO
    END DO

    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          IF( fluxz(i,j,k) >= 0.0) THEN
            cz(i,j,k)=MIN(rplus (i,j,k),rminus(i,j,k-1))
          ELSE
            cz(i,j,k)=MIN(rminus(i,j,k),rplus (i,j,k-1))
          END IF
        END DO
      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Apply pre-limiting factor to the antidiffusive fluxes.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=2,nx-1
          fluxx(i,j,k)=fluxx(i,j,k)*cx(i,j,k)
        END DO
      END DO
    END DO
!
    DO k=1,nz-1
      DO j=2,ny-1
        DO i=1,nx-1
          fluxy(i,j,k)=fluxy(i,j,k)*cy(i,j,k)
        END DO
      END DO
    END DO
!
    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          fluxz(i,j,k)=fluxz(i,j,k)*cz(i,j,k)
        END DO
      END DO
    END DO

  END IF  ! End of 1D prelimiting step

!-----------------------------------------------------------------------
!
!  Calculate total flux directing towards point (i,j,k) pplus
!  then inward flux amplification factor r(i,j,k)+
!  r+=min(1.,q+/p+).
!
!-----------------------------------------------------------------------

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        pplus=(MAX(0.,fluxx(i,j,k))-MIN(0.,fluxx(i+1,j,k))              &
              +MAX(0.,fluxy(i,j,k))-MIN(0.,fluxy(i,j+1,k)))             &
              *mapfct2(i,j)                                             &
              +MAX(0.,fluxz(i,j,k))-MIN(0.,fluxz(i,j,k+1))
        qplus=(MAX(std(i-1,j,k),std(i+1,j,k),std(i,j,k),                &
                   std(i,j-1,k),std(i,j+1,k),std(i,j,k-1),std(i,j,k+1), &
                   sa (i-1,j,k),sa (i+1,j,k),sa (i,j,k),                &
                   sa (i,j-1,k),sa (i,j+1,k),sa (i,j,k-1),sa (i,j,k+1)) &
              -std(i,j,k) )*dxyz * rhostr(i,j,k)
        IF( pplus < 1.0E-20) THEN
          rplus(i,j,k)=0.0
        ELSE
          rplus(i,j,k)=MIN(1.,qplus/pplus )
        END IF

      END DO
    END DO
  END DO

!-----------------------------------------------------------------------
!
!  Calculate total flux directing away from point (i,j,k) pminus
!  then outward flux amplification factor r(i,j,k)-
!  r- = min(1.,q-/p-).
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        pminus=(MAX(0.,fluxx(i+1,j,k))-MIN(0.,fluxx(i,j,k))             &
               +MAX(0.,fluxy(i,j+1,k))-MIN(0.,fluxy(i,j,k)))            &
               *mapfct2(i,j)                                            &
               +MAX(0.,fluxz(i,j,k+1))-MIN(0.,fluxz(i,j,k))
        qminus=-(MIN(std(i-1,j,k),std(i+1,j,k),std(i,j,k),              &
                    std(i,j-1,k),std(i,j+1,k),std(i,j,k-1),std(i,j,k+1), &
                    sa (i-1,j,k),sa (i+1,j,k),sa (i,j,k),               &
                    sa (i,j-1,k),sa (i,j+1,k),sa (i,j,k-1),sa (i,j,k+1)) &
                    -std(i,j,k) )*dxyz * rhostr(i,j,k)

        IF( pminus < 1.e-20 ) THEN
          rminus(i,j,k)=0.0
        ELSE
          rminus(i,j,k)= MIN(1.,qminus/pminus)
        END IF
      END DO
    END DO
  END DO
!
!  Determine the correction coefficient C
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=2,nx-1
        IF( fluxx(i,j,k) >= 0.0) THEN
          cx(i,j,k)=MIN(rplus (i,j,k),rminus(i-1,j,k))
        ELSE
          cx(i,j,k)=MIN(rminus(i,j,k),rplus (i-1,j,k))
        END IF
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=2,ny-1
      DO i=1,nx-1
        IF( fluxy(i,j,k) >= 0.0) THEN
          cy(i,j,k)=MIN(rplus (i,j,k),rminus(i,j-1,k))
        ELSE
          cy(i,j,k)=MIN(rminus(i,j,k),rplus (i,j-1,k))
        END IF
      END DO
    END DO
  END DO

  DO k=2,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        IF( fluxz(i,j,k) >= 0.0) THEN
          cz(i,j,k)=MIN(rplus (i,j,k),rminus(i,j,k-1))
        ELSE
          cz(i,j,k)=MIN(rminus(i,j,k),rplus (i,j,k-1))
        END IF
      END DO
    END DO
  END DO
!
!  Apply final correction factors to the antidiffusive fluxes
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=2,nx-1
        fluxx(i,j,k)=fluxx(i,j,k)*cx(i,j,k)
      END DO
    END DO
  END DO
!
  DO k=1,nz-1
    DO j=2,ny-1
      DO i=1,nx-1
        fluxy(i,j,k)=fluxy(i,j,k)*cy(i,j,k)
      END DO
    END DO
  END DO
!
  DO k=2,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        fluxz(i,j,k)=fluxz(i,j,k)*cz(i,j,k)
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE fctlim