Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_wsm6.F @ 3806

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

Rev	Line
[2759]	1	#if ( RWORDSIZE == 4 )
	2	# define VREC vsrec
	3	# define VSQRT vssqrt
	4	#else
	5	# define VREC vrec
	6	# define VSQRT vsqrt
	7	#endif
	8
	9	MODULE module_mp_wsm6
	10	!
	11	!
	12	REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
	13	REAL, PARAMETER, PRIVATE :: n0r = 8.e6
	14	REAL, PARAMETER, PRIVATE :: n0g = 4.e6
	15	REAL, PARAMETER, PRIVATE :: avtr = 841.9
	16	REAL, PARAMETER, PRIVATE :: bvtr = 0.8
	17	REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
	18	REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
	19	REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
	20	REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
	21	REAL, PARAMETER, PRIVATE :: avts = 11.72
	22	REAL, PARAMETER, PRIVATE :: bvts = .41
	23	REAL, PARAMETER, PRIVATE :: avtg = 330.
	24	REAL, PARAMETER, PRIVATE :: bvtg = 0.8
	25	REAL, PARAMETER, PRIVATE :: deng = 500.
	26	REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
	27	REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
	28	REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
	29	REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
	30	REAL, PARAMETER, PRIVATE :: betai = .6
	31	REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
	32	REAL, PARAMETER, PRIVATE :: dicon = 11.9
	33	REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
	34	REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
	35	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
	36	REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
	37	REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
	38	REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
	39	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
	40	REAL, PARAMETER, PRIVATE :: t40c = 233.16
	41	REAL, PARAMETER, PRIVATE :: eacrc = 1.0
	42	REAL, PARAMETER, PRIVATE :: dens = 100.0
	43	REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! pgaut
	44	REAL, SAVE :: &
	45	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
	46	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
	47	bvtr6,g6pbr, &
	48	precr1,precr2,xm0,xmmax,roqimax,bvts1, &
	49	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
	50	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
	51	pidn0s,xlv1,pacrc, &
	52	bvtg1,bvtg2,bvtg3,bvtg4,g1pbg, &
	53	g3pbg,g4pbg,g5pbgo2,pvtg,pacrg, &
	54	precg1,precg2,pidn0g, &
	55	rslopermax,rslopesmax,rslopegmax, &
	56	rsloperbmax,rslopesbmax,rslopegbmax, &
	57	rsloper2max,rslopes2max,rslopeg2max, &
	58	rsloper3max,rslopes3max,rslopeg3max
	59	CONTAINS
	60	!===================================================================
	61	!
	62	SUBROUTINE wsm6(th, q, qc, qr, qi, qs, qg &
	63	,den, pii, p, delz &
	64	,delt,g, cpd, cpv, rd, rv, t0c &
	65	,ep1, ep2, qmin &
	66	,XLS, XLV0, XLF0, den0, denr &
	67	,cliq,cice,psat &
	68	,rain, rainncv &
	69	,snow, snowncv &
	70	,graupel, graupelncv &
	71	,sr &
	72	,ids,ide, jds,jde, kds,kde &
	73	,ims,ime, jms,jme, kms,kme &
	74	,its,ite, jts,jte, kts,kte &
	75	)
	76	!-------------------------------------------------------------------
	77	IMPLICIT NONE
	78	!-------------------------------------------------------------------
	79	!
	80	! This code is a 6-class GRAUPEL phase microphyiscs scheme (WSM6) of the WRF
	81	! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
	82	! number concentration is a function of temperature, and seperate assumption
	83	! is developed, in which ice crystal number concentration is a function
	84	! of ice amount. A theoretical background of the ice-microphysics and related
	85	! processes in the WSMMPs are described in Hong et al. (2004).
	86	! All production terms in the WSM6 scheme are described in Hong and Lim (2006).
	87	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
	88	!
	89	! WSM6 cloud scheme
	90	!
	91	! Coded by Song-You Hong and Jeong-Ock Jade Lim (Yonsei Univ.)
	92	! Summer 2003
	93	!
	94	! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
	95	! Summer 2004
	96	!
	97	! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
	98	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
	99	! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
	100	! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
	101	! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
	102	! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
	103	!
	104	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
	105	ims,ime, jms,jme, kms,kme , &
	106	its,ite, jts,jte, kts,kte
	107	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	108	INTENT(INOUT) :: &
	109	th, &
	110	q, &
	111	qc, &
	112	qi, &
	113	qr, &
	114	qs, &
	115	qg
	116	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	117	INTENT(IN ) :: &
	118	den, &
	119	pii, &
	120	p, &
	121	delz
	122	REAL, INTENT(IN ) :: delt, &
	123	g, &
	124	rd, &
	125	rv, &
	126	t0c, &
	127	den0, &
	128	cpd, &
	129	cpv, &
	130	ep1, &
	131	ep2, &
	132	qmin, &
	133	XLS, &
	134	XLV0, &
	135	XLF0, &
	136	cliq, &
	137	cice, &
	138	psat, &
	139	denr
	140	REAL, DIMENSION( ims:ime , jms:jme ), &
	141	INTENT(INOUT) :: rain, &
	142	rainncv, &
	143	sr
	144
	145	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
	146	INTENT(INOUT) :: snow, &
	147	snowncv
	148
	149	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
	150	INTENT(INOUT) :: graupel, &
	151	graupelncv
	152	! LOCAL VAR
	153	REAL, DIMENSION( its:ite , kts:kte ) :: t
	154	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
	155	REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs
	156	INTEGER :: i,j,k
	157	!-------------------------------------------------------------------
	158	DO j=jts,jte
	159	DO k=kts,kte
	160	DO i=its,ite
	161	t(i,k)=th(i,k,j)*pii(i,k,j)
	162	qci(i,k,1) = qc(i,k,j)
	163	qci(i,k,2) = qi(i,k,j)
	164	qrs(i,k,1) = qr(i,k,j)
	165	qrs(i,k,2) = qs(i,k,j)
	166	qrs(i,k,3) = qg(i,k,j)
	167	ENDDO
	168	ENDDO
	169
	170	! Sending array starting locations of optional variables may cause
	171	! troubles, so we explicitly change the call.
	172
	173	CALL wsm62D(t, q(ims,kms,j), qci, qrs &
	174	,den(ims,kms,j) &
	175	,p(ims,kms,j), delz(ims,kms,j) &
	176	,delt,g, cpd, cpv, rd, rv, t0c &
	177	,ep1, ep2, qmin &
	178	,XLS, XLV0, XLF0, den0, denr &
	179	,cliq,cice,psat &
	180	,j &
	181	,rain(ims,j),rainncv(ims,j) &
	182	,sr(ims,j) &
	183	,ids,ide, jds,jde, kds,kde &
	184	,ims,ime, jms,jme, kms,kme &
	185	,its,ite, jts,jte, kts,kte &
	186	#if ( EM_CORE == 1 )
	187	,snow(ims,j),snowncv(ims,j) &
	188	,graupel(ims,j),graupelncv(ims,j) &
	189	#endif
	190	)
	191
	192	DO K=kts,kte
	193	DO I=its,ite
	194	th(i,k,j)=t(i,k)/pii(i,k,j)
	195	qc(i,k,j) = qci(i,k,1)
	196	qi(i,k,j) = qci(i,k,2)
	197	qr(i,k,j) = qrs(i,k,1)
	198	qs(i,k,j) = qrs(i,k,2)
	199	qg(i,k,j) = qrs(i,k,3)
	200	ENDDO
	201	ENDDO
	202	ENDDO
	203	END SUBROUTINE wsm6
	204	!===================================================================
	205	!
	206	SUBROUTINE wsm62D(t, q, qci, qrs, den, p, delz &
	207	,delt,g, cpd, cpv, rd, rv, t0c &
	208	,ep1, ep2, qmin &
	209	,XLS, XLV0, XLF0, den0, denr &
	210	,cliq,cice,psat &
	211	,lat &
	212	,rain,rainncv &
	213	,sr &
	214	,ids,ide, jds,jde, kds,kde &
	215	,ims,ime, jms,jme, kms,kme &
	216	,its,ite, jts,jte, kts,kte &
	217	,snow,snowncv &
	218	,graupel,graupelncv &
	219	)
	220	!-------------------------------------------------------------------
	221	IMPLICIT NONE
	222	!-------------------------------------------------------------------
	223	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
	224	ims,ime, jms,jme, kms,kme , &
	225	its,ite, jts,jte, kts,kte, &
	226	lat
	227	REAL, DIMENSION( its:ite , kts:kte ), &
	228	INTENT(INOUT) :: &
	229	t
	230	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
	231	INTENT(INOUT) :: &
	232	qci
	233	REAL, DIMENSION( its:ite , kts:kte, 3 ), &
	234	INTENT(INOUT) :: &
	235	qrs
	236	REAL, DIMENSION( ims:ime , kms:kme ), &
	237	INTENT(INOUT) :: &
	238	q
	239	REAL, DIMENSION( ims:ime , kms:kme ), &
	240	INTENT(IN ) :: &
	241	den, &
	242	p, &
	243	delz
	244	REAL, INTENT(IN ) :: delt, &
	245	g, &
	246	cpd, &
	247	cpv, &
	248	t0c, &
	249	den0, &
	250	rd, &
	251	rv, &
	252	ep1, &
	253	ep2, &
	254	qmin, &
	255	XLS, &
	256	XLV0, &
	257	XLF0, &
	258	cliq, &
	259	cice, &
	260	psat, &
	261	denr
	262	REAL, DIMENSION( ims:ime ), &
	263	INTENT(INOUT) :: rain, &
	264	rainncv, &
	265	sr
	266	REAL, DIMENSION( ims:ime ), OPTIONAL, &
	267	INTENT(INOUT) :: snow, &
	268	snowncv
	269
	270	REAL, DIMENSION( ims:ime ), OPTIONAL, &
	271	INTENT(INOUT) :: graupel, &
	272	graupelncv
	273	! LOCAL VAR
	274	REAL, DIMENSION( its:ite , kts:kte , 3) :: &
	275	rh, qs, rslope, rslope2, rslope3, rslopeb, &
	276	falk, fall, work1
	277	REAL, DIMENSION( its:ite , kts:kte ) :: &
	278	worka
	279	REAL, DIMENSION( its:ite , kts:kte ) :: &
	280	falkc, work1c, work2c, fallc
	281	REAL, DIMENSION( its:ite , kts:kte ) :: &
	282	prevp, psdep, pgdep, praut, psaut, pgaut, &
	283	pracw, psacw, pgacw, pgacr, pgacs, psaci, pgmlt, praci, &
	284	piacr, pracs, psacr, pgaci, pseml, pgeml
	285	REAL, DIMENSION( its:ite , kts:kte) :: qsum
	286	REAL, DIMENSION( its:ite , kts:kte ) :: paacw
	287	REAL, DIMENSION( its:ite , kts:kte ) :: &
	288	pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, denfac, &
	289	xni, pgevp,n0sfac
	290	! variables for optimization
	291	REAL, DIMENSION( its:ite ) :: tvec1
	292	REAL :: temp
	293	INTEGER, DIMENSION( its:ite ) :: mstep, numdt
	294	LOGICAL, DIMENSION( its:ite ) :: flgcld
	295	REAL :: pi, &
	296	cpmcal, xlcal, lamdar, lamdas, lamdag, diffus, &
	297	viscos, xka, venfac, conden, diffac, &
	298	x, y, z, a, b, c, d, e, &
	299	qdt, holdrr, holdrs, holdrg, supcol, supcolt, pvt, &
	300	coeres, supsat, dtcld, xmi, eacrs, satdt, &
	301	qimax, diameter, xni0, roqi0, &
	302	fallsum, fallsum_qsi, fallsum_qg, &
	303	vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
	304	xlwork2, factor, source, value, &
	305	xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
	306	REAL :: vt2ave
	307	REAL :: holdc, holdci
	308	INTEGER :: i, j, k, mstepmax, &
	309	iprt, latd, lond, loop, loops, ifsat, n
	310	! Temporaries used for inlining fpvs function
	311	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
	312	!
	313	!=================================================================
	314	! compute internal functions
	315	!
	316	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
	317	xlcal(x) = xlv0-xlv1*(x-t0c)
	318	!----------------------------------------------------------------
	319	! size distributions: (x=mixing ratio, y=air density):
	320	! valid for mixing ratio > 1.e-9 kg/kg.
	321	!
	322	! Optimizatin : A*B => exp(log(A)(B))
	323	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
	324	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
	325	lamdag(x,y)= sqrt(sqrt(pidn0g/(xy))) ! (pidn0g/(xy))**.25
	326	!
	327	!----------------------------------------------------------------
	328	! diffus: diffusion coefficient of the water vapor
	329	! viscos: kinematic viscosity(m2s-1)
	330	!
	331	diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
	332	viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
	333	xka(x,y) = 1.414e3viscos(x,y)y
	334	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
	335	venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
	336	/sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
	337	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
	338	!
	339	pi = 4. * atan(1.)
	340	!
	341	!
	342	!----------------------------------------------------------------
	343	! paddint 0 for negative values generated by dynamics
	344	!
	345	do k = kts, kte
	346	do i = its, ite
	347	qci(i,k,1) = max(qci(i,k,1),0.0)
	348	qrs(i,k,1) = max(qrs(i,k,1),0.0)
	349	qci(i,k,2) = max(qci(i,k,2),0.0)
	350	qrs(i,k,2) = max(qrs(i,k,2),0.0)
	351	qrs(i,k,3) = max(qrs(i,k,3),0.0)
	352	enddo
	353	enddo
	354	!
	355	!----------------------------------------------------------------
	356	! latent heat for phase changes and heat capacity. neglect the
	357	! changes during microphysical process calculation
	358	! emanuel(1994)
	359	!
	360	do k = kts, kte
	361	do i = its, ite
	362	cpm(i,k) = cpmcal(q(i,k))
	363	xl(i,k) = xlcal(t(i,k))
	364	enddo
	365	enddo
	366	!
	367	!----------------------------------------------------------------
	368	! compute the minor time steps.
	369	!
	370	loops = max(nint(delt/dtcldcr),1)
	371	dtcld = delt/loops
	372	if(delt.le.dtcldcr) dtcld = delt
	373	!
	374	do loop = 1,loops
	375	!
	376	!----------------------------------------------------------------
	377	! initialize the large scale variables
	378	!
	379	do i = its, ite
	380	mstep(i) = 1
	381	flgcld(i) = .true.
	382	enddo
	383	!
	384	! do k = kts, kte
	385	! do i = its, ite
	386	! denfac(i,k) = sqrt(den0/den(i,k))
	387	! enddo
	388	! enddo
	389	do k = kts, kte
	390	CALL VREC( tvec1(its), den(its,k), ite-its+1)
	391	do i = its, ite
	392	tvec1(i) = tvec1(i)*den0
	393	enddo
	394	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
	395	enddo
	396	!
	397	! Inline expansion for fpvs
	398	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	399	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	400	hsub = xls
	401	hvap = xlv0
	402	cvap = cpv
	403	ttp=t0c+0.01
	404	dldt=cvap-cliq
	405	xa=-dldt/rv
	406	xb=xa+hvap/(rv*ttp)
	407	dldti=cvap-cice
	408	xai=-dldti/rv
	409	xbi=xai+hsub/(rv*ttp)
	410	do k = kts, kte
	411	do i = its, ite
	412	tr=ttp/t(i,k)
	413	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	414	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
	415	qs(i,k,1) = max(qs(i,k,1),qmin)
	416	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
	417	tr=ttp/t(i,k)
	418	if(t(i,k).lt.ttp) then
	419	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
	420	else
	421	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	422	endif
	423	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
	424	qs(i,k,2) = max(qs(i,k,2),qmin)
	425	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
	426	enddo
	427	enddo
	428	!
	429	!----------------------------------------------------------------
	430	! initialize the variables for microphysical physics
	431	!
	432	!
	433	do k = kts, kte
	434	do i = its, ite
	435	prevp(i,k) = 0.
	436	psdep(i,k) = 0.
	437	pgdep(i,k) = 0.
	438	praut(i,k) = 0.
	439	psaut(i,k) = 0.
	440	pgaut(i,k) = 0.
	441	pracw(i,k) = 0.
	442	praci(i,k) = 0.
	443	piacr(i,k) = 0.
	444	psaci(i,k) = 0.
	445	psacw(i,k) = 0.
	446	pracs(i,k) = 0.
	447	psacr(i,k) = 0.
	448	pgacw(i,k) = 0.
	449	paacw(i,k) = 0.
	450	pgaci(i,k) = 0.
	451	pgacr(i,k) = 0.
	452	pgacs(i,k) = 0.
	453	pigen(i,k) = 0.
	454	pidep(i,k) = 0.
	455	pcond(i,k) = 0.
	456	psmlt(i,k) = 0.
	457	pgmlt(i,k) = 0.
	458	pseml(i,k) = 0.
	459	pgeml(i,k) = 0.
	460	psevp(i,k) = 0.
	461	pgevp(i,k) = 0.
	462	falk(i,k,1) = 0.
	463	falk(i,k,2) = 0.
	464	falk(i,k,3) = 0.
	465	fall(i,k,1) = 0.
	466	fall(i,k,2) = 0.
	467	fall(i,k,3) = 0.
	468	fallc(i,k) = 0.
	469	falkc(i,k) = 0.
	470	xni(i,k) = 1.e3
	471	enddo
	472	enddo
	473	!
	474	!----------------------------------------------------------------
	475	! compute the fallout term:
	476	! first, vertical terminal velosity for minor loops
	477	!
	478	do k = kts, kte
	479	do i = its, ite
	480	supcol = t0c-t(i,k)
	481	!---------------------------------------------------------------
	482	! n0s: Intercept parameter for snow [m-4] [HDC 6]
	483	!---------------------------------------------------------------
	484	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
	485	if(qrs(i,k,1).le.qcrmin)then
	486	rslope(i,k,1) = rslopermax
	487	rslopeb(i,k,1) = rsloperbmax
	488	rslope2(i,k,1) = rsloper2max
	489	rslope3(i,k,1) = rsloper3max
	490	else
	491	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
	492	rslopeb(i,k,1) = rslope(i,k,1)**bvtr
	493	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
	494	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
	495	endif
	496	if(qrs(i,k,2).le.qcrmin)then
	497	rslope(i,k,2) = rslopesmax
	498	rslopeb(i,k,2) = rslopesbmax
	499	rslope2(i,k,2) = rslopes2max
	500	rslope3(i,k,2) = rslopes3max
	501	else
	502	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
	503	rslopeb(i,k,2) = rslope(i,k,2)**bvts
	504	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
	505	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
	506	endif
	507	if(qrs(i,k,3).le.qcrmin)then
	508	rslope(i,k,3) = rslopegmax
	509	rslopeb(i,k,3) = rslopegbmax
	510	rslope2(i,k,3) = rslopeg2max
	511	rslope3(i,k,3) = rslopeg3max
	512	else
	513	rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
	514	rslopeb(i,k,3) = rslope(i,k,3)**bvtg
	515	rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
	516	rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
	517	endif
	518	!-------------------------------------------------------------
	519	! Ni: ice crystal number concentraiton [HDC 5c]
	520	!-------------------------------------------------------------
	521	! xni(i,k) = min(max(5.38e7*(den(i,k) &
	522	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
	523	temp = (den(i,k)*max(qci(i,k,2),qmin))
	524	temp = sqrt(sqrt(temptemptemp))
	525	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
	526	enddo
	527	enddo
	528	!
	529	mstepmax = 1
	530	numdt = 1
	531	do k = kte, kts, -1
	532	do i = its, ite
	533	work1(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)/delz(i,k)
	534	work1(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)/delz(i,k)
	535	work1(i,k,3) = pvtgrslopeb(i,k,3)denfac(i,k)/delz(i,k)
	536	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
	537	IF ( qsum(i,k) .gt. 1.e-15 ) THEN
	538	worka(i,k) = (work1(i,k,2)qrs(i,k,2) + work1(i,k,3)qrs(i,k,3))/qsum(i,k)
	539	ELSE
	540	worka(i,k) = 0.
	541	ENDIF
	542	numdt(i) = max(nint(max(work1(i,k,1),worka(i,k)) &
	543	*dtcld+.5),1)
	544	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
	545	enddo
	546	enddo
	547	do i = its, ite
	548	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
	549	enddo
	550	!
	551	do n = 1, mstepmax
	552	k = kte
	553	do i = its, ite
	554	if(n.le.mstep(i)) then
	555	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
	556	falk(i,k,2) = den(i,k)qrs(i,k,2)worka(i,k)/mstep(i)
	557	falk(i,k,3) = den(i,k)qrs(i,k,3)worka(i,k)/mstep(i)
	558	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
	559	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
	560	fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
	561	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
	562	qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
	563	qrs(i,k,3) = max(qrs(i,k,3)-falk(i,k,3)*dtcld/den(i,k),0.)
	564	endif
	565	enddo
	566	do k = kte-1, kts, -1
	567	do i = its, ite
	568	if(n.le.mstep(i)) then
	569	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
	570	falk(i,k,2) = den(i,k)qrs(i,k,2)worka(i,k)/mstep(i)
	571	falk(i,k,3) = den(i,k)qrs(i,k,3)worka(i,k)/mstep(i)
	572	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
	573	fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
	574	fall(i,k,3) = fall(i,k,3)+falk(i,k,3)
	575	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
	576	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	577	qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
	578	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	579	qrs(i,k,3) = max(qrs(i,k,3)-(falk(i,k,3)-falk(i,k+1,3) &
	580	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	581	endif
	582	enddo
	583	enddo
	584	do k = kte, kts, -1
	585	do i = its, ite
	586	if(n.le.mstep(i).and.t(i,k).gt.t0c) then
	587	!---------------------------------------------------------------
	588	! psmlt: melting of snow [HL A33] [RH83 A25]
	589	! (T>T0: S->R)
	590	!---------------------------------------------------------------
	591	xlf = xlf0
	592	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
	593	if(qrs(i,k,2).gt.0.) then
	594	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
	595	psmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
	596	n0sfac(i,k)(precs1*rslope2(i,k,2) &
	597	+precs2work2(i,k)coeres)
	598	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
	599	-qrs(i,k,2)/mstep(i)),0.)
	600	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
	601	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
	602	t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
	603	endif
	604	!---------------------------------------------------------------
	605	! pgmlt: melting of graupel [HL A23] [LFO 47]
	606	! (T>T0: G->R)
	607	!---------------------------------------------------------------
	608	if(qrs(i,k,3).gt.0.) then
	609	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
	610	pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf &
	611	(t0c-t(i,k))(precg1*rslope2(i,k,3) &
	612	+precg2work2(i,k)coeres)
	613	pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
	614	-qrs(i,k,3)/mstep(i)),0.)
	615	qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
	616	qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
	617	t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
	618	endif
	619	endif
	620	enddo
	621	enddo
	622	enddo
	623	!---------------------------------------------------------------
	624	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
	625	!---------------------------------------------------------------
	626	mstepmax = 1
	627	mstep = 1
	628	numdt = 1
	629	do k = kte, kts, -1
	630	do i = its, ite
	631	if(qci(i,k,2).le.0.) then
	632	work2c(i,k) = 0.
	633	else
	634	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
	635	! diameter = min(dicon * sqrt(xmi),dimax)
	636	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
	637	work1c(i,k) = 1.49e4diameter*1.31
	638	work2c(i,k) = work1c(i,k)/delz(i,k)
	639	endif
	640	numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
	641	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
	642	enddo
	643	enddo
	644	do i = its, ite
	645	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
	646	enddo
	647	!
	648	do n = 1, mstepmax
	649	k = kte
	650	do i = its, ite
	651	if(n.le.mstep(i)) then
	652	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
	653	holdc = falkc(i,k)
	654	fallc(i,k) = fallc(i,k)+falkc(i,k)
	655	holdci = qci(i,k,2)
	656	qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
	657	endif
	658	enddo
	659	do k = kte-1, kts, -1
	660	do i = its, ite
	661	if(n.le.mstep(i)) then
	662	falkc(i,k) = den(i,k)qci(i,k,2)work2c(i,k)/mstep(i)
	663	holdc = falkc(i,k)
	664	fallc(i,k) = fallc(i,k)+falkc(i,k)
	665	holdci = qci(i,k,2)
	666	qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
	667	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	668	endif
	669	enddo
	670	enddo
	671	enddo
	672	!
	673	!----------------------------------------------------------------
	674	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
	675	!
	676	do i = its, ite
	677	fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
	678	fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
	679	fallsum_qg = fall(i,kts,3)
	680	rainncv(i) = 0.
	681	if(fallsum.gt.0.) then
	682	rainncv(i) = fallsumdelz(i,kts)/denrdtcld*1000.
	683	rain(i) = fallsumdelz(i,kts)/denrdtcld*1000. + rain(i)
	684	endif
	685	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
	686	snowncv(i) = 0.
	687	if(fallsum_qsi.gt.0.) then
	688	snowncv(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000.
	689	snow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i)
	690	endif
	691	ENDIF
	692	IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
	693	graupelncv(i) = 0.
	694	if(fallsum_qg.gt.0.) then
	695	graupelncv(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000.
	696	graupel(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000. + graupel(i)
	697	endif
	698	ENDIF
	699	sr(i) = 0.
	700	if(fallsum.gt.0.)sr(i)=(fallsum_qsidelz(i,kts)/denrdtcld*1000. + &
	701	fallsum_qgdelz(i,kts)/denrdtcld*1000.)/(rainncv(i)+1.e-12)
	702	enddo
	703	!
	704	!---------------------------------------------------------------
	705	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
	706	! (T>T0: I->C)
	707	!---------------------------------------------------------------
	708	do k = kts, kte
	709	do i = its, ite
	710	supcol = t0c-t(i,k)
	711	xlf = xls-xl(i,k)
	712	if(supcol.lt.0.) xlf = xlf0
	713	if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
	714	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
	715	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
	716	qci(i,k,2) = 0.
	717	endif
	718	!---------------------------------------------------------------
	719	! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
	720	! (T<-40C: C->I)
	721	!---------------------------------------------------------------
	722	if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
	723	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
	724	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
	725	qci(i,k,1) = 0.
	726	endif
	727	!---------------------------------------------------------------
	728	! pihtf: heterogeneous freezing of cloud water [HL A44]
	729	! (T0>T>-40C: C->I)
	730	!---------------------------------------------------------------
	731	if(supcol.gt.0..and.qci(i,k,1).gt.qmin) then
	732	! pfrzdtc = min(pfrz1(exp(pfrz2supcol)-1.) &
	733	! den(i,k)/denr/xncrqci(i,k,1)*2dtcld,qci(i,k,1))
	734	supcolt=min(supcol,50.)
	735	pfrzdtc = min(pfrz1(exp(pfrz2supcolt)-1.) &
	736	den(i,k)/denr/xncrqci(i,k,1)qci(i,k,1)dtcld,qci(i,k,1))
	737	qci(i,k,2) = qci(i,k,2) + pfrzdtc
	738	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
	739	qci(i,k,1) = qci(i,k,1)-pfrzdtc
	740	endif
	741	!---------------------------------------------------------------
	742	! pgfrz: freezing of rain water [HL A20] [LFO 45]
	743	! (T<T0, R->G)
	744	!---------------------------------------------------------------
	745	if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
	746	! pfrzdtr = min(20.pi2pfrz1n0rdenr/den(i,k) &
	747	! (exp(pfrz2supcol)-1.)rslope3(i,k,1)*2 &
	748	! rslope(i,k,1)dtcld,qrs(i,k,1))
	749	temp = rslope3(i,k,1)
	750	temp = temptemprslope(i,k,1)
	751	supcolt=min(supcol,50.)
	752	pfrzdtr = min(20.(pipi)pfrz1n0r*denr/den(i,k) &
	753	(exp(pfrz2supcolt)-1.)tempdtcld, &
	754	qrs(i,k,1))
	755	qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
	756	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
	757	qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
	758	endif
	759	enddo
	760	enddo
	761	!
	762	!
	763	!----------------------------------------------------------------
	764	! rsloper: reverse of the slope parameter of the rain(m)
	765	! xka: thermal conductivity of air(jm-1s-1k-1)
	766	! work1: the thermodynamic term in the denominator associated with
	767	! heat conduction and vapor diffusion
	768	! (ry88, y93, h85)
	769	! work2: parameter associated with the ventilation effects(y93)
	770	!
	771	do k = kts, kte
	772	do i = its, ite
	773	if(qrs(i,k,1).le.qcrmin)then
	774	rslope(i,k,1) = rslopermax
	775	rslopeb(i,k,1) = rsloperbmax
	776	rslope2(i,k,1) = rsloper2max
	777	rslope3(i,k,1) = rsloper3max
	778	else
	779	rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
	780	rslopeb(i,k,1) = rslope(i,k,1)**bvtr
	781	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
	782	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
	783	endif
	784	if(qrs(i,k,2).le.qcrmin)then
	785	rslope(i,k,2) = rslopesmax
	786	rslopeb(i,k,2) = rslopesbmax
	787	rslope2(i,k,2) = rslopes2max
	788	rslope3(i,k,2) = rslopes3max
	789	else
	790	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
	791	rslopeb(i,k,2) = rslope(i,k,2)**bvts
	792	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
	793	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
	794	endif
	795	if(qrs(i,k,3).le.qcrmin)then
	796	rslope(i,k,3) = rslopegmax
	797	rslopeb(i,k,3) = rslopegbmax
	798	rslope2(i,k,3) = rslopeg2max
	799	rslope3(i,k,3) = rslopeg3max
	800	else
	801	rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
	802	rslopeb(i,k,3) = rslope(i,k,3)**bvtg
	803	rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
	804	rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
	805	endif
	806	enddo
	807	enddo
	808	!
	809	do k = kts, kte
	810	do i = its, ite
	811	work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
	812	work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
	813	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
	814	enddo
	815	enddo
	816	!
	817	!===============================================================
	818	!
	819	! warm rain processes
	820	!
	821	! - follows the processes in RH83 and LFO except for autoconcersion
	822	!
	823	!===============================================================
	824	!
	825	do k = kts, kte
	826	do i = its, ite
	827	supsat = max(q(i,k),qmin)-qs(i,k,1)
	828	satdt = supsat/dtcld
	829	!---------------------------------------------------------------
	830	! praut: auto conversion rate from cloud to rain [HDC 16]
	831	! (C->R)
	832	!---------------------------------------------------------------
	833	if(qci(i,k,1).gt.qc0) then
	834	praut(i,k) = qck1qci(i,k,1)*(7./3.)
	835	praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
	836	endif
	837	!---------------------------------------------------------------
	838	! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
	839	! (C->R)
	840	!---------------------------------------------------------------
	841	if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
	842	pracw(i,k) = min(pacrrrslope3(i,k,1)rslopeb(i,k,1) &
	843	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
	844	endif
	845	!---------------------------------------------------------------
	846	! prevp: evaporation/condensation rate of rain [HDC 14]
	847	! (V->R or R->V)
	848	!---------------------------------------------------------------
	849	if(qrs(i,k,1).gt.0.) then
	850	coeres = rslope2(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
	851	prevp(i,k) = (rh(i,k,1)-1.)(precr1rslope2(i,k,1) &
	852	+precr2work2(i,k)coeres)/work1(i,k,1)
	853	if(prevp(i,k).lt.0.) then
	854	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
	855	prevp(i,k) = max(prevp(i,k),satdt/2)
	856	else
	857	prevp(i,k) = min(prevp(i,k),satdt/2)
	858	endif
	859	endif
	860	enddo
	861	enddo
	862	!
	863	!===============================================================
	864	!
	865	! cold rain processes
	866	!
	867	! - follows the revised ice microphysics processes in HDC
	868	! - the processes same as in RH83 and RH84 and LFO behave
	869	! following ice crystal hapits defined in HDC, inclduing
	870	! intercept parameter for snow (n0s), ice crystal number
	871	! concentration (ni), ice nuclei number concentration
	872	! (n0i), ice diameter (d)
	873	!
	874	!===============================================================
	875	!
	876	do k = kts, kte
	877	do i = its, ite
	878	supcol = t0c-t(i,k)
	879	supsat = max(q(i,k),qmin)-qs(i,k,2)
	880	satdt = supsat/dtcld
	881	ifsat = 0
	882	!-------------------------------------------------------------
	883	! Ni: ice crystal number concentraiton [HDC 5c]
	884	!-------------------------------------------------------------
	885	! xni(i,k) = min(max(5.38e7*(den(i,k) &
	886	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
	887	temp = (den(i,k)*max(qci(i,k,2),qmin))
	888	temp = sqrt(sqrt(temptemptemp))
	889	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
	890	eacrs = exp(0.07*(-supcol))
	891	!
	892	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
	893	diameter = min(dicon * sqrt(xmi),dimax)
	894	vt2i = 1.49e4diameter*1.31
	895	vt2r=pvtrrslopeb(i,k,1)denfac(i,k)
	896	vt2s=pvtsrslopeb(i,k,2)denfac(i,k)
	897	vt2g=pvtgrslopeb(i,k,3)denfac(i,k)
	898	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
	899	if(qsum(i,k) .gt. 1.e-15) then
	900	vt2ave=(vt2sqrs(i,k,2)+vt2gqrs(i,k,3))/(qsum(i,k))
	901	else
	902	vt2ave=0.
	903	endif
	904	if(supcol.gt.0.and.qci(i,k,2).gt.qmin) then
	905	if(qrs(i,k,1).gt.qcrmin) then
	906	!-------------------------------------------------------------
	907	! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
	908	! (T<T0: I->R)
	909	!-------------------------------------------------------------
	910	acrfac = 2.rslope3(i,k,1)+2.diameter*rslope2(i,k,1) &
	911	+diameter*2rslope(i,k,1)
	912	praci(i,k) = piqci(i,k,2)n0rabs(vt2r-vt2i)acrfac/4.
	913	praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
	914	!-------------------------------------------------------------
	915	! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
	916	! (T<T0: R->S or R->G)
	917	!-------------------------------------------------------------
	918	piacr(i,k) = pi*2avtrn0rdenrxni(i,k)denfac(i,k) &
	919	g6pbrrslope3(i,k,1)*rslope3(i,k,1) &
	920	*rslopeb(i,k,1)/24./den(i,k)
	921	piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
	922	endif
	923	!-------------------------------------------------------------
	924	! psaci: Accretion of cloud ice by snow [HDC 10]
	925	! (T<T0: I->S)
	926	!-------------------------------------------------------------
	927	if(qrs(i,k,2).gt.qcrmin) then
	928	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
	929	+diameter*2rslope(i,k,2)
	930	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
	931	abs(vt2ave-vt2i)acrfac/4.
	932	psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
	933	endif
	934	!-------------------------------------------------------------
	935	! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
	936	! (T<T0: I->G)
	937	!-------------------------------------------------------------
	938	if(qrs(i,k,3).gt.qcrmin) then
	939	egi = exp(0.07*(-supcol))
	940	acrfac = 2.rslope3(i,k,3)+2.diameter*rslope2(i,k,3) &
	941	+diameter*2rslope(i,k,3)
	942	pgaci(i,k) = piegiqci(i,k,2)n0gabs(vt2ave-vt2i)*acrfac/4.
	943	pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
	944	endif
	945	endif
	946	!-------------------------------------------------------------
	947	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
	948	! (T<T0: C->S, and T>=T0: C->R)
	949	!-------------------------------------------------------------
	950	if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
	951	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2) &
	952	rslopeb(i,k,2)qci(i,k,1)*denfac(i,k) &
	953	,qci(i,k,1)/dtcld)
	954	endif
	955	!-------------------------------------------------------------
	956	! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
	957	! (T<T0: C->G, and T>=T0: C->R)
	958	!-------------------------------------------------------------
	959	if(qrs(i,k,3).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
	960	pgacw(i,k) = min(pacrgrslope3(i,k,3)rslopeb(i,k,3) &
	961	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
	962	endif
	963	!-------------------------------------------------------------
	964	! paacw: Accretion of cloud water by averaged snow/graupel
	965	! (T<T0: C->G or S, and T>=T0: C->R)
	966	!-------------------------------------------------------------
	967	if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,3).gt.qcrmin) then
	968	paacw(i,k) = (qrs(i,k,2)psacw(i,k)+qrs(i,k,3)pgacw(i,k))/(qsum(i,k))
	969	endif
	970	!-------------------------------------------------------------
	971	! pracs: Accretion of snow by rain [HL A11] [LFO 27]
	972	! (T<T0: S->G)
	973	!-------------------------------------------------------------
	974	if(qrs(i,k,2).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
	975	if(supcol.gt.0) then
	976	acrfac = 5.rslope3(i,k,2)rslope3(i,k,2)*rslope(i,k,1) &
	977	+2.rslope3(i,k,2)rslope2(i,k,2)*rslope2(i,k,1) &
	978	+.5rslope2(i,k,2)rslope2(i,k,2)*rslope3(i,k,1)
	979	pracs(i,k) = pi*2n0rn0sn0sfac(i,k)*abs(vt2r-vt2ave) &
	980	(dens/den(i,k))acrfac
	981	pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
	982	endif
	983	!-------------------------------------------------------------
	984	! psacr: Accretion of rain by snow [HL A10] [LFO 28]
	985	! (T<T0:R->S or R->G) (T>=T0: enhance melting of snow)
	986	!-------------------------------------------------------------
	987	acrfac = 5.rslope3(i,k,1)rslope3(i,k,1)*rslope(i,k,2) &
	988	+2.rslope3(i,k,1)rslope2(i,k,1)*rslope2(i,k,2) &
	989	+.5rslope2(i,k,1)rslope2(i,k,1)*rslope3(i,k,2)
	990	psacr(i,k) = pi*2n0rn0sn0sfac(i,k)*abs(vt2ave-vt2r) &
	991	(denr/den(i,k))acrfac
	992	psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
	993	endif
	994	!-------------------------------------------------------------
	995	! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
	996	! (T<T0: R->G) (T>=T0: enhance melting of graupel)
	997	!-------------------------------------------------------------
	998	if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,1).gt.qcrmin) then
	999	acrfac = 5.rslope3(i,k,1)rslope3(i,k,1)*rslope(i,k,3) &
	1000	+2.rslope3(i,k,1)rslope2(i,k,1)*rslope2(i,k,3) &
	1001	+.5rslope2(i,k,1)rslope2(i,k,1)*rslope3(i,k,3)
	1002	pgacr(i,k) = pi*2n0rn0gabs(vt2ave-vt2r)*(denr/den(i,k)) &
	1003	*acrfac
	1004	pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
	1005	endif
	1006	!
	1007	!-------------------------------------------------------------
	1008	! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
	1009	! (S->G): This process is eliminated in V3.0 with the
	1010	! new combined snow/graupel fall speeds
	1011	!-------------------------------------------------------------
	1012	! if(qrs(i,k,3).gt.qcrmin.and.qrs(i,k,2).gt.qcrmin) then
	1013	! acrfac = 5.rslope3(i,k,2)rslope3(i,k,2)*rslope(i,k,3) &
	1014	! +2.rslope3(i,k,2)rslope2(i,k,2)*rslope2(i,k,3) &
	1015	! +.5rslope2(i,k,2)rslope2(i,k,2)*rslope3(i,k,3)
	1016	! if(supcol.gt.0) then
	1017	! egs = exp(-0.09*supcol)
	1018	! else
	1019	! egs = 1.
	1020	! endif
	1021	! pgacs(i,k) = pi*2egsn0sn0sfac(i,k)n0gabs(vt2ave-vt2ave) &
	1022	! (dens/den(i,k))acrfac
	1023	pgacs(i,k) = 0.
	1024	! pgacs(i,k) = min(pgacs(i,k),qrs(i,k,2)/dtcld)
	1025	! endif
	1026	if(supcol.le.0) then
	1027	xlf = xlf0
	1028	!-------------------------------------------------------------
	1029	! pseml: Enhanced melting of snow by accretion of water [HL A34]
	1030	! (T>=T0: S->R)
	1031	!-------------------------------------------------------------
	1032	if(qrs(i,k,2).gt.0.) &
	1033	pseml(i,k) = min(max(cliqsupcol(paacw(i,k)+psacr(i,k)) &
	1034	/xlf,-qrs(i,k,2)/dtcld),0.)
	1035	!-------------------------------------------------------------
	1036	! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
	1037	! (T>=T0: G->R)
	1038	!-------------------------------------------------------------
	1039	if(qrs(i,k,3).gt.0.) &
	1040	pgeml(i,k) = min(max(cliqsupcol(paacw(i,k)+pgacr(i,k)) &
	1041	/xlf,-qrs(i,k,3)/dtcld),0.)
	1042	endif
	1043	if(supcol.gt.0) then
	1044	!-------------------------------------------------------------
	1045	! pidep: Deposition/Sublimation rate of ice [HDC 9]
	1046	! (T<T0: V->I or I->V)
	1047	!-------------------------------------------------------------
	1048	if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
	1049	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
	1050	supice = satdt-prevp(i,k)
	1051	if(pidep(i,k).lt.0.) then
	1052	pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
	1053	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
	1054	else
	1055	pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
	1056	endif
	1057	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
	1058	endif
	1059	!-------------------------------------------------------------
	1060	! psdep: deposition/sublimation rate of snow [HDC 14]
	1061	! (T<T0: V->S or S->V)
	1062	!-------------------------------------------------------------
	1063	if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
	1064	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
	1065	psdep(i,k) = (rh(i,k,2)-1.)n0sfac(i,k)(precs1 &
	1066	rslope2(i,k,2)+precs2work2(i,k) &
	1067	*coeres)/work1(i,k,2)
	1068	supice = satdt-prevp(i,k)-pidep(i,k)
	1069	if(psdep(i,k).lt.0.) then
	1070	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
	1071	psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
	1072	else
	1073	psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
	1074	endif
	1075	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
	1076	ifsat = 1
	1077	endif
	1078	!-------------------------------------------------------------
	1079	! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
	1080	! (T<T0: V->G or G->V)
	1081	!-------------------------------------------------------------
	1082	if(qrs(i,k,3).gt.0..and.ifsat.ne.1) then
	1083	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
	1084	pgdep(i,k) = (rh(i,k,2)-1.)(precg1rslope2(i,k,3) &
	1085	+precg2work2(i,k)coeres)/work1(i,k,2)
	1086	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
	1087	if(pgdep(i,k).lt.0.) then
	1088	pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
	1089	pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
	1090	else
	1091	pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
	1092	endif
	1093	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
	1094	abs(satdt)) ifsat = 1
	1095	endif
	1096	!-------------------------------------------------------------
	1097	! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
	1098	! (T<T0: V->I)
	1099	!-------------------------------------------------------------
	1100	if(supsat.gt.0.and.ifsat.ne.1) then
	1101	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
	1102	xni0 = 1.e3exp(0.1supcol)
	1103	roqi0 = 4.92e-11xni0*1.33
	1104	pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
	1105	/dtcld)
	1106	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
	1107	endif
	1108	!
	1109	!-------------------------------------------------------------
	1110	! psaut: conversion(aggregation) of ice to snow [HDC 12]
	1111	! (T<T0: I->S)
	1112	!-------------------------------------------------------------
	1113	if(qci(i,k,2).gt.0.) then
	1114	qimax = roqimax/den(i,k)
	1115	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
	1116	endif
	1117	!
	1118	!-------------------------------------------------------------
	1119	! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
	1120	! (T<T0: S->G)
	1121	!-------------------------------------------------------------
	1122	if(qrs(i,k,2).gt.0.) then
	1123	alpha2 = 1.e-3exp(0.09(-supcol))
	1124	pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)) &
	1125	,qrs(i,k,2)/dtcld)
	1126	endif
	1127	endif
	1128	!
	1129	!-------------------------------------------------------------
	1130	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
	1131	! (T>=T0: S->V)
	1132	!-------------------------------------------------------------
	1133	if(supcol.lt.0.) then
	1134	if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) then
	1135	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
	1136	psevp(i,k) = (rh(i,k,1)-1.)n0sfac(i,k)(precs1 &
	1137	rslope2(i,k,2)+precs2work2(i,k) &
	1138	*coeres)/work1(i,k,1)
	1139	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
	1140	endif
	1141	!-------------------------------------------------------------
	1142	! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
	1143	! (T>=T0: G->V)
	1144	!-------------------------------------------------------------
	1145	if(qrs(i,k,3).gt.0..and.rh(i,k,1).lt.1.) then
	1146	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
	1147	pgevp(i,k) = (rh(i,k,1)-1.)(precg1rslope2(i,k,3) &
	1148	+precg2work2(i,k)coeres)/work1(i,k,1)
	1149	pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
	1150	endif
	1151	endif
	1152	enddo
	1153	enddo
	1154	!
	1155	!
	1156	!----------------------------------------------------------------
	1157	! check mass conservation of generation terms and feedback to the
	1158	! large scale
	1159	!
	1160	do k = kts, kte
	1161	do i = its, ite
	1162	!
	1163	delta2=0.
	1164	delta3=0.
	1165	if(qrs(i,k,1).lt.1.e-4.and.qrs(i,k,2).lt.1.e-4) delta2=1.
	1166	if(qrs(i,k,1).lt.1.e-4) delta3=1.
	1167	if(t(i,k).le.t0c) then
	1168	!
	1169	! cloud water
	1170	!
	1171	value = max(qmin,qci(i,k,1))
	1172	source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
	1173	if (source.gt.value) then
	1174	factor = value/source
	1175	praut(i,k) = praut(i,k)*factor
	1176	pracw(i,k) = pracw(i,k)*factor
	1177	paacw(i,k) = paacw(i,k)*factor
	1178	endif
	1179	!
	1180	! cloud ice
	1181	!
	1182	value = max(qmin,qci(i,k,2))
	1183	source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k) &
	1184	+psaci(i,k)+pgaci(i,k))*dtcld
	1185	if (source.gt.value) then
	1186	factor = value/source
	1187	psaut(i,k) = psaut(i,k)*factor
	1188	pigen(i,k) = pigen(i,k)*factor
	1189	pidep(i,k) = pidep(i,k)*factor
	1190	praci(i,k) = praci(i,k)*factor
	1191	psaci(i,k) = psaci(i,k)*factor
	1192	pgaci(i,k) = pgaci(i,k)*factor
	1193	endif
	1194	!
	1195	! rain
	1196	!
	1197	value = max(qmin,qrs(i,k,1))
	1198	source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) &
	1199	+psacr(i,k)+pgacr(i,k))*dtcld
	1200	if (source.gt.value) then
	1201	factor = value/source
	1202	praut(i,k) = praut(i,k)*factor
	1203	prevp(i,k) = prevp(i,k)*factor
	1204	pracw(i,k) = pracw(i,k)*factor
	1205	piacr(i,k) = piacr(i,k)*factor
	1206	psacr(i,k) = psacr(i,k)*factor
	1207	pgacr(i,k) = pgacr(i,k)*factor
	1208	endif
	1209	!
	1210	! snow
	1211	!
	1212	value = max(qmin,qrs(i,k,2))
	1213	source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) &
	1214	+piacr(i,k)delta3+praci(i,k)delta3 &
	1215	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2 &
	1216	+psaci(i,k)-pgacs(i,k) )*dtcld
	1217	if (source.gt.value) then
	1218	factor = value/source
	1219	psdep(i,k) = psdep(i,k)*factor
	1220	psaut(i,k) = psaut(i,k)*factor
	1221	pgaut(i,k) = pgaut(i,k)*factor
	1222	paacw(i,k) = paacw(i,k)*factor
	1223	piacr(i,k) = piacr(i,k)*factor
	1224	praci(i,k) = praci(i,k)*factor
	1225	psaci(i,k) = psaci(i,k)*factor
	1226	pracs(i,k) = pracs(i,k)*factor
	1227	psacr(i,k) = psacr(i,k)*factor
	1228	pgacs(i,k) = pgacs(i,k)*factor
	1229	endif
	1230	!
	1231	! graupel
	1232	!
	1233	value = max(qmin,qrs(i,k,3))
	1234	source = -(pgdep(i,k)+pgaut(i,k) &
	1235	+piacr(i,k)(1.-delta3)+praci(i,k)(1.-delta3) &
	1236	+psacr(i,k)(1.-delta2)+pracs(i,k)(1.-delta2) &
	1237	+pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
	1238	if (source.gt.value) then
	1239	factor = value/source
	1240	pgdep(i,k) = pgdep(i,k)*factor
	1241	pgaut(i,k) = pgaut(i,k)*factor
	1242	piacr(i,k) = piacr(i,k)*factor
	1243	praci(i,k) = praci(i,k)*factor
	1244	psacr(i,k) = psacr(i,k)*factor
	1245	pracs(i,k) = pracs(i,k)*factor
	1246	paacw(i,k) = paacw(i,k)*factor
	1247	pgaci(i,k) = pgaci(i,k)*factor
	1248	pgacr(i,k) = pgacr(i,k)*factor
	1249	pgacs(i,k) = pgacs(i,k)*factor
	1250	endif
	1251	!
	1252	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k) &
	1253	+pidep(i,k))
	1254	! update
	1255	q(i,k) = q(i,k)+work2(i,k)*dtcld
	1256	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
	1257	+paacw(i,k)+paacw(i,k))*dtcld,0.)
	1258	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
	1259	+prevp(i,k)-piacr(i,k)-pgacr(i,k) &
	1260	-psacr(i,k))*dtcld,0.)
	1261	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
	1262	+psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
	1263	*dtcld,0.)
	1264	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
	1265	-pgaut(i,k)+piacr(i,k)*delta3 &
	1266	+praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
	1267	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2) &
	1268	*dtcld,0.)
	1269	qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
	1270	+piacr(i,k)*(1.-delta3) &
	1271	+praci(i,k)(1.-delta3)+psacr(i,k)(1.-delta2)&
	1272	+pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
	1273	+pgacr(i,k)+pgacs(i,k))*dtcld,0.)
	1274	xlf = xls-xl(i,k)
	1275	xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
	1276	-xl(i,k)prevp(i,k)-xlf(piacr(i,k)+paacw(i,k) &
	1277	+paacw(i,k)+pgacr(i,k)+psacr(i,k))
	1278	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
	1279	else
	1280	!
	1281	! cloud water
	1282	!
	1283	value = max(qmin,qci(i,k,1))
	1284	source=(praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))*dtcld
	1285	if (source.gt.value) then
	1286	factor = value/source
	1287	praut(i,k) = praut(i,k)*factor
	1288	pracw(i,k) = pracw(i,k)*factor
	1289	paacw(i,k) = paacw(i,k)*factor
	1290	endif
	1291	!
	1292	! rain
	1293	!
	1294	value = max(qmin,qrs(i,k,1))
	1295	source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) &
	1296	-pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
	1297	if (source.gt.value) then
	1298	factor = value/source
	1299	praut(i,k) = praut(i,k)*factor
	1300	prevp(i,k) = prevp(i,k)*factor
	1301	pracw(i,k) = pracw(i,k)*factor
	1302	paacw(i,k) = paacw(i,k)*factor
	1303	pseml(i,k) = pseml(i,k)*factor
	1304	pgeml(i,k) = pgeml(i,k)*factor
	1305	endif
	1306	!
	1307	! snow
	1308	!
	1309	value = max(qcrmin,qrs(i,k,2))
	1310	source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
	1311	if (source.gt.value) then
	1312	factor = value/source
	1313	pgacs(i,k) = pgacs(i,k)*factor
	1314	psevp(i,k) = psevp(i,k)*factor
	1315	pseml(i,k) = pseml(i,k)*factor
	1316	endif
	1317	!
	1318	! graupel
	1319	!
	1320	value = max(qcrmin,qrs(i,k,3))
	1321	source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
	1322	if (source.gt.value) then
	1323	factor = value/source
	1324	pgacs(i,k) = pgacs(i,k)*factor
	1325	pgevp(i,k) = pgevp(i,k)*factor
	1326	pgeml(i,k) = pgeml(i,k)*factor
	1327	endif
	1328	work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
	1329	! update
	1330	q(i,k) = q(i,k)+work2(i,k)*dtcld
	1331	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
	1332	+paacw(i,k)+paacw(i,k))*dtcld,0.)
	1333	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
	1334	+prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
	1335	-pgeml(i,k))*dtcld,0.)
	1336	qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
	1337	+pseml(i,k))*dtcld,0.)
	1338	qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
	1339	+pgeml(i,k))*dtcld,0.)
	1340	xlf = xls-xl(i,k)
	1341	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
	1342	-xlf*(pseml(i,k)+pgeml(i,k))
	1343	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
	1344	endif
	1345	enddo
	1346	enddo
	1347	!
	1348	! Inline expansion for fpvs
	1349	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	1350	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	1351	hsub = xls
	1352	hvap = xlv0
	1353	cvap = cpv
	1354	ttp=t0c+0.01
	1355	dldt=cvap-cliq
	1356	xa=-dldt/rv
	1357	xb=xa+hvap/(rv*ttp)
	1358	dldti=cvap-cice
	1359	xai=-dldti/rv
	1360	xbi=xai+hsub/(rv*ttp)
	1361	do k = kts, kte
	1362	do i = its, ite
	1363	tr=ttp/t(i,k)
	1364	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	1365	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
	1366	qs(i,k,1) = max(qs(i,k,1),qmin)
	1367	tr=ttp/t(i,k)
	1368	if(t(i,k).lt.ttp) then
	1369	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
	1370	else
	1371	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	1372	endif
	1373	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
	1374	qs(i,k,2) = max(qs(i,k,2),qmin)
	1375	enddo
	1376	enddo
	1377	!
	1378	!----------------------------------------------------------------
	1379	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
	1380	! if there exists additional water vapor condensated/if
	1381	! evaporation of cloud water is not enough to remove subsaturation
	1382	!
	1383	do k = kts, kte
	1384	do i = its, ite
	1385	work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
	1386	work2(i,k) = qci(i,k,1)+work1(i,k,1)
	1387	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
	1388	if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
	1389	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
	1390	q(i,k) = q(i,k)-pcond(i,k)*dtcld
	1391	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
	1392	t(i,k) = t(i,k)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
	1393	enddo
	1394	enddo
	1395	!
	1396	!
	1397	!----------------------------------------------------------------
	1398	! padding for small values
	1399	!
	1400	do k = kts, kte
	1401	do i = its, ite
	1402	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
	1403	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
	1404	enddo
	1405	enddo
	1406	enddo ! big loops
	1407	END SUBROUTINE wsm62d
	1408	! ...................................................................
	1409	REAL FUNCTION rgmma(x)
	1410	!-------------------------------------------------------------------
	1411	IMPLICIT NONE
	1412	!-------------------------------------------------------------------
	1413	! rgmma function: use infinite product form
	1414	REAL :: euler
	1415	PARAMETER (euler=0.577215664901532)
	1416	REAL :: x, y
	1417	INTEGER :: i
	1418	if(x.eq.1.)then
	1419	rgmma=0.
	1420	else
	1421	rgmma=xexp(eulerx)
	1422	do i=1,10000
	1423	y=float(i)
	1424	rgmma=rgmma(1.000+x/y)exp(-x/y)
	1425	enddo
	1426	rgmma=1./rgmma
	1427	endif
	1428	END FUNCTION rgmma
	1429	!
	1430	!--------------------------------------------------------------------------
	1431	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
	1432	!--------------------------------------------------------------------------
	1433	IMPLICIT NONE
	1434	!--------------------------------------------------------------------------
	1435	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
	1436	xai,xbi,ttp,tr
	1437	INTEGER ice
	1438	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1439	ttp=t0c+0.01
	1440	dldt=cvap-cliq
	1441	xa=-dldt/rv
	1442	xb=xa+hvap/(rv*ttp)
	1443	dldti=cvap-cice
	1444	xai=-dldti/rv
	1445	xbi=xai+hsub/(rv*ttp)
	1446	tr=ttp/t
	1447	if(t.lt.ttp.and.ice.eq.1) then
	1448	fpvs=psat(trxai)exp(xbi*(1.-tr))
	1449	else
	1450	fpvs=psat(trxa)exp(xb*(1.-tr))
	1451	endif
	1452	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1453	END FUNCTION fpvs
	1454	!-------------------------------------------------------------------
	1455	SUBROUTINE wsm6init(den0,denr,dens,cl,cpv,allowed_to_read)
	1456	!-------------------------------------------------------------------
	1457	IMPLICIT NONE
	1458	!-------------------------------------------------------------------
	1459	!.... constants which may not be tunable
	1460	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
	1461	LOGICAL, INTENT(IN) :: allowed_to_read
	1462	REAL :: pi
	1463	!
	1464	pi = 4.*atan(1.)
	1465	xlv1 = cl-cpv
	1466	!
	1467	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
	1468	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
	1469	!
	1470	bvtr1 = 1.+bvtr
	1471	bvtr2 = 2.5+.5*bvtr
	1472	bvtr3 = 3.+bvtr
	1473	bvtr4 = 4.+bvtr
	1474	bvtr6 = 6.+bvtr
	1475	g1pbr = rgmma(bvtr1)
	1476	g3pbr = rgmma(bvtr3)
	1477	g4pbr = rgmma(bvtr4) ! 17.837825
	1478	g6pbr = rgmma(bvtr6)
	1479	g5pbro2 = rgmma(bvtr2) ! 1.8273
	1480	pvtr = avtr*g4pbr/6.
	1481	eacrr = 1.0
	1482	pacrr = pin0ravtrg3pbr.25*eacrr
	1483	precr1 = 2.pin0r*.78
	1484	precr2 = 2.pin0r.31avtr*.5g5pbro2
	1485	xm0 = (di0/dicon)**2
	1486	xmmax = (dimax/dicon)**2
	1487	roqimax = 2.08e22dimax*8
	1488	!
	1489	bvts1 = 1.+bvts
	1490	bvts2 = 2.5+.5*bvts
	1491	bvts3 = 3.+bvts
	1492	bvts4 = 4.+bvts
	1493	g1pbs = rgmma(bvts1) !.8875
	1494	g3pbs = rgmma(bvts3)
	1495	g4pbs = rgmma(bvts4) ! 12.0786
	1496	g5pbso2 = rgmma(bvts2)
	1497	pvts = avts*g4pbs/6.
	1498	pacrs = pin0savtsg3pbs.25
	1499	precs1 = 4.n0s.65
	1500	precs2 = 4.n0s.44avts.5g5pbso2
	1501	pidn0r = pidenrn0r
	1502	pidn0s = pidensn0s
	1503	!
	1504	pacrc = pin0savtsg3pbs.25*eacrc
	1505	!
	1506	bvtg1 = 1.+bvtg
	1507	bvtg2 = 2.5+.5*bvtg
	1508	bvtg3 = 3.+bvtg
	1509	bvtg4 = 4.+bvtg
	1510	g1pbg = rgmma(bvtg1)
	1511	g3pbg = rgmma(bvtg3)
	1512	g4pbg = rgmma(bvtg4)
	1513	pacrg = pin0gavtgg3pbg.25
	1514	g5pbgo2 = rgmma(bvtg2)
	1515	pvtg = avtg*g4pbg/6.
	1516	precg1 = 2.pin0g*.78
	1517	precg2 = 2.pin0g.31avtg*.5g5pbgo2
	1518	pidn0g = pidengn0g
	1519	!
	1520	rslopermax = 1./lamdarmax
	1521	rslopesmax = 1./lamdasmax
	1522	rslopegmax = 1./lamdagmax
	1523	rsloperbmax = rslopermax ** bvtr
	1524	rslopesbmax = rslopesmax ** bvts
	1525	rslopegbmax = rslopegmax ** bvtg
	1526	rsloper2max = rslopermax * rslopermax
	1527	rslopes2max = rslopesmax * rslopesmax
	1528	rslopeg2max = rslopegmax * rslopegmax
	1529	rsloper3max = rsloper2max * rslopermax
	1530	rslopes3max = rslopes2max * rslopesmax
	1531	rslopeg3max = rslopeg2max * rslopegmax
	1532	!
	1533	END SUBROUTINE wsm6init
	1534	END MODULE module_mp_wsm6

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

Download in other formats:

Original Format