phylmd/CV30_routines/cv30_mixing.f

module cv30_mixing_m

  implicit none

contains

  SUBROUTINE cv30_mixing(icb, inb, t, rr, rs, u, v, h, lv, hp, ep, clw, m, &
       sig, ment, qent, uent, vent, nent, sij, elij, ments, qents)

    ! MIXING

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

    use cv30_param_m, only: minorig, nl
    USE dimphy, ONLY: klev, klon
    use suphec_m, only: rcpd, rcpv, rv

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

    ! outputs:
    real ment(klon, klev, klev), qent(klon, klev, klev)
    real uent(klon, klev, klev), vent(klon, klev, klev)
    integer nent(klon, klev)
    real sij(klon, klev, klev), elij(klon, klev, klev)
    real ments(klon, klev, klev), qents(klon, klev, klev)

    ! Local:
    integer ncum, 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(klon), smax(klon), scrit(klon)
    real asum(klon, klev), bsum(klon, klev), csum(klon, klev)
    real wgh
    real zm(klon, klev)
    logical lwork(klon)

    !-------------------------------------------------------------------------

    ncum = size(icb)

    ! INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS

    do j = 1, nl
       do i = 1, ncum
          nent(i, j) = 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
          end do
       end do
    end do

    ment(1:ncum, 1:klev, 1:klev) = 0.0
    sij(1:ncum, 1:klev, 1:klev) = 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) / (rv &
                     * t(il, j) * t(il, j) * rcpd)
                anum = h(il, j) - hp(il, i) + (rcpv - rcpd) * t(il, j) * (rti &
                     - rr(il, j))
                denom = h(il, i) - hp(il, i) + (rcpd - rcpv) * (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, minorig)
                   vent(il, i, j) = sij(il, i, j) * v(il, i) + (1. &
                        - sij(il, i, j)) * v(il, minorig)
                   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 il = 1, ncum
          if ((i >= icb(il)).and.(i <= inb(il)).and.(nent(il, i) == 0)) then
             ment(il, i, i) = m(il, i)
             qent(il, i, i) = rr(il, minorig) - ep(il, i) * clw(il, i)
             uent(il, i, i) = u(il, minorig)
             vent(il, i, i) = v(il, minorig)
             elij(il, i, i) = clw(il, i)
             sij(il, i, i) = 0.0
          end if
       end do
    end do

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

    asum = 0.
    csum = 0.

    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)) &
                  + (rcpv - rcpd) * t(il, i) * (qp - rr(il, i))
             denom = h(il, i) - hp(il, i) + lv(il, i) * (rr(il, i) - qp) &
                  + (rcpd - rcpv) * 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, minorig)
             vent(il, i, i) = v(il, minorig)
             elij(il, i, i) = clw(il, i)
             sij(il, i, i) = 0.0
          endif
       enddo ! il

    end DO

    ! MAF: renormalisation de MENT
    do jm = 1, klev
       do im = 1, klev
          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, klev
       do im = 1, klev
          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, klev
       do im = 1, klev
          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 cv30_mixing

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