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

Philippe 21/11/2019: OpenACC: bl89: add OpenACC directives + bit reproducibility with MNH_BITREP

parent 071a5d29
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 82 additions and 9 deletions
src/MNH/bl89.f90
+ 82
− 9
View file @ 49ae1e0b
......@@ -135,12 +135,16 @@ INTEGER :: IIU,IJU
INTEGER :: J1D ! horizontal loop counter
INTEGER :: JK,JKK,J3RD ! loop counters
INTEGER :: JRR ! moist loop counter
#ifdef MNH_OPENACC
integer :: ji, jj
#endif
REAL :: ZRVORD ! Rv/Rd
REAL :: ZPOTE,ZLWORK1,ZLWORK2
REAL :: ZTEST,ZTEST0,ZTESTM ! test for vectorization
REAL :: Z2SQRT2
!-------------------------------------------------------------------------------
!
!$acc data present( pzz, pdzz, pthvref, pthlm, prm, ptkem, pshear, plm )
if ( mppdb_initialized ) then
!Check all in arrays
......@@ -171,6 +175,11 @@ allocate( zrm (size( prm, 1 ) * size( prm, 2 ), size( prm, 3 ), size
if ( krr > 0 ) &
allocate( zsum (size( prm, 1 ) * size( prm, 2 ), size( prm, 3 ) ) )
!$acc data create ( zvpt, zdeltvpt, zhlvpt, zlwork, zinte, &
!$acc & zzz, zdzz, zg_o_thvref, zthm, ztkem, zlm, zlmdn, &
!$acc & zshear, zsqrt_tke, zrm, zsum )
!$acc kernels
Z2SQRT2=2.*SQRT(2.)
IIU=SIZE(PTKEM,1)
IJU=SIZE(PTKEM,2)
......@@ -202,6 +211,7 @@ IF (CPROGRAM=='AROME ') THEN
END DO
END DO
ELSE
#ifndef MNH_OPENACC
DO JK=1,IKT
ZZZ (:,JK) = RESHAPE(PZZ (:,:,JK),(/ IIU*IJU /) )
ZDZZ (:,JK) = RESHAPE(PDZZ (:,:,JK),(/ IIU*IJU /) )
......@@ -213,15 +223,45 @@ ELSE
ZRM (:,JK,JRR) = RESHAPE(PRM (:,:,JK,JRR),(/ IIU*IJU /) )
END DO
END DO
#else
!$acc loop independent collapse(3)
do jk = 1, ikt
do jj = 1, iju
do ji = 1, iiu
zzz (( jj - 1 ) * iiu + ji, jk ) = pzz (ji, jj, jk)
zdzz (( jj - 1 ) * iiu + ji, jk ) = pdzz (ji, jj, jk)
zthm (( jj - 1 ) * iiu + ji, jk ) = pthlm (ji, jj, jk)
zshear (( jj - 1 ) * iiu + ji, jk ) = pshear (ji, jj, jk)
ztkem (( jj - 1 ) * iiu + ji, jk ) = ptkem (ji, jj, jk)
zg_o_thvref(( jj - 1 ) * iiu + ji, jk ) = xg / pthvref(ji, jj, jk)
end do
end do
end do
!$acc loop independent collapse(4)
do jrr = 1, krr
do jk = 1, ikt
do jj = 1, iju
do ji = 1, iiu
zrm(( jj - 1 ) * iiu + ji, jk, jrr ) = prm(ji, jj, jk, jrr )
end do
end do
end do
end do
#endif
END IF
!
#ifndef MNH_BITREP
ZSQRT_TKE = SQRT(ZTKEM)
#else
zsqrt_tke(:, : ) = Br_pow( ztkem, 0.5 )
#endif
!-------------------------------------------------------------------------------
!
!* 2. Virtual potential temperature on the model grid
! -----------------------------------------------
!
IF(KRR /= 0) THEN
IF( KRR > 0 ) THEN
ZSUM(:,:) = 0.
DO JRR=1,KRR
ZSUM(:,:) = ZSUM(:,:)+ZRM(:,:,JRR)
......@@ -246,9 +286,17 @@ END IF
ZDELTVPT(:,IKTB:IKTE)=ZVPT(:,IKTB:IKTE)-ZVPT(:,IKTB-KKL:IKTE-KKL)
ZDELTVPT(:,KKU)=ZVPT(:,KKU)-ZVPT(:,KKU-KKL)
ZDELTVPT(:,KKA)=0.
#ifndef MNH_OPENACC
WHERE (ABS(ZDELTVPT(:,:))<XLINF)
ZDELTVPT(:,:)=XLINF
END WHERE
#else
do jk = 1, ikt
do ji = 1, iiu * iju
if ( abs( zdeltvpt(ji, jk ) ) < xlinf ) zdeltvpt(ji, jk ) = xlinf
end do
end do
#endif
!
ZHLVPT(:,IKTB:IKTE)= 0.5 * ( ZVPT(:,IKTB:IKTE)+ZVPT(:,IKTB-KKL:IKTE-KKL) )
ZHLVPT(:,KKU)= 0.5 * ( ZVPT(:,KKU)+ZVPT(:,KKU-KKL) )
......@@ -266,7 +314,7 @@ DO JK=IKTB,IKTE
!* 4. mixing length for a downwards displacement
! ------------------------------------------
ZINTE(:)=ZTKEM(:,JK)
ZLWORK=0.
ZLWORK(:)=0.
ZTESTM=1.
DO JKK=JK,IKB,-KKL
IF(ZTESTM > 0.) THEN
......@@ -327,7 +375,7 @@ DO JK=IKTB,IKTE
! -----------------------------------------
!
ZINTE(:)=ZTKEM(:,JK)
ZLWORK=0.
ZLWORK(:)=0.
ZTESTM=1.
!
DO JKK=JK+KKL,IKE,KKL
......@@ -360,19 +408,25 @@ DO JK=IKTB,IKTE
! ( ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK) )
!--------- SHEAR + STABILITY -----------
#ifndef MNH_BITREP
ZLWORK2= ( - ZG_O_THVREF(J1D,JK) *(ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK) ) &
- XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
+ SQRT (ABS( &
#ifndef MNH_BITREP
(XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
+ ( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK))) )**2 &
+ 2. * ZINTE(J1D) * &
( ZG_O_THVREF(J1D,JK)* ZDELTVPT(J1D,JKK)/ZDZZ(J1D,JKK))))) / &
(ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK))
#else
ZLWORK2= ( - ZG_O_THVREF(J1D,JK) *(ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK) ) &
- XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
+ BR_POW (ABS( &
BR_P2(XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
+ ( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK))) ) &
#endif
+ 2. * ZINTE(J1D) * &
( ZG_O_THVREF(J1D,JK)* ZDELTVPT(J1D,JKK)/ZDZZ(J1D,JKK))))) / &
( ZG_O_THVREF(J1D,JK)* ZDELTVPT(J1D,JKK)/ZDZZ(J1D,JKK))), 0.5 ) ) / &
(ZG_O_THVREF(J1D,JK) * ZDELTVPT(J1D,JKK) / ZDZZ(J1D,JKK))
#endif
ZLWORK(J1D)=ZLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2)
ZINTE(J1D) = ZINTE(J1D) - ZPOTE
......@@ -388,9 +442,14 @@ DO JK=IKTB,IKTE
ZLWORK1=MAX(ZLMDN(J1D,JK),1.E-10_MNHREAL)
ZLWORK2=MAX(ZLWORK(J1D),1.E-10_MNHREAL)
ZPOTE = ZLWORK1 / ZLWORK2
#ifndef MNH_BITREP
ZLWORK2=1.d0 + ZPOTE**(2./3.)
ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2))
END DO
#else
ZLWORK2=1.d0 + br_pow(ZPOTE,2./3.)
ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*br_pow(ZLWORK2,0.5))
#endif
END DO
ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)
......@@ -420,15 +479,29 @@ IF (CPROGRAM=='AROME ') THEN
PLM (:,1,JK) = ZLM (:,JK)
END DO
ELSE
#ifndef MNH_OPENACC
DO JK=1,IKT
PLM (:,:,JK) = RESHAPE(ZLM (:,JK), (/ IIU,IJU /) )
END DO
#else
do jk = 1, ikt
do jj = 1, iju
do ji = 1, iiu
plm(ji, jj, jk ) = zlm(( jj - 1 ) * iiu + ji, jk )
end do
end do
end do
#endif
END IF
!$acc end kernels
if ( mppdb_initialized ) then
!Check all out arrays
call Mppdb_check( plm, "Bl89 end:plm" )
end if
!
!$acc end data
!$acc end data
END SUBROUTINE BL89
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