Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_bep_bem.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 159.1 KB

Line
1	MODULE module_sf_bep_bem
2
3	!USE module_model_constants
4	USE module_sf_urban
5	USE module_sf_bem
6
7	! SGClarke 09/11/2008
8	! Access urban_param.tbl values through calling urban_param_init in module_physics_init
9	! for CASE (BEPSCHEME) select sf_urban_physics
10	!
11	! -----------------------------------------------------------------------
12	! Dimension for the array used in the BEP module
13	! -----------------------------------------------------------------------
14
15	integer nurbm ! Maximum number of urban classes
16	parameter (nurbm=3)
17
18	integer ndm ! Maximum number of street directions
19	parameter (ndm=2)
20
21	integer nz_um ! Maximum number of vertical levels in the urban grid
22	parameter(nz_um=13)
23
24	integer ng_u ! Number of grid levels in the ground
25	parameter (ng_u=10)
26	integer nwr_u ! Number of grid levels in the walls or roofs
27	parameter (nwr_u=10)
28
29	integer nf_u !Number of grid levels in the floors (BEM)
30	parameter (nf_u=10)
31
32	integer ngb_u !Number of grid levels in the ground below building (BEM)
33	parameter (ngb_u=10)
34
35	real dz_u ! Urban grid resolution
36	parameter (dz_u=5.)
37
38	integer nbui_max !maximum number of types of buildings in an urban class
39	parameter (nbui_max=4) !must be less or equal than nz_um
40
41	!---------------------------------------------------------------------------------
42	!Parameters of the windows. The glasses of windows are considered without films -
43	!Read the paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour -
44	!of the total solar energy transmittance of windows".Solar Energy Vol.69,No.4, -
45	!pp. 321-329, for more details. -
46	!---------------------------------------------------------------------------------
47	integer p_num !number of panes in the windows (1,2 or 3)
48	parameter (p_num=2)
49	integer q_num !category number for the windows (q_num= 4, standard glasses)
50	parameter(q_num=4) !Possible values 1,2,...,10
51
52
53	! The change of ng_u, nwr_u should be done in agreement with the block data
54	! in the routine "surf_temp"
55	! -----------------------------------------------------------------------
56	! Constant used in the BEP module
57	! -----------------------------------------------------------------------
58
59	real vk ! von Karman constant
60	real g_u ! Gravity acceleration
61	real pi !
62	real r ! Perfect gas constant
63	real cp_u ! Specific heat at constant pressure
64	real rcp_u !
65	real sigma !
66	real p0 ! Reference pressure at the sea level
67	real cdrag ! Drag force constant
68	real latent ! Latent heat of vaporization [J/kg] (used in BEM)
69	parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)
70	parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4,latent=2.45e+06)
71
72	! -----------------------------------------------------------------------
73
74
75
76
77	CONTAINS
78
79	subroutine BEP_BEM(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, &
80	th_phy,rho,p_phy,swdown,glw, &
81	gmt,julday,xlong,xlat, &
82	declin_urb,cosz_urb2d,omg_urb2d, &
83	num_urban_layers, &
84	trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
85	tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
86	tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
87	cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d, &
88	sfwin1_urb3d,sfwin2_urb3d, &
89	sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
90	a_u,a_v,a_t,a_e,b_u,b_v, &
91	b_t,b_e,b_q,dlg,dl_u,sf,vl, &
92	rl_up,rs_abs,emiss,grdflx_urb,qv_phy, &
93	ids,ide, jds,jde, kds,kde, &
94	ims,ime, jms,jme, kms,kme, &
95	its,ite, jts,jte, kts,kte)
96
97	implicit none
98
99	!------------------------------------------------------------------------
100	! Input
101	!------------------------------------------------------------------------
102	INTEGER :: ids,ide, jds,jde, kds,kde, &
103	ims,ime, jms,jme, kms,kme, &
104	its,ite, jts,jte, kts,kte, &
105	itimestep
106
107
108	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: DZ8W
109	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: P_PHY
110	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: RHO
111	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: TH_PHY
112	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: T_PHY
113	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U_PHY
114	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V_PHY
115	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U
116	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V
117	REAL, DIMENSION( ims:ime , jms:jme ) :: GLW
118	REAL, DIMENSION( ims:ime , jms:jme ) :: swdown
119	REAL, DIMENSION( ims:ime, jms:jme ) :: UST
120	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: UTYPE_URB2D
121	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: FRC_URB2D
122	REAL, INTENT(IN ) :: GMT
123	INTEGER, INTENT(IN ) :: JULDAY
124	REAL, DIMENSION( ims:ime, jms:jme ), &
125	INTENT(IN ) :: XLAT, XLONG
126	REAL, INTENT(IN) :: DECLIN_URB
127	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
128	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
129	INTEGER, INTENT(IN ) :: num_urban_layers
130	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
131	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
132	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
133	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
134	!New variables used for BEM
135	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: qv_phy
136	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tlev_urb3d
137	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: qlev_urb3d
138	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
139	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
140	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tglev_urb3d
141	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tflev_urb3d
142	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
143	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
144	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
145	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
146	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
147	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
148	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
149	!End variables
150	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
151	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
152	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
153	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
154
155	real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates
156	REAL, INTENT(IN ):: DT ! Time step
157
158	!
159	!------------------------------------------------------------------------
160	! Output
161	!------------------------------------------------------------------------
162	!
163	! Implicit and explicit components of the source and sink terms at each levels,
164	! the fluxes can be computed as follow: FX = AX + B example: t_fluxes = a_t pt + b_t
165	real a_u(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in X-direction (center)
166	real a_v(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in Y-direction (center)
167	real a_t(ims:ime,kms:kme,jms:jme) ! Implicit component for the temperature
168	real a_e(ims:ime,kms:kme,jms:jme) ! Implicit component for the TKE
169	real b_u(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in X-direction (center)
170	real b_v(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in Y-direction (center)
171	real b_t(ims:ime,kms:kme,jms:jme) ! Explicit component for the temperature
172	real b_e(ims:ime,kms:kme,jms:jme) ! Explicit component for the TKE
173	real b_q(ims:ime,kms:kme,jms:jme) ! Explicit component for the Humidity
174	real dlg(ims:ime,kms:kme,jms:jme) ! Height above ground (L_ground in formula (24) of the BLM paper).
175	real dl_u(ims:ime,kms:kme,jms:jme) ! Length scale (lb in formula (22) ofthe BLM paper).
176	! urban surface and volumes
177	real sf(ims:ime,kms:kme,jms:jme) ! surface of the urban grid cells
178	real vl(ims:ime,kms:kme,jms:jme) ! volume of the urban grid cells
179	! urban fluxes
180	real rl_up(ims:ime,jms:jme) ! upward long wave radiation
181	real rs_abs(ims:ime,jms:jme) ! absorbed short wave radiation
182	real emiss(ims:ime,jms:jme) ! emissivity averaged for urban surfaces
183	real grdflx_urb(ims:ime,jms:jme) ! ground heat flux for urban areas
184	!------------------------------------------------------------------------
185	! Local
186	!------------------------------------------------------------------------
187	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188	! Building parameters
189	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
190	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
191	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
192	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
193	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
194	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
195	real twini_u(nurbm) ! Initial temperature inside the building's wall [K]
196	real trini_u(nurbm) ! Initial temperature inside the building's roof [K]
197	real tgini_u(nurbm) ! Initial road temperature
198	!
199	! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
200	!
201	! Radiation paramters
202	real albg_u(nurbm) ! Albedo of the ground
203	real albw_u(nurbm) ! Albedo of the wall
204	real albr_u(nurbm) ! Albedo of the roof
205	real albwin_u(nurbm) ! Albedo of the windows
206	real emwind_u(nurbm) ! Emissivity of windows
207	real emg_u(nurbm) ! Emissivity of ground
208	real emw_u(nurbm) ! Emissivity of wall
209	real emr_u(nurbm) ! Emissivity of roof
210
211	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
212	! and the short wave radation.
213	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
214	real fwg(nz_um,ndm,nurbm) ! from wall to ground
215	real fgw(nz_um,ndm,nurbm) ! from ground to wall
216	real fsw(nz_um,ndm,nurbm) ! from sky to wall
217	real fws(nz_um,ndm,nurbm) ! from sky to wall
218	real fsg(ndm,nurbm) ! from sky to ground
219
220	! Roughness parameters
221	real z0g_u(nurbm) ! The ground's roughness length
222	real z0r_u(nurbm) ! The roof's roughness length
223
224	! Street parameters
225	integer nd_u(nurbm) ! Number of street direction for each urban class
226	real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
227	real drst_u(ndm,nurbm) ! Street direction
228	real ws_u(ndm,nurbm) ! Street width
229	real bs_u(ndm,nurbm) ! Building width
230	real h_b(nz_um,nurbm) ! Bulding's heights
231	real d_b(nz_um,nurbm) ! Probability that a building has an height h_b
232	real ss_u(nz_um,nurbm) ! Probability that a building has an height equal to z
233	real pb_u(nz_um,nurbm) ! Probability that a building has an height greater or equal to z
234
235
236	! Grid parameters
237	integer nz_u(nurbm) ! Number of layer in the urban grid
238
239	real z_u(nz_um) ! Height of the urban grid levels
240
241	! MT
242	real cop_u(nurbm)
243	real pwin_u(nurbm)
244	real beta_u(nurbm)
245	integer sw_cond_u(nurbm)
246	real time_on_u(nurbm)
247	real time_off_u(nurbm)
248	real targtemp_u(nurbm)
249	real gaptemp_u(nurbm)
250	real targhum_u(nurbm)
251	real gaphum_u(nurbm)
252	real perflo_u(nurbm)
253	real hsesf_u(nurbm)
254	real hsequip(24)
255
256	! 1D array used for the input and output of the routine "urban"
257
258	real z1D(kms:kme) ! vertical coordinates
259	real ua1D(kms:kme) ! wind speed in the x directions
260	real va1D(kms:kme) ! wind speed in the y directions
261	real pt1D(kms:kme) ! potential temperature
262	real da1D(kms:kme) ! air density
263	real pr1D(kms:kme) ! air pressure
264	real pt01D(kms:kme) ! reference potential temperature
265	real zr1D ! zenith angle
266	real deltar1D ! declination of the sun
267	real ah1D ! hour angle (it should come from the radiation routine)
268	real rs1D ! solar radiation
269	real rld1D ! downward flux of the longwave radiation
270
271	real tw1D(2*ndm,nz_um,nwr_u,nbui_max) ! temperature in each layer of the wall
272	real tg1D(ndm,ng_u) ! temperature in each layer of the ground
273	real tr1D(ndm,nz_um,nwr_u) ! temperature in each layer of the roof
274	!
275	!New variable for BEM
276	!
277	real tlev1D(nz_um,nbui_max) ! temperature in each floor and in each different type of building
278	real qlev1D(nz_um,nbui_max) ! specific humidity in each floor and in each different type of building
279	real twlev1D(2*ndm,nz_um,nbui_max) ! temperature in each window in each floor in each different type of building
280	real tglev1D(ndm,ngb_u,nbui_max) ! temperature in each layer of the ground below of a type of building
281	real tflev1D(ndm,nf_u,nz_um-1,nbui_max)! temperature in each layer of the floors in each building
282	real lflev1D(nz_um,nz_um) ! latent heat flux due to the air conditioning systems
283	real sflev1D(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems
284	real lfvlev1D(nz_um,nz_um) ! latent heat flux due to ventilation
285	real sfvlev1D(nz_um,nz_um) ! sensible heat flux due to ventilation
286	real sfwin1D(2*ndm,nz_um,nbui_max) ! sensible heat flux from windows
287	real consumlev1D(nz_um,nz_um) ! consumption due to the air conditioning systems
288	real qv1D(kms:kme) ! specific humidity
289	real meso_urb ! constant to link meso and urban scales [mÃÂ¯2]
290	real d_urb(nz_um)
291	real sf_ac
292	integer ibui,nbui
293	integer nlev(nz_um)
294	!
295	!End new variables
296	!
297
298	real sfw1D(2*ndm,nz_um,nbui_max) ! sensible heat flux from walls
299	real sfg1D(ndm) ! sensible heat flux from ground (road)
300	real sfr1D(ndm,nz_um) ! sensible heat flux from roofs
301	real sf1D(kms:kme) ! surface of the urban grid cells
302	real vl1D(kms:kme) ! volume of the urban grid cells
303	real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
304	real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
305	real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
306	real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
307	real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
308	real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
309	real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
310	real b_ac1D(kms:kme)
311	real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
312	real b_q1D(kms:kme) ! Explicit component of the Humidity sources or sinks
313	real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
314	real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
315	real time_bep
316	! arrays used to collapse indexes
317	integer ind_zwd(nbui_max,nz_um,nwr_u,ndm)
318	integer ind_gd(ng_u,ndm)
319	integer ind_zd(nbui_max,nz_um,ndm)
320	integer ind_zdf(nz_um,ndm)
321	integer ind_zrd(nz_um,nwr_u,ndm)
322	!
323	integer ind_bd(nbui_max,nz_um)
324	integer ind_wd(nbui_max,nz_um,ndm)
325	integer ind_gbd(nbui_max,ngb_u,ndm)
326	integer ind_fbd(nbui_max,nf_u,nz_um-1,ndm)
327	!
328	integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
329	integer it, nint
330	integer iii
331	real tempo
332	logical first
333	character(len=80) :: text
334	data first/.true./
335	save first,time_bep
336
337	save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
338	albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,fww,fwg,fgw,fsw,fws,fsg, &
339	z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
340	nz_u,z_u,albwin_u,emwind_u , &
341	cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, targtemp_u, &
342	gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip
343
344	!------------------------------------------------------------------------
345	! Calculation of the momentum, heat and turbulent kinetic fluxes
346	! produced by builgings
347	!
348	! Reference:
349	! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
350	! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
351	! 261-304
352	!
353	! F. Salamanca and A. Martilli, 2009: 'A new Building Energy Model coupled
354	! with an Urban Canopy Parameterization for urban climate simulations_part II.
355	! Validation with one dimension off-line simulations'. Theor Appl Climatol
356	! DOI 10.1007/s00704-009-0143-8
357	!------------------------------------------------------------------------
358	!prepare the arrays to collapse indexes
359	if(num_urban_layers.lt.nbui_maxnz_umndm*max(nwr_u,ng_u))then
360	write(,)'num_urban_layers too small, please increase to at least ', nbui_maxnz_umndm*max(nwr_u,ng_u)
361	stop
362	endif
363	!
364	!New conditions for BEM
365	!
366	if(num_urban_layers.lt.nbui_max*nz_um)then !limit for indoor temperature and indoor humidity
367	write(,)'num_urban_layers too small, please increase to at least ', nbui_max*nz_um
368	stop
369	endif
370
371	if(num_urban_layers.lt.nbui_maxnz_umndm)then !limit for window temperature
372	write(,)'num_urban_layers too small, please increase to at least ', nbui_maxnz_umndm
373	stop
374	endif
375
376	if(num_urban_layers.lt.nbui_maxndmngb_u)then !limit for ground temperature below a building
377	write(,)'num_urban_layers too small, please increase to at least ', nbui_maxndmngb_u
378	stop
379	endif
380
381	if(num_urban_layers.lt.(nz_um-1)nbui_maxndm*nf_u)then !limit for floor temperature
382	write(,)'num_urban_layers too small, please increase to at least ', nbui_maxndmnf_u*(nz_um-1),num_urban_layers
383	stop
384	endif
385
386	if (ndm.ne.2)then
387	write(,) 'number of directions is not correct',ndm
388	stop
389	endif
390
391	!
392	!End of new conditions
393	!
394	!
395	!Initialize collapse indexes
396	!
397	ind_zwd=0
398	ind_gd=0
399	ind_zd=0
400	ind_zdf=0
401	ind_zrd=0
402	ind_bd=0
403	ind_wd=0
404	ind_gbd=0
405	ind_fbd=0
406	!
407	!End initialization indexes
408	!
409
410	iii=0
411	do ibui=1,nbui_max
412	do iz_u=1,nz_um
413	do iw=1,nwr_u
414	do id=1,ndm
415	iii=iii+1
416	ind_zwd(ibui,iz_u,iw,id)=iii
417	enddo
418	enddo
419	enddo
420	enddo
421
422	iii=0
423	do ig=1,ng_u
424	do id=1,ndm
425	iii=iii+1
426	ind_gd(ig,id)=iii
427	enddo
428	enddo
429
430	iii=0
431	do ibui=1,nbui_max
432	do iz_u=1,nz_um
433	do id=1,ndm
434	iii=iii+1
435	ind_zd(ibui,iz_u,id)=iii
436	enddo
437	enddo
438	enddo
439
440	iii=0
441	do iz_u=1,nz_um
442	do iw=1,nwr_u
443	do id=1,ndm
444	iii=iii+1
445	ind_zrd(iz_u,iw,id)=iii
446	enddo
447	enddo
448	enddo
449
450	!
451	!New indexes for BEM
452	!
453	iii=0
454	do iz_u=1,nz_um
455	do id=1,ndm
456	iii=iii+1
457	ind_zdf(iz_u,id)=iii
458	enddo ! id
459	enddo ! iz_u
460
461	iii=0
462	do ibui=1,nbui_max !Type of building
463	do iz_u=1,nz_um !vertical levels
464	iii=iii+1
465	ind_bd(ibui,iz_u)=iii
466	enddo !iz_u
467	enddo !ibui
468
469
470	iii=0
471	do ibui=1,nbui_max !type of building
472	do iz_u=1,nz_um !vertical levels
473	do id=1,ndm !direction
474	iii=iii+1
475	ind_wd(ibui,iz_u,id)=iii
476	enddo !id
477	enddo !iz_u
478	enddo !ibui
479
480	iii=0
481	do ibui=1,nbui_max!type of building
482	do iw=1,ngb_u !layers in the wall (ground below a building)
483	do id=1,ndm !direction
484	iii=iii+1
485	ind_gbd(ibui,iw,id)=iii
486	enddo !id
487	enddo !iw
488	enddo !ibui
489
490	iii=0
491	do ibui=1,nbui_max !type of building
492	do iw=1,nf_u !layers in the wall (floor)
493	do iz_u=1,nz_um-1 !vertical levels
494	do id=1,ndm !direction
495	iii=iii+1
496	ind_fbd(ibui,iw,iz_u,id)=iii
497	enddo !id
498	enddo !iz_u
499	enddo !iw
500	enddo !ibui
501
502	!
503	!End of new indexes
504	!
505
506	do ix=its,ite
507	do iy=jts,jte
508	z(ix,kts,iy)=0.
509	do iz=kts+1,kte+1
510	z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
511	enddo
512	enddo
513	enddo
514
515	if (first) then ! True only on first call
516	call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
517	twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
518	emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,&
519	cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
520	targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)
521
522	!Initialisation of the urban parameters and calculation of the view factor
523	call icBEP(fww,fwg,fgw,fsw,fws,fsg, &
524	z0g_u,z0r_u, &
525	nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
526	nz_u,z_u)
527
528	first=.false.
529
530	endif ! first
531
532	do ix=its,ite
533	do iy=jts,jte
534	if (FRC_URB2D(ix,iy).gt.0.) then ! Calling BEP only for existing urban classes.
535
536	iurb=UTYPE_URB2D(ix,iy)
537	ibui=0
538	nlev=0
539	nbui=0
540	d_urb=0.
541	do iz=1,nz_um
542	if(ss_u(iz,iurb).gt.0) then
543	ibui=ibui+1
544	nlev(ibui)=iz-1
545	d_urb(ibui)=ss_u(iz,iurb)
546	nbui=ibui
547	endif
548	end do !iz
549	if (nbui.gt.nbui_max) then
550	write (,) 'nbui_max must be increased to',nbui
551	stop
552	endif
553
554	do iz= kts,kte
555	ua1D(iz)=u_phy(ix,iz,iy)
556	va1D(iz)=v_phy(ix,iz,iy)
557	pt1D(iz)=th_phy(ix,iz,iy)
558	da1D(iz)=rho(ix,iz,iy)
559	pr1D(iz)=p_phy(ix,iz,iy)
560	! pt01D(iz)=th_phy(ix,iz,iy)
561	pt01D(iz)=300.
562	z1D(iz)=z(ix,iz,iy)
563	qv1D(iz)=qv_phy(ix,iz,iy)
564	a_u1D(iz)=0.
565	a_v1D(iz)=0.
566	a_t1D(iz)=0.
567	a_e1D(iz)=0.
568	b_u1D(iz)=0.
569	b_v1D(iz)=0.
570	b_t1D(iz)=0.
571	b_ac1D(iz)=0.
572	b_e1D(iz)=0.
573	enddo
574	z1D(kte+1)=z(ix,kte+1,iy)
575
576	do id=1,ndm
577	do iz_u=1,nz_um
578	do iw=1,nwr_u
579	do ibui=1,nbui_max
580	tw1D(2*id-1,iz_u,iw,ibui)=tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
581	tw1D(2*id,iz_u,iw,ibui)=tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)
582	enddo
583	enddo
584	enddo
585	enddo
586
587	do id=1,ndm
588	do ig=1,ng_u
589	tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
590	enddo
591	do iz_u=1,nz_um
592	do ir=1,nwr_u
593	tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)
594	enddo
595	enddo
596	enddo
597	!
598	!Initialize variables for BEM
599	!
600	tlev1D=0. !Indoor temperature
601	qlev1D=0. !Indoor humidity
602
603	twlev1D=0. !Window temperature
604	tglev1D=0. !Ground temperature
605	tflev1D=0. !Floor temperature
606
607	sflev1D=0. !Sensible heat flux from the a.c.
608	lflev1D=0. !latent heat flux from the a.c.
609	consumlev1D=0.!consumption of the a.c.
610	sfvlev1D=0. !Sensible heat flux from natural ventilation
611	lfvlev1D=0. !Latent heat flux from natural ventilation
612	sfwin1D=0. !Sensible heat flux from windows
613	sfw1D=0. !Sensible heat flux from walls
614
615	do iz_u=1,nz_um !vertical levels
616	do ibui=1,nbui_max !Type of building
617	tlev1D(iz_u,ibui)= tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)
618	qlev1D(iz_u,ibui)= qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)
619	enddo !ibui
620	enddo !iz_u
621
622	do id=1,ndm !direction
623	do iz_u=1,nz_um !vertical levels
624	do ibui=1,nbui_max !type of building
625	twlev1D(2*id-1,iz_u,ibui)=tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
626	twlev1D(2*id,iz_u,ibui)=tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
627	sfwin1D(2*id-1,iz_u,ibui)=sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
628	sfwin1D(2*id,iz_u,ibui)=sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)
629	enddo !ibui
630	enddo !iz_u
631	enddo !id
632
633	do id=1,ndm !direction
634	do iw=1,ngb_u !layer in the wall
635	do ibui=1,nbui_max !type of building
636	tglev1D(id,iw,ibui)=tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)
637	enddo !ibui
638	enddo !iw
639	enddo !id
640
641	do id=1,ndm !direction
642	do iw=1,nf_u !layer in the walls
643	do iz_u=1,nz_um-1 !verticals levels
644	do ibui=1,nbui_max !type of building
645	tflev1D(id,iw,iz_u,ibui)=tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)
646	enddo !ibui
647	enddo ! iz_u
648	enddo !iw
649	enddo !id
650
651	!
652	!End initialization for BEM
653	!
654
655	do id=1,ndm
656	do iz=1,nz_um
657	do ibui=1,nbui_max !type of building
658	sfw1D(2*id-1,iz,ibui)=sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)
659	sfw1D(2*id,iz,ibui)=sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)
660	enddo
661	enddo
662	enddo
663
664	do id=1,ndm
665	sfg1D(id)=sfg_urb3d(ix,id,iy)
666	enddo
667
668	do id=1,ndm
669	do iz=1,nz_um
670	sfr1D(id,iz)=sfr_urb3d(ix,ind_zdf(iz,id),iy)
671	enddo
672	enddo
673
674	rs1D=swdown(ix,iy)
675	rld1D=glw(ix,iy)
676
677	zr1D=acos(COSZ_URB2D(ix,iy))
678	deltar1D=DECLIN_URB
679	ah1D=OMG_URB2D(ix,iy)
680
681	call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
682	zr1D,deltar1D,ah1D,rs1D,rld1D, &
683	alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
684	albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, &
685	emwind_u,fww,fwg,fgw,fsw,fws,fsg, &
686	z0g_u,z0r_u, &
687	nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
688	nz_u,z_u, &
689	cop_u,pwin_u,beta_u,sw_cond_u,time_on_u, &
690	time_off_u,targtemp_u,gaptemp_u,targhum_u, &
691	gaphum_u, perflo_u, hsesf_u, hsequip, &
692	tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D, &
693	a_u1D,a_v1D,a_t1D,a_e1D, &
694	b_u1D,b_v1D,b_t1D,b_ac1D,b_e1D,b_q1D, &
695	dlg1D,dl_u1D,sf1D,vl1D,rl_up(ix,iy), &
696	rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy), &
697	qv1D,tlev1D,qlev1D,sflev1D,lflev1D,consumlev1D, &
698	sfvlev1D,lfvlev1D,twlev1D,tglev1D,tflev1D,sfwin1D,&
699	ix,iy)
700
701	do id=1,ndm ! direction
702	do iz=1,nz_um !vertical levels
703	do ibui=1,nbui_max !type of building
704	sfw1_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id-1,iz,ibui)
705	sfw2_urb3d(ix,ind_zd(ibui,iz,id),iy)=sfw1D(2*id,iz,ibui)
706	enddo
707	enddo
708	enddo
709	do id=1,ndm
710	sfg_urb3d(ix,id,iy)=sfg1D(id)
711	enddo
712
713	do id=1,ndm
714	do iz=1,nz_um
715	sfr_urb3d(ix,ind_zdf(iz,id),iy)=sfr1D(id,iz)
716	enddo
717	enddo
718
719	do id=1,ndm
720	do iz_u=1,nz_um
721	do iw=1,nwr_u
722	do ibui=1,nbui_max
723	tw1_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw,ibui)
724	tw2_urb4d(ix,ind_zwd(ibui,iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw,ibui)
725	enddo
726	enddo
727	enddo
728	enddo
729
730	do id=1,ndm
731	do ig=1,ng_u
732	tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
733	enddo
734	do iz_u=1,nz_um
735	do ir=1,nwr_u
736	trb_urb4d(ix,ind_zrd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
737	enddo
738	enddo
739	enddo
740	!
741	!Outputs of BEM
742	!
743	do ibui=1,nbui_max !type of building
744	do iz_u=1,nz_um !vertical levels
745	tlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=tlev1D(iz_u,ibui)
746	qlev_urb3d(ix,ind_bd(ibui,iz_u),iy)=qlev1D(iz_u,ibui)
747	enddo !iz_u
748	enddo !ibui
749	do id=1,ndm !direction
750	do iz_u=1,nz_um !vertical levels
751	do ibui=1,nbui_max !type of building
752	tw1lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id-1,iz_u,ibui)
753	tw2lev_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=twlev1D(2*id,iz_u,ibui)
754	sfwin1_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id-1,iz_u,ibui)
755	sfwin2_urb3d(ix,ind_wd(ibui,iz_u,id),iy)=sfwin1D(2*id,iz_u,ibui)
756	enddo !ibui
757	enddo !iz_u
758	enddo !id
759
760	do id=1,ndm !direction
761	do iw=1,ngb_u !layers in the walls
762	do ibui=1,nbui_max !type of building
763	tglev_urb3d(ix,ind_gbd(ibui,iw,id),iy)=tglev1D(id,iw,ibui)
764	enddo !ibui
765	enddo !iw
766	enddo !id
767
768	do id=1,ndm !direction
769	do iw=1,nf_u !layer in the walls
770	do iz_u=1,nz_um-1 !verticals levels
771	do ibui=1,nbui_max !type of building
772	tflev_urb3d(ix,ind_fbd(ibui,iw,iz_u,id),iy)=tflev1D(id,iw,iz_u,ibui)
773	enddo !ibui
774	enddo !iz_u
775	enddo !iw
776	enddo !id
777	sf_ac_urb3d(ix,iy)=0.
778	lf_ac_urb3d(ix,iy)=0.
779	cm_ac_urb3d(ix,iy)=0.
780	sfvent_urb3d(ix,iy)=0.
781	lfvent_urb3d(ix,iy)=0.
782
783	meso_urb=(1./4.)FRC_URB2D(ix,iy)/((bs_u(1,iurb)+ws_u(1,iurb))bs_u(2,iurb))+ &
784	(1./4.)FRC_URB2D(ix,iy)/((bs_u(2,iurb)+ws_u(2,iurb))bs_u(1,iurb))
785
786
787	ibui=0
788	nlev=0
789	nbui=0
790	d_urb=0.
791	do iz=1,nz_um
792	if(ss_u(iz,iurb).gt.0) then
793	ibui=ibui+1
794	nlev(ibui)=iz-1
795	d_urb(ibui)=ss_u(iz,iurb)
796	nbui=ibui
797	endif
798	end do !iz
799
800	do ibui=1,nbui !type of building
801	do iz_u=1,nlev(ibui) !vertical levels
802	sf_ac_urb3d(ix,iy)=sf_ac_urb3d(ix,iy)+meso_urbd_urb(ibui)sflev1D(iz_u,ibui)
803	lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+meso_urbd_urb(ibui)lflev1D(iz_u,ibui)
804	cm_ac_urb3d(ix,iy)=cm_ac_urb3d(ix,iy)+meso_urbd_urb(ibui)consumlev1D(iz_u,ibui)
805	sfvent_urb3d(ix,iy)=sfvent_urb3d(ix,iy)+meso_urbd_urb(ibui)sfvlev1D(iz_u,ibui)
806	lfvent_urb3d(ix,iy)=lfvent_urb3d(ix,iy)+meso_urbd_urb(ibui)lfvlev1D(iz_u,ibui)
807	enddo !iz_u
808	enddo !ibui
809	!
810	!Add the latent heat exchanged throughout the ventilation in the lf_ac_urb3d output variable.
811	!it is only a print variable
812	!
813	! lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)+lfvent_urb3d(ix,iy)
814	!
815
816	lf_ac_urb3d(ix,iy)=lf_ac_urb3d(ix,iy)-lfvent_urb3d(ix,iy)
817
818	!
819	!End outputs of bem
820	!
821
822	sf_ac=0.
823	do iz= kts,kte
824	sf(ix,iz,iy)=sf1D(iz)
825	vl(ix,iz,iy)=vl1D(iz)
826	a_u(ix,iz,iy)=a_u1D(iz)
827	a_v(ix,iz,iy)=a_v1D(iz)
828	a_t(ix,iz,iy)=a_t1D(iz)
829	a_e(ix,iz,iy)=a_e1D(iz)
830	b_u(ix,iz,iy)=b_u1D(iz)
831	b_v(ix,iz,iy)=b_v1D(iz)
832	b_t(ix,iz,iy)=b_t1D(iz)
833	sf_ac=sf_ac+b_ac1D(iz)da1D(iz)cp_udz8w(ix,iz,iy)vl1D(iz)*FRC_URB2D(ix,iy)
834	b_e(ix,iz,iy)=b_e1D(iz)
835	b_q(ix,iz,iy)=b_q1D(iz)
836	dlg(ix,iz,iy)=dlg1D(iz)
837	dl_u(ix,iz,iy)=dl_u1D(iz)
838	enddo
839	sf(ix,kte+1,iy)=sf1D(kte+1)
840
841	endif ! FRC_URB2D
842
843	enddo ! iy
844	enddo ! ix
845
846	time_bep=time_bep+dt
847
848	return
849	end subroutine BEP_BEM
850
851	! ===6=8===============================================================72
852
853	subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, &
854	zr,deltar,ah,rs,rld, &
855	alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
856	albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,emr_u, &
857	emwind_u,fww,fwg,fgw,fsw,fws,fsg, &
858	z0g_u,z0r_u, &
859	nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
860	nz_u,z_u, &
861	cop_u,pwin_u,beta_u,sw_cond_u,time_on_u, &
862	time_off_u,targtemp_u,gaptemp_u,targhum_u, &
863	gaphum_u, perflo_u, hsesf_u, hsequip, &
864	tw,tg,tr,sfw,sfg,sfr, &
865	a_u,a_v,a_t,a_e, &
866	b_u,b_v,b_t,b_ac,b_e,b_q, &
867	dlg,dl_u,sf,vl,rl_up,rs_abs,emiss,grdflx_urb, &
868	qv,tlev,qlev,sflev,lflev,consumlev, &
869	sfvlev,lfvlev,twlev,tglev,tflev,sfwin,ix,iy)
870
871	! ----------------------------------------------------------------------
872	! This routine computes the effects of buildings on momentum, heat and
873	! TKE (turbulent kinetic energy) sources or sinks and on the mixing length.
874	! It provides momentum, heat and TKE sources or sinks at different levels of a
875	! mesoscale grid defined by the altitude of its cell interfaces "z" and
876	! its number of levels "nz".
877	! The meteorological input parameters (wind, temperature, solar radiation)
878	! are specified on the "mesoscale grid".
879	! The inputs concerning the building and street charateristics are defined
880	! on a "urban grid". The "urban grid" is defined with its number of levels
881	! "nz_u" and its space step "dz_u".
882	! The input parameters are interpolated on the "urban grid". The sources or sinks
883	! are calculated on the "urban grid". Finally the sources or sinks are
884	! interpolated on the "mesoscale grid".
885
886
887	! Mesoscale grid Urban grid Mesoscale grid
888	!
889	! z(4) --- ---
890	! \| \|
891	! \| \|
892	! \| Interpolation Interpolation \|
893	! \| Sources or sinks calculation \|
894	! z(3) --- ---
895	! \| ua ua_u --- uv_a a_u \|
896	! \| va va_u \| uv_b b_u \|
897	! \| pt pt_u --- uh_b a_v \|
898	! z(2) --- \| etc... etc...---
899	! \| z_u(1) --- \|
900	! \| \| \|
901	! z(1) ------------------------------------------------------------
902
903	!
904	! Reference:
905	! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
906	! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
907	! 261-304
908
909	! ----------------------------------------------------------------------
910
911	implicit none
912
913	! ----------------------------------------------------------------------
914	! INPUT:
915	! ----------------------------------------------------------------------
916
917	! Data relative to the "mesoscale grid"
918
919	integer kms,kme,kts,kte
920	real z(kms:kme) ! Altitude above the ground of the cell interfaces.
921	real ua(kms:kme) ! Wind speed in the x direction
922	real va(kms:kme) ! Wind speed in the y direction
923	real pt(kms:kme) ! Potential temperature
924	real da(kms:kme) ! Air density
925	real pr(kms:kme) ! Air pressure
926	real pt0(kms:kme) ! Reference potential temperature (could be equal to "pt")
927	real qv(kms:kme) ! Specific humidity
928	real dt ! Time step
929	real zr ! Zenith angle
930	real deltar ! Declination of the sun
931	real ah ! Hour angle
932	real rs ! Solar radiation
933	real rld ! Downward flux of the longwave radiation
934
935	! Data relative to the "urban grid"
936
937	integer iurb ! Current urban class
938
939	! Building parameters
940	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
941	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
942	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
943	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
944	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
945	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
946
947	! Radiation parameters
948	real albg_u(nurbm) ! Albedo of the ground
949	real albw_u(nurbm) ! Albedo of the wall
950	real albr_u(nurbm) ! Albedo of the roof
951	real albwin_u(nurbm) ! Albedo of the windows
952	real emwind_u(nurbm) ! Emissivity of windows
953	real emg_u(nurbm) ! Emissivity of ground
954	real emw_u(nurbm) ! Emissivity of wall
955	real emr_u(nurbm) ! Emissivity of roof
956
957	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and
958	! short wave radation.
959	! The calculation of these factor is explained in the Appendix A of the BLM paper
960	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
961	real fwg(nz_um,ndm,nurbm) ! from wall to ground
962	real fgw(nz_um,ndm,nurbm) ! from ground to wall
963	real fsw(nz_um,ndm,nurbm) ! from sky to wall
964	real fws(nz_um,ndm,nurbm) ! from wall to sky
965	real fsg(ndm,nurbm) ! from sky to ground
966
967	! Roughness parameters
968	real z0g_u(nurbm) ! The ground's roughness length
969	real z0r_u(nurbm) ! The roof's roughness length
970
971	! Street parameters
972	integer nd_u(nurbm) ! Number of street direction for each urban class
973	real strd_u(ndm,nurbm) ! Street length (set to a greater value then the horizontal length of the cells)
974	real drst_u(ndm,nurbm) ! Street direction
975	real ws_u(ndm,nurbm) ! Street width
976	real bs_u(ndm,nurbm) ! Building width
977	real h_b(nz_um,nurbm) ! Bulding's heights
978	real d_b(nz_um,nurbm) ! The probability that a building has an height "h_b"
979	real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
980	real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
981
982	! Grid parameters
983	integer nz_u(nurbm) ! Number of layer in the urban grid
984	real z_u(nz_um) ! Height of the urban grid levels
985
986	! MT
987	real cop_u(nurbm)
988	real pwin_u(nurbm)
989	real beta_u(nurbm)
990	integer sw_cond_u(nurbm)
991	real time_on_u(nurbm)
992	real time_off_u(nurbm)
993	real targtemp_u(nurbm)
994	real gaptemp_u(nurbm)
995	real targhum_u(nurbm)
996	real gaphum_u(nurbm)
997	real perflo_u(nurbm)
998	real hsesf_u(nurbm)
999	real hsequip(24)
1000
1001	! ----------------------------------------------------------------------
1002	! INPUT-OUTPUT
1003	! ----------------------------------------------------------------------
1004
1005	! Data relative to the "urban grid" which should be stored from the current time step to the next one
1006
1007	real tw(2*ndm,nz_um,nwr_u,nbui_max) ! Temperature in each layer of the wall [K]
1008	real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
1009	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1010	real sfw(2*ndm,nz_um,nbui_max) ! Sensible heat flux from walls
1011	real sfg(ndm) ! Sensible heat flux from ground (road)
1012	real sfr(ndm,nz_um) ! Sensible heat flux from roofs
1013	real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) towards the interior
1014	real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof towards the interior
1015	real gfw(2*ndm,nz_um,nbui_max) ! Heat flux transfered from the surface of the walls towards the interior
1016	! ----------------------------------------------------------------------
1017	! OUTPUT:
1018	! ----------------------------------------------------------------------
1019
1020	! Data relative to the "mesoscale grid"
1021
1022	real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
1023	real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
1024
1025	! Implicit and explicit components of the source and sink terms at each levels,
1026	! the fluxes can be computed as follow: FX = AX + B example: Heat fluxes = a_t pt + b_t
1027	real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
1028	real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
1029	real a_t(kms:kme) ! Implicit component of the heat sources or sinks
1030	real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
1031	real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
1032	real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
1033	real b_t(kms:kme) ! Explicit component of the heat sources or sinks
1034	real b_ac(kms:kme)
1035	real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
1036	real b_q(kms:kme) ! Explicit component of the humidity sources or sinks
1037	real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
1038	real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
1039
1040	! ----------------------------------------------------------------------
1041	! LOCAL:
1042	! ----------------------------------------------------------------------
1043
1044	real dz(kms:kme) ! vertical space steps of the "mesoscale grid"
1045
1046	! Data interpolated from the "mesoscale grid" to the "urban grid"
1047
1048	real ua_u(nz_um) ! Wind speed in the x direction
1049	real va_u(nz_um) ! Wind speed in the y direction
1050	real pt_u(nz_um) ! Potential temperature
1051	real da_u(nz_um) ! Air density
1052	real pt0_u(nz_um) ! Reference potential temperature
1053	real pr_u(nz_um) ! Air pressure
1054	real qv_u(nz_um) !Specific humidity
1055
1056	! Data defining the building and street charateristics
1057
1058	real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
1059
1060	real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
1061	real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
1062	real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
1063
1064	real z0(ndm,nz_um) ! Roughness lengths "profiles"
1065	real ws(ndm) ! Street widths of the current urban class
1066	real bs(ndm) ! Building widths of the current urban class
1067	real strd(ndm) ! Street lengths for the current urban class
1068	real drst(ndm) ! Street directions for the current urban class
1069	real ss(nz_um) ! Probability to have a building with height h
1070	real pb(nz_um) ! Probability to have a building with an height equal
1071
1072	! Solar radiation at each level of the "urban grid"
1073
1074	real rsg(ndm) ! Short wave radiation from the ground
1075	real rsw(2*ndm,nz_um) ! Short wave radiation from the walls
1076	real rlg(ndm) ! Long wave radiation from the ground
1077	real rlw(2*ndm,nz_um) ! Long wave radiation from the walls
1078
1079	! Potential temperature of the surfaces at each level of the "urban grid"
1080
1081	real ptg(ndm) ! Ground potential temperatures
1082	real ptr(ndm,nz_um) ! Roof potential temperatures
1083	real ptw(2*ndm,nz_um,nbui_max) ! Walls potential temperatures
1084
1085
1086	! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
1087	! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
1088	! The fluxes can be computed as follow: Fluxes of X = A*X + B
1089	! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
1090
1091	real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
1092	real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
1093	real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
1094	real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
1095	real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
1096	real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
1097	real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1098	real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1099	real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1100	real tvb_ac(2*ndm,nz_um)
1101	real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1102	real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1103	real qhb_u(ndm,nz_um) ! Humidity Horizontal surfaces, B (explicit) term
1104	real qvb_u(2*ndm,nz_um) ! Humidity Vertical surfaces, B (explicit) term
1105	!
1106	real rs_abs ! solar radiation absorbed by urban surfaces
1107	real rl_up ! longwave radiation emitted by urban surface to the atmosphere
1108	real emiss ! mean emissivity of the urban surface
1109	real grdflx_urb ! ground heat flux
1110	real dt_int ! internal time step
1111	integer nt_int ! number of internal time step
1112	integer iz,id, it_int
1113	integer iw,ix,iy
1114
1115	!---------------------------------------
1116	!New variables uses in BEM
1117	!----------------------------------------
1118
1119	real tmp_u(nz_um) !Air Temperature [K]
1120
1121	real dzw(nwr_u) !Layer sizes in the walls
1122	real dzr(nwr_u) !Layer sizes in the roofs
1123	real dzf(nf_u) !Layer sizes in the floors
1124	real dzgb(ngb_u) !Layer sizes in the ground below the buildings
1125
1126	real csgb(ngb_u) !Specific heat of the ground material below the buildings
1127	!of the current urban class at each ground levels[J m^3 K^-1]
1128	real csf(nf_u) !Specific heat of the floors materials in the buildings
1129	!of the current urban class at each levels[J m^3 K^-1]
1130	real alar(nwr_u+1) ! Roof thermal diffusivity for the current urban class [W/m K]
1131	real alaw(nwr_u+1) ! Walls thermal diffusivity for the current urban class [W/m K]
1132	real alaf(nf_u+1) ! Floor thermal diffusivity at each wall layers [W/m K]
1133	real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall layer [W/m K]
1134
1135	real sfrb(ndm,nbui_max) ! Sensible heat flux from roofs [W/mÃÂ²]
1136	real gfrb(ndm,nbui_max) ! Heat flux flowing inside the roofs [W/mÃÂ²]
1137	real sfwb1D(2*ndm,nz_um) !Sensible heat flux from the walls [W/mÃÂ²]
1138	real sfwin(2*ndm,nz_um,nbui_max)!Sensible heat flux from windows [W/mÃÂ²]
1139	real sfwinb1D(2*ndm,nz_um) !Sensible heat flux from windows [W/mÃÂ²]
1140	real gfwb1D(2*ndm,nz_um) !Heat flux flowing inside the walls [W/mÃÂ²]
1141
1142	real qlev(nz_um,nbui_max) !specific humidity [kg/kg]
1143	real qlevb1D(nz_um) !specific humidity [kg/kg]
1144	real tlev(nz_um,nbui_max) !Indoor temperature [K]
1145	real tlevb1D(nz_um) !Indoor temperature [K]
1146	real twb1D(2*ndm,nwr_u,nz_um) !Wall temperature in BEM [K]
1147	real twlev(2*ndm,nz_um,nbui_max) !Window temperature in BEM [K]
1148	real twlevb1D(2*ndm,nz_um) !Window temperature in BEM [K]
1149	real tglev(ndm,ngb_u,nbui_max) !Ground temperature below a building in BEM [K]
1150	real tglevb1D(ngb_u) !Ground temperature below a building in BEM [K]
1151	real tflev(ndm,nf_u,nz_um-1,nbui_max)!Floor temperature in BEM[K]
1152	real tflevb1D(nf_u,nz_um-1) !Floor temperature in BEM[K]
1153	real trb(ndm,nwr_u,nbui_max) !Roof temperature in BEM [K]
1154	real trb1D(nwr_u) !Roof temperature in BEM [K]
1155
1156	real sflev(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems [W]
1157	real lflev(nz_um,nz_um) ! latent heat flux due to the air conditioning systems [W]
1158	real consumlev(nz_um,nz_um) ! consumption due to the air conditioning systems [W]
1159	real sflev1D(nz_um) ! sensible heat flux due to the air conditioning systems [W]
1160	real lflev1D(nz_um) ! latent heat flux due to the air conditioning systems [W]
1161	real consumlev1D(nz_um) ! consumption due to the air conditioning systems [W]
1162	real sfvlev(nz_um,nz_um) ! sensible heat flux due to ventilation [W]
1163	real lfvlev(nz_um,nz_um) ! latent heat flux due to ventilation [W]
1164	real sfvlev1D(nz_um) ! sensible heat flux due to ventilation [W]
1165	real lfvlev1D(nz_um) ! Latent heat flux due to ventilation [W]
1166
1167	real ptwin(2*ndm,nz_um,nbui_max) ! window potential temperature
1168	real tw_av(2*ndm,nz_um) ! Averaged temperature of the wall surfaces
1169	real twlev_av(2*ndm,nz_um) ! Averaged temperature of the windows
1170	real sfw_av(2*ndm,nz_um) ! Averaged sensible heat from walls
1171	real sfwind_av(2*ndm,nz_um) ! Averaged sensible heat from windows
1172
1173	integer nbui !Total number of different type of buildings in an urban class
1174	integer nlev(nz_um) !Number of levels in each different type of buildings in an urban class
1175	integer ibui,ily
1176	real :: nhourday ! Number of hours from midnight, local time
1177	! ----------------------------------------------------------------------
1178	! END VARIABLES DEFINITIONS
1179	! ----------------------------------------------------------------------
1180
1181	! Fix some usefull parameters for the computation of the sources or sinks
1182	!
1183	!initialize inside param
1184	!
1185	! ss=0.
1186	! pb=0.
1187
1188	do iz=kts,kte
1189	dz(iz)=z(iz+1)-z(iz)
1190	end do
1191	call param(iurb,nz_u(iurb),nd_u(iurb), &
1192	csg_u,csg,alag_u,alag,csr_u,csr, &
1193	alar_u,alar,csw_u,csw,alaw_u,alaw, &
1194	ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
1195	strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb, &
1196	dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)
1197
1198	! Interpolation on the "urban grid"
1199
1200	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,ua,ua_u)
1201	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,va,va_u)
1202	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt,pt_u)
1203	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt0,pt0_u)
1204	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pr,pr_u)
1205	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,da,da_u)
1206	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,qv,qv_u)
1207
1208	! Compute the modification of the radiation due to the buildings
1209
1210	call averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av, &
1211	sfw_av,sfwind_av,sfw,sfwin)
1212	call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws, &
1213	drst,strd,ss,pb, &
1214	tw_av,tg,twlev_av,albg_u(iurb),albw_u(iurb), &
1215	emw_u(iurb),emg_u(iurb),pwin_u(iurb),albwin_u(iurb), &
1216	emwind_u(iurb),fww,fwg,fgw,fsw,fsg, &
1217	zr,deltar,ah, &
1218	rs,rld,rsw,rsg,rlw,rlg)
1219
1220	! calculation of the urban albedo and the upward long wave radiation
1221
1222	call upward_rad(nd_u(iurb),nz_u(iurb),ws,bs,sigma,pb,ss, &
1223	tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg, &
1224	tw_av,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw_av, &
1225	tr,emr_u(iurb),albr_u(iurb),emwind_u(iurb), &
1226	albwin_u(iurb),twlev_av,pwin_u(iurb),sfwind_av,rld,rs,sfr, &
1227	rs_abs,rl_up,emiss,grdflx_urb)
1228
1229	! Compute the surface temperatures
1230
1231	call surf_temp(nd_u(iurb),pr_u,dt, &
1232	rld,rsg,rlg, &
1233	tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg)
1234
1235
1236	! Call the BEM (Building Energy Model) routine
1237
1238	do iz=1,nz_um !Compute the outdoor temperature
1239	tmp_u(iz)=pt_u(iz)(pr_u(iz)/p0)*(rcp_u)
1240	end do
1241
1242	ibui=0
1243	nlev=0
1244	nbui=0
1245
1246	sfrb=0. !Sensible heat flux from roof
1247	gfrb=0. !Heat flux flowing inside the roof
1248	sfwb1D=0. !Sensible heat flux from walls
1249	sfwinb1D=0. !Sensible heat flux from windows
1250	gfwb1D=0. !Heat flux flowing inside the walls[W/mÃÂ²]
1251
1252
1253	twb1D=0. !Wall temperature
1254	twlevb1D=0. !Window temperature
1255	tglevb1D=0. !Ground temperature below a building
1256	tflevb1D=0. !Floor temperature
1257	trb=0. !Roof temperature
1258	trb1D=0. !Roof temperature
1259
1260	qlevb1D=0. !Indoor humidity
1261	tlevb1D=0. !indoor temperature
1262
1263	sflev1D=0. !Sensible heat flux from the a.c.
1264	lflev1D=0. !Latent heat flux from the a.c.
1265	consumlev1D=0.!Consumption from the a.c.
1266	sfvlev1D=0. !Sensible heat flux from the natural ventilation
1267	lfvlev1D=0. !Latent heat flux from natural ventilation
1268
1269	ptw=0. !Wall potential temperature
1270	ptwin=0. !Window potential temperature
1271	ptr=0. !Roof potential temperature
1272
1273	do iz=1,nz_um
1274	if(ss(iz).gt.0) then
1275	ibui=ibui+1
1276	nlev(ibui)=iz-1
1277	nbui=ibui
1278	do id=1,ndm
1279	sfrb(id,ibui)=sfr(id,iz)
1280	do ily=1,nwr_u
1281	trb(id,ily,ibui)=tr(id,iz,ily)
1282	enddo
1283	enddo
1284	endif
1285	end do !iz
1286	!--------------------------------------------------------------------------------
1287	!Loop over BEM -----------------------------------------------------------------
1288	!--------------------------------------------------------------------------------
1289	!--------------------------------------------------------------------------------
1290
1291	nhourday=ah/PI*180./15.+12.
1292	if (nhourday >= 24) nhourday = nhourday - 24
1293	if (nhourday < 0) nhourday = nhourday + 24
1294
1295	do ibui=1,nbui
1296
1297	do iz=1,nz_um
1298	qlevb1D(iz)=qlev(iz,ibui)
1299	tlevb1D(iz)=tlev(iz,ibui)
1300	enddo
1301
1302	do id=1,ndm
1303
1304	do ily=1,nwr_u
1305	trb1D(ily)=trb(id,ily,ibui)
1306	enddo
1307	do ily=1,ngb_u
1308	tglevb1D(ily)=tglev(id,ily,ibui)
1309	enddo
1310
1311	do ily=1,nf_u
1312	do iz=1,nz_um-1
1313	tflevb1D(ily,iz)=tflev(id,ily,iz,ibui)
1314	enddo
1315	enddo
1316
1317	do iz=1,nz_um
1318	sfwinb1D(2id-1,iz)=sfwin(2id-1,iz,ibui)
1319	sfwinb1D(2id,iz)=sfwin(2id,iz,ibui)
1320	enddo
1321
1322	do iz=1,nz_um
1323	do ily=1,nwr_u
1324	twb1D(2id-1,ily,iz)=tw(2id-1,iz,ily,ibui)
1325	twb1D(2id,ily,iz)=tw(2id,iz,ily,ibui)
1326	enddo
1327	sfwb1D(2id-1,iz)=sfw(2id-1,iz,ibui)
1328	sfwb1D(2id,iz)=sfw(2id,iz,ibui)
1329	twlevb1D(2id-1,iz)=twlev(2id-1,iz,ibui)
1330	twlevb1D(2id,iz)=twlev(2id,iz,ibui)
1331	enddo
1332	enddo
1333
1334	call BEM(nz_um,nlev(ibui),nhourday,dt,bs_u(1,iurb), &
1335	bs_u(2,iurb),dz_u,nwr_u,nf_u,nwr_u,ngb_u,sfwb1D,gfwb1D, &
1336	sfwinb1D,sfrb(1,ibui),gfrb(1,ibui), &
1337	latent,sigma,albw_u(iurb),albwin_u(iurb),albr_u(iurb), &
1338	emr_u(iurb),emw_u(iurb),emwind_u(iurb),rsw,rlw,r,cp_u, &
1339	da_u,tmp_u,qv_u,pr_u,rs,rld,dzw,csw,alaw,pwin_u(iurb), &
1340	cop_u(iurb),beta_u(iurb),sw_cond_u(iurb),time_on_u(iurb), &
1341	time_off_u(iurb),targtemp_u(iurb),gaptemp_u(iurb), &
1342	targhum_u(iurb),gaphum_u(iurb),perflo_u(iurb),hsesf_u(iurb), &
1343	hsequip, &
1344	dzf,csf,alaf,dzgb,csgb,alagb,dzr,csr, &
1345	alar,tlevb1D,qlevb1D,twb1D,twlevb1D,tflevb1D,tglevb1D, &
1346	trb1D,sflev1D,lflev1D,consumlev1D,sfvlev1D,lfvlev1D)
1347
1348
1349	!
1350	!Temporal modifications
1351	!
1352	sfrb(2,ibui)=sfrb(1,ibui)
1353	gfrb(2,ibui)=gfrb(1,ibui)
1354	!
1355	!End temporal modifications
1356	!
1357	do iz=1,nz_um
1358	qlev(iz,ibui)=qlevb1D(iz)
1359	tlev(iz,ibui)=tlevb1D(iz)
1360	sflev(iz,ibui)=sflev1D(iz)
1361	lflev(iz,ibui)=lflev1D(iz)
1362	consumlev(iz,ibui)=consumlev1D(iz)
1363	sfvlev(iz,ibui)=sfvlev1D(iz)
1364	lfvlev(iz,ibui)=lfvlev1D(iz)
1365	enddo
1366
1367	do id=1,ndm
1368	do ily=1,nwr_u
1369	trb(id,ily,ibui)=trb1D(ily)
1370	enddo
1371	do ily=1,ngb_u
1372	tglev(id,ily,ibui)=tglevb1D(ily)
1373	enddo
1374
1375	do ily=1,nf_u
1376	do iz=1,nz_um-1
1377	tflev(id,ily,iz,ibui)=tflevb1D(ily,iz)
1378	enddo
1379	enddo
1380
1381
1382	do iz=1,nz_um
1383	do ily=1,nwr_u
1384	tw(2id-1,iz,ily,ibui)=twb1D(2id-1,ily,iz)
1385	tw(2id,iz,ily,ibui)=twb1D(2id,ily,iz)
1386	enddo
1387	gfw(2id-1,iz,ibui)=gfwb1D(2id-1,iz)
1388	gfw(2id,iz,ibui)=gfwb1D(2id,iz)
1389	twlev(2id-1,iz,ibui)=twlevb1D(2id-1,iz)
1390	twlev(2id,iz,ibui)=twlevb1D(2id,iz)
1391	enddo
1392	enddo
1393
1394	enddo !ibui
1395
1396	!-----------------------------------------------------------------------------
1397	!End loop over BEM -----------------------------------------------------------
1398	!-----------------------------------------------------------------------------
1399	!-----------------------------------------------------------------------------
1400
1401	ibui=0
1402	do iz=1,nz_um
1403	if(ss(iz).gt.0) then
1404	ibui=ibui+1
1405	do id=1,ndm
1406	gfr(id,iz)=gfrb(id,ibui)
1407	sfr(id,iz)=sfrb(id,ibui)
1408	do ily=1,nwr_u
1409	tr(id,iz,ily)=trb(id,ily,ibui)
1410	enddo
1411	ptr(id,iz)=tr(id,iz,nwr_u)(pr_u(iz)/p0)*(-rcp_u)
1412	enddo
1413	endif
1414	enddo !iz
1415
1416	!Compute the potential temperature for the vertical surfaces of the buildings
1417
1418	do id=1,ndm
1419	do iz=1,nz_um
1420	do ibui=1,nbui
1421	ptw(2id-1,iz,ibui)=tw(2id-1,iz,nwr_u,ibui)(pr_u(iz)/p0)*(-rcp_u)
1422	ptw(2id,iz,ibui)=tw(2id,iz,nwr_u,ibui)(pr_u(iz)/p0)*(-rcp_u)
1423	ptwin(2id-1,iz,ibui)=twlev(2id-1,iz,ibui)(pr_u(iz)/p0)*(-rcp_u)
1424	ptwin(2id,iz,ibui)=twlev(2id,iz,ibui)(pr_u(iz)/p0)*(-rcp_u)
1425	enddo
1426	enddo
1427	enddo
1428
1429	! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
1430
1431	call buildings(iurb,nd_u(iurb),nz_u(iurb),z0,ua_u,va_u, &
1432	pt_u,pt0_u,ptg,ptr,da_u,ptw,ptwin,pwin_u(iurb),drst, &
1433	uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u,qvb_u,qhb_u, &
1434	uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr, &
1435	sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)
1436
1437	! Calculation of the sensible heat fluxes for the ground, the wall and roof
1438	! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
1439	! where A and B are the implicit and explicit components of the heat sources or sinks.
1440
1441
1442	! Interpolation on the "mesoscale grid"
1443
1444	call urban_meso(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z,dz,z_u,pb,ss,bs,ws,sf, &
1445	vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
1446	uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u, &
1447	a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)
1448
1449
1450	! Calculation of the length scale taking into account the buildings effects
1451
1452	call interp_length(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z_u,z,ss,ws,bs,dlg,dl_u)
1453
1454	return
1455	end subroutine BEP1D
1456
1457	! ===6=8===============================================================72
1458	! ===6=8===============================================================72
1459
1460	subroutine param(iurb,nz,nd, &
1461	csg_u,csg,alag_u,alag,csr_u,csr, &
1462	alar_u,alar,csw_u,csw,alaw_u,alaw, &
1463	ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
1464	strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb, &
1465	dzw,dzr,dzf,csf,alaf,dzgb,csgb,alagb)
1466
1467	! ----------------------------------------------------------------------
1468	! This routine prepare some usefull parameters
1469	! ----------------------------------------------------------------------
1470
1471	implicit none
1472
1473
1474	! ----------------------------------------------------------------------
1475	! INPUT:
1476	! ----------------------------------------------------------------------
1477	integer iurb ! Current urban class
1478	integer nz ! Number of vertical urban levels in the current class
1479	integer nd ! Number of street direction for the current urban class
1480	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
1481	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
1482	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
1483	real bs_u(ndm,nurbm) ! Building width
1484	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
1485	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
1486	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
1487	real drst_u(ndm,nurbm) ! Street direction
1488	real strd_u(ndm,nurbm) ! Street length
1489	real ws_u(ndm,nurbm) ! Street width
1490	real z0g_u(nurbm) ! The ground's roughness length
1491	real z0r_u(nurbm) ! The roof's roughness length
1492	real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
1493	real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
1494
1495
1496	! ----------------------------------------------------------------------
1497	! OUTPUT:
1498	! ----------------------------------------------------------------------
1499	real alag(ng_u) ! Ground thermal diffusivity at each ground levels
1500	real csg(ng_u) ! Specific heat of the ground material at each ground levels
1501	real bs(ndm) ! Building width for the current urban class
1502	real drst(ndm) ! street directions for the current urban class
1503	real strd(ndm) ! Street lengths for the current urban class
1504	real ws(ndm) ! Street widths of the current urban class
1505	real z0(ndm,nz_um) ! Roughness lengths "profiles"
1506	real ss(nz_um) ! Probability to have a building with height h
1507	real pb(nz_um) ! Probability to have a building with an height greater or equal to "z"
1508
1509	!-----------------------------------------------------------------------------
1510	!INPUT/OUTPUT
1511	!-----------------------------------------------------------------------------
1512
1513	real dzw(nwr_u) !Layer sizes in the walls [m]
1514	real dzr(nwr_u) !Layer sizes in the roofs [m]
1515	real dzf(nf_u) !Layer sizes in the floors [m]
1516	real dzgb(ngb_u) !layer sizes in the ground below the buildings [m]
1517
1518	real csr(nwr_u) ! Specific heat of the roof material at each roof levels
1519	real csw(nwr_u) ! Specific heat of the wall material at each wall levels
1520
1521	real csf(nf_u) !Specific heat of the floors materials in the buildings
1522	!of the current urban class [J m^3 K^-1]
1523	real csgb(ngb_u) !Specific heat of the ground material below the buildings
1524	!of the current urban class [J m^3 K^-1]
1525	real alar(nwr_u+1) ! Roof thermal diffusivity at each roof levels [W/ m K]
1526	real alaw(nwr_u+1) ! Wall thermal diffusivity at each wall levels [W/ m K]
1527	real alaf(nf_u+1) ! Floor thermal diffusivity at each wall levels [W/m K]
1528	real alagb(ngb_u+1) ! Ground thermal diffusivity below the building at each wall levels [W/m K]
1529	! ----------------------------------------------------------------------
1530	! LOCAL:
1531	! ----------------------------------------------------------------------
1532	integer id,ig,ir,iw,iz,iflo
1533	! ----------------------------------------------------------------------
1534	! END VARIABLES DEFINITIONS
1535	! ----------------------------------------------------------------------
1536	!Define the layer sizes in the walls
1537
1538	ss=0.
1539	pb=0.
1540	csg=0.
1541	alag=0.
1542	csr=0.
1543	alar=0.
1544	csw=0.
1545	alaw=0.
1546	z0=0.
1547	ws=0.
1548	bs=0.
1549	strd=0.
1550	drst=0.
1551	csgb=0.
1552	alagb=0.
1553	csf=0.
1554	alaf=0.
1555
1556	dzgb=(/0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/)
1557	dzr=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
1558	dzw=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/)
1559	dzf=(/0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02/)
1560
1561	do ig=1,ngb_u
1562	csgb(ig) = csg_u(iurb)
1563	alagb(ig)= csg_u(iurb)*alag_u(iurb)
1564	enddo
1565	alagb(ngb_u+1)= csg_u(iurb)*alag_u(iurb)
1566
1567	do iflo=1,nf_u
1568	csf(iflo) = csw_u(iurb)
1569	alaf(iflo)= csw_u(iurb)*alaw_u(iurb)
1570	enddo
1571	alaf(nf_u+1)= csw_u(iurb)*alaw_u(iurb)
1572
1573	do ir=1,nwr_u
1574	csr(ir) = csr_u(iurb)
1575	alar(ir)= csr_u(iurb)*alar_u(iurb)
1576	enddo
1577	alar(nwr_u+1)= csr_u(iurb)*alar_u(iurb)
1578
1579	do iw=1,nwr_u
1580	csw(iw) = csw_u(iurb)
1581	alaw(iw)= csw_u(iurb)*alaw_u(iurb)
1582	enddo
1583	alaw(nwr_u+1)=csw_u(iurb)*alaw_u(iurb)
1584
1585	!------------------------------------------------------------------------
1586
1587	do iz=1,nz+1
1588	ss(iz)=ss_u(iz,iurb)
1589	pb(iz)=pb_u(iz,iurb)
1590	end do
1591
1592	do ig=1,ng_u
1593	csg(ig)=csg_u(iurb)
1594	alag(ig)=alag_u(iurb)
1595	enddo
1596
1597	do id=1,nd
1598	z0(id,1)=z0g_u(iurb)
1599	do iz=2,nz+1
1600	z0(id,iz)=z0r_u(iurb)
1601	enddo
1602	enddo
1603
1604	do id=1,nd
1605	ws(id)=ws_u(id,iurb)
1606	bs(id)=bs_u(id,iurb)
1607	strd(id)=strd_u(id,iurb)
1608	drst(id)=drst_u(id,iurb)
1609	enddo
1610
1611
1612	return
1613	end subroutine param
1614
1615	! ===6=8===============================================================72
1616	! ===6=8===============================================================72
1617
1618	subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
1619
1620	! ----------------------------------------------------------------------
1621	! This routine interpolate para
1622	! meters from the "mesoscale grid" to
1623	! the "urban grid".
1624	! See p300 Appendix B.1 of the BLM paper.
1625	! ----------------------------------------------------------------------
1626
1627	implicit none
1628
1629	! ----------------------------------------------------------------------
1630	! INPUT:
1631	! ----------------------------------------------------------------------
1632	! Data relative to the "mesoscale grid"
1633	integer kts,kte,kms,kme
1634	real z(kms:kme) ! Altitude of the cell interface
1635	real c(kms:kme) ! Parameter which has to be interpolated
1636	! Data relative to the "urban grid"
1637	integer nz_u ! Number of levels
1638	real z_u(nz_um) ! Altitude of the cell interface
1639	! ----------------------------------------------------------------------
1640	! OUTPUT:
1641	! ----------------------------------------------------------------------
1642	real c_u(nz_um) ! Interpolated paramters in the "urban grid"
1643
1644	! LOCAL:
1645	! ----------------------------------------------------------------------
1646	integer iz_u,iz
1647	real ctot,dz
1648
1649	! ----------------------------------------------------------------------
1650	! END VARIABLES DEFINITIONS
1651	! ----------------------------------------------------------------------
1652
1653	do iz_u=1,nz_u
1654	ctot=0.
1655	do iz=kts,kte
1656	dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
1657	ctot=ctot+c(iz)*dz
1658	enddo
1659	c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
1660	enddo
1661
1662	return
1663	end subroutine interpol
1664
1665	! ===6=8===============================================================72
1666	! ===6=8===============================================================72
1667
1668	subroutine averaging_temp(tw,twlev,ss,pb,tw_av,twlev_av,&
1669	sfw_av,sfwind_av,sfw,sfwin)
1670
1671	implicit none
1672	!
1673	!INPUT VARIABLES
1674	!
1675	real tw(2*ndm,nz_um,nwr_u,nbui_max) ! Temperature in each layer of the wall [K]
1676	real twlev(2*ndm,nz_um,nbui_max) ! Window temperature in BEM [K]
1677	real pb(nz_um) ! Probability to have a building with an height equal or greater h
1678	real ss(nz_um) ! Probability to have a building with height h
1679	real sfw(2*ndm,nz_um,nbui_max) ! Surface fluxes from the walls
1680	real sfwin(2*ndm,nz_um,nbui_max) ! Surface fluxes from the windows
1681	!
1682	!OUTPUT VARIABLES
1683	!
1684	real tw_av(2*ndm,nz_um) ! Averaged temperature of the wall surfaces
1685	real twlev_av(2*ndm,nz_um) ! Averaged temperature of the windows
1686	real sfw_av(2*ndm,nz_um) ! Averaged sensible heat from walls
1687	real sfwind_av(2*ndm,nz_um) ! Averaged sensible heat from windows
1688	!
1689	!LOCAL VARIABLES
1690	!
1691	real d_urb(nz_um)
1692	integer nlev(nz_um)
1693	integer id,iz
1694	integer nbui,ibui
1695	!
1696	!Initialize Variables
1697	!
1698	tw_av=0.
1699	twlev_av=0.
1700	sfw_av=0.
1701	sfwind_av=0.
1702	ibui=0
1703	nbui=0
1704	nlev=0
1705	d_urb=0.
1706
1707	do iz=1,nz_um
1708	if(ss(iz).gt.0) then
1709	ibui=ibui+1
1710	d_urb(ibui)=ss(iz)
1711	nlev(ibui)=iz-1
1712	nbui=ibui
1713	endif
1714	enddo
1715
1716	do id=1,ndm
1717	do iz=1,nz_um-1
1718	if (pb(iz+1).gt.0) then
1719	do ibui=1,nbui
1720	if (iz.le.nlev(ibui)) then
1721	tw_av(2id-1,iz)=tw_av(2id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
1722	tw(2id-1,iz,nwr_u,ibui)*4
1723	tw_av(2id,iz)=tw_av(2id,iz)+(d_urb(ibui)/pb(iz+1))*&
1724	tw(2id,iz,nwr_u,ibui)*4
1725	twlev_av(2id-1,iz)=twlev_av(2id-1,iz)+(d_urb(ibui)/pb(iz+1))*&
1726	twlev(2id-1,iz,ibui)*4
1727	twlev_av(2id,iz)=twlev_av(2id,iz)+(d_urb(ibui)/pb(iz+1))*&
1728	twlev(2id,iz,ibui)*4
1729	sfw_av(2id-1,iz)=sfw_av(2id-1,iz)+(d_urb(ibui)/pb(iz+1))sfw(2id-1,iz,ibui)
1730	sfw_av(2id,iz)=sfw_av(2id,iz)+(d_urb(ibui)/pb(iz+1))sfw(2id,iz,ibui)
1731	sfwind_av(2id-1,iz)=sfwind_av(2id-1,iz)+(d_urb(ibui)/pb(iz+1))sfwin(2id-1,iz,ibui)
1732	sfwind_av(2id,iz)=sfwind_av(2id,iz)+(d_urb(ibui)/pb(iz+1))sfwin(2id,iz,ibui)
1733	endif
1734	enddo
1735	tw_av(2id-1,iz)=tw_av(2id-1,iz)**(1./4.)
1736	tw_av(2id,iz)=tw_av(2id,iz)**(1./4.)
1737	twlev_av(2id-1,iz)=twlev_av(2id-1,iz)**(1./4.)
1738	twlev_av(2id,iz)=twlev_av(2id,iz)**(1./4.)
1739	endif
1740	enddo !iz
1741	enddo !id
1742	return
1743	end subroutine averaging_temp
1744	! ===6=8===============================================================72
1745	! ===6=8===============================================================72
1746
1747	subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, &
1748	tw,tg,twlev,albg,albw,emw,emg,pwin,albwin, &
1749	emwin,fww,fwg,fgw,fsw,fsg, &
1750	zr,deltar,ah, &
1751	rs,rl,rsw,rsg,rlw,rlg)
1752
1753	! ----------------------------------------------------------------------
1754	! This routine computes the modification of the short wave and
1755	! long wave radiation due to the buildings.
1756	! ----------------------------------------------------------------------
1757
1758	implicit none
1759
1760
1761	! ----------------------------------------------------------------------
1762	! INPUT:
1763	! ----------------------------------------------------------------------
1764	integer iurb ! current urban class
1765	integer nd ! Number of street direction for the current urban class
1766	integer nz_u ! Number of layer in the urban grid
1767	real z(nz_um) ! Height of the urban grid levels
1768	real ws(ndm) ! Street widths of the current urban class
1769	real drst(ndm) ! street directions for the current urban class
1770	real strd(ndm) ! Street lengths for the current urban class
1771	real ss(nz_um) ! probability to have a building with height h
1772	real pb(nz_um) ! probability to have a building with an height equal
1773	real tw(2*ndm,nz_um) ! Temperature in each layer of the wall [K]
1774	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1775	real albg ! Albedo of the ground for the current urban class
1776	real albw ! Albedo of the wall for the current urban class
1777	real emg ! Emissivity of ground for the current urban class
1778	real emw ! Emissivity of wall for the current urban class
1779	real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
1780	real fsg(ndm,nurbm) ! View factors from sky to ground
1781	real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
1782	real fws(nz_um,ndm,nurbm) ! View factors from wall to sky
1783	real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
1784	real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
1785	real ah ! Hour angle (it should come from the radiation routine)
1786	real zr ! zenith angle
1787	real deltar ! Declination of the sun
1788	real rs ! solar radiation
1789	real rl ! downward flux of the longwave radiation
1790	!
1791	!New variables BEM
1792	!
1793	real twlev(2*ndm,nz_um) ! Window temperature in BEM [K]
1794	real pwin ! Coverage area fraction of windows in the walls of the buildings
1795	real albwin ! Albedo of the windows for the current urban class
1796	real emwin ! Emissivity of the windows for the current urban class
1797	real alb_av ! Averaged albedo (window and wall)
1798	! ----------------------------------------------------------------------
1799	! OUTPUT:
1800	! ----------------------------------------------------------------------
1801	real rlg(ndm) ! Long wave radiation at the ground
1802	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1803	real rsg(ndm) ! Short wave radiation at the ground
1804	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1805
1806	! ----------------------------------------------------------------------
1807	! LOCAL:
1808	! ----------------------------------------------------------------------
1809
1810	integer id,iz
1811
1812	! Calculation of the shadow effects
1813
1814	call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
1815	rs,rsw,rsg)
1816
1817	! Calculation of the reflection effects
1818	do id=1,nd
1819	call long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev, &
1820	fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
1821
1822	alb_av=pwinalbwin+(1.-pwin)albw
1823
1824	call short_rad(iurb,nz_u,id,alb_av,albg,fwg,fww,fgw,rsg,rsw,pb)
1825
1826	enddo
1827	return
1828	end subroutine modif_rad
1829
1830
1831	! ===6=8===============================================================72
1832	! ===6=8===============================================================72
1833
1834	subroutine surf_temp(nd,pr,dt,rl,rsg,rlg, &
1835	tg,alag,csg,emg,albg,ptg,sfg,gfg)
1836
1837	! ----------------------------------------------------------------------
1838	! Computation of the surface temperatures for walls, ground and roofs
1839	! ----------------------------------------------------------------------
1840
1841	implicit none
1842
1843	! ----------------------------------------------------------------------
1844	! INPUT:
1845	! ----------------------------------------------------------------------
1846
1847	integer nd ! Number of street direction for the current urban class
1848	real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
1849
1850	real albg ! Albedo of the ground for the current urban class
1851
1852	real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
1853
1854	real dt ! Time step
1855	real emg ! Emissivity of ground for the current urban class
1856
1857	real pr(nz_um) ! Air pressure
1858	real rl ! Downward flux of the longwave radiation
1859	real rlg(ndm) ! Long wave radiation at the ground
1860	real rsg(ndm) ! Short wave radiation at the ground
1861	real sfg(ndm) ! Sensible heat flux from ground (road)
1862
1863	real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) toward the interior
1864
1865	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1866
1867
1868
1869	! ----------------------------------------------------------------------
1870	! OUTPUT:
1871	! ----------------------------------------------------------------------
1872	real ptg(ndm) ! Ground potential temperatures
1873
1874
1875	! ----------------------------------------------------------------------
1876	! LOCAL:
1877	! ----------------------------------------------------------------------
1878	integer id,ig,ir,iw,iz
1879
1880	real rtg(ndm) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
1881
1882	real tg_tmp(ng_u)
1883
1884
1885	real dzg_u(ng_u) ! Layer sizes in the ground
1886
1887	data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
1888
1889	! ----------------------------------------------------------------------
1890	! END VARIABLES DEFINITIONS
1891	! ----------------------------------------------------------------------
1892
1893
1894
1895	do id=1,nd
1896
1897	! Calculation for the ground surfaces
1898	do ig=1,ng_u
1899	tg_tmp(ig)=tg(id,ig)
1900	end do
1901	!
1902	call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg, &
1903	rsg(id),rlg(id),pr(1), &
1904	dt,emg,albg, &
1905	rtg(id),sfg(id),gfg(id))
1906	do ig=1,ng_u
1907	tg(id,ig)=tg_tmp(ig)
1908	end do
1909
1910	end do !id
1911
1912	return
1913	end subroutine surf_temp
1914
1915	! ===6=8===============================================================72
1916	! ===6=8===============================================================72
1917
1918	subroutine buildings(iurb,nd,nz,z0,ua_u,va_u,pt_u,pt0_u, &
1919	ptg,ptr,da_u,ptw,ptwin,pwin, &
1920	drst,uva_u,vva_u,uvb_u,vvb_u, &
1921	tva_u,tvb_u,evb_u,qvb_u,qhb_u, &
1922	uhb_u,vhb_u,thb_u,ehb_u,ss,dt,sfw,sfg,sfr, &
1923	sfwin,pb,bs_u,dz_u,sflev,lflev,sfvlev,lfvlev,tvb_ac)
1924
1925	! ----------------------------------------------------------------------
1926	! This routine computes the sources or sinks of the different quantities
1927	! on the urban grid. The actual calculation is done in the subroutines
1928	! called flux_wall, and flux_flat.
1929	! ----------------------------------------------------------------------
1930
1931	implicit none
1932
1933
1934	! ----------------------------------------------------------------------
1935	! INPUT:
1936	! ----------------------------------------------------------------------
1937	integer nd ! Number of street direction for the current urban class
1938	integer nz ! number of vertical space steps
1939	real ua_u(nz_um) ! Wind speed in the x direction on the urban grid
1940	real va_u(nz_um) ! Wind speed in the y direction on the urban grid
1941	real da_u(nz_um) ! air density on the urban grid
1942	real drst(ndm) ! Street directions for the current urban class
1943	real dz
1944	real pt_u(nz_um) ! Potential temperature on the urban grid
1945	real pt0_u(nz_um) ! reference potential temperature on the urban grid
1946	real ptg(ndm) ! Ground potential temperatures
1947	real ptr(ndm,nz_um) ! Roof potential temperatures
1948	real ptw(2*ndm,nz_um,nbui_max) ! Walls potential temperatures
1949	real ss(nz_um) ! probability to have a building with height h
1950	real pb(nz_um)
1951	real z0(ndm,nz_um) ! Roughness lengths "profiles"
1952	real dt ! time step
1953	integer iurb !Urban class
1954
1955	!
1956	!New variables (BEM)
1957	!
1958	real bs_u(ndm,nurbm) ! Building width [m]
1959	real dz_u ! Urban grid resolution
1960	real sflev(nz_um,nz_um) ! sensible heat flux due to the air conditioning systems [W]
1961	real lflev(nz_um,nz_um) ! latent heat flux due to the air conditioning systems [W]
1962	real sfvlev(nz_um,nz_um) ! sensible heat flux due to ventilation [W]
1963	real lfvlev(nz_um,nz_um) ! latent heat flux due to ventilation [W]
1964	real qvb_u(2*ndm,nz_um)
1965	real qhb_u(ndm,nz_um)
1966	real ptwin(2*ndm,nz_um,nbui_max) ! window potential temperature
1967	real pwin
1968	real tvb_ac(2*ndm,nz_um)
1969	! ----------------------------------------------------------------------
1970	! OUTPUT:
1971	! ----------------------------------------------------------------------
1972	! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
1973	! vertical surfaces (walls) and horizontal surfaces (roofs and street)
1974	! The fluxes can be computed as follow: Fluxes of X = A*X + B
1975	! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
1976
1977	real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
1978	real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
1979	real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
1980	real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
1981	real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
1982	real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
1983	real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1984	real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1985	real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1986	real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1987	real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1988	real sfw(2*ndm,nz_um,nbui_max) ! sensible heat flux from walls
1989	real sfwin(2*ndm,nz_um,nbui_max) ! sensible heat flux form windows
1990	real sfr(ndm,nz_um) ! sensible heat flux from roof
1991	real sfg(ndm) ! sensible heat flux from street
1992
1993	! ----------------------------------------------------------------------
1994	! LOCAL:
1995	! ----------------------------------------------------------------------
1996	real d_urb(nz_um)
1997	real uva_tmp
1998	real vva_tmp
1999	real uvb_tmp
2000	real vvb_tmp
2001	real evb_tmp
2002	integer nlev(nz_um)
2003	integer id,iz,ibui,nbui
2004
2005	! ----------------------------------------------------------------------
2006	! END VARIABLES DEFINITIONS
2007	! ----------------------------------------------------------------------
2008	dz=dz_u
2009	ibui=0
2010	nbui=0
2011	nlev=0
2012	d_urb=0.
2013
2014	uva_u=0.
2015	uvb_u=0.
2016	vhb_u=0.
2017	vva_u=0.
2018	vvb_u=0.
2019	thb_u=0.
2020	tva_u=0.
2021	tvb_u=0.
2022	tvb_ac=0.
2023	ehb_u=0.
2024	evb_u=0.
2025	qvb_u=0.
2026	qhb_u=0.
2027
2028	do iz=1,nz_um
2029	if(ss(iz).gt.0) then
2030	ibui=ibui+1
2031	d_urb(ibui)=ss(iz)
2032	nlev(ibui)=iz-1
2033	nbui=ibui
2034	endif
2035	enddo
2036	if (nbui.gt.nbui_max) then
2037	write(,) 'nbui_max must be increased to',nbui
2038	stop
2039	endif
2040	do id=1,nd
2041
2042	! Calculation at the ground surfaces
2043
2044	call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1), &
2045	ptg(id),uhb_u(id,1), &
2046	vhb_u(id,1),sfg(id),ehb_u(id,1),da_u(1))
2047	thb_u(id,1)=- sfg(id)/(da_u(1)*cp_u)
2048
2049	! Calculation at the roof surfaces
2050
2051	do iz=2,nz
2052	if(ss(iz).gt.0)then
2053	call flux_flat(dz,z0(id,iz),ua_u(iz), &
2054	va_u(iz),pt_u(iz),pt0_u(iz), &
2055	ptr(id,iz),uhb_u(id,iz), &
2056	vhb_u(id,iz),sfr(id,iz),ehb_u(id,iz),da_u(iz))
2057	thb_u(id,iz)=- sfr(id,iz)/(da_u(iz)*cp_u)
2058	else
2059	uhb_u(id,iz) = 0.0
2060	vhb_u(id,iz) = 0.0
2061	thb_u(id,iz) = 0.0
2062	ehb_u(id,iz) = 0.0
2063	endif
2064	end do
2065
2066	! Calculation at the wall surfaces
2067
2068	do ibui=1,nbui
2069	do iz=1,nlev(ibui)
2070
2071	call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
2072	ptw(2id-1,iz,ibui),ptwin(2id-1,iz,ibui), &
2073	uva_tmp,vva_tmp, &
2074	uvb_tmp,vvb_tmp, &
2075	sfw(2id-1,iz,ibui),sfwin(2id-1,iz,ibui), &
2076	evb_tmp,drst(id),dt)
2077
2078	if (pb(iz+1).gt.0.) then
2079
2080	uva_u(2id-1,iz)=uva_u(2id-1,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
2081	vva_u(2id-1,iz)=vva_u(2id-1,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
2082	uvb_u(2id-1,iz)=uvb_u(2id-1,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
2083	vvb_u(2id-1,iz)=vvb_u(2id-1,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
2084	evb_u(2id-1,iz)=evb_u(2id-1,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
2085	tvb_u(2id-1,iz)=tvb_u(2id-1,iz)-(d_urb(ibui)/pb(iz+1))* &
2086	(sfw(2id-1,iz,ibui)(1.-pwin)+sfwin(2id-1,iz,ibui)pwin)/ &
2087	da_u(iz)/cp_u-(1./4.)(d_urb(ibui)/pb(iz+1))(sfvlev(iz,ibui)-sflev(iz,ibui))/&
2088	(dz*bs_u(id,iurb))/da_u(iz)/cp_u
2089	tvb_ac(2id-1,iz)=tvb_ac(2id-1,iz)-(1./4.)(d_urb(ibui)/pb(iz+1))(-sflev(iz,ibui))/&
2090	(dz*bs_u(id,iurb))/da_u(iz)/cp_u
2091	qvb_u(2id-1,iz)=qvb_u(2id-1,iz)-(1./4.)(d_urb(ibui)/pb(iz+1))(lfvlev(iz,ibui)-lflev(iz,ibui))/&
2092	(dz*bs_u(id,iurb))/da_u(iz)/latent
2093
2094	endif
2095
2096	call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
2097	ptw(2id,iz,ibui),ptwin(2id,iz,ibui), &
2098	uva_tmp,vva_tmp, &
2099	uvb_tmp,vvb_tmp, &
2100	sfw(2id,iz,ibui),sfwin(2id,iz,ibui), &
2101	evb_tmp,drst(id),dt)
2102
2103	if (pb(iz+1).gt.0.) then
2104
2105	uva_u(2id,iz)=uva_u(2id,iz)+d_urb(ibui)/pb(iz+1)*uva_tmp
2106	vva_u(2id,iz)=vva_u(2id,iz)+d_urb(ibui)/pb(iz+1)*vva_tmp
2107	uvb_u(2id,iz)=uvb_u(2id,iz)+d_urb(ibui)/pb(iz+1)*uvb_tmp
2108	vvb_u(2id,iz)=vvb_u(2id,iz)+d_urb(ibui)/pb(iz+1)*vvb_tmp
2109	evb_u(2id,iz)=evb_u(2id,iz)+d_urb(ibui)/pb(iz+1)*evb_tmp
2110	tvb_u(2id,iz)=tvb_u(2id,iz)-(d_urb(ibui)/pb(iz+1))* &
2111	(sfw(2id,iz,ibui)(1.-pwin)+sfwin(2id,iz,ibui)pwin)/ &
2112	da_u(iz)/cp_u-(1./4.)(d_urb(ibui)/pb(iz+1))(sfvlev(iz,ibui)-sflev(iz,ibui))/&
2113	(dz*bs_u(id,iurb))/da_u(iz)/cp_u
2114	tvb_ac(2id,iz)=tvb_ac(2id,iz)-(1./4.)(d_urb(ibui)/pb(iz+1))(-sflev(iz,ibui))/&
2115	(dz*bs_u(id,iurb))/da_u(iz)/cp_u
2116	qvb_u(2id,iz)=qvb_u(2id,iz)-(1./4.)(d_urb(ibui)/pb(iz+1))(lfvlev(iz,ibui)-lflev(iz,ibui))/&
2117	(dz*bs_u(id,iurb))/da_u(iz)/latent
2118
2119	endif
2120	!
2121	enddo !iz
2122	enddo !ibui
2123
2124	end do !id
2125
2126	return
2127	end subroutine buildings
2128
2129
2130	! ===6=8===============================================================72
2131	! ===6=8===============================================================72
2132
2133	subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl, &
2134	uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
2135	uhb_u,vhb_u,thb_u,ehb_u,qhb_u,qvb_u, &
2136	a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e,b_q,tvb_ac,b_ac)
2137
2138	! ----------------------------------------------------------------------
2139	! This routine interpolates the parameters from the "urban grid" to the
2140	! "mesoscale grid".
2141	! See p300-301 Appendix B.2 of the BLM paper.
2142	! ----------------------------------------------------------------------
2143
2144	implicit none
2145
2146	! ----------------------------------------------------------------------
2147	! INPUT:
2148	! ----------------------------------------------------------------------
2149	! Data relative to the "mesoscale grid"
2150	integer kms,kme,kts,kte
2151	real z(kms:kme) ! Altitude above the ground of the cell interface
2152	real dz(kms:kme) ! Vertical space steps
2153
2154	! Data relative to the "uban grid"
2155	integer nz_u ! Number of layer in the urban grid
2156	integer nd ! Number of street direction for the current urban class
2157	real bs(ndm) ! Building widths of the current urban class
2158	real ws(ndm) ! Street widths of the current urban class
2159	real z_u(nz_um) ! Height of the urban grid levels
2160	real pb(nz_um) ! Probability to have a building with an height equal
2161	real ss(nz_um) ! Probability to have a building with height h
2162	real uhb_u(ndm,nz_um) ! U (x-wind component) Horizontal surfaces, B (explicit) term
2163	real uva_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, A (implicit) term
2164	real uvb_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, B (explicit) term
2165	real vhb_u(ndm,nz_um) ! V (y-wind component) Horizontal surfaces, B (explicit) term
2166	real vva_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, A (implicit) term
2167	real vvb_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, B (explicit) term
2168	real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
2169	real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
2170	real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
2171	real tvb_ac(2*ndm,nz_um)
2172	real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
2173	real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
2174	!
2175	!New variables for BEM
2176	!
2177	real qhb_u(ndm,nz_um)
2178	real qvb_u(2*ndm,nz_um)
2179
2180
2181	! ----------------------------------------------------------------------
2182	! OUTPUT:
2183	! ----------------------------------------------------------------------
2184	! Data relative to the "mesoscale grid"
2185	real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
2186	real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
2187	real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
2188	real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
2189	real a_t(kms:kme) ! Implicit component of the heat sources or sinks
2190	real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
2191	real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
2192	real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
2193	real b_t(kms:kme) ! Explicit component of the heat sources or sinks
2194	real b_ac(kms:kme)
2195	real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
2196	real b_q(kms:kme) ! Explicit component of the humidity sources or sinks
2197	! ----------------------------------------------------------------------
2198	! LOCAL:
2199	! ----------------------------------------------------------------------
2200	real dzz
2201	real fact
2202	integer id,iz,iz_u
2203	real se,sr,st,su,sv,sq
2204	real uet(kms:kme) ! Contribution to TKE due to walls
2205	real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb,vqb,vtb_ac
2206
2207
2208	! ----------------------------------------------------------------------
2209	! END VARIABLES DEFINITIONS
2210	! ----------------------------------------------------------------------
2211
2212	! initialisation
2213
2214	do iz=kts,kte
2215	a_u(iz)=0.
2216	a_v(iz)=0.
2217	a_t(iz)=0.
2218	a_e(iz)=0.
2219	b_u(iz)=0.
2220	b_v(iz)=0.
2221	b_e(iz)=0.
2222	b_t(iz)=0.
2223	b_ac(iz)=0.
2224	uet(iz)=0.
2225	b_q(iz)=0.
2226	end do
2227
2228	! horizontal surfaces
2229	do iz=kts,kte
2230	sf(iz)=0.
2231	vl(iz)=0.
2232	enddo
2233	sf(kte+1)=0.
2234
2235	do id=1,nd
2236	do iz=kts+1,kte+1
2237	sr=0.
2238	do iz_u=2,nz_u
2239	if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
2240	sr=pb(iz_u)
2241	endif
2242	enddo
2243	sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
2244	enddo
2245	enddo
2246
2247	! volume
2248	do iz=kts,kte
2249	do id=1,nd
2250	vtot=0.
2251	do iz_u=1,nz_u
2252	dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
2253	vtot=vtot+pb(iz_u+1)*dzz
2254	enddo
2255	vtot=vtot/(z(iz+1)-z(iz))
2256	vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
2257	enddo
2258	enddo
2259
2260	! horizontal surface impact
2261
2262	do id=1,nd
2263
2264	fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
2265	b_t(kts)=b_t(kts)+thb_u(id,1)*fact
2266	b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
2267	b_v(kts)=b_v(kts)+vhb_u(id,1)*fact
2268	b_e(kts)=b_e(kts)+ehb_u(id,1)fact(z_u(2)-z_u(1))
2269	b_q(kts)=b_q(kts)+qhb_u(id,1)*fact
2270
2271	do iz=kts,kte
2272	st=0.
2273	su=0.
2274	sv=0.
2275	se=0.
2276	sq=0.
2277	do iz_u=2,nz_u
2278	if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
2279	st=st+ss(iz_u)*thb_u(id,iz_u)
2280	su=su+ss(iz_u)*uhb_u(id,iz_u)
2281	sv=sv+ss(iz_u)*vhb_u(id,iz_u)
2282	se=se+ss(iz_u)ehb_u(id,iz_u)(z_u(iz_u+1)-z_u(iz_u))
2283	sq=sq+ss(iz_u)*qhb_u(id,iz_u)
2284	endif
2285	enddo
2286
2287	fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
2288	b_t(iz)=b_t(iz)+st*fact
2289	b_u(iz)=b_u(iz)+su*fact
2290	b_v(iz)=b_v(iz)+sv*fact
2291	b_e(iz)=b_e(iz)+se*fact
2292	b_q(iz)=b_q(iz)+sq*fact
2293	enddo
2294	enddo
2295
2296	! vertical surface impact
2297
2298	do iz=kts,kte
2299	uet(iz)=0.
2300	do id=1,nd
2301	vtb=0.
2302	vtb_ac=0.
2303	vta=0.
2304	vua=0.
2305	vub=0.
2306	vva=0.
2307	vvb=0.
2308	veb=0.
2309	vte=0.
2310	vqb=0.
2311	do iz_u=1,nz_u
2312	dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
2313	fact=dzz/(ws(id)+bs(id))
2314	vtb=vtb+pb(iz_u+1)* &
2315	(tvb_u(2id-1,iz_u)+tvb_u(2id,iz_u))*fact
2316	vtb_ac=vtb_ac+pb(iz_u+1)* &
2317	(tvb_ac(2id-1,iz_u)+tvb_ac(2id,iz_u))*fact
2318	vta=vta+pb(iz_u+1)* &
2319	(tva_u(2id-1,iz_u)+tva_u(2id,iz_u))*fact
2320	vua=vua+pb(iz_u+1)* &
2321	(uva_u(2id-1,iz_u)+uva_u(2id,iz_u))*fact
2322	vva=vva+pb(iz_u+1)* &
2323	(vva_u(2id-1,iz_u)+vva_u(2id,iz_u))*fact
2324	vub=vub+pb(iz_u+1)* &
2325	(uvb_u(2id-1,iz_u)+uvb_u(2id,iz_u))*fact
2326	vvb=vvb+pb(iz_u+1)* &
2327	(vvb_u(2id-1,iz_u)+vvb_u(2id,iz_u))*fact
2328	veb=veb+pb(iz_u+1)* &
2329	(evb_u(2id-1,iz_u)+evb_u(2id,iz_u))*fact
2330	vqb=vqb+pb(iz_u+1)* &
2331	(qvb_u(2id-1,iz_u)+qvb_u(2id,iz_u))*fact
2332	enddo
2333
2334	fact=1./vl(iz)/dz(iz)/nd
2335	b_t(iz)=b_t(iz)+vtb*fact
2336	b_ac(iz)=b_ac(iz)+vtb_ac*fact
2337	a_t(iz)=a_t(iz)+vta*fact
2338	a_u(iz)=a_u(iz)+vua*fact
2339	a_v(iz)=a_v(iz)+vva*fact
2340	b_u(iz)=b_u(iz)+vub*fact
2341	b_v(iz)=b_v(iz)+vvb*fact
2342	b_e(iz)=b_e(iz)+veb*fact
2343	uet(iz)=uet(iz)+vte*fact
2344	b_q(iz)=b_q(iz)+vqb*fact
2345	enddo
2346	enddo
2347
2348
2349
2350	return
2351	end subroutine urban_meso
2352
2353	! ===6=8===============================================================72
2354	! ===6=8===============================================================72
2355
2356	subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
2357	dlg,dl_u)
2358
2359	! ----------------------------------------------------------------------
2360	! Calculation of the length scales
2361	! See p272-274 formula (22) and (24) of the BLM paper
2362	! ----------------------------------------------------------------------
2363
2364	implicit none
2365
2366
2367	! ----------------------------------------------------------------------
2368	! INPUT:
2369	! ----------------------------------------------------------------------
2370	integer kms,kme,kts,kte
2371	real z(kms:kme) ! Altitude above the ground of the cell interface
2372	integer nd ! Number of street direction for the current urban class
2373	integer nz_u ! Number of levels in the "urban grid"
2374	real z_u(nz_um) ! Height of the urban grid levels
2375	real bs(ndm) ! Building widths of the current urban class
2376	real ss(nz_um) ! Probability to have a building with height h
2377	real ws(ndm) ! Street widths of the current urban class
2378
2379
2380	! ----------------------------------------------------------------------
2381	! OUTPUT:
2382	! ----------------------------------------------------------------------
2383	real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
2384	real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
2385
2386	! ----------------------------------------------------------------------
2387	! LOCAL:
2388	! ----------------------------------------------------------------------
2389	real dlgtmp
2390	integer id,iz,iz_u
2391	real sftot
2392	real ulu,ssl
2393
2394	! ----------------------------------------------------------------------
2395	! END VARIABLES DEFINITIONS
2396	! ----------------------------------------------------------------------
2397
2398	do iz=kts,kte
2399	ulu=0.
2400	ssl=0.
2401	do id=1,nd
2402	do iz_u=2,nz_u
2403	if(z_u(iz_u).gt.z(iz))then
2404	ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
2405	ssl=ssl+ss(iz_u)/nd
2406	endif
2407	enddo
2408	enddo
2409
2410	if(ulu.ne.0)then
2411	dl_u(iz)=ssl/ulu
2412	else
2413	dl_u(iz)=0.
2414	endif
2415	enddo
2416
2417
2418	do iz=kts,kte
2419	dlg(iz)=0.
2420	do id=1,nd
2421	sftot=ws(id)
2422	dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
2423	do iz_u=1,nz_u
2424	if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
2425	dlgtmp=dlgtmp+ss(iz_u)*bs(id)/ &
2426	((z(iz)+z(iz+1))/2.-z_u(iz_u))
2427	sftot=sftot+ss(iz_u)*bs(id)
2428	endif
2429	enddo
2430	dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
2431	enddo
2432	dlg(iz)=1./dlg(iz)
2433	enddo
2434
2435	return
2436	end subroutine interp_length
2437
2438	! ===6=8===============================================================72
2439	! ===6=8===============================================================72
2440
2441	subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
2442	rs,rsw,rsg)
2443
2444	! ----------------------------------------------------------------------
2445	! Modification of short wave radiation to take into account
2446	! the shadow produced by the buildings
2447	! ----------------------------------------------------------------------
2448
2449	implicit none
2450
2451	! ----------------------------------------------------------------------
2452	! INPUT:
2453	! ----------------------------------------------------------------------
2454	integer nd ! Number of street direction for the current urban class
2455	integer nz_u ! number of vertical layers defined in the urban grid
2456	real ah ! Hour angle (it should come from the radiation routine)
2457	real deltar ! Declination of the sun
2458	real drst(ndm) ! street directions for the current urban class
2459	real rs ! solar radiation
2460	real ss(nz_um) ! probability to have a building with height h
2461	real pb(nz_um) ! Probability that a building has an height greater or equal to h
2462	real ws(ndm) ! Street width of the current urban class
2463	real z(nz_um) ! Height of the urban grid levels
2464	real zr ! zenith angle
2465
2466	! ----------------------------------------------------------------------
2467	! OUTPUT:
2468	! ----------------------------------------------------------------------
2469	real rsg(ndm) ! Short wave radiation at the ground
2470	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
2471
2472	! ----------------------------------------------------------------------
2473	! LOCAL:
2474	! ----------------------------------------------------------------------
2475	integer id,iz,jz
2476	real aae,aaw,bbb,phix,rd,rtot,wsd
2477
2478	! ----------------------------------------------------------------------
2479	! END VARIABLES DEFINITIONS
2480	! ----------------------------------------------------------------------
2481
2482	if(rs.eq.0.or.sin(zr).eq.1)then
2483	do id=1,nd
2484	rsg(id)=0.
2485	do iz=1,nz_u
2486	rsw(2*id-1,iz)=0.
2487	rsw(2*id,iz)=0.
2488	enddo
2489	enddo
2490	else
2491	!test
2492	if(abs(sin(zr)).gt.1.e-10)then
2493	if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
2494	bbb=pi/2.
2495	elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
2496	bbb=-pi/2.
2497	else
2498	bbb=asin(cos(deltar)*sin(ah)/sin(zr))
2499	endif
2500	else
2501	if(cos(deltar)*sin(ah).ge.0)then
2502	bbb=pi/2.
2503	elseif(cos(deltar)*sin(ah).lt.0)then
2504	bbb=-pi/2.
2505	endif
2506	endif
2507
2508	phix=zr
2509
2510	do id=1,nd
2511
2512	rsg(id)=0.
2513
2514	aae=bbb-drst(id)
2515	aaw=bbb-drst(id)+pi
2516
2517	do iz=1,nz_u
2518	rsw(2*id-1,iz)=0.
2519	rsw(2*id,iz)=0.
2520	if (pb(iz+1).gt.0.) then
2521	do jz=1,nz_u
2522	if(abs(sin(aae)).gt.1.e-10)then
2523	call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae, &
2524	ws(id),rd)
2525	rsw(2id-1,iz)=rsw(2id-1,iz)+rsrdss(jz+1)/pb(iz+1)
2526	endif
2527
2528	if(abs(sin(aaw)).gt.1.e-10)then
2529	call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw, &
2530	ws(id),rd)
2531	rsw(2id,iz)=rsw(2id,iz)+rsrdss(jz+1)/pb(iz+1)
2532	endif
2533	enddo
2534	endif
2535	enddo
2536	if(abs(sin(aae)).gt.1.e-10)then
2537	wsd=abs(ws(id)/sin(aae))
2538
2539	do jz=1,nz_u
2540	rd=max(0.,wsd-z(jz+1)*tan(phix))
2541	rsg(id)=rsg(id)+rsrdss(jz+1)/wsd
2542	enddo
2543	rtot=0.
2544
2545	do iz=1,nz_u
2546	rtot=rtot+(rsw(2id,iz)+rsw(2id-1,iz))* &
2547	(z(iz+1)-z(iz))
2548	enddo
2549	rtot=rtot+rsg(id)*ws(id)
2550	else
2551	rsg(id)=rs
2552	endif
2553
2554	enddo
2555	endif
2556
2557	return
2558	end subroutine shadow_mas
2559
2560	! ===6=8===============================================================72
2561	! ===6=8===============================================================72
2562
2563	subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
2564
2565	! ----------------------------------------------------------------------
2566	! This routine computes the effects of a shadow induced by a building of
2567	! height hu, on a portion of wall between z1 and z2. See equation A10,
2568	! and correction described below formula A11, and figure A1. Basically rd
2569	! is the ratio between the horizontal surface illuminated and the portion
2570	! of wall. Referring to figure A1, multiplying radiation flux density on
2571	! a horizontal surface (rs) by x1-x2 we have the radiation energy per
2572	! unit time. Dividing this by z2-z1, we obtain the radiation flux
2573	! density reaching the portion of the wall between z2 and z1
2574	! (everything is assumed in 2D)
2575	! ----------------------------------------------------------------------
2576
2577	implicit none
2578
2579	! ----------------------------------------------------------------------
2580	! INPUT:
2581	! ----------------------------------------------------------------------
2582	real aa ! Angle between the sun direction and the face of the wall (A12)
2583	real hu ! Height of the building that generates the shadow
2584	real phix ! Solar zenith angle
2585	real ws ! Width of the street
2586	real z1 ! Height of the level z(iz)
2587	real z2 ! Height of the level z(iz+1)
2588
2589	! ----------------------------------------------------------------------
2590	! OUTPUT:
2591	! ----------------------------------------------------------------------
2592	real rd ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A.
2593	! Multiplying rd by rs (radiation flux
2594	! density on a horizontal surface) gives
2595	! the radiation flux density on the
2596	! portion of wall between z1 and z2.
2597	! ----------------------------------------------------------------------
2598	! LOCAL:
2599	! ----------------------------------------------------------------------
2600	real x1,x2 ! x1,x2 see Fig. A1.
2601
2602	! ----------------------------------------------------------------------
2603	! END VARIABLES DEFINITIONS
2604	! ----------------------------------------------------------------------
2605
2606	x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
2607
2608	x2=max((hu-z2)*tan(phix),0.)
2609
2610	rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
2611
2612	return
2613	end subroutine shade_wall
2614
2615	! ===6=8===============================================================72
2616	! ===6=8===============================================================72
2617
2618	subroutine long_rad(iurb,nz_u,id,emw,emg,emwin,pwin,twlev,&
2619	fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
2620
2621	! ----------------------------------------------------------------------
2622	! This routine computes the effects of the reflections of long-wave
2623	! radiation in the street canyon by solving the system
2624	! of 2nz_u+1 eqn. in 2nz_u+1
2625	! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
2626	! The system is solved by solving A X= B,
2627	! with A matrix B vector, and X solution.
2628	! ----------------------------------------------------------------------
2629
2630	implicit none
2631
2632
2633
2634	! ----------------------------------------------------------------------
2635	! INPUT:
2636	! ----------------------------------------------------------------------
2637	real emg ! Emissivity of ground for the current urban class
2638	real emw ! Emissivity of wall for the current urban class
2639	real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
2640	real fsg(ndm,nurbm) ! View factors from sky to ground
2641	real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
2642	real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
2643	real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
2644	integer id ! Current street direction
2645	integer iurb ! Current urban class
2646	integer nz_u ! Number of layer in the urban grid
2647	real pb(nz_um) ! Probability to have a building with an height equal
2648	real rl ! Downward flux of the longwave radiation
2649	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
2650	real tw(2*ndm,nz_um) ! Temperature in each layer of the wall [K]
2651	!
2652	!New Variables for BEM
2653	!
2654	real twlev(2*ndm,nz_um) ! Window temperature in BEM [K]
2655	real emwin ! Emissivity of windows
2656	real pwin ! Coverage area fraction of windows in the walls of the buildings (BEM)
2657
2658
2659	! ----------------------------------------------------------------------
2660	! OUTPUT:
2661	! ----------------------------------------------------------------------
2662	real rlg(ndm) ! Long wave radiation at the ground
2663	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
2664
2665	! ----------------------------------------------------------------------
2666	! LOCAL:
2667	! ----------------------------------------------------------------------
2668	integer i,j
2669	real aaa(2nz_um+1,2nz_um+1) ! terms of the matrix
2670	real bbb(2*nz_um+1) ! terms of the vector
2671
2672	! ----------------------------------------------------------------------
2673	! END VARIABLES DEFINITIONS
2674	! ----------------------------------------------------------------------
2675
2676
2677	! west wall
2678
2679	do i=1,nz_u
2680
2681	do j=1,nz_u
2682	aaa(i,j)=0.
2683	enddo
2684
2685	aaa(i,i)=1.
2686
2687	do j=nz_u+1,2*nz_u
2688	aaa(i,j)=-(1.-emw(1.-pwin)-emwinpwin)* &
2689	fww(j-nz_u,i,id,iurb)*pb(j-nz_u+1)
2690	enddo
2691
2692	aaa(i,2nz_u+1)=-(1.-emg)fgw(i,id,iurb)
2693
2694	bbb(i)=fsw(i,id,iurb)rl+emgfgw(i,id,iurb)sigmatg(id,ng_u)**4
2695	do j=1,nz_u
2696	bbb(i)=bbb(i)+pb(j+1)sigmafww(j,i,id,iurb)* &
2697	(emw(1.-pwin)tw(2id,j)4+emwinpwintwlev(2id,j)**4)+ &
2698	fww(j,i,id,iurb)rl(1.-pb(j+1))
2699	enddo
2700
2701	enddo
2702
2703	! east wall
2704
2705	do i=1+nz_u,2*nz_u
2706
2707	do j=1,nz_u
2708	aaa(i,j)=-(1.-emw(1.-pwin)-emwinpwin)fww(j,i-nz_u,id,iurb)pb(j+1)
2709	enddo
2710
2711	do j=1+nz_u,2*nz_u
2712	aaa(i,j)=0.
2713	enddo
2714
2715	aaa(i,i)=1.
2716
2717	aaa(i,2nz_u+1)=-(1.-emg)fgw(i-nz_u,id,iurb)
2718
2719	bbb(i)=fsw(i-nz_u,id,iurb)*rl+ &
2720	emgfgw(i-nz_u,id,iurb)sigmatg(id,ng_u)*4
2721
2722	do j=1,nz_u
2723	bbb(i)=bbb(i)+pb(j+1)sigmafww(j,i-nz_u,id,iurb)* &
2724	(emw(1.-pwin)tw(2id-1,j)4+emwinpwintwlev(2id-1,j)**4)+&
2725	fww(j,i-nz_u,id,iurb)rl(1.-pb(j+1))
2726	enddo
2727
2728	enddo
2729
2730	! ground
2731	do j=1,nz_u
2732	aaa(2nz_u+1,j)=-(1.-emw(1.-pwin)-emwinpwin) &
2733	fwg(j,id,iurb)*pb(j+1)
2734	enddo
2735
2736	do j=nz_u+1,2*nz_u
2737	aaa(2nz_u+1,j)=-(1.-emw(1.-pwin)-emwinpwin) &
2738	fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
2739	enddo
2740
2741	aaa(2nz_u+1,2nz_u+1)=1.
2742
2743	bbb(2nz_u+1)=fsg(id,iurb)rl
2744
2745	do i=1,nz_u
2746	bbb(2nz_u+1)=bbb(2nz_u+1)+sigmafwg(i,id,iurb)pb(i+1)* &
2747	(emw(1.-pwin)(tw(2id-1,i)4+tw(2id,i)**4)+ &
2748	emwinpwin(twlev(2id-1,i)4+twlev(2id,i)**4))+ &
2749	2.fwg(i,id,iurb)(1.-pb(i+1))*rl
2750	enddo
2751
2752
2753
2754	call gaussj(aaa,2nz_u+1,bbb,2nz_um+1)
2755
2756	do i=1,nz_u
2757	rlw(2*id-1,i)=bbb(i)
2758	enddo
2759
2760	do i=nz_u+1,2*nz_u
2761	rlw(2*id,i-nz_u)=bbb(i)
2762	enddo
2763
2764	rlg(id)=bbb(2*nz_u+1)
2765
2766	return
2767	end subroutine long_rad
2768
2769	! ===6=8===============================================================72
2770	! ===6=8===============================================================72
2771
2772	subroutine short_rad(iurb,nz_u,id,albw, &
2773	albg,fwg,fww,fgw,rsg,rsw,pb)
2774
2775	! ----------------------------------------------------------------------
2776	! This routine computes the effects of the reflections of short-wave
2777	! (solar) radiation in the street canyon by solving the system
2778	! of 2nz_u+1 eqn. in 2nz_u+1
2779	! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
2780	! The system is solved by solving A X= B,
2781	! with A matrix B vector, and X solution.
2782	! ----------------------------------------------------------------------
2783
2784	implicit none
2785
2786
2787
2788	! ----------------------------------------------------------------------
2789	! INPUT:
2790	! ----------------------------------------------------------------------
2791	real albg ! Albedo of the ground for the current urban class
2792	real albw ! Albedo of the wall for the current urban class
2793	real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
2794	real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
2795	real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
2796	integer id ! current street direction
2797	integer iurb ! current urban class
2798	integer nz_u ! Number of layer in the urban grid
2799	real pb(nz_um) ! probability to have a building with an height equal
2800
2801	! ----------------------------------------------------------------------
2802	! OUTPUT:
2803	! ----------------------------------------------------------------------
2804	real rsg(ndm) ! Short wave radiation at the ground
2805	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
2806
2807	! ----------------------------------------------------------------------
2808	! LOCAL:
2809	! ----------------------------------------------------------------------
2810	integer i,j
2811	real aaa(2nz_um+1,2nz_um+1) ! terms of the matrix
2812	real bbb(2*nz_um+1) ! terms of the vector
2813
2814	! ----------------------------------------------------------------------
2815	! END VARIABLES DEFINITIONS
2816	! ----------------------------------------------------------------------
2817
2818
2819	! west wall
2820
2821	do i=1,nz_u
2822	do j=1,nz_u
2823	aaa(i,j)=0.
2824	enddo
2825
2826	aaa(i,i)=1.
2827
2828	do j=nz_u+1,2*nz_u
2829	aaa(i,j)=-albwfww(j-nz_u,i,id,iurb)pb(j-nz_u+1)
2830	enddo
2831
2832	aaa(i,2nz_u+1)=-albgfgw(i,id,iurb)
2833	bbb(i)=rsw(2*id-1,i)
2834
2835	enddo
2836
2837	! east wall
2838
2839	do i=1+nz_u,2*nz_u
2840	do j=1,nz_u
2841	aaa(i,j)=-albwfww(j,i-nz_u,id,iurb)pb(j+1)
2842	enddo
2843
2844	do j=1+nz_u,2*nz_u
2845	aaa(i,j)=0.
2846	enddo
2847
2848	aaa(i,i)=1.
2849	aaa(i,2nz_u+1)=-albgfgw(i-nz_u,id,iurb)
2850	bbb(i)=rsw(2*id,i-nz_u)
2851
2852	enddo
2853
2854	! ground
2855
2856	do j=1,nz_u
2857	aaa(2nz_u+1,j)=-albwfwg(j,id,iurb)*pb(j+1)
2858	enddo
2859
2860	do j=nz_u+1,2*nz_u
2861	aaa(2nz_u+1,j)=-albwfwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
2862	enddo
2863
2864	aaa(2nz_u+1,2nz_u+1)=1.
2865	bbb(2*nz_u+1)=rsg(id)
2866
2867	call gaussj(aaa,2nz_u+1,bbb,2nz_um+1)
2868
2869	do i=1,nz_u
2870	rsw(2*id-1,i)=bbb(i)
2871	enddo
2872
2873	do i=nz_u+1,2*nz_u
2874	rsw(2*id,i-nz_u)=bbb(i)
2875	enddo
2876
2877	rsg(id)=bbb(2*nz_u+1)
2878
2879	return
2880	end subroutine short_rad
2881
2882
2883	! ===6=8===============================================================72
2884	! ===6=8===============================================================72
2885
2886	subroutine gaussj(a,n,b,np)
2887
2888	! ----------------------------------------------------------------------
2889	! This routine solve a linear system of n equations of the form
2890	! A X = B
2891	! where A is a matrix a(i,j)
2892	! B a vector and X the solution
2893	! In output b is replaced by the solution
2894	! ----------------------------------------------------------------------
2895
2896	implicit none
2897
2898	! ----------------------------------------------------------------------
2899	! INPUT:
2900	! ----------------------------------------------------------------------
2901	integer np
2902	real a(np,np)
2903
2904	! ----------------------------------------------------------------------
2905	! OUTPUT:
2906	! ----------------------------------------------------------------------
2907	real b(np)
2908
2909	! ----------------------------------------------------------------------
2910	! LOCAL:
2911	! ----------------------------------------------------------------------
2912	integer nmax
2913	parameter (nmax=150)
2914
2915	real big,dum
2916	integer i,icol,irow
2917	integer j,k,l,ll,n
2918	integer ipiv(nmax)
2919	real pivinv
2920
2921	! ----------------------------------------------------------------------
2922	! END VARIABLES DEFINITIONS
2923	! ----------------------------------------------------------------------
2924
2925	do j=1,n
2926	ipiv(j)=0.
2927	enddo
2928
2929	do i=1,n
2930	big=0.
2931	do j=1,n
2932	if(ipiv(j).ne.1)then
2933	do k=1,n
2934	if(ipiv(k).eq.0)then
2935	if(abs(a(j,k)).ge.big)then
2936	big=abs(a(j,k))
2937	irow=j
2938	icol=k
2939	endif
2940	elseif(ipiv(k).gt.1)then
2941	pause 'singular matrix in gaussj'
2942	endif
2943	enddo
2944	endif
2945	enddo
2946
2947	ipiv(icol)=ipiv(icol)+1
2948
2949	if(irow.ne.icol)then
2950	do l=1,n
2951	dum=a(irow,l)
2952	a(irow,l)=a(icol,l)
2953	a(icol,l)=dum
2954	enddo
2955
2956	dum=b(irow)
2957	b(irow)=b(icol)
2958	b(icol)=dum
2959
2960	endif
2961
2962	if(a(icol,icol).eq.0)pause 'singular matrix in gaussj'
2963
2964	pivinv=1./a(icol,icol)
2965	a(icol,icol)=1
2966
2967	do l=1,n
2968	a(icol,l)=a(icol,l)*pivinv
2969	enddo
2970
2971	b(icol)=b(icol)*pivinv
2972
2973	do ll=1,n
2974	if(ll.ne.icol)then
2975	dum=a(ll,icol)
2976	a(ll,icol)=0.
2977	do l=1,n
2978	a(ll,l)=a(ll,l)-a(icol,l)*dum
2979	enddo
2980
2981	b(ll)=b(ll)-b(icol)*dum
2982
2983	endif
2984	enddo
2985	enddo
2986
2987	return
2988	end subroutine gaussj
2989
2990	! ===6=8===============================================================72
2991	! ===6=8===============================================================72
2992
2993	subroutine soil_temp(nz,dz,temp,pt,ala,cs, &
2994	rs,rl,press,dt,em,alb,rt,sf,gf)
2995
2996	! ----------------------------------------------------------------------
2997	! This routine solves the Fourier diffusion equation for heat in
2998	! the material (wall, roof, or ground). Resolution is done implicitely.
2999	! Boundary conditions are:
3000	! - fixed temperature at the interior
3001	! - energy budget at the surface
3002	! ----------------------------------------------------------------------
3003
3004	implicit none
3005
3006
3007
3008	! ----------------------------------------------------------------------
3009	! INPUT:
3010	! ----------------------------------------------------------------------
3011	integer nz ! Number of layers
3012	real ala(nz) ! Thermal diffusivity in each layers [m^2 s^-1]
3013	real alb ! Albedo of the surface
3014	real cs(nz) ! Specific heat of the material [J m^3 K^-1]
3015	real dt ! Time step
3016	real em ! Emissivity of the surface
3017	real press ! Pressure at ground level
3018	real rl ! Downward flux of the longwave radiation
3019	real rs ! Solar radiation
3020	real sf ! Sensible heat flux at the surface
3021	real temp(nz) ! Temperature in each layer [K]
3022	real dz(nz) ! Layer sizes [m]
3023
3024
3025	! ----------------------------------------------------------------------
3026	! OUTPUT:
3027	! ----------------------------------------------------------------------
3028	real gf ! Heat flux transferred from the surface toward the interior
3029	real pt ! Potential temperature at the surface
3030	real rt ! Total radiation at the surface (solar+incoming long+outgoing long)
3031
3032	! ----------------------------------------------------------------------
3033	! LOCAL:
3034	! ----------------------------------------------------------------------
3035	integer iz
3036	real a(nz,3)
3037	real alpha
3038	real c(nz)
3039	real cddz(nz+2)
3040	real tsig
3041
3042	! ----------------------------------------------------------------------
3043	! END VARIABLES DEFINITIONS
3044	! ----------------------------------------------------------------------
3045
3046	tsig=temp(nz)
3047	alpha=(1.-alb)rs+emrl-emsigma(tsig**4)+sf
3048	! Compute cddz=2*cd/dz
3049
3050	cddz(1)=ala(1)/dz(1)
3051	do iz=2,nz
3052	cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
3053	enddo
3054
3055	a(1,1)=0.
3056	a(1,2)=1.
3057	a(1,3)=0.
3058	c(1)=temp(1)
3059
3060	do iz=2,nz-1
3061	a(iz,1)=-cddz(iz)*dt/dz(iz)
3062	a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
3063	a(iz,3)=-cddz(iz+1)*dt/dz(iz)
3064	c(iz)=temp(iz)
3065	enddo
3066
3067	a(nz,1)=-dt*cddz(nz)/dz(nz)
3068	a(nz,2)=1.+dt*cddz(nz)/dz(nz)
3069	a(nz,3)=0.
3070	c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
3071
3072
3073	call invert(nz,a,c,temp)
3074
3075
3076	pt=temp(nz)(press/1.e+5)*(-rcp_u)
3077
3078	rt=(1.-alb)rs+emrl-emsigma(tsig**4)
3079
3080	gf=(1.-alb)rs+emrl-emsigma(tsig**4)+sf
3081	return
3082	end subroutine soil_temp
3083
3084	! ===6=8===============================================================72
3085	! ===6=8===============================================================72
3086
3087	subroutine invert(n,a,c,x)
3088
3089	! ----------------------------------------------------------------------
3090	! Inversion and resolution of a tridiagonal matrix
3091	! A X = C
3092	! ----------------------------------------------------------------------
3093
3094	implicit none
3095
3096	! ----------------------------------------------------------------------
3097	! INPUT:
3098	! ----------------------------------------------------------------------
3099	integer n
3100	real a(n,3) ! a(*,1) lower diagonal (Ai,i-1)
3101	! a(*,2) principal diagonal (Ai,i)
3102	! a(*,3) upper diagonal (Ai,i+1)
3103	real c(n)
3104
3105	! ----------------------------------------------------------------------
3106	! OUTPUT:
3107	! ----------------------------------------------------------------------
3108	real x(n)
3109
3110	! ----------------------------------------------------------------------
3111	! LOCAL:
3112	! ----------------------------------------------------------------------
3113	integer i
3114
3115	! ----------------------------------------------------------------------
3116	! END VARIABLES DEFINITIONS
3117	! ----------------------------------------------------------------------
3118
3119	do i=n-1,1,-1
3120	c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
3121	a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
3122	enddo
3123
3124	do i=2,n
3125	c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
3126	enddo
3127
3128	do i=1,n
3129	x(i)=c(i)/a(i,2)
3130	enddo
3131
3132	return
3133	end subroutine invert
3134
3135
3136	! ===6=8===============================================================72
3137	! ===6=8===============================================================72
3138
3139	subroutine flux_wall(ua,va,pt,da,ptw,ptwin,uva,vva,uvb,vvb, &
3140	sfw,sfwin,evb,drst,dt)
3141
3142	! ----------------------------------------------------------------------
3143	! This routine computes the surface sources or sinks of momentum, tke,
3144	! and heat from vertical surfaces (walls).
3145	! ----------------------------------------------------------------------
3146	implicit none
3147
3148	! INPUT:
3149	! -----
3150	real drst ! street directions for the current urban class
3151	real da ! air density
3152	real pt ! potential temperature
3153	real ptw ! Walls potential temperatures
3154	real ptwin ! Windows potential temperatures
3155	real ua ! wind speed
3156	real va ! wind speed
3157	real dt !time step
3158
3159	! OUTPUT:
3160	! ------
3161	! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
3162	! vertical surfaces (walls).
3163	! The fluxes can be computed as follow: Fluxes of X = A*X + B
3164	! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
3165
3166	real uva ! U (wind component) Vertical surfaces, A (implicit) term
3167	real uvb ! U (wind component) Vertical surfaces, B (explicit) term
3168	real vva ! V (wind component) Vertical surfaces, A (implicit) term
3169	real vvb ! V (wind component) Vertical surfaces, B (explicit) term
3170	! real tva ! Temperature Vertical surfaces, A (implicit) term
3171	! real tvb ! Temperature Vertical surfaces, B (explicit) term
3172	real evb ! Energy (TKE) Vertical surfaces, B (explicit) term
3173	real sfw ! Surfaces fluxes from the walls
3174	real sfwin ! Surfaces fluxes from the windows
3175
3176	! LOCAL:
3177	! -----
3178	real hc
3179	real hcwin
3180	real u_ort
3181	real vett
3182
3183
3184	! -------------------------
3185	! END VARIABLES DEFINITIONS
3186	! -------------------------
3187
3188	vett=(ua2+va2)**.5
3189
3190	u_ort=abs((cos(drst)ua-sin(drst)va))
3191
3192	uva=-cdragu_ort/2.cos(drst)*cos(drst)
3193	vva=-cdragu_ort/2.sin(drst)*sin(drst)
3194
3195	uvb=cdragu_ort/2.sin(drst)cos(drst)va
3196	vvb=cdragu_ort/2.sin(drst)cos(drst)ua
3197
3198	if (vett.lt.4.88) then
3199	hc=5.678(1.09+0.23(vett/0.3048))
3200	else
3201	hc=5.6780.53((vett/0.3048)**0.78)
3202	endif
3203
3204	if (hc.gt.da*cp_u/dt)then
3205	hc=da*cp_u/dt
3206	endif
3207
3208	if (vett.lt.4.88) then
3209	hcwin=5.678(0.99+0.21(vett/0.3048))
3210	else
3211	hcwin=5.6780.50((vett/0.3048)**0.78)
3212	endif
3213
3214	if (hcwin.gt.da*cp_u/dt) then
3215	hcwin=da*cp_u/dt
3216	endif
3217
3218	! tvb=hc*ptw/da/cp_u
3219	! tva=-hc/da/cp_u
3220	!!!!!!!!!!!!!!!!!!!!
3221	! explicit
3222
3223	sfw=hc*(pt-ptw)
3224	sfwin=hcwin*(pt-ptwin)
3225
3226
3227	evb=cdrag(abs(u_ort)*3.)/2.
3228
3229	return
3230	end subroutine flux_wall
3231
3232	! ===6=8===============================================================72
3233	! ===6=8===============================================================72
3234
3235	subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg, &
3236	uhb,vhb,sf,ehb,da)
3237
3238	! ----------------------------------------------------------------------
3239	! Calculation of the flux at the ground
3240	! Formulation of Louis (Louis, 1979)
3241	! ----------------------------------------------------------------------
3242
3243	implicit none
3244
3245	real dz ! first vertical level
3246	real pt ! potential temperature
3247	real pt0 ! reference potential temperature
3248	real ptg ! ground potential temperature
3249	real ua ! wind speed
3250	real va ! wind speed
3251	real z0 ! Roughness length
3252	real da ! air density
3253
3254
3255
3256	! ----------------------------------------------------------------------
3257	! OUTPUT:
3258	! ----------------------------------------------------------------------
3259	! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
3260	! surfaces (roofs and street)
3261	! The fluxes can be computed as follow: Fluxes of X = B
3262	! Example: Momentum fluxes on horizontal surfaces = uhb_u
3263	real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
3264	real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
3265	! real thb ! Temperature Horizontal surfaces, B (explicit) term
3266	! real tva ! Temperature Vertical surfaces, A (implicit) term
3267	! real tvb ! Temperature Vertical surfaces, B (explicit) term
3268	real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
3269	real sf
3270
3271
3272	! ----------------------------------------------------------------------
3273	! LOCAL:
3274	! ----------------------------------------------------------------------
3275	real aa
3276	real al
3277	real buu
3278	real c
3279	real fbuw
3280	real fbpt
3281	real fh
3282	real fm
3283	real ric
3284	real tstar
3285	real ustar
3286	real utot
3287	real wstar
3288	real zz
3289
3290	real b,cm,ch,rr,tol
3291	parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
3292
3293	! ----------------------------------------------------------------------
3294	! END VARIABLES DEFINITIONS
3295	! ----------------------------------------------------------------------
3296
3297
3298	! computation of the ground temperature
3299
3300	utot=(ua2+va2)**.5
3301
3302
3303	!!!! Louis formulation
3304	!
3305	! compute the bulk Richardson Number
3306
3307	zz=dz/2.
3308
3309
3310	utot=max(utot,0.01)
3311
3312	ric=2.g_uzz(pt-ptg)/((pt+ptg)(utot**2))
3313
3314	aa=vk/log(zz/z0)
3315
3316	! determine the parameters fm and fh for stable, neutral and unstable conditions
3317
3318	if(ric.gt.0)then
3319	fm=1/(1+0.5bric)**2
3320	fh=fm
3321	else
3322	c=bcmaaaa(zz/z0)**.5
3323	fm=1-bric/(1+c(-ric)**.5)
3324	c=c*ch/cm
3325	fh=1-bric/(1+c(-ric)**.5)
3326	endif
3327
3328	fbuw=-aaaautotutotfm
3329	fbpt=-aaaautot(pt-ptg)fh/rr
3330
3331	ustar=(-fbuw)**.5
3332	tstar=-fbpt/ustar
3333
3334	al=(vkg_utstar)/(ptustarustar)
3335
3336	buu=-g_u/pt0ustartstar
3337
3338	uhb=-ustarustarua/utot
3339	vhb=-ustarustarva/utot
3340	sf= ustartstarda*cp_u
3341
3342	! thb= 0.
3343	ehb=buu
3344	!!!!!!!!!!!!!!!
3345
3346	return
3347	end subroutine flux_flat
3348
3349	! ===6=8===============================================================72
3350	! ===6=8===============================================================72
3351
3352	subroutine icBEP (fww,fwg,fgw,fsw,fws,fsg, &
3353	z0g_u,z0r_u, &
3354	nd_u,strd_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
3355	nz_u,z_u)
3356
3357
3358	implicit none
3359
3360
3361	! Building parameters
3362
3363	! Radiation parameters
3364
3365	! Roughness parameters
3366	real z0g_u(nurbm) ! The ground's roughness length
3367	real z0r_u(nurbm) ! The roof's roughness length
3368
3369	! Street parameters
3370	integer nd_u(nurbm) ! Number of street direction for each urban class
3371
3372	real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
3373	real ws_u(ndm,nurbm) ! Street width [m]
3374	real bs_u(ndm,nurbm) ! Building width [m]
3375	real h_b(nz_um,nurbm) ! Bulding's heights [m]
3376	real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
3377	! -----------------------------------------------------------------------
3378	! Output
3379	!------------------------------------------------------------------------
3380
3381
3382
3383	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
3384	! and the short wave radation. They are the part of radiation from a surface
3385	! or from the sky to another surface.
3386	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
3387	real fwg(nz_um,ndm,nurbm) ! from wall to ground
3388	real fgw(nz_um,ndm,nurbm) ! from ground to wall
3389	real fsw(nz_um,ndm,nurbm) ! from sky to wall
3390	real fws(nz_um,ndm,nurbm) ! from wall to sky
3391	real fsg(ndm,nurbm) ! from sky to ground
3392
3393	real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to z
3394	real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to z
3395
3396	! Grid parameters
3397	integer nz_u(nurbm) ! Number of layer in the urban grid
3398	real z_u(nz_um) ! Height of the urban grid levels
3399
3400
3401	! -----------------------------------------------------------------------
3402	! Local
3403	!------------------------------------------------------------------------
3404
3405	integer iz_u,id,ilu,iurb
3406
3407	real dtot
3408	real hbmax
3409
3410	!------------------------------------------------------------------------
3411
3412
3413	! -----------------------------------------------------------------------
3414	! This routine initialise the urban paramters for the BEP module
3415	!------------------------------------------------------------------------
3416	!
3417	!Initialize some variables
3418	!
3419	nz_u=0
3420	z_u=0.
3421	ss_u=0.
3422	pb_u=0.
3423	fww=0.
3424	fwg=0.
3425	fgw=0.
3426	fsw=0.
3427	fws=0.
3428	fsg=0.
3429
3430	! Computation of the urban levels height
3431
3432	z_u(1)=0.
3433
3434	do iz_u=1,nz_um-1
3435	z_u(iz_u+1)=z_u(iz_u)+dz_u
3436	enddo
3437
3438	! Normalisation of the building density
3439
3440	do iurb=1,nurbm
3441	dtot=0.
3442	do ilu=1,nz_um
3443	dtot=dtot+d_b(ilu,iurb)
3444	enddo
3445	do ilu=1,nz_um
3446	d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
3447	enddo
3448	enddo
3449
3450	! Compute the view factors, pb and ss
3451
3452	do iurb=1,nurbm
3453	hbmax=0.
3454	nz_u(iurb)=0
3455	do ilu=1,nz_um
3456	if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
3457	enddo
3458
3459	do iz_u=1,nz_um-1
3460	if(z_u(iz_u+1).gt.hbmax)go to 10
3461	enddo
3462
3463	10 continue
3464	nz_u(iurb)=iz_u+1
3465
3466	do id=1,nd_u(iurb)
3467
3468	call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb), &
3469	z_u,ws_u(id,iurb), &
3470	fww,fwg,fgw,fsg,fsw,fws)
3471
3472	do iz_u=1,nz_u(iurb)
3473	ss_u(iz_u,iurb)=0.
3474	do ilu=1,nz_um
3475	if(z_u(iz_u).le.h_b(ilu,iurb) &
3476	.and.z_u(iz_u+1).gt.h_b(ilu,iurb))then
3477	ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
3478	endif
3479	enddo
3480	enddo
3481
3482	pb_u(1,iurb)=1.
3483	do iz_u=1,nz_u(iurb)
3484	pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
3485	enddo
3486
3487	enddo
3488	end do
3489
3490
3491	return
3492	end subroutine icBEP
3493
3494	! ===6=8===============================================================72
3495	! ===6=8===============================================================72
3496
3497	subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws)
3498
3499	implicit none
3500
3501
3502
3503	! -----------------------------------------------------------------------
3504	! Input
3505	!------------------------------------------------------------------------
3506
3507	integer iurb ! Number of the urban class
3508	integer nz_u ! Number of levels in the urban grid
3509	integer id ! Street direction number
3510	real ws ! Street width
3511	real z(nz_um) ! Height of the urban grid levels
3512	real dxy ! Street lenght
3513
3514
3515	! -----------------------------------------------------------------------
3516	! Output
3517	!------------------------------------------------------------------------
3518
3519	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
3520	! and the short wave radation. They are the part of radiation from a surface
3521	! or from the sky to another surface.
3522
3523	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
3524	real fwg(nz_um,ndm,nurbm) ! from wall to ground
3525	real fgw(nz_um,ndm,nurbm) ! from ground to wall
3526	real fsw(nz_um,ndm,nurbm) ! from sky to wall
3527	real fws(nz_um,ndm,nurbm) ! from wall to sky
3528	real fsg(ndm,nurbm) ! from sky to ground
3529
3530
3531	! -----------------------------------------------------------------------
3532	! Local
3533	!------------------------------------------------------------------------
3534
3535	integer jz,iz
3536
3537	real hut
3538	real f1,f2,f12,f23,f123,ftot
3539	real fprl,fnrm
3540	real a1,a2,a3,a4,a12,a23,a123
3541
3542	! -----------------------------------------------------------------------
3543	! This routine calculates the view factors
3544	!------------------------------------------------------------------------
3545
3546	hut=z(nz_u+1)
3547
3548	do jz=1,nz_u
3549
3550	! radiation from wall to wall
3551
3552	do iz=1,nz_u
3553
3554	call fprls (fprl,dxy,abs(z(jz+1)-z(iz )),ws)
3555	f123=fprl
3556	call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
3557	f23=fprl
3558	call fprls (fprl,dxy,abs(z(jz )-z(iz )),ws)
3559	f12=fprl
3560	call fprls (fprl,dxy,abs(z(jz )-z(iz+1)),ws)
3561	f2 = fprl
3562
3563	a123=dxy*(abs(z(jz+1)-z(iz )))
3564	a12 =dxy*(abs(z(jz )-z(iz )))
3565	a23 =dxy*(abs(z(jz+1)-z(iz+1)))
3566	a1 =dxy*(abs(z(iz+1)-z(iz )))
3567	a2 =dxy*(abs(z(jz )-z(iz+1)))
3568	a3 =dxy*(abs(z(jz+1)-z(jz )))
3569
3570	ftot=0.5(a123f123-a23f23-a12f12+a2*f2)/a1
3571
3572	fww(iz,jz,id,iurb)=ftot*a1/a3
3573
3574	enddo
3575
3576	! radiation from ground to wall
3577
3578	call fnrms (fnrm,z(jz+1),dxy,ws)
3579	f12=fnrm
3580	call fnrms (fnrm,z(jz) ,dxy,ws)
3581	f1=fnrm
3582
3583	a1 = ws*dxy
3584
3585	a12= ws*dxy
3586
3587	a4=(z(jz+1)-z(jz))*dxy
3588
3589	ftot=(a12f12-a12f1)/a1
3590
3591	fgw(jz,id,iurb)=ftot*a1/a4
3592
3593	! radiation from sky to wall
3594
3595	call fnrms(fnrm,hut-z(jz) ,dxy,ws)
3596	f12 = fnrm
3597	call fnrms (fnrm,hut-z(jz+1),dxy,ws)
3598	f1 =fnrm
3599
3600	a1 = ws*dxy
3601
3602	a12= ws*dxy
3603
3604	a4 = (z(jz+1)-z(jz))*dxy
3605
3606	ftot=(a12f12-a12f1)/a1
3607
3608	fsw(jz,id,iurb)=ftot*a1/a4
3609
3610	enddo
3611
3612	! radiation from wall to sky
3613	do iz=1,nz_u
3614	call fnrms(fnrm,ws,dxy,hut-z(iz))
3615	f12=fnrm
3616	call fnrms(fnrm,ws,dxy,hut-z(iz+1))
3617	f1=fnrm
3618	a1 = (z(iz+1)-z(iz))*dxy
3619	a2 = (hut-z(iz+1))*dxy
3620	a12= (hut-z(iz))*dxy
3621	a4 = ws*dxy
3622	ftot=(a12f12-a2f1)/a1
3623	fws(iz,id,iurb)=ftot*a1/a4
3624
3625	enddo
3626	!!!!!!!!!!!!!
3627
3628
3629	do iz=1,nz_u
3630
3631	! radiation from wall to ground
3632
3633	call fnrms (fnrm,ws,dxy,z(iz+1))
3634	f12=fnrm
3635	call fnrms (fnrm,ws,dxy,z(iz ))
3636	f1 =fnrm
3637
3638	a1= (z(iz+1)-z(iz) )*dxy
3639
3640	a2 = z(iz)*dxy
3641	a12= z(iz+1)*dxy
3642	a4 = ws*dxy
3643
3644	ftot=(a12f12-a2f1)/a1
3645
3646	fwg(iz,id,iurb)=ftot*a1/a4
3647
3648	enddo
3649
3650	! radiation from sky to ground
3651
3652	call fprls (fprl,dxy,ws,hut)
3653	fsg(id,iurb)=fprl
3654
3655	return
3656	end subroutine view_factors
3657
3658	! ===6=8===============================================================72
3659	! ===6=8===============================================================72
3660
3661	SUBROUTINE fprls (fprl,a,b,c)
3662
3663	implicit none
3664
3665
3666
3667	real a,b,c
3668	real x,y
3669	real fprl
3670
3671
3672	x=a/c
3673	y=b/c
3674
3675	if(a.eq.0.or.b.eq.0.)then
3676	fprl=0.
3677	else
3678	fprl=log( ( (1.+x*2)(1.+y2)/(1.+x2+y2) ).5)+ &
3679	y((1.+x2).5)atan(y/((1.+x2).5))+ &
3680	x((1.+y2).5)atan(x/((1.+y2).5))- &
3681	yatan(y)-xatan(x)
3682	fprl=fprl2./(pix*y)
3683	endif
3684
3685	return
3686	end subroutine fprls
3687
3688	! ===6=8===============================================================72
3689	! ===6=8===============================================================72
3690
3691	SUBROUTINE fnrms (fnrm,a,b,c)
3692
3693	implicit none
3694
3695
3696
3697	real a,b,c
3698	real x,y,z,a1,a2,a3,a4,a5,a6
3699	real fnrm
3700
3701	x=a/b
3702	y=c/b
3703	z=x2+y2
3704
3705	if(y.eq.0.or.x.eq.0)then
3706	fnrm=0.
3707	else
3708	a1=log( (1.+xx)(1.+y*y)/(1.+z) )
3709	a2=yylog(yy(1.+z)/z/(1.+y*y) )
3710	a3=xxlog(xx(1.+z)/z/(1.+x*x) )
3711	a4=y*atan(1./y)
3712	a5=x*atan(1./x)
3713	a6=sqrt(z)*atan(1./sqrt(z))
3714	fnrm=0.25*(a1+a2+a3)+a4+a5-a6
3715	fnrm=fnrm/(pi*y)
3716	endif
3717
3718	return
3719	end subroutine fnrms
3720	! ===6=8===============================================================72
3721
3722	SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
3723	twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,albwin_u,emg_u,emw_u,&
3724	emr_u,emwind_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b, &
3725	cop_u, pwin_u, beta_u, sw_cond_u, time_on_u, time_off_u, &
3726	targtemp_u, gaptemp_u, targhum_u, gaphum_u, perflo_u, hsesf_u, hsequip)
3727
3728	! initialization routine, where the variables from the table are read
3729
3730	implicit none
3731
3732	integer iurb ! urban class number
3733	! Building parameters
3734	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
3735	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
3736	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
3737	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
3738	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
3739	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
3740	real twini_u(nurbm) ! Temperature inside the buildings behind the wall [K]
3741	real trini_u(nurbm) ! Temperature inside the buildings behind the roof [K]
3742	real tgini_u(nurbm) ! Initial road temperature
3743
3744	! Radiation parameters
3745	real albg_u(nurbm) ! Albedo of the ground
3746	real albw_u(nurbm) ! Albedo of the wall
3747	real albr_u(nurbm) ! Albedo of the roof
3748	real albwin_u(nurbm) ! Albedo of the window
3749	real emg_u(nurbm) ! Emissivity of ground
3750	real emw_u(nurbm) ! Emissivity of wall
3751	real emr_u(nurbm) ! Emissivity of roof
3752	real emwind_u(nurbm) ! Emissivity of windows
3753
3754	! Roughness parameters
3755	real z0g_u(nurbm) ! The ground's roughness length
3756	real z0r_u(nurbm) ! The roof's roughness length
3757
3758	! Street parameters
3759	integer nd_u(nurbm) ! Number of street direction for each urban class
3760
3761	real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
3762	real drst_u(ndm,nurbm) ! Street direction [degree]
3763	real ws_u(ndm,nurbm) ! Street width [m]
3764	real bs_u(ndm,nurbm) ! Building width [m]
3765	real h_b(nz_um,nurbm) ! Bulding's heights [m]
3766	real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
3767
3768	integer i,iu
3769	integer nurb ! number of urban classes used
3770
3771	real, intent(out) :: cop_u(nurbm)
3772	real, intent(out) :: pwin_u(nurbm)
3773	real, intent(out) :: beta_u(nurbm)
3774	integer, intent(out) :: sw_cond_u(nurbm)
3775	real, intent(out) :: time_on_u(nurbm)
3776	real, intent(out) :: time_off_u(nurbm)
3777	real, intent(out) :: targtemp_u(nurbm)
3778	real, intent(out) :: gaptemp_u(nurbm)
3779	real, intent(out) :: targhum_u(nurbm)
3780	real, intent(out) :: gaphum_u(nurbm)
3781	real, intent(out) :: perflo_u(nurbm)
3782	real, intent(out) :: hsesf_u(nurbm)
3783	real, intent(out) :: hsequip(24)
3784
3785	!
3786	!We initialize
3787	!
3788	h_b=0.
3789	d_b=0.
3790
3791	nurb=ICATE
3792	do iu=1,nurb
3793	nd_u(iu)=0
3794	enddo
3795
3796	csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
3797	csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
3798	csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
3799	do i=1,icate
3800	alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3801	alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3802	alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3803	enddo
3804	twini_u=TBLEND_TBL
3805	trini_u=TRLEND_TBL
3806	tgini_u=TGLEND_TBL
3807	albw_u=ALBB_TBL
3808	albr_u=ALBR_TBL
3809	albg_u=ALBG_TBL
3810	emw_u=EPSB_TBL
3811	emr_u=EPSR_TBL
3812	emg_u=EPSG_TBL
3813	z0r_u=Z0R_TBL
3814	z0g_u=Z0G_TBL
3815	nd_u=NUMDIR_TBL
3816	!MT BEM
3817	cop_u = cop_tbl
3818	pwin_u = pwin_tbl
3819	beta_u = beta_tbl
3820	sw_cond_u = sw_cond_tbl
3821	time_on_u = time_on_tbl
3822	time_off_u = time_off_tbl
3823	targtemp_u = targtemp_tbl
3824	gaptemp_u = gaptemp_tbl
3825	targhum_u = targhum_tbl
3826	gaphum_u = gaphum_tbl
3827	perflo_u = perflo_tbl
3828	hsesf_u = hsesf_tbl
3829	hsequip = hsequip_tbl
3830
3831	do iu=1,icate
3832	if(ndm.lt.nd_u(iu))then
3833	write(,)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
3834	write(,)'remember also that num_urban_layers should be equal or greater than nz_umndmnwr-u!'
3835	stop
3836	endif
3837	do i=1,nd_u(iu)
3838	drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
3839	ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
3840	bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
3841	enddo
3842	enddo
3843	do iu=1,ICATE
3844	if(nz_um.lt.numhgt_tbl(iu)+3)then
3845	write(,)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
3846	write(,)'remember also that num_urban_layers should be equal or greater than nz_umndmnwr-u!'
3847	stop
3848	endif
3849	do i=1,NUMHGT_TBL(iu)
3850	h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
3851	d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
3852	enddo
3853	enddo
3854
3855	do i=1,ndm
3856	do iu=1,nurbm
3857	strd_u(i,iu)=100000.
3858	enddo
3859	enddo
3860
3861	do iu=1,nurb
3862	emwind_u(iu)=0.9
3863	call albwindow(albwin_u(iu))
3864	enddo
3865
3866	return
3867	end subroutine init_para
3868	! ===6================================================================72
3869	! ===6================================================================72
3870	subroutine upward_rad(nd_u,nz_u,ws,bs,sigma,pb,ss, &
3871	tg,emg_u,albg_u,rlg,rsg,sfg, &
3872	tw,emw_u,albw_u,rlw,rsw,sfw, &
3873	tr,emr_u,albr_u,emwind,albwind,twlev,pwin, &
3874	sfwind,rld,rs, sfr, &
3875	rs_abs,rl_up,emiss,grdflx_urb)
3876	!
3877	! IN this surboutine we compute the upward longwave flux, and the albedo
3878	! needed for the radiation scheme
3879	!
3880	implicit none
3881
3882	!
3883	!INPUT VARIABLES
3884	!
3885	real rsw(2*ndm,nz_um) ! Short wave radiation at the wall for a given canyon direction [W/m2]
3886	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls for a given canyon direction [W/m2]
3887	real rsg(ndm) ! Short wave radiation at the canyon for a given canyon direction [W/m2]
3888	real rlg(ndm) ! Long wave radiation at the ground for a given canyon direction [W/m2]
3889	real rs ! Short wave radiation at the horizontal surface from the sun [W/mÃÂ²]
3890	real sfw(2*ndm,nz_um) ! Sensible heat flux from walls [W/mÃÂ²]
3891	real sfg(ndm) ! Sensible heat flux from ground (road) [W/mÃÂ²]
3892	real sfr(ndm,nz_um) ! Sensible heat flux from roofs [W/mÃÂ²]
3893	real rld ! Long wave radiation from the sky [W/mÃÂ²]
3894	real albg_u ! albedo of the ground/street
3895	real albw_u ! albedo of the walls
3896	real albr_u ! albedo of the roof
3897	real ws(ndm) ! width of the street
3898	real bs(ndm)
3899	! building size
3900	real pb(nz_um) ! Probability to have a building with an height equal or higher
3901	integer nz_u
3902	real ss(nz_um) ! Probability to have a building of a given height
3903	real sigma
3904	real emg_u ! emissivity of the street
3905	real emw_u ! emissivity of the wall
3906	real emr_u ! emissivity of the roof
3907	real tw(2*ndm,nz_um) ! Temperature in each layer of the wall [K]
3908	real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
3909	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
3910	integer id ! street direction
3911	integer nd_u ! number of street directions
3912	!
3913	!New variables BEM
3914	!
3915	real emwind !Emissivity of the windows
3916	real albwind !Albedo of the windows
3917	real twlev(2*ndm,nz_um) !Averaged Temperature of the windows
3918	real pwin !Coverage area fraction of the windows
3919	real gflwin !Heat stored for the windows
3920	real sfwind(2*ndm,nz_um) !Sensible heat flux from windows [W/mÃÂ²]
3921
3922	!OUTPUT/INPUT
3923	real rs_abs ! absrobed solar radiationfor this street direction
3924	real rl_up ! upward longwave radiation for this street direction
3925	real emiss ! mean emissivity
3926	real grdflx_urb ! ground heat flux
3927	!LOCAL
3928	integer iz,iw
3929	real rl_inc,rl_emit
3930	real gfl
3931	integer ix,iy,iwrong
3932
3933	iwrong=1
3934	do iz=1,nz_u+1
3935	do id=1,nd_u
3936	do iw=1,nwr_u
3937	if(tr(id,iz,iw).lt.100.)then
3938	write(,)'in upward_rad ',iz,id,iw,tr(id,iz,iw)
3939	iwrong=0
3940	endif
3941	enddo
3942	enddo
3943	enddo
3944	if(iwrong.eq.0)stop
3945
3946	rl_up=0.
3947
3948	rs_abs=0.
3949	rl_inc=0.
3950	emiss=0.
3951	rl_emit=0.
3952	grdflx_urb=0.
3953	do id=1,nd_u
3954	rl_emit=rl_emit-( emg_usigma(tg(id,ng_u)*4.)+(1-emg_u)rlg(id))*ws(id)/(ws(id)+bs(id))/nd_u
3955	rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u
3956	rs_abs=rs_abs+(1.-albg_u)rsg(id)ws(id)/(ws(id)+bs(id))/nd_u
3957	gfl=(1.-albg_u)rsg(id)+emg_urlg(id)-emg_usigma(tg(id,ng_u)**4.)+sfg(id)
3958	grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u
3959
3960	do iz=2,nz_u
3961	rl_emit=rl_emit-(emr_usigma(tr(id,iz,nwr_u)*4.)+(1-emr_u)rld)ss(iz)bs(id)/(ws(id)+bs(id))/nd_u
3962	rl_inc=rl_inc+rldss(iz)bs(id)/(ws(id)+bs(id))/nd_u
3963	rs_abs=rs_abs+(1.-albr_u)rsss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
3964	gfl=(1.-albr_u)rs+emr_urld-emr_usigma(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
3965	grdflx_urb=grdflx_urb-gflss(iz)bs(id)/(ws(id)+bs(id))/nd_u
3966	enddo
3967
3968	do iz=1,nz_u
3969
3970	rl_emit=rl_emit-(emw_u(1.-pwin)sigma(tw(2id-1,iz)*4.+tw(2id,iz)**4.)+ &
3971	(emwindpwinsigma(twlev(2id-1,iz)*4.+twlev(2id,iz)**4.))+ &
3972	((1.-emw_u)(1.-pwin)+pwin(1.-emwind))(rlw(2id-1,iz)+rlw(2id,iz))) &
3973	dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
3974
3975	rl_inc=rl_inc+((rlw(2id-1,iz)+rlw(2id,iz)))dz_upb(iz+1)/(ws(id)+bs(id))/nd_u
3976
3977	rs_abs=rs_abs+(((1.-albw_u)(1.-pwin)+(1.-albwind)pwin)(rsw(2id-1,iz)+rsw(2id,iz)))&
3978	dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
3979
3980	gfl=(1.-albw_u)(rsw(2id-1,iz)+rsw(2id,iz)) +emw_u( rlw(2id-1,iz)+rlw(2id,iz) ) &
3981	-emw_usigma( tw(2id-1,iz)4.+tw(2id,iz)*4. )+(sfw(2id-1,iz)+sfw(2*id,iz))
3982
3983	gflwin=(1.-albwind)(rsw(2id-1,iz)+rsw(2id,iz)) +emwind(rlw(2id-1,iz)+rlw(2id,iz)) &
3984	-emwindsigma( twlev(2id-1,iz)4.+twlev(2id,iz)*4.)+(sfwind(2id-1,iz)+sfwind(2*id,iz))
3985
3986
3987	grdflx_urb=grdflx_urb-(gfl(1.-pwin)+pwingflwin)dz_upb(iz+1)/(ws(id)+bs(id))/nd_u
3988
3989	enddo
3990
3991	enddo
3992	emiss=(emg_u+emw_u+emr_u)/3.
3993	rl_up=(rl_inc+rl_emit)-rld
3994
3995
3996	return
3997
3998	END SUBROUTINE upward_rad
3999
4000	!====6=8===============================================================72
4001	!====6=8===============================================================72
4002	! ===6================================================================72
4003	! ===6================================================================72
4004
4005	subroutine albwindow(albwin)
4006
4007	!-------------------------------------------------------------------
4008	implicit none
4009
4010
4011	! -------------------------------------------------------------------
4012	!Based on the
4013	!paper of J.Karlsson and A.Roos(2000):"modelling the angular behaviour
4014	!of the total solar energy transmittance of windows"
4015	!Solar Energy Vol.69,No.4,pp. 321-329.
4016	! -------------------------------------------------------------------
4017	!Input
4018	!-----
4019
4020	!Output
4021	!------
4022	real albwin ! albedo of the window
4023	!Local
4024	!-----
4025	real a,b,c !Polynomial coefficients
4026	real alfa,delta,gama !Polynomial powers
4027	real g0 !transmittance when the angle
4028	!of incidence is normal to the surface.
4029	real asup,ainf
4030	real fonc
4031
4032	!Constants
4033	!--------------------
4034
4035	real epsilon !accuracy of the integration
4036	parameter (epsilon=1.e-07)
4037	real n1,n2 !Index of refraction for glasses and air
4038	parameter(n1=1.,n2=1.5)
4039	integer intg,k
4040	!--------------------------------------------------------------------
4041	if (q_num.eq.0) then
4042	write(,) 'Category parameter of the windows no valid'
4043	stop
4044	endif
4045
4046	g0=4.n1n2/((n1+n2)*(n1+n2))
4047	a=8.
4048	b=0.25/q_num
4049	c=1.-a-b
4050	alfa =5.2 + (0.7*q_num)
4051	delta =2.
4052	gama =(5.26+0.06p_num)+(0.73+0.04p_num)*q_num
4053
4054	intg=1
4055	!----------------------------------------------------------------------
4056
4057
4058	100 asup=0.
4059	ainf=0.
4060
4061	do k=1,intg
4062	call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
4063	asup=asup+(pi/intg)*fonc
4064	enddo
4065
4066	intg=intg+1
4067
4068	do k=1,intg
4069	call foncs(fonc,(pi*k/intg),a,b,c,alfa,delta,gama)
4070	ainf=ainf+(pi/intg)*fonc
4071	enddo
4072
4073	if(abs(asup-ainf).lt.epsilon) then
4074	albwin=1-g0+(g0/2.)*asup
4075	else
4076	goto 100
4077	endif
4078
4079	!----------------------------------------------------------------------
4080	return
4081	end subroutine albwindow
4082	!====================================================================72
4083	!====================================================================72
4084
4085	subroutine foncs(fonc,x,aa,bb,cc,alf,delt,gam)
4086
4087	implicit none
4088	!
4089	real x,aa,bb,cc
4090	real alf,delt,gam
4091	real fonc
4092
4093	fonc=(((aa(xalf))/(pi*alf))+ &
4094	((bb(xdelt))/(pi*delt))+ &
4095	((cc(xgam))/(pigam)))sin(x)
4096
4097	return
4098	end subroutine foncs
4099	!====================================================================72
4100	!====================================================================72
4101	END MODULE module_sf_bep_bem

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format