Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2104)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2127)
@@ -2,25 +2,14 @@
 !
 !
-SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep,                          &
-                          pplay,pplev,pphi,firstcall,                         &
-                          pu,pv,pt,po,                                        &
-                          pduadj,pdvadj,pdtadj,pdoadj,                        &
-                          f0,fm0,entr0,detr0,                                 &
-                          zqta,zqla,ztv,ztva,ztla,zthl,zqsatth,               &
-                          zw2,fraca,                                          &
-                          lmin,lmix,lalim,lmax,                               &
-                          zpspsk,ratqscth,ratqsdiff,                          &
-                          Ale_bl,Alp_bl,lalim_conv,wght_th,                   &
-!!! nrlmd le 10/04/2012
-                          pbl_tke,pctsrf,omega,airephy,                       &
-                          zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,  &
-                          n2,s2,ale_bl_stat,                                  &
-                          therm_tke_max,env_tke_max,                          &
-                          alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,         &
-                          alp_bl_conv,alp_bl_stat)
-!!! fin nrlmd le 10/04/2012
-      
-      
-!==============================================================================
+SUBROUTINE thermcell_main(ngrid,nlay,nq,ptimestep,firstcall,                  &
+                          pplay,pplev,pphi,zpopsk,                            &
+                          pu,pv,pt,pq,                                        &
+                          pduadj,pdvadj,pdtadj,pdqadj,                        &
+                          f0,fm0,entr0,detr0,zw2,fraca,                       &
+                          zqta,zqla,ztv,ztva,zhla,zhl,zqsa,                   &
+                          lmin,lmix,lalim,lmax)
+      
+      
+!===============================================================================
 !   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
 !   Version du 09.02.07
@@ -44,13 +33,21 @@
 !    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
+!    impl = 0 : explicit ; impl = 1 : implicit ; impl =-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"
-!
-!==============================================================================
+! 
+! Major changes 2018-19 (AB alexandre.boissinot@lmd.jussieu.fr)
+!    New detr and entre formulae (no longer alimentation)
+!    lmin can be greater than 1
+!    Mix every tracer (EN COURS)
+!    Old version of thermcell_dq is removed
+!    Alternative version thermcell_dv2 is removed
+! 
+!===============================================================================
       
       USE thermcell_mod
+      USE tracer_h, ONLY: igcm_h2o_vap
       USE print_control_mod, ONLY: lunout, prt_level
       
@@ -58,230 +55,109 @@
       
       
-!==============================================================================
+!===============================================================================
 ! Declaration
-!==============================================================================
-      
-!   inputs:
-!   -------
-      
-      INTEGER itap
-      INTEGER ngrid
-      INTEGER nlay
+!===============================================================================
+      
+!     Inputs:
+!     -------
+      
+      INTEGER ngrid, nlay, nq
       
       REAL ptimestep
-      REAL pplay(ngrid,nlay)
-      REAL pplev(ngrid,nlay+1)
-      REAL pphi(ngrid,nlay)
-      
-      REAL pt(ngrid,nlay)                       ! 
-      REAL pu(ngrid,nlay)                       ! 
-      REAL pv(ngrid,nlay)                       ! 
-      REAL po(ngrid,nlay)                       ! 
+      REAL pplay(ngrid,nlay)                    ! Layer pressure
+      REAL pplev(ngrid,nlay+1)                  ! Level pressure
+      REAL pphi(ngrid,nlay)                     ! Geopotential
+      
+      REAL pu(ngrid,nlay)                       ! Zonal wind
+      REAL pv(ngrid,nlay)                       ! Meridional wind
+      REAL pt(ngrid,nlay)                       ! Temperature
+      REAL pq(ngrid,nlay,nq)                    ! Tracers mass mixing ratio
       
       LOGICAL firstcall
       
-!   outputs:
-!   --------
-      
-      REAL pdtadj(ngrid,nlay)                   ! t convective variations
+!     Outputs:
+!     --------
+      
       REAL pduadj(ngrid,nlay)                   ! u convective variations
       REAL pdvadj(ngrid,nlay)                   ! v convective variations
-      REAL pdoadj(ngrid,nlay)                   ! water convective variations
-      
+      REAL pdtadj(ngrid,nlay)                   ! t convective variations
+      REAL pdqadj(ngrid,nlay,nq)                ! q convective variations
+      
+      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
       REAL fm0(ngrid,nlay+1)                    ! mass flux      (after possible time relaxation)
       REAL entr0(ngrid,nlay)                    ! entrainment    (after possible time relaxation)
       REAL detr0(ngrid,nlay)                    ! detrainment    (after possible time relaxation)
-      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
-      
-!   local:
-!   ------
-      
-      INTEGER, save :: dvimpl=1
-!$OMP THREADPRIVATE(dvimpl)
-      
-      INTEGER, save :: dqimpl=-1
-!$OMP THREADPRIVATE(dqimpl)
-      
-      INTEGER, SAVE :: igout=1
-!$OMP THREADPRIVATE(igout)
-      
-      INTEGER, SAVE :: lunout1=6
-!$OMP THREADPRIVATE(lunout1)
-      
-      INTEGER, SAVE :: lev_out=10
-!$OMP THREADPRIVATE(lev_out)
-      
-      INTEGER ig,k,l,ll,ierr
-      INTEGER lmix_bis(ngrid)
-      INTEGER lmax(ngrid)                       ! 
-      INTEGER lmin(ngrid)                       ! 
-      INTEGER lalim(ngrid)                      ! 
-      INTEGER lmix(ngrid)                       ! 
-      
-      REAL linter(ngrid)
-      REAL zmix(ngrid)
-      REAL zmax(ngrid)
-      REAL zw2(ngrid,nlay+1)
-      REAL zw_est(ngrid,nlay+1)
-      REAL zmax_sec(ngrid)
-      
-      REAL zlay(ngrid,nlay)                     ! layers altitude
-      REAL zlev(ngrid,nlay+1)                   ! levels altitude
-      REAL rho(ngrid,nlay)                      ! layers density
-      REAL rhobarz(ngrid,nlay)                  ! levels density
-      REAL deltaz(ngrid,nlay)                   ! layers heigth
-      REAL masse(ngrid,nlay)                    ! layers mass
-      REAL zpspsk(ngrid,nlay)                   ! Exner function
-      
-      REAL zu(ngrid,nlay)                       ! environment zonal wind
-      REAL zv(ngrid,nlay)                       ! environment meridional wind
-      REAL zo(ngrid,nlay)                       ! environment water tracer
-      REAL zva(ngrid,nlay)                      ! plume zonal wind
-      REAL zua(ngrid,nlay)                      ! plume meridional wind
-      REAL zoa(ngrid,nlay)                      ! plume water tracer
-      
-      REAL zt(ngrid,nlay)                       ! T    environment
-      REAL zh(ngrid,nlay)                       ! T,TP environment
-      REAL zthl(ngrid,nlay)                     ! TP   environment
-      REAL ztv(ngrid,nlay)                      ! TPV  environment ?
-      REAL zl(ngrid,nlay)                       ! ql   environment
-      
-      REAL zta(ngrid,nlay)                      ! 
-      REAL zha(ngrid,nlay)                      ! TRPV plume
-      REAL ztla(ngrid,nlay)                     ! TP   plume
-      REAL ztva(ngrid,nlay)                     ! TRPV plume
-      REAL ztva_est(ngrid,nlay)                 ! TRPV plume (temporary)
+      
+!     Local:
+!     ------
+      
+      INTEGER ig, k, l, iq
+      INTEGER lmax(ngrid)                       ! Highest layer reached by the plume
+      INTEGER lmix(ngrid)                       ! Layer in which plume vertical speed is maximal
+      INTEGER lmin(ngrid)                       ! First unstable layer
+      
+      REAL zlay(ngrid,nlay)                     ! Layers altitudes
+      REAL zlev(ngrid,nlay+1)                   ! Levels altitudes
+      REAL rho(ngrid,nlay)                      ! Layers densities
+      REAL rhobarz(ngrid,nlay)                  ! Levels densities
+      REAL masse(ngrid,nlay)                    ! Layers masses
+      REAL zpopsk(ngrid,nlay)                   ! Exner function
+      
+      REAL zu(ngrid,nlay)                       ! u    environment
+      REAL zv(ngrid,nlay)                       ! v    environment
+      REAL zt(ngrid,nlay)                       ! TR   environment
+      REAL zqt(ngrid,nlay)                      ! qt   environment
+      REAL zql(ngrid,nlay)                      ! ql   environment
+      REAL zhl(ngrid,nlay)                      ! TP   environment
+      REAL ztv(ngrid,nlay)                      ! TRPV environment
+      REAL zqs(ngrid,nlay)                      ! qsat environment
+      
+      REAL zua(ngrid,nlay)                      ! u    plume
+      REAL zva(ngrid,nlay)                      ! v    plume
+      REAL zta(ngrid,nlay)                      ! TR   plume
       REAL zqla(ngrid,nlay)                     ! qv   plume
       REAL zqta(ngrid,nlay)                     ! qt   plume
-      
-      REAL wmax(ngrid)                          ! maximal vertical speed
-      REAL wmax_tmp(ngrid)                      ! temporary maximal vertical speed
-      REAL wmax_sec(ngrid)                      ! maximal vertical speed if dry case
-      
-      REAL fraca(ngrid,nlay+1)                  ! updraft fraction
-      REAL f_star(ngrid,nlay+1)                 ! normalized mass flux
-      REAL entr_star(ngrid,nlay)                ! normalized entrainment
-      REAL detr_star(ngrid,nlay)                ! normalized detrainment
-      REAL alim_star_tot(ngrid)                 ! integrated alimentation
-      REAL alim_star(ngrid,nlay)                ! normalized alimentation
-      REAL alim_star_clos(ngrid,nlay)           ! closure alimentation
-      
-      REAL fm(ngrid,nlay+1)                     ! mass flux
-      REAL entr(ngrid,nlay)                     ! entrainment
-      REAL detr(ngrid,nlay)                     ! detrainment
-      REAL f(ngrid)                             ! mass flux norm
-      
-      REAL zdthladj(ngrid,nlay)                 ! 
-      REAL lambda                               ! time relaxation coefficent
-      
-      REAL zsortie(ngrid,nlay)
-      REAL zsortie1d(ngrid)
-      REAL susqr2pi, Reuler
-      REAL zf
-      REAL zf2
-      REAL thetath2(ngrid,nlay)
-      REAL wth2(ngrid,nlay)
-      REAL wth3(ngrid,nlay)
-      REAL q2(ngrid,nlay)
-! FH : probleme de dimensionnement avec l'allocation dynamique
-!     common/comtherm/thetath2,wth2
-      real wq(ngrid,nlay)
-      real wthl(ngrid,nlay)
-      real wthv(ngrid,nlay)
-      real ratqscth(ngrid,nlay)
-      real var
-      real vardiff
-      real ratqsdiff(ngrid,nlay)
-! niveau de condensation
-      integer nivcon(ngrid)
-      real zcon(ngrid)
-      real CHI
-      real zcon2(ngrid)
-      real pcon(ngrid)
-      real zqsat(ngrid,nlay)
-      real zqsatth(ngrid,nlay)
-      real zlevinter(ngrid)
-      real seuil
-      
-!!! nrlmd le 10/04/2012
-      
-!------Entrees
-      real pbl_tke(ngrid,nlay+1,nbsrf)
-      real pctsrf(ngrid,nbsrf)
-      real omega(ngrid,nlay)
-      real airephy(ngrid)
-!------Sorties
-      real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
-      real therm_tke_max0(ngrid)
-      real env_tke_max0(ngrid)
-      real n2(ngrid),s2(ngrid)
-      real ale_bl_stat(ngrid)
-      real therm_tke_max(ngrid,nlay)
-      real env_tke_max(ngrid,nlay)
-      real alp_bl_det(ngrid)
-      real alp_bl_fluct_m(ngrid)
-      real alp_bl_fluct_tke(ngrid)
-      real alp_bl_conv(ngrid)
-      real alp_bl_stat(ngrid)
-!------Local
-      integer nsrf
-      real rhobarz0(ngrid)                      ! Densite au LCL
-      logical ok_lcl(ngrid)                     ! Existence du LCL des thermiques
-      integer klcl(ngrid)                       ! Niveau du LCL
-      real interp(ngrid)                        ! Coef d'interpolation pour le LCL
-!------Triggering
-      real,parameter :: Su=4e4                  ! Surface unite: celle d'un updraft elementaire
-      real,parameter :: hcoef=1.                ! Coefficient directeur pour le calcul de s2
-      real,parameter :: hmincoef=0.3            ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3)
-      real,parameter :: eps1=0.3                ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
-      real hmin(ngrid)                          ! Ordonnee a l'origine pour le calcul de s2
-      real zmax_moy(ngrid)                      ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
-      real,parameter :: zmax_moy_coef=0.33
-      real depth(ngrid)                         ! Epaisseur moyenne du cumulus
-      real w_max(ngrid)                         ! Vitesse max statistique 
-      real s_max(ngrid)
-!------Closure
-      real pbl_tke_max(ngrid,nlay)              ! Profil de TKE moyenne 
-      real pbl_tke_max0(ngrid)                  ! TKE moyenne au LCL
-      real w_ls(ngrid,nlay)                     ! Vitesse verticale grande echelle (m/s)
-      real,parameter :: coef_m=1.               ! On considere un rendement pour alp_bl_fluct_m
-      real,parameter :: coef_tke=1.             ! On considere un rendement pour alp_bl_fluct_tke
-      
-!!! fin nrlmd le 10/04/2012
-      
-! Nouvelles variables pour la convection
-      integer lalim_conv(ngrid)
-      integer n_int(ngrid)
-      real Ale_bl(ngrid)
-      real Alp_bl(ngrid)
-      real alp_int(ngrid)
-      real dp_int(ngrid),zdp
-      real ale_int(ngrid)
-      real fm_tot(ngrid)
-      real wght_th(ngrid,nlay)
-      
-      CHARACTER*2 str2
-      CHARACTER*10 str10
+      REAL zhla(ngrid,nlay)                     ! TP   plume
+      REAL ztva(ngrid,nlay)                     ! TRPV plume
+      REAL zqsa(ngrid,nlay)                     ! qsat plume
+      
+      REAL zqa(ngrid,nlay,nq)                   ! q    plume (ql=0, qv=qt)
+      
+      REAL linter(ngrid)                        ! Level (continuous) of maximal vertical speed
+      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
+      REAL zmax(ngrid)                          ! Maximal altitude reached by the plume
+      REAL wmax(ngrid)                          ! Maximal vertical speed
+      REAL zw2(ngrid,nlay+1)                    ! Plume vertical speed
+      
+      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
+      
+      REAL f_star(ngrid,nlay+1)                 ! Normalized mass flux
+      REAL entr_star(ngrid,nlay)                ! Normalized entrainment
+      REAL detr_star(ngrid,nlay)                ! Normalized detrainment
+      
+      REAL f(ngrid)                             ! Mass flux norm
+      REAL fm(ngrid,nlay+1)                     ! Mass flux
+      REAL entr(ngrid,nlay)                     ! Entrainment
+      REAL detr(ngrid,nlay)                     ! Detrainment
+      
+      REAL lambda                               ! Time relaxation coefficent
+      
+      REAL zdthladj(ngrid,nlay)                 ! Potential temperature variations
+      REAL dummy(ngrid,nlay)                    ! Dummy argument for thermcell_dq()
       
       CHARACTER (len=20) :: modname='thermcell_main'
       CHARACTER (len=80) :: abort_message
       
-      EXTERNAL SCOPY
-      
-!==============================================================================
+! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
+      INTEGER lalim(ngrid)                      ! Highest alimentation level
+      REAL alim_star(ngrid,nlay)                ! Normalized alimentation
+      REAL alim_star_tot(ngrid)                 ! Integrated alimentation
+      REAL alim_star_clos(ngrid,nlay)           ! Closure alimentation
+      
+!===============================================================================
 ! Initialization
-!==============================================================================
-      
-      seuil = 0.25
+!===============================================================================
       
       IF (firstcall) THEN
-         IF (iflag_thermals==15.or.iflag_thermals==16) THEN
-            dvimpl = 0
-            dqimpl = -1
-         ELSE
-            dvimpl = 1
-            dqimpl = 1
-         ENDIF
-         
          fm0(:,:) = 0.
          entr0(:,:) = 0.
@@ -289,9 +165,30 @@
       ENDIF
       
+      f_star(:,:) = 0.
+      entr_star(:,:) = 0.
+      detr_star(:,:) = 0.
+      
+      f(:) = 0.
+      
       fm(:,:) = 0.
       entr(:,:) = 0.
       detr(:,:) = 0.
-      f(:) = 0.
-      
+      
+      lmax(:) = 1
+      lmix(:) = 1
+      lmin(:) = 1
+      
+      pduadj(:,:) = 0.0
+      pdvadj(:,:) = 0.0
+      pdtadj(:,:) = 0.0
+      pdqadj(:,:,:) = 0.0
+      
+! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
+      alim_star(:,:) = 0.
+      alim_star_tot(:) = 0.
+      alim_star_clos(:,:) = 0.
+      lalim(:) = 1
+      
+! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model!
       DO ig=1,ngrid
          f0(ig) = max(f0(ig), 1.e-2)
@@ -304,33 +201,18 @@
       ENDIF
       
-      wmax_tmp(:) = 0.
-      
-!------------------------------------------------------------------------------
+!===============================================================================
+! Environment settings
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
 ! Calcul de T,q,ql a partir de Tl et qt dans l environnement
-!------------------------------------------------------------------------------
-      
-      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,                        &
-      &                  pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
-      
-!------------------------------------------------------------------------------
-!
-!                       --------------------
-!
-!
-!                       + + + + + + + + + + +
-!
-!
-!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
-!  wh,wt,wo ...
-!
-!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
-!
-!
-!                       --------------------   zlev(1)
-!                       \\\\\\\\\\\\\\\\\\\\
-!
-!------------------------------------------------------------------------------
-! Calcul des altitudes des couches
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+      
+      CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
+      &                  zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
+      
+!-------------------------------------------------------------------------------
+! Levels and layers altitudes
+!-------------------------------------------------------------------------------
       
       DO l=2,nlay
@@ -345,17 +227,9 @@
       ENDDO
       
-!------------------------------------------------------------------------------
-! Calcul de l'epaisseur des couches
-!------------------------------------------------------------------------------
-      
-      DO l=1,nlay
-         deltaz(:,l) = zlev(:,l+1)-zlev(:,l)
-      ENDDO
-      
-!------------------------------------------------------------------------------
-! Calcul des densites
-!------------------------------------------------------------------------------
-      
-      rho(:,:) = pplay(:,:) / (zpspsk(:,:) * RD * ztv(:,:))
+!-------------------------------------------------------------------------------
+! Levels and layers densities
+!-------------------------------------------------------------------------------
+      
+      rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
       
       IF (prt_level.ge.10) THEN
@@ -369,7 +243,7 @@
       ENDDO
       
-!------------------------------------------------------------------------------
-! Calcul de la masse
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! Layers masses
+!-------------------------------------------------------------------------------
       
       DO l=1,nlay
@@ -377,139 +251,105 @@
       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                          |      |
-!    ----- 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
-!  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.
-!------------------------------------------------------------------------------
-      
-      entr_star(:,:) = 0.
-      detr_star(:,:) = 0.
-      alim_star(:,:) = 0.
-      
-      alim_star_tot(:) = 0.
-      
-      lmin(:) = 1
-      
-!==============================================================================
-! Calcul du melange et des variables dans le thermique
-!==============================================================================
-      
-      CALL 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,ztva,                          &
-      &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,      &
-      &    lmin,lev_out,lunout1,igout)
-      
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
-      
-!==============================================================================
+!===============================================================================
+! Explicative schemes
+!===============================================================================
+
+!-------------------------------------------------------------------------------
+! Thermal plume variables
+!-------------------------------------------------------------------------------
+      
+!           top of the model     
+!     ===========================
+!                                
+!     ---------------------------
+!                                        _
+!     ----- F_lmax+1=0 ------zmax         \
+!     lmax                                 |
+!     ------F_lmax>0-------------          |
+!                                          |
+!     ---------------------------          |
+!                                          |
+!     ---------------------------          |
+!                                          |
+!     ------------------wmax,zmix          |
+!     lmix                                 |
+!     ---------------------------          |
+!                                          |
+!     ---------------------------          |
+!                                          | E, D
+!     ---------------------------          |
+!                                          |
+!     --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca
+!         zt, zu, zv, zo, rho              |
+!     ---------------------------          |
+!                                          |
+!     ---------------------------          |
+!                                          |
+!     ---------------------------          |
+!                                          |
+!     ------F_lmin+1>0-----------          |
+!     lmin                                 |
+!     ----- F_lmin=0 ------------        _/
+!                                
+!     ---------------------------
+!                                
+!     ===========================
+!         bottom of the model    
+      
+!-------------------------------------------------------------------------------
+! Zoom on layers k and k-1
+!-------------------------------------------------------------------------------
+      
+!     |     /|\                  |                          |
+!     |----  |  F_k+1 -----------|--------------------------|   level k+1
+!     |      |  w_k+1         |                             |
+!     |                     --|--> D_k                      |
+!     |                       |                             |   layer k
+!     |                    <--|--  E_k                      |
+!     |     /|\               |                             |
+!     |----  |  F_k ----------|-----------------------------|   level k
+!     |      |  w_k        |                                |
+!     |                  --|--> D_k-1                       |
+!     |                    |                                |   layer k-1
+!     |                 <--|--  E_k-1                       |
+!     |     /|\            |                                |
+!     |----  |  F_k-1 -----|--------------------------------|   level k-1
+!            |  w_k-1                                        
+!     0                  fraca                              1
+!      \__________________/ \______________________________/
+!          plume (fraca)          environment (1-fraca)   
+      
+!===============================================================================
+! Thermal plumes computation
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
+! Thermal plumes speeds, fluxes, tracers and temperatures
+!-------------------------------------------------------------------------------
+      
+      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
+      &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
+      &                    detr_star,entr_star,f_star,                        &
+      &                    ztva,zhla,zqla,zqta,zta,                           &
+      &                    zw2,zqsa,lmix,lmin)
+      
+!-------------------------------------------------------------------------------
 ! Thermal plumes characteristics: zmax, zmix, wmax
-!==============================================================================
-      
-      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,            &
-      &                     zlev,lmax,zmax,zmix,wmax,f_star,lev_out)
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : WARNING: zw2 became its square root in thermcell_height!
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
-      
-!==============================================================================
+!-------------------------------------------------------------------------------
+      
+! AB: Careful, zw2 became its square root in thermcell_height!
+      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
+      &                     zlev,lmax,zmax,zmix,wmax,f_star)
+      
+!===============================================================================
 ! Closure and mass fluxes computation
-!==============================================================================
-      
-      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,                  &
-      &                  lalim,lmin,zmax_sec,wmax_sec,lev_out)
-      
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! 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(:,:)
-      
-!------------------------------------------------------------------------------
-! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2)
-!------------------------------------------------------------------------------
-      
-      IF (iflag_thermals_closure.eq.1) THEN
-         CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                &
-         &                      lalim,alim_star_clos,f_star,                  &
-         &                      zmax_sec,wmax_sec,f,lev_out)
-      ELSEIF (iflag_thermals_closure.eq.2) THEN
-         CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                &
-         &                      lalim,alim_star,f_star,                       &
-         &                      zmax,wmax,f,lev_out)
-      ELSE
-         print *, 'ERROR: no closure had been selected!'
-         call abort
-      ENDIF
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
+! Closure
+!-------------------------------------------------------------------------------
+      
+      CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
+      &                      lmax,entr_star,zmax,wmax,f)
       
       IF (tau_thermals>1.) THEN
@@ -520,7 +360,7 @@
       ENDIF
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! Test valable seulement en 1D mais pas genant
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       IF (.not. (f0(1).ge.0.) ) THEN
          abort_message = '.not. (f0(1).ge.0.)'
@@ -529,21 +369,19 @@
       ENDIF
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Mass fluxes
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
-      &                   lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
-      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla,               &
-      &                   lev_out,lunout1,igout)
-      
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
-!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
+!     That is not already done for thermcell_flux.
+      &                   lalim,lmin,lmax,entr_star,detr_star,                &
+      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
+      
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! 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).
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! tau_thermals (typiquement 1800s sur Terre).
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       
       IF (tau_thermals>1.) THEN
@@ -558,7 +396,7 @@
       ENDIF
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Updraft fraction
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO ig=1,ngrid
@@ -577,417 +415,41 @@
       ENDDO
       
-!==============================================================================
+!===============================================================================
 ! Transport vertical
-!==============================================================================
-      
-!------------------------------------------------------------------------------
-! Calcul du transport vertical (de qt et tp)
-!------------------------------------------------------------------------------
-      
-      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
-      &                 zthl,zdthladj,zta,lmin,lev_out)
-      
-      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
-      &                 po,pdoadj,zoa,lmin,lev_out)
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
+! Calcul du transport vertical de la temperature potentielle
+!-------------------------------------------------------------------------------
+      
+      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
+      &                 zhl,zdthladj,dummy,lmin)
       
       DO l=1,nlay
          DO ig=1,ngrid
-           pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l)  
+            pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l)  
          ENDDO
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! Calcul du transport vertical des traceurs
+!-------------------------------------------------------------------------------
+      
+      DO iq=1,nq
+         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,        &
+         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
+      ENDDO
+      
+!-------------------------------------------------------------------------------
 ! Calcul du transport vertical du moment horizontal
-!------------------------------------------------------------------------------
-      
-      IF (dvimpl.eq.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,       &
-         &                  zmax,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
-      ELSE
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! Calcul purement conservatif pour le transport de V=(u,v)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
-         &                 zu,pduadj,zua,lmin,lev_out)
-         
-         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
-         &                 zv,pdvadj,zva,lmin,lev_out)
-      ENDIF
-      
-!==============================================================================
-! Calculs de diagnostiques pour les sorties
-!==============================================================================
-      
-      IF (sorties) THEN
-         
-!------------------------------------------------------------------------------
-! Calcul du niveau de condensation
-!------------------------------------------------------------------------------
-         
-         DO ig=1,ngrid
-            nivcon(ig) = 0
-            zcon(ig) = 0.
-         ENDDO 
-!nouveau calcul
-         do ig=1,ngrid
-            CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
-            pcon(ig) = pplay(ig,1)*(zo(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
-         
-         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
-              write(*,*) 'ERROR: thermal plumes rise the model top'
-              CALL abort
-         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
-!------------------------------------------------------------------------------
-         
-         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'
-         
-         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.-po(ig,l)*1000.)**2
-!test: on calcul q2/po=ratqsc
-               ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(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.-po(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
-         
-         CALL thermcell_alp(ngrid,nlay,ptimestep,                             &
-         &                  pplay,pplev,                                      &
-         &                  fm0,entr0,lmax,                                   &
-         &                  Ale_bl,Alp_bl,lalim_conv,wght_th,                 &
-         &                  zw2,fraca,pcon,                                   &
-         &                  rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,    &
-         &                  pbl_tke,pctsrf,omega,airephy,                     &
-         &                  zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,&
-         &                  n2,s2,ale_bl_stat,                                &
-         &                  therm_tke_max,env_tke_max,                        &
-         &                  alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,       &
-         &                  alp_bl_conv,alp_bl_stat)
-         
-!------------------------------------------------------------------------------
-! Calcul du ratqscdiff
-!------------------------------------------------------------------------------
-      
-         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) / (po(ig,l) * 1000.)   
-            ENDDO
-         ENDDO
-         
-      ENDIF
+!-------------------------------------------------------------------------------
+      
+      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
+      &                 zu,pduadj,zua,lmin)
+      
+      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
+      &                 zv,pdvadj,zva,lmin)
       
       
 RETURN
 END
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,              &
-                       ztva,zqla,f_star,zw2,comment)
-      
-      
-      USE print_control_mod, ONLY: prt_level
-      
-      IMPLICIT NONE
-      
-      
-!==============================================================================
-! Declaration
-!==============================================================================
-      
-!      inputs:
-!      -------
-      
-      INTEGER ngrid
-      INTEGER nlay
-      INTEGER long(ngrid)
-      
-      REAL pplev(ngrid,nlay+1)
-      REAL pplay(ngrid,nlay)
-      REAL ztv(ngrid,nlay)
-      REAL po(ngrid,nlay)
-      REAL ztva(ngrid,nlay)
-      REAL zqla(ngrid,nlay)
-      REAL f_star(ngrid,nlay)
-      REAL zw2(ngrid,nlay)
-      REAL seuil
-      
-      CHARACTER*21 comment
-      
-!      local:
-!      ------
-      
-      INTEGER i, k
-      
-!==============================================================================
-! Test
-!==============================================================================
-      
-      IF (prt_level.ge.1) THEN
-         write(*,*) 'WARNING: in 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*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
-            ENDDO
-        ENDIF
-      ENDDO
-      
-      
-RETURN
-END
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,            &
-                                   rg,pplev,therm_tke_max)
-      
-      
-!==============================================================================
-!   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
-!
-!   Transport de la TKE par le thermique moyen pour la fermeture en ALP
-!   On transporte pbl_tke pour donner therm_tke
-!==============================================================================
-      
-      USE print_control_mod, ONLY: prt_level
-      
-      IMPLICIT NONE
-      
-      
-!==============================================================================
-! Declarations
-!==============================================================================
-      
-      integer ngrid,nlay,nsrf
-      
-      real ptimestep
-      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
-      real entr0(ngrid,nlay),rg
-      real therm_tke_max(ngrid,nlay)
-      real detr0(ngrid,nlay)
-      
-      real masse(ngrid,nlay),fm(ngrid,nlay+1)
-      real entr(ngrid,nlay)
-      real q(ngrid,nlay)
-      
-      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
-      
-      real zzm
-      
-      integer ig,k
-      integer isrf
-      
-!------------------------------------------------------------------------------
-! 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.
-      
-!!! nrlmd le 16/09/2010
-!   calcul de la valeur dans les ascendances
-!      do ig=1,ngrid
-!         qa(ig,1) = q(ig,1)
-!      ENDDO
-      
-      q(:,:)=therm_tke_max(:,:)
-      
-      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.) print *, 'WARNING: qa is negative!'
-!               IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
-            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.) print*,'WARNING: wqd is negative!'
-            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) + ptimestep / masse(ig,k)                    &
-               &       * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k)        &
-               &       - wqd(ig,k) + wqd(ig,k+1))
-            ENDDO
-         ENDDO
-      ENDIF
-      
-      therm_tke_max(:,:) = q(:,:)
-      
-      
-RETURN
-END
-
