Context Navigation

module_sf_bep.F

Last change on this file 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: 126.0 KB

Line
1	MODULE module_sf_bep
2
3	!USE module_model_constants
4	USE module_sf_urban
5
6	! SGClarke 09/11/2008
7	! Access urban_param.tbl values through calling urban_param_init in module_physics_init
8	! for CASE (BEPSCHEME) select sf_urban_physics
9	!
10	! -----------------------------------------------------------------------
11	! Dimension for the array used in the BEP module
12	! -----------------------------------------------------------------------
13
14	integer nurbm ! Maximum number of urban classes
15	parameter (nurbm=3)
16
17	integer ndm ! Maximum number of street directions
18	parameter (ndm=2)
19
20	integer nz_um ! Maximum number of vertical levels in the urban grid
21	parameter(nz_um=13)
22
23	integer ng_u ! Number of grid levels in the ground
24	parameter (ng_u=10)
25	integer nwr_u ! Number of grid levels in the walls or roofs
26	parameter (nwr_u=10)
27
28	real dz_u ! Urban grid resolution
29	parameter (dz_u=5.)
30
31	! The change of ng_u, nwr_u should be done in agreement with the block data
32	! in the routine "surf_temp"
33	! -----------------------------------------------------------------------
34	! Constant used in the BEP module
35	! -----------------------------------------------------------------------
36
37	real vk ! von Karman constant
38	real g_u ! Gravity acceleration
39	real pi !
40	real r ! Perfect gas constant
41	real cp_u ! Specific heat at constant pressure
42	real rcp_u !
43	real sigma !
44	real p0 ! Reference pressure at the sea level
45	real cdrag ! Drag force constant
46	parameter(vk=0.40,g_u=9.81,pi=3.141592653,r=287.,cp_u=1004.)
47	parameter(rcp_u=r/cp_u,sigma=5.67e-08,p0=1.e+5,cdrag=0.4)
48
49	! -----------------------------------------------------------------------
50
51
52
53
54	CONTAINS
55
56	subroutine BEP(FRC_URB2D,UTYPE_URB2D,itimestep,dz8w,dt,u_phy,v_phy, &
57	th_phy,rho,p_phy,swdown,glw, &
58	gmt,julday,xlong,xlat, &
59	declin_urb,cosz_urb2d,omg_urb2d, &
60	num_urban_layers, &
61	trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
62	sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
63	a_u,a_v,a_t,a_e,b_u,b_v, &
64	b_t,b_e,dlg,dl_u,sf,vl, &
65	rl_up,rs_abs,emiss,grdflx_urb, &
66	ids,ide, jds,jde, kds,kde, &
67	ims,ime, jms,jme, kms,kme, &
68	its,ite, jts,jte, kts,kte)
69
70	implicit none
71
72	!------------------------------------------------------------------------
73	! Input
74	!------------------------------------------------------------------------
75	INTEGER :: ids,ide, jds,jde, kds,kde, &
76	ims,ime, jms,jme, kms,kme, &
77	its,ite, jts,jte, kts,kte, &
78	itimestep
79
80
81	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: DZ8W
82	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: P_PHY
83	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: RHO
84	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: TH_PHY
85	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: T_PHY
86	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U_PHY
87	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V_PHY
88	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: U
89	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ):: V
90	REAL, DIMENSION( ims:ime , jms:jme ) :: GLW
91	REAL, DIMENSION( ims:ime , jms:jme ) :: swdown
92	REAL, DIMENSION( ims:ime, jms:jme ) :: UST
93	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: UTYPE_URB2D
94	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: FRC_URB2D
95	REAL, INTENT(IN ) :: GMT
96	INTEGER, INTENT(IN ) :: JULDAY
97	REAL, DIMENSION( ims:ime, jms:jme ), &
98	INTENT(IN ) :: XLAT, XLONG
99	REAL, INTENT(IN) :: DECLIN_URB
100	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
101	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
102	INTEGER, INTENT(IN ) :: num_urban_layers
103	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
104	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
105	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
106	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
107	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
108	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
109	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
110	REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
111	! integer nx,ny,nz ! Number of points in the mesocsale grid
112	real z(ims:ime,kms:kme,jms:jme) ! Vertical coordinates
113	REAL, INTENT(IN ):: DT ! Time step
114	! real zr(ims:ime,jms:jme) ! Solar zenith angle
115	! real deltar(ims:ime,jms:jme) ! Declination of the sun
116	! real ah(ims:ime,jms:jme) ! Hour angle
117	! real rs(ims:ime,jms:jme) ! Solar radiation
118	!------------------------------------------------------------------------
119	! Output
120	!------------------------------------------------------------------------
121	! real tsk(ims:ime,jms:jme) ! Average of surface temperatures (roads and roofs)
122	!
123	! Implicit and explicit components of the source and sink terms at each levels,
124	! the fluxes can be computed as follow: FX = AX + B example: t_fluxes = a_t pt + b_t
125	real a_u(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in X-direction (center)
126	real a_v(ims:ime,kms:kme,jms:jme) ! Implicit component for the momemtum in Y-direction (center)
127	real a_t(ims:ime,kms:kme,jms:jme) ! Implicit component for the temperature
128	real a_e(ims:ime,kms:kme,jms:jme) ! Implicit component for the TKE
129	real b_u(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in X-direction (center)
130	real b_v(ims:ime,kms:kme,jms:jme) ! Explicit component for the momemtum in Y-direction (center)
131	real b_t(ims:ime,kms:kme,jms:jme) ! Explicit component for the temperature
132	real b_e(ims:ime,kms:kme,jms:jme) ! Explicit component for the TKE
133	real dlg(ims:ime,kms:kme,jms:jme) ! Height above ground (L_ground in formula (24) of the BLM paper).
134	real dl_u(ims:ime,kms:kme,jms:jme) ! Length scale (lb in formula (22) ofthe BLM paper).
135	! urban surface and volumes
136	real sf(ims:ime,kms:kme,jms:jme) ! surface of the urban grid cells
137	real vl(ims:ime,kms:kme,jms:jme) ! volume of the urban grid cells
138	! urban fluxes
139	real rl_up(ims:ime,jms:jme) ! upward long wave radiation
140	real rs_abs(ims:ime,jms:jme) ! absorbed short wave radiation
141	real emiss(ims:ime,jms:jme) ! emissivity averaged for urban surfaces
142	real grdflx_urb(ims:ime,jms:jme) ! ground heat flux for urban areas
143	!------------------------------------------------------------------------
144	! Local
145	!------------------------------------------------------------------------
146	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
147	! Building parameters
148	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
149	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
150	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
151	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
152	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
153	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
154	real twini_u(nurbm) ! Initial temperature inside the building's wall [K]
155	real trini_u(nurbm) ! Initial temperature inside the building's roof [K]
156	real tgini_u(nurbm) ! Initial road temperature
157	!
158	! for twini_u, and trini_u the initial value at the deepest level is kept constant during the simulation
159	!
160	! Radiation paramters
161	real albg_u(nurbm) ! Albedo of the ground
162	real albw_u(nurbm) ! Albedo of the wall
163	real albr_u(nurbm) ! Albedo of the roof
164	real emg_u(nurbm) ! Emissivity of ground
165	real emw_u(nurbm) ! Emissivity of wall
166	real emr_u(nurbm) ! Emissivity of roof
167
168	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
169	! and the short wave radation.
170	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
171	real fwg(nz_um,ndm,nurbm) ! from wall to ground
172	real fgw(nz_um,ndm,nurbm) ! from ground to wall
173	real fsw(nz_um,ndm,nurbm) ! from sky to wall
174	real fws(nz_um,ndm,nurbm) ! from sky to wall
175	real fsg(ndm,nurbm) ! from sky to ground
176
177	! Roughness parameters
178	real z0g_u(nurbm) ! The ground's roughness length
179	real z0r_u(nurbm) ! The roof's roughness length
180
181	! Street parameters
182	integer nd_u(nurbm) ! Number of street direction for each urban class
183	real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
184	real drst_u(ndm,nurbm) ! Street direction
185	real ws_u(ndm,nurbm) ! Street width
186	real bs_u(ndm,nurbm) ! Building width
187	real h_b(nz_um,nurbm) ! Bulding's heights
188	real d_b(nz_um,nurbm) ! Probability that a building has an height h_b
189	real ss_u(nz_um,nurbm) ! Probability that a building has an height equal to z
190	real pb_u(nz_um,nurbm) ! Probability that a building has an height greater or equal to z
191
192
193	! Grid parameters
194	integer nz_u(nurbm) ! Number of layer in the urban grid
195
196	real z_u(nz_um) ! Height of the urban grid levels
197
198	! 1D array used for the input and output of the routine "urban"
199
200	real z1D(kms:kme) ! vertical coordinates
201	real ua1D(kms:kme) ! wind speed in the x directions
202	real va1D(kms:kme) ! wind speed in the y directions
203	real pt1D(kms:kme) ! potential temperature
204	real da1D(kms:kme) ! air density
205	real pr1D(kms:kme) ! air pressure
206	real pt01D(kms:kme) ! reference potential temperature
207	real zr1D ! zenith angle
208	real deltar1D ! declination of the sun
209	real ah1D ! hour angle (it should come from the radiation routine)
210	real rs1D ! solar radiation
211	real rld1D ! downward flux of the longwave radiation
212
213
214	real tw1D(2*ndm,nz_um,nwr_u) ! temperature in each layer of the wall
215	real tg1D(ndm,ng_u) ! temperature in each layer of the ground
216	real tr1D(ndm,nz_um,nwr_u) ! temperature in each layer of the roof
217	real sfw1D(2*ndm,nz_um) ! sensible heat flux from walls
218	real sfg1D(ndm) ! sensible heat flux from ground (road)
219	real sfr1D(ndm,nz_um) ! sensible heat flux from roofs
220	real sf1D(kms:kme) ! surface of the urban grid cells
221	real vl1D(kms:kme) ! volume of the urban grid cells
222	real a_u1D(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
223	real a_v1D(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
224	real a_t1D(kms:kme) ! Implicit component of the heat sources or sinks
225	real a_e1D(kms:kme) ! Implicit component of the TKE sources or sinks
226	real b_u1D(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
227	real b_v1D(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
228	real b_t1D(kms:kme) ! Explicit component of the heat sources or sinks
229	real b_e1D(kms:kme) ! Explicit component of the TKE sources or sinks
230	real dlg1D(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
231	real dl_u1D(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper)
232	real tsk1D ! Average of the road surface temperatures
233	real time_bep
234	! arrays used to collapse indexes
235	integer ind_zwd(nz_um,nwr_u,ndm)
236	integer ind_gd(ng_u,ndm)
237	integer ind_zd(nz_um,ndm)
238	!
239	integer ix,iy,iz,iurb,id,iz_u,iw,ig,ir,ix1,iy1,k
240	integer it, nint
241	integer iii
242	real time_h,tempo,shtot
243	logical first
244	character(len=80) :: text
245	data first/.true./
246	save first,time_bep
247
248	save alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
249	albg_u,albw_u,albr_u,emg_u,emw_u,emr_u,fww,fwg,fgw,fsw,fws,fsg, &
250	z0g_u,z0r_u, nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
251	nz_u,z_u
252
253	!------------------------------------------------------------------------
254	! Calculation of the momentum, heat and turbulent kinetic fluxes
255	! produced by builgings
256	!
257	! Reference:
258	! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
259	! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
260	! 261-304
261	!------------------------------------------------------------------------
262	!prepare the arrays to collapse indexes
263
264	if(num_urban_layers.lt.nz_umndmnwr_u)then
265	write(,)'num_urban_layers too small, please increase to at least ', nz_umndmnwr_u
266	stop
267	endif
268	iii=0
269	do iz_u=1,nz_um
270	do iw=1,nwr_u
271	do id=1,ndm
272	iii=iii+1
273	ind_zwd(iz_u,iw,id)=iii
274	enddo
275	enddo
276	enddo
277
278	iii=0
279	do ig=1,ng_u
280	do id=1,ndm
281	iii=iii+1
282	ind_gd(ig,id)=iii
283	enddo
284	enddo
285
286	iii=0
287	do iz_u=1,nz_um
288	do id=1,ndm
289	iii=iii+1
290	ind_zd(iz_u,id)=iii
291	enddo
292	enddo
293	do ix=its,ite
294	do iy=jts,jte
295	z(ix,kts,iy)=0.
296	do iz=kts+1,kte+1
297	z(ix,iz,iy)=z(ix,iz-1,iy)+dz8w(ix,iz-1,iy)
298	enddo
299	enddo
300	enddo
301
302	if (first) then ! True only on first call
303	call init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
304	twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
305	emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
306
307	! Initialisation of the urban parameters and calculation of the view factors
308	call icBEP(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
309	albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
310	fww,fwg,fgw,fsw,fws,fsg, &
311	z0g_u,z0r_u, &
312	nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
313	nz_u,z_u, &
314	twini_u,trini_u)
315
316	first=.false.
317
318	endif ! first
319
320	do ix=its,ite
321	do iy=jts,jte
322	if (FRC_URB2D(ix,iy).gt.0.) then ! Calling BEP only for existing urban classes.
323
324	iurb=UTYPE_URB2D(ix,iy)
325
326	do iz= kts,kte
327	ua1D(iz)=u_phy(ix,iz,iy)
328	va1D(iz)=v_phy(ix,iz,iy)
329	pt1D(iz)=th_phy(ix,iz,iy)
330	da1D(iz)=rho(ix,iz,iy)
331	pr1D(iz)=p_phy(ix,iz,iy)
332	! pt01D(iz)=th_phy(ix,iz,iy)
333	pt01D(iz)=300.
334	z1D(iz)=z(ix,iz,iy)
335	a_u1D(iz)=0.
336	a_v1D(iz)=0.
337	a_t1D(iz)=0.
338	a_e1D(iz)=0.
339	b_u1D(iz)=0.
340	b_v1D(iz)=0.
341	b_t1D(iz)=0.
342	b_e1D(iz)=0.
343	enddo
344	z1D(kte+1)=z(ix,kte+1,iy)
345
346	do id=1,ndm
347	do iz_u=1,nz_um
348	do iw=1,nwr_u
349	! tw1D(2*id-1,iz_u,iw)=tw1_u(ix,iy,ind_zwd(iz_u,iw,id))
350	! tw1D(2*id,iz_u,iw)=tw2_u(ix,iy,ind_zwd(iz_u,iw,id))
351	if(ind_zwd(iz_u,iw,id).gt.num_urban_layers)write(,)'ind_zwd too big w',ind_zwd(iz_u,iw,id)
352	tw1D(2*id-1,iz_u,iw)=tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
353	tw1D(2*id,iz_u,iw)=tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)
354	enddo
355	enddo
356	enddo
357
358	do id=1,ndm
359	do ig=1,ng_u
360	! tg1D(id,ig)=tg_u(ix,iy,ind_gd(ig,id))
361	tg1D(id,ig)=tgb_urb4d(ix,ind_gd(ig,id),iy)
362	enddo
363	do iz_u=1,nz_um
364	do ir=1,nwr_u
365	! tr1D(id,iz_u,ir)=tr_u(ix,iy,ind_zwd(iz_u,ir,id))
366	if(ind_zwd(iz_u,ir,id).gt.num_urban_layers)write(,)'ind_zwd too big r',ind_zwd(iz_u,ir,id)
367	tr1D(id,iz_u,ir)=trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)
368	enddo
369	enddo
370	enddo
371
372	do id=1,ndm
373	do iz=1,nz_um
374	! sfw1D(2*id-1,iz)=sfw1(ix,iy,ind_zd(iz,id))
375	! sfw1D(2*id,iz)=sfw2(ix,iy,ind_zd(iz,id))
376	sfw1D(2*id-1,iz)=sfw1_urb3d(ix,ind_zd(iz,id),iy)
377	sfw1D(2*id,iz)=sfw2_urb3d(ix,ind_zd(iz,id),iy)
378	enddo
379	enddo
380
381	do id=1,ndm
382	! sfg1D(id)=sfg(ix,iy,id)
383	sfg1D(id)=sfg_urb3d(ix,id,iy)
384	enddo
385
386	do id=1,ndm
387	do iz=1,nz_um
388	! sfr1D(id,iz)=sfr(ix,iy,ind_zd(iz,id))
389	sfr1D(id,iz)=sfr_urb3d(ix,ind_zd(iz,id),iy)
390	enddo
391	enddo
392
393
394	rs1D=swdown(ix,iy)
395	rld1D=glw(ix,iy)
396	time_h=(itimestep*dt)/3600.+gmt
397
398	zr1D=acos(COSZ_URB2D(ix,iy))
399	deltar1D=DECLIN_URB
400	ah1D=OMG_URB2D(ix,iy)
401	! call angle(xlong(ix,iy),xlat(ix,iy),julday,time_h,zr1D,deltar1D,ah1D)
402
403	call BEP1D(iurb,kms,kme,kts,kte,z1D,dt,ua1D,va1D,pt1D,da1D,pr1D,pt01D, &
404	zr1D,deltar1D,ah1D,rs1D,rld1D, &
405	alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
406	albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
407	fww,fwg,fgw,fsw,fws,fsg, &
408	z0g_u,z0r_u, &
409	nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
410	nz_u,z_u, &
411	tw1D,tg1D,tr1D,sfw1D,sfg1D,sfr1D, &
412	a_u1D,a_v1D,a_t1D,a_e1D, &
413	b_u1D,b_v1D,b_t1D,b_e1D, &
414	dlg1D,dl_u1D,tsk1D,sf1D,vl1D,rl_up(ix,iy), &
415	rs_abs(ix,iy),emiss(ix,iy),grdflx_urb(ix,iy))
416
417	do id=1,ndm
418	do iz=1,nz_um
419	sfw1_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id-1,iz)
420	sfw2_urb3d(ix,ind_zd(iz,id),iy)=sfw1D(2*id,iz)
421	enddo
422	enddo
423
424	do id=1,ndm
425	sfg_urb3d(ix,id,iy)=sfg1D(id)
426	enddo
427
428	do id=1,ndm
429	do iz=1,nz_um
430	sfr_urb3d(ix,ind_zd(iz,id),iy)=sfr1D(id,iz)
431	enddo
432	enddo
433	!
434	do id=1,ndm
435	do iz_u=1,nz_um
436	do iw=1,nwr_u
437	tw1_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id-1,iz_u,iw)
438	tw2_urb4d(ix,ind_zwd(iz_u,iw,id),iy)=tw1D(2*id,iz_u,iw)
439	enddo
440	enddo
441	enddo
442
443	do id=1,ndm
444	do ig=1,ng_u
445	tgb_urb4d(ix,ind_gd(ig,id),iy)=tg1D(id,ig)
446	enddo
447	do iz_u=1,nz_um
448	do ir=1,nwr_u
449	trb_urb4d(ix,ind_zwd(iz_u,ir,id),iy)=tr1D(id,iz_u,ir)
450	enddo
451	enddo
452	enddo
453
454	do iz= kts,kte
455	sf(ix,iz,iy)=sf1D(iz)
456	vl(ix,iz,iy)=vl1D(iz)
457	a_u(ix,iz,iy)=a_u1D(iz)
458	a_v(ix,iz,iy)=a_v1D(iz)
459	a_t(ix,iz,iy)=a_t1D(iz)
460	a_e(ix,iz,iy)=a_e1D(iz)
461	b_u(ix,iz,iy)=b_u1D(iz)
462	b_v(ix,iz,iy)=b_v1D(iz)
463	b_t(ix,iz,iy)=b_t1D(iz)
464	b_e(ix,iz,iy)=b_e1D(iz)
465	dlg(ix,iz,iy)=dlg1D(iz)
466	dl_u(ix,iz,iy)=dl_u1D(iz)
467	enddo
468	sf(ix,kte+1,iy)=sf1D(kte+1)
469	! tsk(ix,iy)=tsk1D
470	!
471	endif ! FRC_URB2D
472
473	enddo ! iy
474	enddo ! ix
475
476
477	time_bep=time_bep+dt
478
479
480	return
481	end subroutine BEP
482
483	! ===6=8===============================================================72
484
485	subroutine BEP1D(iurb,kms,kme,kts,kte,z,dt,ua,va,pt,da,pr,pt0, &
486	zr,deltar,ah,rs,rld, &
487	alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
488	albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
489	fww,fwg,fgw,fsw,fws,fsg, &
490	z0g_u,z0r_u, &
491	nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
492	nz_u,z_u, &
493	tw,tg,tr,sfw,sfg,sfr, &
494	a_u,a_v,a_t,a_e, &
495	b_u,b_v,b_t,b_e, &
496	dlg,dl_u,tsk,sf,vl,rl_up,rs_abs,emiss,grdflx_urb)
497
498	! ----------------------------------------------------------------------
499	! This routine computes the effects of buildings on momentum, heat and
500	! TKE (turbulent kinetic energy) sources or sinks and on the mixing length.
501	! It provides momentum, heat and TKE sources or sinks at different levels of a
502	! mesoscale grid defined by the altitude of its cell interfaces "z" and
503	! its number of levels "nz".
504	! The meteorological input parameters (wind, temperature, solar radiation)
505	! are specified on the "mesoscale grid".
506	! The inputs concerning the building and street charateristics are defined
507	! on a "urban grid". The "urban grid" is defined with its number of levels
508	! "nz_u" and its space step "dz_u".
509	! The input parameters are interpolated on the "urban grid". The sources or sinks
510	! are calculated on the "urban grid". Finally the sources or sinks are
511	! interpolated on the "mesoscale grid".
512
513
514	! Mesoscale grid Urban grid Mesoscale grid
515	!
516	! z(4) --- ---
517	! \| \|
518	! \| \|
519	! \| Interpolation Interpolation \|
520	! \| Sources or sinks calculation \|
521	! z(3) --- ---
522	! \| ua ua_u --- uv_a a_u \|
523	! \| va va_u \| uv_b b_u \|
524	! \| pt pt_u --- uh_b a_v \|
525	! z(2) --- \| etc... etc...---
526	! \| z_u(1) --- \|
527	! \| \| \|
528	! z(1) ------------------------------------------------------------
529
530	!
531	! Reference:
532	! Martilli, A., Clappier, A., Rotach, M.W.:2002, 'AN URBAN SURFACE EXCHANGE
533	! PARAMETERISATION FOR MESOSCALE MODELS', Boundary-Layer Meteorolgy 104:
534	! 261-304
535
536	! ----------------------------------------------------------------------
537
538	implicit none
539
540
541
542	! ----------------------------------------------------------------------
543	! INPUT:
544	! ----------------------------------------------------------------------
545
546	! Data relative to the "mesoscale grid"
547
548	! integer nz ! Number of vertical levels
549	integer kms,kme,kts,kte
550	real z(kms:kme) ! Altitude above the ground of the cell interfaces.
551	real ua(kms:kme) ! Wind speed in the x direction
552	real va(kms:kme) ! Wind speed in the y direction
553	real pt(kms:kme) ! Potential temperature
554	real da(kms:kme) ! Air density
555	real pr(kms:kme) ! Air pressure
556	real pt0(kms:kme) ! Reference potential temperature (could be equal to "pt")
557	real dt ! Time step
558	real zr ! Zenith angle
559	real deltar ! Declination of the sun
560	real ah ! Hour angle
561	real rs ! Solar radiation
562	real rld ! Downward flux of the longwave radiation
563
564	! Data relative to the "urban grid"
565
566	integer iurb ! Current urban class
567
568	! Building parameters
569	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
570	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
571	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
572	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
573	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
574	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
575
576	! Radiation parameters
577	real albg_u(nurbm) ! Albedo of the ground
578	real albw_u(nurbm) ! Albedo of the wall
579	real albr_u(nurbm) ! Albedo of the roof
580	real emg_u(nurbm) ! Emissivity of ground
581	real emw_u(nurbm) ! Emissivity of wall
582	real emr_u(nurbm) ! Emissivity of roof
583
584	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long and
585	! short wave radation.
586	! The calculation of these factor is explained in the Appendix A of the BLM paper
587	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
588	real fwg(nz_um,ndm,nurbm) ! from wall to ground
589	real fgw(nz_um,ndm,nurbm) ! from ground to wall
590	real fsw(nz_um,ndm,nurbm) ! from sky to wall
591	real fws(nz_um,ndm,nurbm) ! from wall to sky
592	real fsg(ndm,nurbm) ! from sky to ground
593
594	! Roughness parameters
595	real z0g_u(nurbm) ! The ground's roughness length
596	real z0r_u(nurbm) ! The roof's roughness length
597
598	! Street parameters
599	integer nd_u(nurbm) ! Number of street direction for each urban class
600	real strd_u(ndm,nurbm) ! Street length (set to a greater value then the horizontal length of the cells)
601	real drst_u(ndm,nurbm) ! Street direction
602	real ws_u(ndm,nurbm) ! Street width
603	real bs_u(ndm,nurbm) ! Building width
604	real h_b(nz_um,nurbm) ! Bulding's heights
605	real d_b(nz_um,nurbm) ! The probability that a building has an height "h_b"
606	real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
607	real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
608
609	! Grid parameters
610	integer nz_u(nurbm) ! Number of layer in the urban grid
611	! real dz_u ! Urban grid resolution
612	real z_u(nz_um) ! Height of the urban grid levels
613
614
615	! ----------------------------------------------------------------------
616	! INPUT-OUTPUT
617	! ----------------------------------------------------------------------
618
619	! Data relative to the "urban grid" which should be stored from the current time step to the next one
620
621	real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
622	real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
623	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
624	real sfw(2*ndm,nz_um) ! Sensible heat flux from walls
625	real sfg(ndm) ! Sensible heat flux from ground (road)
626	real sfr(ndm,nz_um) ! Sensible heat flux from roofs
627	real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) towards the interior
628	real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof towards the interior
629	real gfw(2*ndm,nz_um) ! Heat flux transfered from the surface of the walls towards the interior
630	! ----------------------------------------------------------------------
631	! OUTPUT:
632	! ----------------------------------------------------------------------
633
634	! Data relative to the "mesoscale grid"
635
636	real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
637	real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
638
639	! Implicit and explicit components of the source and sink terms at each levels,
640	! the fluxes can be computed as follow: FX = AX + B example: Heat fluxes = a_t pt + b_t
641	real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
642	real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
643	real a_t(kms:kme) ! Implicit component of the heat sources or sinks
644	real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
645	real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
646	real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
647	real b_t(kms:kme) ! Explicit component of the heat sources or sinks
648	real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
649	real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
650	real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
651	real tsk ! Average of the road surface temperatures
652
653	! ----------------------------------------------------------------------
654	! LOCAL:
655	! ----------------------------------------------------------------------
656
657	real dz(kms:kme) ! vertical space steps of the "mesoscale grid"
658
659	! Data interpolated from the "mesoscale grid" to the "urban grid"
660
661	real ua_u(nz_um) ! Wind speed in the x direction
662	real va_u(nz_um) ! Wind speed in the y direction
663	real pt_u(nz_um) ! Potential temperature
664	real da_u(nz_um) ! Air density
665	real pt0_u(nz_um) ! Reference potential temperature
666	real pr_u(nz_um) ! Air pressure
667
668	! Data defining the building and street charateristics
669
670	integer nd ! Number of street direction for the current urban class
671
672	real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
673	real alar(nwr_u) ! Roof thermal diffusivity for the current urban class [m^2 s^-1]
674	real alaw(nwr_u) ! Walls thermal diffusivity for the current urban class [m^2 s^-1]
675	real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
676	real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
677	real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
678
679	real z0(ndm,nz_um) ! Roughness lengths "profiles"
680	real ws(ndm) ! Street widths of the current urban class
681	real bs(ndm) ! Building widths of the current urban class
682	real strd(ndm) ! Street lengths for the current urban class
683	real drst(ndm) ! Street directions for the current urban class
684	real ss(nz_um) ! Probability to have a building with height h
685	real pb(nz_um) ! Probability to have a building with an height equal
686
687	! Solar radiation at each level of the "urban grid"
688
689	real rsg(ndm) ! Short wave radiation from the ground
690	real rsw(2*ndm,nz_um) ! Short wave radiation from the walls
691	real rlg(ndm) ! Long wave radiation from the ground
692	real rlw(2*ndm,nz_um) ! Long wave radiation from the walls
693
694	! Potential temperature of the surfaces at each level of the "urban grid"
695
696	real ptg(ndm) ! Ground potential temperatures
697	real ptr(ndm,nz_um) ! Roof potential temperatures
698	real ptw(2*ndm,nz_um) ! Walls potential temperatures
699
700
701	! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
702	! vertical surfaces (walls) ans horizontal surfaces (roofs and street)
703	! The fluxes can be computed as follow: Fluxes of X = A*X + B
704	! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
705
706	real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
707	real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
708	real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
709	real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
710	real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
711	real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
712	real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
713	real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
714	real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
715	real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
716	real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
717
718	!
719	real rs_abs ! solar radiation absorbed by urban surfaces
720	real rl_up ! longwave radiation emitted by urban surface to the atmosphere
721	real emiss ! mean emissivity of the urban surface
722	real grdflx_urb ! ground heat flux
723	real shtot,aaa
724	real dt_int ! internal time step
725	integer nt_int ! number of internal time step
726	integer iz,id, it_int
727	integer iwrong,iw,ix,iy
728
729	! ----------------------------------------------------------------------
730	! END VARIABLES DEFINITIONS
731	! ----------------------------------------------------------------------
732
733	! Fix some usefull parameters for the computation of the sources or sinks
734
735	do iz=kts,kte
736	dz(iz)=z(iz+1)-z(iz)
737	end do
738	call param(iurb,nz_u(iurb),nd_u(iurb), &
739	csg_u,csg,alag_u,alag,csr_u,csr, &
740	alar_u,alar,csw_u,csw,alaw_u,alaw, &
741	ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
742	strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb)
743
744	! Interpolation on the "urban grid"
745	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,ua,ua_u)
746	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,va,va_u)
747	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt,pt_u)
748	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pt0,pt0_u)
749	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,pr,pr_u)
750	call interpol(kms,kme,kts,kte,nz_u(iurb),z,z_u,da,da_u)
751
752
753	! Compute the modification of the radiation due to the buildings
754
755	call modif_rad(iurb,nd_u(iurb),nz_u(iurb),z_u,ws, &
756	drst,strd,ss,pb, &
757	tw,tg,albg_u(iurb),albw_u(iurb), &
758	emw_u(iurb),emg_u(iurb), &
759	fww,fwg,fgw,fsw,fsg, &
760	zr,deltar,ah, &
761	rs,rld,rsw,rsg,rlw,rlg)
762
763	! calculation of the urban albedo and the upward long wave radiation
764	call upward_rad(nd_u(iurb),iurb,nz_u(iurb),ws,bs,sigma,fsw,fsg,pb,ss, &
765	tg,emg_u(iurb),albg_u(iurb),rlg,rsg,sfg, &
766	tw,emw_u(iurb),albw_u(iurb),rlw,rsw,sfw, &
767	tr,emr_u(iurb),albr_u(iurb),rld,rs,sfr, &
768	rs_abs,rl_up,emiss,grdflx_urb)
769
770	! Compute the surface temperatures
771
772
773	call surf_temp(nz_u(iurb),nd_u(iurb),pr_u,dt,ss, &
774	rs,rld,rsg,rlg,rsw,rlw, &
775	tg,alag,csg,emg_u(iurb),albg_u(iurb),ptg,sfg,gfg, &
776	tr,alar,csr,emr_u(iurb),albr_u(iurb),ptr,sfr,gfr, &
777	tw,alaw,csw,emw_u(iurb),albw_u(iurb),ptw,sfw,gfw)
778
779
780	! Compute the implicit and explicit components of the sources or sinks on the "urban grid"
781
782	call buildings(nd_u(iurb),nz_u(iurb),z0,ua_u,va_u, &
783	pt_u,pt0_u,ptg,ptr,da_u,ptw,drst, &
784	uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
785	uhb_u,vhb_u,thb_u,ehb_u,ss,dt)
786
787
788	! Calculation of the sensible heat fluxes for the ground, the wall and roof
789	! Sensible Heat Flux = density * Cp_U * ( A* potential temperature + B )
790	! where A and B are the implicit and explicit components of the heat sources or sinks.
791	!
792	!
793
794	do id=1,nd_u(iurb)
795	sfg(id)=-da_u(1)cp_uthb_u(id,1)
796	do iz=2,nz_u(iurb)
797	sfr(id,iz)=-da_u(iz)cp_uthb_u(id,iz)
798	enddo
799
800	do iz=1,nz_u(iurb)
801	sfw(2id-1,iz)=-da_u(iz)cp_u(tvb_u(2id-1,iz)+ &
802	tva_u(2id-1,iz)pt_u(iz))
803	sfw(2id,iz)=-da_u(iz)cp_u(tvb_u(2id,iz)+ &
804	tva_u(2id,iz)pt_u(iz))
805	enddo
806	enddo
807
808	! calculation of the urban albedo and the upward long wave radiation
809
810	! call upward_rad(nd_u(iurb),iurb,nz_u(iurb),ws,bs,sigma,fsw,fsg,pb,ss, &
811	! tg,emg_u(iurb),albg_u(iurb),rlg,rsg, &
812	! tw,emw_u(iurb),albw_u(iurb),rlw,rsw, &
813	! tr,emr_u(iurb),albr_u(iurb),rld,rs, &
814	! rs_abs,rl_up,emiss)
815
816	! Interpolation on the "mesoscale grid"
817
818	call urban_meso(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z,dz,z_u,pb,ss,bs,ws,sf, &
819	vl,uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
820	uhb_u,vhb_u,thb_u,ehb_u, &
821	a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)
822
823
824	! computation of the mean road temperature tsk (this value could be used
825	! to replace the surface temperature in the radiation routines, if needed).
826
827	! tsk=0.
828	! do id=1,nd_u(iurb)
829	! tsk=tsk+tg(id,ng_u)/nd_u(iurb)
830	! enddo
831
832	! Calculation of the length scale taking into account the buildings effects
833
834	call interp_length(nd_u(iurb),kms,kme,kts,kte,nz_u(iurb),z_u,z,ss,ws,bs,dlg,dl_u)
835
836	return
837	end subroutine BEP1D
838
839	! ===6=8===============================================================72
840	! ===6=8===============================================================72
841
842	subroutine param(iurb,nz,nd, &
843	csg_u,csg,alag_u,alag,csr_u,csr, &
844	alar_u,alar,csw_u,csw,alaw_u,alaw, &
845	ws_u,ws,bs_u,bs,z0g_u,z0r_u,z0, &
846	strd_u,strd,drst_u,drst,ss_u,ss,pb_u,pb)
847
848	! ----------------------------------------------------------------------
849	! This routine prepare some usefull parameters
850	! ----------------------------------------------------------------------
851
852	implicit none
853
854
855	! ----------------------------------------------------------------------
856	! INPUT:
857	! ----------------------------------------------------------------------
858	integer iurb ! Current urban class
859	integer nz ! Number of vertical urban levels in the current class
860	integer nd ! Number of street direction for the current urban class
861	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
862	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
863	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
864	real bs_u(ndm,nurbm) ! Building width
865	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
866	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
867	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
868	real drst_u(ndm,nurbm) ! Street direction
869	real strd_u(ndm,nurbm) ! Street length
870	real ws_u(ndm,nurbm) ! Street width
871	real z0g_u(nurbm) ! The ground's roughness length
872	real z0r_u(nurbm) ! The roof's roughness length
873	real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to "z"
874	real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to "z"
875
876	! ----------------------------------------------------------------------
877	! OUTPUT:
878	! ----------------------------------------------------------------------
879	real alag(ng_u) ! Ground thermal diffusivity at each ground levels
880	real alar(nwr_u) ! Roof thermal diffusivity at each roof levels
881	real alaw(nwr_u) ! Wall thermal diffusivity at each wall levels
882	real csg(ng_u) ! Specific heat of the ground material at each ground levels
883	real csr(nwr_u) ! Specific heat of the roof material at each roof levels
884	real csw(nwr_u) ! Specific heat of the wall material at each wall levels
885	real bs(ndm) ! Building width for the current urban class
886	real drst(ndm) ! street directions for the current urban class
887	real strd(ndm) ! Street lengths for the current urban class
888	real ws(ndm) ! Street widths of the current urban class
889	real z0(ndm,nz_um) ! Roughness lengths "profiles"
890	real ss(nz_um) ! Probability to have a building with height h
891	real pb(nz_um) ! Probability to have a building with an height equal
892
893	! ----------------------------------------------------------------------
894	! LOCAL:
895	! ----------------------------------------------------------------------
896	integer id,ig,ir,iw,iz
897
898	! ----------------------------------------------------------------------
899	! END VARIABLES DEFINITIONS
900	! ----------------------------------------------------------------------
901	!
902	!Initialize the variables
903	!
904	ss=0.
905	pb=0.
906	csg=0.
907	alag=0.
908	csr=0.
909	alar=0.
910	csw=0.
911	alaw=0.
912	z0=0.
913	ws=0.
914	bs=0.
915	strd=0.
916	drst=0.
917
918	do iz=1,nz+1
919	ss(iz)=ss_u(iz,iurb)
920	pb(iz)=pb_u(iz,iurb)
921	end do
922
923	do ig=1,ng_u
924	csg(ig)=csg_u(iurb)
925	alag(ig)=alag_u(iurb)
926	enddo
927
928	do ir=1,nwr_u
929	csr(ir)=csr_u(iurb)
930	alar(ir)=alar_u(iurb)
931	enddo
932
933	do iw=1,nwr_u
934	csw(iw)=csw_u(iurb)
935	alaw(iw)=alaw_u(iurb)
936	enddo
937
938	do id=1,nd
939	z0(id,1)=z0g_u(iurb)
940	do iz=2,nz+1
941	z0(id,iz)=z0r_u(iurb)
942	enddo
943	enddo
944
945	do id=1,nd
946	ws(id)=ws_u(id,iurb)
947	bs(id)=bs_u(id,iurb)
948	strd(id)=strd_u(id,iurb)
949	drst(id)=drst_u(id,iurb)
950	enddo
951
952
953	return
954	end subroutine param
955
956	! ===6=8===============================================================72
957	! ===6=8===============================================================72
958
959	subroutine interpol(kms,kme,kts,kte,nz_u,z,z_u,c,c_u)
960
961	! ----------------------------------------------------------------------
962	! This routine interpolate para
963	! meters from the "mesoscale grid" to
964	! the "urban grid".
965	! See p300 Appendix B.1 of the BLM paper.
966	! ----------------------------------------------------------------------
967
968	implicit none
969
970	! ----------------------------------------------------------------------
971	! INPUT:
972	! ----------------------------------------------------------------------
973	! Data relative to the "mesoscale grid"
974	integer kts,kte,kms,kme
975	real z(kms:kme) ! Altitude of the cell interface
976	real c(kms:kme) ! Parameter which has to be interpolated
977	! Data relative to the "urban grid"
978	integer nz_u ! Number of levels
979	!! real z_u(nz_u+1) ! Altitude of the cell interface
980	real z_u(nz_um) ! Altitude of the cell interface
981	! ----------------------------------------------------------------------
982	! OUTPUT:
983	! ----------------------------------------------------------------------
984	!! real c_u(nz_u) ! Interpolated paramters in the "urban grid"
985	real c_u(nz_um) ! Interpolated paramters in the "urban grid"
986
987	! LOCAL:
988	! ----------------------------------------------------------------------
989	integer iz_u,iz
990	real ctot,dz
991
992	! ----------------------------------------------------------------------
993	! END VARIABLES DEFINITIONS
994	! ----------------------------------------------------------------------
995
996	do iz_u=1,nz_u
997	ctot=0.
998	do iz=kts,kte
999	dz=max(min(z(iz+1),z_u(iz_u+1))-max(z(iz),z_u(iz_u)),0.)
1000	ctot=ctot+c(iz)*dz
1001	enddo
1002	c_u(iz_u)=ctot/(z_u(iz_u+1)-z_u(iz_u))
1003	enddo
1004
1005	return
1006	end subroutine interpol
1007
1008	! ===6=8===============================================================72
1009	! ===6=8===============================================================72
1010
1011	subroutine modif_rad(iurb,nd,nz_u,z,ws,drst,strd,ss,pb, &
1012	tw,tg,albg,albw,emw,emg, &
1013	fww,fwg,fgw,fsw,fsg, &
1014	zr,deltar,ah, &
1015	rs,rl,rsw,rsg,rlw,rlg)
1016
1017	! ----------------------------------------------------------------------
1018	! This routine computes the modification of the short wave and
1019	! long wave radiation due to the buildings.
1020	! ----------------------------------------------------------------------
1021
1022	implicit none
1023
1024
1025	! ----------------------------------------------------------------------
1026	! INPUT:
1027	! ----------------------------------------------------------------------
1028	integer iurb ! current urban class
1029	integer nd ! Number of street direction for the current urban class
1030	integer nz_u ! Number of layer in the urban grid
1031	real z(nz_um) ! Height of the urban grid levels
1032	real ws(ndm) ! Street widths of the current urban class
1033	real drst(ndm) ! street directions for the current urban class
1034	real strd(ndm) ! Street lengths for the current urban class
1035	real ss(nz_um) ! probability to have a building with height h
1036	real pb(nz_um) ! probability to have a building with an height equal
1037	real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
1038	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1039	real albg ! Albedo of the ground for the current urban class
1040	real albw ! Albedo of the wall for the current urban class
1041	real emg ! Emissivity of ground for the current urban class
1042	real emw ! Emissivity of wall for the current urban class
1043	real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
1044	real fsg(ndm,nurbm) ! View factors from sky to ground
1045	real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
1046	real fws(nz_um,ndm,nurbm) ! View factors from wall to sky
1047	real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
1048	real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
1049	real ah ! Hour angle (it should come from the radiation routine)
1050	real zr ! zenith angle
1051	real deltar ! Declination of the sun
1052	real rs ! solar radiation
1053	real rl ! downward flux of the longwave radiation
1054	! ----------------------------------------------------------------------
1055	! OUTPUT:
1056	! ----------------------------------------------------------------------
1057	real rlg(ndm) ! Long wave radiation at the ground
1058	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1059	real rsg(ndm) ! Short wave radiation at the ground
1060	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1061
1062	! ----------------------------------------------------------------------
1063	! LOCAL:
1064	! ----------------------------------------------------------------------
1065
1066	integer id,iz
1067
1068	! Calculation of the shadow effects
1069
1070	call shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
1071	rs,rsw,rsg)
1072
1073	! Calculation of the reflection effects
1074	do id=1,nd
1075	call long_rad(iurb,nz_u,id,emw,emg, &
1076	fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
1077
1078	call short_rad(iurb,nz_u,id,albw,albg,fwg,fww,fgw,rsg,rsw,pb)
1079
1080	enddo
1081
1082	return
1083	end subroutine modif_rad
1084
1085
1086	! ===6=8===============================================================72
1087	! ===6=8===============================================================72
1088
1089	subroutine surf_temp(nz_u,nd,pr,dt,ss,rs,rl,rsg,rlg,rsw,rlw, &
1090	tg,alag,csg,emg,albg,ptg,sfg,gfg, &
1091	tr,alar,csr,emr,albr,ptr,sfr,gfr, &
1092	tw,alaw,csw,emw,albw,ptw,sfw,gfw)
1093
1094
1095	! ----------------------------------------------------------------------
1096	! Computation of the surface temperatures for walls, ground and roofs
1097	! ----------------------------------------------------------------------
1098
1099	implicit none
1100
1101
1102
1103	! ----------------------------------------------------------------------
1104	! INPUT:
1105	! ----------------------------------------------------------------------
1106	integer nz_u ! Number of vertical layers defined in the urban grid
1107	integer nd ! Number of street direction for the current urban class
1108	real alag(ng_u) ! Ground thermal diffusivity for the current urban class [m^2 s^-1]
1109	real alar(nwr_u) ! Roof thermal diffusivity for the current urban class [m^2 s^-1]
1110	real alaw(nwr_u) ! Wall thermal diffusivity for the current urban class [m^2 s^-1]
1111	real albg ! Albedo of the ground for the current urban class
1112	real albr ! Albedo of the roof for the current urban class
1113	real albw ! Albedo of the wall for the current urban class
1114	real csg(ng_u) ! Specific heat of the ground material of the current urban class [J m^3 K^-1]
1115	real csr(nwr_u) ! Specific heat of the roof material for the current urban class [J m^3 K^-1]
1116	real csw(nwr_u) ! Specific heat of the wall material for the current urban class [J m^3 K^-1]
1117	real dt ! Time step
1118	real emg ! Emissivity of ground for the current urban class
1119	real emr ! Emissivity of roof for the current urban class
1120	real emw ! Emissivity of wall for the current urban class
1121	real pr(nz_um) ! Air pressure
1122	real rs ! Solar radiation
1123	real rl ! Downward flux of the longwave radiation
1124	real rlg(ndm) ! Long wave radiation at the ground
1125	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1126	real rsg(ndm) ! Short wave radiation at the ground
1127	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1128	real sfg(ndm) ! Sensible heat flux from ground (road)
1129	real sfr(ndm,nz_um) ! Sensible heat flux from roofs
1130	real sfw(2*ndm,nz_um) ! Sensible heat flux from walls
1131	real gfg(ndm) ! Heat flux transferred from the surface of the ground (road) toward the interior
1132	real gfr(ndm,nz_um) ! Heat flux transferred from the surface of the roof toward the interior
1133	real gfw(2*ndm,nz_um) ! Heat flux transfered from the surface of the walls toward the interior
1134	real ss(nz_um) ! Probability to have a building with height h
1135	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1136	real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
1137	real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
1138
1139
1140	! ----------------------------------------------------------------------
1141	! OUTPUT:
1142	! ----------------------------------------------------------------------
1143	real ptg(ndm) ! Ground potential temperatures
1144	real ptr(ndm,nz_um) ! Roof potential temperatures
1145	real ptw(2*ndm,nz_um) ! Walls potential temperatures
1146
1147	! ----------------------------------------------------------------------
1148	! LOCAL:
1149	! ----------------------------------------------------------------------
1150	integer id,ig,ir,iw,iz
1151
1152	real rtg(ndm) ! Total radiation at ground(road) surface (solar+incoming long+outgoing long)
1153	real rtr(ndm,nz_um) ! Total radiation at roof surface (solar+incoming long+outgoing long)
1154	real rtw(2*ndm,nz_um) ! Radiation at walls surface (solar+incoming long+outgoing long)
1155	real tg_tmp(ng_u)
1156	real tr_tmp(nwr_u)
1157	real tw_tmp(nwr_u)
1158
1159	real dzg_u(ng_u) ! Layer sizes in the ground
1160
1161	real dzr_u(nwr_u) ! Layers sizes in the roof
1162
1163	real dzw_u(nwr_u) ! Layer sizes in the wall
1164
1165
1166	data dzg_u /0.2,0.12,0.08,0.05,0.03,0.02,0.02,0.01,0.005,0.0025/
1167	data dzr_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
1168	data dzw_u /0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.005,0.0025/
1169	! ----------------------------------------------------------------------
1170	! END VARIABLES DEFINITIONS
1171	! ----------------------------------------------------------------------
1172
1173
1174
1175	do id=1,nd
1176
1177	! Calculation for the ground surfaces
1178	do ig=1,ng_u
1179	tg_tmp(ig)=tg(id,ig)
1180	end do
1181	!
1182	call soil_temp(ng_u,dzg_u,tg_tmp,ptg(id),alag,csg, &
1183	rsg(id),rlg(id),pr(1), &
1184	dt,emg,albg, &
1185	rtg(id),sfg(id),gfg(id))
1186	do ig=1,ng_u
1187	tg(id,ig)=tg_tmp(ig)
1188	end do
1189
1190	! Calculation for the roofs surfaces
1191
1192	do iz=2,nz_u
1193
1194	if(ss(iz).gt.0.)then
1195	do ir=1,nwr_u
1196	tr_tmp(ir)=tr(id,iz,ir)
1197	end do
1198
1199	call soil_temp(nwr_u,dzr_u,tr_tmp,ptr(id,iz), &
1200	alar,csr,rs,rl,pr(iz),dt,emr,albr, &
1201	rtr(id,iz),sfr(id,iz),gfr(id,iz))
1202	do ir=1,nwr_u
1203	tr(id,iz,ir)=tr_tmp(ir)
1204	end do
1205
1206	end if
1207
1208	end do !iz
1209
1210	! Calculation for the walls surfaces
1211
1212	do iz=1,nz_u
1213
1214	do iw=1,nwr_u
1215	tw_tmp(iw)=tw(2*id-1,iz,iw)
1216	end do
1217	call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id-1,iz),alaw, &
1218	csw, &
1219	rsw(2id-1,iz),rlw(2id-1,iz), &
1220	pr(iz),dt,emw, &
1221	albw,rtw(2id-1,iz),sfw(2id-1,iz),gfw(2*id-1,iz))
1222
1223	do iw=1,nwr_u
1224	tw(2*id-1,iz,iw)=tw_tmp(iw)
1225	end do
1226
1227	do iw=1,nwr_u
1228	tw_tmp(iw)=tw(2*id,iz,iw)
1229	end do
1230
1231	call soil_temp(nwr_u,dzw_u,tw_tmp,ptw(2*id,iz),alaw, &
1232	csw, &
1233	rsw(2id,iz),rlw(2id,iz), &
1234	pr(iz),dt,emw, &
1235	albw,rtw(2id,iz),sfw(2id,iz),gfw(2*id,iz))
1236	do iw=1,nwr_u
1237	tw(2*id,iz,iw)=tw_tmp(iw)
1238	end do
1239
1240	end do !iz
1241
1242	end do !id
1243
1244	return
1245	end subroutine surf_temp
1246
1247	! ===6=8===============================================================72
1248	! ===6=8===============================================================72
1249
1250	subroutine buildings(nd,nz,z0,ua_u,va_u,pt_u,pt0_u, &
1251	ptg,ptr,da_u,ptw, &
1252	drst,uva_u,vva_u,uvb_u,vvb_u, &
1253	tva_u,tvb_u,evb_u, &
1254	uhb_u,vhb_u,thb_u,ehb_u,ss,dt)
1255
1256	! ----------------------------------------------------------------------
1257	! This routine computes the sources or sinks of the different quantities
1258	! on the urban grid. The actual calculation is done in the subroutines
1259	! called flux_wall, and flux_flat.
1260	! ----------------------------------------------------------------------
1261
1262	implicit none
1263
1264
1265	! ----------------------------------------------------------------------
1266	! INPUT:
1267	! ----------------------------------------------------------------------
1268	integer nd ! Number of street direction for the current urban class
1269	integer nz ! number of vertical space steps
1270	real ua_u(nz_um) ! Wind speed in the x direction on the urban grid
1271	real va_u(nz_um) ! Wind speed in the y direction on the urban grid
1272	real da_u(nz_um) ! air density on the urban grid
1273	real drst(ndm) ! Street directions for the current urban class
1274	real dz
1275	real pt_u(nz_um) ! Potential temperature on the urban grid
1276	real pt0_u(nz_um) ! reference potential temperature on the urban grid
1277	real ptg(ndm) ! Ground potential temperatures
1278	real ptr(ndm,nz_um) ! Roof potential temperatures
1279	real ptw(2*ndm,nz_um) ! Walls potential temperatures
1280	real ss(nz_um) ! probability to have a building with height h
1281	real z0(ndm,nz_um) ! Roughness lengths "profiles"
1282	real dt ! time step
1283
1284
1285	! ----------------------------------------------------------------------
1286	! OUTPUT:
1287	! ----------------------------------------------------------------------
1288	! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
1289	! vertical surfaces (walls) and horizontal surfaces (roofs and street)
1290	! The fluxes can be computed as follow: Fluxes of X = A*X + B
1291	! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
1292
1293	real uhb_u(ndm,nz_um) ! U (wind component) Horizontal surfaces, B (explicit) term
1294	real uva_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, A (implicit) term
1295	real uvb_u(2*ndm,nz_um) ! U (wind component) Vertical surfaces, B (explicit) term
1296	real vhb_u(ndm,nz_um) ! V (wind component) Horizontal surfaces, B (explicit) term
1297	real vva_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, A (implicit) term
1298	real vvb_u(2*ndm,nz_um) ! V (wind component) Vertical surfaces, B (explicit) term
1299	real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1300	real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1301	real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1302	real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1303	real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1304
1305	! ----------------------------------------------------------------------
1306	! LOCAL:
1307	! ----------------------------------------------------------------------
1308	integer id,iz
1309
1310	! ----------------------------------------------------------------------
1311	! END VARIABLES DEFINITIONS
1312	! ----------------------------------------------------------------------
1313	dz=dz_u
1314
1315	do id=1,nd
1316
1317	! Calculation at the ground surfaces
1318	call flux_flat(dz,z0(id,1),ua_u(1),va_u(1),pt_u(1),pt0_u(1), &
1319	ptg(id),uhb_u(id,1), &
1320	vhb_u(id,1),thb_u(id,1),ehb_u(id,1))
1321
1322	! Calculation at the roof surfaces
1323	do iz=2,nz
1324	if(ss(iz).gt.0)then
1325	call flux_flat(dz,z0(id,iz),ua_u(iz), &
1326	va_u(iz),pt_u(iz),pt0_u(iz), &
1327	ptr(id,iz),uhb_u(id,iz), &
1328	vhb_u(id,iz),thb_u(id,iz),ehb_u(id,iz))
1329	else
1330	uhb_u(id,iz) = 0.0
1331	vhb_u(id,iz) = 0.0
1332	thb_u(id,iz) = 0.0
1333	ehb_u(id,iz) = 0.0
1334	endif
1335	end do
1336
1337	! Calculation at the wall surfaces
1338	do iz=1,nz
1339	call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
1340	ptw(2*id-1,iz), &
1341	uva_u(2id-1,iz),vva_u(2id-1,iz), &
1342	uvb_u(2id-1,iz),vvb_u(2id-1,iz), &
1343	tva_u(2id-1,iz),tvb_u(2id-1,iz), &
1344	evb_u(2*id-1,iz),drst(id),dt)
1345
1346	call flux_wall(ua_u(iz),va_u(iz),pt_u(iz),da_u(iz), &
1347	ptw(2*id,iz), &
1348	uva_u(2id,iz),vva_u(2id,iz), &
1349	uvb_u(2id,iz),vvb_u(2id,iz), &
1350	tva_u(2id,iz),tvb_u(2id,iz), &
1351	evb_u(2*id,iz),drst(id),dt)
1352	!
1353
1354	end do
1355
1356	end do
1357
1358	return
1359	end subroutine buildings
1360
1361
1362	! ===6=8===============================================================72
1363	! ===6=8===============================================================72
1364
1365	subroutine urban_meso(nd,kms,kme,kts,kte,nz_u,z,dz,z_u,pb,ss,bs,ws,sf,vl, &
1366	uva_u,vva_u,uvb_u,vvb_u,tva_u,tvb_u,evb_u, &
1367	uhb_u,vhb_u,thb_u,ehb_u, &
1368	a_u,a_v,a_t,a_e,b_u,b_v,b_t,b_e)
1369
1370	! ----------------------------------------------------------------------
1371	! This routine interpolates the parameters from the "urban grid" to the
1372	! "mesoscale grid".
1373	! See p300-301 Appendix B.2 of the BLM paper.
1374	! ----------------------------------------------------------------------
1375
1376	implicit none
1377
1378	! ----------------------------------------------------------------------
1379	! INPUT:
1380	! ----------------------------------------------------------------------
1381	! Data relative to the "mesoscale grid"
1382	integer kms,kme,kts,kte
1383	real z(kms:kme) ! Altitude above the ground of the cell interface
1384	real dz(kms:kme) ! Vertical space steps
1385
1386	! Data relative to the "uban grid"
1387	integer nz_u ! Number of layer in the urban grid
1388	integer nd ! Number of street direction for the current urban class
1389	real bs(ndm) ! Building widths of the current urban class
1390	real ws(ndm) ! Street widths of the current urban class
1391	real z_u(nz_um) ! Height of the urban grid levels
1392	real pb(nz_um) ! Probability to have a building with an height equal
1393	real ss(nz_um) ! Probability to have a building with height h
1394	real uhb_u(ndm,nz_um) ! U (x-wind component) Horizontal surfaces, B (explicit) term
1395	real uva_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, A (implicit) term
1396	real uvb_u(2*ndm,nz_um) ! U (x-wind component) Vertical surfaces, B (explicit) term
1397	real vhb_u(ndm,nz_um) ! V (y-wind component) Horizontal surfaces, B (explicit) term
1398	real vva_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, A (implicit) term
1399	real vvb_u(2*ndm,nz_um) ! V (y-wind component) Vertical surfaces, B (explicit) term
1400	real thb_u(ndm,nz_um) ! Temperature Horizontal surfaces, B (explicit) term
1401	real tva_u(2*ndm,nz_um) ! Temperature Vertical surfaces, A (implicit) term
1402	real tvb_u(2*ndm,nz_um) ! Temperature Vertical surfaces, B (explicit) term
1403	real ehb_u(ndm,nz_um) ! Energy (TKE) Horizontal surfaces, B (explicit) term
1404	real evb_u(2*ndm,nz_um) ! Energy (TKE) Vertical surfaces, B (explicit) term
1405
1406	! ----------------------------------------------------------------------
1407	! OUTPUT:
1408	! ----------------------------------------------------------------------
1409	! Data relative to the "mesoscale grid"
1410	real sf(kms:kme) ! Surface of the "mesoscale grid" cells taking into account the buildings
1411	real vl(kms:kme) ! Volume of the "mesoscale grid" cells taking into account the buildings
1412	real a_u(kms:kme) ! Implicit component of the momentum sources or sinks in the X-direction
1413	real a_v(kms:kme) ! Implicit component of the momentum sources or sinks in the Y-direction
1414	real a_t(kms:kme) ! Implicit component of the heat sources or sinks
1415	real a_e(kms:kme) ! Implicit component of the TKE sources or sinks
1416	real b_u(kms:kme) ! Explicit component of the momentum sources or sinks in the X-direction
1417	real b_v(kms:kme) ! Explicit component of the momentum sources or sinks in the Y-direction
1418	real b_t(kms:kme) ! Explicit component of the heat sources or sinks
1419	real b_e(kms:kme) ! Explicit component of the TKE sources or sinks
1420
1421	! ----------------------------------------------------------------------
1422	! LOCAL:
1423	! ----------------------------------------------------------------------
1424	real dzz
1425	real fact
1426	integer id,iz,iz_u
1427	real se,sr,st,su,sv
1428	real uet(kms:kme) ! Contribution to TKE due to walls
1429	real veb,vta,vtb,vte,vtot,vua,vub,vva,vvb
1430
1431
1432	! ----------------------------------------------------------------------
1433	! END VARIABLES DEFINITIONS
1434	! ----------------------------------------------------------------------
1435
1436	! initialisation
1437
1438	do iz=kts,kte
1439	a_u(iz)=0.
1440	a_v(iz)=0.
1441	a_t(iz)=0.
1442	a_e(iz)=0.
1443	b_u(iz)=0.
1444	b_v(iz)=0.
1445	b_e(iz)=0.
1446	b_t(iz)=0.
1447	uet(iz)=0.
1448	end do
1449
1450	! horizontal surfaces
1451	do iz=kts,kte
1452	sf(iz)=0.
1453	vl(iz)=0.
1454	enddo
1455	sf(kte+1)=0.
1456
1457	do id=1,nd
1458	do iz=kts+1,kte+1
1459	sr=0.
1460	do iz_u=2,nz_u
1461	if(z(iz).lt.z_u(iz_u).and.z(iz).ge.z_u(iz_u-1))then
1462	sr=pb(iz_u)
1463	endif
1464	enddo
1465	sf(iz)=sf(iz)+((ws(id)+(1.-sr)*bs(id))/(ws(id)+bs(id)))/nd
1466	enddo
1467	enddo
1468
1469	! volume
1470	do iz=kts,kte
1471	do id=1,nd
1472	vtot=0.
1473	do iz_u=1,nz_u
1474	dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
1475	vtot=vtot+pb(iz_u+1)*dzz
1476	enddo
1477	vtot=vtot/(z(iz+1)-z(iz))
1478	vl(iz)=vl(iz)+(1.-vtot*bs(id)/(ws(id)+bs(id)))/nd
1479	enddo
1480	enddo
1481
1482	! horizontal surface impact
1483
1484	do id=1,nd
1485
1486	fact=1./vl(kts)/dz(kts)*ws(id)/(ws(id)+bs(id))/nd
1487	b_t(kts)=b_t(kts)+thb_u(id,1)*fact
1488	b_u(kts)=b_u(kts)+uhb_u(id,1)*fact
1489	b_v(kts)=b_v(kts)+vhb_u(id,1)*fact
1490	b_e(kts)=b_e(kts)+ehb_u(id,1)fact(z_u(2)-z_u(1))
1491
1492	do iz=kts,kte
1493	st=0.
1494	su=0.
1495	sv=0.
1496	se=0.
1497	do iz_u=2,nz_u
1498	if(z(iz).le.z_u(iz_u).and.z(iz+1).gt.z_u(iz_u))then
1499	st=st+ss(iz_u)*thb_u(id,iz_u)
1500	su=su+ss(iz_u)*uhb_u(id,iz_u)
1501	sv=sv+ss(iz_u)*vhb_u(id,iz_u)
1502	se=se+ss(iz_u)ehb_u(id,iz_u)(z_u(iz_u+1)-z_u(iz_u))
1503	endif
1504	enddo
1505
1506	fact=bs(id)/(ws(id)+bs(id))/vl(iz)/dz(iz)/nd
1507	b_t(iz)=b_t(iz)+st*fact
1508	b_u(iz)=b_u(iz)+su*fact
1509	b_v(iz)=b_v(iz)+sv*fact
1510	b_e(iz)=b_e(iz)+se*fact
1511	enddo
1512	enddo
1513
1514	! vertical surface impact
1515
1516	do iz=kts,kte
1517	uet(iz)=0.
1518	do id=1,nd
1519	vtb=0.
1520	vta=0.
1521	vua=0.
1522	vub=0.
1523	vva=0.
1524	vvb=0.
1525	veb=0.
1526	vte=0.
1527	do iz_u=1,nz_u
1528	dzz=max(min(z_u(iz_u+1),z(iz+1))-max(z_u(iz_u),z(iz)),0.)
1529	fact=dzz/(ws(id)+bs(id))
1530	vtb=vtb+pb(iz_u+1)* &
1531	(tvb_u(2id-1,iz_u)+tvb_u(2id,iz_u))*fact
1532	vta=vta+pb(iz_u+1)* &
1533	(tva_u(2id-1,iz_u)+tva_u(2id,iz_u))*fact
1534	vua=vua+pb(iz_u+1)* &
1535	(uva_u(2id-1,iz_u)+uva_u(2id,iz_u))*fact
1536	vva=vva+pb(iz_u+1)* &
1537	(vva_u(2id-1,iz_u)+vva_u(2id,iz_u))*fact
1538	vub=vub+pb(iz_u+1)* &
1539	(uvb_u(2id-1,iz_u)+uvb_u(2id,iz_u))*fact
1540	vvb=vvb+pb(iz_u+1)* &
1541	(vvb_u(2id-1,iz_u)+vvb_u(2id,iz_u))*fact
1542	veb=veb+pb(iz_u+1)* &
1543	(evb_u(2id-1,iz_u)+evb_u(2id,iz_u))*fact
1544	enddo
1545
1546	fact=1./vl(iz)/dz(iz)/nd
1547	b_t(iz)=b_t(iz)+vtb*fact
1548	a_t(iz)=a_t(iz)+vta*fact
1549	a_u(iz)=a_u(iz)+vua*fact
1550	a_v(iz)=a_v(iz)+vva*fact
1551	b_u(iz)=b_u(iz)+vub*fact
1552	b_v(iz)=b_v(iz)+vvb*fact
1553	b_e(iz)=b_e(iz)+veb*fact
1554	uet(iz)=uet(iz)+vte*fact
1555	enddo
1556	enddo
1557
1558
1559
1560	return
1561	end subroutine urban_meso
1562
1563	! ===6=8===============================================================72
1564
1565	! ===6=8===============================================================72
1566
1567	subroutine interp_length(nd,kms,kme,kts,kte,nz_u,z_u,z,ss,ws,bs, &
1568	dlg,dl_u)
1569
1570	! ----------------------------------------------------------------------
1571	! Calculation of the length scales
1572	! See p272-274 formula (22) and (24) of the BLM paper
1573	! ----------------------------------------------------------------------
1574
1575	implicit none
1576
1577
1578	! ----------------------------------------------------------------------
1579	! INPUT:
1580	! ----------------------------------------------------------------------
1581	integer kms,kme,kts,kte
1582	real z(kms:kme) ! Altitude above the ground of the cell interface
1583	integer nd ! Number of street direction for the current urban class
1584	integer nz_u ! Number of levels in the "urban grid"
1585	real z_u(nz_um) ! Height of the urban grid levels
1586	real bs(ndm) ! Building widths of the current urban class
1587	real ss(nz_um) ! Probability to have a building with height h
1588	real ws(ndm) ! Street widths of the current urban class
1589
1590
1591	! ----------------------------------------------------------------------
1592	! OUTPUT:
1593	! ----------------------------------------------------------------------
1594	real dlg(kms:kme) ! Height above ground (L_ground in formula (24) of the BLM paper).
1595	real dl_u(kms:kme) ! Length scale (lb in formula (22) ofthe BLM paper).
1596
1597	! ----------------------------------------------------------------------
1598	! LOCAL:
1599	! ----------------------------------------------------------------------
1600	real dlgtmp
1601	integer id,iz,iz_u
1602	real sftot
1603	real ulu,ssl
1604
1605	! ----------------------------------------------------------------------
1606	! END VARIABLES DEFINITIONS
1607	! ----------------------------------------------------------------------
1608
1609	do iz=kts,kte
1610	ulu=0.
1611	ssl=0.
1612	do id=1,nd
1613	do iz_u=2,nz_u
1614	if(z_u(iz_u).gt.z(iz))then
1615	ulu=ulu+ss(iz_u)/z_u(iz_u)/nd
1616	ssl=ssl+ss(iz_u)/nd
1617	endif
1618	enddo
1619	enddo
1620
1621	if(ulu.ne.0)then
1622	dl_u(iz)=ssl/ulu
1623	else
1624	dl_u(iz)=0.
1625	endif
1626	enddo
1627
1628
1629	do iz=kts,kte
1630	dlg(iz)=0.
1631	do id=1,nd
1632	sftot=ws(id)
1633	dlgtmp=ws(id)/((z(iz)+z(iz+1))/2.)
1634	do iz_u=1,nz_u
1635	if((z(iz)+z(iz+1))/2..gt.z_u(iz_u))then
1636	dlgtmp=dlgtmp+ss(iz_u)*bs(id)/ &
1637	((z(iz)+z(iz+1))/2.-z_u(iz_u))
1638	sftot=sftot+ss(iz_u)*bs(id)
1639	endif
1640	enddo
1641	dlg(iz)=dlg(iz)+dlgtmp/sftot/nd
1642	enddo
1643	dlg(iz)=1./dlg(iz)
1644	enddo
1645
1646	return
1647	end subroutine interp_length
1648
1649	! ===6=8===============================================================72
1650	! ===6=8===============================================================72
1651
1652	subroutine shadow_mas(nd,nz_u,zr,deltar,ah,drst,ws,ss,pb,z, &
1653	rs,rsw,rsg)
1654
1655	! ----------------------------------------------------------------------
1656	! Modification of short wave radiation to take into account
1657	! the shadow produced by the buildings
1658	! ----------------------------------------------------------------------
1659
1660	implicit none
1661
1662	! ----------------------------------------------------------------------
1663	! INPUT:
1664	! ----------------------------------------------------------------------
1665	integer nd ! Number of street direction for the current urban class
1666	integer nz_u ! number of vertical layers defined in the urban grid
1667	real ah ! Hour angle (it should come from the radiation routine)
1668	real deltar ! Declination of the sun
1669	real drst(ndm) ! street directions for the current urban class
1670	real rs ! solar radiation
1671	real ss(nz_um) ! probability to have a building with height h
1672	real pb(nz_um) ! Probability that a building has an height greater or equal to h
1673	real ws(ndm) ! Street width of the current urban class
1674	real z(nz_um) ! Height of the urban grid levels
1675	real zr ! zenith angle
1676
1677	! ----------------------------------------------------------------------
1678	! OUTPUT:
1679	! ----------------------------------------------------------------------
1680	real rsg(ndm) ! Short wave radiation at the ground
1681	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
1682
1683	! ----------------------------------------------------------------------
1684	! LOCAL:
1685	! ----------------------------------------------------------------------
1686	integer id,iz,jz
1687	real aae,aaw,bbb,phix,rd,rtot,wsd
1688
1689	! ----------------------------------------------------------------------
1690	! END VARIABLES DEFINITIONS
1691	! ----------------------------------------------------------------------
1692
1693	if(rs.eq.0.or.sin(zr).eq.1)then
1694	do id=1,nd
1695	rsg(id)=0.
1696	do iz=1,nz_u
1697	rsw(2*id-1,iz)=0.
1698	rsw(2*id,iz)=0.
1699	enddo
1700	enddo
1701	else
1702	!test
1703	if(abs(sin(zr)).gt.1.e-10)then
1704	if(cos(deltar)*sin(ah)/sin(zr).ge.1)then
1705	bbb=pi/2.
1706	elseif(cos(deltar)*sin(ah)/sin(zr).le.-1)then
1707	bbb=-pi/2.
1708	else
1709	bbb=asin(cos(deltar)*sin(ah)/sin(zr))
1710	endif
1711	else
1712	if(cos(deltar)*sin(ah).ge.0)then
1713	bbb=pi/2.
1714	elseif(cos(deltar)*sin(ah).lt.0)then
1715	bbb=-pi/2.
1716	endif
1717	endif
1718
1719	phix=zr
1720
1721	do id=1,nd
1722
1723	rsg(id)=0.
1724
1725	aae=bbb-drst(id)
1726	aaw=bbb-drst(id)+pi
1727
1728	do iz=1,nz_u
1729	rsw(2*id-1,iz)=0.
1730	rsw(2*id,iz)=0.
1731	if(pb(iz+1).gt.0.)then
1732	do jz=1,nz_u
1733	if(abs(sin(aae)).gt.1.e-10)then
1734	call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aae, &
1735	ws(id),rd)
1736	rsw(2id-1,iz)=rsw(2id-1,iz)+rsrdss(jz+1)/pb(iz+1)
1737
1738	endif
1739
1740	if(abs(sin(aaw)).gt.1.e-10)then
1741	call shade_wall(z(iz),z(iz+1),z(jz+1),phix,aaw, &
1742	ws(id),rd)
1743	rsw(2id,iz)=rsw(2id,iz)+rsrdss(jz+1)/pb(iz+1)
1744
1745	endif
1746	enddo
1747	endif
1748	enddo
1749	if(abs(sin(aae)).gt.1.e-10)then
1750	wsd=abs(ws(id)/sin(aae))
1751
1752	do jz=1,nz_u
1753	rd=max(0.,wsd-z(jz+1)*tan(phix))
1754	rsg(id)=rsg(id)+rsrdss(jz+1)/wsd
1755	enddo
1756	rtot=0.
1757
1758	do iz=1,nz_u
1759	rtot=rtot+(rsw(2id,iz)+rsw(2id-1,iz))* &
1760	(z(iz+1)-z(iz))
1761	enddo
1762	rtot=rtot+rsg(id)*ws(id)
1763	else
1764	rsg(id)=rs
1765	endif
1766
1767	enddo
1768	endif
1769
1770	return
1771	end subroutine shadow_mas
1772
1773	! ===6=8===============================================================72
1774	! ===6=8===============================================================72
1775
1776	subroutine shade_wall(z1,z2,hu,phix,aa,ws,rd)
1777
1778	! ----------------------------------------------------------------------
1779	! This routine computes the effects of a shadow induced by a building of
1780	! height hu, on a portion of wall between z1 and z2. See equation A10,
1781	! and correction described below formula A11, and figure A1. Basically rd
1782	! is the ratio between the horizontal surface illuminated and the portion
1783	! of wall. Referring to figure A1, multiplying radiation flux density on
1784	! a horizontal surface (rs) by x1-x2 we have the radiation energy per
1785	! unit time. Dividing this by z2-z1, we obtain the radiation flux
1786	! density reaching the portion of the wall between z2 and z1
1787	! (everything is assumed in 2D)
1788	! ----------------------------------------------------------------------
1789
1790	implicit none
1791
1792	! ----------------------------------------------------------------------
1793	! INPUT:
1794	! ----------------------------------------------------------------------
1795	real aa ! Angle between the sun direction and the face of the wall (A12)
1796	real hu ! Height of the building that generates the shadow
1797	real phix ! Solar zenith angle
1798	real ws ! Width of the street
1799	real z1 ! Height of the level z(iz)
1800	real z2 ! Height of the level z(iz+1)
1801
1802	! ----------------------------------------------------------------------
1803	! OUTPUT:
1804	! ----------------------------------------------------------------------
1805	real rd ! Ratio between (x1-x2)/(z2-z1), see Fig. 1A.
1806	! Multiplying rd by rs (radiation flux
1807	! density on a horizontal surface) gives
1808	! the radiation flux density on the
1809	! portion of wall between z1 and z2.
1810	! ----------------------------------------------------------------------
1811	! LOCAL:
1812	! ----------------------------------------------------------------------
1813	real x1,x2 ! x1,x2 see Fig. A1.
1814
1815	! ----------------------------------------------------------------------
1816	! END VARIABLES DEFINITIONS
1817	! ----------------------------------------------------------------------
1818
1819	x1=min((hu-z1)*tan(phix),max(0.,ws/sin(aa)))
1820
1821	x2=max((hu-z2)*tan(phix),0.)
1822
1823	rd=max(0.,sin(aa)*(max(0.,x1-x2))/(z2-z1))
1824
1825	return
1826	end subroutine shade_wall
1827
1828	! ===6=8===============================================================72
1829	! ===6=8===============================================================72
1830
1831	subroutine long_rad(iurb,nz_u,id,emw,emg, &
1832	fwg,fww,fgw,fsw,fsg,tg,tw,rlg,rlw,rl,pb)
1833
1834	! ----------------------------------------------------------------------
1835	! This routine computes the effects of the reflections of long-wave
1836	! radiation in the street canyon by solving the system
1837	! of 2nz_u+1 eqn. in 2nz_u+1
1838	! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
1839	! The system is solved by solving A X= B,
1840	! with A matrix B vector, and X solution.
1841	! ----------------------------------------------------------------------
1842
1843	implicit none
1844
1845
1846
1847	! ----------------------------------------------------------------------
1848	! INPUT:
1849	! ----------------------------------------------------------------------
1850	real emg ! Emissivity of ground for the current urban class
1851	real emw ! Emissivity of wall for the current urban class
1852	real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
1853	real fsg(ndm,nurbm) ! View factors from sky to ground
1854	real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
1855	real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
1856	real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
1857	integer id ! Current street direction
1858	integer iurb ! Current urban class
1859	integer nz_u ! Number of layer in the urban grid
1860	real pb(nz_um) ! Probability to have a building with an height equal
1861	real rl ! Downward flux of the longwave radiation
1862	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
1863	real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
1864
1865
1866	! ----------------------------------------------------------------------
1867	! OUTPUT:
1868	! ----------------------------------------------------------------------
1869	real rlg(ndm) ! Long wave radiation at the ground
1870	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls
1871
1872	! ----------------------------------------------------------------------
1873	! LOCAL:
1874	! ----------------------------------------------------------------------
1875	integer i,j
1876	real aaa(2nz_um+1,2nz_um+1) ! terms of the matrix
1877	real bbb(2*nz_um+1) ! terms of the vector
1878
1879	! ----------------------------------------------------------------------
1880	! END VARIABLES DEFINITIONS
1881	! ----------------------------------------------------------------------
1882
1883
1884	! west wall
1885
1886	do i=1,nz_u
1887
1888	do j=1,nz_u
1889	aaa(i,j)=0.
1890	enddo
1891
1892	aaa(i,i)=1.
1893
1894	do j=nz_u+1,2*nz_u
1895	aaa(i,j)=-(1.-emw)fww(j-nz_u,i,id,iurb)pb(j-nz_u+1)
1896	enddo
1897
1898	!! aaa(i,2nz_u+1)=-(1.-emg)fgw(i,id,iurb)*pb(i+1)
1899	aaa(i,2nz_u+1)=-(1.-emg)fgw(i,id,iurb)
1900
1901	bbb(i)=fsw(i,id,iurb)rl+emgfgw(i,id,iurb)sigmatg(id,ng_u)**4
1902	do j=1,nz_u
1903	bbb(i)=bbb(i)+pb(j+1)emwsigmafww(j,i,id,iurb) &
1904	tw(2id,j,nwr_u)*4+ &
1905	fww(j,i,id,iurb)rl(1.-pb(j+1))
1906	enddo
1907
1908	enddo
1909
1910	! east wall
1911
1912	do i=1+nz_u,2*nz_u
1913
1914	do j=1,nz_u
1915	aaa(i,j)=-(1.-emw)fww(j,i-nz_u,id,iurb)pb(j+1)
1916	enddo
1917
1918	do j=1+nz_u,2*nz_u
1919	aaa(i,j)=0.
1920	enddo
1921
1922	aaa(i,i)=1.
1923
1924	!! aaa(i,2nz_u+1)=-(1.-emg)fgw(i-nz_u,id,iurb)*pb(i-nz_u+1)
1925	aaa(i,2nz_u+1)=-(1.-emg)fgw(i-nz_u,id,iurb)
1926
1927	bbb(i)=fsw(i-nz_u,id,iurb)*rl+ &
1928	emgfgw(i-nz_u,id,iurb)sigmatg(id,ng_u)*4
1929
1930	do j=1,nz_u
1931	bbb(i)=bbb(i)+pb(j+1)emwsigmafww(j,i-nz_u,id,iurb) &
1932	tw(2id-1,j,nwr_u)*4+ &
1933	fww(j,i-nz_u,id,iurb)rl(1.-pb(j+1))
1934	enddo
1935
1936	enddo
1937
1938	! ground
1939	do j=1,nz_u
1940	aaa(2nz_u+1,j)=-(1.-emw)fwg(j,id,iurb)*pb(j+1)
1941	enddo
1942
1943	do j=nz_u+1,2*nz_u
1944	aaa(2nz_u+1,j)=-(1.-emw)fwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
1945	enddo
1946
1947	aaa(2nz_u+1,2nz_u+1)=1.
1948
1949	bbb(2nz_u+1)=fsg(id,iurb)rl
1950
1951	do i=1,nz_u
1952	bbb(2nz_u+1)=bbb(2nz_u+1)+emwsigmafwg(i,id,iurb)pb(i+1) &
1953	(tw(2id-1,i,nwr_u)4+tw(2id,i,nwr_u)**4)+ &
1954	2.fwg(i,id,iurb)(1.-pb(i+1))*rl
1955	enddo
1956
1957
1958
1959	call gaussj(aaa,2nz_u+1,bbb,2nz_um+1)
1960
1961	do i=1,nz_u
1962	rlw(2*id-1,i)=bbb(i)
1963	enddo
1964
1965	do i=nz_u+1,2*nz_u
1966	rlw(2*id,i-nz_u)=bbb(i)
1967	enddo
1968
1969	rlg(id)=bbb(2*nz_u+1)
1970
1971	return
1972	end subroutine long_rad
1973
1974	! ===6=8===============================================================72
1975	! ===6=8===============================================================72
1976
1977	subroutine short_rad(iurb,nz_u,id,albw, &
1978	albg,fwg,fww,fgw,rsg,rsw,pb)
1979
1980	! ----------------------------------------------------------------------
1981	! This routine computes the effects of the reflections of short-wave
1982	! (solar) radiation in the street canyon by solving the system
1983	! of 2nz_u+1 eqn. in 2nz_u+1
1984	! unkonwn defined in A4, A5 and A6 of the paper (pages 295 and 296).
1985	! The system is solved by solving A X= B,
1986	! with A matrix B vector, and X solution.
1987	! ----------------------------------------------------------------------
1988
1989	implicit none
1990
1991
1992
1993	! ----------------------------------------------------------------------
1994	! INPUT:
1995	! ----------------------------------------------------------------------
1996	real albg ! Albedo of the ground for the current urban class
1997	real albw ! Albedo of the wall for the current urban class
1998	real fgw(nz_um,ndm,nurbm) ! View factors from ground to wall
1999	real fwg(nz_um,ndm,nurbm) ! View factors from wall to ground
2000	real fww(nz_um,nz_um,ndm,nurbm) ! View factors from wall to wall
2001	integer id ! current street direction
2002	integer iurb ! current urban class
2003	integer nz_u ! Number of layer in the urban grid
2004	real pb(nz_um) ! probability to have a building with an height equal
2005
2006	! ----------------------------------------------------------------------
2007	! OUTPUT:
2008	! ----------------------------------------------------------------------
2009	real rsg(ndm) ! Short wave radiation at the ground
2010	real rsw(2*ndm,nz_um) ! Short wave radiation at the walls
2011
2012	! ----------------------------------------------------------------------
2013	! LOCAL:
2014	! ----------------------------------------------------------------------
2015	integer i,j
2016	real aaa(2nz_um+1,2nz_um+1) ! terms of the matrix
2017	real bbb(2*nz_um+1) ! terms of the vector
2018
2019	! ----------------------------------------------------------------------
2020	! END VARIABLES DEFINITIONS
2021	! ----------------------------------------------------------------------
2022
2023
2024	! west wall
2025
2026	do i=1,nz_u
2027	do j=1,nz_u
2028	aaa(i,j)=0.
2029	enddo
2030
2031	aaa(i,i)=1.
2032
2033	do j=nz_u+1,2*nz_u
2034	aaa(i,j)=-albwfww(j-nz_u,i,id,iurb)pb(j-nz_u+1)
2035	enddo
2036
2037	aaa(i,2nz_u+1)=-albgfgw(i,id,iurb)
2038	bbb(i)=rsw(2*id-1,i)
2039
2040	enddo
2041
2042	! east wall
2043
2044	do i=1+nz_u,2*nz_u
2045	do j=1,nz_u
2046	aaa(i,j)=-albwfww(j,i-nz_u,id,iurb)pb(j+1)
2047	enddo
2048
2049	do j=1+nz_u,2*nz_u
2050	aaa(i,j)=0.
2051	enddo
2052
2053	aaa(i,i)=1.
2054	aaa(i,2nz_u+1)=-albgfgw(i-nz_u,id,iurb)
2055	bbb(i)=rsw(2*id,i-nz_u)
2056
2057	enddo
2058
2059	! ground
2060
2061	do j=1,nz_u
2062	aaa(2nz_u+1,j)=-albwfwg(j,id,iurb)*pb(j+1)
2063	enddo
2064
2065	do j=nz_u+1,2*nz_u
2066	aaa(2nz_u+1,j)=-albwfwg(j-nz_u,id,iurb)*pb(j-nz_u+1)
2067	enddo
2068
2069	aaa(2nz_u+1,2nz_u+1)=1.
2070	bbb(2*nz_u+1)=rsg(id)
2071
2072	call gaussj(aaa,2nz_u+1,bbb,2nz_um+1)
2073
2074	do i=1,nz_u
2075	rsw(2*id-1,i)=bbb(i)
2076	enddo
2077
2078	do i=nz_u+1,2*nz_u
2079	rsw(2*id,i-nz_u)=bbb(i)
2080	enddo
2081
2082	rsg(id)=bbb(2*nz_u+1)
2083
2084	return
2085	end subroutine short_rad
2086
2087
2088	! ===6=8===============================================================72
2089	! ===6=8===============================================================72
2090
2091	subroutine gaussj(a,n,b,np)
2092
2093	! ----------------------------------------------------------------------
2094	! This routine solve a linear system of n equations of the form
2095	! A X = B
2096	! where A is a matrix a(i,j)
2097	! B a vector and X the solution
2098	! In output b is replaced by the solution
2099	! ----------------------------------------------------------------------
2100
2101	implicit none
2102
2103	! ----------------------------------------------------------------------
2104	! INPUT:
2105	! ----------------------------------------------------------------------
2106	integer np
2107	real a(np,np)
2108
2109	! ----------------------------------------------------------------------
2110	! OUTPUT:
2111	! ----------------------------------------------------------------------
2112	real b(np)
2113
2114	! ----------------------------------------------------------------------
2115	! LOCAL:
2116	! ----------------------------------------------------------------------
2117	integer nmax
2118	parameter (nmax=150)
2119
2120	real big,dum
2121	integer i,icol,irow
2122	integer j,k,l,ll,n
2123	integer ipiv(nmax)
2124	real pivinv
2125
2126	! ----------------------------------------------------------------------
2127	! END VARIABLES DEFINITIONS
2128	! ----------------------------------------------------------------------
2129
2130	do j=1,n
2131	ipiv(j)=0.
2132	enddo
2133
2134	do i=1,n
2135	big=0.
2136	do j=1,n
2137	if(ipiv(j).ne.1)then
2138	do k=1,n
2139	if(ipiv(k).eq.0)then
2140	if(abs(a(j,k)).ge.big)then
2141	big=abs(a(j,k))
2142	irow=j
2143	icol=k
2144	endif
2145	elseif(ipiv(k).gt.1)then
2146	pause 'singular matrix in gaussj'
2147	endif
2148	enddo
2149	endif
2150	enddo
2151
2152	ipiv(icol)=ipiv(icol)+1
2153
2154	if(irow.ne.icol)then
2155	do l=1,n
2156	dum=a(irow,l)
2157	a(irow,l)=a(icol,l)
2158	a(icol,l)=dum
2159	enddo
2160
2161	dum=b(irow)
2162	b(irow)=b(icol)
2163	b(icol)=dum
2164
2165	endif
2166
2167	if(a(icol,icol).eq.0)pause 'singular matrix in gaussj'
2168
2169	pivinv=1./a(icol,icol)
2170	a(icol,icol)=1
2171
2172	do l=1,n
2173	a(icol,l)=a(icol,l)*pivinv
2174	enddo
2175
2176	b(icol)=b(icol)*pivinv
2177
2178	do ll=1,n
2179	if(ll.ne.icol)then
2180	dum=a(ll,icol)
2181	a(ll,icol)=0.
2182	do l=1,n
2183	a(ll,l)=a(ll,l)-a(icol,l)*dum
2184	enddo
2185
2186	b(ll)=b(ll)-b(icol)*dum
2187
2188	endif
2189	enddo
2190	enddo
2191
2192	return
2193	end subroutine gaussj
2194
2195	! ===6=8===============================================================72
2196	! ===6=8===============================================================72
2197
2198	subroutine soil_temp(nz,dz,temp,pt,ala,cs, &
2199	rs,rl,press,dt,em,alb,rt,sf,gf)
2200
2201	! ----------------------------------------------------------------------
2202	! This routine solves the Fourier diffusion equation for heat in
2203	! the material (wall, roof, or ground). Resolution is done implicitely.
2204	! Boundary conditions are:
2205	! - fixed temperature at the interior
2206	! - energy budget at the surface
2207	! ----------------------------------------------------------------------
2208
2209	implicit none
2210
2211
2212
2213	! ----------------------------------------------------------------------
2214	! INPUT:
2215	! ----------------------------------------------------------------------
2216	integer nz ! Number of layers
2217	real ala(nz) ! Thermal diffusivity in each layers [m^2 s^-1]
2218	real alb ! Albedo of the surface
2219	real cs(nz) ! Specific heat of the material [J m^3 K^-1]
2220	real dt ! Time step
2221	real em ! Emissivity of the surface
2222	real press ! Pressure at ground level
2223	real rl ! Downward flux of the longwave radiation
2224	real rs ! Solar radiation
2225	real sf ! Sensible heat flux at the surface
2226	real temp(nz) ! Temperature in each layer [K]
2227	real dz(nz) ! Layer sizes [m]
2228
2229
2230	! ----------------------------------------------------------------------
2231	! OUTPUT:
2232	! ----------------------------------------------------------------------
2233	real gf ! Heat flux transferred from the surface toward the interior
2234	real pt ! Potential temperature at the surface
2235	real rt ! Total radiation at the surface (solar+incoming long+outgoing long)
2236
2237	! ----------------------------------------------------------------------
2238	! LOCAL:
2239	! ----------------------------------------------------------------------
2240	integer iz
2241	real a(nz,3)
2242	real alpha
2243	real c(nz)
2244	real cddz(nz+2)
2245	real tsig
2246
2247	! ----------------------------------------------------------------------
2248	! END VARIABLES DEFINITIONS
2249	! ----------------------------------------------------------------------
2250
2251	tsig=temp(nz)
2252	alpha=(1.-alb)rs+emrl-emsigma(tsig**4)+sf
2253	! Compute cddz=2*cd/dz
2254
2255	cddz(1)=ala(1)/dz(1)
2256	do iz=2,nz
2257	cddz(iz)=2.*ala(iz)/(dz(iz)+dz(iz-1))
2258	enddo
2259	! cddz(nz+1)=ala(nz+1)/dz(nz)
2260
2261	a(1,1)=0.
2262	a(1,2)=1.
2263	a(1,3)=0.
2264	c(1)=temp(1)
2265
2266	do iz=2,nz-1
2267	a(iz,1)=-cddz(iz)*dt/dz(iz)
2268	a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dz(iz)
2269	a(iz,3)=-cddz(iz+1)*dt/dz(iz)
2270	c(iz)=temp(iz)
2271	enddo
2272
2273	a(nz,1)=-dt*cddz(nz)/dz(nz)
2274	a(nz,2)=1.+dt*cddz(nz)/dz(nz)
2275	a(nz,3)=0.
2276	c(nz)=temp(nz)+dt*alpha/cs(nz)/dz(nz)
2277
2278
2279	call invert(nz,a,c,temp)
2280
2281
2282	pt=temp(nz)(press/1.e+5)*(-rcp_u)
2283
2284	rt=(1.-alb)rs+emrl-emsigma(tsig**4)
2285
2286	! gf=-cddz(nz)(temp(nz)-temp(nz-1))cs(nz)
2287	gf=(1.-alb)rs+emrl-emsigma(tsig**4)+sf
2288	return
2289	end subroutine soil_temp
2290
2291	! ===6=8===============================================================72
2292	! ===6=8===============================================================72
2293
2294	subroutine invert(n,a,c,x)
2295
2296	! ----------------------------------------------------------------------
2297	! Inversion and resolution of a tridiagonal matrix
2298	! A X = C
2299	! ----------------------------------------------------------------------
2300
2301	implicit none
2302
2303	! ----------------------------------------------------------------------
2304	! INPUT:
2305	! ----------------------------------------------------------------------
2306	integer n
2307	real a(n,3) ! a(*,1) lower diagonal (Ai,i-1)
2308	! a(*,2) principal diagonal (Ai,i)
2309	! a(*,3) upper diagonal (Ai,i+1)
2310	real c(n)
2311
2312	! ----------------------------------------------------------------------
2313	! OUTPUT:
2314	! ----------------------------------------------------------------------
2315	real x(n)
2316
2317	! ----------------------------------------------------------------------
2318	! LOCAL:
2319	! ----------------------------------------------------------------------
2320	integer i
2321
2322	! ----------------------------------------------------------------------
2323	! END VARIABLES DEFINITIONS
2324	! ----------------------------------------------------------------------
2325
2326	do i=n-1,1,-1
2327	c(i)=c(i)-a(i,3)*c(i+1)/a(i+1,2)
2328	a(i,2)=a(i,2)-a(i,3)*a(i+1,1)/a(i+1,2)
2329	enddo
2330
2331	do i=2,n
2332	c(i)=c(i)-a(i,1)*c(i-1)/a(i-1,2)
2333	enddo
2334
2335	do i=1,n
2336	x(i)=c(i)/a(i,2)
2337	enddo
2338
2339	return
2340	end subroutine invert
2341
2342
2343	! ===6=8===============================================================72
2344	! ===6=8===============================================================72
2345
2346	subroutine flux_wall(ua,va,pt,da,ptw,uva,vva,uvb,vvb, &
2347	tva,tvb,evb,drst,dt)
2348
2349	! ----------------------------------------------------------------------
2350	! This routine computes the surface sources or sinks of momentum, tke,
2351	! and heat from vertical surfaces (walls).
2352	! ----------------------------------------------------------------------
2353
2354	implicit none
2355
2356
2357
2358	! INPUT:
2359	! -----
2360	real drst ! street directions for the current urban class
2361	real da ! air density
2362	real pt ! potential temperature
2363	real ptw ! Walls potential temperatures
2364	real ua ! wind speed
2365	real va ! wind speed
2366
2367	real dt !time step
2368	! OUTPUT:
2369	! ------
2370	! Explicit and implicit component of the momentum, temperature and TKE sources or sinks on
2371	! vertical surfaces (walls).
2372	! The fluxes can be computed as follow: Fluxes of X = A*X + B
2373	! Example: Momentum fluxes on vertical surfaces = uva_u * ua_u + uvb_u
2374	real uva ! U (wind component) Vertical surfaces, A (implicit) term
2375	real uvb ! U (wind component) Vertical surfaces, B (explicit) term
2376	real vva ! V (wind component) Vertical surfaces, A (implicit) term
2377	real vvb ! V (wind component) Vertical surfaces, B (explicit) term
2378	real tva ! Temperature Vertical surfaces, A (implicit) term
2379	real tvb ! Temperature Vertical surfaces, B (explicit) term
2380	real evb ! Energy (TKE) Vertical surfaces, B (explicit) term
2381
2382	! LOCAL:
2383	! -----
2384	real hc
2385	real u_ort
2386	real vett
2387
2388	! -------------------------
2389	! END VARIABLES DEFINITIONS
2390	! -------------------------
2391
2392	vett=(ua2+va2)**.5
2393
2394	u_ort=abs((cos(drst)ua-sin(drst)va))
2395
2396	uva=-cdragu_ort/2.cos(drst)*cos(drst)
2397	vva=-cdragu_ort/2.sin(drst)*sin(drst)
2398
2399	uvb=cdragu_ort/2.sin(drst)cos(drst)va
2400	vvb=cdragu_ort/2.sin(drst)cos(drst)ua
2401
2402	hc=5.678(1.09+0.23(vett/0.3048))
2403
2404	if(hc.gt.da*cp_u/dt)then
2405	hc=da*cp_u/dt
2406	endif
2407
2408	! tvb=hc*ptw/da/cp_u
2409	! tva=-hc/da/cp_u
2410	!!!!!!!!!!!!!!!!!!!!
2411	! explicit
2412	tvb=hcptw/da/cp_u-hc/da/cp_upt !c
2413	tva = 0. !c
2414
2415	evb=cdrag(abs(u_ort)*3.)/2.
2416
2417	return
2418	end subroutine flux_wall
2419
2420	! ===6=8===============================================================72
2421
2422	! ===6=8===============================================================72
2423
2424	subroutine flux_flat(dz,z0,ua,va,pt,pt0,ptg, &
2425	uhb,vhb,thb,ehb)
2426
2427	! ----------------------------------------------------------------------
2428	! Calculation of the flux at the ground
2429	! Formulation of Louis (Louis, 1979)
2430	! ----------------------------------------------------------------------
2431
2432	implicit none
2433
2434
2435
2436	! ----------------------------------------------------------------------
2437	! INPUT:
2438	! ----------------------------------------------------------------------
2439	real dz ! first vertical level
2440	real pt ! potential temperature
2441	real pt0 ! reference potential temperature
2442	real ptg ! ground potential temperature
2443	real ua ! wind speed
2444	real va ! wind speed
2445	real z0 ! Roughness length
2446
2447	! ----------------------------------------------------------------------
2448	! OUTPUT:
2449	! ----------------------------------------------------------------------
2450	! Explicit component of the momentum, temperature and TKE sources or sinks on horizontal
2451	! surfaces (roofs and street)
2452	! The fluxes can be computed as follow: Fluxes of X = B
2453	! Example: Momentum fluxes on horizontal surfaces = uhb_u
2454	real uhb ! U (wind component) Horizontal surfaces, B (explicit) term
2455	real vhb ! V (wind component) Horizontal surfaces, B (explicit) term
2456	real thb ! Temperature Horizontal surfaces, B (explicit) term
2457	real tva ! Temperature Vertical surfaces, A (implicit) term
2458	real tvb ! Temperature Vertical surfaces, B (explicit) term
2459	real ehb ! Energy (TKE) Horizontal surfaces, B (explicit) term
2460
2461
2462	! ----------------------------------------------------------------------
2463	! LOCAL:
2464	! ----------------------------------------------------------------------
2465	real aa
2466	real al
2467	real buu
2468	real c
2469	real fbuw
2470	real fbpt
2471	real fh
2472	real fm
2473	real ric
2474	real tstar
2475	real ustar
2476	real utot
2477	real wstar
2478	real zz
2479
2480	real b,cm,ch,rr,tol
2481	parameter(b=9.4,cm=7.4,ch=5.3,rr=0.74,tol=.001)
2482
2483	! ----------------------------------------------------------------------
2484	! END VARIABLES DEFINITIONS
2485	! ----------------------------------------------------------------------
2486
2487
2488	! computation of the ground temperature
2489
2490	utot=(ua2+va2)**.5
2491
2492
2493	!!!! Louis formulation
2494	!
2495	! compute the bulk Richardson Number
2496
2497	zz=dz/2.
2498
2499	! if(tstar.lt.0.)then
2500	! wstar=(-ustartstarghii/pt)*(1./3.)
2501	! else
2502	! wstar=0.
2503	! endif
2504	!
2505	! if (utot.le.0.7wstar) utot=max(0.7wstar,0.00001)
2506
2507	utot=max(utot,0.01)
2508
2509	ric=2.g_uzz(pt-ptg)/((pt+ptg)(utot**2))
2510
2511	aa=vk/log(zz/z0)
2512
2513	! determine the parameters fm and fh for stable, neutral and unstable conditions
2514
2515	if(ric.gt.0)then
2516	fm=1/(1+0.5bric)**2
2517	fh=fm
2518	else
2519	c=bcmaaaa(zz/z0)**.5
2520	fm=1-bric/(1+c(-ric)**.5)
2521	c=c*ch/cm
2522	fh=1-bric/(1+c(-ric)**.5)
2523	endif
2524
2525	fbuw=-aaaautotutotfm
2526	fbpt=-aaaautot(pt-ptg)fh/rr
2527
2528	ustar=(-fbuw)**.5
2529	tstar=-fbpt/ustar
2530
2531	al=(vkg_utstar)/(ptustarustar)
2532
2533	buu=-g_u/pt0ustartstar
2534
2535	uhb=-ustarustarua/utot
2536	vhb=-ustarustarva/utot
2537	thb=-ustar*tstar
2538	! thb= 0.
2539	ehb=buu
2540	!!!!!!!!!!!!!!!
2541
2542	return
2543	end subroutine flux_flat
2544
2545	! ===6=8===============================================================72
2546	! ===6=8===============================================================72
2547
2548	subroutine icBEP (alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u, &
2549	albg_u,albw_u,albr_u,emg_u,emw_u,emr_u, &
2550	fww,fwg,fgw,fsw,fws,fsg, &
2551	z0g_u,z0r_u, &
2552	nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b,ss_u,pb_u, &
2553	nz_u,z_u, &
2554	twini_u,trini_u)
2555
2556
2557	implicit none
2558
2559
2560	! Building parameters
2561	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
2562	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
2563	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
2564	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
2565	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
2566	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
2567	real twini_u(nurbm) ! Temperature inside the buildings behind the wall [K]
2568	real trini_u(nurbm) ! Temperature inside the buildings behind the roof [K]
2569
2570	! Radiation parameters
2571	real albg_u(nurbm) ! Albedo of the ground
2572	real albw_u(nurbm) ! Albedo of the wall
2573	real albr_u(nurbm) ! Albedo of the roof
2574	real emg_u(nurbm) ! Emissivity of ground
2575	real emw_u(nurbm) ! Emissivity of wall
2576	real emr_u(nurbm) ! Emissivity of roof
2577
2578	! Roughness parameters
2579	real z0g_u(nurbm) ! The ground's roughness length
2580	real z0r_u(nurbm) ! The roof's roughness length
2581
2582	! Street parameters
2583	integer nd_u(nurbm) ! Number of street direction for each urban class
2584
2585	real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
2586	real drst_u(ndm,nurbm) ! Street direction [degree]
2587	real ws_u(ndm,nurbm) ! Street width [m]
2588	real bs_u(ndm,nurbm) ! Building width [m]
2589	real h_b(nz_um,nurbm) ! Bulding's heights [m]
2590	real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
2591	! -----------------------------------------------------------------------
2592	! Output
2593	!------------------------------------------------------------------------
2594
2595
2596
2597	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
2598	! and the short wave radation. They are the part of radiation from a surface
2599	! or from the sky to another surface.
2600	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
2601	real fwg(nz_um,ndm,nurbm) ! from wall to ground
2602	real fgw(nz_um,ndm,nurbm) ! from ground to wall
2603	real fsw(nz_um,ndm,nurbm) ! from sky to wall
2604	real fws(nz_um,ndm,nurbm) ! from wall to sky
2605	real fsg(ndm,nurbm) ! from sky to ground
2606
2607	real ss_u(nz_um,nurbm) ! The probability that a building has an height equal to z
2608	real pb_u(nz_um,nurbm) ! The probability that a building has an height greater or equal to z
2609
2610	! Grid parameters
2611	integer nz_u(nurbm) ! Number of layer in the urban grid
2612	real z_u(nz_um) ! Height of the urban grid levels
2613
2614
2615	! -----------------------------------------------------------------------
2616	! Local
2617	!------------------------------------------------------------------------
2618
2619	integer iz_u,id,ilu,iurb
2620
2621	real dtot
2622	real hbmax
2623
2624	!------------------------------------------------------------------------
2625
2626
2627	! -----------------------------------------------------------------------
2628	! This routine initialise the urban paramters for the BEP module
2629	!------------------------------------------------------------------------
2630	!
2631	!Initialize variables
2632	!
2633	nz_u=0
2634	z_u=0.
2635	ss_u=0.
2636	pb_u=0.
2637	fww=0.
2638	fwg=0.
2639	fgw=0.
2640	fsw=0.
2641	fws=0.
2642	fsg=0.
2643
2644	! Computation of the urban levels height
2645
2646	z_u(1)=0.
2647
2648	do iz_u=1,nz_um-1
2649	z_u(iz_u+1)=z_u(iz_u)+dz_u
2650	enddo
2651
2652	! Normalisation of the building density
2653
2654	do iurb=1,nurbm
2655	dtot=0.
2656	do ilu=1,nz_um
2657	dtot=dtot+d_b(ilu,iurb)
2658	enddo
2659	do ilu=1,nz_um
2660	d_b(ilu,iurb)=d_b(ilu,iurb)/dtot
2661	enddo
2662	enddo
2663
2664	! Compute the view factors, pb and ss
2665
2666	do iurb=1,nurbm
2667	hbmax=0.
2668	nz_u(iurb)=0
2669	do ilu=1,nz_um
2670	if(h_b(ilu,iurb).gt.hbmax)hbmax=h_b(ilu,iurb)
2671	enddo
2672
2673	do iz_u=1,nz_um-1
2674	if(z_u(iz_u+1).gt.hbmax)go to 10
2675	enddo
2676
2677	10 continue
2678	nz_u(iurb)=iz_u+1
2679
2680	do id=1,nd_u(iurb)
2681
2682	call view_factors(iurb,nz_u(iurb),id,strd_u(id,iurb), &
2683	z_u,ws_u(id,iurb), &
2684	fww,fwg,fgw,fsg,fsw,fws)
2685
2686	do iz_u=1,nz_u(iurb)
2687	ss_u(iz_u,iurb)=0.
2688	do ilu=1,nz_um
2689	if(z_u(iz_u).le.h_b(ilu,iurb) &
2690	.and.z_u(iz_u+1).gt.h_b(ilu,iurb))then
2691	ss_u(iz_u,iurb)=ss_u(iz_u,iurb)+d_b(ilu,iurb)
2692	endif
2693	enddo
2694	enddo
2695
2696	pb_u(1,iurb)=1.
2697	do iz_u=1,nz_u(iurb)
2698	pb_u(iz_u+1,iurb)=max(0.,pb_u(iz_u,iurb)-ss_u(iz_u,iurb))
2699	enddo
2700
2701	enddo
2702	end do
2703
2704
2705	return
2706	end subroutine icBEP
2707
2708	! ===6=8===============================================================72
2709	! ===6=8===============================================================72
2710
2711	subroutine view_factors(iurb,nz_u,id,dxy,z,ws,fww,fwg,fgw,fsg,fsw,fws)
2712
2713	implicit none
2714
2715
2716
2717	! -----------------------------------------------------------------------
2718	! Input
2719	!------------------------------------------------------------------------
2720
2721	integer iurb ! Number of the urban class
2722	integer nz_u ! Number of levels in the urban grid
2723	integer id ! Street direction number
2724	real ws ! Street width
2725	real z(nz_um) ! Height of the urban grid levels
2726	real dxy ! Street lenght
2727
2728
2729	! -----------------------------------------------------------------------
2730	! Output
2731	!------------------------------------------------------------------------
2732
2733	! fww,fwg,fgw,fsw,fsg are the view factors used to compute the long wave
2734	! and the short wave radation. They are the part of radiation from a surface
2735	! or from the sky to another surface.
2736
2737	real fww(nz_um,nz_um,ndm,nurbm) ! from wall to wall
2738	real fwg(nz_um,ndm,nurbm) ! from wall to ground
2739	real fgw(nz_um,ndm,nurbm) ! from ground to wall
2740	real fsw(nz_um,ndm,nurbm) ! from sky to wall
2741	real fws(nz_um,ndm,nurbm) ! from wall to sky
2742	real fsg(ndm,nurbm) ! from sky to ground
2743
2744
2745	! -----------------------------------------------------------------------
2746	! Local
2747	!------------------------------------------------------------------------
2748
2749	integer jz,iz
2750
2751	real hut
2752	real f1,f2,f12,f23,f123,ftot
2753	real fprl,fnrm
2754	real a1,a2,a3,a4,a12,a23,a123
2755
2756	! -----------------------------------------------------------------------
2757	! This routine calculates the view factors
2758	!------------------------------------------------------------------------
2759
2760	hut=z(nz_u+1)
2761
2762	do jz=1,nz_u
2763
2764	! radiation from wall to wall
2765
2766	do iz=1,nz_u
2767
2768	call fprls (fprl,dxy,abs(z(jz+1)-z(iz )),ws)
2769	f123=fprl
2770	call fprls (fprl,dxy,abs(z(jz+1)-z(iz+1)),ws)
2771	f23=fprl
2772	call fprls (fprl,dxy,abs(z(jz )-z(iz )),ws)
2773	f12=fprl
2774	call fprls (fprl,dxy,abs(z(jz )-z(iz+1)),ws)
2775	f2 = fprl
2776
2777	a123=dxy*(abs(z(jz+1)-z(iz )))
2778	a12 =dxy*(abs(z(jz )-z(iz )))
2779	a23 =dxy*(abs(z(jz+1)-z(iz+1)))
2780	a1 =dxy*(abs(z(iz+1)-z(iz )))
2781	a2 =dxy*(abs(z(jz )-z(iz+1)))
2782	a3 =dxy*(abs(z(jz+1)-z(jz )))
2783
2784	ftot=0.5(a123f123-a23f23-a12f12+a2*f2)/a1
2785
2786	fww(iz,jz,id,iurb)=ftot*a1/a3
2787
2788	enddo
2789
2790	! radiation from ground to wall
2791
2792	call fnrms (fnrm,z(jz+1),dxy,ws)
2793	f12=fnrm
2794	call fnrms (fnrm,z(jz) ,dxy,ws)
2795	f1=fnrm
2796
2797	a1 = ws*dxy
2798
2799	a12= ws*dxy
2800
2801	a4=(z(jz+1)-z(jz))*dxy
2802
2803	ftot=(a12f12-a12f1)/a1
2804
2805	fgw(jz,id,iurb)=ftot*a1/a4
2806
2807	! radiation from sky to wall
2808
2809	call fnrms(fnrm,hut-z(jz) ,dxy,ws)
2810	f12 = fnrm
2811	call fnrms (fnrm,hut-z(jz+1),dxy,ws)
2812	f1 =fnrm
2813
2814	a1 = ws*dxy
2815
2816	a12= ws*dxy
2817
2818	a4 = (z(jz+1)-z(jz))*dxy
2819
2820	ftot=(a12f12-a12f1)/a1
2821
2822	fsw(jz,id,iurb)=ftot*a1/a4
2823
2824	enddo
2825
2826	! radiation from wall to sky
2827	do iz=1,nz_u
2828	call fnrms(fnrm,ws,dxy,hut-z(iz))
2829	f12=fnrm
2830	call fnrms(fnrm,ws,dxy,hut-z(iz+1))
2831	f1=fnrm
2832	a1 = (z(iz+1)-z(iz))*dxy
2833	a2 = (hut-z(iz+1))*dxy
2834	a12= (hut-z(iz))*dxy
2835	a4 = ws*dxy
2836	ftot=(a12f12-a2f1)/a1
2837	fws(iz,id,iurb)=ftot*a1/a4
2838
2839	enddo
2840	!!!!!!!!!!!!!
2841
2842
2843	do iz=1,nz_u
2844
2845	! radiation from wall to ground
2846
2847	call fnrms (fnrm,ws,dxy,z(iz+1))
2848	f12=fnrm
2849	call fnrms (fnrm,ws,dxy,z(iz ))
2850	f1 =fnrm
2851
2852	a1= (z(iz+1)-z(iz) )*dxy
2853
2854	a2 = z(iz)*dxy
2855	a12= z(iz+1)*dxy
2856	a4 = ws*dxy
2857
2858	ftot=(a12f12-a2f1)/a1
2859
2860	fwg(iz,id,iurb)=ftot*a1/a4
2861
2862	enddo
2863
2864	! radiation from sky to ground
2865
2866	call fprls (fprl,dxy,ws,hut)
2867	fsg(id,iurb)=fprl
2868
2869	return
2870	end subroutine view_factors
2871
2872	! ===6=8===============================================================72
2873	! ===6=8===============================================================72
2874
2875	SUBROUTINE fprls (fprl,a,b,c)
2876
2877	implicit none
2878
2879
2880
2881	real a,b,c
2882	real x,y
2883	real fprl
2884
2885
2886	x=a/c
2887	y=b/c
2888
2889	if(a.eq.0.or.b.eq.0.)then
2890	fprl=0.
2891	else
2892	fprl=log( ( (1.+x*2)(1.+y2)/(1.+x2+y2) ).5)+ &
2893	y((1.+x2).5)atan(y/((1.+x2).5))+ &
2894	x((1.+y2).5)atan(x/((1.+y2).5))- &
2895	yatan(y)-xatan(x)
2896	fprl=fprl2./(pix*y)
2897	endif
2898
2899	return
2900	end subroutine fprls
2901
2902	! ===6=8===============================================================72
2903	! ===6=8===============================================================72
2904
2905	SUBROUTINE fnrms (fnrm,a,b,c)
2906
2907	implicit none
2908
2909
2910
2911	real a,b,c
2912	real x,y,z,a1,a2,a3,a4,a5,a6
2913	real fnrm
2914
2915	x=a/b
2916	y=c/b
2917	z=x2+y2
2918
2919	if(y.eq.0.or.x.eq.0)then
2920	fnrm=0.
2921	else
2922	a1=log( (1.+xx)(1.+y*y)/(1.+z) )
2923	a2=yylog(yy(1.+z)/z/(1.+y*y) )
2924	a3=xxlog(xx(1.+z)/z/(1.+x*x) )
2925	a4=y*atan(1./y)
2926	a5=x*atan(1./x)
2927	a6=sqrt(z)*atan(1./sqrt(z))
2928	fnrm=0.25*(a1+a2+a3)+a4+a5-a6
2929	fnrm=fnrm/(pi*y)
2930	endif
2931
2932	return
2933	end subroutine fnrms
2934	! ===6=8===============================================================72
2935
2936	SUBROUTINE init_para(alag_u,alaw_u,alar_u,csg_u,csw_u,csr_u,&
2937	twini_u,trini_u,tgini_u,albg_u,albw_u,albr_u,emg_u,emw_u,&
2938	emr_u,z0g_u,z0r_u,nd_u,strd_u,drst_u,ws_u,bs_u,h_b,d_b)
2939
2940	! initialization routine, where the variables from the table are read
2941
2942	implicit none
2943
2944	integer iurb ! urban class number
2945	! Building parameters
2946	real alag_u(nurbm) ! Ground thermal diffusivity [m^2 s^-1]
2947	real alaw_u(nurbm) ! Wall thermal diffusivity [m^2 s^-1]
2948	real alar_u(nurbm) ! Roof thermal diffusivity [m^2 s^-1]
2949	real csg_u(nurbm) ! Specific heat of the ground material [J m^3 K^-1]
2950	real csw_u(nurbm) ! Specific heat of the wall material [J m^3 K^-1]
2951	real csr_u(nurbm) ! Specific heat of the roof material [J m^3 K^-1]
2952	real twini_u(nurbm) ! Temperature inside the buildings behind the wall [K]
2953	real trini_u(nurbm) ! Temperature inside the buildings behind the roof [K]
2954	real tgini_u(nurbm) ! Initial road temperature
2955
2956	! Radiation parameters
2957	real albg_u(nurbm) ! Albedo of the ground
2958	real albw_u(nurbm) ! Albedo of the wall
2959	real albr_u(nurbm) ! Albedo of the roof
2960	real emg_u(nurbm) ! Emissivity of ground
2961	real emw_u(nurbm) ! Emissivity of wall
2962	real emr_u(nurbm) ! Emissivity of roof
2963
2964	! Roughness parameters
2965	real z0g_u(nurbm) ! The ground's roughness length
2966	real z0r_u(nurbm) ! The roof's roughness length
2967
2968	! Street parameters
2969	integer nd_u(nurbm) ! Number of street direction for each urban class
2970
2971	real strd_u(ndm,nurbm) ! Street length (fix to greater value to the horizontal length of the cells)
2972	real drst_u(ndm,nurbm) ! Street direction [degree]
2973	real ws_u(ndm,nurbm) ! Street width [m]
2974	real bs_u(ndm,nurbm) ! Building width [m]
2975	real h_b(nz_um,nurbm) ! Bulding's heights [m]
2976	real d_b(nz_um,nurbm) ! The probability that a building has an height h_b
2977
2978	integer i,iu
2979	integer nurb ! number of urban classes used
2980
2981	!
2982	!Initialize some variables
2983	!
2984
2985	h_b=0.
2986	d_b=0.
2987
2988	nurb=ICATE
2989	do iu=1,nurb
2990	nd_u(iu)=0
2991	enddo
2992
2993	csw_u=CAPB_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
2994	csr_u=CAPR_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
2995	csg_u=CAPG_TBL / (( 1.0 / 4.1868 ) * 1.E-6)
2996	do i=1,icate
2997	alaw_u(i)=AKSB_TBL(i) / csw_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
2998	alar_u(i)=AKSR_TBL(i) / csr_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
2999	alag_u(i)=AKSG_TBL(i) / csg_u(i) / (( 1.0 / 4.1868 ) * 1.E-2)
3000	enddo
3001	twini_u=TBLEND_TBL
3002	trini_u=TRLEND_TBL
3003	tgini_u=TGLEND_TBL
3004	albw_u=ALBB_TBL
3005	albr_u=ALBR_TBL
3006	albg_u=ALBG_TBL
3007	emw_u=EPSB_TBL
3008	emr_u=EPSR_TBL
3009	emg_u=EPSG_TBL
3010	z0r_u=Z0R_TBL
3011	z0g_u=Z0G_TBL
3012	nd_u=NUMDIR_TBL
3013	do iu=1,icate
3014	if(ndm.lt.nd_u(iu))then
3015	write(,)'ndm too small in module_sf_bep, please increase to at least ', nd_u(iu)
3016	write(,)'remember also that num_urban_layers should be equal or greater than nz_umndmnwr-u!'
3017	stop
3018	endif
3019	do i=1,nd_u(iu)
3020	drst_u(i,iu)=STREET_DIRECTION_TBL(i,iu) * pi/180.
3021	ws_u(i,iu)=STREET_WIDTH_TBL(i,iu)
3022	bs_u(i,iu)=BUILDING_WIDTH_TBL(i,iu)
3023	enddo
3024	enddo
3025	do iu=1,ICATE
3026	if(nz_um.lt.numhgt_tbl(iu)+3)then
3027	write(,)'nz_um too small in module_sf_bep, please increase to at least ',numhgt_tbl(iu)+3
3028	write(,)'remember also that num_urban_layers should be equal or greater than nz_umndmnwr-u!'
3029	stop
3030	endif
3031	do i=1,NUMHGT_TBL(iu)
3032	h_b(i,iu)=HEIGHT_BIN_TBL(i,iu)
3033	d_b(i,iu)=HPERCENT_BIN_TBL(i,iu)
3034	enddo
3035	enddo
3036
3037	do i=1,ndm
3038	do iu=1,nurbm
3039	strd_u(i,iu)=100000.
3040	enddo
3041	enddo
3042
3043	return
3044	END SUBROUTINE init_para
3045	!==============================================================
3046
3047	!==============================================================
3048	subroutine angle(along,alat,day,realt,zr,deltar,ah)
3049	! ----------------
3050	!
3051	! Computation of the solar angles
3052	! schayes (1982,atm. env. , p1407)
3053	! Inputs
3054	!========================
3055	! along=longitud
3056	! alat=latitude
3057	! day=julian day (from the beginning of the year)
3058	! realt= time GMT in hours
3059	! Outputs
3060	!============================
3061	! zr=solar zenith angle
3062	! deltar=declination angle
3063	! ah=hour angle
3064	!===============================
3065
3066	implicit none
3067	real along,alat, realt, zr, deltar, ah, arg
3068	real rad,om,radh,initt, pii, drad, alongt, cphi, sphi
3069	real c1, c2, c3, s1, s2, s3, delta, rmsr2, cd, sid
3070	real et, ahor, chor, coznt
3071	integer day
3072
3073
3074	data rad,om,radh,initt/0.0174533,0.0172142,0.26179939,0/
3075
3076	zr=0.
3077	deltar=0.
3078	ah=0.
3079
3080	pii = 3.14159265358979312
3081	drad = pii/180.
3082
3083	alongt=along/15.
3084	cphi=cos(alat*drad)
3085	sphi=sin(alat*drad)
3086	!
3087	! declination
3088	!
3089	arg=om*day
3090	c1=cos(arg)
3091	c2=cos(2.*arg)
3092	c3=cos(3.*arg)
3093	s1=sin(arg)
3094	s2=sin(2.*arg)
3095	s3=sin(3.*arg)
3096	delta=0.33281-22.984c1-0.3499c2-0.1398c3+3.7872s1+0.03205s2+0.07187s3
3097	rmsr2=(1./(1.-0.01673c1))*2
3098	deltar=delta*rad
3099	cd=cos(deltar)
3100	sid=sin(deltar)
3101	!
3102	! time equation in hours
3103	!
3104	et=0.0072c1-0.0528c2-0.0012c3-0.1229s1-0.1565s2-0.0041s3
3105	!
3106	!
3107	! hour angle
3108	!
3109
3110	! ifh=0
3111
3112	! ahor=realt-12.+ifh+et+alongt
3113	ahor=realt-12.+et+alongt
3114	ah=ahor*radh
3115	chor=cos(ah)
3116	!
3117	! zenith angle
3118	!
3119	coznt=sphisid+cphicd*chor
3120
3121	zr=acos(coznt)
3122
3123	return
3124
3125	END SUBROUTINE angle
3126	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3127
3128	subroutine upward_rad(nd_u,iurb,nz_u,ws,bs,sigma,fsw,fsg,pb,ss, &
3129	tg,emg_u,albg_u,rlg,rsg,sfg, &
3130	tw,emw_u,albw_u,rlw,rsw,sfw, &
3131	tr,emr_u,albr_u,rld,rs, sfr, &
3132	rs_abs,rl_up,emiss,grdflx_urb)
3133	!
3134	! IN this surboutine we compute the upward longwave flux, and the albedo
3135	! needed for the radiation scheme
3136	!
3137	implicit none
3138
3139	!
3140	!INPUT VARIABLES
3141	!
3142	real rsw(2*ndm,nz_um) ! Short wave radiation at the wall for a given canyon direction [W/m2]
3143	real rlw(2*ndm,nz_um) ! Long wave radiation at the walls for a given canyon direction [W/m2]
3144	real rsg(ndm) ! Short wave radiation at the canyon for a given canyon direction [W/m2]
3145	real rlg(ndm) ! Long wave radiation at the ground for a given canyon direction [W/m2]
3146	real rs ! Short wave radiation at the horizontal surface from the sun [W/mÂ²]
3147	real sfw(2*ndm,nz_um) ! Sensible heat flux from walls [W/mÂ²]
3148	real sfg(ndm) ! Sensible heat flux from ground (road) [W/mÂ²]
3149	real sfr(ndm,nz_um) ! Sensible heat flux from roofs [W/mÂ²]
3150	real rld ! Long wave radiation from the sky [W/mÂ²]
3151	real albg_u ! albedo of the ground/street
3152	real albw_u ! albedo of the walls
3153	real albr_u ! albedo of the roof
3154	real ws(ndm) ! width of the street
3155	real bs(ndm)
3156	! building size
3157	real pb(nz_um) ! Probability to have a building with an height equal or higher
3158	integer nz_u
3159	real ss(nz_um) ! Probability to have a building of a given height
3160	real sigma
3161	real emg_u ! emissivity of the street
3162	real emw_u ! emissivity of the wall
3163	real emr_u ! emissivity of the roof
3164	real fsw(nz_um,ndm,nurbm) ! View factors from sky to wall
3165	real fsg(ndm,nurbm) ! groud to sky view factor
3166	real tw(2*ndm,nz_um,nwr_u) ! Temperature in each layer of the wall [K]
3167	real tr(ndm,nz_um,nwr_u) ! Temperature in each layer of the roof [K]
3168	real tg(ndm,ng_u) ! Temperature in each layer of the ground [K]
3169	integer iurb ! urban class
3170	integer id ! street direction
3171	integer nd_u ! number of street directions
3172	!OUTPUT/INPUT
3173	real rs_abs ! absrobed solar radiationfor this street direction
3174	real rl_up ! upward longwave radiation for this street direction
3175	real emiss ! mean emissivity
3176	real grdflx_urb ! ground heat flux
3177	!LOCAL
3178	integer iz,iw
3179	real rl_inc,rl_emit
3180	real gfl
3181	integer ix,iy,iwrong
3182
3183	iwrong=1
3184	do iz=1,nz_u+1
3185	do id=1,nd_u
3186	do iw=1,nwr_u
3187	if(tr(id,iz,iw).lt.100.)then
3188	write(,)'in upward_rad ',iz,id,iw,tr(id,iz,iw)
3189	iwrong=0
3190	endif
3191	enddo
3192	enddo
3193	enddo
3194	if(iwrong.eq.0)stop
3195
3196	rl_up=0.
3197
3198	rs_abs=0.
3199	rl_inc=0.
3200	emiss=0.
3201	rl_emit=0.
3202	grdflx_urb=0.
3203	do id=1,nd_u
3204	rl_emit=rl_emit-( emg_usigma(tg(id,ng_u)*4.)+(1-emg_u)rlg(id))*ws(id)/(ws(id)+bs(id))/nd_u
3205	rl_inc=rl_inc+rlg(id)*ws(id)/(ws(id)+bs(id))/nd_u
3206	rs_abs=rs_abs+(1.-albg_u)rsg(id)ws(id)/(ws(id)+bs(id))/nd_u
3207	gfl=(1.-albg_u)rsg(id)+emg_urlg(id)-emg_usigma(tg(id,ng_u)**4.)+sfg(id)
3208	grdflx_urb=grdflx_urb-gfl*ws(id)/(ws(id)+bs(id))/nd_u
3209
3210	do iz=2,nz_u
3211	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
3212	rl_inc=rl_inc+rldss(iz)bs(id)/(ws(id)+bs(id))/nd_u
3213	rs_abs=rs_abs+(1.-albr_u)rsss(iz)*bs(id)/(ws(id)+bs(id))/nd_u
3214	gfl=(1.-albr_u)rs+emr_urld-emr_usigma(tr(id,iz,nwr_u)**4.)+sfr(id,iz)
3215	grdflx_urb=grdflx_urb-gflss(iz)bs(id)/(ws(id)+bs(id))/nd_u
3216	enddo
3217
3218	do iz=1,nz_u
3219	rl_emit=rl_emit-(emw_usigma( tw(2id-1,iz,nwr_u)4.+tw(2id,iz,nwr_u)**4. )+ &
3220	(1-emw_u)( rlw(2id-1,iz)+rlw(2id,iz) ) )dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
3221	rl_inc=rl_inc+(( rlw(2id-1,iz)+rlw(2id,iz) ) )dz_upb(iz+1)/(ws(id)+bs(id))/nd_u
3222	rs_abs=rs_abs+((1.-albw_u)( rsw(2id-1,iz)+rsw(2id,iz) ) )dz_u*pb(iz+1)/(ws(id)+bs(id))/nd_u
3223	gfl=(1.-albw_u)(rsw(2id-1,iz)+rsw(2id,iz)) +emw_u( rlw(2id-1,iz)+rlw(2id,iz) ) &
3224	-emw_usigma( tw(2id-1,iz,nwr_u)4.+tw(2id,iz,nwr_u)*4. )+(sfw(2id-1,iz)+sfw(2*id,iz))
3225	grdflx_urb=grdflx_urb-gfldz_upb(iz+1)/(ws(id)+bs(id))/nd_u
3226	enddo
3227
3228	enddo
3229	emiss=(emg_u+emw_u+emr_u)/3.
3230	rl_up=(rl_inc+rl_emit)-rld
3231
3232
3233	return
3234
3235	END SUBROUTINE upward_rad
3236
3237	!====6=8===============================================================72
3238	!====6=8===============================================================72
3239	END MODULE module_sf_bep

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

Download in other formats:

Original Format