MODULE module_sf_bep 2
!USE module_model_constants
USE module_sf_urban
! SGClarke 09/11/2008
! Access urban_param.tbl values through calling urban_param_init in module_physics_init
! for CASE (BEPSCHEME) select sf_urban_physics
!
! -----------------------------------------------------------------------
! Dimension for the array used in the BEP module
! -----------------------------------------------------------------------
integer nurbm ! Maximum number of urban classes
parameter (nurbm=3)
integer ndm ! Maximum number of street directions
parameter (ndm=2)
integer nz_um ! Maximum number of vertical levels in the urban grid
parameter(nz_um=18)
integer ng_u ! Number of grid levels in the ground
parameter (ng_u=10)
integer nwr_u ! Number of grid levels in the walls or roofs
parameter (nwr_u=10)
real dz_u ! Urban grid resolution
parameter (dz_u=5.)
! The change of ng_u, nwr_u should be done in agreement with the block data
! in the routine "surf_temp"
! -----------------------------------------------------------------------
! Constant used in the BEP module
! -----------------------------------------------------------------------
real vk ! von Karman constant
real g_u ! Gravity acceleration
real pi !
real r ! Perfect gas constant
real cp_u ! Specific heat at constant pressure
real rcp_u !
real sigma !
real p0 ! Reference pressure at the sea level
real cdrag ! Drag force constant
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, & 1,6
th_phy,rho,p_phy,swdown,glw, &
gmt,julday,xlong,xlat, &
declin_urb,cosz_urb2d,omg_urb2d, &
num_urban_layers,num_urban_hi, &
trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
a_u,a_v,a_t,a_e,b_u,b_v, &
b_t,b_e,b_q,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
!------------------------------------------------------------------------
! Input
!------------------------------------------------------------------------
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
INTEGER, INTENT(IN ) :: num_urban_hi
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, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lp_urb2d
REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lb_urb2d
REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: hgt_urb2d
! integer nx,ny,nz ! Number of points in the mesocsale grid
real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates
REAL, INTENT(IN ):: DT ! Time step
! real zr(ims:ime,jms:jme) ! Solar zenith angle
! real deltar(ims:ime,jms:jme) ! Declination of the sun
! real ah(ims:ime,jms:jme) ! Hour angle
! real rs(ims:ime,jms:jme) ! Solar radiation
!------------------------------------------------------------------------
! Output
!------------------------------------------------------------------------
! real tsk(ims:ime,jms:jme) ! Average of surface temperatures (roads and roofs)
!
! Implicit and explicit components of the source and sink terms at each levels,
! the fluxes can be computed as follow: FX = A*X + B example: t_fluxes = a_t * pt + b_t
real a_u(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in X-direction (center)
real a_v(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in Y-direction (center)
real a_t(ims:ime,kms:kme,jms:jme) ! Implicit component for the temperature
real a_e(ims:ime,kms:kme,jms:jme) ! Implicit component for the TKE
real b_u(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in X-direction (center)
real b_v(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in Y-direction (center)
real b_t(ims:ime,kms:kme,jms:jme) ! Explicit component for the temperature
real b_e(ims:ime,kms:kme,jms:jme) ! Explicit component for the TKE
real b_q(ims:ime,kms:kme,jms:jme) ! Explicit component for the humidity
real dlg(ims:ime,kms:kme,jms:jme) ! Height above ground (L_ground in formula (24) of the BLM paper).
real dl_u(ims:ime,kms:kme,jms:jme) ! Length scale (lb in formula (22) ofthe BLM paper).
! urban surface and volumes
real sf(ims:ime,kms:kme,jms:jme) ! surface of the urban grid cells
real vl(ims:ime,kms:kme,jms:jme) ! volume of the urban grid cells
! urban fluxes
real rl_up(its:ite,jts:jte) ! upward long wave radiation
real rs_abs(its:ite,jts:jte) ! absorbed short wave radiation
real emiss(its:ite,jts:jte) ! emissivity averaged for urban surfaces
real grdflx_urb(its:ite,jts:jte) ! ground heat flux for urban areas
!------------------------------------------------------------------------
! Local
!------------------------------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real hi_urb(its:ite,1:nz_um,jts:jte) ! Height histograms of buildings
real hi_urb1D(nz_um) ! Height histograms of buildings
real hb_u(nz_um) ! Bulding's heights
real ss_urb(nz_um) ! Probability that a building has an height equal to z
real pb_urb(nz_um) ! Probability that a building has an height greater or equal to z
integer nz_urb(nurbm) ! Number of layer in the urban grid
integer nzurban(nurbm)
! Building parameters
real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
real twini_u(nurbm) ! Initial temperature inside the building's wall [K]
real trini_u(nurbm) ! Initial temperature inside the building's roof [K]
real tgini_u(nurbm) ! Initial road temperature
!
! Building materials
!
real csg(ng_u) ! Specific heat of the ground material [J m^3 K^-1]
real csr(nwr_u) ! Specific heat of the roof material [J m^3 K^-1]
real csw(nwr_u) ! Specific heat of the wall material [J m^3 K^-1]
real alag(ng_u) ! Ground thermal diffusivity [m^2 s^-1]
real alaw(nwr_u) ! Wall thermal diffusivity [m^2 s^-1]
real alar(nwr_u) ! Roof thermal diffusivity [m^2 s^-1]
!
! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
!
! Radiation parameters
real albg_u(nurbm) ! Albedo of the ground
real albw_u(nurbm) ! Albedo of the wall
real albr_u(nurbm) ! Albedo of the roof
real emg_u(nurbm) ! Emissivity of ground
real emw_u(nurbm) ! Emissivity of wall
real emr_u(nurbm) ! Emissivity of roof
! fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
! and the short wave radation.
real fww_u(nz_um,nz_um,ndm,nurbm) ! from wall to wall
real fwg_u(nz_um,ndm,nurbm) ! from wall to ground
real fgw_u(nz_um,ndm,nurbm) ! from ground to wall
real fsw_u(nz_um,ndm,nurbm) ! from sky to wall
real fws_u(nz_um,ndm,nurbm) ! from sky to wall
real fsg_u(ndm,nurbm) ! from sky to ground
! Roughness parameters
real z0g_u(nurbm) ! The ground's roughness length
real z0r_u(nurbm) ! The roof's roughness length
! Roughness parameters
real z0(ndm,nz_um) ! Roughness lengths "profiles"
! Street parameters
integer nd_u(nurbm) ! Number of street direction for each urban class
real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
real drst_u(ndm,nurbm) ! Street direction
real ws_u(ndm,nurbm) ! Street width
real bs_u(ndm,nurbm) ! Building width
real h_b(nz_um,nurbm) ! Bulding's heights
real d_b(nz_um,nurbm) ! Probability that a building has an height h_b
real ss_u(nz_um,nurbm) ! Probability that a building has an height equal to z
real pb_u(nz_um,nurbm) ! Probability that a building has an height greater or equal to z
!
! Street parameters
!
real bs(ndm) ! Building width
real ws(ndm) ! Street width
real drst(ndm) ! street directions
real strd(ndm) ! Street lengths
real ss(nz_um) ! Probability to have a building with height h
real pb(nz_um) ! Probability to have a building with an height equal
! Grid parameters
integer nz_u(nurbm) ! Number of layer in the urban grid
real z_u(nz_um) ! Height of the urban grid levels
! 1D array used for the input and output of the routine "urban"
real z1D(kms:kme) ! vertical coordinates
real ua1D(kms:kme) ! wind speed in the x directions
real va1D(kms:kme) ! wind speed in the y directions
real pt1D(kms:kme) ! potential temperature
real da1D(kms:kme) ! air density
real pr1D(kms:kme) ! air pressure
real pt01D(kms:kme) ! reference potential temperature
real zr1D ! zenith angle
real deltar1D ! declination of the sun
real ah1D ! hour angle (it should come from the radiation routine)
real rs1D ! solar radiation
real rld1D ! downward flux of the longwave radiation
real tw1D(2*ndm,nz_um,nwr_u) ! temperature in each layer of the wall
real tg1D(ndm,ng_u) ! temperature in each layer of the ground
real tr1D(ndm,nz_um,nwr_u) ! temperature in each layer of the roof
real sfw1D(2*ndm,nz_um) ! sensible heat flux from walls
real sfg1D(ndm) ! sensible heat flux from ground (road)
real sfr1D(ndm,nz_um) ! sensible heat flux from roofs
real sf1D(kms:kme) ! surface of the urban grid cells
real vl1D(kms:kme) ! volume of the urban grid cells
real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
real time_bep
! arrays used to collapse indexes
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
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, &
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
!------------------------------------------------------------------------
! Calculation of the momentum, heat and turbulent kinetic fluxes
! produced by builgings
!
! Reference:
! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
! 261-304
!------------------------------------------------------------------------
!prepare the arrays to collapse indexes
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
if (num_urban_hi.ge.nz_um)then
write(*,*)'nz_um too small, please increase to at least ', num_urban_hi+1
stop
endif
do ix=its,ite
do iy=jts,jte
do iz_u=1,nz_um
hi_urb(ix,iz_u,iy)=0.
enddo
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 !iz
do iz_u=1,num_urban_hi
hi_urb(ix,iz_u,iy)= hi_urb2d(ix,iz_u,iy)
enddo !iz_u
enddo
enddo
if (first) then ! True only on first call
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)
! Initialisation of the urban parameters and calculation of the view factors
call icBEP(nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u)
first=.false.
endif ! first
do ix=its,ite
do iy=jts,jte
if (FRC_URB2D(ix,iy).gt.0.) then ! Calling BEP only for existing urban classes.
iurb=UTYPE_URB2D(ix,iy)
hi_urb1D=0.
do iz_u=1,nz_um
hi_urb1D(iz_u)=hi_urb(ix,iz_u,iy)
enddo
call icBEPHI_XY(hb_u,hi_urb1D,ss_urb,pb_urb, &
nz_urb(iurb),z_u)
call param(iurb,nz_u(iurb),nz_urb(iurb),nzurban(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_urb,ss,pb_u, &
pb_urb,pb,lp_urb2d(ix,iy), &
lb_urb2d(ix,iy),hgt_urb2d(ix,iy),FRC_URB2D(ix,iy))
!
!We compute the view factors in the icBEP_XY routine
!
call icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u,fws_u,fsg_u, &
nd_u(iurb),strd,ws,nzurban(iurb),z_u)
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)=th_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
! tw1D(2*id-1,iz_u,iw)=tw1_u(ix,iy,ind_zwd(iz_u,iw,id))
! tw1D(2*id,iz_u,iw)=tw2_u(ix,iy,ind_zwd(iz_u,iw,id))
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)=tg_u(ix,iy,ind_gd(ig,id))
tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
enddo
do iz_u=1,nz_um
do ir=1,nwr_u
! tr1D(id,iz_u,ir)=tr_u(ix,iy,ind_zwd(iz_u,ir,id))
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(ix,iy,ind_zd(iz,id))
! sfw1D(2*id,iz)=sfw2(ix,iy,ind_zd(iz,id))
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(ix,iy,id)
sfg1D(id)=sfg_urb3d(ix,id,iy)
enddo
do id=1,ndm
do iz=1,nz_um
! sfr1D(id,iz)=sfr(ix,iy,ind_zd(iz,id))
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 angle(xlong(ix,iy),xlat(ix,iy),julday,time_h,zr1D,deltar1D,ah1D)
! write(*,*) 'entro en BEP1D'
call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
zr1D,deltar1D,ah1D,rs1D,rld1D, &
alag,alaw,alar,csg,csw,csr, &
albg_u(iurb),albw_u(iurb),albr_u(iurb), &
emg_u(iurb),emw_u(iurb),emr_u(iurb), &
fww_u,fwg_u,fgw_u,fsw_u, &
fws_u,fsg_u,z0, &
nd_u(iurb),strd,drst,ws,bs,ss,pb, &
nzurban(iurb),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,sf1D,vl1D,rl_up(ix,iy), &
rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy))
! write(*,*) 'salgo de BEP1D'
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
sf(ix,kts:kte,iy)=0.
vl(ix,kts:kte,iy)=0.
a_u(ix,kts:kte,iy)=0.
a_v(ix,kts:kte,iy)=0.
a_t(ix,kts:kte,iy)=0.
a_e(ix,kts:kte,iy)=0.
b_u(ix,kts:kte,iy)=0.
b_v(ix,kts:kte,iy)=0.
b_t(ix,kts:kte,iy)=0.
b_e(ix,kts:kte,iy)=0.
b_q(ix,kts:kte,iy)=0.
dlg(ix,kts:kte,iy)=0.
dl_u(ix,kts:kte,iy)=0.
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 ! FRC_URB2D
enddo ! iy
enddo ! ix
time_bep=time_bep+dt
return
end subroutine BEP
! ===6=8===============================================================72
subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, & 2,27
zr,deltar,ah,rs,rld, &
alag,alaw,alar,csg,csw,csr, &
albg,albw,albr,emg,emw,emr, &
fww,fwg,fgw,fsw,fws,fsg,z0, &
ndu,strd,drst,ws,bs,ss,pb, &
nzu,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,sf,vl,rl_up,rs_abs,emiss,grdflx_urb)
! ----------------------------------------------------------------------
! This routine computes the effects of buildings on momentum, heat and
! TKE (turbulent kinetic energy) sources or sinks and on the mixing length.
! It provides momentum, heat and TKE sources or sinks at different levels of a
! mesoscale grid defined by the altitude of its cell interfaces "z" and
! its number of levels "nz".
! The meteorological input parameters (wind, temperature, solar radiation)
! are specified on the "mesoscale grid".
! The inputs concerning the building and street charateristics are defined
! on a "urban grid". The "urban grid" is defined with its number of levels
! "nz_u" and its space step "dz_u".
! The input parameters are interpolated on the "urban grid". The sources or sinks
! are calculated on the "urban grid". Finally the sources or sinks are
! interpolated on the "mesoscale grid".
! Mesoscale grid Urban grid Mesoscale grid
!
! z(4) --- ---
! | |
! | |
! | Interpolation Interpolation |
! | Sources or sinks calculation |
! z(3) --- ---
! | ua ua_u --- uv_a a_u |
! | va va_u | uv_b b_u |
! | pt pt_u --- uh_b a_v |
! z(2) --- | etc... etc...---
! | z_u(1) --- |
! | | |
! z(1) ------------------------------------------------------------
!
! Reference:
! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
! 261-304
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
! integer nz ! Number of vertical levels
integer kms,kme,kts,kte
real z(kms:kme) ! Altitude above the ground of the cell interfaces.
real ua(kms:kme) ! Wind speed in the x direction
real va(kms:kme) ! Wind speed in the y direction
real pt(kms:kme) ! Potential temperature
real da(kms:kme) ! Air density
real pr(kms:kme) ! Air pressure
real pt0(kms:kme) ! Reference potential temperature (could be equal to "pt")
real dt ! Time step
real zr ! Zenith angle
real deltar ! Declination of the sun
real ah ! Hour angle
real rs ! Solar radiation
real rld ! Downward flux of the longwave radiation
! Data relative to the "urban grid"
integer iurb ! Current urban class
! Building parameters
real alag(ng_u) ! Ground thermal diffusivity [m^2 s^-1]
real alaw(nwr_u) ! Wall thermal diffusivity [m^2 s^-1]
real alar(nwr_u) ! Roof thermal diffusivity [m^2 s^-1]
real csg(ng_u) ! Specific heat of the ground material [J m^3 K^-1]
real csw(nwr_u) ! Specific heat of the wall material [J m^3 K^-1]
real csr(nwr_u) ! Specific heat of the roof material [J m^3 K^-1]
! Radiation parameters
real albg ! Albedo of the ground
real albw ! Albedo of the wall
real albr ! Albedo of the roof
real emg ! Emissivity of ground
real emw ! Emissivity of wall
real emr ! Emissivity of roof
! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and
! short wave radation.
! The calculation of these factor is explained in the Appendix A of the BLM paper
real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
real fwg(nz_um,ndm,nurbm) ! from wall to ground
real fgw(nz_um,ndm,nurbm) ! from ground to wall
real fsw(nz_um,ndm,nurbm) ! from sky to wall
real fws(nz_um,ndm,nurbm) ! from wall to sky
real fsg(ndm,nurbm) ! from sky to ground
! Roughness parameters
real z0(ndm,nz_um) ! Roughness lengths "profiles"
! Street parameters
integer ndu ! Number of street direction for each urban class
real strd(ndm) ! Street length (set to a greater value then the horizontal length of the cells)
real drst(ndm) ! Street direction
real ws(ndm) ! Street width
real bs(ndm) ! Building width
real ss(nz_um) ! The probability that a building has an height equal to "z"
real pb(nz_um) ! The probability that a building has an height greater or equal to "z"
! Grid parameters
integer nzu ! Number of layer in the urban grid
real z_u(nz_um) ! Height of the urban grid levels
! ----------------------------------------------------------------------
! INPUT-OUTPUT
! ----------------------------------------------------------------------
! Data relative to the "urban grid" which should be stored from the current time step to the next one
real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
real sfw(2*ndm,nz_um) ! Sensible heat flux from walls
real sfg(ndm) ! Sensible heat flux from ground (road)
real sfr(ndm,nz_um) ! Sensible heat flux from roofs
real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) towards the interior
real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof towards the interior
real gfw(2*ndm,nz_um) ! Heat flux transfered from the surface of the walls towards the interior
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
! Implicit and explicit components of the source and sink terms at each levels,
! the fluxes can be computed as follow: FX = A*X + B example: Heat fluxes = a_t * pt + b_t
real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
real a_t(kms:kme) ! Implicit component of the heat sources or sinks
real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
real b_t(kms:kme) ! Explicit component of the heat sources or sinks
real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
real dz(kms:kme) ! vertical space steps of the "mesoscale grid"
! Data interpolated from the "mesoscale grid" to the "urban grid"
real ua_u(nz_um) ! Wind speed in the x direction
real va_u(nz_um) ! Wind speed in the y direction
real pt_u(nz_um) ! Potential temperature
real da_u(nz_um) ! Air density
real pt0_u(nz_um) ! Reference potential temperature
real pr_u(nz_um) ! Air pressure
! Solar radiation at each level of the "urban grid"
real rsg(ndm) ! Short wave radiation from the ground
real rsw(2*ndm,nz_um) ! Short wave radiation from the walls
real rlg(ndm) ! Long wave radiation from the ground
real rlw(2*ndm,nz_um) ! Long wave radiation from the walls
! Potential temperature of the surfaces at each level of the "urban grid"
real ptg(ndm) ! Ground potential temperatures
real ptr(ndm,nz_um) ! Roof potential temperatures
real ptw(2*ndm,nz_um) ! Walls potential temperatures
! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
! The fluxes can be computed as follow: Fluxes of X = A*X + B
! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
!
real rs_abs ! solar radiation absorbed by urban surfaces
real rl_up ! longwave radiation emitted by urban surface to the atmosphere
real emiss ! mean emissivity of the urban surface
real grdflx_urb ! ground heat flux
integer iz,id
integer iw,ix,iy
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
! Fix some usefull parameters for the computation of the sources or sinks
do iz=kts,kte
dz(iz)=z(iz+1)-z(iz)
end do
! Interpolation on the "urban grid"
call interpol(kms,kme,kts,kte,nzu,z,z_u,ua,ua_u)
call interpol(kms,kme,kts,kte,nzu,z,z_u,va,va_u)
call interpol(kms,kme,kts,kte,nzu,z,z_u,pt,pt_u)
call interpol(kms,kme,kts,kte,nzu,z,z_u,pt0,pt0_u)
call interpol(kms,kme,kts,kte,nzu,z,z_u,pr,pr_u)
call interpol(kms,kme,kts,kte,nzu,z,z_u,da,da_u)
! Compute the modification of the radiation due to the buildings
call modif_rad(iurb,ndu,nzu,z_u,ws, &
drst,strd,ss,pb, &
tw,tg,albg,albw,emw,emg, &
fww,fwg,fgw,fsw,fsg, &
zr,deltar,ah, &
rs,rld,rsw,rsg,rlw,rlg)
! calculation of the urban albedo and the upward long wave radiation
call upward_rad(ndu,nzu,ws,bs, &
sigma,pb,ss, &
tg,emg,albg,rlg,rsg,sfg, &
tw,emw,albw,rlw,rsw,sfw, &
tr,emr,albr,rld,rs,sfr, &
rs_abs,rl_up,emiss,grdflx_urb)
! Compute the surface temperatures
call surf_temp(nzu,ndu,pr_u,dt,ss, &
rs,rld,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)
! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
call buildings(ndu,nzu,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)
! Calculation of the sensible heat fluxes for the ground, the wall and roof
! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
! where A and B are the implicit and explicit components of the heat sources or sinks.
!
!
do id=1,ndu
sfg(id)=-da_u(1)*cp_u*thb_u(id,1)
do iz=2,nzu
sfr(id,iz)=-da_u(iz)*cp_u*thb_u(id,iz)
enddo
do iz=1,nzu
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
! calculation of the urban albedo and the upward long wave radiation
!! call upward_rad(ndu,nzu,ws,bs, &
!! sigma,pb,ss, &
!! tg,emg,albg,rlg,rsg,sfg, &
!! tw,emw,albw,rlw,rsw,sfw, &
!! tr,emr,albr,rld,rs,sfr, &
!! rs_abs,rl_up,emiss,grdflx_urb)
! Interpolation on the "mesoscale grid"
call urban_meso(ndu,kms,kme,kts,kte,nzu,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)
! computation of the mean road temperature tsk (this value could be used
! to replace the surface temperature in the radiation routines, if needed).
! Calculation of the length scale taking into account the buildings effects
call interp_length(ndu,kms,kme,kts,kte,nzu,z_u,z,ss,ws,bs,dlg,dl_u)
return
end subroutine BEP1D
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine param(iurb,nzu,nzurb,nzurban,ndu, & 2
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_urb,ss,pb_u, &
pb_urb,pb,lp_urb,lb_urb,hgt_urb,frc_urb)
! ----------------------------------------------------------------------
! This routine prepare some usefull parameters
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer iurb ! Current urban class
integer nzu ! Number of vertical urban levels in the current class
integer nzurb ! Number of vertical urban levels in the current class
integer ndu ! Number of street direction for the current urban class
real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
real bs_u(ndm,nurbm) ! Building width
real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
real drst_u(ndm,nurbm) ! Street direction
real strd_u(ndm,nurbm) ! Street length
real ws_u(ndm,nurbm) ! Street width
real z0g_u(nurbm) ! The ground's roughness length
real z0r_u(nurbm) ! The roof's roughness length
real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
real ss_urb(nz_um) ! The probability that a building has an height equal to "z"
real pb_urb(nz_um) ! The probability that a building has an height greater or equal to "z"
real lp_urb ! Building plan area density
real lb_urb ! Building surface area to plan area ratio
real hgt_urb ! Average building height weighted by building plan area [m]
real frc_urb ! Urban fraction
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real alag(ng_u) ! Ground thermal diffusivity at each ground levels
real alar(nwr_u) ! Roof thermal diffusivity at each roof levels
real alaw(nwr_u) ! Wall thermal diffusivity at each wall levels
real csg(ng_u) ! Specific heat of the ground material at each ground levels
real csr(nwr_u) ! Specific heat of the roof material at each roof levels
real csw(nwr_u) ! Specific heat of the wall material at each wall levels
real bs(ndm) ! Building width for the current urban class
real drst(ndm) ! street directions for the current urban class
real strd(ndm) ! Street lengths for the current urban class
real ws(ndm) ! Street widths of the current urban class
real z0(ndm,nz_um) ! Roughness lengths "profiles"
real ss(nz_um) ! Probability to have a building with height h
real pb(nz_um) ! Probability to have a building with an height equal
integer nzurban
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer id,ig,ir,iw,iz,ihu
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
!
!Initialize
!
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.
nzurban=0
ihu=0
do iz=1,nz_um
if (ss_urb(iz)/=0.) then
ihu=1
exit
else
continue
endif
enddo
if (ihu==1) then
do iz=1,nzurb+1
ss(iz)=ss_urb(iz)
pb(iz)=pb_urb(iz)
enddo
nzurban=nzurb
else
do iz=1,nzu+1
ss(iz)=ss_u(iz,iurb)
pb(iz)=pb_u(iz,iurb)
end do
nzurban=nzu
endif
do id=1,ndu
z0(id,1)=z0g_u(iurb)
do iz=2,nzurban+1
z0(id,iz)=z0r_u(iurb)
enddo
enddo
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,ndu
strd(id)=strd_u(id,iurb)
drst(id)=drst_u(id,iurb)
enddo
do id=1,ndu
if ((hgt_urb<=0.).OR.(lp_urb<=0.).OR.(lb_urb<=0.)) then
ws(id)=ws_u(id,iurb)
bs(id)=bs_u(id,iurb)
else if ((lp_urb/frc_urb<1.).and.(lp_urb<lb_urb)) then
bs(id)=2.*hgt_urb*lp_urb/(lb_urb-lp_urb)
ws(id)=2.*hgt_urb*lp_urb*((frc_urb/lp_urb)-1.)/(lb_urb-lp_urb)
else
ws(id)=ws_u(id,iurb)
bs(id)=bs_u(id,iurb)
endif
enddo
do id=1,ndu
if ((bs(id)<=1.).OR.(bs(id)>=150.)) then
! write(*,*) 'WARNING, WIDTH OF THE BUILDING WRONG',id,bs(id)
! write(*,*) 'WIDTH OF THE STREET',id,ws(id)
bs(id)=bs_u(id,iurb)
ws(id)=ws_u(id,iurb)
endif
if ((ws(id)<=1.).OR.(ws(id)>=150.)) then
! write(*,*) 'WARNING, WIDTH OF THE STREET WRONG',id,ws(id)
! write(*,*) 'WIDTH OF THE BUILDING',id,bs(id)
bs(id)=bs_u(id,iurb)
ws(id)=ws_u(id,iurb)
endif
enddo
return
end subroutine param
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u) 13
! ----------------------------------------------------------------------
! This routine interpolate para
! meters from the "mesoscale grid" to
! the "urban grid".
! See p300 Appendix B.1 of the BLM paper.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
integer kts,kte,kms,kme
real z(kms:kme) ! Altitude of the cell interface
real c(kms:kme) ! Parameter which has to be interpolated
! Data relative to the "urban grid"
integer nz_u ! Number of levels
!! real z_u(nz_u+1) ! Altitude of the cell interface
real z_u(nz_um) ! Altitude of the cell interface
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
!! real c_u(nz_u) ! Interpolated paramters in the "urban grid"
real c_u(nz_um) ! Interpolated paramters in the "urban grid"
! LOCAL:
! ----------------------------------------------------------------------
integer iz_u,iz
real ctot,dz
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, & 2,6
tw,tg,albg,albw,emw,emg, &
fww,fwg,fgw,fsw,fsg, &
zr,deltar,ah, &
rs,rl,rsw,rsg,rlw,rlg)
! ----------------------------------------------------------------------
! This routine computes the modification of the short wave and
! long wave radiation due to the buildings.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer iurb ! current urban class
integer nd ! Number of street direction for the current urban class
integer nz_u ! Number of layer in the urban grid
real z(nz_um) ! Height of the urban grid levels
real ws(ndm) ! Street widths of the current urban class
real drst(ndm) ! street directions for the current urban class
real strd(ndm) ! Street lengths for the current urban class
real ss(nz_um) ! probability to have a building with height h
real pb(nz_um) ! probability to have a building with an height equal
real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
real albg ! Albedo of the ground for the current urban class
real albw ! Albedo of the wall for the current urban class
real emg ! Emissivity of ground for the current urban class
real emw ! Emissivity of wall for the current urban class
real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
real fsg(ndm,nurbm) ! View factors from sky to ground
real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
real fws(nz_um,ndm,nurbm) ! View factors from wall to sky
real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
real ah ! Hour angle (it should come from the radiation routine)
real zr ! zenith angle
real deltar ! Declination of the sun
real rs ! solar radiation
real rl ! downward flux of the longwave radiation
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real rlg(ndm) ! Long wave radiation at the ground
real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
real rsg(ndm) ! Short wave radiation at the ground
real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer id,iz
! Calculation of the shadow effects
call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
rs,rsw,rsg)
! Calculation of the reflection effects
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine surf_temp(nz_u,nd,pr,dt,ss,rs,rl,rsg,rlg,rsw,rlw, & 2,5
tg,alag,csg,emg,albg,ptg,sfg,gfg, &
tr,alar,csr,emr,albr,ptr,sfr,gfr, &
tw,alaw,csw,emw,albw,ptw,sfw,gfw)
! ----------------------------------------------------------------------
! Computation of the surface temperatures for walls, ground and roofs
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer nz_u ! Number of vertical layers defined in the urban grid
integer nd ! Number of street direction for the current urban class
real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
real alar(nwr_u) ! Roof thermal diffusivity for the current urban class [m^2 s^-1]
real alaw(nwr_u) ! Wall thermal diffusivity for the current urban class [m^2 s^-1]
real albg ! Albedo of the ground for the current urban class
real albr ! Albedo of the roof for the current urban class
real albw ! Albedo of the wall for the current urban class
real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
real dt ! Time step
real emg ! Emissivity of ground for the current urban class
real emr ! Emissivity of roof for the current urban class
real emw ! Emissivity of wall for the current urban class
real pr(nz_um) ! Air pressure
real rs ! Solar radiation
real rl ! Downward flux of the longwave radiation
real rlg(ndm) ! Long wave radiation at the ground
real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
real rsg(ndm) ! Short wave radiation at the ground
real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
real sfg(ndm) ! Sensible heat flux from ground (road)
real sfr(ndm,nz_um) ! Sensible heat flux from roofs
real sfw(2*ndm,nz_um) ! Sensible heat flux from walls
real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) toward the interior
real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof toward the interior
real gfw(2*ndm,nz_um) ! Heat flux transfered from the surface of the walls toward the interior
real ss(nz_um) ! Probability to have a building with height h
real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real ptg(ndm) ! Ground potential temperatures
real ptr(ndm,nz_um) ! Roof potential temperatures
real ptw(2*ndm,nz_um) ! Walls potential temperatures
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer id,ig,ir,iw,iz
real rtg(ndm) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
real rtr(ndm,nz_um) ! Total radiation at roof surface (solar+incoming long+outgoing long)
real rtw(2*ndm,nz_um) ! Radiation at walls surface (solar+incoming long+outgoing long)
real tg_tmp(ng_u)
real tr_tmp(nwr_u)
real tw_tmp(nwr_u)
real dzg_u(ng_u) ! Layer sizes in the ground
real dzr_u(nwr_u) ! Layers sizes in the roof
real dzw_u(nwr_u) ! Layer sizes in the wall
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/
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
do id=1,nd
! Calculation for the ground surfaces
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
! Calculation for the roofs surfaces
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 !iz
! Calculation for the walls surfaces
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 !iz
end do !id
return
end subroutine surf_temp
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine buildings(nd,nz,z0,ua_u,va_u,pt_u,pt0_u, & 2,8
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)
! ----------------------------------------------------------------------
! This routine computes the sources or sinks of the different quantities
! on the urban grid. The actual calculation is done in the subroutines
! called flux_wall, and flux_flat.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer nd ! Number of street direction for the current urban class
integer nz ! number of vertical space steps
real ua_u(nz_um) ! Wind speed in the x direction on the urban grid
real va_u(nz_um) ! Wind speed in the y direction on the urban grid
real da_u(nz_um) ! air density on the urban grid
real drst(ndm) ! Street directions for the current urban class
real dz
real pt_u(nz_um) ! Potential temperature on the urban grid
real pt0_u(nz_um) ! reference potential temperature on the urban grid
real ptg(ndm) ! Ground potential temperatures
real ptr(ndm,nz_um) ! Roof potential temperatures
real ptw(2*ndm,nz_um) ! Walls potential temperatures
real ss(nz_um) ! probability to have a building with height h
real z0(ndm,nz_um) ! Roughness lengths "profiles"
real dt ! time step
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
! vertical surfaces (walls) and horizontal surfaces (roofs and street)
! The fluxes can be computed as follow: Fluxes of X = A*X + B
! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer id,iz
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
dz=dz_u
do id=1,nd
! Calculation at the ground surfaces
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))
! Calculation at the roof surfaces
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
! Calculation at the wall surfaces
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl, & 2
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)
! ----------------------------------------------------------------------
! This routine interpolates the parameters from the "urban grid" to the
! "mesoscale grid".
! See p300-301 Appendix B.2 of the BLM paper.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
integer kms,kme,kts,kte
real z(kms:kme) ! Altitude above the ground of the cell interface
real dz(kms:kme) ! Vertical space steps
! Data relative to the "uban grid"
integer nz_u ! Number of layer in the urban grid
integer nd ! Number of street direction for the current urban class
real bs(ndm) ! Building widths of the current urban class
real ws(ndm) ! Street widths of the current urban class
real z_u(nz_um) ! Height of the urban grid levels
real pb(nz_um) ! Probability to have a building with an height equal
real ss(nz_um) ! Probability to have a building with height h
real uhb_u(ndm,nz_um) ! U (x-wind component) Horizontal surfaces, B (explicit) term
real uva_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, A (implicit) term
real uvb_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, B (explicit) term
real vhb_u(ndm,nz_um) ! V (y-wind component) Horizontal surfaces, B (explicit) term
real vva_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, A (implicit) term
real vvb_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, B (explicit) term
real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Data relative to the "mesoscale grid"
real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
real a_t(kms:kme) ! Implicit component of the heat sources or sinks
real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
real b_t(kms:kme) ! Explicit component of the heat sources or sinks
real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
real dzz
real fact
integer id,iz,iz_u
real se,sr,st,su,sv
real uet(kms:kme) ! Contribution to TKE due to walls
real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
! initialisation
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
! horizontal surfaces
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
! volume
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
! horizontal surface impact
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
! vertical surface impact
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, & 2
dlg,dl_u)
! ----------------------------------------------------------------------
! Calculation of the length scales
! See p272-274 formula (22) and (24) of the BLM paper
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer kms,kme,kts,kte
real z(kms:kme) ! Altitude above the ground of the cell interface
integer nd ! Number of street direction for the current urban class
integer nz_u ! Number of levels in the "urban grid"
real z_u(nz_um) ! Height of the urban grid levels
real bs(ndm) ! Building widths of the current urban class
real ss(nz_um) ! Probability to have a building with height h
real ws(ndm) ! Street widths of the current urban class
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
real dlgtmp
integer id,iz,iz_u
real sftot
real ulu,ssl
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, & 2,4
rs,rsw,rsg)
! ----------------------------------------------------------------------
! Modification of short wave radiation to take into account
! the shadow produced by the buildings
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer nd ! Number of street direction for the current urban class
integer nz_u ! number of vertical layers defined in the urban grid
real ah ! Hour angle (it should come from the radiation routine)
real deltar ! Declination of the sun
real drst(ndm) ! street directions for the current urban class
real rs ! solar radiation
real ss(nz_um) ! probability to have a building with height h
real pb(nz_um) ! Probability that a building has an height greater or equal to h
real ws(ndm) ! Street width of the current urban class
real z(nz_um) ! Height of the urban grid levels
real zr ! zenith angle
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real rsg(ndm) ! Short wave radiation at the ground
real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer id,iz,jz
real aae,aaw,bbb,phix,rd,rtot,wsd
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
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
!test
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd) 4
! ----------------------------------------------------------------------
! This routine computes the effects of a shadow induced by a building of
! height hu, on a portion of wall between z1 and z2. See equation A10,
! and correction described below formula A11, and figure A1. Basically rd
! is the ratio between the horizontal surface illuminated and the portion
! of wall. Referring to figure A1, multiplying radiation flux density on
! a horizontal surface (rs) by x1-x2 we have the radiation energy per
! unit time. Dividing this by z2-z1, we obtain the radiation flux
! density reaching the portion of the wall between z2 and z1
! (everything is assumed in 2D)
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
real aa ! Angle between the sun direction and the face of the wall (A12)
real hu ! Height of the building that generates the shadow
real phix ! Solar zenith angle
real ws ! Width of the street
real z1 ! Height of the level z(iz)
real z2 ! Height of the level z(iz+1)
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real rd ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A.
! Multiplying rd by rs (radiation flux
! density on a horizontal surface) gives
! the radiation flux density on the
! portion of wall between z1 and z2.
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
real x1,x2 ! x1,x2 see Fig. A1.
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine long_rad(iurb,nz_u,id,emw,emg, & 2,2
fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
! ----------------------------------------------------------------------
! This routine computes the effects of the reflections of long-wave
! radiation in the street canyon by solving the system
! of 2*nz_u+1 eqn. in 2*nz_u+1
! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
! The system is solved by solving A X= B,
! with A matrix B vector, and X solution.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
real emg ! Emissivity of ground for the current urban class
real emw ! Emissivity of wall for the current urban class
real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
real fsg(ndm,nurbm) ! View factors from sky to ground
real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
integer id ! Current street direction
integer iurb ! Current urban class
integer nz_u ! Number of layer in the urban grid
real pb(nz_um) ! Probability to have a building with an height equal
real rl ! Downward flux of the longwave radiation
real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real rlg(ndm) ! Long wave radiation at the ground
real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer i,j
real aaa(2*nz_um+1,2*nz_um+1) ! terms of the matrix
real bbb(2*nz_um+1) ! terms of the vector
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
! west wall
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)*pb(i+1)
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
! east wall
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)*pb(i-nz_u+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
! ground
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine short_rad(iurb,nz_u,id,albw, & 2,2
albg,fwg,fww,fgw,rsg,rsw,pb)
! ----------------------------------------------------------------------
! This routine computes the effects of the reflections of short-wave
! (solar) radiation in the street canyon by solving the system
! of 2*nz_u+1 eqn. in 2*nz_u+1
! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
! The system is solved by solving A X= B,
! with A matrix B vector, and X solution.
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
real albg ! Albedo of the ground for the current urban class
real albw ! Albedo of the wall for the current urban class
real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
integer id ! current street direction
integer iurb ! current urban class
integer nz_u ! Number of layer in the urban grid
real pb(nz_um) ! probability to have a building with an height equal
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real rsg(ndm) ! Short wave radiation at the ground
real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer i,j
real aaa(2*nz_um+1,2*nz_um+1) ! terms of the matrix
real bbb(2*nz_um+1) ! terms of the vector
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
! west wall
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
! east wall
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
! ground
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine gaussj(a,n,b,np) 4,4
! ----------------------------------------------------------------------
! This routine solve a linear system of n equations of the form
! A X = B
! where A is a matrix a(i,j)
! B a vector and X the solution
! In output b is replaced by the solution
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer np
real a(np,np)
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real b(np)
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer nmax
parameter (nmax=150)
real big,dum
integer i,icol,irow
integer j,k,l,ll,n
integer ipiv(nmax)
real pivinv
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
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
CALL wrf_error_fatal('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) CALL wrf_error_fatal('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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine soil_temp(nz,dz,temp,pt,ala,cs, & 5,2
rs,rl,press,dt,em,alb,rt,sf,gf)
! ----------------------------------------------------------------------
! This routine solves the Fourier diffusion equation for heat in
! the material (wall, roof, or ground). Resolution is done implicitely.
! Boundary conditions are:
! - fixed temperature at the interior
! - energy budget at the surface
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer nz ! Number of layers
real ala(nz) ! Thermal diffusivity in each layers [m^2 s^-1]
real alb ! Albedo of the surface
real cs(nz) ! Specific heat of the material [J m^3 K^-1]
real dt ! Time step
real em ! Emissivity of the surface
real press ! Pressure at ground level
real rl ! Downward flux of the longwave radiation
real rs ! Solar radiation
real sf ! Sensible heat flux at the surface
real temp(nz) ! Temperature in each layer [K]
real dz(nz) ! Layer sizes [m]
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real gf ! Heat flux transferred from the surface toward the interior
real pt ! Potential temperature at the surface
real rt ! Total radiation at the surface (solar+incoming long+outgoing long)
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer iz
real a(nz,3)
real alpha
real c(nz)
real cddz(nz+2)
real tsig
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
tsig=temp(nz)
alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
! Compute cddz=2*cd/dz
cddz(1)=ala(1)/dz(1)
do iz=2,nz
cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
enddo
! cddz(nz+1)=ala(nz+1)/dz(nz)
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=-cddz(nz)*(temp(nz)-temp(nz-1))*cs(nz)
gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
return
end subroutine soil_temp
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine invert(n,a,c,x) 4
! ----------------------------------------------------------------------
! Inversion and resolution of a tridiagonal matrix
! A X = C
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
integer n
real a(n,3) ! a(*,1) lower diagonal (Ai,i-1)
! a(*,2) principal diagonal (Ai,i)
! a(*,3) upper diagonal (Ai,i+1)
real c(n)
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
real x(n)
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
integer i
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine flux_wall(ua,va,pt,da,ptw,uva,vva,uvb,vvb, & 4
tva,tvb,evb,drst,dt)
! ----------------------------------------------------------------------
! This routine computes the surface sources or sinks of momentum, tke,
! and heat from vertical surfaces (walls).
! ----------------------------------------------------------------------
implicit none
! INPUT:
! -----
real drst ! street directions for the current urban class
real da ! air density
real pt ! potential temperature
real ptw ! Walls potential temperatures
real ua ! wind speed
real va ! wind speed
real dt !time step
! OUTPUT:
! ------
! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
! vertical surfaces (walls).
! The fluxes can be computed as follow: Fluxes of X = A*X + B
! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
real uva ! U (wind component) Vertical surfaces, A (implicit) term
real uvb ! U (wind component) Vertical surfaces, B (explicit) term
real vva ! V (wind component) Vertical surfaces, A (implicit) term
real vvb ! V (wind component) Vertical surfaces, B (explicit) term
real tva ! Temperature Vertical surfaces, A (implicit) term
real tvb ! Temperature Vertical surfaces, B (explicit) term
real evb ! Energy (TKE) Vertical surfaces, B (explicit) term
! LOCAL:
! -----
real hc
real u_ort
real vett
! -------------------------
! END VARIABLES DEFINITIONS
! -------------------------
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
! tva=-hc/da/cp_u
!!!!!!!!!!!!!!!!!!!!
! explicit
tvb=hc*ptw/da/cp_u-hc/da/cp_u*pt !c
tva = 0. !c
evb=cdrag*(abs(u_ort)**3.)/2.
return
end subroutine flux_wall
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg, & 4
uhb,vhb,thb,ehb)
! ----------------------------------------------------------------------
! Calculation of the flux at the ground
! Formulation of Louis (Louis, 1979)
! ----------------------------------------------------------------------
implicit none
! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------
real dz ! first vertical level
real pt ! potential temperature
real pt0 ! reference potential temperature
real ptg ! ground potential temperature
real ua ! wind speed
real va ! wind speed
real z0 ! Roughness length
! ----------------------------------------------------------------------
! OUTPUT:
! ----------------------------------------------------------------------
! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
! surfaces (roofs and street)
! The fluxes can be computed as follow: Fluxes of X = B
! Example: Momentum fluxes on horizontal surfaces = uhb_u
real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
real thb ! Temperature Horizontal surfaces, B (explicit) term
real tva ! Temperature Vertical surfaces, A (implicit) term
real tvb ! Temperature Vertical surfaces, B (explicit) term
real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
! ----------------------------------------------------------------------
! LOCAL:
! ----------------------------------------------------------------------
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)
! ----------------------------------------------------------------------
! END VARIABLES DEFINITIONS
! ----------------------------------------------------------------------
! computation of the ground temperature
utot=(ua**2+va**2)**.5
!!!! Louis formulation
!
! compute the bulk Richardson Number
zz=dz/2.
! if(tstar.lt.0.)then
! wstar=(-ustar*tstar*g*hii/pt)**(1./3.)
! else
! wstar=0.
! endif
!
! if (utot.le.0.7*wstar) utot=max(0.7*wstar,0.00001)
utot=max(utot,0.01)
ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
aa=vk/log(zz/z0)
! determine the parameters fm and fh for stable, neutral and unstable conditions
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
! thb= 0.
ehb=buu
!!!!!!!!!!!!!!!
return
end subroutine flux_flat
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine icBEP (nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u) 2
implicit none
! Street parameters
integer nd_u(nurbm) ! Number of street direction for each urban class
real h_b(nz_um,nurbm) ! Bulding's heights [m]
real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
! -----------------------------------------------------------------------
! Output
!------------------------------------------------------------------------
real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to z
real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to z
! Grid parameters
integer nz_u(nurbm) ! Number of layer in the urban grid
real z_u(nz_um) ! Height of the urban grid levels
! -----------------------------------------------------------------------
! Local
!------------------------------------------------------------------------
integer iz_u,id,ilu,iurb
real dtot
real hbmax
!------------------------------------------------------------------------
! -----------------------------------------------------------------------
! This routine initialise the urban paramters for the BEP module
!------------------------------------------------------------------------
!
!Initialize variables
!
z_u=0.
nz_u=0
ss_u=0.
pb_u=0.
! Computation of the urban levels height
z_u(1)=0.
do iz_u=1,nz_um-1
z_u(iz_u+1)=z_u(iz_u)+dz_u
enddo
! Normalisation of the building density
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
! Compute the view factors, pb and ss
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)
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
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws) 2,26
implicit none
! -----------------------------------------------------------------------
! Input
!------------------------------------------------------------------------
integer iurb ! Number of the urban class
integer nz_u ! Number of levels in the urban grid
integer id ! Street direction number
real ws ! Street width
real z(nz_um) ! Height of the urban grid levels
real dxy ! Street lenght
! -----------------------------------------------------------------------
! Output
!------------------------------------------------------------------------
! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
! and the short wave radation. They are the part of radiation from a surface
! or from the sky to another surface.
real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
real fwg(nz_um,ndm,nurbm) ! from wall to ground
real fgw(nz_um,ndm,nurbm) ! from ground to wall
real fsw(nz_um,ndm,nurbm) ! from sky to wall
real fws(nz_um,ndm,nurbm) ! from wall to sky
real fsg(ndm,nurbm) ! from sky to ground
! -----------------------------------------------------------------------
! Local
!------------------------------------------------------------------------
integer jz,iz
real hut
real f1,f2,f12,f23,f123,ftot
real fprl,fnrm
real a1,a2,a3,a4,a12,a23,a123
! -----------------------------------------------------------------------
! This routine calculates the view factors
!------------------------------------------------------------------------
hut=z(nz_u+1)
do jz=1,nz_u
! radiation from wall to wall
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
! radiation from ground to wall
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
! radiation from sky to wall
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
! radiation from wall to sky
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
! radiation from wall to ground
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
! radiation from sky to ground
call fprls (fprl,dxy,ws,hut)
fsg(id,iurb)=fprl
return
end subroutine view_factors
! ===6=8===============================================================72
! ===6=8===============================================================72
SUBROUTINE fprls (fprl,a,b,c) 10
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
! ===6=8===============================================================72
! ===6=8===============================================================72
SUBROUTINE fnrms (fnrm,a,b,c) 16
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
! ===6=8===============================================================72
SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,& 2,1
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)
! initialization routine, where the variables from the table are read
implicit none
integer iurb ! urban class number
! Building parameters
real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
real twini_u(nurbm) ! Temperature inside the buildings behind the wall [K]
real trini_u(nurbm) ! Temperature inside the buildings behind the roof [K]
real tgini_u(nurbm) ! Initial road temperature
! Radiation parameters
real albg_u(nurbm) ! Albedo of the ground
real albw_u(nurbm) ! Albedo of the wall
real albr_u(nurbm) ! Albedo of the roof
real emg_u(nurbm) ! Emissivity of ground
real emw_u(nurbm) ! Emissivity of wall
real emr_u(nurbm) ! Emissivity of roof
! Roughness parameters
real z0g_u(nurbm) ! The ground's roughness length
real z0r_u(nurbm) ! The roof's roughness length
! Street parameters
integer nd_u(nurbm) ! Number of street direction for each urban class
real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
real drst_u(ndm,nurbm) ! Street direction [degree]
real ws_u(ndm,nurbm) ! Street width [m]
real bs_u(ndm,nurbm) ! Building width [m]
real h_b(nz_um,nurbm) ! Bulding's heights [m]
real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
integer i,iu
integer nurb ! number of urban classes used
!
!Initialize some variables
!
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)
! ----------------
!
! Computation of the solar angles
! schayes (1982,atm. env. , p1407)
! Inputs
!========================
! along=longitud
! alat=latitude
! day=julian day (from the beginning of the year)
! realt= time GMT in hours
! Outputs
!============================
! zr=solar zenith angle
! deltar=declination angle
! ah=hour angle
!===============================
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)
!
! declination
!
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)
!
! time equation in hours
!
et=0.0072*c1-0.0528*c2-0.0012*c3-0.1229*s1-0.1565*s2-0.0041*s3
!
!
! hour angle
!
! ifh=0
! ahor=realt-12.+ifh+et+alongt
ahor=realt-12.+et+alongt
ah=ahor*radh
chor=cos(ah)
!
! zenith angle
!
coznt=sphi*sid+cphi*cd*chor
zr=acos(coznt)
return
END SUBROUTINE angle
!
!====6=8===============================================================72
!====6=8===============================================================72
subroutine upward_rad(ndu,nzu,ws,bs,sigma,pb,ss, & 2
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)
!
! IN this surboutine we compute the upward longwave flux, and the albedo
! needed for the radiation scheme
!
implicit none
!
!INPUT VARIABLES
!
real rsw(2*ndm,nz_um) ! Short wave radiation at the wall for a given canyon direction [W/m2]
real rlw(2*ndm,nz_um) ! Long wave radiation at the walls for a given canyon direction [W/m2]
real rsg(ndm) ! Short wave radiation at the canyon for a given canyon direction [W/m2]
real rlg(ndm) ! Long wave radiation at the ground for a given canyon direction [W/m2]
real rs ! Short wave radiation at the horizontal surface from the sun [W/m²]
real sfw(2*ndm,nz_um) ! Sensible heat flux from walls [W/m²]
real sfg(ndm) ! Sensible heat flux from ground (road) [W/m²]
real sfr(ndm,nz_um) ! Sensible heat flux from roofs [W/m²]
real rld ! Long wave radiation from the sky [W/m²]
real albg_u ! albedo of the ground/street
real albw_u ! albedo of the walls
real albr_u ! albedo of the roof
real ws(ndm) ! width of the street
real bs(ndm)
! building size
real pb(nz_um) ! Probability to have a building with an height equal or higher
integer nzu
real ss(nz_um) ! Probability to have a building of a given height
real sigma
real emg_u ! emissivity of the street
real emw_u ! emissivity of the wall
real emr_u ! emissivity of the roof
real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
integer id ! street direction
integer ndu ! number of street directions
!OUTPUT/INPUT
real rs_abs ! absrobed solar radiationfor this street direction
real rl_up ! upward longwave radiation for this street direction
real emiss ! mean emissivity
real grdflx_urb ! ground heat flux
!LOCAL
integer iz,iw
real rl_inc,rl_emit
real gfl
integer ix,iy,iwrong
iwrong=1
do iz=1,nzu+1
do id=1,ndu
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
if(tw(2*id-1,iz,iw).lt.100.) then
write(*,*)'in upward_rad ',iz,id,iw,tw(2*id-1,iz,iw)
iwrong=0
endif
if(tw(2*id,iz,iw).lt.100.) then
write(*,*)'in upward_rad ',iz,id,iw,tw(2*id,iz,iw)
iwrong=0
endif
enddo
enddo
enddo
do id=1,ndu
do iw=1,ng_u
if(tg(id,iw).lt.100.) then
write(*,*)'in upward_rad ',id,iw,tg(id,iw)
iwrong=0
endif
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,ndu
rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/ndu
rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/ndu
rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/ndu
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))/ndu
do iz=2,nzu
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))/ndu
rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
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))/ndu
enddo
do iz=1,nzu
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))/ndu
rl_inc=rl_inc+(( rlw(2*id-1,iz)+rlw(2*id,iz) ) )*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
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))/ndu
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))/ndu
enddo
enddo
emiss=(emg_u+emw_u+emr_u)/3.
rl_up=(rl_inc+rl_emit)-rld
return
END SUBROUTINE upward_rad
!====6=8===============================================================72
!====6=8===============================================================72
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u, & 2,2
fws_u,fsg_u,ndu,strd,ws,nzu,z_u)
implicit none
! Street parameters
integer ndu ! Number of street direction for each urban class
integer iurb
real strd(ndm) ! Street length (fix to greater value to the horizontal length of the cells)
real ws(ndm) ! Street width [m]
! Grid parameters
integer nzu ! Number of layer in the urban grid
real z_u(nz_um) ! Height of the urban grid levels
! -----------------------------------------------------------------------
! Output
!------------------------------------------------------------------------
! fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
! and the short wave radation. They are the part of radiation from a surface
! or from the sky to another surface.
real fww_u(nz_um,nz_um,ndm,nurbm) ! from wall to wall
real fwg_u(nz_um,ndm,nurbm) ! from wall to ground
real fgw_u(nz_um,ndm,nurbm) ! from ground to wall
real fsw_u(nz_um,ndm,nurbm) ! from sky to wall
real fws_u(nz_um,ndm,nurbm) ! from sky to wall
real fsg_u(ndm,nurbm) ! from sky to ground
! -----------------------------------------------------------------------
! Local
!------------------------------------------------------------------------
integer id
! -----------------------------------------------------------------------
! This routine compute the view factors
!------------------------------------------------------------------------
!
!Initialize
!
fww_u=0.
fwg_u=0.
fgw_u=0.
fsw_u=0.
fws_u=0.
fsg_u=0.
do id=1,ndu
call view_factors(iurb,nzu,id,strd(id),z_u,ws(id), &
fww_u,fwg_u,fgw_u,fsg_u,fsw_u,fws_u)
enddo
return
end subroutine icBEP_XY
! ===6=8===============================================================72
! ===6=8===============================================================72
subroutine icBEPHI_XY(hb_u,hi_urb1D,ss_u,pb_u,nzu,z_u) 2
implicit none
!-----------------------------------------------------------------------
! Inputs
!-----------------------------------------------------------------------
! Street parameters
!
real hi_urb1D(nz_um) ! The probability that a building has an height h_b
!
! Grid parameters
!
real z_u(nz_um) ! Height of the urban grid levels
! -----------------------------------------------------------------------
! Output
!------------------------------------------------------------------------
real ss_u(nz_um) ! The probability that a building has an height equal to z
real pb_u(nz_um) ! The probability that a building has an height greater or equal to z
!
! Grid parameters
!
integer nzu ! Number of layer in the urban grid
! -----------------------------------------------------------------------
! Local
!------------------------------------------------------------------------
real hb_u(nz_um) ! Bulding's heights [m]
integer iz_u,id,ilu
real dtot
real hbmax
!------------------------------------------------------------------------
!Initialize variables
!
nzu=0
ss_u=0.
pb_u=0.
! Normalisation of the building density
dtot=0.
hb_u=0.
do ilu=1,nz_um
dtot=dtot+hi_urb1D(ilu)
enddo
do ilu=1,nz_um
if (hi_urb1D(ilu)<0.) then
! write(*,*) 'WARNING, HI_URB1D(ilu) < 0 IN BEP'
go to 20
endif
enddo
if (dtot.gt.0.) then
continue
else
! write(*,*) 'WARNING, HI_URB1D <= 0 IN BEP'
go to 20
endif
do ilu=1,nz_um
hi_urb1D(ilu)=hi_urb1D(ilu)/dtot
enddo
hb_u(1)=dz_u
do ilu=2,nz_um
hb_u(ilu)=dz_u+hb_u(ilu-1)
enddo
! Compute pb and ss
hbmax=0.
do ilu=1,nz_um
if (hi_urb1D(ilu)>0.and.hi_urb1D(ilu)<=1.) then
hbmax=hb_u(ilu)
endif
enddo
do iz_u=1,nz_um-1
if(z_u(iz_u+1).gt.hbmax)go to 10
enddo
10 continue
nzu=iz_u+1
if ((nzu+1).gt.nz_um) then
write(*,*) 'error, nz_um has to be increased to at least',nzu+1
stop
endif
do iz_u=1,nzu
ss_u(iz_u)=0.
do ilu=1,nz_um
if(z_u(iz_u).le.hb_u(ilu) &
.and.z_u(iz_u+1).gt.hb_u(ilu))then
ss_u(iz_u)=ss_u(iz_u)+hi_urb1D(ilu)
endif
enddo
enddo
pb_u(1)=1.
do iz_u=1,nzu
pb_u(iz_u+1)=max(0.,pb_u(iz_u)-ss_u(iz_u))
enddo
20 continue
return
end subroutine icBEPHI_XY
! ===6=8===============================================================72
! ===6=8===============================================================72
END MODULE module_sf_bep