SUBROUTINE cv3_yield24(nloc, ncum, nd, na, ntra, ok_conserv_q, &
                     icb, inb, delt, &
                     t, rr, t_wake, rr_wake, s_wake, u, v, tra, &
                     gz, p, ph, h, h_wake, hp, lv, lf, cpn, lv_wake, lf_wake, cpn_wake, th, th_wake, &
                     ep, clw, qpreca, m, tp, mpdn, hpdn, rpdn, updn, vpdn, trap, &
                     wt, water, ice, evap, fondue, faci, b, sigd, &
                     ment, qent, hent, iflag_mix, uent, vent, &
                     nent, elij, traent, sig, &
                     tv, tvp, wghti, &
                     iflag, precip, Vprecip, Vprecipi, &     ! jyg: Vprecipi
                     ft, fr, fr_comp, fu, fv, ftra, &                 ! jyg
                     cbmf, upwd, dnwd, dnwd0, ma, mip, &
!!                     tls, tps,                             ! useless . jyg
                     qcondc, wd, &
                     ftd, fqd, qta, qtc, sigt, detrain, tau_cld_cv, coefw_cld_cv)

    USE print_control_mod, ONLY: lunout, prt_level
    USE add_phys_tend_mod, only : fl_cor_ebil

  IMPLICIT NONE

  include "cvthermo.h"
  include "cv3param.h"
  include "cvflag.h"
  include "conema3.h"

!inputs:
      INTEGER, INTENT (IN)                               :: iflag_mix
      INTEGER, INTENT (IN)                               :: ncum, nd, na, ntra, nloc
      LOGICAL, INTENT (IN)                               :: ok_conserv_q
      INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb
      REAL, INTENT (IN)                                  :: delt
      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, u, v
      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t_wake, rr_wake
      REAL, DIMENSION (nloc), INTENT (IN)                :: s_wake
      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)            :: gz, h, h_wake, hp
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: th, tp
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv, cpn, ep, clw
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv_wake, lf_wake, cpn_wake
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: hpdn, rpdn, updn, vpdn
      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: wt
      REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: trap
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: water, evap, b
      REAL, DIMENSION (nloc, na), INTENT (IN)            :: fondue, faci, ice
      REAL, DIMENSION (nloc, na, na), INTENT (IN)        :: qent, uent
      REAL, DIMENSION (nloc, na, na), INTENT (IN)        :: hent
      REAL, DIMENSION (nloc, na, na), INTENT (IN)        :: vent, elij
      INTEGER, DIMENSION (nloc, nd), INTENT (IN)         :: nent
      REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN)  :: traent
      REAL, DIMENSION (nloc, nd), INTENT (IN)            :: tv, tvp, wghti
      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:
      REAL, DIMENSION (nloc, na), INTENT (INOUT)         :: m, mpdn
      REAL, DIMENSION (nloc, na, na), INTENT (INOUT)     :: ment
      INTEGER, DIMENSION (nloc), INTENT (INOUT)          :: iflag
      REAL, DIMENSION (nloc, nd), INTENT (INOUT)         :: sig
      REAL, DIMENSION (nloc), INTENT (INOUT)             :: sigd
!
!outputs:
      REAL, DIMENSION (nloc), INTENT (OUT)               :: precip
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ft, fr, fu, fv , fr_comp
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ftd, fqd
      REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT)     :: ftra
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: upwd, dnwd, ma
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: dnwd0, mip
      REAL, DIMENSION (nloc, nd+1), INTENT (OUT)         :: Vprecip
      REAL, DIMENSION (nloc, nd+1), INTENT (OUT)         :: Vprecipi
!!      REAL tls(nloc, nd), tps(nloc, nd)                    ! useless . jyg
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qcondc                      ! cld
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qtc, sigt                   ! cld
      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: detrain                     ! Louis : pour le calcul de Klein du terme de variance qui dÃ©traine dans lenvironnement
      REAL, DIMENSION (nloc), INTENT (OUT)               :: wd                          ! gust
      REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmf
!
!local variables:
      INTEGER                                            :: i, k, il, n, j, num1
      REAL                                               :: rat, delti
      REAL                                               :: ax, bx, cx, dx, ex
      REAL                                               :: cpinv, cpinv_wake, rdcp, dpinv
      REAL                                               :: sigaq
      REAL, DIMENSION (nloc)                             ::  awat
      REAL, DIMENSION (nloc, nd)                         :: flux_h, flux_q    ! enthalpy and moisture fluxes
                                                                              ! due to unsaturated downdraughts
      REAL, DIMENSION (nloc, nd)                         :: lvcp, lfcp              ! , mke ! unused . jyg
      REAL, DIMENSION (nloc, nd)                         :: lvcp_wake, lfcp_wake    ! , mke ! unused . jyg
      REAL, DIMENSION (nloc)                             :: am, work, ad, amp1
!!      real up1(nloc), dn1(nloc)
      REAL, DIMENSION (nloc, nd, nd)                     :: up1, dn1
!jyg<
      REAL, DIMENSION (nloc, nd)                         :: up_to, up_from
      REAL, DIMENSION (nloc, nd)                         :: dn_to, dn_from
!>jyg
      REAL, DIMENSION (nloc)                             :: asum, bsum, csum, dsum
      REAL, DIMENSION (nloc)                             :: esum, fsum, gsum, hsum
      REAL, DIMENSION (nloc, nd)                         :: th_wake
      REAL, DIMENSION (nloc)                             :: alpha_qpos, alpha_qpos1
      REAL, DIMENSION (nloc, nd)                         :: qcond, nqcond, wa           ! cld
      REAL, DIMENSION (nloc, nd)                         :: siga, sax, mac              ! cld
      REAL, DIMENSION (nloc)                             :: sument
      REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
      REAL, DIMENSION (nloc, nd, nd)                     :: qdet
      REAL sumdq !jyg
! 
! -------------------------------------------------------------

! initialization:

  delti = 1.0/delt
! print*,'cv3_yield initialisation delt', delt

  DO il = 1, ncum
    precip(il) = 0.0
    wd(il) = 0.0 ! gust
  END DO

!   Fluxes are on a staggered grid : loops extend up to nl+1
  DO i = 1, nlp
    DO il = 1, ncum
      Vprecip(il, i) = 0.0
      Vprecipi(il, i) = 0.0                               ! jyg
      upwd(il, i) = 0.0
      dnwd(il, i) = 0.0
      dnwd0(il, i) = 0.0
      mip(il, i) = 0.0
    END DO
  END DO
  DO i = 1, nl
    DO il = 1, ncum
      ft(il, i) = 0.0
      fr(il, i) = 0.0
      fr_comp(il,i) = 0.0
      fu(il, i) = 0.0
      fv(il, i) = 0.0
      ftd(il, i) = 0.0
      fqd(il, i) = 0.0
      qcondc(il, i) = 0.0 ! cld
      qcond(il, i) = 0.0 ! cld
      qtc(il, i) = 0.0 ! cld
      qtment(il, i) = 0.0 ! cld
      sigment(il, i) = 0.0 ! cld
      sigt(il, i) = 0.0 ! cld
      qdet(il,i,:) = 0.0 ! cld
      detrain(il, i) = 0.0 ! cld
      nqcond(il, i) = 0.0 ! cld
    END DO
  END DO
! print*,'cv3_yield initialisation 2'
!AC!      do j=1,ntra
!AC!       do i=1,nd
!AC!        do il=1,ncum
!AC!          ftra(il,i,j)=0.0
!AC!        enddo
!AC!       enddo
!AC!      enddo
! print*,'cv3_yield initialisation 3'
  DO i = 1, nl
    DO il = 1, ncum
      lvcp(il, i) = lv(il, i)/cpn(il, i)
      lfcp(il, i) = lf(il, i)/cpn(il, i)
      lvcp_wake(il, i) = lv_wake(il, i)/cpn_wake(il, i)
      lfcp_wake(il, i) = lf_wake(il, i)/cpn_wake(il, i)
    END DO
  END DO



! ***  calculate surface precipitation in mm/day     ***

  DO il = 1, ncum
    IF (ep(il,inb(il))>=0.0001 .AND. iflag(il)<=1) THEN
      IF (cvflag_ice) THEN
        precip(il) = wt(il, 1)*sigd(il)*(water(il,1)+ice(il,1)) &
                              *86400.*1000./(rowl*grav)
      ELSE
        precip(il) = wt(il, 1)*sigd(il)*water(il, 1) &
                              *86400.*1000./(rowl*grav)
      END IF
    END IF
  END DO
! print*,'cv3_yield apres calcul precip'


! ===  calculate vertical profile of  precipitation in kg/m2/s  ===

  DO i = 1, nl
    DO il = 1, ncum
      IF (ep(il,inb(il))>=0.0001 .AND. i<=inb(il) .AND. iflag(il)<=1) THEN
        IF (cvflag_ice) THEN
          Vprecip(il, i) = wt(il, i)*sigd(il)*(water(il,i)+ice(il,i))/grav
          Vprecipi(il, i) = wt(il, i)*sigd(il)*ice(il,i)/grav                   ! jyg
        ELSE
          Vprecip(il, i) = wt(il, i)*sigd(il)*water(il, i)/grav
          Vprecipi(il, i) = 0.                                                  ! jyg
        END IF
      END IF
    END DO
  END DO


! ***  Calculate downdraft velocity scale    ***
! ***  NE PAS UTILISER POUR L'INSTANT ***

!!      do il=1,ncum
!!        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) &
!!                                       /(sigd(il)*p(il,icb(il)))
!!      enddo


! ***  calculate tendencies of lowest level potential temperature  ***
! ***                      and mixing ratio                        ***

  DO il = 1, ncum
    work(il) = 1.0/(ph(il,1)-ph(il,2))
    cbmf(il) = 0.0
  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
      IF (k>=icb(il)) THEN
        cbmf(il) = cbmf(il) + m(il, k)
      END IF
    END DO
  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
  DO il = 1, ncum
    am(il) = cbmf(il)*wghti(il, 1)
  END DO

!1712 BEGIN
!
!  Vertical fluxes of moisture and enthalpy due to unsaturated downdraughts
  DO il = 1, ncum
   flux_h(il,1) = 0.
   flux_q(il,1) = 0.
  END DO 
  DO j= 2, nl
    DO il = 1, ncum
! Gk = flux_h
! Fk = flux_q
      flux_h(il,j) = mpdn(il,j)*(hpdn(il,j)-h_wake(il,j-1))
      flux_q(il,j) = mpdn(il,j)*(rpdn(il,j)-rr_wake(il,j-1))
    END DO 
  END DO 
!1712 END 

  DO il = 1, ncum
    IF (iflag(il)<=1) THEN
! convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
!JYG  Correction pour conserver l'eau
! cc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2))          !precip
!1712     IF (cvflag_ice) THEN
!1712       ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1) - &
!1712                    lfcp(il, 1)*sigd(il)*evap(il, 1)*faci(il, 1) - &
!1712                    lfcp(il, 1)*sigd(il)*(fondue(il,1)*wt(il,1)) / &
!1712                      (100.*(ph(il,1)-ph(il,2)))                             !precip
!1712     ELSE
!1712       ft(il, 1) = -lvcp(il, 1)*sigd(il)*evap(il, 1)
!1712     END IF
!1712
!1712     ft(il, 1) = ft(il, 1) - 0.009*grav*sigd(il)*mpdn(il, 2)*t_wake(il, 1)*b(il, 1)*work(il)
!1712
!1712     IF (cvflag_ice) THEN
!1712       ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
!1712                                    (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1) + &
!1712                               0.01*sigd(il)*wt(il, 1)*(ci-cpd)*ice(il, 2) * &
!1712                                    (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
!1712     ELSE
!1712       ft(il, 1) = ft(il, 1) + 0.01*sigd(il)*wt(il, 1)*(cl-cpd)*water(il, 2) * &
!1712                                    (t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il, 1)
!1712     END IF
!lfcp=lf/cp
!lvcp=lv/cp
!
! Formule VS code
!Mk = - mpdn(k)
!qk = rpdn(k)
!hk = hpdn(k)
!cpk = cpn_wake(k)
!qtilde(k) = rr_wake(k)
!htilde(k) = h_wake(k)
!sigma_s(k)*Ek(k) = evap(k)
!Sk(k) = fondue(k)*wt(k) ; wt=vitesse terminale de chute des precipitations
!
      ft(il,1)= 1./cpn_wake(il, i)*grav/(100.*(ph(il,1)-ph(il,2)))*(-flux_h(il, 2)) -  &
                sigd(il)*((lvcp_wake(il, 1)+faci(il, 1)*lfcp_wake(il, 1))*evap(il, 1) + &
                          lfcp_wake(il, 1)*fondue(il, 1)*wt(il,1) )

      ftd(il, 1) = ft(il, 1)                                                  ! fin precip

      IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect
!jyg<
        IF (fl_cor_ebil >= 2) THEN
          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
                    ((t(il,2)-t(il,1))*cpn(il,2)+gz(il,2)-gz(il,1))/cpn(il,1)
        ELSE
          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*am(il) * &
                    (t(il,2)-t(il,1)+(gz(il,2)-gz(il,1))/cpn(il,1))
        ENDIF
!>jyg
    END IF ! iflag
  END DO


  DO j = 2, nl
    IF (iflag_mix>0) THEN
      DO il = 1, ncum
! FH WARNING a modifier :
        cpinv = 0.
! cpinv=1.0/cpn(il,1)
        IF (j<=inb(il) .AND. iflag(il)<=1) THEN
          ft(il, 1) = ft(il, 1) + 0.01*grav*work(il)*ment(il, j, 1) * &
                     (hent(il,j,1)-h(il,1)+t(il,1)*(cpv-cpd)*(rr(il,1)-qent(il,j,1)))*cpinv
        END IF ! j
      END DO
    END IF
  END DO
! fin sature


  DO il = 1, ncum
    IF (iflag(il)<=1) THEN
!JYG1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
!1401     fr(il, 1) = 0.01*grav*mpdn(il, 2)*(rpdn(il,2)-rr_wake(il,1))*work(il) + &
!1401                 sigd(il)*evap(il, 1)
          fr(il, 1) = grav/(100.*(ph(il,1)-ph(il,2)))* (-flux_q(il, 2))+ &
                 sigd(il)*evap(il, 1)
!!!                  sigd(il)*0.5*(evap(il,1)+evap(il,2))

      fqd(il, 1) = fr(il, 1) !precip

      fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)        !sature

      fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mpdn(il,2)*(updn(il,2)-u(il,1)) + &
                                                  am(il)*(u(il,2)-u(il,1)))
      fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mpdn(il,2)*(vpdn(il,2)-v(il,1)) + &
                                                  am(il)*(v(il,2)-v(il,1)))
    END IF ! iflag
  END DO ! il


!AC!     do j=1,ntra
!AC!      do il=1,ncum
!AC!       if (iflag(il) .le. 1) then
!AC!       if (cvflag_grav) then
!AC!        ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
!AC!       else
!AC!        ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
!AC!    :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
!AC!    :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
!AC!       endif
!AC!       endif  ! iflag
!AC!      enddo
!AC!     enddo

  DO j = 2, nl
    DO il = 1, ncum
      IF (j<=inb(il) .AND. iflag(il)<=1) THEN
        fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
        fr_comp(il,1) = fr_comp(il,1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
        fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il,j,1)-u(il,1))
        fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il,j,1)-v(il,1))
      END IF ! j
    END DO
  END DO

!AC!      do k=1,ntra
!AC!       do j=2,nl
!AC!        do il=1,ncum
!AC!         if (j.le.inb(il) .and. iflag(il) .le. 1) then
!AC!
!AC!          if (cvflag_grav) then
!AC!           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
!AC!          else
!AC!           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
!AC!     :                *(traent(il,j,1,k)-tra(il,1,k))
!AC!          endif
!AC!
!AC!         endif
!AC!        enddo
!AC!       enddo
!AC!      enddo
! print*,'cv3_yield apres ft'

!jyg<
!-----------------------------------------------------------
           IF (ok_optim_yield) THEN                       !|
!-----------------------------------------------------------
!
!***                                                      ***
!***    Compute convective mass fluxes upwd and dnwd      ***

!
! =================================================
!              upward fluxes                      |
! ------------------------------------------------
!
upwd(:,:) = 0.
up_to(:,:) = 0.
up_from(:,:) = 0.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  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
    IF (i<=inb(il)) THEN
      up_to(il,i) = m(il,i)
    ENDIF
  ENDDO
  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 i = 1, nl
  DO il = 1, ncum
    IF (i<=inb(il)) THEN
      up_from(il,i) = cbmf(il)*wghti(il,i)
    ENDIF
  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


  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
    DO il = 1, ncum
      IF (j<=inb(il)) THEN
!!        dn_to(il,i) = dn_to(il,i) + ment(il,j,i)       !jyg,20220202
        dn_to(il,i) = dn_to(il,i) - ment(il,j,i)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
DO j = 1, nl
  DO i = j+1, nl
    DO il = 1, ncum
      IF (i<=inb(il)) THEN
!!        dn_from(il,i) = dn_from(il,i) + ment(il,i,j)   !jyg,20220202
        dn_from(il,i) = dn_from(il,i) - ment(il,i,j)
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
! The difference between dnwd(il,i) and dnwd(il,i+1) is due to downdrafts ending in layer 
!(i) (theses drafts cross interface (i+1) but not interface(i)) and to downdrafts 
!starting from layer (i) (theses drafts cross interface (i) but not interface(i+1)): 
!
DO i = nl-1, 1, -1
  DO il = 1, ncum
!!    dnwd(il,i) = max(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i)) !jyg,20220202
    dnwd(il,i) = min(0., dnwd(il,i+1) - dn_to(il,i) + dn_from(il,i))
  ENDDO
ENDDO
! =================================================
!
!-----------------------------------------------------------
        ENDIF !(ok_optim_yield)                           !|
!-----------------------------------------------------------
!>jyg

! ***  calculate tendencies of potential temperature and mixing ratio  ***
! ***               at levels above the lowest level                   ***

! ***  first find the net saturated updraft and downdraft mass fluxes  ***
! ***                      through each level                          ***


!jyg<
!!  DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1?
  DO i = 2, nl
!>jyg 

    num1 = 0
    DO il = 1, ncum
      IF (i<=inb(il) .AND. iflag(il)<=1) num1 = num1 + 1
    END DO
    IF (num1<=0) GO TO 500

!
!jyg<
!-----------------------------------------------------------
           IF (ok_optim_yield) THEN                       !|
!-----------------------------------------------------------
DO il = 1, ncum
   amp1(il) = upwd(il,i+1)
   ad(il) = dnwd(il,i)
ENDDO
!-----------------------------------------------------------
        ELSE !(ok_optim_yield)                            !|
!-----------------------------------------------------------
!>jyg
    DO il = 1,ncum
      amp1(il) = 0.
      ad(il) = 0.
    ENDDO

    DO k = 1, nl + 1
      DO il = 1, ncum
        IF (i>=icb(il)) THEN
          IF (k>=i+1 .AND. k<=(inb(il)+1)) THEN
            amp1(il) = amp1(il) + m(il, k)
          END IF
        ELSE
! AMP1 is the part of cbmf taken from layers I and lower
          IF (k<=i) THEN
            amp1(il) = amp1(il) + cbmf(il)*wghti(il, k)
          END IF
        END IF
      END DO
    END DO

    DO j = i + 1, nl + 1         
       DO k = 1, i
          !yor! reverted j and k loops 
          DO il = 1, ncum
!yor!        IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN ! the second condition implies the first !
             IF (j<=(inb(il)+1)) THEN  
                amp1(il) = amp1(il) + ment(il, k, j)
             END IF
          END DO
       END DO
    END DO

    DO k = 1, i - 1
!jyg<
!!      DO j = i, nl + 1 ! newvecto: nl au lieu nl+1?
      DO j = i, nl
!>jyg
        DO il = 1, ncum
!yor!        IF (i<=inb(il) .AND. j<=inb(il)) THEN ! the second condition implies the 1st !
             IF (j<=inb(il)) THEN   
            ad(il) = ad(il) + ment(il, j, k)
          END IF
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------
        ENDIF !(ok_optim_yield)                           !|
!-----------------------------------------------------------
!
!!   print *,'yield, i, amp1, ad', i, amp1(1), ad(1)

    DO il = 1, ncum
      IF (i<=inb(il) .AND. iflag(il)<=1) THEN
        dpinv = 1.0/(ph(il,i)-ph(il,i+1))
        cpinv = 1.0/cpn(il, i)
        cpinv_wake = 1.0/cpn_wake(il, i)

! convect3      if((0.1*dpinv*amp1).ge.delti)iflag(il)=4
        IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto

! precip
! cc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
!140125       IF (cvflag_ice) THEN
!140125         ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i) - &
!140125                      sigd(il)*lfcp(il, i)*evap(il, i)*faci(il, i) - &
!140125                      sigd(il)*lfcp(il, i)*fondue(il, i)*wt(il, i)/(100.*(p(il,i-1)-p(il,i)))
!140125       ELSE
!140125         ft(il, i) = -sigd(il)*lvcp(il, i)*evap(il, i)
!140125       END IF

!140125       rat = cpn(il, i-1)*cpinv

!140125       ft(il, i) = ft(il, i) - 0.009*grav*sigd(il) * &
!140125                    (mpdn(il,i+1)*t_wake(il,i)*b(il,i)-mpdn(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
!140125       IF (cvflag_ice) THEN
!140125         ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
!140125                                      (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv + &
!140125                                 0.01*sigd(il)*wt(il, i)*(ci-cpd)*ice(il, i+1) * &
!140125                                      (t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
!140125       ELSE
!140125         ft(il, i) = ft(il, i) + 0.01*sigd(il)*wt(il, i)*(cl-cpd)*water(il, i+1) * &
!140125                                      (t_wake(il,i+1)-t_wake(il,i))*dpinv* &
!140125           cpinv
!140125       END IF

        ft(il, i) = cpinv_wake*grav*0.01*dpinv*(flux_h(il, i)-flux_h(il, i+1)) - &
                    cpinv_wake*sigd(il)*((lvcp_wake(il, i)+faci(il, i)*lfcp_wake(il, i))*evap(il, i) + &
                    lfcp_wake(il, i)*fondue(il, i)*wt(il,i) )
   
        ftd(il, i) = ft(il, i)
! fin precip

! sature
!jyg<
        IF (fl_cor_ebil >= 2) THEN
          ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
              ( amp1(il)*( (t(il,i+1)-t(il,i))*cpn(il,i+1) + gz(il,i+1)-gz(il,i))*cpinv - &
                ad(il)*( (t(il,i)-t(il,i-1))*cpn(il,i-1) + gz(il,i)-gz(il,i-1))*cpinv)
        ELSE
          ft(il, i) = ft(il, i) + 0.01*grav*dpinv * &
                     (amp1(il)*(t(il,i+1)-t(il,i) + (gz(il,i+1)-gz(il,i))*cpinv) - &
                      ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
        ENDIF
!>jyg


        IF (iflag_mix==0) THEN
          ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i) + &
                                    t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
        END IF
!
! sb: on ne fait pas encore la correction permettant de mieux
! conserver l'eau:
!JYG: correction permettant de mieux conserver l'eau:
! cc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
!140125       fr(il, i) = sigd(il)*evap(il, i) + 0.01*grav*(mpdn(il,i+1)*(rpdn(il,i+1)-rr_wake(il,i)) - &
!140125                                                     mpdn(il,i)*(rpdn(il,i)-rr_wake(il,i-1)))*dpinv
        fr(il, i) = grav*0.01*dpinv*(flux_q(il, i) - flux_q(il, i+1)) + &
                    sigd(il)*evap(il,i) 
        fqd(il, i) = fr(il, i)                                                                     ! precip

        fu(il, i) = 0.01*grav*(mpdn(il,i+1)*(updn(il,i+1)-u(il,i)) - &
                               mpdn(il,i)*(updn(il,i)-u(il,i-1)))*dpinv
        fv(il, i) = 0.01*grav*(mpdn(il,i+1)*(vpdn(il,i+1)-v(il,i)) - &
                               mpdn(il,i)*(vpdn(il,i)-v(il,i-1)))*dpinv


        fr(il, i) = fr(il, i) + 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) - &
                                                 ad(il)*(rr(il,i)-rr(il,i-1)))
        fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i)) - &
                                                 ad(il)*(u(il,i)-u(il,i-1)))
        fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il,i)) - &
                                                 ad(il)*(v(il,i)-v(il,i-1)))

      END IF ! i
    END DO

!AC!      do k=1,ntra
!AC!       do il=1,ncum
!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
!AC!         cpinv=1.0/cpn(il,i)
!AC!         if (cvflag_grav) then
!AC!           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
!AC!         else
!AC!           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
!AC!     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
!AC!     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
!AC!         endif
!AC!        endif
!AC!       enddo
!AC!      enddo

    DO k = 1, i - 1

      DO il = 1, ncum
        awat(il) = elij(il, k, i) - (1.-ep(il,i))*clw(il, i)
        awat(il) = max(awat(il), 0.0)
      END DO

      IF (iflag_mix/=0) THEN
        DO il = 1, ncum
          IF (i<=inb(il) .AND. iflag(il)<=1) THEN
            dpinv = 1.0/(ph(il,i)-ph(il,i+1))
            cpinv = 1.0/cpn(il, i)
            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
                 (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-qent(il,k,i)))*cpinv
!
!
          END IF ! i
        END DO
      END IF

      DO il = 1, ncum
        IF (i<=inb(il) .AND. iflag(il)<=1) THEN
          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
          cpinv = 1.0/cpn(il, i)
          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
                                                       (qent(il,k,i)-awat(il)-rr(il,i))
          fr_comp(il,i) = fr_comp(il,i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-awat(il)-rr(il,i))
          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))

! (saturated updrafts resulting from mixing)                                   ! cld
          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
          qdet(il,k,i) = (qent(il,k,i)-awat(il))                               ! cld Louis : specific humidity in detraining water
          qtment(il, i) = qtment(il, i) + qent(il,k,i)                         ! cld
          nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
        END IF ! i
      END DO
    END DO

!AC!      do j=1,ntra
!AC!       do k=1,i-1
!AC!        do il=1,ncum
!AC!         if (i.le.inb(il) .and. iflag(il) .le. 1) then
!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
!AC!          cpinv=1.0/cpn(il,i)
!AC!          if (cvflag_grav) then
!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
!AC!          else
!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
!AC!     :        *(traent(il,k,i,j)-tra(il,i,j))
!AC!          endif
!AC!         endif
!AC!        enddo
!AC!       enddo
!AC!      enddo

!jyg<
!!    DO k = i, nl + 1
    DO k = i, nl
!>jyg

      IF (iflag_mix/=0) THEN
        DO il = 1, ncum
          IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
            dpinv = 1.0/(ph(il,i)-ph(il,i+1))
            cpinv = 1.0/cpn(il, i)
            ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
                  (hent(il,k,i)-h(il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,k,i)))*cpinv


          END IF ! i
        END DO
      END IF

      DO il = 1, ncum
        IF (i<=inb(il) .AND. k<=inb(il) .AND. iflag(il)<=1) THEN
          dpinv = 1.0/(ph(il,i)-ph(il,i+1))
          cpinv = 1.0/cpn(il, i)

          fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-rr(il,i))
          fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
          fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
        END IF ! i and k
      END DO
    END DO

!AC!      do j=1,ntra
!AC!       do k=i,nl+1
!AC!        do il=1,ncum
!AC!         if (i.le.inb(il) .and. k.le.inb(il)
!AC!     $                .and. iflag(il) .le. 1) then
!AC!          dpinv=1.0/(ph(il,i)-ph(il,i+1))
!AC!          cpinv=1.0/cpn(il,i)
!AC!          if (cvflag_grav) then
!AC!           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
!AC!     :         *(traent(il,k,i,j)-tra(il,i,j))
!AC!          else
!AC!           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
!AC!     :             *(traent(il,k,i,j)-tra(il,i,j))
!AC!          endif
!AC!         endif ! i and k
!AC!        enddo
!AC!       enddo
!AC!      enddo

! sb: interface with the cloud parameterization:                               ! cld

    DO k = i + 1, nl
      DO il = 1, ncum
        IF (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN               ! cld
! (saturated downdrafts resulting from mixing)                                 ! cld
          qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
          qdet(il,k,i) = qent(il,k,i)                                          ! cld Louis : specific humidity in detraining water
          qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
          nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
        END IF ! cld
      END DO ! cld
    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
        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
!jyg<   Bug correction 20180620
!      PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
!!        qtment(il, i) = qent(il,k,i) + qtment(il,i)                            ! cld
        qdet(il,i,i) = qent(il,i,i)                                            ! cld Louis : specific humidity in detraining water
        qtment(il, i) = qent(il,i,i) + qtment(il,i)                            ! cld
!>jyg
        nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
      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
      IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN            ! cld
        qcond(il, i) = qcond(il, i)/nqcond(il, i)                              ! cld
        qtment(il, i) = qtment(il,i)/nqcond(il, i)                             ! cld
      END IF                                                                   ! cld
    END DO

!AC!      do j=1,ntra
!AC!       do il=1,ncum
!AC!        if (i.le.inb(il) .and. iflag(il) .le. 1) then
!AC!         dpinv=1.0/(ph(il,i)-ph(il,i+1))
!AC!         cpinv=1.0/cpn(il,i)
!AC!
!AC!         if (cvflag_grav) then
!AC!          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
!AC!         else
!AC!          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
!AC!     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
!AC!     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
!AC!         endif
!AC!        endif ! i
!AC!       enddo
!AC!      enddo


500 END DO

!JYG<
!Conservation de l'eau
!   sumdq = 0.
!   DO k = 1, nl
!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
!   END DO
!   PRINT *, 'cv3_yield, apres 500, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
!JYG>
! ***   move the detrainment at level inb down to level inb-1   ***
! ***        in such a way as to preserve the vertically        ***
! ***          integrated enthalpy and water tendencies         ***

! Correction bug le 18-03-09
  DO il = 1, ncum
    IF (iflag(il)<=1) THEN
      ax = 0.01*grav*ment(il, inb(il), inb(il))* &
           (hp(il,inb(il))-h(il,inb(il))+t(il,inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il),inb(il))))/ &
                                (cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
      ft(il, inb(il)) = ft(il, inb(il)) - ax
      ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
                              (cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il,inb(il))))

      bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb(il)))/ &
                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
      fr(il, inb(il)) = fr(il, inb(il)) - bx
      fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))

      cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il)))/ &
                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
      fu(il, inb(il)) = fu(il, inb(il)) - cx
      fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))

      dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il)))/ &
                                                 (ph(il,inb(il))-ph(il,inb(il)+1))
      fv(il, inb(il)) = fv(il, inb(il)) - dx
      fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+1))/ &
                                                 (ph(il,inb(il)-1)-ph(il,inb(il)))
    END IF !iflag
  END DO

!JYG<
!Conservation de l'eau
!   sumdq = 0.
!   DO k = 1, nl
!     sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
!   END DO
!   PRINT *, 'cv3_yield, apres 503, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
!JYG>

!AC!      do j=1,ntra
!AC!       do il=1,ncum
!AC!        IF (iflag(il) .le. 1) THEN
!AC!	IF (cvflag_grav) then
!AC!        ex=0.01*grav*ment(il,inb(il),inb(il))
!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
!AC!	else
!AC!        ex=0.1*ment(il,inb(il),inb(il))
!AC!     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
!AC!     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
!AC!        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
!AC!        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
!AC!     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
!AC!     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
!AC!        ENDIF   !cvflag grav
!AC!        ENDIF    !iflag
!AC!       enddo
!AC!      enddo


! ***    homogenize tendencies below cloud base    ***


  DO il = 1, ncum
    asum(il) = 0.0
    bsum(il) = 0.0
    csum(il) = 0.0
    dsum(il) = 0.0
    esum(il) = 0.0
    fsum(il) = 0.0
    gsum(il) = 0.0
    hsum(il) = 0.0
  END DO

!do i=1,nl
!do il=1,ncum
!th_wake(il,i)=t_wake(il,i)*(1000.0/p(il,i))**rdcp
!enddo
!enddo

  DO i = 1, nl
    DO il = 1, ncum
      IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
!jyg  Saturated part : use T profile
        asum(il) = asum(il) + (ft(il,i)-ftd(il,i))*(ph(il,i)-ph(il,i+1))
!jyg<20140311
!Correction pour conserver l eau
        IF (ok_conserv_q) THEN
          bsum(il) = bsum(il) + (fr(il,i)-fqd(il,i))*(ph(il,i)-ph(il,i+1))
          csum(il) = csum(il) + (ph(il,i)-ph(il,i+1))

        ELSE
          bsum(il)=bsum(il)+(fr(il,i)-fqd(il,i))*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
                            (ph(il,i)-ph(il,i+1))
          csum(il)=csum(il)+(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))* &
                            (ph(il,i)-ph(il,i+1))
        ENDIF ! (ok_conserv_q)
!jyg>
        dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i)
!jyg  Unsaturated part : use T_wake profile
        esum(il) = esum(il) + ftd(il, i)*(ph(il,i)-ph(il,i+1))
!jyg<20140311
!Correction pour conserver l eau
        IF (ok_conserv_q) THEN
          fsum(il) = fsum(il) + fqd(il, i)*(ph(il,i)-ph(il,i+1))
          gsum(il) = gsum(il) + (ph(il,i)-ph(il,i+1))
        ELSE
          fsum(il)=fsum(il)+fqd(il,i)*(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
                            (ph(il,i)-ph(il,i+1))
          gsum(il)=gsum(il)+(lv(il,i)+(cl-cpd)*(t_wake(il,i)-t_wake(il,1)))* &
                            (ph(il,i)-ph(il,i+1))
        ENDIF ! (ok_conserv_q)
!jyg>
        hsum(il) = hsum(il) + t_wake(il, i)*(ph(il,i)-ph(il,i+1))/th_wake(il, i)
      END IF
    END DO
  END DO

!!!!      do 700 i=1,icb(il)-1
  IF (ok_homo_tend) THEN
    DO i = 1, nl
      DO il = 1, ncum
        IF (i<=(icb(il)-1) .AND. iflag(il)<=1) THEN
          ftd(il, i) = esum(il)*t_wake(il, i)/(th_wake(il,i)*hsum(il))
          fqd(il, i) = fsum(il)/gsum(il)
          ft(il, i) = ftd(il, i) + asum(il)*t(il, i)/(th(il,i)*dsum(il))
          fr(il, i) = fqd(il, i) + bsum(il)/csum(il)
        END IF
      END DO
    END DO
  ENDIF

!jyg<
!Conservation de l'eau
!!  sumdq = 0.
!!  DO k = 1, nl
!!    sumdq = sumdq + fr(1, k)*100.*(ph(1,k)-ph(1,k+1))/grav
!!  END DO
!!  PRINT *, 'cv3_yield, apres hom, sum(dq), precip, somme ', sumdq, Vprecip(1, 1), sumdq + vprecip(1, 1)
!jyg>


! ***   Check that moisture stays positive. If not, scale tendencies
! in order to ensure moisture positivity
  DO il = 1, ncum
    alpha_qpos(il) = 1.
    IF (iflag(il)<=1) THEN
      IF (fr(il,1)<=0.) THEN
        alpha_qpos(il) = max(alpha_qpos(il), (-delt*fr(il,1))/(s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
      END IF
    END IF
  END DO
  DO i = 2, nl
    DO il = 1, ncum
      IF (iflag(il)<=1) THEN
        IF (fr(il,i)<=0.) THEN
          alpha_qpos1(il) = max(1., (-delt*fr(il,i))/(s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
          IF (alpha_qpos1(il)>=alpha_qpos(il)) alpha_qpos(il) = alpha_qpos1(il)
        END IF
      END IF
    END DO
  END DO
  DO il = 1, ncum
    IF (iflag(il)<=1 .AND. alpha_qpos(il)>1.001) THEN
      alpha_qpos(il) = alpha_qpos(il)*1.1
    END IF
  END DO
!
    IF (prt_level .GE. 5) THEN
      print *,' CV3_YIELD : alpha_qpos ',alpha_qpos(1)
    ENDIF
!
  DO il = 1, ncum
    IF (iflag(il)<=1) THEN
      sigd(il) = sigd(il)/alpha_qpos(il)
      precip(il) = precip(il)/alpha_qpos(il)
      cbmf(il) = cbmf(il)/alpha_qpos(il)
    END IF
  END DO
  DO i = 1, nl
    DO il = 1, ncum
      IF (iflag(il)<=1) THEN
        fr(il, i) = fr(il, i)/alpha_qpos(il)
        ft(il, i) = ft(il, i)/alpha_qpos(il)
        fqd(il, i) = fqd(il, i)/alpha_qpos(il)
        ftd(il, i) = ftd(il, i)/alpha_qpos(il)
        fu(il, i) = fu(il, i)/alpha_qpos(il)
        fv(il, i) = fv(il, i)/alpha_qpos(il)
        m(il, i) = m(il, i)/alpha_qpos(il)
        mpdn(il, i) = mpdn(il, i)/alpha_qpos(il)
        Vprecip(il, i) = Vprecip(il, i)/alpha_qpos(il)
        Vprecipi(il, i) = Vprecipi(il, i)/alpha_qpos(il)                     ! jyg
      END IF
    END DO
  END DO
!jyg<
!-----------------------------------------------------------
           IF (ok_optim_yield) THEN                       !|
!-----------------------------------------------------------
  DO i = 1, nl
    DO il = 1, ncum
      IF (iflag(il)<=1) THEN
        upwd(il, i) = upwd(il, i)/alpha_qpos(il)
        dnwd(il, i) = dnwd(il, i)/alpha_qpos(il)
      END IF
    END DO
  END DO
!-----------------------------------------------------------
        ENDIF !(ok_optim_yield)                           !|
!-----------------------------------------------------------
!>jyg
  DO j = 1, nl !yor! inverted i and j loops
     DO i = 1, nl
      DO il = 1, ncum
        IF (iflag(il)<=1) THEN
          ment(il, i, j) = ment(il, i, j)/alpha_qpos(il)
        END IF
      END DO
    END DO
  END DO

!AC!      DO j = 1,ntra
!AC!      DO i = 1,nl
!AC!       DO il = 1,ncum
!AC!        IF (iflag(il) .le. 1) THEN
!AC!         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
!AC!        ENDIF
!AC!       ENDDO
!AC!      ENDDO
!AC!      ENDDO


! ***           reset counter and return           ***

! Reset counter only for points actually convective (jyg)
! In order take into account the possibility of changing the compression,
! reset m, sig and w0 to zero for non-convecting points.
  DO il = 1, ncum
    IF (iflag(il) < 3) THEN
      sig(il, nd) = 2.0
    ENDIF
  END DO


  DO i = 1, nl
    DO il = 1, ncum
      dnwd0(il, i) = -mpdn(il, i)
    END DO
  END DO
!jyg<  (loops stop at nl)
!!  DO i = nl + 1, nd
!!    DO il = 1, ncum
!!      dnwd0(il, i) = 0.
!!    END DO
!!  END DO
!>jyg


!jyg<
!-----------------------------------------------------------
           IF (.NOT.ok_optim_yield) THEN                  !|
!-----------------------------------------------------------
  DO i = 1, nl
    DO il = 1, ncum
      upwd(il, i) = 0.0
      dnwd(il, i) = 0.0
    END DO
  END DO

!!  DO i = 1, nl                                           ! useless; jyg
!!    DO il = 1, ncum                                      ! useless; jyg
!!      IF (i>=icb(il) .AND. i<=inb(il)) THEN              ! useless; jyg
!!        upwd(il, i) = 0.0                                ! useless; jyg
!!        dnwd(il, i) = 0.0                                ! useless; jyg
!!      END IF                                             ! useless; jyg
!!    END DO                                               ! useless; jyg
!!  END DO                                                 ! useless; jyg

  DO i = 1, nl
    DO k = 1, nl
      DO il = 1, ncum
        up1(il, k, i) = 0.0
        dn1(il, k, i) = 0.0
      END DO
    END DO
  END DO

!yor! commented original
!  DO i = 1, nl
!    DO k = i, nl
!      DO n = 1, i - 1
!        DO il = 1, ncum
!          IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN
!            up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
!            dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
!          END IF
!        END DO
!      END DO
!    END DO
!  END DO
!yor! replaced with
  DO i = 1, nl
    DO k = i, nl
      DO n = 1, i - 1
        DO il = 1, ncum
          IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor ! as i always <= k
             up1(il, k, i) = up1(il, k, i) + ment(il, n, k)
          END IF
        END DO
      END DO
    END DO
  END DO
  DO i = 1, nl
    DO n = 1, i - 1
      DO k = i, nl
        DO il = 1, ncum
          IF (i>=icb(il) .AND. k<=inb(il)) THEN ! yor !  i always <= k
             dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n)
          END IF
        END DO
      END DO
    END DO
  END DO
!yor! end replace 

  DO i = 1, nl
    DO k = 1, nl
      DO il = 1, ncum
        IF (i>=icb(il)) THEN
          IF (k>=i .AND. k<=(inb(il))) THEN
            upwd(il, i) = upwd(il, i) + m(il, k)
          END IF
        ELSE
          IF (k<i) THEN
            upwd(il, i) = upwd(il, i) + cbmf(il)*wghti(il, k)
          END IF
        END IF
! c        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
      END DO
    END DO
  END DO

  DO i = 2, nl
    DO k = i, nl
      DO il = 1, ncum
! test         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
        IF (i<=inb(il) .AND. k<=inb(il)) THEN
          upwd(il, i) = upwd(il, i) + up1(il, k, i)
          dnwd(il, i) = dnwd(il, i) + dn1(il, k, i)
        END IF
! c         print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
      END DO
    END DO
  END DO


!!!!      DO il=1,ncum
!!!!      do i=icb(il),inb(il)
!!!!
!!!!      upwd(il,i)=0.0
!!!!      dnwd(il,i)=0.0
!!!!      do k=i,inb(il)
!!!!      up1=0.0
!!!!      dn1=0.0
!!!!      do n=1,i-1
!!!!      up1=up1+ment(il,n,k)
!!!!      dn1=dn1-ment(il,k,n)
!!!!      enddo
!!!!      upwd(il,i)=upwd(il,i)+m(il,k)+up1
!!!!      dnwd(il,i)=dnwd(il,i)+dn1
!!!!      enddo
!!!!      enddo
!!!!
!!!!      ENDDO

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

! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! determination de la variation de flux ascendant entre
! deux niveau non dilue mip
! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  DO i = 1, nl
    DO il = 1, ncum
      mip(il, i) = m(il, i)
    END DO
  END DO

!jyg<  (loops stop at nl)
!!  DO i = nl + 1, nd
!!    DO il = 1, ncum
!!      mip(il, i) = 0.
!!    END DO
!!  END DO
!>jyg


! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! icb represente de niveau ou se trouve la
! base du nuage , et inb le top du nuage
! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!!  DO i = 1, nd                                  ! unused . jyg
!!    DO il = 1, ncum                             ! unused . jyg
!!      mke(il, i) = upwd(il, i) + dnwd(il, i)    ! unused . jyg
!!    END DO                                      ! unused . jyg
!!  END DO                                        ! unused . jyg

!!  DO i = 1, nd                                                                 ! unused . jyg
!!    DO il = 1, ncum                                                            ! unused . jyg
!!      rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il,i))+rr(il,i)*cpv) ! unused . jyg
!!      tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp                             ! unused . jyg
!!      tps(il, i) = tp(il, i)                                                   ! unused . jyg
!!    END DO                                                                     ! unused . jyg
!!  END DO                                                                       ! unused . jyg


! *** diagnose the in-cloud mixing ratio   ***                       ! cld
! ***           of condensed water         ***                       ! cld
!! cld                                                               
                                                                     
  DO i = 1, nl+1                                                     ! cld
    DO il = 1, ncum                                                  ! cld
      mac(il, i) = 0.0                                               ! cld
      wa(il, i) = 0.0                                                ! cld
      siga(il, i) = 0.0                                              ! cld
      sax(il, i) = 0.0                                               ! cld
    END DO                                                           ! cld
  END DO                                                             ! cld
                                                                     
  DO i = minorig, nl                                                 ! cld
    DO k = i + 1, nl + 1                                             ! cld
      DO il = 1, ncum                                                ! cld
        IF (i<=inb(il) .AND. k<=(inb(il)+1) .AND. iflag(il)<=1) THEN ! cld
          mac(il, i) = mac(il, i) + m(il, k)                         ! cld
        END IF                                                       ! cld
      END DO                                                         ! cld
    END DO                                                           ! cld
  END DO                                                             ! cld

  DO i = 1, nl                                                       ! cld
    DO j = 1, i                                                      ! cld
      DO il = 1, ncum                                                ! cld
        IF (i>=icb(il) .AND. i<=(inb(il)-1) &                        ! cld
            .AND. j>=icb(il) .AND. iflag(il)<=1) THEN                ! cld
          sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) &       ! cld
            *(ph(il,j)-ph(il,j+1))/p(il, j)                          ! cld
        END IF                                                       ! cld
      END DO                                                         ! cld
    END DO                                                           ! cld
  END DO                                                             ! cld

  DO i = 1, nl                                                       ! cld
    DO il = 1, ncum                                                  ! cld
      IF (i>=icb(il) .AND. i<=(inb(il)-1) &                          ! cld
          .AND. sax(il,i)>0.0 .AND. iflag(il)<=1) THEN               ! cld
        wa(il, i) = sqrt(2.*sax(il,i))                               ! cld
      END IF                                                         ! cld
    END DO                                                           ! cld
  END DO  
                                                           ! cld
  DO i = 1, nl  

! 14/01/15 AJ je remets les parties manquantes cf JYG
! Initialize sument to 0

    DO il = 1,ncum
     sument(il) = 0.
    ENDDO

! Sum mixed mass fluxes in sument

    DO k = 1,nl
      DO il = 1,ncum
        IF  (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN   ! cld
          sument(il) =sument(il) + abs(ment(il,k,i))
          detrain(il,i) = detrain(il,i) + abs(ment(il,k,i))*(qdet(il,k,i) - rr(il,i))*(qdet(il,k,i) - rr(il,i)) ! Louis terme de dÃ©trainement dans le bilan de variance
        ENDIF
      ENDDO     ! il
    ENDDO       ! k

! 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
!!        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
        sigaq = siga(il,i)*qta(il,i-1)                               ! cld
      ENDIF

! IM cf. FH 
! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB           
                                                         
      IF (iflag_clw==0) THEN                                         ! cld
        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) &       ! cld
          +(1.-siga(il,i))*qcond(il, i)                              ! cld


        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)*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)*qta(il,i-1)+(1.-siga(il,i))*qtment(il,i) ! cld
!     print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i)  
      		
      ELSE IF (iflag_clw==1) THEN                                    ! cld
        qcondc(il, i) = qcond(il, i)                                 ! cld
        qtc(il,i) = qtment(il,i)                                     ! cld
      END IF                                                         ! cld

    END DO                                                           ! cld
  END DO
! print*,'cv3_yield fin'

  RETURN
END SUBROUTINE cv3_yield24
