diff --git a/src/MNH/bl89.f90 b/src/MNH/bl89.f90
index 2cb05385b235a5d3ce7ec1c420a2c5a8070958e8..e28f1cc0e73223cb604dffbf9368770284d44a53 100644
--- a/src/MNH/bl89.f90
+++ b/src/MNH/bl89.f90
@@ -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