MODULE lmdz_thermcell_plume

  ! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $

CONTAINS

  SUBROUTINE thermcell_plume(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)
    !--------------------------------------------------------------------------
    ! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam

    !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    !   This versions starts from a cleaning of thermcell_plume_6A (2019/01/20)
    !   thermcell_plume_6A is activate for flag_thermas_ed < 10
    !   thermcell_plume_5B for flag_thermas_ed < 20
    !   thermcell_plume for flag_thermals_ed>= 20
    !   Various options are controled by the flag_thermals_ed parameter
    !   = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8
    !   = 21 : the Jam strato-cumulus modif is not activated in detrainment
    !   = 29 : an other way to compute the modified buoyancy (to be tested)
    !--------------------------------------------------------------------------

    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(ngrid)
    integer, intent(out), dimension(ngrid) :: lmix_bis(ngrid)
    real, intent(out), dimension(ngrid) :: linter(ngrid)

    REAL, dimension(ngrid, nlay + 1) :: wa_moy
    REAL, dimension(ngrid, nlay) :: entr, detr
    REAL, dimension(ngrid, nlay) :: ztv_est
    REAL, dimension(ngrid, nlay) :: zqla_est
    REAL, dimension(ngrid, nlay) :: zta_est
    REAL, dimension(ngrid) :: ztemp, zqsat
    REAL zdw2, zdw2bis
    REAL zw2modif
    REAL zw2fact, zw2factbis
    REAL, dimension(ngrid, nlay) :: zeps

    REAL, dimension(ngrid) :: wmaxa

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

    real, dimension(ngrid, nlay) :: zbuoy, gamma, zdqt
    real zdz, zalpha, zw2m
    real, dimension(ngrid, nlay) :: zbuoyjam, zdqtjam
    real 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, 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.


    !-------------------------------------------------------------------------
    ! 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
          !=========================================================================

          !--------------------------------------------------
          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)

          !------------------------------------------------
          !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
          lm = Max(1, l - 2)
          w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2)
        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)
          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)

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

          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))

          if (iflag_thermals_ed == 20) then
            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))
          else
            entr_star(ig, l) = f_star(ig, l) * zdz * (&
                    mix0 * 0.1 / (zalpha + 0.001)               &
                            + zbetalpha * MAX(entr_min, &
                            afact * zbuoy(ig, l) / zw2m - fact_epsilon))
          endif

          ! 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
          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)
          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, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2)
        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
END MODULE lmdz_thermcell_plume
