Context Navigation

module_mp_wsm5.F

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

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

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

Download in other formats:

Original Format