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

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