ppm.f90

!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! Modifications:
!  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
!  P. Wautelet 20/06/2019: OpenACC: correct intent of some dummy variables
!  P. Wautelet 01/07/2019: OpenACC: optimisation of ppm_s0_x/y/z_d for GPU
!  P. Wautelet 18/07/2019: OpenACC: remove use of macros for dif2x/y/z
!-----------------------------------------------------------------
#ifdef MNH_OPENACC
!
! inline shuman with macro
! 
!#define dxf(PDXF,PA) PDXF(1:IIU-1,:,:) = PA(2:IIU,:,:) - PA(1:IIU-1,:,:) ; PDXF(IIU,:,:)    = PDXF(2*JPHEXT,:,:) ! DXF(PDXF,PA)
!#define  dyf(PDYF,PA) PDYF(:,1:IJU-1,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:); PDYF(:,IJU,:) = PDYF(:,2*JPHEXT,:) !   DYF(PDYF,PA)
!!#define dyf(PDYF,PA) PDYF(1:IIU,1:IJU-1,IKB:IKE) = PA(1:IIU,2:IJU,IKB:IKE) - PA(1:IIU,1:IJU-1,IKB:IKE); ! PDYF(1:IIU,IJU,IKB:IKE) = PDYF(1:IIU,2*JPHEXT,IKB:IKE) !   DYF(PDYF,PA)
!#define dzf(PDZF,PA) PDZF(:,:,1:IKU-1) = PA(:,:,2:IKU) - PA(:,:,1:IKU-1) ; PDZF(:,:,IKU) = -999. ! DZF(PDZF,PA)
!
!#define mxm(PMXM,PA) PMXM(2:IIU,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:) !  MXM(PMXM,PA)
!!#define mym(PMYM,PA) PMYM(1:IIU,2:IJU,IKB:IKE) = 0.5*( PA(1:IIU,2:IJU,IKB:IKE)+PA(1:IIU,1:IJU-1,IKB:IKE) ) ; ! PMYM(1:IIU,1,IKB:IKE) = PMYM(1:IIU,IJU-2*JPHEXT+1,IKB:IKE) !  MYM(PMYM,PA)
!#define mzm(PMZM,PA) PMZM(:,:,2:IKU) = 0.5*( PA(:,:,2:IKU)+PA(:,:,1:IKU-1) ) ; PMZM(:,:,1)    = -999. !  MZM(PMZM,PA)
!#define mym(PMYM,PA) PMYM(:,2:IJU,:) = 0.5*( PA(:,2:IJU,:)+PA(:,1:IJU-1,:) ) ; PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:) !  MYM(PMYM,PA)
!
! #define dif2x(DQ,PQ) DQ(IIB:IIE,:,:)=0.5*(PQ(IIB+1:IIE+1,:,:)-PQ(IIB-1:IIE-1,:,:));\
! DQ(IIB-1,:,:)=0.5*(PQ(IIB,:,:)-PQ(IIE-1,:,:));\
! DQ(IIE+1,:,:)=0.5*(PQ(IIB+1,:,:)-PQ(IIE,:,:)) ! DIF2X(DQ,PQ)
!
! #define dif2y(DQ,PQ) DQ(1:IIU,IJB:IJE,IKB:IKE) = 0.5*(PQ(1:IIU,IJB+1:IJE+1,IKB:IKE) - PQ(1:IIU,IJB-1:IJE-1,IKB:IKE)) ; !
! ! DQ(1:IIU,IJB-1,IKB:IKE) = 0.5*(PQ(1:IIU,IJB,IKB:IKE) - PQ(1:IIU,IJE-1,IKB:IKE)) ; \
! DQ(1:IIU,IJE+1,IKB:IKE) = 0.5*(PQ(1:IIU,IJB+1,IKB:IKE) - PQ(1:IIU,IJE,IKB:IKE)) ! DIF2Y(DQ,PQ)
!
! #define dif2z(DQ,PQ) DQ(:,:,IKB:IKE) = 0.5*(PQ(:,:,IKB+1:IKE+1) - PQ(:,:,IKB-1:IKE-1)) ; \
! DQ(:,:,IKB-1) = -DQ(:,:,IKB) ;\
! DQ(:,:,IKE+1) = -DQ(:,:,IKE) ! DIF2Z(DQ,PQ)
!
#endif
!     ###############
      MODULE MODI_PPM
!     ###############
!
INTERFACE
!
#ifndef MNH_OPENACC
FUNCTION PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
               RESULT(PR)
#else
SUBROUTINE PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
#endif
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_01_X
#else
END SUBROUTINE PPM_01_X
#endif
!
!
#ifndef MNH_OPENACC
FUNCTION PPM_01_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP) &
               RESULT(PR)
#else
SUBROUTINE PPM_01_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
#endif
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_01_Y
#else
END SUBROUTINE PPM_01_Y
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_01_Z(KGRID, PSRC, PCR, PRHO, PTSTEP) RESULT(PR)
#else
SUBROUTINE PPM_01_Z(KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
#endif
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
#ifndef MNH_OPENACC
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC    ! variable at t
#else
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
#endif
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_01_Z
#else
END SUBROUTINE PPM_01_Z
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_S0_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
               RESULT(PR)
#else
SUBROUTINE PPM_S0_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP &
                   , PR)
#endif
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_S0_X
#else
END SUBROUTINE PPM_S0_X
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_S0_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP) &
               RESULT(PR)
#else
SUBROUTINE PPM_S0_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PTSTEP &
                   , PR)
#endif
!
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_S0_Y
#else
END SUBROUTINE PPM_S0_Y
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP) &
               RESULT(PR)
#else
SUBROUTINE PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP &
                   , PR)
#endif
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_S0_Z
#else
END SUBROUTINE PPM_S0_Z
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_S1_X(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
                        PTSTEP) RESULT(PR)
#else
SUBROUTINE PPM_S1_X(HLBCX, KGRID, PSRC, PCR, PRHO, PRHOT, &
                        PTSTEP, PR)
#endif
!
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOT ! density at t+dt
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_S1_X
#else
END SUBROUTINE PPM_S1_X
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_S1_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PRHOT, &
                        PTSTEP) RESULT(PR)
#else
SUBROUTINE PPM_S1_Y(HLBCY, KGRID, PSRC, PCR, PRHO, PRHOT, &
                        PTSTEP, PR)
#endif
!
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! X direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOT ! density at t+dt
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_S1_Y
#else
END SUBROUTINE PPM_S1_Y
#endif
!
#ifndef MNH_OPENACC
FUNCTION PPM_S1_Z(KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP) &
                        RESULT(PR)
#else
SUBROUTINE PPM_S1_Z(KGRID, PSRC, PCR, PRHO, PRHOT, PTSTEP, &
                    PR)
#endif
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
#ifndef MNH_OPENACC
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC    ! variable at t
#else
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
#endif
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOT ! density at t+dt
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
#endif
!
#ifndef MNH_OPENACC
END FUNCTION PPM_S1_Z
#else
END SUBROUTINE PPM_S1_Z
#endif
!
END INTERFACE
!
END MODULE MODI_PPM
!
!
!-------------------------------------------------------------------------------
!
#ifdef MNH_OPENACC
!     ########################################################################
!!$      FUNCTION PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
!!$               RESULT(PR)
      SUBROUTINE  PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP, PR)
!     ########################################################################

  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
  USE MODE_MNH_ZWORK, ONLY : IIU,IJU,IKU

  IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR

INTEGER :: IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG

!$acc data present( PSRC, PCR, PRHO, PR )

        CALL  MNH_GET_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)

        CALL  PPM_01_X_D(IIU,IJU,IKU,HLBCX, KGRID, &
                     & PSRC, PCR, PRHO, PTSTEP, PR, &
                     & ZT3D(:,:,:,IZQL),ZT3D(:,:,:,IZQR),ZT3D(:,:,:,IZDQ),ZT3D(:,:,:,IZQ6), &
                     & ZT3D(:,:,:,IZDMQ),ZT3D(:,:,:,IZQL0),ZT3D(:,:,:,IZQR0), ZT3D(:,:,:,IZQ60), &
                     & ZT3D(:,:,:,IZFPOS),ZT3D(:,:,:,IZFNEG) )

        CALL  MNH_REL_ZT3D(IZQL,IZQR,IZDQ,IZQ6,IZDMQ,IZQL0,IZQR0,IZQ60,IZFPOS,IZFNEG)

!$acc end data
!
CONTAINS
!
!     ########################################################################
        SUBROUTINE  PPM_01_X_D(IIU,IJU,IKU,HLBCX, KGRID, &
                     & PSRC, PCR, PRHO, PTSTEP, PR, &
                     & ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG)

!     ########################################################################
#else
!     ########################################################################
      FUNCTION PPM_01_X(HLBCX, KGRID, PSRC, PCR, PRHO, PTSTEP) &
               RESULT(PR)
!     ########################################################################
#endif
!!
!!****  PPM_01_X - PPM_01 fully monotonic PPM advection scheme in X direction
!!                 Colella notation
!!
!!    MODIFICATIONS
!!    -------------
!!    
!!    11.5.2006.  T. Maric - original version
!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
!!      J.Escobar 28/06/2018: limit computation on TAB(:,IJS:IJN,:) to avoid unneeded NaN
!!      J.Escobr  16/07/2018: still NaN pb => reintroduce initialization of temporary local array
!!
!-------------------------------------------------------------------------------
!
USE MODD_CONF

USE MODE_ll
use mode_mppdb
#ifdef MNH_OPENACC
use mode_msg
#endif

#if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
USE MODI_BITREP
#endif
#ifdef MNH_BITREP_OMP
USE MODI_BITREPZ
#endif
USE MODI_GET_HALO
#ifndef MNH_OPENACC
USE MODI_SHUMAN
#else
USE MODI_SHUMAN_DEVICE
#endif
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
#ifdef MNH_OPENACC
INTEGER                        , INTENT(IN) :: IIU,IJU,IKU
#endif
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
!
INTEGER,                INTENT(IN)  :: KGRID   ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: PSRC    ! variable at t
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCR     ! Courant number
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHO  ! density
REAL,                   INTENT(IN)  :: PTSTEP  ! Time step 
!
! output source term
#ifndef MNH_OPENACC
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(IIU,IJU,IKU), INTENT(OUT) :: PR
#endif
!
!*       0.2   Declarations of local variables :
!
INTEGER:: IIB,IJB    ! Begining useful area in x,y,z directions
INTEGER:: IIE,IJE    ! End useful area in x,y,z directions
!
integer :: ji, jj, jk
#ifndef MNH_OPENACC
integer :: iiu, iju, iku
! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL,ZQR
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZDQ,ZQ6
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZDMQ
!
! extra variables for the initial guess of parabolae parameters
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZQL0,ZQR0,ZQ60
!
! advection fluxes
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
!
!BEG JUAN PPM_LL
INTEGER                          :: IJS,IJN
!END JUAN PPM_LL
#else
INTEGER                          :: I,J,K 
!
!!$!
!!$! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
REAL , DIMENSION(IIU,IJU,IKU) :: &
  ZQL,ZQR, ZDQ,ZQ6, ZDMQ &
!!$!
!!$! extra variables for the initial guess of parabolae parameters
  , ZQL0,ZQR0,ZQ60 &
!!$!
!!$! advection fluxes
  , ZFPOS, ZFNEG
!
INTEGER                :: IJS,IJN
#endif
LOGICAL                :: GWEST , GEAST
!-------------------------------------------------------------------------------
!
#ifdef MNH_BITREP_OMP
CALL SBR_FZ(PSRC(:,:,:))
#endif
!
!$acc data present( PSRC, PCR, PRHO, PR , &
!$acc &             ZQL, ZQR, ZDQ, ZQ6, ZDMQ, ZQL0, ZQR0, ZQ60, ZFPOS, ZFNEG )

IF (MPPDB_INITIALIZED) THEN
  !Check all IN arrays
  CALL MPPDB_CHECK(PCR, "PPM_01_X beg:PCR")
  CALL MPPDB_CHECK(PRHO,"PPM_01_X beg:PRHO")
  !Check all INOUT arrays
  CALL MPPDB_CHECK(PSRC,"PPM_01_X beg:PSRC")
END IF
!
!*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
!                 ------------------------------
!
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IJS=IJB
IJN=IJE
!
GWEST = LWEST_ll()
GEAST = LEAST_ll()
!
!BEG JUAN PPM_LL
!
!*              initialise & update halo & halo2 for PSRC
!
#ifndef MNH_OPENACC
iiu = size( PSRC, 1 )
iju = size( PSRC, 2 )
iku = size( PSRC, 3 )

CALL GET_HALO(PSRC, HNAME='PSRC')
!
PR   (:,:,:) = PSRC(:,:,:)
ZQL  (:,:,:) = PSRC(:,:,:)
ZQR  (:,:,:) = PSRC(:,:,:)
ZDQ  (:,:,:) = PSRC(:,:,:)
ZQ6  (:,:,:) = PSRC(:,:,:)
ZDMQ (:,:,:) = PSRC(:,:,:)
ZQL0 (:,:,:) = PSRC(:,:,:)
ZQR0 (:,:,:) = PSRC(:,:,:)
ZQ60 (:,:,:) = PSRC(:,:,:)
ZFPOS(:,:,:) = PSRC(:,:,:)
ZFNEG(:,:,:) = PSRC(:,:,:)
#else
CALL GET_HALO_D(PSRC,HDIR="01_X", HNAME='PSRC')
!
!$acc kernels 
!$mnh_do_concurrent (ji=1:iiu,jj=1:iju,jk=1:iku)
        PR   (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZQL  (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZQR  (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZDQ  (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZQ6  (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZDMQ (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZQL0 (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZQR0 (ji, jj, jk ) = PSRC(ji, jj, jk )
        ZQ60 (ji, jj, jk ) = PSRC(ji, jj, jk )
!$mnh_end_do()
!
#if 0
ZFPOS(:,1:IJS,:)=PSRC(:,1:IJS,:)
ZFNEG(:,1:IJS,:)=PSRC(:,1:IJS,:)
ZFPOS(:,IJN:,:)=PSRC(:,IJN:,:)
ZFNEG(:,IJN:,:)=PSRC(:,IJN:,:)
#else
ZFPOS(:,:,:) = PSRC(:,:,:)
ZFNEG(:,:,:) = PSRC(:,:,:)
#endif
!$acc end kernels
#endif
!
!-------------------------------------------------------------------------------
!
SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
!
!        1.1   CYCLIC BOUNDARY CONDITIONS IN X DIRECTION
!              -----------------------------------------
!
CASE ('CYCL','WALL')          ! In that case one must have HLBCX(1) == HLBCX(2)
#ifdef MNH_OPENACC
  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_X', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
#endif
!
! calculate dmq
   ZDMQ = DIF2X(PSRC)
!
! monotonize the difference followinq eq. 5 in Lin94
!
!BEG JUAN PPM_LL01
!
!  ZDMQ(i) = Fct[ ZDMQ(i),PSRC(i),PSRC(i-1),PSRC(i+1) ]
!
   ZDMQ(IIB:IIE,IJS:IJN,:) = &
        SIGN( (MIN( ABS(ZDMQ(IIB:IIE,IJS:IJN,:)),2.0*(PSRC(IIB:IIE,IJS:IJN,:) - &
        MIN(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:))),    &
        2.0*(MAX(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:)) - &
        PSRC(IIB:IIE,IJS:IJN,:)) )), ZDMQ(IIB:IIE,IJS:IJN,:) )
!
!  WEST BOUND
!
!!$   ZDMQ(IIB-1,:,:) = & 
!!$        SIGN( (MIN( ABS(ZDMQ(IIB-1,:,:)), 2.0*(PSRC(IIB-1,:,:) - &
!!$        MIN(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:))),   &
!!$        2.0*(MAX(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:)) - &
!!$        PSRC(IIB-1,:,:)) )), ZDMQ(IIB-1,:,:) )
!
!  EAST BOUND
!
!!$   ZDMQ(IIE+1,:,:) = &
!!$        SIGN( (MIN( ABS(ZDMQ(IIE+1,:,:)), 2.0*(PSRC(IIE+1,:,:) - &
!!$        MIN(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:))),  &
!!$        2.0*(MAX(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:)) - &
!!$        PSRC(IIE+1,:,:)) )), ZDMQ(IIE+1,:,:) )
!
!  update ZDMQ HALO before next/further  utilisation 
!
   CALL  GET_HALO(ZDMQ, HNAME='ZDMQ')
!
! calculate qL and qR with the modified dmq
!
!  ZQL0(i) = Fct[ PSRC(i),PSRC(i-1),ZDMQ(i),ZDMQ(i-1) ]
!
   ZQL0(IIB:IIE+1,IJS:IJN,:) = 0.5*(PSRC(IIB:IIE+1,IJS:IJN,:) + PSRC(IIB-1:IIE,IJS:IJN,:)) - &
        (ZDMQ(IIB:IIE+1,IJS:IJN,:) - ZDMQ(IIB-1:IIE,IJS:IJN,:))/6.0
!
   CALL  GET_HALO(ZQL0, HNAME='ZQL0')
!  
!  WEST BOUND
!
!!$   ZQL0(IIB-1,:,:) = ZQL0(IIE,:,:) JUAN PPMLL01
!
   ZQR0(IIB-1:IIE,IJS:IJN,:) = ZQL0(IIB:IIE+1,IJS:IJN,:)
!
   CALL  GET_HALO(ZQR0, HNAME='ZQR0')
!
!  EAST BOUND
!
!!$   ZQR0(IIE+1,:,:) = ZQR0(IIB,:,:) JUAN PPMLL01
!
! determine initial coefficients of the parabolae
!
   ZDQ(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:) - ZQL0(:,IJS:IJN,:)
   ZQ60(:,IJS:IJN,:) = 6.0*(PSRC(:,IJS:IJN,:) - 0.5*(ZQL0(:,IJS:IJN,:) + ZQR0(:,IJS:IJN,:)))
!
! initialize final parabolae parameters
!
   ZQL(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:)
   ZQR(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:)
   ZQ6(:,IJS:IJN,:) = ZQ60(:,IJS:IJN,:) 
!
! eliminate over and undershoots and create qL and qR as in Lin96
!
   WHERE ( ZDMQ(:,IJS:IJN,:) == 0.0 )
      ZQL(:,IJS:IJN,:) = PSRC(:,IJS:IJN,:)
      ZQR(:,IJS:IJN,:) = PSRC(:,IJS:IJN,:)
      ZQ6(:,IJS:IJN,:) = 0.0
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) < -(ZDQ(:,IJS:IJN,:))**2 )
#else
   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) < -BR_P2(ZDQ(:,IJS:IJN,:)) )
#endif
      ZQ6(:,IJS:IJN,:) = 3.0*(ZQL0(:,IJS:IJN,:) - PSRC(:,IJS:IJN,:))
      ZQR(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:) - ZQ6(:,IJS:IJN,:)
      ZQL(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:)
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) > (ZDQ(:,IJS:IJN,:))**2 )
#else
   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) > BR_P2(ZDQ(:,IJS:IJN,:)) )
#endif
      ZQ6(:,IJS:IJN,:) = 3.0*(ZQR0(:,IJS:IJN,:) - PSRC(:,IJS:IJN,:))
      ZQL(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:) - ZQ6(:,IJS:IJN,:)
      ZQR(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:)
   END WHERE
!
! recalculate coefficients of the parabolae
!
   ZDQ(:,IJS:IJN,:) = ZQR(:,IJS:IJN,:) - ZQL(:,IJS:IJN,:)
!
! and finally calculate fluxes for the advection
!
!  ZFPOS(i) = Fct[ ZQR(i-1),PCR(i),ZDQ(i-1),ZQ6(i-1) ]
!
   ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*PCR(IIB:IIE+1,IJS:IJN,:) * &
        (ZDQ(IIB-1:IIE,IJS:IJN,:) - (1.0 - 2.0*PCR(IIB:IIE+1,IJS:IJN,:)/3.0)        &
        * ZQ6(IIB-1:IIE,IJS:IJN,:))
!
   CALL GET_HALO(ZFPOS, HNAME='ZFPOS')
!
!  WEST BOUND
!
! PPOSX(IIB-1,:,:) is not important for the calc of advection so 
! we set it to 0
!!$   ZFPOS(IIB-1,:,:) = 0.0 JUANPPMLL01
!
   ZFNEG(:,IJS:IJN,:) = ZQL(:,IJS:IJN,:) - 0.5*PCR(:,IJS:IJN,:) *      &            
        ( ZDQ(:,IJS:IJN,:) + (1.0 + 2.0*PCR(:,IJS:IJN,:)/3.0) * ZQ6(:,IJS:IJN,:) )
!
   CALL GET_HALO(ZFNEG, HNAME='ZFNEG')
!
! advect the actual field in X direction by U*dt
!
#ifndef MNH_OPENACC
   PR = DXF( PCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                             ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
#else
  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_X', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
#endif
   CALL GET_HALO(PR, HNAME='PR')
!
!
!*       1.2    NON-CYCLIC BOUNDARY CONDITIONS IN THE X DIRECTION 
!               -------------------------------------------------
!
CASE('OPEN')
!
! calculate dmq
!
#ifndef  MNH_OPENACC  
   ZDMQ = DIF2X(PSRC)
#else
   CALL DIF2X_DEVICE(ZDMQ,PSRC)
#endif
   
!$acc kernels
!
! overwrite the values on the boundary to get second order difference
! for qL and qR at the boundary
!
!  WEST BOUND
!
  IF (GWEST) THEN
   ZDMQ(IIB-1,IJS:IJN,:) = -ZDMQ(IIB,IJS:IJN,:)
  ENDIF
!
!  EAST BOUND
!
  IF (GEAST) THEN 
   ZDMQ(IIE+1,IJS:IJN,:) = -ZDMQ(IIE,IJS:IJN,:)
  ENDIF
!
! monotonize the difference followinq eq. 5 in Lin94
!
!  ZDMQ(i) = Fct[ ZDMQ(i),PSRC(i),PSRC(i-1),PSRC(i+1) ]
!
!$acc loop independent collapse(3)
  do jk = 1, iku
    do jj = ijs, ijn
      do ji = iib, iie
        ZDMQ(ji, jj, jk ) = SIGN(                                                                  &
            MIN( ABS(ZDMQ(ji, jj, jk )),                                                           &
                 2.0 * ( PSRC(ji, jj, jk )                                                         &
                         - MIN(PSRC(ji - 1, jj, jk ),PSRC(ji, jj, jk ),PSRC(ji + 1, jj, jk )) ),   &
                 2.0 * (-PSRC(ji, jj, jk )                                                         &
                        + MAX(PSRC(ji - 1, jj, jk ),PSRC(ji, jj, jk ),PSRC(ji + 1, jj, jk )) )  ), &
            ZDMQ(ji, jj, jk ) )
      end do
    end do
  end do
!
!  WEST BOUND
!
!!$   ZDMQ(IIB-1,:,:) = & 
!!$        SIGN( (MIN( ABS(ZDMQ(IIB-1,:,:)), 2.0*(PSRC(IIB-1,:,:) - &
!!$        MIN(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:))),   &
!!$        2.0*(MAX(PSRC(IIE-1,:,:),PSRC(IIB-1,:,:),PSRC(IIB,:,:)) - &
!!$        PSRC(IIB-1,:,:)) )), ZDMQ(IIB-1,:,:) )
!
!  EAST BOUND
!
!!$   ZDMQ(IIE+1,:,:) = &
!!$        SIGN( (MIN( ABS(ZDMQ(IIE+1,:,:)), 2.0*(PSRC(IIE+1,:,:) - &
!!$        MIN(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:))),  &
!!$        2.0*(MAX(PSRC(IIE,:,:),PSRC(IIE+1,:,:),PSRC(IIB+1,:,:)) - &
!!$        PSRC(IIE+1,:,:)) )), ZDMQ(IIE+1,:,:) )
!
!
!  update ZDMQ HALO before next/further  utilisation 
!
#ifndef MNH_OPENACC
   CALL  GET_HALO(ZDMQ, HNAME='ZDMQ')
#else
!$acc end kernels
   CALL  GET_HALO_D(ZDMQ, HDIR="01_X", HNAME='ZDMQ')
#endif
!$acc kernels
!
! calculate qL and qR
!
!  ZQL0(i) = Fct[ PSRC(i),PSRC(i-1),ZDMQ(i),ZDMQ(i-1) ]
!
!$acc loop independent collapse(3)
  do jk = 1, iku
    do jj = ijs, ijn
      do ji = iib, iie + 1
        ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji-1, jj, jk ) ) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji-1, jj, jk ) ) / 6.0
      end do
    end do
  end do
!
#ifndef MNH_OPENACC
   CALL  GET_HALO(ZQL0, HNAME='ZQL0')
#else
!$acc end kernels
   CALL  GET_HALO_D(ZQL0,HDIR="01_X", HNAME='ZQL0')
!$acc kernels
#endif
!  
!  WEST BOUND
!
  IF (GWEST) THEN
   ZQL0(IIB-1,IJS:IJN,:) = ZQL0(IIB,IJS:IJN,:)
  ENDIF
!
   ZQR0(IIB-1:IIE,IJS:IJN,:) = ZQL0(IIB:IIE+1,IJS:IJN,:)
!
#ifndef MNH_OPENACC
   CALL  GET_HALO(ZQR0, HNAME='ZQR0')
#else
!$acc end kernels
   CALL  GET_HALO_D(ZQR0, HDIR="01_X", HNAME='ZQR0')
!$acc kernels
#endif
!
!  EAST BOUND
!
  IF (GEAST) THEN
   ZQR0(IIE+1,IJS:IJN,:) = ZQR0(IIE,IJS:IJN,:)
  ENDIF
#ifndef MNH_OPENACC
!
! determine initial coefficients of the parabolae
!
   ZDQ(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:) - ZQL0(:,IJS:IJN,:)
   ZQ60(:,IJS:IJN,:) = 6.0*(PSRC(:,IJS:IJN,:) - 0.5*(ZQL0(:,IJS:IJN,:) + ZQR0(:,IJS:IJN,:)))
!
! initialize final parabolae parameters
!
   ZQL(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:)
   ZQR(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:)
   ZQ6(:,IJS:IJN,:) = ZQ60(:,IJS:IJN,:)
!
! eliminate over and undershoots and create qL and qR as in Lin96
!
   WHERE ( ZDMQ(:,IJS:IJN,:) == 0.0 )
      ZQL(:,IJS:IJN,:) = PSRC(:,IJS:IJN,:)
      ZQR(:,IJS:IJN,:) = PSRC(:,IJS:IJN,:)
      ZQ6(:,IJS:IJN,:) = 0.0
   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) < -ZDQ(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) )
      ZQ6(:,IJS:IJN,:) = 3.0*(ZQL0(:,IJS:IJN,:) - PSRC(:,IJS:IJN,:))
      ZQR(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:) - ZQ6(:,IJS:IJN,:)
      ZQL(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:)
   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) > ZDQ(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) )
      ZQ6(:,IJS:IJN,:) = 3.0*(ZQR0(:,IJS:IJN,:) - PSRC(:,IJS:IJN,:))
      ZQL(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:) - ZQ6(:,IJS:IJN,:)
      ZQR(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:)
   END WHERE
!
! recalculate coefficients of the parabolae
!
   ZDQ(:,IJS:IJN,:) = ZQR(:,IJS:IJN,:) - ZQL(:,IJS:IJN,:)
#else
!$acc loop independent collapse(3)
DO K=1,IKU 
 DO J = IJS,IJN
    ! acc loop vector(24)
   DO I=1,IIU
!
! determine initial coefficients of the parabolae
!
   ZDQ (I,J,K)= ZQR0(I,J,K) - ZQL0(I,J,K)
   ZQ60(I,J,K) = 6.0*(PSRC(I,J,K) - 0.5*(ZQL0(I,J,K) + ZQR0(I,J,K)))
!
! initialize final parabolae parameters
!
   ZQL(I,J,K) = ZQL0(I,J,K)
   ZQR(I,J,K) = ZQR0(I,J,K)
   ZQ6(I,J,K) = ZQ60(I,J,K)
!
! eliminate over and undershoots and create qL and qR as in Lin96
!
   IF ( ZDMQ(I,J,K) == 0.0 ) THEN
      ZQL(I,J,K) = PSRC(I,J,K)
      ZQR(I,J,K) = PSRC(I,J,K)
      ZQ6(I,J,K) = 0.0
   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
      ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K))
      ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K)
      ZQL(I,J,K) = ZQL0(I,J,K)
   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
      ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K))
      ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K)
      ZQR(I,J,K) = ZQR0(I,J,K)
   ENDIF
!
! recalculate coefficients of the parabolae
!
   ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
ENDDO ; ENDDO ; ENDDO
#endif
!
! and finally calculate fluxes for the advection
!
!
!  ZFPOS(i) = Fct[ ZQR(i-1),PCR(i),ZDQ(i-1),ZQ6(i-1) ]
!
!!$   ZFPOS(IIB+1:IIE+1,:,:) = ZQR(IIB:IIE,:,:) - 0.5*PCR(IIB+1:IIE+1,:,:) * &            
!!$        (ZDQ(IIB:IIE,:,:) - (1.0 - 2.0*PCR(IIB+1:IIE+1,:,:)/3.0)          &
!!$        * ZQ6(IIB:IIE,:,:))
!$acc loop independent collapse(3)
  do jk = 1, iku
    do jj = ijs, ijn
      do ji = iib, iie + 1
        ZFPOS(ji, jj, jk ) = ZQR(ji - 1, jj, jk ) - 0.5 * PCR(ji, jj, jk )                  &
                            * ( ZDQ(ji - 1, jj, jk) - (1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0 ) &
                            * ZQ6(ji - 1, jj, jk) )
      end do
    end do
  end do
!
!
#ifndef MNH_OPENACC
   CALL GET_HALO(ZFPOS, HNAME='ZFPOS')
#else
!$acc end kernels
   CALL GET_HALO_D(ZFPOS, HDIR="01_X", HNAME='ZFPOS')
!$acc kernels
#endif
!
!
!  WEST BOUND
!
! advection flux at open boundary when u(IIB) > 0
! 
  IF (GWEST) THEN
   ZFPOS(IIB,IJS:IJN,:) = (PSRC(IIB-1,IJS:IJN,:) - ZQR(IIB-1,IJS:IJN,:))*PCR(IIB,IJS:IJN,:) + &
                    ZQR(IIB-1,IJS:IJN,:)
! PPOSX(IIB-1,:,:) is not important for the calc of advection so 
! we set it to 0
!!$   ZFPOS(IIB-1,:,:) = 0.0
   ENDIF
!
!!$   ZFNEG(IIB-1:IIE,:,:) = ZQL(IIB-1:IIE,:,:) - 0.5*PCR(IIB-1:IIE,:,:) *  &            
!!$        (ZDQ(IIB-1:IIE,:,:) + (1.0 + 2.0*PCR(IIB-1:IIE,:,:)/3.0)         &
!!$        * ZQ6(IIB-1:IIE,:,:))
!$acc loop independent collapse(3)
  do jk = 1, iku
    do jj = ijs, ijn
      do ji = 1, iiu
   ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5*PCR(ji, jj, jk ) *      &
        ( ZDQ(ji, jj, jk ) + (1.0 + 2.0*PCR(ji, jj, jk )/3.0) * ZQ6(ji, jj, jk ) )
      end do
    end do
  end do
!
#ifndef MNH_OPENACC
   CALL GET_HALO(ZFNEG, HNAME='ZFNEG')
#else
!$acc end kernels
   CALL GET_HALO_D(ZFNEG, HDIR="01_X", HNAME='ZFNEG')
!$acc kernels
#endif
!
!  EAST BOUND
!
! advection flux at open boundary when u(IIE+1) < 0
  IF (GEAST) THEN
   ZFNEG(IIE+1,IJS:IJN,:) = (ZQR(IIE,IJS:IJN,:)-PSRC(IIE+1,IJS:IJN,:))*PCR(IIE+1,IJS:IJN,:) + &
                      ZQR(IIE,IJS:IJN,:)
  ENDIF
!
! advect the actual field in X direction by U*dt
!
#ifndef MNH_OPENACC
   PR = DXF( PCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                             ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
   CALL GET_HALO(PR, HNAME='PR')
#else
   !mxm(ZQL,PRHO)
!$acc end kernels
   CALL MXM_DEVICE(PRHO,ZQL)
!$acc kernels
   where ( PCR(:,:,:) > 0. )
     ZQR(:,:,:) = PCR(:,:,:) * ZQL(:,:,:) * ZFPOS(:,:,:)
   elsewhere
     ZQR(:,:,:) = PCR(:,:,:) * ZQL(:,:,:) * ZFNEG(:,:,:)
   end where
    !dxf(PR,ZQR)
!$acc end kernels
   CALL DXF_DEVICE(ZQR,PR)
   CALL GET_HALO_D(PR, HDIR="01_X", HNAME='PR')
#endif
!
END SELECT
!
#ifdef MNH_BITREP_OMP
CALL SBR_FZ(PR(:,:,:))
#endif
!
IF (MPPDB_INITIALIZED) THEN
  !Check all INOUT arrays
  CALL MPPDB_CHECK(PSRC,"PPM_01_X end:PSRC")
  !Check all OUT arrays
  CALL MPPDB_CHECK(PR,"PPM_01_X end:PR")
END IF

!$acc end data

#ifndef MNH_OPENACC
CONTAINS
#else
END SUBROUTINE PPM_01_X_D
#endif
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!     ########################################################################
      FUNCTION DIF2X(PQ) RESULT(DQ)
!     ########################################################################
!!