Context Navigation

cv3p_mixing.f90 @ 5322

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

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

Download in other formats:

Original Format