!
!
!
      SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,          &
      &                           zw2,zlev,lmax,zmax,zmax0,zmix,              &
      &                           wmax,f_star,lev_out)
      
      
!==============================================================================
! thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
!==============================================================================
      
      USE thermcell_mod
      USE print_control_mod, ONLY: prt_level
      
      IMPLICIT NONE
      
      
!==============================================================================
! Declaration
!==============================================================================
      
!   inputs:
!   -------
      
      INTEGER ngrid,nlay
      INTEGER lalim(ngrid)
      INTEGER lmin(ngrid)
      INTEGER lev_out                           ! niveau pour les print
      
      REAL zlev(ngrid,nlay+1)
      REAL f_star(ngrid,nlay+1)
      
!   outputs:
!   --------
      
      INTEGER lmix(ngrid)
      INTEGER lmax(ngrid)
      
      REAL linter(ngrid)
      REAL zmax(ngrid)
      REAL zmax0(ngrid)
      REAL zmix(ngrid)
      REAL wmax(ngrid)
      REAL zw2(ngrid,nlay+1)
      
!   local:
!   ------
      
      INTEGER ig,l
      
      REAL num(ngrid)
      REAL denom(ngrid)
      REAL zlevinter(ngrid)
      REAL zlev2(ngrid,nlay+1)
      
!==============================================================================
! Initialization
!==============================================================================
      
      DO ig=1,ngrid
         lmax(ig) = lalim(ig)
         lmix(ig) = lmin(ig)
      ENDDO
      
      DO ig=1,ngrid
         wmax(ig) = 0.
      ENDDO
      
      DO ig=1,ngrid
         zmax(ig) = 0.
         zlevinter(ig) = zlev(ig,1)
      ENDDO
      
      DO ig=1,ngrid
         zlev2(ig,:) = zlev(ig,:) - zlev(ig,lmin(ig))
      ENDDO
      
!==============================================================================
! Calcul de la hauteur max du thermique
!==============================================================================
      
!------------------------------------------------------------------------------
! Calcul de lmax
!------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         DO l=nlay,lmin(ig)+1,-1
            IF (zw2(ig,l).le.0.) THEN
               lmax(ig) = l - 1
            ELSEIF (f_star(ig,l).le.0.) THEN
               lmax(ig) = l - 1
            ENDIF
         ENDDO
      ENDDO
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! On traite le cas particulier qu'il faudrait eviter ou le thermique
! atteind le haut du modele ...
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      DO ig=1,ngrid
         IF (zw2(ig,nlay).gt.1.e-10) THEN
            print *, 'WARNING: a thermal plume reaches the highest level!'
            lmax(ig)=nlay
         ENDIF
      ENDDO
      
!------------------------------------------------------------------------------
! Calcul de zmax continu via le linter
!------------------------------------------------------------------------------
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : lmin=lmax means the plume is not active and then zw2=0 at each level.
!      Otherwise we have lmax>lmin.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      DO ig=1,ngrid
         l = lmax(ig)
         
         IF (l==lmin(ig)) THEN
            linter(ig) = l
         ELSE
            linter(ig) = l - zw2(ig,l) / (zw2(ig,l+1) - zw2(ig,l))
         ENDIF
      ENDDO
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : zmax is equal to zlevinter-zlev(lmin) because lmin can be different
!      from 1.
!      We have to substract this level because zmax must be the plume height
!      (to derive the good closure), not the plume highest altitude.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      DO  ig=1,ngrid
         zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig)   &
         &             + zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)      &
         &             - zlev(ig,lmax(ig)))
         zmax(ig) = max(zmax(ig), zlevinter(ig) - zlev(ig,lmin(ig)))
         zmax0(ig) = zmax(ig)
      ENDDO
      
!==============================================================================
! Calcul de wmax et du niveau d'inversion.
!==============================================================================
      
!------------------------------------------------------------------------------
! Calcul de lmix et wmax
!------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid
            IF(l.le.lmax(ig) .and. l.gt.lmin(ig)) THEN
               IF (zw2(ig,l).lt.0.) THEN
                  print *, 'ERROR: zw2 has negative value(s)!'
                  CALL abort
               ENDIF
               
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : WARNING zw2 becomes its square root!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
               zw2(ig,l) = sqrt(zw2(ig,l))
               
               IF (zw2(ig,l).gt.wmax(ig)) THEN
                  wmax(ig)  = zw2(ig,l)
                  lmix(ig)  = l
               ENDIF
            ELSE
               zw2(ig,l) = 0.
            ENDIF
         ENDDO
      ENDDO
      
!------------------------------------------------------------------------------
! Calcul de zmix continu (profil parabolique des vitesses)
!------------------------------------------------------------------------------
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : We substract zlev(lmin) in zmix formula because we have to
!      compare it with zmax which is zlev(linter)-zlev(lmin).
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      DO ig=1,ngrid
         IF (lmix(ig).gt.lmin(ig)) THEN
            IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))                        &
            &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))             &
            &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                   &
            &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)  &
            &        THEN
               zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))              &
               &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)    &
               &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                &
               &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))   &
               &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))           &
               &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))          &
               &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                &
               &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))        &
               &        - zlev(ig,lmin(ig))
            ELSE
               zmix(ig) = zlev(ig,lmix(ig)) - zlev(ig,lmin(ig))
               print*,'WARNING: problematic zmix value!'
            ENDIF
         ELSE
            zmix(ig) = 0.
         ENDIF
      ENDDO
      
!==============================================================================
! Verification zmax > zmix
!==============================================================================
      
      DO ig=1,ngrid
         IF ((zmax(ig)-zmix(ig)).lt.0.) THEN
            
            zmix(ig) = 0.9 * zmax(ig)
            
            IF (prt_level.ge.100) THEN
               print *, 'WARNING: we have zmix > zmax!'
            ENDIF
         ENDIF
      ENDDO
      
!------------------------------------------------------------------------------
! Calcul du nouveau lmix correspondant
!------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         DO l=1,nlay
            IF (zmix(ig).ge.zlev2(ig,l).and.zmix(ig).lt.zlev2(ig,l+1)) THEN
               lmix(ig)=l
            ENDIF
         ENDDO
      ENDDO
      
      
RETURN
END