turb.F90

  END IF

  IF( LBUDGET_RC ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:, :, :, 2) )
  IF( LBUDGET_RI ) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:, :, :, 4) )

  IF( LBUDGET_SV )  THEN
    DO JSV = 1, NSV
      CALL BUDGET_STORE_END( 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 

       PTP = PTP + XG / PTHVREF * MZF(PFLXZTHVMF,KKA, KKU, KKL)
       PTPMF=XG / PTHVREF * MZF(PFLXZTHVMF, KKA, KKU, KKL)

!  6.2 TKE evolution equation

IF (.NOT. LHARAT) THEN
!
IF (LBUDGET_TH)  THEN
  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
                                                          & + ZLSOCPEXNM * PRRS(:,:,:,4) )
  ELSE IF ( KRRL >= 1 ) THEN
    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
  ELSE
    CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:, :, :) )
  END IF
END IF

CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH,       &
                   & PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,            &
                   & PTSTEP,PIMPL,ZEXPL,                         &
                   & HTURBLEN,HTURBDIM,                              &
                   & TPFILE,OTURB_DIAG,           &
                   & PTP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,&
                   & TBUDGETS,KBUDGETS,&
                   & PEDR=PEDR)
IF (LBUDGET_TH)  THEN
  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) &
                                                          & + ZLSOCPEXNM * PRRS(:,:,:,4) )
  ELSE IF ( KRRL >= 1 ) THEN
    CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2) )
  ELSE
    CALL BUDGET_STORE_END( 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(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(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(TPFILE,TZFIELD,PRT(:,:,:,1))
   END IF
END IF
!
!* stores value of conservative variables & wind before turbulence tendency
PDRUS_TURB = PRUS - PDRUS_TURB
PDRVS_TURB = PRVS - PDRVS_TURB
PDRTHLS_TURB = PRTHLS - PDRTHLS_TURB
PDRRTS_TURB  = PRRS(:,:,:,1) - PDRRTS_TURB 
PDRSVS_TURB  = PRSVS - PDRSVS_TURB
!----------------------------------------------------------------------------
!
!*      8. RETRIEVE NON-CONSERVATIVE VARIABLES
!          -----------------------------------
!
IF ( KRRL >= 1 ) THEN
  IF ( KRRI >= 1 ) THEN
    PRT(:,:,:,1)  = PRT(:,:,:,1)  - PRT(:,:,:,2)  - PRT(:,:,:,4)
    PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2) - PRRS(:,:,:,4)
    PTHLT(:,:,:)  = PTHLT(:,:,:)  + ZLVOCPEXNM(:,:,:) * PRT(:,:,:,2) &
                                  + ZLSOCPEXNM(:,:,:) * PRT(:,:,:,4)
    PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLVOCPEXNM(:,:,:) * PRRS(:,:,:,2) &
                                  + ZLSOCPEXNM(:,:,:) * PRRS(:,:,:,4)
!
    DEALLOCATE(ZLVOCPEXNM)
    DEALLOCATE(ZLSOCPEXNM)
  ELSE
    PRT(:,:,:,1)  = PRT(:,:,:,1)  - PRT(:,:,:,2) 
    PRRS(:,:,:,1) = PRRS(:,:,:,1) - PRRS(:,:,:,2)
    PTHLT(:,:,:)  = PTHLT(:,:,:)  + ZLOCPEXNM(:,:,:) * PRT(:,:,:,2)
    PRTHLS(:,:,:) = PRTHLS(:,:,:) + ZLOCPEXNM(:,:,:) * PRRS(:,:,:,2)
  END IF
END IF

! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
!CALL SOURCES_NEG_CORRECT(HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
!----------------------------------------------------------------------------
!
!*      9. LES averaged surface fluxes
!          ---------------------------
!
IF (LLES_CALL) THEN
  CALL SECOND_MNH(ZTIME1)
  CALL LES_MEAN_SUBGRID(PSFTH,X_LES_Q0)
  CALL LES_MEAN_SUBGRID(PSFRV,X_LES_E0)
  DO JSV=1,NSV
    CALL LES_MEAN_SUBGRID(PSFSV(:,:,JSV),X_LES_SV0(:,JSV))
  END DO
  CALL LES_MEAN_SUBGRID(PSFU,X_LES_UW0)
  CALL LES_MEAN_SUBGRID(PSFV,X_LES_VW0)
  CALL LES_MEAN_SUBGRID((PSFU*PSFU+PSFV*PSFV)**0.25,X_LES_USTAR)
!----------------------------------------------------------------------------
!
!*     10. LES for 3rd order moments
!          -------------------------
!
  CALL LES_MEAN_SUBGRID(ZMWTH,X_LES_SUBGRID_W2Thl)
  CALL LES_MEAN_SUBGRID(ZMTH2,X_LES_SUBGRID_WThl2)
  IF (KRR>0) THEN
    CALL LES_MEAN_SUBGRID(ZMWR,X_LES_SUBGRID_W2Rt)
    CALL LES_MEAN_SUBGRID(ZMTHR,X_LES_SUBGRID_WThlRt)
    CALL LES_MEAN_SUBGRID(ZMR2,X_LES_SUBGRID_WRt2)
  END IF
!
!----------------------------------------------------------------------------
!
!*     11. LES quantities depending on <w'2> in "1DIM" mode
!          ------------------------------------------------
!
  IF (HTURBDIM=="1DIM") THEN
    CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_U2)
    X_LES_SUBGRID_V2 = X_LES_SUBGRID_U2
    X_LES_SUBGRID_W2 = X_LES_SUBGRID_U2
    CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PTHLT,PDZZ),&
                          KKA, KKU, KKL),X_LES_RES_ddz_Thl_SBG_W2)
    IF (KRR>=1) &
    CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PRT(:,:,:,1),PDZZ),&
                         &KKA, KKU, KKL),X_LES_RES_ddz_Rt_SBG_W2)
    DO JSV=1,NSV
      CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PSVT(:,:,:,JSV),PDZZ), &
                           &KKA, KKU, KKL), X_LES_RES_ddz_Sv_SBG_W2(:,:,:,JSV))
    END DO
  END IF

!----------------------------------------------------------------------------
!
!*     12. LES mixing end dissipative lengths, presso-correlations
!          -------------------------------------------------------
!
  CALL LES_MEAN_SUBGRID(ZLM,X_LES_SUBGRID_LMix)
  CALL LES_MEAN_SUBGRID(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(ZLEPS,X_LES_SUBGRID_WP)
!
  CALL SECOND_MNH(ZTIME2)
  XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
END IF
!
IF(PRESENT(PLEM)) PLEM = ZLM
!----------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
CONTAINS
!
!
!     ##############################################
      SUBROUTINE UPDATE_ROTATE_WIND(PUSLOPE,PVSLOPE)
!     ##############################################
!!
!!****  *UPDATE_ROTATE_WIND* routine to set rotate wind values at the border
!
!!    AUTHOR
!!    ------
!!
!!     P Jabouille   *CNRM METEO-FRANCE
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   24/06/99
!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
USE MODE_ll
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
USE MODD_CONF
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: PUSLOPE,PVSLOPE
! tangential surface fluxes in the axes following the orography
!
!*       0.2   Declarations of local variables :
!
INTEGER             :: IIB,IIE,IJB,IJE,IIU,IJU ! index values for the physical subdomain
TYPE(LIST_ll), POINTER :: TZFIELDS_ll  ! list of fields to exchange
INTEGER                :: IINFO_ll     ! return code of parallel routine
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',0,ZHOOK_HANDLE)
!
!*        1  PROLOGUE
!
NULLIFY(TZFIELDS_ll)
!
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
!
!         2 Update halo if necessary
!
!!$IF (NHALO == 1) THEN
  CALL ADD2DFIELD_ll( TZFIELDS_ll, PUSLOPE, 'UPDATE_ROTATE_WIND::PUSLOPE' )
  CALL ADD2DFIELD_ll( TZFIELDS_ll, PVSLOPE, 'UPDATE_ROTATE_WIND::PVSLOPE' )
  CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
  CALL CLEANLIST_ll(TZFIELDS_ll)
!!$ENDIF
!
!        3 Boundary conditions for non cyclic case
!
IF ( HLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
  PUSLOPE(IIB-1,:)=PUSLOPE(IIB,:)
  PVSLOPE(IIB-1,:)=PVSLOPE(IIB,:)
END IF
IF ( HLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
  PUSLOPE(IIE+1,:)=PUSLOPE(IIE,:)
  PVSLOPE(IIE+1,:)=PVSLOPE(IIE,:)
END IF
IF ( HLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
  PUSLOPE(:,IJB-1)=PUSLOPE(:,IJB)
  PVSLOPE(:,IJB-1)=PVSLOPE(:,IJB)
END IF
IF(  HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
  PUSLOPE(:,IJE+1)=PUSLOPE(:,IJE)
  PVSLOPE(:,IJE+1)=PVSLOPE(:,IJE)
END IF
!
IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',1,ZHOOK_HANDLE)
!
END SUBROUTINE UPDATE_ROTATE_WIND
!
!     ########################################################################
      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
!              ------------
USE MODD_CST
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments 
!
REAL,                   INTENT(IN)    :: PALP,PBETA,PGAM,PLTT,PC
REAL, DIMENSION(:,:,:), INTENT(IN)    :: PT,PEXN,PCP
!
REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLOCPEXN
REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PAMOIST,PATHETA
! 
!*       0.2   Declarations of local variables
!
REAL                :: ZEPS         ! XMV / XMD
REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZRVSAT
REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
!
!-------------------------------------------------------------------------------
!
  REAL(KIND=JPRB) :: ZHOOK_HANDLE
  IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',0,ZHOOK_HANDLE)
  ZEPS = XMV / XMD
!
!*       1.1 Lv/Cph at  t
!
  PLOCPEXN(:,:,:) = ( PLTT + (XCPV-PC) *  (PT(:,:,:)-XTT) ) / PCP(:,:,:)
!
!*      1.2 Saturation vapor pressure at t
!
  ZRVSAT(:,:,:) =  EXP( PALP - PBETA/PT(:,:,:) - PGAM*ALOG( PT(:,:,:) ) )
!
!*      1.3 saturation  mixing ratio at t
!
  ZRVSAT(:,:,:) =  ZRVSAT(:,:,:) * ZEPS / ( PPABST(:,:,:) - ZRVSAT(:,:,:) )
!
!*      1.4 compute the saturation mixing ratio derivative (rvs')
!
  ZDRVSATDT(:,:,:) = ( PBETA / PT(:,:,:)  - PGAM ) / PT(:,:,:)   &
                 * ZRVSAT(:,:,:) * ( 1. + ZRVSAT(:,:,:) / ZEPS )
!
!*      1.5 compute Amoist
!
  PAMOIST(:,:,:)=  0.5 / ( 1.0 + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) )
!
!*      1.6 compute Atheta
!
  PATHETA(:,:,:)= PAMOIST(:,:,:) * PEXN(:,:,:) *                             &
        ( ( ZRVSAT(:,:,:) - PRT(:,:,:,1) ) * PLOCPEXN(:,:,:) /               &
          ( 1. + ZDRVSATDT(:,:,:) * PLOCPEXN(:,:,:) )        *               &
          (                                                                  &
           ZRVSAT(:,:,:) * (1. + ZRVSAT(:,:,:)/ZEPS)                         &
                        * ( -2.*PBETA/PT(:,:,:) + PGAM ) / PT(:,:,:)**2      &
          +ZDRVSATDT(:,:,:) * (1. + 2. * ZRVSAT(:,:,:)/ZEPS)                 &
                        * ( PBETA/PT(:,:,:) - PGAM ) / PT(:,:,:)             &
          )                                                                  &
         - ZDRVSATDT(:,:,:)                                                  &
        )
!
!*      1.7 Lv/Cph/Exner at t-1
!
  PLOCPEXN(:,:,:) = PLOCPEXN(:,:,:) / PEXN(:,:,:)
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE)
END SUBROUTINE COMPUTE_FUNCTION_THERMO
!
!     ####################
      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(:,:,:), INTENT(OUT)   :: PLM
LOGICAL,                INTENT(IN)    :: ODZ
!
!*       0.2   Declarations of local variables
!
REAL                :: ZD           ! distance to the surface
!
!-------------------------------------------------------------------------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:DELT',0,ZHOOK_HANDLE)
IF (ODZ) THEN
  ! Dz is take into account in the computation
  DO JK = IKTB,IKTE ! 1D turbulence scheme
    PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK)
  END DO
  PLM(:,:,KKU) = PLM(:,:,IKE)
  PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
  IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
    IF ( L2D) THEN
      PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) 
    ELSE
      PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
    END IF
  END IF
ELSE
  ! Dz not taken into account in computation to assure invariability with vertical grid mesh
  PLM=1.E10
  IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
    IF ( L2D) THEN
      PLM(:,:,:) = MXF(PDXX(:,:,:))
    ELSE
      PLM(:,:,:) = (MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./2.)
    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 JJ=1,SIZE(PUT,2)
    DO JI=1,SIZE(PUT,1)
      IF (LOCEAN) THEN
        DO JK=IKTE,IKTB,-1
          ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
          IF ( PLM(JI,JJ,JK)>ZD) THEN
            PLM(JI,JJ,JK)=ZD
          ELSE
            EXIT
          ENDIF
       END DO
      ELSE
        DO JK=IKTB,IKTE
          ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))&
          -PZZ(JI,JJ,IKB)) *PDIRCOSZW(JI,JJ)
          IF ( PLM(JI,JJ,JK)>ZD) THEN
            PLM(JI,JJ,JK)=ZD
          ELSE
            EXIT
          ENDIF
        END DO
      ENDIF   
    END DO
  END DO
END IF
!
PLM(:,:,KKA) = PLM(:,:,IKB  )
PLM(:,:,KKU  ) = PLM(:,:,IKE)
!
IF (LHOOK) CALL DR_HOOK('TURB:DELT',1,ZHOOK_HANDLE)
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(:,:,:), INTENT(OUT)   :: PLM
!
!*       0.2   Declarations of local variables
!
REAL                :: ZD           ! distance to the surface
REAL                :: ZVAR         ! Intermediary variable
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) ::   ZWORK2D
!
REAL, DIMENSION(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)) ::     &
            ZDTHLDZ,ZDRTDZ,     &!dtheta_l/dz, drt_dz used for computing the stablity
!                                ! criterion 
            ZETHETA,ZEMOIST             !coef ETHETA and EMOIST
!----------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
!   initialize the mixing length with the mesh grid
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',0,ZHOOK_HANDLE)
! 1D turbulence scheme
PLM(:,:,IKTB:IKTE) = PZZ(:,:,IKTB+KKL:IKTE+KKL) - PZZ(:,:,IKTB:IKTE)
PLM(:,:,KKU) = PLM(:,:,IKE)
PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA)
IF ( HTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
  IF ( L2D) THEN
    PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) )
  ELSE
    PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.)
  END IF
END IF
!   compute a mixing length limited by the stability
!
ZETHETA(:,:,:) = ETHETA(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT)
ZEMOIST(:,:,:) = EMOIST(KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT)
!
IF (KRR>0) THEN
  DO JK = IKTB+1,IKTE-1
    DO JJ=1,SIZE(PUT,2)
      DO JI=1,SIZE(PUT,1)
        ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK    ))/PDZZ(JI,JJ,JK+KKL)+ &
                                (PTHLT(JI,JJ,JK    )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK    ))
        ZDRTDZ(JI,JJ,JK) = 0.5*((PRT(JI,JJ,JK+KKL,1)-PRT(JI,JJ,JK    ,1))/PDZZ(JI,JJ,JK+KKL)+ &
                                (PRT(JI,JJ,JK    ,1)-PRT(JI,JJ,JK-KKL,1))/PDZZ(JI,JJ,JK    ))
        IF (LOCEAN) THEN
          ZVAR=XG*(XALPHAOC*ZDTHLDZ(JI,JJ,JK)-XBETAOC*ZDRTDZ(JI,JJ,JK))
        ELSE
          ZVAR=XG/PTHVREF(JI,JJ,JK)*                                                  &
             (ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)+ZEMOIST(JI,JJ,JK)*ZDRTDZ(JI,JJ,JK))
        END IF
        !
        IF (ZVAR>0.) THEN
          PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
                        0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
        END IF
      END DO
    END DO
  END DO
ELSE! For dry atmos or unsalted ocean runs
  DO JK = IKTB+1,IKTE-1
    DO JJ=1,SIZE(PUT,2)
      DO JI=1,SIZE(PUT,1)
        ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK    ))/PDZZ(JI,JJ,JK+KKL)+ &
                                (PTHLT(JI,JJ,JK    )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK    ))
        IF (LOCEAN) THEN
          ZVAR= XG*XALPHAOC*ZDTHLDZ(JI,JJ,JK)
        ELSE
          ZVAR= XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)
        END IF
!
        IF (ZVAR>0.) THEN
          PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), &
                        0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR)))
        END IF
      END DO
    END DO
  END DO
END IF
!  special case near the surface 
ZDTHLDZ(:,:,IKB)=(PTHLT(:,:,IKB+KKL)-PTHLT(:,:,IKB))/PDZZ(:,:,IKB+KKL)
! For dry simulations
IF (KRR>0) THEN
  ZDRTDZ(:,:,IKB)=(PRT(:,:,IKB+KKL,1)-PRT(:,:,IKB,1))/PDZZ(:,:,IKB+KKL)
ELSE
  ZDRTDZ(:,:,IKB)=0
ENDIF
!
IF (LOCEAN) THEN
  ZWORK2D(:,:)=XG*(XALPHAOC*ZDTHLDZ(:,:,IKB)-XBETAOC*ZDRTDZ(:,:,IKB))
ELSE
  ZWORK2D(:,:)=XG/PTHVREF(:,:,IKB)*                                           &
              (ZETHETA(:,:,IKB)*ZDTHLDZ(:,:,IKB)+ZEMOIST(:,:,IKB)*ZDRTDZ(:,:,IKB))
END IF
WHERE(ZWORK2D(:,:)>0.)
  PLM(:,:,IKB)=MAX(XMNH_EPSILON,MIN( PLM(:,:,IKB),                 &
                    0.76* SQRT(PTKET(:,:,IKB)/ZWORK2D(:,:))))
END WHERE
!
!  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 JJ=1,SIZE(PUT,2)
    DO JI=1,SIZE(PUT,1)
      IF (LOCEAN) THEN
        DO JK=IKTE,IKTB,-1
          ZD=ZALPHA*(PZZ(JI,JJ,IKTE+1)-PZZ(JI,JJ,JK))
          IF ( PLM(JI,JJ,JK)>ZD) THEN
            PLM(JI,JJ,JK)=ZD
          ELSE
            EXIT
          ENDIF
        END DO
      ELSE
        DO JK=IKTB,IKTE
          ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,IKB)) &
            *PDIRCOSZW(JI,JJ)
          IF ( PLM(JI,JJ,JK)>ZD) THEN
            PLM(JI,JJ,JK)=ZD
          ELSE
            EXIT
          ENDIF
        END DO
      ENDIF 
    END DO
  END DO
END IF
!
PLM(:,:,KKA) = PLM(:,:,IKB  )
PLM(:,:,IKE  ) = PLM(:,:,IKE-KKL)
PLM(:,:,KKU  ) = PLM(:,:,KKU-KKL)
!
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',1,ZHOOK_HANDLE)
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
!
REAL :: ZPENTE            ! Slope of the amplification straight line
REAL :: ZCOEF_AMPL_CEI_NUL! Ordonnate at the origin of the
                          ! amplification straight line
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCOEF_AMPL
                          ! Amplification coefficient of the mixing length
                          ! when the instability criterium is verified 
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZLM_CLOUD
                          ! Turbulent mixing length in the clouds
!
!-------------------------------------------------------------------------------
!
!*       1.    INITIALISATION
!              --------------
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',0,ZHOOK_HANDLE)
ZPENTE = ( PCOEF_AMPL_SAT - 1. ) / ( PCEI_MAX - PCEI_MIN ) 
ZCOEF_AMPL_CEI_NUL = 1. - ZPENTE * PCEI_MIN
!
ZCOEF_AMPL(:,:,:) = 1.
!
!*       2.    CALCULATION OF THE AMPLIFICATION COEFFICIENT
!              --------------------------------------------
!
! Saturation
!
WHERE ( PCEI(:,:,:)>=PCEI_MAX ) ZCOEF_AMPL(:,:,:)=PCOEF_AMPL_SAT
!
! Between the min and max limits of CEI index, linear variation of the
! amplification coefficient ZCOEF_AMPL as a function of CEI
!
WHERE ( PCEI(:,:,:) <  PCEI_MAX .AND.                                        &
        PCEI(:,:,:) >  PCEI_MIN      )                                       &
        ZCOEF_AMPL(:,:,:) = ZPENTE * PCEI(:,:,:) + ZCOEF_AMPL_CEI_NUL  
!
!
!*       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(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD)
!
!*         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(TPFILE,TZFIELD,ZLM)
ENDIF
!
! Amplification of the mixing length when the criteria are verified
!
WHERE (ZCOEF_AMPL(:,:,:) /= 1.) ZLM(:,:,:) = ZCOEF_AMPL(:,:,:)*ZLM_CLOUD(:,:,:)
!
! Cloud mixing length in the clouds at the points which do not verified the CEI
!
WHERE (PCEI(:,:,:) == -1.) ZLM(:,:,:) = ZLM_CLOUD(:,:,:)
!
!
!*       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(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(TPFILE,TZFIELD,ZLM_CLOUD)
  !
ENDIF
!
IF (LHOOK) CALL DR_HOOK('TURB:CLOUD_MODIF_LM',1,ZHOOK_HANDLE)
END SUBROUTINE CLOUD_MODIF_LM
!
END SUBROUTINE TURB