bl89.f90

!MNH_LIC Copyright 1997-2021 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

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 

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
!
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,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

allocate( zvpt       (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) )
allocate( zdeltvpt   (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) )
allocate( zhlvpt     (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) )
allocate( zlwork     (size( ptkem, 1 ) * size( ptkem, 2 ) ) )
allocate( zinte      (size( ptkem, 1 ) * size( ptkem, 2 ) ) )
allocate( zzz        (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zdzz       (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zg_o_thvref(size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zthm       (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          ztkem      (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zlm        (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zlmdn      (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zshear     (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) ) , &
          zsqrt_tke  (size( ptkem, 1 ) * size( ptkem, 2 ), size( ptkem, 3 ) )   )
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 ) ) )

!$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)
!
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 (CPROGRAM=='AROME ') THEN
  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
  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
zsqrt_tke(:, : ) = Br_pow( ztkem, 0.5 )
!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
!-------------------------------------------------------------------------------
!
!*       2.    Virtual potential temperature on the model grid
!              -----------------------------------------------
!
IF( KRR > 0 ) THEN
  ZSUM(:,:) = 0.
  DO JRR=1,KRR
    ZSUM(:,:) = ZSUM(:,:)+ZRM(:,:,JRR)
  ENDDO
  ZVPT(:,1:)=ZTHM(:,:) * ( 1. + ZRVORD*ZRM(:,:,1) )  &
                           / ( 1. + ZSUM(:,:) )
ELSE
  ZVPT(:,1:)=ZTHM(:,:)
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.
!!!!!!!!!!!!

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 (CPROGRAM=='AROME ') 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

!$acc end data

END SUBROUTINE BL89