Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 2142)
+++ /trunk/LMDZ.GENERIC/README	(revision 2143)
@@ -1483,2 +1483,7 @@
  but would require to be able to have writediag in 5D)
  EDIT (22/05/19) : To have correct calculations, output is now exp(-tau), you need to postproc it to have tau.
+
+== 20/06/2019 == AB
+- update the thermal plume model (check formulae consistency between thermcell_plume and thercell_closure,
+  compute correctly thermal plume height, fix alimentation computation in the first unstable layer)
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90	(revision 2143)
@@ -3,5 +3,6 @@
 !
 SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
-                             lmax,alim_star,zmax,wmax,f)
+                             lmax,alim_star,zmin,zmax,wmax,f)
+      
       
 !===============================================================================
@@ -13,5 +14,5 @@
 ! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec
 ! 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
+! 3. Vectorisation en mettant les boucles en l a l'exterieur avec des if
 !===============================================================================
       
@@ -35,4 +36,5 @@
       REAL zlev(ngrid,nlay)
       REAL alim_star(ngrid,nlay)
+      REAL zmin(ngrid)
       REAL zmax(ngrid)
       REAL wmax(ngrid)
@@ -51,5 +53,5 @@
       REAL alim_star_tot(ngrid)
       REAL alim_star2(ngrid)
-      REAL zdenom(ngrid)
+      REAL plume_height(ngrid)
       
 !===============================================================================
@@ -64,4 +66,8 @@
       llmax = 1
       
+      DO ig=1,ngrid
+         plume_height(ig) = zmax(ig) - zmin(ig)
+      ENDDO
+      
 !===============================================================================
 ! Closure
@@ -73,5 +79,5 @@
       
       DO ig=1,ngrid
-         IF (lmax(ig).GT.llmax) THEN
+         IF (lmax(ig) > llmax) THEN
             llmax = lmax(ig)
          ENDIF
@@ -86,5 +92,5 @@
             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)))
+               &              / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l)))        ! => intergration is ok because alim_star = a* dz
                alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
             ENDIF
@@ -93,7 +99,7 @@
       
       DO ig=1,ngrid
-         IF (alim_star2(ig).GT.0.) THEN
-            f(ig) = wmax(ig) * alim_star_tot(ig)                              &
-            &     / (max(1.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
+         IF ((alim_star2(ig) > 0.).and.(plume_height(ig) > 0.)) THEN
+            f(ig) = wmax(ig) * alim_star_tot(ig)                              &  ! => normalization is ok
+            &     / (plume_height(ig) * r_aspect_thermals * alim_star2(ig))
          ELSE
             f(ig) = 0.
@@ -101,4 +107,8 @@
       ENDDO
       
+      print *, 'A*2 ', alim_star2(1)
+      print *, 'A*  ', alim_star_tot(1)
+      print *, 'H   ', plume_height(1)
+      print *, 'wmax', wmax(1)
       
 RETURN
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90	(revision 2143)
@@ -53,5 +53,5 @@
 !     ------
       
-      INTEGER ig, l
+      INTEGER ig, l, k
       INTEGER niter, iter
       
@@ -59,5 +59,4 @@
       REAL qold(ngrid,nlay)
       REAL fqa(ngrid,nlay+1)
-      REAL wqd(ngrid,nlay+1)
       REAL zzm
       
@@ -86,12 +85,14 @@
                print *, 'ERROR: entrainment is greater than the layer mass!'
                print *, 'ig,l,entr', ig, l, entr(ig,l)
+               print *, 'lmin', lmin(ig)
                print *, '-------------------------------'
                print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l)
                print *, '-------------------------------'
-               print *, 'fm+', fm(ig,l+1)
-               print *, 'entr,detr', entr(ig,l), detr(ig,l)
-               print *, 'fm ', fm(ig,l)
-               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
-               print *, 'fm-', fm(ig,l-1)
+               DO k=nlay,1,-1
+                  print *, 'fm ', fm(ig,k+1)
+                  print *, 'entr,detr', entr(ig,k), detr(ig,k)
+               ENDDO
+               print *, 'fm ', fm(ig,1)
+               print *, '-------------------------------'
                CALL abort
             ENDIF
@@ -117,14 +118,4 @@
                qa(ig,l) = q(ig,l)
             ENDIF
-            
-!            IF (qa(ig,l).lt.0.) THEN
-!               print *, 'WARNING: qa is negative!'
-!               print *, 'qa', qa(ig,l)
-!            ENDIF
-            
-!            IF (q(ig,l).lt.0.) THEN
-!               print *, 'WARNING: q is negative!'
-!               print *, 'q', q(ig,l)
-!            ENDIF
          ENDDO
       ENDDO
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90	(revision 2143)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90	(revision 2143)
@@ -0,0 +1,194 @@
+!
+!
+!
+SUBROUTINE thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,detr,masse,fraca,       &
+                         zmax,zmin,u,v,du,dv,ua,va)
+     
+      
+!===============================================================================
+!
+!   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
+!
+! Vectorisation, FH : 2010/03/08
+!
+!===============================================================================
+      
+      USE print_control_mod, ONLY: prt_level,lunout
+      
+      IMPLICIT none
+      
+      
+!===============================================================================
+! Declaration
+!===============================================================================
+      
+!     Inputs:
+!     -------
+      
+      INTEGER ngrid, nlay
+      
+      REAL ptimestep
+      REAL masse(ngrid,nlay)
+      REAL fm(ngrid,nlay+1)
+      REAL entr(ngrid,nlay)
+      REAL detr(ngrid,nlay)
+      REAL fraca(ngrid,nlay+1)
+      REAL zmax(ngrid)
+      REAL zmin(ngrid)
+      REAL u(ngrid,nlay)
+      REAL v(ngrid,nlay)
+      
+!     Outputs:
+!     --------
+      
+      REAL dua(ngrid,nlay)
+      REAL dva(ngrid,nlay)
+      REAL ua(ngrid,nlay)
+      REAL va(ngrid,nlay)
+      REAL du(ngrid,nlay)
+      REAL dv(ngrid,nlay)
+      
+!     Local:
+!     ------
+      
+      INTEGER ig, l
+      INTEGER iter
+      
+      REAL qa(ngrid,nlay)
+      REAL zf
+      REAL zf2
+      REAL wvd(ngrid,nlay+1)
+      REAL wud(ngrid,nlay+1)
+      REAL gamma0(ngrid,nlay+1)
+      REAL gamma(ngrid,nlay+1)
+      REAL ue(ngrid,nlay)
+      REAL ve(ngrid,nlay)
+      REAL plume_height(ngrid)
+      
+      LOGICAL ltherm(ngrid,nlay)
+      
+!===============================================================================
+! Initialization
+!===============================================================================
+      
+      DO ig=1,ngrid
+         plume_height(ig) = zmax(ig) - zmin(ig)
+      ENDDO
+      
+      DO ig=1,ngrid
+         ua(ig,1)=u(ig,1)
+         va(ig,1)=v(ig,1)
+         ue(ig,1)=u(ig,1)
+         ve(ig,1)=v(ig,1)
+      ENDDO
+      
+      gamma0(1:ngrid,1)=0.
+      
+      DO l=2,nlay
+         DO ig=1,ngrid
+            ltherm(ig,l) = (fm(ig,l+1) + detr(ig,l)) * ptimestep > 1.e-5 * masse(ig,l)
+            IF (ltherm(ig,l).and.(zmax(ig) > 0.)) THEN
+               gamma0(ig,l) = masse(ig,l) * 0.5 / zmax(ig)                    &
+               &            * sqrt(0.5 * (fraca(ig,l+1) + fraca(ig,l)))
+            ELSE
+               gamma0(ig,l) = 0.
+            ENDIF
+         ENDDO
+      ENDDO
+      
+      gamma(:,:) = 0.
+      
+!===============================================================================
+! 
+!===============================================================================
+      
+      DO l=2,nlay
+         
+         DO ig=1,ngrid
+            IF (ltherm(ig,l)) THEN
+               dua(ig,l)=ua(ig,l-1)-u(ig,l-1)
+               dva(ig,l)=va(ig,l-1)-v(ig,l-1)
+            ELSE
+               ua(ig,l)=u(ig,l)
+               va(ig,l)=v(ig,l)
+               ue(ig,l)=u(ig,l)
+               ve(ig,l)=v(ig,l)
+            ENDIF
+         ENDDO
+         
+         DO iter=1,5
+            DO ig=1,ngrid
+! Calcul prenant en compte la fraction reelle
+               zf = 0.5 * (fraca(ig,l) + fraca(ig,l+1))
+               zf2 = 1./(1.-zf)
+! Calcul avec fraction infiniement petite
+!               zf = 0.
+!               zf2 = 1.
+               
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! FH: La première fois on multiplie le coefficient de freinage par le module du
+!     vent dans la couche en dessous. Mais pourquoi donc ?
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               IF (ltherm(ig,l)) THEN
+!   On choisit une relaxation lineaire.
+!                  gamma(ig,l) = gamma0(ig,l)
+!   On choisit une relaxation quadratique.
+                  gamma(ig,l) = gamma0(ig,l) * sqrt(dua(ig,l)**2 + dva(ig,l)**2)
+                  
+                  ua(ig,l) = (fm(ig,l) * ua(ig,l-1)                           &
+                  &        + (zf2 * entr(ig,l) + gamma(ig,l)) * u(ig,l))      &
+                  &        / (fm(ig,l+1) + detr(ig,l) + entr(ig,l) * zf * zf2 &
+                  &        + gamma(ig,l))
+                  va(ig,l) = (fm(ig,l) * va(ig,l-1)                           &
+                  &        + (zf2 * entr(ig,l) + gamma(ig,l)) * v(ig,l))      &
+                  &        / (fm(ig,l+1) + detr(ig,l) + entr(ig,l) * zf * zf2 &
+                  &        + gamma(ig,l))
+                  dua(ig,l) = ua(ig,l) - u(ig,l)
+                  dva(ig,l) = va(ig,l) - v(ig,l)
+                  ue(ig,l) = (u(ig,l) - zf * ua(ig,l)) * zf2
+                  ve(ig,l) = (v(ig,l) - zf * va(ig,l)) * zf2
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      
+!-------------------------------------------------------------------------------
+! Calcul du flux vertical de moment dans l'environnement
+!-------------------------------------------------------------------------------
+      
+      DO l=2,nlay
+         DO ig=1,ngrid
+            wud(ig,l) = fm(ig,l) * ue(ig,l)
+            wvd(ig,l) = fm(ig,l) * ve(ig,l)
+         ENDDO
+      ENDDO
+      
+      DO ig=1,ngrid
+         wud(ig,1) = 0.
+         wud(ig,nlay+1) = 0.
+         wvd(ig,1) = 0.
+         wvd(ig,nlay+1) = 0.
+      ENDDO
+      
+!===============================================================================
+! Tendencies
+!===============================================================================
+      
+      DO l=1,nlay
+         DO ig=1,ngrid
+            du(ig,l) = ((detr(ig,l) + gamma(ig,l)) * ua(ig,l)                 &
+            &        - (entr(ig,l) + gamma(ig,l)) * ue(ig,l)                  &
+            &        - wud(ig,l) + wud(ig,l+1))                               &
+            &        / masse(ig,l)
+            dv(ig,l) = ((detr(ig,l) + gamma(ig,l)) * va(ig,l)                 &
+            &        - (entr(ig,l) + gamma(ig,l)) * ve(ig,l)                  &
+            &        - wvd(ig,l) + wvd(ig,l+1))                               &
+            &        / masse(ig,l)
+         ENDDO
+      ENDDO
+      
+      
+RETURN
+END
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90	(revision 2143)
@@ -51,6 +51,6 @@
 !     --------
       
-      REAL zqt(ngrid,nlay)                      ! q    environment
-      REAL zql(ngrid,nlay)                      ! q    environment
+      REAL zqt(ngrid,nlay)                      ! qt   environment
+      REAL zql(ngrid,nlay)                      ! ql   environment
       REAL zt(ngrid,nlay)                       ! T    environment
       REAL ztv(ngrid,nlay)                      ! TRPV environment
@@ -76,4 +76,7 @@
       zhl(:,:) = pt(:,:) / zpopsk(:,:)
       
+      zqt(:,:) = pq(:,:,igcm_h2o_vap)
+      zql(:,:) = 0.
+      
 !===============================================================================
 ! Condensation and latent heat release
@@ -81,6 +84,4 @@
       
       IF (water) THEN
-         
-         zqt(:,:) = pq(:,:,igcm_h2o_vap)
          
          DO l=1,nlay
@@ -101,8 +102,6 @@
       ELSE
          
-         zqt(:,:) = 0.
-         zql(:,:) = 0.
          zt(:,:) = pt(:,:)
-         ztv(:,:) = zhl(:,:)
+         ztv(:,:) = pt(:,:) / zpopsk(:,:)
          
       ENDIF
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90	(revision 2143)
@@ -4,9 +4,9 @@
 SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
                           lmin,lmax,entr_star,detr_star,                      &
-                          f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
-      
-      
-!===============================================================================
-!  Purpose: fluxes deduction
+                          f,rhobarz,zlev,zw2,fm,entr,detr)
+      
+      
+!===============================================================================
+!  Purpose: deduction des flux
 !  
 !  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
@@ -39,6 +39,4 @@
       REAL rhobarz(ngrid,nlay)
       REAL f(ngrid)
-      REAL zqla(ngrid,nlay)
-      REAL zmax(ngrid)
       
 !     Outputs:
@@ -53,5 +51,5 @@
       
       INTEGER ig, l, k
-      INTEGER igout, lout                 ! Error grid point number and level
+      INTEGER igout, lout                 ! Error grid point and level
       
       REAL zfm                            ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max
@@ -62,21 +60,23 @@
       REAL ddd                            ! ddd0 - eee
       REAL zzz                            ! Temporary variable set to mass flux
-      REAL fact                           ! Factor between old norm and new norm
-      
-      INTEGER ncorecdetr                  ! detr > fm-   counter
-      INTEGER ncorecentr                  ! entr > e_max counter
-      INTEGER ncorecalpha                 ! fm > zfm     counter
+      REAL fact
+      REAL test
+      
+      INTEGER ncorecentr
+      INTEGER ncorecdetr
+      INTEGER nerrorequa
+      INTEGER ncorecfact
+      INTEGER ncorecalpha
       
       LOGICAL labort_physic
       
-      CHARACTER (len=20) :: modname='thermcell_flux'
-      CHARACTER (len=80) :: abort_message
-      
 !===============================================================================
 ! Initialization
 !===============================================================================
       
-      ncorecdetr  = 0
-      ncorecentr  = 0
+      nerrorequa = 0
+      ncorecentr = 0
+      ncorecdetr = 0
+      ncorecfact = 0
       ncorecalpha = 0
       
@@ -85,6 +85,10 @@
       fm(:,:)   = 0.
       
-!===============================================================================
-! Mass flux, entrainment and detrainment
+      labort_physic = .false.
+      
+      fact = 0.
+      
+!===============================================================================
+! Calcul de l'entrainement, du detrainement et du flux de masse
 !===============================================================================
       
@@ -99,12 +103,12 @@
       
 !-------------------------------------------------------------------------------
-! Mass flux and boundary conditions
+! Calcul du flux de masse
 !-------------------------------------------------------------------------------
       
       DO l=1,nlay
          DO ig=1,ngrid
-            IF (l.lt.lmax(ig) .and. l.ge.lmin(ig)) THEN
+            IF (l < lmax(ig) .and. l >= lmin(ig)) THEN
                fm(ig,l+1) = fm(ig,l) + entr(ig,l) - detr(ig,l)
-            ELSEIF (l.eq.lmax(ig)) THEN
+            ELSEIF (l == lmax(ig)) THEN
                fm(ig,l+1) = 0.
                detr(ig,l) = fm(ig,l) + entr(ig,l)
@@ -118,5 +122,5 @@
       
 !===============================================================================
-! Validity tests and corrections
+! Checking
 !===============================================================================
       
@@ -127,16 +131,6 @@
 !-------------------------------------------------------------------------------
          
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : I remove the correction and replace it by an uncompromising test.
-!      According to the previous derivations, we MUST have positive mass flux
-!      everywhere! Indeed, as soon as fm becomes negative, the plume stops.
-!      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.
-         
-         DO ig=1,ngrid
-            IF (fm(ig,l).lt.0.) THEN
+         DO ig=1,ngrid
+            IF (fm(ig,l) < 0.) THEN
                labort_physic = .true.
                igout = ig
@@ -145,32 +139,10 @@
          ENDDO
          
-         IF (labort_physic) THEN
-            print *, 'ERROR: mass flux has negative value(s)!'
-            print *, 'ig,l,fm',igout, lout, fm(igout,lout)
-            print *, 'lmin,lmax', lmin(igout), lmax(igout)
-            print *, '-------------------------------'
-            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
-         
 !-------------------------------------------------------------------------------
 ! Is entrainment positive ?
 !-------------------------------------------------------------------------------
          
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! 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.
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         
-         DO ig=1,ngrid
-            IF (entr(ig,l).lt.0.) THEN
+         DO ig=1,ngrid
+            IF (entr(ig,l) < 0.) THEN
                labort_physic = .true.
                igout = ig
@@ -179,9 +151,26 @@
          ENDDO
          
+!-------------------------------------------------------------------------------
+! Is detrainment positive ?
+!-------------------------------------------------------------------------------
+         
+         DO ig=1,ngrid
+            IF (detr(ig,l) < 0.) THEN
+               labort_physic = .true.
+               igout = ig
+               lout = l
+            ENDIF
+         ENDDO
+         
+!-------------------------------------------------------------------------------
+! Abort
+!-------------------------------------------------------------------------------
+         
          IF (labort_physic) THEN
-            print *, 'ERROR: entrainment has negative value(s)!'
+            print *, '---------------------------------------------------------'
+            print *, 'ERROR: mass flux has negative value(s)!'
             print *, 'ig,l,entr', igout, lout, entr(igout,lout)
             print *, 'lmin,lmax', lmin(igout), lmax(igout)
-            print *, '-------------------------------'
+            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
             DO k=nlay,1,-1
                print *, 'fm ', fm(igout,k+1)
@@ -189,38 +178,10 @@
             ENDDO
             print *, 'fm ', fm(igout,1)
-            print *, '-------------------------------'
-            abort_message = 'Negative entrainment in thermcell_flux'
-            CALL abort_physic(modname,abort_message,1)
-         ENDIF
-         
-!-------------------------------------------------------------------------------
-! Is detrainment positive ?
-!-------------------------------------------------------------------------------
-         
-         DO ig=1,ngrid
-            IF (detr(ig,l).lt.0.) THEN
-               labort_physic = .true.
-               igout = ig
-               lout = l
-            ENDIF
-         ENDDO
-         
-         IF (labort_physic) THEN
-            print *, 'ERROR: detrainment has negative value(s)!'
-            print *, 'ig,l,detr', igout, lout, detr(igout,lout)
-            print *, 'lmin,lmax', lmin(igout), lmax(igout)
-            print *, '-------------------------------'
-            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
-         
-!-------------------------------------------------------------------------------
-! Is detrainment lesser than incoming mass flux ?
+            print *, '---------------------------------------------------------'
+            CALL abort
+         ENDIF
+         
+!-------------------------------------------------------------------------------
+! Is detrainment lessser than incoming mass flux ?
 !-------------------------------------------------------------------------------
          
@@ -238,23 +199,10 @@
                detr(ig,l) = fm(ig,l)
                entr(ig,l) = fm(ig,l+1)
-               
-               IF (prt_level.ge.10) THEN
-                  print *, 'WARNING: Detrainment is modified in thermcell_flux!'
-                  print *, 'ig,l,lmax', ig, l, lmax(ig)
-               ENDIF
-               
-               IF (prt_level.ge.100) THEN
-                  print *, 'fm+', fm(ig,l+1)
-                  print *, 'entr,detr', entr(ig,l), detr(ig,l)
-                  print *, 'fm-', fm(ig,l)
-               ENDIF
-               
                ncorecdetr  = ncorecdetr + 1
-               
-            ENDIF
-         ENDDO
-         
-!-------------------------------------------------------------------------------
-! Is entrainment mass fraction lesser than fomass_max ?
+            ENDIF
+         ENDDO
+         
+!-------------------------------------------------------------------------------
+! Is entrained mass lesser than fomass_max ?
 !-------------------------------------------------------------------------------
          
@@ -271,6 +219,4 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          
-         labort_physic=.false.
-         
          DO ig=1,ngrid
             eee0 = entr(ig,l)
@@ -278,19 +224,12 @@
             eee  = entr(ig,l) - masse(ig,l) * fomass_max / ptimestep
             ddd  = detr(ig,l) - eee
-            
-            IF (eee.gt.0.) THEN
+            IF (eee > 0.) THEN
                entr(ig,l) = entr(ig,l) - eee
                ncorecentr  = ncorecentr + 1
-               
-               IF (prt_level.ge.10) THEN
-                  print *, 'WARNING: Entrainment is modified in thermcell_flux!'
-                  print *, 'ig,l,lmax', ig, l, lmax(ig)
-               ENDIF
-               
-               IF (ddd.gt.0.) THEN                          ! detr in the current layer is reduced
+               IF (ddd > 0.) THEN
                   detr(ig,l) = ddd
-               ELSEIF (l.eq.lmax(ig)) THEN                  ! detr in the last layer is adjusted
+               ELSEIF (l == lmax(ig)) THEN
                   detr(ig,l) = fm(ig,l) + entr(ig,l)
-               ELSEIF (entr(ig,l+1).gt.ABS(ddd)) THEN       ! detr in the current layer is set to 0 and entr in the layer above is reduced
+               ELSEIF (entr(ig,l+1) > ABS(ddd)) THEN
                   detr(ig,l) = 0.
                   fm(ig,l+1) = fm(ig,l) + entr(ig,l)
@@ -298,91 +237,54 @@
 !               ELSE
 !                  entr(ig,l) = entr(ig,l) + eee
-!                  igout=ig
-!                  lout=l
-!                  labort_physic=.true.
+!                  igout = ig
+!                  lout = l
+!                  labort_physic = .true.
                ELSE
-                  fact = (eee0 - eee) / eee0
+                  fact = max(fact, eee0 / (eee0 - eee))
                   entr(ig,l) = eee0
-                  detr(ig,:) = detr(ig,:) * fact
-                  entr(ig,:) = entr(ig,:) * fact
-                  fm(ig,:)   = fm(ig,:) * fact
-                  
-                  IF (prt_level.ge.1) THEN
-                     print *, 'WARNING: Normalisation is modified in thermcell_flux!'
-                     print *, 'old f, new f :', f(ig), fact * f(ig)
-                  ENDIF
+                  ncorecfact = ncorecfact + 1
                ENDIF
             ENDIF
          ENDDO
          
-!         IF (labort_physic) THEN
-!            print *, 'ERROR: Entrainment is greater than its maximal authorized value!'
-!            print *, '       Nor detrainment neither entrainment can compensate it!'
-!            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
-!            print *, 'lmin,lmax', lmin(igout), lmax(igout)
-!            print *, '--------------------'
-!            print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep
-!            print *, '   fomass_max :', fomass_max
-!            print *, '   masse      :', masse(igout,lout)
-!            print *, '   ptimestep  :', ptimestep
-!            print *, 'norm  :', f(igout)
-!            print *, 'entr* :', entr_star(igout,lout)
-!            print *, '--------------------'
-!            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 = 'Entrainment is too large in thermcell_flux'
-!            CALL abort_physic (modname,abort_message,1)
-!         ENDIF
-         
-!-------------------------------------------------------------------------------
-! Is updraft fraction lesser than alpha_max ?
-!-------------------------------------------------------------------------------
-         
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! FH : Partie a revisiter.
-!      Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1
-!      - soit limiter la hauteur du thermique en considerant que c'est 
-!        la derniere chouche
-!      - soit limiter F a rho w.
-!      Dans le second cas, il faut en fait limiter a un peu moins que ca parce
-!      qu'on a des 1/(1-alpha) un peu plus loin dans thermcell_main et qu'il
-!      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.  
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
+         IF (labort_physic) THEN
+            print *, '---------------------------------------------------------'
+            print *, 'ERROR: Entrainment is greater than maximal authorized value!'
+            print *, '       Nor detrainment neither entrainment can compensate it!'
+            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
+            print *, 'lmin,lmax', lmin(igout), lmax(igout)
+            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
+            print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep
+            print *, '   fomass_max :', fomass_max
+            print *, '   masse      :', masse(igout,lout)
+            print *, '   ptimestep  :', ptimestep
+            print *, 'norm  :', f(igout)
+            print *, 'entr* :', entr_star(igout,lout)
+            print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
+            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 *, '---------------------------------------------------------'
+            CALL abort
+         ENDIF
+         
+!-------------------------------------------------------------------------------
+! Is Updraft fraction lesser than alpha_max ?
+!-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
             zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max
             
-            IF (fm(ig,l+1) .gt. zfm) THEN
+            IF (fm(ig,l+1) > zfm) THEN
                f_old = fm(ig,l+1)
                fm(ig,l+1) = zfm
                detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
-               
+               ncorecalpha = ncorecalpha + 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)
-               ENDIF
-               
-               IF (prt_level.ge.10) THEN
-                  print *, 'Mass flux is modified in thermcell_flux'
-                  print *, 'ig,l,lmax', ig, l, lmax(ig)
-               ENDIF
-               
-               IF (prt_level.ge.100) THEN
-                  print *, 'fm-', fm(ig,l)
-                  print *, 'entr,detr', entr(ig,l), detr(ig,l)
-                  print *, 'fm+', fm(ig,l+1)
-               ENDIF
-               
-               ncorecalpha = ncorecalpha + 1
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+               entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
             ENDIF
          ENDDO
@@ -390,14 +292,29 @@
       ENDDO
       
-!-------------------------------------------------------------------------------
-! Boundary conditions
-!-------------------------------------------------------------------------------
-      
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
-!      Theoretically, we always get lmax >= lmin >= linf > 0
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (fact > 0.) THEN
+         entr(:,:) = entr(:,:) / fact
+         detr(:,:) = entr(:,:) / fact
+         fm(:,:) = entr(:,:) / fact
+      ENDIF
+      
+!-------------------------------------------------------------------------------
+! Is equation df/dz = e - d still verified ?
+!-------------------------------------------------------------------------------
+      
+!      DO l=1,nlay
+!         DO ig=1,ngrid
+!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
+!            IF (test > 1.e-10) THEN
+!               nerrorequa = nerrorequa + 1
+!            ENDIF
+!         ENDDO
+!      ENDDO
+      
+!-------------------------------------------------------------------------------
+! Reset top boundary conditions
+!-------------------------------------------------------------------------------
+      
       DO ig=1,ngrid
-         IF (lmax(ig).gt.0) THEN
+         IF (lmax(ig).GT.0) THEN
             detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig))
             fm(ig,lmax(ig)+1) = 0.
@@ -406,43 +323,36 @@
       ENDDO
       
-!-------------------------------------------------------------------------------
-! Corrections display
-!-------------------------------------------------------------------------------
-      
-      IF (prt_level.GE.1) THEN
-                  
-         IF (ncorecdetr.GE.1) THEN
+!===============================================================================
+! Outputs
+!===============================================================================
+      
+      IF (prt_level > 0) THEN
+         
+         IF (ncorecdetr.ge.1) THEN
             print *, 'WARNING: Detrainment is greater than mass flux!'
-            print *, '         In', ncorecdetr, 'point(s)'
-         ENDIF
-         
-         IF (ncorecentr.GE.1) THEN
+            print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid
+         ENDIF
+         
+         IF (ncorecentr.ge.1) THEN
             print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
-            print *, '         In', ncorecentr, 'point(s)'
-         ENDIF
-         
-         IF (ncorecalpha.GE.1) THEN
+            print *, 'in', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid
+         ENDIF
+         
+         IF (ncorecfact.ge.1) THEN
+            print *, 'WARNING: Entrained mass needs renormalization to be ok!'
+            print *, 'in', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid
+         ENDIF
+         
+!         IF (nerrorequa.ge.1) THEN
+!            print *, 'WARNING: !'
+!            print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid
+!         ENDIF
+         
+         IF (ncorecalpha.ge.1) THEN
             print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
-            print *, '         In', ncorecalpha, 'point(s)'
+            print *, 'in', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid
          ENDIF
          
       ENDIF
-      
-! AB : temporary test added to check the validity of eq. df/dz = e - d
-!      DO l=1,nlay
-!         DO ig=1,ngrid
-!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
-!            IF (test.gt.1.e-10) THEN
-!               print *, 'WARNING: df/dz != entr - detr'
-!               print *, 'ig,l', ig, l
-!               print *, 'fm+  :', fm(ig,l+1)
-!               print *, 'entr,detr', entr(ig,l), detr(ig,l)
-!               print *, 'fm   :', fm(ig,l)
-!               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
-!               print *, 'fm-  :', fm(ig,l-1)
-!               print *, 'err. :', test
-!            ENDIF
-!         ENDDO
-!      ENDDO
       
       
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90	(revision 2143)
@@ -2,6 +2,6 @@
 !
 !
-SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
-                            zlev,lmax,zmax,zmix,wmax,f_star)
+SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
+                            zlev,zmin,zmix,zmax,wmax,f_star)
       
       
@@ -9,7 +9,4 @@
 !  Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
 !===============================================================================
-      
-      USE thermcell_mod
-      USE print_control_mod, ONLY: prt_level
       
       IMPLICIT NONE
@@ -36,8 +33,9 @@
       
       REAL linter(ngrid)
-      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)
+      REAL zmin(ngrid)                          ! Altitude of the plume bottom
+      REAL zmax(ngrid)                          ! Altitude of the plume top
+      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
+      REAL wmax(ngrid)                          ! Maximal vertical speed
+      REAL zw2(ngrid,nlay+1)                    ! Square vertical speed (becomes its square root)
       
 !     Local:
@@ -46,8 +44,5 @@
       INTEGER ig, l
       
-      REAL num(ngrid)
-      REAL denom(ngrid)
       REAL zlevinter(ngrid)
-      REAL zlev2(ngrid,nlay+1)
       
 !===============================================================================
@@ -69,11 +64,17 @@
       ENDDO
       
-      DO ig=1,ngrid
-         zlev2(ig,:) = zlev(ig,:) - zlev(ig,lmin(ig))
-      ENDDO
-      
-!===============================================================================
-! Calcul de la hauteur max du thermique
-!===============================================================================
+!===============================================================================
+! Thermal plume height
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
+! Calcul de zmin
+!-------------------------------------------------------------------------------
+      
+      DO l=1,nlay
+         DO ig=1,ngrid
+            zmin(ig) = zlev(ig,lmin(ig))
+         ENDDO
+      ENDDO
       
 !-------------------------------------------------------------------------------
@@ -83,5 +84,5 @@
       DO ig=1,ngrid
          DO l=nlay,lmin(ig)+1,-1
-            IF (zw2(ig,l).le.0..or.f_star(ig,l).le.0.) THEN
+            IF ((zw2(ig,l) <= 0.).or.(f_star(ig,l) <= 0.)) THEN
                lmax(ig) = l - 1
             ENDIF
@@ -90,9 +91,8 @@
       
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! 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.0.) THEN
+! AB: Problematic case where thermal plume reaches the top of the model...
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      DO ig=1,ngrid
+         IF (zw2(ig,nlay) > 0.) THEN
             print *, 'WARNING: a thermal plume reaches the highest layer!'
             print *, 'ig,l', ig, nlay
@@ -110,9 +110,7 @@
 !      Otherwise we have lmax>lmin.
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      
       DO ig=1,ngrid
          l = lmax(ig)
-         
-         IF (l==lmin(ig)) THEN
+         IF (l == lmin(ig)) THEN
             linter(ig) = l
          ELSE
@@ -121,20 +119,12 @@
       ENDDO
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! 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
-         zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig)   &
-         &             + zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)      &
-         &             - zlev(ig,lmax(ig)))
-         zmax(ig) = max(zmax(ig), zlevinter(ig) - zlev(ig,lmin(ig)))
-      ENDDO
-      
-!===============================================================================
-! Calcul de wmax et du niveau d'inversion.
+      DO ig=1,ngrid
+         zlevinter(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig))          &
+         &             * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig)))
+         zmax(ig) = max(zmax(ig), zlevinter(ig))
+      ENDDO
+      
+!===============================================================================
+! Thermal plume maximal speed and inversion layer
 !===============================================================================
       
@@ -145,5 +135,5 @@
       DO l=1,nlay
          DO ig=1,ngrid
-            IF(l.le.lmax(ig) .and. l.gt.lmin(ig)) THEN
+            IF((l <= lmax(ig)).and.(l > lmin(ig))) THEN
                IF (zw2(ig,l).lt.0.) THEN
                   print *, 'ERROR: zw2 has negative value(s)!'
@@ -158,5 +148,5 @@
                zw2(ig,l) = sqrt(zw2(ig,l))
                
-               IF (zw2(ig,l).gt.wmax(ig)) THEN
+               IF (zw2(ig,l) > wmax(ig)) THEN
                   wmax(ig)  = zw2(ig,l)
                   lmix(ig)  = l
@@ -172,15 +162,10 @@
 !-------------------------------------------------------------------------------
       
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-! 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
-         IF (lmix(ig).gt.lmin(ig)) THEN
+      DO ig=1,ngrid
+         IF (lmix(ig) > lmin(ig)) THEN
             IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))                        &
             &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))             &
             &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                   &
-            &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)  &
+            &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))) > 1e-10)   &
             &        THEN
                zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))              &
@@ -191,8 +176,7 @@
                &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))          &
                &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                &
-               &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))        &
-               &        - zlev(ig,lmin(ig))
+               &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
             ELSE
-               zmix(ig) = zlev(ig,lmix(ig)) - zlev(ig,lmin(ig))
+               zmix(ig) = zlev(ig,lmix(ig))
                print *, 'WARNING: problematic zmix value!'
             ENDIF
@@ -203,12 +187,16 @@
       
 !===============================================================================
-! Verification zmax > zmix
-!===============================================================================
-      
-      DO ig=1,ngrid
-         IF ((zmax(ig)-zmix(ig)).lt.0.) THEN
+! Check consistence between zmax and zmix
+!===============================================================================
+      
+!-------------------------------------------------------------------------------
+! 
+!-------------------------------------------------------------------------------
+      
+      DO ig=1,ngrid
+         IF ((zmax(ig)-zmix(ig)) < 0.) THEN
+            print *, 'WARNING: we have zmix > zmax!'
+            print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig)
             zmix(ig) = 0.9 * zmax(ig)
-            print *, 'WARNING: we have zmix > zmax!'
-            print *, 'zmax,zmix', zmax(ig), zmix(ig)
          ENDIF
       ENDDO
@@ -220,5 +208,5 @@
       DO ig=1,ngrid
          DO l=1,nlay
-            IF (zmix(ig).ge.zlev2(ig,l).and.zmix(ig).lt.zlev2(ig,l+1)) THEN
+            IF ((zmix(ig) >= zlev(ig,l)).and.(zmix(ig) < zlev(ig,l+1))) THEN
                lmix(ig) = l
             ENDIF
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90	(revision 2143)
@@ -49,6 +49,5 @@
       
       USE thermcell_mod
-      USE tracer_h, ONLY: igcm_h2o_vap
-      USE print_control_mod, ONLY: lunout, prt_level
+      USE print_control_mod, ONLY: prt_level
       
       IMPLICIT NONE
@@ -96,4 +95,8 @@
       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
@@ -124,28 +127,20 @@
       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 zmix(ngrid)                          ! Altitude of maximal vertical speed
-      REAL zmax(ngrid)                          ! Plume height
       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
       
 !===============================================================================
@@ -178,14 +173,10 @@
       pdqadj(:,:,:) = 0.0
       
-! 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)
-      ENDDO
-      
-      IF (prt_level.ge.20) then
-         DO ig=1,ngrid
-            print *, 'ig,f0', ig, f0(ig)
-         ENDDO
-      ENDIF
+! 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
       
 !===============================================================================
@@ -212,5 +203,5 @@
       
       DO l=1,nlay
-         zlay(:,l) = pphi(:,l)/RG
+         zlay(:,l) = pphi(:,l) / RG
       ENDDO
       
@@ -222,5 +213,7 @@
       
       IF (prt_level.ge.10) THEN
-         write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)'
+         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
       
@@ -319,7 +312,7 @@
       &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
       &                    detr_star,entr_star,f_star,                        &
-      &                    ztva,zhla,zqla,zqta,zta,                           &
-      &                    zw2,zqsa,lmix,lmin)
-            
+      &                    ztva,zhla,zqla,zqta,zta,zqsa,                      &
+      &                    zw2,lmix,lmin)
+      
 !-------------------------------------------------------------------------------
 ! Thermal plumes characteristics: zmax, zmix, wmax
@@ -327,6 +320,6 @@
       
 ! 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)
+      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
+      &                     zlev,zmin,zmix,zmax,wmax,f_star)
       
 !===============================================================================
@@ -339,5 +332,5 @@
       
       CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
-      &                      lmax,entr_star,zmax,wmax,f)
+      &                      lmax,entr_star,zmin,zmax,wmax,f)
       
       IF (tau_thermals>1.) THEN
@@ -352,7 +345,7 @@
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       IF (.not. (f0(1).ge.0.) ) THEN
-         abort_message = '.not. (f0(1).ge.0.)'
+         print *, 'ERROR: mass flux norm is not positive!'
          print *, 'f0 =', f0(1)
-         CALL abort_physic(modname,abort_message,1)
+         CALL abort
       ENDIF
       
@@ -393,5 +386,5 @@
       DO l=2,nlay
          DO ig=1,ngrid
-            IF (zw2(ig,l).gt.0.) THEN
+            IF (zw2(ig,l) > 0.) THEN
                fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
             ELSE
@@ -433,8 +426,10 @@
       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)
       
+!      CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,fraca,    &
+!      &                  zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva)
       
 RETURN
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90	(revision 2143)
@@ -5,18 +5,18 @@
       
 ! Flags for computations
-                                                                  !  default
-      INTEGER,PARAMETER :: dqimpl                     = 1         !  1        flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme)
+                                                 !  default
+      INTEGER,SAVE :: dqimpl                     !  1        flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme)
       
 ! Physical parameters
       
-      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                      = 1.0       !  0.9       - factor between 0 (e=d) and 1 (rho*fraca=cst)
-      REAL,PARAMETER :: afact                         = 1.        !  2./3.    Buoyancy contribution - factor between 0 and 1
-      REAL,PARAMETER :: fact_epsilon                  = 1.e-4     !  2.e-3    Friction at plume borders - exponential decrease
-      REAL,PARAMETER :: nu                            = 0.000     !           Geometrical contributions to entrainment and detrainment
-      REAL,PARAMETER :: alpha_max                     = 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     Minimal permitted pressure to trigger a thermal plume
+      REAL,SAVE :: r_aspect_thermals             !  1.0      Aspect ratio of the thermals (width / height)
+      REAL,SAVE :: tau_thermals                  !  0.       Relaxation time
+      REAL,SAVE :: betalpha                      !  0.9      - between 0 (e=d) and 1 (rho*fraca=cst)
+      REAL,SAVE :: afact                         !  2./3.    - buoyancy acceleration efficiency, between 0 and 1
+      REAL,SAVE :: fact_epsilon                  !  2.e-3    - turbulence and friction deceleration
+      REAL,SAVE :: nu                            !  0.       Geometrical contributions to entrainment and detrainment
+      REAL,SAVE :: alpha_max                     !  0.7      Maximal permitted updraft fraction
+      REAL,SAVE :: fomass_max                    !  0.5      Maximal permitted outgoing layer mass fraction
+      REAL,SAVE :: pres_limit                    !  1.e5     -
       
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -29,5 +29,5 @@
 !      ngrid. 
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      INTEGER,PARAMETER :: linf                       = 1         !     1
+      INTEGER,SAVE :: linf                       !     1
       
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -35,5 +35,5 @@
 !      layer linf which can be used to force convection to begin in it.
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      REAL,PARAMETER :: d_temp                        = 0.        !     0.
+      REAL,SAVE :: d_temp                        !     0.
       
 ! Physical constants
@@ -52,4 +52,6 @@
       SUBROUTINE init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV)
          
+         USE ioipsl_getin_p_mod, ONLY: getin_p
+         
          IMPLICIT NONE
          
@@ -67,4 +69,59 @@
          RD = r
          
+         print *, 'Plume aspect ratio ?'
+         r_aspect_thermals = 1.
+         call getin_p('r_aspect_thermals', r_aspect_thermals)
+         print *, 'r_aspect_thermals = ', r_aspect_thermals
+         
+         print *, 'Relaxation time ?'
+         tau_thermals = 0.
+         call getin_p('tau_thermals', tau_thermals)
+         print *, 'tau_thermals = ', tau_thermals
+         
+         print *, 'betalpha ?'
+         betalpha = 0.9
+         call getin_p('betalpha', betalpha)
+         print *, 'betalpha = ', betalpha
+         
+         print *, 'Buoyancy acceleration efficiency ?'
+         afact = 2./3.
+         call getin_p('afact', afact)
+         print *, 'afact = ', afact
+         
+         print *, 'Turbulence and friction factor ?'
+         fact_epsilon = 2.e-3
+         call getin_p('fact_epsilon', fact_epsilon)
+         print *, 'fact_epsilon = ', fact_epsilon
+         
+         print *, 'Constant entrainment and detrainment term ?'
+         nu = 0.
+         call getin_p('nu', nu)
+         print *, 'nu = ', nu
+         
+         print *, 'Maximal authorized updraft fraction ?'
+         alpha_max = 0.7
+         call getin_p('alpha_max', alpha_max)
+         print *, 'alpha_max = ', alpha_max
+         
+         print *, 'Maximal authorized entrained mass flux ?'
+         fomass_max = 0.5
+         call getin_p('fomass_max', fomass_max)
+         print *, 'fomass_max = ', fomass_max
+         
+         print *, 'Minimal pressure below which plume can no longer trigger ?'
+         pres_limit = 1.e5
+         call getin_p('pres_limit', pres_limit)
+         print *, 'pres_limit = ', pres_limit
+         
+         print *, 'Deepest layer from which a plume can rise ?'
+         linf = 1
+         call getin_p('linf', linf)
+         print *, 'linf = ', linf
+         
+         print *, 'd_temp ?'
+         d_temp = 0.
+         call getin_p('d_temp', d_temp)
+         print *, 'd_temp = ', d_temp
+         
          RETURN
       END SUBROUTINE init_thermcell_mod
Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90	(revision 2142)
+++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90	(revision 2143)
@@ -5,6 +5,6 @@
                            zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
                            detr_star,entr_star,f_star,                        &
-                           ztva,zhla,zqla,zqta,zta,                           &
-                           zw2,zqsa,lmix,lmin)
+                           ztva,zhla,zqla,zqta,zta,zqsa,                      &
+                           zw2,lmix,lmin)
       
       
@@ -62,9 +62,9 @@
       
       REAL ztva(ngrid,nlay)                     ! TRPV plume (after mixing)
-      REAL zhla(ngrid,nlay)                     ! TP   plume (after mixing)
+      REAL zhla(ngrid,nlay)                     ! TP   plume ?
       REAL zqla(ngrid,nlay)                     ! ql   plume (after mixing)
-      REAL zqta(ngrid,nlay)                     ! qt   plume (after mixing)
+      REAL zqta(ngrid,nlay)                     ! qt   plume ?
       REAL zqsa(ngrid,nlay)                     ! qsat plume (after mixing)
-      REAL zw2(ngrid,nlay+1)                    ! w2   plume (after mixing)
+      REAL zw2(ngrid,nlay+1)                    ! w    plume (after mixing)
       
 !     Local:
@@ -77,5 +77,5 @@
       REAL zta_est(ngrid,nlay)                  ! TR   plume (before mixing)
       REAL zqsa_est(ngrid)                      ! qsat plume (before mixing)
-      REAL zw2_est(ngrid,nlay+1)                ! w2   plume (before mixing)
+      REAL zw2_est(ngrid,nlay+1)                ! w    plume (before mixing)
       
       REAL zta(ngrid,nlay)                      ! TR   plume (after mixing)
@@ -128,4 +128,6 @@
       ztv2(:,linf)     = ztv(:,linf) + d_temp
       
+      active(:) = .false.
+      
 !===============================================================================
 ! First layer computation
@@ -133,22 +135,21 @@
       
       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))
-            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.
+         DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
+            zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
+            zdz = zlev(ig,l+1) - zlev(ig,l)
+            zw2m = afact * zbuoy(ig,l) * zdz / (1. + betalpha)
+!            gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
+!            test = gamma / zw2m - nu
+            test = zbuoy(ig,l)
+            IF (test > 0.) THEN
                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.
+!               entr_star(ig,l) = zdz * f_star(ig,l) * zbetalpha * gamma / zw2m - nu ! Problem because f*(ig,l) = 0
+!               detr_star(ig,l) = f_star(ig,l) * nu                                  ! Problem because f*(ig,l) = 0
+!               f_star(ig,l+1)  = entr_star(ig,l) - detr_star(ig,l)
+               entr_star(ig,l) = 1.
+               f_star(ig,l+1) = 1.
+               zw2_est(ig,l+1) = zw2m * 2.
+               zw2(ig,l+1) = zw2_est(ig,l+1)
             ENDIF
             l = l + 1
@@ -163,21 +164,16 @@
          
 !-------------------------------------------------------------------------------
-! 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
-               active(ig) = .true.
-            ENDIF
-            
-            active(ig) = active(ig)                                           &
-            &          .and. (zw2(ig,l).GT.1.e-10)                            &
-            &          .and. (f_star(ig,l) + entr_star(ig,l)).GT.1.e-10
-         ENDDO
+! Is thermal plume (still) active ?
+!-------------------------------------------------------------------------------
+         
+         DO ig=1,ngrid
+            active(ig) = (active(ig).or.(l == lmin(ig)+1))                    &
+            &           .and.(zw2(ig,l) > 1.e-10)                             &
+            &           .and.(f_star(ig,l) > 1.e-10)
+         ENDDO
+         
+!-------------------------------------------------------------------------------
+! Latent heat release (before mixing)
+!-------------------------------------------------------------------------------
          
          ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
@@ -188,10 +184,10 @@
          
 !-------------------------------------------------------------------------------
-! Vertical speed (before mixing between plume and environment)
+! Vertical speed (before mixing)
 !-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
             IF (active(ig)) THEN
-               zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsa_est(ig))                ! zqla_est set to ql   plume
+               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)
@@ -229,8 +225,8 @@
                
                zdz = zlev(ig,l+1) - zlev(ig,l)
-               zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.                     ! AB: est-ce la bonne vitesse a utiliser ?
+               zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.
                gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
                
-               IF (zw2_est(ig,l).GT.0.) THEN
+               IF (zw2_est(ig,l) > 0.) THEN
                   test = gamma / zw2m - nu
                ELSE
@@ -241,29 +237,14 @@
                ENDIF
                
-               IF (test.gt.0.) THEN
-                  detr_star(ig,l) = zdz * nu
-                  entr_star(ig,l) = zdz * (zbetalpha * gamma / zw2m + nu)
+               IF (test > 0.) THEN
+                  detr_star(ig,l) = zdz * f_star(ig,l) * nu
+                  entr_star(ig,l) = zdz * f_star(ig,l) * (zbetalpha * gamma / zw2m + nu)
                ELSE
-                  detr_star(ig,l) = zdz * (nu - betalpha * gamma / zw2m)
-                  entr_star(ig,l) = zdz * nu
+                  detr_star(ig,l) = zdz * f_star(ig,l) * (nu - betalpha * gamma / zw2m)
+                  entr_star(ig,l) = zdz * f_star(ig,l) * nu
                ENDIF
                
-!               IF (detr_star(ig,l).lt.0.) THEN
-!                  print *, 'WARNING: detrainment is negative!'
-!                  print *, 'l,detr', l, detr_star(ig,l)
-!               ENDIF
-               
-!               IF (entr_star(ig,l).lt.0.) THEN
-!                  print *, 'WARNING: entrainment is negative!'
-!                  print *, 'l,entr', l, entr_star(ig,l)
-!               ENDIF
-               
                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
-!                  print *, 'WARNING: mass flux is negative!'
-!                  print *, 'l,f_star', l+1, f_star(ig,l+1)
-!               ENDIF
-               
             ENDIF
          ENDDO
@@ -273,5 +254,5 @@
 !-------------------------------------------------------------------------------
          
-         activetmp(:) = active(:) .and. (f_star(:,l+1).GT.1.e-10)
+         activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-10)
          
          DO ig=1,ngrid
@@ -286,4 +267,8 @@
          ENDDO
          
+!-------------------------------------------------------------------------------
+! Latent heat release (after mixing)
+!-------------------------------------------------------------------------------
+         
          ztemp(:) = zpopsk(:,l) * zhla(:,l)
          
@@ -295,12 +280,12 @@
          
 !-------------------------------------------------------------------------------
-! Vertical speed (after mixing between plume and environment)
+! Vertical speed (after mixing)
 !-------------------------------------------------------------------------------
          
          DO ig=1,ngrid
             IF (activetmp(ig)) THEN
-               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)
+               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)                         &  ! ztva is set to TR   plume (mixed)
+               &          + RLvCp * zqla(ig,l)
                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))
@@ -329,4 +314,4 @@
       
       
-RETURN 
+RETURN
 END
