diff --git a/src/MNH/ppm.f90 b/src/MNH/ppm.f90
index 5774ff8dd74f7f0a8fbcbf30e60c7d0a0135614b..2085af4a1c141d2c7be514272a260e4646a2024b 100644
--- a/src/MNH/ppm.f90
+++ b/src/MNH/ppm.f90
@@ -7,6 +7,7 @@
! 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 _OPENACC
!
@@ -22,17 +23,17 @@
!#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 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 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)
+! #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
! ###############
@@ -52,11 +53,8 @@ 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
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -85,11 +83,8 @@ 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
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -97,7 +92,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -118,12 +112,9 @@ INTEGER, INTENT(IN) :: KGRID ! C grid localisation
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRC ! variable at t
#else
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC ! variable at t
-!$acc declare present(PSRC)
#endif
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -131,7 +122,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -152,11 +142,8 @@ 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
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -164,7 +151,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -186,12 +172,9 @@ 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
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -199,7 +182,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -219,11 +201,8 @@ SUBROUTINE PPM_S0_Z(KGRID, PSRC, PCR, PRHO, PTSTEP &
INTEGER, INTENT(IN) :: KGRID ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC ! variable at t
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -231,7 +210,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -253,13 +231,9 @@ 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
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHOT ! density at t+dt
-!$acc declare present(PRHOT)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -267,7 +241,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -289,13 +262,9 @@ 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
-!$acc declare present(PSRC)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHOT ! density at t+dt
-!$acc declare present(PRHOT)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -303,7 +272,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -326,14 +294,10 @@ INTEGER, INTENT(IN) :: KGRID ! C grid localisation
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRC ! variable at t
#else
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC ! variable at t
-!$acc declare present(PSRC)
#endif
REAL, DIMENSION(:,:,:), INTENT(IN) :: PCR ! Courant number
-!$acc declare present(PCR)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO ! density
-!$acc declare present(PRHO)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHOT ! density at t+dt
-!$acc declare present(PRHOT)
REAL, INTENT(IN) :: PTSTEP ! Time step
!
! output source term
@@ -341,7 +305,6 @@ REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
#else
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
-!$acc declare present(PR)
#endif
!
#ifndef _OPENACC
@@ -427,23 +390,23 @@ CONTAINS
!!
!-------------------------------------------------------------------------------
!
+USE MODD_CONF
+
USE MODE_ll
use mode_mppdb
#ifdef _OPENACC
use mode_msg
#endif
+
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
+USE MODI_GET_HALO
#ifndef _OPENACC
USE MODI_SHUMAN
#else
USE MODI_SHUMAN_DEVICE
#endif
-USE MODI_GET_HALO
-!
-USE MODD_CONF
-!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
-#ifdef MNH_BITREP
-USE MODI_BITREP
-#endif
!
IMPLICIT NONE
!
@@ -575,10 +538,15 @@ CALL GET_HALO_D(PSRC,HDIR="01_X", HNAME='PSRC')
end do
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
!
@@ -595,11 +563,7 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCX(1) == HLBCX(2)
#endif
!
! calculate dmq
-#ifndef _OPENACC
ZDMQ = DIF2X(PSRC)
-#else
- dif2x(ZDMQ,PSRC)
-#endif
!
! monotonize the difference followinq eq. 5 in Lin94
!
@@ -732,12 +696,8 @@ CASE('OPEN')
!
! calculate dmq
!
-!$acc kernels
-#ifndef _OPENACC
ZDMQ = DIF2X(PSRC)
-#else
- dif2x(ZDMQ,PSRC)
-#endif
+!$acc kernels
!
! overwrite the values on the boundary to get second order difference
! for qL and qR at the boundary
@@ -762,11 +722,11 @@ CASE('OPEN')
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 ) &
+ 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 ) &
+ 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
@@ -1017,6 +977,9 @@ END IF
#ifndef _OPENACC
CONTAINS
+#else
+END SUBROUTINE PPM_01_X_D
+#endif
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
@@ -1048,6 +1011,7 @@ IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+!$acc declare present(PQ, DQ)
!
!* 0.2 Declarations of local variables :
!
@@ -1067,18 +1031,17 @@ IJB=2 ; IJE = SIZE(PQ,2) -1
!
!* 2.0. COMPUTE THE DIFFERENCE
! ----------------------
-!
+!
+!$acc kernels
DQ(IIB:IIE,:,:) = PQ(IIB+1:IIE+1,:,:) - PQ(IIB-1:IIE-1,:,:)
DQ(IIB-1,:,:) = PQ(IIB,:,:) - PQ(IIE-1,:,:)
DQ(IIE+1,:,:) = PQ(IIB+1,:,:) - PQ(IIE,:,:)
-DQ = 0.5*DQ
+DQ = 0.5*DQ
+!$acc end kernels
!
END FUNCTION DIF2X
-#endif
!
#ifdef _OPENACC
-END SUBROUTINE PPM_01_X_D
-
END SUBROUTINE PPM_01_X
#else
END FUNCTION PPM_01_X
@@ -1157,24 +1120,23 @@ CONTAINS
!!
!-------------------------------------------------------------------------------
!
+USE MODD_CONF
+
USE MODE_ll
#ifdef _OPENACC
use mode_msg
#endif
use mode_mppdb
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
+USE MODI_GET_HALO
#ifndef _OPENACC
USE MODI_SHUMAN
#else
USE MODI_SHUMAN_DEVICE
#endif
-USE MODI_GET_HALO
-#ifdef MNH_BITREP
-USE MODI_BITREP
-#endif
-!
-USE MODD_CONF
-!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
!
IMPLICIT NONE
!
@@ -1228,12 +1190,12 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZFPOS, ZFNEG
#else
! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
REAL, DIMENSION(IIU,IJU,IKU) :: &
- ZQL,ZQR , ZDQ,ZQ6 , ZDMQ &
+ ZQL,ZQR , ZDQ,ZQ6 , ZDMQ &
! extra variables for the initial guess of parabolae parameters
, ZQL0,ZQR0,ZQ60 &
! advection fluxes
, ZFPOS, ZFNEG
-!$acc declare present (ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG)
+!$acc declare present (ZQL,ZQR,ZDQ,ZQ6,ZDMQ,ZQL0,ZQR0,ZQ60,ZFPOS,ZFNEG)
!
!JUAN ACC
@@ -1310,10 +1272,15 @@ end do
ZFPOS=PSRC
ZFNEG=PSRC
#else
+#if 0
ZFPOS(1:IIW,:,IKB:IKE)=PSRC(1:IIW,:,IKB:IKE)
ZFNEG(1:IIW,:,IKB:IKE)=PSRC(1:IIW,:,IKB:IKE)
ZFPOS(IIA:,:,IKB:IKE)=PSRC(IIA:,:,IKB:IKE)
ZFNEG(IIA:,:,IKB:IKE)=PSRC(IIA:,:,IKB:IKE)
+#else
+ZFPOS(:,:,:) = PSRC(:,:,:)
+ZFNEG(:,:,:) = PSRC(:,:,:)
+#endif
#endif
!$acc end kernels
!
@@ -1323,22 +1290,16 @@ SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
! ---------------------------------------------
!
CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
-#ifdef _OPENACC
- call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_Y', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
-#endif
!
! calculate dmq
-#ifndef _OPENACC
ZDMQ = DIF2Y(PSRC)
-#else
- dif2y(ZDMQ,PSRC)
-#endif
!
! monotonize the difference followinq eq. 5 in Lin94
!BEG JUAN PPM_LL01
!
! ZDMQ(j) = Fct[ ZDMQ(j),PSRC(j),PSRC(j-1),PSRC(j+1) ]
!
+!$acc kernels
ZDMQ(IIW:IIA,IJB:IJE,:) = &
SIGN( (MIN( ABS(ZDMQ(IIW:IIA,IJB:IJE,:)),2.0*(PSRC(IIW:IIA,IJB:IJE,:) - &
MIN(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:))), &
@@ -1363,26 +1324,53 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
!
! update ZDMQ HALO before next/further utilisation
!
+#ifndef _OPENACC
CALL GET_HALO(ZDMQ, HNAME='ZDMQ')
+#else
+!$acc end kernels
+ CALL GET_HALO_D(ZDMQ,HDIR="01_Y", HNAME='ZDMQ')
+!$acc kernels
+#endif
!
! calculate qL and qR with the modified dmq
!
- ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
- (ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
+! ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
+! (ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
+!$acc loop independent collapse(3)
+ do jk = 1, iku
+ do jj = ijb, ije + 1
+ do ji = iiw, iia
+ ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj-1, jk )) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj-1, jk ) ) / 6.0
+ end do
+ end do
+ end do
+!
+#ifndef _OPENACC
CALL GET_HALO(ZQL0, HNAME='ZQL0')
+#else
+!$acc end kernels
+ CALL GET_HALO_D(ZQL0,HDIR="01_Y", HNAME='ZQL0')
+!$acc kernels
+#endif
!
-!
! SOUTH BOUND
!
!!$ ZQL0(:,IJB-1,:) = ZQL0(:,IJE,:) JUAN PPMLL01
!
ZQR0(IIW:IIA,IJB-1:IJE,:) = ZQL0(IIW:IIA,IJB:IJE+1,:)
!
+#ifndef _OPENACC
CALL GET_HALO(ZQR0, HNAME='ZQR0')
+#else
+!$acc end kernels
+ CALL GET_HALO_D(ZQR0,HDIR="01_Y", HNAME='ZQR0')
+!$acc kernels
+#endif
!
! NORTH BOUND
!
!!$ ZQR0(:,IJE+1,:) = ZQR0(:,IJB,:) JUAN PPMLL01
+#ifndef _OPENACC
!
! determine initial coefficients of the parabolae
!
@@ -1414,6 +1402,47 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
! recalculate coefficients of the parabolae
!
ZDQ(IIW:IIA,:,:) = ZQR(IIW:IIA,:,:) - ZQL(IIW:IIA,:,:)
+#else
+!$acc loop independent collapse(3)
+ DO K=IKB,IKE
+ DO J=IJS,IJN
+ 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
+ ELSE IF ( 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)
+ ELSE IF ( 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)
+ END IF
+ !
+ ! recalculate coefficients of the parabolae
+ !
+ ZDQ(I,J,K) = ZQR(I,J,K) - ZQL(I,J,K)
+ !
+ END DO
+ END DO
+ END DO
+#endif
!
! and finally calculate fluxes for the advection
!
@@ -1423,7 +1452,13 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
(ZDQ(IIW:IIA,IJB-1:IJE,:) - (1.0 - 2.0*PCR(IIW:IIA,IJB:IJE+1,:)/3.0) &
* ZQ6(IIW:IIA,IJB-1:IJE,:))
!
+#ifndef _OPENACC
CALL GET_HALO(ZFPOS, HNAME='ZFPOS')
+#else
+!$acc end kernels
+ CALL GET_HALO_D(ZFPOS,HDIR="01_Y", HNAME='ZFPOS')
+!$acc kernels
+#endif
!
! SOUTH BOUND
!
@@ -1434,7 +1469,13 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
ZFNEG(IIW:IIA,:,:) = ZQL(IIW:IIA,:,:) - 0.5*PCR(IIW:IIA,:,:) * &
( ZDQ(IIW:IIA,:,:) + (1.0 + 2.0*PCR(IIW:IIA,:,:)/3.0) * ZQ6(IIW:IIA,:,:) )
!
+#ifndef _OPENACC
CALL GET_HALO(ZFNEG, HNAME='ZFNEG')
+#else
+!$acc end kernels
+ CALL GET_HALO_D(ZFNEG,HDIR="01_Y", HNAME='ZFNEG')
+!$acc kernels
+#endif
!
! advect the actual field in Y direction by V*dt
!
@@ -1442,8 +1483,28 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
PR = DYF( PCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + &
ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
#else
- call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_Y', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
+!$acc end kernels
+ CALL MYM_DEVICE(PRHO,ZQL)
+!$acc kernels
+!$acc loop independent collapse(3)
+ DO K=IKB,IKE
+ DO J=IJS,IJN
+ DO I=1,IIU
+ if ( PCR(I,J,K) > 0. ) then
+ ZQR(I,J,K) = PCR(I,J,K)* ZQL(I,J,K) * ZFPOS(I,J,K)
+ else
+ ZQR(I,J,K) = PCR(I,J,K)* ZQL(I,J,K) * ZFNEG(I,J,K)
+ end if
+ END DO
+ END DO
+ END DO
+!$acc end kernels
+ CALL DYF_DEVICE(ZQR,PR)
+#endif
+#ifndef _OPENACC
CALL GET_HALO(PR, HNAME='PR')
+#else
+ CALL GET_HALO_D(PR,HDIR="01_Y", HNAME='PR')
#endif
!
!* 2.2 NON-CYCLIC BOUNDARY CONDITIONS IN THE Y DIRECTION
@@ -1451,13 +1512,9 @@ CASE ('CYCL','WALL') ! In that case one must have HLBCY(1) == HLBCY(2)
!
CASE('OPEN')
!
-#ifndef _OPENACC
! calculate dmq
ZDMQ = DIF2Y(PSRC)
-#else
!$acc kernels
- dif2y(ZDMQ,PSRC)
-#endif
! overwrite the values on the boundary to get second order difference
! for qL and qR at the boundary
!
@@ -1732,6 +1789,9 @@ END IF
#ifndef _OPENACC
CONTAINS
+#else
+END SUBROUTINE PPM_01_Y_D
+#endif
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
@@ -1765,6 +1825,7 @@ IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+!$acc declare present(PQ, DQ)
!
!* 0.2 Declarations of local variables :
!
@@ -1785,16 +1846,18 @@ IJB=2 ; IJE = SIZE(PQ,2) -1
!* 2.0. COMPUTE THE DIFFERENCE
! ----------------------
!
+!$acc kernels
DQ(:,IJB:IJE,:) = PQ(:,IJB+1:IJE+1,:) - PQ(:,IJB-1:IJE-1,:)
DQ(:,IJB-1,:) = PQ(:,IJB,:) - PQ(:,IJE-1,:)
DQ(:,IJE+1,:) = PQ(:,IJB+1,:) - PQ(:,IJE,:)
-DQ = 0.5 * DQ
+DQ = 0.5 * DQ
+!$acc end kernels
!
END FUNCTION DIF2Y
-#endif
+! #endif
!
#ifdef _OPENACC
-END SUBROUTINE PPM_01_Y_D
+! END SUBROUTINE PPM_01_Y_D
END SUBROUTINE PPM_01_Y
#else
@@ -1867,6 +1930,7 @@ CONTAINS
!-------------------------------------------------------------------------------
!
USE MODE_ll
+
#ifndef _OPENACC
USE MODI_SHUMAN
#else
@@ -1967,6 +2031,39 @@ iiu = size( PSRC, 1 )
iju = size( PSRC, 2 )
iku = size( PSRC, 3 )
#endif
+
+!$acc kernels
+!$acc loop independent collapse(3)
+ do jk = 1, iku
+ do jj = 1, iju
+ do ji = 1, iiu
+ 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 )
+ end do
+ end do
+end do
+#ifndef _OPENACC
+ZFPOS=PSRC
+ZFNEG=PSRC
+#else
+#if 0
+ZFPOS(1:IIW,:,IKB:IKE)=PSRC(1:IIW,:,IKB:IKE)
+ZFNEG(1:IIW,:,IKB:IKE)=PSRC(1:IIW,:,IKB:IKE)
+ZFPOS(IIA:,:,IKB:IKE)=PSRC(IIA:,:,IKB:IKE)
+ZFNEG(IIA:,:,IKB:IKE)=PSRC(IIA:,:,IKB:IKE)
+#else
+ZFPOS(:,:,:) = PSRC(:,:,:)
+ZFNEG(:,:,:) = PSRC(:,:,:)
+#endif
+#endif
+!$acc end kernels
!
!-------------------------------------------------------------------------------
!
@@ -1974,12 +2071,8 @@ iku = size( PSRC, 3 )
! --------------------------------
!
! calculate dmq
-#ifndef _OPENACC
ZDMQ = DIF2Z(PSRC)
-#else
!$acc kernels
- dif2z(ZDMQ,PSRC)
-#endif
!
! monotonize the difference followinq eq. 5 in Lin94
! use the periodic BC here, it doesn't matter for vertical (hopefully)
@@ -2028,11 +2121,14 @@ do jk = ikb - 1, ike
end do
end do
ZQR0(:,:,IKE+1) = ZQR0(:,:,IKB)
+#ifndef _OPENACC
!
! determine initial coefficients of the parabolae
!
-!$acc loop independent collapse(3)
-do jk = ikb - 1, ike
+!Note: do loop on jk is done from 1 to iku to prevent problems with unitialized value
+! in the next "where"
+! do jk = ikb - 1, ike
+do jk = 1, iku
do jj = 1, iju
do ji = 1, iiu
ZDQ (ji, jj, jk ) = ZQR0(ji, jj, jk ) - ZQL0(ji, jj, jk )
@@ -2040,7 +2136,6 @@ do jk = ikb - 1, ike
end do
end do
end do
-#ifndef _OPENACC
!
! initialize final parabolae parameters
!
@@ -2233,6 +2328,9 @@ END IF
!
#ifndef _OPENACC
CONTAINS
+#else
+END SUBROUTINE PPM_01_Z_D
+#endif
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
@@ -2266,6 +2364,7 @@ IMPLICIT NONE
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PQ
REAL, DIMENSION(SIZE(PQ,1),SIZE(PQ,2),SIZE(PQ,3)) :: DQ
+!$acc declare present(PQ, DQ)
!
!* 0.2 Declarations of local variables :
!
@@ -2285,17 +2384,19 @@ IKE = SIZE(PQ,3) - JPVEXT
!
!* 2.0. COMPUTE THE DIFFERENCE
! ----------------------
-!
+!
+!$acc kernels
DQ(:,:,IKB:IKE) = PQ(:,:,IKB+1:IKE+1) - PQ(:,:,IKB-1:IKE-1)
DQ(:,:,IKB-1) = -DQ(:,:,IKB)
DQ(:,:,IKE+1) = -DQ(:,:,IKE)
DQ = 0.5 * DQ
+!$acc end kernels
!
END FUNCTION DIF2Z
-#endif
+! #endif
!
#ifdef _OPENACC
-END SUBROUTINE PPM_01_Z_D
+! END SUBROUTINE PPM_01_Z_D
!
END SUBROUTINE PPM_01_Z
#else
@@ -2371,21 +2472,21 @@ CONTAINS
!!
!-------------------------------------------------------------------------------
!
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODD_CONF
+
USE MODE_ll
#ifdef _OPENACC
use mode_msg
#endif
+
+USE MODI_GET_HALO
#ifndef _OPENACC
USE MODI_SHUMAN
#else
USE MODI_SHUMAN_DEVICE
#endif
-USE MODI_GET_HALO
!
-USE MODD_CONF
-!BEG JUAN PPM_LL
-USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
-!END JUAN PPM_LL
#ifdef _OPENACC
USE MODD_PARAMETERS, ONLY : JPHEXT
!
@@ -2718,8 +2819,6 @@ CALL GET_HALO(PR, HNAME='PR')
#else
CALL GET_HALO_D(PR, HDIR="S0_X", HNAME='PR')
#endif
-CALL MPPDB_CHECK3DM("PPM::PPM_S0_X OPEN ::PR",PRECISION,PR)
-!
!-------------------------------------------------------------------------------
CALL DEL_HALO2_ll(TZ_PSRC_HALO2_ll)
!
@@ -2805,19 +2904,20 @@ CONTAINS
!!
!-------------------------------------------------------------------------------
!
+USE MODD_CONF
+
USE MODE_ll
#ifdef _OPENACC
use mode_msg
#endif
+
+USE MODI_GET_HALO
#ifndef _OPENACC
USE MODI_SHUMAN
#else
USE MODI_SHUMAN_DEVICE
#endif
-USE MODI_GET_HALO
!
-USE MODD_CONF
-!USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
#ifdef _OPENACC
USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
USE MODD_PARAMETERS, ONLY : JPHEXT
@@ -2911,8 +3011,13 @@ GSOUTH = LSOUTH_ll()
!-------------------------------------------------------------------------------
!
IF ( L2D ) THEN
- PR = PSRC*PRHO
- RETURN
+!$acc kernels
+ PR(:, :, : ) = PSRC(:, :, : ) * PRHO(:, :, : )
+!$acc end kernels
+
+ !Check all OUT arrays
+ CALL MPPDB_CHECK(PR,"PPM_S0_Y end:PR")
+ RETURN
END IF
!
CALL GET_HALO2(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC')
@@ -3151,7 +3256,6 @@ CALL GET_HALO(PR, HNAME='PR')
#else
CALL GET_HALO_D(PR, HDIR="S0_Y", HNAME='PR')
#endif
-CALL MPPDB_CHECK3DM("PPM::PPM_S0_Y OPEN ::PR",PRECISION,PR)
!
CALL DEL_HALO2_ll(TZ_PSRC_HALO2_ll)
!
@@ -3397,11 +3501,6 @@ PR = PSRC * PRHO - &
#ifndef _OPENACC
CALL GET_HALO(PR, HNAME='PR')
#else
-!PW: BUG?: it is necessary to update self before call to get_halo_d
-!and update device after to get correct results on 4 processes and more with OpenACC
-!The reason is not clear, maybe a (compiler PGI 16.4?) bug in the updates of get_halo_d
-!$acc update self(PR)
-!$acc update device(PR)
CALL GET_HALO_D(PR, HNAME='PR')
#endif
IF (MPPDB_INITIALIZED) THEN
@@ -4164,6 +4263,7 @@ END FUNCTION PPM_S1_Y
USE MODE_ll
USE MODE_IO
use mode_msg
+
#ifndef _OPENACC
USE MODI_SHUMAN
#else