! $Id: lmdz_thermcell_plume_6A.f90 5103 2024-07-23 13:29:36Z abarral $


MODULE lmdz_thermcell_plume_6A

CONTAINS

  SUBROUTINE thermcell_plume_6A(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
          , lev_out, lunout1, igout)
    !     &           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)
    !--------------------------------------------------------------------------
    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    !--------------------------------------------------------------------------

    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
    USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell
    USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power
    USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim
    USE lmdz_thermcell_alim, ONLY: thermcell_alim
    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat

    IMPLICIT NONE

    integer, intent(in) :: itap, lev_out, lunout1, igout, ngrid, nlay
    real, intent(in) :: ptimestep
    real, intent(in), dimension(ngrid, nlay) :: ztv
    real, intent(in), dimension(ngrid, nlay) :: zthl
    real, intent(in), dimension(ngrid, nlay) :: po
    real, intent(in), dimension(ngrid, nlay) :: zl
    real, intent(in), dimension(ngrid, nlay) :: rhobarz
    real, intent(in), dimension(ngrid, nlay + 1) :: zlev
    real, intent(in), dimension(ngrid, nlay + 1) :: pplev
    real, intent(in), dimension(ngrid, nlay) :: pphi
    real, intent(in), dimension(ngrid, nlay) :: zpspsk
    real, intent(in), dimension(ngrid) :: f0

    integer, intent(out) :: lalim(ngrid)
    real, intent(out), dimension(ngrid, nlay) :: alim_star
    real, intent(out), dimension(ngrid) :: alim_star_tot
    real, intent(out), dimension(ngrid, nlay) :: detr_star
    real, intent(out), dimension(ngrid, nlay) :: entr_star
    real, intent(out), dimension(ngrid, nlay + 1) :: f_star
    real, intent(out), dimension(ngrid, nlay) :: csc
    real, intent(out), dimension(ngrid, nlay) :: ztva
    real, intent(out), dimension(ngrid, nlay) :: ztla
    real, intent(out), dimension(ngrid, nlay) :: zqla
    real, intent(out), dimension(ngrid, nlay) :: zqta
    real, intent(out), dimension(ngrid, nlay) :: zha
    real, intent(out), dimension(ngrid, nlay + 1) :: zw2
    real, intent(out), dimension(ngrid, nlay + 1) :: w_est
    real, intent(out), dimension(ngrid, nlay) :: ztva_est
    real, intent(out), dimension(ngrid, nlay) :: zqsatth
    integer, intent(out), dimension(ngrid) :: lmix
    integer, intent(out), dimension(ngrid) :: lmix_bis
    real, intent(out), dimension(ngrid) :: linter

    REAL zdw2, zdw2bis
    REAL zw2modif
    REAL zw2fact, zw2factbis
    REAL, dimension(ngrid, nlay) :: zeps

    REAL, dimension(ngrid) :: wmaxa

    INTEGER ig, l, k, lt, it, lm
    integer nbpb

    real, dimension(ngrid, nlay) :: detr
    real, dimension(ngrid, nlay) :: entr
    real, dimension(ngrid, nlay + 1) :: wa_moy
    real, dimension(ngrid, nlay) :: ztv_est
    real, dimension(ngrid) :: ztemp, zqsat
    real, dimension(ngrid, nlay) :: zqla_est
    real, dimension(ngrid, nlay) :: zta_est

    real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt
    real zdz, zalpha, zw2m
    real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam
    real zbuoybis, zdz2, zdz3, lmel, entrbis, zdzbis
    real, dimension(ngrid) :: d_temp
    real ztv1, ztv2, factinv, zinv, zlmel
    real zlmelup, zlmeldwn, zlt, zltdwn, zltup
    real atv1, atv2, btv1, btv2
    real ztv_est1, ztv_est2
    real zcor, zdelta, zcvm5, qlbef
    real zbetalpha, coefzlmel
    real eps
    logical Zsat
    LOGICAL, dimension(ngrid) :: active, activetmp
    REAL fact_gamma, fact_gamma2, fact_epsilon2
    REAL coefc
    REAL, dimension(ngrid, nlay) :: c2

    if (ngrid==1) PRINT*, 'THERMCELL PLUME MODIFIE 2014/07/11'
    Zsat = .FALSE.
    ! Initialisation

    zbetalpha = betalpha / (1. + betalpha)


    ! Initialisations des variables r?elles
    if (1==1) then
      ztva(:, :) = ztv(:, :)
      ztva_est(:, :) = ztva(:, :)
      ztv_est(:, :) = ztv(:, :)
      ztla(:, :) = zthl(:, :)
      zqta(:, :) = po(:, :)
      zqla(:, :) = 0.
      zha(:, :) = ztva(:, :)
    else
      ztva(:, :) = 0.
      ztv_est(:, :) = 0.
      ztva_est(:, :) = 0.
      ztla(:, :) = 0.
      zqta(:, :) = 0.
      zha(:, :) = 0.
    endif

    zqla_est(:, :) = 0.
    zqsatth(:, :) = 0.
    zqla(:, :) = 0.
    detr_star(:, :) = 0.
    entr_star(:, :) = 0.
    alim_star(:, :) = 0.
    alim_star_tot(:) = 0.
    csc(:, :) = 0.
    detr(:, :) = 0.
    entr(:, :) = 0.
    zw2(:, :) = 0.
    zbuoy(:, :) = 0.
    zbuoyjam(:, :) = 0.
    gamma(:, :) = 0.
    zeps(:, :) = 0.
    w_est(:, :) = 0.
    f_star(:, :) = 0.
    wa_moy(:, :) = 0.
    linter(:) = 1.
    !     linter(:)=1.
    ! Initialisation des variables entieres
    lmix(:) = 1
    lmix_bis(:) = 2
    wmaxa(:) = 0.

    ! Initialisation a 0  en cas de sortie dans replay
    zqsat(:) = 0.
    zta_est(:, :) = 0.
    zdqt(:, :) = 0.
    zdqtjam(:, :) = 0.
    c2(:, :) = 0.


    !-------------------------------------------------------------------------
    ! On ne considere comme actif que les colonnes dont les deux premieres
    ! couches sont instables.
    !-------------------------------------------------------------------------

    active(:) = ztv(:, 1)>ztv(:, 2)
    d_temp(:) = 0. ! Pour activer un contraste de temperature a la base
    ! du panache
    !  Cet appel pourrait être fait avant thermcell_plume dans thermcell_main
    CALL thermcell_alim(thermals_flag_alim, ngrid, nlay, ztv, d_temp, zlev, alim_star, lalim)

    !------------------------------------------------------------------------------
    ! Calcul dans la premiere couche
    ! On decide dans cette version que le thermique n'est actif que si la premiere
    ! couche est instable.
    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    ! dans une couche l>1
    !------------------------------------------------------------------------------
    do ig = 1, ngrid
      ! Le panache va prendre au debut les caracteristiques de l'air contenu
      ! dans cette couche.
      if (active(ig)) then
        ztla(ig, 1) = zthl(ig, 1)
        zqta(ig, 1) = po(ig, 1)
        zqla(ig, 1) = zl(ig, 1)
        !cr: attention, prise en compte de f*(1)=1
        f_star(ig, 2) = alim_star(ig, 1)
        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
                & * (zlev(ig, 2) - zlev(ig, 1))  &
                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
        w_est(ig, 2) = zw2(ig, 2)
      endif
    enddo

    !==============================================================================
    !boucle de calcul de la vitesse verticale dans le thermique
    !==============================================================================
    do l = 2, nlay - 1
      !==============================================================================


      ! On decide si le thermique est encore actif ou non
      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
      do ig = 1, ngrid
        active(ig) = active(ig) &
                &                 .and. zw2(ig, l)>1.e-10 &
                &                 .and. f_star(ig, l) + alim_star(ig, l)>1.e-10
      enddo



      !---------------------------------------------------------------------------
      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
      ! sans tenir compte du detrainement et de l'entrainement dans cette
      ! couche
      ! C'est a dire qu'on suppose
      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
      ! avant) a l'alimentation pour avoir un calcul plus propre
      !---------------------------------------------------------------------------

      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
      CALL thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))
      do ig = 1, ngrid
        !       PRINT*,'active',active(ig),ig,l
        if(active(ig)) then
          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
          zta_est(ig, l) = ztva_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))


          !Modif AJAM

          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          lmel = fact_thermals_ed_dz * zlev(ig, l)
          !        lmel=0.09*zlev(ig,l)
          zlmel = zlev(ig, l) + lmel
          zlmelup = zlmel + (zdz / 2)
          zlmeldwn = zlmel - (zdz / 2)

          lt = l + 1
          zlt = zlev(ig, lt)
          zdz3 = zlev(ig, lt + 1) - zlt
          zltdwn = zlt - zdz3 / 2
          zltup = zlt + zdz3 / 2

          !=========================================================================
          ! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus
          !=========================================================================

          !--------------------------------------------------
          if (iflag_thermals_ed<8) then
            !--------------------------------------------------
            !AJ052014: J'ai remplac?? la boucle do par un do while
            ! afin de faire moins de calcul dans la boucle
            !--------------------------------------------------
            do while (zlmelup>zltup)
              lt = lt + 1
              zlt = zlev(ig, lt)
              zdz3 = zlev(ig, lt + 1) - zlt
              zltdwn = zlt - zdz3 / 2
              zltup = zlt + zdz3 / 2
            enddo
            !--------------------------------------------------
            !AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
            ! on cherche o?? se trouve l'altitude d'inversion
            ! en calculant ztv1 (interpolation de la valeur de
            ! theta au niveau lt en utilisant les niveaux lt-1 et
            ! lt-2) et ztv2 (interpolation avec les niveaux lt+1
            ! et lt+2). Si theta r??ellement calcul??e au niveau lt
            ! comprise entre ztv1 et ztv2, alors il y a inversion
            ! et on calcule son altitude zinv en supposant que ztv(lt)
            ! est une combinaison lineaire de ztv1 et ztv2.
            ! Ensuite, on calcule la flottabilite en comparant
            ! la temperature de la couche l a celle de l'air situe
            ! l+lmel plus haut, ce qui necessite de savoir quel fraction
            ! de cet air est au-dessus ou en-dessous de l'inversion
            !--------------------------------------------------
            atv1 = (ztv(ig, lt - 1) - ztv(ig, lt - 2)) / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
            btv1 = (ztv(ig, lt - 2) * zlev(ig, lt - 1) - ztv(ig, lt - 1) * zlev(ig, lt - 2)) &
                    / (zlev(ig, lt - 1) - zlev(ig, lt - 2))
            atv2 = (ztv(ig, lt + 2) - ztv(ig, lt + 1)) / (zlev(ig, lt + 2) - zlev(ig, lt + 1))
            btv2 = (ztv(ig, lt + 1) * zlev(ig, lt + 2) - ztv(ig, lt + 2) * zlev(ig, lt + 1)) &
                    / (zlev(ig, lt + 2) - zlev(ig, lt + 1))

            ztv1 = atv1 * zlt + btv1
            ztv2 = atv2 * zlt + btv2

            if (ztv(ig, lt)>ztv1.and.ztv(ig, lt)<ztv2) then

              !--------------------------------------------------
              !AJ052014: D??calage de zinv qui est entre le haut
              !          et le bas de la couche lt
              !--------------------------------------------------
              factinv = (ztv2 - ztv(ig, lt)) / (ztv2 - ztv1)
              zinv = zltdwn + zdz3 * factinv

              if (zlmeldwn>=zinv) then
                ztv_est(ig, l) = atv2 * zlmel + btv2
                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
                        + (1. - fact_shell) * zbuoy(ig, l)
              elseif (zlmelup>=zinv) then
                ztv_est2 = atv2 * 0.5 * (zlmelup + zinv) + btv2
                ztv_est1 = atv1 * 0.5 * (zinv + zlmeldwn) + btv1
                ztv_est(ig, l) = ((zlmelup - zinv) / zdz) * ztv_est2 + ((zinv - zlmeldwn) / zdz) * ztv_est1

                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zinv) / zdz) * (ztva_est(ig, l) - &
                        ztv_est2) / ztv_est2 + ((zinv - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
                        ztv_est1) / ztv_est1) + (1. - fact_shell) * zbuoy(ig, l)

              else
                ztv_est(ig, l) = atv1 * zlmel + btv1
                zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) &
                        + (1. - fact_shell) * zbuoy(ig, l)
              endif

            else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then

              if (zlmeldwn>zltdwn) then
                zbuoyjam(ig, l) = fact_shell * RG * ((ztva_est(ig, l) - &
                        ztv(ig, lt)) / ztv(ig, lt)) + (1. - fact_shell) * zbuoy(ig, l)
              else
                zbuoyjam(ig, l) = fact_shell * RG * (((zlmelup - zltdwn) / zdz) * (ztva_est(ig, l) - &
                        ztv(ig, lt)) / ztv(ig, lt) + ((zltdwn - zlmeldwn) / zdz) * (ztva_est(ig, l) - &
                        ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + (1. - fact_shell) * zbuoy(ig, l)

              endif

              !          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
              !    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
              !    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
              !         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- &
              !    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
              !     &          po(ig,lt-1))/po(ig,lt-1))
            endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then

          else  !   if (iflag_thermals_ed.lt.8) then
            lt = l + 1
            zlt = zlev(ig, lt)
            zdz2 = zlev(ig, lt) - zlev(ig, l)

            do while (lmel>zdz2)
              lt = lt + 1
              zlt = zlev(ig, lt)
              zdz2 = zlt - zlev(ig, l)
            enddo
            zdz3 = zlev(ig, lt + 1) - zlt
            zltdwn = zlev(ig, lt) - zdz3 / 2
            zlmelup = zlmel + (zdz / 2)
            coefzlmel = Min(1., (zlmelup - zltdwn) / zdz)
            zbuoyjam(ig, l) = 1. * RG * (coefzlmel * (ztva_est(ig, l) - &
                    ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - &
                    ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l)
          endif !   if (iflag_thermals_ed.lt.8) then

          !------------------------------------------------
          !AJAM:nouveau calcul de w?
          !------------------------------------------------
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)

          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
          zdw2 = afact * zbuoy(ig, l) / fact_epsilon
          zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
          !              zdw2bis=0.5*(zdw2+zdw2bis)
          lm = Max(1, l - 2)
          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
          !             w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
          !             w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* &
          !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
          !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
          !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))

          !--------------------------------------------------
          !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
          !--------------------------------------------------
          if (iflag_thermals_ed==8) then
            ! Ancienne version
            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))

            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)

            ! Nouvelle version Arnaud
          else
            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
            !    &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))

            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)

            !             w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* &
            !    &                     (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* &
            !    &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis))



            !            w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact))

            !             w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
            !    &                      (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))

            !             w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)

          endif

          if (iflag_thermals_ed<6) then
            zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
            !              fact_epsilon=0.0005/(zalpha+0.025)**0.5
            !              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
            fact_epsilon = 0.0002 / (zalpha + 0.1)
            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
            zdw2 = afact * zbuoy(ig, l) / fact_epsilon
            zdw2bis = afact * zbuoy(ig, l - 1) / fact_epsilon
            !              w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))

            !              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))

            w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2bis) + zdw2)

          endif
          !--------------------------------------------------
          !AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
          !on fait max(0.0001,.....)
          !--------------------------------------------------

          !             if (w_est(ig,l+1).lt.0.) then
          !               w_est(ig,l+1)=zw2(ig,l)
          !                w_est(ig,l+1)=0.0001
          !             endif

        endif
      enddo


      !-------------------------------------------------
      !calcul des taux d'entrainement et de detrainement
      !-------------------------------------------------

      do ig = 1, ngrid
        if (active(ig)) then

          !          zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1)
          zw2m = w_est(ig, l + 1)
          !          zw2m=zw2(ig,l)
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
          zbuoybis = zbuoy(ig, l)
          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)


          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
          !    &     afact*zbuoybis/zw2m - fact_epsilon )

          !          entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha*  &
          !    &     afact*zbuoybis/zw2m - fact_epsilon )



          !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)

          !=========================================================================
          ! 4. Calcul de l'entrainement et du detrainement
          !=========================================================================

          !          entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
          !    &     afact*zbuoyjam(ig,l)/zw2m - fact_epsilon )
          !          entrbis=entr_star(ig,l)

          if (iflag_thermals_ed<6) then
            fact_epsilon = 0.0002 / (zalpha + 0.1)
          endif

          detr_star(ig, l) = f_star(ig, l) * zdz             &
                  * (mix0 * 0.1 / (zalpha + 0.001)               &
                          + MAX(detr_min, -afact * zbetalpha * zbuoyjam(ig, l) / zw2m   &
                                  + detr_q_coef * (zdqt(ig, l) / zw2m)**detr_q_power))

          !          detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ &
          !    &                          ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1)

          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)

          entr_star(ig, l) = f_star(ig, l) * zdz * (&
                  mix0 * 0.1 / (zalpha + 0.001)               &
                          + zbetalpha * MAX(entr_min, &
                          afact * zbuoyjam(ig, l) / zw2m - fact_epsilon))


          !          entr_star(ig,l)=f_star(ig,l)*zdz* (         &
          !    &       mix0 * 0.1 / (zalpha+0.001)               &
          !    &     + MAX(entr_min,                   &
          !    &     zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon +  &
          !    &     detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power))


          !          entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ &
          !    &                          ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1)

          !          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &
          !    &     afact*zbuoy(ig,l)/zw2m &
          !    &     - 1.*fact_epsilon)


          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
          ! alim_star et 0 sinon
          if (l<lalim(ig)) then
            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
            entr_star(ig, l) = 0.
          endif
          !        if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then
          !          alim_star(ig,l)=entrbis
          !        endif

          !        PRINT*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
          ! Calcul du flux montant normalise
          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
                  - detr_star(ig, l)

        endif
      enddo


      !============================================================================
      ! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique
      !===========================================================================

      activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10
      do ig = 1, ngrid
        if (activetmp(ig)) then
          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))

        endif
      enddo

      ztemp(:) = zpspsk(:, l) * ztla(:, l)
      CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))
      do ig = 1, ngrid
        if (activetmp(ig)) then
          ! on ecrit de maniere conservative (sat ou non)
          !          T = Tl +Lv/Cp ql
          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
          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))
          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          zdzbis = zlev(ig, l) - zlev(ig, l - 1)
          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)
          !!!!!!!          fact_epsilon=0.002
          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
          zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
          zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)
          !              zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) &
          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
          !              lm=Max(1,l-2)
          !              zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) &
          !    &              +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1))
          !            zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2)
          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ &
          !     &                   (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis))
          !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
          !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
          !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
          if (iflag_thermals_ed==8) then
            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
          else
            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)
          endif
          !            zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* &
          !    &                     (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* &
          !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis))

          if (iflag_thermals_ed<6) then
            zalpha = f0(ig) * f_star(ig, l) / sqrt(zw2(ig, l + 1)) / rhobarz(ig, l)
            !           fact_epsilon=0.0005/(zalpha+0.025)**0.5
            !           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
            fact_epsilon = 0.0002 / (zalpha + 0.1)**1
            zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
            zw2factbis = fact_epsilon * 2. * zdzbis / (1. + betalpha)
            zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
            zdw2bis = afact * zbuoy(ig, l - 1) / (fact_epsilon)

            !            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
            !    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
            !    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
            !            zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact))
            zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2bis) + zdw2)

          endif

        endif
      enddo

      if (prt_level>=20) PRINT*, 'coucou calcul detr 460: ig, l', ig, l

      !===========================================================================
      ! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
      !===========================================================================

      nbpb = 0
      do ig = 1, ngrid
        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
          !               stop'On tombe sur le cas particulier de thermcell_dry'
          !               PRINT*,'On tombe sur le cas particulier de thermcell_plume'
          nbpb = nbpb + 1
          zw2(ig, l + 1) = 0.
          linter(ig) = l + 1
        endif

        if (zw2(ig, l + 1)<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.
          !+CR:04/05/12:correction calcul linter pour calcul de zmax continu
        elseif (f_star(ig, l + 1)<0.) then
          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
          zw2(ig, l + 1) = 0.
          !fin CR:04/05/12
        endif

        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))

        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
          !on rajoute le calcul de lmix_bis
          if (zqla(ig, l)<1.e-10) then
            lmix_bis(ig) = l + 1
          endif
          lmix(ig) = l + 1
          wmaxa(ig) = wa_moy(ig, l + 1)
        endif
      enddo

      if (nbpb>0) then
        PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
      endif

      !=========================================================================
      ! FIN DE LA BOUCLE VERTICALE
    enddo
    !=========================================================================

    !on recalcule alim_star_tot
    do ig = 1, ngrid
      alim_star_tot(ig) = 0.
    enddo
    do ig = 1, ngrid
      do l = 1, lalim(ig) - 1
        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
      enddo
    enddo

    if (prt_level>=20) PRINT*, 'coucou calcul detr 470: ig, l', ig, l
    RETURN
  end


  SUBROUTINE thermcell_plume_5B(itap, ngrid, nlay, ptimestep, ztv, zthl, po, zl, rhobarz, &
          zlev, pplev, pphi, zpspsk, alim_star, alim_star_tot, &
          lalim, f0, detr_star, entr_star, f_star, csc, ztva, &
          ztla, zqla, zqta, zha, zw2, w_est, ztva_est, zqsatth, lmix, lmix_bis, linter &
          , lev_out, lunout1, igout)
    !&           ,lev_out,lunout1,igout,zbuoy,zbuoyjam)

    !--------------------------------------------------------------------------
    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    ! Version conforme a l'article de Rio et al. 2010.
    ! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin
    !--------------------------------------------------------------------------

    USE lmdz_thermcell_ini, ONLY: prt_level, fact_thermals_ed_dz, iflag_thermals_ed, RLvCP, RETV, RG
    USE lmdz_thermcell_qsat, ONLY: thermcell_qsat
    IMPLICIT NONE

    INTEGER itap
    INTEGER lunout1, igout
    INTEGER ngrid, nlay
    REAL ptimestep
    REAL ztv(ngrid, nlay)
    REAL zthl(ngrid, nlay)
    REAL, INTENT(IN) :: po(ngrid, nlay)
    REAL zl(ngrid, nlay)
    REAL rhobarz(ngrid, nlay)
    REAL zlev(ngrid, nlay + 1)
    REAL pplev(ngrid, nlay + 1)
    REAL pphi(ngrid, nlay)
    REAL zpspsk(ngrid, nlay)
    REAL alim_star(ngrid, nlay)
    REAL f0(ngrid)
    INTEGER lalim(ngrid)
    integer lev_out                           ! niveau pour les print
    integer nbpb

    real alim_star_tot(ngrid)

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

    REAL detr_star(ngrid, nlay)
    REAL coefc
    REAL entr_star(ngrid, nlay)
    REAL detr(ngrid, nlay)
    REAL entr(ngrid, nlay)

    REAL csc(ngrid, nlay)

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

    REAL ztva_est(ngrid, nlay)
    REAL zqla_est(ngrid, nlay)
    REAL zqsatth(ngrid, nlay)
    REAL zta_est(ngrid, nlay)
    REAL zbuoyjam(ngrid, nlay)
    REAL ztemp(ngrid), zqsat(ngrid)
    REAL zdw2
    REAL zw2modif
    REAL zw2fact
    REAL zeps(ngrid, nlay)

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

    INTEGER ig, l, k

    real zdz, zbuoy(ngrid, nlay), zalpha, gamma(ngrid, nlay), zdqt(ngrid, nlay), zw2m
    real zbuoybis
    real zcor, zdelta, zcvm5, qlbef, zdz2
    real betalpha, zbetalpha
    real eps, afact
    logical Zsat
    LOGICAL active(ngrid), activetmp(ngrid)
    REAL fact_gamma, fact_epsilon, fact_gamma2, fact_epsilon2
    REAL c2(ngrid, nlay)
    Zsat = .FALSE.
    ! Initialisation

    fact_epsilon = 0.002
    betalpha = 0.9
    afact = 2. / 3.

    zbetalpha = betalpha / (1. + betalpha)


    ! Initialisations des variables reeles
    if (1==1) then
      ztva(:, :) = ztv(:, :)
      ztva_est(:, :) = ztva(:, :)
      ztla(:, :) = zthl(:, :)
      zqta(:, :) = po(:, :)
      zha(:, :) = ztva(:, :)
    else
      ztva(:, :) = 0.
      ztva_est(:, :) = 0.
      ztla(:, :) = 0.
      zqta(:, :) = 0.
      zha(:, :) = 0.
    endif

    zqla_est(:, :) = 0.
    zqsatth(:, :) = 0.
    zqla(:, :) = 0.
    detr_star(:, :) = 0.
    entr_star(:, :) = 0.
    alim_star(:, :) = 0.
    alim_star_tot(:) = 0.
    csc(:, :) = 0.
    detr(:, :) = 0.
    entr(:, :) = 0.
    zw2(:, :) = 0.
    zbuoy(:, :) = 0.
    zbuoyjam(:, :) = 0.
    gamma(:, :) = 0.
    zeps(:, :) = 0.
    w_est(:, :) = 0.
    f_star(:, :) = 0.
    wa_moy(:, :) = 0.
    linter(:) = 1.
    !     linter(:)=1.
    ! Initialisation des variables entieres
    lmix(:) = 1
    lmix_bis(:) = 2
    wmaxa(:) = 0.
    lalim(:) = 1


    !-------------------------------------------------------------------------
    ! On ne considere comme actif que les colonnes dont les deux premieres
    ! couches sont instables.
    !-------------------------------------------------------------------------
    active(:) = ztv(:, 1)>ztv(:, 2)

    !-------------------------------------------------------------------------
    ! Definition de l'alimentation
    !-------------------------------------------------------------------------
    do l = 1, nlay - 1
      do ig = 1, ngrid
        if (ztv(ig, l)> ztv(ig, l + 1) .and. ztv(ig, 1)>=ztv(ig, l)) then
          alim_star(ig, l) = MAX((ztv(ig, l) - ztv(ig, l + 1)), 0.)  &
                  * sqrt(zlev(ig, l + 1))
          lalim(ig) = l + 1
          alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
        endif
      enddo
    enddo
    do l = 1, nlay
      do ig = 1, ngrid
        if (alim_star_tot(ig) > 1.e-10) then
          alim_star(ig, l) = alim_star(ig, l) / alim_star_tot(ig)
        endif
      enddo
    enddo
    alim_star_tot(:) = 1.



    !------------------------------------------------------------------------------
    ! Calcul dans la premiere couche
    ! On decide dans cette version que le thermique n'est actif que si la premiere
    ! couche est instable.
    ! Pourrait etre change si on veut que le thermiques puisse se d??clencher
    ! dans une couche l>1
    !------------------------------------------------------------------------------
    do ig = 1, ngrid
      ! Le panache va prendre au debut les caracteristiques de l'air contenu
      ! dans cette couche.
      if (active(ig)) then
        ztla(ig, 1) = zthl(ig, 1)
        zqta(ig, 1) = po(ig, 1)
        zqla(ig, 1) = zl(ig, 1)
        !cr: attention, prise en compte de f*(1)=1
        f_star(ig, 2) = alim_star(ig, 1)
        zw2(ig, 2) = 2. * RG * (ztv(ig, 1) - ztv(ig, 2)) / ztv(ig, 2)  &
                & * (zlev(ig, 2) - zlev(ig, 1))  &
                & * 0.4 * pphi(ig, 1) / (pphi(ig, 2) - pphi(ig, 1))
        w_est(ig, 2) = zw2(ig, 2)
      endif
    enddo

    !==============================================================================
    !boucle de calcul de la vitesse verticale dans le thermique
    !==============================================================================
    do l = 2, nlay - 1
      !==============================================================================


      ! On decide si le thermique est encore actif ou non
      ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
      do ig = 1, ngrid
        active(ig) = active(ig) &
                &                 .and. zw2(ig, l)>1.e-10 &
                &                 .and. f_star(ig, l) + alim_star(ig, l)>1.e-10
      enddo



      !---------------------------------------------------------------------------
      ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
      ! sans tenir compte du detrainement et de l'entrainement dans cette
      ! couche
      ! C'est a dire qu'on suppose
      ! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1)
      ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
      ! avant) a l'alimentation pour avoir un calcul plus propre
      !---------------------------------------------------------------------------

      ztemp(:) = zpspsk(:, l) * ztla(:, l - 1)
      CALL thermcell_qsat(ngrid, active, pplev(:, l), ztemp, zqta(:, l - 1), zqsat(:))

      do ig = 1, ngrid
        !       PRINT*,'active',active(ig),ig,l
        if(active(ig)) then
          zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig))
          ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l)
          zta_est(ig, l) = ztva_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))

          !------------------------------------------------
          !AJAM:nouveau calcul de w?
          !------------------------------------------------
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)

          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
          zdw2 = (afact) * zbuoy(ig, l) / (fact_epsilon)
          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)

          if (w_est(ig, l + 1)<0.) then
            w_est(ig, l + 1) = zw2(ig, l)
          endif
        endif
      enddo


      !-------------------------------------------------
      !calcul des taux d'entrainement et de detrainement
      !-------------------------------------------------

      do ig = 1, ngrid
        if (active(ig)) then

          zw2m = max(0.5 * (w_est(ig, l) + w_est(ig, l + 1)), 0.1)
          zw2m = w_est(ig, l + 1)
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          zbuoy(ig, l) = RG * (ztva_est(ig, l) - ztv(ig, l)) / ztv(ig, l)
          !          zbuoybis=zbuoy(ig,l)+RG*0.1/300.
          zbuoybis = zbuoy(ig, l)
          zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l)
          zdqt(ig, l) = max(zqta(ig, l - 1) - po(ig, l), 0.) / po(ig, l)

          entr_star(ig, l) = f_star(ig, l) * zdz * zbetalpha * MAX(0., &
                  afact * zbuoybis / zw2m - fact_epsilon)

          detr_star(ig, l) = f_star(ig, l) * zdz                        &
                  * MAX(1.e-3, -afact * zbetalpha * zbuoy(ig, l) / zw2m          &
                          + 0.012 * (zdqt(ig, l) / zw2m)**0.5)

          ! En dessous de lalim, on prend le max de alim_star et entr_star pour
          ! alim_star et 0 sinon
          if (l<lalim(ig)) then
            alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l))
            entr_star(ig, l) = 0.
          endif

          ! Calcul du flux montant normalise
          f_star(ig, l + 1) = f_star(ig, l) + alim_star(ig, l) + entr_star(ig, l)  &
                  - detr_star(ig, l)

        endif
      enddo


      !----------------------------------------------------------------------------
      !calcul de la vitesse verticale en melangeant Tl et qt du thermique
      !---------------------------------------------------------------------------
      activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10
      do ig = 1, ngrid
        if (activetmp(ig)) then
          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))

        endif
      enddo

      ztemp(:) = zpspsk(:, l) * ztla(:, l)
      CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l))

      do ig = 1, ngrid
        if (activetmp(ig)) then
          ! on ecrit de maniere conservative (sat ou non)
          !          T = Tl +Lv/Cp ql
          zqla(ig, l) = max(0., zqta(ig, l) - zqsatth(ig, l))
          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))
          zbuoy(ig, l) = RG * (ztva(ig, l) - ztv(ig, l)) / ztv(ig, l)
          zdz = zlev(ig, l + 1) - zlev(ig, l)
          zeps(ig, l) = (entr_star(ig, l) + alim_star(ig, l)) / (f_star(ig, l) * zdz)

          zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
          zdw2 = afact * zbuoy(ig, l) / (fact_epsilon)
          zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
        endif
      enddo

      if (prt_level>=20) PRINT*, 'coucou calcul detr 460: ig, l', ig, l

      !---------------------------------------------------------------------------
      !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
      !---------------------------------------------------------------------------

      nbpb = 0
      do ig = 1, ngrid
        if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then
          !               stop'On tombe sur le cas particulier de thermcell_dry'
          !               PRINT*,'On tombe sur le cas particulier de thermcell_plume'
          nbpb = nbpb + 1
          zw2(ig, l + 1) = 0.
          linter(ig) = l + 1
        endif

        if (zw2(ig, l + 1)<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.
        elseif (f_star(ig, l + 1)<0.) then
          linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l))  &
                  - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l))
          !           PRINT*,"linter plume", linter(ig)
          zw2(ig, l + 1) = 0.
        endif

        wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1))

        if (wa_moy(ig, l + 1)>wmaxa(ig)) then
          !   lmix est le niveau de la couche ou w (wa_moy) est maximum
          !on rajoute le calcul de lmix_bis
          if (zqla(ig, l)<1.e-10) then
            lmix_bis(ig) = l + 1
          endif
          lmix(ig) = l + 1
          wmaxa(ig) = wa_moy(ig, l + 1)
        endif
      enddo

      if (nbpb>0) then
        PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume'
      endif

      !=========================================================================
      ! FIN DE LA BOUCLE VERTICALE
    enddo
    !=========================================================================

    !on recalcule alim_star_tot
    do ig = 1, ngrid
      alim_star_tot(ig) = 0.
    enddo
    do ig = 1, ngrid
      do l = 1, lalim(ig) - 1
        alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig, l)
      enddo
    enddo

    if (prt_level>=20) PRINT*, 'coucou calcul detr 470: ig, l', ig, l
  end
END MODULE lmdz_thermcell_plume_6A
