Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_ysu.F @ 3553

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

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

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

Download in other formats:

Original Format