     MODULE MODULE_BL_MFSHCONVPBL




USE MODULE_MODEL_CONSTANTS

REAL,PARAMETER :: XG      = 9.80665

REAL,PARAMETER :: XP00= 1.E5               ! Reference pressure
!
!
REAL,PARAMETER :: XMD= 28.9644E-3
REAL,PARAMETER :: XMV= 18.0153E-3            ! Molar mass of dry air and molar mass of vapor
REAL,PARAMETER :: XRD=R_D
REAL,PARAMETER :: XRV=R_V        ! Gaz constant for dry air, gaz constant for vapor
REAL,PARAMETER :: XCPD=7.* XRD /2.
REAL,PARAMETER :: XCPV=4.* XRV          ! Cpd (dry air), Cpv (vapor)
REAL,PARAMETER :: XCL= 4.218E+3
REAL,PARAMETER :: XCI= 2.106E+3      ! Cl (liquid), Ci (ice)
REAL,PARAMETER :: XTT= 273.16        ! Triple point temperature
REAL,PARAMETER :: XLVTT=2.5008E+6              ! Vaporization heat constant
REAL,PARAMETER :: XLSTT=2.8345E+6   ! Sublimation heat constant
                               !  temperature
REAL,PARAMETER :: XGAMW = (XCL - XCPV) / XRV! Constants for saturation vapor
REAL,PARAMETER :: XBETAW= (XLVTT/XRV) + (XGAMW * XTT)
!The use of intrinsics in an initialization expressions is a F2003 feature.
!For backward compatibility, hard coded Log(644.11) & Log(XTT) here
!REAL,PARAMETER :: XALPW= LOG(611.14) + (XBETAW /XTT) + (XGAMW *LOG(XTT))
REAL,PARAMETER :: LOG_611_14 = 6.415326
REAL,PARAMETER :: LOG_XTT    = 5.610058
REAL,PARAMETER :: XALPW= LOG_611_14 + (XBETAW /XTT) + (XGAMW *LOG_XTT)
                                !  pressure  function
REAL,PARAMETER :: XGAMI  = (XCI - XCPV) / XRV
REAL,PARAMETER :: XBETAI = (XLSTT/XRV) + (XGAMI * XTT)
!REAL,PARAMETER :: XALPI  = LOG(611.14) + (XBETAI /XTT) + (XGAMI *LOG(XTT))
REAL,PARAMETER :: XALPI  = LOG_611_14 + (XBETAI /XTT) + (XGAMI *LOG_XTT)
REAL,PARAMETER :: XLINI = 0.32


REAL, PARAMETER :: XALP_PERT   = 0.3  ! coefficient for the perturbation of
                   ! theta_l and r_t at the first level of
                   ! the updraft
REAL, PARAMETER ::XABUO       = 1.   ! coefficient of the buoyancy term in the w_up equation
REAL, PARAMETER ::XBENTR      = 1.   ! coefficient of the entrainment term in thew_up equation
REAL, PARAMETER ::XBDETR      = 0.   ! coefficient of the detrainment term in thew_up equation


REAL, PARAMETER ::XCMF        = 0.065! coefficient for the mass flux at the firstlevel 0.065
                                ! of the updraft (closure) XCMF        = 0.065
REAL, PARAMETER ::XENTR_DRY   = 0.55 ! coefficient for entrainment in dry part XENTR_DRY   = 0.55
REAL, PARAMETER ::XDETR_DRY   = 10.  ! coefficient for detrainment in dry part XDETR_DRY   = 10.
REAL, PARAMETER ::XDETR_LUP   = 1.0   !  coefficient for detrainment in dry part XDETR_LUP   = 1.

REAL, PARAMETER ::XENTR_MF    = 0.035! entrainment constant (m/Pa) = 0.2 (m)  XENTR_MF    = 0.035
REAL, PARAMETER ::XCRAD_MF    = 50.  ! cloud radius in cloudy part
REAL, PARAMETER ::XKCF_MF     = 2.75 ! coefficient for cloud fraction
REAL, PARAMETER ::XKRC_MF     = 1.   ! coefficient for convective rc
REAL, PARAMETER ::XTAUSIGMF   = 600.
REAL, PARAMETER ::XPRES_UV    = 0.5  ! coefficient for pressure term in wind mixing
!
REAL, PARAMETER ::XFRAC_UP_MAX= 0.33 ! maximum Updraft fraction


!

     CONTAINS


      SUBROUTINE MFSHCONVPBL (DT,STEPBL,HT,DZ           &
                    ,RHO,PMID,PINT,TH,EXNER             &
                    ,QV, QC, U, V                       &
                    ,HFX, QFX, TKE                      &
                    ,RUBLTEN,RVBLTEN,RTHBLTEN           &
                    ,RQVBLTEN,RQCBLTEN                  &
                    ,IDS,IDE,JDS,JDE,KDS,KDE            &
                    ,IMS,IME,JMS,JME,KMS,KME            &
                    ,ITS,ITE,JTS,JTE,KTS,KTE,KRR        &
                    ,MASSFLUX_EDKF, ENTR_EDKF, DETR_EDKF & 
                    ,THL_UP, THV_UP, RT_UP, RV_UP       &
                    ,RC_UP, U_UP, V_UP, FRAC_UP, RC_MF  &
                    ,WTHV,PLM_BL89 )                   
 
      IMPLICIT NONE
!
!----------------------------------------------------------------------
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
     &                     ,IMS,IME,JMS,JME,KMS,KME                    &
     &                     ,ITS,ITE,JTS,JTE,KTS,KTE
!
      INTEGER,INTENT(IN) :: KRR

      INTEGER,INTENT(IN) :: STEPBL

      REAL,INTENT(IN) :: DT
!
      REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT, HFX, QFX
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ          &
     &                                                     ,EXNER       &
     &                                                     ,PMID,PINT   &
     &                                                     ,QV,QC,RHO   &
     &                                                     ,TH,U,V,TKE   
!
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: RQCBLTEN,RQVBLTEN       &
     &                                        ,RTHBLTEN                                &
     &                                        ,RUBLTEN,RVBLTEN        
  
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),OPTIONAL,INTENT(OUT) ::              &
     &                                         MASSFLUX_EDKF, ENTR_EDKF, DETR_EDKF &
     &                                        ,THL_UP, THV_UP, RT_UP, RV_UP        &
     &                                        ,RC_UP, U_UP, V_UP, FRAC_UP

      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),OPTIONAL,INTENT(INOUT) ::            &
     &                                         RC_MF

      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) ::   WTHV
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) ::   PLM_BL89

!
!Local declaration

INTEGER   :: KRRL      ! number of liquid water var.
INTEGER   :: KRRI      ! number of ice water var.
LOGICAL   :: OMIXUV    ! True if mixing of momentum

REAL   :: PIMPL_MF     ! degre of implicitness
REAL   ::  PTSTEP   ! Dynamical timestep 
REAL   ::  PTSTEP_MET! Timestep for meteorological variables                        

REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)   :: PZZ       !  Height at the flux point

REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)   :: PZZM      !  Height at the mass point

REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)   :: PDZZ      ! depth between mass levels 
 
REAL, DIMENSION(ITS:ITE,JTS:JTE)   ::  PSFTH,PSFRV
                                            ! normal surface fluxes of theta,rv
!
!    prognostic variables at t- deltat
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PPABSM      ! Pressure at mass point
!REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PPABSF    ! Pressure at flux point
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PEXNM       ! Exner function at t-dt
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PRHODREF    ! dry density of the
                                                     ! reference state
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PRHODJ    ! dry density of the

REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PTKEM       ! TKE
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PUM,PVM     ! momentum

!   thermodynamical variables which are transformed in conservative var.
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     ::  PTHM       ! pot. temp. = PTHLM in turb.f90
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE,KRR) ::  PRM        ! water species 

REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     ::  PRUS,PRVS,PRTHS 
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE,KRR) ::  PRRS

! For diagnostic output
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PEMF, PENTR, PDETR
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PTHL_UP, PRT_UP, PRV_UP, PRC_UP
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PU_UP, PV_UP, PTHV_UP, PFRAC_UP                   
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PRC_MF
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: WTHV_MF
REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PLM_MF

INTEGER :: I,J,K  ! loop variables
      
! Transform WRF Variable to input of mass flux scheme

DO J=JTS,JTE
   DO K=KTS,KTE
      DO I=ITS,ITE
        IF (K==KTS) PZZ(I,J,K)=0.
        PEMF(I,J,K)=0.
        PENTR(I,J,K)=0.
        PDETR(I,J,K)=0.
        PTHL_UP(I,J,K)=0.
        PTHV_UP(I,J,K)=0.
        PRT_UP(I,J,K)=0.
        PRV_UP(I,J,K)=0.
        PRC_UP(I,J,K)=0.
        PU_UP(I,J,K)=0.
        PV_UP(I,J,K)=0.
        PFRAC_UP(I,J,K)=0.
        WTHV_MF(I,J,K)=0.
        PTHM(I,J,K)=TH(I,K,J)
        PTKEM(I,J,K)=TKE(I,K,J)
        PRM(I,J,K,1)=QV(I,K,J)-RC_MF(I,K,J)
        PRM(I,J,K,2)=RC_MF(I,K,J)
        PUM(I,J,K)=U(I,K,J)
        PVM(I,J,K)=V(I,K,J)
        PRHODREF(I,J,K)=RHO(I,K,J)/(1.+QV(I,K,J))
        PEXNM(I,J,K)=EXNER(I,K,J)
        PPABSM(I,J,K)=PMID(I,K,J)
        IF (K/=KTE) THEN
            PZZ(I,J,K+1)=PZZ(I,J,K)+DZ(I,K,J)
            PZZM(I,J,K)=0.5*(PZZ(I,J,K+1)+PZZ(I,J,K)) ! z at mass point
        ELSE
            PZZM(I,J,K)=PZZ(I,J,K)+0.5*DZ(I,K-1,J) ! z at mass point
        ENDIF
        IF (K==KTS) THEN
           PDZZ(I,J,K)=2*(PZZM(I,J,K))      
        ELSE
           PDZZ(I,J,K)=PZZM(I,J,K)-PZZM(I,J,K-1)      
        ENDIF
     
        PRHODJ(I,J,K)=PRHODREF(I,J,K)*DZ(I,K,J)
  

      ENDDO
   ENDDO
ENDDO

! fill the kte+1 level
PTHM(:,:,KTE)=PTHM(:,:,KTE-1)
PTKEM(:,:,KTE)=PTKEM(:,:,KTE-1)
PRM(:,:,KTE,1)=PRM(:,:,KTE-1,1)
PRM(:,:,KTE,2)=PRM(:,:,KTE-1,2)
PUM(:,:,KTE)=PUM(:,:,KTE-1)
PVM(:,:,KTE)=PVM(:,:,KTE-1)
PRHODREF(:,:,KTE)=PRHODREF(:,:,KTE-1)
PEXNM(:,:,KTE)=PEXNM(:,:,KTE-1)
PPABSM(:,:,KTE)=PPABSM(:,:,KTE-1)
PRHODJ(:,:,KTE)=PRHODJ(:,:,KTE-1)

PSFTH(:,:)=HFX(ITS:ITE,JTS:JTE)/(PRHODREF(:,:,KTS)*XCPD)
PSFRV(:,:)=QFX(ITS:ITE,JTS:JTE)/(PRHODREF(:,:,KTS))



! Assign some variables

OMIXUV=.FALSE.
KRRL=1 !Qc is managed
KRRI=0 !Qi not
PIMPL_MF=0.
PTSTEP=DT*STEPBL
PTSTEP_MET=PTSTEP

      CALL MFSHCONVPBL_CORE(KRR,KRRL,KRRI,                            &
                OMIXUV,                                               &
                PIMPL_MF,PTSTEP,PTSTEP_MET,                           &
                PDZZ, PZZ,                                            &
                PRHODJ, PRHODREF,                                     &
                PPABSM, PEXNM,                                        &
                PSFTH,PSFRV,                                          &
                PTHM,PRM,PUM,PVM,PTKEM,                               &
                PRTHS,PRRS,PRUS,PRVS,PEMF, PENTR, PDETR,              &
                PTHL_UP, PRT_UP, PRV_UP, PRC_UP,                      &
                PU_UP, PV_UP, PTHV_UP, PFRAC_UP, PRC_MF, WTHV_MF,PLM_MF  )


DO J=JTS,JTE
   DO K=KTS,KTE
      DO I=ITS,ITE
         RQCBLTEN(I,K,J)=PRRS(I,J,K,2)       
         RQVBLTEN(I,K,J)=PRRS(I,J,K,1)       
         RTHBLTEN(I,K,J)=PRTHS(I,J,K)                
         RUBLTEN(I,K,J)=PRUS(I,J,K)
         RVBLTEN(I,K,J)=PRVS(I,J,K)
         WTHV(I,K,J)=WTHV_MF(I,J,K)
         PLM_BL89(I,K,J)=PLM_MF(I,J,K)
      ENDDO
   ENDDO
ENDDO

      IF ( PRESENT(MASSFLUX_EDKF) ) THEN
         DO J=JTS,JTE
            DO K=KTS,KTE
               DO I=ITS,ITE
                  MASSFLUX_EDKF(I,K,J)=PEMF(I,J,K)
                  ENTR_EDKF(I,K,J)=PENTR(I,J,K)
                  DETR_EDKF(I,K,J)=PDETR(I,J,K)
                  THL_UP(I,K,J)=PTHL_UP(I,J,K)
                  THV_UP(I,K,J)=PTHV_UP(I,J,K)
                  RT_UP(I,K,J)=PRT_UP(I,J,K)
                  RV_UP(I,K,J)=PRV_UP(I,J,K)
                  RC_UP(I,K,J)=PRC_UP(I,J,K)
                  U_UP(I,K,J)=PU_UP(I,J,K)
                  V_UP(I,K,J)=PV_UP(I,J,K)
                  FRAC_UP(I,K,J)=PFRAC_UP(I,J,K)
                  RC_MF(I,K,J)=PRC_MF(I,J,K)
               ENDDO
            ENDDO
         ENDDO
      ENDIF


END SUBROUTINE MFSHCONVPBL

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   WRAPPER from WRF to MASS FLUX SCHEME 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
  
      SUBROUTINE MFSHCONVPBL_CORE(KRR,KRRL,KRRI,                            &
                OMIXUV,                       &
                PIMPL_MF,PTSTEP,PTSTEP_MET,                 &
                PDZZ, PZZ,                                            &
                PRHODJ, PRHODREF,                                     &
                PPABSM, PEXNM,                                        &
                PSFTH,PSFRV,                                          &
                PTHM,PRM,PUM,PVM,PTKEM,                               &
                PRTHS,PRRS,PRUS,PRVS,  PEMF, PENTR, PDETR,            &
                PTHL_UP, PRT_UP, PRV_UP, PRC_UP,                      &
                PU_UP, PV_UP, PTHV_UP, PFRAC_UP,  PRC_MF,             &
                PFLXZTHVMF,PLM       )
!!
!!****  *MFSHCONVPBL_CORE* - Interfacing routine
!!

!! --------------------------------------------------------------------------
!


IMPLICIT NONE

INTEGER,                INTENT(IN)   :: KRR          ! number of moist var.
INTEGER,                INTENT(IN)   :: KRRL         ! number of liquid water var.
INTEGER,                INTENT(IN)   :: KRRI         ! number of ice water var.
LOGICAL,                INTENT(IN)   :: OMIXUV    ! True if mixing of momentum
             
REAL,                   INTENT(IN)   :: PIMPL_MF     ! degre of implicitness
REAL,                 INTENT(IN)     ::  PTSTEP   ! Dynamical timestep 
REAL,                 INTENT(IN)     ::  PTSTEP_MET! Timestep for meteorological variables                        
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PZZ         ! Height of flux point
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PDZZ        ! Metric coefficients
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PRHODJ      ! dry density * Grid size
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PRHODREF    ! dry density of the
                                                     ! reference state
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PEXNM       ! Exner function at t-dt

REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTH,PSFRV ! normal surface fluxes of theta and Rv 
REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHM        ! Theta at t-dt
REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM         ! water var. at t-dt
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUM,PVM     ! wind components at t-dt
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! tke at t-dt


REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PRUS,PRVS,PRTHS ! Meso-NH sources
REAL, DIMENSION(:,:,:,:), INTENT(OUT) ::  PRRS 

REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PEMF, PENTR, PDETR
REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PTHL_UP, PRT_UP, PRV_UP, PRC_UP 
REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PU_UP, PV_UP, PTHV_UP, PFRAC_UP
REAL, DIMENSION(:,:,:), INTENT(INOUT) ::  PRC_MF 
REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PFLXZTHVMF 
REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PLM
  

!
!                     0.2  Declaration of local variables
!
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2),SIZE(PTHM,3)) ::     &
          ZEXN,ZCPH,                                      &
          PRV,PRL,PTH,                                    &
          ZTM,                                            &  ! Temperature at t-dt
          ZLVOCPEXN,                                      &  !
          ZCF_MF,                                         & 
          ZLSOCPEXN,                                      &  !
          ZAMOIST,                                        &  !
          ZATHETA,                                        &  !
          ZTHLM,                                          &  !
          ZRTM,                                           &  !
          ZTHVM,ZTHVREF,ZUMM,ZVMM,                        &  !
          ZRI_UP,ZW_UP,                                   & ! 
          ZEMF_O_RHODREF,                                 & ! entrainment/detrainment
          ZTHLDT,ZRTDT,ZUDT,ZVDT,                         & ! tendencies
          ZFLXZTHMF,ZFLXZRMF,ZFLXZUMF,ZFLXZVMF   ! fluxes


INTEGER,DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2))  :: IKLCL,IKETL,IKCTL 
                    
INTEGER :: IKU, IKB, IKE  
INTEGER :: JI,JJ,JK,JSV                          ! Loop counters
INTEGER :: IRESP    ! error code


!------------------------------------------------------------------------

!!! 1. Initialisation

! Internal Domain
IKU=SIZE(PTHM,3)
IKB=1              ! Modif WRF JP
IKE=IKU-1


! number of scalar var

  ZUMM=PUM   !Modif WRF JP
  ZVMM=PVM   !Modif WRF JP

! Thermodynamics functions

  CALL COMPUTE_FUNCTION_THERMO_MF( KRR,KRRL,KRRI,                   &
                                   PTHM,PRM,PEXNM,PPABSM,           &
                                   ZTM,ZLVOCPEXN,ZLSOCPEXN,         &
                                   ZAMOIST,ZATHETA                  )



! Conservative variables at t-dt

  CALL THL_RT_FROM_TH_R_MF( KRR,KRRL,KRRI,                         &
                                  PTHM, PRM, ZLVOCPEXN, ZLSOCPEXN,       &
                                  ZTHLM, ZRTM                            )



! Virtual potential temperature at t-dt

ZTHVM(:,:,:) = PTHM(:,:,:)*((1.+XRV / XRD *PRM(:,:,:,1))/(1.+ZRTM(:,:,:))) 


       ZTHVREF=XG/ZTHVM
       CALL BL89(PZZ,PDZZ,ZTHVREF,ZTHLM,KRR, &
                 PRM,PTKEM,PLM)

!!! 2. Compute updraft


      CALL UPDRAFT_SOPE (KRR,KRRL,KRRI,OMIXUV,                     &
                         PZZ,PDZZ,PSFTH,PSFRV,PPABSM,PRHODREF,     &
                         PTKEM,PTHM,PRM,ZTHLM,ZRTM,ZUMM,ZVMM,      &
                         PTHL_UP,PRT_UP,PRV_UP,PU_UP,PV_UP,        &
                         PRC_UP,ZRI_UP,PTHV_UP,ZW_UP,PFRAC_UP,PEMF,&
                         PDETR,PENTR,IKLCL,IKETL,IKCTL )

!!! 3. Compute fluxes of conservative variables and their divergence = tendency

       ZEMF_O_RHODREF=PEMF/PRHODREF
       CALL MF_TURB(OMIXUV, PIMPL_MF, PTSTEP,PTSTEP_MET,       &
                PDZZ, PRHODJ, ZTHLM,ZTHVM,ZRTM,ZUMM,ZVMM,           &
                ZTHLDT,ZRTDT,ZUDT,ZVDT,                            &
                ZEMF_O_RHODREF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP,&
                ZFLXZTHMF,PFLXZTHVMF,ZFLXZRMF,ZFLXZUMF,ZFLXZVMF   )

               
!!! 5. Compute diagnostic convective cloud fraction and content
! ! !  ONLY liquid cloud implemented (yet)

      CALL COMPUTE_MF_CLOUD(KRRL,ZTHLM,PRC_UP,PFRAC_UP,PDZZ,IKLCL,  &
                            PRC_MF,ZCF_MF                      )


!!! 6. Compute tendency terms for pronostic variables

ZEXN(:,:,:)=(PPABSM(:,:,:)/XP00) ** (XRD/XCPD)
!


PRV(:,:,:)=PRM(:,:,:,1)-PRC_MF(:,:,:)
PRL(:,:,:)=PRC_MF(:,:,:)
      !       2.1 Cph
ZCPH(:,:,:)=XCPD+ XCPV * PRV(:,:,:)+ XCL * PRL(:,:,:) 
      
PTH(:,:,:)=(ZTHLM(:,:,:)+ZTHLDT(:,:,:))+(XLVTT/(ZCPH*ZEXN(:,:,:))*PRL(:,:,:))          

PRTHS(:,:,:)  = ZTHLDT(:,:,:)
PRTHS(:,:,:)  = (PTH(:,:,:)-PTHM(:,:,:))/PTSTEP_MET
PRRS(:,:,:,2) = (PRC_MF-PRM(:,:,:,2))/PTSTEP_MET
PRRS(:,:,:,1) = ZRTDT(:,:,:)-PRRS(:,:,:,2)

PRTHS(:,:,:)  = ZTHLDT(:,:,:)
PRRS(:,:,:,1) = ZRTDT(:,:,:)
PRRS(:,:,:,2) = 0
PRUS(:,:,:)   = ZUDT(:,:,:)
PRVS(:,:,:)   = ZVDT(:,:,:) 
  

END SUBROUTINE MFSHCONVPBL_CORE 


!     ###################################################################
      SUBROUTINE COMPUTE_BL89_ML(PDZZ2D, &
             PTKEM2D,PG_O_THVREF2D,PVPT,KK,OUPORDN,PLWORK)
!     ###################################################################
!!
!!     COMPUTE_BL89_ML routine to:
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!
!              ------------
!USE MODD_CST
!
!
IMPLICIT NONE
!
!          0.1 arguments
!
REAL, DIMENSION(:,:),   INTENT(IN)  :: PDZZ2D
REAL, DIMENSION(:,:),   INTENT(IN)  :: PTKEM2D
REAL, DIMENSION(:,:),   INTENT(IN)  :: PG_O_THVREF2D
REAL, DIMENSION(:,:),   INTENT(IN)  :: PVPT
INTEGER,                INTENT(IN)  :: KK
LOGICAL,                INTENT(IN)  :: OUPORDN
REAL, DIMENSION(:),     INTENT(OUT) :: PLWORK

!          0.2 Local variable
!
REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZLWORK1,ZLWORK2 ! Temporary mixing length
REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZINTE,ZPOTE     ! TKE and potential energy
                                                    ! between 2 levels
INTEGER :: IKB,IKE
!
REAL, DIMENSION(SIZE(PTKEM2D,1),SIZE(PTKEM2D,2)) :: ZDELTVPT,ZHLVPT                                
                      !Virtual Potential Temp at Half level and DeltaThv between
                      !2 levels
REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZTH! Potential Temp 

INTEGER :: IIJU,IKU             !Internal Domain
INTEGER :: J1D                  !horizontal loop counter
INTEGER :: JKK                  !loop counters
INTEGER :: JRR                  !moist loop counter
INTEGER :: JIJK                 !loop counters 
REAL    :: ZTEST,ZTEST0,ZTESTM  !test for vectorization
!-------------------------------------------------------------------------------------
!
!*       1.    INITIALISATION
!              --------------
IIJU=SIZE(PTKEM2D,1)
!
IKB = 1 !Modif WRF JP
IKE = SIZE(PTKEM2D,2)-1 !Modif WRF JP
IKU = SIZE(PTKEM2D,2)

ZDELTVPT(:,2:IKU)=PVPT(:,2:IKU)-PVPT(:,1:IKU-1)
ZDELTVPT(:,1)=0.

! to prevent division by zero
WHERE (ABS(ZDELTVPT(:,:))<1.E-10)
  ZDELTVPT(:,:)=1.E-10
END WHERE
!
ZHLVPT(:,2:IKU)= 0.5 * ( PVPT(:,2:IKU)+PVPT(:,1:IKU-1) )
ZHLVPT(:,1)    = PVPT(:,1)

!
!
!
!*       2.    CALCULATION OF THE UPWARD MIXING LENGTH
!              ---------------------------------------
!

IF (OUPORDN.EQV..TRUE.) THEN 
 ZINTE(:)=PTKEM2D(:,KK)
 PLWORK=0.
 ZLWORK1=0.
 ZLWORK2=0.
 ZTESTM=1.
 ZTH(:)=PVPT(:,KK)
 DO JKK=KK+1,IKE
    IF(ZTESTM > 0.) THEN
      ZTESTM=0
      DO J1D=1,IIJU
        ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
        ZPOTE(J1D) = ZTEST0*(PG_O_THVREF2D(J1D,KK)      *      &
            (ZHLVPT(J1D,JKK) - ZTH(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_THVREF2D(J1D,KK) *                     &
            (  PVPT(J1D,JKK-1) - ZTH(J1D) )                              &
          + SQRT (ABS(                                                       &
            ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK-1) - ZTH(J1D)) )**2  &
            + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK)                        &
                 * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) ))    ) /             &
        ( PG_O_THVREF2D(J1D,KK) * 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 
 ZINTE(:)=PTKEM2D(:,KK)
 PLWORK=0.
 ZLWORK1=0.
 ZLWORK2=0.
 ZTESTM=1.
 ZTH(:)=PVPT(:,KK)
 DO JKK=KK,IKB,-1 
    IF(ZTESTM > 0.) THEN
      ZTESTM=0
      DO J1D=1,IIJU
        ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
        ZPOTE(J1D) = -ZTEST0*(PG_O_THVREF2D(J1D,KK)      *      &
            (ZHLVPT(J1D,JKK) - ZTH(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(J1D)=        ( + PG_O_THVREF2D(J1D,KK) *                     &
            (  PVPT(J1D,JKK) - ZTH(J1D) )                              &
          + SQRT (ABS(                                                       &
            ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK) - ZTH(J1D)) )**2  &
            + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK)                        &
                 * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) ))    ) /             &
        ( PG_O_THVREF2D(J1D,KK) * 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
!

!
!         #############################################################
          SUBROUTINE COMPUTE_ENTR_DETR(OTEST,OTESTLCL,&
                            HFRAC_ICE,PFRAC_ICE,KK,PPABSM,PZZ,PDZZ,&
                            PTHVM,PTHLM,PRTM,PW_UP2,&
                            PTHL_UP,PRT_UP,PLUP,&
                            PENTR,PDETR,PBUO_INTEG)
!         #############################################################

!!
!!***COMPUTE_ENTR_DETR* - calculates caracteristics of the updraft or downdraft
!!                       using model of the EDMF scheme 
!!
!!
!! --------------------------------------------------------------------------
!



IMPLICIT NONE
!
!                         
!*                    1.1  Declaration of Arguments
!
!

LOGICAL,DIMENSION(:),INTENT(INOUT)  :: OTEST ! test to see if updraft is running
LOGICAL,DIMENSION(:),INTENT(INOUT)  :: OTESTLCL !test of condensation 
CHARACTER*1                         :: HFRAC_ICE ! frac_ice can be compute using
                                              ! Temperature (T) or prescribed
                                              ! (Y)
REAL, DIMENSION(:), INTENT(INOUT)   :: PFRAC_ICE ! if frac_ice is prescribed
INTEGER,            INTENT(IN)      :: KK  ! level where E and D are computed     
!
!    prognostic variables at t- deltat
!
REAL, DIMENSION(:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
REAL, DIMENSION(:,:),   INTENT(IN) ::  PZZ       !  Height at the flux point
REAL, DIMENSION(:,:),   INTENT(IN) ::  PDZZ       !  metrics coefficient
REAL, DIMENSION(:,:),   INTENT(IN) ::  PTHVM      ! ThetaV environment 

!
!   thermodynamical variables which are transformed in conservative var.
!
REAL, DIMENSION(:),   INTENT(IN)     ::  PTHLM     ! Thetal
REAL, DIMENSION(:),   INTENT(IN)     ::  PRTM      ! total mixing ratio 
REAL, DIMENSION(:,:), INTENT(INOUT)  ::  PW_UP2    ! Vertical velocity^2
REAL, DIMENSION(:),   INTENT(IN)     ::  PTHL_UP,PRT_UP  ! updraft properties
REAL, DIMENSION(:),   INTENT(IN)     ::  PLUP      ! LUP compute from the ground
REAL, DIMENSION(:),   INTENT(INOUT)  ::  PENTR     ! Mass flux entrainment of the updraft
REAL, DIMENSION(:),   INTENT(INOUT)  ::  PDETR     ! Mass flux detrainment of the updraft
REAL, DIMENSION(:),   INTENT(INOUT)  ::  PBUO_INTEG! Integral Buoyancy
!
!
!                       1.2  Declaration of local variables
!
!

! Variables for cloudy part

REAL, DIMENSION(SIZE(PTHLM))   :: ZKIC           ! fraction of env. mass in the muxtures
REAL, DIMENSION(SIZE(PTHLM))   :: ZEPSI,ZDELTA   ! factor entrainment detrainment
REAL, DIMENSION(SIZE(PTHLM))   :: ZEPSI_CLOUD    ! factor entrainment detrainment
REAL, DIMENSION(SIZE(PTHLM))   :: ZCOEFFMF_CLOUD ! factor for compputing entr. detr.

REAL, DIMENSION(SIZE(PTHLM))   :: ZMIXTHL,ZMIXRT ! Thetal and rt in the mixtures
!
REAL, DIMENSION(SIZE(PTHLM))   :: ZTHMIX,ZTHVMIX ! Theta and Thetav  of mixtures
REAL, DIMENSION(SIZE(PTHLM))   :: ZRVMIX,ZRCMIX,ZRIMIX ! mixing ratios in mixtures

REAL, DIMENSION(SIZE(PTHLM))   :: ZTHMIX_F2     ! Theta and Thetav  of mixtures
REAL, DIMENSION(SIZE(PTHLM))   :: ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2 ! mixing ratios in mixtures

REAL, DIMENSION(SIZE(PTHLM))   :: ZTHV_UP       ! thvup at mass point kk


REAL, DIMENSION(SIZE(PTHLM))   :: ZTHVMIX_1,ZTHVMIX_2 ! Theta and Thetav  of mixtures


! Variables for dry part

REAL, DIMENSION(SIZE(PTHLM))   :: ZBUO_INTEG,&         ! Temporary integral Buoyancy
                                  ZDZ_HALF,&           ! half-DeltaZ between 2 flux points
                                  ZDZ_STOP,&           ! Exact Height of the LCL 
                                  ZTHV_MINUS_HALF,&    ! Thv at flux point(kk)  
                                  ZTHV_PLUS_HALF,&     ! Thv at flux point(kk+1)
                                  ZCOEFF_MINUS_HALF,&  ! Variation of Thv between mass points kk-1 and kk
                                  ZCOEFF_PLUS_HALF,&   ! Variation of Thv between mass points kk and kk+1
                                  ZCOTHVU_MINUS_HALF,& ! Variation of Thvup between flux point kk and mass point kk
                                  ZCOTHVU_PLUS_HALF,&  ! Variation of Thvup between mass point kk and flux point kk+1    
                                  ZW2_HALF             ! w**2 at mass point KK  

REAL, DIMENSION(SIZE(PTHLM))   :: ZCOPRE_MINUS_HALF,&  ! Variation of pressure between mass points kk-1 and kk
                                  ZCOPRE_PLUS_HALF,&   ! Variation of pressure between mass points kk and kk+1
                                  ZPRE_MINUS_HALF,&    ! pressure at flux point kk
                                  ZPRE_PLUS_HALF,&     ! pressure at flux point kk+1
                                  ZTHV_UP_F1,&         ! thv_up at flux point kk
                                  ZTHV_UP_F2           ! thv_up at flux point kk+1
REAL, DIMENSION(SIZE(PTHLM))   :: ZCOEFF_QSAT,&        ! variation of Qsat at the transition between dry part and cloudy part
                                  ZRC_ORD,&            ! 
                                  ZPART_DRY            ! part of dry part at the transition level
!
REAL, DIMENSION(SIZE(PTHLM))   :: ZQVSAT ! QV at saturation

REAL, DIMENSION(SIZE(PTHLM))   :: PT ! temperature to compute saturation vapour mixing ratio 

REAL, DIMENSION(SIZE(PTHVM,1),SIZE(PTHVM,2))   ::ZG_O_THVREF
!
LOGICAL, DIMENSION(SIZE(OTEST,1)) :: GTEST_LOCAL_LCL,& ! true if LCL found between flux point KK and mass point KK  
                                     GTEST_LOCAL_LCL2  ! true if LCL found between mass point KK and flux point KK+1
!
REAL     :: ZRDORV       ! RD/RV
REAL     :: ZRVORD       ! RV/RD
INTEGER  :: ILON, ITEST, IKB         !


!----------------------------------------------------------------------------------
                        
!                1.3 Initialisation
!                ------------------
  IKB=1 ! modif WRF JP    

  
  ZRDORV   = XRD / XRV   !=0.622
  ZRVORD   = XRV / XRD   !=1.607
  ZG_O_THVREF=XG/PTHVM
  
  ZCOEFF_QSAT=0.
  ZRC_ORD=0.
  ZPART_DRY=1.
  GTEST_LOCAL_LCL=.FALSE.
  ZDZ_HALF(:) = (PZZ(:,KK+1)-PZZ(:,KK))/2.
  ZDZ_STOP(:) = ZDZ_HALF(:)
 
  ZKIC(:)=0.1  ! starting value for critical mixed fraction for CLoudy Part


!                Computation of KIC
!                ---------------------

!        2.1    Compute critical mixed fraction by estimating unknown  
!               T^mix r_c^mix and r_i^mix from thl^mix and r_t^mix
!               We determine the zero crossing of the linear curve
!               evaluating the derivative using ZMIXF=0.1.
!               -----------------------------------------------------

  ZMIXTHL(:) = ZKIC(:) * PTHLM(:)+(1. - ZKIC(:))*PTHL_UP(:)
  ZMIXRT(:)  = ZKIC(:) * PRTM(:)+(1. - ZKIC(:))*PRT_UP(:)

! MIXTURE FOR CLOUDY PART
  !  Compute pressure at flux level KK and at flux Level KK+1  


IF (KK==IKB) THEN !MODIF WRF JP
  ZCOPRE_MINUS_HALF(:) = 0. !MODIF WRF JP
ELSE!MODIF WRF JP
  ZCOPRE_MINUS_HALF(:) = ((PPABSM(:,KK)-PPABSM(:,KK-1))/PDZZ(:,KK))
ENDIF!MODIF WRF JP   
  ZCOPRE_PLUS_HALF(:) = ((PPABSM(:,KK+1)-PPABSM(:,KK))/PDZZ(:,KK+1))
  
IF (KK==IKB) THEN !MODIF WRF JP   
  ZPRE_MINUS_HALF(:)= PPABSM(:,KK)
ELSE!MODIF WRF JP
  ZPRE_MINUS_HALF(:)= ZCOPRE_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-1))+PPABSM(:,KK-1)
ENDIF!MODIF WRF JP   
  ZPRE_PLUS_HALF(:) = ZCOPRE_PLUS_HALF*0.5*(PZZ(:,KK+1)-PZZ(:,KK))+PPABSM(:,KK)  


  !  Compute non cons. var. of mixture at the mass level
  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
               PPABSM(:,KK),ZMIXTHL,ZMIXRT,&
               ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
             

  ! Compute theta_v of mixture at mass level KK for KF90        
  ZTHVMIX_1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))

  !  Compute non cons. var. of mixture at the flux level KK+1
  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
               ZPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
               ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
             

  ! compute theta_v of mixture at the flux level KK+1 for KF90       
  ZTHVMIX_2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))


!        2.1    Compute critical mixed fraction by estimating unknown  
!               T^mix r_c^mix and r_i^mix from thl^mix and r_t^mix
!               We determine the zero crossing of the linear curve
!               evaluating the derivative using ZMIXF=0.1.
!               -----------------------------------------------------


! THV_UP FOR DRY PART
  ! Compute theta_v of updraft at flux level KK                 
  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
               ZPRE_MINUS_HALF,PTHL_UP,PRT_UP,&
               ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
  ZTHV_UP_F1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
 
  ! Compute theta_v of updraft at mass level KK                 
  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
               PPABSM(:,KK),PTHL_UP,PRT_UP,&
               ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)            
  ZTHV_UP(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))

  ! Compute theta_v of updraft at flux level KK+1                   
  CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
               ZPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
               ZTHMIX_F2,ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2)
  ZTHV_UP_F2(:) = ZTHMIX_F2(:)*(1.+ZRVORD*ZRVMIX_F2(:))/(1.+PRT_UP(:))
 
   
!
!*   2.2     Compute final values for entr. and detr. 
!            ----------------------------------------
!
! Dry PART  

 ! Computation of integral entrainment and detrainment between flux level KK
 ! and mass level KK 

  WHERE ((ZRCMIX(:)>0.).AND.(.NOT.OTESTLCL)) 
! If rc is found between flux level KK and mass level KK
! a part of dry entrainment/detrainment is defined
! the exact height of LCL is also determined
     ZCOEFF_QSAT(:) = (ZRCMIX_F2(:) - ZRCMIX(:))/ ZDZ_HALF(:)
     ZRC_ORD(:) = ZRCMIX(:) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
     ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))     
     ZPART_DRY(:) = ZDZ_STOP / (PZZ(:,KK+1)-PZZ(:,KK))
     GTEST_LOCAL_LCL(:)=.TRUE.
  
  ENDWHERE

IF (KK==IKB) THEN !MODIF WRF JP   
  ZCOEFF_MINUS_HALF = 0.
ELSE!MODIF WRF JP
  ZCOEFF_MINUS_HALF = ((PTHVM(:,KK)-PTHVM(:,KK-1))/PDZZ(:,KK))
ENDIF!MODIF WRF JP   

  ZCOEFF_PLUS_HALF  = ((PTHVM(:,KK+1)-PTHVM(:,KK))/PDZZ(:,KK+1))
  
  ZCOTHVU_MINUS_HALF = (ZTHV_UP(:)-ZTHV_UP_F1(:))/ZDZ_HALF(:)
  ZCOTHVU_PLUS_HALF  = (ZTHV_UP_F2(:)-ZTHV_UP(:))/ZDZ_HALF(:)

IF (KK==IKB) THEN !MODIF WRF JP   
  ZTHV_MINUS_HALF = PTHVM(:,KK)
ELSE!MODIF WRF JP
  ZTHV_MINUS_HALF = ZCOEFF_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-1))+PTHVM(:,KK-1)
ENDIF!MODIF WRF JP   

  ZTHV_PLUS_HALF  = ZCOEFF_PLUS_HALF*0.5*(PZZ(:,KK+1)-PZZ(:,KK))+ ZTHV_MINUS_HALF  ! BUG JP 16082010
 
  ! Integral Buoyancy between flux level KK and mass level KK 
          
  PBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
              (0.5*( ZCOTHVU_MINUS_HALF - ZCOEFF_MINUS_HALF)*ZDZ_HALF(:) &
                - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) ) 


 
  WHERE ((OTEST).AND.(.NOT.OTESTLCL))
     PENTR=0.
     PDETR=0.

 
     ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
                (0.5 * (  - ZCOEFF_MINUS_HALF)* ZDZ_STOP(:) &
                  - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) ) 
           
  
     WHERE (ZBUO_INTEG(:)>=0.)
         PENTR = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
                 LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(:,KK))* &
                 ZBUO_INTEG)
         PDETR = 0.
    
         ZW2_HALF = PW_UP2(:,KK) +  2*(XABUO-XBENTR*XENTR_DRY)*(ZBUO_INTEG)
     ELSEWHERE
         PENTR = 0.
         PDETR = 0.5/(XABUO)*&
                 LOG(1.+ (2.*(XABUO)/PW_UP2(:,KK))* &
                 MAX(0.,-ZBUO_INTEG))

         ZW2_HALF = PW_UP2(:,KK) +  2*(XABUO)*(ZBUO_INTEG)
     ENDWHERE
     
 ENDWHERE

 
 ZDZ_STOP(:) = ZDZ_HALF(:)
  
! total Integral Buoyancy between flux level KK and flux level KK+1 

 PBUO_INTEG = PBUO_INTEG + ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
                (0.5*(ZCOTHVU_PLUS_HALF - ZCOEFF_PLUS_HALF)* ZDZ_HALF(:) - & 
                PTHVM(:,KK) + ZTHV_UP(:) ) 
              
               
 WHERE ((((ZRCMIX_F2(:)>0.).AND.(ZRCMIX(:)<=0.)).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
 ! If rc is found between mass level KK and flux level KK+1
 ! a part of dry entrainment is defined
 ! the exact height of LCL is also determined
     PT(:)=ZTHMIX_F2(:)*((PPABSM(:,KK+1)/XP00)**(XRD/XCPD)) 
     ZQVSAT(:)=EXP( XALPW - XBETAW/PT(:) - XGAMW*LOG(PT(:))  )
     ZQVSAT(:)=XRD/XRV*ZQVSAT(:)/PPABSM(:,KK+1)   &
                   / (1.+(XRD/XRV-1.)*ZQVSAT(:)/PPABSM(:,KK+1))   
     ZCOEFF_QSAT(:) = (PRT_UP(:) - ZQVSAT(:) - &
                ZRCMIX(:))/ (0.5* (PZZ(:,KK+2)-PZZ(:,KK+1)))
     ZRC_ORD(:) = ZRCMIX_F2(:) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
     ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))     
     ZPART_DRY(:) = 0.5+ZDZ_STOP  / (PZZ(:,KK+1)-PZZ(:,KK))
     GTEST_LOCAL_LCL2(:)=.TRUE.     
 ENDWHERE


 WHERE (((OTEST).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
   
     ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
                (0.5*( - ZCOEFF_PLUS_HALF)* ZDZ_STOP(:)&
                - PTHVM(:,KK) + ZTHV_UP(:) )

             

     WHERE (ZW2_HALF>0.)
        WHERE (ZBUO_INTEG(:)>=0.)
           PENTR = PENTR + 0.5/(XABUO-XBENTR*XENTR_DRY)* &
                LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/ZW2_HALF(:)) * ZBUO_INTEG)
          
           PDETR = PDETR
              
        ELSEWHERE
          PENTR = PENTR
          
          PDETR = PDETR + 0.5/(XABUO)* &
                LOG(1.+ (2.*(XABUO)/ZW2_HALF(:)) * &
                MAX(-ZBUO_INTEG,0.))
        ENDWHERE     
     ELSEWHERE
        ! if w**2<0 the updraft is stopped 
           OTEST=.FALSE.
           PENTR = PENTR 
           PDETR = PDETR 
     ENDWHERE
              
 ENDWHERE
 PENTR = XENTR_DRY*PENTR/(PZZ(:,KK+1)-PZZ(:,KK))    
 PDETR = XDETR_DRY*PDETR/(PZZ(:,KK+1)-PZZ(:,KK))
 PDETR = MAX(ZPART_DRY(:)*XDETR_LUP/(PLUP-0.5*(PZZ(:,KK)+PZZ(:,KK+1))),PDETR)


  

! compute final value of critical mixed fraction using theta_v
! of mixture, grid-scale and updraft in cloud

 WHERE ((OTEST).AND.(OTESTLCL))
     ZKIC(:) = MAX(0.,ZTHV_UP(:)-PTHVM(:,KK))*ZKIC(:) /  &  
                 (ZTHV_UP(:)-ZTHVMIX_1(:)+1.E-10)
                       
     ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
    
     ZEPSI(:) = ZKIC(:) **2.
     ZDELTA(:) = (1.-ZKIC(:))**2.
     ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
     ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF         
     PENTR(:) = ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
     PDETR(:) = ZCOEFFMF_CLOUD(:)*ZDELTA(:)
 ENDWHERE
 
! compute final value of critical mixed fraction using theta_v
! of mixture, grid-scale and updraft in cloud


 WHERE (((OTEST).AND.(.NOT.(OTESTLCL))).AND.((GTEST_LOCAL_LCL(:).OR.GTEST_LOCAL_LCL2(:))))
     ZKIC(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF)*ZKIC(:) /  &  
                       (ZTHV_UP_F2(:)-ZTHVMIX_2(:)+1.E-10)                      
     ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
     ZEPSI(:) = ZKIC(:) **2.
     ZDELTA(:) = (1.-ZKIC(:))**2.
     ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
     ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF     
     PENTR(:) = PENTR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
     PDETR(:) = PDETR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZDELTA(:)
 ENDWHERE
 


END SUBROUTINE COMPUTE_ENTR_DETR  

                                       
!     #################################################################
      SUBROUTINE COMPUTE_FUNCTION_THERMO_MF( KRR,KRRL,KRRI,                  &
                                       PTH, PR, PEXN, PPABS,                 &
                                       PT,PLVOCPEXN,PLSOCPEXN,               &
                                       PAMOIST,PATHETA                       )
!     #################################################################
!
!!
!!****  *COMPUTE_FUNCTION_THERMO_MF* - 
!!
!!
!! --------------------------------------------------------------------------
!       

!
IMPLICIT NONE
!
!
!*      0.1  declarations of arguments
!
INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.

REAL, DIMENSION(:,:,:), INTENT(IN)   :: PTH      ! theta
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR       ! water species
REAL, DIMENSION(:,:,:)  , INTENT(IN) :: PPABS,PEXN    ! pressure, Exner funct.

REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PT      ! temperature
REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLVOCPEXN,PLSOCPEXN     ! L/(cp*Pi)

REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PAMOIST,PATHETA
!
!-------------------------------------------------------------------------------
! 
!*       0.2   Declarations of local variables
!
REAL                :: ZEPS         ! XMV / XMD
REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)) ::     &
          ZCP,                        &  ! Cp 
          ZE,                         &  ! Saturation mixing ratio
          ZDEDT,                      &  ! Saturation mixing ratio derivative
          ZFRAC_ICE,                  &  ! Ice fraction
          ZAMOIST_W,                  &  ! Coefficients for s = f (Thetal,Rnp)
          ZATHETA_W,                  &  !
          ZAMOIST_I,                  &  !
          ZATHETA_I                      !

INTEGER             :: JRR
!
!-------------------------------------------------------------------------------
!
  ZEPS = XMV / XMD

  PLVOCPEXN(:,:,:) = 0.
  PLSOCPEXN(:,:,:) = 0.
  PAMOIST(:,:,:) = 0.
  PATHETA(:,:,:) = 0.
  ZFRAC_ICE(:,:,:) = 0.0

!
!*       Cph
!
ZCP=XCPD

IF (KRR > 0) ZCP(:,:,:) = ZCP(:,:,:) + XCPV * PR(:,:,:,1)

DO JRR = 2,1+KRRL  ! loop on the liquid components  
   ZCP(:,:,:)  = ZCP(:,:,:) + XCL * PR(:,:,:,JRR)
END DO

DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components   
  ZCP(:,:,:)  = ZCP(:,:,:)  + XCI * PR(:,:,:,JRR)
END DO

!*      Temperature
!
PT(:,:,:) =  PTH(:,:,:) * PEXN(:,:,:)
!
!
!! Liquid water
!
IF ( KRRL >= 1 ) THEN 
!
!*       Lv/Cph 
!
  PLVOCPEXN(:,:,:) = (XLVTT + (XCPV-XCL) *  (PT(:,:,:)-XTT) ) / ZCP(:,:,:)
!
!*      Saturation vapor pressure with respect to water
!
  ZE(:,:,:) =  EXP( XALPW - XBETAW/PT(:,:,:) - XGAMW*ALOG( PT(:,:,:) ) )
!
!*      Saturation  mixing ratio with respect to water
!
  ZE(:,:,:) =  ZE(:,:,:) * ZEPS / ( PPABS(:,:,:) - ZE(:,:,:) )
!
!*      Compute the saturation mixing ratio derivative (rvs')
!
  ZDEDT(:,:,:) = ( XBETAW / PT(:,:,:)  - XGAMW ) / PT(:,:,:)   &
                 * ZE(:,:,:) * ( 1. + ZE(:,:,:) / ZEPS )
!
!*      Compute Amoist
!
  ZAMOIST_W(:,:,:)=  0.5 / ( 1.0 + ZDEDT(:,:,:) * PLVOCPEXN(:,:,:) )
!
!*      Compute Atheta
!
  ZATHETA_W(:,:,:)= ZAMOIST_W(:,:,:) * PEXN(:,:,:) *                       &
        ( ( ZE(:,:,:) - PR(:,:,:,1) ) * PLVOCPEXN(:,:,:) /                  &
          ( 1. + ZDEDT(:,:,:) * PLVOCPEXN(:,:,:) )           *              &
          (                                                                 &
           ZE(:,:,:) * (1. + ZE(:,:,:)/ZEPS)                                &
                        * ( -2.*XBETAW/PT(:,:,:) + XGAMW ) / PT(:,:,:)**2   &
          +ZDEDT(:,:,:) * (1. + 2. * ZE(:,:,:)/ZEPS)                        &
                        * ( XBETAW/PT(:,:,:) - XGAMW ) / PT(:,:,:)          &
          )                                                                 &
         - ZDEDT(:,:,:)                                                     &
        )

!
!! Solid water
!
  IF ( KRRI >= 1 ) THEN 

!
!*       Fraction of ice
!
    WHERE(PR(:,:,:,2)+PR(:,:,:,4)>0.0)
      ZFRAC_ICE(:,:,:) = PR(:,:,:,4) / ( PR(:,:,:,2)+PR(:,:,:,4) )
    END WHERE
!
!*       Ls/Cph 
!
    PLSOCPEXN(:,:,:) = (XLSTT + (XCPV-XCI) *  (PT(:,:,:)-XTT) ) / ZCP(:,:,:)
!
!*      Saturation vapor pressure with respect to ice
!
    ZE(:,:,:) =  EXP( XALPI - XBETAI/PT(:,:,:) - XGAMI*ALOG( PT(:,:,:) ) )
!
!*      Saturation  mixing ratio with respect to ice
!
    ZE(:,:,:) =  ZE(:,:,:) * ZEPS / ( PPABS(:,:,:) - ZE(:,:,:) )
!
!*      Compute the saturation mixing ratio derivative (rvs')
!
    ZDEDT(:,:,:) = ( XBETAI / PT(:,:,:)  - XGAMI ) / PT(:,:,:)   &
                   * ZE(:,:,:) * ( 1. + ZE(:,:,:) / ZEPS )
!
!*      Compute Amoist
!
    ZAMOIST_I(:,:,:)=  0.5 / ( 1.0 + ZDEDT(:,:,:) * PLSOCPEXN(:,:,:) )
!
!*      Compute Atheta
!
    ZATHETA_I(:,:,:)= ZAMOIST_I(:,:,:) * PEXN(:,:,:) *                     &
        ( ( ZE(:,:,:) - PR(:,:,:,1) ) * PLSOCPEXN(:,:,:) /                  &
          ( 1. + ZDEDT(:,:,:) * PLSOCPEXN(:,:,:) )           *              &
          (                                                                 &
           ZE(:,:,:) * (1. + ZE(:,:,:)/ZEPS)                                &
                        * ( -2.*XBETAI/PT(:,:,:) + XGAMI ) / PT(:,:,:)**2   &
          +ZDEDT(:,:,:) * (1. + 2. * ZE(:,:,:)/ZEPS)                        &
                        * ( XBETAI/PT(:,:,:) - XGAMI ) / PT(:,:,:)          &
          )                                                                 &
         - ZDEDT(:,:,:)                                                     &
        )


  ENDIF

  PAMOIST(:,:,:) = (1.0-ZFRAC_ICE(:,:,:))*ZAMOIST_W(:,:,:) &
                         +ZFRAC_ICE(:,:,:) *ZAMOIST_I(:,:,:)
  PATHETA(:,:,:) = (1.0-ZFRAC_ICE(:,:,:))*ZATHETA_W(:,:,:) &
                         +ZFRAC_ICE(:,:,:) *ZATHETA_I(:,:,:)

!
!*      Lv/Cph/Exner and  Ls/Cph/Exner
!
  PLVOCPEXN(:,:,:) = PLVOCPEXN(:,:,:) / PEXN(:,:,:)
  PLSOCPEXN(:,:,:) = PLSOCPEXN(:,:,:) / PEXN(:,:,:)
!
ENDIF

END SUBROUTINE COMPUTE_FUNCTION_THERMO_MF
!
!     #################################################################
      SUBROUTINE COMPUTE_UPDRAFT(OMIXUV,PZZ,PDZZ,KK,              &
                                 PSFTH,PSFRV,                     &
                                 PPABSM,PRHODREF,PUM,PVM, PTKEM,  &
                                 PTHM,PRVM,PRCM,PRIM,PTHLM,PRTM,  &
                                 PTHL_UP,PRT_UP,                  &
                                 PRV_UP,PRC_UP,PRI_UP,PTHV_UP,    &
                                 PW_UP,PU_UP, PV_UP,              &
                                 PFRAC_UP,PEMF,PDETR,PENTR,       &
                                 KKLCL,KKETL,KKCTL)
!     #################################################################
!!
!!****  *COMPUTE_UPDRAFT* - calculates caracteristics of the updraft 
!!                         
!! --------------------------------------------------------------------------



IMPLICIT NONE

!*                    1.1  Declaration of Arguments
!
!
!
LOGICAL,                INTENT(IN) :: OMIXUV    ! True if mixing of momentum
REAL, DIMENSION(:,:), INTENT(IN)   :: PZZ       !  Height at the flux point
REAL, DIMENSION(:,:), INTENT(IN)   :: PDZZ      !  Metrics coefficient
 
INTEGER,              INTENT(IN)   ::  KK        
REAL, DIMENSION(:),   INTENT(IN)   ::  PSFTH,PSFRV
! normal surface fluxes of theta,rv,(u,v) parallel to the orography
!
REAL, DIMENSION(:,:),   INTENT(IN) ::  PPABSM     ! Pressure at t-dt
REAL, DIMENSION(:,:),   INTENT(IN) ::  PRHODREF   ! dry density of the
                                                  ! reference state
REAL, DIMENSION(:,:),   INTENT(IN) ::  PUM        ! u mean wind
REAL, DIMENSION(:,:),   INTENT(IN) ::  PVM        ! v mean wind
REAL, DIMENSION(:,:),   INTENT(IN) ::  PTKEM      ! TKE at t-dt
!
REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTHM           ! liquid pot. temp. at t-dt
REAL, DIMENSION(:,:),   INTENT(IN)   ::  PRVM,PRCM,PRIM ! vapor mixing ratio at t-dt
REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTHLM,PRTM     ! cons. var. at t-dt


REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PTHL_UP,PRT_UP   ! updraft properties
REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PU_UP, PV_UP     ! updraft wind components
REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PRV_UP,PRC_UP, & ! updraft rv, rc
                                         PRI_UP,PTHV_UP,& ! updraft ri, THv
                                         PW_UP,PFRAC_UP   ! updraft w, fraction

                                         
REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PEMF,PDETR,PENTR ! Mass_flux,
                                                          ! detrainment,entrainment
INTEGER, DIMENSION(:),  INTENT(OUT)  ::  KKLCL,KKETL,KKCTL! LCL, ETL, CTL                                     
!                       1.2  Declaration of local variables
!
!
! Mean environment variables at t-dt at flux point
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::    &
                        ZTHM_F,ZRVM_F,ZRCM_F,ZRIM_F   ! Theta,rv,rl,ri of
                                                      ! updraft environnement
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::    &
                        ZRTM_F, ZTHLM_F, ZTKEM_F,&    ! rt, thetal,TKE,pressure,
                        ZUM_F,ZVM_F,ZRHO_F,      &    ! density,momentum
                        ZPRES_F,ZTHVM_F,ZTHVM,   &    ! interpolated at the flux point
                        ZG_O_THVREF,             &    ! g*ThetaV ref
                        ZW_UP2                        ! w**2  of the updraft

                        
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::  &
                        ZTH_UP,                  &    ! updraft THETA 
                        ZFRAC_ICE                     ! Ice fraction

REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::  ZCOEF  ! diminution coefficient for too high clouds 
                        
REAL, DIMENSION(SIZE(PSFTH,1) )            ::  ZWTHVSURF  ! Surface w'thetav'
CHARACTER(LEN=1)                           ::  YFRAC_ICE  ! Ice Fraction
                                                          ! precribed or computed 

REAL  :: ZRDORV       ! RD/RV
REAL  :: ZRVORD       ! RV/RD


REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3
REAL, DIMENSION(SIZE(PTHM,1)) :: ZBUO_INTEG           ! Integrated Buoyancy 

REAL, DIMENSION(SIZE(PTHM,1)) :: ZLUP         ! Upward Mixing length from the ground
REAL, DIMENSION(SIZE(PTHM,1)) :: ZDEPTH       ! Deepness limit for cloud

INTEGER  :: IKB,IKE            ! index value for the Beginning and the End
                               ! of the physical domain for the mass points
INTEGER  :: IKU           ! array size in k
INTEGER  :: JK,JI,JSV          ! loop counters

LOGICAL, DIMENSION(SIZE(PTHM,1)) ::  GTEST,GTESTLCL,GTESTETL
                               ! Test if the ascent continue, if LCL or ETL is reached
LOGICAL                          ::  GLMIX 
                               ! To choose upward or downward mixing length
LOGICAL, DIMENSION(SIZE(PTHM,1))              :: GWORK1
LOGICAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: GWORK2

INTEGER  :: ITEST
REAL  :: ZTMAX,ZRMAX  ! control value
! Thresholds for the  perturbation of
! theta_l and r_t at the first level of the updraft
ZTMAX=2.0
ZRMAX=1.E-3
!------------------------------------------------------------------------

!                     INITIALISATION

!                 Local variables, internal domain
! Internal Domain
IKU=SIZE(PTHM,2)
IKB=1        ! modif WRF JP
IKE=IKU-1      ! modif WRF JP

! Initialisation of intersesting Level :LCL,ETL,CTL
KKLCL(:)=IKE
KKETL(:)=IKE
KKCTL(:)=IKE

! 
! Initialisation
PRV_UP(:,:)=0.
PRC_UP(:,:)=0.

PW_UP(:,:)=0.
PEMF(:,:)=0.
PDETR(:,:)=0.
PENTR(:,:)=0.
ZTH_UP(:,:)=0.
PFRAC_UP(:,:)=0.
PTHV_UP(:,:)=0.

!no ice cloud coded yet 
PRI_UP(:,:)=0.
ZFRAC_ICE(:,:)=0.
YFRAC_ICE='T'

ZBUO_INTEG=0.
! Initialisation of the constants   
ZRDORV   = XRD / XRV   !=0.622
ZRVORD   = (XRV / XRD) 


! Initialisation of environment variables at t-dt

! variables at flux level


ZTHM_F (:,IKB+1:IKU) = 0.5*(PTHM(:,IKB:IKU-1)+PTHM(:,IKB+1:IKU))
ZTHLM_F(:,IKB+1:IKU) = 0.5*(PTHLM(:,IKB:IKU-1)+PTHLM(:,IKB+1:IKU))
ZRTM_F (:,IKB+1:IKU) = 0.5*(PRTM(:,IKB:IKU-1)+PRTM(:,IKB+1:IKU))
ZTKEM_F(:,IKB+1:IKU) = 0.5*(PTKEM(:,IKB:IKU-1)+PTKEM(:,IKB+1:IKU))
ZPRES_F(:,IKB+1:IKU) = 0.5*(PPABSM(:,IKB:IKU-1)+PPABSM(:,IKB+1:IKU))
ZRHO_F (:,IKB+1:IKU) = 0.5*(PRHODREF(:,IKB:IKU-1)+PRHODREF(:,IKB+1:IKU))
ZRVM_F (:,IKB+1:IKU) = 0.5*(PRVM(:,IKB:IKU-1)+PRVM(:,IKB+1:IKU))
ZUM_F  (:,IKB+1:IKU) = 0.5*(PUM(:,IKB:IKU-1)+PUM(:,IKB+1:IKU))
ZVM_F  (:,IKB+1:IKU) = 0.5*(PVM(:,IKB:IKU-1)+PVM(:,IKB+1:IKU))



ZTHM_F (:,IKB) = PTHM(:,IKB) 
ZTHLM_F(:,IKB) = PTHLM(:,IKB)
ZRTM_F (:,IKB) = PRTM (:,IKB)
ZTKEM_F(:,IKB) = PTKEM(:,IKB)
ZPRES_F(:,IKB) = PPABSM(:,IKB)
ZRHO_F(:,IKB)  = PRHODREF(:,IKB)
ZRVM_F(:,IKB)  = PRVM(:,IKB)
ZUM_F(:,IKB)   = PUM(:,IKB)
ZVM_F(:,IKB)   = PVM(:,IKB)


! thetav at mass and flux levels 
ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
ZTHVM(:,:)=PTHM(:,:)*((1.+ZRVORD*PRVM(:,:))/(1.+PRTM(:,:)))


!                     
!          Initialisation of updraft characteristics 
PTHL_UP(:,:)=ZTHLM_F(:,:)
PRT_UP(:,:)=ZRTM_F(:,:)
ZW_UP2(:,:)=0.
PTHV_UP(:,:)=ZTHVM_F(:,:)
PU_UP(:,:)=ZUM_F(:,:)
PV_UP(:,:)=ZVM_F(:,:)


! Computation or initialisation of updraft characteristics at the KK level
! thetal_up,rt_up,thetaV_up, w²,Buoyancy term and mass flux (PEMF)

! 03/2009
!PTHL_UP(:,KK)= ZTHLM_F(:,KK)+(PSFTH(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT
!PRT_UP(:,KK) = ZRTM_F(:,KK)+(PSFRV(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT
PTHL_UP(:,KK)= ZTHLM_F(:,KK)+MAX(0.,MIN(ZTMAX,(PSFTH(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT))
PRT_UP(:,KK) = ZRTM_F(:,KK)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT)) 



ZW_UP2(:,KK) = MAX(0.0001,(2./3.)*ZTKEM_F(:,KK))


! Computation of non conservative variable for the KK level of the updraft
! (all or nothing ajustement)
CALL TH_R_FROM_THL_RT_2D(YFRAC_ICE,ZFRAC_ICE(:,KK:KK),ZPRES_F(:,KK:KK), &
             PTHL_UP(:,KK:KK),PRT_UP(:,KK:KK),ZTH_UP(:,KK:KK), &
             PRV_UP(:,KK:KK),PRC_UP(:,KK:KK),PRI_UP(:,KK:KK))

! compute updraft thevav and buoyancy term at KK level             
PTHV_UP(:,KK) = ZTH_UP(:,KK)*((1+ZRVORD*PRV_UP(:,KK))/(1+PRT_UP(:,KK))) 
                                                            
! Closure assumption for mass flux at KK level                                                  
ZG_O_THVREF=XG/ZTHVM_F

! compute L_up
GLMIX=.TRUE.
ZTKEM_F(:,KK)=0.

CALL COMPUTE_BL89_ML(PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KK,GLMIX,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)

! Compute Buoyancy flux at the ground
ZWTHVSURF(:) = (ZTHVM_F(:,IKB)/ZTHM_F(:,IKB))*PSFTH(:)+     &
                (0.61*ZTHM_F(:,IKB))*PSFRV(:)

! Mass flux at KK level (updraft triggered if PSFTH>0.)
WHERE (ZWTHVSURF(:)>0.)
  PEMF(:,KK) = XCMF * ZRHO_F(:,KK) * ((ZG_O_THVREF(:,KK))*ZWTHVSURF*ZLUP)**(1./3.)
  PFRAC_UP(:,KK)=MIN(PEMF(:,KK)/(SQRT(ZW_UP2(:,KK))*ZRHO_F(:,KK)),XFRAC_UP_MAX)
  ZW_UP2(:,KK)=(PEMF(:,KK)/(PFRAC_UP(:,KK)*ZRHO_F(:,KK)))**2
  GTEST(:)=.TRUE.
ELSEWHERE
  PEMF(:,KK) =0.
  GTEST(:)=.FALSE.
ENDWHERE


!--------------------------------------------------------------------------

!                        3. Iteration
!                           ---------
!
! If GTEST = T the updraft starts from the KK level and stops when GTEST becomes F
!
!
GTESTLCL(:)=.FALSE.
GTESTETL(:)=.FALSE.

!       Loop on vertical level

DO JK=KK,IKE-1
! IF the updraft top is reached for all column, stop the loop on levels
ITEST=COUNT(GTEST)
IF (ITEST==0) CYCLE
!       Computation of entrainment and detrainment with KF90
!       parameterization in clouds and LR01 in subcloud layer


! to find the LCL (check if JK is LCL or not)

  WHERE(GTEST)
  WHERE ((PRC_UP(:,JK)>0.).AND.(.NOT.(GTESTLCL)))
      KKLCL(:) = JK           
      GTESTLCL(:)=.TRUE.
  ENDWHERE
  ENDWHERE


! COMPUTE PENTR and PDETR at mass level JK

 CALL COMPUTE_ENTR_DETR(GTEST,GTESTLCL,YFRAC_ICE,ZFRAC_ICE(:,JK),JK,&
                       PPABSM(:,:),PZZ(:,:),PDZZ(:,:),ZTHVM(:,:),&
                       PTHLM(:,JK),PRTM(:,JK),ZW_UP2(:,:),   &
                       PTHL_UP(:,JK),PRT_UP(:,JK),ZLUP(:), &
                       PENTR(:,JK),PDETR(:,JK),ZBUO_INTEG)
  

  IF (JK==KK) THEN
       PDETR(:,JK)=0.
  ENDIF   

 
!       Computation of updraft characteristics at level JK+1
  WHERE(GTEST)
    ZMIX1(:)=0.5*(PZZ(:,JK+1)-PZZ(:,JK))*(PENTR(:,JK)-PDETR(:,JK))
    PEMF(:,JK+1)=PEMF(:,JK)*EXP(2*ZMIX1(:))
  ENDWHERE
  
  

! stop the updraft if MF becomes negative
  WHERE (GTEST.AND.(PEMF(:,JK+1)<=0.))
    PEMF(:,JK+1)=0.
    GTEST(:)=.FALSE.
    KKCTL(:) = JK+1         
  ENDWHERE


! If the updraft did not stop, compute cons updraft characteritics at jk+1
  WHERE(GTEST)     
    ZMIX2(:) = (PZZ(:,JK+1)-PZZ(:,JK))*PENTR(:,JK) !&
    ZMIX3(:) = (PZZ(:,JK+1)-PZZ(:,JK))*PDETR(:,JK) !&                   
                
    PTHL_UP(:,JK+1) = (PTHL_UP(:,JK)*(1.-0.5*ZMIX2(:)) + PTHLM(:,JK)*ZMIX2(:)) &
                          /(1.+0.5*ZMIX2(:))   
    PRT_UP(:,JK+1) = (PRT_UP (:,JK)*(1.-0.5*ZMIX2(:)) + PRTM(:,JK)*ZMIX2(:))   &
                          /(1.+0.5*ZMIX2(:))
  ENDWHERE
  

  IF(OMIXUV) THEN
    WHERE(GTEST) 
      PU_UP(:,JK+1) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
                        0.5*XPRES_UV*(PZZ(:,JK+1)-PZZ(:,JK))*&
                        ((PUM(:,JK+1)-PUM(:,JK))/PDZZ(:,JK+1)+&
                         (PUM(:,JK)-PUM(:,JK-1))/PDZZ(:,JK))        )   &
                        /(1+0.5*ZMIX2(:))
      PV_UP(:,JK+1) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
                        0.5*XPRES_UV*(PZZ(:,JK+1)-PZZ(:,JK))*&
                        ((PVM(:,JK+1)-PVM(:,JK))/PDZZ(:,JK+1)+&
                         (PVM(:,JK)-PVM(:,JK-1))/PDZZ(:,JK))    )   &
                        /(1+0.5*ZMIX2(:))
    ENDWHERE
  ENDIF
  
! Compute non cons. var. at level JK+1
  CALL TH_R_FROM_THL_RT_2D(YFRAC_ICE,ZFRAC_ICE(:,JK+1:JK+1),ZPRES_F(:,JK+1:JK+1), &
          PTHL_UP(:,JK+1:JK+1),PRT_UP(:,JK+1:JK+1),ZTH_UP(:,JK+1:JK+1), &
          PRV_UP(:,JK+1:JK+1),PRC_UP(:,JK+1:JK+1),PRI_UP(:,JK+1:JK+1))
  

! Compute the updraft theta_v, buoyancy and w**2 for level JK+1   
  WHERE(GTEST)
      PTHV_UP(:,JK+1) = ZTH_UP(:,JK+1)*((1+ZRVORD*PRV_UP(:,JK+1))/(1+PRT_UP(:,JK+1)))
   WHERE (.NOT.(GTESTLCL)) 
      WHERE (ZBUO_INTEG(:)>0.)
        ZW_UP2(:,JK+1)  = ZW_UP2(:,JK) + 2.*(XABUO-XBENTR*XENTR_DRY)* ZBUO_INTEG(:)              
      ENDWHERE
      WHERE (ZBUO_INTEG(:)<=0.)
        ZW_UP2(:,JK+1)  = ZW_UP2(:,JK) + 2.*XABUO* ZBUO_INTEG(:)
      ENDWHERE      
   ENDWHERE      
   WHERE (GTESTLCL)
      ZW_UP2(:,JK+1)  = ZW_UP2(:,JK)*(1.-(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:)))&
            /(1.+(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:))) &
            +2.*(XABUO)*ZBUO_INTEG/(1.+(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:)))
   ENDWHERE
 ENDWHERE


! Test if the updraft has reach the ETL
  GTESTETL(:)=.FALSE.
  WHERE (GTEST.AND.(ZBUO_INTEG(:)<=0.))
      KKETL(:) = JK+1           
      GTESTETL(:)=.TRUE.
  ENDWHERE

! Test is we have reached the top of the updraft
  WHERE (GTEST.AND.((ZW_UP2(:,JK+1)<=0.).OR.(PEMF(:,JK+1)<=0.)))
      ZW_UP2(:,JK+1)=0.
      PEMF(:,JK+1)=0.
      GTEST(:)=.FALSE.
      PTHL_UP(:,JK+1)=ZTHLM_F(:,JK+1)
      PRT_UP(:,JK+1)=ZRTM_F(:,JK+1)
      PRC_UP(:,JK+1)=0.
      PRV_UP(:,JK+1)=0.
      PTHV_UP(:,JK+1)=ZTHVM_F(:,JK+1)
      PFRAC_UP(:,JK+1)=0.
      KKCTL(:)=JK+1
  ENDWHERE
 
! compute frac_up at JK+1
  WHERE (GTEST)
      PFRAC_UP(:,JK+1)=PEMF(:,JK+1)/(SQRT(ZW_UP2(:,JK+1))*ZRHO_F(:,JK+1))
  ENDWHERE

! Updraft fraction must be smaller than XFRAC_UP_MAX
  WHERE (GTEST)
      PFRAC_UP(:,JK+1)=MIN(XFRAC_UP_MAX,PFRAC_UP(:,JK+1))
  ENDWHERE


! When cloudy and non-buoyant, updraft fraction must decrease
  
  WHERE ((GTEST.AND.GTESTETL).AND.GTESTLCL)
    PFRAC_UP(:,JK+1)=MIN(PFRAC_UP(:,JK+1),PFRAC_UP(:,JK))
  ENDWHERE

! Mass flux is updated with the new updraft fraction
  
  PEMF(:,JK+1)=PFRAC_UP(:,JK+1)*SQRT(ZW_UP2(:,JK+1))*ZRHO_F(:,JK+1)

ENDDO


PW_UP(:,:)=SQRT(ZW_UP2(:,:))

PEMF(:,KK) =0.
PEMF(:,IKU)=0. ! modif WRF JP
PFRAC_UP(:,IKU)=0.


DO JI=1,SIZE(PTHM,1) 
   ZDEPTH(JI) = (PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
END DO

GWORK1(:)= (GTESTLCL(:) .AND. (ZDEPTH(:) > 3000.) )
GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=IKU )
ZCOEF(:,:) = SPREAD( (1.-(ZDEPTH(:)-3000.)/1000.), DIM=2, NCOPIES=IKU)
ZCOEF=MIN(MAX(ZCOEF,0.),1.)

WHERE (GWORK2) 
   PEMF(:,:)     = PEMF(:,:)     * ZCOEF(:,:)
   PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
ENDWHERE


END SUBROUTINE COMPUTE_UPDRAFT


                
!     #################################################################
      SUBROUTINE MF_TURB(OMIXUV, PIMPL, PTSTEP,                       &
                PTSTEP_MET,  PDZZ,                          &
                PRHODJ,                                               &
                PTHLM,PTHVM,PRTM,PUM,PVM,                             &
                PTHLDT,PRTDT,PUDT,PVDT,                               &
                PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP,              &
                PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF       )

!     #################################################################
!
!
!!****  *MF_TURB* - computes the MF_turbulent source terms for the prognostic
!!                  variables. 
!!

!! --------------------------------------------------------------------------

!
!
IMPLICIT NONE
!
!
!*      0.1  declarations of arguments
!
!
LOGICAL,                INTENT(IN)   :: OMIXUV      ! True if mixing of momentum
REAL,                   INTENT(IN)   :: PIMPL       ! degree of implicitness
REAL,                 INTENT(IN)     ::  PTSTEP   ! Dynamical timestep 
REAL,                 INTENT(IN)     ::  PTSTEP_MET! Timestep for meteorological variables                        
!
REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDZZ        ! metric coefficients

REAL, DIMENSION(:,:,:), INTENT(IN)   :: PRHODJ      ! dry density * Grid size

!   Conservative var. at t-dt
REAL, DIMENSION(:,:,:), INTENT(IN) ::  PTHLM        ! conservative pot. temp.
REAL, DIMENSION(:,:,:), INTENT(IN) ::  PRTM         ! water var.  where 
!  Virtual potential temperature at t-dt
REAL, DIMENSION(:,:,:), INTENT(IN) ::  PTHVM 
!  Momentum at t-dt
REAL, DIMENSION(:,:,:), INTENT(IN) ::  PUM
REAL, DIMENSION(:,:,:), INTENT(IN) ::  PVM
!
! Tendencies of conservative variables
REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PTHLDT

REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PRTDT 
! Tendencies of momentum
REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PUDT
REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PVDT
! Tendencies of scalar variables


! Updraft characteritics
REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP
! Fluxes
REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF

!
!
!
!-------------------------------------------------------------------------------
!
!       0.2  declaration of local variables
!

REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: ZVARS

!
!----------------------------------------------------------------------------
!
!*      1.PRELIMINARIES
!         -------------
!
!

!
PFLXZTHMF = 0.
PFLXZRMF = 0.
PFLXZTHVMF = 0.
PFLXZUMF = 0.
PFLXZVMF = 0.
PTHLDT = 0.
PRTDT = 0.
PUDT = 0.
PVDT = 0.

!
!----------------------------------------------------------------------------
!
!*      2. COMPUTE THE MEAN FLUX OF CONSERVATIVE VARIABLES at time t-dt
!          (equation (3) of Soares et al)
!          + THE MEAN FLUX OF THETA_V (buoyancy flux)
!          -----------------------------------------------
!   ( Resulting fluxes are in flux level (w-point) as PEMF and PTHL_UP )
!

PFLXZTHMF(:,:,:) = PEMF(:,:,:)*(PTHL_UP(:,:,:)-MZM(PTHLM(:,:,:)))

PFLXZRMF(:,:,:) =  PEMF(:,:,:)*(PRT_UP(:,:,:)-MZM(PRTM(:,:,:)))

PFLXZTHVMF(:,:,:) = PEMF(:,:,:)*(PTHV_UP(:,:,:)-MZM(PTHVM(:,:,:)))

PFLXZTHVMF(:,:,:) = 9.81/PTHVM(:,:,:)* PEMF(:,:,:)*(PTHV_UP(:,:,:)-MZM(PTHVM(:,:,:))) !JP

IF (OMIXUV) THEN
  PFLXZUMF(:,:,:) =  PEMF(:,:,:)*(PU_UP(:,:,:)-MZM(PUM(:,:,:)))
  PFLXZVMF(:,:,:) =  PEMF(:,:,:)*(PV_UP(:,:,:)-MZM(PVM(:,:,:)))
ENDIF
!
!
!----------------------------------------------------------------------------
!
!*      3. COMPUTE TENDENCIES OF CONSERVATIVE VARIABLES (or treated as such...)
!          (implicit formulation)
!          --------------------------------------------
!

!
!
! 3.1 Compute the tendency for the conservative potential temperature
!     (PDZZ and flux in w-point and PRHODJ is mass point, result in mass point)
!
CALL TRIDIAG_MASSFLUX(PTHLM,PFLXZTHMF,-PEMF,PTSTEP_MET,PIMPL,  &
                      PDZZ,PRHODJ,ZVARS )
! compute new flux
PFLXZTHMF(:,:,:) = PEMF(:,:,:)*(PTHL_UP(:,:,:)-MZM(ZVARS(:,:,:)))

!!! compute THL tendency
!
PTHLDT(:,:,:)= (ZVARS(:,:,:)-PTHLM(:,:,:))/PTSTEP_MET

!
! 3.2 Compute the tendency for the conservative mixing ratio
!
CALL TRIDIAG_MASSFLUX(PRTM(:,:,:),PFLXZRMF,-PEMF,PTSTEP_MET,PIMPL,  &
                                 PDZZ,PRHODJ,ZVARS )
! compute new flux
PFLXZRMF(:,:,:) =  PEMF(:,:,:)*(PRT_UP(:,:,:)-MZM(ZVARS(:,:,:)))

!!! compute RT tendency
PRTDT(:,:,:) = (ZVARS(:,:,:)-PRTM(:,:,:))/PTSTEP_MET
!
! 3.3 Compute the tendency for the (non conservative but treated as it)  mixing ratio
!
!CALL TRIDIAG_MASSFLUX(PTHVM(:,:,:),PFLXZTHVMF,-PEMF,PTSTEP,PIMPL,  &
!                                PDZZ,PRHODJ,ZVARS )
! compute new flux
!PFLXZTHVMF(:,:,:) =  PEMF(:,:,:)*(PTHV_UP(:,:,:)-MZM(ZVARS(:,:,:)))



IF (OMIXUV) THEN
!
! 3.3 Compute the tendency for the (non conservative but treated as it) zonal momentum
!     (PDZZ and flux in w-point and PRHODJ is mass point, result in mass point)
!

 CALL TRIDIAG_MASSFLUX(PUM,PFLXZUMF,-PEMF,PTSTEP,PIMPL,  &
                                 PDZZ,PRHODJ,ZVARS )
! compute new flux
PFLXZUMF(:,:,:) = PEMF(:,:,:)*(PU_UP(:,:,:)-MZM(ZVARS(:,:,:)))

! compute U tendency
  PUDT(:,:,:)= (ZVARS(:,:,:)-PUM(:,:,:))/PTSTEP_MET

!
!
! 3.4 Compute the tendency for the (non conservative but treated as it for the time beiing) 
!                                  meridian momentum
!     (PDZZ and flux in w-point and PRHODJ is mass point, result in mass point)
!
  CALL TRIDIAG_MASSFLUX(PVM,PFLXZVMF,-PEMF,PTSTEP,PIMPL,  &
                                 PDZZ,PRHODJ,ZVARS )
! compute new flux
PFLXZVMF(:,:,:) = PEMF(:,:,:)*(PV_UP(:,:,:)-MZM(ZVARS(:,:,:)))

! compute V tendency
  PVDT(:,:,:)= (ZVARS(:,:,:)-PVM(:,:,:))/PTSTEP_MET

ENDIF



!
END SUBROUTINE MF_TURB   

FUNCTION MZM(PA)  RESULT(PMZM)


!
IMPLICIT NONE
!
!*       0.1   Declarations of argument and result
!              ------------------------------------
!

REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass localization
REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM   ! result at flux localization

!        0.2   Declarations of local variables
!              -------------------------------
!

INTEGER :: JK             ! Loop index in z direction
INTEGER :: IKU            ! upper bound in z direction of PA
!
IKU = SIZE(PA,3)

!DO JK=2,IKU MODIF WRF JP
DO JK=2,IKU
   PMZM(:,:,JK)=0.5*(PA(:,:,JK)+PA(:,:,JK-1))
ENDDO

PMZM(:,:,1)=PA(:,:,2)


END FUNCTION MZM


!
!
!     #################################################################
      SUBROUTINE TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,PP,             &
                                  PTHL, PRT, PTH, PRV, PRL, PRI       )
!     #################################################################
!
!
!!****  *TH_R_FROM_THL_RT_1D* - computes the non-conservative variables
!!                          from conservative variables
!!

!!
!! --------------------------------------------------------------------------
!       

!
IMPLICIT NONE
!
!
!*      0.1  declarations of arguments
!
CHARACTER*1           , INTENT(IN) :: HFRAC_ICE
REAL, DIMENSION(:), INTENT(INOUT) :: PFRAC_ICE
REAL, DIMENSION(:), INTENT(IN) :: PP          ! Pressure
REAL, DIMENSION(:), INTENT(IN) :: PTHL    ! thetal (or thetal tendency) to
                                             ! transform into th (or tendency)
REAL, DIMENSION(:),INTENT(IN)  :: PRT    ! Total mixing ratios (or tendency) to
                                             ! transform into rv,rc and ri
                                             ! (or its tendency)
REAL, DIMENSION(:), INTENT(OUT):: PTH    ! th (or th_l tendency)
REAL, DIMENSION(:), INTENT(OUT):: PRV    ! vapor mixing ratio (or its tendency)
REAL, DIMENSION(:), INTENT(OUT):: PRL    ! vapor mixing ratio (or its tendency)
REAL, DIMENSION(:), INTENT(OUT):: PRI    ! vapor mixing ratio (or its tendency)

!
!-------------------------------------------------------------------------------
!
!       0.2  declaration of local variables
!
INTEGER                       :: II
! Loop control
REAL                          :: ZCOEF
REAL, DIMENSION(SIZE(PP,1))   :: ZEXN,ZFOES,ZQVSAT
REAL, DIMENSION(SIZE(PTHL,1)) :: ZRVSAT,ZCPH,ZRLTEMP  
REAL, DIMENSION(SIZE(PTHL,1)) :: ZT,ZLOVCPEXN,ZLOSCPEXN
!----------------------------------------------------------------------------
!
!*      1 Initialisation 
!         --------------
!
!
ZCOEF = 0.8
PRL(:)=0.
PRI(:)=0.
ZRLTEMP(:)=0.
PRV(:)=PRT(:)
PTH(:)=PTHL(:)
ZEXN(:)=(PP(:)/XP00) ** (XRD/XCPD)
!
!       2 Iteration
!         ---------

    DO II=1,20
      ZT(:)=PTH(:)*ZEXN(:)

      WHERE (ZT(:) > 273.15)
! warm cloud
      ZFOES(:) = EXP( XALPW - XBETAW/ZT(:) - XGAMW*LOG(ZT(:))  )
      ZQVSAT(:) = XRD/XRV*ZFOES(:)/PP(:)   &
                   / (1.+(XRD/XRV-1.)*ZFOES(:)/PP(:))

      ZRVSAT(:) = (1-ZCOEF)*ZQVSAT(:)*(1+PRT(:))+(ZCOEF)*PRV(:)
!      CALL COMPUTE_FRAC_ICE_1D(HFRAC_ICE,PFRAC_ICE,PP,ZT) 
      PFRAC_ICE(:)=0.
      ZRLTEMP(:)=MAX(0.,PRV(:)-ZRVSAT(:))
      PRV(:)=PRV(:)-ZRLTEMP(:)
      PRL(:)=PRL(:)+PRI(:)+ZRLTEMP(:)
      PRI(:)    = PFRAC_ICE(:)    * (PRL(:))
      PRL(:)    = (1-PFRAC_ICE(:))* (PRT(:) - PRV(:))
      !       2.1 Cph
      ZCPH(:)=XCPD+ XCPV * PRV(:)+ XCL * PRL(:) + XCI * PRI(:)
      
      !       2.2 L/Cp/EXN
      ZLOVCPEXN(:) = (XLVTT + (XCPV-XCL) * (ZT(:)-XTT)) &
                        /(ZCPH*ZEXN(:))
      ZLOSCPEXN(:) = (XLSTT + (XCPV-XCI) * (ZT(:)-XTT)) &
                        /(ZCPH*ZEXN(:))
      PTH(:)=PTHL(:)+ZLOVCPEXN*PRL(:)+ZLOSCPEXN(:)*PRI(:)          

      ELSEWHERE
             ! cold shallow cloud not yet coded 
             ! probleme also for the highest level of fire case
        PRL(:)=0.
        PRI(:)=0.
        PRV(:)=PRT(:)
        PTH(:)=PTHL(:)
      ENDWHERE 
    ENDDO
   
END SUBROUTINE TH_R_FROM_THL_RT_1D

              
!     #################################################################
      SUBROUTINE TH_R_FROM_THL_RT_2D(HFRAC_ICE,PFRAC_ICE,PP,             &
                                  PTHL, PRT, PTH, PRV, PRL, PRI       )
!     #################################################################
!
!
!!****  *TH_R_FROM_THL_RT_2D* - computes the non-conservative variables
!!                          from conservative variables
!!
!!
!!
!! --------------------------------------------------------------------------


!
IMPLICIT NONE
!
!
!*      0.1  declarations of arguments
!
CHARACTER*1         , INTENT(IN) :: HFRAC_ICE
REAL, DIMENSION(:,:), INTENT(INOUT) :: PFRAC_ICE
REAL, DIMENSION(:,:), INTENT(IN) :: PP          ! Pressure
REAL, DIMENSION(:,:), INTENT(IN) :: PTHL    ! thetal (or thetal tendency) to
                                             ! transform into th (or tendency)
REAL, DIMENSION(:,:),INTENT(IN)  :: PRT    ! Total mixing ratios (or tendency) to
                                             ! transform into rv,rc and ri
                                             ! (or its tendency)
REAL, DIMENSION(:,:), INTENT(OUT):: PTH    ! th (or th_l tendency)
REAL, DIMENSION(:,:), INTENT(OUT):: PRV    ! vapor mixing ratio (or its tendency)
REAL, DIMENSION(:,:), INTENT(OUT):: PRL    ! vapor mixing ratio (or its tendency)
REAL, DIMENSION(:,:), INTENT(OUT):: PRI    ! vapor mixing ratio (or its tendency)

!
!-------------------------------------------------------------------------------
!
!       0.2  declaration of local variables
!
INTEGER                :: II
! Loop control
REAL                                         :: ZCOEF
REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2))       :: ZEXN,ZFOES,ZQVSAT
REAL, DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2))   :: ZRVSAT,ZCPH,ZRLTEMP  
REAL, DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2))   :: ZT,ZLOVCPEXN,ZLOSCPEXN
!----------------------------------------------------------------------------
!
!*      1 Initialisation 
!         --------------
!
!
ZCOEF = 0.8
PRL(:,:)=0.
PRI(:,:)=0.
ZRLTEMP(:,:)=0.
PRV(:,:)=PRT(:,:)
PTH(:,:)=PTHL(:,:)
ZEXN(:,:)=(PP(:,:)/XP00) ** (XRD/XCPD)
!
!       2 Iteration
!         ---------

    DO II=1,20
      ZT(:,:)=PTH(:,:)*ZEXN(:,:)
      WHERE (ZT(:,:) > 273.15)
! warm cloud

      ZFOES(:,:) = EXP( XALPW - XBETAW/ZT(:,:) - XGAMW*LOG(ZT(:,:))  )
      ZQVSAT(:,:) = XRD/XRV*ZFOES(:,:)/PP(:,:)   &
                   / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:))
      ZRVSAT(:,:) = (1-ZCOEF)*ZQVSAT(:,:)*(1+PRT(:,:))+(ZCOEF)*PRV(:,:)

!      CALL COMPUTE_FRAC_ICE_2D(HFRAC_ICE,PFRAC_ICE,PP,ZT) 
      PFRAC_ICE(:,:) = 0.
      ZRLTEMP(:,:)=MAX(0.,PRV(:,:)-ZRVSAT(:,:))
      PRV(:,:)=PRV(:,:)-ZRLTEMP(:,:)
      PRL(:,:)=PRL(:,:)+PRI(:,:)+ZRLTEMP(:,:)
      PRI(:,:)    = PFRAC_ICE(:,:)    * (PRL(:,:))
      PRL(:,:)    = (1-PFRAC_ICE(:,:))* (PRT(:,:) - PRV(:,:))
      !       2.1 Cph
        ZCPH(:,:)=XCPD+ XCPV * PRV(:,:)+ XCL * PRL(:,:) + XCI * PRI(:,:)
      
        !       2.2 L/Cp/EXN
        ZLOVCPEXN(:,:) = (XLVTT + (XCPV-XCL) * (ZT(:,:)-XTT)) &
                        /(ZCPH*ZEXN(:,:))
        ZLOSCPEXN(:,:) = (XLSTT + (XCPV-XCI) * (ZT(:,:)-XTT)) &
                        /(ZCPH*ZEXN(:,:))
        PTH(:,:)=PTHL(:,:)+ZLOVCPEXN*PRL(:,:)+ZLOSCPEXN(:,:)*PRI(:,:)          
      ELSEWHERE
        ! cold shallow cloud not yet coded 
        ! probleme also for the highest level of fire case
        PRL(:,:)=0.
        PRI(:,:)=0.
        PRV(:,:)=PRT(:,:)
        PTH(:,:)=PTHL(:,:)
       ENDWHERE   
    ENDDO
   
END SUBROUTINE TH_R_FROM_THL_RT_2D


                                       
!     #################################################################
      SUBROUTINE THL_RT_FROM_TH_R_MF( KRR,KRRL,KRRI,                  &
                                       PTH, PR, PLVOCPEXN, PLSOCPEXN, &
                                       PTHL, PRT                      )
!     #################################################################
!
!!
!!****  *THL_RT_FROM_TH_R* - computes the conservative variables THL and RT
!!                           from TH and the non precipitating water species
!!
!! --------------------------------------------------------------------------
!       
!*      0. DECLARATIONS
!          ------------
!
!USE MODD_CST
!
IMPLICIT NONE
!
!
!*      0.1  declarations of arguments
!
INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.

REAL, DIMENSION(:,:,:), INTENT(IN)   :: PTH      ! theta
REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR       ! water species
REAL, DIMENSION(:,:,:), INTENT(IN)   :: PLVOCPEXN, PLSOCPEXN      ! L/(cp*Pi)

REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHL     ! th_l
REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PRT      ! total non precip. water
!
!-------------------------------------------------------------------------------
!
!       0.2  declaration of local variables
!

!----------------------------------------------------------------------------

!----------------------------------------------------------------------------
!
!

IF ( KRRL >= 1 ) THEN
  IF ( KRRI >= 1 ) THEN 
    ! Rnp 
    PRT(:,:,:)  = PR(:,:,:,1)  + PR(:,:,:,2)  + PR(:,:,:,4)
    ! Theta_l 
    PTHL(:,:,:)  = PTH(:,:,:)  - PLVOCPEXN(:,:,:) * PR(:,:,:,2) &
                               - PLSOCPEXN(:,:,:) * PR(:,:,:,4)
  ELSE
    ! Rnp
    PRT(:,:,:)  = PR(:,:,:,1)  + PR(:,:,:,2) 
    ! Theta_l
    PTHL(:,:,:) = PTH(:,:,:)  - PLVOCPEXN(:,:,:) * PR(:,:,:,2)
  END IF
ELSE
    ! Rnp = rv
PRT(:,:,:)  = PR(:,:,:,1)
    ! Theta_l = Theta
PTHL(:,:,:) = PTH(:,:,:)
END IF

END SUBROUTINE THL_RT_FROM_TH_R_MF


!      #################################################
       SUBROUTINE TRIDIAG_MASSFLUX(PVARM,PF,PDFDT,PTSTEP,PIMPL,  &
                                 PDZZ,PRHODJ,PVARP             )
!      #################################################
!
!
!!****   *TRIDIAG_MASSFLUX* - routine to solve a time implicit scheme
!!
!!
!! ---------------------------------------------------------------------
!
!*       0. DECLARATIONS
!
!
!
IMPLICIT NONE
!
!
!*       0.1 declarations of arguments
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM   ! variable at t-1      at mass point
REAL, DIMENSION(:,:,:), INTENT(IN) :: PF      ! flux in dT/dt=-dF/dz at flux point
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDFDT   ! dF/dT                at flux point
REAL,                   INTENT(IN) :: PTSTEP  ! Double time step
REAL,                   INTENT(IN) :: PIMPL   ! implicit weight
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ    ! Dz                   at flux point
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ  ! (dry rho)*J          at mass point
!
REAL, DIMENSION(:,:,:), INTENT(OUT):: PVARP   ! variable at t+1      at mass point
!
!
!*       0.2 declarations of local variables
!
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZRHODJ_DFDT_O_DZ
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZMZM_RHODJ
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZA, ZB, ZC
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZY ,ZGAM 
                                         ! RHS of the equation, 3D work array
REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))                :: ZBET
                                         ! 2D work array
INTEGER                              :: JK            ! loop counter
INTEGER                              :: IKB,IKE       ! inner vertical limits
!
! ---------------------------------------------------------------------------
!                                              
!*      1.  Preliminaries
!           -------------
!
IKB=1
IKE=SIZE(PVARM,3)-1 
!
ZMZM_RHODJ  = MZM(PRHODJ)
ZRHODJ_DFDT_O_DZ = ZMZM_RHODJ*PDFDT/PDZZ
!
ZA=0.
ZB=0.
ZC=0.
ZY=0.
!
!
!*      2.  COMPUTE THE RIGHT HAND SIDE
!           ---------------------------
!
ZY(:,:,IKB) = PRHODJ(:,:,IKB)*PVARM(:,:,IKB)/PTSTEP             &
    - ZMZM_RHODJ(:,:,IKB+1) * PF(:,:,IKB+1)/PDZZ(:,:,IKB+1)     &
    + ZMZM_RHODJ(:,:,IKB  ) * PF(:,:,IKB  )/PDZZ(:,:,IKB  )     &
    + ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL * PVARM(:,:,IKB+1)    &
    + ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL * PVARM(:,:,IKB  )
!
DO JK=IKB+1,IKE-1
  ZY(:,:,JK) = PRHODJ(:,:,JK)*PVARM(:,:,JK)/PTSTEP          &
    - ZMZM_RHODJ(:,:,JK+1) * PF(:,:,JK+1)/PDZZ(:,:,JK+1)    &
    + ZMZM_RHODJ(:,:,JK  ) * PF(:,:,JK  )/PDZZ(:,:,JK  )    &
    + ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL * PVARM(:,:,JK+1)  &
    + ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL * PVARM(:,:,JK  )  &
    - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL * PVARM(:,:,JK  )  &
    - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL * PVARM(:,:,JK-1)
END DO
! 
ZY(:,:,IKE) = PRHODJ(:,:,IKE)*PVARM(:,:,IKE)/PTSTEP         &
    - ZMZM_RHODJ(:,:,IKE+1) * PF(:,:,IKE+1)/PDZZ(:,:,IKE+1) &
    + ZMZM_RHODJ(:,:,IKE  ) * PF(:,:,IKE  )/PDZZ(:,:,IKE  ) &
    - ZRHODJ_DFDT_O_DZ(:,:,IKE ) * 0.5*PIMPL * PVARM(:,:,IKE  ) &
    - ZRHODJ_DFDT_O_DZ(:,:,IKE ) * 0.5*PIMPL * PVARM(:,:,IKE-1)
!
!
!*       3.  INVERSION OF THE TRIDIAGONAL SYSTEM
!            -----------------------------------
!
IF ( PIMPL > 1.E-10 ) THEN
!
!*       3.1 arrays A, B, C
!            --------------
!
  ZB(:,:,IKB) =   PRHODJ(:,:,IKB)/PTSTEP                   &
                + ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL
  ZC(:,:,IKB) =   ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL

  DO JK=IKB+1,IKE-1
    ZA(:,:,JK) = - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL
    ZB(:,:,JK) =   PRHODJ(:,:,JK)/PTSTEP                   &
                 + ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL &
                 - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL
    ZC(:,:,JK) =   ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL
  END DO

  ZA(:,:,IKE) = - ZRHODJ_DFDT_O_DZ(:,:,IKE  ) * 0.5*PIMPL
  ZB(:,:,IKE) =   PRHODJ(:,:,IKE)/PTSTEP                   &
                - ZRHODJ_DFDT_O_DZ(:,:,IKE  ) * 0.5*PIMPL
!
!*       3.2 going up
!            --------
!
  ZBET(:,:) = ZB(:,:,IKB)  ! bet = b(ikb)
  PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)

  !
  DO JK = IKB+1,IKE-1
    ZGAM(:,:,JK) = ZC(:,:,JK-1) / ZBET(:,:)  
                                                    ! gam(k) = c(k-1) / bet
    ZBET(:,:)    = ZB(:,:,JK) - ZA(:,:,JK) * ZGAM(:,:,JK)
                                                    ! bet = b(k) - a(k)* gam(k)  
    PVARP(:,:,JK)= ( ZY(:,:,JK) - ZA(:,:,JK) * PVARP(:,:,JK-1) ) / ZBET(:,:)
                                        ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
  END DO 
  ! special treatment for the last level
  ZGAM(:,:,IKE) = ZC(:,:,IKE-1) / ZBET(:,:) 
                                                    ! gam(k) = c(k-1) / bet
  ZBET(:,:)     = ZB(:,:,IKE) - ZA(:,:,IKE) * ZGAM(:,:,IKE)
                                                    ! bet = b(k) - a(k)* gam(k)  
  PVARP(:,:,IKE)= ( ZY(:,:,IKE) - ZA(:,:,IKE) * PVARP(:,:,IKE-1) ) / ZBET(:,:)
                                       ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
!
!*       3.3 going down
!            ----------
!
  DO JK = IKE-1,IKB,-1
    PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+1) * PVARP(:,:,JK+1)
  END DO
!
!
ELSE
!!! EXPLICIT FORMULATION
! 
  PVARP(:,:,IKB:IKE) = ZY(:,:,IKB:IKE) * PTSTEP / PRHODJ(:,:,IKB:IKE)
!
END IF 
!
!
!*       4.  FILL THE UPPER AND LOWER EXTERNAL VALUES
!            ----------------------------------------
!
!PVARP(:,:,IKB-1)=PVARP(:,:,IKB) MODIF WRF JP
PVARP(:,:,IKE+1)=PVARP(:,:,IKE)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE TRIDIAG_MASSFLUX

!     
!     #################################################################
      SUBROUTINE UPDRAFT_SOPE(KRR,KRRL,KRRI,OMIXUV,                   & 
                         PZZ,PDZZ,PSFTH,PSFRV,PPABSM,PRHODREF,        &
                         PTKEM,PTHM,PRM,PTHLM,PRTM,PUM,PVM,           &
                         PTHL_UP,PRT_UP,PRV_UP,PU_UP,PV_UP,           &
                         PRC_UP,PRI_UP,PTHV_UP,PW_UP,PFRAC_UP,PEMF,   &
                         PDETR,PENTR,KKLCL,KKETL,KKCTL )
!     #################################################################
!!
!!****  *UPDRAFT_SOPE* - Interfacing routine
!!
!! --------------------------------------------------------------------------
!
!*      0. DECLARATIONS
!          ------------
!

IMPLICIT NONE

!*                    1.1  Declaration of Arguments
!
!
INTEGER,                INTENT(IN)   :: KRR       ! number of moist var.
INTEGER,                INTENT(IN)   :: KRRL      ! number of liquid water var.
INTEGER,                INTENT(IN)   :: KRRI      ! number of ice water var.
LOGICAL,                INTENT(IN)   :: OMIXUV    ! True if mixing of momentum

REAL, DIMENSION(:,:,:), INTENT(IN)   :: PZZ       !  Height at the flux point

REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDZZ      ! depth between mass levels 
 
REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTH,PSFRV
                                            ! normal surface fluxes of theta,rv
!
!    prognostic variables at t- deltat
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PRHODREF    ! dry density of the
                                                     ! reference state
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! TKE
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUM,PVM     ! momentum

!
!   thermodynamical variables which are transformed in conservative var.
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PTHM       ! pot. temp. = PTHLM in turb.f90
REAL, DIMENSION(:,:,:,:), INTENT(IN)   ::  PRM        ! water species 
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PTHLM,PRTM !cons. var.  
REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PTHL_UP,PRT_UP  ! updraft properties
REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PRV_UP,PRC_UP,PRI_UP,&!Thl,Rt,Rv,Rc,Ri
                                           PW_UP,PFRAC_UP,PEMF, &!w,Updraft Fraction, Mass Flux
                                           PDETR,PENTR,PTHV_UP, &!entrainment, detrainment, ThV
                                           PU_UP, PV_UP          !updraft wind component
                                           
INTEGER, DIMENSION(:,:),  INTENT(OUT)  ::  KKLCL,KKETL,KKCTL     !index for LCL,ETL,CTL       
!
!
!
!                       1.2  Declaration of local variables

INTEGER :: IKB

!
! Variables to transform 3D fields in 2D fields
REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZPABSM,ZRHODREF
       
REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZZZ,ZDZZ

REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZTHLM,ZRTM,&
                                                              ZTHM,ZTKEM,&
                                                              ZUM,ZVM
                                                              
REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZRVM,ZRCM,ZRIM

REAL, DIMENSION(SIZE(PSFTH,1)*SIZE(PSFTH,2))    ::  ZSFTH,ZSFRV

REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3))  :: ZTHL_UP,ZRT_UP,&
                                      ZRV_UP,ZRC_UP, ZRI_UP, &
                                      ZW_UP,ZU_UP,ZV_UP,ZTHV_UP,           &
                                      ZFRAC_UP,ZEMF_UP,ZENTR_UP,ZDETR_UP

REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2))    :: ZFRAC_GRID
INTEGER, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2)) :: JKETL,JKCTL,JKLCL

INTEGER :: IIU,IJU,IKU! Limit of the grid
INTEGER :: J1D        ! horizontal loop counter
INTEGER :: JK,JKK,JSV ! loop counters
INTEGER :: JRR        ! moist loop counter
REAL    :: ZRVORD     ! Rv/Rd (=1/0.622 cf glossaire)

!------------------------------------------------------------------------

!                     2. INITIALISATION

!                     2.1 Local variables, internal domain

IIU=SIZE(PTKEM,1)
IJU=SIZE(PTKEM,2)
!
IKB = 1              ! Modif WRF JP
IKU = SIZE(PTKEM,3)

ZRVORD = XRV / XRD


DO JK=IKB,IKU
  ZZZ    (:,JK) = RESHAPE(PZZ    (:,:,JK),(/ IIU*IJU /) )
  ZDZZ   (:,JK) = RESHAPE(PDZZ    (:,:,JK),(/ IIU*IJU /) )  
  ZTHM   (:,JK) = RESHAPE(PTHM   (:,:,JK),(/ IIU*IJU /) )
  ZTKEM  (:,JK) = RESHAPE(PTKEM  (:,:,JK),(/ IIU*IJU /) )
  ZPABSM (:,JK) = RESHAPE(PPABSM (:,:,JK),(/ IIU*IJU /) )
  ZRHODREF(:,JK) = RESHAPE(PRHODREF(:,:,JK),(/ IIU*IJU /) )  
  ZRVM   (:,JK) = RESHAPE(PRM    (:,:,JK,1),(/ IIU*IJU /) ) 
  ZTHLM  (:,JK) = RESHAPE(PTHLM   (:,:,JK),(/ IIU*IJU /) )
  ZRTM   (:,JK) = RESHAPE(PRTM    (:,:,JK),(/ IIU*IJU /) )
  ZUM    (:,JK) = RESHAPE(PUM    (:,:,JK),(/ IIU*IJU /) )
  ZVM    (:,JK) = RESHAPE(PVM    (:,:,JK),(/ IIU*IJU /) )
END DO

IF (KRRL>1) THEN
DO JK=1,IKU
  ZRCM   (:,JK) = RESHAPE(PRM    (:,:,JK,2),(/ IIU*IJU /) ) 
END DO
ELSE
  ZRCM   (:,:) =0.
ENDIF
IF (KRRI>1) THEN
DO JK=1,IKU
  ZRIM   (:,JK) = RESHAPE(PRM    (:,:,JK,4),(/ IIU*IJU /) ) 
END DO
ELSE
  ZRIM   (:,:) =0.
ENDIF

ZSFTH(:)=RESHAPE(PSFTH(:,:),(/ IIU*IJU /) )
ZSFRV(:)=RESHAPE(PSFRV(:,:),(/ IIU*IJU /) )

!Updraft begins at level 1 (Modif WRF)
JK=IKB 

!        6.2 compute properties of the updraft
CALL COMPUTE_UPDRAFT(OMIXUV,ZZZ,ZDZZ,JK,      &
                     ZSFTH,ZSFRV,ZPABSM,ZRHODREF,ZUM,ZVM,ZTKEM, &
                     ZTHM,ZRVM,ZRCM,ZRIM,ZTHLM,ZRTM,       &
                     ZTHL_UP,ZRT_UP,ZRV_UP,ZRC_UP,ZRI_UP,&
                     ZTHV_UP,ZW_UP,ZU_UP,ZV_UP, &
                     ZFRAC_UP,ZEMF_UP,&
                     ZDETR_UP,ZENTR_UP,&
                     JKLCL,JKETL,JKCTL)


PTHL_UP(:,:,:)= RESHAPE(ZTHL_UP(:,:), (/ IIU,IJU,IKU /))
PRT_UP(:,:,:)=RESHAPE(ZRT_UP(:,:), (/ IIU,IJU,IKU /) )
PRV_UP(:,:,:)=RESHAPE(ZRV_UP(:,:), (/ IIU,IJU,IKU /) )
PRC_UP(:,:,:)=RESHAPE(ZRC_UP(:,:), (/ IIU,IJU,IKU /) )
PRI_UP(:,:,:)=RESHAPE(ZRI_UP(:,:), (/ IIU,IJU,IKU /) )
PW_UP(:,:,:)=RESHAPE(ZW_UP(:,:), (/ IIU,IJU,IKU /) )
PU_UP(:,:,:)=RESHAPE(ZU_UP(:,:), (/ IIU,IJU,IKU /) )
PV_UP(:,:,:)=RESHAPE(ZV_UP(:,:), (/ IIU,IJU,IKU /) )
PEMF(:,:,:)=RESHAPE(ZEMF_UP(:,:), (/ IIU,IJU,IKU /) )
PDETR(:,:,:)=RESHAPE(ZDETR_UP(:,:), (/ IIU,IJU,IKU /) )
PENTR(:,:,:)=RESHAPE(ZENTR_UP(:,:), (/ IIU,IJU,IKU /) )
PTHV_UP(:,:,:)=RESHAPE(ZTHV_UP(:,:), (/ IIU,IJU,IKU /) )
KKETL(:,:)=RESHAPE(JKETL(:),(/ IIU,IJU/) )
KKCTL(:,:)=RESHAPE(JKCTL(:),(/ IIU,IJU/) )
KKLCL(:,:)=RESHAPE(JKLCL(:),(/ IIU,IJU/) )
PFRAC_UP(:,:,:)=RESHAPE(ZFRAC_UP(:,:),(/ IIU,IJU,IKU /) )

END SUBROUTINE UPDRAFT_SOPE

!     #################################################################
      SUBROUTINE COMPUTE_MF_CLOUD(KRRL, PTHLM,PRC_UP, PFRAC_UP, PDZZ, KKLCL,&
                                  PRC_MF,PCF_MF                     )
!     #################################################################
!!
!!****  *COMPUTE_MF_CLOUD* - 
!!       compute diagnostic subgrid cumulus cloud caracteristics
!!
!!    PURPOSE
!!    -------
!!****  The purpose of this routine is to compute the cloud fraction and 
!!      the mean cloud content associated with clouds described by the 
!!      mass flux scheme
!!
!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!      
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!     REFERENCE
!!     ---------
!!
!!
!!     AUTHOR
!!     ------
!! --------------------------------------------------------------------------
!
!*      0. DECLARATIONS
!          ------------
!


IMPLICIT NONE

!*                    1.1  Declaration of Arguments
!
!
!
INTEGER,                  INTENT(IN)   ::  KRRL         ! number of liquid water var.
                                                        ! scheme
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PTHLM ! updraft characteritics

REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRC_UP ! updraft characteritics
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PFRAC_UP          ! Updraft Fraction

REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PDZZ
INTEGER, DIMENSION(:,:),  INTENT(IN)   ::  KKLCL          ! index of updraft condensation level
REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PRC_MF, PCF_MF ! cloud content and
                                                          ! cloud fraction for MF scheme

!
!                       1.2  Declaration of local variables
!
REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::  ZFLXZ
INTEGER  :: IKU, IKB, IKE, JI,JJ,JK
!------------------------------------------------------------------------

!                     1. INITIALISATION

!                     2.1 Internal domain

! Internal Domain
IKU=SIZE(PTHLM,3)
IKB=1
IKE=IKU-1

PCF_MF = 0.
PRC_MF = 0.


!
!                    Direct cloud scheme
!        This scheme may be activated only if the selected updraft model
!        gives the updraft fraction as an output
!        

! attention, les variables de l'updraft sont en niveaux flux
! ils faut les passer aux niveaux masse pour calculer PRC_MF et PCF_MF
DO JI=1,SIZE(PCF_MF,1)
DO JJ=1,SIZE(PCF_MF,2)

DO JK=KKLCL(JI,JJ),IKE

PCF_MF(JI,JJ,JK ) = XKCF_MF *0.5* (       &
  &    PFRAC_UP(JI,JJ,JK) +  PFRAC_UP(JI,JJ,JK+1) )
PRC_MF(JI,JJ,JK)  = 0.5* XKCF_MF * ( PFRAC_UP(JI,JJ,JK)*PRC_UP(JI,JJ,JK)  &
                         + PFRAC_UP(JI,JJ,JK+1)*PRC_UP(JI,JJ,JK+1) )


END DO
END DO
END DO


END SUBROUTINE COMPUTE_MF_CLOUD


!     #################################################################
      SUBROUTINE mfshconvpblinit(massflux_EDKF, entr_EDKF, detr_EDKF & 
                                     ,thl_up, thv_up, rt_up       &
                                     ,rv_up, rc_up, u_up, v_up    &
                                     ,frac_up,RESTART,ALLOWED_TO_READ,     &
     &                      IDS,IDE,JDS,JDE,KDS,KDE,                    &
     &                      IMS,IME,JMS,JME,KMS,KME,                    &
     &                      ITS,ITE,JTS,JTE,KTS,KTE                 )
!     #################################################################

      IMPLICIT NONE

      LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
     &                      IMS,IME,JMS,JME,KMS,KME,                    &
     &                      ITS,ITE,JTS,JTE,KTS,KTE

      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT),OPTIONAL ::             &
     &                                        MASSFLUX_EDKF, ENTR_EDKF, DETR_EDKF &
     &                                       ,THL_UP, THV_UP, RT_UP, RV_UP        &
     &                                       ,RC_UP, U_UP, V_UP, FRAC_UP 

      INTEGER :: I,J,K,ITF,JTF,KTF

!---------------------------------------------
      JTF=MIN0(JTE,JDE-1)
      KTF=MIN0(KTE,KDE-1)
      ITF=MIN0(ITE,IDE-1)

!     IF(.NOT.RESTART)THEN
      IF( PRESENT (MASSFLUX_EDKF) ) THEN
        DO J=JTS,JTF
        DO K=KTS,KTF
        DO I=ITS,ITF
          MASSFLUX_EDKF(I,K,J)=0.
          ENTR_EDKF(I,K,J)=0.
          DETR_EDKF(I,K,J)=0.
          THL_UP(I,K,J)=0.
          THV_UP(I,K,J)=0.
          RT_UP(I,K,J)=0.
          RV_UP(I,K,J)=0.
          RC_UP(I,K,J)=0.
          U_UP(I,K,J)=0.
          V_UP(I,K,J)=0.
          FRAC_UP(I,K,J)=0.
        ENDDO
        ENDDO
        ENDDO
      ENDIF

END SUBROUTINE mfshconvpblinit


!     #########################################################
      SUBROUTINE BL89(PZZ,PDZZ,PTHVREF,PTHLM,KRR, &
                 PRM,PTKEM,PLM)
!     #########################################################
!
!!****  *BL89* -
!!
!!    PURPOSE
!!    -------
!!    This routine computes the mixing length from Bougeault-Lacarrere 89
!!    formula.
!!
!!**  METHOD
!!    ------
!!
!*       0.1   Declaration of arguments
!              ------------------------
!
INTEGER,                INTENT(IN)   :: KRR

REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PZZ
REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PDZZ
REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTHVREF
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! TKE
!   thermodynamical variables PTHLM=Theta at the begining
REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTHLM     ! conservative pot. temp.
REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM       ! water var.                 
REAL, DIMENSION(:,:,:),   INTENT(OUT) :: PLM       ! Mixing length

!*       0.2   Declaration of local variables
!              ------------------------------
!
INTEGER :: IKB,IKE

REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZRTM
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)) ::  ZLWORK
!           ! downwards then upwards vertical displacement,
REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZZZ,ZDZZ,&
                                                              ZG_O_THVREF, &
                                                              ZTHM,ZTKEM,ZLM,&
                                                              ZLMUP,ZLMDN,&
                                                              ZLMTEST
!           ! 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

REAL   :: ZPOTE,ZLWORK1,ZLWORK2

INTEGER :: IIU,IJU,IKU,IICE
INTEGER :: J1D        ! horizontal loop counter
INTEGER :: JK,JKK     ! loop counters
INTEGER :: JRR        ! moist loop counter
REAL    :: ZRVORD     ! Rv/Rd (=1/0.622 cf glossaire)
LOGICAL :: GUPORDN    !Test for computation of upward or downward mixing length
REAL    :: Z2SQRT2
!-------------------------------------------------------------------------------
!
Z2SQRT2=2.*SQRT(2.)
IIU=SIZE(PTKEM,1)
IJU=SIZE(PTKEM,2)
!
IKB =  1
IKE = SIZE(PTKEM,3)-1
IKU = SIZE(PTKEM,3)
ZRVORD = XRV / XRD
!-------------------------------------------------------------------------------
!
!*       1.    pack the horizontal dimensions into one
!              ---------------------------------------
!
DO JK=1,IKU
  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
!-------------------------------------------------------------------------------
!
!*       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
!
!-------------------------------------------------------------------------------
!
!*       3.  loop on model levels
!            --------------------
!
DO JK=IKB,IKE
!-------------------------------------------------------------------------------
!
!*       4.  mixing length for a downwards displacement
!            ------------------------------------------
   GUPORDN=.FALSE.
   CALL COMPUTE_BL89_ML(ZDZZ,ZTKEM,ZG_O_THVREF,ZVPT,JK,GUPORDN,ZLWORK)



!-------------------------------------------------------------------------------
!
!*       5.  intermediate storage of the final mixing length
!            -----------------------------------------------
!
   ZLMDN(:,JK)=MIN(ZLWORK(:),0.5*(ZZZ(:,JK)+ZZZ(:,JK+1))-ZZZ(:,IKB))
!
!-------------------------------------------------------------------------------
!
!*       6.  mixing length for an upwards displacement
!            -----------------------------------------
 
   GUPORDN=.TRUE.
   CALL COMPUTE_BL89_ML(ZDZZ,ZTKEM,ZG_O_THVREF,ZVPT,JK,GUPORDN,ZLWORK)

   ZLMUP(:,JK)=ZLWORK(:)

!-------------------------------------------------------------------------------
! 
  DO J1D=1,IIU*IJU
    ZLWORK1=MAX(ZLMDN(J1D,JK),1.E-10)
    ZLWORK2=MAX(ZLMUP(J1D,JK),1.E-10)
    ZPOTE = ZLWORK1 / ZLWORK2
    ZLWORK2=1.d0 + ZPOTE**(2./3.)
    ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2))
  END DO

ZLM(:,JK)=( 0.5* (MAX(ZLMDN(:,JK),1.E-10)**(-2./3.)+MAX(ZLMUP(:,JK),1.E-10)**(-2./3.)) )**(-1.5)
ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)

!*       8.  end of the loop on the vertical levels
!            --------------------------------------
!
END DO
!
!-------------------------------------------------------------------------------
!
!*       9.  boundaries
!            ----------
!
ZLM(:,IKE)  =ZLM(:,IKE-1)
ZLM(:,IKE+1)=ZLM(:,IKE-1)
ZLMUP(:,IKE)  =ZLMUP(:,IKE-1)
ZLMUP(:,IKE+1)=ZLMUP(:,IKE-1)
ZLMDN(:,IKE)  =ZLMDN(:,IKE-1)
ZLMDN(:,IKE+1)=ZLMDN(:,IKE-1)
!
  DO JK=1,IKU
    PLM  (:,:,JK)   = RESHAPE(ZLM  (:,JK), (/ IIU,IJU /) )
  END DO

END SUBROUTINE BL89
                        
END  MODULE MODULE_BL_MFSHCONVPBL

