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


      USE IOIPSL, ONLY : getin
      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"
#include "thermcell.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 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 alpha,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

      REAL, SAVE :: iflag_cloudth_vert, iflag_cloudth_vert_omp=0


      LOGICAL, SAVE :: first=.true.





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Astuce pour gérer deux versions de cloudth en attendant
! de converger sur une version nouvelle
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IF (first) THEN
     !$OMP MASTER
     CALL getin('iflag_cloudth_vert',iflag_cloudth_vert_omp)
     !$OMP END MASTER
     !$OMP BARRIER
     iflag_cloudth_vert=iflag_cloudth_vert_omp
      first=.false.
      ENDIF
       IF (iflag_cloudth_vert==1) THEN
       CALL cloudth_vert(ngrid,klev,ind2,  &
     &           ztv,po,zqta,fraca, & 
     &           qcloud,ctot,zpspsk,paprs,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).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.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)   
!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*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)
!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
     

!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      if (ctot(ind1,ind2).lt.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).lt.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



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


      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"
#include "thermcell.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 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 alpha,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





!------------------------------------------------------------------------------
! 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).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.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)   
!      ctot(ind1,ind2)=alpha*cth(ind1,ind2)+(1.-1.*alpha)*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)
!      qltot(ind1,ind2)=alpha*qlth(ind1,ind2)+(1.-1.*alpha)*qlenv(ind1,ind2)
     

!      print*,senv,sth,sigma1s,sigma2s,fraca(ind1,ind2),'senv et sth et sig1 et sig2 et alpha'


!-------------------------------------------------------------------------------
!  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)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.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).lt.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