diff --git a/src/common/turb/bl89.f90 b/src/common/turb/bl89.f90
new file mode 100644
index 0000000000000000000000000000000000000000..303d92d262407a571fbe894e1da0274d1626abf3
--- /dev/null
+++ b/src/common/turb/bl89.f90
@@ -0,0 +1,335 @@
+! ######spl
+ SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,PLM)
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! #########################################################
+!
+!!**** *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
+!
+!
+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(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZVPT ! Virtual Potential Temp at half levels
+REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZDELTVPT
+ ! Increment of Virtual Potential Temp between two following levels
+REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZHLVPT
+ ! Virtual Potential Temp at half levels
+REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2)) :: ZLWORK,ZINTE
+! ! downwards then upwards vertical displacement,
+! ! residual internal energy,
+! ! residual potential energy
+REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZZZ,ZDZZ, &
+ ZG_O_THVREF, &
+ ZTHM,ZTKEM,ZLM, &
+ ZLMDN,ZSHEAR, &
+ ZSQRT_TKE
+! ! input and output arrays packed according one horizontal coord.
+REAL, DIMENSION(SIZE(PRM,1)*SIZE(PRM,2),SIZE(PRM,3),SIZE(PRM,4)) :: ZRM
+! ! input array packed according one horizontal coord.
+REAL, DIMENSION(SIZE(PRM,1)*SIZE(PRM,2),SIZE(PRM,3)) :: ZSUM ! to replace SUM function
+!
+INTEGER :: IIU,IJU
+INTEGER :: J1D ! horizontal loop counter
+INTEGER :: JK,JKK ! loop counters
+INTEGER :: JRR ! moist loop counter
+REAL :: ZRVORD ! Rv/Rd
+REAL :: ZPOTE,ZLWORK1,ZLWORK2
+REAL :: ZTEST,ZTEST0,ZTESTM ! test for vectorization
+REAL :: Z2SQRT2,ZUSRBL89,ZBL89EXP
+!-------------------------------------------------------------------------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('BL89',0,ZHOOK_HANDLE)
+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
+ 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
+END IF
+!
+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))
+ZUSRBL89 = 1./ZBL89EXP
+!-------------------------------------------------------------------------------
+!
+!* 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.
+WHERE (ABS(ZDELTVPT(:,:))<XLINF)
+ ZDELTVPT(:,:)=XLINF
+END WHERE
+!
+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)
+
+ !--------- SHEAR + STABILITY -----------
+ ZLWORK2 = (ZG_O_THVREF(J1D,JK) *(ZVPT(J1D,JKK) - 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) - ZVPT(J1D,JK)) ))**2.0 + &
+ 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
+ !--------- 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)
+ !--------- SHEAR + STABILITY -----------
+ 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))
+ 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
+ ZLWORK2=1.d0 + ZPOTE**ZBL89EXP
+ ZLM(J1D,JK) = ZLWORK1*(2./ZLWORK2)**ZUSRBL89
+ 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
+ DO JK=1,IKT
+ PLM (:,:,JK) = RESHAPE(ZLM (:,JK), (/ IIU,IJU /) )
+ END DO
+END IF
+
+!
+IF (LHOOK) CALL DR_HOOK('BL89',1,ZHOOK_HANDLE)
+END SUBROUTINE BL89
diff --git a/src/common/turb/compute_bl89_ml.f90 b/src/common/turb/compute_bl89_ml.f90
new file mode 100644
index 0000000000000000000000000000000000000000..5a75011c8ed80ee7d7a3b364136a9ecff36168ea
--- /dev/null
+++ b/src/common/turb/compute_bl89_ml.f90
@@ -0,0 +1,216 @@
+! ######spl
+ SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
+ PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PSHEAR,PLWORK)
+
+ USE PARKIND1, ONLY : JPRB
+ USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+! ###################################################################
+!!
+!! COMPUTE_BL89_ML routine to:
+!! 1/ compute upward or downward mixing length with BL89 formulation
+!!
+!! AUTHOR
+!! ------
+!! J. PERGAUD
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 19/01/06
+!! S. Riette Jan 2012: support for both order of vertical levels and cleaning
+!! R.Honnert Oct 2016 : Update with AROME
+!! Q.Rodier 01/2019 : support RM17 mixing length as in bl89.f90
+! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+! ------------
+!
+!!!!!!!!!!!!
+!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!! WARNING !!
+!!!!!!!!!!!!
+!!!!!!!!!!!!
+!Any modification done to this routine must be copied in bl89.f90.
+!This routine was inlined in bl89 for numerical performance reasons
+!but algorithm must remain the same.
+!!!!!!!!!!!!
+!
+USE MODD_CTURB
+USE MODD_PARAMETERS, ONLY: JPVEXT
+!
+USE MODE_MSG
+!
+USE MODI_SHUMAN_MF
+!
+IMPLICIT NONE
+!
+! 0.1 arguments
+!
+INTEGER, INTENT(IN) :: KKA ! near ground array index
+INTEGER, INTENT(IN) :: KKB ! near ground physical index
+INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
+INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
+INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ2D ! height difference between two mass levels
+REAL, DIMENSION(:), INTENT(IN) :: PTKEM_DEP ! TKE to consume
+REAL, DIMENSION(:), INTENT(IN) :: PG_O_THVREF ! g/ThetaVRef at the departure point
+REAL, DIMENSION(:,:), INTENT(IN) :: PVPT ! ThetaV on mass levels
+INTEGER, INTENT(IN) :: KK ! index of departure level
+LOGICAL, INTENT(IN) :: OUPORDN ! switch to compute upward (true) or
+ ! downward (false) mixing length
+LOGICAL, INTENT(IN) :: OFLUX ! Computation must be done from flux level
+REAL, DIMENSION(:), INTENT(OUT) :: PLWORK ! Resulting mixing length
+REAL, DIMENSION(:,:), INTENT(IN) :: PSHEAR ! vertical wind shear for RM17 mixing length
+
+! 0.2 Local variable
+!
+REAL, DIMENSION(SIZE(PVPT,1)) :: ZLWORK1,ZLWORK2 ! Temporary mixing length
+REAL, DIMENSION(SIZE(PVPT,1)) :: ZINTE,ZPOTE ! TKE and potential energy
+ ! between 2 levels
+REAL, DIMENSION(SIZE(PVPT,1)) :: ZVPT_DEP ! Thetav on departure point
+!
+REAL, DIMENSION(SIZE(PVPT,1),SIZE(PVPT,2)) :: ZDELTVPT,ZHLVPT
+ !Virtual Potential Temp at Half level and DeltaThv between
+ !2 mass levels
+
+INTEGER :: IIJU !Internal Domain
+INTEGER :: J1D !horizontal loop counter
+INTEGER :: JKK !loop counters
+REAL :: ZTEST,ZTEST0,ZTESTM !test for vectorization
+!-------------------------------------------------------------------------------------
+!
+!* 1. INITIALISATION
+! --------------
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('COMPUTE_BL89_ML',0,ZHOOK_HANDLE)
+IIJU=SIZE(PVPT,1)
+!
+ZDELTVPT(:,:)=DZM_MF(KKA,KKU,KKL,PVPT(:,:))
+ZDELTVPT(:,KKA)=0.
+WHERE (ABS(ZDELTVPT(:,:))<XLINF)
+ ZDELTVPT(:,:)=XLINF
+END WHERE
+!
+ZHLVPT(:,:)=MZM_MF(KKA,KKU,KKL,PVPT(:,:))
+!
+!We consider that gradient between mass levels KKB and KKB+KKL is the same as
+!the gradient between flux level KKB and mass level KKB
+ZDELTVPT(:,KKB)=PDZZ2D(:,KKB)*ZDELTVPT(:,KKB+KKL)/PDZZ2D(:,KKB+KKL)
+ZHLVPT(:,KKB)=PVPT(:,KKB)-ZDELTVPT(:,KKB)*0.5
+!
+!
+!
+!* 2. CALCULATION OF THE UPWARD MIXING LENGTH
+! ---------------------------------------
+!
+
+IF (OUPORDN.EQV..TRUE.) THEN
+ ZINTE(:)=PTKEM_DEP(:)
+ PLWORK=0.
+ ZTESTM=1.
+ IF(OFLUX)THEN
+ ZVPT_DEP(:)=ZHLVPT(:,KK) ! departure point is on flux level
+ !We must compute what happens between flux level KK and mass level KK
+ DO J1D=1,IIJU
+ ZTEST0=0.5+SIGN(0.5,ZINTE(J1D)) ! test if there's energy to consume
+ ! Energy consumed if parcel cross the entire layer
+ ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
+ (0.5*(ZHLVPT(J1D,KK)+ PVPT(J1D,KK)) - ZVPT_DEP(J1D)) + &
+ XRM17*PSHEAR(J1D,KK)*SQRT(ABS(PTKEM_DEP(J1D)))) * &
+ PDZZ2D(J1D,KK)*0.5
+ ! Test if it rests some energy to consume
+ ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
+ ! Length travelled by parcel if it rests energy to consume
+ ZLWORK1(J1D)=PDZZ2D(J1D,KK)*0.5
+ ! Lenght travelled by parcel to nullify energy
+ ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
+ ( ZHLVPT(J1D,KK) - ZVPT_DEP(J1D) ) &
+ - XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))) &
+ + SQRT (ABS( &
+ (XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))) + &
+ PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
+ * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )) ) / &
+ ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )
+ ! Effective length travelled by parcel
+ PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
+ (1-ZTEST)*ZLWORK2(J1D))
+ ! Rest of energy to consume
+ ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
+ ENDDO
+ ELSE
+ ZVPT_DEP(:)=PVPT(:,KK) ! departure point is on mass level
+ ENDIF
+
+ DO JKK=KK+KKL,KKE,KKL
+ IF(ZTESTM > 0.) THEN
+ ZTESTM=0
+ DO J1D=1,IIJU
+ ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
+ ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
+ (ZHLVPT(J1D,JKK) - ZVPT_DEP(J1D)) &
+ + XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))))* PDZZ2D(J1D,JKK)
+ ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
+ ZTESTM=ZTESTM+ZTEST0
+ ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
+ !ZLWORK2 jump of the last reached level
+ ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
+ ( PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D) ) &
+ - XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ + SQRT (ABS( &
+ (XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) + &
+ PG_O_THVREF(J1D) * (PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
+ * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
+ ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
+ !
+ PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
+ (1-ZTEST)*ZLWORK2(J1D))
+ ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
+ END DO
+ ENDIF
+ END DO
+ENDIF
+!!
+!* 2. CALCULATION OF THE DOWNWARD MIXING LENGTH
+! ---------------------------------------
+!
+
+IF (OUPORDN.EQV..FALSE.) THEN
+ IF(OFLUX) CALL PRINT_MSG(NVERB_FATAL,'GEN','COMPUTE_BL89_ML','OFLUX option not coded for downward mixing length')
+ ZINTE(:)=PTKEM_DEP(:)
+ PLWORK=0.
+ ZTESTM=1.
+ DO JKK=KK,KKB,-KKL
+ IF(ZTESTM > 0.) THEN
+ ZTESTM=0
+ DO J1D=1,IIJU
+ ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
+ ZPOTE(J1D) = ZTEST0*(-PG_O_THVREF(J1D) * &
+ (ZHLVPT(J1D,JKK) - PVPT(J1D,KK)) &
+ + XRM17*PSHEAR(J1D,JKK)*SQRT(ABS(PTKEM_DEP(J1D))))* PDZZ2D(J1D,JKK)
+ ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
+ ZTESTM=ZTESTM+ZTEST0
+ ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
+ ZLWORK2(J1D)= ( + PG_O_THVREF(J1D) * &
+ ( PVPT(J1D,JKK) - PVPT(J1D,KK) ) &
+ -XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) &
+ + SQRT (ABS( &
+ (XRM17*PSHEAR(J1D,JKK)*sqrt(abs(PTKEM_DEP(J1D))) - &
+ PG_O_THVREF(J1D) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
+ * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
+ ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
+ !
+ PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
+ (1-ZTEST)*ZLWORK2(J1D))
+ ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
+ END DO
+ ENDIF
+ END DO
+ENDIF
+
+IF (LHOOK) CALL DR_HOOK('COMPUTE_BL89_ML',1,ZHOOK_HANDLE)
+END SUBROUTINE COMPUTE_BL89_ML