Context Navigation

cv3p_mixing.F90 @ 2226

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

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

Download in other formats:

Original Format