Index: LMDZ6/branches/Ocean_skin/libf/phylmd/cv3_routines.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/cv3_routines.F90	(revision 3428)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/cv3_routines.F90	(revision 3605)
@@ -35,4 +35,5 @@
 
   include "cv3param.h"
+  include "cvflag.h"
   include "conema3.h"
 
@@ -125,4 +126,16 @@
      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)
+     keepbug_ice_frac = .TRUE.
+     CALL getin_p('keepbug_ice_frac', keepbug_ice_frac)
 
     WRITE (*, *) 't_top_max=', t_top_max
@@ -144,4 +157,11 @@
     WRITE (*, *) 'elcrit=', elcrit
     WRITE (*, *) 'tlcrit=', tlcrit
+    WRITE (*, *) 'ejectliq=', ejectliq
+    WRITE (*, *) 'ejectice=', ejectice
+    WRITE (*, *) 'cvflag_prec_eject =', cvflag_prec_eject 
+    WRITE (*, *) 'qsat_depends_on_qt =', qsat_depends_on_qt 
+    WRITE (*, *) 'adiab_ascent_mass_flux_depends_on_ejectliq =', adiab_ascent_mass_flux_depends_on_ejectliq
+    WRITE (*, *) 'keepbug_ice_frac =', keepbug_ice_frac 
+
     first = .FALSE.
   END IF ! (first)
@@ -170,4 +190,5 @@
 
   include "cv3param.h"
+  include "cvflag.h"
 
 !inputs:
@@ -236,5 +257,6 @@
 ! 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)
+!!      lf(i, k) = lf0 - clmci*(t(i,k)-273.15)   ! erreur de signe !!
+      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)
@@ -289,4 +311,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
 
@@ -516,4 +539,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
 
 ! -------------------------------------------------------------------
@@ -1105,5 +1132,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
@@ -1153,22 +1181,45 @@
   REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ep, sigp, hp
   REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: buoy
-  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac
+  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:
   INTEGER i, j, k
-  REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
+  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, qi(nloc, nd)
-  REAL by, defrac, pden, tbis
-  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
-  LOGICAL lcape(nloc)
-  INTEGER iposit(nloc)
-  REAL fracg
-  REAL deltap
+  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
+  REAL                                               :: deltap
+  REAL, SAVE                                         :: Tx, Tm
+  DATA Tx/263.15/, Tm/243.15/
+!$OMP THREADPRIVATE(Tx, Tm)
+  REAL                                               :: aa, bb, dd, ddelta, discr
+  REAL                                               :: ff, fp
+  REAL                                               :: coefx, coefm, Zx, Zm, Ux, U, Um
 
   IF (prt_level >= 10) THEN
-    print *,'cv3_undilute2.0. t(1,k), q(1,k), qs(1,k) ', &
-                        (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
+    print *,'cv3_undilute2.0. icvflag_Tpa, t(1,k), q(1,k), qs(1,k) ', &
+                        icvflag_Tpa, (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
   ENDIF
+  smallestreal=tiny(smallestreal)
 
 ! =====================================================================
@@ -1181,4 +1232,5 @@
     END DO
   END DO
+
 
 ! =====================================================================
@@ -1197,9 +1249,298 @@
              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
+
+!jyg<
+! =====================================================================
+! --- SET THE THE FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
+! =====================================================================
+  DO k = 1, nl
+    DO i = 1, ncum
+      ep(i, k) = 0.0
+      sigp(i, k) = spfac
+    END DO
+  END DO
+!>jyg
+!
 
 ! ***  Find lifted parcel quantities above cloud base    ***
 
-
+!----------------------------------------------------------------------------
+!
+  IF (icvflag_Tpa == 2) THEN
+!
+!----------------------------------------------------------------------------
+!
+    DO k = minorig + 1, nl
+      DO i = 1,ncum
+        tp(i,k) = t(i,k)
+      ENDDO
+!!      alv = lv0 - clmcpv*(t(i,k)-273.15)
+!!      alf = lf0 + clmci*(t(i,k)-273.15)
+!!      als = alf + alv
+      DO j = 1,4
+        DO i = 1, ncum
+! ori	    if(k.ge.(icb(i)+1))then
+          IF (k>=(icbs(i)+1)) THEN                                ! convect3
+            tg = tp(i, k)
+            IF (tg .gt. Tx) THEN
+              es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
+              qg = eps*es/(p(i,k)-es*(1.-eps))
+            ELSE
+              esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
+              qg = eps*esi/(p(i,k)-esi*(1.-eps))
+            ENDIF
+! Ice fraction
+            ff = 0.
+            fp = 1./(Tx - Tm)
+            IF (tg < Tx) THEN
+              IF (tg > Tm) THEN
+                ff = (Tx - tg)*fp
+              ELSE
+                ff = 1.
+              ENDIF ! (tg > Tm)
+            ENDIF ! (tg < Tx)
+! Intermediate variables
+            aa = cpd + (cl-cpd)*qnk(i) + lv(i,k)*lv(i,k)*qg/(rrv*tg*tg)
+            ahg = (cpd + (cl-cpd)*qnk(i))*tg + lv(i,k)*qg - &
+                  lf(i,k)*ff*(qnk(i) - qg) + gz(i,k)
+            dd = lf(i,k)*lv(i,k)*qg/(rrv*tg*tg)
+            ddelta = lf(i,k)*(qnk(i) - qg)
+            bb = aa + ddelta*fp + dd*fp*(Tx-tg)
+! Compute Zx and Zm
+            coefx = aa
+            coefm = aa + dd
+            IF (tg .gt. Tx) THEN
+              Zx = ahg            + coefx*(Tx - tg)
+              Zm = ahg - ddelta   + coefm*(Tm - tg)
+            ELSE
+              IF (tg .gt. Tm) THEN
+                Zx = ahg          + (coefx +fp*ddelta)*(Tx - Tg)
+                Zm = ahg          + (coefm +fp*ddelta)*(Tm - Tg)
+              ELSE
+                Zx = ahg + ddelta + coefx*(Tx - tg)
+                Zm = ahg          + coefm*(Tm - tg)
+              ENDIF ! (tg .gt. Tm)
+            ENDIF ! (tg .gt. Tx)
+! Compute the masks Um, U, Ux
+            Um = (sign(1., Zm-ah0(i))+1.)/2.
+            Ux = (sign(1., ah0(i)-Zx)+1.)/2.
+            U = (1. - Um)*(1. - Ux)
+! Compute the updated parcell temperature Tp : 3 cases depending on tg value
+            IF (tg .gt. Tx) THEN
+              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tx-tg))
+              Tp(i,k) = tg + &
+                  Um*  (ah0(i) - ahg + ddelta)           /(aa + dd) + &
+                  U *2*(ah0(i) - ahg + ddelta*fp*(Tx-tg))/(bb + sqrt(discr)) + &
+                  Ux*  (ah0(i) - ahg)                    /aa
+            ELSEIF (tg .gt. Tm) THEN
+              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg)
+              Tp(i,k) = tg + &
+                  Um*  (ah0(i) - ahg + ddelta*fp*(tg-Tm))/(aa + dd) + &
+                  U *2*(ah0(i) - ahg)                    /(bb + sqrt(discr)) + &
+                  Ux*  (ah0(i) - ahg + ddelta*fp*(tg-Tx))/aa
+            ELSE
+              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tm-tg))
+              Tp(i,k) = tg + &
+                  Um*  (ah0(i) - ahg)                    /(aa + dd) + &
+                  U *2*(ah0(i) - ahg + ddelta*fp*(Tm-tg))/(bb + sqrt(discr)) + &
+                  Ux*  (ah0(i) - ahg - ddelta)           /aa
+            ENDIF ! (tg .gt. Tx)
+!
+!!     print *,' j, k, Um, U, Ux, aa, bb, discr, dd, ddelta ', j, k, Um, U, Ux, aa, bb, discr, dd, ddelta
+!!     print *,' j, k, ah0(i), ahg, tg, qg, tp(i,k), ff ', j, k, ah0(i), ahg, tg, qg, tp(i,k), ff
+          END IF ! (k>=(icbs(i)+1))
+        END DO ! i = 1, ncum
+      END DO ! j = 1,4
+      DO i = 1, ncum
+        IF (k>=(icbs(i)+1)) THEN                                ! convect3
+          tg = tp(i, k)
+          IF (tg .gt. Tx) THEN
+            es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
+            qg = eps*es/(p(i,k)-es*(1.-eps))
+          ELSE
+            esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
+            qg = eps*esi/(p(i,k)-esi*(1.-eps))
+          ENDIF
+          clw(i, k) = qnk(i) - qg
+          clw(i, k) = max(0.0, clw(i,k))
+          tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
+! print*,tvp(i,k),'tvp'
+          IF (clw(i,k)<1.E-11) THEN
+            tp(i, k) = tv(i, k)
+            tvp(i, k) = tv(i, k)
+          END IF ! (clw(i,k)<1.E-11)
+        END IF ! (k>=(icbs(i)+1))
+      END DO ! i = 1, ncum
+    END DO ! k = minorig + 1, nl
+!----------------------------------------------------------------------------
+!
+  ELSE IF (icvflag_Tpa == 1) THEN  ! (icvflag_Tpa == 2)
+!
+!----------------------------------------------------------------------------
+!
+    DO k = minorig + 1, nl
+      DO i = 1,ncum
+        tp(i,k) = t(i,k)
+      ENDDO
+!!      alv = lv0 - clmcpv*(t(i,k)-273.15)
+!!      alf = lf0 + clmci*(t(i,k)-273.15)
+!!      als = alf + alv
+      DO j = 1,4
+        DO i = 1, ncum
+! ori	    if(k.ge.(icb(i)+1))then
+          IF (k>=(icbs(i)+1)) THEN                                ! convect3
+            tg = tp(i, k)
+            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))
+              dqgdT = lv(i,k)*qg/(rrv*tg*tg)
+            ELSE
+              esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
+              qg = eps*esi/(p(i,k)-esi*(1.-eps))
+              dqgdT = (lv(i,k)+lf(i,k))*qg/(rrv*tg*tg)
+            ENDIF
+            IF (qsat_depends_on_qt) THEN
+              dqgdT = dqgdT*(1.-qta(i,k-1))/(1.-qg)**2
+              qg = qg*(1.-qta(i,k-1))/(1.-qg)            
+            ENDIF
+            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)*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
+          END IF ! (k>=(icbs(i)+1))
+        END DO ! i = 1, ncum
+      END DO ! j = 1,4
+      DO i = 1, ncum
+        IF (k>=(icbs(i)+1)) THEN                                ! convect3
+          tg = tp(i, k)
+          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))
+          ELSE
+            esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
+            qg = eps*esi/(p(i,k)-esi*(1.-eps))
+          ENDIF
+          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.+qhsat(i,k)/eps-qta(i,k-1)))
+! print*,tvp(i,k),'tvp'
+          IF (clw(i,k)<1.E-11) THEN
+            tp(i, k) = tv(i, k)
+            tvp(i, k) = tv(i, k)
+          END IF ! (clw(i,k)<1.E-11)
+        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 IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
+!
+!----------------------------------------------------------------------------
+!
   DO k = minorig + 1, nl
     DO i = 1, ncum
@@ -1358,23 +1699,14 @@
   END DO
 
-  IF (prt_level >= 10) THEN
-    print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
-                          (k, tp(1,k), tvp(1,k), k = 1,nl)
-  ENDIF
-
+!----------------------------------------------------------------------------
+!
+  ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0)
+!
+!----------------------------------------------------------------------------
+!
 ! =====================================================================
-! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
-! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
+! --- SET THE PRECIPITATION EFFICIENCIES 
 ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
 ! =====================================================================
-!
-!jyg<
-  DO k = 1, nl
-    DO i = 1, ncum
-      ep(i, k) = 0.0
-      sigp(i, k) = spfac
-    END DO
-  END DO
-!>jyg
 !
   IF (flag_epkeorig/=1) THEN
@@ -1413,4 +1745,10 @@
     END DO
   END IF
+!
+!   =========================================================================
+  IF (prt_level >= 10) THEN
+    print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
+                          (k, tp(1,k), tvp(1,k), k = 1,nl)
+  ENDIF
 !
 ! =====================================================================
@@ -1648,24 +1986,37 @@
   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
         IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
-          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)
+!jyg< frac computation moved to beginning of cv3_undilute2.
+!     kept here for compatibility test with CMip6 version
+          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
+  ENDIF  ! (cvflag_prec_eject) ELSE
+    DO k = minorig + 1, nl
       DO i = 1, ncum
-        IF (k<icb(i)) THEN
-          frac(i,k) = frac(i,icb(i))
+        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
+  ELSE   ! (cvflag_ice)
 !
     DO k = minorig + 1, nl
@@ -2350,9 +2701,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
@@ -2372,7 +2723,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
@@ -2387,5 +2743,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
@@ -2395,8 +2752,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
@@ -2405,8 +2763,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)
@@ -2415,4 +2775,6 @@
 ! ------------------------------------------------------
 IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
+
+smallestreal=tiny(smallestreal)
 
 ! =============================
@@ -2434,4 +2796,5 @@
 !! RomP >>>
 wdtrainA(:,:) = 0.
+wdtrainS(:,:) = 0.
 wdtrainM(:,:) = 0.
 !! RomP <<<
@@ -2489,4 +2852,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
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -2513,15 +2911,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
@@ -2531,11 +2926,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
@@ -2543,4 +2934,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   ***
@@ -2548,4 +2970,36 @@
 
 
+    IF (cvflag_ice) THEN                                                                                !!jygprl
+      IF (cvflag_prec_eject) THEN
+        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
+      ELSE  ! (cvflag_prec_eject)
+        DO il = 1, ncum                                                                                   !!jygprl
+          IF (i<=inb(il) .AND. lwork(il)) THEN                                                            !!jygprl
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+            IF (keepbug_ice_frac) THEN
+              frac(il, i) = frac_s(il, i)
+!       Ice fraction computed again here as a function of the temperature seen by unsaturated downdraughts
+!       (i.e. the cold pool temperature) for compatibility with earlier versions.
+              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)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+            ELSE  ! (keepbug_ice_frac)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+              frac(il, i) = frac_s(il, i)
+              fraci(il, i) = frac(il, i)                                                                    !!jygprl
+            ENDIF  ! (keepbug_ice_frac)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+          END IF                                                                                          !!jygprl
+        END DO                                                                                            !!jygprl
+      ENDIF  ! (cvflag_prec_eject)
+    END IF                                                                                              !!jygprl
+
+
     DO il = 1, ncum
       IF (i<=inb(il) .AND. lwork(il)) THEN
@@ -2553,29 +3007,9 @@
         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
-
-          IF (cvflag_ice) THEN
-!CR:tmax_fonte_cv: T for which ice is totally melted (used to be 275.15)
-            thaw = (t(il,i)-273.15)/(tmax_fonte_cv-273.15)
-            thaw = min(max(thaw,0.0), 1.0)
-            frac(il, i) = frac(il, i)*(1.-thaw)
-          ELSE
-            CONTINUE
-          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) = 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)
         rp(il, i) = max(rp(il,i), 0.0)
         rp(il, i) = amin1(rp(il,i), rs(il,i))
@@ -2998,4 +3432,5 @@
 
   RETURN
+
 END SUBROUTINE cv3_unsat
 
@@ -3004,5 +3439,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, &
@@ -3014,5 +3449,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
@@ -3054,5 +3489,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:
@@ -3083,4 +3520,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
@@ -3100,5 +3538,4 @@
       REAL, DIMENSION (nloc)                             :: sument
       REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
-      REAL, DIMENSION (nloc)                             :: qnk
       REAL sumdq !jyg
 ! 
@@ -3211,4 +3648,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
@@ -3219,4 +3685,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
@@ -3355,14 +3841,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
@@ -3387,6 +3931,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
@@ -3410,8 +3953,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
@@ -3424,6 +3975,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
@@ -3749,5 +4298,10 @@
     END DO ! cld
 
+!ym BIG Warning : it seems that the k loop is missing !!!
+!ym Strong advice to check this
+!ym add a k loop temporary 
+
 ! (particular case: no detraining level is found)                              ! cld
+! Verif merge Dynamico<<<<<<< .working
     DO il = 1, ncum                                                            ! cld
       IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
@@ -3761,4 +4315,15 @@
       END IF                                                                   ! cld
     END DO                                                                     ! cld
+! Verif merge Dynamico =======
+! Verif merge Dynamico     DO k = i + 1, nl
+! Verif merge Dynamico       DO il = 1, ncum        !ym k loop added                                    ! cld
+! Verif merge Dynamico         IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
+! Verif merge Dynamico           qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
+! Verif merge Dynamico           qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
+! Verif merge Dynamico           nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
+! Verif merge Dynamico         END IF                                                                   ! cld
+! Verif merge Dynamico       END DO
+! Verif merge Dynamico     ENDDO                                                                     ! cld
+! Verif merge Dynamico >>>>>>> .merge-right.r3413
 
     DO il = 1, ncum                                                            ! cld
@@ -4181,4 +4746,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)                      !|
@@ -4205,33 +4801,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
@@ -4320,9 +4887,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 
@@ -4336,9 +4910,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)  
       		
@@ -4629,6 +5204,6 @@
       do k=1,nl
         do i=1,ncum
-	  hp(i,k)=h(i,k)
-	enddo
+          hp(i,k)=h(i,k)
+        enddo
       enddo
 
