Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_lin.F @ 6

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

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

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

Download in other formats:

Original Format