Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wdm6.F @ 14

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

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

Download in other formats:

Original Format