Index: /LMDZ6/trunk/libf/phylmd/concvl.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/concvl.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/concvl.F90	(revision 3496)
@@ -13,8 +13,8 @@
 !RomP >>>
 !!     .             da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
-                  da, phi, mp, phii, d1a, dam, sij, clw, elij, &     ! RomP
+                  da, phi, mp, phii, d1a, dam, sij, qta, clw, elij, &! RomP
                   dd_t, dd_q, lalim_conv, wght_th, &                 ! RomP
                   evap, ep, epmlmMm, eplaMm, &                       ! RomP
-                  wdtrainA, wdtrainM, wght, qtc, sigt, &
+                  wdtrainA, wdtrainS, wdtrainM, wght, qtc, sigt, &
                   tau_cld_cv, coefw_cld_cv, &                           ! RomP+RL, AJ
 !RomP <<<
@@ -80,4 +80,5 @@
 ! eplaMm-----output-R
 ! wdtrainA---output-R
+! wdtrainS---output-R
 ! wdtrainM---output-R
 ! wght-------output-R
@@ -134,4 +135,5 @@
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: d1a, dam
   REAL, DIMENSION(klon,klev,klev),INTENT(OUT)   :: sij, elij
+  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: qta
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: clw
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: dd_t, dd_q
@@ -139,5 +141,5 @@
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: eplaMm
   REAL, DIMENSION(klon,klev,klev), INTENT(OUT)  :: epmlmMm
-  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: wdtrainA, wdtrainM
+  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: wdtrainA, wdtrainS, wdtrainM
 ! RomP <<<
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: wght                       !RL
@@ -437,6 +439,6 @@
 !!                   evap,ep,epmlmMm,eplaMm,                              ! RomP
                     da, phi, mp, phii, d1a, dam, sij, wght, &           ! RomP+RL
-                    clw, elij, evap, ep, epmlmMm, eplaMm, &             ! RomP+RL
-                    wdtrainA, wdtrainM, qtc, sigt, &
+                    qta, clw, elij, evap, ep, epmlmMm, eplaMm, &        ! RomP+RL
+                    wdtrainA, wdtrainS, wdtrainM, qtc, sigt, &
                     tau_cld_cv, coefw_cld_cv, &                         ! RomP,AJ
 !AC!+!RomP+jyg
Index: /LMDZ6/trunk/libf/phylmd/cv3_mixscale.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cv3_mixscale.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cv3_mixscale.F90	(revision 3496)
@@ -13,15 +13,20 @@
   include "cv3param.h"
 
-  INTEGER nloc, ncum, na
+!inputs:
+  INTEGER, INTENT (IN)                               :: ncum, na, nloc
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: m
+!input/outputs:
+  REAL, DIMENSION (nloc, na, na), INTENT (INOUT)     :: ment
+
+!local variables:
   INTEGER i, j, il
-  REAL ment(nloc, na, na), m(nloc, na)
 
-  DO j = 1, nl
-    DO i = 1, nl
-      DO il = 1, ncum
-        ment(il, i, j) = m(il, i)*ment(il, i, j)
+    DO j = 1, nl
+      DO i = 1, nl
+        DO il = 1, ncum
+          ment(il, i, j) = m(il, i)*ment(il, i, j)
+        END DO
       END DO
     END DO
-  END DO
 
 
Index: /LMDZ6/trunk/libf/phylmd/cv3_routines.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cv3_routines.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cv3_routines.F90	(revision 3496)
@@ -35,4 +35,5 @@
 
   include "cv3param.h"
+  include "cvflag.h"
   include "conema3.h"
 
@@ -125,6 +126,14 @@
      tlcrit=-55.0
      CALL getin_p('tlcrit',tlcrit)
+     ejectliq=0.
+     CALL getin_p('ejectliq',ejectliq)
+     ejectice=0.
+     CALL getin_p('ejectice',ejectice)
+     cvflag_prec_eject = .FALSE.
+     CALL getin_p('cvflag_prec_eject',cvflag_prec_eject)
      qsat_depends_on_qt = .FALSE.
      CALL getin_p('qsat_depends_on_qt',qsat_depends_on_qt)
+     adiab_ascent_mass_flux_depends_on_ejectliq = .FALSE.
+     CALL getin_p('adiab_ascent_mass_flux_depends_on_ejectliq',adiab_ascent_mass_flux_depends_on_ejectliq)
 
     WRITE (*, *) 't_top_max=', t_top_max
@@ -172,4 +181,5 @@
 
   include "cv3param.h"
+  include "cvflag.h"
 
 !inputs:
@@ -292,4 +302,5 @@
   USE mod_phys_lmdz_transfert_para, ONLY : bcast
   USE add_phys_tend_mod, ONLY: fl_cor_ebil
+  USE print_control_mod, ONLY: prt_level
   IMPLICIT NONE
 
@@ -519,4 +530,8 @@
     END DO
   ENDIF
+  IF (prt_level .GE. 10) THEN
+    print *,'cv3_feed : iflag(1), pfeed(1), plcl(1), wghti(1,k) ', &
+                        iflag(1), pfeed(1), plcl(1), (wghti(1,k),k=1,10)
+  ENDIF
 
 ! -------------------------------------------------------------------
@@ -1108,5 +1123,6 @@
                          tnk, qnk, gznk, hnk, t, q, qs, gz, &
                          p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
-                         inb, tp, tvp, clw, hp, ep, sigp, buoy, frac)
+                         inb, tp, tvp, clw, hp, ep, sigp, buoy, &
+                         frac_a, frac_s, qpreca, qta)
   USE print_control_mod, ONLY: prt_level
   IMPLICIT NONE
@@ -1156,5 +1172,7 @@
   REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ep, sigp, hp
   REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: buoy
-  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac   ! ice fraction
+  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac_a, frac_s
+  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qpreca
+  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qta
 
 !local variables:
@@ -1162,12 +1180,22 @@
   REAL smallestreal
   REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
+  REAL                                               :: phinu2p
   REAL als
   REAL                                               :: qsat_new, snew
   REAL, DIMENSION (nloc,nd)                          :: qi
+  REAL, DIMENSION (nloc,nd)                          :: ha    ! moist static energy of adiabatic ascents 
+                                                              ! taking into account precip ejection
+  REAL, DIMENSION (nloc,nd)                          :: hla   ! liquid water static energy of adiabatic ascents 
+                                                              ! taking into account precip ejection
+  REAL, DIMENSION (nloc,nd)                          :: qcld  ! specific cloud water
   REAL, DIMENSION (nloc,nd)                          :: qhsat    ! specific humidity at saturation
   REAL, DIMENSION (nloc,nd)                          :: dqhsatdT ! dqhsat/dT
+  REAL, DIMENSION (nloc,nd)                          :: frac  ! ice fraction function of envt temperature
+  REAL, DIMENSION (nloc,nd)                          :: qps   ! specific solid precipitation
+  REAL, DIMENSION (nloc,nd)                          :: qpl   ! specific liquid precipitation
   REAL, DIMENSION (nloc)                             :: ah0, cape, capem, byp
   LOGICAL, DIMENSION (nloc)                          :: lcape
   INTEGER, DIMENSION (nloc)                          :: iposit
+  REAL                                               :: denomm1
   REAL                                               :: by, defrac, pden, tbis
   REAL                                               :: fracg
@@ -1196,4 +1224,5 @@
   END DO
 
+
 ! =====================================================================
 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
@@ -1211,8 +1240,47 @@
              qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i)
   END DO
+!
+!  Ice fraction
+!
+  IF (cvflag_ice) THEN
+    DO k = minorig, nl
+      DO i = 1, ncum
+          frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
+          frac(i, k) = min(max(frac(i,k),0.0), 1.0)
+      END DO
+    END DO
+! Below cloud base, set ice fraction to cloud base value
+    DO k = 1, nl
+      DO i = 1, ncum
+        IF (k<icb(i)) THEN
+          frac(i,k) = frac(i,icb(i))
+        END IF
+      END DO
+    END DO
+  ELSE
+    DO k = 1, nl
+      DO i = 1, ncum
+          frac(i,k) = 0.
+      END DO
+    END DO
+  ENDIF ! (cvflag_ice)
+
 
   DO k = minorig, nl
     DO i = 1,ncum
+      ha(i,k) = ah0(i)
+      hla(i,k) = hnk(i)
+      qta(i,k) = qnk(i)
+      qpreca(i,k) = 0.
+      frac_a(i,k) = 0.
+      frac_s(i,k) = frac(i,k)
+      qpl(i,k) = 0.
+      qps(i,k) = 0.
       qhsat(i,k) = qs(i,k)
+      qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.)
+      IF (k <= icb(i)+1) THEN
+        qhsat(i,k) = qnk(i)-clw(i,k)
+        qcld(i,k) = clw(i,k)
+      ENDIF 
     ENDDO
   ENDDO
@@ -1347,21 +1415,4 @@
 !----------------------------------------------------------------------------
 !
-!  Ice fraction
-!
-    DO k = minorig + 1, nl
-      DO i = 1, ncum
-          frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
-          frac(i, k) = min(max(frac(i,k),0.0), 1.0)
-      END DO
-    END DO
-! Below cloud base, set ice fraction to cloud base value
-    DO k = 1, nl
-      DO i = 1, ncum
-        IF (k<icb(i)) THEN
-          frac(i,k) = frac(i,icb(i))
-        END IF
-      END DO
-    END DO
-
     DO k = minorig + 1, nl
       DO i = 1,ncum
@@ -1376,5 +1427,5 @@
           IF (k>=(icbs(i)+1)) THEN                                ! convect3
             tg = tp(i, k)
-            IF (tg .gt. Tx) THEN
+            IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
               es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
               qg = eps*es/(p(i,k)-es*(1.-eps))
@@ -1386,11 +1437,11 @@
             ENDIF
             IF (qsat_depends_on_qt) THEN
-              dqgdT = dqgdT*(1.-qnk(i))/(1.-qg)**2
-              qg = qg*(1.-qnk(i))/(1.-qg)            
+              dqgdT = dqgdT*(1.-qta(i,k-1))/(1.-qg)**2
+              qg = qg*(1.-qta(i,k-1))/(1.-qg)            
             ENDIF
-            ahg = (cpd + (cl-cpd)*qnk(i))*tg + lv(i,k)*qg - &
-                  lf(i,k)*frac(i,k)*(qnk(i) - qg) + gz(i,k)
+            ahg = (cpd + (cl-cpd)*qta(i,k-1))*tg + lv(i,k)*qg - &
+                  lf(i,k)*frac(i,k)*(qta(i,k-1) - qg) + gz(i,k)
             Tp(i,k) = tg + (ah0(i) - ahg)/ &
-                    (cpd + (cl-cpd)*qnk(i) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT)
+                    (cpd + (cl-cpd)*qta(i,k-1) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT)
 !!   print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', &
 !!                                 k, Tp(i,k), ah0(i), ahg
@@ -1401,5 +1452,5 @@
         IF (k>=(icbs(i)+1)) THEN                                ! convect3
           tg = tp(i, k)
-          IF (tg .gt. Tx) THEN
+          IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN
             es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
             qg = eps*es/(p(i,k)-es*(1.-eps))
@@ -1408,10 +1459,15 @@
             qg = eps*esi/(p(i,k)-esi*(1.-eps))
           ENDIF
-            IF (qsat_depends_on_qt) THEN
-              qg = qg*(1.-qnk(i))/(1.-qg)            
-            ENDIF
-          clw(i, k) = qnk(i) - qg
+          IF (qsat_depends_on_qt) THEN
+            qg = qg*(1.-qta(i,k-1))/(1.-qg)            
+          ENDIF
+          qhsat(i,k) = qg
+        END IF ! (k>=(icbs(i)+1))
+      END DO ! i = 1, ncum
+      DO i = 1, ncum
+        IF (k>=(icbs(i)+1)) THEN                                ! convect3
+          clw(i, k) = qta(i,k-1) - qhsat(i,k)
           clw(i, k) = max(0.0, clw(i,k))
-          tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
+          tvp(i, k) = max(0., tp(i,k)*(1.+qhsat(i,k)/eps-qta(i,k-1)))
 ! print*,tvp(i,k),'tvp'
           IF (clw(i,k)<1.E-11) THEN
@@ -1421,9 +1477,58 @@
         END IF ! (k>=(icbs(i)+1))
       END DO ! i = 1, ncum
+!
+      IF (cvflag_prec_eject) THEN
+        DO i = 1, ncum
+          IF (k>=(icbs(i)+1)) THEN                                ! convect3
+!  Specific precipitation (liquid and solid) and ice content 
+!  before ejection of precipitation                                                     !!jygprl
+            elacrit = elcrit*min(max(1.-(tp(i,k)-T0)/Tlcrit, 0.), 1.)                   !!jygprl
+!!!!            qcld(i,k) = min(clw(i,k), elacrit)                                          !!jygprl
+            qcld(i,k) = min(clw(i,k), elacrit*(1.-qta(i,k-1))/(1.-elacrit))             !!jygprl
+            phinu2p = qhsat(i,k-1) + qcld(i,k-1) - (qhsat(i,k) + qcld(i,k))                     !!jygprl
+            qpl(i,k) = qpl(i,k-1) + (1.-frac(i,k))*phinu2p                            !!jygprl
+            qps(i,k) = qps(i,k-1) + frac(i,k)     *phinu2p                            !!jygprl
+            qi(i,k) = (1.-ejectliq)*clw(i,k)*frac(i,k) + &                            !!jygprl
+                     ejectliq*(qps(i,k-1) + frac(i,k)*(phinu2p+qcld(i,k)))            !!jygprl
+!!
+!  =====================================================================================
+!  Ejection of precipitation from adiabatic ascents if requested (cvflag_prec_eject=True):
+!  Compute the steps of total water (qta), of moist static energy (ha), of specific 
+!  precipitation (qpl and qps) and of specific cloud water (qcld) associated with precipitation
+!   ejection.
+!  =====================================================================================
+!  
+!   Verif
+            qpreca(i,k) = ejectliq*qpl(i,k) + ejectice*qps(i,k)                                   !!jygprl
+            frac_a(i,k) = ejectice*qps(i,k)/max(qpreca(i,k),smallestreal)                         !!jygprl
+            frac_s(i,k) = (1.-ejectliq)*frac(i,k) + &                                             !!jygprl
+               ejectliq*(1. - (qpl(i,k)+(1.-frac(i,k))*qcld(i,k))/max(clw(i,k),smallestreal))     !!jygprl
+!          
+            denomm1 = 1./(1. - qpreca(i,k))
+!          
+            qta(i,k) = qta(i,k-1) - &
+                      qpreca(i,k)*(1.-qta(i,k-1))*denomm1
+            ha(i,k)  = ha(i,k-1) + &
+                      ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cl-cpd)*tp(i,k) + &
+                                  lv(i,k)*qhsat(i,k) - lf(i,k)*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
+                        lf(i,k)*ejectice*qps(i,k))*denomm1
+            hla(i,k) = hla(i,k-1) + &
+                      ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cpv-cpd)*tp(i,k) - &
+                                  lv(i,k)*((1.-frac_s(i,k))*qcld(i,k)+qpl(i,k)) - &
+                                  (lv(i,k)+lf(i,k))*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + &
+                        lv(i,k)*ejectliq*qpl(i,k) + (lv(i,k)+lf(i,k))*ejectice*qps(i,k))*denomm1
+            qpl(i,k) = qpl(i,k)*(1.-ejectliq)*denomm1
+            qps(i,k) = qps(i,k)*(1.-ejectice)*denomm1
+            qcld(i,k) = qcld(i,k)*denomm1
+            qhsat(i,k) = qhsat(i,k)*(1.-qta(i,k))/(1.-qta(i,k-1))
+         END IF ! (k>=(icbs(i)+1))
+        END DO ! i = 1, ncum
+      ENDIF  ! (cvflag_prec_eject)
+! 
     END DO ! k = minorig + 1, nl
 !
 !----------------------------------------------------------------------------
 !
-  ELSE ! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
+  ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
 !
 !----------------------------------------------------------------------------
@@ -1587,5 +1692,5 @@
 !----------------------------------------------------------------------------
 !
-  ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE
+  ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0)
 !
 !----------------------------------------------------------------------------
@@ -1872,4 +1977,14 @@
   IF (cvflag_ice) THEN
 !
+  IF (cvflag_prec_eject) THEN
+!!    DO k = minorig + 1, nl
+!!      DO i = 1, ncum
+!!        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
+!!          frac_s(i,k) = qi(i,k)/max(clw(i,k),smallestreal)    
+!!          frac_s(i,k) = 1. - (qpl(i,k)+(1.-frac_s(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)    
+!!        END IF
+!!      END DO
+!!    END DO
+  ELSE    ! (cvflag_prec_eject)
     DO k = minorig + 1, nl
       DO i = 1, ncum
@@ -1877,19 +1992,20 @@
 !jyg< frac computation moved to beginning of cv3_undilute2.
 !     kept here for compatibility test with CMip6 version
-          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)
+          frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
+          frac_s(i, k) = min(max(frac_s(i,k),0.0), 1.0)
         END IF
       END DO
     END DO
-! Below cloud base, set ice fraction to cloud base value
-!!    DO k = 1, nl
-!!      DO i = 1, ncum
-!!        IF (k<icb(i)) THEN
-!!          frac(i,k) = frac(i,icb(i))
-!!        END IF
-!!      END DO
-!!    END DO
+  ENDIF  ! (cvflag_prec_eject) ELSE
+    DO k = minorig + 1, nl
+      DO i = 1, ncum
+        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
+!!          hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* &     !!jygprl
+!!                              ep(i, k)*clw(i, k)                                    !!jygprl
+          hp(i, k) = hla(i,k-1) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* &   !!jygprl
+                              ep(i, k)*clw(i, k)                                      !!jygprl
+        END IF
+      END DO
+    END DO
 !
   ELSE   ! (cvflag_ice)
@@ -2576,9 +2692,9 @@
 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, &
+                     th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , &                       !!jygprl
                      m, ment, elij, delt, plcl, coef_clos, &
                      mp, rp, up, vp, trap, wt, water, evap, fondue, ice, &
                      faci, b, sigd, &
-                     wdtrainA, wdtrainM)                                      ! RomP
+                     wdtrainA, wdtrainS, wdtrainM)                                      ! RomP
   USE print_control_mod, ONLY: prt_level, lunout
   IMPLICIT NONE
@@ -2598,7 +2714,12 @@
   REAL, DIMENSION (nloc, na), INTENT (IN)            :: gz
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
-  REAL tra(nloc, nd, ntra)
-  REAL p(nloc, nd), ph(nloc, nd+1)
-  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, sigp, clw
+  REAL, DIMENSION (nloc, nd, ntra), INTENT(IN)       :: tra
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: p
+  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, sigp, clw   !adiab ascent shedding
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac_s          !ice fraction in adiab ascent shedding !!jygprl
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: qpreca          !adiab ascent precip                   !!jygprl
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac_a          !ice fraction in adiab ascent precip   !!jygprl
+  REAL, DIMENSION (nloc, na), INTENT (IN)            :: qta             !adiab ascent specific total water     !!jygprl
   REAL, DIMENSION (nloc, na), INTENT (IN)            :: th, tv, lv, cpn
   REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
@@ -2613,5 +2734,6 @@
   REAL, DIMENSION (nloc, na), INTENT (OUT)           :: mp, rp, up, vp
   REAL, DIMENSION (nloc, na), INTENT (OUT)           :: water, evap, wt
-  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: ice, fondue, faci
+  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: ice, fondue
+  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: faci            ! ice fraction in precipitation
   REAL, DIMENSION (nloc, na, ntra), INTENT (OUT)     :: trap
   REAL, DIMENSION (nloc, na), INTENT (OUT)           :: b
@@ -2621,8 +2743,9 @@
 ! Distinction des wdtrain
 ! Pa = wdtrainA     Pm = wdtrainM
-  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: wdtrainA, wdtrainM
+  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: wdtrainA, wdtrainS, wdtrainM
 
 !local variables
   INTEGER i, j, k, il, num1, ndp1
+  REAL smallestreal
   REAL tinv, delti, coef
   REAL awat, afac, afac1, afac2, bfac
@@ -2631,8 +2754,10 @@
   REAL ampmax, thaw
   REAL tevap(nloc)
-  REAL lvcp(nloc, na), lfcp(nloc, na)
-  REAL h(nloc, na), hm(nloc, na)
-  REAL frac(nloc, na)
-  REAL fraci(nloc, na), prec(nloc, na)
+  REAL, DIMENSION (nloc, na)      :: lvcp, lfcp
+  REAL, DIMENSION (nloc, na)      :: h, hm
+  REAL, DIMENSION (nloc, na)      :: ma
+  REAL, DIMENSION (nloc, na)      :: frac          ! ice fraction in precipitation source
+  REAL, DIMENSION (nloc, na)      :: fraci         ! provisionnal ice fraction in precipitation
+  REAL, DIMENSION (nloc, na)      :: prec
   REAL wdtrain(nloc)
   LOGICAL lwork(nloc), mplus(nloc)
@@ -2641,4 +2766,6 @@
 ! ------------------------------------------------------
 IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
+
+smallestreal=tiny(smallestreal)
 
 ! =============================
@@ -2660,4 +2787,5 @@
 !! RomP >>>
 wdtrainA(:,:) = 0.
+wdtrainS(:,:) = 0.
 wdtrainM(:,:) = 0.
 !! RomP <<<
@@ -2715,4 +2843,39 @@
   END DO
 
+!
+! Get adiabatic ascent mass flux 
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!! Warning : this option leads to water conservation violation
+!!!           Expert only
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    DO il = 1, ncum
+      ma(il, nlp) = 0.
+      ma(il, 1)   = 0.
+    END DO
+
+  DO i = nl, 2, -1
+      DO il = 1, ncum
+        ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i)
+      END DO
+  END DO
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    DO il = 1, ncum
+      ma(il, nlp) = 0.
+      ma(il, 1)   = 0.
+    END DO
+
+  DO i = nl, 2, -1
+      DO il = 1, ncum
+        ma(il, i) = ma(il, i+1) + m(il, i)
+      END DO
+  END DO
+
+  ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -2739,15 +2902,12 @@
 ! ***              calculate detrained precipitation             ***
 
-    DO il = 1, ncum
-      IF (i<=inb(il) .AND. lwork(il)) THEN
-        IF (cvflag_grav) THEN
-          wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)
-          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
-        END IF
-      END IF
-    END DO
+
+    DO il = 1, ncum                                                   
+      IF (i<=inb(il) .AND. lwork(il)) THEN                            
+        wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i)            
+        wdtrainS(il, i) = wdtrain(il)/grav                                            !   Ps   jyg
+!!        wdtrainA(il, i) = wdtrain(il)/grav                                          !   Ps   RomP
+      END IF                                                          
+    END DO                                                            
 
     IF (i>1) THEN
@@ -2757,11 +2917,7 @@
             awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
             awat = max(awat, 0.0)
-            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
-            ELSE
-              wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i)
-              wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i)   !   Pm  RomP
-            END IF
+            wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i)
+            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i)    !   Pm  jyg
+!!            wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i)  !   Pm  RomP
           END IF
         END DO
@@ -2769,4 +2925,35 @@
     END IF
 
+    IF (cvflag_prec_eject) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!! Warning : this option leads to water conservation violation
+!!!           Expert only
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+          IF ( i > 1) THEN
+            DO il = 1, ncum
+              IF (i<=inb(il) .AND. lwork(il)) THEN
+                wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1))    !   Pa   jygprl
+                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
+              END IF
+            END DO
+          ENDIF  ! ( i > 1)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+          IF ( i > 1) THEN
+            DO il = 1, ncum
+              IF (i<=inb(il) .AND. lwork(il)) THEN
+                wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))                        !   Pa   jygprl
+                wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i)
+              END IF
+            END DO
+          ENDIF  ! ( i > 1)
+
+      ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+    ENDIF  ! (cvflag_prec_eject)
+
 
 ! ***    find rain water and evaporation using provisional   ***
@@ -2774,4 +2961,14 @@
 
 
+    IF (cvflag_ice) THEN                                                                                !!jygprl
+      DO il = 1, ncum                                                                                   !!jygprl
+        IF (i<=inb(il) .AND. lwork(il)) THEN                                                            !!jygprl
+          frac(il, i) = (frac_a(il,i)*wdtrainA(il,i)+frac_s(il,i)*(wdtrainS(il,i)+wdtrainM(il,i))) / &  !!jygprl
+                        max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal)                  !!jygprl
+          fraci(il, i) = frac(il, i)                                                                    !!jygprl
+        END IF                                                                                          !!jygprl
+      END DO                                                                                            !!jygprl
+    END IF                                                                                              !!jygprl
+
     DO il = 1, ncum
       IF (i<=inb(il) .AND. lwork(il)) THEN
@@ -2779,11 +2976,4 @@
         wt(il, i) = 45.0
 
-        IF (cvflag_ice) THEN
-          frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)
-          frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)
-          fraci(il, inb(il)) = frac(il, inb(il))
-        ELSE
-          CONTINUE
-        END IF
 
         IF (i<inb(il)) THEN
@@ -2802,6 +2992,6 @@
           rp(il, i) = 0.5*(rp(il,i)+rr(il,i))
         END IF
-        fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
-        fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
+!!        fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)
+!!        fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)
         rp(il, i) = max(rp(il,i), 0.0)
         rp(il, i) = amin1(rp(il,i), rs(il,i))
@@ -3230,5 +3420,5 @@
                      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, &
+                     ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, &
                      wt, water, ice, evap, fondue, faci, b, sigd, &
                      ment, qent, hent, iflag_mix, uent, vent, &
@@ -3240,5 +3430,5 @@
 !!                     tls, tps,                             ! useless . jyg
                      qcondc, wd, &
-                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
+                     ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)
 
     USE print_control_mod, ONLY: lunout, prt_level
@@ -3280,5 +3470,7 @@
       REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN)  :: traent
       REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, wghti
-      REAL,INTENT(IN)                                    :: tau_cld_cv, coefw_cld_cv
+      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
+      REAL, DIMENSION (nloc, na),INTENT(IN)              :: qpreca
+      REAL, INTENT(IN)                                   :: tau_cld_cv, coefw_cld_cv
 !
 !input/output:
@@ -3309,4 +3501,5 @@
       REAL                                               :: ax, bx, cx, dx, ex
       REAL                                               :: cpinv, rdcp, dpinv
+      REAL                                               :: sigaq
       REAL, DIMENSION (nloc)                             ::  awat
       REAL, DIMENSION (nloc, nd)                         :: lvcp, lfcp              ! , mke ! unused . jyg
@@ -3326,5 +3519,4 @@
       REAL, DIMENSION (nloc)                             :: sument
       REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
-      REAL, DIMENSION (nloc)                             :: qnk
       REAL sumdq !jyg
 ! 
@@ -3437,4 +3629,33 @@
   END DO
 
+! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf"
+!-----------------------------------------------------------------
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!! Warning : this option leads to water conservation violation
+!!!           Expert only
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  DO il = 1, ncum
+    ma(il, nlp) = 0.
+    ma(il, 1)   = 0.
+  END DO
+  DO k = nl, 2, -1
+    DO il = 1, ncum
+      ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k)
+      cbmf(il) = max(cbmf(il), ma(il,k))
+    END DO
+  END DO
+  DO k = 2,nl
+    DO il = 1, ncum
+      IF (k <icb(il)) THEN
+        ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
+      ENDIF
+    END DO
+  END DO
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Line kept for compatibility with earlier versions
   DO k = 2, nl
     DO il = 1, ncum
@@ -3445,4 +3666,24 @@
   END DO
 
+  DO il = 1, ncum
+    ma(il, nlp) = 0.
+    ma(il, 1)   = 0.
+  END DO
+  DO k = nl, 2, -1
+    DO il = 1, ncum
+      ma(il, k) = ma(il, k+1) + m(il, k)
+    END DO
+  END DO
+  DO k = 2,nl
+    DO il = 1, ncum
+      IF (k <icb(il)) THEN
+        ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il)
+      ENDIF
+    END DO
+  END DO
+
+  ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
 !    print*,'cv3_yield avant ft'
 ! am is the part of cbmf taken from the first level
@@ -3581,14 +3822,72 @@
 !***    Compute convective mass fluxes upwd and dnwd      ***
 
+!
+! =================================================
+!              upward fluxes                      |
+! ------------------------------------------------
+!
 upwd(:,:) = 0.
 up_to(:,:) = 0.
 up_from(:,:) = 0.
-dnwd(:,:) = 0.
-dn_to(:,:) = 0.
-dn_from(:,:) = 0.
-!
-! =================================================
-!              upward fluxes                      |
-! ------------------------------------------------
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation 
+!! is taken into account. 
+!! WARNING : in the present version, taking into account the mass-flux decrease due to 
+!! precipitation ejection leads to water conservation violation.
+!
+! - Upward mass flux of mixed draughts
+!---------------------------------------
+DO i = 2, nl
+  DO j = 1, i-1
+    DO il = 1, ncum
+      IF (i<=inb(il)) THEN
+        up_to(il,i) = up_to(il,i) + ment(il,j,i)
+      ENDIF
+    ENDDO
+  ENDDO
+ENDDO
+!
+DO j = 3, nl
+  DO i = 2, j-1
+    DO il = 1, ncum
+      IF (j<=inb(il)) THEN
+        up_from(il,i) = up_from(il,i) + ment(il,i,j)
+      ENDIF
+    ENDDO
+  ENDDO
+ENDDO
+!
+! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer 
+!(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting 
+!from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)): 
+! 
+DO i = 2, nlp
+  DO il = 1, ncum
+    IF (i<=inb(il)+1) THEN
+      upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1))
+    ENDIF
+  ENDDO
+ENDDO
+!
+! - Total upward mass flux
+!---------------------------
+DO i = 2, nlp
+  DO il = 1, ncum
+    IF (i<=inb(il)+1) THEN
+      upwd(il,i) = upwd(il,i) + ma(il,i)
+    ENDIF
+  ENDDO
+ENDDO
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! The decrease of the adiabatic ascent mass flux due to ejection of precipitation 
+!! is not taken into account. 
+!
+! - Upward mass flux
+!-------------------
 DO i = 2, nl
   DO il = 1, ncum
@@ -3613,6 +3912,5 @@
   ENDDO
 ENDDO
-!!DO i = 2, nl
-!!  DO j = i+1, nl          !! Permuter les boucles i et j
+!
 DO j = 3, nl
   DO i = 2, j-1
@@ -3636,8 +3934,16 @@
   ENDDO
 ENDDO
+
+
+  ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 !
 ! =================================================
 !              downward fluxes                    |
 ! ------------------------------------------------
+dnwd(:,:) = 0.
+dn_to(:,:) = 0.
+dn_from(:,:) = 0.
 DO i = 1, nl
   DO j = i+1, nl
@@ -3650,6 +3956,4 @@
 ENDDO
 !
-!!DO i = 2, nl
-!!  DO j = 1, i-1          !! Permuter les boucles i et j 
 DO j = 1, nl
   DO i = j+1, nl
@@ -4423,4 +4727,35 @@
 !!!!
 !!!!      ENDDO
+
+!!  DO i = 1, nlp
+!!    DO il = 1, ncum
+!!      ma(il, i) = 0
+!!    END DO
+!!  END DO
+!!
+!!  DO i = 1, nl
+!!    DO j = i, nl
+!!      DO il = 1, ncum
+!!        ma(il, i) = ma(il, i) + m(il, j)
+!!      END DO
+!!    END DO
+!!  END DO
+
+!jyg<  (loops stop at nl)
+!!  DO i = nl + 1, nd
+!!    DO il = 1, ncum
+!!      ma(il, i) = 0.
+!!    END DO
+!!  END DO
+!>jyg
+
+!!  DO i = 1, nl
+!!    DO il = 1, ncum
+!!      IF (i<=(icb(il)-1)) THEN
+!!        ma(il, i) = 0
+!!      END IF
+!!    END DO
+!!  END DO
+
 !-----------------------------------------------------------
         ENDIF !(.NOT.ok_optim_yield)                      !|
@@ -4447,33 +4782,4 @@
 !>jyg
 
-  DO i = 1, nlp
-    DO il = 1, ncum
-      ma(il, i) = 0
-    END DO
-  END DO
-
-  DO i = 1, nl
-    DO j = i, nl
-      DO il = 1, ncum
-        ma(il, i) = ma(il, i) + m(il, j)
-      END DO
-    END DO
-  END DO
-
-!jyg<  (loops stop at nl)
-!!  DO i = nl + 1, nd
-!!    DO il = 1, ncum
-!!      ma(il, i) = 0.
-!!    END DO
-!!  END DO
-!>jyg
-
-  DO i = 1, nl
-    DO il = 1, ncum
-      IF (i<=(icb(il)-1)) THEN
-        ma(il, i) = 0
-      END IF
-    END DO
-  END DO
 
 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -4562,9 +4868,16 @@
 ! 14/01/15 AJ delta n'a rien à faire là...                                                 
     DO il = 1, ncum                                                  ! cld
-      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
+!!      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
+!!        siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) &          ! cld
+!!        *rrd*tvp(il, i)/p(il, i)/100.                                ! cld
+!!
+!!      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
+      sigaq = 0.
+      IF (wa(il,i)>0.0 .AND. iflag(il)<=1)  THEN                     ! cld
         siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) &          ! cld
-        *rrd*tvp(il, i)/p(il, i)/100.                                ! cld
-
-      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
+                     *rrd*tvp(il, i)/p(il, i)/100.                   ! cld
+        siga(il, i) = min(siga(il,i), 1.0)                           ! cld
+        sigaq = siga(il,i)*qta(il,i-1)                               ! cld
+      ENDIF
 
 ! IM cf. FH 
@@ -4578,9 +4891,10 @@
         sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1))    ! cld
         sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i))  ! cld
-        qtc(il, i) = (siga(il,i)*qnk(il)+sigment(il,i)*qtment(il,i)) & ! cld
+!!        qtc(il, i) = (siga(il,i)*qta(il,i-1)+sigment(il,i)*qtment(il,i)) & ! cld
+        qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld
                      /(siga(il,i)+sigment(il,i))                     ! cld
         sigt(il,i) = sigment(il, i) + siga(il, i)
 
-!        qtc(il, i) = siga(il,i)*qnk(il)+(1.-siga(il,i))*qtment(il,i) ! cld
+!        qtc(il, i) = siga(il,i)*qta(il,i-1)+(1.-siga(il,i))*qtment(il,i) ! cld
 !     print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)  
       		
@@ -4871,6 +5185,6 @@
       do k=1,nl
         do i=1,ncum
-	  hp(i,k)=h(i,k)
-	enddo
+          hp(i,k)=h(i,k)
+        enddo
       enddo
 
Index: /LMDZ6/trunk/libf/phylmd/cv3a_uncompress.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cv3a_uncompress.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cv3a_uncompress.F90	(revision 3496)
@@ -10,6 +10,6 @@
                            asupmaxmin, &
                            da, phi, mp, phi2, d1a, dam, sigij, &                ! RomP+AC+jyg
-                           clw, elij, evap, ep, epmlmMm, eplaMm, &              ! RomP
-                           wdtrainA, wdtrainM, &                                ! RomP
+                           qta, clw, elij, evap, ep, epmlmMm, eplaMm, &         ! RomP+jyg
+                           wdtrainA, wdtrainS, wdtrainM, &                      ! RomP
                            qtc, sigt,          &
                            epmax_diag, & ! epmax_cape
@@ -24,6 +24,6 @@
                            asupmaxmin1, &
                            da1, phi1, mp1, phi21, d1a1, dam1, sigij1, &         ! RomP+AC+jyg
-                           clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &        ! RomP
-                           wdtrainA1, wdtrainM1, &                              ! RomP
+                           qta1, clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &  ! RomP+jyg
+                           wdtrainA1, wdtrainS1, wdtrainM1, &                   ! RomP
                            qtc1, sigt1, &
                            epmax_diag1) ! epmax_cape
@@ -75,4 +75,5 @@
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: d1a, dam               !RomP
   REAL, DIMENSION (nloc, nd, nd), INTENT (IN)        :: sigij                  !RomP
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta                    !jyg
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: clw                    !RomP
   REAL, DIMENSION (nloc, nd, nd), INTENT (IN)        :: elij                   !RomP
@@ -81,5 +82,5 @@
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: eplamM                 !RomP+jyg
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qtc, sigt              !RomP
-  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: wdtrainA, wdtrainM     !RomP
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: wdtrainA, wdtrainS, wdtrainM     !RomP
 
   ! outputs:
@@ -111,4 +112,5 @@
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: d1a1, dam1 !RomP       !RomP
   REAL, DIMENSION (len, nd, nd), INTENT (OUT)        :: sigij1                 !RomP
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: qta1                   !jyg
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: clw1                   !RomP
   REAL, DIMENSION (len, nd, nd), INTENT (OUT)        :: elij1                  !RomP
@@ -117,5 +119,5 @@
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: eplamM1                !RomP+jyg
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: qtc1, sigt1            !RomP
-  REAL, DIMENSION (len, nd), INTENT (OUT)            :: wdtrainA1, wdtrainM1   !RomP
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: wdtrainA1, wdtrainS1, wdtrainM1   !RomP
 
 
@@ -175,4 +177,5 @@
         d1a1(idcum(i), k) = d1a(i, k) !RomP
         dam1(idcum(i), k) = dam(i, k) !RomP
+        qta1(idcum(i), k) = qta(i, k) !jyg
         clw1(idcum(i), k) = clw(i, k) !RomP
         evap1(idcum(i), k) = evap(i, k) !RomP
@@ -180,4 +183,5 @@
         eplamM1(idcum(i), k) = eplamM(i, k) !RomP+jyg
         wdtrainA1(idcum(i), k) = wdtrainA(i, k) !RomP
+        wdtrainS1(idcum(i), k) = wdtrainS(i, k) !RomP
         wdtrainM1(idcum(i), k) = wdtrainM(i, k) !RomP
         qtc1(idcum(i), k) = qtc(i, k)
@@ -300,4 +304,5 @@
       d1a1(:, 1:nl) = d1a(:, 1:nl)            !RomP
       dam1(:, 1:nl) = dam(:, 1:nl)            !RomP
+      qta1(:, 1:nl) = qta(:, 1:nl)            !jyg
       clw1(:, 1:nl) = clw(:, 1:nl)            !RomP
       evap1(:, 1:nl) = evap(:, 1:nl)          !RomP
@@ -305,4 +310,5 @@
       eplamM1(:, 1:nl) = eplamM(:, 1:nl)       !RomP+jyg
       wdtrainA1(:, 1:nl) = wdtrainA(:, 1:nl)  !RomP
+      wdtrainS1(:, 1:nl) = wdtrainS(:, 1:nl)  !RomP
       wdtrainM1(:, 1:nl) = wdtrainM(:, 1:nl)  !RomP
       qtc1(:, 1:nl) = qtc(:, 1:nl)
Index: /LMDZ6/trunk/libf/phylmd/cv3p2_closure.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cv3p2_closure.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cv3p2_closure.F90	(revision 3496)
@@ -23,4 +23,5 @@
   include "cvthermo.h"
   include "cv3param.h"
+  include "cvflag.h"
   include "YOMCST2.h"
   include "YOMCST.h"
Index: /LMDZ6/trunk/libf/phylmd/cv3p_mixing.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cv3p_mixing.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cv3p_mixing.F90	(revision 3496)
@@ -1,4 +1,4 @@
 SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
-                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qnk, &
+                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
                        unk, vnk, hp, tv, tvp, ep, clw, sig, &
                        Ment, Qent, hent, uent, vent, nent, &
@@ -29,5 +29,6 @@
   INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
-  REAL, DIMENSION (nloc), INTENT (IN)                :: qnk, unk, vnk
+  REAL, DIMENSION (nloc), INTENT (IN)                :: unk, vnk
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
   REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
@@ -173,5 +174,6 @@
                            .AND. (j<=inb(il))) THEN
 
-            rti = qnk(il) - ep(il, i)*clw(il, i)
+!!            rti = qnk(il) - ep(il, i)*clw(il, i)
+            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
             bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
 !jyg(from aj)<
@@ -219,4 +221,8 @@
             Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
             Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
+          ELSE IF (j > i) THEN
+            IF (prt_level >= 10) THEN
+              print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
+            ENDIF
           END IF ! new
         END DO
@@ -248,5 +254,6 @@
 !!!       Ment(il,i,i)=m(il,i)
         Ment(il, i, i) = 1.
-        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
+!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
+        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
         uent(il, i, i) = unk(il)
         vent(il, i, i) = vnk(il)
@@ -332,5 +339,6 @@
       IF (i>=icb(il) .AND. i<=inb(il)) THEN
         lwork(il) = (nent(il,i)/=0)
-        rti = qnk(il) - ep(il, i)*clw(il, i)
+!!        rti = qnk(il) - ep(il, i)*clw(il, i)
+        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
 !jyg<
         IF (cvflag_ice) THEN
@@ -462,5 +470,6 @@
             lwork(il)) THEN
           IF (Sij(il,i,j)>0.0) THEN
-            rti = qnk(il) - ep(il, i)*clw(il, i)
+!!            rti = qnk(il) - ep(il, i)*clw(il, i)
+            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
             Qmixmax(il) = Qmix(Sjmax(il))
             Qmixmin(il) = Qmix(Sjmin(il))
@@ -590,5 +599,6 @@
               lwork(il)) THEN
             IF (Sij(il,i,j)>0.0) THEN
-              rti = qnk(il) - ep(il, i)*clw(il, i)
+!!              rti = qnk(il) - ep(il, i)*clw(il, i)
+              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
 !!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
               Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
@@ -659,5 +669,6 @@
 ! cc        Ment(il,i,i)=m(il,i)
         Ment(il, i, i) = 1.
-        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
+!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
+        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
         uent(il, i, i) = unk(il)
         vent(il, i, i) = vnk(il)
Index: /LMDZ6/trunk/libf/phylmd/cv3param.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cv3param.h	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cv3param.h	(revision 3496)
@@ -7,13 +7,7 @@
 !------------------------------------------------------------
 
-      logical ok_homo_tend
-      logical ok_optim_yield
-      logical ok_entrain
-      logical ok_convstop
-      logical ok_intermittent
-      logical qsat_depends_on_qt
+      integer flag_epKEorig, flag_wb
+      integer cv_flag_feed
       integer noff, minorig, nl, nlp, nlm
-      integer cv_flag_feed
-      integer flag_epKEorig,flag_wb
       real sigdz, spfac
       real pbcrit, ptcrit
@@ -28,4 +22,6 @@
       real delta
       real betad
+      real ejectliq
+      real ejectice
 
       COMMON /cv3param/ sigdz, spfac &
@@ -40,12 +36,8 @@
                       ,wbmax &
                       ,delta, betad  &
+                      ,ejectliq, ejectice &
                       ,flag_epKEorig &
                       ,flag_wb, cv_flag_feed &
-                      ,noff, minorig, nl, nlp, nlm  &
-                      ,ok_convstop, ok_intermittent &
-                      ,ok_optim_yield &
-                      ,ok_entrain &
-                      ,ok_homo_tend &
-                      ,qsat_depends_on_qt
+                      ,noff, minorig, nl, nlp, nlm
 !$OMP THREADPRIVATE(/cv3param/)
 
Index: /LMDZ6/trunk/libf/phylmd/cva_driver.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cva_driver.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cva_driver.F90	(revision 3496)
@@ -25,6 +25,6 @@
 !!                      elij1,evap1,ep1,epmlmMm1,eplaMm1, &                ! RomP
                       da1, phi1, mp1, phi21, d1a1, dam1, sigij1, wghti1, & ! RomP, RL
-                      clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &        ! RomP, RL
-                      wdtrainA1, wdtrainM1, qtc1, sigt1, tau_cld_cv, &
+                      qta1, clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &  ! RomP, RL
+                      wdtrainA1, wdtrainS1, wdtrainM1, qtc1, sigt1, tau_cld_cv, &     !!jygprl
                       coefw_cld_cv, &                                      ! RomP, AJ
                       epmax_diag1)  ! epmax_cape
@@ -124,5 +124,7 @@
 !                                      of dimension ND, defined at same grid levels as T, Q, QS and P.
 
-! wdtrainA1     Real           Output   precipitation detrained from adiabatic draught;
+! wdtrainA1     Real           Output   precipitation ejected from adiabatic draught;
+!                                         should be used in tracer transport (cvltr)
+! wdtrainS1     Real           Output   precipitation detrained from shedding of adiabatic draught;
 !                                         used in tracer transport (cvltr)
 ! wdtrainM1     Real           Output   precipitation detrained from mixed draughts;
@@ -248,5 +250,5 @@
 
 ! RomP >>>
-  REAL, DIMENSION (len, nd), INTENT (OUT)            :: wdtrainA1, wdtrainM1 ! precipitation sources (extensive)
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: wdtrainA1, wdtrainS1, wdtrainM1 ! precipitation sources (extensive)
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: mp1  ! unsat. mass flux (staggered grid)
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: da1  ! detrained mass flux of adiab. asc. air (extensive)
@@ -258,4 +260,5 @@
   REAL, DIMENSION (len, nd, nd), INTENT (OUT)        :: sigij1 ! mass fraction of env. air in mixed draughts (intensive) 
   REAL, DIMENSION (len, nd, nd), INTENT (OUT)        :: elij1! cond. water per unit mass of mixed draughts (intensive) 
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: qta1 ! total water per unit mass of the adiab. asc. (intensive)
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: clw1 ! cond. water per unit mass of the adiab. asc. (intensive) 
 !JYG,RL
@@ -467,4 +470,5 @@
   REAL tv_wake(nloc, nd)
   REAL clw(nloc, nd)
+  REAL, DIMENSION(nloc, nd)    :: qta, qpreca                       !!jygprl
   REAL dph(nloc, nd)
   REAL pbase(nloc), buoybase(nloc), th(nloc, nd)
@@ -477,4 +481,5 @@
   REAL cin(nloc)
   REAL m(nloc, nd)
+  REAL mm(nloc, nd)
   REAL ment(nloc, nd, nd), sigij(nloc, nd, nd)
   REAL qent(nloc, nd, nd)
@@ -494,5 +499,5 @@
   REAL, DIMENSION(len,nd)     :: wt, water, evap
   REAL, DIMENSION(len,nd)     :: ice, fondue, b
-  REAL, DIMENSION(len,nd)     :: frac, faci
+  REAL, DIMENSION(len,nd)     :: frac_a, frac_s, faci               !!jygprl
   REAL ft(nloc, nd), fq(nloc, nd)
   REAL ftd(nloc, nd), fqd(nloc, nd)
@@ -523,5 +528,5 @@
  
 ! RomP >>>
-  REAL wdtrainA(nloc, nd), wdtrainM(nloc, nd)
+  REAL wdtrainA(nloc, nd), wdtrainS(nloc, nd), wdtrainM(nloc, nd)   !!jygprl
   REAL da(len, nd), phi(len, nd, nd)
   REAL epmlmMm(nloc, nd, nd), eplaMm(nloc, nd)
@@ -615,5 +620,5 @@
   tvp(:, :) = 0. !ym missing init, need to have a look by developpers
   tv(:, :) = 0. !ym missing init, need to have a look by developpers
-   
+
   DO il = 1, len
     cin1(il) = -100000.
@@ -636,4 +641,5 @@
   qtc1(:, :) = 0.
   wdtrainA1(:, :) = 0.
+  wdtrainS1(:, :) = 0.
   wdtrainM1(:, :) = 0.
   da1(:, :) = 0.
@@ -646,4 +652,5 @@
   sigij1(:, :, :) = 0.
   elij1(:, :, :) = 0.
+  qta1(:,:) = 0.
   clw1(:,:) = 0.
   wghti1(:,:) = 0.
@@ -906,5 +913,5 @@
                          p, ph, h, tv, lv, lf, pbase, buoybase, plcl, &
                          inb, tp, tvp, clw, hp, ep, sigp, buoy, &
-                         frac)
+                         frac_a, frac_s, qpreca, qta)                        !!jygprl
     END IF
 
@@ -915,5 +922,5 @@
                         tnk, qnk, gznk, t, q, qs, gz, &
                         p, dph, h, tv, lv, &
-                        inb, inbis, tp, tvp, clw, hp, ep, sigp, frac)
+                        inb, inbis, tp, tvp, clw, hp, ep, sigp, frac_s)
     END IF
 
@@ -923,5 +930,5 @@
              PRINT *, 'cva_driver -> cv3_epmax_cape'
     call cv3_epmax_fn_cape(nloc,ncum,nd &
-                , ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac &
+                , ep,hp,icb,inb,clw,nk,t,h,hnk,lv,lf,frac_s &
                 , pbase, p, ph, tv, buoy, sig, w0,iflag &
                 , epmax_diag)
@@ -941,5 +948,6 @@
              PRINT *, 'cva_driver -> cv3p_mixing'
         CALL cv3p_mixing(nloc, ncum, nd, nd, ntra, icb, nk, inb, &           ! na->nd
-                         ph, t, q, qs, u, v, tra, h, lv, lf, frac, qnk, &
+!!                         ph, t, q, qs, u, v, tra, h, lv, lf, frac, qnk, &
+                         ph, t, q, qs, u, v, tra, h, lv, lf, frac_s, qta, &      !!jygprl
                          unk, vnk, hp, tv, tvp, ep, clw, sig, &
                          ment, qent, hent, uent, vent, nent, &
@@ -1021,10 +1029,16 @@
              PRINT *, 'cva_driver -> cv3_mixing'
         CALL cv3_mixing(nloc, ncum, nd, nd, ntra, icb, nk, inb, &             ! na->nd
-                        ph, t, q, qs, u, v, tra, h, lv, lf, frac, qnk, &
+                        ph, t, q, qs, u, v, tra, h, lv, lf, frac_s, qnk, &
                         unk, vnk, hp, tv, tvp, ep, clw, m, sig, &
                         ment, qent, uent, vent, nent, sigij, elij, ments, qents, traent)
         CALL zilch(hent, nloc*nd*nd)
       ELSE
-        CALL cv3_mixscale(nloc, ncum, nd, ment, m)
+!!jyg:  Essais absurde pour voir
+!!        mm(:,1) = 0.
+!!        DO  i = 2,nd
+!!          mm(:,i) = m(:,i)*(1.-qta(:,i-1))
+!!        ENDDO
+        mm(:,:) = m(:,:)
+        CALL cv3_mixscale(nloc, ncum, nd, ment, mm)
         IF (debut) THEN
           PRINT *, ' cv3_mixscale-> '
@@ -1062,9 +1076,10 @@
                      t_wake, q_wake, qs_wake, gz, u, v, tra, p, ph, &
                      th_wake, tv_wake, lv_wake, lf_wake, cpn_wake, &
-                     ep, sigp, clw, &
+                     ep, sigp, clw, frac_s, qpreca, frac_a, qta, &                    !!jygprl
                      m, ment, elij, delt, plcl, coef_clos, &
                      mp, qp, up, vp, trap, wt, water, evap, fondue, ice, &
                      faci, b, sigd, &
-                     wdtrainA, wdtrainM)                                       ! RomP
+!!                     wdtrainA, wdtrainM)                                       ! RomP
+                     wdtrainA, wdtrainS, wdtrainM)                               !!jygprl
 !
       IF (prt_level >= 10) THEN
@@ -1075,8 +1090,8 @@
            evap(igout,k), fondue(igout,k)
         ENDDO
-        Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainM '
+        Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainS, wdtrainM '     !!jygprl
         DO k = 1,nd
-        write (6, '(i4,2(1x,e13.6))'), &
-           k, wdtrainA(igout,k), wdtrainM(igout,k)
+        write (6, '(i4,3(1x,e13.6))'), &
+           k, wdtrainA(igout,k), wdtrainS(igout,k), wdtrainM(igout,k)            !!jygprl
         ENDDO
       ENDIF
@@ -1112,5 +1127,5 @@
                      t, q, t_wake, q_wake, s_wake, u, v, tra, &
                      gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, &
-                     ep, clw, m, tp, mp, qp, up, vp, trap, &
+                     ep, clw, qpreca, m, tp, mp, qp, up, vp, trap, &
                      wt, water, ice, evap, fondue, faci, b, sigd, &
                      ment, qent, hent, iflag_mix, uent, vent, &
@@ -1121,5 +1136,8 @@
 !!                     tls, tps, &                            ! useless . jyg
                      qcondc, wd, &
-                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
+!!                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
+                     ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)         !!jygprl
+!
+!         Test conseravtion de l'eau
 !
       IF (debut) THEN
@@ -1142,5 +1160,5 @@
                      t, q, u, v, &
                      gz, p, ph, h, hp, lv, cpn, &
-                     ep, clw, frac, m, mp, qp, up, vp, &
+                     ep, clw, frac_s, m, mp, qp, up, vp, &
                      wt, water, evap, &
                      ment, qent, uent, vent, nent, elij, &
@@ -1187,6 +1205,6 @@
                            asupmaxmin, &
                            da, phi, mp, phi2, d1a, dam, sigij, &         ! RomP
-                           clw, elij, evap, ep, epmlmMm, eplaMm, &       ! RomP
-                           wdtrainA, wdtrainM, &                         ! RomP
+                           qta, clw, elij, evap, ep, epmlmMm, eplaMm, &  ! RomP
+                           wdtrainA, wdtrainS, wdtrainM, &                         ! RomP
                            qtc, sigt, epmax_diag, & ! epmax_cape
                            iflag1, kbas1, ktop1, &
@@ -1199,7 +1217,7 @@
                            Plim11, plim21, asupmax1, supmax01, &
                            asupmaxmin1, &
-                           da1, phi1, mp1, phi21, d1a1, dam1, sigij1,  & ! RomP
-                           clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, & ! RomP
-                           wdtrainA1, wdtrainM1,                       & ! RomP
+                           da1, phi1, mp1, phi21, d1a1, dam1, sigij1,  &       ! RomP
+                           qta1, clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, & ! RomP
+                           wdtrainA1, wdtrainS1, wdtrainM1,                  & ! RomP
                            qtc1, sigt1, epmax_diag1) ! epmax_cape
 !   
Index: /LMDZ6/trunk/libf/phylmd/cvflag.h
===================================================================
--- /LMDZ6/trunk/libf/phylmd/cvflag.h	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/cvflag.h	(revision 3496)
@@ -4,6 +4,22 @@
       logical cvflag_grav
       logical cvflag_ice
+      logical ok_optim_yield
+      logical ok_entrain
+      logical ok_homo_tend
+      logical ok_convstop
+      logical ok_intermittent
+      logical cvflag_prec_eject
+      logical qsat_depends_on_qt
+      logical adiab_ascent_mass_flux_depends_on_ejectliq
       integer icvflag_Tpa
 
-      COMMON /cvflag/ icvflag_Tpa, cvflag_grav, cvflag_ice
+      COMMON /cvflag/ icvflag_Tpa, &
+                      cvflag_grav, cvflag_ice, &
+                      ok_optim_yield, &
+                      ok_entrain, &
+                      ok_homo_tend, &
+                      ok_convstop, ok_intermittent, &
+                      cvflag_prec_eject, &
+                      qsat_depends_on_qt, &
+                      adiab_ascent_mass_flux_depends_on_ejectliq
 !$OMP THREADPRIVATE(/cvflag/)
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 3496)
@@ -416,6 +416,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: pmflxr, pmflxs
 !$OMP THREADPRIVATE(pmflxr, pmflxs)
-      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: wdtrainA, wdtrainM
-!$OMP THREADPRIVATE(wdtrainA, wdtrainM)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: wdtrainA, wdtrainS, wdtrainM
+!$OMP THREADPRIVATE(wdtrainA, wdtrainS, wdtrainM)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: da, mp
 !$OMP THREADPRIVATE(da, mp)
@@ -428,4 +428,6 @@
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: ev
 !$OMP THREADPRIVATE(ev)
+      REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: qtaa
+!$OMP THREADPRIVATE(qtaa)
       REAL,ALLOCATABLE,SAVE,DIMENSION(:,:)  :: clw
 !$OMP THREADPRIVATE(clw)
@@ -791,5 +793,5 @@
 !  Deep convective variables used in phytrac
       ALLOCATE(pmflxr(klon, klev+1), pmflxs(klon, klev+1))
-      ALLOCATE(wdtrainA(klon,klev),wdtrainM(klon,klev))
+      ALLOCATE(wdtrainA(klon,klev),wdtrainS(klon,klev),wdtrainM(klon,klev))
       ALLOCATE(dnwd(klon, klev), upwd(klon, klev) )
       ALLOCATE(ep(klon,klev))                          ! epmax_cape
@@ -801,4 +803,5 @@
       ALLOCATE(ev(klon,klev) )
       ALLOCATE(elij(klon,klev,klev) )
+      ALLOCATE(qtaa(klon,klev) )
       ALLOCATE(clw(klon,klev) )
       ALLOCATE(epmlmMm(klon,klev,klev), eplaMm(klon,klev) )
@@ -1082,5 +1085,5 @@
 
       DEALLOCATE(pmflxr, pmflxs)
-      DEALLOCATE(wdtrainA, wdtrainM)
+      DEALLOCATE(wdtrainA, wdtrainS, wdtrainM)
       DEALLOCATE(upwd, dnwd)
       DEALLOCATE(ep)
@@ -1092,4 +1095,5 @@
       DEALLOCATE(ev )
       DEALLOCATE(elij )
+      DEALLOCATE(qtaa )
       DEALLOCATE(clw )
       DEALLOCATE(epmlmMm, eplaMm )
Index: /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 3496)
@@ -1526,6 +1526,8 @@
   TYPE(ctrl_out), SAVE :: o_clwcon = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'clwcon', 'Convective Cloud Liquid water content', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_Mipsh = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'Mipsh', 'mass flux shed from adiab. ascents', 'kg/m2/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_Ma = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
-    'Ma', 'undilute adiab updraft', 'kg/m2/s', (/ ('', i=1, 10) /))
+    'Ma', 'undilute adiab updraft mass flux', 'kg/m2/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_dnwd = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'dnwd', 'saturated downdraft', 'kg/m2/s', (/ ('', i=1, 10) /))
@@ -1591,8 +1593,14 @@
   TYPE(ctrl_out), SAVE :: o_wdtrainA = ctrl_out((/ 4, 5, 10,  4, 10, 10, 11, 11, 11, 11 /), &
     'wdtrainA', 'precipitation from AA', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_wdtrainS = ctrl_out((/ 4, 5, 10,  4, 10, 10, 11, 11, 11, 11 /), &
+    'wdtrainS', 'precipitation from shedding of AA', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_wdtrainM = ctrl_out((/ 4, 5, 10,  4, 10, 10, 11, 11, 11, 11 /), &
     'wdtrainM', 'precipitation from mixture', '-', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_Vprecip = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'Vprecip', 'precipitation vertical profile', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_qtaa = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'qtaa', 'specific total water in adiabatic ascents', 'kg/kg', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_clwaa = ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
+    'Clwaa', 'specific condensed water in adiabatic ascents', 'kg/kg', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_ftd = ctrl_out((/ 4, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'ftd', 'tend temp due aux descentes precip', '-', (/ ('', i=1, 10) /))
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 3496)
@@ -73,5 +73,6 @@
          o_uwat, o_vwat, &
          o_ptop, o_fbase, o_plcl, o_plfc, &
-         o_wbeff, o_convoccur, o_cape_max, o_upwd, o_ep,o_epmax_diag, o_Ma, &
+         o_wbeff, o_convoccur, o_cape_max, o_upwd, o_ep,o_epmax_diag, &
+         o_Mipsh, o_Ma, &
          o_dnwd, o_dnwd0, o_ftime_deepcv, o_ftime_con, o_mc, &
          o_prw, o_prlw, o_prsw, o_s_pblh, o_s_pblt, o_s_lcl, &
@@ -88,5 +89,6 @@
          o_wake_s, o_wake_deltat, o_wake_deltaq, &
          o_wake_omg, o_dtwak, o_dqwak, o_dqwak2d, o_Vprecip, &
-         o_ftd, o_fqd, o_wdtrainA, o_wdtrainM, &
+         o_qtaa, o_Clwaa, &
+         o_ftd, o_fqd, o_wdtrainA, o_wdtrainS, o_wdtrainM, &
          o_n2, o_s2, o_proba_notrig, &
          o_random_notrig, o_ale_bl_stat, &
@@ -232,5 +234,5 @@
          delta_tsurf, &
          wstar, cape, ema_pcb, ema_pct, &
-         ema_cbmf, Ma, fm_therm, ale_bl, alp_bl, ale, &
+         ema_cbmf, Mipsh, Ma, fm_therm, ale_bl, alp_bl, ale, &
          alp, cin, wake_pe, wake_dens, wake_s, wake_deltat, &
          wake_deltaq, ftd, fqd, ale_bl_trig, albsol1, &
@@ -273,6 +275,6 @@
          kh         ,kh_x       ,kh_w       , &
          cv_gen, wake_h, &
-         wake_omg, d_t_wake, d_q_wake, Vprecip, &
-         wdtrainA, wdtrainM, n2, s2, proba_notrig, &
+         wake_omg, d_t_wake, d_q_wake, Vprecip, qtaa, Clw, &
+         wdtrainA, wdtrainS, wdtrainM, n2, s2, proba_notrig, &
          random_notrig, &
          alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
@@ -1058,4 +1060,5 @@
 
           CALL histwrite_phy(o_upwd, upwd)
+          CALL histwrite_phy(o_Mipsh, Mipsh)
           CALL histwrite_phy(o_Ma, Ma)
           CALL histwrite_phy(o_dnwd, dnwd)
@@ -1242,5 +1245,8 @@
           ! etendue a iflag_con=3 (jyg)
           CALL histwrite_phy(o_Vprecip, Vprecip)
+          CALL histwrite_phy(o_qtaa, qtaa)
+          CALL histwrite_phy(o_clwaa, clw)
           CALL histwrite_phy(o_wdtrainA, wdtrainA)
+          CALL histwrite_phy(o_wdtrainS, wdtrainS)
           CALL histwrite_phy(o_wdtrainM, wdtrainM)
        ENDIF !(iflag_con.EQ.3.or.iflag_con.EQ.30)
Index: /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 3496)
@@ -202,5 +202,7 @@
       REAL,ALLOCATABLE,SAVE :: ema_pcb(:), ema_pct(:)
 !$OMP THREADPRIVATE(ema_pcb,ema_pct)
-      REAL,ALLOCATABLE,SAVE :: Ma(:,:)        ! undilute upward mass flux
+      REAL,ALLOCATABLE,SAVE :: Mipsh(:,:)     ! mass flux shed from  adiab. ascents
+!$OMP THREADPRIVATE(Mipsh)
+      REAL,ALLOCATABLE,SAVE :: Ma(:,:)       ! undilute upward mass flux
 !$OMP THREADPRIVATE(Ma)
       REAL,ALLOCATABLE,SAVE :: qcondc(:,:)    ! in-cld water content from convect
@@ -541,4 +543,5 @@
       ALLOCATE(ema_pcb(klon), ema_pct(klon))
 !
+      ALLOCATE(Mipsh(klon,klev))
       ALLOCATE(Ma(klon,klev))
       ALLOCATE(qcondc(klon,klev))
@@ -699,5 +702,5 @@
       deallocate(ema_cbmf)
       deallocate(ema_pcb, ema_pct)
-      deallocate(Ma, qcondc)
+      deallocate(Mipsh, Ma, qcondc)
       deallocate(wd, sigd)
       deallocate(cin, ALE, ALP)
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3495)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 3496)
@@ -172,5 +172,5 @@
        !    Deep convective variables used in phytrac
        pmflxr, pmflxs,  &
-       wdtrainA, wdtrainM,  &
+       wdtrainA, wdtrainS, wdtrainM,  &
        upwd, dnwd, &
        ep,  &
@@ -182,4 +182,5 @@
        ev, &
        elij, &
+       qtaa, &
        clw, &
        epmlmMm, eplaMm, &
@@ -559,5 +560,6 @@
     ! Variables li\'ees \`a la poche froide (jyg)
 
-    REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
+!!    REAL mipsh(klon,klev)  ! mass flux shed by the adiab ascent at each level
+!!      Moved to phys_state_var_mod
     !
     REAL wape_prescr, fip_prescr
@@ -2629,4 +2631,5 @@
     pmflxs(:,:) = 0.
     wdtrainA(:,:) = 0.
+    wdtrainS(:,:) = 0.
     wdtrainM(:,:) = 0.
     upwd(:,:) = 0.
@@ -2644,4 +2647,5 @@
     elij(:,:,:)=0.
     ev(:,:)=0.
+    qtaa(:,:)=0.
     clw(:,:)=0.
     sij(:,:,:)=0.
@@ -2780,13 +2784,13 @@
                rain_con, snow_con, ibas_con, itop_con, sigd, &
                ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
-               Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
+               Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
                pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
                                 ! RomP >>>
                                 !!     .        pmflxr,pmflxs,da,phi,mp,
                                 !!     .        ftd,fqd,lalim_conv,wght_th)
-               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
+               pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
                ftd,fqd,lalim_conv,wght_th, &
                ev, ep,epmlmMm,eplaMm, &
-               wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
+               wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
                tau_cld_cv,coefw_cld_cv,epmax_diag)
 
