Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_wsm3.F @ 3553

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

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

Download in other formats:

Original Format