Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_ncloud5.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
Property svn:executable set to ``*
File size: 37.8 KB

Line
1	MODULE module_mp_ncloud5
2
3	REAL, PARAMETER, PRIVATE :: dtcldcr = 240.
4	INTEGER, PARAMETER, PRIVATE :: mstepmax = 100
5
6	REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
7	REAL, PARAMETER, PRIVATE :: n0r = 8.e6
8	REAL, PARAMETER, PRIVATE :: avtr = 841.9
9	REAL, PARAMETER, PRIVATE :: bvtr = 0.8
10	REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
11	REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
12	REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
13
14	REAL, PARAMETER, PRIVATE :: avts = 16.2
15	REAL, PARAMETER, PRIVATE :: bvts = .527
16	REAL, PARAMETER, PRIVATE :: xncmax = 1.e8
17	REAL, PARAMETER, PRIVATE :: n0smax = 1.e9
18	REAL, PARAMETER, PRIVATE :: betai = .6
19	REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
20	REAL, PARAMETER, PRIVATE :: dicon = 16.3
21	REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6*.8
22	REAL, PARAMETER, PRIVATE :: dimax = 400.e-6
23	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
24	REAL, PARAMETER, PRIVATE :: alpha = 1./8.18 ! .122 exponen factor for n0s
25	REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
26	REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
27	! REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e15
28	REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e5
29	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-6
30
31	REAL, PARAMETER, PRIVATE :: t40c = 233.16
32	REAL, PARAMETER, PRIVATE :: eacrc = 1.0
33
34	REAL, SAVE :: &
35	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
36	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
37	precr1,precr2,xm0,xmmax,bvts1, &
38	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
39	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
40	pidn0s,xlv1,pacrc
41
42	CONTAINS
43
44	!===================================================================
45	!
46	SUBROUTINE ncloud5(th, q, qc, qr, qi, qs &
47	,w, den, pii, p, delz &
48	,delt,g, cpd, cpv, rd, rv, t0c &
49	,ep1, ep2, qmin &
50	,XLS, XLV0, XLF0, den0, denr &
51	,cliq,cice,psat &
52	,rain, rainncv &
53	,ids,ide, jds,jde, kds,kde &
54	,ims,ime, jms,jme, kms,kme &
55	,its,ite, jts,jte, kts,kte &
56	)
57	!-------------------------------------------------------------------
58	IMPLICIT NONE
59	!-------------------------------------------------------------------
60	!
61	!Coded by Song-You Hong (NCEP) and implemented by Shuhua Chen (NCAR)
62	!
63	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
64	ims,ime, jms,jme, kms,kme , &
65	its,ite, jts,jte, kts,kte
66
67	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
68	INTENT(INOUT) :: &
69	th, &
70	q, &
71	qc, &
72	qi, &
73	qr, &
74	qs
75
76	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
77	INTENT(IN ) :: w, &
78	den, &
79	pii, &
80	p, &
81	delz
82
83
84	REAL, INTENT(IN ) :: delt, &
85	g, &
86	rd, &
87	rv, &
88	t0c, &
89	den0, &
90	cpd, &
91	cpv, &
92	ep1, &
93	ep2, &
94	qmin, &
95	xls, &
96	xlv0, &
97	xlf0, &
98	cliq, &
99	cice, &
100	psat, &
101	denr
102	REAL, DIMENSION( ims:ime , jms:jme ), &
103	INTENT(INOUT) :: rain, &
104	rainncv
105
106	! LOCAL VAR
107
108	REAL, DIMENSION( its:ite , kts:kte ) :: t
109	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
110	INTEGER :: i,j,k
111
112	!-------------------------------------------------------------------
113	DO J=jts,jte
114
115	DO k=kts,kte
116	DO i=its,ite
117	t(i,k)=th(i,k,j)*pii(i,k,j)
118	qci(i,k,1) = qc(i,k,j)
119	qci(i,k,2) = qi(i,k,j)
120	qrs(i,k,1) = qr(i,k,j)
121	qrs(i,k,2) = qs(i,k,j)
122	ENDDO
123	ENDDO
124
125	CALL ncloud52D(t, q(ims,kms,j), qci, qrs &
126	,w(ims,kms,j), den(ims,kms,j) &
127	,p(ims,kms,j), delz(ims,kms,j) &
128	,delt,g, cpd, cpv, rd, rv, t0c &
129	,ep1, ep2, qmin &
130	,XLS, XLV0, XLF0, den0, denr &
131	,cliq,cice,psat &
132	,j &
133	,rain(ims,j),rainncv(ims,j) &
134	,ids,ide, jds,jde, kds,kde &
135	,ims,ime, jms,jme, kms,kme &
136	,its,ite, jts,jte, kts,kte &
137	)
138
139	DO K=kts,kte
140	DO I=its,ite
141	th(i,k,j)=t(i,k)/pii(i,k,j)
142	qc(i,k,j) = qci(i,k,1)
143	qi(i,k,j) = qci(i,k,2)
144	qr(i,k,j) = qrs(i,k,1)
145	qs(i,k,j) = qrs(i,k,2)
146	ENDDO
147	ENDDO
148
149	ENDDO
150
151	END SUBROUTINE ncloud5
152
153	!===================================================================
154	!
155	SUBROUTINE ncloud52D(t, q, qci, qrs,w, den, p, delz &
156	,delt,g, cpd, cpv, rd, rv, t0c &
157	,ep1, ep2, qmin &
158	,XLS, XLV0, XLF0, den0, denr &
159	,cliq,cice,psat &
160	,lat &
161	,rain,rainncv &
162	,ids,ide, jds,jde, kds,kde &
163	,ims,ime, jms,jme, kms,kme &
164	,its,ite, jts,jte, kts,kte &
165	)
166	!-------------------------------------------------------------------
167	IMPLICIT NONE
168	!-------------------------------------------------------------------
169	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
170	ims,ime, jms,jme, kms,kme , &
171	its,ite, jts,jte, kts,kte, &
172	lat
173
174	REAL, DIMENSION( its:ite , kts:kte ), &
175	INTENT(INOUT) :: &
176	t
177	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
178	INTENT(INOUT) :: &
179	qci, qrs
180
181	REAL, DIMENSION( ims:ime , kms:kme ), &
182	INTENT(INOUT) :: &
183	q
184
185	REAL, DIMENSION( ims:ime , kms:kme ), &
186	INTENT(IN ) :: p, &
187	w, &
188	den, &
189	delz
190
191	REAL, INTENT(IN ) :: delt, &
192	g, &
193	cpd, &
194	cpv, &
195	t0c, &
196	den0, &
197	rd, &
198	rv, &
199	ep1, &
200	ep2, &
201	qmin, &
202	xls, &
203	xlv0, &
204	xlf0, &
205	cliq, &
206	cice, &
207	psat, &
208	denr
209	REAL, DIMENSION( ims:ime ), &
210	INTENT(INOUT) :: rain, &
211	rainncv
212
213	! LOCAL VAR
214
215	INTEGER, PARAMETER :: iun = 84
216
217	REAL, DIMENSION( its:ite , kts:kte , 2) :: &
218	rh, qs, slope, slope2, slopeb, &
219	paut, pres, falk, fall, work1
220	REAL, DIMENSION( its:ite , kts:kte ) :: &
221	falkc, work1c, work2c, fallc
222	REAL, DIMENSION( its:ite , kts:kte, 3 ) :: &
223	pacr
224	REAL, DIMENSION( its:ite , kts:kte ) :: &
225	pgen, pisd, pcon, xl, cpm, work2, q1, t1, denfac, &
226	pgens, pauts, pacrss, pisds, press, pcons, psml, psev
227
228	INTEGER, DIMENSION( its:ite ) :: mstep
229	LOGICAL, DIMENSION( its:ite ) :: flgcld
230
231	REAL :: n0sfac, pi, &
232	cpmcal, xlcal, lamdar, lamdas, diffus, &
233	viscos, xka, venfac, conden, diffac, &
234	x, y, z, a, b, c, d, e, &
235	qdt, holdrr, holdrs, supcol, &
236	coeres, supsat, dtcld, xmi, eacrs, satdt, xnc, &
237	fallsum, xlwork2, factor, source, value, &
238	xlf,pfrzdtc,pfrzdtr
239
240	INTEGER :: i,j,k, &
241	iprt, latd, lond, loop, loops, ifsat, n, numdt
242
243	!
244	!=================================================================
245	! compute internal functions
246	!
247	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
248	xlcal(x) = xlv0-xlv1*(x-t0c)
249	! tvcal(x,y) = x+xep1max(y,qmin)
250	!----------------------------------------------------------------
251	! size distributions: (x=mixing ratio, y=air density):
252	! valid for mixing ratio > 1.e-30 kg/kg.
253	! otherwise use uniform distribution value (1.e15)
254	!
255	lamdar(x,y)=(pidn0r/(xy))*.25
256	lamdas(x,y,z)=(pidn0sz/(xy))**.25
257	!
258	!----------------------------------------------------------------
259	! diffus: diffusion coefficient of the water vapor
260	! viscos: kinematic viscosity(m2s-1)
261	!
262	diffus(x,y) = 8.794e-5x*1.81/y
263	viscos(x,y) = 1.496e-6x*1.5/(x+120.)/y
264	xka(x,y) = 1.414e3viscos(x,y)y
265	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
266	venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
267	/viscos(b,c)*(.5)(den0/c)**0.25
268	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
269	!
270	pi = 4. * atan(1.)
271	!
272	!=================================================================
273	! set iprt = 0 for no unit fort.84 output
274	!
275	iprt = 1
276	if(iprt.eq.1) then
277	qdt = delt * 1000.
278	latd = (jte+jts)/2 + 1
279	lond = (ite+its)/2 + 1
280	else
281	latd = jts
282	lond = its
283	endif
284	!
285	!----------------------------------------------------------------
286	! paddint 0 for negative values generated by dynamics
287	!
288	do k = kts, kte
289	do i = its, ite
290	qci(i,k,1) = max(qci(i,k,1),0.0)
291	qrs(i,k,1) = max(qrs(i,k,1),0.0)
292	qci(i,k,2) = max(qci(i,k,2),0.0)
293	qrs(i,k,2) = max(qrs(i,k,2),0.0)
294	enddo
295	enddo
296	!
297	!----------------------------------------------------------------
298	! latent heat for phase changes and heat capacity. neglect the
299	! changes during microphysical process calculation
300	! emanuel(1994)
301	!
302	do k = kts, kte
303	do i = its, ite
304	cpm(i,k) = cpmcal(q(i,k))
305	xl(i,k) = xlcal(t(i,k))
306	enddo
307	enddo
308	if(lat.eq.latd) then
309	i = lond
310	do k = kts, kte
311	press(i,k) = 0.
312	pauts(i,k) = 0.
313	pacrss(i,k)= 0.
314	pgens(i,k) = 0.
315	pisds(i,k) = 0.
316	pcons(i,k) = 0.
317	t1(i,k) = t(i,k)
318	q1(i,k) = q(i,k)
319	enddo
320	endif
321	!
322	!----------------------------------------------------------------
323	! compute the minor time steps.
324	!
325	loops = max(nint(delt/dtcldcr),1)
326	dtcld = delt/loops
327	if(delt.le.dtcldcr) dtcld = delt
328	!
329	do loop = 1,loops
330	!
331	do i = its, ite
332	mstep(i) = 1
333	flgcld(i) = .true.
334	enddo
335	!
336	do k = kts, kte
337	do i = its, ite
338	denfac(i,k) = sqrt(den0/den(i,k))
339	enddo
340	enddo
341	!
342	do k = kts, kte
343	do i = its, ite
344	qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
345	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
346	qs(i,k,1) = max(qs(i,k,1),qmin)
347	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
348	qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
349	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
350	qs(i,k,2) = max(qs(i,k,2),qmin)
351	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
352	! if(lat.eq.latd.and.i.eq.lond) write(iun,700) qci(i,k,1)*1000., &
353	! qrs(i,k,1)1000.,qci(i,k,2)1000.,qrs(i,k,2)*1000.
354	!700 format(1x,'before computation', 4f10.4)
355	enddo
356	enddo
357	!
358	!
359	!----------------------------------------------------------------
360	! initialize the variables for microphysical physics
361	!
362	!
363	do k = kts, kte
364	do i = its, ite
365	pres(i,k,1) = 0.
366	pres(i,k,2) = 0.
367	paut(i,k,1) = 0.
368	paut(i,k,2) = 0.
369	pacr(i,k,1) = 0.
370	pacr(i,k,2) = 0.
371	pacr(i,k,3) = 0.
372	pgen(i,k) = 0.
373	pisd(i,k) = 0.
374	pcon(i,k) = 0.
375	psml(i,k) = 0.
376	psev(i,k) = 0.
377	falk(i,k,1) = 0.
378	falk(i,k,2) = 0.
379	fall(i,k,1) = 0.
380	fall(i,k,2) = 0.
381	fallc(i,k) = 0.
382	falkc(i,k) = 0.
383	enddo
384	enddo
385	!
386	!----------------------------------------------------------------
387	! sloper: the slope parameter of the rain(m-1)
388	! xka: thermal conductivity of air(jm-1s-1k-1)
389	! work1: the thermodynamic term in the denominator associated with
390	! heat conduction and vapor diffusion
391	! (ry88, y93, h85)
392	! work2: parameter associated with the ventilation effects(y93)
393	!
394	do k = kts, kte
395	do i = its, ite
396	if(qrs(i,k,1).le.qcrmin)then
397	slope(i,k,1) = lamdarmax
398	slopeb(i,k,1) = slope(i,k,1)**bvtr
399	else
400	slope(i,k,1) = lamdar(qrs(i,k,1),den(i,k))
401	slopeb(i,k,1) = slope(i,k,1)**bvtr
402	endif
403	if(qrs(i,k,2).le.qcrmin)then
404	slope(i,k,2) = lamdarmax
405	slopeb(i,k,2) = slope(i,k,2)**bvts
406	else
407	supcol = t0c-t(i,k)
408	n0sfac = min(exp(alpha*supcol),n0smax)
409	slope(i,k,2) = lamdas(qrs(i,k,2),den(i,k),n0sfac)
410	slopeb(i,k,2) = slope(i,k,2)**bvts
411	endif
412	slope2(i,k,1) = slope(i,k,1)*slope(i,k,1)
413	slope2(i,k,2) = slope(i,k,2)*slope(i,k,2)
414	enddo
415	enddo
416	!
417	do k = kts, kte
418	do i = its, ite
419	work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
420	work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
421	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
422	enddo
423	enddo
424	!
425	! instantaneous melting of cloud ice
426	!
427	do k = kts, kte
428	do i = its, ite
429	supcol = t0c-t(i,k)
430	xlf = xls-xl(i,k)
431	if(supcol.lt.0..and.qci(i,k,2).gt.qcrmin) then
432	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
433	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
434	qci(i,k,2) = 0.
435	! if(lat.eq.latd.and.i.eq.lond) write(iun,607) k,-supcol, &
436	! qci(i,k,2)1000.,xlf/cpm(i,k)qci(i,k,2)
437	endif
438	!607 format(1x,'k = ',i3,' t = ',f5.1,' qi melting = ',f10.6, &
439	! ' del t = ',f10.6)
440	!
441	! homogeneous freezing of cloud water below -40c
442	!
443	if(supcol.gt.40..and.qci(i,k,1).gt.qcrmin) then
444	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
445	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
446	qci(i,k,1) = 0.
447	! if(lat.eq.latd.and.i.eq.lond) write(iun,608) k,-supcol, &
448	! qci(i,k,1)1000.,xlf/cpm(i,k)qci(i,k,1)
449	endif
450	!608 format(1x,'k = ',i3,' t = ',f5.1,' qc home freezing = ',f10.6, &
451	! ' del t = ',f10.6)
452	!
453	! heterogeneous freezing of cloud water
454	!
455	if(supcol.gt.0..and.qci(i,k,1).gt.qcrmin) then
456	pfrzdtc = min(pfrz1exp(pfrz2supcol-1.) &
457	den(i,k)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
458	qci(i,k,2) = qci(i,k,2) + pfrzdtc
459	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
460	qci(i,k,1) = qci(i,k,1)-pfrzdtc
461	! if(lat.eq.latd.and.i.eq.lond) write(iun,609) k,-supcol, &
462	! pfrzdtc1000.,xlf/cpm(i,k)pfrzdtc
463	endif
464	!609 format(1x,'k = ',i3,' t = ',f5.1,' qc hete freezing = ',f10.6, &
465	! ' del t = ',f10.6)
466	!
467	! freezing of rain water
468	!
469	if(supcol.gt.0..and.qrs(i,k,1).gt.qcrmin) then
470	pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k) &
471	exp(pfrz2supcol-1.)slope(i,k,1)(-7)dtcld, &
472	qrs(i,k,1))
473	qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
474	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
475	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
476	! if(lat.eq.latd.and.i.eq.lond) write(iun,610) k,-supcol, &
477	! pfrzdtr1000.,xlf/cpm(i,k)pfrzdtr
478	!610 format(1x,'k = ',i3,' t = ',f5.1,' qr freezing = ',f10.6, &
479	! ' del t = ',f10.6)
480	endif
481	enddo
482	enddo
483	!
484	! warm rain process
485	! paut: auto conversion rate from cloud to rain (kgkg-1s-1)
486	! pacr: accretion rate of rain by cloud(lin83)
487	! pres: evaporation/condensation rate of rain(rh83)
488	!
489	do k = kts, kte
490	do i = its, ite
491	supsat = max(q(i,k),qmin)-qs(i,k,1)
492	satdt = supsat/dtcld
493	if(qci(i,k,1).gt.qc0) then
494	paut(i,k,1) = qck1qci(i,k,1)*(7./3.)
495	paut(i,k,1) = min(paut(i,k,1),qci(i,k,1)/dtcld)
496	endif
497	if(qrs(i,k,1).gt.qcrmin) then
498	if(qci(i,k,1).gt.qcrmin) pacr(i,k,1) &
499	= min(pacrr/slope2(i,k,1)/slope(i,k,1)/slopeb(i,k,1) &
500	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
501	coeres = slope2(i,k,1)sqrt(slope(i,k,1)slopeb(i,k,1))
502	pres(i,k,1) = (rh(i,k,1)-1.)*(precr1/slope2(i,k,1) &
503	+precr2*work2(i,k)/coeres)/work1(i,k,1)
504	if(pres(i,k,1).lt.0.) then
505	pres(i,k,1) = max(pres(i,k,1),-qrs(i,k,1)/dtcld)
506	pres(i,k,1) = max(pres(i,k,1),satdt/2)
507	else
508	pres(i,k,1) = min(pres(i,k,1),satdt/2)
509	pres(i,k,1) = min(pres(i,k,1),qrs(i,k,1)/dtcld)
510	endif
511	endif
512	enddo
513	enddo
514	!
515	!----------------------------------------------------------------
516	! cold rain process
517	! paut: conversion(aggregation) of ice to snow(kgkg-1s-1)(rh83)
518	! pgen: generation(nucleation) of ice from vapor(kgkg-1s-1)(rh83)
519	! pacr: accretion rate of snow by ice(lin83)
520	! pisd: deposition/sublimation rate of ice(rh83)
521	! pres: deposition/sublimation rate of snow(lin83)
522	!
523	do k = kts, kte
524	do i = its, ite
525	supcol = t0c-t(i,k)
526	supsat = max(q(i,k),qmin)-qs(i,k,2)
527	satdt = supsat/dtcld
528	ifsat = 0
529	n0sfac = min(exp(alpha*supcol),n0smax)
530	xnc = min(xn0 * exp(betai*supcol)/den(i,k),xncmax)
531	!
532	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qcrmin) then
533	eacrs = exp(0.025*(-supcol))
534	pacr(i,k,2) = min(pacrsn0sfaceacrs/slope2(i,k,2)/slope(i,k,2) &
535	/slopeb(i,k,2)qci(i,k,2)denfac(i,k),qci(i,k,2)/dtcld)
536	endif
537	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qcrmin) then
538	pacr(i,k,3) = min(pacrc/slope2(i,k,2)/slope(i,k,2) &
539	/slopeb(i,k,2)qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
540	endif
541	!
542	if(qci(i,k,2).gt.qcrmin) then
543	xmi = qci(i,k,2)*xnc
544	if(ifsat.ne.1) pisd(i,k) = 4.diconsqrt(xmi)*den(i,k) &
545	*(rh(i,k,2)-1.)/work1(i,k,2)
546	if(pisd(i,k).lt.0.) then
547	pisd(i,k) = max(pisd(i,k),satdt)
548	pisd(i,k) = max(pisd(i,k),-qci(i,k,2)/dtcld)
549	else
550	pisd(i,k) = min(pisd(i,k),satdt)
551	endif
552	if(abs(pisd(i,k)).gt.abs(satdt)) ifsat = 1
553	endif
554	!
555	if(qrs(i,k,2).gt.qcrmin.and.ifsat.ne.1) then
556	coeres = slope2(i,k,2)sqrt(slope(i,k,2)slopeb(i,k,2))
557	if(ifsat.ne.1) pres(i,k,2) = (rh(i,k,2)-1.)*(precs1 &
558	/slope2(i,k,2)+precs2*work2(i,k) &
559	/coeres)/work1(i,k,2)
560	if(pres(i,k,2).lt.0.) then
561	pres(i,k,2) = max(pres(i,k,2),-qrs(i,k,2)/dtcld)
562	pres(i,k,2) = max(pres(i,k,2),satdt/2)
563	else
564	pres(i,k,2) = min(pres(i,k,2),satdt/2)
565	pres(i,k,2) = min(pres(i,k,2),qrs(i,k,2)/dtcld)
566	endif
567	if(abs(pisd(i,k)+pres(i,k,2)).ge.abs(satdt)) ifsat = 1
568	endif
569	!
570	if(supsat.gt.0.and.ifsat.ne.1) then
571	pgen(i,k) = max(0.,(xm0*xnc-max(qci(i,k,2),0.))/dtcld)
572	pgen(i,k) = min(pgen(i,k),satdt)
573	endif
574	!
575	if(qci(i,k,2).gt.qcrmin) paut(i,k,2) &
576	= max(0.,(qci(i,k,2)-xmmax*xnc)/dtcld)
577	!
578	if(t(i,k).gt.t0c) then
579	xlf = xls - xl(i,k)
580	if(qrs(i,k,2).gt.qcrmin) then
581	psml(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
582	(precs1/slope2(i,k,2)+precs2work2(i,k)/coeres)
583	psml(i,k) = min(max(psml(i,k),-qrs(i,k,2)/dtcld),0.)
584	endif
585	if(qrs(i,k,2).gt.qcrmin.and.rh(i,k,1).lt.1.) &
586	psev(i,k) = pres(i,k,2)*work1(i,k,2)/work1(i,k,1)
587	psev(i,k) = min(max(psev(i,k),-qrs(i,k,2)/dtcld),0.)
588	endif
589	enddo
590	enddo
591	!
592	!----------------------------------------------------------------
593	! check mass conservation of generation terms and feedback to the
594	! large scale
595	!
596	do k = kts, kte
597	do i = its, ite
598	if(t(i,k).le.t0c) then
599	!
600	! cloud water
601	!
602	value = max(qcrmin,qci(i,k,1))
603	source = (paut(i,k,1)+pacr(i,k,1)+pacr(i,k,3))*dtcld
604	if (source.gt.value) then
605	factor = value/source
606	paut(i,k,1) = paut(i,k,1)*factor
607	pacr(i,k,1) = pacr(i,k,1)*factor
608	pacr(i,k,3) = pacr(i,k,3)*factor
609	endif
610	!
611	! cloud ice
612	!
613	value = max(qcrmin,qci(i,k,2))
614	source = (paut(i,k,2)+pacr(i,k,2)-pgen(i,k)-pisd(i,k))*dtcld
615	if (source.gt.value) then
616	factor = value/source
617	paut(i,k,2) = paut(i,k,2)*factor
618	pacr(i,k,2) = pacr(i,k,2)*factor
619	pgen(i,k) = pgen(i,k)*factor
620	pisd(i,k) = pisd(i,k)*factor
621	endif
622	work2(i,k)=-(pres(i,k,1)+pres(i,k,2)+pgen(i,k)+pisd(i,k))
623	! update
624	q(i,k) = q(i,k)+work2(i,k)*dtcld
625	qci(i,k,1) = max(qci(i,k,1)-(paut(i,k,1)+pacr(i,k,1) &
626	+pacr(i,k,3))*dtcld,0.)
627	qrs(i,k,1) = max(qrs(i,k,1)+(paut(i,k,1)+pacr(i,k,1) &
628	+pres(i,k,1))*dtcld,0.)
629	qci(i,k,2) = max(qci(i,k,2)-(paut(i,k,2)+pacr(i,k,2) &
630	-pgen(i,k)-pisd(i,k))*dtcld,0.)
631	qrs(i,k,2) = max(qrs(i,k,2)+(pres(i,k,2)+paut(i,k,2) &
632	+pacr(i,k,2)+pacr(i,k,3))*dtcld,0.)
633	xlf = xls-xl(i,k)
634	xlwork2 = -xls*(pres(i,k,2)+pisd(i,k)+pgen(i,k)) &
635	-xl(i,k)*pres(i,k,1) &
636	-xlf*pacr(i,k,3)
637	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
638	else
639	!
640	! cloud water
641	!
642	value = max(qcrmin,qci(i,k,1))
643	source=(paut(i,k,1)+pacr(i,k,1)+pacr(i,k,3))*dtcld
644	if (source.gt.value) then
645	factor = value/source
646	paut(i,k,1) = paut(i,k,1)*factor
647	pacr(i,k,1) = pacr(i,k,1)*factor
648	pacr(i,k,3) = pacr(i,k,3)*factor
649	endif
650	!
651	! snow
652	!
653	value = max(qcrmin,qrs(i,k,2))
654	source=(-psev(i,k)-psml(i,k))*dtcld
655	if (source.gt.value) then
656	factor = value/source
657	psev(i,k) = psev(i,k)*factor
658	psml(i,k) = psml(i,k)*factor
659	endif
660	work2(i,k)=-(pres(i,k,1)+psev(i,k))
661	! update
662	q(i,k) = q(i,k)+work2(i,k)*dtcld
663	qci(i,k,1) = max(qci(i,k,1)-(paut(i,k,1)+pacr(i,k,1) &
664	+pacr(i,k,3))*dtcld,0.)
665	qrs(i,k,1) = max(qrs(i,k,1)+(paut(i,k,1)+pacr(i,k,1) &
666	+pres(i,k,1) -psml(i,k)+pacr(i,k,3))*dtcld,0.)
667	qrs(i,k,2) = max(qrs(i,k,2)+(psml(i,k)+psev(i,k))*dtcld,0.)
668	xlf = xls-xl(i,k)
669	xlwork2 = -xlf(psml(i,k))-xl(i,k)(pres(i,k,1)+psev(i,k))
670	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
671	endif
672	enddo
673	enddo
674	!
675	do k = kts, kte
676	do i = its, ite
677	qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
678	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
679	qs(i,k,1) = max(qs(i,k,1),qmin)
680	qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
681	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
682	qs(i,k,2) = max(qs(i,k,2),qmin)
683	enddo
684	enddo
685	!
686	!----------------------------------------------------------------
687	! condensational/evaporational rate of cloud water if there exists
688	! additional water vapor condensated/if evaporation of cloud water
689	! is not enough to remove subsaturation.
690	! use fall bariable for this process(pcon)
691	!
692	! if(lat.eq.latd) write(iun,603)
693	do k = kts, kte
694	do i = its, ite
695	work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
696	work2(i,k) = qci(i,k,1)+work1(i,k,1)
697	pcon(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
698	if(qci(i,k,1).gt.qcrmin.and.work1(i,k,1).lt.0.and.t(i,k).gt.t0c) &
699	pcon(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
700	q(i,k) = q(i,k)-pcon(i,k)*dtcld
701	qci(i,k,1) = max(qci(i,k,1)+pcon(i,k)*dtcld,0.)
702	t(i,k) = t(i,k)+pcon(i,k)xl(i,k)/cpm(i,k)dtcld
703	!
704	if(lat.eq.latd.and.i.eq.lond) then
705	pgens(i,k) = pgens(i,k)+pgen(i,k)
706	pcons(i,k) = pcons(i,k)+pcon(i,k)
707	pisds(i,k) = pisds(i,k)+pisd(i,k)
708	pacrss(i,k) = pacrss(i,k)+pacr(i,k,1)+pacr(i,k,2)+pacr(i,k,3)
709	press(i,k) = press(i,k)+pres(i,k,1)+pres(i,k,2)
710	pauts(i,k) = pauts(i,k)+paut(i,k,1)+paut(i,k,2)
711	! write(iun,604) k,p(i,k)/100., &
712	! t(i,k)-t0c,t(i,k)-t1(i,k),q(i,k)*1000., &
713	! (q(i,k)-q1(i,k))1000.,rh(i,k,2)100.,pgens(i,k)*qdt, &
714	! pcons(i,k)qdt,pisds(i,k)qdt,pauts(i,k)qdt,pacrss(i,k)qdt, &
715	! press(i,k)qdt,qci(i,k,1)1000.,qrs(i,k,1)*1000., &
716	! qci(i,k,2)1000.,qrs(i,k,2)1000.
717	endif
718	enddo
719	enddo
720	!603 format(1x,' k',' p', &
721	! ' t',' delt',' q',' delq',' rh', &
722	! ' pgen',' pcon',' pisd',' paut',' pacr',' pres', &
723	! ' qc',' qr',' qi',' qs')
724	!604 format(1x,i3,f6.0,4f5.1,f5.0,10f5.2)
725	!
726	!----------------------------------------------------------------
727	! compute the fallout term:
728	! first, vertical terminal velosity for minor loops
729	!
730	do k = kts, kte
731	do i = its, ite
732	if(qrs(i,k,1).le.qcrmin)then
733	slope(i,k,1) = lamdarmax
734	slopeb(i,k,1) = slope(i,k,1)**bvtr
735	else
736	slope(i,k,1) = lamdar(qrs(i,k,1),den(i,k))
737	slopeb(i,k,1) = slope(i,k,1)**bvtr
738	endif
739	if(qrs(i,k,2).le.qcrmin)then
740	slope(i,k,2) = lamdarmax
741	slopeb(i,k,2) = slope(i,k,2)**bvts
742	else
743	supcol = t0c-t(i,k)
744	n0sfac = min(exp(alpha*supcol),n0smax)
745	slope(i,k,2) = lamdas(qrs(i,k,2),den(i,k),n0sfac)
746	slopeb(i,k,2) = slope(i,k,2)**bvts
747	endif
748	slope2(i,k,1) = slope(i,k,1)*slope(i,k,1)
749	slope2(i,k,2) = slope(i,k,2)*slope(i,k,2)
750	enddo
751	enddo
752	!
753	do i = its, ite
754	do k = kte, kts, -1
755	work1(i,k,1) = pvtr/slopeb(i,k,1)*denfac(i,k)/delz(i,k)
756	work1(i,k,2) = pvts/slopeb(i,k,2)*denfac(i,k)/delz(i,k)
757	if(qrs(i,k,1).le.qcrmin) work1(i,k,1) = 0.
758	if(qrs(i,k,2).le.qcrmin) work1(i,k,2) = 0.
759	numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
760	if(qci(i,k,2).le.qmin) then
761	work2c(i,k) = 0.
762	else
763	work1c(i,k) = 3.29(den(i,k)qci(i,k,2))**0.16
764	work2c(i,k) = work1c(i,k)/delz(i,k)
765	endif
766	numdt = max(nint(work2c(i,k)*dtcld+.5),numdt)
767	if(numdt.ge.mstep(i)) mstep(i) = numdt
768	enddo
769	mstep(i) = min(mstep(i),mstepmax)
770	enddo
771	!
772	! if(lat.eq.latd) write(iun,605)
773	do n = 1,mstepmax
774	k = kte
775	do i = its, ite
776	if(n.le.mstep(i)) then
777	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
778	falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)/mstep(i)
779	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
780	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
781	holdrr = qrs(i,k,1)
782	holdrs = qrs(i,k,2)
783	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
784	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
785	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
786	fallc(i,k) = fallc(i,k)+falkc(i,k)
787	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
788	!
789	! if(lat.eq.latd.and.i.eq.lond) &
790	! write(iun,606) k,p(i,k)/100., &
791	! t(i,k)-t0c,q(i,k)1000.,rh(i,k,2)100.,w(i,k),work1(i,k,1) &
792	! delz(i,k), work1(i,k,2)delz(i,k),holdrr1000.,holdrs1000., &
793	! qrs(i,k,1)1000.,qrs(i,k,2)1000.,n
794	endif
795	enddo
796	do k = kte-1, kts, -1
797	do i = its, ite
798	if(n.le.mstep(i)) then
799	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
800	falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)/mstep(i)
801	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
802	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
803	holdrr = qrs(i,k,1)
804	holdrs = qrs(i,k,2)
805	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1) &
806	-falk(i,k+1,1)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
807	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2) &
808	-falk(i,k+1,2)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
809	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
810	fallc(i,k) = fallc(i,k)+falkc(i,k)
811	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k) &
812	-falkc(i,k+1)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
813	!
814	! if(lat.eq.latd.and.i.eq.lond) &
815	! write(iun,606) k,p(i,k)/100., &
816	! t(i,k)-t0c,q(i,k)1000.,rh(i,k,2)100.,w(i,k),work1(i,k,1) &
817	! delz(i,k), work1(i,k,2)delz(i,k),holdrr1000.,holdrs1000., &
818	! qrs(i,k,1)1000.,qrs(i,k,2)1000.,n
819	endif
820	enddo
821	enddo
822	enddo
823	!605 format(1x,' k',' p',' t',' q',' rh',' w', &
824	! ' vtr',' vts',' qri',' qsi',' qrf',' qsf',' mstep')
825	!606 format(1x,i3,f6.0,2f5.1,f5.0,f6.2,6f6.2,i5)
826	!
827	!----------------------------------------------------------------
828	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
829	!
830	do i = its, ite
831	fallsum = fall(i,1,1)+fall(i,1,2)
832	if(fallsum.gt.0.) then
833	rainncv(i) = fallsumdelz(i,1)/denrdtcld*1000.
834	rain(i) = fallsumdelz(i,1)/denrdtcld*1000. &
835	+ rain(i)
836	endif
837	enddo
838	!
839	! if(lat.eq.latd) write(iun,601) latd,lond,loop,rain(lond)
840	! 601 format(1x,' ncloud5 lat lon loop : rain(mm) ',3i6,f20.2)
841	!
842	enddo ! big loops
843
844	END SUBROUTINE ncloud52d
845	! ...................................................................
846	REAL FUNCTION rgmma(x)
847	!-------------------------------------------------------------------
848	IMPLICIT NONE
849	!-------------------------------------------------------------------
850	! rgmma function: use infinite product form
851	REAL :: euler
852	PARAMETER (euler=0.577215664901532)
853	REAL :: x, y
854	INTEGER :: i
855
856	if(x.eq.1.)then
857	rgmma=0.
858	else
859	rgmma=xexp(eulerx)
860	do i=1,10000
861	y=float(i)
862	rgmma=rgmma(1.000+x/y)exp(-x/y)
863	enddo
864	rgmma=1./rgmma
865	endif
866	end function rgmma
867	!
868	!--------------------------------------------------------------------------
869	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
870	!--------------------------------------------------------------------------
871	IMPLICIT NONE
872	!--------------------------------------------------------------------------
873	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
874	xai,xbi,ttp,tr
875	INTEGER ice
876	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
877	ttp=t0c+0.01
878	dldt=cvap-cliq
879	xa=-dldt/rv
880	xb=xa+hvap/(rv*ttp)
881	dldti=cvap-cice
882	xai=-dldti/rv
883	xbi=xai+hsub/(rv*ttp)
884	tr=ttp/t
885	if(t.lt.ttp.and.ice.eq.1) then
886	fpvs=psat(trxai)exp(xbi*(1.-tr))
887	else
888	fpvs=psat(trxa)exp(xb*(1.-tr))
889	endif
890	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
891	END FUNCTION fpvs
892
893	!-------------------------------------------------------------------
894	SUBROUTINE ncloud5init(den0,denr,dens,cl,cpv,allowed_to_read)
895	!-------------------------------------------------------------------
896	IMPLICIT NONE
897	!-------------------------------------------------------------------
898	!.... constants which may not be tunable
899
900	real, intent(in) :: den0,denr,dens,cl,cpv
901	logical, intent(in) :: allowed_to_read
902	real :: pi
903
904	pi = 4.*atan(1.)
905	xlv1 = cl-cpv
906
907	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
908	qck1 = .1049.8peaut/(xncrdenr)*(1./3.)/xmyu ! 7.03
909	bvtr1 = 1.+bvtr
910	bvtr2 = 2.5+.5*bvtr
911	bvtr3 = 3.+bvtr
912	bvtr4 = 4.+bvtr
913	g1pbr = rgmma(bvtr1)
914	g3pbr = rgmma(bvtr3)
915	g4pbr = rgmma(bvtr4) ! 17.837825
916	g5pbro2 = rgmma(bvtr2) ! 1.8273
917	pvtr = avtr*g4pbr/6.
918	eacrr = 1.0
919	pacrr = pin0ravtrg3pbr.25*eacrr
920	precr1 = 2.pin0r*.78
921	precr2 = 2.pin0r.31avtr*.5g5pbro2
922	xm0 = (di0/dicon)**2
923	xmmax = (dimax/dicon)**2
924	!
925	bvts1 = 1.+bvts
926	bvts2 = 2.5+.5*bvts
927	bvts3 = 3.+bvts
928	bvts4 = 4.+bvts
929	g1pbs = rgmma(bvts1) !.8875
930	g3pbs = rgmma(bvts3)
931	g4pbs = rgmma(bvts4) ! 12.0786
932	g5pbso2 = rgmma(bvts2)
933	pvts = avts*g4pbs/6.
934	pacrs = pin0savtsg3pbs.25
935	precs1 = 4.n0s.65
936	precs2 = 4.n0s.44avts.5g5pbso2
937	pidn0r = pidenrn0r
938	pidn0s = pidensn0s
939	!
940	pacrc = pin0savtsg3pbs.25*eacrc
941	!
942	END SUBROUTINE ncloud5init
943
944	END MODULE module_mp_ncloud5
945

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

Download in other formats:

Original Format