Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm5.F @ 3431

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

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

Download in other formats:

Original Format