!
!
!
SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
                            zlev,lmax,zmax,zmix,wmax,f_star)
      
      
!===============================================================================
!  Purpose: 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 lmin(ngrid)
      
      REAL zlev(ngrid,nlay+1)
      REAL f_star(ngrid,nlay+1)
      
!     Outputs:
!     --------
      
      INTEGER lmix(ngrid)
      INTEGER lmax(ngrid)
      
      REAL linter(ngrid)
      REAL zmax(ngrid)                          ! maximal vertical speed
      REAL zmix(ngrid)                          ! altitude of maximal vertical speed
      REAL wmax(ngrid)                          ! maximal vertical speed
      REAL zw2(ngrid,nlay+1)                    ! square vertical speed (becomes its squere root)
      
!     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) = lmin(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..or.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.0.) THEN
            print *, 'WARNING: a thermal plume reaches the highest layer!'
            print *, 'ig,l', ig, nlay
            print *, 'zw2', zw2(ig,nlay)
            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 now 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)))
      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)!'
                  print *, 'ig,l', ig, l
                  print *, 'zw2', zw2(ig,l)
                  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)
            print *, 'WARNING: we have zmix > zmax!'
            print *, 'zmax,zmix', zmax(ig), zmix(ig)
         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