Forked from
Méso-NH / Méso-NH code
2183 commits behind, 1063 commits ahead of the upstream repository.
-
ESCOBAR MUNOZ Juan authoredJuan&Naima 16/03/2023:MNH/advec_weno_k_3_aux.f90,bl89.f90,ice4_sedimentation_split.f90, Cray OpenACC opt , add present_cr to avoid recurrence in kernels
ESCOBAR MUNOZ Juan authoredJuan&Naima 16/03/2023:MNH/advec_weno_k_3_aux.f90,bl89.f90,ice4_sedimentation_split.f90, Cray OpenACC opt , add present_cr to avoid recurrence in kernels
bl89.f90 20.39 KiB
!MNH_LIC Copyright 1997-2022 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! ################
MODULE MODI_BL89
! ################
INTERFACE
SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM)
!
INTEGER, INTENT(IN) :: KKA
INTEGER, INTENT(IN) :: KKU
INTEGER, INTENT(IN) :: KKL
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
INTEGER, INTENT(IN) :: KRR
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSHEAR
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM
END SUBROUTINE BL89
END INTERFACE
END MODULE MODI_BL89
!
! #########################################################
SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM)
! #########################################################
!
!!**** *BL89* -
!!
!! PURPOSE
!! -------
!! This routine computes the mixing length from Bougeault-Lacarrere 89
!! formula.
!!
!!** METHOD
!! ------
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!!
!! REFERENCE
!! ---------
!!
!! Book 2
!!
!! AUTHOR
!! ------
!!
!! J. Cuxart INM and Meteo-France
!!
!! MODIFICATIONS
!! -------------
!! Original 27/04/97 (V. Masson) separation from turb.f90
!! and optimization
!! 06/01/98 (V. Masson and P. Jabouille) optimization
!! 15/03/99 (V. Masson) new lup ldown averaging
!! 21/02/01 (P. Jabouille) improve vectorization
!! 2012-02 (Y. Seity) add possibility to run with
!! reversed vertical levels
!! Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
!! 01/2019 (Q. Rodier) support for RM17 mixing length!! 03/2021 (JL Redelsperger) Ocean model case
!! 06/2021 (P. Marquet) correction of exponent on final length according to Lemarié et al. 2021
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CONF, ONLY: CPROGRAM
USE MODD_CST
USE MODD_CTURB
USE MODD_DYN_n, ONLY: LOCEAN
USE MODD_PARAMETERS
use modd_precision, only: MNHREAL
#ifdef MNH_OPENACC
USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
#endif
use mode_mppdb
#ifdef MNH_BITREP
USE MODI_BITREP
#endif
!
!
IMPLICIT NONE
!
!* 0.1 Declaration of arguments
! ------------------------
!
INTEGER, INTENT(IN) :: KKA !near ground array index
INTEGER, INTENT(IN) :: KKU !uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL !vert. levels type 1=MNH -1=ARO
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM ! conservative pot. temp.
INTEGER, INTENT(IN) :: KRR
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! water var.
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSHEAR
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM ! Mixing length
! thermodynamical variables PTHLM=Theta at the begining
!
!* 0.2 Declaration of local variables
! ------------------------------
!
INTEGER :: IKB,IKE
INTEGER :: IKT ! array size in k direction
INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
#ifndef MNH_OPENACC
real, dimension(:,:), allocatable :: ZVPT ! Virtual Potential Temp at half levels
real, dimension(:,:), allocatable :: ZDELTVPT
! Increment of Virtual Potential Temp between two following levels
real, dimension(:,:), allocatable :: ZHLVPT
! Virtual Potential Temp at half levels
real, dimension(:), allocatable :: ZLWORK,ZINTE
! ! downwards then upwards vertical displacement,
! ! residual internal energy,
! ! residual potential energy
real, dimension(:,:), allocatable :: ZZZ,ZDZZ, &
ZG_O_THVREF, &
ZTHM,ZTKEM,ZLM, &
ZLMDN,ZSHEAR, &
ZSQRT_TKE
! ! input and output arrays packed according one horizontal coord.
real, dimension(:,:,:), allocatable :: ZRM
! ! input array packed according one horizontal coord.
real, dimension(:,:), allocatable :: ZSUM ! to replace SUM function
#elsereal, dimension(:,:), pointer, contiguous :: ZVPT ! Virtual Potential Temp at half levels
real, dimension(:,:), pointer, contiguous :: ZDELTVPT
! Increment of Virtual Potential Temp between two following levels
real, dimension(:,:), pointer, contiguous :: ZHLVPT
! Virtual Potential Temp at half levels
real, dimension(:), pointer, contiguous :: ZLWORK,ZINTE
! ! downwards then upwards vertical displacement,
! ! residual internal energy,
! ! residual potential energy
real, dimension(:,:), pointer, contiguous :: ZZZ,ZDZZ, &
ZG_O_THVREF, &
ZTHM,ZTKEM,ZLM, &
ZLMDN,ZSHEAR, &
ZSQRT_TKE
! ! input and output arrays packed according one horizontal coord.
real, dimension(:,:,:), pointer, contiguous :: ZRM
! ! input array packed according one horizontal coord.
real, dimension(:,:), pointer, contiguous :: ZSUM ! to replace SUM function
#endif
!
INTEGER :: IIU, IJU
INTEGER :: J1D ! horizontal loop counter
INTEGER :: JK,JKK,J3RD ! loop counters
INTEGER :: JRR ! moist loop counter
integer :: ji, jj
LOGICAL :: LAROME
REAL :: ZRVORD ! Rv/Rd
REAL :: ZPOTE,ZLWORK1,ZLWORK2
REAL :: ZTEST,ZTEST0,ZTESTM ! test for vectorization
REAL :: Z2SQRT2,ZUSRBL89,ZBL89EXP
!-------------------------------------------------------------------------------
!$acc data present( pzz, pdzz, pthvref, pthlm, prm, ptkem, pshear, plm )
if ( mppdb_initialized ) then
!Check all in arrays
call Mppdb_check( pzz, "Bl89 beg:pzz" )
call Mppdb_check( pdzz, "Bl89 beg:pdzz" )
call Mppdb_check( pthvref, "Bl89 beg:pthvref" )
call Mppdb_check( pthlm, "Bl89 beg:pthlm" )
call Mppdb_check( prm, "Bl89 beg:prm" )
call Mppdb_check( ptkem, "Bl89 beg:ptkem" )
call Mppdb_check( pshear, "Bl89 beg:pshear" )
end if
IF ( CPROGRAM == 'AROME' ) THEN
LAROME = .TRUE.
ELSE
LAROME = .FALSE.
END IF
IIU = SIZE( PTKEM, 1 )
IJU = SIZE( PTKEM, 2 )
IKT = SIZE( PTKEM, 3 )
#ifndef MNH_OPENACC
allocate( zvpt (IIU * IJU, IKT ) )
allocate( zdeltvpt (IIU * IJU, IKT ) )
allocate( zhlvpt (IIU * IJU, IKT ) )
allocate( zlwork (IIU * IJU ) )
allocate( zinte (IIU * IJU ) )
allocate( zzz (IIU * IJU, IKT ) , &
zdzz (IIU * IJU, IKT ) , &
zg_o_thvref(IIU * IJU, IKT ) , &
zthm (IIU * IJU, IKT ) , &
ztkem (IIU * IJU, IKT ) , &
zlm (IIU * IJU, IKT ) , &
zlmdn (IIU * IJU, IKT ) , &
zshear (IIU * IJU, IKT ) , &
zsqrt_tke (IIU * IJU, IKT ) )allocate( zrm (size( prm, 1 ) * size( prm, 2 ), size( prm, 3 ), size( prm, 4 ) ) )
if ( krr > 0 ) &
allocate( zsum (size( prm, 1 ) * size( prm, 2 ), size( prm, 3 ) ) )
#else
!Pin positions in the pools of MNH memory
CALL MNH_MEM_POSITION_PIN()
CALL MNH_MEM_GET( zvpt, IIU * IJU, IKT )
CALL MNH_MEM_GET( zdeltvpt, IIU * IJU, IKT )
CALL MNH_MEM_GET( zhlvpt, IIU * IJU, IKT )
CALL MNH_MEM_GET( zlwork, IIU * IJU )
CALL MNH_MEM_GET( zinte, IIU * IJU )
CALL MNH_MEM_GET( zzz, IIU * IJU, IKT )
CALL MNH_MEM_GET( zdzz, IIU * IJU, IKT )
CALL MNH_MEM_GET( zg_o_thvref, IIU * IJU, IKT )
CALL MNH_MEM_GET( zthm, IIU * IJU, IKT )
CALL MNH_MEM_GET( ztkem, IIU * IJU, IKT )
CALL MNH_MEM_GET( zlm, IIU * IJU, IKT )
CALL MNH_MEM_GET( zlmdn, IIU * IJU, IKT )
CALL MNH_MEM_GET( zshear, IIU * IJU, IKT )
CALL MNH_MEM_GET( zsqrt_tke, IIU * IJU, IKT )
CALL MNH_MEM_GET( zrm, size( prm, 1 ) * size( prm, 2 ), size( prm, 3 ), size( prm, 4 ) )
if ( krr > 0 ) &
CALL MNH_MEM_GET( zsum, size( prm, 1 ) * size( prm, 2 ), size( prm, 3 ) )
!$acc data present ( zvpt, zdeltvpt, zhlvpt, zlwork, zinte, &
!$acc & zzz, zdzz, zg_o_thvref, zthm, ztkem, zlm, zlmdn, &
!$acc & zshear, zsqrt_tke, zrm, zsum )
#endif
!$acc kernels
Z2SQRT2=2.*SQRT(2.)
!
IKB=KKA+JPVEXT_TURB*KKL
IKE=KKU-JPVEXT_TURB*KKL
IKTB = JPVEXT_TURB + 1
IKT = SIZE(PTKEM,3)
IKTE = IKT-JPVEXT_TURB
ZRVORD = XRV / XRD
!
!-------------------------------------------------------------------------------
!
!* 1. pack the horizontal dimensions into one
! ---------------------------------------
!
IF ( LAROME ) THEN
!$acc loop independent
DO JK=1,IKT
ZZZ (:,JK) = PZZ (:,1,JK)
ZDZZ (:,JK) = PDZZ (:,1,JK)
ZTHM (:,JK) = PTHLM (:,1,JK)
ZTKEM (:,JK) = PTKEM (:,1,JK)
ZG_O_THVREF(:,JK) = XG/PTHVREF(:,1,JK)
END DO
!$acc loop independent collapse(2)
DO JK=1,IKT
DO JRR=1,KRR
ZRM (:,JK,JRR) = PRM (:,1,JK,JRR)
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 /) )
ZTHM (:,JK) = RESHAPE(PTHLM (:,:,JK),(/ IIU*IJU /) )
ZSHEAR (:,JK) = RESHAPE(PSHEAR (:,:,JK),(/ IIU*IJU /) )
ZTKEM (:,JK) = RESHAPE(PTKEM (:,:,JK),(/ IIU*IJU /) )
ZG_O_THVREF(:,JK) = RESHAPE(XG/PTHVREF(:,:,JK),(/ IIU*IJU /) ) IF (LOCEAN) ZG_O_THVREF(:,JK) = XG * XALPHAOC
DO JRR=1,KRR
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)
!ZBL89EXP is defined here because (and not in ini_cturb) because XCED is defined in read_exseg (depending on BL89/RM17)
ZBL89EXP = LOG(16.)/(4.*LOG(XKARMAN)+LOG(XCED)-3.*LOG(XCMFS))
#else
#ifdef MNH_COMPILER_NVHPC
!$acc loop independent collapse(2)
#endif
do concurrent( ji = 1 : iiu, jj = 1 : iju )
zsqrt_tke(ji, jj) = Br_pow( ztkem(ji,jj), 0.5 )
end do
!ZBL89EXP is defined here because (and not in ini_cturb) because XCED is defined in read_exseg (depending on BL89/RM17)
ZBL89EXP = Br_log( 16. ) / ( 4. * Br_log( XKARMAN )+ Br_log( XCED ) - 3. * Br_log( XCMFS) )
#endif
ZUSRBL89 = 1./ZBL89EXP
!$acc end kernels
!-------------------------------------------------------------------------------
!
!* 2. Virtual potential temperature on the model grid
! -----------------------------------------------
!
IF( KRR > 0 ) THEN
!$acc kernels
ZSUM(:,:) = 0.
DO JRR=1,KRR
ZSUM(:,:) = ZSUM(:,:)+ZRM(:,:,JRR)
ENDDO
ZVPT(:,:)=ZTHM(:,:) * ( 1. + ZRVORD*ZRM(:,:,1) ) &
/ ( 1. + ZSUM(:,:) )
!$acc end kernels
ELSE
!$acc kernels
ZVPT(:,:)=ZTHM(:,:)
!$acc end kernels
END IF
!
!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!! WARNING !!
!!!!!!!!!!!!
!!!!!!!!!!!!
!Any modification done to the following lines and to the sections 4 and
!6 must be copied in compute_bl89_ml routine.
!We do not call directly this routine for numerical performance reasons
!but algorithm must remain the same.
!!!!!!!!!!!!
!$acc kernels present_cr(ZTEST0,ZPOTE,ZTEST,ZTESTM,ZLWORK1,ZLWORK2,ZLWORK,ZINTE) &
!$acc & present_cr(ZTEST0,ZPOTE,ZTEST,ZTESTM,ZLWORK1,ZLWORK2,ZLWORK,ZINTE) &
!$acc & present_cr(ZLWORK1,ZLWORK2,ZPOTE,ZLWORK2,ZLM)
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) )
ZHLVPT(:,KKA) = ZVPT(:,KKA)
!-------------------------------------------------------------------------------
!
!* 3. loop on model levels
! --------------------
DO JK=IKTB,IKTE
!
!-------------------------------------------------------------------------------
!
!
!* 4. mixing length for a downwards displacement
! ------------------------------------------
ZINTE(:)=ZTKEM(:,JK)
ZLWORK(:)=0.
ZTESTM=1.
DO JKK=JK,IKB,-KKL
IF(ZTESTM > 0.) THEN
ZTESTM=0.
DO J1D=1,IIU*IJU
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
!--------- SHEAR + STABILITY -----------
ZPOTE = ZTEST0* &
(-ZG_O_THVREF(J1D,JK)*(ZHLVPT(J1D,JKK)-ZVPT(J1D,JK)) &
+ XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
)*ZDZZ(J1D,JKK)
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE)
ZTESTM=ZTESTM+ZTEST0
ZLWORK1=ZDZZ(J1D,JKK)
!-------- ORIGINAL -------------
! ZLWORK2= ( + ZG_O_THVREF(J1D,JK) * &
! ( ZVPT(J1D,JKK) - ZVPT(J1D,JK) ) &
! + SQRT (ABS( &
! ( ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK) - 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))
!--------- SHEAR + STABILITY -----------
ZLWORK2 = (ZG_O_THVREF(J1D,JK) *(ZVPT(J1D,JKK) - ZVPT(J1D,JK)) &
-XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
#ifndef MNH_BITREP
+ SQRT(ABS( (XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
+ ( -ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK) - ZVPT(J1D,JK)) ))**2.0 + &
#else
+ SQRT(ABS( BR_P2(XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
+ ( -ZG_O_THVREF(J1D,JK) * (ZVPT(J1D,JKK) - 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))
ZLWORK(J1D)=ZLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2)
ZINTE(J1D) = ZINTE(J1D) - ZPOTE
END DO
ENDIF
END DO
!-------------------------------------------------------------------------------
!
!* 5. intermediate storage of the final mixing length
! -----------------------------------------------
!
ZLMDN(:,JK)=MIN(ZLWORK(:),0.5*(ZZZ(:,JK)+ZZZ(:,JK+KKL))-ZZZ(:,IKB))
!
!-------------------------------------------------------------------------------
!
!* 6. mixing length for an upwards displacement
! -----------------------------------------
!
ZINTE(:)=ZTKEM(:,JK)
ZLWORK(:)=0.
ZTESTM=1.
!
DO JKK=JK+KKL,IKE,KKL
IF(ZTESTM > 0.) THEN
ZTESTM=0.
DO J1D=1,IIU*IJU
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
!-------- ORIGINAL -------------
!ZPOTE = ZTEST0*ZG_O_THVREF(J1D,JK) * &
! (ZHLVPT(J1D,JKK) - ZVPT(J1D,JK) ) *ZDZZ(J1D,JKK)
!--------- SHEAR + STABILITY -----------
ZPOTE = ZTEST0* &
(ZG_O_THVREF(J1D,JK)*(ZHLVPT(J1D,JKK)-ZVPT(J1D,JK)) &
+XRM17*ZSHEAR(J1D,JKK)*ZSQRT_TKE(J1D,JK) &
)*ZDZZ(J1D,JKK)
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE)
ZTESTM=ZTESTM+ZTEST0
ZLWORK1=ZDZZ(J1D,JKK)
!-------- ORIGINAL -------------
! ZLWORK2= ( - ZG_O_THVREF(J1D,JK) * &
! ( ZVPT(J1D,JKK-KKL) - ZVPT(J1D,JK) ) &
! + SQRT (ABS( &
! ( 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) )
!--------- 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( &
(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))) ) &
+ 2. * ZINTE(J1D) * &
( 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
END DO
ENDIF
END DO
!
!-------------------------------------------------------------------------------
!
!* 7. final mixing length
!
DO J1D=1,IIU*IJU
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**ZBL89EXP
ZLM(J1D,JK) = ZLWORK1*(2./ZLWORK2)**ZUSRBL89
#else
ZLWORK2=1.d0 + br_pow( ZPOTE, ZBL89EXP )
ZLM(J1D,JK) = ZLWORK1 * br_pow( 2. / ZLWORK2, ZUSRBL89 )
#endif
END DO
ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)
!
!
!* 8. end of the loop on the vertical levels
! --------------------------------------
!
END DO
!
!-------------------------------------------------------------------------------
!
!* 9. boundaries
! ----------
!
ZLM(:,KKA)=ZLM(:,IKB)
ZLM(:,IKE)=ZLM(:,IKE-KKL)
ZLM(:,KKU)=ZLM(:,IKE-KKL)
!
!-------------------------------------------------------------------------------
!
!* 10. retrieve output array in model coordinates
! ------------------------------------------
!
IF ( LAROME ) THEN
DO JK=1,IKT
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
#ifdef MNH_OPENACC
!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
CALL MNH_MEM_RELEASE()
#endif
!$acc end data
END SUBROUTINE BL89