Context Navigation

module_mp_lin.F

Last change on this file 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: 83.8 KB

Line
1	!WRF:MODEL_LAYER:PHYSICS
2	!
3
4	MODULE module_mp_lin
5
6	USE module_wrf_error
7	!
8	REAL , PARAMETER, PRIVATE :: RH = 1.0
9	! REAL , PARAMETER, PRIVATE :: episp0 = 0.622*611.21
10	REAL , PARAMETER, PRIVATE :: xnor = 8.0e6
11	REAL , PARAMETER, PRIVATE :: xnos = 3.0e6
12
13	! Lin
14	! REAL , PARAMETER, PRIVATE :: xnog = 4.0e4
15	! REAL , PARAMETER, PRIVATE :: rhograul = 917.
16
17	! Hobbs
18	REAL , PARAMETER, PRIVATE :: xnog = 4.0e6
19	REAL , PARAMETER, PRIVATE :: rhograul = 400.
20
21	!
22	REAL , PARAMETER, PRIVATE :: &
23	qi0 = 1.0e-3, ql0 = 7.0e-4, qs0 = 6.0E-4, &
24	xmi50 = 4.8e-10, xmi40 = 2.46e-10, &
25	constb = 0.8, constd = 0.25, &
26	o6 = 1./6., cdrag = 0.6, &
27	avisc = 1.49628e-6, adiffwv = 8.7602e-5, &
28	axka = 1.4132e3, di50 = 1.0e-4, xmi = 4.19e-13, &
29	cw = 4.187e3, vf1s = 0.78, vf2s = 0.31, &
30	xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5, &
31	ci = 2.093e3
32	CONTAINS
33
34	!-------------------------------------------------------------------
35	! Lin et al., 1983, JAM, 1065-1092, and
36	! Rutledge and Hobbs, 1984, JAS, 2949-2972
37	!-------------------------------------------------------------------
38	SUBROUTINE lin_et_al(th &
39	,qv, ql, qr &
40	,qi, qs &
41	,rho, pii, p &
42	,dt_in &
43	,z,ht, dz8w &
44	,grav, cp, Rair, rvapor &
45	,XLS, XLV, XLF, rhowater, rhosnow &
46	,EP2,SVP1,SVP2,SVP3,SVPT0 &
47	, RAINNC, RAINNCV &
48	,ids,ide, jds,jde, kds,kde &
49	,ims,ime, jms,jme, kms,kme &
50	,its,ite, jts,jte, kts,kte &
51	! Optional
52	,qlsink, precr, preci, precs, precg &
53	, F_QG,F_QNDROP &
54	, qg, qndrop &
55	)
56	!-------------------------------------------------------------------
57	IMPLICIT NONE
58	!-------------------------------------------------------------------
59	!
60	! Shuhua 12/17/99
61	!
62	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
63	ims,ime, jms,jme, kms,kme , &
64	its,ite, jts,jte, kts,kte
65
66	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
67	INTENT(INOUT) :: &
68	th, &
69	qv, &
70	ql, &
71	qr
72
73	!
74	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
75	INTENT(IN ) :: &
76	rho, &
77	pii, &
78	p, &
79	dz8w
80
81	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
82	INTENT(IN ) :: z
83
84
85
86	REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
87
88	REAL, INTENT(IN ) :: dt_in, &
89	grav, &
90	Rair, &
91	rvapor, &
92	cp, &
93	XLS, &
94	XLV, &
95	XLF, &
96	rhowater, &
97	rhosnow
98
99	REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
100
101	REAL, DIMENSION( ims:ime , jms:jme ), &
102	INTENT(INOUT) :: RAINNC, &
103	RAINNCV
104
105	! Optional
106
107	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
108	OPTIONAL, &
109	INTENT(INOUT) :: &
110	qi, &
111	qs, &
112	qg, &
113	qndrop
114
115	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
116	OPTIONAL, INTENT(OUT ) :: &
117	qlsink, & ! cloud water conversion to rain (/s)
118	precr, & ! rain precipitation rate at all levels (kg/m2/s)
119	preci, & ! ice precipitation rate at all levels (kg/m2/s)
120	precs, & ! snow precipitation rate at all levels (kg/m2/s)
121	precg ! graupel precipitation rate at all levels (kg/m2/s)
122
123	LOGICAL, INTENT(IN), OPTIONAL :: F_QG, F_QNDROP
124
125	! LOCAL VAR
126
127	INTEGER :: min_q, max_q
128
129	REAL, DIMENSION( its:ite , jts:jte ) &
130	:: rain, snow, graupel,ice
131
132	REAL, DIMENSION( kts:kte ) :: qvz, qlz, qrz, &
133	qiz, qsz, qgz, &
134	thz, &
135	tothz, rhoz, &
136	orhoz, sqrhoz, &
137	prez, zz, &
138	precrz, preciz, precsz, precgz, &
139	qndropz, &
140	dzw, preclw
141
142	LOGICAL :: flag_qg, flag_qndrop
143	!
144	REAL :: dt, pptrain, pptsnow, pptgraul, rhoe_s, &
145	gindex, pptice
146	real :: qndropconst
147
148	INTEGER :: i,j,k
149	!
150	flag_qg = .false.
151	flag_qndrop = .false.
152	IF ( PRESENT ( f_qg ) ) flag_qg = f_qg
153	IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
154	!
155	dt=dt_in
156	rhoe_s=1.29
157	qndropconst=100.e6 !sg
158	gindex=1.0
159
160	IF (.not.flag_qg) gindex=0.
161
162	j_loop: DO j = jts, jte
163	i_loop: DO i = its, ite
164	!
165	!- write data from 3-D to 1-D
166	!
167	DO k = kts, kte
168	qvz(k)=qv(i,k,j)
169	qlz(k)=ql(i,k,j)
170	qrz(k)=qr(i,k,j)
171	thz(k)=th(i,k,j)
172	!
173	rhoz(k)=rho(i,k,j)
174	orhoz(k)=1./rhoz(k)
175	prez(k)=p(i,k,j)
176	sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
177	tothz(k)=pii(i,k,j)
178	zz(k)=z(i,k,j)
179	dzw(k)=dz8w(i,k,j)
180	END DO
181
182	IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
183	DO k = kts, kte
184	qndropz(k)=qndrop(i,k,j)
185	ENDDO
186	ELSE
187	DO k = kts, kte
188	qndropz(k)=qndropconst
189	ENDDO
190	ENDIF
191
192	DO k = kts, kte
193	qiz(k)=qi(i,k,j)
194	qsz(k)=qs(i,k,j)
195	ENDDO
196
197	IF ( flag_qg .AND. PRESENT( qg ) ) THEN
198	DO k = kts, kte
199	qgz(k)=qg(i,k,j)
200	ENDDO
201	ELSE
202	DO k = kts, kte
203	qgz(k)=0.
204	ENDDO
205	ENDIF
206	!
207	pptrain=0.
208	pptsnow=0.
209	pptgraul=0.
210	pptice=0.
211	CALL clphy1d( dt, qvz, qlz, qrz, qiz, qsz, qgz, &
212	qndropz,flag_qndrop, &
213	thz, tothz, rhoz, orhoz, sqrhoz, &
214	prez, zz, dzw, ht(I,J), preclw, &
215	precrz, preciz, precsz, precgz, &
216	pptrain, pptsnow, pptgraul, pptice, &
217	grav, cp, Rair, rvapor, gindex, &
218	XLS, XLV, XLF, rhowater, rhosnow, &
219	EP2,SVP1,SVP2,SVP3,SVPT0, &
220	kts, kte, i, j )
221
222	!
223	! Precipitation from cloud microphysics -- only for one time step
224	!
225	! unit is transferred from m to mm
226
227	!
228	rain(i,j)=pptrain
229	snow(i,j)=pptsnow
230	graupel(i,j)=pptgraul
231	ice(i,j)=pptice
232	!
233	RAINNCV(i,j)= pptrain + pptsnow + pptgraul + pptice
234	RAINNC(i,j)=RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
235
236	!
237	!- update data from 1-D back to 3-D
238	!
239	!
240	IF ( present(qlsink) .and. present(precr) ) THEN !sg beg
241	DO k = kts, kte
242	if(ql(i,k,j)>1.e-20) then
243	qlsink(i,k,j)=-preclw(k)/ql(i,k,j)
244	else
245	qlsink(i,k,j)=0.
246	endif
247	precr(i,k,j)=precrz(k)
248	END DO
249	END IF !sg end
250
251	DO k = kts, kte
252	qv(i,k,j)=qvz(k)
253	ql(i,k,j)=qlz(k)
254	qr(i,k,j)=qrz(k)
255	th(i,k,j)=thz(k)
256	END DO
257	!
258	IF ( flag_qndrop .AND. PRESENT( qndrop ) ) THEN !sg beg
259	DO k = kts, kte
260	qndrop(i,k,j)=qndropz(k)
261	ENDDO
262	END IF !sg end
263
264	DO k = kts, kte
265	qi(i,k,j)=qiz(k)
266	qs(i,k,j)=qsz(k)
267	ENDDO
268
269	IF ( present(preci) ) THEN !sg beg
270	DO k = kts, kte
271	preci(i,k,j)=preciz(k)
272	ENDDO
273	END IF
274
275	IF ( present(precs) ) THEN
276	DO k = kts, kte
277	precs(i,k,j)=precsz(k)
278	ENDDO
279	END IF !sg end
280
281	IF ( flag_qg .AND. PRESENT( qg ) ) THEN
282	DO k = kts, kte
283	qg(i,k,j)=qgz(k)
284	ENDDO
285
286	IF ( present(precg) ) THEN !sg beg
287	DO k = kts, kte
288	precg(i,k,j)=precgz(k)
289	ENDDO !sg end
290	END IF
291	ELSE !sg beg
292	IF ( present(precg) ) precg(i,:,j)=0. !sg end
293	ENDIF
294	!
295	ENDDO i_loop
296	ENDDO j_loop
297
298	END SUBROUTINE lin_et_al
299
300
301	!-----------------------------------------------------------------------
302	SUBROUTINE clphy1d(dt, qvz, qlz, qrz, qiz, qsz, qgz, &
303	qndropz,flag_qndrop, &
304	thz, tothz, rho, orho, sqrho, &
305	prez, zz, dzw, zsfc, preclw, &
306	precrz, preciz, precsz, precgz, &
307	pptrain, pptsnow, pptgraul, &
308	pptice, grav, cp, Rair, rvapor, gindex, &
309	XLS, XLV, XLF, rhowater, rhosnow, &
310	EP2,SVP1,SVP2,SVP3,SVPT0, &
311	kts, kte, i, j )
312	!-----------------------------------------------------------------------
313	IMPLICIT NONE
314	!-----------------------------------------------------------------------
315	! This program handles the vertical 1-D cloud micphysics
316	!-----------------------------------------------------------------------
317	! avisc: constant in empirical formular for dynamic viscosity of air
318	! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
319	! adiffwv: constant in empirical formular for diffusivity of water
320	! vapor in air
321	! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
322	! axka: constant in empirical formular for thermal conductivity of air
323	! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
324	! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
325	! xmi50: mass of a 50 micron ice crystal
326	! = 4.8e-10 [kg] =4.8e-7 [g]
327	! xmi40: mass of a 40 micron ice crystal
328	! = 2.46e-10 [kg] = 2.46e-7 [g]
329	! di50: diameter of a 50 micro (radius) ice crystal
330	! =1.0e-4 [m]
331	! xmi: mass of one cloud ice crystal
332	! =4.19e-13 [kg] = 4.19e-10 [g]
333	! oxmi=1.0/xmi
334	!
335	! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
336	! Hsie et al.(1980) and Rutledge and Hobbs(1983) )
337	! bni=0.5 [K-1]
338	! xmnin: mass of a natural ice nucleus
339	! = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
340	! Hsie et al. (1980)
341	! = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
342	! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
343	! consta: constant in empirical formular for terminal
344	! velocity of raindrop
345	! =2115.0 [cm*(1-b)/s] = 2115.00.01(1-b) [m(1-b)/s]
346	! constb: constant in empirical formular for terminal
347	! velocity of raindrop
348	! =0.8
349	! xnor: intercept parameter of the raindrop size distribution
350	! = 0.08 cm-4 = 8.0e6 m-4
351	! araut: time sacle for autoconversion of cloud water to raindrops
352	! =1.0e-3 [s-1]
353	! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
354	! vf1r: ventilation factors for rain =0.78
355	! vf2r: ventilation factors for rain =0.31
356	! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
357	! constc: constant in empirical formular for terminal
358	! velocity of snow
359	! =152.93 [cm*(1-d)/s] = 152.930.01(1-d) [m(1-d)/s]
360	! constd: constant in empirical formular for terminal
361	! velocity of snow
362	! =0.25
363	! xnos: intercept parameter of the snowflake size distribution
364	! vf1s: ventilation factors for snow =0.78
365	! vf2s: ventilation factors for snow =0.31
366	!
367	!----------------------------------------------------------------------
368
369	INTEGER, INTENT(IN ) :: kts, kte, i, j
370
371	REAL, DIMENSION( kts:kte ), &
372	INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, &
373	qndropz, &
374	qgz, thz
375
376	REAL, DIMENSION( kts:kte ), &
377	INTENT(IN ) :: tothz, rho, orho, sqrho, &
378	prez, zz, dzw
379
380	REAL, INTENT(IN ) :: dt, grav, cp, Rair, rvapor, &
381	XLS, XLV, XLF, rhowater, &
382	rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0
383
384	REAL, DIMENSION( kts:kte ), INTENT(OUT) :: preclw, &
385	precrz, preciz, precsz, precgz
386
387	REAL, INTENT(INOUT) :: pptrain, pptsnow, pptgraul, pptice
388
389	REAL, INTENT(IN ) :: zsfc
390	logical, intent(in) :: flag_qndrop !sg
391
392	! local vars
393
394	REAL :: obp4, bp3, bp5, bp6, odp4, &
395	dp3, dp5, dp5o2
396
397
398	! temperary vars
399
400	REAL :: tmp, tmp0, tmp1, tmp2,tmp3, &
401	tmp4,delta2,delta3, delta4, &
402	tmpa,tmpb,tmpc,tmpd,alpha1, &
403	qic, abi,abr, abg, odtberg, &
404	vti50,eiw,eri,esi,esr, esw, &
405	erw,delrs,term0,term1,araut, &
406	constg2, vf1r, vf2r,alpha2, &
407	Ap, Bp, egw, egi, egs, egr, &
408	constg, gdelta4, g1sdelt4, &
409	factor, tmp_r, tmp_s,tmp_g, &
410	qlpqi, rsat, a1, a2, xnin
411
412	INTEGER :: k
413	!
414	REAL, DIMENSION( kts:kte ) :: oprez, tem, temcc, theiz, qswz, &
415	qsiz, qvoqswz, qvoqsiz, qvzodt, &
416	qlzodt, qizodt, qszodt, qrzodt, &
417	qgzodt
418
419	REAL, DIMENSION( kts:kte ) :: psnow, psaut, psfw, psfi, praci, &
420	piacr, psaci, psacw, psdep, pssub, &
421	pracs, psacr, psmlt, psmltevp, &
422	prain, praut, pracw, prevp, pvapor, &
423	pclw, pladj, pcli, pimlt, pihom, &
424	pidw, piadj, pgraupel, pgaut, &
425	pgfr, pgacw, pgaci, pgacr, pgacs, &
426	pgacip,pgacrp,pgacsp,pgwet, pdry, &
427	pgsub, pgdep, pgmlt, pgmltevp, &
428	qschg, qgchg
429	!
430
431	REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, &
432	schmidt, xka
433
434	REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, &
435	vtrold, vtsold, vtgold, vtiold, &
436	xlambdar, xlambdas, xlambdag, &
437	olambdar, olambdas, olambdag
438
439	REAL :: episp0k, dtb, odtb, pi, pio4, &
440	pio6, oxLf, xLvocp, xLfocp, consta, &
441	constc, ocdrag, gambp4, gamdp4, &
442	gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
443	gambp6, gam3pt5, gam2pt75, gambp5o2,&
444	gamdp5o2, cwoxlf, ocp, xni50, es
445	!
446	REAL :: qvmin=1.e-20
447	REAL :: gindex
448	REAL :: temc1,save1,save2,xni50mx
449
450	! for terminal velocity flux
451
452	INTEGER :: min_q, max_q
453	REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
454	LOGICAL :: notlast
455	!
456
457	!sg: begin
458	! liqconc = liquid water content in gcm^-3
459	! capn = droplet number concentration cm^-3
460	! dis = relative dispersion (dimensionless) between 0.2 and 1.
461	! Written by Yangang Liu based on Liu et al., GRL 32, 2005.
462	! Autoconversion rate P = P0*T
463	! p0 = rate function
464	! kappa = constant in Long kernel
465	! beta = Condensation rate constant
466	! xc = Normalized critical mass
467	! ***********************************************************
468	real liqconc, dis, beta, kappa, p0, xc, capn,rhocgs
469	if(flag_qndrop)then
470	dis = 0.5 ! droplet dispersion, set to 0.5 per SG 8-Nov-2006
471	! Give empirical constants
472	kappa=1.1d10
473	! Calculate Condensation rate constant
474	beta = (1.0d0+3.0d0dis2)(1.0d0+4.0d0dis2) &
475	(1.0d0+5.0d0dis2)/((1.0d0+dis2)(1.0d0+2.0d0dis*2))
476	endif
477	!sg: end
478
479	dtb=dt
480	odtb=1./dtb
481	pi=acos(-1.)
482	pio4=acos(-1.)/4.
483	pio6=acos(-1.)/6.
484	ocp=1./cp
485	oxLf=1./xLf
486	xLvocp=xLv/cp
487	xLfocp=xLf/cp
488	consta=2115.00.01*(1-constb)
489	constc=152.930.01*(1-constd)
490	ocdrag=1./Cdrag
491	! episp0k=RH*episp0
492	episp0k=RHep21000.*svp1
493	!
494	gambp4=ggamma(constb+4.)
495	gamdp4=ggamma(constd+4.)
496	gam4pt5=ggamma(4.5)
497	Cpor=cp/Rair
498	oxmi=1.0/xmi
499	gambp3=ggamma(constb+3.)
500	gamdp3=ggamma(constd+3.)
501	gambp6=ggamma(constb+6)
502	gam3pt5=ggamma(3.5)
503	gam2pt75=ggamma(2.75)
504	gambp5o2=ggamma((constb+5.)/2.)
505	gamdp5o2=ggamma((constd+5.)/2.)
506	cwoxlf=cw/xlf
507	delta2=0.
508	delta3=0.
509	delta4=0.
510	!
511	!-----------------------------------------------------------------------
512	! oprez 1./prez ( prez : pressure)
513	! qsw saturated mixing ratio on water surface
514	! qsi saturated mixing ratio on ice surface
515	! episp0k RHesaturated pressure at 273.15 K
516	! qvoqsw qv/qsw
517	! qvoqsi qv/qsi
518	! qvzodt qv/dt
519	! qlzodt ql/dt
520	! qizodt qi/dt
521	! qszodt qs/dt
522	! qrzodt qr/dt
523	! qgzodt qg/dt
524	!
525	! temcc temperature in dregee C
526	!
527
528	obp4=1.0/(constb+4.0)
529	bp3=constb+3.0
530	bp5=constb+5.0
531	bp6=constb+6.0
532	odp4=1.0/(constd+4.0)
533	dp3=constd+3.0
534	dp5=constd+5.0
535	dp5o2=0.5*(constd+5.0)
536	!
537	do k=kts,kte
538	oprez(k)=1./prez(k)
539	enddo
540
541	do k=kts,kte
542	qlz(k)=amax1( 0.0,qlz(k) )
543	qiz(k)=amax1( 0.0,qiz(k) )
544	qvz(k)=amax1( qvmin,qvz(k) )
545	qsz(k)=amax1( 0.0,qsz(k) )
546	qrz(k)=amax1( 0.0,qrz(k) )
547	qgz(k)=amax1( 0.0,qgz(k) )
548	qndropz(k)=amax1( 0.0,qndropz(k) ) !sg
549	!
550	tem(k)=thz(k)*tothz(k)
551	temcc(k)=tem(k)-273.15
552	!
553	! qswz(k)=episp0koprez(k) &
554	! exp( svp2*temcc(k)/(tem(k)-svp3) )
555	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
556	qswz(k)=ep2*es/(prez(k)-es)
557	if (tem(k) .lt. 273.15 ) then
558	! qsiz(k)=episp0koprez(k) &
559	! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
560	es=1000.svp1exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
561	qsiz(k)=ep2*es/(prez(k)-es)
562	if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
563	else
564	qsiz(k)=qswz(k)
565	endif
566	!
567	qvoqswz(k)=qvz(k)/qswz(k)
568	qvoqsiz(k)=qvz(k)/qsiz(k)
569	qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
570	qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
571	qizodt(k)=amax1( 0.0,odtb*qiz(k) )
572	qszodt(k)=amax1( 0.0,odtb*qsz(k) )
573	qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
574	qgzodt(k)=amax1( 0.0,odtb*qgz(k) )
575
576	theiz(k)=thz(k)+(xlvocpqvz(k)-xlfocpqiz(k))/tothz(k)
577	enddo
578
579
580	!
581	!
582	!-----------------------------------------------------------------------
583	! In this simple stable cloud parameterization scheme, only five
584	! forms of water substance (water vapor, cloud water, cloud ice,
585	! rain and snow are considered. The prognostic variables are total
586	! water (qp),cloud water (ql), and cloud ice (qi). Rain and snow are
587	! diagnosed following Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7).
588	! the micro physics are based on (1) Hsie et al.,1980, JAM, 950-977 ;
589	! (2) Lin et al., 1983, JAM, 1065-1092 ; (3) Rutledge and Hobbs, 1983,
590	! JAS, 1185-1206 ; (4) Rutledge and Hobbs, 1984, JAS, 2949-2972.
591	!-----------------------------------------------------------------------
592	!
593	! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
594	! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
595	! xnor: intercept parameter of the raindrop size distribution
596	! = 0.08 cm-4 = 8.0e6 m-4
597	! xnos: intercept parameter of the snowflake size distribution
598	! = 0.03 cm-4 = 3.0e6 m-4
599	! xnog: intercept parameter of the graupel size distribution
600	! = 4.0e-4 cm-4 = 4.0e4 m-4
601	! consta: constant in empirical formular for terminal
602	! velocity of raindrop
603	! =2115.0 [cm*(1-b)/s] = 2115.00.01(1-b) [m(1-b)/s]
604	! constb: constant in empirical formular for terminal
605	! velocity of raindrop
606	! =0.8
607	! constc: constant in empirical formular for terminal
608	! velocity of snow
609	! =152.93 [cm*(1-d)/s] = 152.930.01(1-d) [m(1-d)/s]
610	! constd: constant in empirical formular for terminal
611	! velocity of snow
612	! =0.25
613	! avisc: constant in empirical formular for dynamic viscosity of air
614	! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
615	! adiffwv: constant in empirical formular for diffusivity of water
616	! vapor in air
617	! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
618	! axka: constant in empirical formular for thermal conductivity of air
619	! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
620	! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
621	! = 1.0e-3 g/g = 1.0e-3 kg/gk
622	! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
623	! = 2.0e-3 g/g = 2.0e-3 kg/gk
624	! qs0: mixing ratio threshold for snow aggregation
625	! = 6.0e-4 g/g = 6.0e-4 kg/gk
626	! xmi50: mass of a 50 micron ice crystal
627	! = 4.8e-10 [kg] =4.8e-7 [g]
628	! xmi40: mass of a 40 micron ice crystal
629	! = 2.46e-10 [kg] = 2.46e-7 [g]
630	! di50: diameter of a 50 micro (radius) ice crystal
631	! =1.0e-4 [m]
632	! xmi: mass of one cloud ice crystal
633	! =4.19e-13 [kg] = 4.19e-10 [g]
634	! oxmi=1.0/xmi
635	!
636
637
638	! if gindex=1.0 include graupel
639	! if gindex=0. no graupel
640	!
641	!
642	do k=kts,kte
643	psnow(k)=0.0
644	psaut(k)=0.0
645	psfw(k)=0.0
646	psfi(k)=0.0
647	praci(k)=0.0
648	piacr(k)=0.0
649	psaci(k)=0.0
650	psacw(k)=0.0
651	psdep(k)=0.0
652	pssub(k)=0.0
653	pracs(k)=0.0
654	psacr(k)=0.0
655	psmlt(k)=0.0
656	psmltevp(k)=0.0
657	!
658	prain(k)=0.0
659	praut(k)=0.0
660	pracw(k)=0.0
661	prevp(k)=0.0
662	!
663	pvapor(k)=0.0
664	!
665	pclw(k)=0.0
666	preclw(k)=0.0 !sg
667	pladj(k)=0.0
668	!
669	pcli(k)=0.0
670	pimlt(k)=0.0
671	pihom(k)=0.0
672	pidw(k)=0.0
673	piadj(k)=0.0
674	enddo
675
676	!
677	!!! graupel
678	!
679	do k=kts,kte
680	pgraupel(k)=0.0
681	pgaut(k)=0.0
682	pgfr(k)=0.0
683	pgacw(k)=0.0
684	pgaci(k)=0.0
685	pgacr(k)=0.0
686	pgacs(k)=0.0
687	pgacip(k)=0.0
688	pgacrP(k)=0.0
689	pgacsp(k)=0.0
690	pgwet(k)=0.0
691	pdry(k)=0.0
692	pgsub(k)=0.0
693	pgdep(k)=0.0
694	pgmlt(k)=0.0
695	pgmltevp(k)=0.0
696	qschg(k)=0.
697	qgchg(k)=0.
698	enddo
699	!
700	!
701	! Set rs0=episp0*oprez(k)
702	! episp0=e*saturated pressure at 273.15 K
703	! e = 0.622
704	!
705	DO k=kts,kte
706	rs0(k)=ep21000.svp1/(prez(k)-1000.*svp1)
707	END DO
708	!
709	!***********************************************************************
710	! Calculate precipitation fluxes due to terminal velocities.
711	!***********************************************************************
712	!
713	!- Calculate termianl velocity (vt?) of precipitation q?z
714	!- Find maximum vt? to determine the small delta t
715	!
716	!-- rain
717	!
718	t_del_tv=0.
719	del_tv=dtb
720	notlast=.true.
721	DO while (notlast)
722	!
723	min_q=kte
724	max_q=kts-1
725	!
726	do k=kts,kte-1
727	if (qrz(k) .gt. 1.0e-8) then
728	min_q=min0(min_q,k)
729	max_q=max0(max_q,k)
730	tmp1=sqrt(pirhowaterxnor/rho(k)/qrz(k))
731	tmp1=sqrt(tmp1)
732	vtrold(k)=o6constagambp4sqrho(k)/tmp1*constb
733	if (k .eq. 1) then
734	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
735	else
736	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
737	endif
738	else
739	vtrold(k)=0.
740	endif
741	enddo
742
743	if (max_q .ge. min_q) then
744	!
745	!- Check if the summation of the small delta t >= big delta t
746	! (t_del_tv) (del_tv) (dtb)
747
748	t_del_tv=t_del_tv+del_tv
749	!
750	if ( t_del_tv .ge. dtb ) then
751	notlast=.false.
752	del_tv=dtb+del_tv-t_del_tv
753	endif
754	!
755	fluxin=0.
756	do k=max_q,min_q,-1
757	fluxout=rho(k)vtrold(k)qrz(k)
758	flux=(fluxin-fluxout)/rho(k)/dzw(k)
759	tmpqrz=qrz(k)
760	qrz(k)=qrz(k)+del_tv*flux
761	fluxin=fluxout
762	enddo
763	if (min_q .eq. 1) then
764	pptrain=pptrain+fluxin*del_tv
765	else
766	qrz(min_q-1)=qrz(min_q-1)+del_tv* &
767	fluxin/rho(min_q-1)/dzw(min_q-1)
768	endif
769	!
770	else
771	notlast=.false.
772	endif
773	ENDDO
774
775	!
776	!-- snow
777	!
778	t_del_tv=0.
779	del_tv=dtb
780	notlast=.true.
781
782	DO while (notlast)
783	!
784	min_q=kte
785	max_q=kts-1
786	!
787	do k=kts,kte-1
788	if (qsz(k) .gt. 1.0e-8) then
789	min_q=min0(min_q,k)
790	max_q=max0(max_q,k)
791	tmp1=sqrt(pirhosnowxnos/rho(k)/qsz(k))
792	tmp1=sqrt(tmp1)
793	vtsold(k)=o6constcgamdp4sqrho(k)/tmp1*constd
794	if (k .eq. 1) then
795	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
796	else
797	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
798	endif
799	else
800	vtsold(k)=0.
801	endif
802	enddo
803
804	if (max_q .ge. min_q) then
805	!
806	!
807	!- Check if the summation of the small delta t >= big delta t
808	! (t_del_tv) (del_tv) (dtb)
809
810	t_del_tv=t_del_tv+del_tv
811
812	if ( t_del_tv .ge. dtb ) then
813	notlast=.false.
814	del_tv=dtb+del_tv-t_del_tv
815	endif
816	!
817	fluxin=0.
818	do k=max_q,min_q,-1
819	fluxout=rho(k)vtsold(k)qsz(k)
820	flux=(fluxin-fluxout)/rho(k)/dzw(k)
821	qsz(k)=qsz(k)+del_tv*flux
822	qsz(k)=amax1(0.,qsz(k))
823	fluxin=fluxout
824	enddo
825	if (min_q .eq. 1) then
826	pptsnow=pptsnow+fluxin*del_tv
827	else
828	qsz(min_q-1)=qsz(min_q-1)+del_tv* &
829	fluxin/rho(min_q-1)/dzw(min_q-1)
830	endif
831	!
832	else
833	notlast=.false.
834	endif
835
836	ENDDO
837	!
838	!-- grauupel
839	!
840	t_del_tv=0.
841	del_tv=dtb
842	notlast=.true.
843	!
844	DO while (notlast)
845	!
846	min_q=kte
847	max_q=kts-1
848	!
849	do k=kts,kte-1
850	if (qgz(k) .gt. 1.0e-8) then
851	min_q=min0(min_q,k)
852	max_q=max0(max_q,k)
853	tmp1=sqrt(pirhograulxnog/rho(k)/qgz(k))
854	tmp1=sqrt(tmp1)
855	term0=sqrt(4.gravrhograul*0.33334/rho(k)/cdrag)
856	vtgold(k)=o6gam4pt5term0*sqrt(1./tmp1)
857	if (k .eq. 1) then
858	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtgold(k))
859	else
860	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtgold(k))
861	endif
862	else
863	vtgold(k)=0.
864	endif
865	enddo
866
867	if (max_q .ge. min_q) then
868	!
869	!
870	!- Check if the summation of the small delta t >= big delta t
871	! (t_del_tv) (del_tv) (dtb)
872
873	t_del_tv=t_del_tv+del_tv
874
875	if ( t_del_tv .ge. dtb ) then
876	notlast=.false.
877	del_tv=dtb+del_tv-t_del_tv
878	endif
879
880	!
881	fluxin=0.
882	do k=max_q,min_q,-1
883	fluxout=rho(k)vtgold(k)qgz(k)
884	flux=(fluxin-fluxout)/rho(k)/dzw(k)
885	qgz(k)=qgz(k)+del_tv*flux
886	qgz(k)=amax1(0.,qgz(k))
887	fluxin=fluxout
888	enddo
889	if (min_q .eq. 1) then
890	pptgraul=pptgraul+fluxin*del_tv
891	else
892	qgz(min_q-1)=qgz(min_q-1)+del_tv* &
893	fluxin/rho(min_q-1)/dzw(min_q-1)
894	endif
895	!
896	else
897	notlast=.false.
898	endif
899	!
900	ENDDO
901
902	!
903	!-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL
904	!
905	t_del_tv=0.
906	del_tv=dtb
907	notlast=.true.
908	!
909	DO while (notlast)
910	!
911	min_q=kte
912	max_q=kts-1
913	!
914	do k=kts,kte-1
915	if (qiz(k) .gt. 1.0e-8) then
916	min_q=min0(min_q,k)
917	max_q=max0(max_q,k)
918	vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16 ! Heymsfield and Donner
919	if (k .eq. 1) then
920	del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
921	else
922	del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
923	endif
924	else
925	vtiold(k)=0.
926	endif
927	enddo
928
929	if (max_q .ge. min_q) then
930	!
931	!
932	!- Check if the summation of the small delta t >= big delta t
933	! (t_del_tv) (del_tv) (dtb)
934
935	t_del_tv=t_del_tv+del_tv
936
937	if ( t_del_tv .ge. dtb ) then
938	notlast=.false.
939	del_tv=dtb+del_tv-t_del_tv
940	endif
941
942	fluxin=0.
943	do k=max_q,min_q,-1
944	fluxout=rho(k)vtiold(k)qiz(k)
945	flux=(fluxin-fluxout)/rho(k)/dzw(k)
946	qiz(k)=qiz(k)+del_tv*flux
947	qiz(k)=amax1(0.,qiz(k))
948	fluxin=fluxout
949	enddo
950	if (min_q .eq. 1) then
951	pptice=pptice+fluxin*del_tv
952	else
953	qiz(min_q-1)=qiz(min_q-1)+del_tv* &
954	fluxin/rho(min_q-1)/dzw(min_q-1)
955	endif
956	!
957	else
958	notlast=.false.
959	endif
960	!
961	ENDDO
962	do k=kts,kte-1 !sg beg
963	precrz(k)=rho(k)vtrold(k)qrz(k)
964	preciz(k)=rho(k)vtiold(k)qiz(k)
965	precsz(k)=rho(k)vtsold(k)qsz(k)
966	precgz(k)=rho(k)vtgold(k)qgz(k)
967	enddo !sg end
968	precrz(kte)=0. !wig - top level never set for vtXold vars
969	preciz(kte)=0. !wig
970	precsz(kte)=0. !wig
971	precgz(kte)=0. !wig
972
973
974	! Microphpysics processes
975	!
976	DO 2000 k=kts,kte
977	!
978	!***********************************************************************
979	!*** diagnose mixing ratios (qrz,qsz), terminal ***
980	!*** velocities (vtr,vts), and slope parameters in size ***
981	!*** distribution(xlambdar,xlambdas) of rain and snow ***
982	!*** follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7) ***
983	!***********************************************************************
984	!
985	!**** assuming no cloud water can exist in the top two levels due to
986	!**** radiation consideration
987	!
988	!! if
989	!! unsaturated,
990	!! no cloud water, rain, ice, snow and graupel
991	!! then
992	!! skip these processes and jump to line 2000
993	!
994	!
995	tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)+qgz(k)*gindex
996	if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k) &
997	.and. tmp .eq. 0.0 ) go to 2000
998
999	!! calculate terminal velocity of rain
1000	!
1001	if (qrz(k) .gt. 1.0e-8) then
1002	tmp1=sqrt(pirhowaterxnor*orho(k)/qrz(k))
1003	xlambdar(k)=sqrt(tmp1)
1004	olambdar(k)=1.0/xlambdar(k)
1005	vtrold(k)=o6constagambp4sqrho(k)olambdar(k)**constb
1006	else
1007	vtrold(k)=0.
1008	olambdar(k)=0.
1009	endif
1010	!
1011	! if (qrz(k) .gt. 1.0e-12) then
1012	if (qrz(k) .gt. 1.0e-8) then
1013	tmp1=sqrt(pirhowaterxnor*orho(k)/qrz(k))
1014	xlambdar(k)=sqrt(tmp1)
1015	olambdar(k)=1.0/xlambdar(k)
1016	vtr(k)=o6constagambp4sqrho(k)olambdar(k)**constb
1017	else
1018	vtr(k)=0.
1019	olambdar(k)=0.
1020	endif
1021	!
1022	!! calculate terminal velocity of snow
1023	!
1024	if (qsz(k) .gt. 1.0e-8) then
1025	tmp1=sqrt(pirhosnowxnos*orho(k)/qsz(k))
1026	xlambdas(k)=sqrt(tmp1)
1027	olambdas(k)=1.0/xlambdas(k)
1028	vtsold(k)=o6constcgamdp4sqrho(k)olambdas(k)**constd
1029	else
1030	vtsold(k)=0.
1031	olambdas(k)=0.
1032	endif
1033	!
1034	! if (qsz(k) .gt. 1.0e-12) then
1035	if (qsz(k) .gt. 1.0e-8) then
1036	tmp1=sqrt(pirhosnowxnos*orho(k)/qsz(k))
1037	xlambdas(k)=sqrt(tmp1)
1038	olambdas(k)=1.0/xlambdas(k)
1039	vts(k)=o6constcgamdp4sqrho(k)olambdas(k)**constd
1040	else
1041	vts(k)=0.
1042	olambdas(k)=0.
1043	endif
1044	!
1045	!! calculate terminal velocity of graupel
1046	!
1047	if (qgz(k) .gt. 1.0e-8) then
1048	tmp1=sqrt( pirhograulxnog*orho(k)/qgz(k))
1049	xlambdag(k)=sqrt(tmp1)
1050	olambdag(k)=1.0/xlambdag(k)
1051	term0=sqrt(4.gravrhograul0.33334orho(k)*ocdrag)
1052	vtgold(k)=o6gam4pt5term0*sqrt(olambdag(k))
1053	else
1054	vtgold(k)=0.
1055	olambdag(k)=0.
1056	endif
1057	!
1058	! if (qgz(k) .gt. 1.0e-12) then
1059	if (qgz(k) .gt. 1.0e-8) then
1060	tmp1=sqrt( pirhograulxnog*orho(k)/qgz(k))
1061	xlambdag(k)=sqrt(tmp1)
1062	olambdag(k)=1.0/xlambdag(k)
1063	term0=sqrt(4.gravrhograul0.33334orho(k)*ocdrag)
1064	vtg(k)=o6gam4pt5term0*sqrt(olambdag(k))
1065	else
1066	vtg(k)=0.
1067	olambdag(k)=0.
1068	endif
1069	!
1070	!***********************************************************************
1071	!*** compute viscosity,difusivity,thermal conductivity, and ****
1072	!*** Schmidt number ****
1073	!***********************************************************************
1074	!c------------------------------------------------------------------
1075	!c viscmu: dynamic viscosity of air kg/m/s
1076	!c visc: kinematic viscosity of air = viscmu/rho (m2/s)
1077	!c avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
1078	!c viscmu=1.718e-5 kg/m/s in RH
1079	!c diffwv: Diffusivity of water vapor in air
1080	!c adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
1081	!c = 8.7602 (8.794 in MM5) gcm/s3
1082	!c diffwv(k)=2.26e-5 m2/s
1083	!c schmidt: Schmidt number=visc/diffwv
1084	!c xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
1085	!c xka(k)=2.43e-2 J/m/s/K in RH
1086	!c axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
1087	!c------------------------------------------------------------------
1088
1089	viscmu(k)=avisctem(k)*1.5/(tem(k)+120.0)
1090	visc(k)=viscmu(k)*orho(k)
1091	diffwv(k)=adiffwvtem(k)1.81oprez(k)
1092	schmidt(k)=visc(k)/diffwv(k)
1093	xka(k)=axka*viscmu(k)
1094
1095	if (tem(k) .lt. 273.15) then
1096
1097	!
1098	!***********************************************************************
1099	!******* snow production processes for T < 0 C ********
1100	!***********************************************************************
1101	!c
1102	!c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
1103	!c! psaut=alpha1*(qi-qi0)
1104	!c! alpha1=1.0e-3exp(0.025(T-T0))
1105	!c
1106	! alpha1=1.0e-3exp( 0.025temcc(k) )
1107
1108	alpha1=1.0e-3exp( 0.025temcc(k) )
1109	!
1110	if(temcc(k) .lt. -20.0) then
1111	tmp1=-7.6+4.0exp( -0.2443e-3(abs(temcc(k))-20)**2.455 )
1112	qic=1.0e-3exp(tmp1)orho(k)
1113	else
1114	qic=qi0
1115	end if
1116	!testing
1117	! tmp1=amax1( 0.0,alpha1*(qiz(k)-qic) )
1118	! psaut(k)=amin1( tmp1,qizodt(k) )
1119
1120	tmp1=odtb(qiz(k)-qic)(1.0-exp(-alpha1*dtb))
1121	psaut(k)=amax1( 0.0,tmp1 )
1122
1123	!c
1124	!c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
1125	!c this process only considered when -31 C < T < 0 C
1126	!c Lin (33) and Hsie (17)
1127	!c
1128	!c!
1129	!c! parama1 and parama2 functions must be user supplied
1130	!c!
1131
1132	! testing
1133	if( qlz(k) .gt. 1.0e-10 ) then
1134	temc1=amax1(-30.99,temcc(k))
1135	! print*,'temc1',temc1,qlz(k)
1136	a1=parama1( temc1 )
1137	a2=parama2( temc1 )
1138	tmp1=1.0-a2
1139	!! change unit from cgs to mks
1140	a1=a10.001*tmp1
1141	!c! dtberg is the time needed for a crystal to grow from 40 to 50 um
1142	!c ! odtberg=1.0/dtberg
1143	odtberg=(a1tmp1)/(xmi50tmp1-xmi40*tmp1)
1144	!
1145	!c! compute terminal velocity of a 50 micron ice cystal
1146	!
1147	vti50=constcdi50constdsqrho(k)
1148	!
1149	eiw=1.0
1150	save1=a1xmi50*a2
1151	save2=0.25pieiwrho(k)di50di50vti50
1152	!
1153	tmp2=( save1 + save2*qlz(k) )
1154	!
1155	!! maximum number of 50 micron crystals limited by the amount
1156	!! of supercool water
1157	!
1158	xni50mx=qlzodt(k)/tmp2
1159	!
1160	!! number of 50 micron crystals produced
1161	!
1162	!
1163	xni50=qiz(k)( 1.0-exp(-dtbodtberg) )/xmi50
1164	xni50=amin1(xni50,xni50mx)
1165	!
1166	tmp3=odtbtmp2/save2( 1.0-exp(-save2xni50dtb) )
1167	psfw(k)=amin1( tmp3,qlzodt(k) )
1168	!testing
1169	! psfw(k)=0.
1170
1171	!0915 if( temcc(k).gt.-30.99 ) then
1172	!0915 a1=parama1( temcc(k) )
1173	!0915 a2=parama2( temcc(k) )
1174	!0915 tmp1=1.0-a2
1175	!! change unit from cgs to mks
1176	!0915 a1=a10.001*tmp1
1177
1178	!c! dtberg is the time needed for a crystal to grow from 40 to 50 um
1179	!c! odtberg=1.0/dtberg
1180	!0915 odtberg=(a1tmp1)/(xmi50tmp1-xmi40*tmp1)
1181
1182	!c! number of 50 micron crystals produced
1183	!0915 xni50=qiz(k)dtbodtberg/xmi50
1184
1185	!c! need to calculate the terminal velocity of a 50 micron
1186	!c! ice cystal
1187	!0915 vti50=constcdi50constdsqrho(k)
1188	!0915 eiw=1.0
1189	!0915 tmp2=xni50( a1xmi50**a2 + &
1190	!0915 0.25qlz(k)pieiwrho(k)di50di50*vti50 )
1191	!0915 psfw(k)=amin1( tmp2,qlzodt(k) )
1192	!0915 psfw(k)=0.
1193	!c
1194	!c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
1195	!c this process only considered when -31 C < T < 0 C
1196	!c
1197	tmp1=xni50*xmi50-psfw(k)
1198	psfi(k)=amin1(tmp1,qizodt(k))
1199	! testing
1200	! psfi(k)=0.
1201	end if
1202	!
1203
1204	!0915 tmp1=qiz(k)*odtberg
1205	!0915 psfi(k)=amin1(tmp1,qizodt(k))
1206	! testing
1207	!0915 psfi(k)=0.
1208	!0915 end if
1209	!
1210	if(qrz(k) .le. 0.0) go to 1000
1211	!
1212	! Processes (4) and (5) only need when qrz > 0.0
1213	!
1214	!c
1215	!c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
1216	!c may produce snow or graupel
1217	!c
1218	eri=1.0
1219	!0915 tmp1=qiz(k)pio4erixnorconsta*sqrho(k)
1220	!0915 tmp2=tmp1gambp3olambdar(k)**bp3
1221	!0915 praci(k)=amin1( tmp2,qizodt(k) )
1222
1223	save1=pio4erixnorconstasqrho(k)
1224	tmp1=save1gambp3olambdar(k)**bp3
1225	praci(k)=qizodt(k)( 1.0-exp(-tmp1dtb) )
1226
1227	!c
1228	!c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
1229	!c
1230	!0915 tmp2=tmp1rho(k)pio6rhowatergambp6oxmi &
1231	!0915 olambdar(k)**bp6
1232	!0915 piacr(k)=amin1( tmp2,qrzodt(k) )
1233
1234	tmp2=qiz(k)save1rho(k)pio6rhowatergambp6oxmi* &
1235	olambdar(k)**bp6
1236	piacr(k)=amin1( tmp2,qrzodt(k) )
1237
1238	!
1239	1000 continue
1240	!
1241	if(qsz(k) .le. 0.0) go to 1200
1242	!
1243	! Compute the following processes only when qsz > 0.0
1244	!
1245	!c
1246	!c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
1247	!c
1248	esi=exp( 0.025*temcc(k) )
1249	save1=pio4xnosconstcgamdp3sqrho(k)* &
1250	olambdas(k)**dp3
1251	tmp1=esi*save1
1252	psaci(k)=qizodt(k)( 1.0-exp(-tmp1dtb) )
1253
1254	!0915 tmp1=pio4xnosconstcgamdp3sqrho(k)* &
1255	!0915 olambdas(k)**dp3
1256	!0915 tmp2=qiz(k)esitmp1
1257	!0915 psaci(k)=amin1( tmp2,qizodt(k) )
1258	!c
1259	!c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1260	!c
1261	esw=1.0
1262	tmp1=esw*save1
1263	psacw(k)=qlzodt(K)( 1.0-exp(-tmp1dtb) )
1264
1265	!0915 tmp2=qlz(k)eswtmp1
1266	!0915 psacw(k)=amin1( tmp2,qlzodt(k) )
1267	!c
1268	!c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
1269	!c includes consideration of ventilation effect
1270	!c
1271	!c abi=2pi(Si-1)/rho/(A"+B")
1272	!c
1273	tmpa=rvaporxka(k)tem(k)*tem(k)
1274	tmpb=xlsxlsrho(k)qsiz(k)diffwv(k)
1275	tmpc=tmpaqsiz(k)diffwv(k)
1276	abi=2.0pi(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1277	!
1278	!c vf1s,vf2s=ventilation factors for snow
1279	!c vf1s=0.78,vf2s=0.31 in LIN
1280	!
1281	tmp1=constcsqrho(k)olambdas(k)**dp5/visc(k)
1282	tmp2=abixnos( vf1solambdas(k)olambdas(k)+ &
1283	vf2sschmidt(k)0.33334gamdp5o2*sqrt(tmp1) )
1284	tmp3=odtb*( qvz(k)-qsiz(k) )
1285	!
1286	if( tmp2 .le. 0.0) then
1287	tmp2=amax1( tmp2,tmp3)
1288	pssub(k)=amax1( tmp2,-qszodt(k) )
1289	psdep(k)=0.0
1290	else
1291	psdep(k)=amin1( tmp2,tmp3 )
1292	pssub(k)=0.0
1293	end if
1294
1295	!0915 psdep(k)=amax1(0.0,tmp2)
1296	!0915 pssub(k)=amin1(0.0,tmp2)
1297	!0915 pssub(k)=amax1( pssub(k),-qszodt(k) )
1298	!
1299	if(qrz(k) .le. 0.0) go to 1200
1300	!
1301	! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
1302	!
1303	!c
1304	!c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
1305	!c
1306	esr=1.0
1307	tmpa=olambdar(k)*olambdar(k)
1308	tmpb=olambdas(k)*olambdas(k)
1309	tmpc=olambdar(k)*olambdas(k)
1310	tmp1=pipiesrxnorxnosabs( vtr(k)-vts(k) )orho(k)
1311	tmp2=tmpbtmpbolambdar(k)(5.0tmpb+2.0tmpc+0.5tmpa)
1312	tmp3=tmp1rhosnowtmp2
1313	pracs(k)=amin1( tmp3,qszodt(k) )
1314	!c
1315	!c (10) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1316	!c
1317	tmp3=tmpatmpaolambdas(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1318	tmp4=tmp1rhowatertmp3
1319	psacr(k)=amin1( tmp4,qrzodt(k) )
1320	!
1321	1200 continue
1322	!
1323	else
1324	!
1325	!***********************************************************************
1326	!******* snow production processes for T > 0 C ********
1327	!***********************************************************************
1328	!
1329	if (qsz(k) .le. 0.0) go to 1400
1330	!c
1331	!c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1332	!c
1333	esw=1.0
1334
1335	tmp1=eswpio4xnosconstcgamdp3sqrho(k) &
1336	olambdas(k)**dp3
1337	psacw(k)=qlzodt(k)( 1.0-exp(-tmp1dtb) )
1338
1339	!0915 tmp1=pio4xnosconstcgamdp3sqrho(k)* &
1340	!0915 olambdas(k)**dp3
1341	!0915 tmp2=qlz(k)eswtmp1
1342	!0915 psacw(k)=amin1( tmp2,qlzodt(k) )
1343	!c
1344	!c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1345	!c
1346	esr=1.0
1347	tmpa=olambdar(k)*olambdar(k)
1348	tmpb=olambdas(k)*olambdas(k)
1349	tmpc=olambdar(k)*olambdas(k)
1350	tmp1=pipiesrxnorxnosabs( vtr(k)-vts(k) )orho(k)
1351	tmp2=tmpatmpaolambdas(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1352	tmp3=tmp1rhowatertmp2
1353	psacr(k)=amin1( tmp3,qrzodt(k) )
1354	!c
1355	!c (3) MELTING OF SNOW (Psmlt): Lin (32)
1356	!c Psmlt is negative value
1357	!
1358	delrs=rs0(k)-qvz(k)
1359	term1=2.0piorho(k)( xlvdiffwv(k)rho(k)delrs- &
1360	xka(k)*temcc(k) )
1361	tmp1=constcsqrho(k)olambdas(k)**dp5/visc(k)
1362	tmp2=xnos( vf1solambdas(k)*olambdas(k)+ &
1363	vf2sschmidt(k)0.33334gamdp5o2*sqrt(tmp1) )
1364	tmp3=term1oxlftmp2-cwoxlftemcc(k)( psacw(k)+psacr(k) )
1365	tmp4=amin1(0.0,tmp3)
1366	psmlt(k)=amax1( tmp4,-qszodt(k) )
1367	!c
1368	!c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
1369	!c but use Lin et al. coefficience
1370	!c Psmltevp is a negative value
1371	!c
1372	tmpa=rvaporxka(k)tem(k)*tem(k)
1373	tmpb=xlvxlvrho(k)qswz(k)diffwv(k)
1374	tmpc=tmpaqswz(k)diffwv(k)
1375	tmpd=amin1( 0.0,(qvoqswz(k)-0.90)qswz(k)odtb )
1376
1377	! abr=2.0pi(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1378
1379	abr=2.0pi(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1380	!
1381	!**** allow evaporation to occur when RH less than 90%
1382	!**** here not using 100% because the evaporation cooling
1383	!**** of temperature is not taking into account yet; hence,
1384	!**** the qsw value is a little bit larger. This will avoid
1385	!**** evaporation can generate cloud.
1386	!
1387	!c vf1s,vf2s=ventilation factors for snow
1388	!c vf1s=0.78,vf2s=0.31 in LIN
1389	!
1390	tmp1=constcsqrho(k)olambdas(k)**dp5/visc(k)
1391	tmp2=abrxnos( vf1solambdas(k)olambdas(k)+ &
1392	vf2sschmidt(k)0.33334gamdp5o2*sqrt(tmp1) )
1393	tmp3=amin1(0.0,tmp2)
1394	tmp3=amax1( tmp3,tmpd )
1395	psmltevp(k)=amax1( tmp3,-qszodt(k) )
1396	1400 continue
1397	!
1398	end if
1399
1400	!***********************************************************************
1401	!******* rain production processes ********
1402	!***********************************************************************
1403	!
1404	!c
1405	!c (1) AUTOCONVERSION OF RAIN (Praut): RH
1406	!sg: begin
1407	if(flag_qndrop)then
1408	if( qndropz(k) >= 1. ) then
1409	! Liu et al. autoconversion scheme
1410	rhocgs=rho(k)*1.e-3
1411	liqconc=rhocgs*qlz(k)
1412	capn=rhocgs*qndropz(k)
1413	! rate function
1414	if(liqconc.gt.1.e-10)then
1415	p0=kappabeta/capn(liqconcliqconcliqconc)
1416	xc=9.7d-17capnsqrt(capn)/(liqconc*liqconc)
1417	! Calculate autoconversion rate (g/g/s)
1418	if(xc.lt.10.)then
1419	praut(k)=p0/rhocgs0.5d0(xcxc+2xc+2.0d0)* &
1420	(1.0d0+xc)dexp(-2.0d0xc)
1421	endif
1422	endif
1423	endif
1424	else
1425	!sg: end
1426	!c araut=afa*rho
1427	!c afa=0.001 Rate coefficient for autoconvergence
1428	!c
1429	!c araut=1.0e-3
1430	!c
1431	araut=0.001
1432	!testing
1433	! tmp1=amax1( 0.0,araut*(qlz(k)-ql0) )
1434	! praut(k)=amin1( tmp1,qlzodt(k) )
1435	tmp1=odtb(qlz(k)-ql0)( 1.0-exp(-araut*dtb) )
1436	praut(k)=amax1( 0.0,tmp1 )
1437	endif !sg
1438
1439	!c
1440	!c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
1441	!c
1442	erw=1.0
1443	! tmp1=qlz(k)pio4erwxnorconsta*sqrho(k)
1444	! tmp2=tmp1gambp3olambdar(k)**bp3
1445	! pracw(k)=amin1( tmp2,qlzodt(k) )
1446
1447	tmp1=pio4erwxnorconstasqrho(k)* &
1448	gambp3olambdar(k)*bp3
1449	pracw(k)=qlzodt(k)( 1.0-exp(-tmp1dtb) )
1450
1451	!c
1452	!c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
1453	!c Prevp is negative value
1454	!c
1455	!c Sw=qvoqsw : saturation ratio
1456	!c
1457	tmpa=rvaporxka(k)tem(k)*tem(k)
1458	tmpb=xlvxlvrho(k)qswz(k)diffwv(k)
1459	tmpc=tmpaqswz(k)diffwv(k)
1460	tmpd=amin1(0.0,(qvoqswz(k)-0.90)qswz(k)odtb)
1461	!
1462	! abr=2.0pi(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1463
1464	abr=2.0pi(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1465	!
1466	!c vf1r,vf2r=ventilation factors for rain
1467	!c vf1r=0.78,vf2r=0.31 in RH, LIN and MM5
1468	!
1469	vf1r=0.78
1470	vf2r=0.31
1471	tmp1=constasqrho(k)olambdar(k)**bp5/visc(k)
1472	tmp2=abrxnor( vf1rolambdar(k)olambdar(k)+ &
1473	vf2rschmidt(k)0.33334gambp5o2*sqrt(tmp1) )
1474	tmp3=amin1( 0.0,tmp2 )
1475	tmp3=amax1( tmp3,tmpd )
1476	prevp(k)=amax1( tmp3,-qrzodt(k) )
1477
1478	!
1479	! if(iout .gt. 0) write(20,*)'tmp1,tmp2,tmp3=',tmp1,tmp2,tmp3
1480	! if(iout .gt. 0) write(20,*)'qlz,qiz,qrz=',qlz(k),qiz(k),qrz(k)
1481	! if(iout .gt. 0) write(20,*)'tem,qsz,qvz=',tem(k),qsz(k),qvz(k)
1482
1483
1484
1485	! if (gindex .eq. 0.) goto 900
1486	!
1487	if (tem(k) .lt. 273.15) then
1488	!
1489	!
1490	!-- graupel
1491	!***********************************************************************
1492	!******* graupel production processes for T < 0 C ********
1493	!***********************************************************************
1494	!c
1495	!c (1) AUTOCONVERSION OF SNOW TO FORM GRAUPEL (Pgaut): Lin (37)
1496	!c pgaut=alpha2*(qsz-qs0)
1497	!c qs0=6.0E-4
1498	!c alpha2=1.0e-3exp(0.09temcc(k)) Lin (38)
1499	!
1500	alpha2=1.0e-3exp(0.09temcc(k))
1501	!
1502
1503	! testing
1504	! tmp1=alpha2*(qsz(k)-qs0)
1505	! tmp1=amax1(0.0,tmp1)
1506	! pgaut(k)=amin1( tmp1,qszodt(k) )
1507
1508	tmp1=odtb(qsz(k)-qs0)(1.0-exp(-alpha2*dtb))
1509	pgaut(k)=amax1( 0.0,tmp1 )
1510
1511	!c
1512	!c (2) FREEZING OF RAIN TO FORM GRAUPEL (Pgfr): Lin (45)
1513	!c positive value
1514	!c Constant in Bigg freezing Aplume=Ap=0.66 /k
1515	!c Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
1516	!
1517
1518	if (qrz(k) .gt. 1.e-8 ) then
1519	Bp=100.
1520	Ap=0.66
1521	tmp1=olambdar(k)olambdar(k)olambdar(k)
1522	tmp2=20.pipiBpxnorrhowaterorho(k)* &
1523	(exp(-Aptemcc(k))-1.0)tmp1tmp1olambdar(k)
1524	Pgfr(k)=amin1( tmp2,qrzodt(k) )
1525	else
1526	Pgfr(k)=0
1527	endif
1528
1529	!c
1530	!c if (qgz(k) = 0.0) skip the other step below about graupel
1531	!c
1532	if (qgz(k) .eq. 0.0) goto 4000
1533
1534	!c
1535	!c Comparing Pgwet(wet process) and Pdry(dry process),
1536	!c we will pick up the small one.
1537	!c
1538
1539	!c ---------------
1540	!c \| dry processes \|
1541	!c ---------------
1542	!c
1543	!c (3) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1544	!c egw=1.0
1545	!c Cdrag=0.6 drag coefficients for hairstone
1546	!c constg=sqrt(4.gravrhograul0.33334orho(k)/Cdrag)
1547	!c
1548	egw=1.0
1549	constg=sqrt(4.gravrhograul0.33334orho(k)*oCdrag)
1550	tmp1=pio4xnoggam3pt5constgolambdag(k)**3.5
1551	tmp2=qlz(k)egwtmp1
1552	Pgacw(k)=amin1( tmp2,qlzodt(k) )
1553	!c
1554	!c (4) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgaci): Lin (41)
1555	!c egi=1. for wet growth
1556	!c egi=0.1 for dry growth
1557	!c
1558	egi=0.1
1559	tmp2=qiz(k)egitmp1
1560	pgaci(k)=amin1( tmp2,qizodt(k) )
1561
1562
1563	!c
1564	!c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1565	!c Compute processes (6) only when qsz > 0.0 and qgz > 0.0
1566	!c
1567	egs=exp(0.09*temcc(k))
1568	tmpa=olambdas(k)*olambdas(k)
1569	tmpb=olambdag(k)*olambdag(k)
1570	tmpc=olambdas(k)*olambdag(k)
1571	tmp1=pipixnosxnogabs( vts(k)-vtg(k) )*orho(k)
1572	tmp2=tmpatmpaolambdag(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1573	tmp3=tmp1egsrhosnow*tmp2
1574	Pgacs(k)=amin1( tmp3,qszodt(k) )
1575
1576
1577	!c
1578	!c (6) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1579	!c Compute processes (6) only when qrz > 0.0 and qgz > 0.0
1580	!c egr=1.
1581	!c
1582	egr=1.
1583	tmpa=olambdar(k)*olambdar(k)
1584	tmpb=olambdag(k)*olambdag(k)
1585	tmpc=olambdar(k)*olambdag(k)
1586	tmp1=pipixnorxnogabs( vtr(k)-vtg(k) )*orho(k)
1587	tmp2=tmpatmpaolambdag(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1588	tmp3=tmp1egrrhowater*tmp2
1589	pgacr(k)=amin1( tmp3,qrzodt(k) )
1590
1591	!c
1592	!c (7) Calculate total dry process effect Pdry(k)
1593	!c
1594	Pdry(k)=Pgacw(k)+pgaci(k)+Pgacs(k)+pgacr(k)
1595
1596	!c ---------------
1597	!c \| wet processes \|
1598	!c ---------------
1599	!c
1600	!c (3) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgacip): Lin (41)
1601	!c egi=1. for wet growth
1602	!c egi=0.1 for dry growth
1603	!c
1604	tmp2=10.*pgaci(k)
1605	pgacip(k)=amin1( tmp2,qizodt(k) )
1606
1607	!c
1608	!c (4) ACCRETION OF SNOW BY GRAUPEL ((Pgacsp) : Lin (29)
1609	!c Compute processes (6) only when qsz > 0.0 and qgz > 0.0
1610	!c egs=exp(0.09*(tem(k)-273.15)) when T < 273.15 k
1611	!c
1612	tmp3=Pgacs(k)*1.0/egs
1613	Pgacsp(k)=amin1( tmp3,qszodt(k) )
1614
1615	!c
1616	!c (5) WET GROWTH OF GRAUPEL (Pgwet) : Lin (43)
1617	!c may involve Pgacs or Pgaci and
1618	!c must include PPgacw or Pgacr, or both.
1619	!c ( The amount of Pgacw which is not able
1620	!c to freeze is shed to rain. )
1621	IF(temcc(k).gt.-40.)THEN
1622
1623	term0=constgolambdag(k)*5.5/visc(k)
1624
1625	!c
1626	!c vf1s,vf2s=ventilation factors for graupel
1627	!c vf1s=0.78,vf2s=0.31 in LIN
1628	!c Cdrag=0.6 drag coefficient for hairstone
1629	!c constg2=vf1solambdag(k)olambdag(k)+
1630	!c vf2sschmidt(k)0.33334gam2pt75*sqrt(term0)
1631
1632	delrs=rs0(k)-qvz(k)
1633	tmp0=1./(xlf+cw*temcc(k))
1634	tmp1=2.pixnog(rho(k)xlvdiffwv(k)delrs-xka(k)* &
1635	temcc(k))orho(k)tmp0
1636	constg2=vf1solambdag(k)olambdag(k)+ &
1637	vf2sschmidt(k)0.33334gam2pt75*sqrt(term0)
1638	tmp3=tmp1constg2+(Pgacip(k)+Pgacsp(k)) &
1639	(1-Citemcc(k)tmp0)
1640	tmp3=amax1(0.0,tmp3)
1641	Pgwet(k)=amax1(tmp3,qlzodt(k)+qszodt(k)+qizodt(k) )
1642
1643	!c
1644	!c Comparing Pgwet(wet process) and Pdry(dry process),
1645	!c we will apply the small one.
1646	!c if dry processes then delta4=1.0
1647	!c if wet processes then delta4=0.0
1648	!
1649	if ( Pdry(k) .lt. Pgwet(k) ) then
1650	delta4=1.0
1651	else
1652	delta4=0.0
1653	endif
1654	ELSE
1655	delta4=1.0
1656	ENDIF
1657
1658	!c
1659	!c
1660	!c (6) Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
1661	!c if Pgacrp(k) > 0. then some of the rain is frozen to hail
1662	!c if Pgacrp(k) < 0. then some of the cloud water collected
1663	!c by the hail is unable to freeze and is
1664	!c shed as rain.
1665	!c
1666	Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
1667
1668	!c
1669	!c (8) DEPOSITION/SUBLIMATION OF GRAUPEL (Pgdep/Pgsub): Lin (46)
1670	!c includes ventilation effect
1671	!c constg=sqrt(4.gravrhograul0.33334orho(k)/Cdrag)
1672	!c constg2=vf1solambdag(k)olambdag(k)+
1673	!c vf2sschmidt(k)0.33334gam2pt75*constg
1674	!c
1675	!c abg=2pi(Si-1)/rho/(A"+B")
1676	!c
1677	tmpa=rvaporxka(k)tem(k)*tem(k)
1678	tmpb=xlsxlsrho(k)qsiz(k)diffwv(k)
1679	tmpc=tmpaqsiz(k)diffwv(k)
1680	abg=2.0pi(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1681	!c
1682	!c vf1s,vf2s=ventilation factors for graupel
1683	!c vf1s=0.78,vf2s=0.31 in LIN
1684	!c Cdrag=0.6 drag coefficient for hairstone
1685	!c
1686	term0=constgolambdag(k)*5.5/visc(k)
1687	constg2=vf1solambdag(k)olambdag(k)+ &
1688	vf2sschmidt(k)0.33334gam2pt75*sqrt(term0)
1689	tmp2=abgxnogconstg2
1690	pgdep(k)=amax1(0.0,tmp2)
1691	pgsub(k)=amin1(0.0,tmp2)
1692	pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
1693
1694	4000 continue
1695	else
1696	!
1697	!***********************************************************************
1698	!******* graupel production processes for T > 0 C ********
1699	!***********************************************************************
1700	!
1701	!c
1702	!c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1703	!c egw=1.0
1704	!c Cdrag=0.6 drag coefficients for hairstone
1705	!c constg=sqrt(4.gravrhograul0.33334orho(k)/Cdrag)
1706
1707	egw=1.0
1708	constg=sqrt(4.gravrhograul0.33334orho(k)*oCdrag)
1709	tmp1=pio4xnoggam3pt5constgolambdag(k)**3.5
1710	tmp2=qlz(k)egwtmp1
1711	Pgacw(k)=amin1( tmp2,qlzodt(k) )
1712
1713	!c
1714	!c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1715	!c Compute processes (5) only when qrz > 0.0 and qgz > 0.0
1716	!c egr=1.
1717	!c
1718	egr=1.
1719	tmpa=olambdar(k)*olambdar(k)
1720	tmpb=olambdag(k)*olambdag(k)
1721	tmpc=olambdar(k)*olambdag(k)
1722	tmp1=pipixnorxnogabs( vtr(k)-vtg(k) )*orho(k)
1723	tmp2=tmpatmpaolambdag(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1724	tmp3=tmp1egrrhowater*tmp2
1725	pgacr(k)=amin1( tmp3,qrzodt(k) )
1726
1727
1728	!c
1729	!c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
1730	!c Pgmlt is negative value
1731	!c constg=sqrt(4.gravrhograul0.33334orho(k)/Cdrag)
1732	!c constg2=vf1solambdag(k)olambdag(k)+
1733	!c vf2sschmidt(k)0.33334gam2pt75*constg
1734	!c Cdrag=0.6 drag coefficients for hairstone
1735	!
1736	delrs=rs0(k)-qvz(k)
1737	term1=2.0piorho(k)( xlvdiffwv(k)rho(k)delrs- &
1738	xka(k)*temcc(k) )
1739	term0=sqrt(4.gravrhograul0.33334orho(k)*ocdrag) &
1740	olambdag(k)*5.5/visc(k)
1741
1742	constg2=vf1solambdag(k)olambdag(k)+ &
1743	vf2sschmidt(k)0.33334gam2pt75*sqrt(term0)
1744	tmp2=xnog*constg2
1745	tmp3=term1oxlftmp2-cwoxlftemcc(k)( pgacw(k)+pgacr(k) )
1746	tmp4=amin1(0.0,tmp3)
1747	pgmlt(k)=amax1( tmp4,-qgzodt(k) )
1748
1749
1750	!c
1751	!c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
1752	!c but use Lin et al. coefficience
1753	!c Pgmltevp is a negative value
1754	!c abg=2.0pi(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1755	!c
1756	tmpa=rvaporxka(k)tem(k)*tem(k)
1757	tmpb=xlvxlvrho(k)qswz(k)diffwv(k)
1758	tmpc=tmpaqswz(k)diffwv(k)
1759	tmpd=amin1( 0.0,(qvoqswz(k)-0.90)qswz(k)odtb )
1760
1761	!c
1762	!c abg=2pi(Si-1)/rho/(A"+B")
1763	!c
1764	abg=2.0pi(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1765	!
1766	!**** allow evaporation to occur when RH less than 90%
1767	!**** here not using 100% because the evaporation cooling
1768	!**** of temperature is not taking into account yet; hence,
1769	!**** the qgw value is a little bit larger. This will avoid
1770	!**** evaporation can generate cloud.
1771	!
1772	!c vf1s,vf2s=ventilation factors for snow
1773	!c vf1s=0.78,vf2s=0.31 in LIN
1774	!c constg=sqrt(4.gravrhograul0.33334orho(k)/Cdrag)
1775	!c constg2=vf1solambdag(k)olambdag(k)+
1776	!c vf2sschmidt(k)0.33334gam2pt75*constg
1777	!
1778	tmp2=abgxnogconstg2
1779	tmp3=amin1(0.0,tmp2)
1780	tmp3=amax1( tmp3,tmpd )
1781	pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
1782
1783	!c
1784	!c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1785	!c Compute processes (3) only when qsz > 0.0 and qgz > 0.0
1786	!c egs=1.0
1787	!c
1788	egs=1.
1789	tmpa=olambdas(k)*olambdas(k)
1790	tmpb=olambdag(k)*olambdag(k)
1791	tmpc=olambdas(k)*olambdag(k)
1792	tmp1=pipixnosxnogabs( vts(k)-vtg(k) )*orho(k)
1793	tmp2=tmpatmpaolambdag(k)(5.0tmpa+2.0tmpc+0.5tmpb)
1794	tmp3=tmp1egsrhosnow*tmp2
1795	Pgacs(k)=amin1( tmp3,qszodt(k) )
1796
1797	endif
1798
1799
1800	!
1801	900 continue
1802
1803	!cc
1804	!c
1805	!c**********************************************************************
1806	!c*** combine all processes together and avoid negative ***
1807	!c***** water substances
1808	!***********************************************************************
1809	!c
1810	if ( temcc(k) .lt. 0.0) then
1811	!,delta4,1.-delta4
1812	!c
1813	!c gdelta4=gindex*delta4
1814	!c g1sdelt4=gindex*(1.-delta4)
1815	!c
1816	gdelta4=gindex*delta4
1817	g1sdelt4=gindex*(1.-delta4)
1818	!c
1819	!c combined water vapor depletions
1820	!c
1821	!cc graupel
1822	tmp=psdep(k)+pgdep(k)*gindex
1823	if ( tmp .gt. qvzodt(k) ) then
1824	factor=qvzodt(k)/tmp
1825	psdep(k)=psdep(k)*factor
1826	pgdep(k)=pgdep(k)factorgindex
1827	end if
1828	!c
1829	!c combined cloud water depletions
1830	!c
1831	tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
1832	if ( tmp .gt. qlzodt(k) ) then
1833	factor=qlzodt(k)/tmp
1834	praut(k)=praut(k)*factor
1835	psacw(k)=psacw(k)*factor
1836	psfw(k)=psfw(k)*factor
1837	pracw(k)=pracw(k)*factor
1838	!cc graupel
1839	Pgacw(k)=Pgacw(k)factorgindex
1840	end if
1841	!c
1842	!c combined cloud ice depletions
1843	!c
1844	tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
1845	+Pgacip(k)*g1sdelt4
1846	if (tmp .gt. qizodt(k) ) then
1847	factor=qizodt(k)/tmp
1848	psaut(k)=psaut(k)*factor
1849	psaci(k)=psaci(k)*factor
1850	praci(k)=praci(k)*factor
1851	psfi(k)=psfi(k)*factor
1852	!cc graupel
1853	Pgaci(k)=Pgaci(k)factorgdelta4
1854	Pgacip(k)=Pgacip(k)factorg1sdelt4
1855	endif
1856	!c
1857	!c combined all rain processes
1858	!c
1859	tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1860	+Pgfr(k)gindex+Pgacr(k)gdelta4 &
1861	+Pgacrp(k)*g1sdelt4
1862	if (tmp_r .gt. qrzodt(k) ) then
1863	factor=qrzodt(k)/tmp_r
1864	piacr(k)=piacr(k)*factor
1865	psacr(k)=psacr(k)*factor
1866	prevp(k)=prevp(k)*factor
1867	!cc graupel
1868	Pgfr(k)=Pgfr(k)factorgindex
1869	Pgacr(k)=Pgacr(k)factorgdelta4
1870	Pgacrp(k)=Pgacrp(k)factorg1sdelt4
1871	endif
1872
1873	!c
1874	!c if qrz < 1.0E-4 and qsz < 1.0E-4 then delta2=1.
1875	!c (all Pracs and Psacr become to snow)
1876	!c if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
1877	!c (all Pracs and Psacr become to graupel)
1878	!c
1879	if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
1880	delta2=1.0
1881	else
1882	delta2=0.0
1883	endif
1884	!
1885	!cc graupel
1886
1887	!c
1888	!c if qrz(k) < 1.0e-4 then delta3=1. means praci(k) --> qs
1889	!c piacr(k) --> qs
1890	!c if qrz(k) > 1.0e-4 then delta3=0. means praci(k) --> qg
1891	!c piacr(k) --> qg : Lin (20)
1892
1893	if (qrz(k) .lt. 1.0e-4) then
1894	delta3=1.0
1895	else
1896	delta3=0.0
1897	endif
1898	!
1899	!c
1900	!c if gindex = 0.(no graupel) then delta2=1.0
1901	!c delta3=1.0
1902	!c
1903	if (gindex .eq. 0.) then
1904	delta2=1.0
1905	delta3=1.0
1906	endif
1907	!
1908	!c
1909	!c combined all snow processes
1910	!c
1911	tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1912	psfi(k)+praci(k)delta3+piacr(k)delta3+ &
1913	psdep(k))+Pgaut(k)gindex+Pgacs(k)gdelta4+ &
1914	Pgacsp(k)g1sdelt4+Pracs(k)(1.-delta2)- &
1915	Psacr(k)*delta2
1916	if ( tmp_s .gt. qszodt(k) ) then
1917	factor=qszodt(k)/tmp_s
1918	pssub(k)=pssub(k)*factor
1919	Pracs(k)=Pracs(k)*factor
1920	!cc graupel
1921	Pgaut(k)=Pgaut(k)factorgindex
1922	Pgacs(k)=Pgacs(k)factorgdelta4
1923	Pgacsp(k)=Pgacsp(k)factorg1sdelt4
1924	endif
1925
1926	!cc graupel
1927	!
1928
1929	! if (gindex .eq. 0.) goto 998
1930	!c
1931	!c combined all graupel processes
1932	!c
1933	tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)delta4-Pgaci(k)delta4 &
1934	-Pgacr(k)delta4-Pgacs(k)delta4 &
1935	-pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1936	-psacr(k)(1-delta2)-Pracs(k)(1-delta2) &
1937	-praci(k)(1-delta3)-piacr(k)(1-delta3)
1938	if (tmp_g .gt. qgzodt(k)) then
1939	factor=qgzodt(k)/tmp_g
1940	pgsub(k)=pgsub(k)*factor
1941	endif
1942
1943	998 continue
1944	!c
1945	!c calculate new water substances, thetae, tem, and qvsbar
1946	!c
1947
1948	!cc graupel
1949	pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
1950	-pgdep(k)*gindex
1951	qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
1952	pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
1953	if(flag_qndrop)then
1954	if( qlz(k) > 1e-20 ) &
1955	qndropz(k)=amax1( 0.0,qndropz(k)+dtbpclw(k)qndropz(k)/qlz(k) ) !sg
1956	endif
1957	qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1958	pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
1959	-Pgacip(k)*g1sdelt4
1960	qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1961	tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1962	+Pgfr(k)gindex+Pgacr(k)gdelta4 &
1963	+Pgacrp(k)*g1sdelt4
1964	232 format(i2,1x,6(f9.3,1x))
1965	prain(k)=-tmp_r
1966	qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1967	tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1968	psfi(k)+praci(k)delta3+piacr(k)delta3+ &
1969	psdep(k))+Pgaut(k)gindex+Pgacs(k)gdelta4+ &
1970	Pgacsp(k)g1sdelt4+Pracs(k)(1.-delta2)- &
1971	Psacr(k)*delta2
1972	psnow(k)=-tmp_s
1973	qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1974	qschg(k)=qschg(k)+psnow(k)
1975	qschg(k)=psnow(k)
1976	!cc graupel
1977	tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)delta4-Pgaci(k)delta4 &
1978	-Pgacr(k)delta4-Pgacs(k)delta4 &
1979	-pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1980	-psacr(k)(1-delta2)-Pracs(k)(1-delta2) &
1981	-praci(k)(1-delta3)-piacr(k)(1-delta3)
1982	252 format(i2,1x,6(f12.9,1x))
1983	262 format(i2,1x,7(f12.9,1x))
1984	pgraupel(k)=-tmp_g
1985	pgraupel(k)=pgraupel(k)*gindex
1986	qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
1987	! qgchg(k)=qgchg(k)+pgraupel(k)
1988	qgchg(k)=pgraupel(k)
1989	qgz(k)=qgz(k)*gindex
1990
1991	tmp=ocp/tothz(k)xLf(qschg(k)+qgchg(k))
1992	theiz(k)=theiz(k)+dtb*tmp
1993	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
1994	tem(k)=thz(k)*tothz(k)
1995
1996	temcc(k)=tem(k)-273.15
1997
1998	if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
1999	qlpqi=qlz(k)+qiz(k)
2000	if ( qlpqi .eq. 0.0 ) then
2001	qvsbar(k)=qsiz(k)
2002	else
2003	qvsbar(k)=( qiz(k)qsiz(k)+qlz(k)qswz(k) )/qlpqi
2004	endif
2005
2006	!
2007	else
2008	!c
2009	!c combined cloud water depletions
2010	!c
2011	tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
2012	if ( tmp .gt. qlzodt(k) ) then
2013	factor=qlzodt(k)/tmp
2014	praut(k)=praut(k)*factor
2015	psacw(k)=psacw(k)*factor
2016	pracw(k)=pracw(k)*factor
2017	!cc graupel
2018	pgacw(k)=pgacw(k)factorgindex
2019	end if
2020	!c
2021	!c combined all snow processes
2022	!c
2023	tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2024	if (tmp_s .gt. qszodt(k) ) then
2025	factor=qszodt(k)/tmp_s
2026	psmlt(k)=psmlt(k)*factor
2027	psmltevp(k)=psmltevp(k)*factor
2028	!cc graupel
2029	Pgacs(k)=Pgacs(k)factorgindex
2030	endif
2031
2032	!c
2033	!c
2034	!cc graupel
2035	!c
2036	! if (gindex .eq. 0.) goto 997
2037
2038	!c
2039	!c combined all graupel processes
2040	!c
2041	tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2042	if (tmp_g .gt. qgzodt(k)) then
2043	factor=qgzodt(k)/tmp_g
2044	pgmltevp(k)=pgmltevp(k)*factor
2045	pgmlt(k)=pgmlt(k)*factor
2046	endif
2047	!c
2048	997 continue
2049
2050	!c
2051	!c combined all rain processes
2052	!c
2053	tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2054	+pgmlt(k)gindex-pgacw(k)gindex
2055	if (tmp_r .gt. qrzodt(k) ) then
2056	factor=qrzodt(k)/tmp_r
2057	prevp(k)=prevp(k)*factor
2058	endif
2059	!c
2060	!c
2061	!c calculate new water substances and thetae
2062	!c
2063
2064
2065	pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex
2066	qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
2067	pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
2068	if(flag_qndrop)then
2069	if( qlz(k) > 1e-20 ) &
2070	qndropz(k)=amax1( 0.0,qndropz(k)+dtbpclw(k)qndropz(k)/qlz(k) ) !sg
2071	endif
2072	qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
2073	pcli(k)=0.0
2074	qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
2075	tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2076	+pgmlt(k)gindex-pgacw(k)gindex
2077	242 format(i2,1x,7(f9.6,1x))
2078	prain(k)=-tmp_r
2079	tmpqrz=qrz(k)
2080	qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
2081	tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2082	psnow(k)=-tmp_s
2083	qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
2084	! qschg(k)=qschg(k)+psnow(k)
2085	qschg(k)=psnow(k)
2086	!cc graupel
2087
2088	tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2089	! write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
2090	272 format(i2,1x,3(f12.9,1x))
2091	pgraupel(k)=-tmp_g*gindex
2092	qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
2093	! qgchg(k)=qgchg(k)+pgraupel(k)
2094	qgchg(k)=pgraupel(k)
2095	qgz(k)=qgz(k)*gindex
2096	!
2097	tmp=ocp/tothz(k)xLf(qschg(k)+qgchg(k))
2098	theiz(k)=theiz(k)+dtb*tmp
2099	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
2100
2101	tem(k)=thz(k)*tothz(k)
2102	temcc(k)=tem(k)-273.15
2103	! qswz(k)=episp0koprez(k) &
2104	! exp( svp2*temcc(k)/(tem(k)-svp3) )
2105	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
2106	qswz(k)=ep2*es/(prez(k)-es)
2107	qsiz(k)=qswz(k)
2108	qvsbar(k)=qswz(k)
2109	!
2110	end if
2111	preclw(k)=pclw(k) !sg
2112
2113	!
2114	!***********************************************************************
2115	!******** saturation adjustment ********
2116	!***********************************************************************
2117	!
2118	! allow supersaturation exits linearly from 0% at 500 mb to 50%
2119	! above 300 mb
2120	! 5.0e-5=1.0/(500mb-300mb)
2121	!
2122	rsat=1.0+0.5(50000.0-prez(k))5.0e-5
2123	rsat=amax1(1.0,rsat)
2124	rsat=amin1(1.5,rsat)
2125	rsat=1.0
2126	if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
2127
2128	!c
2129	!c unsaturated
2130	!c
2131	qvz(k)=qvz(k)+qlz(k)+qiz(k)
2132	qlz(k)=0.0
2133	qiz(k)=0.0
2134
2135	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
2136	tem(k)=thz(k)*tothz(k)
2137	temcc(k)=tem(k)-273.15
2138
2139	go to 1800
2140	!
2141	else
2142	!c
2143	!c saturated
2144	!c
2145	!
2146	pladj(k)=qlz(k)
2147	piadj(k)=qiz(k)
2148	!
2149
2150	CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2151	k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
2152
2153	!
2154	pladj(k)=odtb*(qlz(k)-pladj(k))
2155	piadj(k)=odtb*(qiz(k)-piadj(k))
2156	!
2157	pclw(k)=pclw(k)+pladj(k)
2158	pcli(k)=pcli(k)+piadj(k)
2159	pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
2160	!
2161	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
2162	tem(k)=thz(k)*tothz(k)
2163
2164	temcc(k)=tem(k)-273.15
2165
2166	! qswz(k)=episp0koprez(k) &
2167	! exp( svp2*temcc(k)/(tem(k)-svp3) )
2168	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
2169	qswz(k)=ep2*es/(prez(k)-es)
2170	if (tem(k) .lt. 273.15 ) then
2171	! qsiz(k)=episp0koprez(k) &
2172	! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2173	es=1000.svp1exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2174	qsiz(k)=ep2*es/(prez(k)-es)
2175	if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2176	else
2177	qsiz(k)=qswz(k)
2178	endif
2179	qlpqi=qlz(k)+qiz(k)
2180	if ( qlpqi .eq. 0.0 ) then
2181	qvsbar(k)=qsiz(k)
2182	else
2183	qvsbar(k)=( qiz(k)qsiz(k)+qlz(k)qswz(k) )/qlpqi
2184	endif
2185
2186	end if
2187
2188	!
2189	!***********************************************************************
2190	!*** melting and freezing of cloud ice and cloud water ***
2191	!***********************************************************************
2192	qlpqi=qlz(k)+qiz(k)
2193	if(qlpqi .le. 0.0) go to 1800
2194	!
2195	!c
2196	!c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
2197	!c
2198	if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
2199	!c
2200	!c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
2201	!c
2202	if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
2203	!c
2204	!c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
2205	!c this process only considered when -31 C < T < 0 C
2206	!c
2207	if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
2208	!c!
2209	!c! parama1 and parama2 functions must be user supplied
2210	!c!
2211	a1=parama1( temcc(k) )
2212	a2=parama2( temcc(k) )
2213	!! change unit from cgs to mks
2214	a1=a10.001*(1.0-a2)
2215	xnin=xni0exp(-bnitemcc(k))
2216	pidw(k)=xninorho(k)(a1xmnin*a2)
2217	end if
2218	!
2219	pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
2220	pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
2221	qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
2222	qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
2223
2224	!
2225	CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2226	k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
2227
2228	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
2229	tem(k)=thz(k)*tothz(k)
2230
2231	temcc(k)=tem(k)-273.15
2232
2233	! qswz(k)=episp0koprez(k) &
2234	! exp( svp2*temcc(k)/(tem(k)-svp3) )
2235	es=1000.svp1exp( svp2*temcc(k)/(tem(k)-svp3) )
2236	qswz(k)=ep2*es/(prez(k)-es)
2237
2238	if (tem(k) .lt. 273.15 ) then
2239	! qsiz(k)=episp0koprez(k) &
2240	! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2241	es=1000.svp1exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2242	qsiz(k)=ep2*es/(prez(k)-es)
2243	if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2244	else
2245	qsiz(k)=qswz(k)
2246	endif
2247	qlpqi=qlz(k)+qiz(k)
2248	if ( qlpqi .eq. 0.0 ) then
2249	qvsbar(k)=qsiz(k)
2250	else
2251	qvsbar(k)=( qiz(k)qsiz(k)+qlz(k)qswz(k) )/qlpqi
2252	endif
2253
2254	1800 continue
2255	!
2256	!***********************************************************************
2257	!******** integrate the productions of rain and snow ********
2258	!***********************************************************************
2259	!c
2260
2261	2000 continue
2262
2263
2264	!---------------------------------------------------------------------
2265
2266	!
2267	!***********************************************************************
2268	!**** Write terms in cloud physics to time series dataset ***
2269	!***********************************************************************
2270	!
2271	! open(unit=24,form='formatted',status='new',
2272	! & file='cloud.dat')
2273
2274	!9030 format(10e12.6)
2275
2276	! write(24,*)'tmp'
2277	! write(24,9030) (tem(k),k=kts+1,kte)
2278	! write(24,*)'qiz'
2279	! write(24,9030) (qiz(k),k=kts+1,kte)
2280	! write(24,*)'qsz'
2281	! write(24,9030) (qsz(k),k=kts+1,kte)
2282	! write(24,*)'qrz'
2283	! write(24,9030) (qrz(k),k=kts+1,kte)
2284	! write(24,*)'qgz'
2285	! write(24,9030) (qgz(k),k=kts+1,kte)
2286	! write(24,*)'qvoqsw'
2287	! write(24,9030) (qvoqswz(k),k=kts+1,kte)
2288	! write(24,*)'qvoqsi'
2289	! write(24,9030) (qvoqsiz(k),k=kts+1,kte)
2290	! write(24,*)'vtr'
2291	! write(24,9030) (vtr(k),k=kts+1,kte)
2292	! write(24,*)'vts'
2293	! write(24,9030) (vts(k),k=kts+1,kte)
2294	! write(24,*)'vtg'
2295	! write(24,9030) (vtg(k),k=kts+1,kte)
2296	! write(24,*)'pclw'
2297	! write(24,9030) (pclw(k),k=kts+1,kte)
2298	! write(24,*)'pvapor'
2299	! write(24,9030) (pvapor(k),k=kts+1,kte)
2300	! write(24,*)'pcli'
2301	! write(24,9030) (pcli(k),k=kts+1,kte)
2302	! write(24,*)'pimlt'
2303	! write(24,9030) (pimlt(k),k=kts+1,kte)
2304	! write(24,*)'pihom'
2305	! write(24,9030) (pihom(k),k=kts+1,kte)
2306	! write(24,*)'pidw'
2307	! write(24,9030) (pidw(k),k=kts+1,kte)
2308	! write(24,*)'prain'
2309	! write(24,9030) (prain(k),k=kts+1,kte)
2310	! write(24,*)'praut'
2311	! write(24,9030) (praut(k),k=kts+1,kte)
2312	! write(24,*)'pracw'
2313	! write(24,9030) (pracw(k),k=kts+1,kte)
2314	! write(24,*)'prevp'
2315	! write(24,9030) (prevp(k),k=kts+1,kte)
2316	! write(24,*)'psnow'
2317	! write(24,9030) (psnow(k),k=kts+1,kte)
2318	! write(24,*)'psaut'
2319	! write(24,9030) (psaut(k),k=kts+1,kte)
2320	! write(24,*)'psfw'
2321	! write(24,9030) (psfw(k),k=kts+1,kte)
2322	! write(24,*)'psfi'
2323	! write(24,9030) (psfi(k),k=kts+1,kte)
2324	! write(24,*)'praci'
2325	! write(24,9030) (praci(k),k=kts+1,kte)
2326	! write(24,*)'piacr'
2327	! write(24,9030) (piacr(k),k=kts+1,kte)
2328	! write(24,*)'psaci'
2329	! write(24,9030) (psaci(k),k=kts+1,kte)
2330	! write(24,*)'psacw'
2331	! write(24,9030) (psacw(k),k=kts+1,kte)
2332	! write(24,*)'psdep'
2333	! write(24,9030) (psdep(k),k=kts+1,kte)
2334	! write(24,*)'pssub'
2335	! write(24,9030) (pssub(k),k=kts+1,kte)
2336	! write(24,*)'pracs'
2337	! write(24,9030) (pracs(k),k=kts+1,kte)
2338	! write(24,*)'psacr'
2339	! write(24,9030) (psacr(k),k=kts+1,kte)
2340	! write(24,*)'psmlt'
2341	! write(24,9030) (psmlt(k),k=kts+1,kte)
2342	! write(24,*)'psmltevp'
2343	! write(24,9030) (psmltevp(k),k=kts+1,kte)
2344	! write(24,*)'pladj'
2345	! write(24,9030) (pladj(k),k=kts+1,kte)
2346	! write(24,*)'piadj'
2347	! write(24,9030) (piadj(k),k=kts+1,kte)
2348	! write(24,*)'pgraupel'
2349	! write(24,9030) (pgraupel(k),k=kts+1,kte)
2350	! write(24,*)'pgaut'
2351	! write(24,9030) (pgaut(k),k=kts+1,kte)
2352	! write(24,*)'pgfr'
2353	! write(24,9030) (pgfr(k),k=kts+1,kte)
2354	! write(24,*)'pgacw'
2355	! write(24,9030) (pgacw(k),k=kts+1,kte)
2356	! write(24,*)'pgaci'
2357	! write(24,9030) (pgaci(k),k=kts+1,kte)
2358	! write(24,*)'pgacr'
2359	! write(24,9030) (pgacr(k),k=kts+1,kte)
2360	! write(24,*)'pgacs'
2361	! write(24,9030) (pgacs(k),k=kts+1,kte)
2362	! write(24,*)'pgacip'
2363	! write(24,9030) (pgacip(k),k=kts+1,kte)
2364	! write(24,*)'pgacrP'
2365	! write(24,9030) (pgacrP(k),k=kts+1,kte)
2366	! write(24,*)'pgacsp'
2367	! write(24,9030) (pgacsp(k),k=kts+1,kte)
2368	! write(24,*)'pgwet'
2369	! write(24,9030) (pgwet(k),k=kts+1,kte)
2370	! write(24,*)'pdry'
2371	! write(24,9030) (pdry(k),k=kts+1,kte)
2372	! write(24,*)'pgsub'
2373	! write(24,9030) (pgsub(k),k=kts+1,kte)
2374	! write(24,*)'pgdep'
2375	! write(24,9030) (pgdep(k),k=kts+1,kte)
2376	! write(24,*)'pgmlt'
2377	! write(24,9030) (pgmlt(k),k=kts+1,kte)
2378	! write(24,*)'pgmltevp'
2379	! write(24,9030) (pgmltevp(k),k=kts+1,kte)
2380
2381
2382
2383	!**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
2384	!
2385	do k=kts+1,kte
2386	if ( qvz(k) .lt. qvmin ) then
2387	qlz(k)=0.0
2388	qiz(k)=0.0
2389	qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
2390	end if
2391	enddo
2392	!
2393	END SUBROUTINE clphy1d
2394
2395
2396	!---------------------------------------------------------------------
2397	! SATURATED ADJUSTMENT
2398	!---------------------------------------------------------------------
2399	SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, &
2400	kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
2401	!---------------------------------------------------------------------
2402	IMPLICIT NONE
2403	!---------------------------------------------------------------------
2404	! This program use Newton's method for finding saturated temperature
2405	! and saturation mixing ratio.
2406	!
2407	! In this saturation adjustment scheme we assume
2408	! (1) the saturation mixing ratio is the mass weighted average of
2409	! saturation values over liquid water (qsw), and ice (qsi)
2410	! following Lord et al., 1984 and Tao, 1989
2411	!
2412	! (2) the percentage of cloud liquid and cloud ice will
2413	! be fixed during the saturation calculation
2414	!---------------------------------------------------------------------
2415	!
2416
2417	INTEGER, INTENT(IN ) :: kts, kte, k
2418
2419	REAL, DIMENSION( kts:kte ), &
2420	INTENT(INOUT) :: qvz, qlz, qiz
2421	!
2422	REAL, DIMENSION( kts:kte ), &
2423	INTENT(IN ) :: prez, theiz, tothz
2424
2425	REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
2426	REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
2427
2428	! LOCAL VARS
2429
2430	INTEGER :: n
2431
2432	REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
2433	qswz, qvsbar
2434
2435	REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
2436	denom1, denom2, dqvsbar, ftsat, dftsat, qpz, &
2437	gindex, es
2438	!
2439	!---------------------------------------------------------------------
2440
2441	thz(k)=theiz(k)-(xLvocpqvz(k)-xLfocpqiz(k))/tothz(k)
2442
2443	tem(k)=tothz(k)*thz(k)
2444	if (tem(k) .gt. 273.15) then
2445	! qsat=episp0k/prez(k)* &
2446	! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
2447	es=1000.svp1exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
2448	qsat=ep2*es/(prez(k)-es)
2449	else
2450	qsat=episp0k/prez(k)* &
2451	exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
2452	end if
2453	qpz=qvz(k)+qlz(k)+qiz(k)
2454	if (qpz .lt. qsat) then
2455	qvz(k)=qpz
2456	qiz(k)=0.0
2457	qlz(k)=0.0
2458	go to 400
2459	end if
2460	qlpqi=qlz(k)+qiz(k)
2461	if( qlpqi .ge. 1.0e-5) then
2462	ratql=qlz(k)/qlpqi
2463	ratqi=qiz(k)/qlpqi
2464	else
2465	t0=273.15
2466	! t1=233.15
2467	t1=248.15
2468	tmp1=( t0-tem(k) )/(t0-t1)
2469	tmp1=amin1(1.0,tmp1)
2470	tmp1=amax1(0.0,tmp1)
2471	ratqi=tmp1
2472	ratql=1.0-tmp1
2473	end if
2474	!
2475	!
2476	!-- saturation mixing ratios over water and ice
2477	!-- at the outset we will follow Bolton 1980 MWR for
2478	!-- the water and Murray JAS 1967 for the ice
2479	!
2480	!-- dqvsbar=d(qvsbar)/dT
2481	!-- ftsat=F(Tsat)
2482	!-- dftsat=d(F(T))/dT
2483	!
2484	! First guess of tsat
2485
2486	tsat=tem(k)
2487	absft=1.0
2488	!
2489	do 200 n=1,20
2490	denom1=1.0/(tsat-svp3)
2491	denom2=1.0/(tsat-7.66)
2492	! qswz(k)=episp0k/prez(k)* &
2493	! exp( svp2denom1(tsat-273.15) )
2494	es=1000.svp1exp( svp2denom1(tsat-svpt0) )
2495	qswz(k)=ep2*es/(prez(k)-es)
2496	if (tem(k) .lt. 273.15) then
2497	! qsiz(k)=episp0k/prez(k)* &
2498	! exp( 21.8745584denom2(tsat-273.15) )
2499	es=1000.svp1exp( 21.8745584denom2(tsat-273.15) )
2500	qsiz(k)=ep2*es/(prez(k)-es)
2501	if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
2502	else
2503	qsiz(k)=qswz(k)
2504	endif
2505	qvsbar(k)=ratqlqswz(k)+ratqiqsiz(k)
2506	!
2507	! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
2508	if( absft .lt. 0.01 ) go to 300
2509	!
2510	dqvsbar=ratqlqswz(k)svp2243.5denom1*denom1+ &
2511	ratqiqsiz(k)21.8745584265.5denom2*denom2
2512	ftsat=tsat+(xlvocp+ratqixlfocp)qvsbar(k)- &
2513	tothz(k)theiz(k)-xlfocpratqi*(qvz(k)+qlz(k)+qiz(k))
2514	dftsat=1.0+(xlvocp+ratqixlfocp)dqvsbar
2515	tsat=tsat-ftsat/dftsat
2516	absft=abs(ftsat)
2517
2518	200 continue
2519	9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
2520	!
2521	300 continue
2522	if( qpz .gt. qvsbar(k) ) then
2523	qvz(k)=qvsbar(k)
2524	qiz(k)=ratqi*( qpz-qvz(k) )
2525	qlz(k)=ratql*( qpz-qvz(k) )
2526	else
2527	qvz(k)=qpz
2528	qiz(k)=0.0
2529	qlz(k)=0.0
2530	end if
2531	400 continue
2532
2533	END SUBROUTINE satadj
2534
2535
2536	!----------------------------------------------------------------
2537	REAL FUNCTION parama1(temp)
2538	!----------------------------------------------------------------
2539	IMPLICIT NONE
2540	!----------------------------------------------------------------
2541	! This program calculate the parameter for crystal growth rate
2542	! in Bergeron process
2543	!----------------------------------------------------------------
2544
2545	REAL, INTENT (IN ) :: temp
2546	REAL, DIMENSION(32) :: a1
2547	INTEGER :: i1, i1p1
2548	REAL :: ratio
2549
2550	data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
2551	0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
2552	0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
2553	0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
2554	0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
2555	0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
2556	0.5333e-6,0.4834e-6/
2557
2558	i1=int(-temp)+1
2559	i1p1=i1+1
2560	ratio=-(temp)-float(i1-1)
2561	parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
2562
2563	END FUNCTION parama1
2564
2565	!----------------------------------------------------------------
2566	REAL FUNCTION parama2(temp)
2567	!----------------------------------------------------------------
2568	IMPLICIT NONE
2569	!----------------------------------------------------------------
2570	! This program calculate the parameter for crystal growth rate
2571	! in Bergeron process
2572	!----------------------------------------------------------------
2573
2574	REAL, INTENT (IN ) :: temp
2575	REAL, DIMENSION(32) :: a2
2576	INTEGER :: i1, i1p1
2577	REAL :: ratio
2578
2579	data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
2580	0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
2581	0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
2582	0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
2583	0.4433,0.4413,0.4382,0.4361/
2584	i1=int(-temp)+1
2585	i1p1=i1+1
2586	ratio=-(temp)-float(i1-1)
2587	parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
2588
2589	END FUNCTION parama2
2590
2591	!----------------------------------------------------------------
2592	REAL FUNCTION ggamma(X)
2593	!----------------------------------------------------------------
2594	IMPLICIT NONE
2595	!----------------------------------------------------------------
2596	REAL, INTENT(IN ) :: x
2597	REAL, DIMENSION(8) :: B
2598	INTEGER ::j, K1
2599	REAL ::PF, G1TO2 ,TEMP
2600
2601	DATA B/-.577191652,.988205891,-.897056937,.918206857, &
2602	-.756704078,.482199394,-.193527818,.035868343/
2603
2604	PF=1.
2605	TEMP=X
2606	DO 10 J=1,200
2607	IF (TEMP .LE. 2) GO TO 20
2608	TEMP=TEMP-1.
2609	10 PF=PF*TEMP
2610	100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
2611	WRITE(wrf_err_message,100)X
2612	CALL wrf_error_fatal(wrf_err_message)
2613	20 G1TO2=1.
2614	TEMP=TEMP - 1.
2615	DO 30 K1=1,8
2616	30 G1TO2=G1TO2 + B(K1)TEMP*K1
2617	ggamma=PF*G1TO2
2618
2619	END FUNCTION ggamma
2620
2621	!----------------------------------------------------------------
2622
2623	END MODULE module_mp_lin
2624

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

Download in other formats:

Original Format