MODULE lmdz_thermcell_main
! $Id: lmdz_thermcell_main.F90 4843 2024-03-04 17:58:03Z crio $
!
! A REGARDER !!!!!!!!!!!!!!!!!
! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)
! ATTENTION : dans thermcell_env, on condense potentiellement de l'eau. Mais comme on ne mélange pas l'eau liquide supposant qu'il n'y en n'a pas, c'est potentiellement un souci
CONTAINS

      subroutine thermcell_main(itap,ngrid,nlay,ptimestep  &
     &                  ,pplay,pplev,pphi,debut  &
     &                  ,puwind,pvwind,ptemp,p_o,ptemp_env, po_env  &
     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
     &                  ,ratqscth,ratqsdiff,zqsatth  &
     &                  ,zmax0, f0,zw2,fraca,ztv &
     &                  ,zpspsk,ztla,zthl,ztva &
     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,zcong &
#ifdef ISO         
     &      ,xtpo,xtpdoadj &
#endif         
     &   )



! USE necessaires pour les lignes importees de thermcell_env
      USE lmdz_thermcell_ini, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level
      USE lmdz_thermcell_ini, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals
      USE lmdz_thermcell_ini, ONLY: iflag_thermals_down,fact_thermals_down
      USE lmdz_thermcell_ini, ONLY: iflag_thermals_tenv
      USE lmdz_thermcell_ini, ONLY: RD,RG

      USE lmdz_thermcell_down, ONLY: thermcell_updown_dq
      USE lmdz_thermcell_closure, ONLY: thermcell_closure
      USE lmdz_thermcell_dq, ONLY: thermcell_dq
      USE lmdz_thermcell_dry, ONLY: thermcell_dry
      USE lmdz_thermcell_dv2, ONLY: thermcell_dv2
      USE lmdz_thermcell_env, ONLY: thermcell_env
      USE lmdz_thermcell_flux2, ONLY: thermcell_flux2
      USE lmdz_thermcell_height, ONLY: thermcell_height
      USE lmdz_thermcell_plume, ONLY: thermcell_plume
      USE lmdz_thermcell_plume_6A, ONLY: thermcell_plume_6A,thermcell_plume_5B

! USE necessaires pour les lignes importees de thermcell_env
   USE lmdz_thermcell_ini, ONLY : RLvCp,RKAPPA,RETV
   USE lmdz_thermcell_qsat, ONLY : thermcell_qsat


#ifdef ISO
  USE infotrac_phy, ONLY : ntiso
#ifdef ISOVERIF
  USE isotopes_mod, ONLY : iso_eau,iso_HDO
  USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
        iso_verif_aberrant_encadre
#endif
#endif


      IMPLICIT NONE

!=======================================================================
!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
!   Version du 09.02.07
!   Calcul du transport vertical dans la couche limite en presence
!   de "thermiques" explicitement representes avec processus nuageux
!
!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
!
!   le thermique est suppose homogene et dissipe par melange avec
!   son environnement. la longueur l_mix controle l'efficacite du
!   melange
!
!   Le calcul du transport des differentes especes se fait en prenant
!   en compte:
!     1. un flux de masse montant
!     2. un flux de masse descendant
!     3. un entrainement
!     4. un detrainement
!
! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
!    Introduction of an implicit computation of vertical advection in
!    the environment of thermal plumes in thermcell_dq
!    impl =     0 : explicit, 1 : implicit, -1 : old version
!    controled by iflag_thermals =
!       15, 16 run with impl=-1 : numerical convergence with NPv3
!       17, 18 run with impl=1  : more stable
!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
!
! Using
!    abort_physic 
!    iso_verif_aberrant_encadre
!    iso_verif_egalite
!    test_ltherm
!    thermcell_closure
!    thermcell_dq
!    thermcell_dry
!    thermcell_dv2
!    thermcell_env
!    thermcell_flux2
!    thermcell_height
!    thermcell_plume
!    thermcell_plume_5B
!    thermcell_plume_6A
!
!=======================================================================


!-----------------------------------------------------------------------
!   declarations:
!   -------------


!   arguments:
!   ----------
      integer, intent(in) :: itap,ngrid,nlay
      real, intent(in) ::  ptimestep
      real, intent(in), dimension(ngrid,nlay)    :: ptemp,puwind,pvwind,pplay,pphi,ptemp_env,po_env
! ATTENTION : zpspsk est inout et out mais c'est pas forcement pour de bonnes raisons (FH, 2023)
      real, intent(in), dimension(ngrid,nlay)    :: p_o
      real, intent(out), dimension(ngrid,nlay)    :: zpspsk
      real, intent(in), dimension(ngrid,nlay+1)  :: pplev
      integer, intent(out), dimension(ngrid) :: lmax
      real, intent(out), dimension(ngrid,nlay)   :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0
      real, intent(out), dimension(ngrid,nlay)   :: ztla,zqla,zqta,zqsatth,zthl
      real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca
      real, intent(inout), dimension(ngrid) :: zmax0,f0
      real, intent(out), dimension(ngrid,nlay) :: ztva,ztv
      logical, intent(in) :: debut
      real,intent(out), dimension(ngrid,nlay) :: ratqscth,ratqsdiff

      real, intent(out), dimension(ngrid) :: pcon
      real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3
      real, intent(out), dimension(ngrid) :: wmax_sec
      integer,intent(out), dimension(ngrid) :: lalim
      real, intent(out), dimension(ngrid,nlay+1) :: fm
      real, intent(out), dimension(ngrid,nlay) :: alim_star
      real, intent(out), dimension(ngrid) :: zmax,zcong

!   local:
!   ------


      integer,save :: igout=1
!$OMP THREADPRIVATE(igout)
      integer,save :: lunout1=6
!$OMP THREADPRIVATE(lunout1)
      integer,save :: lev_out=10
!$OMP THREADPRIVATE(lev_out)

      real lambda, zf,zf2,var,vardiff,CHI
      integer ig,k,l,ierr,ll
      logical sorties
      real, dimension(ngrid) :: linter,zmix, zmax_sec,lintercong
      integer,dimension(ngrid) :: lmin,lmix,lmix_bis,nivcon, lcong
      real, dimension(ngrid,nlay) :: ztva_est
      real, dimension(ngrid,nlay) :: deltaz,zlay,zdthladj,zu,zv,z_o,zl,zva,zua,z_oa
      real, dimension(ngrid,nlay) :: ztemp_env ! temperarure liquide de l'environnement
      real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2
      real, dimension(ngrid,nlay) :: rho,masse
      real, dimension(ngrid,nlay+1) :: zw_est,zlev
      real, dimension(ngrid) :: wmax,wmax_tmp
      real, dimension(ngrid,nlay+1) :: f_star
      real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos
      real, dimension(ngrid,nlay) :: zqsat,csc
      real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f
      real, dimension(ngrid,nlay) :: entrdn,detrdn
      logical, dimension(ngrid,nlay) :: mask

      character (len=20) :: modname='thermcell_main'
      character (len=80) :: abort_message


#ifdef ISO
      REAL xtpo(ntiso,ngrid,nlay),xtpdoadj(ntiso,ngrid,nlay)
      REAL xtzo(ntiso,ngrid,nlay)
      REAL xtpdoadj_tmp(ngrid,nlay)
      REAL xtpo_tmp(ngrid,nlay)
      REAL xtzo_tmp(ngrid,nlay)
      integer ixt
#endif

!

!-----------------------------------------------------------------------
!   initialisation:
!   ---------------
!
   fm=0. ; entr=0. ; detr=0.

      if (prt_level.ge.1) print*,'thermcell_main V4'

       sorties=.true.
      IF(ngrid.NE.ngrid) THEN
         PRINT*
         PRINT*,'STOP dans convadj'
         PRINT*,'ngrid    =',ngrid
         PRINT*,'ngrid  =',ngrid
      ENDIF
!
!print*,'thermcell_main debut'
!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
     do ig=1,ngrid
         f0(ig)=max(f0(ig),1.e-2)
         zmax0(ig)=max(zmax0(ig),40.)
!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
     enddo

      if (prt_level.ge.20) then
       do ig=1,ngrid
          print*,'th_main ig f0',ig,f0(ig)
       enddo
      endif

!-----------------------------------------------------------------------
! Calcul de T,q,ql a partir de Tl et qT dans l environnement
!   --------------------------------------------------------------------
!
          ! On condense l'eau liquide si besoin.
          ! En fait on arrive ici d'habitude (jusque 6A) après réévaporation
          ! Dans une nouvelle mouture, on passe les profiles 
          ! avant la couche limite : iflag_thermals_tenv=1
          !     dés le début de la physique : iflag_thermals_tenv=2
          ! Mais même pour 2) on ne veut sans doute pas réévaporer.
          ! On veut comparer thetav dans le thermique, après condensation,
          ! avec le theta_v effectif de l'environnement.

      if (iflag_thermals_tenv - 10 * ( iflag_thermals_tenv / 10 ) == 0) then

          CALL thermcell_env(ngrid,nlay,p_o,ptemp_env,puwind,pvwind,pplay,  &
         &           pplev,z_o,ztemp_env,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lcong,lintercong,lev_out)

      else

        ! Chantier en cours : ne pas effacer (Fredho). 15 septembre 2023
        ! Dans la version originale de thermcell_env, on condense l'eau de l'environnement
        ! pour calculer une temperature potentielle liquide.
        ! On en déduit un Theta v.

 ! ... 
        ! contenu de thermcell_env
        !    SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
        ! &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
        ! contenu thermcell_env : call thermcell_qsat(ngrid*nlay,mask,pplev,pt,po,pqsat)
        ! contenu thermcell_env : do ll=1,nlay
        ! contenu thermcell_env :    do ig=1,ngrid
        ! contenu thermcell_env :       zl(ig,ll) = max(0.,po(ig,ll)-pqsat(ig,ll))
        ! contenu thermcell_env :       zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll)         !   T = Tl + Lv/Cp ql
        ! contenu thermcell_env :       zo(ig,ll) = po(ig,ll)-zl(ig,ll)
        ! contenu thermcell_env :    enddo
        ! contenu thermcell_env : enddo
        ! contenu thermcell_env : do ll=1,nlay
        ! contenu thermcell_env :    do ig=1,ngrid
        ! contenu thermcell_env :        zpspsk(ig,ll)=(pplay(ig,ll)/100000.)**RKAPPA
        ! contenu thermcell_env :        zu(ig,ll)=pu(ig,ll)
        ! contenu thermcell_env :        zv(ig,ll)=pv(ig,ll)
        ! contenu thermcell_env :          ztv(ig,ll)=zh(ig,ll)/zpspsk(ig,ll)
        ! contenu thermcell_env :          ztv(ig,ll)=ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll))
        ! contenu thermcell_env :          zthl(ig,ll)=pt(ig,ll)/zpspsk(ig,ll)
        ! contenu thermcell_env :    enddo
        ! contenu thermcell_env : enddo

        do l=1,nlay
            do ig=1,ngrid
                 zl(ig,l)=0.
                 zu(ig,l)=puwind(ig,l)
                 zv(ig,l)=pvwind(ig,l)
                 ztemp_env(ig,l)=ptemp_env(ig,l)
                 zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA
                 ztv(ig,l)=ztemp_env(ig,l)/zpspsk(ig,l)
                 ztv(ig,l)=ztv(ig,l)*(1.+RETV*po_env(ig,l))
                 zthl(ig,l)=ptemp(ig,l)/zpspsk(ig,l)
                 mask(ig,l)=.true.
            enddo
        enddo
        call thermcell_qsat(ngrid*nlay,mask,pplev,ptemp_env,p_o,zqsat)
          
      endif
       
      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'

!------------------------------------------------------------------------
!                       --------------------
!
!
!                       + + + + + + + + + + +
!
!
!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
!  wh,wt,wo ...
!
!                       + + + + + + + + + + +  zh,zu,zv,z_o,rho
!
!
!                       --------------------   zlev(1)
!                       \\\\\\\\\\\\\\\\\\\\
!
!

!-----------------------------------------------------------------------
!   Calcul des altitudes des couches
!-----------------------------------------------------------------------

      do l=2,nlay
         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
      enddo
      zlev(:,1)=0.
      zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
      do l=1,nlay
         zlay(:,l)=pphi(:,l)/RG
      enddo
      do l=1,nlay
         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
      enddo

!-----------------------------------------------------------------------
!   Calcul des densites et masses
!-----------------------------------------------------------------------

      rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
      if (prt_level.ge.10) write(lunout,*) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
      rhobarz(:,1)=rho(:,1)
      do l=2,nlay
         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
      enddo
      do l=1,nlay
         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
      enddo
      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'

!------------------------------------------------------------------
!
!             /|\
!    --------  |  F_k+1 -------   
!                              ----> D_k
!             /|\              <---- E_k , A_k
!    --------  |  F_k --------- 
!                              ----> D_k-1
!                              <---- E_k-1 , A_k-1
!
!
!
!
!
!    ---------------------------
!
!    ----- F_lmax+1=0 ----------         \
!            lmax     (zmax)              |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------          |
!                                         |  E
!    ---------------------------          |  D
!                                         |
!    ---------------------------          |
!                                         |
!    ---------------------------  \       |
!            lalim                 |      |
!    ---------------------------   |      |
!                                  |      |
!    ---------------------------   |      |
!                                  | A    |
!    ---------------------------   |      |
!                                  |      |
!    ---------------------------   |      |
!    lmin  (=1 pour le moment)     |      |
!    ----- F_lmin=0 ------------  /      /
!
!    ---------------------------
!    //////////////////////////
!
!
!=============================================================================
!  Calculs initiaux ne faisant pas intervenir les changements de phase
!=============================================================================

!------------------------------------------------------------------
!  1. alim_star est le profil vertical de l'alimentation a la base du
!     panache thermique, calcule a partir de la flotabilite de l'air sec
!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
!------------------------------------------------------------------
!
      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
      lmin=1

!-----------------------------------------------------------------------------
!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
!     panache sec conservatif (e=d=0) alimente selon alim_star 
!     Il s'agit d'un calcul de type CAPE
!     zmax_sec est utilise pour determiner la geometrie du thermique.
!------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
!calcul du melange et des variables dans le thermique
!--------------------------------------------------------------------------------
!
      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out

!=====================================================================
! Old version of thermcell_plume in thermcell_plume_6A.F90
! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
! to the 5B and 6A versions used for CMIP5 and CMIP6.
! The latest was previously named thermcellV1_plume.
! The new thermcell_plume is a clean version (removing obsolete
! options) of thermcell_plume_6A.
! The 3 versions are controled by
! flag_thermals_ed <= 9 thermcell_plume_6A
!                  <= 19 thermcell_plume_5B
!                  else thermcell_plume (default 20 for convergence with 6A)
! Fredho
!=====================================================================

      if (iflag_thermals_ed<=9) then
!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'
         CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,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,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     &    ,lev_out,lunout1,igout)

      elseif (iflag_thermals_ed<=19) then
!        print*,'THERM RIO et al 2010, version d Arnaud'
         CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,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,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     &    ,lev_out,lunout1,igout)
      else
         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,p_o,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,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     &    ,lev_out,lunout1,igout)
      endif

      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out

      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')

      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
      if (prt_level.ge.10) then
         write(lunout1,*) 'Dans thermcell_main 2'
         write(lunout1,*) 'lmin ',lmin(igout)
         write(lunout1,*) 'lalim ',lalim(igout)
         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
      endif

!-------------------------------------------------------------------------------
! Calcul des caracteristiques du thermique:zmax,zmix,wmax
!-------------------------------------------------------------------------------
!
      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lcong,lintercong,lmix,zw2,  &
     &           zlev,lmax,zmax,zmax0,zmix,wmax,zcong)
! Attention, w2 est transforme en sa racine carree dans cette routine
! Le probleme vient du fait que linter et lmix sont souvent egaux a 1.
      wmax_tmp=0.
      do  l=1,nlay
         wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
      enddo
!     print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax



      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
      call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
      call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')

      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'

!-------------------------------------------------------------------------------
! Fermeture,determination de f
!-------------------------------------------------------------------------------
!
!
      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
    &                      lalim,lmin,zmax_sec,wmax_sec)

 
call test_ltherm(ngrid,nlay,pplay,lmin,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')

      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
      if (prt_level.ge.10) then
         write(lunout1,*) 'Dans thermcell_main 1b'
         write(lunout1,*) 'lmin ',lmin(igout)
         write(lunout1,*) 'lalim ',lalim(igout)
         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
     &    ,l=1,lalim(igout)+4)
      endif




! Choix de la fonction d'alimentation utilisee pour la fermeture.
! Apparemment sans importance
      alim_star_clos(:,:)=alim_star(:,:)
      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
!
!CR Appel de la fermeture seche 
      if (iflag_thermals_closure.eq.1) then

     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
    &   zlev,lalim,alim_star_clos,zmax_sec,wmax_sec,f)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Appel avec les zmax et wmax tenant compte de la condensation
! Semble moins bien marcher
     else if (iflag_thermals_closure.eq.2) then

     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
    &   zlev,lalim,alim_star,zmax,wmax,f)


     endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'

      if (tau_thermals>1.) then
         lambda=exp(-ptimestep/tau_thermals)
         f0=(1.-lambda)*f+lambda*f0
      else
         f0=f
      endif

! Test valable seulement en 1D mais pas genant
      if (.not. (f0(1).ge.0.) ) then
              abort_message = '.not. (f0(1).ge.0.)'
              CALL abort_physic (modname,abort_message,1)
      endif

!-------------------------------------------------------------------------------
!deduction des flux

      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
     &       lalim,lmax,alim_star,  &
     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
     &       detr,zqla,lev_out,lunout1,igout)

!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)

      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
      call test_ltherm(ngrid,nlay,pplay,lalim,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
      call test_ltherm(ngrid,nlay,pplay,lmax ,ztv,p_o,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')

!------------------------------------------------------------------
!   On ne prend pas directement les profils issus des calculs precedents
!   mais on s'autorise genereusement une relaxation vers ceci avec
!   une constante de temps tau_thermals (typiquement 1800s).
!------------------------------------------------------------------

      if (tau_thermals>1.) then
         lambda=exp(-ptimestep/tau_thermals)
         fm0=(1.-lambda)*fm+lambda*fm0
         entr0=(1.-lambda)*entr+lambda*entr0
         detr0=(1.-lambda)*detr+lambda*detr0
      else
         fm0=fm
         entr0=entr
         detr0=detr
      endif

!------------------------------------------------------------------
! Calcul de la fraction de l'ascendance
!------------------------------------------------------------------
      do ig=1,ngrid
         fraca(ig,1)=0.
         fraca(ig,nlay+1)=0.
      enddo
      do l=2,nlay
         do ig=1,ngrid
            if (zw2(ig,l).gt.1.e-10) then
            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
            else
            fraca(ig,l)=0.
            endif
         enddo
      enddo
     
!c------------------------------------------------------------------
!   calcul du transport vertical
!------------------------------------------------------------------
      IF (iflag_thermals_down .GT. 0) THEN
        if (debut) print*,'WARNING !!! routine thermcell_down en cours de developpement'
        entrdn=fact_thermals_down*detr0
        detrdn=fact_thermals_down*entr0 
        ! we want to transport potential temperature, total water and momentum
        CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zthl,zdthladj)
        CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,p_o,pdoadj)
        CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zu,pduadj)
        CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,entrdn,detrdn,masse,zv,pdvadj)
      ELSE
      !--------------------------------------------------------------

        ! Temperature potentielle liquide effectivement mélangée par les thermiques
        do ll=1,nlay
           do ig=1,ngrid
              zthl(ig,ll)=ptemp(ig,ll)/zpspsk(ig,ll)
           enddo
        enddo
        call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
        &                    zthl,zdthladj,zta,lev_out)

        do ll=1,nlay
           do ig=1,ngrid
              z_o(ig,ll)=p_o(ig,ll)
           enddo
        enddo
        call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
        &                   z_o,pdoadj,z_oa,lev_out)

#ifdef ISO
        ! C Risi: on utilise directement la meme routine
        do ixt=1,ntiso
          do ll=1,nlay
            DO ig=1,ngrid
                xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll)
                xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll)
            enddo
          enddo
          call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
     &                   xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out)
          do ll=1,nlay
            DO ig=1,ngrid
                xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll)
            enddo
          enddo
        enddo
#endif

#ifdef ISO      
#ifdef ISOVERIF
      DO  ll=1,nlay
        DO ig=1,ngrid
          if (iso_eau.gt.0) then
              call iso_verif_egalite(xtpo(iso_eau,ig,ll), &
     &          p_o(ig,ll),'thermcell_main 594')
              call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), &
     &          pdoadj(ig,ll),'thermcell_main 596')
          endif
          if (iso_HDO.gt.0) then
              call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) &
     &           /p_o(ig,ll),'thermcell_main 610')
          endif
        enddo
      enddo !DO  ll=1,nlay 
      write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq'
#endif      
#endif


!------------------------------------------------------------------
!  calcul du transport vertical du moment horizontal
!------------------------------------------------------------------

!IM 090508  
      if (dvdq == 0 ) then

! Calcul du transport de V tenant compte d'echange par gradient
! de pression horizontal avec l'environnement

         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
!    &    ,fraca*dvdq,zmax &
     &    ,fraca,zmax &
     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)

      else

! calcul purement conservatif pour le transport de V
         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
     &    ,zu,pduadj,zua,lev_out)
         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
     &    ,zv,pdvadj,zva,lev_out)

      endif
    ENDIF

!     print*,'13 OK convect8'
      do l=1,nlay
         do ig=1,ngrid
           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)  
         enddo
      enddo

      if (prt_level.ge.1) print*,'14 OK convect8'
!------------------------------------------------------------------
!   Calculs de diagnostiques pour les sorties
!------------------------------------------------------------------
!calcul de fraca pour les sorties
      
      if (sorties) then
      if (prt_level.ge.1) print*,'14a OK convect8'
! calcul du niveau de condensation
! initialisation
      do ig=1,ngrid
         nivcon(ig)=0
         zcon(ig)=0.
      enddo 
!nouveau calcul
      do ig=1,ngrid
      ! WARNING !!! verifier que c'est bien ztemp_env qu'on veut là
      CHI=ztemp_env(ig,1)/(1669.0-122.0*z_o(ig,1)/zqsat(ig,1)-ztemp_env(ig,1))
      pcon(ig)=pplay(ig,1)*(z_o(ig,1)/zqsat(ig,1))**CHI
      enddo
!IM   do k=1,nlay
      do k=1,nlay-1
         do ig=1,ngrid
         if ((pcon(ig).le.pplay(ig,k))  &
     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
         endif
         enddo
      enddo
!IM
      ierr=0
      do ig=1,ngrid
        if (pcon(ig).le.pplay(ig,nlay)) then 
           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
           ierr=1
        endif
      enddo
!      if (ierr==1) then
!           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
!           CALL abort_physic (modname,abort_message,1)
!      endif

      if (prt_level.ge.1) print*,'14b OK convect8'
      do k=nlay,1,-1
         do ig=1,ngrid
            if (zqla(ig,k).gt.1e-10) then
               nivcon(ig)=k
               zcon(ig)=zlev(ig,k)
            endif
         enddo
      enddo
      if (prt_level.ge.1) print*,'14c OK convect8'
!calcul des moments
!initialisation
      do l=1,nlay
         do ig=1,ngrid
            q2(ig,l)=0.
            wth2(ig,l)=0.
            wth3(ig,l)=0.
            ratqscth(ig,l)=0.
            ratqsdiff(ig,l)=0.
         enddo
      enddo      
      if (prt_level.ge.1) print*,'14d OK convect8'
      if (prt_level.ge.10)write(lunout,*)                                &
    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
      do l=1,nlay
         do ig=1,ngrid
            zf=fraca(ig,l)
            zf2=zf/(1.-zf)
!
            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
            if(zw2(ig,l).gt.1.e-10) then
             wth2(ig,l)=zf2*(zw2(ig,l))**2
            else
             wth2(ig,l)=0.
            endif
            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
            q2(ig,l)=zf2*(zqta(ig,l)*1000.-p_o(ig,l)*1000.)**2
!test: on calcul q2/p_o=ratqsc
            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(p_o(ig,l)*1000.))
         enddo
      enddo
!calcul des flux: q, thetal et thetav
      do l=1,nlay
         do ig=1,ngrid
      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-p_o(ig,l)*1000.)
      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
         enddo
      enddo

!calcul du ratqscdiff
      if (prt_level.ge.1) print*,'14e OK convect8'
      var=0.
      vardiff=0.
      ratqsdiff(:,:)=0.

      do l=1,nlay
         do ig=1,ngrid
            if (l<=lalim(ig)) then
            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
            endif
         enddo
      enddo

      if (prt_level.ge.1) print*,'14f OK convect8'

      do l=1,nlay
         do ig=1,ngrid
            if (l<=lalim(ig)) then
               zf=fraca(ig,l)
               zf2=zf/(1.-zf)
               vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
            endif
         enddo
      enddo

      if (prt_level.ge.1) print*,'14g OK convect8'
         do l=1,nlay
            do ig=1,ngrid
               ratqsdiff(ig,l)=sqrt(vardiff)/(p_o(ig,l)*1000.)   
            enddo
         enddo 
      endif

      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'

!print*,'thermcell_main fin'
 RETURN
      end subroutine thermcell_main

!=============================================================================
!/////////////////////////////////////////////////////////////////////////////
!=============================================================================
      subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,p_o,ztva, &  ! in
    &            zqla,f_star,zw2,comment)                          ! in
!=============================================================================
      USE lmdz_thermcell_ini, ONLY: prt_level
      IMPLICIT NONE

      integer i, k, ngrid,nlay
      real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,p_o,ztva,zqla
      real, intent(in), dimension(ngrid,nlay) :: f_star,zw2
      integer, intent(in), dimension(ngrid) :: long
      real seuil
      character*21 comment

      seuil=0.25

      if (prt_level.ge.1) THEN
       print*,'WARNING !!! TEST ',comment
      endif
      return

!  test sur la hauteur des thermiques ...
         do i=1,ngrid
!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
           if (prt_level.ge.10) then
               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
               do k=1,nlay
                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*p_o(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
               enddo
           endif
         enddo


      return
      end

! nrlmd le 10/04/2012   Transport de la TKE par le thermique moyen pour la fermeture en ALP 
!                       On transporte pbl_tke pour donner therm_tke
!                       Copie conforme de la subroutine DTKE dans physiq.F ecrite par Frederic Hourdin

!=======================================================================
!///////////////////////////////////////////////////////////////////////
!=======================================================================

      subroutine thermcell_tke_transport( &
     &     ngrid,nlay,ptimestep,fm0,entr0,rg,pplev,  &   ! in
     &     therm_tke_max)                                ! out
      USE lmdz_thermcell_ini, ONLY: prt_level
      implicit none

!=======================================================================
!
!   Calcul du transport verticale dans la couche limite en presence
!   de "thermiques" explicitement representes
!   calcul du dq/dt une fois qu'on connait les ascendances
!
!=======================================================================

      integer ngrid,nlay

      real, intent(in) :: ptimestep
      real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev
      real, intent(in), dimension(ngrid,nlay) :: entr0
      real, intent(in) :: rg
      real, intent(out), dimension(ngrid,nlay) :: therm_tke_max

      real detr0(ngrid,nlay)
      real masse0(ngrid,nlay)
      real masse(ngrid,nlay),fm(ngrid,nlay+1)
      real entr(ngrid,nlay)
      real q(ngrid,nlay)
      integer lev_out                           ! niveau pour les print

      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
      integer ig,k


      lev_out=0


      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'

!   calcul du detrainement
      do k=1,nlay
         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
      enddo


! Decalage vertical des entrainements et detrainements.
      masse(:,1)=0.5*masse0(:,1)
      entr(:,1)=0.5*entr0(:,1)
      detr(:,1)=0.5*detr0(:,1)
      fm(:,1)=0.
      do k=1,nlay-1
         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
      enddo
      fm(:,nlay+1)=0.


   q(:,:)=therm_tke_max(:,:)
!!! nrlmd le 16/09/2010
      do ig=1,ngrid
         qa(ig,1)=q(ig,1)
      enddo
!!!

    if (1==1) then
      do k=2,nlay
         do ig=1,ngrid
            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     &         1.e-5*masse(ig,k)) then
         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     &         /(fm(ig,k+1)+detr(ig,k))
            else
               qa(ig,k)=q(ig,k)
            endif
            if (qa(ig,k).lt.0.) then
!               print*,'qa<0!!!'
            endif
            if (q(ig,k).lt.0.) then
!               print*,'q<0!!!'
            endif
         enddo
      enddo

! Calcul du flux subsident

      do k=2,nlay
         do ig=1,ngrid
            wqd(ig,k)=fm(ig,k)*q(ig,k)
            if (wqd(ig,k).lt.0.) then
!               print*,'wqd<0!!!'
            endif
         enddo
      enddo
      do ig=1,ngrid
         wqd(ig,1)=0.
         wqd(ig,nlay+1)=0.
      enddo

! Calcul des tendances
      do k=1,nlay
         do ig=1,ngrid
            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
     &               -wqd(ig,k)+wqd(ig,k+1))  &
     &               *ptimestep/masse(ig,k)
         enddo
      enddo

 endif

   therm_tke_max(:,:)=q(:,:)

      return
!!! fin nrlmd le 10/04/2012
     end

END MODULE lmdz_thermcell_main
