Context Navigation

cv3p_mixing.F90 @ 2408

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

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

Download in other formats:

Original Format