Index: LMDZ4/trunk/libf/phytherm/thermcell_flux.F90
===================================================================
--- LMDZ4/trunk/libf/phytherm/thermcell_flux.F90	(revision 814)
+++ LMDZ4/trunk/libf/phytherm/thermcell_flux.F90	(revision 852)
@@ -28,9 +28,10 @@
 
       integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
-      integer ncorecfm4,ncorecfm5,ncorecfm6
+      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
       
 
       REAL entr(ngrid,klev),detr(ngrid,klev)
       REAL fm(ngrid,klev+1)
+      REAL zfm
 
       integer igout
@@ -40,8 +41,9 @@
       REAL f_old,ddd0,eee0,ddd,eee,zzz
 
-      REAL fracmax
-      save fracmax
-
-      fracmax=0.5
+      REAL fomass_max,alphamax
+      save fomass_max,alphamax
+
+      fomass_max=0.5
+      alphamax=0.7
 
       ncorecfm1=0
@@ -51,4 +53,6 @@
       ncorecfm5=0
       ncorecfm6=0
+      ncorecfm7=0
+      ncorecfm8=0
       ncorecalpha=0
 
@@ -56,4 +60,15 @@
       fm(:,:)=0.
       
+      if (lev_out.ge.10) then
+         write(lunout,*) 'Dans thermcell_flux 0'
+         write(lunout,*) 'flux base ',f(igout)
+         write(lunout,*) 'lmax ',lmax(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) 'ig= ',igout
+         write(lunout,*) ' l E*    A*     D*  '
+         write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
+     &    ,l=1,lmax(igout))
+      endif
+
 
 !-------------------------------------------------------------------------
@@ -66,5 +81,13 @@
             if (l.le.lmax(ig)) then
                if (entr_star(ig,l).gt.1.) then
-                    print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
+                    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)
@@ -72,12 +95,4 @@
                     stop
                endif
-            else
-               if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0) then
-                    print*,'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
@@ -92,4 +107,15 @@
          detr(:,l)=f(:)*detr_star(:,l)
       enddo
+
+      if (lev_out.ge.10) then
+         write(lunout,*) 'Dans thermcell_flux 1'
+         write(lunout,*) 'flux base ',f(igout)
+         write(lunout,*) 'lmax ',lmax(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) 'ig= ',igout
+         write(lunout,*) ' l   E    D     W2'
+         write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
+     &    ,zw2(igout,l+1),l=1,lmax(igout))
+      endif
 
       fm(:,1)=0.
@@ -107,12 +133,22 @@
       enddo
 
-      if (lev_out.ge.10) then
-         write(lunout,*) 'Dans thermcell_flux 1'
-         write(lunout,*) 'flux base ',f(igout)
-         write(lunout,*) 'lmax ',lmax(igout)
-         write(lunout,*) 'lalim ',lalim(igout)
-         write(lunout,'(i6,i4,3e15.5)') (ig,l,entr(igout,l),detr(igout,l) &
-     &    ,fm(igout,l+1),l=1,lmax(igout))
-      endif
+
+
+! 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 (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'2  ')
 
 !-------------------------------------------------------------------------
@@ -130,4 +166,8 @@
          enddo
       enddo
+
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'3  ')
 
 !-------------------------------------------------------------------------
@@ -151,4 +191,8 @@
       enddo
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'4  ')
+
 
 !-------------------------------------------------------------------------
@@ -167,4 +211,8 @@
       enddo
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'5  ')
+
 !-------------------------------------------------------------------------
 !detr ne peut pas etre superieur a fm
@@ -172,4 +220,5 @@
 
       if(1.eq.1) then
+
       do l=1,klev
          do ig=1,ngrid
@@ -181,6 +230,23 @@
                ncorecfm6=ncorecfm6+1
                detr(ig,l)=fm(ig,l)
-               entr(ig,l)=fm(ig,l+1)
-            endif
+!              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.
+               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)
@@ -194,4 +260,8 @@
 
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'6  ')
+
 !-------------------------------------------------------------------------
 !fm ne peut pas etre negatif
@@ -207,5 +277,5 @@
             endif
             if (detr(ig,l).lt.0.) then
-               print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
+               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)
@@ -215,19 +285,41 @@
      enddo
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'7  ')
+
 !-----------------------------------------------------------------------
 !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.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+     if(1.eq.0) then
+
      do l=1,klev
         do ig=1,ngrid
            if (zw2(ig,l+1).gt.1.e-10) then
-           if ((((fm(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.  &
-     &      1.)) 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)=rhobarz(ig,l+1)*zw2(ig,l+1)
-              zw2(ig,l+1)=0.
-              zqla(ig,l+1)=0.
+              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))
+!             lmax(ig)=l+1
+!             zmax(ig)=zlev(ig,lmax(ig))
 !             print*,'alpha>1',l+1,lmax(ig)
               ncorecalpha=ncorecalpha+1
@@ -237,4 +329,9 @@
      enddo
 !
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'8  ')
+
+     endif
 
 !-----------------------------------------------------------------------
@@ -248,5 +345,5 @@
             eee0=entr(ig,l)
             ddd0=detr(ig,l)
-            eee=entr(ig,l)-masse(ig,l)*fracmax/ptimestep
+            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
             ddd=detr(ig,l)-eee
             if (eee.gt.0.) then
@@ -272,6 +369,6 @@
                          print*,'detr',detr(ig,l)
                          print*,'masse',masse(ig,l)
-                         print*,'fracmax',fracmax
-                         print*,'masse(ig,l)*fracmax/ptimestep',masse(ig,l)*fracmax/ptimestep
+                         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)
@@ -309,8 +406,57 @@
     &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
     &     ncorecfm6,'x fm6', &
+    &     ncorecfm7,'x fm7', &
+    &     ncorecfm8,'x fm8', &
     &     ncorecalpha,'x alpha'
       endif
 
+      if (lev_out.ge.10) &
+    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,'fin')
+
 
       return
       end
+
+      subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
+    &    ptimestep,masse,entr,detr,fm,descr)
+
+     implicit none
+
+      integer ngrid,klev,lunout,igout,l,lm
+
+      integer lmax(klev),lalim(klev)
+      real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
+      real fm(ngrid,klev+1),f(ngrid)
+
+      character*3 descr
+
+      lm=lmax(igout)+5
+      if(lm.gt.klev) lm=klev
+
+      print*,'Impression jusque lm=',lm
+
+         write(lunout,*) 'Dans thermcell_flux '//descr
+         write(lunout,*) 'flux base ',f(igout)
+         write(lunout,*) 'lmax ',lmax(igout)
+         write(lunout,*) 'lalim ',lalim(igout)
+         write(lunout,*) 'ig= ',igout
+         write(lunout,'(a3,4a14)') 'l','M','E','D','F'
+         write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, &
+     &     entr(igout,l),detr(igout,l) &
+     &    ,fm(igout,l+1),l=1,lm)
+
+
+      do l=lmax(igout)+1,klev
+          if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then
+          print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout)
+          print*,'entr(igout,l)',entr(igout,l)
+          print*,'detr(igout,l)',detr(igout,l)
+          print*,'fm(igout,l)',fm(igout,l)
+          stop
+          endif
+      enddo
+
+      return
+      end
+
