!
!
!
SUBROUTINE thermcell_alim(ngrid,nlay,ztv2,zlev,alim_star,lalim,lmin)
      
      
!==============================================================================
! FH : 2015/11/06
! thermcell_alim: calcule la distribution verticale de l'alimentation 
! laterale a la base des panaches thermiques
!==============================================================================
      
      IMPLICIT NONE
      
      
!==============================================================================
! Declaration
!==============================================================================
      
!      inputs:
!      -------
      
      INTEGER, INTENT(IN) :: ngrid
      INTEGER, INTENT(IN) :: nlay
      INTEGER, INTENT(IN) :: lmin(ngrid)           ! plume initial level
      
      REAL, INTENT(IN) :: ztv2(ngrid,nlay)         ! Virtual potential temperature
      REAL, INTENT(IN) :: zlev(ngrid,nlay+1)       ! levels altitude
      
!      outputs:
!      --------
      
      INTEGER, INTENT(OUT) :: lalim(ngrid)         ! alimentation maximal level
      
      REAL, INTENT(OUT) :: alim_star(ngrid,nlay)   ! Normalized alimentation
      
!      local:
!      ------
      
      INTEGER :: ig, l
      
      REAL :: alim_star_tot(ngrid)                 ! integrated alimentation
      
!==============================================================================
! Initialization
!==============================================================================
      
      lalim(:) = 1
      
      alim_star_tot(:) = 0.
      
!==============================================================================
! Alimentation computation
!==============================================================================
      
      DO ig=1,ngrid
         DO l=lmin(ig),nlay-1
            IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(ztv2(ig,lmin(ig)).ge.ztv2(ig,l))) THEN
               alim_star(ig,l) = MAX( (ztv2(ig,l) - ztv2(ig,l+1)), 0.)  &
               &               * sqrt(zlev(ig,l+1))
               lalim(ig) = l + 1
               alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
            ELSE
               alim_star(ig,l) = 0.
            ENDIF
         ENDDO
      ENDDO
      
!------------------------------------------------------------------------------
! Alimentation normalization
!------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid 
            IF (alim_star_tot(ig) > 1.e-10 ) THEN
               alim_star(ig,l) = alim_star(ig,l) / alim_star_tot(ig)
            ENDIF
         ENDDO
      ENDDO
      
      alim_star_tot(:) = 1.
      
      
RETURN
END