Context Navigation

cv3p_mixing.F90 @ 2380

Last change on this file since 2380 was 2226, checked in by Laurent Fairhead, 10 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