Index: LMDZ4/trunk/libf/phylmd/calltherm.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/calltherm.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/calltherm.F90	(revision 1026)
@@ -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)
+     &       zmax0,f0,zw2,fraca)
 
       USE dimphy
@@ -39,4 +39,5 @@
       REAL wght_th(klon,klev)
       INTEGER lalim_conv(klon)
+      REAL zw2(klon,klev+1),fraca(klon,klev+1)
 
 !FH Update Thermiques
@@ -50,4 +51,5 @@
       real fmc_therm(klon,klev+1),zqasc(klon,klev)
       real zqla(klon,klev)
+      real zqta(klon,klev)
       real wmax_sec(klon)
       real zmax_sec(klon)
@@ -187,19 +189,4 @@
      &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
      &      ,tau_thermals)
-!            CALL calcul_sec_entr(klon,klev,zdt
-!     s      ,pplay,paprs,pphi,zlev,debut
-!     s      ,u_seri,v_seri,t_seri,q_seri                
-!     s      ,zmax_sec,wmax_sec,zw_sec,lmix_sec
-!     s      ,r_aspect_thermals,l_mix_thermals,w2di_thermals
-!     s      ,tau_thermals,3)
-!           CALL thermcell_pluie_detr(klon,klev,zdt  &
-!    &      ,pplay,paprs,pphi,zlev,debut  &
-!    &      ,u_seri,v_seri,t_seri,q_seri  &
-!    &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-!    &      ,zfm_therm,zentr_therm,zqla,lmax  &
-!    &      ,zmax_sec,wmax_sec,zw_sec,lmix_sec  &
-!    &      ,ratqscth,ratqsdiff,zqsatth  &
-!    &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
-!    &      ,tau_thermals)
           else if (iflag_thermals.ge.13) then
             CALL thermcell_main(itap,klon,klev,zdt  &
@@ -207,9 +194,9 @@
      &      ,u_seri,v_seri,t_seri,q_seri  &
      &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
-     &      ,zfm_therm,zentr_therm,zdetr_therm,zqla,lmax  &
+     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
      &      ,ratqscth,ratqsdiff,zqsatth  &
      &      ,r_aspect_thermals,l_mix_thermals  &
      &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
-     &      ,zmax0,f0)
+     &      ,zmax0,f0,zw2,fraca)
          endif
 
@@ -231,4 +218,6 @@
             entr_therm(:,k)=entr_therm(:,k)  &
      &       +zentr_therm(:,k)*fact(:)
+            detr_therm(:,k)=detr_therm(:,k)  &
+     &       +zdetr_therm(:,k)*fact(:)
       ENDDO
        fm_therm(:,klev+1)=0.
Index: LMDZ4/trunk/libf/phylmd/thermcell.h
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell.h	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell.h	(revision 1026)
@@ -3,4 +3,5 @@
       integer w2di_thermals,isplit
       integer iflag_coupl,iflag_clos,iflag_wake
+      integer iflag_thermals_ed,iflag_thermals_optflux
 
       common/ctherm1/iflag_thermals,nsplit_thermals
@@ -8,3 +9,5 @@
       common/ctherm3/w2di_thermals
       common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
+      common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux
+
 !$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm3/,/ctherm4/)
Index: LMDZ4/trunk/libf/phylmd/thermcell_closure.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell_closure.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell_closure.F90	(revision 1026)
@@ -1,4 +1,4 @@
       SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
+     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
 
 !-------------------------------------------------------------------------
@@ -8,4 +8,5 @@
 
 #include "iniprint.h"
+#include "thermcell.h"
       INTEGER ngrid,nlay
       INTEGER ig,k       
@@ -15,4 +16,5 @@
       INTEGER lalim(ngrid)
       REAL alim_star(ngrid,nlay)
+      REAL alim_star_tot(ngrid)
       REAL rho(ngrid,nlay)
       REAL zlev(ngrid,nlay)
@@ -48,13 +50,15 @@
                 stop
              endif
-             if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then 
-             f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
+             if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then 
+             f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
      &             *alim_star2(ig))
 !            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
 !    &                     zmax_sec(ig))*wmax_sec(ig))
+             print*,'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
              else
-             f(ig)=wmax(ig)/zdenom
+             f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
 !            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
 !     &                     zmax(ig))*wmax(ig))
+             print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
              endif
          endif
Index: LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90	(revision 1026)
@@ -12,4 +12,5 @@
       IMPLICIT NONE
 #include "iniprint.h"
+#include "thermcell.h"
       
       INTEGER ig,l
@@ -45,5 +46,4 @@
       REAL fomass_max,alphamax
       save fomass_max,alphamax
-!$OMP THREADPRIVATE(fomass_max,alphamax)
 
       fomass_max=0.5
@@ -202,5 +202,5 @@
 !Test sur fraca croissant
 !-------------------------------------------------------------------------
-      if (1.eq.1) then 
+      if (iflag_thermals_optflux==0) then 
 !     do l=1,klev
          do ig=1,ngrid
@@ -228,5 +228,5 @@
 !test sur flux de masse croissant
 !-------------------------------------------------------------------------
-      if (1.eq.1) then
+      if (iflag_thermals_optflux==0) then
 !     do l=1,klev
          do ig=1,ngrid
Index: LMDZ4/trunk/libf/phylmd/thermcell_height.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell_height.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell_height.F90	(revision 1026)
@@ -7,4 +7,5 @@
       IMPLICIT NONE
 #include "iniprint.h"
+#include "thermcell.h"
 
       INTEGER ig,l
@@ -23,4 +24,6 @@
       REAL zmax0(ngrid)
       REAL zmix(ngrid)
+      REAL num(ngrid)
+      REAL denom(ngrid)
 
       REAL zlevinter(ngrid)
@@ -70,4 +73,24 @@
          zlevinter(ig)=zlev(ig,1)
       enddo
+
+      if (iflag_thermals_ed.ge.1) then
+
+         num(:)=0.
+         denom(:)=0.
+         do ig=1,ngrid
+          do l=1,nlay
+             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+          enddo
+       enddo
+       do ig=1,ngrid
+       if (denom(ig).gt.1.e-10) then
+          zmax(ig)=2.*num(ig)/denom(ig)
+          zmax0(ig)=zmax(ig)
+       endif 
+       enddo
+
+       else
+
       do  ig=1,ngrid
 ! calcul de zlevinter
@@ -80,4 +103,8 @@
        zmax0(ig)=zmax(ig)
       enddo
+
+
+      endif
+!endif iflag_thermals_ed
 !
 ! def de  zmix continu (profil parabolique des vitesses)
Index: LMDZ4/trunk/libf/phylmd/thermcell_init.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell_init.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell_init.F90	(revision 1026)
@@ -7,4 +7,5 @@
       IMPLICIT NONE
 #include "iniprint.h"
+#include "thermcell.h"
 
       INTEGER l,ig
@@ -34,4 +35,24 @@
          lalim(ig)=1
       enddo
+
+      if (iflag_thermals_ed.ge.1) then
+!si la première couche est instable, on declenche un thermique
+         do ig=1,ngrid
+            if (ztv(ig,1).gt.ztv(ig,2)) then
+               lmin(ig)=1
+               lalim(ig)=2
+               alim_star(ig,1)=1.
+               alim_star_tot(ig)=alim_star(ig,1)
+               print*,'init',alim_star(ig,1),alim_star_tot(ig)
+            else
+                lmin(ig)=1
+                lalim(ig)=1
+                alim_star(ig,1)=0.
+                alim_star_tot(ig)=0. 
+            endif
+         enddo
+ 
+         else
+!else iflag_thermals_ed=0 ancienne def de l alim 
 
 !on ne considere que les premieres couches instables
@@ -131,4 +152,11 @@
       enddo
        
+!on remet alim_star_tot a 1
+      do ig=1,ngrid 
+         alim_star_tot(ig)=1.
+      enddo
+
+      endif
+!endif iflag_thermals_ed
       return 
       end  
Index: LMDZ4/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell_main.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell_main.F90	(revision 1026)
@@ -6,11 +6,12 @@
      &                  ,pu,pv,pt,po  &
      &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
-     &                  ,fm0,entr0,detr0,zqla,lmax  &
+     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
      &                  ,ratqscth,ratqsdiff,zqsatth  &
      &                  ,r_aspect,l_mix,tau_thermals &
      &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
-     &                  ,zmax0, f0)
+     &                  ,zmax0, f0,zw2,fraca)
 
       USE dimphy
+      USE comgeomphy , ONLY:rlond,rlatd
       IMPLICIT NONE
 
@@ -41,5 +42,4 @@
 
 #include "dimensions.h"
-!#include "dimphy.h"
 #include "YOMCST.h"
 #include "YOETHF.h"
@@ -82,9 +82,11 @@
       INTEGER lmax(klon),lmin(klon),lalim(klon)
       INTEGER lmix(klon)
+      INTEGER lmix_bis(klon)
       real linter(klon)
       real zmix(klon)
-      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev)
+      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1)
+!      real fraca(klon,klev)
+
       real zmax_sec(klon)
-      real w_est(klon,klev+1)
 !on garde le zmax du pas de temps precedent
       real zmax0(klon)
@@ -118,9 +120,4 @@
       real vardiff
       real ratqsdiff(klon,klev)
-      integer isplit,nsplit
-      parameter (nsplit=10)
-      data isplit/0/
-      save isplit
-!$OMP THREADPRIVATE(isplit)
 
       logical sorties
@@ -183,4 +180,20 @@
          entr0=0.
          detr0=0.
+
+
+! #define wrgrads_thermcell
+#undef wrgrads_thermcell
+#ifdef wrgrads_thermcell
+! Initialisation des sorties grads pour les thermiques.
+! Pour l'instant en 1D sur le point igout.
+! Utilise par thermcell_out3d.h
+         str10='therm'
+         call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
+     &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
+     &   ptimestep,str10,'therm ')
+#endif
+
+
+
       endif
 
@@ -409,7 +422,8 @@
 !IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
       CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
+     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
-     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
+     &           ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
+     &            ,lev_out,lunout1,igout)
       if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
 
@@ -445,7 +459,15 @@
 ! Fermeture,determination de f
 !-------------------------------------------------------------------------------
+!
+!avant closure: on redéfinit lalim, alim_star_tot et alim_star
+!       do ig=1,klon
+!       do l=2,lalim(ig)
+!       alim_star(ig,l)=entr_star(ig,l)
+!       entr_star(ig,l)=0.
+!       enddo
+!       enddo
 
       CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
+     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
 
       if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
@@ -757,17 +779,6 @@
 !     print*,'15 OK convect8'
 
-      isplit=isplit+1       
-
-
-#ifdef und
-      if (prt_level.ge.1) print*,'thermcell_main sorties 1D'
-#include "thermcell_out1d.h"
-#endif
-
-
-#define troisD
-#undef troisD
       if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
-#ifdef troisD
+#ifdef wrgrads_thermcell
 #include "thermcell_out3d.h"
 #endif
Index: LMDZ4/trunk/libf/phylmd/thermcell_plume.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/thermcell_plume.F90	(revision 1023)
+++ LMDZ4/trunk/libf/phylmd/thermcell_plume.F90	(revision 1026)
@@ -1,6 +1,7 @@
       SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
+     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
-     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
+     &           ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
+     &           ,lev_out,lunout1,igout)
 
 !--------------------------------------------------------------------------
@@ -14,7 +15,8 @@
 #include "FCTTRE.h"
 #include "iniprint.h"
+#include "thermcell.h"
 
       INTEGER itap
-
+      INTEGER lunout1,igout
       INTEGER ngrid,klev
       REAL ptimestep
@@ -36,4 +38,6 @@
       integer lev_out                           ! niveau pour les print
       real zcon2(ngrid)
+    
+      real alim_star_tot(ngrid)
 
       REAL ztva(ngrid,klev)
@@ -67,4 +71,5 @@
       REAL zqla_est(ngrid,klev)
       REAL zqsatth(ngrid,klev)
+      REAL zta_est(ngrid,klev)
 
       REAL linter(ngrid)
@@ -81,4 +86,6 @@
       PARAMETER (DDT0=.01)
       logical Zsat
+      REAL fact_gamma,fact_epsilon
+      REAL c2(ngrid,klev)
 
       Zsat=.false.
@@ -86,4 +93,15 @@
       RLvCp = RLVTT/RCPD
      
+      if (iflag_thermals_ed==0) then
+         fact_gamma=1.
+         fact_epsilon=1.
+      else if (iflag_thermals_ed==1)  then
+         fact_gamma=1.
+         fact_epsilon=1.
+      else if (iflag_thermals_ed==2)  then
+         fact_gamma=1.
+         fact_epsilon=2.
+      endif
+
       do l=1,klev
          do ig=1,ngrid
@@ -94,4 +112,10 @@
       enddo
 
+!CR: attention test couche alim
+!     do l=2,klev
+!     do ig=1,ngrid
+!        alim_star(ig,l)=0.
+!     enddo
+!     enddo
 !AM:initialisations du thermique
       do k=1,klev
@@ -141,4 +165,5 @@
          linter(ig)=1.
          lmix(ig)=1
+         lmix_bis(ig)=2
          wmaxa(ig)=0.
       enddo
@@ -212,4 +237,5 @@
      &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
      &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
+!     &      *1./3. &
      &      *(zlev(ig,3)-zlev(ig,2))
        endif
@@ -253,4 +279,5 @@
         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
         ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
+        zta_est(ig,l)=ztva_est(ig,l)
         ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
      &      -zqla_est(ig,l))-zqla_est(ig,l))
@@ -260,4 +287,5 @@
      &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
      &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+!     &                   *1./3. &
      &                   *(zlev(ig,l+1)-zlev(ig,l))
              if (w_est(ig,l+1).lt.0.) then
@@ -268,4 +296,6 @@
 !=======================
 
+!CR:on vire les modifs
+         if (iflag_thermals_ed==0) then
 
 ! Modifications du calcul du detrainement.
@@ -523,4 +553,86 @@
         if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
 
+       endif  ! iflag_thermals_ed==0
+
+!CR:nvlle def de entr_star et detr_star
+      if (iflag_thermals_ed>=1) then
+!      if (l.lt.lalim(ig)) then
+!      if (l.lt.2) then 
+!        entr_star(ig,l)=0.
+!        detr_star(ig,l)=0.
+!      else
+!      if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then 
+!         entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+!      else
+!         entr_star(ig,l)=  &
+!     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
+!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
+!     &                *(zlev(ig,l+1)-zlev(ig,l))
+
+ 
+         entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &          
+     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
+     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)   &
+     &                *(zlev(ig,l+1)-zlev(ig,l))) &
+     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+
+        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
+            alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
+            lalim(ig)=lmix_bis(ig)
+            print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
+        endif
+
+        if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
+!        c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
+         c2(ig,l)=0.001
+         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
+     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
+     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
+     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
+     &                *(zlev(ig,l+1)-zlev(ig,l)))                    &
+     &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+
+       else
+!         c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
+          c2(ig,l)=0.003
+
+         detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
+     &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
+     &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
+     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
+     &                *(zlev(ig,l+1)-zlev(ig,l))) &
+     &                +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+       endif
+         
+           
+!        detr_star(ig,l)=detr_star(ig,l)*3.
+!        if (l.lt.lalim(ig)) then
+!          entr_star(ig,l)=0.
+!        endif
+!        if (l.lt.2) then
+!          entr_star(ig,l)=0.
+!          detr_star(ig,l)=0.
+!        endif
+
+
+!      endif 
+!      else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
+!      entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1))        &
+!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
+!     &                *(zlev(ig,l+1)-zlev(ig,l))
+!      detr_star(ig,l)=0.002*f_star(ig,l)                         &
+!     &                *(zlev(ig,l+1)-zlev(ig,l))
+!      else
+!      entr_star(ig,l)=0.001*f_star(ig,l)                         &
+!     &                *(zlev(ig,l+1)-zlev(ig,l))
+!      detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1))       &
+!     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))       &
+!     &                *(zlev(ig,l+1)-zlev(ig,l))                      &
+!     &                +0.002*f_star(ig,l)                             &
+!     &                *(zlev(ig,l+1)-zlev(ig,l))
+!      endif
+
+      endif   ! iflag_thermals_ed==1
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! FH inutile si on conserve comme on l'a fait plus haut entr=detr
@@ -533,4 +645,5 @@
 !
 !prise en compte du detrainement et de l entrainement dans le calcul du flux
+
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
      &              -detr_star(ig,l)
@@ -597,8 +710,34 @@
            zw2(ig,l+1)=zw2(ig,l)*  &
      &                 ((f_star(ig,l))**2)  &
-     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
+!  Tests de Catherine
+!     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
+     &      /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
      &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+     &                 *fact_gamma &
      &                 *(zlev(ig,l+1)-zlev(ig,l))
-               
+!prise en compte des forces de pression que qd flottabilité<0
+!              zw2(ig,l+1)=zw2(ig,l)*  &
+!     &            1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &        
+!     &                 (f_star(ig,l))**2 &
+!     &                 /(f_star(ig,l)+entr_star(ig,l))**2+ &
+!     &                 (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+  &        
+!     &                 /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
+!     &                 /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
+!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+!     &                 *1./3. &
+!     &                 *(zlev(ig,l+1)-zlev(ig,l))
+          
+!        write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
+!     &              -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
+!     &               2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
+
+ 
+!             zw2(ig,l+1)=zw2(ig,l)*  &
+!     &                 (2.-2.*entr_star(ig,l)/f_star(ig,l)) &  
+!     &                 -zw2(ig,l-1)+  &        
+!     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+!     &                 *1./3. &
+!     &                 *(zlev(ig,l+1)-zlev(ig,l))             
+
             endif
         endif
@@ -614,5 +753,5 @@
             endif
 
-
+!        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
         if (zw2(ig,l+1).lt.0.) then 
            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
@@ -625,4 +764,8 @@
         if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
+!on rajoute le calcul de lmix_bis
+            if (zqla(ig,l).lt.1.e-10) then
+               lmix_bis(ig)=l+1
+            endif
             lmix(ig)=l+1
             wmaxa(ig)=wa_moy(ig,l+1)
@@ -631,17 +774,26 @@
       enddo
 
+!on remplace a* par e* ds premiere couche
+!      if (iflag_thermals_ed.ge.1) then
+!       do ig=1,ngrid
+!       do l=2,klev
+!          if (l.lt.lalim(ig)) then
+!             alim_star(ig,l)=entr_star(ig,l)
+!          endif
+!       enddo
+!       enddo
+!       do ig=1,ngrid
+!          lalim(ig)=lmix_bis(ig)
+!       enddo
+!      endif
+       if (iflag_thermals_ed.ge.1) then
+          do ig=1,ngrid
+             do l=2,lalim(ig)
+                alim_star(ig,l)=entr_star(ig,l)
+                entr_star(ig,l)=0.
+             enddo
+           enddo
+       endif
         if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
-#ifdef troisD
-      call wrgradsfi(1,klev,zqla_est,'ql_es     ','ql_es     ')
-      call wrgradsfi(1,klev,entr_star1,'e_st1     ','e_st1     ')
-      call wrgradsfi(1,klev,entr_star2,'e_st2     ','e_st2     ')
-      call wrgradsfi(1,klev,detr_stara,'d_sta     ','d_sta     ')
-      call wrgradsfi(1,klev,detr_starb,'d_stb     ','d_stb     ')
-      call wrgradsfi(1,klev,detr_starc,'d_stc     ','d_stc     ')
-      call wrgradsfi(1,klev,zqla0     ,'zqla0     ','zqla0     ')
-      call wrgradsfi(1,klev,detr_star0,'d_st0     ','d_st0     ')
-      call wrgradsfi(1,klev,detr_star1,'d_st1     ','d_st1     ')
-      call wrgradsfi(1,klev,detr_star2,'d_st2     ','d_st2     ')
-#endif
 
 !     print*,'thermcell_plume OK'
