Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_lin.F @ 3667

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

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

Download in other formats:

Original Format