Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_mp_ncloud3.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
Property svn:executable set to ``*
File size: 31.2 KB

Line
1	MODULE module_mp_ncloud3
2
3	REAL, PARAMETER, PRIVATE :: dtcldcr = 240.
4	INTEGER, PARAMETER, PRIVATE :: mstepmax = 100
5
6	REAL, PARAMETER, PRIVATE :: n0r = 8.e6
7	REAL, PARAMETER, PRIVATE :: avtr = 841.9
8	REAL, PARAMETER, PRIVATE :: bvtr = 0.8
9	REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
10	REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
11	REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
12	REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
13
14	REAL, PARAMETER, PRIVATE :: avts = 16.2
15	REAL, PARAMETER, PRIVATE :: bvts = .527
16	REAL, PARAMETER, PRIVATE :: xncmax = 1.e8
17	REAL, PARAMETER, PRIVATE :: n0smax = 1.e9
18	REAL, PARAMETER, PRIVATE :: betai = .6
19	REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
20	REAL, PARAMETER, PRIVATE :: dicon = 16.3
21	REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6*.8
22	REAL, PARAMETER, PRIVATE :: dimax = 400.e-6
23	REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
24	REAL, PARAMETER, PRIVATE :: alpha = 1./8.18 ! .122 exponen factor for n0s
25	! REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e15
26	REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e5
27	REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-6
28
29	REAL, SAVE :: &
30	qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
31	g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
32	precr1,precr2,xm0,xmmax,bvts1, &
33	bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
34	g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
35	pidn0s,xlv1
36
37	CONTAINS
38
39	!===================================================================
40	!
41	SUBROUTINE ncloud3(th, q, qci, qrs &
42	, w, den, pii, p, delz &
43	, delt,g, cpd, cpv, rd, rv, t0c &
44	, ep1, ep2, qmin &
45	, XLS, XLV0, XLF0, den0, denr &
46	, cliq,cice,psat &
47	, rain, rainncv &
48	, ids,ide, jds,jde, kds,kde &
49	, ims,ime, jms,jme, kms,kme &
50	, its,ite, jts,jte, kts,kte &
51	)
52
53	!-------------------------------------------------------------------
54	IMPLICIT NONE
55	!-------------------------------------------------------------------
56	!
57	!Coded by Song-You Hong (NCEP) and implemented by Shuhua Chen (NCAR)
58	!
59	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
60	ims,ime, jms,jme, kms,kme , &
61	its,ite, jts,jte, kts,kte
62
63	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
64	INTENT(INOUT) :: &
65	th, &
66	q, &
67	qci, &
68	qrs
69
70	REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
71	INTENT(IN ) :: w, &
72	den, &
73	pii, &
74	p, &
75	delz
76
77	REAL, INTENT(IN ) :: delt, &
78	g, &
79	rd, &
80	rv, &
81	T0c, &
82	den0, &
83	cpd, &
84	cpv, &
85	ep1, &
86	ep2, &
87	qmin, &
88	XLS, &
89	XLV0, &
90	XLF0, &
91	cliq, &
92	cice, &
93	psat, &
94	denr
95	REAL, DIMENSION( ims:ime , jms:jme ), &
96	INTENT(INOUT) :: rain, &
97	rainncv
98
99	! LOCAL VAR
100
101	REAL, DIMENSION( its:ite , kts:kte ) :: t
102	INTEGER :: i,j,k
103
104	!-------------------------------------------------------------------
105	DO J=jts,jte
106
107	DO K=kts,kte
108	DO I=its,ite
109	t(i,k)=th(i,k,j)*pii(i,k,j)
110	ENDDO
111	ENDDO
112
113	CALL ncloud32D(t, q(ims,kms,j), qci(ims,kms,j) &
114	,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j) &
115	,p(ims,kms,j), delz(ims,kms,j) &
116	, delt,g, cpd, cpv, rd, rv, t0c &
117	,ep1, ep2, qmin &
118	,XLS, XLV0, XLF0, den0, denr &
119	,cliq,cice,psat &
120	,J &
121	,rain(ims,j), rainncv(ims,j) &
122	,ids,ide, jds,jde, kds,kde &
123	,ims,ime, jms,jme, kms,kme &
124	,its,ite, jts,jte, kts,kte &
125	)
126
127	DO K=kts,kte
128	DO I=its,ite
129	th(i,k,j)=t(i,k)/pii(i,k,j)
130	ENDDO
131	ENDDO
132
133	ENDDO
134
135	END SUBROUTINE ncloud3
136
137	!===================================================================
138	!
139	SUBROUTINE ncloud32D(t, q, qci, qrs,w, den, p, delz &
140	,delt,g, cpd, cpv, rd, rv, t0c &
141	,ep1, ep2, qmin &
142	,XLS, XLV0, XLF0, den0, denr &
143	,cliq,cice,psat &
144	,lat &
145	,rain, rainncv &
146	,ids,ide, jds,jde, kds,kde &
147	,ims,ime, jms,jme, kms,kme &
148	,its,ite, jts,jte, kts,kte &
149	)
150	!-------------------------------------------------------------------
151	IMPLICIT NONE
152	!-------------------------------------------------------------------
153	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
154	ims,ime, jms,jme, kms,kme , &
155	its,ite, jts,jte, kts,kte, &
156	lat
157
158	REAL, DIMENSION( its:ite , kts:kte ), &
159	INTENT(INOUT) :: &
160	t
161
162	REAL, DIMENSION( ims:ime , kms:kme ), &
163	INTENT(INOUT) :: &
164	q, &
165	qci, &
166	qrs
167
168	REAL, DIMENSION( ims:ime , kms:kme ), &
169	INTENT(IN ) :: w, &
170	den, &
171	p, &
172	delz
173
174	REAL, INTENT(IN ) :: delt, &
175	g, &
176	cpd, &
177	cpv, &
178	T0c, &
179	den0, &
180	rd, &
181	rv, &
182	ep1, &
183	ep2, &
184	qmin, &
185	XLS, &
186	XLV0, &
187	XLF0, &
188	cliq, &
189	cice, &
190	psat, &
191	denr
192	REAL, DIMENSION( ims:ime ), &
193	INTENT(INOUT) :: rain, &
194	rainncv
195
196	! LOCAL VAR
197
198	INTEGER, PARAMETER :: iun = 84
199
200	REAL, DIMENSION( its:ite , kts:kte ) :: &
201	rh, qs, denfac, slope, slope2, slopeb, &
202	pgen, paut, pacr, pisd, pres, pcon, fall, falk, &
203	xl, cpm, work1, work2, q1, t1, &
204	pgens, pauts, pacrss, pisds, press, pcons
205
206	REAL, DIMENSION( its:ite , kts:kte ) :: &
207	falkc, work1c, work2c, fallc
208
209	INTEGER, DIMENSION( its:ite ) :: mstep
210	LOGICAL, DIMENSION( its:ite ) :: flgcld
211
212	REAL :: n0sfac, pi, &
213	cpmcal, xlcal, tvcal, lamdar, lamdas, diffus, &
214	viscos, xka, venfac, conden, diffac, &
215	x, y, z, a, b, c, d, e, &
216	qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
217	snomlt, hold, holdrs, facqci, supcol, coeres, &
218	supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, xnc
219
220	INTEGER :: i,j,k, &
221	iprt, latd, lond, loop, loops, ifsat, kk, n, numdt
222
223	!
224	!=================================================================
225	! compute internal functions
226	!
227	cpmcal(x) = cpd(1.-max(x,qmin))+max(x,qmin)cpv
228	xlcal(x) = xlv0-xlv1*(x-t0c)
229	tvcal(x,y) = x+xep1max(y,qmin)
230	!----------------------------------------------------------------
231	! size distributions: (x=mixing ratio, y=air density):
232	! valid for mixing ratio > 1.e-30 kg/kg.
233	! otherwise use uniform distribution value (1.e15)
234	!
235	lamdar(x,y)=(pidn0r/(xy))*.25
236	lamdas(x,y,z)=(pidn0sz/(xy))**.25
237	!
238	!----------------------------------------------------------------
239	! diffus: diffusion coefficient of the water vapor
240	! viscos: kinematic viscosity(m2s-1)
241	!
242	diffus(x,y) = 8.794e-5x*1.81/y
243	viscos(x,y) = 1.496e-6x*1.5/(x+120.)/y
244	xka(x,y) = 1.414e3viscos(x,y)y
245	diffac(a,b,c,d,e) = daa/(xka(c,d)rvcc)+1./(ediffus(c,b))
246	venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
247	/viscos(b,c)*(.5)(den0/c)**0.25
248	conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+dd/(rve)c/(aa))
249	!
250	pi = 4. * atan(1.)
251	!
252	!=================================================================
253	! set iprt = 0 for no unit fort.84 output
254	!
255	! iprt = 0
256	! if(iprt.eq.1) then
257	! qdt = delt * 1000.
258	! latd = jts
259	! lond = its
260	! else
261	! latd = 0
262	! lond = 0
263	! endif
264	!
265	!----------------------------------------------------------------
266	! paddint 0 for negative values generated by dynamics
267	!
268	do k = kts, kte
269	do i = its, ite
270	qci(i,k) = max(qci(i,k),0.0)
271	qrs(i,k) = max(qrs(i,k),0.0)
272	enddo
273	enddo
274	!
275	!----------------------------------------------------------------
276	! latent heat for phase changes and heat capacity. neglect the
277	! changes during microphysical process calculation
278	! emanuel(1994)
279	!
280	do k = kts, kte
281	do i = its, ite
282	cpm(i,k) = cpmcal(q(i,k))
283	xl(i,k) = xlcal(t(i,k))
284	enddo
285	enddo
286	!
287	!----------------------------------------------------------------
288	! compute the minor time steps.
289	!
290	loops = max(nint(delt/dtcldcr),1)
291	dtcld = delt/loops
292	if(delt.le.dtcldcr) dtcld = delt
293	!
294	do loop = 1,loops
295	!
296	!----------------------------------------------------------------
297	! initialize the large scale variables
298	!
299	do i = its, ite
300	mstep(i) = 1
301	flgcld(i) = .true.
302	enddo
303	!
304	do k = kts, kte
305	do i = its, ite
306	work1(i,k) = tvcal(t(i,k),q(i,k))
307	denfac(i,k) = sqrt(den0/den(i,k))
308	enddo
309	enddo
310	!
311	do k = kts, kte
312	do i = its, ite
313	qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
314	qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
315	qs(i,k) = max(qs(i,k),qmin)
316	rh(i,k) = max(q(i,k) / qs(i,k),qmin)
317	enddo
318	enddo
319	!
320	!----------------------------------------------------------------
321	! initialize the variables for microphysical physics
322	!
323	! if(lat.eq.latd) then
324	! i = lond
325	! print*,'lat',latd,lat,i
326
327	! do k = kts, kte
328	! press(i,k) = 0.
329	! pauts(i,k) = 0.
330	! pacrss(i,k)= 0.
331	! pgens(i,k) = 0.
332	! pisds(i,k) = 0.
333	! pcons(i,k) = 0.
334	! t1(i,k) = t(i,k)
335	! q1(i,k) = q(i,k)
336	! enddo
337	! endif
338	!
339	do k = kts, kte
340	do i = its, ite
341	pres(i,k) = 0.
342	paut(i,k) = 0.
343	pacr(i,k) = 0.
344	pgen(i,k) = 0.
345	pisd(i,k) = 0.
346	pcon(i,k) = 0.
347	fall(i,k) = 0.
348	falk(i,k) = 0.
349	fallc(i,k) = 0.
350	falkc(i,k) = 0.
351	enddo
352	enddo
353	!
354	!----------------------------------------------------------------
355	! sloper: the slope parameter of the rain(m-1)
356	! xka: thermal conductivity of air(jm-1s-1k-1)
357	! work1: the thermodynamic term in the denominator associated with
358	! heat conduction and vapor diffusion
359	! (ry88, y93, h85)
360	! work2: parameter associated with the ventilation effects(y93)
361	!
362	do k = kts, kte
363	do i = its, ite
364	if(qrs(i,k).le.qcrmin)then
365	slope(i,k) = lamdarmax
366	slopeb(i,k) = slope(i,k)**bvtr
367	else
368	if(t(i,k).ge.t0c) then
369	slope(i,k) = lamdar(qrs(i,k),den(i,k))
370	slopeb(i,k) = slope(i,k)**bvtr
371	else
372	supcol = t0c-t(i,k)
373	n0sfac = min(exp(alpha*supcol),n0smax)
374	slope(i,k) = lamdas(qrs(i,k),den(i,k),n0sfac)
375	slopeb(i,k) = slope(i,k)**bvts
376	endif
377	endif
378	slope2(i,k) = slope(i,k)*slope(i,k)
379	enddo
380	enddo
381	!
382	do k = kts, kte
383	do i = its, ite
384	if(t(i,k).ge.t0c) then
385	work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k))
386	else
387	work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k))
388	endif
389	work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
390	enddo
391	enddo
392	!
393	do k = kts, kte
394	do i = its, ite
395	supsat = max(q(i,k),qmin)-qs(i,k)
396	satdt = supsat/dtcld
397	if(t(i,k).ge.t0c) then
398	!
399	!----------------------------------------------------------------
400	! warm rain process
401	! paut: auto conversion rate from cloud to rain (kgkg-1s-1)(kessler)
402	! pacr: accretion rate of rain by cloud(lin83)
403	! pres: evaporation/condensation rate of rain(rh83)
404	!
405	if(qci(i,k).gt.qc0) then
406	paut(i,k) = qck1qci(i,k)*(7./3.)
407	paut(i,k) = min(paut(i,k),qci(i,k)/dtcld)
408	endif
409	!
410	if(qrs(i,k).gt.qcrmin) then
411	if(qci(i,k).gt.qcrmin) &
412	pacr(i,k) = min(pacrr/slope2(i,k)/slope(i,k)/slopeb(i,k) &
413	qci(i,k)denfac(i,k),qci(i,k)/dtcld)
414	coeres = slope2(i,k)sqrt(slope(i,k)slopeb(i,k))
415	pres(i,k) = (rh(i,k)-1.)*(precr1/slope2(i,k) &
416	+precr2*work2(i,k)/coeres)/work1(i,k)
417	if(pres(i,k).lt.0.) then
418	pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
419	pres(i,k) = max(pres(i,k),satdt/2)
420	else
421	pres(i,k) = min(pres(i,k),qrs(i,k)/dtcld)
422	pres(i,k) = min(pres(i,k),satdt/2)
423	endif
424	endif
425	else
426	!
427	!----------------------------------------------------------------
428	! cold rain process
429	! paut: conversion(aggregation) of ice to snow(kgkg-1s-1)(rh83)
430	! pgen: generation(nucleation) of ice from vapor(kgkg-1s-1)(rh83)
431	! pacr: accretion rate of snow by ice(lin83)
432	! pisd: deposition/sublimation rate of ice(rh83)
433	! pres: deposition/sublimation rate of snow(lin83)
434	!
435	supcol = t0c-t(i,k)
436	ifsat = 0
437	n0sfac = min(exp(alpha*supcol),n0smax)
438	xnc = min(xn0 * exp(betai*supcol)/den(i,k),xncmax)
439	!
440	if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qcrmin) then
441	eacrs = exp(0.025*(-supcol))
442	pacr(i,k) = pacrsn0sfaceacrs/slope2(i,k)/slope(i,k) &
443	/slopeb(i,k)qci(i,k)denfac(i,k)
444	endif
445	!
446	if(qci(i,k).gt.qcrmin) then
447	xmi = qci(i,k)*xnc
448	pisd(i,k) = 4.diconsqrt(xmi)den(i,k)(rh(i,k)-1.) &
449	/work1(i,k)
450	if(pisd(i,k).lt.0.) then
451	pisd(i,k) = max(pisd(i,k),satdt/2)
452	pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld)
453	else
454	pisd(i,k) = min(pisd(i,k),satdt/2)
455	endif
456	if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
457	endif
458	!
459	if(qrs(i,k).gt.qcrmin.and.ifsat.ne.1) then
460	coeres = slope2(i,k)sqrt(slope(i,k)slopeb(i,k))
461	pres(i,k) = (rh(i,k)-1.)n0sfac(precs1/slope2(i,k) &
462	+precs2*work2(i,k)/coeres)/work1(i,k)
463	if(pres(i,k).lt.0.) then
464	pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
465	pres(i,k) = max(pres(i,k),satdt/2)
466	else
467	pres(i,k) = min(pres(i,k),satdt/2)
468	pres(i,k) = min(pres(i,k),qrs(i,k)/dtcld)
469	endif
470	if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
471	endif
472	!
473	if(supsat.gt.0.and.ifsat.ne.1) then
474	pgen(i,k) = max(0.,(xm0*xnc-max(qci(i,k),0.))/dtcld)
475	pgen(i,k) = min(pgen(i,k),satdt)
476	endif
477	!
478	if(qci(i,k).gt.qcrmin) paut(i,k) &
479	= max(0.,(qci(i,k)-xmmax*xnc)/dtcld)
480	endif
481	enddo
482	enddo
483	!
484	!----------------------------------------------------------------
485	! check mass conservation of generation terms and feedback to the
486	! large scale
487	!
488	do k = kts, kte
489	do i = its, ite
490	qciik = max(qcrmin,qci(i,k))
491	delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
492	if(delqci.ge.qciik) then
493	facqci = qciik/delqci
494	paut(i,k) = paut(i,k)*facqci
495	pacr(i,k) = pacr(i,k)*facqci
496	pgen(i,k) = pgen(i,k)*facqci
497	pisd(i,k) = pisd(i,k)*facqci
498	endif
499	qik = max(qcrmin,q(i,k))
500	delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
501	if(delq.ge.qik) then
502	facq = qik/delq
503	pres(i,k) = pres(i,k)*facq
504	pgen(i,k) = pgen(i,k)*facq
505	pisd(i,k) = pisd(i,k)*facq
506	endif
507	work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
508	q(i,k) = q(i,k)+work2(i,k)*dtcld
509	qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k) &
510	-pisd(i,k))*dtcld,0.)
511	qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k) &
512	+pres(i,k))*dtcld,0.)
513	if(t(i,k).lt.t0c) then
514	t(i,k) = t(i,k)-xlswork2(i,k)/cpm(i,k)dtcld
515	else
516	t(i,k) = t(i,k)-xl(i,k)work2(i,k)/cpm(i,k)dtcld
517	endif
518	enddo
519	enddo
520	!
521	do k = kts, kte
522	do i = its, ite
523	qs(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
524	qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
525	qs(i,k) = max(qs(i,k),qmin)
526	denfac(i,k) = sqrt(den0/den(i,k))
527	enddo
528	enddo
529	!
530	!----------------------------------------------------------------
531	! condensational/evaporational rate of cloud water if there exists
532	! additional water vapor condensated/if evaporation of cloud water
533	! is not enough to remove subsaturation.
534	! use fall bariable for this process(pcon)
535	!
536	! if(lat.eq.latd) write(iun,603)
537	do k = kts, kte
538	do i = its, ite
539	work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k))
540	work2(i,k) = qci(i,k)+work1(i,k)
541	pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld
542	if(qci(i,k).gt.qcrmin.and.work1(i,k).lt.0.and.t(i,k).gt.t0c) &
543	pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld
544	q(i,k) = q(i,k)-pcon(i,k)*dtcld
545	qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.)
546	t(i,k) = t(i,k)+pcon(i,k)xl(i,k)/cpm(i,k)dtcld
547
548	!
549	! if(lat.eq.latd.and.i.eq.lond) then
550	! pgens(i,k) = pgens(i,k)+pgen(i,k)
551	! pcons(i,k) = pcons(i,k)+pcon(i,k)
552	! pisds(i,k) = pisds(i,k)+pisd(i,k)
553	! pacrss(i,k) = pacrss(i,k)+pacr(i,k)
554	! press(i,k) = press(i,k)+pres(i,k)
555	! pauts(i,k) = pauts(i,k)+paut(i,k)
556	! write(iun,604) k,p(i,k)/100., &
557	! t(i,k)-t0c,t(i,k)-t1(i,k),q(i,k)*1000., &
558	! (q(i,k)-q1(i,k))1000.,rh(i,k)100.,pgens(i,k)*qdt, &
559	! pcons(i,k)qdt,pisds(i,k)qdt,pauts(i,k)qdt,pacrss(i,k)qdt, &
560	! press(i,k)qdt,qci(i,k)1000.,qrs(i,k)*1000.
561	! endif
562
563	enddo
564	enddo
565	603 format(1x,' k',' p', &
566	' t',' delt',' q',' delq',' rh', &
567	' pgen',' pcon',' pisd',' paut',' pacr',' pres', &
568	' qci',' qrs')
569	604 format(1x,i3,f6.0,4f5.1,f5.0,8f5.2)
570	!
571	!----------------------------------------------------------------
572	! compute the fallout term:
573	! first, vertical terminal velosity for minor loops
574	!
575	do k = kts, kte
576	do i = its, ite
577	denfac(i,k) = sqrt(den0/den(i,k))
578	enddo
579	enddo
580	!
581	do k = kts, kte
582	do i = its, ite
583	if(qrs(i,k).le.qcrmin)then
584	slope(i,k) = lamdarmax
585	slopeb(i,k) = slope(i,k)**bvtr
586	else
587	if(t(i,k).ge.t0c) then
588	slope(i,k) = lamdar(qrs(i,k),den(i,k))
589	slopeb(i,k) = slope(i,k)**bvtr
590	else
591	supcol = t0c-t(i,k)
592	n0sfac = min(exp(alpha*supcol),n0smax)
593	slope(i,k) = lamdas(qrs(i,k),den(i,k),n0sfac)
594	slopeb(i,k) = slope(i,k)**bvts
595	endif
596	endif
597	slope2(i,k) = slope(i,k)*slope(i,k)
598	enddo
599	enddo
600	!
601	do i = its, ite
602	do k = kte, kts, -1
603	if(t(i,k).lt.t0c) then
604	pvt = pvts
605	else
606	pvt = pvtr
607	endif
608	work1(i,k) = pvt/slopeb(i,k)*denfac(i,k)
609	work2(i,k) = work1(i,k)/delz(i,k)
610	if(qrs(i,k).le.qcrmin) work2(i,k) = 0.
611	numdt = max(nint(work2(i,k)*dtcld+.5),1)
612	if(t(i,k).lt.t0c.and.qci(i,k).gt.qmin) then
613	work1c(i,k) = 3.29(den(i,k)qci(i,k))**0.16
614	else
615	work1c(i,k) = 0.
616	endif
617	if(qci(i,k).le.qmin) then
618	work2c(i,k) = 0.
619	else
620	work2c(i,k) = work1c(i,k)/delz(i,k)
621	endif
622	numdt = max(nint(work2c(i,k)*dtcld+.5),numdt)
623	if(numdt.ge.mstep(i)) mstep(i) = numdt
624	enddo
625	mstep(i) = min(mstep(i),mstepmax)
626	enddo
627	!
628	! if(lat.eq.latd) write(iun,605)
629	do n = 1,mstepmax
630	k = kte
631	do i = its, ite
632	if(n.le.mstep(i)) then
633	falk(i,k) = den(i,k)qrs(i,k)work2(i,k)/mstep(i)
634	hold = falk(i,k)
635	fall(i,k) = fall(i,k)+falk(i,k)
636	holdrs = qrs(i,k)
637	qrs(i,k) = max(qrs(i,k)-falk(i,k)*dtcld/den(i,k),0.)
638	falkc(i,k) = den(i,k)qci(i,k)work2c(i,k)/mstep(i)
639	fallc(i,k) = fallc(i,k)+falkc(i,k)
640	qci(i,k) = max(qci(i,k)-falkc(i,k)*dtcld/den(i,k),0.)
641	endif
642	enddo
643	do k = kte-1, kts, -1
644	do i = its, ite
645	if(n.le.mstep(i)) then
646	falk(i,k) = den(i,k)qrs(i,k)work2(i,k)/mstep(i)
647	hold = falk(i,k)
648	fall(i,k) = fall(i,k)+falk(i,k)
649	holdrs = qrs(i,k)
650	qrs(i,k) = max(qrs(i,k)-(falk(i,k) &
651	-falk(i,k+1)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
652	falkc(i,k) = den(i,k)qci(i,k)work2c(i,k)/mstep(i)
653	fallc(i,k) = fallc(i,k)+falkc(i,k)
654	qci(i,k) = max(qci(i,k)-(falkc(i,k) &
655	-falkc(i,k+1)delz(i,k+1)/delz(i,k))dtcld/den(i,k),0.)
656	endif
657	enddo
658	enddo
659	enddo
660	605 format(1x,' k',' p',' t',' q',' rh',' w', &
661	' vt',' falk',' falt',' qrsi',' qrsf',' mstep')
662	606 format(1x,i3,f6.0,2f5.1,f5.0,f6.2,5f6.2,i5)
663	!
664	!----------------------------------------------------------------
665	! compute the freezing/melting term.
666	! freezing occurs one layer above the melting level
667	!
668	do i = its, ite
669	mstep(i) = 0
670	enddo
671	do k = kts, kte
672	!
673	do i = its, ite
674	if(t(i,k).ge.t0c) then
675	mstep(i) = k
676	endif
677	enddo
678	enddo
679	!
680	do i = its, ite
681	if(mstep(i).ne.0.and.w(i,mstep(i)).gt.0.) then
682	work1(i,1) = float(mstep(i) + 1)
683	work1(i,2) = float(mstep(i))
684	else
685	work1(i,1) = float(mstep(i))
686	work1(i,2) = float(mstep(i))
687	endif
688	enddo
689	!
690	do i = its, ite
691	k = nint(work1(i,1))
692	kk = nint(work1(i,2))
693	if(k*kk.ge.1) then
694	qrsci = qrs(i,k) + qci(i,k)
695	if(qrsci.gt.qcrmin.or.fall(i,kk).gt.0.) then
696	frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld), &
697	qrsci/dtcld)
698	snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld),qrs(i,k)/dtcld)
699	if(k.eq.kk) then
700	t(i,k) = t(i,k) - xlf0/cpm(i,k)(frzmlt+snomlt)dtcld
701	else
702	t(i,k) = t(i,k) - xlf0/cpm(i,k)frzmltdtcld
703	t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)snomltdtcld
704	endif
705
706	! if(lat.eq.latd.and.i.eq.lond) write(iun,608) k,t(i,k)-t0c, &
707	! w(i,k),frzmltqdt,snomltqdt
708
709	endif
710	endif
711	enddo
712	608 format(1x,'k = ',i3,' t = ',f5.1,' w = ',f6.2,' frz/mlt = ',f5.1, &
713	' snomlt = ',f5.1)
714	!
715	!----------------------------------------------------------------
716	! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
717	!
718	do i = its, ite
719	if(fall(i,1).gt.0.) then
720	rainncv(i) = fall(i,1)delz(i,1)/denrdtcld*1000.
721	rain(i) = fall(i,1)delz(i,1)/denrdtcld*1000. &
722	+ rain(i)
723	endif
724	enddo
725	!
726	! if(lat.eq.latd) write(iun,601) latd,lond,loop,rain(lond)
727	601 format(1x,' ncloud3 lat lon loop : rain(mm) ',3i6,f20.2)
728	!
729	enddo ! big loops
730
731	END SUBROUTINE ncloud32D
732	! ...................................................................
733	real function rgmma(x)
734	!-------------------------------------------------------------------
735	IMPLICIT NONE
736	!-------------------------------------------------------------------
737	! rgmma function: use infinite product form
738	REAL :: euler
739	parameter (euler=0.577215664901532)
740	REAL :: x, y
741	INTEGER :: i
742
743	if(x.eq.1.)then
744	rgmma=0.
745	else
746	rgmma=xexp(eulerx)
747	do i=1,10000
748	y=float(i)
749	rgmma=rgmma(1.000+x/y)exp(-x/y)
750	enddo
751	rgmma=1./rgmma
752	endif
753	END FUNCTION rgmma
754	!
755	!--------------------------------------------------------------------------
756	real function fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
757	!--------------------------------------------------------------------------
758	IMPLICIT NONE
759	!--------------------------------------------------------------------------
760	real t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
761	xai,xbi,ttp,tr
762	INTEGER ice
763	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
764	ttp=t0c+0.01
765	dldt=cvap-cliq
766	xa=-dldt/rv
767	xb=xa+hvap/(rv*ttp)
768	dldti=cvap-cice
769	xai=-dldti/rv
770	xbi=xai+hsub/(rv*ttp)
771	tr=ttp/t
772	if(t.lt.ttp.and.ice.eq.1) then
773	fpvs=psat(trxai)exp(xbi*(1.-tr))
774	else
775	fpvs=psat(trxa)exp(xb*(1.-tr))
776	endif
777	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
778	END FUNCTION fpvs
779
780	!-------------------------------------------------------------------
781	SUBROUTINE ncloud3init(den0,denr,dens,cl,cpv,allowed_to_read)
782	!-------------------------------------------------------------------
783	IMPLICIT NONE
784	!-------------------------------------------------------------------
785	!.... constants which may not be tunable
786
787	REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
788	LOGICAL,INTENT(IN) :: allowed_to_read
789	REAL :: pi
790
791	pi = 4.*atan(1.)
792	xlv1 = cl-cpv
793
794	qc0 = 4./3.pidenrr03xncr/den0 ! 0.419e-3 -- .61e-3
795	qck1 = .1049.8peaut/(xncrdenr)*(1./3.)/xmyu ! 7.03
796	bvtr1 = 1.+bvtr
797	bvtr2 = 2.5+.5*bvtr
798	bvtr3 = 3.+bvtr
799	bvtr4 = 4.+bvtr
800	g1pbr = rgmma(bvtr1)
801	g3pbr = rgmma(bvtr3)
802	g4pbr = rgmma(bvtr4) ! 17.837825
803	g5pbro2 = rgmma(bvtr2) ! 1.8273
804	pvtr = avtr*g4pbr/6.
805	eacrr = 1.0
806	pacrr = pin0ravtrg3pbr.25*eacrr
807	precr1 = 2.pin0r*.78
808	precr2 = 2.pin0r.31avtr*.5g5pbro2
809	xm0 = (di0/dicon)**2
810	xmmax = (dimax/dicon)**2
811	!
812	bvts1 = 1.+bvts
813	bvts2 = 2.5+.5*bvts
814	bvts3 = 3.+bvts
815	bvts4 = 4.+bvts
816	g1pbs = rgmma(bvts1) !.8875
817	g3pbs = rgmma(bvts3)
818	g4pbs = rgmma(bvts4) ! 12.0786
819	g5pbso2 = rgmma(bvts2)
820	pvts = avts*g4pbs/6.
821	pacrs = pin0savtsg3pbs.25
822	precs1 = 4.n0s.65
823	precs2 = 4.n0s.44avts.5g5pbso2
824	pidn0r = pidenrn0r
825	pidn0s = pidensn0s
826	!
827	END SUBROUTINE ncloud3init
828
829	END MODULE module_mp_ncloud3
830

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

Download in other formats:

Original Format