Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_bl_ysu.F @ 1382

Last change on this file since 1382 was 1, checked in by lfita, 11 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 50.1 KB

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

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

Download in other formats:

Original Format