Index: /LMDZ5/trunk/libf/phylmd/cv3_routines.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/cv3_routines.F	(revision 1493)
+++ /LMDZ5/trunk/libf/phylmd/cv3_routines.F	(revision 1494)
@@ -38,7 +38,4 @@
       CHARACTER (LEN=20) :: modname='cv3_param'
       CHARACTER (LEN=80) :: abort_message
-      LOGICAL, SAVE :: FIRST=.TRUE.
-!$OMP THREADPRIVATE(FIRST)
-
 
 c noff: integer limit for convection (nd-noff)
@@ -103,11 +100,9 @@
       CLOSE(99)
 9999  Continue
-      if (first) then
-        WRITE(*,*)'dpbase=',dpbase
-        WRITE(*,*)'pbcrit=',pbcrit
-        WRITE(*,*)'ptcrit=',ptcrit
-        WRITE(*,*)'sigdz=',sigdz
-        WRITE(*,*)'spfac=',spfac
-      endif
+      WRITE(*,*)'dpbase=',dpbase
+      WRITE(*,*)'pbcrit=',pbcrit
+      WRITE(*,*)'ptcrit=',ptcrit
+      WRITE(*,*)'sigdz=',sigdz
+      WRITE(*,*)'spfac=',spfac
 
       return
@@ -1055,5 +1050,5 @@
            pden=ptcrit-pbcrit
            ep(i,k)=(plcl(i)-p(i,k)-pbcrit)/pden*epmax
-           ep(i,k)=amax1(ep(i,k),0.0)
+           ep(i,k)=max(ep(i,k),0.0)
            ep(i,k)=amin1(ep(i,k),epmax)
            sigp(i,k)=spfac
@@ -1157,12 +1152,12 @@
       enddo
 
-      do 530 i=1,ncum
-       do 535 k=1,nl-1
+      do 535 k=1,nl-1
+       do 530 i=1,ncum
         if ((k.ge.iposit(i)).and.(buoy(i,k).lt.dtovsh)) then
          inb(i)=MIN(inb(i),k)
         endif
- 535   continue
- 530  continue
-
+ 530   continue
+ 535  continue
+c
 c-- end convect3
 
@@ -1354,5 +1349,5 @@
         if (k.le.icb(i))then
          sig(i,k)=beta*sig(i,k)-2.*alpha*buoy(i,icb(i))*buoy(i,icb(i))
-         sig(i,k)=amax1(sig(i,k),0.0)
+         sig(i,k)=max(sig(i,k),0.0)
          w0(i,k)=beta*w0(i,k)
         endif
@@ -1364,5 +1359,5 @@
 c!            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
 c!     1              abs(buoy(inb))
-c!            sig(i)=amax1(sig(i),0.0)
+c!            sig(i)=max(sig(i),0.0)
 c!            w0(i)=beta*w0(i)
 c!   85    continue
@@ -1371,5 +1366,5 @@
 c!      do 87 i=1,icb
 c!         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
-c!         sig(i)=amax1(sig(i),0.0)
+c!         sig(i)=max(sig(i),0.0)
 c!         w0(i)=beta*w0(i)
 c!   87 continue
@@ -1436,5 +1431,5 @@
 
          sig(i,k)=beta*sig(i,k)+alpha*dtmin(i,k)*ABS(dtmin(i,k))
-         sig(i,k)=amax1(sig(i,k),0.0)
+         sig(i,k)=max(sig(i,k),0.0)
          sig(i,k)=amin1(sig(i,k),0.01)
          fac=AMIN1(((dtcrit-dtmin(i,k))/dtcrit),1.0)
@@ -1496,5 +1491,5 @@
 c!         dcape=rrd*buoy(i-1)*deltap/p(i-1)
 c!         dlnp=deltap/p(i-1)
-c!         cape=amax1(0.0,cape)
+c!         cape=max(0.0,cape)
 c!         sigold=sig(i)
 
@@ -1505,5 +1500,5 @@
 
 c!         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
-c!         sig(i)=amax1(sig(i),0.0)
+c!         sig(i)=max(sig(i),0.0)
 c!         sig(i)=amin1(sig(i),0.01)
 c!         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
@@ -1654,9 +1649,9 @@
 c!!!      end do
           elij(il,i,j)=altem
-          elij(il,i,j)=amax1(0.0,elij(il,i,j))
+          elij(il,i,j)=max(0.0,elij(il,i,j))
           ment(il,i,j)=m(il,i)/(1.-sij(il,i,j))
           nent(il,i)=nent(il,i)+1
          end if
-         sij(il,i,j)=amax1(0.0,sij(il,i,j))
+         sij(il,i,j)=max(0.0,sij(il,i,j))
          sij(il,i,j)=amin1(1.0,sij(il,i,j))
          endif ! new
@@ -1778,17 +1773,17 @@
         wgh=1.0
         if(j.gt.i)then
-         sjmax=amax1(sij(il,i,j+1),smax(il))
+         sjmax=max(sij(il,i,j+1),smax(il))
          sjmax=amin1(sjmax,scrit(il))
-         smax(il)=amax1(sij(il,i,j),smax(il))
-         sjmin=amax1(sij(il,i,j-1),smax(il))
+         smax(il)=max(sij(il,i,j),smax(il))
+         sjmin=max(sij(il,i,j-1),smax(il))
          sjmin=amin1(sjmin,scrit(il))
          if(sij(il,i,j).lt.(smax(il)-1.0e-16))wgh=0.0
          smid=amin1(sij(il,i,j),scrit(il))
         else
-         sjmax=amax1(sij(il,i,j+1),scrit(il))
-         smid=amax1(sij(il,i,j),scrit(il))
+         sjmax=max(sij(il,i,j+1),scrit(il))
+         smid=max(sij(il,i,j),scrit(il))
          sjmin=0.0
          if(j.gt.1)sjmin=sij(il,i,j-1)
-         sjmin=amax1(sjmin,scrit(il))
+         sjmin=max(sjmin,scrit(il))
         endif
         delp=abs(sjmax-smid)
@@ -1804,5 +1799,5 @@
       do il=1,ncum
        if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
-        asij(il)=amax1(1.0e-16,asij(il))
+        asij(il)=max(1.0e-16,asij(il))
         asij(il)=1.0/asij(il)
         asum(il,i)=0.0
@@ -1834,5 +1829,5 @@
       do il=1,ncum
        if (i.ge.icb(il).and.i.le.inb(il).and.lwork(il)) then
-        bsum(il,i)=amax1(bsum(il,i),1.0e-16)
+        bsum(il,i)=max(bsum(il,i),1.0e-16)
         bsum(il,i)=1.0/bsum(il,i)
        endif
@@ -1951,11 +1946,12 @@
       real tinv, delti
       real awat, afac, afac1, afac2, bfac
-      real pr1, pr2, sigt, b6, c6, revap, tevap, delth
+      real pr1, pr2, sigt, b6, c6, revap, delth
       real amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
       real ampmax
+      real tevap(nloc)
       real lvcp(nloc,na)
       real h(nloc,na),hm(nloc,na)
       real wdtrain(nloc)
-      logical lwork(nloc)
+      logical lwork(nloc),mplus(nloc)
 
 
@@ -1993,9 +1989,8 @@
 
         do il=1,ncum
-          lwork(il)=.TRUE.
-          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
+!!          lwork(il)=.TRUE.
+!!          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
+          lwork(il)= ep(il,inb(il)) .ge. 0.0001
         enddo
-
-        call zilch(wdtrain,ncum)
 
 c   ***  Set the fractionnal area sigd of precipitating downdraughts
@@ -2004,4 +1999,11 @@
         enddo
 
+
+c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c
+c    ***                    begin downdraft loop                    ***
+c
+c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c
         DO 400 i=nl+1,1,-1
 
@@ -2012,13 +2014,10 @@
         if (num1.le.0) goto 400
 
+        call zilch(wdtrain,ncum)
+
 c
 c   ***  integrate liquid water equation to find condensed water   ***
 c   ***                and condensed water flux                    ***
 c
-
-c
-c    ***                    begin downdraft loop                    ***
-c
-
 c
 c    ***              calculate detrained precipitation             ***
@@ -2039,5 +2038,5 @@
           if (i.le.inb(il) .and. lwork(il)) then
            awat=elij(il,j,i)-(1.-ep(il,i))*clw(il,i)
-           awat=amax1(awat,0.0)
+           awat=max(awat,0.0)
            if (cvflag_grav) then
             wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i)
@@ -2055,6 +2054,5 @@
 c
 
-      do 999 il=1,ncum
-
+      do 995 il=1,ncum
        if (i.le.inb(il) .and. lwork(il)) then
 
@@ -2066,5 +2064,5 @@
        rp(il,i)=0.5*(rp(il,i)+rr(il,i))
       endif
-      rp(il,i)=amax1(rp(il,i),0.0)
+      rp(il,i)=max(rp(il,i),0.0)
       rp(il,i)=amin1(rp(il,i),rs(il,i))
       rp(il,inb(il))=rr(il,inb(il))
@@ -2077,5 +2075,5 @@
        rp(il,i-1)=0.5*(rp(il,i-1)+rr(il,i-1))
        rp(il,i-1)=amin1(rp(il,i-1),rs(il,i-1))
-       rp(il,i-1)=amax1(rp(il,i-1),0.0)
+       rp(il,i-1)=max(rp(il,i-1),0.0)
        afac1=p(il,i)*(rs(il,i)-rp(il,i))/(1.0e4+2000.0*p(il,i)*rs(il,i))
        afac2=p(il,i-1)*(rs(il,i-1)-rp(il,i-1))
@@ -2084,5 +2082,5 @@
       endif
       if(i.eq.inb(il))afac=0.0
-      afac=amax1(afac,0.0)
+      afac=max(afac,0.0)
       bfac=1./(sigd(il)*wt(il,i))
 c
@@ -2148,19 +2146,28 @@
      :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
 c
+       endif !(i.le.inb(il) .and. lwork(il))
+995   Continue
+c----------------------------------------------------------------
+c
 ccc
 c    ***  calculate precipitating downdraft mass flux under     ***
 c    ***              hydrostatic approximation                 ***
 c
-      if (i.ne.1) then
-
-      tevap=amax1(0.0,evap(il,i))
-      delth=amax1(0.001,(th(il,i)-th(il,i-1)))
+      Do 996 il = 1,ncum
+       if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then
+c
+      tevap(il)=max(0.0,evap(il,i))
+      delth=max(0.001,(th(il,i)-th(il,i-1)))
       if (cvflag_grav) then
-       mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap
+       mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap(il)
      :              *(p(il,i-1)-p(il,i))/delth
       else
-       mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap
+       mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap(il)
      :         *(p(il,i-1)-p(il,i))/delth
       endif
+c
+       endif !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
+996   Continue
+c----------------------------------------------------------------
 c
 c    ***           if hydrostatic assumption fails,             ***
@@ -2168,7 +2175,11 @@
 c    ***  and mass flux from two simultaneous differential eqns ***
 c
+      Do 997 il = 1,ncum
+       if (i.le.inb(il) .and. lwork(il) .and. i.ne.1) then
+c
       amfac=sigd(il)*sigd(il)*70.0*ph(il,i)*(p(il,i-1)-p(il,i))
      :          *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
       amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
+c
       if(amp2.gt.(0.1*amfac))then
        xf=100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1))
@@ -2177,5 +2188,5 @@
        af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv
        bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
-     :            +50.*(p(il,i-1)-p(il,i))*xf*tevap
+     :            +50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
        fac2=1.0
        if(bf.lt.0.0)fac2=-1.0
@@ -2193,5 +2204,5 @@
         mp(il,i)=mp(il,i+1)*tinv+2.*sqrt(af*tinv)*cos(d*tinv)
        endif
-       mp(il,i)=amax1(0.0,mp(il,i))
+       mp(il,i)=max(0.0,mp(il,i))
 
        if (cvflag_grav) then
@@ -2199,16 +2210,17 @@
 C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
 C Et il faut bien revoir les facteurs 100.
-        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
+        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
      2   /(mp(il,i)+sigd(il)*0.1)
      3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
      : *sigd(il)*th(il,i))
        else
-        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
+        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap(il)
      2   /(mp(il,i)+sigd(il)*0.1)
      3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
      : *sigd(il)*th(il,i))
        endif
-       b(il,i-1)=amax1(b(il,i-1),0.0)
-      endif
+       b(il,i-1)=max(b(il,i-1),0.0)
+c
+      endif !(amp2.gt.(0.1*amfac))
 c
 c   ***         limit magnitude of mp(i) to meet cfl condition      ***
@@ -2216,6 +2228,6 @@
       ampmax=2.0*(ph(il,i)-ph(il,i+1))*delti
       amp2=2.0*(ph(il,i-1)-ph(il,i))*delti
-      ampmax=amin1(ampmax,amp2)
-      mp(il,i)=amin1(mp(il,i),ampmax)
+      ampmax=min(ampmax,amp2)
+      mp(il,i)=min(mp(il,i),ampmax)
 c
 c    ***      force mp to decrease linearly to zero                 ***
@@ -2231,15 +2243,22 @@
       endif
 
-360   continue
-      endif ! i.eq.1
+       endif ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
+997   Continue
+c----------------------------------------------------------------
 c
 c    ***       find mixing ratio of precipitating downdraft     ***
 c
-
-      if (i.ne.inb(il)) then
-
+      Do il = 1,ncum
+       if (i.lt.inb(il) .and. lwork(il)) then
+        mplus(il) = mp(il,i).gt.mp(il,i+1)
+       endif ! (i.lt.inb(il) .and. lwork(il))
+      enddo
+c
+      Do 999 il = 1,ncum
+       if (i.lt.inb(il) .and. lwork(il)) then
+c
       rp(il,i)=rr(il,i)
 
-      if(mp(il,i).gt.mp(il,i+1))then
+      if(mplus(il))then
 
        if (cvflag_grav) then
@@ -2258,11 +2277,5 @@
       vp(il,i)=vp(il,i)/mp(il,i)
 
-      do j=1,ntra
-      trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
-     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
-      trap(il,i,j)=trap(il,i,j)/mp(il,i)
-      end do
-
-      else
+      else ! if (mplus(il))
 
        if(mp(il,i+1).gt.1.0e-16)then
@@ -2278,19 +2291,41 @@
        up(il,i)=up(il,i+1)
        vp(il,i)=vp(il,i+1)
-
-       do j=1,ntra
-       trap(il,i,j)=trap(il,i+1,j)
-       end do
-
-       endif
-      endif
+       endif ! (mp(il,i+1).gt.1.0e-16)
+      endif ! (mplus(il)) else if (.not.mplus(il))
+c
       rp(il,i)=amin1(rp(il,i),rs(il,i))
-      rp(il,i)=amax1(rp(il,i),0.0)
-
-      endif
-      endif
+      rp(il,i)=max(rp(il,i),0.0)
+
+      endif ! (i.lt.inb(il) .and. lwork(il))
 999   continue
+c----------------------------------------------------------------
+c
+c    ***       find tracer concentrations in precipitating downdraft     ***
+c
+      do j=1,ntra
+       do il = 1,ncum
+       if (i.lt.inb(il) .and. lwork(il)) then
+c
+         if(mplus(il))then
+          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
+     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
+          trap(il,i,j)=trap(il,i,j)/mp(il,i)
+         else ! if (mplus(il))
+          if(mp(il,i+1).gt.1.0e-16)then
+           trap(il,i,j)=trap(il,i+1,j)
+          endif
+         endif ! (mplus(il)) else if (.not.mplus(il))
+c
+        endif ! (i.lt.inb(il) .and. lwork(il))
+       enddo
+      end do
 
 400   continue
+c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c
+c    ***                    end of downdraft loop                    ***
+c
+c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+c
 
        return
@@ -2814,5 +2849,5 @@
        do il = 1,ncum
         awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
-        awat(il)=amax1(awat(il),0.0)
+        awat(il)=max(awat(il),0.0)
        enddo
 c
Index: /LMDZ5/trunk/libf/phylmd/cv3p_mixing.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/cv3p_mixing.F	(revision 1493)
+++ /LMDZ5/trunk/libf/phylmd/cv3p_mixing.F	(revision 1494)
@@ -53,6 +53,5 @@
       real Qmixmin(nloc), Rmixmin(nloc), SQmRmin(nloc)
       real signhpmh(nloc)
-      real Sx, Scrit2
-      integer Jx
+      real Sx(nloc), Scrit2
       real smid(nloc), sjmin(nloc), sjmax(nloc)
       real Sbef(nloc), Sup(nloc), Smin(nloc)
@@ -240,5 +239,7 @@
       enddo
 
-      DO 789 i=minorig+1,nl
+c---------------------------------------------------------------
+      DO 789 i=minorig+1,nl     !Loop on origin level "i"
+c---------------------------------------------------------------
 
       num1=0
@@ -248,9 +249,32 @@
       if (num1.le.0) goto 789
 
+c
+cjyg1    Find maximum of SIJ for J>I, if any.
+c
+       Sx(:) = 0.
+c
+       DO il = 1,ncum
+        IF ( i.ge.icb(il) .and. i.le.inb(il) ) THEN
+         Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
+         Sbef(il) = max(0.,signhpmh(il))
+        ENDIF
+       ENDDO
+
+       DO j = i+1,nl
+        DO il = 1,ncum
+         IF ( i.ge.icb(il) .and. i.le.inb(il)
+     :         .and. j.le.inb(il) ) THEN
+          IF (Sbef(il) .LT. Sij(il,i,j)) THEN
+            Sx(il) = max(Sij(il,i,j),Sx(il))
+          ENDIF
+          Sbef(il) = Sij(il,i,j)
+         ENDIF
+        ENDDO
+       ENDDO
+c
 
       do 781 il=1,ncum
        if ( i.ge.icb(il) .and. i.le.inb(il) ) then
         lwork(il)=(nent(il,i).ne.0)
-        Signhpmh(il) = sign(1.,hp(il,i)-h(il,i))
         qp=qnk(il)-ep(il,i)*clw(il,i)
         anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
@@ -262,20 +286,9 @@
         alt=qp-rs(il,i)+scrit(il)*(rr(il,i)-qp)
 c
-cjyg1    Find maximum of SIJ for J>I, if any, and new critical value Scrit2
+cjyg1    Find new critical value Scrit2
 c        such that : Sij > Scrit2  => mixed draught will detrain at J<I
 c                    Sij < Scrit2  => mixed draught will detrain at J>I
 c
-       Sx = 0.
-       Jx = 0.
-       Sbef(il) = max(0.,signhpmh(il))
-       DO j = i+1,inb(il)
-         IF (Sbef(il) .LT. Sij(il,i,j)) THEN
-           Sx = max(Sij(il,i,j),Sx)
-           Jx = J
-         ENDIF
-         Sbef(il) = Sij(il,i,j)
-       ENDDO
-c
-       Scrit2 = min(Scrit(il),Sx)*max(0.,-signhpmh(il))
+       Scrit2 = min(Scrit(il),Sx(il))*max(0.,-signhpmh(il))
      :         +Scrit(il)*max(0.,signhpmh(il))
 c
@@ -293,5 +306,7 @@
 781   continue
 
-      do 175 j=minorig,nl
+c---------------------------------------------------------------
+      do 175 j=minorig,nl   !Loop on destination level "j"
+c---------------------------------------------------------------
 
       num2=0
Index: /LMDZ5/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 1493)
+++ /LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 1494)
@@ -80,5 +80,5 @@
 !$OMP THREADPRIVATE(lev_out)
 
-      INTEGER ig,k,l,ll
+      INTEGER ig,k,l,ll,ierr
       real zsortie1d(klon)
       INTEGER lmax(klon),lmin(klon),lalim(klon)
@@ -233,7 +233,4 @@
 !     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
      do ig=1,klon
-      if (prt_level.ge.20) then
-       print*,'th_main ig f0',ig,f0(ig)
-      endif
          f0(ig)=max(f0(ig),1.e-2)
          zmax0(ig)=max(zmax0(ig),40.)
@@ -241,4 +238,9 @@
      enddo
 
+      if (prt_level.ge.20) then
+       do ig=1,ngrid
+          print*,'th_main ig f0',ig,f0(ig)
+       enddo
+      endif
 !-----------------------------------------------------------------------
 ! Calcul de T,q,ql a partir de Tl et qT dans l environnement
@@ -290,9 +292,6 @@
 !-----------------------------------------------------------------------
 
-      do l=1,nlay
-         rho(:,l)=pplay(:,l)/(zpspsk(:,l)*RD*ztv(:,l))
-      enddo
-
-!IM
+     rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
+
      if (prt_level.ge.10)write(lunout,*)                                &
     &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
@@ -619,11 +618,16 @@
       enddo
 !IM
+      ierr=0
       do ig=1,ngrid
         if (pcon(ig).le.pplay(ig,nlay)) then 
            zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
+           ierr=1
+        endif
+      enddo
+      if (ierr==1) then
            abort_message = 'thermcellV0_main: les thermiques vont trop haut '
            CALL abort_gcm (modname,abort_message,1)
-        endif
-      enddo
+      endif
+
       if (prt_level.ge.1) print*,'14b OK convect8'
       do k=nlay,1,-1
@@ -655,7 +659,4 @@
             zf2=zf/(1.-zf)
 !
-      if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
-!
-      if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
             thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
             if(zw2(ig,l).gt.1.e-10) then
@@ -664,8 +665,6 @@
              wth2(ig,l)=0.
             endif
-!           print*,'wth2=',wth2(ig,l)
             wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
      &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
-      if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
             q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
 !test: on calcul q2/po=ratqsc
@@ -682,4 +681,12 @@
       enddo
 !
+      if (prt_level.ge.10) then
+         ig=igout
+         do l=1,nlay
+            print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
+            print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
+         enddo
+      endif
+
 !      print*,'avant calcul ale et alp' 
 !calcul de ALE et ALP pour la convection
@@ -705,36 +712,25 @@
 !initialisations
 !      print*,'ponderation'
-      do ig=1,ngrid
-           fm_tot(ig)=0.
-      enddo
-       do ig=1,ngrid
-        do k=1,klev
-           wght_th(ig,k)=1.
-        enddo
-       enddo
-       do ig=1,ngrid
-!         lalim_conv(ig)=lmix_bis(ig)
-!la hauteur de la couche alim_conv = hauteur couche alim_therm
-         lalim_conv(ig)=lalim(ig)
-!         zentr(ig)=zlev(ig,lalim(ig))
-      enddo
-      do ig=1,ngrid
-        do k=1,lalim_conv(ig)
-           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
-        enddo
-      enddo
-      do ig=1,ngrid
-        do k=1,lalim_conv(ig)
-           if (fm_tot(ig).gt.1.e-10) then
-!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
-           endif
-!on pondere chaque couche par a*
-             if (alim_star(ig,k).gt.1.e-10) then
-             wght_th(ig,k)=alim_star(ig,k)
-             else
-             wght_th(ig,k)=1.
-             endif
-        enddo
-      enddo
+
+      fm_tot(:)=0.
+      wght_th(:,:)=1.
+      lalim_conv(:)=lalim(:)
+
+      do k=1,klev
+         do ig=1,ngrid
+            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
+         enddo
+      enddo
+
+! assez bizarre car, si on est dans la couche d'alim et que alim_star et
+! plus petit que 1.e-10, on prend wght_th=1.
+      do k=1,klev
+         do ig=1,ngrid
+            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
+               wght_th(ig,k)=alim_star(ig,k)
+            endif
+         enddo
+      enddo
+
 !      print*,'apres wght_th'
 !test pour prolonger la convection
@@ -748,4 +744,5 @@
       enddo
 
+
 !calcul du ratqscdiff
       if (prt_level.ge.1) print*,'14e OK convect8'
@@ -753,20 +750,25 @@
       vardiff=0.
       ratqsdiff(:,:)=0.
-      do ig=1,ngrid
-         do l=1,lalim(ig)
+
+      do l=1,klev
+         do ig=1,ngrid
+            if (l<=lalim(ig)) then
             var=var+alim_star(ig,l)*zqta(ig,l)*1000.
-         enddo
-      enddo
+            endif
+         enddo
+      enddo
+
       if (prt_level.ge.1) print*,'14f OK convect8'
-      do ig=1,ngrid
-          do l=1,lalim(ig)
-          zf=fraca(ig,l)
-          zf2=zf/(1.-zf)
-          vardiff=vardiff+alim_star(ig,l)  &
-     &           *(zqta(ig,l)*1000.-var)**2
-!         ratqsdiff=ratqsdiff+alim_star(ig,l)*
-!     s          (zqta(ig,l)*1000.-po(ig,l)*1000.)**2
-          enddo
-      enddo
+
+      do l=1,klev
+         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
@@ -779,5 +781,5 @@
 !
 !ecriture des fichiers sortie
-!     print*,'15 OK convect8'
+!     print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
 
 #ifdef wrgrads_thermcell
