Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm3_accel.F @ 354

Last change on this file since 354 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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

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

Download in other formats:

Original Format