Context Navigation

cv3p_mixing.f90 @ 5279

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

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

Download in other formats:

Original Format