Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_wsm5.F @ 13

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 68.5 KB

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

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

Download in other formats:

Original Format