Context Navigation

cv3p_mixing.F90 @ 5443

Last change on this file since 5443 was 2478, checked in by fhourdin, 9 years ago
Correction d'un bug de non initialisation d'un tableau qui faisait planter le modèle en haute résolution sur 8 thread. Bug identifié par Yann Meurdesoif.
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	!ym
146	Sigij(1:ncum, 1:nd, 1:nd) = 0.0
147	!ym
148
149	!jyg! DO k = 1, ntra
150	!jyg! DO j = 1, nd ! instead nlp
151	!jyg! DO i = 1, nd ! instead nlp
152	!jyg! DO il = 1, ncum
153	!jyg! traent(il, i, j, k) = tra(il, j, k)
154	!jyg! END DO
155	!jyg! END DO
156	!jyg! END DO
157	!jyg! END DO
158
159	! =====================================================================
160	! --- CALCULATE ENTRAINED AIR MASS FLUX (Ment), TOTAL WATER MIXING
161	! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
162	! --- FRACTION (Sij)
163	! =====================================================================
164
165	DO i = minorig + 1, nl
166
167	DO j = minorig, nl
168	DO il = 1, ncum
169	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)-1)) &
170	.AND. (j<=inb(il))) THEN
171
172	rti = qnk(il) - ep(il, i)*clw(il, i)
173	bf2 = 1. + lv(il, j)lv(il, j)rs(il, j)/(rrvt(il,j)t(il,j)*cpd)
174	!jyg(from aj)<
175	IF (cvflag_ice) THEN
176	! print*,cvflag_ice,'cvflag_ice dans do 700'
177	IF (t(il,j)<=263.15) THEN
178	bf2 = 1. + (lf(il,j)+lv(il,j))(lv(il,j)+frac(il,j) &
179	lf(il,j))rs(il, j)/(rrvt(il,j)t(il,j)cpd)
180	END IF
181	END IF
182	!>jyg
183	anum = h(il, j) - hp(il, i) + (cpv-cpd)t(il, j)(rti-rr(il,j))
184	denom = h(il, i) - hp(il, i) + (cpd-cpv)(rr(il,i)-rti)t(il, j)
185	dei = denom
186	IF (abs(dei)<0.01) dei = 0.01
187	Sij(il, i, j) = anum/dei
188	Sij(il, i, i) = 1.0
189	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
190	altem = altem/bf2
191	cwat = clw(il, j)*(1.-ep(il,j))
192	stemp = Sij(il, i, j)
193	IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
194	!jyg(from aj)<
195	IF (cvflag_ice) THEN
196	anum = anum - (lv(il,j)+frac(il,j)lf(il,j))(rti-rs(il,j)-cwat*bf2)
197	denom = denom + (lv(il,j)+frac(il,j)lf(il,j))(rr(il,i)-rti)
198	ELSE
199	anum = anum - lv(il, j)(rti-rs(il,j)-cwatbf2)
200	denom = denom + lv(il, j)*(rr(il,i)-rti)
201	END IF
202	!>jyg
203	IF (abs(denom)<0.01) denom = 0.01
204	Sij(il, i, j) = anum/denom
205	altem = Sij(il, i, j)rr(il, i) + (1.-Sij(il,i,j))rti - rs(il, j)
206	altem = altem - (bf2-1.)*cwat
207	END IF
208	IF (Sij(il,i,j)>0.0) THEN
209	!!! Ment(il,i,j)=m(il,i)
210	Ment(il, i, j) = 1.
211	elij(il, i, j) = altem
212	elij(il, i, j) = amax1(0.0, elij(il,i,j))
213	nent(il, i) = nent(il, i) + 1
214	END IF
215
216	Sij(il, i, j) = amax1(0.0, Sij(il,i,j))
217	Sij(il, i, j) = amin1(1.0, Sij(il,i,j))
218	END IF ! new
219	END DO
220	END DO
221
222	!jygdebug<
223	IF (prt_level >= 10) THEN
224	print *,'cv3p_mixing i, nent(i), icb, inb ',i, nent(igout,i), icb(igout), inb(igout)
225	IF (nent(igout,i) .gt. 0) THEN
226	print *,'i,(j,Sij(i,j),j=icb-1,inb) ',i,(j,Sij(igout,i,j),j=icb(igout)-1,inb(igout))
227	ENDIF
228	ENDIF
229	!>jygdebug
230
231	! * if no air can entrain at level i assume that updraft detrains *
232	! * at that level and calculate detrained air flux and properties *
233
234
235	! @ do 170 i=icb(il),inb(il)
236
237	DO il = 1, ncum
238	IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN
239	! @ if(nent(il,i).eq.0)then
240	!!! Ment(il,i,i)=m(il,i)
241	Ment(il, i, i) = 1.
242	Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
243	uent(il, i, i) = unk(il)
244	vent(il, i, i) = vnk(il)
245	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
246	Sij(il, i, i) = 0.0
247	END IF
248	END DO
249	END DO
250
251	!jyg! DO j = 1, ntra
252	!jyg! DO i = minorig + 1, nl
253	!jyg! DO il = 1, ncum
254	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. nent(il,i)==0) THEN
255	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
256	!jyg! END IF
257	!jyg! END DO
258	!jyg! END DO
259	!jyg! END DO
260
261	DO j = minorig, nl
262	DO i = minorig, nl
263	DO il = 1, ncum
264	IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
265	(i>=icb(il)) .AND. (i<=inb(il))) THEN
266	Sigij(il, i, j) = Sij(il, i, j)
267	END IF
268	END DO
269	END DO
270	END DO
271	! @ enddo
272
273	! @170 continue
274
275	! =====================================================================
276	! --- NORMALIZE ENTRAINED AIR MASS FLUXES
277	! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING
278	! =====================================================================
279
280	CALL zilch(csum, nloc*nd)
281
282	DO il = 1, ncum
283	lwork(il) = .FALSE.
284	END DO
285
286	! ---------------------------------------------------------------
287	DO i = minorig + 1, nl !Loop on origin level "i"
288	! ---------------------------------------------------------------
289
290	num1 = 0
291	DO il = 1, ncum
292	IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1
293	END DO
294	IF (num1<=0) GO TO 789
295
296
297	!JYG1 Find maximum of SIJ for J>I, if any.
298
299	Sx(:) = 0.
300
301	DO il = 1, ncum
302	IF (i>=icb(il) .AND. i<=inb(il)) THEN
303	signhpmh(il) = sign(1., hp(il,i)-h(il,i))
304	Sbef(il) = max(0., signhpmh(il))
305	END IF
306	END DO
307
308	DO j = i + 1, nl
309	DO il = 1, ncum
310	IF (i>=icb(il) .AND. i<=inb(il) .AND. j<=inb(il)) THEN
311	IF (Sbef(il)<Sij(il,i,j)) THEN
312	Sx(il) = max(Sij(il,i,j), Sx(il))
313	END IF
314	Sbef(il) = Sij(il, i, j)
315	END IF
316	END DO
317	END DO
318
319
320	DO il = 1, ncum
321	IF (i>=icb(il) .AND. i<=inb(il)) THEN
322	lwork(il) = (nent(il,i)/=0)
323	rti = qnk(il) - ep(il, i)*clw(il, i)
324	!jyg<
325	IF (cvflag_ice) THEN
326
327	anum = h(il, i) - hp(il, i) - (lv(il,i)+frac(il,i)lf(il,i)) &
328	(rti-rs(il,i)) + (cpv-cpd)t(il, i)(rti-rr(il,i))
329	denom = h(il, i) - hp(il, i) + (lv(il,i)+frac(il,i)lf(il,i)) &
330	(rr(il,i)-rti) + (cpd-cpv)t(il, i)(rr(il,i)-rti)
331	ELSE
332
333	anum = h(il, i) - hp(il, i) - lv(il, i)*(rti-rs(il,i)) + &
334	(cpv-cpd)t(il, i)(rti-rr(il,i))
335	denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-rti) + &
336	(cpd-cpv)t(il, i)(rr(il,i)-rti)
337	END IF
338	!>jyg
339	IF (abs(denom)<0.01) denom = 0.01
340	Scrit(il) = min(anum/denom, 1.)
341	alt = rti - rs(il, i) + Scrit(il)*(rr(il,i)-rti)
342
343	!JYG1 Find new critical value Scrit2
344	! such that : Sij > Scrit2 => mixed draught will detrain at J<I
345	! Sij < Scrit2 => mixed draught will detrain at J>I
346
347	Scrit2 = min(Scrit(il), Sx(il))*max(0., -signhpmh(il)) + &
348	Scrit(il)*max(0., signhpmh(il))
349
350	Scrit(il) = Scrit2
351
352	!JYG Correction pour la nouvelle logique; la correction pour ALT
353	! est un peu au hazard
354	IF (Scrit(il)<=0.0) Scrit(il) = 0.0
355	IF (alt<=0.0) Scrit(il) = 1.0
356
357	smax(il) = 0.0
358	ASij(il) = 0.0
359	sup(il) = 0. ! upper S-value reached by descending draughts
360	END IF
361	END DO
362
363	! ---------------------------------------------------------------
364	DO j = minorig, nl !Loop on destination level "j"
365	! ---------------------------------------------------------------
366
367	num2 = 0
368	DO il = 1, ncum
369	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
370	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
371	lwork(il)) num2 = num2 + 1
372	END DO
373	IF (num2<=0) GO TO 175
374
375	! -----------------------------------------------
376	IF (j>i) THEN
377	! -----------------------------------------------
378	DO il = 1, ncum
379	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
380	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
381	lwork(il)) THEN
382	IF (Sij(il,i,j)>0.0) THEN
383	Smid(il) = min(Sij(il,i,j), Scrit(il))
384	Sjmax(il) = Smid(il)
385	Sjmin(il) = Smid(il)
386	IF (Smid(il)<smin(il) .AND. Sij(il,i,j+1)<Smid(il)) THEN
387	smin(il) = Smid(il)
388	Sjmax(il) = min((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j), Scrit(il))
389	Sjmin(il) = max((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
390	Sjmin(il) = min(Sjmin(il), Scrit(il))
391	Sbef(il) = Sij(il, i, j)
392	END IF
393	END IF
394	END IF
395	END DO
396	! -----------------------------------------------
397	ELSE IF (j==i) THEN
398	! -----------------------------------------------
399	DO il = 1, ncum
400	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
401	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
402	lwork(il)) THEN
403	IF (Sij(il,i,j)>0.0) THEN
404	Smid(il) = 1.
405	Sjmin(il) = max((Sij(il,i,j-1)+Smid(il))/2., Scrit(il))*max(0., -signhpmh(il)) + &
406	min((Sij(il,i,j+1)+Smid(il))/2., Scrit(il))*max(0., signhpmh(il))
407	Sjmin(il) = max(Sjmin(il), sup(il))
408	Sjmax(il) = 1.
409
410	! - preparation des variables Scrit, Smin et Sbef pour la partie j>i
411	Scrit(il) = min(Sjmin(il), Sjmax(il), Scrit(il))
412
413	smin(il) = 1.
414	Sbef(il) = max(0., signhpmh(il))
415	supmax(il, i) = sign(Scrit(il), -signhpmh(il))
416	END IF
417	END IF
418	END DO
419	! -----------------------------------------------
420	ELSE IF (j<i) THEN
421	! -----------------------------------------------
422	DO il = 1, ncum
423	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
424	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
425	lwork(il)) THEN
426	IF (Sij(il,i,j)>0.0) THEN
427	Smid(il) = max(Sij(il,i,j), Scrit(il))
428	Sjmax(il) = Smid(il)
429	Sjmin(il) = Smid(il)
430	IF (Smid(il)>smax(il) .AND. Sij(il,i,j+1)>Smid(il)) THEN
431	smax(il) = Smid(il)
432	Sjmax(il) = max((Sij(il,i,j+1)+Sij(il,i,j))/2., Sij(il,i,j))
433	Sjmax(il) = max(Sjmax(il), Scrit(il))
434	Sjmin(il) = min((Sbef(il)+Sij(il,i,j))/2., Sij(il,i,j))
435	Sjmin(il) = max(Sjmin(il), Scrit(il))
436	Sbef(il) = Sij(il, i, j)
437	END IF
438	IF (abs(Sjmin(il)-Sjmax(il))>1.E-10) &
439	sup(il) = max(Sjmin(il), Sjmax(il), sup(il))
440	END IF
441	END IF
442	END DO
443	! -----------------------------------------------
444	END IF
445	! -----------------------------------------------
446
447
448	DO il = 1, ncum
449	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
450	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
451	lwork(il)) THEN
452	IF (Sij(il,i,j)>0.0) THEN
453	rti = qnk(il) - ep(il, i)*clw(il, i)
454	Qmixmax(il) = Qmix(Sjmax(il))
455	Qmixmin(il) = Qmix(Sjmin(il))
456	Rmixmax(il) = Rmix(Sjmax(il))
457	Rmixmin(il) = Rmix(Sjmin(il))
458	sqmrmax(il) = Sjmax(il)*Qmix(Sjmax(il)) - Rmix(Sjmax(il))
459	sqmrmin(il) = Sjmin(il)*Qmix(Sjmin(il)) - Rmix(Sjmin(il))
460
461	Ment(il, i, j) = abs(Qmixmax(il)-Qmixmin(il))*Ment(il, i, j)
462
463	! Sigij(i,j) is the 'true' mixing fraction of mixture Ment(i,j)
464	IF (abs(Qmixmax(il)-Qmixmin(il))>1.E-10) THEN
465	Sigij(il, i, j) = (sqmrmax(il)-sqmrmin(il))/(Qmixmax(il)-Qmixmin(il))
466	ELSE
467	Sigij(il, i, j) = 0.
468	END IF
469
470	! -- Compute Qent, uent, vent according to the true mixing fraction
471	Qent(il, i, j) = (1.-Sigij(il,i,j))rti + Sigij(il, i, j)rr(il, i)
472	uent(il, i, j) = (1.-Sigij(il,i,j))unk(il) + Sigij(il, i, j)u(il, i)
473	vent(il, i, j) = (1.-Sigij(il,i,j))vnk(il) + Sigij(il, i, j)v(il, i)
474
475	! -- Compute liquid water static energy of mixed draughts
476	! IF (j .GT. i) THEN
477	! awat=elij(il,i,j)-(1.-ep(il,j))*clw(il,j)
478	! awat=amax1(awat,0.0)
479	! ELSE
480	! awat = 0.
481	! ENDIF
482	! Hent(il,i,j) = (1.-Sigij(il,i,j))*HP(il,i)
483	! : + Sigij(il,i,j)*H(il,i)
484	! : + (LV(il,j)+(cpd-cpv)t(il,j))awat
485	!IM 301008 beg
486	hent(il, i, j) = (1.-Sigij(il,i,j))hp(il, i) + Sigij(il, i, j)h(il, i)
487
488	!jyg<
489	! elij(il, i, j) = Qent(il, i, j) - rs(il, j)
490	! elij(il, i, j) = elij(il, i, j) + &
491	! ((h(il,j)-hent(il,i,j))rs(il,j)lv(il,j) / &
492	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
493	! elij(il, i, j) = elij(il, i, j) / &
494	! (1.+lv(il,j)lv(il,j)rs(il,j) / &
495	! ((cpd(1.-Qent(il,i,j))+Qent(il,i,j)cpv)rrvt(il,j)*t(il,j)))
496	!
497	! Computation of condensate amount Elij, taking into account the ice fraction frac
498	! Warning : the same saturation humidity rs is used over both liquid water and ice; this
499	! should be corrected.
500	!
501	! Heat capacity of mixed draught
502	cpm = cpd+Qent(il,i,j)*(cpv-cpd)
503	!
504	IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
505	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
506	elij(il, i, j) = elij(il, i, j) + &
507	(h(il,j)-hent(il,i,j)+(cpv-cpd)(Qent(il,i,j)-rr(il,j))t(il,j))* &
508	rs(il,j)lv(il,j) / (cpmrrvt(il,j)t(il,j))
509	elij(il, i, j) = elij(il, i, j) / &
510	(1.+(lv(il,j)+frac(il,j)lf(il,j))lv(il,j)*rs(il,j) / &
511	(cpmrrvt(il,j)*t(il,j)))
512	ELSE
513	elij(il, i, j) = Qent(il, i, j) - rs(il, j)
514	elij(il, i, j) = elij(il, i, j) + &
515	(h(il,j)-hent(il,i,j)+(cpv-cpd)(Qent(il,i,j)-rr(il,j))t(il,j))* &
516	rs(il,j)lv(il,j) / (cpmrrvt(il,j)t(il,j))
517	elij(il, i, j) = elij(il, i, j) / &
518	(1.+lv(il,j)lv(il,j)rs(il,j) / &
519	(cpmrrvt(il,j)*t(il,j)))
520	ENDIF
521	!>jyg
522	elij(il, i, j) = max(elij(il,i,j), 0.)
523
524	elij(il, i, j) = min(elij(il,i,j), Qent(il,i,j))
525
526	IF (j>i) THEN
527	awat = elij(il, i, j) - (1.-ep(il,j))*clw(il, j)
528	awat = amax1(awat, 0.0)
529	ELSE
530	awat = 0.
531	END IF
532
533	! print ,h(il,j)-hent(il,i,j),LV(il,j)rs(il,j)/(cpdrrvt(il,j)*
534	! : t(il,j))
535
536	!jyg<
537	! hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)t(il,j))awat
538	! Mixed draught temperature at level j
539	IF (cvflag_ice .and. frac(il,j) .gt. 0.) THEN
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)+frac(il,j)lf(il,j)+(cpd-cpv)Tm)*awat
542	ELSE
543	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))
544	hent(il, i, j) = hent(il, i, j) + (lv(il,j)+(cpd-cpv)Tm)awat
545	ENDIF
546	!>jyg
547
548	!IM 301008 end
549
550	! print *,'mix : i,j,hent(il,i,j),Sigij(il,i,j) ',
551	! : i,j,hent(il,i,j),Sigij(il,i,j)
552
553	! -- ASij is the integral of P(F) over the relevant F interval
554	ASij(il) = ASij(il) + abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
555	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
556
557	END IF
558	END IF
559	END DO
560	!jyg! DO k = 1, ntra
561	!jyg! DO il = 1, ncum
562	!jyg! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
563	!jyg! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
564	!jyg! lwork(il)) THEN
565	!jyg! IF (Sij(il,i,j)>0.0) THEN
566	!jyg! traent(il, i, j, k) = Sigij(il, i, j)*tra(il, i, k) + &
567	!jyg! (1.-Sigij(il,i,j))*tra(il, nk(il), k)
568	!jyg! END IF
569	!jyg! END IF
570	!jyg! END DO
571	!jyg! END DO
572
573	! -- If I=J (detrainement and entrainement at the same level), then only the
574	! -- adiabatic ascent part of the mixture is considered
575	IF (i==j) THEN
576	DO il = 1, ncum
577	IF (i>=icb(il) .AND. i<=inb(il) .AND. &
578	j>=(icb(il)-1) .AND. j<=inb(il) .AND. &
579	lwork(il)) THEN
580	IF (Sij(il,i,j)>0.0) THEN
581	rti = qnk(il) - ep(il, i)*clw(il, i)
582	!!! Ment(il,i,i) = m(il,i)abs(Qmixmax(il)(1.-Sjmax(il))
583	Ment(il, i, i) = abs(Qmixmax(il)*(1.-Sjmax(il))+Rmixmax(il) - &
584	Qmixmin(il)*(1.-Sjmin(il))-Rmixmin(il))
585	Qent(il, i, i) = rti
586	uent(il, i, i) = unk(il)
587	vent(il, i, i) = vnk(il)
588	hent(il, i, i) = hp(il, i)
589	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
590	Sigij(il, i, i) = 0.
591	END IF
592	END IF
593	END DO
594	!jyg! DO k = 1, ntra
595	!jyg! DO il = 1, ncum
596	!jyg! IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. &
597	!jyg! (j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. &
598	!jyg! lwork(il)) THEN
599	!jyg! IF (Sij(il,i,j)>0.0) THEN
600	!jyg! traent(il, i, i, k) = tra(il, nk(il), k)
601	!jyg! END IF
602	!jyg! END IF
603	!jyg! END DO
604	!jyg! END DO
605
606	END IF
607
608	! ---------------------------------------------------------------
609	175 END DO ! End loop on destination level "j"
610	! ---------------------------------------------------------------
611
612	DO il = 1, ncum
613	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN
614	ASij(il) = amax1(1.0E-16, ASij(il))
615	!jyg+lluis<
616	!! ASij(il) = 1.0/ASij(il)
617	ASij_inv(il) = 1.0/ASij(il)
618	! IF the F-interval spanned by possible mixtures is less than 0.01, no mixing occurs
619	IF (ASij_inv(il) > 100.) ASij_inv(il) = 0.
620	!>jyg+lluis
621	csum(il, i) = 0.0
622	END IF
623	END DO
624
625	DO j = minorig, nl
626	DO il = 1, ncum
627	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
628	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
629	!jyg Ment(il, i, j) = Ment(il, i, j)*ASij(il)
630	Ment(il, i, j) = Ment(il, i, j)*ASij_inv(il)
631	END IF
632	END DO
633	END DO
634
635	DO j = minorig, nl
636	DO il = 1, ncum
637	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. &
638	j>=(icb(il)-1) .AND. j<=inb(il)) THEN
639	csum(il, i) = csum(il, i) + Ment(il, i, j)
640	END IF
641	END DO
642	END DO
643
644	DO il = 1, ncum
645	IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
646	! cc : .and. csum(il,i).lt.m(il,i) ) then
647	nent(il, i) = 0
648	! cc Ment(il,i,i)=m(il,i)
649	Ment(il, i, i) = 1.
650	Qent(il, i, i) = qnk(il) - ep(il, i)*clw(il, i)
651	uent(il, i, i) = unk(il)
652	vent(il, i, i) = vnk(il)
653	elij(il, i, i) = clw(il, i)*(1.-ep(il,i))
654	Sij(il, i, i) = 0.0
655	END IF
656	END DO ! il
657
658	!jyg! DO j = 1, ntra
659	!jyg! DO il = 1, ncum
660	!jyg! IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. csum(il,i)<1.) THEN
661	!jyg!! cc : .and. csum(il,i).lt.m(il,i) ) then
662	!jyg! traent(il, i, i, j) = tra(il, nk(il), j)
663	!jyg! END IF
664	!jyg! END DO
665	!jyg! END DO
666
667	! ---------------------------------------------------------------
668	789 END DO ! End loop on origin level "i"
669	! ---------------------------------------------------------------
670
671
672	RETURN
673	END SUBROUTINE cv3p_mixing
674

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

Download in other formats:

Original Format