SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
     &                  lalim,lmin,alim_star,alim_star_tot,lev_out)

!----------------------------------------------------------------------
!thermcell_init: calcul du profil d alimentation du thermique
!----------------------------------------------------------------------
      IMPLICIT NONE
#include "iniprint.h"
#include "thermcell.h"

      INTEGER l,ig
!arguments d entree
      INTEGER ngrid,nlay
      REAL ztv(ngrid,nlay)
      REAL zlay(ngrid,nlay)
      REAL zlev(ngrid,nlay+1)
!arguments de sortie
      INTEGER lalim(ngrid)
      INTEGER lmin(ngrid)
      REAL alim_star(ngrid,nlay)
      REAL alim_star_tot(ngrid)
      integer lev_out                           ! niveau pour les print
      
      REAL zzalim(ngrid)
!CR: ponderation entrainement des couches instables
!def des alim_star tels que alim=f*alim_star      

      do l=1,nlay
         do ig=1,ngrid 
            alim_star(ig,l)=0.
         enddo
      enddo
! determination de la longueur de la couche d entrainement
      do ig=1,ngrid
         lalim(ig)=1
      enddo

      if (iflag_thermals_ed.ge.1) then
!si la première couche est instable, on declenche un thermique
         do ig=1,ngrid
            if (ztv(ig,1).gt.ztv(ig,2)) then
               lmin(ig)=1
               lalim(ig)=2
               alim_star(ig,1)=1.
               alim_star_tot(ig)=alim_star(ig,1)
               print*,'init',alim_star(ig,1),alim_star_tot(ig)
            else
                lmin(ig)=1
                lalim(ig)=1
                alim_star(ig,1)=0.
                alim_star_tot(ig)=0. 
            endif
         enddo
 
         else
!else iflag_thermals_ed=0 ancienne def de l alim 

!on ne considere que les premieres couches instables
      do l=nlay-2,1,-1
         do ig=1,ngrid
            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
     &          ztv(ig,l+1).le.ztv(ig,l+2)) then
               lalim(ig)=l+1
            endif
          enddo
      enddo

! determination du lmin: couche d ou provient le thermique

      do ig=1,ngrid
! FH initialisation de lmin a nlay plutot que 1.
!        lmin(ig)=nlay
         lmin(ig)=1
      enddo
      do l=nlay,2,-1
         do ig=1,ngrid
            if (ztv(ig,l-1).gt.ztv(ig,l)) then
               lmin(ig)=l-1
            endif
         enddo
      enddo
!
      zzalim(:)=0.
      do l=1,nlay-1
         do ig=1,ngrid 
             if (l<lalim(ig)) then
                zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
             endif
          enddo
      enddo
      do ig=1,ngrid
          if (lalim(ig)>1) then
             zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
          else
             zzalim(ig)=zlay(ig,1)
          endif
      enddo

      if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))

! definition de l'entrainement des couches
      if (1.eq.1) then
      do l=1,nlay-1
         do ig=1,ngrid 
            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
!def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
             alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     &                       *sqrt(zlev(ig,l+1)) 
            endif
         enddo
      enddo
      else
      do l=1,nlay-1
         do ig=1,ngrid
            if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
     &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
             alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
     &        *(zlev(ig,l+1)-zlev(ig,l))
            endif
         enddo
      enddo
      endif
      
! pas de thermique si couche 1 stable
      do ig=1,ngrid
!CRnouveau test
        if (alim_star(ig,1).lt.1.e-10) then 
            do l=1,nlay
                alim_star(ig,l)=0.
            enddo
            lmin(ig)=1
         endif
      enddo 
! calcul de l alimentation totale
      do ig=1,ngrid
         alim_star_tot(ig)=0.
      enddo
      do l=1,nlay
         do ig=1,ngrid
            alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
         enddo
      enddo
!
! Calcul entrainement normalise
      do l=1,nlay 
         do ig=1,ngrid 
            if (alim_star_tot(ig).gt.1.e-10) then
               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
            endif
         enddo
      enddo
       
!on remet alim_star_tot a 1
      do ig=1,ngrid 
         alim_star_tot(ig)=1.
      enddo

      endif
!endif iflag_thermals_ed
      return 
      end