phylmd/CV30_routines/cv30_mixing.f

module cv3_mixing_m

  implicit none

contains

  SUBROUTINE cv3_mixing(nloc, ncum, nd, na, icb, nk, inb, t, rr, rs, u, v, h, &
       lv, hp, ep, clw, m, sig, ment, qent, uent, vent, nent, sij, elij, &
       ments, qents)
    use cv3_param_m
    use cvthermo

    !---------------------------------------------------------------------
    ! a faire:
    ! - changer rr(il, 1) -> qnk(il)
    ! - vectorisation de la partie normalisation des flux (do 789)
    !---------------------------------------------------------------------

    ! inputs:
    integer, intent(in):: ncum, nd, na, nloc
    integer icb(nloc), inb(nloc), nk(nloc)
    real sig(nloc, nd)
    real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
    real u(nloc, nd), v(nloc, nd)
    real lv(nloc, na), h(nloc, na), hp(nloc, na)
    real ep(nloc, na), clw(nloc, na)
    real m(nloc, na) ! input of convect3

    ! outputs:
    real ment(nloc, na, na), qent(nloc, na, na)
    real uent(nloc, na, na), vent(nloc, na, na)
    real sij(nloc, na, na), elij(nloc, na, na)
    real ments(nloc, nd, nd), qents(nloc, nd, nd)
    real sigij(nloc, nd, nd)
    integer nent(nloc, nd)

    ! local variables:
    integer i, j, k, il, im, jm
    integer num1, num2
    real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
    real alt, smid, sjmin, sjmax, delp, delm
    real asij(nloc), smax(nloc), scrit(nloc)
    real asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
    real wgh
    real zm(nloc, na)
    logical lwork(nloc)

    ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS

    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

    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
             !ym ment(i, k, j)=0.0
             !ym sij(i, k, j)=0.0
          end do
       end do
    end do

    !ym
    ment(1:ncum, 1:nd, 1:nd)=0.0
    sij(1:ncum, 1:nd, 1:nd)=0.0

    zm(:, :)=0.

    ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
    ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
    ! --- FRACTION (sij)

    do i=minorig+1, nl

       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=rr(il, 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)
                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
                   anum=anum-lv(il, j)*(rti-rs(il, j)-cwat*bf2)
                   denom=denom+lv(il, j)*(rr(il, i)-rti)
                   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.and.sij(il, i, j) < 0.95)then
                   qent(il, i, j)=sij(il, i, j)*rr(il, i)+(1.-sij(il, i, j))*rti
                   uent(il, i, j)=sij(il, i, j)*u(il, i)+(1.-sij(il, i, j))*u(il, nk(il))
                   vent(il, i, j)=sij(il, i, j)*v(il, i)+(1.-sij(il, i, j))*v(il, nk(il))
                   elij(il, i, j)=altem
                   elij(il, i, j)=amax1(0.0, elij(il, i, j))
                   ment(il, i, j)=m(il, i)/(1.-sij(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))
             endif ! new
          end do
       end do

       ! *** 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) == 0)then
             ment(il, i, i)=m(il, i)
             qent(il, i, i)=rr(il, nk(il))-ep(il, i)*clw(il, i)
             uent(il, i, i)=u(il, nk(il))
             vent(il, i, i)=v(il, nk(il))
             elij(il, i, i)=clw(il, i)
             !MAF sij(il, i, i)=1.0
             sij(il, i, i)=0.0
          end if
       end do
    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)
             endif
          end do
       end do
    end do
    !@ enddo

    !@170 continue

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

    call zilch(asum, nloc*nd)
    call zilch(csum, nloc*nd)
    call zilch(csum, nloc*nd)

    do il=1, ncum
       lwork(il) = .FALSE.
    enddo

    DO i=minorig+1, nl

       num1=0
       do il=1, ncum
          if (i >= icb(il) .and. i <= inb(il)) num1=num1+1
       enddo
       if (num1 <= 0) cycle

       do il=1, ncum
          if (i >= icb(il) .and. i <= inb(il)) then
             lwork(il)=(nent(il, i) /= 0)
             qp=rr(il, 1)-ep(il, i)*clw(il, i)
             anum=h(il, i)-hp(il, i)-lv(il, i)*(qp-rs(il, i)) &
                  +(cpv-cpd)*t(il, i)*(qp-rr(il, i))
             denom=h(il, i)-hp(il, i)+lv(il, i)*(rr(il, i)-qp) &
                  +(cpd-cpv)*t(il, i)*(rr(il, i)-qp)
             if(abs(denom) < 0.01)denom=0.01
             scrit(il)=anum/denom
             alt=qp-rs(il, i)+scrit(il)*(rr(il, i)-qp)
             if(scrit(il) <= 0.0.or.alt <= 0.0)scrit(il)=1.0
             smax(il)=0.0
             asij(il)=0.0
          endif
       end do

       do j=nl, minorig, -1

          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
          enddo
          if (num2 <= 0) cycle

          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) > 1.0e-16.and.sij(il, i, j) < 0.95)then
                   wgh=1.0
                   if(j > i)then
                      sjmax=amax1(sij(il, i, j+1), smax(il))
                      sjmax=amin1(sjmax, scrit(il))
                      smax(il)=amax1(sij(il, i, j), smax(il))
                      sjmin=amax1(sij(il, i, j-1), smax(il))
                      sjmin=amin1(sjmin, scrit(il))
                      if(sij(il, i, j) < (smax(il)-1.0e-16))wgh=0.0
                      smid=amin1(sij(il, i, j), scrit(il))
                   else
                      sjmax=amax1(sij(il, i, j+1), scrit(il))
                      smid=amax1(sij(il, i, j), scrit(il))
                      sjmin=0.0
                      if(j > 1)sjmin=sij(il, i, j-1)
                      sjmin=amax1(sjmin, scrit(il))
                   endif
                   delp=abs(sjmax-smid)
                   delm=abs(sjmin-smid)
                   asij(il)=asij(il)+wgh*(delp+delm)
                   ment(il, i, j)=ment(il, i, j)*(delp+delm)*wgh
                endif
             endif
          end do

       end do

       do il=1, ncum
          if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
             asij(il)=amax1(1.0e-16, asij(il))
             asij(il)=1.0/asij(il)
             asum(il, i)=0.0
             bsum(il, i)=0.0
             csum(il, i)=0.0
          endif
       enddo

       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
                ment(il, i, j)=ment(il, i, j)*asij(il)
             endif
          enddo
       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
                asum(il, i)=asum(il, i)+ment(il, i, j)
                ment(il, i, j)=ment(il, i, j)*sig(il, j)
                bsum(il, i)=bsum(il, i)+ment(il, i, j)
             endif
          enddo
       end do

       do il=1, ncum
          if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
             bsum(il, i)=amax1(bsum(il, i), 1.0e-16)
             bsum(il, i)=1.0/bsum(il, i)
          endif
       enddo

       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
                ment(il, i, j)=ment(il, i, j)*asum(il, i)*bsum(il, i)
             endif
          enddo
       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)
             endif
          enddo
       end do

       do il=1, ncum
          if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
               .and. csum(il, i) < m(il, i)) then
             nent(il, i)=0
             ment(il, i, i)=m(il, i)
             qent(il, i, i)=rr(il, 1)-ep(il, i)*clw(il, i)
             uent(il, i, i)=u(il, nk(il))
             vent(il, i, i)=v(il, nk(il))
             elij(il, i, i)=clw(il, i)
             !MAF sij(il, i, i)=1.0
             sij(il, i, i)=0.0
          endif
       enddo ! il

    end DO

    ! MAF: renormalisation de MENT
    do jm=1, nd
       do im=1, nd
          do il=1, ncum
             zm(il, im)=zm(il, im)+(1.-sij(il, im, jm))*ment(il, im, jm)
          end do
       end do
    end do

    do jm=1, nd
       do im=1, nd
          do il=1, ncum
             if(zm(il, im) /= 0.) then
                ment(il, im, jm)=ment(il, im, jm)*m(il, im)/zm(il, im)
             endif
          end do
       end do
    end do

    do jm=1, nd
       do im=1, nd
          do il=1, ncum
             qents(il, im, jm)=qent(il, im, jm)
             ments(il, im, jm)=ment(il, im, jm)
          end do
       enddo
    enddo

  end SUBROUTINE cv3_mixing

end module cv3_mixing_m
1	guez	97	module cv3_mixing_m
2	guez	47
3	guez	97	implicit none
4	guez	47
5	guez	97	contains
6	guez	47
7	guez	145	SUBROUTINE cv3_mixing(nloc, ncum, nd, na, icb, nk, inb, t, rr, rs, u, v, h, &
8			lv, hp, ep, clw, m, sig, ment, qent, uent, vent, nent, sij, elij, &
9			ments, qents)
10	guez	97	use cv3_param_m
11			use cvthermo
12	guez	47
13	guez	97	!---------------------------------------------------------------------
14			! a faire:
15	guez	145	! - changer rr(il, 1) -> qnk(il)
16			! - vectorisation de la partie normalisation des flux (do 789)
17	guez	97	!---------------------------------------------------------------------
18	guez	47
19	guez	97	! inputs:
20			integer, intent(in):: ncum, nd, na, nloc
21			integer icb(nloc), inb(nloc), nk(nloc)
22	guez	145	real sig(nloc, nd)
23			real t(nloc, nd), rr(nloc, nd), rs(nloc, nd)
24			real u(nloc, nd), v(nloc, nd)
25			real lv(nloc, na), h(nloc, na), hp(nloc, na)
26			real ep(nloc, na), clw(nloc, na)
27			real m(nloc, na) ! input of convect3
28	guez	47
29	guez	97	! outputs:
30	guez	145	real ment(nloc, na, na), qent(nloc, na, na)
31			real uent(nloc, na, na), vent(nloc, na, na)
32			real sij(nloc, na, na), elij(nloc, na, na)
33			real ments(nloc, nd, nd), qents(nloc, nd, nd)
34			real sigij(nloc, nd, nd)
35			integer nent(nloc, nd)
36	guez	47
37	guez	97	! local variables:
38			integer i, j, k, il, im, jm
39			integer num1, num2
40			real rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
41			real alt, smid, sjmin, sjmax, delp, delm
42			real asij(nloc), smax(nloc), scrit(nloc)
43	guez	145	real asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd)
44	guez	97	real wgh
45	guez	145	real zm(nloc, na)
46	guez	97	logical lwork(nloc)
47
48			! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
49
50	guez	145	do j=1, nl
51			do i=1, ncum
52			nent(i, j)=0
53	guez	97	! in convect3, m is computed in cv3_closure
54	guez	145	! ori m(i, 1)=0.0
55	guez	97	end do
56			end do
57	guez	47
58	guez	145	do j=1, nl
59			do k=1, nl
60			do i=1, ncum
61			qent(i, k, j)=rr(i, j)
62			uent(i, k, j)=u(i, j)
63			vent(i, k, j)=v(i, j)
64			elij(i, k, j)=0.0
65			!ym ment(i, k, j)=0.0
66			!ym sij(i, k, j)=0.0
67	guez	97	end do
68			end do
69			end do
70	guez	47
71	guez	97	!ym
72	guez	145	ment(1:ncum, 1:nd, 1:nd)=0.0
73			sij(1:ncum, 1:nd, 1:nd)=0.0
74	guez	47
75	guez	145	zm(:, :)=0.
76	guez	47
77	guez	97	! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
78			! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
79			! --- FRACTION (sij)
80	guez	47
81	guez	145	do i=minorig+1, nl
82	guez	47
83	guez	145	do j=minorig, nl
84			do il=1, ncum
85			if((i >= icb(il)).and.(i <= inb(il)).and. &
86			(j >= (icb(il)-1)).and.(j <= inb(il)))then
87	guez	47
88	guez	145	rti=rr(il, 1)-ep(il, i)*clw(il, i)
89			bf2=1.+lv(il, j)lv(il, j)rs(il, j)/(rrvt(il, j)t(il, j)*cpd)
90			anum=h(il, j)-hp(il, i)+(cpv-cpd)t(il, j)(rti-rr(il, j))
91			denom=h(il, i)-hp(il, i)+(cpd-cpv)(rr(il, i)-rti)t(il, j)
92	guez	97	dei=denom
93	guez	145	if(abs(dei) < 0.01)dei=0.01
94			sij(il, i, j)=anum/dei
95			sij(il, i, i)=1.0
96			altem=sij(il, i, j)rr(il, i)+(1.-sij(il, i, j))rti-rs(il, j)
97	guez	97	altem=altem/bf2
98	guez	145	cwat=clw(il, j)*(1.-ep(il, j))
99			stemp=sij(il, i, j)
100			if((stemp < 0.0.or.stemp > 1.0.or.altem > cwat) &
101			.and.j > i)then
102			anum=anum-lv(il, j)(rti-rs(il, j)-cwatbf2)
103			denom=denom+lv(il, j)*(rr(il, i)-rti)
104			if(abs(denom) < 0.01)denom=0.01
105			sij(il, i, j)=anum/denom
106			altem=sij(il, i, j)rr(il, i)+(1.-sij(il, i, j))rti-rs(il, j)
107	guez	97	altem=altem-(bf2-1.)*cwat
108			end if
109	guez	145	if(sij(il, i, j) > 0.0.and.sij(il, i, j) < 0.95)then
110			qent(il, i, j)=sij(il, i, j)rr(il, i)+(1.-sij(il, i, j))rti
111			uent(il, i, j)=sij(il, i, j)u(il, i)+(1.-sij(il, i, j))u(il, nk(il))
112			vent(il, i, j)=sij(il, i, j)v(il, i)+(1.-sij(il, i, j))v(il, nk(il))
113			elij(il, i, j)=altem
114			elij(il, i, j)=amax1(0.0, elij(il, i, j))
115			ment(il, i, j)=m(il, i)/(1.-sij(il, i, j))
116			nent(il, i)=nent(il, i)+1
117	guez	97	end if
118	guez	145	sij(il, i, j)=amax1(0.0, sij(il, i, j))
119			sij(il, i, j)=amin1(1.0, sij(il, i, j))
120	guez	97	endif ! new
121			end do
122			end do
123	guez	47
124	guez	145	! * if no air can entrain at level i assume that updraft detrains *
125			! * at that level and calculate detrained air flux and properties *
126	guez	47
127	guez	145	!@ do 170 i=icb(il), inb(il)
128	guez	47
129	guez	145	do il=1, ncum
130			if ((i >= icb(il)).and.(i <= inb(il)).and.(nent(il, i) == 0)) then
131			!@ if(nent(il, i) == 0)then
132			ment(il, i, i)=m(il, i)
133			qent(il, i, i)=rr(il, nk(il))-ep(il, i)*clw(il, i)
134			uent(il, i, i)=u(il, nk(il))
135			vent(il, i, i)=v(il, nk(il))
136			elij(il, i, i)=clw(il, i)
137			!MAF sij(il, i, i)=1.0
138			sij(il, i, i)=0.0
139	guez	97	end if
140			end do
141			end do
142	guez	47
143	guez	145	do j=minorig, nl
144			do i=minorig, nl
145			do il=1, ncum
146			if ((j >= (icb(il)-1)).and.(j <= inb(il)) &
147			.and.(i >= icb(il)).and.(i <= inb(il)))then
148			sigij(il, i, j)=sij(il, i, j)
149	guez	97	endif
150			end do
151			end do
152			end do
153	guez	145	!@ enddo
154	guez	47
155	guez	145	!@170 continue
156	guez	47
157	guez	145	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
158			! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
159	guez	47
160	guez	145	call zilch(asum, nloc*nd)
161			call zilch(csum, nloc*nd)
162			call zilch(csum, nloc*nd)
163	guez	47
164	guez	145	do il=1, ncum
165	guez	47	lwork(il) = .FALSE.
166	guez	97	enddo
167	guez	47
168	guez	145	DO i=minorig+1, nl
169	guez	47
170	guez	97	num1=0
171	guez	145	do il=1, ncum
172			if (i >= icb(il) .and. i <= inb(il)) num1=num1+1
173	guez	97	enddo
174	guez	145	if (num1 <= 0) cycle
175	guez	47
176	guez	145	do il=1, ncum
177			if (i >= icb(il) .and. i <= inb(il)) then
178			lwork(il)=(nent(il, i) /= 0)
179			qp=rr(il, 1)-ep(il, i)*clw(il, i)
180			anum=h(il, i)-hp(il, i)-lv(il, i)*(qp-rs(il, i)) &
181			+(cpv-cpd)t(il, i)(qp-rr(il, i))
182			denom=h(il, i)-hp(il, i)+lv(il, i)*(rr(il, i)-qp) &
183			+(cpd-cpv)t(il, i)(rr(il, i)-qp)
184			if(abs(denom) < 0.01)denom=0.01
185	guez	97	scrit(il)=anum/denom
186	guez	145	alt=qp-rs(il, i)+scrit(il)*(rr(il, i)-qp)
187			if(scrit(il) <= 0.0.or.alt <= 0.0)scrit(il)=1.0
188	guez	97	smax(il)=0.0
189			asij(il)=0.0
190			endif
191			end do
192	guez	47
193	guez	145	do j=nl, minorig, -1
194	guez	47
195	guez	97	num2=0
196	guez	145	do il=1, ncum
197			if (i >= icb(il) .and. i <= inb(il) .and. &
198			j >= (icb(il)-1) .and. j <= inb(il) &
199			.and. lwork(il)) num2=num2+1
200	guez	97	enddo
201	guez	145	if (num2 <= 0) cycle
202	guez	47
203	guez	145	do il=1, ncum
204			if (i >= icb(il) .and. i <= inb(il) .and. &
205			j >= (icb(il)-1) .and. j <= inb(il) &
206			.and. lwork(il)) then
207	guez	47
208	guez	145	if(sij(il, i, j) > 1.0e-16.and.sij(il, i, j) < 0.95)then
209	guez	97	wgh=1.0
210	guez	145	if(j > i)then
211			sjmax=amax1(sij(il, i, j+1), smax(il))
212			sjmax=amin1(sjmax, scrit(il))
213			smax(il)=amax1(sij(il, i, j), smax(il))
214			sjmin=amax1(sij(il, i, j-1), smax(il))
215			sjmin=amin1(sjmin, scrit(il))
216			if(sij(il, i, j) < (smax(il)-1.0e-16))wgh=0.0
217			smid=amin1(sij(il, i, j), scrit(il))
218	guez	97	else
219	guez	145	sjmax=amax1(sij(il, i, j+1), scrit(il))
220			smid=amax1(sij(il, i, j), scrit(il))
221	guez	97	sjmin=0.0
222	guez	145	if(j > 1)sjmin=sij(il, i, j-1)
223			sjmin=amax1(sjmin, scrit(il))
224	guez	97	endif
225			delp=abs(sjmax-smid)
226			delm=abs(sjmin-smid)
227			asij(il)=asij(il)+wgh*(delp+delm)
228	guez	145	ment(il, i, j)=ment(il, i, j)(delp+delm)wgh
229	guez	97	endif
230			endif
231			end do
232	guez	47
233	guez	97	end do
234	guez	47
235	guez	145	do il=1, ncum
236			if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
237			asij(il)=amax1(1.0e-16, asij(il))
238	guez	97	asij(il)=1.0/asij(il)
239	guez	145	asum(il, i)=0.0
240			bsum(il, i)=0.0
241			csum(il, i)=0.0
242	guez	97	endif
243	guez	47	enddo
244
245	guez	145	do j=minorig, nl
246			do il=1, ncum
247			if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
248			.and. j >= (icb(il)-1) .and. j <= inb(il)) then
249			ment(il, i, j)=ment(il, i, j)*asij(il)
250	guez	97	endif
251			enddo
252			end do
253	guez	47
254	guez	145	do j=minorig, nl
255			do il=1, ncum
256			if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
257			.and. j >= (icb(il)-1) .and. j <= inb(il)) then
258			asum(il, i)=asum(il, i)+ment(il, i, j)
259			ment(il, i, j)=ment(il, i, j)*sig(il, j)
260			bsum(il, i)=bsum(il, i)+ment(il, i, j)
261	guez	97	endif
262			enddo
263			end do
264	guez	47
265	guez	145	do il=1, ncum
266			if (i >= icb(il).and.i <= inb(il).and.lwork(il)) then
267			bsum(il, i)=amax1(bsum(il, i), 1.0e-16)
268			bsum(il, i)=1.0/bsum(il, i)
269	guez	97	endif
270	guez	47	enddo
271
272	guez	145	do j=minorig, nl
273			do il=1, ncum
274			if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
275			.and. j >= (icb(il)-1) .and. j <= inb(il)) then
276			ment(il, i, j)=ment(il, i, j)asum(il, i)bsum(il, i)
277	guez	97	endif
278			enddo
279			end do
280
281	guez	145	do j=minorig, nl
282			do il=1, ncum
283			if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
284			.and. j >= (icb(il)-1) .and. j <= inb(il)) then
285			csum(il, i)=csum(il, i)+ment(il, i, j)
286	guez	97	endif
287			enddo
288			end do
289
290	guez	145	do il=1, ncum
291			if (i >= icb(il) .and. i <= inb(il) .and. lwork(il) &
292			.and. csum(il, i) < m(il, i)) then
293			nent(il, i)=0
294			ment(il, i, i)=m(il, i)
295			qent(il, i, i)=rr(il, 1)-ep(il, i)*clw(il, i)
296			uent(il, i, i)=u(il, nk(il))
297			vent(il, i, i)=v(il, nk(il))
298			elij(il, i, i)=clw(il, i)
299			!MAF sij(il, i, i)=1.0
300			sij(il, i, i)=0.0
301	guez	97	endif
302			enddo ! il
303	guez	47
304	guez	97	end DO
305	guez	145
306	guez	97	! MAF: renormalisation de MENT
307	guez	145	do jm=1, nd
308			do im=1, nd
309			do il=1, ncum
310			zm(il, im)=zm(il, im)+(1.-sij(il, im, jm))*ment(il, im, jm)
311	guez	97	end do
312			end do
313			end do
314	guez	145
315			do jm=1, nd
316			do im=1, nd
317			do il=1, ncum
318			if(zm(il, im) /= 0.) then
319			ment(il, im, jm)=ment(il, im, jm)*m(il, im)/zm(il, im)
320	guez	97	endif
321			end do
322	guez	47	end do
323	guez	97	end do
324	guez	145
325			do jm=1, nd
326			do im=1, nd
327			do il=1, ncum
328			qents(il, im, jm)=qent(il, im, jm)
329			ments(il, im, jm)=ment(il, im, jm)
330	guez	97	end do
331	guez	47	enddo
332	guez	97	enddo
333	guez	47
334	guez	97	end SUBROUTINE cv3_mixing
335
336			end module cv3_mixing_m