SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
     &   zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)

!-------------------------------------------------------------------------
!thermcell_closure: fermeture, determination de f
!-------------------------------------------------------------------------
      IMPLICIT NONE

      INTEGER ngrid,nlay
      INTEGER ig,k       
      REAL r_aspect,ptimestep
      integer lev_out                           ! niveau pour les print

      INTEGER lalim(ngrid)
      REAL alim_star(ngrid,nlay)
      REAL rho(ngrid,nlay)
      REAL zlev(ngrid,nlay)
      REAL zmax(ngrid),zmax_sec(ngrid)
      REAL wmax(ngrid),wmax_sec(ngrid)
      real zdenom

      REAL alim_star2(ngrid)

      REAL f(ngrid)
      REAL f0(ngrid)

      do ig=1,ngrid
         alim_star2(ig)=0.
      enddo
      do ig=1,ngrid
         if (alim_star(ig,1).LT.1.e-10) then
            f(ig)=0.
         else   
             do k=1,lalim(ig)
                alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
             enddo
             zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
             if (zdenom<1.e-14) then
                print*,'ig=',ig
                print*,'alim_star2',alim_star2(ig)
                print*,'zmax',zmax(ig)
                print*,'r_aspect',r_aspect
                print*,'zdenom',zdenom
                print*,'alim_star',alim_star(ig,:)
                print*,'zmax_sec',zmax_sec(ig)
                print*,'wmax_sec',wmax_sec(ig)
                stop
             endif
             if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then 
             f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
     &             *alim_star2(ig))
             f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
     &                     zmax_sec(ig))*wmax_sec(ig))
             else
             f(ig)=wmax(ig)/zdenom
            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
     &                     zmax(ig))*wmax(ig))
             endif
         endif
         f0(ig)=f(ig)
      enddo
      if (lev_out.ge.1) print*,'apres fermeture'

!
      return
      end