!
!
!
SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
                          lmin,lmax,entr_star,detr_star,                      &
                          f,rhobarz,zlev,zw2,fm,entr,detr)
      
      
!===============================================================================
!  Purpose: deduction des flux
!  
!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
!  
!===============================================================================
      
      USE print_control_mod, ONLY: prt_level
      USE thermcell_mod, ONLY: fomass_max, alpha_max
      
      IMPLICIT NONE
      
      
!===============================================================================
! Declaration
!===============================================================================
      
!     Inputs:
!     ------- 
      
      INTEGER ngrid, nlay
      INTEGER lmin(ngrid)
      INTEGER lmax(ngrid)
      
      REAL entr_star(ngrid,nlay)
      REAL detr_star(ngrid,nlay)
      REAL zw2(ngrid,nlay+1)
      REAL zlev(ngrid,nlay+1)
      REAL masse(ngrid,nlay)
      REAL ptimestep
      REAL rhobarz(ngrid,nlay)
      REAL f(ngrid)
      
!     Outputs:
!     --------      
      
      REAL entr(ngrid,nlay)
      REAL detr(ngrid,nlay)
      REAL fm(ngrid,nlay+1)
      
!     Local:
!     ------
      
      INTEGER ig, l, k
      INTEGER igout, lout                 ! Error grid point and level
      
      REAL zfm                            ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max
      REAL f_old                          ! Save initial value of mass flux
      REAL eee0                           ! Save initial value of entrainment
      REAL ddd0                           ! Save initial value of detrainment
      REAL eee                            ! eee0 - layer mass * maximal authorized mass fraction / dt
      REAL ddd                            ! ddd0 - eee
      REAL zzz                            ! Temporary variable set to mass flux
      REAL fact
      REAL test
      
      INTEGER ncorecentr
      INTEGER ncorecdetr
      INTEGER nerrorequa
      INTEGER ncorecfact
      INTEGER ncorecalpha
      
      LOGICAL labort_physic
      
!===============================================================================
! Initialization
!===============================================================================
      
      nerrorequa = 0
      ncorecentr = 0
      ncorecdetr = 0
      ncorecfact = 0
      ncorecalpha = 0
      
      entr(:,:) = 0.
      detr(:,:) = 0.
      fm(:,:)   = 0.
      
      labort_physic = .false.
      
      fact = 0.
      
!===============================================================================
! Calcul de l'entrainement, du detrainement et du flux de masse
!===============================================================================
      
!-------------------------------------------------------------------------------
! Multiplication par la norme issue de la relation de fermeture
!-------------------------------------------------------------------------------
      
      DO l=1,nlay
         entr(:,l) = f(:) * entr_star(:,l)
         detr(:,l) = f(:) * detr_star(:,l)
      ENDDO
      
!-------------------------------------------------------------------------------
! Calcul du flux de masse
!-------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid
            IF (l < lmax(ig) .and. l >= lmin(ig)) THEN
               fm(ig,l+1) = fm(ig,l) + entr(ig,l) - detr(ig,l)
            ELSEIF (l == lmax(ig)) THEN
               fm(ig,l+1) = 0.
               detr(ig,l) = fm(ig,l) + entr(ig,l)
            ELSE
               fm(ig,l+1) = 0.
               entr(ig,l) = 0.
               detr(ig,l) = 0.
            ENDIF
         ENDDO
      ENDDO
      
!===============================================================================
! Checking
!===============================================================================
      
      DO l=1,nlay
         
!-------------------------------------------------------------------------------
! Is incoming mass flux positive ?
!-------------------------------------------------------------------------------
         
         DO ig=1,ngrid
            IF (fm(ig,l) < 0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
!-------------------------------------------------------------------------------
! Is entrainment positive ?
!-------------------------------------------------------------------------------
         
         DO ig=1,ngrid
            IF (entr(ig,l) < 0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
!-------------------------------------------------------------------------------
! Is detrainment positive ?
!-------------------------------------------------------------------------------
         
         DO ig=1,ngrid
            IF (detr(ig,l) < 0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
!-------------------------------------------------------------------------------
! Abort
!-------------------------------------------------------------------------------
         
         IF (labort_physic) THEN
            print *, '---------------------------------------------------------'
            print *, 'ERROR: mass flux has negative value(s)!'
            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
            print *, 'lmin,lmax', lmin(igout), lmax(igout)
            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
            DO k=nlay,1,-1
               print *, 'fm ', fm(igout,k+1)
               print *, 'entr,detr', entr(igout,k), detr(igout,k)
            ENDDO
            print *, 'fm ', fm(igout,1)
            print *, '---------------------------------------------------------'
            CALL abort
         ENDIF
         
!-------------------------------------------------------------------------------
! Is detrainment lessser than incoming mass flux ?
!-------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : Even if fm has no negative value, it can be lesser than detr.
!      That's not suitable because when we will mix the plume with the
!      environment, it will detrain more mass than it is physically able to do.
!      When it occures, that imply that entr + fm is greater than detr,
!      otherwise we get a negative outgoing mass flux (cf. above).
!      That is why we correct the detrainment as follows.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         DO ig=1,ngrid
            IF (detr(ig,l).gt.fm(ig,l)) THEN
               detr(ig,l) = fm(ig,l)
               entr(ig,l) = fm(ig,l+1)
               ncorecdetr  = ncorecdetr + 1
            ENDIF
         ENDDO
         
!-------------------------------------------------------------------------------
! Is entrained mass lesser than fomass_max ?
!-------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : Entrainment is bigger than the maximal authorized value.
!      If we consider that the excess entrainement is in fact plume air which
!      is not detrained then we compensate it by decreasing detr.
!      If it's not enough, we can increase entr in the layer above and decrease
!      the outgoing mass flux in the current layer.
!      If it's still unsufficient, code will abort (now commented).
!      Else we reset entr to its intial value and we renormalize entrainment,
!      detrainment and mass flux profiles such as we do not exceed the maximal
!      authorized entrained mass.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         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 > 0.) THEN
               entr(ig,l) = entr(ig,l) - eee
               ncorecentr  = ncorecentr + 1
               IF (ddd > 0.) THEN
                  detr(ig,l) = ddd
               ELSEIF (l == lmax(ig)) THEN
                  detr(ig,l) = fm(ig,l) + entr(ig,l)
               ELSEIF (entr(ig,l+1) > ABS(ddd)) THEN
                  detr(ig,l) = 0.
                  fm(ig,l+1) = fm(ig,l) + entr(ig,l)
                  entr(ig,l+1) = entr(ig,l+1) + ddd
!               ELSE
!                  entr(ig,l) = entr(ig,l) + eee
!                  igout = ig
!                  lout = l
!                  labort_physic = .true.
               ELSE
                  fact = max(fact, eee0 / (eee0 - eee))
                  entr(ig,l) = eee0
                  ncorecfact = ncorecfact + 1
               ENDIF
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, '---------------------------------------------------------'
            print *, 'ERROR: Entrainment is greater than maximal authorized value!'
            print *, '       Nor detrainment neither entrainment can compensate it!'
            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
            print *, 'lmin,lmax', lmin(igout), lmax(igout)
            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
            print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep
            print *, '   fomass_max :', fomass_max
            print *, '   masse      :', masse(igout,lout)
            print *, '   ptimestep  :', ptimestep
            print *, 'norm  :', f(igout)
            print *, 'entr* :', entr_star(igout,lout)
            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
            DO k=nlay,1,-1
               print *, 'fm ', fm(igout,k+1)
               print *, 'entr,detr', entr(igout,k), detr(igout,k)
            ENDDO
            print *, 'fm ', fm(igout,1)
            print *, '---------------------------------------------------------'
            CALL abort
         ENDIF
         
!-------------------------------------------------------------------------------
! Is Updraft fraction lesser than alpha_max ?
!-------------------------------------------------------------------------------
         
         DO ig=1,ngrid
            zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max
            
            IF (fm(ig,l+1) > zfm) THEN
               f_old = fm(ig,l+1)
               fm(ig,l+1) = zfm
               detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
               ncorecalpha = ncorecalpha + 1
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : The previous change doesn't observe the equation df/dz = e - d.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
               entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
            ENDIF
         ENDDO
         
      ENDDO
      
      IF (fact > 0.) THEN
         entr(:,:) = entr(:,:) / fact
         detr(:,:) = detr(:,:) / fact
         fm(:,:) = fm(:,:) / fact
      ENDIF
      
!-------------------------------------------------------------------------------
! Is equation df/dz = e - d still verified ?
!-------------------------------------------------------------------------------
      
!      DO l=1,nlay
!         DO ig=1,ngrid
!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
!            IF (test > 1.e-10) THEN
!               nerrorequa = nerrorequa + 1
!            ENDIF
!         ENDDO
!      ENDDO
      
!-------------------------------------------------------------------------------
! Reset top boundary conditions
!-------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         IF (lmax(ig).GT.0) THEN
            detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig))
            fm(ig,lmax(ig)+1) = 0.
            entr(ig,lmax(ig)) = 0.
         ENDIF
      ENDDO
      
!===============================================================================
! Outputs
!===============================================================================
      
      IF (prt_level > 0) THEN
         
         IF (ncorecdetr.ge.1) THEN
            print *, 'WARNING: Detrainment is greater than mass flux!'
            print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid
         ENDIF
         
         IF (ncorecentr.ge.1) THEN
            print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
            print *, 'in', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid
         ENDIF
         
         IF (ncorecfact.ge.1) THEN
            print *, 'WARNING: Entrained mass needs renormalization to be ok!'
            print *, 'in', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid
         ENDIF
         
!         IF (nerrorequa.ge.1) THEN
!            print *, 'WARNING: !'
!            print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid
!         ENDIF
         
         IF (ncorecalpha.ge.1) THEN
            print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
            print *, 'in', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid
         ENDIF
         
      ENDIF
      
      
RETURN
END