Index: trunk/LMDZ.GENERIC/libf/phystd/convadj.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/convadj.F	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/convadj.F	(revision 2127)
@@ -197,5 +197,5 @@
             l2 = l2 + 1
             IF (l2 .GT. nlay) EXIT
-            IF ((zhc(i, l2).LT.zhc(i, l2-1)).and.(l.GT.lmax(i))) THEN
+            IF ((zhc(i, l2).LT.zhc(i, l2-1)).and.(l2.GT.lmax(i))) THEN
  
 !     l2 is the highest level of the unstable column
@@ -226,5 +226,5 @@
                 down = .false.
                 IF (l1 .ne. 1) then    !--  and then
-                  IF ((zhmc.LT.zhc(i, l1-1)).and.(l.GT.lmax(i))) then
+                  IF ((zhmc.LT.zhc(i, l1-1)).and.(l1.GT.lmax(i))) then
                     down = .true.
                   END IF
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 2127)
@@ -1190,5 +1190,6 @@
       
       IF (calltherm) THEN
-! AB : We need to evaporate ice before calling thermcell_main
+         
+! AB: We need to evaporate ice before calling thermcell_main.
          IF (water) THEN
             CALL evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,zqtherm,zttherm)
@@ -1198,24 +1199,13 @@
          ENDIF
          
-! AB:  If a plume stops, the parametrization never look above if somewhere the atmosphere is still unstable!
-!      Maybe it's a good idea to call convective adjustment too.
-         CALL thermcell_main(icount, ngrid, nlayer, ptimestep,                   &
-                             pplay, pplev, pphi, firstcall,                      &
-                             pu, pv, zttherm, zqtherm(:,:,igcm_h2o_vap),         &
-                             zdutherm, zdvtherm, zdttherm, zdotherm,             &
-                             f0, fm0, entr0, detr0,                              &
+! AB: If a plume stops, the parametrization never look above if somewhere the atmosphere is still unstable!
+!     Maybe it's a good idea to call convective adjustment too.
+         CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall,            &
+                             pplay, pplev, pphi, zpopsk,                         &
+                             pu, pv, zttherm, zqtherm,                           &
+                             zdutherm, zdvtherm, zdttherm, zdqtherm,             &
+                             f0, fm0, entr0, detr0, zw2, fraca,                  &
                              zqta, zqla, ztv, ztva, ztla, zthl, zqsatth,         &
-                             zw2, fraca,                                         &
-                             lmin, lmix, lalim, lmax,                            &
-                             zpopsk, ratqscth, ratqsdiff,                        &
-! AB : Next variables are only used for diagnoses
-                             Ale_bl,Alp_bl,lalim_conv,wght_th,                   &
-                             pbl_tke,pctsrf,omega_therm,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)
+                             lmin, lmix, lalim, lmax)
          
          DO ig=1,ngrid
@@ -1228,11 +1218,11 @@
          pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer)
          
-         IF(tracer) THEN
-            pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) + zdotherm(1:ngrid,1:nlayer)
+         IF (tracer) THEN
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
          ENDIF
          
       ELSE
          
-         lmax(1:ngrid) = 0
+         lmax(1:ngrid) = 1
          
       ENDIF ! end of 'calltherm'
@@ -2358,7 +2348,6 @@
          CALL send_xios_field('entr',entr0)
          CALL send_xios_field('detr',detr0)
-         CALL send_xios_field('dt_plm',zdttherm)
-         CALL send_xios_field('pmax',pmax)
       ENDIF
+      
       IF (water) THEN
          CALL send_xios_field('h2o_vap',zq(:,:,igcm_h2o_vap))
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90	(revision 2127)
@@ -2,10 +2,9 @@
 !
 !
-      SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,             &
-      &                            lalim,alim_star,f_star,                    &
-      &                            zmax,wmax,f,lev_out)
+SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
+                             lmax,alim_star,zmax,wmax,f)
       
-!==============================================================================
-!thermcell_closure: fermeture, determination de f
+!===============================================================================
+!  Purpose: fermeture, determination de f
 !
 ! Modification 7 septembre 2009
@@ -15,5 +14,5 @@
 ! l'idee etant que le choix se fasse a l'appel de thermcell_closure
 ! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
-!==============================================================================
+!===============================================================================
       
       USE thermcell_mod
@@ -22,32 +21,30 @@
       
       
-!=============================================================================
+!===============================================================================
 ! Declaration
-!=============================================================================
+!===============================================================================
       
-!      inputs:
-!      -------
+!     Inputs:
+!     -------
       
-      INTEGER ngrid,nlay
-      INTEGER lalim(ngrid)
-      INTEGER lev_out                           ! niveau pour les print
+      INTEGER ngrid, nlay
+      INTEGER lmax(ngrid)
       
-      REAL alim_star(ngrid,nlay)
-      REAL f_star(ngrid,nlay+1)
+      REAL ptimestep
       REAL rho(ngrid,nlay)
       REAL zlev(ngrid,nlay)
+      REAL alim_star(ngrid,nlay)
       REAL zmax(ngrid)
       REAL wmax(ngrid)
-      REAL ptimestep
       
-!      outputs:
-!      --------
+!     Outputs:
+!     --------
       
       REAL f(ngrid)
       
-!      local:
-!      ------
+!     Local:
+!     ------
       
-      INTEGER ig,k 
+      INTEGER ig, l 
       INTEGER llmax
       
@@ -56,7 +53,7 @@
       REAL zdenom(ngrid)
       
-!==============================================================================
+!===============================================================================
 ! Initialization
-!==============================================================================
+!===============================================================================
       
       alim_star2(:) = 0.
@@ -67,28 +64,28 @@
       llmax = 1
       
-!==============================================================================
+!===============================================================================
 ! Closure
-!==============================================================================
+!===============================================================================
       
-!------------------------------------------------------------------------------
-! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! Indice vertical max atteint par les thermiques sur le domaine
+!-------------------------------------------------------------------------------
       
       DO ig=1,ngrid
-         IF (lalim(ig)>llmax) THEN
-            llmax = lalim(ig)
+         IF (lmax(ig).GT.llmax) THEN
+            llmax = lmax(ig)
          ENDIF
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Calcul des integrales sur la verticale de alim_star et de alim_star^2/(rho dz)
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
-      DO k=1,llmax-1
+      DO l=1,llmax-1
          DO ig=1,ngrid
-            IF (k<lalim(ig)) THEN
-               alim_star2(ig) = alim_star2(ig) + alim_star(ig,k)**2           &
-               &              / (rho(ig,k) * (zlev(ig,k+1) - zlev(ig,k)))
-               alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,k)
+            IF (l < lmax(ig)) THEN
+               alim_star2(ig) = alim_star2(ig) + alim_star(ig,l)**2           &
+               &              / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l)))
+               alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
             ENDIF
          ENDDO
@@ -96,10 +93,13 @@
       
       DO ig=1,ngrid
-         IF (alim_star2(ig)>1.e-10) THEN
+         IF (alim_star2(ig).GT.0.) THEN
             f(ig) = wmax(ig) * alim_star_tot(ig)                              &
-            &     / (max(500.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
+            &     / (max(1.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
+         ELSE
+            f(ig) = 0.
          ENDIF
       ENDDO
       
+      
 RETURN
 END
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90	(revision 2127)
@@ -2,39 +2,38 @@
 !
 !
-SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse,         &
-                        q,dq,qa,lmin,lev_out)
+SUBROUTINE thermcell_dq(ngrid,nlay,ptimestep,fm,entr,detr,masse,              &
+                        q,dq,qa,lmin)
       
       
-!==============================================================================
-!  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
-!
+!===============================================================================
+!  Purpose: 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
+!  
 !  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
-!
-!  Modif 2018/11/05 (AB alexandre.boissinot@lmd.jussieu.fr)
+!  Introduction of an implicit computation of vertical advection in the environ-
+!     ment of thermal plumes in thermcell_dq
 !  
-!
-!==============================================================================
+!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
+!     dqimpl = 1 : implicit scheme
+!     dqimpl = 0 : explicit scheme
+!  
+!===============================================================================
       
       USE print_control_mod, ONLY: prt_level
+      USe thermcell_mod, ONLY: dqimpl
       
       IMPLICIT NONE
       
       
-!==============================================================================
+!===============================================================================
 ! Declaration
-!==============================================================================
+!===============================================================================
       
-!      inputs:
-!      -------
+!     Inputs:
+!     -------
       
       INTEGER ngrid, nlay
-      INTEGER impl
       INTEGER lmin(ngrid)
-      INTEGER lev_out                           ! niveau pour les print
       
       REAL ptimestep
@@ -44,13 +43,13 @@
       REAL detr(ngrid,nlay)
       REAL q(ngrid,nlay)
+      
+!     Outputs:
+!     --------
+      
+      REAL dq(ngrid,nlay)
       REAL qa(ngrid,nlay)
       
-!      outputs:
-!      --------
-      
-      REAL dq(ngrid,nlay)
-      
-!      local:
-!      ------
+!     Local:
+!     ------
       
       INTEGER ig, l
@@ -63,25 +62,17 @@
       REAL zzm
       
-      CHARACTER (LEN=20) :: modname
-      CHARACTER (LEN=80) :: abort_message
+!===============================================================================
+! Initialization
+!===============================================================================
       
-!==============================================================================
-! Initialization
-!==============================================================================
+      qold(:,:) = q(:,:)
       
-      qold = q
+!===============================================================================
+! Tracer variation computation
+!===============================================================================
       
-      modname = 'thermcell_dq'
-      
-      IF (impl.lt.0) THEN
-         print *, 'OLD EXPLICIT SCHEME IS USED'
-         CALL thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,               &
-         &                   masse,q,dq,qa,lev_out)
-         RETURN
-      ENDIF
-      
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! CFL criterion computation for advection in downdraft
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       cfl = 0.
@@ -103,13 +94,12 @@
                print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
                print *, 'fm-', fm(ig,l-1)
-               abort_message = 'entr dt > m, 1st'
-               CALL abort_physic(modname,abort_message,1)
+               CALL abort
             ENDIF
          ENDDO
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Computation of tracer concentrations in the ascending plume
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO ig=1,ngrid
@@ -121,10 +111,10 @@
       DO ig=1,ngrid
          DO l=lmin(ig)+1,nlay
-            if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then
+            IF ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-6*masse(ig,l)) THEN
                qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l))      &
                &        / (fm(ig,l+1) + detr(ig,l))
-            else
+            ELSE
                qa(ig,l) = q(ig,l)
-            endif
+            ENDIF
             
 !            IF (qa(ig,l).lt.0.) THEN
@@ -140,7 +130,7 @@
       ENDDO
       
-!------------------------------------------------------------------------------
-! Plume vertical flux of q
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! Plume vertical flux of tracer
+!-------------------------------------------------------------------------------
       
       DO l=2,nlay-1
@@ -151,24 +141,29 @@
       fqa(:,nlay) = 0.
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Trace species evolution
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
-      IF (impl==0) THEN
-         do l=1,nlay-1
+      IF (dqimpl==0) THEN
+         DO l=1,nlay-1
             q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l)       &
             &      + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
-         enddo
-      ELSE
-         do l=nlay-1,1,-1
+         ENDDO
+      ELSEIF (dqimpl==1) THEN
+         DO l=nlay-1,1,-1
             q(:,l) = ( q(:,l) + ptimestep / masse(:,l)                        &
             &      * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) )       &
             &      / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
          ENDDO
+      ELSE
+         print *, 'ERROR: No corresponding scheme for mixing computations!'
+         print *, '       dqimpl must be equal to 1, 0 or -1 but'
+         print *, 'dqimpl =', dqimpl
+         call abort
       ENDIF
       
-!==============================================================================
+!===============================================================================
 ! Tendencies
-!==============================================================================
+!===============================================================================
       
       DO l=1,nlay
@@ -182,192 +177,2 @@
 RETURN
 END
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
-Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse,            &
-                          q,dq,qa,lev_out)
-      
-      
-!==============================================================================
-!
-!   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
-!
-!==============================================================================
-      
-      USE print_control_mod, ONLY: prt_level
-      
-      implicit none
-      
-      
-!==============================================================================
-! Declaration
-!==============================================================================
-      
-      integer ngrid,nlay,impl
-
-      real ptimestep
-      real masse(ngrid,nlay),fm(ngrid,nlay+1)
-      real entr(ngrid,nlay)
-      real q(ngrid,nlay)
-      real dq(ngrid,nlay)
-      integer lev_out                           ! niveau pour les print
-
-      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
-
-      real zzm
-
-      integer ig,l
-      real cfl
-
-      real qold(ngrid,nlay)
-      real ztimestep
-      integer niter,iter
-      CHARACTER (LEN=20) :: modname='thermcell_dq'
-      CHARACTER (LEN=80) :: abort_message
-
-
-
-! Calcul du critere CFL pour l'advection dans la subsidence
-      cfl = 0.
-      do l=1,nlay
-         do ig=1,ngrid
-            zzm=masse(ig,l)/ptimestep
-            cfl=max(cfl,fm(ig,l)/zzm)
-            if (entr(ig,l).gt.zzm) then
-               print*,'entr*dt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l)
-               abort_message = 'entr dt > m, 2nd'
-               CALL abort_physic (modname,abort_message,1)
-            endif
-         enddo
-      enddo
-      
-!IM 090508     print*,'CFL CFL CFL CFL ',cfl
-      
-#undef CFL
-#ifdef CFL
-! On subdivise le calcul en niter pas de temps.
-      niter=int(cfl)+1
-#else
-      niter=1
-#endif
-      
-      ztimestep=ptimestep/niter
-      qold=q
-      
-      do iter=1,niter
-         if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
-         
-!   calcul du detrainement
-         do l=1,nlay
-            do ig=1,ngrid
-               detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l)
-!               print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l)
-!test
-               if (detr(ig,l).lt.0.) then
-                  entr(ig,l)=entr(ig,l)-detr(ig,l)
-                  detr(ig,l)=0.
-!                  print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l),
-!     s            'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l)
-               endif
-               
-!               if (fm(ig,l+1).lt.0.) then
-!                  print*,'fm2<0!!!'
-!               endif
-               
-!               if (entr(ig,l).lt.0.) then
-!                  print*,'entr2<0!!!'
-!               endif
-            enddo
-         enddo
-         
-! calcul de la valeur dans les ascendances
-        do ig=1,ngrid
-            qa(ig,1)=q(ig,1)
-         enddo
-         
-         do l=2,nlay
-            do ig=1,ngrid
-               if ((fm(ig,l+1)+detr(ig,l))*ztimestep.gt.1.e-5*masse(ig,l)) then
-                  qa(ig,l) = (fm(ig,l)*qa(ig,l-1)+entr(ig,l)*q(ig,l))         &
-                  &        / (fm(ig,l+1)+detr(ig,l))
-               else
-                  qa(ig,l)=q(ig,l)
-               endif
-               
-               if (qa(ig,l).lt.0.) then
-!                  print*,'qa<0!!!'
-               endif
-               
-               if (q(ig,l).lt.0.) then
-!                  print*,'q<0!!!'
-               endif
-            enddo
-         enddo
-         
-! Calcul du flux subsident
-         
-         do l=2,nlay
-            do ig=1,ngrid
-#undef centre
-#ifdef centre
-                wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l))
-#else
-
-#define plusqueun
-#ifdef plusqueun
-! Schema avec advection sur plus qu'une maille.
-               zzm=masse(ig,l)/ztimestep
-               if (fm(ig,l)>zzm) then
-                 wqd(ig,l)=zzm*q(ig,l)+(fm(ig,l)-zzm)*q(ig,l+1)
-               else
-                  wqd(ig,l)=fm(ig,l)*q(ig,l)
-               endif
-#else
-               wqd(ig,l)=fm(ig,l)*q(ig,l)
-#endif
-#endif
-
-               if (wqd(ig,l).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 l=1,nlay
-            do ig=1,ngrid
-               q(ig,l)=q(ig,l)+(detr(ig,l)*qa(ig,l)-entr(ig,l)*q(ig,l)           &
-               &        -wqd(ig,l)+wqd(ig,l+1))                                  &
-               &        *ztimestep/masse(ig,l)
-!               if (dq(ig,l).lt.0.) then
-!                  print*,'dq<0!!!'
-!               endif
-           enddo
-         enddo
-      enddo
-      
-! Calcul des tendances
-      do l=1,nlay
-         do ig=1,ngrid
-            dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
-            q(ig,l) = qold(ig,l)
-         enddo
-      enddo
-      
-      
-return
-end
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90	(revision 2127)
@@ -2,107 +2,107 @@
 !
 !
-      SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
-      &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
-
-!==============================================================================
-! thermcell_env: calcule les caracteristiques de l environnement
-!                necessaires au calcul des proprietes dans le thermique
-!==============================================================================
-
+SUBROUTINE thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
+                         zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
+      
+      
+!===============================================================================
+!  Purpose: calcul des caracteristiques de l'environnement necessaires au calcul
+!           des proprietes dans le thermique.
+!  
+!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
+!  
+!  Nota Bene
+!     ql   means "non-gaseous water mass mixing ratio" (liquid and solid)
+!     qv   means "vapor mass mixing ratio"
+!     qt   means "total water mass mixing ratio"
+!     TP   means "potential temperature"
+!     TRPV means "virtual potential temperature with latent heat release"  
+!     TPV  means "virtual potential temperature"
+!     TR   means "temperature with latent heat release"
+!===============================================================================
+      
       USE print_control_mod, ONLY: prt_level
       USE thermcell_mod, ONLY: RKAPPA
       USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
-      USE planete_mod, ONLY: preff
+      USE tracer_h, ONLY: igcm_h2o_vap, igcm_h2o_ice
+      USE callkeys_mod, ONLY: water
       
-      IMPLICIT NONE   
+      IMPLICIT NONE
       
-!==============================================================================
+      
+!===============================================================================
 ! Declaration
-!==============================================================================
+!===============================================================================
       
-!      inputs:
-!      -------
+!     Inputs:
+!     -------
       
-      INTEGER ngrid
-      INTEGER nlay
-      INTEGER lev_out                           ! niveau pour les print
+      INTEGER ngrid, nlay, nq
       
-      REAL po(ngrid,nlay)                       ! large scale water
-      REAL pt(ngrid,nlay)                       ! large scale temperature
-      REAL pu(ngrid,nlay)                       ! large scale zonal wind
-      REAL pv(ngrid,nlay)                       ! large scale meridional wind
-      REAL pplay(ngrid,nlay)                    ! layers pressure
-      REAL pplev(ngrid,nlay+1)                  ! levels pressure
+      REAL pq(ngrid,nlay,nq)                    ! Large scale water
+      REAL pt(ngrid,nlay)                       ! Large scale temperature
+      REAL pu(ngrid,nlay)                       ! Large scale zonal wind
+      REAL pv(ngrid,nlay)                       ! Large scale meridional wind
+      REAL pplay(ngrid,nlay)                    ! Layers pressure
+      REAL pplev(ngrid,nlay+1)                  ! Levels pressure
+      REAL zpopsk(ngrid,nlay)                   ! Exner function
       
-!      outputs:
-!      --------
+!     Outputs:
+!     --------
       
-      REAL zo(ngrid,nlay)                       ! qt  environment
-      REAL zl(ngrid,nlay)                       ! ql  environment
-      REAL zh(ngrid,nlay)                       ! T   environment
-      REAL ztv(ngrid,nlay)                      ! TV  environment
-      REAL zthl(ngrid,nlay)                     ! TPV environment
-      REAL zu(ngrid,nlay)                       ! u   environment
-      REAL zv(ngrid,nlay)                       ! v   environment
+      REAL zqt(ngrid,nlay)                      ! q    environment
+      REAL zql(ngrid,nlay)                      ! q    environment
+      REAL zt(ngrid,nlay)                       ! T    environment
+      REAL ztv(ngrid,nlay)                      ! TRPV environment
+      REAL zhl(ngrid,nlay)                      ! TP   environment
+      REAL zu(ngrid,nlay)                       ! u    environment
+      REAL zv(ngrid,nlay)                       ! v    environment
+      REAL zqs(ngrid,nlay)                      ! qsat environment
       
-      REAL zpspsk(ngrid,nlay)                   ! Exner function
-      REAL pqsat(ngrid,nlay)                    ! saturation vapor pressure
+!     Local:
+!     ------
       
-!      local:
-!      ------
+      INTEGER ig, l
       
-      INTEGER ig, ll
+      REAL psat                                 ! Dummy argument for Psat_water()
       
-      REAL dummy
+!===============================================================================
+! Initialization
+!===============================================================================
       
-!==============================================================================
-! Initialization
-!==============================================================================
+      zu(:,:) = pu(:,:)
+      zv(:,:) = pv(:,:)
       
-!------------------------------------------------------------------------------
-! Calcul des caracteristiques de l'environnement
-!------------------------------------------------------------------------------
-      DO ll=1,nlay
-         DO ig=1,ngrid
-            zo(ig,ll) = po(ig,ll)
-            zl(ig,ll) = 0.
-            zh(ig,ll) = pt(ig,ll)
+      zhl(:,:) = pt(:,:) / zpopsk(:,:)
+      
+      zqt(:,:) = pq(:,:,igcm_h2o_vap)
+      
+!===============================================================================
+! Condensation and latent heat release
+!===============================================================================
+      
+      IF (water) THEN
+         
+         DO l=1,nlay
+            DO ig=1,ngrid
+               CALL Psat_water(pt(ig,l), pplev(ig,l), psat, zqs(ig,l))
+            ENDDO
          ENDDO
-      ENDDO
-      
-!==============================================================================
-! Computation of qsat and condensation
-!==============================================================================
-
-      DO ll=1,nlay
-         DO ig=1,ngrid
-            CALL Psat_water(pt(ig,ll), pplev(ig,ll), dummy, pqsat(ig,ll))
+         
+         DO l=1,nlay
+            DO ig=1,ngrid
+               zql(ig,l) = max(0.,pq(ig,l,igcm_h2o_vap) - zqs(ig,l))
+               zt(ig,l) = pt(ig,l) + RLvCp * zql(ig,l)
+               ztv(ig,l) = zt(ig,l) / zpopsk(ig,l)                            &
+               &         * (1. + RETV * (zqt(ig,l)-zql(ig,l)) - zql(ig,l))
+            ENDDO
          ENDDO
-      ENDDO
-      
-      DO ll=1,nlay
-         DO ig=1,ngrid
-            zl(ig,ll) = max(0.,po(ig,ll) - pqsat(ig,ll))                      ! zl   set to ql   env
-            zh(ig,ll) = pt(ig,ll) + RLvCp * zl(ig,ll)                         ! zh   set to TR   env
-            zo(ig,ll) = po(ig,ll) - zl(ig,ll)                                 ! zo   set to qv   env
-         ENDDO
-      ENDDO
-      
-      DO ll=1,nlay
-         DO ig=1,ngrid
-             zpspsk(ig,ll) = (pplay(ig,ll)/preff)**RKAPPA
-             zu(ig,ll) = pu(ig,ll)
-             zv(ig,ll) = pv(ig,ll)
-             
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : WARNING, zh is no longer the potentiel temperature.
-!      It is now the temperature. How awful!
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            ztv(ig,ll) = zh(ig,ll)/zpspsk(ig,ll)                              ! ztv  set to TRP  env
-            ztv(ig,ll) = ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll))           ! ztv  set to TRPV env
-            zthl(ig,ll) = pt(ig,ll)/zpspsk(ig,ll)                             ! zthl set to TP   env
-            
-         ENDDO
-      ENDDO
+         
+      ELSE
+         
+         zt(:,:) = pt(:,:)
+         ztv(:,:) = pt(:,:) / zpopsk(:,:)
+         
+      ENDIF
       
       
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2127)
@@ -3,12 +3,14 @@
 !
 SUBROUTINE 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)
-      
-      
-!==============================================================================
-! thermcell_flux: deduction des flux
-!==============================================================================
+                          lalim,lmin,lmax,entr_star,detr_star,                &
+                          f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
+      
+      
+!===============================================================================
+!  Purpose: deduction des flux
+!  
+!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
+!  
+!===============================================================================
       
       USE print_control_mod, ONLY: prt_level
@@ -18,21 +20,18 @@
       
       
-!==============================================================================
+!===============================================================================
 ! Declaration
-!==============================================================================
-      
-!      inputs:
-!      ------- 
-      
-      INTEGER ngrid
-      INTEGER nlay
-      INTEGER igout
-      INTEGER lev_out
-      INTEGER lunout1
+!===============================================================================
+      
+!     Inputs:
+!     ------- 
+      
+      INTEGER ngrid, nlay
       INTEGER lmin(ngrid)
       INTEGER lmax(ngrid)
       INTEGER lalim(ngrid)
       
-      REAL alim_star(ngrid,nlay)
+! AB : I remove alimentation, which is in fact included in entr_star. EN COURS
+!      REAL alim_star(ngrid,nlay)
       REAL entr_star(ngrid,nlay)
       REAL detr_star(ngrid,nlay)
@@ -46,6 +45,6 @@
       REAL zmax(ngrid)
       
-!      outputs:
-!      --------      
+!     Outputs:
+!     --------      
       
       REAL entr(ngrid,nlay)
@@ -53,17 +52,17 @@
       REAL fm(ngrid,nlay+1)
       
-!      local:
-!      ------
-      
-      INTEGER ig,l,k
-      INTEGER lout
-      
-      REAL zfm                            ! mass flux such as updraft fraction is equal to its maximal authorized value alphamax
-      REAL f_old                          ! save initial value of mass flux
-      REAL eee0                           ! save initial value of entrainment
-      REAL ddd0                           ! save initial value of detrainment
+!     Local:
+!     ------
+      
+      INTEGER ig, l, k
+      INTEGER igout, lout                 ! Error grid point number and level
+      
+      REAL zfm                            ! Mass flux such as updraft fraction is equal to its maximal authorized value alphamax
+      REAL f_old                          ! Save initial value of mass flux
+      REAL eee0                           ! Save initial value of entrainment
+      REAL ddd0                           ! Save initial value of detrainment
       REAL eee                            ! eee0 - layer mass * maximal authorized mass fraction / dt
       REAL ddd                            ! ddd0 - eee
-      REAL zzz                            ! temporary variable set to mass flux
+      REAL zzz                            ! Temporary variable set to mass flux
       REAL fact
       REAL test
@@ -83,7 +82,7 @@
       CHARACTER (len=80) :: abort_message
       
-!==============================================================================
+!===============================================================================
 ! Initialization
-!==============================================================================
+!===============================================================================
       
       ncorecfm1 = 0
@@ -100,32 +99,34 @@
       fm(:,:)   = 0.
       
-      IF (prt_level.ge.10) THEN
-         write(lunout1,*) 'Dans thermcell_flux 0'
-         write(lunout1,*) 'flux base ', f(igout)
-         write(lunout1,*) 'lmax ', lmax(igout)
-         write(lunout1,*) 'lalim ', lalim(igout)
-         write(lunout1,*) 'ig= ', igout
-         write(lunout1,*) ' l E*    A*     D*  '
-         write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),                  &
-         &                             alim_star(igout,l),detr_star(igout,l), &
-         &                             l=1,lmax(igout))
-      ENDIF
-      
-!==============================================================================
+!===============================================================================
 ! Calcul de l'entrainement, du detrainement et du flux de masse
-!==============================================================================
-      
-!------------------------------------------------------------------------------
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
 ! Multiplication par la norme issue de la relation de fermeture
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO l=1,nlay
-         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
+! AB : I remove alimentation, which is in fact included in entr_star. EN COURS
+!         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
+         entr(:,l) = f(:) * entr_star(:,l)
          detr(:,l) = f(:) * detr_star(:,l)
       ENDDO
       
-!------------------------------------------------------------------------------
+! AB : temporary test
+      DO l=1,nlay
+         DO ig=1,ngrid
+            IF (detr(ig,l)<0.) THEN
+               print *, 'ERROR: detr is negative in thermcell_flux!'
+               print *, 'ig,l', ig, l
+               print *, 'detr_star', detr_star(ig,l)
+               print *, 'detr', detr(ig,l)
+            ENDIF
+         ENDDO
+      ENDDO
+      
+!-------------------------------------------------------------------------------
 ! Calcul du flux de masse
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO l=1,nlay
@@ -144,15 +145,15 @@
       ENDDO
       
-!==============================================================================
+!===============================================================================
 ! Tests de validite
-!==============================================================================
+!===============================================================================
       
       DO l=1,nlay
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Verification de la positivite du flux de masse entrant
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : I remove the correction and replace it by an uncompromising test.
 !      According to the previous derivations, we MUST have positive mass flux
@@ -160,5 +161,5 @@
 !      Then the only value which can be negative is the mass flux at the top of
 !      the plume. However, this value was set to zero a few lines above.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          labort_physic=.false.
@@ -177,25 +178,24 @@
             print *, 'lmin,lmax', lmin(igout), lmax(igout)
             print *, '-------------------------------'
-            print *, 'fm+', fm(igout,lout+1)
-            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
-            print *, 'fm ', fm(igout,lout)
-            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
-            print *, 'fm-', fm(igout,lout-1)
+            DO k=nlay,1,-1
+               print *, 'fm ', fm(igout,k+1)
+               print *, 'entr,detr', entr(igout,k), detr(igout,k)
+            ENDDO
+            print *, 'fm-', fm(igout,1)
+            print *, '-------------------------------'
             abort_message = 'Negative incoming fm in thermcell_flux'
             CALL abort_physic(modname,abort_message,1)
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test entrainment positif
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : According to the previous derivations, we MUST have positive entrainment
 !      and detrainment everywhere!
 !      Indeed, they are set to max between zero and a computed value.
 !      Then tests are uncompromising.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         
-         labort_physic=.false.
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          DO ig=1,ngrid
@@ -212,16 +212,17 @@
             print *, 'lmin,lmax', lmin(igout), lmax(igout)
             print *, '-------------------------------'
-            print *, 'fm+', fm(igout,lout+1)
-            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
-            print *, 'fm ', fm(igout,lout)
-            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
-            print *, 'fm-', fm(igout,lout-1)
+            DO k=nlay,1,-1
+               print *, 'fm ', fm(igout,k+1)
+               print *, 'entr,detr', entr(igout,k), detr(igout,k)
+            ENDDO
+            print *, 'fm-', fm(igout,1)
+            print *, '-------------------------------'
             abort_message = 'Negative entrainment in thermcell_flux'
             CALL abort_physic(modname,abort_message,1)
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test detrainment positif
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
@@ -238,16 +239,17 @@
             print *, 'lmin,lmax', lmin(igout), lmax(igout)
             print *, '-------------------------------'
-            print *, 'fm+', fm(igout,lout+1)
-            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
-            print *, 'fm ', fm(igout,lout)
-            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
-            print *, 'fm-', fm(igout,lout-1)
+            DO k=nlay,1,-1
+               print *, 'fm ', fm(igout,k+1)
+               print *, 'entr,detr', entr(igout,k), detr(igout,k)
+            ENDDO
+            print *, 'fm-', fm(igout,1)
+            print *, '-------------------------------'
             abort_message = 'Negative detrainment in thermcell_flux'
             CALL abort_physic(modname,abort_message,1)
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test sur fraca constante ou croissante au-dessus de lalim
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
          
 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
@@ -267,7 +269,7 @@
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test sur flux de masse constant ou decroissant au-dessus de lalim
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
          
 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
@@ -283,9 +285,9 @@
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test detrainment inferieur au flux de masse
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : Even if fm has no negative value, it can be lesser than detr.
 !      That's not suitable because when we will mix the plume with the
@@ -294,5 +296,5 @@
 !      otherwise we get a negative outgoing mass flux (cf. above).
 !      That is why we correct the detrainment as follows.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          DO ig=1,ngrid
@@ -314,5 +316,5 @@
                ncorecfm6  = ncorecfm6 + 1
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le
 ! detrainement est plus fort que le flux de masse, on stope le thermique.
@@ -320,7 +322,7 @@
 !
 ! AB : Do we have to stop the plume here?
-! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are
-!      set to zero above this level!
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB : Attention, if lmax is modified, be sure that fm, zw2, entr, detr are set
+!      to zero above this level!
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               IF (l.gt.lalim(ig)) THEN
 !                  lmax(ig)=l
@@ -331,5 +333,5 @@
 !               ENDIF
             ENDIF
-!            
+            
 !            IF (l.gt.lmax(ig)) THEN
 !               detr(ig,l) = 0.
@@ -360,9 +362,9 @@
 !         ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test fraction masse entrainee inferieure a fomass_max
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : Entrainment is bigger than the maximal authorized value.
 !      If we consider that the excess entrainement is in fact plume air which
@@ -374,5 +376,5 @@
 !      detrainment and mass flux profiles such as we do not exceed the maximal
 !      authorized entrained mass.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
          labort_physic=.false.
@@ -432,5 +434,5 @@
             print *, '   ptimestep  :', ptimestep
             print *, 'norm  :', f(igout)
-            print *, 'alim* :', alim_star(igout,lout)
+!            print *, 'alim* :', alim_star(igout,lout)
             print *, 'entr* :', entr_star(igout,lout)
             print *, '--------------------'
@@ -447,9 +449,9 @@
          ENDIF
          
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Test fraction couverte inferieure a alphamax
-!------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! FH : Partie a revisiter.
 !      Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1
@@ -461,8 +463,9 @@
 !      semble de toutes facons deraisonable d'avoir des fractions de 1...
 !      Ci-dessous, et dans l'etat actuel, le premier des deux if est sans doute inutile.  
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
          
          DO ig=1,ngrid
             zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alphamax
+            
             IF (fm(ig,l+1) .gt. zfm) THEN
                f_old = fm(ig,l+1)
@@ -470,10 +473,10 @@
                detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : The previous change doesn't observe the equation df/dz = e - d.
 !      To avoid this issue, what is better to do? I choose to increase the
 !      entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and
 !      there are never any problems.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                IF (l.lt.lmax(ig)) THEN
                   entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
@@ -498,12 +501,12 @@
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! We check if fm, entr and detr vanish correctly in layer lmax
-!------------------------------------------------------------------------------
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+      
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
 !      Theoretically, we always get lmax >= lmin >= linf > 0
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       DO ig=1,ngrid
          IF (lmax(ig).gt.0) THEN
@@ -514,7 +517,7 @@
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Impression des bidouilles utilisees pour corriger les problemes
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       IF (prt_level.ge.1) THEN
@@ -553,7 +556,7 @@
       ENDIF
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : temporary test added to check the validity of eq. df/dz = e - d
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !      DO l=1,nlay
 !         DO ig=1,ngrid
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90	(revision 2127)
@@ -2,12 +2,11 @@
 !
 !
-      SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,          &
-      &                           zw2,zlev,lmax,zmax,zmix,                    &
-      &                           wmax,f_star,lev_out)
-      
-      
-!==============================================================================
-! thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
-!==============================================================================
+SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
+                            zlev,lmax,zmax,zmix,wmax,f_star)
+      
+      
+!===============================================================================
+!  Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
+!===============================================================================
       
       USE thermcell_mod
@@ -17,21 +16,19 @@
       
       
-!==============================================================================
+!===============================================================================
 ! Declaration
-!==============================================================================
-      
-!   inputs:
-!   -------
-      
-      INTEGER ngrid,nlay
-      INTEGER lalim(ngrid)
+!===============================================================================
+      
+!     Inputs:
+!     -------
+      
+      INTEGER ngrid, nlay
       INTEGER lmin(ngrid)
-      INTEGER lev_out                           ! niveau pour les print
       
       REAL zlev(ngrid,nlay+1)
       REAL f_star(ngrid,nlay+1)
       
-!   outputs:
-!   --------
+!     Outputs:
+!     --------
       
       INTEGER lmix(ngrid)
@@ -39,13 +36,13 @@
       
       REAL linter(ngrid)
-      REAL zmax(ngrid)
-      REAL zmix(ngrid)
-      REAL wmax(ngrid)
-      REAL zw2(ngrid,nlay+1)
-      
-!   local:
-!   ------
-      
-      INTEGER ig,l
+      REAL zmax(ngrid)                          ! maximal vertical speed
+      REAL zmix(ngrid)                          ! altitude of maximal vertical speed
+      REAL wmax(ngrid)                          ! maximal vertical speed
+      REAL zw2(ngrid,nlay+1)                    ! square vertical speed (becomes its squere root)
+      
+!     Local:
+!     ------
+      
+      INTEGER ig, l
       
       REAL num(ngrid)
@@ -54,10 +51,10 @@
       REAL zlev2(ngrid,nlay+1)
       
-!==============================================================================
+!===============================================================================
 ! Initialization
-!==============================================================================
-      
-      DO ig=1,ngrid
-         lmax(ig) = lalim(ig)
+!===============================================================================
+      
+      DO ig=1,ngrid
+         lmax(ig) = lmin(ig)
          lmix(ig) = lmin(ig)
       ENDDO
@@ -76,42 +73,41 @@
       ENDDO
       
-!==============================================================================
+!===============================================================================
 ! Calcul de la hauteur max du thermique
-!==============================================================================
-      
-!------------------------------------------------------------------------------
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
 ! Calcul de lmax
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO ig=1,ngrid
          DO l=nlay,lmin(ig)+1,-1
-            IF (zw2(ig,l).le.0.) THEN
+            IF (zw2(ig,l).le.0..or.f_star(ig,l).le.0.) THEN
                lmax(ig) = l - 1
-            ELSEIF (f_star(ig,l).le.0.) THEN
-               lmax(ig) = l - 1
             ENDIF
          ENDDO
       ENDDO
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! On traite le cas particulier qu'il faudrait eviter ou le thermique
 ! atteind le haut du modele ...
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
-      DO ig=1,ngrid
-         IF (zw2(ig,nlay).gt.1.e-10) THEN
-            print *, 'WARNING: a thermal plume reaches the highest level!'
-            lmax(ig)=nlay
-         ENDIF
-      ENDDO
-      
-!------------------------------------------------------------------------------
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO ig=1,ngrid
+         IF (zw2(ig,nlay).gt.0.) THEN
+            print *, 'WARNING: a thermal plume reaches the highest layer!'
+            print *, 'ig,l', ig, nlay
+            print *, 'zw2', zw2(ig,nlay)
+            lmax(ig) = nlay
+         ENDIF
+      ENDDO
+      
+!-------------------------------------------------------------------------------
 ! Calcul de zmax continu via le linter
-!------------------------------------------------------------------------------
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+      
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : lmin=lmax means the plume is not active and then zw2=0 at each level.
 !      Otherwise we have lmax>lmin.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       
       DO ig=1,ngrid
@@ -125,12 +121,12 @@
       ENDDO
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : zmax is equal to zlevinter-zlev(lmin) because lmin can be different
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB : zmax is now equal to zlevinter-zlev(lmin) because lmin can be different
 !      from 1.
-!      We have to substract this level because zmax must be the plume height
-!      (to derive the good closure), not the plume highest altitude.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
-      DO  ig=1,ngrid
+!      We have to substract this level because zmax must be the plume height (to
+!      derive the good closure), not the plume highest altitude.
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      
+      DO ig=1,ngrid
          zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig)   &
          &             + zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)      &
@@ -139,11 +135,11 @@
       ENDDO
       
-!==============================================================================
+!===============================================================================
 ! Calcul de wmax et du niveau d'inversion.
-!==============================================================================
-      
-!------------------------------------------------------------------------------
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
 ! Calcul de lmix et wmax
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO l=1,nlay
@@ -152,10 +148,12 @@
                IF (zw2(ig,l).lt.0.) THEN
                   print *, 'ERROR: zw2 has negative value(s)!'
+                  print *, 'ig,l', ig, l
+                  print *, 'zw2', zw2(ig,l)
                   CALL abort
                ENDIF
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : WARNING zw2 becomes its square root!
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                zw2(ig,l) = sqrt(zw2(ig,l))
                
@@ -170,12 +168,12 @@
       ENDDO
       
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
 ! Calcul de zmix continu (profil parabolique des vitesses)
-!------------------------------------------------------------------------------
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+      
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! AB : We substract zlev(lmin) in zmix formula because we have to
 !      compare it with zmax which is zlev(linter)-zlev(lmin).
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       
       DO ig=1,ngrid
@@ -197,5 +195,5 @@
             ELSE
                zmix(ig) = zlev(ig,lmix(ig)) - zlev(ig,lmin(ig))
-               print*,'WARNING: problematic zmix value!'
+               print *, 'WARNING: problematic zmix value!'
             ENDIF
          ELSE
@@ -204,27 +202,24 @@
       ENDDO
       
-!==============================================================================
+!===============================================================================
 ! Verification zmax > zmix
-!==============================================================================
+!===============================================================================
       
       DO ig=1,ngrid
          IF ((zmax(ig)-zmix(ig)).lt.0.) THEN
-            
             zmix(ig) = 0.9 * zmax(ig)
-            
-            IF (prt_level.ge.100) THEN
-               print *, 'WARNING: we have zmix > zmax!'
-            ENDIF
-         ENDIF
-      ENDDO
-      
-!------------------------------------------------------------------------------
+            print *, 'WARNING: we have zmix > zmax!'
+            print *, 'zmax,zmix', zmax(ig), zmix(ig)
+         ENDIF
+      ENDDO
+      
+!-------------------------------------------------------------------------------
 ! Calcul du nouveau lmix correspondant
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
       
       DO ig=1,ngrid
          DO l=1,nlay
             IF (zmix(ig).ge.zlev2(ig,l).and.zmix(ig).lt.zlev2(ig,l+1)) THEN
-               lmix(ig)=l
+               lmix(ig) = l
             ENDIF
          ENDDO
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2121)
+++ 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
-
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90	(revision 2127)
@@ -6,31 +6,18 @@
 ! Flags for computations
                                                                   !  default
-      INTEGER,PARAMETER :: iflag_thermals_optflux     = 1         !  0        !
-      INTEGER,PARAMETER :: iflag_thermals_closure     = 2         !  2        !
-      INTEGER,PARAMETER :: iflag_thermals             = 18        !  18       !
-      
-! Flags for (terrestrial) diagnoses
-      
-      LOGICAL,PARAMETER :: sorties                    = .false.   !  false
-      INTEGER,PARAMETER :: iflag_trig_bl              = 1         !  1
-      INTEGER,PARAMETER :: iflag_clos_bl              = 1         !  1
-      INTEGER,PARAMETER :: iflag_coupl                = 5         !  5
+      INTEGER,PARAMETER :: iflag_thermals_optflux     = 1         !  0        -
+      INTEGER,PARAMETER :: dqimpl                     = 1         !  1        flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme)
       
 ! Physical parameters
       
-      REAL,PARAMETER :: fact_thermals_ed_dz           = 0.007     !  0.007    !
       REAL,PARAMETER :: r_aspect_thermals             = 1.0       !           Aspect ratio of the thermals (width / height)
       REAL,PARAMETER :: tau_thermals                  = 0.        !  0.       Relaxation time
-      REAL,PARAMETER :: betalpha                      = 0.9       !  0.9      !
-      REAL,PARAMETER :: afact                         = 2./3.     !  2./3.    !
-      REAL,PARAMETER :: fact_epsilon                  = 0.000     !  0.002    !
-      REAL,PARAMETER :: detr_q_power                  = 0.5       !  0.5      !
-      REAL,PARAMETER :: detr_q_coef                   = 0.012     !  0.012    !
-      REAL,PARAMETER :: mix0                          = 0.        !  0.       !
-      REAL,PARAMETER :: detr_min                      = 1.d-5     !  1.e-5    Minimal detrainment value
-      REAL,PARAMETER :: entr_min                      = 1.d-5     !  1.e-5    Maximal detrainment value
-      REAL,PARAMETER :: alphamax                      = 0.7       !           Maximal permitted updraft fraction
-      REAL,PARAMETER :: fomass_max                    = 0.5       !           Maximal permitted outgoing layer mass fraction
-      REAL,PARAMETER :: pres_limit                    = 1.e5      !  1.e5     !
+      REAL,PARAMETER :: betalpha                      = 1.0       !  0.9      - between 0 (e=d) and 1 (rho*fraca=cst)
+      REAL,PARAMETER :: afact                         = 1.        !  2./3.    - buoyancy contribution, between 0 and 1
+      REAL,PARAMETER :: fact_epsilon                  = 1.e-4     !  2.e-3    - friction
+      REAL,PARAMETER :: nu                            = 0.000     !           Geometrical contributions to entrainment and detrainment
+      REAL,PARAMETER :: alphamax                      = 0.7       !  0.7      Maximal permitted updraft fraction
+      REAL,PARAMETER :: fomass_max                    = 0.5       !  0.5      Maximal permitted outgoing layer mass fraction
+      REAL,PARAMETER :: pres_limit                    = 2.e5      !  1.e5     -
       
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -43,5 +30,5 @@
 !      ngrid. 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      INTEGER,PARAMETER :: linf                       = 2         !     1
+      INTEGER,PARAMETER :: linf                       = 1         !     1
       
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -50,8 +37,4 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       REAL,PARAMETER :: d_temp                        = 0.        !     0.
-      
-! Parameters for diagnoses
-      
-      REAL,PARAMETER :: alp_bl_k                      = 0.5       !     0.5
       
 ! Physical constants
@@ -64,10 +47,4 @@
       
 !$OMP THREADPRIVATE(RTT, RG, RKAPPA, RPI, RD)
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : Parameters needed only for a loop in thermcell_alp (diagnoses).
-!      Maybe to be removed.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      INTEGER,PARAMETER :: nbsrf                      = 1
       
       
Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90	(revision 2121)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90	(revision 2127)
@@ -2,24 +2,27 @@
 !
 !
-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,ztva,ztla,zqla,zqta,zha,          &
-                           zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,          &
-                           lmin,lev_out,lunout1,igout)
-      
-      
-!==============================================================================
-! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
-! AB : ql means "liquid water mass mixing ratio"
-!      qt means "total water mass mixing ratio"
-!      TP means "potential temperature"
-!      TRPV means "virtual potential temperature with latent heat release"  
-!      TPV means "virtual potential temperature"
-!      TR means "temperature with latent heat release"
-!==============================================================================
+SUBROUTINE 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)
+      
+      
+!===============================================================================
+!  Purpose: calcule les valeurs de qt, thetal et w dans l ascendance
+!  
+!  Nota Bene
+!     ql   means "non-gaseous water mass mixing ratio" (liquid and solid)
+!     qv   means "vapor mass mixing ratio"
+!     qt   means "total water mass mixing ratio"
+!     TP   means "potential temperature"
+!     TRPV means "virtual potential temperature with latent heat release"  
+!     TPV  means "virtual potential temperature"
+!     TR   means "temperature with latent heat release"
+!===============================================================================
       
       USE print_control_mod, ONLY: prt_level
       USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
+      USE tracer_h, ONLY: igcm_h2o_vap
       USE thermcell_mod
       
@@ -27,28 +30,24 @@
       
       
-!==============================================================================
+!===============================================================================
 ! Declaration
-!==============================================================================
+!===============================================================================
       
 !     Inputs:
 !     -------
       
-      INTEGER itap
-      INTEGER ngrid
-      INTEGER nlay
-      INTEGER lunout1
-      INTEGER igout
-      INTEGER lev_out                           ! niveau pour les print
-      
-      REAL ptimestep                            ! time step
+      INTEGER ngrid, nlay, nq
+      
+      REAL ptimestep
+      REAL rhobarz(ngrid,nlay)                  ! Levels density
+      REAL zlev(ngrid,nlay+1)                   ! Levels altitude
+      REAL pplev(ngrid,nlay+1)                  ! Levels pressure
+      REAL pphi(ngrid,nlay)                     ! Geopotential
+      REAL zpopsk(ngrid,nlay)                   ! Exner function
+      
       REAL ztv(ngrid,nlay)                      ! TRPV environment
-      REAL zthl(ngrid,nlay)                     ! TP   environment
-      REAL po(ngrid,nlay)                       ! qt   environment
-      REAL zl(ngrid,nlay)                       ! ql   environment
-      REAL rhobarz(ngrid,nlay)                  ! levels density
-      REAL zlev(ngrid,nlay+1)                   ! levels altitude
-      REAL pplev(ngrid,nlay+1)                  ! levels pressure
-      REAL pphi(ngrid,nlay)                     ! geopotential
-      REAL zpspsk(ngrid,nlay)                   ! Exner function
+      REAL zhl(ngrid,nlay)                      ! TP   environment
+      REAL zqt(ngrid,nlay)                      ! qt   environment
+      REAL zql(ngrid,nlay)                      ! ql   environment
       
 !     Outputs:
@@ -56,12 +55,5 @@
       
       INTEGER lmin(ngrid)                       ! plume base level (first unstable level)
-      INTEGER lalim(ngrid)                      ! higher alimentation level
       INTEGER lmix(ngrid)                       ! maximum vertical speed level
-      INTEGER lmix_bis(ngrid)                   ! maximum vertical speed level (modified)
-      
-      REAL alim_star(ngrid,nlay)                ! normalized alimentation
-      REAL alim_star_tot(ngrid)                 ! integrated alimentation
-      
-      REAL f0(ngrid)                            ! previous time step mass flux norm
       
       REAL detr_star(ngrid,nlay)                ! normalized detrainment
@@ -70,14 +62,9 @@
       
       REAL ztva(ngrid,nlay)                     ! TRPV plume (after mixing)
-      REAL ztva_est(ngrid,nlay)                 ! TRPV plume (before mixing)
-      REAL ztla(ngrid,nlay)                     ! TP   plume
+      REAL zhla(ngrid,nlay)                     ! TP   plume ?
       REAL zqla(ngrid,nlay)                     ! ql   plume (after mixing)
-      REAL zqta(ngrid,nlay)                     ! qt   plume
-      REAL zha(ngrid,nlay)                      ! TRPV plume
-      
-      REAL w_est(ngrid,nlay+1)                  ! updraft square vertical speed (before mixing)
-      REAL zw2(ngrid,nlay+1)                    ! updraft square vertical speed (after mixing)
-      
-      REAL zqsatth(ngrid,nlay)                  ! saturation vapor pressure (after mixing)
+      REAL zqta(ngrid,nlay)                     ! qt   plume ?
+      REAL zqsa(ngrid,nlay)                     ! qsat plume (after mixing)
+      REAL zw2(ngrid,nlay+1)                    ! w    plume (after mixing)
       
 !     Local:
@@ -85,101 +72,83 @@
       
       INTEGER ig, l, k
-      INTEGER lt
-      INTEGER lm
-      
+      
+      REAL ztva_est(ngrid,nlay)                 ! TRPV plume (before mixing)
       REAL zqla_est(ngrid,nlay)                 ! ql   plume (before mixing)
       REAL zta_est(ngrid,nlay)                  ! TR   plume (before mixing)
-      REAL zbuoy(ngrid,nlay)                    ! plume buoyancy
-      REAL zbuoyjam(ngrid,nlay)                 ! plume buoyancy (modified)
-      
-      REAL ztemp(ngrid)                         ! temperature for saturation vapor pressure computation in plume
-      REAL zqsat(ngrid)                         ! saturation vapor pressure (before mixing)
-      REAL zdz                                  ! layers height
-      REAL ztv2(ngrid,nlay)                     ! ztv + d_temp * Dirac(l=lmin)
-      
-      REAL zalpha                               !
-      REAL zdqt(ngrid,nlay)                     !
+      REAL zqsa_est(ngrid)                      ! qsat plume (before mixing)
+      REAL zw2_est(ngrid,nlay+1)                ! w    plume (before mixing)
+      
+      REAL zta(ngrid,nlay)                      ! TR   plume (after mixing)
+      
+      REAL zbuoy(ngrid,nlay)                    ! Plume buoyancy
+      REAL ztemp(ngrid)                         ! Temperature for saturation vapor pressure computation in plume
+      REAL zdz                                  ! Layers heights
+      REAL ztv2(ngrid,nlay)                     ! ztv + d_temp * Dirac(l=linf)
+      
       REAL zbetalpha                            !
       REAL zdw2                                 ! 
       REAL zdw2bis                              ! 
       REAL zw2fact                              ! 
-      REAL zw2factbis                           ! 
-      REAL zw2m                                 ! 
-      REAL zdzbis                               !
-      REAL coefzlmel                            !
-      REAL zdz2                                 !
-      REAL zdz3                                 !
-      REAL lmel                                 !
-      REAL zlmel                                !
-      REAL zlmelup                              !
-      REAL zlmeldwn                             !
-      REAL zlt                                  !
-      REAL zltdwn                               !
-      REAL zltup                                ! useless here
-      
-      REAL psat                                 ! dummy argument for Psat_water()
-      
-      LOGICAL active(ngrid)                     ! if the plume is active at ig,l (speed and incoming mass flux > 0 or l=lmin)
-      LOGICAL activetmp(ngrid)                  ! if the plus is active at ig,l (active=true and outgoing mass flux > 0)
-      LOGICAL, SAVE :: first = .true.           ! if it is the first time step
-      
-!$OMP THREADPRIVATE(first)
-      
-!==============================================================================
+      REAL zw2m                                 ! Average vertical velocity between two successive levels
+      REAL gamma                                ! Plume acceleration term (to compute vertical velocity)
+      REAL test                                 ! Test to know how to compute entrainment and detrainment
+       
+      REAL psat                                 ! Dummy argument for Psat_water()
+      
+      LOGICAL active(ngrid)                     ! If the plume is active at ig (speed and incoming mass flux > 0 or l=lmin)
+      LOGICAL activetmp(ngrid)                  ! If the plume is active at ig (active=true and outgoing mass flux > 0)
+      
+!===============================================================================
 ! Initialization
-!==============================================================================
+!===============================================================================
       
       zbetalpha = betalpha / (1. + betalpha)
       
-      ztva(:,:)        = ztv(:,:)               ! ztva     is set to the virtual potential temperature without latent heat release
-      ztva_est(:,:)    = ztva(:,:)              ! ztva_est is set to the virtual potential temperature without latent heat release
-      ztla(:,:)        = zthl(:,:)              ! ztla     is set to the potential temperature
-      zqta(:,:)        = po(:,:)                ! zqta     is set to qt
-      zqla(:,:)        = 0.                     ! zqla     is set to ql
-      zqla_est(:,:)    = 0.                     ! zqla_est is set to ql
-      zha(:,:)         = ztva(:,:)              ! zha      is set to the plume virtual potential temperature without latent heat release
-      
-      zqsat(:)         = 0.
-      zqsatth(:,:)     = 0.
-      
-      w_est(:,:)       = 0.
+      ztva(:,:)        = ztv(:,:)                                                ! ztva     is set to TPV environment
+      zhla(:,:)        = zhl(:,:)                                                ! zhla     is set to TP  environment
+      zqta(:,:)        = zqt(:,:)                                                ! zqta     is set to qt  environment
+      zqla(:,:)        = zql(:,:)                                                ! zqla     is set to ql  environment
+      
+      zqsa_est(:)      = 0.
+      zqsa(:,:)        = 0.
+      
+      zw2_est(:,:)     = 0.
       zw2(:,:)         = 0.
       
       zbuoy(:,:)       = 0.
-      zbuoyjam(:,:)    = 0.
       
       f_star(:,:)      = 0.
       detr_star(:,:)   = 0.
       entr_star(:,:)   = 0.
-      alim_star(:,:)   = 0.
-      alim_star_tot(:) = 0.
       
       lmix(:)          = 1
-      lmix_bis(:)      = 2
-      lalim(:)         = 1
-      lmin(:)          = linf
+      lmin(:)          = 1
       
       ztv2(:,:)        = ztv(:,:)
       ztv2(:,linf)     = ztv(:,linf) + d_temp
       
-!==============================================================================
-! 0. Calcul de l'alimentation
-!==============================================================================
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : Convective plumes can go off from every layer above the linf-th and
-!      where pressure is lesser than pres_limit (cf. thermcell_plume).
-!      The second constraint is added to avoid the parametrization occurs too
-!      high when the low atmosphere is stable.
-!      However, once there is a triggered plume, it can rise as high as its
-!      velocity allows it (it can overshoot).
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!===============================================================================
+! First layer computation
+!===============================================================================
+      
       DO ig=1,ngrid
          active(ig) = .false.
          l = linf
-         DO WHILE ((.not.active(ig)) .and. pplev(ig,l+1).gt.pres_limit .and. l.lt.nlay)
-            IF (ztv2(ig,l).gt.ztv2(ig,l+1)) THEN
+         DO WHILE (.not.active(ig).and.(pplev(ig,l+1).GT.pres_limit).and.(l.LT.nlay))
+            zdz = (zlev(ig,l+1) - zlev(ig,l))
+            zw2(ig,l+1) = 2. * afact * RG * zdz                               &  ! Do we have to divide by 1+betalpha (consider entrainment) ?
+            &           * (ztv2(ig,l) - ztv2(ig,l+1))                         &
+            &           / (ztv2(ig,l+1) * (1. + betalpha))
+            zw2m = zw2(ig,l+1) / 2.
+            entr_star(ig,l) = zdz * zbetalpha * (afact * RG * (ztv2(ig,l)     &
+            &               - ztv2(ig,l+1)) / ztv2(ig,l+1) / zw2m - fact_epsilon)
+            IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(entr_star(ig,l).GT.0.)) THEN
                active(ig) = .true.
                lmin(ig) = l
+               f_star(ig,l+1) = entr_star(ig,l)
+               zw2_est(ig,l+1) = zw2(ig,l+1)
+            ELSE
+               zw2(ig,l+1) = 0.
+               entr_star(ig,l) = 0.
             ENDIF
             l = l + 1
@@ -187,46 +156,19 @@
       ENDDO
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : On pourrait n'appeler thermcell_alim que si la plume est active
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      CALL thermcell_alim(ngrid,nlay,ztv2,zlev,alim_star,lalim,lmin)
-      
-!==============================================================================
-! 1. Calcul dans la premiere couche
-!==============================================================================
-      
-      DO ig=1,ngrid
-         IF (active(ig)) THEN
-            
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : plume takes the environment features for every layer below lmin.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            DO l=1,lmin(ig)
-               ztla(ig,l) = zthl(ig,l)
-               zqta(ig,l) = po(ig,l)
-               zqla(ig,l) = zl(ig,l)
-            ENDDO
-            
-            l = lmin(ig)
-            f_star(ig,l+1) = alim_star(ig,l)
-            
-            zw2(ig,l+1) = 2. * RG * (zlev(ig,l+1) - zlev(ig,l))               &
-            &                * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
-            
-            w_est(ig,l+1) = zw2(ig,l+1)
-         ENDIF
-      ENDDO
-      
-!==============================================================================
-! 2. Boucle de calcul de la vitesse verticale dans le thermique
-!==============================================================================
+!===============================================================================
+! Thermal plumes computations
+!===============================================================================
       
       DO l=2,nlay-1
          
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : we decide here if the plume is still active or not. When the plume's
-!      first level is reached, we set active to "true". Otherwise, it is given
-!      by zw2, f_star, alim_star and entr_star.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!-------------------------------------------------------------------------------
+! Check if thermal plume is (still) active
+!-------------------------------------------------------------------------------
+         
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB: we decide here if the plume is still active or not. When the plume's
+!     first level is reached, we set active to "true". Otherwise, it is given
+!     by zw2, f_star and entr_star.
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          DO ig=1,ngrid
             IF (l==lmin(ig)+1) THEN
@@ -235,90 +177,51 @@
             
             active(ig) = active(ig)                                           &
-            &            .and. zw2(ig,l)>1.e-10                               &
-            &            .and. f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)>1.e-10
-         ENDDO
-         
-         ztemp(:) = zpspsk(:,l) * ztla(:,l-1)
-         
-         DO ig=1,ngrid
-            CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsat(ig))
-         ENDDO
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : we compute thermodynamical values and speed in the plume in the layer l
-!      without mixing with environment.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            &          .and. (zw2(ig,l).GT.1.e-10)                            &
+            &          .and. (f_star(ig,l) + entr_star(ig,l)).GT.1.e-10
+         ENDDO
+         
+         ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
+         
+         DO ig=1,ngrid
+            CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
+         ENDDO
+         
+!-------------------------------------------------------------------------------
+! Vertical speed (before mixing between plume and environment)
+!-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
             IF (active(ig)) THEN
-               zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig))                   ! zqla_est set to ql   plume
-               ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)   ! ztva_est set to TR   plume
-               zta_est(ig,l)  = ztva_est(ig,l)                                   ! zta_est  set to TR   plume
-               ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)                      ! ztva_est set to TRP  plume
-               ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)         &  ! ztva_est set to TRPV plume
-               &              - zqla_est(ig,l))-zqla_est(ig,l))
-               
-               zbuoy(ig,l) = RG * (ztva_est(ig,l)-ztv(ig,l)) / ztv(ig,l)
+               zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsa_est(ig))                ! zqla_est set to ql   plume
+               zta_est(ig,l)  = zhla(ig,l-1) * zpopsk(ig,l)                   &  ! zta_est  set to TR   plume
+               &              + RLvCp * zqla_est(ig,l)
+               ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l)                  &  ! ztva_est set to TRPV plume
+               &              * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
+               
+               zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
                zdz = zlev(ig,l+1) - zlev(ig,l)
                
-!==============================================================================
-! 3. Calcul de la flotabilite modifiee par melange avec l'air au dessus
-!==============================================================================
-               
-               lmel = fact_thermals_ed_dz * zlev(ig,l)
-               zlmel = zlev(ig,l) + lmel
-               lt = l + 1
-               zlt = zlev(ig,lt)
-               zdz2 = zlev(ig,lt) - zlev(ig,l)
-               
-               DO while (lmel.gt.zdz2)
-                  lt   = lt + 1
-                  zlt  = zlev(ig,lt)
-                  zdz2 = zlev(ig,lt) - zlev(ig,l)
-               ENDDO
-               
-!               IF (lt-l.gt.1) THEN
-!                  print *, 'WARNING: lt is greater than l+1!'
-!                  print *, 'l,lt', l, lt
-!               ENDIF
-               
-               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)
-               
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : initial formulae
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB: initial formulae
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
 !               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
 !               zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
-!               w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
-!               w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : own derivation for w_est (Rio 2010 formula with e=d=0)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!               zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2)
+!               zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB: own derivation for zw2_est (Rio et al. 2010)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               zw2fact = 2. * fact_epsilon * zdz
 !               zdw2 = 2. * afact * zbuoy(ig,l) * zdz
                zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
                zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
-               w_est(ig,l+1) = Max(0., exp(-zw2fact) * w_est(ig,l) + zdw2)
-               
-!               IF (w_est(ig,l+1).le.0.) THEN
-!                  print *, 'WARNING: w_est is negative!'
-!                  print *, 'l,w_est', l+1, w_est(ig,l+1)
-!               ENDIF
-            ENDIF
-         ENDDO
-         
-!==============================================================================
-! 4. Calcul de l'entrainement et du detrainement
-!==============================================================================
+               zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2)
+            ENDIF
+         ENDDO
+         
+!-------------------------------------------------------------------------------
+! Mass flux, entrainment and detrainment
+!-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
@@ -326,44 +229,23 @@
                
                zdz = zlev(ig,l+1) - zlev(ig,l)
-               
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : The next test is added to avoid divisions by zero when w_est vanishes.
-!      Indeed, entr and detr computed here are of no importance because w_est
-!      <= 0 means it will be the last layer reached by the plume and then they
-!      will be reset in thermcell_flux.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               IF (w_est(ig,l+1).eq.0.) THEN
-                  zw2m = 1.
-                  zalpha = 0.
+               zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.                     ! AB: est-ce la bonne vitesse a utiliser ?
+               gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
+               
+               IF (zw2_est(ig,l).GT.0.) THEN
+                  test = gamma / zw2m - nu
                ELSE
-                  zw2m = w_est(ig,l+1)
-                  zalpha = f0(ig) * f_star(ig,l) / sqrt(w_est(ig,l+1)) / rhobarz(ig,l)
+                  print *, 'ERROR: zw2_est is negative while plume is active!'
+                  print *, 'ig,l', ig, l
+                  print *, 'zw2_est', zw2_est(ig,l)
+                  call abort
                ENDIF
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : The next test is added to avoid a division by zero if there is no water
-!      in the environment.
-!      In the case where there is no water in the env. but water in the plume
-!      (ascending from depth) we set the effect on detrainment equal to zero
-!      but at the next time step, po will be positive thanks to the mixing and
-!      then the physical effect of the water gradient will be taken on board.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-               IF (po(ig,l).lt.1.e-6) THEN
-!                  print *, 'WARNING: po=0 in layer',l,'!'
-!                  print *, 'po,zqta', po(ig,l), zqta(ig,l-1)
-                  zdqt(ig,l) = 0.0
+               IF (test.gt.0.) THEN
+                  detr_star(ig,l) = zdz * nu
+                  entr_star(ig,l) = zdz * (zbetalpha * gamma / zw2m + nu)
                ELSE
-                  zdqt(ig,l) = max(zqta(ig,l-1)-po(ig,l),0.) / po(ig,l)
+                  detr_star(ig,l) = zdz * (nu - betalpha * gamma / zw2m)
+                  entr_star(ig,l) = zdz * nu
                ENDIF
-               
-!------------------------------------------------------------------------------
-! Detrainment
-!------------------------------------------------------------------------------
-               
-               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 (detr_star(ig,l).lt.0.) THEN
@@ -372,14 +254,4 @@
 !               ENDIF
                
-!------------------------------------------------------------------------------
-! Entrainment
-!------------------------------------------------------------------------------
-               
-               entr_star(ig,l) = f_star(ig,l) * zdz * (                       &
-               &                 mix0 * 0.1 / (zalpha+0.001)                  &
-               &               + MAX(entr_min,                                &
-               &                 zbetalpha * afact * zbuoyjam(ig,l) / zw2m    &
-               &               - zbetalpha * fact_epsilon) )
-               
 !               IF (entr_star(ig,l).lt.0.) THEN
 !                  print *, 'WARNING: entrainment is negative!'
@@ -387,19 +259,5 @@
 !               ENDIF
                
-!------------------------------------------------------------------------------
-! Alimentation and entrainment
-!------------------------------------------------------------------------------
-               
-               IF (l.lt.lalim(ig)) THEN
-                  alim_star(ig,l) = max(alim_star(ig,l),entr_star(ig,l))
-                  entr_star(ig,l) = 0.
-               ENDIF
-               
-!------------------------------------------------------------------------------
-! Mass flux
-!------------------------------------------------------------------------------
-               
-               f_star(ig,l+1) = f_star(ig,l) + alim_star(ig,l)                &
-               &              + entr_star(ig,l) - detr_star(ig,l)
+               f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
                
 !               IF (f_star(ig,l+1).le.0.) THEN
@@ -411,58 +269,53 @@
          ENDDO
          
-!==============================================================================
-! 5. Calcul de la vitesse verticale en melangeant Tl et qt du thermique
-!==============================================================================
-         
-         activetmp(:) = active(:) .and. f_star(:,l+1)>1.e-10
-         
-!------------------------------------------------------------------------------
-! Calcul du melange avec l'environnement
-!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+! Mixing between thermal plume and environment
+!-------------------------------------------------------------------------------
+         
+         activetmp(:) = active(:) .and. (f_star(:,l+1).GT.1.e-10)
          
          DO ig=1,ngrid
             IF (activetmp(ig)) THEN
-               ztla(ig,l) = (f_star(ig,l) * ztla(ig,l-1)                      &  ! ztla is set to TP in plume (mixed)
-               &          + (alim_star(ig,l) + entr_star(ig,l)) * zthl(ig,l)) &
+               zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1)                      &  ! zhla is set to TP in plume (mixed)
+               &          + entr_star(ig,l) * zhl(ig,l))                      &
                &          / (f_star(ig,l+1) + detr_star(ig,l))
                zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) +                    &  ! zqta is set to qt in plume (mixed)
-               &          + (alim_star(ig,l) + entr_star(ig,l)) * po(ig,l))   &
+               &          + entr_star(ig,l) * zqt(ig,l))                      &
                &          / (f_star(ig,l+1) + detr_star(ig,l))
             ENDIF
          ENDDO
          
-         ztemp(:) = zpspsk(:,l) * ztla(:,l)
+         ztemp(:) = zpopsk(:,l) * zhla(:,l)
          
          DO ig=1,ngrid
             IF (activetmp(ig)) THEN
-               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsatth(ig,l))
-            ENDIF
-         ENDDO
-         
-!------------------------------------------------------------------------------
-! Calcul de la vitesse verticale zw2 apres le melange
-!------------------------------------------------------------------------------
+               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
+            ENDIF
+         ENDDO
+         
+!-------------------------------------------------------------------------------
+! Vertical speed (after mixing between plume and environment)
+!-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
             IF (activetmp(ig)) THEN
-               zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l))                     ! zqla is set to ql   plume (mixed)
-               ztva(ig,l) = ztla(ig,l) * zpspsk(ig,l)+RLvCp*zqla(ig,l)           ! ztva is set to TR   plume (mixed)
-               ztva(ig,l) = ztva(ig,l) / zpspsk(ig,l)                            ! ztva is set to TRP  plume (mixed)
-               zha(ig,l)  = ztva(ig,l)                                           ! zha  is set to TRP  plume (mixed)
-               ztva(ig,l) = ztva(ig,l) * (1. + RETV*(zqta(ig,l)-zqla(ig,l))   &  ! ztva is set to TRPV plume (mixed)
-               &                             - zqla(ig,l))
+               zqla(ig,l) = max(0.,zqta(ig,l)-zqsa(ig,l))                        ! zqla is set to ql   plume (mixed)
+               
+               zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) + RLvCp * zqla(ig,l)        ! ztva is set to TR   plume (mixed)
+               ztva(ig,l) = zta(ig,l) / zpopsk(ig,l)                          &  ! ztva is set to TRPV plume (mixed)           
+               &          * (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)
                
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : initial formula
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB: initial formula
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
 !               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
 !               zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : own derivation for zw2 (Rio 2010 formula)
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! AB: own derivation for zw2 (Rio et al. 2010)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 !               zw2fact = 2. * (fact_epsilon * zdz  + entr_star(ig,l) / f_star(ig,l))
 !               zdw2 = 2. * afact * zbuoy(ig,l) * zdz
@@ -470,26 +323,7 @@
                zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
                zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)
-               
-!               IF (zw2(ig,l+1).le.0.) THEN
-!                  print *, 'WARNING: zw2 is negative!'
-!                  print *, 'l,zw2', l+1, zw2(ig,l+1)
-!               ENDIF
-            ENDIF
-         ENDDO
-         
-      ENDDO
-      
-!==============================================================================
-! 6. New computation of 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
+            ENDIF
+         ENDDO
+         
       ENDDO
       
