Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm5.F @ 3547

Last change on this file since 3547 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 51.1 KB

Line
1	#if ( RWORDSIZE == 4 )
2	# define VREC vsrec
3	# define VSQRT vssqrt
4	#else
5	# define VREC vrec
6	# define VSQRT vsqrt
7	#endif
8
9	!Including inline expansion statistical function
10	MODULE module_mp_wsm5
11	!
12	!
13	REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
14	REAL, PARAMETER, PRIVATE :: n0r = 8.e6
15	REAL, PARAMETER, PRIVATE :: avtr = 841.9
16	REAL, PARAMETER, PRIVATE :: bvtr = 0.8
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
22	REAL, PARAMETER, PRIVATE :: bvts = .41
23	REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
24	REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
25	REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
26	REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
27	REAL, PARAMETER, PRIVATE :: betai = .6
28	REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
29	REAL, PARAMETER, PRIVATE :: dicon = 11.9
30	REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
31	REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
32	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
33	REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
34	REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
35	REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
36	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
37	REAL, PARAMETER, PRIVATE :: t40c = 233.16
38	REAL, PARAMETER, PRIVATE :: eacrc = 1.0
39	REAL, SAVE :: &
40	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
41	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
42	precr1,precr2,xm0,xmmax,roqimax,bvts1, &
43	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
44	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
45	pidn0s,xlv1,pacrc, &
46	rslopermax,rslopesmax,rslopegmax, &
47	rsloperbmax,rslopesbmax,rslopegbmax, &
48	rsloper2max,rslopes2max,rslopeg2max, &
49	rsloper3max,rslopes3max,rslopeg3max
50	!
51	! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
52	!
53	CONTAINS
54	!===================================================================
55	!
56	SUBROUTINE wsm5(th, q, qc, qr, qi, qs &
57	,den, pii, p, delz &
58	,delt,g, cpd, cpv, rd, rv, t0c &
59	,ep1, ep2, qmin &
60	,XLS, XLV0, XLF0, den0, denr &
61	,cliq,cice,psat &
62	,rain, rainncv &
63	,snow, snowncv &
64	,sr &
65	,ids,ide, jds,jde, kds,kde &
66	,ims,ime, jms,jme, kms,kme &
67	,its,ite, jts,jte, kts,kte &
68	)
69	!-------------------------------------------------------------------
70	IMPLICIT NONE
71	!-------------------------------------------------------------------
72	!
73	! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
74	! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
75	! number concentration is a function of temperature, and seperate assumption
76	! is developed, in which ice crystal number concentration is a function
77	! of ice amount. A theoretical background of the ice-microphysics and related
78	! processes in the WSMMPs are described in Hong et al. (2004).
79	! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
80	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
81	!
82	! WSM5 cloud scheme
83	!
84	! Coded by Song-You Hong (Yonsei Univ.)
85	! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
86	! Summer 2002
87	!
88	! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
89	! Summer 2003
90	!
91	! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
92	! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
93	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
94	!
95	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
96	ims,ime, jms,jme, kms,kme , &
97	its,ite, jts,jte, kts,kte
98	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
99	INTENT(INOUT) :: &
100	th, &
101	q, &
102	qc, &
103	qi, &
104	qr, &
105	qs
106	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
107	INTENT(IN ) :: &
108	den, &
109	pii, &
110	p, &
111	delz
112	REAL, INTENT(IN ) :: delt, &
113	g, &
114	rd, &
115	rv, &
116	t0c, &
117	den0, &
118	cpd, &
119	cpv, &
120	ep1, &
121	ep2, &
122	qmin, &
123	XLS, &
124	XLV0, &
125	XLF0, &
126	cliq, &
127	cice, &
128	psat, &
129	denr
130	REAL, DIMENSION( ims:ime , jms:jme ), &
131	INTENT(INOUT) :: rain, &
132	rainncv, &
133	sr
134
135	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
136	INTENT(INOUT) :: snow, &
137	snowncv
138
139	! LOCAL VAR
140	REAL, DIMENSION( its:ite , kts:kte ) :: t
141	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
142	INTEGER :: i,j,k
143	!-------------------------------------------------------------------
144	DO j=jts,jte
145	DO k=kts,kte
146	DO i=its,ite
147	t(i,k)=th(i,k,j)*pii(i,k,j)
148	qci(i,k,1) = qc(i,k,j)
149	qci(i,k,2) = qi(i,k,j)
150	qrs(i,k,1) = qr(i,k,j)
151	qrs(i,k,2) = qs(i,k,j)
152	ENDDO
153	ENDDO
154	CALL wsm52D(t, q(ims,kms,j), qci, qrs &
155	,den(ims,kms,j) &
156	,p(ims,kms,j), delz(ims,kms,j) &
157	,delt,g, cpd, cpv, rd, rv, t0c &
158	,ep1, ep2, qmin &
159	,XLS, XLV0, XLF0, den0, denr &
160	,cliq,cice,psat &
161	,j &
162	,rain(ims,j),rainncv(ims,j) &
163	,sr(ims,j) &
164	,ids,ide, jds,jde, kds,kde &
165	,ims,ime, jms,jme, kms,kme &
166	,its,ite, jts,jte, kts,kte &
167	,snow(ims,j),snowncv(ims,j) &
168	)
169	DO K=kts,kte
170	DO I=its,ite
171	th(i,k,j)=t(i,k)/pii(i,k,j)
172	qc(i,k,j) = qci(i,k,1)
173	qi(i,k,j) = qci(i,k,2)
174	qr(i,k,j) = qrs(i,k,1)
175	qs(i,k,j) = qrs(i,k,2)
176	ENDDO
177	ENDDO
178	ENDDO
179	END SUBROUTINE wsm5
180	!===================================================================
181	!
182	SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz &
183	,delt,g, cpd, cpv, rd, rv, t0c &
184	,ep1, ep2, qmin &
185	,XLS, XLV0, XLF0, den0, denr &
186	,cliq,cice,psat &
187	,lat &
188	,rain,rainncv &
189	,sr &
190	,ids,ide, jds,jde, kds,kde &
191	,ims,ime, jms,jme, kms,kme &
192	,its,ite, jts,jte, kts,kte &
193	,snow,snowncv &
194	)
195	!-------------------------------------------------------------------
196	IMPLICIT NONE
197	!-------------------------------------------------------------------
198	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
199	ims,ime, jms,jme, kms,kme , &
200	its,ite, jts,jte, kts,kte, &
201	lat
202	REAL, DIMENSION( its:ite , kts:kte ), &
203	INTENT(INOUT) :: &
204	t
205	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
206	INTENT(INOUT) :: &
207	qci, &
208	qrs
209	REAL, DIMENSION( ims:ime , kms:kme ), &
210	INTENT(INOUT) :: &
211	q
212	REAL, DIMENSION( ims:ime , kms:kme ), &
213	INTENT(IN ) :: &
214	den, &
215	p, &
216	delz
217	REAL, INTENT(IN ) :: delt, &
218	g, &
219	cpd, &
220	cpv, &
221	t0c, &
222	den0, &
223	rd, &
224	rv, &
225	ep1, &
226	ep2, &
227	qmin, &
228	XLS, &
229	XLV0, &
230	XLF0, &
231	cliq, &
232	cice, &
233	psat, &
234	denr
235	REAL, DIMENSION( ims:ime ), &
236	INTENT(INOUT) :: rain, &
237	rainncv, &
238	sr
239
240	REAL, DIMENSION( ims:ime ), OPTIONAL, &
241	INTENT(INOUT) :: snow, &
242	snowncv
243
244	! LOCAL VAR
245	REAL, DIMENSION( its:ite , kts:kte , 2) :: &
246	rh, qs, rslope, rslope2, rslope3, rslopeb, &
247	falk, fall, work1
248	REAL, DIMENSION( its:ite , kts:kte ) :: &
249	falkc, work1c, work2c, fallc
250	REAL, DIMENSION( its:ite , kts:kte ) :: &
251	praut, psaut, prevp, psdep, pracw, psaci, psacw, &
252	pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, denfac, xni,&
253	n0sfac
254	! variables for optimization
255	REAL, DIMENSION( its:ite ) :: tvec1
256	INTEGER, DIMENSION( its:ite ) :: mstep, numdt
257	REAL, DIMENSION(its:ite) :: rmstep
258	REAL dtcldden, rdelz, rdtcld
259	LOGICAL, DIMENSION( its:ite ) :: flgcld
260	REAL :: pi, &
261	cpmcal, xlcal, lamdar, lamdas, diffus, &
262	viscos, xka, venfac, conden, diffac, &
263	x, y, z, a, b, c, d, e, &
264	qdt, holdrr, holdrs, supcol, pvt, &
265	coeres, supsat, dtcld, xmi, eacrs, satdt, &
266	vt2i,vt2s,acrfac, &
267	qimax, diameter, xni0, roqi0, &
268	fallsum, fallsum_qsi, xlwork2, factor, source, &
269	value, xlf, pfrzdtc, pfrzdtr, supice
270	REAL :: temp
271	REAL :: holdc, holdci
272	INTEGER :: i, j, k, mstepmax, &
273	iprt, latd, lond, loop, loops, ifsat, n
274	! Temporaries used for inlining fpvs function
275	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
276	!
277	!=================================================================
278	! compute internal functions
279	!
280	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
281	xlcal(x) = xlv0-xlv1*(x-t0c)
282	!----------------------------------------------------------------
283	! size distributions: (x=mixing ratio, y=air density):
284	! valid for mixing ratio > 1.e-9 kg/kg.
285	!
286	! Optimizatin : A*B => exp(log(A)(B))
287	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
288	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
289	!
290	!----------------------------------------------------------------
291	! diffus: diffusion coefficient of the water vapor
292	! viscos: kinematic viscosity(m2s-1)
293	! diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
294	! viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
295	! xka(x,y) = 1.414e3viscos(x,y)y
296	! diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
297	! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
298	! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
299	! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
300	!
301	!
302	pi = 4. * atan(1.)
303	!
304	!----------------------------------------------------------------
305	! paddint 0 for negative values generated by dynamics
306	!
307	do k = kts, kte
308	do i = its, ite
309	qci(i,k,1) = max(qci(i,k,1),0.0)
310	qrs(i,k,1) = max(qrs(i,k,1),0.0)
311	qci(i,k,2) = max(qci(i,k,2),0.0)
312	qrs(i,k,2) = max(qrs(i,k,2),0.0)
313	enddo
314	enddo
315	!
316	!----------------------------------------------------------------
317	! latent heat for phase changes and heat capacity. neglect the
318	! changes during microphysical process calculation
319	! emanuel(1994)
320	!
321	do k = kts, kte
322	do i = its, ite
323	cpm(i,k) = cpmcal(q(i,k))
324	xl(i,k) = xlcal(t(i,k))
325	enddo
326	enddo
327	!
328	!----------------------------------------------------------------
329	! compute the minor time steps.
330	!
331	loops = max(nint(delt/dtcldcr),1)
332	dtcld = delt/loops
333	if(delt.le.dtcldcr) dtcld = delt
334	!
335	do loop = 1,loops
336	!
337	!----------------------------------------------------------------
338	! initialize the large scale variables
339	!
340	do i = its, ite
341	mstep(i) = 1
342	flgcld(i) = .true.
343	enddo
344	!
345	! do k = kts, kte
346	! do i = its, ite
347	! denfac(i,k) = sqrt(den0/den(i,k))
348	! enddo
349	! enddo
350	do k = kts, kte
351	CALL VREC( tvec1(its), den(its,k), ite-its+1)
352	do i = its, ite
353	tvec1(i) = tvec1(i)*den0
354	enddo
355	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
356	enddo
357	!
358	! Inline expansion for fpvs
359	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
360	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
361	hsub = xls
362	hvap = xlv0
363	cvap = cpv
364	ttp=t0c+0.01
365	dldt=cvap-cliq
366	xa=-dldt/rv
367	xb=xa+hvap/(rv*ttp)
368	dldti=cvap-cice
369	xai=-dldti/rv
370	xbi=xai+hsub/(rv*ttp)
371	do k = kts, kte
372	do i = its, ite
373	tr=ttp/t(i,k)
374	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
375	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
376	qs(i,k,1) = max(qs(i,k,1),qmin)
377	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
378	tr=ttp/t(i,k)
379	if(t(i,k).lt.ttp) then
380	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
381	else
382	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
383	endif
384	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
385	qs(i,k,2) = max(qs(i,k,2),qmin)
386	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
387	enddo
388	enddo
389	!
390	!----------------------------------------------------------------
391	! initialize the variables for microphysical physics
392	!
393	!
394	do k = kts, kte
395	do i = its, ite
396	prevp(i,k) = 0.
397	psdep(i,k) = 0.
398	praut(i,k) = 0.
399	psaut(i,k) = 0.
400	pracw(i,k) = 0.
401	psaci(i,k) = 0.
402	psacw(i,k) = 0.
403	pigen(i,k) = 0.
404	pidep(i,k) = 0.
405	pcond(i,k) = 0.
406	psmlt(i,k) = 0.
407	psevp(i,k) = 0.
408	falk(i,k,1) = 0.
409	falk(i,k,2) = 0.
410	fall(i,k,1) = 0.
411	fall(i,k,2) = 0.
412	fallc(i,k) = 0.
413	falkc(i,k) = 0.
414	xni(i,k) = 1.e3
415	enddo
416	enddo
417	!
418	!----------------------------------------------------------------
419	! compute the fallout term:
420	! first, vertical terminal velosity for minor loops
421	!
422	do k = kts, kte
423	do i = its, ite
424	supcol = t0c-t(i,k)
425	!---------------------------------------------------------------
426	! n0s: Intercept parameter for snow [m-4] [HDC 6]
427	!---------------------------------------------------------------
428	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
429	if(qrs(i,k,1).le.qcrmin)then
430	rslope(i,k,1) = rslopermax
431	rslopeb(i,k,1) = rsloperbmax
432	rslope2(i,k,1) = rsloper2max
433	rslope3(i,k,1) = rsloper3max
434	else
435	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
436	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
437	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
438	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
439	endif
440	if(qrs(i,k,2).le.qcrmin)then
441	rslope(i,k,2) = rslopesmax
442	rslopeb(i,k,2) = rslopesbmax
443	rslope2(i,k,2) = rslopes2max
444	rslope3(i,k,2) = rslopes3max
445	else
446	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
447	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
448	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
449	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
450	endif
451	!-------------------------------------------------------------
452	! Ni: ice crystal number concentraiton [HDC 5c]
453	!-------------------------------------------------------------
454	! xni(i,k) = min(max(5.38e7*(den(i,k) &
455	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
456	temp = (den(i,k)*max(qci(i,k,2),qmin))
457	temp = sqrt(sqrt(temptemptemp))
458	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
459	enddo
460	enddo
461	!
462	mstepmax = 1
463	numdt = 1
464	do k = kte, kts, -1
465	do i = its, ite
466	work1(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)/delz(i,k)
467	work1(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)/delz(i,k)
468	numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
469	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
470	enddo
471	enddo
472	do i = its, ite
473	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
474	rmstep(i) = 1./mstep(i)
475	enddo
476	!
477	do n = 1, mstepmax
478	k = kte
479	do i = its, ite
480	if(n.le.mstep(i)) then
481	! falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
482	! falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)/mstep(i)
483	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)*rmstep(i)
484	falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)*rmstep(i)
485	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
486	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
487	! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
488	! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
489	dtcldden = dtcld/den(i,k)
490	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
491	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
492	endif
493	enddo
494	do k = kte-1, kts, -1
495	do i = its, ite
496	if(n.le.mstep(i)) then
497	! falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
498	! falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)/mstep(i)
499	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)*rmstep(i)
500	falk(i,k,2) = den(i,k)qrs(i,k,2)work1(i,k,2)*rmstep(i)
501	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
502	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
503	! qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
504	! delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
505	! qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
506	! delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
507	dtcldden = dtcld/den(i,k)
508	rdelz = 1./delz(i,k)
509	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
510	delz(i,k+1)rdelz)*dtcldden,0.)
511	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
512	delz(i,k+1)rdelz)*dtcldden,0.)
513	endif
514	enddo
515	enddo
516	do k = kte, kts, -1
517	do i = its, ite
518	if(n.le.mstep(i)) then
519	if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
520	!----------------------------------------------------------------
521	! psmlt: melting of snow [HL A33] [RH83 A25]
522	! (T>T0: S->R)
523	!----------------------------------------------------------------
524	xlf = xlf0
525	! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
526	work2(i,k)= (exp(log(((1.496e-6((t(i,k))sqrt(t(i,k))) &
527	/((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
528	exp(log(t(i,k))(1.81))/p(i,k)))) &
529	((.3333333)))/sqrt((1.496e-6((t(i,k)) &
530	*sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
531	*sqrt(sqrt(den0/(den(i,k)))))
532	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
533	! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
534	! n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
535	! work2(i,k)coeres)
536	psmlt(i,k) = &
537	(1.414e3(1.496e-6 ((t(i,k))sqrt(t(i,k))) /((t(i,k))+120.)/(den(i,k)) )(den(i,k)))&
538	/xlf(t0c-t(i,k))pi/2. &
539	n0sfac(i,k)(precs1*rslope2(i,k,2)+precs2 &
540	work2(i,k)coeres)
541	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
542	-qrs(i,k,2)/mstep(i)),0.)
543	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
544	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
545	t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
546	endif
547	endif
548	enddo
549	enddo
550	enddo
551	!---------------------------------------------------------------
552	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
553	!---------------------------------------------------------------
554	mstepmax = 1
555	mstep = 1
556	numdt = 1
557	do k = kte, kts, -1
558	do i = its, ite
559	if(qci(i,k,2).le.0.) then
560	work2c(i,k) = 0.
561	else
562	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
563	! diameter = min(dicon * sqrt(xmi),dimax)
564	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
565	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
566	work2c(i,k) = work1c(i,k)/delz(i,k)
567	endif
568	numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
569	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
570	enddo
571	enddo
572	do i = its, ite
573	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
574	enddo
575	!
576	do n = 1, mstepmax
577	k = kte
578	do i = its, ite
579	if(n.le.mstep(i)) then
580	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
581	holdc = falkc(i,k)
582	fallc(i,k) = fallc(i,k)+falkc(i,k)
583	holdci = qci(i,k,2)
584	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
585	endif
586	enddo
587	do k = kte-1, kts, -1
588	do i = its, ite
589	if(n.le.mstep(i)) then
590	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
591	holdc = falkc(i,k)
592	fallc(i,k) = fallc(i,k)+falkc(i,k)
593	holdci = qci(i,k,2)
594	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
595	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
596	endif
597	enddo
598	enddo
599	enddo
600	!
601	!
602	!----------------------------------------------------------------
603	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
604	!
605	do i = its, ite
606	fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
607	fallsum_qsi = fall(i,1,2)+fallc(i,1)
608	rainncv(i) = 0.
609	if(fallsum.gt.0.) then
610	rainncv(i) = fallsumdelz(i,1)/denrdtcld*1000.
611	rain(i) = fallsumdelz(i,1)/denrdtcld*1000. + rain(i)
612	endif
613	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
614	snowncv(i) = 0.
615	if(fallsum_qsi.gt.0.) then
616	snowncv(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000.
617	snow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i)
618	endif
619	ENDIF
620	sr(i) = 0.
621	if(fallsum.gt.0.)sr(i)=fallsum_qsidelz(i,kts)/denrdtcld*1000./(rainncv(i)+1.e-12)
622	enddo
623	!
624	!---------------------------------------------------------------
625	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
626	! (T>T0: I->C)
627	!---------------------------------------------------------------
628	do k = kts, kte
629	do i = its, ite
630	supcol = t0c-t(i,k)
631	xlf = xls-xl(i,k)
632	if(supcol.lt.0.) xlf = xlf0
633	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
634	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
635	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
636	qci(i,k,2) = 0.
637	endif
638	!---------------------------------------------------------------
639	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
640	! (T<-40C: C->I)
641	!---------------------------------------------------------------
642	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
643	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
644	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
645	qci(i,k,1) = 0.
646	endif
647	!---------------------------------------------------------------
648	! pihtf: heterogeneous freezing of cloud water [HL A44]
649	! (T0>T>-40C: C->I)
650	!---------------------------------------------------------------
651	if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
652	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
653	! den(i,k)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
654	pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
655	den(i,k)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
656	qci(i,k,2) = qci(i,k,2) + pfrzdtc
657	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
658	qci(i,k,1) = qci(i,k,1)-pfrzdtc
659	endif
660	!---------------------------------------------------------------
661	! psfrz: freezing of rain water [HL A20] [LFO 45]
662	! (T<T0, R->S)
663	!---------------------------------------------------------------
664	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
665	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k) &
666	! (exp(pfrz2supcol)-1.)rslope(i,k,1)7dtcld, &
667	! qrs(i,k,1))
668	temp = rslope(i,k,1)
669	temp = temptemptemptemptemptemptemp
670	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k) &
671	(exp(pfrz2supcol)-1.)tempdtcld, &
672	qrs(i,k,1))
673	qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
674	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
675	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
676	endif
677	enddo
678	enddo
679	!
680	!----------------------------------------------------------------
681	! rsloper: reverse of the slope parameter of the rain(m)
682	! xka: thermal conductivity of air(jm-1s-1k-1)
683	! work1: the thermodynamic term in the denominator associated with
684	! heat conduction and vapor diffusion
685	! (ry88, y93, h85)
686	! work2: parameter associated with the ventilation effects(y93)
687	!
688	do k = kts, kte
689	do i = its, ite
690	if(qrs(i,k,1).le.qcrmin)then
691	rslope(i,k,1) = rslopermax
692	rslopeb(i,k,1) = rsloperbmax
693	rslope2(i,k,1) = rsloper2max
694	rslope3(i,k,1) = rsloper3max
695	else
696	! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
697	rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
698	rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
699	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
700	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
701	endif
702	if(qrs(i,k,2).le.qcrmin)then
703	rslope(i,k,2) = rslopesmax
704	rslopeb(i,k,2) = rslopesbmax
705	rslope2(i,k,2) = rslopes2max
706	rslope3(i,k,2) = rslopes3max
707	else
708	! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
709	rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s(n0sfac(i,k))/((qrs(i,k,2))(den(i,k))))))
710	rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
711	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
712	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
713	endif
714	enddo
715	enddo
716	!
717	do k = kts, kte
718	do i = its, ite
719	! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
720	work1(i,k,1) = &
721	((((den(i,k))(xl(i,k))(xl(i,k))) * ((t(i,k))+120.) * (den(i,k))) &
722	/ &
723	( 1.414e3 * (1.496e-6 * ((t(i,k))sqrt(t(i,k)))) (den(i,k)) * &
724	(rv(t(i,k))(t(i,k))))) &
725	+ &
726	p(i,k) / ( (qs(i,k,1)) * ( 8.794e-5 * exp(log(t(i,k))*(1.81)) ) )
727	! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
728	work1(i,k,2) = &
729	( &
730	(((den(i,k))(xls)(xls))((t(i,k))+120.)(den(i,k))) &
731	/ &
732	( &
733	1.414e3 * (1.496e-6 * ((t(i,k))sqrt(t(i,k)))) (den(i,k)) * &
734	(rv(t(i,k))(t(i,k))) &
735	) &
736	+ &
737	p(i,k) &
738	/ &
739	( qs(i,k,2) * (8.794e-5 * exp(log(t(i,k))*(1.81)))) &
740	)
741	! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
742	work2(i,k) = &
743	( &
744	exp(.3333333*log( &
745	((1.496e-6 * ((t(i,k))sqrt(t(i,k))))p(i,k)) &
746	/ &
747	(((t(i,k))+120.)den(i,k)(8.794e-5 * exp(log(t(i,k))*(1.81)))) &
748	)) &
749	* &
750	sqrt(sqrt(den0/(den(i,k)))) &
751	) &
752	/ &
753	sqrt( &
754	(1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
755	/ &
756	( &
757	((t(i,k))+120.) * den(i,k) &
758	) &
759	)
760	ENDDO
761	ENDDO
762	!
763	!===============================================================
764	!
765	! warm rain processes
766	!
767	! - follows the processes in RH83 and LFO except for autoconcersion
768	!
769	!===============================================================
770	!
771	do k = kts, kte
772	do i = its, ite
773	supsat = max(q(i,k),qmin)-qs(i,k,1)
774	satdt = supsat/dtcld
775	!---------------------------------------------------------------
776	! praut: auto conversion rate from cloud to rain [HDC 16]
777	! (C->R)
778	!---------------------------------------------------------------
779	if(qci(i,k,1).gt.qc0) then
780	praut(i,k) = qck1exp(log(qci(i,k,1))((7./3.)))
781	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
782	endif
783	!---------------------------------------------------------------
784	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
785	! (C->R)
786	!---------------------------------------------------------------
787	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
788	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
789	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
790	endif
791	!---------------------------------------------------------------
792	! prevp: evaporation/condensation rate of rain [HDC 14]
793	! (V->R or R->V)
794	!---------------------------------------------------------------
795	if(qrs(i,k,1).gt.0.) then
796	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
797	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
798	+precr2work2(i,k)coeres)/work1(i,k,1)
799	if(prevp(i,k).lt.0.) then
800	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
801	prevp(i,k) = max(prevp(i,k),satdt/2)
802	else
803	prevp(i,k) = min(prevp(i,k),satdt/2)
804	endif
805	endif
806	enddo
807	enddo
808	!
809	!===============================================================
810	!
811	! cold rain processes
812	!
813	! - follows the revised ice microphysics processes in HDC
814	! - the processes same as in RH83 and RH84 and LFO behave
815	! following ice crystal hapits defined in HDC, inclduing
816	! intercept parameter for snow (n0s), ice crystal number
817	! concentration (ni), ice nuclei number concentration
818	! (n0i), ice diameter (d)
819	!
820	!===============================================================
821	!
822	rdtcld = 1./dtcld
823	do k = kts, kte
824	do i = its, ite
825	supcol = t0c-t(i,k)
826	supsat = max(q(i,k),qmin)-qs(i,k,2)
827	satdt = supsat/dtcld
828	ifsat = 0
829	!-------------------------------------------------------------
830	! Ni: ice crystal number concentraiton [HDC 5c]
831	!-------------------------------------------------------------
832	! xni(i,k) = min(max(5.38e7*(den(i,k) &
833	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
834	temp = (den(i,k)*max(qci(i,k,2),qmin))
835	temp = sqrt(sqrt(temptemptemp))
836	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
837	eacrs = exp(0.07*(-supcol))
838	!
839	if(supcol.gt.0) then
840	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
841	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
842	diameter = min(dicon * sqrt(xmi),dimax)
843	vt2i = 1.49e4diameter*1.31
844	vt2s = pvtsrslopeb(i,k,2)denfac(i,k)
845	!-------------------------------------------------------------
846	! psaci: Accretion of cloud ice by rain [HDC 10]
847	! (T<T0: I->S)
848	!-------------------------------------------------------------
849	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
850	+diameter*2rslope(i,k,2)
851	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
852	abs(vt2s-vt2i)acrfac/4.
853	endif
854	!-------------------------------------------------------------
855	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
856	! (T<T0: C->S, and T>=T0: C->R)
857	!-------------------------------------------------------------
858	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
859	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2) &
860	rslopeb(i,k,2)qci(i,k,1)*denfac(i,k) &
861	! ,qci(i,k,1)/dtcld)
862	,qci(i,k,1)*rdtcld)
863	endif
864	!-------------------------------------------------------------
865	! pidep: Deposition/Sublimation rate of ice [HDC 9]
866	! (T<T0: V->I or I->V)
867	!-------------------------------------------------------------
868	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
869	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
870	diameter = dicon * sqrt(xmi)
871	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
872	supice = satdt-prevp(i,k)
873	if(pidep(i,k).lt.0.) then
874	! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
875	! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
876	pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
877	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
878	else
879	! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
880	pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
881	endif
882	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
883	endif
884	endif
885	!-------------------------------------------------------------
886	! psdep: deposition/sublimation rate of snow [HDC 14]
887	! (V->S or S->V)
888	!-------------------------------------------------------------
889	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
890	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
891	psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
892	(precs1rslope2(i,k,2)+precs2 &
893	work2(i,k)coeres)/work1(i,k,2)
894	supice = satdt-prevp(i,k)-pidep(i,k)
895	if(psdep(i,k).lt.0.) then
896	! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
897	! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
898	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
899	psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
900	else
901	! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
902	psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
903	endif
904	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
905	ifsat = 1
906	endif
907	!-------------------------------------------------------------
908	! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
909	! (T<T0: V->I)
910	!-------------------------------------------------------------
911	if(supcol.gt.0) then
912	if(supsat.gt.0.and.ifsat.ne.1) then
913	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
914	xni0 = 1.e3exp(0.1supcol)
915	roqi0 = 4.92e-11exp(log(xni0)(1.33))
916	pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
917	! /dtcld)
918	*rdtcld)
919	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
920	endif
921	!
922	!-------------------------------------------------------------
923	! psaut: conversion(aggregation) of ice to snow [HDC 12]
924	! (T<T0: I->S)
925	!-------------------------------------------------------------
926	if(qci(i,k,2).gt.0.) then
927	qimax = roqimax/den(i,k)
928	! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
929	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
930	endif
931	endif
932	!-------------------------------------------------------------
933	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
934	! (T>T0: S->V)
935	!-------------------------------------------------------------
936	if(supcol.lt.0.) then
937	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
938	psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
939	! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
940	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
941	endif
942	enddo
943	enddo
944	!
945	!
946	!----------------------------------------------------------------
947	! check mass conservation of generation terms and feedback to the
948	! large scale
949	!
950	do k = kts, kte
951	do i = its, ite
952	if(t(i,k).le.t0c) then
953	!
954	! cloud water
955	!
956	value = max(qmin,qci(i,k,1))
957	source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
958	if (source.gt.value) then
959	factor = value/source
960	praut(i,k) = praut(i,k)*factor
961	pracw(i,k) = pracw(i,k)*factor
962	psacw(i,k) = psacw(i,k)*factor
963	endif
964	!
965	! cloud ice
966	!
967	value = max(qmin,qci(i,k,2))
968	source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
969	if (source.gt.value) then
970	factor = value/source
971	psaut(i,k) = psaut(i,k)*factor
972	psaci(i,k) = psaci(i,k)*factor
973	pigen(i,k) = pigen(i,k)*factor
974	pidep(i,k) = pidep(i,k)*factor
975	endif
976	!
977	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
978	! update
979	q(i,k) = q(i,k)+work2(i,k)*dtcld
980	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
981	+psacw(i,k))*dtcld,0.)
982	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
983	+prevp(i,k))*dtcld,0.)
984	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
985	-pigen(i,k)-pidep(i,k))*dtcld,0.)
986	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
987	+psaci(i,k)+psacw(i,k))*dtcld,0.)
988	xlf = xls-xl(i,k)
989	xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
990	-xl(i,k)prevp(i,k)-xlfpsacw(i,k)
991	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
992	else
993	!
994	! cloud water
995	!
996	value = max(qmin,qci(i,k,1))
997	source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
998	if (source.gt.value) then
999	factor = value/source
1000	praut(i,k) = praut(i,k)*factor
1001	pracw(i,k) = pracw(i,k)*factor
1002	psacw(i,k) = psacw(i,k)*factor
1003	endif
1004	!
1005	! snow
1006	!
1007	value = max(qcrmin,qrs(i,k,2))
1008	source=(-psevp(i,k))*dtcld
1009	if (source.gt.value) then
1010	factor = value/source
1011	psevp(i,k) = psevp(i,k)*factor
1012	endif
1013	work2(i,k)=-(prevp(i,k)+psevp(i,k))
1014	! update
1015	q(i,k) = q(i,k)+work2(i,k)*dtcld
1016	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1017	+psacw(i,k))*dtcld,0.)
1018	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1019	+prevp(i,k) +psacw(i,k))*dtcld,0.)
1020	qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1021	xlf = xls-xl(i,k)
1022	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1023	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1024	endif
1025	enddo
1026	enddo
1027	!
1028	! Inline expansion for fpvs
1029	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1030	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1031	hsub = xls
1032	hvap = xlv0
1033	cvap = cpv
1034	ttp=t0c+0.01
1035	dldt=cvap-cliq
1036	xa=-dldt/rv
1037	xb=xa+hvap/(rv*ttp)
1038	dldti=cvap-cice
1039	xai=-dldti/rv
1040	xbi=xai+hsub/(rv*ttp)
1041	do k = kts, kte
1042	do i = its, ite
1043	tr=ttp/t(i,k)
1044	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
1045	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1046	qs(i,k,1) = max(qs(i,k,1),qmin)
1047	tr=ttp/t(i,k)
1048	if(t(i,k).lt.ttp) then
1049	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
1050	else
1051	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
1052	endif
1053	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1054	qs(i,k,2) = max(qs(i,k,2),qmin)
1055	enddo
1056	enddo
1057	!
1058	!----------------------------------------------------------------
1059	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1060	! if there exists additional water vapor condensated/if
1061	! evaporation of cloud water is not enough to remove subsaturation
1062	!
1063	do k = kts, kte
1064	do i = its, ite
1065	! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1066	work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/ &
1067	(1.+(xl(i,k))(xl(i,k))/(rv(cpm(i,k)))(qs(i,k,1))/((t(i,k))(t(i,k)))))
1068	work2(i,k) = qci(i,k,1)+work1(i,k,1)
1069	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1070	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1071	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1072	q(i,k) = q(i,k)-pcond(i,k)*dtcld
1073	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1074	t(i,k) = t(i,k)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
1075	enddo
1076	enddo
1077	!
1078	!
1079	!----------------------------------------------------------------
1080	! padding for small values
1081	!
1082	do k = kts, kte
1083	do i = its, ite
1084	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1085	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1086	enddo
1087	enddo
1088	enddo ! big loops
1089	END SUBROUTINE wsm52d
1090	! ...................................................................
1091	REAL FUNCTION rgmma(x)
1092	!-------------------------------------------------------------------
1093	IMPLICIT NONE
1094	!-------------------------------------------------------------------
1095	! rgmma function: use infinite product form
1096	REAL :: euler
1097	PARAMETER (euler=0.577215664901532)
1098	REAL :: x, y
1099	INTEGER :: i
1100	if(x.eq.1.)then
1101	rgmma=0.
1102	else
1103	rgmma=xexp(eulerx)
1104	do i=1,10000
1105	y=float(i)
1106	rgmma=rgmma(1.000+x/y)exp(-x/y)
1107	enddo
1108	rgmma=1./rgmma
1109	endif
1110	END FUNCTION rgmma
1111	!
1112	!--------------------------------------------------------------------------
1113	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1114	!--------------------------------------------------------------------------
1115	IMPLICIT NONE
1116	!--------------------------------------------------------------------------
1117	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1118	xai,xbi,ttp,tr
1119	INTEGER ice
1120	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1121	ttp=t0c+0.01
1122	dldt=cvap-cliq
1123	xa=-dldt/rv
1124	xb=xa+hvap/(rv*ttp)
1125	dldti=cvap-cice
1126	xai=-dldti/rv
1127	xbi=xai+hsub/(rv*ttp)
1128	tr=ttp/t
1129	if(t.lt.ttp.and.ice.eq.1) then
1130	fpvs=psatexp(log(tr)(xai))exp(xbi(1.-tr))
1131	else
1132	fpvs=psatexp(log(tr)(xa))exp(xb(1.-tr))
1133	endif
1134	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1135	END FUNCTION fpvs
1136	!-------------------------------------------------------------------
1137	SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
1138	!-------------------------------------------------------------------
1139	IMPLICIT NONE
1140	!-------------------------------------------------------------------
1141	!.... constants which may not be tunable
1142	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
1143	LOGICAL, INTENT(IN) :: allowed_to_read
1144	REAL :: pi
1145	!
1146	pi = 4.*atan(1.)
1147	xlv1 = cl-cpv
1148	!
1149	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
1150	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
1151	!
1152	bvtr1 = 1.+bvtr
1153	bvtr2 = 2.5+.5*bvtr
1154	bvtr3 = 3.+bvtr
1155	bvtr4 = 4.+bvtr
1156	g1pbr = rgmma(bvtr1)
1157	g3pbr = rgmma(bvtr3)
1158	g4pbr = rgmma(bvtr4) ! 17.837825
1159	g5pbro2 = rgmma(bvtr2) ! 1.8273
1160	pvtr = avtr*g4pbr/6.
1161	eacrr = 1.0
1162	pacrr = pin0ravtrg3pbr.25*eacrr
1163	precr1 = 2.pin0r*.78
1164	precr2 = 2.pin0r.31avtr*.5g5pbro2
1165	xm0 = (di0/dicon)**2
1166	xmmax = (dimax/dicon)**2
1167	roqimax = 2.08e22dimax*8
1168	!
1169	bvts1 = 1.+bvts
1170	bvts2 = 2.5+.5*bvts
1171	bvts3 = 3.+bvts
1172	bvts4 = 4.+bvts
1173	g1pbs = rgmma(bvts1) !.8875
1174	g3pbs = rgmma(bvts3)
1175	g4pbs = rgmma(bvts4) ! 12.0786
1176	g5pbso2 = rgmma(bvts2)
1177	pvts = avts*g4pbs/6.
1178	pacrs = pin0savtsg3pbs.25
1179	precs1 = 4.n0s.65
1180	precs2 = 4.n0s.44avts.5g5pbso2
1181	pidn0r = pidenrn0r
1182	pidn0s = pidensn0s
1183	pacrc = pin0savtsg3pbs.25*eacrc
1184	!
1185	rslopermax = 1./lamdarmax
1186	rslopesmax = 1./lamdasmax
1187	rsloperbmax = rslopermax ** bvtr
1188	rslopesbmax = rslopesmax ** bvts
1189	rsloper2max = rslopermax * rslopermax
1190	rslopes2max = rslopesmax * rslopesmax
1191	rsloper3max = rsloper2max * rslopermax
1192	rslopes3max = rslopes2max * rslopesmax
1193	!
1194	END SUBROUTINE wsm5init
1195	END MODULE module_mp_wsm5

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

Download in other formats:

Original Format