Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm6.F @ 3567

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

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

Download in other formats:

Original Format