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