Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_mp_wsm3.F @ 3094

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

Rev	Line
[2759]	1	#if ( RWORDSIZE == 4 )
	2	# define VREC vsrec
	3	# define VSQRT vssqrt
	4	#else
	5	# define VREC vrec
	6	# define VSQRT vsqrt
	7	#endif
	8
	9	MODULE module_mp_wsm3
	10	!
	11	!
	12	REAL, PARAMETER, PRIVATE :: dtcldcr = 120.
	13	REAL, PARAMETER, PRIVATE :: n0r = 8.e6
	14	REAL, PARAMETER, PRIVATE :: avtr = 841.9
	15	REAL, PARAMETER, PRIVATE :: bvtr = 0.8
	16	REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
	17	REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
	18	REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
	19	REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
	20	REAL, PARAMETER, PRIVATE :: avts = 11.72
	21	REAL, PARAMETER, PRIVATE :: bvts = .41
	22	REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! t=-90C unlimited
	23	REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4
	24	REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5
	25	REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4
	26	REAL, PARAMETER, PRIVATE :: betai = .6
	27	REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
	28	REAL, PARAMETER, PRIVATE :: dicon = 11.9
	29	REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6
	30	REAL, PARAMETER, PRIVATE :: dimax = 500.e-6
	31	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
	32	REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
	33	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9
	34	REAL, SAVE :: &
	35	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
	36	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
	37	precr1,precr2,xm0,xmmax,roqimax,bvts1, &
	38	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
	39	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
	40	pidn0s,xlv1, &
	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	!
	69	! This code is a 3-class simple ice microphyiscs scheme (WSM3) of the WRF
	70	! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
	71	! number concentration is a function of temperature, and seperate assumption
	72	! is developed, in which ice crystal number concentration is a function
	73	! of ice amount. A theoretical background of the ice-microphysics and related
	74	! processes in the WSMMPs are described in Hong et al. (2004).
	75	! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
	76	! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
	77	!
	78	! WSM3 cloud scheme
	79	!
	80	! Coded by Song-You Hong (Yonsei Univ.)
	81	! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
	82	! Summer 2002
	83	!
	84	! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
	85	! Summer 2003
	86	!
	87	! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
	88	! Dudhia (D89, 1989) J. Atmos. Sci.
	89	! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
	90	!
	91	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
	92	ims,ime, jms,jme, kms,kme , &
	93	its,ite, jts,jte, kts,kte
	94	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	95	INTENT(INOUT) :: &
	96	th, &
	97	q, &
	98	qci, &
	99	qrs
	100	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
	101	INTENT(IN ) :: w, &
	102	den, &
	103	pii, &
	104	p, &
	105	delz
	106	REAL, INTENT(IN ) :: delt, &
	107	g, &
	108	rd, &
	109	rv, &
	110	t0c, &
	111	den0, &
	112	cpd, &
	113	cpv, &
	114	ep1, &
	115	ep2, &
	116	qmin, &
	117	XLS, &
	118	XLV0, &
	119	XLF0, &
	120	cliq, &
	121	cice, &
	122	psat, &
	123	denr
	124	REAL, DIMENSION( ims:ime , jms:jme ), &
	125	INTENT(INOUT) :: rain, &
	126	rainncv, &
	127	sr
	128
	129	REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
	130	INTENT(INOUT) :: snow, &
	131	snowncv
	132
	133	! LOCAL VAR
	134	REAL, DIMENSION( its:ite , kts:kte ) :: t
	135	INTEGER :: i,j,k
	136	!-------------------------------------------------------------------
	137	DO j=jts,jte
	138	DO k=kts,kte
	139	DO i=its,ite
	140	t(i,k)=th(i,k,j)*pii(i,k,j)
	141	ENDDO
	142	ENDDO
	143	CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j) &
	144	,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j) &
	145	,p(ims,kms,j), delz(ims,kms,j) &
	146	,delt,g, cpd, cpv, rd, rv, t0c &
	147	,ep1, ep2, qmin &
	148	,XLS, XLV0, XLF0, den0, denr &
	149	,cliq,cice,psat &
	150	,j &
	151	,rain(ims,j), rainncv(ims,j) &
	152	,sr(ims,j) &
	153	,ids,ide, jds,jde, kds,kde &
	154	,ims,ime, jms,jme, kms,kme &
	155	,its,ite, jts,jte, kts,kte &
	156	,snow(ims,j),snowncv(ims,j) &
	157	)
	158	DO K=kts,kte
	159	DO I=its,ite
	160	th(i,k,j)=t(i,k)/pii(i,k,j)
	161	ENDDO
	162	ENDDO
	163	ENDDO
	164	END SUBROUTINE wsm3
	165	!===================================================================
	166	!
	167	SUBROUTINE wsm32D(t, q, qci, qrs,w, den, p, delz &
	168	,delt,g, cpd, cpv, rd, rv, t0c &
	169	,ep1, ep2, qmin &
	170	,XLS, XLV0, XLF0, den0, denr &
	171	,cliq,cice,psat &
	172	,lat &
	173	,rain, rainncv &
	174	,sr &
	175	,ids,ide, jds,jde, kds,kde &
	176	,ims,ime, jms,jme, kms,kme &
	177	,its,ite, jts,jte, kts,kte &
	178	,snow,snowncv &
	179	)
	180	!-------------------------------------------------------------------
	181	IMPLICIT NONE
	182	!-------------------------------------------------------------------
	183	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
	184	ims,ime, jms,jme, kms,kme , &
	185	its,ite, jts,jte, kts,kte, &
	186	lat
	187	REAL, DIMENSION( its:ite , kts:kte ), &
	188	INTENT(INOUT) :: &
	189	t
	190	REAL, DIMENSION( ims:ime , kms:kme ), &
	191	INTENT(INOUT) :: &
	192	q, &
	193	qci, &
	194	qrs
	195	REAL, DIMENSION( ims:ime , kms:kme ), &
	196	INTENT(IN ) :: w, &
	197	den, &
	198	p, &
	199	delz
	200	REAL, INTENT(IN ) :: delt, &
	201	g, &
	202	cpd, &
	203	cpv, &
	204	t0c, &
	205	den0, &
	206	rd, &
	207	rv, &
	208	ep1, &
	209	ep2, &
	210	qmin, &
	211	XLS, &
	212	XLV0, &
	213	XLF0, &
	214	cliq, &
	215	cice, &
	216	psat, &
	217	denr
	218	REAL, DIMENSION( ims:ime ), &
	219	INTENT(INOUT) :: rain, &
	220	rainncv, &
	221	sr
	222
	223	REAL, DIMENSION( ims:ime ), OPTIONAL, &
	224	INTENT(INOUT) :: snow, &
	225	snowncv
	226	! LOCAL VAR
	227	REAL, DIMENSION( its:ite , kts:kte ) :: &
	228	rh, qs, denfac, rslope, rslope2, rslope3, rslopeb, &
	229	pgen, paut, pacr, pisd, pres, pcon, fall, falk, &
	230	xl, cpm, work1, work2, xni, qs0, n0sfac
	231	INTEGER, DIMENSION( its:ite ) :: kwork1, kwork2
	232	REAL, DIMENSION( its:ite , kts:kte ) :: &
	233	falkc, work1c, work2c, fallc
	234	! variables for optimization
	235	REAL, DIMENSION( its:ite ) :: tvec1
	236	INTEGER, DIMENSION( its:ite ) :: mstep, numdt
	237	LOGICAL, DIMENSION( its:ite ) :: flgcld
	238	REAL :: pi, &
	239	cpmcal, xlcal, lamdar, lamdas, diffus, &
	240	viscos, xka, venfac, conden, diffac, &
	241	x, y, z, a, b, c, d, e, &
	242	fallsum, fallsum_qsi, vt2i,vt2s,acrfac, &
	243	qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
	244	snomlt, hold, holdrs, facqci, supcol, coeres, &
	245	supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, &
	246	qimax, diameter, xni0, roqi0, supice
	247	REAL :: holdc, holdci
	248	INTEGER :: i, j, k, mstepmax, &
	249	iprt, latd, lond, loop, loops, ifsat, kk, n
	250	! Temporaries used for inlining fpvs function
	251	REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
	252	!
	253	!=================================================================
	254	! compute internal functions
	255	!
	256	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
	257	xlcal(x) = xlv0-xlv1*(x-t0c)
	258	!----------------------------------------------------------------
	259	! size distributions: (x=mixing ratio, y=air density):
	260	! valid for mixing ratio > 1.e-9 kg/kg.
	261	!
	262	! Optimizatin : A*B => exp(log(A)(B))
	263	lamdar(x,y)= sqrt(sqrt(pidn0r/(xy))) ! (pidn0r/(xy))**.25
	264	lamdas(x,y,z)= sqrt(sqrt(pidn0sz/(xy))) ! (pidn0sz/(xy))**.25
	265	!
	266	!----------------------------------------------------------------
	267	! diffus: diffusion coefficient of the water vapor
	268	! viscos: kinematic viscosity(m2s-1)
	269	!
	270	diffus(x,y) = 8.794e-5 * exp(log(x)(1.81)) / y ! 8.794e-5x**1.81/y
	271	viscos(x,y) = 1.496e-6 * (xsqrt(x)) /(x+120.)/y ! 1.496e-6x**1.5/(x+120.)/y
	272	xka(x,y) = 1.414e3viscos(x,y)y
	273	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
	274	! venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
	275	! /viscos(b,c)*(.5)(den0/c)**0.25
	276	venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
	277	/sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
	278	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
	279	!
	280	pi = 4. * atan(1.)
	281	!
	282	!----------------------------------------------------------------
	283	! paddint 0 for negative values generated by dynamics
	284	!
	285	do k = kts, kte
	286	do i = its, ite
	287	qci(i,k) = max(qci(i,k),0.0)
	288	qrs(i,k) = max(qrs(i,k),0.0)
	289	enddo
	290	enddo
	291	!
	292	!----------------------------------------------------------------
	293	! latent heat for phase changes and heat capacity. neglect the
	294	! changes during microphysical process calculation
	295	! emanuel(1994)
	296	!
	297	do k = kts, kte
	298	do i = its, ite
	299	cpm(i,k) = cpmcal(q(i,k))
	300	xl(i,k) = xlcal(t(i,k))
	301	enddo
	302	enddo
	303	!
	304	!----------------------------------------------------------------
	305	! compute the minor time steps.
	306	!
	307	loops = max(nint(delt/dtcldcr),1)
	308	dtcld = delt/loops
	309	if(delt.le.dtcldcr) dtcld = delt
	310	!
	311	do loop = 1,loops
	312	!
	313	!----------------------------------------------------------------
	314	! initialize the large scale variables
	315	!
	316	do i = its, ite
	317	mstep(i) = 1
	318	flgcld(i) = .true.
	319	enddo
	320	!
	321	! do k = kts, kte
	322	! do i = its, ite
	323	! denfac(i,k) = sqrt(den0/den(i,k))
	324	! enddo
	325	! enddo
	326	do k = kts, kte
	327	CALL VREC( tvec1(its), den(its,k), ite-its+1)
	328	do i = its, ite
	329	tvec1(i) = tvec1(i)*den0
	330	enddo
	331	CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
	332	enddo
	333	!
	334	! Inline expansion for fpvs
	335	! qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	336	! qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
	337	cvap = cpv
	338	hvap=xlv0
	339	hsub=xls
	340	ttp=t0c+0.01
	341	dldt=cvap-cliq
	342	xa=-dldt/rv
	343	xb=xa+hvap/(rv*ttp)
	344	dldti=cvap-cice
	345	xai=-dldti/rv
	346	xbi=xai+hsub/(rv*ttp)
	347	do k = kts, kte
	348	do i = its, ite
	349	! tr=ttp/t(i,k)
	350	! if(t(i,k).lt.ttp) then
	351	! qs(i,k) =psat(trxai)exp(xbi*(1.-tr))
	352	! else
	353	! qs(i,k) =psat(trxa)exp(xb*(1.-tr))
	354	! endif
	355	! qs0(i,k) =psat(trxa)exp(xb*(1.-tr))
	356	tr=ttp/t(i,k)
	357	if(t(i,k).lt.ttp) then
	358	qs(i,k) =psat(exp(log(tr)(xai)))exp(xbi(1.-tr))
	359	else
	360	qs(i,k) =psat(exp(log(tr)(xa)))exp(xb(1.-tr))
	361	endif
	362	qs0(i,k) =psat(exp(log(tr)(xa)))exp(xb(1.-tr))
	363	qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
	364	qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
	365	qs(i,k) = max(qs(i,k),qmin)
	366	rh(i,k) = max(q(i,k) / qs(i,k),qmin)
	367	enddo
	368	enddo
	369	!
	370	!----------------------------------------------------------------
	371	! initialize the variables for microphysical physics
	372	!
	373	!
	374	do k = kts, kte
	375	do i = its, ite
	376	pres(i,k) = 0.
	377	paut(i,k) = 0.
	378	pacr(i,k) = 0.
	379	pgen(i,k) = 0.
	380	pisd(i,k) = 0.
	381	pcon(i,k) = 0.
	382	fall(i,k) = 0.
	383	falk(i,k) = 0.
	384	fallc(i,k) = 0.
	385	falkc(i,k) = 0.
	386	xni(i,k) = 1.e3
	387	enddo
	388	enddo
	389	!
	390	!----------------------------------------------------------------
	391	! compute the fallout term:
	392	! first, vertical terminal velosity for minor loops
	393	!---------------------------------------------------------------
	394	! n0s: Intercept parameter for snow [m-4] [HDC 6]
	395	!---------------------------------------------------------------
	396	do k = kts, kte
	397	do i = its, ite
	398	supcol = t0c-t(i,k)
	399	n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
	400	if(t(i,k).ge.t0c) then
	401	if(qrs(i,k).le.qcrmin)then
	402	rslope(i,k) = rslopermax
	403	rslopeb(i,k) = rsloperbmax
	404	rslope2(i,k) = rsloper2max
	405	rslope3(i,k) = rsloper3max
	406	else
	407	rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
	408	! rslopeb(i,k) = rslope(i,k)**bvtr
	409	rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
	410	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	411	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	412	endif
	413	else
	414	if(qrs(i,k).le.qcrmin)then
	415	rslope(i,k) = rslopesmax
	416	rslopeb(i,k) = rslopesbmax
	417	rslope2(i,k) = rslopes2max
	418	rslope3(i,k) = rslopes3max
	419	else
	420	rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
	421	! rslopeb(i,k) = rslope(i,k)**bvts
	422	rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
	423	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	424	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	425	endif
	426	endif
	427	!-------------------------------------------------------------
	428	! Ni: ice crystal number concentraiton [HDC 5c]
	429	!-------------------------------------------------------------
	430	! xni(i,k) = min(max(5.38e7*(den(i,k) &
	431	! max(qci(i,k),qmin))*0.75,1.e3),1.e6)
	432	xni(i,k) = min(max(5.38e7exp(log((den(i,k)max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
	433	enddo
	434	enddo
	435	!
	436	mstepmax = 1
	437	numdt = 1
	438	do k = kte, kts, -1
	439	do i = its, ite
	440	if(t(i,k).lt.t0c) then
	441	pvt = pvts
	442	else
	443	pvt = pvtr
	444	endif
	445	work1(i,k) = pvtrslopeb(i,k)denfac(i,k)
	446	work2(i,k) = work1(i,k)/delz(i,k)
	447	numdt(i) = max(nint(work2(i,k)*dtcld+.5),1)
	448	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
	449	enddo
	450	enddo
	451	do i = its, ite
	452	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
	453	enddo
	454	!
	455	do n = 1, mstepmax
	456	k = kte
	457	do i = its, ite
	458	if(n.le.mstep(i)) then
	459	falk(i,k) = den(i,k)qrs(i,k)work2(i,k)/mstep(i)
	460	hold = falk(i,k)
	461	fall(i,k) = fall(i,k)+falk(i,k)
	462	holdrs = qrs(i,k)
	463	qrs(i,k) = max(qrs(i,k)-falk(i,k)*dtcld/den(i,k),0.)
	464	endif
	465	enddo
	466	do k = kte-1, kts, -1
	467	do i = its, ite
	468	if(n.le.mstep(i)) then
	469	falk(i,k) = den(i,k)qrs(i,k)work2(i,k)/mstep(i)
	470	hold = falk(i,k)
	471	fall(i,k) = fall(i,k)+falk(i,k)
	472	holdrs = qrs(i,k)
	473	qrs(i,k) = max(qrs(i,k)-(falk(i,k) &
	474	-falk(i,k+1)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	475	endif
	476	enddo
	477	enddo
	478	enddo
	479	!---------------------------------------------------------------
	480	! Vice [ms-1] : fallout of ice crystal [HDC 5a]
	481	!---------------------------------------------------------------
	482	mstepmax = 1
	483	mstep = 1
	484	numdt = 1
	485	do k = kte, kts, -1
	486	do i = its, ite
	487	if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then
	488	xmi = den(i,k)*qci(i,k)/xni(i,k)
	489	! diameter = dicon * sqrt(xmi)
	490	! work1c(i,k) = 1.49e4diameter*1.31
	491	diameter = max(dicon * sqrt(xmi), 1.e-25)
	492	work1c(i,k) = 1.49e4exp(log(diameter)(1.31))
	493	else
	494	work1c(i,k) = 0.
	495	endif
	496	if(qci(i,k).le.0.) then
	497	work2c(i,k) = 0.
	498	else
	499	work2c(i,k) = work1c(i,k)/delz(i,k)
	500	endif
	501	numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
	502	if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
	503	enddo
	504	enddo
	505	do i = its, ite
	506	if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
	507	enddo
	508	!
	509	do n = 1, mstepmax
	510	k = kte
	511	do i = its, ite
	512	if (n.le.mstep(i)) then
	513	falkc(i,k) = den(i,k)qci(i,k)work2c(i,k)/mstep(i)
	514	holdc = falkc(i,k)
	515	fallc(i,k) = fallc(i,k)+falkc(i,k)
	516	holdci = qci(i,k)
	517	qci(i,k) = max(qci(i,k)-falkc(i,k)*dtcld/den(i,k),0.)
	518	endif
	519	enddo
	520	do k = kte-1, kts, -1
	521	do i = its, ite
	522	if (n.le.mstep(i)) then
	523	falkc(i,k) = den(i,k)qci(i,k)work2c(i,k)/mstep(i)
	524	holdc = falkc(i,k)
	525	fallc(i,k) = fallc(i,k)+falkc(i,k)
	526	holdci = qci(i,k)
	527	qci(i,k) = max(qci(i,k)-(falkc(i,k) &
	528	-falkc(i,k+1)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
	529	endif
	530	enddo
	531	enddo
	532	enddo
	533	!
	534	!----------------------------------------------------------------
	535	! compute the freezing/melting term. [D89 B16-B17]
	536	! freezing occurs one layer above the melting level
	537	!
	538	do i = its, ite
	539	mstep(i) = 0
	540	enddo
	541	do k = kts, kte
	542	do i = its, ite
	543	if(t(i,k).ge.t0c) then
	544	mstep(i) = k
	545	endif
	546	enddo
	547	enddo
	548	!
	549	do i = its, ite
	550	kwork2(i) = mstep(i)
	551	kwork1(i) = mstep(i)
	552	if(mstep(i).ne.0) then
	553	if (w(i,mstep(i)).gt.0.) then
	554	kwork1(i) = mstep(i) + 1
	555	endif
	556	endif
	557	enddo
	558	!
	559	do i = its, ite
	560	k = kwork1(i)
	561	kk = kwork2(i)
	562	if(k*kk.ge.1) then
	563	qrsci = qrs(i,k) + qci(i,k)
	564	if(qrsci.gt.0..or.fall(i,kk).gt.0.) then
	565	frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld), &
	566	qrsci/dtcld)
	567	snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld), &
	568	qrs(i,k)/dtcld)
	569	if(k.eq.kk) then
	570	t(i,k) = t(i,k) - xlf0/cpm(i,k)(frzmlt+snomlt)dtcld
	571	else
	572	t(i,k) = t(i,k) - xlf0/cpm(i,k)frzmltdtcld
	573	t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)snomltdtcld
	574	endif
	575	endif
	576	endif
	577	enddo
	578	!
	579	!----------------------------------------------------------------
	580	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
	581	!
	582	do i = its, ite
	583	fallsum = fall(i,1)
	584	fallsum_qsi = 0.
	585	if((t0c-t(i,1)).gt.0) then
	586	fallsum = fallsum+fallc(i,1)
	587	fallsum_qsi = fall(i,1)+fallc(i,1)
	588	endif
	589	rainncv(i) = 0.
	590	if(fallsum.gt.0.) then
	591	rainncv(i) = fallsumdelz(i,1)/denrdtcld*1000.
	592	rain(i) = fallsumdelz(i,1)/denrdtcld*1000. &
	593	+ rain(i)
	594	endif
	595	IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
	596	snowncv(i) = 0.
	597	if(fallsum_qsi.gt.0.) then
	598	snowncv(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000.
	599	snow(i) = fallsum_qsidelz(i,kts)/denrdtcld*1000. + snow(i)
	600	endif
	601	ENDIF
	602	sr(i) = 0.
	603	if(fallsum.gt.0.)sr(i)=fallsum_qsidelz(i,kts)/denrdtcld*1000./(rainncv(i)+1.e-12)
	604	enddo
	605	!
	606	!----------------------------------------------------------------
	607	! rsloper: reverse of the slope parameter of the rain(m)
	608	! xka: thermal conductivity of air(jm-1s-1k-1)
	609	! work1: the thermodynamic term in the denominator associated with
	610	! heat conduction and vapor diffusion
	611	! (ry88, y93, h85)
	612	! work2: parameter associated with the ventilation effects(y93)
	613	!
	614	do k = kts, kte
	615	do i = its, ite
	616	if(t(i,k).ge.t0c) then
	617	if(qrs(i,k).le.qcrmin)then
	618	rslope(i,k) = rslopermax
	619	rslopeb(i,k) = rsloperbmax
	620	rslope2(i,k) = rsloper2max
	621	rslope3(i,k) = rsloper3max
	622	else
	623	rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
	624	rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
	625	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	626	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	627	endif
	628	else
	629	if(qrs(i,k).le.qcrmin)then
	630	rslope(i,k) = rslopesmax
	631	rslopeb(i,k) = rslopesbmax
	632	rslope2(i,k) = rslopes2max
	633	rslope3(i,k) = rslopes3max
	634	else
	635	rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
	636	rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
	637	rslope2(i,k) = rslope(i,k)*rslope(i,k)
	638	rslope3(i,k) = rslope2(i,k)*rslope(i,k)
	639	endif
	640	endif
	641	enddo
	642	enddo
	643	!
	644	do k = kts, kte
	645	do i = its, ite
	646	if(t(i,k).ge.t0c) then
	647	work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k))
	648	else
	649	work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k))
	650	endif
	651	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
	652	enddo
	653	enddo
	654	!
	655	do k = kts, kte
	656	do i = its, ite
	657	supsat = max(q(i,k),qmin)-qs(i,k)
	658	satdt = supsat/dtcld
	659	if(t(i,k).ge.t0c) then
	660	!
	661	!===============================================================
	662	!
	663	! warm rain processes
	664	!
	665	! - follows the processes in RH83 and LFO except for autoconcersion
	666	!
	667	!===============================================================
	668	!---------------------------------------------------------------
	669	! praut: auto conversion rate from cloud to rain [HDC 16]
	670	! (C->R)
	671	!---------------------------------------------------------------
	672	if(qci(i,k).gt.qc0) then
	673	! paut(i,k) = qck1qci(i,k)*(7./3.)
	674	paut(i,k) = qck1exp(log(qci(i,k))((7./3.)))
	675	paut(i,k) = min(paut(i,k),qci(i,k)/dtcld)
	676	endif
	677	!---------------------------------------------------------------
	678	! pracw: accretion of cloud water by rain [HL A40] [D89 B15]
	679	! (C->R)
	680	!---------------------------------------------------------------
	681	if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
	682	pacr(i,k) = min(pacrrrslope3(i,k)rslopeb(i,k) &
	683	qci(i,k)denfac(i,k),qci(i,k)/dtcld)
	684	endif
	685	!---------------------------------------------------------------
	686	! prevp: evaporation/condensation rate of rain [HDC 14]
	687	! (V->R or R->V)
	688	!---------------------------------------------------------------
	689	if(qrs(i,k).gt.0.) then
	690	coeres = rslope2(i,k)sqrt(rslope(i,k)rslopeb(i,k))
	691	pres(i,k) = (rh(i,k)-1.)(precr1rslope2(i,k) &
	692	+precr2work2(i,k)coeres)/work1(i,k)
	693	if(pres(i,k).lt.0.) then
	694	pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
	695	pres(i,k) = max(pres(i,k),satdt/2)
	696	else
	697	pres(i,k) = min(pres(i,k),satdt/2)
	698	endif
	699	endif
	700	else
	701	!
	702	!===============================================================
	703	!
	704	! cold rain processes
	705	!
	706	! - follows the revised ice microphysics processes in HDC
	707	! - the processes same as in RH83 and LFO behave
	708	! following ice crystal hapits defined in HDC, inclduing
	709	! intercept parameter for snow (n0s), ice crystal number
	710	! concentration (ni), ice nuclei number concentration
	711	! (n0i), ice diameter (d)
	712	!
	713	!===============================================================
	714	!
	715	supcol = t0c-t(i,k)
	716	ifsat = 0
	717	!-------------------------------------------------------------
	718	! Ni: ice crystal number concentraiton [HDC 5c]
	719	!-------------------------------------------------------------
	720	! xni(i,k) = min(max(5.38e7*(den(i,k) &
	721	! max(qci(i,k),qmin))*0.75,1.e3),1.e6)
	722	xni(i,k) = min(max(5.38e7exp(log((den(i,k)max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
	723	eacrs = exp(0.07*(-supcol))
	724	!
	725	if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
	726	xmi = den(i,k)*qci(i,k)/xni(i,k)
	727	diameter = min(dicon * sqrt(xmi),dimax)
	728	vt2i = 1.49e4diameter*1.31
	729	vt2s = pvtsrslopeb(i,k)denfac(i,k)
	730	!-------------------------------------------------------------
	731	! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
	732	! (T<T0: I->R)
	733	!-------------------------------------------------------------
	734	acrfac = 2.rslope3(i,k)+2.diameter*rslope2(i,k) &
	735	+diameter*2rslope(i,k)
	736	pacr(i,k) = min(piqci(i,k)eacrsn0sn0sfac(i,k) &
	737	abs(vt2s-vt2i)acrfac/4.,qci(i,k)/dtcld)
	738	endif
	739	!-------------------------------------------------------------
	740	! pidep: Deposition/Sublimation rate of ice [HDC 9]
	741	! (T<T0: V->I or I->V)
	742	!-------------------------------------------------------------
	743	if(qci(i,k).gt.0.) then
	744	xmi = den(i,k)*qci(i,k)/xni(i,k)
	745	diameter = dicon * sqrt(xmi)
	746	pisd(i,k) = 4.diameterxni(i,k)*(rh(i,k)-1.)/work1(i,k)
	747	if(pisd(i,k).lt.0.) then
	748	pisd(i,k) = max(pisd(i,k),satdt/2)
	749	pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld)
	750	else
	751	pisd(i,k) = min(pisd(i,k),satdt/2)
	752	endif
	753	if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
	754	endif
	755	!-------------------------------------------------------------
	756	! psdep: deposition/sublimation rate of snow [HDC 14]
	757	! (V->S or S->V)
	758	!-------------------------------------------------------------
	759	if(qrs(i,k).gt.0..and.ifsat.ne.1) then
	760	coeres = rslope2(i,k)sqrt(rslope(i,k)rslopeb(i,k))
	761	pres(i,k) = (rh(i,k)-1.)n0sfac(i,k)(precs1*rslope2(i,k) &
	762	+precs2work2(i,k)coeres)/work1(i,k)
	763	supice = satdt-pisd(i,k)
	764	if(pres(i,k).lt.0.) then
	765	pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
	766	pres(i,k) = max(max(pres(i,k),satdt/2),supice)
	767	else
	768	pres(i,k) = min(min(pres(i,k),satdt/2),supice)
	769	endif
	770	if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
	771	endif
	772	!-------------------------------------------------------------
	773	! pigen: generation(nucleation) of ice from vapor [HDC 7-8]
	774	! (T<T0: V->I)
	775	!-------------------------------------------------------------
	776	if(supsat.gt.0.and.ifsat.ne.1) then
	777	supice = satdt-pisd(i,k)-pres(i,k)
	778	xni0 = 1.e3exp(0.1supcol)
	779	! roqi0 = 4.92e-11xni0*1.33
	780	roqi0 = 4.92e-11exp(log(xni0)(1.33))
	781	pgen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k),0.))/dtcld)
	782	pgen(i,k) = min(min(pgen(i,k),satdt),supice)
	783	endif
	784	!-------------------------------------------------------------
	785	! psaut: conversion(aggregation) of ice to snow [HDC 12]
	786	! (T<T0: I->S)
	787	!-------------------------------------------------------------
	788	if(qci(i,k).gt.0.) then
	789	qimax = roqimax/den(i,k)
	790	paut(i,k) = max(0.,(qci(i,k)-qimax)/dtcld)
	791	endif
	792	endif
	793	enddo
	794	enddo
	795	!
	796	!----------------------------------------------------------------
	797	! check mass conservation of generation terms and feedback to the
	798	! large scale
	799	!
	800	do k = kts, kte
	801	do i = its, ite
	802	qciik = max(qmin,qci(i,k))
	803	delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
	804	if(delqci.ge.qciik) then
	805	facqci = qciik/delqci
	806	paut(i,k) = paut(i,k)*facqci
	807	pacr(i,k) = pacr(i,k)*facqci
	808	pgen(i,k) = pgen(i,k)*facqci
	809	pisd(i,k) = pisd(i,k)*facqci
	810	endif
	811	qik = max(qmin,q(i,k))
	812	delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
	813	if(delq.ge.qik) then
	814	facq = qik/delq
	815	pres(i,k) = pres(i,k)*facq
	816	pgen(i,k) = pgen(i,k)*facq
	817	pisd(i,k) = pisd(i,k)*facq
	818	endif
	819	work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
	820	q(i,k) = q(i,k)+work2(i,k)*dtcld
	821	qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k) &
	822	-pisd(i,k))*dtcld,0.)
	823	qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k) &
	824	+pres(i,k))*dtcld,0.)
	825	if(t(i,k).lt.t0c) then
	826	t(i,k) = t(i,k)-xlswork2(i,k)/cpm(i,k)dtcld
	827	else
	828	t(i,k) = t(i,k)-xl(i,k)work2(i,k)/cpm(i,k)dtcld
	829	endif
	830	enddo
	831	enddo
	832	!
	833	cvap = cpv
	834	hvap = xlv0
	835	hsub = xls
	836	ttp=t0c+0.01
	837	dldt=cvap-cliq
	838	xa=-dldt/rv
	839	xb=xa+hvap/(rv*ttp)
	840	dldti=cvap-cice
	841	xai=-dldti/rv
	842	xbi=xai+hsub/(rv*ttp)
	843	do k = kts, kte
	844	do i = its, ite
	845	tr=ttp/t(i,k)
	846	! qs(i,k)=psat(trxa)exp(xb*(1.-tr))
	847	qs(i,k)=psat(exp(log(tr)(xa)))exp(xb(1.-tr))
	848	qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
	849	qs(i,k) = max(qs(i,k),qmin)
	850	denfac(i,k) = sqrt(den0/den(i,k))
	851	enddo
	852	enddo
	853	!
	854	!----------------------------------------------------------------
	855	! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
	856	! if there exists additional water vapor condensated/if
	857	! evaporation of cloud water is not enough to remove subsaturation
	858	!
	859	do k = kts, kte
	860	do i = its, ite
	861	work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k))
	862	work2(i,k) = qci(i,k)+work1(i,k)
	863	pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld
	864	if(qci(i,k).gt.0..and.work1(i,k).lt.0.and.t(i,k).gt.t0c) &
	865	pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld
	866	q(i,k) = q(i,k)-pcon(i,k)*dtcld
	867	qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.)
	868	t(i,k) = t(i,k)+pcon(i,k)xl(i,k)/cpm(i,k)dtcld
	869	enddo
	870	enddo
	871	!
	872	!----------------------------------------------------------------
	873	! padding for small values
	874	!
	875	do k = kts, kte
	876	do i = its, ite
	877	if(qci(i,k).le.qmin) qci(i,k) = 0.0
	878	enddo
	879	enddo
	880	!
	881	enddo ! big loops
	882	END SUBROUTINE wsm32D
	883	! ...................................................................
	884	REAL FUNCTION rgmma(x)
	885	!-------------------------------------------------------------------
	886	IMPLICIT NONE
	887	!-------------------------------------------------------------------
	888	! rgmma function: use infinite product form
	889	REAL :: euler
	890	PARAMETER (euler=0.577215664901532)
	891	REAL :: x, y
	892	INTEGER :: i
	893	if(x.eq.1.)then
	894	rgmma=0.
	895	else
	896	rgmma=xexp(eulerx)
	897	do i=1,10000
	898	y=float(i)
	899	rgmma=rgmma(1.000+x/y)exp(-x/y)
	900	enddo
	901	rgmma=1./rgmma
	902	endif
	903	END FUNCTION rgmma
	904	!
	905	!--------------------------------------------------------------------------
	906	REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
	907	!--------------------------------------------------------------------------
	908	IMPLICIT NONE
	909	!--------------------------------------------------------------------------
	910	REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
	911	xai,xbi,ttp,tr
	912	INTEGER ice
	913	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	914	ttp=t0c+0.01
	915	dldt=cvap-cliq
	916	xa=-dldt/rv
	917	xb=xa+hvap/(rv*ttp)
	918	dldti=cvap-cice
	919	xai=-dldti/rv
	920	xbi=xai+hsub/(rv*ttp)
	921	tr=ttp/t
	922	if(t.lt.ttp.and.ice.eq.1) then
	923	fpvs=psat(trxai)exp(xbi*(1.-tr))
	924	else
	925	fpvs=psat(trxa)exp(xb*(1.-tr))
	926	endif
	927	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	928	END FUNCTION fpvs
	929	!-------------------------------------------------------------------
	930	SUBROUTINE wsm3init(den0,denr,dens,cl,cpv,allowed_to_read)
	931	!-------------------------------------------------------------------
	932	IMPLICIT NONE
	933	!-------------------------------------------------------------------
	934	!.... constants which may not be tunable
	935	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
	936	LOGICAL, INTENT(IN) :: allowed_to_read
	937	REAL :: pi
	938	!
	939	pi = 4.*atan(1.)
	940	xlv1 = cl-cpv
	941	!
	942	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
	943	qck1 = .1049.8peaut/(xncrdenr)(1./3.)/xmyuden0**(4./3.) ! 7.03
	944	!
	945	bvtr1 = 1.+bvtr
	946	bvtr2 = 2.5+.5*bvtr
	947	bvtr3 = 3.+bvtr
	948	bvtr4 = 4.+bvtr
	949	g1pbr = rgmma(bvtr1)
	950	g3pbr = rgmma(bvtr3)
	951	g4pbr = rgmma(bvtr4) ! 17.837825
	952	g5pbro2 = rgmma(bvtr2) ! 1.8273
	953	pvtr = avtr*g4pbr/6.
	954	eacrr = 1.0
	955	pacrr = pin0ravtrg3pbr.25*eacrr
	956	precr1 = 2.pin0r*.78
	957	precr2 = 2.pin0r.31avtr*.5g5pbro2
	958	xm0 = (di0/dicon)**2
	959	xmmax = (dimax/dicon)**2
	960	roqimax = 2.08e22dimax*8
	961	!
	962	bvts1 = 1.+bvts
	963	bvts2 = 2.5+.5*bvts
	964	bvts3 = 3.+bvts
	965	bvts4 = 4.+bvts
	966	g1pbs = rgmma(bvts1) !.8875
	967	g3pbs = rgmma(bvts3)
	968	g4pbs = rgmma(bvts4) ! 12.0786
	969	g5pbso2 = rgmma(bvts2)
	970	pvts = avts*g4pbs/6.
	971	pacrs = pin0savtsg3pbs.25
	972	precs1 = 4.n0s.65
	973	precs2 = 4.n0s.44avts.5g5pbso2
	974	pidn0r = pidenrn0r
	975	pidn0s = pidensn0s
	976	!
	977	rslopermax = 1./lamdarmax
	978	rslopesmax = 1./lamdasmax
	979	rsloperbmax = rslopermax ** bvtr
	980	rslopesbmax = rslopesmax ** bvts
	981	rsloper2max = rslopermax * rslopermax
	982	rslopes2max = rslopesmax * rslopesmax
	983	rsloper3max = rsloper2max * rslopermax
	984	rslopes3max = rslopes2max * rslopesmax
	985	!
	986	END SUBROUTINE wsm3init
	987	END MODULE module_mp_wsm3

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

Download in other formats:

Original Format