<p><b>laura@ucar.edu</b> 2011-01-13 16:30:40 -0700 (Thu, 13 Jan 2011)</p><p>One more for Noah Land Surface Model<br>
</p><hr noshade><pre><font color="gray">Added: branches/atmos_physics/src/core_physics/physics_wrf/module_sf_bep_bem.F
===================================================================
--- branches/atmos_physics/src/core_physics/physics_wrf/module_sf_bep_bem.F (rev 0)
+++ branches/atmos_physics/src/core_physics/physics_wrf/module_sf_bep_bem.F 2011-01-13 23:30:40 UTC (rev 686)
@@ -0,0 +1,4101 @@
+MODULE module_sf_bep_bem
+
+!USE module_model_constants
+ USE module_sf_urban
+ USE module_sf_bem
+
+! 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=13)
+
+ 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)
+
+ integer nf_u !Number of grid levels in the floors (BEM)
+ parameter (nf_u=10)
+
+ integer ngb_u !Number of grid levels in the ground below building (BEM)
+ parameter (ngb_u=10)
+
+ real dz_u ! Urban grid resolution
+ parameter (dz_u=5.)
+
+ integer nbui_max !maximum number of types of buildings in an urban class
+ parameter (nbui_max=4) !must be less or equal than nz_um
+
+!---------------------------------------------------------------------------------
+!Parameters of the windows. The glasses of windows are considered without films -
+!Read the paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour -
+!of the total solar energy transmittance of windows".Solar Energy Vol.69,No.4, -
+!pp. 321-329, for more details. -
+!---------------------------------------------------------------------------------
+ integer p_num !number of panes in the windows (1,2 or 3)
+ parameter (p_num=2)
+ integer q_num !category number for the windows (q_num= 4, standard glasses)
+ parameter(q_num=4) !Possible values 1,2,...,10
+
+
+! 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
+ real latent ! Latent heat of vaporization [J/kg] (used in BEM)
+ 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,latent=2.45e+06)
+
+! -----------------------------------------------------------------------
+
+
+
+
+ CONTAINS
+
+ subroutine BEP_BEM(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, &
+ th_phy,rho,p_phy,swdown,glw, &
+ gmt,julday,xlong,xlat, &
+ declin_urb,cosz_urb2d,omg_urb2d, &
+ num_urban_layers, &
+ trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
+ tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
+ tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
+ cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d, &
+ sfwin1_urb3d,sfwin2_urb3d, &
+ sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
+ 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,qv_phy, &
+ 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
+ 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
+!New variables used for BEM
+ REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: qv_phy
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tlev_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: qlev_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tglev_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tflev_urb3d
+ REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
+ REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
+ REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
+ REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
+ REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
+!End variables
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
+ REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
+
+ real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates
+ REAL, INTENT(IN ):: DT ! Time step
+
+!
+!------------------------------------------------------------------------
+! Output
+!------------------------------------------------------------------------
+!
+! 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(ims:ime,jms:jme) ! upward long wave radiation
+ real rs_abs(ims:ime,jms:jme) ! absorbed short wave radiation
+ real emiss(ims:ime,jms:jme) ! emissivity averaged for urban surfaces
+ real grdflx_urb(ims:ime,jms:jme) ! ground heat flux for urban areas
+!------------------------------------------------------------------------
+! Local
+!------------------------------------------------------------------------
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! 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
+!
+! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
+!
+! Radiation paramters
+ 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 albwin_u(nurbm) ! Albedo of the windows
+ real emwind_u(nurbm) ! Emissivity of windows
+ real emg_u(nurbm) ! Emissivity of ground
+ real emw_u(nurbm) ! Emissivity of wall
+ real emr_u(nurbm) ! Emissivity of roof
+
+! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
+! and the short wave radation.
+ 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 sky to wall
+ real fsg(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
+
+! 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
+
+
+! Grid parameters
+ integer nz_u(nurbm) ! Number of layer in the urban grid
+
+ real z_u(nz_um) ! Height of the urban grid levels
+
+! MT
+ real cop_u(nurbm)
+ real pwin_u(nurbm)
+ real beta_u(nurbm)
+ integer sw_cond_u(nurbm)
+ real time_on_u(nurbm)
+ real time_off_u(nurbm)
+ real targtemp_u(nurbm)
+ real gaptemp_u(nurbm)
+ real targhum_u(nurbm)
+ real gaphum_u(nurbm)
+ real perflo_u(nurbm)
+ real hsesf_u(nurbm)
+ real hsequip(24)
+
+! 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,nbui_max) ! 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
+!
+!New variable for BEM
+!
+ real tlev1D(nz_um,nbui_max) ! temperature in each floor and in each different type of building
+ real qlev1D(nz_um,nbui_max) ! specific humidity in each floor and in each different type of building
+ real twlev1D(2*ndm,nz_um,nbui_max) ! temperature in each window in each floor in each different type of building
+ real tglev1D(ndm,ngb_u,nbui_max) ! temperature in each layer of the ground below of a type of building
+ real tflev1D(ndm,nf_u,nz_um-1,nbui_max)! temperature in each layer of the floors in each building
+ real lflev1D(nz_um,nz_um) ! latent heat flux due to the air conditioning systems
+ real sflev1D(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems
+ real lfvlev1D(nz_um,nz_um) ! latent heat flux due to ventilation
+ real sfvlev1D(nz_um,nz_um) ! sensible heat flux due to ventilation
+ real sfwin1D(2*ndm,nz_um,nbui_max) ! sensible heat flux from windows
+ real consumlev1D(nz_um,nz_um) ! consumption due to the air conditioning systems
+ real qv1D(kms:kme) ! specific humidity
+ real meso_urb ! constant to link meso and urban scales [m¯2]
+ real d_urb(nz_um)
+ real sf_ac
+ integer ibui,nbui
+ integer nlev(nz_um)
+!
+!End new variables
+!
+
+ real sfw1D(2*ndm,nz_um,nbui_max) ! 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_ac1D(kms:kme)
+ real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
+ real b_q1D(kms:kme) ! Explicit component of the Humidity 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(nbui_max,nz_um,nwr_u,ndm)
+ integer ind_gd(ng_u,ndm)
+ integer ind_zd(nbui_max,nz_um,ndm)
+ integer ind_zdf(nz_um,ndm)
+ integer ind_zrd(nz_um,nwr_u,ndm)
+!
+ integer ind_bd(nbui_max,nz_um)
+ integer ind_wd(nbui_max,nz_um,ndm)
+ integer ind_gbd(nbui_max,ngb_u,ndm)
+ integer ind_fbd(nbui_max,nf_u,nz_um-1,ndm)
+!
+ integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
+ integer it, nint
+ integer iii
+ real 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,fww,fwg,fgw,fsw,fws,fsg, &
+ z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
+ nz_u,z_u,albwin_u,emwind_u , &
+ cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, targtemp_u, &
+ gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip
+
+!------------------------------------------------------------------------
+! 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
+!
+! F. Salamanca and A. Martilli, 2009: 'A new Building Energy Model coupled
+! with an Urban Canopy Parameterization for urban climate simulations_part II.
+! Validation with one dimension off-line simulations'. Theor Appl Climatol
+! DOI 10.1007/s00704-009-0143-8
+!------------------------------------------------------------------------
+!prepare the arrays to collapse indexes
+ if(num_urban_layers.lt.nbui_max*nz_um*ndm*max(nwr_u,ng_u))then
+ write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um*ndm*max(nwr_u,ng_u)
+ stop
+ endif
+!
+!New conditions for BEM
+!
+ if(num_urban_layers.lt.nbui_max*nz_um)then !limit for indoor temperature and indoor humidity
+ write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um
+ stop
+ endif
+
+ if(num_urban_layers.lt.nbui_max*nz_um*ndm)then !limit for window temperature
+ write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um*ndm
+ stop
+ endif
+
+ if(num_urban_layers.lt.nbui_max*ndm*ngb_u)then !limit for ground temperature below a building
+ write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*ndm*ngb_u
+ stop
+ endif
+
+ if(num_urban_layers.lt.(nz_um-1)*nbui_max*ndm*nf_u)then !limit for floor temperature
+ write(*,*)'num_urban_layers too small, please increase to at least ', nbui_max*ndm*nf_u*(nz_um-1),num_urban_layers
+ stop
+ endif
+
+ if (ndm.ne.2)then
+ write(*,*) 'number of directions is not correct',ndm
+ stop
+ endif
+
+!
+!End of new conditions
+!
+!
+!Initialize collapse indexes
+!
+ ind_zwd=0
+ ind_gd=0
+ ind_zd=0
+ ind_zdf=0
+ ind_zrd=0
+ ind_bd=0
+ ind_wd=0
+ ind_gbd=0
+ ind_fbd=0
+!
+!End initialization indexes
+!
+
+ iii=0
+ do ibui=1,nbui_max
+ do iz_u=1,nz_um
+ do iw=1,nwr_u
+ do id=1,ndm
+ iii=iii+1
+ ind_zwd(ibui,iz_u,iw,id)=iii
+ enddo
+ 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 ibui=1,nbui_max
+ do iz_u=1,nz_um
+ do id=1,ndm
+ iii=iii+1
+ ind_zd(ibui,iz_u,id)=iii
+ enddo
+ enddo
+ enddo
+
+ iii=0
+ do iz_u=1,nz_um
+ do iw=1,nwr_u
+ do id=1,ndm
+ iii=iii+1
+ ind_zrd(iz_u,iw,id)=iii
+ enddo
+ enddo
+ enddo
+
+!
+!New indexes for BEM
+!
+ iii=0
+ do iz_u=1,nz_um
+ do id=1,ndm
+ iii=iii+1
+ ind_zdf(iz_u,id)=iii
+ enddo ! id
+ enddo ! iz_u
+
+ iii=0
+ do ibui=1,nbui_max !Type of building
+ do iz_u=1,nz_um !vertical levels
+ iii=iii+1
+ ind_bd(ibui,iz_u)=iii
+ enddo !iz_u
+ enddo !ibui
+
+
+ iii=0
+ do ibui=1,nbui_max !type of building
+ do iz_u=1,nz_um !vertical levels
+ do id=1,ndm !direction
+ iii=iii+1
+ ind_wd(ibui,iz_u,id)=iii
+ enddo !id
+ enddo !iz_u
+ enddo !ibui
+
+ iii=0
+ do ibui=1,nbui_max!type of building
+ do iw=1,ngb_u !layers in the wall (ground below a building)
+ do id=1,ndm !direction
+ iii=iii+1
+ ind_gbd(ibui,iw,id)=iii
+ enddo !id
+ enddo !iw
+ enddo !ibui
+
+ iii=0
+ do ibui=1,nbui_max !type of building
+ do iw=1,nf_u !layers in the wall (floor)
+ do iz_u=1,nz_um-1 !vertical levels
+ do id=1,ndm !direction
+ iii=iii+1
+ ind_fbd(ibui,iw,iz_u,id)=iii
+ enddo !id
+ enddo !iz_u
+ enddo !iw
+ enddo !ibui
+
+!
+!End of new indexes
+!
+
+ do ix=its,ite
+ do iy=jts,jte
+ z(ix,kts,iy)=0.
+ do iz=kts+1,kte+1
+ z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
+ enddo
+ enddo
+ enddo
+
+ if (first) then ! 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,albwin_u,emg_u,emw_u,&
+ emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,&
+ cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
+ targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)
+
+!Initialisation of the urban parameters and calculation of the view factor
+ call icBEP(fww,fwg,fgw,fsw,fws,fsg, &
+ z0g_u,z0r_u, &
+ nd_u,strd_u,ws_u,bs_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)
+ ibui=0
+ nlev=0
+ nbui=0
+ d_urb=0.
+ do iz=1,nz_um
+ if(ss_u(iz,iurb).gt.0) then
+ ibui=ibui+1
+ nlev(ibui)=iz-1
+ d_urb(ibui)=ss_u(iz,iurb)
+ nbui=ibui
+ endif
+ end do !iz
+ if (nbui.gt.nbui_max) then
+ write (*,*) 'nbui_max must be increased to',nbui
+ stop
+ endif
+
+ 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)
+ qv1D(iz)=qv_phy(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_ac1D(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
+ do ibui=1,nbui_max
+ tw1D(2*id-1,iz_u,iw,ibui)=tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
+ tw1D(2*id,iz_u,iw,ibui)=tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do id=1,ndm
+ do ig=1,ng_u
+ tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
+ enddo
+ do iz_u=1,nz_um
+ do ir=1,nwr_u
+ tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)
+ enddo
+ enddo
+ enddo
+!
+!Initialize variables for BEM
+!
+ tlev1D=0. !Indoor temperature
+ qlev1D=0. !Indoor humidity
+
+ twlev1D=0. !Window temperature
+ tglev1D=0. !Ground temperature
+ tflev1D=0. !Floor temperature
+
+ sflev1D=0. !Sensible heat flux from the a.c.
+ lflev1D=0. !latent heat flux from the a.c.
+ consumlev1D=0.!consumption of the a.c.
+ sfvlev1D=0. !Sensible heat flux from natural ventilation
+ lfvlev1D=0. !Latent heat flux from natural ventilation
+ sfwin1D=0. !Sensible heat flux from windows
+ sfw1D=0. !Sensible heat flux from walls
+
+ do iz_u=1,nz_um !vertical levels
+ do ibui=1,nbui_max !Type of building
+ tlev1D(iz_u,ibui)= tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)
+ qlev1D(iz_u,ibui)= qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)
+ enddo !ibui
+ enddo !iz_u
+
+ do id=1,ndm !direction
+ do iz_u=1,nz_um !vertical levels
+ do ibui=1,nbui_max !type of building
+ twlev1D(2*id-1,iz_u,ibui)=tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
+ twlev1D(2*id,iz_u,ibui)=tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
+ sfwin1D(2*id-1,iz_u,ibui)=sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
+ sfwin1D(2*id,iz_u,ibui)=sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
+ enddo !ibui
+ enddo !iz_u
+ enddo !id
+
+ do id=1,ndm !direction
+ do iw=1,ngb_u !layer in the wall
+ do ibui=1,nbui_max !type of building
+ tglev1D(id,iw,ibui)=tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)
+ enddo !ibui
+ enddo !iw
+ enddo !id
+
+ do id=1,ndm !direction
+ do iw=1,nf_u !layer in the walls
+ do iz_u=1,nz_um-1 !verticals levels
+ do ibui=1,nbui_max !type of building
+ tflev1D(id,iw,iz_u,ibui)=tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)
+ enddo !ibui
+ enddo ! iz_u
+ enddo !iw
+ enddo !id
+
+!
+!End initialization for BEM
+!
+
+ do id=1,ndm
+ do iz=1,nz_um
+ do ibui=1,nbui_max !type of building
+ sfw1D(2*id-1,iz,ibui)=sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)
+ sfw1D(2*id,iz,ibui)=sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)
+ enddo
+ enddo
+ enddo
+
+ do id=1,ndm
+ sfg1D(id)=sfg_urb3d(ix,id,iy)
+ enddo
+
+ do id=1,ndm
+ do iz=1,nz_um
+ sfr1D(id,iz)=sfr_urb3d(ix,ind_zdf(iz,id),iy)
+ enddo
+ enddo
+
+ rs1D=swdown(ix,iy)
+ rld1D=glw(ix,iy)
+
+ zr1D=acos(COSZ_URB2D(ix,iy))
+ deltar1D=DECLIN_URB
+ ah1D=OMG_URB2D(ix,iy)
+
+ call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
+ zr1D,deltar1D,ah1D,rs1D,rld1D, &
+ alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
+ albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, &
+ emwind_u,fww,fwg,fgw,fsw,fws,fsg, &
+ z0g_u,z0r_u, &
+ nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
+ nz_u,z_u, &
+ cop_u,pwin_u,beta_u,sw_cond_u,time_on_u, &
+ time_off_u,targtemp_u,gaptemp_u,targhum_u, &
+ gaphum_u, perflo_u, hsesf_u, hsequip, &
+ tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D, &
+ a_u1D,a_v1D,a_t1D,a_e1D, &
+ b_u1D,b_v1D,b_t1D,b_ac1D,b_e1D,b_q1D, &
+ dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy), &
+ rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy), &
+ qv1D,tlev1D,qlev1D,sflev1D,lflev1D,consumlev1D, &
+ sfvlev1D,lfvlev1D,twlev1D,tglev1D,tflev1D,sfwin1D,&
+ ix,iy)
+
+ do id=1,ndm ! direction
+ do iz=1,nz_um !vertical levels
+ do ibui=1,nbui_max !type of building
+ sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id-1,iz,ibui)
+ sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id,iz,ibui)
+ enddo
+ 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_zdf(iz,id),iy)=sfr1D(id,iz)
+ enddo
+ enddo
+
+ do id=1,ndm
+ do iz_u=1,nz_um
+ do iw=1,nwr_u
+ do ibui=1,nbui_max
+ tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw,ibui)
+ tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw,ibui)
+ enddo
+ 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_zrd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
+ enddo
+ enddo
+ enddo
+!
+!Outputs of BEM
+!
+ do ibui=1,nbui_max !type of building
+ do iz_u=1,nz_um !vertical levels
+ tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=tlev1D(iz_u,ibui)
+ qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=qlev1D(iz_u,ibui)
+ enddo !iz_u
+ enddo !ibui
+ do id=1,ndm !direction
+ do iz_u=1,nz_um !vertical levels
+ do ibui=1,nbui_max !type of building
+ tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id-1,iz_u,ibui)
+ tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id,iz_u,ibui)
+ sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id-1,iz_u,ibui)
+ sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id,iz_u,ibui)
+ enddo !ibui
+ enddo !iz_u
+ enddo !id
+
+ do id=1,ndm !direction
+ do iw=1,ngb_u !layers in the walls
+ do ibui=1,nbui_max !type of building
+ tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)=tglev1D(id,iw,ibui)
+ enddo !ibui
+ enddo !iw
+ enddo !id
+
+ do id=1,ndm !direction
+ do iw=1,nf_u !layer in the walls
+ do iz_u=1,nz_um-1 !verticals levels
+ do ibui=1,nbui_max !type of building
+ tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)=tflev1D(id,iw,iz_u,ibui)
+ enddo !ibui
+ enddo !iz_u
+ enddo !iw
+ enddo !id
+ sf_ac_urb3d(ix,iy)=0.
+ lf_ac_urb3d(ix,iy)=0.
+ cm_ac_urb3d(ix,iy)=0.
+ sfvent_urb3d(ix,iy)=0.
+ lfvent_urb3d(ix,iy)=0.
+
+ meso_urb=(1./4.)*FRC_URB2D(ix,iy)/((bs_u(1,iurb)+ws_u(1,iurb))*bs_u(2,iurb))+ &
+ (1./4.)*FRC_URB2D(ix,iy)/((bs_u(2,iurb)+ws_u(2,iurb))*bs_u(1,iurb))
+
+
+ ibui=0
+ nlev=0
+ nbui=0
+ d_urb=0.
+ do iz=1,nz_um
+ if(ss_u(iz,iurb).gt.0) then
+ ibui=ibui+1
+ nlev(ibui)=iz-1
+ d_urb(ibui)=ss_u(iz,iurb)
+ nbui=ibui
+ endif
+ end do !iz
+
+ do ibui=1,nbui !type of building
+ do iz_u=1,nlev(ibui) !vertical levels
+ sf_ac_urb3d(ix,iy)=sf_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*sflev1D(iz_u,ibui)
+ lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*lflev1D(iz_u,ibui)
+ cm_ac_urb3d(ix,iy)=cm_ac_urb3d(ix,iy)+meso_urb*d_urb(ibui)*consumlev1D(iz_u,ibui)
+ sfvent_urb3d(ix,iy)=sfvent_urb3d(ix,iy)+meso_urb*d_urb(ibui)*sfvlev1D(iz_u,ibui)
+ lfvent_urb3d(ix,iy)=lfvent_urb3d(ix,iy)+meso_urb*d_urb(ibui)*lfvlev1D(iz_u,ibui)
+ enddo !iz_u
+ enddo !ibui
+!
+!Add the latent heat exchanged throughout the ventilation in the lf_ac_urb3d output variable.
+!it is only a print variable
+!
+! lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+lfvent_urb3d(ix,iy)
+!
+
+ lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)-lfvent_urb3d(ix,iy)
+
+!
+!End outputs of bem
+!
+
+ sf_ac=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)
+ sf_ac=sf_ac+b_ac1D(iz)*da1D(iz)*cp_u*dz8w(ix,iz,iy)*vl1D(iz)*FRC_URB2D(ix,iy)
+ b_e(ix,iz,iy)=b_e1D(iz)
+ b_q(ix,iz,iy)=b_q1D(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_BEM
+
+! ===6=8===============================================================72
+
+ subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, &
+ zr,deltar,ah,rs,rld, &
+ alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
+ albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, &
+ emwind_u,fww,fwg,fgw,fsw,fws,fsg, &
+ z0g_u,z0r_u, &
+ nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
+ nz_u,z_u, &
+ cop_u,pwin_u,beta_u,sw_cond_u,time_on_u, &
+ time_off_u,targtemp_u,gaptemp_u,targhum_u, &
+ gaphum_u, perflo_u, hsesf_u, hsequip, &
+ tw,tg,tr,sfw,sfg,sfr, &
+ a_u,a_v,a_t,a_e, &
+ b_u,b_v,b_t,b_ac,b_e,b_q, &
+ dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb, &
+ qv,tlev,qlev,sflev,lflev,consumlev, &
+ sfvlev,lfvlev,twlev,tglev,tflev,sfwin,ix,iy)
+
+! ----------------------------------------------------------------------
+! 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 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 qv(kms:kme) ! Specific humidity
+ 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_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]
+
+! 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 albwin_u(nurbm) ! Albedo of the windows
+ real emwind_u(nurbm) ! Emissivity of windows
+ real emg_u(nurbm) ! Emissivity of ground
+ real emw_u(nurbm) ! Emissivity of wall
+ real emr_u(nurbm) ! 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 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 (set to a greater value then 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) ! The probability that a building has an height "h_b"
+ 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
+
+! MT
+ real cop_u(nurbm)
+ real pwin_u(nurbm)
+ real beta_u(nurbm)
+ integer sw_cond_u(nurbm)
+ real time_on_u(nurbm)
+ real time_off_u(nurbm)
+ real targtemp_u(nurbm)
+ real gaptemp_u(nurbm)
+ real targhum_u(nurbm)
+ real gaphum_u(nurbm)
+ real perflo_u(nurbm)
+ real hsesf_u(nurbm)
+ real hsequip(24)
+
+! ----------------------------------------------------------------------
+! 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,nbui_max) ! 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,nbui_max) ! 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,nbui_max) ! 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_ac(kms:kme)
+ real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
+ real b_q(kms:kme) ! Explicit component of the humidity 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
+ real qv_u(nz_um) !Specific humidity
+
+! Data defining the building and street charateristics
+
+ real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
+
+ 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 z0(ndm,nz_um) ! Roughness lengths "profiles"
+ real ws(ndm) ! Street widths of the current urban class
+ real bs(ndm) ! Building widths of the current urban class
+ real strd(ndm) ! Street lengths for the current urban class
+ real drst(ndm) ! Street directions 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
+
+! 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,nbui_max) ! 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 tvb_ac(2*ndm,nz_um)
+ 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 qhb_u(ndm,nz_um) ! Humidity Horizontal surfaces, B (explicit) term
+ real qvb_u(2*ndm,nz_um) ! Humidity 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
+ real dt_int ! internal time step
+ integer nt_int ! number of internal time step
+ integer iz,id, it_int
+ integer iw,ix,iy
+
+!---------------------------------------
+!New variables uses in BEM
+!----------------------------------------
+
+ real tmp_u(nz_um) !Air Temperature [K]
+
+ real dzw(nwr_u) !Layer sizes in the walls
+ real dzr(nwr_u) !Layer sizes in the roofs
+ real dzf(nf_u) !Layer sizes in the floors
+ real dzgb(ngb_u) !Layer sizes in the ground below the buildings
+
+ real csgb(ngb_u) !Specific heat of the ground material below the buildings
+ !of the current urban class at each ground levels[J m^3 K^-1]
+ real csf(nf_u) !Specific heat of the floors materials in the buildings
+ !of the current urban class at each levels[J m^3 K^-1]
+ real alar(nwr_u+1) ! Roof thermal diffusivity for the current urban class [W/m K]
+ real alaw(nwr_u+1) ! Walls thermal diffusivity for the current urban class [W/m K]
+ real alaf(nf_u+1) ! Floor thermal diffusivity at each wall layers [W/m K]
+ real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall layer [W/m K]
+
+ real sfrb(ndm,nbui_max) ! Sensible heat flux from roofs [W/m²]
+ real gfrb(ndm,nbui_max) ! Heat flux flowing inside the roofs [W/m²]
+ real sfwb1D(2*ndm,nz_um) !Sensible heat flux from the walls [W/m²]
+ real sfwin(2*ndm,nz_um,nbui_max)!Sensible heat flux from windows [W/m²]
+ real sfwinb1D(2*ndm,nz_um) !Sensible heat flux from windows [W/m²]
+ real gfwb1D(2*ndm,nz_um) !Heat flux flowing inside the walls [W/m²]
+
+ real qlev(nz_um,nbui_max) !specific humidity [kg/kg]
+ real qlevb1D(nz_um) !specific humidity [kg/kg]
+ real tlev(nz_um,nbui_max) !Indoor temperature [K]
+ real tlevb1D(nz_um) !Indoor temperature [K]
+ real twb1D(2*ndm,nwr_u,nz_um) !Wall temperature in BEM [K]
+ real twlev(2*ndm,nz_um,nbui_max) !Window temperature in BEM [K]
+ real twlevb1D(2*ndm,nz_um) !Window temperature in BEM [K]
+ real tglev(ndm,ngb_u,nbui_max) !Ground temperature below a building in BEM [K]
+ real tglevb1D(ngb_u) !Ground temperature below a building in BEM [K]
+ real tflev(ndm,nf_u,nz_um-1,nbui_max)!Floor temperature in BEM[K]
+ real tflevb1D(nf_u,nz_um-1) !Floor temperature in BEM[K]
+ real trb(ndm,nwr_u,nbui_max) !Roof temperature in BEM [K]
+ real trb1D(nwr_u) !Roof temperature in BEM [K]
+
+ real sflev(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems [W]
+ real lflev(nz_um,nz_um) ! latent heat flux due to the air conditioning systems [W]
+ real consumlev(nz_um,nz_um) ! consumption due to the air conditioning systems [W]
+ real sflev1D(nz_um) ! sensible heat flux due to the air conditioning systems [W]
+ real lflev1D(nz_um) ! latent heat flux due to the air conditioning systems [W]
+ real consumlev1D(nz_um) ! consumption due to the air conditioning systems [W]
+ real sfvlev(nz_um,nz_um) ! sensible heat flux due to ventilation [W]
+ real lfvlev(nz_um,nz_um) ! latent heat flux due to ventilation [W]
+ real sfvlev1D(nz_um) ! sensible heat flux due to ventilation [W]
+ real lfvlev1D(nz_um) ! Latent heat flux due to ventilation [W]
+
+ real ptwin(2*ndm,nz_um,nbui_max) ! window potential temperature
+ real tw_av(2*ndm,nz_um) ! Averaged temperature of the wall surfaces
+ real twlev_av(2*ndm,nz_um) ! Averaged temperature of the windows
+ real sfw_av(2*ndm,nz_um) ! Averaged sensible heat from walls
+ real sfwind_av(2*ndm,nz_um) ! Averaged sensible heat from windows
+
+ integer nbui !Total number of different type of buildings in an urban class
+ integer nlev(nz_um) !Number of levels in each different type of buildings in an urban class
+ integer ibui,ily
+ real :: nhourday ! Number of hours from midnight, local time
+! ----------------------------------------------------------------------
+! END VARIABLES DEFINITIONS
+! ----------------------------------------------------------------------
+
+! Fix some usefull parameters for the computation of the sources or sinks
+!
+!initialize inside param
+!
+! ss=0.
+! pb=0.
+
+ do iz=kts,kte
+ dz(iz)=z(iz+1)-z(iz)
+ end do
+ call param(iurb,nz_u(iurb),nd_u(iurb), &
+ csg_u,csg,alag_u,alag,csr_u,csr, &
+ alar_u,alar,csw_u,csw,alaw_u,alaw, &
+ ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
+ strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb, &
+ dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)
+
+! Interpolation on the "urban grid"
+
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,ua,ua_u)
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,va,va_u)
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt,pt_u)
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt0,pt0_u)
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pr,pr_u)
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,da,da_u)
+ call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,qv,qv_u)
+
+! Compute the modification of the radiation due to the buildings
+
+ call averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av, &
+ sfw_av,sfwind_av,sfw,sfwin)
+ call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws, &
+ drst,strd,ss,pb, &
+ tw_av,tg,twlev_av,albg_u(iurb),albw_u(iurb), &
+ emw_u(iurb),emg_u(iurb),pwin_u(iurb),albwin_u(iurb), &
+ emwind_u(iurb),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(nd_u(iurb),nz_u(iurb),ws,bs,sigma,pb,ss, &
+ tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg, &
+ tw_av,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw_av, &
+ tr,emr_u(iurb),albr_u(iurb),emwind_u(iurb), &
+ albwin_u(iurb),twlev_av,pwin_u(iurb),sfwind_av,rld,rs,sfr, &
+ rs_abs,rl_up,emiss,grdflx_urb)
+
+! Compute the surface temperatures
+
+ call surf_temp(nd_u(iurb),pr_u,dt, &
+ rld,rsg,rlg, &
+ tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg)
+
+
+! Call the BEM (Building Energy Model) routine
+
+ do iz=1,nz_um !Compute the outdoor temperature
+ tmp_u(iz)=pt_u(iz)*(pr_u(iz)/p0)**(rcp_u)
+ end do
+
+ ibui=0
+ nlev=0
+ nbui=0
+
+ sfrb=0. !Sensible heat flux from roof
+ gfrb=0. !Heat flux flowing inside the roof
+ sfwb1D=0. !Sensible heat flux from walls
+ sfwinb1D=0. !Sensible heat flux from windows
+ gfwb1D=0. !Heat flux flowing inside the walls[W/m²]
+
+
+ twb1D=0. !Wall temperature
+ twlevb1D=0. !Window temperature
+ tglevb1D=0. !Ground temperature below a building
+ tflevb1D=0. !Floor temperature
+ trb=0. !Roof temperature
+ trb1D=0. !Roof temperature
+
+ qlevb1D=0. !Indoor humidity
+ tlevb1D=0. !indoor temperature
+
+ sflev1D=0. !Sensible heat flux from the a.c.
+ lflev1D=0. !Latent heat flux from the a.c.
+ consumlev1D=0.!Consumption from the a.c.
+ sfvlev1D=0. !Sensible heat flux from the natural ventilation
+ lfvlev1D=0. !Latent heat flux from natural ventilation
+
+ ptw=0. !Wall potential temperature
+ ptwin=0. !Window potential temperature
+ ptr=0. !Roof potential temperature
+
+ do iz=1,nz_um
+ if(ss(iz).gt.0) then
+ ibui=ibui+1
+ nlev(ibui)=iz-1
+ nbui=ibui
+ do id=1,ndm
+ sfrb(id,ibui)=sfr(id,iz)
+ do ily=1,nwr_u
+ trb(id,ily,ibui)=tr(id,iz,ily)
+ enddo
+ enddo
+ endif
+ end do !iz
+!--------------------------------------------------------------------------------
+!Loop over BEM -----------------------------------------------------------------
+!--------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------
+
+ nhourday=ah/PI*180./15.+12.
+ if (nhourday >= 24) nhourday = nhourday - 24
+ if (nhourday < 0) nhourday = nhourday + 24
+
+ do ibui=1,nbui
+
+ do iz=1,nz_um
+ qlevb1D(iz)=qlev(iz,ibui)
+ tlevb1D(iz)=tlev(iz,ibui)
+ enddo
+
+ do id=1,ndm
+
+ do ily=1,nwr_u
+ trb1D(ily)=trb(id,ily,ibui)
+ enddo
+ do ily=1,ngb_u
+ tglevb1D(ily)=tglev(id,ily,ibui)
+ enddo
+
+ do ily=1,nf_u
+ do iz=1,nz_um-1
+ tflevb1D(ily,iz)=tflev(id,ily,iz,ibui)
+ enddo
+ enddo
+
+ do iz=1,nz_um
+ sfwinb1D(2*id-1,iz)=sfwin(2*id-1,iz,ibui)
+ sfwinb1D(2*id,iz)=sfwin(2*id,iz,ibui)
+ enddo
+
+ do iz=1,nz_um
+ do ily=1,nwr_u
+ twb1D(2*id-1,ily,iz)=tw(2*id-1,iz,ily,ibui)
+ twb1D(2*id,ily,iz)=tw(2*id,iz,ily,ibui)
+ enddo
+ sfwb1D(2*id-1,iz)=sfw(2*id-1,iz,ibui)
+ sfwb1D(2*id,iz)=sfw(2*id,iz,ibui)
+ twlevb1D(2*id-1,iz)=twlev(2*id-1,iz,ibui)
+ twlevb1D(2*id,iz)=twlev(2*id,iz,ibui)
+ enddo
+ enddo
+
+ call BEM(nz_um,nlev(ibui),nhourday,dt,bs_u(1,iurb), &
+ bs_u(2,iurb),dz_u,nwr_u,nf_u,nwr_u,ngb_u,sfwb1D,gfwb1D, &
+ sfwinb1D,sfrb(1,ibui),gfrb(1,ibui), &
+ latent,sigma,albw_u(iurb),albwin_u(iurb),albr_u(iurb), &
+ emr_u(iurb),emw_u(iurb),emwind_u(iurb),rsw,rlw,r,cp_u, &
+ da_u,tmp_u,qv_u,pr_u,rs,rld,dzw,csw,alaw,pwin_u(iurb), &
+ cop_u(iurb),beta_u(iurb),sw_cond_u(iurb),time_on_u(iurb), &
+ time_off_u(iurb),targtemp_u(iurb),gaptemp_u(iurb), &
+ targhum_u(iurb),gaphum_u(iurb),perflo_u(iurb),hsesf_u(iurb), &
+ hsequip, &
+ dzf,csf,alaf,dzgb,csgb,alagb,dzr,csr, &
+ alar,tlevb1D,qlevb1D,twb1D,twlevb1D,tflevb1D,tglevb1D, &
+ trb1D,sflev1D,lflev1D,consumlev1D,sfvlev1D,lfvlev1D)
+
+
+!
+!Temporal modifications
+!
+ sfrb(2,ibui)=sfrb(1,ibui)
+ gfrb(2,ibui)=gfrb(1,ibui)
+!
+!End temporal modifications
+!
+ do iz=1,nz_um
+ qlev(iz,ibui)=qlevb1D(iz)
+ tlev(iz,ibui)=tlevb1D(iz)
+ sflev(iz,ibui)=sflev1D(iz)
+ lflev(iz,ibui)=lflev1D(iz)
+ consumlev(iz,ibui)=consumlev1D(iz)
+ sfvlev(iz,ibui)=sfvlev1D(iz)
+ lfvlev(iz,ibui)=lfvlev1D(iz)
+ enddo
+
+ do id=1,ndm
+ do ily=1,nwr_u
+ trb(id,ily,ibui)=trb1D(ily)
+ enddo
+ do ily=1,ngb_u
+ tglev(id,ily,ibui)=tglevb1D(ily)
+ enddo
+
+ do ily=1,nf_u
+ do iz=1,nz_um-1
+ tflev(id,ily,iz,ibui)=tflevb1D(ily,iz)
+ enddo
+ enddo
+
+
+ do iz=1,nz_um
+ do ily=1,nwr_u
+ tw(2*id-1,iz,ily,ibui)=twb1D(2*id-1,ily,iz)
+ tw(2*id,iz,ily,ibui)=twb1D(2*id,ily,iz)
+ enddo
+ gfw(2*id-1,iz,ibui)=gfwb1D(2*id-1,iz)
+ gfw(2*id,iz,ibui)=gfwb1D(2*id,iz)
+ twlev(2*id-1,iz,ibui)=twlevb1D(2*id-1,iz)
+ twlev(2*id,iz,ibui)=twlevb1D(2*id,iz)
+ enddo
+ enddo
+
+ enddo !ibui
+
+!-----------------------------------------------------------------------------
+!End loop over BEM -----------------------------------------------------------
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+
+ ibui=0
+ do iz=1,nz_um
+ if(ss(iz).gt.0) then
+ ibui=ibui+1
+ do id=1,ndm
+ gfr(id,iz)=gfrb(id,ibui)
+ sfr(id,iz)=sfrb(id,ibui)
+ do ily=1,nwr_u
+ tr(id,iz,ily)=trb(id,ily,ibui)
+ enddo
+ ptr(id,iz)=tr(id,iz,nwr_u)*(pr_u(iz)/p0)**(-rcp_u)
+ enddo
+ endif
+ enddo !iz
+
+!Compute the potential temperature for the vertical surfaces of the buildings
+
+ do id=1,ndm
+ do iz=1,nz_um
+ do ibui=1,nbui
+ ptw(2*id-1,iz,ibui)=tw(2*id-1,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u)
+ ptw(2*id,iz,ibui)=tw(2*id,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u)
+ ptwin(2*id-1,iz,ibui)=twlev(2*id-1,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u)
+ ptwin(2*id,iz,ibui)=twlev(2*id,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u)
+ enddo
+ enddo
+ enddo
+
+! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
+
+ call buildings(iurb,nd_u(iurb),nz_u(iurb),z0,ua_u,va_u, &
+ pt_u,pt0_u,ptg,ptr,da_u,ptw,ptwin,pwin_u(iurb),drst, &
+ uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,qvb_u,qhb_u, &
+ uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr, &
+ sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)
+
+! 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.
+
+
+! Interpolation on the "mesoscale grid"
+
+ call urban_meso(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z,dz,z_u,pb,ss,bs,ws,sf, &
+ vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
+ uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u, &
+ a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)
+
+
+! Calculation of the length scale taking into account the buildings effects
+
+ call interp_length(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z_u,z,ss,ws,bs,dlg,dl_u)
+
+ return
+ end subroutine BEP1D
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine param(iurb,nz,nd, &
+ csg_u,csg,alag_u,alag,csr_u,csr, &
+ alar_u,alar,csw_u,csw,alaw_u,alaw, &
+ ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
+ strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb, &
+ dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)
+
+! ----------------------------------------------------------------------
+! This routine prepare some usefull parameters
+! ----------------------------------------------------------------------
+
+ implicit none
+
+
+! ----------------------------------------------------------------------
+! INPUT:
+! ----------------------------------------------------------------------
+ integer iurb ! Current urban class
+ integer nz ! Number of vertical urban levels in the current class
+ integer nd ! 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"
+
+
+! ----------------------------------------------------------------------
+! OUTPUT:
+! ----------------------------------------------------------------------
+ real alag(ng_u) ! Ground thermal diffusivity at each ground levels
+ real csg(ng_u) ! Specific heat of the ground material at each ground 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 greater or equal to "z"
+
+!-----------------------------------------------------------------------------
+!INPUT/OUTPUT
+!-----------------------------------------------------------------------------
+
+ real dzw(nwr_u) !Layer sizes in the walls [m]
+ real dzr(nwr_u) !Layer sizes in the roofs [m]
+ real dzf(nf_u) !Layer sizes in the floors [m]
+ real dzgb(ngb_u) !layer sizes in the ground below the buildings [m]
+
+ 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 csf(nf_u) !Specific heat of the floors materials in the buildings
+ !of the current urban class [J m^3 K^-1]
+ real csgb(ngb_u) !Specific heat of the ground material below the buildings
+ !of the current urban class [J m^3 K^-1]
+ real alar(nwr_u+1) ! Roof thermal diffusivity at each roof levels [W/ m K]
+ real alaw(nwr_u+1) ! Wall thermal diffusivity at each wall levels [W/ m K]
+ real alaf(nf_u+1) ! Floor thermal diffusivity at each wall levels [W/m K]
+ real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall levels [W/m K]
+! ----------------------------------------------------------------------
+! LOCAL:
+! ----------------------------------------------------------------------
+ integer id,ig,ir,iw,iz,iflo
+! ----------------------------------------------------------------------
+! END VARIABLES DEFINITIONS
+! ----------------------------------------------------------------------
+!Define the layer sizes in the walls
+
+ 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.
+ csgb=0.
+ alagb=0.
+ csf=0.
+ alaf=0.
+
+ dzgb=(/0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/)
+ dzr=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
+ dzw=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
+ dzf=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02/)
+
+ do ig=1,ngb_u
+ csgb(ig) = csg_u(iurb)
+ alagb(ig)= csg_u(iurb)*alag_u(iurb)
+ enddo
+ alagb(ngb_u+1)= csg_u(iurb)*alag_u(iurb)
+
+ do iflo=1,nf_u
+ csf(iflo) = csw_u(iurb)
+ alaf(iflo)= csw_u(iurb)*alaw_u(iurb)
+ enddo
+ alaf(nf_u+1)= csw_u(iurb)*alaw_u(iurb)
+
+ do ir=1,nwr_u
+ csr(ir) = csr_u(iurb)
+ alar(ir)= csr_u(iurb)*alar_u(iurb)
+ enddo
+ alar(nwr_u+1)= csr_u(iurb)*alar_u(iurb)
+
+ do iw=1,nwr_u
+ csw(iw) = csw_u(iurb)
+ alaw(iw)= csw_u(iurb)*alaw_u(iurb)
+ enddo
+ alaw(nwr_u+1)=csw_u(iurb)*alaw_u(iurb)
+
+!------------------------------------------------------------------------
+
+ do iz=1,nz+1
+ ss(iz)=ss_u(iz,iurb)
+ pb(iz)=pb_u(iz,iurb)
+ end do
+
+ do ig=1,ng_u
+ csg(ig)=csg_u(iurb)
+ alag(ig)=alag_u(iurb)
+ enddo
+
+ do id=1,nd
+ z0(id,1)=z0g_u(iurb)
+ do iz=2,nz+1
+ z0(id,iz)=z0r_u(iurb)
+ enddo
+ enddo
+
+ do id=1,nd
+ ws(id)=ws_u(id,iurb)
+ bs(id)=bs_u(id,iurb)
+ strd(id)=strd_u(id,iurb)
+ drst(id)=drst_u(id,iurb)
+ enddo
+
+
+ return
+ end subroutine param
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
+
+! ----------------------------------------------------------------------
+! 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_um) ! Altitude of the cell interface
+! ----------------------------------------------------------------------
+! OUTPUT:
+! ----------------------------------------------------------------------
+ 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 averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av,&
+ sfw_av,sfwind_av,sfw,sfwin)
+
+ implicit none
+!
+!INPUT VARIABLES
+!
+ real tw(2*ndm,nz_um,nwr_u,nbui_max) ! Temperature in each layer of the wall [K]
+ real twlev(2*ndm,nz_um,nbui_max) ! Window temperature in BEM [K]
+ real pb(nz_um) ! Probability to have a building with an height equal or greater h
+ real ss(nz_um) ! Probability to have a building with height h
+ real sfw(2*ndm,nz_um,nbui_max) ! Surface fluxes from the walls
+ real sfwin(2*ndm,nz_um,nbui_max) ! Surface fluxes from the windows
+!
+!OUTPUT VARIABLES
+!
+ real tw_av(2*ndm,nz_um) ! Averaged temperature of the wall surfaces
+ real twlev_av(2*ndm,nz_um) ! Averaged temperature of the windows
+ real sfw_av(2*ndm,nz_um) ! Averaged sensible heat from walls
+ real sfwind_av(2*ndm,nz_um) ! Averaged sensible heat from windows
+!
+!LOCAL VARIABLES
+!
+ real d_urb(nz_um)
+ integer nlev(nz_um)
+ integer id,iz
+ integer nbui,ibui
+!
+!Initialize Variables
+!
+ tw_av=0.
+ twlev_av=0.
+ sfw_av=0.
+ sfwind_av=0.
+ ibui=0
+ nbui=0
+ nlev=0
+ d_urb=0.
+
+ do iz=1,nz_um
+ if(ss(iz).gt.0) then
+ ibui=ibui+1
+ d_urb(ibui)=ss(iz)
+ nlev(ibui)=iz-1
+ nbui=ibui
+ endif
+ enddo
+
+ do id=1,ndm
+ do iz=1,nz_um-1
+ if (pb(iz+1).gt.0) then
+ do ibui=1,nbui
+ if (iz.le.nlev(ibui)) then
+ tw_av(2*id-1,iz)=tw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
+ tw(2*id-1,iz,nwr_u,ibui)**4
+ tw_av(2*id,iz)=tw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
+ tw(2*id,iz,nwr_u,ibui)**4
+ twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
+ twlev(2*id-1,iz,ibui)**4
+ twlev_av(2*id,iz)=twlev_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
+ twlev(2*id,iz,ibui)**4
+ sfw_av(2*id-1,iz)=sfw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id-1,iz,ibui)
+ sfw_av(2*id,iz)=sfw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id,iz,ibui)
+ sfwind_av(2*id-1,iz)=sfwind_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id-1,iz,ibui)
+ sfwind_av(2*id,iz)=sfwind_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id,iz,ibui)
+ endif
+ enddo
+ tw_av(2*id-1,iz)=tw_av(2*id-1,iz)**(1./4.)
+ tw_av(2*id,iz)=tw_av(2*id,iz)**(1./4.)
+ twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)**(1./4.)
+ twlev_av(2*id,iz)=twlev_av(2*id,iz)**(1./4.)
+ endif
+ enddo !iz
+ enddo !id
+ return
+ end subroutine averaging_temp
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, &
+ tw,tg,twlev,albg,albw,emw,emg,pwin,albwin, &
+ emwin,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) ! 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
+!
+!New variables BEM
+!
+ real twlev(2*ndm,nz_um) ! Window temperature in BEM [K]
+ real pwin ! Coverage area fraction of windows in the walls of the buildings
+ real albwin ! Albedo of the windows for the current urban class
+ real emwin ! Emissivity of the windows for the current urban class
+ real alb_av ! Averaged albedo (window and wall)
+! ----------------------------------------------------------------------
+! 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,emwin,pwin,twlev, &
+ fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
+
+ alb_av=pwin*albwin+(1.-pwin)*albw
+
+ call short_rad(iurb,nz_u,id,alb_av,albg,fwg,fww,fgw,rsg,rsw,pb)
+
+ enddo
+ return
+ end subroutine modif_rad
+
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine surf_temp(nd,pr,dt,rl,rsg,rlg, &
+ tg,alag,csg,emg,albg,ptg,sfg,gfg)
+
+! ----------------------------------------------------------------------
+! Computation of the surface temperatures for walls, ground and roofs
+! ----------------------------------------------------------------------
+
+ implicit none
+
+! ----------------------------------------------------------------------
+! INPUT:
+! ----------------------------------------------------------------------
+
+ 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 albg ! Albedo of the ground 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 dt ! Time step
+ real emg ! Emissivity of ground for the current urban class
+
+ real pr(nz_um) ! Air pressure
+ real rl ! Downward flux of the longwave radiation
+ real rlg(ndm) ! Long wave radiation at the ground
+ real rsg(ndm) ! Short wave radiation at the ground
+ real sfg(ndm) ! Sensible heat flux from ground (road)
+
+ real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) toward the interior
+
+ real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
+
+
+
+! ----------------------------------------------------------------------
+! OUTPUT:
+! ----------------------------------------------------------------------
+ real ptg(ndm) ! Ground potential temperatures
+
+
+! ----------------------------------------------------------------------
+! LOCAL:
+! ----------------------------------------------------------------------
+ integer id,ig,ir,iw,iz
+
+ real rtg(ndm) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
+
+ real tg_tmp(ng_u)
+
+
+ real dzg_u(ng_u) ! Layer sizes in the ground
+
+ data dzg_u /0.2,0.12,0.08,0.05,0.03,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
+
+ end do !id
+
+ return
+ end subroutine surf_temp
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine buildings(iurb,nd,nz,z0,ua_u,va_u,pt_u,pt0_u, &
+ ptg,ptr,da_u,ptw,ptwin,pwin, &
+ drst,uva_u,vva_u,uvb_u,vvb_u, &
+ tva_u,tvb_u,evb_u,qvb_u,qhb_u, &
+ uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr, &
+ sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)
+
+! ----------------------------------------------------------------------
+! 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,nbui_max) ! Walls potential temperatures
+ real ss(nz_um) ! probability to have a building with height h
+ real pb(nz_um)
+ real z0(ndm,nz_um) ! Roughness lengths "profiles"
+ real dt ! time step
+ integer iurb !Urban class
+
+!
+!New variables (BEM)
+!
+ real bs_u(ndm,nurbm) ! Building width [m]
+ real dz_u ! Urban grid resolution
+ real sflev(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems [W]
+ real lflev(nz_um,nz_um) ! latent heat flux due to the air conditioning systems [W]
+ real sfvlev(nz_um,nz_um) ! sensible heat flux due to ventilation [W]
+ real lfvlev(nz_um,nz_um) ! latent heat flux due to ventilation [W]
+ real qvb_u(2*ndm,nz_um)
+ real qhb_u(ndm,nz_um)
+ real ptwin(2*ndm,nz_um,nbui_max) ! window potential temperature
+ real pwin
+ real tvb_ac(2*ndm,nz_um)
+! ----------------------------------------------------------------------
+! 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
+ real sfw(2*ndm,nz_um,nbui_max) ! sensible heat flux from walls
+ real sfwin(2*ndm,nz_um,nbui_max) ! sensible heat flux form windows
+ real sfr(ndm,nz_um) ! sensible heat flux from roof
+ real sfg(ndm) ! sensible heat flux from street
+
+! ----------------------------------------------------------------------
+! LOCAL:
+! ----------------------------------------------------------------------
+ real d_urb(nz_um)
+ real uva_tmp
+ real vva_tmp
+ real uvb_tmp
+ real vvb_tmp
+ real evb_tmp
+ integer nlev(nz_um)
+ integer id,iz,ibui,nbui
+
+! ----------------------------------------------------------------------
+! END VARIABLES DEFINITIONS
+! ----------------------------------------------------------------------
+ dz=dz_u
+ ibui=0
+ nbui=0
+ nlev=0
+ d_urb=0.
+
+ uva_u=0.
+ uvb_u=0.
+ vhb_u=0.
+ vva_u=0.
+ vvb_u=0.
+ thb_u=0.
+ tva_u=0.
+ tvb_u=0.
+ tvb_ac=0.
+ ehb_u=0.
+ evb_u=0.
+ qvb_u=0.
+ qhb_u=0.
+
+ do iz=1,nz_um
+ if(ss(iz).gt.0) then
+ ibui=ibui+1
+ d_urb(ibui)=ss(iz)
+ nlev(ibui)=iz-1
+ nbui=ibui
+ endif
+ enddo
+ if (nbui.gt.nbui_max) then
+ write(*,*) 'nbui_max must be increased to',nbui
+ stop
+ endif
+ 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),sfg(id),ehb_u(id,1),da_u(1))
+ thb_u(id,1)=- sfg(id)/(da_u(1)*cp_u)
+
+! 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),sfr(id,iz),ehb_u(id,iz),da_u(iz))
+ thb_u(id,iz)=- sfr(id,iz)/(da_u(iz)*cp_u)
+ 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 ibui=1,nbui
+ do iz=1,nlev(ibui)
+
+ call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
+ ptw(2*id-1,iz,ibui),ptwin(2*id-1,iz,ibui), &
+ uva_tmp,vva_tmp, &
+ uvb_tmp,vvb_tmp, &
+ sfw(2*id-1,iz,ibui),sfwin(2*id-1,iz,ibui), &
+ evb_tmp,drst(id),dt)
+
+ if (pb(iz+1).gt.0.) then
+
+ uva_u(2*id-1,iz)=uva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
+ vva_u(2*id-1,iz)=vva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
+ uvb_u(2*id-1,iz)=uvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
+ vvb_u(2*id-1,iz)=vvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
+ evb_u(2*id-1,iz)=evb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
+ tvb_u(2*id-1,iz)=tvb_u(2*id-1,iz)-(d_urb(ibui)/pb(iz+1))* &
+ (sfw(2*id-1,iz,ibui)*(1.-pwin)+sfwin(2*id-1,iz,ibui)*pwin)/ &
+ da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
+ (dz*bs_u(id,iurb))/da_u(iz)/cp_u
+ tvb_ac(2*id-1,iz)=tvb_ac(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
+ (dz*bs_u(id,iurb))/da_u(iz)/cp_u
+ qvb_u(2*id-1,iz)=qvb_u(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
+ (dz*bs_u(id,iurb))/da_u(iz)/latent
+
+ endif
+
+ call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
+ ptw(2*id,iz,ibui),ptwin(2*id,iz,ibui), &
+ uva_tmp,vva_tmp, &
+ uvb_tmp,vvb_tmp, &
+ sfw(2*id,iz,ibui),sfwin(2*id,iz,ibui), &
+ evb_tmp,drst(id),dt)
+
+ if (pb(iz+1).gt.0.) then
+
+ uva_u(2*id,iz)=uva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
+ vva_u(2*id,iz)=vva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
+ uvb_u(2*id,iz)=uvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
+ vvb_u(2*id,iz)=vvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
+ evb_u(2*id,iz)=evb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
+ tvb_u(2*id,iz)=tvb_u(2*id,iz)-(d_urb(ibui)/pb(iz+1))* &
+ (sfw(2*id,iz,ibui)*(1.-pwin)+sfwin(2*id,iz,ibui)*pwin)/ &
+ da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
+ (dz*bs_u(id,iurb))/da_u(iz)/cp_u
+ tvb_ac(2*id,iz)=tvb_ac(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
+ (dz*bs_u(id,iurb))/da_u(iz)/cp_u
+ qvb_u(2*id,iz)=qvb_u(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
+ (dz*bs_u(id,iurb))/da_u(iz)/latent
+
+ endif
+!
+ enddo !iz
+ enddo !ibui
+
+ end do !id
+
+ 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, &
+ uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
+ uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u, &
+ a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)
+
+! ----------------------------------------------------------------------
+! 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 tvb_ac(2*ndm,nz_um)
+ 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
+!
+!New variables for BEM
+!
+ real qhb_u(ndm,nz_um)
+ real qvb_u(2*ndm,nz_um)
+
+
+! ----------------------------------------------------------------------
+! 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_ac(kms:kme)
+ real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
+ real b_q(kms:kme) ! Explicit component of the humidity sources or sinks
+! ----------------------------------------------------------------------
+! LOCAL:
+! ----------------------------------------------------------------------
+ real dzz
+ real fact
+ integer id,iz,iz_u
+ real se,sr,st,su,sv,sq
+ real uet(kms:kme) ! Contribution to TKE due to walls
+ real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb,vqb,vtb_ac
+
+
+! ----------------------------------------------------------------------
+! 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.
+ b_ac(iz)=0.
+ uet(iz)=0.
+ b_q(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))
+ b_q(kts)=b_q(kts)+qhb_u(id,1)*fact
+
+ do iz=kts,kte
+ st=0.
+ su=0.
+ sv=0.
+ se=0.
+ sq=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))
+ sq=sq+ss(iz_u)*qhb_u(id,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
+ b_q(iz)=b_q(iz)+sq*fact
+ enddo
+ enddo
+
+! vertical surface impact
+
+ do iz=kts,kte
+ uet(iz)=0.
+ do id=1,nd
+ vtb=0.
+ vtb_ac=0.
+ vta=0.
+ vua=0.
+ vub=0.
+ vva=0.
+ vvb=0.
+ veb=0.
+ vte=0.
+ vqb=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
+ vtb_ac=vtb_ac+pb(iz_u+1)* &
+ (tvb_ac(2*id-1,iz_u)+tvb_ac(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
+ vqb=vqb+pb(iz_u+1)* &
+ (qvb_u(2*id-1,iz_u)+qvb_u(2*id,iz_u))*fact
+ enddo
+
+ fact=1./vl(iz)/dz(iz)/nd
+ b_t(iz)=b_t(iz)+vtb*fact
+ b_ac(iz)=b_ac(iz)+vtb_ac*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
+ b_q(iz)=b_q(iz)+vqb*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, &
+ 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, &
+ 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)
+
+! ----------------------------------------------------------------------
+! 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,emwin,pwin,twlev,&
+ 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) ! Temperature in each layer of the wall [K]
+!
+!New Variables for BEM
+!
+ real twlev(2*ndm,nz_um) ! Window temperature in BEM [K]
+ real emwin ! Emissivity of windows
+ real pwin ! Coverage area fraction of windows in the walls of the buildings (BEM)
+
+
+! ----------------------------------------------------------------------
+! 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*(1.-pwin)-emwin*pwin)* &
+ fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
+ enddo
+
+ aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
+
+ bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg(id,ng_u)**4
+ do j=1,nz_u
+ bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i,id,iurb)* &
+ (emw*(1.-pwin)*tw(2*id,j)**4+emwin*pwin*twlev(2*id,j)**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*(1.-pwin)-emwin*pwin)*fww(j,i-nz_u,id,iurb)*pb(j+1)
+ enddo
+
+ do j=1+nz_u,2*nz_u
+ aaa(i,j)=0.
+ enddo
+
+ aaa(i,i)=1.
+
+ aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
+
+ bbb(i)=fsw(i-nz_u,id,iurb)*rl+ &
+ emg*fgw(i-nz_u,id,iurb)*sigma*tg(id,ng_u)**4
+
+ do j=1,nz_u
+ bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i-nz_u,id,iurb)* &
+ (emw*(1.-pwin)*tw(2*id-1,j)**4+emwin*pwin*twlev(2*id-1,j)**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*(1.-pwin)-emwin*pwin)* &
+ fwg(j,id,iurb)*pb(j+1)
+ enddo
+
+ do j=nz_u+1,2*nz_u
+ aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
+ 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)+sigma*fwg(i,id,iurb)*pb(i+1)* &
+ (emw*(1.-pwin)*(tw(2*id-1,i)**4+tw(2*id,i)**4)+ &
+ emwin*pwin*(twlev(2*id-1,i)**4+twlev(2*id,i)**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, &
+ 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)
+
+! ----------------------------------------------------------------------
+! 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
+ pause 'singular matrix in gaussj'
+ endif
+ enddo
+ endif
+ enddo
+
+ ipiv(icol)=ipiv(icol)+1
+
+ if(irow.ne.icol)then
+ do l=1,n
+ dum=a(irow,l)
+ a(irow,l)=a(icol,l)
+ a(icol,l)=dum
+ enddo
+
+ dum=b(irow)
+ b(irow)=b(icol)
+ b(icol)=dum
+
+ endif
+
+ if(a(icol,icol).eq.0)pause 'singular matrix in gaussj'
+
+ pivinv=1./a(icol,icol)
+ a(icol,icol)=1
+
+ do l=1,n
+ a(icol,l)=a(icol,l)*pivinv
+ enddo
+
+ b(icol)=b(icol)*pivinv
+
+ do ll=1,n
+ if(ll.ne.icol)then
+ dum=a(ll,icol)
+ a(ll,icol)=0.
+ do l=1,n
+ a(ll,l)=a(ll,l)-a(icol,l)*dum
+ enddo
+
+ b(ll)=b(ll)-b(icol)*dum
+
+ endif
+ enddo
+ enddo
+
+ return
+ end subroutine gaussj
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine soil_temp(nz,dz,temp,pt,ala,cs, &
+ 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
+
+ a(1,1)=0.
+ a(1,2)=1.
+ a(1,3)=0.
+ c(1)=temp(1)
+
+ do iz=2,nz-1
+ a(iz,1)=-cddz(iz)*dt/dz(iz)
+ a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
+ a(iz,3)=-cddz(iz+1)*dt/dz(iz)
+ c(iz)=temp(iz)
+ enddo
+
+ a(nz,1)=-dt*cddz(nz)/dz(nz)
+ a(nz,2)=1.+dt*cddz(nz)/dz(nz)
+ a(nz,3)=0.
+ c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
+
+
+ call invert(nz,a,c,temp)
+
+
+ pt=temp(nz)*(press/1.e+5)**(-rcp_u)
+
+ rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)
+
+ gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf
+ return
+ end subroutine soil_temp
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine invert(n,a,c,x)
+
+! ----------------------------------------------------------------------
+! 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,ptwin,uva,vva,uvb,vvb, &
+ sfw,sfwin,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 ptwin ! Windows 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
+ real sfw ! Surfaces fluxes from the walls
+ real sfwin ! Surfaces fluxes from the windows
+
+! LOCAL:
+! -----
+ real hc
+ real hcwin
+ 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
+
+ if (vett.lt.4.88) then
+ hc=5.678*(1.09+0.23*(vett/0.3048))
+ else
+ hc=5.678*0.53*((vett/0.3048)**0.78)
+ endif
+
+ if (hc.gt.da*cp_u/dt)then
+ hc=da*cp_u/dt
+ endif
+
+ if (vett.lt.4.88) then
+ hcwin=5.678*(0.99+0.21*(vett/0.3048))
+ else
+ hcwin=5.678*0.50*((vett/0.3048)**0.78)
+ endif
+
+ if (hcwin.gt.da*cp_u/dt) then
+ hcwin=da*cp_u/dt
+ endif
+
+! tvb=hc*ptw/da/cp_u
+! tva=-hc/da/cp_u
+!!!!!!!!!!!!!!!!!!!!
+! explicit
+
+ sfw=hc*(pt-ptw)
+ sfwin=hcwin*(pt-ptwin)
+
+
+ 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, &
+ uhb,vhb,sf,ehb,da)
+
+! ----------------------------------------------------------------------
+! Calculation of the flux at the ground
+! Formulation of Louis (Louis, 1979)
+! ----------------------------------------------------------------------
+
+ implicit none
+
+ 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
+ real da ! air density
+
+
+
+! ----------------------------------------------------------------------
+! 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
+ real sf
+
+
+! ----------------------------------------------------------------------
+! 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.
+
+
+ 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
+ sf= ustar*tstar*da*cp_u
+
+! thb= 0.
+ ehb=buu
+!!!!!!!!!!!!!!!
+
+ return
+ end subroutine flux_flat
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine icBEP (fww,fwg,fgw,fsw,fws,fsg, &
+ z0g_u,z0r_u, &
+ nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
+ nz_u,z_u)
+
+
+ implicit none
+
+
+! Building parameters
+
+! Radiation parameters
+
+! 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 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
+! -----------------------------------------------------------------------
+! 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
+
+ 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 some variables
+!
+ nz_u=0
+ z_u=0.
+ ss_u=0.
+ pb_u=0.
+ fww=0.
+ fwg=0.
+ fgw=0.
+ fsw=0.
+ fws=0.
+ fsg=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)
+
+ call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb), &
+ z_u,ws_u(id,iurb), &
+ fww,fwg,fgw,fsg,fsw,fws)
+
+ do iz_u=1,nz_u(iurb)
+ ss_u(iz_u,iurb)=0.
+ do ilu=1,nz_um
+ if(z_u(iz_u).le.h_b(ilu,iurb) &
+ .and.z_u(iz_u+1).gt.h_b(ilu,iurb))then
+ ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
+ endif
+ enddo
+ enddo
+
+ pb_u(1,iurb)=1.
+ do iz_u=1,nz_u(iurb)
+ pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
+ enddo
+
+ enddo
+ end do
+
+
+ return
+ end subroutine icBEP
+
+! ===6=8===============================================================72
+! ===6=8===============================================================72
+
+ subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws)
+
+ 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)
+
+ 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)
+
+ 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,&
+ twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
+ emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b, &
+ cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
+ targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)
+
+! 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 albwin_u(nurbm) ! Albedo of the window
+ real emg_u(nurbm) ! Emissivity of ground
+ real emw_u(nurbm) ! Emissivity of wall
+ real emr_u(nurbm) ! Emissivity of roof
+ real emwind_u(nurbm) ! Emissivity of windows
+
+! 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
+
+ real, intent(out) :: cop_u(nurbm)
+ real, intent(out) :: pwin_u(nurbm)
+ real, intent(out) :: beta_u(nurbm)
+ integer, intent(out) :: sw_cond_u(nurbm)
+ real, intent(out) :: time_on_u(nurbm)
+ real, intent(out) :: time_off_u(nurbm)
+ real, intent(out) :: targtemp_u(nurbm)
+ real, intent(out) :: gaptemp_u(nurbm)
+ real, intent(out) :: targhum_u(nurbm)
+ real, intent(out) :: gaphum_u(nurbm)
+ real, intent(out) :: perflo_u(nurbm)
+ real, intent(out) :: hsesf_u(nurbm)
+ real, intent(out) :: hsequip(24)
+
+!
+!We initialize
+!
+ 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
+!MT BEM
+ cop_u = cop_tbl
+ pwin_u = pwin_tbl
+ beta_u = beta_tbl
+ sw_cond_u = sw_cond_tbl
+ time_on_u = time_on_tbl
+ time_off_u = time_off_tbl
+ targtemp_u = targtemp_tbl
+ gaptemp_u = gaptemp_tbl
+ targhum_u = targhum_tbl
+ gaphum_u = gaphum_tbl
+ perflo_u = perflo_tbl
+ hsesf_u = hsesf_tbl
+ hsequip = hsequip_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
+
+ do iu=1,nurb
+ emwind_u(iu)=0.9
+ call albwindow(albwin_u(iu))
+ enddo
+
+ return
+ end subroutine init_para
+! ===6================================================================72
+! ===6================================================================72
+ subroutine upward_rad(nd_u,nz_u,ws,bs,sigma,pb,ss, &
+ tg,emg_u,albg_u,rlg,rsg,sfg, &
+ tw,emw_u,albw_u,rlw,rsw,sfw, &
+ tr,emr_u,albr_u,emwind,albwind,twlev,pwin, &
+ sfwind,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 nz_u
+ 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) ! 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 nd_u ! number of street directions
+!
+!New variables BEM
+!
+ real emwind !Emissivity of the windows
+ real albwind !Albedo of the windows
+ real twlev(2*ndm,nz_um) !Averaged Temperature of the windows
+ real pwin !Coverage area fraction of the windows
+ real gflwin !Heat stored for the windows
+ real sfwind(2*ndm,nz_um) !Sensible heat flux from windows [W/m²]
+
+!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,nz_u+1
+ do id=1,nd_u
+ do iw=1,nwr_u
+ if(tr(id,iz,iw).lt.100.)then
+ write(*,*)'in upward_rad ',iz,id,iw,tr(id,iz,iw)
+ iwrong=0
+ endif
+ enddo
+ enddo
+ enddo
+ if(iwrong.eq.0)stop
+
+ rl_up=0.
+
+ rs_abs=0.
+ rl_inc=0.
+ emiss=0.
+ rl_emit=0.
+ grdflx_urb=0.
+ do id=1,nd_u
+ rl_emit=rl_emit-( emg_u*sigma*(tg(id,ng_u)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/nd_u
+ rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u
+ rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/nd_u
+ gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg(id,ng_u)**4.)+sfg(id)
+ grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u
+
+ do iz=2,nz_u
+ rl_emit=rl_emit-(emr_u*sigma*(tr(id,iz,nwr_u)**4.)+(1-emr_u)*rld)*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
+ rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
+ rs_abs=rs_abs+(1.-albr_u)*rs*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
+ gfl=(1.-albr_u)*rs+emr_u*rld-emr_u*sigma*(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
+ grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
+ enddo
+
+ do iz=1,nz_u
+
+ rl_emit=rl_emit-(emw_u*(1.-pwin)*sigma*(tw(2*id-1,iz)**4.+tw(2*id,iz)**4.)+ &
+ (emwind*pwin*sigma*(twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.))+ &
+ ((1.-emw_u)*(1.-pwin)+pwin*(1.-emwind))*(rlw(2*id-1,iz)+rlw(2*id,iz)))* &
+ dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
+
+ rl_inc=rl_inc+((rlw(2*id-1,iz)+rlw(2*id,iz)))*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
+
+ rs_abs=rs_abs+(((1.-albw_u)*(1.-pwin)+(1.-albwind)*pwin)*(rsw(2*id-1,iz)+rsw(2*id,iz)))*&
+ dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
+
+ gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) ) &
+ -emw_u*sigma*( tw(2*id-1,iz)**4.+tw(2*id,iz)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))
+
+ gflwin=(1.-albwind)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emwind*(rlw(2*id-1,iz)+rlw(2*id,iz)) &
+ -emwind*sigma*( twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.)+(sfwind(2*id-1,iz)+sfwind(2*id,iz))
+
+
+ grdflx_urb=grdflx_urb-(gfl*(1.-pwin)+pwin*gflwin)*dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
+
+ enddo
+
+ enddo
+ emiss=(emg_u+emw_u+emr_u)/3.
+ rl_up=(rl_inc+rl_emit)-rld
+
+
+ return
+
+ END SUBROUTINE upward_rad
+
+!====6=8===============================================================72
+!====6=8===============================================================72
+! ===6================================================================72
+! ===6================================================================72
+
+ subroutine albwindow(albwin)
+
+!-------------------------------------------------------------------
+ implicit none
+
+
+! -------------------------------------------------------------------
+!Based on the
+!paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour
+!of the total solar energy transmittance of windows"
+!Solar Energy Vol.69,No.4,pp. 321-329.
+! -------------------------------------------------------------------
+!Input
+!-----
+
+!Output
+!------
+ real albwin ! albedo of the window
+!Local
+!-----
+ real a,b,c !Polynomial coefficients
+ real alfa,delta,gama !Polynomial powers
+ real g0 !transmittance when the angle
+ !of incidence is normal to the surface.
+ real asup,ainf
+ real fonc
+
+!Constants
+!--------------------
+
+ real epsilon !accuracy of the integration
+ parameter (epsilon=1.e-07)
+ real n1,n2 !Index of refraction for glasses and air
+ parameter(n1=1.,n2=1.5)
+ integer intg,k
+!--------------------------------------------------------------------
+ if (q_num.eq.0) then
+ write(*,*) 'Category parameter of the windows no valid'
+ stop
+ endif
+
+ g0=4.*n1*n2/((n1+n2)*(n1+n2))
+ a=8.
+ b=0.25/q_num
+ c=1.-a-b
+ alfa =5.2 + (0.7*q_num)
+ delta =2.
+ gama =(5.26+0.06*p_num)+(0.73+0.04*p_num)*q_num
+
+ intg=1
+!----------------------------------------------------------------------
+
+
+100 asup=0.
+ ainf=0.
+
+ do k=1,intg
+ call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
+ asup=asup+(pi/intg)*fonc
+ enddo
+
+ intg=intg+1
+
+ do k=1,intg
+ call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
+ ainf=ainf+(pi/intg)*fonc
+ enddo
+
+ if(abs(asup-ainf).lt.epsilon) then
+ albwin=1-g0+(g0/2.)*asup
+ else
+ goto 100
+ endif
+
+!----------------------------------------------------------------------
+ return
+ end subroutine albwindow
+!====================================================================72
+!====================================================================72
+
+ subroutine foncs(fonc,x,aa,bb,cc,alf,delt,gam)
+
+ implicit none
+!
+ real x,aa,bb,cc
+ real alf,delt,gam
+ real fonc
+
+ fonc=(((aa*(x**alf))/(pi**alf))+ &
+ ((bb*(x**delt))/(pi**delt))+ &
+ ((cc*(x**gam))/(pi**gam)))*sin(x)
+
+ return
+ end subroutine foncs
+!====================================================================72
+!====================================================================72
+END MODULE module_sf_bep_bem
</font>
</pre>