Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90	(revision 1392)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90	(revision 1399)
@@ -8,5 +8,5 @@
      &      ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth,  &
      &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
-     &       zmax0,f0,zw2,fraca)
+     &       zmax0,f0,zw2,fraca,ztv,zpspsk,ztla,zthl)
 
       USE dimphy
@@ -49,4 +49,8 @@
       real zqla(klon,klev)
       real zqta(klon,klev)
+      real ztv(klon,klev)
+      real zpspsk(klon,klev)
+      real ztla(klon,klev)
+      real zthl(klon,klev)
       real wmax_sec(klon)
       real zmax_sec(klon)
@@ -214,5 +218,6 @@
      &      ,r_aspect_thermals,l_mix_thermals &
      &      ,tau_thermals,iflag_thermals_ed,Ale,Alp,lalim_conv,wght_th &
-     &      ,zmax0,f0,zw2,fraca)
+     &      ,zmax0,f0,zw2,fraca,ztv,zpspsk &
+     &      ,ztla,zthl)
            if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
          else
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90	(revision 1399)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cloudth.F90	(revision 1399)
@@ -0,0 +1,263 @@
+       SUBROUTINE cloudth(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 "iniprint.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
+
+
+
+
+
+!      print*,ngrid,klev,ind1,ind2,ztv(ind1,ind2),po(ind1),zqta(ind1,ind2), &
+!     &       fraca(ind1,ind2),zpspsk(ind1,ind2),paprs(ind1,ind2),ztla(ind1,ind2),zthl(ind1,ind2), &
+!     &       'verif'
+
+
+!      LOGICAL active(ngrid)   
+      
+!-----------------------------------------------------------------------------------------------------------------
+! 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
+!-----------------------------------------------------------------------------------------------------------------
+!      alpha=0.5*(fraca(ind1,ind2)+fraca(ind1,ind2-1))
+!      sigma1s=(Max(1.2*fraca(ind1,ind2)*(sth-senv)**2,(senv/100)**2))**0.5
+!      sigma2s=((0.021/(fraca(ind1,ind2)+0.0)**0.5)*(sth-senv)**2+(sth/100)**2)**0.5 
+!      sigma1s=(1.5**0.5)*(fraca(ind1,ind2)**0.65)*(sth-senv)+0.000045
+!      sigma2s=0.1265*(sth-senv)/(fraca(ind1,ind2)+0.0)**0.35+0.000045     
+
+      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.00003
+      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003  
+
+!      sigma1s=(1.5**0.5)*(alpha**0.65)*(sth-senv)+0.000045
+!      sigma2s=0.1265*(sth-senv)/alpha**0.35+0.000045
+      
+!      sigma1s=(2.8**0.5)*(0.1**0.7)*(sth-senv)+0.00002
+!      sigma2s=((0.126/(0.1)**0.3)*(sth-senv))+0.00002
+
+
+ 
+!-----------------------------------------------------------------------------------------------------------------
+! 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)) 
+      
+!      if (zqenv(ind1).gt.0.005) then
+      sigma1s=0.005*zqenv(ind1)
+!      else 
+!      sigma1s=0.0005
+!      endif
+
+!      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
+
+
+
+!      sigma1s=0.00003
+      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
+
+
+
+
+
+                                                                            
+
+
+
+
+
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F	(revision 1392)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/fisrtilp.F	(revision 1399)
@@ -7,5 +7,6 @@
      s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
      s                   frac_impa, frac_nucl,
-     s                   prfl, psfl, rhcl)
+     s                   prfl, psfl, rhcl, zqta, fraca,
+     s                   ztv, zpspsk, ztla, zthl, iflag_cldcon)
 
 c
@@ -41,5 +42,12 @@
       REAL snow(klon) ! neige (mm/s)
       REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
-      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
+      REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) 
+      REAL ztv(klon,klev)
+      REAL zqta(klon,klev),fraca(klon,klev) 
+      REAL sigma1(klon,klev),sigma2(klon,klev)
+      REAL qltot(klon,klev),ctot(klon,klev)
+      REAL zpspsk(klon,klev),ztla(klon,klev)
+      REAL zthl(klon,klev)
+
 cAA
 c Coeffients de fraction lessivee : pour OFF-LINE
@@ -63,5 +71,6 @@
 
       INTEGER ninter ! sous-intervals pour la precipitation
-      INTEGER ncoreczq
+      INTEGER ncoreczq  
+      INTEGER iflag_cldcon
       PARAMETER (ninter=5)
       LOGICAL evap_prec ! evaporation de la pluie
@@ -72,5 +81,6 @@
       real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
       real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
-      real erf
+      real erf   
+      REAL qcloud(klon)
 c
       LOGICAL cpartiel ! condensation partielle
@@ -82,5 +92,5 @@
 c
       INTEGER i, k, n, kk
-      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
+      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
       REAL zrfl(klon), zrfln(klon), zqev, zqevt
       REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
@@ -130,4 +140,5 @@
       zdelq=0.0
       
+      print*,'CLOUDTH4 A. JAM'
       IF (appel1er) THEN
 c
@@ -322,4 +333,5 @@
 c de l'eau condensee:
 c
+
       IF (cpartiel) THEN
 
@@ -351,6 +363,21 @@
                 zq(i)=1.e-15
               endif
-           enddo
-           do i=1,klon
+           enddo 
+
+              if (iflag_cldcon.eq.5) then
+
+                 call cloudth(klon,klev,k,ztv,
+     .           zq,zqta,fraca,
+     .           qcloud,ctot,zpspsk,paprs,ztla,zthl,
+     .           ratqs,zqs,t)
+
+                 do i=1,klon
+                 rneb(i,k)=ctot(i,k)
+                 zqn(i)=qcloud(i)
+                 enddo
+
+              else
+
+            do i=1,klon
             zpdf_sig(i)=ratqs(i,k)*zq(i)
             zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
@@ -372,5 +399,7 @@
             endif
             
-           enddo
+           enddo 
+
+         endif ! iflag_cldcon
 
         endif ! iflag_pdf
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1392)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1399)
@@ -669,5 +669,9 @@
 
       REAL zw2(klon,klev+1)
-      REAL fraca(klon,klev+1)
+      REAL fraca(klon,klev+1)        
+      REAL ztv(klon,klev) 
+      REAL zpspsk(klon,klev)
+      REAL ztla(klon,klev) 
+      REAL zthl(klon,klev)
 
 c Variables locales pour la couche limite (al1):
@@ -2440,5 +2444,6 @@
      s      ,ratqsdiff,zqsatth
 con rajoute ale et alp, et les caracteristiques de la couche alim
-     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca)
+     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca
+     s      ,ztv,zpspsk,ztla,zthl)
 
 ! ----------------------------------------------------------------------
@@ -2702,5 +2707,6 @@
      .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
      .           frac_impa, frac_nucl,
-     .           prfl, psfl, rhcl)
+     .           prfl, psfl, rhcl, 
+     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
 
       WHERE (rain_lsc < 0) rain_lsc = 0.
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90	(revision 1392)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90	(revision 1399)
@@ -10,5 +10,6 @@
      &                  ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed &
      &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
-     &                  ,zmax0, f0,zw2,fraca)
+     &                  ,zmax0, f0,zw2,fraca,ztv &
+     &                  ,zpspsk,ztla,zthl)
 
       USE dimphy
@@ -389,5 +390,5 @@
 !      print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
       if (iflag_thermals_ed<=9) then
-!         print*,'THERM NOVUELLE/NOUVELLE/ANCIENNE'
+!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'
          CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
      &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
@@ -396,5 +397,5 @@
      &    ,lev_out,lunout1,igout)
 
-      elseif (iflag_thermals_ed<=19) then
+      elseif (iflag_thermals_ed>9) then
 !        print*,'THERM RIO et al 2010, version d Arnaud'
          CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90	(revision 1392)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90	(revision 1399)
@@ -64,5 +64,4 @@
       REAL zqsatth(ngrid,klev)
       REAL zta_est(ngrid,klev)
-      REAL ztemp(ngrid),zqsat(ngrid)
       REAL zdw2
       REAL zw2modif
@@ -76,6 +75,8 @@
       INTEGER ig,l,k
 
-      real zdz,zfact,zbuoy,zalpha
+      real zdz,zfact,zbuoy,zalpha,zdrag
       real zcor,zdelta,zcvm5,qlbef
+      real Tbef,qsatbef
+      real dqsat_dT,DT,num,denom
       REAL REPS,RLvCp,DDT0
       PARAMETER (DDT0=.01)
@@ -84,4 +85,5 @@
       REAL fact_gamma,fact_epsilon,fact_gamma2
       REAL c2(ngrid,klev)
+      REAL a1,m
 
       REAL zw2fact,expa
@@ -90,32 +92,16 @@
       RLvCp = RLVTT/RCPD
      
-      if (iflag_thermals_ed==0) then
-         fact_gamma=1.
-         fact_epsilon=1.
-      else if (iflag_thermals_ed==1)  then
-! Valeurs au moment de la premiere soumission des papiers
-         fact_gamma=1.
+      
          fact_epsilon=0.002
-         fact_gamma2=0.6
-! Suggestions des Fleurs, Septembre 2009
-         fact_epsilon=0.015
-!test cr
-!         fact_epsilon=0.002
+         a1=2./3.
          fact_gamma=0.9
-         fact_gamma2=0.7
-
-      else if (iflag_thermals_ed==2)  then
-         fact_gamma=1.
-         fact_epsilon=2.
-      else if (iflag_thermals_ed==3)  then
-         fact_gamma=3./4.
-         fact_epsilon=3.
-      endif
-
-!      write(lunout,*)'THERM 31H '
-
-      print*,'THERMCELL_PLUME OPTIMISE V0 '
+         zfact=fact_gamma/(1+fact_gamma)
+         fact_gamma2=zfact
+         expa=0.
+      
+      print*,'THERM 31H '
+
 ! Initialisations des variables reeles
-if (1==0) then
+if (1==1) then
       ztva(:,:)=ztv(:,:)
       ztva_est(:,:)=ztva(:,:)
@@ -168,5 +154,5 @@
                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
      &                       *sqrt(zlev(ig,l+1)) 
-               lalim(ig)=l+1
+               lalim(:)=l+1
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
             endif
@@ -187,5 +173,5 @@
 ! On decide dans cette version que le thermique n'est actif que si la premiere
 ! couche est instable.
-! Pourrait etre change si on veut que le thermiques puisse se déclencher
+! Pourrait etre change si on veut que le thermiques puisse se dÃ©clencher
 ! dans une couche l>1
 !------------------------------------------------------------------------------
@@ -239,17 +225,13 @@
 ! sans tenir compte du detrainement et de l'entrainement dans cette
 ! couche
-! C'est a dire qu'on suppose 
-! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
 ! avant) a l'alimentation pour avoir un calcul plus propre
 !---------------------------------------------------------------------------
 
-   ztemp(:)=zpspsk(:,l)*ztla(:,l-1)
-   call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat)
-
+     call thermcell_condens(ngrid,active, &
+&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
 
     do ig=1,ngrid
         if(active(ig)) then
-        zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig))
         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
         zta_est(ig,l)=ztva_est(ig,l)
@@ -258,4 +240,6 @@
      &      -zqla_est(ig,l))-zqla_est(ig,l))
 
+         if (1.eq.0) then 
+!calcul de w_est sans prendre en compte le drag 
              w_est(ig,l+1)=zw2(ig,l)*  &
      &                   ((f_star(ig,l))**2)  &
@@ -263,4 +247,16 @@
      &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
      &                   *(zlev(ig,l+1)-zlev(ig,l))
+         else
+
+           zdz=zlev(ig,l+1)-zlev(ig,l)
+           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
+           zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
+           zdrag=fact_epsilon/(zalpha**expa)
+           zw2fact=zbuoy/zdrag*a1
+           w_est(ig,l+1)=(w_est(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
+      &    +zw2fact
+
+         endif
+
              if (w_est(ig,l+1).lt.0.) then
                 w_est(ig,l+1)=zw2(ig,l)
@@ -277,65 +273,34 @@
           zdz=zlev(ig,l+1)-zlev(ig,l)
           zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
-          zfact=fact_gamma/(1.+fact_gamma)
  
 ! estimation de la fraction couverte par les thermiques
           zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
 
-!calcul de la soumission papier
-          if (1.eq.1) then
-             fact_epsilon=0.0007
-!             zfact=0.9/(1.+0.9)
-             zfact=0.3
-             fact_gamma=0.7
-             fact_gamma2=0.6
-             expa=0.25
+!calcul de la soumission papier 
 ! Calcul  du taux d'entrainement entr_star (epsilon)
            entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
-     &     zbuoy/w_est(ig,l+1) )&
-!- fact_epsilon/zalpha**0.25  ) &
-     &     +0.000 )
-
-!           entr_star(ig,l)=f_star(ig,l)*zdz * (  1./3 * MAX(0.,  &     
-!     &     zbuoy/w_est(ig,l+1) - 1./zalpha**0.25  ) &
-!     &     +0.000 )
-! Calcul du taux de detrainement detr_star (delta)
-!           if (zqla_est(ig,l).lt.1.e-10) then
+     &     a1*zbuoy/w_est(ig,l+1) &
+     &     - fact_epsilon/zalpha**expa  ) &
+     &     +0. )
+           
+!calcul du taux de detrainment (delta)
 !           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
-!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
-!     &     +0.0006 )
-!           else
+!     &      MAX(1.e-3, &
+!     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
+!     &      +0.01*(max(zqta(ig,l-1)-po(ig,l),0.)/(po(ig,l))/(w_est(ig,l+1)))**0.5    &    
+!     &     +0. ))
+
+          m=0.5
+
+          detr_star(ig,l)=1.*f_star(ig,l)*zdz *                    &
+    &     MAX(5.e-4,-fact_gamma2*a1*(1./w_est(ig,l+1))*((1.-(1.-m)/(1.+70*zqta(ig,l-1)))*zbuoy        &
+    &     -40*(1.-m)*(max(zqta(ig,l-1)-po(ig,l),0.))/(1.+70*zqta(ig,l-1)) )   )
+
 !           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
-!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
-!     &     +0.002 )
-!           endif
-           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
-     &     fact_gamma2 * MAX(0., &
-!fact_epsilon/zalpha**0.25
-     &      -zbuoy/w_est(ig,l+1) )       &
-!     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
-!test
-     &     +  0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &    
-     &     +0.0000 )
-          else
-
-! Calcul  du taux d'entrainement entr_star (epsilon)
-           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
-     &     zbuoy/w_est(ig,l+1) - fact_epsilon  ) &
-     &     +0.0000 )
-
-! Calcul du taux de detrainement detr_star (delta)
-           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
-     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
-     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
-     &     +0.0000 )
-
-          endif
-!endif choix du calcul de E* et D*
-
-!cr test
-!           entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)     
-
-! Prise en compte de la fraction
-!      detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
+!     &      MAX(0.0, &
+!     &      -fact_gamma2*a1*zbuoy/w_est(ig,l+1)        &
+!     &      +20*(max(zqta(ig,l-1)-po(ig,l),0.))**1*(zalpha/w_est(ig,l+1))**0.5    &    
+!     &     +0. ))
+
 
 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
@@ -346,4 +311,8 @@
         endif
 
+!attention test
+!        if (detr_star(ig,l).gt.(f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l))) then       
+!            detr_star(ig,l)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
+!        endif
 ! Calcul du flux montant normalise
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
@@ -370,6 +339,6 @@
     enddo
 
-   ztemp(:)=zpspsk(:,l)*ztla(:,l)
-   call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l))
+   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
+
 
    do ig=1,ngrid
@@ -377,5 +346,4 @@
         if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
 ! on ecrit de maniere conservative (sat ou non)
-           zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l))
 !          T = Tl +Lv/Cp ql
            ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
@@ -387,4 +355,5 @@
 
 !on ecrit zqsat 
+           zqsatth(ig,l)=qsatbef  
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -407,25 +376,12 @@
            zw2(ig,l+1)=zw2modif+zdw2
 else
-! Tentative de reecriture de l'equation de w2. A reprendre ...
-!           zdw2=2*zdz*zbuoy
-!           zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon))
-!!!!!      zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
-!formulation Arnaud
-!           zw2fact=zbuoy*zalpha**expa/fact_epsilon
-!           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) &
-!      &    +zw2fact
-           if (zbuoy.gt.1.e-10) then
-           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact)
-           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
+           zdrag=fact_epsilon/(zalpha**expa)
+           zw2fact=zbuoy/zdrag*a1
+           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*zdrag/(1+fact_gamma)*zdz) &
       &    +zw2fact
-           else
-           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma)
-           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
-      &    +zw2fact
-
-           endif
+
 
 endif
-!           zw2(ig,l+1)=zw2modif+zdw2
+
       endif
    enddo
@@ -440,6 +396,5 @@
             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
 !               stop'On tombe sur le cas particulier de thermcell_dry'
-                write(lunout,*)                                         &
-     &          'On tombe sur le cas particulier de thermcell_plume'
+                print*,'On tombe sur le cas particulier de thermcell_plume'
                 zw2(ig,l+1)=0.
                 linter(ig)=l+1
@@ -487,5 +442,4 @@
       return 
       end
-
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
