Index: /LMDZ4/trunk/libf/phylmd/add_phys_tend.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/add_phys_tend.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/add_phys_tend.F90	(revision 972)
@@ -18,4 +18,5 @@
 use phys_state_var_mod
 IMPLICIT none
+#include "iniprint.h"
 
 ! Arguments :
@@ -32,12 +33,20 @@
 INTEGER jadrs(klon*klev), jbad
 INTEGER jqadrs(klon*klev), jqbad
+INTEGER kadrs(klon*klev)
+INTEGER kqadrs(klon*klev)
 
 integer debug_level
-
+logical, save :: first=.true.
+INTEGER, SAVE :: itap
 !======================================================================
 ! Initialisations
 
 debug_level=10
-
+     if (first) then
+        itap=0
+        first=.false.
+     endif
+! Incrementer le compteur de la physique
+     itap   = itap + 1
 !======================================================================
 ! Ajout des tendances sur le vent et l'eau liquide
@@ -62,8 +71,10 @@
             jbad = jbad + 1
             jadrs(jbad) = i
+            kadrs(jbad) = k
             ENDIF
             IF ( zq<0. .or. zq>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
             jqbad = jqbad + 1
             jqadrs(jqbad) = i
+            kqadrs(jqbad) = k
             ENDIF
             t_seri(i,k)=zt
@@ -97,5 +108,12 @@
          print*,'l    T     dT       Q     dQ    '
          DO k = 1, klev
-            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
+           zq=q_seri(i,k)+zdq(i,k)
+           if (zq.lt.1.e-15) then
+              if (q_seri(i,k).lt.1.e-15) then
+!              print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
+               q_seri(i,k)=1.e-15
+               zdq(i,k)=(1.e-15-q_seri(i,k))
+              endif
+           endif
 !           zq=q_seri(i,k)+zdq(i,k)
 !           if (zq.lt.1.e-15) then
@@ -103,9 +121,65 @@
 !           endif
          ENDDO
-         call print_debug_phys(i,debug_level,text)
       ENDDO
 ENDIF
 !
 
+!IM ajout memes tests pour reverifier les jbad, jqbad beg
+      jbad=0
+      jqbad=0
+      DO k = 1, klev
+         DO i = 1, klon
+            IF ( t_seri(i,k)>370. .or. t_seri(i,k)<130. .or. abs(zdt(i,k))>50. ) then
+            jbad = jbad + 1
+            jadrs(jbad) = i
+!           if(prt_level.ge.10) THEN 
+!             print*,'cas2 i k t_seri zdt',i,k,t_seri(i,k),zdt(i,k)
+!           endif
+            ENDIF
+            IF ( q_seri(i,k)<0. .or. q_seri(i,k)>0.1 .or. abs(zdq(i,k))>1.e-2 ) then
+            jqbad = jqbad + 1
+            jqadrs(jqbad) = i
+            kqadrs(jqbad) = k
+!           if(prt_level.ge.10) THEN 
+!             print*,'cas2 i k q_seri zdq',i,k,q_seri(i,k),zdq(i,k)
+!           endif
+            ENDIF
+         ENDDO
+      ENDDO
+IF (jbad .GT. 0) THEN
+      DO j = 1, jbad
+         i=jadrs(j)
+         k=kadrs(j)
+         print*,'PLANTAGE2 POUR LE POINT i itap rlon rlat txt jbad zdt t',i,itap,rlon(i),rlat(i),text,jbad, &
+       &        zdt(i,k),t_seri(i,k)-zdt(i,k)
+!        if(prt_level.ge.10) THEN 
+         if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN 
+          print*,'l    T     dT       Q     dQ    '
+          DO k = 1, klev
+             write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
+          ENDDO
+          call print_debug_phys(i,debug_level,text)
+         endif
+      ENDDO
+ENDIF
+!
+IF (jqbad .GT. 0) THEN
+      DO j = 1, jqbad
+         i=jqadrs(j)
+         k=kqadrs(j)
+         print*,'WARNING  : EAU2 POUR LE POINT i itap rlon rlat txt jqbad zdq q zdql ql',i,itap,rlon(i),rlat(i),text,jqbad,&
+       &        zdq(i,k), q_seri(i,k)-zdq(i,k), zdql(i,k), ql_seri(i,k)-zdql(i,k)
+!        if(prt_level.ge.10) THEN 
+         if(prt_level.ge.10.and.itap.GE.229.and.i.EQ.3027) THEN 
+          print*,'l    T     dT       Q     dQ    '
+          DO k = 1, klev
+            write(*,'(i3,2f14.4,2e14.2)') k,t_seri(i,k),zdt(i,k),q_seri(i,k),zdq(i,k)
+          ENDDO
+          call print_debug_phys(i,debug_level,text)
+         endif
+      ENDDO
+ENDIF
+
+      CALL hgardfou(t_seri,ftsol,text)
       RETURN
       END
Index: /LMDZ4/trunk/libf/phylmd/cv3a_compress.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/cv3a_compress.F	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/cv3a_compress.F	(revision 972)
@@ -123,5 +123,5 @@
 
       if (nn.ne.ncum) then
-         print*,'strange! nn not equal to ncum: ',nn,ncum
+         print*,'WARNING nn not equal to ncum: ',nn,ncum
          stop
       endif
@@ -150,4 +150,9 @@
  150  continue
 
+      if (nn.ne.ncum) then
+         print*,'WARNING nn not equal to ncum: ',nn,ncum
+         stop
+      endif
+
       RETURN
       END
Index: /LMDZ4/trunk/libf/phylmd/hgardfou.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/hgardfou.F	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/hgardfou.F	(revision 972)
@@ -51,6 +51,6 @@
            ok = .FALSE.
            DO i = 1, jbad
-             PRINT *,'i,k,temperature =',jadrs(i),k,zt(jadrs(i)),
-     $       rlon(jadrs(i)),rlat(jadrs(i))
+            PRINT *,'i,k,temperature rlon rlat=',jadrs(i),k,zt(jadrs(i))
+     $      ,rlon(jadrs(i)),rlat(jadrs(i))
            ENDDO
          ENDIF
@@ -70,6 +70,6 @@
            ok = .FALSE.
            DO i = 1, jbad
-             PRINT *,'i,k,temperature =',jadrs(i),k,zt(jadrs(i)),
-     $       rlon(jadrs(i)),rlat(jadrs(i))
+            PRINT *,'i,k,temperature rlon rlat=',jadrs(i),k,zt(jadrs(i))
+     $      ,rlon(jadrs(i)),rlat(jadrs(i))
            ENDDO
          ENDIF
Index: /LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90	(revision 972)
@@ -880,6 +880,6 @@
 !****************************************************************************************
 ! H and Q
-       print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
-       print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
+!      print *,'pbl_surface: ok_flux_surf=',ok_flux_surf
+!      print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
        if (ok_flux_surf) then
           y_flux_t1(:) =  fsens
Index: /LMDZ4/trunk/libf/phylmd/thermcell.h
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell.h	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell.h	(revision 972)
@@ -1,9 +1,9 @@
       integer iflag_thermals,nsplit_thermals
-      real r_aspect_thermals,l_mix_thermals,tho_thermals
+      real r_aspect_thermals,l_mix_thermals,tau_thermals
       integer w2di_thermals,isplit
       integer iflag_coupl,iflag_clos,iflag_wake
 
       common/ctherm1/iflag_thermals,nsplit_thermals
-      common/ctherm2/r_aspect_thermals,l_mix_thermals,tho_thermals
+      common/ctherm2/r_aspect_thermals,l_mix_thermals,tau_thermals
       common/ctherm3/w2di_thermals
       common/ctherm4/iflag_coupl,iflag_clos,iflag_wake
Index: /LMDZ4/trunk/libf/phylmd/thermcell_closure.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_closure.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_closure.F90	(revision 972)
@@ -1,4 +1,4 @@
       SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
+     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
 
 !-------------------------------------------------------------------------
@@ -24,5 +24,4 @@
 
       REAL f(ngrid)
-      REAL f0(ngrid)
 
       do ig=1,ngrid
@@ -52,13 +51,13 @@
              f(ig)=wmax_sec(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))
+!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
+!    &                     zmax_sec(ig))*wmax_sec(ig))
              else
              f(ig)=wmax(ig)/zdenom
-            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
-     &                     zmax(ig))*wmax(ig))
+!            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
+!     &                     zmax(ig))*wmax(ig))
              endif
          endif
-         f0(ig)=f(ig)
+!         f0(ig)=f(ig)
       enddo
       if (prt_level.ge.1) print*,'apres fermeture'
Index: /LMDZ4/trunk/libf/phylmd/thermcell_dq.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_dq.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_dq.F90	(revision 972)
@@ -23,10 +23,45 @@
       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
 
+      real zzm
+
       integer ig,k
+      real cfl
+
+      real qold(ngrid,nlay)
+      real ztimestep
+      integer niter,iter
+
+
+
+! Calcul du critere CFL pour l'advection dans la subsidence
+      do k=1,nlay
+         do ig=1,ngrid
+            zzm=masse(ig,k)/ptimestep
+            cfl=max(cfl,fm(ig,k)/zzm)
+            if (entr(ig,k).gt.zzm) then
+               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
+               stop
+            endif
+         enddo
+      enddo
+
+!IM 090508     print*,'CFL CFL CFL CFL ',cfl
+
+#undef CFL
+#ifdef CFL
+! On subdivise le calcul en niter pas de temps.
+      niter=int(cfl)+1
+#else
+      niter=1
+#endif
+
+      ztimestep=ptimestep/niter
+      qold=q
+
+
+do iter=1,niter
+      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
 
 !   calcul du detrainement
-
-      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
-
       do k=1,nlay
          do ig=1,ngrid
@@ -56,5 +91,5 @@
       do k=2,nlay
          do ig=1,ngrid
-            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
+            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
      &         1.e-5*masse(ig,k)) then
          qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
@@ -72,8 +107,27 @@
       enddo
 
+! Calcul du flux subsident
+
       do k=2,nlay
          do ig=1,ngrid
-!             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
+#undef centre
+#ifdef centre
+             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
+#else
+
+#define plusqueun
+#ifdef plusqueun
+! Schema avec advection sur plus qu'une maille.
+            zzm=masse(ig,k)/ztimestep
+            if (fm(ig,k)>zzm) then
+               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
+            else
+               wqd(ig,k)=fm(ig,k)*q(ig,k)
+            endif
+#else
             wqd(ig,k)=fm(ig,k)*q(ig,k)
+#endif
+#endif
+
             if (wqd(ig,k).lt.0.) then
 !               print*,'wqd<0!!!'
@@ -86,9 +140,11 @@
       enddo
      
+
+! Calcul des tendances
       do k=1,nlay
          do ig=1,ngrid
-            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
+            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
      &               -wqd(ig,k)+wqd(ig,k+1))  &
-     &               /masse(ig,k)
+     &               *ztimestep/masse(ig,k)
 !            if (dq(ig,k).lt.0.) then
 !               print*,'dq<0!!!'
@@ -97,4 +153,16 @@
       enddo
 
+
+enddo
+
+
+! Calcul des tendances
+      do k=1,nlay
+         do ig=1,ngrid
+            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
+            q(ig,k)=qold(ig,k)
+         enddo
+      enddo
+
       return
       end
Index: /LMDZ4/trunk/libf/phylmd/thermcell_dv2.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_dv2.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_dv2.F90	(revision 972)
@@ -54,4 +54,6 @@
       enddo
 
+      print*,'WARNING on initialise gamma(1:ngrid,1)=0.'
+      gamma(1:ngrid,1)=0.
       do k=2,nlay
          do ig=1,ngrid
@@ -60,8 +62,24 @@
 !   On itère sur la valeur du coeff de freinage.
 !              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
+!IM 060508 beg
+!             if(0.5*(fraca(ig,k+1)+fraca(ig,k)).LT.0.) THEN
+!              print*,'th_dv2 ig k fraca(:,k) fraca(:k+1)', &
+!    &         ig,k,fraca(ig,k),fraca(ig,k+1)
+!             endif
+!             if(larga(ig).EQ.0.) THEN
+!              print*,'th_dv2 ig larga=0.',ig
+!             endif
+              if(larga(ig).GT.0.) THEN
+!IM 060508 end
                gamma0=masse(ig,k)  &
      &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
      &         *0.5/larga(ig)  &
      &         *1.
+!IM 060508 beg
+              else 
+               if(prt_level.GE.10) print*,'WARNING cas ELSE on initialise gamma0=0.'
+               gamma0=0.
+              endif !(larga(ig).GT.0.) THEN
+!IM 060508 end
 !    s         *0.5
 !              gamma0=0.
@@ -117,4 +135,9 @@
       do k=1,nlay
          do ig=1,ngrid
+!IM
+         if(prt_level.GE.10) print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
+     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
+     &   masse(ig,k)
+!
             du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
      &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
Index: /LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90	(revision 972)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_flux2.F90	(revision 972)
@@ -0,0 +1,458 @@
+      SUBROUTINE thermcell_flux2(ngrid,klev,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)
+
+
+!---------------------------------------------------------------------------
+!thermcell_flux: deduction des flux
+!---------------------------------------------------------------------------
+
+      IMPLICIT NONE
+#include "iniprint.h"
+      
+      INTEGER ig,l
+      INTEGER ngrid,klev
+      
+      REAL alim_star(ngrid,klev),entr_star(ngrid,klev)
+      REAL detr_star(ngrid,klev)
+      REAL zw2(ngrid,klev+1)
+      REAL zlev(ngrid,klev+1)
+      REAL masse(ngrid,klev)
+      REAL ptimestep
+      REAL rhobarz(ngrid,klev)
+      REAL f(ngrid)
+      INTEGER lmax(ngrid)
+      INTEGER lalim(ngrid)
+      REAL zqla(ngrid,klev)
+      REAL zmax(ngrid)
+
+      integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
+      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
+      
+
+      REAL entr(ngrid,klev),detr(ngrid,klev)
+      REAL fm(ngrid,klev+1)
+      REAL zfm
+
+      integer igout
+      integer lev_out
+      integer lunout1
+
+      REAL f_old,ddd0,eee0,ddd,eee,zzz
+
+      REAL fomass_max,alphamax
+      save fomass_max,alphamax
+
+      fomass_max=0.5
+      alphamax=0.7
+
+      ncorecfm1=0
+      ncorecfm2=0
+      ncorecfm3=0
+      ncorecfm4=0
+      ncorecfm5=0
+      ncorecfm6=0
+      ncorecfm7=0
+      ncorecfm8=0
+      ncorecalpha=0
+
+!initialisation
+      fm(:,:)=0.
+      
+      if (prt_level.ge.10) then
+         write(lunout1,*) 'Dans thermcell_flux 0'
+         write(lunout1,*) 'flux base ',f(igout)
+         write(lunout1,*) 'lmax ',lmax(igout)
+         write(lunout1,*) 'lalim ',lalim(igout)
+         write(lunout1,*) 'ig= ',igout
+         write(lunout1,*) ' l E*    A*     D*  '
+         write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
+     &    ,l=1,lmax(igout))
+      endif
+
+
+!-------------------------------------------------------------------------
+! Verification de la nullite des entrainement et detrainement au dessus
+! de lmax(ig)
+!-------------------------------------------------------------------------
+
+      do l=1,klev
+         do ig=1,ngrid
+            if (l.le.lmax(ig)) then
+               if (entr_star(ig,l).gt.1.) then
+                    print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
+                    print*,'entr_star(ig,l)',entr_star(ig,l)
+                    print*,'alim_star(ig,l)',alim_star(ig,l)
+                    print*,'detr_star(ig,l)',detr_star(ig,l)
+!                   stop
+               endif
+            else
+               if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
+                    print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
+                    print*,'entr_star(ig,l)',entr_star(ig,l)
+                    print*,'alim_star(ig,l)',alim_star(ig,l)
+                    print*,'detr_star(ig,l)',detr_star(ig,l)
+                    stop
+               endif
+            endif
+         enddo
+      enddo
+
+!-------------------------------------------------------------------------
+! Multiplication par le flux de masse issu de la femreture
+!-------------------------------------------------------------------------
+
+      do l=1,klev
+         entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
+         detr(:,l)=f(:)*detr_star(:,l)
+      enddo
+
+      if (prt_level.ge.10) then
+         write(lunout1,*) 'Dans thermcell_flux 1'
+         write(lunout1,*) 'flux base ',f(igout)
+         write(lunout1,*) 'lmax ',lmax(igout)
+         write(lunout1,*) 'lalim ',lalim(igout)
+         write(lunout1,*) 'ig= ',igout
+         write(lunout1,*) ' l   E    D     W2'
+         write(lunout1,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
+     &    ,zw2(igout,l+1),l=1,lmax(igout))
+      endif
+
+      fm(:,1)=0.
+      do l=1,klev
+         do ig=1,ngrid
+            if (l.lt.lmax(ig)) then
+               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
+            elseif(l.eq.lmax(ig)) then
+               fm(ig,l+1)=0.
+               detr(ig,l)=fm(ig,l)+entr(ig,l)
+            else
+               fm(ig,l+1)=0.
+            endif
+         enddo
+      enddo
+
+
+
+! Test provisoire : pour comprendre pourquoi on corrige plein de fois 
+! le cas fm6, on commence par regarder une premiere fois avant les
+! autres corrections.
+
+      do l=1,klev
+         do ig=1,ngrid
+            if (detr(ig,l).gt.fm(ig,l)) then
+               ncorecfm8=ncorecfm8+1
+!              igout=ig
+            endif
+         enddo
+      enddo
+
+!      if (prt_level.ge.10) &
+!    &    call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
+!    &    ptimestep,masse,entr,detr,fm,'2  ')
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH Version en cours de test;
+! par rapport a thermcell_flux, on fait une grande boucle sur "l"
+! et on modifie le flux avec tous les contrôles appliques d'affilee
+! pour la meme couche
+! Momentanement, on duplique le calcule du flux pour pouvoir comparer
+! les flux avant et apres modif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      do l=1,klev
+
+         do ig=1,ngrid
+            if (l.lt.lmax(ig)) then
+               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
+            elseif(l.eq.lmax(ig)) then
+               fm(ig,l+1)=0.
+               detr(ig,l)=fm(ig,l)+entr(ig,l)
+            else
+               fm(ig,l+1)=0.
+            endif
+         enddo
+
+
+!-------------------------------------------------------------------------
+! Verification de la positivite des flux de masse
+!-------------------------------------------------------------------------
+
+!     do l=1,klev
+         do ig=1,ngrid
+            if (fm(ig,l+1).lt.0.) then
+!              print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
+                ncorecfm1=ncorecfm1+1
+               fm(ig,l+1)=fm(ig,l)
+               detr(ig,l)=entr(ig,l)
+            endif
+         enddo
+!     enddo
+
+      if (prt_level.ge.10) &
+     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
+
+!-------------------------------------------------------------------------
+!Test sur fraca croissant
+!-------------------------------------------------------------------------
+      if (1.eq.1) then 
+!     do l=1,klev
+         do ig=1,ngrid
+          if (l.ge.lalim(ig).and.l.le.lmax(ig) &
+     &    .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
+!  zzz est le flux en l+1 a frac constant
+             zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)  &
+     &                          /(rhobarz(ig,l)*zw2(ig,l))
+             if (fm(ig,l+1).gt.zzz) then
+                detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
+                fm(ig,l+1)=zzz
+                ncorecfm4=ncorecfm4+1
+             endif
+          endif
+        enddo
+!     enddo
+      endif
+
+      if (prt_level.ge.10) &
+     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
+
+
+!-------------------------------------------------------------------------
+!test sur flux de masse croissant
+!-------------------------------------------------------------------------
+      if (1.eq.1) then
+!     do l=1,klev
+         do ig=1,ngrid
+            if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
+              f_old=fm(ig,l+1)
+              fm(ig,l+1)=fm(ig,l)
+              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
+               ncorecfm5=ncorecfm5+1
+            endif
+         enddo
+!     enddo
+      endif
+
+      if (prt_level.ge.10) &
+     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
+
+!fin 1.eq.0
+!-------------------------------------------------------------------------
+!detr ne peut pas etre superieur a fm
+!-------------------------------------------------------------------------
+
+      if(1.eq.1) then
+
+!     do l=1,klev
+         do ig=1,ngrid
+            if (entr(ig,l)<0.) then
+               print*,'N1 ig,l,entr',ig,l,entr(ig,l)
+               stop'entr negatif'
+            endif
+            if (detr(ig,l).gt.fm(ig,l)) then
+               ncorecfm6=ncorecfm6+1
+               detr(ig,l)=fm(ig,l)
+               entr(ig,l)=fm(ig,l+1)
+
+! Dans le cas ou on est au dessus de la couche d'alimentation et que le
+! detrainement est plus fort que le flux de masse, on stope le thermique.
+!test:on commente
+!               if (l.gt.lalim(ig)) then
+!                  lmax(ig)=l
+!                  fm(ig,l+1)=0.
+!                  entr(ig,l)=0.
+!               else
+!                  ncorecfm7=ncorecfm7+1
+!               endif
+            endif
+
+            if(l.gt.lmax(ig)) then
+               detr(ig,l)=0.
+               fm(ig,l+1)=0.
+               entr(ig,l)=0.
+            endif
+
+            if (entr(ig,l).lt.0.) then
+               print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
+               print*,'entr(ig,l)',entr(ig,l)
+               print*,'fm(ig,l)',fm(ig,l)
+               stop'probleme dans thermcell flux'
+            endif
+         enddo
+!     enddo
+      endif
+
+
+      if (prt_level.ge.10) &
+     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
+
+!-------------------------------------------------------------------------
+!fm ne peut pas etre negatif
+!-------------------------------------------------------------------------
+
+!     do l=1,klev
+         do ig=1,ngrid
+            if (fm(ig,l+1).lt.0.) then
+               detr(ig,l)=detr(ig,l)+fm(ig,l+1)
+               fm(ig,l+1)=0.
+!              print*,'fm2<0',l+1,lmax(ig)
+               ncorecfm2=ncorecfm2+1
+            endif
+            if (detr(ig,l).lt.0.) then
+               print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
+               print*,'detr(ig,l)',detr(ig,l)
+               print*,'fm(ig,l)',fm(ig,l)
+               stop'probleme dans thermcell flux'
+            endif
+        enddo
+!    enddo
+
+      if (prt_level.ge.10) &
+     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
+
+!-----------------------------------------------------------------------
+!la fraction couverte ne peut pas etre superieure a 1            
+!-----------------------------------------------------------------------
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH Partie a revisiter.
+! Il semble qu'etaient codees ici deux optiques dans le cas
+! F/ (rho *w) > 1
+! soit limiter la hauteur du thermique en considerant que c'est 
+! la derniere chouche, soit limiter F a rho w.
+! Dans le second cas, il faut en fait limiter a un peu moins
+! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
+! dans thermcell_main et qu'il semble de toutes facons deraisonable
+! d'avoir des fractions de 1..
+! Ci dessous, et dans l'etat actuel, le premier des  deux if est
+! sans doute inutile.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!    do l=1,klev
+        do ig=1,ngrid
+           if (zw2(ig,l+1).gt.1.e-10) then
+           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
+           if ( fm(ig,l+1) .gt. zfm) then
+              f_old=fm(ig,l+1)
+              fm(ig,l+1)=zfm
+!             zw2(ig,l+1)=0.
+!             zqla(ig,l+1)=0.
+              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
+!             lmax(ig)=l+1
+!             zmax(ig)=zlev(ig,lmax(ig))
+!             print*,'alpha>1',l+1,lmax(ig)
+              ncorecalpha=ncorecalpha+1
+           endif
+           endif
+        enddo
+!    enddo
+!
+
+
+      if (prt_level.ge.10) &
+     &   write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l),fm(igout,l+1)
+
+! Fin de la grande boucle sur les niveaux verticaux
+      enddo
+
+!      if (prt_level.ge.10) &
+!    &    call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
+!    &    ptimestep,masse,entr,detr,fm,'8  ')
+
+
+!-----------------------------------------------------------------------
+! On fait en sorte que la quantite totale d'air entraine dans le 
+! panache ne soit pas trop grande comparee a la masse de la maille
+!-----------------------------------------------------------------------
+
+      if (1.eq.1) then
+      do l=1,klev-1
+         do ig=1,ngrid
+            eee0=entr(ig,l)
+            ddd0=detr(ig,l)
+            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
+            ddd=detr(ig,l)-eee
+            if (eee.gt.0.) then
+                ncorecfm3=ncorecfm3+1
+                entr(ig,l)=entr(ig,l)-eee
+                if ( ddd.gt.0.) then
+!   l'entrainement est trop fort mais l'exces peut etre compense par une
+!   diminution du detrainement)
+                   detr(ig,l)=ddd
+                else
+!   l'entrainement est trop fort mais l'exces doit etre compense en partie
+!   par un entrainement plus fort dans la couche superieure
+                   if(l.eq.lmax(ig)) then
+                      detr(ig,l)=fm(ig,l)+entr(ig,l)
+                   else
+                      if(l.ge.lmax(ig).and.0.eq.1) then
+                         print*,'ig,l',ig,l
+                         print*,'eee0',eee0
+                         print*,'ddd0',ddd0
+                         print*,'eee',eee
+                         print*,'ddd',ddd
+                         print*,'entr',entr(ig,l)
+                         print*,'detr',detr(ig,l)
+                         print*,'masse',masse(ig,l)
+                         print*,'fomass_max',fomass_max
+                         print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
+                         print*,'ptimestep',ptimestep
+                         print*,'lmax(ig)',lmax(ig)
+                         print*,'fm(ig,l+1)',fm(ig,l+1)
+                         print*,'fm(ig,l)',fm(ig,l)
+                         stop'probleme dans thermcell_flux'
+                      endif
+                      entr(ig,l+1)=entr(ig,l+1)-ddd
+                      detr(ig,l)=0.
+                      fm(ig,l+1)=fm(ig,l)+entr(ig,l)
+                      detr(ig,l)=0.
+                   endif
+                endif
+            endif
+         enddo
+      enddo
+      endif
+!                  
+!              ddd=detr(ig)-entre
+!on s assure que tout s annule bien en zmax
+      do ig=1,ngrid
+         fm(ig,lmax(ig)+1)=0.
+         entr(ig,lmax(ig))=0.
+         detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
+      enddo
+
+!-----------------------------------------------------------------------
+! Impression du nombre de bidouilles qui ont ete necessaires
+!-----------------------------------------------------------------------
+
+!IM 090508 beg
+!     if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
+!
+!         print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
+!   &     ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
+!   &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
+!   &     ncorecfm6,'x fm6', &
+!   &     ncorecfm7,'x fm7', &
+!   &     ncorecfm8,'x fm8', &
+!   &     ncorecalpha,'x alpha'
+!     endif
+!IM 090508 end
+
+!      if (prt_level.ge.10) &
+!    &    call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
+!    &    ptimestep,masse,entr,detr,fm,'fin')
+
+
+      return
+      end
Index: /LMDZ4/trunk/libf/phylmd/thermcell_init.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_init.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_init.F90	(revision 972)
@@ -1,3 +1,3 @@
-      SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlev,  &
+      SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
      &                  lalim,lmin,alim_star,alim_star_tot,lev_out)
 
@@ -12,5 +12,6 @@
       INTEGER ngrid,nlay
       REAL ztv(ngrid,nlay)
-      REAL zlev(ngrid,nlay)
+      REAL zlay(ngrid,nlay)
+      REAL zlev(ngrid,nlay+1)
 !arguments de sortie
       INTEGER lalim(ngrid)
@@ -20,4 +21,5 @@
       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      
@@ -58,5 +60,24 @@
       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 
@@ -69,4 +90,15 @@
          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
Index: /LMDZ4/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_main.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_main.F90	(revision 972)
@@ -2,15 +2,15 @@
 ! $Header$
 !
-      SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
+      SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep  &
      &                  ,pplay,pplev,pphi,debut  &
      &                  ,pu,pv,pt,po  &
      &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
-     &                  ,fm0,entr0,zqla,lmax  &
+     &                  ,fm0,entr0,detr0,zqla,lmax  &
      &                  ,ratqscth,ratqsdiff,zqsatth  &
-     &                  ,r_aspect,l_mix,w2di,tho &
+     &                  ,r_aspect,l_mix,tau_thermals &
      &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
      &                  ,zmax0, f0)
 
-      use dimphy
+      USE dimphy
       IMPLICIT NONE
 
@@ -50,5 +50,9 @@
 !   ----------
 
-      INTEGER ngrid,nlay,w2di,tho
+!IM 140508
+      INTEGER itap
+
+      INTEGER ngrid,nlay,w2di
+      real tau_thermals
       real ptimestep,l_mix,r_aspect
       REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
@@ -61,4 +65,8 @@
 !   local:
 !   ------
+
+      integer icount
+      data icount/0/
+      save icount
 
       integer,save :: igout=1
@@ -79,7 +87,9 @@
 !FH/IM     save zmax0
 
+      real lambda
+
       real zlev(klon,klev+1),zlay(klon,klev)
       real deltaz(klon,klev)
-      REAL zh(klon,klev),zdhadj(klon,klev)
+      REAL zh(klon,klev)
       real zthl(klon,klev),zdthladj(klon,klev)
       REAL ztv(klon,klev)
@@ -97,5 +107,6 @@
       real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
       real q2(klon,klev)
-!      common/comtherm/thetath2,wth2
+! FH probleme de dimensionnement avec l'allocation dynamique
+!     common/comtherm/thetath2,wth2
     
       real ratqscth(klon,klev)
@@ -109,12 +120,11 @@
 
       logical sorties
-      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
+      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),detr(klon,klev)
-      real detr0(klon,klev)
-      real fm(klon,klev+1),entr(klon,klev)
+      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)
@@ -164,4 +174,22 @@
       seuil=0.25
 
+      if (debut)  then
+         fm0=0.
+         entr0=0.
+         detr0=0.
+      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'
 
@@ -176,13 +204,22 @@
 !Initialisation
 !
-      do ig=1,klon      
+!    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.
+!     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
+!     enddo 
+!    ENDIF !(1.eq.0) THEN
+     print*,'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
 
 !-----------------------------------------------------------------------
@@ -238,4 +275,8 @@
          rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
       enddo
+
+!IM
+      print*,'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
+      rhobarz(:,1)=rho(:,1)
 
       do l=2,nlay
@@ -309,5 +350,5 @@
 !
       entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
-      CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
+      CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
      &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
 
@@ -318,9 +359,9 @@
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
       if (prt_level.ge.10) then
-         write(lunout,*) 'Dans thermcell_main 1'
-         write(lunout,*) 'lmin ',lmin(igout)
-         write(lunout,*) 'lalim ',lalim(igout)
-         write(lunout,*) ' ig l alim_star thetav'
-         write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
+         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
@@ -346,9 +387,9 @@
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
       if (prt_level.ge.10) then
-         write(lunout,*) 'Dans thermcell_main 1b'
-         write(lunout,*) 'lmin ',lmin(igout)
-         write(lunout,*) 'lalim ',lalim(igout)
-         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
-         write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
+         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
@@ -360,8 +401,12 @@
 !--------------------------------------------------------------------------------
 !
-      CALL thermcell_plume(ngrid,nlay,ztv,zthl,po,zl,rhobarz,  &
+      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,  &
      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
      &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
+      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  ')
@@ -369,9 +414,9 @@
       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
       if (prt_level.ge.10) then
-         write(lunout,*) 'Dans thermcell_main 2'
-         write(lunout,*) 'lmin ',lmin(igout)
-         write(lunout,*) 'lalim ',lalim(igout)
-         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
-         write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
+         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
@@ -397,7 +442,19 @@
 
       CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
-     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
+     &   zlev,lalim,alim_star,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
 
 !-------------------------------------------------------------------------------
@@ -405,8 +462,9 @@
 !-------------------------------------------------------------------------------
 
-      CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
+      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
      &       lalim,lmax,alim_star,  &
      &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
-     &       detr,zqla,zmax,lev_out,lunout,igout)
+     &       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'
@@ -414,11 +472,15 @@
       call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
 
-!c------------------------------------------------------------------
-!   calcul du transport vertical
-!------------------------------------------------------------------
-
-      if (w2di.eq.1) then
-         fm0=fm0+ptimestep*(fm-fm0)/float(tho)
-         entr0=entr0+ptimestep*(entr-entr0)/float(tho)
+!------------------------------------------------------------------
+!   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
@@ -426,4 +488,8 @@
          detr0=detr
       endif
+
+!c------------------------------------------------------------------
+!   calcul du transport vertical
+!------------------------------------------------------------------
 
       call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
@@ -453,5 +519,8 @@
 !------------------------------------------------------------------
 
+!IM 090508  
       if (1.eq.1) then
+!IM 070508 vers. _dq       
+!     if (1.eq.0) then
 
 
@@ -461,5 +530,6 @@
          call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
      &    ,fraca,zmax  &
-     &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
+     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
+!IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
       else
 
@@ -478,5 +548,5 @@
       enddo
 
-!     print*,'14 OK convect8'
+      if (prt_level.ge.1) print*,'14 OK convect8'
 !------------------------------------------------------------------
 !   Calculs de diagnostiques pour les sorties
@@ -485,5 +555,5 @@
       
       if (sorties) then
-!      print*,'14a OK convect8'
+      if (prt_level.ge.1) print*,'14a OK convect8'
 ! calcul du niveau de condensation
 ! initialisation
@@ -505,5 +575,5 @@
          enddo
       enddo
-!     print*,'14b OK convect8'
+      if (prt_level.ge.1) print*,'14b OK convect8'
       do k=nlay,1,-1
          do ig=1,ngrid
@@ -514,5 +584,5 @@
          enddo
       enddo
-!     print*,'14c OK convect8'
+      if (prt_level.ge.1) print*,'14c OK convect8'
 !calcul des moments
 !initialisation
@@ -526,14 +596,24 @@
          enddo
       enddo      
-!     print*,'14d OK convect8'
+      if (prt_level.ge.1) print*,'14d OK convect8'
+      print*,'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
-            wth2(ig,l)=zf2*(zw2(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
@@ -548,10 +628,12 @@
       n_int(ig)=0
       enddo
-      do ig=1,ngrid
-!      do l=nivcon(ig),lmax(ig)
-      do l=1,lmax(ig)
-      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
+!
+      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
@@ -638,5 +720,5 @@
 
 !calcul du ratqscdiff
-!     print*,'14e OK convect8'
+      if (prt_level.ge.1) print*,'14e OK convect8'
       var=0.
       vardiff=0.
@@ -647,5 +729,5 @@
          enddo
       enddo
-!     print*,'14f OK convect8'
+      if (prt_level.ge.1) print*,'14f OK convect8'
       do ig=1,ngrid
           do l=1,lalim(ig)
@@ -658,5 +740,5 @@
           enddo
       enddo
-!     print*,'14g OK convect8'
+      if (prt_level.ge.1) print*,'14g OK convect8'
       do l=1,nlay
          do ig=1,ngrid
@@ -679,5 +761,6 @@
 
 
-! #define troisD
+#define troisD
+#undef troisD
       if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
 #ifdef troisD
@@ -689,4 +772,5 @@
       if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
 
+!     if(icount.eq.501) stop'au pas 301 dans thermcell_main'
       return
       end
@@ -717,5 +801,6 @@
 !  test sur la hauteur des thermiques ...
          do i=1,klon
-            if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
+!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'
@@ -724,5 +809,5 @@
                enddo
 !              stop
-            endif
+           endif
          enddo
 
Index: /LMDZ4/trunk/libf/phylmd/thermcell_plume.F90
===================================================================
--- /LMDZ4/trunk/libf/phylmd/thermcell_plume.F90	(revision 971)
+++ /LMDZ4/trunk/libf/phylmd/thermcell_plume.F90	(revision 972)
@@ -1,3 +1,3 @@
-      SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz,  &
+      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
      &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
@@ -15,5 +15,8 @@
 #include "iniprint.h"
 
+      INTEGER itap
+
       INTEGER ngrid,klev
+      REAL ptimestep
       REAL ztv(ngrid,klev)
       REAL zthl(ngrid,klev)
@@ -32,13 +35,25 @@
       INTEGER lalim(ngrid)
       integer lev_out                           ! niveau pour les print
+      real zcon2(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)
@@ -55,4 +70,5 @@
       REAL linter(ngrid)
       INTEGER lmix(ngrid)
+      INTEGER lmix_bis(ngrid)
       REAL    wmaxa(ngrid)
 
@@ -85,4 +101,9 @@
             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 
@@ -91,4 +112,15 @@
            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.
@@ -117,11 +149,23 @@
       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))  &
@@ -129,4 +173,6 @@
                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
@@ -147,5 +193,19 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+
+
+! 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) &
@@ -153,4 +213,5 @@
      &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
      &      *(zlev(ig,3)-zlev(ig,2))
+       endif
 
 
@@ -160,4 +221,5 @@
 !
 !test:estimation de ztva_new_est sans entrainement
+
                Tbef=ztla(ig,l-1)*zpspsk(ig,l)
                zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
@@ -204,14 +266,43 @@
 !
 !calcul du detrainement
+!=======================
+
+
+! 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)))
-       if (prt_level.ge.20) print*,'coucou calcul detr 1'
+             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))*  &
@@ -222,14 +313,108 @@
      &       *((zmax_sec(ig)-zlev(ig,l))/  &
      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
-        if (prt_level.ge.20) print*,'coucou calcul detr 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))
-        if (prt_level.ge.20) print*,'coucou calcul detr 3'
+             detr_starc(ig,l)=detr_star(ig,l)
+
+        if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
              
           endif
-          detr_star(ig,l)=detr_star(ig,l)/f0(ig)
+
+#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
@@ -237,4 +422,39 @@
 ! 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)
@@ -249,9 +469,25 @@
 ! 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
 
@@ -261,4 +497,5 @@
 
            detr_star(ig,l)=f_star(ig,l)
+#ifdef pre_int1c
            if (l.gt.lalim(ig)+1) then
                entr_star(ig,l)=0.
@@ -269,12 +506,20 @@
                linter(ig)=l+1
             else
-               entr_star(ig,l)=detr_star(ig,l)
+               entr_star(ig,l)=0.4*detr_star(ig,l)
             endif
+#else
+           entr_star(ig,l)=0.4*detr_star(ig,l)
+#endif
         endif
-
-      else
+#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
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -292,4 +537,5 @@
 
 !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 
 !----------------------------------------------------------------------------
@@ -336,4 +582,5 @@
               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
@@ -356,4 +603,5 @@
             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 
@@ -383,4 +631,19 @@
       enddo
 
+        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
+#ifdef troisD
+      call wrgradsfi(1,klev,zqla_est,'ql_es     ','ql_es     ')
+      call wrgradsfi(1,klev,entr_star1,'e_st1     ','e_st1     ')
+      call wrgradsfi(1,klev,entr_star2,'e_st2     ','e_st2     ')
+      call wrgradsfi(1,klev,detr_stara,'d_sta     ','d_sta     ')
+      call wrgradsfi(1,klev,detr_starb,'d_stb     ','d_stb     ')
+      call wrgradsfi(1,klev,detr_starc,'d_stc     ','d_stc     ')
+      call wrgradsfi(1,klev,zqla0     ,'zqla0     ','zqla0     ')
+      call wrgradsfi(1,klev,detr_star0,'d_st0     ','d_st0     ')
+      call wrgradsfi(1,klev,detr_star1,'d_st1     ','d_st1     ')
+      call wrgradsfi(1,klev,detr_star2,'d_st2     ','d_st2     ')
+#endif
+
+!     print*,'thermcell_plume OK'
 
       return 
