Index: trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2144)
+++ trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2177)
@@ -42,7 +42,5 @@
 !    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
+!    Mix every tracer
 ! 
 !===============================================================================
@@ -61,30 +59,36 @@
 !     -------
       
-      INTEGER ngrid, nlay, nq
-      
-      REAL ptimestep
-      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
+      INTEGER, INTENT(in) :: ngrid
+      INTEGER, INTENT(in) :: nlay
+      INTEGER, INTENT(in) :: nq
+      
+      REAL, INTENT(in) :: ptimestep
+      REAL, INTENT(in) :: pplay(ngrid,nlay)           ! Layer pressure
+      REAL, INTENT(in) :: pplev(ngrid,nlay+1)         ! Level pressure
+      REAL, INTENT(in) :: pphi(ngrid,nlay)            ! Geopotential
+      
+      REAL, INTENT(in) :: pu(ngrid,nlay)              ! Zonal wind
+      REAL, INTENT(in) :: pv(ngrid,nlay)              ! Meridional wind
+      REAL, INTENT(in) :: pt(ngrid,nlay)              ! Temperature
+      REAL, INTENT(in) :: pq(ngrid,nlay,nq)           ! Tracers mass mixing ratio
+      
+      LOGICAL, INTENT(in) :: firstcall
       
 !     Outputs:
 !     --------
       
-      REAL pduadj(ngrid,nlay)                   ! u convective variations
-      REAL pdvadj(ngrid,nlay)                   ! v 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)
+      INTEGER, INTENT(out) :: lmax(ngrid)             ! Highest layer reached by the plume
+      INTEGER, INTENT(out) :: lmix(ngrid)             ! Layer in which plume vertical speed is maximal
+      INTEGER, INTENT(out) :: lmin(ngrid)             ! First unstable layer
+      
+      REAL, INTENT(out) :: pduadj(ngrid,nlay)         ! u convective variations
+      REAL, INTENT(out) :: pdvadj(ngrid,nlay)         ! v convective variations
+      REAL, INTENT(out) :: pdtadj(ngrid,nlay)         ! t convective variations
+      REAL, INTENT(out) :: pdqadj(ngrid,nlay,nq)      ! q convective variations
+      
+      REAL, INTENT(inout) :: f0(ngrid)                ! mass flux norm (after possible time relaxation)
+      REAL, INTENT(inout) :: fm0(ngrid,nlay+1)        ! mass flux      (after possible time relaxation)
+      REAL, INTENT(inout) :: entr0(ngrid,nlay)        ! entrainment    (after possible time relaxation)
+      REAL, INTENT(inout) :: detr0(ngrid,nlay)        ! detrainment    (after possible time relaxation)
       
 !     Local:
@@ -92,55 +96,50 @@
       
       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 zmix(ngrid)                          ! Altitude of maximal vertical speed
-      REAL zmax(ngrid)                          ! Maximal altitudes where plumes are active
-      REAL zmin(ngrid)                          ! Minimal altitudes where plumes are active
-      
-      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 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 f_star(ngrid,nlay+1)                 ! Normalized mass flux
-      REAL entr_star(ngrid,nlay)                ! Normalized entrainment (E* = e* dz)
-      REAL detr_star(ngrid,nlay)                ! Normalized detrainment (D* = d* dz)
-      
-      REAL fm(ngrid,nlay+1)                     ! Mass flux
-      REAL entr(ngrid,nlay)                     ! Entrainment (E = e dz)
-      REAL detr(ngrid,nlay)                     ! Detrainment (D = d dz)
-      
-      REAL f(ngrid)                             ! Mass flux norm
-      REAL lambda                               ! Time relaxation coefficent
-      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
-      REAL linter(ngrid)                        ! Level (continuous) of maximal vertical speed
-      REAL wmax(ngrid)                          ! Maximal vertical speed
-      REAL zw2(ngrid,nlay+1)                    ! Plume vertical speed
-      REAL zdthladj(ngrid,nlay)                 ! Potential temperature variations
-      REAL dummy(ngrid,nlay)                    ! Dummy argument for thermcell_dq()
+      
+      REAL zmix(ngrid)                                ! Altitude of maximal vertical speed
+      REAL zmax(ngrid)                                ! Maximal altitudes where plumes are active
+      REAL zmin(ngrid)                                ! Minimal altitudes where plumes are active
+      
+      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 zqla(ngrid,nlay)                           ! qv   plume
+      REAL zqta(ngrid,nlay)                           ! qt   plume
+      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 f_star(ngrid,nlay+1)                       ! Normalized mass flux
+      REAL entr_star(ngrid,nlay)                      ! Normalized entrainment (E* = e* dz)
+      REAL detr_star(ngrid,nlay)                      ! Normalized detrainment (D* = d* dz)
+      
+      REAL fm(ngrid,nlay+1)                           ! Mass flux
+      REAL entr(ngrid,nlay)                           ! Entrainment (E = e dz)
+      REAL detr(ngrid,nlay)                           ! Detrainment (D = d dz)
+      
+      REAL f(ngrid)                                   ! Mass flux norm
+      REAL lambda                                     ! Time relaxation coefficent
+      REAL fraca(ngrid,nlay+1)                        ! Updraft fraction
+      REAL wmax(ngrid)                                ! Maximal vertical speed
+      REAL zw2(ngrid,nlay+1)                          ! Plume vertical speed
+      REAL zdthladj(ngrid,nlay)                       ! Potential temperature variations
+      REAL dummy(ngrid,nlay)                          ! Dummy argument for thermcell_dq()
       
 !===============================================================================
@@ -154,17 +153,8 @@
       ENDIF
       
-      f_star(:,:) = 0.
-      entr_star(:,:) = 0.
-      detr_star(:,:) = 0.
-      
-      f(:) = 0.
-      
-      fm(:,:) = 0.
-      entr(:,:) = 0.
-      detr(:,:) = 0.
-      
-      lmax(:) = 1
-      lmix(:) = 1
-      lmin(:) = 1
+      DO ig=1,ngrid
+! AB: Minimal f0 value is set to 0. (instead of 1.e-2 in Earth version)
+         f0(ig) = MAX(f0(ig), 0.)
+      ENDDO
       
       pduadj(:,:) = 0.0
@@ -173,10 +163,5 @@
       pdqadj(:,:,:) = 0.0
       
-      DO ig=1,ngrid
-! AB: Careful: Hard-coded value from Earth version!
-!         f0(ig) = max(f0(ig), 1.e-2)
-! AB: No pescribed minimal value for f0
-         f0(ig) = max(f0(ig), 0.)
-      ENDDO
+      zdthladj(:,:) = 0.0
       
 !===============================================================================
@@ -212,11 +197,8 @@
       rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
       
+      rhobarz(:,1) = rho(:,1)
       IF (prt_level.ge.10) THEN
          print *, 'WARNING: density in the first layer is equal to density at the first level!'
-         print *, 'rhobarz(:,1)', rhobarz(:,1)
-         print *, 'rho(:,1)    ', rho(:,1)
-      ENDIF
-      
-      rhobarz(:,1) = rho(:,1)
+      ENDIF
       
       DO l=2,nlay
@@ -244,32 +226,32 @@
 !                                
 !     ---------------------------
-!                                        _
-!     ----- F_lmax+1=0 ------zmax         \
-!     lmax                                 |
-!     ------F_lmax>0-------------          |
-!                                          |
-!     ---------------------------          |
-!                                          |
-!     ---------------------------          |
-!                                          |
-!     ------------------wmax,zmix          |
-!     lmix                                 |
-!     ---------------------------          |
-!                                          |
-!     ---------------------------          |
-!                                          | E, D
-!     ---------------------------          |
-!                                          |
+!                                       _
+!     ----- 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 ------------        _/
+!         zt, zu, zv, zo, rho             |
+!     ---------------------------         |
+!                                         |
+!     ---------------------------         |
+!                                         |
+!     ---------------------------         |
+!                                         |
+!     ------F_lmin+1>0-----------         |
+!     lmin                                |
+!     ----- F_lmin=0 ------------       _/
 !                                
 !     ---------------------------
@@ -309,9 +291,9 @@
 !-------------------------------------------------------------------------------
       
-      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
-      &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
+      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,                           &
+      &                    ztv,zhl,zqt,zql,zlev,pplev,zpopsk,                 &
       &                    detr_star,entr_star,f_star,                        &
-      &                    ztva,zhla,zqla,zqta,zta,zqsa,                      &
-      &                    zw2,lmix,lmin)
+      &                    ztva,zhla,zqta,zqla,zqsa,                          &
+      &                    zw2,lmin)
       
 !-------------------------------------------------------------------------------
@@ -320,6 +302,6 @@
       
 ! AB: Careful, zw2 became its square root in thermcell_height!
-      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
-      &                     zlev,zmin,zmix,zmax,wmax,f_star)
+      CALL thermcell_height(ngrid,nlay,lmin,lmix,lmax,zlev,                   &
+      &                     zmin,zmix,zmax,zw2,wmax,f_star)
       
 !===============================================================================
@@ -341,7 +323,5 @@
       ENDIF
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! Test valable seulement en 1D mais pas genant
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! FH: Test valable seulement en 1D mais pas genant
       IF (.not. (f0(1).ge.0.) ) THEN
          print *, 'ERROR: mass flux norm is not positive!'
@@ -403,5 +383,5 @@
       
       CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
-      &                 zhl,zdthladj,dummy,lmin)
+      &                 zhl,zdthladj,dummy,lmin,lmax)
       
       DO l=1,nlay
@@ -417,5 +397,5 @@
       DO iq=1,nq
          CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,        &
-         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
+         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin,lmax)
       ENDDO
       
@@ -429,7 +409,7 @@
       ELSE
          CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
-         &                 zu,pduadj,zua,lmin)
+         &                 zu,pduadj,zua,lmin,lmax)
          CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
-         &                 zv,pdvadj,zva,lmin)
+         &                 zv,pdvadj,zva,lmin,lmax)
       ENDIF
       
