Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_ysu.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
File size: 42.6 KB

Line
1	!wrf:model_layer:physics
2	!
3	!
4	!
5	module module_bl_ysu
6	contains
7	!
8	!-------------------------------------------------------------------
9	!
10	subroutine ysu(u3d,v3d,th3d,t3d,qv3d,qc3d,qi3d,p3d,pi3d, &
11	rublten,rvblten,rthblten, &
12	rqvblten,rqcblten,rqiblten, &
13	cp,g,rovcp,rd,rovg, &
14	dz8w,z,xlv,rv,psfc, &
15	znt,ust,zol,hol,hpbl,psim,psih, &
16	xland,hfx,qfx,tsk,gz1oz0,wspd,br, &
17	dt,dtmin,kpbl2d, &
18	svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt,&
19	exch_h, &
20	ids,ide, jds,jde, kds,kde, &
21	ims,ime, jms,jme, kms,kme, &
22	its,ite, jts,jte, kts,kte, &
23	!optional
24	regime )
25	!-------------------------------------------------------------------
26	implicit none
27	!-------------------------------------------------------------------
28	!-- u3d 3d u-velocity interpolated to theta points (m/s)
29	!-- v3d 3d v-velocity interpolated to theta points (m/s)
30	!-- th3d 3d potential temperature (k)
31	!-- t3d temperature (k)
32	!-- qv3d 3d water vapor mixing ratio (kg/kg)
33	!-- qc3d 3d cloud mixing ratio (kg/kg)
34	!-- qi3d 3d ice mixing ratio (kg/kg)
35	! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)
36	!-- p3d 3d pressure (pa)
37	!-- pi3d 3d exner function (dimensionless)
38	!-- rr3d 3d dry air density (kg/m^3)
39	!-- rublten u tendency due to
40	! pbl parameterization (m/s/s)
41	!-- rvblten v tendency due to
42	! pbl parameterization (m/s/s)
43	!-- rthblten theta tendency due to
44	! pbl parameterization (K/s)
45	!-- rqvblten qv tendency due to
46	! pbl parameterization (kg/kg/s)
47	!-- rqcblten qc tendency due to
48	! pbl parameterization (kg/kg/s)
49	!-- rqiblten qi tendency due to
50	! pbl parameterization (kg/kg/s)
51	!-- cp heat capacity at constant pressure for dry air (j/kg/k)
52	!-- g acceleration due to gravity (m/s^2)
53	!-- rovcp r/cp
54	!-- rd gas constant for dry air (j/kg/k)
55	!-- rovg r/g
56	!-- dz8w dz between full levels (m)
57	!-- z height above sea level (m)
58	!-- xlv latent heat of vaporization (j/kg)
59	!-- rv gas constant for water vapor (j/kg/k)
60	!-- psfc pressure at the surface (pa)
61	!-- znt roughness length (m)
62	!-- ust u* in similarity theory (m/s)
63	!-- zol z/l height over monin-obukhov length
64	!-- hol pbl height over monin-obukhov length
65	!-- hpbl pbl height (m)
66	!-- regime flag indicating pbl regime (stable, unstable, etc.)
67	!-- psim similarity stability function for momentum
68	!-- psih similarity stability function for heat
69	!-- xland land mask (1 for land, 2 for water)
70	!-- hfx upward heat flux at the surface (w/m^2)
71	!-- qfx upward moisture flux at the surface (kg/m^2/s)
72	!-- tsk surface temperature (k)
73	!-- gz1oz0 log(z/z0) where z0 is roughness length
74	!-- wspd wind speed at lowest model level (m/s)
75	!-- br bulk richardson number in surface layer
76	!-- dt time step (s)
77	!-- dtmin time step (minute)
78	!-- rvovrd r_v divided by r_d (dimensionless)
79	!-- svp1 constant for saturation vapor pressure (kpa)
80	!-- svp2 constant for saturation vapor pressure (dimensionless)
81	!-- svp3 constant for saturation vapor pressure (k)
82	!-- svpt0 constant for saturation vapor pressure (k)
83	!-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
84	!-- ep2 constant for specific humidity calculation
85	!-- karman von karman constant
86	!-- eomeg angular velocity of earths rotation (rad/s)
87	!-- stbolt stefan-boltzmann constant (w/m^2/k^4)
88	!-- ids start index for i in domain
89	!-- ide end index for i in domain
90	!-- jds start index for j in domain
91	!-- jde end index for j in domain
92	!-- kds start index for k in domain
93	!-- kde end index for k in domain
94	!-- ims start index for i in memory
95	!-- ime end index for i in memory
96	!-- jms start index for j in memory
97	!-- jme end index for j in memory
98	!-- kms start index for k in memory
99	!-- kme end index for k in memory
100	!-- its start index for i in tile
101	!-- ite end index for i in tile
102	!-- jts start index for j in tile
103	!-- jte end index for j in tile
104	!-- kts start index for k in tile
105	!-- kte end index for k in tile
106	!-------------------------------------------------------------------
107	!
108	integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
109	ims,ime, jms,jme, kms,kme, &
110	its,ite, jts,jte, kts,kte
111	!
112	real, intent(in ) :: dt,dtmin,cp,g,rovcp, &
113	rovg,rd,xlv,rv
114	!
115	real, intent(in ) :: svp1,svp2,svp3,svpt0
116	real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt
117	!
118	real, dimension( ims:ime, kms:kme, jms:jme ) , &
119	intent(in ) :: qv3d, &
120	qc3d, &
121	qi3d, &
122	p3d, &
123	pi3d, &
124	th3d, &
125	t3d, &
126	dz8w, &
127	z
128	!
129	real, dimension( ims:ime, kms:kme, jms:jme ) , &
130	intent(inout) :: rublten, &
131	rvblten, &
132	rthblten, &
133	rqvblten, &
134	rqcblten
135	real, dimension( ims:ime, kms:kme, jms:jme ) , &
136	intent(inout) :: exch_h
137	!
138	real, dimension( ims:ime, jms:jme ) , &
139	intent(in ) :: xland, &
140	hfx, &
141	qfx, &
142	psim, &
143	psih, &
144	gz1oz0, &
145	br, &
146	psfc, &
147	tsk
148	!
149	real, dimension( ims:ime, jms:jme ) , &
150	intent(inout) :: hol, &
151	ust, &
152	hpbl, &
153	znt, &
154	wspd, &
155	zol
156	!
157	real, dimension( ims:ime, kms:kme, jms:jme ) , &
158	intent(in ) :: u3d, &
159	v3d
160	!
161	integer, dimension( ims:ime, jms:jme ) , &
162	intent(out ) :: kpbl2d
163	!optional
164	real, dimension( ims:ime, jms:jme ) , &
165	optional , &
166	intent(inout) :: regime
167	!
168	real, dimension( ims:ime, kms:kme, jms:jme ) , &
169	optional , &
170	intent(inout) :: rqiblten
171	!
172	!local
173	integer :: i,j,k
174	!
175	do j = jts,jte
176	call ysu2d(J=j,ux=u3d(ims,kms,j),vx=v3d(ims,kms,j) &
177	,tx=t3d(ims,kms,j) &
178	,qx=qv3d(ims,kms,j),qcx=qc3d(ims,kms,j) &
179	,qix=qi3d(ims,kms,j) &
180	,p2d=p3d(ims,kms,j),utnp=rublten(ims,kms,j) &
181	,vtnp=rvblten(ims,kms,j) &
182	,ttnp=rthblten(ims,kms,j),qtnp=rqvblten(ims,kms,j) &
183	,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
184	,cp=cp,g=g,rovcp=rovcp,rd=rd,rovg=rovg &
185	,dz8w2d=dz8w(ims,kms,j),z2d=z(ims,kms,j) &
186	,xlv=xlv,rv=rv &
187	,psfcpa=psfc(ims,j),znt=znt(ims,j),ust=ust(ims,j) &
188	,zol=zol(ims,j),hol=hol(ims,j),hpbl=hpbl(ims,j) &
189	,regime=regime(ims,j),psim=psim(ims,j) &
190	,psih=psih(ims,j),xland=xland(ims,j) &
191	,hfx=hfx(ims,j),qfx=qfx(ims,j) &
192	,tsk=tsk(ims,j),gz1oz0=gz1oz0(ims,j) &
193	,wspd=wspd(ims,j),br=br(ims,j) &
194	,dt=dt,dtmin=dtmin,kpbl1d=kpbl2d(ims,j) &
195	,svp1=svp1,svp2=svp2,svp3=svp3,svpt0=svpt0 &
196	,ep1=ep1,ep2=ep2,karman=karman,eomeg=eomeg &
197	,stbolt=stbolt &
198	,exch_hx=exch_h(ims,kms,j) &
199	,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
200	,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
201	,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
202	do k = kts,kte
203	do i = its,ite
204	rthblten(i,k,j) = rthblten(i,k,j)/pi3d(i,k,j)
205	enddo
206	enddo
207	enddo
208	!
209	end subroutine ysu
210	!
211	!-------------------------------------------------------------------
212	!
213	subroutine ysu2d(j,ux,vx,tx,qx,qcx,qix,p2d, &
214	utnp,vtnp,ttnp, &
215	qtnp,qctnp,qitnp, &
216	cp,g,rovcp,rd,rovg, &
217	dz8w2d,z2d,xlv,rv,psfcpa, &
218	znt,ust,zol,hol,hpbl,psim,psih, &
219	xland,hfx,qfx,tsk,gz1oz0,wspd,br, &
220	dt,dtmin,kpbl1d, &
221	svp1,svp2,svp3,svpt0,ep1,ep2,karman,eomeg,stbolt,&
222	exch_hx, &
223	ids,ide, jds,jde, kds,kde, &
224	ims,ime, jms,jme, kms,kme, &
225	its,ite, jts,jte, kts,kte, &
226	!optional
227	regime )
228	!-------------------------------------------------------------------
229	implicit none
230	!-------------------------------------------------------------------
231	!
232	! this code is a revised vertical diffusion package ("ysupbl")
233	! with a nonlocal turbulent mixing in the pbl after "mrfpbl".
234	! the ysupbl is based on the study of noh et al(2003) and
235	! accumulated realism of the behavior of the troen and mahrt
236	! (1986) concept implemented by hong and pan(1996). the major
237	! ingredient of the ysupbl is the inclusion of an explicit
238	! treatment of the entrainment processes at the entrainment layer.
239	! this routine uses an implicit approach for vertical flux
240	! divergence and does not require "miter" timesteps.
241	! it includes vertical diffusion in the stable atmosphere
242	! and moist vertical diffusion in clouds.
243	! surface fluxes calculated as in hirpbl.
244	! 5-layer soil model option required in slab due to long timestep
245	!
246	! mrfpbl:
247	! coded by song-you hong (ncep), implemented by jimy dudhia (ncar)
248	! fall 1996
249	!
250	! ysupbl:
251	! coded by song-you hong (yonsei university) and implemented by
252	! song-you hong (yonsei university) and jimy dudhia (ncar)
253	! summer 2002
254	!
255	! references:
256	!
257	! hong and pan (1996), mon. wea. rev.
258	! hong, noh, and dudhia (2004), mon. wea. rev, to be submitted
259	! dudhia and hong (2004), mon. wea. rev, to be submitted
260	! noh, chun, hong, and raasch (2003), boundary layer met.
261	! troen and mahrt (1986), boundary layer met.
262	!
263	!-------------------------------------------------------------------
264	!
265	real rlam,prmin,prmax,xkzmin,xkzmax,rimin,brcr, &
266	bfac,pfac,sfcfrac,ckz,zfmin,aphi5,aphi16,gamcrt, &
267	gamcrq,xka
268	!
269	parameter (xkzmin = 0.01,xkzmax = 1000.,rimin = -100.)
270	integer ncloud
271	parameter (ncloud = 3)
272	real afac, phifac, qmin
273	real d1,d2,d3
274	real h1,h2
275	parameter (rlam = 30.,prmin = 0.25,prmax = 4.)
276	parameter (brcr = 0.0,bfac = 6.8,pfac = 2.0,sfcfrac = 0.1)
277	parameter (afac = 6.8,phifac = 8.,qmin=1.e-2)
278	parameter (d1 = 0.02, d2 = 0.05, d3 = 0.001)
279	parameter (h1 = 0.33333335, h2 = 0.6666667)
280	parameter (ckz = 0.001,zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.)
281	parameter (gamcrt = 3.,gamcrq = 2.e-3)
282	parameter (xka = 2.4e-5)
283	!wig 16-sep-2005: Parameter added to set a minimum pbl height. To disable
284	! this option, set hpblmin=0.
285	real, parameter :: hpblmin=0.
286	!
287	integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
288	ims,ime, jms,jme, kms,kme, &
289	its,ite, jts,jte, kts,kte, j
290	!
291	real, intent(in ) :: dt,dtmin,cp,g,rovcp, &
292	rovg,rd,xlv,rv
293	!
294	real, intent(in ) :: svp1,svp2,svp3,svpt0
295	real, intent(in ) :: ep1,ep2,karman,eomeg,stbolt
296	!
297	real, dimension( ims:ime, kms:kme ), &
298	intent(in) :: dz8w2d, &
299	z2d
300	!
301	real, dimension( ims:ime, kms:kme ) , &
302	intent(in ) :: tx, &
303	qx, &
304	qcx, &
305	qix, &
306	p2d
307	!
308	real, dimension( ims:ime, kms:kme ) , &
309	intent(inout) :: utnp, &
310	vtnp, &
311	ttnp, &
312	qtnp, &
313	qctnp, &
314	qitnp
315	!
316	real, dimension( ims:ime ) , &
317	intent(inout) :: hol, &
318	ust, &
319	hpbl, &
320	znt
321	real, dimension( ims:ime ) , &
322	intent(in ) :: xland, &
323	hfx, &
324	qfx
325	!
326	real, dimension( ims:ime ), intent(inout) :: wspd
327	real, dimension( ims:ime ), intent(in ) :: br
328	!
329	real, dimension( ims:ime ), intent(in ) :: psim, &
330	psih
331	real, dimension( ims:ime ), intent(in ) :: gz1oz0
332	!
333	real, dimension( ims:ime ), intent(in ) :: psfcpa
334	real, dimension( ims:ime ), intent(in ) :: tsk
335	real, dimension( ims:ime ), intent(inout) :: zol
336	integer, dimension( ims:ime ), intent(out ) :: kpbl1d
337	!
338	!
339	real, dimension( ims:ime, kms:kme ) , &
340	intent(in ) :: ux, &
341	vx
342	!optional
343	real, dimension( ims:ime ) , &
344	optional , &
345	intent(inout) :: regime
346	!
347	! local vars
348	!
349	real, dimension( its:ite, kts:kte+1 ) :: zq
350	!
351	real, dimension( its:ite, kts:kte ) :: &
352	thx,thvx, &
353	dzq,dza, &
354	za, &
355	uxs,vxs, &
356	thxs,qxs, &
357	qcxs,qixs
358	!
359	real, dimension( its:ite ) :: qixsv,rhox, &
360	govrth, &
361	thxsv, &
362	uxsv,vxsv, &
363	qxsv,qcxsv, &
364	qgh,tgdsa,ps
365	!
366	real, dimension( its:ite ) :: &
367	zl1,thermal, &
368	wscale,hgamt, &
369	hgamq,brdn, &
370	brup,phim, &
371	phih,cpm, &
372	dusfc,dvsfc, &
373	dtsfc,dqsfc, &
374	thgb,prpbl, &
375	wspd1
376	!
377	real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, &
378	f1,f2, &
379	r1,r2, &
380	ad,au, &
381	cu, &
382	al, &
383	zfac, &
384	scr4
385	!jdf added exch_hx
386	real, dimension( ims:ime, kms:kme ) , &
387	intent(inout) :: exch_hx
388	real, dimension( its:ite, kts:kte, ncloud) :: r3,f3
389	!
390	logical, dimension( its:ite ) :: pblflg, &
391	sfcflg, &
392	stable
393	!
394	integer :: n,i,k,l,nzol,imvdif,ic
395	integer :: klpbl
396	!
397	!wig 16-sep-2005: Added kmin for setting minimum pbl height
398	integer, dimension( its:ite ) :: kpbl, &
399	kmin
400	!
401	real :: zoln,x,y,thcon,tvcon,e1,dtstep
402	real :: zl,tskv,dthvdz,dthvm,vconv,rzol
403	real :: dtthx,psix,dtg,psiq,ustm
404	real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum
405	real :: xkzo,ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
406	real :: brint,dtodsd,rdz,dsdzt,dsdzq,dsdz2,ttend,qtend
407	real :: utend,vtend,qctend,qitend,tgc,dtodsu,govrthv
408	!
409	real, dimension( its:ite, kts:kte ) :: wscalek, &
410	xkzml,xkzhl, &
411	zfacent,entfac
412	real, dimension( its:ite ) :: ust3, &
413	wstar3,wstar, &
414	hgamu,hgamv, &
415	wm2, we, &
416	bfxpbl, &
417	hfxpbl,qfxpbl, &
418	ufxpbl,vfxpbl, &
419	delta,dthvx
420	real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
421	dsdzu,dsdzv,wm3,dthx,dqx
422	!
423	!----------------------------------------------------------------------
424	!
425	klpbl = kte
426	!
427	!-- imvdif imvdif = 1 for moist adiabat vertical diffusion
428	imvdif = 1
429	!
430	!----convert ground temperature to potential temperature:
431	!
432	do i = its,ite
433	tgdsa(i) = tsk(i)
434	ps(i) = psfcpa(i)/1000. ! ps psfc cmb
435	thgb(i) = tsk(i)(100./ps(i))*rovcp
436	enddo
437	!
438	! scr4(i,k) store virtual temperature.
439	!
440	do k = kts,kte
441	do i = its,ite
442	thcon = (100000./p2d(i,k))**rovcp
443	thx(i,k) = tx(i,k)*thcon
444	scr4(i,k) = tx(i,k)
445	thvx(i,k) = thx(i,k)
446	enddo
447	enddo
448	!
449	do i = its,ite
450	qgh(i) = 0.
451	cpm(i) = cp
452	enddo
453	!
454	do k = kts,kte
455	do i = its,ite
456	tvcon = (1.+ep1*qx(i,k))
457	thvx(i,k) = thx(i,k)*tvcon
458	scr4(i,k) = tx(i,k)*tvcon
459	enddo
460	enddo
461	!
462	do i = its,ite
463	e1 = svp1exp(svp2(tgdsa(i)-svpt0)/(tgdsa(i)-svp3))
464	qgh(i) = ep2*e1/(ps(i)-e1)
465	cpm(i) = cp(1.+0.8qx(i,1))
466	enddo
467	!
468	!-----compute the height of full- and half-sigma levels above ground
469	! level, and the layer thicknesses.
470	!
471	do i = its,ite
472	zq(i,1) = 0.
473	rhox(i) = ps(i)1000./(rdscr4(i,1))
474	enddo
475	!
476	do k = kts,kte
477	do i = its,ite
478	zq(i,k+1) = dz8w2d(i,k)+zq(i,k)
479	enddo
480	enddo
481	!
482	do k = kts,kte
483	do i = its,ite
484	za(i,k) = 0.5*(zq(i,k)+zq(i,k+1))
485	dzq(i,k) = zq(i,k+1)-zq(i,k)
486	enddo
487	enddo
488	!
489	do i = its,ite
490	dza(i,1) = za(i,1)
491	enddo
492	!
493	do k = kts+1,kte
494	do i = its,ite
495	dza(i,k) = za(i,k)-za(i,k-1)
496	enddo
497	enddo
498	!
499	do i = its,ite
500	govrth(i) = g/thx(i,1)
501	enddo
502	!
503	!wig 16-sep-2005: Determine the k level associated with hpblmin.
504	if( hpblmin >= 1. ) then
505	do i = its,ite
506	kmin(i) = 1
507	do k = kte-1,kts,-1
508	if( zq(i,k) <= hpblmin ) then
509	kmin(i) = k
510	exit
511	end if
512	end do
513	end do
514	end if
515	!
516	!-----initialize vertical tendencies and
517	!
518	do i = its,ite
519	do k = kts,kte
520	utnp(i,k) = 0.
521	vtnp(i,k) = 0.
522	ttnp(i,k) = 0.
523	enddo
524	enddo
525	!
526	do k = kts,kte
527	do i = its,ite
528	qtnp(i,k) = 0.
529	enddo
530	enddo
531	!
532	do k = kts,kte
533	do i = its,ite
534	qctnp(i,k) = 0.
535	qitnp(i,k) = 0.
536	enddo
537	enddo
538	!
539	do i = its,ite
540	wspd1(i) = sqrt(ux(i,1)ux(i,1)+vx(i,1)vx(i,1))+1.e-9
541	enddo
542	!
543	!---- compute vertical diffusion
544	!
545	! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
546	! compute preliminary variables
547	!
548	dtstep = dt
549	dt2 = 2.*dtstep
550	rdt = 1./dt2
551	!
552	do i = its,ite
553	bfxpbl(i) = 0.0
554	hfxpbl(i) = 0.0
555	qfxpbl(i) = 0.0
556	ufxpbl(i) = 0.0
557	vfxpbl(i) = 0.0
558	hgamu(i) = 0.0
559	hgamv(i) = 0.0
560	delta(i) = 0.0
561	enddo
562	!
563	do k = kts,klpbl
564	do i = its,ite
565	wscalek(i,k) = 0.0
566	enddo
567	enddo
568	!
569	do k = kts,klpbl
570	do i = its,ite
571	zfac(i,k) = 0.0
572	enddo
573	enddo
574	!
575	do i = its,ite
576	hgamt(i) = 0.
577	hgamq(i) = 0.
578	wscale(i) = 0.
579	kpbl(i) = 1
580	hpbl(i) = zq(i,1)
581	zl1(i) = za(i,1)
582	thermal(i)= thvx(i,1)
583	pblflg(i) = .true.
584	sfcflg(i) = .true.
585	!wig 12-aug-2004: Turn off check for sfc stability if using a minimum
586	! pbl height > 0.
587	if(br(i).gt.0.0 .and. hpblmin<=1.) sfcflg(i) = .false.
588	enddo
589	!
590	! compute the first guess of pbl height
591	!
592	do i = its,ite
593	stable(i) = .false.
594	brup(i) = br(i)
595	enddo
596	!
597	do k = 2,klpbl
598	do i = its,ite
599	if(.not.stable(i))then
600	brdn(i) = brup(i)
601	spdk2 = max(ux(i,k)2+vx(i,k)2,1.)
602	brup(i) = (thvx(i,k)-thermal(i))(gza(i,k)/thvx(i,1))/spdk2
603	kpbl(i) = k
604	stable(i) = brup(i).gt.brcr
605	endif
606	enddo
607	enddo
608	!
609	do i = its,ite
610	k = kpbl(i)
611	if(brdn(i).ge.brcr)then
612	brint = 0.
613	elseif(brup(i).le.brcr)then
614	brint = 1.
615	else
616	brint = (brcr-brdn(i))/(brup(i)-brdn(i))
617	endif
618	hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
619
620	!wig 16-sep-2005: rig a minimum PBL heigt
621	if( hpblmin >= 1. .and. hpbl(i).lt.hpblmin) then
622	kpbl(i) = kmin(i)
623	hpbl(i) = hpblmin
624	else
625	if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
626	end if
627
628	if(kpbl(i).le.1) pblflg(i) = .false.
629	enddo
630	!
631	do i = its,ite
632	fm = gz1oz0(i)-psim(i)
633	fh = gz1oz0(i)-psih(i)
634	hol(i) = max(br(i)fmfm/fh,rimin)
635	if(sfcflg(i))then
636	hol(i) = min(hol(i),-zfmin)
637	else
638	hol(i) = max(hol(i),zfmin)
639	endif
640	!
641	hol1 = hol(i)hpbl(i)/zl1(i)sfcfrac
642	hol(i) = -hol(i)*hpbl(i)/zl1(i)
643	if(sfcflg(i))then
644	phim(i) = (1.-aphi16hol1)*(-1./4.)
645	phih(i) = (1.-aphi16hol1)*(-1./2.)
646	bfx0 = max(hfx(i)/rhox(i)/cpm(i) &
647	+ep1thx(i,1)qfx(i)/rhox(i),0.)
648	hfx0 = max(hfx(i)/rhox(i)/cpm(i),0.)
649	qfx0 = max(ep1thx(i,1)qfx(i)/rhox(i),0.)
650	wstar3(i) = (govrth(i)bfx0hpbl(i))
651	wstar(i) = (wstar3(i))**h1
652	else
653	phim(i) = (1.+aphi5*hol1)
654	phih(i) = phim(i)
655	wstar(i) = 0.
656	wstar3(i) = 0.
657	endif
658	ust3(i) = ust(i)**3.
659	wscale(i) = (ust3(i)+phifackarmanwstar3(i)0.5)*h1
660	wscale(i) = min(wscale(i),ust(i)*aphi16)
661	wscale(i) = max(wscale(i),ust(i)/aphi5)
662	enddo
663	!
664	! compute the surface variables for pbl height estimation
665	! under unstable conditions
666	!
667	do i = its,ite
668	if(sfcflg(i))then
669	gamfac = bfac/rhox(i)/wscale(i)
670	hgamt(i) = min(gamfac*hfx(i)/cpm(i),gamcrt)
671	hgamq(i) = min(gamfac*qfx(i),gamcrq)
672	vpert = (hgamt(i)+ep1thx(i,1)hgamq(i))/bfac*afac
673	thermal(i) = thermal(i)+max(vpert,0.)
674	hgamt(i) = max(hgamt(i),0.0)
675	hgamq(i) = 0.0
676	brint = -15.9ust(i)ust(i)/wspd(i)*wstar3(i) &
677	/(wscale(i)**4.)
678	hgamu(i) = brint*ux(i,1)
679	hgamv(i) = brint*vx(i,1)
680	else
681	pblflg(i) = .false.
682	endif
683	enddo
684	!
685	! enhance the pbl height by considering the thermal
686	!
687	do i = its,ite
688	if(pblflg(i))then
689	kpbl(i) = 1
690	hpbl(i) = zq(i,1)
691	endif
692	enddo
693	!
694	do i = its,ite
695	if(pblflg(i))then
696	stable(i) = .false.
697	brup(i) = br(i)
698	endif
699	enddo
700	!
701	do k = 2,klpbl
702	do i = its,ite
703	if(.not.stable(i).and.pblflg(i))then
704	brdn(i) = brup(i)
705	spdk2 = max((ux(i,k)2+vx(i,k)2),1.)
706	brup(i) = (thvx(i,k)-thermal(i))(gza(i,k)/thvx(i,1))/spdk2
707	kpbl(i) = k
708	stable(i) = brup(i).gt.brcr
709	endif
710	enddo
711	enddo
712	!
713	do i = its,ite
714	if(pblflg(i))then
715	k = kpbl(i)
716	if(brdn(i).ge.brcr)then
717	brint = 0.
718	elseif(brup(i).le.brcr)then
719	brint = 1.
720	else
721	brint = (brcr-brdn(i))/(brup(i)-brdn(i))
722	endif
723	hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
724
725	!wig 16-sep-2005: rig a minimum PBL heigt
726	if(hpblmin >= 1. .and. hpbl(i).lt.hpblmin) then
727	kpbl(i) = kmin(i)
728	hpbl(i) = hpblmin
729	else
730	if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
731	end if
732
733	if(kpbl(i).le.1) pblflg(i) = .false.
734	endif
735	enddo
736	!
737	! estimate the entrainment parameters
738	!
739	do i = its,ite
740	if(pblflg(i)) then
741	k = kpbl(i) - 1
742	ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k))+(vx(i,k+1)- &
743	vx(i,k))(vx(i,k+1)-vx(i,k)))/(dza(i,k+1)dza(i,k+1))+ &
744	1.e-9
745	govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
746	ri = govrthv(thvx(i,k+1)-thvx(i,k))/(ssdza(i,k+1))
747	if(imvdif.eq.1)then
748	if((qcx(i,k)+qix(i,k)).gt.0.01e-3.and.(qcx(i,k+1)+ &
749	qix(i,k+1)).gt.0.01e-3)then
750	! in cloud
751	qmean = 0.5*(qx(i,k)+qx(i,k+1))
752	tmean = 0.5*(tx(i,k)+tx(i,k+1))
753	alph = xlv*qmean/rd/tmean
754	chi = xlvxlvqmean/cp/rv/tmean/tmean
755	ri = (1.+alph)(ri-gg/ss/tmean/cp*((chi-alph)/(1.+chi)))
756	endif
757	endif
758	prpbl(i) = 1.0
759	wm3 = wstar3(i) + 5. * ust3(i)
760	wm2(i) = wm3**h2
761	bfxpbl(i) = -0.15thvx(i,1)/gwm3/hpbl(i)
762	dthvx(i) = max(thvx(i,k+1)-thvx(i,k),qmin)
763	dthx = max(thx(i,k+1)-thx(i,k),qmin)
764	dqx = min(qx(i,k+1)-qx(i,k),0.0)
765	we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
766	hfxpbl(i) = we(i)*dthx
767	qfxpbl(i) = we(i)*dqx
768	!
769	dux = ux(i,k+1)-ux(i,k)
770	dvx = vx(i,k+1)-vx(i,k)
771	if(dux.gt.qmin) then
772	ufxpbl(i) = max(prpbl(i)we(i)dux,-ust(i)*ust(i))
773	elseif(dux.lt.-qmin) then
774	ufxpbl(i) = min(prpbl(i)we(i)dux,ust(i)*ust(i))
775	else
776	ufxpbl(i) = 0.0
777	endif
778	if(dvx.gt.qmin) then
779	vfxpbl(i) = max(prpbl(i)we(i)dvx,-ust(i)*ust(i))
780	elseif(dvx.lt.-qmin) then
781	vfxpbl(i) = min(prpbl(i)we(i)dvx,ust(i)*ust(i))
782	else
783	vfxpbl(i) = 0.0
784	endif
785	delb = govrth(i)d3hpbl(i)
786	delta(i) = min(d1hpbl(i) + d2wm2(i)/delb,100.)
787	endif
788	enddo
789	!
790	do k = kts,klpbl
791	do i = its,ite
792	if(pblflg(i).and.k.ge.kpbl(i))then
793	entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
794	else
795	entfac(i,k) = 1.e30
796	endif
797	enddo
798	enddo
799	!
800	! compute diffusion coefficients below pbl
801	!
802	do k = kts,klpbl
803	do i = its,ite
804	if(k.lt.kpbl(i)) then
805	zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i)) &
806	/(hpbl(i)-zl1(i))),zfmin),1.)
807	xkzo = ckz*dza(i,k+1)
808	zfacent(i,k) = (1.-zfac(i,k))**3.
809	prnumfac = -3.(max(zq(i,k+1)-sfcfrachpbl(i),0.))**2. &
810	/hpbl(i)**2.
811	prnum = (phih(i)/phim(i)+bfackarmansfcfrac)
812	prnum = 1. + (prnum-1.)*exp(prnumfac)
813	prnum = min(prnum,prmax)
814	prnum = max(prnum,prmin)
815	wscalek(i,k) = (ust3(i)+phifackarmanwstar3(i) &
816	(1.-zfac(i,k)))*h1
817	xkzm(i,k) = xkzo+wscalek(i,k)karmanzq(i,k+1) &
818	zfac(i,k)*pfac
819	xkzh(i,k) = xkzm(i,k)/prnum
820	xkzm(i,k) = min(xkzm(i,k),xkzmax)
821	xkzm(i,k) = max(xkzm(i,k),xkzmin)
822	xkzh(i,k) = min(xkzh(i,k),xkzmax)
823	xkzh(i,k) = max(xkzh(i,k),xkzmin)
824	!wig 16-sep-2005: This code is unnecessary. exch_hx will get overwritten
825	! with the same values further down in this routine
826	! exch_hx(i,k) = xkzh(i,k)
827	endif
828	enddo
829	enddo
830	!
831	! compute diffusion coefficients over pbl (free atmosphere)
832	!
833	do k = kts,kte-1
834	do i = its,ite
835	xkzo = ckz*dza(i,k+1)
836	if(k.ge.kpbl(i)) then
837	ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k))+(vx(i,k+1)- &
838	vx(i,k))(vx(i,k+1)-vx(i,k)))/(dza(i,k+1)dza(i,k+1))+ &
839	1.e-9
840	govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
841	ri = govrthv(thvx(i,k+1)-thvx(i,k))/(ssdza(i,k+1))
842	if(imvdif.eq.1)then
843	if((qcx(i,k)+qix(i,k)).gt.0.01e-3.and.(qcx(i,k+1)+ &
844	qix(i,k+1)).gt.0.01e-3)then
845	! in cloud
846	qmean = 0.5*(qx(i,k)+qx(i,k+1))
847	tmean = 0.5*(tx(i,k)+tx(i,k+1))
848	alph = xlv*qmean/rd/tmean
849	chi = xlvxlvqmean/cp/rv/tmean/tmean
850	ri = (1.+alph)(ri-gg/ss/tmean/cp*((chi-alph)/(1.+chi)))
851	endif
852	endif
853	zk = karman*zq(i,k+1)
854	rl2 = (zkrlam/(rlam+zk))*2
855	dk = rl2*sqrt(ss)
856	if(ri.lt.0.)then
857	! unstable regime
858	sri = sqrt(-ri)
859	xkzm(i,k) = xkzo+dk(1+8.(-ri)/(1+1.746*sri))
860	xkzh(i,k) = xkzo+dk(1+8.(-ri)/(1+1.286*sri))
861	else
862	! stable regime
863	xkzh(i,k) = xkzo+dk/(1+5.ri)*2
864	prnum = 1.0+2.1*ri
865	prnum = min(prnum,prmax)
866	xkzm(i,k) = (xkzh(i,k)-xkzo)*prnum+xkzo
867	endif
868	!
869	xkzm(i,k) = min(xkzm(i,k),xkzmax)
870	xkzm(i,k) = max(xkzm(i,k),xkzmin)
871	xkzh(i,k) = min(xkzh(i,k),xkzmax)
872	xkzh(i,k) = max(xkzh(i,k),xkzmin)
873	xkzml(i,k) = xkzm(i,k)
874	xkzhl(i,k) = xkzh(i,k)
875	!wig 16-sep-2005: This code is unnecessary. exch_hx will get overwritten
876	! further down in this routine after adjustments are made
877	! for entrainment.
878	! exch_hx(i,k) = xkzh(i,k)
879	endif
880	enddo
881	enddo
882	!
883	! compute tridiagonal matrix elements for heat and moisture, and clouds
884	!
885	do i = its,ite
886	do k = kts,kte
887	au(i,k) = 0.
888	al(i,k) = 0.
889	ad(i,k) = 0.
890	f1(i,k) = 0.
891	enddo
892	enddo
893	!
894	do ic = 1,ncloud
895	do i = its,ite
896	do k = kts,kte
897	f3(i,k,ic) = 0.
898	enddo
899	enddo
900	enddo
901	!
902	do i = its,ite
903	ad(i,1) = 1.
904	f1(i,1) = thx(i,1)+hfx(i)/(rhox(i)cpm(i))/zq(i,2)dt2
905	f3(i,1,1) = qx(i,1)+qfx(i)/(rhox(i))/zq(i,2)*dt2
906	enddo
907	!
908	if(ncloud.ge.2) then
909	do ic = 2,ncloud
910	do i = its,ite
911	if(ic.eq.2) then
912	f3(i,1,ic) = qcx(i,1)
913	elseif(ic.eq.3) then
914	f3(i,1,ic) = qix(i,1)
915	endif
916	enddo
917	enddo
918	endif
919	!
920	do k = kts,kte-1
921	do i = its,ite
922	dtodsd = dt2/dz8w2d(i,k)
923	dtodsu = dt2/dz8w2d(i,k+1)
924	rdz = 1./dza(i,k+1)
925	if(pblflg(i).and.k.lt.kpbl(i)) then
926	dsdzt = xkzh(i,k)*(-hgamt(i)/hpbl(i) &
927	-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
928	dsdzq = xkzh(i,k)*(-hgamq(i)/hpbl(i) &
929	-qfxpbl(i)*zfacent(i,k)/xkzh(i,k))
930	f1(i,k) = f1(i,k)+dtodsd*dsdzt
931	f1(i,k+1) = thx(i,k+1)-dtodsu*dsdzt
932	f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
933	f3(i,k+1,1) = qx(i,k+1)-dtodsu*dsdzq
934	elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
935	xkzh(i,k) = -we(i)dza(i,kpbl(i))exp(-entfac(i,k))
936	xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
937	xkzh(i,k) = min(xkzh(i,k),xkzmax)
938	xkzh(i,k) = max(xkzh(i,k),xkzmin)
939	f1(i,k+1) = thx(i,k+1)
940	f3(i,k+1,1) = qx(i,k+1)
941	else
942	f1(i,k+1) = thx(i,k+1)
943	f3(i,k+1,1) = qx(i,k+1)
944	endif
945	dsdz2 = xkzh(i,k)*rdz
946	au(i,k) = -dtodsd*dsdz2
947	al(i,k) = -dtodsu*dsdz2
948	ad(i,k) = ad(i,k)-au(i,k)
949	ad(i,k+1) = 1.-al(i,k)
950	exch_hx(i,k) = xkzh(i,k)
951	enddo
952	enddo
953	!
954	if(ncloud.ge.2) then
955	do ic = 2,ncloud
956	do k = kts,kte-1
957	do i = its,ite
958	if(ic.eq.2) then
959	f3(i,k+1,ic) = qcx(i,k+1)
960	elseif(ic.eq.3) then
961	f3(i,k+1,ic) = qix(i,k+1)
962	endif
963	enddo
964	enddo
965	enddo
966	endif
967
968	! copies here to avoid duplicate input args for tridin
969
970	do k = kts,kte
971	do i = its,ite
972	cu(i,k) = au(i,k)
973	r1(i,k) = f1(i,k)
974	enddo
975	enddo
976	!
977	do ic = 1,ncloud
978	do k = kts,kte
979	do i = its,ite
980	r3(i,k,ic) = f3(i,k,ic)
981	enddo
982	enddo
983	enddo
984	!
985	!
986	! solve tridiagonal problem for heat and moisture, and clouds
987	!
988	call tridin(al,ad,cu,r1,r3,au,f1,f3, &
989	its,ite,kts,kte,ncloud )
990	!
991	! recover tendencies of heat and moisture
992	!
993	do k = kte,kts,-1
994	do i = its,ite
995	ttend = (f1(i,k)-thx(i,k))rdt(tx(i,k)/thx(i,k))
996	qtend = (f3(i,k,1)-qx(i,k))*rdt
997	ttnp(i,k) = ttnp(i,k)+ttend
998	qtnp(i,k) = qtnp(i,k)+qtend
999	enddo
1000	enddo
1001	!
1002	if(ncloud.ge.2) then
1003	do ic = 2,ncloud
1004	do k = kte,kts,-1
1005	do i = its,ite
1006	if(ic.eq.2) then
1007	qctend = (f3(i,k,ic)-qcx(i,k))*rdt
1008	qctnp(i,k) = qctnp(i,k)+qctend
1009	elseif(ic.eq.3) then
1010	qitend = (f3(i,k,ic)-qix(i,k))*rdt
1011	qitnp(i,k) = qitnp(i,k)+qitend
1012	endif
1013	enddo
1014	enddo
1015	enddo
1016	endif
1017	!
1018	! compute tridiagonal matrix elements for momentum
1019	!
1020	do i = its,ite
1021	do k = kts,kte
1022	au(i,k) = 0.
1023	al(i,k) = 0.
1024	ad(i,k) = 0.
1025	f1(i,k) = 0.
1026	f2(i,k) = 0.
1027	enddo
1028	enddo
1029	!
1030	do i = its,ite
1031	ad(i,1) = 1.
1032	f1(i,1) = ux(i,1)-ux(i,1)/wspd1(i)ust(i)ust(i)/zq(i,2)*dt2 &
1033	(wspd1(i)/wspd(i))*2
1034	f2(i,1) = vx(i,1)-vx(i,1)/wspd1(i)ust(i)ust(i)/zq(i,2)*dt2 &
1035	(wspd1(i)/wspd(i))*2
1036	enddo
1037	!
1038	do k = kts,kte-1
1039	do i = its,ite
1040	dtodsd = dt2/dz8w2d(i,k)
1041	dtodsu = dt2/dz8w2d(i,k+1)
1042	rdz = 1./dza(i,k+1)
1043	if(pblflg(i).and.k.lt.kpbl(i))then
1044	dsdzu=xkzm(i,k)*(-hgamu(i)/hpbl(i) &
1045	-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1046	dsdzv=xkzm(i,k)*(-hgamv(i)/hpbl(i) &
1047	-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1048	f1(i,k) = f1(i,k)+dtodsd*dsdzu
1049	f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1050	f2(i,k) = f2(i,k)+dtodsd*dsdzv
1051	f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1052	elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1053	xkzm(i,k) = prpbl(i)*xkzh(i,k)
1054	xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1055	xkzm(i,k)=min(xkzm(i,k),xkzmax)
1056	xkzm(i,k)=max(xkzm(i,k),xkzmin)
1057	f1(i,k+1)=ux(i,k+1)
1058	f2(i,k+1)=vx(i,k+1)
1059	else
1060	f1(i,k+1) = ux(i,k+1)
1061	f2(i,k+1) = vx(i,k+1)
1062	endif
1063	dsdz2 = xkzm(i,k)*rdz
1064	au(i,k) = -dtodsd*dsdz2
1065	al(i,k) = -dtodsu*dsdz2
1066	ad(i,k) = ad(i,k)-au(i,k)
1067	ad(i,k+1) = 1.-al(i,k)
1068	enddo
1069	enddo
1070
1071	! copies here to avoid duplicate input args for tridin
1072
1073	do k = kts,kte
1074	do i = its,ite
1075	cu(i,k) = au(i,k)
1076	r1(i,k) = f1(i,k)
1077	r2(i,k) = f2(i,k)
1078	enddo
1079	enddo
1080	!
1081	!
1082	! solve tridiagonal problem for momentum
1083	!
1084	call tridin(al,ad,cu,r1,r2,au,f1,f2, &
1085	its,ite,kts,kte,1 )
1086	!
1087	! recover tendencies of momentum
1088	!
1089	do k = kte,kts,-1
1090	do i = its,ite
1091	utend = (f1(i,k)-ux(i,k))*rdt
1092	vtend = (f2(i,k)-vx(i,k))*rdt
1093	utnp(i,k) = utnp(i,k)+utend
1094	vtnp(i,k) = vtnp(i,k)+vtend
1095	enddo
1096	enddo
1097	!
1098	!---- end of vertical diffusion
1099	!
1100	do i = its,ite
1101	kpbl1d(i) = kpbl(i)
1102	enddo
1103	!
1104	end subroutine ysu2d
1105	!
1106	subroutine tridin(cl,cm,cu,r1,r2,au,f1,f2, &
1107	its,ite,kts,kte,nt )
1108	!----------------------------------------------------------------
1109	implicit none
1110	!----------------------------------------------------------------
1111	!
1112	integer, intent(in ) :: its,ite, kts,kte, nt
1113	!
1114	real, dimension( its:ite, kts+1:kte+1 ) , &
1115	intent(in ) :: cl
1116	!
1117	real, dimension( its:ite, kts:kte ) , &
1118	intent(in ) :: cm, &
1119	r1
1120	real, dimension( its:ite, kts:kte,nt ) , &
1121	intent(in ) :: r2
1122	!
1123	real, dimension( its:ite, kts:kte ) , &
1124	intent(inout) :: au, &
1125	cu, &
1126	f1
1127	real, dimension( its:ite, kts:kte,nt ) , &
1128	intent(inout) :: f2
1129	!
1130	real :: fk
1131	integer :: i,k,l,n,it
1132	!
1133	!----------------------------------------------------------------
1134	!
1135	l = ite
1136	n = kte
1137	!
1138	do i = its,l
1139	fk = 1./cm(i,1)
1140	au(i,1) = fk*cu(i,1)
1141	f1(i,1) = fk*r1(i,1)
1142	enddo
1143	do it = 1,nt
1144	do i = its,l
1145	fk = 1./cm(i,1)
1146	f2(i,1,it) = fk*r2(i,1,it)
1147	enddo
1148	enddo
1149	do k = kts+1,n-1
1150	do i = its,l
1151	fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1152	au(i,k) = fk*cu(i,k)
1153	f1(i,k) = fk(r1(i,k)-cl(i,k)f1(i,k-1))
1154	enddo
1155	enddo
1156	do it = 1,nt
1157	do k = kts+1,n-1
1158	do i = its,l
1159	fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1160	f2(i,k,it) = fk(r2(i,k,it)-cl(i,k)f2(i,k-1,it))
1161	enddo
1162	enddo
1163	enddo
1164	do i = its,l
1165	fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1166	f1(i,n) = fk(r1(i,n)-cl(i,n)f1(i,n-1))
1167	enddo
1168	do it = 1,nt
1169	do i = its,l
1170	fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1171	f2(i,n,it) = fk(r2(i,n,it)-cl(i,n)f2(i,n-1,it))
1172	enddo
1173	enddo
1174	do k = n-1,kts,-1
1175	do i = its,l
1176	f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
1177	enddo
1178	enddo
1179	do it = 1,nt
1180	do k = n-1,kts,-1
1181	do i = its,l
1182	f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1183	enddo
1184	enddo
1185	enddo
1186	!
1187	end subroutine tridin
1188	!
1189	subroutine ysuinit(rublten,rvblten,rthblten,rqvblten, &
1190	rqcblten,rqiblten,p_qi,p_first_scalar, &
1191	restart, allowed_to_read, &
1192	ids, ide, jds, jde, kds, kde, &
1193	ims, ime, jms, jme, kms, kme, &
1194	its, ite, jts, jte, kts, kte )
1195	!-------------------------------------------------------------------
1196	implicit none
1197	!-------------------------------------------------------------------
1198	!
1199	logical , intent(in) :: restart, allowed_to_read
1200	integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
1201	ims, ime, jms, jme, kms, kme, &
1202	its, ite, jts, jte, kts, kte
1203	integer , intent(in) :: p_qi,p_first_scalar
1204	real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &
1205	rublten, &
1206	rvblten, &
1207	rthblten, &
1208	rqvblten, &
1209	rqcblten, &
1210	rqiblten
1211	integer :: i, j, k, itf, jtf, ktf
1212	!
1213	jtf = min0(jte,jde-1)
1214	ktf = min0(kte,kde-1)
1215	itf = min0(ite,ide-1)
1216	!
1217	if(.not.restart)then
1218	do j = jts,jtf
1219	do k = kts,ktf
1220	do i = its,itf
1221	rublten(i,k,j) = 0.
1222	rvblten(i,k,j) = 0.
1223	rthblten(i,k,j) = 0.
1224	rqvblten(i,k,j) = 0.
1225	rqcblten(i,k,j) = 0.
1226	enddo
1227	enddo
1228	enddo
1229	endif
1230	!
1231	if (p_qi .ge. p_first_scalar .and. .not.restart) then
1232	do j = jts,jtf
1233	do k = kts,ktf
1234	do i = its,itf
1235	rqiblten(i,k,j) = 0.
1236	enddo
1237	enddo
1238	enddo
1239	endif
1240	!
1241	end subroutine ysuinit
1242	!-------------------------------------------------------------------
1243	end module module_bl_ysu

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

Download in other formats:

Original Format