Context Navigation

cv3p_mixing.F90 @ 5440

Last change on this file since 5440 was 2905, checked in by musat, 8 years ago
Bug correction on ok_entrain making an MPI error
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: 24.8 KB

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

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

Download in other formats:

Original Format