SUBROUTINE thermcell_flux1(ngrid,klev,ptimestep,masse, &
     &       lentr,lmax,alim_star,  &
     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,alim,entr,  &
     &       detr,zqla,zmax,lev_out,lunout,igout)


!---------------------------------------------------------------------------
!thermcell_flux: deduction des flux
!---------------------------------------------------------------------------

      IMPLICIT NONE
      
      INTEGER ig,k,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 lentr(ngrid)
      REAL zqla(ngrid,klev)
      REAL zmax(ngrid)


     integer lev_out,lunout,igout


      integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
      

      REAL alim(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
      REAL fm(ngrid,klev+1)

      REAL f_old,ddd,eee

      REAL fracmax
      save fracmax

      fracmax=0.5

      ncorecfm1=0
      ncorecfm2=0
      ncorecfm3=0
      ncorecalpha=0

!initialisation
      do ig=1,ngrid
         do k=1,klev+1
            fm(ig,k)=0.
         enddo
      enddo
      
! Calcul de l alimentation
      do ig=1,ngrid 
         do k=1,klev
            alim(ig,k)=f(ig)*alim_star(ig,k)
         enddo
      enddo
!CR:test pour entrainer moins que la masse
!       do ig=1,ngrid
!          do l=1,lentr(ig)
!             if ((alim(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
!                alim(ig,l+1)=alim(ig,l+1)+alim(ig,l)
!     s                       -0.9*masse(ig,l)/ptimestep
!                alim(ig,l)=0.9*masse(ig,l)/ptimestep
!             endif
!          enddo
!       enddo
! calcul du détrainement et de l entrainement
       do ig=1,ngrid
          do k=1,klev
             detr(ig,k)=f(ig)*detr_star(ig,k)
          enddo
          do k=1,klev
             entr(ig,k)=f(ig)*entr_star(ig,k)
          enddo
       enddo

      if (lev_out.ge.10) then
         write(lunout,*) 'Dans thermcell_flux 1'
         write(lunout,*) 'flux base ',f(igout)
         write(lunout,*) 'lmax ',lmax(igout)
      endif

!
! Calcul des flux
      do ig=1,ngrid
         do l=1,lmax(ig)
!calcul du flux
            fm(ig,l+1)=fm(ig,l)+alim(ig,l)+entr(ig,l)-detr(ig,l)

      if (lev_out.ge.10.and.ig.eq.igout) then
         write(lunout,'(i6,i4,3e15.5)') ig,l,entr(igout,l)+alim(igout,l),detr(igout,l) &
     &    ,fm(igout,l+1)
      endif


            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)=alim(ig,l)+entr(ig,l)
            endif
!verification des valeurs du flux
!Test sur fraca croissant
          if ((zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
          if ( (l.ge.lentr(ig)).and.  &
     &       ((fm(ig,l+1)/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.  &
     &       (fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))))) then
             f_old=fm(ig,l+1)
             fm(ig,l+1)=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)  &
     &                          /(rhobarz(ig,l)*zw2(ig,l))
              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) 
          endif
          endif
!test sur flux de masse croissant
        if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lentr(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)
       endif
!detr ne peut pas etre superieur a fm
       if (detr(ig,l).gt.fm(ig,l)) then
               detr(ig,l)=fm(ig,l)
               entr(ig,l)=fm(ig,l+1)-alim(ig,l)
       endif
!fm ne peut pas etre negatif
       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
!la fraction couverte ne peut pas etre superieure a 1            
       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
          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.
          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
!

!-----------------------------------------------------------------------
! 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.0) then
      do l=1,klev-1
         do ig=1,ngrid
            eee=entr(ig,l)+alim(ig,l)-masse(ig,l)*fracmax/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.ge.lmax(ig)) stop'probleme dans thermcell_flux'
                   entr(ig,l+1)=entr(ig,l+1)-ddd
                   detr(ig,l)=0.
                   fm(ig,l+1)=fm(ig,l)+entr(ig,l)+alim(ig,l)
                   detr(ig,l)=0.
                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))  &
     &                     +alim(ig,lmax(ig))
      enddo

!-----------------------------------------------------------------------
! Impression du nombre de bidouilles qui ont ete necessaires
!-----------------------------------------------------------------------

      if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecalpha > 0 ) then

          print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',ncorecfm2 &
    &   ,'x fm2',ncorecfm3,'x fm3 et', ncorecalpha,'x alpha'
      endif

      stop
      return
      end