!
!
!
      SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,             &
      &                            lalim,alim_star,f_star,                    &
      &                            zmax,wmax,f,lev_out)
      
!==============================================================================
!thermcell_closure: fermeture, determination de f
!
! Modification 7 septembre 2009
! 1. On enleve alim_star_tot des arguments pour le recalculer et etre ainis
! coherent avec l'integrale au numerateur.
! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec
! l'idee etant que le choix se fasse a l'appel de thermcell_closure
! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
!==============================================================================
      
      USE thermcell_mod
      
      IMPLICIT NONE
      
      
!=============================================================================
! Declaration
!=============================================================================
      
!      inputs:
!      -------
      
      INTEGER ngrid,nlay
      INTEGER lalim(ngrid)
      INTEGER lev_out                           ! niveau pour les print
      
      REAL alim_star(ngrid,nlay)
      REAL f_star(ngrid,nlay+1)
      REAL rho(ngrid,nlay)
      REAL zlev(ngrid,nlay)
      REAL zmax(ngrid)
      REAL wmax(ngrid)
      REAL ptimestep
      
!      outputs:
!      --------
      
      REAL f(ngrid)
      
!      local:
!      ------
      
      INTEGER ig,k 
      INTEGER llmax
      
      REAL alim_star_tot(ngrid)
      REAL alim_star2(ngrid)
      REAL zdenom(ngrid)
      
!==============================================================================
! Initialization
!==============================================================================
      
      alim_star2(:) = 0.
      alim_star_tot(:) = 0.
      
      f(:) = 0.
      
      llmax = 1
      
!==============================================================================
! Closure
!==============================================================================
      
!------------------------------------------------------------------------------
! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
!------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         IF (lalim(ig)>llmax) THEN
            llmax = lalim(ig)
         ENDIF
      ENDDO
      
!------------------------------------------------------------------------------
! Calcul des integrales sur la verticale de alim_star et de alim_star^2/(rho dz)
!------------------------------------------------------------------------------
      
      DO k=1,llmax-1
         DO ig=1,ngrid
            IF (k<lalim(ig)) THEN
               alim_star2(ig) = alim_star2(ig) + alim_star(ig,k)**2           &
               &              / (rho(ig,k) * (zlev(ig,k+1) - zlev(ig,k)))
               alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,k)
            ENDIF
         ENDDO
      ENDDO
      
      DO ig=1,ngrid
         IF (alim_star2(ig)>1.e-10) THEN
            f(ig) = wmax(ig) * alim_star_tot(ig)                              &
            &     / (max(500.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
         ENDIF
      ENDDO
      
RETURN
END