diff --git a/src/arome/turb/bl89.F90 b/src/arome/turb/bl89.F90
deleted file mode 100644
index b4191c6151748166ce025e6077ddc437a0ff8423..0000000000000000000000000000000000000000
--- a/src/arome/turb/bl89.F90
+++ /dev/null
@@ -1,310 +0,0 @@
-! ######spl
- SUBROUTINE BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,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
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_CONF, ONLY: CPROGRAM
-USE MODD_CST
-USE MODD_CTURB
-USE MODD_PARAMETERS
-!
-!
-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(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
-! ! 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
-!-------------------------------------------------------------------------------
-!
-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 /) )
- ZTKEM (:,JK) = RESHAPE(PTKEM (:,:,JK),(/ IIU*IJU /) )
- ZG_O_THVREF(:,JK) = RESHAPE(XG/PTHVREF(:,:,JK),(/ IIU*IJU /) )
- DO JRR=1,KRR
- ZRM (:,JK,JRR) = RESHAPE(PRM (:,:,JK,JRR),(/ IIU*IJU /) )
- END DO
- END DO
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!* 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))
- ZPOTE = -ZTEST0*ZG_O_THVREF(J1D,JK)* &
- (ZHLVPT(J1D,JKK)-ZVPT(J1D,JK) )*ZDZZ(J1D,JKK)
- ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE)
- ZTESTM=ZTESTM+ZTEST0
- ZLWORK1=ZDZZ(J1D,JKK)
- 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))
- 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))
- ZPOTE = ZTEST0*ZG_O_THVREF(J1D,JK) * &
- (ZHLVPT(J1D,JKK) - ZVPT(J1D,JK) ) *ZDZZ(J1D,JKK)
- ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE)
- ZTESTM=ZTESTM+ZTEST0
- ZLWORK1=ZDZZ(J1D,JKK)
- 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) )
- 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)
- ZLWORK2=MAX(ZLWORK(J1D),1.E-10)
- ZPOTE = ZLWORK1 / ZLWORK2
- ZLWORK2=1.d0 + ZPOTE**(2./3.)
- ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2))
- 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/arome/turb/compute_bl89_ml.F90 b/src/arome/turb/compute_bl89_ml.F90
deleted file mode 100644
index 696447465f2c90a2e43847a416b519c5725999ff..0000000000000000000000000000000000000000
--- a/src/arome/turb/compute_bl89_ml.F90
+++ /dev/null
@@ -1,205 +0,0 @@
-! ######spl
- SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
- PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,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
-!!
-!-------------------------------------------------------------------------------
-!
-!* 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 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
-
-! 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))) * &
- 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) ) &
- + SQRT (ABS( &
- ( 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))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
- 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) ) &
- + SQRT (ABS( &
- ( 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) THEN
- WRITE(*,*) ' STOP'
- WRITE(*,*) ' OFLUX OPTION NOT CODED FOR DOWNWARD MIXING LENGTH'
- CALL ABORT
- STOP
- ENDIF
- 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))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
- 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) ) &
- + SQRT (ABS( &
- ( 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
diff --git a/src/arome/turb/sbl_depth.F90 b/src/arome/turb/sbl_depth.F90
deleted file mode 100644
index f1c4c22c2e9f72130858f3622d87030f9672cb05..0000000000000000000000000000000000000000
--- a/src/arome/turb/sbl_depth.F90
+++ /dev/null
@@ -1,121 +0,0 @@
-! ######spl
- SUBROUTINE SBL_DEPTH(KKB,KKE,PZZ,PFLXU,PFLXV,PWTHV,PLMO,PSBL_DEPTH)
- USE PARKIND1, ONLY : JPRB
- USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-! #################################################################
-!
-!
-!!**** *SBL_DEPTH* - computes SBL depth
-!!
-!! PURPOSE
-!! -------
-!
-!!** METHOD
-!! ------
-!!
-!! SBL is defined as the layer where momentum flux is equal to XSBL_FRAC of its surface value
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!! V. Masson * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original nov. 2005
-!!
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_PARAMETERS, ONLY : XUNDEF
-USE MODD_CTURB, ONLY : XFTOP_O_FSURF, XSBL_O_BL
-!
-USE MODI_BL_DEPTH_DIAG
-!
-IMPLICIT NONE
-!
-!
-!* 0.1 declarations of arguments
-!
-INTEGER, INTENT(IN) :: KKB ! first physical level
-INTEGER, INTENT(IN) :: KKE ! upper physical level
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitude of flux levels
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXU ! u'w'
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXV ! v'w'
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWTHV ! buoyancy flux
-REAL, DIMENSION(:,:), INTENT(IN) :: PLMO ! Monin-Obukhov length
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PSBL_DEPTH! boundary layer height
-!
-!-------------------------------------------------------------------------------
-!
-! 0.2 declaration of local variables
-!
-!
-INTEGER :: JLOOP ! loop counter
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZQ0 ! surface buoyancy flux
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZWU ! surface friction u'w'
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZWV ! surface friction v'w'
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZUSTAR2 ! surface friction
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZSBL_DYN ! SBL wih dynamical criteria
-REAL, DIMENSION(SIZE(PFLXU,1),SIZE(PFLXU,2),SIZE(PFLXU,3)) :: ZWIND
- ! intermediate wind for SBL calculation
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZSBL_THER! SBL wih thermal criteria
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZA ! ponderation coefficient
-!----------------------------------------------------------------------------
-!
-!* initialisations
-!
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK('SBL_DEPTH',0,ZHOOK_HANDLE)
-ZWU (:,:) = PFLXU(:,:,KKB)
-ZWV (:,:) = PFLXV(:,:,KKB)
-ZQ0 (:,:) = PWTHV(:,:,KKB)
-!
-ZUSTAR2(:,:) = SQRT(ZWU**2+ZWV**2)
-!
-!----------------------------------------------------------------------------
-!
-!* BL and SBL diagnosed with friction criteria
-!
-ZWIND=SQRT(PFLXU**2+PFLXV**2)
-ZSBL_DYN = XSBL_O_BL * BL_DEPTH_DIAG(KKB,KKE,ZUSTAR2,PZZ(:,:,KKB),ZWIND,PZZ,XFTOP_O_FSURF)
-!
-!----------------------------------------------------------------------------
-!
-!* BL and SBL diagnosed with buoyancy flux criteria
-!
-ZSBL_THER= XSBL_O_BL * BL_DEPTH_DIAG(KKB,KKE,ZQ0,PZZ(:,:,KKB),PWTHV,PZZ,XFTOP_O_FSURF)
-!
-!----------------------------------------------------------------------------
-!
-!* SBL depth
-!
-PSBL_DEPTH = 0.
-WHERE (ZSBL_THER> 0. .AND. ZSBL_DYN> 0.) PSBL_DEPTH = MIN(ZSBL_THER(:,:),ZSBL_DYN(:,:))
-WHERE (ZSBL_THER> 0. .AND. ZSBL_DYN==0.) PSBL_DEPTH = ZSBL_THER(:,:)
-WHERE (ZSBL_THER==0. .AND. ZSBL_DYN> 0.) PSBL_DEPTH = ZSBL_THER(:,:)
-!
-DO JLOOP=1,5
- WHERE (PLMO(:,:)/=XUNDEF .AND. ABS(PLMO(:,:))>=0.01 )
- ZA = TANH(2.*PSBL_DEPTH/PLMO)**2
- PSBL_DEPTH = 0.2 * PSBL_DEPTH + 0.8 * ((1.-ZA) * ZSBL_DYN + ZA * ZSBL_THER )
- END WHERE
-END DO
-WHERE (ABS(PLMO(:,:))<=0.01 ) PSBL_DEPTH = ZSBL_DYN
-WHERE (PLMO(:,:)==XUNDEF) PSBL_DEPTH = ZSBL_THER
-!
-!----------------------------------------------------------------------------
-IF (LHOOK) CALL DR_HOOK('SBL_DEPTH',1,ZHOOK_HANDLE)
-END SUBROUTINE SBL_DEPTH
diff --git a/src/mesonh/turb/bl89.f90 b/src/mesonh/turb/bl89.f90
deleted file mode 100644
index 8d9fe3e369828bca9115ce04c6ea4b01cfb18570..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/bl89.f90
+++ /dev/null
@@ -1,396 +0,0 @@
-!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
-!
-!
-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,J3RD ! 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
-!-------------------------------------------------------------------------------
-!
-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)
-
- !-------- 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) &
- + 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
- 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 -----------
- 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
-
-!
-END SUBROUTINE BL89
diff --git a/src/mesonh/turb/compute_bl89_ml.f90 b/src/mesonh/turb/compute_bl89_ml.f90
deleted file mode 100644
index 20c9a078db5e550dbf3bb03cd2fcc6b13ddc89a5..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/compute_bl89_ml.f90
+++ /dev/null
@@ -1,250 +0,0 @@
-!MNH_LIC Copyright 2006-2019 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_COMPUTE_BL89_ML
-! ###########################
-
-INTERFACE
-
-! ###################################################################
- SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
- PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PSHEAR,PLWORK)
-! ###################################################################
-
-!* 1.1 Declaration of 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
-
-END SUBROUTINE COMPUTE_BL89_ML
-
-END INTERFACE
-!
-END MODULE MODI_COMPUTE_BL89_ML
-! ######spl
- SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
- PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PSHEAR,PLWORK)
-
-! ###################################################################
-!!
-!! 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
-! --------------
-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
-
-END SUBROUTINE COMPUTE_BL89_ML
diff --git a/src/mesonh/turb/sbl_depth.f90 b/src/mesonh/turb/sbl_depth.f90
deleted file mode 100644
index e83d8f784d4b1f4b84e755fe227fff2c588f24fc..0000000000000000000000000000000000000000
--- a/src/mesonh/turb/sbl_depth.f90
+++ /dev/null
@@ -1,145 +0,0 @@
-!MNH_LIC Copyright 1994-2014 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_SBL_DEPTH
-! ################
-!
-INTERFACE
-!
- SUBROUTINE SBL_DEPTH(KKB,KKE,PZZ,PFLXU,PFLXV,PWTHV,PLMO,PSBL_DEPTH)
-!
-INTEGER, INTENT(IN) :: KKB ! first physical level
-INTEGER, INTENT(IN) :: KKE ! upper physical level
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitude of flux levels
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXU ! u'w'
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXV ! v'w'
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWTHV ! buoyancy flux
-REAL, DIMENSION(:,:), INTENT(IN) :: PLMO ! Monin-Obukhov length
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PSBL_DEPTH! boundary layer height
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE SBL_DEPTH
-!
-END INTERFACE
-!
-END MODULE MODI_SBL_DEPTH
-!
-! #################################################################
- SUBROUTINE SBL_DEPTH(KKB,KKE,PZZ,PFLXU,PFLXV,PWTHV,PLMO,PSBL_DEPTH)
-! #################################################################
-!
-!
-!!**** *SBL_DEPTH* - computes SBL depth
-!!
-!! PURPOSE
-!! -------
-!
-!!** METHOD
-!! ------
-!!
-!! SBL is defined as the layer where momentum flux is equal to XSBL_FRAC of its surface value
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!! V. Masson * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original nov. 2005
-!! 26/02/2020 T.Nagel Correction of SBL depth computation in neutral stratification
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_PARAMETERS, ONLY : XUNDEF
-USE MODD_CTURB, ONLY : XFTOP_O_FSURF, XSBL_O_BL
-!
-USE MODI_BL_DEPTH_DIAG
-!
-IMPLICIT NONE
-!
-!
-!* 0.1 declarations of arguments
-!
-INTEGER, INTENT(IN) :: KKB ! first physical level
-INTEGER, INTENT(IN) :: KKE ! upper physical level
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! altitude of flux levels
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXU ! u'w'
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFLXV ! v'w'
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWTHV ! buoyancy flux
-REAL, DIMENSION(:,:), INTENT(IN) :: PLMO ! Monin-Obukhov length
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PSBL_DEPTH! boundary layer height
-!
-!-------------------------------------------------------------------------------
-!
-! 0.2 declaration of local variables
-!
-!
-INTEGER :: JLOOP ! loop counter
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZQ0 ! surface buoyancy flux
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZWU ! surface friction u'w'
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZWV ! surface friction v'w'
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZUSTAR2 ! surface friction
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZSBL_DYN ! SBL wih dynamical criteria
-REAL, DIMENSION(SIZE(PFLXU,1),SIZE(PFLXU,2),SIZE(PFLXU,3)) :: ZWIND
- ! intermediate wind for SBL calculation
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZSBL_THER! SBL wih thermal criteria
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2)) :: ZA ! ponderation coefficient
-!----------------------------------------------------------------------------
-!
-!* initialisations
-!
-!
-ZWU (:,:) = PFLXU(:,:,KKB)
-ZWV (:,:) = PFLXV(:,:,KKB)
-ZQ0 (:,:) = PWTHV(:,:,KKB)
-!
-ZUSTAR2(:,:) = SQRT(ZWU**2+ZWV**2)
-!
-!----------------------------------------------------------------------------
-!
-!* BL and SBL diagnosed with friction criteria
-!
-ZWIND=SQRT(PFLXU**2+PFLXV**2)
-ZSBL_DYN = XSBL_O_BL * BL_DEPTH_DIAG(KKB,KKE,ZUSTAR2,PZZ(:,:,KKB),ZWIND,PZZ,XFTOP_O_FSURF)
-!
-!----------------------------------------------------------------------------
-!
-!* BL and SBL diagnosed with buoyancy flux criteria
-!
-ZSBL_THER= XSBL_O_BL * BL_DEPTH_DIAG(KKB,KKE,ZQ0,PZZ(:,:,KKB),PWTHV,PZZ,XFTOP_O_FSURF)
-!
-!----------------------------------------------------------------------------
-!
-!* SBL depth
-!
-PSBL_DEPTH = 0.
-WHERE (ZSBL_THER> 0. .AND. ZSBL_DYN> 0.) PSBL_DEPTH = MIN(ZSBL_THER(:,:),ZSBL_DYN(:,:))
-WHERE (ZSBL_THER> 0. .AND. ZSBL_DYN==0.) PSBL_DEPTH = ZSBL_THER(:,:)
-WHERE (ZSBL_THER==0. .AND. ZSBL_DYN> 0.) PSBL_DEPTH = ZSBL_DYN(:,:)
-!
-DO JLOOP=1,5
- WHERE (PLMO(:,:)/=XUNDEF .AND. ABS(PLMO(:,:))>=0.01 )
- ZA = TANH(2.*PSBL_DEPTH/PLMO)**2
- PSBL_DEPTH = 0.2 * PSBL_DEPTH + 0.8 * ((1.-ZA) * ZSBL_DYN + ZA * ZSBL_THER )
- END WHERE
-END DO
-WHERE (ABS(PLMO(:,:))<=0.01 ) PSBL_DEPTH = ZSBL_THER
-WHERE (PLMO(:,:)==XUNDEF) PSBL_DEPTH = ZSBL_DYN
-!
-!----------------------------------------------------------------------------
-END SUBROUTINE SBL_DEPTH