phylmd/CV30_routines/cv30_mixing.f

module cv30_mixing_m

  implicit none

contains

  SUBROUTINE cv30_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 cv30_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, intent(in):: 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 cv30_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 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

    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)) &
                  +(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 cv30_mixing

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