Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm5_accel.F @ 1

Last change on this file since 1 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: 130.2 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	!Including inline expansion statistical function
10	MODULE module_mp_wsm5
11	!
12	!
13	REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
14	REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
15	REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
16	REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
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 ! a constant for terminal velocity of snow
22	REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
23	REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
24	REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
25	REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
26	REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
27	REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
28	REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
29	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
30	REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
31	REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
32	REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
33	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
34	REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
35	REAL, SAVE :: &
36	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
37	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
38	precr1,precr2,xmmax,roqimax,bvts1, &
39	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
40	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
41	pidn0s,xlv1,pacrc, &
42	rslopermax,rslopesmax,rslopegmax, &
43	rsloperbmax,rslopesbmax,rslopegbmax, &
44	rsloper2max,rslopes2max,rslopeg2max, &
45	rsloper3max,rslopes3max,rslopeg3max
46	!
47	! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
48	!
49	CONTAINS
50	!===================================================================
51	!
52	SUBROUTINE wsm5(th, q, qc, qr, qi, qs &
53	,den, pii, p, delz &
54	,delt,g, cpd, cpv, rd, rv, t0c &
55	,ep1, ep2, qmin &
56	,XLS, XLV0, XLF0, den0, denr &
57	,cliq,cice,psat &
58	,rain, rainncv &
59	,snow, snowncv &
60	,sr &
61	,ids,ide, jds,jde, kds,kde &
62	,ims,ime, jms,jme, kms,kme &
63	,its,ite, jts,jte, kts,kte &
64	)
65	#ifdef _OPENMP
66	use omp_lib
67	#endif
68	!-------------------------------------------------------------------
69	IMPLICIT NONE
70	!-------------------------------------------------------------------
71	!
72	! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
73	! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
74	! number concentration is a function of temperature, and seperate assumption
75	! is developed, in which ice crystal number concentration is a function
76	! of ice amount. A theoretical background of the ice-microphysics and related
77	! processes in the WSMMPs are described in Hong et al. (2004).
78	! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
79	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
80	!
81	! WSM5 cloud scheme
82	!
83	! Coded by Song-You Hong (Yonsei Univ.)
84	! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
85	! Summer 2002
86	!
87	! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
88	! Summer 2003
89	!
90	! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
91	! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
92	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
93	!
94	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
95	ims,ime, jms,jme, kms,kme , &
96	its,ite, jts,jte, kts,kte
97	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
98	INTENT(INOUT) :: &
99	th, &
100	q, &
101	qc, &
102	qi, &
103	qr, &
104	qs
105	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
106	INTENT(IN ) :: &
107	den, &
108	pii, &
109	p, &
110	delz
111	REAL, INTENT(IN ) :: delt, &
112	g, &
113	rd, &
114	rv, &
115	t0c, &
116	den0, &
117	cpd, &
118	cpv, &
119	ep1, &
120	ep2, &
121	qmin, &
122	XLS, &
123	XLV0, &
124	XLF0, &
125	cliq, &
126	cice, &
127	psat, &
128	denr
129	REAL, DIMENSION( ims:ime , jms:jme ), &
130	INTENT(INOUT) :: rain, &
131	rainncv, &
132	sr
133
134	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
135	INTENT(INOUT) :: snow, &
136	snowncv
137
138	! LOCAL VAR
139	REAL, DIMENSION( its:ite , kts:kte ) :: t
140	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
141	CHARACTER*256 :: emess
142	INTEGER :: mkx_test
143	INTEGER :: i,j,k
144
145	#ifdef _ACCEL_PROF
146	INTEGER :: l
147	real*8 wsm3_t(8,256), wsm5_t(8,256), t1, t2
148	common /wsm_times/ wsm3_t(8,256), wsm5_t(8,256)
149	#endif
150
151
152	!-------------------------------------------------------------------
153
154	#ifdef _ACCEL_PROF
155	call cpu_time(t1)
156	#endif
157
158	#ifndef RUN_ON_GPU
159
160	#ifdef _ACCEL
161
162	! Need to send th, pii, qc, qi, qr, qs
163	! Don't send t
164
165	CALL wsm52D(th, pii, q, qc, qr, qi, qs &
166	,den &
167	,p, delz &
168	,delt,g, cpd, cpv, rd, rv, t0c &
169	,ep1, ep2, qmin &
170	,XLS, XLV0, XLF0, den0, denr &
171	,cliq,cice,psat &
172	,rain,rainncv &
173	,sr &
174	,ids,ide, jds,jde, kds,kde &
175	,ims,ime, jms,jme, kms,kme &
176	,its,ite, jts,jte, kts,kte &
177	,snow,snowncv &
178	)
179
180	#else
181
182	DO j=jts,jte
183	DO k=kts,kte
184	DO i=its,ite
185	t(i,k)=th(i,k,j)*pii(i,k,j)
186	qci(i,k,1) = qc(i,k,j)
187	qci(i,k,2) = qi(i,k,j)
188	qrs(i,k,1) = qr(i,k,j)
189	qrs(i,k,2) = qs(i,k,j)
190	ENDDO
191	ENDDO
192
193	! Sending array starting locations of optional variables may cause
194	! troubles, so we explicitly change the call.
195
196	CALL wsm52D(t, q(ims,kms,j), qci, qrs &
197	,den(ims,kms,j) &
198	,p(ims,kms,j), delz(ims,kms,j) &
199	,delt,g, cpd, cpv, rd, rv, t0c &
200	,ep1, ep2, qmin &
201	,XLS, XLV0, XLF0, den0, denr &
202	,cliq,cice,psat &
203	,j &
204	,rain(ims,j),rainncv(ims,j) &
205	,sr(ims,j) &
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	)
211
212	DO K=kts,kte
213	DO I=its,ite
214	th(i,k,j)=t(i,k)/pii(i,k,j)
215	qc(i,k,j) = qci(i,k,1)
216	qi(i,k,j) = qci(i,k,2)
217	qr(i,k,j) = qrs(i,k,1)
218	qs(i,k,j) = qrs(i,k,2)
219	ENDDO
220	ENDDO
221	ENDDO
222	#endif
223	#else
224	CALL get_wsm5_gpu_levels ( mkx_test )
225	IF ( mkx_test .LT. kte ) THEN
226	WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ', &
227	mkx_test,' < ',kte
228	CALL wrf_error_fatal(emess)
229	ENDIF
230	CALL wsm5_host ( &
231	th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte) &
232	,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte) &
233	,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte) &
234	,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte) &
235	,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte) &
236	,delt &
237	,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte) &
238	,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte) &
239	,sr(its:ite,jts:jte) &
240	,its, ite, jts, jte, kts, kte &
241	,its, ite, jts, jte, kts, kte &
242	,its, ite, jts, jte, kts, kte &
243	)
244	#endif
245
246
247	#ifdef _ACCEL_PROF
248	call cpu_time(t2)
249	#ifdef _OPENMP
250	l = omp_get_thread_num() + 1
251	#else
252	l = 1
253	#endif
254	wsm5_t(1,l) = wsm5_t(1,l) + (t2 - t1)
255	#endif
256
257	END SUBROUTINE wsm5
258
259
260	#ifdef _ACCEL
261
262	!===================================================================
263	!
264	SUBROUTINE wsm52D(th, pii, q, qc, qr, qi, qqs, den, p, delz &
265	,delt,g, cpd, cpv, rd, rv, t0c &
266	,ep1, ep2, qmin &
267	,XLS, XLV0, XLF0, den0, denr &
268	,cliq,cice,psat &
269	,rain,rainncv &
270	,sr &
271	,ids,ide, jds,jde, kds,kde &
272	,ims,ime, jms,jme, kms,kme &
273	,its,ite, jts,jte, kts,kte &
274	,snow,snowncv &
275	)
276	!-------------------------------------------------------------------
277	IMPLICIT NONE
278	!-------------------------------------------------------------------
279	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
280	ims,ime, jms,jme, kms,kme , &
281	its,ite, jts,jte, kts,kte
282	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
283	INTENT(INOUT) :: &
284	th
285	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
286	INTENT(IN) :: &
287	pii
288	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
289	INTENT(INOUT) :: &
290	qc, &
291	qr, &
292	qi, &
293	qqs
294	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
295	INTENT(INOUT) :: &
296	q
297	REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
298	INTENT(IN ) :: &
299	den, &
300	p, &
301	delz
302	REAL, INTENT(IN ) :: delt, &
303	g, &
304	cpd, &
305	cpv, &
306	t0c, &
307	den0, &
308	rd, &
309	rv, &
310	ep1, &
311	ep2, &
312	qmin, &
313	XLS, &
314	XLV0, &
315	XLF0, &
316	cliq, &
317	cice, &
318	psat, &
319	denr
320	REAL, DIMENSION( ims:ime, jms:jme ), &
321	INTENT(INOUT) :: rain, &
322	rainncv, &
323	sr
324
325	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
326	INTENT(INOUT) :: snow, &
327	snowncv
328
329	! LOCAL VAR
330	REAL, DIMENSION( its:ite , kts:kte , 2) :: &
331	rh, &
332	qs, &
333	rslope, &
334	rslope2, &
335	rslope3, &
336	rslopeb, &
337	falk, &
338	fall, &
339	work1
340	REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: &
341	t
342	REAL, DIMENSION( its:ite , kts:kte , 2 ) :: &
343	qci, &
344	qrs
345	REAL, DIMENSION( its:ite , kts:kte ) :: &
346	falkc, &
347	fallc, &
348	xl, &
349	cpm, &
350	denfac, &
351	xni, &
352	n0sfac, &
353	work2, &
354	work1c, &
355	work2c
356
357	REAL, DIMENSION( its:ite , kts:kte ) :: &
358	pigen, &
359	pidep, &
360	psdep, &
361	praut, &
362	psaut, &
363	prevp, &
364	psevp, &
365	pracw, &
366	psacw, &
367	psaci, &
368	pcond, &
369	psmlt
370	INTEGER :: &
371	mstep, &
372	numdt
373	REAL :: rmstep
374	REAL dtcldden, rdelz, rdtcld
375	LOGICAL :: flgcld
376
377	#define WSM_NO_CONDITIONAL_IN_VECTOR
378	#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
379	REAL :: xal, xbl
380	#endif
381
382	REAL :: pi, &
383	cpmcal, xlcal, lamdar, lamdas, diffus, &
384	viscos, xka, venfac, conden, diffac, &
385	x, y, z, a, b, c, d, e, &
386	qdt, holdrr, holdrs, supcol, supcolt, pvt, &
387	coeres, supsat, dtcld, xmi, eacrs, satdt, &
388	vt2i,vt2s,acrfac, &
389	qimax, diameter, xni0, roqi0, &
390	fallsum, fallsum_qsi, xlwork2, factor, source, &
391	value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
392	! variables for optimization
393	REAL, DIMENSION( its:ite ) :: tvec1
394	REAL :: temp
395	INTEGER :: i, j, k, &
396	iprt, latd, lond, loop, loops, ifsat, n
397	! Temporaries used for inlining fpvs function
398	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
399	REAL :: logtr
400	!
401	!=================================================================
402	! compute internal functions
403	!
404	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
405	xlcal(x) = xlv0-xlv1*(x-t0c)
406	!----------------------------------------------------------------
407	! size distributions: (x=mixing ratio, y=air density):
408	! valid for mixing ratio > 1.e-9 kg/kg.
409	!
410	! Optimizatin : A*B => exp(log(A)(B))
411	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
412	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
413	!
414	!----------------------------------------------------------------
415	! diffus: diffusion coefficient of the water vapor
416	! viscos: kinematic viscosity(m2s-1)
417	! diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
418	! viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
419	! xka(x,y) = 1.414e3viscos(x,y)y
420	! diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
421	! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
422	! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
423	! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
424	!
425	!
426	pi = 4. * atan(1.)
427	!
428	!----------------------------------------------------------------
429	! paddint 0 for negative values generated by dynamics
430	!
431
432	!
433	! Moved outside of accelerator region
434	!
435	loops = max(nint(delt/dtcldcr),1)
436	dtcld = delt/loops
437	if(delt.le.dtcldcr) dtcld = delt
438	!
439
440	!....!$acc local(t) &
441
442	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
443
444	!$acc region &
445	!$acc local(t) &
446	!$acc copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
447	!$acc copyout(snowncv(:,:),rainncv(:,:),sr(:,:)) &
448	!$acc copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
449	!$acc copy(th(:,:,:),q(:,:,:),snow(:,:),rain(:,:))
450	!$acc do &
451	!$acc private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
452	!$acc private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
453	!$acc private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
454	!$acc private(praut,psaut,prevp,psevp) &
455	!$acc private(pracw,psacw,psaci,pcond,psmlt) &
456	!$acc parallel
457	do j = jts, jte
458	!$acc do &
459	!$acc private(numdt,mstep) &
460	!$acc kernel vector
461	do i = its, ite
462	do k = kts, kte
463	t(i,k,j)=th(i,k,j)*pii(i,k,j)
464	qci(i,k,1) = max(qc(i,k,j),0.0)
465	qci(i,k,2) = max(qi(i,k,j),0.0)
466	qrs(i,k,1) = max(qr(i,k,j),0.0)
467	qrs(i,k,2) = max(qqs(i,k,j),0.0)
468	enddo
469	!
470	!----------------------------------------------------------------
471	! latent heat for phase changes and heat capacity. neglect the
472	! changes during microphysical process calculation
473	! emanuel(1994)
474	!
475	do k = kts, kte
476	cpm(i,k) = cpmcal(q(i,k,j))
477	xl(i,k) = xlcal(t(i,k,j))
478	enddo
479	!
480	!----------------------------------------------------------------
481	! compute the minor time steps.
482	!
483	! loops = max(nint(delt/dtcldcr),1)
484	! dtcld = delt/loops
485	! if(delt.le.dtcldcr) dtcld = delt
486	!
487	do loop = 1,loops
488	!
489	!----------------------------------------------------------------
490	! initialize the large scale variables
491	!
492	mstep = 1
493	flgcld = .true.
494	!
495	do k = kts, kte
496	denfac(i,k) = sqrt(den0/den(i,k,j))
497	enddo
498	! do k = kts, kte
499	! CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
500	! do i = its, ite
501	! tvec1(i) = tvec1(i)*den0
502	! enddo
503	! CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
504	! enddo
505	!
506	! Inline expansion for fpvs
507	! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
508	! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
509	hsub = xls
510	hvap = xlv0
511	cvap = cpv
512	ttp=t0c+0.01
513	dldt=cvap-cliq
514	xa=-dldt/rv
515	xb=xa+hvap/(rv*ttp)
516	dldti=cvap-cice
517	xai=-dldti/rv
518	xbi=xai+hsub/(rv*ttp)
519
520	! this is for compilers where the conditional inhibits vectorization
521	#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
522	do k = kts, kte
523	if(t(i,k,j).lt.ttp) then
524	xal = xai
525	xbl = xbi
526	else
527	xal = xa
528	xbl = xb
529	endif
530	tr=ttp/t(i,k,j)
531	logtr=log(tr)
532	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
533	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
534	qs(i,k,1) = max(qs(i,k,1),qmin)
535	rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
536	qs(i,k,2)=psatexp(logtr(xal)+xbl*(1.-tr))
537	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
538	qs(i,k,2) = max(qs(i,k,2),qmin)
539	rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
540	enddo
541	#else
542	do k = kts, kte
543	tr=ttp/t(i,k,j)
544	logtr=log(tr)
545	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
546	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
547	qs(i,k,1) = max(qs(i,k,1),qmin)
548	rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
549	if(t(i,k,j).lt.ttp) then
550	qs(i,k,2)=psatexp(logtr(xai)+xbi*(1.-tr))
551	else
552	qs(i,k,2)=psatexp(logtr(xa)+xb*(1.-tr))
553	endif
554	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
555	qs(i,k,2) = max(qs(i,k,2),qmin)
556	rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
557	enddo
558	#endif
559
560	!
561	!----------------------------------------------------------------
562	! initialize the variables for microphysical physics
563	!
564	!
565	do k = kts, kte
566	prevp(i,k) = 0.
567	psdep(i,k) = 0.
568	praut(i,k) = 0.
569	psaut(i,k) = 0.
570	pracw(i,k) = 0.
571	psaci(i,k) = 0.
572	psacw(i,k) = 0.
573	pigen(i,k) = 0.
574	pidep(i,k) = 0.
575	pcond(i,k) = 0.
576	psmlt(i,k) = 0.
577	psevp(i,k) = 0.
578	falk(i,k,1) = 0.
579	falk(i,k,2) = 0.
580	fall(i,k,1) = 0.
581	fall(i,k,2) = 0.
582	fallc(i,k) = 0.
583	falkc(i,k) = 0.
584	xni(i,k) = 1.e3
585	enddo
586	!
587	!----------------------------------------------------------------
588	! compute the fallout term:
589	! first, vertical terminal velosity for minor loops
590	!
591	do k = kts, kte
592	supcol = t0c-t(i,k,j)
593	!---------------------------------------------------------------
594	! n0s: Intercept parameter for snow [m-4] [HDC 6]
595	!---------------------------------------------------------------
596	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
597	if(qrs(i,k,1).le.qcrmin)then
598	rslope(i,k,1) = rslopermax
599	rslopeb(i,k,1) = rsloperbmax
600	rslope2(i,k,1) = rsloper2max
601	rslope3(i,k,1) = rsloper3max
602	else
603	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
604	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
605	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
606	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
607	endif
608	if(qrs(i,k,2).le.qcrmin)then
609	rslope(i,k,2) = rslopesmax
610	rslopeb(i,k,2) = rslopesbmax
611	rslope2(i,k,2) = rslopes2max
612	rslope3(i,k,2) = rslopes3max
613	else
614	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
615	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
616	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
617	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
618	endif
619	!-------------------------------------------------------------
620	! Ni: ice crystal number concentraiton [HDC 5c]
621	!-------------------------------------------------------------
622	! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
623	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
624	temp = (den(i,k,j)*max(qci(i,k,2),qmin))
625	temp = sqrt(sqrt(temptemptemp))
626	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
627	enddo
628	!
629	numdt = 1
630	do k = kte, kts, -1
631	work1(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)/delz(i,k,j)
632	work1(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)/delz(i,k,j)
633	numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
634	if(numdt.ge.mstep) mstep = numdt
635	enddo
636	rmstep = 1./mstep
637	!
638	do n = 1, mstep
639	k = kte
640	! falk(i,k,1) = den(i,k,j)qrs(i,k,1)work1(i,k,1)/mstep(i)
641	! falk(i,k,2) = den(i,k,j)qrs(i,k,2)work1(i,k,2)/mstep(i)
642	falk(i,k,1) = den(i,k,j)qrs(i,k,1)work1(i,k,1)*rmstep
643	falk(i,k,2) = den(i,k,j)qrs(i,k,2)work1(i,k,2)*rmstep
644	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
645	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
646	! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
647	! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
648	dtcldden = dtcld/den(i,k,j)
649	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
650	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
651	! endif
652	do k = kte-1, kts, -1
653	falk(i,k,1) = den(i,k,j)qrs(i,k,1)work1(i,k,1)*rmstep
654	falk(i,k,2) = den(i,k,j)qrs(i,k,2)work1(i,k,2)*rmstep
655	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
656	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
657	dtcldden = dtcld/den(i,k,j)
658	rdelz = 1./delz(i,k,j)
659	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
660	delz(i,k+1,j)rdelz)*dtcldden,0.)
661	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
662	delz(i,k+1,j)rdelz)*dtcldden,0.)
663	enddo
664	do k = kte, kts, -1
665	if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
666	!----------------------------------------------------------------
667	! psmlt: melting of snow [HL A33] [RH83 A25]
668	! (T>T0: S->R)
669	!----------------------------------------------------------------
670	xlf = xlf0
671	! work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
672	work2(i,k)= (exp(log(((1.496e-6((t(i,k,j))sqrt(t(i,k,j))) &
673	/((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5 &
674	exp(log(t(i,k,j))(1.81))/p(i,k,j)))) &
675	((.3333333)))/sqrt((1.496e-6((t(i,k,j)) &
676	*sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j)))) &
677	*sqrt(sqrt(den0/(den(i,k,j)))))
678	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
679	! psmlt(i,k) = xka(t(i,k,j),den(i,k,j))/xlf(t0c-t(i,k,j))pi/2. &
680	! n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
681	! work2(i,k)coeres)
682	psmlt(i,k) = (1.414e3(1.496e-6((t(i,k,j))*sqrt(t(i,k,j))) &
683	/((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j))) &
684	/xlf(t0c-t(i,k,j))pi/2. &
685	n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
686	work2(i,k)coeres)
687	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep, &
688	-qrs(i,k,2)/mstep),0.)
689	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
690	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
691	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
692	endif
693	enddo
694	enddo
695
696
697	!---------------------------------------------------------------
698	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
699	!---------------------------------------------------------------
700	mstep = 1
701	numdt = 1
702	do k = kte, kts, -1
703	if(qci(i,k,2).le.0.) then
704	work2c(i,k) = 0.
705	else
706	xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
707	! diameter = min(dicon * sqrt(xmi),dimax)
708	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
709	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
710	work2c(i,k) = work1c(i,k)/delz(i,k,j)
711	endif
712	numdt = max(nint(work2c(i,k)*dtcld+.5),1)
713	if(numdt.ge.mstep) mstep = numdt
714	enddo
715	!
716	do n = 1, mstep
717	k = kte
718	falkc(i,k) = den(i,k,j)qci(i,k,2)work2c(i,k)/mstep
719	holdc = falkc(i,k)
720	fallc(i,k) = fallc(i,k)+falkc(i,k)
721	holdci = qci(i,k,2)
722	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
723	! endif
724	do k = kte-1, kts, -1
725	falkc(i,k) = den(i,k,j)qci(i,k,2)work2c(i,k)/mstep
726	holdc = falkc(i,k)
727	fallc(i,k) = fallc(i,k)+falkc(i,k)
728	holdci = qci(i,k,2)
729	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
730	delz(i,k+1,j)/delz(i,k,j))dtcld/den(i,k,j),0.)
731	! endif
732	enddo
733	enddo
734	!
735	!
736	!----------------------------------------------------------------
737	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
738	!
739	fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
740	fallsum_qsi = fall(i,1,2)+fallc(i,1)
741	rainncv(i,j) = 0.
742	if(fallsum.gt.0.) then
743	rainncv(i,j) = fallsumdelz(i,1,j)/denrdtcld*1000.
744	rain(i,j) = fallsumdelz(i,1,j)/denrdtcld*1000. + rain(i,j)
745	endif
746	snowncv(i,j) = 0.
747	if(fallsum_qsi.gt.0.) then
748	snowncv(i,j) = fallsum_qsidelz(i,kts,j)/denrdtcld*1000.
749	snow(i,j) = fallsum_qsidelz(i,kts,j)/denrdtcld*1000. + snow(i,j)
750	endif
751	sr(i,j) = 0.
752	if(fallsum.gt.0.)sr(i,j)=fallsum_qsidelz(i,kts,j)/denrdtcld*1000. &
753	/(rainncv(i,j)+1.e-12)
754	!
755	!---------------------------------------------------------------
756	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
757	! (T>T0: I->C)
758	!---------------------------------------------------------------
759	do k = kts, kte
760	supcol = t0c-t(i,k,j)
761	xlf = xls-xl(i,k)
762	if(supcol.lt.0.) xlf = xlf0
763	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
764	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
765	t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
766	qci(i,k,2) = 0.
767	endif
768	!---------------------------------------------------------------
769	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
770	! (T<-40C: C->I)
771	!---------------------------------------------------------------
772	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
773	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
774	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
775	qci(i,k,1) = 0.
776	endif
777	!---------------------------------------------------------------
778	! pihtf: heterogeneous freezing of cloud water [HL A44]
779	! (T0>T>-40C: C->I)
780	!---------------------------------------------------------------
781	if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
782	supcolt=min(supcol,50.)
783	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
784	! den(i,k,j)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
785	pfrzdtc = min(pfrz1(exp(pfrz2supcolt)-1.) &
786	den(i,k,j)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
787	qci(i,k,2) = qci(i,k,2) + pfrzdtc
788	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
789	qci(i,k,1) = qci(i,k,1)-pfrzdtc
790	endif
791	!---------------------------------------------------------------
792	! psfrz: freezing of rain water [HL A20] [LFO 45]
793	! (T<T0, R->S)
794	!---------------------------------------------------------------
795	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
796	supcolt=min(supcol,50.)
797	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k,j) &
798	! (exp(pfrz2supcol)-1.)rslope(i,k,1)7dtcld, &
799	! qrs(i,k,1))
800	temp = rslope(i,k,1)
801	temp = temptemptemptemptemptemptemp
802	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k,j) &
803	(exp(pfrz2supcolt)-1.)tempdtcld, &
804	qrs(i,k,1))
805	qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
806	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
807	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
808	endif
809	enddo
810	!
811	!----------------------------------------------------------------
812	! rsloper: reverse of the slope parameter of the rain(m)
813	! xka: thermal conductivity of air(jm-1s-1k-1)
814	! work1: the thermodynamic term in the denominator associated with
815	! heat conduction and vapor diffusion
816	! (ry88, y93, h85)
817	! work2: parameter associated with the ventilation effects(y93)
818	!
819	do k = kts, kte
820	if(qrs(i,k,1).le.qcrmin)then
821	rslope(i,k,1) = rslopermax
822	rslopeb(i,k,1) = rsloperbmax
823	rslope2(i,k,1) = rsloper2max
824	rslope3(i,k,1) = rsloper3max
825	else
826	! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
827	rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
828	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
829	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
830	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
831	endif
832	if(qrs(i,k,2).le.qcrmin)then
833	rslope(i,k,2) = rslopesmax
834	rslopeb(i,k,2) = rslopesbmax
835	rslope2(i,k,2) = rslopes2max
836	rslope3(i,k,2) = rslopes3max
837	else
838	! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
839	rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2)) &
840	*(den(i,k,j))))))
841	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
842	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
843	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
844	endif
845	enddo
846	!
847	do k = kts, kte
848	! work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
849	work1(i,k,1) = ((((den(i,k,j))(xl(i,k))(xl(i,k)))*((t(i,k,j))+120.) &
850	(den(i,k,j)))/(1.414e3(1.496e-6((t(i,k,j))sqrt(t(i,k,j))))&
851	(den(i,k,j))(rv(t(i,k,j))(t(i,k,j))))) &
852	+ p(i,k,j)/((qs(i,k,1))(8.794e-5exp(log(t(i,k,j))*(1.81))))
853	! work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
854	work1(i,k,2) = ((((den(i,k,j))(xls)(xls))((t(i,k,j))+120.)(den(i,k,j)))&
855	/(1.414e3(1.496e-6((t(i,k,j))sqrt(t(i,k,j))))(den(i,k,j)) &
856	(rv(t(i,k,j))*(t(i,k,j)))) &
857	+ p(i,k,j)/(qs(i,k,2)(8.794e-5exp(log(t(i,k,j))*(1.81)))))
858	! work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
859	work2(i,k) = (exp(.3333333log(((1.496e-6 ((t(i,k,j))*sqrt(t(i,k,j)))) &
860	p(i,k,j))/(((t(i,k,j))+120.)den(i,k,j)*(8.794e-5 &
861	exp(log(t(i,k,j))(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
862	/sqrt((1.496e-6((t(i,k,j))sqrt(t(i,k,j)))) &
863	/(((t(i,k,j))+120.)*den(i,k,j)))
864	enddo
865	!
866	!===============================================================
867	!
868	! warm rain processes
869	!
870	! - follows the processes in RH83 and LFO except for autoconcersion
871	!
872	!===============================================================
873	!
874	do k = kts, kte
875	supsat = max(q(i,k,j),qmin)-qs(i,k,1)
876	satdt = supsat/dtcld
877	!---------------------------------------------------------------
878	! praut: auto conversion rate from cloud to rain [HDC 16]
879	! (C->R)
880	!---------------------------------------------------------------
881	if(qci(i,k,1).gt.qc0) then
882	praut(i,k) = qck1exp(log(qci(i,k,1))((7./3.)))
883	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
884	endif
885	!---------------------------------------------------------------
886	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
887	! (C->R)
888	!---------------------------------------------------------------
889	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
890	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
891	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
892	endif
893	!---------------------------------------------------------------
894	! prevp: evaporation/condensation rate of rain [HDC 14]
895	! (V->R or R->V)
896	!---------------------------------------------------------------
897	if(qrs(i,k,1).gt.0.) then
898	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
899	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
900	+precr2work2(i,k)coeres)/work1(i,k,1)
901	if(prevp(i,k).lt.0.) then
902	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
903	prevp(i,k) = max(prevp(i,k),satdt/2)
904	else
905	prevp(i,k) = min(prevp(i,k),satdt/2)
906	endif
907	endif
908	enddo
909	!
910	!===============================================================
911	!
912	! cold rain processes
913	!
914	! - follows the revised ice microphysics processes in HDC
915	! - the processes same as in RH83 and RH84 and LFO behave
916	! following ice crystal hapits defined in HDC, inclduing
917	! intercept parameter for snow (n0s), ice crystal number
918	! concentration (ni), ice nuclei number concentration
919	! (n0i), ice diameter (d)
920	!
921	!===============================================================
922	!
923	rdtcld = 1./dtcld
924	do k = kts, kte
925	supcol = t0c-t(i,k,j)
926	supsat = max(q(i,k,j),qmin)-qs(i,k,2)
927	satdt = supsat/dtcld
928	ifsat = 0
929	!-------------------------------------------------------------
930	! Ni: ice crystal number concentraiton [HDC 5c]
931	!-------------------------------------------------------------
932	! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
933	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
934	temp = (den(i,k,j)*max(qci(i,k,2),qmin))
935	temp = sqrt(sqrt(temptemptemp))
936	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
937	eacrs = exp(0.07*(-supcol))
938	!
939	if(supcol.gt.0) then
940	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
941	xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
942	diameter = min(dicon * sqrt(xmi),dimax)
943	vt2i = 1.49e4diameter*1.31
944	vt2s = pvtsrslopeb(i,k,2)denfac(i,k)
945	!-------------------------------------------------------------
946	! psaci: Accretion of cloud ice by rain [HDC 10]
947	! (T<T0: I->S)
948	!-------------------------------------------------------------
949	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
950	+diameter*2rslope(i,k,2)
951	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
952	abs(vt2s-vt2i)acrfac/4.
953	endif
954	endif
955	!-------------------------------------------------------------
956	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
957	! (T<T0: C->S, and T>=T0: C->R)
958	!-------------------------------------------------------------
959	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
960	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2) &
961	rslopeb(i,k,2)qci(i,k,1)*denfac(i,k) &
962	! ,qci(i,k,1)/dtcld)
963	,qci(i,k,1)*rdtcld)
964	endif
965	if(supcol .gt. 0) then
966	!-------------------------------------------------------------
967	! pidep: Deposition/Sublimation rate of ice [HDC 9]
968	! (T<T0: V->I or I->V)
969	!-------------------------------------------------------------
970	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
971	xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
972	diameter = dicon * sqrt(xmi)
973	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
974	supice = satdt-prevp(i,k)
975	if(pidep(i,k).lt.0.) then
976	! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
977	! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
978	pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
979	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
980	else
981	! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
982	pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
983	endif
984	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
985	endif
986	!-------------------------------------------------------------
987	! psdep: deposition/sublimation rate of snow [HDC 14]
988	! (V->S or S->V)
989	!-------------------------------------------------------------
990	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
991	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
992	psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
993	(precs1rslope2(i,k,2)+precs2 &
994	work2(i,k)coeres)/work1(i,k,2)
995	supice = satdt-prevp(i,k)-pidep(i,k)
996	if(psdep(i,k).lt.0.) then
997	! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
998	! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
999	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1000	psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1001	else
1002	! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1003	psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1004	endif
1005	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1006	ifsat = 1
1007	endif
1008	!-------------------------------------------------------------
1009	! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1010	! (T<T0: V->I)
1011	!-------------------------------------------------------------
1012	if(supsat.gt.0.and.ifsat.ne.1) then
1013	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1014	xni0 = 1.e3exp(0.1supcol)
1015	roqi0 = 4.92e-11exp(log(xni0)(1.33))
1016	pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.)) &
1017	! /dtcld)
1018	*rdtcld)
1019	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1020	endif
1021	!
1022	!-------------------------------------------------------------
1023	! psaut: conversion(aggregation) of ice to snow [HDC 12]
1024	! (T<T0: I->S)
1025	!-------------------------------------------------------------
1026	if(qci(i,k,2).gt.0.) then
1027	qimax = roqimax/den(i,k,j)
1028	! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1029	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1030	endif
1031	endif
1032	!-------------------------------------------------------------
1033	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1034	! (T>T0: S->V)
1035	!-------------------------------------------------------------
1036	if(supcol.lt.0.) then
1037	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
1038	psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1039	! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1040	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1041	endif
1042	enddo
1043	!
1044	!
1045	!----------------------------------------------------------------
1046	! check mass conservation of generation terms and feedback to the
1047	! large scale
1048	!
1049	do k = kts, kte
1050	if(t(i,k,j).le.t0c) then
1051	!
1052	! cloud water
1053	!
1054	value = max(qmin,qci(i,k,1))
1055	source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1056	if (source.gt.value) then
1057	factor = value/source
1058	praut(i,k) = praut(i,k)*factor
1059	pracw(i,k) = pracw(i,k)*factor
1060	psacw(i,k) = psacw(i,k)*factor
1061	endif
1062	!
1063	! cloud ice
1064	!
1065	value = max(qmin,qci(i,k,2))
1066	source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1067	if (source.gt.value) then
1068	factor = value/source
1069	psaut(i,k) = psaut(i,k)*factor
1070	psaci(i,k) = psaci(i,k)*factor
1071	pigen(i,k) = pigen(i,k)*factor
1072	pidep(i,k) = pidep(i,k)*factor
1073	endif
1074	!
1075	! rain
1076	!
1077	!
1078	value = max(qmin,qrs(i,k,1))
1079	source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1080	if (source.gt.value) then
1081	factor = value/source
1082	praut(i,k) = praut(i,k)*factor
1083	pracw(i,k) = pracw(i,k)*factor
1084	prevp(i,k) = prevp(i,k)*factor
1085	endif
1086	!
1087	! snow
1088	!
1089	value = max(qmin,qrs(i,k,2))
1090	source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1091	if (source.gt.value) then
1092	factor = value/source
1093	psdep(i,k) = psdep(i,k)*factor
1094	psaut(i,k) = psaut(i,k)*factor
1095	psaci(i,k) = psaci(i,k)*factor
1096	psacw(i,k) = psacw(i,k)*factor
1097	endif
1098	!
1099	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1100	! update
1101	q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1102	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1103	+psacw(i,k))*dtcld,0.)
1104	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1105	+prevp(i,k))*dtcld,0.)
1106	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
1107	-pigen(i,k)-pidep(i,k))*dtcld,0.)
1108	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
1109	+psaci(i,k)+psacw(i,k))*dtcld,0.)
1110	xlf = xls-xl(i,k)
1111	xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1112	-xl(i,k)prevp(i,k)-xlfpsacw(i,k)
1113	t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1114	else
1115	!
1116	! cloud water
1117	!
1118	value = max(qmin,qci(i,k,1))
1119	source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1120	if (source.gt.value) then
1121	factor = value/source
1122	praut(i,k) = praut(i,k)*factor
1123	pracw(i,k) = pracw(i,k)*factor
1124	psacw(i,k) = psacw(i,k)*factor
1125	endif
1126	!
1127	! rain
1128	!
1129	value = max(qmin,qrs(i,k,1))
1130	source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1131	if (source.gt.value) then
1132	factor = value/source
1133	praut(i,k) = praut(i,k)*factor
1134	pracw(i,k) = pracw(i,k)*factor
1135	prevp(i,k) = prevp(i,k)*factor
1136	psacw(i,k) = psacw(i,k)*factor
1137	endif
1138	!
1139	! snow
1140	!
1141	value = max(qcrmin,qrs(i,k,2))
1142	source=(-psevp(i,k))*dtcld
1143	if (source.gt.value) then
1144	factor = value/source
1145	psevp(i,k) = psevp(i,k)*factor
1146	endif
1147	work2(i,k)=-(prevp(i,k)+psevp(i,k))
1148	! update
1149	q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1150	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1151	+psacw(i,k))*dtcld,0.)
1152	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1153	+prevp(i,k) +psacw(i,k))*dtcld,0.)
1154	qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1155	xlf = xls-xl(i,k)
1156	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1157	t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1158	endif
1159	enddo
1160	!
1161	! Inline expansion for fpvs
1162	! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1163	! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1164	hsub = xls
1165	hvap = xlv0
1166	cvap = cpv
1167	ttp=t0c+0.01
1168	dldt=cvap-cliq
1169	xa=-dldt/rv
1170	xb=xa+hvap/(rv*ttp)
1171	dldti=cvap-cice
1172	xai=-dldti/rv
1173	xbi=xai+hsub/(rv*ttp)
1174	do k = kts, kte
1175	tr=ttp/t(i,k,j)
1176	logtr = log(tr)
1177	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
1178	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1179	qs(i,k,1) = max(qs(i,k,1),qmin)
1180	enddo
1181	!
1182	!----------------------------------------------------------------
1183	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1184	! if there exists additional water vapor condensated/if
1185	! evaporation of cloud water is not enough to remove subsaturation
1186	!
1187	do k = kts, kte
1188	! work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
1189	work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1190	(xl(i,k))/(rv(cpm(i,k)))*(qs(i,k,1)) &
1191	/((t(i,k,j))*(t(i,k,j)))))
1192	work2(i,k) = qci(i,k,1)+work1(i,k,1)
1193	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
1194	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1195	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1196	q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
1197	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1198	t(i,k,j) = t(i,k,j)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
1199	enddo
1200	!
1201	!
1202	!----------------------------------------------------------------
1203	! padding for small values
1204	!
1205	do k = kts, kte
1206	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1207	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1208	enddo
1209	enddo ! big loops
1210
1211	DO K=kts,kte
1212	th(i,k,j)=t(i,k,j)/pii(i,k,j)
1213	qc(i,k,j) = qci(i,k,1)
1214	qi(i,k,j) = qci(i,k,2)
1215	qr(i,k,j) = qrs(i,k,1)
1216	qqs(i,k,j) = qrs(i,k,2)
1217	ENDDO
1218
1219
1220
1221	ENDDO ! i loop
1222	enddo ! j loop
1223	!$acc end region
1224
1225	ELSE
1226
1227	!
1228	! Moved outside of accelerator region
1229	!
1230	loops = max(nint(delt/dtcldcr),1)
1231	dtcld = delt/loops
1232	if(delt.le.dtcldcr) dtcld = delt
1233
1234	!$acc region &
1235	!$acc local(t) &
1236	!$acc copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
1237	!$acc copyout(rainncv(:,:),sr(:,:)) &
1238	!$acc copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
1239	!$acc copy(th(:,:,:),q(:,:,:),rain(:,:))
1240	!$acc do &
1241	!$acc private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
1242	!$acc private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
1243	!$acc private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
1244	!$acc private(praut,psaut,prevp,psevp) &
1245	!$acc private(pracw,psacw,psaci,pcond,psmlt) &
1246	!$acc parallel
1247	do j = jts, jte
1248	!$acc do &
1249	!$acc private(numdt,mstep) &
1250	!$acc kernel vector
1251	do i = its, ite
1252	do k = kts, kte
1253	t(i,k,j)=th(i,k,j)*pii(i,k,j)
1254	qci(i,k,1) = max(qc(i,k,j),0.0)
1255	qci(i,k,2) = max(qi(i,k,j),0.0)
1256	qrs(i,k,1) = max(qr(i,k,j),0.0)
1257	qrs(i,k,2) = max(qqs(i,k,j),0.0)
1258	enddo
1259	!
1260	!----------------------------------------------------------------
1261	! latent heat for phase changes and heat capacity. neglect the
1262	! changes during microphysical process calculation
1263	! emanuel(1994)
1264	!
1265	do k = kts, kte
1266	cpm(i,k) = cpmcal(q(i,k,j))
1267	xl(i,k) = xlcal(t(i,k,j))
1268	enddo
1269	!
1270	!----------------------------------------------------------------
1271	! compute the minor time steps.
1272	!
1273	! loops = max(nint(delt/dtcldcr),1)
1274	! dtcld = delt/loops
1275	! if(delt.le.dtcldcr) dtcld = delt
1276	!
1277	do loop = 1,loops
1278	!
1279	!----------------------------------------------------------------
1280	! initialize the large scale variables
1281	!
1282	mstep = 1
1283	flgcld = .true.
1284	!
1285	do k = kts, kte
1286	denfac(i,k) = sqrt(den0/den(i,k,j))
1287	enddo
1288	! do k = kts, kte
1289	! CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
1290	! do i = its, ite
1291	! tvec1(i) = tvec1(i)*den0
1292	! enddo
1293	! CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
1294	! enddo
1295	!
1296	! Inline expansion for fpvs
1297	! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1298	! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1299	hsub = xls
1300	hvap = xlv0
1301	cvap = cpv
1302	ttp=t0c+0.01
1303	dldt=cvap-cliq
1304	xa=-dldt/rv
1305	xb=xa+hvap/(rv*ttp)
1306	dldti=cvap-cice
1307	xai=-dldti/rv
1308	xbi=xai+hsub/(rv*ttp)
1309
1310	! this is for compilers where the conditional inhibits vectorization
1311	#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
1312	do k = kts, kte
1313	if(t(i,k,j).lt.ttp) then
1314	xal = xai
1315	xbl = xbi
1316	else
1317	xal = xa
1318	xbl = xb
1319	endif
1320	tr=ttp/t(i,k,j)
1321	logtr=log(tr)
1322	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
1323	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1324	qs(i,k,1) = max(qs(i,k,1),qmin)
1325	rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
1326	qs(i,k,2)=psatexp(logtr(xal)+xbl*(1.-tr))
1327	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
1328	qs(i,k,2) = max(qs(i,k,2),qmin)
1329	rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
1330	enddo
1331	#else
1332	do k = kts, kte
1333	tr=ttp/t(i,k,j)
1334	logtr=log(tr)
1335	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
1336	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1337	qs(i,k,1) = max(qs(i,k,1),qmin)
1338	rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
1339	if(t(i,k,j).lt.ttp) then
1340	qs(i,k,2)=psatexp(logtr(xai)+xbi*(1.-tr))
1341	else
1342	qs(i,k,2)=psatexp(logtr(xa)+xb*(1.-tr))
1343	endif
1344	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
1345	qs(i,k,2) = max(qs(i,k,2),qmin)
1346	rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
1347	enddo
1348	#endif
1349	!
1350	!----------------------------------------------------------------
1351	! initialize the variables for microphysical physics
1352	!
1353	!
1354	do k = kts, kte
1355	prevp(i,k) = 0.
1356	psdep(i,k) = 0.
1357	praut(i,k) = 0.
1358	psaut(i,k) = 0.
1359	pracw(i,k) = 0.
1360	psaci(i,k) = 0.
1361	psacw(i,k) = 0.
1362	pigen(i,k) = 0.
1363	pidep(i,k) = 0.
1364	pcond(i,k) = 0.
1365	psmlt(i,k) = 0.
1366	psevp(i,k) = 0.
1367	falk(i,k,1) = 0.
1368	falk(i,k,2) = 0.
1369	fall(i,k,1) = 0.
1370	fall(i,k,2) = 0.
1371	fallc(i,k) = 0.
1372	falkc(i,k) = 0.
1373	xni(i,k) = 1.e3
1374	enddo
1375	!
1376	!----------------------------------------------------------------
1377	! compute the fallout term:
1378	! first, vertical terminal velosity for minor loops
1379	!
1380	do k = kts, kte
1381	supcol = t0c-t(i,k,j)
1382	!---------------------------------------------------------------
1383	! n0s: Intercept parameter for snow [m-4] [HDC 6]
1384	!---------------------------------------------------------------
1385	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1386	if(qrs(i,k,1).le.qcrmin)then
1387	rslope(i,k,1) = rslopermax
1388	rslopeb(i,k,1) = rsloperbmax
1389	rslope2(i,k,1) = rsloper2max
1390	rslope3(i,k,1) = rsloper3max
1391	else
1392	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
1393	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1394	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1395	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1396	endif
1397	if(qrs(i,k,2).le.qcrmin)then
1398	rslope(i,k,2) = rslopesmax
1399	rslopeb(i,k,2) = rslopesbmax
1400	rslope2(i,k,2) = rslopes2max
1401	rslope3(i,k,2) = rslopes3max
1402	else
1403	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
1404	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1405	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1406	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1407	endif
1408	!-------------------------------------------------------------
1409	! Ni: ice crystal number concentraiton [HDC 5c]
1410	!-------------------------------------------------------------
1411	! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
1412	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
1413	temp = (den(i,k,j)*max(qci(i,k,2),qmin))
1414	temp = sqrt(sqrt(temptemptemp))
1415	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1416	enddo
1417	!
1418	numdt = 1
1419	do k = kte, kts, -1
1420	work1(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)/delz(i,k,j)
1421	work1(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)/delz(i,k,j)
1422	numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
1423	if(numdt.ge.mstep) mstep = numdt
1424	enddo
1425	rmstep = 1./mstep
1426	!
1427	do n = 1, mstep
1428	k = kte
1429	! falk(i,k,1) = den(i,k,j)qrs(i,k,1)work1(i,k,1)/mstep(i)
1430	! falk(i,k,2) = den(i,k,j)qrs(i,k,2)work1(i,k,2)/mstep(i)
1431	falk(i,k,1) = den(i,k,j)qrs(i,k,1)work1(i,k,1)*rmstep
1432	falk(i,k,2) = den(i,k,j)qrs(i,k,2)work1(i,k,2)*rmstep
1433	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
1434	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
1435	! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
1436	! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
1437	dtcldden = dtcld/den(i,k,j)
1438	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
1439	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
1440	! endif
1441	do k = kte-1, kts, -1
1442	falk(i,k,1) = den(i,k,j)qrs(i,k,1)work1(i,k,1)*rmstep
1443	falk(i,k,2) = den(i,k,j)qrs(i,k,2)work1(i,k,2)*rmstep
1444	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
1445	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
1446	dtcldden = dtcld/den(i,k,j)
1447	rdelz = 1./delz(i,k,j)
1448	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
1449	delz(i,k+1,j)rdelz)*dtcldden,0.)
1450	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
1451	delz(i,k+1,j)rdelz)*dtcldden,0.)
1452	enddo
1453	do k = kte, kts, -1
1454	if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
1455	!----------------------------------------------------------------
1456	! psmlt: melting of snow [HL A33] [RH83 A25]
1457	! (T>T0: S->R)
1458	!----------------------------------------------------------------
1459	xlf = xlf0
1460	! work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
1461	work2(i,k)= (exp(log(((1.496e-6((t(i,k,j))sqrt(t(i,k,j))) &
1462	/((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5 &
1463	exp(log(t(i,k,j))(1.81))/p(i,k,j)))) &
1464	((.3333333)))/sqrt((1.496e-6((t(i,k,j)) &
1465	*sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j)))) &
1466	*sqrt(sqrt(den0/(den(i,k,j)))))
1467	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
1468	! psmlt(i,k) = xka(t(i,k,j),den(i,k,j))/xlf(t0c-t(i,k,j))pi/2. &
1469	! n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
1470	! work2(i,k)coeres)
1471	psmlt(i,k) = (1.414e3(1.496e-6((t(i,k,j))*sqrt(t(i,k,j))) &
1472	/((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j))) &
1473	/xlf(t0c-t(i,k,j))pi/2. &
1474	n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
1475	work2(i,k)coeres)
1476	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep, &
1477	-qrs(i,k,2)/mstep),0.)
1478	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
1479	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
1480	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
1481	endif
1482	enddo
1483	enddo
1484	!---------------------------------------------------------------
1485	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
1486	!---------------------------------------------------------------
1487	mstep = 1
1488	numdt = 1
1489	do k = kte, kts, -1
1490	if(qci(i,k,2).le.0.) then
1491	work2c(i,k) = 0.
1492	else
1493	xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1494	! diameter = min(dicon * sqrt(xmi),dimax)
1495	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
1496	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
1497	work2c(i,k) = work1c(i,k)/delz(i,k,j)
1498	endif
1499	numdt = max(nint(work2c(i,k)*dtcld+.5),1)
1500	if(numdt.ge.mstep) mstep = numdt
1501	enddo
1502	!
1503	do n = 1, mstep
1504	k = kte
1505	falkc(i,k) = den(i,k,j)qci(i,k,2)work2c(i,k)/mstep
1506	holdc = falkc(i,k)
1507	fallc(i,k) = fallc(i,k)+falkc(i,k)
1508	holdci = qci(i,k,2)
1509	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
1510	do k = kte-1, kts, -1
1511	falkc(i,k) = den(i,k,j)qci(i,k,2)work2c(i,k)/mstep
1512	holdc = falkc(i,k)
1513	fallc(i,k) = fallc(i,k)+falkc(i,k)
1514	holdci = qci(i,k,2)
1515	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
1516	delz(i,k+1,j)/delz(i,k,j))dtcld/den(i,k,j),0.)
1517	enddo
1518	enddo
1519	!
1520	!
1521	!----------------------------------------------------------------
1522	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
1523	!
1524	fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
1525	fallsum_qsi = fall(i,1,2)+fallc(i,1)
1526	rainncv(i,j) = 0.
1527	if(fallsum.gt.0.) then
1528	rainncv(i,j) = fallsumdelz(i,1,j)/denrdtcld*1000.
1529	rain(i,j) = fallsumdelz(i,1,j)/denrdtcld*1000. + rain(i,j)
1530	endif
1531	sr(i,j) = 0.
1532	if(fallsum.gt.0.)sr(i,j)=fallsum_qsidelz(i,kts,j)/denrdtcld*1000. &
1533	/(rainncv(i,j)+1.e-12)
1534	!
1535	!---------------------------------------------------------------
1536	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
1537	! (T>T0: I->C)
1538	!---------------------------------------------------------------
1539	do k = kts, kte
1540	supcol = t0c-t(i,k,j)
1541	xlf = xls-xl(i,k)
1542	if(supcol.lt.0.) xlf = xlf0
1543	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
1544	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
1545	t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
1546	qci(i,k,2) = 0.
1547	endif
1548	!---------------------------------------------------------------
1549	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
1550	! (T<-40C: C->I)
1551	!---------------------------------------------------------------
1552	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
1553	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
1554	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
1555	qci(i,k,1) = 0.
1556	endif
1557	!---------------------------------------------------------------
1558	! pihtf: heterogeneous freezing of cloud water [HL A44]
1559	! (T0>T>-40C: C->I)
1560	!---------------------------------------------------------------
1561	if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
1562	supcolt=min(supcol,50.)
1563	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
1564	! den(i,k,j)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
1565	pfrzdtc = min(pfrz1(exp(pfrz2supcolt)-1.) &
1566	den(i,k,j)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
1567	qci(i,k,2) = qci(i,k,2) + pfrzdtc
1568	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
1569	qci(i,k,1) = qci(i,k,1)-pfrzdtc
1570	endif
1571	!---------------------------------------------------------------
1572	! psfrz: freezing of rain water [HL A20] [LFO 45]
1573	! (T<T0, R->S)
1574	!---------------------------------------------------------------
1575	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
1576	supcolt=min(supcol,50.)
1577	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k,j) &
1578	! (exp(pfrz2supcol)-1.)rslope(i,k,1)7dtcld, &
1579	! qrs(i,k,1))
1580	temp = rslope(i,k,1)
1581	temp = temptemptemptemptemptemptemp
1582	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k,j) &
1583	(exp(pfrz2supcolt)-1.)tempdtcld, &
1584	qrs(i,k,1))
1585	qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
1586	t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
1587	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
1588	endif
1589	enddo
1590	!
1591	!----------------------------------------------------------------
1592	! rsloper: reverse of the slope parameter of the rain(m)
1593	! xka: thermal conductivity of air(jm-1s-1k-1)
1594	! work1: the thermodynamic term in the denominator associated with
1595	! heat conduction and vapor diffusion
1596	! (ry88, y93, h85)
1597	! work2: parameter associated with the ventilation effects(y93)
1598	!
1599	do k = kts, kte
1600	if(qrs(i,k,1).le.qcrmin)then
1601	rslope(i,k,1) = rslopermax
1602	rslopeb(i,k,1) = rsloperbmax
1603	rslope2(i,k,1) = rsloper2max
1604	rslope3(i,k,1) = rsloper3max
1605	else
1606	! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
1607	rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
1608	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1609	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1610	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1611	endif
1612	if(qrs(i,k,2).le.qcrmin)then
1613	rslope(i,k,2) = rslopesmax
1614	rslopeb(i,k,2) = rslopesbmax
1615	rslope2(i,k,2) = rslopes2max
1616	rslope3(i,k,2) = rslopes3max
1617	else
1618	! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
1619	rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2)) &
1620	*(den(i,k,j))))))
1621	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1622	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1623	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1624	endif
1625	enddo
1626	!
1627	do k = kts, kte
1628	! work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
1629	work1(i,k,1) = ((((den(i,k,j))(xl(i,k))(xl(i,k)))*((t(i,k,j))+120.) &
1630	(den(i,k,j)))/(1.414e3(1.496e-6((t(i,k,j))sqrt(t(i,k,j))))&
1631	(den(i,k,j))(rv(t(i,k,j))(t(i,k,j))))) &
1632	+ p(i,k,j)/((qs(i,k,1))(8.794e-5exp(log(t(i,k,j))*(1.81))))
1633	! work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
1634	work1(i,k,2) = ((((den(i,k,j))(xls)(xls))((t(i,k,j))+120.)(den(i,k,j)))&
1635	/(1.414e3(1.496e-6((t(i,k,j))sqrt(t(i,k,j))))(den(i,k,j)) &
1636	(rv(t(i,k,j))*(t(i,k,j)))) &
1637	+ p(i,k,j)/(qs(i,k,2)(8.794e-5exp(log(t(i,k,j))*(1.81)))))
1638	! work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
1639	work2(i,k) = (exp(.3333333log(((1.496e-6 ((t(i,k,j))*sqrt(t(i,k,j)))) &
1640	p(i,k,j))/(((t(i,k,j))+120.)den(i,k,j)*(8.794e-5 &
1641	exp(log(t(i,k,j))(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
1642	/sqrt((1.496e-6((t(i,k,j))sqrt(t(i,k,j)))) &
1643	/(((t(i,k,j))+120.)*den(i,k,j)))
1644	enddo
1645	!
1646	!===============================================================
1647	!
1648	! warm rain processes
1649	!
1650	! - follows the processes in RH83 and LFO except for autoconcersion
1651	!
1652	!===============================================================
1653	!
1654	do k = kts, kte
1655	supsat = max(q(i,k,j),qmin)-qs(i,k,1)
1656	satdt = supsat/dtcld
1657	!---------------------------------------------------------------
1658	! praut: auto conversion rate from cloud to rain [HDC 16]
1659	! (C->R)
1660	!---------------------------------------------------------------
1661	if(qci(i,k,1).gt.qc0) then
1662	praut(i,k) = qck1exp(log(qci(i,k,1))((7./3.)))
1663	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
1664	endif
1665	!---------------------------------------------------------------
1666	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
1667	! (C->R)
1668	!---------------------------------------------------------------
1669	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
1670	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
1671	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
1672	endif
1673	!---------------------------------------------------------------
1674	! prevp: evaporation/condensation rate of rain [HDC 14]
1675	! (V->R or R->V)
1676	!---------------------------------------------------------------
1677	if(qrs(i,k,1).gt.0.) then
1678	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
1679	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
1680	+precr2work2(i,k)coeres)/work1(i,k,1)
1681	if(prevp(i,k).lt.0.) then
1682	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1683	prevp(i,k) = max(prevp(i,k),satdt/2)
1684	else
1685	prevp(i,k) = min(prevp(i,k),satdt/2)
1686	endif
1687	endif
1688	enddo
1689	!
1690	!===============================================================
1691	!
1692	! cold rain processes
1693	!
1694	! - follows the revised ice microphysics processes in HDC
1695	! - the processes same as in RH83 and RH84 and LFO behave
1696	! following ice crystal hapits defined in HDC, inclduing
1697	! intercept parameter for snow (n0s), ice crystal number
1698	! concentration (ni), ice nuclei number concentration
1699	! (n0i), ice diameter (d)
1700	!
1701	!===============================================================
1702	!
1703	rdtcld = 1./dtcld
1704	do k = kts, kte
1705	supcol = t0c-t(i,k,j)
1706	supsat = max(q(i,k,j),qmin)-qs(i,k,2)
1707	satdt = supsat/dtcld
1708	ifsat = 0
1709	!-------------------------------------------------------------
1710	! Ni: ice crystal number concentraiton [HDC 5c]
1711	!-------------------------------------------------------------
1712	! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
1713	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
1714	temp = (den(i,k,j)*max(qci(i,k,2),qmin))
1715	temp = sqrt(sqrt(temptemptemp))
1716	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1717	eacrs = exp(0.07*(-supcol))
1718	!
1719	if(supcol.gt.0) then
1720	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
1721	xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1722	diameter = min(dicon * sqrt(xmi),dimax)
1723	vt2i = 1.49e4diameter*1.31
1724	vt2s = pvtsrslopeb(i,k,2)denfac(i,k)
1725	!-------------------------------------------------------------
1726	! psaci: Accretion of cloud ice by rain [HDC 10]
1727	! (T<T0: I->S)
1728	!-------------------------------------------------------------
1729	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
1730	+diameter*2rslope(i,k,2)
1731	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
1732	abs(vt2s-vt2i)acrfac/4.
1733	endif
1734	endif
1735	!-------------------------------------------------------------
1736	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
1737	! (T<T0: C->S, and T>=T0: C->R)
1738	!-------------------------------------------------------------
1739	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
1740	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2) &
1741	rslopeb(i,k,2)qci(i,k,1)*denfac(i,k) &
1742	! ,qci(i,k,1)/dtcld)
1743	,qci(i,k,1)*rdtcld)
1744	endif
1745	if(supcol .gt. 0) then
1746	!-------------------------------------------------------------
1747	! pidep: Deposition/Sublimation rate of ice [HDC 9]
1748	! (T<T0: V->I or I->V)
1749	!-------------------------------------------------------------
1750	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
1751	xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1752	diameter = dicon * sqrt(xmi)
1753	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1754	supice = satdt-prevp(i,k)
1755	if(pidep(i,k).lt.0.) then
1756	! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1757	! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1758	pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
1759	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
1760	else
1761	! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1762	pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
1763	endif
1764	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1765	endif
1766	!-------------------------------------------------------------
1767	! psdep: deposition/sublimation rate of snow [HDC 14]
1768	! (V->S or S->V)
1769	!-------------------------------------------------------------
1770	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1771	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
1772	psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
1773	(precs1rslope2(i,k,2)+precs2 &
1774	work2(i,k)coeres)/work1(i,k,2)
1775	supice = satdt-prevp(i,k)-pidep(i,k)
1776	if(psdep(i,k).lt.0.) then
1777	! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1778	! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1779	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1780	psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1781	else
1782	! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1783	psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1784	endif
1785	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1786	ifsat = 1
1787	endif
1788	!-------------------------------------------------------------
1789	! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1790	! (T<T0: V->I)
1791	!-------------------------------------------------------------
1792	if(supsat.gt.0.and.ifsat.ne.1) then
1793	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1794	xni0 = 1.e3exp(0.1supcol)
1795	roqi0 = 4.92e-11exp(log(xni0)(1.33))
1796	pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.)) &
1797	! /dtcld)
1798	*rdtcld)
1799	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1800	endif
1801	!
1802	!-------------------------------------------------------------
1803	! psaut: conversion(aggregation) of ice to snow [HDC 12]
1804	! (T<T0: I->S)
1805	!-------------------------------------------------------------
1806	if(qci(i,k,2).gt.0.) then
1807	qimax = roqimax/den(i,k,j)
1808	! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1809	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1810	endif
1811	endif
1812	!-------------------------------------------------------------
1813	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1814	! (T>T0: S->V)
1815	!-------------------------------------------------------------
1816	if(supcol.lt.0.) then
1817	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
1818	psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1819	! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1820	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1821	endif
1822	enddo
1823	!
1824	!
1825	!----------------------------------------------------------------
1826	! check mass conservation of generation terms and feedback to the
1827	! large scale
1828	!
1829	do k = kts, kte
1830	if(t(i,k,j).le.t0c) then
1831	!
1832	! cloud water
1833	!
1834	value = max(qmin,qci(i,k,1))
1835	source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1836	if (source.gt.value) then
1837	factor = value/source
1838	praut(i,k) = praut(i,k)*factor
1839	pracw(i,k) = pracw(i,k)*factor
1840	psacw(i,k) = psacw(i,k)*factor
1841	endif
1842	!
1843	! cloud ice
1844	!
1845	value = max(qmin,qci(i,k,2))
1846	source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1847	if (source.gt.value) then
1848	factor = value/source
1849	psaut(i,k) = psaut(i,k)*factor
1850	psaci(i,k) = psaci(i,k)*factor
1851	pigen(i,k) = pigen(i,k)*factor
1852	pidep(i,k) = pidep(i,k)*factor
1853	endif
1854	!
1855	! rain
1856	!
1857	!
1858	value = max(qmin,qrs(i,k,1))
1859	source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1860	if (source.gt.value) then
1861	factor = value/source
1862	praut(i,k) = praut(i,k)*factor
1863	pracw(i,k) = pracw(i,k)*factor
1864	prevp(i,k) = prevp(i,k)*factor
1865	endif
1866	!
1867	! snow
1868	!
1869	value = max(qmin,qrs(i,k,2))
1870	source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1871	if (source.gt.value) then
1872	factor = value/source
1873	psdep(i,k) = psdep(i,k)*factor
1874	psaut(i,k) = psaut(i,k)*factor
1875	psaci(i,k) = psaci(i,k)*factor
1876	psacw(i,k) = psacw(i,k)*factor
1877	endif
1878	!
1879	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1880	! update
1881	q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1882	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1883	+psacw(i,k))*dtcld,0.)
1884	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1885	+prevp(i,k))*dtcld,0.)
1886	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
1887	-pigen(i,k)-pidep(i,k))*dtcld,0.)
1888	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
1889	+psaci(i,k)+psacw(i,k))*dtcld,0.)
1890	xlf = xls-xl(i,k)
1891	xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1892	-xl(i,k)prevp(i,k)-xlfpsacw(i,k)
1893	t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1894	else
1895	!
1896	! cloud water
1897	!
1898	value = max(qmin,qci(i,k,1))
1899	source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1900	if (source.gt.value) then
1901	factor = value/source
1902	praut(i,k) = praut(i,k)*factor
1903	pracw(i,k) = pracw(i,k)*factor
1904	psacw(i,k) = psacw(i,k)*factor
1905	endif
1906	!
1907	! rain
1908	!
1909	value = max(qmin,qrs(i,k,1))
1910	source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1911	if (source.gt.value) then
1912	factor = value/source
1913	praut(i,k) = praut(i,k)*factor
1914	pracw(i,k) = pracw(i,k)*factor
1915	prevp(i,k) = prevp(i,k)*factor
1916	psacw(i,k) = psacw(i,k)*factor
1917	endif
1918	!
1919	! snow
1920	!
1921	value = max(qcrmin,qrs(i,k,2))
1922	source=(-psevp(i,k))*dtcld
1923	if (source.gt.value) then
1924	factor = value/source
1925	psevp(i,k) = psevp(i,k)*factor
1926	endif
1927	work2(i,k)=-(prevp(i,k)+psevp(i,k))
1928	! update
1929	q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1930	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1931	+psacw(i,k))*dtcld,0.)
1932	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1933	+prevp(i,k) +psacw(i,k))*dtcld,0.)
1934	qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1935	xlf = xls-xl(i,k)
1936	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1937	t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1938	endif
1939	enddo
1940	!
1941	! Inline expansion for fpvs
1942	! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1943	! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1944	hsub = xls
1945	hvap = xlv0
1946	cvap = cpv
1947	ttp=t0c+0.01
1948	dldt=cvap-cliq
1949	xa=-dldt/rv
1950	xb=xa+hvap/(rv*ttp)
1951	dldti=cvap-cice
1952	xai=-dldti/rv
1953	xbi=xai+hsub/(rv*ttp)
1954	do k = kts, kte
1955	tr=ttp/t(i,k,j)
1956	logtr = log(tr)
1957	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
1958	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1959	qs(i,k,1) = max(qs(i,k,1),qmin)
1960	enddo
1961	!
1962	!----------------------------------------------------------------
1963	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1964	! if there exists additional water vapor condensated/if
1965	! evaporation of cloud water is not enough to remove subsaturation
1966	!
1967	do k = kts, kte
1968	! work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
1969	work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1970	(xl(i,k))/(rv(cpm(i,k)))*(qs(i,k,1)) &
1971	/((t(i,k,j))*(t(i,k,j)))))
1972	work2(i,k) = qci(i,k,1)+work1(i,k,1)
1973	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
1974	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1975	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1976	q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
1977	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1978	t(i,k,j) = t(i,k,j)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
1979	enddo
1980	!
1981	!
1982	!----------------------------------------------------------------
1983	! padding for small values
1984	!
1985	do k = kts, kte
1986	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1987	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1988	enddo
1989	enddo ! big loops
1990
1991	DO K=kts,kte
1992	th(i,k,j)=t(i,k,j)/pii(i,k,j)
1993	qc(i,k,j) = qci(i,k,1)
1994	qi(i,k,j) = qci(i,k,2)
1995	qr(i,k,j) = qrs(i,k,1)
1996	qqs(i,k,j) = qrs(i,k,2)
1997	ENDDO
1998
1999	ENDDO ! i loop
2000	enddo ! j loop
2001	!$acc end region
2002
2003	ENDIF
2004
2005	END SUBROUTINE wsm52d
2006
2007
2008	#else
2009
2010	!===================================================================
2011	!
2012	SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz &
2013	,delt,g, cpd, cpv, rd, rv, t0c &
2014	,ep1, ep2, qmin &
2015	,XLS, XLV0, XLF0, den0, denr &
2016	,cliq,cice,psat &
2017	,lat &
2018	,rain,rainncv &
2019	,sr &
2020	,ids,ide, jds,jde, kds,kde &
2021	,ims,ime, jms,jme, kms,kme &
2022	,its,ite, jts,jte, kts,kte &
2023	,snow,snowncv &
2024	)
2025	!-------------------------------------------------------------------
2026	IMPLICIT NONE
2027	!-------------------------------------------------------------------
2028	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
2029	ims,ime, jms,jme, kms,kme , &
2030	its,ite, jts,jte, kts,kte, &
2031	lat
2032	REAL, DIMENSION( its:ite , kts:kte ), &
2033	INTENT(INOUT) :: &
2034	t
2035	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
2036	INTENT(INOUT) :: &
2037	qci, &
2038	qrs
2039	REAL, DIMENSION( ims:ime , kms:kme ), &
2040	INTENT(INOUT) :: &
2041	q
2042	REAL, DIMENSION( ims:ime , kms:kme ), &
2043	INTENT(IN ) :: &
2044	den, &
2045	p, &
2046	delz
2047	REAL, INTENT(IN ) :: delt, &
2048	g, &
2049	cpd, &
2050	cpv, &
2051	t0c, &
2052	den0, &
2053	rd, &
2054	rv, &
2055	ep1, &
2056	ep2, &
2057	qmin, &
2058	XLS, &
2059	XLV0, &
2060	XLF0, &
2061	cliq, &
2062	cice, &
2063	psat, &
2064	denr
2065	REAL, DIMENSION( ims:ime ), &
2066	INTENT(INOUT) :: rain, &
2067	rainncv, &
2068	sr
2069
2070	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
2071	INTENT(INOUT) :: snow, &
2072	snowncv
2073
2074	! LOCAL VAR
2075	REAL, DIMENSION( its:ite , kts:kte , 2) :: &
2076	rh, &
2077	qs, &
2078	rslope, &
2079	rslope2, &
2080	rslope3, &
2081	rslopeb, &
2082	falk, &
2083	fall, &
2084	work1
2085
2086	REAL, DIMENSION( its:ite , kts:kte ) :: &
2087	falkc, &
2088	fallc, &
2089	xl, &
2090	cpm, &
2091	denfac, &
2092	xni, &
2093	n0sfac, &
2094	work2, &
2095	work1c, &
2096	work2c
2097
2098	REAL, DIMENSION( its:ite , kts:kte ) :: &
2099	pigen, &
2100	pidep, &
2101	psdep, &
2102	praut, &
2103	psaut, &
2104	prevp, &
2105	psevp, &
2106	pracw, &
2107	psacw, &
2108	psaci, &
2109	pcond, &
2110	psmlt
2111	INTEGER, DIMENSION( its:ite ) :: &
2112	mstep, &
2113	numdt
2114	REAL, DIMENSION(its:ite) :: rmstep
2115	REAL dtcldden, rdelz, rdtcld
2116	LOGICAL, DIMENSION( its:ite ) :: flgcld
2117
2118	#define WSM_NO_CONDITIONAL_IN_VECTOR
2119	#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
2120	REAL, DIMENSION(its:ite) :: xal, xbl
2121	#endif
2122
2123	REAL :: pi, &
2124	cpmcal, xlcal, lamdar, lamdas, diffus, &
2125	viscos, xka, venfac, conden, diffac, &
2126	x, y, z, a, b, c, d, e, &
2127	qdt, holdrr, holdrs, supcol, supcolt, pvt, &
2128	coeres, supsat, dtcld, xmi, eacrs, satdt, &
2129	vt2i,vt2s,acrfac, &
2130	qimax, diameter, xni0, roqi0, &
2131	fallsum, fallsum_qsi, xlwork2, factor, source, &
2132	value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
2133	! variables for optimization
2134	REAL, DIMENSION( its:ite ) :: tvec1
2135	REAL :: temp
2136	INTEGER :: i, j, k, mstepmax, &
2137	iprt, latd, lond, loop, loops, ifsat, n
2138	! Temporaries used for inlining fpvs function
2139	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
2140	REAL :: logtr
2141	!
2142	!=================================================================
2143	! compute internal functions
2144	!
2145	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
2146	xlcal(x) = xlv0-xlv1*(x-t0c)
2147	!----------------------------------------------------------------
2148	! size distributions: (x=mixing ratio, y=air density):
2149	! valid for mixing ratio > 1.e-9 kg/kg.
2150	!
2151	! Optimizatin : A*B => exp(log(A)(B))
2152	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
2153	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
2154	!
2155	!----------------------------------------------------------------
2156	! diffus: diffusion coefficient of the water vapor
2157	! viscos: kinematic viscosity(m2s-1)
2158	! diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
2159	! viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
2160	! xka(x,y) = 1.414e3viscos(x,y)y
2161	! diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
2162	! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
2163	! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
2164	! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
2165	!
2166	!
2167	pi = 4. * atan(1.)
2168	!
2169	!----------------------------------------------------------------
2170	! paddint 0 for negative values generated by dynamics
2171	!
2172	do k = kts, kte
2173	do i = its, ite
2174	qci(i,k,1) = max(qci(i,k,1),0.0)
2175	qrs(i,k,1) = max(qrs(i,k,1),0.0)
2176	qci(i,k,2) = max(qci(i,k,2),0.0)
2177	qrs(i,k,2) = max(qrs(i,k,2),0.0)
2178	enddo
2179	enddo
2180	!
2181	!----------------------------------------------------------------
2182	! latent heat for phase changes and heat capacity. neglect the
2183	! changes during microphysical process calculation
2184	! emanuel(1994)
2185	!
2186	do k = kts, kte
2187	do i = its, ite
2188	cpm(i,k) = cpmcal(q(i,k))
2189	xl(i,k) = xlcal(t(i,k))
2190	enddo
2191	enddo
2192	!
2193	!----------------------------------------------------------------
2194	! compute the minor time steps.
2195	!
2196	loops = max(nint(delt/dtcldcr),1)
2197	dtcld = delt/loops
2198	if(delt.le.dtcldcr) dtcld = delt
2199	!
2200	do loop = 1,loops
2201	!
2202	!----------------------------------------------------------------
2203	! initialize the large scale variables
2204	!
2205	do i = its, ite
2206	mstep(i) = 1
2207	flgcld(i) = .true.
2208	enddo
2209	!
2210	! do k = kts, kte
2211	! do i = its, ite
2212	! denfac(i,k) = sqrt(den0/den(i,k))
2213	! enddo
2214	! enddo
2215	do k = kts, kte
2216	CALL VREC( tvec1(its), den(its,k), ite-its+1)
2217	do i = its, ite
2218	tvec1(i) = tvec1(i)*den0
2219	enddo
2220	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
2221	enddo
2222	!
2223	! Inline expansion for fpvs
2224	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2225	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2226	hsub = xls
2227	hvap = xlv0
2228	cvap = cpv
2229	ttp=t0c+0.01
2230	dldt=cvap-cliq
2231	xa=-dldt/rv
2232	xb=xa+hvap/(rv*ttp)
2233	dldti=cvap-cice
2234	xai=-dldti/rv
2235	xbi=xai+hsub/(rv*ttp)
2236
2237	! this is for compilers where the conditional inhibits vectorization
2238	#ifdef WSM_NO_CONDITIONAL_IN_VECTOR
2239	do k = kts, kte
2240	do i = its, ite
2241	if(t(i,k).lt.ttp) then
2242	xal(i) = xai
2243	xbl(i) = xbi
2244	else
2245	xal(i) = xa
2246	xbl(i) = xb
2247	endif
2248	enddo
2249	do i = its, ite
2250	tr=ttp/t(i,k)
2251	logtr=log(tr)
2252	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
2253	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2254	qs(i,k,1) = max(qs(i,k,1),qmin)
2255	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2256	qs(i,k,2)=psatexp(logtr(xal(i))+xbl(i)*(1.-tr))
2257	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2258	qs(i,k,2) = max(qs(i,k,2),qmin)
2259	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
2260	enddo
2261	enddo
2262	#else
2263	do k = kts, kte
2264	do i = its, ite
2265	tr=ttp/t(i,k)
2266	logtr=log(tr)
2267	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
2268	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2269	qs(i,k,1) = max(qs(i,k,1),qmin)
2270	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2271	if(t(i,k).lt.ttp) then
2272	qs(i,k,2)=psatexp(logtr(xai)+xbi*(1.-tr))
2273	else
2274	qs(i,k,2)=psatexp(logtr(xa)+xb*(1.-tr))
2275	endif
2276	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2277	qs(i,k,2) = max(qs(i,k,2),qmin)
2278	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
2279	enddo
2280	enddo
2281	#endif
2282	!
2283	!----------------------------------------------------------------
2284	! initialize the variables for microphysical physics
2285	!
2286	!
2287	do k = kts, kte
2288	do i = its, ite
2289	prevp(i,k) = 0.
2290	psdep(i,k) = 0.
2291	praut(i,k) = 0.
2292	psaut(i,k) = 0.
2293	pracw(i,k) = 0.
2294	psaci(i,k) = 0.
2295	psacw(i,k) = 0.
2296	pigen(i,k) = 0.
2297	pidep(i,k) = 0.
2298	pcond(i,k) = 0.
2299	psmlt(i,k) = 0.
2300	psevp(i,k) = 0.
2301	falk(i,k,1) = 0.
2302	falk(i,k,2) = 0.
2303	fall(i,k,1) = 0.
2304	fall(i,k,2) = 0.
2305	fallc(i,k) = 0.
2306	falkc(i,k) = 0.
2307	xni(i,k) = 1.e3
2308	enddo
2309	enddo
2310	!
2311	!----------------------------------------------------------------
2312	! compute the fallout term:
2313	! first, vertical terminal velosity for minor loops
2314	!
2315	do k = kts, kte
2316	do i = its, ite
2317	supcol = t0c-t(i,k)
2318	!---------------------------------------------------------------
2319	! n0s: Intercept parameter for snow [m-4] [HDC 6]
2320	!---------------------------------------------------------------
2321	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2322	if(qrs(i,k,1).le.qcrmin)then
2323	rslope(i,k,1) = rslopermax
2324	rslopeb(i,k,1) = rsloperbmax
2325	rslope2(i,k,1) = rsloper2max
2326	rslope3(i,k,1) = rsloper3max
2327	else
2328	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
2329	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
2330	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2331	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2332	endif
2333	if(qrs(i,k,2).le.qcrmin)then
2334	rslope(i,k,2) = rslopesmax
2335	rslopeb(i,k,2) = rslopesbmax
2336	rslope2(i,k,2) = rslopes2max
2337	rslope3(i,k,2) = rslopes3max
2338	else
2339	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2340	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
2341	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2342	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2343	endif
2344	!-------------------------------------------------------------
2345	! Ni: ice crystal number concentraiton [HDC 5c]
2346	!-------------------------------------------------------------
2347	! xni(i,k) = min(max(5.38e7*(den(i,k) &
2348	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
2349	temp = (den(i,k)*max(qci(i,k,2),qmin))
2350	temp = sqrt(sqrt(temptemptemp))
2351	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2352	enddo
2353	enddo
2354	!
2355	mstepmax = 1
2356	numdt = 1
2357	do k = kte, kts, -1
2358	do i = its, ite
2359	work1(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)/delz(i,k)
2360	work1(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)/delz(i,k)
2361	numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
2362	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
2363	enddo
2364	enddo
2365	do i = its, ite
2366	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
2367	rmstep(i) = 1./mstep(i)
2368	enddo
2369	!
2370	do n = 1, mstepmax
2371	k = kte
2372	do i = its, ite
2373	if(n.le.mstep(i)) then
2374	! falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
2375	! falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)/mstep(i)
2376	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)*rmstep(i)
2377	falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)*rmstep(i)
2378	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
2379	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
2380	! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
2381	! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
2382	dtcldden = dtcld/den(i,k)
2383	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
2384	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
2385	endif
2386	enddo
2387	do k = kte-1, kts, -1
2388	do i = its, ite
2389	if(n.le.mstep(i)) then
2390	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)*rmstep(i)
2391	falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)*rmstep(i)
2392	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
2393	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
2394	dtcldden = dtcld/den(i,k)
2395	rdelz = 1./delz(i,k)
2396	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
2397	delz(i,k+1)rdelz)*dtcldden,0.)
2398	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
2399	delz(i,k+1)rdelz)*dtcldden,0.)
2400	endif
2401	enddo
2402	enddo
2403	do k = kte, kts, -1
2404	do i = its, ite
2405	if(n.le.mstep(i)) then
2406	if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
2407	!----------------------------------------------------------------
2408	! psmlt: melting of snow [HL A33] [RH83 A25]
2409	! (T>T0: S->R)
2410	!----------------------------------------------------------------
2411	xlf = xlf0
2412	! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
2413	work2(i,k)= (exp(log(((1.496e-6((t(i,k))sqrt(t(i,k))) &
2414	/((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
2415	exp(log(t(i,k))(1.81))/p(i,k)))) &
2416	((.3333333)))/sqrt((1.496e-6((t(i,k)) &
2417	*sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
2418	*sqrt(sqrt(den0/(den(i,k)))))
2419	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
2420	! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
2421	! n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
2422	! work2(i,k)coeres)
2423	psmlt(i,k) = (1.414e3(1.496e-6((t(i,k))*sqrt(t(i,k))) &
2424	/((t(i,k))+120.)/(den(i,k)) )*(den(i,k))) &
2425	/xlf(t0c-t(i,k))pi/2. &
2426	n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
2427	work2(i,k)coeres)
2428	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
2429	-qrs(i,k,2)/mstep(i)),0.)
2430	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
2431	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
2432	t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
2433	endif
2434	endif
2435	enddo
2436	enddo
2437	enddo
2438	!---------------------------------------------------------------
2439	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
2440	!---------------------------------------------------------------
2441	mstepmax = 1
2442	mstep = 1
2443	numdt = 1
2444	do k = kte, kts, -1
2445	do i = its, ite
2446	if(qci(i,k,2).le.0.) then
2447	work2c(i,k) = 0.
2448	else
2449	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2450	! diameter = min(dicon * sqrt(xmi),dimax)
2451	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
2452	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
2453	work2c(i,k) = work1c(i,k)/delz(i,k)
2454	endif
2455	numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
2456	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
2457	enddo
2458	enddo
2459	do i = its, ite
2460	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
2461	enddo
2462	!
2463	do n = 1, mstepmax
2464	k = kte
2465	do i = its, ite
2466	if(n.le.mstep(i)) then
2467	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
2468	holdc = falkc(i,k)
2469	fallc(i,k) = fallc(i,k)+falkc(i,k)
2470	holdci = qci(i,k,2)
2471	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
2472	endif
2473	enddo
2474	do k = kte-1, kts, -1
2475	do i = its, ite
2476	if(n.le.mstep(i)) then
2477	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
2478	holdc = falkc(i,k)
2479	fallc(i,k) = fallc(i,k)+falkc(i,k)
2480	holdci = qci(i,k,2)
2481	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
2482	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
2483	endif
2484	enddo
2485	enddo
2486	enddo
2487	!
2488	!
2489	!----------------------------------------------------------------
2490	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
2491	!
2492	do i = its, ite
2493	fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
2494	fallsum_qsi = fall(i,1,2)+fallc(i,1)
2495	rainncv(i) = 0.
2496	if(fallsum.gt.0.) then
2497	rainncv(i) = fallsumdelz(i,1)/denrdtcld*1000.
2498	rain(i) = fallsumdelz(i,1)/denrdtcld*1000. + rain(i)
2499	endif
2500	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
2501	snowncv(i,lat) = 0.
2502	if(fallsum_qsi.gt.0.) then
2503	snowncv(i,lat) = fallsum_qsidelz(i,kts)/denrdtcld*1000.
2504	snow(i,lat) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i,lat)
2505	endif
2506	ENDIF
2507	sr(i) = 0.
2508	if(fallsum.gt.0.)sr(i)=fallsum_qsidelz(i,kts)/denrdtcld*1000. &
2509	/(rainncv(i)+1.e-12)
2510	enddo
2511	!
2512	!---------------------------------------------------------------
2513	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
2514	! (T>T0: I->C)
2515	!---------------------------------------------------------------
2516	do k = kts, kte
2517	do i = its, ite
2518	supcol = t0c-t(i,k)
2519	xlf = xls-xl(i,k)
2520	if(supcol.lt.0.) xlf = xlf0
2521	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
2522	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
2523	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
2524	qci(i,k,2) = 0.
2525	endif
2526	!---------------------------------------------------------------
2527	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
2528	! (T<-40C: C->I)
2529	!---------------------------------------------------------------
2530	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
2531	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
2532	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
2533	qci(i,k,1) = 0.
2534	endif
2535	!---------------------------------------------------------------
2536	! pihtf: heterogeneous freezing of cloud water [HL A44]
2537	! (T0>T>-40C: C->I)
2538	!---------------------------------------------------------------
2539	if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
2540	supcolt=min(supcol,50.)
2541	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
2542	! den(i,k)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
2543	pfrzdtc = min(pfrz1(exp(pfrz2supcolt)-1.) &
2544	den(i,k)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
2545	qci(i,k,2) = qci(i,k,2) + pfrzdtc
2546	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
2547	qci(i,k,1) = qci(i,k,1)-pfrzdtc
2548	endif
2549	!---------------------------------------------------------------
2550	! psfrz: freezing of rain water [HL A20] [LFO 45]
2551	! (T<T0, R->S)
2552	!---------------------------------------------------------------
2553	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
2554	supcolt=min(supcol,50.)
2555	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k) &
2556	! (exp(pfrz2supcol)-1.)rslope(i,k,1)7dtcld, &
2557	! qrs(i,k,1))
2558	temp = rslope(i,k,1)
2559	temp = temptemptemptemptemptemptemp
2560	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k) &
2561	(exp(pfrz2supcolt)-1.)tempdtcld, &
2562	qrs(i,k,1))
2563	qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
2564	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
2565	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
2566	endif
2567	enddo
2568	enddo
2569	!
2570	!----------------------------------------------------------------
2571	! rsloper: reverse of the slope parameter of the rain(m)
2572	! xka: thermal conductivity of air(jm-1s-1k-1)
2573	! work1: the thermodynamic term in the denominator associated with
2574	! heat conduction and vapor diffusion
2575	! (ry88, y93, h85)
2576	! work2: parameter associated with the ventilation effects(y93)
2577	!
2578	do k = kts, kte
2579	do i = its, ite
2580	if(qrs(i,k,1).le.qcrmin)then
2581	rslope(i,k,1) = rslopermax
2582	rslopeb(i,k,1) = rsloperbmax
2583	rslope2(i,k,1) = rsloper2max
2584	rslope3(i,k,1) = rsloper3max
2585	else
2586	! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
2587	rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
2588	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
2589	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2590	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2591	endif
2592	if(qrs(i,k,2).le.qcrmin)then
2593	rslope(i,k,2) = rslopesmax
2594	rslopeb(i,k,2) = rslopesbmax
2595	rslope2(i,k,2) = rslopes2max
2596	rslope3(i,k,2) = rslopes3max
2597	else
2598	! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2599	rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2)) &
2600	*(den(i,k))))))
2601	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
2602	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2603	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2604	endif
2605	enddo
2606	enddo
2607	!
2608	do k = kts, kte
2609	do i = its, ite
2610	! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
2611	work1(i,k,1) = ((((den(i,k))(xl(i,k))(xl(i,k)))*((t(i,k))+120.) &
2612	(den(i,k)))/(1.414e3(1.496e-6((t(i,k))sqrt(t(i,k))))&
2613	(den(i,k))(rv(t(i,k))(t(i,k))))) &
2614	+ p(i,k)/((qs(i,k,1))(8.794e-5exp(log(t(i,k))*(1.81))))
2615	! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
2616	work1(i,k,2) = ((((den(i,k))(xls)(xls))((t(i,k))+120.)(den(i,k)))&
2617	/(1.414e3(1.496e-6((t(i,k))sqrt(t(i,k))))(den(i,k)) &
2618	(rv(t(i,k))*(t(i,k)))) &
2619	+ p(i,k)/(qs(i,k,2)(8.794e-5exp(log(t(i,k))*(1.81)))))
2620	! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
2621	work2(i,k) = (exp(.3333333log(((1.496e-6 ((t(i,k))*sqrt(t(i,k)))) &
2622	p(i,k))/(((t(i,k))+120.)den(i,k)*(8.794e-5 &
2623	exp(log(t(i,k))(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
2624	/sqrt((1.496e-6((t(i,k))sqrt(t(i,k)))) &
2625	/(((t(i,k))+120.)*den(i,k)))
2626	enddo
2627	enddo
2628	!
2629	!===============================================================
2630	!
2631	! warm rain processes
2632	!
2633	! - follows the processes in RH83 and LFO except for autoconcersion
2634	!
2635	!===============================================================
2636	!
2637	do k = kts, kte
2638	do i = its, ite
2639	supsat = max(q(i,k),qmin)-qs(i,k,1)
2640	satdt = supsat/dtcld
2641	!---------------------------------------------------------------
2642	! praut: auto conversion rate from cloud to rain [HDC 16]
2643	! (C->R)
2644	!---------------------------------------------------------------
2645	if(qci(i,k,1).gt.qc0) then
2646	praut(i,k) = qck1exp(log(qci(i,k,1))((7./3.)))
2647	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
2648	endif
2649	!---------------------------------------------------------------
2650	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
2651	! (C->R)
2652	!---------------------------------------------------------------
2653	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
2654	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
2655	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
2656	endif
2657	!---------------------------------------------------------------
2658	! prevp: evaporation/condensation rate of rain [HDC 14]
2659	! (V->R or R->V)
2660	!---------------------------------------------------------------
2661	if(qrs(i,k,1).gt.0.) then
2662	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
2663	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
2664	+precr2work2(i,k)coeres)/work1(i,k,1)
2665	if(prevp(i,k).lt.0.) then
2666	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
2667	prevp(i,k) = max(prevp(i,k),satdt/2)
2668	else
2669	prevp(i,k) = min(prevp(i,k),satdt/2)
2670	endif
2671	endif
2672	enddo
2673	enddo
2674	!
2675	!===============================================================
2676	!
2677	! cold rain processes
2678	!
2679	! - follows the revised ice microphysics processes in HDC
2680	! - the processes same as in RH83 and RH84 and LFO behave
2681	! following ice crystal hapits defined in HDC, inclduing
2682	! intercept parameter for snow (n0s), ice crystal number
2683	! concentration (ni), ice nuclei number concentration
2684	! (n0i), ice diameter (d)
2685	!
2686	!===============================================================
2687	!
2688	rdtcld = 1./dtcld
2689	do k = kts, kte
2690	do i = its, ite
2691	supcol = t0c-t(i,k)
2692	supsat = max(q(i,k),qmin)-qs(i,k,2)
2693	satdt = supsat/dtcld
2694	ifsat = 0
2695	!-------------------------------------------------------------
2696	! Ni: ice crystal number concentraiton [HDC 5c]
2697	!-------------------------------------------------------------
2698	! xni(i,k) = min(max(5.38e7*(den(i,k) &
2699	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
2700	temp = (den(i,k)*max(qci(i,k,2),qmin))
2701	temp = sqrt(sqrt(temptemptemp))
2702	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2703	eacrs = exp(0.07*(-supcol))
2704	!
2705	if(supcol.gt.0) then
2706	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
2707	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2708	diameter = min(dicon * sqrt(xmi),dimax)
2709	vt2i = 1.49e4diameter*1.31
2710	vt2s = pvtsrslopeb(i,k,2)denfac(i,k)
2711	!-------------------------------------------------------------
2712	! psaci: Accretion of cloud ice by rain [HDC 10]
2713	! (T<T0: I->S)
2714	!-------------------------------------------------------------
2715	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
2716	+diameter*2rslope(i,k,2)
2717	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
2718	abs(vt2s-vt2i)acrfac/4.
2719	endif
2720	endif
2721	!-------------------------------------------------------------
2722	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
2723	! (T<T0: C->S, and T>=T0: C->R)
2724	!-------------------------------------------------------------
2725	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
2726	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2) &
2727	rslopeb(i,k,2)qci(i,k,1)*denfac(i,k) &
2728	! ,qci(i,k,1)/dtcld)
2729	,qci(i,k,1)*rdtcld)
2730	endif
2731	if(supcol .gt. 0) then
2732	!-------------------------------------------------------------
2733	! pidep: Deposition/Sublimation rate of ice [HDC 9]
2734	! (T<T0: V->I or I->V)
2735	!-------------------------------------------------------------
2736	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
2737	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2738	diameter = dicon * sqrt(xmi)
2739	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
2740	supice = satdt-prevp(i,k)
2741	if(pidep(i,k).lt.0.) then
2742	! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
2743	! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
2744	pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
2745	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
2746	else
2747	! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
2748	pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
2749	endif
2750	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
2751	endif
2752	!-------------------------------------------------------------
2753	! psdep: deposition/sublimation rate of snow [HDC 14]
2754	! (V->S or S->V)
2755	!-------------------------------------------------------------
2756	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
2757	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
2758	psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
2759	(precs1rslope2(i,k,2)+precs2 &
2760	work2(i,k)coeres)/work1(i,k,2)
2761	supice = satdt-prevp(i,k)-pidep(i,k)
2762	if(psdep(i,k).lt.0.) then
2763	! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
2764	! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
2765	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
2766	psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
2767	else
2768	! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
2769	psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
2770	endif
2771	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
2772	ifsat = 1
2773	endif
2774	!-------------------------------------------------------------
2775	! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
2776	! (T<T0: V->I)
2777	!-------------------------------------------------------------
2778	if(supsat.gt.0.and.ifsat.ne.1) then
2779	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
2780	xni0 = 1.e3exp(0.1supcol)
2781	roqi0 = 4.92e-11exp(log(xni0)(1.33))
2782	pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
2783	! /dtcld)
2784	*rdtcld)
2785	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
2786	endif
2787	!
2788	!-------------------------------------------------------------
2789	! psaut: conversion(aggregation) of ice to snow [HDC 12]
2790	! (T<T0: I->S)
2791	!-------------------------------------------------------------
2792	if(qci(i,k,2).gt.0.) then
2793	qimax = roqimax/den(i,k)
2794	! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
2795	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
2796	endif
2797	endif
2798	!-------------------------------------------------------------
2799	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
2800	! (T>T0: S->V)
2801	!-------------------------------------------------------------
2802	if(supcol.lt.0.) then
2803	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
2804	psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
2805	! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
2806	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
2807	endif
2808	enddo
2809	enddo
2810	!
2811	!
2812	!----------------------------------------------------------------
2813	! check mass conservation of generation terms and feedback to the
2814	! large scale
2815	!
2816	do k = kts, kte
2817	do i = its, ite
2818	if(t(i,k).le.t0c) then
2819	!
2820	! cloud water
2821	!
2822	value = max(qmin,qci(i,k,1))
2823	source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
2824	if (source.gt.value) then
2825	factor = value/source
2826	praut(i,k) = praut(i,k)*factor
2827	pracw(i,k) = pracw(i,k)*factor
2828	psacw(i,k) = psacw(i,k)*factor
2829	endif
2830	!
2831	! cloud ice
2832	!
2833	value = max(qmin,qci(i,k,2))
2834	source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
2835	if (source.gt.value) then
2836	factor = value/source
2837	psaut(i,k) = psaut(i,k)*factor
2838	psaci(i,k) = psaci(i,k)*factor
2839	pigen(i,k) = pigen(i,k)*factor
2840	pidep(i,k) = pidep(i,k)*factor
2841	endif
2842	!
2843	! rain
2844	!
2845	!
2846	value = max(qmin,qrs(i,k,1))
2847	source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
2848	if (source.gt.value) then
2849	factor = value/source
2850	praut(i,k) = praut(i,k)*factor
2851	pracw(i,k) = pracw(i,k)*factor
2852	prevp(i,k) = prevp(i,k)*factor
2853	endif
2854	!
2855	! snow
2856	!
2857	value = max(qmin,qrs(i,k,2))
2858	source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
2859	if (source.gt.value) then
2860	factor = value/source
2861	psdep(i,k) = psdep(i,k)*factor
2862	psaut(i,k) = psaut(i,k)*factor
2863	psaci(i,k) = psaci(i,k)*factor
2864	psacw(i,k) = psacw(i,k)*factor
2865	endif
2866	!
2867	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
2868	! update
2869	q(i,k) = q(i,k)+work2(i,k)*dtcld
2870	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
2871	+psacw(i,k))*dtcld,0.)
2872	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
2873	+prevp(i,k))*dtcld,0.)
2874	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
2875	-pigen(i,k)-pidep(i,k))*dtcld,0.)
2876	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
2877	+psaci(i,k)+psacw(i,k))*dtcld,0.)
2878	xlf = xls-xl(i,k)
2879	xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
2880	-xl(i,k)prevp(i,k)-xlfpsacw(i,k)
2881	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2882	else
2883	!
2884	! cloud water
2885	!
2886	value = max(qmin,qci(i,k,1))
2887	source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
2888	if (source.gt.value) then
2889	factor = value/source
2890	praut(i,k) = praut(i,k)*factor
2891	pracw(i,k) = pracw(i,k)*factor
2892	psacw(i,k) = psacw(i,k)*factor
2893	endif
2894	!
2895	! rain
2896	!
2897	value = max(qmin,qrs(i,k,1))
2898	source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
2899	if (source.gt.value) then
2900	factor = value/source
2901	praut(i,k) = praut(i,k)*factor
2902	pracw(i,k) = pracw(i,k)*factor
2903	prevp(i,k) = prevp(i,k)*factor
2904	psacw(i,k) = psacw(i,k)*factor
2905	endif
2906	!
2907	! snow
2908	!
2909	value = max(qcrmin,qrs(i,k,2))
2910	source=(-psevp(i,k))*dtcld
2911	if (source.gt.value) then
2912	factor = value/source
2913	psevp(i,k) = psevp(i,k)*factor
2914	endif
2915	work2(i,k)=-(prevp(i,k)+psevp(i,k))
2916	! update
2917	q(i,k) = q(i,k)+work2(i,k)*dtcld
2918	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
2919	+psacw(i,k))*dtcld,0.)
2920	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
2921	+prevp(i,k) +psacw(i,k))*dtcld,0.)
2922	qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
2923	xlf = xls-xl(i,k)
2924	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
2925	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2926	endif
2927	enddo
2928	enddo
2929	!
2930	! Inline expansion for fpvs
2931	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2932	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2933	hsub = xls
2934	hvap = xlv0
2935	cvap = cpv
2936	ttp=t0c+0.01
2937	dldt=cvap-cliq
2938	xa=-dldt/rv
2939	xb=xa+hvap/(rv*ttp)
2940	dldti=cvap-cice
2941	xai=-dldti/rv
2942	xbi=xai+hsub/(rv*ttp)
2943	do k = kts, kte
2944	do i = its, ite
2945	tr=ttp/t(i,k)
2946	logtr = log(tr)
2947	qs(i,k,1)=psatexp(logtr(xa)+xb*(1.-tr))
2948	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2949	qs(i,k,1) = max(qs(i,k,1),qmin)
2950	enddo
2951	enddo
2952	!
2953	!----------------------------------------------------------------
2954	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
2955	! if there exists additional water vapor condensated/if
2956	! evaporation of cloud water is not enough to remove subsaturation
2957	!
2958	do k = kts, kte
2959	do i = its, ite
2960	! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
2961	work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
2962	(xl(i,k))/(rv(cpm(i,k)))*(qs(i,k,1)) &
2963	/((t(i,k))*(t(i,k)))))
2964	work2(i,k) = qci(i,k,1)+work1(i,k,1)
2965	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
2966	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
2967	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
2968	q(i,k) = q(i,k)-pcond(i,k)*dtcld
2969	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
2970	t(i,k) = t(i,k)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
2971	enddo
2972	enddo
2973	!
2974	!
2975	!----------------------------------------------------------------
2976	! padding for small values
2977	!
2978	do k = kts, kte
2979	do i = its, ite
2980	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
2981	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
2982	enddo
2983	enddo
2984	enddo ! big loops
2985	END SUBROUTINE wsm52d
2986
2987	#endif
2988
2989	! ...................................................................
2990	REAL FUNCTION rgmma(x)
2991	!-------------------------------------------------------------------
2992	IMPLICIT NONE
2993	!-------------------------------------------------------------------
2994	! rgmma function: use infinite product form
2995	REAL :: euler
2996	PARAMETER (euler=0.577215664901532)
2997	REAL :: x, y
2998	INTEGER :: i
2999	if(x.eq.1.)then
3000	rgmma=0.
3001	else
3002	rgmma=xexp(eulerx)
3003	do i=1,10000
3004	y=float(i)
3005	rgmma=rgmma(1.000+x/y)exp(-x/y)
3006	enddo
3007	rgmma=1./rgmma
3008	endif
3009	END FUNCTION rgmma
3010	!
3011	!--------------------------------------------------------------------------
3012	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
3013	!--------------------------------------------------------------------------
3014	IMPLICIT NONE
3015	!--------------------------------------------------------------------------
3016	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
3017	xai,xbi,ttp,tr
3018	INTEGER ice
3019	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3020	ttp=t0c+0.01
3021	dldt=cvap-cliq
3022	xa=-dldt/rv
3023	xb=xa+hvap/(rv*ttp)
3024	dldti=cvap-cice
3025	xai=-dldti/rv
3026	xbi=xai+hsub/(rv*ttp)
3027	tr=ttp/t
3028	if(t.lt.ttp.and.ice.eq.1) then
3029	fpvs=psatexp(log(tr)(xai))exp(xbi(1.-tr))
3030	else
3031	fpvs=psatexp(log(tr)(xa))exp(xb(1.-tr))
3032	endif
3033	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3034	END FUNCTION fpvs
3035	!-------------------------------------------------------------------
3036	SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
3037	!-------------------------------------------------------------------
3038	IMPLICIT NONE
3039	!-------------------------------------------------------------------
3040	!.... constants which may not be tunable
3041	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
3042	LOGICAL, INTENT(IN) :: allowed_to_read
3043	REAL :: pi
3044	!
3045	pi = 4.*atan(1.)
3046	xlv1 = cl-cpv
3047	!
3048	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
3049	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
3050	!
3051	bvtr1 = 1.+bvtr
3052	bvtr2 = 2.5+.5*bvtr
3053	bvtr3 = 3.+bvtr
3054	bvtr4 = 4.+bvtr
3055	g1pbr = rgmma(bvtr1)
3056	g3pbr = rgmma(bvtr3)
3057	g4pbr = rgmma(bvtr4) ! 17.837825
3058	g5pbro2 = rgmma(bvtr2) ! 1.8273
3059	pvtr = avtr*g4pbr/6.
3060	eacrr = 1.0
3061	pacrr = pin0ravtrg3pbr.25*eacrr
3062	precr1 = 2.pin0r*.78
3063	precr2 = 2.pin0r.31avtr*.5g5pbro2
3064	xmmax = (dimax/dicon)**2
3065	roqimax = 2.08e22dimax*8
3066	!
3067	bvts1 = 1.+bvts
3068	bvts2 = 2.5+.5*bvts
3069	bvts3 = 3.+bvts
3070	bvts4 = 4.+bvts
3071	g1pbs = rgmma(bvts1) !.8875
3072	g3pbs = rgmma(bvts3)
3073	g4pbs = rgmma(bvts4) ! 12.0786
3074	g5pbso2 = rgmma(bvts2)
3075	pvts = avts*g4pbs/6.
3076	pacrs = pin0savtsg3pbs.25
3077	precs1 = 4.n0s.65
3078	precs2 = 4.n0s.44avts.5g5pbso2
3079	pidn0r = pidenrn0r
3080	pidn0s = pidensn0s
3081	pacrc = pin0savtsg3pbs.25*eacrc
3082	!
3083	rslopermax = 1./lamdarmax
3084	rslopesmax = 1./lamdasmax
3085	rsloperbmax = rslopermax ** bvtr
3086	rslopesbmax = rslopesmax ** bvts
3087	rsloper2max = rslopermax * rslopermax
3088	rslopes2max = rslopesmax * rslopesmax
3089	rsloper3max = rsloper2max * rslopermax
3090	rslopes3max = rslopes2max * rslopesmax
3091	!
3092	END SUBROUTINE wsm5init
3093	END MODULE module_mp_wsm5

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

Download in other formats:

Original Format