Juan 15/04/2022:condensation.f90, Cray GPU , add JKPK,JKMK to integrate main...

Juan 15/04/2022:condensation.f90, Cray GPU , add JKPK,JKMK to integrate main JK loop(seq) in 3D do concurrent , mush faster !!!

src/MNH/condensation.f90

@@ -229,6 +229,9 @@ REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZINCP
 LOGICAL :: GPRESENT_PLV, GPRESENT_PLS, GPRESENT_PCPH
 LOGICAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: GWORK
 CHARACTER(LEN=4) :: YLAMBDA3 !Necessary to workaround NVHPC bug (version 21.7 if OpenACC enabled)
+LOGICAL :: GPRESENT_PHLC_HCF, GPRESENT_PHLC_HRC, GPRESENT_PHLI_HCF, GPRESENT_PHLI_HRI
+!
+INTEGER,    DIMENSION(:), POINTER, CONTIGUOUS :: JKPK,JKMK
 !
 !*       0.3  Definition of constants :
 !
@@ -252,6 +255,14 @@ REAL, DIMENSION(-22:11),PARAMETER :: ZSRC_1D =(/                         &
 !
 !-------------------------------------------------------------------------------
 !
+GPRESENT_PLV  = PRESENT(PLV)
+GPRESENT_PLS  = PRESENT(PLS)
+GPRESENT_PCPH = PRESENT(PCPH)
+GPRESENT_PHLC_HCF = PRESENT(PHLC_HCF)
+GPRESENT_PHLC_HRC = PRESENT(PHLC_HRC)
+GPRESENT_PHLI_HCF = PRESENT(PHLI_HCF)
+GPRESENT_PHLI_HRI = PRESENT(PHLI_HRI)
+!
 IF (MPPDB_INITIALIZED) THEN
   !Check all IN arrays
   CALL MPPDB_CHECK3D(PPABS,"CONDENSATION beg:PPABS",PRECISION)
@@ -260,9 +271,9 @@ IF (MPPDB_INITIALIZED) THEN
   CALL MPPDB_CHECK3D(PRG,"CONDENSATION beg:PRG",PRECISION)
   CALL MPPDB_CHECK3D(PSIGS,"CONDENSATION beg:PSIGS",PRECISION)
   CALL MPPDB_CHECK3D(PMFCONV,"CONDENSATION beg:PMFCONV",PRECISION)
-  IF (PRESENT(PLV)) CALL MPPDB_CHECK3D(PLV,"CONDENSATION beg:PLV",PRECISION)
-  IF (PRESENT(PLS)) CALL MPPDB_CHECK3D(PLS,"CONDENSATION beg:PLS",PRECISION)
-  IF (PRESENT(PCPH)) CALL MPPDB_CHECK3D(PCPH,"CONDENSATION beg:PCPH",PRECISION)
+  IF (GPRESENT_PLV) CALL MPPDB_CHECK3D(PLV,"CONDENSATION beg:PLV",PRECISION)
+  IF (GPRESENT_PLS) CALL MPPDB_CHECK3D(PLS,"CONDENSATION beg:PLS",PRECISION)
+  IF (GPRESENT_PCPH) CALL MPPDB_CHECK3D(PCPH,"CONDENSATION beg:PCPH",PRECISION)
   !Check all INOUT arrays
   CALL MPPDB_CHECK3D(PT,"CONDENSATION beg:PT",PRECISION)
   CALL MPPDB_CHECK3D(PRV,"CONDENSATION beg:PRV",PRECISION)
@@ -321,6 +332,8 @@ allocate( INQ1P(kiu, kju, kku )
 allocate( ZINCP(kiu, kju, kku )
 allocate( ZRCOLDP(kiu, kju, kku )
 allocate( ZRIOLDP(kiu, kju, kku )
+allocate( JKPK(kku) )
+allocate( JKMK(kku) )
 
 #else
 !Pin positions in the pools of MNH memory
@@ -371,23 +384,20 @@ CALL MNH_MEM_GET( INQ1P, kiu, kju, kku )
 CALL MNH_MEM_GET( ZINCP, kiu, kju, kku )
 CALL MNH_MEM_GET( ZRCOLDP, kiu, kju, kku )
 CALL MNH_MEM_GET( ZRIOLDP, kiu, kju, kku )
+CALL MNH_MEM_GET( JKPK, kku )
+CALL MNH_MEM_GET( JKMK, kku )
 
 !$acc data present( PPABS, PZZ, PT, PRV, PRC, PRI, PRS, PRG, PSIGS, PMFCONV, PCLDFR, PSIGRC, &
 !$acc &             ztlk, zrt, zl, zfrac, itpl, ztmin, zlv, zls,zcpd, gwork,&
 !$acc &             ZZZP,JKPP,ZPVP,ZQSLP,ZPIVP,ZQSIP,ZLVSP,ZAHP,ZAP,ZBP,ZSBARP,ZSIGMAP,&
 !$acc &             DZZP,ZDRWP,ZDTLP,ZLLP,ZSIG_CONVP,ZQ1P,ZGCONDP,ZGAUVP,ZCONDP,ZAUTCP,&
-!$acc &             ZGAUTCP,ZGAUCP,ZCRIAUTIP,ZAUTIP,ZGAUTIP,ZGAUIP,INQ1P,ZINCP,ZRCOLDP,ZRIOLDP)
+!$acc &             ZGAUTCP,ZGAUCP,ZCRIAUTIP,ZAUTIP,ZGAUTIP,ZGAUIP,INQ1P,ZINCP,ZRCOLDP,ZRIOLDP,&
+!$acc &             JKPK,JKMK )
 #endif
 
 IKTB=1+JPVEXT
 IKTE=KKU-JPVEXT
 !
-GPRESENT_PLV  = .FALSE.
-GPRESENT_PLS  = .FALSE.
-GPRESENT_PCPH = .FALSE.
-IF (PRESENT(PLV))  GPRESENT_PLV  = .TRUE.
-IF (PRESENT(PLS))  GPRESENT_PLS  = .TRUE.
-IF (PRESENT(PCPH)) GPRESENT_PCPH = .TRUE.
 !
 !$acc kernels
 PCLDFR(:,:,:) = 0. ! Initialize values
@@ -437,8 +447,9 @@ IF ( .NOT. OSIGMAS ) THEN
               - ZLV(KIB:KIE,KJB:KJE,IKTB:IKTE)*PRC(KIB:KIE,KJB:KJE,IKTB:IKTE)/ZCPD(KIB:KIE,KJB:KJE,IKTB:IKTE) &
               - ZLS(KIB:KIE,KJB:KJE,IKTB:IKTE)*PRI(KIB:KIE,KJB:KJE,IKTB:IKTE)/ZCPD(KIB:KIE,KJB:KJE,IKTB:IKTE)
   ! Determine tropopause/inversion  height from minimum temperature
-  ITPL(KIB:KIE,KJB:KJE)  = KIB+1
+  ITPL(KIB:KIE,KJB:KJE)  = KKB+KKL
   ZTMIN(KIB:KIE,KJB:KJE) = 400.
+  !$acc loop seq
   DO JK = IKTB+1,IKTE-1
      WHERE ( PT(KIB:KIE,KJB:KJE,JK) < ZTMIN(KIB:KIE,KJB:KJE) ) 
           ZTMIN(KIB:KIE,KJB:KJE) = PT(KIB:KIE,KJB:KJE,JK)
@@ -447,6 +458,7 @@ IF ( .NOT. OSIGMAS ) THEN
   END DO
   ! Set the mixing length scale
   ZL(KIB:KIE,KJB:KJE,KKB) = 20.
+  !$acc loop seq
   DO JK = KKB+KKL,KKE,KKL
         ! free troposphere
         ZL(KIB:KIE,KJB:KJE,JK) = ZL0
@@ -483,16 +495,14 @@ IF (OUSERI) CALL COMPUTE_FRAC_ICE(HFRAC_ICE, ZFRAC, PT)
 IF (OUSERI) CALL COMPUTE_FRAC_ICE3D_DEVICE(HFRAC_ICE, ZFRAC, PT)
 #endif
 !
-!acc kernels
-DO JK=IKTB,IKTE
-!PW: note: 10x faster to put the kernels zone inside the JK loop (NVHPC 21.9, NVHPC 22.2) even if indenpendent forced
 !$acc kernels
-  JKP=MAX(MIN(JK+KKL,IKTE),IKTB)
-  JKM=MAX(MIN(JK-KKL,IKTE),IKTB)
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-  DO CONCURRENT (JI=KIB:KIE,JJ=KJB:KJE)
+!$acc_nv loop independent
+DO  CONCURRENT ( JK=IKTB:IKTE )
+  JKPK(JK)=MAX(MIN(JK+KKL,IKTE),IKTB)
+  JKMK(JK)=MAX(MIN(JK-KKL,IKTE),IKTB)
+END DO
+!$acc_nv loop independent collapse(3)
+DO CONCURRENT (JI=KIB:KIE,JJ=KJB:KJE,JK=IKTB:IKTE)
       ! latent heats
       ! saturated water vapor mixing ratio over liquid water
 #ifndef MNH_BITREP
@@ -541,9 +551,9 @@ DO JK=IKTB,IKTE
         END IF
       ELSE
         ! parameterize Sigma_s with first_order closure
-        DZZP(JI,JJ,JK)    =  PZZ(JI,JJ,JKP) - PZZ(JI,JJ,JKM)
-        ZDRWP(JI,JJ,JK)   =  ZRT(JI,JJ,JKP) - ZRT(JI,JJ,JKM)
-        ZDTLP(JI,JJ,JK)   =  ZTLK(JI,JJ,JKP) - ZTLK(JI,JJ,JKM) + XG/ZCPD(JI,JJ,JK) * DZZP(JI,JJ,JK)
+        DZZP(JI,JJ,JK)    =  PZZ(JI,JJ,JKPK(JK)) - PZZ(JI,JJ,JKMK(JK))
+        ZDRWP(JI,JJ,JK)   =  ZRT(JI,JJ,JKPK(JK)) - ZRT(JI,JJ,JKMK(JK))
+        ZDTLP(JI,JJ,JK)   =  ZTLK(JI,JJ,JKPK(JK)) - ZTLK(JI,JJ,JKMK(JK)) + XG/ZCPD(JI,JJ,JK) * DZZP(JI,JJ,JK)
         ZLLP(JI,JJ,JK) = ZL(JI,JJ,JK)
         ! standard deviation due to convection
         ZSIG_CONVP(JI,JJ,JK) =0.
@@ -582,7 +592,7 @@ DO JK=IKTB,IKTE
         PSIGRC(JI,JJ,JK) = PCLDFR(JI,JJ,JK)
 
         !Computation warm/cold Cloud Fraction and content in high water content part
-        IF(PRESENT(PHLC_HCF) .AND. PRESENT(PHLC_HRC))THEN
+        IF(GPRESENT_PHLC_HCF .AND. GPRESENT_PHLC_HRC)THEN
           IF(1-ZFRAC(JI,JJ,JK) > 1.E-20)THEN
              ZAUTCP(JI,JJ,JK) = (ZSBARP(JI,JJ,JK) - XCRIAUTC/(PRHODREF(JI,JJ,JK)*(1-ZFRAC(JI,JJ,JK)))) &
                                 /ZSIGMAP(JI,JJ,JK)
@@ -600,7 +610,7 @@ DO JK=IKTB,IKTE
           ENDIF
         ENDIF
 
-        IF(PRESENT(PHLI_HCF) .AND. PRESENT(PHLI_HRI))THEN
+        IF(GPRESENT_PHLI_HCF .AND. GPRESENT_PHLI_HRI)THEN
           IF(ZFRAC(JI,JJ,JK) > 1.E-20)THEN
             ZCRIAUTIP(JI,JJ,JK)=MIN(XCRIAUTI,10**(XACRIAUTI*(PT(JI,JJ,JK)-XTT)+XBCRIAUTI))
             ZAUTIP(JI,JJ,JK) = (ZSBARP(JI,JJ,JK) - ZCRIAUTIP(JI,JJ,JK)/ZFRAC(JI,JJ,JK))/ZSIGMAP(JI,JJ,JK)
@@ -648,11 +658,11 @@ DO JK=IKTB,IKTE
 
         PSIGRC(JI,JJ,JK) =  MIN(1.,(1.-ZINCP(JI,JJ,JK))*ZSRC_1D(INQ1P(JI,JJ,JK))+ZINCP(JI,JJ,JK)*ZSRC_1D(INQ1P(JI,JJ,JK)+1))
 
-        IF(PRESENT(PHLC_HCF) .AND. PRESENT(PHLC_HRC))THEN
+        IF(GPRESENT_PHLC_HCF .AND. GPRESENT_PHLC_HRC)THEN
           PHLC_HCF(JI,JJ,JK)=0.
           PHLC_HRC(JI,JJ,JK)=0.
         ENDIF
-        IF(PRESENT(PHLI_HCF) .AND. PRESENT(PHLI_HRI))THEN
+        IF(GPRESENT_PHLI_HCF .AND. GPRESENT_PHLI_HRI)THEN
           PHLI_HCF(JI,JJ,JK)=0.
           PHLI_HRI(JI,JJ,JK)=0.
         ENDIF
@@ -689,10 +699,8 @@ DO JK=IKTB,IKTE
         PSIGRC(JI,JJ,JK) = PSIGRC(JI,JJ,JK)* MIN( 3. , MAX(1.,1.-ZQ1P(JI,JJ,JK)) )
       ENDIF
 
-   END DO
-!$acc end kernels
 END DO
-!acc end kernels
+!$acc end kernels
 
 !$acc end data
 
@@ -701,6 +709,7 @@ deallocate( ztlk, zrt, zl, zfrac, itpl, ztmin, zlv, zls,zcpd, gwork )
 deallocate( ZZZP, JKPP, ZPVP, ZQSLP, ZPIVP, ZQSIP, ZLVSP, ZAHP, ZAP, ZBP, ZSBARP, ZSIGMAP )
 deallocate( DZZP, ZDRWP, ZDTLP, ZLLP, ZSIG_CONVP, ZQ1P, ZGCONDP, ZGAUVP, ZCONDP, ZAUTCP )
 deallocate( ZGAUTCP, ZGAUCP, ZCRIAUTIP, ZAUTIP, ZGAUTIP, ZGAUIP, INQ1P, ZINCP, ZRCOLDP, ZRIOLDP )
+deallocate( JKPK,JKMK )
 #else
 !Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
 CALL MNH_MEM_RELEASE()