Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90	(revision 1294)
@@ -45,4 +45,5 @@
 !********************************************************
 !     declarations
+      LOGICAL flag_bidouille_stratocu
       real fmc_therm(klon,klev+1),zqasc(klon,klev)
       real zqla(klon,klev)
@@ -187,6 +188,6 @@
      &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
      &      ,tau_thermals)
-          else if (iflag_thermals.ge.13) then
-            CALL thermcell_main(itap,klon,klev,zdt  &
+          else if (iflag_thermals==13.or.iflag_thermals==14) then
+            CALL thermcellV0_main(itap,klon,klev,zdt  &
      &      ,pplay,paprs,pphi,debut  &
      &      ,u_seri,v_seri,t_seri,q_seri  &
@@ -197,10 +198,24 @@
      &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
      &      ,zmax0,f0,zw2,fraca)
+          else if (iflag_thermals==15.or.iflag_thermals==16) then
+            CALL thermcell_main(itap,klon,klev,zdt  &
+     &      ,pplay,paprs,pphi,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,zdetr_therm,zqasc,zqla,lmax  &
+     &      ,ratqscth,ratqsdiff,zqsatth  &
+     &      ,r_aspect_thermals,l_mix_thermals  &
+     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
+     &      ,zmax0,f0,zw2,fraca)
+           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
+         else
+            STOP'Cas des thermiques non prevu'
          endif
 
+       flag_bidouille_stratocu=iflag_thermals.lt.14.or.iflag_thermals.lt.16
 
       fact(:)=0.
       DO i=1,klon
-       logexpr1(i)=iflag_thermals.lt.14.or.weak_inversion(i).gt.0.5
+       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
        IF(logexpr1(i)) fact(i)=1./float(nsplit_thermals)
       ENDDO
@@ -235,4 +250,5 @@
             qmemoire(:,:)=q_seri(:,:)
             q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
+           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
 
        DO i=1,klon
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1294)
@@ -2548,5 +2548,5 @@
 
 !   les ratqs sont une combinaison de ratqss et ratqsc
-       print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
+       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
 
          if (tau_ratqs>1.e-10) then
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcellV0_main.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcellV0_main.F90	(revision 1294)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcellV0_main.F90	(revision 1294)
@@ -0,0 +1,2088 @@
+!
+! $Header$
+!
+      SUBROUTINE thermcellV0_main(itap,ngrid,nlay,ptimestep  &
+     &                  ,pplay,pplev,pphi,debut  &
+     &                  ,pu,pv,pt,po  &
+     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
+     &                  ,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,zw2,fraca)
+
+      USE dimphy
+      USE comgeomphy , ONLY:rlond,rlatd
+      IMPLICIT NONE
+
+!=======================================================================
+!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
+!   Version du 09.02.07
+!   Calcul du transport vertical dans la couche limite en presence
+!   de "thermiques" explicitement representes avec processus nuageux
+!
+!   Réécriture à partir d'un listing papier à Habas, le 14/02/00
+!
+!   le thermique est supposé homogène et dissipé par mélange avec
+!   son environnement. la longueur l_mix contrôle l'efficacité du
+!   mélange
+!
+!   Le calcul du transport des différentes espèces se fait en prenant
+!   en compte:
+!     1. un flux de masse montant
+!     2. un flux de masse descendant
+!     3. un entrainement
+!     4. un detrainement
+!
+!=======================================================================
+
+!-----------------------------------------------------------------------
+!   declarations:
+!   -------------
+
+#include "dimensions.h"
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "iniprint.h"
+
+!   arguments:
+!   ----------
+
+!IM 140508
+      INTEGER itap
+
+      INTEGER ngrid,nlay,w2di
+      real tau_thermals
+      real ptimestep,l_mix,r_aspect
+      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
+      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
+      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
+      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
+      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
+      real pphi(ngrid,nlay)
+
+!   local:
+!   ------
+
+      integer icount
+      data icount/0/
+      save icount
+!$OMP THREADPRIVATE(icount)
+
+      integer,save :: igout=1
+!$OMP THREADPRIVATE(igout)
+      integer,save :: lunout1=6
+!$OMP THREADPRIVATE(lunout1)
+      integer,save :: lev_out=10
+!$OMP THREADPRIVATE(lev_out)
+
+      INTEGER ig,k,l,ll
+      real zsortie1d(klon)
+      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),zw_est(klon,klev+1)
+!      real fraca(klon,klev)
+
+      real zmax_sec(klon)
+!on garde le zmax du pas de temps precedent
+      real zmax0(klon)
+!FH/IM     save zmax0
+
+      real lambda
+
+      real zlev(klon,klev+1),zlay(klon,klev)
+      real deltaz(klon,klev)
+      REAL zh(klon,klev)
+      real zthl(klon,klev),zdthladj(klon,klev)
+      REAL ztv(klon,klev)
+      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
+      real zl(klon,klev)
+      real zsortie(klon,klev)
+      real zva(klon,klev)
+      real zua(klon,klev)
+      real zoa(klon,klev)
+
+      real zta(klon,klev)
+      real zha(klon,klev)
+      real fraca(klon,klev+1)
+      real zf,zf2
+      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
+      real q2(klon,klev)
+! FH probleme de dimensionnement avec l'allocation dynamique
+!     common/comtherm/thetath2,wth2
+    
+      real ratqscth(klon,klev)
+      real var
+      real vardiff
+      real ratqsdiff(klon,klev)
+
+      logical sorties
+      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
+      real zpspsk(klon,klev)
+
+      real wmax(klon)
+      real wmax_sec(klon)
+      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
+      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
+
+      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
+!niveau de condensation
+      integer nivcon(klon)
+      real zcon(klon)
+      REAL CHI
+      real zcon2(klon)
+      real pcon(klon)
+      real zqsat(klon,klev)
+      real zqsatth(klon,klev) 
+
+      real f_star(klon,klev+1),entr_star(klon,klev)
+      real detr_star(klon,klev)
+      real alim_star_tot(klon),alim_star2(klon)
+      real alim_star(klon,klev)
+      real f(klon), f0(klon)
+!FH/IM     save f0
+      real zlevinter(klon)
+      logical debut
+       real seuil
+
+!
+      !nouvelles variables pour la convection
+      real Ale_bl(klon)
+      real Alp_bl(klon)
+      real alp_int(klon)
+      real ale_int(klon)
+      integer n_int(klon)
+      real fm_tot(klon)
+      real wght_th(klon,klev)
+      integer lalim_conv(klon)
+!v1d     logical therm
+!v1d     save therm
+
+      character*2 str2
+      character*10 str10
+
+      EXTERNAL SCOPY
+!
+
+!-----------------------------------------------------------------------
+!   initialisation:
+!   ---------------
+!
+
+      seuil=0.25
+
+      if (debut)  then
+         fm0=0.
+         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
+
+      fm=0. ; entr=0. ; detr=0.
+
+      icount=icount+1
+
+!IM 090508 beg
+!print*,'====================================================================='
+!print*,'====================================================================='
+!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
+!print*,'====================================================================='
+!print*,'====================================================================='
+!IM 090508 end
+
+      if (prt_level.ge.1) print*,'thermcell_main V4'
+
+       sorties=.true.
+      IF(ngrid.NE.klon) THEN
+         PRINT*
+         PRINT*,'STOP dans convadj'
+         PRINT*,'ngrid    =',ngrid
+         PRINT*,'klon  =',klon
+      ENDIF
+!
+!Initialisation
+!
+!    IF (1.eq.0) THEN
+!     do ig=1,klon      
+!FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
+!     if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
+!           f0(ig)=1.e-5
+!           zmax0(ig)=40.
+!v1d        therm=.false.
+!     endif
+!     enddo 
+!    ENDIF !(1.eq.0) THEN
+     if (prt_level.ge.10)write(lunout,*)                                &
+    &     'WARNING thermcell_main f0=max(f0,1.e-2)'
+     do ig=1,klon
+      if (prt_level.ge.20) then
+       print*,'th_main ig f0',ig,f0(ig)
+      endif
+         f0(ig)=max(f0(ig),1.e-2)
+!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
+     enddo
+
+!-----------------------------------------------------------------------
+! Calcul de T,q,ql a partir de Tl et qT dans l environnement
+!   --------------------------------------------------------------------
+!
+      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
+     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
+       
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
+
+!------------------------------------------------------------------------
+!                       --------------------
+!
+!
+!                       + + + + + + + + + + +
+!
+!
+!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
+!  wh,wt,wo ...
+!
+!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
+!
+!
+!                       --------------------   zlev(1)
+!                       \\\\\\\\\\\\\\\\\\\\
+!
+!
+
+!-----------------------------------------------------------------------
+!   Calcul des altitudes des couches
+!-----------------------------------------------------------------------
+
+      do l=2,nlay
+         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
+      enddo
+         zlev(:,1)=0.
+         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
+      do l=1,nlay
+         zlay(:,l)=pphi(:,l)/RG
+      enddo
+!calcul de l epaisseur des couches
+      do l=1,nlay
+         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
+      enddo
+
+!     print*,'2 OK convect8'
+!-----------------------------------------------------------------------
+!   Calcul des densites
+!-----------------------------------------------------------------------
+
+      do l=1,nlay
+         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
+      enddo
+
+!IM
+     if (prt_level.ge.10)write(lunout,*)                                &
+    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
+      rhobarz(:,1)=rho(:,1)
+
+      do l=2,nlay
+         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
+      enddo
+
+!calcul de la masse
+      do l=1,nlay
+         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
+      enddo
+
+      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
+
+!------------------------------------------------------------------
+!
+!             /|\
+!    --------  |  F_k+1 -------   
+!                              ----> D_k
+!             /|\              <---- E_k , A_k
+!    --------  |  F_k --------- 
+!                              ----> D_k-1
+!                              <---- E_k-1 , A_k-1
+!
+!
+!
+!
+!
+!    ---------------------------
+!
+!    ----- F_lmax+1=0 ----------         \
+!            lmax     (zmax)              |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------          |
+!                                         |  E
+!    ---------------------------          |  D
+!                                         |
+!    ---------------------------          |
+!                                         |
+!    ---------------------------  \       |
+!            lalim                 |      |
+!    ---------------------------   |      |
+!                                  |      |
+!    ---------------------------   |      |
+!                                  | A    |
+!    ---------------------------   |      |
+!                                  |      |
+!    ---------------------------   |      |
+!    lmin  (=1 pour le moment)     |      |
+!    ----- F_lmin=0 ------------  /      /
+!
+!    ---------------------------
+!    //////////////////////////
+!
+!
+!=============================================================================
+!  Calculs initiaux ne faisant pas intervenir les changements de phase
+!=============================================================================
+
+!------------------------------------------------------------------
+!  1. alim_star est le profil vertical de l'alimentation à la base du
+!     panache thermique, calculé à partir de la flotabilité de l'air sec
+!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
+!------------------------------------------------------------------
+!
+      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
+      CALL thermcellV0_init(ngrid,nlay,ztv,zlay,zlev,  &
+     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
+
+call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
+call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
+
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
+      if (prt_level.ge.10) then
+         write(lunout1,*) 'Dans thermcell_main 1'
+         write(lunout1,*) 'lmin ',lmin(igout)
+         write(lunout1,*) 'lalim ',lalim(igout)
+         write(lunout1,*) ' ig l alim_star thetav'
+         write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
+     &   ,ztv(igout,l),l=1,lalim(igout)+4)
+      endif
+
+!v1d      do ig=1,klon
+!v1d     if (alim_star(ig,1).gt.1.e-10) then
+!v1d     therm=.true.
+!v1d     endif
+!v1d      enddo
+!-----------------------------------------------------------------------------
+!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
+!     panache sec conservatif (e=d=0) alimente selon alim_star 
+!     Il s'agit d'un calcul de type CAPE
+!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
+!------------------------------------------------------------------------------
+!
+      CALL thermcellV0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
+     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
+
+call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
+call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
+      if (prt_level.ge.10) then
+         write(lunout1,*) 'Dans thermcell_main 1b'
+         write(lunout1,*) 'lmin ',lmin(igout)
+         write(lunout1,*) 'lalim ',lalim(igout)
+         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
+         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
+     &    ,l=1,lalim(igout)+4)
+      endif
+
+
+
+!---------------------------------------------------------------------------------
+!calcul du melange et des variables dans le thermique
+!--------------------------------------------------------------------------------
+!
+      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
+!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
+      CALL thermcellV0_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
+     &           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,zw_est,zqsatth,lmix,lmix_bis,linter &
+     &            ,lev_out,lunout1,igout)
+      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
+
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
+      if (prt_level.ge.10) then
+         write(lunout1,*) 'Dans thermcell_main 2'
+         write(lunout1,*) 'lmin ',lmin(igout)
+         write(lunout1,*) 'lalim ',lalim(igout)
+         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
+         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
+     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
+      endif
+
+!-------------------------------------------------------------------------------
+! Calcul des caracteristiques du thermique:zmax,zmix,wmax
+!-------------------------------------------------------------------------------
+!
+      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
+     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
+
+
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
+
+!-------------------------------------------------------------------------------
+! 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 thermcellV0_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
+     &   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'
+
+      if (tau_thermals>1.) then
+         lambda=exp(-ptimestep/tau_thermals)
+         f0=(1.-lambda)*f+lambda*f0
+      else
+         f0=f
+      endif
+
+! Test valable seulement en 1D mais pas genant
+      if (.not. (f0(1).ge.0.) ) then
+           stop 'Dans thermcell_main'
+      endif
+
+!-------------------------------------------------------------------------------
+!deduction des flux
+!-------------------------------------------------------------------------------
+
+      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
+     &       lalim,lmax,alim_star,  &
+     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
+     &       detr,zqla,lev_out,lunout1,igout)
+!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
+      call testV0_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
+
+!------------------------------------------------------------------
+!   On ne prend pas directement les profils issus des calculs precedents
+!   mais on s'autorise genereusement une relaxation vers ceci avec
+!   une constante de temps tau_thermals (typiquement 1800s).
+!------------------------------------------------------------------
+
+      if (tau_thermals>1.) then
+         lambda=exp(-ptimestep/tau_thermals)
+         fm0=(1.-lambda)*fm+lambda*fm0
+         entr0=(1.-lambda)*entr+lambda*entr0
+!        detr0=(1.-lambda)*detr+lambda*detr0
+      else
+         fm0=fm
+         entr0=entr
+         detr0=detr
+      endif
+
+!c------------------------------------------------------------------
+!   calcul du transport vertical
+!------------------------------------------------------------------
+
+      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
+     &                    zthl,zdthladj,zta,lev_out)
+      call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
+     &                   po,pdoadj,zoa,lev_out)
+
+!------------------------------------------------------------------
+! Calcul de la fraction de l'ascendance
+!------------------------------------------------------------------
+      do ig=1,klon
+         fraca(ig,1)=0.
+         fraca(ig,nlay+1)=0.
+      enddo
+      do l=2,nlay
+         do ig=1,klon
+            if (zw2(ig,l).gt.1.e-10) then
+            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
+            else
+            fraca(ig,l)=0.
+            endif
+         enddo
+      enddo
+     
+!------------------------------------------------------------------
+!  calcul du transport vertical du moment horizontal
+!------------------------------------------------------------------
+
+!IM 090508  
+      if (1.eq.1) then
+!IM 070508 vers. _dq       
+!     if (1.eq.0) then
+
+
+! Calcul du transport de V tenant compte d'echange par gradient
+! de pression horizontal avec l'environnement
+
+         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
+     &    ,fraca,zmax  &
+     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
+!IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
+      else
+
+! calcul purement conservatif pour le transport de V
+         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
+     &    ,zu,pduadj,zua,lev_out)
+         call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
+     &    ,zv,pdvadj,zva,lev_out)
+      endif
+
+!     print*,'13 OK convect8'
+      do l=1,nlay
+         do ig=1,ngrid
+           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)  
+         enddo
+      enddo
+
+      if (prt_level.ge.1) print*,'14 OK convect8'
+!------------------------------------------------------------------
+!   Calculs de diagnostiques pour les sorties
+!------------------------------------------------------------------
+!calcul de fraca pour les sorties
+      
+      if (sorties) then
+      if (prt_level.ge.1) print*,'14a OK convect8'
+! calcul du niveau de condensation
+! initialisation
+      do ig=1,ngrid
+         nivcon(ig)=0
+         zcon(ig)=0.
+      enddo 
+!nouveau calcul
+      do ig=1,ngrid
+      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
+      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
+      enddo
+      do k=1,nlay
+         do ig=1,ngrid
+         if ((pcon(ig).le.pplay(ig,k))  &
+     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
+            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
+         endif
+         enddo
+      enddo
+      if (prt_level.ge.1) print*,'14b OK convect8'
+      do k=nlay,1,-1
+         do ig=1,ngrid
+            if (zqla(ig,k).gt.1e-10) then
+               nivcon(ig)=k
+               zcon(ig)=zlev(ig,k)
+            endif
+         enddo
+      enddo
+      if (prt_level.ge.1) print*,'14c OK convect8'
+!calcul des moments
+!initialisation
+      do l=1,nlay
+         do ig=1,ngrid
+            q2(ig,l)=0.
+            wth2(ig,l)=0.
+            wth3(ig,l)=0.
+            ratqscth(ig,l)=0.
+            ratqsdiff(ig,l)=0.
+         enddo
+      enddo      
+      if (prt_level.ge.1) print*,'14d OK convect8'
+      if (prt_level.ge.10)write(lunout,*)                                &
+    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
+      do l=1,nlay
+         do ig=1,ngrid
+            zf=fraca(ig,l)
+            zf2=zf/(1.-zf)
+!
+      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
+!
+      if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
+            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
+            if(zw2(ig,l).gt.1.e-10) then
+             wth2(ig,l)=zf2*(zw2(ig,l))**2
+            else
+             wth2(ig,l)=0.
+            endif
+!           print*,'wth2=',wth2(ig,l)
+            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
+     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
+      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
+            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
+!test: on calcul q2/po=ratqsc
+            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
+         enddo
+      enddo
+!calcul de ale_bl et alp_bl
+!pour le calcul d'une valeur intégrée entre la surface et lmax
+      do ig=1,ngrid
+      alp_int(ig)=0.
+      ale_int(ig)=0.
+      n_int(ig)=0
+      enddo
+!
+      do l=1,nlay
+      do ig=1,ngrid
+       if(l.LE.lmax(ig)) THEN
+        alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
+        ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
+        n_int(ig)=n_int(ig)+1
+       endif
+      enddo
+      enddo
+!      print*,'avant calcul ale et alp' 
+!calcul de ALE et ALP pour la convection
+      do ig=1,ngrid
+!      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
+!          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
+!      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig)) 
+!     &           *0.1
+!valeur integree de alp_bl * 0.5:
+       if (n_int(ig).gt.0) then
+       Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
+!       if (Alp_bl(ig).lt.0.) then
+!       Alp_bl(ig)=0.
+       endif
+!       endif
+!         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
+!     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
+!            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
+!      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
+!      if (nivcon(ig).eq.1) then
+!       Ale_bl(ig)=0.
+!       else
+!valeur max de ale_bl:
+       Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2 
+!     & /2.
+!     & *0.1
+!        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2 
+!       if (n_int(ig).gt.0) then
+!       Ale_bl(ig)=ale_int(ig)/n_int(ig)
+!        Ale_bl(ig)=4.
+!       endif
+!       endif
+!            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
+!          Ale_bl(ig)=wth2(ig,nivcon(ig))
+!          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
+      enddo
+!test:calcul de la ponderation des couches pour KE
+!initialisations
+!      print*,'ponderation'
+      do ig=1,ngrid
+           fm_tot(ig)=0.
+      enddo
+       do ig=1,ngrid
+        do k=1,klev
+           wght_th(ig,k)=1.
+        enddo
+       enddo
+       do ig=1,ngrid
+!         lalim_conv(ig)=lmix_bis(ig)
+!la hauteur de la couche alim_conv = hauteur couche alim_therm
+         lalim_conv(ig)=lalim(ig)
+!         zentr(ig)=zlev(ig,lalim(ig))
+      enddo
+      do ig=1,ngrid
+        do k=1,lalim_conv(ig)
+           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
+        enddo
+      enddo
+      do ig=1,ngrid
+        do k=1,lalim_conv(ig)
+           if (fm_tot(ig).gt.1.e-10) then
+!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
+           endif
+!on pondere chaque couche par a*
+             if (alim_star(ig,k).gt.1.e-10) then
+             wght_th(ig,k)=alim_star(ig,k)
+             else
+             wght_th(ig,k)=1.
+             endif
+        enddo
+      enddo
+!      print*,'apres wght_th'
+!test pour prolonger la convection
+      do ig=1,ngrid
+!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
+      if ((alim_star(ig,1).lt.1.e-10)) then
+      lalim_conv(ig)=1
+      wght_th(ig,1)=1.
+!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
+      endif
+      enddo
+
+!calcul du ratqscdiff
+      if (prt_level.ge.1) print*,'14e OK convect8'
+      var=0.
+      vardiff=0.
+      ratqsdiff(:,:)=0.
+      do ig=1,ngrid
+         do l=1,lalim(ig)
+            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
+         enddo
+      enddo
+      if (prt_level.ge.1) print*,'14f OK convect8'
+      do ig=1,ngrid
+          do l=1,lalim(ig)
+          zf=fraca(ig,l)
+          zf2=zf/(1.-zf)
+          vardiff=vardiff+alim_star(ig,l)  &
+     &           *(zqta(ig,l)*1000.-var)**2
+!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
+!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
+          enddo
+      enddo
+      if (prt_level.ge.1) print*,'14g OK convect8'
+      do l=1,nlay
+         do ig=1,ngrid
+            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
+!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
+         enddo
+      enddo 
+!--------------------------------------------------------------------    
+!
+!ecriture des fichiers sortie
+!     print*,'15 OK convect8'
+
+      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
+#ifdef wrgrads_thermcell
+#include "thermcell_out3d.h"
+#endif
+
+      endif
+
+      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
+
+!     if(icount.eq.501) stop'au pas 301 dans thermcell_main'
+      return
+      end
+
+!-----------------------------------------------------------------------------
+
+      subroutine testV0_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
+      IMPLICIT NONE
+#include "iniprint.h"
+
+      integer i, k, klon,klev
+      real pplev(klon,klev+1),pplay(klon,klev)
+      real ztv(klon,klev)
+      real po(klon,klev)
+      real ztva(klon,klev)
+      real zqla(klon,klev)
+      real f_star(klon,klev)
+      real zw2(klon,klev)
+      integer long(klon)
+      real seuil
+      character*21 comment
+
+      if (prt_level.ge.1) THEN
+       print*,'WARNING !!! TEST ',comment
+      endif
+      return
+
+!  test sur la hauteur des thermiques ...
+         do i=1,klon
+!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
+           if (prt_level.ge.10) then
+               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
+               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
+               do k=1,klev
+                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
+               enddo
+!              stop
+           endif
+         enddo
+
+
+      return
+      end
+
+!==============================================================================
+      SUBROUTINE thermcellV0_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
+     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
+
+!-------------------------------------------------------------------------
+!thermcell_closure: fermeture, determination de f
+!-------------------------------------------------------------------------
+      IMPLICIT NONE
+
+#include "iniprint.h"
+#include "thermcell.h"
+      INTEGER ngrid,nlay
+      INTEGER ig,k       
+      REAL r_aspect,ptimestep
+      integer lev_out                           ! niveau pour les print
+
+      INTEGER lalim(ngrid)
+      REAL alim_star(ngrid,nlay)
+      REAL alim_star_tot(ngrid)
+      REAL rho(ngrid,nlay)
+      REAL zlev(ngrid,nlay)
+      REAL zmax(ngrid),zmax_sec(ngrid)
+      REAL wmax(ngrid),wmax_sec(ngrid)
+      real zdenom
+
+      REAL alim_star2(ngrid)
+
+      REAL f(ngrid)
+
+      do ig=1,ngrid
+         alim_star2(ig)=0.
+      enddo
+      do ig=1,ngrid
+         if (alim_star(ig,1).LT.1.e-10) then
+            f(ig)=0.
+         else   
+             do k=1,lalim(ig)
+                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
+     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
+             enddo
+             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
+             if (zdenom<1.e-14) then
+                print*,'ig=',ig
+                print*,'alim_star2',alim_star2(ig)
+                print*,'zmax',zmax(ig)
+                print*,'r_aspect',r_aspect
+                print*,'zdenom',zdenom
+                print*,'alim_star',alim_star(ig,:)
+                print*,'zmax_sec',zmax_sec(ig)
+                print*,'wmax_sec',wmax_sec(ig)
+                stop
+             endif
+             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))
+             if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
+             else
+             f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
+!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
+!     &                     zmax(ig))*wmax(ig))
+             if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
+             endif
+         endif
+!         f0(ig)=f(ig)
+      enddo
+      if (prt_level.ge.1) print*,'apres fermeture'
+
+!
+      return
+      end
+!==============================================================================
+      SUBROUTINE thermcellV0_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
+     &           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,w_est,zqsatth,lmix,lmix_bis,linter &
+     &           ,lev_out,lunout1,igout)
+
+!--------------------------------------------------------------------------
+!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
+!--------------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+#include "iniprint.h"
+#include "thermcell.h"
+
+      INTEGER itap
+      INTEGER lunout1,igout
+      INTEGER ngrid,klev
+      REAL ptimestep
+      REAL ztv(ngrid,klev)
+      REAL zthl(ngrid,klev)
+      REAL po(ngrid,klev)
+      REAL zl(ngrid,klev)
+      REAL rhobarz(ngrid,klev)
+      REAL zlev(ngrid,klev+1)
+      REAL pplev(ngrid,klev+1)
+      REAL pphi(ngrid,klev)
+      REAL zpspsk(ngrid,klev)
+      REAL alim_star(ngrid,klev)
+      REAL zmax_sec(ngrid)
+      REAL f0(ngrid)
+      REAL l_mix
+      REAL r_aspect
+      INTEGER lalim(ngrid)
+      integer lev_out                           ! niveau pour les print
+      real zcon2(ngrid)
+    
+      real alim_star_tot(ngrid)
+
+      REAL ztva(ngrid,klev)
+      REAL ztla(ngrid,klev)
+      REAL zqla(ngrid,klev)
+      REAL zqla0(ngrid,klev)
+      REAL zqta(ngrid,klev)
+      REAL zha(ngrid,klev)
+
+      REAL detr_star(ngrid,klev)
+      REAL coefc
+      REAL detr_stara(ngrid,klev)
+      REAL detr_starb(ngrid,klev)
+      REAL detr_starc(ngrid,klev)
+      REAL detr_star0(ngrid,klev)
+      REAL detr_star1(ngrid,klev)
+      REAL detr_star2(ngrid,klev)
+
+      REAL entr_star(ngrid,klev)
+      REAL entr_star1(ngrid,klev)
+      REAL entr_star2(ngrid,klev)
+      REAL detr(ngrid,klev)
+      REAL entr(ngrid,klev)
+
+      REAL zw2(ngrid,klev+1)
+      REAL w_est(ngrid,klev+1)
+      REAL f_star(ngrid,klev+1)
+      REAL wa_moy(ngrid,klev+1)
+
+      REAL ztva_est(ngrid,klev)
+      REAL zqla_est(ngrid,klev)
+      REAL zqsatth(ngrid,klev)
+      REAL zta_est(ngrid,klev)
+
+      REAL linter(ngrid)
+      INTEGER lmix(ngrid)
+      INTEGER lmix_bis(ngrid)
+      REAL    wmaxa(ngrid)
+
+      INTEGER ig,l,k
+
+      real zcor,zdelta,zcvm5,qlbef
+      real Tbef,qsatbef
+      real dqsat_dT,DT,num,denom
+      REAL REPS,RLvCp,DDT0
+      PARAMETER (DDT0=.01)
+      logical Zsat
+      REAL fact_gamma,fact_epsilon
+      REAL c2(ngrid,klev)
+
+      Zsat=.false.
+! Initialisation
+      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
+            zqla_est(ig,l)=0.
+            ztva_est(ig,l)=ztva(ig,l)
+            zqsatth(ig,l)=0.
+         enddo
+      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
+         do ig=1,ngrid
+            ztva(ig,k)=ztv(ig,k)
+            ztla(ig,k)=zthl(ig,k)
+            zqla(ig,k)=0.
+            zqta(ig,k)=po(ig,k)
+!
+            ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
+            ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
+            zha(ig,k) = ztva(ig,k)
+!
+         enddo
+      enddo 
+      do k=1,klev
+        do ig=1,ngrid
+           detr_star(ig,k)=0.
+           entr_star(ig,k)=0.
+
+           detr_stara(ig,k)=0.
+           detr_starb(ig,k)=0.
+           detr_starc(ig,k)=0.
+           detr_star0(ig,k)=0.
+           zqla0(ig,k)=0.
+           detr_star1(ig,k)=0.
+           detr_star2(ig,k)=0.
+           entr_star1(ig,k)=0.
+           entr_star2(ig,k)=0.
+
+           detr(ig,k)=0.
+           entr(ig,k)=0.
+        enddo
+      enddo
+      if (prt_level.ge.1) print*,'7 OK convect8'
+      do k=1,klev+1
+         do ig=1,ngrid
+            zw2(ig,k)=0.
+            w_est(ig,k)=0.
+            f_star(ig,k)=0.
+            wa_moy(ig,k)=0.
+         enddo
+      enddo
+
+      if (prt_level.ge.1) print*,'8 OK convect8'
+      do ig=1,ngrid
+         linter(ig)=1.
+         lmix(ig)=1
+         lmix_bis(ig)=2
+         wmaxa(ig)=0.
+      enddo
+
+!-----------------------------------------------------------------------------------
+!boucle de calcul de la vitesse verticale dans le thermique
+!-----------------------------------------------------------------------------------
+      do l=1,klev-1
+         do ig=1,ngrid
+
+
+
+! Calcul dans la premiere couche active du thermique (ce qu'on teste
+! en disant que la couche est instable et que w2 en bas de la couche
+! est nulle.
+
+            if (ztv(ig,l).gt.ztv(ig,l+1)  &
+     &         .and.alim_star(ig,l).gt.1.e-10  &
+     &         .and.zw2(ig,l).lt.1e-10) then
+
+
+! Le panache va prendre au debut les caracteristiques de l'air contenu
+! dans cette couche.
+               ztla(ig,l)=zthl(ig,l) 
+               zqta(ig,l)=po(ig,l)
+               zqla(ig,l)=zl(ig,l)
+               f_star(ig,l+1)=alim_star(ig,l)
+
+               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
+     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
+     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
+               w_est(ig,l+1)=zw2(ig,l+1)
+!
+
+
+            else if ((zw2(ig,l).ge.1e-10).and.  &
+     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
+!estimation du detrainement a partir de la geometrie du pas precedent
+!tests sur la definition du detr
+!calcul de detr_star et entr_star
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH le test miraculeux de Catherine ? Le bout du tunel ?
+!               w_est(ig,3)=zw2(ig,2)*  &
+!    &                   ((f_star(ig,2))**2)  &
+!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
+!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
+!    &                   *(zlev(ig,3)-zlev(ig,2))
+!     if (l.gt.2) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+
+! Premier calcul de la vitesse verticale a partir de la temperature
+! potentielle virtuelle
+
+! FH CESTQUOI CA ????
+#define int1d2
+!#undef int1d2
+#ifdef int1d2
+      if (l.ge.2) then
+#else
+      if (l.gt.2) then
+#endif
+
+      if (1.eq.1) then
+          w_est(ig,3)=zw2(ig,2)* &
+     &      ((f_star(ig,2))**2) &
+     &      /(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
+
+
+!---------------------------------------------------------------------------
+!calcul de l entrainement et du detrainement lateral
+!---------------------------------------------------------------------------
+!
+!test:estimation de ztva_new_est sans entrainement
+
+               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
+               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
+               qsatbef=MIN(0.5,qsatbef)
+               zcor=1./(1.-retv*qsatbef)
+               qsatbef=qsatbef*zcor
+               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
+               if (Zsat) then
+               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
+               DT = 0.5*RLvCp*qlbef
+               do while (abs(DT).gt.DDT0)
+                 Tbef=Tbef+DT
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
+                 qsatbef=MIN(0.5,qsatbef)
+                 zcor=1./(1.-retv*qsatbef)
+                 qsatbef=qsatbef*zcor
+                 qlbef=zqta(ig,l-1)-qsatbef
+
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
+                 zcor=1./(1.-retv*qsatbef)
+                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
+                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
+                 denom=1.+RLvCp*dqsat_dT
+                 DT=num/denom
+               enddo
+                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef) 
+               endif
+        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))
+
+             w_est(ig,l+1)=zw2(ig,l)*  &
+     &                   ((f_star(ig,l))**2)  &
+     &                   /(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
+                w_est(ig,l+1)=zw2(ig,l)
+             endif
+!
+!calcul du detrainement
+!=======================
+
+!CR:on vire les modifs
+         if (iflag_thermals_ed==0) then
+
+! Modifications du calcul du detrainement.
+! Dans la version de la these de Catherine, on passe brusquement
+! de la version seche a la version nuageuse pour le detrainement
+! ce qui peut occasioner des oscillations.
+! dans la nouvelle version, on commence par calculer un detrainement sec.
+! Puis un autre en cas de nuages.
+! Puis on combine les deux lineairement en fonction de la quantite d'eau.
+
+#define int1d3
+!#undef int1d3
+#define RIO_TH
+#ifdef RIO_TH
+!1. Cas non nuageux
+! 1.1 on est sous le zmax_sec et w croit
+          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
+     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
+#ifdef int1d3
+     &       (zqla_est(ig,l).lt.1.e-10)) then 
+#else
+     &       (zqla(ig,l-1).lt.1.e-10)) then 
+#endif
+             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
+     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
+     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
+     &       /(r_aspect*zmax_sec(ig)))
+             detr_stara(ig,l)=detr_star(ig,l)
+
+       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
+
+! 1.2 on est sous le zmax_sec et w decroit
+          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
+#ifdef int1d3
+     &            (zqla_est(ig,l).lt.1.e-10)) then
+#else
+     &            (zqla(ig,l-1).lt.1.e-10)) then
+#endif
+             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
+     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
+     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
+     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
+     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
+     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
+     &       *((zmax_sec(ig)-zlev(ig,l))/  &
+     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
+             detr_starb(ig,l)=detr_star(ig,l)
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
+
+          else
+
+! 1.3 dans les autres cas
+             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
+     &                      *(zlev(ig,l+1)-zlev(ig,l))
+             detr_starc(ig,l)=detr_star(ig,l)
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
+             
+          endif
+
+#else
+
+! 1.1 on est sous le zmax_sec et w croit
+          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
+     &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
+             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
+     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
+     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
+     &       /(r_aspect*zmax_sec(ig)))
+
+       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
+
+! 1.2 on est sous le zmax_sec et w decroit
+          else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
+             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
+     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
+     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
+     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
+     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
+     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
+     &       *((zmax_sec(ig)-zlev(ig,l))/  &
+     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
+       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
+
+          else
+             detr_star=0.
+          endif
+
+! 1.3 dans les autres cas
+          detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
+     &                      *(zlev(ig,l+1)-zlev(ig,l))
+
+          coefc=min(zqla(ig,l-1)/1.e-3,1.)
+          if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
+          coefc=1.
+! il semble qu'il soit important de baser le calcul sur
+! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
+          detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
+
+#endif
+
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
+!IM 730508 beg
+!        if(itap.GE.7200) THEN
+!         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
+!        endif
+!IM 730508 end
+         
+         zqla0(ig,l)=zqla_est(ig,l)
+         detr_star0(ig,l)=detr_star(ig,l)
+!IM 060508 beg
+!         if(detr_star(ig,l).GT.1.) THEN
+!          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
+!   &      ,detr_starc(ig,l),coefc
+!         endif
+!IM 060508 end
+!IM 160508 beg
+!IM 160508       IF (f0(ig).NE.0.) THEN
+           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
+!IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
+!IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
+!IM 160508       ELSE
+!IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
+!IM 160508       ENDIF
+!IM 160508 end
+!IM 060508 beg
+!        if(detr_star(ig,l).GT.1.) THEN
+!         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
+!   &     float(1)/f0(ig)
+!        endif
+!IM 060508 end
+        if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
+!
+!calcul de entr_star
+
+! #undef test2
+! #ifdef test2
+! La version test2 destabilise beaucoup le modele.
+! Il semble donc que ca aide d'avoir un entrainement important sous
+! le nuage.
+!         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
+!          entr_star(ig,l)=0.4*detr_star(ig,l)
+!         else
+!          entr_star(ig,l)=0.
+!         endif
+! #else
+!
+! Deplacement du calcul de entr_star pour eviter d'avoir aussi
+! entr_star > fstar.
+! Redeplacer suite a la transformation du cas detr>f
+! FH
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
+#define int1d
+!FH 070508 #define int1d4
+!#undef int1d4
+! L'option int1d4 correspond au choix dans le cas ou le detrainement
+! devient trop grand.
+
+#ifdef int1d
+
+#ifdef int1d4
+#else
+       detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
+!FH 070508 plus
+       detr_star(ig,l)=min(detr_star(ig,l),1.)
+#endif
+
+       entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
+#ifdef int1d4
+! Si le detrainement excede le flux en bas + l'entrainement, le thermique
+! doit disparaitre.
+       if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
+          detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
+          f_star(ig,l+1)=0.
+          linter(ig)=l+1
+          zw2(ig,l+1)=-1.e-10
+       endif
+#endif
+
+
+#else
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
+        if(l.gt.lalim(ig)) then
+         entr_star(ig,l)=0.4*detr_star(ig,l)
+        else
+
+! FH :
+! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
+! en haut de la couche d'alimentation.
+! A remettre en questoin a la premiere occasion mais ca peut aider a 
+! ecrire un code robuste.
+! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
+! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
+! d'alimentation, ce qui n'est pas forcement heureux.
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
+#undef pre_int1c
+#ifdef pre_int1c
+         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
+         detr_star(ig,l)=entr_star(ig,l)
+#else
+         entr_star(ig,l)=0.
+#endif
+
+        endif
+
+#endif
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
+        entr_star1(ig,l)=entr_star(ig,l)
+        detr_star1(ig,l)=detr_star(ig,l)
+!
+
+#ifdef int1d
+#else
+        if (detr_star(ig,l).gt.f_star(ig,l)) then
+
+!  Ce test est là entre autres parce qu'on passe par des valeurs
+!  delirantes de detr_star.
+!  ca vaut sans doute le coup de verifier pourquoi.
+
+           detr_star(ig,l)=f_star(ig,l)
+#ifdef pre_int1c
+           if (l.gt.lalim(ig)+1) then
+               entr_star(ig,l)=0.
+               alim_star(ig,l)=0.
+! FH ajout pour forcer a stoper le thermique juste sous le sommet
+! de la couche (voir calcul de finter)
+               zw2(ig,l+1)=-1.e-10
+               linter(ig)=l+1
+            else
+               entr_star(ig,l)=0.4*detr_star(ig,l)
+            endif
+#else
+           entr_star(ig,l)=0.4*detr_star(ig,l)
+#endif
+        endif
+#endif
+
+      else !l > 2
+         detr_star(ig,l)=0.
+         entr_star(ig,l)=0.
+      endif
+
+        entr_star2(ig,l)=entr_star(ig,l)
+        detr_star2(ig,l)=detr_star(ig,l)
+        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)
+            if(prt_level.GE.10) 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
+! dans la couche d'alimentation
+!pas d entrainement dans la couche alim
+!      if ((l.le.lalim(ig))) then
+!           entr_star(ig,l)=0.
+!      endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!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)
+
+!test sur le signe de f_star
+        if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
+       if (f_star(ig,l+1).gt.1.e-10) then 
+!----------------------------------------------------------------------------
+!calcul de la vitesse verticale en melangeant Tl et qt du thermique
+!---------------------------------------------------------------------------
+!
+       Zsat=.false.
+       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
+     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
+     &            /(f_star(ig,l+1)+detr_star(ig,l))
+!
+       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
+     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
+     &            /(f_star(ig,l+1)+detr_star(ig,l))
+!  
+               Tbef=ztla(ig,l)*zpspsk(ig,l)
+               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
+               qsatbef=MIN(0.5,qsatbef)
+               zcor=1./(1.-retv*qsatbef)
+               qsatbef=qsatbef*zcor
+               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
+               if (Zsat) then
+               qlbef=max(0.,zqta(ig,l)-qsatbef)
+               DT = 0.5*RLvCp*qlbef
+               do while (abs(DT).gt.DDT0)
+                 Tbef=Tbef+DT
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
+                 qsatbef=MIN(0.5,qsatbef)
+                 zcor=1./(1.-retv*qsatbef)
+                 qsatbef=qsatbef*zcor
+                 qlbef=zqta(ig,l)-qsatbef
+
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
+                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
+                 zcor=1./(1.-retv*qsatbef)
+                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
+                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
+                 denom=1.+RLvCp*dqsat_dT
+                 DT=num/denom
+              enddo
+                 zqla(ig,l) = max(0.,qlbef) 
+              endif
+!    
+        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
+! on ecrit de maniere conservative (sat ou non)
+!          T = Tl +Lv/Cp ql
+           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
+           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
+!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
+           zha(ig,l) = ztva(ig,l)
+           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
+     &              -zqla(ig,l))-zqla(ig,l))
+
+!on ecrit zqsat 
+           zqsatth(ig,l)=qsatbef  
+!calcul de vitesse
+           zw2(ig,l+1)=zw2(ig,l)*  &
+     &                 ((f_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
+        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
+!
+!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
+
+            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'
+                print*,'On tombe sur le cas particulier de thermcell_plume'
+                zw2(ig,l+1)=0.
+                linter(ig)=l+1
+            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))  &
+     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
+           zw2(ig,l+1)=0.
+        endif
+
+           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 
+
+        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)
+        endif
+        enddo
+      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
+
+!     print*,'thermcell_plume OK'
+
+      return 
+      end
+!==============================================================================
+       SUBROUTINE thermcellV0_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
+     &                            lalim,lmin,zmax,wmax,lev_out)
+
+!--------------------------------------------------------------------------
+!thermcell_dry: calcul de zmax et wmax du thermique sec
+!--------------------------------------------------------------------------
+       IMPLICIT NONE
+#include "YOMCST.h"       
+#include "iniprint.h"
+       INTEGER l,ig
+
+       INTEGER ngrid,nlay
+       REAL zlev(ngrid,nlay+1)
+       REAL pphi(ngrid,nlay)
+       REAl ztv(ngrid,nlay)
+       REAL alim_star(ngrid,nlay)
+       INTEGER lalim(ngrid)
+      integer lev_out                           ! niveau pour les print
+
+       REAL zmax(ngrid)
+       REAL wmax(ngrid)
+
+!variables locales
+       REAL zw2(ngrid,nlay+1)
+       REAL f_star(ngrid,nlay+1)
+       REAL ztva(ngrid,nlay+1)
+       REAL wmaxa(ngrid)
+       REAL wa_moy(ngrid,nlay+1)
+       REAL linter(ngrid),zlevinter(ngrid)
+       INTEGER lmix(ngrid),lmax(ngrid),lmin(ngrid)
+
+!initialisations
+       do ig=1,ngrid
+          do l=1,nlay+1
+             zw2(ig,l)=0.
+             wa_moy(ig,l)=0.
+          enddo
+       enddo
+       do ig=1,ngrid
+          do l=1,nlay
+             ztva(ig,l)=ztv(ig,l)
+          enddo
+       enddo
+       do ig=1,ngrid
+          wmax(ig)=0.
+          wmaxa(ig)=0.
+       enddo
+!calcul de la vitesse a partir de la CAPE en melangeant thetav
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! A eliminer
+! Ce if complique etait fait pour reperer la premiere couche instable
+! Ici, c'est lmin.
+!
+!       do l=1,nlay-2
+!         do ig=1,ngrid
+!            if (ztv(ig,l).gt.ztv(ig,l+1)  &
+!     &         .and.alim_star(ig,l).gt.1.e-10  &
+!     &         .and.zw2(ig,l).lt.1e-10) then
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+! Calcul des F^*, integrale verticale de E^*
+       f_star(:,1)=0.
+       do l=1,nlay
+          f_star(:,l+1)=f_star(:,l)+alim_star(:,l)
+       enddo
+
+! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
+       linter(:)=0.
+
+! couche la plus haute concernee par le thermique. 
+       lmax(:)=1
+
+! Le niveau linter est une variable continue qui se trouve dans la couche
+! lmax
+
+       do l=1,nlay-2
+         do ig=1,ngrid
+            if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
+
+!------------------------------------------------------------------------
+!  Calcul de la vitesse en haut de la premiere couche instable.
+!  Premiere couche du panache thermique
+!------------------------------------------------------------------------
+               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
+     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
+     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
+
+!------------------------------------------------------------------------
+! Tant que la vitesse en bas de la couche et la somme du flux de masse
+! et de l'entrainement (c'est a dire le flux de masse en haut) sont
+! positifs, on calcul
+! 1. le flux de masse en haut  f_star(ig,l+1)
+! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
+! 3. la vitesse au carré en haut zw2(ig,l+1)
+!------------------------------------------------------------------------
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  A eliminer : dans cette version, si zw2 est > 0 on a un therique.
+!  et donc, au dessus, f_star(ig,l+1) est forcement suffisamment 
+!  grand puisque on n'a pas de detrainement.
+!  f_star est une fonction croissante.
+!  c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
+!           else if ((zw2(ig,l).ge.1e-10).and.  &
+!    &               (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
+!              f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+            else if (zw2(ig,l).ge.1e-10) then
+
+               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l)  &
+     &                    *ztv(ig,l))/f_star(ig,l+1)
+               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+  &
+     &                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+     &                     *(zlev(ig,l+1)-zlev(ig,l))
+            endif
+! determination de zmax continu par interpolation lineaire
+!------------------------------------------------------------------------
+
+            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'
+!               print*,'On tombe sur le cas particulier de thermcell_dry'
+                zw2(ig,l+1)=0.
+                linter(ig)=l+1
+                lmax(ig)=l
+            endif
+
+            if (zw2(ig,l+1).lt.0.) then
+               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
+     &           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
+               zw2(ig,l+1)=0.
+               lmax(ig)=l
+            endif
+
+               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
+
+            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
+!   lmix est le niveau de la couche ou w (wa_moy) est maximum
+               lmix(ig)=l+1
+               wmaxa(ig)=wa_moy(ig,l+1)
+            endif
+         enddo
+      enddo
+       if (prt_level.ge.1) print*,'fin calcul zw2'
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! A eliminer :
+! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
+! Calcul de la couche correspondant a la hauteur du thermique
+!      do ig=1,ngrid
+!         lmax(ig)=lalim(ig)
+!      enddo
+!      do ig=1,ngrid
+!         do l=nlay,lalim(ig)+1,-1
+!            if (zw2(ig,l).le.1.e-10) then
+!               lmax(ig)=l-1
+!            endif
+!         enddo
+!      enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!    
+! Determination de zw2 max
+      do ig=1,ngrid
+         wmax(ig)=0.
+      enddo
+
+      do l=1,nlay
+         do ig=1,ngrid
+            if (l.le.lmax(ig)) then
+                zw2(ig,l)=sqrt(zw2(ig,l))
+                wmax(ig)=max(wmax(ig),zw2(ig,l))
+            else
+                 zw2(ig,l)=0.
+            endif
+          enddo
+      enddo
+
+!   Longueur caracteristique correspondant a la hauteur des thermiques.
+      do  ig=1,ngrid
+         zmax(ig)=0.
+         zlevinter(ig)=zlev(ig,1)
+      enddo
+      do  ig=1,ngrid
+! calcul de zlevinter
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH A eliminer
+! Simplification
+!          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
+!     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
+!     &    -zlev(ig,lmax(ig)))
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+          zlevinter(ig)=zlev(ig,lmax(ig)) + &
+     &    (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
+           zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
+      enddo
+
+! Verification que lalim<=lmax
+      do ig=1,ngrid
+         if(lalim(ig)>lmax(ig)) then
+           if ( prt_level > 1 ) THEN
+            print*,'WARNING thermcell_dry ig=',ig,'  lalim=',lalim(ig),'  lmax(ig)=',lmax(ig)
+           endif
+           lmax(ig)=lalim(ig)
+         endif
+      enddo
+      
+      RETURN
+      END
+!==============================================================================
+      SUBROUTINE thermcellV0_init(ngrid,nlay,ztv,zlay,zlev,  &
+     &                  lalim,lmin,alim_star,alim_star_tot,lev_out)
+
+!----------------------------------------------------------------------
+!thermcell_init: calcul du profil d alimentation du thermique
+!----------------------------------------------------------------------
+      IMPLICIT NONE
+#include "iniprint.h"
+#include "thermcell.h"
+
+      INTEGER l,ig
+!arguments d entree
+      INTEGER ngrid,nlay
+      REAL ztv(ngrid,nlay)
+      REAL zlay(ngrid,nlay)
+      REAL zlev(ngrid,nlay+1)
+!arguments de sortie
+      INTEGER lalim(ngrid)
+      INTEGER lmin(ngrid)
+      REAL alim_star(ngrid,nlay)
+      REAL alim_star_tot(ngrid)
+      integer lev_out                           ! niveau pour les print
+      
+      REAL zzalim(ngrid)
+!CR: ponderation entrainement des couches instables
+!def des alim_star tels que alim=f*alim_star      
+
+      do l=1,nlay
+         do ig=1,ngrid 
+            alim_star(ig,l)=0.
+         enddo
+      enddo
+! determination de la longueur de la couche d entrainement
+      do ig=1,ngrid
+         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)
+               if(prt_level.GE.10) 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
+      do l=nlay-2,1,-1
+         do ig=1,ngrid
+            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
+     &          ztv(ig,l+1).le.ztv(ig,l+2)) then
+               lalim(ig)=l+1
+            endif
+          enddo
+      enddo
+
+! determination du lmin: couche d ou provient le thermique
+
+      do ig=1,ngrid
+! FH initialisation de lmin a nlay plutot que 1.
+!        lmin(ig)=nlay
+         lmin(ig)=1
+      enddo
+      do l=nlay,2,-1
+         do ig=1,ngrid
+            if (ztv(ig,l-1).gt.ztv(ig,l)) then
+               lmin(ig)=l-1
+            endif
+         enddo
+      enddo
+!
+      zzalim(:)=0.
+      do l=1,nlay-1
+         do ig=1,ngrid 
+             if (l<lalim(ig)) then
+                zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
+             endif
+          enddo
+      enddo
+      do ig=1,ngrid
+          if (lalim(ig)>1) then
+             zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
+          else
+             zzalim(ig)=zlay(ig,1)
+          endif
+      enddo
+
+      if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))
+
+! definition de l'entrainement des couches
+      if (1.eq.1) then
+      do l=1,nlay-1
+         do ig=1,ngrid 
+            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
+     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
+!def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
+             alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
+     &                       *sqrt(zlev(ig,l+1)) 
+            endif
+         enddo
+      enddo
+      else
+      do l=1,nlay-1
+         do ig=1,ngrid
+            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
+     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
+             alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
+     &        *(zlev(ig,l+1)-zlev(ig,l))
+            endif
+         enddo
+      enddo
+      endif
+      
+! pas de thermique si couche 1 stable
+      do ig=1,ngrid
+!CRnouveau test
+        if (alim_star(ig,1).lt.1.e-10) then 
+            do l=1,nlay
+                alim_star(ig,l)=0.
+            enddo
+            lmin(ig)=1
+         endif
+      enddo 
+! calcul de l alimentation totale
+      do ig=1,ngrid
+         alim_star_tot(ig)=0.
+      enddo
+      do l=1,nlay
+         do ig=1,ngrid
+            alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
+         enddo
+      enddo
+!
+! Calcul entrainement normalise
+      do l=1,nlay 
+         do ig=1,ngrid 
+            if (alim_star_tot(ig).gt.1.e-10) then
+               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
+            endif
+         enddo
+      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/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_closure.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_closure.F90	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_closure.F90	(revision 1294)
@@ -1,7 +1,17 @@
+!
+! $Header$
+!
       SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
+     &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
 
 !-------------------------------------------------------------------------
 !thermcell_closure: fermeture, determination de f
+!
+! Modification 7 septembre 2009
+! 1. On enleve alim_star_tot des arguments pour le recalculer et etre ainis
+! coherent avec l'integrale au numerateur.
+! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec
+! l'idee etant que le choix se fasse a l'appel de thermcell_closure
+! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
 !-------------------------------------------------------------------------
       IMPLICIT NONE
@@ -9,63 +19,64 @@
 #include "iniprint.h"
 #include "thermcell.h"
-      INTEGER ngrid,nlay
-      INTEGER ig,k       
-      REAL r_aspect,ptimestep
-      integer lev_out                           ! niveau pour les print
+INTEGER ngrid,nlay
+INTEGER ig,k       
+REAL r_aspect,ptimestep
+integer lev_out                           ! niveau pour les print
 
-      INTEGER lalim(ngrid)
-      REAL alim_star(ngrid,nlay)
-      REAL alim_star_tot(ngrid)
-      REAL rho(ngrid,nlay)
-      REAL zlev(ngrid,nlay)
-      REAL zmax(ngrid),zmax_sec(ngrid)
-      REAL wmax(ngrid),wmax_sec(ngrid)
-      real zdenom
+INTEGER lalim(ngrid)
+REAL alim_star(ngrid,nlay)
+REAL f_star(ngrid,nlay+1)
+REAL rho(ngrid,nlay)
+REAL zlev(ngrid,nlay)
+REAL zmax(ngrid)
+REAL wmax(ngrid)
+REAL zdenom(ngrid)
+REAL alim_star2(ngrid)
+REAL f(ngrid)
 
-      REAL alim_star2(ngrid)
+REAL alim_star_tot(ngrid)
+INTEGER llmax
 
-      REAL f(ngrid)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+print*,'THERMCELL CLOSURE 26E'
 
-      do ig=1,ngrid
-         alim_star2(ig)=0.
-      enddo
-      do ig=1,ngrid
-         if (alim_star(ig,1).LT.1.e-10) then
-            f(ig)=0.
-         else   
-             do k=1,lalim(ig)
-                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
-     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
-             enddo
-             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
-             if (zdenom<1.e-14) then
-                print*,'ig=',ig
-                print*,'alim_star2',alim_star2(ig)
-                print*,'zmax',zmax(ig)
-                print*,'r_aspect',r_aspect
-                print*,'zdenom',zdenom
-                print*,'alim_star',alim_star(ig,:)
-                print*,'zmax_sec',zmax_sec(ig)
-                print*,'wmax_sec',wmax_sec(ig)
-                stop
-             endif
-             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))
-             if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
-             else
-             f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
-!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
-!     &                     zmax(ig))*wmax(ig))
-             if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
-             endif
-         endif
-!         f0(ig)=f(ig)
-      enddo
-      if (prt_level.ge.1) print*,'apres fermeture'
+alim_star2(:)=0.
+alim_star_tot(:)=0.
+f(:)=0.
 
-!
+! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
+llmax=1
+do ig=1,ngrid
+   if (lalim(ig)>llmax) llmax=lalim(ig)
+enddo
+
+
+! Calcul des integrales sur la verticale de alim_star et de
+!   alim_star^2/(rho dz)
+do k=1,llmax-1
+   do ig=1,ngrid
+      if (k<lalim(ig)) then
+         alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
+&                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
+         alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,k)
+      endif
+   enddo
+enddo
+
+
+do ig=1,ngrid
+   if (alim_star2(ig)>1.e-10) then
+      f(ig)=wmax(ig)*alim_star_tot(ig)/  &
+&     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
+   endif
+enddo
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! TESTS POUR UNE NOUVELLE FERMETURE DANS LAQUELLE ALIM_STAR NE SERAIT
+! PAS NORMALISE
+!           f(ig)=f(ig)*f_star(ig,2)/(f_star(ig,lalim(ig)))
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
       return
       end
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_condens.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_condens.F90	(revision 1294)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_condens.F90	(revision 1294)
@@ -0,0 +1,81 @@
+subroutine thermcell_condens(klon,active,zpspsk,pplev,ztla,zqta,zqla)
+implicit none
+
+#include "YOMCST.h"
+#include "YOETHF.h"
+#include "FCTTRE.h"
+
+
+!====================================================================
+! DECLARATIONS
+!====================================================================
+
+! Arguments
+INTEGER klon
+REAL zpspsk(klon),pplev(klon)
+REAL ztla(klon),zqta(klon),zqla(klon)
+LOGICAL active(klon)
+
+! Variables locales
+INTEGER ig,iter
+REAL Tbef(klon),DT(klon)
+REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dT
+logical Zsat
+REAL RLvCp
+REAL, SAVE :: DDT0=.01
+LOGICAL afaire(klon),tout_converge
+
+!====================================================================
+! INITIALISATIONS
+!====================================================================
+
+RLvCp = RLVTT/RCPD
+tout_converge=.false.
+afaire(:)=.false.
+DT(:)=0.
+
+
+!====================================================================
+! Routine a vectoriser en copiant active dans converge et en mettant
+! la boucle sur les iterations a l'exterieur est en mettant
+! converge= false des que la convergence est atteinte.
+!====================================================================
+
+do ig=1,klon
+   if (active(ig)) then
+               Tbef(ig)=ztla(ig)*zpspsk(ig)
+               zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
+               qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
+               qsatbef=MIN(0.5,qsatbef)
+               zcor=1./(1.-retv*qsatbef)
+               qsatbef=qsatbef*zcor
+               qlbef=max(0.,zqta(ig)-qsatbef)
+               DT(ig) = 0.5*RLvCp*qlbef
+     endif
+enddo
+
+do iter=1,10
+    afaire(:)=abs(DT(:)).gt.DDT0
+    do ig=1,klon
+               if (afaire(ig)) then
+                 Tbef(ig)=Tbef(ig)+DT(ig)
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
+                 qsatbef= R2ES * FOEEW(Tbef(ig),zdelta)/pplev(ig)
+                 qsatbef=MIN(0.5,qsatbef)
+                 zcor=1./(1.-retv*qsatbef)
+                 qsatbef=qsatbef*zcor
+                 qlbef=zqta(ig)-qsatbef
+                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef(ig)))
+                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
+                 zcor=1./(1.-retv*qsatbef)
+                 dqsat_dT=FOEDE(Tbef(ig),zdelta,zcvm5,qsatbef,zcor)
+                 num=-Tbef(ig)+ztla(ig)*zpspsk(ig)+RLvCp*qlbef
+                 denom=1.+RLvCp*dqsat_dT
+                 zqla(ig) = max(0.,zqta(ig)-qsatbef) 
+                 DT(ig)=num/denom
+               endif
+    enddo
+enddo
+
+return
+end
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dry.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dry.F90	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dry.F90	(revision 1294)
@@ -1,2 +1,5 @@
+!
+! $Header$
+!
        SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
      &                            lalim,lmin,zmax,wmax,lev_out)
@@ -4,5 +7,12 @@
 !--------------------------------------------------------------------------
 !thermcell_dry: calcul de zmax et wmax du thermique sec
+! Calcul de la vitesse maximum et de la hauteur maximum pour un panache
+! ascendant avec une fonction d'alimentation alim_star et sans changement 
+! de phase.
+! Le calcul pourrait etre sans doute simplifier.
+! La temperature potentielle virtuelle dans la panache ascendant est
+! la temperature potentielle virtuelle pondÃ©rÃ©e par alim_star.
 !--------------------------------------------------------------------------
+
        IMPLICIT NONE
 #include "YOMCST.h"       
@@ -48,16 +58,4 @@
 !calcul de la vitesse a partir de la CAPE en melangeant thetav
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! A eliminer
-! Ce if complique etait fait pour reperer la premiere couche instable
-! Ici, c'est lmin.
-!
-!       do l=1,nlay-2
-!         do ig=1,ngrid
-!            if (ztv(ig,l).gt.ztv(ig,l+1)  &
-!     &         .and.alim_star(ig,l).gt.1.e-10  &
-!     &         .and.zw2(ig,l).lt.1e-10) then
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
 
 ! Calcul des F^*, integrale verticale de E^*
@@ -84,4 +82,5 @@
 !  Premiere couche du panache thermique
 !------------------------------------------------------------------------
+
                zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
      &                     *(zlev(ig,l+1)-zlev(ig,l))  &
@@ -96,15 +95,4 @@
 ! 3. la vitesse au carré en haut zw2(ig,l+1)
 !------------------------------------------------------------------------
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!  A eliminer : dans cette version, si zw2 est > 0 on a un therique.
-!  et donc, au dessus, f_star(ig,l+1) est forcement suffisamment 
-!  grand puisque on n'a pas de detrainement.
-!  f_star est une fonction croissante.
-!  c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
-!           else if ((zw2(ig,l).ge.1e-10).and.  &
-!    &               (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
-!              f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
             else if (zw2(ig,l).ge.1e-10) then
@@ -145,21 +133,4 @@
        if (prt_level.ge.1) print*,'fin calcul zw2'
 !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! A eliminer :
-! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
-! Calcul de la couche correspondant a la hauteur du thermique
-!      do ig=1,ngrid
-!         lmax(ig)=lalim(ig)
-!      enddo
-!      do ig=1,ngrid
-!         do l=nlay,lalim(ig)+1,-1
-!            if (zw2(ig,l).le.1.e-10) then
-!               lmax(ig)=l-1
-!            endif
-!         enddo
-!      enddo
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!    
 ! Determination de zw2 max
       do ig=1,ngrid
@@ -185,13 +156,4 @@
       do  ig=1,ngrid
 ! calcul de zlevinter
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH A eliminer
-! Simplification
-!          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
-!     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
-!     &    -zlev(ig,lmax(ig)))
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
           zlevinter(ig)=zlev(ig,lmax(ig)) + &
      &    (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
@@ -199,14 +161,4 @@
       enddo
 
-! Verification que lalim<=lmax
-      do ig=1,ngrid
-         if(lalim(ig)>lmax(ig)) then
-           if ( prt_level > 1 ) THEN
-            print*,'WARNING thermcell_dry ig=',ig,'  lalim=',lalim(ig),'  lmax(ig)=',lmax(ig)
-           endif
-           lmax(ig)=lalim(ig)
-         endif
-      enddo
-      
       RETURN
       END
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_init.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_init.F90	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_init.F90	(revision 1294)
@@ -1,2 +1,5 @@
+!
+! $Header$
+!
       SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
      &                  lalim,lmin,alim_star,alim_star_tot,lev_out)
@@ -26,137 +29,31 @@
 !def des alim_star tels que alim=f*alim_star      
 
-      do l=1,nlay
-         do ig=1,ngrid 
-            alim_star(ig,l)=0.
-         enddo
-      enddo
-! determination de la longueur de la couche d entrainement
-      do ig=1,ngrid
-         lalim(ig)=1
-      enddo
 
-      if (iflag_thermals_ed.ge.1) then
-!si la première couche est instable, on declenche un thermique
+      write(lunout,*)'THERM INIT V20C '
+
+      alim_star_tot(:)=0.
+      alim_star(:,:)=0.
+      lmin(:)=1
+      lalim(:)=1
+
+      do l=1,nlay-1
          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)
-               if(prt_level.GE.10) 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
-      do l=nlay-2,1,-1
-         do ig=1,ngrid
-            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
-     &          ztv(ig,l+1).le.ztv(ig,l+2)) then
-               lalim(ig)=l+1
-            endif
-          enddo
-      enddo
-
-! determination du lmin: couche d ou provient le thermique
-
-      do ig=1,ngrid
-! FH initialisation de lmin a nlay plutot que 1.
-!        lmin(ig)=nlay
-         lmin(ig)=1
-      enddo
-      do l=nlay,2,-1
-         do ig=1,ngrid
-            if (ztv(ig,l-1).gt.ztv(ig,l)) then
-               lmin(ig)=l-1
+            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
+               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
+     &                       *sqrt(zlev(ig,l+1)) 
+               lalim(:)=l+1
+               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
             endif
          enddo
       enddo
-!
-      zzalim(:)=0.
-      do l=1,nlay-1
+      do l=1,nlay
          do ig=1,ngrid 
-             if (l<lalim(ig)) then
-                zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
-             endif
-          enddo
-      enddo
-      do ig=1,ngrid
-          if (lalim(ig)>1) then
-             zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
-          else
-             zzalim(ig)=zlay(ig,1)
-          endif
-      enddo
-
-      if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))
-
-! definition de l'entrainement des couches
-      if (1.eq.1) then
-      do l=1,nlay-1
-         do ig=1,ngrid 
-            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
-     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
-!def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
-             alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
-     &                       *sqrt(zlev(ig,l+1)) 
-            endif
-         enddo
-      enddo
-      else
-      do l=1,nlay-1
-         do ig=1,ngrid
-            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
-     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
-             alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
-     &        *(zlev(ig,l+1)-zlev(ig,l))
-            endif
-         enddo
-      enddo
-      endif
-      
-! pas de thermique si couche 1 stable
-      do ig=1,ngrid
-!CRnouveau test
-        if (alim_star(ig,1).lt.1.e-10) then 
-            do l=1,nlay
-                alim_star(ig,l)=0.
-            enddo
-            lmin(ig)=1
-         endif
-      enddo 
-! calcul de l alimentation totale
-      do ig=1,ngrid
-         alim_star_tot(ig)=0.
-      enddo
-      do l=1,nlay
-         do ig=1,ngrid
-            alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
-         enddo
-      enddo
-!
-! Calcul entrainement normalise
-      do l=1,nlay 
-         do ig=1,ngrid 
-            if (alim_star_tot(ig).gt.1.e-10) then
+            if (alim_star_tot(ig) > 1.e-10 ) then
                alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
             endif
          enddo
       enddo
-       
-!on remet alim_star_tot a 1
-      do ig=1,ngrid 
-         alim_star_tot(ig)=1.
-      enddo
+      alim_star_tot(:)=1.
 
-      endif
-!endif iflag_thermals_ed
-      return 
+      return
       end  
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90	(revision 1294)
@@ -22,11 +22,11 @@
 !   de "thermiques" explicitement representes avec processus nuageux
 !
-!   Réécriture à partir d'un listing papier à Habas, le 14/02/00
-!
-!   le thermique est supposé homogène et dissipé par mélange avec
-!   son environnement. la longueur l_mix contrôle l'efficacité du
-!   mélange
-!
-!   Le calcul du transport des différentes espèces se fait en prenant
+!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
+!
+!   le thermique est suppose homogene et dissipe par melange avec
+!   son environnement. la longueur l_mix controle l'efficacite du
+!   melange
+!
+!   Le calcul du transport des differentes especes se fait en prenant
 !   en compte:
 !     1. un flux de masse montant
@@ -85,5 +85,5 @@
       real linter(klon)
       real zmix(klon)
-      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1)
+      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
 !      real fraca(klon,klev)
 
@@ -115,4 +115,7 @@
 ! FH probleme de dimensionnement avec l'allocation dynamique
 !     common/comtherm/thetath2,wth2
+      real wq(klon,klev)
+      real wthl(klon,klev)
+      real wthv(klon,klev)
     
       real ratqscth(klon,klev)
@@ -142,6 +145,7 @@
       real f_star(klon,klev+1),entr_star(klon,klev)
       real detr_star(klon,klev)
-      real alim_star_tot(klon),alim_star2(klon)
+      real alim_star_tot(klon)
       real alim_star(klon,klev)
+      real alim_star_clos(klon,klev)
       real f(klon), f0(klon)
 !FH/IM     save f0
@@ -149,4 +153,5 @@
       logical debut
        real seuil
+      real csc(klon,klev)
 
 !
@@ -220,18 +225,5 @@
       ENDIF
 !
-!Initialisation
-!
-!    IF (1.eq.0) THEN
-!     do ig=1,klon      
-!FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
-!     if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
-!           f0(ig)=1.e-5
-!           zmax0(ig)=40.
-!v1d        therm=.false.
-!     endif
-!     enddo 
-!    ENDIF !(1.eq.0) THEN
-     if (prt_level.ge.10)write(lunout,*)                                &
-    &     'WARNING thermcell_main f0=max(f0,1.e-2)'
+     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
      do ig=1,klon
       if (prt_level.ge.20) then
@@ -239,4 +231,5 @@
       endif
          f0(ig)=max(f0(ig),1.e-2)
+         zmax0(ig)=max(zmax0(ig),40.)
 !IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
      enddo
@@ -364,41 +357,68 @@
 
 !------------------------------------------------------------------
-!  1. alim_star est le profil vertical de l'alimentation à la base du
-!     panache thermique, calculé à partir de la flotabilité de l'air sec
+!  1. alim_star est le profil vertical de l'alimentation a la base du
+!     panache thermique, calcule a partir de la flotabilite de l'air sec
 !  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
 !------------------------------------------------------------------
 !
       entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
-      CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
-     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
-
-call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
-call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
-
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
-      if (prt_level.ge.10) then
-         write(lunout1,*) 'Dans thermcell_main 1'
-         write(lunout1,*) 'lmin ',lmin(igout)
-         write(lunout1,*) 'lalim ',lalim(igout)
-         write(lunout1,*) ' ig l alim_star thetav'
-         write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
-     &   ,ztv(igout,l),l=1,lalim(igout)+4)
-      endif
-
-!v1d      do ig=1,klon
-!v1d     if (alim_star(ig,1).gt.1.e-10) then
-!v1d     therm=.true.
-!v1d     endif
-!v1d      enddo
+      lmin=1
+
 !-----------------------------------------------------------------------------
 !  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
 !     panache sec conservatif (e=d=0) alimente selon alim_star 
 !     Il s'agit d'un calcul de type CAPE
-!     zmax_sec est utilisé pour déterminer la géométrie du thermique.
+!     zmax_sec est utilise pour determiner la geometrie du thermique.
 !------------------------------------------------------------------------------
-!
+!---------------------------------------------------------------------------------
+!calcul du melange et des variables dans le thermique
+!--------------------------------------------------------------------------------
+!
+      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
+!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,alim_star,alim_star_tot,  &
+     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
+     &           ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
+     &            ,lev_out,lunout1,igout)
+      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
+
+      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
+      if (prt_level.ge.10) then
+         write(lunout1,*) 'Dans thermcell_main 2'
+         write(lunout1,*) 'lmin ',lmin(igout)
+         write(lunout1,*) 'lalim ',lalim(igout)
+         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
+         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
+     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
+      endif
+
+!-------------------------------------------------------------------------------
+! Calcul des caracteristiques du thermique:zmax,zmix,wmax
+!-------------------------------------------------------------------------------
+!
+      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
+     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
+
+
+
+      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
+      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
+
+      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
+
+!-------------------------------------------------------------------------------
+! Fermeture,determination de f
+!-------------------------------------------------------------------------------
+!
+!
+      write(lunout,*)'THERM NOUVEAU RIO2009 31B'
       CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
-     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
+    &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
 
 call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
@@ -417,59 +437,21 @@
 
 
-!---------------------------------------------------------------------------------
-!calcul du melange et des variables dans le thermique
-!--------------------------------------------------------------------------------
-!
-      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
-!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,alim_star_tot,  &
-     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
-     &           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
-
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
-      if (prt_level.ge.10) then
-         write(lunout1,*) 'Dans thermcell_main 2'
-         write(lunout1,*) 'lmin ',lmin(igout)
-         write(lunout1,*) 'lalim ',lalim(igout)
-         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
-         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
-     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
-      endif
-
-!-------------------------------------------------------------------------------
-! Calcul des caracteristiques du thermique:zmax,zmix,wmax
-!-------------------------------------------------------------------------------
-!
-      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
-     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
-
-
-      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
-      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
-
-      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
-
-!-------------------------------------------------------------------------------
-! 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
-
+print*,'THERM 26JJJ'
+
+! Choix de la fonction d'alimentation utilisee pour la fermeture.
+! Apparemment sans importance
+      alim_star_clos(:,:)=alim_star(:,:)
+      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
+
+! Appel avec la version seche
       CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
+     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Appel avec les zmax et wmax tenant compte de la condensation
+! Semble moins bien marcher
+!     CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
+!    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
@@ -484,5 +466,5 @@
 ! Test valable seulement en 1D mais pas genant
       if (.not. (f0(1).ge.0.) ) then
-           stop 'Dans thermcell_main'
+           stop'Dans thermcell_main'
       endif
 
@@ -511,5 +493,5 @@
          fm0=(1.-lambda)*fm+lambda*fm0
          entr0=(1.-lambda)*entr+lambda*entr0
-!        detr0=(1.-lambda)*detr+lambda*detr0
+         detr0=(1.-lambda)*detr+lambda*detr0
       else
          fm0=fm
@@ -560,5 +542,5 @@
      &    ,fraca,zmax  &
      &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
-!IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
+
       else
 
@@ -636,5 +618,5 @@
 !
       if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
-            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
+            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
             if(zw2(ig,l).gt.1.e-10) then
              wth2(ig,l)=zf2*(zw2(ig,l))**2
@@ -651,6 +633,14 @@
          enddo
       enddo
+!calcul des flux: q, thetal et thetav
+      do l=1,nlay
+         do ig=1,ngrid
+      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
+      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
+      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
+         enddo
+      enddo
 !calcul de ale_bl et alp_bl
-!pour le calcul d'une valeur intégrée entre la surface et lmax
+!pour le calcul d'une valeur integree entre la surface et lmax
       do ig=1,ngrid
       alp_int(ig)=0.
@@ -782,6 +772,6 @@
 !     print*,'15 OK convect8'
 
+#ifdef wrgrads_thermcell
       if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
-#ifdef wrgrads_thermcell
 #include "thermcell_out3d.h"
 #endif
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_out3d.h
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_out3d.h	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_out3d.h	(revision 1294)
@@ -27,4 +27,9 @@
          call wrgradsfi(1,nlay,q2(igout,1:klev),'q2       ','q2       ')
 !
+!
+         call wrgradsfi(1,nlay,wthl(igout,1:klev),'wthl       ','wthl       ')
+         call wrgradsfi(1,nlay,wthv(igout,1:klev),'wthv       ','wthv       ')
+         call wrgradsfi(1,nlay,wq(igout,1:klev),'wq       ','wq       ')
+         
          call wrgradsfi(1,nlay,ztva(igout,1:klev),'ztva      ','ztva      ')
          call wrgradsfi(1,nlay,ztv(igout,1:klev),'ztv       ','ztv       ')
@@ -53,5 +58,4 @@
       call wrgradsfi(1,1,f(igout),'f      ','f      ')
       call wrgradsfi(1,1,alim_star_tot(igout),'a_s_t      ','a_s_t      ')
-      call wrgradsfi(1,1,alim_star2(igout),'a_2      ','a_2      ')
       call wrgradsfi(1,1,zmax(igout),'zmax      ','zmax      ')
       call wrgradsfi(1,1,zmax_sec(igout),'z_sec      ','z_sec      ')
Index: /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90
===================================================================
--- /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90	(revision 1293)
+++ /LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90	(revision 1294)
@@ -1,6 +1,9 @@
+!
+! $Header$
+!
       SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
-     &           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,w_est,zqsatth,lmix,lmix_bis,linter &
+     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
+     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
+     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
      &           ,lev_out,lunout1,igout)
 
@@ -31,8 +34,5 @@
       REAL zpspsk(ngrid,klev)
       REAL alim_star(ngrid,klev)
-      REAL zmax_sec(ngrid)
       REAL f0(ngrid)
-      REAL l_mix
-      REAL r_aspect
       INTEGER lalim(ngrid)
       integer lev_out                           ! niveau pour les print
@@ -44,5 +44,4 @@
       REAL ztla(ngrid,klev)
       REAL zqla(ngrid,klev)
-      REAL zqla0(ngrid,klev)
       REAL zqta(ngrid,klev)
       REAL zha(ngrid,klev)
@@ -50,16 +49,9 @@
       REAL detr_star(ngrid,klev)
       REAL coefc
-      REAL detr_stara(ngrid,klev)
-      REAL detr_starb(ngrid,klev)
-      REAL detr_starc(ngrid,klev)
-      REAL detr_star0(ngrid,klev)
-      REAL detr_star1(ngrid,klev)
-      REAL detr_star2(ngrid,klev)
-
       REAL entr_star(ngrid,klev)
-      REAL entr_star1(ngrid,klev)
-      REAL entr_star2(ngrid,klev)
       REAL detr(ngrid,klev)
       REAL entr(ngrid,klev)
+
+      REAL csc(ngrid,klev)
 
       REAL zw2(ngrid,klev+1)
@@ -72,4 +64,7 @@
       REAL zqsatth(ngrid,klev)
       REAL zta_est(ngrid,klev)
+      REAL zdw2
+      REAL zw2modif
+      REAL zeps
 
       REAL linter(ngrid)
@@ -80,4 +75,5 @@
       INTEGER ig,l,k
 
+      real zdz,zfact,zbuoy,zalpha
       real zcor,zdelta,zcvm5,qlbef
       real Tbef,qsatbef
@@ -86,7 +82,9 @@
       PARAMETER (DDT0=.01)
       logical Zsat
-      REAL fact_gamma,fact_epsilon
+      LOGICAL active(ngrid),activetmp(ngrid)
+      REAL fact_gamma,fact_epsilon,fact_gamma2
       REAL c2(ngrid,klev)
 
+      REAL zw2fact,expa
       Zsat=.false.
 ! Initialisation
@@ -97,124 +95,130 @@
          fact_epsilon=1.
       else if (iflag_thermals_ed==1)  then
+! Valeurs au moment de la premiere soumission des papiers
          fact_gamma=1.
-         fact_epsilon=1.
+         fact_epsilon=0.002
+         fact_gamma2=0.6
+! Suggestions des Fleurs, Septembre 2009
+         fact_epsilon=0.015
+!test cr
+!         fact_epsilon=0.002
+         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
 
-      do l=1,klev
+      write(lunout,*)'THERM 31H '
+
+! Initialisations des variables reeles
+if (1==0) then
+      ztva(:,:)=ztv(:,:)
+      ztva_est(:,:)=ztva(:,:)
+      ztla(:,:)=zthl(:,:)
+      zqta(:,:)=po(:,:)
+      zha(:,:) = ztva(:,:)
+else
+      ztva(:,:)=0.
+      ztva_est(:,:)=0.
+      ztla(:,:)=0.
+      zqta(:,:)=0.
+      zha(:,:) =0.
+endif
+
+      zqla_est(:,:)=0.
+      zqsatth(:,:)=0.
+      zqla(:,:)=0.
+      detr_star(:,:)=0.
+      entr_star(:,:)=0.
+      alim_star(:,:)=0.
+      alim_star_tot(:)=0.
+      csc(:,:)=0.
+      detr(:,:)=0.
+      entr(:,:)=0.
+      zw2(:,:)=0.
+      w_est(:,:)=0.
+      f_star(:,:)=0.
+      wa_moy(:,:)=0.
+      linter(:)=1.
+      linter(:)=1.
+
+! Initialisation des variables entieres
+      lmix(:)=1
+      lmix_bis(:)=2
+      wmaxa(:)=0.
+      lalim(:)=1
+
+!-------------------------------------------------------------------------
+! On ne considere comme actif que les colonnes dont les deux premieres
+! couches sont instables.
+!-------------------------------------------------------------------------
+      active(:)=ztv(:,1)>ztv(:,2)
+
+!-------------------------------------------------------------------------
+! Definition de l'alimentation a l'origine dans thermcell_init
+!-------------------------------------------------------------------------
+      do l=1,klev-1
          do ig=1,ngrid
-            zqla_est(ig,l)=0.
-            ztva_est(ig,l)=ztva(ig,l)
-            zqsatth(ig,l)=0.
+            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
+               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
+     &                       *sqrt(zlev(ig,l+1)) 
+               lalim(:)=l+1
+               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
+            endif
          enddo
       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
-         do ig=1,ngrid
-            ztva(ig,k)=ztv(ig,k)
-            ztla(ig,k)=zthl(ig,k)
-            zqla(ig,k)=0.
-            zqta(ig,k)=po(ig,k)
-!
-            ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
-            ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
-            zha(ig,k) = ztva(ig,k)
-!
-         enddo
-      enddo 
-      do k=1,klev
-        do ig=1,ngrid
-           detr_star(ig,k)=0.
-           entr_star(ig,k)=0.
-
-           detr_stara(ig,k)=0.
-           detr_starb(ig,k)=0.
-           detr_starc(ig,k)=0.
-           detr_star0(ig,k)=0.
-           zqla0(ig,k)=0.
-           detr_star1(ig,k)=0.
-           detr_star2(ig,k)=0.
-           entr_star1(ig,k)=0.
-           entr_star2(ig,k)=0.
-
-           detr(ig,k)=0.
-           entr(ig,k)=0.
-        enddo
-      enddo
-      if (prt_level.ge.1) print*,'7 OK convect8'
-      do k=1,klev+1
-         do ig=1,ngrid
-            zw2(ig,k)=0.
-            w_est(ig,k)=0.
-            f_star(ig,k)=0.
-            wa_moy(ig,k)=0.
+      do l=1,klev
+         do ig=1,ngrid 
+            if (alim_star_tot(ig) > 1.e-10 ) then
+               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
+            endif
          enddo
       enddo
-
-      if (prt_level.ge.1) print*,'8 OK convect8'
-      do ig=1,ngrid
-         linter(ig)=1.
-         lmix(ig)=1
-         lmix_bis(ig)=2
-         wmaxa(ig)=0.
-      enddo
-
-!-----------------------------------------------------------------------------------
-!boucle de calcul de la vitesse verticale dans le thermique
-!-----------------------------------------------------------------------------------
-      do l=1,klev-1
-         do ig=1,ngrid
-
-
-
-! Calcul dans la premiere couche active du thermique (ce qu'on teste
-! en disant que la couche est instable et que w2 en bas de la couche
-! est nulle.
-
-            if (ztv(ig,l).gt.ztv(ig,l+1)  &
-     &         .and.alim_star(ig,l).gt.1.e-10  &
-     &         .and.zw2(ig,l).lt.1e-10) then
-
-
+      alim_star_tot(:)=1.
+
+
+!------------------------------------------------------------------------------
+! Calcul dans la premiere couche
+! 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
+! dans une couche l>1
+!------------------------------------------------------------------------------
+do ig=1,ngrid
 ! Le panache va prendre au debut les caracteristiques de l'air contenu
 ! dans cette couche.
-               ztla(ig,l)=zthl(ig,l) 
-               zqta(ig,l)=po(ig,l)
-               zqla(ig,l)=zl(ig,l)
-               f_star(ig,l+1)=alim_star(ig,l)
-
-               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
-     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
-     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
-               w_est(ig,l+1)=zw2(ig,l+1)
+    if (active(ig)) then
+    ztla(ig,1)=zthl(ig,1) 
+    zqta(ig,1)=po(ig,1)
+    zqla(ig,1)=zl(ig,1)
+!cr: attention, prise en compte de f*(1)=1
+    f_star(ig,2)=alim_star(ig,1)
+    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
+&                     *(zlev(ig,2)-zlev(ig,1))  &
+&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
+    w_est(ig,2)=zw2(ig,2)
+    endif
+enddo
 !
 
-
-            else if ((zw2(ig,l).ge.1e-10).and.  &
-     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
-!estimation du detrainement a partir de la geometrie du pas precedent
-!tests sur la definition du detr
-!calcul de detr_star et entr_star
-
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH le test miraculeux de Catherine ? Le bout du tunel ?
-!               w_est(ig,3)=zw2(ig,2)*  &
-!    &                   ((f_star(ig,2))**2)  &
-!    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
-!    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
-!    &                   *(zlev(ig,3)-zlev(ig,2))
-!     if (l.gt.2) then
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!==============================================================================
+!boucle de calcul de la vitesse verticale dans le thermique
+!==============================================================================
+do l=2,klev-1
+!==============================================================================
+
+
+! On decide si le thermique est encore actif ou non
+! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
+    do ig=1,ngrid
+       active(ig)=active(ig) &
+&                 .and. zw2(ig,l)>1.e-10 &
+&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
+    enddo
 
 
@@ -222,62 +226,29 @@
 ! Premier calcul de la vitesse verticale a partir de la temperature
 ! potentielle virtuelle
-
-! FH CESTQUOI CA ????
-#define int1d2
-!#undef int1d2
-#ifdef int1d2
-      if (l.ge.2) then
-#else
-      if (l.gt.2) then
-#endif
-
-      if (1.eq.1) then
-          w_est(ig,3)=zw2(ig,2)* &
-     &      ((f_star(ig,2))**2) &
-     &      /(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
+!     if (1.eq.1) then
+!         w_est(ig,3)=zw2(ig,2)* &
+!    &      ((f_star(ig,2))**2) &
+!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
+!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
+!    &      *(zlev(ig,3)-zlev(ig,2))
+!      endif
 
 
 !---------------------------------------------------------------------------
-!calcul de l entrainement et du detrainement lateral
+! calcul des proprietes thermodynamiques et de la vitesse de la couche l
+! sans tenir compte du detrainement et de l'entrainement dans cette
+! couche
+! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
+! avant) a l'alimentation pour avoir un calcul plus propre
 !---------------------------------------------------------------------------
-!
-!test:estimation de ztva_new_est sans entrainement
-
-               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
-               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
-               qsatbef=MIN(0.5,qsatbef)
-               zcor=1./(1.-retv*qsatbef)
-               qsatbef=qsatbef*zcor
-               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
-               if (Zsat) then
-               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
-               DT = 0.5*RLvCp*qlbef
-               do while (abs(DT).gt.DDT0)
-                 Tbef=Tbef+DT
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
-                 qsatbef=MIN(0.5,qsatbef)
-                 zcor=1./(1.-retv*qsatbef)
-                 qsatbef=qsatbef*zcor
-                 qlbef=zqta(ig,l-1)-qsatbef
-
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
-                 zcor=1./(1.-retv*qsatbef)
-                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
-                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
-                 denom=1.+RLvCp*dqsat_dT
-                 DT=num/denom
-               enddo
-                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef) 
-               endif
+
+     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
         ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
+        zta_est(ig,l)=ztva_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))
@@ -287,412 +258,117 @@
      &                   /(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
                 w_est(ig,l+1)=zw2(ig,l)
              endif
-!
-!calcul du detrainement
-!=======================
-
-!CR:on vire les modifs
-         if (iflag_thermals_ed==0) then
-
-! Modifications du calcul du detrainement.
-! Dans la version de la these de Catherine, on passe brusquement
-! de la version seche a la version nuageuse pour le detrainement
-! ce qui peut occasioner des oscillations.
-! dans la nouvelle version, on commence par calculer un detrainement sec.
-! Puis un autre en cas de nuages.
-! Puis on combine les deux lineairement en fonction de la quantite d'eau.
-
-#define int1d3
-!#undef int1d3
-#define RIO_TH
-#ifdef RIO_TH
-!1. Cas non nuageux
-! 1.1 on est sous le zmax_sec et w croit
-          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
-     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
-#ifdef int1d3
-     &       (zqla_est(ig,l).lt.1.e-10)) then 
-#else
-     &       (zqla(ig,l-1).lt.1.e-10)) then 
-#endif
-             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
-     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
-     &       /(r_aspect*zmax_sec(ig)))
-             detr_stara(ig,l)=detr_star(ig,l)
-
-       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
-
-! 1.2 on est sous le zmax_sec et w decroit
-          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
-#ifdef int1d3
-     &            (zqla_est(ig,l).lt.1.e-10)) then
-#else
-     &            (zqla(ig,l-1).lt.1.e-10)) then
-#endif
-             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
-     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
-     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
-     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
-     &       *((zmax_sec(ig)-zlev(ig,l))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
-             detr_starb(ig,l)=detr_star(ig,l)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
-
+       endif
+    enddo
+
+!-------------------------------------------------
+!calcul des taux d'entrainement et de detrainement
+!-------------------------------------------------
+
+     do ig=1,ngrid
+        if (active(ig)) then
+          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  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
+!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
+!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
+!     &     +0.0006 )
+!           else
+!           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
 
-! 1.3 dans les autres cas
-             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
-     &                      *(zlev(ig,l+1)-zlev(ig,l))
-             detr_starc(ig,l)=detr_star(ig,l)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
-             
+! 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
-
-#else
-
-! 1.1 on est sous le zmax_sec et w croit
-          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
-     &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
-             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
-     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
-     &       /(r_aspect*zmax_sec(ig)))
-
-       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
-
-! 1.2 on est sous le zmax_sec et w decroit
-          else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
-             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
-     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
-     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
-     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
-     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
-     &       *((zmax_sec(ig)-zlev(ig,l))/  &
-     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
-       if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
-
-          else
-             detr_star=0.
-          endif
-
-! 1.3 dans les autres cas
-          detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
-     &                      *(zlev(ig,l+1)-zlev(ig,l))
-
-          coefc=min(zqla(ig,l-1)/1.e-3,1.)
-          if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
-          coefc=1.
-! il semble qu'il soit important de baser le calcul sur
-! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
-          detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
-
-#endif
-
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
-!IM 730508 beg
-!        if(itap.GE.7200) THEN
-!         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
-!        endif
-!IM 730508 end
-         
-         zqla0(ig,l)=zqla_est(ig,l)
-         detr_star0(ig,l)=detr_star(ig,l)
-!IM 060508 beg
-!         if(detr_star(ig,l).GT.1.) THEN
-!          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
-!   &      ,detr_starc(ig,l),coefc
-!         endif
-!IM 060508 end
-!IM 160508 beg
-!IM 160508       IF (f0(ig).NE.0.) THEN
-           detr_star(ig,l)=detr_star(ig,l)/f0(ig)
-!IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
-!IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
-!IM 160508       ELSE
-!IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
-!IM 160508       ENDIF
-!IM 160508 end
-!IM 060508 beg
-!        if(detr_star(ig,l).GT.1.) THEN
-!         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
-!   &     float(1)/f0(ig)
-!        endif
-!IM 060508 end
-        if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
-!
-!calcul de entr_star
-
-! #undef test2
-! #ifdef test2
-! La version test2 destabilise beaucoup le modele.
-! Il semble donc que ca aide d'avoir un entrainement important sous
-! le nuage.
-!         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
-!          entr_star(ig,l)=0.4*detr_star(ig,l)
-!         else
-!          entr_star(ig,l)=0.
-!         endif
-! #else
-!
-! Deplacement du calcul de entr_star pour eviter d'avoir aussi
-! entr_star > fstar.
-! Redeplacer suite a la transformation du cas detr>f
-! FH
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
-#define int1d
-!FH 070508 #define int1d4
-!#undef int1d4
-! L'option int1d4 correspond au choix dans le cas ou le detrainement
-! devient trop grand.
-
-#ifdef int1d
-
-#ifdef int1d4
-#else
-       detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
-!FH 070508 plus
-       detr_star(ig,l)=min(detr_star(ig,l),1.)
-#endif
-
-       entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
-#ifdef int1d4
-! Si le detrainement excede le flux en bas + l'entrainement, le thermique
-! doit disparaitre.
-       if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
-          detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
-          f_star(ig,l+1)=0.
-          linter(ig)=l+1
-          zw2(ig,l+1)=-1.e-10
-       endif
-#endif
-
-
-#else
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
-        if(l.gt.lalim(ig)) then
-         entr_star(ig,l)=0.4*detr_star(ig,l)
-        else
-
-! FH :
-! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
-! en haut de la couche d'alimentation.
-! A remettre en questoin a la premiere occasion mais ca peut aider a 
-! ecrire un code robuste.
-! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
-! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
-! d'alimentation, ce qui n'est pas forcement heureux.
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
-#undef pre_int1c
-#ifdef pre_int1c
-         entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
-         detr_star(ig,l)=entr_star(ig,l)
-#else
-         entr_star(ig,l)=0.
-#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))
+
+! En dessous de lalim, on prend le max de alim_star et entr_star pour
+! alim_star et 0 sinon
+        if (l.lt.lalim(ig)) then
+          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
+          entr_star(ig,l)=0.
         endif
 
-#endif
-
-        if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
-        entr_star1(ig,l)=entr_star(ig,l)
-        detr_star1(ig,l)=detr_star(ig,l)
-!
-
-#ifdef int1d
-#else
-        if (detr_star(ig,l).gt.f_star(ig,l)) then
-
-!  Ce test est là entre autres parce qu'on passe par des valeurs
-!  delirantes de detr_star.
-!  ca vaut sans doute le coup de verifier pourquoi.
-
-           detr_star(ig,l)=f_star(ig,l)
-#ifdef pre_int1c
-           if (l.gt.lalim(ig)+1) then
-               entr_star(ig,l)=0.
-               alim_star(ig,l)=0.
-! FH ajout pour forcer a stoper le thermique juste sous le sommet
-! de la couche (voir calcul de finter)
-               zw2(ig,l+1)=-1.e-10
-               linter(ig)=l+1
-            else
-               entr_star(ig,l)=0.4*detr_star(ig,l)
-            endif
-#else
-           entr_star(ig,l)=0.4*detr_star(ig,l)
-#endif
-        endif
-#endif
-
-      else !l > 2
-         detr_star(ig,l)=0.
-         entr_star(ig,l)=0.
-      endif
-
-        entr_star2(ig,l)=entr_star(ig,l)
-        detr_star2(ig,l)=detr_star(ig,l)
-        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)
-            if(prt_level.GE.10) 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
-! dans la couche d'alimentation
-!pas d entrainement dans la couche alim
-!      if ((l.le.lalim(ig))) then
-!           entr_star(ig,l)=0.
-!      endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-!prise en compte du detrainement et de l entrainement dans le calcul du flux
-
+! Calcul du flux montant normalise
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
      &              -detr_star(ig,l)
 
-!test sur le signe de f_star
-        if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
-       if (f_star(ig,l+1).gt.1.e-10) then 
+      endif
+   enddo
+
 !----------------------------------------------------------------------------
 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
 !---------------------------------------------------------------------------
-!
-       Zsat=.false.
-       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
+   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
+   do ig=1,ngrid
+       if (activetmp(ig)) then 
+           Zsat=.false.
+           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
      &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
      &            /(f_star(ig,l+1)+detr_star(ig,l))
-!
-       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
+           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
      &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
      &            /(f_star(ig,l+1)+detr_star(ig,l))
-!  
-               Tbef=ztla(ig,l)*zpspsk(ig,l)
-               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
-               qsatbef=MIN(0.5,qsatbef)
-               zcor=1./(1.-retv*qsatbef)
-               qsatbef=qsatbef*zcor
-               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
-               if (Zsat) then
-               qlbef=max(0.,zqta(ig,l)-qsatbef)
-               DT = 0.5*RLvCp*qlbef
-               do while (abs(DT).gt.DDT0)
-                 Tbef=Tbef+DT
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
-                 qsatbef=MIN(0.5,qsatbef)
-                 zcor=1./(1.-retv*qsatbef)
-                 qsatbef=qsatbef*zcor
-                 qlbef=zqta(ig,l)-qsatbef
-
-                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
-                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
-                 zcor=1./(1.-retv*qsatbef)
-                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
-                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
-                 denom=1.+RLvCp*dqsat_dT
-                 DT=num/denom
-              enddo
-                 zqla(ig,l) = max(0.,qlbef) 
-              endif
-!    
+
+        endif
+    enddo
+
+   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
+
+
+   do ig=1,ngrid
+      if (activetmp(ig)) then
         if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
 ! on ecrit de maniere conservative (sat ou non)
@@ -707,43 +383,55 @@
 !on ecrit zqsat 
            zqsatth(ig,l)=qsatbef  
-!calcul de vitesse
-           zw2(ig,l+1)=zw2(ig,l)*  &
-     &                 ((f_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. &
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!          zw2(ig,l+1)=&
+!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
+!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
+!     &                 *1./(1.+fact_gamma) &
 !     &                 *(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
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! La meme en plus modulaire :
+           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
+           zdz=zlev(ig,l+1)-zlev(ig,l)
+
+
+           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
+
+if (1==0) then
+           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
+           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
+           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) &
+      &    +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
+
         if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
 !
+!---------------------------------------------------------------------------
 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
-
+!---------------------------------------------------------------------------
+
+   do ig=1,ngrid
             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'
@@ -753,5 +441,4 @@
             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))  &
@@ -771,28 +458,22 @@
             wmaxa(ig)=wa_moy(ig,l+1)
         endif
-        enddo
+   enddo
+
+!=========================================================================
+! FIN DE LA BOUCLE VERTICALE
       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
+!=========================================================================
+
+!on recalcule alim_star_tot
+       do ig=1,ngrid
+          alim_star_tot(ig)=0.
+       enddo
+       do ig=1,ngrid
+          do l=1,lalim(ig)-1
+          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
+          enddo
+       enddo
+       
+
         if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
 
