Context Navigation

cv3p_mixing.F90 @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1997:2055 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 20.1 KB

Line
1	SUBROUTINE cv3p_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, &
2	ph, t, rr, rs, u, v, tra, h, lv, qnk, &
3	unk, vnk, hp, tv, tvp, ep, clw, sig, &
4	Ment, Qent, hent, uent, vent, nent, &
5	Sigij, elij, supmax, Ments, Qents, traent)
6	! **************************************************************
7	! *
8	! CV3P_MIXING : compute mixed draught properties and, *
9	! within a scaling factor, mixed draught *
10	! mass fluxes. *
11	! written by : VTJ Philips,JY Grandpeix, 21/05/2003, 09.14.15*
12	! modified by : *
13	! **************************************************************
14
15	IMPLICIT NONE
16
17	include "cvthermo.h"
18	include "cv3param.h"
19	include "YOMCST2.h"
20
21	!inputs:
22	INTEGER, INTENT (IN) :: ncum, nd, na
23	INTEGER, INTENT (IN) :: ntra, nloc
24	INTEGER, DIMENSION (nloc), INTENT (IN) :: icb, inb, nk
25	REAL, DIMENSION (nloc, nd), INTENT (IN) :: sig
26	REAL, DIMENSION (nloc), INTENT (IN) :: qnk, unk, vnk
27	REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph
28	REAL, DIMENSION (nloc, nd), INTENT (IN) :: t, rr, rs
29	REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v
30	REAL, DIMENSION (nloc, nd, ntra), INTENT (IN) :: tra ! input of convect3
31	REAL, DIMENSION (nloc, na), INTENT (IN) :: lv
32	REAL, DIMENSION (nloc, na), INTENT (IN) :: h !liquid water static energy of environMent
33	REAL, DIMENSION (nloc, na), INTENT (IN) :: hp !liquid water static energy of air shed from adiab. asc.
34	REAL, DIMENSION (nloc, na), INTENT (IN) :: tv, tvp
35	REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, clw
36
37	!outputs:
38	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Ment, Qent
39	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: uent, vent
40	REAL, DIMENSION (nloc, na, na), INTENT (OUT) :: Sigij, elij
41	REAL, DIMENSION (nloc, na), INTENT (OUT) :: supmax(nloc, na) ! Highest mixing fraction of mixed
42	! updraughts with the sign of (h-hp)
43	REAL, DIMENSION (nloc, nd, nd, ntra), INTENT (OUT) :: traent
44	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: Ments, Qents
45	REAL, DIMENSION (nloc, nd, nd), INTENT (OUT) :: hent
46	INTEGER, DIMENSION (nloc, nd), INTENT (OUT) :: nent
47
48	!local variables:
49	INTEGER i, j, k, il, im, jm
50	INTEGER num1, num2
51	REAL :: rti, bf2, anum, denom, dei, altem, cwat, stemp, qp
52	REAL :: alt, delp, delm
53	REAL, DIMENSION (nloc) :: Qmixmax, Rmixmax, sqmrmax
54	REAL, DIMENSION (nloc) :: Qmixmin, Rmixmin, sqmrmin
55	REAL, DIMENSION (nloc) :: signhpmh
56	REAL, DIMENSION (nloc) :: Sx
57	REAL :: Scrit2
58	REAL, DIMENSION (nloc) :: Smid, Sjmin, Sjmax
59	REAL, DIMENSION (nloc) :: Sbef, sup, smin
60	REAL, DIMENSION (nloc) :: ASij, smax, Scrit
61	REAL, DIMENSION (nloc, nd, nd) :: Sij
62	REAL, DIMENSION (nloc, nd) :: csum
63	REAL :: awat
64	LOGICAL, DIMENSION (nloc) :: lwork
65
66	REAL amxupcrit, df, ff
67	INTEGER nstep
68
69	! -- Mixing probability distribution functions
70
71	REAL Qcoef1, Qcoef2, QFF, QFFF, Qmix, Rmix, Qmix1, Rmix1, Qmix2, Rmix2, F
72
73	Qcoef1(F) = tanh(F/gammas)
74	Qcoef2(F) = (tanh(F/gammas)+gammas*log(cosh((1.-F)/gammas)/cosh(F/gammas)))
75	QFF(F) = max(min(F,1.), 0.)
76	QFFf(F) = min(QFF(F), scut)
77	Qmix1(F) = (tanh((QFF(F)-Fmax)/gammas)+Qcoef1max)/Qcoef2max
78	Rmix1(F) = (gammaslog(cosh((QFF(F)-Fmax)/gammas))+QFF(F)Qcoef1max)/Qcoef2max
79	Qmix2(F) = -log(1.-QFFf(F))/scut
80	Rmix2(F) = (QFFf(F)+(1.-QFF(F))*log(1.-QFFf(F)))/scut
81	Qmix(F) = qqa1Qmix1(F) + qqa2Qmix2(F)
82	Rmix(F) = qqa1Rmix1(F) + qqa2Rmix2(F)
83
84	INTEGER, SAVE :: ifrst
85	DATA ifrst/0/
86	!$OMP THREADPRIVATE(ifrst)
87
88
89	! =====================================================================
90	! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
91	! =====================================================================
92
93	! -- Initialize mixing PDF coefficients
94	IF (ifrst==0) THEN
95	ifrst = 1
96	Qcoef1max = Qcoef1(Fmax)
97	Qcoef2max = Qcoef2(Fmax)
98
99	END IF
100
101
102	! ori do 360 i=1,ncum*nlp
103	DO j = 1, nl
104	DO i = 1, ncum
105	nent(i, j) = 0
106	! in convect3, m is computed in cv3_closure
107	! ori m(i,1)=0.0
108	END DO
109	END DO
110
111	! ori do 400 k=1,nlp
112	! ori do 390 j=1,nlp
113	DO j = 1, nl
114	DO k = 1, nl
115	DO i = 1, ncum
116	Qent(i, k, j) = rr(i, j)
117	uent(i, k, j) = u(i, j)
118	vent(i, k, j) = v(i, j)
119	elij(i, k, j) = 0.0
120	hent(i, k, j) = 0.0
121	!AC! Ment(i,k,j)=0.0
122	!AC! Sij(i,k,j)=0.0
123	END DO
124	END DO
125	END DO
126
127	!AC!
128	Ment(1:ncum, 1:nd, 1:nd) = 0.0
129	Sij(1:ncum, 1:nd, 1:nd) = 0.0
130	!AC!
131
132	DO k = 1, ntra
133	DO j = 1, nd ! instead nlp
134	DO i = 1, nd ! instead nlp
135	DO il = 1, ncum
136	traent(il, i, j, k) = tra(il, j, k)
137	END DO
138	END DO
139	END DO
140	END DO
141
142	! =====================================================================
143	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
144	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
145	! --- FRACTION (Sij)
146	! =====================================================================
147
148	DO i = minorig + 1, nl
149
150	DO j = minorig, nl
151	DO il = 1, ncum
152	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
153	.AND. (j<=inb(il))) THEN
154
155	rti = qnk(il) - ep(il, i)*clw(il, i)
156	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
157	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
158	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
159	dei = denom
160	IF (abs(dei)<0.01) dei = 0.01
161	Sij(il, i, j) = anum/dei
162	Sij(il, i, i) = 1.0
163	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
164	altem = altem/bf2
165	cwat = clw(il, j)*(1.-ep(il,j))
166	stemp = Sij(il, i, j)
167	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
168	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
169	denom = denom + lv(il, j)*(rr(il,i)-rti)
170	IF (abs(denom)<0.01) denom = 0.01
171	Sij(il, i, j) = anum/denom
172	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
173	altem = altem - (bf2-1.)*cwat
174	END IF
175	IF (Sij(il,i,j)>0.0) THEN
176	!!! Ment(il,i,j)=m(il,i)
177	Ment(il, i, j) = 1.
178	elij(il, i, j) = altem
179	elij(il, i, j) = amax1(0.0, elij(il,i,j))
180	nent(il, i) = nent(il, i) + 1
181	END IF
182
183	Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
184	Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
185	END IF ! new
186	END DO
187	END DO
188
189
190	! * if no air can entrain at level i assume that updraft detrains *
191	! * at that level and calculate detrained air flux and properties *
192
193
194	! @ do 170 i=icb(il),inb(il)
195
196	DO il = 1, ncum
197	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
198	! @ if(nent(il,i).eq.0)then
199	!!! Ment(il,i,i)=m(il,i)
200	Ment(il, i, i) = 1.
201	Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
202	uent(il, i, i) = unk(il)
203	vent(il, i, i) = vnk(il)
204	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
205	Sij(il, i, i) = 0.0
206	END IF
207	END DO
208	END DO
209
210	DO j = 1, ntra
211	DO i = minorig + 1, nl
212	DO il = 1, ncum
213	IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
214	traent(il, i, i, j) = tra(il, nk(il), j)
215	END IF
216	END DO
217	END DO
218	END DO
219
220	DO j = minorig, nl
221	DO i = minorig, nl
222	DO il = 1, ncum
223	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
224	(i>=icb(il)) .AND. (i<=inb(il))) THEN
225	Sigij(il, i, j) = Sij(il, i, j)
226	END IF
227	END DO
228	END DO
229	END DO
230	! @ enddo
231
232	! @170 continue
233
234	! =====================================================================
235	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
236	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
237	! =====================================================================
238
239	CALL zilch(csum, nloc*nd)
240
241	DO il = 1, ncum
242	lwork(il) = .FALSE.
243	END DO
244
245	! ---------------------------------------------------------------
246	DO i = minorig + 1, nl !Loop on origin level "i"
247	! ---------------------------------------------------------------
248
249	num1 = 0
250	DO il = 1, ncum
251	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
252	END DO
253	IF (num1<=0) GO TO 789
254
255
256	!JYG1 Find maximum of SIJ for J>I, if any.
257
258	Sx(:) = 0.
259
260	DO il = 1, ncum
261	IF (i>=icb(il) .AND. i<=inb(il)) THEN
262	signhpmh(il) = sign(1., hp(il,i)-h(il,i))
263	Sbef(il) = max(0., signhpmh(il))
264	END IF
265	END DO
266
267	DO j = i + 1, nl
268	DO il = 1, ncum
269	IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
270	IF (Sbef(il)<Sij(il,i,j)) THEN
271	Sx(il) = max(Sij(il,i,j), Sx(il))
272	END IF
273	Sbef(il) = Sij(il, i, j)
274	END IF
275	END DO
276	END DO
277
278
279	DO il = 1, ncum
280	IF (i>=icb(il) .AND. i<=inb(il)) THEN
281	lwork(il) = (nent(il,i)/=0)
282	rti = qnk(il) - ep(il, i)*clw(il, i)
283	anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
284	(cpv-cpd)t(il, i)(rti-rr(il,i))
285	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
286	(cpd-cpv)t(il, i)(rr(il,i)-rti)
287	IF (abs(denom)<0.01) denom = 0.01
288	Scrit(il) = min(anum/denom, 1.)
289	alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
290
291	!JYG1 Find new critical value Scrit2
292	! such that : Sij > Scrit2 => mixed draught will detrain at J<I
293	! Sij < Scrit2 => mixed draught will detrain at J>I
294
295	Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
296	Scrit(il)*max(0., signhpmh(il))
297
298	Scrit(il) = Scrit2
299
300	!JYG Correction pour la nouvelle logique; la correction pour ALT
301	! est un peu au hazard
302	IF (Scrit(il)<=0.0) Scrit(il) = 0.0
303	IF (alt<=0.0) Scrit(il) = 1.0
304
305	smax(il) = 0.0
306	ASij(il) = 0.0
307	sup(il) = 0. ! upper S-value reached by descending draughts
308	END IF
309	END DO
310
311	! ---------------------------------------------------------------
312	DO j = minorig, nl !Loop on destination level "j"
313	! ---------------------------------------------------------------
314
315	num2 = 0
316	DO il = 1, ncum
317	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
318	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
319	lwork(il)) num2 = num2 + 1
320	END DO
321	IF (num2<=0) GO TO 175
322
323	! -----------------------------------------------
324	IF (j>i) THEN
325	! -----------------------------------------------
326	DO il = 1, ncum
327	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
328	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
329	lwork(il)) THEN
330	IF (Sij(il,i,j)>0.0) THEN
331	Smid(il) = min(Sij(il,i,j), Scrit(il))
332	Sjmax(il) = Smid(il)
333	Sjmin(il) = Smid(il)
334	IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
335	smin(il) = Smid(il)
336	Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
337	Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
338	Sjmin(il) = min(Sjmin(il), Scrit(il))
339	Sbef(il) = Sij(il, i, j)
340	END IF
341	END IF
342	END IF
343	END DO
344	! -----------------------------------------------
345	ELSE IF (j==i) THEN
346	! -----------------------------------------------
347	DO il = 1, ncum
348	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
349	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
350	lwork(il)) THEN
351	IF (Sij(il,i,j)>0.0) THEN
352	Smid(il) = 1.
353	Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
354	min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
355	Sjmin(il) = max(Sjmin(il), sup(il))
356	Sjmax(il) = 1.
357
358	! - preparation des variables Scrit, Smin et Sbef pour la partie j>i
359	Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
360
361	smin(il) = 1.
362	Sbef(il) = max(0., signhpmh(il))
363	supmax(il, i) = sign(Scrit(il), -signhpmh(il))
364	END IF
365	END IF
366	END DO
367	! -----------------------------------------------
368	ELSE IF (j<i) THEN
369	! -----------------------------------------------
370	DO il = 1, ncum
371	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
372	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
373	lwork(il)) THEN
374	IF (Sij(il,i,j)>0.0) THEN
375	Smid(il) = max(Sij(il,i,j), Scrit(il))
376	Sjmax(il) = Smid(il)
377	Sjmin(il) = Smid(il)
378	IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
379	smax(il) = Smid(il)
380	Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
381	Sjmax(il) = max(Sjmax(il), Scrit(il))
382	Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
383	Sjmin(il) = max(Sjmin(il), Scrit(il))
384	Sbef(il) = Sij(il, i, j)
385	END IF
386	IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
387	sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
388	END IF
389	END IF
390	END DO
391	! -----------------------------------------------
392	END IF
393	! -----------------------------------------------
394
395
396	DO il = 1, ncum
397	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
398	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
399	lwork(il)) THEN
400	IF (Sij(il,i,j)>0.0) THEN
401	rti = qnk(il) - ep(il, i)*clw(il, i)
402	Qmixmax(il) = Qmix(Sjmax(il))
403	Qmixmin(il) = Qmix(Sjmin(il))
404	Rmixmax(il) = Rmix(Sjmax(il))
405	Rmixmin(il) = Rmix(Sjmin(il))
406	sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
407	sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
408
409	Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
410
411	! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
412	IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
413	Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
414	ELSE
415	Sigij(il, i, j) = 0.
416	END IF
417
418	! -- Compute Qent, uent, vent according to the true mixing fraction
419	Qent(il, i, j) = (1.-Sigij(il,i,j))rti + Sigij(il, i, j)rr(il, i)
420	uent(il, i, j) = (1.-Sigij(il,i,j))unk(il) + Sigij(il, i, j)u(il, i)
421	vent(il, i, j) = (1.-Sigij(il,i,j))vnk(il) + Sigij(il, i, j)v(il, i)
422
423	! -- Compute liquid water static energy of mixed draughts
424	! IF (j .GT. i) THEN
425	! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
426	! awat=amax1(awat,0.0)
427	! ELSE
428	! awat = 0.
429	! ENDIF
430	! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
431	! : + Sigij(il,i,j)*H(il,i)
432	! : + (LV(il,j)+(cpd-cpv)t(il,j))awat
433	!IM 301008 beg
434	hent(il, i, j) = (1.-Sigij(il,i,j))hp(il, i) + Sigij(il, i, j)h(il, i)
435
436	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
437	elij(il, i, j) = elij(il, i, j) + &
438	((h(il,j)-hent(il,i,j))rs(il,j)lv(il,j) / &
439	((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
440	elij(il, i, j) = elij(il, i, j) / &
441	(1.+lv(il,j)lv(il,j)rs(il,j) / &
442	((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
443
444	elij(il, i, j) = max(elij(il,i,j), 0.)
445
446	elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
447
448	IF (j>i) THEN
449	awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
450	awat = amax1(awat, 0.0)
451	ELSE
452	awat = 0.
453	END IF
454
455	! print ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
456	! : t(il,j))
457
458	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)t(il,j))awat
459	!IM 301008 end
460
461	! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
462	! : i,j,hent(il,i,j),Sigij(il,i,j)
463
464	! -- ASij is the integral of P(F) over the relevant F interval
465	ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
466	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
467
468	END IF
469	END IF
470	END DO
471	DO k = 1, ntra
472	DO il = 1, ncum
473	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
474	(j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
475	lwork(il)) THEN
476	IF (Sij(il,i,j)>0.0) THEN
477	traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
478	(1.-Sigij(il,i,j))*tra(il, nk(il), k)
479	END IF
480	END IF
481	END DO
482	END DO
483
484	! -- If I=J (detrainement and entrainement at the same level), then only the
485	! -- adiabatic ascent part of the mixture is considered
486	IF (i==j) THEN
487	DO il = 1, ncum
488	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
489	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
490	lwork(il)) THEN
491	IF (Sij(il,i,j)>0.0) THEN
492	rti = qnk(il) - ep(il, i)*clw(il, i)
493	!!! Ment(il,i,i) = m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
494	Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
495	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
496	Qent(il, i, i) = rti
497	uent(il, i, i) = unk(il)
498	vent(il, i, i) = vnk(il)
499	hent(il, i, i) = hp(il, i)
500	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
501	Sigij(il, i, i) = 0.
502	END IF
503	END IF
504	END DO
505	DO k = 1, ntra
506	DO il = 1, ncum
507	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
508	(j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
509	lwork(il)) THEN
510	IF (Sij(il,i,j)>0.0) THEN
511	traent(il, i, i, k) = tra(il, nk(il), k)
512	END IF
513	END IF
514	END DO
515	END DO
516
517	END IF
518
519	! ---------------------------------------------------------------
520	175 END DO ! End loop on destination level "j"
521	! ---------------------------------------------------------------
522
523	DO il = 1, ncum
524	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
525	ASij(il) = amax1(1.0E-16, ASij(il))
526	ASij(il) = 1.0/ASij(il)
527	csum(il, i) = 0.0
528	END IF
529	END DO
530
531	DO j = minorig, nl
532	DO il = 1, ncum
533	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
534	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
535	Ment(il, i, j) = Ment(il, i, j)*ASij(il)
536	END IF
537	END DO
538	END DO
539
540	DO j = minorig, nl
541	DO il = 1, ncum
542	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
543	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
544	csum(il, i) = csum(il, i) + Ment(il, i, j)
545	END IF
546	END DO
547	END DO
548
549	DO il = 1, ncum
550	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
551	! cc : .and. csum(il,i).lt.m(il,i) ) then
552	nent(il, i) = 0
553	! cc Ment(il,i,i)=m(il,i)
554	Ment(il, i, i) = 1.
555	Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
556	uent(il, i, i) = unk(il)
557	vent(il, i, i) = vnk(il)
558	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
559	Sij(il, i, i) = 0.0
560	END IF
561	END DO ! il
562
563	DO j = 1, ntra
564	DO il = 1, ncum
565	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
566	! cc : .and. csum(il,i).lt.m(il,i) ) then
567	traent(il, i, i, j) = tra(il, nk(il), j)
568	END IF
569	END DO
570	END DO
571
572	! ---------------------------------------------------------------
573	789 END DO ! End loop on origin level "i"
574	! ---------------------------------------------------------------
575
576
577	RETURN
578	END SUBROUTINE cv3p_mixing
579

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format