From 71e0d56f80b8c031173e12a6364aa447d76d2b1f Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Thu, 1 Aug 2019 10:13:11 +0200
Subject: [PATCH] Philippe 01/08/2019: OpenACC: implementation of slow_terms
---
src/MNH/slow_terms.f90 | 228 +++++++++++++++++++++++++++++++----------
1 file changed, 175 insertions(+), 53 deletions(-)
diff --git a/src/MNH/slow_terms.f90 b/src/MNH/slow_terms.f90
index 3699b5af8..5291742d1 100644
--- a/src/MNH/slow_terms.f90
+++ b/src/MNH/slow_terms.f90
@@ -156,12 +156,16 @@ END MODULE MODI_SLOW_TERMS
!* 0. DECLARATIONS
! ------------
!
-USE MODD_PARAMETERS
-USE MODD_CLOUDPAR
-USE MODD_CST
-USE MODD_CONF
-USE MODD_BUDGET
-!
+USE MODD_BUDGET, only: LBUDGET_RC, LBUDGET_RR, LBUDGET_RV, LBUDGET_TH
+USE MODD_CLOUDPAR, only: XC1RC, XC2RC, XC1RE, XC2RE, XCEXRA, XCEXRE, XCEXRS, XCEXVT, XCRA, XCRS, XDIVA, XTHCO
+USE MODD_CST, only: XALPW, XBETAW, XGAMW, XCL, XCPD, XCPV, XLVTT, XMD, XMV, XP00, XRD, XRHOLW, XRV, XTT
+USE MODD_PARAMETERS, only: JPVEXT
+
+use mode_mppdb
+
+#ifdef MNH_BITREP
+use modi_bitrep
+#endif
USE MODI_BUDGET
!
IMPLICIT NONE
@@ -207,14 +211,49 @@ INTEGER :: IKE ! the microphysical sources have to be computed
!
REAL :: ZTSPLITR ! Small time step for rain sedimentation
!
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) &
- :: ZT,ZW,ZW1,ZW2,ZW3,ZEXNT,ZDZZ ! Work arrays
-LOGICAL, DIMENSION(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3)) :: G3D
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZT,ZW,ZW1,ZW2,ZW3,ZEXNT,ZDZZ ! Work arrays
+LOGICAL, DIMENSION(:,:,:), ALLOCATABLE :: G3D
!
-INTEGER , DIMENSION(:), ALLOCATABLE :: I1,I2 ! index for packed array
INTEGER :: JI,JJ,IC,JL ! loop control for packed array
-!
+
+!$acc declare present( PZZ, PRHODJ, PRHODREF, PCLDFR, PTHT, PRVT, PRCT, PRRT, &
+!$acc & PPABST, PTHS, PRVS, PRCS, PRRS, PINPRR, PINPRR3D, PEVAP3D )
+
+!$acc declare create( ZT, ZW, ZW1, ZW2, ZW3, ZEXNT, ZDZZ, G3D )
+
+!$acc declare copyin( XC1RC, XC2RC, XC1RE, XC2RE, XCEXRA, XCEXRE, XCEXRS, XCEXVT, XCRA, XCRS, XDIVA, XTHCO, &
+!$acc & XALPW, XBETAW, XGAMW, XCL, XCPD, XCPV, XLVTT, XMD, XMV, XP00, XRD, XRHOLW, XRV, XTT)
+
!-------------------------------------------------------------------------------
+IF (MPPDB_INITIALIZED) THEN
+ !Check all IN arrays
+ CALL MPPDB_CHECK(PZZ, "SLOW_TERMS beg:PZZ")
+ CALL MPPDB_CHECK(PRHODJ, "SLOW_TERMS beg:PRHODJ")
+ CALL MPPDB_CHECK(PRHODREF,"SLOW_TERMS beg:PRHODREF")
+ CALL MPPDB_CHECK(PCLDFR, "SLOW_TERMS beg:PCLDFR")
+ CALL MPPDB_CHECK(PTHT, "SLOW_TERMS beg:PTHT")
+ CALL MPPDB_CHECK(PRVT, "SLOW_TERMS beg:PRVT")
+ CALL MPPDB_CHECK(PRCT, "SLOW_TERMS beg:PRCT")
+ CALL MPPDB_CHECK(PRRT, "SLOW_TERMS begPRRT:")
+ CALL MPPDB_CHECK(PPABST, "SLOW_TERMS beg:PPABST")
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PTHS, "SLOW_TERMS beg:PTHS")
+ CALL MPPDB_CHECK(PRVS, "SLOW_TERMS beg:PRVS")
+ CALL MPPDB_CHECK(PRCS, "SLOW_TERMS beg:PRCS")
+ CALL MPPDB_CHECK(PRRS, "SLOW_TERMS beg:PRRS")
+ CALL MPPDB_CHECK(PINPRR, "SLOW_TERMS beg:PINPRR")
+ CALL MPPDB_CHECK(PINPRR3D,"SLOW_TERMS beg:PINPRR3D")
+ CALL MPPDB_CHECK(PEVAP3D, "SLOW_TERMS beg:PEVAP3D")
+END IF
+
+allocate( zt ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( zw ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( zw1 ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( zw2 ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( zw3 ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( zexnt( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( zdzz ( size(pzz,1), size(pzz,2), size(pzz,3) ) )
+allocate( g3d ( size(prhodj,1), size(prhodj,2), size(prhodj,3) ) )
!
!* 1. COMPUTE THE LOOP BOUNDS AND EXNER FUNCTION
! ------------------------------------------
@@ -222,10 +261,10 @@ INTEGER :: JI,JJ,IC,JL ! loop control for packed ar
IKB=1+JPVEXT
IKE=SIZE(PZZ,3) - JPVEXT
! compute Delta z at the mass point
+!$acc kernels
DO JK = 1,IKE
ZDZZ(:,:,JK) = PZZ(:,:,JK+1)-PZZ(:,:,JK)
END DO
-
!
!-------------------------------------------------------------------------------
!
@@ -241,30 +280,50 @@ ZW1(:,:,:) = PRRS(:,:,:) * PTSTEP
ZW2(:,:,:) = 0.
ZW3(:,:,:) = 0.
!
+#ifndef _OPENACC
G3D(:,:,:)=.FALSE.
DO JK=IKE,1,-1
WHERE( ZW1(:,:,JK)>0. .OR. G3D(:,:,JK+1) )
G3D(:,:,JK)=.TRUE.
END WHERE
END DO
+#else
+g3d(:, :, : ) = .false.
+do jk = ike, 1, -1
+ do jj = 1, size(pzz,2)
+ do ji = 1, size(pzz,1)
+ if ( zw1(ji, jj, jk ) >0. .or. g3d(ji, jj, jk+1 ) ) then
+ g3d(ji, jj, jk ) = .true.
+ end if
+ end do
+ end do
+end do
+#endif
!
WHERE (G3D(:,:,:))
+#ifndef MNH_BITREP
ZW3(:,:,:) = PRHODREF(:,:,:) ** (XCEXRS-XCEXVT)
+#else
+ ZW3(:,:,:) = BR_POW( PRHODREF(:,:,:), XCEXRS - XCEXVT )
+#endif
END WHERE
!
WHERE (ZW1(:,:,IKE+1)>0.)
+#ifndef MNH_BITREP
ZW2(:,:,IKE+1) = XCRS &
* ZW1(:,:,IKE+1) ** XCEXRS &
* PRHODREF(:,:,IKE+1) ** (XCEXRS-XCEXVT)
+#else
+ ZW2(:,:,IKE+1) = XCRS &
+ * BR_POW( ZW1(:,:,IKE+1), XCEXRS ) &
+ * BR_POW( PRHODREF(:,:,IKE+1), XCEXRS - XCEXVT )
+#endif
END WHERE
!
!
!* 2.2 small time step integration
! ---------------------------
!
-ALLOCATE (I1(SIZE(ZW1)))
-ALLOCATE (I2(SIZE(ZW1)))
-!
DO JN=1,KSPLITR
!
!* 2.2.1 test where computations should be made and pack arrays
@@ -272,37 +331,31 @@ DO JN=1,KSPLITR
!
DO JK=IKE,IKB,-1
- !
- IC = 0
+!$acc loop collapse(2) independent
DO JJ = 1,SIZE( ZW1,2)
DO JI = 1,SIZE( ZW1,1)
IF ( ( ZW1(JI,JJ,JK+1)>0. ) .AND. ( ZW1(JI,JJ,JK)>0. ) ) THEN
-!!$ IF ( (ZW1(JI,JJ,JK)+ZW1(JI,JJ,JK+1)>0.) ) THEN
- IC = IC +1
- I1(IC) = JI
- I2(IC) = JJ
- ENDIF
- ENDDO
- ENDDO
!
!* 2.2.2 compute the rain flux
! ---------------------
!
!
- DO JL=1,IC
- ZW2(I1(JL),I2(JL),JK) = XCRS &
- * ZW1(I1(JL),I2(JL),JK) ** XCEXRS &
- * ZW3(I1(JL),I2(JL),JK)
+#ifndef MNH_BITREP
+ ZW2(JI,JJ,JK) = XCRS * ( ZW1(JI,JJ,JK) ** XCEXRS ) * ZW3(JI,JJ,JK)
+#else
+ ZW2(JI,JJ,JK) = XCRS * BR_POW( ZW1(JI,JJ,JK), XCEXRS ) * ZW3(JI,JJ,JK)
+#endif
!
!
!* 2.2.3 update the rain tendency with the small timestep
! ------------------------------------------------
-!
- ZW1(I1(JL),I2(JL),JK) = ZW1(I1(JL),I2(JL),JK) + ZTSPLITR * &
- ( ZW2(I1(JL),I2(JL),JK+1)-ZW2(I1(JL),I2(JL),JK) ) / &
- ( ZDZZ(I1(JL),I2(JL),JK) * PRHODREF(I1(JL),I2(JL),JK) )
- ENDDO
!
+ ZW1(JI,JJ,JK) = ZW1(JI,JJ,JK) + ZTSPLITR * &
+ ( ZW2(JI,JJ,JK+1)-ZW2(JI,JJ,JK) ) / &
+ ( ZDZZ(JI,JJ,JK) * PRHODREF(JI,JJ,JK) )
+ END IF
+ END DO
+ ENDDO
ENDDO
!
!* 2.2.4 compute the explicit accumulated precipitations
@@ -315,16 +368,17 @@ DO JN=1,KSPLITR
!
END DO
!
-DEALLOCATE (I1)
-DEALLOCATE (I2)
-!
!* 2.4 update the rain tendency
!
PRRS(:,:,:) = ZW1(:,:,:) / PTSTEP
-!
+!$acc end kernels
+!
!* 2.5 budget storage
!
-IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'SEDI_BU_RRR')
+IF (LBUDGET_RR) THEN
+!$acc update self(PRRS)
+ CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'SEDI_BU_RRR')
+END IF
!
!-------------------------------------------------------------------------------
!
@@ -335,21 +389,35 @@ IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'SEDI_BU_RRR')
!
!* 3.1 compute the accretion and update the tendencies
!
-WHERE ( (PRCT(:,:,:)>0.0) .AND. (PRRT(:,:,:)>0.0) &
- .AND. (PRCS(:,:,:)>0.0) )
+!$acc kernels
+G3D(:,:,:) = PRCT(:,:,:)>0.0 .AND. PRRT(:,:,:)>0.0 .AND. PRCS(:,:,:)>0.0
+WHERE ( G3D(:,:,:) )
+#ifndef MNH_BITREP
ZW(:,:,:) = XCRA * PRCT(:,:,:) &
* PRRT(:,:,:) ** XCEXRA &
* PRHODREF(:,:,:) ** (XCEXRA - XCEXVT)
+#else
+ ZW(:,:,:) = XCRA * PRCT(:,:,:) &
+ * BR_POW( PRRT(:,:,:), XCEXRA ) &
+ * BR_POW( PRHODREF(:,:,:), XCEXRA - XCEXVT )
+#endif
ZW(:,:,:) = MIN ( ZW(:,:,:) , PRCS(:,:,:) )
!
PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
END WHERE
+!$acc end kernels
!
!* 3.2 budget storage
!
-IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'ACCR_BU_RRC')
-IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'ACCR_BU_RRR')
+IF (LBUDGET_RC) THEN
+!$acc update self(PRCS)
+ CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'ACCR_BU_RRC')
+END IF
+IF (LBUDGET_RR) THEN
+!$acc update self(PRRS)
+ CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'ACCR_BU_RRR')
+END IF
!
!-------------------------------------------------------------------------------
!
@@ -360,15 +428,18 @@ IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'ACCR_BU_RRR')
!
!* 4.1 compute the autoconversion and update the tendencies
!
+!$acc kernels
IF ( HSUBG_AUCV == 'CLFR' ) THEN
- WHERE ( (PRCT(:,:,:)>0.0) .AND. (PRCS(:,:,:)>0.0) .AND. (PCLDFR(:,:,:)>0.0) )
+ G3D(:,:,:) = PRCT(:,:,:)>0.0 .AND. PRCS(:,:,:)>0.0 .AND. PCLDFR(:,:,:)>0.0
+ WHERE ( G3D(:,:,:) )
ZW(:,:,:) = XC1RC * MAX(PRCT(:,:,:) /(PCLDFR(:,:,:)) - XC2RC / PRHODREF(:,:,:),0.)
ZW(:,:,:) = MIN ( (ZW(:,:,:)* PCLDFR(:,:,:)), PRCS(:,:,:) )
PRCS(:,:,:) = PRCS(:,:,:) - ZW(:,:,:)
PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
END WHERE
ELSE
- WHERE ( (PRCT(:,:,:)>0.0) .AND. (PRCS(:,:,:)>0.0) )
+ G3D(:,:,:) = PRCT(:,:,:)>0.0 .AND. PRCS(:,:,:)>0.0
+ WHERE ( G3D(:,:,:) )
ZW(:,:,:) = XC1RC * MAX ( PRCT(:,:,:) - XC2RC / PRHODREF(:,:,:), 0. )
ZW(:,:,:) = MIN ( ZW(:,:,:) , PRCS(:,:,:) )
!
@@ -376,23 +447,36 @@ ELSE
PRRS(:,:,:) = PRRS(:,:,:) + ZW(:,:,:)
END WHERE
END IF
+!$acc end kernels
!
!* 4.2 budget storage
!
-IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'AUTO_BU_RRC')
-IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'AUTO_BU_RRR')
+IF (LBUDGET_RC) THEN
+!$acc update self(PRCS)
+ CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'AUTO_BU_RRC')
+END IF
+IF (LBUDGET_RR) THEN
+!$acc update self(PRRS)
+ CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'AUTO_BU_RRR')
+END IF
!
!-------------------------------------------------------------------------------
!
!* 5. COMPUTES THE RAIN EVAPORATION (RE) SOURCE
! -----------------------------------------
!
+!$acc kernels
PEVAP3D(:,:,:)=0.
-WHERE ( (PRRT(:,:,:)>0.0) .AND. (PRCT(:,:,:)==0.0) )
+G3D(:,:,:) = PRRT(:,:,:)>0.0 .AND. PRCT(:,:,:)==0.0
+WHERE ( G3D(:,:,:) )
!
!* 5.1 compute the Exner function
!
+#ifndef MNH_BITREP
ZEXNT(:,:,:) = (PPABST(:,:,:)/XP00)**(XRD/XCPD)
+#else
+ ZEXNT(:,:,:) = BR_POW( PPABST(:,:,:)/XP00, XRD/XCPD )
+#endif
!
!* 5.2 compute the temperature
!
@@ -400,7 +484,11 @@ WHERE ( (PRRT(:,:,:)>0.0) .AND. (PRCT(:,:,:)==0.0) )
!
!* 5.3 compute the saturation vapor pressure
!
+#ifndef MNH_BITREP
ZW1(:,:,:) = EXP( XALPW - XBETAW/ZT(:,:,:) - XGAMW*ALOG(ZT(:,:,:) ) )
+#else
+ ZW1(:,:,:) = BR_EXP( XALPW - XBETAW/ZT(:,:,:) - XGAMW*BR_LOG(ZT(:,:,:) ) )
+#endif
!
!* 5.4 compute the saturation mixing ratio
!
@@ -413,14 +501,25 @@ WHERE ( (PRRT(:,:,:)>0.0) .AND. (PRCT(:,:,:)==0.0) )
!
!* 5.6 compute the source
!
- ZW(:,:,:) = MAX( 1. - PRVT(:,:,:)/ZW2(:,:,:) , 0.0 ) * &
+#ifndef MNH_BITREP
+ ZW(:,:,:) = MAX( 1. - PRVT(:,:,:)/ZW2(:,:,:) , 0.0 ) * &
( XC1RE * SQRT( PRRT(:,:,:)/PRHODREF(:,:,:) ) &
+XC2RE * PRRT(:,:,:)**XCEXRE &
* PRHODREF(:,:,:)**(XCEXRE-0.5*XCEXVT-1.) &
) / &
( XRV*ZT(:,:,:)/(ZW1(:,:,:)*XDIVA) &
+ZW3(:,:,:) ** 2 / (XTHCO*XRV*ZT(:,:,:)**2) &
- )
+ )
+#else
+ ZW(:,:,:) = MAX( 1. - PRVT(:,:,:)/ZW2(:,:,:) , 0.0 ) * &
+ ( XC1RE * BR_POW( PRRT(:,:,:)/PRHODREF(:,:,:), 0.5 ) &
+ +XC2RE * BR_POW( PRRT(:,:,:), XCEXRE ) &
+ * BR_POW( PRHODREF(:,:,:), XCEXRE-0.5*XCEXVT-1. ) &
+ ) / &
+ ( XRV*ZT(:,:,:)/(ZW1(:,:,:)*XDIVA) &
+ +ZW3(:,:,:)*ZW3(:,:,:) / (XTHCO*XRV*ZT(:,:,:)*ZT(:,:,:)) &
+ )
+#endif
ZW(:,:,:) = MIN ( ZW(:,:,:) , PRRS(:,:,:) )
!
!* 5.7 update the total tendencies
@@ -433,22 +532,45 @@ WHERE ( (PRRT(:,:,:)>0.0) .AND. (PRCT(:,:,:)==0.0) )
+ XCL * (PRCT(:,:,:) + PRRT(:,:,:)) ) )
PEVAP3D(:,:,:)=ZW(:,:,:)
END WHERE
+!$acc end kernels
!
!* 5.8 budget storage
!
-IF (LBUDGET_RV) CALL BUDGET (PRVS(:,:,:)*PRHODJ(:,:,:),6,'REVA_BU_RRV')
-IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'REVA_BU_RRR')
-IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'REVA_BU_RTH')
+IF (LBUDGET_RV) THEN
+!$acc update self(PRVS)
+ CALL BUDGET (PRVS(:,:,:)*PRHODJ(:,:,:),6,'REVA_BU_RRV')
+END IF
+IF (LBUDGET_RR) THEN
+!$acc update self(PRRS)
+ CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'REVA_BU_RRR')
+END IF
+IF (LBUDGET_TH) THEN
+!$acc update self(PTHS)
+ CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'REVA_BU_RTH')
+END IF
!
!-------------------------------------------------------------------------------
!
!* 6. REMOVES NON-PHYSICAL LOW VALUES
! -------------------------------
!
-WHERE (PRRS(:,:,:)<1.E-16)
+!$acc kernels
+G3D(:,:,:) = PRRS(:,:,:) < 1.E-16
+WHERE ( G3D(:,:,:) )
PRRS(:,:,:)=0.
END WHERE
+!$acc end kernels
!-------------------------------------------------------------------------------
!
+IF (MPPDB_INITIALIZED) THEN
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PTHS, "SLOW_TERMS end:PTHS")
+ CALL MPPDB_CHECK(PRVS, "SLOW_TERMS end:PRVS")
+ CALL MPPDB_CHECK(PRCS, "SLOW_TERMS end:PRCS")
+ CALL MPPDB_CHECK(PRRS, "SLOW_TERMS end:PRRS")
+ CALL MPPDB_CHECK(PINPRR, "SLOW_TERMS end:PINPRR")
+ CALL MPPDB_CHECK(PINPRR3D,"SLOW_TERMS end:PINPRR3D")
+ CALL MPPDB_CHECK(PEVAP3D, "SLOW_TERMS end:PEVAP3D")
+END IF
!
END SUBROUTINE SLOW_TERMS
--
GitLab