!
!
!
SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
                            zlev,zmin,zmix,zmax,wmax,f_star)
      
      
!===============================================================================
!  Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
!===============================================================================
      
      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 zmin(ngrid)                          ! Altitude of the plume bottom
      REAL zmax(ngrid)                          ! Altitude of the plume top
      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
      REAL wmax(ngrid)                          ! Maximal vertical speed
      REAL zw2(ngrid,nlay+1)                    ! Square vertical speed (becomes its square root)
      
!     Local:
!     ------
      
      INTEGER ig, l
      
      REAL zlevinter(ngrid)
      
!===============================================================================
! 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
      
!===============================================================================
! Thermal plume height
!===============================================================================
      
!-------------------------------------------------------------------------------
! Calcul de zmin
!-------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid
            zmin(ig) = zlev(ig,lmin(ig))
         ENDDO
      ENDDO
      
!-------------------------------------------------------------------------------
! Calcul de lmax
!-------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         DO l=nlay,lmin(ig)+1,-1
            IF ((zw2(ig,l) <= 0.).or.(f_star(ig,l) <= 0.)) THEN
               lmax(ig) = l - 1
            ENDIF
         ENDDO
      ENDDO
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB: Problematic case where thermal plume reaches the top of the model...
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO ig=1,ngrid
         IF (zw2(ig,nlay) > 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
      
      DO ig=1,ngrid
         zlevinter(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig))          &
         &             * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig)))
         zmax(ig) = max(zmax(ig), zlevinter(ig))
      ENDDO
      
!===============================================================================
! Thermal plume maximal speed and inversion layer
!===============================================================================
      
!-------------------------------------------------------------------------------
! Calcul de lmix et wmax
!-------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid
            IF((l <= lmax(ig)).and.(l > 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) > 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)
!-------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         IF (lmix(ig) > 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))))) > 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))))))
            ELSE
               zmix(ig) = zlev(ig,lmix(ig))
               print *, 'WARNING: problematic zmix value!'
            ENDIF
         ELSE
            zmix(ig) = 0.
         ENDIF
      ENDDO
      
!===============================================================================
! Check consistence between zmax and zmix
!===============================================================================
      
!-------------------------------------------------------------------------------
! 
!-------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         IF ((zmax(ig)-zmix(ig)) < 0.) THEN
            print *, 'WARNING: we have zmix > zmax!'
            print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig)
            zmix(ig) = 0.9 * zmax(ig)
         ENDIF
      ENDDO
      
!-------------------------------------------------------------------------------
! Calcul du nouveau lmix correspondant
!-------------------------------------------------------------------------------
      
      DO ig=1,ngrid
         DO l=1,nlay
            IF ((zmix(ig) >= zlev(ig,l)).and.(zmix(ig) < zlev(ig,l+1))) THEN
               lmix(ig) = l
            ENDIF
         ENDDO
      ENDDO
      
      
RETURN
END