Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
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

Rev	Line
[2759]	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