MODULE lmdz_thermcell_height
CONTAINS

      SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lcong,lintercong,lmix,  &
             zw2,zlev,lmax,zmax,zmax0,zmix,wmax,zcong)
      IMPLICIT NONE

!-----------------------------------------------------------------------------
!thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
!-----------------------------------------------------------------------------

! arguments

! Entree
      INTEGER, INTENT(IN) :: ngrid,nlay
      REAL, INTENT(IN), DIMENSION(ngrid) :: linter,lintercong
      REAL, INTENT(IN), DIMENSION(ngrid,nlay+1) :: zlev
! Sortie
      REAL, INTENT(OUT), DIMENSION(ngrid) :: wmax,zmax,zmax0,zmix,zcong
      INTEGER, INTENT(OUT), DIMENSION(ngrid) :: lmax,lcong
! Les deux
     INTEGER, INTENT(INOUT), DIMENSION(ngrid) :: lmix,lalim,lmin
     REAL, INTENT(INOUT), DIMENSION(ngrid,nlay+1) :: zw2

! local
     REAL, DIMENSION(ngrid) :: num,denom,zlevinter,zlevintercong
     INTEGER ig,l

!calcul de la hauteur max du thermique
      do ig=1,ngrid
         lmax(ig)=lalim(ig)
      enddo
      do ig=1,ngrid
         do l=nlay,lalim(ig)+1,-1
            IF (zw2(ig,l)<=1.e-10) THEN
               lmax(ig)=l-1
            endif
         enddo
      enddo

! On traite le cas particulier qu'il faudrait éviter ou le thermique
! atteind le haut du modele ...
      do ig=1,ngrid
      IF ( zw2(ig,nlay) > 1.e-10 ) THEN
          PRINT*,'WARNING !!!!! W2 thermiques non nul derniere couche '
          lmax(ig)=nlay
      endif
      enddo

! pas de thermique si couche 1 stable
      do ig=1,ngrid
         IF (lmin(ig)>1) THEN
             lmax(ig)=1
             lmin(ig)=1
             lalim(ig)=1
         endif
      enddo 

! Determination de zw2 max
      do ig=1,ngrid
         wmax(ig)=0.
      enddo

      do l=1,nlay
         do ig=1,ngrid
            IF (l<=lmax(ig)) THEN
                IF (zw2(ig,l)<0.)THEN
                  PRINT*,'pb2 zw2<0'
                endif
                zw2(ig,l)=sqrt(zw2(ig,l))
                wmax(ig)=max(wmax(ig),zw2(ig,l))
            else
                 zw2(ig,l)=0.
            endif
          enddo
      enddo

!   Longueur caracteristique correspondant a la hauteur des thermiques.
      do  ig=1,ngrid
         zmax(ig)=0.
         zlevinter(ig)=zlev(ig,1)
      enddo

!     if (iflag_thermals_ed.ge.1) THEN
      IF (1==0) THEN
!CR:date de quand le calcul du zmax continu etait buggue 
         num(:)=0.
         denom(:)=0.
         do ig=1,ngrid
          do l=1,nlay
             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
          enddo
       enddo
       do ig=1,ngrid
       IF (denom(ig)>1.e-10) THEN
          zmax(ig)=2.*num(ig)/denom(ig)
          zmax0(ig)=zmax(ig)
       endif 
       enddo
 
      else
!CR:Calcul de zmax continu via le linter      
      do  ig=1,ngrid
! calcul de zlevinter
          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)))
!pour le cas ou on prend tjs lmin=1
!       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
       zmax0(ig)=zmax(ig)
      enddo

!CR:calcul du zcong
      do  ig=1,ngrid
! calcul de zlevintercong
         zlevintercong(ig)=(zlev(ig,lcong(ig)+1)-zlev(ig,lcong(ig)))*  &
      lintercong(ig)+zlev(ig,lcong(ig))-lcong(ig)*(zlev(ig,lcong(ig)+1)  &
      -zlev(ig,lcong(ig)))
         zcong(ig)=zlevintercong(ig)-zlev(ig,1)
!         PRINT*,"calcul zcong",lcong(ig),lintercong(ig),zlevintercong(ig),zcong(ig)
      enddo

      endif
!endif iflag_thermals_ed

! def de  zmix continu (profil parabolique des vitesses)
      do ig=1,ngrid
           IF (lmix(ig)>1) THEN
! test 
              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*,'pb zmix'
              endif
          else 
              zmix(ig)=0.
          endif
!test
         IF ((zmax(ig)-zmix(ig))<=0.) THEN
            zmix(ig)=0.9*zmax(ig)
!            PRINT*,'pb zmix>zmax'
         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
END MODULE lmdz_thermcell_height