Skip to content
Snippets Groups Projects
Commit 71e0d56f authored by WAUTELET Philippe's avatar WAUTELET Philippe
Browse files

Philippe 01/08/2019: OpenACC: implementation of slow_terms

parent 6fa32f36
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 175 additions and 53 deletions
src/MNH/slow_terms.f90
+ 175
− 53
View file @ 71e0d56f
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme