turb.F90

          PDP,PTP,PSIGS,PWTH,PWRC,PWSV,                  &
          PSSTFL, PSSTFL_C, PSSRFL_C,PSSUFL_C,PSSVFL_C,  &
          PSSUFL,PSSVFL                                  )

IF( BUCONF%LBUDGET_U ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_U), 'VTURB', PRUS(:,:) )
IF( BUCONF%LBUDGET_V ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_V), 'VTURB', PRVS(:,:) )
IF( BUCONF%LBUDGET_W ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_W), 'VTURB', PRWS(:,:) )

IF( BUCONF%LBUDGET_TH ) THEN
  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) &
                                                                          + ZLSOCPEXNM(:,:) * PRRS(:,:, 4) )
  ELSE IF( KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2) )
  ELSE
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:,:) )
  END IF
END IF

IF( BUCONF%LBUDGET_RV ) THEN
  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4) )
  ELSE IF( KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:,:, 1) - PRRS(:,:, 2) )
  ELSE
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:,:, 1) )
  END IF
END IF

IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'VTURB', PRRS(:,:, 2) )
IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'VTURB', PRRS(:,:, 4) )

IF( BUCONF%LBUDGET_SV )  THEN
  DO JSV = 1, KSV
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'VTURB', PRSVS(:,:, JSV) )
  END DO
END IF
!
!Les budgets des termes horizontaux de la turb sont présents dans AROME
! alors que ces termes ne sont pas calculés
#ifdef REPRO48
#else
IF( HTURBDIM == '3DIM' ) THEN
#endif
  IF( BUCONF%LBUDGET_U  ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_U ), 'HTURB', PRUS  (:,:) )
  IF( BUCONF%LBUDGET_V  ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_V ), 'HTURB', PRVS  (:,:) )
  IF( BUCONF%LBUDGET_W  ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_W ), 'HTURB', PRWS  (:,:) )

  IF(BUCONF%LBUDGET_TH)  THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) &
                                                                             + ZLSOCPEXNM(:,:) * PRRS(:,:, 4) )
    ELSE IF( KRRL >= 1 ) THEN
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2) )
    ELSE
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RV ) THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4) )
    ELSE IF( KRRL >= 1 ) THEN
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) - PRRS(:,:, 2) )
    ELSE
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:,:, 2) )
  IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:,:, 4) )

  IF( BUCONF%LBUDGET_SV )  THEN
    DO JSV = 1, KSV
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:,:, JSV) )
    END DO
  END IF
!à supprimer une fois le précédent ifdef REPRO48 validé
#ifdef REPRO48
#else
    CALL TURB_HOR_SPLT(D,CST,CSTURB,                           &
          KSPLIT, KRR, KRRL, KRRI, KSV, PTSTEP,HLBCX,HLBCY,    &
          OTURB_FLX,OSUBG_COND,OOCEAN,OCOMPUTE_SRC,            &
          TPFILE,                                              &
          PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,                        &
          PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,                       &
          PCOSSLOPE,PSINSLOPE,                                 &
          PRHODJ,PTHVREF,                                      &
          PSFTH,PSFRV,PSFSV,                                   &
          ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU22M,ZTAU33M,             &
          PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,          &
          PTKET,ZLM,ZLEPS,                                     &
          ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE,           &
          PDP,PTP,PSIGS,                                       &
          ZTRH,                                                &
          PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS                     )
#endif
  IF( BUCONF%LBUDGET_U ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_U), 'HTURB', PRUS(:,:) )
  IF( BUCONF%LBUDGET_V ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_V), 'HTURB', PRVS(:,:) )
  IF( BUCONF%LBUDGET_W ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_W), 'HTURB', PRWS(:,:) )

  IF( BUCONF%LBUDGET_TH ) THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) &
                                                                            + ZLSOCPEXNM(:,:) * PRRS(:,:, 4) )
    ELSE IF( KRRL >= 1 ) THEN
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2) )
    ELSE
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RV ) THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4) )
    ELSE IF( KRRL >= 1 ) THEN
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) - PRRS(:,:, 2) )
    ELSE
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:,:, 2) )
  IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:,:, 4) )

  IF( BUCONF%LBUDGET_SV )  THEN
    DO JSV = 1, KSV
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:,:, JSV) )
    END DO
  END IF
#ifdef REPRO48
#else
END IF
#endif
!----------------------------------------------------------------------------
!
!*      6. EVOLUTION OF THE TKE AND ITS DISSIPATION
!          ----------------------------------------
!
!  6.1 Contribution of mass-flux in the TKE buoyancy production if
!      cloud computation is not statistical
CALL MZF_PHY(D,PFLXZTHVMF,ZWORK1)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
PTP(IIJB:IIJE,1:D%NKT) = PTP(IIJB:IIJE,1:D%NKT) &
                             + CST%XG / PTHVREF(IIJB:IIJE,1:D%NKT) * ZWORK1(IIJB:IIJE,1:D%NKT)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)

IF(PRESENT(PTPMF))  THEN
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  PTPMF(IIJB:IIJE,1:D%NKT)=CST%XG / PTHVREF(IIJB:IIJE,1:D%NKT) * ZWORK1(IIJB:IIJE,1:D%NKT)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
END IF
!  6.2 TKE evolution equation

IF (.NOT. OHARAT) THEN
!
IF (BUCONF%LBUDGET_TH)  THEN
  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:)+ ZLVOCPEXNM(:,:) * PRRS(:,:,2) &
                                                          & + ZLSOCPEXNM(:,:) * PRRS(:,:,4) )
  ELSE IF ( KRRL >= 1 ) THEN
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:,2) )
  ELSE
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:) )
  END IF
END IF
!
IF(PRESENT(PRTKEMS)) THEN
  ZRTKEMS(:,:)=PRTKEMS(:,:)
ELSE
  ZRTKEMS(:,:)=0.
END IF
!
CALL TKE_EPS_SOURCES(D,CST,CSTURB,BUCONF,HPROGRAM,                      &
                   & KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH,                      &
                   & PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,               &
                   & PTSTEP,PIMPL,ZEXPL,                                &
                   & HTURBLEN,HTURBDIM,                                 &
                   & TPFILE,OTURB_DIAG,OLES_CALL,ODIAG_IN_RUN,          &
                   & PTP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,ZRTKEMS,&
                   & TBUDGETS,KBUDGETS, PEDR=PEDR, PTR=PTR,PDISS=PDISS, &
                   & PCURRENT_TKE_DISS=PCURRENT_TKE_DISS                )
IF (BUCONF%LBUDGET_TH)  THEN
  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:)+ ZLVOCPEXNM(:,:) * PRRS(:,:,2) &
                                                          & + ZLSOCPEXNM(:,:) * PRRS(:,:,4) )
  ELSE IF ( KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:)+ ZLOCPEXNM(:,:) * PRRS(:,:,2) )
  ELSE
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:) )
  END IF
END IF
!
ENDIF
!
!----------------------------------------------------------------------------
!
!*      7. STORES SOME INFORMATIONS RELATED TO THE TURBULENCE SCHEME
!          ---------------------------------------------------------
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
!
! stores the mixing length
!
  TZFIELD%CMNHNAME   = 'LM'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'LM'
  TZFIELD%CUNITS     = 'm'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'Mixing length'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZLM)
!
  IF (KRR /= 0) THEN
!
! stores the conservative potential temperature
!
    TZFIELD%CMNHNAME   = 'THLM'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'THLM'
    TZFIELD%CUNITS     = 'K'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = 'Conservative potential temperature'
    TZFIELD%NGRID      = 1
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 3
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,PTHLT)
!
! stores the conservative mixing ratio
!
    TZFIELD%CMNHNAME   = 'RNPM'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'RNPM'
    TZFIELD%CUNITS     = 'kg kg-1'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = 'Conservative mixing ratio'
    TZFIELD%NGRID      = 1
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 3
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,PRT(:,:,1))
   END IF
END IF
!
!* stores value of conservative variables & wind before turbulence tendency (AROME only)
IF(PRESENT(PDRUS_TURB)) THEN
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  PDRUS_TURB(IIJB:IIJE,1:D%NKT)   = PRUS(IIJB:IIJE,1:D%NKT) - PDRUS_TURB(IIJB:IIJE,1:D%NKT)
  PDRVS_TURB(IIJB:IIJE,1:D%NKT)   = PRVS(IIJB:IIJE,1:D%NKT) - PDRVS_TURB(IIJB:IIJE,1:D%NKT)
  PDRTHLS_TURB(IIJB:IIJE,1:D%NKT) = PRTHLS(IIJB:IIJE,1:D%NKT) - PDRTHLS_TURB(IIJB:IIJE,1:D%NKT)
  PDRRTS_TURB(IIJB:IIJE,1:D%NKT)  = PRRS(IIJB:IIJE,1:D%NKT,1) - PDRRTS_TURB(IIJB:IIJE,1:D%NKT)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT,JSV=1:KSV)  
  PDRSVS_TURB(IIJB:IIJE,1:D%NKT,:)  = PRSVS(IIJB:IIJE,1:D%NKT,:) - PDRSVS_TURB(IIJB:IIJE,1:D%NKT,:)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT,JSV=1:KSV)
END IF
!----------------------------------------------------------------------------
!
!*      8. RETRIEVE NON-CONSERVATIVE VARIABLES
!          -----------------------------------
!
IF ( KRRL >= 1 ) THEN
  IF ( KRRI >= 1 ) THEN
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    PRT(IIJB:IIJE,1:D%NKT,1)  = PRT(IIJB:IIJE,1:D%NKT,1)  - PRT(IIJB:IIJE,1:D%NKT,2)  &
                                    - PRT(IIJB:IIJE,1:D%NKT,4)
    PRRS(IIJB:IIJE,1:D%NKT,1) = PRRS(IIJB:IIJE,1:D%NKT,1) - PRRS(IIJB:IIJE,1:D%NKT,2) &
                                    - PRRS(IIJB:IIJE,1:D%NKT,4)
    PTHLT(IIJB:IIJE,1:D%NKT)  = PTHLT(IIJB:IIJE,1:D%NKT)  + ZLVOCPEXNM(IIJB:IIJE,1:D%NKT) &
                                    * PRT(IIJB:IIJE,1:D%NKT,2) &
                                    + ZLSOCPEXNM(IIJB:IIJE,1:D%NKT) * PRT(IIJB:IIJE,1:D%NKT,4)
    PRTHLS(IIJB:IIJE,1:D%NKT) = PRTHLS(IIJB:IIJE,1:D%NKT) + ZLVOCPEXNM(IIJB:IIJE,1:D%NKT) &
                                    * PRRS(IIJB:IIJE,1:D%NKT,2) &
                                    + ZLSOCPEXNM(IIJB:IIJE,1:D%NKT) * PRRS(IIJB:IIJE,1:D%NKT,4)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
  ELSE
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    PRT(IIJB:IIJE,1:D%NKT,1)  = PRT(IIJB:IIJE,1:D%NKT,1)  - PRT(IIJB:IIJE,1:D%NKT,2)
    PRRS(IIJB:IIJE,1:D%NKT,1) = PRRS(IIJB:IIJE,1:D%NKT,1) - PRRS(IIJB:IIJE,1:D%NKT,2)
    PTHLT(IIJB:IIJE,1:D%NKT)  = PTHLT(IIJB:IIJE,1:D%NKT)  + ZLOCPEXNM(IIJB:IIJE,1:D%NKT) &
                                    * PRT(IIJB:IIJE,1:D%NKT,2)
    PRTHLS(IIJB:IIJE,1:D%NKT) = PRTHLS(IIJB:IIJE,1:D%NKT) + ZLOCPEXNM(IIJB:IIJE,1:D%NKT) &
                                    * PRRS(IIJB:IIJE,1:D%NKT,2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  END IF
END IF

! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
CALL SOURCES_NEG_CORRECT_PHY(D,KSV,HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
!----------------------------------------------------------------------------
!
!*      9. LES averaged surface fluxes
!          ---------------------------
!
IF (OLES_CALL) THEN
  CALL SECOND_MNH(ZTIME1)
  CALL LES_MEAN_SUBGRID_PHY(D,PSFTH,X_LES_Q0)
  CALL LES_MEAN_SUBGRID_PHY(D,PSFRV,X_LES_E0)
  DO JSV=1,KSV
    CALL LES_MEAN_SUBGRID_PHY(D,PSFSV(:,JSV),X_LES_SV0(:,JSV))
  END DO
  CALL LES_MEAN_SUBGRID_PHY(D,PSFU,X_LES_UW0)
  CALL LES_MEAN_SUBGRID_PHY(D,PSFV,X_LES_VW0)
  !
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  ZWORK2D(IIJB:IIJE) = (PSFU(IIJB:IIJE)*PSFU(IIJB:IIJE)+PSFV(IIJB:IIJE)*PSFV(IIJB:IIJE))**0.25
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
  CALL LES_MEAN_SUBGRID_PHY(D,ZWORK2D,X_LES_USTAR)
!----------------------------------------------------------------------------
!
!*     10. LES for 3rd order moments
!          -------------------------
!
  CALL LES_MEAN_SUBGRID_PHY(D,ZMWTH,X_LES_SUBGRID_W2Thl)
  CALL LES_MEAN_SUBGRID_PHY(D,ZMTH2,X_LES_SUBGRID_WThl2)
  IF (KRR>0) THEN
    CALL LES_MEAN_SUBGRID_PHY(D,ZMWR,X_LES_SUBGRID_W2Rt)
    CALL LES_MEAN_SUBGRID_PHY(D,ZMTHR,X_LES_SUBGRID_WThlRt)
    CALL LES_MEAN_SUBGRID_PHY(D,ZMR2,X_LES_SUBGRID_WRt2)
  END IF
!
!----------------------------------------------------------------------------
!
!*     11. LES quantities depending on <w'2> in "1DIM" mode
!          ------------------------------------------------
!
  IF (HTURBDIM=="1DIM") THEN
    !
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    ZWORK1(IIJB:IIJE,1:D%NKT) = 2./3.*PTKET(IIJB:IIJE,1:D%NKT)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    CALL LES_MEAN_SUBGRID_PHY(D,ZWORK1,X_LES_SUBGRID_U2)
    X_LES_SUBGRID_V2(:,:,:) = X_LES_SUBGRID_U2(:,:,:)
    X_LES_SUBGRID_W2(:,:,:) = X_LES_SUBGRID_U2(:,:,:)
    !
    CALL GZ_M_W_PHY(D,PTHLT,PDZZ,ZWORK1)
    CALL MZF_PHY(D,ZWORK1,ZWORK2)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    ZWORK2(IIJB:IIJE,1:D%NKT)  = 2./3.*PTKET(IIJB:IIJE,1:D%NKT) *ZWORK2(IIJB:IIJE,1:D%NKT)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    CALL LES_MEAN_SUBGRID_PHY(D,ZWORK2,X_LES_RES_ddz_Thl_SBG_W2)
    !
    IF (KRR>=1) THEN
      CALL GZ_M_W_PHY(D,PRT(:,:,1),PDZZ,ZWORK1)
      CALL MZF_PHY(D,ZWORK1,ZWORK2)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      ZWORK2(IIJB:IIJE,1:D%NKT)  = 2./3.*PTKET(IIJB:IIJE,1:D%NKT) *ZWORK2(IIJB:IIJE,1:D%NKT)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      CALL LES_MEAN_SUBGRID_PHY(D,ZWORK2,X_LES_RES_ddz_Rt_SBG_W2)
    END IF
    DO JSV=1,KSV
      CALL GZ_M_W_PHY(D,PSVT(:,:,JSV),PDZZ,ZWORK1)
      CALL MZF_PHY(D,ZWORK1,ZWORK2)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      ZWORK2(IIJB:IIJE,1:D%NKT)  = 2./3.*PTKET(IIJB:IIJE,1:D%NKT) *ZWORK2(IIJB:IIJE,1:D%NKT)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      CALL LES_MEAN_SUBGRID_PHY(D,ZWORK2, X_LES_RES_ddz_Sv_SBG_W2(:,:,:,JSV))
    END DO
  END IF

!----------------------------------------------------------------------------
!
!*     12. LES mixing end dissipative lengths, presso-correlations
!          -------------------------------------------------------
!
  CALL LES_MEAN_SUBGRID_PHY(D,ZLM,X_LES_SUBGRID_LMix)
  CALL LES_MEAN_SUBGRID_PHY(D,ZLEPS,X_LES_SUBGRID_LDiss)
!
!* presso-correlations for subgrid Tke are equal to zero.
!
  ZLEPS(:,:) = 0. !ZLEPS is used as a work array (not used anymore)
  CALL LES_MEAN_SUBGRID_PHY(D,ZLEPS,X_LES_SUBGRID_WP)
!
  CALL SECOND_MNH(ZTIME2)
  XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
IF(PRESENT(PLEM)) PLEM(IIJB:IIJE,IKTB:IKTE) = ZLM(IIJB:IIJE,IKTB:IKTE)
!----------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
CONTAINS
!
!     ########################################################################
      SUBROUTINE COMPUTE_FUNCTION_THERMO(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
                                         PLOCPEXN,PAMOIST,PATHETA            )
!     ########################################################################
!!
!!****  *COMPUTE_FUNCTION_THERMO* routine to compute several thermo functions
!
!!    AUTHOR
!!    ------
!!
!!     JP Pinty      *LA*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   24/02/03
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments
!
REAL,                   INTENT(IN)    :: PALP,PBETA,PGAM,PLTT,PC
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN)    :: PT,PEXN,PCP
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLOCPEXN
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PAMOIST,PATHETA
!
!-------------------------------------------------------------------------------
!
  IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',0,ZHOOK_HANDLE2)
  ZEPS = CST%XMV / CST%XMD
!
!*       1.1 Lv/Cph at  t
!
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  PLOCPEXN(IIJB:IIJE,1:D%NKT) = ( PLTT + (CST%XCPV-PC) *  (PT(IIJB:IIJE,1:D%NKT)-CST%XTT) ) &
                                     / PCP(IIJB:IIJE,1:D%NKT)
!
!*      1.2 Saturation vapor pressure at t
!
  ZRVSAT(IIJB:IIJE,1:D%NKT) =  EXP( PALP - PBETA/PT(IIJB:IIJE,1:D%NKT) - PGAM*ALOG( PT(IIJB:IIJE,1:D%NKT) ) )
!
!*      1.3 saturation  mixing ratio at t
!
  ZRVSAT(IIJB:IIJE,1:D%NKT) =  ZRVSAT(IIJB:IIJE,1:D%NKT) &
                                    * ZEPS / ( PPABST(IIJB:IIJE,1:D%NKT) - ZRVSAT(IIJB:IIJE,1:D%NKT) )
!
!*      1.4 compute the saturation mixing ratio derivative (rvs')
!
  ZDRVSATDT(IIJB:IIJE,1:D%NKT) = ( PBETA / PT(IIJB:IIJE,1:D%NKT)  - PGAM ) / PT(IIJB:IIJE,1:D%NKT)   &
                 * ZRVSAT(IIJB:IIJE,1:D%NKT) * ( 1. + ZRVSAT(IIJB:IIJE,1:D%NKT) / ZEPS )
!
!*      1.5 compute Amoist
!
  PAMOIST(IIJB:IIJE,1:D%NKT)=  0.5 / ( 1.0 + ZDRVSATDT(IIJB:IIJE,1:D%NKT) * PLOCPEXN(IIJB:IIJE,1:D%NKT) )
!
!*      1.6 compute Atheta
!
  PATHETA(IIJB:IIJE,1:D%NKT)= PAMOIST(IIJB:IIJE,1:D%NKT) * PEXN(IIJB:IIJE,1:D%NKT) *               &
        ( ( ZRVSAT(IIJB:IIJE,1:D%NKT) - PRT(IIJB:IIJE,1:D%NKT,1) ) * PLOCPEXN(IIJB:IIJE,1:D%NKT) / &
          ( 1. + ZDRVSATDT(IIJB:IIJE,1:D%NKT) * PLOCPEXN(IIJB:IIJE,1:D%NKT) )        *               &
          (                                                                  &
           ZRVSAT(IIJB:IIJE,1:D%NKT) * (1. + ZRVSAT(IIJB:IIJE,1:D%NKT)/ZEPS)                         &
                        * ( -2.*PBETA/PT(IIJB:IIJE,1:D%NKT) + PGAM ) / PT(IIJB:IIJE,1:D%NKT)**2      &
          +ZDRVSATDT(IIJB:IIJE,1:D%NKT) * (1. + 2. * ZRVSAT(IIJB:IIJE,1:D%NKT)/ZEPS)                 &
                        * ( PBETA/PT(IIJB:IIJE,1:D%NKT) - PGAM ) / PT(IIJB:IIJE,1:D%NKT)             &
          )                                                                  &
         - ZDRVSATDT(IIJB:IIJE,1:D%NKT)                                                  &
        )
!
!*      1.7 Lv/Cph/Exner at t-1
!
  PLOCPEXN(IIJB:IIJE,1:D%NKT) = PLOCPEXN(IIJB:IIJE,1:D%NKT) / PEXN(IIJB:IIJE,1:D%NKT)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE2)
END SUBROUTINE COMPUTE_FUNCTION_THERMO

!     ########################################################################
      SUBROUTINE COMPUTE_FUNCTION_THERMO_NEW_STAT(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
                                         PLOCPEXN,PAMOIST,PATHETA            )
!     ########################################################################
!!
!!****  *COMPUTE_FUNCTION_THERMO* routine to compute several thermo functions
!
!!    AUTHOR
!!    ------
!!
!!     JP Pinty      *LA*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   24/02/03
!!     Modified: Wim de Rooy 06-02-2019
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
USE MODD_CST
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments
!
REAL, INTENT(IN)                      :: PALP,PBETA,PGAM,PLTT,PC
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN)    :: PT,PEXN,PCP
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLOCPEXN
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PAMOIST,PATHETA
!
!-------------------------------------------------------------------------------
!
  IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO_NEW_STAT',0,ZHOOK_HANDLE2)
  ZEPS = CST%XMV / CST%XMD
!
!*       1.1 Lv/Cph at  t
!
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  PLOCPEXN(IIJB:IIJE,1:D%NKT) = ( PLTT + (CST%XCPV-PC) *  (PT(IIJB:IIJE,1:D%NKT)-CST%XTT) ) / PCP(IIJB:IIJE,1:D%NKT)
!
!*      1.2 Saturation vapor pressure at t
!
  ZRVSAT(IIJB:IIJE,1:D%NKT) =  EXP( PALP - PBETA/PT(IIJB:IIJE,1:D%NKT) - PGAM*ALOG( PT(IIJB:IIJE,1:D%NKT) ) )
!
!*      1.3 saturation  mixing ratio at t
!
  ZRVSAT(IIJB:IIJE,1:D%NKT) =  ZRVSAT(IIJB:IIJE,1:D%NKT) * ZEPS / ( PPABST(IIJB:IIJE,1:D%NKT) - ZRVSAT(IIJB:IIJE,1:D%NKT) )
!
!*      1.4 compute the saturation mixing ratio derivative (rvs')
!
  ZDRVSATDT(IIJB:IIJE,1:D%NKT) = ( PBETA / PT(IIJB:IIJE,1:D%NKT)  - PGAM ) / PT(IIJB:IIJE,1:D%NKT)   &
                 * ZRVSAT(IIJB:IIJE,1:D%NKT) * ( 1. + ZRVSAT(IIJB:IIJE,1:D%NKT) / ZEPS )
!
!*      1.5 compute Amoist
!
  PAMOIST(IIJB:IIJE,1:D%NKT)=  1.0 / ( 1.0 + ZDRVSATDT(IIJB:IIJE,1:D%NKT) * PLOCPEXN(IIJB:IIJE,1:D%NKT) )
!
!*      1.6 compute Atheta
!
  PATHETA(IIJB:IIJE,1:D%NKT)= PAMOIST(IIJB:IIJE,1:D%NKT) * PEXN(IIJB:IIJE,1:D%NKT) * ZDRVSATDT(IIJB:IIJE,1:D%NKT)
!
!*      1.7 Lv/Cph/Exner at t-1
!
  PLOCPEXN(IIJB:IIJE,1:D%NKT) = PLOCPEXN(IIJB:IIJE,1:D%NKT) / PEXN(IIJB:IIJE,1:D%NKT)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO_NEW_STAT',1,ZHOOK_HANDLE2)
END SUBROUTINE COMPUTE_FUNCTION_THERMO_NEW_STAT

!
!     ####################
      SUBROUTINE DELT(PLM,ODZ)
!     ####################
!!
!!****  *DELT* routine to compute mixing length for DELT case
!
!!    AUTHOR
!!    ------
!!
!!     M Tomasini      *Meteo-France
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   01/05
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!*       0.1   Declarations of dummy arguments
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLM
LOGICAL,                INTENT(IN)    :: ODZ
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB:DELT',0,ZHOOK_HANDLE2)
!
CALL MXF_PHY(D,PDXX,ZWORK1)
IF (.NOT. O2D) THEN
  CALL MYF_PHY(D,PDYY,ZWORK2)
END IF
!
IF (ODZ) THEN
  ! Dz is take into account in the computation
  DO JK = IKTB,IKTE ! 1D turbulence scheme
    !$mnh_expand_array(JIJ=IIJB:IIJE)
    PLM(IIJB:IIJE,JK) = PZZ(IIJB:IIJE,JK+D%NKL) - PZZ(IIJB:IIJE,JK)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE)
  END DO
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  PLM(IIJB:IIJE,D%NKU) = PLM(IIJB:IIJE,IKE)
  PLM(IIJB:IIJE,D%NKA) = PZZ(IIJB:IIJE,IKB) - PZZ(IIJB:IIJE,D%NKA)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
  IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
    IF ( O2D) THEN
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      PLM(IIJB:IIJE,1:D%NKT) = SQRT( PLM(IIJB:IIJE,1:D%NKT)*ZWORK1(IIJB:IIJE,1:D%NKT) )
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    ELSE
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      PLM(IIJB:IIJE,1:D%NKT) = (PLM(IIJB:IIJE,1:D%NKT)*ZWORK1(IIJB:IIJE,1:D%NKT) &
                                   * ZWORK2(IIJB:IIJE,1:D%NKT) ) ** (1./3.)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    END IF
  END IF
ELSE
  ! Dz not taken into account in computation to assure invariability with vertical grid mesh
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  PLM(IIJB:IIJE,1:D%NKT)=1.E10
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
    IF ( O2D) THEN
      PLM(:,:) = ZWORK1(:,:)
    ELSE
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
      PLM(IIJB:IIJE,1:D%NKT) = (ZWORK1(IIJB:IIJE,1:D%NKT)*ZWORK2(IIJB:IIJE,1:D%NKT) ) ** (1./2.)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    END IF
  END IF
END IF
!
!  mixing length limited by the distance normal to the surface
!  (with the same factor as for BL89)
!
IF (.NOT. ORMC01) THEN
  ZALPHA=0.5**(-1.5)
  !
  DO JIJ=IIJB,IIJE
    IF (OOCEAN) THEN
      DO JK=IKTE,IKTB,-1
        ZD=ZALPHA*(PZZ(JIJ,IKTE+1)-PZZ(JIJ,JK))
        IF ( PLM(JIJ,JK)>ZD) THEN
          PLM(JIJ,JK)=ZD
        ELSE
          EXIT
        ENDIF
      END DO
    ELSE
      DO JK=IKTB,IKTE
        ZD=ZALPHA*(0.5*(PZZ(JIJ,JK)+PZZ(JIJ,JK+D%NKL))&
        -PZZ(JIJ,IKB)) *PDIRCOSZW(JIJ)
        IF ( PLM(JIJ,JK)>ZD) THEN
          PLM(JIJ,JK)=ZD
        ELSE
          EXIT
        ENDIF
      END DO
    ENDIF
  END DO
END IF
!
!$mnh_expand_array(JIJ=IIJB:IIJE)
PLM(IIJB:IIJE,D%NKA) = PLM(IIJB:IIJE,IKB)
PLM(IIJB:IIJE,D%NKU) = PLM(IIJB:IIJE,IKE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
IF (LHOOK) CALL DR_HOOK('TURB:DELT',1,ZHOOK_HANDLE2)
END SUBROUTINE DELT
!
!     ####################
      SUBROUTINE DEAR(PLM)
!     ####################
!!
!!****  *DEAR* routine to compute mixing length for DEARdorff case
!
!!    AUTHOR
!!    ------
!!
!!     M Tomasini      *Meteo-France
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   01/05
!!      I.Sandu (Sept.2006) : Modification of the stability criterion
!!                            (theta_v -> theta_l)
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!*       0.1   Declarations of dummy arguments
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLM
!
!-------------------------------------------------------------------------------
!
!   initialize the mixing length with the mesh grid
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',0,ZHOOK_HANDLE2)
IF ( HTURBDIM /= '1DIM' ) THEN
  CALL MXF_PHY(D,PDXX,ZWORK1)
  IF (.NOT. O2D) THEN
    CALL MYF_PHY(D,PDYY,ZWORK2)
  END IF
END IF
! 1D turbulence scheme
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=IKTB:IKTE)
PLM(IIJB:IIJE,IKTB:IKTE) = PZZ(IIJB:IIJE,D%NKL+IKTB:IKTE+D%NKL) - PZZ(IIJB:IIJE,IKTB:IKTE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=IKTB:IKTE)
!$mnh_expand_array(JIJ=IIJB:IIJE)
PLM(IIJB:IIJE,D%NKU) = PLM(IIJB:IIJE,IKE)
PLM(IIJB:IIJE,D%NKA) = PZZ(IIJB:IIJE,IKB) - PZZ(IIJB:IIJE,D%NKA)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
  IF ( O2D) THEN
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    PLM(IIJB:IIJE,1:D%NKT) = SQRT( PLM(IIJB:IIJE,1:D%NKT)*ZWORK1(IIJB:IIJE,1:D%NKT) )
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  ELSE
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
    PLM(IIJB:IIJE,1:D%NKT) = (PLM(IIJB:IIJE,1:D%NKT)*ZWORK1(IIJB:IIJE,1:D%NKT) &
                                 * ZWORK2(IIJB:IIJE,1:D%NKT) ) ** (1./3.)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
  END IF
END IF
!   compute a mixing length limited by the stability
!
CALL ETHETA(D,CST,KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT,OOCEAN,OCOMPUTE_SRC,ZETHETA)
CALL EMOIST(D,CST,KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT,OOCEAN,ZEMOIST)
!
IF (KRR>0) THEN
  DO JK = IKTB+1,IKTE-1
    DO JIJ=IIJB,IIJE
      ZDTHLDZ(JIJ,JK)= 0.5*((PTHLT(JIJ,JK+D%NKL)-PTHLT(JIJ,JK    ))/PDZZ(JIJ,JK+D%NKL)+ &
                              (PTHLT(JIJ,JK    )-PTHLT(JIJ,JK-D%NKL))/PDZZ(JIJ,JK    ))
      ZDRTDZ(JIJ,JK) = 0.5*((PRT(JIJ,JK+D%NKL,1)-PRT(JIJ,JK    ,1))/PDZZ(JIJ,JK+D%NKL)+ &
                              (PRT(JIJ,JK    ,1)-PRT(JIJ,JK-D%NKL,1))/PDZZ(JIJ,JK    ))
      IF (OOCEAN) THEN
        ZVAR=CST%XG*(CST%XALPHAOC*ZDTHLDZ(JIJ,JK)-CST%XBETAOC*ZDRTDZ(JIJ,JK))
      ELSE
        ZVAR=CST%XG/PTHVREF(JIJ,JK)*                                                  &
           (ZETHETA(JIJ,JK)*ZDTHLDZ(JIJ,JK)+ZEMOIST(JIJ,JK)*ZDRTDZ(JIJ,JK))
      END IF
      !
      IF (ZVAR>0.) THEN
        PLM(JIJ,JK)=MAX(CST%XMNH_EPSILON,MIN(PLM(JIJ,JK), &
                      0.76* SQRT(PTKET(JIJ,JK)/ZVAR)))
      END IF
    END DO
  END DO
ELSE! For dry atmos or unsalted ocean runs
  DO JK = IKTB+1,IKTE-1
    DO JIJ=IIJB,IIJE
      ZDTHLDZ(JIJ,JK)= 0.5*((PTHLT(JIJ,JK+D%NKL)-PTHLT(JIJ,JK    ))/PDZZ(JIJ,JK+D%NKL)+ &
                              (PTHLT(JIJ,JK    )-PTHLT(JIJ,JK-D%NKL))/PDZZ(JIJ,JK    ))
      IF (OOCEAN) THEN
        ZVAR= CST%XG*CST%XALPHAOC*ZDTHLDZ(JIJ,JK)
      ELSE
        ZVAR= CST%XG/PTHVREF(JIJ,JK)*ZETHETA(JIJ,JK)*ZDTHLDZ(JIJ,JK)
      END IF
!
      IF (ZVAR>0.) THEN
        PLM(JIJ,JK)=MAX(CST%XMNH_EPSILON,MIN(PLM(JIJ,JK), &
                      0.76* SQRT(PTKET(JIJ,JK)/ZVAR)))
      END IF
    END DO
  END DO
END IF
!  special case near the surface
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZDTHLDZ(IIJB:IIJE,IKB)=(PTHLT(IIJB:IIJE,IKB+D%NKL)-PTHLT(IIJB:IIJE,IKB))/PDZZ(IIJB:IIJE,IKB+D%NKL)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
! For dry simulations
IF (KRR>0) THEN
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  ZDRTDZ(IIJB:IIJE,IKB)=(PRT(IIJB:IIJE,IKB+D%NKL,1)-PRT(IIJB:IIJE,IKB,1))/PDZZ(IIJB:IIJE,IKB+D%NKL)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
ELSE
  ZDRTDZ(:,IKB)=0
ENDIF
!
IF (OOCEAN) THEN
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  ZWORK2D(IIJB:IIJE)=CST%XG*(CST%XALPHAOC*ZDTHLDZ(IIJB:IIJE,IKB)-CST%XBETAOC*ZDRTDZ(IIJB:IIJE,IKB))
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
ELSE
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  ZWORK2D(IIJB:IIJE)=CST%XG/PTHVREF(IIJB:IIJE,IKB)*                                           &
              (ZETHETA(IIJB:IIJE,IKB)*ZDTHLDZ(IIJB:IIJE,IKB)+ZEMOIST(IIJB:IIJE,IKB)*ZDRTDZ(IIJB:IIJE,IKB))
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
!$mnh_expand_where(JIJ=IIJB:IIJE)
WHERE(ZWORK2D(IIJB:IIJE)>0.)
  PLM(IIJB:IIJE,IKB)=MAX(CST%XMNH_EPSILON,MIN( PLM(IIJB:IIJE,IKB),                 &
                    0.76* SQRT(PTKET(IIJB:IIJE,IKB)/ZWORK2D(IIJB:IIJE))))
END WHERE
!$mnh_end_expand_where(JIJ=IIJB:IIJE)
!
!  mixing length limited by the distance normal to the surface (with the same factor as for BL89)
!
IF (.NOT. ORMC01) THEN
  ZALPHA=0.5**(-1.5)
  !
  DO JIJ=IIJB,IIJE
    IF (OOCEAN) THEN
      DO JK=IKTE,IKTB,-1
        ZD=ZALPHA*(PZZ(JIJ,IKTE+1)-PZZ(JIJ,JK))
        IF ( PLM(JIJ,JK)>ZD) THEN
          PLM(JIJ,JK)=ZD
        ELSE
          EXIT
        ENDIF
      END DO
    ELSE
      DO JK=IKTB,IKTE
        ZD=ZALPHA*(0.5*(PZZ(JIJ,JK)+PZZ(JIJ,JK+D%NKL))-PZZ(JIJ,IKB)) &
          *PDIRCOSZW(JIJ)
        IF ( PLM(JIJ,JK)>ZD) THEN
          PLM(JIJ,JK)=ZD
        ELSE
          EXIT
        ENDIF
      END DO
    ENDIF
  END DO
END IF
!
!$mnh_expand_array(JIJ=IIJB:IIJE)
PLM(IIJB:IIJE,D%NKA) = PLM(IIJB:IIJE,IKB)
PLM(IIJB:IIJE,IKE) = PLM(IIJB:IIJE,IKE-D%NKL)
PLM(IIJB:IIJE,D%NKU) = PLM(IIJB:IIJE,D%NKU-D%NKL)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',1,ZHOOK_HANDLE2)
END SUBROUTINE DEAR
!
!     #########################
      SUBROUTINE CLOUD_MODIF_LM
!     #########################
!!
!!*****CLOUD_MODIF_LM routine to:
!!       1/ change the mixing length in the clouds
!!       2/ emphasize the mixing length in the cloud
!!           by the coefficient ZCOEF_AMPL calculated here
!!             when the CEI index is above ZCEI_MIN.
!!
!!
!!      ZCOEF_AMPL ^
!!                 |
!!                 |
!!  ZCOEF_AMPL_SAT -                       ---------- Saturation
!!    (XDUMMY1)    |                      -
!!                 |                     -
!!                 |                    -
!!                 |                   -
!!                 |                  - Amplification
!!                 |                 - straight
!!                 |                - line
!!                 |               -
!!                 |              -
!!                 |             -
!!                 |            -
!!                 |           -
!!               1 ------------
!!                 |
!!                 |
!!               0 -----------|------------|----------> PCEI
!!                 0      ZCEI_MIN     ZCEI_MAX
!!                        (XDUMMY2)    (XDUMMY3)
!!
!!
!!
!!    AUTHOR
!!    ------
!!     M. Tomasini   *CNRM METEO-FRANCE
!!
!!    MODIFICATIONS
!!    -------------
!!     Original   09/07/04
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
IMPLICIT NONE
!
!-------------------------------------------------------------------------------
!
!*       1.    INITIALISATION
!              --------------
!
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',0,ZHOOK_HANDLE2)
ZPENTE = ( PCOEF_AMPL_SAT - 1. ) / ( PCEI_MAX - PCEI_MIN )
ZCOEF_AMPL_CEI_NUL = 1. - ZPENTE * PCEI_MIN
!
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
ZCOEF_AMPL(IIJB:IIJE,1:D%NKT) = 1.
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
!*       2.    CALCULATION OF THE AMPLIFICATION COEFFICIENT
!              --------------------------------------------
!
! Saturation
!
!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
WHERE ( PCEI(IIJB:IIJE,1:D%NKT)>=PCEI_MAX ) 
  ZCOEF_AMPL(IIJB:IIJE,1:D%NKT)=PCOEF_AMPL_SAT
END WHERE
!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
! Between the min and max limits of CEI index, linear variation of the
! amplification coefficient ZCOEF_AMPL as a function of CEI
!
!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
WHERE ( PCEI(IIJB:IIJE,1:D%NKT) <  PCEI_MAX .AND. PCEI(IIJB:IIJE,1:D%NKT) >  PCEI_MIN)
  ZCOEF_AMPL(IIJB:IIJE,1:D%NKT) = ZPENTE * PCEI(IIJB:IIJE,1:D%NKT) + ZCOEF_AMPL_CEI_NUL
END WHERE
!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
!
!*       3.    CALCULATION OF THE MIXING LENGTH IN CLOUDS
!              ------------------------------------------
!
IF (HTURBLEN_CL == HTURBLEN) THEN
  ZLM_CLOUD(:,:) = ZLM(:,:)
ELSE
  SELECT CASE (HTURBLEN_CL)
!
!*         3.1 BL89 mixing length
!           ------------------
  CASE ('BL89','RM17','ADAP')
    ZSHEAR(:,:)=0.
    CALL BL89(D,CST,CSTURB,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD,OOCEAN,HPROGRAM)
!
!*         3.2 Delta mixing length
!           -------------------
  CASE ('DELT')
    CALL DELT(ZLM_CLOUD,ODZ=.TRUE.)
!
!*         3.3 Deardorff mixing length
!           -----------------------
  CASE ('DEAR')
    CALL DEAR(ZLM_CLOUD)
!
  END SELECT
ENDIF
!
!*       4.    MODIFICATION OF THE MIXING LENGTH IN THE CLOUDS
!              -----------------------------------------------
!
! Impression before modification of the mixing length
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
  TZFIELD%CMNHNAME   = 'LM_CLEAR_SKY'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'LM_CLEAR_SKY'
  TZFIELD%CUNITS     = 'm'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'X_Y_Z_LM CLEAR SKY'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZLM)
ENDIF
!
! Amplification of the mixing length when the criteria are verified
!
!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
WHERE (ZCOEF_AMPL(IIJB:IIJE,1:D%NKT) /= 1.) 
  ZLM(IIJB:IIJE,1:D%NKT) = ZCOEF_AMPL(IIJB:IIJE,1:D%NKT)*ZLM_CLOUD(IIJB:IIJE,1:D%NKT)
END WHERE
!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
! Cloud mixing length in the clouds at the points which do not verified the CEI
!
!$mnh_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
WHERE (PCEI(IIJB:IIJE,1:D%NKT) == -1.)
  ZLM(IIJB:IIJE,1:D%NKT) = ZLM_CLOUD(IIJB:IIJE,1:D%NKT)
END WHERE
!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:D%NKT)
!
!
!*       5.    IMPRESSION
!              ----------
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
  TZFIELD%CMNHNAME   = 'COEF_AMPL'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'COEF_AMPL'
  TZFIELD%CUNITS     = '1'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'X_Y_Z_COEF AMPL'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZCOEF_AMPL)
  !
  TZFIELD%CMNHNAME   = 'LM_CLOUD'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'LM_CLOUD'
  TZFIELD%CUNITS     = 'm'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'X_Y_Z_LM CLOUD'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZLM_CLOUD)
  !
ENDIF
!
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',1,ZHOOK_HANDLE2)
END SUBROUTINE CLOUD_MODIF_LM
!