diff --git a/src/ZSOLVER/gdiv.f90 b/src/ZSOLVER/gdiv.f90
index b22065908c6e897d61706999fd1374582060f1d4..d548d9fa49aca7695193f55b0d6d2022574124b8 100644
--- a/src/ZSOLVER/gdiv.f90
+++ b/src/ZSOLVER/gdiv.f90
@@ -106,6 +106,10 @@ USE MODI_CONTRAV
!
USE MODE_ll
!
+#ifdef MNH_OPENACC
+USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D
+#endif
+!
IMPLICIT NONE
!
!* 0.1 declarations of arguments
@@ -145,6 +149,14 @@ INTEGER :: IKE ! indice K for the last inner mass point along z
!
INTEGER :: JI,JJ,JK ! loop indexes
!
+#ifdef MNH_OPENACC
+INTEGER :: IIU,IJU,IKU
+!
+REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZTMP1,ZTMP2
+INTEGER :: IZTMP1,IZTMP2
+#endif
+!
+LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH
!
!-------------------------------------------------------------------------------
!
@@ -156,6 +168,20 @@ CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IKB=1+JPVEXT
IKE=SIZE(PU,3) - JPVEXT
!
+#ifdef MNH_OPENACC
+IIU=SIZE(PU,1)
+IJU=SIZE(PU,2)
+IKU=SIZE(PU,3)
+!
+IZTMP1 = MNH_ALLOCATE_ZT3D( ZTMP1,IIU,IJU,IKU )
+IZTMP2 = MNH_ALLOCATE_ZT3D( ZTMP2,IIU,IJU,IKU )
+#endif
+!
+GWEST = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() )
+GEAST = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() )
+GSOUTH = ( .NOT. L2D .AND. HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() )
+GNORTH = ( .NOT. L2D .AND. HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() )
+!
!-------------------------------------------------------------------------------
!
!* 2. COMPUTE THE CONTRAVARIANT COMPONENTS
@@ -164,22 +190,31 @@ IKE=SIZE(PU,3) - JPVEXT
!* 2.1 prepare the boundary conditions
!
!
- PU(:,:,IKB-1)=PU(:,:,IKB)
- PU(:,:,IKE+1)=PU(:,:,IKE)
- PV(:,:,IKB-1)=PV(:,:,IKB)
- PV(:,:,IKE+1)=PV(:,:,IKE)
-
+!$acc kernels
+ DO CONCURRENT ( JI=1:IIU,JJ=1:IJU )
+ PU(JI,JJ,IKB-1)=PU(JI,JJ,IKB)
+ PU(JI,JJ,IKE+1)=PU(JI,JJ,IKE)
+ PV(JI,JJ,IKB-1)=PV(JI,JJ,IKB)
+ PV(JI,JJ,IKE+1)=PV(JI,JJ,IKE)
+ END DO
+!$acc end kernels
!
!
!* 2.1 compute the contravariant components
!
+#ifndef MNH_OPENACC
CALL CONTRAV(HLBCX,HLBCY,PU,PV,PW,PDXX,PDYY,PDZZ,PDZX,PDZY,ZUC,ZVC,ZWC,4)
+#else
+CALL CONTRAV_DEVICE(HLBCX,HLBCY,PU,PV,PW,PDXX,PDYY,PDZZ,PDZX,PDZY,ZUC,ZVC,ZWC,4, &
+ ZTMP1,ZTMP2,ODATA_ON_DEVICE=.TRUE.)
+#endif
!
!-------------------------------------------------------------------------------
!
!* 3. COMPUTE THE DIVERGENCE
! ----------------------
!
+!$acc kernels
PGDIV=0. !usefull for the four corners and halo zones
!
Z1(IIB:IIE,:,:)=ZUC(IIB+1:IIE+1,:,:)-ZUC(IIB:IIE,:,:)
@@ -190,7 +225,8 @@ PGDIV(IIB:IIE,IJB:IJE,IKB:IKE)= Z1(IIB:IIE,IJB:IJE,IKB:IKE) + &
Z2(IIB:IIE,IJB:IJE,IKB:IKE) + &
Z3(IIB:IIE,IJB:IJE,IKB:IKE)
! only the divergences computed
- ! in the inner mass points are meaningful
+ ! in the inner mass points are meaningful
+!$acc end kernels
!
!-------------------------------------------------------------------------------
!
@@ -201,12 +237,14 @@ PGDIV(IIB:IIE,IJB:IJE,IKB:IKE)= Z1(IIB:IIE,IJB:IJE,IKB:IKE) + &
!
! we set the divergence equal to the vertical contravariant component above
! and under the physical domain
+!$acc kernels async
DO JJ=IJB,IJE
DO JI=IIB,IIE
PGDIV(JI,JJ,IKB-1)=ZWC(JI,JJ,IKB)
PGDIV(JI,JJ,IKE+1)=ZWC(JI,JJ,IKE+1)
END DO
END DO
+!$acc end kernels
!
!* 4.2 set divergence at the lateral boundaries
!
@@ -214,64 +252,89 @@ END DO
! the right and the left of the physical domain in both horizontal directions
! for non-periodic cases
!
-IF(HLBCX(1) /= 'CYCL' .AND. LWEST_ll()) THEN
- DO JK=IKB,IKE
- DO JJ=IJB,IJE
- PGDIV(IIB-1,JJ,JK)=ZUC(IIB,JJ,JK)
- END DO
- END DO
+IF( GWEST ) THEN
+ !$acc kernels async
+ DO JK=IKB,IKE
+ DO JJ=IJB,IJE
+ PGDIV(IIB-1,JJ,JK)=ZUC(IIB,JJ,JK)
+ END DO
+ END DO
+ !$acc end kernels
END IF
!
-IF(HLBCX(2) /= 'CYCL' .AND. LEAST_ll()) THEN
- DO JK=IKB,IKE
- DO JJ=IJB,IJE
- PGDIV(IIE+1,JJ,JK)=ZUC(IIE+1,JJ,JK)
- END DO
- END DO
+IF( GEAST ) THEN
+ !$acc kernels async
+ DO JK=IKB,IKE
+ DO JJ=IJB,IJE
+ PGDIV(IIE+1,JJ,JK)=ZUC(IIE+1,JJ,JK)
+ END DO
+ END DO
+ !$acc end kernels
END IF
!
!
-IF (.NOT. L2D .AND. HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll()) THEN
- DO JK=IKB,IKE
- DO JI=IIB,IIE
- PGDIV(JI,IJB-1,JK)=ZVC(JI,IJB,JK)
- END DO
- END DO
+IF ( GSOUTH ) THEN
+ !$acc kernels async
+ DO JK=IKB,IKE
+ DO JI=IIB,IIE
+ PGDIV(JI,IJB-1,JK)=ZVC(JI,IJB,JK)
+ END DO
+ END DO
+ !$acc end kernels
END IF
!
-IF (.NOT. L2D .AND. HLBCY(2) /= 'CYCL' .AND. LNORTH_ll()) THEN
- DO JK=IKB,IKE
- DO JI=IIB,IIE
- PGDIV(JI,IJE+1,JK)=ZVC(JI,IJE+1,JK)
- END DO
- END DO
+IF ( GNORTH ) THEN
+ !$acc kernels async
+ DO JK=IKB,IKE
+ DO JI=IIB,IIE
+ PGDIV(JI,IJE+1,JK)=ZVC(JI,IJE+1,JK)
+ END DO
+ END DO
+ !$acc end kernels
END IF
!
+! wait on GPU for all boundary condition update
+!$acc wait
+!
!* 4.3 set divergence at the corner points
!
! it is the following of the condition of copy the horizontal component
! under the bottom of the model
!
-IF(HLBCX(1) /= 'CYCL' .AND. LWEST_ll()) THEN
- PGDIV(IIB-1,IJB:IJE,IKB-1)=PGDIV(IIB-1,IJB:IJE,IKB)
- PGDIV(IIB-1,IJB:IJE,IKE+1)=PGDIV(IIB-1,IJB:IJE,IKE)
+IF( GWEST ) THEN
+ !$acc kernels async
+ PGDIV(IIB-1,IJB:IJE,IKB-1)=PGDIV(IIB-1,IJB:IJE,IKB)
+ PGDIV(IIB-1,IJB:IJE,IKE+1)=PGDIV(IIB-1,IJB:IJE,IKE)
+ !$acc end kernels
END IF
!
-IF (HLBCX(2) /= 'CYCL' .AND. LEAST_ll()) THEN
- PGDIV(IIE+1,IJB:IJE,IKB-1)=PGDIV(IIE+1,IJB:IJE,IKB)
- PGDIV(IIE+1,IJB:IJE,IKE+1)=PGDIV(IIE+1,IJB:IJE,IKE)
+IF ( GEAST ) THEN
+ !$acc kernels async
+ PGDIV(IIE+1,IJB:IJE,IKB-1)=PGDIV(IIE+1,IJB:IJE,IKB)
+ PGDIV(IIE+1,IJB:IJE,IKE+1)=PGDIV(IIE+1,IJB:IJE,IKE)
+ !$acc end kernels
END IF
!
-IF (.NOT. L2D .AND. HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll()) THEN
- PGDIV(IIB:IIE,IJB-1,IKB-1)=PGDIV(IIB:IIE,IJB-1,IKB)
- PGDIV(IIB:IIE,IJB-1,IKE+1)=PGDIV(IIB:IIE,IJB-1,IKE)
+IF ( GSOUTH ) THEN
+ !$acc kernels async
+ PGDIV(IIB:IIE,IJB-1,IKB-1)=PGDIV(IIB:IIE,IJB-1,IKB)
+ PGDIV(IIB:IIE,IJB-1,IKE+1)=PGDIV(IIB:IIE,IJB-1,IKE)
+ !$acc end kernels
END IF
!
-IF (.NOT. L2D .AND. HLBCY(2) /= 'CYCL' .AND. LNORTH_ll()) THEN
- PGDIV(IIB:IIE,IJE+1,IKB-1)=PGDIV(IIB:IIE,IJE+1,IKB)
- PGDIV(IIB:IIE,IJE+1,IKE+1)=PGDIV(IIB:IIE,IJE+1,IKE)
+IF ( GNORTH ) THEN
+ !$acc kernels async
+ PGDIV(IIB:IIE,IJE+1,IKB-1)=PGDIV(IIB:IIE,IJE+1,IKB)
+ PGDIV(IIB:IIE,IJE+1,IKE+1)=PGDIV(IIB:IIE,IJE+1,IKE)
+ !$acc end kernels
END IF
!
+! wait on GPU for all corner update
+!$acc wait
+!
+#ifdef MNH_OPENACC
+CALL MNH_REL_ZT3D ( IZTMP1,IZTMP2 )
+#endif
!
!-------------------------------------------------------------------------------
!
diff --git a/src/ZSOLVER/pressurez.f90 b/src/ZSOLVER/pressurez.f90
index c9b2b9db971f01c05b40f5c809a24c8f18a47d18..06210d56f144c701f472cd212e5e2213779c0701 100644
--- a/src/ZSOLVER/pressurez.f90
+++ b/src/ZSOLVER/pressurez.f90
@@ -267,6 +267,15 @@ USE MODI_MASS_LEAK
USE MODI_P_ABS
USE MODI_RICHARDSON
USE MODI_SHUMAN
+#ifdef MNH_OPENACC
+USE MODI_SHUMAN_DEVICE
+USE MODI_GET_HALO
+#endif
+!
+#ifdef MNH_BITREP
+USE MODI_BITREP
+USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D
+#endif
!
IMPLICIT NONE
!
@@ -393,6 +402,14 @@ TYPE(LIST_ll), POINTER :: TZFIELDS_ll, TZFIELDS2_ll ! list of fields to exchang
INTEGER :: IIB_I,IIE_I,IJB_I,IJE_I
INTEGER :: IIMAX_ll,IJMAX_ll
!
+REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZPRHODJ,ZMXM_PRHODJ,ZMZM_PRHODJ,ZGZ_M_W
+INTEGER :: IZPRHODJ,IZMXM_PRHODJ,IZMZM_PRHODJ,IZGZ_M_W
+REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZMYM_PRHODJ
+INTEGER :: IZMYM_PRHODJ
+!
+!
+LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH
+LOGICAL :: GSOUTH2D,GNORTH2D,GPRVREF0
!
!------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
@@ -416,10 +433,27 @@ IKB= 1+JPVEXT
IKU= SIZE(PPABST,3)
IKE= IKU - JPVEXT
!
+GWEST = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() )
+GEAST = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() )
+GSOUTH = ( .NOT. L2D .AND. HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() )
+GNORTH = ( .NOT. L2D .AND. HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() )
+GSOUTH2D = ( L2D .AND. LSOUTH_ll() )
+GNORTH2D = ( L2D .AND. LNORTH_ll() )
+!
+GPRVREF0 = ( SIZE(PRVREF,1) == 0 )
+!
+IZPRHODJ = MNH_ALLOCATE_ZT3D( ZPRHODJ,IIU,IJU,IKU )
+IZMXM_PRHODJ = MNH_ALLOCATE_ZT3D( ZMXM_PRHODJ,IIU,IJU,IKU )
+IZMZM_PRHODJ = MNH_ALLOCATE_ZT3D( ZMZM_PRHODJ,IIU,IJU,IKU )
+IZGZ_M_W = MNH_ALLOCATE_ZT3D( ZGZ_M_W,IIU,IJU,IKU )
+IZMYM_PRHODJ = MNH_ALLOCATE_ZT3D( ZMYM_PRHODJ,IIU,IJU,IKU )
+!
+!$acc kernels
ZPABS_S(:,:) = 0.
ZPABS_N(:,:) = 0.
ZPABS_E(:,:) = 0.
ZPABS_W(:,:) = 0.
+!$acc end kernels
!
!
!-------------------------------------------------------------------------------
@@ -454,15 +488,27 @@ CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
!
! The non-homogenous Neuman problem is transformed in an homogenous Neuman
! problem in the non-periodic cases
-IF (HLBCX(1) /= 'CYCL') THEN
- IF (LWEST_ll()) ZDV_SOURCE(IIB-1,:,:) = 0.
- IF (LEAST_ll()) ZDV_SOURCE(IIE+1,:,:) = 0.
-ENDIF
-!
-IF (.NOT. L2D .AND. HLBCY(1) /= 'CYCL') THEN
- IF (LSOUTH_ll()) ZDV_SOURCE(:,IJB-1,:) = 0.
- IF (LNORTH_ll()) ZDV_SOURCE(:,IJE+1,:) = 0.
-ENDIF
+IF ( GWEST ) THEN
+ !$acc kernels async
+ ZDV_SOURCE(IIB-1,:,:) = 0.
+ !$acc end kernels
+END IF
+IF ( GEAST ) THEN
+ !$acc kernels async
+ ZDV_SOURCE(IIE+1,:,:) = 0.
+ !$acc end kernels
+END IF
+IF ( GSOUTH ) THEN
+ !$acc kernels async
+ ZDV_SOURCE(:,IJB-1,:) = 0.
+ !$acc end kernels
+END IF
+IF ( GNORTH ) THEN
+ !$acc kernels async
+ ZDV_SOURCE(:,IJE+1,:) = 0.
+ !$acc end kernels
+END IF
+!$acc wait
!
!-------------------------------------------------------------------------------
!
@@ -474,6 +520,7 @@ ENDIF
! -------------------------------------------------------
!
IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
+ !$acc kernels
IF(KRR > 0) THEN
!
! compute the ratio : 1 + total water mass / dry air mass
@@ -489,8 +536,15 @@ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
! compute the virtual potential temperature when water is absent
ZTHETAV(:,:,:) = PTHT(:,:,:)
END IF
+ !$acc end kernels
!
+#ifndef MNH_OPENACC
ZPHIT(:,:,:)=(PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)
+#else
+ !$acc kernels
+ ZPHIT(:,:,:)=BR_POW((PPABST(:,:,:)/XP00),(XRD/XCPD))-PEXNREF(:,:,:)
+ !$acc end kernels
+#endif
!
ELSEIF(CEQNSYS=='LHE') THEN
ZPHIT(:,:,:)= ((PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) &
@@ -553,26 +607,48 @@ END IF
!* 6. ADD THE PRESSURE GRADIENT TO THE SOURCES
! ----------------------------------------
!
-IF ( HLBCX(1) /= 'CYCL' ) THEN
- IF(LWEST_ll()) ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB-1)
- IF(LEAST_ll()) ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB-1)
-ENDIF
-IF ( HLBCY(1) /= 'CYCL' ) THEN
- IF (LSOUTH_ll()) ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB-1)
- IF (LNORTH_ll()) ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB-1)
-ENDIF
-!
-IF ( L2D ) THEN
- IF (LSOUTH_ll()) ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)
- IF (LNORTH_ll()) ZPHIT(:,IJE+1,:) = ZPHIT(:,IJB,:)
+IF ( GWEST ) THEN
+ !$acc kernels async
+ ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB-1)
+ !$acc end kernels
+END IF
+IF ( GEAST ) THEN
+ !$acc kernels async
+ ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB-1)
+ !$acc end kernels
+END IF
+IF ( GSOUTH ) THEN
+ !$acc kernels async
+ ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB-1)
+ !$acc end kernels
END IF
+IF ( GNORTH ) THEN
+ !$acc kernels async
+ ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB-1)
+ !$acc end kernels
+END IF
+IF ( GSOUTH2D ) THEN
+ !$acc kernels async
+ ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)
+ !$acc end kernels
+END IF
+IF ( GNORTH2D ) THEN
+ !$acc kernels async
+ ZPHIT(:,IJE+1,:) = ZPHIT(:,IJB,:)
+ !$acc end kernels
+END IF
+!$acc wait
!
+#ifndef MNH_OPENACC
ZDV_SOURCE = GX_M_U(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX)
+#else
+CALL GX_M_U_DEVICE(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX,ZDV_SOURCE)
+#endif
!
-IF ( HLBCX(1) /= 'CYCL' ) THEN
- IF(LWEST_ll()) THEN
+IF ( GWEST ) THEN
!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
+ !$acc kernels
DO JK=2,IKU-1
ZDV_SOURCE(IIB,:,JK)= &
(ZPHIT(IIB,:,JK) - ZPHIT(IIB-1,:,JK) - 0.5 * ( &
@@ -580,10 +656,12 @@ IF ( HLBCX(1) /= 'CYCL' ) THEN
+PDZX(IIB,:,JK+1) * (ZPHIT(IIB,:,JK+1)-ZPHIT(IIB,:,JK)) / PDZZ(IIB,:,JK+1) &
) &
) / PDXX(IIB,:,JK)
- END DO
- ENDIF
+ END DO
+ !$acc end kernels
+ENDIF
!
- IF(LEAST_ll()) THEN
+IF( GEAST ) THEN
+ !$acc kernels
DO JK=2,IKU-1
ZDV_SOURCE(IIE+1,:,JK)= &
(ZPHIT(IIE+1,:,JK) - ZPHIT(IIE+1-1,:,JK) - 0.5 * ( &
@@ -593,27 +671,41 @@ IF ( HLBCX(1) /= 'CYCL' ) THEN
/ PDZZ(IIE+1-1,:,JK+1) &
) &
) / PDXX(IIE+1,:,JK)
- END DO
- END IF
+ END DO
+ !$acc end kernels
END IF
!
CALL MPPDB_CHECK3DM("before MXM PRESSUREZ :PRU/V/WS",PRECISION,PRUS,PRVS,PRWS)
-IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- PRUS = PRUS - MXM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
- PRWS = PRWS - MZM(PRHODJ * XCPD * ZTHETAV) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
+IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
+!!$ PRUS = PRUS - MXM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
+!!$ PRWS = PRWS - MZM(PRHODJ * XCPD * ZTHETAV) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
+ !$acc kernels
+ ZPRHODJ = PRHODJ * XCPD * ZTHETAV
+ !$acc end kernels
+ CALL MXM_DEVICE(ZPRHODJ, ZMXM_PRHODJ)
+ CALL MZM_DEVICE(ZPRHODJ, ZMZM_PRHODJ)
+ CALL GZ_M_W_DEVICE(1,IKU,1,ZPHIT,PDZZ,ZGZ_M_W)
+ !$acc kernels
+ PRUS = PRUS - ZMXM_PRHODJ * ZDV_SOURCE
+ PRWS = PRWS - ZMZM_PRHODJ * ZGZ_M_W
+ !$acc end kernels
ELSEIF(CEQNSYS=='LHE') THEN
PRUS = PRUS - MXM(PRHODJ) * ZDV_SOURCE
PRWS = PRWS - MZM(PRHODJ) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
END IF
!
IF(.NOT. L2D) THEN
+ !
+#ifndef MNH_OPENACC
+ ZDV_SOURCE = GY_M_V(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY)
+#else
+ CALL GY_M_V_DEVICE(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY,ZDV_SOURCE)
+#endif
!
- ZDV_SOURCE = GY_M_V(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY)
-!
- IF ( HLBCY(1) /= 'CYCL' ) THEN
- IF (LSOUTH_ll()) THEN
+ IF ( GSOUTH ) THEN
!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
+ !$acc kernels async
DO JK=2,IKU-1
ZDV_SOURCE(:,IJB,JK)= &
(ZPHIT(:,IJB,JK) - ZPHIT(:,IJB-1,JK) - 0.5 * ( &
@@ -621,10 +713,12 @@ IF(.NOT. L2D) THEN
+PDZY(:,IJB,JK+1) * (ZPHIT(:,IJB,JK+1)-ZPHIT(:,IJB,JK)) / PDZZ(:,IJB,JK+1) &
) &
) / PDYY(:,IJB,JK)
- END DO
+ END DO
+ !$acc end kernels
END IF
!
- IF (LNORTH_ll()) THEN
+ IF ( GNORTH ) THEN
+ !$acc kernels async
DO JK=2,IKU-1
ZDV_SOURCE(:,IJE+1,JK)= &
(ZPHIT(:,IJE+1,JK) - ZPHIT(:,IJE+1-1,JK) - 0.5 * ( &
@@ -634,13 +728,18 @@ IF(.NOT. L2D) THEN
/ PDZZ(:,IJE+1-1,JK+1) &
) &
) / PDYY(:,IJE+1,JK)
- END DO
- END IF
+ END DO
+ !$acc end kernels
END IF
+!$acc wait
!
CALL MPPDB_CHECK3DM("before MYM PRESSUREZ :PRU/V/WS",PRECISION,PRUS,PRVS,PRWS)
IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
+!!$ PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
+ CALL MYM_DEVICE(ZPRHODJ,ZMYM_PRHODJ)
+ !$acc kernels
+ PRVS = PRVS - ZMYM_PRHODJ * ZDV_SOURCE
+ !$acc end kernels
ELSEIF(CEQNSYS=='LHE') THEN
PRVS = PRVS - MYM(PRHODJ) * ZDV_SOURCE
END IF
@@ -649,10 +748,12 @@ END IF
!! same boundary conditions as in gdiv ... !! (provisory coding)
!! (necessary when NVERB=1)
!!
+ !$acc kernels
PRUS(:,:,IKB-1)=PRUS(:,:,IKB)
PRUS(:,:,IKE+1)=PRUS(:,:,IKE)
PRVS(:,:,IKB-1)=PRVS(:,:,IKB)
PRVS(:,:,IKE+1)=PRVS(:,:,IKE)
+ !$acc end kernels
!
NULLIFY(TZFIELDS2_ll)
CALL ADD3DFIELD_ll( TZFIELDS2_ll, PRUS, 'PRESSUREZ::PRUS' )
@@ -665,10 +766,14 @@ CALL CLEANLIST_ll(TZFIELDS2_ll)
CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
!
IF ( CEQNSYS=='DUR' ) THEN
- IF ( SIZE(PRVREF,1) == 0 ) THEN
+ IF ( GPRVREF0 ) THEN
+ !$acc kernels
ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF
- ELSE
+ !$acc end kernels
+ ELSE
+ !$acc kernels
ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF*(1.+PRVREF)
+ !$acc end kernels
END IF
ELSEIF( CEQNSYS=='MAE' .OR. CEQNSYS=='LHE' ) THEN
ZDV_SOURCE=ZDV_SOURCE/PRHODJ*PRHODREF
@@ -729,61 +834,83 @@ IF ((ZMAX_ll > 1.E-12) .AND. KTCOUNT >0 ) THEN
PRVREF, PEXNREF, ZPHIT )
!
IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD)
+#ifndef MNH_OPENACC
+ PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD)
+#else
+ !$acc kernels
+ PPABST(:,:,:)=XP00*BR_POW((ZPHIT+PEXNREF),(XCPD/XRD))
+ !$acc end kernels
+#endif
ELSEIF(CEQNSYS=='LHE') THEN
PPABST(:,:,:)=XP00*(ZPHIT/(XCPD*PTHVREF)+PEXNREF)**(XCPD/XRD)
ENDIF
!
- IF( HLBCX(1) == 'CYCL' ) THEN
- IF (LWEST_ll()) THEN
- ZPABS_W(:,:)= PPABST(IIB,:,:)
- END IF
-!
- IF (LEAST_ll()) THEN
- ZPABS_E(:,:)= PPABST(IIE+1,:,:)
- END IF
-!
+ IF( GWEST ) THEN
+ !$acc kernels async
+ ZPABS_W(:,:)= PPABST(IIB,:,:)
+ !$acc end kernels
END IF
!
- IF( HLBCY(1) == 'CYCL' ) THEN
- IF (LSOUTH_ll()) THEN
- ZPABS_S(:,:)= PPABST(:,IJB,:)
- END IF
-!
- IF (LNORTH_ll()) THEN
- ZPABS_N(:,:)= PPABST(:,IJE+1,:)
- END IF
+ IF ( GEAST ) THEN
+ !$acc kernels async
+ ZPABS_E(:,:)= PPABST(IIE+1,:,:)
+ !$acc end kernels
+ END IF
!
+ IF( GSOUTH ) THEN
+ !$acc kernels async
+ ZPABS_S(:,:)= PPABST(:,IJB,:)
+ !$acc end kernels
END IF
!
+ IF ( GNORTH ) THEN
+ !$acc kernels async
+ ZPABS_N(:,:)= PPABST(:,IJE+1,:)
+ !$acc end kernels
+ END IF
+ !
+ !$acc wait
+ !
+#ifndef MNH_OPENACC
CALL ADD3DFIELD_ll( TZFIELDS_ll, PPABST, 'PRESSUREZ::PPABST' )
CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS_ll)
-!
- IF( HLBCX(1) == 'CYCL' ) THEN
- IF (LWEST_ll()) THEN
- PPABST(IIB,:,:) = ZPABS_W(:,:)
- END IF
-!
- IF (LEAST_ll()) THEN
- PPABST(IIE+1,:,:) = ZPABS_E(:,:)
- END IF
-!
+#else
+ CALL GET_HALO_D( PPABST,HNAME='PRESSUREZ::PPABST' )
+#endif
+
+!
+ IF( GWEST ) THEN
+ !$acc kernels async
+ PPABST(IIB,:,:) = ZPABS_W(:,:)
+ !$acc end kernels
END IF
!
- IF( HLBCY(1) == 'CYCL' ) THEN
- IF (LSOUTH_ll()) THEN
- PPABST(:,IJB,:) = ZPABS_S(:,:)
- END IF
+ IF ( GEAST ) THEN
+ !$acc kernels async
+ PPABST(IIE+1,:,:) = ZPABS_E(:,:)
+ !$acc end kernels
+ END IF
!
- IF (LNORTH_ll()) THEN
- PPABST(:,IJE+1,:) = ZPABS_N(:,:)
- END IF
+ IF( GSOUTH ) THEN
+ !$acc kernels async
+ PPABST(:,IJB,:) = ZPABS_S(:,:)
+ !$acc end kernels
+ END IF
!
+ IF ( GNORTH ) THEN
+ !$acc kernels async
+ PPABST(:,IJE+1,:) = ZPABS_N(:,:)
+ !$acc end kernels
END IF
!
+!$acc wait
+!
END IF
!
+#ifdef MNH_OPENACC
+CALL MNH_REL_ZT3D ( IZPRHODJ,IZMXM_PRHODJ,IZMZM_PRHODJ,IZGZ_M_W,IZMYM_PRHODJ )
+#endif
!-------------------------------------------------------------------------------
!
END SUBROUTINE PRESSUREZ
diff --git a/src/ZSOLVER/qlap.f90 b/src/ZSOLVER/qlap.f90
index d4da491d97bc6c716267088375a88a11fbeac7cc..5f76e50d9b5f77db9717497307716ce2d7e2fef0 100644
--- a/src/ZSOLVER/qlap.f90
+++ b/src/ZSOLVER/qlap.f90
@@ -33,6 +33,33 @@ REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PQLAP ! final divergence
!
END FUNCTION QLAP
!
+!
+#ifdef MNH_OPENACC
+SUBROUTINE QLAP_DEVICE(PQLAP,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY)
+!
+IMPLICIT NONE
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQLAP ! final divergence
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at time t
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components
+!
+
+!
+END SUBROUTINE QLAP_DEVICE
+#endif
END INTERFACE
!
END MODULE MODI_QLAP
@@ -262,3 +289,280 @@ CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP)
!-------------------------------------------------------------------------------
!
END FUNCTION QLAP
+
+#ifdef MNH_OPENACC
+! #########################################################################
+ SUBROUTINE QLAP_DEVICE(PQLAP,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY)
+! #########################################################################
+!
+!!**** *QLAP * - compute the complete quasi-laplacien QLAP of a field P
+!!
+!! PURPOSE
+!! -------
+! This function computes the quasi-laplacien QLAP of the scalar field P
+! localized at a mass point, with non-vanishing orography.
+! The result is localized at a mass point and defined by:
+! for Durran and MAE anelastic equations
+! ( ( GX_M_U (PY) ) )
+! PQLAP = GDIV ( rho * CPd * Thetav * J ( GX_M_V (PY) ) )
+! ( ( GX_M_W (PY) ) )
+! or for Lipps and Hemler
+! ( ( GX_M_U (PY) ) )
+! PQLAP = GDIV ( rho * J ( GX_M_V (PY) ) )
+! ( ( GX_M_W (PY) ) )
+! Where GX_M_.. are the cartesian components of the gradient of PY and
+! GDIV is the operator acting on a vector AA:
+!
+! GDIV ( AA ) = J * divergence (1/J AA )
+!
+!!** METHOD
+!! ------
+!! First, we compute the gradients along x, y , z of the P field. The
+!! result is multiplied by rhod * CPd * Thetav or rhod depending on the
+!! chosen anelastic system where the suffixes indicate
+!! d dry and v for virtual.
+!! Then, the pseudo-divergence ( J * DIV (1/J o ) ) is computed by the
+!! subroutine GDIV. The result is localized at a mass point.
+!!
+!! EXTERNAL
+!! --------
+!! Function GX_M_U : compute the gradient along x
+!! Function GY_M_V : compute the gradient along y
+!! Function GZ_M_W : compute the gradient along z
+!! FUNCTION MXM: compute an average in the x direction for a variable
+!! at a mass localization
+!! FUNCTION MYM: compute an average in the y direction for a variable
+!! at a mass localization
+!! FUNCTION MZM: compute an average in the z direction for a variable
+!! at a mass localization
+!! Subroutine GDIV : compute J times the divergence of 1/J times a vector
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: JPHEXT, JPVEXT
+!! Module MODD_CONF: L2D,CEQNSYS
+!! Module MODD_CST : XCPD
+!!
+!! REFERENCE
+!! ---------
+!! Pressure solver documentation ( Scientific documentation )
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 11/07/94
+!! Modification 16/03/95 change the argument list of the gradient
+!! 14/01/97 New anelastic equation ( Stein )
+!! 17/12/97 include the case of non-vanishing orography
+!! at the lbc ( Stein )
+!! 06/12 V.Masson : update_halo due to CONTRAV changes
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODE_ll
+!
+USE MODD_PARAMETERS
+USE MODD_CONF
+USE MODD_CST
+USE MODI_GDIV
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN
+USE MODI_SHUMAN_DEVICE
+!
+USE MODE_MPPDB
+!
+USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D
+!
+USE MODI_GET_HALO
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of arguments
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQLAP ! final divergence
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at time t
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components
+!
+!* 0.2 declarations of local variables
+!
+REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZU ! rho*J*gradient along x
+!
+REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZV ! rho*J*gradient along y
+!
+REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZW ! rho*J*gradient along z
+!
+INTEGER :: IZU,IZV,IZW
+!
+INTEGER :: IIU,IJU,IKU ! I,J,K array sizes
+INTEGER :: JK,JJ,JI ! vertical loop index
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+INTEGER :: IINFO_ll
+INTEGER :: IIB,IIE,IJB,IJE
+!
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMXM,ZMYM,ZMZM,ZRHODJ,ZGZMW
+INTEGER :: IZMXM,IZMYM,IZMZM,IZRHODJ,IZGZMW
+!
+LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. COMPUTE THE GRADIENT COMPONENTS
+! -------------------------------
+!
+CALL GET_DIM_EXT_ll('B',IIU,IJU)
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IKU=SIZE(PY,3)
+!
+GWEST = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() )
+GEAST = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() )
+GSOUTH = ( HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() )
+GNORTH = ( HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() )
+!
+IZU = MNH_ALLOCATE_ZT3D( ZU,IIU,IJU,IKU )
+IZV = MNH_ALLOCATE_ZT3D( ZV,IIU,IJU,IKU )
+IZW = MNH_ALLOCATE_ZT3D( ZW,IIU,IJU,IKU )
+!
+IZMXM = MNH_ALLOCATE_ZT3D( ZMXM,IIU,IJU,IKU )
+IZMYM = MNH_ALLOCATE_ZT3D( ZMYM,IIU,IJU,IKU )
+IZMZM = MNH_ALLOCATE_ZT3D( ZMZM,IIU,IJU,IKU )
+IZRHODJ = MNH_ALLOCATE_ZT3D( ZRHODJ ,IIU,IJU,IKU )
+IZGZMW = MNH_ALLOCATE_ZT3D( ZGZMW ,IIU,IJU,IKU )
+!
+CALL GX_M_U_DEVICE(1,IKU,1,PY,PDXX,PDZZ,PDZX,ZU)
+CALL MPPDB_CHECK3D(ZU,'QLAP::ZU',PRECISION)
+!
+IF ( GWEST ) THEN
+ !$acc kernels async
+ DO JK=2,IKU-1
+ DO JJ=1,IJU
+ ZU(IIB,JJ,JK)= (PY(IIB,JJ,JK) - PY(IIB-1,JJ,JK) - 0.5 * ( &
+ PDZX(IIB,JJ,JK) * (PY(IIB,JJ,JK)-PY(IIB,JJ,JK-1)) / PDZZ(IIB,JJ,JK) &
+ +PDZX(IIB,JJ,JK+1) * (PY(IIB,JJ,JK+1)-PY(IIB,JJ,JK)) / PDZZ(IIB,JJ,JK+1) &
+ ) ) / PDXX(IIB,JJ,JK)
+ END DO
+ END DO
+ !$acc end kernels
+END IF
+!
+IF ( GEAST ) THEN
+ !$acc kernels async
+ DO JK=2,IKU-1
+ DO JJ=1,IJU
+ ZU(IIE+1,JJ,JK)= (PY(IIE+1,JJ,JK) - PY(IIE+1-1,JJ,JK) - 0.5 * ( &
+ PDZX(IIE+1,JJ,JK) * (PY(IIE+1-1,JJ,JK)-PY(IIE+1-1,JJ,JK-1)) / PDZZ(IIE+1-1,JJ,JK) &
+ +PDZX(IIE+1,JJ,JK+1) * (PY(IIE+1-1,JJ,JK+1)-PY(IIE+1-1,JJ,JK)) / PDZZ(IIE+1-1,JJ,JK+1)&
+ ) ) / PDXX(IIE+1,JJ,JK)
+ END DO
+ END DO
+ !$acc end kernels
+END IF
+!$acc wait
+!
+CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/W',PRECISION)
+CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/E',PRECISION)
+!
+IF(.NOT. L2D) THEN
+!
+ CALL GY_M_V_DEVICE(1,IKU,1,PY,PDYY,PDZZ,PDZY,ZV)
+ CALL MPPDB_CHECK3D(ZV,'QLAP::ZV',PRECISION)
+!
+ IF ( GSOUTH ) THEN
+ !$acc kernels async
+ DO JK=2,IKU-1
+ DO JI=1,IIU
+ ZV(JI,IJB,JK)= (PY(JI,IJB,JK) - PY(JI,IJB-1,JK) - 0.5 * ( &
+ PDZY(JI,IJB,JK) * (PY(JI,IJB,JK)-PY(JI,IJB,JK-1)) / PDZZ(JI,IJB,JK) &
+ +PDZY(JI,IJB,JK+1) * (PY(JI,IJB,JK+1)-PY(JI,IJB,JK)) / PDZZ(JI,IJB,JK+1) &
+ ) ) / PDYY(JI,IJB,JK)
+ END DO
+ END DO
+ !$acc end kernels
+ END IF
+
+
+ IF ( GNORTH ) THEN
+ !$acc kernels async
+ DO JK=2,IKU-1
+ DO JI=1,IIU
+ ZV(JI,IJE+1,JK)= (PY(JI,IJE+1,JK) - PY(JI,IJE+1-1,JK) - 0.5 * ( &
+ PDZY(JI,IJE+1,JK) * (PY(JI,IJE+1-1,JK)-PY(JI,IJE+1-1,JK-1)) &
+ / PDZZ(JI,IJE+1-1,JK) &
+ +PDZY(JI,IJE+1,JK+1) * (PY(JI,IJE+1-1,JK+1)-PY(JI,IJE+1-1,JK)) &
+ / PDZZ(JI,IJE+1-1,JK+1)&
+ ) ) / PDYY(JI,IJE+1,JK)
+ END DO
+ END DO
+ !$acc end kernels
+ END IF
+ !$acc wait
+ CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/S',PRECISION)
+ CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/N',PRECISION)
+!
+ELSE
+ ZV=0.
+ENDIF
+!
+IF ( CEQNSYS == 'DUR' .OR. CEQNSYS == 'MAE' ) THEN
+ !$acc kernels
+ ZRHODJ = PRHODJ * XCPD * PTHETAV
+ !$acc end kernels
+ CALL MXM_DEVICE (ZRHODJ, ZMXM )
+ !$acc kernels
+ ZU = ZMXM * ZU
+ !$acc end kernels
+ IF(.NOT. L2D) THEN
+ CALL MYM_DEVICE (ZRHODJ, ZMYM )
+ !$acc kernels
+ ZV = ZMYM * ZV
+ !$acc end kernels
+ END IF
+ CALL MZM_DEVICE (ZRHODJ, ZMZM )
+ CALL GZ_M_W_DEVICE(1,IKU,1,PY,PDZZ,ZGZMW)
+ !$acc kernels
+ ZW = ZMZM * ZGZMW
+ !$acc end kernels
+ELSEIF ( CEQNSYS == 'LHE' ) THEN
+ ZU = MXM(PRHODJ) * ZU
+ IF(.NOT. L2D) THEN
+ ZV = MYM(PRHODJ) * ZV
+ ENDIF
+ ZW = MZM(PRHODJ) * GZ_M_W(1,IKU,1,PY,PDZZ)
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTE THE DIVERGENCE
+! ----------------------
+!
+CALL GET_HALO_D(ZU,HNAME='QLAP::ZU')
+CALL GET_HALO_D(ZV,HNAME='QLAP::ZV')
+CALL GET_HALO_D(ZW,HNAME='QLAP::ZW')
+!
+CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP)
+!
+CALL MNH_REL_ZT3D ( IZU,IZV,IZW, IZMXM,IZMYM,IZMZM,IZRHODJ,IZGZMW )
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE QLAP_DEVICE
+#endif
diff --git a/src/ZSOLVER/zsolver.f90 b/src/ZSOLVER/zsolver.f90
index 67c35779264fbbf9c63b31a89c940a7c27ad8c41..779a6b5685912b69249552d79b93337d17edd1bf 100644
--- a/src/ZSOLVER/zsolver.f90
+++ b/src/ZSOLVER/zsolver.f90
@@ -227,8 +227,13 @@ REAL :: ZDOT_DELTA ! dot product of ZDELTA by itself
!
!
!* 1.1 compute the vector: r^(0) = Q(PHI) - Y
-!
+!
+#ifndef MNH_OPENACC
ZRESIDUE = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY
+#else
+CALL QLAP_DEVICE(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI)
+ZRESIDUE = ZRESIDUE - PY
+#endif
!
!* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y )
!
@@ -240,7 +245,11 @@ CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
!
!* 1.3 compute the vector: delta^(0) = Q ( p^(0) )
!
+#ifndef MNH_OPENACC
ZDELTA = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP)
+#else
+CALL QLAP_DEVICE(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP)
+#endif
!
!-------------------------------------------------------------------------------
!
@@ -256,15 +265,19 @@ DO JM = 1,KITR
!
!* 2.2 update the pressure function PHI
!
+ !$acc kernels
PPHI = PPHI + ZLAMBDA * ZP
+ !$acc end kernels
!
!
IF( JM == KITR ) EXIT
!
!
!* 2.3 update the residual error: r
-!
+ !
+ !$acc kernels
ZRESIDUE = ZRESIDUE + ZLAMBDA * ZDELTA
+ !$acc end kernels
!
!* 2.4 compute the vector: q = F^(-1)*( Q(PHI) - Y )
!
@@ -276,7 +289,11 @@ DO JM = 1,KITR
!
!* 2.5 compute the auxiliary field: ksi = Q ( q )
!
+#ifndef MNH_OPENACC
ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ)
+#else
+ CALL QLAP_DEVICE(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ)
+#endif
! -1
!* 2.6 compute the step ALPHA
!
@@ -284,8 +301,10 @@ DO JM = 1,KITR
!
!* 2.7 update p and DELTA
!
+ !$acc kernels
ZP = ZQ + ZALPHA * ZP
ZDELTA = ZKSI + ZALPHA * ZDELTA
+ !$acc end kernels
!
END DO ! end of the loop for the iterative solver
!