Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm6.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: 92.6 KB

Line
1	#if ( RWORDSIZE == 4 )
2	# define VREC vsrec
3	# define VSQRT vssqrt
4	#else
5	# define VREC vrec
6	# define VSQRT vsqrt
7	#endif
8
9	MODULE module_mp_wsm6
10	!
11	!
12	REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
13	REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
14	REAL, PARAMETER, PRIVATE :: n0g = 4.e6 ! intercept parameter graupel
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 :: avtg = 330. ! a constant for terminal velocity of graupel
24	REAL, PARAMETER, PRIVATE :: bvtg = 0.8 ! a constant for terminal velocity of graupel
25	REAL, PARAMETER, PRIVATE :: deng = 500. ! density of graupel
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, PARAMETER, PRIVATE :: dens = 100.0 ! Density of snow
39	REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! threshold amount for aggretion to occur
40	REAL, SAVE :: &
41	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
42	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
43	bvtr6,g6pbr, &
44	precr1,precr2,roqimax,bvts1, &
45	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
46	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
47	pidn0s,xlv1,pacrc,pi, &
48	bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
49	g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
50	precg1,precg2,pidn0g, &
51	rslopermax,rslopesmax,rslopegmax, &
52	rsloperbmax,rslopesbmax,rslopegbmax, &
53	rsloper2max,rslopes2max,rslopeg2max, &
54	rsloper3max,rslopes3max,rslopeg3max
55	CONTAINS
56	!===================================================================
57	!
58	SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg &
59	,den, pii, p, delz &
60	,delt,g, cpd, cpv, rd, rv, t0c &
61	,ep1, ep2, qmin &
62	,XLS, XLV0, XLF0, den0, denr &
63	,cliq,cice,psat &
64	,rain, rainncv &
65	,snow, snowncv &
66	,sr &
67	,graupel, graupelncv &
68	,ids,ide, jds,jde, kds,kde &
69	,ims,ime, jms,jme, kms,kme &
70	,its,ite, jts,jte, kts,kte &
71	)
72	!-------------------------------------------------------------------
73	IMPLICIT NONE
74	!-------------------------------------------------------------------
75	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
76	ims,ime, jms,jme, kms,kme , &
77	its,ite, jts,jte, kts,kte
78	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
79	INTENT(INOUT) :: &
80	th, &
81	q, &
82	qc, &
83	qi, &
84	qr, &
85	qs, &
86	qg
87	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
88	INTENT(IN ) :: &
89	den, &
90	pii, &
91	p, &
92	delz
93	REAL, INTENT(IN ) :: delt, &
94	g, &
95	rd, &
96	rv, &
97	t0c, &
98	den0, &
99	cpd, &
100	cpv, &
101	ep1, &
102	ep2, &
103	qmin, &
104	XLS, &
105	XLV0, &
106	XLF0, &
107	cliq, &
108	cice, &
109	psat, &
110	denr
111	REAL, DIMENSION( ims:ime , jms:jme ), &
112	INTENT(INOUT) :: rain, &
113	rainncv, &
114	sr
115	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
116	INTENT(INOUT) :: snow, &
117	snowncv
118	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
119	INTENT(INOUT) :: graupel, &
120	graupelncv
121	! LOCAL VAR
122	REAL, DIMENSION( its:ite , kts:kte ) :: t
123	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
124	REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
125	INTEGER :: i,j,k
126	!-------------------------------------------------------------------
127	DO j=jts,jte
128	DO k=kts,kte
129	DO i=its,ite
130	t(i,k)=th(i,k,j)*pii(i,k,j)
131	qci(i,k,1) = qc(i,k,j)
132	qci(i,k,2) = qi(i,k,j)
133	qrs(i,k,1) = qr(i,k,j)
134	qrs(i,k,2) = qs(i,k,j)
135	qrs(i,k,3) = qg(i,k,j)
136	ENDDO
137	ENDDO
138	! Sending array starting locations of optional variables may cause
139	! troubles, so we explicitly change the call.
140	CALL wsm62D(t, q(ims,kms,j), qci, qrs &
141	,den(ims,kms,j) &
142	,p(ims,kms,j), delz(ims,kms,j) &
143	,delt,g, cpd, cpv, rd, rv, t0c &
144	,ep1, ep2, qmin &
145	,XLS, XLV0, XLF0, den0, denr &
146	,cliq,cice,psat &
147	,j &
148	,rain(ims,j),rainncv(ims,j) &
149	,sr(ims,j) &
150	,ids,ide, jds,jde, kds,kde &
151	,ims,ime, jms,jme, kms,kme &
152	,its,ite, jts,jte, kts,kte &
153	,snow,snowncv &
154	,graupel,graupelncv &
155	)
156	DO K=kts,kte
157	DO I=its,ite
158	th(i,k,j)=t(i,k)/pii(i,k,j)
159	qc(i,k,j) = qci(i,k,1)
160	qi(i,k,j) = qci(i,k,2)
161	qr(i,k,j) = qrs(i,k,1)
162	qs(i,k,j) = qrs(i,k,2)
163	qg(i,k,j) = qrs(i,k,3)
164	ENDDO
165	ENDDO
166	ENDDO
167	END SUBROUTINE wsm6
168	!===================================================================
169	!
170	SUBROUTINE wsm62D(t, q &
171	,qci, qrs, den, p, delz &
172	,delt,g, cpd, cpv, rd, rv, t0c &
173	,ep1, ep2, qmin &
174	,XLS, XLV0, XLF0, den0, denr &
175	,cliq,cice,psat &
176	,lat &
177	,rain,rainncv &
178	,sr &
179	,ids,ide, jds,jde, kds,kde &
180	,ims,ime, jms,jme, kms,kme &
181	,its,ite, jts,jte, kts,kte &
182	,snow,snowncv &
183	,graupel,graupelncv &
184	)
185	!-------------------------------------------------------------------
186	IMPLICIT NONE
187	!-------------------------------------------------------------------
188	!
189	! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the
190	! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
191	! number concentration is a function of temperature, and seperate assumption
192	! is developed, in which ice crystal number concentration is a function
193	! of ice amount. A theoretical background of the ice-microphysics and related
194	! processes in the WSMMPs are described in Hong et al. (2004).
195	! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
196	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
197	!
198	! WSM6 cloud scheme
199	!
200	! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
201	! Summer 2003
202	!
203	! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
204	! Summer 2004
205	!
206	! History : semi-lagrangian scheme sedimentation(JH), and clean up
207	! Hong, August 2009
208	!
209	! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
210	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
211	! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
212	! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
213	! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
214	! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
215	! Juang and Hong (JH, 2010) Mon. Wea. Rev.
216	!
217	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
218	ims,ime, jms,jme, kms,kme , &
219	its,ite, jts,jte, kts,kte, &
220	lat
221	REAL, DIMENSION( its:ite , kts:kte ), &
222	INTENT(INOUT) :: &
223	t
224	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
225	INTENT(INOUT) :: &
226	qci
227	REAL, DIMENSION( its:ite , kts:kte, 3 ), &
228	INTENT(INOUT) :: &
229	qrs
230	REAL, DIMENSION( ims:ime , kms:kme ), &
231	INTENT(INOUT) :: &
232	q
233	REAL, DIMENSION( ims:ime , kms:kme ), &
234	INTENT(IN ) :: &
235	den, &
236	p, &
237	delz
238	REAL, INTENT(IN ) :: delt, &
239	g, &
240	cpd, &
241	cpv, &
242	t0c, &
243	den0, &
244	rd, &
245	rv, &
246	ep1, &
247	ep2, &
248	qmin, &
249	XLS, &
250	XLV0, &
251	XLF0, &
252	cliq, &
253	cice, &
254	psat, &
255	denr
256	REAL, DIMENSION( ims:ime ), &
257	INTENT(INOUT) :: rain, &
258	rainncv, &
259	sr
260	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
261	INTENT(INOUT) :: snow, &
262	snowncv
263	REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
264	INTENT(INOUT) :: graupel, &
265	graupelncv
266	! LOCAL VAR
267	REAL, DIMENSION( its:ite , kts:kte , 3) :: &
268	rh, &
269	qs, &
270	rslope, &
271	rslope2, &
272	rslope3, &
273	rslopeb, &
274	qrs_tmp, &
275	falk, &
276	fall, &
277	work1
278	REAL, DIMENSION( its:ite , kts:kte ) :: &
279	fallc, &
280	falkc, &
281	work1c, &
282	work2c, &
283	workr, &
284	worka
285	REAL, DIMENSION( its:ite , kts:kte ) :: &
286	den_tmp, &
287	delz_tmp
288	REAL, DIMENSION( its:ite , kts:kte ) :: &
289	pigen, &
290	pidep, &
291	pcond, &
292	prevp, &
293	psevp, &
294	pgevp, &
295	psdep, &
296	pgdep, &
297	praut, &
298	psaut, &
299	pgaut, &
300	piacr, &
301	pracw, &
302	praci, &
303	pracs, &
304	psacw, &
305	psaci, &
306	psacr, &
307	pgacw, &
308	pgaci, &
309	pgacr, &
310	pgacs, &
311	paacw, &
312	psmlt, &
313	pgmlt, &
314	pseml, &
315	pgeml
316	REAL, DIMENSION( its:ite , kts:kte ) :: &
317	qsum, &
318	xl, &
319	cpm, &
320	work2, &
321	denfac, &
322	xni, &
323	denqrs1, &
324	denqrs2, &
325	denqrs3, &
326	denqci, &
327	n0sfac
328	REAL, DIMENSION( its:ite ) :: delqrs1, &
329	delqrs2, &
330	delqrs3, &
331	delqi
332	REAL, DIMENSION( its:ite ) :: tstepsnow, &
333	tstepgraup
334	INTEGER, DIMENSION( its:ite ) :: mstep, &
335	numdt
336	LOGICAL, DIMENSION( its:ite ) :: flgcld
337	REAL :: &
338	cpmcal, xlcal, diffus, &
339	viscos, xka, venfac, conden, diffac, &
340	x, y, z, a, b, c, d, e, &
341	qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
342	coeres, supsat, dtcld, xmi, eacrs, satdt, &
343	qimax, diameter, xni0, roqi0, &
344	fallsum, fallsum_qsi, fallsum_qg, &
345	vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
346	xlwork2, factor, source, value, &
347	xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
348	REAL :: vt2ave
349	REAL :: holdc, holdci
350	INTEGER :: i, j, k, mstepmax, &
351	iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
352	! Temporaries used for inlining fpvs function
353	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
354	! variables for optimization
355	REAL, DIMENSION( its:ite ) :: tvec1
356	REAL :: temp
357	!
358	!=================================================================
359	! compute internal functions
360	!
361	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
362	xlcal(x) = xlv0-xlv1*(x-t0c)
363	!----------------------------------------------------------------
364	! diffus: diffusion coefficient of the water vapor
365	! viscos: kinematic viscosity(m2s-1)
366	! Optimizatin : A*B => exp(log(A)(B))
367	!
368	diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
369	viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
370	xka(x,y) = 1.414e3viscos(x,y)y
371	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
372	venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
373	/sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
374	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
375	!
376	!
377	idim = ite-its+1
378	kdim = kte-kts+1
379	!
380	!----------------------------------------------------------------
381	! paddint 0 for negative values generated by dynamics
382	!
383	do k = kts, kte
384	do i = its, ite
385	qci(i,k,1) = max(qci(i,k,1),0.0)
386	qrs(i,k,1) = max(qrs(i,k,1),0.0)
387	qci(i,k,2) = max(qci(i,k,2),0.0)
388	qrs(i,k,2) = max(qrs(i,k,2),0.0)
389	qrs(i,k,3) = max(qrs(i,k,3),0.0)
390	enddo
391	enddo
392	!
393	!----------------------------------------------------------------
394	! latent heat for phase changes and heat capacity. neglect the
395	! changes during microphysical process calculation
396	! emanuel(1994)
397	!
398	do k = kts, kte
399	do i = its, ite
400	cpm(i,k) = cpmcal(q(i,k))
401	xl(i,k) = xlcal(t(i,k))
402	enddo
403	enddo
404	do k = kts, kte
405	do i = its, ite
406	delz_tmp(i,k) = delz(i,k)
407	den_tmp(i,k) = den(i,k)
408	enddo
409	enddo
410	!
411	!----------------------------------------------------------------
412	! initialize the surface rain, snow, graupel
413	!
414	do i = its, ite
415	rainncv(i) = 0.
416	if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i,lat) = 0.
417	if(PRESENT (graupelncv) .AND. PRESENT (graupel)) graupelncv(i,lat) = 0.
418	sr(i) = 0.
419	! new local array to catch step snow and graupel
420	tstepsnow(i) = 0.
421	tstepgraup(i) = 0.
422	enddo
423	!
424	!----------------------------------------------------------------
425	! compute the minor time steps.
426	!
427	loops = max(nint(delt/dtcldcr),1)
428	dtcld = delt/loops
429	if(delt.le.dtcldcr) dtcld = delt
430	!
431	do loop = 1,loops
432	!
433	!----------------------------------------------------------------
434	! initialize the large scale variables
435	!
436	do i = its, ite
437	mstep(i) = 1
438	flgcld(i) = .true.
439	enddo
440	!
441	! do k = kts, kte
442	! do i = its, ite
443	! denfac(i,k) = sqrt(den0/den(i,k))
444	! enddo
445	! enddo
446	do k = kts, kte
447	CALL VREC( tvec1(its), den(its,k), ite-its+1)
448	do i = its, ite
449	tvec1(i) = tvec1(i)*den0
450	enddo
451	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
452	enddo
453	!
454	! Inline expansion for fpvs
455	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
456	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
457	hsub = xls
458	hvap = xlv0
459	cvap = cpv
460	ttp=t0c+0.01
461	dldt=cvap-cliq
462	xa=-dldt/rv
463	xb=xa+hvap/(rv*ttp)
464	dldti=cvap-cice
465	xai=-dldti/rv
466	xbi=xai+hsub/(rv*ttp)
467	do k = kts, kte
468	do i = its, ite
469	tr=ttp/t(i,k)
470	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
471	qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
472	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
473	qs(i,k,1) = max(qs(i,k,1),qmin)
474	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
475	tr=ttp/t(i,k)
476	if(t(i,k).lt.ttp) then
477	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
478	else
479	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
480	endif
481	qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
482	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
483	qs(i,k,2) = max(qs(i,k,2),qmin)
484	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
485	enddo
486	enddo
487	!
488	!----------------------------------------------------------------
489	! initialize the variables for microphysical physics
490	!
491	!
492	do k = kts, kte
493	do i = its, ite
494	prevp(i,k) = 0.
495	psdep(i,k) = 0.
496	pgdep(i,k) = 0.
497	praut(i,k) = 0.
498	psaut(i,k) = 0.
499	pgaut(i,k) = 0.
500	pracw(i,k) = 0.
501	praci(i,k) = 0.
502	piacr(i,k) = 0.
503	psaci(i,k) = 0.
504	psacw(i,k) = 0.
505	pracs(i,k) = 0.
506	psacr(i,k) = 0.
507	pgacw(i,k) = 0.
508	paacw(i,k) = 0.
509	pgaci(i,k) = 0.
510	pgacr(i,k) = 0.
511	pgacs(i,k) = 0.
512	pigen(i,k) = 0.
513	pidep(i,k) = 0.
514	pcond(i,k) = 0.
515	psmlt(i,k) = 0.
516	pgmlt(i,k) = 0.
517	pseml(i,k) = 0.
518	pgeml(i,k) = 0.
519	psevp(i,k) = 0.
520	pgevp(i,k) = 0.
521	falk(i,k,1) = 0.
522	falk(i,k,2) = 0.
523	falk(i,k,3) = 0.
524	fall(i,k,1) = 0.
525	fall(i,k,2) = 0.
526	fall(i,k,3) = 0.
527	fallc(i,k) = 0.
528	falkc(i,k) = 0.
529	xni(i,k) = 1.e3
530	enddo
531	enddo
532	!-------------------------------------------------------------
533	! Ni: ice crystal number concentraiton [HDC 5c]
534	!-------------------------------------------------------------
535	do k = kts, kte
536	do i = its, ite
537	temp = (den(i,k)*max(qci(i,k,2),qmin))
538	temp = sqrt(sqrt(temptemptemp))
539	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
540	enddo
541	enddo
542	!
543	!----------------------------------------------------------------
544	! compute the fallout term:
545	! first, vertical terminal velosity for minor loops
546	!----------------------------------------------------------------
547	do k = kts, kte
548	do i = its, ite
549	qrs_tmp(i,k,1) = qrs(i,k,1)
550	qrs_tmp(i,k,2) = qrs(i,k,2)
551	qrs_tmp(i,k,3) = qrs(i,k,3)
552	enddo
553	enddo
554	call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
555	work1,its,ite,kts,kte)
556	!
557	do k = kte, kts, -1
558	do i = its, ite
559	workr(i,k) = work1(i,k,1)
560	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
561	IF ( qsum(i,k) .gt. 1.e-15 ) THEN
562	worka(i,k) = (work1(i,k,2)qrs(i,k,2) + work1(i,k,3)qrs(i,k,3)) &
563	/qsum(i,k)
564	ELSE
565	worka(i,k) = 0.
566	ENDIF
567	denqrs1(i,k) = den(i,k)*qrs(i,k,1)
568	denqrs2(i,k) = den(i,k)*qrs(i,k,2)
569	denqrs3(i,k) = den(i,k)*qrs(i,k,3)
570	if(qrs(i,k,1).le.0.0) workr(i,k) = 0.0
571	enddo
572	enddo
573	call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
574	delqrs1,dtcld,1,1)
575	call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
576	denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
577	do k = kts, kte
578	do i = its, ite
579	qrs(i,k,1) = max(denqrs1(i,k)/den(i,k),0.)
580	qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
581	qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
582	fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
583	fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
584	fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
585	enddo
586	enddo
587	do i = its, ite
588	fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
589	fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
590	fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
591	enddo
592	do k = kts, kte
593	do i = its, ite
594	qrs_tmp(i,k,1) = qrs(i,k,1)
595	qrs_tmp(i,k,2) = qrs(i,k,2)
596	qrs_tmp(i,k,3) = qrs(i,k,3)
597	enddo
598	enddo
599	call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
600	work1,its,ite,kts,kte)
601	!
602	do k = kte, kts, -1
603	do i = its, ite
604	supcol = t0c-t(i,k)
605	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
606	if(t(i,k).gt.t0c) then
607	!---------------------------------------------------------------
608	! psmlt: melting of snow [HL A33] [RH83 A25]
609	! (T>T0: S->R)
610	!---------------------------------------------------------------
611	xlf = xlf0
612	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
613	if(qrs(i,k,2).gt.0.) then
614	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
615	psmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
616	n0sfac(i,k)(precs1*rslope2(i,k,2) &
617	+precs2work2(i,k)coeres)
618	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
619	-qrs(i,k,2)/mstep(i)),0.)
620	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
621	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
622	t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
623	endif
624	!---------------------------------------------------------------
625	! pgmlt: melting of graupel [HL A23] [LFO 47]
626	! (T>T0: G->R)
627	!---------------------------------------------------------------
628	if(qrs(i,k,3).gt.0.) then
629	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
630	pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
631	(t0c-t(i,k))(precg1*rslope2(i,k,3) &
632	+precg2work2(i,k)coeres)
633	pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
634	-qrs(i,k,3)/mstep(i)),0.)
635	qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
636	qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
637	t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
638	endif
639	endif
640	enddo
641	enddo
642	!---------------------------------------------------------------
643	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
644	!---------------------------------------------------------------
645	do k = kte, kts, -1
646	do i = its, ite
647	if(qci(i,k,2).le.0.) then
648	work1c(i,k) = 0.
649	else
650	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
651	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
652	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
653	endif
654	enddo
655	enddo
656	!
657	! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
658	!
659	do k = kte, kts, -1
660	do i = its, ite
661	denqci(i,k) = den(i,k)*qci(i,k,2)
662	enddo
663	enddo
664	call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
665	delqi,dtcld,1,0)
666	do k = kts, kte
667	do i = its, ite
668	qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
669	enddo
670	enddo
671	do i = its, ite
672	fallc(i,1) = delqi(i)/delz(i,1)/dtcld
673	enddo
674	!
675	!----------------------------------------------------------------
676	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
677	!
678	do i = its, ite
679	fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
680	fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
681	fallsum_qg = fall(i,kts,3)
682	if(fallsum.gt.0.) then
683	rainncv(i) = fallsumdelz(i,kts)/denrdtcld*1000. + rainncv(i)
684	rain(i) = fallsumdelz(i,kts)/denrdtcld*1000. + rain(i)
685	endif
686	if(fallsum_qsi.gt.0.) then
687	tstepsnow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. &
688	+tstepsnow(i)
689	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
690	snowncv(i,lat) = fallsum_qsidelz(i,kts)/denrdtcld*1000. &
691	+snowncv(i,lat)
692	snow(i,lat) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i,lat)
693	ENDIF
694	endif
695	if(fallsum_qg.gt.0.) then
696	tstepgraup(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. &
697	+tstepgraup(i)
698	IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
699	graupelncv(i,lat) = fallsum_qgdelz(i,kts)/denrdtcld*1000. &
700	+ graupelncv(i,lat)
701	graupel(i,lat) = fallsum_qgdelz(i,kts)/denrdtcld*1000. + graupel(i,lat)
702	ENDIF
703	endif
704	! if(fallsum.gt.0.)sr(i)=(snowncv(i,lat) + graupelncv(i,lat))/(rainncv(i)+1.e-12)
705	if(fallsum.gt.0.)sr(i)=(tstepsnow(i) + tstepgraup(i))/(rainncv(i)+1.e-12)
706	enddo
707	!
708	!---------------------------------------------------------------
709	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
710	! (T>T0: I->C)
711	!---------------------------------------------------------------
712	do k = kts, kte
713	do i = its, ite
714	supcol = t0c-t(i,k)
715	xlf = xls-xl(i,k)
716	if(supcol.lt.0.) xlf = xlf0
717	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
718	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
719	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
720	qci(i,k,2) = 0.
721	endif
722	!---------------------------------------------------------------
723	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
724	! (T<-40C: C->I)
725	!---------------------------------------------------------------
726	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
727	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
728	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
729	qci(i,k,1) = 0.
730	endif
731	!---------------------------------------------------------------
732	! pihtf: heterogeneous freezing of cloud water [HL A44]
733	! (T0>T>-40C: C->I)
734	!---------------------------------------------------------------
735	if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
736	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
737	! den(i,k)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
738	supcolt=min(supcol,50.)
739	pfrzdtc = min(pfrz1(exp(pfrz2supcolt)-1.) &
740	den(i,k)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
741	qci(i,k,2) = qci(i,k,2) + pfrzdtc
742	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
743	qci(i,k,1) = qci(i,k,1)-pfrzdtc
744	endif
745	!---------------------------------------------------------------
746	! pgfrz: freezing of rain water [HL A20] [LFO 45]
747	! (T<T0, R->G)
748	!---------------------------------------------------------------
749	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
750	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k) &
751	! (exp(pfrz2supcol)-1.)rslope3(i,k,1)*2 &
752	! rslope(i,k,1)dtcld,qrs(i,k,1))
753	temp = rslope3(i,k,1)
754	temp = temptemprslope(i,k,1)
755	supcolt=min(supcol,50.)
756	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k) &
757	(exp(pfrz2supcolt)-1.)tempdtcld, &
758	qrs(i,k,1))
759	qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
760	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
761	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
762	endif
763	enddo
764	enddo
765	!
766	!
767	!----------------------------------------------------------------
768	! update the slope parameters for microphysics computation
769	!
770	do k = kts, kte
771	do i = its, ite
772	qrs_tmp(i,k,1) = qrs(i,k,1)
773	qrs_tmp(i,k,2) = qrs(i,k,2)
774	qrs_tmp(i,k,3) = qrs(i,k,3)
775	enddo
776	enddo
777	call slope_wsm6(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
778	work1,its,ite,kts,kte)
779	!------------------------------------------------------------------
780	! work1: the thermodynamic term in the denominator associated with
781	! heat conduction and vapor diffusion
782	! (ry88, y93, h85)
783	! work2: parameter associated with the ventilation effects(y93)
784	!
785	do k = kts, kte
786	do i = its, ite
787	work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
788	work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
789	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
790	enddo
791	enddo
792	!
793	!===============================================================
794	!
795	! warm rain processes
796	!
797	! - follows the processes in RH83 and LFO except for autoconcersion
798	!
799	!===============================================================
800	!
801	do k = kts, kte
802	do i = its, ite
803	supsat = max(q(i,k),qmin)-qs(i,k,1)
804	satdt = supsat/dtcld
805	!---------------------------------------------------------------
806	! praut: auto conversion rate from cloud to rain [HDC 16]
807	! (C->R)
808	!---------------------------------------------------------------
809	if(qci(i,k,1).gt.qc0) then
810	praut(i,k) = qck1qci(i,k,1)*(7./3.)
811	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
812	endif
813	!---------------------------------------------------------------
814	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
815	! (C->R)
816	!---------------------------------------------------------------
817	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
818	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
819	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
820	endif
821	!---------------------------------------------------------------
822	! prevp: evaporation/condensation rate of rain [HDC 14]
823	! (V->R or R->V)
824	!---------------------------------------------------------------
825	if(qrs(i,k,1).gt.0.) then
826	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
827	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
828	+precr2work2(i,k)coeres)/work1(i,k,1)
829	if(prevp(i,k).lt.0.) then
830	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
831	prevp(i,k) = max(prevp(i,k),satdt/2)
832	else
833	prevp(i,k) = min(prevp(i,k),satdt/2)
834	endif
835	endif
836	enddo
837	enddo
838	!
839	!===============================================================
840	!
841	! cold rain processes
842	!
843	! - follows the revised ice microphysics processes in HDC
844	! - the processes same as in RH83 and RH84 and LFO behave
845	! following ice crystal hapits defined in HDC, inclduing
846	! intercept parameter for snow (n0s), ice crystal number
847	! concentration (ni), ice nuclei number concentration
848	! (n0i), ice diameter (d)
849	!
850	!===============================================================
851	!
852	do k = kts, kte
853	do i = its, ite
854	supcol = t0c-t(i,k)
855	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
856	supsat = max(q(i,k),qmin)-qs(i,k,2)
857	satdt = supsat/dtcld
858	ifsat = 0
859	!-------------------------------------------------------------
860	! Ni: ice crystal number concentraiton [HDC 5c]
861	!-------------------------------------------------------------
862	! xni(i,k) = min(max(5.38e7*(den(i,k) &
863	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
864	temp = (den(i,k)*max(qci(i,k,2),qmin))
865	temp = sqrt(sqrt(temptemptemp))
866	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
867	eacrs = exp(0.07*(-supcol))
868	!
869	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
870	diameter = min(dicon * sqrt(xmi),dimax)
871	vt2i = 1.49e4diameter*1.31
872	vt2r=pvtrrslopeb(i,k,1)denfac(i,k)
873	vt2s=pvtsrslopeb(i,k,2)denfac(i,k)
874	vt2g=pvtgrslopeb(i,k,3)denfac(i,k)
875	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
876	if(qsum(i,k) .gt. 1.e-15) then
877	vt2ave=(vt2sqrs(i,k,2)+vt2gqrs(i,k,3))/(qsum(i,k))
878	else
879	vt2ave=0.
880	endif
881	if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
882	if(qrs(i,k,1).gt.qcrmin) then
883	!-------------------------------------------------------------
884	! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
885	! (T<T0: I->R)
886	!-------------------------------------------------------------
887	acrfac = 2.rslope3(i,k,1)+2.diameter*rslope2(i,k,1) &
888	+diameter*2rslope(i,k,1)
889	praci(i,k) = piqci(i,k,2)n0rabs(vt2r-vt2i)acrfac/4.
890	praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
891	!-------------------------------------------------------------
892	! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
893	! (T<T0: R->S or R->G)
894	!-------------------------------------------------------------
895	piacr(i,k) = pi*2avtrn0rdenrxni(i,k)denfac(i,k) &
896	g6pbrrslope3(i,k,1)*rslope3(i,k,1) &
897	*rslopeb(i,k,1)/24./den(i,k)
898	piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
899	endif
900	!-------------------------------------------------------------
901	! psaci: Accretion of cloud ice by snow [HDC 10]
902	! (T<T0: I->S)
903	!-------------------------------------------------------------
904	if(qrs(i,k,2).gt.qcrmin) then
905	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
906	+diameter*2rslope(i,k,2)
907	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
908	abs(vt2ave-vt2i)acrfac/4.
909	psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
910	endif
911	!-------------------------------------------------------------
912	! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
913	! (T<T0: I->G)
914	!-------------------------------------------------------------
915	if(qrs(i,k,3).gt.qcrmin) then
916	egi = exp(0.07*(-supcol))
917	acrfac = 2.rslope3(i,k,3)+2.diameter*rslope2(i,k,3) &
918	+diameter*2rslope(i,k,3)
919	pgaci(i,k) = piegiqci(i,k,2)n0gabs(vt2ave-vt2i)*acrfac/4.
920	pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
921	endif
922	endif
923	!-------------------------------------------------------------
924	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
925	! (T<T0: C->S, and T>=T0: C->R)
926	!-------------------------------------------------------------
927	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
928	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2)*rslopeb(i,k,2) &
929	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
930	endif
931	!-------------------------------------------------------------
932	! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
933	! (T<T0: C->G, and T>=T0: C->R)
934	!-------------------------------------------------------------
935	if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
936	pgacw(i,k) = min(pacrgrslope3(i,k,3)rslopeb(i,k,3) &
937	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
938	endif
939	!-------------------------------------------------------------
940	! paacw: Accretion of cloud water by averaged snow/graupel
941	! (T<T0: C->G or S, and T>=T0: C->R)
942	!-------------------------------------------------------------
943	if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,3).gt.qcrmin) then
944	paacw(i,k) = (qrs(i,k,2)psacw(i,k)+qrs(i,k,3)pgacw(i,k)) &
945	/(qsum(i,k))
946	endif
947	!-------------------------------------------------------------
948	! pracs: Accretion of snow by rain [HL A11] [LFO 27]
949	! (T<T0: S->G)
950	!-------------------------------------------------------------
951	if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
952	if(supcol.gt.0) then
953	acrfac = 5.rslope3(i,k,2)rslope3(i,k,2)*rslope(i,k,1) &
954	+2.rslope3(i,k,2)rslope2(i,k,2)*rslope2(i,k,1) &
955	+.5rslope2(i,k,2)rslope2(i,k,2)*rslope3(i,k,1)
956	pracs(i,k) = pi*2n0rn0sn0sfac(i,k)*abs(vt2r-vt2ave) &
957	(dens/den(i,k))acrfac
958	pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
959	endif
960	!-------------------------------------------------------------
961	! psacr: Accretion of rain by snow [HL A10] [LFO 28]
962	! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
963	!-------------------------------------------------------------
964	acrfac = 5.rslope3(i,k,1)rslope3(i,k,1)*rslope(i,k,2) &
965	+2.rslope3(i,k,1)rslope2(i,k,1)*rslope2(i,k,2) &
966	+.5rslope2(i,k,1)rslope2(i,k,1)*rslope3(i,k,2)
967	psacr(i,k) = pi*2n0rn0sn0sfac(i,k)*abs(vt2ave-vt2r) &
968	(denr/den(i,k))acrfac
969	psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
970	endif
971	!-------------------------------------------------------------
972	! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
973	! (T<T0: R->G) (T>=T0: enhance melting of graupel)
974	!-------------------------------------------------------------
975	if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
976	acrfac = 5.rslope3(i,k,1)rslope3(i,k,1)*rslope(i,k,3) &
977	+2.rslope3(i,k,1)rslope2(i,k,1)*rslope2(i,k,3) &
978	+.5rslope2(i,k,1)rslope2(i,k,1)*rslope3(i,k,3)
979	pgacr(i,k) = pi*2n0rn0gabs(vt2ave-vt2r)*(denr/den(i,k)) &
980	*acrfac
981	pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
982	endif
983	!
984	!-------------------------------------------------------------
985	! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
986	! (S->G): This process is eliminated in V3.0 with the
987	! new combined snow/graupel fall speeds
988	!-------------------------------------------------------------
989	if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
990	pgacs(i,k) = 0.
991	endif
992	if(supcol.le.0) then
993	xlf = xlf0
994	!-------------------------------------------------------------
995	! pseml: Enhanced melting of snow by accretion of water [HL A34]
996	! (T>=T0: S->R)
997	!-------------------------------------------------------------
998	if(qrs(i,k,2).gt.0.) &
999	pseml(i,k) = min(max(cliqsupcol(paacw(i,k)+psacr(i,k)) &
1000	/xlf,-qrs(i,k,2)/dtcld),0.)
1001	!-------------------------------------------------------------
1002	! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1003	! (T>=T0: G->R)
1004	!-------------------------------------------------------------
1005	if(qrs(i,k,3).gt.0.) &
1006	pgeml(i,k) = min(max(cliqsupcol(paacw(i,k)+pgacr(i,k)) &
1007	/xlf,-qrs(i,k,3)/dtcld),0.)
1008	endif
1009	if(supcol.gt.0) then
1010	!-------------------------------------------------------------
1011	! pidep: Deposition/Sublimation rate of ice [HDC 9]
1012	! (T<T0: V->I or I->V)
1013	!-------------------------------------------------------------
1014	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
1015	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1016	supice = satdt-prevp(i,k)
1017	if(pidep(i,k).lt.0.) then
1018	pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1019	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1020	else
1021	pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1022	endif
1023	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1024	endif
1025	!-------------------------------------------------------------
1026	! psdep: deposition/sublimation rate of snow [HDC 14]
1027	! (T<T0: V->S or S->V)
1028	!-------------------------------------------------------------
1029	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1030	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
1031	psdep(i,k) = (rh(i,k,2)-1.)n0sfac(i,k)(precs1*rslope2(i,k,2) &
1032	+ precs2work2(i,k)coeres)/work1(i,k,2)
1033	supice = satdt-prevp(i,k)-pidep(i,k)
1034	if(psdep(i,k).lt.0.) then
1035	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1036	psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1037	else
1038	psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1039	endif
1040	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1041	ifsat = 1
1042	endif
1043	!-------------------------------------------------------------
1044	! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1045	! (T<T0: V->G or G->V)
1046	!-------------------------------------------------------------
1047	if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
1048	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
1049	pgdep(i,k) = (rh(i,k,2)-1.)(precg1rslope2(i,k,3) &
1050	+precg2work2(i,k)coeres)/work1(i,k,2)
1051	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1052	if(pgdep(i,k).lt.0.) then
1053	pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
1054	pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1055	else
1056	pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1057	endif
1058	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1059	abs(satdt)) ifsat = 1
1060	endif
1061	!-------------------------------------------------------------
1062	! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1063	! (T<T0: V->I)
1064	!-------------------------------------------------------------
1065	if(supsat.gt.0.and.ifsat.ne.1) then
1066	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1067	xni0 = 1.e3exp(0.1supcol)
1068	roqi0 = 4.92e-11xni0*1.33
1069	pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
1070	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1071	endif
1072	!
1073	!-------------------------------------------------------------
1074	! psaut: conversion(aggregation) of ice to snow [HDC 12]
1075	! (T<T0: I->S)
1076	!-------------------------------------------------------------
1077	if(qci(i,k,2).gt.0.) then
1078	qimax = roqimax/den(i,k)
1079	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1080	endif
1081	!
1082	!-------------------------------------------------------------
1083	! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1084	! (T<T0: S->G)
1085	!-------------------------------------------------------------
1086	if(qrs(i,k,2).gt.0.) then
1087	alpha2 = 1.e-3exp(0.09(-supcol))
1088	pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
1089	endif
1090	endif
1091	!
1092	!-------------------------------------------------------------
1093	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1094	! (T>=T0: S->V)
1095	!-------------------------------------------------------------
1096	if(supcol.lt.0.) then
1097	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
1098	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
1099	psevp(i,k) = (rh(i,k,1)-1.)n0sfac(i,k)(precs1 &
1100	rslope2(i,k,2)+precs2work2(i,k) &
1101	*coeres)/work1(i,k,1)
1102	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1103	endif
1104	!-------------------------------------------------------------
1105	! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
1106	! (T>=T0: G->V)
1107	!-------------------------------------------------------------
1108	if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
1109	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
1110	pgevp(i,k) = (rh(i,k,1)-1.)(precg1rslope2(i,k,3) &
1111	+precg2work2(i,k)coeres)/work1(i,k,1)
1112	pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
1113	endif
1114	endif
1115	enddo
1116	enddo
1117	!
1118	!
1119	!----------------------------------------------------------------
1120	! check mass conservation of generation terms and feedback to the
1121	! large scale
1122	!
1123	do k = kts, kte
1124	do i = its, ite
1125	!
1126	delta2=0.
1127	delta3=0.
1128	if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
1129	if(qrs(i,k,1).lt.1.e-4) delta3=1.
1130	if(t(i,k).le.t0c) then
1131	!
1132	! cloud water
1133	!
1134	value = max(qmin,qci(i,k,1))
1135	source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1136	if (source.gt.value) then
1137	factor = value/source
1138	praut(i,k) = praut(i,k)*factor
1139	pracw(i,k) = pracw(i,k)*factor
1140	paacw(i,k) = paacw(i,k)*factor
1141	endif
1142	!
1143	! cloud ice
1144	!
1145	value = max(qmin,qci(i,k,2))
1146	source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
1147	+pgaci(i,k))*dtcld
1148	if (source.gt.value) then
1149	factor = value/source
1150	psaut(i,k) = psaut(i,k)*factor
1151	pigen(i,k) = pigen(i,k)*factor
1152	pidep(i,k) = pidep(i,k)*factor
1153	praci(i,k) = praci(i,k)*factor
1154	psaci(i,k) = psaci(i,k)*factor
1155	pgaci(i,k) = pgaci(i,k)*factor
1156	endif
1157	!
1158	! rain
1159	!
1160	value = max(qmin,qrs(i,k,1))
1161	source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k)+psacr(i,k) &
1162	+pgacr(i,k))*dtcld
1163	if (source.gt.value) then
1164	factor = value/source
1165	praut(i,k) = praut(i,k)*factor
1166	prevp(i,k) = prevp(i,k)*factor
1167	pracw(i,k) = pracw(i,k)*factor
1168	piacr(i,k) = piacr(i,k)*factor
1169	psacr(i,k) = psacr(i,k)*factor
1170	pgacr(i,k) = pgacr(i,k)*factor
1171	endif
1172	!
1173	! snow
1174	!
1175	value = max(qmin,qrs(i,k,2))
1176	source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k)+piacr(i,k) &
1177	delta3+praci(i,k)delta3-pracs(i,k)*(1.-delta2) &
1178	+psacr(i,k)delta2+psaci(i,k)-pgacs(i,k) )dtcld
1179	if (source.gt.value) then
1180	factor = value/source
1181	psdep(i,k) = psdep(i,k)*factor
1182	psaut(i,k) = psaut(i,k)*factor
1183	pgaut(i,k) = pgaut(i,k)*factor
1184	paacw(i,k) = paacw(i,k)*factor
1185	piacr(i,k) = piacr(i,k)*factor
1186	praci(i,k) = praci(i,k)*factor
1187	psaci(i,k) = psaci(i,k)*factor
1188	pracs(i,k) = pracs(i,k)*factor
1189	psacr(i,k) = psacr(i,k)*factor
1190	pgacs(i,k) = pgacs(i,k)*factor
1191	endif
1192	!
1193	! graupel
1194	!
1195	value = max(qmin,qrs(i,k,3))
1196	source = -(pgdep(i,k)+pgaut(i,k) &
1197	+piacr(i,k)(1.-delta3)+praci(i,k)(1.-delta3) &
1198	+psacr(i,k)(1.-delta2)+pracs(i,k)(1.-delta2) &
1199	+pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1200	if (source.gt.value) then
1201	factor = value/source
1202	pgdep(i,k) = pgdep(i,k)*factor
1203	pgaut(i,k) = pgaut(i,k)*factor
1204	piacr(i,k) = piacr(i,k)*factor
1205	praci(i,k) = praci(i,k)*factor
1206	psacr(i,k) = psacr(i,k)*factor
1207	pracs(i,k) = pracs(i,k)*factor
1208	paacw(i,k) = paacw(i,k)*factor
1209	pgaci(i,k) = pgaci(i,k)*factor
1210	pgacr(i,k) = pgacr(i,k)*factor
1211	pgacs(i,k) = pgacs(i,k)*factor
1212	endif
1213	!
1214	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
1215	! update
1216	q(i,k) = q(i,k)+work2(i,k)*dtcld
1217	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1218	+paacw(i,k)+paacw(i,k))*dtcld,0.)
1219	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1220	+prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1221	-psacr(i,k))*dtcld,0.)
1222	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
1223	+psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1224	*dtcld,0.)
1225	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1226	-pgaut(i,k)+piacr(i,k)*delta3 &
1227	+praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1228	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2) &
1229	*dtcld,0.)
1230	qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
1231	+piacr(i,k)*(1.-delta3) &
1232	+praci(i,k)(1.-delta3)+psacr(i,k)(1.-delta2) &
1233	+pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1234	+pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1235	xlf = xls-xl(i,k)
1236	xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1237	-xl(i,k)prevp(i,k)-xlf(piacr(i,k)+paacw(i,k) &
1238	+paacw(i,k)+pgacr(i,k)+psacr(i,k))
1239	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1240	else
1241	!
1242	! cloud water
1243	!
1244	value = max(qmin,qci(i,k,1))
1245	source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
1246	if (source.gt.value) then
1247	factor = value/source
1248	praut(i,k) = praut(i,k)*factor
1249	pracw(i,k) = pracw(i,k)*factor
1250	paacw(i,k) = paacw(i,k)*factor
1251	endif
1252	!
1253	! rain
1254	!
1255	value = max(qmin,qrs(i,k,1))
1256	source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k)-pracw(i,k) &
1257	-paacw(i,k)-prevp(i,k))*dtcld
1258	if (source.gt.value) then
1259	factor = value/source
1260	praut(i,k) = praut(i,k)*factor
1261	prevp(i,k) = prevp(i,k)*factor
1262	pracw(i,k) = pracw(i,k)*factor
1263	paacw(i,k) = paacw(i,k)*factor
1264	pseml(i,k) = pseml(i,k)*factor
1265	pgeml(i,k) = pgeml(i,k)*factor
1266	endif
1267	!
1268	! snow
1269	!
1270	value = max(qcrmin,qrs(i,k,2))
1271	source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1272	if (source.gt.value) then
1273	factor = value/source
1274	pgacs(i,k) = pgacs(i,k)*factor
1275	psevp(i,k) = psevp(i,k)*factor
1276	pseml(i,k) = pseml(i,k)*factor
1277	endif
1278	!
1279	! graupel
1280	!
1281	value = max(qcrmin,qrs(i,k,3))
1282	source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1283	if (source.gt.value) then
1284	factor = value/source
1285	pgacs(i,k) = pgacs(i,k)*factor
1286	pgevp(i,k) = pgevp(i,k)*factor
1287	pgeml(i,k) = pgeml(i,k)*factor
1288	endif
1289	work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1290	! update
1291	q(i,k) = q(i,k)+work2(i,k)*dtcld
1292	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1293	+paacw(i,k)+paacw(i,k))*dtcld,0.)
1294	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1295	+prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1296	-pgeml(i,k))*dtcld,0.)
1297	qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
1298	+pseml(i,k))*dtcld,0.)
1299	qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
1300	+pgeml(i,k))*dtcld,0.)
1301	xlf = xls-xl(i,k)
1302	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1303	-xlf*(pseml(i,k)+pgeml(i,k))
1304	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1305	endif
1306	enddo
1307	enddo
1308	!
1309	! Inline expansion for fpvs
1310	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1311	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1312	hsub = xls
1313	hvap = xlv0
1314	cvap = cpv
1315	ttp=t0c+0.01
1316	dldt=cvap-cliq
1317	xa=-dldt/rv
1318	xb=xa+hvap/(rv*ttp)
1319	dldti=cvap-cice
1320	xai=-dldti/rv
1321	xbi=xai+hsub/(rv*ttp)
1322	do k = kts, kte
1323	do i = its, ite
1324	tr=ttp/t(i,k)
1325	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
1326	qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
1327	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1328	qs(i,k,1) = max(qs(i,k,1),qmin)
1329	tr=ttp/t(i,k)
1330	if(t(i,k).lt.ttp) then
1331	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
1332	else
1333	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
1334	endif
1335	qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
1336	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1337	qs(i,k,2) = max(qs(i,k,2),qmin)
1338	enddo
1339	enddo
1340	!
1341	!----------------------------------------------------------------
1342	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1343	! if there exists additional water vapor condensated/if
1344	! evaporation of cloud water is not enough to remove subsaturation
1345	!
1346	do k = kts, kte
1347	do i = its, ite
1348	work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1349	work2(i,k) = qci(i,k,1)+work1(i,k,1)
1350	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1351	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1352	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1353	q(i,k) = q(i,k)-pcond(i,k)*dtcld
1354	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1355	t(i,k) = t(i,k)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
1356	enddo
1357	enddo
1358	!
1359	!
1360	!----------------------------------------------------------------
1361	! padding for small values
1362	!
1363	do k = kts, kte
1364	do i = its, ite
1365	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1366	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1367	enddo
1368	enddo
1369	enddo ! big loops
1370	END SUBROUTINE wsm62d
1371	! ...................................................................
1372	REAL FUNCTION rgmma(x)
1373	!-------------------------------------------------------------------
1374	IMPLICIT NONE
1375	!-------------------------------------------------------------------
1376	! rgmma function: use infinite product form
1377	REAL :: euler
1378	PARAMETER (euler=0.577215664901532)
1379	REAL :: x, y
1380	INTEGER :: i
1381	if(x.eq.1.)then
1382	rgmma=0.
1383	else
1384	rgmma=xexp(eulerx)
1385	do i=1,10000
1386	y=float(i)
1387	rgmma=rgmma(1.000+x/y)exp(-x/y)
1388	enddo
1389	rgmma=1./rgmma
1390	endif
1391	END FUNCTION rgmma
1392	!
1393	!--------------------------------------------------------------------------
1394	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1395	!--------------------------------------------------------------------------
1396	IMPLICIT NONE
1397	!--------------------------------------------------------------------------
1398	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1399	xai,xbi,ttp,tr
1400	INTEGER ice
1401	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1402	ttp=t0c+0.01
1403	dldt=cvap-cliq
1404	xa=-dldt/rv
1405	xb=xa+hvap/(rv*ttp)
1406	dldti=cvap-cice
1407	xai=-dldti/rv
1408	xbi=xai+hsub/(rv*ttp)
1409	tr=ttp/t
1410	if(t.lt.ttp.and.ice.eq.1) then
1411	fpvs=psat(trxai)exp(xbi*(1.-tr))
1412	else
1413	fpvs=psat(trxa)exp(xb*(1.-tr))
1414	endif
1415	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1416	END FUNCTION fpvs
1417	!-------------------------------------------------------------------
1418	SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,allowed_to_read)
1419	!-------------------------------------------------------------------
1420	IMPLICIT NONE
1421	!-------------------------------------------------------------------
1422	!.... constants which may not be tunable
1423	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
1424	LOGICAL, INTENT(IN) :: allowed_to_read
1425	!
1426	pi = 4.*atan(1.)
1427	xlv1 = cl-cpv
1428	!
1429	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
1430	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
1431	!
1432	bvtr1 = 1.+bvtr
1433	bvtr2 = 2.5+.5*bvtr
1434	bvtr3 = 3.+bvtr
1435	bvtr4 = 4.+bvtr
1436	bvtr6 = 6.+bvtr
1437	g1pbr = rgmma(bvtr1)
1438	g3pbr = rgmma(bvtr3)
1439	g4pbr = rgmma(bvtr4) ! 17.837825
1440	g6pbr = rgmma(bvtr6)
1441	g5pbro2 = rgmma(bvtr2) ! 1.8273
1442	pvtr = avtr*g4pbr/6.
1443	eacrr = 1.0
1444	pacrr = pin0ravtrg3pbr.25*eacrr
1445	precr1 = 2.pin0r*.78
1446	precr2 = 2.pin0r.31avtr*.5g5pbro2
1447	roqimax = 2.08e22dimax*8
1448	!
1449	bvts1 = 1.+bvts
1450	bvts2 = 2.5+.5*bvts
1451	bvts3 = 3.+bvts
1452	bvts4 = 4.+bvts
1453	g1pbs = rgmma(bvts1) !.8875
1454	g3pbs = rgmma(bvts3)
1455	g4pbs = rgmma(bvts4) ! 12.0786
1456	g5pbso2 = rgmma(bvts2)
1457	pvts = avts*g4pbs/6.
1458	pacrs = pin0savtsg3pbs.25
1459	precs1 = 4.n0s.65
1460	precs2 = 4.n0s.44avts.5g5pbso2
1461	pidn0r = pidenrn0r
1462	pidn0s = pidensn0s
1463	!
1464	pacrc = pin0savtsg3pbs.25*eacrc
1465	!
1466	bvtg1 = 1.+bvtg
1467	bvtg2 = 2.5+.5*bvtg
1468	bvtg3 = 3.+bvtg
1469	bvtg4 = 4.+bvtg
1470	g1pbg = rgmma(bvtg1)
1471	g3pbg = rgmma(bvtg3)
1472	g4pbg = rgmma(bvtg4)
1473	pacrg = pin0gavtgg3pbg.25
1474	g5pbgo2 = rgmma(bvtg2)
1475	pvtg = avtg*g4pbg/6.
1476	precg1 = 2.pin0g*.78
1477	precg2 = 2.pin0g.31avtg*.5g5pbgo2
1478	pidn0g = pidengn0g
1479	!
1480	rslopermax = 1./lamdarmax
1481	rslopesmax = 1./lamdasmax
1482	rslopegmax = 1./lamdagmax
1483	rsloperbmax = rslopermax ** bvtr
1484	rslopesbmax = rslopesmax ** bvts
1485	rslopegbmax = rslopegmax ** bvtg
1486	rsloper2max = rslopermax * rslopermax
1487	rslopes2max = rslopesmax * rslopesmax
1488	rslopeg2max = rslopegmax * rslopegmax
1489	rsloper3max = rsloper2max * rslopermax
1490	rslopes3max = rslopes2max * rslopesmax
1491	rslopeg3max = rslopeg2max * rslopegmax
1492	!
1493	END SUBROUTINE wsm6init
1494	!------------------------------------------------------------------------------
1495	subroutine slope_wsm6(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1496	vt,its,ite,kts,kte)
1497	IMPLICIT NONE
1498	INTEGER :: its,ite, jts,jte, kts,kte
1499	REAL, DIMENSION( its:ite , kts:kte,3) :: &
1500	qrs, &
1501	rslope, &
1502	rslopeb, &
1503	rslope2, &
1504	rslope3, &
1505	vt
1506	REAL, DIMENSION( its:ite , kts:kte) :: &
1507	den, &
1508	denfac, &
1509	t
1510	REAL, PARAMETER :: t0c = 273.15
1511	REAL, DIMENSION( its:ite , kts:kte ) :: &
1512	n0sfac
1513	REAL :: lamdar, lamdas, lamdag, x, y, z, supcol
1514	integer :: i, j, k
1515	!----------------------------------------------------------------
1516	! size distributions: (x=mixing ratio, y=air density):
1517	! valid for mixing ratio > 1.e-9 kg/kg.
1518	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
1519	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
1520	lamdag(x,y)= sqrt(sqrt(pidn0g/(xy))) ! (pidn0g/(xy))**.25
1521	!
1522	do k = kts, kte
1523	do i = its, ite
1524	supcol = t0c-t(i,k)
1525	!---------------------------------------------------------------
1526	! n0s: Intercept parameter for snow [m-4] [HDC 6]
1527	!---------------------------------------------------------------
1528	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1529	if(qrs(i,k,1).le.qcrmin)then
1530	rslope(i,k,1) = rslopermax
1531	rslopeb(i,k,1) = rsloperbmax
1532	rslope2(i,k,1) = rsloper2max
1533	rslope3(i,k,1) = rsloper3max
1534	else
1535	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
1536	rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1537	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1538	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1539	endif
1540	if(qrs(i,k,2).le.qcrmin)then
1541	rslope(i,k,2) = rslopesmax
1542	rslopeb(i,k,2) = rslopesbmax
1543	rslope2(i,k,2) = rslopes2max
1544	rslope3(i,k,2) = rslopes3max
1545	else
1546	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1547	rslopeb(i,k,2) = rslope(i,k,2)**bvts
1548	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1549	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1550	endif
1551	if(qrs(i,k,3).le.qcrmin)then
1552	rslope(i,k,3) = rslopegmax
1553	rslopeb(i,k,3) = rslopegbmax
1554	rslope2(i,k,3) = rslopeg2max
1555	rslope3(i,k,3) = rslopeg3max
1556	else
1557	rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
1558	rslopeb(i,k,3) = rslope(i,k,3)**bvtg
1559	rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
1560	rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
1561	endif
1562	vt(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)
1563	vt(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)
1564	vt(i,k,3) = pvtgrslopeb(i,k,3)denfac(i,k)
1565	if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1566	if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1567	if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
1568	enddo
1569	enddo
1570	END subroutine slope_wsm6
1571	!-----------------------------------------------------------------------------
1572	subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1573	vt,its,ite,kts,kte)
1574	IMPLICIT NONE
1575	INTEGER :: its,ite, jts,jte, kts,kte
1576	REAL, DIMENSION( its:ite , kts:kte) :: &
1577	qrs, &
1578	rslope, &
1579	rslopeb, &
1580	rslope2, &
1581	rslope3, &
1582	vt, &
1583	den, &
1584	denfac, &
1585	t
1586	REAL, PARAMETER :: t0c = 273.15
1587	REAL, DIMENSION( its:ite , kts:kte ) :: &
1588	n0sfac
1589	REAL :: lamdar, x, y, z, supcol
1590	integer :: i, j, k
1591	!----------------------------------------------------------------
1592	! size distributions: (x=mixing ratio, y=air density):
1593	! valid for mixing ratio > 1.e-9 kg/kg.
1594	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
1595	!
1596	do k = kts, kte
1597	do i = its, ite
1598	if(qrs(i,k).le.qcrmin)then
1599	rslope(i,k) = rslopermax
1600	rslopeb(i,k) = rsloperbmax
1601	rslope2(i,k) = rsloper2max
1602	rslope3(i,k) = rsloper3max
1603	else
1604	rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
1605	rslopeb(i,k) = rslope(i,k)**bvtr
1606	rslope2(i,k) = rslope(i,k)*rslope(i,k)
1607	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1608	endif
1609	vt(i,k) = pvtrrslopeb(i,k)denfac(i,k)
1610	if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1611	enddo
1612	enddo
1613	END subroutine slope_rain
1614	!------------------------------------------------------------------------------
1615	subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1616	vt,its,ite,kts,kte)
1617	IMPLICIT NONE
1618	INTEGER :: its,ite, jts,jte, kts,kte
1619	REAL, DIMENSION( its:ite , kts:kte) :: &
1620	qrs, &
1621	rslope, &
1622	rslopeb, &
1623	rslope2, &
1624	rslope3, &
1625	vt, &
1626	den, &
1627	denfac, &
1628	t
1629	REAL, PARAMETER :: t0c = 273.15
1630	REAL, DIMENSION( its:ite , kts:kte ) :: &
1631	n0sfac
1632	REAL :: lamdas, x, y, z, supcol
1633	integer :: i, j, k
1634	!----------------------------------------------------------------
1635	! size distributions: (x=mixing ratio, y=air density):
1636	! valid for mixing ratio > 1.e-9 kg/kg.
1637	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
1638	!
1639	do k = kts, kte
1640	do i = its, ite
1641	supcol = t0c-t(i,k)
1642	!---------------------------------------------------------------
1643	! n0s: Intercept parameter for snow [m-4] [HDC 6]
1644	!---------------------------------------------------------------
1645	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1646	if(qrs(i,k).le.qcrmin)then
1647	rslope(i,k) = rslopesmax
1648	rslopeb(i,k) = rslopesbmax
1649	rslope2(i,k) = rslopes2max
1650	rslope3(i,k) = rslopes3max
1651	else
1652	rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1653	rslopeb(i,k) = rslope(i,k)**bvts
1654	rslope2(i,k) = rslope(i,k)*rslope(i,k)
1655	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1656	endif
1657	vt(i,k) = pvtsrslopeb(i,k)denfac(i,k)
1658	if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1659	enddo
1660	enddo
1661	END subroutine slope_snow
1662	!----------------------------------------------------------------------------------
1663	subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1664	vt,its,ite,kts,kte)
1665	IMPLICIT NONE
1666	INTEGER :: its,ite, jts,jte, kts,kte
1667	REAL, DIMENSION( its:ite , kts:kte) :: &
1668	qrs, &
1669	rslope, &
1670	rslopeb, &
1671	rslope2, &
1672	rslope3, &
1673	vt, &
1674	den, &
1675	denfac, &
1676	t
1677	REAL, PARAMETER :: t0c = 273.15
1678	REAL, DIMENSION( its:ite , kts:kte ) :: &
1679	n0sfac
1680	REAL :: lamdag, x, y, z, supcol
1681	integer :: i, j, k
1682	!----------------------------------------------------------------
1683	! size distributions: (x=mixing ratio, y=air density):
1684	! valid for mixing ratio > 1.e-9 kg/kg.
1685	lamdag(x,y)= sqrt(sqrt(pidn0g/(xy))) ! (pidn0g/(xy))**.25
1686	!
1687	do k = kts, kte
1688	do i = its, ite
1689	!---------------------------------------------------------------
1690	! n0s: Intercept parameter for snow [m-4] [HDC 6]
1691	!---------------------------------------------------------------
1692	if(qrs(i,k).le.qcrmin)then
1693	rslope(i,k) = rslopegmax
1694	rslopeb(i,k) = rslopegbmax
1695	rslope2(i,k) = rslopeg2max
1696	rslope3(i,k) = rslopeg3max
1697	else
1698	rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
1699	rslopeb(i,k) = rslope(i,k)**bvtg
1700	rslope2(i,k) = rslope(i,k)*rslope(i,k)
1701	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1702	endif
1703	vt(i,k) = pvtgrslopeb(i,k)denfac(i,k)
1704	if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1705	enddo
1706	enddo
1707	END subroutine slope_graup
1708	!---------------------------------------------------------------------------------
1709	!-------------------------------------------------------------------
1710	SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1711	!-------------------------------------------------------------------
1712	!
1713	! for non-iteration semi-Lagrangain forward advection for cloud
1714	! with mass conservation and positive definite advection
1715	! 2nd order interpolation with monotonic piecewise linear method
1716	! this routine is under assumption of decfl < 1 for semi_Lagrangian
1717	!
1718	! dzl depth of model layer in meter
1719	! wwl terminal velocity at model layer m/s
1720	! rql cloud density*mixing ration
1721	! precip precipitation
1722	! dt time step
1723	! id kind of precip: 0 test case; 1 raindrop
1724	! iter how many time to guess mean terminal velocity: 0 pure forward.
1725	! 0 : use departure wind for advection
1726	! 1 : use mean wind for advection
1727	! > 1 : use mean wind after iter-1 iterations
1728	!
1729	! author: hann-ming henry juang <henry.juang@noaa.gov>
1730	! implemented by song-you hong
1731	!
1732	implicit none
1733	integer im,km,id
1734	real dt
1735	real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1736	real denl(im,km),denfacl(im,km),tkl(im,km)
1737	!
1738	integer i,k,n,m,kk,kb,kt,iter
1739	real tl,tl2,qql,dql,qqd
1740	real th,th2,qqh,dqh
1741	real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1742	real allold, allnew, zz, dzamin, cflmax, decfl
1743	real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1744	real den(km), denfac(km), tk(km)
1745	real wi(km+1), zi(km+1), za(km+1)
1746	real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1747	real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1748	!
1749	precip(:) = 0.0
1750	!
1751	i_loop : do i=1,im
1752	! -----------------------------------
1753	dz(:) = dzl(i,:)
1754	qq(:) = rql(i,:)
1755	ww(:) = wwl(i,:)
1756	den(:) = denl(i,:)
1757	denfac(:) = denfacl(i,:)
1758	tk(:) = tkl(i,:)
1759	! skip for no precipitation for all layers
1760	allold = 0.0
1761	do k=1,km
1762	allold = allold + qq(k)
1763	enddo
1764	if(allold.le.0.0) then
1765	cycle i_loop
1766	endif
1767	!
1768	! compute interface values
1769	zi(1)=0.0
1770	do k=1,km
1771	zi(k+1) = zi(k)+dz(k)
1772	enddo
1773	!
1774	! save departure wind
1775	wd(:) = ww(:)
1776	n=1
1777	100 continue
1778	! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1779	! 2nd order interpolation to get wi
1780	wi(1) = ww(1)
1781	wi(km+1) = ww(km)
1782	do k=2,km
1783	wi(k) = (ww(k)dz(k-1)+ww(k-1)dz(k))/(dz(k-1)+dz(k))
1784	enddo
1785	! 3rd order interpolation to get wi
1786	fa1 = 9./16.
1787	fa2 = 1./16.
1788	wi(1) = ww(1)
1789	wi(2) = 0.5*(ww(2)+ww(1))
1790	do k=3,km-1
1791	wi(k) = fa1(ww(k)+ww(k-1))-fa2(ww(k+1)+ww(k-2))
1792	enddo
1793	wi(km) = 0.5*(ww(km)+ww(km-1))
1794	wi(km+1) = ww(km)
1795	!
1796	! terminate of top of raingroup
1797	do k=2,km
1798	if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1799	enddo
1800	!
1801	! diffusivity of wi
1802	con1 = 0.05
1803	do k=km,1,-1
1804	decfl = (wi(k+1)-wi(k))*dt/dz(k)
1805	if( decfl .gt. con1 ) then
1806	wi(k) = wi(k+1) - con1*dz(k)/dt
1807	endif
1808	enddo
1809	! compute arrival point
1810	do k=1,km+1
1811	za(k) = zi(k) - wi(k)*dt
1812	enddo
1813	!
1814	do k=1,km
1815	dza(k) = za(k+1)-za(k)
1816	enddo
1817	dza(km+1) = zi(km+1) - za(km+1)
1818	!
1819	! computer deformation at arrival point
1820	do k=1,km
1821	qa(k) = qq(k)*dz(k)/dza(k)
1822	qr(k) = qa(k)/den(k)
1823	enddo
1824	qa(km+1) = 0.0
1825	! call maxmin(km,1,qa,' arrival points ')
1826	!
1827	! compute arrival terminal velocity, and estimate mean terminal velocity
1828	! then back to use mean terminal velocity
1829	if( n.le.iter ) then
1830	call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1831	if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1832	do k=1,km
1833	!#ifdef DEBUG
1834	! print,' slope_wsm3 ',qr(k)1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1835	!#endif
1836	! mean wind is average of departure and new arrival winds
1837	ww(k) = 0.5* ( wd(k)+wa(k) )
1838	enddo
1839	was(:) = wa(:)
1840	n=n+1
1841	go to 100
1842	endif
1843	!
1844	! estimate values at arrival cell interface with monotone
1845	do k=2,km
1846	dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1847	dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1848	if( dip*dim.le.0.0 ) then
1849	qmi(k)=qa(k)
1850	qpi(k)=qa(k)
1851	else
1852	qpi(k)=qa(k)+0.5(dip+dim)dza(k)
1853	qmi(k)=2.0*qa(k)-qpi(k)
1854	if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1855	qpi(k) = qa(k)
1856	qmi(k) = qa(k)
1857	endif
1858	endif
1859	enddo
1860	qpi(1)=qa(1)
1861	qmi(1)=qa(1)
1862	qmi(km+1)=qa(km+1)
1863	qpi(km+1)=qa(km+1)
1864	!
1865	! interpolation to regular point
1866	qn = 0.0
1867	kb=1
1868	kt=1
1869	intp : do k=1,km
1870	kb=max(kb-1,1)
1871	kt=max(kt-1,1)
1872	! find kb and kt
1873	if( zi(k).ge.za(km+1) ) then
1874	exit intp
1875	else
1876	find_kb : do kk=kb,km
1877	if( zi(k).le.za(kk+1) ) then
1878	kb = kk
1879	exit find_kb
1880	else
1881	cycle find_kb
1882	endif
1883	enddo find_kb
1884	find_kt : do kk=kt,km
1885	if( zi(k+1).le.za(kk) ) then
1886	kt = kk
1887	exit find_kt
1888	else
1889	cycle find_kt
1890	endif
1891	enddo find_kt
1892	kt = kt - 1
1893	! compute q with piecewise constant method
1894	if( kt.eq.kb ) then
1895	tl=(zi(k)-za(kb))/dza(kb)
1896	th=(zi(k+1)-za(kb))/dza(kb)
1897	tl2=tl*tl
1898	th2=th*th
1899	qqd=0.5*(qpi(kb)-qmi(kb))
1900	qqh=qqdth2+qmi(kb)th
1901	qql=qqdtl2+qmi(kb)tl
1902	qn(k) = (qqh-qql)/(th-tl)
1903	else if( kt.gt.kb ) then
1904	tl=(zi(k)-za(kb))/dza(kb)
1905	tl2=tl*tl
1906	qqd=0.5*(qpi(kb)-qmi(kb))
1907	qql=qqdtl2+qmi(kb)tl
1908	dql = qa(kb)-qql
1909	zsum = (1.-tl)*dza(kb)
1910	qsum = dql*dza(kb)
1911	if( kt-kb.gt.1 ) then
1912	do m=kb+1,kt-1
1913	zsum = zsum + dza(m)
1914	qsum = qsum + qa(m) * dza(m)
1915	enddo
1916	endif
1917	th=(zi(k+1)-za(kt))/dza(kt)
1918	th2=th*th
1919	qqd=0.5*(qpi(kt)-qmi(kt))
1920	dqh=qqdth2+qmi(kt)th
1921	zsum = zsum + th*dza(kt)
1922	qsum = qsum + dqh*dza(kt)
1923	qn(k) = qsum/zsum
1924	endif
1925	cycle intp
1926	endif
1927	!
1928	enddo intp
1929	!
1930	! rain out
1931	sum_precip: do k=1,km
1932	if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1933	precip(i) = precip(i) + qa(k)*dza(k)
1934	cycle sum_precip
1935	else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1936	precip(i) = precip(i) + qa(k)*(0.0-za(k))
1937	exit sum_precip
1938	endif
1939	exit sum_precip
1940	enddo sum_precip
1941	!
1942	! replace the new values
1943	rql(i,:) = qn(:)
1944	!
1945	! ----------------------------------
1946	enddo i_loop
1947	!
1948	END SUBROUTINE nislfv_rain_plm
1949	!-------------------------------------------------------------------
1950	SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
1951	!-------------------------------------------------------------------
1952	!
1953	! for non-iteration semi-Lagrangain forward advection for cloud
1954	! with mass conservation and positive definite advection
1955	! 2nd order interpolation with monotonic piecewise linear method
1956	! this routine is under assumption of decfl < 1 for semi_Lagrangian
1957	!
1958	! dzl depth of model layer in meter
1959	! wwl terminal velocity at model layer m/s
1960	! rql cloud density*mixing ration
1961	! precip precipitation
1962	! dt time step
1963	! id kind of precip: 0 test case; 1 raindrop
1964	! iter how many time to guess mean terminal velocity: 0 pure forward.
1965	! 0 : use departure wind for advection
1966	! 1 : use mean wind for advection
1967	! > 1 : use mean wind after iter-1 iterations
1968	!
1969	! author: hann-ming henry juang <henry.juang@noaa.gov>
1970	! implemented by song-you hong
1971	!
1972	implicit none
1973	integer im,km,id
1974	real dt
1975	real dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
1976	real denl(im,km),denfacl(im,km),tkl(im,km)
1977	!
1978	integer i,k,n,m,kk,kb,kt,iter,ist
1979	real tl,tl2,qql,dql,qqd
1980	real th,th2,qqh,dqh
1981	real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1982	real allold, allnew, zz, dzamin, cflmax, decfl
1983	real dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
1984	real den(km), denfac(km), tk(km)
1985	real wi(km+1), zi(km+1), za(km+1)
1986	real qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1987	real dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
1988	!
1989	precip(:) = 0.0
1990	precip1(:) = 0.0
1991	precip2(:) = 0.0
1992	!
1993	i_loop : do i=1,im
1994	! -----------------------------------
1995	dz(:) = dzl(i,:)
1996	qq(:) = rql(i,:)
1997	qq2(:) = rql2(i,:)
1998	ww(:) = wwl(i,:)
1999	den(:) = denl(i,:)
2000	denfac(:) = denfacl(i,:)
2001	tk(:) = tkl(i,:)
2002	! skip for no precipitation for all layers
2003	allold = 0.0
2004	do k=1,km
2005	allold = allold + qq(k)
2006	enddo
2007	if(allold.le.0.0) then
2008	cycle i_loop
2009	endif
2010	!
2011	! compute interface values
2012	zi(1)=0.0
2013	do k=1,km
2014	zi(k+1) = zi(k)+dz(k)
2015	enddo
2016	!
2017	! save departure wind
2018	wd(:) = ww(:)
2019	n=1
2020	100 continue
2021	! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2022	! 2nd order interpolation to get wi
2023	wi(1) = ww(1)
2024	wi(km+1) = ww(km)
2025	do k=2,km
2026	wi(k) = (ww(k)dz(k-1)+ww(k-1)dz(k))/(dz(k-1)+dz(k))
2027	enddo
2028	! 3rd order interpolation to get wi
2029	fa1 = 9./16.
2030	fa2 = 1./16.
2031	wi(1) = ww(1)
2032	wi(2) = 0.5*(ww(2)+ww(1))
2033	do k=3,km-1
2034	wi(k) = fa1(ww(k)+ww(k-1))-fa2(ww(k+1)+ww(k-2))
2035	enddo
2036	wi(km) = 0.5*(ww(km)+ww(km-1))
2037	wi(km+1) = ww(km)
2038	!
2039	! terminate of top of raingroup
2040	do k=2,km
2041	if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2042	enddo
2043	!
2044	! diffusivity of wi
2045	con1 = 0.05
2046	do k=km,1,-1
2047	decfl = (wi(k+1)-wi(k))*dt/dz(k)
2048	if( decfl .gt. con1 ) then
2049	wi(k) = wi(k+1) - con1*dz(k)/dt
2050	endif
2051	enddo
2052	! compute arrival point
2053	do k=1,km+1
2054	za(k) = zi(k) - wi(k)*dt
2055	enddo
2056	!
2057	do k=1,km
2058	dza(k) = za(k+1)-za(k)
2059	enddo
2060	dza(km+1) = zi(km+1) - za(km+1)
2061	!
2062	! computer deformation at arrival point
2063	do k=1,km
2064	qa(k) = qq(k)*dz(k)/dza(k)
2065	qa2(k) = qq2(k)*dz(k)/dza(k)
2066	qr(k) = qa(k)/den(k)
2067	qr2(k) = qa2(k)/den(k)
2068	enddo
2069	qa(km+1) = 0.0
2070	qa2(km+1) = 0.0
2071	! call maxmin(km,1,qa,' arrival points ')
2072	!
2073	! compute arrival terminal velocity, and estimate mean terminal velocity
2074	! then back to use mean terminal velocity
2075	if( n.le.iter ) then
2076	call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2077	call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
2078	do k = 1, km
2079	tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
2080	IF ( tmp(k) .gt. 1.e-15 ) THEN
2081	wa(k) = (wa(k)qr(k) + wa2(k)qr2(k))/tmp(k)
2082	ELSE
2083	wa(k) = 0.
2084	ENDIF
2085	enddo
2086	if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2087	do k=1,km
2088	!#ifdef DEBUG
2089	! print,' slope_wsm3 ',qr(k)1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
2090	! ww(k),wa(k)
2091	!#endif
2092	! mean wind is average of departure and new arrival winds
2093	ww(k) = 0.5* ( wd(k)+wa(k) )
2094	enddo
2095	was(:) = wa(:)
2096	n=n+1
2097	go to 100
2098	endif
2099	ist_loop : do ist = 1, 2
2100	if (ist.eq.2) then
2101	qa(:) = qa2(:)
2102	endif
2103	!
2104	precip(i) = 0.
2105	!
2106	! estimate values at arrival cell interface with monotone
2107	do k=2,km
2108	dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2109	dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2110	if( dip*dim.le.0.0 ) then
2111	qmi(k)=qa(k)
2112	qpi(k)=qa(k)
2113	else
2114	qpi(k)=qa(k)+0.5(dip+dim)dza(k)
2115	qmi(k)=2.0*qa(k)-qpi(k)
2116	if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2117	qpi(k) = qa(k)
2118	qmi(k) = qa(k)
2119	endif
2120	endif
2121	enddo
2122	qpi(1)=qa(1)
2123	qmi(1)=qa(1)
2124	qmi(km+1)=qa(km+1)
2125	qpi(km+1)=qa(km+1)
2126	!
2127	! interpolation to regular point
2128	qn = 0.0
2129	kb=1
2130	kt=1
2131	intp : do k=1,km
2132	kb=max(kb-1,1)
2133	kt=max(kt-1,1)
2134	! find kb and kt
2135	if( zi(k).ge.za(km+1) ) then
2136	exit intp
2137	else
2138	find_kb : do kk=kb,km
2139	if( zi(k).le.za(kk+1) ) then
2140	kb = kk
2141	exit find_kb
2142	else
2143	cycle find_kb
2144	endif
2145	enddo find_kb
2146	find_kt : do kk=kt,km
2147	if( zi(k+1).le.za(kk) ) then
2148	kt = kk
2149	exit find_kt
2150	else
2151	cycle find_kt
2152	endif
2153	enddo find_kt
2154	kt = kt - 1
2155	! compute q with piecewise constant method
2156	if( kt.eq.kb ) then
2157	tl=(zi(k)-za(kb))/dza(kb)
2158	th=(zi(k+1)-za(kb))/dza(kb)
2159	tl2=tl*tl
2160	th2=th*th
2161	qqd=0.5*(qpi(kb)-qmi(kb))
2162	qqh=qqdth2+qmi(kb)th
2163	qql=qqdtl2+qmi(kb)tl
2164	qn(k) = (qqh-qql)/(th-tl)
2165	else if( kt.gt.kb ) then
2166	tl=(zi(k)-za(kb))/dza(kb)
2167	tl2=tl*tl
2168	qqd=0.5*(qpi(kb)-qmi(kb))
2169	qql=qqdtl2+qmi(kb)tl
2170	dql = qa(kb)-qql
2171	zsum = (1.-tl)*dza(kb)
2172	qsum = dql*dza(kb)
2173	if( kt-kb.gt.1 ) then
2174	do m=kb+1,kt-1
2175	zsum = zsum + dza(m)
2176	qsum = qsum + qa(m) * dza(m)
2177	enddo
2178	endif
2179	th=(zi(k+1)-za(kt))/dza(kt)
2180	th2=th*th
2181	qqd=0.5*(qpi(kt)-qmi(kt))
2182	dqh=qqdth2+qmi(kt)th
2183	zsum = zsum + th*dza(kt)
2184	qsum = qsum + dqh*dza(kt)
2185	qn(k) = qsum/zsum
2186	endif
2187	cycle intp
2188	endif
2189	!
2190	enddo intp
2191	!
2192	! rain out
2193	sum_precip: do k=1,km
2194	if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2195	precip(i) = precip(i) + qa(k)*dza(k)
2196	cycle sum_precip
2197	else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2198	precip(i) = precip(i) + qa(k)*(0.0-za(k))
2199	exit sum_precip
2200	endif
2201	exit sum_precip
2202	enddo sum_precip
2203	!
2204	! replace the new values
2205	if(ist.eq.1) then
2206	rql(i,:) = qn(:)
2207	precip1(i) = precip(i)
2208	else
2209	rql2(i,:) = qn(:)
2210	precip2(i) = precip(i)
2211	endif
2212	enddo ist_loop
2213	!
2214	! ----------------------------------
2215	enddo i_loop
2216	!
2217	END SUBROUTINE nislfv_rain_plm6
2218	END MODULE module_mp_wsm6

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

Download in other formats:

Original Format