Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm6.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 65.0 KB

Line
1	#if ( RWORDSIZE == 4 )
2	# define VREC vsrec
3	# define VSQRT vssqrt
4	#else
5	# define VREC vrec
6	# define VSQRT vsqrt
7	#endif
8
9	MODULE module_mp_wsm6
10	!
11	!
12	REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
13	REAL, PARAMETER, PRIVATE :: n0r = 8.e6
14	REAL, PARAMETER, PRIVATE :: n0g = 4.e6
15	REAL, PARAMETER, PRIVATE :: avtr = 841.9
16	REAL, PARAMETER, PRIVATE :: bvtr = 0.8
17	REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
18	REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
19	REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
20	REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
21	REAL, PARAMETER, PRIVATE :: avts = 11.72
22	REAL, PARAMETER, PRIVATE :: bvts = .41
23	REAL, PARAMETER, PRIVATE :: avtg = 330.
24	REAL, PARAMETER, PRIVATE :: bvtg = 0.8
25	REAL, PARAMETER, PRIVATE :: deng = 500.
26	REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
27	REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
28	REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
29	REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
30	REAL, PARAMETER, PRIVATE :: betai = .6
31	REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
32	REAL, PARAMETER, PRIVATE :: dicon = 11.9
33	REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
34	REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
35	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
36	REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
37	REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
38	REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
39	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
40	REAL, PARAMETER, PRIVATE :: t40c = 233.16
41	REAL, PARAMETER, PRIVATE :: eacrc = 1.0
42	REAL, PARAMETER, PRIVATE :: dens = 100.0
43	REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! pgaut
44	REAL, SAVE :: &
45	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
46	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
47	bvtr6,g6pbr, &
48	precr1,precr2,xm0,xmmax,roqimax,bvts1, &
49	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
50	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
51	pidn0s,xlv1,pacrc, &
52	bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
53	g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
54	precg1,precg2,pidn0g, &
55	rslopermax,rslopesmax,rslopegmax, &
56	rsloperbmax,rslopesbmax,rslopegbmax, &
57	rsloper2max,rslopes2max,rslopeg2max, &
58	rsloper3max,rslopes3max,rslopeg3max
59	CONTAINS
60	!===================================================================
61	!
62	SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg &
63	,den, pii, p, delz &
64	,delt,g, cpd, cpv, rd, rv, t0c &
65	,ep1, ep2, qmin &
66	,XLS, XLV0, XLF0, den0, denr &
67	,cliq,cice,psat &
68	,rain, rainncv &
69	,snow, snowncv &
70	,graupel, graupelncv &
71	,sr &
72	,ids,ide, jds,jde, kds,kde &
73	,ims,ime, jms,jme, kms,kme &
74	,its,ite, jts,jte, kts,kte &
75	)
76	!-------------------------------------------------------------------
77	IMPLICIT NONE
78	!-------------------------------------------------------------------
79	!
80	! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the WRF
81	! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
82	! number concentration is a function of temperature, and seperate assumption
83	! is developed, in which ice crystal number concentration is a function
84	! of ice amount. A theoretical background of the ice-microphysics and related
85	! processes in the WSMMPs are described in Hong et al. (2004).
86	! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
87	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
88	!
89	! WSM6 cloud scheme
90	!
91	! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
92	! Summer 2003
93	!
94	! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
95	! Summer 2004
96	!
97	! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
98	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
99	! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
100	! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
101	! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
102	! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
103	!
104	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
105	ims,ime, jms,jme, kms,kme , &
106	its,ite, jts,jte, kts,kte
107	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
108	INTENT(INOUT) :: &
109	th, &
110	q, &
111	qc, &
112	qi, &
113	qr, &
114	qs, &
115	qg
116	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
117	INTENT(IN ) :: &
118	den, &
119	pii, &
120	p, &
121	delz
122	REAL, INTENT(IN ) :: delt, &
123	g, &
124	rd, &
125	rv, &
126	t0c, &
127	den0, &
128	cpd, &
129	cpv, &
130	ep1, &
131	ep2, &
132	qmin, &
133	XLS, &
134	XLV0, &
135	XLF0, &
136	cliq, &
137	cice, &
138	psat, &
139	denr
140	REAL, DIMENSION( ims:ime , jms:jme ), &
141	INTENT(INOUT) :: rain, &
142	rainncv, &
143	sr
144
145	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
146	INTENT(INOUT) :: snow, &
147	snowncv
148
149	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
150	INTENT(INOUT) :: graupel, &
151	graupelncv
152	! LOCAL VAR
153	REAL, DIMENSION( its:ite , kts:kte ) :: t
154	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
155	REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
156	INTEGER :: i,j,k
157	!-------------------------------------------------------------------
158	DO j=jts,jte
159	DO k=kts,kte
160	DO i=its,ite
161	t(i,k)=th(i,k,j)*pii(i,k,j)
162	qci(i,k,1) = qc(i,k,j)
163	qci(i,k,2) = qi(i,k,j)
164	qrs(i,k,1) = qr(i,k,j)
165	qrs(i,k,2) = qs(i,k,j)
166	qrs(i,k,3) = qg(i,k,j)
167	ENDDO
168	ENDDO
169
170	! Sending array starting locations of optional variables may cause
171	! troubles, so we explicitly change the call.
172
173	CALL wsm62D(t, q(ims,kms,j), qci, qrs &
174	,den(ims,kms,j) &
175	,p(ims,kms,j), delz(ims,kms,j) &
176	,delt,g, cpd, cpv, rd, rv, t0c &
177	,ep1, ep2, qmin &
178	,XLS, XLV0, XLF0, den0, denr &
179	,cliq,cice,psat &
180	,j &
181	,rain(ims,j),rainncv(ims,j) &
182	,sr(ims,j) &
183	,ids,ide, jds,jde, kds,kde &
184	,ims,ime, jms,jme, kms,kme &
185	,its,ite, jts,jte, kts,kte &
186	#if ( EM_CORE == 1 )
187	,snow(ims,j),snowncv(ims,j) &
188	,graupel(ims,j),graupelncv(ims,j) &
189	#endif
190	)
191
192	DO K=kts,kte
193	DO I=its,ite
194	th(i,k,j)=t(i,k)/pii(i,k,j)
195	qc(i,k,j) = qci(i,k,1)
196	qi(i,k,j) = qci(i,k,2)
197	qr(i,k,j) = qrs(i,k,1)
198	qs(i,k,j) = qrs(i,k,2)
199	qg(i,k,j) = qrs(i,k,3)
200	ENDDO
201	ENDDO
202	ENDDO
203	END SUBROUTINE wsm6
204	!===================================================================
205	!
206	SUBROUTINE wsm62D(t, q, qci, qrs, den, p, delz &
207	,delt,g, cpd, cpv, rd, rv, t0c &
208	,ep1, ep2, qmin &
209	,XLS, XLV0, XLF0, den0, denr &
210	,cliq,cice,psat &
211	,lat &
212	,rain,rainncv &
213	,sr &
214	,ids,ide, jds,jde, kds,kde &
215	,ims,ime, jms,jme, kms,kme &
216	,its,ite, jts,jte, kts,kte &
217	,snow,snowncv &
218	,graupel,graupelncv &
219	)
220	!-------------------------------------------------------------------
221	IMPLICIT NONE
222	!-------------------------------------------------------------------
223	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
224	ims,ime, jms,jme, kms,kme , &
225	its,ite, jts,jte, kts,kte, &
226	lat
227	REAL, DIMENSION( its:ite , kts:kte ), &
228	INTENT(INOUT) :: &
229	t
230	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
231	INTENT(INOUT) :: &
232	qci
233	REAL, DIMENSION( its:ite , kts:kte, 3 ), &
234	INTENT(INOUT) :: &
235	qrs
236	REAL, DIMENSION( ims:ime , kms:kme ), &
237	INTENT(INOUT) :: &
238	q
239	REAL, DIMENSION( ims:ime , kms:kme ), &
240	INTENT(IN ) :: &
241	den, &
242	p, &
243	delz
244	REAL, INTENT(IN ) :: delt, &
245	g, &
246	cpd, &
247	cpv, &
248	t0c, &
249	den0, &
250	rd, &
251	rv, &
252	ep1, &
253	ep2, &
254	qmin, &
255	XLS, &
256	XLV0, &
257	XLF0, &
258	cliq, &
259	cice, &
260	psat, &
261	denr
262	REAL, DIMENSION( ims:ime ), &
263	INTENT(INOUT) :: rain, &
264	rainncv, &
265	sr
266	REAL, DIMENSION( ims:ime ), OPTIONAL, &
267	INTENT(INOUT) :: snow, &
268	snowncv
269
270	REAL, DIMENSION( ims:ime ), OPTIONAL, &
271	INTENT(INOUT) :: graupel, &
272	graupelncv
273	! LOCAL VAR
274	REAL, DIMENSION( its:ite , kts:kte , 3) :: &
275	rh, qs, rslope, rslope2, rslope3, rslopeb, &
276	falk, fall, work1
277	REAL, DIMENSION( its:ite , kts:kte ) :: &
278	worka
279	REAL, DIMENSION( its:ite , kts:kte ) :: &
280	falkc, work1c, work2c, fallc
281	REAL, DIMENSION( its:ite , kts:kte ) :: &
282	prevp, psdep, pgdep, praut, psaut, pgaut, &
283	pracw, psacw, pgacw, pgacr, pgacs, psaci, pgmlt, praci, &
284	piacr, pracs, psacr, pgaci, pseml, pgeml
285	REAL, DIMENSION( its:ite , kts:kte) :: qsum
286	REAL, DIMENSION( its:ite , kts:kte ) :: paacw
287	REAL, DIMENSION( its:ite , kts:kte ) :: &
288	pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, denfac, &
289	xni, pgevp,n0sfac
290	! variables for optimization
291	REAL, DIMENSION( its:ite ) :: tvec1
292	REAL :: temp
293	INTEGER, DIMENSION( its:ite ) :: mstep, numdt
294	LOGICAL, DIMENSION( its:ite ) :: flgcld
295	REAL :: pi, &
296	cpmcal, xlcal, lamdar, lamdas, lamdag, diffus, &
297	viscos, xka, venfac, conden, diffac, &
298	x, y, z, a, b, c, d, e, &
299	qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
300	coeres, supsat, dtcld, xmi, eacrs, satdt, &
301	qimax, diameter, xni0, roqi0, &
302	fallsum, fallsum_qsi, fallsum_qg, &
303	vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
304	xlwork2, factor, source, value, &
305	xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
306	REAL :: vt2ave
307	REAL :: holdc, holdci
308	INTEGER :: i, j, k, mstepmax, &
309	iprt, latd, lond, loop, loops, ifsat, n
310	! Temporaries used for inlining fpvs function
311	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
312	!
313	!=================================================================
314	! compute internal functions
315	!
316	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
317	xlcal(x) = xlv0-xlv1*(x-t0c)
318	!----------------------------------------------------------------
319	! size distributions: (x=mixing ratio, y=air density):
320	! valid for mixing ratio > 1.e-9 kg/kg.
321	!
322	! Optimizatin : A*B => exp(log(A)(B))
323	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
324	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
325	lamdag(x,y)= sqrt(sqrt(pidn0g/(xy))) ! (pidn0g/(xy))**.25
326	!
327	!----------------------------------------------------------------
328	! diffus: diffusion coefficient of the water vapor
329	! viscos: kinematic viscosity(m2s-1)
330	!
331	diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
332	viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
333	xka(x,y) = 1.414e3viscos(x,y)y
334	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
335	venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
336	/sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
337	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
338	!
339	pi = 4. * atan(1.)
340	!
341	!
342	!----------------------------------------------------------------
343	! paddint 0 for negative values generated by dynamics
344	!
345	do k = kts, kte
346	do i = its, ite
347	qci(i,k,1) = max(qci(i,k,1),0.0)
348	qrs(i,k,1) = max(qrs(i,k,1),0.0)
349	qci(i,k,2) = max(qci(i,k,2),0.0)
350	qrs(i,k,2) = max(qrs(i,k,2),0.0)
351	qrs(i,k,3) = max(qrs(i,k,3),0.0)
352	enddo
353	enddo
354	!
355	!----------------------------------------------------------------
356	! latent heat for phase changes and heat capacity. neglect the
357	! changes during microphysical process calculation
358	! emanuel(1994)
359	!
360	do k = kts, kte
361	do i = its, ite
362	cpm(i,k) = cpmcal(q(i,k))
363	xl(i,k) = xlcal(t(i,k))
364	enddo
365	enddo
366	!
367	!----------------------------------------------------------------
368	! compute the minor time steps.
369	!
370	loops = max(nint(delt/dtcldcr),1)
371	dtcld = delt/loops
372	if(delt.le.dtcldcr) dtcld = delt
373	!
374	do loop = 1,loops
375	!
376	!----------------------------------------------------------------
377	! initialize the large scale variables
378	!
379	do i = its, ite
380	mstep(i) = 1
381	flgcld(i) = .true.
382	enddo
383	!
384	! do k = kts, kte
385	! do i = its, ite
386	! denfac(i,k) = sqrt(den0/den(i,k))
387	! enddo
388	! enddo
389	do k = kts, kte
390	CALL VREC( tvec1(its), den(its,k), ite-its+1)
391	do i = its, ite
392	tvec1(i) = tvec1(i)*den0
393	enddo
394	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
395	enddo
396	!
397	! Inline expansion for fpvs
398	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
399	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
400	hsub = xls
401	hvap = xlv0
402	cvap = cpv
403	ttp=t0c+0.01
404	dldt=cvap-cliq
405	xa=-dldt/rv
406	xb=xa+hvap/(rv*ttp)
407	dldti=cvap-cice
408	xai=-dldti/rv
409	xbi=xai+hsub/(rv*ttp)
410	do k = kts, kte
411	do i = its, ite
412	tr=ttp/t(i,k)
413	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
414	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
415	qs(i,k,1) = max(qs(i,k,1),qmin)
416	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
417	tr=ttp/t(i,k)
418	if(t(i,k).lt.ttp) then
419	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
420	else
421	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
422	endif
423	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
424	qs(i,k,2) = max(qs(i,k,2),qmin)
425	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
426	enddo
427	enddo
428	!
429	!----------------------------------------------------------------
430	! initialize the variables for microphysical physics
431	!
432	!
433	do k = kts, kte
434	do i = its, ite
435	prevp(i,k) = 0.
436	psdep(i,k) = 0.
437	pgdep(i,k) = 0.
438	praut(i,k) = 0.
439	psaut(i,k) = 0.
440	pgaut(i,k) = 0.
441	pracw(i,k) = 0.
442	praci(i,k) = 0.
443	piacr(i,k) = 0.
444	psaci(i,k) = 0.
445	psacw(i,k) = 0.
446	pracs(i,k) = 0.
447	psacr(i,k) = 0.
448	pgacw(i,k) = 0.
449	paacw(i,k) = 0.
450	pgaci(i,k) = 0.
451	pgacr(i,k) = 0.
452	pgacs(i,k) = 0.
453	pigen(i,k) = 0.
454	pidep(i,k) = 0.
455	pcond(i,k) = 0.
456	psmlt(i,k) = 0.
457	pgmlt(i,k) = 0.
458	pseml(i,k) = 0.
459	pgeml(i,k) = 0.
460	psevp(i,k) = 0.
461	pgevp(i,k) = 0.
462	falk(i,k,1) = 0.
463	falk(i,k,2) = 0.
464	falk(i,k,3) = 0.
465	fall(i,k,1) = 0.
466	fall(i,k,2) = 0.
467	fall(i,k,3) = 0.
468	fallc(i,k) = 0.
469	falkc(i,k) = 0.
470	xni(i,k) = 1.e3
471	enddo
472	enddo
473	!
474	!----------------------------------------------------------------
475	! compute the fallout term:
476	! first, vertical terminal velosity for minor loops
477	!
478	do k = kts, kte
479	do i = its, ite
480	supcol = t0c-t(i,k)
481	!---------------------------------------------------------------
482	! n0s: Intercept parameter for snow [m-4] [HDC 6]
483	!---------------------------------------------------------------
484	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
485	if(qrs(i,k,1).le.qcrmin)then
486	rslope(i,k,1) = rslopermax
487	rslopeb(i,k,1) = rsloperbmax
488	rslope2(i,k,1) = rsloper2max
489	rslope3(i,k,1) = rsloper3max
490	else
491	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
492	rslopeb(i,k,1) = rslope(i,k,1)**bvtr
493	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
494	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
495	endif
496	if(qrs(i,k,2).le.qcrmin)then
497	rslope(i,k,2) = rslopesmax
498	rslopeb(i,k,2) = rslopesbmax
499	rslope2(i,k,2) = rslopes2max
500	rslope3(i,k,2) = rslopes3max
501	else
502	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
503	rslopeb(i,k,2) = rslope(i,k,2)**bvts
504	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
505	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
506	endif
507	if(qrs(i,k,3).le.qcrmin)then
508	rslope(i,k,3) = rslopegmax
509	rslopeb(i,k,3) = rslopegbmax
510	rslope2(i,k,3) = rslopeg2max
511	rslope3(i,k,3) = rslopeg3max
512	else
513	rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
514	rslopeb(i,k,3) = rslope(i,k,3)**bvtg
515	rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
516	rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
517	endif
518	!-------------------------------------------------------------
519	! Ni: ice crystal number concentraiton [HDC 5c]
520	!-------------------------------------------------------------
521	! xni(i,k) = min(max(5.38e7*(den(i,k) &
522	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
523	temp = (den(i,k)*max(qci(i,k,2),qmin))
524	temp = sqrt(sqrt(temptemptemp))
525	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
526	enddo
527	enddo
528	!
529	mstepmax = 1
530	numdt = 1
531	do k = kte, kts, -1
532	do i = its, ite
533	work1(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)/delz(i,k)
534	work1(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)/delz(i,k)
535	work1(i,k,3) = pvtgrslopeb(i,k,3)denfac(i,k)/delz(i,k)
536	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
537	IF ( qsum(i,k) .gt. 1.e-15 ) THEN
538	worka(i,k) = (work1(i,k,2)qrs(i,k,2) + work1(i,k,3)qrs(i,k,3))/qsum(i,k)
539	ELSE
540	worka(i,k) = 0.
541	ENDIF
542	numdt(i) = max(nint(max(work1(i,k,1),worka(i,k)) &
543	*dtcld+.5),1)
544	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
545	enddo
546	enddo
547	do i = its, ite
548	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
549	enddo
550	!
551	do n = 1, mstepmax
552	k = kte
553	do i = its, ite
554	if(n.le.mstep(i)) then
555	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
556	falk(i,k,2) = den(i,k)qrs(i,k,2)worka(i,k)/mstep(i)
557	falk(i,k,3) = den(i,k)qrs(i,k,3)worka(i,k)/mstep(i)
558	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
559	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
560	fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
561	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
562	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
563	qrs(i,k,3) = max(qrs(i,k,3)-falk(i,k,3)*dtcld/den(i,k),0.)
564	endif
565	enddo
566	do k = kte-1, kts, -1
567	do i = its, ite
568	if(n.le.mstep(i)) then
569	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
570	falk(i,k,2) = den(i,k)qrs(i,k,2)worka(i,k)/mstep(i)
571	falk(i,k,3) = den(i,k)qrs(i,k,3)worka(i,k)/mstep(i)
572	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
573	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
574	fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
575	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
576	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
577	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
578	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
579	qrs(i,k,3) = max(qrs(i,k,3)-(falk(i,k,3)-falk(i,k+1,3) &
580	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
581	endif
582	enddo
583	enddo
584	do k = kte, kts, -1
585	do i = its, ite
586	if(n.le.mstep(i).and.t(i,k).gt.t0c) then
587	!---------------------------------------------------------------
588	! psmlt: melting of snow [HL A33] [RH83 A25]
589	! (T>T0: S->R)
590	!---------------------------------------------------------------
591	xlf = xlf0
592	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
593	if(qrs(i,k,2).gt.0.) then
594	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
595	psmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
596	n0sfac(i,k)(precs1*rslope2(i,k,2) &
597	+precs2work2(i,k)coeres)
598	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
599	-qrs(i,k,2)/mstep(i)),0.)
600	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
601	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
602	t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
603	endif
604	!---------------------------------------------------------------
605	! pgmlt: melting of graupel [HL A23] [LFO 47]
606	! (T>T0: G->R)
607	!---------------------------------------------------------------
608	if(qrs(i,k,3).gt.0.) then
609	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
610	pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
611	(t0c-t(i,k))(precg1*rslope2(i,k,3) &
612	+precg2work2(i,k)coeres)
613	pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
614	-qrs(i,k,3)/mstep(i)),0.)
615	qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
616	qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
617	t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
618	endif
619	endif
620	enddo
621	enddo
622	enddo
623	!---------------------------------------------------------------
624	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
625	!---------------------------------------------------------------
626	mstepmax = 1
627	mstep = 1
628	numdt = 1
629	do k = kte, kts, -1
630	do i = its, ite
631	if(qci(i,k,2).le.0.) then
632	work2c(i,k) = 0.
633	else
634	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
635	! diameter = min(dicon * sqrt(xmi),dimax)
636	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
637	work1c(i,k) = 1.49e4diameter*1.31
638	work2c(i,k) = work1c(i,k)/delz(i,k)
639	endif
640	numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
641	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
642	enddo
643	enddo
644	do i = its, ite
645	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
646	enddo
647	!
648	do n = 1, mstepmax
649	k = kte
650	do i = its, ite
651	if(n.le.mstep(i)) then
652	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
653	holdc = falkc(i,k)
654	fallc(i,k) = fallc(i,k)+falkc(i,k)
655	holdci = qci(i,k,2)
656	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
657	endif
658	enddo
659	do k = kte-1, kts, -1
660	do i = its, ite
661	if(n.le.mstep(i)) then
662	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
663	holdc = falkc(i,k)
664	fallc(i,k) = fallc(i,k)+falkc(i,k)
665	holdci = qci(i,k,2)
666	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
667	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
668	endif
669	enddo
670	enddo
671	enddo
672	!
673	!----------------------------------------------------------------
674	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
675	!
676	do i = its, ite
677	fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
678	fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
679	fallsum_qg = fall(i,kts,3)
680	rainncv(i) = 0.
681	if(fallsum.gt.0.) then
682	rainncv(i) = fallsumdelz(i,kts)/denrdtcld*1000.
683	rain(i) = fallsumdelz(i,kts)/denrdtcld*1000. + rain(i)
684	endif
685	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
686	snowncv(i) = 0.
687	if(fallsum_qsi.gt.0.) then
688	snowncv(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000.
689	snow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i)
690	endif
691	ENDIF
692	IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
693	graupelncv(i) = 0.
694	if(fallsum_qg.gt.0.) then
695	graupelncv(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000.
696	graupel(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000. + graupel(i)
697	endif
698	ENDIF
699	sr(i) = 0.
700	if(fallsum.gt.0.)sr(i)=(fallsum_qsidelz(i,kts)/denrdtcld*1000. + &
701	fallsum_qgdelz(i,kts)/denrdtcld*1000.)/(rainncv(i)+1.e-12)
702	enddo
703	!
704	!---------------------------------------------------------------
705	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
706	! (T>T0: I->C)
707	!---------------------------------------------------------------
708	do k = kts, kte
709	do i = its, ite
710	supcol = t0c-t(i,k)
711	xlf = xls-xl(i,k)
712	if(supcol.lt.0.) xlf = xlf0
713	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
714	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
715	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
716	qci(i,k,2) = 0.
717	endif
718	!---------------------------------------------------------------
719	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
720	! (T<-40C: C->I)
721	!---------------------------------------------------------------
722	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
723	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
724	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
725	qci(i,k,1) = 0.
726	endif
727	!---------------------------------------------------------------
728	! pihtf: heterogeneous freezing of cloud water [HL A44]
729	! (T0>T>-40C: C->I)
730	!---------------------------------------------------------------
731	if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
732	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
733	! den(i,k)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
734	supcolt=min(supcol,50.)
735	pfrzdtc = min(pfrz1(exp(pfrz2supcolt)-1.) &
736	den(i,k)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
737	qci(i,k,2) = qci(i,k,2) + pfrzdtc
738	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
739	qci(i,k,1) = qci(i,k,1)-pfrzdtc
740	endif
741	!---------------------------------------------------------------
742	! pgfrz: freezing of rain water [HL A20] [LFO 45]
743	! (T<T0, R->G)
744	!---------------------------------------------------------------
745	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
746	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k) &
747	! (exp(pfrz2supcol)-1.)rslope3(i,k,1)*2 &
748	! rslope(i,k,1)dtcld,qrs(i,k,1))
749	temp = rslope3(i,k,1)
750	temp = temptemprslope(i,k,1)
751	supcolt=min(supcol,50.)
752	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k) &
753	(exp(pfrz2supcolt)-1.)tempdtcld, &
754	qrs(i,k,1))
755	qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
756	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
757	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
758	endif
759	enddo
760	enddo
761	!
762	!
763	!----------------------------------------------------------------
764	! rsloper: reverse of the slope parameter of the rain(m)
765	! xka: thermal conductivity of air(jm-1s-1k-1)
766	! work1: the thermodynamic term in the denominator associated with
767	! heat conduction and vapor diffusion
768	! (ry88, y93, h85)
769	! work2: parameter associated with the ventilation effects(y93)
770	!
771	do k = kts, kte
772	do i = its, ite
773	if(qrs(i,k,1).le.qcrmin)then
774	rslope(i,k,1) = rslopermax
775	rslopeb(i,k,1) = rsloperbmax
776	rslope2(i,k,1) = rsloper2max
777	rslope3(i,k,1) = rsloper3max
778	else
779	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
780	rslopeb(i,k,1) = rslope(i,k,1)**bvtr
781	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
782	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
783	endif
784	if(qrs(i,k,2).le.qcrmin)then
785	rslope(i,k,2) = rslopesmax
786	rslopeb(i,k,2) = rslopesbmax
787	rslope2(i,k,2) = rslopes2max
788	rslope3(i,k,2) = rslopes3max
789	else
790	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
791	rslopeb(i,k,2) = rslope(i,k,2)**bvts
792	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
793	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
794	endif
795	if(qrs(i,k,3).le.qcrmin)then
796	rslope(i,k,3) = rslopegmax
797	rslopeb(i,k,3) = rslopegbmax
798	rslope2(i,k,3) = rslopeg2max
799	rslope3(i,k,3) = rslopeg3max
800	else
801	rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
802	rslopeb(i,k,3) = rslope(i,k,3)**bvtg
803	rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
804	rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
805	endif
806	enddo
807	enddo
808	!
809	do k = kts, kte
810	do i = its, ite
811	work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
812	work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
813	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
814	enddo
815	enddo
816	!
817	!===============================================================
818	!
819	! warm rain processes
820	!
821	! - follows the processes in RH83 and LFO except for autoconcersion
822	!
823	!===============================================================
824	!
825	do k = kts, kte
826	do i = its, ite
827	supsat = max(q(i,k),qmin)-qs(i,k,1)
828	satdt = supsat/dtcld
829	!---------------------------------------------------------------
830	! praut: auto conversion rate from cloud to rain [HDC 16]
831	! (C->R)
832	!---------------------------------------------------------------
833	if(qci(i,k,1).gt.qc0) then
834	praut(i,k) = qck1qci(i,k,1)*(7./3.)
835	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
836	endif
837	!---------------------------------------------------------------
838	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
839	! (C->R)
840	!---------------------------------------------------------------
841	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
842	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
843	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
844	endif
845	!---------------------------------------------------------------
846	! prevp: evaporation/condensation rate of rain [HDC 14]
847	! (V->R or R->V)
848	!---------------------------------------------------------------
849	if(qrs(i,k,1).gt.0.) then
850	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
851	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
852	+precr2work2(i,k)coeres)/work1(i,k,1)
853	if(prevp(i,k).lt.0.) then
854	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
855	prevp(i,k) = max(prevp(i,k),satdt/2)
856	else
857	prevp(i,k) = min(prevp(i,k),satdt/2)
858	endif
859	endif
860	enddo
861	enddo
862	!
863	!===============================================================
864	!
865	! cold rain processes
866	!
867	! - follows the revised ice microphysics processes in HDC
868	! - the processes same as in RH83 and RH84 and LFO behave
869	! following ice crystal hapits defined in HDC, inclduing
870	! intercept parameter for snow (n0s), ice crystal number
871	! concentration (ni), ice nuclei number concentration
872	! (n0i), ice diameter (d)
873	!
874	!===============================================================
875	!
876	do k = kts, kte
877	do i = its, ite
878	supcol = t0c-t(i,k)
879	supsat = max(q(i,k),qmin)-qs(i,k,2)
880	satdt = supsat/dtcld
881	ifsat = 0
882	!-------------------------------------------------------------
883	! Ni: ice crystal number concentraiton [HDC 5c]
884	!-------------------------------------------------------------
885	! xni(i,k) = min(max(5.38e7*(den(i,k) &
886	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
887	temp = (den(i,k)*max(qci(i,k,2),qmin))
888	temp = sqrt(sqrt(temptemptemp))
889	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
890	eacrs = exp(0.07*(-supcol))
891	!
892	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
893	diameter = min(dicon * sqrt(xmi),dimax)
894	vt2i = 1.49e4diameter*1.31
895	vt2r=pvtrrslopeb(i,k,1)denfac(i,k)
896	vt2s=pvtsrslopeb(i,k,2)denfac(i,k)
897	vt2g=pvtgrslopeb(i,k,3)denfac(i,k)
898	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
899	if(qsum(i,k) .gt. 1.e-15) then
900	vt2ave=(vt2sqrs(i,k,2)+vt2gqrs(i,k,3))/(qsum(i,k))
901	else
902	vt2ave=0.
903	endif
904	if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
905	if(qrs(i,k,1).gt.qcrmin) then
906	!-------------------------------------------------------------
907	! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
908	! (T<T0: I->R)
909	!-------------------------------------------------------------
910	acrfac = 2.rslope3(i,k,1)+2.diameter*rslope2(i,k,1) &
911	+diameter*2rslope(i,k,1)
912	praci(i,k) = piqci(i,k,2)n0rabs(vt2r-vt2i)acrfac/4.
913	praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
914	!-------------------------------------------------------------
915	! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
916	! (T<T0: R->S or R->G)
917	!-------------------------------------------------------------
918	piacr(i,k) = pi*2avtrn0rdenrxni(i,k)denfac(i,k) &
919	g6pbrrslope3(i,k,1)*rslope3(i,k,1) &
920	*rslopeb(i,k,1)/24./den(i,k)
921	piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
922	endif
923	!-------------------------------------------------------------
924	! psaci: Accretion of cloud ice by snow [HDC 10]
925	! (T<T0: I->S)
926	!-------------------------------------------------------------
927	if(qrs(i,k,2).gt.qcrmin) then
928	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
929	+diameter*2rslope(i,k,2)
930	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
931	abs(vt2ave-vt2i)acrfac/4.
932	psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
933	endif
934	!-------------------------------------------------------------
935	! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
936	! (T<T0: I->G)
937	!-------------------------------------------------------------
938	if(qrs(i,k,3).gt.qcrmin) then
939	egi = exp(0.07*(-supcol))
940	acrfac = 2.rslope3(i,k,3)+2.diameter*rslope2(i,k,3) &
941	+diameter*2rslope(i,k,3)
942	pgaci(i,k) = piegiqci(i,k,2)n0gabs(vt2ave-vt2i)*acrfac/4.
943	pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
944	endif
945	endif
946	!-------------------------------------------------------------
947	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
948	! (T<T0: C->S, and T>=T0: C->R)
949	!-------------------------------------------------------------
950	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
951	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2) &
952	rslopeb(i,k,2)qci(i,k,1)*denfac(i,k) &
953	,qci(i,k,1)/dtcld)
954	endif
955	!-------------------------------------------------------------
956	! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
957	! (T<T0: C->G, and T>=T0: C->R)
958	!-------------------------------------------------------------
959	if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
960	pgacw(i,k) = min(pacrgrslope3(i,k,3)rslopeb(i,k,3) &
961	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
962	endif
963	!-------------------------------------------------------------
964	! paacw: Accretion of cloud water by averaged snow/graupel
965	! (T<T0: C->G or S, and T>=T0: C->R)
966	!-------------------------------------------------------------
967	if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,3).gt.qcrmin) then
968	paacw(i,k) = (qrs(i,k,2)psacw(i,k)+qrs(i,k,3)pgacw(i,k))/(qsum(i,k))
969	endif
970	!-------------------------------------------------------------
971	! pracs: Accretion of snow by rain [HL A11] [LFO 27]
972	! (T<T0: S->G)
973	!-------------------------------------------------------------
974	if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
975	if(supcol.gt.0) then
976	acrfac = 5.rslope3(i,k,2)rslope3(i,k,2)*rslope(i,k,1) &
977	+2.rslope3(i,k,2)rslope2(i,k,2)*rslope2(i,k,1) &
978	+.5rslope2(i,k,2)rslope2(i,k,2)*rslope3(i,k,1)
979	pracs(i,k) = pi*2n0rn0sn0sfac(i,k)*abs(vt2r-vt2ave) &
980	(dens/den(i,k))acrfac
981	pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
982	endif
983	!-------------------------------------------------------------
984	! psacr: Accretion of rain by snow [HL A10] [LFO 28]
985	! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
986	!-------------------------------------------------------------
987	acrfac = 5.rslope3(i,k,1)rslope3(i,k,1)*rslope(i,k,2) &
988	+2.rslope3(i,k,1)rslope2(i,k,1)*rslope2(i,k,2) &
989	+.5rslope2(i,k,1)rslope2(i,k,1)*rslope3(i,k,2)
990	psacr(i,k) = pi*2n0rn0sn0sfac(i,k)*abs(vt2ave-vt2r) &
991	(denr/den(i,k))acrfac
992	psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
993	endif
994	!-------------------------------------------------------------
995	! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
996	! (T<T0: R->G) (T>=T0: enhance melting of graupel)
997	!-------------------------------------------------------------
998	if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
999	acrfac = 5.rslope3(i,k,1)rslope3(i,k,1)*rslope(i,k,3) &
1000	+2.rslope3(i,k,1)rslope2(i,k,1)*rslope2(i,k,3) &
1001	+.5rslope2(i,k,1)rslope2(i,k,1)*rslope3(i,k,3)
1002	pgacr(i,k) = pi*2n0rn0gabs(vt2ave-vt2r)*(denr/den(i,k)) &
1003	*acrfac
1004	pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
1005	endif
1006	!
1007	!-------------------------------------------------------------
1008	! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
1009	! (S->G): This process is eliminated in V3.0 with the
1010	! new combined snow/graupel fall speeds
1011	!-------------------------------------------------------------
1012	! if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
1013	! acrfac = 5.rslope3(i,k,2)rslope3(i,k,2)*rslope(i,k,3) &
1014	! +2.rslope3(i,k,2)rslope2(i,k,2)*rslope2(i,k,3) &
1015	! +.5rslope2(i,k,2)rslope2(i,k,2)*rslope3(i,k,3)
1016	! if(supcol.gt.0) then
1017	! egs = exp(-0.09*supcol)
1018	! else
1019	! egs = 1.
1020	! endif
1021	! pgacs(i,k) = pi*2egsn0sn0sfac(i,k)n0gabs(vt2ave-vt2ave) &
1022	! (dens/den(i,k))acrfac
1023	pgacs(i,k) = 0.
1024	! pgacs(i,k) = min(pgacs(i,k),qrs(i,k,2)/dtcld)
1025	! endif
1026	if(supcol.le.0) then
1027	xlf = xlf0
1028	!-------------------------------------------------------------
1029	! pseml: Enhanced melting of snow by accretion of water [HL A34]
1030	! (T>=T0: S->R)
1031	!-------------------------------------------------------------
1032	if(qrs(i,k,2).gt.0.) &
1033	pseml(i,k) = min(max(cliqsupcol(paacw(i,k)+psacr(i,k)) &
1034	/xlf,-qrs(i,k,2)/dtcld),0.)
1035	!-------------------------------------------------------------
1036	! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1037	! (T>=T0: G->R)
1038	!-------------------------------------------------------------
1039	if(qrs(i,k,3).gt.0.) &
1040	pgeml(i,k) = min(max(cliqsupcol(paacw(i,k)+pgacr(i,k)) &
1041	/xlf,-qrs(i,k,3)/dtcld),0.)
1042	endif
1043	if(supcol.gt.0) then
1044	!-------------------------------------------------------------
1045	! pidep: Deposition/Sublimation rate of ice [HDC 9]
1046	! (T<T0: V->I or I->V)
1047	!-------------------------------------------------------------
1048	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
1049	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1050	supice = satdt-prevp(i,k)
1051	if(pidep(i,k).lt.0.) then
1052	pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1053	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1054	else
1055	pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1056	endif
1057	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1058	endif
1059	!-------------------------------------------------------------
1060	! psdep: deposition/sublimation rate of snow [HDC 14]
1061	! (T<T0: V->S or S->V)
1062	!-------------------------------------------------------------
1063	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1064	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
1065	psdep(i,k) = (rh(i,k,2)-1.)n0sfac(i,k)(precs1 &
1066	rslope2(i,k,2)+precs2work2(i,k) &
1067	*coeres)/work1(i,k,2)
1068	supice = satdt-prevp(i,k)-pidep(i,k)
1069	if(psdep(i,k).lt.0.) then
1070	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1071	psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1072	else
1073	psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1074	endif
1075	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1076	ifsat = 1
1077	endif
1078	!-------------------------------------------------------------
1079	! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1080	! (T<T0: V->G or G->V)
1081	!-------------------------------------------------------------
1082	if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
1083	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
1084	pgdep(i,k) = (rh(i,k,2)-1.)(precg1rslope2(i,k,3) &
1085	+precg2work2(i,k)coeres)/work1(i,k,2)
1086	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1087	if(pgdep(i,k).lt.0.) then
1088	pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
1089	pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1090	else
1091	pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1092	endif
1093	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1094	abs(satdt)) ifsat = 1
1095	endif
1096	!-------------------------------------------------------------
1097	! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1098	! (T<T0: V->I)
1099	!-------------------------------------------------------------
1100	if(supsat.gt.0.and.ifsat.ne.1) then
1101	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1102	xni0 = 1.e3exp(0.1supcol)
1103	roqi0 = 4.92e-11xni0*1.33
1104	pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
1105	/dtcld)
1106	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1107	endif
1108	!
1109	!-------------------------------------------------------------
1110	! psaut: conversion(aggregation) of ice to snow [HDC 12]
1111	! (T<T0: I->S)
1112	!-------------------------------------------------------------
1113	if(qci(i,k,2).gt.0.) then
1114	qimax = roqimax/den(i,k)
1115	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1116	endif
1117	!
1118	!-------------------------------------------------------------
1119	! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1120	! (T<T0: S->G)
1121	!-------------------------------------------------------------
1122	if(qrs(i,k,2).gt.0.) then
1123	alpha2 = 1.e-3exp(0.09(-supcol))
1124	pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)) &
1125	,qrs(i,k,2)/dtcld)
1126	endif
1127	endif
1128	!
1129	!-------------------------------------------------------------
1130	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1131	! (T>=T0: S->V)
1132	!-------------------------------------------------------------
1133	if(supcol.lt.0.) then
1134	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
1135	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
1136	psevp(i,k) = (rh(i,k,1)-1.)n0sfac(i,k)(precs1 &
1137	rslope2(i,k,2)+precs2work2(i,k) &
1138	*coeres)/work1(i,k,1)
1139	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1140	endif
1141	!-------------------------------------------------------------
1142	! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
1143	! (T>=T0: G->V)
1144	!-------------------------------------------------------------
1145	if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
1146	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
1147	pgevp(i,k) = (rh(i,k,1)-1.)(precg1rslope2(i,k,3) &
1148	+precg2work2(i,k)coeres)/work1(i,k,1)
1149	pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
1150	endif
1151	endif
1152	enddo
1153	enddo
1154	!
1155	!
1156	!----------------------------------------------------------------
1157	! check mass conservation of generation terms and feedback to the
1158	! large scale
1159	!
1160	do k = kts, kte
1161	do i = its, ite
1162	!
1163	delta2=0.
1164	delta3=0.
1165	if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
1166	if(qrs(i,k,1).lt.1.e-4) delta3=1.
1167	if(t(i,k).le.t0c) then
1168	!
1169	! cloud water
1170	!
1171	value = max(qmin,qci(i,k,1))
1172	source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1173	if (source.gt.value) then
1174	factor = value/source
1175	praut(i,k) = praut(i,k)*factor
1176	pracw(i,k) = pracw(i,k)*factor
1177	paacw(i,k) = paacw(i,k)*factor
1178	endif
1179	!
1180	! cloud ice
1181	!
1182	value = max(qmin,qci(i,k,2))
1183	source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k) &
1184	+psaci(i,k)+pgaci(i,k))*dtcld
1185	if (source.gt.value) then
1186	factor = value/source
1187	psaut(i,k) = psaut(i,k)*factor
1188	pigen(i,k) = pigen(i,k)*factor
1189	pidep(i,k) = pidep(i,k)*factor
1190	praci(i,k) = praci(i,k)*factor
1191	psaci(i,k) = psaci(i,k)*factor
1192	pgaci(i,k) = pgaci(i,k)*factor
1193	endif
1194	!
1195	! rain
1196	!
1197	value = max(qmin,qrs(i,k,1))
1198	source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) &
1199	+psacr(i,k)+pgacr(i,k))*dtcld
1200	if (source.gt.value) then
1201	factor = value/source
1202	praut(i,k) = praut(i,k)*factor
1203	prevp(i,k) = prevp(i,k)*factor
1204	pracw(i,k) = pracw(i,k)*factor
1205	piacr(i,k) = piacr(i,k)*factor
1206	psacr(i,k) = psacr(i,k)*factor
1207	pgacr(i,k) = pgacr(i,k)*factor
1208	endif
1209	!
1210	! snow
1211	!
1212	value = max(qmin,qrs(i,k,2))
1213	source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) &
1214	+piacr(i,k)delta3+praci(i,k)delta3 &
1215	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2 &
1216	+psaci(i,k)-pgacs(i,k) )*dtcld
1217	if (source.gt.value) then
1218	factor = value/source
1219	psdep(i,k) = psdep(i,k)*factor
1220	psaut(i,k) = psaut(i,k)*factor
1221	pgaut(i,k) = pgaut(i,k)*factor
1222	paacw(i,k) = paacw(i,k)*factor
1223	piacr(i,k) = piacr(i,k)*factor
1224	praci(i,k) = praci(i,k)*factor
1225	psaci(i,k) = psaci(i,k)*factor
1226	pracs(i,k) = pracs(i,k)*factor
1227	psacr(i,k) = psacr(i,k)*factor
1228	pgacs(i,k) = pgacs(i,k)*factor
1229	endif
1230	!
1231	! graupel
1232	!
1233	value = max(qmin,qrs(i,k,3))
1234	source = -(pgdep(i,k)+pgaut(i,k) &
1235	+piacr(i,k)(1.-delta3)+praci(i,k)(1.-delta3) &
1236	+psacr(i,k)(1.-delta2)+pracs(i,k)(1.-delta2) &
1237	+pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1238	if (source.gt.value) then
1239	factor = value/source
1240	pgdep(i,k) = pgdep(i,k)*factor
1241	pgaut(i,k) = pgaut(i,k)*factor
1242	piacr(i,k) = piacr(i,k)*factor
1243	praci(i,k) = praci(i,k)*factor
1244	psacr(i,k) = psacr(i,k)*factor
1245	pracs(i,k) = pracs(i,k)*factor
1246	paacw(i,k) = paacw(i,k)*factor
1247	pgaci(i,k) = pgaci(i,k)*factor
1248	pgacr(i,k) = pgacr(i,k)*factor
1249	pgacs(i,k) = pgacs(i,k)*factor
1250	endif
1251	!
1252	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k) &
1253	+pidep(i,k))
1254	! update
1255	q(i,k) = q(i,k)+work2(i,k)*dtcld
1256	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1257	+paacw(i,k)+paacw(i,k))*dtcld,0.)
1258	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1259	+prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1260	-psacr(i,k))*dtcld,0.)
1261	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
1262	+psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1263	*dtcld,0.)
1264	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1265	-pgaut(i,k)+piacr(i,k)*delta3 &
1266	+praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1267	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2) &
1268	*dtcld,0.)
1269	qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
1270	+piacr(i,k)*(1.-delta3) &
1271	+praci(i,k)(1.-delta3)+psacr(i,k)(1.-delta2)&
1272	+pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1273	+pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1274	xlf = xls-xl(i,k)
1275	xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1276	-xl(i,k)prevp(i,k)-xlf(piacr(i,k)+paacw(i,k) &
1277	+paacw(i,k)+pgacr(i,k)+psacr(i,k))
1278	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1279	else
1280	!
1281	! cloud water
1282	!
1283	value = max(qmin,qci(i,k,1))
1284	source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1285	if (source.gt.value) then
1286	factor = value/source
1287	praut(i,k) = praut(i,k)*factor
1288	pracw(i,k) = pracw(i,k)*factor
1289	paacw(i,k) = paacw(i,k)*factor
1290	endif
1291	!
1292	! rain
1293	!
1294	value = max(qmin,qrs(i,k,1))
1295	source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) &
1296	-pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
1297	if (source.gt.value) then
1298	factor = value/source
1299	praut(i,k) = praut(i,k)*factor
1300	prevp(i,k) = prevp(i,k)*factor
1301	pracw(i,k) = pracw(i,k)*factor
1302	paacw(i,k) = paacw(i,k)*factor
1303	pseml(i,k) = pseml(i,k)*factor
1304	pgeml(i,k) = pgeml(i,k)*factor
1305	endif
1306	!
1307	! snow
1308	!
1309	value = max(qcrmin,qrs(i,k,2))
1310	source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1311	if (source.gt.value) then
1312	factor = value/source
1313	pgacs(i,k) = pgacs(i,k)*factor
1314	psevp(i,k) = psevp(i,k)*factor
1315	pseml(i,k) = pseml(i,k)*factor
1316	endif
1317	!
1318	! graupel
1319	!
1320	value = max(qcrmin,qrs(i,k,3))
1321	source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1322	if (source.gt.value) then
1323	factor = value/source
1324	pgacs(i,k) = pgacs(i,k)*factor
1325	pgevp(i,k) = pgevp(i,k)*factor
1326	pgeml(i,k) = pgeml(i,k)*factor
1327	endif
1328	work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1329	! update
1330	q(i,k) = q(i,k)+work2(i,k)*dtcld
1331	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1332	+paacw(i,k)+paacw(i,k))*dtcld,0.)
1333	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1334	+prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1335	-pgeml(i,k))*dtcld,0.)
1336	qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
1337	+pseml(i,k))*dtcld,0.)
1338	qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
1339	+pgeml(i,k))*dtcld,0.)
1340	xlf = xls-xl(i,k)
1341	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1342	-xlf*(pseml(i,k)+pgeml(i,k))
1343	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1344	endif
1345	enddo
1346	enddo
1347	!
1348	! Inline expansion for fpvs
1349	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1350	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1351	hsub = xls
1352	hvap = xlv0
1353	cvap = cpv
1354	ttp=t0c+0.01
1355	dldt=cvap-cliq
1356	xa=-dldt/rv
1357	xb=xa+hvap/(rv*ttp)
1358	dldti=cvap-cice
1359	xai=-dldti/rv
1360	xbi=xai+hsub/(rv*ttp)
1361	do k = kts, kte
1362	do i = its, ite
1363	tr=ttp/t(i,k)
1364	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
1365	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1366	qs(i,k,1) = max(qs(i,k,1),qmin)
1367	tr=ttp/t(i,k)
1368	if(t(i,k).lt.ttp) then
1369	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
1370	else
1371	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
1372	endif
1373	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1374	qs(i,k,2) = max(qs(i,k,2),qmin)
1375	enddo
1376	enddo
1377	!
1378	!----------------------------------------------------------------
1379	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1380	! if there exists additional water vapor condensated/if
1381	! evaporation of cloud water is not enough to remove subsaturation
1382	!
1383	do k = kts, kte
1384	do i = its, ite
1385	work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1386	work2(i,k) = qci(i,k,1)+work1(i,k,1)
1387	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1388	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1389	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1390	q(i,k) = q(i,k)-pcond(i,k)*dtcld
1391	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1392	t(i,k) = t(i,k)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
1393	enddo
1394	enddo
1395	!
1396	!
1397	!----------------------------------------------------------------
1398	! padding for small values
1399	!
1400	do k = kts, kte
1401	do i = its, ite
1402	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1403	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1404	enddo
1405	enddo
1406	enddo ! big loops
1407	END SUBROUTINE wsm62d
1408	! ...................................................................
1409	REAL FUNCTION rgmma(x)
1410	!-------------------------------------------------------------------
1411	IMPLICIT NONE
1412	!-------------------------------------------------------------------
1413	! rgmma function: use infinite product form
1414	REAL :: euler
1415	PARAMETER (euler=0.577215664901532)
1416	REAL :: x, y
1417	INTEGER :: i
1418	if(x.eq.1.)then
1419	rgmma=0.
1420	else
1421	rgmma=xexp(eulerx)
1422	do i=1,10000
1423	y=float(i)
1424	rgmma=rgmma(1.000+x/y)exp(-x/y)
1425	enddo
1426	rgmma=1./rgmma
1427	endif
1428	END FUNCTION rgmma
1429	!
1430	!--------------------------------------------------------------------------
1431	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1432	!--------------------------------------------------------------------------
1433	IMPLICIT NONE
1434	!--------------------------------------------------------------------------
1435	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1436	xai,xbi,ttp,tr
1437	INTEGER ice
1438	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1439	ttp=t0c+0.01
1440	dldt=cvap-cliq
1441	xa=-dldt/rv
1442	xb=xa+hvap/(rv*ttp)
1443	dldti=cvap-cice
1444	xai=-dldti/rv
1445	xbi=xai+hsub/(rv*ttp)
1446	tr=ttp/t
1447	if(t.lt.ttp.and.ice.eq.1) then
1448	fpvs=psat(trxai)exp(xbi*(1.-tr))
1449	else
1450	fpvs=psat(trxa)exp(xb*(1.-tr))
1451	endif
1452	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1453	END FUNCTION fpvs
1454	!-------------------------------------------------------------------
1455	SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,allowed_to_read)
1456	!-------------------------------------------------------------------
1457	IMPLICIT NONE
1458	!-------------------------------------------------------------------
1459	!.... constants which may not be tunable
1460	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
1461	LOGICAL, INTENT(IN) :: allowed_to_read
1462	REAL :: pi
1463	!
1464	pi = 4.*atan(1.)
1465	xlv1 = cl-cpv
1466	!
1467	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
1468	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
1469	!
1470	bvtr1 = 1.+bvtr
1471	bvtr2 = 2.5+.5*bvtr
1472	bvtr3 = 3.+bvtr
1473	bvtr4 = 4.+bvtr
1474	bvtr6 = 6.+bvtr
1475	g1pbr = rgmma(bvtr1)
1476	g3pbr = rgmma(bvtr3)
1477	g4pbr = rgmma(bvtr4) ! 17.837825
1478	g6pbr = rgmma(bvtr6)
1479	g5pbro2 = rgmma(bvtr2) ! 1.8273
1480	pvtr = avtr*g4pbr/6.
1481	eacrr = 1.0
1482	pacrr = pin0ravtrg3pbr.25*eacrr
1483	precr1 = 2.pin0r*.78
1484	precr2 = 2.pin0r.31avtr*.5g5pbro2
1485	xm0 = (di0/dicon)**2
1486	xmmax = (dimax/dicon)**2
1487	roqimax = 2.08e22dimax*8
1488	!
1489	bvts1 = 1.+bvts
1490	bvts2 = 2.5+.5*bvts
1491	bvts3 = 3.+bvts
1492	bvts4 = 4.+bvts
1493	g1pbs = rgmma(bvts1) !.8875
1494	g3pbs = rgmma(bvts3)
1495	g4pbs = rgmma(bvts4) ! 12.0786
1496	g5pbso2 = rgmma(bvts2)
1497	pvts = avts*g4pbs/6.
1498	pacrs = pin0savtsg3pbs.25
1499	precs1 = 4.n0s.65
1500	precs2 = 4.n0s.44avts.5g5pbso2
1501	pidn0r = pidenrn0r
1502	pidn0s = pidensn0s
1503	!
1504	pacrc = pin0savtsg3pbs.25*eacrc
1505	!
1506	bvtg1 = 1.+bvtg
1507	bvtg2 = 2.5+.5*bvtg
1508	bvtg3 = 3.+bvtg
1509	bvtg4 = 4.+bvtg
1510	g1pbg = rgmma(bvtg1)
1511	g3pbg = rgmma(bvtg3)
1512	g4pbg = rgmma(bvtg4)
1513	pacrg = pin0gavtgg3pbg.25
1514	g5pbgo2 = rgmma(bvtg2)
1515	pvtg = avtg*g4pbg/6.
1516	precg1 = 2.pin0g*.78
1517	precg2 = 2.pin0g.31avtg*.5g5pbgo2
1518	pidn0g = pidengn0g
1519	!
1520	rslopermax = 1./lamdarmax
1521	rslopesmax = 1./lamdasmax
1522	rslopegmax = 1./lamdagmax
1523	rsloperbmax = rslopermax ** bvtr
1524	rslopesbmax = rslopesmax ** bvts
1525	rslopegbmax = rslopegmax ** bvtg
1526	rsloper2max = rslopermax * rslopermax
1527	rslopes2max = rslopesmax * rslopesmax
1528	rslopeg2max = rslopegmax * rslopegmax
1529	rsloper3max = rsloper2max * rslopermax
1530	rslopes3max = rslopes2max * rslopesmax
1531	rslopeg3max = rslopeg2max * rslopegmax
1532	!
1533	END SUBROUTINE wsm6init
1534	END MODULE module_mp_wsm6

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

Download in other formats:

Original Format