SUBROUTINE cv3_unsat24(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, frac_s, qpreca, frac_a, qta , &                       !!jygprl
                     m, ment, elij, delt, plcl, coef_clos, &
                     mpdn, hpdn, rpdn, updn, vpdn, trap, wt, water, evap, fondue, ice, &
                     faci, b, sigd, &
                     wdtrainA, wdtrainS, wdtrainM)                                      ! RomP
  USE print_control_mod, ONLY: prt_level, lunout
  IMPLICIT NONE


  include "cvthermo.h"
  include "cv3param.h"
  include "cvflag.h"
  include "nuage.h"

!inputs:
  INTEGER, INTENT (IN)                               :: ncum, nd, na, ntra, nloc
  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb
  REAL, INTENT(IN)                                   :: delt
  REAL, DIMENSION (nloc), INTENT (IN)                :: plcl
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: gz
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
  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
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: m
  REAL, DIMENSION (nloc, na, na), INTENT (IN)        :: ment, elij
  REAL, DIMENSION (nloc), INTENT (IN)                :: coef_clos

!input/output
  INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag(nloc)

!outputs:
  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: mpdn, hpdn, rpdn, updn, vpdn
  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: water, evap, wt
  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
  REAL, DIMENSION (nloc), INTENT (OUT)               :: sigd
! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
! de l ascendance adiabatique et des flux melanges Pa et Pm.
! Distinction des wdtrain
! Pa = wdtrainA     Pm = wdtrainM
  REAL, 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
  REAL pr1, pr2, sigt, b6, c6, d6, e6, f6, revap, delth
  REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf
  REAL ampmax, thaw
  REAL tevap(nloc)
  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)
!
!!jyg!    Pour nouveau calcul des flux de masse
!!jyg  REAL, DIMENSION (nloc, na)      :: mp_diag, rp_diag, hp_diag
!!jyg  REAL, DIMENSION (nloc, na)      :: hp


! ------------------------------------------------------
IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1)
                                                                  
print *,'XXXunsat24 nl, inb(1) ', nl, inb(1)

smallestreal=tiny(smallestreal)

! =============================
! --- INITIALIZE OUTPUT ARRAYS 
! =============================
!  (loops up to nl+1)
mp_diag(:,:) = 0.
!!jyg rp_diag(:,:) = 0.
!!jyg hp_diag(:,:) = 0.
hpdn(:,:) = 0.
mpdn(:,:) = 0.
rpdn(:,:) = 0.
updn(:,:) = 0.
vpdn(:,:) = 0.
water(:,:) = 0.
evap(:,:) = 0.
wt(:,:) = 0.
ice(:,:) = 0.
fondue(:,:) = 0.
faci(:,:) = 0.
b(:,:) = 0.
sigd(:) = 0.
!! RomP >>>
wdtrainA(:,:) = 0.
wdtrainS(:,:) = 0.
wdtrainM(:,:) = 0.
!! RomP <<<

  DO i = 1, nlp
    DO il = 1, ncum
      rpdn(il, i) = rr(il, i)
      updn(il, i) = u(il, i)
      vpdn(il, i) = v(il, i)
      wt(il, i) = 0.001
    END DO
  END DO

! ***  Set the fractionnal area sigd of precipitating downdraughts
  DO il = 1, ncum
    sigd(il) = sigdz*coef_clos(il)
  END DO

! =====================================================================
! --- INITIALIZE VARIOUS ARRAYS AND PARAMETERS USED IN THE COMPUTATIONS
! =====================================================================
!  (loops up to nl+1)

  delti = 1./delt
  tinv = 1./3.

  DO i = 1, nlp
    DO il = 1, ncum
      frac(il, i) = 0.0
      fraci(il, i) = 0.0
      prec(il, i) = 0.0
      lvcp(il, i) = lv(il, i)/cpn(il, i)
      lfcp(il, i) = lf(il, i)/cpn(il, i)
    END DO
  END DO

!AC!        do k=1,ntra
!AC!         do i=1,nd
!AC!          do il=1,ncum
!AC!           trap(il,i,k)=tra(il,i,k)
!AC!          enddo
!AC!         enddo
!AC!        enddo

! ***  check whether ep(inb)=0, if so, skip precipitating    ***
! ***             downdraft calculation                      ***


  DO il = 1, ncum
!!          lwork(il)=.TRUE.
!!          if(ep(il,inb(il)).lt.0.0001)lwork(il)=.FALSE.
!jyg<
!!    lwork(il) = ep(il, inb(il)) >= 0.0001
    lwork(il) = ep(il, inb(il)) >= 0.0001 .AND. iflag(il) <= 2
  END DO

  print *,'XXXunsat24 lwork(1) ', lwork(1)

!
! 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! ***                    begin downdraft loop                    ***
!
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  DO i = nl + 1, 1, -1

    num1 = 0
    DO il = 1, ncum
      IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1
    END DO

  print *,'XXXunsat24 i, num1 ', i, num1

    IF (num1<=0) GO TO 400

    CALL zilch(wdtrain, ncum)


! ***  integrate liquid water equation to find condensed water   ***
! ***                and condensed water flux                    ***
!
!
! ***              calculate detrained precipitation             ***


    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
      DO j = 1, i - 1
        DO il = 1, ncum
          IF (i<=inb(il) .AND. lwork(il)) THEN
            awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i)
            awat = max(awat, 0.0)
            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
      END DO
    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   ***
! ***              estimates of rp(i)and rp(i-1)             ***


    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

        wt(il, i) = 45.0

        IF (i<inb(il)) THEN
          rpdn(il, i) = rpdn(il, i+1) + &
                      (cpd*(t(il,i+1)-t(il,i))+gz(il,i+1)-gz(il,i))/lv(il, i)
          rpdn(il, i) = 0.5*(rpdn(il,i)+rr(il,i))
        END IF
        rpdn(il, i) = max(rpdn(il,i), 0.0)
        rpdn(il, i) = amin1(rpdn(il,i), rs(il,i))
        rpdn(il, inb(il)) = rr(il, inb(il))

        IF (i==1) THEN
          afac = p(il, 1)*(rs(il,1)-rpdn(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
          IF (cvflag_ice) THEN
            afac1 = p(il, i)*(rs(il,1)-rpdn(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1))
          END IF
        ELSE
          rpdn(il, i-1) = rpdn(il, i) + (cpd*(t(il,i)-t(il,i-1))+gz(il,i)-gz(il,i-1))/lv(il, i)
          rpdn(il, i-1) = 0.5*(rpdn(il,i-1)+rr(il,i-1))
          rpdn(il, i-1) = amin1(rpdn(il,i-1), rs(il,i-1))
          rpdn(il, i-1) = max(rpdn(il,i-1), 0.0)
          afac1 = p(il, i)*(rs(il,i)-rpdn(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i))
          afac2 = p(il, i-1)*(rs(il,i-1)-rpdn(il,i-1))/(1.0E4+2000.0*p(il,i-1)*rs(il,i-1))
          afac = 0.5*(afac1+afac2)
        END IF
        IF (i==inb(il)) afac = 0.0
        afac = max(afac, 0.0)
        bfac = 1./(sigd(il)*wt(il,i))

!
    IF (prt_level >= 20) THEN
      Print*, 'cv3_unsat after provisional rpdn estimate: rpdn, afac, bfac ', &
          i, rpdn(1, i), afac,bfac
    ENDIF
!
!JYG1
! cc        sigt=1.0
! cc        if(i.ge.icb)sigt=sigp(i)
! prise en compte de la variation progressive de sigt dans
! les couches icb et icb-1:
! pour plcl<ph(i+1), pr1=0 & pr2=1
! pour plcl>ph(i),   pr1=1 & pr2=0
! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval
! sur le nuage, et pr2 est la proportion sous la base du
! nuage.
        pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1))
        pr1 = max(0., min(1.,pr1))
        pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1))
        pr2 = max(0., min(1.,pr2))
        sigt = sigp(il, i)*pr1 + pr2
!JYG2

!JYG----
!    b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
!    c6 = water(il,i+1) + wdtrain(il)*bfac
!    c6 = prec(il,i+1) + wdtrain(il)*bfac
!    revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
!    evap(il,i)=sigt*afac*revap
!    water(il,i)=revap*revap
!    prec(il,i)=revap*revap
!!        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ', &
!!                 i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
!!---end jyg---

! --------retour à la formulation originale d'Emanuel.
        IF (cvflag_ice) THEN

!   b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
!   c6=prec(il,i+1)+bfac*wdtrain(il) &
!       -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
!   if(c6.gt.0.0)then
!   revap=0.5*(-b6+sqrt(b6*b6+4.*c6))

!JAM  Attention: evap=sigt*E
!    Modification: evap devient l'évaporation en milieu de couche
!    car nécessaire dans cv3_yield
!    Du coup, il faut modifier pas mal d'équations...
!    et l'expression de afac qui devient afac1
!    revap=sqrt((prec(i+1)+prec(i))/2)

          b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac1
          c6 = prec(il, i+1) + 0.5*bfac*wdtrain(il)
! print *,'bfac,sigd(il),sigt,afac1 ',bfac,sigd(il),sigt,afac1
! print *,'prec(il,i+1),wdtrain(il) ',prec(il,i+1),wdtrain(il)
! print *,'b6,c6,b6*b6+4.*c6 ',b6,c6,b6*b6+4.*c6
          IF (c6>b6*b6+1.E-20) THEN
            revap = 2.*c6/(b6+sqrt(b6*b6+4.*c6))
          ELSE
            revap = (-b6+sqrt(b6*b6+4.*c6))/2.
          END IF
          prec(il, i) = max(0., 2.*revap*revap-prec(il,i+1))
! print*,prec(il,i),'neige'

!JYG    Dans sa formulation originale, Emanuel calcule l'evaporation par:
! c             evap(il,i)=sigt*afac*revap
! ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
! Ici,l'evaporation evap est simplement calculee par l'equation de
! conservation.
! prec(il,i)=revap*revap
! else
!JYG----   Correction : si c6 <= 0, water(il,i)=0.
! prec(il,i)=0.
! endif

!JYG---   Dans tous les cas, evaporation = [tt ce qui entre dans la couche i]
! moins [tt ce qui sort de la couche i]
! print *, 'evap avec ice'
          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(prec(il,i+1)-prec(il,i))) / &
                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)
!
    IF (prt_level >= 20) THEN
      Print*, 'cv3_unsat after evap computation: wdtrain, sigd, wt, prec(i+1),prec(i) ', &
          i, wdtrain(1), sigd(1), wt(1,i), prec(1,i+1),prec(1,i)
    ENDIF
!

!jyg<
          d6 = prec(il,i)-prec(il,i+1)

!!          d6 = bfac*wdtrain(il) - 100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
!!          e6 = bfac*wdtrain(il)
!!          f6 = -100.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i)
!>jyg
!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)
!jyg<
          water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
          ice(il, i)   = ice(il, i+1)   + fraci(il, i)*d6
          water(il, i) = min(prec(il,i), max(water(il,i), 0.))
          ice(il, i)   = min(prec(il,i), max(ice(il,i),   0.))

!!          water(il, i) = water(il, i+1) + (1-fraci(il,i))*d6
!!          water(il, i) = max(water(il,i), 0.)
!!          ice(il, i) = ice(il, i+1) + fraci(il, i)*d6
!!          ice(il, i) = max(ice(il,i), 0.)
!>jyg
          fondue(il, i) = ice(il, i)*thaw
          water(il, i) = water(il, i) + fondue(il, i)
          ice(il, i) = ice(il, i) - fondue(il, i)

          IF (water(il,i)+ice(il,i)<1.E-30) THEN
            faci(il, i) = 0.
          ELSE
            faci(il, i) = ice(il, i)/(water(il,i)+ice(il,i))
          END IF

!           water(il,i)=water(il,i+1)+(1.-fraci(il,i))*e6+(1.-faci(il,i))*f6
!           water(il,i)=max(water(il,i),0.)
!           ice(il,i)=ice(il,i+1)+fraci(il,i)*e6+faci(il,i)*f6
!           ice(il,i)=max(ice(il,i),0.)
!           fondue(il,i)=ice(il,i)*thaw
!           water(il,i)=water(il,i)+fondue(il,i)
!           ice(il,i)=ice(il,i)-fondue(il,i)
            
!           if((water(il,i)+ice(il,i)).lt.1.e-30)then
!             faci(il,i)=0.
!           else
!             faci(il,i)=ice(il,i)/(water(il,i)+ice(il,i))
!           endif

        ELSE
          b6 = bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
          c6 = water(il, i+1) + bfac*wdtrain(il) - &
               50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il, i+1)
          IF (c6>0.0) THEN
            revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
            water(il, i) = revap*revap
          ELSE
            water(il, i) = 0.
          END IF
! print *, 'evap sans ice'
          evap(il, i) = (wdtrain(il)+sigd(il)*wt(il,i)*(water(il,i+1)-water(il,i)))/ &
                        (sigd(il)*(ph(il,i)-ph(il,i+1))*100.)

        END IF
      END IF !(i.le.inb(il) .and. lwork(il))
    END DO
! ----------------------------------------------------------------
!  ------------------ IM + JYG + NR -------------------------------
! |    Precipitating downdraft mass flux computation               |
!  ----------------------------------------------------------------
!
    IF (i /= 1 .AND. i /= nl+1) THEN
  print *,'XXXunsat24 avant mf i, (mpdn(1,j), j=32,41) ', i, (mpdn(1,j), j=32,41)
      CALL cv3_unsat_mf(nloc, ncum, inb, lwork, i,     &
                        t(:,i), rr(:,i),               &
                        cpn(:,i),                      & 
                        tv(:,i), tv(:,i-1), tv(:,i+1), &
                        gz(:,i),                  &
                        p(:,i), p(:,i-1), p(:,i+1), ph(:,i), ph(:,i+1),  &
                        lv(:,i), lf(:,i),   &
                        sigd,  &
                        faci(:,i),  &
                        evap(:,i), fondue(:,i), wt(:,i),  &
!!                        mp_diag(:,i+1),   &
                        mpdn(:,i+1),   &
!!                        rp_diag(:,i+1), hp_diag(:,i+1),   &
                        rpdn(:,i+1), hpdn(:,i+1),   &
!IM121224               mp_diag(:,i),   &
!IM121224               rp_diag(:,i), hp_diag(:,i)  &
                        mpdn(:,i),   &
                        rpdn(:,i), hpdn(:,i)  &
                        )
  print *,'XXXunsat24 apres mf i, (mpdn(1,j), j=32,41) ', i, (mpdn(1,j), j=32,41)
  print *, ' '
    ENDIF ! (i /= 1)

! cc
! ***  calculate precipitating downdraft mass flux under     ***
! ***              hydrostatic approximation                 ***

    DO il = 1, ncum
      IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN

        tevap(il) = max(0.0, evap(il,i))
        delth = max(0.001, (th(il,i)-th(il,i-1)))
!IM       IF (cvflag_ice) THEN
!IM         IF (cvflag_grav) THEN
!IM121224           mp(il, i) = 100.*ginv*(lvcp(il,i)*sigd(il)*tevap(il)* &
!IM121224                                              (p(il,i-1)-p(il,i))/delth + &
!IM121224                                  lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
!IM121224                                              (p(il,i-1)-p(il,i))/delth + &
!IM121224                                  lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
!IM121224                                              (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))
!IM         ELSE
!IM121224           mp(il, i) = 10.*(lvcp(il,i)*sigd(il)*tevap(il)* &
!IM121224                                               (p(il,i-1)-p(il,i))/delth + &
!IM121224                            lfcp(il,i)*sigd(il)*faci(il,i)*tevap(il)* &
!IM121224                                               (p(il,i-1)-p(il,i))/delth + &
!IM121224                            lfcp(il,i)*sigd(il)*wt(il,i)/100.*fondue(il,i)* &
!IM121224                                               (p(il,i-1)-p(il,i))/delth/(ph(il,i)-ph(il,i+1)))

!IM         END IF
!IM       ELSE
!IM         IF (cvflag_grav) THEN
!IM121224           mp(il, i) = 100.*ginv*lvcp(il, i)*sigd(il)*tevap(il)* &
!IM121224                                               (p(il,i-1)-p(il,i))/delth
!IM         ELSE
!IM121224           mp(il, i) = 10.*lvcp(il, i)*sigd(il)*tevap(il)* &
!IM121224                                               (p(il,i-1)-p(il,i))/delth
!IM         END IF

!IM       END IF

      END IF !(i.le.inb(il) .and. lwork(il) .and. i.ne.1)
      IF (prt_level .GE. 20) THEN
        PRINT *,'cv3_unsat, mpdn hydrostatic ', i, mpdn(il,i)
      ENDIF
    END DO
! ----------------------------------------------------------------

! ***           if hydrostatic assumption fails,             ***
! ***   solve cubic difference equation for downdraft theta  ***
! ***  and mass flux from two simultaneous differential eqns ***

    DO il = 1, ncum
      IF (i<=inb(il) .AND. lwork(il) .AND. i/=1) THEN

!IMJYG110225        amfac = sigd(il)*sigd(il)*70.0*ph(il, i)*(p(il,i-1)-p(il,i))* &
!IMJYG110225                         (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
!IMJYG110225        amp2 = abs(mpdn(il,i+1)*mpdn(il,i+1)-mpdn(il,i)*mpdn(il,i))
!IMJYG110225
!IMJYG110225        IF (amp2>(0.1*amfac)) THEN
!IMJYG110225          xf = 100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1))
!IMJYG110225          tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i) / &
!IMJYG110225                              (lvcp(il,i)*sigd(il)*th(il,i))
!IMJYG110225          af = xf*tf + mpdn(il, i+1)*mpdn(il, i+1)*tinv
!IMJYG110225
!IMJYG110225          IF (cvflag_ice) THEN
!IMJYG110225            bf = 2.*(tinv*mpdn(il,i+1))**3 + tinv*mpdn(il, i+1)*xf*tf + &
!IMJYG110225                 50.*(p(il,i-1)-p(il,i))*xf*(tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
!IMJYG110225                (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i)/(ph(il,i)-ph(il,i+1)))
!IMJYG110225          ELSE
!IMJYG110225
!IMJYG110225            bf = 2.*(tinv*mpdn(il,i+1))**3 + tinv*mpdn(il, i+1)*xf*tf + &
!IMJYG110225                                           50.*(p(il,i-1)-p(il,i))*xf*tevap(il)
!IMJYG110225          END IF
!IMJYG110225
!IMJYG110225          fac2 = 1.0
!IMJYG110225          IF (bf<0.0) fac2 = -1.0
!IMJYG110225          bf = abs(bf)
!IMJYG110225          ur = 0.25*bf*bf - af*af*af*tinv*tinv*tinv
!IMJYG110225          IF (ur>=0.0) THEN
!IMJYG110225            sru = sqrt(ur)
!IMJYG110225            fac = 1.0
!IMJYG110225            IF ((0.5*bf-sru)<0.0) fac = -1.0
!IMJYG110225            mpdn(il, i) = mpdn(il, i+1)*tinv + (0.5*bf+sru)**tinv + &
!IMJYG110225                                           fac*(abs(0.5*bf-sru))**tinv
!IMJYG110225          ELSE
!IMJYG110225            d = atan(2.*sqrt(-ur)/(bf+1.0E-28))
!IMJYG110225            IF (fac2<0.0) d = 3.14159 - d
!IMJYG110225            mpdn(il, i) = mpdn(il, i+1)*tinv + 2.*sqrt(af*tinv)*cos(d*tinv)
!IMJYG110225          END IF
!IMJYG110225          mpdn(il, i) = max(0.0, mpdn(il,i))
!IMJYG110225          IF (prt_level .GE. 20) THEN
!IMJYG110225            PRINT *,'cv3_unsat, mp cubic ', i, mpdn(il,i)
!IMJYG110225          ENDIF
!IMJYG110225
!IMJYG110225          IF (cvflag_ice) THEN
!IMJYG110225            IF (cvflag_grav) THEN
!IMJYG110225!JYG : il y a vraisemblablement une erreur dans la ligne 2 suivante:
!IMJYG110225! il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
!IMJYG110225! Et il faut bien revoir les facteurs 100.
!IMJYG110225              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))* &
!IMJYG110225                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
!IMJYG110225                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
!IMJYG110225                           (ph(il,i)-ph(il,i+1))) / &
!IMJYG110225                           (mpdn(il,i)+sigd(il)*0.1) - &
!IMJYG110225                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
!IMJYG110225                           (lvcp(il,i)*sigd(il)*th(il,i))
!IMJYG110225            ELSE
!IMJYG110225              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*&
!IMJYG110225                           (tevap(il)*(1.+(lf(il,i)/lv(il,i))*faci(il,i)) + &
!IMJYG110225                           (lf(il,i)/lv(il,i))*wt(il,i)/100.*fondue(il,i) / &
!IMJYG110225                           (ph(il,i)-ph(il,i+1))) / &
!IMJYG110225                           (mpdn(il,i)+sigd(il)*0.1) - &
!IMJYG110225                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
!IMJYG110225                           (lvcp(il,i)*sigd(il)*th(il,i))
!IMJYG110225            END IF
!IMJYG110225          ELSE
!IMJYG110225            IF (cvflag_grav) THEN
!IMJYG110225              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
!IMJYG110225                           (mpdn(il,i)+sigd(il)*0.1) - &
!IMJYG110225                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
!IMJYG110225                           (lvcp(il,i)*sigd(il)*th(il,i))
!IMJYG110225            ELSE
!IMJYG110225              b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap(il) / &
!IMJYG110225                           (mpdn(il,i)+sigd(il)*0.1) - &
!IMJYG110225                           10.0*(th(il,i)-th(il,i-1))*t(il, i) / &
!IMJYG110225                           (lvcp(il,i)*sigd(il)*th(il,i))
!IMJYG110225            END IF
!IMJYG110225          END IF
!IMJYG110225          b(il, i-1) = max(b(il,i-1), 0.0)
!IMJYG110225
!IMJYG110225        END IF !(amp2.gt.(0.1*amfac))
!IMJYG110225
!jyg<    This part shifted 10 lines farther
!!! ***         limit magnitude of mp(i) to meet cfl condition      ***
!!
!!        ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
!!        amp2 = 2.0*(ph(il,i-1)-ph(il,i))*delti
!!        ampmax = min(ampmax, amp2)
!!        mp(il, i) = min(mp(il,i), ampmax)
!>jyg

! ***      force mp to decrease linearly to zero                 ***
! ***       between cloud base and the surface                   ***


! c      if(p(il,i).gt.p(il,icb(il)))then
! c       mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
! c      endif
!IM121224 BEGIN ???
        IF (ph(il,i)>0.9*plcl(il)) THEN
          mpdn(il, i) = mpdn(il, i)*(ph(il,1)-ph(il,i))/(ph(il,1)-0.9*plcl(il))
        END IF
!IM121224 END ???

!jyg<    Shifted part
! ***         limit magnitude of mp(i) to meet cfl condition      ***

!IM121224 BEGIN ???
        ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti
        amp2 = 2.0*(ph(il,i-1)-ph(il,i))*delti
        ampmax = min(ampmax, amp2)
        mpdn(il, i) = min(mpdn(il,i), ampmax)
!IM121224 END ???
!>jyg

      END IF ! (i.le.inb(il) .and. lwork(il) .and. i.ne.1)
    END DO
! ----------------------------------------------------------------
!
    IF (prt_level >= 20) THEN
      Print*, 'cv3_unsat after mpdn computation: mpdn, b(i), b(i-1) ', &
          i, mpdn(1, i), b(1,i), b(1,max(i-1,1))
    ENDIF
!

! ***       find mixing ratio of precipitating downdraft     ***

    DO il = 1, ncum
      IF (i<inb(il) .AND. lwork(il)) THEN
        mplus(il) = mpdn(il, i) > mpdn(il, i+1)
      END IF ! (i.lt.inb(il) .and. lwork(il))
    END DO

    DO il = 1, ncum
      IF (i<inb(il) .AND. lwork(il)) THEN

!IM121224  ???
!IM       rp(il, i) = rr(il, i)
!IM       hp(il, i) = cpn(il, i)*t(il,i)+gz(il,i)
!IM121224  ???

        IF (mplus(il)) THEN

!IM121224         IF (cvflag_grav) THEN
!IM121224           rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
!IM121224             100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
!IM121224         ELSE
!IM121224           rp(il, i) = rp(il, i+1)*mp(il, i+1) + rr(il, i)*(mp(il,i)-mp(il,i+1)) + &
!IM121224             5.*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
!IM121224         END IF
!IM121224         rp(il, i) = rp(il, i)/mp(il, i)
!IM121224         hp(il, i) = hp(il, i+1)*mp(il, i+1) + (cpn(il, i)*t(il,i)+gz(il,i))*(mp(il,i)-mp(il,i+1)) - &
!IM121224           lv(il,i)*100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))*(evap(il,i+1)+evap(il,i))
!IM121224         hp(il, i) = p(il, i)/mp(il, i)
          updn(il, i) = updn(il, i+1)*mpdn(il, i+1) + u(il, i)*(mpdn(il,i)-mpdn(il,i+1))
          updn(il, i) = updn(il, i)/mpdn(il, i)
          vpdn(il, i) = vpdn(il, i+1)*mpdn(il, i+1) + v(il, i)*(mpdn(il,i)-mpdn(il,i+1))
          vpdn(il, i) = vpdn(il, i)/mpdn(il, i)

        ELSE ! if (mplus(il))

          IF (mpdn(il,i+1)>1.0E-16) THEN
!IM121224           IF (cvflag_grav) THEN
!IM121224             rp(il, i) = rp(il,i+1) + 100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
!IM121224                                      (evap(il,i+1)+evap(il,i))/mp(il,i+1)
!IM121224           ELSE
!IM121224             rp(il, i) = rp(il,i+1) + 5.*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
!IM121224                                      (evap(il,i+1)+evap(il,i))/mp(il, i+1)
!IM121224           END IF
!IM121224           hp(il, i) = hp(il,i+1) - lv(il,i)*100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1)) * &
!IM121224                                    (evap(il,i+1)+evap(il,i))/mp(il,i+1)
            updn(il, i) = updn(il, i+1)
            vpdn(il, i) = vpdn(il, i+1)
          END IF ! (mp(il,i+1).gt.1.0e-16)
        END IF ! (mplus(il)) else if (.not.mplus(il))

!IM121224 BEGIN ??? 
        rpdn(il, i) = amin1(rpdn(il,i), rs(il,i))
        rpdn(il, i) = max(rpdn(il,i), 0.0)
!IM121224 END ???
      END IF ! (i.lt.inb(il) .and. lwork(il))
    END DO
! ----------------------------------------------------------------

! ***       find tracer concentrations in precipitating downdraft     ***

!AC!      do j=1,ntra
!AC!       do il = 1,ncum
!AC!       if (i.lt.inb(il) .and. lwork(il)) then
!AC!c
!AC!         if(mplus(il))then
!AC!          trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
!AC!     :              +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
!AC!          trap(il,i,j)=trap(il,i,j)/mp(il,i)
!AC!         else ! if (mplus(il))
!AC!          if(mp(il,i+1).gt.1.0e-16)then
!AC!           trap(il,i,j)=trap(il,i+1,j)
!AC!          endif
!AC!         endif ! (mplus(il)) else if (.not.mplus(il))
!AC!c
!AC!        endif ! (i.lt.inb(il) .and. lwork(il))
!AC!       enddo
!AC!      end do

400 END DO
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

! ***                    end of downdraft loop                    ***

! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


  RETURN

END SUBROUTINE cv3_unsat24
