!
!
!
      SUBROUTINE thermcell_alim(flag,ngrid,klev,ztv,zlev,alim_star,           &
                                lalim,lmin)
      
      
      USE thermcell_mod, ONLY: linf, d_temp
      
      IMPLICIT NONE
      
!==============================================================================
! FH : 2015/11/06
! thermcell_alim: calcule la distribution verticale de l'alimentation 
! laterale a la base des panaches thermiques
!==============================================================================
      
      
!==============================================================================
! Declaration
!==============================================================================
      
!      inputs:
!      -------
      
      INTEGER, INTENT(IN) :: ngrid
      INTEGER, INTENT(IN) :: klev
      INTEGER, INTENT(IN) :: lmin(ngrid)        ! plume initial level
      INTEGER, INTENT(IN) :: flag               ! useless here
      
      REAL, INTENT(IN) :: ztv(ngrid,klev)       ! Large scale virtual potential temperature
! Virtual potential temperature offset in linf layer
      REAL, INTENT(IN) :: zlev(ngrid,klev+1)    ! levels altitude
      
!      outputs:
!      --------
      
      INTEGER, INTENT(OUT) :: lalim(ngrid)      ! alimentation maximal level
      
      REAL, INTENT(OUT) :: alim_star(ngrid,klev)
      
!      local:
!      ------
      
      INTEGER :: ig, l
      
      REAL :: alim_star_tot(ngrid)              ! integrated alimentation
      REAL :: ztv2(ngrid,klev)                  ! ztv + d_temp
      
!==============================================================================
! Initialization
!==============================================================================
      
      lalim(:) = 1
      
      alim_star_tot(:) = 0.
      
      ztv2(:,:) = ztv(:,:)
      ztv2(:,linf) = ztv2(:,linf) + d_temp
      
!==============================================================================
! Alimentation computation
!==============================================================================
      
      DO l=lmin(ig),klev-1
         DO ig=1,ngrid
            IF ((ztv2(ig,l)>ztv2(ig,l+1)).and.(ztv2(ig,lmin(ig))>=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,klev
         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