| blob | f249700a08eb72fa2c844e95da7229d4dc2c432f |
1 MODULE module_sf_bep_bem
2 #if defined(mpas)
3 use mpas_atmphys_utilities, only: physics_error_fatal
4 #define FATAL_ERROR(M) call physics_error_fatal( M )
5 #else
6 use module_wrf_error
7 #define FATAL_ERROR(M) call wrf_error_fatal( M )
8 !USE module_model_constants
9 #endif
10 USE module_sf_urban
11 USE module_sf_bem
12 USE module_bep_bem_helper, ONLY: nurbm
14 ! SGClarke 09/11/2008
15 ! Access urban_param.tbl values through calling urban_param_init in module_physics_init
16 ! for CASE (BEPSCHEME) select sf_urban_physics
17 !
18 ! -----------------------------------------------------------------------
19 ! Dimension for the array used in the BEP module
20 ! -----------------------------------------------------------------------
22 integer nurbmax ! Maximum number of urban classes
23 parameter (nurbmax=11)
25 integer ndm ! Maximum number of street directions
26 parameter (ndm=2)
28 integer nz_um ! Maximum number of vertical levels in the urban grid
29 parameter(nz_um=18)
31 integer ng_u ! Number of grid levels in the ground
32 parameter (ng_u=10)
34 integer ngr_u ! Number of grid levels in green roof
35 parameter (ngr_u=10)
37 integer nwr_u ! Number of grid levels in the walls or roofs
38 parameter (nwr_u=10)
40 integer nf_u !Number of grid levels in the floors (BEM)
41 parameter (nf_u=10)
43 integer ngb_u !Number of grid levels in the ground below building (BEM)
44 parameter (ngb_u=10)
46 real dz_u ! Urban grid resolution
47 parameter (dz_u=5.)
49 integer nbui_max !maximum number of types of buildings in an urban class
50 parameter (nbui_max=15) !must be less or equal than nz_um
53 real h_water
54 parameter(h_water=0.0009722) !mm of irrigation per hour
56 !---------------------------------------------------------------------------------
57 !Parameters of the windows. The glasses of windows are considered without films -
58 !Read the paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour -
59 !of the total solar energy transmittance of windows".Solar Energy Vol.69,No.4, -
60 !pp. 321-329, for more details. -
61 !---------------------------------------------------------------------------------
62 integer p_num !number of panes in the windows (1,2 or 3)
63 parameter (p_num=2)
64 integer q_num !category number for the windows (q_num= 4, standard glasses)
65 parameter(q_num=4) !Possible values 1,2,...,10
67 ! The change of ng_u, nwr_u should be done in agreement with the block data
68 ! in the routine "surf_temp"
69 ! -----------------------------------------------------------------------
70 ! Constant used in the BEP module
71 ! -----------------------------------------------------------------------
73 real vk ! von Karman constant
74 real g_u ! Gravity acceleration
75 real pi !
76 real r ! Perfect gas constant
77 real cp_u ! Specific heat at constant pressure
78 real rcp_u !
79 real sigma !
80 real p0 ! Reference pressure at the sea level
81 real latent ! Latent heat of vaporization [J/kg] (used in BEM)
82 real dgmax ! Maximum ground water holding capacity (mm)
83 real drmax ! Maximum ground roof holding capacity (mm)
85 parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)
86 parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,latent=2.45e+06,dgmax=1.,drmax=1.)
88 ! -----------------------------------------------------------------------
93 CONTAINS
95 subroutine BEP_BEM(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, &
96 th_phy,rho,p_phy,swdown,glw, &
97 gmt,julday,xlong,xlat, &
98 declin_urb,cosz_urb2d,omg_urb2d, &
99 num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &
100 urban_map_zd, urban_map_zdf, urban_map_bd, urban_map_wd, &
101 urban_map_gbd, urban_map_fbd, &
102 urban_map_zgrd, num_urban_hi, &
103 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
104 tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
105 tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
106 cm_ac_urb3d, &
107 sfvent_urb3d,lfvent_urb3d, &
108 sfwin1_urb3d,sfwin2_urb3d, &
109 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
110 ep_pv_urb3d,t_pv_urb3d, &
111 trv_urb4d,qr_urb4d,qgr_urb3d,tgr_urb3d, &
112 drain_urb4d,draingr_urb3d, &
113 sfrv_urb3d,lfrv_urb3d, &
114 dgr_urb3d,dg_urb3d, &
115 lfr_urb3d,lfg_urb3d,rainbl,swddir,swddif, &
116 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
117 a_u,a_v,a_t,a_e,b_u,b_v, &
118 b_t,b_e,b_q,dlg,dl_u,sf,vl, &
119 rl_up,rs_abs,emiss,grdflx_urb,qv_phy, &
120 ids,ide, jds,jde, kds,kde, &
121 ims,ime, jms,jme, kms,kme, &
122 its,ite, jts,jte, kts,kte)
124 implicit none
126 !------------------------------------------------------------------------
127 ! Input
128 !------------------------------------------------------------------------
129 INTEGER :: ids,ide, jds,jde, kds,kde, &
130 ims,ime, jms,jme, kms,kme, &
131 its,ite, jts,jte, kts,kte, &
132 itimestep
135 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: DZ8W
136 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: P_PHY
137 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: RHO
138 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: TH_PHY
139 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: T_PHY
140 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U_PHY
141 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V_PHY
142 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U
143 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V
144 REAL, DIMENSION( ims:ime , jms:jme ) :: GLW
145 REAL, DIMENSION( ims:ime , jms:jme ) :: swdown
146 REAL, DIMENSION( ims:ime , jms:jme ) :: swddir
147 REAL, DIMENSION( ims:ime , jms:jme ) :: swddif
148 REAL, DIMENSION( ims:ime, jms:jme ) :: UST
149 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: UTYPE_URB2D
150 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: FRC_URB2D
151 REAL, INTENT(IN ) :: GMT
152 INTEGER, INTENT(IN ) :: JULDAY
153 REAL, DIMENSION( ims:ime, jms:jme ), &
154 INTENT(IN ) :: XLAT, XLONG
155 REAL, INTENT(IN) :: DECLIN_URB
156 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
157 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
158 INTEGER, INTENT(IN ) :: urban_map_zrd
159 INTEGER, INTENT(IN ) :: urban_map_zwd
160 INTEGER, INTENT(IN ) :: urban_map_gd
161 INTEGER, INTENT(IN ) :: urban_map_zd
162 INTEGER, INTENT(IN ) :: urban_map_zdf
163 INTEGER, INTENT(IN ) :: urban_map_bd
164 INTEGER, INTENT(IN ) :: urban_map_wd
165 INTEGER, INTENT(IN ) :: urban_map_gbd
166 INTEGER, INTENT(IN ) :: urban_map_fbd
167 INTEGER, INTENT(IN ) :: num_urban_ndm
168 INTEGER, INTENT(IN) :: num_urban_hi
169 INTEGER , INTENT(IN) :: urban_map_zgrd
170 REAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
171 REAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
172 REAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
173 REAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
174 REAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ), INTENT(INOUT) :: trv_urb4d
175 REAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ), INTENT(INOUT) :: qr_urb4d
176 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: qgr_urb3d
177 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tgr_urb3d
178 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: drain_urb4d
179 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: rainbl
180 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: draingr_urb3d
181 !New variables used for BEM
182 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: qv_phy
183 REAL, DIMENSION( ims:ime, 1:urban_map_bd, jms:jme ), INTENT(INOUT) :: tlev_urb3d
184 REAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: qlev_urb3d
185 REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
186 REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
187 REAL, DIMENSION( ims:ime, 1:urban_map_gbd, jms:jme ), INTENT(INOUT) :: tglev_urb3d
188 REAL, DIMENSION( ims:ime, 1:urban_map_fbd, jms:jme ), INTENT(INOUT) :: tflev_urb3d
189 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
190 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
191 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
192 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ep_pv_urb3d
193 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: t_pv_urb3d
194 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
195 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
196 REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
197 REAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
199 !End variables
200 REAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
201 REAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
202 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
203 REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
204 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfrv_urb3d
205 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: lfrv_urb3d
206 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: dgr_urb3d !GRZ
207 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: dg_urb3d !GRZ
208 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfr_urb3d !GRZ
209 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: lfg_urb3d !G
211 REAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
212 REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lp_urb2d
213 REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: lb_urb2d
214 REAL, DIMENSION( ims:ime,jms:jme), INTENT(IN) :: hgt_urb2d
216 real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates
217 REAL, INTENT(IN ):: DT ! Time step
219 !------------------------------------------------------------------------
220 ! Output
221 !------------------------------------------------------------------------
222 !
223 ! Implicit and explicit components of the source and sink terms at each levels,
224 ! the fluxes can be computed as follow: FX = A*X + B example: t_fluxes = a_t * pt + b_t
225 real a_u(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in X-direction (center)
226 real a_v(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in Y-direction (center)
227 real a_t(ims:ime,kms:kme,jms:jme) ! Implicit component for the temperature
228 real a_e(ims:ime,kms:kme,jms:jme) ! Implicit component for the TKE
229 real b_u(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in X-direction (center)
230 real b_v(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in Y-direction (center)
231 real b_t(ims:ime,kms:kme,jms:jme) ! Explicit component for the temperature
232 real b_e(ims:ime,kms:kme,jms:jme) ! Explicit component for the TKE
233 real b_q(ims:ime,kms:kme,jms:jme) ! Explicit component for the Humidity
234 real dlg(ims:ime,kms:kme,jms:jme) ! Height above ground (L_ground in formula (24) of the BLM paper).
235 real dl_u(ims:ime,kms:kme,jms:jme) ! Length scale (lb in formula (22) ofthe BLM paper).
236 ! urban surface and volumes
237 real sf(ims:ime,kms:kme,jms:jme) ! surface of the urban grid cells
238 real vl(ims:ime,kms:kme,jms:jme) ! volume of the urban grid cells
239 ! urban fluxes
240 real rl_up(its:ite,jts:jte) ! upward long wave radiation
241 real rs_abs(its:ite,jts:jte) ! absorbed short wave radiation
242 real emiss(its:ite,jts:jte) ! emissivity averaged for urban surfaces
243 real grdflx_urb(its:ite,jts:jte) ! ground heat flux for urban areas
244 !------------------------------------------------------------------------
245 ! Local
246 !------------------------------------------------------------------------
247 real hi_urb(its:ite,1:nz_um,jts:jte) ! Height histograms of buildings
248 real hi_urb1D(nz_um) ! Height histograms of buildings
249 real ss_urb(nz_um,nurbmax) ! Probability that a building has an height equal to z
250 real pb_urb(nz_um) ! Probability that a building has an height greater or equal to z
251 real hb_u(nz_um) ! Bulding's heights
252 integer nz_urb(nurbmax) ! Number of layer in the urban grid
253 integer nzurban(nurbmax)
255 ! Building parameters
256 real alag_u(nurbmax) ! Ground thermal diffusivity [m^2 s^-1]
257 real alaw_u(nurbmax) ! Wall thermal diffusivity [m^2 s^-1]
258 real alar_u(nurbmax) ! Roof thermal diffusivity [m^2 s^-1]
259 real csg_u(nurbmax) ! Specific heat of the ground material [J m^3 K^-1]
260 real csw_u(nurbmax) ! Specific heat of the wall material [J m^3 K^-1]
261 real csr_u(nurbmax) ! Specific heat of the roof material [J m^3 K^-1]
262 real twini_u(nurbmax) ! Initial temperature inside the building's wall [K]
263 real trini_u(nurbmax) ! Initial temperature inside the building's roof [K]
264 real tgini_u(nurbmax) ! Initial road temperature
266 !
267 ! Building materials
268 !
270 real csg(ng_u) ! Specific heat of the ground material [J m^3 K^-1]
271 real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
272 real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
273 real csgb(ngb_u) ! Specific heat of the ground material below the buildings at each ground levels[J m^3 K^-1]
274 real csf(nf_u) ! Specific heat of the floors materials in the buildings at each levels[J m^3 K^-1]
275 real alar(nwr_u+1) ! Roof thermal diffusivity for the current urban class [W/m K]
276 real alaw(nwr_u+1) ! Walls thermal diffusivity for the current urban class [W/m K]
277 real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
278 real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall layer [W/m K]
279 real alaf(nf_u+1) ! Floor thermal diffusivity at each wall layers [W/m K]
280 real dzr(nwr_u) ! Layer sizes in the roofs [m]
281 real dzf(nf_u) ! Layer sizes in the floors[m]
282 real dzw(nwr_u) ! Layer sizes in the walls [m]
283 real dzgb(ngb_u) ! Layer sizes in the ground below the buildings [m]
285 !
286 !New street and radiation parameters
289 real bs(ndm) ! Building width for the current urban class
290 real ws(ndm) ! Street widths of the current urban class
291 real strd(ndm) ! Street lengths for the current urban class
292 real drst(ndm) ! street directions for the current urban class
293 real ss(nz_um) ! Probability to have a building with height h
294 real pb(nz_um) ! Probability to have a building with an height equal
295 real HFGR_D(nz_um)
296 !New roughness and buildings parameters
297 !
298 real z0(ndm,nz_um) ! Roughness lengths "profiles"
299 real bs_urb(ndm,nurbmax) ! Building width
300 real ws_urb(ndm,nurbmax) ! Street width
302 !
303 ! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
304 !
305 ! Radiation paramters
306 real albg_u(nurbmax) ! Albedo of the ground
307 real albw_u(nurbmax) ! Albedo of the wall
308 real albr_u(nurbmax) ! Albedo of the roof
309 real albwin_u(nurbmax) ! Albedo of the windows
310 real emwind_u(nurbmax) ! Emissivity of windows
311 real emg_u(nurbmax) ! Emissivity of ground
312 real emw_u(nurbmax) ! Emissivity of wall
313 real emr_u(nurbmax) ! Emissivity of roof
314 real gr_frac_roof_u(nurbmax)
315 real pv_frac_roof_u(nurbmax)
316 integer gr_flag_u
317 integer gr_type_u
319 ! fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
320 ! and the short wave radiation.
321 real fww_u(nz_um,nz_um,ndm,nurbmax) ! from wall to wall
322 real fwg_u(nz_um,ndm,nurbmax) ! from wall to ground
323 real fgw_u(nz_um,ndm,nurbmax) ! from ground to wall
324 real fsw_u(nz_um,ndm,nurbmax) ! from sky to wall
325 real fws_u(nz_um,ndm,nurbmax) ! from sky to wall
326 real fsg_u(ndm,nurbmax) ! from sky to ground
328 ! Roughness parameters
329 real z0g_u(nurbmax) ! The ground's roughness length
330 real z0r_u(nurbmax) ! The roof's roughness length
332 ! Street parameters
333 integer nd_u(nurbmax) ! Number of street direction for each urban class
334 real strd_u(ndm,nurbmax) ! Street length (fix to greater value to the horizontal length of the cells)
335 real drst_u(ndm,nurbmax) ! Street direction
336 real ws_u(ndm,nurbmax) ! Street width
337 real bs_u(ndm,nurbmax) ! Building width
338 real h_b(nz_um,nurbmax) ! Bulding's heights
339 real d_b(nz_um,nurbmax) ! Probability that a building has an height h_b
340 real ss_u(nz_um,nurbmax)! Probability that a building has an height equal to z
341 real pb_u(nz_um,nurbmax)! Probability that a building has an height greater or equal to z
344 ! Grid parameters
345 integer nz_u(nurbmax) ! Number of layer in the urban grid
347 real z_u(nz_um) ! Height of the urban grid levels
348 !FS
349 real cop_u(nurbmax)
350 real bldac_frc_u(nurbmax)
351 real cooled_frc_u(nurbmax)
352 real pwin_u(nurbmax)
353 real beta_u(nurbmax)
354 integer sw_cond_u(nurbmax)
355 real time_on_u(nurbmax)
356 real time_off_u(nurbmax)
357 real targtemp_u(nurbmax)
358 real gaptemp_u(nurbmax)
359 real targhum_u(nurbmax)
360 real gaphum_u(nurbmax)
361 real perflo_u(nurbmax)
362 real hsesf_u(nurbmax)
363 real hsequip(24)
364 real irho(24)
365 ! 1D array used for the input and output of the routine "urban"
367 real z1D(kms:kme) ! vertical coordinates
368 real ua1D(kms:kme) ! wind speed in the x directions
369 real va1D(kms:kme) ! wind speed in the y directions
370 real pt1D(kms:kme) ! potential temperature
371 real da1D(kms:kme) ! air density
372 real pr1D(kms:kme) ! air pressure
373 real pt01D(kms:kme) ! reference potential temperature
374 real zr1D ! zenith angle
375 real deltar1D ! declination of the sun
376 real ah1D ! hour angle (it should come from the radiation routine)
377 real rs1D ! solar radiation
378 real rld1D ! downward flux of the longwave radiation
379 real swddir1D
380 real swddif1D ! short wave diffuse solar radiation _gl
384 real tw1D(2*ndm,nz_um,nwr_u,nbui_max) ! temperature in each layer of the wall
385 real tg1D(ndm,ng_u) ! temperature in each layer of the ground
386 real tr1D(ndm,nz_um,nwr_u) ! temperature in each layer of the roof
387 real trv1D(ndm,nz_um,ngr_u) ! temperature in each layer of the GREEN roof
388 real qr1D(ndm,nz_um,ngr_u) ! humidity in each layer of the GREEN roof
390 !
391 !New variable for BEM
392 !
393 real tlev1D(nz_um,nbui_max) ! temperature in each floor and in each different type of building
394 real qlev1D(nz_um,nbui_max) ! specific humidity in each floor and in each different type of building
395 real twlev1D(2*ndm,nz_um,nbui_max) ! temperature in each window in each floor in each different type of building
396 real tglev1D(ndm,ngb_u,nbui_max) ! temperature in each layer of the ground below of a type of building
397 real tflev1D(ndm,nf_u,nz_um-1,nbui_max)! temperature in each layer of the floors in each building
398 real lflev1D(nz_um,nz_um) ! latent heat flux due to the air conditioning systems
399 real sflev1D(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems
400 real lfvlev1D(nz_um,nz_um) ! latent heat flux due to ventilation
401 real sfvlev1D(nz_um,nz_um) ! sensible heat flux due to ventilation
402 real sfwin1D(2*ndm,nz_um,nbui_max) ! sensible heat flux from windows
403 real consumlev1D(nz_um,nz_um) ! consumption due to the air conditioning systems
404 real eppvlev1D(nz_um) ! electricity production of PV panels
405 real tair1D(nz_um)
406 real tpvlev1D(ndm,nz_um)
407 real qv1D(kms:kme) ! specific humidity
408 real meso_urb ! constant to link meso and urban scales [m-2]
409 real meso_urb_ac
410 real roof_frac ! Surface fraction occupied by roof
411 real d_urb(nz_um)
412 real sf_ac
413 integer ibui,nbui
414 integer nlev(nz_um)
416 !
417 !End new variables
418 !
420 real sfw1D(2*ndm,nz_um,nbui_max) ! sensible heat flux from walls
421 real sfg1D(ndm) ! sensible heat flux from ground (road)
422 real sfr1D(ndm,nz_um) ! sensible heat flux from roofs
423 real sfrpv1D(ndm,nz_um)
425 real tpv1D(nbui_max)
426 real sfr_indoor1D(nbui_max)
427 real sfrv1D(ndm,nz_um) ! sensible heat flux from roofs
428 real lfrv1D(ndm,nz_um) ! latent heat flux from roofs
429 real dg1D(ndm) ! water depth from ground
430 real dgr1D(ndm,nz_um) ! water depth from roofs
431 real lfg1D(ndm) ! latent heat flux from ground (road)
432 real lfr1D(ndm,nz_um) ! latent heat flux from roofs
433 real drain1D(ndm,nz_um) ! sensible heat flux from roofs
434 real sf1D(kms:kme) ! surface of the urban grid cells
435 real vl1D(kms:kme) ! volume of the urban grid cells
436 real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
437 real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
438 real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
439 real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
440 real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
441 real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
442 real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
443 real b_ac1D(kms:kme)
444 real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
445 real b_q1D(kms:kme) ! Explicit component of the Humidity sources or sinks
446 real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
447 real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
448 real gfr1D(ndm,nz_um)
449 real time_bep
450 ! arrays used to collapse indexes
451 integer ind_zwd(nbui_max,nz_um,nwr_u,ndm)
452 integer ind_gd(ng_u,ndm)
453 integer ind_zd(nbui_max,nz_um,ndm)
454 integer ind_zdf(nz_um,ndm)
455 integer ind_zrd(nz_um,nwr_u,ndm)
456 integer ind_grd(nz_um,ngr_u,ndm)
457 !
458 integer ind_bd(nbui_max,nz_um)
459 integer ind_wd(nbui_max,nz_um,ndm)
460 integer ind_gbd(nbui_max,ngb_u,ndm)
461 integer ind_fbd(nbui_max,nf_u,nz_um-1,ndm)
463 integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
464 integer it, nint
465 integer iii
466 logical first
467 character(len=80) :: text
468 data first/.true./
469 save first,time_bep
471 save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
472 albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
473 z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
474 nz_u,z_u,albwin_u,emwind_u,cop_u,pwin_u,beta_u,sw_cond_u, &
475 bldac_frc_u,cooled_frc_u, &
476 time_on_u,time_off_u,targtemp_u,gaptemp_u,targhum_u,gaphum_u, &
477 perflo_u,gr_frac_roof_u, &
478 pv_frac_roof_u,hsesf_u,hsequip,irho,gr_flag_u,gr_type_u
480 !------------------------------------------------------------------------
481 ! Calculation of the momentum, heat and turbulent kinetic fluxes
482 ! produced by buildings
483 !
484 ! References:
485 ! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
486 ! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
487 ! 261-304
488 !
489 ! F. Salamanca and A. Martilli, 2009: 'A new Building Energy Model coupled
490 ! with an Urban Canopy Parameterization for urban climate simulations - part II.
491 ! Validation with one dimension off-line simulations'. Theor Appl Climatol
492 ! DOI 10.1007/s00704-009-0143-8
493 !------------------------------------------------------------------------
494 !
495 !prepare the arrays to collapse indexes
498 !
499 if(urban_map_zwd.lt.nbui_max*nz_um*ndm*max(nwr_u,ng_u))then
500 write(*,*)'urban_map_zwd too small, please increase to at least ', nbui_max*nz_um*ndm*max(nwr_u,ng_u)
501 stop
502 endif
503 !
504 !New conditions for BEM
505 !
506 if(urban_map_bd.lt.nbui_max*nz_um)then !limit for indoor temperature and indoor humidity
507 write(*,*)'urban_map_bd too small, please increase to at least ', nbui_max*nz_um
508 stop
509 endif
511 if(urban_map_wd.lt.nbui_max*nz_um*ndm)then !limit for window temperature
512 write(*,*)'urban_map_wd too small, please increase to at least ', nbui_max*nz_um*ndm
513 stop
514 endif
516 if(urban_map_gbd.lt.nbui_max*ndm*ngb_u)then !limit for ground temperature below a building
517 write(*,*)'urban_map_gbd too small, please increase to at least ', nbui_max*ndm*ngb_u
518 stop
519 endif
521 if(urban_map_fbd.lt.(nz_um-1)*nbui_max*ndm*nf_u)then !limit for floor temperature
522 write(*,*)'urban_map_fbd too small, please increase to at least ', nbui_max*ndm*nf_u*(nz_um-1)
523 stop
524 endif
526 if (ndm.ne.2)then
527 write(*,*) 'number of directions is not correct',ndm
528 stop
529 endif
531 !End of new conditions
532 !
533 !
534 !Initialize collapse indexes
535 !
536 ind_zwd=0
537 ind_gd=0
538 ind_zd=0
539 ind_zdf=0
540 ind_zrd=0
541 ind_grd=0
542 ind_bd=0
543 ind_wd=0
544 ind_gbd=0
545 ind_fbd=0
546 !
547 !End initialization indexes
548 !
550 iii=0
551 do ibui=1,nbui_max
552 do iz_u=1,nz_um
553 do iw=1,nwr_u
554 do id=1,ndm
555 iii=iii+1
556 ind_zwd(ibui,iz_u,iw,id)=iii
557 enddo
558 enddo
559 enddo
560 enddo
562 iii=0
563 do ig=1,ng_u
564 do id=1,ndm
565 iii=iii+1
566 ind_gd(ig,id)=iii
567 enddo
568 enddo
570 iii=0
571 do ibui=1,nbui_max
572 do iz_u=1,nz_um
573 do id=1,ndm
574 iii=iii+1
575 ind_zd(ibui,iz_u,id)=iii
576 enddo
577 enddo
578 enddo
580 iii=0
581 do iz_u=1,nz_um
582 do iw=1,nwr_u
583 do id=1,ndm
584 iii=iii+1
585 ind_zrd(iz_u,iw,id)=iii
586 enddo
587 enddo
588 enddo
590 iii=0
591 do iz_u=1,nz_um
592 do iw=1,ngr_u
593 do id=1,ndm
594 iii=iii+1
595 ind_grd(iz_u,iw,id)=iii
596 enddo
597 enddo
598 enddo
600 !
601 !New indexes for BEM
603 iii=0
604 do iz_u=1,nz_um
605 do id=1,ndm
606 iii=iii+1
607 ind_zdf(iz_u,id)=iii
608 enddo ! id
609 enddo ! iz_u
611 iii=0
612 do ibui=1,nbui_max !Type of building
613 do iz_u=1,nz_um !vertical levels
614 iii=iii+1
615 ind_bd(ibui,iz_u)=iii
616 enddo !iz_u
617 enddo !ibui
621 iii=0
622 do ibui=1,nbui_max !type of building
623 do iz_u=1,nz_um !vertical levels
624 do id=1,ndm !direction
625 iii=iii+1
626 ind_wd(ibui,iz_u,id)=iii
627 enddo !id
628 enddo !iz_u
629 enddo !ibui
631 iii=0
632 do ibui=1,nbui_max!type of building
633 do iw=1,ngb_u !layers in the wall (ground below a building)
634 do id=1,ndm !direction
635 iii=iii+1
636 ind_gbd(ibui,iw,id)=iii
637 enddo !id
638 enddo !iw
639 enddo !ibui
641 iii=0
642 do ibui=1,nbui_max !type of building
643 do iw=1,nf_u !layers in the wall (floor)
644 do iz_u=1,nz_um-1 !vertical levels
645 do id=1,ndm !direction
646 iii=iii+1
647 ind_fbd(ibui,iw,iz_u,id)=iii
648 enddo !id
649 enddo !iz_u
650 enddo !iw
651 enddo !ibui
654 !End of new indexes
656 if (num_urban_hi.ge.nz_um)then
657 write(*,*)'nz_um too small, please increase to at least ', num_urban_hi+1
658 stop
659 endif
661 do ix=its,ite
662 do iy=jts,jte
663 do iz_u=1,nz_um
664 hi_urb(ix,iz_u,iy)=0.
665 enddo
666 enddo
667 enddo
669 do ix=its,ite
670 do iy=jts,jte
671 z(ix,kts,iy)=0.
672 do iz=kts+1,kte+1
673 z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
674 enddo
675 iii=0
676 do iz_u=1,num_urban_hi
677 hi_urb(ix,iz_u,iy)= hi_urb2d(ix,iz_u,iy)
678 if (hi_urb(ix,iz_u,iy)/=0.) then
679 iii=iii+1
680 endif
681 enddo !iz_u
682 if (iii.gt.nbui_max) then
683 write(*,*) 'nbui_max too small, please increase to at least ',iii
684 stop
685 endif
686 enddo
687 enddo
690 if (first) then ! True only on first call
692 call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
693 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
694 emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b, &
695 cop_u,pwin_u,beta_u,sw_cond_u,time_on_u,time_off_u,targtemp_u, &
696 bldac_frc_u,cooled_frc_u, &
697 gaptemp_u,targhum_u,gaphum_u,perflo_u, &
698 gr_frac_roof_u,pv_frac_roof_u, &
699 hsesf_u,hsequip,irho,gr_flag_u,gr_type_u)
701 !Initialisation of the urban parameters and calculation of the view factor
702 call icBEP(nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u)
704 first=.false.
706 endif ! first
708 do ix=its,ite
709 do iy=jts,jte
710 if (FRC_URB2D(ix,iy).gt.0.) then ! Calling BEP only for existing urban classes.
712 iurb=UTYPE_URB2D(ix,iy)
714 hi_urb1D=0.
715 do iz_u=1,nz_um
716 hi_urb1D(iz_u)=hi_urb(ix,iz_u,iy)
718 enddo
720 call icBEPHI_XY(iurb,hb_u,hi_urb1D,ss_urb,pb_urb, &
721 nz_urb(iurb),z_u)
723 call param(iurb,nz_u(iurb),nz_urb(iurb),nzurban(iurb), &
724 nd_u(iurb),csg_u,csg,alag_u,alag,csr_u,csr, &
725 alar_u,alar,csw_u,csw,alaw_u,alaw, &
726 ws_u,ws_urb,ws,bs_u,bs_urb,bs,z0g_u,z0r_u,z0, &
727 strd_u,strd,drst_u,drst,ss_u,ss_urb,ss,pb_u, &
728 pb_urb,pb,dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb, &
729 lp_urb2d(ix,iy),lb_urb2d(ix,iy), &
730 hgt_urb2d(ix,iy),FRC_URB2D(ix,iy))
732 !
733 !We compute the view factors in the icBEP_XY routine
734 !
736 call icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u,fws_u,fsg_u, &
737 nd_u(iurb),strd,ws,nzurban(iurb),z_u)
739 ibui=0
740 nlev=0
741 nbui=0
742 d_urb=0.
743 do iz=1,nz_um
744 if(ss_urb(iz,iurb).gt.0) then
745 ibui=ibui+1
746 nlev(ibui)=iz-1
747 d_urb(ibui)=ss_urb(iz,iurb)
748 nbui=ibui
749 endif
750 end do !iz
752 if (nbui.gt.nbui_max) then
753 write (*,*) 'nbui_max must be increased to',nbui
754 stop
755 endif
759 do iz= kts,kte
760 ua1D(iz)=u_phy(ix,iz,iy)
761 va1D(iz)=v_phy(ix,iz,iy)
762 pt1D(iz)=th_phy(ix,iz,iy)
763 da1D(iz)=rho(ix,iz,iy)
764 pr1D(iz)=p_phy(ix,iz,iy)
765 pt01D(iz)=300.
766 z1D(iz)=z(ix,iz,iy)
767 qv1D(iz)=qv_phy(ix,iz,iy)
768 a_u1D(iz)=0.
769 a_v1D(iz)=0.
770 a_t1D(iz)=0.
771 a_e1D(iz)=0.
772 b_u1D(iz)=0.
773 b_v1D(iz)=0.
774 b_t1D(iz)=0.
775 b_ac1D(iz)=0.
776 b_e1D(iz)=0.
777 enddo
778 z1D(kte+1)=z(ix,kte+1,iy)
782 do id=1,ndm
783 do iz_u=1,nz_um
784 do iw=1,nwr_u
785 do ibui=1,nbui_max
786 tw1D(2*id-1,iz_u,iw,ibui)=tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
787 tw1D(2*id,iz_u,iw,ibui)=tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
788 enddo
789 enddo
790 enddo
791 enddo
793 do id=1,ndm
794 do ig=1,ng_u
795 tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
796 enddo
798 do iz_u=1,nz_um
799 do ir=1,nwr_u
800 tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)
801 enddo
802 do ir=1,ngr_u
803 if(gr_flag_u.eq.1)then
804 trv1D(id,iz_u,ir)=trv_urb4d(ix,ind_grd(iz_u,ir,id),iy)
805 qr1D(id,iz_u,ir)=qr_urb4d(ix,ind_grd(iz_u,ir,id),iy)
806 else
807 trv1D(id,iz_u,ir)=0.
808 qr1D(id,iz_u,ir)=0.
809 endif
810 enddo
811 enddo
812 enddo
816 !Initialize variables for BEM
818 tlev1D=0. !Indoor temperature
819 qlev1D=0. !Indoor humidity
821 twlev1D=0. !Window temperature
822 tglev1D=0. !Ground temperature
823 tflev1D=0. !Floor temperature
825 sflev1D=0. !Sensible heat flux from the a.c.
826 lflev1D=0. !latent heat flux from the a.c.
827 consumlev1D=0.!consumption of the a.c.
828 eppvlev1D=0. !electricity production of PV panels
829 tpvlev1D=0.
830 sfvlev1D=0. !Sensible heat flux from natural ventilation
831 lfvlev1D=0. !Latent heat flux from natural ventilation
832 sfwin1D=0. !Sensible heat flux from windows
833 sfw1D=0. !Sensible heat flux from walls
835 do iz_u=1,nz_um !vertical levels
836 do ibui=1,nbui_max !Type of building
837 tlev1D(iz_u,ibui)= tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)
838 qlev1D(iz_u,ibui)= qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)
839 enddo !ibui
840 enddo !iz_u
844 do id=1,ndm !direction
845 do iz_u=1,nz_um !vertical levels
846 do ibui=1,nbui_max !type of building
847 twlev1D(2*id-1,iz_u,ibui)=tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
848 twlev1D(2*id,iz_u,ibui)=tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
849 sfwin1D(2*id-1,iz_u,ibui)=sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
850 sfwin1D(2*id,iz_u,ibui)=sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
851 enddo !ibui
852 enddo !iz_u
853 enddo !id
855 do id=1,ndm !direction
856 do iw=1,ngb_u !layer in the wall
857 do ibui=1,nbui_max !type of building
858 tglev1D(id,iw,ibui)=tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)
859 enddo !ibui
860 enddo !iw
861 enddo !id
863 do id=1,ndm !direction
864 do iw=1,nf_u !layer in the walls
865 do iz_u=1,nz_um-1 !verticals levels
866 do ibui=1,nbui_max !type of building
867 tflev1D(id,iw,iz_u,ibui)=tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)
869 enddo !ibui
870 enddo ! iz_u
871 enddo !iw
872 enddo !id
874 !
875 !End initialization for BEM
876 !
878 do id=1,ndm
879 do iz=1,nz_um
880 do ibui=1,nbui_max !type of building
881 !! sfw1D(2*id-1,iz)=sfw1(ix,iy,ind_zd(iz,id))
882 !! sfw1D(2*id,iz)=sfw2(ix,iy,ind_zd(iz,id))
883 sfw1D(2*id-1,iz,ibui)=sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)
884 sfw1D(2*id,iz,ibui)=sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)
885 enddo
886 enddo
887 enddo
889 do id=1,ndm
890 sfg1D(id)=sfg_urb3d(ix,id,iy)
891 lfg1D(id)=lfg_urb3d(ix,id,iy)
892 dg1D(id)=dg_urb3d(ix,id,iy)
894 enddo
896 do id=1,ndm
897 do iz=1,nz_um
898 tpvlev1D(id,iz)=t_pv_urb3d(ix,ind_zdf(iz,id),iy)
899 sfr1D(id,iz)=sfr_urb3d(ix,ind_zdf(iz,id),iy)
900 lfr1D(id,iz)=lfr_urb3d(ix,ind_zdf(iz,id),iy)
901 dgr1D(id,iz)=dgr_urb3d(ix,ind_zdf(iz,id),iy)
902 if(gr_flag_u.eq.1)then
903 sfrv1D(id,iz)=sfrv_urb3d(ix,ind_zdf(iz,id),iy)
904 lfrv1D(id,iz)=lfrv_urb3d(ix,ind_zdf(iz,id),iy)
905 drain1D(id,iz)=drain_urb4d(ix,ind_zdf(iz,id),iy)
906 else
907 sfrv1D(id,iz)=0.
908 lfrv1D(id,iz)=0.
909 drain1D(id,iz)=0.
910 endif
911 enddo
912 enddo
916 rs1D=swdown(ix,iy)
917 rld1D=glw(ix,iy)
918 swddir1D=swddir(ix,iy) !_gl
919 swddif1D=swddif(ix,iy) !_gl
920 zr1D=acos(COSZ_URB2D(ix,iy))
921 deltar1D=DECLIN_URB
922 ah1D=OMG_URB2D(ix,iy)
925 call BEP1D(itimestep,ix,iy,iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
926 zr1D,deltar1D,ah1D,rs1D,rld1D,alagb, &
927 alag,alaw,alar,alaf,csgb,csg,csw,csr,csf, &
928 dzr,dzf,dzw,dzgb,xlat(ix,iy),swddir1D,swddif1D, &
929 albg_u(iurb),albw_u(iurb),albr_u(iurb), &
930 albwin_u(iurb),emg_u(iurb),emw_u(iurb), &
931 emr_u(iurb),emwind_u(iurb),fww_u,fwg_u, &
932 fgw_u,fsw_u,fws_u,fsg_u,z0, &
933 nd_u(iurb),strd,drst,ws,bs_urb,bs,ss,pb, &
934 nzurban(iurb),z_u,cop_u,pwin_u,beta_u, &
935 sw_cond_u,time_on_u,time_off_u,targtemp_u, &
936 gaptemp_u,targhum_u,gaphum_u,perflo_u, &
937 gr_frac_roof_u(iurb),pv_frac_roof_u(iurb), &
938 hsesf_u,hsequip,irho,gr_flag_u,gr_type_u, &
939 tw1D,tg1D,tr1D,trv1D,sfw1D,sfg1D,sfr1D, &
940 sfrv1D,lfrv1D, &
941 dgr1D,dg1D,lfr1D,lfg1D, &
942 drain1D,rainbl(ix,iy),qr1D, &
943 a_u1D,a_v1D,a_t1D,a_e1D, &
944 b_u1D,b_v1D,b_t1D,b_ac1D,b_e1D,b_q1D, &
945 dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy), &
946 rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy), &
947 qv1D,tlev1D,qlev1D,sflev1D,lflev1D,consumlev1D, &
948 eppvlev1D,tpvlev1D,sfvlev1D,lfvlev1D,twlev1D,tglev1D,tflev1D,sfwin1D,tair1D,sfr_indoor1D,sfrpv1D,gfr1D)
950 do ibui=1,nbui_max !type of building
951 do iz=1,nz_um !vertical levels
952 do id=1,ndm ! direction
953 sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id-1,iz,ibui)
954 sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id,iz,ibui)
955 enddo
956 enddo
957 enddo
959 do id=1,ndm
960 sfg_urb3d(ix,id,iy)=sfg1D(id)
961 lfg_urb3d(ix,id,iy)=lfg1D(id)
962 dg_urb3d(ix,id,iy)=dg1D(id)
963 enddo
965 do id=1,ndm
966 do iz=1,nz_um
967 t_pv_urb3d(ix,ind_zdf(iz,id),iy)=tpvlev1D(id,iz)
968 sfr_urb3d(ix,ind_zdf(iz,id),iy)=sfr1D(id,iz)
969 dgr_urb3d(ix,ind_zdf(iz,id),iy)=dgr1D(id,iz)
970 lfr_urb3d(ix,ind_zdf(iz,id),iy)=lfr1D(id,iz)
971 if(gr_flag_u.eq.1)then
972 sfrv_urb3d(ix,ind_zdf(iz,id),iy)=sfrv1D(id,iz)
973 lfrv_urb3d(ix,ind_zdf(iz,id),iy)=lfrv1D(id,iz)
974 drain_urb4d(ix,ind_zdf(iz,id),iy)=drain1D(id,iz)
975 endif
976 enddo
977 enddo
979 do ibui=1,nbui_max
980 do iz_u=1,nz_um
981 do iw=1,nwr_u
982 do id=1,ndm
983 tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw,ibui)
984 tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw,ibui)
985 enddo
986 enddo
987 enddo
988 enddo
991 do id=1,ndm
992 do ig=1,ng_u
994 tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
995 enddo
996 do iz_u=1,nz_um
997 do ir=1,nwr_u
998 trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
999 enddo
1000 if(gr_flag_u.eq.1)then
1001 do ir=1,ngr_u
1002 trv_urb4d(ix,ind_grd(iz_u,ir,id),iy)=trv1D(id,iz_u,ir)
1003 qr_urb4d(ix,ind_grd(iz_u,ir,id),iy)=qr1D(id,iz_u,ir)
1004 enddo
1005 endif
1006 enddo
1007 enddo
1008 !
1010 !
1011 !Outputs of BEM
1012 !
1014 do ibui=1,nbui_max !type of building
1015 do iz_u=1,nz_um !vertical levels
1016 tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=tlev1D(iz_u,ibui)
1017 qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=qlev1D(iz_u,ibui)
1018 enddo !iz_u
1019 enddo !ibui
1021 do ibui=1,nbui_max !type of building
1022 do iz_u=1,nz_um !vertical levels
1023 do id=1,ndm !direction
1024 tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id-1,iz_u,ibui)
1025 tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id,iz_u,ibui)
1026 sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id-1,iz_u,ibui)
1027 sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id,iz_u,ibui)
1028 enddo !id
1029 enddo !iz_u
1030 enddo !ibui
1032 do ibui=1,nbui_max !type of building
1033 do iw=1,ngb_u !layers in the walls
1034 do id=1,ndm !direction
1035 tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)=tglev1D(id,iw,ibui)
1036 enddo !id
1037 enddo !iw
1038 enddo !ibui
1040 do ibui=1,nbui_max !type of building
1041 do iw=1,nf_u !layer in the walls
1042 do iz_u=1,nz_um-1 !verticals levels
1043 do id=1,ndm
1044 tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)=tflev1D(id,iw,iz_u,ibui)
1045 enddo !ibui
1046 enddo ! iz_u
1047 enddo !iw
1048 enddo !id
1052 sf_ac_urb3d(ix,iy)=0.
1053 lf_ac_urb3d(ix,iy)=0.
1054 cm_ac_urb3d(ix,iy)=0.
1055 ep_pv_urb3d(ix,iy)=0.
1056 sfvent_urb3d(ix,iy)=0.
1057 lfvent_urb3d(ix,iy)=0.
1058 draingr_urb3d(ix,iy)=0.
1059 qgr_urb3d(ix,iy)=0.
1060 tgr_urb3d(ix,iy)=0.
1061 meso_urb=(1./4.)*FRC_URB2D(ix,iy)/((bs_urb(1,iurb)+ws_urb(1,iurb))*bs_urb(2,iurb))+ &
1062 (1./4.)*FRC_URB2D(ix,iy)/((bs_urb(2,iurb)+ws_urb(2,iurb))*bs_urb(1,iurb))
1063 meso_urb_ac=meso_urb*bldac_frc_u(iurb)*cooled_frc_u(iurb)
1064 roof_frac=FRC_URB2D(ix,iy)*bs_urb(1,iurb)/(bs_urb(1,iurb)+ws_urb(1,iurb))
1065 ibui=0
1066 nlev=0
1067 nbui=0
1068 d_urb=0.
1069 do iz=1,nz_um
1070 if(ss_urb(iz,iurb).gt.0) then
1071 ibui=ibui+1
1072 nlev(ibui)=iz-1
1073 d_urb(ibui)=ss_urb(iz,iurb)
1074 nbui=ibui
1075 endif
1076 end do !iz
1081 do ibui=1,nbui !type of building
1082 ep_pv_urb3d(ix,iy)=ep_pv_urb3d(ix,iy)+meso_urb_ac*d_urb(ibui)*eppvlev1D(ibui)
1083 do iz_u=1,nlev(ibui) !vertical levels
1084 sf_ac_urb3d(ix,iy)=sf_ac_urb3d(ix,iy)+meso_urb_ac*d_urb(ibui)*sflev1D(iz_u,ibui)
1085 lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+meso_urb_ac*d_urb(ibui)*lflev1D(iz_u,ibui)
1086 cm_ac_urb3d(ix,iy)=cm_ac_urb3d(ix,iy)+meso_urb_ac*d_urb(ibui)*consumlev1D(iz_u,ibui)
1087 !if(consumlev1D(iz_u,ibui).gt.0.)then
1088 !print*,'IX',ix,'IY',iy,'IZ_U',iz_u,'IBUI',ibui,'CONSUM',consumlev1D(iz_u,ibui),'D_URB',d_urb(ibui),'MESO_URB',meso_urb_ac
1090 !endif
1091 sfvent_urb3d(ix,iy)=sfvent_urb3d(ix,iy)+meso_urb_ac*d_urb(ibui)*sfvlev1D(iz_u,ibui)
1092 lfvent_urb3d(ix,iy)=lfvent_urb3d(ix,iy)+meso_urb_ac*d_urb(ibui)*lfvlev1D(iz_u,ibui)
1093 enddo !iz_u
1094 enddo !ibui
1098 if(gr_flag_u.eq.1)then
1099 do id=1,ndm
1100 do iz=2,nz_um-1
1101 draingr_urb3d(ix,iy)=draingr_urb3d(ix,iy)+d_urb(iz-1)*roof_frac*drain1D(id,iz)*1000
1102 do ig=1,ngr_u
1103 qgr_urb3d(ix,iy)=qgr_urb3d(ix,iy)+qr1D(id,iz,ig)/ndm/(nz_um-2)/ngr_u
1104 tgr_urb3d(ix,iy)=tgr_urb3d(ix,iy)+trv1D(id,iz,ig)/ndm/(nz_um-2)/ngr_u
1106 enddo
1107 enddo
1108 enddo
1109 endif
1110 !End outputs of bem
1111 !
1113 sf_ac=0.
1114 sf(ix,kts:kte,iy)=0.
1115 vl(ix,kts:kte,iy)=0.
1116 a_u(ix,kts:kte,iy)=0.
1117 a_v(ix,kts:kte,iy)=0.
1118 a_t(ix,kts:kte,iy)=0.
1119 a_e(ix,kts:kte,iy)=0.
1120 b_u(ix,kts:kte,iy)=0.
1121 b_v(ix,kts:kte,iy)=0.
1122 b_t(ix,kts:kte,iy)=0.
1123 b_e(ix,kts:kte,iy)=0.
1124 b_q(ix,kts:kte,iy)=0.
1125 dlg(ix,kts:kte,iy)=0.
1126 dl_u(ix,kts:kte,iy)=0.
1128 do iz= kts,kte
1129 sf(ix,iz,iy)=sf1D(iz)
1130 vl(ix,iz,iy)=vl1D(iz)
1131 a_u(ix,iz,iy)=a_u1D(iz)
1132 a_v(ix,iz,iy)=a_v1D(iz)
1133 a_t(ix,iz,iy)=a_t1D(iz)
1134 a_e(ix,iz,iy)=a_e1D(iz)
1135 b_u(ix,iz,iy)=b_u1D(iz)
1136 b_v(ix,iz,iy)=b_v1D(iz)
1137 b_t(ix,iz,iy)=b_t1D(iz)
1138 sf_ac=sf_ac+b_ac1D(iz)*da1D(iz)*cp_u*dz8w(ix,iz,iy)*vl1D(iz)*FRC_URB2D(ix,iy)
1139 b_e(ix,iz,iy)=b_e1D(iz)
1140 b_q(ix,iz,iy)=b_q1D(iz)
1141 dlg(ix,iz,iy)=dlg1D(iz)
1142 dl_u(ix,iz,iy)=dl_u1D(iz)
1143 enddo
1144 sf(ix,kte+1,iy)=sf1D(kte+1)
1146 endif ! FRC_URB2D
1149 enddo ! iy
1150 enddo ! ix
1153 time_bep=time_bep+dt
1155 ! print*, 'ss_urb', ss_urb
1156 ! print*, 'pb_urb', pb_urb
1157 ! print*, 'nz_urb', nz_urb
1158 ! print*, 'd_urb', d_urb
1161 return
1162 end subroutine BEP_BEM
1165 ! ===6=8===============================================================72
1167 subroutine BEP1D(itimestep,ix,iy,iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, &
1168 zr,deltar,ah,rs,rld,alagb, &
1169 alag,alaw,alar,alaf,csgb,csg,csw,csr,csf, &
1170 dzr,dzf,dzw,dzgb,xlat,swddir,swddif, &
1171 albg,albw,albr,albwin,emg,emw,emr, &
1172 emwind,fww,fwg,fgw,fsw,fws,fsg,z0, &
1173 ndu,strd,drst,ws,bs_u,bs,ss,pb, &
1174 nzu,z_u,cop_u,pwin_u,beta_u,sw_cond_u, &
1175 time_on_u,time_off_u,targtemp_u, &
1176 gaptemp_u,targhum_u,gaphum_u,perflo_u, &
1177 gr_frac_roof,pv_frac_roof, &
1178 hsesf_u,hsequip,irho,gr_flag,gr_type, &
1179 tw,tg,tr,trv,sfw,sfg,sfr, &
1180 sfrv,lfrv,dgr,dg,lfr,lfg,drain,rainbl,qr, &
1181 a_u,a_v,a_t,a_e, &
1182 b_u,b_v,b_t,b_ac,b_e,b_q, &
1183 dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb, &
1184 qv,tlev,qlev,sflev,lflev,consumlev, &
1185 eppvlev,tpvlev,sfvlev,lfvlev,twlev,tglev,tflev,sfwin,tmp_u,sfr_indoor,sfrpv,gfr)
1186 ! print*,'SFR_AFT',sfr(id,iz)
1190 implicit none
1192 ! ----------------------------------------------------------------------
1193 ! INPUT:
1194 ! ----------------------------------------------------------------------
1196 ! Data relative to the "mesoscale grid"
1198 !! integer nz ! Number of vertical levels
1199 integer kms,kme,kts,kte,ix,iy,itimestep
1200 real z(kms:kme) ! Altitude above the ground of the cell interfaces.
1201 real ua(kms:kme) ! Wind speed in the x direction
1202 real va(kms:kme) ! Wind speed in the y direction
1203 real pt(kms:kme) ! Potential temperature
1204 real da(kms:kme) ! Air density
1205 real pr(kms:kme) ! Air pressure
1206 real pt0(kms:kme) ! Reference potential temperature (could be equal to "pt")
1207 real qv(kms:kme) ! Specific humidity
1208 real dt ! Time step
1209 real zr ! Zenith angle
1210 real deltar ! Declination of the sun
1211 real ah ! Hour angle
1212 real rs ! Solar radiation
1213 real rld ! Downward flux of the longwave radiation
1214 real xlat ! Latitude
1215 real swddir ! short wave direct solar radiation !_gl
1216 real swddif ! short wave diffuse solar radiation !_gl
1218 ! Data relative to the "urban grid"
1220 integer iurb ! Current urban class
1222 ! Radiation parameters
1223 real albg ! Albedo of the ground
1224 real albw ! Albedo of the wall
1225 real albr ! Albedo of the roof
1226 real albwin ! Albedo of the windows
1227 real emwind ! Emissivity of windows
1228 real emg ! Emissivity of ground
1229 real emw ! Emissivity of wall
1230 real emr ! Emissivity of roof
1233 ! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and
1234 ! short wave radation.
1235 ! The calculation of these factor is explained in the Appendix A of the BLM paper
1236 real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
1237 real fwg(nz_um,ndm,nurbm) ! from wall to ground
1238 real fgw(nz_um,ndm,nurbm) ! from ground to wall
1239 real fsw(nz_um,ndm,nurbm) ! from sky to wall
1240 real fws(nz_um,ndm,nurbm) ! from wall to sky
1241 real fsg(ndm,nurbm) ! from sky to ground
1243 ! Street parameters
1244 integer ndu ! Number of street direction for each urban class
1245 real bs_u(ndm,nurbm) ! Building width
1247 ! Grid parameters
1248 integer nzu ! Number of layer in the urban grid
1249 real z_u(nz_um) ! Height of the urban grid levels
1250 !FS
1251 real cop_u(nurbm)
1252 real pwin_u(nurbm)
1253 real beta_u(nurbm)
1254 integer sw_cond_u(nurbm)
1255 real time_on_u(nurbm)
1256 real time_off_u(nurbm)
1257 real targtemp_u(nurbm)
1258 real gaptemp_u(nurbm)
1259 real targhum_u(nurbm)
1260 real gaphum_u(nurbm)
1261 real perflo_u(nurbm)
1262 real hsesf_u(nurbm)
1263 real hsequip(24)
1264 real irho(24)
1265 real gr_frac_roof
1266 real pv_frac_roof
1267 integer gr_flag
1268 integer gr_type
1269 real tpv(nbui_max)
1270 real sfpv(nbui_max)
1271 real sfr_indoor(nbui_max)
1272 ! ----------------------------------------------------------------------
1273 ! INPUT-OUTPUT
1274 ! ----------------------------------------------------------------------
1276 ! Data relative to the "urban grid" which should be stored from the current time step to the next one
1278 real tw(2*ndm,nz_um,nwr_u,nbui_max) ! Temperature in each layer of the wall [K]
1279 real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
1280 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1281 real trv(ndm,nz_um,ngr_u) ! Temperature in each layer of the green roof [K]
1282 real sfw(2*ndm,nz_um,nbui_max) ! Sensible heat flux from walls
1283 real sfg(ndm) ! Sensible heat flux from ground (road)
1284 real sfr(ndm,nz_um) ! Sensible heat flux from roofs
1285 real sfrv(ndm,nz_um) ! Sensible heat flux from green roofs
1286 real lfrv(ndm,nz_um) ! Latent heat flux from green roofs
1287 real dg(ndm) ! water depth ground (road)
1288 real dgr(ndm,nz_um) ! water depth roofs
1289 real lfr(ndm,nz_um) ! Latent heat flux from roofs
1290 real lfg(ndm) ! Latent heat flux from ground (road)
1291 real drain(ndm,nz_um) ! Green roof drainage
1292 real rainbl ! Rainfall
1293 real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) towards the interior
1294 real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof towards the interior
1295 real gfw(2*ndm,nz_um,nbui_max) ! Heat flux transfered from the surface of the walls towards the interior
1296 real qr(ndm,nz_um,ngr_u) ! Green Roof soil moisture
1298 ! ----------------------------------------------------------------------
1299 ! OUTPUT:
1300 ! ----------------------------------------------------------------------
1303 ! Data relative to the "mesoscale grid"
1305 real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
1306 real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
1308 ! Implicit and explicit components of the source and sink terms at each levels,
1309 ! the fluxes can be computed as follow: FX = A*X + B example: Heat fluxes = a_t * pt + b_t
1310 real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
1311 real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
1312 real a_t(kms:kme) ! Implicit component of the heat sources or sinks
1313 real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
1314 real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
1315 real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
1316 real b_t(kms:kme) ! Explicit component of the heat sources or sinks
1317 real b_ac(kms:kme)
1318 real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
1319 real b_q(kms:kme) ! Explicit component of the humidity sources or sinks
1320 real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
1321 real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
1322 ! ----------------------------------------------------------------------
1323 ! LOCAL:
1324 ! ----------------------------------------------------------------------
1326 real dz(kms:kme) ! vertical space steps of the "mesoscale grid"
1327 ! Data interpolated from the "mesoscale grid" to the "urban grid"
1329 real ua_u(nz_um) ! Wind speed in the x direction
1330 real va_u(nz_um) ! Wind speed in the y direction
1331 real pt_u(nz_um) ! Potential temperature
1332 real da_u(nz_um) ! Air density
1333 real pt0_u(nz_um) ! Reference potential temperature
1334 real pr_u(nz_um) ! Air pressure
1335 real qv_u(nz_um) !Specific humidity
1337 ! Data defining the building and street charateristics
1339 real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
1341 real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
1342 real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
1343 real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
1345 real z0(ndm,nz_um) ! Roughness lengths "profiles"
1346 real ws(ndm) ! Street widths of the current urban class
1347 real bs(ndm) ! Building widths of the current urban class
1348 real strd(ndm) ! Street lengths for the current urban class
1349 real drst(ndm) ! Street directions for the current urban class
1350 real ss(nz_um) ! Probability to have a building with height h
1351 real pb(nz_um) ! Probability to have a building with an height equal
1352 real cdrag(nz_um)
1353 real alp
1355 ! Solar radiation at each level of the "urban grid"
1357 real rsg(ndm) ! Short wave radiation from the ground
1358 real rsw(2*ndm,nz_um) ! Short wave radiation from the walls
1359 real rsd(2*ndm,nz_um) ! Direct Short wave radiation received by the walls
1360 real rlg(ndm) ! Long wave radiation from the ground
1361 real rlw(2*ndm,nz_um) ! Long wave radiation from the walls
1363 ! Potential temperature of the surfaces at each level of the "urban grid"
1365 real ptg(ndm) ! Ground potential temperatures
1366 real ptr(ndm,nz_um) ! Roof potential temperatures
1367 real ptrv(ndm,nz_um) ! Roof potential temperatures
1368 real ptw(2*ndm,nz_um,nbui_max) ! Walls potential temperatures
1370 real tg_av(ndm)
1371 ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
1372 ! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
1373 ! The fluxes can be computed as follow: Fluxes of X = A*X + B
1374 ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
1376 real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
1377 real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
1378 real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
1379 real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
1380 real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
1381 real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
1382 real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1383 real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1384 real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1387 real tvb_ac(2*ndm,nz_um)
1388 real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1389 real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1390 real qhb_u(ndm,nz_um) ! Humidity Horizontal surfaces, B (explicit) term
1391 real qvb_u(2*ndm,nz_um) ! Humidity Vertical surfaces, B (explicit) term
1392 !
1393 real rs_abs ! solar radiation absorbed by urban surfaces
1394 real rl_up ! longwave radiation emitted by urban surface to the atmosphere
1395 real emiss ! mean emissivity of the urban surface
1396 real grdflx_urb ! ground heat flux
1397 real dt_int ! internal time step
1398 integer nt_int ! number of internal time step
1399 integer iz,id, it_int,it
1400 integer iw
1402 !---------------------------------------
1403 !New variables uses in BEM
1404 !----------------------------------------
1406 real tmp_u(nz_um) !Air Temperature [K]
1408 real dzw(nwr_u) !Layer sizes in the walls
1409 real dzr(nwr_u) !Layer sizes in the roofs
1410 real dzf(nf_u) !Layer sizes in the floors
1411 real dzgb(ngb_u) !Layer sizes in the ground below the buildings
1413 real csgb(ngb_u) !Specific heat of the ground material below the buildings
1415 real csf(nf_u) !Specific heat of the floors materials in the buildings
1416 !of the current urban class at each levels[J m^3 K^-1]
1417 real alar(nwr_u+1) ! Roof thermal diffusivity for the current urban class [W/m K]
1418 real alaw(nwr_u+1) ! Walls thermal diffusivity for the current urban class [W/m K]
1419 real alaf(nf_u+1) ! Floor thermal diffusivity at each wall layers [W/m K]
1420 real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall layer [W/m K]
1422 real sfrb(ndm,nbui_max) ! Sensible heat flux from roofs [W/m2]
1423 real sfrbpv(ndm,nbui_max) ! Sensible heat flux from PV panels [W/m2]
1424 real sfrpv(ndm,nz_um) ! Sensible heat flux from PV panels [W/m2]
1425 real sfrvb(ndm,nbui_max) ! Sensible heat flux from roofs [W/m2]
1426 real lfrvb(ndm,nbui_max) ! Sensible heat flux from roofs [W/m2]
1427 real lfrb(ndm,nbui_max) ! Sensible heat flux from roofs [W/m2]
1429 real gfrb(ndm,nbui_max) ! Heat flux flowing inside the roofs [W/m2]
1430 real sfwb1D(2*ndm,nz_um) !Sensible heat flux from the walls [W/m2]
1431 real sfwin(2*ndm,nz_um,nbui_max)!Sensible heat flux from windows [W/m2]
1432 real sfwinb1D(2*ndm,nz_um) !Sensible heat flux from windows [W/m2]
1433 real gfwb1D(2*ndm,nz_um) !Heat flux flowing inside the walls [W/m2]
1435 real qlev(nz_um,nbui_max) !specific humidity [kg/kg]
1436 real qlevb1D(nz_um) !specific humidity [kg/kg]
1437 real tlev(nz_um,nbui_max) !Indoor temperature [K]
1438 real tlevb1D(nz_um) !Indoor temperature [K]
1439 real twb1D(2*ndm,nwr_u,nz_um) !Wall temperature in BEM [K]
1440 real twlev(2*ndm,nz_um,nbui_max) !Window temperature in BEM [K]
1441 real twlevb1D(2*ndm,nz_um) !Window temperature in BEM [K]
1442 real tglev(ndm,ngb_u,nbui_max) !Ground temperature below a building in BEM [K]
1443 real tglevb1D(ngb_u) !Ground temperature below a building in BEM [K]
1444 real tflev(ndm,nf_u,nz_um-1,nbui_max)!Floor temperature in BEM[K]
1445 real tflevb1D(nf_u,nz_um-1) !Floor temperature in BEM[K]
1446 real trb(ndm,nwr_u,nbui_max) !Roof temperature in BEM [K]
1447 real trvb(ndm,ngr_u,nbui_max) !Roof temperature in BEM [K]
1448 real trb1D(nwr_u)
1450 real sflev(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems [W]
1451 real lflev(nz_um,nz_um) ! latent heat flux due to the air conditioning systems [W]
1452 real consumlev(nz_um,nz_um) ! consumption due to the air conditioning systems [W]
1453 real sflev1D(nz_um) ! sensible heat flux due to the air conditioning systems [W]
1454 real lflev1D(nz_um) ! latent heat flux due to the air conditioning systems [W]
1455 real consumlev1D(nz_um) ! consumption due to the air conditioning systems [W]
1456 real eppvlev(nz_um) ! Electricity production of PV panels [W]
1457 real tpvlev(ndm,nz_um)
1458 real tpvlevb(ndm,nbui_max) ! Sensible heat flux from roofs [W/m2]
1459 real sfvlev(nz_um,nz_um) ! sensible heat flux due to ventilation [W]
1460 real lfvlev(nz_um,nz_um) ! latent heat flux due to ventilation [W]
1461 real sfvlev1D(nz_um) ! sensible heat flux due to ventilation [W]
1462 real lfvlev1D(nz_um) ! Latent heat flux due to ventilation [W]
1464 real ptwin(2*ndm,nz_um,nbui_max) ! window potential temperature
1465 real tw_av(2*ndm,nz_um) ! Averaged temperature of the wall surfaces
1466 real twlev_av(2*ndm,nz_um) ! Averaged temperature of the windows
1467 real sfw_av(2*ndm,nz_um) ! Averaged sensible heat from walls
1468 real sfwind_av(2*ndm,nz_um) ! Averaged sensible heat from windows
1469 integer flag_pvp
1470 integer nbui !Total number of different type of buildings in an urban class
1471 integer nlev(nz_um) !Number of levels in each different type of buildings in an urban class
1472 integer ibui,ily
1473 real :: nhourday ! Number of hours from midnight, local time
1474 real :: st4,gamma,fp,lmr,smr,prova
1475 real hfgr(ndm,nz_um)!heat flux green roof
1476 real hfgrb(ndm,nbui_max)
1477 real irri_per_ts
1478 real irri_now
1479 real tr_av(ndm,nz_um)
1480 real tr_avb(ndm,nbui_max)
1481 real sfr_avb(ndm,nbui_max)
1482 ! ----------------------------------------------------------------------
1483 ! END VARIABLES DEFINITIONS
1484 ! ----------------------------------------------------------------------
1486 ! Fix some usefull parameters for the computation of the sources or sinks
1487 !
1488 !initialize the variables inside the param routine
1490 nhourday=ah/PI*180./15.+12.
1491 if (nhourday >= 24) nhourday = nhourday - 24
1492 if (nhourday < 0) nhourday = nhourday + 24
1495 if(sum(irho).gt.0)then
1496 irri_per_ts=h_water/sum(irho)
1497 else
1498 irri_per_ts=0.
1499 endif
1501 if(irho(int(nhourday)+1).ne.0)then
1502 irri_now=irri_per_ts
1503 else
1504 irri_now=0.
1505 endif
1507 do iz=kts,kte
1508 dz(iz)=z(iz+1)-z(iz)
1509 end do
1510 ! Interpolation on the "urban grid"
1511 call interpol(kms,kme,kts,kte,nzu,z,z_u,ua,ua_u)
1512 call interpol(kms,kme,kts,kte,nzu,z,z_u,va,va_u)
1513 call interpol(kms,kme,kts,kte,nzu,z,z_u,pt,pt_u)
1514 call interpol(kms,kme,kts,kte,nzu,z,z_u,pt0,pt0_u)
1515 call interpol(kms,kme,kts,kte,nzu,z,z_u,pr,pr_u)
1516 call interpol(kms,kme,kts,kte,nzu,z,z_u,da,da_u)
1517 call interpol(kms,kme,kts,kte,nzu,z,z_u,qv,qv_u)
1518 ! Compute the modification of the radiation due to the buildings
1521 call averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av, &
1522 sfw_av,sfwind_av,sfw,sfwin)
1524 do id=1,ndu
1525 tg_av(id)=tg(id,ng_u)
1526 do iz=1,nz_um
1528 tr_av(id,iz)=((1-gr_frac_roof)*tr(id,iz,nwr_u)**4.+ &
1529 gr_frac_roof*trv(id,iz,ngr_u)**4.)**(1./4.)
1531 enddo
1532 enddo
1538 call modif_rad(iurb,ndu,nzu,z_u,ws, &
1539 drst,strd,ss,pb, &
1540 tw_av,tg_av,twlev_av,albg,albw, &
1541 emw,emg,pwin_u(iurb),albwin, &
1542 emwind,fww,fwg,fgw,fsw,fsg, &
1543 zr,deltar,ah,xlat,swddir,swddif, & !_gl
1544 rs,rld,rsw,rsd,rsg,rlw,rlg)
1549 ! calculation of the urban albedo and the upward long wave radiation
1552 call upward_rad(ndu,nzu,ws,bs,sigma,pb,ss, &
1553 tg_av,emg,albg,rlg,rsg,sfg,lfg, &
1554 tw_av,emw,albw,rlw,rsw,sfw_av, &
1555 tr_av,emr,albr,emwind, &
1556 albwin,twlev_av,pwin_u(iurb),sfwind_av,rld,rs,sfr,sfrv,lfr,lfrv, &
1557 rs_abs,rl_up,emiss,grdflx_urb,gr_frac_roof,tpvlev,pv_frac_roof)
1559 do id=1,ndu
1560 if(dg(id).le.dgmax) then
1561 dg(id)=dg(id)+(rainbl+(lfg(id)*dt)/latent)
1562 endif
1563 if (dg(id).lt.0) then
1564 dg(id)=0
1565 endif
1566 if (dg(id).gt.dgmax) then
1567 dg(id)=dgmax
1568 endif
1569 do iz=2,nz_um
1570 if(dgr(id,iz).le.drmax) then
1571 dgr(id,iz)=dgr(id,iz)+(rainbl+(lfr(id,iz)*dt)/latent)
1572 endif
1573 if (dgr(id,iz).lt.0) then
1574 dgr(id,iz)=0
1575 endif
1576 if (dgr(id,iz).gt.drmax) then
1577 dgr(id,iz)=drmax
1578 endif
1579 enddo
1580 enddo !id
1585 call surf_temp(ndu,pr_u,dt, &
1586 rld,rsg,rlg, &
1587 tg,alag,csg,emg,albg,ptg,sfg,lfg,gfg)
1588 if(gr_flag.eq.1)then
1589 if(gr_frac_roof.gt.0.)then
1590 hfgr=0.
1591 call roof_temp_veg(ndu,pr_u,dt, &
1592 rld,rs, &
1593 trv,ptrv,sfrv,lfrv,gfr,qr,rainbl,drain,hfgr,tr,alar(5),dzr(5),csr(5),nzu,irri_now,gr_type,pv_frac_roof,tpvlev)
1595 endif
1596 endif
1600 do iz=1,nz_um !Compute the outdoor temperature
1601 tmp_u(iz)=pt_u(iz)*(pr_u(iz)/p0)**(rcp_u)
1602 end do
1604 ibui=0
1605 nlev=0
1606 nbui=0
1607 hfgrb=0.
1608 sfrb=0. !Sensible heat flux from roof
1609 sfrbpv=0. !Sensible heat flux from PV panels
1610 sfrpv=0. !Sensible heat flux from PV panels
1611 lfrvb=0.
1612 lfrb=0.
1613 sfrvb=0.
1614 gfrb=0. !Heat flux flowing inside the roof
1615 sfwb1D=0. !Sensible heat flux from walls
1616 sfwinb1D=0. !Sensible heat flux from windows
1617 gfwb1D=0. !Heat flux flowing inside the walls[W/m2]
1620 twb1D=0. !Wall temperature
1621 twlevb1D=0. !Window temperature
1622 tglevb1D=0. !Ground temperature below a building
1623 tflevb1D=0. !Floor temperature
1624 trvb=0.
1625 trb=0. !Roof temperature
1626 trb1D=0. !Roof temperature
1627 tr_avb=0.
1628 qlevb1D=0. !Indoor humidity
1629 tlevb1D=0. !indoor temperature
1631 sflev1D=0. !Sensible heat flux from the a.c.
1632 lflev1D=0. !Latent heat flux from the a.c.
1633 consumlev1D=0.!Consumption from the a.c.
1634 tpvlevb=0.
1635 eppvlev=0.
1636 sfvlev1D=0. !Sensible heat flux from the natural ventilation
1637 lfvlev1D=0. !Latent heat flux from natural ventilation
1638 ptw=0. !Wall potential temperature
1639 ptwin=0. !Window potential temperature
1640 ptr=0. !Roof potential temperature
1642 do iz=1,nz_um
1643 if(ss(iz).gt.0) then
1644 ibui=ibui+1
1645 nlev(ibui)=iz-1
1646 nbui=ibui
1647 do id=1,ndm
1648 tr_avb(id,ibui)=tr_av(id,iz)
1649 tpvlevb(id,ibui)=tpvlev(id,iz)
1650 hfgrb(id,ibui)=hfgr(id,iz)
1651 sfrb(id,ibui)=sfr(id,iz)
1652 sfrvb(id,ibui)=sfrv(id,iz)
1653 lfrvb(id,ibui)=lfrv(id,iz)
1654 lfrb(id,ibui)=lfr(id,iz)
1655 sfr_avb(id,ibui)=(1-gr_frac_roof)*sfr(id,iz)+gr_frac_roof*(sfrv(id,iz))
1656 do ily=1,nwr_u
1657 trb(id,ily,ibui)=tr(id,iz,ily)
1658 enddo
1659 do ily=1,ngr_u
1660 trvb(id,ily,ibui)=trv(id,iz,ily)
1661 enddo
1663 enddo
1664 endif
1665 end do !iz
1667 !--------------------------------------------------------------------------------
1668 !Loop over BEM -----------------------------------------------------------------
1669 !--------------------------------------------------------------------------------
1670 !--------------------------------------------------------------------------------
1672 do ibui=1,nbui
1673 do iz=1,nz_um
1674 qlevb1D(iz)=qlev(iz,ibui)
1675 tlevb1D(iz)=tlev(iz,ibui)
1676 enddo
1678 do id=1,ndm
1680 do ily=1,nwr_u
1681 trb1D(ily)=trb(id,ily,ibui)
1682 enddo
1683 do ily=1,ngb_u
1684 tglevb1D(ily)=tglev(id,ily,ibui)
1685 enddo
1687 do ily=1,nf_u
1688 do iz=1,nz_um-1
1689 tflevb1D(ily,iz)=tflev(id,ily,iz,ibui)
1690 enddo
1691 enddo
1693 do iz=1,nz_um
1694 sfwinb1D(2*id-1,iz)=sfwin(2*id-1,iz,ibui)
1695 sfwinb1D(2*id,iz)=sfwin(2*id,iz,ibui)
1696 enddo
1698 do iz=1,nz_um
1699 do ily=1,nwr_u
1700 twb1D(2*id-1,ily,iz)=tw(2*id-1,iz,ily,ibui)
1701 twb1D(2*id,ily,iz)=tw(2*id,iz,ily,ibui)
1702 enddo
1703 sfwb1D(2*id-1,iz)=sfw(2*id-1,iz,ibui)
1704 sfwb1D(2*id,iz)=sfw(2*id,iz,ibui)
1705 twlevb1D(2*id-1,iz)=twlev(2*id-1,iz,ibui)
1706 twlevb1D(2*id,iz)=twlev(2*id,iz,ibui)
1707 enddo
1708 enddo
1710 !print*,'HFGR_BEFORE_CALLING_BEM',hfgr(nlev(ibui))
1712 call BEM(nz_um,nlev(ibui),nhourday,dt,bs_u(1,iurb), &
1713 bs_u(2,iurb),dz_u,nwr_u,nf_u,nwr_u,ngb_u,sfwb1D,gfwb1D, &
1714 sfwinb1D,sfr_avb(1,ibui),lfrb(1,ibui),gfrb(1,ibui), &
1715 sfrbpv(1,ibui), &
1716 latent,sigma,albw,albwin,albr, &
1717 emr,emw,emwind,rsw,rlw,r,cp_u, &
1718 da_u,tmp_u,qv_u,pr_u,rs,swddif,rld,dzw,csw,alaw,pwin_u(iurb), &
1719 cop_u(iurb),beta_u(iurb),sw_cond_u(iurb),time_on_u(iurb), &
1720 time_off_u(iurb),targtemp_u(iurb),gaptemp_u(iurb), &
1721 targhum_u(iurb),gaphum_u(iurb),perflo_u(iurb), &
1722 gr_frac_roof,pv_frac_roof,gr_flag, &
1723 ua_u,va_u, &
1724 hsesf_u(iurb),hsequip, &
1725 dzf,csf,alaf,dzgb,csgb,alagb,dzr,csr, &
1726 alar,tlevb1D,qlevb1D,twb1D,twlevb1D,tflevb1D,tglevb1D, &
1727 trb1D,sflev1D,lflev1D,consumlev1D,eppvlev(ibui), &
1728 tpvlevb(1,ibui), &
1729 sfvlev1D,lfvlev1D,hfgrb(1,ibui),tr_avb(1,ibui), &
1730 tpv(ibui),sfpv(ibui),sfr_indoor(ibui))
1733 !
1734 !Temporal modifications
1735 !
1736 tpvlevb(2,ibui)=tpvlevb(1,ibui)
1737 sfrb(2,ibui)=sfrb(1,ibui)
1738 sfrvb(2,ibui)=sfrvb(1,ibui)
1739 lfrvb(2,ibui)=lfrvb(1,ibui)
1740 lfrb(2,ibui)=lfrb(1,ibui)
1741 sfrbpv(2,ibui)=sfrbpv(1,ibui)
1742 gfrb(2,ibui)=gfrb(1,ibui)
1743 hfgrb(2,ibui)=hfgrb(1,ibui)
1744 !End temporal modifications
1745 !
1746 do iz=1,nz_um
1747 qlev(iz,ibui)=qlevb1D(iz)
1748 tlev(iz,ibui)=tlevb1D(iz)
1749 sflev(iz,ibui)=sflev1D(iz)
1750 lflev(iz,ibui)=lflev1D(iz)
1751 consumlev(iz,ibui)=consumlev1D(iz)
1752 sfvlev(iz,ibui)=sfvlev1D(iz)
1753 lfvlev(iz,ibui)=lfvlev1D(iz)
1754 enddo
1756 do id=1,ndm
1757 do ily=1,nwr_u
1758 trb(id,ily,ibui)=trb1D(ily)
1759 enddo
1760 do ily=1,ngb_u
1761 tglev(id,ily,ibui)=tglevb1D(ily)
1762 enddo
1764 do ily=1,nf_u
1765 do iz=1,nz_um-1
1766 tflev(id,ily,iz,ibui)=tflevb1D(ily,iz)
1767 enddo
1768 enddo
1771 do iz=1,nz_um
1772 do ily=1,nwr_u
1773 tw(2*id-1,iz,ily,ibui)=twb1D(2*id-1,ily,iz)
1774 tw(2*id,iz,ily,ibui)=twb1D(2*id,ily,iz)
1775 enddo
1776 gfw(2*id-1,iz,ibui)=gfwb1D(2*id-1,iz)
1777 gfw(2*id,iz,ibui)=gfwb1D(2*id,iz)
1778 twlev(2*id-1,iz,ibui)=twlevb1D(2*id-1,iz)
1779 twlev(2*id,iz,ibui)=twlevb1D(2*id,iz)
1780 enddo
1781 enddo
1783 enddo !ibui
1785 !-----------------------------------------------------------------------------
1786 !End loop over BEM -----------------------------------------------------------
1787 !-----------------------------------------------------------------------------
1788 !-----------------------------------------------------------------------------
1790 ibui=0
1792 do iz=1,nzu!nz_um
1794 if(ss(iz).gt.0) then
1795 ibui=ibui+1
1796 do id=1,ndm
1797 gfr(id,iz)=gfrb(id,ibui)
1798 tpvlev(id,iz)=tpvlevb(id,ibui)
1799 sfr(id,iz)=sfrb(id,ibui)
1800 hfgr(id,iz)=hfgrb(id,ibui)
1801 sfrpv(id,iz)=-sfrbpv(id,ibui)
1802 lfr(id,iz)=lfrb(id,ibui)
1803 do ily=1,nwr_u
1804 tr(id,iz,ily)=trb(id,ily,ibui)
1805 enddo
1806 ptr(id,iz)=tr(id,iz,nwr_u)*(pr_u(iz)/p0)**(-rcp_u)
1807 enddo
1808 endif
1809 enddo !iz
1811 !Compute the potential temperature for the vertical surfaces of the buildings
1813 do id=1,ndm
1814 do iz=1,nzu!nz_um
1815 do ibui=1,nbui
1816 ptw(2*id-1,iz,ibui)=tw(2*id-1,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u)
1817 ptw(2*id,iz,ibui)=tw(2*id,iz,nwr_u,ibui)*(pr_u(iz)/p0)**(-rcp_u)
1818 ptwin(2*id-1,iz,ibui)=twlev(2*id-1,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u)
1819 ptwin(2*id,iz,ibui)=twlev(2*id,iz,ibui)*(pr_u(iz)/p0)**(-rcp_u)
1821 enddo
1822 enddo
1823 enddo
1824 !NEW CDRAG!
1825 do iz=1,nz_um
1826 alp=0.
1827 do id=1,ndu
1828 alp=alp+bs(id)/(ws(id)+bs(id))*pb(iz)
1829 enddo
1830 alp=alp/ndu
1831 if(alp.lt.0.29)then
1832 cdrag(iz)=3.32*alp**0.47
1833 else
1834 cdrag(iz)=1.85
1835 endif
1836 enddo
1840 ! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
1842 call buildings(iurb,ndu,nzu,z0,cdrag,ua_u,va_u, &
1843 pt_u,pt0_u,ptg,ptr,ptrv,da_u,qv_u,pr_u,tmp_u,ptw,ptwin,pwin_u(iurb),drst, &
1844 uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,qvb_u,qhb_u, &
1845 uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr,sfrpv,sfrv,lfrv, &
1846 dgr,dg,lfr,lfg, &
1847 sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac,ix,iy,rsg,rs,qr,gr_frac_roof, &
1848 pv_frac_roof,gr_flag,gr_type)
1853 ! Calculation of the sensible heat fluxes for the ground, the wall and roof
1854 ! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
1855 ! where A and B are the implicit and explicit components of the heat sources or sinks.
1857 ! Interpolation on the "mesoscale grid"
1859 call urban_meso(ndu,kms,kme,kts,kte,nzu,z,dz,z_u,pb,ss,bs,ws,sf, &
1860 vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
1861 uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u, &
1862 a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)
1865 ! Calculation of the length scale taking into account the buildings effects
1867 call interp_length(ndu,kms,kme,kts,kte,nzu,z_u,z,ss,ws,bs,dlg,dl_u)
1869 return
1870 end subroutine BEP1D
1872 ! ===6=8===============================================================72
1873 ! ===6=8===============================================================72
1875 subroutine param(iurb,nzu,nzurb,nzurban,ndu, &
1876 csg_u,csg,alag_u,alag,csr_u,csr, &
1877 alar_u,alar,csw_u,csw,alaw_u,alaw, &
1878 ws_u,ws_urb,ws,bs_u,bs_urb,bs,z0g_u,z0r_u,z0, &
1879 strd_u,strd,drst_u,drst,ss_u,ss_urb,ss,pb_u, &
1880 pb_urb,pb,dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb,&
1881 lp_urb,lb_urb,hgt_urb,frc_urb)
1883 ! ----------------------------------------------------------------------
1884 ! This routine prepare some usefull parameters
1885 ! ----------------------------------------------------------------------
1887 implicit none
1890 ! ----------------------------------------------------------------------
1891 ! INPUT:
1892 ! ----------------------------------------------------------------------
1893 integer iurb ! Current urban class
1894 integer nzu ! Number of vertical urban levels in the current class
1895 integer ndu ! Number of street direction for the current urban class
1896 integer nzurb ! Number of vertical urban levels in the current class
1897 real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
1898 real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
1899 real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
1900 real bs_u(ndm,nurbm) ! Building width
1901 real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
1902 real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
1903 real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
1904 real drst_u(ndm,nurbm) ! Street direction
1905 real strd_u(ndm,nurbm) ! Street length
1906 real ws_u(ndm,nurbm) ! Street width
1907 real z0g_u(nurbm) ! The ground's roughness length
1908 real z0r_u(nurbm) ! The roof's roughness length
1909 real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
1910 real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
1911 real lp_urb ! Building plan area density
1912 real lb_urb ! Building surface area to plan area ratio
1913 real hgt_urb ! Average building height weighted by building plan area [m]
1914 real frc_urb ! Urban fraction
1916 ! ----------------------------------------------------------------------
1917 ! OUTPUT:
1918 ! ----------------------------------------------------------------------
1919 real alag(ng_u) ! Ground thermal diffusivity at each ground levels
1920 real csg(ng_u) ! Specific heat of the ground material at each ground levels
1921 real bs(ndm) ! Building width for the current urban class
1922 real drst(ndm) ! street directions for the current urban class
1923 real strd(ndm) ! Street lengths for the current urban class
1924 real ws(ndm) ! Street widths of the current urban class
1925 real z0(ndm,nz_um) ! Roughness lengths "profiles"
1926 real ss(nz_um) ! Probability to have a building with height h
1927 real pb(nz_um) ! Probability to have a building with an height greater or equal to "z"
1928 integer nzurban
1930 !-----------------------------------------------------------------------------
1931 !INPUT/OUTPUT
1932 !-----------------------------------------------------------------------------
1934 real dzw(nwr_u) !Layer sizes in the walls [m]
1935 real dzr(nwr_u) !Layer sizes in the roofs [m]
1936 real dzf(nf_u) !Layer sizes in the floors [m]
1937 real dzgb(ngb_u) !layer sizes in the ground below the buildings [m]
1939 real csr(nwr_u) ! Specific heat of the roof material at each roof levels
1940 real csw(nwr_u) ! Specific heat of the wall material at each wall levels
1942 real csf(nf_u) !Specific heat of the floors materials in the buildings
1943 !of the current urban class [J m^3 K^-1]
1944 real csgb(ngb_u) !Specific heat of the ground material below the buildings
1945 !of the current urban class [J m^3 K^-1]
1946 real alar(nwr_u+1) ! Roof thermal diffusivity at each roof levels [W/ m K]
1947 real alaw(nwr_u+1) ! Wall thermal diffusivity at each wall levels [W/ m K]
1948 real alaf(nf_u+1) ! Floor thermal diffusivity at each wall levels [W/m K]
1949 real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall levels [W/m K]
1950 real bs_urb(ndm,nurbm) ! Building width
1951 real ws_urb(ndm,nurbm) ! Street width
1952 real ss_urb(nz_um,nurbm) ! The probability that a building has an height equal to "z"
1953 real pb_urb(nz_um) ! Probability that a building has an height greater or equal to z
1954 ! ----------------------------------------------------------------------
1955 ! LOCAL:
1956 ! ----------------------------------------------------------------------
1957 integer id,ig,ir,iw,iz,iflo,ihu
1958 ! ----------------------------------------------------------------------
1959 ! END VARIABLES DEFINITIONS
1960 ! ----------------------------------------------------------------------
1961 !
1962 !Initialize variables
1963 !
1964 ss=0.
1965 pb=0.
1966 csg=0.
1967 alag=0.
1968 csgb=0.
1969 alagb=0.
1970 csf=0.
1971 alaf=0.
1972 csr=0.
1973 alar=0.
1974 csw=0.
1975 alaw=0.
1976 z0=0.
1977 ws=0.
1978 bs=0.
1979 bs_urb=0.
1980 ws_urb=0.
1981 strd=0.
1982 drst=0.
1983 nzurban=0
1985 !Define the layer sizes in the walls
1987 dzgb=(/0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/)
1988 dzr=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
1989 dzw=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
1990 dzf=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02/)
1992 ihu=0
1994 do iz=1,nz_um
1995 if (ss_urb(iz,iurb)/=0.) then
1996 ihu=1
1997 exit
1998 else
1999 continue
2000 endif
2001 enddo
2003 if (ihu==1) then
2004 do iz=1,nzurb+1
2005 ss(iz)=ss_urb(iz,iurb)
2006 pb(iz)=pb_urb(iz)
2007 enddo
2008 nzurban=nzurb
2009 else
2010 do iz=1,nzu+1
2011 ss(iz)=ss_u(iz,iurb)
2012 pb(iz)=pb_u(iz,iurb)
2013 ss_urb(iz,iurb)=ss_u(iz,iurb)
2014 pb_urb(iz)=pb_u(iz,iurb)
2015 end do
2016 nzurban=nzu
2017 endif
2019 do ig=1,ngb_u
2020 csgb(ig) = csg_u(iurb)
2021 alagb(ig)= csg_u(iurb)*alag_u(iurb)
2022 enddo
2023 alagb(ngb_u+1)= csg_u(iurb)*alag_u(iurb)
2025 do iflo=1,nf_u
2026 csf(iflo) = csw_u(iurb)
2027 alaf(iflo)= csw_u(iurb)*alaw_u(iurb)
2028 enddo
2029 alaf(nf_u+1)= csw_u(iurb)*alaw_u(iurb)
2031 do ir=1,nwr_u
2032 csr(ir) = csr_u(iurb)
2033 alar(ir)= csr_u(iurb)*alar_u(iurb)
2034 enddo
2035 alar(nwr_u+1)= csr_u(iurb)*alar_u(iurb)
2037 do iw=1,nwr_u
2038 csw(iw) = csw_u(iurb)
2039 alaw(iw)= csw_u(iurb)*alaw_u(iurb)
2040 enddo
2041 alaw(nwr_u+1)=csw_u(iurb)*alaw_u(iurb)
2043 !------------------------------------------------------------------------
2045 do ig=1,ng_u
2046 csg(ig)=csg_u(iurb)
2047 alag(ig)=alag_u(iurb)
2048 enddo
2050 do id=1,ndu
2051 z0(id,1)=z0g_u(iurb)
2052 do iz=2,nzurban+1
2053 z0(id,iz)=z0r_u(iurb)
2054 enddo
2055 enddo
2057 do id=1,ndu
2058 strd(id)=strd_u(id,iurb)
2059 drst(id)=drst_u(id,iurb)
2060 enddo
2062 do id=1,ndu
2063 if ((hgt_urb<=0.).OR.(lp_urb<=0.).OR.(lb_urb<=0.)) then
2064 ws(id)=ws_u(id,iurb)
2065 bs(id)=bs_u(id,iurb)
2066 bs_urb(id,iurb)=bs_u(id,iurb)
2067 ws_urb(id,iurb)=ws_u(id,iurb)
2068 else if ((lp_urb/frc_urb<1.).and.(lp_urb<lb_urb)) then
2069 bs(id)=2.*hgt_urb*lp_urb/(lb_urb-lp_urb)
2070 ws(id)=2.*hgt_urb*lp_urb*((frc_urb/lp_urb)-1.)/(lb_urb-lp_urb)
2071 bs_urb(id,iurb)=bs(id)
2072 ws_urb(id,iurb)=ws(id)
2073 else
2074 ws(id)=ws_u(id,iurb)
2075 bs(id)=bs_u(id,iurb)
2076 bs_urb(id,iurb)=bs_u(id,iurb)
2077 ws_urb(id,iurb)=ws_u(id,iurb)
2078 endif
2079 enddo
2080 do id=1,ndu
2081 if ((bs(id)<=1.).OR.(bs(id)>=150.)) then
2082 bs(id)=bs_u(id,iurb)
2083 ws(id)=ws_u(id,iurb)
2084 bs_urb(id,iurb)=bs_u(id,iurb)
2085 ws_urb(id,iurb)=ws_u(id,iurb)
2086 endif
2087 if ((ws(id)<=1.).OR.(ws(id)>=150.)) then
2088 ws(id)=ws_u(id,iurb)
2089 bs(id)=bs_u(id,iurb)
2090 bs_urb(id,iurb)=bs_u(id,iurb)
2091 ws_urb(id,iurb)=ws_u(id,iurb)
2092 endif
2093 enddo
2094 return
2095 end subroutine param
2097 ! ===6=8===============================================================72
2098 ! ===6=8===============================================================72
2100 subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
2102 ! ----------------------------------------------------------------------
2103 ! This routine interpolate para
2104 ! meters from the "mesoscale grid" to
2105 ! the "urban grid".
2106 ! See p300 Appendix B.1 of the BLM paper.
2107 ! ----------------------------------------------------------------------
2109 implicit none
2111 ! ----------------------------------------------------------------------
2112 ! INPUT:
2113 ! ----------------------------------------------------------------------
2114 ! Data relative to the "mesoscale grid"
2115 integer kts,kte,kms,kme
2116 real z(kms:kme) ! Altitude of the cell interface
2117 real c(kms:kme) ! Parameter which has to be interpolated
2118 ! Data relative to the "urban grid"
2119 integer nz_u ! Number of levels
2120 !! real z_u(nz_u+1) ! Altitude of the cell interface
2121 real z_u(nz_um) ! Altitude of the cell interface
2123 ! ----------------------------------------------------------------------
2124 ! OUTPUT:
2125 ! ----------------------------------------------------------------------
2126 !! real c_u(nz_u) ! Interpolated paramters in the "urban grid"
2127 real c_u(nz_um) ! Interpolated paramters in the "urban grid"
2129 ! LOCAL:
2130 ! ----------------------------------------------------------------------
2131 integer iz_u,iz
2132 real ctot,dz
2134 ! ----------------------------------------------------------------------
2135 ! END VARIABLES DEFINITIONS
2136 ! ----------------------------------------------------------------------
2138 do iz_u=1,nz_u
2139 ctot=0.
2140 do iz=kts,kte
2141 dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
2142 ctot=ctot+c(iz)*dz
2143 enddo
2144 c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
2145 enddo
2147 return
2148 end subroutine interpol
2150 ! ===6=8===============================================================72
2151 ! ===6=8===============================================================72
2153 subroutine averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av, &
2154 sfw_av,sfwind_av,sfw,sfwin)
2156 implicit none
2157 !
2158 !INPUT VARIABLES
2159 !
2160 real tw(2*ndm,nz_um,nwr_u,nbui_max) ! Temperature in each layer of the wall [K]
2161 real twlev(2*ndm,nz_um,nbui_max) ! Window temperature in BEM [K]
2162 real pb(nz_um) ! Probability to have a building with an height equal or greater h
2163 real ss(nz_um) ! Probability to have a building with height h
2164 real sfw(2*ndm,nz_um,nbui_max) ! Surface fluxes from the walls
2165 real sfwin(2*ndm,nz_um,nbui_max) ! Surface fluxes from the windows
2166 !
2167 !OUTPUT VARIABLES
2168 !
2169 real tw_av(2*ndm,nz_um) ! Averaged temperature of the wall surfaces
2170 real twlev_av(2*ndm,nz_um) ! Averaged temperature of the windows
2171 real sfw_av(2*ndm,nz_um) ! Averaged sensible heat from walls
2172 real sfwind_av(2*ndm,nz_um) ! Averaged sensible heat from windows
2173 !
2174 !LOCAL VARIABLES
2175 !
2176 real d_urb(nz_um)
2177 integer nlev(nz_um)
2178 integer id,iz
2179 integer nbui,ibui
2180 !
2181 !Initialize Variables
2182 !
2183 tw_av=0.
2184 twlev_av=0.
2185 sfw_av=0.
2186 sfwind_av=0.
2187 ibui=0
2188 nbui=0
2189 nlev=0
2190 d_urb=0.
2192 do iz=1,nz_um
2193 if(ss(iz).gt.0) then
2194 ibui=ibui+1
2195 d_urb(ibui)=ss(iz)
2196 nlev(ibui)=iz-1
2197 nbui=ibui
2198 endif
2199 enddo
2201 do id=1,ndm
2202 do iz=1,nz_um-1
2203 if (pb(iz+1).gt.0) then
2204 do ibui=1,nbui
2205 if (iz.le.nlev(ibui)) then
2206 tw_av(2*id-1,iz)=tw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
2207 tw(2*id-1,iz,nwr_u,ibui)**4
2208 tw_av(2*id,iz)=tw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
2209 tw(2*id,iz,nwr_u,ibui)**4
2210 twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
2211 twlev(2*id-1,iz,ibui)**4
2212 twlev_av(2*id,iz)=twlev_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*&
2213 twlev(2*id,iz,ibui)**4
2214 sfw_av(2*id-1,iz)=sfw_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id-1,iz,ibui)
2215 sfw_av(2*id,iz)=sfw_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfw(2*id,iz,ibui)
2216 sfwind_av(2*id-1,iz)=sfwind_av(2*id-1,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id-1,iz,ibui)
2217 sfwind_av(2*id,iz)=sfwind_av(2*id,iz)+(d_urb(ibui)/pb(iz+1))*sfwin(2*id,iz,ibui)
2218 endif
2219 enddo
2220 tw_av(2*id-1,iz)=tw_av(2*id-1,iz)**(1./4.)
2221 tw_av(2*id,iz)=tw_av(2*id,iz)**(1./4.)
2222 twlev_av(2*id-1,iz)=twlev_av(2*id-1,iz)**(1./4.)
2223 twlev_av(2*id,iz)=twlev_av(2*id,iz)**(1./4.)
2224 endif
2225 enddo !iz
2226 enddo !id
2227 return
2228 end subroutine averaging_temp
2229 ! ===6=8===============================================================72
2230 ! ===6=8===============================================================72
2232 subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, &
2233 tw,tg_av,twlev,albg,albw,emw,emg,pwin,albwin, &
2234 emwin,fww,fwg,fgw,fsw,fsg, &
2235 zr,deltar,ah,xlat,swddir,swddif, &
2236 rs,rl,rsw,rsd,rsg,rlw,rlg)
2238 ! ----------------------------------------------------------------------
2239 ! This routine computes the modification of the short wave and
2240 ! long wave radiation due to the buildings.
2241 ! ----------------------------------------------------------------------
2243 implicit none
2246 ! ----------------------------------------------------------------------
2247 ! INPUT:
2248 ! ----------------------------------------------------------------------
2249 integer iurb ! current urban class
2250 integer nd ! Number of street direction for the current urban class
2251 integer nz_u ! Number of layer in the urban grid
2252 real z(nz_um) ! Height of the urban grid levels
2253 real ws(ndm) ! Street widths of the current urban class
2254 real drst(ndm) ! street directions for the current urban class
2255 real strd(ndm) ! Street lengths for the current urban class
2256 real ss(nz_um) ! probability to have a building with height h
2257 real pb(nz_um) ! probability to have a building with an height equal
2258 real tw(2*ndm,nz_um) ! Temperature in each layer of the wall [K]
2259 real tg_av(ndm) ! Temperature in each layer of the ground [K]
2260 real albg ! Albedo of the ground for the current urban class
2261 real albw ! Albedo of the wall for the current urban class
2262 real emg ! Emissivity of ground for the current urban class
2263 real emw ! Emissivity of wall for the current urban class
2264 real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
2265 real fsg(ndm,nurbm) ! View factors from sky to ground
2266 real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
2267 real fws(nz_um,ndm,nurbm) ! View factors from wall to sky
2268 real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
2269 real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
2270 real ah ! Hour angle (it should come from the radiation routine)
2271 real zr ! zenith angle
2272 real deltar ! Declination of the sun
2273 real rs ! solar radiation
2274 real rl ! downward flux of the longwave radiation
2275 real xlat ! latitudine
2276 real swddir ! short wave direct solar radiation _gl
2277 real swddif ! short wave diffuse solar radiation _gl
2279 !
2280 !New variables BEM
2281 !
2282 real twlev(2*ndm,nz_um) ! Window temperature in BEM [K]
2283 real pwin ! Coverage area fraction of windows in the walls of the buildings
2284 real albwin ! Albedo of the windows for the current urban class
2285 real emwin ! Emissivity of the windows for the current urban class
2286 real alb_av ! Averaged albedo (window and wall)
2287 ! ----------------------------------------------------------------------
2288 ! OUTPUT:
2289 ! ----------------------------------------------------------------------
2290 real rlg(ndm) ! Long wave radiation at the ground
2291 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
2292 real rsg(ndm) ! Short wave radiation at the ground
2293 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
2294 real rsd(2*ndm,nz_um) ! Direct Short wave radiation at the walls
2296 ! ----------------------------------------------------------------------
2297 ! LOCAL:
2298 ! ----------------------------------------------------------------------
2300 integer id,iz
2302 ! Calculation of the shadow effects
2304 call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
2305 swddir,rsw,rsg,xlat)
2306 rsd=rsw
2308 ! Calculation of the reflection effects
2309 do id=1,nd
2310 call long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev, &
2311 fwg,fww,fgw,fsw,fsg,tg_av,tw,rlg,rlw,rl,pb)
2313 alb_av=pwin*albwin+(1.-pwin)*albw
2315 call short_rad_dd(iurb,nz_u,id,alb_av, &
2316 albg,swddif,fwg,fww,fgw,fsw,fsg,rsg,rsw,pb)
2319 enddo
2320 return
2321 end subroutine modif_rad
2325 ! ===6=8===============================================================72
2326 ! ===6=8===============================================================72
2328 subroutine surf_temp(nd,pr,dt,rl,rsg,rlg, &
2329 tg,alag,csg,emg,albg,ptg,sfg,lfg,gfg)
2331 ! ----------------------------------------------------------------------
2332 ! Computation of the surface temperatures for walls, ground and roofs
2333 ! ----------------------------------------------------------------------
2335 implicit none
2337 ! ----------------------------------------------------------------------
2338 ! INPUT:
2339 ! ----------------------------------------------------------------------
2341 integer nd ! Number of street direction for the current urban class
2342 real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
2344 real albg ! Albedo of the ground for the current urban class
2346 real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
2348 real dt ! Time step
2349 real emg ! Emissivity of ground for the current urban class
2351 real pr(nz_um) ! Air pressure
2353 real rl ! Downward flux of the longwave radiation
2354 real rlg(ndm) ! Long wave radiation at the ground
2356 real rsg(ndm) ! Short wave radiation at the ground
2358 real sfg(ndm) ! Sensible heat flux from ground (road)
2360 real lfg(ndm) ! Latent heat flux from ground (road)
2362 real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) toward the interior
2364 real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
2366 ! ----------------------------------------------------------------------
2367 ! OUTPUT:
2368 ! ----------------------------------------------------------------------
2369 real ptg(ndm) ! Ground potential temperatures
2371 ! ----------------------------------------------------------------------
2372 ! LOCAL:
2373 ! ----------------------------------------------------------------------
2374 integer id,ig,ir,iw,iz
2376 real rtg(ndm) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
2378 real tg_tmp(ng_u)
2380 real dzg_u(ng_u) ! Layer sizes in the ground
2382 data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
2385 ! ----------------------------------------------------------------------
2386 ! END VARIABLES DEFINITIONS
2387 ! ----------------------------------------------------------------------
2391 do id=1,nd
2393 ! Calculation for the ground surfaces
2394 do ig=1,ng_u
2395 tg_tmp(ig)=tg(id,ig)
2396 end do
2397 !
2398 ! print*,'alag','cs',alag(1),csg(1)
2400 call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg, &
2401 rsg(id),rlg(id),pr(1), &
2402 dt,emg,albg, &
2403 rtg(id),sfg(id),lfg(id),gfg(id))
2405 do ig=1,ng_u
2406 tg(id,ig)=tg_tmp(ig)
2407 end do
2409 end do !id
2411 return
2412 end subroutine surf_temp
2415 ! ===6=8===============================================================72
2416 ! ===6=8===============================================================72
2419 subroutine roof_temp_veg(nd,pr,dt,rl,rsr, &
2420 trv,ptrv,sfrv,lfrv,gfr,qr,rainbl,drain,hfgroof,tr,alar,dzr,csr,nzu,irri_now,gr_type,pv_frac_roof,tpvlev)
2422 ! ----------------------------------------------------------------------
2423 ! Computation of the surface temperatures for walls, ground and roofs
2424 ! ----------------------------------------------------------------------
2426 implicit none
2428 ! ----------------------------------------------------------------------
2429 ! INPUT:
2430 ! ----------------------------------------------------------------------
2431 real rainbl
2432 integer nd ! Number of street direction for the current urban class
2434 integer nzu ! Number of urban layers
2435 real irho(24) ! Which hour of irrigation\
2438 real alar ! Roof thermal diffusivity for the current urban class [m^2 s^-1]
2439 real pv_frac_roof
2440 real csr
2442 real dzr ! Layer sizes in the roofs [m]
2444 real dt ! Time step
2446 real pr(nz_um) ! Air pressure
2448 real rl ! Downward flux of the longwave radiation
2450 real rsr ! Short wave radiation at the ground
2452 real tpvlev(ndm,nz_um)
2454 real sfrv(ndm,nz_um) ! Sensible heat flux from ground (road)
2456 real lfrv(ndm,nz_um) ! Latent heat flux from ground (road)
2458 real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the ground (road) toward the interior
2460 real trv(ndm,nz_um,ngr_u) ! Temperature in each layer of the green roof [K]
2462 real qr(ndm,nz_um,ngr_u) ! Humidity in each layer of the green roof
2464 real tr(ndm,nz_um,nwr_u) !Roof temperature in BEM [K]
2466 ! ----------------------------------------------------------------------
2467 ! OUTPUT:
2468 ! ----------------------------------------------------------------------
2469 real ptrv(ndm,nz_um) ! Ground potential temperatures
2471 real hfgroof(ndm,nz_um)
2472 ! ----------------------------------------------------------------------
2473 ! LOCAL:
2474 ! ----------------------------------------------------------------------
2475 integer id,ig,ir,iw,iz
2477 real alagr(ngr_u) ! Green Roof thermal diffusivity for the current urban class [m^2 s^-1]
2479 real rtr(ndm,nz_um) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
2481 real tr_tmp(ngr_u)
2483 real qr_tmp(ngr_u)
2484 real qr_tmp_old(ngr_u)
2485 real dzgr_u(ngr_u) ! Layer sizes in the green roof
2486 !MODIFICA
2487 data dzgr_u /0.1,0.003,0.06,0.003,0.05,0.04,0.02,0.0125,0.005,0.0025/
2488 real cs(ngr_u) ! Specific heat of the ground material
2489 real cw
2490 parameter(cw=4.295e6)
2491 real s(ngr_u)
2492 real d(ngr_u)
2493 real k(ngr_u)
2494 real qr_m ! mean soil moisture between layers
2495 real qrmax(ngr_u)
2496 real smax(ngr_u)
2497 real kmax(ngr_u)
2498 real b(ngr_u)
2499 real cd(ngr_u)
2500 real csa(4)
2501 real ka(4)
2502 real qref
2503 parameter(qref=0.37)
2504 data qrmax /0.0,0.0,0.0,0.0,0.439,0.37,0.37,0.37,0.37,0.37/
2505 data smax /0,0,0,0,-0.01,-0.1,-0.1,-0.1,-0.1,-0.1/
2506 data kmax /0,0,0,0,3.32e-3,2.162e-3,2.162e-3,2.162e-3,2.162e-3,2.162e-3/
2507 data b /0,0,0,0,2.7,3.9,3.9,3.9,3.9,3.9/
2508 data cd /0,0,0,0,331500,1.342e6,1.342e6,1.342e6,1.342e6,1.342e6/
2509 data csa /7.5e4,2.1e6,4.48e4,2.1e6/
2510 data ka /0.035,0.7,0.024,0.7/
2511 real em_gr(1)
2512 real alb_gr(1)
2513 real irri_now
2514 integer gr_type
2515 real drain(ndm,nz_um)
2516 ! ----------------------------------------------------------------------
2517 ! END VARIABLES DEFINITIONS
2519 if(gr_type.eq.1)then
2520 em_gr=0.95
2521 alb_gr=0.3
2522 elseif(gr_type.eq.2)then
2523 em_gr=0.83
2524 alb_gr=0.154
2525 endif
2528 do iz=2,nzu
2530 do id=1,nd
2535 ! Calculation for the ground surfaces
2537 do ig=1,ngr_u
2538 tr_tmp(ig)=trv(id,iz,ig)
2539 qr(id,iz,ig) = max(qr(id,iz,ig),1e-6) !cenlin, 11/4/2020
2540 qr_tmp(ig)=qr(id,iz,ig)
2541 qr_tmp_old(ig)=qr(id,iz,ig)
2543 if(ig.le.4) then
2545 cs(ig)=csa(ig)
2546 alagr(ig)=ka(ig)/csa(ig)
2548 else
2551 if (ig.gt.5) then
2552 qr_m=(qr(id,iz,ig)*dzgr_u(ig-1)+qr(id,iz,ig-1)*dzgr_u(ig))/(dzgr_u(ig)+dzgr_u(ig-1))
2553 else
2554 qr_m=qr(id,iz,ig)
2555 endif
2556 cs(ig)=(1-qr_m)*cd(ig)+qr_m*cw
2557 s(ig)=smax(ig)*(qrmax(ig)/qr_m)**b(ig)
2558 k(ig)=kmax(ig)*(qr_m/qrmax(ig))**(2*b(ig)+3)
2559 d(ig)=-b(ig)*kmax(ig)*smax(ig)*((qr_m/qrmax(ig))**(b(ig)+3))/qr_m
2560 if (log10(abs(s(ig))).le.5.1) then
2561 alagr(ig)=exp(-(log10(abs(s(ig)))+2.7))*4.186e2/cs(ig)
2562 endif
2563 if (log10(abs(s(ig))).gt.5.1) then
2564 alagr(ig)=0.00041*4.186e2/cs(ig)
2565 endif
2567 endif
2569 end do
2570 hfgroof(id,iz)=(alar/csr+alagr(1))*(tr_tmp(1)-tr(id,iz,5))/(dzr+dzgr_u(1))
2572 call soil_temp_veg(hfgroof(id,iz),ngr_u,dzgr_u,tr_tmp,ptrv(id,iz),alagr,cs, &
2573 rsr,rl,pr(iz), &
2574 dt,em_gr(1),alb_gr(1), &
2575 rtr(id,iz),sfrv(id,iz),lfrv(id,iz),gfr(id,iz),pv_frac_roof,tpvlev(id,iz))
2576 do ig=1,ngr_u
2577 trv(id,iz,ig)=tr_tmp(ig)
2578 end do
2579 drain(id,iz)=kmax(5)*(qr(id,iz,5)/qrmax(5))**(2*b(5)+3)
2580 call soil_moist(ngr_u,dzgr_u,qr_tmp,dt,lfrv(id,iz),d,k,rainbl,drain(id,iz),irri_now)
2583 do ig=1,ngr_u
2584 ! qr(id,iz,ig)=min(qr_tmp(ig),qrmax(ig))
2585 qr(id,iz,ig)=max(min(qr_tmp(ig),qrmax(ig)),1e-6) !cenlin,11/4/2020
2586 end do
2588 end do !id
2589 end do !iz
2591 return
2592 end subroutine roof_temp_veg
2594 ! ===6=8===============================================================72
2595 ! ===6=8===============================================================72
2597 subroutine buildings(iurb,nd,nz,z0,cdrag,ua_u,va_u,pt_u,pt0_u, &
2598 ptg,ptr,ptrv,da_u,qv_u,pr_u,tmp_u,ptw,ptwin,pwin, &
2599 drst,uva_u,vva_u,uvb_u,vvb_u, &
2600 tva_u,tvb_u,evb_u,qvb_u,qhb_u, &
2601 uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr,sfrpv,sfrv,lfrv, &
2602 dgr,dg,lfr,lfg, &
2603 sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac,ix,iy,rsg,rs,qr,gr_frac_roof, &
2604 pv_frac_roof,gr_flag,gr_type)
2606 ! ----------------------------------------------------------------------
2607 ! This routine computes the sources or sinks of the different quantities
2608 ! on the urban grid. The actual calculation is done in the subroutines
2609 ! called flux_wall, and flux_flat.
2610 ! ----------------------------------------------------------------------
2612 implicit none
2615 ! ----------------------------------------------------------------------
2616 ! INPUT:
2617 ! ----------------------------------------------------------------------
2618 integer nd ! Number of street direction for the current urban class
2619 integer ix,iy
2620 integer nz ! number of vertical space steps
2621 real ua_u(nz_um) ! Wind speed in the x direction on the urban grid
2622 real va_u(nz_um) ! Wind speed in the y direction on the urban grid
2623 real da_u(nz_um) ! air density on the urban grid
2624 real qv_u(nz_um) ! specific humidity on the urban grid
2625 real pr_u(nz_um) ! pressure on the urban grid
2626 real tmp_u(nz_um) ! temperaure on the urban grid
2627 real drst(ndm) ! Street directions for the current urban class
2628 real dz
2629 real pt_u(nz_um) ! Potential temperature on the urban grid
2630 real pt0_u(nz_um) ! reference potential temperature on the urban grid
2631 real ptg(ndm) ! Ground potential temperatures
2632 real ptr(ndm,nz_um) ! Roof potential temperatures
2633 real ptrv(ndm,nz_um) ! Green Roof potential temperatures
2634 real ptw(2*ndm,nz_um,nbui_max) ! Walls potential temperatures
2635 real ss(nz_um) ! probability to have a building with height h
2636 real pb(nz_um)
2637 real cdrag(nz_um)
2638 real z0(ndm,nz_um) ! Roughness lengths "profiles"
2639 real dt ! time step
2640 integer iurb !Urban class
2641 real rsg(ndm) ! Solar Radiation
2642 real rs ! Solar Radiation
2643 real qr(ndm,nz_um,ngr_u) ! Ground Soil Moisture
2644 real trv(ndm,nz_um,ngr_u) ! Ground Soil Moisture
2645 real roof_frac
2646 real road_frac
2647 !
2648 !New variables (BEM)
2649 !
2650 real bs_u(ndm,nurbm) ! Building width [m]
2651 real dz_u ! Urban grid resolution
2652 real sflev(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems [W]
2653 real lflev(nz_um,nz_um) ! latent heat flux due to the air conditioning systems [W]
2654 real sfvlev(nz_um,nz_um) ! sensible heat flux due to ventilation [W]
2655 real lfvlev(nz_um,nz_um) ! latent heat flux due to ventilation [W]
2656 real qvb_u(2*ndm,nz_um)
2657 real qhb_u(ndm,nz_um)
2658 real ptwin(2*ndm,nz_um,nbui_max) ! window potential temperature
2659 real pwin
2660 real tvb_ac(2*ndm,nz_um)
2661 real gr_frac_roof
2662 real pv_frac_roof
2663 integer gr_flag,gr_type
2665 ! ----------------------------------------------------------------------
2666 ! OUTPUT:
2667 ! ----------------------------------------------------------------------
2668 ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
2669 ! vertical surfaces (walls) and horizontal surfaces (roofs and street)
2670 ! The fluxes can be computed as follow: Fluxes of X = A*X + B
2671 ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
2673 real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
2674 real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
2675 real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
2676 real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
2677 real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
2678 real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
2679 real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
2680 real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
2681 real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
2682 real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
2683 real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
2684 real uhb(2*ndm,nz_um)
2685 real vhb(2*ndm,nz_um)
2686 real ehb(2*ndm,nz_um)
2687 real sfw(2*ndm,nz_um,nbui_max) ! sensible heat flux from walls
2688 real sfwin(2*ndm,nz_um,nbui_max) ! sensible heat flux form windows
2689 real sfr(ndm,nz_um) ! sensible heat flux from roof
2690 real sfrv(ndm,nz_um) ! sensible heat flux from roof
2691 real lfrv(ndm,nz_um) ! Latent heat flux from roof
2692 real dgr(ndm,nz_um) ! sensible heat flux from roof
2693 real dg(ndm)
2694 real lfr(ndm,nz_um) ! Latent heat flux from roof
2695 real lfg(ndm) ! Latent heat flux from street
2696 real sfrpv(ndm,nz_um) ! sensible heat flux from PV panels
2697 real sfg(ndm) ! sensible heat flux from street
2700 ! ----------------------------------------------------------------------
2701 ! LOCAL:
2702 ! ----------------------------------------------------------------------
2703 real d_urb(nz_um)
2704 real uva_tmp
2705 real vva_tmp
2706 real uvb_tmp
2707 real vvb_tmp
2708 real evb_tmp
2709 integer nlev(nz_um)
2710 integer id,iz,ibui,nbui,il
2711 real wfg !Ground water pool fraction
2712 real wfr !Roof water pool fraction
2713 real uhbv(2*ndm,nz_um)
2714 real vhbv(2*ndm,nz_um)
2715 real ehbv(2*ndm,nz_um)
2716 real z0v !Vegetation roughness
2717 parameter(z0v=0.01)
2718 real resg
2719 real rsveg
2720 real f1,f2,f3,f4
2721 integer rsv(2)
2722 real qr_tmp(ngr_u)
2723 data rsv /0,1/
2724 real fh,ric,utot
2725 !------------------------------------------------------------------
2726 ! END VARIABLES DEFINITIONS
2727 ! ----------------------------------------------------------------------
2728 dz=dz_u
2729 ibui=0
2730 nbui=0
2731 nlev=0
2732 d_urb=0.
2734 uva_u=0.
2735 uvb_u=0.
2736 vhb_u=0.
2737 vva_u=0.
2738 vvb_u=0.
2739 thb_u=0.
2740 tva_u=0.
2741 tvb_u=0.
2742 tvb_ac=0.
2743 ehb_u=0.
2744 evb_u=0.
2745 qvb_u=0.
2746 qhb_u=0.
2748 uhb=0.
2749 vhb=0.
2750 ehb=0.
2751 uhbv=0.
2752 vhbv=0.
2753 ehbv=0.
2756 do iz=1,nz_um
2757 if(ss(iz).gt.0)then
2758 ibui=ibui+1
2759 d_urb(ibui)=ss(iz)
2760 nlev(ibui)=iz-1
2761 nbui=ibui
2762 endif
2763 enddo
2765 ! Calculation at the ground surfaces
2766 do id=1,nd
2768 call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1), &
2769 ptg(id),qv_u(1),uhb(id,1), &
2770 vhb(id,1),sfg(id),lfg(id),ehb(id,1),da_u(1),pr_u(1))
2771 if(dg(id).gt.0)then
2772 wfg=dg(id)/dgmax
2773 lfg(id)=-da_u(1)*latent*(-(wfg*lfg(id))/(da_u(1)*latent))
2774 else
2775 qhb_u(id,1)=0.
2776 lfg(id)=0.
2777 endif
2778 thb_u(id,1)=-(sfg(id))/(da_u(1)*cp_u)
2779 vhb_u(id,1)=vhb(id,1)
2780 uhb_u(id,1)=uhb(id,1)
2781 ehb_u(id,1)=ehb(id,1)
2782 qhb_u(id,1)=-lfg(id)/(da_u(1)*latent)
2783 do iz=2,nz
2784 if(ss(iz).gt.0)then
2786 call flux_flat(dz,z0(id,iz),ua_u(iz),&
2787 va_u(iz),pt_u(iz),pt0_u(iz), &
2788 ptr(id,iz),qv_u(iz),uhb(id,iz), &
2789 vhb(id,iz),sfr(id,iz),lfr(id,iz),ehb(id,iz),da_u(iz),pr_u(iz))
2790 if(dgr(id,iz).gt.0)then
2791 wfr=dgr(id,iz)/drmax
2792 lfr(id,iz)=-da_u(iz)*latent*(-(wfr*lfr(id,iz))/(da_u(iz)*latent))
2793 else
2794 lfr(id,iz)=0.
2795 endif
2796 if(gr_flag.eq.1.and.gr_frac_roof.gt.0.)then
2797 do il=1,ngr_u
2798 qr_tmp(il)=qr(id,iz,il)
2799 enddo
2800 call flux_flat_roof(dz,z0v,ua_u(iz),va_u(iz),pt_u(iz),pt0_u(iz), &
2801 ptrv(id,iz),uhbv(id,iz), &
2802 vhbv(id,iz),sfrv(id,iz),lfrv(id,iz),ehbv(id,iz),da_u(iz),qv_u(iz),pr_u(iz),rs,qr_tmp,resg,rsveg,f1,f2,f3,f4,gr_type,pv_frac_roof)
2803 sfr(id,iz)=sfr(id,iz)+pv_frac_roof*sfrpv(id,iz)
2804 thb_u(id,iz)=-((1.-gr_frac_roof)*sfr(id,iz)+gr_frac_roof*sfrv(id,iz))/(da_u(iz)*cp_u)
2805 vhb_u(id,iz)=(1.-gr_frac_roof)*vhb(id,iz)+gr_frac_roof*vhbv(id,iz)
2806 uhb_u(id,iz)=(1.-gr_frac_roof)*uhb(id,iz)+gr_frac_roof*uhbv(id,iz)
2807 ehb_u(id,iz)=(1.-gr_frac_roof)*ehb(id,iz)+gr_frac_roof*ehbv(id,iz)
2808 qhb_u(id,iz)=-(gr_frac_roof*lfrv(id,iz)+(1.-gr_frac_roof)*lfr(id,iz))/(da_u(iz)*latent)
2809 sfr(id,iz)=sfr(id,iz)-pv_frac_roof*sfrpv(id,iz)
2810 else
2811 sfr(id,iz)=sfr(id,iz)+pv_frac_roof*sfrpv(id,iz)
2812 thb_u(id,iz)=-sfr(id,iz)/(da_u(iz)*cp_u)
2813 vhb_u(id,iz)=vhb(id,iz)
2814 uhb_u(id,iz)=uhb(id,iz)
2815 ehb_u(id,iz)=ehb(id,iz)
2816 qhb_u(id,iz)=-lfr(id,iz)/(da_u(iz)*latent)
2817 sfr(id,iz)=sfr(id,iz)-pv_frac_roof*sfrpv(id,iz)
2818 endif
2819 else
2820 uhb_u(id,iz) = 0.0
2821 vhb_u(id,iz) = 0.0
2822 thb_u(id,iz) = 0.0
2823 ehb_u(id,iz) = 0.0
2824 qhb_u(id,iz) = 0.0
2825 endif
2826 enddo
2830 ! Calculation at the wall surfaces
2832 do ibui=1,nbui
2833 do iz=1,nlev(ibui)
2835 call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
2836 ptw(2*id-1,iz,ibui),ptwin(2*id-1,iz,ibui), &
2837 uva_tmp,vva_tmp, &
2838 uvb_tmp,vvb_tmp, &
2839 sfw(2*id-1,iz,ibui),sfwin(2*id-1,iz,ibui), &
2840 evb_tmp,drst(id),dt,cdrag(iz))
2842 if (pb(iz+1).gt.0.) then
2844 uva_u(2*id-1,iz)=uva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
2845 vva_u(2*id-1,iz)=vva_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
2846 uvb_u(2*id-1,iz)=uvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
2847 vvb_u(2*id-1,iz)=vvb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
2848 evb_u(2*id-1,iz)=evb_u(2*id-1,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
2849 tvb_u(2*id-1,iz)=tvb_u(2*id-1,iz)-(d_urb(ibui)/pb(iz+1))* &
2850 (sfw(2*id-1,iz,ibui)*(1.-pwin)+sfwin(2*id-1,iz,ibui)*pwin)/ &
2851 da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
2852 (dz*bs_u(id,iurb))/da_u(iz)/cp_u
2853 tvb_ac(2*id-1,iz)=tvb_ac(2*id-1,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
2854 (dz*bs_u(id,iurb))/da_u(iz)/cp_u
2855 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))/&
2856 (dz*bs_u(id,iurb))/da_u(iz)/latent
2858 endif
2860 call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
2861 ptw(2*id,iz,ibui),ptwin(2*id,iz,ibui), &
2862 uva_tmp,vva_tmp, &
2863 uvb_tmp,vvb_tmp, &
2864 sfw(2*id,iz,ibui),sfwin(2*id,iz,ibui), &
2865 evb_tmp,drst(id),dt,cdrag(iz))
2867 if (pb(iz+1).gt.0.) then
2869 uva_u(2*id,iz)=uva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
2870 vva_u(2*id,iz)=vva_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
2871 uvb_u(2*id,iz)=uvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
2872 vvb_u(2*id,iz)=vvb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
2873 evb_u(2*id,iz)=evb_u(2*id,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
2874 tvb_u(2*id,iz)=tvb_u(2*id,iz)-(d_urb(ibui)/pb(iz+1))* &
2875 (sfw(2*id,iz,ibui)*(1.-pwin)+sfwin(2*id,iz,ibui)*pwin)/ &
2876 da_u(iz)/cp_u-(1./4.)*(d_urb(ibui)/pb(iz+1))*(sfvlev(iz,ibui)-sflev(iz,ibui))/&
2877 (dz*bs_u(id,iurb))/da_u(iz)/cp_u
2878 tvb_ac(2*id,iz)=tvb_ac(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(-sflev(iz,ibui))/&
2879 (dz*bs_u(id,iurb))/da_u(iz)/cp_u
2880 qvb_u(2*id,iz)=qvb_u(2*id,iz)-(1./4.)*(d_urb(ibui)/pb(iz+1))*(lfvlev(iz,ibui)-lflev(iz,ibui))/&
2881 (dz*bs_u(id,iurb))/da_u(iz)/latent
2883 endif
2884 !
2885 enddo !iz
2886 enddo !ibui
2888 end do !id
2890 return
2891 end subroutine buildings
2894 ! ===6=8===============================================================72
2895 ! ===6=8===============================================================72
2897 subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl, &
2898 uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
2899 uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u, &
2900 a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)
2902 ! ----------------------------------------------------------------------
2903 ! This routine interpolates the parameters from the "urban grid" to the
2904 ! "mesoscale grid".
2905 ! See p300-301 Appendix B.2 of the BLM paper.
2906 ! ----------------------------------------------------------------------
2908 implicit none
2910 ! ----------------------------------------------------------------------
2911 ! INPUT:
2912 ! ----------------------------------------------------------------------
2913 ! Data relative to the "mesoscale grid"
2914 integer kms,kme,kts,kte
2915 real z(kms:kme) ! Altitude above the ground of the cell interface
2916 real dz(kms:kme) ! Vertical space steps
2918 ! Data relative to the "uban grid"
2919 integer nz_u ! Number of layer in the urban grid
2920 integer nd ! Number of street direction for the current urban class
2921 real bs(ndm) ! Building widths of the current urban class
2922 real ws(ndm) ! Street widths of the current urban class
2923 real z_u(nz_um) ! Height of the urban grid levels
2924 real pb(nz_um) ! Probability to have a building with an height equal
2925 real ss(nz_um) ! Probability to have a building with height h
2926 real uhb_u(ndm,nz_um) ! U (x-wind component) Horizontal surfaces, B (explicit) term
2927 real uva_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, A (implicit) term
2928 real uvb_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, B (explicit) term
2929 real vhb_u(ndm,nz_um) ! V (y-wind component) Horizontal surfaces, B (explicit) term
2930 real vva_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, A (implicit) term
2931 real vvb_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, B (explicit) term
2932 real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
2933 real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
2934 real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
2935 real tvb_ac(2*ndm,nz_um)
2936 real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
2937 real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
2938 !
2939 !New variables for BEM
2940 !
2941 real qhb_u(ndm,nz_um)
2942 real qvb_u(2*ndm,nz_um)
2945 ! ----------------------------------------------------------------------
2946 ! OUTPUT:
2947 ! ----------------------------------------------------------------------
2948 ! Data relative to the "mesoscale grid"
2949 real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
2950 real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
2951 real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
2952 real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
2953 real a_t(kms:kme) ! Implicit component of the heat sources or sinks
2954 real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
2955 real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
2956 real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
2957 real b_t(kms:kme) ! Explicit component of the heat sources or sinks
2958 real b_ac(kms:kme)
2959 real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
2960 real b_q(kms:kme) ! Explicit component of the humidity sources or sinks
2961 ! ----------------------------------------------------------------------
2962 ! LOCAL:
2963 ! ----------------------------------------------------------------------
2964 real dzz
2965 real fact
2966 integer id,iz,iz_u
2967 real se,sr,st,su,sv,sq
2968 real uet(kms:kme) ! Contribution to TKE due to walls
2969 real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb,vqb,vtb_ac
2972 ! ----------------------------------------------------------------------
2973 ! END VARIABLES DEFINITIONS
2974 ! ----------------------------------------------------------------------
2976 ! initialisation
2978 do iz=kts,kte
2979 a_u(iz)=0.
2980 a_v(iz)=0.
2981 a_t(iz)=0.
2982 a_e(iz)=0.
2983 b_u(iz)=0.
2984 b_v(iz)=0.
2985 b_e(iz)=0.
2986 b_t(iz)=0.
2987 b_ac(iz)=0.
2988 uet(iz)=0.
2989 b_q(iz)=0.
2990 end do
2992 ! horizontal surfaces
2993 do iz=kts,kte
2994 sf(iz)=0.
2995 vl(iz)=0.
2996 enddo
2997 sf(kte+1)=0.
2999 do id=1,nd
3000 do iz=kts+1,kte+1
3001 sr=0.
3002 do iz_u=2,nz_u
3003 if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
3004 sr=pb(iz_u)
3005 endif
3006 enddo
3007 sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
3008 enddo
3009 enddo
3011 ! volume
3012 do iz=kts,kte
3013 do id=1,nd
3014 vtot=0.
3015 do iz_u=1,nz_u
3016 dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
3017 vtot=vtot+pb(iz_u+1)*dzz
3018 enddo
3019 vtot=vtot/(z(iz+1)-z(iz))
3020 vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
3021 enddo
3022 enddo
3024 ! horizontal surface impact
3026 do id=1,nd
3028 fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
3029 b_t(kts)=b_t(kts)+thb_u(id,1)*fact
3030 b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
3031 b_v(kts)=b_v(kts)+vhb_u(id,1)*fact
3032 b_e(kts)=b_e(kts)+ehb_u(id,1)*fact*(z_u(2)-z_u(1))
3033 b_q(kts)=b_q(kts)+qhb_u(id,1)*fact
3035 do iz=kts,kte
3036 st=0.
3037 su=0.
3038 sv=0.
3039 se=0.
3040 sq=0.
3041 do iz_u=2,nz_u
3042 if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
3043 st=st+ss(iz_u)*thb_u(id,iz_u)
3044 su=su+ss(iz_u)*uhb_u(id,iz_u)
3045 sv=sv+ss(iz_u)*vhb_u(id,iz_u)
3046 se=se+ss(iz_u)*ehb_u(id,iz_u)*(z_u(iz_u+1)-z_u(iz_u))
3047 sq=sq+ss(iz_u)*qhb_u(id,iz_u)
3048 endif
3049 enddo
3051 fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
3052 b_t(iz)=b_t(iz)+st*fact
3053 b_u(iz)=b_u(iz)+su*fact
3054 b_v(iz)=b_v(iz)+sv*fact
3055 b_e(iz)=b_e(iz)+se*fact
3056 b_q(iz)=b_q(iz)+sq*fact
3057 enddo
3058 enddo
3060 ! vertical surface impact
3062 do iz=kts,kte
3063 uet(iz)=0.
3064 do id=1,nd
3065 vtb=0.
3066 vtb_ac=0.
3067 vta=0.
3068 vua=0.
3069 vub=0.
3070 vva=0.
3071 vvb=0.
3072 veb=0.
3073 vte=0.
3074 vqb=0.
3075 do iz_u=1,nz_u
3076 dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
3077 fact=dzz/(ws(id)+bs(id))
3078 vtb=vtb+pb(iz_u+1)* &
3079 (tvb_u(2*id-1,iz_u)+tvb_u(2*id,iz_u))*fact
3080 vtb_ac=vtb_ac+pb(iz_u+1)* &
3081 (tvb_ac(2*id-1,iz_u)+tvb_ac(2*id,iz_u))*fact
3082 vta=vta+pb(iz_u+1)* &
3083 (tva_u(2*id-1,iz_u)+tva_u(2*id,iz_u))*fact
3084 vua=vua+pb(iz_u+1)* &
3085 (uva_u(2*id-1,iz_u)+uva_u(2*id,iz_u))*fact
3086 vva=vva+pb(iz_u+1)* &
3087 (vva_u(2*id-1,iz_u)+vva_u(2*id,iz_u))*fact
3088 vub=vub+pb(iz_u+1)* &
3089 (uvb_u(2*id-1,iz_u)+uvb_u(2*id,iz_u))*fact
3090 vvb=vvb+pb(iz_u+1)* &
3091 (vvb_u(2*id-1,iz_u)+vvb_u(2*id,iz_u))*fact
3092 veb=veb+pb(iz_u+1)* &
3093 (evb_u(2*id-1,iz_u)+evb_u(2*id,iz_u))*fact
3094 vqb=vqb+pb(iz_u+1)* &
3095 (qvb_u(2*id-1,iz_u)+qvb_u(2*id,iz_u))*fact
3096 enddo
3098 fact=1./vl(iz)/dz(iz)/nd
3099 b_t(iz)=b_t(iz)+vtb*fact
3100 b_ac(iz)=b_ac(iz)+vtb_ac*fact
3101 a_t(iz)=a_t(iz)+vta*fact
3102 a_u(iz)=a_u(iz)+vua*fact
3103 a_v(iz)=a_v(iz)+vva*fact
3104 b_u(iz)=b_u(iz)+vub*fact
3105 b_v(iz)=b_v(iz)+vvb*fact
3106 b_e(iz)=b_e(iz)+veb*fact
3107 uet(iz)=uet(iz)+vte*fact
3108 b_q(iz)=b_q(iz)+vqb*fact
3109 enddo
3110 enddo
3114 return
3115 end subroutine urban_meso
3117 ! ===6=8===============================================================72
3118 ! ===6=8===============================================================72
3120 subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
3121 dlg,dl_u)
3123 ! ----------------------------------------------------------------------
3124 ! Calculation of the length scales
3125 ! See p272-274 formula (22) and (24) of the BLM paper
3126 ! ----------------------------------------------------------------------
3128 implicit none
3131 ! ----------------------------------------------------------------------
3132 ! INPUT:
3133 ! ----------------------------------------------------------------------
3134 integer kms,kme,kts,kte
3135 real z(kms:kme) ! Altitude above the ground of the cell interface
3136 integer nd ! Number of street direction for the current urban class
3137 integer nz_u ! Number of levels in the "urban grid"
3138 real z_u(nz_um) ! Height of the urban grid levels
3139 real bs(ndm) ! Building widths of the current urban class
3140 real ss(nz_um) ! Probability to have a building with height h
3141 real ws(ndm) ! Street widths of the current urban class
3144 ! ----------------------------------------------------------------------
3145 ! OUTPUT:
3146 ! ----------------------------------------------------------------------
3147 real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
3148 real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
3150 ! ----------------------------------------------------------------------
3151 ! LOCAL:
3152 ! ----------------------------------------------------------------------
3153 real dlgtmp
3154 integer id,iz,iz_u
3155 real sftot
3156 real ulu,ssl
3158 ! ----------------------------------------------------------------------
3159 ! END VARIABLES DEFINITIONS
3160 ! ----------------------------------------------------------------------
3162 do iz=kts,kte
3163 ulu=0.
3164 ssl=0.
3165 do id=1,nd
3166 do iz_u=2,nz_u
3167 if(z_u(iz_u).gt.z(iz))then
3168 ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
3169 ssl=ssl+ss(iz_u)/nd
3170 endif
3171 enddo
3172 enddo
3174 if(ulu.ne.0)then
3175 dl_u(iz)=ssl/ulu
3176 else
3177 dl_u(iz)=0.
3178 endif
3179 enddo
3182 do iz=kts,kte
3183 dlg(iz)=0.
3184 do id=1,nd
3185 sftot=ws(id)
3186 dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
3187 do iz_u=1,nz_u
3188 if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
3189 dlgtmp=dlgtmp+ss(iz_u)*bs(id)/ &
3190 ((z(iz)+z(iz+1))/2.-z_u(iz_u))
3191 sftot=sftot+ss(iz_u)*bs(id)
3192 endif
3193 enddo
3194 dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
3195 enddo
3196 dlg(iz)=1./dlg(iz)
3197 enddo
3199 return
3200 end subroutine interp_length
3202 ! ===6=8===============================================================72
3203 ! ===6=8===============================================================72
3205 subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
3206 swddir,rsw,rsg,xlat)
3208 ! ----------------------------------------------------------------------
3209 ! Modification of short wave radiation to take into account
3210 ! the shadow produced by the buildings
3211 ! ----------------------------------------------------------------------
3213 implicit none
3215 ! ----------------------------------------------------------------------
3216 ! INPUT:
3217 ! ----------------------------------------------------------------------
3218 integer nd ! Number of street direction for the current urban class
3219 integer nz_u ! number of vertical layers defined in the urban grid
3220 real ah ! Hour angle (it should come from the radiation routine)
3221 real deltar ! Declination of the sun
3222 real drst(ndm) ! street directions for the current urban class
3223 real swddir ! solar radiation
3224 real ss(nz_um) ! probability to have a building with height h
3225 real pb(nz_um) ! Probability that a building has an height greater or equal to h
3226 real ws(ndm) ! Street width of the current urban class
3227 real z(nz_um) ! Height of the urban grid levels
3228 real zr ! zenith angle
3229 real xlat
3230 real xlat_r
3231 ! ----------------------------------------------------------------------
3232 ! OUTPUT:
3233 ! ----------------------------------------------------------------------
3234 real rsg(ndm) ! Short wave radiation at the ground
3235 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
3237 ! ----------------------------------------------------------------------
3238 ! LOCAL:
3239 ! ----------------------------------------------------------------------
3240 integer id,iz,jz
3241 real aae,aaw,bbb,phix,rd,rtot,wsd
3243 ! ----------------------------------------------------------------------
3244 ! END VARIABLES DEFINITIONS
3245 ! ----------------------------------------------------------------------
3247 xlat_r=xlat*pi/180
3249 if(swddir.eq.0.or.sin(zr).eq.1)then
3250 do id=1,nd
3251 rsg(id)=0.
3252 do iz=1,nz_u
3253 rsw(2*id-1,iz)=0.
3254 rsw(2*id,iz)=0.
3255 enddo
3256 enddo
3257 else
3258 !test
3260 if(abs(sin(zr)).gt.1.e-10)then
3261 if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
3262 bbb=pi/2.
3263 elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
3264 bbb=-pi/2.
3265 else
3266 bbb=asin(cos(deltar)*sin(ah)/sin(zr)) !
3267 if(sin(deltar).lt.(cos(zr)*sin(xlat_r)))then !
3268 bbb=pi-bbb !
3269 endif
3270 endif
3271 else
3272 if(cos(deltar)*sin(ah).ge.0)then
3273 bbb=pi/2.
3274 elseif(cos(deltar)*sin(ah).lt.0)then
3275 bbb=-pi/2.
3276 endif
3277 endif
3278 phix=zr
3280 do id=1,nd
3282 rsg(id)=0.
3284 aae=bbb-drst(id)
3285 aaw=bbb-drst(id)+pi
3287 do iz=1,nz_u
3288 rsw(2*id-1,iz)=0.
3289 rsw(2*id,iz)=0.
3290 if(pb(iz+1).gt.0.)then
3291 do jz=1,nz_u
3292 if(abs(sin(aae)).gt.1.e-10)then
3293 call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae, &
3294 ws(id),rd)
3295 rsw(2*id-1,iz)=rsw(2*id-1,iz)+swddir*rd*ss(jz+1)/pb(iz+1)
3296 endif
3298 if(abs(sin(aaw)).gt.1.e-10)then
3299 call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw, &
3300 ws(id),rd)
3301 rsw(2*id,iz)=rsw(2*id,iz)+swddir*rd*ss(jz+1)/pb(iz+1)
3302 endif
3303 enddo
3304 endif
3305 enddo
3306 if(abs(sin(aae)).gt.1.e-10)then
3307 wsd=abs(ws(id)/sin(aae))
3309 do jz=1,nz_u
3310 rd=max(0.,wsd-z(jz+1)*tan(phix))
3311 rsg(id)=rsg(id)+swddir*rd*ss(jz+1)/wsd
3312 enddo
3313 rtot=0.
3315 do iz=1,nz_u
3316 rtot=rtot+(rsw(2*id,iz)+rsw(2*id-1,iz))* &
3317 (z(iz+1)-z(iz))
3318 enddo
3319 rtot=rtot+rsg(id)*ws(id)
3320 else
3321 rsg(id)=swddir
3322 endif
3326 enddo
3327 endif
3329 return
3330 end subroutine shadow_mas
3336 ! ===6=8===============================================================72
3337 ! ===6=8===============================================================72
3339 subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
3341 ! ----------------------------------------------------------------------
3342 ! This routine computes the effects of a shadow induced by a building of
3343 ! height hu, on a portion of wall between z1 and z2. See equation A10,
3344 ! and correction described below formula A11, and figure A1. Basically rd
3345 ! is the ratio between the horizontal surface illuminated and the portion
3346 ! of wall. Referring to figure A1, multiplying radiation flux density on
3347 ! a horizontal surface (rs) by x1-x2 we have the radiation energy per
3348 ! unit time. Dividing this by z2-z1, we obtain the radiation flux
3349 ! density reaching the portion of the wall between z2 and z1
3350 ! (everything is assumed in 2D)
3351 ! ----------------------------------------------------------------------
3353 implicit none
3355 ! ----------------------------------------------------------------------
3356 ! INPUT:
3357 ! ----------------------------------------------------------------------
3358 real aa ! Angle between the sun direction and the face of the wall (A12)
3359 real hu ! Height of the building that generates the shadow
3360 real phix ! Solar zenith angle
3361 real ws ! Width of the street
3362 real z1 ! Height of the level z(iz)
3363 real z2 ! Height of the level z(iz+1)
3365 ! ----------------------------------------------------------------------
3366 ! OUTPUT:
3367 ! ----------------------------------------------------------------------
3368 real rd ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A.
3369 ! Multiplying rd by rs (radiation flux
3370 ! density on a horizontal surface) gives
3371 ! the radiation flux density on the
3372 ! portion of wall between z1 and z2.
3373 ! ----------------------------------------------------------------------
3374 ! LOCAL:
3375 ! ----------------------------------------------------------------------
3376 real x1,x2 ! x1,x2 see Fig. A1.
3378 ! ----------------------------------------------------------------------
3379 ! END VARIABLES DEFINITIONS
3380 ! ----------------------------------------------------------------------
3382 x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
3384 x2=max((hu-z2)*tan(phix),0.)
3386 rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
3388 return
3389 end subroutine shade_wall
3391 ! ===6=8===============================================================72
3392 ! ===6=8===============================================================72
3394 subroutine long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev,&
3395 fwg,fww,fgw,fsw,fsg,tg_av,tw,rlg,rlw,rl,pb)
3397 ! ----------------------------------------------------------------------
3398 ! This routine computes the effects of the reflections of long-wave
3399 ! radiation in the street canyon by solving the system
3400 ! of 2*nz_u+1 eqn. in 2*nz_u+1
3401 ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
3402 ! The system is solved by solving A X= B,
3403 ! with A matrix B vector, and X solution.
3404 ! ----------------------------------------------------------------------
3406 implicit none
3410 ! ----------------------------------------------------------------------
3411 ! INPUT:
3412 ! ----------------------------------------------------------------------
3413 real emg ! Emissivity of ground for the current urban class
3414 real emw ! Emissivity of wall for the current urban class
3415 real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
3416 real fsg(ndm,nurbm) ! View factors from sky to ground
3417 real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
3418 real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
3419 real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
3420 integer id ! Current street direction
3421 integer iurb ! Current urban class
3422 integer nz_u ! Number of layer in the urban grid
3423 real pb(nz_um) ! Probability to have a building with an height equal
3424 real rl ! Downward flux of the longwave radiation
3425 real tg_av(ndm) ! Temperature in each layer of the ground [K]
3426 real tw(2*ndm,nz_um) ! Temperature in each layer of the wall [K]
3427 !
3428 !New Variables for BEM
3429 !
3430 real twlev(2*ndm,nz_um) ! Window temperature in BEM [K]
3431 real emwin ! Emissivity of windows
3432 real pwin ! Coverage area fraction of windows in the walls of the buildings (BEM)
3435 ! ----------------------------------------------------------------------
3436 ! OUTPUT:
3437 ! ----------------------------------------------------------------------
3438 real rlg(ndm) ! Long wave radiation at the ground
3439 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
3441 ! ----------------------------------------------------------------------
3442 ! LOCAL:
3443 ! ----------------------------------------------------------------------
3444 integer i,j
3445 real aaa(2*nz_um+1,2*nz_um+1) ! terms of the matrix
3446 real bbb(2*nz_um+1) ! terms of the vector
3448 ! ----------------------------------------------------------------------
3449 ! END VARIABLES DEFINITIONS
3450 ! ----------------------------------------------------------------------
3453 ! west wall
3455 do i=1,nz_u
3457 do j=1,nz_u
3458 aaa(i,j)=0.
3459 enddo
3461 aaa(i,i)=1.
3463 do j=nz_u+1,2*nz_u
3464 aaa(i,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
3465 fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
3466 enddo
3468 !! aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)*pb(i+1)
3469 aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i,id,iurb)
3471 bbb(i)=fsw(i,id,iurb)*rl+emg*fgw(i,id,iurb)*sigma*tg_av(id)**4
3472 do j=1,nz_u
3473 bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i,id,iurb)* &
3474 (emw*(1.-pwin)*tw(2*id,j)**4+emwin*pwin*twlev(2*id,j)**4)+ &
3475 fww(j,i,id,iurb)*rl*(1.-pb(j+1))
3476 enddo
3478 enddo
3480 ! east wall
3482 do i=1+nz_u,2*nz_u
3484 do j=1,nz_u
3485 aaa(i,j)=-(1.-emw*(1.-pwin)-emwin*pwin)*fww(j,i-nz_u,id,iurb)*pb(j+1)
3486 enddo
3488 do j=1+nz_u,2*nz_u
3489 aaa(i,j)=0.
3490 enddo
3492 aaa(i,i)=1.
3494 !! aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)*pb(i-nz_u+1)
3495 aaa(i,2*nz_u+1)=-(1.-emg)*fgw(i-nz_u,id,iurb)
3497 bbb(i)=fsw(i-nz_u,id,iurb)*rl+ &
3498 emg*fgw(i-nz_u,id,iurb)*sigma*tg_av(id)**4
3500 do j=1,nz_u
3501 bbb(i)=bbb(i)+pb(j+1)*sigma*fww(j,i-nz_u,id,iurb)* &
3502 (emw*(1.-pwin)*tw(2*id-1,j)**4+emwin*pwin*twlev(2*id-1,j)**4)+&
3503 fww(j,i-nz_u,id,iurb)*rl*(1.-pb(j+1))
3504 enddo
3506 enddo
3508 ! ground
3509 do j=1,nz_u
3510 aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
3511 fwg(j,id,iurb)*pb(j+1)
3512 enddo
3514 do j=nz_u+1,2*nz_u
3515 aaa(2*nz_u+1,j)=-(1.-emw*(1.-pwin)-emwin*pwin)* &
3516 fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
3517 enddo
3519 aaa(2*nz_u+1,2*nz_u+1)=1.
3521 bbb(2*nz_u+1)=fsg(id,iurb)*rl
3523 do i=1,nz_u
3524 bbb(2*nz_u+1)=bbb(2*nz_u+1)+sigma*fwg(i,id,iurb)*pb(i+1)* &
3525 (emw*(1.-pwin)*(tw(2*id-1,i)**4+tw(2*id,i)**4)+ &
3526 emwin*pwin*(twlev(2*id-1,i)**4+twlev(2*id,i)**4))+ &
3527 2.*fwg(i,id,iurb)*(1.-pb(i+1))*rl
3528 enddo
3532 call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
3534 do i=1,nz_u
3535 rlw(2*id-1,i)=bbb(i)
3536 enddo
3538 do i=nz_u+1,2*nz_u
3539 rlw(2*id,i-nz_u)=bbb(i)
3540 enddo
3542 rlg(id)=bbb(2*nz_u+1)
3544 return
3545 end subroutine long_rad
3547 ! ===6=8===============================================================72
3548 ! ===6=8===============================================================72
3551 subroutine short_rad_dd(iurb,nz_u,id,albw, &
3552 albg,rsdif,fwg,fww,fgw,fsw,fsg,rsg,rsw,pb)
3554 ! ----------------------------------------------------------------------
3555 ! This routine computes the effects of the reflections of short-wave
3556 ! (solar) radiation in the street canyon by solving the system
3557 ! of 2*nz_u+1 eqn. in 2*nz_u+1
3558 ! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
3559 ! The system is solved by solving A X= B,
3560 ! with A matrix B vector, and X solution.
3561 ! ----------------------------------------------------------------------
3563 implicit none
3567 ! ----------------------------------------------------------------------
3568 ! INPUT:
3569 ! ----------------------------------------------------------------------
3570 real albg ! Albedo of the ground for the current urban class
3571 real albw ! Albedo of the wall for the current urban class
3572 real rsdif ! diffused short wave radiation
3573 real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
3574 real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
3575 real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
3576 real fsg(ndm,nurbm) ! View factors from sky to ground
3577 real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
3578 integer id ! current street direction
3579 integer iurb ! current urban class
3580 integer nz_u ! Number of layer in the urban grid
3581 real pb(nz_um) ! probability to have a building with an height equal
3583 ! ----------------------------------------------------------------------
3584 ! OUTPUT:
3585 ! ----------------------------------------------------------------------
3586 real rsg(ndm) ! Short wave radiation at the ground
3587 real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
3589 ! ----------------------------------------------------------------------
3590 ! LOCAL:
3591 ! ----------------------------------------------------------------------
3592 integer i,j
3593 real aaa(2*nz_um+1,2*nz_um+1) ! terms of the matrix
3594 real bbb(2*nz_um+1) ! terms of the vector
3596 ! ----------------------------------------------------------------------
3597 ! END VARIABLES DEFINITIONS
3598 ! ----------------------------------------------------------------------
3601 ! west wall
3604 do i=1,nz_u
3605 do j=1,nz_u
3606 aaa(i,j)=0.
3607 enddo
3609 aaa(i,i)=1.
3611 do j=nz_u+1,2*nz_u
3612 aaa(i,j)=-albw*fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
3613 enddo
3615 aaa(i,2*nz_u+1)=-albg*fgw(i,id,iurb)
3616 bbb(i)=rsw(2*id-1,i)+fsw(i,id,iurb)*rsdif
3618 enddo
3620 ! east wall
3621 do i=1+nz_u,2*nz_u
3622 do j=1,nz_u
3623 aaa(i,j)=-albw*fww(j,i-nz_u,id,iurb)*pb(j+1)
3624 enddo
3626 do j=1+nz_u,2*nz_u
3627 aaa(i,j)=0.
3628 enddo
3630 aaa(i,i)=1.
3631 aaa(i,2*nz_u+1)=-albg*fgw(i-nz_u,id,iurb)
3632 bbb(i)=rsw(2*id,i-nz_u)+fsw(i-nz_u,id,iurb)*rsdif
3634 enddo
3637 ! ground
3638 do j=1,nz_u
3639 aaa(2*nz_u+1,j)=-albw*fwg(j,id,iurb)*pb(j+1)
3640 enddo
3642 do j=nz_u+1,2*nz_u
3643 aaa(2*nz_u+1,j)=-albw*fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
3644 enddo
3646 aaa(2*nz_u+1,2*nz_u+1)=1.
3647 bbb(2*nz_u+1)=rsg(id)+fsg(id,iurb)*rsdif
3649 call gaussj(aaa,2*nz_u+1,bbb,2*nz_um+1)
3651 do i=1,nz_u
3652 rsw(2*id-1,i)=bbb(i)
3653 enddo
3655 do i=nz_u+1,2*nz_u
3656 rsw(2*id,i-nz_u)=bbb(i)
3657 enddo
3659 rsg(id)=bbb(2*nz_u+1)
3662 return
3663 end subroutine short_rad_dd
3667 ! ===6=8===============================================================72
3668 ! ===6=8===============================================================72
3670 subroutine gaussj(a,n,b,np)
3672 ! ----------------------------------------------------------------------
3673 ! This routine solve a linear system of n equations of the form
3674 ! A X = B
3675 ! where A is a matrix a(i,j)
3676 ! B a vector and X the solution
3677 ! In output b is replaced by the solution
3678 ! ----------------------------------------------------------------------
3680 implicit none
3682 ! ----------------------------------------------------------------------
3683 ! INPUT:
3684 ! ----------------------------------------------------------------------
3685 integer np
3686 real a(np,np)
3688 ! ----------------------------------------------------------------------
3689 ! OUTPUT:
3690 ! ----------------------------------------------------------------------
3691 real b(np)
3693 ! ----------------------------------------------------------------------
3694 ! LOCAL:
3695 ! ----------------------------------------------------------------------
3696 integer nmax
3697 parameter (nmax=150)
3699 real big,dum
3700 integer i,icol,irow
3701 integer j,k,l,ll,n
3702 integer ipiv(nmax)
3703 real pivinv
3705 ! ----------------------------------------------------------------------
3706 ! END VARIABLES DEFINITIONS
3707 ! ----------------------------------------------------------------------
3709 do j=1,n
3710 ipiv(j)=0.
3711 enddo
3713 do i=1,n
3714 big=0.
3715 do j=1,n
3716 if(ipiv(j).ne.1)then
3717 do k=1,n
3718 if(ipiv(k).eq.0)then
3719 if(abs(a(j,k)).ge.big)then
3720 big=abs(a(j,k))
3721 irow=j
3722 icol=k
3723 endif
3724 elseif(ipiv(k).gt.1)then
3725 FATAL_ERROR('singular matrix in gaussj')
3726 endif
3727 enddo
3728 endif
3729 enddo
3731 ipiv(icol)=ipiv(icol)+1
3733 if(irow.ne.icol)then
3734 do l=1,n
3735 dum=a(irow,l)
3736 a(irow,l)=a(icol,l)
3737 a(icol,l)=dum
3738 enddo
3740 dum=b(irow)
3741 b(irow)=b(icol)
3742 b(icol)=dum
3744 endif
3746 if(a(icol,icol).eq.0) FATAL_ERROR('singular matrix in gaussj')
3748 pivinv=1./a(icol,icol)
3749 a(icol,icol)=1
3751 do l=1,n
3752 a(icol,l)=a(icol,l)*pivinv
3753 enddo
3755 b(icol)=b(icol)*pivinv
3757 do ll=1,n
3758 if(ll.ne.icol)then
3759 dum=a(ll,icol)
3760 a(ll,icol)=0.
3761 do l=1,n
3762 a(ll,l)=a(ll,l)-a(icol,l)*dum
3763 enddo
3765 b(ll)=b(ll)-b(icol)*dum
3767 endif
3768 enddo
3769 enddo
3771 return
3772 end subroutine gaussj
3776 ! ===6=8===============================================================72
3777 ! ===6=8===============================================================72
3779 subroutine soil_moist(nz,dz,qv,dt,lf,d,k,rainbl,drain,irri_now)
3781 ! ----------------------------------------------------------------------
3782 ! This routine solves the Fourier diffusion equation for heat in
3783 ! the material (wall, roof, or ground). Resolution is done implicitely.
3784 ! Boundary conditions are:
3785 ! - fixed temperature at the interior
3786 ! - energy budget at the surface
3787 ! ----------------------------------------------------------------------
3789 implicit none
3793 ! ----------------------------------------------------------------------
3794 ! INPUT:
3795 ! ----------------------------------------------------------------------
3796 integer nz ! Number of layers
3797 real dt ! Time step
3798 real lf ! Latent heat flux at the surface
3799 real qv(nz) ! Moisture in each layer [K]
3800 real dz(nz) ! Layer sizes [m]
3801 real rainbl ! Rainfall [mm]
3802 real d(nz) ! Soil water diffusivity
3803 real k(nz) ! Hydraulic conductivity
3804 real gr ! Dummy variable
3805 real drain
3806 real irri_now
3807 ! ----------------------------------------------------------------------
3808 ! OUTPUT:
3809 ! ----------------------------------------------------------------------
3812 ! ----------------------------------------------------------------------
3813 ! LOCAL:
3814 ! ----------------------------------------------------------------------
3815 integer iz
3816 real a(nz,3)
3817 real alpha
3818 real c(nz)
3819 real cddz(nz+2)
3820 real dw !water density Kg/m3
3821 parameter(dw=1000.)
3822 !----------------------------------------------------------------------
3823 ! END VARIABLES DEFINITIONS
3824 ! ----------------------------------------------------------------------
3826 alpha=rainbl/(dw*dt)+lf/latent/dw+irri_now/dw
3827 cddz(1)=0.
3828 do iz=2,nz
3829 cddz(iz)=2.*d(iz)/(dz(iz)+dz(iz-1))
3830 enddo
3831 do iz=1,4
3832 a(iz,1)=0.
3833 a(iz,2)=1.
3834 a(iz,3)=0.
3835 c(iz)=qv(iz)
3836 enddo
3837 do iz=6,nz-1
3838 a(iz,1)=-cddz(iz)*dt/dz(iz)
3839 a(iz,2)=1.+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
3840 a(iz,3)=-cddz(iz+1)*dt/dz(iz)
3841 c(iz)=qv(iz)+dt*(k(iz+1)-k(iz))/dz(iz)
3842 enddo
3843 a(5,1)=0.
3844 a(5,2)=1.+dt*(cddz(5+1))/dz(5)
3845 a(5,3)=-cddz(5+1)*dt/dz(5)
3846 c(5)=qv(5)+dt*(k(5+1)-drain)/dz(5)
3849 a(nz,1)=-dt*cddz(nz)/dz(nz)
3850 a(nz,2)=1.+dt*cddz(nz)/dz(nz)
3851 a(nz,3)=0.
3852 c(nz)=qv(nz)+dt*alpha/dz(nz)-dt*k(nz-1)/dz(nz)
3854 call invert(nz,a,c,qv)
3856 return
3857 end subroutine soil_moist
3859 ! ===6=8===============================================================72
3860 ! ===6=8===============================================================72
3863 ! ===6=8===============================================================72
3864 ! ===6=8===============================================================72
3866 subroutine soil_temp_veg(heflro,nz,dz,temp,pt,ala,cs, &
3867 rs,rl,press,dt,em,alb,rt,sf,lf,gf,pv_frac_roof,tpv)
3869 ! ----------------------------------------------------------------------
3870 ! This routine solves the Fourier diffusion equation for heat in
3871 ! the material (wall, roof, or ground). Resolution is done implicitely.
3872 ! Boundary conditions are:
3873 ! - fixed temperature at the interior
3874 ! - energy budget at the surface
3875 ! ----------------------------------------------------------------------
3877 implicit none
3881 ! ----------------------------------------------------------------------
3882 ! INPUT:
3883 ! ----------------------------------------------------------------------
3884 integer nz ! Number of layers
3885 real ala(nz) ! Thermal diffusivity in each layers [m^2 s^-1]
3886 real alb ! Albedo of the surface
3887 real cs(nz) ! Specific heat of the material [J m^3 K^-1]
3888 real dt ! Time step
3889 real em ! Emissivity of the surface
3890 real press ! Pressure at ground level
3891 real rl ! Downward flux of the longwave radiation
3892 real rs ! Solar radiation
3893 real sf ! Sensible heat flux at the surface
3894 real lf ! Latent heat flux at the surface
3895 real temp(nz) ! Temperature in each layer [K]
3896 real dz(nz) ! Layer sizes [m]
3897 real heflro ! Heat flux between roof and green roof
3898 real rs_eff
3899 real rl_eff
3900 real tpv
3901 real pv_frac_roof
3902 ! ----------------------------------------------------------------------
3903 ! OUTPUT:
3904 ! ----------------------------------------------------------------------
3905 real gf ! Heat flux transferred from the surface toward the interior
3906 real pt ! Potential temperature at the surface
3907 real rt ! Total radiation at the surface (solar+incoming long+outgoing long)
3909 ! ----------------------------------------------------------------------
3910 ! LOCAL:
3911 ! ----------------------------------------------------------------------
3912 integer iz
3913 real a(nz,3)
3914 real alpha
3915 real c(nz)
3916 real cddz(nz+2)
3917 real tsig
3919 ! ----------------------------------------------------------------------
3920 ! END VARIABLES DEFINITIONS
3921 ! ----------------------------------------------------------------------
3922 if(pv_frac_roof.gt.0)then
3923 rl_eff=(1-pv_frac_roof)*em*rl+em*sigma*tpv**4*pv_frac_roof
3924 rs_eff=(1.-pv_frac_roof)*rs
3925 else
3926 rl_eff=em*rl
3927 rs_eff=rs
3928 endif
3929 tsig=temp(nz)
3930 alpha=(1.-alb)*rs_eff+rl_eff-em*sigma*(tsig**4)+sf+lf
3931 cddz(1)=ala(1)/dz(1)
3932 do iz=2,nz
3933 cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
3934 enddo
3936 a(1,1)=0.
3937 a(1,2)=1.
3938 a(1,3)=0.
3939 c(1)=temp(1)-heflro*dt/dz(1)
3940 do iz=2,nz-1
3941 a(iz,1)=-cddz(iz)*dt/dz(iz)
3942 a(iz,2)=1.+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
3943 a(iz,3)=-cddz(iz+1)*dt/dz(iz)
3944 c(iz)=temp(iz)
3945 enddo
3946 a(nz,1)=-dt*cddz(nz)/dz(nz)
3947 a(nz,2)=1.+dt*cddz(nz)/dz(nz)
3948 a(nz,3)=0.
3949 c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
3951 call invert(nz,a,c,temp)
3953 pt=temp(nz)*(press/1.e+5)**(-rcp_u)
3955 rt=(1.-alb)*rs_eff+rl_eff-em*sigma*(tsig**4.)
3957 gf=(1.-alb)*rs_eff+rl_eff-em*sigma*(tsig**4.)+sf
3958 return
3959 end subroutine soil_temp_veg
3961 ! ===6=8===============================================================72
3962 ! ===6=8===============================================================72
3964 subroutine soil_temp(nz,dz,temp,pt,ala,cs, &
3965 rs,rl,press,dt,em,alb,rt,sf,lf,gf)
3967 ! ----------------------------------------------------------------------
3968 ! This routine solves the Fourier diffusion equation for heat in
3969 ! the material (wall, roof, or ground). Resolution is done implicitely.
3970 ! Boundary conditions are:
3971 ! - fixed temperature at the interior
3972 ! - energy budget at the surface
3973 ! ----------------------------------------------------------------------
3975 implicit none
3979 ! ----------------------------------------------------------------------
3980 ! INPUT:
3981 ! ----------------------------------------------------------------------
3982 integer nz ! Number of layers
3983 real ala(nz) ! Thermal diffusivity in each layers [m^2 s^-1]
3984 real alb ! Albedo of the surface
3985 real cs(nz) ! Specific heat of the material [J m^3 K^-1]
3986 real dt ! Time step
3987 real em ! Emissivity of the surface
3988 real press ! Pressure at ground level
3989 real rl ! Downward flux of the longwave radiation
3990 real rs ! Solar radiation
3991 real sf ! Sensible heat flux at the surface
3992 real lf ! Latent heat flux at the surface
3993 real temp(nz) ! Temperature in each layer [K]
3994 real dz(nz) ! Layer sizes [m]
3996 ! ----------------------------------------------------------------------
3997 ! OUTPUT:
3998 ! ----------------------------------------------------------------------
3999 real gf ! Heat flux transferred from the surface toward the interior
4000 real pt ! Potential temperature at the surface
4001 real rt ! Total radiation at the surface (solar+incoming long+outgoing long)
4003 ! ----------------------------------------------------------------------
4004 ! LOCAL:
4005 ! ----------------------------------------------------------------------
4006 integer iz
4007 real a(nz,3)
4008 real alpha
4009 real c(nz)
4010 real cddz(nz+2)
4011 real tsig
4013 ! ----------------------------------------------------------------------
4014 ! END VARIABLES DEFINITIONS
4015 ! ----------------------------------------------------------------------
4017 tsig=temp(nz)
4018 alpha=(1.-alb)*rs+em*rl-em*sigma*(tsig**4)+sf+lf
4019 ! Compute cddz=2*cd/dz
4020 cddz(1)=ala(1)/dz(1)
4021 do iz=2,nz
4022 cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
4023 enddo
4025 a(1,1)=0.
4026 a(1,2)=1.
4027 a(1,3)=0.
4028 c(1)=temp(1)
4029 do iz=2,nz-1
4030 a(iz,1)=-cddz(iz)*dt/dz(iz)
4031 a(iz,2)=1.+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
4032 a(iz,3)=-cddz(iz+1)*dt/dz(iz)
4033 c(iz)=temp(iz)
4034 enddo
4035 a(nz,1)=-dt*cddz(nz)/dz(nz)
4036 a(nz,2)=1.+dt*cddz(nz)/dz(nz)
4037 a(nz,3)=0.
4038 c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
4041 call invert(nz,a,c,temp)
4043 pt=temp(nz)*(press/1.e+5)**(-rcp_u)
4045 rt=(1.-alb)*rs+em*rl-em*sigma*(tsig**4.)
4047 gf=(1.-alb)*rs+em*rl-em*sigma*(tsig**4.)+sf
4048 return
4049 end subroutine soil_temp
4052 ! ===6=8===============================================================72
4053 ! ===6=8===============================================================72
4056 subroutine invert(n,a,c,x)
4058 ! ----------------------------------------------------------------------
4059 ! Inversion and resolution of a tridiagonal matrix
4060 ! A X = C
4061 ! ----------------------------------------------------------------------
4063 implicit none
4065 ! ----------------------------------------------------------------------
4066 ! INPUT:
4067 ! ----------------------------------------------------------------------
4068 integer n
4069 real a(n,3) ! a(*,1) lower diagonal (Ai,i-1)
4070 ! a(*,2) principal diagonal (Ai,i)
4071 ! a(*,3) upper diagonal (Ai,i+1)
4072 real c(n)
4074 ! ----------------------------------------------------------------------
4075 ! OUTPUT:
4076 ! ----------------------------------------------------------------------
4077 real x(n)
4079 ! ----------------------------------------------------------------------
4080 ! LOCAL:
4081 ! ----------------------------------------------------------------------
4082 integer i
4084 ! ----------------------------------------------------------------------
4085 ! END VARIABLES DEFINITIONS
4086 ! ----------------------------------------------------------------------
4088 do i=n-1,1,-1
4089 c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
4090 a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
4091 enddo
4093 do i=2,n
4094 c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
4095 enddo
4097 do i=1,n
4098 x(i)=c(i)/a(i,2)
4099 enddo
4101 return
4102 end subroutine invert
4105 ! ===6=8===============================================================72
4106 ! ===6=8===============================================================72
4108 subroutine flux_wall(ua,va,pt,da,ptw,ptwin,uva,vva,uvb,vvb, &
4109 sfw,sfwin,evb,drst,dt,cdrag)
4111 ! ----------------------------------------------------------------------
4112 ! This routine computes the surface sources or sinks of momentum, tke,
4113 ! and heat from vertical surfaces (walls).
4114 ! ----------------------------------------------------------------------
4115 implicit none
4117 ! INPUT:
4118 ! -----
4119 real drst ! street directions for the current urban class
4120 real da ! air density
4121 real pt ! potential temperature
4122 real ptw ! Walls potential temperatures
4123 real ptwin ! Windows potential temperatures
4124 real ua ! wind speed
4125 real va ! wind speed
4126 real dt !time step
4127 real cdrag
4128 ! OUTPUT:
4129 ! ------
4130 ! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
4131 ! vertical surfaces (walls).
4132 ! The fluxes can be computed as follow: Fluxes of X = A*X + B
4133 ! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
4135 real uva ! U (wind component) Vertical surfaces, A (implicit) term
4136 real uvb ! U (wind component) Vertical surfaces, B (explicit) term
4137 real vva ! V (wind component) Vertical surfaces, A (implicit) term
4138 real vvb ! V (wind component) Vertical surfaces, B (explicit) term
4139 real tva ! Temperature Vertical surfaces, A (implicit) term
4140 real tvb ! Temperature Vertical surfaces, B (explicit) term
4141 real evb ! Energy (TKE) Vertical surfaces, B (explicit) term
4142 real sfw ! Surfaces fluxes from the walls
4143 real sfwin ! Surfaces fluxes from the windows
4145 ! LOCAL:
4146 ! -----
4147 real hc
4148 real hcwin
4149 real u_ort
4150 real vett
4153 ! -------------------------
4154 ! END VARIABLES DEFINITIONS
4155 ! -------------------------
4157 vett=(ua**2+va**2)**.5
4159 u_ort=abs((cos(drst)*ua-sin(drst)*va))
4161 uva=-cdrag*u_ort/2.*cos(drst)*cos(drst)
4162 vva=-cdrag*u_ort/2.*sin(drst)*sin(drst)
4164 uvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*va
4165 vvb=cdrag*u_ort/2.*sin(drst)*cos(drst)*ua
4167 if (vett.lt.4.88) then
4168 hc=5.678*(1.09+0.23*(vett/0.3048))
4169 else
4170 hc=5.678*0.53*((vett/0.3048)**0.78)
4171 endif
4173 if (hc.gt.da*cp_u/dt)then
4174 hc=da*cp_u/dt
4175 endif
4177 if (vett.lt.4.88) then
4178 hcwin=5.678*(0.99+0.21*(vett/0.3048))
4179 else
4180 hcwin=5.678*0.50*((vett/0.3048)**0.78)
4181 endif
4183 if (hcwin.gt.da*cp_u/dt) then
4184 hcwin=da*cp_u/dt
4185 endif
4187 ! tvb=hc*ptw/da/cp_u
4188 ! tva=-hc/da/cp_u
4189 !!!!!!!!!!!!!!!!!!!!
4190 ! explicit
4192 sfw=hc*(pt-ptw)
4193 sfwin=hcwin*(pt-ptwin)
4196 evb=cdrag*(abs(u_ort)**3.)/2.
4198 return
4199 end subroutine flux_wall
4201 ! ===6=8===============================================================72
4202 ! ===6=8===============================================================72
4204 subroutine flux_flat_ground(dz,z0,ua,va,pt,pt0,ptg, &
4205 uhb,vhb,sf,ehb,da,qv,pr,rsg,qg,resg,rsveg,f1,f2,f3,f4,fh,ric,utot,gr_type)
4207 ! ----------------------------------------------------------------------
4208 ! Calculation of the flux at the ground
4209 ! Formulation of Louis (Louis, 1979)
4210 ! ----------------------------------------------------------------------
4212 implicit none
4214 real dz ! first vertical level
4215 real pt ! potential temperature
4216 real pt0 ! reference potential temperature
4217 real ptg ! ground potential temperature
4218 real ua ! wind speed
4219 real va ! wind speed
4220 real z0 ! Roughness length
4221 real da ! air density
4222 real qv ! specific humidity
4223 real pr ! pressure
4224 real rsg ! solar radiation
4225 real qg(ng_u) ! Ground Soil Moisture
4229 ! ----------------------------------------------------------------------
4230 ! OUTPUT:
4231 ! ----------------------------------------------------------------------
4232 ! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
4233 ! surfaces (roofs and street)
4234 ! The fluxes can be computed as follow: Fluxes of X = B
4235 ! Example: Momentum fluxes on horizontal surfaces = uhb_u
4236 real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
4237 real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
4238 ! real thb ! Temperature Horizontal surfaces, B (explicit) term
4239 real tva ! Temperature Vertical surfaces, A (implicit) term
4240 real tvb ! Temperature Vertical surfaces, B (explicit) term
4241 real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
4242 real sf
4243 real lf
4245 ! ----------------------------------------------------------------------
4246 ! LOCAL:
4247 ! ----------------------------------------------------------------------
4248 real aa,ah
4249 real z0t
4250 real al
4251 real buu
4252 real c
4253 real fbuw
4254 real fbpt
4255 real fh
4256 real fm
4257 real ric
4258 real tstar
4259 real qstar
4260 real ustar
4261 real utot
4262 real wstar
4263 real zz
4264 real qvsg,qvs,es,esa,fbqq
4265 real b,cm,ch,rr,tol
4266 parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
4268 real f
4269 real f1
4270 real f2
4271 real f3
4272 real f4
4273 real ta ! surface air temperature
4274 real tmp ! ground temperature
4275 real rsveg ! Stomatal resistance
4276 real resg
4277 real lai ! leaf area index
4278 real sdlim ! radiation limit at which photosyntesis start W/m2
4279 parameter(sdlim=100.)
4280 real rsmin ! Minimum stomatal resistance
4281 real rsmax ! Maximun stomatal resistance
4282 real qw
4283 parameter(qw=0.06)
4284 real qref
4285 parameter(qref=0.37)
4286 real hs
4287 parameter(hs=36.35)
4289 real dzg_u(ng_u) ! Layer sizes in the ground
4291 data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
4293 real gx,dzg_tot
4294 integer gr_type,iz
4295 ! ----------------------------------------------------------------------
4296 ! END VARIABLES DEFINITIONS
4297 ! ----------------------------------------------------------------------
4298 z0t=z0/10.
4299 if(gr_type.eq.1)then
4300 rsmin=40.
4301 rsmax=5000.
4302 lai=2.
4303 elseif(gr_type.eq.2)then
4304 rsmin=150.
4305 rsmax=5000.
4306 lai=3.
4307 endif
4308 ! computation of the ground temperature
4310 utot=(ua**2.+va**2.)**.5
4313 !!!! Louis formulation
4314 !
4315 ! compute the bulk Richardson Number
4317 zz=dz/2.
4319 ! if(tstar.lt.0.)then
4320 ! wstar=(-ustar*tstar*g*hii/pt)**(1./3.)
4321 ! else
4322 ! wstar=0.
4323 ! endif
4324 !
4325 ! if (utot.le.0.7*wstar) utot=max(0.7*wstar,0.00001)
4327 utot=max(utot,0.01)
4329 ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
4331 aa=vk/log(zz/z0)
4332 ah=vk/log(zz/z0t)
4334 ! determine the parameters fm and fh for stable, neutral and unstable conditions
4336 if(ric.gt.0)then
4337 fm=1/(1+0.5*b*ric)**2.
4338 fh=fm
4339 else
4340 c=b*cm*aa*aa*(zz/z0)**.5
4341 fm=1-b*ric/(1+c*(-ric)**.5)
4342 c=b*cm*aa*ah*(zz/z0t)**.5
4343 c=c*ch/cm
4344 fh=1-b*ric/(1+c*(-ric)**.5)
4345 endif
4347 fbuw=-aa*aa*utot*utot*fm
4348 fbpt=-aa*ah*utot*(pt-ptg)*fh/rr
4349 tmp=ptg*(pr/p0)**(rcp_u)-273.15
4350 es=6.11*(10.**(tmp*7.5/(237.7+tmp)))
4351 qvsg=0.62197*es/(0.01*pr-0.378*es)
4354 f=0.55*rsg/sdlim*2./lai
4356 f1=(f+rsmin/rsmax)/(1.+f)
4358 ta=pt*(pr/p0)**(rcp_u)-273.15
4359 esa=6.11*(10**(ta*7.5/(237.7+ta)))
4360 qvs=0.62197*esa/(0.01*pr-0.378*esa)
4362 f2= 1./(1.+hs*(qvs-qv))
4363 f3=1.-0.0016*(25.-ta)**2.
4364 f4=0.
4365 dzg_tot=0.
4366 do iz=1,ng_u
4367 gx=(qg(iz)-qw)/(qref-qw)
4368 if (gx.gt.1)gx=1.
4369 if (gx.lt.0)gx=0.
4370 f4=f4+gx*dzg_u(iz)
4371 dzg_tot=dzg_tot+dzg_u(iz)
4372 enddo
4373 f4=f4/dzg_tot
4375 rsveg=min(rsmin/max(lai*f1*f2*f3*f4,1e-9),rsmax)
4376 resg= rr/(aa*aa*utot*fh)
4379 fbqq=-(qv-qvsg)/(resg+rsveg)
4382 ustar=(-fbuw)**.5
4383 tstar=-fbpt/ustar
4384 qstar=-fbqq/ustar
4386 al=(vk*g_u*tstar)/(pt*ustar*ustar)
4388 buu=-g_u/pt0*ustar*tstar
4390 uhb=-ustar*ustar*ua/utot
4391 vhb=-ustar*ustar*va/utot
4392 sf= ustar*tstar*da*cp_u
4393 lf= ustar*qstar*da*latent
4395 ! thb= 0.
4396 ehb=buu
4397 !!!!!!!!!!!!!!!
4399 return
4400 end subroutine flux_flat_ground
4402 ! ===6=8===============================================================72
4403 ! ===6=8===============================================================72
4404 subroutine flux_flat_roof(dz,z0,ua,va,pt,pt0,ptg, &
4405 uhb,vhb,sf,lf,ehb,da,qv,pr,rsg,qr,resg,rsveg,f1,f2,f3,f4,gr_type,pv_frac_roof)
4407 ! ----------------------------------------------------------------------
4408 ! Calculation of the flux at the ground
4409 ! Formulation of Louis (Louis, 1979)
4410 ! ----------------------------------------------------------------------
4412 implicit none
4414 real dz ! first vertical level
4415 real pt ! potential temperature
4416 real pt0 ! reference potential temperature
4417 real ptg ! ground potential temperature
4418 real ua ! wind speed
4419 real va ! wind speed
4420 real z0 ! Roughness length
4421 real da ! air density
4422 real qv ! specific humidity
4423 real pr ! pressure
4424 real rsg ! solar radiation
4425 real qr(ngr_u) ! Ground Soil Moisture
4426 real pv_frac_roof
4427 real rs_eff
4429 ! ----------------------------------------------------------------------
4430 ! OUTPUT:
4431 ! ----------------------------------------------------------------------
4432 ! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
4433 ! surfaces (roofs and street)
4434 ! The fluxes can be computed as follow: Fluxes of X = B
4435 ! Example: Momentum fluxes on horizontal surfaces = uhb_u
4436 real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
4437 real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
4438 ! real thb ! Temperature Horizontal surfaces, B (explicit) term
4439 real tva ! Temperature Vertical surfaces, A (implicit) term
4440 real tvb ! Temperature Vertical surfaces, B (explicit) term
4441 real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
4442 real sf
4443 real lf
4445 ! ----------------------------------------------------------------------
4446 ! LOCAL:
4447 ! ----------------------------------------------------------------------
4448 real aa,ah
4449 real al
4450 real buu
4451 real c
4452 real fbuw
4453 real fbpt
4454 real fh
4455 real fm
4456 real ric
4457 real tstar
4458 real qstar
4459 real ustar
4460 real utot
4461 real wstar
4462 real zz
4463 real z0t
4464 real qvsg,qvs,es,esa,fbqq
4465 real b,cm,ch,rr,tol
4466 parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
4468 real f
4469 real f1
4470 real f2
4471 real f3
4472 real f4
4473 real ta ! surface air temperature
4474 real tmp ! ground temperature
4475 real rsveg ! Stomatal resistance
4476 real resg
4477 real lai ! leaft area index
4478 real sdlim ! radiation limit at which photosyntesis start W/m2
4479 parameter(sdlim=100.)
4480 real rsmin
4481 real rsmax ! Maximun stomatal resistance
4482 real qw ! Wilting point
4483 parameter(qw=0.06)
4484 real qref ! Field capacity
4485 parameter(qref=0.37)
4486 real hs
4487 parameter(hs=36.35)
4489 real dzgr_u(ngr_u) ! Layer sizes in the ground
4491 data dzgr_u /0.1,0.003,0.06,0.003,0.05,0.04,0.02,0.0125,0.005,0.0025/
4493 real gx,dzgr_tot
4494 integer gr_type,iz
4495 ! ----------------------------------------------------------------------
4496 ! END VARIABLES DEFINITIONS
4498 ! ----------------------------------------------------------------------
4499 z0t=z0/10.
4500 if(gr_type.eq.1)then
4501 rsmin=40.
4502 rsmax=5000.
4503 lai=2.
4504 elseif(gr_type.eq.2)then
4505 rsmin=150.
4506 rsmax=5000.
4507 lai=3.
4508 endif
4509 rs_eff=(1-pv_frac_roof)*rsg
4510 ! computation of the ground temperature
4512 utot=(ua**2.+va**2.)**.5
4514 !!!! Louis formulation
4515 !
4516 ! compute the bulk Richardson Number
4518 zz=dz/2.
4521 utot=max(utot,0.01)
4523 ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
4525 aa=vk/log(zz/z0)
4526 ah=vk/log(zz/z0t)
4528 if(ric.gt.0.)then
4529 fm=1./(1.+0.5*b*ric)**2.
4530 fh=fm
4531 else
4532 c=b*cm*aa*aa*(zz/z0)**.5
4533 fm=1.-b*ric/(1.+c*(-ric)**.5)
4534 c=b*cm*aa*ah*(zz/z0t)**.5
4535 c=c*ch/cm
4536 fh=1.-b*ric/(1+c*(-ric)**.5)
4537 endif
4539 fbuw=-aa*aa*utot*utot*fm
4540 fbpt=-aa*ah*utot*(pt-ptg)*fh/rr
4541 tmp=ptg*(pr/p0)**(rcp_u)-273.15
4542 es=6.11*(10.**(tmp*7.5/(237.7+tmp)))
4543 qvsg=0.62197*es/(0.01*pr-0.378*es)
4546 f=0.55*rs_eff/sdlim*2./lai
4548 f1=(f+rsmin/rsmax)/(1.+f)
4550 ta=pt*(pr/p0)**(rcp_u)-273.15
4551 esa=6.11*(10**(ta*7.5/(237.7+ta)))
4552 qvs=0.62197*esa/(0.01*pr-0.378*esa)
4554 f2= 1./(1.+hs*(qvs-qv))
4555 f3=1.-0.0016*(25.-ta)**2.
4556 f4=0.
4557 dzgr_tot=0.
4558 do iz=5,ngr_u
4559 gx=(qr(iz)-qw)/(qref-qw)
4560 if (gx.gt.1)gx=1.
4561 if (gx.lt.0)gx=0.
4562 f4=f4+gx*dzgr_u(iz)
4563 dzgr_tot=dzgr_tot+dzgr_u(iz)
4564 enddo
4565 f4=f4/dzgr_tot
4567 rsveg=min(rsmin/max(lai*f1*f2*f3*f4,1e-9),rsmax)
4570 resg= rr/(aa*aa*utot*fh)
4573 fbqq=-(qv-qvsg)/(resg+rsveg)
4575 ustar=(-fbuw)**.5
4576 tstar=-fbpt/ustar
4577 qstar=-fbqq/ustar
4579 al=(vk*g_u*tstar)/(pt*ustar*ustar)
4581 buu=-g_u/pt0*ustar*tstar
4583 uhb=-ustar*ustar*ua/utot
4584 vhb=-ustar*ustar*va/utot
4585 sf= ustar*tstar*da*cp_u
4586 lf= ustar*qstar*da*latent
4588 ehb=buu
4589 end subroutine flux_flat_roof
4591 !!!!!!!===============================
4593 ! ===6=8===============================================================72
4594 ! ===6=8===============================================================72
4596 subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg,qv, &
4597 uhb,vhb,sf,lf,ehb,da,pr)
4599 ! ----------------------------------------------------------------------
4600 ! Calculation of the flux at the ground
4601 ! Formulation of Louis (Louis, 1979)
4602 ! ----------------------------------------------------------------------
4604 implicit none
4605 real pr
4606 real dz ! first vertical level
4607 real pt ! potential temperature
4608 real pt0 ! reference potential temperature
4609 real ptg ! ground potential temperature
4610 real ua ! wind speed
4611 real va ! wind speed
4612 real z0 ! Roughness length
4613 real da ! air density
4614 real qv
4615 ! ----------------------------------------------------------------------
4616 ! OUTPUT:
4617 ! ----------------------------------------------------------------------
4618 ! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
4619 ! surfaces (roofs and street)
4620 ! The fluxes can be computed as follow: Fluxes of X = B
4621 ! Example: Momentum fluxes on horizontal surfaces = uhb_u
4622 real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
4623 real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
4624 ! real thb ! Temperature Horizontal surfaces, B (explicit) term
4625 real tva ! Temperature Vertical surfaces, A (implicit) term
4626 real tvb ! Temperature Vertical surfaces, B (explicit) term
4627 real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
4628 real sf
4629 real lf
4631 ! ----------------------------------------------------------------------
4632 ! LOCAL:
4633 ! ----------------------------------------------------------------------
4634 real aa
4635 real al
4636 real buu
4637 real c
4638 real fbuw
4639 real fbpt
4640 real fh
4641 real fm
4642 real ric
4643 real tstar
4644 real ustar
4645 real qstar
4646 real utot
4647 real wstar
4648 real zz
4649 real qvsg,qvs,es,esa,fbqq,tmp,resg
4650 real b,cm,ch,rr,tol
4651 parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
4653 ! ----------------------------------------------------------------------
4654 ! END VARIABLES DEFINITIONS
4655 ! ----------------------------------------------------------------------
4658 ! computation of the ground temperature
4660 utot=(ua**2+va**2)**.5
4663 !!!! Louis formulation
4664 !
4665 ! compute the bulk Richardson Number
4667 zz=dz/2.
4670 utot=max(utot,0.01)
4672 ric=2.*g_u*zz*(pt-ptg)/((pt+ptg)*(utot**2))
4674 aa=vk/log(zz/z0)
4678 tmp=ptg*(pr/(1.e+5))**(rcp_u)-273.15
4679 es=6.11*(10**(tmp*7.5/(237.7+tmp)))
4680 qvsg=0.62197*es/(0.01*pr-0.378*es)
4684 ! determine the parameters fm and fh for stable, neutral and unstable conditions
4686 if(ric.gt.0.)then
4687 fm=1./(1.+0.5*b*ric)**2
4688 fh=fm
4689 else
4690 c=b*cm*aa*aa*(zz/z0)**.5
4691 fm=1.-b*ric/(1.+c*(-ric)**.5)
4692 c=c*ch/cm
4693 fh=1.-b*ric/(1.+c*(-ric)**.5)
4694 endif
4696 resg= rr/(aa*aa*utot*fh)
4697 fbuw=-aa*aa*utot*utot*fm
4698 fbpt=-aa*aa*utot*(pt-ptg)*fh/rr
4699 fbqq=-(qv-qvsg)/(resg)
4701 ustar=(-fbuw)**.5
4702 tstar=-fbpt/ustar
4703 qstar=-fbqq/ustar
4704 al=(vk*g_u*tstar)/(pt*ustar*ustar)
4706 buu=-g_u/pt0*ustar*tstar
4708 uhb=-ustar*ustar*ua/utot
4709 vhb=-ustar*ustar*va/utot
4710 sf= ustar*tstar*da*cp_u
4711 lf= ustar*qstar*da*latent
4712 ehb=buu
4713 !!!!!!!!!!!!!!!
4715 return
4716 end subroutine flux_flat
4717 !!!!!!!!!!!!!================!!!!!!!!!!!!!!!!!!!
4718 ! ===6=8===============================================================72
4719 ! ===6=8===============================================================72
4721 subroutine icBEP (nd_u,h_b,d_b,ss_u,pb_u,nz_u,z_u)
4723 implicit none
4725 ! Street parameters
4726 integer nd_u(nurbm) ! Number of street direction for each urban class
4727 real h_b(nz_um,nurbm) ! Bulding's heights [m]
4728 real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
4729 ! -----------------------------------------------------------------------
4730 ! Output
4731 !------------------------------------------------------------------------
4733 real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to z
4734 real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to z
4736 ! Grid parameters
4737 integer nz_u(nurbm) ! Number of layer in the urban grid
4738 real z_u(nz_um) ! Height of the urban grid levels
4741 ! -----------------------------------------------------------------------
4742 ! Local
4743 !------------------------------------------------------------------------
4745 integer iz_u,id,ilu,iurb
4747 real dtot
4748 real hbmax
4750 ! -----------------------------------------------------------------------
4751 ! This routine initialise the urban paramters for the BEP module
4752 !------------------------------------------------------------------------
4753 !
4754 !Initialize variables
4755 !
4756 !
4757 nz_u=0
4758 z_u=0.
4759 ss_u=0.
4760 pb_u=0.
4762 ! Computation of the urban levels height
4764 z_u(1)=0.
4766 do iz_u=1,nz_um-1
4767 z_u(iz_u+1)=z_u(iz_u)+dz_u
4768 enddo
4770 ! Normalisation of the building density
4772 do iurb=1,nurbm
4773 dtot=0.
4774 do ilu=1,nz_um
4775 dtot=dtot+d_b(ilu,iurb)
4776 enddo
4777 do ilu=1,nz_um
4778 d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
4779 enddo
4780 enddo
4782 ! Compute the view factors, pb and ss
4784 do iurb=1,nurbm
4785 hbmax=0.
4786 nz_u(iurb)=0
4787 do ilu=1,nz_um
4788 if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
4789 enddo
4791 do iz_u=1,nz_um-1
4792 if(z_u(iz_u+1).gt.hbmax)go to 10
4793 enddo
4795 10 continue
4796 nz_u(iurb)=iz_u+1
4798 do id=1,nd_u(iurb)
4800 do iz_u=1,nz_u(iurb)
4801 ss_u(iz_u,iurb)=0.
4802 do ilu=1,nz_um
4803 if(z_u(iz_u).le.h_b(ilu,iurb) &
4804 .and.z_u(iz_u+1).gt.h_b(ilu,iurb))then
4805 ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
4806 endif
4807 enddo
4808 enddo
4810 pb_u(1,iurb)=1.
4811 do iz_u=1,nz_u(iurb)
4812 pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
4813 enddo
4815 enddo
4816 end do
4819 return
4820 end subroutine icBEP
4822 ! ===6=8===============================================================72
4823 ! ===6=8===============================================================72
4826 subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws)
4828 implicit none
4832 ! -----------------------------------------------------------------------
4833 ! Input
4834 !------------------------------------------------------------------------
4836 integer iurb ! Number of the urban class
4837 integer nz_u ! Number of levels in the urban grid
4838 integer id ! Street direction number
4839 real ws ! Street width
4840 real z(nz_um) ! Height of the urban grid levels
4841 real dxy ! Street lenght
4844 ! -----------------------------------------------------------------------
4845 ! Output
4846 !------------------------------------------------------------------------
4848 ! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
4849 ! and the short wave radation. They are the part of radiation from a surface
4850 ! or from the sky to another surface.
4852 real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
4853 real fwg(nz_um,ndm,nurbm) ! from wall to ground
4854 real fgw(nz_um,ndm,nurbm) ! from ground to wall
4855 real fsw(nz_um,ndm,nurbm) ! from sky to wall
4856 real fws(nz_um,ndm,nurbm) ! from wall to sky
4857 real fsg(ndm,nurbm) ! from sky to ground
4860 ! -----------------------------------------------------------------------
4861 ! Local
4862 !------------------------------------------------------------------------
4864 integer jz,iz
4866 real hut
4867 real f1,f2,f12,f23,f123,ftot
4868 real fprl,fnrm
4869 real a1,a2,a3,a4,a12,a23,a123
4871 ! -----------------------------------------------------------------------
4872 ! This routine calculates the view factors
4873 !------------------------------------------------------------------------
4875 hut=z(nz_u+1)
4877 do jz=1,nz_u
4879 ! radiation from wall to wall
4881 do iz=1,nz_u
4883 call fprls (fprl,dxy,abs(z(jz+1)-z(iz )),ws)
4884 f123=fprl
4885 call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
4886 f23=fprl
4887 call fprls (fprl,dxy,abs(z(jz )-z(iz )),ws)
4888 f12=fprl
4889 call fprls (fprl,dxy,abs(z(jz )-z(iz+1)),ws)
4890 f2 = fprl
4892 a123=dxy*(abs(z(jz+1)-z(iz )))
4893 a12 =dxy*(abs(z(jz )-z(iz )))
4894 a23 =dxy*(abs(z(jz+1)-z(iz+1)))
4895 a1 =dxy*(abs(z(iz+1)-z(iz )))
4896 a2 =dxy*(abs(z(jz )-z(iz+1)))
4897 a3 =dxy*(abs(z(jz+1)-z(jz )))
4899 ftot=0.5*(a123*f123-a23*f23-a12*f12+a2*f2)/a1
4901 fww(iz,jz,id,iurb)=ftot*a1/a3
4903 enddo
4905 ! radiation from ground to wall
4907 call fnrms (fnrm,z(jz+1),dxy,ws)
4908 f12=fnrm
4909 call fnrms (fnrm,z(jz) ,dxy,ws)
4910 f1=fnrm
4912 a1 = ws*dxy
4914 a12= ws*dxy
4916 a4=(z(jz+1)-z(jz))*dxy
4918 ftot=(a12*f12-a12*f1)/a1
4920 fgw(jz,id,iurb)=ftot*a1/a4
4922 ! radiation from sky to wall
4924 call fnrms(fnrm,hut-z(jz) ,dxy,ws)
4925 f12 = fnrm
4926 call fnrms (fnrm,hut-z(jz+1),dxy,ws)
4927 f1 =fnrm
4929 a1 = ws*dxy
4931 a12= ws*dxy
4933 a4 = (z(jz+1)-z(jz))*dxy
4935 ftot=(a12*f12-a12*f1)/a1
4937 fsw(jz,id,iurb)=ftot*a1/a4
4939 enddo
4941 ! radiation from wall to sky
4942 do iz=1,nz_u
4943 call fnrms(fnrm,ws,dxy,hut-z(iz))
4944 f12=fnrm
4945 call fnrms(fnrm,ws,dxy,hut-z(iz+1))
4946 f1=fnrm
4947 a1 = (z(iz+1)-z(iz))*dxy
4948 a2 = (hut-z(iz+1))*dxy
4949 a12= (hut-z(iz))*dxy
4950 a4 = ws*dxy
4951 ftot=(a12*f12-a2*f1)/a1
4952 fws(iz,id,iurb)=ftot*a1/a4
4954 enddo
4955 !!!!!!!!!!!!!
4958 do iz=1,nz_u
4960 ! radiation from wall to ground
4962 call fnrms (fnrm,ws,dxy,z(iz+1))
4963 f12=fnrm
4964 call fnrms (fnrm,ws,dxy,z(iz ))
4965 f1 =fnrm
4967 a1= (z(iz+1)-z(iz) )*dxy
4969 a2 = z(iz)*dxy
4970 a12= z(iz+1)*dxy
4971 a4 = ws*dxy
4973 ftot=(a12*f12-a2*f1)/a1
4975 fwg(iz,id,iurb)=ftot*a1/a4
4977 enddo
4979 ! radiation from sky to ground
4981 call fprls (fprl,dxy,ws,hut)
4982 fsg(id,iurb)=fprl
4984 return
4985 end subroutine view_factors
4987 ! ===6=8===============================================================72
4988 ! ===6=8===============================================================72
4991 SUBROUTINE fprls (fprl,a,b,c)
4993 implicit none
4997 real a,b,c
4998 real x,y
4999 real fprl
5002 x=a/c
5003 y=b/c
5005 if(a.eq.0.or.b.eq.0.)then
5006 fprl=0.
5007 else
5008 fprl=log( ( (1.+x**2)*(1.+y**2)/(1.+x**2+y**2) )**.5)+ &
5009 y*((1.+x**2)**.5)*atan(y/((1.+x**2)**.5))+ &
5010 x*((1.+y**2)**.5)*atan(x/((1.+y**2)**.5))- &
5011 y*atan(y)-x*atan(x)
5012 fprl=fprl*2./(pi*x*y)
5013 endif
5015 return
5016 end subroutine fprls
5018 ! ===6=8===============================================================72
5019 ! ===6=8===============================================================72
5021 SUBROUTINE fnrms (fnrm,a,b,c)
5023 implicit none
5027 real a,b,c
5028 real x,y,z,a1,a2,a3,a4,a5,a6
5029 real fnrm
5031 x=a/b
5032 y=c/b
5033 z=x**2+y**2
5035 if(y.eq.0.or.x.eq.0)then
5036 fnrm=0.
5037 else
5038 a1=log( (1.+x*x)*(1.+y*y)/(1.+z) )
5039 a2=y*y*log(y*y*(1.+z)/z/(1.+y*y) )
5040 a3=x*x*log(x*x*(1.+z)/z/(1.+x*x) )
5041 a4=y*atan(1./y)
5042 a5=x*atan(1./x)
5043 a6=sqrt(z)*atan(1./sqrt(z))
5044 fnrm=0.25*(a1+a2+a3)+a4+a5-a6
5045 fnrm=fnrm/(pi*y)
5046 endif
5048 return
5049 end subroutine fnrms
5050 ! ===6=8===============================================================72
5052 SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
5053 twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
5054 emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b, &
5055 cop_u,pwin_u,beta_u,sw_cond_u,time_on_u,time_off_u,targtemp_u, &
5056 bldac_frc_u,cooled_frc_u, &
5057 gaptemp_u, targhum_u,gaphum_u,perflo_u, &
5058 gr_frac_roof_u,pv_frac_roof_u, &
5059 hsesf_u,hsequip,irho,gr_flag_u,gr_type_u)
5062 ! initialization routine, where the variables from the table are read
5064 implicit none
5065 integer iurb ! urban class number
5066 ! Building parameters
5067 real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
5068 real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
5069 real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
5070 real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
5071 real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
5072 real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
5073 real twini_u(nurbm) ! Temperature inside the buildings behind the wall [K]
5074 real trini_u(nurbm) ! Temperature inside the buildings behind the roof [K]
5075 real tgini_u(nurbm) ! Initial road temperature
5077 ! Radiation parameters
5078 real albg_u(nurbm) ! Albedo of the ground
5079 real albw_u(nurbm) ! Albedo of the wall
5080 real albr_u(nurbm) ! Albedo of the roof
5081 real albwin_u(nurbm) ! Albedo of the window
5082 real emg_u(nurbm) ! Emissivity of ground
5083 real emw_u(nurbm) ! Emissivity of wall
5084 real emr_u(nurbm) ! Emissivity of roof
5085 real emwind_u(nurbm) ! Emissivity of windows
5087 ! Roughness parameters
5088 real z0g_u(nurbm) ! The ground's roughness length
5089 real z0r_u(nurbm) ! The roof's roughness length
5091 ! Street parameters
5092 integer nd_u(nurbm) ! Number of street direction for each urban class
5094 real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
5095 real drst_u(ndm,nurbm) ! Street direction [degree]
5096 real ws_u(ndm,nurbm) ! Street width [m]
5097 real bs_u(ndm,nurbm) ! Building width [m]
5098 real h_b(nz_um,nurbm) ! Bulding's heights [m]
5099 real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
5101 integer i,iu
5102 integer nurb ! number of urban classes used
5103 real, intent(out) :: bldac_frc_u(nurbm)
5104 real, intent(out) :: cooled_frc_u(nurbm)
5105 real, intent(out) :: cop_u(nurbm)
5106 real, intent(out) :: pwin_u(nurbm)
5107 real, intent(out) :: beta_u(nurbm)
5108 integer, intent(out) :: sw_cond_u(nurbm)
5109 real, intent(out) :: time_on_u(nurbm)
5110 real, intent(out) :: time_off_u(nurbm)
5111 real, intent(out) :: targtemp_u(nurbm)
5112 real, intent(out) :: gaptemp_u(nurbm)
5113 real, intent(out) :: targhum_u(nurbm)
5114 real, intent(out) :: gaphum_u(nurbm)
5115 real, intent(out) :: perflo_u(nurbm)
5116 real, intent(out) :: gr_frac_roof_u(nurbm)
5117 real, intent(out) :: pv_frac_roof_u(nurbm)
5118 real, intent(out) :: hsesf_u(nurbm)
5119 real, intent(out) :: hsequip(24)
5120 real, intent(out) :: irho(24)
5121 integer, intent(out) :: gr_flag_u,gr_type_u
5122 !
5123 !Initialize some variables
5124 !
5126 h_b=0.
5127 d_b=0.
5129 nurb=ICATE
5130 do iu=1,nurb
5131 nd_u(iu)=0
5132 enddo
5134 csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
5135 csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
5136 csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
5137 do i=1,icate
5138 alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
5139 alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
5140 alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
5141 enddo
5142 twini_u=TBLEND_TBL
5143 trini_u=TRLEND_TBL
5144 tgini_u=TGLEND_TBL
5145 albw_u=ALBB_TBL
5146 albr_u=ALBR_TBL
5147 albg_u=ALBG_TBL
5148 emw_u=EPSB_TBL
5149 emr_u=EPSR_TBL
5150 emg_u=EPSG_TBL
5151 z0r_u=Z0R_TBL
5152 z0g_u=Z0G_TBL
5153 nd_u=NUMDIR_TBL
5154 !FS
5155 ! print*, 'g alla call', gr_frac_roof_u(iurb)
5156 bldac_frc_u = bldac_frc_tbl
5157 cooled_frc_u = cooled_frc_tbl
5158 cop_u = cop_tbl
5159 pwin_u = pwin_tbl
5160 beta_u = beta_tbl
5161 sw_cond_u = sw_cond_tbl
5162 time_on_u = time_on_tbl
5163 time_off_u = time_off_tbl
5164 targtemp_u = targtemp_tbl
5165 gaptemp_u = gaptemp_tbl
5166 targhum_u = targhum_tbl
5167 gaphum_u = gaphum_tbl
5168 perflo_u = perflo_tbl
5169 gr_frac_roof_u =gr_frac_roof_tbl
5170 gr_flag_u=gr_flag_tbl
5171 pv_frac_roof_u = pv_frac_roof_tbl
5172 hsesf_u = hsesf_tbl
5173 hsequip = hsequip_tbl
5174 irho=irho_tbl
5175 gr_type_u=gr_type_tbl
5176 do iu=1,icate
5177 if(ndm.lt.nd_u(iu))then
5178 write(*,*)'ndm too small in module_sf_bep_bem, please increase to at least ', nd_u(iu)
5179 write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
5180 stop
5181 endif
5182 do i=1,nd_u(iu)
5183 drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
5184 ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
5185 bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
5186 enddo
5187 enddo
5188 do iu=1,ICATE
5189 if(nz_um.lt.numhgt_tbl(iu)+3)then
5190 write(*,*)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
5191 write(*,*)'remember also that num_urban_layers should be equal or greater than nz_um*ndm*nwr-u!'
5192 stop
5193 endif
5194 do i=1,NUMHGT_TBL(iu)
5195 h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
5196 d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
5197 enddo
5198 enddo
5200 do i=1,ndm
5201 do iu=1,nurbm
5202 strd_u(i,iu)=100000.
5203 enddo
5204 enddo
5206 do iu=1,nurb
5207 emwind_u(iu)=0.9
5208 call albwindow(albwin_u(iu))
5209 enddo
5211 return
5212 end subroutine init_para
5213 !==============================================================
5214 !==============================================================
5215 !====6=8===============================================================72
5216 !====6=8===============================================================72
5218 subroutine upward_rad(ndu,nzu,ws,bs,sigma,pb,ss, &
5219 tg_av,emg_u,albg_u,rlg,rsg,sfg,lfg, &
5220 tw,emw_u,albw_u,rlw,rsw,sfw, &
5221 tr_av,emr_u,albr_u,emwind,albwind,twlev,pwin, &
5222 sfwind,rld,rs, sfr,sfrv,lfr,lfrv, &
5223 rs_abs,rl_up,emiss,grdflx_urb,gr_frac_roof,tpvlev,pv_frac_roof)
5224 !
5225 ! IN this surboutine we compute the upward longwave flux, and the albedo
5226 ! needed for the radiation scheme
5227 !
5228 implicit none
5230 !
5231 !INPUT VARIABLES
5232 !
5233 real rsw(2*ndm,nz_um) ! Short wave radiation at the wall for a given canyon direction [W/m2]
5234 real rlw(2*ndm,nz_um) ! Long wave radiation at the walls for a given canyon direction [W/m2]
5235 real rsg(ndm) ! Short wave radiation at the canyon for a given canyon direction [W/m2]
5236 real rlg(ndm) ! Long wave radiation at the ground for a given canyon direction [W/m2]
5237 real rs ! Short wave radiation at the horizontal surface from the sun [W/m2]
5238 real sfw(2*ndm,nz_um) ! Sensible heat flux from walls [W/m2]
5239 real sfg(ndm) ! Sensible heat flux from ground (road) [W/m2]
5240 real lfg(ndm)
5241 real sfr(ndm,nz_um) ! Sensible heat flux from roofs [W/m2]
5242 real lfr(ndm,nz_um)
5243 real lfrv(ndm,nz_um)
5244 real sfrv(ndm,nz_um)
5245 real gr_frac_roof
5246 real rld ! Long wave radiation from the sky [W/m2]
5247 real albg_u ! albedo of the ground/street
5248 real albw_u ! albedo of the walls
5249 real albr_u ! albedo of the roof
5250 real ws(ndm) ! width of the street
5251 real bs(ndm)
5252 ! building size
5253 real pb(nz_um) ! Probability to have a building with an height equal or higher
5254 integer nzu
5255 real ss(nz_um) ! Probability to have a building of a given height
5256 real sigma
5257 real emg_u ! emissivity of the street
5258 real emw_u ! emissivity of the wall
5259 real emr_u ! emissivity of the roof
5260 real tw(2*ndm,nz_um) ! Temperature in each layer of the wall [K]
5261 real tr_av(ndm,nz_um) ! Temperature in each layer of the roof [K]
5262 real tpvlev(ndm,nz_um)
5263 real pv_frac_roof
5264 real tg_av(ndm) ! Temperature in each layer of the ground [K]
5265 integer id ! street direction
5266 integer ndu ! number of street directions
5267 !
5268 !New variables BEM
5269 !
5270 real emwind !Emissivity of the windows
5271 real albwind !Albedo of the windows
5272 real twlev(2*ndm,nz_um) !Averaged Temperature of the windows
5273 real pwin !Coverage area fraction of the windows
5274 real gflwin !Heat stored for the windows
5275 real sfwind(2*ndm,nz_um) !Sensible heat flux from windows [W/m2]
5277 !OUTPUT/INPUT
5278 real rs_abs ! absrobed solar radiationfor this street direction
5279 real rl_up ! upward longwave radiation for this street direction
5280 real emiss ! mean emissivity
5281 real grdflx_urb ! ground heat flux
5282 !LOCAL
5283 integer iz,iw
5284 real rl_inc,rl_emit
5285 real gfl
5286 integer ix,iy,iwrong
5288 iwrong=1
5289 do iz=1,nzu+1
5290 do id=1,ndu
5291 if(tr_av(id,iz).lt.100.)then
5292 write(203,*) tr_av(id,iz)
5293 write(*,*)'in upward_rad ',iz,id,iw,tr_av(id,iz)
5294 iwrong=0
5295 endif
5296 enddo
5297 enddo
5298 if(iwrong.eq.0)stop
5300 rl_up=0.
5302 rs_abs=0.
5303 rl_inc=0.
5304 emiss=0.
5305 rl_emit=0.
5306 grdflx_urb=0.
5307 do id=1,ndu
5308 rl_emit=rl_emit-( emg_u*sigma*(tg_av(id)**4.)+(1-emg_u)*rlg(id))*ws(id)/(ws(id)+bs(id))/ndu
5309 rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/ndu
5310 rs_abs=rs_abs+(1.-albg_u)*rsg(id)*ws(id)/(ws(id)+bs(id))/ndu
5311 gfl=(1.-albg_u)*rsg(id)+emg_u*rlg(id)-emg_u*sigma*(tg_av(id)**4.)+sfg(id)+lfg(id)
5312 grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/ndu
5314 do iz=2,nzu
5315 rl_emit=rl_emit-(emr_u*sigma*(1.-pv_frac_roof)*tr_av(id,iz)**4.+0.79*sigma*pv_frac_roof*tpvlev(id,iz)**4+ &
5316 (1-emr_u)*rld*(1.-pv_frac_roof)+(1-0.79)*pv_frac_roof*rld)*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
5317 rl_inc=rl_inc+rld*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
5318 rs_abs=rs_abs+((1.-albr_u)*rs*(1.-pv_frac_roof)+(1.-0.11)*rs*pv_frac_roof)*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
5319 gfl=(1.-albr_u)*rs*(1-pv_frac_roof)+emr_u*rld*(1-pv_frac_roof)+pv_frac_roof*emr_u*sigma*tpvlev(id,iz)**4 &
5320 -emr_u*sigma*(tr_av(id,iz)**4.)+(1-gr_frac_roof)*sfr(id,iz)+(sfrv(id,iz)+lfrv(id,iz))*gr_frac_roof+(1.-gr_frac_roof)*lfr(id,iz)
5321 grdflx_urb=grdflx_urb-gfl*ss(iz)*bs(id)/(ws(id)+bs(id))/ndu
5322 enddo
5324 do iz=1,nzu
5326 rl_emit=rl_emit-(emw_u*(1.-pwin)*sigma*(tw(2*id-1,iz)**4.+tw(2*id,iz)**4.)+ &
5327 (emwind*pwin*sigma*(twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.))+ &
5328 ((1.-emw_u)*(1.-pwin)+pwin*(1.-emwind))*(rlw(2*id-1,iz)+rlw(2*id,iz)))* &
5329 dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
5331 rl_inc=rl_inc+((rlw(2*id-1,iz)+rlw(2*id,iz)))*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
5333 rs_abs=rs_abs+(((1.-albw_u)*(1.-pwin)+(1.-albwind)*pwin)*(rsw(2*id-1,iz)+rsw(2*id,iz)))*&
5334 dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
5336 gfl=(1.-albw_u)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emw_u*( rlw(2*id-1,iz)+rlw(2*id,iz) ) &
5337 -emw_u*sigma*( tw(2*id-1,iz)**4.+tw(2*id,iz)**4. )+(sfw(2*id-1,iz)+sfw(2*id,iz))
5339 gflwin=(1.-albwind)*(rsw(2*id-1,iz)+rsw(2*id,iz)) +emwind*(rlw(2*id-1,iz)+rlw(2*id,iz)) &
5340 -emwind*sigma*( twlev(2*id-1,iz)**4.+twlev(2*id,iz)**4.)+(sfwind(2*id-1,iz)+sfwind(2*id,iz))
5343 grdflx_urb=grdflx_urb-(gfl*(1.-pwin)+pwin*gflwin)*dz_u*pb(iz+1)/(ws(id)+bs(id))/ndu
5345 enddo
5347 enddo
5348 emiss=(emg_u+emw_u+emr_u)/3.
5349 rl_up=(rl_inc+rl_emit)-rld
5352 return
5354 END SUBROUTINE upward_rad
5356 !====6=8===============================================================72
5357 !====6=8===============================================================72
5358 ! ===6================================================================72
5359 ! ===6================================================================72
5361 subroutine albwindow(albwin)
5363 !-------------------------------------------------------------------
5364 implicit none
5367 ! -------------------------------------------------------------------
5368 !Based on the
5369 !paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour
5370 !of the total solar energy transmittance of windows"
5371 !Solar Energy Vol.69,No.4,pp. 321-329.
5372 ! -------------------------------------------------------------------
5373 !Input
5374 !-----
5376 !Output
5377 !------
5378 real albwin ! albedo of the window
5379 !Local
5380 !-----
5381 real a,b,c !Polynomial coefficients
5382 real alfa,delta,gama !Polynomial powers
5383 real g0 !transmittance when the angle
5384 !of incidence is normal to the surface.
5385 real asup,ainf
5386 real fonc
5388 !Constants
5389 !--------------------
5391 real epsilon !accuracy of the integration
5392 parameter (epsilon=1.e-07)
5393 real n1,n2 !Index of refraction for glasses and air
5394 parameter(n1=1.,n2=1.5)
5395 integer intg,k
5396 !--------------------------------------------------------------------
5397 if (q_num.eq.0) then
5398 write(*,*) 'Category parameter of the windows no valid'
5399 stop
5400 endif
5402 g0=4.*n1*n2/((n1+n2)*(n1+n2))
5403 a=8.
5404 b=0.25/q_num
5405 c=1.-a-b
5406 alfa =5.2 + (0.7*q_num)
5407 delta =2.
5408 gama =(5.26+0.06*p_num)+(0.73+0.04*p_num)*q_num
5410 intg=1
5411 !----------------------------------------------------------------------
5414 100 asup=0.
5415 ainf=0.
5417 do k=1,intg
5418 call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
5419 asup=asup+(pi/intg)*fonc
5420 enddo
5422 intg=intg+1
5424 do k=1,intg
5425 call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
5426 ainf=ainf+(pi/intg)*fonc
5427 enddo
5429 if(abs(asup-ainf).lt.epsilon) then
5430 albwin=1-g0+(g0/2.)*asup
5431 else
5432 goto 100
5433 endif
5435 !----------------------------------------------------------------------
5436 return
5437 end subroutine albwindow
5438 !====================================================================72
5439 !====================================================================72
5441 subroutine foncs(fonc,x,aa,bb,cc,alf,delt,gam)
5443 implicit none
5444 !
5445 real x,aa,bb,cc
5446 real alf,delt,gam
5447 real fonc
5449 fonc=(((aa*(x**alf))/(pi**alf))+ &
5450 ((bb*(x**delt))/(pi**delt))+ &
5451 ((cc*(x**gam))/(pi**gam)))*sin(x)
5453 return
5454 end subroutine foncs
5455 !====================================================================72
5456 !====================================================================72
5458 subroutine icBEP_XY(iurb,fww_u,fwg_u,fgw_u,fsw_u, &
5459 fws_u,fsg_u,ndu,strd,ws,nzu,z_u)
5461 implicit none
5463 ! Street parameters
5464 integer ndu ! Number of street direction for each urban class
5465 integer iurb
5467 real strd(ndm) ! Street length (fix to greater value to the horizontal length of the cells)
5468 real ws(ndm) ! Street width [m]
5470 ! Grid parameters
5471 integer nzu ! Number of layer in the urban grid
5472 real z_u(nz_um) ! Height of the urban grid levels
5473 ! -----------------------------------------------------------------------
5474 ! Output
5475 !------------------------------------------------------------------------
5477 ! fww_u,fwg_u,fgw_u,fsw_u,fsg_u are the view factors used to compute the long wave
5478 ! and the short wave radation. They are the part of radiation from a surface
5479 ! or from the sky to another surface.
5481 real fww_u(nz_um,nz_um,ndm,nurbm) ! from wall to wall
5482 real fwg_u(nz_um,ndm,nurbm) ! from wall to ground
5483 real fgw_u(nz_um,ndm,nurbm) ! from ground to wall
5484 real fsw_u(nz_um,ndm,nurbm) ! from sky to wall
5485 real fws_u(nz_um,ndm,nurbm) ! from sky to wall
5486 real fsg_u(ndm,nurbm) ! from sky to ground
5488 ! -----------------------------------------------------------------------
5489 ! Local
5490 !------------------------------------------------------------------------
5492 integer id
5494 ! -----------------------------------------------------------------------
5495 ! This routine compute the view factors
5496 !------------------------------------------------------------------------
5497 !
5498 !Initialize
5499 !
5500 fww_u=0.
5501 fwg_u=0.
5502 fgw_u=0.
5503 fsw_u=0.
5504 fws_u=0.
5505 fsg_u=0.
5507 do id=1,ndu
5509 call view_factors(iurb,nzu,id,strd(id),z_u,ws(id), &
5510 fww_u,fwg_u,fgw_u,fsg_u,fsw_u,fws_u)
5512 enddo
5513 return
5514 end subroutine icBEP_XY
5515 !====================================================================72
5516 !====================================================================72
5517 subroutine icBEPHI_XY(iurb,hb_u,hi_urb1D,ss_u,pb_u,nzu,z_u)
5519 implicit none
5520 !-----------------------------------------------------------------------
5521 ! Inputs
5522 !-----------------------------------------------------------------------
5523 ! Street parameters
5524 !
5525 real hi_urb1D(nz_um) ! The probability that a building has an height h_b
5526 integer iurb ! Number of the urban class
5527 !
5528 ! Grid parameters
5529 !
5530 real z_u(nz_um) ! Height of the urban grid levels
5531 ! -----------------------------------------------------------------------
5532 ! Output
5533 !------------------------------------------------------------------------
5535 real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to z
5536 real pb_u(nz_um) ! The probability that a building has an height greater or equal to z
5537 !
5538 ! Grid parameters
5539 !
5540 integer nzu ! Number of layer in the urban grid
5542 ! -----------------------------------------------------------------------
5543 ! Local
5544 !------------------------------------------------------------------------
5545 real hb_u(nz_um) ! Bulding's heights [m]
5546 integer iz_u,id,ilu
5548 real dtot
5549 real hbmax
5551 !------------------------------------------------------------------------
5553 !Initialize variables
5554 !
5556 nzu=0
5557 ss_u=0.
5558 pb_u=0.
5560 ! Normalisation of the building density
5562 dtot=0.
5563 hb_u=0.
5565 do ilu=1,nz_um
5566 dtot=dtot+hi_urb1D(ilu)
5567 enddo
5569 do ilu=1,nz_um
5570 if (hi_urb1D(ilu)<0.) then
5571 ! write(*,*) 'WARNING, HI_URB1D(ilu) < 0 IN BEP_BEM'
5572 go to 20
5573 endif
5574 enddo
5576 if (dtot.gt.0.) then
5577 continue
5578 else
5579 ! write(*,*) 'WARNING, HI_URB1D <= 0 IN BEP_BEM'
5580 go to 20
5581 endif
5583 do ilu=1,nz_um
5584 hi_urb1D(ilu)=hi_urb1D(ilu)/dtot
5585 enddo
5587 hb_u(1)=dz_u
5588 do ilu=2,nz_um
5589 hb_u(ilu)=dz_u+hb_u(ilu-1)
5590 enddo
5593 ! Compute pb and ss
5596 hbmax=0.
5598 do ilu=1,nz_um
5599 if (hi_urb1D(ilu)>0.and.hi_urb1D(ilu)<=1.) then
5600 hbmax=hb_u(ilu)
5601 endif
5602 enddo
5604 do iz_u=1,nz_um-1
5605 if(z_u(iz_u+1).gt.hbmax)go to 10
5606 enddo
5608 10 continue
5610 nzu=iz_u+1
5612 if ((nzu+1).gt.nz_um) then
5613 write(*,*) 'error, nz_um has to be increased to at least',nzu+1
5614 stop
5615 endif
5617 do iz_u=1,nzu
5618 ss_u(iz_u,iurb)=0.
5619 do ilu=1,nz_um
5620 if(z_u(iz_u).le.hb_u(ilu) &
5621 .and.z_u(iz_u+1).gt.hb_u(ilu))then
5622 ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+hi_urb1D(ilu)
5623 endif
5624 enddo
5625 enddo
5627 pb_u(1)=1.
5628 do iz_u=1,nzu
5629 pb_u(iz_u+1)=max(0.,pb_u(iz_u)-ss_u(iz_u,iurb))
5630 enddo
5632 20 continue
5633 return
5634 end subroutine icBEPHI_XY
5635 !====================================================================72
5636 !====================================================================72
5637 END MODULE module_sf_bep_bem
5639 ! ===6=8===============================================================72
5640 ! ===6=8===============================================================72
5642 FUNCTION bep_bem_nurbm () RESULT (bep_bem_val_nurbm)
5643 USE module_sf_bep_bem
5644 IMPLICIT NONE
5645 INTEGER :: bep_bem_val_nurbm
5646 bep_bem_val_nurbm = nurbm
5647 END FUNCTION bep_bem_nurbm
5649 FUNCTION bep_bem_ndm () RESULT (bep_bem_val_ndm)
5650 USE module_sf_bep_bem
5651 IMPLICIT NONE
5652 INTEGER :: bep_bem_val_ndm
5653 bep_bem_val_ndm = ndm
5654 END FUNCTION bep_bem_ndm
5656 FUNCTION bep_bem_nz_um () RESULT (bep_bem_val_nz_um)
5657 USE module_sf_bep_bem
5658 IMPLICIT NONE
5659 INTEGER :: bep_bem_val_nz_um
5660 bep_bem_val_nz_um = nz_um
5661 END FUNCTION bep_bem_nz_um
5663 FUNCTION bep_bem_ng_u () RESULT (bep_bem_val_ng_u)
5664 USE module_sf_bep_bem
5665 IMPLICIT NONE
5666 INTEGER :: bep_bem_val_ng_u
5667 bep_bem_val_ng_u = ng_u
5668 END FUNCTION bep_bem_ng_u
5670 FUNCTION bep_bem_nwr_u () RESULT (bep_bem_val_nwr_u)
5671 USE module_sf_bep_bem
5672 IMPLICIT NONE
5673 INTEGER :: bep_bem_val_nwr_u
5674 bep_bem_val_nwr_u = nwr_u
5675 END FUNCTION bep_bem_nwr_u
5677 FUNCTION bep_bem_nf_u () RESULT (bep_bem_val_nf_u)
5678 USE module_sf_bep_bem
5679 IMPLICIT NONE
5680 INTEGER :: bep_bem_val_nf_u
5681 bep_bem_val_nf_u = nf_u
5682 END FUNCTION bep_bem_nf_u
5685 FUNCTION bep_bem_ngb_u () RESULT (bep_bem_val_ngb_u)
5686 USE module_sf_bep_bem
5687 IMPLICIT NONE
5688 INTEGER :: bep_bem_val_ngb_u
5689 bep_bem_val_ngb_u = ngb_u
5690 END FUNCTION bep_bem_ngb_u
5692 FUNCTION bep_bem_nbui_max () RESULT (bep_bem_val_nbui_max)
5693 USE module_sf_bep_bem
5694 IMPLICIT NONE
5695 INTEGER :: bep_bem_val_nbui_max
5696 bep_bem_val_nbui_max = nbui_max
5697 END FUNCTION bep_bem_nbui_max
5700 FUNCTION bep_bem_ngr_u () RESULT (bep_bem_val_ngr_u)
5701 USE module_sf_bep_bem
5702 IMPLICIT NONE
5703 INTEGER :: bep_bem_val_ngr_u
5704 bep_bem_val_ngr_u = ngr_u
5705 END FUNCTION bep_bem_ngr_u