Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_sbu_ylin.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 58.0 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2
3	!--- The code is based on Lin and Colle (A New Bulk Microphysical Scheme
4	! that Includes Riming Intensity and Temperature Dependent Ice Characteristics, 2011, MWR)
5	! and Lin et al. (Parameterization of riming intensity and its impact on ice fall speed using ARM data, 2011, MWR)
6	!--- NOTE: 1) Prognose variables are: qi,PI(precipitating ice, qs, which includes snow, partially rimed snow and graupel),qw,qr
7	!--- 2) Sedimentation flux is based on Prudue Lin scheme
8	!--- 2) PI has varying properties depending on riming intensity (Ri, diagnosed currently following Lin et al. (2011, MWR) and T
9	!--- 3) Autoconverion is based on Liu and Daum (2004)
10	!--- 4) PI size distribution assuming Gamma distribution, but mu_s=0 (Exponential) currently
11	!--- 5) No density dependent fall speed since the V-D is derived using Best number approach, which already includes density effect
12	!--- 6) Future work will include radar equivalent reflectivity using the new PI property (A-D, M-D, N(D)). If you use RIP for reflectivity
13	!--- computation, please note that snow is (1-Ri)qs and graupel is Riqs. Otherwise, reflectivity will be underestimated.
14	!--- 7) The Liu and Daum autoconverion is quite sensitive on Nt_c. For mixed-phase cloud and marine environment, Nt_c of 10 or 20 is suggested.
15	!--- default value is 10E.6. Change accordingly for your use.
16	!--- 8) Eq.7 and 8 are not in SI units and need to be converted in the code. the
17	! paper treats the units in Eq.7 and 8 as cgs, and so need 1e-2^(2-ba) in
18	! the code, and that would give the plots in the paper. However, there is
19	! large uncertainty with this parameter, and one could argue that the units
20	! for these equations could be mm-g-s instead, which would mean 1e-3^(2-ba)
21	! in the code. This increases the snow fallspeed and gives an even
22	! better comparison of aa and ba with obs in paper.
23
24
25	MODULE module_mp_sbu_ylin
26	USE module_wrf_error
27	!
28	!..Parameters user might change based on their need
29	REAL, PARAMETER, PRIVATE :: RH = 1.0
30	REAL, PARAMETER, PRIVATE :: xnor = 8.0e6
31	REAL, PARAMETER, PRIVATE :: Nt_c = 10.E6
32	!..Water vapor and air gas constants at constant pressure
33	REAL, PARAMETER, PRIVATE :: Rvapor = 461.5
34	REAL, PARAMETER, PRIVATE :: oRv = 1./Rvapor
35	REAL, PARAMETER, PRIVATE :: Rair = 287.04
36	REAL, PARAMETER, PRIVATE :: Cp = 1004.0
37	REAL, PARAMETER, PRIVATE :: grav = 9.81
38	REAL, PARAMETER, PRIVATE :: rhowater = 1000.0
39	REAL, PARAMETER, PRIVATE :: rhosnow = 100.0
40
41	REAL, PARAMETER, PRIVATE :: SVP1=0.6112
42	REAL, PARAMETER, PRIVATE :: SVP2=17.67
43	REAL, PARAMETER, PRIVATE :: SVP3=29.65
44	REAL, PARAMETER, PRIVATE :: SVPT0=273.15
45	REAL, PARAMETER, PRIVATE :: EP1=Rvapor/Rair-1.
46	REAL, PARAMETER, PRIVATE :: EP2=Rair/Rvapor
47	!..Enthalpy of sublimation, vaporization, and fusion at 0C.
48	REAL, PARAMETER, PRIVATE :: XLS = 2.834E6
49	REAL, PARAMETER, PRIVATE :: XLV = 2.5E6
50	REAL, PARAMETER, PRIVATE :: XLF = XLS - XLV
51
52	!
53	REAL, PARAMETER, PRIVATE :: &
54	qi0 = 1.0e-3, & !--- ice aggregation to snow threshold
55	xmi50 = 4.8e-10, xmi40 = 2.46e-10, &
56	xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5, &
57	di50 = 1.0e-4, xmi = 4.19e-13, & !--- parameters used in BF process
58	bv_r = 0.8, bv_i = 0.25, &
59	o6 = 1./6., cdrag = 0.6, &
60	avisc = 1.49628e-6, adiffwv = 8.7602e-5, &
61	axka = 1.4132e3, cw = 4.187e3, ci = 2.093e3
62	CONTAINS
63
64	!-------------------------------------------------------------------
65	! Lin et al., 1983, JAM, 1065-1092, and
66	! Rutledge and Hobbs, 1984, JAS, 2949-2972
67	!-------------------------------------------------------------------
68	SUBROUTINE sbu_ylin(th &
69	,qv, ql, qr &
70	,qi, qs, Ri3D &
71	,rho, pii, p &
72	,dt_in &
73	,z,ht, dz8w &
74	,RAINNC, RAINNCV &
75	,ids,ide, jds,jde, kds,kde &
76	,ims,ime, jms,jme, kms,kme &
77	,its,ite, jts,jte, kts,kte &
78	)
79	!-------------------------------------------------------------------
80	IMPLICIT NONE
81	!-------------------------------------------------------------------
82	!
83	!
84	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
85	ims,ime, jms,jme, kms,kme , &
86	its,ite, jts,jte, kts,kte
87
88	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
89	INTENT(INOUT) :: &
90	th, &
91	qv, &
92	qi,ql, &
93	qs,qr
94
95	! YLIN
96	! Adding RI3D as a variable to the interface
97	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
98	INTENT(INOUT) :: Ri3D
99
100
101	!
102	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
103	INTENT(IN ) :: &
104	rho, &
105	pii, &
106	z,p, &
107	dz8w
108
109
110	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
111	REAL, INTENT(IN ) :: dt_in
112
113	REAL, DIMENSION( ims:ime , jms:jme ), &
114	INTENT(INOUT) :: RAINNC, &
115	RAINNCV
116
117	! LOCAL VAR
118
119	INTEGER :: min_q, max_q
120
121	REAL, DIMENSION( its:ite , jts:jte ) &
122	:: rain, snow,ice
123
124	REAL, DIMENSION( kts:kte ) :: qvz, qlz, qrz, &
125	qiz, qsz, qgz, &
126	thz, &
127	tothz, rhoz, &
128	orhoz, sqrhoz, &
129	prez, zz, &
130	dzw
131
132	! Added vertical profile of Ri (riz) as a variable
133	REAL, DIMENSION( kts:kte ) :: riz
134
135	!
136	REAL :: dt, pptice, pptrain, pptsnow, pptgraul, rhoe_s
137
138	INTEGER :: i,j,k
139
140	dt=dt_in
141	rhoe_s=1.29
142
143	j_loop: DO j = jts, jte
144	i_loop: DO i = its, ite
145	!
146	!- write data from 3-D to 1-D
147	!
148	DO k = kts, kte
149	qvz(k)=qv(i,k,j)
150	qlz(k)=ql(i,k,j)
151	qrz(k)=qr(i,k,j)
152	qiz(k)=qi(i,k,j)
153	qsz(k)=qs(i,k,j)
154	thz(k)=th(i,k,j)
155	rhoz(k)=rho(i,k,j)
156	orhoz(k)=1./rhoz(k)
157	prez(k)=p(i,k,j)
158	! sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
159	! no density dependence of fall speed as Note #5, you can turn it on to increase fall speed at low pressure.
160	sqrhoz(k)=1.0
161	tothz(k)=pii(i,k,j)
162	zz(k)=z(i,k,j)
163	dzw(k)=dz8w(i,k,j)
164	END DO
165	!
166	pptrain=0.
167	pptsnow=0.
168	pptice =0.
169
170	! CALL wrf_debug ( 100 , 'microphysics_driver: calling clphy1d_ylin' )
171
172	CALL clphy1d_ylin( dt, qvz, qlz, qrz, qiz, qsz, &
173	thz, tothz, rhoz, orhoz, sqrhoz, &
174	prez, zz, dzw, ht(I,J), &
175	pptrain, pptsnow, pptice, &
176	kts, kte, i, j, riz )
177
178	!
179	! Precipitation from cloud microphysics -- only for one time step
180	!
181	! unit is transferred from m to mm
182	!
183	rain(i,j)= pptrain
184	snow(i,j)= pptsnow
185	ice(i,j) = pptice
186	!
187	RAINNCV(i,j)= pptrain + pptsnow + pptice
188	RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptice
189
190	!
191	!- update data from 1-D back to 3-D
192	!
193	DO k = kts, kte
194	qv(i,k,j)=qvz(k)
195	ql(i,k,j)=qlz(k)
196	qr(i,k,j)=qrz(k)
197	th(i,k,j)=thz(k)
198	qi(i,k,j)=qiz(k)
199	qs(i,k,j)=qsz(k)
200	ri3d(i,k,j)=riz(k)
201	END DO
202	!
203	ENDDO i_loop
204	ENDDO j_loop
205
206	END SUBROUTINE sbu_ylin
207
208
209	!-----------------------------------------------------------------------
210	SUBROUTINE clphy1d_ylin(dt, qvz, qlz, qrz, qiz, qsz, &
211	thz, tothz, rho, orho, sqrho, &
212	prez, zz, dzw, zsfc, pptrain, pptsnow,pptice, &
213	kts, kte, i, j,riz )
214	!-----------------------------------------------------------------------
215	IMPLICIT NONE
216	!-----------------------------------------------------------------------
217	! This program handles the vertical 1-D cloud micphysics
218	!-----------------------------------------------------------------------
219	! avisc: constant in empirical formular for dynamic viscosity of air
220	! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
221	! adiffwv: constant in empirical formular for diffusivity of water
222	! vapor in air
223	! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
224	! axka: constant in empirical formular for thermal conductivity of air
225	! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
226	! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
227	! xmi50: mass of a 50 micron ice crystal
228	! = 4.8e-10 [kg] =4.8e-7 [g]
229	! xmi40: mass of a 40 micron ice crystal
230	! = 2.46e-10 [kg] = 2.46e-7 [g]
231	! di50: diameter of a 50 micro (radius) ice crystal
232	! =1.0e-4 [m]
233	! xmi: mass of one cloud ice crystal
234	! =4.19e-13 [kg] = 4.19e-10 [g]
235	! oxmi=1.0/xmi
236	!
237	! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
238	! Hsie et al.(1980) and Rutledge and Hobbs(1983) )
239	! bni=0.5 [K-1]
240	! xmnin: mass of a natural ice nucleus
241	! = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
242	! Hsie et al. (1980)
243	! = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
244
245	! av_r: av_r in empirical formular for terminal
246	! velocity of raindrop
247	! =2115.0 [cm*(1-b)/s] = 2115.00.01(1-b) [m(1-b)/s]
248	! bv_r: bv_r in empirical formular for terminal
249	! velocity of raindrop
250	! =0.8
251	! av_i: av_i in empirical formular for terminal
252	! velocity of snow
253	! =152.93 [cm*(1-d)/s] = 152.930.01(1-d) [m(1-d)/s]
254	! bv_i: bv_i in empirical formular for terminal
255	! velocity of snow
256	! =0.25
257	! vf1r: ventilation factors for rain =0.78
258	! vf2r: ventilation factors for rain =0.31
259	! vf1s: ventilation factors for snow =0.65
260	! vf2s: ventilation factors for snow =0.44
261	!
262	!----------------------------------------------------------------------
263
264	INTEGER, INTENT(IN ) :: kts, kte, i, j
265
266	REAL, DIMENSION( kts:kte ), &
267	INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, &
268	thz
269
270	REAL, DIMENSION( kts:kte ), &
271	INTENT(IN ) :: tothz, rho, orho, sqrho, &
272	prez, zz, dzw
273
274
275	REAL, INTENT(INOUT) :: pptrain, pptsnow, pptice
276
277	REAL, INTENT(IN ) :: dt, zsfc
278
279	! local vars
280
281	REAL :: obp4, bp3, bp5, bp6, odp4, &
282	dp3, dp5, dp5o2
283
284	! temperary vars
285
286	REAL :: tmp, tmp0, tmp1, tmp2,tmp3, &
287	tmp4, tmpa,tmpb,tmpc,tmpd,alpha1, &
288	qic, abi,abr, abg, odtberg, &
289	vti50,eiw,eri,esi,esr, esw, &
290	erw,delrs,term0,term1, &
291	Ap, Bp, &
292	factor, tmp_r, tmp_s,tmp_g, &
293	qlpqi, rsat, a1, a2, xnin
294
295	!
296	REAL, DIMENSION( kts:kte ) :: oprez, tem, temcc, theiz, qswz, &
297	qsiz, qvoqswz, qvoqsiz, qvzodt, &
298	qlzodt, qizodt, qszodt, qrzodt
299
300	!--- microphysical processes
301
302	REAL, DIMENSION( kts:kte ) :: psnow, psaut, psfw, psfi, praci, &
303	piacr, psaci, psacw, psdep, pssub, &
304	pracs, psacr, psmlt, psmltevp, &
305	prain, praut, pracw, prevp, pvapor, &
306	pclw, pladj, pcli, pimlt, pihom, &
307	pidw, piadj, pgfr, &
308	qschg
309	!
310
311	REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, &
312	schmidt, xka
313	!---- new snow parameters
314
315	REAL, DIMENSION( kts:kte ):: ab_s,ab_r,ab_riming,lamc
316	REAL, DIMENSION( kts:kte ):: cap_s !---- capacitance of snow
317
318	REAL, PARAMETER :: vf1s = 0.65, vf2s = 0.44, vf1r =0.78 , vf2r = 0.31
319
320	REAL, PARAMETER :: am_c1=0.004, am_c2= 6e-5, am_c3=0.15
321	REAL, PARAMETER :: bm_c1=1.85, bm_c2= 0.003, bm_c3=1.25
322	REAL, PARAMETER :: aa_c1=1.28, aa_c2= -0.012, aa_c3=-0.6
323	REAL, PARAMETER :: ba_c1=1.5, ba_c2= 0.0075, ba_c3=0.5
324
325	REAL, PARAMETER :: best_a=1.08 , best_b = 0.499
326	REAL, DIMENSION(kts:kte):: am_s,bm_s,av_s,bv_s,Ri,N0_s,tmp_ss,lams
327	REAL, DIMENSION(kts:kte):: aa_s,ba_s,tmp_sa
328	REAL, PARAMETER :: mu_s=0.,mu_i=0.,mu_r=0.
329
330	REAL :: tc0, disp, Dc_liu, eta, mu_c, R6c !--- for Liu's autoconversion
331
332	! Adding variable Riz, which will duplicate Ri but be a copy passed upward
333
334	REAL, DIMENSION(kts:kte) :: Riz
335
336	REAL, DIMENSION( kts:kte ) :: vtr, vts, &
337	vtrold, vtsold, vtiold, &
338	xlambdar, xlambdas, &
339	olambdar, olambdas
340
341	REAL :: episp0k, dtb, odtb, pi, pio4, &
342	pio6, oxLf, xLvocp, xLfocp, av_r, &
343	av_i, ocdrag, gambp4, gamdp4, &
344	gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
345	gambp6, gam3pt5, gam2pt75, gambp5o2,&
346	gamdp5o2, cwoxlf, ocp, xni50, es
347	!
348	REAL :: qvmin=1.e-20
349	REAL :: temc1,save1,save2,xni50mx
350
351	! for terminal velocity flux
352
353	INTEGER :: min_q, max_q, max_ri_k, k
354	REAL :: max_ri
355	REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
356	LOGICAL :: notlast
357	!
358	mu_c = AMIN1(15., (1000.E6/Nt_c + 2.))
359	R6c = 10.0E-6 !---- 10 micron, threshold radius of cloud droplet
360	dtb=dt
361	odtb=1./dtb
362	pi =acos(-1.)
363	pio4=acos(-1.)/4.
364	pio6=acos(-1.)/6.
365	ocp=1./cp
366	oxLf=1./xLf
367	xLvocp=xLv/cp
368	xLfocp=xLf/cp
369	Cpor=cp/Rair
370	oxmi=1.0/xmi
371	cwoxlf=cw/xlf
372	av_r=2115.00.01*(1-bv_r)
373	av_i=152.930.01*(1-bv_i)
374	ocdrag=1./Cdrag
375	episp0k=RHep21000.*svp1
376	!
377	gambp4=ggamma(bv_r+4.)
378	gamdp4=ggamma(bv_i+4.)
379	gambp3=ggamma(bv_r+3.)
380	gambp6=ggamma(bv_r+6)
381	gambp5o2=ggamma((bv_r+5.)/2.)
382	gamdp5o2=ggamma((bv_i+5.)/2.)
383	!
384	! oprez 1./prez ( prez : pressure)
385	! qsw saturated mixing ratio on water surface
386	! qsi saturated mixing ratio on ice surface
387	! episp0k RHesaturated pressure at 273.15 K = 611.2 hPa (Rogers and Yau 1989)
388	! qvoqsw qv/qsw
389	! qvoqsi qv/qsi
390	! qvzodt qv/dt
391	! qlzodt ql/dt
392	! qizodt qi/dt
393	! qszodt qs/dt
394	! qrzodt qr/dt
395	! temcc temperature in dregee C
396	!
397
398	obp4=1.0/(bv_r+4.0)
399	bp3=bv_r+3.0
400	bp5=bv_r+5.0
401	bp6=bv_r+6.0
402	odp4=1.0/(bv_i+4.0)
403	dp3=bv_i+3.0
404	dp5=bv_i+5.0
405	dp5o2=0.5*(bv_i+5.0)
406	!
407	do k=kts,kte
408	oprez(k)=1./prez(k)
409	qlz(k)=amax1( 0.0,qlz(k) )
410	qiz(k)=amax1( 0.0,qiz(k) )
411	qvz(k)=amax1( qvmin,qvz(k) )
412	qsz(k)=amax1( 0.0,qsz(k) )
413	qrz(k)=amax1( 0.0,qrz(k) )
414	tem(k)=thz(k)*tothz(k)
415	temcc(k)=tem(k)-273.15
416	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) ) !--- RY89 Eq(2.17)
417	qswz(k)=ep2*es/(prez(k)-es)
418	if (tem(k) .lt. 273.15 ) then
419	es=1000.svp1exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
420	qsiz(k)=ep2*es/(prez(k)-es)
421	if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
422	else
423	qsiz(k)=qswz(k)
424	endif
425	!
426	qvoqswz(k)=qvz(k)/qswz(k)
427	qvoqsiz(k)=qvz(k)/qsiz(k)
428	qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
429	qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
430	qizodt(k)=amax1( 0.0,odtb*qiz(k) )
431	qszodt(k)=amax1( 0.0,odtb*qsz(k) )
432	qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
433	theiz(k)=thz(k)+(xlvocpqvz(k)-xlfocpqiz(k))/tothz(k)
434	enddo
435
436	do k=kts,kte
437
438	psnow(k)=0.0
439	psaut(k)=0.0
440	psfw(k)=0.0
441	psfi(k)=0.0
442	praci(k)=0.0
443	piacr(k)=0.0
444	psaci(k)=0.0
445	psacw(k)=0.0
446	psdep(k)=0.0
447	pssub(k)=0.0
448	pracs(k)=0.0
449	psacr(k)=0.0
450	psmlt(k)=0.0
451	psmltevp(k)=0.0
452
453	prain(k)=0.0
454	praut(k)=0.0
455	pracw(k)=0.0
456	prevp(k)=0.0
457	pgfr(k)=0.0
458
459	pvapor(k)=0.0
460
461	pclw(k)=0.0
462	pladj(k)=0.0
463
464	pcli(k)=0.0
465	pimlt(k)=0.0
466	pihom(k)=0.0
467	pidw(k)=0.0
468	piadj(k)=0.0
469
470	qschg(k)=0.
471
472	enddo
473
474	!***********************************************************************
475	!*** compute viscosity,difusivity,thermal conductivity, and ****
476	!*** Schmidt number ****
477	!***********************************************************************
478	!c------------------------------------------------------------------
479	!c viscmu: dynamic viscosity of air kg/m/s
480	!c visc: kinematic viscosity of air = viscmu/rho (m2/s)
481	!c avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
482	!c viscmu=1.718e-5 kg/m/s in RH
483	!c diffwv: Diffusivity of water vapor in air
484	!c adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
485	!c = 8.7602 (8.794 in MM5) gcm/s3
486	!c diffwv(k)=2.26e-5 m2/s
487	!c schmidt: Schmidt number=visc/diffwv
488	!c xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
489	!c xka(k)=2.43e-2 J/m/s/K in RH
490	!c axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
491	!c------------------------------------------------------------------
492	DO k=kts,kte
493	viscmu(k)=avisctem(k)*1.5/(tem(k)+120.0)
494	visc(k)=viscmu(k)*orho(k)
495	diffwv(k)=adiffwvtem(k)1.81oprez(k)
496	schmidt(k)=visc(k)/diffwv(k)
497	xka(k)=axka*viscmu(k)
498	rs0(k)=ep21000.svp1/(prez(k)-1000.*svp1)
499	END DO
500	!
501	! ---- YLIN, set snow variables
502	!
503	!---- A+B in depositional growth, the first try just take from Rogers and Yau(1989)
504	! ab_s(k) = lsublsuborv/(tcond(k)*temp(k))+&
505	! rvtemp(k)/(diffu(k)qvsi(k))
506
507	do k = kts, kte
508	tc0 = tem(k)-273.15
509	if (rho(k)qlz(k) .gt. 1e-5 .AND. rho(k)qsz(k) .gt. 1e-5) then
510	Ri(k) = 1.0/(1.0+6e-5/(rho(k)*1.170qlz(k)qsz(k)*0.170))
511	else
512	Ri(k) = 0
513	endif
514	enddo
515	!
516	!--- make sure Ri does not decrease downward in a column
517	!
518	max_ri_k = MAXLOC(Ri,dim=1)
519	max_ri = MAXVAL(Ri)
520
521	do k = kts, max_ri_k
522	Ri(k) = max_ri
523	enddo
524
525	!--- YLIN, get PI properties
526	do k = kts, kte
527	Ri(k) = AMAX1(0.,AMIN1(Ri(k),1.0))
528	! Store the value of Ri(k) as Riz(k)
529	Riz(k) = Ri(k)
530
531	cap_s(k)= 0.25*(1+Ri(k))
532	tc0 = AMIN1(-0.1, tem(k)-273.15)
533	N0_s(k) = amin1(2.0E8, 2.0E6exp(-0.12tc0))
534	am_s(k) = am_c1+am_c2tc0+am_c3Ri(k)*Ri(k) !--- Heymsfield 2007
535	am_s(k) = AMAX1(0.000023,am_s(k)) !--- use the a_min in table 1 of Heymsfield
536	bm_s(k) = bm_c1+bm_c2tc0+bm_c3Ri(k)
537	bm_s(k) = AMIN1(bm_s(k),3.0) !---- capped by 3
538	!--- converting from cgs to SI unit
539	am_s(k) = 10*(2bm_s(k)-3.0)*am_s(k)
540	aa_s(k) = aa_c1 + aa_c2tc0 + aa_c3Ri(k)
541	ba_s(k) = ba_c1 + ba_c2tc0 + ba_c3Ri(k)
542	!--- convert from mm.g.s to SI unit (this will give larger PI fall speed than in the paper)
543	aa_s(k) = (1e-3)*(2.0-ba_s(k))aa_s(k)
544	!---- get v from Mitchell 1996
545	av_s(k) = best_aviscmu(k)(2gravam_s(k)/rho(k)/aa_s(k)/(viscmu(k)2))best_b
546	bv_s(k) = best_b*(bm_s(k)-ba_s(k)+2)-1
547
548	tmp_ss(k)= bm_s(k)+mu_s+1
549	tmp_sa(k)= ba_s(k)+mu_s+1
550
551	enddo
552
553	!
554	!***********************************************************************
555	! Calculate precipitation fluxes due to terminal velocities.
556	!***********************************************************************
557	!
558	!- Calculate termianl velocity (vt?) of precipitation q?z
559	!- Find maximum vt? to determine the small delta t
560	!
561	!-- rain
562	!
563	! CALL wrf_debug ( 100 , 'module_ylin, start precip fluxes' )
564
565	t_del_tv=0.
566	del_tv=dtb
567	notlast=.true.
568
569	DO while (notlast)
570	!
571	min_q=kte
572	max_q=kts-1
573	!
574	do k=kts,kte-1
575	if (qrz(k) .gt. 1.0e-8) then
576	min_q=min0(min_q,k)
577	max_q=max0(max_q,k)
578	tmp1=sqrt(pirhowaterxnor/rho(k)/qrz(k))
579	tmp1=sqrt(tmp1)
580	vtrold(k)=o6av_rgambp4sqrho(k)/tmp1*bv_r
581	if (k .eq. 1) then
582	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
583	else
584	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
585	endif
586	else
587	vtrold(k)=0.
588	endif
589	enddo
590
591	if (max_q .ge. min_q) then
592	!
593	!- Check if the summation of the small delta t >= big delta t
594	! (t_del_tv) (del_tv) (dtb)
595
596	t_del_tv=t_del_tv+del_tv
597	!
598	if ( t_del_tv .ge. dtb ) then
599	notlast=.false.
600	del_tv=dtb+del_tv-t_del_tv
601	endif
602	!
603	fluxin=0.
604	do k=max_q,min_q,-1
605	fluxout=rho(k)vtrold(k)qrz(k)
606	flux=(fluxin-fluxout)/rho(k)/dzw(k)
607	tmpqrz=qrz(k)
608	qrz(k)=qrz(k)+del_tv*flux
609	fluxin=fluxout
610	enddo
611	if (min_q .eq. 1) then
612	pptrain=pptrain+fluxin*del_tv
613	else
614	qrz(min_q-1)=qrz(min_q-1)+del_tv* &
615	fluxin/rho(min_q-1)/dzw(min_q-1)
616	endif
617	!
618	else
619	notlast=.false.
620	endif
621	ENDDO
622
623	!
624	!-- snow
625	!
626	t_del_tv=0.
627	del_tv=dtb
628	notlast=.true.
629
630	DO while (notlast)
631	!
632	min_q=kte
633	max_q=kts-1
634	!
635	do k=kts,kte-1
636	if (qsz(k) .gt. 1.0e-8) then
637	min_q=min0(min_q,k)
638	max_q=max0(max_q,k)
639
640	tmp1= (am_s(k)N0_s(k)ggamma(tmp_ss(k))*orho(k)/qsz(k))&
641	**(1./tmp_ss(k))
642
643	vtsold(k)= sqrho(k)av_s(k)ggamma(bv_s(k)+tmp_ss(k))/ &
644	ggamma(tmp_ss(k))/(tmp1**bv_s(k))
645
646	if (k .eq. 1) then
647	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
648	else
649	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
650	endif
651	else
652	vtsold(k)=0.
653	endif
654	enddo
655
656	if (max_q .ge. min_q) then
657	!
658	!
659	!- Check if the summation of the small delta t >= big delta t
660	! (t_del_tv) (del_tv) (dtb)
661
662	t_del_tv=t_del_tv+del_tv
663
664	if ( t_del_tv .ge. dtb ) then
665	notlast=.false.
666	del_tv=dtb+del_tv-t_del_tv
667	endif
668	!
669	fluxin=0.
670	do k=max_q,min_q,-1
671	fluxout=rho(k)vtsold(k)qsz(k)
672	flux=(fluxin-fluxout)/rho(k)/dzw(k)
673	qsz(k)=qsz(k)+del_tv*flux
674	qsz(k)=amax1(0.,qsz(k))
675	fluxin=fluxout
676	enddo
677	if (min_q .eq. 1) then
678	pptsnow=pptsnow+fluxin*del_tv
679	else
680	qsz(min_q-1)=qsz(min_q-1)+del_tv* &
681	fluxin/rho(min_q-1)/dzw(min_q-1)
682	endif
683	!
684	else
685	notlast=.false.
686	endif
687
688	ENDDO
689
690	!
691	!-- cloud ice (03/21/02) using Heymsfield and Donner (1990) Vi=3.29*qi^0.16
692	!
693	t_del_tv=0.
694	del_tv=dtb
695	notlast=.true.
696	!
697	DO while (notlast)
698	!
699	min_q=kte
700	max_q=kts-1
701	!
702	do k=kts,kte-1
703	if (qiz(k) .gt. 1.0e-8) then
704	min_q=min0(min_q,k)
705	max_q=max0(max_q,k)
706	vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16 ! Heymsfield and Donner
707	if (k .eq. 1) then
708	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
709	else
710	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
711	endif
712	else
713	vtiold(k)=0.
714	endif
715	enddo
716
717	if (max_q .ge. min_q) then
718	!
719	!- Check if the summation of the small delta t >= big delta t
720	! (t_del_tv) (del_tv) (dtb)
721
722	t_del_tv=t_del_tv+del_tv
723
724	if ( t_del_tv .ge. dtb ) then
725	notlast=.false.
726	del_tv=dtb+del_tv-t_del_tv
727	endif
728
729	fluxin=0.
730	do k=max_q,min_q,-1
731	fluxout=rho(k)vtiold(k)qiz(k)
732	flux=(fluxin-fluxout)/rho(k)/dzw(k)
733	qiz(k)=qiz(k)+del_tv*flux
734	qiz(k)=amax1(0.,qiz(k))
735	fluxin=fluxout
736	enddo
737	if (min_q .eq. 1) then
738	pptice=pptice+fluxin*del_tv
739	else
740	qiz(min_q-1)=qiz(min_q-1)+del_tv* &
741	fluxin/rho(min_q-1)/dzw(min_q-1)
742	endif
743	!
744	else
745	notlast=.false.
746	endif
747	!
748	ENDDO
749
750	! CALL wrf_debug ( 100 , 'module_ylin: end precip flux' )
751
752	! Microphpysics processes
753
754	DO 2000 k=kts,kte
755	!
756	qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
757	qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
758	qizodt(k)=amax1( 0.0,odtb*qiz(k) )
759	qszodt(k)=amax1( 0.0,odtb*qsz(k) )
760	qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
761
762	!***********************************************************************
763	!*** diagnose mixing ratios (qrz,qsz), terminal ***
764	!*** velocities (vtr,vts), and slope parameters in size ***
765	!*** distribution(xlambdar,xlambdas) of rain and snow ***
766	!*** follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7) ***
767	!***********************************************************************
768	!
769	!**** assuming no cloud water can exist in the top two levels due to
770	!**** radiation consideration
771	!
772	!! if
773	!! unsaturated,
774	!! no cloud water, rain, ice, snow
775	!! then
776	!! skip these processes and jump to line 2000
777	!
778	!
779	tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)
780	if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k) &
781	.and. tmp .eq. 0.0 ) go to 2000
782	!
783	!! calculate terminal velocity of rain
784	!
785	if (qrz(k) .gt. 1.0e-8) then
786	tmp1=sqrt(pirhowaterxnor*orho(k)/qrz(k))
787	xlambdar(k)=sqrt(tmp1)
788	olambdar(k)=1.0/xlambdar(k)
789	vtrold(k)=o6av_rgambp4sqrho(k)olambdar(k)**bv_r
790	else
791	vtrold(k)=0.
792	olambdar(k)=0.
793	endif
794	!
795	if (qrz(k) .gt. 1.0e-8) then
796	tmp1=sqrt(pirhowaterxnor*orho(k)/qrz(k))
797	xlambdar(k)=sqrt(tmp1)
798	olambdar(k)=1.0/xlambdar(k)
799	vtr(k)=o6av_rgambp4sqrho(k)olambdar(k)**bv_r
800	else
801	vtr(k)=0.
802	olambdar(k)=0.
803	endif
804	!
805	!! calculate terminal velocity of snow
806	!
807	if (qsz(k) .gt. 1.0e-8) then
808	tmp1= (am_s(k)N0_s(k)ggamma(tmp_ss(k))*orho(k)/qsz(k))&
809	**(1./tmp_ss(k))
810	olambdas(k)=1.0/tmp1
811	vtsold(k)= sqrho(k)av_s(k)ggamma(bv_s(k)+tmp_ss(k))/ &
812	ggamma(tmp_ss(k))/(tmp1**bv_s(k))
813
814	else
815	vtsold(k)=0.
816	olambdas(k)=0.
817	endif
818	!
819	if (qsz(k) .gt. 1.0e-8) then
820	tmp1= (am_s(k)N0_s(k)ggamma(tmp_ss(k))*orho(k)/qsz(k))&
821	**(1./tmp_ss(k))
822	olambdas(k)=1.0/tmp1
823	vts(k)= sqrho(k)av_s(k)ggamma(bv_s(k)+tmp_ss(k))/ &
824	ggamma(tmp_ss(k))/(tmp1**bv_s(k))
825
826	else
827	vts(k)=0.
828	olambdas(k)=0.
829	endif
830
831	!---------- start of snow/ice processes below freezing
832
833	if (tem(k) .lt. 273.15) then
834
835	!
836	!***********************************************************************
837	!******* snow production processes for T < 0 C ********
838	!***********************************************************************
839	!c
840	!c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
841	!c! psaut=alpha1*(qi-qi0)
842	!c! alpha1=1.0e-3exp(0.025(T-T0))
843	!c
844	alpha1=1.0e-3exp( 0.025temcc(k) )
845	!
846	if(temcc(k) .lt. -20.0) then
847	tmp1=-7.6+4.0exp( -0.2443e-3(abs(temcc(k))-20)**2.455 )
848	qic=1.0e-3exp(tmp1)orho(k)
849	else
850	qic=qi0
851	end if
852
853	tmp1=odtb(qiz(k)-qic)(1.0-exp(-alpha1*dtb))
854	psaut(k)=amax1( 0.0,tmp1 )
855
856	!c
857	!c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
858	!c this process only considered when -31 C < T < 0 C
859	!c Lin (33) and Hsie (17)
860	!c
861	!c!
862	!c! parama1 and parama2 functions must be user supplied
863	!c!
864
865	if( qlz(k) .gt. 1.0e-10 ) then
866	temc1=amax1(-30.99,temcc(k))
867	a1=parama1( temc1 )
868	a2=parama2( temc1 )
869	tmp1=1.0-a2
870	!! change unit from cgs to mks
871	a1=a10.001*tmp1
872	!! dtberg is the time needed for a crystal to grow from 40 to 50 um
873	!! odtberg=1.0/dtberg
874	odtberg=(a1tmp1)/(xmi50tmp1-xmi40*tmp1)
875	!
876	!! compute terminal velocity of a 50 micron ice cystal
877	!
878	vti50=av_idi50bv_isqrho(k)
879	!
880	eiw=1.0
881	save1=a1xmi50*a2
882	save2=0.25pieiwrho(k)di50di50vti50
883	!
884	tmp2=( save1 + save2*qlz(k) )
885	!
886	!! maximum number of 50 micron crystals limited by the amount
887	!! of supercool water
888	!
889	xni50mx=qlzodt(k)/tmp2
890	!
891	!! number of 50 micron crystals produced
892	!
893	xni50=qiz(k)( 1.0-exp(-dtbodtberg) )/xmi50
894	xni50=amin1(xni50,xni50mx)
895	!
896	tmp3=odtbtmp2/save2( 1.0-exp(-save2xni50dtb) )
897	psfw(k)=amin1( tmp3,qlzodt(k) )
898	!c
899	!c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
900	!c this process only considered when -31 C < T < 0 C
901	!c
902	tmp1=xni50*xmi50-psfw(k)
903	psfi(k)=amin1(tmp1,qizodt(k))
904	end if
905	!
906	!
907	if(qrz(k) .le. 0.0) go to 1000
908	!
909	! Processes (4) and (5) only need when qrz > 0.0
910	!
911	!c
912	!c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
913	!c produce PI
914	!c
915	eri=1.0
916	save1=pio4erixnorav_rsqrho(k)
917	tmp1=save1gambp3olambdar(k)**bp3
918	praci(k)=qizodt(k)( 1.0-exp(-tmp1dtb) )
919
920	!c
921	!c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
922	!c
923	tmp2=qiz(k)save1rho(k)pio6rhowatergambp6oxmi* &
924	olambdar(k)**bp6
925	piacr(k)=amin1( tmp2,qrzodt(k) )
926
927	!
928	1000 continue
929	!
930	if(qsz(k) .le. 0.0) go to 1200
931	!
932	! Compute the following processes only when qsz > 0.0
933	!
934	!c
935	!c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
936	!c
937	esi=exp( 0.025*temcc(k) )
938	save1 = aa_s(k)sqrho(k)N0_s(k)* &
939	ggamma(bv_s(k)+tmp_sa(k))olambdas(k)*(bv_s(k)+tmp_sa(k))
940
941	tmp1=esi*save1
942	psaci(k)=qizodt(k)( 1.0-exp(-tmp1dtb) )
943
944	!c
945	!c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
946	!c
947	esw=1.0
948	tmp1=esw*save1
949	psacw(k)=qlzodt(K)( 1.0-exp(-tmp1dtb) )
950
951	!c
952	!c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
953	!c includes consideration of ventilation effect
954	!c
955	tmpa=rvaporxka(k)tem(k)*tem(k)
956	tmpb=xlsxlsrho(k)qsiz(k)diffwv(k)
957	tmpc=tmpaqsiz(k)diffwv(k)
958	abi=4.0picap_s(k)(qvoqsiz(k)-1.0)tmpc/(tmpa+tmpb)
959	tmp1=av_s(k)sqrho(k)olambdas(k)*(5+bv_s(k)+2mu_s)/visc(k)
960
961	!---- YLIN, here there is some approximation assuming mu_s =1, so gamma(2)=1, etc.
962
963	tmp2= abiN0_s(k)( vf1solambdas(k)olambdas(k)+ &
964	vf2sschmidt(k)0.33334 &
965	ggamma(2.5+0.5bv_s(k)+mu_s)sqrt(tmp1) )
966
967	tmp3=odtb*( qvz(k)-qsiz(k) )
968	!
969	if( tmp2 .le. 0.0) then
970	tmp2=amax1( tmp2,tmp3)
971	pssub(k)=amax1( tmp2,-qszodt(k) )
972	psdep(k)=0.0
973	else
974	psdep(k)=amin1( tmp2,tmp3 )
975	pssub(k)=0.0
976	end if
977
978	!
979	if(qrz(k) .le. 0.0) go to 1200
980	!
981	! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
982	! these two terms need to be refined in the future, they should be equal
983	!c
984	!c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
985	!c
986	esr=1.0
987	tmpa=olambdar(k)*olambdar(k)
988	tmpb=olambdas(k)*olambdas(k)
989	tmpc=olambdar(k)*olambdas(k)
990	tmp1=pipiesrxnorN0_s(k)abs( vtr(k)-vts(k) )orho(k)
991	tmp2=tmpbtmpbolambdar(k)(5.0tmpb+2.0tmpc+0.5tmpa)
992	tmp3=tmp1rhosnowtmp2
993	pracs(k)=amin1( tmp3,qszodt(k) )
994	pracs(k)=0.0
995	!c
996	!c (10) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
997	!c
998	tmp3=tmpatmpaolambdas(k)(5.0tmpa+2.0tmpc+0.5tmpb)
999	tmp4=tmp1rhowatertmp3
1000	psacr(k)=amin1( tmp4,qrzodt(k) )
1001	!
1002	!c
1003	!c (2) FREEZING OF RAIN TO FORM GRAUPEL (pgfr): Lin (45), added to PI
1004	!c positive value
1005	!c Constant in Bigg freezing Aplume=Ap=0.66 /k
1006	!c Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
1007	!
1008
1009	if (qrz(k) .gt. 1.e-8 ) then
1010	Bp=100.
1011	Ap=0.66
1012	tmp1=olambdar(k)olambdar(k)olambdar(k)
1013	tmp2=20.pipiBpxnorrhowaterorho(k)* &
1014	(exp(-Aptemcc(k))-1.0)tmp1tmp1olambdar(k)
1015	pgfr(k)=amin1( tmp2,qrzodt(k) )
1016	else
1017	pgfr(k)=0
1018	endif
1019
1020	1200 continue
1021	!
1022
1023	else
1024
1025	!
1026	!***********************************************************************
1027	!******* snow production processes for T > 0 C ********
1028	!***********************************************************************
1029	!
1030	if (qsz(k) .le. 0.0) go to 1400
1031	!c
1032	!c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1033	!c
1034	esw=1.0
1035
1036	save1 =aa_s(k)sqrho(k)N0_s(k)* &
1037	ggamma(bv_s(k)+tmp_sa(k))olambdas(k)*(bv_s(k)+tmp_sa(k))
1038
1039	tmp1=esw*save1
1040	psacw(k)=qlzodt(k)( 1.0-exp(-tmp1dtb) )
1041
1042	!c
1043	!c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1044	!c
1045	esr=1.0
1046	tmpa=olambdar(k)*olambdar(k)
1047	tmpb=olambdas(k)*olambdas(k)
1048	tmpc=olambdar(k)*olambdas(k)
1049	tmp1=pipiesrxnorN0_s(k)abs( vtr(k)-vts(k) )orho(k)
1050	tmp2=tmpatmpaolambdas(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1051	tmp3=tmp1rhowatertmp2
1052	psacr(k)=amin1( tmp3,qrzodt(k) )
1053	!c
1054	!c (3) MELTING OF SNOW (Psmlt): Lin (32)
1055	!c Psmlt is negative value
1056	!
1057	delrs=rs0(k)-qvz(k)
1058	term1=2.0piorho(k)( xlvdiffwv(k)rho(k)delrs- &
1059	xka(k)*temcc(k) )
1060	tmp1= av_s(k)sqrho(k)olambdas(k)*(5+bv_s(k)+2mu_s)/visc(k)
1061	tmp2= N0_s(k)( vf1solambdas(k)*olambdas(k)+ &
1062	vf2sschmidt(k)0.33334 &
1063	ggamma(2.5+0.5bv_s(k)+mu_s)sqrt(tmp1) )
1064	tmp3=term1oxlftmp2-cwoxlftemcc(k)( psacw(k)+psacr(k) )
1065	tmp4=amin1(0.0,tmp3)
1066	psmlt(k)=amax1( tmp4,-qszodt(k) )
1067	!c
1068	!c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
1069	!c but use Lin et al. coefficience
1070	!c Psmltevp is a negative value
1071	!c
1072	tmpa=rvaporxka(k)tem(k)*tem(k)
1073	tmpb=xlvxlvrho(k)qswz(k)diffwv(k)
1074	tmpc=tmpaqswz(k)diffwv(k)
1075	tmpd=amin1( 0.0,(qvoqswz(k)-0.90)qswz(k)odtb )
1076
1077	abr=2.0pi(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1078	!
1079	!**** allow evaporation to occur when RH less than 90%
1080	!**** here not using 100% because the evaporation cooling
1081	!**** of temperature is not taking into account yet; hence,
1082	!**** the qsw value is a little bit larger. This will avoid
1083	!**** evaporation can generate cloud.
1084	!
1085	tmp1=av_s(k)sqrho(k)olambdas(k)*(5+bv_s(k)+2mu_s)/visc(k)
1086	tmp2= N0_s(k)( vf1solambdas(k)*olambdas(k)+ &
1087	vf2sschmidt(k)0.33334 &
1088	ggamma(2.5+0.5bv_s(k)+mu_s)sqrt(tmp1) )
1089	tmp3=amin1(0.0,tmp2)
1090	tmp3=amax1( tmp3,tmpd )
1091	psmltevp(k)=amax1( tmp3,-qszodt(k) )
1092	1400 continue
1093	!
1094	end if !---- end of snow/ice processes
1095
1096	! CALL wrf_debug ( 100 , 'module_ylin: finish ice/snow processes' )
1097
1098
1099	!***********************************************************************
1100	!******* rain production processes ********
1101	!***********************************************************************
1102
1103	!c
1104	!c (1) AUTOCONVERSION OF RAIN (Praut): using Liu and Daum (2004)
1105	!c
1106
1107	!---- YLIN, autoconversion use Liu and Daum (2004), unit = g cm-3 s-1, in the scheme kg/kg s-1, so
1108
1109	if (qlz(k) .gt. 1e-6) then
1110	lamc(k) = (Nt_crhowaterpiggamma(4.+mu_c)/(6.rho(k)qlz(k))/ & !--- N(D) = N0D^muexp(-lamcD)
1111	ggamma(1+mu_c))**0.3333
1112	Dc_liu = (ggamma(6+1+mu_c)/ggamma(1+mu_c))**(1./6.)/lamc(k) !----- R6 in m
1113
1114	if (Dc_liu .gt. R6c) then
1115	disp = 1./(mu_c+1.) !--- square of relative dispersion
1116	eta = (0.75/pi/(1e-3rhowater))21.9e11((1+3disp)(1+4disp)*&
1117	(1+5disp)/(1+disp)/(1+2disp))
1118	praut(k) = eta(1e-3rho(k)qlz(k))3/(1e-6Nt_c) !--- g cm-3 s-1
1119	praut(k) = praut(k)/(1e-3*rho(k)) !--- kg kg-1 s-1
1120	else
1121	praut(k) = 0.0
1122	endif
1123	else
1124	praut(k) = 0.0
1125	endif
1126
1127	!c
1128	!c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
1129	!c
1130	erw=1.0
1131
1132	tmp1=pio4erwxnorav_rsqrho(k)* &
1133	gambp3olambdar(k)*bp3
1134	pracw(k)=qlzodt(k)( 1.0-exp(-tmp1dtb) )
1135
1136	!c
1137	!c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
1138	!c Prevp is negative value
1139	!c
1140	!c Sw=qvoqsw : saturation ratio
1141	!c
1142	tmpa=rvaporxka(k)tem(k)*tem(k)
1143	tmpb=xlvxlvrho(k)qswz(k)diffwv(k)
1144	tmpc=tmpaqswz(k)diffwv(k)
1145	tmpd=amin1(0.0,(qvoqswz(k)-0.90)qswz(k)odtb)
1146
1147	abr=2.0pi(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1148	tmp1=av_rsqrho(k)olambdar(k)**bp5/visc(k)
1149	tmp2=abrxnor( vf1rolambdar(k)olambdar(k)+ &
1150	vf2rschmidt(k)0.33334gambp5o2*sqrt(tmp1) )
1151	tmp3=amin1( 0.0,tmp2 )
1152	tmp3=amax1( tmp3,tmpd )
1153	prevp(k)=amax1( tmp3,-qrzodt(k) )
1154
1155	! CALL wrf_debug ( 100 , 'module_ylin: finish rain processes' )
1156
1157	!c
1158	!c**********************************************************************
1159	!c*** combine all processes together and avoid negative ***
1160	!c***** water substances
1161	!***********************************************************************
1162	!c
1163	if ( temcc(k) .lt. 0.0) then
1164	!c
1165	!c combined water vapor depletions
1166	!c
1167	tmp=psdep(k)
1168	if ( tmp .gt. qvzodt(k) ) then
1169	factor=qvzodt(k)/tmp
1170	psdep(k)=psdep(k)*factor
1171	end if
1172	!c
1173	!c combined cloud water depletions
1174	!c
1175	tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)
1176	if ( tmp .gt. qlzodt(k) ) then
1177	factor=qlzodt(k)/tmp
1178	praut(k)=praut(k)*factor
1179	psacw(k)=psacw(k)*factor
1180	psfw(k)=psfw(k)*factor
1181	pracw(k)=pracw(k)*factor
1182	end if
1183	!c
1184	!c combined cloud ice depletions
1185	!c
1186	tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)
1187	if (tmp .gt. qizodt(k) ) then
1188	factor=qizodt(k)/tmp
1189	psaut(k)=psaut(k)*factor
1190	psaci(k)=psaci(k)*factor
1191	praci(k)=praci(k)*factor
1192	psfi(k)=psfi(k)*factor
1193	endif
1194	!c
1195	!c combined all rain processes
1196	!c
1197	tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k)+pgfr(k)
1198	if (tmp_r .gt. qrzodt(k) ) then
1199	factor=qrzodt(k)/tmp_r
1200	piacr(k)=piacr(k)*factor
1201	psacr(k)=psacr(k)*factor
1202	prevp(k)=prevp(k)*factor
1203	pgfr(k)=pgfr(k)*factor
1204	endif
1205	!c
1206	!c combined all snow processes
1207	!c
1208	tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+pgfr(k)+ &
1209	psfi(k)+praci(k)+piacr(k)+ &
1210	psdep(k)+psacr(k)-pracs(k))
1211	if ( tmp_s .gt. qszodt(k) ) then
1212	factor=qszodt(k)/tmp_s
1213	pssub(k)=pssub(k)*factor
1214	Pracs(k)=Pracs(k)*factor
1215	endif
1216
1217	!c
1218	!c calculate new water substances, thetae, tem, and qvsbar
1219	!c
1220
1221	pvapor(k)=-pssub(k)-psdep(k)-prevp(k)
1222	qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
1223	pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)
1224	qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1225	pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)
1226	qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1227	tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k)+pgfr(k)-pracs(k)
1228	prain(k)=-tmp_r
1229	qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1230	tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+pgfr(k)+ &
1231	psfi(k)+praci(k)+piacr(k)+ &
1232	psdep(k)+psacr(k)-pracs(k))
1233	psnow(k)=-tmp_s
1234	qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1235
1236	qschg(k)=qschg(k)+psnow(k)
1237	qschg(k)=psnow(k)
1238
1239	tmp=ocp/tothz(k)xLfqschg(k)
1240	theiz(k)=theiz(k)+dtb*tmp
1241	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1242	tem(k)=thz(k)*tothz(k)
1243
1244	temcc(k)=tem(k)-273.15
1245
1246	if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
1247	qlpqi=qlz(k)+qiz(k)
1248	if ( qlpqi .eq. 0.0 ) then
1249	qvsbar(k)=qsiz(k)
1250	else
1251	qvsbar(k)=( qiz(k)qsiz(k)+qlz(k)qswz(k) )/qlpqi
1252	endif
1253
1254	!
1255	else !>0 C
1256	!c
1257	!c combined cloud water depletions
1258	!c
1259	tmp=praut(k)+psacw(k)+pracw(k)
1260	if ( tmp .gt. qlzodt(k) ) then
1261	factor=qlzodt(k)/tmp
1262	praut(k)=praut(k)*factor
1263	psacw(k)=psacw(k)*factor
1264	pracw(k)=pracw(k)*factor
1265	end if
1266	!c
1267	!c combined all snow processes
1268	!c
1269	tmp_s=-(psmlt(k)+psmltevp(k))
1270	if (tmp_s .gt. qszodt(k) ) then
1271	factor=qszodt(k)/tmp_s
1272	psmlt(k)=psmlt(k)*factor
1273	psmltevp(k)=psmltevp(k)*factor
1274	endif
1275	!c
1276	!c combined all rain processes
1277	!c
1278	tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k))
1279	if (tmp_r .gt. qrzodt(k) ) then
1280	factor=qrzodt(k)/tmp_r
1281	prevp(k)=prevp(k)*factor
1282	endif
1283	!c
1284	!c calculate new water substances and thetae
1285	!c
1286	pvapor(k)=-psmltevp(k)-prevp(k)
1287	qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
1288	pclw(k)=-praut(k)-pracw(k)-psacw(k)
1289	qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1290	pcli(k)=0.0
1291	qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1292	tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k))
1293	prain(k)=-tmp_r
1294	tmpqrz=qrz(k)
1295	qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1296	tmp_s=-(psmlt(k)+psmltevp(k))
1297	psnow(k)=-tmp_s
1298	qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1299	qschg(k)=psnow(k)
1300	!
1301	tmp=ocp/tothz(k)xLfqschg(k)
1302	theiz(k)=theiz(k)+dtb*tmp
1303	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1304
1305	tem(k)=thz(k)*tothz(k)
1306	temcc(k)=tem(k)-273.15
1307	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
1308	qswz(k)=ep2*es/(prez(k)-es)
1309	qsiz(k)=qswz(k)
1310	qvsbar(k)=qswz(k)
1311	!
1312	end if
1313	! CALL wrf_debug ( 100 , 'module_ylin: finish sum of all processes' )
1314
1315	!
1316	!***********************************************************************
1317	!******** saturation adjustment ********
1318	!***********************************************************************
1319	!
1320	! allow supersaturation exits linearly from 0% at 500 mb to 50%
1321	! above 300 mb
1322	! 5.0e-5=1.0/(500mb-300mb)
1323	!
1324	rsat=1.0
1325	if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
1326
1327	!c
1328	!c unsaturated
1329	!c
1330	qvz(k)=qvz(k)+qlz(k)+qiz(k)
1331	qlz(k)=0.0
1332	qiz(k)=0.0
1333
1334	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1335	tem(k)=thz(k)*tothz(k)
1336	temcc(k)=tem(k)-273.15
1337
1338	go to 1800
1339	!
1340	else
1341	!c
1342	!c saturated
1343	!c
1344	pladj(k)=qlz(k)
1345	piadj(k)=qiz(k)
1346	!
1347
1348	CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
1349	k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
1350
1351	!
1352	pladj(k)=odtb*(qlz(k)-pladj(k))
1353	piadj(k)=odtb*(qiz(k)-piadj(k))
1354	!
1355	pclw(k)=pclw(k)+pladj(k)
1356	pcli(k)=pcli(k)+piadj(k)
1357	pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
1358	!
1359	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1360	tem(k)=thz(k)*tothz(k)
1361
1362	temcc(k)=tem(k)-273.15
1363
1364	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
1365	qswz(k)=ep2*es/(prez(k)-es)
1366	if (tem(k) .lt. 273.15 ) then
1367	es=1000.svp1exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
1368	qsiz(k)=ep2*es/(prez(k)-es)
1369	if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
1370	else
1371	qsiz(k)=qswz(k)
1372	endif
1373	qlpqi=qlz(k)+qiz(k)
1374	if ( qlpqi .eq. 0.0 ) then
1375	qvsbar(k)=qsiz(k)
1376	else
1377	qvsbar(k)=( qiz(k)qsiz(k)+qlz(k)qswz(k) )/qlpqi
1378	endif
1379
1380	end if
1381
1382	!
1383	!***********************************************************************
1384	!*** melting and freezing of cloud ice and cloud water ***
1385	!***********************************************************************
1386	qlpqi=qlz(k)+qiz(k)
1387	if(qlpqi .le. 0.0) go to 1800
1388	!
1389	!c
1390	!c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
1391	!c
1392	if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
1393	!c
1394	!c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
1395	!c
1396	if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
1397	!c
1398	!c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
1399	!c this process only considered when -31 C < T < 0 C
1400	!c
1401	if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
1402	!c!
1403	!c! parama1 and parama2 functions must be user supplied
1404	!c!
1405	a1=parama1( temcc(k) )
1406	a2=parama2( temcc(k) )
1407	!! change unit from cgs to mks
1408	a1=a10.001*(1.0-a2)
1409	xnin=xni0exp(-bnitemcc(k))
1410	pidw(k)=xninorho(k)(a1xmnin*a2)
1411	end if
1412	!
1413	pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
1414	pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
1415	qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
1416	qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
1417
1418	!
1419	CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
1420	k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
1421
1422	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1423	tem(k)=thz(k)*tothz(k)
1424
1425	temcc(k)=tem(k)-273.15
1426
1427	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
1428	qswz(k)=ep2*es/(prez(k)-es)
1429
1430	if (tem(k) .lt. 273.15 ) then
1431	es=1000.svp1exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
1432	qsiz(k)=ep2*es/(prez(k)-es)
1433	if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
1434	else
1435	qsiz(k)=qswz(k)
1436	endif
1437	qlpqi=qlz(k)+qiz(k)
1438
1439	if ( qlpqi .eq. 0.0 ) then
1440	qvsbar(k)=qsiz(k)
1441	else
1442	qvsbar(k)=( qiz(k)qsiz(k)+qlz(k)qswz(k) )/qlpqi
1443	endif
1444
1445	1800 continue
1446	!
1447	!***********************************************************************
1448	!******** integrate the productions of rain and snow ********
1449	!***********************************************************************
1450	!
1451	2000 continue
1452
1453	!
1454	!**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
1455	!
1456	do k=kts+1,kte
1457	if ( qvz(k) .lt. qvmin ) then
1458	qlz(k)=0.0
1459	qiz(k)=0.0
1460	qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
1461	end if
1462	enddo
1463	!
1464
1465	! CALL wrf_debug ( 100 , 'module_ylin: finish saturation adjustment' )
1466
1467	END SUBROUTINE clphy1d_ylin
1468
1469
1470
1471
1472	!---------------------------------------------------------------------
1473	! SATURATED ADJUSTMENT
1474	!---------------------------------------------------------------------
1475	SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, &
1476	kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
1477	!---------------------------------------------------------------------
1478	IMPLICIT NONE
1479	!---------------------------------------------------------------------
1480	! This program use Newton's method for finding saturated temperature
1481	! and saturation mixing ratio.
1482	!
1483	! In this saturation adjustment scheme we assume
1484	! (1) the saturation mixing ratio is the mass weighted average of
1485	! saturation values over liquid water (qsw), and ice (qsi)
1486	! following Lord et al., 1984 and Tao, 1989
1487	!
1488	! (2) the percentage of cloud liquid and cloud ice will
1489	! be fixed during the saturation calculation
1490	!---------------------------------------------------------------------
1491	!
1492
1493	INTEGER, INTENT(IN ) :: kts, kte, k
1494
1495	REAL, DIMENSION( kts:kte ), &
1496	INTENT(INOUT) :: qvz, qlz, qiz
1497	!
1498	REAL, DIMENSION( kts:kte ), &
1499	INTENT(IN ) :: prez, theiz, tothz
1500
1501	REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
1502	REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
1503
1504	! LOCAL VARS
1505
1506	INTEGER :: n
1507
1508	REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
1509	qswz, qvsbar
1510
1511	REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
1512	denom1, denom2, dqvsbar, ftsat, dftsat, qpz,es
1513	!
1514	!---------------------------------------------------------------------
1515
1516	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1517
1518	tem(k)=tothz(k)*thz(k)
1519	if (tem(k) .gt. 273.15) then
1520	! qsat=episp0k/prez(k)* &
1521	! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
1522	es=1000.svp1exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
1523	qsat=ep2*es/(prez(k)-es)
1524	else
1525	qsat=episp0k/prez(k)* &
1526	exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
1527	end if
1528	qpz=qvz(k)+qlz(k)+qiz(k)
1529	if (qpz .lt. qsat) then
1530	qvz(k)=qpz
1531	qiz(k)=0.0
1532	qlz(k)=0.0
1533	go to 400
1534	end if
1535	qlpqi=qlz(k)+qiz(k)
1536	if( qlpqi .ge. 1.0e-5) then
1537	ratql=qlz(k)/qlpqi
1538	ratqi=qiz(k)/qlpqi
1539	else
1540	t0=273.15
1541	! t1=233.15
1542	t1=248.15
1543	tmp1=( t0-tem(k) )/(t0-t1)
1544	tmp1=amin1(1.0,tmp1)
1545	tmp1=amax1(0.0,tmp1)
1546	ratqi=tmp1
1547	ratql=1.0-tmp1
1548	end if
1549	!
1550	!
1551	!-- saturation mixing ratios over water and ice
1552	!-- at the outset we will follow Bolton 1980 MWR for
1553	!-- the water and Murray JAS 1967 for the ice
1554	!
1555	!-- dqvsbar=d(qvsbar)/dT
1556	!-- ftsat=F(Tsat)
1557	!-- dftsat=d(F(T))/dT
1558	!
1559	! First guess of tsat
1560
1561	tsat=tem(k)
1562	absft=1.0
1563	!
1564	do 200 n=1,20
1565	denom1=1.0/(tsat-svp3)
1566	denom2=1.0/(tsat-7.66)
1567	! qswz(k)=episp0k/prez(k)* &
1568	! exp( svp2denom1(tsat-273.15) )
1569	es=1000.svp1exp( svp2denom1(tsat-svpt0) )
1570	qswz(k)=ep2*es/(prez(k)-es)
1571	if (tem(k) .lt. 273.15) then
1572	! qsiz(k)=episp0k/prez(k)* &
1573	! exp( 21.8745584denom2(tsat-273.15) )
1574	es=1000.svp1exp( 21.8745584denom2(tsat-273.15) )
1575	qsiz(k)=ep2*es/(prez(k)-es)
1576	if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
1577	else
1578	qsiz(k)=qswz(k)
1579	endif
1580	qvsbar(k)=ratqlqswz(k)+ratqiqsiz(k)
1581	!
1582	! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
1583	if( absft .lt. 0.01 ) go to 300
1584	!
1585	dqvsbar=ratqlqswz(k)svp2243.5denom1*denom1+ &
1586	ratqiqsiz(k)21.8745584265.5denom2*denom2
1587	ftsat=tsat+(xlvocp+ratqixlfocp)qvsbar(k)- &
1588	tothz(k)theiz(k)-xlfocpratqi*(qvz(k)+qlz(k)+qiz(k))
1589	dftsat=1.0+(xlvocp+ratqixlfocp)dqvsbar
1590	tsat=tsat-ftsat/dftsat
1591	absft=abs(ftsat)
1592
1593	200 continue
1594	9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
1595	300 continue
1596
1597	if( qpz .gt. qvsbar(k) ) then
1598	qvz(k)=qvsbar(k)
1599	qiz(k)=ratqi*( qpz-qvz(k) )
1600	qlz(k)=ratql*( qpz-qvz(k) )
1601	else
1602	qvz(k)=qpz
1603	qiz(k)=0.0
1604	qlz(k)=0.0
1605	end if
1606	400 continue
1607
1608	END SUBROUTINE satadj
1609
1610
1611	!----------------------------------------------------------------
1612	REAL FUNCTION parama1(temp)
1613	!----------------------------------------------------------------
1614	IMPLICIT NONE
1615	!----------------------------------------------------------------
1616	! This program calculate the parameter for crystal growth rate
1617	! in Bergeron process
1618	!----------------------------------------------------------------
1619
1620	REAL, INTENT (IN ) :: temp
1621	REAL, DIMENSION(32) :: a1
1622	INTEGER :: i1, i1p1
1623	REAL :: ratio
1624
1625	data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
1626	0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
1627	0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
1628	0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
1629	0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
1630	0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
1631	0.5333e-6,0.4834e-6/
1632
1633	i1=int(-temp)+1
1634	i1p1=i1+1
1635	ratio=-(temp)-float(i1-1)
1636	parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
1637
1638	END FUNCTION parama1
1639
1640	!----------------------------------------------------------------
1641	REAL FUNCTION parama2(temp)
1642	!----------------------------------------------------------------
1643	IMPLICIT NONE
1644	!----------------------------------------------------------------
1645	! This program calculate the parameter for crystal growth rate
1646	! in Bergeron process
1647	!----------------------------------------------------------------
1648
1649	REAL, INTENT (IN ) :: temp
1650	REAL, DIMENSION(32) :: a2
1651	INTEGER :: i1, i1p1
1652	REAL :: ratio
1653
1654	data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
1655	0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
1656	0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
1657	0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
1658	0.4433,0.4413,0.4382,0.4361/
1659	i1=int(-temp)+1
1660	i1p1=i1+1
1661	ratio=-(temp)-float(i1-1)
1662	parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
1663
1664	END FUNCTION parama2
1665
1666	!+---+-----------------------------------------------------------------+
1667	! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
1668	! A FUNCTION OF TEMPERATURE AND PRESSURE
1669	!
1670	REAL FUNCTION RSLF(P,T)
1671
1672	IMPLICIT NONE
1673	REAL, INTENT(IN):: P, T
1674	REAL:: ESL,X
1675	REAL, PARAMETER:: C0= .611583699E03
1676	REAL, PARAMETER:: C1= .444606896E02
1677	REAL, PARAMETER:: C2= .143177157E01
1678	REAL, PARAMETER:: C3= .264224321E-1
1679	REAL, PARAMETER:: C4= .299291081E-3
1680	REAL, PARAMETER:: C5= .203154182E-5
1681	REAL, PARAMETER:: C6= .702620698E-8
1682	REAL, PARAMETER:: C7= .379534310E-11
1683	REAL, PARAMETER:: C8=-.321582393E-13
1684
1685	X=MAX(-80.,T-273.16)
1686
1687	! ESL=612.2EXP(17.67X/(T-29.65))
1688	ESL=C0+X(C1+X(C2+X(C3+X(C4+X(C5+X(C6+X(C7+XC8)))))))
1689	RSLF=.622*ESL/(P-ESL)
1690
1691	END FUNCTION RSLF
1692	!
1693	!+---+-----------------------------------------------------------------+
1694	! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
1695	! FUNCTION OF TEMPERATURE AND PRESSURE
1696	!
1697	REAL FUNCTION RSIF(P,T)
1698
1699	IMPLICIT NONE
1700	REAL, INTENT(IN):: P, T
1701	REAL:: ESI,X
1702	REAL, PARAMETER:: C0= .609868993E03
1703	REAL, PARAMETER:: C1= .499320233E02
1704	REAL, PARAMETER:: C2= .184672631E01
1705	REAL, PARAMETER:: C3= .402737184E-1
1706	REAL, PARAMETER:: C4= .565392987E-3
1707	REAL, PARAMETER:: C5= .521693933E-5
1708	REAL, PARAMETER:: C6= .307839583E-7
1709	REAL, PARAMETER:: C7= .105785160E-9
1710	REAL, PARAMETER:: C8= .161444444E-12
1711
1712	X=MAX(-80.,T-273.16)
1713	ESI=C0+X(C1+X(C2+X(C3+X(C4+X(C5+X(C6+X(C7+XC8)))))))
1714	RSIF=.622*ESI/(P-ESI)
1715
1716	END FUNCTION RSIF
1717	!+---+-----------------------------------------------------------------+
1718
1719	!----------------------------------------------------------------
1720	REAL FUNCTION ggamma(X)
1721	!----------------------------------------------------------------
1722	IMPLICIT NONE
1723	!----------------------------------------------------------------
1724	REAL, INTENT(IN ) :: x
1725	REAL, DIMENSION(8) :: B
1726	INTEGER ::j, K1
1727	REAL ::PF, G1TO2 ,TEMP
1728
1729	DATA B/-.577191652,.988205891,-.897056937,.918206857, &
1730	-.756704078,.482199394,-.193527818,.035868343/
1731
1732	PF=1.
1733	TEMP=X
1734	DO 10 J=1,200
1735	IF (TEMP .LE. 2) GO TO 20
1736	TEMP=TEMP-1.
1737	10 PF=PF*TEMP
1738	! 100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
1739	! WRITE(wrf_err_message,100)X
1740	! CALL wrf_error_fatal(wrf_err_message)
1741	20 G1TO2=1.
1742	TEMP=TEMP - 1.
1743	DO 30 K1=1,8
1744	30 G1TO2=G1TO2 + B(K1)TEMP*K1
1745	ggamma=PF*G1TO2
1746
1747	END FUNCTION ggamma
1748
1749	!----------------------------------------------------------------
1750
1751	END MODULE module_mp_sbu_ylin
1752

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

Download in other formats:

Original Format