bl89.f90

!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
#else
real, 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( 'BL89' )

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 * iju, jk = 1 : ikt )
  zsqrt_tke(ji, jk) = Br_pow( ztkem(ji,jk), 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(:,:) = ZRM(:,:,1)
  DO JRR = 2, KRR
    !$mnh_do_concurrent ( JI=1:IIU,JJ=1:IJU )
    ZSUM(JI,JJ) = ZSUM(JI,JJ) + ZRM(JI,JJ,JRR)
    !$mnh_end_do()
  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)
!$acc end kernels
!-------------------------------------------------------------------------------
!
!*       3.  loop on model levels
!            --------------------
DO JK=IKTB,IKTE
!$acc kernels
!
!-------------------------------------------------------------------------------
!
!

!*       4.  mixing length for a downwards displacement
!            ------------------------------------------
  ZINTE(:)=ZTKEM(:,JK)
  ZLWORK(:)=0.
  ZTESTM=1.
!$acc loop seq
  DO JKK=JK,IKB,-KKL  
    IF(ZTESTM > 0.) THEN
      ZTESTM=0.
      !$mnh_do_concurrent( 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
      !$mnh_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.
!
!$acc loop seq
  DO JKK=JK+KKL,IKE,KKL 
    IF(ZTESTM > 0.) THEN
      ZTESTM=0.
      !$mnh_do_concurrent( 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
      !$mnh_end_do()
    ENDIF
  END DO
!
!-------------------------------------------------------------------------------
!
!*       7.  final mixing length
!
  !$mnh_do_concurrent( 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
  !$mnh_end_do()

  ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)
!
!
!*       8.  end of the loop on the vertical levels
!            --------------------------------------
!
!$acc end kernels
END DO
!
!-------------------------------------------------------------------------------
!
!*       9.  boundaries
!            ----------
!
!$acc kernels
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( 'BL89')
#endif

!$acc end data

END SUBROUTINE BL89