Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_wdm6.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 110.9 KB

Rev	Line
[1]	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_wdm6
	10	!
	11	!
	12	!
	13	REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
	14	REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
	15	REAL, PARAMETER, PRIVATE :: n0g = 4.e6 ! intercept parameter graupel
	16	REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
	17	REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
	18	REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
	19	REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
	20	REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
	21	REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
	22	REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
	23	REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
	24	REAL, PARAMETER, PRIVATE :: avtg = 330. ! a constant for terminal velocity of graupel
	25	REAL, PARAMETER, PRIVATE :: bvtg = 0.8 ! a constant for terminal velocity of graupel
	26	REAL, PARAMETER, PRIVATE :: deng = 500. ! density of graupel
	27	REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
	28	REAL, PARAMETER, PRIVATE :: lamdacmax = 1.e10 ! limited maximum value for slope parameter of cloud water
	29	REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e8 ! limited maximum value for slope parameter of rain
	30	REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
	31	REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
	32	REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
	33	REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
	34	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
	35	REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
	36	REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
	37	REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
	38	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
	39	REAL, PARAMETER, PRIVATE :: ncmin = 1.e1 ! minimum value for Nc
	40	REAL, PARAMETER, PRIVATE :: nrmin = 1.e-2 ! minimum value for Nr
	41	REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
	42	REAL, PARAMETER, PRIVATE :: dens = 100.0 ! Density of snow
	43	REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! threshold amount for aggretion to occur
	44	!
	45	REAL, PARAMETER, PRIVATE :: satmax = 1.0048 ! maximum saturation value for CCN activation
	46	! 1.008 for maritime /1.0048 for conti
	47	REAL, PARAMETER, PRIVATE :: actk = 0.6 ! parameter for the CCN activation
	48	REAL, PARAMETER, PRIVATE :: actr = 1.5 ! radius of activated CCN drops
	49	REAL, PARAMETER, PRIVATE :: ncrk1 = 3.03e3 ! Long's collection kernel coefficient
	50	REAL, PARAMETER, PRIVATE :: ncrk2 = 2.59e15 ! Long's collection kernel coefficient
	51	REAL, PARAMETER, PRIVATE :: di100 = 1.e-4 ! parameter related with accretion and collection of cloud drops
	52	REAL, PARAMETER, PRIVATE :: di600 = 6.e-4 ! parameter related with accretion and collection of cloud drops
	53	REAL, PARAMETER, PRIVATE :: di2000 = 2000.e-6 ! parameter related with accretion and collection of cloud drops
	54	REAL, PARAMETER, PRIVATE :: di82 = 82.e-6 ! dimater related with raindrops evaporation
	55	REAL, PARAMETER, PRIVATE :: di15 = 15.e-6 ! auto conversion takes place beyond this diameter
	56	!
	57	REAL, SAVE :: &
	58	qc0,qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4,bvtr5, &
	59	bvtr6,bvtr7, bvtr2o5,bvtr3o5, &
	60	g1pbr,g2pbr,g3pbr,g4pbr,g5pbr,g6pbr,g7pbr, &
	61	g5pbro2,g7pbro2,pi, &
	62	pvtr,pvtrn,eacrr,pacrr,pidn0r,pidnr, &
	63	precr1,precr2,xmmax,roqimax,bvts1,bvts2, &
	64	bvts3,bvts4,g1pbs,g3pbs,g4pbs,g5pbso2, &
	65	pvts,pacrs,precs1,precs2,pidn0s,xlv1,pacrc, &
	66	bvtg1,bvtg2,bvtg3,bvtg4,g1pbg,g3pbg,g4pbg, &
	67	g5pbgo2,pvtg,pacrg,precg1,precg2,pidn0g, &
	68	rslopecmax,rslopec2max,rslopec3max, &
	69	rslopermax,rslopesmax,rslopegmax, &
	70	rsloperbmax,rslopesbmax,rslopegbmax, &
	71	rsloper2max,rslopes2max,rslopeg2max, &
	72	rsloper3max,rslopes3max,rslopeg3max
	73	CONTAINS
	74	!===================================================================
	75	!
	76	SUBROUTINE wdm6(th, q, qc, qr, qi, qs, qg, &
	77	nn, nc, nr, &
	78	den, pii, p, delz, &
	79	delt,g, cpd, cpv, ccn0, rd, rv, t0c, &
	80	ep1, ep2, qmin, &
	81	XLS, XLV0, XLF0, den0, denr, &
	82	cliq,cice,psat, &
	83	rain, rainncv, &
	84	snow, snowncv, &
	85	sr, &
	86	graupel, graupelncv, &
	87	itimestep, &
	88	ids,ide, jds,jde, kds,kde, &
	89	ims,ime, jms,jme, kms,kme, &
	90	its,ite, jts,jte, kts,kte &
	91	)
	92	!-------------------------------------------------------------------
	93	IMPLICIT NONE
	94	!-------------------------------------------------------------------
	95	!
	96	! This code is a WRF double-moment 6-class GRAUPEL phase
	97	! microphyiscs scheme (WDM6). The WDM microphysics scheme predicts
	98	! number concentrations for warm rain species including clouds and
	99	! rain. cloud condensation nuclei (CCN) is also predicted.
	100	! The cold rain species including ice, snow, graupel follow the
	101	! WRF single-moment 6-class microphysics (WSM6, Hong and Lim 2006)
	102	! in which theoretical background for WSM ice phase microphysics is
	103	! based on Hong et al. (2004). A new mixed-phase terminal velocity
	104	! for precipitating ice is introduced in WSM6 (Dudhia et al. 2008).
	105	! The WDM scheme is described in Lim and Hong (2010).
	106	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
	107	!
	108	! WDM6 cloud scheme
	109	!
	110	! Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
	111	!
	112	! Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
	113	!
	114	! Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev.
	115	! Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
	116	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
	117	! Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
	118	! Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
	119	! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
	120	!
	121	! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
	122	! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
	123	! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
	124	!
	125	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
	126	ims,ime, jms,jme, kms,kme , &
	127	its,ite, jts,jte, kts,kte
	128	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	129	INTENT(INOUT) :: &
	130	th, &
	131	q, &
	132	qc, &
	133	qi, &
	134	qr, &
	135	qs, &
	136	qg, &
	137	nn, &
	138	nc, &
	139	nr
	140	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	141	INTENT(IN ) :: &
	142	den, &
	143	pii, &
	144	p, &
	145	delz
	146	REAL, INTENT(IN ) :: delt, &
	147	g, &
	148	rd, &
	149	rv, &
	150	t0c, &
	151	den0, &
	152	cpd, &
	153	cpv, &
	154	ccn0, &
	155	ep1, &
	156	ep2, &
	157	qmin, &
	158	XLS, &
	159	XLV0, &
	160	XLF0, &
	161	cliq, &
	162	cice, &
	163	psat, &
	164	denr
	165	INTEGER, INTENT(IN ) :: itimestep
	166	REAL, DIMENSION( ims:ime , jms:jme ), &
	167	INTENT(INOUT) :: rain, &
	168	rainncv, &
	169	sr
	170	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
	171	INTENT(INOUT) :: snow, &
	172	snowncv
	173	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
	174	INTENT(INOUT) :: graupel, &
	175	graupelncv
	176	! LOCAL VAR
	177	REAL, DIMENSION( its:ite , kts:kte ) :: t
	178	REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
	179	REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs, ncr
	180	INTEGER :: i,j,k
	181	!-------------------------------------------------------------------
	182	IF (itimestep .eq. 1) THEN
	183	DO j=jms,jme
	184	DO k=kms,kme
	185	DO i=ims,ime
	186	nn(i,k,j) = ccn0
	187	ENDDO
	188	ENDDO
	189	ENDDO
	190	ENDIF
	191	!
	192	DO j=jts,jte
	193	DO k=kts,kte
	194	DO i=its,ite
	195	t(i,k)=th(i,k,j)*pii(i,k,j)
	196	qci(i,k,1) = qc(i,k,j)
	197	qci(i,k,2) = qi(i,k,j)
	198	qrs(i,k,1) = qr(i,k,j)
	199	qrs(i,k,2) = qs(i,k,j)
	200	qrs(i,k,3) = qg(i,k,j)
	201	ncr(i,k,1) = nn(i,k,j)
	202	ncr(i,k,2) = nc(i,k,j)
	203	ncr(i,k,3) = nr(i,k,j)
	204	ENDDO
	205	ENDDO
	206	! Sending array starting locations of optional variables may cause
	207	! troubles, so we explicitly change the call.
	208	CALL wdm62D(t, q(ims,kms,j), qci, qrs, ncr &
	209	,den(ims,kms,j) &
	210	,p(ims,kms,j), delz(ims,kms,j) &
	211	,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
	212	,ep1, ep2, qmin &
	213	,XLS, XLV0, XLF0, den0, denr &
	214	,cliq,cice,psat &
	215	,j &
	216	,rain(ims,j),rainncv(ims,j) &
	217	,sr(ims,j) &
	218	,ids,ide, jds,jde, kds,kde &
	219	,ims,ime, jms,jme, kms,kme &
	220	,its,ite, jts,jte, kts,kte &
	221	,snow(ims,j),snowncv(ims,j) &
	222	,graupel(ims,j),graupelncv(ims,j) &
	223	)
	224	DO K=kts,kte
	225	DO I=its,ite
	226	th(i,k,j)=t(i,k)/pii(i,k,j)
	227	qc(i,k,j) = qci(i,k,1)
	228	qi(i,k,j) = qci(i,k,2)
	229	qr(i,k,j) = qrs(i,k,1)
	230	qs(i,k,j) = qrs(i,k,2)
	231	qg(i,k,j) = qrs(i,k,3)
	232	nn(i,k,j) = ncr(i,k,1)
	233	nc(i,k,j) = ncr(i,k,2)
	234	nr(i,k,j) = ncr(i,k,3)
	235	ENDDO
	236	ENDDO
	237	ENDDO
	238	END SUBROUTINE wdm6
	239	!===================================================================
	240	!
	241	SUBROUTINE wdm62D(t, q, qci, qrs, ncr, den, p, delz &
	242	,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
	243	,ep1, ep2, qmin &
	244	,XLS, XLV0, XLF0, den0, denr &
	245	,cliq,cice,psat &
	246	,lat &
	247	,rain,rainncv &
	248	,sr &
	249	,ids,ide, jds,jde, kds,kde &
	250	,ims,ime, jms,jme, kms,kme &
	251	,its,ite, jts,jte, kts,kte &
	252	,snow,snowncv &
	253	,graupel,graupelncv &
	254	)
	255	!-------------------------------------------------------------------
	256	IMPLICIT NONE
	257	!-------------------------------------------------------------------
	258	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
	259	ims,ime, jms,jme, kms,kme , &
	260	its,ite, jts,jte, kts,kte , &
	261	lat
	262	REAL, DIMENSION( its:ite , kts:kte ), &
	263	INTENT(INOUT) :: &
	264	t
	265	REAL, DIMENSION( its:ite , kts:kte, 2 ), &
	266	INTENT(INOUT) :: &
	267	qci
	268	REAL, DIMENSION( its:ite , kts:kte, 3 ), &
	269	INTENT(INOUT) :: &
	270	qrs, &
	271	ncr
	272	REAL, DIMENSION( ims:ime , kms:kme ), &
	273	INTENT(INOUT) :: &
	274	q
	275	REAL, DIMENSION( ims:ime , kms:kme ), &
	276	INTENT(IN ) :: &
	277	den, &
	278	p, &
	279	delz
	280	REAL, INTENT(IN ) :: delt, &
	281	g, &
	282	cpd, &
	283	cpv, &
	284	ccn0, &
	285	t0c, &
	286	den0, &
	287	rd, &
	288	rv, &
	289	ep1, &
	290	ep2, &
	291	qmin, &
	292	XLS, &
	293	XLV0, &
	294	XLF0, &
	295	cliq, &
	296	cice, &
	297	psat, &
	298	denr
	299	REAL, DIMENSION( ims:ime ), &
	300	INTENT(INOUT) :: rain, &
	301	rainncv, &
	302	sr
	303	REAL, DIMENSION( ims:ime ), OPTIONAL, &
	304	INTENT(INOUT) :: snow, &
	305	snowncv
	306	REAL, DIMENSION( ims:ime ), OPTIONAL, &
	307	INTENT(INOUT) :: graupel, &
	308	graupelncv
	309	! LOCAL VAR
	310	REAL, DIMENSION( its:ite , kts:kte , 3) :: &
	311	rh, qs, rslope, rslope2, rslope3, rslopeb, &
	312	falk, fall, work1, qrs_tmp
	313	REAL, DIMENSION( its:ite , kts:kte ) :: &
	314	rslopec, rslopec2,rslopec3
	315	REAL, DIMENSION( its:ite , kts:kte, 2) :: &
	316	avedia
	317	REAL, DIMENSION( its:ite , kts:kte ) :: &
	318	workn,falln,falkn
	319	REAL, DIMENSION( its:ite , kts:kte ) :: &
	320	worka,workr
	321	REAL, DIMENSION( its:ite , kts:kte ) :: &
	322	den_tmp, delz_tmp, ncr_tmp
	323	REAL, DIMENSION( its:ite , kts:kte ) :: &
	324	falkc, work1c, work2c, fallc
	325	REAL, DIMENSION( its:ite , kts:kte ) :: &
	326	pcact, prevp, psdep, pgdep, praut, psaut, pgaut, &
	327	pracw, psacw, pgacw, pgacr, pgacs, psaci, pgmlt, praci, &
	328	piacr, pracs, psacr, pgaci, pseml, pgeml
	329	REAL, DIMENSION( its:ite , kts:kte ) :: paacw
	330	REAL, DIMENSION( its:ite , kts:kte ) :: &
	331	nraut, nracw, ncevp, nccol, nrcol, &
	332	nsacw, ngacw, niacr, nsacr, ngacr, naacw, &
	333	nseml, ngeml, ncact
	334	REAL, DIMENSION( its:ite , kts:kte ) :: &
	335	pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, &
	336	denfac, xni, pgevp,n0sfac, qsum, &
	337	denqrs1, denqr1, denqrs2, denqrs3, denncr3, denqci
	338	REAL, DIMENSION( its:ite ) :: &
	339	delqrs1, delqrs2, delqrs3, delncr3, delqi
	340	REAL, DIMENSION( its:ite ) :: tstepsnow, tstepgraup
	341	REAL :: gfac, sfac
	342	! variables for optimization
	343	REAL, DIMENSION( its:ite ) :: tvec1
	344	REAL :: temp
	345	INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
	346	INTEGER, DIMENSION( its:ite ) :: mstep, numdt
	347	LOGICAL, DIMENSION( its:ite ) :: flgcld
	348	REAL :: &
	349	cpmcal, xlcal, lamdac, &
	350	diffus, &
	351	viscos, xka, venfac, conden, diffac, &
	352	x, y, z, a, b, c, d, e, &
	353	ndt, qdt, holdrr, holdrs, holdrg, supcol, supcolt, &
	354	pvt, coeres, supsat, dtcld, xmi, eacrs, satdt, &
	355	qimax, diameter, xni0, roqi0, &
	356	fallsum, fallsum_qsi, fallsum_qg, &
	357	vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
	358	xlwork2, factor, source, value, coecol, &
	359	nfrzdtr, nfrzdtc, &
	360	taucon, lencon, lenconcr, &
	361	xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
	362	REAL :: vt2ave
	363	REAL :: holdc, holdci
	364	!
	365	INTEGER :: i, j, k, mstepmax, &
	366	iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
	367	! Temporaries used for inlining fpvs function
	368	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
	369	!
	370	!=================================================================
	371	! compute internal functions
	372	!
	373	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
	374	xlcal(x) = xlv0-xlv1*(x-t0c)
	375	!----------------------------------------------------------------
	376	! size distributions: (x=mixing ratio, y=air density):
	377	! valid for mixing ratio > 1.e-9 kg/kg.
	378	!
	379	! Optimizatin : A*B => exp(log(A)(B))
	380	lamdac(x,y,z)= exp(log(((pidncz)/(xy)))*((.33333333)))
	381	!----------------------------------------------------------------
	382	! diffus: diffusion coefficient of the water vapor
	383	! viscos: kinematic viscosity(m2s-1)
	384	!
	385	diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
	386	viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
	387	xka(x,y) = 1.414e3viscos(x,y)y
	388	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
	389	venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
	390	/sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
	391	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
	392	!
	393	idim = ite-its+1
	394	kdim = kte-kts+1
	395	!
	396	!----------------------------------------------------------------
	397	! paddint 0 for negative values generated by dynamics
	398	!
	399	do k = kts, kte
	400	do i = its, ite
	401	qci(i,k,1) = max(qci(i,k,1),0.0)
	402	qrs(i,k,1) = max(qrs(i,k,1),0.0)
	403	qci(i,k,2) = max(qci(i,k,2),0.0)
	404	qrs(i,k,2) = max(qrs(i,k,2),0.0)
	405	qrs(i,k,3) = max(qrs(i,k,3),0.0)
	406	ncr(i,k,1) = max(ncr(i,k,1),0.0)
	407	ncr(i,k,2) = max(ncr(i,k,2),0.0)
	408	ncr(i,k,3) = max(ncr(i,k,3),0.0)
	409	enddo
	410	enddo
	411	!
	412	!----------------------------------------------------------------
	413	! latent heat for phase changes and heat capacity. neglect the
	414	! changes during microphysical process calculation
	415	! emanuel(1994)
	416	!
	417	do k = kts, kte
	418	do i = its, ite
	419	cpm(i,k) = cpmcal(q(i,k))
	420	xl(i,k) = xlcal(t(i,k))
	421	enddo
	422	enddo
	423	do k = kts, kte
	424	do i = its, ite
	425	delz_tmp(i,k) = delz(i,k)
	426	den_tmp(i,k) = den(i,k)
	427	enddo
	428	enddo
	429	!
	430	!----------------------------------------------------------------
	431	! initialize the surface rain, snow, graupel
	432	!
	433	do i = its, ite
	434	rainncv(i) = 0.
	435	if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i) = 0.
	436	if(PRESENT (graupelncv) .AND. PRESENT (graupel)) graupelncv(i) = 0.
	437	sr(i) = 0.
	438	! new local array to catch step snow and graupel
	439	tstepsnow(i) = 0.
	440	tstepgraup(i) = 0.
	441	enddo
	442	!
	443	!----------------------------------------------------------------
	444	! compute the minor time steps.
	445	!
	446	loops = max(nint(delt/dtcldcr),1)
	447	dtcld = delt/loops
	448	if(delt.le.dtcldcr) dtcld = delt
	449	!
	450	do loop = 1,loops
	451	!
	452	!----------------------------------------------------------------
	453	! initialize the large scale variables
	454	!
	455	do i = its, ite
	456	mstep(i) = 1
	457	mnstep(i) = 1
	458	flgcld(i) = .true.
	459	enddo
	460	!
	461	do k = kts, kte
	462	CALL VREC( tvec1(its), den(its,k), ite-its+1)
	463	do i = its, ite
	464	tvec1(i) = tvec1(i)*den0
	465	enddo
	466	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
	467	enddo
	468	!
	469	! Inline expansion for fpvs
	470	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	471	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	472	hsub = xls
	473	hvap = xlv0
	474	cvap = cpv
	475	ttp=t0c+0.01
	476	dldt=cvap-cliq
	477	xa=-dldt/rv
	478	xb=xa+hvap/(rv*ttp)
	479	dldti=cvap-cice
	480	xai=-dldti/rv
	481	xbi=xai+hsub/(rv*ttp)
	482	do k = kts, kte
	483	do i = its, ite
	484	tr=ttp/t(i,k)
	485	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	486	qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
	487	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
	488	qs(i,k,1) = max(qs(i,k,1),qmin)
	489	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
	490	tr=ttp/t(i,k)
	491	if(t(i,k).lt.ttp) then
	492	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
	493	else
	494	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	495	endif
	496	qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
	497	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
	498	qs(i,k,2) = max(qs(i,k,2),qmin)
	499	rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
	500	enddo
	501	enddo
	502	!
	503	!----------------------------------------------------------------
	504	! initialize the variables for microphysical physics
	505	!
	506	!
	507	do k = kts, kte
	508	do i = its, ite
	509	prevp(i,k) = 0.
	510	psdep(i,k) = 0.
	511	pgdep(i,k) = 0.
	512	praut(i,k) = 0.
	513	psaut(i,k) = 0.
	514	pgaut(i,k) = 0.
	515	pracw(i,k) = 0.
	516	praci(i,k) = 0.
	517	piacr(i,k) = 0.
	518	psaci(i,k) = 0.
	519	psacw(i,k) = 0.
	520	pracs(i,k) = 0.
	521	psacr(i,k) = 0.
	522	pgacw(i,k) = 0.
	523	paacw(i,k) = 0.
	524	pgaci(i,k) = 0.
	525	pgacr(i,k) = 0.
	526	pgacs(i,k) = 0.
	527	pigen(i,k) = 0.
	528	pidep(i,k) = 0.
	529	pcond(i,k) = 0.
	530	psmlt(i,k) = 0.
	531	pgmlt(i,k) = 0.
	532	pseml(i,k) = 0.
	533	pgeml(i,k) = 0.
	534	psevp(i,k) = 0.
	535	pgevp(i,k) = 0.
	536	pcact(i,k) = 0.
	537	falk(i,k,1) = 0.
	538	falk(i,k,2) = 0.
	539	falk(i,k,3) = 0.
	540	fall(i,k,1) = 0.
	541	fall(i,k,2) = 0.
	542	fall(i,k,3) = 0.
	543	fallc(i,k) = 0.
	544	falkc(i,k) = 0.
	545	falln(i,k) =0.
	546	falkn(i,k) =0.
	547	xni(i,k) = 1.e3
	548	nsacw(i,k) = 0.
	549	ngacw(i,k) = 0.
	550	naacw(i,k) = 0.
	551	niacr(i,k) = 0.
	552	nsacr(i,k) = 0.
	553	ngacr(i,k) = 0.
	554	nseml(i,k) = 0.
	555	ngeml(i,k) = 0.
	556	nracw(i,k) = 0.
	557	nccol(i,k) = 0.
	558	nrcol(i,k) = 0.
	559	ncact(i,k) = 0.
	560	nraut(i,k) = 0.
	561	ncevp(i,k) = 0.
	562	enddo
	563	enddo
	564	do k = kts, kte
	565	do i = its, ite
	566	if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin ) then
	567	rslopec(i,k) = rslopecmax
	568	rslopec2(i,k) = rslopec2max
	569	rslopec3(i,k) = rslopec3max
	570	else
	571	rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
	572	rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
	573	rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
	574	endif
	575	!-------------------------------------------------------------
	576	! Ni: ice crystal number concentraiton [HDC 5c]
	577	!-------------------------------------------------------------
	578	temp = (den(i,k)*max(qci(i,k,2),qmin))
	579	temp = sqrt(sqrt(temptemptemp))
	580	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
	581	enddo
	582	enddo
	583	!----------------------------------------------------------------
	584	! compute the fallout term:
	585	! first, vertical terminal velosity for minor loops
	586	!----------------------------------------------------------------
	587	do k = kts, kte
	588	do i = its, ite
	589	qrs_tmp(i,k,1) = qrs(i,k,1)
	590	qrs_tmp(i,k,2) = qrs(i,k,2)
	591	qrs_tmp(i,k,3) = qrs(i,k,3)
	592	ncr_tmp(i,k) = ncr(i,k,3)
	593	enddo
	594	enddo
	595	call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
	596	rslope3,work1,workn,its,ite,kts,kte)
	597	!
	598	! vt update for qr and nr
	599	mstepmax = 1
	600	numdt = 1
	601	do k = kte, kts, -1
	602	do i = its, ite
	603	work1(i,k,1) = work1(i,k,1)/delz(i,k)
	604	workn(i,k) = workn(i,k)/delz(i,k)
	605	numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
	606	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
	607	enddo
	608	enddo
	609	do i = its, ite
	610	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
	611	enddo
	612	!
	613	do n = 1, mstepmax
	614	k = kte
	615	do i = its, ite
	616	if(n.le.mstep(i)) then
	617	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
	618	falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
	619	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
	620	falln(i,k) = falln(i,k)+falkn(i,k)
	621	qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
	622	ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
	623	endif
	624	enddo
	625	do k = kte-1, kts, -1
	626	do i = its, ite
	627	if(n.le.mstep(i)) then
	628	falk(i,k,1) = den(i,k)qrs(i,k,1)work1(i,k,1)/mstep(i)
	629	falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
	630	fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
	631	falln(i,k) = falln(i,k)+falkn(i,k)
	632	qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
	633	delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	634	ncr(i,k,3) = max(ncr(i,k,3)-(falkn(i,k)-falkn(i,k+1)*delz(i,k+1) &
	635	/delz(i,k))*dtcld,0.)
	636	endif
	637	enddo
	638	enddo
	639	do k = kts, kte
	640	do i = its, ite
	641	qrs_tmp(i,k,1) = qrs(i,k,1)
	642	ncr_tmp(i,k) = ncr(i,k,3)
	643	enddo
	644	enddo
	645	call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
	646	rslope3,work1,workn,its,ite,kts,kte)
	647	do k = kte, kts, -1
	648	do i = its, ite
	649	work1(i,k,1) = work1(i,k,1)/delz(i,k)
	650	workn(i,k) = workn(i,k)/delz(i,k)
	651	enddo
	652	enddo
	653	enddo
	654	! for semi
	655	do k = kte, kts, -1
	656	do i = its, ite
	657	qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
	658	if(qsum(i,k) .gt. 1.e-15 ) then
	659	worka(i,k) = (work1(i,k,2)qrs(i,k,2) + work1(i,k,3)qrs(i,k,3)) &
	660	/qsum(i,k)
	661	else
	662	worka(i,k) = 0.
	663	endif
	664	denqrs2(i,k) = den(i,k)*qrs(i,k,2)
	665	denqrs3(i,k) = den(i,k)*qrs(i,k,3)
	666	enddo
	667	enddo
	668	call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
	669	denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
	670	do k = kts, kte
	671	do i = its, ite
	672	qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
	673	qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
	674	fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
	675	fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
	676	enddo
	677	enddo
	678	do i = its, ite
	679	fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
	680	fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
	681	enddo
	682	do k = kts, kte
	683	do i = its, ite
	684	qrs_tmp(i,k,1) = qrs(i,k,1)
	685	qrs_tmp(i,k,2) = qrs(i,k,2)
	686	qrs_tmp(i,k,3) = qrs(i,k,3)
	687	ncr_tmp(i,k) = ncr(i,k,3)
	688	enddo
	689	enddo
	690	call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
	691	rslope3,work1,workn,its,ite,kts,kte)
	692	!
	693	do k = kte, kts, -1
	694	do i = its, ite
	695	supcol = t0c-t(i,k)
	696	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
	697	if(t(i,k).gt.t0c) then
	698	!---------------------------------------------------------------
	699	! psmlt: melting of snow [HL A33] [RH83 A25]
	700	! (T>T0: QS->QR)
	701	!---------------------------------------------------------------
	702	xlf = xlf0
	703	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
	704	if(qrs(i,k,2).gt.0.) then
	705	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
	706	psmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))pi/2. &
	707	n0sfac(i,k)(precs1*rslope2(i,k,2) &
	708	+precs2work2(i,k)coeres)
	709	psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2) &
	710	/mstep(i)),0.)
	711	qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
	712	qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
	713	!-------------------------------------------------------------------
	714	! nsmlt: melting of snow [LH A27]
	715	! (T>T0: ->NR)
	716	!-------------------------------------------------------------------
	717	if(qrs(i,k,2).gt.qcrmin) then
	718	sfac = rslope(i,k,2)n0sn0sfac(i,k)/qrs(i,k,2)
	719	ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
	720	endif
	721	t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
	722	endif
	723	!---------------------------------------------------------------
	724	! pgmlt: melting of graupel [HL A23] [LFO 47]
	725	! (T>T0: QG->QR)
	726	!---------------------------------------------------------------
	727	if(qrs(i,k,3).gt.0.) then
	728	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
	729	pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf(t0c-t(i,k))(precg1 &
	730	rslope2(i,k,3) + precg2work2(i,k)*coeres)
	731	pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
	732	-qrs(i,k,3)/mstep(i)),0.)
	733	qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
	734	qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
	735	!-------------------------------------------------------------------
	736	! ngmlt: melting of graupel [LH A28]
	737	! (T>T0: ->NR)
	738	!-------------------------------------------------------------------
	739	if(qrs(i,k,3).gt.qcrmin) then
	740	gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
	741	ncr(i,k,3) = ncr(i,k,3) - gfac*pgmlt(i,k)
	742	endif
	743	t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
	744	endif
	745	endif
	746	enddo
	747	enddo
	748	!---------------------------------------------------------------
	749	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
	750	!---------------------------------------------------------------
	751	do k = kte, kts, -1
	752	do i = its, ite
	753	if(qci(i,k,2).le.0.) then
	754	work1c(i,k) = 0.
	755	else
	756	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
	757	diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
	758	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
	759	endif
	760	enddo
	761	enddo
	762	!
	763	! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
	764	!
	765	do k = kte, kts, -1
	766	do i = its, ite
	767	denqci(i,k) = den(i,k)*qci(i,k,2)
	768	enddo
	769	enddo
	770	call nislfv_rain_plmr(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,denqci, &
	771	delqi,dtcld,1,0,0)
	772	do k = kts, kte
	773	do i = its, ite
	774	qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
	775	enddo
	776	enddo
	777	do i = its, ite
	778	fallc(i,1) = delqi(i)/delz(i,1)/dtcld
	779	enddo
	780	!
	781	!----------------------------------------------------------------
	782	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
	783	!
	784	do i = its, ite
	785	fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
	786	fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
	787	fallsum_qg = fall(i,kts,3)
	788	if(fallsum.gt.0.) then
	789	rainncv(i) = fallsumdelz(i,kts)/denrdtcld*1000. + rainncv(i)
	790	rain(i) = fallsumdelz(i,kts)/denrdtcld*1000. + rain(i)
	791	endif
	792	if(fallsum_qsi.gt.0.) then
	793	tstepsnow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + tstepsnow(i)
	794	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
	795	snowncv(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snowncv(i)
	796	snow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i)
	797	ENDIF
	798	endif
	799	if(fallsum_qg.gt.0.) then
	800	tstepgraup(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000. &
	801	+ tstepgraup(i)
	802	IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
	803	graupelncv(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000. &
	804	+ graupelncv(i)
	805	graupel(i) = fallsum_qgdelz(i,kts)/denrdtcld*1000. + graupel(i)
	806	ENDIF
	807	endif
	808	! if(fallsum.gt.0.) sr(i) = (snowncv(i) + graupelncv(i)) &
	809	if(fallsum.gt.0.) sr(i) = (tstepsnow(i) + tstepgraup(i)) &
	810	/(rainncv(i)+1.e-12)
	811	enddo
	812	!
	813	!---------------------------------------------------------------
	814	! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
	815	! (T>T0: QI->QC)
	816	!---------------------------------------------------------------
	817	do k = kts, kte
	818	do i = its, ite
	819	supcol = t0c-t(i,k)
	820	xlf = xls-xl(i,k)
	821	if(supcol.lt.0.) xlf = xlf0
	822	if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
	823	qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
	824	!---------------------------------------------------------------
	825	! nimlt: instantaneous melting of cloud ice [LH A18]
	826	! (T>T0: ->NC)
	827	!--------------------------------------------------------------
	828	ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
	829	t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
	830	qci(i,k,2) = 0.
	831	endif
	832	!---------------------------------------------------------------
	833	! pihmf: homogeneous of cloud water below -40c [HL A45]
	834	! (T<-40C: QC->QI)
	835	!---------------------------------------------------------------
	836	if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
	837	qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
	838	!---------------------------------------------------------------
	839	! nihmf: homogeneous of cloud water below -40c [LH A17]
	840	! (T<-40C: NC->)
	841	!---------------------------------------------------------------
	842	if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
	843	t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
	844	qci(i,k,1) = 0.
	845	endif
	846	!---------------------------------------------------------------
	847	! pihtf: heterogeneous of cloud water [HL A44]
	848	! (T0>T>-40C: QC->QI)
	849	!---------------------------------------------------------------
	850	if(supcol.gt.0. .and. qci(i,k,1).gt.qmin) then
	851	supcolt=min(supcol,70.)
	852	pfrzdtc = min(pipipfrz1(exp(pfrz2supcolt)-1.)*denr/den(i,k) &
	853	ncr(i,k,2)rslopec3(i,k)rslopec3(i,k)/18.dtcld &
	854	,qci(i,k,1))
	855	!---------------------------------------------------------------
	856	! nihtf: heterogeneous of cloud water [LH A16]
	857	! (T0>T>-40C: NC->)
	858	!---------------------------------------------------------------
	859	if(ncr(i,k,2).gt.ncmin) then
	860	nfrzdtc = min(pipfrz1(exp(pfrz2supcolt)-1.)ncr(i,k,2) &
	861	rslopec3(i,k)/6.dtcld,ncr(i,k,2))
	862	ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
	863	endif
	864	qci(i,k,2) = qci(i,k,2) + pfrzdtc
	865	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
	866	qci(i,k,1) = qci(i,k,1)-pfrzdtc
	867	endif
	868	!---------------------------------------------------------------
	869	! pgfrz: freezing of rain water [HL A20] [LFO 45]
	870	! (T<T0, QR->QG)
	871	!---------------------------------------------------------------
	872	if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
	873	supcolt=min(supcol,70.)
	874	pfrzdtr = min(140.(pipi)pfrz1ncr(i,k,3)*denr/den(i,k) &
	875	(exp(pfrz2supcolt)-1.)rslope3(i,k,1)rslope3(i,k,1) &
	876	*dtcld,qrs(i,k,1))
	877	!---------------------------------------------------------------
	878	! ngfrz: freezing of rain water [LH A26]
	879	! (T<T0, NR-> )
	880	!---------------------------------------------------------------
	881	if(ncr(i,k,3).gt.nrmin) then
	882	nfrzdtr = min(4.pipfrz1ncr(i,k,3)(exp(pfrz2*supcolt)-1.) &
	883	rslope3(i,k,1)dtcld, ncr(i,k,3))
	884	ncr(i,k,3) = ncr(i,k,3) - nfrzdtr
	885	endif
	886	qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
	887	t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
	888	qrs(i,k,1) = qrs(i,k,1) - pfrzdtr
	889	endif
	890	enddo
	891	enddo
	892	!
	893	do k = kts, kte
	894	do i = its, ite
	895	ncr(i,k,2) = max(ncr(i,k,2),0.0)
	896	ncr(i,k,3) = max(ncr(i,k,3),0.0)
	897	enddo
	898	enddo
	899	!
	900	!----------------------------------------------------------------
	901	! update the slope parameters for microphysics computation
	902	!
	903	do k = kts, kte
	904	do i = its, ite
	905	qrs_tmp(i,k,1) = qrs(i,k,1)
	906	qrs_tmp(i,k,2) = qrs(i,k,2)
	907	qrs_tmp(i,k,3) = qrs(i,k,3)
	908	ncr_tmp(i,k) = ncr(i,k,3)
	909	enddo
	910	enddo
	911	call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
	912	rslope3,work1,workn,its,ite,kts,kte)
	913	do k = kts, kte
	914	do i = its, ite
	915	!-----------------------------------------------------------------
	916	! compute the mean-volume drop diameter [LH A10]
	917	! for raindrop distribution
	918	!-----------------------------------------------------------------
	919	avedia(i,k,2) = rslope(i,k,1)((24.)*(.3333333))
	920	!
	921	if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
	922	rslopec(i,k) = rslopecmax
	923	rslopec2(i,k) = rslopec2max
	924	rslopec3(i,k) = rslopec3max
	925	else
	926	rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
	927	rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
	928	rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
	929	endif
	930	!--------------------------------------------------------------------
	931	! compute the mean-volume drop diameter [LH A7]
	932	! for cloud-droplet distribution
	933	!--------------------------------------------------------------------
	934	avedia(i,k,1) = rslopec(i,k)
	935	enddo
	936	enddo
	937	!
	938	do k = kts, kte
	939	do i = its, ite
	940	work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
	941	work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
	942	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
	943	enddo
	944	enddo
	945	!
	946	!===============================================================
	947	!
	948	! warm rain processes
	949	!
	950	! - follows the double-moment processes in Lim and Hong
	951	!
	952	!===============================================================
	953	!
	954	do k = kts, kte
	955	do i = its, ite
	956	supsat = max(q(i,k),qmin)-qs(i,k,1)
	957	satdt = supsat/dtcld
	958	!---------------------------------------------------------------
	959	! praut: auto conversion rate from cloud to rain [LH 9] [CP 17]
	960	! (QC->QR)
	961	!--------------------------------------------------------------
	962	lencon = 2.7e-2den(i,k)qci(i,k,1)(1.e20/16.rslopec2(i,k) &
	963	*rslopec2(i,k)-0.4)
	964	lenconcr = max(1.2*lencon, qcrmin)
	965	if(avedia(i,k,1).gt.di15) then
	966	taucon = 3.7/den(i,k)/qci(i,k,1)/(0.5e6*rslopec(i,k)-7.5)
	967	praut(i,k) = lencon/taucon
	968	praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
	969	!---------------------------------------------------------------
	970	! nraut: auto conversion rate from cloud to rain [LH A6] [CP 18 & 19]
	971	! (NC->NR)
	972	!---------------------------------------------------------------
	973	nraut(i,k) = 3.5e9den(i,k)praut(i,k)
	974	if(qrs(i,k,1).gt.lenconcr) &
	975	nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
	976	nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
	977	endif
	978	!---------------------------------------------------------------
	979	! pracw: accretion of cloud water by rain [LH 10] [CP 22 & 23]
	980	! (QC->QR)
	981	! nracw: accretion of cloud water by rain [LH A9]
	982	! (NC->)
	983	!---------------------------------------------------------------
	984	if(qrs(i,k,1).ge.lenconcr) then
	985	if(avedia(i,k,2).ge.di100) then
	986	nracw(i,k) = min(ncrk1ncr(i,k,2)ncr(i,k,3)*(rslopec3(i,k) &
	987	+ 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
	988	pracw(i,k) = min(pi/6.(denr/den(i,k))ncrk1*ncr(i,k,2) &
	989	ncr(i,k,3)rslopec3(i,k)(2.rslopec3(i,k) &
	990	+ 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
	991	else
	992	nracw(i,k) = min(ncrk2ncr(i,k,2)ncr(i,k,3)(2.rslopec3(i,k) &
	993	rslopec3(i,k)+5040.rslope3(i,k,1) &
	994	*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
	995	pracw(i,k) = min(pi/6.(denr/den(i,k))ncrk2*ncr(i,k,2) &
	996	ncr(i,k,3)rslopec3(i,k)(6.rslopec3(i,k) &
	997	rslopec3(i,k)+5040.rslope3(i,k,1)*rslope3(i,k,1)) &
	998	,qci(i,k,1)/dtcld)
	999	endif
	1000	endif
	1001	!----------------------------------------------------------------
	1002	! nccol: self collection of cloud water [LH A8] [CP 24 & 25]
	1003	! (NC->)
	1004	!----------------------------------------------------------------
	1005	if(avedia(i,k,1).ge.di100) then
	1006	nccol(i,k) = ncrk1ncr(i,k,2)ncr(i,k,2)*rslopec3(i,k)
	1007	else
	1008	nccol(i,k) = 2.ncrk2ncr(i,k,2)ncr(i,k,2)rslopec3(i,k) &
	1009	*rslopec3(i,k)
	1010	endif
	1011	!----------------------------------------------------------------
	1012	! nrcol: self collection of rain-drops and break-up [LH A21] [CP 24 & 25]
	1013	! (NR->)
	1014	!----------------------------------------------------------------
	1015	if(qrs(i,k,1).ge.lenconcr) then
	1016	if(avedia(i,k,2).lt.di100) then
	1017	nrcol(i,k) = 5040.ncrk2ncr(i,k,3)ncr(i,k,3)rslope3(i,k,1) &
	1018	*rslope3(i,k,1)
	1019	elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
	1020	nrcol(i,k) = 24.ncrk1ncr(i,k,3)ncr(i,k,3)rslope3(i,k,1)
	1021	elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
	1022	coecol = -2.5e3*(avedia(i,k,2)-di600)
	1023	nrcol(i,k) = 24.exp(coecol)ncrk1ncr(i,k,3)ncr(i,k,3) &
	1024	*rslope3(i,k,1)
	1025	else
	1026	nrcol(i,k) = 0.
	1027	endif
	1028	endif
	1029	!---------------------------------------------------------------
	1030	! prevp: evaporation/condensation rate of rain [HL A41]
	1031	! (QV->QR or QR->QV)
	1032	!---------------------------------------------------------------
	1033	if(qrs(i,k,1).gt.0.) then
	1034	coeres = rslope(i,k,1)sqrt(rslope(i,k,1)rslopeb(i,k,1))
	1035	prevp(i,k) = (rh(i,k,1)-1.)ncr(i,k,3)(precr1*rslope(i,k,1) &
	1036	+ precr2work2(i,k)coeres)/work1(i,k,1)
	1037	if(prevp(i,k).lt.0.) then
	1038	prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
	1039	prevp(i,k) = max(prevp(i,k),satdt/2)
	1040	!----------------------------------------------------------------
	1041	! Nrevp: evaporation/condensation rate of rain [LH A14]
	1042	! (NR->NCCN)
	1043	!----------------------------------------------------------------
	1044	if(prevp(i,k).eq.-qrs(i,k,1)/dtcld) then
	1045	ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,3)
	1046	ncr(i,k,3) = 0.
	1047	endif
	1048	else
	1049	!
	1050	prevp(i,k) = min(prevp(i,k),satdt/2)
	1051	endif
	1052	endif
	1053	enddo
	1054	enddo
	1055	!
	1056	!===============================================================
	1057	!
	1058	! cold rain processes
	1059	!
	1060	! - follows the revised ice microphysics processes in HDC
	1061	! - the processes same as in RH83 and RH84 and LFO behave
	1062	! following ice crystal hapits defined in HDC, inclduing
	1063	! intercept parameter for snow (n0s), ice crystal number
	1064	! concentration (ni), ice nuclei number concentration
	1065	! (n0i), ice diameter (d)
	1066	!
	1067	!===============================================================
	1068	!
	1069	do k = kts, kte
	1070	do i = its, ite
	1071	supcol = t0c-t(i,k)
	1072	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
	1073	supsat = max(q(i,k),qmin)-qs(i,k,2)
	1074	satdt = supsat/dtcld
	1075	ifsat = 0
	1076	!-------------------------------------------------------------
	1077	! Ni: ice crystal number concentraiton [HDC 5c]
	1078	!-------------------------------------------------------------
	1079	! xni(i,k) = min(max(5.38e7*(den(i,k) &
	1080	! max(qci(i,k,2),qmin))*0.75,1.e3),1.e6)
	1081	temp = (den(i,k)*max(qci(i,k,2),qmin))
	1082	temp = sqrt(sqrt(temptemptemp))
	1083	xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
	1084	eacrs = exp(0.07*(-supcol))
	1085	!
	1086	xmi = den(i,k)*qci(i,k,2)/xni(i,k)
	1087	diameter = min(dicon * sqrt(xmi),dimax)
	1088	vt2i = 1.49e4diameter*1.31
	1089	vt2r=pvtrrslopeb(i,k,1)denfac(i,k)
	1090	vt2s=pvtsrslopeb(i,k,2)denfac(i,k)
	1091	vt2g=pvtgrslopeb(i,k,3)denfac(i,k)
	1092	qsum(i,k) = max((qrs(i,k,2)+qrs(i,k,3)),1.e-15)
	1093	if(qsum(i,k) .gt. 1.e-15) then
	1094	vt2ave=(vt2sqrs(i,k,2)+vt2gqrs(i,k,3))/(qsum(i,k))
	1095	else
	1096	vt2ave=0.
	1097	endif
	1098	if(supcol.gt.0. .and. qci(i,k,2).gt.qmin) then
	1099	if(qrs(i,k,1).gt.qcrmin) then
	1100	!-------------------------------------------------------------
	1101	! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
	1102	! (T<T0: QI->QR)
	1103	!-------------------------------------------------------------
	1104	acrfac = 6.rslope2(i,k,1)+4.diameterrslope(i,k,1) + diameter*2
	1105	praci(i,k) = piqci(i,k,2)ncr(i,k,3)abs(vt2r-vt2i)acrfac/4.
	1106	praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
	1107	!-------------------------------------------------------------
	1108	! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
	1109	! (T<T0: QR->QS or QR->QG)
	1110	!-------------------------------------------------------------
	1111	piacr(i,k) = pipiavtrncr(i,k,3)denrxni(i,k)denfac(i,k) &
	1112	g7pbrrslope3(i,k,1)rslope2(i,k,1)rslopeb(i,k,1) &
	1113	/24./den(i,k)
	1114	piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
	1115	endif
	1116	!-------------------------------------------------------------
	1117	! niacr: Accretion of rain by cloud ice [LH A25]
	1118	! (T<T0: NR->)
	1119	!-------------------------------------------------------------
	1120	if(ncr(i,k,3).gt.nrmin) then
	1121	niacr(i,k) = piavtrncr(i,k,3)xni(i,k)denfac(i,k)*g4pbr &
	1122	rslope2(i,k,1)rslopeb(i,k,1)/4.
	1123	niacr(i,k) = min(niacr(i,k),ncr(i,k,3)/dtcld)
	1124	endif
	1125	!-------------------------------------------------------------
	1126	! psaci: Accretion of cloud ice by snow [HDC 10]
	1127	! (T<T0: QI->QS)
	1128	!-------------------------------------------------------------
	1129	if(qrs(i,k,2).gt.qcrmin) then
	1130	acrfac = 2.rslope3(i,k,2)+2.diameter*rslope2(i,k,2) &
	1131	+ diameter*2rslope(i,k,2)
	1132	psaci(i,k) = piqci(i,k,2)eacrsn0sn0sfac(i,k) &
	1133	abs(vt2ave-vt2i)acrfac/4.
	1134	psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
	1135	endif
	1136	!-------------------------------------------------------------
	1137	! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
	1138	! (T<T0: QI->QG)
	1139	!-------------------------------------------------------------
	1140	if(qrs(i,k,3).gt.qcrmin) then
	1141	egi = exp(0.07*(-supcol))
	1142	acrfac = 2.rslope3(i,k,3)+2.diameter*rslope2(i,k,3) &
	1143	+ diameter*2rslope(i,k,3)
	1144	pgaci(i,k) = piegiqci(i,k,2)n0gabs(vt2ave-vt2i)*acrfac/4.
	1145	pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
	1146	endif
	1147	endif
	1148	!-------------------------------------------------------------
	1149	! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
	1150	! (T<T0: QC->QS, and T>=T0: QC->QR)
	1151	!-------------------------------------------------------------
	1152	if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
	1153	psacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2)*rslopeb(i,k,2) &
	1154	qci(i,k,1)denfac(i,k),qci(i,k,1)/dtcld)
	1155	endif
	1156	!-------------------------------------------------------------
	1157	! nsacw: Accretion of cloud water by snow [LH A12]
	1158	! (NC ->)
	1159	!-------------------------------------------------------------
	1160	if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
	1161	nsacw(i,k) = min(pacrcn0sfac(i,k)rslope3(i,k,2)*rslopeb(i,k,2) &
	1162	ncr(i,k,2)denfac(i,k),ncr(i,k,2)/dtcld)
	1163	endif
	1164	!-------------------------------------------------------------
	1165	! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
	1166	! (T<T0: QC->QG, and T>=T0: QC->QR)
	1167	!-------------------------------------------------------------
	1168	if(qrs(i,k,3).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
	1169	pgacw(i,k) = min(pacrgrslope3(i,k,3)rslopeb(i,k,3)*qci(i,k,1) &
	1170	*denfac(i,k),qci(i,k,1)/dtcld)
	1171	endif
	1172	!-------------------------------------------------------------
	1173	! ngacw: Accretion of cloud water by graupel [LH A13]
	1174	! (NC->
	1175	!-------------------------------------------------------------
	1176	if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
	1177	ngacw(i,k) = min(pacrgrslope3(i,k,3)rslopeb(i,k,3)*ncr(i,k,2) &
	1178	*denfac(i,k),ncr(i,k,2)/dtcld)
	1179	endif
	1180	!-------------------------------------------------------------
	1181	! paacw: Accretion of cloud water by averaged snow/graupel
	1182	! (T<T0: QC->QG or QS, and T>=T0: QC->QR)
	1183	!-------------------------------------------------------------
	1184	if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,3).gt.qcrmin) then
	1185	paacw(i,k) = (qrs(i,k,2)psacw(i,k)+qrs(i,k,3)pgacw(i,k))/(qsum(i,k))
	1186	!-------------------------------------------------------------
	1187	! naacw: Accretion of cloud water by averaged snow/graupel
	1188	! (Nc->)
	1189	!-------------------------------------------------------------
	1190	naacw(i,k) = (qrs(i,k,2)nsacw(i,k)+qrs(i,k,3)ngacw(i,k))/(qsum(i,k))
	1191	endif
	1192	!-------------------------------------------------------------
	1193	! pracs: Accretion of snow by rain [HL A11] [LFO 27]
	1194	! (T<T0: QS->QG)
	1195	!-------------------------------------------------------------
	1196	if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
	1197	if(supcol.gt.0) then
	1198	acrfac = 5.rslope3(i,k,2)rslope3(i,k,2) &
	1199	+ 4.rslope3(i,k,2)rslope2(i,k,2)*rslope(i,k,1) &
	1200	+ 1.5rslope2(i,k,2)rslope2(i,k,2)*rslope2(i,k,1)
	1201	pracs(i,k) = pipincr(i,k,3)n0sn0sfac(i,k)*abs(vt2r-vt2ave) &
	1202	(dens/den(i,k))acrfac
	1203	pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
	1204	endif
	1205	!-------------------------------------------------------------
	1206	! psacr: Accretion of rain by snow [HL A10] [LFO 28]
	1207	! (T<T0:QR->QS or QR->QG) (T>=T0: enhance melting of snow)
	1208	!-------------------------------------------------------------
	1209	acrfac = 30.rslope3(i,k,1)rslope2(i,k,1)*rslope(i,k,2) &
	1210	+ 5.rslope2(i,k,1)rslope2(i,k,1)*rslope2(i,k,2) &
	1211	+ 2.rslope3(i,k,1)rslope3(i,k,2)
	1212	psacr(i,k) = pipincr(i,k,3)n0sn0sfac(i,k)*abs(vt2ave-vt2r) &
	1213	(denr/den(i,k))acrfac
	1214	psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
	1215	endif
	1216	if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
	1217	!-------------------------------------------------------------
	1218	! nsacr: Accretion of rain by snow [LH A23]
	1219	! (T<T0:NR->)
	1220	!-------------------------------------------------------------
	1221	acrfac = 1.5rslope2(i,k,1)rslope(i,k,2) &
	1222	+ 1.0rslope(i,k,1)rslope2(i,k,2)+.5*rslope3(i,k,2)
	1223	nsacr(i,k) = pincr(i,k,3)n0sn0sfac(i,k)abs(vt2ave-vt2r) &
	1224	*acrfac
	1225	nsacr(i,k) = min(nsacr(i,k),ncr(i,k,3)/dtcld)
	1226	endif
	1227	!-------------------------------------------------------------
	1228	! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
	1229	! (T<T0: QR->QG) (T>=T0: enhance melting of graupel)
	1230	!-------------------------------------------------------------
	1231	if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
	1232	acrfac = 30.rslope3(i,k,1)rslope2(i,k,1)*rslope(i,k,3) &
	1233	+ 5.rslope2(i,k,1)rslope2(i,k,1)*rslope2(i,k,3) &
	1234	+ 2.rslope3(i,k,1)rslope3(i,k,3)
	1235	pgacr(i,k) = pipincr(i,k,3)n0gabs(vt2ave-vt2r)*(denr/den(i,k)) &
	1236	*acrfac
	1237	pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
	1238	endif
	1239	!-------------------------------------------------------------
	1240	! ngacr: Accretion of rain by graupel [LH A24]
	1241	! (T<T0: NR->)
	1242	!-------------------------------------------------------------
	1243	if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
	1244	acrfac = 1.5rslope2(i,k,1)rslope(i,k,3) &
	1245	+ 1.0rslope(i,k,1)rslope2(i,k,3) + .5*rslope3(i,k,3)
	1246	ngacr(i,k) = pincr(i,k,3)n0gabs(vt2ave-vt2r)acrfac
	1247	ngacr(i,k) = min(ngacr(i,k),ncr(i,k,3)/dtcld)
	1248	endif
	1249	!
	1250	!-------------------------------------------------------------
	1251	! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
	1252	! (QS->QG) : This process is eliminated in V3.0 with the
	1253	! new combined snow/graupel fall speeds
	1254	!-------------------------------------------------------------
	1255	if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,2).gt.qcrmin) then
	1256	pgacs(i,k) = 0.
	1257	endif
	1258	if(supcol.le.0) then
	1259	xlf = xlf0
	1260	!-------------------------------------------------------------
	1261	! pseml: Enhanced melting of snow by accretion of water [HL A34]
	1262	! (T>=T0: QS->QR)
	1263	!-------------------------------------------------------------
	1264	if(qrs(i,k,2).gt.0.) &
	1265	pseml(i,k) = min(max(cliqsupcol(paacw(i,k)+psacr(i,k)) &
	1266	/xlf,-qrs(i,k,2)/dtcld),0.)
	1267	!--------------------------------------------------------------
	1268	! nseml: Enhanced melting of snow by accretion of water [LH A29]
	1269	! (T>=T0: ->NR)
	1270	!--------------------------------------------------------------
	1271	if (qrs(i,k,2).gt.qcrmin) then
	1272	sfac = rslope(i,k,2)n0sn0sfac(i,k)/qrs(i,k,2)
	1273	nseml(i,k) = -sfac*pseml(i,k)
	1274	endif
	1275	!-------------------------------------------------------------
	1276	! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
	1277	! (T>=T0: QG->QR)
	1278	!-------------------------------------------------------------
	1279	if(qrs(i,k,3).gt.0.) &
	1280	pgeml(i,k) = min(max(cliqsupcol(paacw(i,k)+pgacr(i,k))/xlf &
	1281	,-qrs(i,k,3)/dtcld),0.)
	1282	!--------------------------------------------------------------
	1283	! ngeml: Enhanced melting of graupel by accretion of water [LH A30]
	1284	! (T>=T0: -> NR)
	1285	!--------------------------------------------------------------
	1286	if (qrs(i,k,3).gt.qcrmin) then
	1287	gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
	1288	ngeml(i,k) = -gfac*pgeml(i,k)
	1289	endif
	1290	endif
	1291	if(supcol.gt.0) then
	1292	!-------------------------------------------------------------
	1293	! pidep: Deposition/Sublimation rate of ice [HDC 9]
	1294	! (T<T0: QV->QI or QI->QV)
	1295	!-------------------------------------------------------------
	1296	if(qci(i,k,2).gt.0. .and. ifsat.ne.1) then
	1297	pidep(i,k) = 4.diameterxni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
	1298	supice = satdt-prevp(i,k)
	1299	if(pidep(i,k).lt.0.) then
	1300	pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
	1301	pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
	1302	else
	1303	pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
	1304	endif
	1305	if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
	1306	endif
	1307	!-------------------------------------------------------------
	1308	! psdep: deposition/sublimation rate of snow [HDC 14]
	1309	! (T<T0: QV->QS or QS->QV)
	1310	!-------------------------------------------------------------
	1311	if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then
	1312	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
	1313	psdep(i,k) = (rh(i,k,2)-1.)n0sfac(i,k)(precs1*rslope2(i,k,2) &
	1314	+ precs2work2(i,k)coeres)/work1(i,k,2)
	1315	supice = satdt-prevp(i,k)-pidep(i,k)
	1316	if(psdep(i,k).lt.0.) then
	1317	psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
	1318	psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
	1319	else
	1320	psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
	1321	endif
	1322	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
	1323	endif
	1324	!-------------------------------------------------------------
	1325	! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
	1326	! (T<T0: QV->QG or QG->QV)
	1327	!-------------------------------------------------------------
	1328	if(qrs(i,k,3).gt.0. .and. ifsat.ne.1) then
	1329	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
	1330	pgdep(i,k) = (rh(i,k,2)-1.)(precg1rslope2(i,k,3) &
	1331	+ precg2work2(i,k)coeres)/work1(i,k,2)
	1332	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
	1333	if(pgdep(i,k).lt.0.) then
	1334	pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
	1335	pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
	1336	else
	1337	pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
	1338	endif
	1339	if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
	1340	abs(satdt)) ifsat = 1
	1341	endif
	1342	!-------------------------------------------------------------
	1343	! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
	1344	! (T<T0: QV->QI)
	1345	!-------------------------------------------------------------
	1346	if(supsat.gt.0. .and. ifsat.ne.1) then
	1347	supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
	1348	xni0 = 1.e3exp(0.1supcol)
	1349	roqi0 = 4.92e-11xni0*1.33
	1350	pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
	1351	pigen(i,k) = min(min(pigen(i,k),satdt),supice)
	1352	endif
	1353	!
	1354	!-------------------------------------------------------------
	1355	! psaut: conversion(aggregation) of ice to snow [HDC 12]
	1356	! (T<T0: QI->QS)
	1357	!-------------------------------------------------------------
	1358	if(qci(i,k,2).gt.0.) then
	1359	qimax = roqimax/den(i,k)
	1360	psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
	1361	endif
	1362	!
	1363	!-------------------------------------------------------------
	1364	! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
	1365	! (T<T0: QS->QG)
	1366	!-------------------------------------------------------------
	1367	if(qrs(i,k,2).gt.0.) then
	1368	alpha2 = 1.e-3exp(0.09(-supcol))
	1369	pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
	1370	endif
	1371	endif
	1372	!
	1373	!-------------------------------------------------------------
	1374	! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
	1375	! (T>=T0: QS->QV)
	1376	!-------------------------------------------------------------
	1377	if(supcol.lt.0.) then
	1378	if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.) then
	1379	coeres = rslope2(i,k,2)sqrt(rslope(i,k,2)rslopeb(i,k,2))
	1380	psevp(i,k) = (rh(i,k,1)-1.)n0sfac(i,k)(precs1*rslope2(i,k,2) &
	1381	+precs2work2(i,k)coeres)/work1(i,k,1)
	1382	psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
	1383	endif
	1384	!-------------------------------------------------------------
	1385	! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
	1386	! (T>=T0: QG->QV)
	1387	!-------------------------------------------------------------
	1388	if(qrs(i,k,3).gt.0. .and. rh(i,k,1).lt.1.) then
	1389	coeres = rslope2(i,k,3)sqrt(rslope(i,k,3)rslopeb(i,k,3))
	1390	pgevp(i,k) = (rh(i,k,1)-1.)(precg1rslope2(i,k,3) &
	1391	+ precg2work2(i,k)coeres)/work1(i,k,1)
	1392	pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
	1393	endif
	1394	endif
	1395	enddo
	1396	enddo
	1397	!
	1398	!
	1399	!----------------------------------------------------------------
	1400	! check mass conservation of generation terms and feedback to the
	1401	! large scale
	1402	!
	1403	do k = kts, kte
	1404	do i = its, ite
	1405	!
	1406	delta2=0.
	1407	delta3=0.
	1408	if(qrs(i,k,1).lt.1.e-4 .and. qrs(i,k,2).lt.1.e-4) delta2=1.
	1409	if(qrs(i,k,1).lt.1.e-4) delta3=1.
	1410	if(t(i,k).le.t0c) then
	1411	!
	1412	! cloud water
	1413	!
	1414	value = max(qmin,qci(i,k,1))
	1415	source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k))&
	1416	*dtcld
	1417	if (source.gt.value) then
	1418	factor = value/source
	1419	praut(i,k) = praut(i,k)*factor
	1420	pracw(i,k) = pracw(i,k)*factor
	1421	paacw(i,k) = paacw(i,k)*factor
	1422	endif
	1423	!
	1424	! cloud ice
	1425	!
	1426	value = max(qmin,qci(i,k,2))
	1427	source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
	1428	+pgaci(i,k))*dtcld
	1429	if (source.gt.value) then
	1430	factor = value/source
	1431	psaut(i,k) = psaut(i,k)*factor
	1432	pigen(i,k) = pigen(i,k)*factor
	1433	pidep(i,k) = pidep(i,k)*factor
	1434	praci(i,k) = praci(i,k)*factor
	1435	psaci(i,k) = psaci(i,k)*factor
	1436	pgaci(i,k) = pgaci(i,k)*factor
	1437	endif
	1438	!
	1439	! rain
	1440	!
	1441	value = max(qmin,qrs(i,k,1))
	1442	source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) &
	1443	+psacr(i,k)+pgacr(i,k))*dtcld
	1444	if (source.gt.value) then
	1445	factor = value/source
	1446	praut(i,k) = praut(i,k)*factor
	1447	prevp(i,k) = prevp(i,k)*factor
	1448	pracw(i,k) = pracw(i,k)*factor
	1449	piacr(i,k) = piacr(i,k)*factor
	1450	psacr(i,k) = psacr(i,k)*factor
	1451	pgacr(i,k) = pgacr(i,k)*factor
	1452	endif
	1453	!
	1454	! snow
	1455	!
	1456	value = max(qmin,qrs(i,k,2))
	1457	source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) &
	1458	+piacr(i,k)delta3+praci(i,k)delta3 &
	1459	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2 &
	1460	+psaci(i,k)-pgacs(i,k) )*dtcld
	1461	if (source.gt.value) then
	1462	factor = value/source
	1463	psdep(i,k) = psdep(i,k)*factor
	1464	psaut(i,k) = psaut(i,k)*factor
	1465	pgaut(i,k) = pgaut(i,k)*factor
	1466	paacw(i,k) = paacw(i,k)*factor
	1467	piacr(i,k) = piacr(i,k)*factor
	1468	praci(i,k) = praci(i,k)*factor
	1469	psaci(i,k) = psaci(i,k)*factor
	1470	pracs(i,k) = pracs(i,k)*factor
	1471	psacr(i,k) = psacr(i,k)*factor
	1472	pgacs(i,k) = pgacs(i,k)*factor
	1473	endif
	1474	!
	1475	! graupel
	1476	!
	1477	value = max(qmin,qrs(i,k,3))
	1478	source = -(pgdep(i,k)+pgaut(i,k) &
	1479	+piacr(i,k)(1.-delta3)+praci(i,k)(1.-delta3) &
	1480	+psacr(i,k)(1.-delta2)+pracs(i,k)(1.-delta2) &
	1481	+pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
	1482	if (source.gt.value) then
	1483	factor = value/source
	1484	pgdep(i,k) = pgdep(i,k)*factor
	1485	pgaut(i,k) = pgaut(i,k)*factor
	1486	piacr(i,k) = piacr(i,k)*factor
	1487	praci(i,k) = praci(i,k)*factor
	1488	psacr(i,k) = psacr(i,k)*factor
	1489	pracs(i,k) = pracs(i,k)*factor
	1490	paacw(i,k) = paacw(i,k)*factor
	1491	pgaci(i,k) = pgaci(i,k)*factor
	1492	pgacr(i,k) = pgacr(i,k)*factor
	1493	pgacs(i,k) = pgacs(i,k)*factor
	1494	endif
	1495	!
	1496	! cloud
	1497	!
	1498	value = max(ncmin,ncr(i,k,2))
	1499	source = (nraut(i,k)+nccol(i,k)+nracw(i,k) &
	1500	+naacw(i,k)+naacw(i,k))*dtcld
	1501	if (source.gt.value) then
	1502	factor = value/source
	1503	nraut(i,k) = nraut(i,k)*factor
	1504	nccol(i,k) = nccol(i,k)*factor
	1505	nracw(i,k) = nracw(i,k)*factor
	1506	naacw(i,k) = naacw(i,k)*factor
	1507	endif
	1508	!
	1509	! rain
	1510	!
	1511	value = max(nrmin,ncr(i,k,3))
	1512	source = (-nraut(i,k)+nrcol(i,k)+niacr(i,k)+nsacr(i,k)+ngacr(i,k) &
	1513	)*dtcld
	1514	if (source.gt.value) then
	1515	factor = value/source
	1516	nraut(i,k) = nraut(i,k)*factor
	1517	nrcol(i,k) = nrcol(i,k)*factor
	1518	niacr(i,k) = niacr(i,k)*factor
	1519	nsacr(i,k) = nsacr(i,k)*factor
	1520	ngacr(i,k) = ngacr(i,k)*factor
	1521	endif
	1522	!
	1523	work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
	1524	! update
	1525	q(i,k) = q(i,k)+work2(i,k)*dtcld
	1526	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
	1527	+paacw(i,k)+paacw(i,k))*dtcld,0.)
	1528	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
	1529	+prevp(i,k)-piacr(i,k)-pgacr(i,k) &
	1530	-psacr(i,k))*dtcld,0.)
	1531	qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
	1532	+psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
	1533	*dtcld,0.)
	1534	qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
	1535	-pgaut(i,k)+piacr(i,k)*delta3 &
	1536	+praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
	1537	-pracs(i,k)(1.-delta2)+psacr(i,k)delta2) &
	1538	*dtcld,0.)
	1539	qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
	1540	+piacr(i,k)*(1.-delta3) &
	1541	+praci(i,k)(1.-delta3)+psacr(i,k)(1.-delta2) &
	1542	+pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
	1543	+pgacr(i,k)+pgacs(i,k))*dtcld,0.)
	1544	ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
	1545	-naacw(i,k)-naacw(i,k))*dtcld,0.)
	1546	ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-niacr(i,k) &
	1547	-nsacr(i,k)-ngacr(i,k))*dtcld,0.)
	1548	xlf = xls-xl(i,k)
	1549	xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
	1550	-xl(i,k)prevp(i,k)-xlf(piacr(i,k)+paacw(i,k) &
	1551	+paacw(i,k)+pgacr(i,k)+psacr(i,k))
	1552	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
	1553	else
	1554	!
	1555	! cloud water
	1556	!
	1557	value = max(qmin,qci(i,k,1))
	1558	source= (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)) &
	1559	*dtcld
	1560	if (source.gt.value) then
	1561	factor = value/source
	1562	praut(i,k) = praut(i,k)*factor
	1563	pracw(i,k) = pracw(i,k)*factor
	1564	paacw(i,k) = paacw(i,k)*factor
	1565	endif
	1566	!
	1567	! rain
	1568	!
	1569	value = max(qmin,qrs(i,k,1))
	1570	source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) &
	1571	-pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
	1572	if (source.gt.value) then
	1573	factor = value/source
	1574	praut(i,k) = praut(i,k)*factor
	1575	prevp(i,k) = prevp(i,k)*factor
	1576	pracw(i,k) = pracw(i,k)*factor
	1577	paacw(i,k) = paacw(i,k)*factor
	1578	pseml(i,k) = pseml(i,k)*factor
	1579	pgeml(i,k) = pgeml(i,k)*factor
	1580	endif
	1581	!
	1582	! snow
	1583	!
	1584	value = max(qcrmin,qrs(i,k,2))
	1585	source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
	1586	if (source.gt.value) then
	1587	factor = value/source
	1588	pgacs(i,k) = pgacs(i,k)*factor
	1589	psevp(i,k) = psevp(i,k)*factor
	1590	pseml(i,k) = pseml(i,k)*factor
	1591	endif
	1592	!
	1593	! graupel
	1594	!
	1595	value = max(qcrmin,qrs(i,k,3))
	1596	source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
	1597	if (source.gt.value) then
	1598	factor = value/source
	1599	pgacs(i,k) = pgacs(i,k)*factor
	1600	pgevp(i,k) = pgevp(i,k)*factor
	1601	pgeml(i,k) = pgeml(i,k)*factor
	1602	endif
	1603	!
	1604	! cloud
	1605	!
	1606	value = max(ncmin,ncr(i,k,2))
	1607	source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+naacw(i,k) &
	1608	+naacw(i,k))*dtcld
	1609	if (source.gt.value) then
	1610	factor = value/source
	1611	nraut(i,k) = nraut(i,k)*factor
	1612	nccol(i,k) = nccol(i,k)*factor
	1613	nracw(i,k) = nracw(i,k)*factor
	1614	naacw(i,k) = naacw(i,k)*factor
	1615	endif
	1616	!
	1617	! rain
	1618	!
	1619	value = max(nrmin,ncr(i,k,3))
	1620	source = (-nraut(i,k)+nrcol(i,k)-nseml(i,k)-ngeml(i,k) &
	1621	)*dtcld
	1622	if (source.gt.value) then
	1623	factor = value/source
	1624	nraut(i,k) = nraut(i,k)*factor
	1625	nrcol(i,k) = nrcol(i,k)*factor
	1626	nseml(i,k) = nseml(i,k)*factor
	1627	ngeml(i,k) = ngeml(i,k)*factor
	1628	endif
	1629	!
	1630	work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
	1631	! update
	1632	q(i,k) = q(i,k)+work2(i,k)*dtcld
	1633	qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
	1634	+paacw(i,k)+paacw(i,k))*dtcld,0.)
	1635	qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
	1636	+prevp(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
	1637	-pgeml(i,k))*dtcld,0.)
	1638	qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
	1639	+pseml(i,k))*dtcld,0.)
	1640	qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
	1641	+pgeml(i,k))*dtcld,0.)
	1642	ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
	1643	-naacw(i,k)-naacw(i,k))*dtcld,0.)
	1644	ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)+nseml(i,k) &
	1645	+ngeml(i,k))*dtcld,0.)
	1646	xlf = xls-xl(i,k)
	1647	xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
	1648	-xlf*(pseml(i,k)+pgeml(i,k))
	1649	t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
	1650	endif
	1651	enddo
	1652	enddo
	1653	!
	1654	! Inline expansion for fpvs
	1655	! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	1656	! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	1657	hsub = xls
	1658	hvap = xlv0
	1659	cvap = cpv
	1660	ttp=t0c+0.01
	1661	dldt=cvap-cliq
	1662	xa=-dldt/rv
	1663	xb=xa+hvap/(rv*ttp)
	1664	dldti=cvap-cice
	1665	xai=-dldti/rv
	1666	xbi=xai+hsub/(rv*ttp)
	1667	do k = kts, kte
	1668	do i = its, ite
	1669	tr=ttp/t(i,k)
	1670	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	1671	qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
	1672	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
	1673	qs(i,k,1) = max(qs(i,k,1),qmin)
	1674	tr=ttp/t(i,k)
	1675	if(t(i,k).lt.ttp) then
	1676	qs(i,k,2)=psatexp(log(tr)(xai))exp(xbi(1.-tr))
	1677	else
	1678	qs(i,k,2)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	1679	endif
	1680	qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
	1681	qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
	1682	qs(i,k,2) = max(qs(i,k,2),qmin)
	1683	rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
	1684	enddo
	1685	enddo
	1686	!
	1687	call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
	1688	rslope3,work1,workn,its,ite,kts,kte)
	1689	do k = kts, kte
	1690	do i = its, ite
	1691	!-----------------------------------------------------------------
	1692	! re-compute the mean-volume drop diameter [LH A10]
	1693	! for raindrop distribution
	1694	!-----------------------------------------------------------------
	1695	avedia(i,k,2) = rslope(i,k,1)((24.)*(.3333333))
	1696	!----------------------------------------------------------------
	1697	! Nrevp_s: evaporation/condensation rate of rain [LH A14]
	1698	! (NR->NC)
	1699	!----------------------------------------------------------------
	1700	if(avedia(i,k,2).le.di82) then
	1701	ncr(i,k,2) = ncr(i,k,2)+ncr(i,k,3)
	1702	ncr(i,k,3) = 0.
	1703	!----------------------------------------------------------------
	1704	! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
	1705	! (QR->QC)
	1706	!----------------------------------------------------------------
	1707	qci(i,k,1) = qci(i,k,1)+qrs(i,k,1)
	1708	qrs(i,k,1) = 0.
	1709	endif
	1710	enddo
	1711	enddo
	1712	!
	1713	do k = kts, kte
	1714	do i = its, ite
	1715	!---------------------------------------------------------------
	1716	! rate of change of cloud drop concentration due to CCN activation
	1717	! pcact: QV -> QC [LH 8] [KK 14]
	1718	! ncact: NCCN -> NC [LH A2] [KK 12]
	1719	!---------------------------------------------------------------
	1720	if(rh(i,k,1).gt.1.) then
	1721	ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2)) &
	1722	min(1.,(rh(i,k,1)/satmax)*actk) - ncr(i,k,2)))/dtcld
	1723	ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
	1724	pcact(i,k) = min(4.pidenr(actr1.E-6)*3ncact(i,k)/ &
	1725	(3.*den(i,k)),max(q(i,k),0.)/dtcld)
	1726	q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
	1727	qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
	1728	ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
	1729	ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
	1730	t(i,k) = t(i,k)+pcact(i,k)xl(i,k)/cpm(i,k)dtcld
	1731	endif
	1732	!---------------------------------------------------------------
	1733	! pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
	1734	! if there exists additional water vapor condensated/if
	1735	! evaporation of cloud water is not enough to remove subsaturation
	1736	! (QV->QC or QC->QV)
	1737	!---------------------------------------------------------------
	1738	tr=ttp/t(i,k)
	1739	qs(i,k,1)=psatexp(log(tr)(xa))exp(xb(1.-tr))
	1740	qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
	1741	qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
	1742	qs(i,k,1) = max(qs(i,k,1),qmin)
	1743	work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
	1744	work2(i,k) = qci(i,k,1)+work1(i,k,1)
	1745	pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
	1746	if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.) &
	1747	pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
	1748	!----------------------------------------------------------------
	1749	! ncevp: evpration of Cloud number concentration [LH A3]
	1750	! (NC->NCCN)
	1751	!----------------------------------------------------------------
	1752	if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
	1753	ncr(i,k,2) = 0.
	1754	ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
	1755	endif
	1756	!
	1757	q(i,k) = q(i,k)-pcond(i,k)*dtcld
	1758	qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
	1759	t(i,k) = t(i,k)+pcond(i,k)xl(i,k)/cpm(i,k)dtcld
	1760	enddo
	1761	enddo
	1762	!
	1763	!----------------------------------------------------------------
	1764	! padding for small values
	1765	!
	1766	do k = kts, kte
	1767	do i = its, ite
	1768	if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
	1769	if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
	1770	enddo
	1771	enddo
	1772	enddo ! big loops
	1773	END SUBROUTINE wdm62d
	1774	! ...................................................................
	1775	REAL FUNCTION rgmma(x)
	1776	!-------------------------------------------------------------------
	1777	IMPLICIT NONE
	1778	!-------------------------------------------------------------------
	1779	! rgmma function: use infinite product form
	1780	REAL :: euler
	1781	PARAMETER (euler=0.577215664901532)
	1782	REAL :: x, y
	1783	INTEGER :: i
	1784	if(x.eq.1.)then
	1785	rgmma=0.
	1786	else
	1787	rgmma=xexp(eulerx)
	1788	do i=1,10000
	1789	y=float(i)
	1790	rgmma=rgmma(1.000+x/y)exp(-x/y)
	1791	enddo
	1792	rgmma=1./rgmma
	1793	endif
	1794	END FUNCTION rgmma
	1795	!
	1796	!--------------------------------------------------------------------------
	1797	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
	1798	!--------------------------------------------------------------------------
	1799	IMPLICIT NONE
	1800	!--------------------------------------------------------------------------
	1801	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
	1802	xai,xbi,ttp,tr
	1803	INTEGER ice
	1804	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1805	ttp=t0c+0.01
	1806	dldt=cvap-cliq
	1807	xa=-dldt/rv
	1808	xb=xa+hvap/(rv*ttp)
	1809	dldti=cvap-cice
	1810	xai=-dldti/rv
	1811	xbi=xai+hsub/(rv*ttp)
	1812	tr=ttp/t
	1813	if(t.lt.ttp .and. ice.eq.1) then
	1814	fpvs=psat(trxai)exp(xbi*(1.-tr))
	1815	else
	1816	fpvs=psat(trxa)exp(xb*(1.-tr))
	1817	endif
	1818	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1819	END FUNCTION fpvs
	1820	!-------------------------------------------------------------------
	1821	SUBROUTINE wdm6init(den0,denr,dens,cl,cpv, ccn0, allowed_to_read)
	1822	!-------------------------------------------------------------------
	1823	IMPLICIT NONE
	1824	!-------------------------------------------------------------------
	1825	!.... constants which may not be tunable
	1826	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
	1827	LOGICAL, INTENT(IN) :: allowed_to_read
	1828	!
	1829	pi = 4.*atan(1.)
	1830	xlv1 = cl-cpv
	1831	!
	1832	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
	1833	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
	1834	pidnc = pi*denr/6.
	1835	!
	1836	bvtr1 = 1.+bvtr
	1837	bvtr2 = 2.+bvtr
	1838	bvtr3 = 3.+bvtr
	1839	bvtr4 = 4.+bvtr
	1840	bvtr5 = 5.+bvtr
	1841	bvtr6 = 6.+bvtr
	1842	bvtr7 = 7.+bvtr
	1843	bvtr2o5 = 2.5+.5*bvtr
	1844	bvtr3o5 = 3.5+.5*bvtr
	1845	g1pbr = rgmma(bvtr1)
	1846	g2pbr = rgmma(bvtr2)
	1847	g3pbr = rgmma(bvtr3)
	1848	g4pbr = rgmma(bvtr4) ! 17.837825
	1849	g5pbr = rgmma(bvtr5)
	1850	g6pbr = rgmma(bvtr6)
	1851	g7pbr = rgmma(bvtr7)
	1852	g5pbro2 = rgmma(bvtr2o5)
	1853	g7pbro2 = rgmma(bvtr3o5)
	1854	pvtr = avtr*g5pbr/24.
	1855	pvtrn = avtr*g2pbr
	1856	eacrr = 1.0
	1857	pacrr = pin0ravtrg3pbr.25*eacrr
	1858	precr1 = 2.pi1.56
	1859	precr2 = 2.pi.31avtr.5g7pbro2
	1860	pidn0r = pidenrn0r
	1861	pidnr = 4.pidenr
	1862	!
	1863	xmmax = (dimax/dicon)**2
	1864	roqimax = 2.08e22dimax*8
	1865	!
	1866	bvts1 = 1.+bvts
	1867	bvts2 = 2.5+.5*bvts
	1868	bvts3 = 3.+bvts
	1869	bvts4 = 4.+bvts
	1870	g1pbs = rgmma(bvts1) !.8875
	1871	g3pbs = rgmma(bvts3)
	1872	g4pbs = rgmma(bvts4) ! 12.0786
	1873	g5pbso2 = rgmma(bvts2)
	1874	pvts = avts*g4pbs/6.
	1875	pacrs = pin0savtsg3pbs.25
	1876	precs1 = 4.n0s.65
	1877	precs2 = 4.n0s.44avts.5g5pbso2
	1878	pidn0s = pidensn0s
	1879	!
	1880	pacrc = pin0savtsg3pbs.25*eacrc
	1881	!
	1882	bvtg1 = 1.+bvtg
	1883	bvtg2 = 2.5+.5*bvtg
	1884	bvtg3 = 3.+bvtg
	1885	bvtg4 = 4.+bvtg
	1886	g1pbg = rgmma(bvtg1)
	1887	g3pbg = rgmma(bvtg3)
	1888	g4pbg = rgmma(bvtg4)
	1889	g5pbgo2 = rgmma(bvtg2)
	1890	pacrg = pin0gavtgg3pbg.25
	1891	pvtg = avtg*g4pbg/6.
	1892	precg1 = 2.pin0g*.78
	1893	precg2 = 2.pin0g.31avtg*.5g5pbgo2
	1894	pidn0g = pidengn0g
	1895	!
	1896	rslopecmax = 1./lamdacmax
	1897	rslopermax = 1./lamdarmax
	1898	rslopesmax = 1./lamdasmax
	1899	rslopegmax = 1./lamdagmax
	1900	rsloperbmax = rslopermax ** bvtr
	1901	rslopesbmax = rslopesmax ** bvts
	1902	rslopegbmax = rslopegmax ** bvtg
	1903	rslopec2max = rslopecmax * rslopecmax
	1904	rsloper2max = rslopermax * rslopermax
	1905	rslopes2max = rslopesmax * rslopesmax
	1906	rslopeg2max = rslopegmax * rslopegmax
	1907	rslopec3max = rslopec2max * rslopecmax
	1908	rsloper3max = rsloper2max * rslopermax
	1909	rslopes3max = rslopes2max * rslopesmax
	1910	rslopeg3max = rslopeg2max * rslopegmax
	1911	!
	1912	END SUBROUTINE wdm6init
	1913	!------------------------------------------------------------------------------
	1914	subroutine slope_wdm6(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
	1915	vt,vtn,its,ite,kts,kte)
	1916	IMPLICIT NONE
	1917	INTEGER :: its,ite, jts,jte, kts,kte
	1918	REAL, DIMENSION( its:ite , kts:kte,3) :: &
	1919	qrs, &
	1920	rslope, &
	1921	rslopeb, &
	1922	rslope2, &
	1923	rslope3, &
	1924	vt
	1925	REAL, DIMENSION( its:ite , kts:kte) :: &
	1926	ncr, &
	1927	vtn, &
	1928	den, &
	1929	denfac, &
	1930	t
	1931	REAL, PARAMETER :: t0c = 273.15
	1932	REAL, DIMENSION( its:ite , kts:kte ) :: &
	1933	n0sfac
	1934	REAL :: lamdar, lamdas, lamdag, x, y, z, supcol
	1935	integer :: i, j, k
	1936	!----------------------------------------------------------------
	1937	! size distributions: (x=mixing ratio, y=air density):
	1938	! valid for mixing ratio > 1.e-9 kg/kg.
	1939	!
	1940	! Optimizatin : A*B => exp(log(A)(B))
	1941	lamdar(x,y,z)= exp(log(((pidnrz)/(xy)))*((.33333333)))
	1942	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
	1943	lamdag(x,y)= sqrt(sqrt(pidn0g/(xy))) ! (pidn0g/(xy))**.25
	1944	!
	1945	do k = kts, kte
	1946	do i = its, ite
	1947	supcol = t0c-t(i,k)
	1948	!---------------------------------------------------------------
	1949	! n0s: Intercept parameter for snow [m-4] [HDC 6]
	1950	!---------------------------------------------------------------
	1951	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
	1952	if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
	1953	rslope(i,k,1) = rslopermax
	1954	rslopeb(i,k,1) = rsloperbmax
	1955	rslope2(i,k,1) = rsloper2max
	1956	rslope3(i,k,1) = rsloper3max
	1957	else
	1958	rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
	1959	rslopeb(i,k,1) = rslope(i,k,1)**bvtr
	1960	rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
	1961	rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
	1962	endif
	1963	if(qrs(i,k,2).le.qcrmin) then
	1964	rslope(i,k,2) = rslopesmax
	1965	rslopeb(i,k,2) = rslopesbmax
	1966	rslope2(i,k,2) = rslopes2max
	1967	rslope3(i,k,2) = rslopes3max
	1968	else
	1969	rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
	1970	rslopeb(i,k,2) = rslope(i,k,2)**bvts
	1971	rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
	1972	rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
	1973	endif
	1974	if(qrs(i,k,3).le.qcrmin) then
	1975	rslope(i,k,3) = rslopegmax
	1976	rslopeb(i,k,3) = rslopegbmax
	1977	rslope2(i,k,3) = rslopeg2max
	1978	rslope3(i,k,3) = rslopeg3max
	1979	else
	1980	rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
	1981	rslopeb(i,k,3) = rslope(i,k,3)**bvtg
	1982	rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
	1983	rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
	1984	endif
	1985	vt(i,k,1) = pvtrrslopeb(i,k,1)denfac(i,k)
	1986	vt(i,k,2) = pvtsrslopeb(i,k,2)denfac(i,k)
	1987	vt(i,k,3) = pvtgrslopeb(i,k,3)denfac(i,k)
	1988	vtn(i,k) = pvtrnrslopeb(i,k,1)denfac(i,k)
	1989	if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
	1990	if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
	1991	if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
	1992	if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
	1993	enddo
	1994	enddo
	1995	END subroutine slope_wdm6
	1996	!-----------------------------------------------------------------------------
	1997	subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
	1998	vt,vtn,its,ite,kts,kte)
	1999	IMPLICIT NONE
	2000	INTEGER :: its,ite, jts,jte, kts,kte
	2001	REAL, DIMENSION( its:ite , kts:kte) :: &
	2002	qrs, &
	2003	ncr, &
	2004	rslope, &
	2005	rslopeb, &
	2006	rslope2, &
	2007	rslope3, &
	2008	vt, &
	2009	vtn, &
	2010	den, &
	2011	denfac, &
	2012	t
	2013	REAL, PARAMETER :: t0c = 273.15
	2014	REAL, DIMENSION( its:ite , kts:kte ) :: &
	2015	n0sfac
	2016	REAL :: lamdar, x, y, z, supcol
	2017	integer :: i, j, k
	2018	!----------------------------------------------------------------
	2019	! size distributions: (x=mixing ratio, y=air density):
	2020	! valid for mixing ratio > 1.e-9 kg/kg.
	2021	lamdar(x,y,z)= exp(log(((pidnrz)/(xy)))*((.33333333)))
	2022	!
	2023	do k = kts, kte
	2024	do i = its, ite
	2025	if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
	2026	rslope(i,k) = rslopermax
	2027	rslopeb(i,k) = rsloperbmax
	2028	rslope2(i,k) = rsloper2max
	2029	rslope3(i,k) = rsloper3max
	2030	else
	2031	rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
	2032	rslopeb(i,k) = rslope(i,k)**bvtr
	2033	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	2034	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	2035	endif
	2036	vt(i,k) = pvtrrslopeb(i,k)denfac(i,k)
	2037	vtn(i,k) = pvtrnrslopeb(i,k)denfac(i,k)
	2038	if(qrs(i,k).le.0.0) vt(i,k) = 0.0
	2039	if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
	2040	enddo
	2041	enddo
	2042	END subroutine slope_rain
	2043	!------------------------------------------------------------------------------
	2044	subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
	2045	vt,its,ite,kts,kte)
	2046	IMPLICIT NONE
	2047	INTEGER :: its,ite, jts,jte, kts,kte
	2048	REAL, DIMENSION( its:ite , kts:kte) :: &
	2049	qrs, &
	2050	rslope, &
	2051	rslopeb, &
	2052	rslope2, &
	2053	rslope3, &
	2054	vt, &
	2055	den, &
	2056	denfac, &
	2057	t
	2058	REAL, PARAMETER :: t0c = 273.15
	2059	REAL, DIMENSION( its:ite , kts:kte ) :: &
	2060	n0sfac
	2061	REAL :: lamdas, x, y, z, supcol
	2062	integer :: i, j, k
	2063	!----------------------------------------------------------------
	2064	! size distributions: (x=mixing ratio, y=air density):
	2065	! valid for mixing ratio > 1.e-9 kg/kg.
	2066	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
	2067	!
	2068	do k = kts, kte
	2069	do i = its, ite
	2070	supcol = t0c-t(i,k)
	2071	!---------------------------------------------------------------
	2072	! n0s: Intercept parameter for snow [m-4] [HDC 6]
	2073	!---------------------------------------------------------------
	2074	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
	2075	if(qrs(i,k).le.qcrmin)then
	2076	rslope(i,k) = rslopesmax
	2077	rslopeb(i,k) = rslopesbmax
	2078	rslope2(i,k) = rslopes2max
	2079	rslope3(i,k) = rslopes3max
	2080	else
	2081	rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
	2082	rslopeb(i,k) = rslope(i,k)**bvts
	2083	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	2084	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	2085	endif
	2086	vt(i,k) = pvtsrslopeb(i,k)denfac(i,k)
	2087	if(qrs(i,k).le.0.0) vt(i,k) = 0.0
	2088	enddo
	2089	enddo
	2090	END subroutine slope_snow
	2091	!----------------------------------------------------------------------------------
	2092	subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
	2093	vt,its,ite,kts,kte)
	2094	IMPLICIT NONE
	2095	INTEGER :: its,ite, jts,jte, kts,kte
	2096	REAL, DIMENSION( its:ite , kts:kte) :: &
	2097	qrs, &
	2098	rslope, &
	2099	rslopeb, &
	2100	rslope2, &
	2101	rslope3, &
	2102	vt, &
	2103	den, &
	2104	denfac, &
	2105	t
	2106	REAL, PARAMETER :: t0c = 273.15
	2107	REAL, DIMENSION( its:ite , kts:kte ) :: &
	2108	n0sfac
	2109	REAL :: lamdag, x, y, z, supcol
	2110	integer :: i, j, k
	2111	!----------------------------------------------------------------
	2112	! size distributions: (x=mixing ratio, y=air density):
	2113	! valid for mixing ratio > 1.e-9 kg/kg.
	2114	lamdag(x,y)= sqrt(sqrt(pidn0g/(xy))) ! (pidn0g/(xy))**.25
	2115	!
	2116	do k = kts, kte
	2117	do i = its, ite
	2118	!---------------------------------------------------------------
	2119	! n0s: Intercept parameter for snow [m-4] [HDC 6]
	2120	!---------------------------------------------------------------
	2121	if(qrs(i,k).le.qcrmin)then
	2122	rslope(i,k) = rslopegmax
	2123	rslopeb(i,k) = rslopegbmax
	2124	rslope2(i,k) = rslopeg2max
	2125	rslope3(i,k) = rslopeg3max
	2126	else
	2127	rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
	2128	rslopeb(i,k) = rslope(i,k)**bvtg
	2129	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	2130	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	2131	endif
	2132	vt(i,k) = pvtgrslopeb(i,k)denfac(i,k)
	2133	if(qrs(i,k).le.0.0) vt(i,k) = 0.0
	2134	enddo
	2135	enddo
	2136	END subroutine slope_graup
	2137	!---------------------------------------------------------------------------------
	2138	!-------------------------------------------------------------------
	2139	SUBROUTINE nislfv_rain_plmr(im,km,denl,denfacl,tkl,dzl,wwl,rql,rnl,precip,dt,id,iter,rid)
	2140	!-------------------------------------------------------------------
	2141	!
	2142	! for non-iteration semi-Lagrangain forward advection for cloud
	2143	! with mass conservation and positive definite advection
	2144	! 2nd order interpolation with monotonic piecewise linear method
	2145	! this routine is under assumption of decfl < 1 for semi_Lagrangian
	2146	!
	2147	! dzl depth of model layer in meter
	2148	! wwl terminal velocity at model layer m/s
	2149	! rql cloud density*mixing ration
	2150	! precip precipitation
	2151	! dt time step
	2152	! id kind of precip: 0 test case; 1 raindrop
	2153	! iter how many time to guess mean terminal velocity: 0 pure forward.
	2154	! 0 : use departure wind for advection
	2155	! 1 : use mean wind for advection
	2156	! > 1 : use mean wind after iter-1 iterations
	2157	! rid : 1 for number 0 for mixing ratio
	2158	!
	2159	! author: hann-ming henry juang <henry.juang@noaa.gov>
	2160	! implemented by song-you hong
	2161	!
	2162	implicit none
	2163	integer im,km,id
	2164	real dt
	2165	real dzl(im,km),wwl(im,km),rql(im,km),rnl(im,km),precip(im)
	2166	real denl(im,km),denfacl(im,km),tkl(im,km)
	2167	!
	2168	integer i,k,n,m,kk,kb,kt,iter,rid
	2169	real tl,tl2,qql,dql,qqd
	2170	real th,th2,qqh,dqh
	2171	real zsum,qsum,dim,dip,c1,con1,fa1,fa2
	2172	real allold, allnew, zz, dzamin, cflmax, decfl
	2173	real dz(km), ww(km), qq(km), nr(km), wd(km), wa(km), wa2(km), was(km)
	2174	real den(km), denfac(km), tk(km)
	2175	real wi(km+1), zi(km+1), za(km+1)
	2176	real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
	2177	real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
	2178	!
	2179	precip(:) = 0.0
	2180	!
	2181	i_loop : do i=1,im
	2182	! -----------------------------------
	2183	dz(:) = dzl(i,:)
	2184	qq(:) = rql(i,:)
	2185	nr(:) = rnl(i,:)
	2186	if(rid .eq. 1) nr(:) = rnl(i,:)/denl(i,:)
	2187	ww(:) = wwl(i,:)
	2188	den(:) = denl(i,:)
	2189	denfac(:) = denfacl(i,:)
	2190	tk(:) = tkl(i,:)
	2191	! skip for no precipitation for all layers
	2192	allold = 0.0
	2193	do k=1,km
	2194	allold = allold + qq(k)
	2195	enddo
	2196	if(allold.le.0.0) then
	2197	cycle i_loop
	2198	endif
	2199	!
	2200	! compute interface values
	2201	zi(1)=0.0
	2202	do k=1,km
	2203	zi(k+1) = zi(k)+dz(k)
	2204	enddo
	2205	!
	2206	! save departure wind
	2207	wd(:) = ww(:)
	2208	n=1
	2209	100 continue
	2210	! plm is 2nd order, we can use 2nd order wi or 3rd order wi
	2211	! 2nd order interpolation to get wi
	2212	wi(1) = ww(1)
	2213	wi(km+1) = ww(km)
	2214	do k=2,km
	2215	wi(k) = (ww(k)dz(k-1)+ww(k-1)dz(k))/(dz(k-1)+dz(k))
	2216	enddo
	2217	! 3rd order interpolation to get wi
	2218	fa1 = 9./16.
	2219	fa2 = 1./16.
	2220	wi(1) = ww(1)
	2221	wi(2) = 0.5*(ww(2)+ww(1))
	2222	do k=3,km-1
	2223	wi(k) = fa1(ww(k)+ww(k-1))-fa2(ww(k+1)+ww(k-2))
	2224	enddo
	2225	wi(km) = 0.5*(ww(km)+ww(km-1))
	2226	wi(km+1) = ww(km)
	2227	!
	2228	! terminate of top of raingroup
	2229	do k=2,km
	2230	if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
	2231	enddo
	2232	!
	2233	! diffusivity of wi
	2234	con1 = 0.05
	2235	do k=km,1,-1
	2236	decfl = (wi(k+1)-wi(k))*dt/dz(k)
	2237	if( decfl .gt. con1 ) then
	2238	wi(k) = wi(k+1) - con1*dz(k)/dt
	2239	endif
	2240	enddo
	2241	! compute arrival point
	2242	do k=1,km+1
	2243	za(k) = zi(k) - wi(k)*dt
	2244	enddo
	2245	!
	2246	do k=1,km
	2247	dza(k) = za(k+1)-za(k)
	2248	enddo
	2249	dza(km+1) = zi(km+1) - za(km+1)
	2250	!
	2251	! computer deformation at arrival point
	2252	do k=1,km
	2253	qa(k) = qq(k)*dz(k)/dza(k)
	2254	qr(k) = qa(k)/den(k)
	2255	if(rid .eq. 1) qr(k) = qa(K)
	2256	enddo
	2257	qa(km+1) = 0.0
	2258	! call maxmin(km,1,qa,' arrival points ')
	2259	!
	2260	! compute arrival terminal velocity, and estimate mean terminal velocity
	2261	! then back to use mean terminal velocity
	2262	if( n.le.iter ) then
	2263	if(rid.eq.1) then
	2264	call slope_rain(nr,qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
	2265	else
	2266	call slope_rain(qr,nr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
	2267	endif
	2268	if(rid.eq.1) wa(:) = wa2(:)
	2269	if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
	2270	do k=1,km
	2271	!#ifdef DEBUG
	2272	! print,' slope_wsm3 ',qr(k)1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
	2273	!#endif
	2274	! mean wind is average of departure and new arrival winds
	2275	ww(k) = 0.5* ( wd(k)+wa(k) )
	2276	enddo
	2277	was(:) = wa(:)
	2278	n=n+1
	2279	go to 100
	2280	endif
	2281	!
	2282	! estimate values at arrival cell interface with monotone
	2283	do k=2,km
	2284	dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
	2285	dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
	2286	if( dip*dim.le.0.0 ) then
	2287	qmi(k)=qa(k)
	2288	qpi(k)=qa(k)
	2289	else
	2290	qpi(k)=qa(k)+0.5(dip+dim)dza(k)
	2291	qmi(k)=2.0*qa(k)-qpi(k)
	2292	if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
	2293	qpi(k) = qa(k)
	2294	qmi(k) = qa(k)
	2295	endif
	2296	endif
	2297	enddo
	2298	qpi(1)=qa(1)
	2299	qmi(1)=qa(1)
	2300	qmi(km+1)=qa(km+1)
	2301	qpi(km+1)=qa(km+1)
	2302	!
	2303	! interpolation to regular point
	2304	qn = 0.0
	2305	kb=1
	2306	kt=1
	2307	intp : do k=1,km
	2308	kb=max(kb-1,1)
	2309	kt=max(kt-1,1)
	2310	! find kb and kt
	2311	if( zi(k).ge.za(km+1) ) then
	2312	exit intp
	2313	else
	2314	find_kb : do kk=kb,km
	2315	if( zi(k).le.za(kk+1) ) then
	2316	kb = kk
	2317	exit find_kb
	2318	else
	2319	cycle find_kb
	2320	endif
	2321	enddo find_kb
	2322	find_kt : do kk=kt,km
	2323	if( zi(k+1).le.za(kk) ) then
	2324	kt = kk
	2325	exit find_kt
	2326	else
	2327	cycle find_kt
	2328	endif
	2329	enddo find_kt
	2330	kt = kt - 1
	2331	! compute q with piecewise constant method
	2332	if( kt.eq.kb ) then
	2333	tl=(zi(k)-za(kb))/dza(kb)
	2334	th=(zi(k+1)-za(kb))/dza(kb)
	2335	tl2=tl*tl
	2336	th2=th*th
	2337	qqd=0.5*(qpi(kb)-qmi(kb))
	2338	qqh=qqdth2+qmi(kb)th
	2339	qql=qqdtl2+qmi(kb)tl
	2340	qn(k) = (qqh-qql)/(th-tl)
	2341	else if( kt.gt.kb ) then
	2342	tl=(zi(k)-za(kb))/dza(kb)
	2343	tl2=tl*tl
	2344	qqd=0.5*(qpi(kb)-qmi(kb))
	2345	qql=qqdtl2+qmi(kb)tl
	2346	dql = qa(kb)-qql
	2347	zsum = (1.-tl)*dza(kb)
	2348	qsum = dql*dza(kb)
	2349	if( kt-kb.gt.1 ) then
	2350	do m=kb+1,kt-1
	2351	zsum = zsum + dza(m)
	2352	qsum = qsum + qa(m) * dza(m)
	2353	enddo
	2354	endif
	2355	th=(zi(k+1)-za(kt))/dza(kt)
	2356	th2=th*th
	2357	qqd=0.5*(qpi(kt)-qmi(kt))
	2358	dqh=qqdth2+qmi(kt)th
	2359	zsum = zsum + th*dza(kt)
	2360	qsum = qsum + dqh*dza(kt)
	2361	qn(k) = qsum/zsum
	2362	endif
	2363	cycle intp
	2364	endif
	2365	!
	2366	enddo intp
	2367	!
	2368	! rain out
	2369	sum_precip: do k=1,km
	2370	if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
	2371	precip(i) = precip(i) + qa(k)*dza(k)
	2372	cycle sum_precip
	2373	else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
	2374	precip(i) = precip(i) + qa(k)*(0.0-za(k))
	2375	exit sum_precip
	2376	endif
	2377	exit sum_precip
	2378	enddo sum_precip
	2379	!
	2380	! replace the new values
	2381	rql(i,:) = qn(:)
	2382	!
	2383	! ----------------------------------
	2384	enddo i_loop
	2385	!
	2386	END SUBROUTINE nislfv_rain_plmr
	2387	!-------------------------------------------------------------------
	2388	SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
	2389	!-------------------------------------------------------------------
	2390	!
	2391	! for non-iteration semi-Lagrangain forward advection for cloud
	2392	! with mass conservation and positive definite advection
	2393	! 2nd order interpolation with monotonic piecewise linear method
	2394	! this routine is under assumption of decfl < 1 for semi_Lagrangian
	2395	!
	2396	! dzl depth of model layer in meter
	2397	! wwl terminal velocity at model layer m/s
	2398	! rql cloud density*mixing ration
	2399	! precip precipitation
	2400	! dt time step
	2401	! id kind of precip: 0 test case; 1 raindrop
	2402	! iter how many time to guess mean terminal velocity: 0 pure forward.
	2403	! 0 : use departure wind for advection
	2404	! 1 : use mean wind for advection
	2405	! > 1 : use mean wind after iter-1 iterations
	2406	!
	2407	! author: hann-ming henry juang <henry.juang@noaa.gov>
	2408	! implemented by song-you hong
	2409	!
	2410	implicit none
	2411	integer im,km,id
	2412	real dt
	2413	real dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
	2414	real denl(im,km),denfacl(im,km),tkl(im,km)
	2415	!
	2416	integer i,k,n,m,kk,kb,kt,iter,ist
	2417	real tl,tl2,qql,dql,qqd
	2418	real th,th2,qqh,dqh
	2419	real zsum,qsum,dim,dip,c1,con1,fa1,fa2
	2420	real allold, allnew, zz, dzamin, cflmax, decfl
	2421	real dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
	2422	real den(km), denfac(km), tk(km)
	2423	real wi(km+1), zi(km+1), za(km+1)
	2424	real qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
	2425	real dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
	2426	!
	2427	precip(:) = 0.0
	2428	precip1(:) = 0.0
	2429	precip2(:) = 0.0
	2430	!
	2431	i_loop : do i=1,im
	2432	! -----------------------------------
	2433	dz(:) = dzl(i,:)
	2434	qq(:) = rql(i,:)
	2435	qq2(:) = rql2(i,:)
	2436	ww(:) = wwl(i,:)
	2437	den(:) = denl(i,:)
	2438	denfac(:) = denfacl(i,:)
	2439	tk(:) = tkl(i,:)
	2440	! skip for no precipitation for all layers
	2441	allold = 0.0
	2442	do k=1,km
	2443	allold = allold + qq(k)
	2444	enddo
	2445	if(allold.le.0.0) then
	2446	cycle i_loop
	2447	endif
	2448	!
	2449	! compute interface values
	2450	zi(1)=0.0
	2451	do k=1,km
	2452	zi(k+1) = zi(k)+dz(k)
	2453	enddo
	2454	!
	2455	! save departure wind
	2456	wd(:) = ww(:)
	2457	n=1
	2458	100 continue
	2459	! plm is 2nd order, we can use 2nd order wi or 3rd order wi
	2460	! 2nd order interpolation to get wi
	2461	wi(1) = ww(1)
	2462	wi(km+1) = ww(km)
	2463	do k=2,km
	2464	wi(k) = (ww(k)dz(k-1)+ww(k-1)dz(k))/(dz(k-1)+dz(k))
	2465	enddo
	2466	! 3rd order interpolation to get wi
	2467	fa1 = 9./16.
	2468	fa2 = 1./16.
	2469	wi(1) = ww(1)
	2470	wi(2) = 0.5*(ww(2)+ww(1))
	2471	do k=3,km-1
	2472	wi(k) = fa1(ww(k)+ww(k-1))-fa2(ww(k+1)+ww(k-2))
	2473	enddo
	2474	wi(km) = 0.5*(ww(km)+ww(km-1))
	2475	wi(km+1) = ww(km)
	2476	!
	2477	! terminate of top of raingroup
	2478	do k=2,km
	2479	if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
	2480	enddo
	2481	!
	2482	! diffusivity of wi
	2483	con1 = 0.05
	2484	do k=km,1,-1
	2485	decfl = (wi(k+1)-wi(k))*dt/dz(k)
	2486	if( decfl .gt. con1 ) then
	2487	wi(k) = wi(k+1) - con1*dz(k)/dt
	2488	endif
	2489	enddo
	2490	! compute arrival point
	2491	do k=1,km+1
	2492	za(k) = zi(k) - wi(k)*dt
	2493	enddo
	2494	!
	2495	do k=1,km
	2496	dza(k) = za(k+1)-za(k)
	2497	enddo
	2498	dza(km+1) = zi(km+1) - za(km+1)
	2499	!
	2500	! computer deformation at arrival point
	2501	do k=1,km
	2502	qa(k) = qq(k)*dz(k)/dza(k)
	2503	qa2(k) = qq2(k)*dz(k)/dza(k)
	2504	qr(k) = qa(k)/den(k)
	2505	qr2(k) = qa2(k)/den(k)
	2506	enddo
	2507	qa(km+1) = 0.0
	2508	qa2(km+1) = 0.0
	2509	! call maxmin(km,1,qa,' arrival points ')
	2510	!
	2511	! compute arrival terminal velocity, and estimate mean terminal velocity
	2512	! then back to use mean terminal velocity
	2513	if( n.le.iter ) then
	2514	call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
	2515	call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
	2516	do k = 1, km
	2517	tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
	2518	IF ( tmp(k) .gt. 1.e-15 ) THEN
	2519	wa(k) = (wa(k)qr(k) + wa2(k)qr2(k))/tmp(k)
	2520	ELSE
	2521	wa(k) = 0.
	2522	ENDIF
	2523	enddo
	2524	if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
	2525	do k=1,km
	2526	!#ifdef DEBUG
	2527	! print,' slope_wsm3 ',qr(k)1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
	2528	! ww(k),wa(k)
	2529	!#endif
	2530	! mean wind is average of departure and new arrival winds
	2531	ww(k) = 0.5* ( wd(k)+wa(k) )
	2532	enddo
	2533	was(:) = wa(:)
	2534	n=n+1
	2535	go to 100
	2536	endif
	2537	ist_loop : do ist = 1, 2
	2538	if (ist.eq.2) then
	2539	qa(:) = qa2(:)
	2540	endif
	2541	!
	2542	precip(i) = 0.
	2543	!
	2544	! estimate values at arrival cell interface with monotone
	2545	do k=2,km
	2546	dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
	2547	dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
	2548	if( dip*dim.le.0.0 ) then
	2549	qmi(k)=qa(k)
	2550	qpi(k)=qa(k)
	2551	else
	2552	qpi(k)=qa(k)+0.5(dip+dim)dza(k)
	2553	qmi(k)=2.0*qa(k)-qpi(k)
	2554	if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
	2555	qpi(k) = qa(k)
	2556	qmi(k) = qa(k)
	2557	endif
	2558	endif
	2559	enddo
	2560	qpi(1)=qa(1)
	2561	qmi(1)=qa(1)
	2562	qmi(km+1)=qa(km+1)
	2563	qpi(km+1)=qa(km+1)
	2564	!
	2565	! interpolation to regular point
	2566	qn = 0.0
	2567	kb=1
	2568	kt=1
	2569	intp : do k=1,km
	2570	kb=max(kb-1,1)
	2571	kt=max(kt-1,1)
	2572	! find kb and kt
	2573	if( zi(k).ge.za(km+1) ) then
	2574	exit intp
	2575	else
	2576	find_kb : do kk=kb,km
	2577	if( zi(k).le.za(kk+1) ) then
	2578	kb = kk
	2579	exit find_kb
	2580	else
	2581	cycle find_kb
	2582	endif
	2583	enddo find_kb
	2584	find_kt : do kk=kt,km
	2585	if( zi(k+1).le.za(kk) ) then
	2586	kt = kk
	2587	exit find_kt
	2588	else
	2589	cycle find_kt
	2590	endif
	2591	enddo find_kt
	2592	kt = kt - 1
	2593	! compute q with piecewise constant method
	2594	if( kt.eq.kb ) then
	2595	tl=(zi(k)-za(kb))/dza(kb)
	2596	th=(zi(k+1)-za(kb))/dza(kb)
	2597	tl2=tl*tl
	2598	th2=th*th
	2599	qqd=0.5*(qpi(kb)-qmi(kb))
	2600	qqh=qqdth2+qmi(kb)th
	2601	qql=qqdtl2+qmi(kb)tl
	2602	qn(k) = (qqh-qql)/(th-tl)
	2603	else if( kt.gt.kb ) then
	2604	tl=(zi(k)-za(kb))/dza(kb)
	2605	tl2=tl*tl
	2606	qqd=0.5*(qpi(kb)-qmi(kb))
	2607	qql=qqdtl2+qmi(kb)tl
	2608	dql = qa(kb)-qql
	2609	zsum = (1.-tl)*dza(kb)
	2610	qsum = dql*dza(kb)
	2611	if( kt-kb.gt.1 ) then
	2612	do m=kb+1,kt-1
	2613	zsum = zsum + dza(m)
	2614	qsum = qsum + qa(m) * dza(m)
	2615	enddo
	2616	endif
	2617	th=(zi(k+1)-za(kt))/dza(kt)
	2618	th2=th*th
	2619	qqd=0.5*(qpi(kt)-qmi(kt))
	2620	dqh=qqdth2+qmi(kt)th
	2621	zsum = zsum + th*dza(kt)
	2622	qsum = qsum + dqh*dza(kt)
	2623	qn(k) = qsum/zsum
	2624	endif
	2625	cycle intp
	2626	endif
	2627	!
	2628	enddo intp
	2629	!
	2630	! rain out
	2631	sum_precip: do k=1,km
	2632	if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
	2633	precip(i) = precip(i) + qa(k)*dza(k)
	2634	cycle sum_precip
	2635	else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
	2636	precip(i) = precip(i) + qa(k)*(0.0-za(k))
	2637	exit sum_precip
	2638	endif
	2639	exit sum_precip
	2640	enddo sum_precip
	2641	!
	2642	! replace the new values
	2643	if(ist.eq.1) then
	2644	rql(i,:) = qn(:)
	2645	precip1(i) = precip(i)
	2646	else
	2647	rql2(i,:) = qn(:)
	2648	precip2(i) = precip(i)
	2649	endif
	2650	enddo ist_loop
	2651	!
	2652	! ----------------------------------
	2653	enddo i_loop
	2654	!
	2655	END SUBROUTINE nislfv_rain_plm6
	2656	END MODULE module_mp_wdm6

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

Download in other formats:

Original Format