Context Navigation

cv3p_mixing.F90 @ 4131

Last change on this file since 4131 was 3496, checked in by jyg, 6 years ago

Implementation of the ejection of liquid precipitation from the adiabatic ascents.
New flags:
+cvflag_prec_eject: logical

n -> old code, y -> new code

+ejectliq: real; possible values 0. & 1.

-> no liquid precipitation is ejected
-> all liquid precipitation is ejected

+ejectice: real; any value between 0. and 1.

fraction of solid precipitation ejected at each level

Note that the adiabatic ascent mass flux decrease due to precipitation ejection is not taken into account.

Attempts to do it led to water conservation violation.

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 print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
16	USE ioipsl_getin_p_mod, ONLY: getin_p
17	USE add_phys_tend_mod, ONLY: fl_cor_ebil
18
19	IMPLICIT NONE
20
21	include "cvthermo.h"
22	include "cv3param.h"
23	include "YOMCST2.h"
24	include "cvflag.h"
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