Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm5.F @ 2295

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 68.5 KB

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

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

Download in other formats:

Original Format