Index: LMDZ5/trunk/libf/phylmd/cv3_routines.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cv3_routines.F90	(revision 1992)
+++ LMDZ5/trunk/libf/phylmd/cv3_routines.F90	(revision 2007)
@@ -7,26 +7,26 @@
   IMPLICIT NONE
 
-  ! ------------------------------------------------------------
-  ! Set parameters for convectL for iflag_con = 3
-  ! ------------------------------------------------------------
-
-
-  ! ***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
-  ! ***      PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO     ***
-  ! ***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
-  ! ***            EFFICIENCY IS ASSUMED TO BE UNITY            ***
-  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
-  ! ***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
-  ! ***                        OF CLOUD                         ***
-
-  ! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
-  ! ***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
-  ! ***                 APPROACH TO QUASI-EQUILIBRIUM           ***
-  ! ***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
-  ! ***           (BETA MUST BE LESS THAN OR EQUAL TO 1)        ***
-
-  ! ***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
-  ! ***                 APPROACH TO QUASI-EQUILIBRIUM           ***
-  ! ***                     IT MUST BE LESS THAN 0              ***
+!------------------------------------------------------------
+!Set parameters for convectL for iflag_con = 3
+!------------------------------------------------------------
+
+
+!***  PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE ***
+!***      PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO     ***
+!***  PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. ***
+!***            EFFICIENCY IS ASSUMED TO BE UNITY            ***
+!***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
+!***  SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE ***
+!***                        OF CLOUD                         ***
+
+![TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA]
+!***    ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF ***
+!***                 APPROACH TO QUASI-EQUILIBRIUM           ***
+!***    (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) ***
+!***           (BETA MUST BE LESS THAN OR EQUAL TO 1)        ***
+
+!***    DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE ***
+!***                 APPROACH TO QUASI-EQUILIBRIUM           ***
+!***                     IT MUST BE LESS THAN 0              ***
 
   include "cv3param.h"
@@ -41,10 +41,10 @@
 
   LOGICAL, SAVE :: first = .TRUE.
-  !$OMP THREADPRIVATE(first)
-
-  ! noff: integer limit for convection (nd-noff)
-  ! minorig: First level of convection
-
-  ! -- limit levels for convection:
+!$OMP THREADPRIVATE(first)
+
+!glb  noff: integer limit for convection (nd-noff)
+! minorig: First level of convection
+
+! -- limit levels for convection:
 
   noff = 1
@@ -56,41 +56,41 @@
   IF (first) THEN
 
-    ! -- "microphysical" parameters:
+!$OMP MASTER
+! -- "microphysical" parameters:
     sigdz = 0.01
     spfac = 0.15
     pbcrit = 150.0
     ptcrit = 500.0
-    ! IM beg: ajout fis. reglage ep
+! IM beg: ajout fis. reglage ep
     flag_epkeorig = 1
     elcrit = 0.0003
     tlcrit = -55.0
-    ! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
+! IM lu dans physiq.def via conf_phys.F90     epmax  = 0.993
 
     omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
 
-    ! -- misc:
+! -- misc:
 
     dtovsh = -0.2 ! dT for overshoot
     dpbase = -40. ! definition cloud base (400m above LCL)
-    ! cc      dttrig = 5.   ! (loose) condition for triggering
+! cc      dttrig = 5.   ! (loose) condition for triggering
     dttrig = 10. ! (loose) condition for triggering
     flag_wb = 1
     wbmax = 6. ! (m/s) adiab updraught speed at LFC (used in cv3p1_closure)
 
-    ! -- rate of approach to quasi-equilibrium:
+! -- rate of approach to quasi-equilibrium:
 
     dtcrit = -2.0
     tau = 8000.
 
-    ! -- interface cloud parameterization:
+! -- interface cloud parameterization:
 
     delta = 0.01 ! cld
 
-    ! -- interface with boundary-layer (gust factor): (sb)
+! -- interface with boundary-layer (gust factor): (sb)
 
     betad = 10.0 ! original value (from convect 4.3)
 
-    OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', &
-      ERR=9999)
+    OPEN (99, FILE='conv_param.data', STATUS='old', FORM='formatted', ERR=9999)
     READ (99, *, END=9998) dpbase
     READ (99, *, END=9998) pbcrit
@@ -113,5 +113,5 @@
     WRITE (*, *) 'wbmax =', wbmax
 
-    ! IM Lecture du fichier ep_param.data
+! IM Lecture du fichier ep_param.data
     OPEN (79, FILE='ep_param.data', STATUS='old', FORM='formatted', ERR=7999)
     READ (79, *, END=7998) flag_epkeorig
@@ -124,41 +124,43 @@
     WRITE (*, *) 'elcrit=', elcrit
     WRITE (*, *) 'tlcrit=', tlcrit
-    ! IM end: ajout fis. reglage ep
+! IM end: ajout fis. reglage ep
 
     first = .FALSE.
-  END IF
+!$OMP END MASTER
+
+  END IF ! (first)
 
   beta = 1.0 - delt/tau
   alpha1 = 1.5E-3
-  ! jyg    Correction bug alpha
+!JYG    Correction bug alpha
   alpha1 = alpha1*1.5
   alpha = alpha1*delt/tau
-  ! jyg    Bug
-  ! cc increase alpha to compensate W decrease:
-  ! c      alpha  = alpha*1.5
+!JYG    Bug
+! cc increase alpha to compensate W decrease:
+! c      alpha  = alpha*1.5
 
   RETURN
 END SUBROUTINE cv3_param
 
-SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, lv, lf, cpn, tv, gz, h, hm, &
-    th)
+SUBROUTINE cv3_prelim(len, nd, ndp1, t, q, p, ph, &
+                      lv, lf, cpn, tv, gz, h, hm, th)
   IMPLICIT NONE
 
-  ! =====================================================================
-  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
-  ! "ori": from convect4.3 (vectorized)
-  ! "convect3": to be exactly consistent with convect3
-  ! =====================================================================
-
-  ! inputs:
+! =====================================================================
+! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
+! "ori": from convect4.3 (vectorized)
+! "convect3": to be exactly consistent with convect3
+! =====================================================================
+
+! inputs:
   INTEGER len, nd, ndp1
   REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
 
-  ! outputs:
+! outputs:
   REAL lv(len, nd), lf(len, nd), cpn(len, nd), tv(len, nd)
   REAL gz(len, nd), h(len, nd), hm(len, nd)
   REAL th(len, nd)
 
-  ! local variables:
+! local variables:
   INTEGER k, i
   REAL rdcp
@@ -170,15 +172,15 @@
 
 
-  ! ori      do 110 k=1,nlp
-  ! abderr     do 110 k=1,nl ! convect3
+! ori      do 110 k=1,nlp
+! abderr     do 110 k=1,nl ! convect3
   DO k = 1, nlp
 
     DO i = 1, len
-      ! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
+! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
       lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
       lf(i, k) = lf0 - clmci*(t(i,k)-273.15)
       cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
       cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
-      ! ori          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
+! ori          tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1)
       tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k))
       rdcp = (rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i, k)
@@ -187,28 +189,28 @@
   END DO
 
-  ! gz = phi at the full levels (same as p).
+! gz = phi at the full levels (same as p).
 
   DO i = 1, len
     gz(i, 1) = 0.0
   END DO
-  ! ori      do 140 k=2,nlp
+! ori      do 140 k=2,nlp
   DO k = 2, nl ! convect3
     DO i = 1, len
-      tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3
-      tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3
-      gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy) & !convect3
-        *(p(i,k-1)-p(i,k))/ph(i, k) !convect3
-
-      ! c        print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
-
-      ! ori         gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
-      ! ori    &         *(p(i,k-1)-p(i,k))/ph(i,k)
-    END DO
-  END DO
-
-  ! h  = phi + cpT (dry static energy).
-  ! hm = phi + cp(T-Tbase)+Lq
-
-  ! ori      do 170 k=1,nlp
+      tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k))         !convect3
+      tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1))   !convect3
+      gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy)* & !convect3
+                 (p(i,k-1)-p(i,k))/ph(i, k)        !convect3
+
+! c        print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
+
+! ori         gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
+! ori    &         *(p(i,k-1)-p(i,k))/ph(i,k)
+    END DO
+  END DO
+
+! h  = phi + cpT (dry static energy).
+! hm = phi + cp(T-Tbase)+Lq
+
+! ori      do 170 k=1,nlp
   DO k = 1, nl ! convect3
     DO i = 1, len
@@ -221,28 +223,31 @@
 END SUBROUTINE cv3_prelim
 
-SUBROUTINE cv3_feed(len, nd, t, q, u, v, p, ph, hm, gz, p1feed, p2feed, wght, &
-    wghti, tnk, thnk, qnk, qsnk, unk, vnk, cpnk, hnk, nk, icb, icbmax, iflag, &
-    gznk, plcl)
+SUBROUTINE cv3_feed(len, nd, ok_conserv_q, &
+                    t, q, u, v, p, ph, hm, gz, &
+                    p1feed, p2feed, wght, &
+                    wghti, tnk, thnk, qnk, qsnk, unk, vnk, &
+                    cpnk, hnk, nk, icb, icbmax, iflag, gznk, plcl)
   IMPLICIT NONE
 
-  ! ================================================================
-  ! Purpose: CONVECTIVE FEED
-
-  ! Main differences with cv_feed:
-  ! - ph added in input
-  ! - here, nk(i)=minorig
-  ! - icb defined differently (plcl compared with ph instead of p)
-
-  ! Main differences with convect3:
-  ! - we do not compute dplcldt and dplcldr of CLIFT anymore
-  ! - values iflag different (but tests identical)
-  ! - A,B explicitely defined (!...)
-  ! ================================================================
+! ================================================================
+! Purpose: CONVECTIVE FEED
+
+! Main differences with cv_feed:
+! - ph added in input
+! - here, nk(i)=minorig
+! - icb defined differently (plcl compared with ph instead of p)
+
+! Main differences with convect3:
+! - we do not compute dplcldt and dplcldr of CLIFT anymore
+! - values iflag different (but tests identical)
+! - A,B explicitely defined (!...)
+! ================================================================
 
   include "cv3param.h"
   include "cvthermo.h"
 
-  ! inputs:
+!inputs:
   INTEGER len, nd
+  LOGICAL ok_conserv_q
   REAL t(len, nd), q(len, nd), p(len, nd)
   REAL u(len, nd), v(len, nd)
@@ -250,11 +255,11 @@
   REAL ph(len, nd+1)
   REAL p1feed(len)
-  ! ,  wght(len)
+! ,  wght(len)
   REAL wght(nd)
-  ! input-output
+!input-output
   REAL p2feed(len)
-  ! outputs:
+!outputs:
   INTEGER iflag(len), nk(len), icb(len), icbmax
-  ! real   wghti(len)
+!   real   wghti(len)
   REAL wghti(len, nd)
   REAL tnk(len), thnk(len), qnk(len), qsnk(len)
@@ -263,5 +268,5 @@
   REAL plcl(len)
 
-  ! local variables:
+!local variables:
   INTEGER i, k, iter, niter
   INTEGER ihmin(len)
@@ -269,10 +274,34 @@
   REAL pup(len), plo(len), pfeed(len)
   REAL plclup(len), plcllo(len), plclfeed(len)
+  REAL pfeedmin(len)
   REAL posit(len)
   LOGICAL nocond(len)
 
-  ! -------------------------------------------------------------------
-  ! --- Origin level of ascending parcels for convect3:
-  ! -------------------------------------------------------------------
+!jyg20140217<
+  INTEGER iostat
+  LOGICAL, SAVE :: first
+  LOGICAL, SAVE :: ok_new_feed
+  REAL, SAVE :: dp_lcl_feed
+!$OMP THREADPRIVATE (first,ok_new_feed,dp_lcl_feed)
+  DATA first/.TRUE./
+  DATA dp_lcl_feed/2./
+
+  IF (first) THEN
+!$OMP MASTER
+    ok_new_feed = ok_conserv_q
+    OPEN (98, FILE='cv3feed_param.data', STATUS='old', FORM='formatted', IOSTAT=iostat)
+    IF (iostat==0) THEN
+      READ (98, *, END=998) ok_new_feed
+998   CONTINUE
+      CLOSE (98)
+    END IF
+    PRINT *, ' ok_new_feed: ', ok_new_feed
+    first = .FALSE.
+!$OMP END MASTER
+  END IF
+!jyg>
+! -------------------------------------------------------------------
+! --- Origin level of ascending parcels for convect3:
+! -------------------------------------------------------------------
 
   DO i = 1, len
@@ -281,26 +310,28 @@
   END DO
 
-  ! -------------------------------------------------------------------
-  ! --- Adjust feeding layer thickness so that lifting up to the top of
-  ! --- the feeding layer does not induce condensation (i.e. so that
-  ! --- plcl < p2feed).
-  ! --- Method : iterative secant method.
-  ! -------------------------------------------------------------------
-
-  ! 1- First bracketing of the solution : ph(nk+1), p2feed
-
-  ! 1.a- LCL associated to p2feed
+! -------------------------------------------------------------------
+! --- Adjust feeding layer thickness so that lifting up to the top of
+! --- the feeding layer does not induce condensation (i.e. so that
+! --- plcl < p2feed).
+! --- Method : iterative secant method.
+! -------------------------------------------------------------------
+
+! 1- First bracketing of the solution : ph(nk+1), p2feed
+
+! 1.a- LCL associated with p2feed
   DO i = 1, len
     pup(i) = p2feed(i)
   END DO
-  CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, t, q, u, v, wght, &
-    wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
-  ! 1.b- LCL associated to ph(nk+1)
+  CALL cv3_vertmix(len, nd, iflag, p1feed, pup, p, ph, &
+                   t, q, u, v, wght, &
+                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclup)
+! 1.b- LCL associated with ph(nk+1)
   DO i = 1, len
     plo(i) = ph(i, nk(i)+1)
   END DO
-  CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, t, q, u, v, wght, &
-    wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
-  ! 2- Iterations
+  CALL cv3_vertmix(len, nd, iflag, p1feed, plo, p, ph, &
+                   t, q, u, v, wght, &
+                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plcllo)
+! 2- Iterations
   niter = 5
   DO iter = 1, niter
@@ -314,17 +345,55 @@
         pfeed(i) = pup(i)
       ELSE
-        pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+plo(i)*(plclup(i)-pup(i)))/ &
-          (plo(i)-plcllo(i)+plclup(i)-pup(i))
+!JYG20140217<
+        IF (ok_new_feed) THEN
+          pfeed(i) = (pup(i)*(plo(i)-plcllo(i)-dp_lcl_feed)+  &
+                      plo(i)*(plclup(i)-pup(i)+dp_lcl_feed))/ &
+                     (plo(i)-plcllo(i)+plclup(i)-pup(i))
+        ELSE
+          pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+  &
+                      plo(i)*(plclup(i)-pup(i)))/ &
+                     (plo(i)-plcllo(i)+plclup(i)-pup(i))
+        END IF
+!JYG>
       END IF
     END DO
-    CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, t, q, u, v, wght, &
-      wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
+!jyg20140217<
+! For the last iteration, make sure that the top of the feeding layer
+! and LCL are not in the same layer:
+    IF (ok_new_feed) THEN
+      IF (iter==niter) THEN
+        DO k = minorig, nd
+          DO i = 1, len
+            IF (ph(i,k)>=plclfeed(i)) pfeedmin(i) = ph(i, k)
+          END DO
+        END DO
+        DO i = 1, len
+          pfeed(i) = max(pfeedmin(i), pfeed(i))
+        END DO
+      END IF
+    END IF
+!jyg>
+
+    CALL cv3_vertmix(len, nd, iflag, p1feed, pfeed, p, ph, &
+                   t, q, u, v, wght, &
+                   wghti, nk, tnk, thnk, qnk, qsnk, unk, vnk, plclfeed)
+!jyg20140217<
+    IF (ok_new_feed) THEN
+      DO i = 1, len
+        posit(i) = (sign(1.,plclfeed(i)-pfeed(i)+dp_lcl_feed)+1.)*0.5
+        IF (plclfeed(i)-pfeed(i)+dp_lcl_feed==0.) posit(i) = 1.
+      END DO
+    ELSE
+      DO i = 1, len
+        posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
+        IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
+      END DO
+    END IF
+!jyg>
     DO i = 1, len
-      posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
-      IF (plclfeed(i)==pfeed(i)) posit(i) = 1.
-      ! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
-      ! -               => pup=pfeed
-      ! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
-      ! -               => plo=pfeed
+! - posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
+! -               => pup=pfeed
+! - posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
+! -               => plo=pfeed
       pup(i) = posit(i)*pfeed(i) + (1.-posit(i))*pup(i)
       plo(i) = (1.-posit(i))*pfeed(i) + posit(i)*plo(i)
@@ -343,30 +412,30 @@
   END DO
 
-  ! -------------------------------------------------------------------
-  ! --- Check whether parcel level temperature and specific humidity
-  ! --- are reasonable
-  ! -------------------------------------------------------------------
+! -------------------------------------------------------------------
+! --- Check whether parcel level temperature and specific humidity
+! --- are reasonable
+! -------------------------------------------------------------------
   DO i = 1, len
     IF (((tnk(i)<250.0) .OR. (qnk(i)<=0.0)) .AND. (iflag(i)==0)) iflag(i) = 7
   END DO
 
-  ! -------------------------------------------------------------------
-  ! --- Calculate first level above lcl (=icb)
-  ! -------------------------------------------------------------------
-
-  ! @      do 270 i=1,len
-  ! @       icb(i)=nlm
-  ! @ 270  continue
-  ! @c
-  ! @      do 290 k=minorig,nl
-  ! @        do 280 i=1,len
-  ! @          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
-  ! @     &    icb(i)=min(icb(i),k)
-  ! @ 280    continue
-  ! @ 290  continue
-  ! @c
-  ! @      do 300 i=1,len
-  ! @        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
-  ! @ 300  continue
+! -------------------------------------------------------------------
+! --- Calculate first level above lcl (=icb)
+! -------------------------------------------------------------------
+
+!@      do 270 i=1,len
+!@       icb(i)=nlm
+!@ 270  continue
+!@c
+!@      do 290 k=minorig,nl
+!@        do 280 i=1,len
+!@          if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i)))
+!@     &    icb(i)=min(icb(i),k)
+!@ 280    continue
+!@ 290  continue
+!@c
+!@      do 300 i=1,len
+!@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
+!@ 300  continue
 
   DO i = 1, len
@@ -374,7 +443,7 @@
   END DO
 
-  ! la modification consiste a comparer plcl a ph et non a p:
-  ! icb est defini par :  ph(icb)<plcl<ph(icb-1)
-  ! @      do 290 k=minorig,nl
+! la modification consiste a comparer plcl a ph et non a p:
+! icb est defini par :  ph(icb)<plcl<ph(icb-1)
+!@      do 290 k=minorig,nl
   DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2
     DO i = 1, len
@@ -384,10 +453,10 @@
 
 
-  ! print*,'icb dans cv3_feed '
-  ! write(*,'(64i2)') icb(2:len-1)
-  ! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
+! print*,'icb dans cv3_feed '
+! write(*,'(64i2)') icb(2:len-1)
+! call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
 
   DO i = 1, len
-    ! @        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
+!@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
     IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9
   END DO
@@ -397,10 +466,10 @@
   END DO
 
-  ! Compute icbmax.
+! Compute icbmax.
 
   icbmax = 2
   DO i = 1, len
-    ! !        icbmax=max(icbmax,icb(i))
-    IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02
+!!        icbmax=max(icbmax,icb(i))
+    IF (iflag(i)<7) icbmax = max(icbmax, icb(i))     ! sb Jun7th02
   END DO
 
@@ -409,24 +478,24 @@
 
 SUBROUTINE cv3_undilute1(len, nd, t, qs, gz, plcl, p, icb, tnk, qnk, gznk, &
-    tp, tvp, clw, icbs)
+                         tp, tvp, clw, icbs)
   IMPLICIT NONE
 
-  ! ----------------------------------------------------------------
-  ! Equivalent de TLIFT entre NK et ICB+1 inclus
-
-  ! Differences with convect4:
-  ! - specify plcl in input
-  ! - icbs is the first level above LCL (may differ from icb)
-  ! - in the iterations, used x(icbs) instead x(icb)
-  ! - many minor differences in the iterations
-  ! - tvp is computed in only one time
-  ! - icbs: first level above Plcl (IMIN de TLIFT) in output
-  ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
-  ! ----------------------------------------------------------------
+! ----------------------------------------------------------------
+! Equivalent de TLIFT entre NK et ICB+1 inclus
+
+! Differences with convect4:
+!    - specify plcl in input
+!    - icbs is the first level above LCL (may differ from icb)
+!    - in the iterations, used x(icbs) instead x(icb)
+!    - many minor differences in the iterations
+!    - tvp is computed in only one time
+!    - icbs: first level above Plcl (IMIN de TLIFT) in output
+!    - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1)
+! ----------------------------------------------------------------
 
   include "cvthermo.h"
   include "cv3param.h"
 
-  ! inputs:
+! inputs:
   INTEGER len, nd
   INTEGER icb(len)
@@ -436,8 +505,8 @@
   REAL plcl(len) ! convect3
 
-  ! outputs:
+! outputs:
   REAL tp(len, nd), tvp(len, nd), clw(len, nd)
 
-  ! local variables:
+! local variables:
   INTEGER i, k
   INTEGER icb1(len), icbs(len), icbsmax2 ! convect3
@@ -448,160 +517,159 @@
   REAL cpinv(len) ! convect3
 
-  ! -------------------------------------------------------------------
-  ! --- Calculates the lifted parcel virtual temperature at nk,
-  ! --- the actual temperature, and the adiabatic
-  ! --- liquid water content. The procedure is to solve the equation.
-  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
-  ! -------------------------------------------------------------------
-
-
-  ! ***  Calculate certain parcel quantities, including static energy   ***
+! -------------------------------------------------------------------
+! --- Calculates the lifted parcel virtual temperature at nk,
+! --- the actual temperature, and the adiabatic
+! --- liquid water content. The procedure is to solve the equation.
+!     cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
+! -------------------------------------------------------------------
+
+
+! ***  Calculate certain parcel quantities, including static energy   ***
 
   DO i = 1, len
-    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
-      273.15)) + gznk(i)
+    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
     cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
     cpinv(i) = 1./cpp(i)
   END DO
 
-  ! ***   Calculate lifted parcel quantities below cloud base   ***
+! ***   Calculate lifted parcel quantities below cloud base   ***
+
+  DO i = 1, len                                           !convect3
+    icb1(i) = max(icb(i), 2)                              !convect3
+    icb1(i) = min(icb(i), nl)                             !convect3
+! if icb is below LCL, start loop at ICB+1:
+! (icbs est le premier niveau au-dessus du LCL)
+    icbs(i) = icb1(i)                                     !convect3
+    IF (plcl(i)<p(i,icb1(i))) THEN
+      icbs(i) = min(icbs(i)+1, nl)                        !convect3
+    END IF
+  END DO                                                  !convect3
 
   DO i = 1, len !convect3
-    icb1(i) = max(icb(i), 2) !convect3
-    icb1(i) = min(icb(i), nl) !convect3
-    ! if icb is below LCL, start loop at ICB+1:
-    ! (icbs est le premier niveau au-dessus du LCL)
-    icbs(i) = icb1(i) !convect3
-    IF (plcl(i)<p(i,icb1(i))) THEN
-      icbs(i) = min(icbs(i)+1, nl) !convect3
-    END IF
+    ticb(i) = t(i, icbs(i))                               !convect3
+    gzicb(i) = gz(i, icbs(i))                             !convect3
+    qsicb(i) = qs(i, icbs(i))                             !convect3
   END DO !convect3
 
-  DO i = 1, len !convect3
-    ticb(i) = t(i, icbs(i)) !convect3
-    gzicb(i) = gz(i, icbs(i)) !convect3
-    qsicb(i) = qs(i, icbs(i)) !convect3
-  END DO !convect3
-
-
-  ! Re-compute icbsmax (icbsmax2):        !convect3
-  ! !convect3
-  icbsmax2 = 2 !convect3
-  DO i = 1, len !convect3
-    icbsmax2 = max(icbsmax2, icbs(i)) !convect3
-  END DO !convect3
-
-  ! initialization outputs:
-
-  DO k = 1, icbsmax2 ! convect3
-    DO i = 1, len ! convect3
-      tp(i, k) = 0.0 ! convect3
-      tvp(i, k) = 0.0 ! convect3
-      clw(i, k) = 0.0 ! convect3
-    END DO ! convect3
-  END DO ! convect3
-
-  ! tp and tvp below cloud base:
+
+! Re-compute icbsmax (icbsmax2):                          !convect3
+!                                                         !convect3
+  icbsmax2 = 2                                            !convect3
+  DO i = 1, len                                           !convect3
+    icbsmax2 = max(icbsmax2, icbs(i))                     !convect3
+  END DO                                                  !convect3
+
+! initialization outputs:
+
+  DO k = 1, icbsmax2                                      ! convect3
+    DO i = 1, len                                         ! convect3
+      tp(i, k) = 0.0                                      ! convect3
+      tvp(i, k) = 0.0                                     ! convect3
+      clw(i, k) = 0.0                                     ! convect3
+    END DO                                                ! convect3
+  END DO                                                  ! convect3
+
+! tp and tvp below cloud base:
 
   DO k = minorig, icbsmax2 - 1
     DO i = 1, len
       tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i)
-      tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3)
-    END DO
-  END DO
-
-  ! ***  Find lifted parcel quantities above cloud base    ***
+      tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i))        !whole thing (convect3)
+    END DO
+  END DO
+
+! ***  Find lifted parcel quantities above cloud base    ***
 
   DO i = 1, len
     tg = ticb(i)
-    ! ori         qg=qs(i,icb(i))
+! ori         qg=qs(i,icb(i))
     qg = qsicb(i) ! convect3
-    ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
+! debug         alv=lv0-clmcpv*(ticb(i)-t0)
     alv = lv0 - clmcpv*(ticb(i)-273.15)
 
-    ! First iteration.
-
-    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
-    s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
-      +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
+! First iteration.
+
+! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
+    s = cpd*(1.-qnk(i)) + cl*qnk(i) + &                   ! convect3
+        alv*alv*qg/(rrv*ticb(i)*ticb(i))                  ! convect3
     s = 1./s
-    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
+! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
     tg = tg + s*(ah0(i)-ahg)
-    ! ori          tg=max(tg,35.0)
-    ! debug          tc=tg-t0
+! ori          tg=max(tg,35.0)
+! debug          tc=tg-t0
     tc = tg - 273.15
     denom = 243.5 + tc
     denom = max(denom, 1.0) ! convect3
-    ! ori          if(tc.ge.0.0)then
+! ori          if(tc.ge.0.0)then
     es = 6.112*exp(17.67*tc/denom)
-    ! ori          else
-    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-    ! ori          endif
-    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
+! ori          else
+! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+! ori          endif
+! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
 
-    ! Second iteration.
-
-
-    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
-    ! ori          s=1./s
-    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
+! Second iteration.
+
+
+! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
+! ori          s=1./s
+! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
     ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
     tg = tg + s*(ah0(i)-ahg)
-    ! ori          tg=max(tg,35.0)
-    ! debug          tc=tg-t0
+! ori          tg=max(tg,35.0)
+! debug          tc=tg-t0
     tc = tg - 273.15
     denom = 243.5 + tc
-    denom = max(denom, 1.0) ! convect3
-    ! ori          if(tc.ge.0.0)then
+    denom = max(denom, 1.0)                               ! convect3
+! ori          if(tc.ge.0.0)then
     es = 6.112*exp(17.67*tc/denom)
-    ! ori          else
-    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-    ! ori          end if
-    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
+! ori          else
+! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+! ori          end if
+! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     qg = eps*es/(p(i,icbs(i))-es*(1.-eps))
 
     alv = lv0 - clmcpv*(ticb(i)-273.15)
 
-    ! ori c approximation here:
-    ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
-    ! ori     &   -gz(i,icb(i))-alv*qg)/cpd
-
-    ! convect3: no approximation:
+! ori c approximation here:
+! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
+! ori     &   -gz(i,icb(i))-alv*qg)/cpd
+
+! convect3: no approximation:
     tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i))
 
-    ! ori         clw(i,icb(i))=qnk(i)-qg
-    ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
+! ori         clw(i,icb(i))=qnk(i)-qg
+! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
     clw(i, icbs(i)) = qnk(i) - qg
     clw(i, icbs(i)) = max(0.0, clw(i,icbs(i)))
 
     rg = qg/(1.-qnk(i))
-    ! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
-    ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
-    tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing
-
-  END DO
-
-  ! ori      do 380 k=minorig,icbsmax2
-  ! ori       do 370 i=1,len
-  ! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
-  ! ori 370   continue
-  ! ori 380  continue
-
-
-  ! -- The following is only for convect3:
-
-  ! * icbs is the first level above the LCL:
-  ! if plcl<p(icb), then icbs=icb+1
-  ! if plcl>p(icb), then icbs=icb
-
-  ! * the routine above computes tvp from minorig to icbs (included).
-
-  ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
-  ! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
-
-  ! * therefore, in the case icbs=icb, we compute tvp at level icb+1
-  ! (tvp at other levels will be computed in cv3_undilute2.F)
+! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
+! convect3: (qg utilise au lieu du vrai mixing ratio rg)
+    tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i))   !whole thing
+
+  END DO
+
+! ori      do 380 k=minorig,icbsmax2
+! ori       do 370 i=1,len
+! ori         tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i)
+! ori 370   continue
+! ori 380  continue
+
+
+! -- The following is only for convect3:
+
+! * icbs is the first level above the LCL:
+! if plcl<p(icb), then icbs=icb+1
+! if plcl>p(icb), then icbs=icb
+
+! * the routine above computes tvp from minorig to icbs (included).
+
+! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1)
+! must be known. This is the case if icbs=icb+1, but not if icbs=icb.
+
+! * therefore, in the case icbs=icb, we compute tvp at level icb+1
+! (tvp at other levels will be computed in cv3_undilute2.F)
 
 
@@ -615,68 +683,68 @@
     tg = ticb(i)
     qg = qsicb(i) ! convect3
-    ! debug         alv=lv0-clmcpv*(ticb(i)-t0)
+! debug         alv=lv0-clmcpv*(ticb(i)-t0)
     alv = lv0 - clmcpv*(ticb(i)-273.15)
 
-    ! First iteration.
-
-    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
-    s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
-      +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3
+! First iteration.
+
+! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
+    s = cpd*(1.-qnk(i)) + cl*qnk(i) &                         ! convect3
+      +alv*alv*qg/(rrv*ticb(i)*ticb(i))                       ! convect3
     s = 1./s
-    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
-    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
+! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
+    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i)     ! convect3
     tg = tg + s*(ah0(i)-ahg)
-    ! ori          tg=max(tg,35.0)
-    ! debug          tc=tg-t0
+! ori          tg=max(tg,35.0)
+! debug          tc=tg-t0
     tc = tg - 273.15
     denom = 243.5 + tc
-    denom = max(denom, 1.0) ! convect3
-    ! ori          if(tc.ge.0.0)then
+    denom = max(denom, 1.0)                                   ! convect3
+! ori          if(tc.ge.0.0)then
     es = 6.112*exp(17.67*tc/denom)
-    ! ori          else
-    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-    ! ori          endif
-    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
+! ori          else
+! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+! ori          endif
+! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
 
-    ! Second iteration.
-
-
-    ! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
-    ! ori          s=1./s
-    ! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
-    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3
+! Second iteration.
+
+
+! ori          s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i))
+! ori          s=1./s
+! ori          ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i)
+    ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i)     ! convect3
     tg = tg + s*(ah0(i)-ahg)
-    ! ori          tg=max(tg,35.0)
-    ! debug          tc=tg-t0
+! ori          tg=max(tg,35.0)
+! debug          tc=tg-t0
     tc = tg - 273.15
     denom = 243.5 + tc
-    denom = max(denom, 1.0) ! convect3
-    ! ori          if(tc.ge.0.0)then
+    denom = max(denom, 1.0)                                   ! convect3
+! ori          if(tc.ge.0.0)then
     es = 6.112*exp(17.67*tc/denom)
-    ! ori          else
-    ! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-    ! ori          end if
-    ! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
+! ori          else
+! ori           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+! ori          end if
+! ori          qg=eps*es/(p(i,icb(i))-es*(1.-eps))
     qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps))
 
     alv = lv0 - clmcpv*(ticb(i)-273.15)
 
-    ! ori c approximation here:
-    ! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
-    ! ori     &   -gz(i,icb(i))-alv*qg)/cpd
-
-    ! convect3: no approximation:
+! ori c approximation here:
+! ori         tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i)
+! ori     &   -gz(i,icb(i))-alv*qg)/cpd
+
+! convect3: no approximation:
     tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i))
 
-    ! ori         clw(i,icb(i))=qnk(i)-qg
-    ! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
+! ori         clw(i,icb(i))=qnk(i)-qg
+! ori         clw(i,icb(i))=max(0.0,clw(i,icb(i)))
     clw(i, icb(i)+1) = qnk(i) - qg
     clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1))
 
     rg = qg/(1.-qnk(i))
-    ! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
-    ! convect3: (qg utilise au lieu du vrai mixing ratio rg)
-    tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing
+! ori         tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi)
+! convect3: (qg utilise au lieu du vrai mixing ratio rg)
+    tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i))     !whole thing
 
   END DO
@@ -685,26 +753,26 @@
 END SUBROUTINE cv3_undilute1
 
-SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, pbase, &
-    buoybase, iflag, sig, w0)
+SUBROUTINE cv3_trigger(len, nd, icb, plcl, p, th, tv, tvp, thnk, &
+                       pbase, buoybase, iflag, sig, w0)
   IMPLICIT NONE
 
-  ! -------------------------------------------------------------------
-  ! --- TRIGGERING
-
-  ! - computes the cloud base
-  ! - triggering (crude in this version)
-  ! - relaxation of sig and w0 when no convection
-
-  ! Caution1: if no convection, we set iflag=4
-  ! (it used to be 0 in convect3)
-
-  ! Caution2: at this stage, tvp (and thus buoy) are know up
-  ! through icb only!
-  ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
-  ! -------------------------------------------------------------------
+! -------------------------------------------------------------------
+! --- TRIGGERING
+
+! - computes the cloud base
+! - triggering (crude in this version)
+! - relaxation of sig and w0 when no convection
+
+! Caution1: if no convection, we set iflag=4
+! (it used to be 0 in convect3)
+
+! Caution2: at this stage, tvp (and thus buoy) are know up
+! through icb only!
+! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy
+! -------------------------------------------------------------------
 
   include "cv3param.h"
 
-  ! input:
+! input:
   INTEGER len, nd
   INTEGER icb(len)
@@ -713,57 +781,55 @@
   REAL thnk(len)
 
-  ! output:
+! output:
   REAL pbase(len), buoybase(len)
 
-  ! input AND output:
+! input AND output:
   INTEGER iflag(len)
   REAL sig(len, nd), w0(len, nd)
 
-  ! local variables:
+! local variables:
   INTEGER i, k
   REAL tvpbase, tvbase, tdif, ath, ath1
 
 
-  ! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
+! ***   set cloud base buoyancy at (plcl+dpbase) level buoyancy
 
   DO i = 1, len
     pbase(i) = plcl(i) + dpbase
-    tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
-      (p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( &
-      p(i,icb(i))-p(i,icb(i)+1))
-    tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ &
-      (p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p &
-      (i,icb(i))-p(i,icb(i)+1))
+    tvpbase = tvp(i, icb(i))  *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
+              tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))  /(p(i,icb(i))-p(i,icb(i)+1))
+    tvbase = tv(i, icb(i))  *(pbase(i)-p(i,icb(i)+1))/(p(i,icb(i))-p(i,icb(i)+1)) + &
+             tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))  /(p(i,icb(i))-p(i,icb(i)+1))
     buoybase(i) = tvpbase - tvbase
   END DO
 
 
-  ! ***   make sure that column is dry adiabatic between the surface  ***
-  ! ***    and cloud base, and that lifted air is positively buoyant  ***
-  ! ***                         at cloud base                         ***
-  ! ***       if not, return to calling program after resetting       ***
-  ! ***                        sig(i) and w0(i)                       ***
-
-
-  ! oct3      do 200 i=1,len
-  ! oct3
-  ! oct3       tdif = buoybase(i)
-  ! oct3       ath1 = th(i,1)
-  ! oct3       ath  = th(i,icb(i)-1) - dttrig
-  ! oct3
-  ! oct3       if (tdif.lt.dtcrit .or. ath.gt.ath1) then
-  ! oct3         do 60 k=1,nl
-  ! oct3            sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
-  ! oct3            sig(i,k) = AMAX1(sig(i,k),0.0)
-  ! oct3            w0(i,k)  = beta*w0(i,k)
-  ! oct3   60    continue
-  ! oct3         iflag(i)=4 ! pour version vectorisee
-  ! oct3c convect3         iflag(i)=0
-  ! oct3cccc         return
-  ! oct3       endif
-  ! oct3
-  ! oct3200   continue
-
-  ! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
+! ***   make sure that column is dry adiabatic between the surface  ***
+! ***    and cloud base, and that lifted air is positively buoyant  ***
+! ***                         at cloud base                         ***
+! ***       if not, return to calling program after resetting       ***
+! ***                        sig(i) and w0(i)                       ***
+
+
+! oct3      do 200 i=1,len
+! oct3
+! oct3       tdif = buoybase(i)
+! oct3       ath1 = th(i,1)
+! oct3       ath  = th(i,icb(i)-1) - dttrig
+! oct3
+! oct3       if (tdif.lt.dtcrit .or. ath.gt.ath1) then
+! oct3         do 60 k=1,nl
+! oct3            sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif
+! oct3            sig(i,k) = AMAX1(sig(i,k),0.0)
+! oct3            w0(i,k)  = beta*w0(i,k)
+! oct3   60    continue
+! oct3         iflag(i)=4 ! pour version vectorisee
+! oct3c convect3         iflag(i)=0
+! oct3cccc         return
+! oct3       endif
+! oct3
+! oct3200   continue
+
+! -- oct3: on reecrit la boucle 200 (pour la vectorisation)
 
   DO k = 1, nl
@@ -779,5 +845,5 @@
         w0(i, k) = beta*w0(i, k)
         iflag(i) = 4 ! pour version vectorisee
-        ! convect3         iflag(i)=0
+! convect3         iflag(i)=0
       END IF
 
@@ -785,14 +851,22 @@
   END DO
 
-  ! fin oct3 --
+! fin oct3 --
 
   RETURN
 END SUBROUTINE cv3_trigger
 
-SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, &
-    plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, &
-    th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, &
-    iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, &
-    v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0)
+SUBROUTINE cv3_compress(len, nloc, ncum, nd, ntra, &
+                        iflag1, nk1, icb1, icbs1, &
+                        plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, &
+                        t1, q1, qs1, u1, v1, gz1, th1, &
+                        tra1, &
+                        h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
+                        sig1, w01, &
+                        iflag, nk, icb, icbs, &
+                        plcl, tnk, qnk, gznk, pbase, buoybase, &
+                        t, q, qs, u, v, gz, th, &
+                        tra, &
+                        h, lv, cpn, p, ph, tv, tp, tvp, clw, &
+                        sig, w0)
   IMPLICIT NONE
 
@@ -800,5 +874,5 @@
   include 'iniprint.h'
 
-  ! inputs:
+!inputs:
   INTEGER len, ncum, nd, ntra, nloc
   INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len)
@@ -813,6 +887,6 @@
   REAL tra1(len, nd, ntra)
 
-  ! outputs:
-  ! en fait, on a nloc=len pour l'instant (cf cv_driver)
+!outputs:
+! en fait, on a nloc=len pour l'instant (cf cv_driver)
   INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc)
   REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
@@ -826,5 +900,5 @@
   REAL tra(nloc, nd, ntra)
 
-  ! local variables:
+!local variables:
   INTEGER i, k, nn, j
 
@@ -859,16 +933,16 @@
   END DO
 
-  ! AC!      do 121 j=1,ntra
-  ! AC!ccccc      do 111 k=1,nl+1
-  ! AC!      do 111 k=1,nd
-  ! AC!       nn=0
-  ! AC!      do 101 i=1,len
-  ! AC!      if(iflag1(i).eq.0)then
-  ! AC!       nn=nn+1
-  ! AC!       tra(nn,k,j)=tra1(i,k,j)
-  ! AC!      endif
-  ! AC! 101  continue
-  ! AC! 111  continue
-  ! AC! 121  continue
+!AC!      do 121 j=1,ntra
+!AC!ccccc      do 111 k=1,nl+1
+!AC!      do 111 k=1,nd
+!AC!       nn=0
+!AC!      do 101 i=1,len
+!AC!      if(iflag1(i).eq.0)then
+!AC!       nn=nn+1
+!AC!       tra(nn,k,j)=tra1(i,k,j)
+!AC!      endif
+!AC! 101  continue
+!AC! 111  continue
+!AC! 121  continue
 
   IF (nn/=ncum) THEN
@@ -902,7 +976,7 @@
 
 
-  ! JAM--------------------------------------------------------------------
-  ! Calcul de la quantité d'eau sous forme de glace
-  ! --------------------------------------------------------------------
+!JAM--------------------------------------------------------------------
+! Calcul de la quantité d'eau sous forme de glace
+! --------------------------------------------------------------------
   REAL qi(len, nl)
   REAL t(len, nl), clw(len, nl)
@@ -922,5 +996,5 @@
         END IF
       END IF
-      ! print*,t(i,k),qi(i,k),'temp,testglace'
+! print*,t(i,k),qi(i,k),'temp,testglace'
     END DO
   END DO
@@ -930,26 +1004,27 @@
 END SUBROUTINE icefrac
 
-SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, hnk, &
-    t, q, qs, gz, p, h, tv, lv, lf, pbase, buoybase, plcl, inb, tp, tvp, clw, &
-    hp, ep, sigp, buoy, frac)
+SUBROUTINE cv3_undilute2(nloc, ncum, nd, icb, icbs, nk, &
+                         tnk, qnk, gznk, hnk, t, q, qs, gz, &
+                         p, h, tv, lv, lf, pbase, buoybase, plcl, &
+                         inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
   IMPLICIT NONE
 
-  ! ---------------------------------------------------------------------
-  ! Purpose:
-  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
-  ! &
-  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
-  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
-  ! &
-  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
-
-  ! Main differences convect3/convect4:
-  ! - icbs (input) is the first level above LCL (may differ from icb)
-  ! - many minor differences in the iterations
-  ! - condensed water not removed from tvp in convect3
-  ! - vertical profile of buoyancy computed here (use of buoybase)
-  ! - the determination of inb is different
-  ! - no inb1, only inb in output
-  ! ---------------------------------------------------------------------
+! ---------------------------------------------------------------------
+! Purpose:
+! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
+! &
+! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
+! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
+! &
+! FIND THE LEVEL OF NEUTRAL BUOYANCY
+
+! Main differences convect3/convect4:
+!   - icbs (input) is the first level above LCL (may differ from icb)
+!   - many minor differences in the iterations
+!   - condensed water not removed from tvp in convect3
+!   - vertical profile of buoyancy computed here (use of buoybase)
+!   - the determination of inb is different
+!   - no inb1, only inb in output
+! ---------------------------------------------------------------------
 
   include "cvthermo.h"
@@ -958,5 +1033,5 @@
   include "cvflag.h"
 
-  ! inputs:
+!inputs:
   INTEGER ncum, nd, nloc, j
   INTEGER icb(nloc), icbs(nloc), nk(nloc)
@@ -968,5 +1043,5 @@
   REAL pbase(nloc), buoybase(nloc), plcl(nloc)
 
-  ! outputs:
+!outputs:
   INTEGER inb(nloc)
   REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
@@ -974,5 +1049,5 @@
   REAL buoy(nloc, nd)
 
-  ! local variables:
+!local variables:
   INTEGER i, k
   REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
@@ -986,7 +1061,7 @@
   REAL fracg
 
-  ! =====================================================================
-  ! --- SOME INITIALIZATIONS
-  ! =====================================================================
+! =====================================================================
+! --- SOME INITIALIZATIONS
+! =====================================================================
 
   DO k = 1, nl
@@ -998,84 +1073,83 @@
   END DO
 
-  ! =====================================================================
-  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
-  ! =====================================================================
-
-  ! ---       The procedure is to solve the equation.
-  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
-
-  ! ***  Calculate certain parcel quantities, including static energy   ***
+! =====================================================================
+! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
+! =====================================================================
+
+! ---       The procedure is to solve the equation.
+!                cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
+
+! ***  Calculate certain parcel quantities, including static energy   ***
 
 
   DO i = 1, ncum
-    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) & ! debug     &
-                                                  ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
-      +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
-  END DO
-
-
-  ! ***  Find lifted parcel quantities above cloud base    ***
+    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i)+ & 
+! debug          qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i)
+             qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
+  END DO
+
+
+! ***  Find lifted parcel quantities above cloud base    ***
 
 
   DO k = minorig + 1, nl
     DO i = 1, ncum
-      ! ori	    if(k.ge.(icb(i)+1))then
-      IF (k>=(icbs(i)+1)) THEN ! convect3
+! ori	    if(k.ge.(icb(i)+1))then
+      IF (k>=(icbs(i)+1)) THEN                                ! convect3
         tg = t(i, k)
         qg = qs(i, k)
-        ! debug	      alv=lv0-clmcpv*(t(i,k)-t0)
+! debug	      alv=lv0-clmcpv*(t(i,k)-t0)
         alv = lv0 - clmcpv*(t(i,k)-273.15)
 
-        ! First iteration.
-
-        ! ori	       s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
-        s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3
-          +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3
+! First iteration.
+
+! ori	       s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
+        s = cpd*(1.-qnk(i)) + cl*qnk(i) + &                   ! convect3
+            alv*alv*qg/(rrv*t(i,k)*t(i,k))                    ! convect3
         s = 1./s
-        ! ori	       ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
+! ori	       ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
         ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
         tg = tg + s*(ah0(i)-ahg)
-        ! ori	       tg=max(tg,35.0)
-        ! debug	       tc=tg-t0
+! ori	       tg=max(tg,35.0)
+! debug	       tc=tg-t0
         tc = tg - 273.15
         denom = 243.5 + tc
-        denom = max(denom, 1.0) ! convect3
-        ! ori	       if(tc.ge.0.0)then
+        denom = max(denom, 1.0)                               ! convect3
+! ori	       if(tc.ge.0.0)then
         es = 6.112*exp(17.67*tc/denom)
-        ! ori	       else
-        ! ori			es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-        ! ori	       endif
+! ori	       else
+! ori			es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+! ori	       endif
         qg = eps*es/(p(i,k)-es*(1.-eps))
 
-        ! Second iteration.
-
-        ! ori	       s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
-        ! ori	       s=1./s
-        ! ori	       ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
+! Second iteration.
+
+! ori	       s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k))
+! ori	       s=1./s
+! ori	       ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k)
         ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3
         tg = tg + s*(ah0(i)-ahg)
-        ! ori	       tg=max(tg,35.0)
-        ! debug	       tc=tg-t0
+! ori	       tg=max(tg,35.0)
+! debug	       tc=tg-t0
         tc = tg - 273.15
         denom = 243.5 + tc
-        denom = max(denom, 1.0) ! convect3
-        ! ori	       if(tc.ge.0.0)then
+        denom = max(denom, 1.0)                               ! convect3
+! ori	       if(tc.ge.0.0)then
         es = 6.112*exp(17.67*tc/denom)
-        ! ori	       else
-        ! ori			es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
-        ! ori	       endif
+! ori	       else
+! ori			es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
+! ori	       endif
         qg = eps*es/(p(i,k)-es*(1.-eps))
 
-        ! debug	       alv=lv0-clmcpv*(t(i,k)-t0)
+! debug	       alv=lv0-clmcpv*(t(i,k)-t0)
         alv = lv0 - clmcpv*(t(i,k)-273.15)
-        ! print*,'cpd dans convect2 ',cpd
-        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
-        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
-
-        ! ori c approximation here:
-        ! ori
-        ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
-
-        ! convect3: no approximation:
+! print*,'cpd dans convect2 ',cpd
+! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
+! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
+
+! ori c approximation here:
+! ori        tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
+
+! convect3: no approximation:
         IF (cvflag_ice) THEN
           tp(i, k) = max(0., (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)))
@@ -1087,6 +1161,6 @@
         clw(i, k) = max(0.0, clw(i,k))
         rg = qg/(1.-qnk(i))
-        ! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
-        ! convect3: (qg utilise au lieu du vrai mixing ratio rg):
+! ori               tvp(i,k)=tp(i,k)*(1.+rg*epsi)
+! convect3: (qg utilise au lieu du vrai mixing ratio rg):
         tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing
         IF (cvflag_ice) THEN
@@ -1099,6 +1173,6 @@
 
       IF (cvflag_ice) THEN
-        ! CR:attention boucle en klon dans Icefrac
-        ! Call Icefrac(t,clw,qi,nl,nloc)
+!CR:attention boucle en klon dans Icefrac
+! Call Icefrac(t,clw,qi,nl,nloc)
         IF (t(i,k)>263.15) THEN
           qi(i, k) = 0.
@@ -1111,29 +1185,29 @@
           END IF
         END IF
-        ! CR: fin test
+!CR: fin test
         IF (t(i,k)<263.15) THEN
-          ! CR: on commente les calculs d'Arnaud car division par zero
-          ! nouveau calcul propose par JYG
-          ! alv=lv0-clmcpv*(t(i,k)-273.15)
-          ! alf=lf0-clmci*(t(i,k)-273.15)
-          ! tg=tp(i,k)
-          ! tc=tp(i,k)-273.15
-          ! denom=243.5+tc
-          ! do j=1,3
-          ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-          ! il faudra que esi vienne en argument de la convection
-          ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-          ! tbis=t(i,k)+(tp(i,k)-tg)
-          ! esi=exp(23.33086-(6111.72784/tbis)
-          ! :               +0.15215*log(tbis))
-          ! qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
-          ! snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/
-          ! :                               (rrv*tbis*tbis)
-          ! snew=1./snew
-          ! print*,esi,qsat_new,snew,'esi,qsat,snew'
-          ! tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
-          ! print*,k,tp(i,k),qnk(i),'avec glace'
-          ! print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
-          ! enddo
+!CR: on commente les calculs d'Arnaud car division par zero
+! nouveau calcul propose par JYG
+!       alv=lv0-clmcpv*(t(i,k)-273.15)
+!       alf=lf0-clmci*(t(i,k)-273.15)
+!       tg=tp(i,k)
+!       tc=tp(i,k)-273.15
+!       denom=243.5+tc
+!       do j=1,3
+! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! il faudra que esi vienne en argument de la convection
+! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!        tbis=t(i,k)+(tp(i,k)-tg)
+!        esi=exp(23.33086-(6111.72784/tbis) + &
+!                       0.15215*log(tbis))
+!        qsat_new=eps*esi/(p(i,k)-esi*(1.-eps))
+!        snew=cpd*(1.-qnk(i))+cl*qnk(i)+alv*alv*qsat_new/ &
+!                                       (rrv*tbis*tbis)
+!        snew=1./snew
+!        print*,esi,qsat_new,snew,'esi,qsat,snew'
+!        tp(i,k)=tg+(alf*qi(i,k)+alv*qg*(1.-(esi/es)))*snew
+!        print*,k,tp(i,k),qnk(i),'avec glace'
+!        print*,'tpNAN',tg,alf,qi(i,k),alv,qg,esi,es,snew
+!       enddo
 
           alv = lv0 - clmcpv*(t(i,k)-273.15)
@@ -1145,19 +1219,20 @@
             esi = exp(23.33086-(6111.72784/tp(i,k))+0.15215*log(tp(i,k)))
             qsat_new = eps*esi/(p(i,k)-esi*(1.-eps))
-            snew = cpd*(1.-qnk(i)) + cl*qnk(i) + alv*als*qsat_new/(rrv*tp(i,k &
-              )*tp(i,k))
+            snew = cpd*(1.-qnk(i)) + cl*qnk(i) + alv*als*qsat_new/ &
+                                                 (rrv*tp(i,k)*tp(i,k))
             snew = 1./snew
-            ! c             print*,esi,qsat_new,snew,'esi,qsat,snew'
-            tp(i, k) = tp(i, k) + ((cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i, &
-              k))+alv*(qg-qsat_new)+alf*qi(i,k))*snew
-            ! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k),
-            ! :             'k,tp,q,qt,qi avec glace'
+! c             print*,esi,qsat_new,snew,'esi,qsat,snew'
+            tp(i, k) = tp(i, k) + &
+                       ((cpd*(1.-qnk(i))+cl*qnk(i))*(tg-tp(i,k)) + &
+                        alv*(qg-qsat_new)+alf*qi(i,k))*snew
+! print*,k,tp(i,k),qsat_new,qnk(i),qi(i,k), &
+!              'k,tp,q,qt,qi avec glace'
           END DO
 
-          ! CR:reprise du code AJ
+!CR:reprise du code AJ
           clw(i, k) = qnk(i) - qsat_new
           clw(i, k) = max(0.0, clw(i,k))
           tvp(i, k) = max(0., tp(i,k)*(1.+qsat_new/eps-qnk(i)))
-          ! print*,tvp(i,k),'tvp'
+! print*,tvp(i,k),'tvp'
         END IF
         IF (clw(i,k)<1.E-11) THEN
@@ -1170,9 +1245,9 @@
   END DO
 
-  ! =====================================================================
-  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
-  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
-  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
-  ! =====================================================================
+! =====================================================================
+! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
+! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
+! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
+! =====================================================================
 
   IF (flag_epkeorig/=1) THEN
@@ -1205,37 +1280,37 @@
     END DO
   END IF
-  ! =====================================================================
-  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
-  ! --- VIRTUAL TEMPERATURE
-  ! =====================================================================
-
-  ! dans convect3, tvp est calcule en une seule fois, et sans retirer
-  ! l'eau condensee (~> reversible CAPE)
-
-  ! ori      do 340 k=minorig+1,nl
-  ! ori        do 330 i=1,ncum
-  ! ori        if(k.ge.(icb(i)+1))then
-  ! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
-  ! oric         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
-  ! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
-  ! ori        endif
-  ! ori 330    continue
-  ! ori 340  continue
-
-  ! ori      do 350 i=1,ncum
-  ! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
-  ! ori 350  continue
-
-  DO i = 1, ncum ! convect3
-    tp(i, nlp) = tp(i, nl) ! convect3
-  END DO ! convect3
-
-  ! =====================================================================
-  ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
-  ! =====================================================================
-
-  ! -- this is for convect3 only:
-
-  ! first estimate of buoyancy:
+! =====================================================================
+! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
+! --- VIRTUAL TEMPERATURE
+! =====================================================================
+
+! dans convect3, tvp est calcule en une seule fois, et sans retirer
+! l'eau condensee (~> reversible CAPE)
+
+! ori      do 340 k=minorig+1,nl
+! ori        do 330 i=1,ncum
+! ori        if(k.ge.(icb(i)+1))then
+! ori          tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
+! oric         print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
+! oric         print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
+! ori        endif
+! ori 330    continue
+! ori 340  continue
+
+! ori      do 350 i=1,ncum
+! ori       tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd
+! ori 350  continue
+
+  DO i = 1, ncum                                           ! convect3
+    tp(i, nlp) = tp(i, nl)                                 ! convect3
+  END DO                                                   ! convect3
+
+! =====================================================================
+! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
+! =====================================================================
+
+! -- this is for convect3 only:
+
+! first estimate of buoyancy:
 
   DO i = 1, ncum
@@ -1245,6 +1320,6 @@
   END DO
 
-  ! set buoyancy=buoybase for all levels below base
-  ! for safety, set buoy(icb)=buoybase
+! set buoyancy=buoybase for all levels below base
+! for safety, set buoy(icb)=buoybase
 
   DO i = 1, ncum
@@ -1254,16 +1329,16 @@
       END IF
     END DO
-    ! buoy(icb(i),k)=buoybase(i)
+!    buoy(icb(i),k)=buoybase(i)
     buoy(i, icb(i)) = buoybase(i)
   END DO
 
-  ! -- end convect3
-
-  ! =====================================================================
-  ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
-  ! --- LEVEL OF NEUTRAL BUOYANCY
-  ! =====================================================================
-
-  ! -- this is for convect3 only:
+! -- end convect3
+
+! =====================================================================
+! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S
+! --- LEVEL OF NEUTRAL BUOYANCY
+! =====================================================================
+
+! -- this is for convect3 only:
 
   DO i = 1, ncum
@@ -1273,5 +1348,5 @@
 
 
-  ! --    iposit(i) = first level, above icb, with positive buoyancy
+! --    iposit(i) = first level, above icb, with positive buoyancy
   DO k = 1, nl - 1
     DO i = 1, ncum
@@ -1296,99 +1371,99 @@
   END DO
 
-  ! -- end convect3
-
-  ! ori      do 510 i=1,ncum
-  ! ori        cape(i)=0.0
-  ! ori        capem(i)=0.0
-  ! ori        inb(i)=icb(i)+1
-  ! ori        inb1(i)=inb(i)
-  ! ori 510  continue
-
-  ! Originial Code
-
-  ! do 530 k=minorig+1,nl-1
-  ! do 520 i=1,ncum
-  ! if(k.ge.(icb(i)+1))then
-  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
-  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
-  ! cape(i)=cape(i)+by
-  ! if(by.ge.0.0)inb1(i)=k+1
-  ! if(cape(i).gt.0.0)then
-  ! inb(i)=k+1
-  ! capem(i)=cape(i)
-  ! endif
-  ! endif
-  ! 520    continue
-  ! 530  continue
-  ! do 540 i=1,ncum
-  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
-  ! cape(i)=capem(i)+byp
-  ! defrac=capem(i)-cape(i)
-  ! defrac=max(defrac,0.001)
-  ! frac(i)=-cape(i)/defrac
-  ! frac(i)=min(frac(i),1.0)
-  ! frac(i)=max(frac(i),0.0)
-  ! 540   continue
-
-  ! K Emanuel fix
-
-  ! call zilch(byp,ncum)
-  ! do 530 k=minorig+1,nl-1
-  ! do 520 i=1,ncum
-  ! if(k.ge.(icb(i)+1))then
-  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
-  ! cape(i)=cape(i)+by
-  ! if(by.ge.0.0)inb1(i)=k+1
-  ! if(cape(i).gt.0.0)then
-  ! inb(i)=k+1
-  ! capem(i)=cape(i)
-  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
-  ! endif
-  ! endif
-  ! 520    continue
-  ! 530  continue
-  ! do 540 i=1,ncum
-  ! inb(i)=max(inb(i),inb1(i))
-  ! cape(i)=capem(i)+byp(i)
-  ! defrac=capem(i)-cape(i)
-  ! defrac=max(defrac,0.001)
-  ! frac(i)=-cape(i)/defrac
-  ! frac(i)=min(frac(i),1.0)
-  ! frac(i)=max(frac(i),0.0)
-  ! 540   continue
-
-  ! J Teixeira fix
-
-  ! ori      call zilch(byp,ncum)
-  ! ori      do 515 i=1,ncum
-  ! ori        lcape(i)=.true.
-  ! ori 515  continue
-  ! ori      do 530 k=minorig+1,nl-1
-  ! ori        do 520 i=1,ncum
-  ! ori          if(cape(i).lt.0.0)lcape(i)=.false.
-  ! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
-  ! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
-  ! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
-  ! ori            cape(i)=cape(i)+by
-  ! ori            if(by.ge.0.0)inb1(i)=k+1
-  ! ori            if(cape(i).gt.0.0)then
-  ! ori              inb(i)=k+1
-  ! ori              capem(i)=cape(i)
-  ! ori            endif
-  ! ori          endif
-  ! ori 520    continue
-  ! ori 530  continue
-  ! ori      do 540 i=1,ncum
-  ! ori          cape(i)=capem(i)+byp(i)
-  ! ori          defrac=capem(i)-cape(i)
-  ! ori          defrac=max(defrac,0.001)
-  ! ori          frac(i)=-cape(i)/defrac
-  ! ori          frac(i)=min(frac(i),1.0)
-  ! ori          frac(i)=max(frac(i),0.0)
-  ! ori 540  continue
-
-  ! =====================================================================
-  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
-  ! =====================================================================
+! -- end convect3
+
+! ori      do 510 i=1,ncum
+! ori        cape(i)=0.0
+! ori        capem(i)=0.0
+! ori        inb(i)=icb(i)+1
+! ori        inb1(i)=inb(i)
+! ori 510  continue
+
+! Originial Code
+
+!    do 530 k=minorig+1,nl-1
+!     do 520 i=1,ncum
+!      if(k.ge.(icb(i)+1))then
+!       by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
+!       byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
+!       cape(i)=cape(i)+by
+!       if(by.ge.0.0)inb1(i)=k+1
+!       if(cape(i).gt.0.0)then
+!        inb(i)=k+1
+!        capem(i)=cape(i)
+!       endif
+!      endif
+!520    continue
+!530  continue
+!    do 540 i=1,ncum
+!     byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
+!     cape(i)=capem(i)+byp
+!     defrac=capem(i)-cape(i)
+!     defrac=max(defrac,0.001)
+!     frac(i)=-cape(i)/defrac
+!     frac(i)=min(frac(i),1.0)
+!     frac(i)=max(frac(i),0.0)
+!540   continue
+
+!    K Emanuel fix
+
+!    call zilch(byp,ncum)
+!    do 530 k=minorig+1,nl-1
+!     do 520 i=1,ncum
+!      if(k.ge.(icb(i)+1))then
+!       by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
+!       cape(i)=cape(i)+by
+!       if(by.ge.0.0)inb1(i)=k+1
+!       if(cape(i).gt.0.0)then
+!        inb(i)=k+1
+!        capem(i)=cape(i)
+!        byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
+!       endif
+!      endif
+!520    continue
+!530  continue
+!    do 540 i=1,ncum
+!     inb(i)=max(inb(i),inb1(i))
+!     cape(i)=capem(i)+byp(i)
+!     defrac=capem(i)-cape(i)
+!     defrac=max(defrac,0.001)
+!     frac(i)=-cape(i)/defrac
+!     frac(i)=min(frac(i),1.0)
+!     frac(i)=max(frac(i),0.0)
+!540   continue
+
+! J Teixeira fix
+
+! ori      call zilch(byp,ncum)
+! ori      do 515 i=1,ncum
+! ori        lcape(i)=.true.
+! ori 515  continue
+! ori      do 530 k=minorig+1,nl-1
+! ori        do 520 i=1,ncum
+! ori          if(cape(i).lt.0.0)lcape(i)=.false.
+! ori          if((k.ge.(icb(i)+1)).and.lcape(i))then
+! ori            by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
+! ori            byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
+! ori            cape(i)=cape(i)+by
+! ori            if(by.ge.0.0)inb1(i)=k+1
+! ori            if(cape(i).gt.0.0)then
+! ori              inb(i)=k+1
+! ori              capem(i)=cape(i)
+! ori            endif
+! ori          endif
+! ori 520    continue
+! ori 530  continue
+! ori      do 540 i=1,ncum
+! ori          cape(i)=capem(i)+byp(i)
+! ori          defrac=capem(i)-cape(i)
+! ori          defrac=max(defrac,0.001)
+! ori          frac(i)=-cape(i)/defrac
+! ori          frac(i)=min(frac(i),1.0)
+! ori          frac(i)=max(frac(i),0.0)
+! ori 540  continue
+
+! =====================================================================
+! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
+! =====================================================================
 
   DO k = 1, nd
@@ -1405,6 +1480,6 @@
           frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
           frac(i, k) = min(max(frac(i,k),0.0), 1.0)
-          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))*ep &
-            (i, k)*clw(i, k)
+          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))* &
+                              ep(i, k)*clw(i, k)
 
         ELSE
@@ -1419,18 +1494,19 @@
 END SUBROUTINE cv3_undilute2
 
-SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, sig, &
-    w0, cape, m, iflag)
+SUBROUTINE cv3_closure(nloc, ncum, nd, icb, inb, &
+                       pbase, p, ph, tv, buoy, &
+                       sig, w0, cape, m, iflag)
   IMPLICIT NONE
 
-  ! ===================================================================
-  ! ---  CLOSURE OF CONVECT3
-
-  ! vectorization: S. Bony
-  ! ===================================================================
+! ===================================================================
+! ---  CLOSURE OF CONVECT3
+!
+! vectorization: S. Bony
+! ===================================================================
 
   include "cvthermo.h"
   include "cv3param.h"
 
-  ! input:
+!input:
   INTEGER ncum, nd, nloc
   INTEGER icb(nloc), inb(nloc)
@@ -1439,13 +1515,13 @@
   REAL tv(nloc, nd), buoy(nloc, nd)
 
-  ! input/output:
+!input/output:
   REAL sig(nloc, nd), w0(nloc, nd)
   INTEGER iflag(nloc)
 
-  ! output:
+!output:
   REAL cape(nloc)
   REAL m(nloc, nd)
 
-  ! local variables:
+!local variables:
   INTEGER i, j, k, icbmax
   REAL deltap, fac, w, amu
@@ -1454,7 +1530,7 @@
 
 
-  ! -------------------------------------------------------
-  ! -- Initialization
-  ! -------------------------------------------------------
+! -------------------------------------------------------
+! -- Initialization
+! -------------------------------------------------------
 
   DO k = 1, nl
@@ -1464,15 +1540,15 @@
   END DO
 
-  ! -------------------------------------------------------
-  ! -- Reset sig(i) and w0(i) for i>inb and i<icb
-  ! -------------------------------------------------------
-
-  ! update sig and w0 above LNB:
+! -------------------------------------------------------
+! -- Reset sig(i) and w0(i) for i>inb and i<icb
+! -------------------------------------------------------
+
+! update sig and w0 above LNB:
 
   DO k = 1, nl - 1
     DO i = 1, ncum
       IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN
-        sig(i, k) = beta*sig(i, k) + 2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb( &
-          i)))
+        sig(i, k) = beta*sig(i, k) + &
+                    2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb(i)))
         sig(i, k) = amax1(sig(i,k), 0.0)
         w0(i, k) = beta*w0(i, k)
@@ -1481,5 +1557,5 @@
   END DO
 
-  ! compute icbmax:
+! compute icbmax:
 
   icbmax = 2
@@ -1488,10 +1564,11 @@
   END DO
 
-  ! update sig and w0 below cloud base:
+! update sig and w0 below cloud base:
 
   DO k = 1, icbmax
     DO i = 1, ncum
       IF (k<=icb(i)) THEN
-        sig(i, k) = beta*sig(i, k) - 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
+        sig(i, k) = beta*sig(i, k) - &
+                    2.*alpha*buoy(i, icb(i))*buoy(i, icb(i))
         sig(i, k) = max(sig(i,k), 0.0)
         w0(i, k) = beta*w0(i, k)
@@ -1500,23 +1577,23 @@
   END DO
 
-  ! !      if(inb.lt.(nl-1))then
-  ! !         do 85 i=inb+1,nl-1
-  ! !            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
-  ! !     1              abs(buoy(inb))
-  ! !            sig(i)=max(sig(i),0.0)
-  ! !            w0(i)=beta*w0(i)
-  ! !   85    continue
-  ! !      end if
-
-  ! !      do 87 i=1,icb
-  ! !         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
-  ! !         sig(i)=max(sig(i),0.0)
-  ! !         w0(i)=beta*w0(i)
-  ! !   87 continue
-
-  ! -------------------------------------------------------------
-  ! -- Reset fractional areas of updrafts and w0 at initial time
-  ! -- and after 10 time steps of no convection
-  ! -------------------------------------------------------------
+!!      if(inb.lt.(nl-1))then
+!!         do 85 i=inb+1,nl-1
+!!            sig(i)=beta*sig(i)+2.*alpha*buoy(inb)*
+!!     1              abs(buoy(inb))
+!!            sig(i)=max(sig(i),0.0)
+!!            w0(i)=beta*w0(i)
+!!   85    continue
+!!      end if
+
+!!      do 87 i=1,icb
+!!         sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb)
+!!         sig(i)=max(sig(i),0.0)
+!!         w0(i)=beta*w0(i)
+!!   87 continue
+
+! -------------------------------------------------------------
+! -- Reset fractional areas of updrafts and w0 at initial time
+! -- and after 10 time steps of no convection
+! -------------------------------------------------------------
 
   DO k = 1, nl - 1
@@ -1529,9 +1606,9 @@
   END DO
 
-  ! -------------------------------------------------------------
-  ! -- Calculate convective available potential energy (cape),
-  ! -- vertical velocity (w), fractional area covered by
-  ! -- undilute updraft (sig), and updraft mass flux (m)
-  ! -------------------------------------------------------------
+! -------------------------------------------------------------
+! -- Calculate convective available potential energy (cape),
+! -- vertical velocity (w), fractional area covered by
+! -- undilute updraft (sig), and updraft mass flux (m)
+! -------------------------------------------------------------
 
   DO i = 1, ncum
@@ -1539,5 +1616,5 @@
   END DO
 
-  ! compute dtmin (minimum buoyancy between ICB and given level k):
+! compute dtmin (minimum buoyancy between ICB and given level k):
 
   DO i = 1, ncum
@@ -1550,6 +1627,5 @@
     DO k = 1, nl
       DO j = minorig, nl
-        IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- &
-            1))) THEN
+        IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k-1))) THEN
           dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j))
         END IF
@@ -1558,5 +1634,5 @@
   END DO
 
-  ! the interval on which cape is computed starts at pbase :
+! the interval on which cape is computed starts at pbase :
 
   DO k = 1, nl
@@ -1570,8 +1646,8 @@
         sigold(i, k) = sig(i, k)
 
-        ! dtmin(i,k)=100.0
-        ! do 97 j=icb(i),k-1 ! mauvaise vectorisation
-        ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
-        ! 97     continue
+! dtmin(i,k)=100.0
+! do 97 j=icb(i),k-1 ! mauvaise vectorisation
+! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j))
+! 97     continue
 
         sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k))
@@ -1590,85 +1666,83 @@
   DO i = 1, ncum
     w0(i, icb(i)) = 0.5*w0(i, icb(i)+1)
-    m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/ &
-      (ph(i,icb(i)+1)-ph(i,icb(i)+2))
+    m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/(ph(i,icb(i)+1)-ph(i,icb(i)+2))
     sig(i, icb(i)) = sig(i, icb(i)+1)
     sig(i, icb(i)-1) = sig(i, icb(i))
   END DO
 
-  ! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
-  ! ccc    cloud base mass flux is exceedingly small and is decreasing (i.e.
-  ! if
-  ! ccc    the final mass flux (cbmflast) is greater than the target mass
-  ! flux
-  ! ccc    (cbmf) ??).
-  ! cc
-  ! c      do i = 1,ncum
-  ! c       cbmflast(i) = 0.
-  ! c      enddo
-  ! cc
-  ! c      do k= 1,nl
-  ! c       do i = 1,ncum
-  ! c        IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
-  ! c         cbmflast(i) = cbmflast(i)+M(i,k)
-  ! c        ENDIF
-  ! c       enddo
-  ! c      enddo
-  ! cc
-  ! c      do i = 1,ncum
-  ! c       IF (cbmflast(i) .lt. 1.e-6) THEN
-  ! c         iflag(i) = 3
-  ! c       ENDIF
-  ! c      enddo
-  ! cc
-  ! c      do k= 1,nl
-  ! c       do i = 1,ncum
-  ! c        IF (iflag(i) .ge. 3) THEN
-  ! c         M(i,k) = 0.
-  ! c         sig(i,k) = 0.
-  ! c         w0(i,k) = 0.
-  ! c        ENDIF
-  ! c       enddo
-  ! c      enddo
-  ! cc
-  ! !      cape=0.0
-  ! !      do 98 i=icb+1,inb
-  ! !         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
-  ! !         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
-  ! !         dcape=rrd*buoy(i-1)*deltap/p(i-1)
-  ! !         dlnp=deltap/p(i-1)
-  ! !         cape=max(0.0,cape)
-  ! !         sigold=sig(i)
-
-  ! !         dtmin=100.0
-  ! !         do 97 j=icb,i-1
-  ! !            dtmin=amin1(dtmin,buoy(j))
-  ! !   97    continue
-
-  ! !         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
-  ! !         sig(i)=max(sig(i),0.0)
-  ! !         sig(i)=amin1(sig(i),0.01)
-  ! !         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
-  ! !         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
-  ! !         amu=0.5*(sig(i)+sigold)*w
-  ! !         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
-  ! !         w0(i)=w
-  ! !   98 continue
-  ! !      w0(icb)=0.5*w0(icb+1)
-  ! !      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
-  ! !      sig(icb)=sig(icb+1)
-  ! !      sig(icb-1)=sig(icb)
+! ccc 3. Compute final cloud base mass flux and set iflag to 3 if
+! ccc    cloud base mass flux is exceedingly small and is decreasing (i.e. if
+! ccc    the final mass flux (cbmflast) is greater than the target mass flux
+! ccc    (cbmf) ??).
+! cc
+! c      do i = 1,ncum
+! c       cbmflast(i) = 0.
+! c      enddo
+! cc
+! c      do k= 1,nl
+! c       do i = 1,ncum
+! c        IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
+! c         cbmflast(i) = cbmflast(i)+M(i,k)
+! c        ENDIF
+! c       enddo
+! c      enddo
+! cc
+! c      do i = 1,ncum
+! c       IF (cbmflast(i) .lt. 1.e-6) THEN
+! c         iflag(i) = 3
+! c       ENDIF
+! c      enddo
+! cc
+! c      do k= 1,nl
+! c       do i = 1,ncum
+! c        IF (iflag(i) .ge. 3) THEN
+! c         M(i,k) = 0.
+! c         sig(i,k) = 0.
+! c         w0(i,k) = 0.
+! c        ENDIF
+! c       enddo
+! c      enddo
+! cc
+!!      cape=0.0
+!!      do 98 i=icb+1,inb
+!!         deltap = min(pbase,ph(i-1))-min(pbase,ph(i))
+!!         cape=cape+rrd*buoy(i-1)*deltap/p(i-1)
+!!         dcape=rrd*buoy(i-1)*deltap/p(i-1)
+!!         dlnp=deltap/p(i-1)
+!!         cape=max(0.0,cape)
+!!         sigold=sig(i)
+
+!!         dtmin=100.0
+!!         do 97 j=icb,i-1
+!!            dtmin=amin1(dtmin,buoy(j))
+!!   97    continue
+
+!!         sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin)
+!!         sig(i)=max(sig(i),0.0)
+!!         sig(i)=amin1(sig(i),0.01)
+!!         fac=amin1(((dtcrit-dtmin)/dtcrit),1.0)
+!!         w=(1.-beta)*fac*sqrt(cape)+beta*w0(i)
+!!         amu=0.5*(sig(i)+sigold)*w
+!!         m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i)
+!!         w0(i)=w
+!!   98 continue
+!!      w0(icb)=0.5*w0(icb+1)
+!!      m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2))
+!!      sig(icb)=sig(icb+1)
+!!      sig(icb-1)=sig(icb)
 
   RETURN
 END SUBROUTINE cv3_closure
 
-SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, &
-    u, v, tra, h, lv, lf, frac, qnk, unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
-    ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
+SUBROUTINE cv3_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
+                      ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
+                      unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
+                      ment, qent, uent, vent, nent, sij, elij, ments, qents, traent)
   IMPLICIT NONE
 
-  ! ---------------------------------------------------------------------
-  ! a faire:
-  ! - vectorisation de la partie normalisation des flux (do 789...)
-  ! ---------------------------------------------------------------------
+! ---------------------------------------------------------------------
+! a faire:
+! - vectorisation de la partie normalisation des flux (do 789...)
+! ---------------------------------------------------------------------
 
   include "cvthermo.h"
@@ -1676,5 +1750,5 @@
   include "cvflag.h"
 
-  ! inputs:
+!inputs:
   INTEGER ncum, nd, na, ntra, nloc
   INTEGER icb(nloc), inb(nloc), nk(nloc)
@@ -1690,5 +1764,5 @@
   REAL m(nloc, na) ! input of convect3
 
-  ! outputs:
+!outputs:
   REAL ment(nloc, na, na), qent(nloc, na, na)
   REAL uent(nloc, na, na), vent(nloc, na, na)
@@ -1699,5 +1773,5 @@
   INTEGER nent(nloc, nd)
 
-  ! local variables:
+!local variables:
   INTEGER i, j, k, il, im, jm
   INTEGER num1, num2
@@ -1710,19 +1784,19 @@
   LOGICAL lwork(nloc)
 
-  ! =====================================================================
-  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
-  ! =====================================================================
-
-  ! ori        do 360 i=1,ncum*nlp
+! =====================================================================
+! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
+! =====================================================================
+
+! ori        do 360 i=1,ncum*nlp
   DO j = 1, nl
     DO i = 1, ncum
       nent(i, j) = 0
-      ! in convect3, m is computed in cv3_closure
-      ! ori          m(i,1)=0.0
-    END DO
-  END DO
-
-  ! ori      do 400 k=1,nlp
-  ! ori       do 390 j=1,nlp
+! in convect3, m is computed in cv3_closure
+! ori          m(i,1)=0.0
+    END DO
+  END DO
+
+! ori      do 400 k=1,nlp
+! ori       do 390 j=1,nlp
   DO j = 1, nl
     DO k = 1, nl
@@ -1732,30 +1806,30 @@
         vent(i, k, j) = v(i, j)
         elij(i, k, j) = 0.0
-        ! ym            ment(i,k,j)=0.0
-        ! ym            sij(i,k,j)=0.0
+!ym            ment(i,k,j)=0.0
+!ym            sij(i,k,j)=0.0
       END DO
     END DO
   END DO
 
-  ! ym
+!ym
   ment(1:ncum, 1:nd, 1:nd) = 0.0
   sij(1:ncum, 1:nd, 1:nd) = 0.0
 
-  ! AC!      do k=1,ntra
-  ! AC!       do j=1,nd  ! instead nlp
-  ! AC!        do i=1,nd ! instead nlp
-  ! AC!         do il=1,ncum
-  ! AC!            traent(il,i,j,k)=tra(il,j,k)
-  ! AC!         enddo
-  ! AC!        enddo
-  ! AC!       enddo
-  ! AC!      enddo
+!AC!      do k=1,ntra
+!AC!       do j=1,nd  ! instead nlp
+!AC!        do i=1,nd ! instead nlp
+!AC!         do il=1,ncum
+!AC!            traent(il,i,j,k)=tra(il,j,k)
+!AC!         enddo
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
   zm(:, :) = 0.
 
-  ! =====================================================================
-  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
-  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
-  ! --- FRACTION (sij)
-  ! =====================================================================
+! =====================================================================
+! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
+! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
+! --- FRACTION (sij)
+! =====================================================================
 
   DO i = minorig + 1, nl
@@ -1763,6 +1837,5 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- &
-            1)) .AND. (j<=inb(il))) THEN
+        IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) .AND. (j<=inb(il))) THEN
 
           rti = qnk(il) - ep(il, i)*clw(il, i)
@@ -1771,8 +1844,8 @@
 
           IF (cvflag_ice) THEN
-            ! print*,cvflag_ice,'cvflag_ice dans do 700'
+! print*,cvflag_ice,'cvflag_ice dans do 700'
             IF (t(il,j)<=263.15) THEN
-              bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)*lf(il,j))* &
-                rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
+              bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
+                   lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
             END IF
           END IF
@@ -1791,6 +1864,5 @@
 
             IF (cvflag_ice) THEN
-              anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat &
-                *bf2)
+              anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
               denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
             ELSE
@@ -1801,18 +1873,15 @@
             IF (abs(denom)<0.01) denom = 0.01
             sij(il, i, j) = anum/denom
-            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - &
-              rs(il, j)
+            altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j)
             altem = altem - (bf2-1.)*cwat
           END IF
           IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN
             qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti
-            uent(il, i, j) = sij(il, i, j)*u(il, i) + &
-              (1.-sij(il,i,j))*unk(il)
-            vent(il, i, j) = sij(il, i, j)*v(il, i) + &
-              (1.-sij(il,i,j))*vnk(il)
-            ! !!!      do k=1,ntra
-            ! !!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
-            ! !!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
-            ! !!!      end do
+            uent(il, i, j) = sij(il, i, j)*u(il, i) + (1.-sij(il,i,j))*unk(il)
+            vent(il, i, j) = sij(il, i, j)*v(il, i) + (1.-sij(il,i,j))*vnk(il)
+!!!!      do k=1,ntra
+!!!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
+!!!!     :      +(1.-sij(il,i,j))*tra(il,nk(il),k)
+!!!!      end do
             elij(il, i, j) = altem
             elij(il, i, j) = max(0.0, elij(il,i,j))
@@ -1826,28 +1895,26 @@
     END DO
 
-    ! AC!       do k=1,ntra
-    ! AC!        do j=minorig,nl
-    ! AC!         do il=1,ncum
-    ! AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
-    ! AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
-    ! AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
-    ! AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
-    ! AC!          endif
-    ! AC!         enddo
-    ! AC!        enddo
-    ! AC!       enddo
-
-
-    ! ***   if no air can entrain at level i assume that updraft detrains
-    ! ***
-    ! ***   at that level and calculate detrained air flux and properties
-    ! ***
-
-
-    ! @      do 170 i=icb(il),inb(il)
+!AC!       do k=1,ntra
+!AC!        do j=minorig,nl
+!AC!         do il=1,ncum
+!AC!          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
+!AC!     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
+!AC!            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
+!AC!     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
+!AC!          endif
+!AC!         enddo
+!AC!        enddo
+!AC!       enddo
+
+
+! ***   if no air can entrain at level i assume that updraft detrains  ***
+! ***   at that level and calculate detrained air flux and properties  ***
+
+
+! @      do 170 i=icb(il),inb(il)
 
     DO il = 1, ncum
       IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
-        ! @      if(nent(il,i).eq.0)then
+! @      if(nent(il,i).eq.0)then
         ment(il, i, i) = m(il, i)
         qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
@@ -1855,5 +1922,5 @@
         vent(il, i, i) = vnk(il)
         elij(il, i, i) = clw(il, i)
-        ! MAF      sij(il,i,i)=1.0
+! MAF      sij(il,i,i)=1.0
         sij(il, i, i) = 0.0
       END IF
@@ -1861,20 +1928,18 @@
   END DO
 
-  ! AC!      do j=1,ntra
-  ! AC!       do i=minorig+1,nl
-  ! AC!        do il=1,ncum
-  ! AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0)
-  ! then
-  ! AC!          traent(il,i,i,j)=tra(il,nk(il),j)
-  ! AC!         endif
-  ! AC!        enddo
-  ! AC!       enddo
-  ! AC!      enddo
+!AC!      do j=1,ntra
+!AC!       do i=minorig+1,nl
+!AC!        do il=1,ncum
+!AC!         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
+!AC!          traent(il,i,i,j)=tra(il,nk(il),j)
+!AC!         endif
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 
   DO j = minorig, nl
     DO i = minorig, nl
       DO il = 1, ncum
-        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= &
-            inb(il))) THEN
+        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<=inb(il))) THEN
           sigij(il, i, j) = sij(il, i, j)
         END IF
@@ -1882,12 +1947,12 @@
     END DO
   END DO
-  ! @      enddo
-
-  ! @170   continue
-
-  ! =====================================================================
-  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
-  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
-  ! =====================================================================
+! @      enddo
+
+! @170   continue
+
+! =====================================================================
+! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
+! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
+! =====================================================================
 
   CALL zilch(asum, nloc*nd)
@@ -1915,14 +1980,14 @@
         IF (cvflag_ice) THEN
 
-          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))*(qp-rs &
-            (il,i)) + (cpv-cpd)*t(il, i)*(qp-rr(il,i))
-          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))*(rr( &
-            il,i)-qp) + (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
+          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
+                       (qp-rs(il,i)) + (cpv-cpd)*t(il, i)*(qp-rr(il,i))
+          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
+                       (rr(il,i)-qp) + (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
         ELSE
 
           anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + &
-            (cpv-cpd)*t(il, i)*(qp-rr(il,i))
+                       (cpv-cpd)*t(il, i)*(qp-rr(il,i))
           denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + &
-            (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
+                       (cpd-cpv)*t(il, i)*(rr(il,i)-qp)
         END IF
 
@@ -1940,12 +2005,14 @@
       num2 = 0
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-          il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+            lwork(il)) num2 = num2 + 1
       END DO
       IF (num2<=0) GO TO 175
 
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
+            lwork(il)) THEN
 
           IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN
@@ -1988,6 +2055,6 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il)) THEN
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
           ment(il, i, j) = ment(il, i, j)*asij(il)
         END IF
@@ -1997,6 +2064,6 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il)) THEN
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
           asum(il, i) = asum(il, i) + ment(il, i, j)
           ment(il, i, j) = ment(il, i, j)*sig(il, j)
@@ -2015,6 +2082,6 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il)) THEN
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
           ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i)
         END IF
@@ -2024,6 +2091,6 @@
     DO j = minorig, nl
       DO il = 1, ncum
-        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( &
-            il)-1) .AND. j<=inb(il)) THEN
+        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
+            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
           csum(il, i) = csum(il, i) + ment(il, i, j)
         END IF
@@ -2040,20 +2107,20 @@
         vent(il, i, i) = vnk(il)
         elij(il, i, i) = clw(il, i)
-        ! MAF        sij(il,i,i)=1.0
+! MAF        sij(il,i,i)=1.0
         sij(il, i, i) = 0.0
       END IF
     END DO ! il
 
-    ! AC!      do j=1,ntra
-    ! AC!       do il=1,ncum
-    ! AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
-    ! AC!     :     .and. csum(il,i).lt.m(il,i) ) then
-    ! AC!         traent(il,i,i,j)=tra(il,nk(il),j)
-    ! AC!        endif
-    ! AC!       enddo
-    ! AC!      enddo
+!AC!      do j=1,ntra
+!AC!       do il=1,ncum
+!AC!        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
+!AC!     :     .and. csum(il,i).lt.m(il,i) ) then
+!AC!         traent(il,i,i,j)=tra(il,nk(il),j)
+!AC!        endif
+!AC!       enddo
+!AC!      enddo
 789 END DO
 
-  ! MAF: renormalisation de MENT
+! MAF: renormalisation de MENT
   CALL zilch(zm, nloc*na)
   DO jm = 1, nd
@@ -2087,8 +2154,11 @@
 END SUBROUTINE cv3_mixing
 
-SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, t, rr, rs, &
-    gz, u, v, tra, p, ph, th, tv, lv, lf, cpn, ep, sigp, clw, m, ment, elij, &
-    delt, plcl, coef_clos, mp, rp, up, vp, trap, wt, water, evap, fondue, &
-    ice, faci, b, sigd, wdtraina, wdtrainm) ! RomP
+SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, &
+                     t, rr, rs, gz, u, v, tra, p, ph, &
+                     th, tv, lv, lf, cpn, ep, sigp, clw, &
+                     m, ment, elij, delt, plcl, coef_clos, &
+                     mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
+                     faci, b, sigd, &
+                     wdtrainA, wdtrainM)                                      ! RomP
   IMPLICIT NONE
 
@@ -2098,5 +2168,5 @@
   include "cvflag.h"
 
-  ! inputs:
+!inputs:
   INTEGER ncum, nd, na, ntra, nloc
   INTEGER icb(nloc), inb(nloc)
@@ -2112,8 +2182,8 @@
   REAL coef_clos(nloc)
 
-  ! input/output
+!input/output
   INTEGER iflag(nloc)
 
-  ! outputs:
+!outputs:
   REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na)
   REAL water(nloc, na), evap(nloc, na), wt(nloc, na)
@@ -2121,11 +2191,11 @@
   REAL trap(nloc, na, ntra)
   REAL b(nloc, na), sigd(nloc)
-  ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
-  ! lascendance adiabatique et des flux melanges Pa et Pm.
-  ! Distinction des wdtrain
-  ! Pa = wdtrainA     Pm = wdtrainM
-  REAL wdtraina(nloc, na), wdtrainm(nloc, na)
-
-  ! local variables
+! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
+! de l ascendance adiabatique et des flux melanges Pa et Pm.
+! Distinction des wdtrain
+! Pa = wdtrainA     Pm = wdtrainM
+  REAL wdtrainA(nloc, na), wdtrainM(nloc, na)
+
+!local variables
   INTEGER i, j, k, il, num1, ndp1
   REAL tinv, delti, coef
@@ -2143,5 +2213,5 @@
 
 
-  ! ------------------------------------------------------
+! ------------------------------------------------------
 
   delti = 1./delt
@@ -2170,31 +2240,31 @@
     END DO
   END DO
-  ! AC!        do k=1,ntra
-  ! AC!         do i=1,nd
-  ! AC!          do il=1,ncum
-  ! AC!           trap(il,i,k)=tra(il,i,k)
-  ! AC!          enddo
-  ! AC!         enddo
-  ! AC!        enddo
-  ! ! RomP >>>
+!AC!        do k=1,ntra
+!AC!         do i=1,nd
+!AC!          do il=1,ncum
+!AC!           trap(il,i,k)=tra(il,i,k)
+!AC!          enddo
+!AC!         enddo
+!AC!        enddo
+!! RomP >>>
   DO i = 1, nd
     DO il = 1, ncum
-      wdtraina(il, i) = 0.0
-      wdtrainm(il, i) = 0.0
-    END DO
-  END DO
-  ! ! RomP <<<
-
-  ! ***  check whether ep(inb)=0, if so, skip precipitating    ***
-  ! ***             downdraft calculation                      ***
+      wdtrainA(il, i) = 0.0
+      wdtrainM(il, i) = 0.0
+    END DO
+  END DO
+!! RomP <<<
+
+! ***  check whether ep(inb)=0, if so, skip precipitating    ***
+! ***             downdraft calculation                      ***
 
 
   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)) >= 0.0001
   END DO
 
-  ! ***  Set the fractionnal area sigd of precipitating downdraughts
+! ***  Set the fractionnal area sigd of precipitating downdraughts
   DO il = 1, ncum
     sigd(il) = sigdz*coef_clos(il)
@@ -2202,9 +2272,9 @@
 
 
-  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-  ! ***                    begin downdraft loop                    ***
-
-  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!
+! ***                    begin downdraft loop                    ***
+!
+! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   DO i = nl + 1, 1, -1
@@ -2219,9 +2289,9 @@
 
 
-    ! ***  integrate liquid water equation to find condensed water   ***
-    ! ***                and condensed water flux                    ***
-
-
-    ! ***              calculate detrained precipitation             ***
+! ***  integrate liquid water equation to find condensed water   ***
+! ***                and condensed water flux                    ***
+!
+!
+! ***              calculate detrained precipitation             ***
 
     DO il = 1, ncum
@@ -2229,8 +2299,8 @@
         IF (cvflag_grav) THEN
           wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
-          wdtraina(il, i) = wdtrain(il)/grav !   Pa   RomP
+          wdtrainA(il, i) = wdtrain(il)/grav                        !   Pa   RomP
         ELSE
           wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i)
-          wdtraina(il, i) = wdtrain(il)/10. !   Pa   RomP
+          wdtrainA(il, i) = wdtrain(il)/10.                         !   Pa   RomP
         END IF
       END IF
@@ -2245,8 +2315,8 @@
             IF (cvflag_grav) THEN
               wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
-              wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) !   Pm  RomP
+              wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i)  !   Pm  RomP
             ELSE
               wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
-              wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) !   Pm  RomP
+              wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i)   !   Pm  RomP
             END IF
           END IF
@@ -2256,6 +2326,6 @@
 
 
-    ! ***    find rain water and evaporation using provisional   ***
-    ! ***              estimates of rp(i)and rp(i-1)             ***
+! ***    find rain water and evaporation using provisional   ***
+! ***              estimates of rp(i)and rp(i-1)             ***
 
 
@@ -2283,6 +2353,6 @@
           END IF
 
-          rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, &
-            i))+gz(il,i+1)-gz(il,i))/lv(il, i)
+          rp(il, i) = rp(il, i+1) + &
+                      (cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
           rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
         END IF
@@ -2296,17 +2366,13 @@
           afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
           IF (cvflag_ice) THEN
-            afac1 = p(il, i)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il, &
-              1))
+            afac1 = p(il, i)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
           END IF
         ELSE
-          rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, &
-            i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
+          rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
           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) = 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))/ &
-            (1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
+          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))/(1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
           afac = 0.5*(afac1+afac2)
         END IF
@@ -2315,14 +2381,14 @@
         bfac = 1./(sigd(il)*wt(il,i))
 
-        ! jyg1
-        ! cc        sigt=1.0
-        ! cc        if(i.ge.icb)sigt=sigp(i)
-        ! prise en compte de la variation progressive de sigt dans
-        ! les couches icb et icb-1:
-        ! pour plcl<ph(i+1), pr1=0 & pr2=1
-        ! pour plcl>ph(i),   pr1=1 & pr2=0
-        ! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
-        ! sur le nuage, et pr2 est la proportion sous la base du
-        ! nuage.
+!JYG1
+! cc        sigt=1.0
+! cc        if(i.ge.icb)sigt=sigp(i)
+! prise en compte de la variation progressive de sigt dans
+! les couches icb et icb-1:
+! pour plcl<ph(i+1), pr1=0 & pr2=1
+! pour plcl>ph(i),   pr1=1 & pr2=0
+! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
+! sur le nuage, et pr2 est la proportion sous la base du
+! nuage.
         pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
         pr1 = max(0., min(1.,pr1))
@@ -2330,40 +2396,39 @@
         pr2 = max(0., min(1.,pr2))
         sigt = sigp(il, i)*pr1 + pr2
-        ! jyg2
-
-        ! jyg----
-        ! b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
-        ! c6 = water(il,i+1) + wdtrain(il)*bfac
-        ! c6 = prec(il,i+1) + wdtrain(il)*bfac
-        ! revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
-        ! evap(il,i)=sigt*afac*revap
-        ! water(il,i)=revap*revap
-        ! prec(il,i)=revap*revap
-        ! c        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
-        ! ',
-        ! c     $            i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
-        ! c---end jyg---
-
-        ! --------retour à la formulation originale d'Emanuel.
+!JYG2
+
+!JYG----
+!    b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
+!    c6 = water(il,i+1) + wdtrain(il)*bfac
+!    c6 = prec(il,i+1) + wdtrain(il)*bfac
+!    revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
+!    evap(il,i)=sigt*afac*revap
+!    water(il,i)=revap*revap
+!    prec(il,i)=revap*revap
+!!        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
+!!                 i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
+!!---end jyg---
+
+! --------retour à la formulation originale d'Emanuel.
         IF (cvflag_ice) THEN
 
-          ! b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
-          ! c6=prec(il,i+1)+bfac*wdtrain(il)
-          ! :    -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
-          ! if(c6.gt.0.0)then
-          ! revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
-
-          ! JAM  Attention: evap=sigt*E
-          ! Modification: evap devient l'évaporation en milieu de couche
-          ! car nécessaire dans cv3_yield
-          ! Du coup, il faut modifier pas mal d'équations...
-          ! et l'expression de afac qui devient afac1
-          ! revap=sqrt((prec(i+1)+prec(i))/2)
+!   b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
+!   c6=prec(il,i+1)+bfac*wdtrain(il) &
+!       -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
+!   if(c6.gt.0.0)then
+!   revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
+
+!JAM  Attention: evap=sigt*E
+!    Modification: evap devient l'évaporation en milieu de couche
+!    car nécessaire dans cv3_yield
+!    Du coup, il faut modifier pas mal d'équations...
+!    et l'expression de afac qui devient afac1
+!    revap=sqrt((prec(i+1)+prec(i))/2)
 
           b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac1
           c6 = prec(il, i+1) + 0.5*bfac*wdtrain(il)
-          ! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
-          ! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
-          ! print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
+! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
+! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
+! print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
           IF (c6>b6*b6+1.E-20) THEN
             revap = 2.*c6/(b6+sqrt(b6*b6+4.*c6))
@@ -2372,28 +2437,24 @@
           END IF
           prec(il, i) = max(0., 2.*revap*revap-prec(il,i+1))
-          ! print*,prec(il,i),'neige'
-
-          ! jyg    Dans sa formulation originale, Emanuel calcule
-          ! l'evaporation par:
-          ! c             evap(il,i)=sigt*afac*revap
-          ! ce qui n'est pas correct. Dans cv_routines, la formulation a été
-          ! modifiee.
-          ! Ici,l'evaporation evap est simplement calculee par l'equation de
-          ! conservation.
-          ! prec(il,i)=revap*revap
-          ! else
-          ! jyg----   Correction : si c6 <= 0, water(il,i)=0.
-          ! prec(il,i)=0.
-          ! endif
-
-          ! jyg---   Dans tous les cas, evaporation = [tt ce qui entre dans
-          ! la couche i]
-          ! moins [tt ce qui sort de la couche i]
-          ! print *, 'evap avec ice'
-          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il, &
-            i)))/(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
-
-          d6 = bfac*wdtrain(il) - 100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))* &
-            evap(il, i)
+! print*,prec(il,i),'neige'
+
+!JYG    Dans sa formulation originale, Emanuel calcule l'evaporation par:
+! c             evap(il,i)=sigt*afac*revap
+! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
+! Ici,l'evaporation evap est simplement calculee par l'equation de
+! conservation.
+! prec(il,i)=revap*revap
+! else
+!JYG----   Correction : si c6 <= 0, water(il,i)=0.
+! prec(il,i)=0.
+! endif
+
+!JYG---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
+! moins [tt ce qui sort de la couche i]
+! print *, 'evap avec ice'
+          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il,i))) / &
+                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
+
+          d6 = bfac*wdtrain(il) - 100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
           e6 = bfac*wdtrain(il)
           f6 = -100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
@@ -2415,22 +2476,22 @@
           END IF
 
-          ! water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
-          ! water(il,i)=max(water(il,i),0.)
-          ! ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
-          ! ice(il,i)=max(ice(il,i),0.)
-          ! fondue(il,i)=ice(il,i)*thaw
-          ! water(il,i)=water(il,i)+fondue(il,i)
-          ! ice(il,i)=ice(il,i)-fondue(il,i)
-
-          ! if((water(il,i)+ice(il,i)).lt.1.e-30)then
-          ! faci(il,i)=0.
-          ! else
-          ! faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
-          ! endif
+!           water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
+!           water(il,i)=max(water(il,i),0.)
+!           ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
+!           ice(il,i)=max(ice(il,i),0.)
+!           fondue(il,i)=ice(il,i)*thaw
+!           water(il,i)=water(il,i)+fondue(il,i)
+!           ice(il,i)=ice(il,i)-fondue(il,i)
+            
+!           if((water(il,i)+ice(il,i)).lt.1.e-30)then
+!             faci(il,i)=0.
+!           else
+!             faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
+!           endif
 
         ELSE
           b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
-          c6 = water(il, i+1) + bfac*wdtrain(il) - 50.*sigd(il)*bfac*(ph(il,i &
-            )-ph(il,i+1))*evap(il, i+1)
+          c6 = water(il, i+1) + bfac*wdtrain(il) - &
+               50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i+1)
           IF (c6>0.0) THEN
             revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
@@ -2439,16 +2500,16 @@
             water(il, i) = 0.
           END IF
-          ! print *, 'evap sans ice'
-          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(water(il, &
-            i+1)-water(il,i)))/(sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
+! print *, 'evap sans ice'
+          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i)))/ &
+                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
 
         END IF
       END IF !(i.le.inb(il) .and. lwork(il))
     END DO
-    ! ----------------------------------------------------------------
-
-    ! cc
-    ! ***  calculate precipitating downdraft mass flux under     ***
-    ! ***              hydrostatic approximation                 ***
+! ----------------------------------------------------------------
+
+! cc
+! ***  calculate precipitating downdraft mass flux under     ***
+! ***              hydrostatic approximation                 ***
 
     DO il = 1, ncum
@@ -2459,13 +2520,17 @@
         IF (cvflag_ice) THEN
           IF (cvflag_grav) THEN
-            mp(il, i) = 100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)*(p(il, &
-              i-1)-p(il,i))/delth+lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)*(p &
-              (il,i-1)-p(il,i))/delth+lfcp(il,i)*sigd(il)*wt(il,i)/100.* &
-              fondue(il,i)*(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
+            mp(il, i) = 100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)* &
+                                               (p(il,i-1)-p(il,i))/delth + &
+                                   lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
+                                               (p(il,i-1)-p(il,i))/delth + &
+                                   lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
+                                               (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
           ELSE
-            mp(il, i) = 10.*(lvcp(il,i)*sigd(il)*tevap(il)*(p(il,i-1)-p(il, &
-              i))/delth+lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)*(p(il, &
-              i-1)-p(il,i))/delth+lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il &
-              ,i)*(p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
+            mp(il, i) = 10.*(lvcp(il,i)*sigd(il)*tevap(il)* &
+                                                (p(il,i-1)-p(il,i))/delth + &
+                             lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
+                                                (p(il,i-1)-p(il,i))/delth + &
+                             lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
+                                                (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
 
           END IF
@@ -2473,8 +2538,8 @@
           IF (cvflag_grav) THEN
             mp(il, i) = 100.*ginv*lvcp(il, i)*sigd(il)*tevap(il)* &
-              (p(il,i-1)-p(il,i))/delth
+                                                (p(il,i-1)-p(il,i))/delth
           ELSE
             mp(il, i) = 10.*lvcp(il, i)*sigd(il)*tevap(il)* &
-              (p(il,i-1)-p(il,i))/delth
+                                                (p(il,i-1)-p(il,i))/delth
           END IF
 
@@ -2483,9 +2548,9 @@
       END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
     END DO
-    ! ----------------------------------------------------------------
-
-    ! ***           if hydrostatic assumption fails,             ***
-    ! ***   solve cubic difference equation for downdraft theta  ***
-    ! ***  and mass flux from two simultaneous differential eqns ***
+! ----------------------------------------------------------------
+
+! ***           if hydrostatic assumption fails,             ***
+! ***   solve cubic difference equation for downdraft theta  ***
+! ***  and mass flux from two simultaneous differential eqns ***
 
     DO il = 1, ncum
@@ -2493,22 +2558,21 @@
 
         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))
+                         (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))
 
         IF (amp2>(0.1*amfac)) THEN
           xf = 100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1))
-          tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i)*sigd &
-            (il)*th(il,i))
+          tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i) / &
+                              (lvcp(il,i)*sigd(il)*th(il,i))
           af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv
 
           IF (cvflag_ice) THEN
             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(il)*(1.+(lf(il,i)/lv(il,i))* &
-              faci(il,i))+(lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/(ph( &
-              il,i)-ph(il,i+1)))
+                 50.*(p(il,i-1)-p(il,i))*xf*(tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
+                (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il,i+1)))
           ELSE
 
             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(il)
+                                           50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
           END IF
 
@@ -2522,5 +2586,5 @@
             IF ((0.5*bf-sru)<0.0) fac = -1.0
             mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
-              fac*(abs(0.5*bf-sru))**tinv
+                                           fac*(abs(0.5*bf-sru))**tinv
           ELSE
             d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
@@ -2532,30 +2596,34 @@
           IF (cvflag_ice) THEN
             IF (cvflag_grav) THEN
-              ! jyg : il y a vraisemblablement une erreur dans la ligne 2
-              ! suivante:
-              ! il faut diviser par (mp(il,i)*sigd(il)*grav) et non par
-              ! (mp(il,i)+sigd(il)*0.1).
-              ! 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(il)*( &
-                1.+(lf(il,i)/lv(il,i))*faci(il,i))+(lf(il,i)/lv(il, &
-                i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il, &
-                i+1)))/(mp(il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t &
-                (il, i)/(lvcp(il,i)*sigd(il)*th(il,i))
+!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
+! il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
+! 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(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
+                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
+                           (ph(il,i)-ph(il,i+1))) / &
+                           (mp(il,i)+sigd(il)*0.1) - &
+                           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(il)*( &
-                1.+(lf(il,i)/lv(il,i))*faci(il,i))+(lf(il,i)/lv(il, &
-                i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il, &
-                i+1)))/(mp(il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t &
-                (il, i)/(lvcp(il,i)*sigd(il)*th(il,i))
+              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*&
+                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
+                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
+                           (ph(il,i)-ph(il,i+1))) / &
+                           (mp(il,i)+sigd(il)*0.1) - &
+                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
+                           (lvcp(il,i)*sigd(il)*th(il,i))
             END IF
           ELSE
             IF (cvflag_grav) THEN
-              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il)/(mp &
-                (il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/( &
-                lvcp(il,i)*sigd(il)*th(il,i))
+              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
+                           (mp(il,i)+sigd(il)*0.1) - &
+                           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(il)/(mp &
-                (il,i)+sigd(il)*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/( &
-                lvcp(il,i)*sigd(il)*th(il,i))
+              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
+                           (mp(il,i)+sigd(il)*0.1) - &
+                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
+                           (lvcp(il,i)*sigd(il)*th(il,i))
             END IF
           END IF
@@ -2564,5 +2632,5 @@
         END IF !(amp2.gt.(0.1*amfac))
 
-        ! ***         limit magnitude of mp(i) to meet cfl condition      ***
+! ***         limit magnitude of mp(i) to meet cfl condition      ***
 
         ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
@@ -2571,12 +2639,11 @@
         mp(il, i) = min(mp(il,i), ampmax)
 
-        ! ***      force mp to decrease linearly to zero                 ***
-        ! ***       between cloud base and the surface                   ***
-
-
-        ! c      if(p(il,i).gt.p(il,icb(il)))then
-        ! c
-        ! mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
-        ! c      endif
+! ***      force mp to decrease linearly to zero                 ***
+! ***       between cloud base and the surface                   ***
+
+
+! c      if(p(il,i).gt.p(il,icb(il)))then
+! c       mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
+! c      endif
         IF (ph(il,i)>0.9*plcl(il)) THEN
           mp(il, i) = mp(il, i)*(ph(il,1)-ph(il,i))/(ph(il,1)-0.9*plcl(il))
@@ -2585,7 +2652,7 @@
       END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
     END DO
-    ! ----------------------------------------------------------------
-
-    ! ***       find mixing ratio of precipitating downdraft     ***
+! ----------------------------------------------------------------
+
+! ***       find mixing ratio of precipitating downdraft     ***
 
     DO il = 1, ncum
@@ -2603,18 +2670,14 @@
 
           IF (cvflag_grav) THEN
-            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i &
-              +1)) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+ &
-              1)+evap(il,i))
+            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
+              100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
           ELSE
-            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i &
-              +1)) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il, &
-              i))
+            rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
+              5.*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
           END IF
           rp(il, i) = rp(il, i)/mp(il, i)
-          up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+1) &
-            )
+          up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+1))
           up(il, i) = up(il, i)/mp(il, i)
-          vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+1) &
-            )
+          vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+1))
           vp(il, i) = vp(il, i)/mp(il, i)
 
@@ -2623,9 +2686,9 @@
           IF (mp(il,i+1)>1.0E-16) THEN
             IF (cvflag_grav) THEN
-              rp(il, i) = rp(il, i+1) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph( &
-                il,i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1)
+              rp(il, i) = rp(il,i+1) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
+                                       (evap(il,i+1)+evap(il,i))/mp(il,i+1)
             ELSE
-              rp(il, i) = rp(il, i+1) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1))*( &
-                evap(il,i+1)+evap(il,i))/mp(il, i+1)
+              rp(il, i) = rp(il,i+1) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
+                                       (evap(il,i+1)+evap(il,i))/mp(il, i+1)
             END IF
             up(il, i) = up(il, i+1)
@@ -2639,32 +2702,32 @@
       END IF ! (i.lt.inb(il) .and. lwork(il))
     END DO
-    ! ----------------------------------------------------------------
-
-    ! ***       find tracer concentrations in precipitating downdraft     ***
-
-    ! AC!      do j=1,ntra
-    ! AC!       do il = 1,ncum
-    ! AC!       if (i.lt.inb(il) .and. lwork(il)) then
-    ! AC!c
-    ! AC!         if(mplus(il))then
-    ! AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
-    ! AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
-    ! AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
-    ! AC!         else ! if (mplus(il))
-    ! AC!          if(mp(il,i+1).gt.1.0e-16)then
-    ! AC!           trap(il,i,j)=trap(il,i+1,j)
-    ! AC!          endif
-    ! AC!         endif ! (mplus(il)) else if (.not.mplus(il))
-    ! AC!c
-    ! AC!        endif ! (i.lt.inb(il) .and. lwork(il))
-    ! AC!       enddo
-    ! AC!      end do
+! ----------------------------------------------------------------
+
+! ***       find tracer concentrations in precipitating downdraft     ***
+
+!AC!      do j=1,ntra
+!AC!       do il = 1,ncum
+!AC!       if (i.lt.inb(il) .and. lwork(il)) then
+!AC!c
+!AC!         if(mplus(il))then
+!AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
+!AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
+!AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
+!AC!         else ! if (mplus(il))
+!AC!          if(mp(il,i+1).gt.1.0e-16)then
+!AC!           trap(il,i,j)=trap(il,i+1,j)
+!AC!          endif
+!AC!         endif ! (mplus(il)) else if (.not.mplus(il))
+!AC!c
+!AC!        endif ! (i.lt.inb(il) .and. lwork(il))
+!AC!       enddo
+!AC!      end do
 
 400 END DO
-  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-  ! ***                    end of downdraft loop                    ***
-
-  ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+! ***                    end of downdraft loop                    ***
+
+! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
@@ -2672,10 +2735,17 @@
 END SUBROUTINE cv3_unsat
 
-SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, t_wake, &
-    rr_wake, s_wake, u, v, tra, gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
-    ep, clw, m, tp, mp, rp, up, vp, trap, wt, water, ice, evap, fondue, faci, &
-    b, sigd, ment, qent, hent, iflag_mix, uent, vent, nent, elij, traent, &
-    sig, tv, tvp, wghti, iflag, precip, vprecip, ft, fr, fu, fv, ftra, cbmf, &
-    upwd, dnwd, dnwd0, ma, mip, tls, tps, qcondc, wd, ftd, fqd)
+SUBROUTINE cv3_yield(nloc, ncum, nd, na, ntra, ok_conserv_q, &
+                     icb, inb, delt, &
+                     t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
+                     gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
+                     ep, clw, m, tp, mp, rp, up, vp, trap, &
+                     wt, water, ice, evap, fondue, faci, b, sigd, &
+                     ment, qent, hent, iflag_mix, uent, vent, &
+                     nent, elij, traent, sig, &
+                     tv, tvp, wghti, &
+                     iflag, precip, Vprecip, ft, fr, fu, fv, ftra, &
+                     cbmf, upwd, dnwd, dnwd0, ma, mip, &
+                     tls, tps, qcondc, wd, &
+                     ftd, fqd)
 
   IMPLICIT NONE
@@ -2686,77 +2756,75 @@
   include "conema3.h"
 
-  ! inputs:
-  ! print*,'cv3_yield apres include'
-  INTEGER iflag_mix
-  INTEGER ncum, nd, na, ntra, nloc
-  INTEGER icb(nloc), inb(nloc)
-  REAL delt
-  REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
-  REAL t_wake(nloc, nd), rr_wake(nloc, nd)
-  REAL s_wake(nloc)
-  REAL tra(nloc, nd, ntra), sig(nloc, nd)
-  REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
-  REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
-  REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
-  REAL lf(nloc, na)
-  REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
-  REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
-  REAL water(nloc, na), evap(nloc, na), b(nloc, na), sigd(nloc)
-  REAL fondue(nloc, na), faci(nloc, na), ice(nloc, na)
-  REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
-  REAL hent(nloc, na, na)
-  ! IM bug   real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
-  REAL vent(nloc, na, na), elij(nloc, na, na)
-  INTEGER nent(nloc, nd)
-  REAL traent(nloc, na, na, ntra)
-  REAL tv(nloc, nd), tvp(nloc, nd), wghti(nloc, nd)
-  ! print*,'cv3_yield declarations 1'
-  ! input/output:
-  INTEGER iflag(nloc)
-
-  ! outputs:
-  REAL precip(nloc)
-  REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
-  REAL ftd(nloc, nd), fqd(nloc, nd)
-  REAL ftra(nloc, nd, ntra)
-  REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
-  REAL dnwd0(nloc, nd), mip(nloc, nd)
-  REAL vprecip(nloc, nd+1)
-  REAL tls(nloc, nd), tps(nloc, nd)
-  REAL qcondc(nloc, nd) ! cld
-  REAL wd(nloc) ! gust
-  REAL cbmf(nloc)
-  ! print*,'cv3_yield declarations 2'
-  ! local variables:
-  INTEGER i, k, il, n, j, num1
-  REAL rat, delti
-  REAL ax, bx, cx, dx, ex
-  REAL cpinv, rdcp, dpinv
-  REAL awat(nloc)
-  REAL lvcp(nloc, na), lfcp(nloc, na), mke(nloc, na)
-  REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
-  ! !!      real up1(nloc), dn1(nloc)
-  REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
-  REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
-  REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
-  REAL th_wake(nloc, nd)
-  REAL alpha_qpos(nloc), alpha_qpos1(nloc)
-  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
-  REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
-
-  ! print*,'cv3_yield declarations 3'
-  ! -------------------------------------------------------------
-
-  ! initialization:
+!inputs:
+      INTEGER iflag_mix
+      INTEGER ncum, nd, na, ntra, nloc
+      LOGICAL ok_conserv_q
+      INTEGER icb(nloc), inb(nloc)
+      REAL delt
+      REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd)
+      REAL t_wake(nloc, nd), rr_wake(nloc, nd)
+      REAL s_wake(nloc)
+      REAL tra(nloc, nd, ntra), sig(nloc, nd)
+      REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na)
+      REAL th(nloc, na), p(nloc, nd), tp(nloc, na)
+      REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na)
+      REAL lf(nloc, na)
+      REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na)
+      REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra)
+      REAL water(nloc, na), evap(nloc, na), b(nloc, na), sigd(nloc)
+      REAL fondue(nloc, na), faci(nloc, na), ice(nloc, na)
+      REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na)
+      REAL hent(nloc, na, na)
+!IM bug   real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
+      REAL vent(nloc, na, na), elij(nloc, na, na)
+      INTEGER nent(nloc, nd)
+      REAL traent(nloc, na, na, ntra)
+      REAL tv(nloc, nd), tvp(nloc, nd), wghti(nloc, nd)
+!
+!input/output:
+      INTEGER iflag(nloc)
+!
+!outputs:
+      REAL precip(nloc)
+      REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd)
+      REAL ftd(nloc, nd), fqd(nloc, nd)
+      REAL ftra(nloc, nd, ntra)
+      REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd)
+      REAL dnwd0(nloc, nd), mip(nloc, nd)
+      REAL Vprecip(nloc, nd+1)
+      REAL tls(nloc, nd), tps(nloc, nd)
+      REAL qcondc(nloc, nd) ! cld
+      REAL wd(nloc) ! gust
+      REAL cbmf(nloc)
+!
+!local variables:
+      INTEGER i, k, il, n, j, num1
+      REAL rat, delti
+      REAL ax, bx, cx, dx, ex
+      REAL cpinv, rdcp, dpinv
+      REAL awat(nloc)
+      REAL lvcp(nloc, na), lfcp(nloc, na), mke(nloc, na)
+      REAL am(nloc), work(nloc), ad(nloc), amp1(nloc)
+!!      real up1(nloc), dn1(nloc)
+      REAL up1(nloc, nd, nd), dn1(nloc, nd, nd)
+      REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc)
+      REAL esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
+      REAL th_wake(nloc, nd)
+      REAL alpha_qpos(nloc), alpha_qpos1(nloc)
+      REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
+      REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
+
+      REAL sumdq !jyg
+! 
+! -------------------------------------------------------------
+
+! initialization:
 
   delti = 1.0/delt
-  ! print*,'cv3_yield initialisation delt', delt
-  ! precip,Vprecip,ft,fr,fu,fv,ftra
-  ! :                    ,cbmf,upwd,dnwd,dnwd0,ma,mip
-  ! :                    ,tls,tps,qcondc,wd
-  ! :                    ,ftd,fqd  )
+! print*,'cv3_yield initialisation delt', delt
+!
   DO il = 1, ncum
     precip(il) = 0.0
-    vprecip(il, nd+1) = 0.0
+    Vprecip(il, nd+1) = 0.0
     wd(il) = 0.0 ! gust
   END DO
@@ -2764,5 +2832,5 @@
   DO i = 1, nd
     DO il = 1, ncum
-      vprecip(il, i) = 0.0
+      Vprecip(il, i) = 0.0
       ft(il, i) = 0.0
       fr(il, i) = 0.0
@@ -2780,13 +2848,13 @@
     END DO
   END DO
-  ! print*,'cv3_yield initialisation 2'
-  ! AC!      do j=1,ntra
-  ! AC!       do i=1,nd
-  ! AC!        do il=1,ncum
-  ! AC!          ftra(il,i,j)=0.0
-  ! AC!        enddo
-  ! AC!       enddo
-  ! AC!      enddo
-  ! print*,'cv3_yield initialisation 3'
+! print*,'cv3_yield initialisation 2'
+!AC!      do j=1,ntra
+!AC!       do i=1,nd
+!AC!        do il=1,ncum
+!AC!          ftra(il,i,j)=0.0
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
+! print*,'cv3_yield initialisation 3'
   DO i = 1, nl
     DO il = 1, ncum
@@ -2798,29 +2866,21 @@
 
 
-  ! ***  calculate surface precipitation in mm/day     ***
+! ***  calculate surface precipitation in mm/day     ***
 
   DO il = 1, ncum
     IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
       IF (cvflag_ice) THEN
-        IF (cvflag_grav) THEN
-          precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1))*86400.* &
-            1000./(rowl*grav)
-        ELSE
-          precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1))*8640.
-        END IF
+        precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1)) &
+                              *86400.*1000./(rowl*grav)
       ELSE
-        IF (cvflag_grav) THEN
-          precip(il) = wt(il, 1)*sigd(il)*water(il, 1)*86400.*1000./ &
-            (rowl*grav)
-        ELSE
-          precip(il) = wt(il, 1)*sigd(il)*water(il, 1)*8640.
-        END IF
+        precip(il) = wt(il, 1)*sigd(il)*water(il, 1) &
+                              *86400.*1000./(rowl*grav)
       END IF
     END IF
   END DO
-  ! print*,'cv3_yield apres calcul precip'
-
-
-  ! ===  calculate vertical profile of  precipitation in kg/m2/s  ===
+! print*,'cv3_yield apres calcul precip'
+
+
+! ===  calculate vertical profile of  precipitation in kg/m2/s  ===
 
   DO i = 1, nl
@@ -2828,15 +2888,7 @@
       IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
         IF (cvflag_ice) THEN
-          IF (cvflag_grav) THEN
-            vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/grav
-          ELSE
-            vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/10.
-          END IF
+          Vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/grav
         ELSE
-          IF (cvflag_grav) THEN
-            vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/grav
-          ELSE
-            vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/10.
-          END IF
+          Vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/grav
         END IF
       END IF
@@ -2845,15 +2897,15 @@
 
 
-  ! ***  Calculate downdraft velocity scale    ***
-  ! ***  NE PAS UTILISER POUR L'INSTANT ***
-
-  ! !      do il=1,ncum
-  ! !        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
-  ! !     :                                  /(sigd(il)*p(il,icb(il)))
-  ! !      enddo
-
-
-  ! ***  calculate tendencies of lowest level potential temperature  ***
-  ! ***                      and mixing ratio                        ***
+! ***  Calculate downdraft velocity scale    ***
+! ***  NE PAS UTILISER POUR L'INSTANT ***
+
+!!      do il=1,ncum
+!!        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) &
+!!                                       /(sigd(il)*p(il,icb(il)))
+!!      enddo
+
+
+! ***  calculate tendencies of lowest level potential temperature  ***
+! ***                      and mixing ratio                        ***
 
   DO il = 1, ncum
@@ -2870,6 +2922,6 @@
   END DO
 
-  ! print*,'cv3_yield avant ft'
-  ! AM is the part of cbmf taken from the first level
+!    print*,'cv3_yield avant ft'
+! am is the part of cbmf taken from the first level
   DO il = 1, ncum
     am(il) = cbmf(il)*wghti(il, 1)
@@ -2878,46 +2930,33 @@
   DO il = 1, ncum
     IF (iflag(il)<=1) THEN
-      ! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
-      ! jyg  Correction pour conserver l'eau
-      ! cc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2))
-      ! !precip
+! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
+!JYG  Correction pour conserver l'eau
+! cc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2))          !precip
       IF (cvflag_ice) THEN
         ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1) - &
-          lfcp(il, 1)*sigd(il)*evap(il, 1)*faci(il, 1) - &
-          lfcp(il, 1)*sigd(il)*(fondue(il,1)*wt(il,1))/(100.*(ph(il,1)-ph(il, &
-          2))) !precip
+                     lfcp(il, 1)*sigd(il)*evap(il, 1)*faci(il, 1) - &
+                     lfcp(il, 1)*sigd(il)*(fondue(il,1)*wt(il,1)) / &
+                       (100.*(ph(il,1)-ph(il,2)))                             !precip
       ELSE
         ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1)
       END IF
 
-      IF (cvflag_grav) THEN
-        ft(il, 1) = ft(il, 1) - 0.009*grav*sigd(il)*mp(il, 2)*t_wake(il, 1)*b &
-          (il, 1)*work(il)
+      ft(il, 1) = ft(il, 1) - 0.009*grav*sigd(il)*mp(il, 2)*t_wake(il, 1)*b(il, 1)*work(il)
+
+      IF (cvflag_ice) THEN
+        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
+                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
+                                0.01*sigd(il)*wt(il, 1)*(ci-cpd)*ice(il, 2) * &
+                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
       ELSE
-        ft(il, 1) = ft(il, 1) - 0.09*sigd(il)*mp(il, 2)*t_wake(il, 1)*b(il, 1 &
-          )*work(il)
+        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
+                                     (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
       END IF
 
-      IF (cvflag_ice) THEN
-        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) &
-          *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
-          0.01*sigd(il)*wt(il, 1)*(ci-cpd)*ice(il, 2)* &
-          (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
-      ELSE
-        ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) &
-          *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
-      END IF
-
-      ftd(il, 1) = ft(il, 1) ! fin precip
-
-      IF (cvflag_grav) THEN !sature
-        IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
-        ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)+( &
-          gz(il,2)-gz(il,1))/cpn(il,1))
-      ELSE
-        IF ((0.1*work(il)*am(il))>=delti) iflag(il) = 1 !consistency vect
-        ft(il, 1) = ft(il, 1) + 0.1*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il, &
-          2)-gz(il,1))/cpn(il,1))
-      END IF
+      ftd(il, 1) = ft(il, 1)                                                  ! fin precip
+
+      IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
+      ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
+                                   (t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
     END IF ! iflag
   END DO
@@ -2927,117 +2966,86 @@
     IF (iflag_mix>0) THEN
       DO il = 1, ncum
-        ! FH WARNING a modifier :
+! FH WARNING a modifier :
         cpinv = 0.
-        ! cpinv=1.0/cpn(il,1)
+! cpinv=1.0/cpn(il,1)
         IF (j<=inb(il) .AND. iflag(il)<=1) THEN
-          IF (cvflag_grav) THEN
-            ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(hent( &
-              il,j,1)-h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j, &
-              1)))*cpinv
-          ELSE
-            ft(il, 1) = ft(il, 1) + 0.1*work(il)*ment(il, j, 1)*(hent(il,j,1) &
-              -h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j,1)))*cpinv
-          END IF ! cvflag_grav
+          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*ment(il, j, 1) * &
+                     (hent(il,j,1)-h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j,1)))*cpinv
         END IF ! j
       END DO
     END IF
   END DO
-    ! fin sature
+! fin sature
 
 
   DO il = 1, ncum
     IF (iflag(il)<=1) THEN
-      IF (cvflag_grav) THEN
-        ! jyg1  Correction pour mieux conserver l'eau (conformite avec
-        ! CONVECT4.3)
-        fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
-          sigd(il)*evap(il, 1)
-        ! cc     :          +sigd(il)*0.5*(evap(il,1)+evap(il,2))
-
-        fqd(il, 1) = fr(il, 1) !precip
-
-        fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il) !sature
-
-        fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il, &
-          1))+am(il)*(u(il,2)-u(il,1)))
-        fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
-          1))+am(il)*(v(il,2)-v(il,1)))
-      ELSE ! cvflag_grav
-        fr(il, 1) = 0.1*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
-          sigd(il)*evap(il, 1)
-        ! cc     :          +sigd(il)*0.5*(evap(il,1)+evap(il,2))
-        fqd(il, 1) = fr(il, 1) !precip
-        fr(il, 1) = fr(il, 1) + 0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
-        fu(il, 1) = fu(il, 1) + 0.1*work(il)*(mp(il,2)*(up(il,2)-u(il, &
-          1))+am(il)*(u(il,2)-u(il,1)))
-        fv(il, 1) = fv(il, 1) + 0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il, &
-          1))+am(il)*(v(il,2)-v(il,1)))
-      END IF ! cvflag_grav
+!JYG1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
+      fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr_wake(il,1))*work(il) + &
+                  sigd(il)*evap(il, 1)
+!!!                  sigd(il)*0.5*(evap(il,1)+evap(il,2))
+
+      fqd(il, 1) = fr(il, 1) !precip
+
+      fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)        !sature
+
+      fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1)) + &
+                                                  am(il)*(u(il,2)-u(il,1)))
+      fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il,1)) + &
+                                                  am(il)*(v(il,2)-v(il,1)))
     END IF ! iflag
   END DO ! il
 
 
-  ! AC!     do j=1,ntra
-  ! AC!      do il=1,ncum
-  ! AC!       if (iflag(il) .le. 1) then
-  ! AC!       if (cvflag_grav) then
-  ! AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
-  ! AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
-  ! AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
-  ! AC!       else
-  ! AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
-  ! AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
-  ! AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
-  ! AC!       endif
-  ! AC!       endif  ! iflag
-  ! AC!      enddo
-  ! AC!     enddo
+!AC!     do j=1,ntra
+!AC!      do il=1,ncum
+!AC!       if (iflag(il) .le. 1) then
+!AC!       if (cvflag_grav) then
+!AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
+!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
+!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
+!AC!       else
+!AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
+!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
+!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
+!AC!       endif
+!AC!       endif  ! iflag
+!AC!      enddo
+!AC!     enddo
 
   DO j = 2, nl
     DO il = 1, ncum
       IF (j<=inb(il) .AND. iflag(il)<=1) THEN
-        IF (cvflag_grav) THEN
-          fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il, &
-            j,1)-rr(il,1))
-          fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il, &
-            j,1)-u(il,1))
-          fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il, &
-            j,1)-v(il,1))
-        ELSE ! cvflag_grav
-          fr(il, 1) = fr(il, 1) + 0.1*work(il)*ment(il, j, 1)*(qent(il,j,1)- &
-            rr(il,1))
-          fu(il, 1) = fu(il, 1) + 0.1*work(il)*ment(il, j, 1)*(uent(il,j,1)-u &
-            (il,1))
-          fv(il, 1) = fv(il, 1) + 0.1*work(il)*ment(il, j, 1)*(vent(il,j,1)-v &
-            (il,1)) ! fin sature
-        END IF ! cvflag_grav
+        fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
+        fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il,j,1)-u(il,1))
+        fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il,j,1)-v(il,1))
       END IF ! j
     END DO
   END DO
 
-  ! AC!      do k=1,ntra
-  ! AC!       do j=2,nl
-  ! AC!        do il=1,ncum
-  ! AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
-  ! AC!
-  ! AC!          if (cvflag_grav) then
-  ! AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
-  ! AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
-  ! AC!          else
-  ! AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
-  ! AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
-  ! AC!          endif
-  ! AC!
-  ! AC!         endif
-  ! AC!        enddo
-  ! AC!       enddo
-  ! AC!      enddo
-  ! print*,'cv3_yield apres ft'
-
-  ! ***  calculate tendencies of potential temperature and mixing ratio  ***
-  ! ***               at levels above the lowest level                   ***
-
-  ! ***  first find the net saturated updraft and downdraft mass fluxes  ***
-  ! ***                      through each level                          ***
+!AC!      do k=1,ntra
+!AC!       do j=2,nl
+!AC!        do il=1,ncum
+!AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!
+!AC!          if (cvflag_grav) then
+!AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
+!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
+!AC!          else
+!AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
+!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
+!AC!          endif
+!AC!
+!AC!         endif
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
+! print*,'cv3_yield apres ft'
+
+! ***  calculate tendencies of potential temperature and mixing ratio  ***
+! ***               at levels above the lowest level                   ***
+
+! ***  first find the net saturated updraft and downdraft mass fluxes  ***
+! ***                      through each level                          ***
 
 
@@ -3060,5 +3068,5 @@
           END IF
         ELSE
-          ! AMP1 is the part of cbmf taken from layers I and lower
+! AMP1 is the part of cbmf taken from layers I and lower
           IF (k<=i) THEN
             amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
@@ -3093,19 +3101,13 @@
         cpinv = 1.0/cpn(il, i)
 
-        ! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
-        IF (cvflag_grav) THEN
-          IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
-        ELSE
-          IF ((0.1*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
-        END IF
-
-          ! precip
-        ! cc       ft(il,i)=
-        ! -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
+! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
+        IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto
+
+! precip
+! cc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
         IF (cvflag_ice) THEN
           ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i) - &
-            sigd(il)*lfcp(il, i)*evap(il, i)*faci(il, i) - &
-            sigd(il)*lfcp(il, i)*fondue(il, i)*wt(il, i)/(100.*(p(il, &
-            i-1)-p(il,i)))
+                       sigd(il)*lfcp(il, i)*evap(il, i)*faci(il, i) - &
+                       sigd(il)*lfcp(il, i)*fondue(il, i)*wt(il, i)/(100.*(p(il,i-1)-p(il,i)))
         ELSE
           ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i)
@@ -3114,123 +3116,74 @@
         rat = cpn(il, i-1)*cpinv
 
-        IF (cvflag_grav) THEN
-          ft(il, i) = ft(il, i) - 0.009*grav*sigd(il)*(mp(il,i+1)*t_wake(il,i &
-            )*b(il,i)-mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
-          IF (cvflag_ice) THEN
-            ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
-              , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
-              0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1)* &
-              (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
-          ELSE
-            ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
-              , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
-          END IF
-
-          ftd(il, i) = ft(il, i)
-            ! fin precip
-
-            ! sature
-          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
-            i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
-            i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
-
-
-          IF (iflag_mix==0) THEN
-            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)- &
-              h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
-          END IF
-
-        ELSE ! cvflag_grav
-          ft(il, i) = ft(il, i) - 0.09*sigd(il)*(mp(il,i+1)*t_wake(il,i)*b(il &
-            ,i)-mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
-
-          IF (cvflag_ice) THEN
-            ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
-              , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
-              0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1)* &
-              (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
-          ELSE
-            ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il &
-              , i+1)*(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
-          END IF
-
-          ftd(il, i) = ft(il, i)
-            ! fin precip
-
-            ! sature
-          ft(il, i) = ft(il, i) + 0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il, &
-            i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, &
-            i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
-
-
-          IF (iflag_mix==0) THEN
-            ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i &
-              )+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
-          END IF
-        END IF ! cvflag_grav
-
-
-        IF (cvflag_grav) THEN
-          ! sb: on ne fait pas encore la correction permettant de mieux
-          ! conserver l'eau:
-          ! jyg: correction permettant de mieux conserver l'eau:
-          ! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
-          fr(il, i) = sigd(il)*evap(il, i) + 0.01*grav*(mp(il,i+1)*(rp(il, &
-            i+1)-rr_wake(il,i))-mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
-          fqd(il, i) = fr(il, i) ! precip
-
-          fu(il, i) = 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il, &
-            i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
-          fv(il, i) = 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il, &
-            i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
-        ELSE ! cvflag_grav
-          ! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
-          fr(il, i) = sigd(il)*evap(il, i) + 0.1*(mp(il,i+1)*(rp(il, &
-            i+1)-rr_wake(il,i))-mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
-          fqd(il, i) = fr(il, i) ! precip
-
-          fu(il, i) = 0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i))-mp(il,i)*(up(il, &
-            i)-u(il,i-1)))*dpinv
-          fv(il, i) = 0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i))-mp(il,i)*(vp(il, &
-            i)-v(il,i-1)))*dpinv
-        END IF ! cvflag_grav
-
-
-        IF (cvflag_grav) THEN
-          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*(amp1(il)*(rr(il, &
-            i+1)-rr(il,i))-ad(il)*(rr(il,i)-rr(il,i-1)))
-          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
-            i))-ad(il)*(u(il,i)-u(il,i-1)))
-          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
-            i))-ad(il)*(v(il,i)-v(il,i-1)))
-        ELSE ! cvflag_grav
-          fr(il, i) = fr(il, i) + 0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, &
-            i))-ad(il)*(rr(il,i)-rr(il,i-1)))
-          fu(il, i) = fu(il, i) + 0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il, &
-            i))-ad(il)*(u(il,i)-u(il,i-1)))
-          fv(il, i) = fv(il, i) + 0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il, &
-            i))-ad(il)*(v(il,i)-v(il,i-1)))
-        END IF ! cvflag_grav
+        ft(il, i) = ft(il, i) - 0.009*grav*sigd(il) * &
+                     (mp(il,i+1)*t_wake(il,i)*b(il,i)-mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
+        IF (cvflag_ice) THEN
+          ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
+                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
+                                  0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1) * &
+                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
+        ELSE
+          ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
+                                       (t_wake(il,i+1)-t_wake(il,i))*dpinv* &
+            cpinv
+        END IF
+
+        ftd(il, i) = ft(il, i)
+! fin precip
+
+! sature
+        ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
+                     (amp1(il)*(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))*cpinv) - &
+                      ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
+
+
+        IF (iflag_mix==0) THEN
+          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i) + &
+                                    t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
+        END IF
+
+
+
+! sb: on ne fait pas encore la correction permettant de mieux
+! conserver l'eau:
+!JYG: correction permettant de mieux conserver l'eau:
+! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
+        fr(il, i) = sigd(il)*evap(il, i) + 0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i)) - &
+                                                      mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
+        fqd(il, i) = fr(il, i)                                                                     ! precip
+
+        fu(il, i) = 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i)) - &
+                               mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
+        fv(il, i) = 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i)) - &
+                               mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
+
+
+        fr(il, i) = fr(il, i) + 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) - &
+                                                 ad(il)*(rr(il,i)-rr(il,i-1)))
+        fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) - &
+                                                 ad(il)*(u(il,i)-u(il,i-1)))
+        fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) - &
+                                                 ad(il)*(v(il,i)-v(il,i-1)))
 
       END IF ! i
     END DO
 
-    ! AC!      do k=1,ntra
-    ! AC!       do il=1,ncum
-    ! AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
-    ! AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
-    ! AC!         cpinv=1.0/cpn(il,i)
-    ! AC!         if (cvflag_grav) then
-    ! AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
-    ! AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
-    ! AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
-    ! AC!         else
-    ! AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
-    ! AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
-    ! AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
-    ! AC!         endif
-    ! AC!        endif
-    ! AC!       enddo
-    ! AC!      enddo
+!AC!      do k=1,ntra
+!AC!       do il=1,ncum
+!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!         cpinv=1.0/cpn(il,i)
+!AC!         if (cvflag_grav) then
+!AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
+!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
+!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
+!AC!         else
+!AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
+!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
+!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
+!AC!         endif
+!AC!        endif
+!AC!       enddo
+!AC!      enddo
 
     DO k = 1, i - 1
@@ -3246,16 +3199,8 @@
             dpinv = 1.0/(ph(il,i)-ph(il,i+1))
             cpinv = 1.0/cpn(il, i)
-            IF (cvflag_grav) THEN
-              ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(hent(il &
-                ,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k, &
-                i)))*cpinv
-
-
-
-            ELSE
-              ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, k, i)*(hent(il,k,i)- &
-                h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k, &
-                i)))*cpinv
-            END IF !cvflag_grav
+            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
+                 (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
+!
+!
           END IF ! i
         END DO
@@ -3266,22 +3211,11 @@
           dpinv = 1.0/(ph(il,i)-ph(il,i+1))
           cpinv = 1.0/cpn(il, i)
-          IF (cvflag_grav) THEN
-            fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
-              ,i)-awat(il)-rr(il,i))
-            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
-              ,i)-u(il,i))
-            fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
-              ,i)-v(il,i))
-          ELSE ! cvflag_grav
-            fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)- &
-              awat(il)-rr(il,i))
-            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
-              ,i)-u(il,i))
-            fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
-              il,i))
-          END IF ! cvflag_grav
-
-          ! (saturated updrafts resulting from mixing)        ! cld
-          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il)) ! cld
+          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
+                                                       (qent(il,k,i)-awat(il)-rr(il,i))
+          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
+          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
+
+! (saturated updrafts resulting from mixing)                                   ! cld
+          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
           nqcond(il, i) = nqcond(il, i) + 1. ! cld
         END IF ! i
@@ -3289,21 +3223,21 @@
     END DO
 
-    ! AC!      do j=1,ntra
-    ! AC!       do k=1,i-1
-    ! AC!        do il=1,ncum
-    ! AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
-    ! AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
-    ! AC!          cpinv=1.0/cpn(il,i)
-    ! AC!          if (cvflag_grav) then
-    ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
-    ! AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
-    ! AC!          else
-    ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
-    ! AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
-    ! AC!          endif
-    ! AC!         endif
-    ! AC!        enddo
-    ! AC!       enddo
-    ! AC!      enddo
+!AC!      do j=1,ntra
+!AC!       do k=1,i-1
+!AC!        do il=1,ncum
+!AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!          cpinv=1.0/cpn(il,i)
+!AC!          if (cvflag_grav) then
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
+!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          else
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
+!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          endif
+!AC!         endif
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
 
     DO k = i, nl + 1
@@ -3314,14 +3248,8 @@
             dpinv = 1.0/(ph(il,i)-ph(il,i+1))
             cpinv = 1.0/cpn(il, i)
-            IF (cvflag_grav) THEN
-              ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(hent(il &
-                ,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k, &
-                i)))*cpinv
-
-
-            ELSE
-              ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, k, i)*(hent(il,k,i)- &
-                h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k,i)))*cpinv
-            END IF !cvflag_grav
+            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
+                  (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k,i)))*cpinv
+
+
           END IF ! i
         END DO
@@ -3333,49 +3261,37 @@
           cpinv = 1.0/cpn(il, i)
 
-          IF (cvflag_grav) THEN
-            fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k &
-              ,i)-rr(il,i))
-            fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k &
-              ,i)-u(il,i))
-            fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k &
-              ,i)-v(il,i))
-          ELSE ! cvflag_grav
-            fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)-rr &
-              (il,i))
-            fu(il, i) = fu(il, i) + 0.1*dpinv*ment(il, k, i)*(uent(il,k,i)-u( &
-              il,i))
-            fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( &
-              il,i))
-          END IF ! cvflag_grav
+          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-rr(il,i))
+          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
+          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
         END IF ! i and k
       END DO
     END DO
 
-    ! AC!      do j=1,ntra
-    ! AC!       do k=i,nl+1
-    ! AC!        do il=1,ncum
-    ! AC!         if (i.le.inb(il) .and. k.le.inb(il)
-    ! AC!     $                .and. iflag(il) .le. 1) then
-    ! AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
-    ! AC!          cpinv=1.0/cpn(il,i)
-    ! AC!          if (cvflag_grav) then
-    ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
-    ! AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
-    ! AC!          else
-    ! AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
-    ! AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
-    ! AC!          endif
-    ! AC!         endif ! i and k
-    ! AC!        enddo
-    ! AC!       enddo
-    ! AC!      enddo
-
-    ! sb: interface with the cloud parameterization:          ! cld
+!AC!      do j=1,ntra
+!AC!       do k=i,nl+1
+!AC!        do il=1,ncum
+!AC!         if (i.le.inb(il) .and. k.le.inb(il)
+!AC!     $                .and. iflag(il) .le. 1) then
+!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!          cpinv=1.0/cpn(il,i)
+!AC!          if (cvflag_grav) then
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
+!AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          else
+!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
+!AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
+!AC!          endif
+!AC!         endif ! i and k
+!AC!        enddo
+!AC!       enddo
+!AC!      enddo
+
+! sb: interface with the cloud parameterization:                               ! cld
 
     DO k = i + 1, nl
       DO il = 1, ncum
-        IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN ! cld
-          ! (saturated downdrafts resulting from mixing)            ! cld
-          qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld
+        IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN               ! cld
+! (saturated downdrafts resulting from mixing)                                 ! cld
+          qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
           nqcond(il, i) = nqcond(il, i) + 1. ! cld
         END IF ! cld
@@ -3383,130 +3299,117 @@
     END DO ! cld
 
-    ! (particular case: no detraining level is found)         ! cld
-    DO il = 1, ncum ! cld
-      IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN ! cld
-        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld
-        nqcond(il, i) = nqcond(il, i) + 1. ! cld
-      END IF ! cld
-    END DO ! cld
-
-    DO il = 1, ncum ! cld
-      IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN ! cld
-        qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld
-      END IF ! cld
-    END DO
-
-    ! AC!      do j=1,ntra
-    ! AC!       do il=1,ncum
-    ! AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
-    ! AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
-    ! AC!         cpinv=1.0/cpn(il,i)
-    ! AC!
-    ! AC!         if (cvflag_grav) then
-    ! AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
-    ! AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
-    ! AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
-    ! AC!         else
-    ! AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
-    ! AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
-    ! AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
-    ! AC!         endif
-    ! AC!        endif ! i
-    ! AC!       enddo
-    ! AC!      enddo
+! (particular case: no detraining level is found)                              ! cld
+    DO il = 1, ncum                                                            ! cld
+      IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
+        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
+        nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
+      END IF                                                                   ! cld
+    END DO                                                                     ! cld
+
+    DO il = 1, ncum                                                            ! cld
+      IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN            ! cld
+        qcond(il, i) = qcond(il, i)/nqcond(il, i)                              ! cld
+      END IF                                                                   ! cld
+    END DO
+
+!AC!      do j=1,ntra
+!AC!       do il=1,ncum
+!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
+!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
+!AC!         cpinv=1.0/cpn(il,i)
+!AC!
+!AC!         if (cvflag_grav) then
+!AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
+!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
+!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
+!AC!         else
+!AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
+!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
+!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
+!AC!         endif
+!AC!        endif ! i
+!AC!       enddo
+!AC!      enddo
 
 
 500 END DO
 
-
-  ! ***   move the detrainment at level inb down to level inb-1   ***
-  ! ***        in such a way as to preserve the vertically        ***
-  ! ***          integrated enthalpy and water tendencies         ***
-
-  ! Correction bug le 18-03-09
+!JYG<
+!Conservation de l'eau
+!   sumdq = 0.
+!   DO k = 1, nl
+!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
+!   END DO
+!   PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
+!JYG>
+! ***   move the detrainment at level inb down to level inb-1   ***
+! ***        in such a way as to preserve the vertically        ***
+! ***          integrated enthalpy and water tendencies         ***
+
+! Correction bug le 18-03-09
   DO il = 1, ncum
     IF (iflag(il)<=1) THEN
-      IF (cvflag_grav) THEN
-        ax = 0.01*grav*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il &
-          ))+t(il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
-          inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
-        ft(il, inb(il)) = ft(il, inb(il)) - ax
-        ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il, &
-          inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il, &
-          inb(il)-1)-ph(il,inb(il))))
-
-        bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))- &
-          rr(il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
-        fr(il, inb(il)) = fr(il, inb(il)) - bx
-        fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb( &
-          il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
-
-        cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u &
-          (il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
-        fu(il, inb(il)) = fu(il, inb(il)) - cx
-        fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb( &
-          il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
-
-        dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v &
-          (il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
-        fv(il, inb(il)) = fv(il, inb(il)) - dx
-        fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb( &
-          il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
-      ELSE
-        ax = 0.1*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t( &
-          il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), &
-          inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
-        ft(il, inb(il)) = ft(il, inb(il)) - ax
-        ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il, &
-          inb(il))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il, &
-          inb(il)-1)-ph(il,inb(il))))
-
-        bx = 0.1*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il, &
-          inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
-        fr(il, inb(il)) = fr(il, inb(il)) - bx
-        fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb( &
-          il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
-
-        cx = 0.1*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il, &
-          inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
-        fu(il, inb(il)) = fu(il, inb(il)) - cx
-        fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb( &
-          il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
-
-        dx = 0.1*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il, &
-          inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1))
-        fv(il, inb(il)) = fv(il, inb(il)) - dx
-        fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb( &
-          il)+1))/(ph(il,inb(il)-1)-ph(il,inb(il)))
-      END IF
+      ax = 0.01*grav*ment(il, inb(il), inb(il))* &
+           (hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
+                                (cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
+      ft(il, inb(il)) = ft(il, inb(il)) - ax
+      ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
+                              (cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))
+
+      bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
+                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
+      fr(il, inb(il)) = fr(il, inb(il)) - bx
+      fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
+                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
+
+      cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
+                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
+      fu(il, inb(il)) = fu(il, inb(il)) - cx
+      fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
+                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
+
+      dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
+                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
+      fv(il, inb(il)) = fv(il, inb(il)) - dx
+      fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
+                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
     END IF !iflag
   END DO
 
-  ! AC!      do j=1,ntra
-  ! AC!       do il=1,ncum
-  ! AC!        IF (iflag(il) .le. 1) THEN
-  ! AC!	IF (cvflag_grav) then
-  ! AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
-  ! AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
-  ! AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
-  ! AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
-  ! AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
-  ! AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
-  ! AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
-  ! AC!	else
-  ! AC!        ex=0.1*ment(il,inb(il),inb(il))
-  ! AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
-  ! AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
-  ! AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
-  ! AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
-  ! AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
-  ! AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
-  ! AC!        ENDIF   !cvflag grav
-  ! AC!        ENDIF    !iflag
-  ! AC!       enddo
-  ! AC!      enddo
-
-
-  ! ***    homogenize tendencies below cloud base    ***
+!JYG<
+!Conservation de l'eau
+!   sumdq = 0.
+!   DO k = 1, nl
+!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
+!   END DO
+!   PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
+!JYG>
+
+!AC!      do j=1,ntra
+!AC!       do il=1,ncum
+!AC!        IF (iflag(il) .le. 1) THEN
+!AC!	IF (cvflag_grav) then
+!AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
+!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
+!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
+!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
+!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
+!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
+!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
+!AC!	else
+!AC!        ex=0.1*ment(il,inb(il),inb(il))
+!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
+!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
+!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
+!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
+!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
+!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
+!AC!        ENDIF   !cvflag grav
+!AC!        ENDIF    !iflag
+!AC!       enddo
+!AC!      enddo
+
+
+! ***    homogenize tendencies below cloud base    ***
 
 
@@ -3522,33 +3425,49 @@
   END DO
 
-  ! do i=1,nl
-  ! do il=1,ncum
-  ! th_wake(il,i)=t_wake(il,i)*(1000.0/p(il,i))**rdcp
-  ! enddo
-  ! enddo
+!do i=1,nl
+!do il=1,ncum
+!th_wake(il,i)=t_wake(il,i)*(1000.0/p(il,i))**rdcp
+!enddo
+!enddo
 
   DO i = 1, nl
     DO il = 1, ncum
       IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
-        ! jyg  Saturated part : use T profile
+!jyg  Saturated part : use T profile
         asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
-        bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(lv(il,i)+(cl-cpd)*(t(il, &
-          i)-t(il,1)))*(ph(il,i)-ph(il,i+1))
-        csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, &
-          1)))*(ph(il,i)-ph(il,i+1))
+!jyg<20140311
+!Correction pour conserver l eau
+        IF (ok_conserv_q) THEN
+          bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
+          csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))
+
+        ELSE
+          bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
+                            (ph(il,i)-ph(il,i+1))
+          csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
+                            (ph(il,i)-ph(il,i+1))
+        ENDIF ! (ok_conserv_q)
+!jyg>
         dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
-        ! jyg  Unsaturated part : use T_wake profile
+!jyg  Unsaturated part : use T_wake profile
         esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
-        fsum(il) = fsum(il) + fqd(il, i)*(lv(il,i)+(cl-cpd)*(t_wake(il, &
-          i)-t_wake(il,1)))*(ph(il,i)-ph(il,i+1))
-        gsum(il) = gsum(il) + (lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il, &
-          1)))*(ph(il,i)-ph(il,i+1))
-        hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, &
-          i)
+!jyg<20140311
+!Correction pour conserver l eau
+        IF (ok_conserv_q) THEN
+          fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
+          gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
+        ELSE
+          fsum(il)=fsum(il)+fqd(il,i)*(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
+                            (ph(il,i)-ph(il,i+1))
+          gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
+                            (ph(il,i)-ph(il,i+1))
+        ENDIF ! (ok_conserv_q)
+!jyg>
+        hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
       END IF
     END DO
   END DO
 
-  ! !!!      do 700 i=1,icb(il)-1
+!!!!      do 700 i=1,icb(il)-1
   DO i = 1, nl
     DO il = 1, ncum
@@ -3562,13 +3481,21 @@
   END DO
 
-
-  ! ***   Check that moisture stays positive. If not, scale tendencies
-  ! in order to ensure moisture positivity
+!jyg<
+!Conservation de l'eau
+!!  sumdq = 0.
+!!  DO k = 1, nl
+!!    sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
+!!  END DO
+!!  PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
+!jyg>
+
+
+! ***   Check that moisture stays positive. If not, scale tendencies
+! in order to ensure moisture positivity
   DO il = 1, ncum
     alpha_qpos(il) = 1.
     IF (iflag(il)<=1) THEN
       IF (fr(il,1)<=0.) THEN
-        alpha_qpos(il) = max(alpha_qpos(il), (-delt*fr(il, &
-          1))/(s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
+        alpha_qpos(il) = max(alpha_qpos(il), (-delt*fr(il,1))/(s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
       END IF
     END IF
@@ -3578,8 +3505,6 @@
       IF (iflag(il)<=1) THEN
         IF (fr(il,i)<=0.) THEN
-          alpha_qpos1(il) = max(1., (-delt*fr(il,i))/(s_wake(il)*rr_wake(il, &
-            i)+(1.-s_wake(il))*rr(il,i)))
-          IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) &
-            = alpha_qpos1(il)
+          alpha_qpos1(il) = max(1., (-delt*fr(il,i))/(s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
+          IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
         END IF
       END IF
@@ -3608,5 +3533,5 @@
         m(il, i) = m(il, i)/alpha_qpos(il)
         mp(il, i) = mp(il, i)/alpha_qpos(il)
-        vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
+        Vprecip(il, i) = vprecip(il, i)/alpha_qpos(il)
       END IF
     END DO
@@ -3622,16 +3547,16 @@
   END DO
 
-  ! AC!      DO j = 1,ntra
-  ! AC!      DO i = 1,nl
-  ! AC!       DO il = 1,ncum
-  ! AC!        IF (iflag(il) .le. 1) THEN
-  ! AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
-  ! AC!        ENDIF
-  ! AC!       ENDDO
-  ! AC!      ENDDO
-  ! AC!      ENDDO
-
-
-  ! ***           reset counter and return           ***
+!AC!      DO j = 1,ntra
+!AC!      DO i = 1,nl
+!AC!       DO il = 1,ncum
+!AC!        IF (iflag(il) .le. 1) THEN
+!AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
+!AC!        ENDIF
+!AC!       ENDDO
+!AC!      ENDDO
+!AC!      ENDDO
+
+
+! ***           reset counter and return           ***
 
   DO il = 1, ncum
@@ -3702,5 +3627,5 @@
           END IF
         END IF
-        ! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
+! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
       END DO
     END DO
@@ -3710,12 +3635,10 @@
     DO k = i, nl
       DO il = 1, ncum
-        ! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il))
-        ! then
+! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
         IF (i<=inb(il) .AND. k<=inb(il)) THEN
           upwd(il, i) = upwd(il, i) + up1(il, k, i)
           dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
         END IF
-        ! c         print
-        ! *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
+! c         print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
       END DO
     END DO
@@ -3723,27 +3646,27 @@
 
 
-  ! !!!      DO il=1,ncum
-  ! !!!      do i=icb(il),inb(il)
-  ! !!!
-  ! !!!      upwd(il,i)=0.0
-  ! !!!      dnwd(il,i)=0.0
-  ! !!!      do k=i,inb(il)
-  ! !!!      up1=0.0
-  ! !!!      dn1=0.0
-  ! !!!      do n=1,i-1
-  ! !!!      up1=up1+ment(il,n,k)
-  ! !!!      dn1=dn1-ment(il,k,n)
-  ! !!!      enddo
-  ! !!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
-  ! !!!      dnwd(il,i)=dnwd(il,i)+dn1
-  ! !!!      enddo
-  ! !!!      enddo
-  ! !!!
-  ! !!!      ENDDO
-
-  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-  ! determination de la variation de flux ascendant entre
-  ! deux niveau non dilue mip
-  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+!!!!      DO il=1,ncum
+!!!!      do i=icb(il),inb(il)
+!!!!
+!!!!      upwd(il,i)=0.0
+!!!!      dnwd(il,i)=0.0
+!!!!      do k=i,inb(il)
+!!!!      up1=0.0
+!!!!      dn1=0.0
+!!!!      do n=1,i-1
+!!!!      up1=up1+ment(il,n,k)
+!!!!      dn1=dn1-ment(il,k,n)
+!!!!      enddo
+!!!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
+!!!!      dnwd(il,i)=dnwd(il,i)+dn1
+!!!!      enddo
+!!!!      enddo
+!!!!
+!!!!      ENDDO
+
+! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! determination de la variation de flux ascendant entre
+! deux niveau non dilue mip
+! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
   DO i = 1, nl
@@ -3787,8 +3710,8 @@
   END DO
 
-  ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-  ! icb represente de niveau ou se trouve la
-  ! base du nuage , et inb le top du nuage
-  ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! icb represente de niveau ou se trouve la
+! base du nuage , et inb le top du nuage
+! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
   DO i = 1, nd
@@ -3800,6 +3723,5 @@
   DO i = 1, nd
     DO il = 1, ncum
-      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il, &
-        i))+rr(il,i)*cpv)
+      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il,i))+rr(il,i)*cpv)
       tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp
       tps(il, i) = tp(il, i)
@@ -3808,77 +3730,79 @@
 
 
-  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
-  ! ***           of condensed water         ***            ! cld
-  ! ! cld
-
-  DO i = 1, nd ! cld
-    DO il = 1, ncum ! cld
-      mac(il, i) = 0.0 ! cld
-      wa(il, i) = 0.0 ! cld
-      siga(il, i) = 0.0 ! cld
-      sax(il, i) = 0.0 ! cld
-    END DO ! cld
-  END DO ! cld
-
-  DO i = minorig, nl ! cld
-    DO k = i + 1, nl + 1 ! cld
-      DO il = 1, ncum ! cld
+! *** diagnose the in-cloud mixing ratio   ***                       ! cld
+! ***           of condensed water         ***                       ! cld
+!! cld                                                               
+                                                                     
+  DO i = 1, nd                                                       ! cld
+    DO il = 1, ncum                                                  ! cld
+      mac(il, i) = 0.0                                               ! cld
+      wa(il, i) = 0.0                                                ! cld
+      siga(il, i) = 0.0                                              ! cld
+      sax(il, i) = 0.0                                               ! cld
+    END DO                                                           ! cld
+  END DO                                                             ! cld
+                                                                     
+  DO i = minorig, nl                                                 ! cld
+    DO k = i + 1, nl + 1                                             ! cld
+      DO il = 1, ncum                                                ! cld
         IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
-          mac(il, i) = mac(il, i) + m(il, k) ! cld
-        END IF ! cld
-      END DO ! cld
-    END DO ! cld
-  END DO ! cld
-
-  DO i = 1, nl ! cld
-    DO j = 1, i ! cld
-      DO il = 1, ncum ! cld
-        IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
-            .AND. j>=icb(il) .AND. iflag(il)<=1) THEN ! cld
-          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld
-            *(ph(il,j)-ph(il,j+1))/p(il, j) ! cld
-        END IF ! cld
-      END DO ! cld
-    END DO ! cld
-  END DO ! cld
-
-  DO i = 1, nl ! cld
-    DO il = 1, ncum ! cld
-      IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld
-          .AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld
-        wa(il, i) = sqrt(2.*sax(il,i)) ! cld
-      END IF ! cld
-    END DO ! cld
-  END DO ! cld
-
-  DO i = 1, nl ! cld
-    DO il = 1, ncum ! cld
-      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld
-        siga(il, i) = mac(il, i)/wa(il, i) & ! cld
-        *rrd*tvp(il, i)/p(il, i)/100./delta ! cld
-      siga(il, i) = min(siga(il,i), 1.0) ! cld
-      ! IM cf. FH
-      IF (iflag_clw==0) THEN
-        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) & ! cld
-          +(1.-siga(il,i))*qcond(il, i) ! cld
-      ELSE IF (iflag_clw==1) THEN
-        qcondc(il, i) = qcond(il, i) ! cld
-      END IF
-
-    END DO ! cld
-  END DO
-  ! print*,'cv3_yield fin'
-    ! cld
+          mac(il, i) = mac(il, i) + m(il, k)                         ! cld
+        END IF                                                       ! cld
+      END DO                                                         ! cld
+    END DO                                                           ! cld
+  END DO                                                             ! cld
+
+  DO i = 1, nl                                                       ! cld
+    DO j = 1, i                                                      ! cld
+      DO il = 1, ncum                                                ! cld
+        IF (i>=icb(il) .AND. i<=(inb(il)-1) &                        ! cld
+            .AND. j>=icb(il) .AND. iflag(il)<=1) THEN                ! cld
+          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) &       ! cld
+            *(ph(il,j)-ph(il,j+1))/p(il, j)                          ! cld
+        END IF                                                       ! cld
+      END DO                                                         ! cld
+    END DO                                                           ! cld
+  END DO                                                             ! cld
+
+  DO i = 1, nl                                                       ! cld
+    DO il = 1, ncum                                                  ! cld
+      IF (i>=icb(il) .AND. i<=(inb(il)-1) &                          ! cld
+          .AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN               ! cld
+        wa(il, i) = sqrt(2.*sax(il,i))                               ! cld
+      END IF                                                         ! cld
+    END DO                                                           ! cld
+  END DO                                                             ! cld
+
+  DO i = 1, nl                                                       ! cld
+    DO il = 1, ncum                                                  ! cld
+      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
+        siga(il, i) = mac(il, i)/wa(il, i) &                         ! cld
+        *rrd*tvp(il, i)/p(il, i)/100./delta                          ! cld
+      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
+! IM cf. FH                                                          
+      IF (iflag_clw==0) THEN                                         ! cld
+        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) &       ! cld
+          +(1.-siga(il,i))*qcond(il, i)                              ! cld
+      ELSE IF (iflag_clw==1) THEN                                    ! cld
+        qcondc(il, i) = qcond(il, i)                                 ! cld
+      END IF                                                         ! cld
+
+    END DO                                                           ! cld
+  END DO
+! print*,'cv3_yield fin'
+
   RETURN
 END SUBROUTINE cv3_yield
 
-! AC! et !RomP >>>
-SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, ment, sigij, da, phi, phi2, &
-    d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb)
+!AC! et !RomP >>>
+SUBROUTINE cv3_tracer(nloc, len, ncum, nd, na, &
+                      ment, sigij, da, phi, phi2, d1a, dam, &
+                      ep, Vprecip, elij, clw, epmlmMm, eplaMm, &
+                      icb, inb)
   IMPLICIT NONE
 
   include "cv3param.h"
 
-  ! inputs:
+!inputs:
   INTEGER ncum, nd, na, nloc, len
   REAL ment(nloc, na, na), sigij(nloc, na, na)
@@ -3886,23 +3810,23 @@
   REAL ep(nloc, na)
   INTEGER icb(nloc), inb(nloc)
-  REAL vprecip(nloc, nd+1)
-  ! ouputs:
+  REAL Vprecip(nloc, nd+1)
+!ouputs:
   REAL da(nloc, na), phi(nloc, na, na)
   REAL phi2(nloc, na, na)
   REAL d1a(nloc, na), dam(nloc, na)
-  REAL epmlmmm(nloc, na, na), eplamm(nloc, na)
-  ! variables pour tracer dans precip de l'AA et des mel
-  ! local variables:
+  REAL epmlmMm(nloc, na, na), eplaMm(nloc, na)
+! variables pour tracer dans precip de l'AA et des mel
+!local variables:
   INTEGER i, j, k
   REAL epm(nloc, na, na)
 
-  ! variables d'Emanuel : du second indice au troisieme
-  ! --->    tab(i,k,j) -> de l origine k a l arrivee j
-  ! ment, sigij, elij
-  ! variables personnelles : du troisieme au second indice
-  ! --->    tab(i,j,k) -> de k a j
-  ! phi, phi2
-
-  ! initialisations
+! variables d'Emanuel : du second indice au troisieme
+! --->    tab(i,k,j) -> de l origine k a l arrivee j
+! ment, sigij, elij
+! variables personnelles : du troisieme au second indice
+! --->    tab(i,j,k) -> de k a j
+! phi, phi2
+
+! initialisations
 
   da(:, :) = 0.
@@ -3910,20 +3834,20 @@
   dam(:, :) = 0.
   epm(:, :, :) = 0.
-  eplamm(:, :) = 0.
-  epmlmmm(:, :, :) = 0.
+  eplaMm(:, :) = 0.
+  epmlmMm(:, :, :) = 0.
   phi(:, :, :) = 0.
   phi2(:, :, :) = 0.
 
-  ! fraction deau condensee dans les melanges convertie en precip : epm
-  ! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
+! fraction deau condensee dans les melanges convertie en precip : epm
+! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
   DO j = 1, na
     DO k = 1, na
       DO i = 1, ncum
-        IF (k>=icb(i) .AND. k<=inb(i) .AND. & ! !jyg     &
-                                              ! j.ge.k.and.j.le.inb(i)) then
-          ! !jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
+        IF (k>=icb(i) .AND. k<=inb(i) .AND. & 
+!!jyg              j.ge.k.and.j.le.inb(i)) then
+!!jyg             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
             j>k .AND. j<=inb(i)) THEN
           epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16)
-          ! !
+!!
           epm(i, j, k) = max(epm(i,j,k), 0.0)
         END IF
@@ -3937,6 +3861,6 @@
       DO i = 1, ncum
         IF (k>=icb(i) .AND. k<=inb(i)) THEN
-          eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- &
-            sigij(i,j,k))
+          eplaMm(i, j) = eplamm(i, j) + &
+                         ep(i, j)*clw(i, j)*ment(i, j, k)*(1.-sigij(i,j,k))
         END IF
       END DO
@@ -3948,5 +3872,5 @@
       DO i = 1, ncum
         IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN
-          epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
+          epmlmMm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j)
         END IF
       END DO
@@ -3954,5 +3878,5 @@
   END DO
 
-  ! matrices pour calculer la tendance des concentrations dans cvltr.F90
+! matrices pour calculer la tendance des concentrations dans cvltr.F90
   DO j = 1, na
     DO k = 1, na
@@ -3962,7 +3886,5 @@
         d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sigij(i,k,j))
         IF (k<=j) THEN
-          dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, k, j)*(1.-ep(i,k))*(1. &
-            -sigij(i,k,j))
-
+          dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, k, j)*(1.-ep(i,k))*(1.-sigij(i,k,j))
           phi2(i, j, k) = phi(i, j, k)*epm(i, j, k)
         END IF
@@ -3973,15 +3895,20 @@
   RETURN
 END SUBROUTINE cv3_tracer
-! AC! et !RomP <<<
-
-SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, &
-    sig, w0, ft, fq, fu, fv, ftra, ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
-    iflag1, precip1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, ma1, upwd1, dnwd1, &
-    dnwd01, qcondc1, wd1, cape1)
+!AC! et !RomP <<<
+
+SUBROUTINE cv3_uncompress(nloc, len, ncum, nd, ntra, idcum, &
+                          iflag, &
+                          precip, sig, w0, &
+                          ft, fq, fu, fv, ftra, &
+                          Ma, upwd, dnwd, dnwd0, qcondc, wd, cape, &
+                          iflag1, &
+                          precip1, sig1, w01, &
+                          ft1, fq1, fu1, fv1, ftra1, &
+                          Ma1, upwd1, dnwd1, dnwd01, qcondc1, wd1, cape1)
   IMPLICIT NONE
 
   include "cv3param.h"
 
-  ! inputs:
+!inputs:
   INTEGER len, ncum, nd, ntra, nloc
   INTEGER idcum(nloc)
@@ -3996,5 +3923,5 @@
   REAL wd(nloc), cape(nloc)
 
-  ! outputs:
+!outputs:
   INTEGER iflag1(len)
   REAL precip1(len)
@@ -4007,5 +3934,5 @@
   REAL wd1(nloc), cape1(nloc)
 
-  ! local variables:
+!local variables:
   INTEGER i, k, j
 
@@ -4038,12 +3965,13 @@
 
 
-  ! AC!        do 2100 j=1,ntra
-  ! AC!c oct3         do 2110 k=1,nl
-  ! AC!         do 2110 k=1,nd ! oct3
-  ! AC!          do 2120 i=1,ncum
-  ! AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
-  ! AC! 2120     continue
-  ! AC! 2110    continue
-  ! AC! 2100   continue
+!AC!        do 2100 j=1,ntra
+!AC!c oct3         do 2110 k=1,nl
+!AC!         do 2110 k=1,nd ! oct3
+!AC!          do 2120 i=1,ncum
+!AC!            ftra1(idcum(i),k,j)=ftra(i,k,j)
+!AC! 2120     continue
+!AC! 2110    continue
+!AC! 2100   continue
+!
   RETURN
 END SUBROUTINE cv3_uncompress
