MODULE lmdz_cloudth


  IMPLICIT NONE

CONTAINS

       SUBROUTINE cloudth(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,t, &
             cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)


      use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs

      IMPLICIT NONE


!===========================================================================
! Auteur : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================


      INCLUDE "YOMCST.h"
      INCLUDE "YOETHF.h"
      INCLUDE "FCTTRE.h"

      INTEGER itap,ind1,ind2
      INTEGER ngrid,klev,klon,l,ig 
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
      
      REAL ztv(ngrid,klev)
      REAL po(ngrid)
      REAL zqenv(ngrid)   
      REAL zqta(ngrid,klev)
          
      REAL fraca(ngrid,klev+1)
      REAL zpspsk(ngrid,klev)
      REAL paprs(ngrid,klev+1)
      REAL pplay(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zthl(ngrid,klev)

      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      
      REAL sigma1(ngrid,klev)
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)  
      REAL cenv(ngrid,klev)   
      REAL ctot(ngrid,klev)
      REAL rneb(ngrid,klev)
      REAL t(ngrid,klev)
      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
      REAL rdd,cppd,Lv
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,xth,xenv
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
      
      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL zqs(ngrid), qcloud(ngrid)
      REAL erf




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Gestion de deux versions de cloudth
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      IF (iflag_cloudth_vert>=1) THEN
      CALL cloudth_vert(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,t)
      RETURN
      ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!-------------------------------------------------------------------------------
! Initialisation des variables r?elles
!-------------------------------------------------------------------------------
      sigma1(:,:)=0.
      sigma2(:,:)=0.
      qlth(:,:)=0.
      qlenv(:,:)=0.  
      qltot(:,:)=0.
      rneb(:,:)=0.
      qcloud(:)=0.
      cth(:,:)=0.
      cenv(:,:)=0.
      ctot(:,:)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)



!-------------------------------------------------------------------------------
! Calcul de la fraction du thermique et des ?cart-types des distributions
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))


!      zqenv(ind1)=po(ind1)
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef




      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 




      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
      ath=1./(1.+(alth*Lv/cppd))
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
      
      

!------------------------------------------------------------------------------
! Calcul des ?cart-types pour s
!------------------------------------------------------------------------------

!      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) 
!       if (paprs(ind1,ind2).gt.90000) then
!       ratqs(ind1,ind2)=0.002
!       else 
!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
!       endif
       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
!       sigma1s=ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
 
!------------------------------------------------------------------------------
! Calcul de l'eau condens?e et de la couverture nuageuse
!------------------------------------------------------------------------------
      sqrt2pi=sqrt(2.*pi)
      xth=sth/(sqrt(2.)*sigma2s)
      xenv=senv/(sqrt(2.)*sigma1s)
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      if (ctot(ind1,ind2)<1.e-10) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else   
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)

      endif                           
      
          
!     PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'


      else  ! gaussienne environnement seule
      
      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
      

!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
      

      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)

      sqrt2pi=sqrt(2.*pi)
      xenv=senv/(sqrt(2.)*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
      
      if (ctot(ind1,ind2)<1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else   
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)

      endif    
 
 
 
 
 
 
      endif   
      enddo
     
      return
!     end
END SUBROUTINE cloudth



!===========================================================================
     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,t)

!===========================================================================
! Auteur : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================


      use lmdz_cloudth_ini, only: iflag_cloudth_vert, vert_alpha

      IMPLICIT NONE

      INCLUDE "YOMCST.h"
      INCLUDE "YOETHF.h"
      INCLUDE "FCTTRE.h"
      
      INTEGER itap,ind1,ind2
      INTEGER ngrid,klev,klon,l,ig 
      
      REAL ztv(ngrid,klev)
      REAL po(ngrid)
      REAL zqenv(ngrid)   
      REAL zqta(ngrid,klev)
          
      REAL fraca(ngrid,klev+1)
      REAL zpspsk(ngrid,klev)
      REAL paprs(ngrid,klev+1)
      REAL pplay(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zthl(ngrid,klev)

      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      
      REAL sigma1(ngrid,klev)                                                         
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)  
      REAL cenv(ngrid,klev)   
      REAL ctot(ngrid,klev)
      REAL rneb(ngrid,klev)
      REAL t(ngrid,klev)                                                                  
      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
      REAL rdd,cppd,Lv,sqrt2,sqrtpi
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,xth,xenv
      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
      ! (J Jouhaud, JL Dufresne, JB Madeleine)
      
      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL zqs(ngrid), qcloud(ngrid)
      REAL erf

!------------------------------------------------------------------------------
! Initialisation des variables r?elles
!------------------------------------------------------------------------------
      sigma1(:,:)=0.
      sigma2(:,:)=0.
      qlth(:,:)=0.
      qlenv(:,:)=0.  
      qltot(:,:)=0.
      rneb(:,:)=0.
      qcloud(:)=0.
      cth(:,:)=0.
      cenv(:,:)=0.
      ctot(:,:)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(pi)

!-------------------------------------------------------------------------------
! Calcul de la fraction du thermique et des ?cart-types des distributions
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))


!      zqenv(ind1)=po(ind1)
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef




      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 




      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
      ath=1./(1.+(alth*Lv/cppd))
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2)) 
      
      

!------------------------------------------------------------------------------
! Calcul des ?cart-types pour s
!------------------------------------------------------------------------------

      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
!       if (paprs(ind1,ind2).gt.90000) then
!       ratqs(ind1,ind2)=0.002
!       else 
!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
!       endif
!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
!       sigma1s=ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
 
!------------------------------------------------------------------------------
! Calcul de l'eau condens?e et de la couverture nuageuse
!------------------------------------------------------------------------------
      sqrt2pi=sqrt(2.*pi)
      xth=sth/(sqrt(2.)*sigma2s)
      xenv=senv/(sqrt(2.)*sigma1s)
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     
       IF (iflag_cloudth_vert == 1) THEN
!-------------------------------------------------------------------------------
!  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
!-------------------------------------------------------------------------------
!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
!      deltasenv=aenv*0.01*po(ind1)
!     deltasth=ath*0.01*zqta(ind1,ind2)    
      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
      
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))

      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlenv(ind1,ind2)=IntJ 
!      PRINT*, qlenv(ind1,ind2),'VERIF EAU'


      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
!      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
!      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlth(ind1,ind2)=IntJ
!      PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      ELSE IF (iflag_cloudth_vert == 2) THEN

!-------------------------------------------------------------------------------
!  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
!-------------------------------------------------------------------------------
!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
      deltasenv=aenv*vert_alpha*sigma1s
      deltasth=ath*vert_alpha*sigma2s
      
      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
!     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
!     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
      
      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)

      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))

!      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
!      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
!      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))

      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlenv(ind1,ind2)=IntJ 
!      PRINT*, qlenv(ind1,ind2),'VERIF EAU'

      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
      
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
!      qlth(ind1,ind2)=IntJ
!      PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
      



      ENDIF ! of if (iflag_cloudth_vert==1 or 2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      if (cenv(ind1,ind2)<1.e-10.or.cth(ind1,ind2)<1.e-10) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else 
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) 

      endif  
                       
      
          
!     PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'


      else  ! gaussienne environnement seule
      
      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
      

!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2)) 
      

      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)

      sqrt2pi=sqrt(2.*pi)
      xenv=senv/(sqrt(2.)*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
      
      if (ctot(ind1,ind2)<1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 

      else   
                
      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)

      endif    
 
 
 
 
 
 
      endif   
      enddo
     
      return
!     end
END SUBROUTINE cloudth_vert




       SUBROUTINE cloudth_v3(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,t, &
             cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

      use lmdz_cloudth_ini, only: iflag_cloudth_vert

      IMPLICIT NONE


!===========================================================================
! Author : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================


      INCLUDE "YOMCST.h"
      INCLUDE "YOETHF.h"
      INCLUDE "FCTTRE.h"

      integer, intent(in) :: ind2
      integer, intent(in) :: ngrid,klev
      
      real, dimension(ngrid,klev), intent(in) :: ztv
      real, dimension(ngrid), intent(in) :: po
      real, dimension(ngrid,klev), intent(in) :: zqta
      real, dimension(ngrid,klev+1), intent(in) :: fraca
      real, dimension(ngrid), intent(out) :: qcloud
      real, dimension(ngrid,klev), intent(out) :: ctot
      real, dimension(ngrid,klev), intent(out) :: ctot_vol
      real, dimension(ngrid,klev), intent(in) :: zpspsk
      real, dimension(ngrid,klev+1), intent(in) :: paprs
      real, dimension(ngrid,klev), intent(in) :: pplay
      real, dimension(ngrid,klev), intent(in) :: ztla
      real, dimension(ngrid,klev), intent(inout) :: zthl
      real, dimension(ngrid,klev), intent(in) :: ratqs
      real, dimension(ngrid), intent(in) :: zqs
      real, dimension(ngrid,klev), intent(in) :: t
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv


      REAL zqenv(ngrid)   
      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      REAL sigma1(ngrid,klev)                                                         
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)
      REAL cenv(ngrid,klev)   
      REAL cth_vol(ngrid,klev)
      REAL cenv_vol(ngrid,klev)
      REAL rneb(ngrid,klev)      
      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
      REAL rdd,cppd,Lv
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL erf


      INTEGER :: ind1,l, ig

      IF (iflag_cloudth_vert>=1) THEN
      CALL cloudth_vert_v3(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,t, &
             cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
      RETURN
      ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!-------------------------------------------------------------------------------
! Initialisation des variables r?elles
!-------------------------------------------------------------------------------
      sigma1(:,:)=0.
      sigma2(:,:)=0.
      qlth(:,:)=0.
      qlenv(:,:)=0.  
      qltot(:,:)=0.
      rneb(:,:)=0.
      qcloud(:)=0.
      cth(:,:)=0.
      cenv(:,:)=0.
      ctot(:,:)=0.
      cth_vol(:,:)=0.
      cenv_vol(:,:)=0.
      ctot_vol(:,:)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(pi)


!-------------------------------------------------------------------------------
! Cloud fraction in the thermals and standard deviation of the PDFs
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))

      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef


      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84

!po = qt de l'environnement ET des thermique
!zqenv = qt environnement
!zqsatenv = qsat environnement
!zthl = Tl environnement


      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
      
!zqta = qt thermals
!zqsatth = qsat thermals
!ztla = Tl thermals

!------------------------------------------------------------------------------
! s standard deviations
!------------------------------------------------------------------------------

!     tests
!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
!     final option
      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
 
!------------------------------------------------------------------------------
! Condensed water and cloud cover
!------------------------------------------------------------------------------
      xth=sth/(sqrt2*sigma2s)
      xenv=senv/(sqrt2*sigma1s)
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
      ctot_vol(ind1,ind2)=ctot(ind1,ind2)

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      if (ctot(ind1,ind2)<1.e-10) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2) 
      else
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
      endif

      else  ! Environnement only, follow the if l.110
      
      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef

!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)  
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
      
      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)

      xenv=senv/(sqrt2*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))

      if (ctot(ind1,ind2)<1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)
      else   
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
      endif


      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
      enddo       ! from the loop on ngrid l.108
      return
!     end
END SUBROUTINE cloudth_v3



!===========================================================================
     SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,t, &
             cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

!===========================================================================
! Auteur : Arnaud Octavio Jam (LMD/CNRS)
! Date : 25 Mai 2010
! Objet : calcule les valeurs de qc et rneb dans les thermiques
!===========================================================================

      use lmdz_cloudth_ini, ONLY: iflag_cloudth_vert,iflag_ratqs
      use lmdz_cloudth_ini, ONLY: vert_alpha,vert_alpha_th, sigma1s_factor, sigma1s_power , sigma2s_factor , sigma2s_power , cloudth_ratqsmin , iflag_cloudth_vert_noratqs

      IMPLICIT NONE



      INCLUDE "YOMCST.h"
      INCLUDE "YOETHF.h"
      INCLUDE "FCTTRE.h"
      
      INTEGER itap,ind1,ind2
      INTEGER ngrid,klev,klon,l,ig 
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
      
      REAL ztv(ngrid,klev)
      REAL po(ngrid)
      REAL zqenv(ngrid)   
      REAL zqta(ngrid,klev)
          
      REAL fraca(ngrid,klev+1)
      REAL zpspsk(ngrid,klev)
      REAL paprs(ngrid,klev+1)
      REAL pplay(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zthl(ngrid,klev)

      REAL zqsatth(ngrid,klev)
      REAL zqsatenv(ngrid,klev)
      
      REAL sigma1(ngrid,klev)                                                         
      REAL sigma2(ngrid,klev)
      REAL qlth(ngrid,klev)
      REAL qlenv(ngrid,klev)
      REAL qltot(ngrid,klev) 
      REAL cth(ngrid,klev)
      REAL cenv(ngrid,klev)   
      REAL ctot(ngrid,klev)
      REAL cth_vol(ngrid,klev)
      REAL cenv_vol(ngrid,klev)
      REAL ctot_vol(ngrid,klev)
      REAL rneb(ngrid,klev)
      REAL t(ngrid,klev)                                                                  
      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
      REAL rdd,cppd,Lv
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
      REAL Tbef,zdelta,qsatbef,zcor
      REAL qlbef  
      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
      ! (J Jouhaud, JL Dufresne, JB Madeleine)

      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid) 
      REAL zqs(ngrid), qcloud(ngrid)
      REAL erf

      REAL rhodz(ngrid,klev)
      REAL zrho(ngrid,klev)
      REAL dz(ngrid,klev)

      DO ind1 = 1, ngrid
        !Layer calculation
        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg !kg/m2
        zrho(ind1,ind2) = pplay(ind1,ind2)/t(ind1,ind2)/rd !kg/m3
        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2) !m : epaisseur de la couche en metre
      END DO

!------------------------------------------------------------------------------
! Initialize
!------------------------------------------------------------------------------

      sigma1(:,:)=0.
      sigma2(:,:)=0.
      qlth(:,:)=0.
      qlenv(:,:)=0.  
      qltot(:,:)=0.
      rneb(:,:)=0.
      qcloud(:)=0.
      cth(:,:)=0.
      cenv(:,:)=0.
      ctot(:,:)=0.
      cth_vol(:,:)=0.
      cenv_vol(:,:)=0.
      ctot_vol(:,:)=0.
      qsatmmussig1=0.
      qsatmmussig2=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2pi=sqrt(2.*pi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(pi)



!-------------------------------------------------------------------------------
! Calcul de la fraction du thermique et des ecart-types des distributions
!-------------------------------------------------------------------------------                 
      do ind1=1,ngrid

      if ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) then !Thermal and environnement

      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv


      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef


      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84

!zqenv = qt environnement
!zqsatenv = qsat environnement
!zthl = Tl environnement


      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
            
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
      
      
!zqta = qt thermals
!zqsatth = qsat thermals
!ztla = Tl thermals
!------------------------------------------------------------------------------
! s standard deviation
!------------------------------------------------------------------------------

      sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
                  (1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
!     sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5
      IF (cloudth_ratqsmin>0.) THEN
         sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
      ELSE
         sigma1s_ratqs = ratqs(ind1,ind2)*po(ind1)
      ENDIF
      sigma1s = sigma1s_fraca + sigma1s_ratqs
      IF (iflag_ratqs==11) then
         sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
      ENDIF
      sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
!      tests
!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) 
!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
!       if (paprs(ind1,ind2).gt.90000) then
!       ratqs(ind1,ind2)=0.002
!       else 
!       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
!       endif
!       sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
!       sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
!       sigma1s=ratqs(ind1,ind2)*po(ind1)
!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
 
       IF (iflag_cloudth_vert == 1) THEN
!-------------------------------------------------------------------------------
!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
!-------------------------------------------------------------------------------

      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)

      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xth1=(sth-deltasth)/(sqrt(2.)*sigma2s)
      xth2=(sth+deltasth)/(sqrt(2.)*sigma2s)
      coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
      coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
      
      cth(ind1,ind2)=0.5*(1.+1.*erf(xth2))
      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv2)) 
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   

      ! Environment
      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
      IntI2=coeffqlenv*xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
      IntI3=coeffqlenv*0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))

      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3

      ! Thermal
      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      ELSE IF (iflag_cloudth_vert >= 3) THEN
      IF (iflag_cloudth_vert < 5) THEN
!-------------------------------------------------------------------------------
!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
!-------------------------------------------------------------------------------
!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
!      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
!      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
!      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
      IF (iflag_cloudth_vert == 3) THEN
        deltasenv=aenv*vert_alpha*sigma1s
        deltasth=ath*vert_alpha_th*sigma2s
      ELSE IF (iflag_cloudth_vert == 4) THEN
        IF (iflag_cloudth_vert_noratqs == 1) THEN
          deltasenv=vert_alpha*max(sigma1s_fraca,1e-10)
          deltasth=vert_alpha_th*sigma2s
        ELSE
          deltasenv=vert_alpha*sigma1s
          deltasth=vert_alpha_th*sigma2s
        ENDIF
      ENDIF
      
      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
      exp_xenv1 = exp(-1.*xenv1**2)
      exp_xenv2 = exp(-1.*xenv2**2)
      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
      exp_xth1 = exp(-1.*xth1**2)
      exp_xth2 = exp(-1.*xth2**2)
      
      !CF_surfacique
      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
      cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))
      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 


      !CF_volumique & eau condense
      !environnement
      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
      IntJ_CF=0.5*(1.-1.*erf(xenv2))
      if (deltasenv < 1.e-10) then
      qlenv(ind1,ind2)=IntJ
      cenv_vol(ind1,ind2)=IntJ_CF
      else
      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
      IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
      IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
      cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
      endif

      !thermique
      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
      IntJ_CF=0.5*(1.-1.*erf(xth2))
      if (deltasth < 1.e-10) then
      qlth(ind1,ind2)=IntJ
      cth_vol(ind1,ind2)=IntJ_CF
      else
      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
      IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
      IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
      cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
      endif

      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)

      ELSE IF (iflag_cloudth_vert == 5) THEN
         sigma1s=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
              /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
              +ratqs(ind1,ind2)*po(ind1) !Environment
      sigma2s=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.02)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)                   !Thermals
      !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
      !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) 
      xth=sth/(sqrt(2.)*sigma2s)
      xenv=senv/(sqrt(2.)*sigma1s)

      !Volumique
      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
      !print *,'jeanjean_CV=',ctot_vol(ind1,ind2)

      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth_vol(ind1,ind2))
      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv_vol(ind1,ind2))  
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)

      !Surfacique
      !Neggers
      !beta=0.0044
      !inverse_rho=1.+beta*dz(ind1,ind2)
      !print *,'jeanjean : beta=',beta
      !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
      !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
      !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)

      !Brooks
      a_Brooks=0.6694
      b_Brooks=0.1882
      A_Maj_Brooks=0.1635 !-- sans shear
      !A_Maj_Brooks=0.17   !-- ARM LES
      !A_Maj_Brooks=0.18   !-- RICO LES
      !A_Maj_Brooks=0.19   !-- BOMEX LES
      Dx_Brooks=200000.
      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
      !print *,'jeanjean_f=',f_Brooks

      cth(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cth_vol(ind1,ind2),1.)))- 1.))
      cenv(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(cenv_vol(ind1,ind2),1.)))- 1.))
      ctot(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
      !print *,'JJ_ctot_1',ctot(ind1,ind2)





      ENDIF ! of if (iflag_cloudth_vert<5)
      ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)

!      if (ctot(ind1,ind2).lt.1.e-10) then
      if (cenv(ind1,ind2)<1.e-10.or.cth(ind1,ind2)<1.e-10) then
      ctot(ind1,ind2)=0.
      ctot_vol(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)

      else
                
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
!    &             +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) 

      endif  

      else  ! gaussienne environnement seule
      

      zqenv(ind1)=po(ind1)
      Tbef=t(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
      

!      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)
      aenv=1./(1.+(alenv*Lv/cppd))
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
      sth=0.
      

      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
      sigma2s=0.

      sqrt2pi=sqrt(2.*pi)
      xenv=senv/(sqrt(2.)*sigma1s)
      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=ctot(ind1,ind2)
      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
      
      if (ctot(ind1,ind2)<1.e-3) then
      ctot(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)

      else   
                
!      ctot(ind1,ind2)=ctot(ind1,ind2) 
      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)

      endif  
 



      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
      ! Outputs used to check the PDFs
      cloudth_senv(ind1,ind2) = senv
      cloudth_sth(ind1,ind2) = sth
      cloudth_sigmaenv(ind1,ind2) = sigma1s
      cloudth_sigmath(ind1,ind2) = sigma2s

      enddo       ! from the loop on ngrid l.333
      return
!     end
END SUBROUTINE cloudth_vert_v3

       SUBROUTINE cloudth_v6(ngrid,klev,ind2,  &
             ztv,po,zqta,fraca, &
             qcloud,ctot_surf,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
             ratqs,zqs,T, &
             cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)

      use lmdz_cloudth_ini, only: iflag_cloudth_vert

      IMPLICIT NONE


      INCLUDE "YOMCST.h"
      INCLUDE "YOETHF.h"
      INCLUDE "FCTTRE.h"


        !Domain variables
      INTEGER ngrid !indice Max lat-lon
      INTEGER klev  !indice Max alt
      real, dimension(ngrid,klev), intent(out) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv
      INTEGER ind1  !indice in [1:ngrid]
      INTEGER ind2  !indice in [1:klev]
        !thermal plume fraction
      REAL fraca(ngrid,klev+1)   !thermal plumes fraction in the gridbox
        !temperatures
      REAL T(ngrid,klev)       !temperature
      REAL zpspsk(ngrid,klev)  !factor (p/p0)**kappa (used for potential variables)
      REAL ztv(ngrid,klev)     !potential temperature (voir thermcell_env.F90)
      REAL ztla(ngrid,klev)    !liquid temperature in the thermals (Tl_th)
      REAL zthl(ngrid,klev)    !liquid temperature in the environment (Tl_env)
        !pressure
      REAL paprs(ngrid,klev+1)   !pressure at the interface of levels
      REAL pplay(ngrid,klev)     !pressure at the middle of the level
        !humidity
      REAL ratqs(ngrid,klev)   !width of the total water subgrid-scale distribution
      REAL po(ngrid)           !total water (qt)
      REAL zqenv(ngrid)        !total water in the environment (qt_env)
      REAL zqta(ngrid,klev)    !total water in the thermals (qt_th)
      REAL zqsatth(ngrid,klev)   !water saturation level in the thermals (q_sat_th)
      REAL zqsatenv(ngrid,klev)  !water saturation level in the environment (q_sat_env)
      REAL qlth(ngrid,klev)    !condensed water in the thermals
      REAL qlenv(ngrid,klev)   !condensed water in the environment
      REAL qltot(ngrid,klev)   !condensed water in the gridbox
        !cloud fractions
      REAL cth_vol(ngrid,klev)   !cloud fraction by volume in the thermals
      REAL cenv_vol(ngrid,klev)  !cloud fraction by volume in the environment
      REAL ctot_vol(ngrid,klev)  !cloud fraction by volume in the gridbox
      REAL cth_surf(ngrid,klev)  !cloud fraction by surface in the thermals
      REAL cenv_surf(ngrid,klev) !cloud fraction by surface in the environment  
      REAL ctot_surf(ngrid,klev) !cloud fraction by surface in the gridbox
        !PDF of saturation deficit variables
      REAL rdd,cppd,Lv
      REAL Tbef,zdelta,qsatbef,zcor
      REAL alth,alenv,ath,aenv
      REAL sth,senv              !saturation deficits in the thermals and environment
      REAL sigma_env,sigma_th    !standard deviations of the biGaussian PDF
        !cloud fraction variables
      REAL xth,xenv
      REAL inverse_rho,beta                                  !Neggers et al. (2011) method
      REAL a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks !Brooks et al. (2005) method
        !Incloud total water variables
      REAL zqs(ngrid)    !q_sat
      REAL qcloud(ngrid) !eau totale dans le nuage
        !Some arithmetic variables
      REAL erf,pi,sqrt2,sqrt2pi
        !Depth of the layer
      REAL dz(ngrid,klev)    !epaisseur de la couche en metre
      REAL rhodz(ngrid,klev)
      REAL zrho(ngrid,klev)
      DO ind1 = 1, ngrid
        rhodz(ind1,ind2) = (paprs(ind1,ind2)-paprs(ind1,ind2+1))/rg ![kg/m2]
        zrho(ind1,ind2) = pplay(ind1,ind2)/T(ind1,ind2)/rd          ![kg/m3]
        dz(ind1,ind2) = rhodz(ind1,ind2)/zrho(ind1,ind2)            ![m]
      END DO

!------------------------------------------------------------------------------
! Initialization
!------------------------------------------------------------------------------
      qlth(:,:)=0.
      qlenv(:,:)=0.  
      qltot(:,:)=0.
      cth_vol(:,:)=0.
      cenv_vol(:,:)=0.
      ctot_vol(:,:)=0.
      cth_surf(:,:)=0.
      cenv_surf(:,:)=0.
      ctot_surf(:,:)=0.
      qcloud(:)=0.
      rdd=287.04
      cppd=1005.7
      pi=3.14159 
      Lv=2.5e6
      sqrt2=sqrt(2.)
      sqrt2pi=sqrt(2.*pi)


      DO ind1=1,ngrid
!-------------------------------------------------------------------------------
!Both thermal and environment in the gridbox
!-------------------------------------------------------------------------------
      IF ((ztv(ind1,1)>ztv(ind1,2)).and.(fraca(ind1,ind2)>1.e-10)) THEN
        !--------------------------------------------
        !calcul de qsat_env
        !--------------------------------------------
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
        !--------------------------------------------
        !calcul de s_env
        !--------------------------------------------
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam
        !--------------------------------------------
        !calcul de qsat_th
        !--------------------------------------------
      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatth(ind1,ind2)=qsatbef
        !--------------------------------------------
        !calcul de s_th  
        !--------------------------------------------
      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84 these Arnaud Jam
      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84 these Arnaud Jam
      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84 these Arnaud Jam
        !--------------------------------------------
        !calcul standard deviations bi-Gaussian PDF
        !--------------------------------------------
      sigma_th=(0.03218+0.000092655*dz(ind1,ind2))/((fraca(ind1,ind2)+0.01)**0.5)*(((sth-senv)**2)**0.5)+0.002*zqta(ind1,ind2)
      sigma_env=(0.71794+0.000498239*dz(ind1,ind2))*(fraca(ind1,ind2)**0.5) &
           /(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5) &
           +ratqs(ind1,ind2)*po(ind1)
      xth=sth/(sqrt2*sigma_th)
      xenv=senv/(sqrt2*sigma_env)
        !--------------------------------------------
        !Cloud fraction by volume CF_vol
        !--------------------------------------------
      cth_vol(ind1,ind2)=0.5*(1.+1.*erf(xth))
      cenv_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
      ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
        !-------------------------------------------- 
        !Condensed water qc
        !--------------------------------------------
      qlth(ind1,ind2)=sigma_th*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth_vol(ind1,ind2))
      qlenv(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv_vol(ind1,ind2))  
      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
        !--------------------------------------------
        !Cloud fraction by surface CF_surf
        !--------------------------------------------
        !Method Neggers et al. (2011) : ok for cumulus clouds only
      !beta=0.0044 (Jouhaud et al.2018)
      !inverse_rho=1.+beta*dz(ind1,ind2)
      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
        !Method Brooks et al. (2005) : ok for all types of clouds
      a_Brooks=0.6694
      b_Brooks=0.1882
      A_Maj_Brooks=0.1635 !-- sans dependence au cisaillement de vent
      Dx_Brooks=200000.   !-- si l'on considere des mailles de 200km de cote
      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
        !--------------------------------------------
        !Incloud Condensed water qcloud
        !--------------------------------------------
      if (ctot_surf(ind1,ind2) < 1.e-10) then
      ctot_vol(ind1,ind2)=0.
      ctot_surf(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)
      else
      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqs(ind1)
      endif



!-------------------------------------------------------------------------------
!Environment only in the gridbox
!-------------------------------------------------------------------------------
      ELSE
        !--------------------------------------------
        !calcul de qsat_env
        !--------------------------------------------
      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
      qsatbef=MIN(0.5,qsatbef)
      zcor=1./(1.-retv*qsatbef)
      qsatbef=qsatbef*zcor
      zqsatenv(ind1,ind2)=qsatbef
        !--------------------------------------------
        !calcul de s_env
        !--------------------------------------------
      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84 these Arnaud Jam
      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84 these Arnaud Jam
      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84 these Arnaud Jam 
        !--------------------------------------------
        !calcul standard deviations Gaussian PDF
        !--------------------------------------------
      zqenv(ind1)=po(ind1)
      sigma_env=ratqs(ind1,ind2)*zqenv(ind1)
      xenv=senv/(sqrt2*sigma_env)
        !--------------------------------------------
        !Cloud fraction by volume CF_vol
        !--------------------------------------------
      ctot_vol(ind1,ind2)=0.5*(1.+1.*erf(xenv))
        !--------------------------------------------
        !Condensed water qc
        !--------------------------------------------
      qltot(ind1,ind2)=sigma_env*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*ctot_vol(ind1,ind2))
        !--------------------------------------------
        !Cloud fraction by surface CF_surf
        !--------------------------------------------
        !Method Neggers et al. (2011) : ok for cumulus clouds only
      !beta=0.0044 (Jouhaud et al.2018)
      !inverse_rho=1.+beta*dz(ind1,ind2)
      !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
        !Method Brooks et al. (2005) : ok for all types of clouds
      a_Brooks=0.6694
      b_Brooks=0.1882
      A_Maj_Brooks=0.1635 !-- sans dependence au shear
      Dx_Brooks=200000.
      f_Brooks=A_Maj_Brooks*(dz(ind1,ind2)**(a_Brooks))*(Dx_Brooks**(-b_Brooks))
      ctot_surf(ind1,ind2)=1./(1.+exp(-1.*f_Brooks)*((1./max(1.e-15,min(ctot_vol(ind1,ind2),1.)))- 1.))
        !--------------------------------------------
        !Incloud Condensed water qcloud
        !--------------------------------------------
      if (ctot_surf(ind1,ind2) < 1.e-8) then
      ctot_vol(ind1,ind2)=0.
      ctot_surf(ind1,ind2)=0.
      qcloud(ind1)=zqsatenv(ind1,ind2)
      else
      qcloud(ind1)=qltot(ind1,ind2)/ctot_vol(ind1,ind2)+zqsatenv(ind1,ind2)
      endif


      END IF  ! From the separation (thermal/envrionnement) et (environnement only)

      ! Outputs used to check the PDFs
      cloudth_senv(ind1,ind2) = senv
      cloudth_sth(ind1,ind2) = sth
      cloudth_sigmaenv(ind1,ind2) = sigma_env
      cloudth_sigmath(ind1,ind2) = sigma_th

      END DO  ! From the loop on ngrid
      return

END SUBROUTINE cloudth_v6





!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE cloudth_mpc(klon,klev,ind2,mpc_bl_points,                        &
&           temp,qt,qt_th,frac_th,zpspsk,paprsup, paprsdn,play,thetal_th,   &
&           ratqs,qcloud,qincloud,icefrac,ctot,ctot_vol,                    &
&           cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Author : Etienne Vignon (LMDZ/CNRS)
! Date: April 2024
! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      use lmdz_cloudth_ini, only: iflag_cloudth_vert,iflag_ratqs
      use lmdz_cloudth_ini, only: C_mpc ,Ni,C_cap,Ei,d_top ,vert_alpha, vert_alpha_th ,sigma1s_factor,sigma1s_power,sigma2s_factor,sigma2s_power,cloudth_ratqsmin ,iflag_cloudth_vert_noratqs
      use lmdz_lscp_tools,  only: CALC_QSAT_ECMWF
      use lmdz_lscp_ini,    only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th

      IMPLICIT NONE


!------------------------------------------------------------------------------
! Declaration
!------------------------------------------------------------------------------

! INPUT/OUTPUT

      INTEGER, INTENT(IN)                         :: klon,klev,ind2
      

      REAL, DIMENSION(klon),      INTENT(IN)      ::  temp          ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt            ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  qt_th         ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  thetal_th     ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
      REAL, DIMENSION(klon),      INTENT(IN)      ::  frac_th       ! Fraction of the mesh covered by thermals [0-1] 
      REAL, DIMENSION(klon),      INTENT(IN)      ::  zpspsk        ! Exner potential
      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsup       ! Pressure at top layer interface [Pa]
      REAL, DIMENSION(klon),      INTENT(IN)      ::  paprsdn       ! Pressure at bottom layer interface [Pa]
      REAL, DIMENSION(klon),      INTENT(IN)      ::  play          ! Pressure of layers [Pa]
      REAL, DIMENSION(klon),      INTENT(IN)      ::  ratqs         ! Parameter that determines the width of the total water distrib.


      INTEGER, DIMENSION(klon,klev), INTENT(INOUT) :: mpc_bl_points  ! grid points with convective (thermals) mixed phase clouds
      
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot         ! Cloud fraction [0-1]
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  ctot_vol     ! Volume cloud fraction [0-1]
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qcloud       ! In cloud total water content [kg/kg] 
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  qincloud     ! In cloud condensed water content [kg/kg] 
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  icefrac      ! Fraction of ice in clouds [0-1]
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sth  ! mean saturation deficit in thermals
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_senv ! mean saturation deficit in environment
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmath ! std of saturation deficit in thermals
      REAL, DIMENSION(klon),      INTENT(OUT)      ::  cloudth_sigmaenv ! std of saturation deficit in environment


! LOCAL VARIABLES

      INTEGER itap,ind1,l,ig,iter,k 
      INTEGER iflag_topthermals, niter

      REAL qcth(klon)
      REAL qcenv(klon)
      REAL qctot(klon) 
      REAL cth(klon)
      REAL cenv(klon)   
      REAL cth_vol(klon)
      REAL cenv_vol(klon)
      REAL qt_env(klon), thetal_env(klon)
      REAL icefracenv(klon), icefracth(klon)
      REAL sqrtpi,sqrt2,sqrt2pi
      REAL alth,alenv,ath,aenv
      REAL sth,senv,sigma1s,sigma2s,sigma1s_fraca,sigma1s_ratqs
      REAL inverse_rho,beta,a_Brooks,b_Brooks,A_Maj_Brooks,Dx_Brooks,f_Brooks
      REAL xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
      REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth
      REAL zdelta,qsatbef,zcor
      REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
      REAL qlbef
      REAL dqsatenv(klon), dqsatth(klon)
      REAL erf
      REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
      REAL zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) 
      REAL rhodz(klon)
      REAL rho(klon)
      REAL dz(klon)
      REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)  
      REAL alenvl, aenvl
      REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
      REAL senvi, senvl, qbase, sbase, qliqth, qiceth
      REAL qmax, ttarget, stmp, cout, coutref, fraci
      REAL maxi, mini, pas

      ! Modifty the saturation deficit PDF in thermals 
      ! in the presence of ice crystals
      CHARACTER (len = 80) :: abort_message
      CHARACTER (len = 20) :: routname = 'cloudth_mpc'


!------------------------------------------------------------------------------
! Initialisation
!------------------------------------------------------------------------------


! Few initial checksS

      DO k = 1,klev
      DO ind1 = 1, klon
        rhodz(ind1) = (paprsdn(ind1)-paprsup(ind1))/rg !kg/m2
        zrho(ind1)  = play(ind1)/temp(ind1)/rd !kg/m3
        dz(ind1)    = rhodz(ind1)/zrho(ind1) !m : epaisseur de la couche en metre
      END DO
      END DO


      icefracth(:)=0.
      icefracenv(:)=0.
      sqrt2pi=sqrt(2.*rpi)
      sqrt2=sqrt(2.)
      sqrtpi=sqrt(rpi)
      icefrac(:)=0.

!-------------------------------------------------------------------------------
! Identify grid points with potential mixed-phase conditions
!-------------------------------------------------------------------------------


      DO ind1=1,klon
            IF ((temp(ind1) < RTT) .AND. (temp(ind1) > temp_nowater) &
            .AND. (ind2<=klev-2)  &
            .AND. (frac_th(ind1)>min_frac_th_cld)) THEN
                mpc_bl_points(ind1,ind2)=1
            ELSE
                mpc_bl_points(ind1,ind2)=0
            ENDIF
      ENDDO


!-------------------------------------------------------------------------------
! Thermal fraction calculation and standard deviation of the distribution
!-------------------------------------------------------------------------------  

! initialisations and calculation of temperature, humidity and saturation specific humidity

DO ind1=1,klon
  
   rhodz(ind1)   = (paprsdn(ind1)-paprsup(ind1))/rg  ! kg/m2
   rho(ind1)     = play(ind1)/temp(ind1)/rd          ! kg/m3
   dz(ind1)      = rhodz(ind1)/rho(ind1)             ! m : epaisseur de la couche en metre
   Tbefenv(ind1) = temp(ind1)
   thetal_env(ind1) = Tbefenv(ind1)/zpspsk(ind1)
   Tbefth(ind1)  = thetal_th(ind1)*zpspsk(ind1)
   rhoth(ind1)   = play(ind1)/Tbefth(ind1)/RD
   qt_env(ind1)  = (qt(ind1)-frac_th(ind1)*qt_th(ind1))/(1.-frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv 

ENDDO

! Calculation of saturation specific humidity
CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,1,.FALSE.,qslenv,dqsatenv)
CALL CALC_QSAT_ECMWF(klon,Tbefenv,qt_env,play,RTT,2,.FALSE.,qsienv,dqsatenv)
CALL CALC_QSAT_ECMWF(klon,Tbefth,qt_th,play,RTT,1,.FALSE.,qslth,dqsatth)


DO ind1=1,klon


    IF (frac_th(ind1)>min_frac_th_cld) THEN !Thermal and environnement

! unlike in the other cloudth routine, 
! We consider distributions of the saturation deficit WRT liquid 
! at positive AND negative celcius temperature
! subsequently, cloud fraction corresponds to the part of the pdf corresponding
! to superstauration with respect to liquid WHATEVER the temperature

! Environment:


        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)     
        aenv=1./(1.+(alenv*RLVTT/rcpd))                             
        senv=aenv*(qt_env(ind1)-qslenv(ind1))                           
     

! Thermals:


        alth=(0.622*RLVTT*qslth(ind1))/(rd*thetal_th(ind1)**2)       
        ath=1./(1.+(alth*RLVTT/rcpd))                                                         
        sth=ath*(qt_th(ind1)-qslth(ind1))                      


!       IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC


! Standard deviation of the distributions

           sigma1s_fraca = (sigma1s_factor**0.5)*(frac_th(ind1)**sigma1s_power) / &
                  (1-frac_th(ind1))*((sth-senv)**2)**0.5

           IF (cloudth_ratqsmin>0.) THEN
             sigma1s_ratqs = cloudth_ratqsmin*qt(ind1)
           ELSE
             sigma1s_ratqs = ratqs(ind1)*qt(ind1)
           ENDIF
 
           sigma1s = sigma1s_fraca + sigma1s_ratqs
           IF (iflag_ratqs==11) then
              sigma1s = ratqs(ind1)*qt(ind1)*aenv
           ENDIF
           sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((frac_th(ind1)+0.02)**sigma2s_power))+0.002*qt_th(ind1)


           deltasenv=aenv*vert_alpha*sigma1s
           deltasth=ath*vert_alpha_th*sigma2s

           xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
           xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
           exp_xenv1 = exp(-1.*xenv1**2)
           exp_xenv2 = exp(-1.*xenv2**2)
           xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
           xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
           exp_xth1 = exp(-1.*xth1**2)
           exp_xth2 = exp(-1.*xth2**2)
      
!surface CF

           cth(ind1)=0.5*(1.-1.*erf(xth1))
           cenv(ind1)=0.5*(1.-1.*erf(xenv1))
           ctot(ind1)=frac_th(ind1)*cth(ind1)+(1.-1.*frac_th(ind1))*cenv(ind1) 


!volume CF, condensed water and ice fraction

            !environnement

            IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
            IntJ_CF=0.5*(1.-1.*erf(xenv2))

            IF (deltasenv < 1.e-10) THEN
              qcenv(ind1)=IntJ
              cenv_vol(ind1)=IntJ_CF
            ELSE
              IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
              IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
              qcenv(ind1)=IntJ+IntI1+IntI2+IntI3
              cenv_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
              IF (Tbefenv(ind1) < temp_nowater) THEN
                  ! freeze all droplets in cirrus temperature regime
                  icefracenv(ind1)=1. 
              ENDIF
            ENDIF
             


            !thermals

            IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
            IntJ_CF=0.5*(1.-1.*erf(xth2))
     
            IF (deltasth < 1.e-10) THEN
              qcth(ind1)=IntJ
              cth_vol(ind1)=IntJ_CF
            ELSE
              IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
              IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
              IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
              qcth(ind1)=IntJ+IntI1+IntI2+IntI3
              cth_vol(ind1)=IntJ_CF+IntI1_CF+IntI3_CF
              IF (Tbefth(ind1) < temp_nowater) THEN
                  ! freeze all droplets in cirrus temperature regime
                  icefracth(ind1)=1.
              ENDIF

            ENDIF

            qctot(ind1)=frac_th(ind1)*qcth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)
            ctot_vol(ind1)=frac_th(ind1)*cth_vol(ind1)+(1.-1.*frac_th(ind1))*cenv_vol(ind1)

            IF (cenv(ind1)<min_neb_th.and.cth(ind1)<min_neb_th) THEN
                ctot(ind1)=0.
                ctot_vol(ind1)=0.
                qcloud(ind1)=qslenv(ind1)
                qincloud(ind1)=0.
                icefrac(ind1)=0.
            ELSE               
                qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
                qincloud(ind1)=qctot(ind1)/ctot(ind1)
                IF (qctot(ind1) > 0) THEN
                  icefrac(ind1)=(frac_th(ind1)*qcth(ind1)*icefracth(ind1)+(1.-1.*frac_th(ind1))*qcenv(ind1)*icefracenv(ind1))/qctot(ind1)
                  icefrac(ind1)=max(min(1.,icefrac(ind1)), 0.)
                ENDIF
            ENDIF 


!        ELSE ! mpc_bl_points>0

    ELSE  ! gaussian for environment only

      
        alenv=(0.622*RLVTT*qslenv(ind1))/(rd*thetal_env(ind1)**2)
        aenv=1./(1.+(alenv*RLVTT/rcpd))
        senv=aenv*(qt_env(ind1)-qslenv(ind1))
        sth=0.
        sigma1s=ratqs(ind1)*qt_env(ind1)
        sigma2s=0.

        xenv=senv/(sqrt(2.)*sigma1s)
        cenv(ind1)=0.5*(1.-1.*erf(xenv))
        ctot(ind1)=0.5*(1.+1.*erf(xenv))
        ctot_vol(ind1)=ctot(ind1)
        qctot(ind1)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1))


        IF (ctot(ind1)<min_neb_th) THEN
          ctot(ind1)=0.
          qcloud(ind1)=qslenv(ind1)
          qincloud(ind1)=0.
        ELSE 
          qcloud(ind1)=qctot(ind1)/ctot(ind1)+qslenv(ind1)
          qincloud(ind1)=MAX(qctot(ind1)/ctot(ind1),0.)
        ENDIF
 

    ENDIF       ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492

    ! Outputs used to check the PDFs
    cloudth_senv(ind1) = senv
    cloudth_sth(ind1) = sth
    cloudth_sigmaenv(ind1) = sigma1s
    cloudth_sigmath(ind1) = sigma2s


    ENDDO       !loop on klon


RETURN


END SUBROUTINE cloudth_mpc



!        ELSE ! mpc_bl_points>0

!            ! Treat boundary layer mixed phase clouds

!            ! thermals
!            !=========

!           ! ice phase
!           !...........

!            qiceth=0.
!            deltazlev_mpc=dz(ind1,:)
!            temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
!            pres_mpc=pplay(ind1,:)
!            fraca_mpc=fraca(ind1,:)
!            snowf_mpc=snowflux(ind1,:)
!            iflag_topthermals=0
!            IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0))  THEN
!                iflag_topthermals = 1
!            ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & 
!                    .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
!                iflag_topthermals = 2
!            ELSE 
!                iflag_topthermals = 0
!            ENDIF

!            CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
!                                   qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))

!            ! qmax calculation
!            sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
!            deltasth=athl*vert_alpha_th*sigma2s
!            xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
!            xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
!            exp_xth1 = exp(-1.*xth1**2)
!            exp_xth2 = exp(-1.*xth2**2)
!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
!            IntJ_CF=0.5*(1.-1.*erf(xth2))
!            IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
!            IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
!            IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
!            IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
!            IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
!            qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)


!            ! Liquid phase
!            !................
!            ! We account for the effect of ice crystals in thermals on sthl
!            ! and on the width of the distribution

!            sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2))  &
!                + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1))  

!            sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) &
!                 /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
!                 +0.002*zqta(ind1,ind2)
!            deltasthc=athl*vert_alpha_th*sigma2sc


!            xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
!            xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)           
!            exp_xth1 = exp(-1.*xth1**2)
!            exp_xth2 = exp(-1.*xth2**2)
!            IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
!            IntJ_CF=0.5*(1.-1.*erf(xth2))
!            IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1))
!            IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
!            IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2)
!            IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc)
!            IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc)
!            qliqth=IntJ+IntI1+IntI2+IntI3

!            qlth(ind1,ind2)=MAX(0.,qliqth)

!            ! Condensed water

!            qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)


!            ! consistency with subgrid distribution

!             IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
!                 fraci=qith(ind1,ind2)/qcth(ind1,ind2)
!                 qcth(ind1,ind2)=qmax
!                 qith(ind1,ind2)=fraci*qmax
!                 qlth(ind1,ind2)=(1.-fraci)*qmax
!             ENDIF

!            ! Cloud Fraction
!            !...............
!            ! calculation of qbase which is the value of the water vapor within mixed phase clouds
!            ! such that the total water in cloud = qbase+qliqth+qiceth
!            ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
!            ! sbase and qbase calculation (note that sbase is wrt liq so negative)
!            ! look for an approximate solution with iteration

!            ttarget=qcth(ind1,ind2)
!            mini= athl*(qsith(ind1,ind2)-qslth(ind1))
!            maxi= 0. !athl*(3.*sqrt(sigma2s))
!            niter=20
!            pas=(maxi-mini)/niter
!            stmp=mini
!            sbase=stmp
!            coutref=1.E6
!            DO iter=1,niter
!                cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) &
!                     + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
!               IF (cout .LT. coutref) THEN
!                     sbase=stmp
!                     coutref=cout
!                ELSE
!                     stmp=stmp+pas
!                ENDIF
!            ENDDO
!            qbase=MAX(0., sbase/athl+qslth(ind1))

!            ! surface cloud fraction in thermals
!            cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
!            cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)


!            !volume cloud fraction in thermals
!            !to be checked
!            xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
!            xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)           
!            exp_xth1 = exp(-1.*xth1**2)
!            exp_xth2 = exp(-1.*xth2**2)

!            IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
!            IntJ_CF=0.5*(1.-1.*erf(xth2))

!            IF (deltasth .LT. 1.e-10) THEN
!              cth_vol(ind1,ind2)=IntJ_CF
!            ELSE
!              IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
!              IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
!              IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
!              IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth)
!              IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth)
!              cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
!            ENDIF
!              cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)



!            ! Environment
!            !=============
!            ! In the environment/downdrafts, only liquid clouds
!            ! See Shupe et al. 2008, JAS

!            ! standard deviation of the distribution in the environment
!            sth=sthl
!            senv=senvl
!            sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / &
!                          &                (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5
!            ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
!            ! in the environement

!            sigma1s_ratqs=1E-10
!            IF (cloudth_ratqsmin>0.) THEN
!                sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
!            ENDIF

!            sigma1s = sigma1s_fraca + sigma1s_ratqs
!            IF (iflag_ratqs.eq.11) then
!               sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
!            ENDIF
!            IF (iflag_ratqs.eq.11) then
!              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
!            ENDIF
!            deltasenv=aenvl*vert_alpha*sigma1s
!            xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
!            xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
!            exp_xenv1 = exp(-1.*xenv1**2)
!            exp_xenv2 = exp(-1.*xenv2**2)

!            !surface CF
!            cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1))

!            !volume CF and condensed water
!            IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
!            IntJ_CF=0.5*(1.-1.*erf(xenv2))

!            IF (deltasenv .LT. 1.e-10) THEN
!              qcenv(ind1,ind2)=IntJ
!              cenv_vol(ind1,ind2)=IntJ_CF
!            ELSE
!              IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
!              IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
!              IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
!              IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv)
!              IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv)
!              qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
!              cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
!            ENDIF

!            qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
!            cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)



!            ! Thermals + environment
!            ! ======================= 
!            ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) 
!            qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)
!            ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2)
!            IF (qcth(ind1,ind2) .GT. 0) THEN
!               icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
!                    /(fraca(ind1,ind2)*qcth(ind1,ind2) &
!                    +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2))
!                icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
!            ELSE
!                icefrac(ind1,ind2)=0.
!            ENDIF

!            IF (cenv(ind1,ind2).LT.1.e-10.or.cth(ind1,ind2).LT.1.e-10) THEN
!                ctot(ind1,ind2)=0.
!                ctot_vol(ind1,ind2)=0.
!                qincloud(ind1)=0.
!                qcloud(ind1)=zqsatenv(ind1)
!            ELSE               
!                qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
!                qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
!                            +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
!            ENDIF 

!        ENDIF ! mpc_bl_points



!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp,pres,qth,qsith,qlth,deltazlev,vith,fraca,qith)

! parameterization of ice for boundary
! layer mixed-phase clouds assuming a stationary system

! Note that vapor deposition on ice crystals and riming of liquid droplets
! depend on the ice number concentration Ni
! One could assume that Ni depends on qi, e.g.,  Ni=beta*(qi*rho)**xi
! and use values from Hong et al. 2004, MWR for instance
! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962 
! One could also think of a more complex expression of Ni;
! function of qi, T, the concentration in aerosols or INP ..
! Here we prefer fixing Ni to a tuning parameter
! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
! in Mioche et al. 2017


! References:
!------------
! This parameterization is thoroughly described in Vignon et al. 

! More specifically
! for the Water vapor deposition process:

! Rotstayn, L. D., 1997: A physically based scheme for the treat-
! ment of stratiform cloudfs and precipitation in large-scale
! models. I: Description and evaluation of the microphysical
! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.

!  Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
!  stratiform cloud microphysics scheme in the NCAR Com-
!  munity Atmosphere Model (CAM3). Part I: Description and
!  numerical tests. J. Climate, 21, 3642–3659

! for the Riming process:

! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
! scale structure and organization of clouds and precipitation in
! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206

! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
! forecasts of winter precipitation using an improved bulk
! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
! forecasts of winter precipitation using an improved bulk
! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542

! For the formation of clouds by thermals:

! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
! the Atmospheric Sciences, 65, 407–425.

! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3



! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
!=============================================================================

    use phys_state_var_mod, ONLY: fm_therm, detr_therm, entr_therm

    IMPLICIT none

    INCLUDE "YOMCST.h"

    INTEGER, INTENT(IN) :: ind1,ind2, klev           ! horizontal and vertical indices and dimensions
    INTEGER, INTENT(IN) :: iflag_topthermals         ! uppermost layer of thermals ?
    REAL, INTENT(IN)    :: Ni                        ! ice number concentration [m-3]
    REAL, INTENT(IN)    :: Ei                        ! ice-droplet collision efficiency
    REAL, INTENT(IN)    :: C_cap                     ! ice crystal capacitance
    REAL, INTENT(IN)    :: d_top                     ! cloud-top ice crystal mixing parameter
    REAL,  DIMENSION(klev), INTENT(IN) :: temp       ! temperature [K] within thermals
    REAL,  DIMENSION(klev), INTENT(IN) :: pres       ! pressure [Pa]
    REAL,  DIMENSION(klev), INTENT(IN) :: qth        ! mean specific water content in thermals [kg/kg]
    REAL,  DIMENSION(klev), INTENT(IN) :: qsith        ! saturation specific humidity wrt ice in thermals [kg/kg]
    REAL,  DIMENSION(klev), INTENT(IN) :: qlth       ! condensed liquid water in thermals, approximated value [kg/kg]
    REAL,  DIMENSION(klev), INTENT(IN) :: deltazlev  ! layer thickness [m]
    REAL,  DIMENSION(klev), INTENT(IN) :: vith       ! ice crystal fall velocity [m/s]
    REAL,  DIMENSION(klev+1), INTENT(IN) :: fraca      ! fraction of the mesh covered by thermals
    REAL,  DIMENSION(klev), INTENT(INOUT) :: qith       ! condensed ice water , thermals [kg/kg]


    INTEGER ind2p1,ind2p2
    REAL rho(klev)
    REAL unsurtaudet, unsurtaustardep, unsurtaurim
    REAL qi, AA, BB, Ka, Dv, rhoi
    REAL p0, t0, fp1, fp2
    REAL alpha, flux_term
    REAL det_term, precip_term, rim_term, dep_term


    ind2p1=ind2+1
    ind2p2=ind2+2

    rho=pres/temp/RD  ! air density kg/m3

    Ka=2.4e-2      ! thermal conductivity of the air, SI
    p0=101325.0    ! ref pressure
    T0=273.15      ! ref temp
    rhoi=500.0     ! cloud ice density following Reisner et al. 1998
    alpha=700.     ! fallvelocity param


    IF (iflag_topthermals > 0) THEN ! uppermost thermals levels

    Dv=0.0001*0.211*(p0/pres(ind2))*((temp(ind2)/T0)**1.94) ! water vapor diffusivity in air, SI

    ! Detrainment term:

    unsurtaudet=detr_therm(ind1,ind2)/rho(ind2)/deltazlev(ind2)
 

    ! Deposition term
    AA=RLSTT/Ka/temp(ind2)*(RLSTT/RV/temp(ind2)-1.)
    BB=1./(rho(ind2)*Dv*qsith(ind2))
    unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2)-qsith(ind2))/qsith(ind2) &
                    *4.*RPI/(AA+BB)*(6.*rho(ind2)/rhoi/RPI/Gamma(4.))**(0.33)

    ! Riming term  neglected at this level
    !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4)

    qi=fraca(ind2)*unsurtaustardep/MAX((d_top*unsurtaudet),1E-12)
    qi=MAX(qi,0.)**(3./2.)

    ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2

    Dv=0.0001*0.211*(p0/pres(ind2p1))*((temp(ind2p1)/T0)**1.94) ! water vapor diffusivity in air, SI

    ! Detrainment term:

    unsurtaudet=detr_therm(ind1,ind2p1)/rho(ind2p1)/deltazlev(ind2p1)
    det_term=-unsurtaudet*qith(ind2p1)*rho(ind2p1)
    
    
    ! Deposition term

    AA=RLSTT/Ka/temp(ind2p1)*(RLSTT/RV/temp(ind2p1)-1.)
    BB=1./(rho(ind2p1)*Dv*qsith(ind2p1))
    unsurtaustardep=C_cap*(Ni**0.66)*(qth(ind2p1)-qsith(ind2p1)) &
         /qsith(ind2p1)*4.*RPI/(AA+BB) &
         *(6.*rho(ind2p1)/rhoi/RPI/Gamma(4.))**(0.33)
    dep_term=rho(ind2p1)*fraca(ind2p1)*(qith(ind2p1)**0.33)*unsurtaustardep
 
    ! Riming term

    unsurtaurim=rho(ind2p1)*alpha*3./rhoi/2.*Ei*qlth(ind2p1)*((p0/pres(ind2p1))**0.4)
    rim_term=rho(ind2p1)*fraca(ind2p1)*qith(ind2p1)*unsurtaurim

    ! Precip term

    ! We assume that there is no solid precipitation outside thermals 
    ! so the precipitation flux within thermals is equal to the precipitation flux
    ! at mesh-scale divided by  thermals fraction
    
    fp2=0.
    fp1=0.
    IF (fraca(ind2p1) > 0.) THEN
    fp2=-qith(ind2p2)*rho(ind2p2)*vith(ind2p2)*fraca(ind2p2)! flux defined positive upward
    fp1=-qith(ind2p1)*rho(ind2p1)*vith(ind2p1)*fraca(ind2p1)
    ENDIF

    precip_term=-1./deltazlev(ind2p1)*(fp2-fp1)

    ! Calculation in a top-to-bottom loop

    IF (fm_therm(ind1,ind2p1) > 0.) THEN
          qi= 1./fm_therm(ind1,ind2p1)* &
              (deltazlev(ind2p1)*(-rim_term-dep_term-det_term-precip_term) + &
              fm_therm(ind1,ind2p2)*(qith(ind2p1)))
    END IF

    ENDIF ! top thermals

    qith(ind2)=MAX(0.,qi)

    RETURN

END SUBROUTINE ICE_MPC_BL_CLOUDS

END MODULE lmdz_cloudth
