Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_wsm3.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: 56.1 KB

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

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

Download in other formats:

Original Format