phylmd/CV30_routines/cv30_mixing.f90

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(:) ! (ncum) {2 <= icb <= nl - 3}

    integer, intent(in):: inb(:) ! (ncum)
    ! first model level above the level of neutral buoyancy of the
    ! parcel (1 <= inb <= nl - 1)

    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, intent(out):: nent(:, 2:) ! (ncum, 2:nl - 1)
    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

    nent = 0

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