SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz,  &
     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
     &           lentr,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)

!--------------------------------------------------------------------------
!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
!--------------------------------------------------------------------------

      IMPLICIT NONE

#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"

      INTEGER ngrid,klev
      REAL ztv(ngrid,klev)
      REAL zthl(ngrid,klev)
      REAL po(ngrid,klev)
      REAL zl(ngrid,klev)
      REAL rhobarz(ngrid,klev)
      REAL zlev(ngrid,klev+1)
      REAL pplev(ngrid,klev+1)
      REAL pphi(ngrid,klev)
      REAL zpspsk(ngrid,klev)
      REAL alim_star(ngrid,klev)
      REAL zmax_sec(ngrid)
      REAL f0(ngrid)
      REAL l_mix
      REAL r_aspect
      INTEGER lentr(ngrid)
      integer lev_out                           ! niveau pour les print

      REAL ztva(ngrid,klev)
      REAL ztla(ngrid,klev)
      REAL zqla(ngrid,klev)
      REAL zqta(ngrid,klev)
      REAL zha(ngrid,klev)

      REAL detr_star(ngrid,klev)
      REAL entr_star(ngrid,klev)
      REAL detr(ngrid,klev)
      REAL entr(ngrid,klev)

      REAL zw2(ngrid,klev+1)
      REAL w_est(ngrid,klev+1)
      REAL f_star(ngrid,klev+1)
      REAL wa_moy(ngrid,klev+1)

      REAL ztva_est(ngrid,klev)
      REAL zqla_est(ngrid,klev)
      REAL zqsatth(ngrid,klev)

      REAL linter(ngrid)
      INTEGER lmix(ngrid)
      REAL    wmaxa(ngrid)

      INTEGER ig,l,k

      real zcor,zdelta,zcvm5,qlbef
      real Tbef,qsatbef
      real dqsat_dT,DT,num,denom
      REAL REPS,RLvCp,DDT0
      PARAMETER (DDT0=.01)
      logical Zsat

      Zsat=.false.
! Initialisation
      RLvCp = RLVTT/RCPD
     
      do l=1,klev
         do ig=1,ngrid
            zqla_est(ig,l)=0.
            ztva_est(ig,l)=ztva(ig,l)
            zqsatth(ig,l)=0.
         enddo
      enddo

!AM:initialisations du thermique
      do k=1,klev
         do ig=1,ngrid
            ztva(ig,k)=ztv(ig,k)
            ztla(ig,k)=zthl(ig,k)
            zqla(ig,k)=0.
            zqta(ig,k)=po(ig,k)
         enddo
      enddo 
      do k=1,klev
        do ig=1,ngrid
           detr_star(ig,k)=0.
           entr_star(ig,k)=0.
           detr(ig,k)=0.
           entr(ig,k)=0.
        enddo
      enddo
      if (lev_out.ge.1) print*,'7 OK convect8'
      do k=1,klev+1
         do ig=1,ngrid
            zw2(ig,k)=0.
            w_est(ig,k)=0.
            f_star(ig,k)=0.
            wa_moy(ig,k)=0.
         enddo
      enddo

      if (lev_out.ge.1) print*,'8 OK convect8'
      do ig=1,ngrid
         linter(ig)=1.
         lmix(ig)=1
         wmaxa(ig)=0.
      enddo

!-----------------------------------------------------------------------------------
!boucle de calcul de la vitesse verticale dans le thermique
!-----------------------------------------------------------------------------------
      do l=1,klev-1
         do ig=1,ngrid
            if (ztv(ig,l).gt.ztv(ig,l+1)  &
     &         .and.alim_star(ig,l).gt.1.e-10  &
     &         .and.zw2(ig,l).lt.1e-10) then
               ztla(ig,l)=zthl(ig,l) 
               zqta(ig,l)=po(ig,l)
               zqla(ig,l)=zl(ig,l)
               f_star(ig,l+1)=alim_star(ig,l)
               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
     &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
               w_est(ig,l+1)=zw2(ig,l+1)
!
            else if ((zw2(ig,l).ge.1e-10).and.  &
     &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
!estimation du detrainement a partir de la geometrie du pas precedent
!tests sur la definition du detr
!calcul de detr_star et entr_star
                w_est(ig,3)=zw2(ig,2)*  &
     &                   ((f_star(ig,2))**2)  &
     &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
     &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     &                   *(zlev(ig,3)-zlev(ig,2))
      if (l.gt.2) then
!---------------------------------------------------------------------------
!calcul de l entrainement et du detrainement lateral
!---------------------------------------------------------------------------
!
!test:estimation de ztva_new_est sans entrainement
               Tbef=ztla(ig,l-1)*zpspsk(ig,l)
               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
               qsatbef=MIN(0.5,qsatbef)
               zcor=1./(1.-retv*qsatbef)
               qsatbef=qsatbef*zcor
               Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
               if (Zsat) then
               qlbef=max(0.,zqta(ig,l-1)-qsatbef)
               DT = 0.5*RLvCp*qlbef
               do while (abs(DT).gt.DDT0)
                 Tbef=Tbef+DT
                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
                 qsatbef=MIN(0.5,qsatbef)
                 zcor=1./(1.-retv*qsatbef)
                 qsatbef=qsatbef*zcor
                 qlbef=zqta(ig,l-1)-qsatbef

                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
                 zcor=1./(1.-retv*qsatbef)
                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
                 num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
                 denom=1.+RLvCp*dqsat_dT
                 DT=num/denom
               enddo
                 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef) 
               endif
        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
     &      -zqla_est(ig,l))-zqla_est(ig,l))

             w_est(ig,l+1)=zw2(ig,l)*  &
     &                   ((f_star(ig,l))**2)  &
     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
     &                   *(zlev(ig,l+1)-zlev(ig,l))
             if (w_est(ig,l+1).lt.0.) then
                w_est(ig,l+1)=zw2(ig,l)
             endif
!
!calcul du detrainement
          if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
     &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
     &       (zqla(ig,l-1).lt.1.e-10)) then 
             detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
     &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
     &       /(r_aspect*zmax_sec(ig)))
       if (lev_out.ge.20) print*,'coucou calcul detr 1'
          else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
     &            (zqla(ig,l-1).lt.1.e-10)) then
             detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
     &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
     &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
     &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
     &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
     &       *((zmax_sec(ig)-zlev(ig,l))/  &
     &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
        if (lev_out.ge.20) print*,'coucou calcul detr 2'
          else
             detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
     &                      *(zlev(ig,l+1)-zlev(ig,l))
        if (lev_out.ge.20) print*,'coucou calcul detr 3'
             
          endif
          detr_star(ig,l)=detr_star(ig,l)/f0(ig)
!
!calcul de entr_star
!
        if (detr_star(ig,l).gt.f_star(ig,l)) then
           detr_star(ig,l)=f_star(ig,l)
!a decommenter ou pas?
!           entr_star(ig,l)=0.
        endif

! Deplacement du calcul de entr_star pour eviter d'avoir aussi
! entr_star > fstar.
! FH
          entr_star(ig,l)=0.4*detr_star(ig,l)
!
      else
         detr_star(ig,l)=0.
         entr_star(ig,l)=0.
      endif
!pas d entrainement dans la couche alim
      if ((l.lt.lentr(ig))) then
           entr_star(ig,l)=0.
      endif
!
!prise en compte du detrainement et de l entrainement dans le calcul du flux
      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
     &              -detr_star(ig,l)

!test sur le signe de f_star
       if (f_star(ig,l+1).gt.1.e-10) then 
!----------------------------------------------------------------------------
!calcul de la vitesse verticale en melangeant Tl et qt du thermique
!---------------------------------------------------------------------------
!
       Zsat=.false.
       ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
     &            /(f_star(ig,l+1)+detr_star(ig,l))
!
       zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
     &            /(f_star(ig,l+1)+detr_star(ig,l))
!  
               Tbef=ztla(ig,l)*zpspsk(ig,l)
               zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
               qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
               qsatbef=MIN(0.5,qsatbef)
               zcor=1./(1.-retv*qsatbef)
               qsatbef=qsatbef*zcor
               Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
               if (Zsat) then
               qlbef=max(0.,zqta(ig,l)-qsatbef)
               DT = 0.5*RLvCp*qlbef
               do while (abs(DT).gt.DDT0)
                 Tbef=Tbef+DT
                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
                 qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
                 qsatbef=MIN(0.5,qsatbef)
                 zcor=1./(1.-retv*qsatbef)
                 qsatbef=qsatbef*zcor
                 qlbef=zqta(ig,l)-qsatbef

                 zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
                 zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
                 zcor=1./(1.-retv*qsatbef)
                 dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
                 num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
                 denom=1.+RLvCp*dqsat_dT
                 DT=num/denom
              enddo
                 zqla(ig,l) = max(0.,qlbef) 
              endif
!    
! on ecrit de maniere conservative (sat ou non)
!          T = Tl +Lv/Cp ql
           ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l)
           ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l)
!on rajoute le calcul de zha pour diagnostiques (temp potentielle)
           zha(ig,l) = ztva(ig,l)
           ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l)  &
     &              -zqla(ig,l))-zqla(ig,l))

!on ecrit zqsat 
           zqsatth(ig,l)=qsatbef  
!calcul de vitesse
           zw2(ig,l+1)=zw2(ig,l)*  &
     &                 ((f_star(ig,l))**2)  &
     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
     &                 *(zlev(ig,l+1)-zlev(ig,l))
               
            endif
        endif
!
!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 
        if (zw2(ig,l+1).lt.0.) then 
           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
           zw2(ig,l+1)=0.
        else
           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 
        endif
        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
!   lmix est le niveau de la couche ou w (wa_moy) est maximum
            lmix(ig)=l+1
            wmaxa(ig)=wa_moy(ig,l+1)
        endif
        enddo
      enddo


      return 
      end