MODULE module_sf_bep
USE module_sf_urban
integer nurbm
parameter (nurbm=3)
integer ndm
parameter (ndm=2)
integer nz_um
parameter(nz_um=13)
integer ng_u
parameter (ng_u=10)
integer nwr_u
parameter (nwr_u=10)
real dz_u
parameter (dz_u=5.)
real vk
real g_u
real pi
real r
real cp_u
real rcp_u
real sigma
real p0
real cdrag
parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)
parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4)
CONTAINS
subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, &
th_phy,rho,p_phy,swdown,glw, &
gmt,julday,xlong,xlat, &
declin_urb,cosz_urb2d,omg_urb2d, &
num_urban_layers, &
trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
a_u,a_v,a_t,a_e,b_u,b_v, &
b_t,b_e,dlg,dl_u,sf,vl, &
rl_up,rs_abs,emiss,grdflx_urb, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
implicit none
INTEGER :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
itimestep
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: DZ8W
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: P_PHY
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: RHO
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: TH_PHY
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: T_PHY
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U_PHY
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V_PHY
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V
REAL, DIMENSION( ims:ime , jms:jme ) :: GLW
REAL, DIMENSION( ims:ime , jms:jme ) :: swdown
REAL, DIMENSION( ims:ime, jms:jme ) :: UST
INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: UTYPE_URB2D
REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: FRC_URB2D
REAL, INTENT(IN ) :: GMT
INTEGER, INTENT(IN ) :: JULDAY
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN ) :: XLAT, XLONG
REAL, INTENT(IN) :: DECLIN_URB
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
INTEGER, INTENT(IN ) :: num_urban_layers
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
real z(ims:ime,kms:kme,jms:jme)
REAL, INTENT(IN ):: DT
real a_u(ims:ime,kms:kme,jms:jme)
real a_v(ims:ime,kms:kme,jms:jme)
real a_t(ims:ime,kms:kme,jms:jme)
real a_e(ims:ime,kms:kme,jms:jme)
real b_u(ims:ime,kms:kme,jms:jme)
real b_v(ims:ime,kms:kme,jms:jme)
real b_t(ims:ime,kms:kme,jms:jme)
real b_e(ims:ime,kms:kme,jms:jme)
real dlg(ims:ime,kms:kme,jms:jme)
real dl_u(ims:ime,kms:kme,jms:jme)
real sf(ims:ime,kms:kme,jms:jme)
real vl(ims:ime,kms:kme,jms:jme)
real rl_up(ims:ime,jms:jme)
real rs_abs(ims:ime,jms:jme)
real emiss(ims:ime,jms:jme)
real grdflx_urb(ims:ime,jms:jme)
real alag_u(nurbm)
real alaw_u(nurbm)
real alar_u(nurbm)
real csg_u(nurbm)
real csw_u(nurbm)
real csr_u(nurbm)
real twini_u(nurbm)
real trini_u(nurbm)
real tgini_u(nurbm)
real albg_u(nurbm)
real albw_u(nurbm)
real albr_u(nurbm)
real emg_u(nurbm)
real emw_u(nurbm)
real emr_u(nurbm)
real fww(nz_um,nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fgw(nz_um,ndm,nurbm)
real fsw(nz_um,ndm,nurbm)
real fws(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
real z0g_u(nurbm)
real z0r_u(nurbm)
integer nd_u(nurbm)
real strd_u(ndm,nurbm)
real drst_u(ndm,nurbm)
real ws_u(ndm,nurbm)
real bs_u(ndm,nurbm)
real h_b(nz_um,nurbm)
real d_b(nz_um,nurbm)
real ss_u(nz_um,nurbm)
real pb_u(nz_um,nurbm)
integer nz_u(nurbm)
real z_u(nz_um)
real z1D(kms:kme)
real ua1D(kms:kme)
real va1D(kms:kme)
real pt1D(kms:kme)
real da1D(kms:kme)
real pr1D(kms:kme)
real pt01D(kms:kme)
real zr1D
real deltar1D
real ah1D
real rs1D
real rld1D
real tw1D(2*ndm,nz_um,nwr_u)
real tg1D(ndm,ng_u)
real tr1D(ndm,nz_um,nwr_u)
real sfw1D(2*ndm,nz_um)
real sfg1D(ndm)
real sfr1D(ndm,nz_um)
real sf1D(kms:kme)
real vl1D(kms:kme)
real a_u1D(kms:kme)
real a_v1D(kms:kme)
real a_t1D(kms:kme)
real a_e1D(kms:kme)
real b_u1D(kms:kme)
real b_v1D(kms:kme)
real b_t1D(kms:kme)
real b_e1D(kms:kme)
real dlg1D(kms:kme)
real dl_u1D(kms:kme)
real tsk1D
real time_bep
integer ind_zwd(nz_um,nwr_u,ndm)
integer ind_gd(ng_u,ndm)
integer ind_zd(nz_um,ndm)
integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
integer it, nint
integer iii
real time_h,tempo,shtot
logical first
character(len=80) :: text
data first/.true./
save first,time_bep
save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,fww,fwg,fgw,fsw,fws,fsg, &
z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
nz_u,z_u
if(num_urban_layers.lt.nz_um*ndm*nwr_u)then
write(*,*)'num_urban_layers too small, please increase to at least ', nz_um*ndm*nwr_u
stop
endif
iii=0
do iz_u=1,nz_um
do iw=1,nwr_u
do id=1,ndm
iii=iii+1
ind_zwd(iz_u,iw,id)=iii
enddo
enddo
enddo
iii=0
do ig=1,ng_u
do id=1,ndm
iii=iii+1
ind_gd(ig,id)=iii
enddo
enddo
iii=0
do iz_u=1,nz_um
do id=1,ndm
iii=iii+1
ind_zd(iz_u,id)=iii
enddo
enddo
do ix=its,ite
do iy=jts,jte
z(ix,kts,iy)=0.
do iz=kts+1,kte+1
z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
enddo
enddo
enddo
if (first) then
call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
call icBEP(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
fww,fwg,fgw,fsw,fws,fsg, &
z0g_u,z0r_u, &
nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
nz_u,z_u, &
twini_u,trini_u)
first=.false.
endif
do ix=its,ite
do iy=jts,jte
if (FRC_URB2D(ix,iy).gt.0.) then
iurb=UTYPE_URB2D(ix,iy)
do iz= kts,kte
ua1D(iz)=u_phy(ix,iz,iy)
va1D(iz)=v_phy(ix,iz,iy)
pt1D(iz)=th_phy(ix,iz,iy)
da1D(iz)=rho(ix,iz,iy)
pr1D(iz)=p_phy(ix,iz,iy)
pt01D(iz)=300.
z1D(iz)=z(ix,iz,iy)
a_u1D(iz)=0.
a_v1D(iz)=0.
a_t1D(iz)=0.
a_e1D(iz)=0.
b_u1D(iz)=0.
b_v1D(iz)=0.
b_t1D(iz)=0.
b_e1D(iz)=0.
enddo
z1D(kte+1)=z(ix,kte+1,iy)
do id=1,ndm
do iz_u=1,nz_um
do iw=1,nwr_u
if(ind_zwd(iz_u,iw,id).gt.num_urban_layers)write(*,*)'ind_zwd too big w',ind_zwd(iz_u,iw,id)
tw1D(2*id-1,iz_u,iw)=tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
tw1D(2*id,iz_u,iw)=tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
enddo
enddo
enddo
do id=1,ndm
do ig=1,ng_u
tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
enddo
do iz_u=1,nz_um
do ir=1,nwr_u
if(ind_zwd(iz_u,ir,id).gt.num_urban_layers)write(*,*)'ind_zwd too big r',ind_zwd(iz_u,ir,id)
tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)
enddo
enddo
enddo
do id=1,ndm
do iz=1,nz_um
sfw1D(2*id-1,iz)=sfw1_urb3d(ix,ind_zd(iz,id),iy)
sfw1D(2*id,iz)=sfw2_urb3d(ix,ind_zd(iz,id),iy)
enddo
enddo
do id=1,ndm
sfg1D(id)=sfg_urb3d(ix,id,iy)
enddo
do id=1,ndm
do iz=1,nz_um
sfr1D(id,iz)=sfr_urb3d(ix,ind_zd(iz,id),iy)
enddo
enddo
rs1D=swdown(ix,iy)
rld1D=glw(ix,iy)
time_h=(itimestep*dt)/3600.+gmt
zr1D=acos(COSZ_URB2D(ix,iy))
deltar1D=DECLIN_URB
ah1D=OMG_URB2D(ix,iy)
call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
zr1D,deltar1D,ah1D,rs1D,rld1D, &
alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
fww,fwg,fgw,fsw,fws,fsg, &
z0g_u,z0r_u, &
nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
nz_u,z_u, &
tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D, &
a_u1D,a_v1D,a_t1D,a_e1D, &
b_u1D,b_v1D,b_t1D,b_e1D, &
dlg1D,dl_u1D,tsk1D,sf1D,vl1D,rl_up(ix,iy), &
rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy))
do id=1,ndm
do iz=1,nz_um
sfw1_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id-1,iz)
sfw2_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id,iz)
enddo
enddo
do id=1,ndm
sfg_urb3d(ix,id,iy)=sfg1D(id)
enddo
do id=1,ndm
do iz=1,nz_um
sfr_urb3d(ix,ind_zd(iz,id),iy)=sfr1D(id,iz)
enddo
enddo
do id=1,ndm
do iz_u=1,nz_um
do iw=1,nwr_u
tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw)
tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw)
enddo
enddo
enddo
do id=1,ndm
do ig=1,ng_u
tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
enddo
do iz_u=1,nz_um
do ir=1,nwr_u
trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
enddo
enddo
enddo
do iz= kts,kte
sf(ix,iz,iy)=sf1D(iz)
vl(ix,iz,iy)=vl1D(iz)
a_u(ix,iz,iy)=a_u1D(iz)
a_v(ix,iz,iy)=a_v1D(iz)
a_t(ix,iz,iy)=a_t1D(iz)
a_e(ix,iz,iy)=a_e1D(iz)
b_u(ix,iz,iy)=b_u1D(iz)
b_v(ix,iz,iy)=b_v1D(iz)
b_t(ix,iz,iy)=b_t1D(iz)
b_e(ix,iz,iy)=b_e1D(iz)
dlg(ix,iz,iy)=dlg1D(iz)
dl_u(ix,iz,iy)=dl_u1D(iz)
enddo
sf(ix,kte+1,iy)=sf1D(kte+1)
endif
enddo
enddo
time_bep=time_bep+dt
return
end subroutine BEP
subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, &
zr,deltar,ah,rs,rld, &
alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
fww,fwg,fgw,fsw,fws,fsg, &
z0g_u,z0r_u, &
nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
nz_u,z_u, &
tw,tg,tr,sfw,sfg,sfr, &
a_u,a_v,a_t,a_e, &
b_u,b_v,b_t,b_e, &
dlg,dl_u,tsk,sf,vl,rl_up,rs_abs,emiss,grdflx_urb)
implicit none
integer kms,kme,kts,kte
real z(kms:kme)
real ua(kms:kme)
real va(kms:kme)
real pt(kms:kme)
real da(kms:kme)
real pr(kms:kme)
real pt0(kms:kme)
real dt
real zr
real deltar
real ah
real rs
real rld
integer iurb
real alag_u(nurbm)
real alaw_u(nurbm)
real alar_u(nurbm)
real csg_u(nurbm)
real csw_u(nurbm)
real csr_u(nurbm)
real albg_u(nurbm)
real albw_u(nurbm)
real albr_u(nurbm)
real emg_u(nurbm)
real emw_u(nurbm)
real emr_u(nurbm)
real fww(nz_um,nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fgw(nz_um,ndm,nurbm)
real fsw(nz_um,ndm,nurbm)
real fws(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
real z0g_u(nurbm)
real z0r_u(nurbm)
integer nd_u(nurbm)
real strd_u(ndm,nurbm)
real drst_u(ndm,nurbm)
real ws_u(ndm,nurbm)
real bs_u(ndm,nurbm)
real h_b(nz_um,nurbm)
real d_b(nz_um,nurbm)
real ss_u(nz_um,nurbm)
real pb_u(nz_um,nurbm)
integer nz_u(nurbm)
real z_u(nz_um)
real tw(2*ndm,nz_um,nwr_u)
real tr(ndm,nz_um,nwr_u)
real tg(ndm,ng_u)
real sfw(2*ndm,nz_um)
real sfg(ndm)
real sfr(ndm,nz_um)
real gfg(ndm)
real gfr(ndm,nz_um)
real gfw(2*ndm,nz_um)
real sf(kms:kme)
real vl(kms:kme)
real a_u(kms:kme)
real a_v(kms:kme)
real a_t(kms:kme)
real a_e(kms:kme)
real b_u(kms:kme)
real b_v(kms:kme)
real b_t(kms:kme)
real b_e(kms:kme)
real dlg(kms:kme)
real dl_u(kms:kme)
real tsk
real dz(kms:kme)
real ua_u(nz_um)
real va_u(nz_um)
real pt_u(nz_um)
real da_u(nz_um)
real pt0_u(nz_um)
real pr_u(nz_um)
integer nd
real alag(ng_u)
real alar(nwr_u)
real alaw(nwr_u)
real csg(ng_u)
real csr(nwr_u)
real csw(nwr_u)
real z0(ndm,nz_um)
real ws(ndm)
real bs(ndm)
real strd(ndm)
real drst(ndm)
real ss(nz_um)
real pb(nz_um)
real rsg(ndm)
real rsw(2*ndm,nz_um)
real rlg(ndm)
real rlw(2*ndm,nz_um)
real ptg(ndm)
real ptr(ndm,nz_um)
real ptw(2*ndm,nz_um)
real uhb_u(ndm,nz_um)
real uva_u(2*ndm,nz_um)
real uvb_u(2*ndm,nz_um)
real vhb_u(ndm,nz_um)
real vva_u(2*ndm,nz_um)
real vvb_u(2*ndm,nz_um)
real thb_u(ndm,nz_um)
real tva_u(2*ndm,nz_um)
real tvb_u(2*ndm,nz_um)
real ehb_u(ndm,nz_um)
real evb_u(2*ndm,nz_um)
real rs_abs
real rl_up
real emiss
real grdflx_urb
real shtot,aaa
real dt_int
integer nt_int
integer iz,id, it_int
integer iwrong,iw,ix,iy
do iz=kts,kte
dz(iz)=z(iz+1)-z(iz)
end do
call param(iurb,nz_u(iurb),nd_u(iurb), &
csg_u,csg,alag_u,alag,csr_u,csr, &
alar_u,alar,csw_u,csw,alaw_u,alaw, &
ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb)
call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,ua,ua_u)
call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,va,va_u)
call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt,pt_u)
call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt0,pt0_u)
call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pr,pr_u)
call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,da,da_u)
call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws, &
drst,strd,ss,pb, &
tw,tg,albg_u(iurb),albw_u(iurb), &
emw_u(iurb),emg_u(iurb), &
fww,fwg,fgw,fsw,fsg, &
zr,deltar,ah, &
rs,rld,rsw,rsg,rlw,rlg)
call upward_rad(nd_u(iurb),iurb,nz_u(iurb),ws,bs,sigma,fsw,fsg,pb,ss, &
tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg, &
tw,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw, &
tr,emr_u(iurb),albr_u(iurb),rld,rs,sfr, &
rs_abs,rl_up,emiss,grdflx_urb)
call surf_temp(nz_u(iurb),nd_u(iurb),pr_u,dt,ss, &
rs,rld,rsg,rlg,rsw,rlw, &
tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg, &
tr,alar,csr,emr_u(iurb),albr_u(iurb),ptr,sfr,gfr, &
tw,alaw,csw,emw_u(iurb),albw_u(iurb),ptw,sfw,gfw)
call buildings(nd_u(iurb),nz_u(iurb),z0,ua_u,va_u, &
pt_u,pt0_u,ptg,ptr,da_u,ptw,drst, &
uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
uhb_u,vhb_u,thb_u,ehb_u,ss,dt)
do id=1,nd_u(iurb)
sfg(id)=-da_u(1)*cp_u*thb_u(id,1)
do iz=2,nz_u(iurb)
sfr(id,iz)=-da_u(iz)*cp_u*thb_u(id,iz)
enddo
do iz=1,nz_u(iurb)
sfw(2*id-1,iz)=-da_u(iz)*cp_u*(tvb_u(2*id-1,iz)+ &
tva_u(2*id-1,iz)*pt_u(iz))
sfw(2*id,iz)=-da_u(iz)*cp_u*(tvb_u(2*id,iz)+ &
tva_u(2*id,iz)*pt_u(iz))
enddo
enddo
call urban_meso(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z,dz,z_u,pb,ss,bs,ws,sf, &
vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
uhb_u,vhb_u,thb_u,ehb_u, &
a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)
call interp_length(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z_u,z,ss,ws,bs,dlg,dl_u)
return
end subroutine BEP1D
subroutine param(iurb,nz,nd, &
csg_u,csg,alag_u,alag,csr_u,csr, &
alar_u,alar,csw_u,csw,alaw_u,alaw, &
ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb)
implicit none
integer iurb
integer nz
integer nd
real alag_u(nurbm)
real alar_u(nurbm)
real alaw_u(nurbm)
real bs_u(ndm,nurbm)
real csg_u(nurbm)
real csr_u(nurbm)
real csw_u(nurbm)
real drst_u(ndm,nurbm)
real strd_u(ndm,nurbm)
real ws_u(ndm,nurbm)
real z0g_u(nurbm)
real z0r_u(nurbm)
real ss_u(nz_um,nurbm)
real pb_u(nz_um,nurbm)
real alag(ng_u)
real alar(nwr_u)
real alaw(nwr_u)
real csg(ng_u)
real csr(nwr_u)
real csw(nwr_u)
real bs(ndm)
real drst(ndm)
real strd(ndm)
real ws(ndm)
real z0(ndm,nz_um)
real ss(nz_um)
real pb(nz_um)
integer id,ig,ir,iw,iz
ss=0.
pb=0.
csg=0.
alag=0.
csr=0.
alar=0.
csw=0.
alaw=0.
z0=0.
ws=0.
bs=0.
strd=0.
drst=0.
do iz=1,nz+1
ss(iz)=ss_u(iz,iurb)
pb(iz)=pb_u(iz,iurb)
end do
do ig=1,ng_u
csg(ig)=csg_u(iurb)
alag(ig)=alag_u(iurb)
enddo
do ir=1,nwr_u
csr(ir)=csr_u(iurb)
alar(ir)=alar_u(iurb)
enddo
do iw=1,nwr_u
csw(iw)=csw_u(iurb)
alaw(iw)=alaw_u(iurb)
enddo
do id=1,nd
z0(id,1)=z0g_u(iurb)
do iz=2,nz+1
z0(id,iz)=z0r_u(iurb)
enddo
enddo
do id=1,nd
ws(id)=ws_u(id,iurb)
bs(id)=bs_u(id,iurb)
strd(id)=strd_u(id,iurb)
drst(id)=drst_u(id,iurb)
enddo
return
end subroutine param
subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
implicit none
integer kts,kte,kms,kme
real z(kms:kme)
real c(kms:kme)
integer nz_u
real z_u(nz_um)
real c_u(nz_um)
integer iz_u,iz
real ctot,dz
do iz_u=1,nz_u
ctot=0.
do iz=kts,kte
dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
ctot=ctot+c(iz)*dz
enddo
c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
enddo
return
end subroutine interpol
subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, &
tw,tg,albg,albw,emw,emg, &
fww,fwg,fgw,fsw,fsg, &
zr,deltar,ah, &
rs,rl,rsw,rsg,rlw,rlg)
implicit none
integer iurb
integer nd
integer nz_u
real z(nz_um)
real ws(ndm)
real drst(ndm)
real strd(ndm)
real ss(nz_um)
real pb(nz_um)
real tw(2*ndm,nz_um,nwr_u)
real tg(ndm,ng_u)
real albg
real albw
real emg
real emw
real fgw(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
real fsw(nz_um,ndm,nurbm)
real fws(nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fww(nz_um,nz_um,ndm,nurbm)
real ah
real zr
real deltar
real rs
real rl
real rlg(ndm)
real rlw(2*ndm,nz_um)
real rsg(ndm)
real rsw(2*ndm,nz_um)
integer id,iz
call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
rs,rsw,rsg)
do id=1,nd
call long_rad(iurb,nz_u,id,emw,emg, &
fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
call short_rad(iurb,nz_u,id,albw,albg,fwg,fww,fgw,rsg,rsw,pb)
enddo
return
end subroutine modif_rad
subroutine surf_temp(nz_u,nd,pr,dt,ss,rs,rl,rsg,rlg,rsw,rlw, &
tg,alag,csg,emg,albg,ptg,sfg,gfg, &
tr,alar,csr,emr,albr,ptr,sfr,gfr, &
tw,alaw,csw,emw,albw,ptw,sfw,gfw)
implicit none
integer nz_u
integer nd
real alag(ng_u)
real alar(nwr_u)
real alaw(nwr_u)
real albg
real albr
real albw
real csg(ng_u)
real csr(nwr_u)
real csw(nwr_u)
real dt
real emg
real emr
real emw
real pr(nz_um)
real rs
real rl
real rlg(ndm)
real rlw(2*ndm,nz_um)
real rsg(ndm)
real rsw(2*ndm,nz_um)
real sfg(ndm)
real sfr(ndm,nz_um)
real sfw(2*ndm,nz_um)
real gfg(ndm)
real gfr(ndm,nz_um)
real gfw(2*ndm,nz_um)
real ss(nz_um)
real tg(ndm,ng_u)
real tr(ndm,nz_um,nwr_u)
real tw(2*ndm,nz_um,nwr_u)
real ptg(ndm)
real ptr(ndm,nz_um)
real ptw(2*ndm,nz_um)
integer id,ig,ir,iw,iz
real rtg(ndm)
real rtr(ndm,nz_um)
real rtw(2*ndm,nz_um)
real tg_tmp(ng_u)
real tr_tmp(nwr_u)
real tw_tmp(nwr_u)
real dzg_u(ng_u)
real dzr_u(nwr_u)
real dzw_u(nwr_u)
data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
data dzr_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
data dzw_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
do id=1,nd
do ig=1,ng_u
tg_tmp(ig)=tg(id,ig)
end do
call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg, &
rsg(id),rlg(id),pr(1), &
dt,emg,albg, &
rtg(id),sfg(id),gfg(id))
do ig=1,ng_u
tg(id,ig)=tg_tmp(ig)
end do
do iz=2,nz_u
if(ss(iz).gt.0.)then
do ir=1,nwr_u
tr_tmp(ir)=tr(id,iz,ir)
end do
call soil_temp(nwr_u,dzr_u,tr_tmp,ptr(id,iz), &
alar,csr,rs,rl,pr(iz),dt,emr,albr, &
rtr(id,iz),sfr(id,iz),gfr(id,iz))
do ir=1,nwr_u
tr(id,iz,ir)=tr_tmp(ir)
end do
end if
end do
do iz=1,nz_u
do iw=1,nwr_u
tw_tmp(iw)=tw(2*id-1,iz,iw)
end do
call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id-1,iz),alaw, &
csw, &
rsw(2*id-1,iz),rlw(2*id-1,iz), &
pr(iz),dt,emw, &
albw,rtw(2*id-1,iz),sfw(2*id-1,iz),gfw(2*id-1,iz))
do iw=1,nwr_u
tw(2*id-1,iz,iw)=tw_tmp(iw)
end do
do iw=1,nwr_u
tw_tmp(iw)=tw(2*id,iz,iw)
end do
call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id,iz),alaw, &
csw, &
rsw(2*id,iz),rlw(2*id,iz), &
pr(iz),dt,emw, &
albw,rtw(2*id,iz),sfw(2*id,iz),gfw(2*id,iz))
do iw=1,nwr_u
tw(2*id,iz,iw)=tw_tmp(iw)
end do
end do
end do
return
end subroutine surf_temp
subroutine buildings(nd,nz,z0,ua_u,va_u,pt_u,pt0_u, &
ptg,ptr,da_u,ptw, &
drst,uva_u,vva_u,uvb_u,vvb_u, &
tva_u,tvb_u,evb_u, &
uhb_u,vhb_u,thb_u,ehb_u,ss,dt)
implicit none
integer nd
integer nz
real ua_u(nz_um)
real va_u(nz_um)
real da_u(nz_um)
real drst(ndm)
real dz
real pt_u(nz_um)
real pt0_u(nz_um)
real ptg(ndm)
real ptr(ndm,nz_um)
real ptw(2*ndm,nz_um)
real ss(nz_um)
real z0(ndm,nz_um)
real dt
real uhb_u(ndm,nz_um)
real uva_u(2*ndm,nz_um)
real uvb_u(2*ndm,nz_um)
real vhb_u(ndm,nz_um)
real vva_u(2*ndm,nz_um)
real vvb_u(2*ndm,nz_um)
real thb_u(ndm,nz_um)
real tva_u(2*ndm,nz_um)
real tvb_u(2*ndm,nz_um)
real ehb_u(ndm,nz_um)
real evb_u(2*ndm,nz_um)
integer id,iz
dz=dz_u
do id=1,nd
call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1), &
ptg(id),uhb_u(id,1), &
vhb_u(id,1),thb_u(id,1),ehb_u(id,1))
do iz=2,nz
if(ss(iz).gt.0)then
call flux_flat(dz,z0(id,iz),ua_u(iz), &
va_u(iz),pt_u(iz),pt0_u(iz), &
ptr(id,iz),uhb_u(id,iz), &
vhb_u(id,iz),thb_u(id,iz),ehb_u(id,iz))
else
uhb_u(id,iz) = 0.0
vhb_u(id,iz) = 0.0
thb_u(id,iz) = 0.0
ehb_u(id,iz) = 0.0
endif
end do
do iz=1,nz
call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
ptw(2*id-1,iz), &
uva_u(2*id-1,iz),vva_u(2*id-1,iz), &
uvb_u(2*id-1,iz),vvb_u(2*id-1,iz), &
tva_u(2*id-1,iz),tvb_u(2*id-1,iz), &
evb_u(2*id-1,iz),drst(id),dt)
call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
ptw(2*id,iz), &
uva_u(2*id,iz),vva_u(2*id,iz), &
uvb_u(2*id,iz),vvb_u(2*id,iz), &
tva_u(2*id,iz),tvb_u(2*id,iz), &
evb_u(2*id,iz),drst(id),dt)
end do
end do
return
end subroutine buildings
subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl, &
uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
uhb_u,vhb_u,thb_u,ehb_u, &
a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)
implicit none
integer kms,kme,kts,kte
real z(kms:kme)
real dz(kms:kme)
integer nz_u
integer nd
real bs(ndm)
real ws(ndm)
real z_u(nz_um)
real pb(nz_um)
real ss(nz_um)
real uhb_u(ndm,nz_um)
real uva_u(2*ndm,nz_um)
real uvb_u(2*ndm,nz_um)
real vhb_u(ndm,nz_um)
real vva_u(2*ndm,nz_um)
real vvb_u(2*ndm,nz_um)
real thb_u(ndm,nz_um)
real tva_u(2*ndm,nz_um)
real tvb_u(2*ndm,nz_um)
real ehb_u(ndm,nz_um)
real evb_u(2*ndm,nz_um)
real sf(kms:kme)
real vl(kms:kme)
real a_u(kms:kme)
real a_v(kms:kme)
real a_t(kms:kme)
real a_e(kms:kme)
real b_u(kms:kme)
real b_v(kms:kme)
real b_t(kms:kme)
real b_e(kms:kme)
real dzz
real fact
integer id,iz,iz_u
real se,sr,st,su,sv
real uet(kms:kme)
real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb
do iz=kts,kte
a_u(iz)=0.
a_v(iz)=0.
a_t(iz)=0.
a_e(iz)=0.
b_u(iz)=0.
b_v(iz)=0.
b_e(iz)=0.
b_t(iz)=0.
uet(iz)=0.
end do
do iz=kts,kte
sf(iz)=0.
vl(iz)=0.
enddo
sf(kte+1)=0.
do id=1,nd
do iz=kts+1,kte+1
sr=0.
do iz_u=2,nz_u
if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
sr=pb(iz_u)
endif
enddo
sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
enddo
enddo
do iz=kts,kte
do id=1,nd
vtot=0.
do iz_u=1,nz_u
dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
vtot=vtot+pb(iz_u+1)*dzz
enddo
vtot=vtot/(z(iz+1)-z(iz))
vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
enddo
enddo
do id=1,nd
fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
b_t(kts)=b_t(kts)+thb_u(id,1)*fact
b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
b_v(kts)=b_v(kts)+vhb_u(id,1)*fact
b_e(kts)=b_e(kts)+ehb_u(id,1)*fact*(z_u(2)-z_u(1))
do iz=kts,kte
st=0.
su=0.
sv=0.
se=0.
do iz_u=2,nz_u
if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
st=st+ss(iz_u)*thb_u(id,iz_u)
su=su+ss(iz_u)*uhb_u(id,iz_u)
sv=sv+ss(iz_u)*vhb_u(id,iz_u)
se=se+ss(iz_u)*ehb_u(id,iz_u)*(z_u(iz_u+1)-z_u(iz_u))
endif
enddo
fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
b_t(iz)=b_t(iz)+st*fact
b_u(iz)=b_u(iz)+su*fact
b_v(iz)=b_v(iz)+sv*fact
b_e(iz)=b_e(iz)+se*fact
enddo
enddo
do iz=kts,kte
uet(iz)=0.
do id=1,nd
vtb=0.
vta=0.
vua=0.
vub=0.
vva=0.
vvb=0.
veb=0.
vte=0.
do iz_u=1,nz_u
dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
fact=dzz/(ws(id)+bs(id))
vtb=vtb+pb(iz_u+1)* &
(tvb_u(2*id-1,iz_u)+tvb_u(2*id,iz_u))*fact
vta=vta+pb(iz_u+1)* &
(tva_u(2*id-1,iz_u)+tva_u(2*id,iz_u))*fact
vua=vua+pb(iz_u+1)* &
(uva_u(2*id-1,iz_u)+uva_u(2*id,iz_u))*fact
vva=vva+pb(iz_u+1)* &
(vva_u(2*id-1,iz_u)+vva_u(2*id,iz_u))*fact
vub=vub+pb(iz_u+1)* &
(uvb_u(2*id-1,iz_u)+uvb_u(2*id,iz_u))*fact
vvb=vvb+pb(iz_u+1)* &
(vvb_u(2*id-1,iz_u)+vvb_u(2*id,iz_u))*fact
veb=veb+pb(iz_u+1)* &
(evb_u(2*id-1,iz_u)+evb_u(2*id,iz_u))*fact
enddo
fact=1./vl(iz)/dz(iz)/nd
b_t(iz)=b_t(iz)+vtb*fact
a_t(iz)=a_t(iz)+vta*fact
a_u(iz)=a_u(iz)+vua*fact
a_v(iz)=a_v(iz)+vva*fact
b_u(iz)=b_u(iz)+vub*fact
b_v(iz)=b_v(iz)+vvb*fact
b_e(iz)=b_e(iz)+veb*fact
uet(iz)=uet(iz)+vte*fact
enddo
enddo
return
end subroutine urban_meso
subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
dlg,dl_u)
implicit none
integer kms,kme,kts,kte
real z(kms:kme)
integer nd
integer nz_u
real z_u(nz_um)
real bs(ndm)
real ss(nz_um)
real ws(ndm)
real dlg(kms:kme)
real dl_u(kms:kme)
real dlgtmp
integer id,iz,iz_u
real sftot
real ulu,ssl
do iz=kts,kte
ulu=0.
ssl=0.
do id=1,nd
do iz_u=2,nz_u
if(z_u(iz_u).gt.z(iz))then
ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
ssl=ssl+ss(iz_u)/nd
endif
enddo
enddo
if(ulu.ne.0)then
dl_u(iz)=ssl/ulu
else
dl_u(iz)=0.
endif
enddo
do iz=kts,kte
dlg(iz)=0.
do id=1,nd
sftot=ws(id)
dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
do iz_u=1,nz_u
if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
dlgtmp=dlgtmp+ss(iz_u)*bs(id)/ &
((z(iz)+z(iz+1))/2.-z_u(iz_u))
sftot=sftot+ss(iz_u)*bs(id)
endif
enddo
dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
enddo
dlg(iz)=1./dlg(iz)
enddo
return
end subroutine interp_length
subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
rs,rsw,rsg)
implicit none
integer nd
integer nz_u
real ah
real deltar
real drst(ndm)
real rs
real ss(nz_um)
real pb(nz_um)
real ws(ndm)
real z(nz_um)
real zr
real rsg(ndm)
real rsw(2*ndm,nz_um)
integer id,iz,jz
real aae,aaw,bbb,phix,rd,rtot,wsd
if(rs.eq.0.or.sin(zr).eq.1)then
do id=1,nd
rsg(id)=0.
do iz=1,nz_u
rsw(2*id-1,iz)=0.
rsw(2*id,iz)=0.
enddo
enddo
else
if(abs(sin(zr)).gt.1.e-10)then
if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
bbb=pi/2.
elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
bbb=-pi/2.
else
bbb=asin(cos(deltar)*sin(ah)/sin(zr))
endif
else
if(cos(deltar)*sin(ah).ge.0)then
bbb=pi/2.
elseif(cos(deltar)*sin(ah).lt.0)then
bbb=-pi/2.
endif
endif
phix=zr
do id=1,nd
rsg(id)=0.
aae=bbb-drst(id)
aaw=bbb-drst(id)+pi
do iz=1,nz_u
rsw(2*id-1,iz)=0.
rsw(2*id,iz)=0.
if(pb(iz+1).gt.0.)then
do jz=1,nz_u
if(abs(sin(aae)).gt.1.e-10)then
call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae, &
ws(id),rd)
rsw(2*id-1,iz)=rsw(2*id-1,iz)+rs*rd*ss(jz+1)/pb(iz+1)
endif
if(abs(sin(aaw)).gt.1.e-10)then
call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw, &
ws(id),rd)
rsw(2*id,iz)=rsw(2*id,iz)+rs*rd*ss(jz+1)/pb(iz+1)
endif
enddo
endif
enddo
if(abs(sin(aae)).gt.1.e-10)then
wsd=abs(ws(id)/sin(aae))
do jz=1,nz_u
rd=max(0.,wsd-z(jz+1)*tan(phix))
rsg(id)=rsg(id)+rs*rd*ss(jz+1)/wsd
enddo
rtot=0.
do iz=1,nz_u
rtot=rtot+(rsw(2*id,iz)+rsw(2*id-1,iz))* &
(z(iz+1)-z(iz))
enddo
rtot=rtot+rsg(id)*ws(id)
else
rsg(id)=rs
endif
enddo
endif
return
end subroutine shadow_mas
subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
implicit none
real aa
real hu
real phix
real ws
real z1
real z2
real rd
real x1,x2
x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
x2=max((hu-z2)*tan(phix),0.)
rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
return
end subroutine shade_wall
subroutine long_rad(iurb,nz_u,id,emw,emg, &
fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
implicit none
real emg
real emw
real fgw(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
real fsw(nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fww(nz_um,nz_um,ndm,nurbm)
integer id
integer iurb
integer nz_u
real pb(nz_um)
real rl
real tg(ndm,ng_u)
real tw(2*ndm,nz_um,nwr_u)
real rlg(ndm)
real rlw(2*ndm,nz_um)
integer i,j
real aaa(2*nz_um+1,2*nz_um+1)
real bbb(2*nz_um+1)
do i=1,nz_u
do j=1,nz_u
aaa(i,j)=0.
enddo
aaa(i,i)=1.
do j=nz_u+1,2*nz_u
aaa(i,j)=-(1.-emw)*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
enddo
aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg(id,ng_u)**4
do j=1,nz_u
bbb(i)=bbb(i)+pb(j+1)*emw*sigma*fww(j,i,id,iurb)* &
tw(2*id,j,nwr_u)**4+ &
fww(j,i,id,iurb)*rl*(1.-pb(j+1))
enddo
enddo
do i=1+nz_u,2*nz_u
do j=1,nz_u
aaa(i,j)=-(1.-emw)*fww(j,i-nz_u,id,iurb)*pb(j+1)
enddo
do j=1+nz_u,2*nz_u
aaa(i,j)=0.
enddo
aaa(i,i)=1.
aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
bbb(i)=fsw(i-nz_u,id,iurb)*rl+ &
emg*fgw(i-nz_u,id,iurb)*sigma*tg(id,ng_u)**4
do j=1,nz_u
bbb(i)=bbb(i)+pb(j+1)*emw*sigma*fww(j,i-nz_u,id,iurb)* &
tw(2*id-1,j,nwr_u)**4+ &
fww(j,i-nz_u,id,iurb)*rl*(1.-pb(j+1))
enddo
enddo
do j=1,nz_u
aaa(2*nz_u+1,j)=-(1.-emw)*fwg(j,id,iurb)*pb(j+1)
enddo
do j=nz_u+1,2*nz_u
aaa(2*nz_u+1,j)=-(1.-emw)*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
enddo
aaa(2*nz_u+1,2*nz_u+1)=1.
bbb(2*nz_u+1)=fsg(id,iurb)*rl
do i=1,nz_u
bbb(2*nz_u+1)=bbb(2*nz_u+1)+emw*sigma*fwg(i,id,iurb)*pb(i+1)* &
(tw(2*id-1,i,nwr_u)**4+tw(2*id,i,nwr_u)**4)+ &
2.*fwg(i,id,iurb)*(1.-pb(i+1))*rl
enddo
call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
do i=1,nz_u
rlw(2*id-1,i)=bbb(i)
enddo
do i=nz_u+1,2*nz_u
rlw(2*id,i-nz_u)=bbb(i)
enddo
rlg(id)=bbb(2*nz_u+1)
return
end subroutine long_rad
subroutine short_rad(iurb,nz_u,id,albw, &
albg,fwg,fww,fgw,rsg,rsw,pb)
implicit none
real albg
real albw
real fgw(nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fww(nz_um,nz_um,ndm,nurbm)
integer id
integer iurb
integer nz_u
real pb(nz_um)
real rsg(ndm)
real rsw(2*ndm,nz_um)
integer i,j
real aaa(2*nz_um+1,2*nz_um+1)
real bbb(2*nz_um+1)
do i=1,nz_u
do j=1,nz_u
aaa(i,j)=0.
enddo
aaa(i,i)=1.
do j=nz_u+1,2*nz_u
aaa(i,j)=-albw*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
enddo
aaa(i,2*nz_u+1)=-albg*fgw(i,id,iurb)
bbb(i)=rsw(2*id-1,i)
enddo
do i=1+nz_u,2*nz_u
do j=1,nz_u
aaa(i,j)=-albw*fww(j,i-nz_u,id,iurb)*pb(j+1)
enddo
do j=1+nz_u,2*nz_u
aaa(i,j)=0.
enddo
aaa(i,i)=1.
aaa(i,2*nz_u+1)=-albg*fgw(i-nz_u,id,iurb)
bbb(i)=rsw(2*id,i-nz_u)
enddo
do j=1,nz_u
aaa(2*nz_u+1,j)=-albw*fwg(j,id,iurb)*pb(j+1)
enddo
do j=nz_u+1,2*nz_u
aaa(2*nz_u+1,j)=-albw*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
enddo
aaa(2*nz_u+1,2*nz_u+1)=1.
bbb(2*nz_u+1)=rsg(id)
call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
do i=1,nz_u
rsw(2*id-1,i)=bbb(i)
enddo
do i=nz_u+1,2*nz_u
rsw(2*id,i-nz_u)=bbb(i)
enddo
rsg(id)=bbb(2*nz_u+1)
return
end subroutine short_rad
subroutine gaussj(a,n,b,np)
implicit none
integer np
real a(np,np)
real b(np)
integer nmax
parameter (nmax=150)
real big,dum
integer i,icol,irow
integer j,k,l,ll,n
integer ipiv(nmax)
real pivinv
do j=1,n
ipiv(j)=0.
enddo
do i=1,n
big=0.
do j=1,n
if(ipiv(j).ne.1)then
do k=1,n
if(ipiv(k).eq.0)then
if(abs(a(j,k)).ge.big)then
big=abs(a(j,k))
irow=j
icol=k
endif
elseif(ipiv(k).gt.1)then
pause 'singular matrix in gaussj'
endif
enddo
endif
enddo
ipiv(icol)=ipiv(icol)+1
if(irow.ne.icol)then
do l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
enddo
dum=b(irow)
b(irow)=b(icol)
b(icol)=dum
endif
if(a(icol,icol).eq.0)pause 'singular matrix in gaussj'
pivinv=1./a(icol,icol)
a(icol,icol)=1
do l=1,n
a(icol,l)=a(icol,l)*pivinv
enddo
b(icol)=b(icol)*pivinv
do ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.
do l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
enddo
b(ll)=b(ll)-b(icol)*dum
endif
enddo
enddo
return
end subroutine gaussj
subroutine soil_temp(nz,dz,temp,pt,ala,cs, &
rs,rl,press,dt,em,alb,rt,sf,gf)
implicit none
integer nz
real ala(nz)
real alb
real cs(nz)
real dt
real em
real press
real rl
real rs
real sf
real temp(nz)
real dz(nz)
real gf
real pt
real rt
integer iz
real a(nz,3)
real alpha
real c(nz)
real cddz(nz+2)
real tsig
tsig=temp(nz)
alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
cddz(1)=ala(1)/dz(1)
do iz=2,nz
cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
enddo
a(1,1)=0.
a(1,2)=1.
a(1,3)=0.
c(1)=temp(1)
do iz=2,nz-1
a(iz,1)=-cddz(iz)*dt/dz(iz)
a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
a(iz,3)=-cddz(iz+1)*dt/dz(iz)
c(iz)=temp(iz)
enddo
a(nz,1)=-dt*cddz(nz)/dz(nz)
a(nz,2)=1.+dt*cddz(nz)/dz(nz)
a(nz,3)=0.
c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
call invert(nz,a,c,temp)
pt=temp(nz)*(press/1.e+5)**(-rcp_u)
rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)
gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
return
end subroutine soil_temp
subroutine invert(n,a,c,x)
implicit none
integer n
real a(n,3)
real c(n)
real x(n)
integer i
do i=n-1,1,-1
c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
enddo
do i=2,n
c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
enddo
do i=1,n
x(i)=c(i)/a(i,2)
enddo
return
end subroutine invert
subroutine flux_wall(ua,va,pt,da,ptw,uva,vva,uvb,vvb, &
tva,tvb,evb,drst,dt)
implicit none
real drst
real da
real pt
real ptw
real ua
real va
real dt
real uva
real uvb
real vva
real vvb
real tva
real tvb
real evb
real hc
real u_ort
real vett
vett=(ua**2+va**2)**.5
u_ort=abs((cos(drst)*ua-sin(drst)*va))
uva=-cdrag*u_ort/2.*cos(drst)*cos(drst)
vva=-cdrag*u_ort/2.*sin(drst)*sin(drst)
uvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*va
vvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*ua
hc=5.678*(1.09+0.23*(vett/0.3048))
if(hc.gt.da*cp_u/dt)then
hc=da*cp_u/dt
endif
tvb=hc*ptw/da/cp_u-hc/da/cp_u*pt
tva = 0.
evb=cdrag*(abs(u_ort)**3.)/2.
return
end subroutine flux_wall
subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg, &
uhb,vhb,thb,ehb)
implicit none
real dz
real pt
real pt0
real ptg
real ua
real va
real z0
real uhb
real vhb
real thb
real tva
real tvb
real ehb
real aa
real al
real buu
real c
real fbuw
real fbpt
real fh
real fm
real ric
real tstar
real ustar
real utot
real wstar
real zz
real b,cm,ch,rr,tol
parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
utot=(ua**2+va**2)**.5
zz=dz/2.
utot=max(utot,0.01)
ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
aa=vk/log(zz/z0)
if(ric.gt.0)then
fm=1/(1+0.5*b*ric)**2
fh=fm
else
c=b*cm*aa*aa*(zz/z0)**.5
fm=1-b*ric/(1+c*(-ric)**.5)
c=c*ch/cm
fh=1-b*ric/(1+c*(-ric)**.5)
endif
fbuw=-aa*aa*utot*utot*fm
fbpt=-aa*aa*utot*(pt-ptg)*fh/rr
ustar=(-fbuw)**.5
tstar=-fbpt/ustar
al=(vk*g_u*tstar)/(pt*ustar*ustar)
buu=-g_u/pt0*ustar*tstar
uhb=-ustar*ustar*ua/utot
vhb=-ustar*ustar*va/utot
thb=-ustar*tstar
ehb=buu
return
end subroutine flux_flat
subroutine icBEP (alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
fww,fwg,fgw,fsw,fws,fsg, &
z0g_u,z0r_u, &
nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
nz_u,z_u, &
twini_u,trini_u)
implicit none
real alag_u(nurbm)
real alaw_u(nurbm)
real alar_u(nurbm)
real csg_u(nurbm)
real csw_u(nurbm)
real csr_u(nurbm)
real twini_u(nurbm)
real trini_u(nurbm)
real albg_u(nurbm)
real albw_u(nurbm)
real albr_u(nurbm)
real emg_u(nurbm)
real emw_u(nurbm)
real emr_u(nurbm)
real z0g_u(nurbm)
real z0r_u(nurbm)
integer nd_u(nurbm)
real strd_u(ndm,nurbm)
real drst_u(ndm,nurbm)
real ws_u(ndm,nurbm)
real bs_u(ndm,nurbm)
real h_b(nz_um,nurbm)
real d_b(nz_um,nurbm)
real fww(nz_um,nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fgw(nz_um,ndm,nurbm)
real fsw(nz_um,ndm,nurbm)
real fws(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
real ss_u(nz_um,nurbm)
real pb_u(nz_um,nurbm)
integer nz_u(nurbm)
real z_u(nz_um)
integer iz_u,id,ilu,iurb
real dtot
real hbmax
nz_u=0
z_u=0.
ss_u=0.
pb_u=0.
fww=0.
fwg=0.
fgw=0.
fsw=0.
fws=0.
fsg=0.
z_u(1)=0.
do iz_u=1,nz_um-1
z_u(iz_u+1)=z_u(iz_u)+dz_u
enddo
do iurb=1,nurbm
dtot=0.
do ilu=1,nz_um
dtot=dtot+d_b(ilu,iurb)
enddo
do ilu=1,nz_um
d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
enddo
enddo
do iurb=1,nurbm
hbmax=0.
nz_u(iurb)=0
do ilu=1,nz_um
if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
enddo
do iz_u=1,nz_um-1
if(z_u(iz_u+1).gt.hbmax)go to 10
enddo
10 continue
nz_u(iurb)=iz_u+1
do id=1,nd_u(iurb)
call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb), &
z_u,ws_u(id,iurb), &
fww,fwg,fgw,fsg,fsw,fws)
do iz_u=1,nz_u(iurb)
ss_u(iz_u,iurb)=0.
do ilu=1,nz_um
if(z_u(iz_u).le.h_b(ilu,iurb) &
.and.z_u(iz_u+1).gt.h_b(ilu,iurb))then
ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
endif
enddo
enddo
pb_u(1,iurb)=1.
do iz_u=1,nz_u(iurb)
pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
enddo
enddo
end do
return
end subroutine icBEP
subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws)
implicit none
integer iurb
integer nz_u
integer id
real ws
real z(nz_um)
real dxy
real fww(nz_um,nz_um,ndm,nurbm)
real fwg(nz_um,ndm,nurbm)
real fgw(nz_um,ndm,nurbm)
real fsw(nz_um,ndm,nurbm)
real fws(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
integer jz,iz
real hut
real f1,f2,f12,f23,f123,ftot
real fprl,fnrm
real a1,a2,a3,a4,a12,a23,a123
hut=z(nz_u+1)
do jz=1,nz_u
do iz=1,nz_u
call fprls (fprl,dxy,abs(z(jz+1)-z(iz )),ws)
f123=fprl
call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
f23=fprl
call fprls (fprl,dxy,abs(z(jz )-z(iz )),ws)
f12=fprl
call fprls (fprl,dxy,abs(z(jz )-z(iz+1)),ws)
f2 = fprl
a123=dxy*(abs(z(jz+1)-z(iz )))
a12 =dxy*(abs(z(jz )-z(iz )))
a23 =dxy*(abs(z(jz+1)-z(iz+1)))
a1 =dxy*(abs(z(iz+1)-z(iz )))
a2 =dxy*(abs(z(jz )-z(iz+1)))
a3 =dxy*(abs(z(jz+1)-z(jz )))
ftot=0.5*(a123*f123-a23*f23-a12*f12+a2*f2)/a1
fww(iz,jz,id,iurb)=ftot*a1/a3
enddo
call fnrms (fnrm,z(jz+1),dxy,ws)
f12=fnrm
call fnrms (fnrm,z(jz) ,dxy,ws)
f1=fnrm
a1 = ws*dxy
a12= ws*dxy
a4=(z(jz+1)-z(jz))*dxy
ftot=(a12*f12-a12*f1)/a1
fgw(jz,id,iurb)=ftot*a1/a4
call fnrms(fnrm,hut-z(jz) ,dxy,ws)
f12 = fnrm
call fnrms (fnrm,hut-z(jz+1),dxy,ws)
f1 =fnrm
a1 = ws*dxy
a12= ws*dxy
a4 = (z(jz+1)-z(jz))*dxy
ftot=(a12*f12-a12*f1)/a1
fsw(jz,id,iurb)=ftot*a1/a4
enddo
do iz=1,nz_u
call fnrms(fnrm,ws,dxy,hut-z(iz))
f12=fnrm
call fnrms(fnrm,ws,dxy,hut-z(iz+1))
f1=fnrm
a1 = (z(iz+1)-z(iz))*dxy
a2 = (hut-z(iz+1))*dxy
a12= (hut-z(iz))*dxy
a4 = ws*dxy
ftot=(a12*f12-a2*f1)/a1
fws(iz,id,iurb)=ftot*a1/a4
enddo
do iz=1,nz_u
call fnrms (fnrm,ws,dxy,z(iz+1))
f12=fnrm
call fnrms (fnrm,ws,dxy,z(iz ))
f1 =fnrm
a1= (z(iz+1)-z(iz) )*dxy
a2 = z(iz)*dxy
a12= z(iz+1)*dxy
a4 = ws*dxy
ftot=(a12*f12-a2*f1)/a1
fwg(iz,id,iurb)=ftot*a1/a4
enddo
call fprls (fprl,dxy,ws,hut)
fsg(id,iurb)=fprl
return
end subroutine view_factors
SUBROUTINE fprls (fprl,a,b,c)
implicit none
real a,b,c
real x,y
real fprl
x=a/c
y=b/c
if(a.eq.0.or.b.eq.0.)then
fprl=0.
else
fprl=log( ( (1.+x**2)*(1.+y**2)/(1.+x**2+y**2) )**.5)+ &
y*((1.+x**2)**.5)*atan(y/((1.+x**2)**.5))+ &
x*((1.+y**2)**.5)*atan(x/((1.+y**2)**.5))- &
y*atan(y)-x*atan(x)
fprl=fprl*2./(pi*x*y)
endif
return
end subroutine fprls
SUBROUTINE fnrms (fnrm,a,b,c)
implicit none
real a,b,c
real x,y,z,a1,a2,a3,a4,a5,a6
real fnrm
x=a/b
y=c/b
z=x**2+y**2
if(y.eq.0.or.x.eq.0)then
fnrm=0.
else
a1=log( (1.+x*x)*(1.+y*y)/(1.+z) )
a2=y*y*log(y*y*(1.+z)/z/(1.+y*y) )
a3=x*x*log(x*x*(1.+z)/z/(1.+x*x) )
a4=y*atan(1./y)
a5=x*atan(1./x)
a6=sqrt(z)*atan(1./sqrt(z))
fnrm=0.25*(a1+a2+a3)+a4+a5-a6
fnrm=fnrm/(pi*y)
endif
return
end subroutine fnrms
SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
implicit none
integer iurb
real alag_u(nurbm)
real alaw_u(nurbm)
real alar_u(nurbm)
real csg_u(nurbm)
real csw_u(nurbm)
real csr_u(nurbm)
real twini_u(nurbm)
real trini_u(nurbm)
real tgini_u(nurbm)
real albg_u(nurbm)
real albw_u(nurbm)
real albr_u(nurbm)
real emg_u(nurbm)
real emw_u(nurbm)
real emr_u(nurbm)
real z0g_u(nurbm)
real z0r_u(nurbm)
integer nd_u(nurbm)
real strd_u(ndm,nurbm)
real drst_u(ndm,nurbm)
real ws_u(ndm,nurbm)
real bs_u(ndm,nurbm)
real h_b(nz_um,nurbm)
real d_b(nz_um,nurbm)
integer i,iu
integer nurb
h_b=0.
d_b=0.
nurb=ICATE
do iu=1,nurb
nd_u(iu)=0
enddo
csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
do i=1,icate
alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
enddo
twini_u=TBLEND_TBL
trini_u=TRLEND_TBL
tgini_u=TGLEND_TBL
albw_u=ALBB_TBL
albr_u=ALBR_TBL
albg_u=ALBG_TBL
emw_u=EPSB_TBL
emr_u=EPSR_TBL
emg_u=EPSG_TBL
z0r_u=Z0R_TBL
z0g_u=Z0G_TBL
nd_u=NUMDIR_TBL
do iu=1,icate
if(ndm.lt.nd_u(iu))then
write(*,*)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
stop
endif
do i=1,nd_u(iu)
drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
enddo
enddo
do iu=1,ICATE
if(nz_um.lt.numhgt_tbl(iu)+3)then
write(*,*)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
stop
endif
do i=1,NUMHGT_TBL(iu)
h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
enddo
enddo
do i=1,ndm
do iu=1,nurbm
strd_u(i,iu)=100000.
enddo
enddo
return
END SUBROUTINE init_para
subroutine angle(along,alat,day,realt,zr,deltar,ah)
implicit none
real along,alat, realt, zr, deltar, ah, arg
real rad,om,radh,initt, pii, drad, alongt, cphi, sphi
real c1, c2, c3, s1, s2, s3, delta, rmsr2, cd, sid
real et, ahor, chor, coznt
integer day
data rad,om,radh,initt/0.0174533,0.0172142,0.26179939,0/
zr=0.
deltar=0.
ah=0.
pii = 3.14159265358979312
drad = pii/180.
alongt=along/15.
cphi=cos(alat*drad)
sphi=sin(alat*drad)
arg=om*day
c1=cos(arg)
c2=cos(2.*arg)
c3=cos(3.*arg)
s1=sin(arg)
s2=sin(2.*arg)
s3=sin(3.*arg)
delta=0.33281-22.984*c1-0.3499*c2-0.1398*c3+3.7872*s1+0.03205*s2+0.07187*s3
rmsr2=(1./(1.-0.01673*c1))**2
deltar=delta*rad
cd=cos(deltar)
sid=sin(deltar)
et=0.0072*c1-0.0528*c2-0.0012*c3-0.1229*s1-0.1565*s2-0.0041*s3
ahor=realt-12.+et+alongt
ah=ahor*radh
chor=cos(ah)
coznt=sphi*sid+cphi*cd*chor
zr=acos(coznt)
return
END SUBROUTINE angle
subroutine upward_rad(nd_u,iurb,nz_u,ws,bs,sigma,fsw,fsg,pb,ss, &
tg,emg_u,albg_u,rlg,rsg,sfg, &
tw,emw_u,albw_u,rlw,rsw,sfw, &
tr,emr_u,albr_u,rld,rs, sfr, &
rs_abs,rl_up,emiss,grdflx_urb)
implicit none
real rsw(2*ndm,nz_um)
real rlw(2*ndm,nz_um)
real rsg(ndm)
real rlg(ndm)
real rs
real sfw(2*ndm,nz_um)
real sfg(ndm)
real sfr(ndm,nz_um)
real rld
real albg_u
real albw_u
real albr_u
real ws(ndm)
real bs(ndm)
real pb(nz_um)
integer nz_u
real ss(nz_um)
real sigma
real emg_u
real emw_u
real emr_u
real fsw(nz_um,ndm,nurbm)
real fsg(ndm,nurbm)
real tw(2*ndm,nz_um,nwr_u)
real tr(ndm,nz_um,nwr_u)
real tg(ndm,ng_u)
integer iurb
integer id
integer nd_u
real rs_abs
real rl_up
real emiss
real grdflx_urb
integer iz,iw
real rl_inc,rl_emit
real gfl
integer ix,iy,iwrong
iwrong=1
do iz=1,nz_u+1
do id=1,nd_u
do iw=1,nwr_u
if(tr(id,iz,iw).lt.100.)then
write(*,*)'in upward_rad ',iz,id,iw,tr(id,iz,iw)
iwrong=0
endif
enddo
enddo
enddo
if(iwrong.eq.0)stop
rl_up=0.
rs_abs=0.
rl_inc=0.
emiss=0.
rl_emit=0.
grdflx_urb=0.
do id=1,nd_u
rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/nd_u
rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u
rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/nd_u
gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg(id,ng_u)**4.)+sfg(id)
grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u
do iz=2,nz_u
rl_emit=rl_emit-(emr_u*sigma*(tr(id,iz,nwr_u)**4.)+(1-emr_u)*rld)*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
gfl=(1.-albr_u)*rs+emr_u*rld-emr_u*sigma*(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
enddo
do iz=1,nz_u
rl_emit=rl_emit-(emw_u*sigma*( tw(2*id-1,iz,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+ &
(1-emw_u)*( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
rl_inc=rl_inc+(( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
rs_abs=rs_abs+((1.-albw_u)*( rsw(2*id-1,iz)+rsw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) ) &
-emw_u*sigma*( tw(2*id-1,iz,nwr_u)**4.+tw(2*id,iz,nwr_u)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))
grdflx_urb=grdflx_urb-gfl*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
enddo
enddo
emiss=(emg_u+emw_u+emr_u)/3.
rl_up=(rl_inc+rl_emit)-rld
return
END SUBROUTINE upward_rad
END MODULE module_sf_bep