SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
                       ph, t, rr, rs, u, v, tra, h, lv, lf, frac, qta, &
                       unk, vnk, hp, tv, tvp, ep, clw, sig, &
                       Ment, Qent, hent, uent, vent, nent, &
                       Sigij, elij, supmax, Ments, Qents, traent &
#ifdef ISO
                         ,xt,xtta,xtclw,xtent,xtelij  &
#endif         
                         )
! **************************************************************
! *
! CV3P_MIXING : compute mixed draught properties and,         *
! within a scaling factor, mixed draught        *
! mass fluxes.                                  *
! written by  : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
! modified by :                                               *
! **************************************************************

  USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
  USE ioipsl_getin_p_mod, ONLY: getin_p
  USE add_phys_tend_mod, ONLY: fl_cor_ebil
#ifdef ISO
  USE infotrac_phy, ONLY: ntraciso=>ntiso
  USE isotopes_mod, ONLY: pxtmelt,pxtice
  USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall
#ifdef ISOTRAC
  USE infotrac_phy, ONLY: niso
  USE isotrac_mod, ONLY: option_cond,izone_cond, &
        option_tmin,option_traceurs,seuil_tag_tmin, &
        index_zone,index_iso
  USE isotrac_routines_mod, ONLY: iso_recolorise_condensation
  USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac
#endif
#ifdef ISOVERIF
       USE isotopes_verif_mod
!, ONLY: iso_verif_positif_nostop, iso_verif_egalite_nostop, &
!          errmax,errmaxrel,deltalim
       USE isotopes_mod, ONLY: ridicule,iso_eau,iso_hdo
#endif
#endif
  IMPLICIT NONE

  include "cvthermo.h"
  include "cv3param.h"
  include "YOMCST2.h"
  include "cvflag.h"

!inputs:
  INTEGER, INTENT (IN)                               :: ncum, nd, na
  INTEGER, INTENT (IN)                               :: ntra, nloc
  INTEGER, DIMENSION (nloc), INTENT (IN)             :: icb, inb, nk
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: sig
  REAL, DIMENSION (nloc), INTENT (IN)                :: unk, vnk
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qta
  REAL, DIMENSION (nloc, nd+1), INTENT (IN)          :: ph
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: t, rr, rs
  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: u, v
  REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: tra ! input of convect3
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lv
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: lf
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: frac !ice fraction in condensate
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: h  !liquid water static energy of environment
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: hp !liquid water static energy of air shed from adiab. asc.
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: tv, tvp
  REAL, DIMENSION (nloc, na), INTENT (IN)            :: ep, clw
#ifdef ISO
  REAL, DIMENSION (ntraciso,nloc, na), INTENT (IN)            :: xtclw
  REAL, DIMENSION (ntraciso,nloc, nd), INTENT (IN)            :: xt
  REAL, DIMENSION (ntraciso,nloc, nd), INTENT (IN)            :: xtta
#endif

!outputs:
  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Ment, Qent
  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: uent, vent
  REAL, DIMENSION (nloc, na, na), INTENT (OUT)       :: Sigij, elij
  REAL, DIMENSION (nloc, na), INTENT (OUT)           :: supmax           ! Highest mixing fraction of mixed
                                                                         ! updraughts with the sign of (h-hp)
  REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: Ments, Qents
  REAL, DIMENSION (nloc, nd, nd), INTENT (OUT)       :: hent
  INTEGER, DIMENSION (nloc, nd), INTENT (OUT)        :: nent
#ifdef ISO
  REAL, DIMENSION (ntraciso,nloc, na, na), INTENT (OUT)       :: xtent
  REAL, DIMENSION (ntraciso,nloc, na, na), INTENT (OUT)       :: xtelij
#endif


!local variables:
  INTEGER i, j, k, il, im, jm
  INTEGER num1, num2
  REAL                               :: rti, bf2, anum, denom, dei, altem, cwat, stemp
  REAL                               :: alt, delp, delm
  REAL, DIMENSION (nloc)             :: Qmixmax, Rmixmax, sqmrmax
  REAL, DIMENSION (nloc)             :: Qmixmin, Rmixmin, sqmrmin
  REAL, DIMENSION (nloc)             :: signhpmh
  REAL, DIMENSION (nloc)             :: Sx
  REAL                               :: Scrit2
  REAL, DIMENSION (nloc)             :: Smid, Sjmin, Sjmax
  REAL, DIMENSION (nloc)             :: Sbef, sup, smin
  REAL, DIMENSION (nloc)             :: ASij, ASij_inv, smax, Scrit
  REAL, DIMENSION (nloc, nd, nd)     :: Sij
  REAL, DIMENSION (nloc, nd)         :: csum
  REAL                               :: awat
  REAL                               :: cpm        !Mixed draught heat capacity
  REAL                               :: Tm         !Mixed draught temperature
  LOGICAL, DIMENSION (nloc)          :: lwork

  REAL amxupcrit, df, ff
  INTEGER nstep

  INTEGER,SAVE                                       :: igout=1
!$OMP THREADPRIVATE(igout)

! --   Mixing probability distribution functions

  REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F

  Qcoef1(F) = tanh(F/gammas)
  Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
  QFF(F) = max(min(F,1.), 0.)
  QFFf(F) = min(QFF(F), scut)
  Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
  Rmix1(F) = (gammas*log(cosh((QFF(F)-Fmax)/gammas))+QFF(F)*Qcoef1max)/Qcoef2max
  Qmix2(F) = -log(1.-QFFf(F))/scut
  Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
  Qmix(F) = qqa1*Qmix1(F) + qqa2*Qmix2(F)
  Rmix(F) = qqa1*Rmix1(F) + qqa2*Rmix2(F)

  INTEGER, SAVE :: ifrst
  DATA ifrst/0/
!$OMP THREADPRIVATE(ifrst)

#ifdef ISO
      integer ixt
      real xtrti(ntraciso)
      real xtres(ntraciso)
      real zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc)
      real qent_diag(nloc,nd,nd)
      real xtent_diag(ntraciso,nloc,nd,nd)
#endif   


! =====================================================================
! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
! =====================================================================

! -- Initialize mixing PDF coefficients
  IF (ifrst==0) THEN
    ifrst = 1
    Qcoef1max = Qcoef1(Fmax)
    Qcoef2max = Qcoef2(Fmax)
!<jyg
   PRINT*, 'fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max ', &
            fmax, gammas, qqa1, qqa2, Qcoef1max, Qcoef2max
!>jyg

  END IF


! ori        do 360 i=1,ncum*nlp
  DO j = 1, nl
    DO i = 1, ncum
      nent(i, j) = 0
! in convect3, m is computed in cv3_closure
! ori          m(i,1)=0.0
    END DO
  END DO

#ifdef ISOVERIF
        Qent(:,:,:)=0.0
        xtent(:,:,:,:)=0.0
#endif

! ori      do 400 k=1,nlp
! ori       do 390 j=1,nlp
  DO j = 1, nl
    DO k = 1, nl
      DO i = 1, ncum
        Qent(i, k, j) = rr(i, j)
        uent(i, k, j) = u(i, j)
        vent(i, k, j) = v(i, j) 
        elij(i, k, j) = 0.0
        hent(i, k, j) = 0.0
!AC!            Ment(i,k,j)=0.0
!AC!            Sij(i,k,j)=0.0
#ifdef ISO
        xtelij(:,i, k, j) = 0.0
        xtent(:,i, k, j) = xt(:,i, j)
#endif
      END DO
    END DO
  END DO
#ifdef ISO
        ! on initialise mieux par sécurité:
        elij(:, :, :)=0.0
        xtelij(:,:, :, :)=0.0
#endif

!AC!
  Ment(1:ncum, 1:nd, 1:nd) = 0.0
  Sij(1:ncum, 1:nd, 1:nd) = 0.0
!AC!
!ym
  Sigij(1:ncum, 1:nd, 1:nd) = 0.0
!ym 

!jyg!  DO k = 1, ntra
!jyg!    DO j = 1, nd ! instead nlp
!jyg!      DO i = 1, nd ! instead nlp
!jyg!        DO il = 1, ncum
!jyg!          traent(il, i, j, k) = tra(il, j, k)
!jyg!        END DO
!jyg!      END DO
!jyg!    END DO
!jyg!  END DO

! =====================================================================
! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
! --- FRACTION (Sij)
! =====================================================================

  DO i = minorig + 1, nl

    IF (ok_entrain) THEN
      DO j = minorig, nl
        DO il = 1, ncum
          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
                           .AND. (j<=inb(il))) THEN

!!            rti = qnk(il) - ep(il, i)*clw(il, i)
            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
            bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
!jyg(from aj)<
            IF (cvflag_ice) THEN
! PRINT*,cvflag_ice,'cvflag_ice dans do 700'
              IF (t(il,j)<=263.15) THEN
                bf2 = 1. + (lf(il,j)+lv(il,j))*(lv(il,j)+frac(il,j)* &
                     lf(il,j))*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd)
              END IF
            END IF ! IF (cvflag_ice) THEN
!>jyg
            anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j))
            denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j)
            dei = denom
            IF (abs(dei)<0.01) dei = 0.01
            Sij(il, i, j) = anum/dei
            Sij(il, i, i) = 1.0
            altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
            altem = altem/bf2
            cwat = clw(il, j)*(1.-ep(il,j))
            stemp = Sij(il, i, j)
            IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
!jyg(from aj)<
              IF (cvflag_ice) THEN
                anum = anum - (lv(il,j)+frac(il,j)*lf(il,j))*(rti-rs(il,j)-cwat*bf2)
                denom = denom + (lv(il,j)+frac(il,j)*lf(il,j))*(rr(il,i)-rti)
              ELSE
                anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2)
                denom = denom + lv(il, j)*(rr(il,i)-rti)
              END IF
!>jyg
              IF (abs(denom)<0.01) denom = 0.01
              Sij(il, i, j) = anum/denom
              altem = Sij(il, i, j)*rr(il, i) + (1.-Sij(il,i,j))*rti - rs(il, j)
              altem = altem - (bf2-1.)*cwat
            END IF
            IF (Sij(il,i,j)>0.0) THEN
!!!                 Ment(il,i,j)=m(il,i)
              Ment(il, i, j) = 1.
              elij(il, i, j) = altem
              elij(il, i, j) = amax1(0.0, elij(il,i,j))
              nent(il, i) = nent(il, i) + 1
            END IF

            Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
            Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
          ELSE IF (j > i) THEN
            IF (prt_level >= 10) THEN
              print *,'cv3p_mixing i, j, Sij given by the no-precip eq. ', i, j, Sij(il,i,j)
            ENDIF
          END IF ! new ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
#ifdef ISO
         if(sij(il,i,j).gt.0.0)then
               qent_diag(il,i,j)=sij(il,i,j)*rr(il,i) &
                     +(1.-sij(il,i,j))*rti
           ! question: pourquoi qent n'est-il pas recalculé ici?
         endif !if(sij(il,i,j).gt.0.0)then
#endif
        END DO ! il

#ifdef ISO
       do il=1,ncum
         zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
         zfice(il) = MIN(MAX(zfice(il),0.0),1.0)

         if( (i.ge.icb(il)).and.(i.le.inb(il)).and. &
           (j.ge.(icb(il)-1)).and.(j.le.inb(il))) then
          if (sij(il,i,j).gt.0.0) then  
           do ixt=1,ntraciso
            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i) 
            xtent_diag(ixt,il,i,j)=sij(il,i,j)*xt(ixt,il,i) &
                     +(1.-sij(il,i,j))*xtrti(ixt)   
           enddo ! do ixt=1,ntraciso 
#ifdef ISOVERIF 
          if (iso_eau.gt.0) then
          CALL iso_verif_egalite(xtta(iso_eau,il,i-1), &
                qta(il,i-1),'cv3p_mixing 256a')
          CALL iso_verif_egalite(xt(iso_eau,il,i), &
                rr(il,i),'cv3p_mixing 256b')
          CALL iso_verif_egalite(xtclw(iso_eau,il,i), &
                clw(il,i),'cv3p_mixing 256c')
          CALL iso_verif_egalite(xtent_diag(iso_eau,il,i,j), &
                qent_diag(il,i,j),'cv3p_mixing 256d')    
          endif ! if (iso_eau.gt.0) then      
#ifdef ISOTRAC
          CALL iso_verif_traceur(xtrti,'cv3p_mixing 1968')
#endif           
#endif           
           if (qent_diag(il,i,j).lt.elij(il,i,j)) then  
             do ixt=1,ntraciso
               xtent_diag(ixt,il,i,j)=  xtent_diag(ixt,il,i,j) &
                / qent_diag(il,i,j)*elij(il,i,j)    
             enddo ! do ixt=1,ntraciso
             qent_diag(il,i,j)=elij(il,i,j)
           endif !if (qent_diag(il,i,j).lt.elij(il,i,j)) then

#ifdef ISOVERIFIDRIS
          CALL iso_verif_noNaN(qent_diag(il,i,j),'cv3p_routines 261')
#endif
#ifdef ISOVERIF
          if (iso_eau.gt.0) then
           if (iso_verif_egalite_nostop(xtent_diag(iso_eau,il,i,j), &
                qent_diag(il,i,j),'cv3p_mixing 258').eq.1) then
             write(*,*) 'il,i,j=',il,i,j
             write(*,*) 'qent_diag(il,i,j),elij(il,i,j)=', &
                qent_diag(il,i,j),elij(il,i,j)
             write(*,*) 'sij(il,i,j)=',sij(il,i,j)
             write(*,*) 'xt=',xt(iso_eau,il,i),rr(il,i)
             write(*,*) 'xtclw=',xtclw(iso_eau,il,i),clw(il,i)
             write(*,*) 'ep(il,i)=',ep(il,i)
             stop
           endif
          endif
          CALL iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
                'cv3p_routines 262')
#endif
         else !if (sij(il,i,j).gt.0.0) then 
            ! ajout le 4 mai 2012: initialiser au cas où
            elij(il,i,j)=0.0 
            qent_diag(il,i,j)=rr(il,i)
            do ixt=1,ntraciso
              xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
            enddo !do ixt=1,ntraciso
         endif  !if(sij(il,i,j).gt.0.0)then  
        else  !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
            ! ajout le 4 mai 2012: initialiser au cas où
            elij(il,i,j)=0.0 
            qent_diag(il,i,j)=rr(il,i)
            do ixt=1,ntraciso
              xtent_diag(ixt,il,i,j)=xt(ixt,il,i)
            enddo !do ixt=1,ntraciso
        endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
#ifdef ISOVERIF
#ifdef VERIFNEGATIF
        CALL iso_verif_positif(qent_diag(il,i,j), &
                'cv3p_mixing 303: qt<0')
#endif      
        CALL iso_verif_positif(qent_diag(il,i,j)-elij(il,i,j), &
                'cv3p_mixing 305: qt<elij')
#endif       
       enddo  !do il=1,ncum 

       !write(*,*) 'cv3p_mixing 265: CALL condiso'
       CALL condiso_liq_ice_vectall(xtent_diag(1,1,i,j), &
                qent_diag(1,i,j),elij(1,i,j), &
                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
#ifdef ISOTRAC
        CALL condiso_liq_ice_vectall_trac(xtent_diag(1,1,i,j), &
                qent_diag(1,i,j),elij(1,i,j), &
                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)  
#ifdef ISOVERIF
        do il=1,ncum
          CALL iso_verif_traceur(xt(1,il,i),'cv3p_mixing 1967')
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1969')
        enddo !do il=1,ncum
#endif            
#endif    
        do il=1,ncum
         do ixt = 1, ntraciso
          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
         enddo !do ixt = 1, ntraciso
        enddo !do il=1,ncum

#ifdef ISOTRAC    
!        write(*,*) 'cv3p_mixing tmp 1987,option_traceurs=',
!     :           option_traceurs
        if (option_tmin.ge.1) then
        do il=1,ncum    
        ! colorier la vapeur résiduelle selon température de
        ! condensation, et le condensat en un tag spécifique
          if ((elij(il,i,j).gt.0.0).and. &
                (qent_diag(il,i,j).gt.0.0)) then  
            if (option_traceurs.eq.17) then        
             CALL iso_recolorise_condensation(qent_diag(il,i,j), &
                elij(il,i,j), &
                xtent_diag(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
                0.0,xtres, &
                seuil_tag_tmin)
            else !if (option_traceurs.eq.17) then
!             write(*,*) 'cv3 2002: il,i,j  =',il,i,j   
             CALL iso_recolorise_condensation(qent_diag(il,i,j), &
                elij(il,i,j), &
                xtent_diag(1,il,i,j),xtelij(1,il,i,j),rs(il,j), &
                0.0,xtres,seuil_tag_tmin)
            endif !if (option_traceurs.eq.17) then
            do ixt=1+niso,ntraciso
               xtent_diag(ixt,il,i,j)=xtres(ixt)
            enddo     
          endif !if (cond.gt.0.0) then
        enddo !do il=1,ncum
#ifdef ISOVERIF
        do il=1,ncum 
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_mixing 1996')
          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3p_mixing 1997')
          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
                'cv3p_mixing 2042')
        enddo !do il=1,ncum 
#endif        
        endif !if (option_tmin.ge.1) then       
#endif

! fractionation:
#ifdef ISOVERIF                 
        do il=1,ncum
        if ((i.ge.icb(il)).and.(i.le.inb(il)).and. &
           (j.ge.(icb(il)-1)).and.(j.le.inb(il))) then
        if (sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.scut) then  
        if (iso_eau.gt.0) then
          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,j), &
             qent(il,i,j),'cv3p_mixing 1889',errmax,errmaxrel)   
          CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,j), &
             elij(il,i,j),'cv3p_mixing 1890',errmax,errmaxrel) 
        endif
        if (iso_hdo.gt.0) then    
          CALL iso_verif_aberrant_choix(xt(iso_HDO,il,i),rr(il,i), &
                 ridicule,deltalim,'cv3p_mixing 1997')         
          CALL iso_verif_aberrant_choix( &
                 xtent(iso_HDO,il,i,j),qent(il,i,j), &
                 ridicule,deltalim,'cv3p_mixing 1931')
          CALL iso_verif_aberrant_choix( &
                 xtelij(iso_HDO,il,i,j),elij(il,i,j), &
                 ridicule,deltalim_snow,'cv3p_mixing 1993')
        endif !if (iso_hdo.gt.0) then
#ifdef ISOTRAC  
!        write(*,*) 'cv3p_mixing tmp 2039 il=',il
           CALL iso_verif_traceur(xtent(1,il,i,j), &
                        'cv3p_mixing 2031')
           CALL iso_verif_traceur(xtelij(1,il,i,j), &
                        'cv3p_mixing 2033')
#endif        

        endif !if(sij(il,i,j).gt.0.0)then  
        endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
        enddo !do il=1,ncum
#endif
!        write(*,*) 'cv3p_routine tmp 1984: cond=',elij(il,i,j)                
#endif

      END DO ! j
    ELSE  ! (ok_entrain)
      DO il = 1,ncum
        nent(il,i) = 0
      ENDDO
    ENDIF ! (ok_entrain)

!jygdebug<
    IF (prt_level >= 10) THEN
      print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
      IF (nent(igout,i) > 0) THEN
        print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
      ENDIF
    ENDIF
!>jygdebug

! ***   if no air can entrain at level i assume that updraft detrains  ***
! ***   at that level and calculate detrained air flux and properties  ***


! @      do 170 i=icb(il),inb(il)

    DO il = 1, ncum
      IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
! @      if(nent(il,i).eq.0)then
!!!       Ment(il,i,i)=m(il,i)
        Ment(il, i, i) = 1.
!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
        uent(il, i, i) = unk(il)
        vent(il, i, i) = vnk(il)
        IF (fl_cor_ebil >= 2) THEN
          hent(il, i, i) = hp(il,i)
        ENDIF
        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
        Sij(il, i, i) = 0.0


#ifdef ISO
      do ixt = 1, ntraciso
       xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
       xtelij(ixt,il,i,i)=xtclw(ixt,il,i) ! rq: ne sera pas utilise ensuite
      enddo ! do ixt = 1, ntraciso 
#ifdef ISOVERIF
        if (iso_eau.gt.0) then
          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
             qent(il,i,i),'cv3p_mixing 414',errmax,errmaxrel)
        endif
#endif
#ifdef ISOTRAC          
        if (option_tmin.ge.1) then
        ! colorier la vapeur résiduelle selon température de
        ! condensation, et le condensat en un tag spécifique
!        write(*,*) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=',
!     :            il,i,j,xtent(:,il,i,j)
          if ((elij(il,i,i).gt.0.0).and.(qent(il,i,i).gt.0.0)) then 
            if (option_traceurs.eq.17) then
             CALL iso_recolorise_condensation(qent(il,i,i), &
                elij(il,i,i), &
                xt(1,il,nk(il)),xtclw(1,il,i),t(il,i),ep(il,i), &
                xtres, &
                seuil_tag_tmin)
            else !if (option_traceurs.eq.17) then
             CALL iso_recolorise_condensation(qent(il,i,i), &
                elij(il,i,i), &
                xt(1,il,nk(il)),xtclw(1,il,i),rs(il,i),ep(il,i), &
                xtres, &
                seuil_tag_tmin)
            endif !if (option_traceurs.eq.17) then
            do ixt=1+niso,ntraciso
              xtent(ixt,il,i,i)=xtres(ixt)
            enddo
#ifdef ISOVERIF            
            do ixt=1,niso
            CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
                'cv3p_mixing 2102',errmax,errmaxrel)
            CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
                'cv3p_mixing 2154')
            enddo
#endif            
          endif !if (cond.gt.0.0) then
          
#ifdef ISOVERIF          
          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
                qent(il,i,i),'cv3p_mixing 2103',errmax,errmaxrel)
          CALL iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2095')
          CALL iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2096')
#endif        
        endif !if (option_tmin.ge.1) then    
#endif
! endif ISOTRAC
#endif
! endif ISO

      END IF ! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
    END DO ! il
  END DO ! i = minorig + 1, nl

!jyg!  DO j = 1, ntra
!jyg!    DO i = minorig + 1, nl
!jyg!      DO il = 1, ncum
!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
!jyg!        END IF
!jyg!      END DO
!jyg!    END DO
!jyg!  END DO

  DO j = minorig, nl
    DO i = minorig, nl
      DO il = 1, ncum
        IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
            (i>=icb(il)) .AND. (i<=inb(il))) THEN
          Sigij(il, i, j) = Sij(il, i, j)
        END IF
      END DO
    END DO
  END DO
! @      enddo

! @170   continue

! =====================================================================
! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
! =====================================================================

  CALL zilch(csum, nloc*nd)

  DO il = 1, ncum
    lwork(il) = .FALSE.
  END DO

! ---------------------------------------------------------------
  DO i = minorig + 1, nl      !Loop on origin level "i"
! ---------------------------------------------------------------

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


!JYG1    Find maximum of SIJ for J>I, if any.

    Sx(:) = 0.

    DO il = 1, ncum
      IF (i>=icb(il) .AND. i<=inb(il)) THEN
        signhpmh(il) = sign(1., hp(il,i)-h(il,i))
        Sbef(il) = max(0., signhpmh(il))
      END IF
    END DO

    DO j = i + 1, nl
      DO il = 1, ncum
        IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
          IF (Sbef(il)<Sij(il,i,j)) THEN
            Sx(il) = max(Sij(il,i,j), Sx(il))
          END IF
          Sbef(il) = Sij(il, i, j)
        END IF
      END DO
    END DO


    DO il = 1, ncum
      IF (i>=icb(il) .AND. i<=inb(il)) THEN
        lwork(il) = (nent(il,i)/=0)
!!        rti = qnk(il) - ep(il, i)*clw(il, i)
        rti = qta(il,i-1) - ep(il, i)*clw(il, i)
!jyg<
        IF (cvflag_ice) THEN

          anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)*lf(il,i))* &
                       (rti-rs(il,i)) + (cpv-cpd)*t(il, i)*(rti-rr(il,i))
          denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)*lf(il,i))* &
                       (rr(il,i)-rti) + (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
        ELSE

          anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
                       (cpv-cpd)*t(il, i)*(rti-rr(il,i))
          denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
                       (cpd-cpv)*t(il, i)*(rr(il,i)-rti)
        END IF
!>jyg
        IF (abs(denom)<0.01) denom = 0.01
        Scrit(il) = min(anum/denom, 1.)
        alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)

!JYG1    Find new critical value Scrit2
!         such that : Sij > Scrit2  => mixed draught will detrain at J<I
!                     Sij < Scrit2  => mixed draught will detrain at J>I

        Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
                 Scrit(il)*max(0., signhpmh(il))

        Scrit(il) = Scrit2

!JYG    Correction pour la nouvelle logique; la correction pour ALT
! est un peu au hazard
        IF (Scrit(il)<=0.0) Scrit(il) = 0.0
        IF (alt<=0.0) Scrit(il) = 1.0

        smax(il) = 0.0
        ASij(il) = 0.0
        sup(il) = 0.      ! upper S-value reached by descending draughts
      END IF
    END DO

! ---------------------------------------------------------------
    DO j = minorig, nl         !Loop on destination level "j"
! ---------------------------------------------------------------

      num2 = 0
      DO il = 1, ncum
        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
            lwork(il)) num2 = num2 + 1
      END DO
      IF (num2<=0) GO TO 175

! -----------------------------------------------
      IF (j>i) THEN
! -----------------------------------------------
        DO il = 1, ncum
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
              lwork(il)) THEN
            IF (Sij(il,i,j)>0.0) THEN
              Smid(il) = min(Sij(il,i,j), Scrit(il))
              Sjmax(il) = Smid(il)
              Sjmin(il) = Smid(il)
              IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
                smin(il) = Smid(il)
                Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
                Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
                Sjmin(il) = min(Sjmin(il), Scrit(il))
                Sbef(il) = Sij(il, i, j)
              END IF
            END IF
          END IF
        END DO
! -----------------------------------------------
      ELSE IF (j==i) THEN
! -----------------------------------------------
        DO il = 1, ncum
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
              lwork(il)) THEN
            IF (Sij(il,i,j)>0.0) THEN
              Smid(il) = 1.
              Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
                          min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
              Sjmin(il) = max(Sjmin(il), sup(il))
              Sjmax(il) = 1.

! -             preparation des variables Scrit, Smin et Sbef pour la partie j>i
              Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))

              smin(il) = 1.
              Sbef(il) = max(0., signhpmh(il))
              supmax(il, i) = sign(Scrit(il), -signhpmh(il))
            END IF
          END IF
        END DO
! -----------------------------------------------
      ELSE IF (j<i) THEN
! -----------------------------------------------
        DO il = 1, ncum
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
              lwork(il)) THEN
            IF (Sij(il,i,j)>0.0) THEN
              Smid(il) = max(Sij(il,i,j), Scrit(il))
              Sjmax(il) = Smid(il)
              Sjmin(il) = Smid(il)
              IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
                smax(il) = Smid(il)
                Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
                Sjmax(il) = max(Sjmax(il), Scrit(il))
                Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
                Sjmin(il) = max(Sjmin(il), Scrit(il))
                Sbef(il) = Sij(il, i, j)
              END IF
              IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
                             sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
            END IF
          END IF
        END DO
! -----------------------------------------------
      END IF
! -----------------------------------------------


      DO il = 1, ncum
        IF (i>=icb(il) .AND. i<=inb(il) .AND. &
            j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
            lwork(il)) THEN
          IF (Sij(il,i,j)>0.0) THEN
!!            rti = qnk(il) - ep(il, i)*clw(il, i)
            rti = qta(il,i-1) - ep(il, i)*clw(il, i)
            Qmixmax(il) = Qmix(Sjmax(il))
            Qmixmin(il) = Qmix(Sjmin(il))
            Rmixmax(il) = Rmix(Sjmax(il))
            Rmixmin(il) = Rmix(Sjmin(il))
            sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
            sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))

            Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)

! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
            IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
              Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
            ELSE
              Sigij(il, i, j) = 0.
            END IF

! --    Compute Qent, uent, vent according to the true mixing fraction
            Qent(il, i, j) = (1.-Sigij(il,i,j))*rti     + Sigij(il, i, j)*rr(il, i)
            uent(il, i, j) = (1.-Sigij(il,i,j))*unk(il) + Sigij(il, i, j)*u(il, i)
            vent(il, i, j) = (1.-Sigij(il,i,j))*vnk(il) + Sigij(il, i, j)*v(il, i)

! --     Compute liquid water static energy of mixed draughts
!    IF (j .GT. i) THEN
!      awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
!      awat=amax1(awat,0.0)
!    ELSE
!      awat = 0.
!    ENDIF
!    Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
!    :         + Sigij(il,i,j)*H(il,i)
!    :         + (LV(il,j)+(cpd-cpv)*t(il,j))*awat
!IM 301008 beg
            hent(il, i, j) = (1.-Sigij(il,i,j))*hp(il, i) + Sigij(il, i, j)*h(il, i)

!jyg<
!            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
!            elij(il, i, j) = elij(il, i, j) + &
!                             ((h(il,j)-hent(il,i,j))*rs(il,j)*lv(il,j) / &
!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))
!            elij(il, i, j) = elij(il, i, j) / &
!                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
!                              ((cpd*(1.-Qent(il,i,j))+Qent(il,i,j)*cpv)*rrv*t(il,j)*t(il,j)))

!       Computation of condensate amount Elij, taking into account the ice fraction frac
!       Warning : the same saturation humidity rs is used over both liquid water and ice; this
!                 should be corrected.

!  Heat capacity of mixed draught
    cpm = cpd+Qent(il,i,j)*(cpv-cpd)

    IF (cvflag_ice .and. frac(il,j) > 0.) THEN
            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
            elij(il, i, j) = elij(il, i, j) + &
                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
            elij(il, i, j) = elij(il, i, j) / &
                             (1.+(lv(il,j)+frac(il,j)*lf(il,j))*lv(il,j)*rs(il,j) / &
                              (cpm*rrv*t(il,j)*t(il,j)))
    ELSE
            elij(il, i, j) = Qent(il, i, j) - rs(il, j)
            elij(il, i, j) = elij(il, i, j) + &
                             (h(il,j)-hent(il,i,j)+(cpv-cpd)*(Qent(il,i,j)-rr(il,j))*t(il,j))* &
                             rs(il,j)*lv(il,j) / (cpm*rrv*t(il,j)*t(il,j))
            elij(il, i, j) = elij(il, i, j) / &
                             (1.+lv(il,j)*lv(il,j)*rs(il,j) / &
                              (cpm*rrv*t(il,j)*t(il,j)))
    ENDIF
!>jyg
            elij(il, i, j) = max(elij(il,i,j), 0.)

            elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))

            IF (j>i) THEN
              awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
              awat = amax1(awat, 0.0)
            ELSE
              awat = 0.
            END IF

! print *,h(il,j)-hent(il,i,j),LV(il,j)*rs(il,j)/(cpd*rrv*t(il,j)*
! :         t(il,j))

!jyg<
!            hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*t(il,j))*awat
! Mixed draught temperature at level j
    IF (cvflag_ice .and. frac(il,j) > 0.) THEN
          Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))*rrv*t(il,j)*t(il,j)/(lv(il,j)*rs(il,j))
          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+frac(il,j)*lf(il,j)+(cpd-cpv)*Tm)*awat
    ELSE
          Tm = t(il,j) + (Qent(il,i,j)-elij(il,i,j)-rs(il,j))*rrv*t(il,j)*t(il,j)/(lv(il,j)*rs(il,j))
          hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)*Tm)*awat
    ENDIF
!>jyg

!IM 301008 end

! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
! :               i,j,hent(il,i,j),Sigij(il,i,j)

! --      ASij is the integral of P(F) over the relevant F interval
            ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
                                      Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))

          END IF
        END IF
      END DO ! il


              
#ifdef ISO
      do il=1,ncum
          zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
          zfice(il) = MIN(MAX(zfice(il),0.0),1.0)  
       if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
           j.ge.(icb(il)-1) .and. j.le.inb(il) &
           .and. lwork(il) ) then
        if(sij(il,i,j).gt.0.0)then

          do ixt=1,ntraciso           
            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)         
            xtent(ixt,il,i,j) = (1.-Sigij(il,i,j))*xtrti(ixt) &
                    + Sigij(il,i,j)*xt(ixt,il,i)
          enddo !do ixt = 1, ntraciso      
#ifdef ISOVERIF 
#ifdef ISOTRAC          
          CALL iso_verif_traceur(xtrti,'cv3p_routines 714')
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
#endif          
          if (iso_eau.gt.0) then
            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
                qent(il,i,j),'cv3p_routines 718')    
          endif
#endif              
         endif  !if(sij(il,i,j).gt.0.0)then  
        endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
       enddo  !do il=1,ncum
   
#ifdef ISOVERIFIDRIS
        do il=1,ncum       
         CALL iso_verif_noNaN(qent(il,i,j),'cv3p_routines 771')
        enddo !do il=1,ncum   
#endif     
#ifdef ISOVERIF
        if (iso_eau.gt.0) then
         do il=1,ncum   
            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
                qent(il,i,j),'cv3p_routines 744')   
         enddo !do il=1,ncum   
        endif    
        do il=1,ncum       
         CALL iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
                'cv3p_routines 749')
        enddo !do il=1,ncum   
        !write(*,*) 'cv3p_mixing 788: CALL condiso_liq_ice_vectall'
#endif     
      
       CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
                elij(1,i,j), &
                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
#ifdef ISOTRAC
        CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
                elij(1,i,j), &
                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)
#endif    
        do il=1,ncum
        if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
           j.ge.(icb(il)-1) .and. j.le.inb(il) &
           .and. lwork(il) ) then
        if(sij(il,i,j).gt.0.0)then
         do ixt = 1, ntraciso
          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
         enddo !do ixt = 1, ntraciso
        endif !if(sij(il,i,j).gt.0.0)then
        endif !if ( i.ge.icb(il) .and. i.le.inb(il) .and.
        enddo !do il=1,ncum

#ifdef ISOVERIF  
        if (iso_eau.gt.0) then
            do il=1,ncum
            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
                qent(il,i,j),'cv3p_routines 744') 
            CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
                elij(il,i,j),'cv3p_routines 746') 
            enddo !do il=1,ncum   
        endif
#endif  
!#ifdef ISOVERIF  

#ifdef ISOTRAC    
!        write(*,*) 'cv3_routines tmp 1987,option_traceurs=',
!     :           option_traceurs
        if (option_tmin.ge.1) then
        do il=1,ncum    
!        write(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
!     :           'tcond(il),rs(il,j)=',
!     :            il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
        ! colorier la vapeur résiduelle selon température de
        ! condensation, et le condensat en un tag spécifique
          if ((elij(il,i,j).gt.0.0).and.(qent(il,i,j).gt.0.0)) then  
            if (option_traceurs.eq.17) then        
             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
                xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
                0.0,xtres, &
                seuil_tag_tmin)
            else !if (option_traceurs.eq.17) then
!             write(*,*) 'cv3 2002: il,i,j  =',il,i,j   
             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
                xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
                seuil_tag_tmin)
            endif !if (option_traceurs.eq.17) then
            do ixt=1+niso,ntraciso
               xtent(ixt,il,i,j)=xtres(ixt)
            enddo     
          endif !if (cond.gt.0.0) then
        enddo !do il=1,ncum
#ifdef ISOVERIF
        do il=1,ncum 
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
                'cv3_routines 2042')
        enddo !do il=1,ncum 
#endif        
        endif !if (option_tmin.ge.1) then       
#endif
! endif ISOTRAC

#endif
!!#ifdef ISO

!jyg!      DO k = 1, ntra
!jyg!        DO il = 1, ncum
!jyg!          IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
!jyg!              (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
!jyg!              lwork(il)) THEN
!jyg!            IF (Sij(il,i,j)>0.0) THEN
!jyg!              traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
!jyg!                                    (1.-Sigij(il,i,j))*tra(il, nk(il), k)
!jyg!            END IF
!jyg!          END IF
!jyg!        END DO
!jyg!      END DO

! --    If I=J (detrainement and entrainement at the same level), then only the
! --    adiabatic ascent part of the mixture is considered
      IF (i==j) THEN
        DO il = 1, ncum
          IF (i>=icb(il) .AND. i<=inb(il) .AND. &
              j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
              lwork(il)) THEN
            IF (Sij(il,i,j)>0.0) THEN
!!              rti = qnk(il) - ep(il, i)*clw(il, i)
              rti = qta(il,i-1) - ep(il, i)*clw(il, i)
!!!             Ment(il,i,i) = m(il,i)*abs(Qmixmax(il)*(1.-Sjmax(il))
              Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
                                   Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
              Qent(il, i, i) = rti
              uent(il, i, i) = unk(il)
              vent(il, i, i) = vnk(il)
              hent(il, i, i) = hp(il, i)
              elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
              Sigij(il, i, i) = 0.
            END IF
          END IF
        END DO


#ifdef ISO
      IF (I .EQ. J) THEN
      do il=1,ncum
          zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice)
          zfice(il) = MIN(MAX(zfice(il),0.0),1.0)  
       if ( i.ge.icb(il) .and. i.le.inb(il) .and. &
           j.ge.(icb(il)-1) .and. j.le.inb(il) &
           .and. lwork(il) ) then
        if(sij(il,i,j).gt.0.0)then
          do ixt=1,ntraciso           
            xtrti(ixt)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)         
            xtent(ixt,il,i,j) = xtrti(ixt)
          enddo !do ixt = 1, ntraciso
#ifdef ISOVERIFIDRIS
          CALL iso_verif_noNaN(qent(il,i,j),'cv3p_routines 923')
#endif 
#ifdef ISOVERIF  
#ifdef ISOTRAC          
          CALL iso_verif_traceur(xtrti,'cv3p_routines 714')
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 714')
#endif          
          if (iso_eau.gt.0) then
            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
                qent(il,i,j),'cv3p_routines 718')    
          endif          
          CALL iso_verif_positif(qent(il,i,j)-elij(il,i,j), &
                'cv3p_routines 887')
#endif
         endif  !if(sij(il,i,j).gt.0.0)then  
        endif !if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
       enddo  !do il=1,ncum
      
       !write(*,*) 'cv3p_mixing 808: CALL condiso'
       CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), &
                elij(1,i,j), &
                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum)

       do il=1,ncum
         do ixt = 1, ntraciso
          xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il)
         enddo !do ixt = 1, ntraciso
#ifdef ISOVERIF
         if (iso_eau.gt.0) then
           CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
                elij(il,i,j),'cv3p_routines 1074')
         endif ! if (iso_eau.gt.0) then
#endif
        enddo !do il=1,ncum
        
#ifdef ISOTRAC
        CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), &
                elij(1,i,j), &
                t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum) 
#ifdef ISOVERIF
        do il=1,ncum
          CALL iso_verif_traceur(xt(1,il,i),'cv3p_routines 1967')
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3p_routines 1969')
        enddo !do il=1,ncum
#endif   
        
       if (option_tmin.ge.1) then
        do il=1,ncum    
!        write(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),',
!     :           'tcond(il),rs(il,j)=',
!     :            il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j)
        ! colorier la vapeur résiduelle selon température de
        ! condensation, et le condensat en un tag spécifique
          if ((elij(il,i,j).gt.0.0).and.(qent(il,i,j).gt.0.0)) then  
            if (option_traceurs.eq.17) then        
             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
                xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), &
                0.0,xtres, &
                seuil_tag_tmin)
            else !if (option_traceurs.eq.17) then
!             write(*,*) 'cv3 2002: il,i,j  =',il,i,j   
             CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), &
                xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, &
                seuil_tag_tmin)
            endif !if (option_traceurs.eq.17) then
            do ixt=1+niso,ntraciso
               xtent(ixt,il,i,j)=xtres(ixt)
            enddo     
          endif !if (cond.gt.0.0) then
        enddo !do il=1,ncum
#ifdef ISOVERIF
        do il=1,ncum 
          CALL iso_verif_traceur(xtent(1,il,i,j),'cv3_routines 1996')
          CALL iso_verif_traceur(xtelij(1,il,i,j),'cv3_routines 1997')
          CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
                'cv3_routines 2042')
        enddo !do il=1,ncum 
#endif        
        endif !if (option_tmin.ge.1) then
!#ifdef ISOVERIF        
#endif
!#ifdef ISOTRAC        
#ifdef ISOVERIF        
        if (iso_eau.gt.0) then
            do il=1,ncum
            CALL iso_verif_egalite(xtent(iso_eau,il,i,j), &
                qent(il,i,j),'cv3p_routines 883') 
            CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
                elij(il,i,j),'cv3p_routines 885')  
            enddo !do il=1,ncum  
        endif !if (iso_eau.gt.0) then
#endif  
!#ifdef ISOVERIF           
        endif !IF (I .EQ. J) THEN
#endif  
!#ifdef ISO 

!jyg!        DO k = 1, ntra
!jyg!          DO il = 1, ncum
!jyg!            IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
!jyg!                (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
!jyg!                lwork(il)) THEN
!jyg!              IF (Sij(il,i,j)>0.0) THEN
!jyg!                traent(il, i, i, k) = tra(il, nk(il), k)
!jyg!              END IF
!jyg!            END IF
!jyg!          END DO
!jyg!        END DO

      END IF

! ---------------------------------------------------------------
175 END DO        ! End loop on destination level "j"
! ---------------------------------------------------------------

    DO il = 1, ncum
      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
        ASij(il) = amax1(1.0E-16, ASij(il))
!jyg+lluis<
!!        ASij(il) = 1.0/ASij(il)
        ASij_inv(il) = 1.0/ASij(il)
!   IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
        IF (ASij_inv(il) > 100.)  ASij_inv(il) = 0.
!>jyg+lluis
        csum(il, i) = 0.0
      END IF
    END DO

    DO j = minorig, nl
      DO il = 1, ncum
        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
!jyg          Ment(il, i, j) = Ment(il, i, j)*ASij(il)
          Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
        END IF
      END DO
    END DO

    DO j = minorig, nl
      DO il = 1, ncum
        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
            j>=(icb(il)-1) .AND. j<=inb(il)) THEN
          csum(il, i) = csum(il, i) + Ment(il, i, j)
        END IF
      END DO
    END DO

    DO il = 1, ncum
      IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
! cc     :     .and. csum(il,i).lt.m(il,i) ) then
        nent(il, i) = 0
! cc        Ment(il,i,i)=m(il,i)
        Ment(il, i, i) = 1.
!!        Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
        Qent(il, i, i) = qta(il,i-1) - ep(il, i)*clw(il, i)
        uent(il, i, i) = unk(il)
        vent(il, i, i) = vnk(il)
        elij(il, i, i) = clw(il, i)*(1.-ep(il,i))

#ifdef ISO
      do ixt = 1, ntraciso
        xtent(ixt,il,i,i)=xtta(ixt,il,i-1)-ep(il,i)*xtclw(ixt,il,i)
        xtelij(ixt,il,i,i)=xtclw(ixt,il,i)*(1.-ep(il,i))
      enddo
#ifdef ISOVERIF
      if (iso_eau.gt.0) then
        CALL iso_verif_egalite(xtelij(iso_eau,il,i,j), &
                elij(il,i,j),'cv3p_routines 1204')
      endif ! if (iso_eau.gt.0) then
#endif
#endif

#ifdef ISOTRAC          
        if (option_tmin.ge.1) then
        ! colorier la vapeur résiduelle selon température de
        ! condensation, et le condensat en un tag spécifique
!        write(*,*) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=',
!     :            il,i,j,xtent(:,il,i,j)
          if ((elij(il,i,i).gt.0.0).and.(qent(il,i,i).gt.0.0)) then 
            if (option_traceurs.eq.17) then          
              CALL iso_recolorise_condensation(qent(il,i,i), &
                elij(il,i,i), &
                xt(1,il,1),xtclw(1,il,i),t(il,i),ep(il,i), &
                xtres, &
                seuil_tag_tmin)
            else !if (option_traceurs.eq.17) then
              CALL iso_recolorise_condensation(qent(il,i,i), &
                elij(il,i,i), &
                xt(1,il,1),xtclw(1,il,i),rs(il,i),ep(il,i), &
                xtres, &
                seuil_tag_tmin)
            endif ! if (option_traceurs.eq.17) then
            do ixt=1+niso,ntraciso
              xtent(ixt,il,i,i)=xtres(ixt)
            enddo  
#ifdef ISOVERIF               
            do ixt=1,niso
              CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), &
                'cv3p_mixing 2318',errmax,errmaxrel)
              CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), &
                'cv3p_mixing 2383')
            enddo
#endif               
          endif !if (cond.gt.0.0) then
#ifdef ISOVERIF          
          CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), &
                qent(il,i,i),'cv3p_mixing 2321',errmax,errmaxrel)
          CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,i), &
                elij(il,i,i),'cv3p_mixing 2321b',errmax,errmaxrel)
          CALL iso_verif_traceur(xtent(1,il,i,i),'cv3p_mixing 2322')
          CALL iso_verif_traceur(xtelij(1,il,i,i),'cv3p_mixing 2323')
#endif        
        endif !if (option_tmin.ge.1) then
#endif
        IF (fl_cor_ebil .GE. 2) THEN
          hent(il, i, i) = hp(il,i)
          Sigij(il, i, i) = 0.0
        ELSE
          Sij(il, i, i) = 0.0
        ENDIF
      END IF
    END DO ! il

!jyg!    DO j = 1, ntra
!jyg!      DO il = 1, ncum
!jyg!        IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
!jyg!! cc     :     .and. csum(il,i).lt.m(il,i) ) then
!jyg!          traent(il, i, i, j) = tra(il, nk(il), j)
!jyg!        END IF
!jyg!      END DO
!jyg!    END DO

! ---------------------------------------------------------------
789 END DO              ! End loop on origin level "i"
! ---------------------------------------------------------------

#ifdef ISO
#ifdef ISOVERIF
!       write(*,*) 'cv3p_mixing 2540: ', &
!        'verif finale en sortant de cv3p_mixing'
!       write(*,*) 'qent,xtent(1,1,1)=',qent(1,1,1),xtent(iso_eau,1,1,1)
       do im = 1, nd
       do jm = 1, nd
        do il = 1, ncum
          if (iso_eau.gt.0) then
            CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
              elij(il,im,jm),'cv3p_mixing 2110',errmax,errmaxrel)
            CALL iso_verif_egalite_choix(xtent(iso_eau,il,im,jm), &
              qent(il,im,jm),'cv3p_mixing 2112',errmax,errmaxrel)
          endif !if (iso_eau>0) then
#ifdef ISOTRAC
        CALL iso_verif_traceur_justmass(xtelij(1,il,im,jm), &
                       'cv3p_routine 2250')
#endif            
        enddo !do il = 1, nloc
       enddo !do jm = 1, klev
       enddo !do im = 1, klev
#endif
#endif  

#ifdef ISO
#ifdef ISOTRAC
        ! seulement à la fin on taggue le condensat
        if (option_cond.ge.1) then
         do im = 1, nd
         do jm = 1, nd
         do il = 1, ncum    
           ! colorier le condensat en un tag spécifique
           do ixt=niso+1,ntraciso
             if (index_zone(ixt).eq.izone_cond) then
                xtelij(ixt,il,im,jm)=xtelij(index_iso(ixt),il,im,jm)
             else !if (index_zone(ixt).eq.izone_cond) then
                xtelij(ixt,il,im,jm)=0.0
             endif !if (index_zone(ixt).eq.izone_cond) then
           enddo !do ixt=1,ntraciso      
#ifdef ISOVERIF
        CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), &
                elij(il,im,jm),'cv3p_mixing 2408',errmax,errmaxrel)
        CALL iso_verif_traceur(xtelij(1,il,im,jm), &
               'cv3p_mixing 358')
#endif     
         enddo !do il = 1, ncum   
         enddo !do jm = 1, nd
         enddo !do im = 1, nd
        ! 17 juin 2020: on supprime ci-dessous car normalement xtclw a déjà été
        ! recolorié ds cv3_undilute
!         do im = 1, nd
!         do il = 1, ncum    
!           ! colorier le condensat en un tag spécifique
!           do ixt=niso+1,ntraciso
!             if (index_zone(ixt).eq.izone_cond) then
!                xtclw(ixt,il,im)=xtclw(index_iso(ixt),il,im)
!             else !if (index_zone(ixt).eq.izone_cond) then
!                xtclw(ixt,il,im)=0.0
!             endif !if (index_zone(ixt).eq.izone_cond) then
!           enddo !do ixt=1,ntraciso      
!#ifdef ISOVERIF
!        CALL iso_verif_egalite_choix(xtclw(iso_eau,il,im), &
!                clw(il,im),'cv3p_mixing 2427',errmax,errmaxrel)
!        CALL iso_verif_traceur(xtclw(1,il,im), &
!               'cv3p_mixing 358')
!        if (iso_verif_positif_nostop(xtclw(itZonIso( &
!                izone_cond,iso_eau),i,k)-xtclw(iso_eau,i,k) &
!                ,'cv3p_mixing 909').eq.1) then
!               write(*,*) 'i,k=',i,k
!               write(*,*) 'xtclw=',xtclw(:,i,k)
!               write(*,*) 'niso,ntraciso,index_zone,izone_cond=', &
!                  niso,ntraciso,index_zone,izone_cond      
!               stop
!         endif !if (iso_verif_positif_nostop(xtclw(itZonIso(
!#endif             
!         enddo !do il = 1, ncum   
!         enddo !do im = 1, nd
!         write(*,*) 'xtclw(:,1,2)=',xtclw(:,1,2)
        endif !if (option_tmin.eq.1) then
#endif
#endif         


  RETURN
END SUBROUTINE cv3p_mixing

