!
!
!
SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
                          lmin,lmax,entr_star,detr_star,                      &
                          f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
      
      
!===============================================================================
!  Purpose: fluxes deduction
!  
!  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)
      REAL zqla(ngrid,nlay)
      REAL zmax(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 number 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                           ! Factor between old norm and new norm
      
      INTEGER ncorecdetr                  ! detr > fm-   counter
      INTEGER ncorecentr                  ! entr > e_max counter
      INTEGER ncorecalpha                 ! fm > zfm     counter
      
      LOGICAL labort_physic
      
      CHARACTER (len=20) :: modname='thermcell_flux'
      CHARACTER (len=80) :: abort_message
      
!===============================================================================
! Initialization
!===============================================================================
      
      ncorecdetr  = 0
      ncorecentr  = 0
      ncorecalpha = 0
      
      entr(:,:) = 0.
      detr(:,:) = 0.
      fm(:,:)   = 0.
      
!===============================================================================
! Mass flux, entrainment and detrainment
!===============================================================================
      
!-------------------------------------------------------------------------------
! 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
      
!-------------------------------------------------------------------------------
! Mass flux and boundary conditions
!-------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid
            IF (l.lt.lmax(ig) .and. l.ge.lmin(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.
               entr(ig,l) = 0.
               detr(ig,l) = 0.
            ENDIF
         ENDDO
      ENDDO
      
!===============================================================================
! Validity tests and corrections
!===============================================================================
      
      DO l=1,nlay
         
!-------------------------------------------------------------------------------
! Is incoming mass flux positive ?
!-------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : I remove the correction and replace it by an uncompromising test.
!      According to the previous derivations, we MUST have positive mass flux
!      everywhere! Indeed, as soon as fm becomes negative, the plume stops.
!      Then the only value which can be negative is the mass flux at the top of
!      the plume. However, this value was set to zero a few lines above.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         labort_physic=.false.
         
         DO ig=1,ngrid
            IF (fm(ig,l).lt.0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: mass flux has negative value(s)!'
            print *, 'ig,l,fm',igout, lout, fm(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 *, '-------------------------------'
            abort_message = 'Negative incoming fm in thermcell_flux'
            CALL abort_physic(modname,abort_message,1)
         ENDIF
         
!-------------------------------------------------------------------------------
! Is entrainment positive ?
!-------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : According to the previous derivations, we MUST have positive entrainment
!      and detrainment everywhere!
!      Indeed, they are set to max between zero and a computed value.
!      Then tests are uncompromising.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         DO ig=1,ngrid
            IF (entr(ig,l).lt.0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: entrainment 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 *, '-------------------------------'
            abort_message = 'Negative entrainment in thermcell_flux'
            CALL abort_physic(modname,abort_message,1)
         ENDIF
         
!-------------------------------------------------------------------------------
! Is detrainment positive ?
!-------------------------------------------------------------------------------
         
         DO ig=1,ngrid
            IF (detr(ig,l).lt.0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: detrainment has negative value(s)!'
            print *, 'ig,l,detr', igout, lout, detr(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 *, '-------------------------------'
            abort_message = 'Negative detrainment in thermcell_flux'
            CALL abort_physic(modname,abort_message,1)
         ENDIF
         
!-------------------------------------------------------------------------------
! Is detrainment lesser 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)
               
               IF (prt_level.ge.10) THEN
                  print *, 'WARNING: Detrainment is modified in thermcell_flux!'
                  print *, 'ig,l,lmax', ig, l, lmax(ig)
               ENDIF
               
               IF (prt_level.ge.100) THEN
                  print *, 'fm+', fm(ig,l+1)
                  print *, 'entr,detr', entr(ig,l), detr(ig,l)
                  print *, 'fm-', fm(ig,l)
               ENDIF
               
               ncorecdetr  = ncorecdetr + 1
               
            ENDIF
         ENDDO
         
!-------------------------------------------------------------------------------
! Is entrainment mass fraction 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.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         labort_physic=.false.
         
         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
               entr(ig,l) = entr(ig,l) - eee
               ncorecentr  = ncorecentr + 1
               
               IF (prt_level.ge.10) THEN
                  print *, 'WARNING: Entrainment is modified in thermcell_flux!'
                  print *, 'ig,l,lmax', ig, l, lmax(ig)
               ENDIF
               
               IF (ddd.gt.0.) THEN                          ! detr in the current layer is reduced
                  detr(ig,l) = ddd
               ELSEIF (l.eq.lmax(ig)) THEN                  ! detr in the last layer is adjusted
                  detr(ig,l) = fm(ig,l) + entr(ig,l)
               ELSEIF (entr(ig,l+1).gt.ABS(ddd)) THEN       ! detr in the current layer is set to 0 and entr in the layer above is reduced
                  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 = (eee0 - eee) / eee0
                  entr(ig,l) = eee0
                  detr(ig,:) = detr(ig,:) * fact
                  entr(ig,:) = entr(ig,:) * fact
                  fm(ig,:)   = fm(ig,:) * fact
                  
                  IF (prt_level.ge.1) THEN
                     print *, 'WARNING: Normalisation is modified in thermcell_flux!'
                     print *, 'old f, new f :', f(ig), fact * f(ig)
                  ENDIF
               ENDIF
            ENDIF
         ENDDO
         
!         IF (labort_physic) THEN
!            print *, 'ERROR: Entrainment is greater than its 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 *, '-------------------------------'
!            abort_message = 'Entrainment is too large in thermcell_flux'
!            CALL abort_physic (modname,abort_message,1)
!         ENDIF
         
!-------------------------------------------------------------------------------
! Is updraft fraction lesser than alpha_max ?
!-------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 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 ig=1,ngrid
            zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max
            
            IF (fm(ig,l+1) .gt. 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)
               
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : The previous change doesn't observe the equation df/dz = e - d.
!      To avoid this issue, what is better to do? I choose to increase the
!      entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and
!      there are never any problems.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
               IF (l.lt.lmax(ig)) THEN
                  entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
               ENDIF
               
               IF (prt_level.ge.10) THEN
                  print *, 'Mass flux is modified in thermcell_flux'
                  print *, 'ig,l,lmax', ig, l, lmax(ig)
               ENDIF
               
               IF (prt_level.ge.100) THEN
                  print *, 'fm-', fm(ig,l)
                  print *, 'entr,detr', entr(ig,l), detr(ig,l)
                  print *, 'fm+', fm(ig,l+1)
               ENDIF
               
               ncorecalpha = ncorecalpha + 1
            ENDIF
         ENDDO
         
      ENDDO
      
!-------------------------------------------------------------------------------
! Boundary conditions
!-------------------------------------------------------------------------------
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
!      Theoretically, we always get lmax >= lmin >= linf > 0
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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
      
!-------------------------------------------------------------------------------
! Corrections display
!-------------------------------------------------------------------------------
      
      IF (prt_level.GE.1) THEN
                  
         IF (ncorecdetr.GE.1) THEN
            print *, 'WARNING: Detrainment is greater than mass flux!'
            print *, '         In', ncorecdetr, 'point(s)'
         ENDIF
         
         IF (ncorecentr.GE.1) THEN
            print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
            print *, '         In', ncorecentr, 'point(s)'
         ENDIF
         
         IF (ncorecalpha.GE.1) THEN
            print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
            print *, '         In', ncorecalpha, 'point(s)'
         ENDIF
         
      ENDIF
      
! AB : temporary test added to check the validity of eq. df/dz = e - d
!      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.gt.1.e-10) THEN
!               print *, 'WARNING: df/dz != entr - detr'
!               print *, 'ig,l', ig, l
!               print *, 'fm+  :', fm(ig,l+1)
!               print *, 'entr,detr', entr(ig,l), detr(ig,l)
!               print *, 'fm   :', fm(ig,l)
!               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
!               print *, 'fm-  :', fm(ig,l-1)
!               print *, 'err. :', test
!            ENDIF
!         ENDDO
!      ENDDO
      
      
RETURN
END