Context Navigation

module_mp_lin.F

Last change on this file was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 83.8 KB

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

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

Download in other formats:

Original Format