Context Navigation

cv3p_mixing.F90 @ 5872

Last change on this file since 5872 was 2478, checked in by fhourdin, 10 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