Context Navigation

cv3p_mixing.F90 @ 2415

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

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

Download in other formats:

Original Format