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


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

      IMPLICIT NONE
      
      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
      

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

      integer igout
      integer lev_out
      integer lunout

      REAL f_old,ddd0,eee0,ddd,eee,zzz

      REAL fracmax
      save fracmax

      fracmax=0.5

      ncorecfm1=0
      ncorecfm2=0
      ncorecfm3=0
      ncorecfm4=0
      ncorecfm5=0
      ncorecfm6=0
      ncorecalpha=0

!initialisation
      fm(:,:)=0.
      

!-------------------------------------------------------------------------
! 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*,'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*,'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

      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

      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

!-------------------------------------------------------------------------
! 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

!-------------------------------------------------------------------------
!Test sur fraca croissant
!-------------------------------------------------------------------------

      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


!-------------------------------------------------------------------------
!test sur flux de masse croissant
!-------------------------------------------------------------------------

      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

!-------------------------------------------------------------------------
!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)
            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


!-------------------------------------------------------------------------
!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*,'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

!-----------------------------------------------------------------------
!la fraction couverte ne peut pas etre superieure a 1            
!-----------------------------------------------------------------------
     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
              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.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)*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.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*,'fracmax',fracmax
                         print*,'masse(ig,l)*fracmax/ptimestep',masse(ig,l)*fracmax/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
!-----------------------------------------------------------------------

      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', &
    &     ncorecalpha,'x alpha'
      endif


      return
      end