Home

Context Navigation

urban_surface_mod.f90 @ 2766

Last change on this file since 2766 was 2766, checked in by kanani, 3 years ago
Removal of chem directive, plus minor changes
Property svn:keywords set to `Id`
File size: 410.3 KB

Line
1	!> @file urban_surface_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2015-2018 Czech Technical University in Prague
18	! Copyright 1997-2018 Leibniz Universitaet Hannover
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! ------------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: urban_surface_mod.f90 2766 2018-01-22 17:17:47Z kanani $
28	! Removed redundant commas, added some blanks
29	!
30	! 2765 2018-01-22 11:34:58Z maronga
31	! Major bugfix in calculation of f_shf. Adjustment of roughness lengths in
32	! building_pars
33	!
34	! 2750 2018-01-15 16:26:51Z knoop
35	! Move flag plant canopy to modules
36	!
37	! 2737 2018-01-11 14:58:11Z kanani
38	! Removed unused variables t_surf_whole...
39	!
40	! 2735 2018-01-11 12:01:27Z suehring
41	! resistances are saved in surface attributes
42	!
43	! 2723 2018-01-05 09:27:03Z maronga
44	! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
45	!
46	! 2720 2018-01-02 16:27:15Z kanani
47	! Correction of comment
48	!
49	! 2718 2018-01-02 08:49:38Z maronga
50	! Corrected "Former revisions" section
51	!
52	! 2705 2017-12-18 11:26:23Z maronga
53	! Changes from last commit documented
54	!
55	! 2703 2017-12-15 20:12:38Z maronga
56	! Workaround for calculation of r_a
57	!
58	! 2696 2017-12-14 17:12:51Z kanani
59	! - Change in file header (GPL part)
60	! - Bugfix in calculation of pt_surface and related fluxes. (BM)
61	! - Do not write surface temperatures onto pt array as this might cause
62	! problems with nesting. (MS)
63	! - Revised calculation of pt1 (now done in surface_layer_fluxes).
64	! Bugfix, f_shf_window and f_shf_green were not set at vertical surface
65	! elements. (MS)
66	! - merged with branch ebsolver
67	! green building surfaces do not evaporate yet
68	! properties of green wall layers and window layers are taken from wall layers
69	! this input data is missing. (RvT)
70	! - Merged with branch radiation (developed by Mohamed Salim)
71	! - Revised initialization. (MS)
72	! - Rename emiss_surf into emissivity, roughness_wall into z0, albedo_surf into
73	! albedo. (MS)
74	! - Move first call of usm_radiatin from usm_init to init_3d_model
75	! - fixed problem with near surface temperature
76	! - added near surface temperature t_surf_10cm_h(m), t_surf_10cm_v(l)%t(m)
77	! - does not work with temp profile including stability, ol
78	! t_surf_10cm = pt1 now
79	! - merged with 2357 bugfix, error message for nopointer version
80	! - added indoor model coupling with wall heat flux
81	! - added green substrate/ dry vegetation layer for buildings
82	! - merged with 2232 new surface-type structure
83	! - added transmissivity of window tiles
84	! - added MOSAIK tile approach for 3 different surfaces (RvT)
85	!
86	! 2583 2017-10-26 13:58:38Z knoop
87	! Bugfix: reverted MPI_Win_allocate_cptr introduction in last commit
88	!
89	! 2582 2017-10-26 13:19:46Z hellstea
90	! Workaround for gnufortran compiler added in usm_calc_svf. CALL MPI_Win_allocate is
91	! replaced by CALL MPI_Win_allocate_cptr if defined ( __gnufortran ).
92	!
93	! 2544 2017-10-13 18:09:32Z maronga
94	! Date and time quantities are now read from date_and_time_mod. Solar constant is
95	! read from radiation_model_mod
96	!
97	! 2516 2017-10-04 11:03:04Z suehring
98	! Remove tabs
99	!
100	! 2514 2017-10-04 09:52:37Z suehring
101	! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
102	! no output of ghost layer data
103	!
104	! 2350 2017-08-15 11:48:26Z kanani
105	! Bugfix and error message for nopointer version.
106	! Additional "! defined(__nopointer)" as workaround to enable compilation of
107	! nopointer version.
108	!
109	! 2318 2017-07-20 17:27:44Z suehring
110	! Get topography top index via Function call
111	!
112	! 2317 2017-07-20 17:27:19Z suehring
113	! Bugfix: adjust output of shf. Added support for spinups
114	!
115	! 2287 2017-06-15 16:46:30Z suehring
116	! Bugfix in determination topography-top index
117	!
118	! 2269 2017-06-09 11:57:32Z suehring
119	! Enable restart runs with different number of PEs
120	! Bugfixes nopointer branch
121	!
122	! 2258 2017-06-08 07:55:13Z suehring
123	! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
124	!
125	! 2233 2017-05-30 18:08:54Z suehring
126	!
127	! 2232 2017-05-30 17:47:52Z suehring
128	! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
129	! insteat, heat fluxes are directly applied in diffusion_s.
130	!
131	! 2213 2017-04-24 15:10:35Z kanani
132	! Removal of output quantities usm_lad and usm_canopy_hr
133	!
134	! 2209 2017-04-19 09:34:46Z kanani
135	! cpp switch __mpi3 removed,
136	! minor formatting,
137	! small bugfix for division by zero (Krc)
138	!
139	! 2113 2017-01-12 13:40:46Z kanani
140	! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
141	!
142	! 2071 2016-11-17 11:22:14Z maronga
143	! Small bugfix (Resler)
144	!
145	! 2031 2016-10-21 15:11:58Z knoop
146	! renamed variable rho to rho_ocean
147	!
148	! 2024 2016-10-12 16:42:37Z kanani
149	! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
150	! optimization of MPI-RMA operations,
151	! declaration of pcbl as integer,
152	! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
153	! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
154	! use of new control parameter varnamelength,
155	! added output variables usm_rad_ressw, usm_rad_reslw,
156	! minor formatting changes,
157	! minor optimizations.
158	!
159	! 2011 2016-09-19 17:29:57Z kanani
160	! Major reformatting according to PALM coding standard (comments, blanks,
161	! alphabetical ordering, etc.),
162	! removed debug_prints,
163	! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
164	! defined in MODULE control_parameters (modules.f90) to avoid circular
165	! dependencies,
166	! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
167	!
168	! 2007 2016-08-24 15:47:17Z kanani
169	! Initial revision
170	!
171	!
172	! Description:
173	! ------------
174	! 2016/6/9 - Initial version of the USM (Urban Surface Model)
175	! authors: Jaroslav Resler, Pavel Krc
176	! (Czech Technical University in Prague and Institute of
177	! Computer Science of the Czech Academy of Sciences, Prague)
178	! with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
179	! partly inspired by PALM LSM (B. Maronga)
180	! parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
181	!> Module for Urban Surface Model (USM)
182	!> The module includes:
183	!> 1. radiation model with direct/diffuse radiation, shading, reflections
184	!> and integration with plant canopy
185	!> 2. wall and wall surface model
186	!> 3. surface layer energy balance
187	!> 4. anthropogenic heat (only from transportation so far)
188	!> 5. necessary auxiliary subroutines (reading inputs, writing outputs,
189	!> restart simulations, ...)
190	!> It also make use of standard radiation and integrates it into
191	!> urban surface model.
192	!>
193	!> Further work:
194	!> -------------
195	!> 1. Reduce number of shape view factors by merging factors for distant surfaces
196	!> under shallow angles. Idea: Iteratively select the smallest shape view
197	!> factor by value (among all sources and targets) which has a similarly
198	!> oriented source neighbor (or near enough) SVF and merge them by adding
199	!> value of the smaller SVF to the larger one and deleting the smaller one.
200	!> This will allow for better scaling at higher resolutions.
201	!>
202	!> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity
203	!> from surfaces that are visible from local surfaces (i.e. there is a SVF
204	!> where target is local). To do that, radiosity will be exchanged after each
205	!> reflection step using MPI_Alltoall instead of current MPI_Allgather.
206	!>
207	!> 3. Temporarily large values of surface heat flux can be observed, up to
208	!> 1.2 Km/s, which seem to be not realistic.
209	!>
210	!> @todo Output of _av variables in case of restarts
211	!> @todo Revise flux conversion in energy-balance solver
212	!> @todo Bugfixing in nopointer branch
213	!> @todo Check optimizations for RMA operations
214	!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
215	!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
216	!> factor 3 between min and max time
217	!> @todo Move setting of flag indoor_model to indoor_model_mod once available
218	!> @todo Check divisions in wtend (etc.) calculations for possible division
219	!> by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
220	!> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
221	!------------------------------------------------------------------------------!
222	MODULE urban_surface_mod
223
224	#if ! defined( __nopointer )
225	USE arrays_3d, &
226	ONLY: zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
227	#endif
228
229	USE cloud_parameters, &
230	ONLY: cp, r_d
231
232	USE constants, &
233	ONLY: pi
234
235	USE control_parameters, &
236	ONLY: coupling_start_time, dz, topography, dt_3d, &
237	intermediate_timestep_count, initializing_actions, &
238	intermediate_timestep_count_max, simulated_time, end_time, &
239	timestep_scheme, tsc, coupling_char, io_blocks, io_group, &
240	message_string, time_since_reference_point, surface_pressure, &
241	g, pt_surface, large_scale_forcing, lsf_surf, spinup, &
242	spinup_pt_mean, spinup_time, time_do3d, dt_do3d, &
243	average_count_3d, varnamelength, urban_surface, kappa, &
244	plant_canopy
245
246	USE cpulog, &
247	ONLY: cpu_log, log_point, log_point_s
248
249	USE date_and_time_mod, &
250	ONLY: d_seconds_year, day_of_year_init, time_utc_init
251
252	USE grid_variables, &
253	ONLY: dx, dy, ddx, ddy, ddx2, ddy2
254
255	USE indices, &
256	ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, &
257	nysg, nzb, nzt, nbgp, wall_flags_0
258
259	USE, INTRINSIC :: iso_c_binding
260
261	USE kinds
262
263	USE pegrid
264
265	USE plant_canopy_model_mod, &
266	ONLY: pc_heating_rate, usm_lad_rma
267
268	USE radiation_model_mod, &
269	ONLY: albedo_type, radiation, calc_zenith, zenith, &
270	rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, &
271	sigma_sb, solar_constant, sun_direction, sun_dir_lat, &
272	sun_dir_lon, &
273	force_radiation_call, surfinsw, surfinlw, surfinswdir, &
274	surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
275	surfinl, surfinlwdif, energy_balance_surf_h, &
276	energy_balance_surf_v, rad_sw_in_dir, rad_sw_in_diff, &
277	rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
278	surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, pcbl, npcbl, startenergy, &
279	endenergy, nenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, &
280	inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,id, &
281	iz, iy, ix, idir, jdir, kdir, startborder, endborder, nsurf_type, nzub, nzut, &
282	isky, inorth_b,idown_a, isouth_b, ieast_b, iwest_b, nzu, pch, nsurf, &
283	iup_a, inorth_a, isouth_a, ieast_a, iwest_a, idsvf, ndsvf, &
284	idcsf, ndcsf, kdcsf, pct, startland, endland, startwall, endwall
285
286	USE statistics, &
287	ONLY: hom, statistic_regions
288
289	USE surface_mod
290
291
292	IMPLICIT NONE
293
294
295	!-- configuration parameters (they can be setup in PALM config)
296	LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used
297	LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used
298	LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls
299	LOGICAL :: indoor_model = .FALSE. !< whether to use the indoor model
300
301
302	INTEGER(iwp) :: building_type = 1 !< default building type (preleminary setting)
303	INTEGER(iwp) :: land_category = 2 !< default category for land surface
304	INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone
305	INTEGER(iwp) :: pedestrant_category = 2 !< default category for wall surface in pedestrian zone
306	INTEGER(iwp) :: roof_category = 2 !< default category for root surface
307
308	!
309	!-- Indices of input attributes for (above) ground floor level
310	INTEGER(iwp) :: ind_alb_wall = 38 !< index in input list for albedo_type of wall fraction
311	INTEGER(iwp) :: ind_alb_green = 39 !< index in input list for albedo_type of green fraction
312	INTEGER(iwp) :: ind_alb_win = 40 !< index in input list for albedo_type of window fraction
313	INTEGER(iwp) :: ind_emis_wall_agfl = 14 !< index in input list for wall emissivity, above ground floor level
314	INTEGER(iwp) :: ind_emis_wall_gfl = 32 !< index in input list for wall emissivity, ground floor level
315	INTEGER(iwp) :: ind_emis_green_agfl = 15 !< index in input list for green emissivity, above ground floor level
316	INTEGER(iwp) :: ind_emis_green_gfl = 33 !< index in input list for green emissivity, ground floor level
317	INTEGER(iwp) :: ind_emis_win_agfl = 16 !< index in input list for window emissivity, above ground floor level
318	INTEGER(iwp) :: ind_emis_win_gfl = 34 !< index in input list for window emissivity, ground floor level
319	INTEGER(iwp) :: ind_green_frac_w_agfl = 2 !< index in input list for green fraction on wall, above ground floor level
320	INTEGER(iwp) :: ind_green_frac_w_gfl = 23 !< index in input list for green fraction on wall, ground floor level
321	INTEGER(iwp) :: ind_green_frac_r_agfl = 3 !< index in input list for green fraction on roof, above ground floor level
322	INTEGER(iwp) :: ind_green_frac_r_gfl = 24 !< index in input list for green fraction on roof, ground floor level
323	INTEGER(iwp) :: ind_hc1_agfl = 6 !< index in input list for heat capacity at first wall layer, above ground floor level
324	INTEGER(iwp) :: ind_hc1_gfl = 26 !< index in input list for heat capacity at first wall layer, ground floor level
325	INTEGER(iwp) :: ind_hc2_agfl = 7 !< index in input list for heat capacity at second wall layer, above ground floor level
326	INTEGER(iwp) :: ind_hc2_gfl = 27 !< index in input list for heat capacity at second wall layer, ground floor level
327	INTEGER(iwp) :: ind_hc3_agfl = 8 !< index in input list for heat capacity at third wall layer, above ground floor level
328	INTEGER(iwp) :: ind_hc3_gfl = 28 !< index in input list for heat capacity at third wall layer, ground floor level
329	INTEGER(iwp) :: ind_gflh = 20 !< index in input list for ground floor level height
330	INTEGER(iwp) :: ind_lai_r_agfl = 4 !< index in input list for LAI on roof, above ground floor level
331	INTEGER(iwp) :: ind_lai_r_gfl = 4 !< index in input list for LAI on roof, ground floor level
332	INTEGER(iwp) :: ind_lai_w_agfl = 5 !< index in input list for LAI on wall, above ground floor level
333	INTEGER(iwp) :: ind_lai_w_gfl = 25 !< index in input list for LAI on wall, ground floor level
334	INTEGER(iwp) :: ind_tc1_agfl = 9 !< index in input list for thermal conductivity at first wall layer, above ground floor level
335	INTEGER(iwp) :: ind_tc1_gfl = 29 !< index in input list for thermal conductivity at first wall layer, ground floor level
336	INTEGER(iwp) :: ind_tc2_agfl = 10 !< index in input list for thermal conductivity at second wall layer, above ground floor level
337	INTEGER(iwp) :: ind_tc2_gfl = 30 !< index in input list for thermal conductivity at second wall layer, ground floor level
338	INTEGER(iwp) :: ind_tc3_agfl = 11 !< index in input list for thermal conductivity at third wall layer, above ground floor level
339	INTEGER(iwp) :: ind_tc3_gfl = 31 !< index in input list for thermal conductivity at third wall layer, ground floor level
340	INTEGER(iwp) :: ind_thick_1 = 41 !< index for wall layer thickness - 1st layer
341	INTEGER(iwp) :: ind_thick_2 = 42 !< index for wall layer thickness - 2nd layer
342	INTEGER(iwp) :: ind_thick_3 = 43 !< index for wall layer thickness - 3rd layer
343	INTEGER(iwp) :: ind_thick_4 = 44 !< index for wall layer thickness - 4th layer
344	INTEGER(iwp) :: ind_trans_agfl = 17 !< index in input list for window transmissivity, above ground floor level
345	INTEGER(iwp) :: ind_trans_gfl = 35 !< index in input list for window transmissivity, ground floor level
346	INTEGER(iwp) :: ind_wall_frac_agfl = 0 !< index in input list for wall fraction, above ground floor level
347	INTEGER(iwp) :: ind_wall_frac_gfl = 21 !< index in input list for wall fraction, ground floor level
348	INTEGER(iwp) :: ind_win_frac_agfl = 1 !< index in input list for window fraction, above ground floor level
349	INTEGER(iwp) :: ind_win_frac_gfl = 22 !< index in input list for window fraction, ground floor level
350	INTEGER(iwp) :: ind_z0_agfl = 18 !< index in input list for z0, above ground floor level
351	INTEGER(iwp) :: ind_z0_gfl = 36 !< index in input list for z0, ground floor level
352	INTEGER(iwp) :: ind_z0qh_agfl = 19 !< index in input list for z0h / z0q, above ground floor level
353	INTEGER(iwp) :: ind_z0qh_gfl = 37 !< index in input list for z0h / z0q, ground floor level
354
355
356	REAL(wp) :: roof_height_limit = 4._wp !< height for distinguish between land surfaces and roofs
357	REAL(wp) :: ground_floor_level = 4.0_wp !< default ground floor level
358	REAL(wp) :: ra_horiz_coef = 5.0_wp !< mysterious coefficient for correction of overestimation
359	!< of r_a for horizontal surfaces -> TODO
360
361	CHARACTER(37), DIMENSION(0:6), PARAMETER :: building_type_name = (/ &
362	'user-defined ', & ! 0
363	'residential - 1950 ', & ! 1
364	'residential 1951 - 2000 ', & ! 2
365	'residential 2001 - ', & ! 3
366	'office - 1950 ', & ! 4
367	'office 1951 - 2000 ', & ! 5
368	'office 2001 - ' & ! 6
369	/)
370	!
371	!-- building parameters, 4 different types
372	!-- 0 - wall fraction, 1- window fraction, 2 - green fraction on wall, 3- green fraction
373	!-- at roof, 4 - lai of green fraction at roof, 5 - lai of green fraction at wall,
374	!-- 6 - heat capacity of wall layer 1, 7 - heat capacity of wall layer 2,
375	!-- 8 - heat capacity of wall layer 3, 9 - thermal conductivity of wall layer 1,
376	!-- 10 - thermal conductivity of wall layer 2, 11 - thermal conductivity of wall layer 3,
377	!-- 12 - indoor target summer temperature ( K ), 13 - indoor target winter temperature (K),
378	!-- 14 - emissivity of wall fraction, 15 - emissivity of green fraction, 16 - emissivity of window fraction,
379	!-- 17 - transmissivity of window fraction, 18 - z0, 19 - z0h/z0q, 20 - ground floor height,
380	!-- 21 - ground floor wall fraction, 22 - ground floor window fraction, 23 ground floor green fraction,
381	!-- 24 - ground floor green fraction on roof, 25 - ground floor lai of green fraction,
382	!-- 26 - ground floor heat capacity of wall layer 1, 27 - ground floor heat capacity of wall layer 1,
383	!-- 28 - ground floor heat capacity of wall layer 3, 29 - ground floor thermal conductivity of wall layer 1,
384	!-- 30 - ground floor thermal conductivity of wall layer 2, 31 - ground floor thermal conductivity of wall layer 3,
385	!-- 32 - ground floor emissivity of wall fraction, 33 - ground floor emissivity of green fraction,
386	!-- 34 - ground floor emissivity of window fraction, 35 - ground floor transmissivity of window fraction,
387	!-- 36 - ground floor z0, 37 - ground floor z0h/z0q, 38 - albedo type wall fraction
388	!-- 39 - albedo type green fraction, 40 - albedo type window fraction
389	!-- 41 - wall layer thickness - 1st layer, 42 - wall layer thickness - 2nd layer,
390	!-- 43 - wall layer thickness - 3rd layer, 44 - wall layer thickness - 4th layer,
391	!-- 45 - heat capacity of the wall surface, 46 - heat conductivity
392	!-- Please note, only preleminary dummy values so far!
393	REAL(wp), DIMENSION(0:46,1:6), PARAMETER :: building_pars = RESHAPE( (/ &
394	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5
395	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11
396	296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17
397	0.001_wp, 0.0001_wp, 4.0_wp, & !parameter 18-20
398	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25
399	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28
400	0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31
401	0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37
402	24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40
403	0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44
404	20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 1
405	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5
406	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11
407	296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17
408	0.001_wp, 0.0001_wp, 4.0_wp, & !parameter 18-20
409	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25
410	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28
411	0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31
412	0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37
413	24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40
414	0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44
415	20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 2
416	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5
417	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11
418	296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17
419	0.001_wp, 0.0001_wp, 4.0_wp, & !parameter 18-20
420	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25
421	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28
422	0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31
423	0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37
424	24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40
425	0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44
426	20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 3
427	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5
428	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11
429	296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17
430	0.01_wp, 0.001_wp, 4.0_wp, & !parameter 18-20
431	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25
432	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28
433	0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31
434	0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37
435	24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40
436	0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44
437	20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 4
438	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5
439	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11
440	296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17
441	0.001_wp, 0.0001_wp, 4.0_wp, & !parameter 18-20
442	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25
443	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28
444	0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31
445	0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37
446	24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40
447	0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44
448	20000.0_wp, 10.0_wp, & !parameter 45-46 - end of type 5
449	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, & !parameter 0-5
450	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp, & !parameter 6-11
451	296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp, & !parameter 12-17
452	0.001_wp, 0.0001_wp, 4.0_wp, & !parameter 18-20
453	1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, & !parameter 21-25
454	1000000.0_wp, 1000000.0_wp, 1000000.0_wp, & !parameter 26-28
455	0.3_wp, 0.3_wp, 0.3_wp, & !parameter 29-31
456	0.4_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp, & !parameter 32-37
457	24.0_wp, 24.0_wp, 24.0_wp, & !parameter 38-40
458	0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp, & !parameter 41-44
459	20000.0_wp, 10.0_wp & !parameter 45-46 - end of type 6
460	/), &
461	(/47, 6/) )
462
463	!
464	!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
465	TYPE t_surf_vertical
466	REAL(wp), DIMENSION(:), ALLOCATABLE :: t
467	END TYPE t_surf_vertical
468	!
469	!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
470	TYPE t_wall_vertical
471	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t
472	END TYPE t_wall_vertical
473
474
475	!-- arrays for time averages
476	!-- Attention: the variable rad_net_av is also used in the 3d field variable in radiation_model_mod.f90. It may be better to rename it
477	REAL(wp), DIMENSION(:), ALLOCATABLE :: rad_net_av !< average of rad_net_l
478	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
479	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
480	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
481	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
482	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
483	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
484	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
485	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
486	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
487	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
488	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
489	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfhf_av !< average of total radiation flux incoming to minus outgoing from local surface
490	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_av !< average of wghf_eb
491	REAL(wp), DIMENSION(:), ALLOCATABLE :: wshf_eb_av !< average of wshf_eb
492	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_wall_av !< Average of t_wall
493	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_green_av !< average of wghf_eb_green
494	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_green_av !< Average of t_green
495	REAL(wp), DIMENSION(:), ALLOCATABLE :: wghf_eb_window_av !< average of wghf_eb_window
496	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t_window_av !< Average of t_window
497
498
499	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
500	!-- anthropogenic heat sources
501	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
502	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2)
503	REAL(wp), DIMENSION(:), ALLOCATABLE :: aheatprof !< diurnal profile of anthropogenic heat
504
505	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
506	!-- wall surface model
507	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
508	!-- wall surface model constants
509	INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched)
510	INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched)
511	INTEGER(iwp), PARAMETER :: nzw = 4 !< number of wall layers (fixed for now)
512
513	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
514	!< normalized soil, wall and roof layer depths (m/m)
515	! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.33_wp, 0.66_wp, 1.0_wp /)
516	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
517	! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.33_wp, 0.66_wp, 1.0_wp /)
518	! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
519	!< normalized window layer depths (m/m)
520	! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
521	!< normalized green layer depths (m/m)
522	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
523	! REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /)
524
525
526	REAL(wp) :: wall_inner_temperature = 296.0_wp !< temperature of the inner wall surface (~23 degrees C) (K)
527	REAL(wp) :: roof_inner_temperature = 296.0_wp !< temperature of the inner roof surface (~23 degrees C) (K)
528	REAL(wp) :: soil_inner_temperature = 283.0_wp !< temperature of the deep soil (~10 degrees C) (K)
529	REAL(wp) :: window_inner_temperature = 296.0_wp !< temperature of the inner window surface (~23 degrees C) (K)
530
531	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
532	!-- surface and material model variables for walls, ground, roofs
533	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
534	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m)
535	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_window !< normalized window layer depths (m)
536	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_green !< normalized green layer depths (m)
537
538	#if defined( __nopointer )
539	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h !< wall surface temperature (K) at horizontal walls
540	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_p !< progn. wall surface temperature (K) at horizontal walls
541	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h !< window surface temperature (K) at horizontal walls
542	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_p !< progn. window surface temperature (K) at horizontal walls
543	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h !< green surface temperature (K) at horizontal walls
544	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_p !< progn. green surface temperature (K) at horizontal walls
545	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h !< near surface temperature (10cm) (K) at horizontal walls
546	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h_p !< progn. near surface temperature (10cm) (K) at horizontal walls
547	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v
548	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_p
549	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v
550	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_p
551	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v
552	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_p
553	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v
554	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v_p
555	#else
556	REAL(wp), DIMENSION(:), POINTER :: t_surf_h
557	REAL(wp), DIMENSION(:), POINTER :: t_surf_h_p
558	REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h
559	REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h_p
560	REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h
561	REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h_p
562	REAL(wp), DIMENSION(:), POINTER :: t_surf_10cm_h
563	REAL(wp), DIMENSION(:), POINTER :: t_surf_10cm_h_p
564
565	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_1
566	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_h_2
567	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_1
568	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_2
569	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_1
570	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_2
571	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h_1
572	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_h_2
573
574	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v
575	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_v_p
576	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v
577	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v_p
578	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v
579	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v_p
580	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_10cm_v
581	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_10cm_v_p
582
583	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_1
584	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_v_2
585	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_1
586	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_2
587	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_1
588	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_2
589	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v_1
590	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_10cm_v_2
591
592	#endif
593	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_av !< average of wall surface temperature (K)
594	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_av !< average of window surface temperature (K)
595	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_av !< average of green wall surface temperature (K)
596	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_10cm_av !< average of whole wall surface temperature (K)
597
598	!-- Temporal tendencies for time stepping
599	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_m !< surface temperature tendency of wall (K)
600	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_window_m !< surface temperature tendency of window (K)
601	REAL(wp), DIMENSION(:), ALLOCATABLE :: tt_surface_green_m !< surface temperature tendency of green wall (K)
602
603	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
604	!-- Energy balance variables
605	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
606	!-- parameters of the land, roof and wall surfaces
607
608	#if defined( __nopointer )
609	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h !< Wall temperature (K)
610	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av !< Average of t_wall
611	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_p !< Prog. wall temperature (K)
612	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h !< Window temperature (K)
613	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_av !< Average of t_window
614	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_p !< Prog. window temperature (K)
615	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h !< Green temperature (K)
616	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_av !< Average of t_green
617	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_p !< Prog. green temperature (K)
618
619	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v !< Wall temperature (K)
620	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av !< Average of t_wall
621	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_p !< Prog. wall temperature (K)
622	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v !< Window temperature (K)
623	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_av !< Average of t_window
624	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_p !< Prog. window temperature (K)
625	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v !< Green temperature (K)
626	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_av !< Average of t_green
627	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_p !< Prog. green temperature (K)
628	#else
629	REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p
630	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_av, t_wall_h_1, t_wall_h_2
631	REAL(wp), DIMENSION(:,:), POINTER :: t_window_h, t_window_h_p
632	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_av, t_window_h_1, t_window_h_2
633	REAL(wp), DIMENSION(:,:), POINTER :: t_green_h, t_green_h_p
634	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_av, t_green_h_1, t_green_h_2
635
636	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p
637	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_av, t_wall_v_1, t_wall_v_2
638	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_window_v, t_window_v_p
639	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_av, t_window_v_1, t_window_v_2
640	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_green_v, t_green_v_p
641	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_av, t_green_v_1, t_green_v_2
642	#endif
643
644	!-- Wall temporal tendencies for time stepping
645	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_wall_m !< t_wall prognostic array
646	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_window_m !< t_window prognostic array
647	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_green_m !< t_green prognostic array
648
649	!-- Surface and material parameters classes (surface_type)
650	!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
651	INTEGER(iwp) :: n_surface_types !< number of the wall type categories
652	INTEGER(iwp), PARAMETER :: n_surface_params = 8 !< number of parameters for each type of the wall
653	INTEGER(iwp), PARAMETER :: ialbedo = 1 !< albedo of the surface
654	INTEGER(iwp), PARAMETER :: iemiss = 2 !< emissivity of the surface
655	INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity Î»S between air and surface ( W mâ2 Kâ1 )
656	INTEGER(iwp), PARAMETER :: irough = 4 !< roughness relative to concrete
657	INTEGER(iwp), PARAMETER :: icsurf = 5 !< Surface skin layer heat capacity (J mâ2 Kâ1 )
658	INTEGER(iwp), PARAMETER :: ithick = 6 !< thickness of the surface (wall, roof, land) ( m )
659	INTEGER(iwp), PARAMETER :: irhoC = 7 !< volumetric heat capacity rho*C of the material ( J mâ3 Kâ1 )
660	INTEGER(iwp), PARAMETER :: ilambdah = 8 !< thermal conductivity Î»H of the wall (W mâ1 Kâ1 )
661	CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports)
662	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types
663	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: surface_params !< parameters of wall types
664
665
666	!-- interfaces of subroutines accessed from outside of this module
667	INTERFACE usm_boundary_condition
668	MODULE PROCEDURE usm_boundary_condition
669	END INTERFACE usm_boundary_condition
670
671	INTERFACE usm_check_data_output
672	MODULE PROCEDURE usm_check_data_output
673	END INTERFACE usm_check_data_output
674
675	INTERFACE usm_check_parameters
676	MODULE PROCEDURE usm_check_parameters
677	END INTERFACE usm_check_parameters
678
679	INTERFACE usm_data_output_3d
680	MODULE PROCEDURE usm_data_output_3d
681	END INTERFACE usm_data_output_3d
682
683	INTERFACE usm_define_netcdf_grid
684	MODULE PROCEDURE usm_define_netcdf_grid
685	END INTERFACE usm_define_netcdf_grid
686
687	INTERFACE usm_init_urban_surface
688	MODULE PROCEDURE usm_init_urban_surface
689	END INTERFACE usm_init_urban_surface
690
691	INTERFACE usm_material_heat_model
692	MODULE PROCEDURE usm_material_heat_model
693	END INTERFACE usm_material_heat_model
694
695	INTERFACE usm_green_heat_model
696	MODULE PROCEDURE usm_green_heat_model
697	END INTERFACE usm_green_heat_model
698
699	INTERFACE usm_parin
700	MODULE PROCEDURE usm_parin
701	END INTERFACE usm_parin
702
703	INTERFACE usm_temperature_near_surface
704	MODULE PROCEDURE usm_temperature_near_surface
705	END INTERFACE usm_temperature_near_surface
706
707	INTERFACE usm_read_restart_data
708	MODULE PROCEDURE usm_read_restart_data
709	END INTERFACE usm_read_restart_data
710
711	INTERFACE usm_surface_energy_balance
712	MODULE PROCEDURE usm_surface_energy_balance
713	END INTERFACE usm_surface_energy_balance
714
715	INTERFACE usm_swap_timelevel
716	MODULE PROCEDURE usm_swap_timelevel
717	END INTERFACE usm_swap_timelevel
718
719	INTERFACE usm_write_restart_data
720	MODULE PROCEDURE usm_write_restart_data
721	END INTERFACE usm_write_restart_data
722
723	INTERFACE usm_allocate_surface
724	MODULE PROCEDURE usm_allocate_surface
725	END INTERFACE usm_allocate_surface
726
727	INTERFACE usm_average_3d_data
728	MODULE PROCEDURE usm_average_3d_data
729	END INTERFACE usm_average_3d_data
730
731
732	SAVE
733
734	PRIVATE
735
736	!-- Public functions
737	PUBLIC usm_boundary_condition, usm_check_parameters, usm_init_urban_surface,&
738	usm_read_restart_data, &
739	usm_surface_energy_balance, usm_material_heat_model, &
740	usm_swap_timelevel, usm_check_data_output, usm_average_3d_data, &
741	usm_data_output_3d, usm_define_netcdf_grid, usm_parin, &
742	usm_write_restart_data, usm_allocate_surface
743
744	!-- Public parameters, constants and initial values
745	PUBLIC usm_anthropogenic_heat, usm_material_model, ra_horiz_coef, &
746	usm_green_heat_model, usm_temperature_near_surface
747
748
749
750	CONTAINS
751
752	!------------------------------------------------------------------------------!
753	! Description:
754	! ------------
755	!> This subroutine creates the necessary indices of the urban surfaces
756	!> and plant canopy and it allocates the needed arrays for USM
757	!------------------------------------------------------------------------------!
758	SUBROUTINE usm_allocate_surface
759
760	IMPLICIT NONE
761
762	INTEGER(iwp) :: l
763
764	!
765	!-- Allocate radiation arrays which are part of the new data type.
766	!-- For horizontal surfaces.
767	ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns) )
768	ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
769	!
770	!-- For vertical surfaces
771	DO l = 0, 3
772	ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) )
773	ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
774	ENDDO
775
776	!-- Wall surface model
777	!-- allocate arrays for wall surface model and define pointers
778
779	!-- allocate array of wall types and wall parameters
780	ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
781	DO l = 0, 3
782	ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
783	ENDDO
784	!
785	!-- Allocate albedo_type and albedo. Each surface element
786	!-- has 3 values, 0: wall fraction, 1: green fraction, 2: window fraction.
787	ALLOCATE( surf_usm_h%albedo_type(0:2,1:surf_usm_h%ns) )
788	ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) )
789	surf_usm_h%albedo_type = albedo_type
790	DO l = 0, 3
791	ALLOCATE( surf_usm_v(l)%albedo_type(0:2,1:surf_usm_v(l)%ns) )
792	ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
793	surf_usm_v(l)%albedo_type = albedo_type
794	ENDDO
795
796
797	!
798	!-- Allocate indoor target temperature for summer and winter
799	ALLOCATE( surf_usm_h%target_temp_summer(1:surf_usm_h%ns) )
800	ALLOCATE( surf_usm_h%target_temp_winter(1:surf_usm_h%ns) )
801	DO l = 0, 3
802	ALLOCATE( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) )
803	ALLOCATE( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) )
804	ENDDO
805	!
806	!-- Allocate flag indicating ground floor level surface elements
807	ALLOCATE ( surf_usm_h%ground_level(1:surf_usm_h%ns) )
808	DO l = 0, 3
809	ALLOCATE( surf_usm_v(l)%ground_level(1:surf_usm_v(l)%ns) )
810	ENDDO
811	!
812	!-- Allocate arrays for relative surface fraction.
813	!-- 0 - wall fraction, 1 - green fraction, 2 - window fraction
814	ALLOCATE( surf_usm_h%frac(0:2,1:surf_usm_h%ns) )
815	surf_usm_h%frac = 0.0_wp
816	DO l = 0, 3
817	ALLOCATE( surf_usm_v(l)%frac(0:2,1:surf_usm_v(l)%ns) )
818	surf_usm_v(l)%frac = 0.0_wp
819	ENDDO
820
821	!-- wall and roof surface parameters. First for horizontal surfaces
822	ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) )
823	ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) )
824	ALLOCATE ( surf_usm_h%lambda_surf_window(1:surf_usm_h%ns) )
825	ALLOCATE ( surf_usm_h%lambda_surf_green(1:surf_usm_h%ns) )
826	ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) )
827	ALLOCATE ( surf_usm_h%c_surface_window(1:surf_usm_h%ns) )
828	ALLOCATE ( surf_usm_h%c_surface_green(1:surf_usm_h%ns) )
829	ALLOCATE ( surf_usm_h%transmissivity(1:surf_usm_h%ns) )
830	ALLOCATE ( surf_usm_h%lai(1:surf_usm_h%ns) )
831	ALLOCATE ( surf_usm_h%emissivity(0:2,1:surf_usm_h%ns) )
832	ALLOCATE ( surf_usm_h%r_a(1:surf_usm_h%ns) )
833	ALLOCATE ( surf_usm_h%r_a_green(1:surf_usm_h%ns) )
834	ALLOCATE ( surf_usm_h%r_a_window(1:surf_usm_h%ns) )
835
836	!
837	!-- For vertical surfaces.
838	DO l = 0, 3
839	ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) )
840	ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) )
841	ALLOCATE ( surf_usm_v(l)%lambda_surf_window(1:surf_usm_v(l)%ns) )
842	ALLOCATE ( surf_usm_v(l)%c_surface_window(1:surf_usm_v(l)%ns) )
843	ALLOCATE ( surf_usm_v(l)%lambda_surf_green(1:surf_usm_v(l)%ns) )
844	ALLOCATE ( surf_usm_v(l)%c_surface_green(1:surf_usm_v(l)%ns) )
845	ALLOCATE ( surf_usm_v(l)%transmissivity(1:surf_usm_v(l)%ns) )
846	ALLOCATE ( surf_usm_v(l)%lai(1:surf_usm_v(l)%ns) )
847	ALLOCATE ( surf_usm_v(l)%emissivity(0:2,1:surf_usm_v(l)%ns) )
848	ALLOCATE ( surf_usm_v(l)%r_a(1:surf_usm_v(l)%ns) )
849	ALLOCATE ( surf_usm_v(l)%r_a_green(1:surf_usm_v(l)%ns) )
850	ALLOCATE ( surf_usm_v(l)%r_a_window(1:surf_usm_v(l)%ns) )
851	ENDDO
852
853	!
854	!-- allocate wall and roof material parameters. First for horizontal surfaces
855	ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) )
856	ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns) )
857	ALLOCATE ( surf_usm_h%thickness_green(1:surf_usm_h%ns) )
858	ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
859	ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
860	ALLOCATE ( surf_usm_h%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
861	ALLOCATE ( surf_usm_h%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
862	ALLOCATE ( surf_usm_h%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
863	ALLOCATE ( surf_usm_h%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
864
865	!
866	!-- For vertical surfaces.
867	DO l = 0, 3
868	ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) )
869	ALLOCATE ( surf_usm_v(l)%thickness_window(1:surf_usm_v(l)%ns) )
870	ALLOCATE ( surf_usm_v(l)%thickness_green(1:surf_usm_v(l)%ns) )
871	ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
872	ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
873	ALLOCATE ( surf_usm_v(l)%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
874	ALLOCATE ( surf_usm_v(l)%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
875	ALLOCATE ( surf_usm_v(l)%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
876	ALLOCATE ( surf_usm_v(l)%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
877	ENDDO
878
879	!-- allocate wall and roof layers sizes. For horizontal surfaces.
880	ALLOCATE ( zwn(nzb_wall:nzt_wall) )
881	ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
882	ALLOCATE ( zwn_window(nzb_wall:nzt_wall) )
883	ALLOCATE ( surf_usm_h%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
884	ALLOCATE ( zwn_green(nzb_wall:nzt_wall) )
885	ALLOCATE ( surf_usm_h%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
886	ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
887	ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
888	ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
889	ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
890	ALLOCATE ( surf_usm_h%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
891	ALLOCATE ( surf_usm_h%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
892	ALLOCATE ( surf_usm_h%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
893	ALLOCATE ( surf_usm_h%zw_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
894	ALLOCATE ( surf_usm_h%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
895	ALLOCATE ( surf_usm_h%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
896	ALLOCATE ( surf_usm_h%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
897	ALLOCATE ( surf_usm_h%zw_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
898	!
899	!-- For vertical surfaces.
900	DO l = 0, 3
901	ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
902	ALLOCATE ( surf_usm_v(l)%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
903	ALLOCATE ( surf_usm_v(l)%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
904	ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
905	ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
906	ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
907	ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
908	ALLOCATE ( surf_usm_v(l)%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
909	ALLOCATE ( surf_usm_v(l)%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
910	ALLOCATE ( surf_usm_v(l)%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
911	ALLOCATE ( surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
912	ALLOCATE ( surf_usm_v(l)%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
913	ALLOCATE ( surf_usm_v(l)%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
914	ALLOCATE ( surf_usm_v(l)%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
915	ALLOCATE ( surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
916	ENDDO
917
918	!-- allocate wall and roof temperature arrays, for horizontal walls
919	#if defined( __nopointer )
920	IF ( .NOT. ALLOCATED( t_surf_h ) ) &
921	ALLOCATE ( t_surf_h(1:surf_usm_h%ns) )
922	IF ( .NOT. ALLOCATED( t_surf_h_p ) ) &
923	ALLOCATE ( t_surf_h_p(1:surf_usm_h%ns) )
924	IF ( .NOT. ALLOCATED( t_wall_h ) ) &
925	ALLOCATE ( t_wall_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
926	IF ( .NOT. ALLOCATED( t_wall_h_p ) ) &
927	ALLOCATE ( t_wall_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
928	IF ( .NOT. ALLOCATED( t_surf_window_h ) ) &
929	ALLOCATE ( t_surf_window_h(1:surf_usm_h%ns) )
930	IF ( .NOT. ALLOCATED( t_surf_window_h_p ) ) &
931	ALLOCATE ( t_surf_window_h_p(1:surf_usm_h%ns) )
932	IF ( .NOT. ALLOCATED( t_window_h ) ) &
933	ALLOCATE ( t_window_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
934	IF ( .NOT. ALLOCATED( t_window_h_p ) ) &
935	ALLOCATE ( t_window_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
936	IF ( .NOT. ALLOCATED( t_surf_green_h ) ) &
937	ALLOCATE ( t_surf_green_h(1:surf_usm_h%ns) )
938	IF ( .NOT. ALLOCATED( t_surf_green_h_p ) ) &
939	ALLOCATE ( t_surf_green_h_p(1:surf_usm_h%ns) )
940	IF ( .NOT. ALLOCATED( t_green_h ) ) &
941	ALLOCATE ( t_green_h(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
942	IF ( .NOT. ALLOCATED( t_green_h_p ) ) &
943	ALLOCATE ( t_green_h_p(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
944	IF ( .NOT. ALLOCATED( t_surf_10cm_h ) ) &
945	ALLOCATE ( t_surf_10cm_h(1:surf_usm_h%ns) )
946	IF ( .NOT. ALLOCATED( t_surf_10cm_h_p ) ) &
947	ALLOCATE ( t_surf_10cm_h_p(1:surf_usm_h%ns) )
948	#else
949	!
950	!-- Allocate if required. Note, in case of restarts, some of these arrays
951	!-- might be already allocated.
952	IF ( .NOT. ALLOCATED( t_surf_h_1 ) ) &
953	ALLOCATE ( t_surf_h_1(1:surf_usm_h%ns) )
954	IF ( .NOT. ALLOCATED( t_surf_h_2 ) ) &
955	ALLOCATE ( t_surf_h_2(1:surf_usm_h%ns) )
956	IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) &
957	ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
958	IF ( .NOT. ALLOCATED( t_wall_h_2 ) ) &
959	ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
960	IF ( .NOT. ALLOCATED( t_surf_window_h_1 ) ) &
961	ALLOCATE ( t_surf_window_h_1(1:surf_usm_h%ns) )
962	IF ( .NOT. ALLOCATED( t_surf_window_h_2 ) ) &
963	ALLOCATE ( t_surf_window_h_2(1:surf_usm_h%ns) )
964	IF ( .NOT. ALLOCATED( t_window_h_1 ) ) &
965	ALLOCATE ( t_window_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
966	IF ( .NOT. ALLOCATED( t_window_h_2 ) ) &
967	ALLOCATE ( t_window_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
968	IF ( .NOT. ALLOCATED( t_surf_green_h_1 ) ) &
969	ALLOCATE ( t_surf_green_h_1(1:surf_usm_h%ns) )
970	IF ( .NOT. ALLOCATED( t_surf_green_h_2 ) ) &
971	ALLOCATE ( t_surf_green_h_2(1:surf_usm_h%ns) )
972	IF ( .NOT. ALLOCATED( t_green_h_1 ) ) &
973	ALLOCATE ( t_green_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
974	IF ( .NOT. ALLOCATED( t_green_h_2 ) ) &
975	ALLOCATE ( t_green_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
976	IF ( .NOT. ALLOCATED( t_surf_10cm_h_1 ) ) &
977	ALLOCATE ( t_surf_10cm_h_1(1:surf_usm_h%ns) )
978	IF ( .NOT. ALLOCATED( t_surf_10cm_h_2 ) ) &
979	ALLOCATE ( t_surf_10cm_h_2(1:surf_usm_h%ns) )
980	!
981	!-- initial assignment of the pointers
982	t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2
983	t_window_h => t_window_h_1; t_window_h_p => t_window_h_2
984	t_green_h => t_green_h_1; t_green_h_p => t_green_h_2
985	t_surf_h => t_surf_h_1; t_surf_h_p => t_surf_h_2
986	t_surf_window_h => t_surf_window_h_1; t_surf_window_h_p => t_surf_window_h_2
987	t_surf_green_h => t_surf_green_h_1; t_surf_green_h_p => t_surf_green_h_2
988	t_surf_10cm_h => t_surf_10cm_h_1; t_surf_10cm_h_p => t_surf_10cm_h_2
989
990	#endif
991
992	!-- allocate wall and roof temperature arrays, for vertical walls if required
993	#if defined( __nopointer )
994	DO l = 0, 3
995	IF ( .NOT. ALLOCATED( t_surf_v(l)%t ) ) &
996	ALLOCATE ( t_surf_v(l)%t(1:surf_usm_v(l)%ns) )
997	IF ( .NOT. ALLOCATED( t_surf_v_p(l)%t ) ) &
998	ALLOCATE ( t_surf_v_p(l)%t(1:surf_usm_v(l)%ns) )
999	IF ( .NOT. ALLOCATED( t_wall_v(l)%t ) ) &
1000	ALLOCATE ( t_wall_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1001	IF ( .NOT. ALLOCATED( t_wall_v_p(l)%t ) ) &
1002	ALLOCATE ( t_wall_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1003	IF ( .NOT. ALLOCATED( t_surf_window_v(l)%t ) ) &
1004	ALLOCATE ( t_surf_window_v(l)%t(1:surf_usm_v(l)%ns) )
1005	IF ( .NOT. ALLOCATED( t_surf_window_v_p(l)%t ) ) &
1006	ALLOCATE ( t_surf_window_v_p(l)%t(1:surf_usm_v(l)%ns) )
1007	IF ( .NOT. ALLOCATED( t_window_v(l)%t ) ) &
1008	ALLOCATE ( t_window_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1009	IF ( .NOT. ALLOCATED( t_window_v_p(l)%t ) ) &
1010	ALLOCATE ( t_window_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1011	IF ( .NOT. ALLOCATED( t_green_v(l)%t ) ) &
1012	ALLOCATE ( t_green_v(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1013	IF ( .NOT. ALLOCATED( t_green_v_p(l)%t ) ) &
1014	ALLOCATE ( t_green_v_p(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1015	IF ( .NOT. ALLOCATED( t_surf_green_v(l)%t ) ) &
1016	ALLOCATE ( t_surf_green_v(l)%t(1:surf_usm_v(l)%ns) )
1017	IF ( .NOT. ALLOCATED( t_surf_green_v_p(l)%t ) ) &
1018	ALLOCATE ( t_surf_green_v_p(l)%t(1:surf_usm_v(l)%ns) )
1019	IF ( .NOT. ALLOCATED( t_surf_10cm_v(l)%t ) ) &
1020	ALLOCATE ( t_surf_10cm_v(l)%t(1:surf_usm_v(l)%ns) )
1021	IF ( .NOT. ALLOCATED( t_surf_10cm_v_p(l)%t ) ) &
1022	ALLOCATE ( t_surf_10cm_v_p(l)%t(1:surf_usm_v(l)%ns) )
1023	ENDDO
1024	#else
1025	!
1026	!-- Allocate if required. Note, in case of restarts, some of these arrays
1027	!-- might be already allocated.
1028	DO l = 0, 3
1029	IF ( .NOT. ALLOCATED( t_surf_v_1(l)%t ) ) &
1030	ALLOCATE ( t_surf_v_1(l)%t(1:surf_usm_v(l)%ns) )
1031	IF ( .NOT. ALLOCATED( t_surf_v_2(l)%t ) ) &
1032	ALLOCATE ( t_surf_v_2(l)%t(1:surf_usm_v(l)%ns) )
1033	IF ( .NOT. ALLOCATED( t_wall_v_1(l)%t ) ) &
1034	ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1035	IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) ) &
1036	ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1037	IF ( .NOT. ALLOCATED( t_surf_window_v_1(l)%t ) ) &
1038	ALLOCATE ( t_surf_window_v_1(l)%t(1:surf_usm_v(l)%ns) )
1039	IF ( .NOT. ALLOCATED( t_surf_window_v_2(l)%t ) ) &
1040	ALLOCATE ( t_surf_window_v_2(l)%t(1:surf_usm_v(l)%ns) )
1041	IF ( .NOT. ALLOCATED( t_window_v_1(l)%t ) ) &
1042	ALLOCATE ( t_window_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1043	IF ( .NOT. ALLOCATED( t_window_v_2(l)%t ) ) &
1044	ALLOCATE ( t_window_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1045	IF ( .NOT. ALLOCATED( t_surf_green_v_1(l)%t ) ) &
1046	ALLOCATE ( t_surf_green_v_1(l)%t(1:surf_usm_v(l)%ns) )
1047	IF ( .NOT. ALLOCATED( t_surf_green_v_2(l)%t ) ) &
1048	ALLOCATE ( t_surf_green_v_2(l)%t(1:surf_usm_v(l)%ns) )
1049	IF ( .NOT. ALLOCATED( t_green_v_1(l)%t ) ) &
1050	ALLOCATE ( t_green_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1051	IF ( .NOT. ALLOCATED( t_green_v_2(l)%t ) ) &
1052	ALLOCATE ( t_green_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1053	IF ( .NOT. ALLOCATED( t_surf_10cm_v_1(l)%t ) ) &
1054	ALLOCATE ( t_surf_10cm_v_1(l)%t(1:surf_usm_v(l)%ns) )
1055	IF ( .NOT. ALLOCATED( t_surf_10cm_v_2(l)%t ) ) &
1056	ALLOCATE ( t_surf_10cm_v_2(l)%t(1:surf_usm_v(l)%ns) )
1057	ENDDO
1058	!
1059	!-- initial assignment of the pointers
1060	t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2
1061	t_surf_v => t_surf_v_1; t_surf_v_p => t_surf_v_2
1062	t_window_v => t_window_v_1; t_window_v_p => t_window_v_2
1063	t_green_v => t_green_v_1; t_green_v_p => t_green_v_2
1064	t_surf_window_v => t_surf_window_v_1; t_surf_window_v_p => t_surf_window_v_2
1065	t_surf_green_v => t_surf_green_v_1; t_surf_green_v_p => t_surf_green_v_2
1066	t_surf_10cm_v => t_surf_10cm_v_1; t_surf_10cm_v_p => t_surf_10cm_v_2
1067
1068	#endif
1069	!
1070	!-- Allocate intermediate timestep arrays. For horizontal surfaces.
1071	ALLOCATE ( surf_usm_h%tt_surface_m(1:surf_usm_h%ns) )
1072	ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1073	ALLOCATE ( surf_usm_h%tt_surface_window_m(1:surf_usm_h%ns) )
1074	ALLOCATE ( surf_usm_h%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1075	ALLOCATE ( surf_usm_h%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1076	ALLOCATE ( surf_usm_h%tt_surface_green_m(1:surf_usm_h%ns) )
1077
1078	!
1079	!-- Set inital values for prognostic quantities
1080	IF ( ALLOCATED( surf_usm_h%tt_surface_m ) ) surf_usm_h%tt_surface_m = 0.0_wp
1081	IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp
1082	IF ( ALLOCATED( surf_usm_h%tt_surface_window_m ) ) surf_usm_h%tt_surface_window_m = 0.0_wp
1083	IF ( ALLOCATED( surf_usm_h%tt_window_m ) ) surf_usm_h%tt_window_m = 0.0_wp
1084	IF ( ALLOCATED( surf_usm_h%tt_green_m ) ) surf_usm_h%tt_green_m = 0.0_wp
1085	IF ( ALLOCATED( surf_usm_h%tt_surface_green_m ) ) surf_usm_h%tt_surface_green_m = 0.0_wp
1086	!
1087	!-- Now, for vertical surfaces
1088	DO l = 0, 3
1089	ALLOCATE ( surf_usm_v(l)%tt_surface_m(1:surf_usm_v(l)%ns) )
1090	ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1091	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_m ) ) surf_usm_v(l)%tt_surface_m = 0.0_wp
1092	IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp
1093	ALLOCATE ( surf_usm_v(l)%tt_surface_window_m(1:surf_usm_v(l)%ns) )
1094	ALLOCATE ( surf_usm_v(l)%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1095	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_window_m ) ) surf_usm_v(l)%tt_surface_window_m = 0.0_wp
1096	IF ( ALLOCATED( surf_usm_v(l)%tt_window_m ) ) surf_usm_v(l)%tt_window_m = 0.0_wp
1097	ALLOCATE ( surf_usm_v(l)%tt_surface_green_m(1:surf_usm_v(l)%ns) )
1098	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_green_m ) ) surf_usm_v(l)%tt_surface_green_m = 0.0_wp
1099	ALLOCATE ( surf_usm_v(l)%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1100	IF ( ALLOCATED( surf_usm_v(l)%tt_green_m ) ) surf_usm_v(l)%tt_green_m = 0.0_wp
1101	ENDDO
1102
1103	!-- allocate wall heat flux output array and set initial values. For horizontal surfaces
1104	! ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns) ) !can be removed
1105	ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
1106	ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
1107	ALLOCATE ( surf_usm_h%wghf_eb_window(1:surf_usm_h%ns) )
1108	ALLOCATE ( surf_usm_h%wghf_eb_green(1:surf_usm_h%ns) )
1109	ALLOCATE ( surf_usm_h%iwghf_eb(1:surf_usm_h%ns) )
1110	ALLOCATE ( surf_usm_h%iwghf_eb_window(1:surf_usm_h%ns) )
1111	IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp
1112	IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp
1113	IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp
1114	IF ( ALLOCATED( surf_usm_h%wghf_eb_window ) ) surf_usm_h%wghf_eb_window = 0.0_wp
1115	IF ( ALLOCATED( surf_usm_h%wghf_eb_green ) ) surf_usm_h%wghf_eb_green = 0.0_wp
1116	IF ( ALLOCATED( surf_usm_h%iwghf_eb ) ) surf_usm_h%iwghf_eb = 0.0_wp
1117	IF ( ALLOCATED( surf_usm_h%iwghf_eb_window ) ) surf_usm_h%iwghf_eb_window = 0.0_wp
1118	!
1119	!-- Now, for vertical surfaces
1120	DO l = 0, 3
1121	! ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns) ) ! can be removed
1122	ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
1123	ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
1124	ALLOCATE ( surf_usm_v(l)%wghf_eb_window(1:surf_usm_v(l)%ns) )
1125	ALLOCATE ( surf_usm_v(l)%wghf_eb_green(1:surf_usm_v(l)%ns) )
1126	ALLOCATE ( surf_usm_v(l)%iwghf_eb(1:surf_usm_v(l)%ns) )
1127	ALLOCATE ( surf_usm_v(l)%iwghf_eb_window(1:surf_usm_v(l)%ns) )
1128	IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp
1129	IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp
1130	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp
1131	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_window ) ) surf_usm_v(l)%wghf_eb_window = 0.0_wp
1132	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_green ) ) surf_usm_v(l)%wghf_eb_green = 0.0_wp
1133	IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb ) ) surf_usm_v(l)%iwghf_eb = 0.0_wp
1134	IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) ) surf_usm_v(l)%iwghf_eb_window = 0.0_wp
1135	ENDDO
1136
1137	END SUBROUTINE usm_allocate_surface
1138
1139
1140	!------------------------------------------------------------------------------!
1141	! Description:
1142	! ------------
1143	!> Sum up and time-average urban surface output quantities as well as allocate
1144	!> the array necessary for storing the average.
1145	!------------------------------------------------------------------------------!
1146	SUBROUTINE usm_average_3d_data( mode, variable )
1147
1148	IMPLICIT NONE
1149
1150	CHARACTER (len=*), INTENT(IN) :: mode
1151	CHARACTER (len=*), INTENT(IN) :: variable
1152
1153	INTEGER(iwp) :: i, j, k, l, m, ids, iwl,istat
1154	CHARACTER (len=varnamelength) :: var, surfid
1155	INTEGER(iwp), PARAMETER :: nd = 5
1156	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
1157
1158	!-- find the real name of the variable
1159	var = TRIM(variable)
1160	DO i = 0, nd-1
1161	k = len(TRIM(var))
1162	j = len(TRIM(dirname(i)))
1163	IF ( var(k-j+1:k) == dirname(i) ) THEN
1164	ids = i
1165	var = var(:k-j)
1166	EXIT
1167	ENDIF
1168	ENDDO
1169	IF ( ids == -1 ) THEN
1170	var = TRIM(variable)
1171	ENDIF
1172	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
1173	!-- wall layers
1174	READ(var(12:12), '(I1)', iostat=istat ) iwl
1175	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1176	var = var(1:10)
1177	ELSE
1178	!-- wrong wall layer index
1179	RETURN
1180	ENDIF
1181	ENDIF
1182	IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN
1183	!-- wall layers
1184	READ(var(14:14), '(I1)', iostat=istat ) iwl
1185	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1186	var = var(1:12)
1187	ELSE
1188	!-- wrong window layer index
1189	RETURN
1190	ENDIF
1191	ENDIF
1192	IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN
1193	!-- wall layers
1194	READ(var(13:13), '(I1)', iostat=istat ) iwl
1195	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1196	var = var(1:11)
1197	ELSE
1198	!-- wrong green layer index
1199	RETURN
1200	ENDIF
1201	ENDIF
1202
1203	IF ( mode == 'allocate' ) THEN
1204
1205	SELECT CASE ( TRIM( var ) )
1206
1207	CASE ( 'usm_rad_net' )
1208	!-- array of complete radiation balance
1209	IF ( .NOT. ALLOCATED(surf_usm_h%rad_net_av) ) THEN
1210	ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
1211	surf_usm_h%rad_net_av = 0.0_wp
1212	ENDIF
1213	DO l = 0, 3
1214	IF ( .NOT. ALLOCATED(surf_usm_v(l)%rad_net_av) ) THEN
1215	ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
1216	surf_usm_v(l)%rad_net_av = 0.0_wp
1217	ENDIF
1218	ENDDO
1219
1220	CASE ( 'usm_rad_insw' )
1221	!-- array of sw radiation falling to surface after i-th reflection
1222	IF ( .NOT. ALLOCATED(surf_usm_h%surfinsw_av) ) THEN
1223	ALLOCATE( surf_usm_h%surfinsw_av(1:surf_usm_h%ns) )
1224	surf_usm_h%surfinsw_av = 0.0_wp
1225	ENDIF
1226	DO l = 0, 3
1227	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinsw_av) ) THEN
1228	ALLOCATE( surf_usm_v(l)%surfinsw_av(1:surf_usm_v(l)%ns) )
1229	surf_usm_v(l)%surfinsw_av = 0.0_wp
1230	ENDIF
1231	ENDDO
1232
1233	CASE ( 'usm_rad_inlw' )
1234	!-- array of lw radiation falling to surface after i-th reflection
1235	IF ( .NOT. ALLOCATED(surf_usm_h%surfinlw_av) ) THEN
1236	ALLOCATE( surf_usm_h%surfinlw_av(1:surf_usm_h%ns) )
1237	surf_usm_h%surfinlw_av = 0.0_wp
1238	ENDIF
1239	DO l = 0, 3
1240	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfinlw_av) ) THEN
1241	ALLOCATE( surf_usm_v(l)%surfinlw_av(1:surf_usm_v(l)%ns) )
1242	surf_usm_v(l)%surfinlw_av = 0.0_wp
1243	ENDIF
1244	ENDDO
1245
1246	CASE ( 'usm_rad_inswdir' )
1247	!-- array of direct sw radiation falling to surface from sun
1248	IF ( .NOT. ALLOCATED(surfinswdir_av) ) THEN
1249	ALLOCATE( surfinswdir_av(startenergy:endenergy) )
1250	surfinswdir_av = 0.0_wp
1251	ENDIF
1252
1253	CASE ( 'usm_rad_inswdif' )
1254	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1255	IF ( .NOT. ALLOCATED(surfinswdif_av) ) THEN
1256	ALLOCATE( surfinswdif_av(startenergy:endenergy) )
1257	surfinswdif_av = 0.0_wp
1258	ENDIF
1259
1260	CASE ( 'usm_rad_inswref' )
1261	!-- array of sw radiation falling to surface from reflections
1262	IF ( .NOT. ALLOCATED(surfinswref_av) ) THEN
1263	ALLOCATE( surfinswref_av(startenergy:endenergy) )
1264	surfinswref_av = 0.0_wp
1265	ENDIF
1266
1267	CASE ( 'usm_rad_inlwdif' )
1268	!-- array of sw radiation falling to surface after i-th reflection
1269	IF ( .NOT. ALLOCATED(surfinlwdif_av) ) THEN
1270	ALLOCATE( surfinlwdif_av(startenergy:endenergy) )
1271	surfinlwdif_av = 0.0_wp
1272	ENDIF
1273
1274	CASE ( 'usm_rad_inlwref' )
1275	!-- array of lw radiation falling to surface from reflections
1276	IF ( .NOT. ALLOCATED(surfinlwref_av) ) THEN
1277	ALLOCATE( surfinlwref_av(startenergy:endenergy) )
1278	surfinlwref_av = 0.0_wp
1279	ENDIF
1280
1281	CASE ( 'usm_rad_outsw' )
1282	!-- array of sw radiation emitted from surface after i-th reflection
1283	IF ( .NOT. ALLOCATED(surfoutsw_av) ) THEN
1284	ALLOCATE( surfoutsw_av(startenergy:endenergy) )
1285	surfoutsw_av = 0.0_wp
1286	ENDIF
1287
1288	CASE ( 'usm_rad_outlw' )
1289	!-- array of lw radiation emitted from surface after i-th reflection
1290	IF ( .NOT. ALLOCATED(surfoutlw_av) ) THEN
1291	ALLOCATE( surfoutlw_av(startenergy:endenergy) )
1292	surfoutlw_av = 0.0_wp
1293	ENDIF
1294	CASE ( 'usm_rad_ressw' )
1295	!-- array of residua of sw radiation absorbed in surface after last reflection
1296	IF ( .NOT. ALLOCATED(surfins_av) ) THEN
1297	ALLOCATE( surfins_av(startenergy:endenergy) )
1298	surfins_av = 0.0_wp
1299	ENDIF
1300
1301	CASE ( 'usm_rad_reslw' )
1302	!-- array of residua of lw radiation absorbed in surface after last reflection
1303	IF ( .NOT. ALLOCATED(surfinl_av) ) THEN
1304	ALLOCATE( surfinl_av(startenergy:endenergy) )
1305	surfinl_av = 0.0_wp
1306	ENDIF
1307
1308	CASE ( 'usm_rad_hf' )
1309	!-- array of heat flux from radiation for surfaces after i-th reflection
1310	IF ( .NOT. ALLOCATED(surf_usm_h%surfhf_av) ) THEN
1311	ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
1312	surf_usm_h%surfhf_av = 0.0_wp
1313	ENDIF
1314	DO l = 0, 3
1315	IF ( .NOT. ALLOCATED(surf_usm_v(l)%surfhf_av) ) THEN
1316	ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
1317	surf_usm_v(l)%surfhf_av = 0.0_wp
1318	ENDIF
1319	ENDDO
1320
1321	CASE ( 'usm_wshf' )
1322	!-- array of sensible heat flux from surfaces
1323	!-- land surfaces
1324	IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN
1325	ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1326	surf_usm_h%wshf_eb_av = 0.0_wp
1327	ENDIF
1328	DO l = 0, 3
1329	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN
1330	ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1331	surf_usm_v(l)%wshf_eb_av = 0.0_wp
1332	ENDIF
1333	ENDDO
1334
1335	CASE ( 'usm_wghf' )
1336	!-- array of heat flux from ground (wall, roof, land)
1337	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN
1338	ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1339	surf_usm_h%wghf_eb_av = 0.0_wp
1340	ENDIF
1341	DO l = 0, 3
1342	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN
1343	ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1344	surf_usm_v(l)%wghf_eb_av = 0.0_wp
1345	ENDIF
1346	ENDDO
1347
1348	CASE ( 'usm_wghf_window' )
1349	!-- array of heat flux from window ground (wall, roof, land)
1350	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN
1351	ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
1352	surf_usm_h%wghf_eb_window_av = 0.0_wp
1353	ENDIF
1354	DO l = 0, 3
1355	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) ) THEN
1356	ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
1357	surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
1358	ENDIF
1359	ENDDO
1360
1361	CASE ( 'usm_wghf_green' )
1362	!-- array of heat flux from green ground (wall, roof, land)
1363	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN
1364	ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
1365	surf_usm_h%wghf_eb_green_av = 0.0_wp
1366	ENDIF
1367	DO l = 0, 3
1368	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) ) THEN
1369	ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
1370	surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
1371	ENDIF
1372	ENDDO
1373
1374	CASE ( 'usm_iwghf' )
1375	!-- array of heat flux from indoor ground (wall, roof, land)
1376	IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN
1377	ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
1378	surf_usm_h%iwghf_eb_av = 0.0_wp
1379	ENDIF
1380	DO l = 0, 3
1381	IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_av) ) THEN
1382	ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
1383	surf_usm_v(l)%iwghf_eb_av = 0.0_wp
1384	ENDIF
1385	ENDDO
1386
1387	CASE ( 'usm_iwghf_window' )
1388	!-- array of heat flux from indoor window ground (wall, roof, land)
1389	IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN
1390	ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
1391	surf_usm_h%iwghf_eb_window_av = 0.0_wp
1392	ENDIF
1393	DO l = 0, 3
1394	IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) ) THEN
1395	ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
1396	surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
1397	ENDIF
1398	ENDDO
1399
1400	CASE ( 'usm_t_surf' )
1401	!-- surface temperature for surfaces
1402	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_av) ) THEN
1403	ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
1404	surf_usm_h%t_surf_av = 0.0_wp
1405	ENDIF
1406	DO l = 0, 3
1407	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_av) ) THEN
1408	ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
1409	surf_usm_v(l)%t_surf_av = 0.0_wp
1410	ENDIF
1411	ENDDO
1412
1413	CASE ( 'usm_t_surf_window' )
1414	!-- surface temperature for window surfaces
1415	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN
1416	ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
1417	surf_usm_h%t_surf_window_av = 0.0_wp
1418	ENDIF
1419	DO l = 0, 3
1420	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_window_av) ) THEN
1421	ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
1422	surf_usm_v(l)%t_surf_window_av = 0.0_wp
1423	ENDIF
1424	ENDDO
1425
1426	CASE ( 'usm_t_surf_green' )
1427	!-- surface temperature for green surfaces
1428	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN
1429	ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
1430	surf_usm_h%t_surf_green_av = 0.0_wp
1431	ENDIF
1432	DO l = 0, 3
1433	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_green_av) ) THEN
1434	ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
1435	surf_usm_v(l)%t_surf_green_av = 0.0_wp
1436	ENDIF
1437	ENDDO
1438
1439	CASE ( 'usm_t_surf_10cm' )
1440	!-- near surface temperature for whole surfaces
1441	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_10cm_av) ) THEN
1442	ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
1443	surf_usm_h%t_surf_10cm_av = 0.0_wp
1444	ENDIF
1445	DO l = 0, 3
1446	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) ) THEN
1447	ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) )
1448	surf_usm_v(l)%t_surf_10cm_av = 0.0_wp
1449	ENDIF
1450	ENDDO
1451
1452	CASE ( 'usm_t_wall' )
1453	!-- wall temperature for iwl layer of walls and land
1454	IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN
1455	ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1456	surf_usm_h%t_wall_av = 0.0_wp
1457	ENDIF
1458	DO l = 0, 3
1459	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN
1460	ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1461	surf_usm_v(l)%t_wall_av = 0.0_wp
1462	ENDIF
1463	ENDDO
1464
1465	CASE ( 'usm_t_window' )
1466	!-- window temperature for iwl layer of walls and land
1467	IF ( .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN
1468	ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1469	surf_usm_h%t_window_av = 0.0_wp
1470	ENDIF
1471	DO l = 0, 3
1472	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_window_av) ) THEN
1473	ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1474	surf_usm_v(l)%t_window_av = 0.0_wp
1475	ENDIF
1476	ENDDO
1477
1478	CASE ( 'usm_t_green' )
1479	!-- green temperature for iwl layer of walls and land
1480	IF ( .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN
1481	ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1482	surf_usm_h%t_green_av = 0.0_wp
1483	ENDIF
1484	DO l = 0, 3
1485	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_green_av) ) THEN
1486	ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1487	surf_usm_v(l)%t_green_av = 0.0_wp
1488	ENDIF
1489	ENDDO
1490
1491	CASE DEFAULT
1492	CONTINUE
1493
1494	END SELECT
1495
1496	ELSEIF ( mode == 'sum' ) THEN
1497
1498	SELECT CASE ( TRIM( var ) )
1499
1500	CASE ( 'usm_rad_net' )
1501	!-- array of complete radiation balance
1502	DO m = 1, surf_usm_h%ns
1503	surf_usm_h%rad_net_av(m) = &
1504	surf_usm_h%rad_net_av(m) + &
1505	surf_usm_h%rad_net_l(m)
1506	ENDDO
1507	DO l = 0, 3
1508	DO m = 1, surf_usm_v(l)%ns
1509	surf_usm_v(l)%rad_net_av(m) = &
1510	surf_usm_v(l)%rad_net_av(m) + &
1511	surf_usm_v(l)%rad_net_l(m)
1512	ENDDO
1513	ENDDO
1514
1515	CASE ( 'usm_rad_insw' )
1516	!-- array of sw radiation falling to surface after i-th reflection
1517	DO l = startenergy, endenergy
1518	IF ( surfl(id,l) == ids ) THEN
1519	surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
1520	ENDIF
1521	ENDDO
1522
1523	CASE ( 'usm_rad_inlw' )
1524	!-- array of lw radiation falling to surface after i-th reflection
1525	DO l = startenergy, endenergy
1526	IF ( surfl(id,l) == ids ) THEN
1527	surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
1528	ENDIF
1529	ENDDO
1530
1531	CASE ( 'usm_rad_inswdir' )
1532	!-- array of direct sw radiation falling to surface from sun
1533	DO l = startenergy, endenergy
1534	IF ( surfl(id,l) == ids ) THEN
1535	surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
1536	ENDIF
1537	ENDDO
1538
1539	CASE ( 'usm_rad_inswdif' )
1540	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1541	DO l = startenergy, endenergy
1542	IF ( surfl(id,l) == ids ) THEN
1543	surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
1544	ENDIF
1545	ENDDO
1546
1547	CASE ( 'usm_rad_inswref' )
1548	!-- array of sw radiation falling to surface from reflections
1549	DO l = startenergy, endenergy
1550	IF ( surfl(id,l) == ids ) THEN
1551	surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1552	surfinswdir(l) - surfinswdif(l)
1553	ENDIF
1554	ENDDO
1555
1556
1557	CASE ( 'usm_rad_inlwdif' )
1558	!-- array of sw radiation falling to surface after i-th reflection
1559	DO l = startenergy, endenergy
1560	IF ( surfl(id,l) == ids ) THEN
1561	surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1562	surfinswdir(l) - surfinswdif(l)
1563	ENDIF
1564	ENDDO
1565	!
1566	CASE ( 'usm_rad_inlwref' )
1567	!-- array of lw radiation falling to surface from reflections
1568	DO l = startenergy, endenergy
1569	IF ( surfl(id,l) == ids ) THEN
1570	surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
1571	ENDIF
1572	ENDDO
1573
1574	CASE ( 'usm_rad_outsw' )
1575	!-- array of sw radiation emitted from surface after i-th reflection
1576	DO l = startenergy, endenergy
1577	IF ( surfl(id,l) == ids ) THEN
1578	surfinlwref_av(l) = surfinlwref_av(l) + &
1579	surfinlw(l) - surfinlwdif(l)
1580	ENDIF
1581	ENDDO
1582
1583	CASE ( 'usm_rad_outlw' )
1584	!-- array of lw radiation emitted from surface after i-th reflection
1585	DO l = startenergy, endenergy
1586	IF ( surfl(id,l) == ids ) THEN
1587	surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
1588	ENDIF
1589	ENDDO
1590
1591	CASE ( 'usm_rad_ressw' )
1592	!-- array of residua of sw radiation absorbed in surface after last reflection
1593	DO l = startenergy, endenergy
1594	IF ( surfl(id,l) == ids ) THEN
1595	surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
1596	ENDIF
1597	ENDDO
1598
1599	CASE ( 'usm_rad_reslw' )
1600	!-- array of residua of lw radiation absorbed in surface after last reflection
1601	DO l = startenergy, endenergy
1602	IF ( surfl(id,l) == ids ) THEN
1603	surfins_av(l) = surfins_av(l) + surfins(l)
1604	ENDIF
1605	ENDDO
1606
1607	CASE ( 'usm_rad_hf' )
1608	!-- array of heat flux from radiation for surfaces after i-th reflection
1609	DO m = 1, surf_usm_h%ns
1610	surf_usm_h%surfhf_av(m) = &
1611	surf_usm_h%surfhf_av(m) + &
1612	surf_usm_h%surfhf(m)
1613	ENDDO
1614	DO l = 0, 3
1615	DO m = 1, surf_usm_v(l)%ns
1616	surf_usm_v(l)%surfhf_av(m) = &
1617	surf_usm_v(l)%surfhf_av(m) + &
1618	surf_usm_v(l)%surfhf(m)
1619	ENDDO
1620	ENDDO
1621
1622	CASE ( 'usm_wshf' )
1623	!-- array of sensible heat flux from surfaces (land, roof, wall)
1624	DO m = 1, surf_usm_h%ns
1625	surf_usm_h%wshf_eb_av(m) = &
1626	surf_usm_h%wshf_eb_av(m) + &
1627	surf_usm_h%wshf_eb(m)
1628	ENDDO
1629	DO l = 0, 3
1630	DO m = 1, surf_usm_v(l)%ns
1631	surf_usm_v(l)%wshf_eb_av(m) = &
1632	surf_usm_v(l)%wshf_eb_av(m) + &
1633	surf_usm_v(l)%wshf_eb(m)
1634	ENDDO
1635	ENDDO
1636
1637	CASE ( 'usm_wghf' )
1638	!-- array of heat flux from ground (wall, roof, land)
1639	DO m = 1, surf_usm_h%ns
1640	surf_usm_h%wghf_eb_av(m) = &
1641	surf_usm_h%wghf_eb_av(m) + &
1642	surf_usm_h%wghf_eb(m)
1643	ENDDO
1644	DO l = 0, 3
1645	DO m = 1, surf_usm_v(l)%ns
1646	surf_usm_v(l)%wghf_eb_av(m) = &
1647	surf_usm_v(l)%wghf_eb_av(m) + &
1648	surf_usm_v(l)%wghf_eb(m)
1649	ENDDO
1650	ENDDO
1651
1652	CASE ( 'usm_wghf_window' )
1653	!-- array of heat flux from window ground (wall, roof, land)
1654	DO m = 1, surf_usm_h%ns
1655	surf_usm_h%wghf_eb_window_av(m) = &
1656	surf_usm_h%wghf_eb_window_av(m) + &
1657	surf_usm_h%wghf_eb_window(m)
1658	ENDDO
1659	DO l = 0, 3
1660	DO m = 1, surf_usm_v(l)%ns
1661	surf_usm_v(l)%wghf_eb_window_av(m) = &
1662	surf_usm_v(l)%wghf_eb_window_av(m) + &
1663	surf_usm_v(l)%wghf_eb_window(m)
1664	ENDDO
1665	ENDDO
1666
1667	CASE ( 'usm_wghf_green' )
1668	!-- array of heat flux from green ground (wall, roof, land)
1669	DO m = 1, surf_usm_h%ns
1670	surf_usm_h%wghf_eb_green_av(m) = &
1671	surf_usm_h%wghf_eb_green_av(m) + &
1672	surf_usm_h%wghf_eb_green(m)
1673	ENDDO
1674	DO l = 0, 3
1675	DO m = 1, surf_usm_v(l)%ns
1676	surf_usm_v(l)%wghf_eb_green_av(m) = &
1677	surf_usm_v(l)%wghf_eb_green_av(m) + &
1678	surf_usm_v(l)%wghf_eb_green(m)
1679	ENDDO
1680	ENDDO
1681
1682	CASE ( 'usm_iwghf' )
1683	!-- array of heat flux from indoor ground (wall, roof, land)
1684	DO m = 1, surf_usm_h%ns
1685	surf_usm_h%iwghf_eb_av(m) = &
1686	surf_usm_h%iwghf_eb_av(m) + &
1687	surf_usm_h%iwghf_eb(m)
1688	ENDDO
1689	DO l = 0, 3
1690	DO m = 1, surf_usm_v(l)%ns
1691	surf_usm_v(l)%iwghf_eb_av(m) = &
1692	surf_usm_v(l)%iwghf_eb_av(m) + &
1693	surf_usm_v(l)%iwghf_eb(m)
1694	ENDDO
1695	ENDDO
1696
1697	CASE ( 'usm_iwghf_window' )
1698	!-- array of heat flux from indoor window ground (wall, roof, land)
1699	DO m = 1, surf_usm_h%ns
1700	surf_usm_h%iwghf_eb_window_av(m) = &
1701	surf_usm_h%iwghf_eb_window_av(m) + &
1702	surf_usm_h%iwghf_eb_window(m)
1703	ENDDO
1704	DO l = 0, 3
1705	DO m = 1, surf_usm_v(l)%ns
1706	surf_usm_v(l)%iwghf_eb_window_av(m) = &
1707	surf_usm_v(l)%iwghf_eb_window_av(m) + &
1708	surf_usm_v(l)%iwghf_eb_window(m)
1709	ENDDO
1710	ENDDO
1711
1712	CASE ( 'usm_t_surf' )
1713	!-- surface temperature for surfaces
1714	DO m = 1, surf_usm_h%ns
1715	surf_usm_h%t_surf_av(m) = &
1716	surf_usm_h%t_surf_av(m) + &
1717	t_surf_h(m)
1718	ENDDO
1719	DO l = 0, 3
1720	DO m = 1, surf_usm_v(l)%ns
1721	surf_usm_v(l)%t_surf_av(m) = &
1722	surf_usm_v(l)%t_surf_av(m) + &
1723	t_surf_v(l)%t(m)
1724	ENDDO
1725	ENDDO
1726
1727	CASE ( 'usm_t_surf_window' )
1728	!-- surface temperature for window surfaces
1729	DO m = 1, surf_usm_h%ns
1730	surf_usm_h%t_surf_window_av(m) = &
1731	surf_usm_h%t_surf_window_av(m) + &
1732	t_surf_window_h(m)
1733	ENDDO
1734	DO l = 0, 3
1735	DO m = 1, surf_usm_v(l)%ns
1736	surf_usm_v(l)%t_surf_window_av(m) = &
1737	surf_usm_v(l)%t_surf_window_av(m) + &
1738	t_surf_window_v(l)%t(m)
1739	ENDDO
1740	ENDDO
1741
1742	CASE ( 'usm_t_surf_green' )
1743	!-- surface temperature for green surfaces
1744	DO m = 1, surf_usm_h%ns
1745	surf_usm_h%t_surf_green_av(m) = &
1746	surf_usm_h%t_surf_green_av(m) + &
1747	t_surf_green_h(m)
1748	ENDDO
1749	DO l = 0, 3
1750	DO m = 1, surf_usm_v(l)%ns
1751	surf_usm_v(l)%t_surf_green_av(m) = &
1752	surf_usm_v(l)%t_surf_green_av(m) + &
1753	t_surf_green_v(l)%t(m)
1754	ENDDO
1755	ENDDO
1756
1757	CASE ( 'usm_t_surf_10cm' )
1758	!-- near surface temperature for whole surfaces
1759	DO m = 1, surf_usm_h%ns
1760	surf_usm_h%t_surf_10cm_av(m) = &
1761	surf_usm_h%t_surf_10cm_av(m) + &
1762	t_surf_10cm_h(m)
1763	ENDDO
1764	DO l = 0, 3
1765	DO m = 1, surf_usm_v(l)%ns
1766	surf_usm_v(l)%t_surf_10cm_av(m) = &
1767	surf_usm_v(l)%t_surf_10cm_av(m) + &
1768	t_surf_10cm_v(l)%t(m)
1769	ENDDO
1770	ENDDO
1771
1772
1773	CASE ( 'usm_t_wall' )
1774	!-- wall temperature for iwl layer of walls and land
1775	DO m = 1, surf_usm_h%ns
1776	surf_usm_h%t_wall_av(iwl,m) = &
1777	surf_usm_h%t_wall_av(iwl,m) + &
1778	t_wall_h(iwl,m)
1779	ENDDO
1780	DO l = 0, 3
1781	DO m = 1, surf_usm_v(l)%ns
1782	surf_usm_v(l)%t_wall_av(iwl,m) = &
1783	surf_usm_v(l)%t_wall_av(iwl,m) + &
1784	t_wall_v(l)%t(iwl,m)
1785	ENDDO
1786	ENDDO
1787
1788	CASE ( 'usm_t_window' )
1789	!-- window temperature for iwl layer of walls and land
1790	DO m = 1, surf_usm_h%ns
1791	surf_usm_h%t_window_av(iwl,m) = &
1792	surf_usm_h%t_window_av(iwl,m) + &
1793	t_window_h(iwl,m)
1794	ENDDO
1795	DO l = 0, 3
1796	DO m = 1, surf_usm_v(l)%ns
1797	surf_usm_v(l)%t_window_av(iwl,m) = &
1798	surf_usm_v(l)%t_window_av(iwl,m) + &
1799	t_window_v(l)%t(iwl,m)
1800	ENDDO
1801	ENDDO
1802
1803	CASE ( 'usm_t_green' )
1804	!-- green temperature for iwl layer of walls and land
1805	DO m = 1, surf_usm_h%ns
1806	surf_usm_h%t_green_av(iwl,m) = &
1807	surf_usm_h%t_green_av(iwl,m) + &
1808	t_green_h(iwl,m)
1809	ENDDO
1810	DO l = 0, 3
1811	DO m = 1, surf_usm_v(l)%ns
1812	surf_usm_v(l)%t_green_av(iwl,m) = &
1813	surf_usm_v(l)%t_green_av(iwl,m) + &
1814	t_green_v(l)%t(iwl,m)
1815	ENDDO
1816	ENDDO
1817
1818	CASE DEFAULT
1819	CONTINUE
1820
1821	END SELECT
1822
1823	ELSEIF ( mode == 'average' ) THEN
1824
1825	SELECT CASE ( TRIM( var ) )
1826
1827	CASE ( 'usm_rad_net' )
1828	!-- array of complete radiation balance
1829	DO m = 1, surf_usm_h%ns
1830	surf_usm_h%rad_net_av(m) = &
1831	surf_usm_h%rad_net_av(m) / &
1832	REAL( average_count_3d, kind=wp )
1833	ENDDO
1834	DO l = 0, 3
1835	DO m = 1, surf_usm_v(l)%ns
1836	surf_usm_v(l)%rad_net_av(m) = &
1837	surf_usm_v(l)%rad_net_av(m) / &
1838	REAL( average_count_3d, kind=wp )
1839	ENDDO
1840	ENDDO
1841
1842	CASE ( 'usm_rad_insw' )
1843	!-- array of sw radiation falling to surface after i-th reflection
1844	DO l = startenergy, endenergy
1845	IF ( surfl(id,l) == ids ) THEN
1846	surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
1847	ENDIF
1848	ENDDO
1849
1850	CASE ( 'usm_rad_inlw' )
1851	!-- array of lw radiation falling to surface after i-th reflection
1852	DO l = startenergy, endenergy
1853	IF ( surfl(id,l) == ids ) THEN
1854	surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
1855	ENDIF
1856	ENDDO
1857
1858	CASE ( 'usm_rad_inswdir' )
1859	!-- array of direct sw radiation falling to surface from sun
1860	DO l = startenergy, endenergy
1861	IF ( surfl(id,l) == ids ) THEN
1862	surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
1863	ENDIF
1864	ENDDO
1865
1866	CASE ( 'usm_rad_inswdif' )
1867	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
1868	DO l = startenergy, endenergy
1869	IF ( surfl(id,l) == ids ) THEN
1870	surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
1871	ENDIF
1872	ENDDO
1873
1874	CASE ( 'usm_rad_inswref' )
1875	!-- array of sw radiation falling to surface from reflections
1876	DO l = startenergy, endenergy
1877	IF ( surfl(id,l) == ids ) THEN
1878	surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
1879	ENDIF
1880	ENDDO
1881
1882	CASE ( 'usm_rad_inlwdif' )
1883	!-- array of sw radiation falling to surface after i-th reflection
1884	DO l = startenergy, endenergy
1885	IF ( surfl(id,l) == ids ) THEN
1886	surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
1887	ENDIF
1888	ENDDO
1889
1890	CASE ( 'usm_rad_inlwref' )
1891	!-- array of lw radiation falling to surface from reflections
1892	DO l = startenergy, endenergy
1893	IF ( surfl(id,l) == ids ) THEN
1894	surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
1895	ENDIF
1896	ENDDO
1897
1898	CASE ( 'usm_rad_outsw' )
1899	!-- array of sw radiation emitted from surface after i-th reflection
1900	DO l = startenergy, endenergy
1901	IF ( surfl(id,l) == ids ) THEN
1902	surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
1903	ENDIF
1904	ENDDO
1905
1906	CASE ( 'usm_rad_outlw' )
1907	!-- array of lw radiation emitted from surface after i-th reflection
1908	DO l = startenergy, endenergy
1909	IF ( surfl(id,l) == ids ) THEN
1910	surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
1911	ENDIF
1912	ENDDO
1913
1914	CASE ( 'usm_rad_ressw' )
1915	!-- array of residua of sw radiation absorbed in surface after last reflection
1916	DO l = startenergy, endenergy
1917	IF ( surfl(id,l) == ids ) THEN
1918	surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
1919	ENDIF
1920	ENDDO
1921
1922	CASE ( 'usm_rad_reslw' )
1923	!-- array of residua of lw radiation absorbed in surface after last reflection
1924	DO l = startenergy, endenergy
1925	IF ( surfl(id,l) == ids ) THEN
1926	surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
1927	ENDIF
1928	ENDDO
1929
1930	CASE ( 'usm_rad_hf' )
1931	!-- array of heat flux from radiation for surfaces after i-th reflection
1932	DO m = 1, surf_usm_h%ns
1933	surf_usm_h%surfhf_av(m) = &
1934	surf_usm_h%surfhf_av(m) / &
1935	REAL( average_count_3d, kind=wp )
1936	ENDDO
1937	DO l = 0, 3
1938	DO m = 1, surf_usm_v(l)%ns
1939	surf_usm_v(l)%surfhf_av(m) = &
1940	surf_usm_v(l)%surfhf_av(m) / &
1941	REAL( average_count_3d, kind=wp )
1942	ENDDO
1943	ENDDO
1944
1945	CASE ( 'usm_wshf' )
1946	!-- array of sensible heat flux from surfaces (land, roof, wall)
1947	DO m = 1, surf_usm_h%ns
1948	surf_usm_h%wshf_eb_av(m) = &
1949	surf_usm_h%wshf_eb_av(m) / &
1950	REAL( average_count_3d, kind=wp )
1951	ENDDO
1952	DO l = 0, 3
1953	DO m = 1, surf_usm_v(l)%ns
1954	surf_usm_v(l)%wshf_eb_av(m) = &
1955	surf_usm_v(l)%wshf_eb_av(m) / &
1956	REAL( average_count_3d, kind=wp )
1957	ENDDO
1958	ENDDO
1959
1960	CASE ( 'usm_wghf' )
1961	!-- array of heat flux from ground (wall, roof, land)
1962	DO m = 1, surf_usm_h%ns
1963	surf_usm_h%wghf_eb_av(m) = &
1964	surf_usm_h%wghf_eb_av(m) / &
1965	REAL( average_count_3d, kind=wp )
1966	ENDDO
1967	DO l = 0, 3
1968	DO m = 1, surf_usm_v(l)%ns
1969	surf_usm_v(l)%wghf_eb_av(m) = &
1970	surf_usm_v(l)%wghf_eb_av(m) / &
1971	REAL( average_count_3d, kind=wp )
1972	ENDDO
1973	ENDDO
1974
1975	CASE ( 'usm_wghf_window' )
1976	!-- array of heat flux from window ground (wall, roof, land)
1977	DO m = 1, surf_usm_h%ns
1978	surf_usm_h%wghf_eb_window_av(m) = &
1979	surf_usm_h%wghf_eb_window_av(m) / &
1980	REAL( average_count_3d, kind=wp )
1981	ENDDO
1982	DO l = 0, 3
1983	DO m = 1, surf_usm_v(l)%ns
1984	surf_usm_v(l)%wghf_eb_window_av(m) = &
1985	surf_usm_v(l)%wghf_eb_window_av(m) / &
1986	REAL( average_count_3d, kind=wp )
1987	ENDDO
1988	ENDDO
1989
1990	CASE ( 'usm_wghf_green' )
1991	!-- array of heat flux from green ground (wall, roof, land)
1992	DO m = 1, surf_usm_h%ns
1993	surf_usm_h%wghf_eb_green_av(m) = &
1994	surf_usm_h%wghf_eb_green_av(m) / &
1995	REAL( average_count_3d, kind=wp )
1996	ENDDO
1997	DO l = 0, 3
1998	DO m = 1, surf_usm_v(l)%ns
1999	surf_usm_v(l)%wghf_eb_green_av(m) = &
2000	surf_usm_v(l)%wghf_eb_green_av(m) / &
2001	REAL( average_count_3d, kind=wp )
2002	ENDDO
2003	ENDDO
2004
2005	CASE ( 'usm_iwghf' )
2006	!-- array of heat flux from indoor ground (wall, roof, land)
2007	DO m = 1, surf_usm_h%ns
2008	surf_usm_h%iwghf_eb_av(m) = &
2009	surf_usm_h%iwghf_eb_av(m) / &
2010	REAL( average_count_3d, kind=wp )
2011	ENDDO
2012	DO l = 0, 3
2013	DO m = 1, surf_usm_v(l)%ns
2014	surf_usm_v(l)%iwghf_eb_av(m) = &
2015	surf_usm_v(l)%iwghf_eb_av(m) / &
2016	REAL( average_count_3d, kind=wp )
2017	ENDDO
2018	ENDDO
2019
2020	CASE ( 'usm_iwghf_window' )
2021	!-- array of heat flux from indoor window ground (wall, roof, land)
2022	DO m = 1, surf_usm_h%ns
2023	surf_usm_h%iwghf_eb_window_av(m) = &
2024	surf_usm_h%iwghf_eb_window_av(m) / &
2025	REAL( average_count_3d, kind=wp )
2026	ENDDO
2027	DO l = 0, 3
2028	DO m = 1, surf_usm_v(l)%ns
2029	surf_usm_v(l)%iwghf_eb_window_av(m) = &
2030	surf_usm_v(l)%iwghf_eb_window_av(m) / &
2031	REAL( average_count_3d, kind=wp )
2032	ENDDO
2033	ENDDO
2034
2035	CASE ( 'usm_t_surf' )
2036	!-- surface temperature for surfaces
2037	DO m = 1, surf_usm_h%ns
2038	surf_usm_h%t_surf_av(m) = &
2039	surf_usm_h%t_surf_av(m) / &
2040	REAL( average_count_3d, kind=wp )
2041	ENDDO
2042	DO l = 0, 3
2043	DO m = 1, surf_usm_v(l)%ns
2044	surf_usm_v(l)%t_surf_av(m) = &
2045	surf_usm_v(l)%t_surf_av(m) / &
2046	REAL( average_count_3d, kind=wp )
2047	ENDDO
2048	ENDDO
2049
2050	CASE ( 'usm_t_surf_window' )
2051	!-- surface temperature for window surfaces
2052	DO m = 1, surf_usm_h%ns
2053	surf_usm_h%t_surf_window_av(m) = &
2054	surf_usm_h%t_surf_window_av(m) / &
2055	REAL( average_count_3d, kind=wp )
2056	ENDDO
2057	DO l = 0, 3
2058	DO m = 1, surf_usm_v(l)%ns
2059	surf_usm_v(l)%t_surf_window_av(m) = &
2060	surf_usm_v(l)%t_surf_window_av(m) / &
2061	REAL( average_count_3d, kind=wp )
2062	ENDDO
2063	ENDDO
2064
2065	CASE ( 'usm_t_surf_green' )
2066	!-- surface temperature for green surfaces
2067	DO m = 1, surf_usm_h%ns
2068	surf_usm_h%t_surf_green_av(m) = &
2069	surf_usm_h%t_surf_green_av(m) / &
2070	REAL( average_count_3d, kind=wp )
2071	ENDDO
2072	DO l = 0, 3
2073	DO m = 1, surf_usm_v(l)%ns
2074	surf_usm_v(l)%t_surf_green_av(m) = &
2075	surf_usm_v(l)%t_surf_green_av(m) / &
2076	REAL( average_count_3d, kind=wp )
2077	ENDDO
2078	ENDDO
2079
2080	CASE ( 'usm_t_surf_10cm' )
2081	!-- near surface temperature for whole surfaces
2082	DO m = 1, surf_usm_h%ns
2083	surf_usm_h%t_surf_10cm_av(m) = &
2084	surf_usm_h%t_surf_10cm_av(m) / &
2085	REAL( average_count_3d, kind=wp )
2086	ENDDO
2087	DO l = 0, 3
2088	DO m = 1, surf_usm_v(l)%ns
2089	surf_usm_v(l)%t_surf_10cm_av(m) = &
2090	surf_usm_v(l)%t_surf_10cm_av(m) / &
2091	REAL( average_count_3d, kind=wp )
2092	ENDDO
2093	ENDDO
2094
2095	CASE ( 'usm_t_wall' )
2096	!-- wall temperature for iwl layer of walls and land
2097	DO m = 1, surf_usm_h%ns
2098	surf_usm_h%t_wall_av(iwl,m) = &
2099	surf_usm_h%t_wall_av(iwl,m) / &
2100	REAL( average_count_3d, kind=wp )
2101	ENDDO
2102	DO l = 0, 3
2103	DO m = 1, surf_usm_v(l)%ns
2104	surf_usm_v(l)%t_wall_av(iwl,m) = &
2105	surf_usm_v(l)%t_wall_av(iwl,m) / &
2106	REAL( average_count_3d, kind=wp )
2107	ENDDO
2108	ENDDO
2109
2110	CASE ( 'usm_t_window' )
2111	!-- window temperature for iwl layer of walls and land
2112	DO m = 1, surf_usm_h%ns
2113	surf_usm_h%t_window_av(iwl,m) = &
2114	surf_usm_h%t_window_av(iwl,m) / &
2115	REAL( average_count_3d, kind=wp )
2116	ENDDO
2117	DO l = 0, 3
2118	DO m = 1, surf_usm_v(l)%ns
2119	surf_usm_v(l)%t_window_av(iwl,m) = &
2120	surf_usm_v(l)%t_window_av(iwl,m) / &
2121	REAL( average_count_3d, kind=wp )
2122	ENDDO
2123	ENDDO
2124
2125	CASE ( 'usm_t_green' )
2126	!-- green temperature for iwl layer of walls and land
2127	DO m = 1, surf_usm_h%ns
2128	surf_usm_h%t_green_av(iwl,m) = &
2129	surf_usm_h%t_green_av(iwl,m) / &
2130	REAL( average_count_3d, kind=wp )
2131	ENDDO
2132	DO l = 0, 3
2133	DO m = 1, surf_usm_v(l)%ns
2134	surf_usm_v(l)%t_green_av(iwl,m) = &
2135	surf_usm_v(l)%t_green_av(iwl,m) / &
2136	REAL( average_count_3d, kind=wp )
2137	ENDDO
2138	ENDDO
2139
2140
2141	END SELECT
2142
2143	ENDIF
2144
2145	END SUBROUTINE usm_average_3d_data
2146
2147
2148
2149	!------------------------------------------------------------------------------!
2150	! Description:
2151	! ------------
2152	!> Set internal Neumann boundary condition at outer soil grid points
2153	!> for temperature and humidity.
2154	!------------------------------------------------------------------------------!
2155	SUBROUTINE usm_boundary_condition
2156
2157	IMPLICIT NONE
2158
2159	INTEGER(iwp) :: i !< grid index x-direction
2160	INTEGER(iwp) :: ioff !< offset index x-direction indicating location of soil grid point
2161	INTEGER(iwp) :: j !< grid index y-direction
2162	INTEGER(iwp) :: joff !< offset index x-direction indicating location of soil grid point
2163	INTEGER(iwp) :: k !< grid index z-direction
2164	INTEGER(iwp) :: koff !< offset index x-direction indicating location of soil grid point
2165	INTEGER(iwp) :: l !< running index surface-orientation
2166	INTEGER(iwp) :: m !< running index surface elements
2167
2168	koff = surf_usm_h%koff
2169	DO m = 1, surf_usm_h%ns
2170	i = surf_usm_h%i(m)
2171	j = surf_usm_h%j(m)
2172	k = surf_usm_h%k(m)
2173	pt(k+koff,j,i) = pt(k,j,i)
2174	ENDDO
2175
2176	DO l = 0, 3
2177	ioff = surf_usm_v(l)%ioff
2178	joff = surf_usm_v(l)%joff
2179	DO m = 1, surf_usm_v(l)%ns
2180	i = surf_usm_v(l)%i(m)
2181	j = surf_usm_v(l)%j(m)
2182	k = surf_usm_v(l)%k(m)
2183	pt(k,j+joff,i+ioff) = pt(k,j,i)
2184	ENDDO
2185	ENDDO
2186
2187	END SUBROUTINE usm_boundary_condition
2188
2189
2190	!------------------------------------------------------------------------------!
2191	!
2192	! Description:
2193	! ------------
2194	!> Subroutine checks variables and assigns units.
2195	!> It is called out from subroutine check_parameters.
2196	!------------------------------------------------------------------------------!
2197	SUBROUTINE usm_check_data_output( variable, unit )
2198
2199	IMPLICIT NONE
2200
2201	CHARACTER (len=*),INTENT(IN) :: variable !:
2202	CHARACTER (len=*),INTENT(OUT) :: unit !:
2203
2204	CHARACTER (len=varnamelength) :: var
2205
2206	var = TRIM(variable)
2207	IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. &
2208	var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. &
2209	var(1:16) == 'usm_rad_inswdif_' .OR. var(1:16) == 'usm_rad_inswref_' .OR. &
2210	var(1:16) == 'usm_rad_inlwdif_' .OR. var(1:16) == 'usm_rad_inlwref_' .OR. &
2211	var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. &
2212	var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. &
2213	var(1:11) == 'usm_rad_hf_' .OR. &
2214	var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
2215	var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. &
2216	var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' ) THEN
2217	unit = 'W/m2'
2218	ELSE IF ( var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. &
2219	var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR. &
2220	var(1:16) == 'usm_t_surf_green' .OR. &
2221	var(1:11) == 'usm_t_green' .OR. &
2222	var(1:15) == 'usm_t_surf_10cm') THEN
2223	unit = 'K'
2224	ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. &
2225	var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. &
2226	var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis') THEN
2227	unit = '1'
2228	ELSE
2229	unit = 'illegal'
2230	ENDIF
2231
2232	END SUBROUTINE usm_check_data_output
2233
2234
2235	!------------------------------------------------------------------------------!
2236	! Description:
2237	! ------------
2238	!> Check parameters routine for urban surface model
2239	!------------------------------------------------------------------------------!
2240	SUBROUTINE usm_check_parameters
2241
2242	USE control_parameters, &
2243	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing, &
2244	lsf_surf, topography
2245
2246	!
2247	!-- Dirichlet boundary conditions are required as the surface fluxes are
2248	!-- calculated from the temperature/humidity gradients in the urban surface
2249	!-- model
2250	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
2251	message_string = 'urban surface model requires setting of '// &
2252	'bc_pt_b = "dirichlet" and '// &
2253	'bc_q_b = "dirichlet"'
2254	CALL message( 'check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
2255	ENDIF
2256
2257	IF ( .NOT. constant_flux_layer ) THEN
2258	message_string = 'urban surface model requires '// &
2259	'constant_flux_layer = .T.'
2260	CALL message( 'check_parameters', 'PA0591', 1, 2, 0, 6, 0 )
2261	ENDIF
2262	!
2263	!-- Surface forcing has to be disabled for LSF in case of enabled
2264	!-- urban surface module
2265	IF ( large_scale_forcing ) THEN
2266	lsf_surf = .FALSE.
2267	ENDIF
2268	!
2269	!-- Topography
2270	IF ( topography == 'flat' ) THEN
2271	message_string = 'topography /= "flat" is required '// &
2272	'when using the urban surface model'
2273	CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2274	ENDIF
2275
2276	END SUBROUTINE usm_check_parameters
2277
2278
2279	!------------------------------------------------------------------------------!
2280	!
2281	! Description:
2282	! ------------
2283	!> Output of the 3D-arrays in netCDF and/or AVS format
2284	!> for variables of urban_surface model.
2285	!> It resorts the urban surface module output quantities from surf style
2286	!> indexing into temporary 3D array with indices (i,j,k).
2287	!> It is called from subroutine data_output_3d.
2288	!------------------------------------------------------------------------------!
2289	SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2290
2291	IMPLICIT NONE
2292
2293	INTEGER(iwp), INTENT(IN) :: av !<
2294	CHARACTER (len=*), INTENT(IN) :: variable !<
2295	INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0)
2296	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
2297	LOGICAL, INTENT(OUT) :: found !<
2298	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< sp - it has to correspond to module data_output_3d
2299	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: temp_pf !< temp array for urban surface output procedure
2300
2301	CHARACTER (len=varnamelength) :: var, surfid
2302	INTEGER(iwp), PARAMETER :: nd = 5
2303	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
2304	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
2305	INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart
2306	INTEGER(iwp), DIMENSION(0:nd-1) :: dirend
2307	INTEGER(iwp) :: ids,isurf,isvf,isurfs,isurflt
2308	INTEGER(iwp) :: is,js,ks,i,j,k,iwl,istat, l, m
2309	INTEGER(iwp) :: k_topo !< topography top index
2310
2311	dirstart = (/ startland, startwall, startwall, startwall, startwall /)
2312	dirend = (/ endland, endwall, endwall, endwall, endwall /)
2313
2314	found = .TRUE.
2315	temp_pf = -1._wp
2316
2317	ids = -1
2318	var = TRIM(variable)
2319	DO i = 0, nd-1
2320	k = len(TRIM(var))
2321	j = len(TRIM(dirname(i)))
2322	IF ( var(k-j+1:k) == dirname(i) ) THEN
2323	ids = i
2324	var = var(:k-j)
2325	EXIT
2326	ENDIF
2327	ENDDO
2328	IF ( ids == -1 ) THEN
2329	var = TRIM(variable)
2330	ENDIF
2331	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
2332	!-- wall layers
2333	READ(var(12:12), '(I1)', iostat=istat ) iwl
2334	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2335	var = var(1:10)
2336	ENDIF
2337	ENDIF
2338	IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN
2339	!-- window layers
2340	READ(var(14:14), '(I1)', iostat=istat ) iwl
2341	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2342	var = var(1:12)
2343	ENDIF
2344	ENDIF
2345	IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN
2346	!-- green layers
2347	READ(var(13:13), '(I1)', iostat=istat ) iwl
2348	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2349	var = var(1:11)
2350	ENDIF
2351	ENDIF
2352	IF ( (var(1:8) == 'usm_svf_' .OR. var(1:8) == 'usm_dif_') .AND. len(TRIM(var)) >= 13 ) THEN
2353	!-- svf values to particular surface
2354	surfid = var(9:)
2355	i = index(surfid,'_')
2356	j = index(surfid(i+1:),'_')
2357	READ(surfid(1:i-1),*, iostat=istat ) is
2358	IF ( istat == 0 ) THEN
2359	READ(surfid(i+1:i+j-1),*, iostat=istat ) js
2360	ENDIF
2361	IF ( istat == 0 ) THEN
2362	READ(surfid(i+j+1:),*, iostat=istat ) ks
2363	ENDIF
2364	IF ( istat == 0 ) THEN
2365	var = var(1:7)
2366	ENDIF
2367	ENDIF
2368
2369	SELECT CASE ( TRIM(var) )
2370
2371	CASE ( 'usm_surfz' )
2372	!-- array of lw radiation falling to local surface after i-th reflection
2373	DO m = 1, surf_usm_h%ns
2374	i = surf_usm_h%i(m)
2375	j = surf_usm_h%j(m)
2376	k = surf_usm_h%k(m)
2377	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
2378	ENDDO
2379	DO l = 0, 3
2380	DO m = 1, surf_usm_v(l)%ns
2381	i = surf_usm_v(l)%i(m)
2382	j = surf_usm_v(l)%j(m)
2383	k = surf_usm_v(l)%k(m)
2384	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
2385	ENDDO
2386	ENDDO
2387
2388	CASE ( 'usm_surfcat' )
2389	!-- surface category
2390	DO m = 1, surf_usm_h%ns
2391	i = surf_usm_h%i(m)
2392	j = surf_usm_h%j(m)
2393	k = surf_usm_h%k(m)
2394	temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2395	ENDDO
2396	DO l = 0, 3
2397	DO m = 1, surf_usm_v(l)%ns
2398	i = surf_usm_v(l)%i(m)
2399	j = surf_usm_v(l)%j(m)
2400	k = surf_usm_v(l)%k(m)
2401	temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
2402	ENDDO
2403	ENDDO
2404
2405	CASE ( 'usm_surfalb' )
2406	!-- surface albedo, weighted average
2407	DO m = 1, surf_usm_h%ns
2408	i = surf_usm_h%i(m)
2409	j = surf_usm_h%j(m)
2410	k = surf_usm_h%k(m)
2411	temp_pf(k,j,i) = surf_usm_h%frac(0,m) * &
2412	surf_usm_h%albedo(0,m) + &
2413	surf_usm_h%frac(1,m) * &
2414	surf_usm_h%albedo(1,m) + &
2415	surf_usm_h%frac(2,m) * &
2416	surf_usm_h%albedo(2,m)
2417	ENDDO
2418	DO l = 0, 3
2419	DO m = 1, surf_usm_v(l)%ns
2420	i = surf_usm_v(l)%i(m)
2421	j = surf_usm_v(l)%j(m)
2422	k = surf_usm_v(l)%k(m)
2423	temp_pf(k,j,i) = surf_usm_v(l)%frac(0,m) * &
2424	surf_usm_v(l)%albedo(0,m) + &
2425	surf_usm_v(l)%frac(1,m) * &
2426	surf_usm_v(l)%albedo(1,m) + &
2427	surf_usm_v(l)%frac(2,m) * &
2428	surf_usm_v(l)%albedo(2,m)
2429	ENDDO
2430	ENDDO
2431
2432	CASE ( 'usm_surfemis' )
2433	!-- surface emissivity, weighted average
2434	DO m = 1, surf_usm_h%ns
2435	i = surf_usm_h%i(m)
2436	j = surf_usm_h%j(m)
2437	k = surf_usm_h%k(m)
2438	temp_pf(k,j,i) = surf_usm_h%frac(0,m) * &
2439	surf_usm_h%emissivity(0,m) + &
2440	surf_usm_h%frac(1,m) * &
2441	surf_usm_h%emissivity(1,m) + &
2442	surf_usm_h%frac(2,m) * &
2443	surf_usm_h%emissivity(2,m)
2444	ENDDO
2445	DO l = 0, 3
2446	DO m = 1, surf_usm_v(l)%ns
2447	i = surf_usm_v(l)%i(m)
2448	j = surf_usm_v(l)%j(m)
2449	k = surf_usm_v(l)%k(m)
2450	temp_pf(k,j,i) = surf_usm_v(l)%frac(0,m) * &
2451	surf_usm_v(l)%emissivity(0,m) + &
2452	surf_usm_v(l)%frac(1,m) * &
2453	surf_usm_v(l)%emissivity(1,m) + &
2454	surf_usm_v(l)%frac(2,m) * &
2455	surf_usm_v(l)%emissivity(2,m)
2456	ENDDO
2457	ENDDO
2458
2459	CASE ( 'usm_surfwintrans' )
2460	!-- transmissivity window tiles
2461	DO m = 1, surf_usm_h%ns
2462	i = surf_usm_h%i(m)
2463	j = surf_usm_h%j(m)
2464	k = surf_usm_h%k(m)
2465	temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
2466	ENDDO
2467	DO l = 0, 3
2468	DO m = 1, surf_usm_v(l)%ns
2469	i = surf_usm_v(l)%i(m)
2470	j = surf_usm_v(l)%j(m)
2471	k = surf_usm_v(l)%k(m)
2472	temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
2473	ENDDO
2474
2475	ENDDO
2476
2477	!
2478	!-- Not adjusted so far
2479	CASE ( 'usm_svf', 'usm_dif' )
2480	!-- shape view factors or iradiance factors to selected surface
2481	IF ( TRIM(var)=='usm_svf' ) THEN
2482	k = 1
2483	ELSE
2484	k = 2
2485	ENDIF
2486	DO isvf = 1, nsvfl
2487	isurflt = svfsurf(1, isvf)
2488	isurfs = svfsurf(2, isvf)
2489
2490	IF ( surf(ix,isurfs) == is .AND. surf(iy,isurfs) == js .AND. &
2491	surf(iz,isurfs) == ks .AND. surf(id,isurfs) == ids ) THEN
2492	!-- correct source surface
2493	temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
2494	ENDIF
2495	ENDDO
2496
2497	CASE ( 'usm_rad_net' )
2498	!-- array of complete radiation balance
2499	IF ( av == 0 ) THEN
2500	DO m = 1, surf_usm_h%ns
2501	i = surf_usm_h%i(m)
2502	j = surf_usm_h%j(m)
2503	k = surf_usm_h%k(m)
2504	temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
2505	ENDDO
2506	DO l = 0, 3
2507	DO m = 1, surf_usm_v(l)%ns
2508	i = surf_usm_v(l)%i(m)
2509	j = surf_usm_v(l)%j(m)
2510	k = surf_usm_v(l)%k(m)
2511	temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
2512	ENDDO
2513	ENDDO
2514	ELSE
2515	DO m = 1, surf_usm_h%ns
2516	i = surf_usm_h%i(m)
2517	j = surf_usm_h%j(m)
2518	k = surf_usm_h%k(m)
2519	temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
2520	ENDDO
2521	DO l = 0, 3
2522	DO m = 1, surf_usm_v(l)%ns
2523	i = surf_usm_v(l)%i(m)
2524	j = surf_usm_v(l)%j(m)
2525	k = surf_usm_v(l)%k(m)
2526	temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
2527	ENDDO
2528	ENDDO
2529	ENDIF
2530
2531	CASE ( 'usm_rad_insw' )
2532	!-- array of sw radiation falling to surface after i-th reflection
2533	DO isurf = dirstart(ids), dirend(ids)
2534	IF ( surfl(id,isurf) == ids ) THEN
2535	IF ( av == 0 ) THEN
2536	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
2537	ELSE
2538	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
2539	ENDIF
2540	ENDIF
2541	ENDDO
2542
2543	CASE ( 'usm_rad_inlw' )
2544	!-- array of lw radiation falling to surface after i-th reflection
2545	DO isurf = dirstart(ids), dirend(ids)
2546	IF ( surfl(id,isurf) == ids ) THEN
2547	IF ( av == 0 ) THEN
2548	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
2549	ELSE
2550	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
2551	ENDIF
2552	ENDIF
2553	ENDDO
2554
2555	CASE ( 'usm_rad_inswdir' )
2556	!-- array of direct sw radiation falling to surface from sun
2557	DO isurf = dirstart(ids), dirend(ids)
2558	IF ( surfl(id,isurf) == ids ) THEN
2559	IF ( av == 0 ) THEN
2560	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
2561	ELSE
2562	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
2563	ENDIF
2564	ENDIF
2565	ENDDO
2566
2567	CASE ( 'usm_rad_inswdif' )
2568	!-- array of difusion sw radiation falling to surface from sky and borders of the domain
2569	DO isurf = dirstart(ids), dirend(ids)
2570	IF ( surfl(id,isurf) == ids ) THEN
2571	IF ( av == 0 ) THEN
2572	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
2573	ELSE
2574	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
2575	ENDIF
2576	ENDIF
2577	ENDDO
2578
2579	CASE ( 'usm_rad_inswref' )
2580	!-- array of sw radiation falling to surface from reflections
2581	DO isurf = dirstart(ids), dirend(ids)
2582	IF ( surfl(id,isurf) == ids ) THEN
2583	IF ( av == 0 ) THEN
2584	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
2585	surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
2586	ELSE
2587	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
2588	ENDIF
2589	ENDIF
2590	ENDDO
2591
2592	CASE ( 'usm_rad_inlwref' )
2593	!-- array of lw radiation falling to surface from reflections
2594	DO isurf = dirstart(ids), dirend(ids)
2595	IF ( surfl(id,isurf) == ids ) THEN
2596	IF ( av == 0 ) THEN
2597	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
2598	ELSE
2599	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
2600	ENDIF
2601	ENDIF
2602	ENDDO
2603
2604	CASE ( 'usm_rad_outsw' )
2605	!-- array of sw radiation emitted from surface after i-th reflection
2606	DO isurf = dirstart(ids), dirend(ids)
2607	IF ( surfl(id,isurf) == ids ) THEN
2608	IF ( av == 0 ) THEN
2609	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
2610	ELSE
2611	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
2612	ENDIF
2613	ENDIF
2614	ENDDO
2615
2616	CASE ( 'usm_rad_outlw' )
2617	!-- array of lw radiation emitted from surface after i-th reflection
2618	DO isurf = dirstart(ids), dirend(ids)
2619	IF ( surfl(id,isurf) == ids ) THEN
2620	IF ( av == 0 ) THEN
2621	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
2622	ELSE
2623	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
2624	ENDIF
2625	ENDIF
2626	ENDDO
2627
2628	CASE ( 'usm_rad_ressw' )
2629	!-- average of array of residua of sw radiation absorbed in surface after last reflection
2630	DO isurf = dirstart(ids), dirend(ids)
2631	IF ( surfl(id,isurf) == ids ) THEN
2632	IF ( av == 0 ) THEN
2633	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
2634	ELSE
2635	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
2636	ENDIF
2637	ENDIF
2638	ENDDO
2639
2640	CASE ( 'usm_rad_reslw' )
2641	!-- average of array of residua of lw radiation absorbed in surface after last reflection
2642	DO isurf = dirstart(ids), dirend(ids)
2643	IF ( surfl(id,isurf) == ids ) THEN
2644	IF ( av == 0 ) THEN
2645	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
2646	ELSE
2647	temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
2648	ENDIF
2649	ENDIF
2650	ENDDO
2651
2652	CASE ( 'usm_rad_hf' )
2653	!-- array of heat flux from radiation for surfaces after all reflections
2654	IF ( av == 0 ) THEN
2655	DO m = 1, surf_usm_h%ns
2656	i = surf_usm_h%i(m)
2657	j = surf_usm_h%j(m)
2658	k = surf_usm_h%k(m)
2659	temp_pf(k,j,i) = surf_usm_h%surfhf(m)
2660	ENDDO
2661	DO l = 0, 3
2662	DO m = 1, surf_usm_v(l)%ns
2663	i = surf_usm_v(l)%i(m)
2664	j = surf_usm_v(l)%j(m)
2665	k = surf_usm_v(l)%k(m)
2666	temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
2667	ENDDO
2668	ENDDO
2669	ELSE
2670	DO m = 1, surf_usm_h%ns
2671	i = surf_usm_h%i(m)
2672	j = surf_usm_h%j(m)
2673	k = surf_usm_h%k(m)
2674	temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
2675	ENDDO
2676	DO l = 0, 3
2677	DO m = 1, surf_usm_v(l)%ns
2678	i = surf_usm_v(l)%i(m)
2679	j = surf_usm_v(l)%j(m)
2680	k = surf_usm_v(l)%k(m)
2681	temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
2682	ENDDO
2683	ENDDO
2684	ENDIF
2685
2686	CASE ( 'usm_wshf' )
2687	!-- array of sensible heat flux from surfaces
2688	IF ( av == 0 ) THEN
2689	DO m = 1, surf_usm_h%ns
2690	i = surf_usm_h%i(m)
2691	j = surf_usm_h%j(m)
2692	k = surf_usm_h%k(m)
2693	temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2694	ENDDO
2695	DO l = 0, 3
2696	DO m = 1, surf_usm_v(l)%ns
2697	i = surf_usm_v(l)%i(m)
2698	j = surf_usm_v(l)%j(m)
2699	k = surf_usm_v(l)%k(m)
2700	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2701	ENDDO
2702	ENDDO
2703	ELSE
2704	DO m = 1, surf_usm_h%ns
2705	i = surf_usm_h%i(m)
2706	j = surf_usm_h%j(m)
2707	k = surf_usm_h%k(m)
2708	temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2709	ENDDO
2710	DO l = 0, 3
2711	DO m = 1, surf_usm_v(l)%ns
2712	i = surf_usm_v(l)%i(m)
2713	j = surf_usm_v(l)%j(m)
2714	k = surf_usm_v(l)%k(m)
2715	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2716	ENDDO
2717	ENDDO
2718	ENDIF
2719
2720
2721	CASE ( 'usm_wghf' )
2722	!-- array of heat flux from ground (land, wall, roof)
2723	IF ( av == 0 ) THEN
2724	DO m = 1, surf_usm_h%ns
2725	i = surf_usm_h%i(m)
2726	j = surf_usm_h%j(m)
2727	k = surf_usm_h%k(m)
2728	temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2729	ENDDO
2730	DO l = 0, 3
2731	DO m = 1, surf_usm_v(l)%ns
2732	i = surf_usm_v(l)%i(m)
2733	j = surf_usm_v(l)%j(m)
2734	k = surf_usm_v(l)%k(m)
2735	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2736	ENDDO
2737	ENDDO
2738	ELSE
2739	DO m = 1, surf_usm_h%ns
2740	i = surf_usm_h%i(m)
2741	j = surf_usm_h%j(m)
2742	k = surf_usm_h%k(m)
2743	temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2744	ENDDO
2745	DO l = 0, 3
2746	DO m = 1, surf_usm_v(l)%ns
2747	i = surf_usm_v(l)%i(m)
2748	j = surf_usm_v(l)%j(m)
2749	k = surf_usm_v(l)%k(m)
2750	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2751	ENDDO
2752	ENDDO
2753	ENDIF
2754
2755	CASE ( 'usm_wghf_window' )
2756	!-- array of heat flux from window ground (land, wall, roof)
2757
2758	IF ( av == 0 ) THEN
2759	DO m = 1, surf_usm_h%ns
2760	i = surf_usm_h%i(m)
2761	j = surf_usm_h%j(m)
2762	k = surf_usm_h%k(m)
2763	temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
2764	ENDDO
2765	DO l = 0, 3
2766	DO m = 1, surf_usm_v(l)%ns
2767	i = surf_usm_v(l)%i(m)
2768	j = surf_usm_v(l)%j(m)
2769	k = surf_usm_v(l)%k(m)
2770	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m)
2771	ENDDO
2772	ENDDO
2773	ELSE
2774	DO m = 1, surf_usm_h%ns
2775	i = surf_usm_h%i(m)
2776	j = surf_usm_h%j(m)
2777	k = surf_usm_h%k(m)
2778	temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
2779	ENDDO
2780	DO l = 0, 3
2781	DO m = 1, surf_usm_v(l)%ns
2782	i = surf_usm_v(l)%i(m)
2783	j = surf_usm_v(l)%j(m)
2784	k = surf_usm_v(l)%k(m)
2785	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
2786	ENDDO
2787	ENDDO
2788	ENDIF
2789
2790	CASE ( 'usm_wghf_green' )
2791	!-- array of heat flux from green ground (land, wall, roof)
2792
2793	IF ( av == 0 ) THEN
2794	DO m = 1, surf_usm_h%ns
2795	i = surf_usm_h%i(m)
2796	j = surf_usm_h%j(m)
2797	k = surf_usm_h%k(m)
2798	temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
2799	ENDDO
2800	DO l = 0, 3
2801	DO m = 1, surf_usm_v(l)%ns
2802	i = surf_usm_v(l)%i(m)
2803	j = surf_usm_v(l)%j(m)
2804	k = surf_usm_v(l)%k(m)
2805	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
2806	ENDDO
2807	ENDDO
2808	ELSE
2809	DO m = 1, surf_usm_h%ns
2810	i = surf_usm_h%i(m)
2811	j = surf_usm_h%j(m)
2812	k = surf_usm_h%k(m)
2813	temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
2814	ENDDO
2815	DO l = 0, 3
2816	DO m = 1, surf_usm_v(l)%ns
2817	i = surf_usm_v(l)%i(m)
2818	j = surf_usm_v(l)%j(m)
2819	k = surf_usm_v(l)%k(m)
2820	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m)
2821	ENDDO
2822	ENDDO
2823	ENDIF
2824
2825	CASE ( 'usm_iwghf' )
2826	!-- array of heat flux from indoor ground (land, wall, roof)
2827	IF ( av == 0 ) THEN
2828	DO m = 1, surf_usm_h%ns
2829	i = surf_usm_h%i(m)
2830	j = surf_usm_h%j(m)
2831	k = surf_usm_h%k(m)
2832	temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
2833	ENDDO
2834	DO l = 0, 3
2835	DO m = 1, surf_usm_v(l)%ns
2836	i = surf_usm_v(l)%i(m)
2837	j = surf_usm_v(l)%j(m)
2838	k = surf_usm_v(l)%k(m)
2839	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m)
2840	ENDDO
2841	ENDDO
2842	ELSE
2843	DO m = 1, surf_usm_h%ns
2844	i = surf_usm_h%i(m)
2845	j = surf_usm_h%j(m)
2846	k = surf_usm_h%k(m)
2847	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
2848	ENDDO
2849	DO l = 0, 3
2850	DO m = 1, surf_usm_v(l)%ns
2851	i = surf_usm_v(l)%i(m)
2852	j = surf_usm_v(l)%j(m)
2853	k = surf_usm_v(l)%k(m)
2854	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m)
2855	ENDDO
2856	ENDDO
2857	ENDIF
2858
2859	CASE ( 'usm_iwghf_window' )
2860	!-- array of heat flux from indoor window ground (land, wall, roof)
2861
2862	IF ( av == 0 ) THEN
2863	DO m = 1, surf_usm_h%ns
2864	i = surf_usm_h%i(m)
2865	j = surf_usm_h%j(m)
2866	k = surf_usm_h%k(m)
2867	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
2868	ENDDO
2869	DO l = 0, 3
2870	DO m = 1, surf_usm_v(l)%ns
2871	i = surf_usm_v(l)%i(m)
2872	j = surf_usm_v(l)%j(m)
2873	k = surf_usm_v(l)%k(m)
2874	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m)
2875	ENDDO
2876	ENDDO
2877	ELSE
2878	DO m = 1, surf_usm_h%ns
2879	i = surf_usm_h%i(m)
2880	j = surf_usm_h%j(m)
2881	k = surf_usm_h%k(m)
2882	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
2883	ENDDO
2884	DO l = 0, 3
2885	DO m = 1, surf_usm_v(l)%ns
2886	i = surf_usm_v(l)%i(m)
2887	j = surf_usm_v(l)%j(m)
2888	k = surf_usm_v(l)%k(m)
2889	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m)
2890	ENDDO
2891	ENDDO
2892	ENDIF
2893
2894	CASE ( 'usm_t_surf' )
2895	!-- surface temperature for surfaces
2896	IF ( av == 0 ) THEN
2897	DO m = 1, surf_usm_h%ns
2898	i = surf_usm_h%i(m)
2899	j = surf_usm_h%j(m)
2900	k = surf_usm_h%k(m)
2901	temp_pf(k,j,i) = t_surf_h(m)
2902	ENDDO
2903	DO l = 0, 3
2904	DO m = 1, surf_usm_v(l)%ns
2905	i = surf_usm_v(l)%i(m)
2906	j = surf_usm_v(l)%j(m)
2907	k = surf_usm_v(l)%k(m)
2908	temp_pf(k,j,i) = t_surf_v(l)%t(m)
2909	ENDDO
2910	ENDDO
2911	ELSE
2912	DO m = 1, surf_usm_h%ns
2913	i = surf_usm_h%i(m)
2914	j = surf_usm_h%j(m)
2915	k = surf_usm_h%k(m)
2916	temp_pf(k,j,i) = surf_usm_h%t_surf_av(m)
2917	ENDDO
2918	DO l = 0, 3
2919	DO m = 1, surf_usm_v(l)%ns
2920	i = surf_usm_v(l)%i(m)
2921	j = surf_usm_v(l)%j(m)
2922	k = surf_usm_v(l)%k(m)
2923	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_av(m)
2924	ENDDO
2925	ENDDO
2926	ENDIF
2927
2928	CASE ( 'usm_t_surf_window' )
2929	!-- surface temperature for window surfaces
2930
2931	IF ( av == 0 ) THEN
2932	DO m = 1, surf_usm_h%ns
2933	i = surf_usm_h%i(m)
2934	j = surf_usm_h%j(m)
2935	k = surf_usm_h%k(m)
2936	temp_pf(k,j,i) = t_surf_window_h(m)
2937	ENDDO
2938	DO l = 0, 3
2939	DO m = 1, surf_usm_v(l)%ns
2940	i = surf_usm_v(l)%i(m)
2941	j = surf_usm_v(l)%j(m)
2942	k = surf_usm_v(l)%k(m)
2943	temp_pf(k,j,i) = t_surf_window_v(l)%t(m)
2944	ENDDO
2945	ENDDO
2946
2947	ELSE
2948	DO m = 1, surf_usm_h%ns
2949	i = surf_usm_h%i(m)
2950	j = surf_usm_h%j(m)
2951	k = surf_usm_h%k(m)
2952	temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
2953	ENDDO
2954	DO l = 0, 3
2955	DO m = 1, surf_usm_v(l)%ns
2956	i = surf_usm_v(l)%i(m)
2957	j = surf_usm_v(l)%j(m)
2958	k = surf_usm_v(l)%k(m)
2959	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_window_av(m)
2960	ENDDO
2961
2962	ENDDO
2963
2964	ENDIF
2965
2966	CASE ( 'usm_t_surf_green' )
2967	!-- surface temperature for green surfaces
2968
2969	IF ( av == 0 ) THEN
2970	DO m = 1, surf_usm_h%ns
2971	i = surf_usm_h%i(m)
2972	j = surf_usm_h%j(m)
2973	k = surf_usm_h%k(m)
2974	temp_pf(k,j,i) = t_surf_green_h(m)
2975	ENDDO
2976	DO l = 0, 3
2977	DO m = 1, surf_usm_v(l)%ns
2978	i = surf_usm_v(l)%i(m)
2979	j = surf_usm_v(l)%j(m)
2980	k = surf_usm_v(l)%k(m)
2981	temp_pf(k,j,i) = t_surf_green_v(l)%t(m)
2982	ENDDO
2983	ENDDO
2984
2985	ELSE
2986	DO m = 1, surf_usm_h%ns
2987	i = surf_usm_h%i(m)
2988	j = surf_usm_h%j(m)
2989	k = surf_usm_h%k(m)
2990	temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
2991	ENDDO
2992	DO l = 0, 3
2993	DO m = 1, surf_usm_v(l)%ns
2994	i = surf_usm_v(l)%i(m)
2995	j = surf_usm_v(l)%j(m)
2996	k = surf_usm_v(l)%k(m)
2997	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_green_av(m)
2998	ENDDO
2999
3000	ENDDO
3001
3002	ENDIF
3003
3004	CASE ( 'usm_t_surf_10cm' )
3005	!-- near surface temperature for whole surfaces
3006
3007	IF ( av == 0 ) THEN
3008	DO m = 1, surf_usm_h%ns
3009	i = surf_usm_h%i(m)
3010	j = surf_usm_h%j(m)
3011	k = surf_usm_h%k(m)
3012	temp_pf(k,j,i) = t_surf_10cm_h(m)
3013	ENDDO
3014	DO l = 0, 3
3015	DO m = 1, surf_usm_v(l)%ns
3016	i = surf_usm_v(l)%i(m)
3017	j = surf_usm_v(l)%j(m)
3018	k = surf_usm_v(l)%k(m)
3019	temp_pf(k,j,i) = t_surf_10cm_v(l)%t(m)
3020	ENDDO
3021	ENDDO
3022
3023	ELSE
3024	DO m = 1, surf_usm_h%ns
3025	i = surf_usm_h%i(m)
3026	j = surf_usm_h%j(m)
3027	k = surf_usm_h%k(m)
3028	temp_pf(k,j,i) = surf_usm_h%t_surf_10cm_av(m)
3029	ENDDO
3030	DO l = 0, 3
3031	DO m = 1, surf_usm_v(l)%ns
3032	i = surf_usm_v(l)%i(m)
3033	j = surf_usm_v(l)%j(m)
3034	k = surf_usm_v(l)%k(m)
3035	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_10cm_av(m)
3036	ENDDO
3037
3038	ENDDO
3039
3040	ENDIF
3041
3042
3043	CASE ( 'usm_t_wall' )
3044	!-- wall temperature for iwl layer of walls and land
3045	IF ( av == 0 ) THEN
3046	DO m = 1, surf_usm_h%ns
3047	i = surf_usm_h%i(m)
3048	j = surf_usm_h%j(m)
3049	k = surf_usm_h%k(m)
3050	temp_pf(k,j,i) = t_wall_h(iwl,m)
3051	ENDDO
3052	DO l = 0, 3
3053	DO m = 1, surf_usm_v(l)%ns
3054	i = surf_usm_v(l)%i(m)
3055	j = surf_usm_v(l)%j(m)
3056	k = surf_usm_v(l)%k(m)
3057	temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
3058	ENDDO
3059	ENDDO
3060	ELSE
3061	DO m = 1, surf_usm_h%ns
3062	i = surf_usm_h%i(m)
3063	j = surf_usm_h%j(m)
3064	k = surf_usm_h%k(m)
3065	temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
3066	ENDDO
3067	DO l = 0, 3
3068	DO m = 1, surf_usm_v(l)%ns
3069	i = surf_usm_v(l)%i(m)
3070	j = surf_usm_v(l)%j(m)
3071	k = surf_usm_v(l)%k(m)
3072	temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
3073	ENDDO
3074	ENDDO
3075	ENDIF
3076
3077	CASE ( 'usm_t_window' )
3078	!-- window temperature for iwl layer of walls and land
3079	IF ( av == 0 ) THEN
3080	DO m = 1, surf_usm_h%ns
3081	i = surf_usm_h%i(m)
3082	j = surf_usm_h%j(m)
3083	k = surf_usm_h%k(m)
3084	temp_pf(k,j,i) = t_window_h(iwl,m)
3085	ENDDO
3086	DO l = 0, 3
3087	DO m = 1, surf_usm_v(l)%ns
3088	i = surf_usm_v(l)%i(m)
3089	j = surf_usm_v(l)%j(m)
3090	k = surf_usm_v(l)%k(m)
3091	temp_pf(k,j,i) = t_window_v(l)%t(iwl,m)
3092	ENDDO
3093	ENDDO
3094	ELSE
3095	DO m = 1, surf_usm_h%ns
3096	i = surf_usm_h%i(m)
3097	j = surf_usm_h%j(m)
3098	k = surf_usm_h%k(m)
3099	temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
3100	ENDDO
3101	DO l = 0, 3
3102	DO m = 1, surf_usm_v(l)%ns
3103	i = surf_usm_v(l)%i(m)
3104	j = surf_usm_v(l)%j(m)
3105	k = surf_usm_v(l)%k(m)
3106	temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m)
3107	ENDDO
3108	ENDDO
3109	ENDIF
3110
3111	CASE ( 'usm_t_green' )
3112	!-- green temperature for iwl layer of walls and land
3113	IF ( av == 0 ) THEN
3114	DO m = 1, surf_usm_h%ns
3115	i = surf_usm_h%i(m)
3116	j = surf_usm_h%j(m)
3117	k = surf_usm_h%k(m)
3118	temp_pf(k,j,i) = t_green_h(iwl,m)
3119	ENDDO
3120	DO l = 0, 3
3121	DO m = 1, surf_usm_v(l)%ns
3122	i = surf_usm_v(l)%i(m)
3123	j = surf_usm_v(l)%j(m)
3124	k = surf_usm_v(l)%k(m)
3125	temp_pf(k,j,i) = t_green_v(l)%t(iwl,m)
3126	ENDDO
3127	ENDDO
3128	ELSE
3129	DO m = 1, surf_usm_h%ns
3130	i = surf_usm_h%i(m)
3131	j = surf_usm_h%j(m)
3132	k = surf_usm_h%k(m)
3133	temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
3134	ENDDO
3135	DO l = 0, 3
3136	DO m = 1, surf_usm_v(l)%ns
3137	i = surf_usm_v(l)%i(m)
3138	j = surf_usm_v(l)%j(m)
3139	k = surf_usm_v(l)%k(m)
3140	temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m)
3141	ENDDO
3142	ENDDO
3143	ENDIF
3144
3145
3146	CASE DEFAULT
3147	found = .FALSE.
3148
3149	END SELECT
3150
3151	!
3152	!-- Rearrange dimensions for NetCDF output
3153	DO j = nys, nyn
3154	DO i = nxl, nxr
3155	DO k = nzb_do, nzt_do
3156	local_pf(i,j,k) = temp_pf(k,j,i)
3157	ENDDO
3158	ENDDO
3159	ENDDO
3160
3161	END SUBROUTINE usm_data_output_3d
3162
3163
3164	!------------------------------------------------------------------------------!
3165	!
3166	! Description:
3167	! ------------
3168	!> Soubroutine defines appropriate grid for netcdf variables.
3169	!> It is called out from subroutine netcdf.
3170	!------------------------------------------------------------------------------!
3171	SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
3172
3173	IMPLICIT NONE
3174
3175	CHARACTER (len=*), INTENT(IN) :: variable !<
3176	LOGICAL, INTENT(OUT) :: found !<
3177	CHARACTER (len=*), INTENT(OUT) :: grid_x !<
3178	CHARACTER (len=*), INTENT(OUT) :: grid_y !<
3179	CHARACTER (len=*), INTENT(OUT) :: grid_z !<
3180
3181	CHARACTER (len=varnamelength) :: var
3182
3183	var = TRIM(variable)
3184	IF ( var(1:12) == 'usm_rad_net_' .OR. var(1:13) == 'usm_rad_insw_' .OR. &
3185	var(1:13) == 'usm_rad_inlw_' .OR. var(1:16) == 'usm_rad_inswdir_' .OR. &
3186	var(1:16) == 'usm_rad_inswdif_' .OR. var(1:16) == 'usm_rad_inswref_' .OR. &
3187	var(1:16) == 'usm_rad_inlwdif_' .OR. var(1:16) == 'usm_rad_inlwref_' .OR. &
3188	var(1:14) == 'usm_rad_outsw_' .OR. var(1:14) == 'usm_rad_outlw_' .OR. &
3189	var(1:14) == 'usm_rad_ressw_' .OR. var(1:14) == 'usm_rad_reslw_' .OR. &
3190	var(1:11) == 'usm_rad_hf_' .OR. &
3191	var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
3192	var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. &
3193	var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. &
3194	var(1:10) == 'usm_t_surf' .OR. var(1:10) == 'usm_t_wall' .OR. &
3195	var(1:17) == 'usm_t_surf_window' .OR. var(1:12) == 'usm_t_window' .OR. &
3196	var(1:16) == 'usm_t_surf_green' .OR. &
3197	var(1:15) == 'usm_t_surf_10cm' .OR. &
3198	var(1:9) == 'usm_surfz' .OR. var(1:7) == 'usm_svf' .OR. &
3199	var(1:7) == 'usm_dif' .OR. var(1:11) == 'usm_surfcat' .OR. &
3200	var(1:11) == 'usm_surfalb' .OR. var(1:12) == 'usm_surfemis' .OR. &
3201	var(1:16) == 'usm_surfwintrans' ) THEN
3202
3203	found = .TRUE.
3204	grid_x = 'x'
3205	grid_y = 'y'
3206	grid_z = 'zu'
3207	ELSE
3208	found = .FALSE.
3209	grid_x = 'none'
3210	grid_y = 'none'
3211	grid_z = 'none'
3212	ENDIF
3213
3214	END SUBROUTINE usm_define_netcdf_grid
3215
3216
3217	!------------------------------------------------------------------------------!
3218	! Description:
3219	! ------------
3220	!> Initialization of the wall surface model
3221	!------------------------------------------------------------------------------!
3222	SUBROUTINE usm_init_material_model
3223
3224	IMPLICIT NONE
3225
3226	INTEGER(iwp) :: k, l, m !< running indices
3227
3228	CALL location_message( ' initialization of wall surface model', .TRUE. )
3229
3230	!-- Calculate wall grid spacings.
3231	!-- Temperature is defined at the center of the wall layers,
3232	!-- whereas gradients/fluxes are defined at the edges (_stag)
3233	!-- apply for all particular surface grids. First for horizontal surfaces
3234	DO m = 1, surf_usm_h%ns
3235
3236	surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
3237	DO k = nzb_wall+1, nzt_wall
3238	surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) - &
3239	surf_usm_h%zw(k-1,m)
3240	ENDDO
3241	surf_usm_h%dz_window(nzb_wall,m) = surf_usm_h%zw_window(nzb_wall,m)
3242	DO k = nzb_wall+1, nzt_wall
3243	surf_usm_h%dz_window(k,m) = surf_usm_h%zw_window(k,m) - &
3244	surf_usm_h%zw_window(k-1,m)
3245	ENDDO
3246	surf_usm_h%dz_green(nzb_wall,m) = surf_usm_h%zw_green(nzb_wall,m)
3247	DO k = nzb_wall+1, nzt_wall
3248	surf_usm_h%dz_green(k,m) = surf_usm_h%zw_green(k,m) - &
3249	surf_usm_h%zw_green(k-1,m)
3250	ENDDO
3251
3252	surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
3253
3254	DO k = nzb_wall, nzt_wall-1
3255	surf_usm_h%dz_wall_stag(k,m) = 0.5 * ( &
3256	surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
3257	ENDDO
3258	surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
3259
3260	surf_usm_h%dz_window(nzt_wall+1,m) = surf_usm_h%dz_window(nzt_wall,m)
3261
3262	DO k = nzb_wall, nzt_wall-1
3263	surf_usm_h%dz_window_stag(k,m) = 0.5 * ( &
3264	surf_usm_h%dz_window(k+1,m) + surf_usm_h%dz_window(k,m) )
3265	ENDDO
3266	surf_usm_h%dz_window_stag(nzt_wall,m) = surf_usm_h%dz_window(nzt_wall,m)
3267
3268	surf_usm_h%dz_green(nzt_wall+1,m) = surf_usm_h%dz_green(nzt_wall,m)
3269
3270	DO k = nzb_wall, nzt_wall-1
3271	surf_usm_h%dz_green_stag(k,m) = 0.5 * ( &
3272	surf_usm_h%dz_green(k+1,m) + surf_usm_h%dz_green(k,m) )
3273	ENDDO
3274	surf_usm_h%dz_green_stag(nzt_wall,m) = surf_usm_h%dz_green(nzt_wall,m)
3275	ENDDO
3276	surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall
3277	surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag
3278	surf_usm_h%ddz_window = 1.0_wp / surf_usm_h%dz_window
3279	surf_usm_h%ddz_window_stag = 1.0_wp / surf_usm_h%dz_window_stag
3280	surf_usm_h%ddz_green = 1.0_wp / surf_usm_h%dz_green
3281	surf_usm_h%ddz_green_stag = 1.0_wp / surf_usm_h%dz_green_stag
3282	!
3283	!-- For vertical surfaces
3284	DO l = 0, 3
3285	DO m = 1, surf_usm_v(l)%ns
3286	surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
3287	DO k = nzb_wall+1, nzt_wall
3288	surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - &
3289	surf_usm_v(l)%zw(k-1,m)
3290	ENDDO
3291	surf_usm_v(l)%dz_window(nzb_wall,m) = surf_usm_v(l)%zw_window(nzb_wall,m)
3292	DO k = nzb_wall+1, nzt_wall
3293	surf_usm_v(l)%dz_window(k,m) = surf_usm_v(l)%zw_window(k,m) - &
3294	surf_usm_v(l)%zw_window(k-1,m)
3295	ENDDO
3296	surf_usm_v(l)%dz_green(nzb_wall,m) = surf_usm_v(l)%zw_green(nzb_wall,m)
3297	DO k = nzb_wall+1, nzt_wall
3298	surf_usm_v(l)%dz_green(k,m) = surf_usm_v(l)%zw_green(k,m) - &
3299	surf_usm_v(l)%zw_green(k-1,m)
3300	ENDDO
3301
3302	surf_usm_v(l)%dz_wall(nzt_wall+1,m) = &
3303	surf_usm_v(l)%dz_wall(nzt_wall,m)
3304
3305	DO k = nzb_wall, nzt_wall-1
3306	surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * ( &
3307	surf_usm_v(l)%dz_wall(k+1,m) + &
3308	surf_usm_v(l)%dz_wall(k,m) )
3309	ENDDO
3310	surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = &
3311	surf_usm_v(l)%dz_wall(nzt_wall,m)
3312	surf_usm_v(l)%dz_window(nzt_wall+1,m) = &
3313	surf_usm_v(l)%dz_window(nzt_wall,m)
3314
3315	DO k = nzb_wall, nzt_wall-1
3316	surf_usm_v(l)%dz_window_stag(k,m) = 0.5 * ( &
3317	surf_usm_v(l)%dz_window(k+1,m) + &
3318	surf_usm_v(l)%dz_window(k,m) )
3319	ENDDO
3320	surf_usm_v(l)%dz_window_stag(nzt_wall,m) = &
3321	surf_usm_v(l)%dz_window(nzt_wall,m)
3322	surf_usm_v(l)%dz_green(nzt_wall+1,m) = &
3323	surf_usm_v(l)%dz_green(nzt_wall,m)
3324
3325	DO k = nzb_wall, nzt_wall-1
3326	surf_usm_v(l)%dz_green_stag(k,m) = 0.5 * ( &
3327	surf_usm_v(l)%dz_green(k+1,m) + &
3328	surf_usm_v(l)%dz_green(k,m) )
3329	ENDDO
3330	surf_usm_v(l)%dz_green_stag(nzt_wall,m) = &
3331	surf_usm_v(l)%dz_green(nzt_wall,m)
3332	ENDDO
3333	surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall
3334	surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag
3335	surf_usm_v(l)%ddz_window = 1.0_wp / surf_usm_v(l)%dz_window
3336	surf_usm_v(l)%ddz_window_stag = 1.0_wp / surf_usm_v(l)%dz_window_stag
3337	surf_usm_v(l)%ddz_green = 1.0_wp / surf_usm_v(l)%dz_green
3338	surf_usm_v(l)%ddz_green_stag = 1.0_wp / surf_usm_v(l)%dz_green_stag
3339	ENDDO
3340
3341
3342	CALL location_message( ' wall structures filed out', .TRUE. )
3343
3344	CALL location_message( ' initialization of wall surface model finished', .TRUE. )
3345
3346	END SUBROUTINE usm_init_material_model
3347
3348
3349	!------------------------------------------------------------------------------!
3350	! Description:
3351	! ------------
3352	!> Initialization of the urban surface model
3353	!------------------------------------------------------------------------------!
3354	SUBROUTINE usm_init_urban_surface
3355
3356	USE arrays_3d, &
3357	ONLY: zw
3358
3359	USE netcdf_data_input_mod, &
3360	ONLY: building_pars_f, building_type_f, terrain_height_f
3361
3362	IMPLICIT NONE
3363
3364	INTEGER(iwp) :: i !< loop index x-dirction
3365	INTEGER(iwp) :: ind_emis_wall !< index in input list for wall emissivity
3366	INTEGER(iwp) :: ind_emis_green !< index in input list for green emissivity
3367	INTEGER(iwp) :: ind_emis_win !< index in input list for window emissivity
3368	INTEGER(iwp) :: ind_green_frac_w !< index in input list for green fraction on wall
3369	INTEGER(iwp) :: ind_green_frac_r !< index in input list for green fraction on roof
3370	INTEGER(iwp) :: ind_hc1 !< index in input list for heat capacity at first wall layer
3371	INTEGER(iwp) :: ind_hc2 !< index in input list for heat capacity at second wall layer
3372	INTEGER(iwp) :: ind_hc3 !< index in input list for heat capacity at third wall layer
3373	INTEGER(iwp) :: ind_lai_r !< index in input list for LAI on roof
3374	INTEGER(iwp) :: ind_lai_w !< index in input list for LAI on wall
3375	INTEGER(iwp) :: ind_tc1 !< index in input list for thermal conductivity at first wall layer
3376	INTEGER(iwp) :: ind_tc2 !< index in input list for thermal conductivity at second wall layer
3377	INTEGER(iwp) :: ind_tc3 !< index in input list for thermal conductivity at third wall layer
3378	INTEGER(iwp) :: ind_trans !< index in input list for window transmissivity
3379	INTEGER(iwp) :: ind_wall_frac !< index in input list for wall fraction
3380	INTEGER(iwp) :: ind_win_frac !< index in input list for window fraction
3381	INTEGER(iwp) :: ind_z0 !< index in input list for z0
3382	INTEGER(iwp) :: ind_z0qh !< index in input list for z0h / z0q
3383	INTEGER(iwp) :: j !< loop index y-dirction
3384	INTEGER(iwp) :: k !< loop index z-dirction
3385	INTEGER(iwp) :: l !< loop index surface orientation
3386	INTEGER(iwp) :: m !< loop index surface element
3387	INTEGER(iwp) :: st !< dummy
3388
3389	REAL(wp) :: c, d, tin, twin, exn
3390	REAL(wp) :: ground_floor_level_l !< local height of ground floor level
3391	REAL(wp) :: z_agl !< height above ground
3392
3393	!
3394	!-- NOPOINTER version not implemented yet
3395	#if defined( __nopointer )
3396	message_string = 'The urban surface module only runs with POINTER version'
3397	CALL message( 'urban_surface_mod', 'PA0452', 1, 2, 0, 6, 0 )
3398	#endif
3399
3400	CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
3401	!-- surface forcing have to be disabled for LSF
3402	!-- in case of enabled urban surface module
3403	IF ( large_scale_forcing ) THEN
3404	lsf_surf = .FALSE.
3405	ENDIF
3406
3407	!
3408	!-- Flag surface elements belonging to the ground floor level. Therefore,
3409	!-- use terrain height array from file, if available. This flag is later used
3410	!-- to control initialization of surface attributes.
3411	surf_usm_h%ground_level = .FALSE.
3412	DO m = 1, surf_usm_h%ns
3413	i = surf_usm_h%i(m)
3414	j = surf_usm_h%j(m)
3415	k = surf_usm_h%k(m)
3416	!
3417	!-- Get local ground level. If no ground level is given in input file,
3418	!-- use default value.
3419	ground_floor_level_l = ground_floor_level
3420	IF ( building_pars_f%from_file ) THEN
3421	IF ( building_pars_f%pars_xy(ind_gflh,j,i) /= &
3422	building_pars_f%fill ) &
3423	ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)
3424	ENDIF
3425	!
3426	!-- Determine height of surface element above ground level
3427	IF ( terrain_height_f%from_file ) THEN
3428	z_agl = zw(k) - terrain_height_f%var(j,i)
3429	ELSE
3430	z_agl = zw(k)
3431	ENDIF
3432	!
3433	!-- Set flag for ground level
3434	IF ( z_agl <= ground_floor_level_l ) &
3435	surf_usm_h%ground_level(m) = .TRUE.
3436	ENDDO
3437
3438	DO l = 0, 3
3439	surf_usm_v(l)%ground_level = .FALSE.
3440	DO m = 1, surf_usm_v(l)%ns
3441	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
3442	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
3443	k = surf_usm_v(l)%k(m)
3444	!
3445	!-- Get local ground level. If no ground level is given in input file,
3446	!-- use default value.
3447	ground_floor_level_l = ground_floor_level
3448	IF ( building_pars_f%from_file ) THEN
3449	IF ( building_pars_f%pars_xy(ind_gflh,j,i) /= &
3450	building_pars_f%fill ) &
3451	ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)
3452	ENDIF
3453	!
3454	!-- Determine height of surface element above ground level. Please
3455	!-- note, height of surface element is determined with respect to
3456	!-- its height of the adjoing atmospheric grid point.
3457	IF ( terrain_height_f%from_file ) THEN
3458	z_agl = zw(k) - terrain_height_f%var(j-surf_usm_v(l)%joff, &
3459	i-surf_usm_v(l)%ioff)
3460	ELSE
3461	z_agl = zw(k)
3462	ENDIF
3463	!
3464	!-- Set flag for ground level
3465	IF ( z_agl <= ground_floor_level_l ) &
3466	surf_usm_v(l)%ground_level(m) = .TRUE.
3467
3468	ENDDO
3469	ENDDO
3470	!
3471	!-- Initialize urban-type surface attribute. According to initialization in
3472	!-- land-surface model, follow a 3-level approach.
3473	!-- Level 1 - initialization via default attributes
3474	DO m = 1, surf_usm_h%ns
3475	!
3476	!-- Now, all horizontal surfaces are roof surfaces (?)
3477	surf_usm_h%isroof_surf(m) = .TRUE.
3478	surf_usm_h%surface_types(m) = roof_category !< default category for root surface
3479	!
3480	!-- In order to distinguish between ground floor level and
3481	!-- above-ground-floor level surfaces, set input indices.
3482	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3483	surf_usm_h%ground_level(m) )
3484	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3485	surf_usm_h%ground_level(m) )
3486	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3487	surf_usm_h%ground_level(m) )
3488	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3489	surf_usm_h%ground_level(m) )
3490	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3491	surf_usm_h%ground_level(m) )
3492	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
3493	surf_usm_h%ground_level(m) )
3494	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
3495	surf_usm_h%ground_level(m) )
3496	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
3497	surf_usm_h%ground_level(m) )
3498	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
3499	surf_usm_h%ground_level(m) )
3500	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
3501	surf_usm_h%ground_level(m) )
3502	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
3503	surf_usm_h%ground_level(m) )
3504	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
3505	surf_usm_h%ground_level(m) )
3506	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
3507	surf_usm_h%ground_level(m) )
3508	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
3509	surf_usm_h%ground_level(m) )
3510	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
3511	surf_usm_h%ground_level(m) )
3512	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
3513	surf_usm_h%ground_level(m) )
3514	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3515	surf_usm_h%ground_level(m) )
3516	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3517	surf_usm_h%ground_level(m) )
3518	!
3519	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3520	surf_usm_h%frac(0,m) = building_pars(ind_wall_frac,building_type)
3521	surf_usm_h%frac(1,m) = building_pars(ind_green_frac_r,building_type)
3522	surf_usm_h%frac(2,m) = building_pars(ind_win_frac,building_type)
3523	surf_usm_h%lai(m) = building_pars(ind_green_frac_r,building_type)
3524
3525	surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,building_type)
3526	surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3527	surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3528	surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
3529	surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1,building_type)
3530	surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
3531	surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3532	surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3533	surf_usm_h%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,building_type)
3534	surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3535	surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3536	surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
3537	surf_usm_h%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,building_type)
3538	surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
3539	surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3540	surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3541	surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,building_type)
3542	surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3543	surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3544	surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
3545	surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,building_type)
3546	surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
3547	surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3548	surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3549
3550	surf_usm_h%target_temp_summer(m) = building_pars(12,building_type)
3551	surf_usm_h%target_temp_winter(m) = building_pars(13,building_type)
3552	!
3553	!-- emissivity of wall-, green- and window fraction
3554	surf_usm_h%emissivity(0,m) = building_pars(ind_emis_wall,building_type)
3555	surf_usm_h%emissivity(1,m) = building_pars(ind_emis_green,building_type)
3556	surf_usm_h%emissivity(2,m) = building_pars(ind_emis_win,building_type)
3557
3558	surf_usm_h%transmissivity(m) = building_pars(ind_trans,building_type)
3559
3560	surf_usm_h%z0(m) = building_pars(ind_z0,building_type)
3561	surf_usm_h%z0h(m) = building_pars(ind_z0qh,building_type)
3562	surf_usm_h%z0q(m) = building_pars(ind_z0qh,building_type)
3563	!
3564	!-- albedo type for wall fraction, green fraction, window fraction
3565	surf_usm_h%albedo_type(0,m) = INT( building_pars(ind_alb_wall,building_type) )
3566	surf_usm_h%albedo_type(1,m) = INT( building_pars(ind_alb_green,building_type) )
3567	surf_usm_h%albedo_type(2,m) = INT( building_pars(ind_alb_win,building_type) )
3568
3569	surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1,building_type)
3570	surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
3571	surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
3572	surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
3573
3574	surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1,building_type)
3575	surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
3576	surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
3577	surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
3578
3579	surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1,building_type)
3580	surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
3581	surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
3582	surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
3583
3584	surf_usm_h%c_surface(m) = building_pars(45,building_type)
3585	surf_usm_h%lambda_surf(m) = building_pars(46,building_type)
3586	surf_usm_h%c_surface_green(m) = building_pars(45,building_type)
3587	surf_usm_h%lambda_surf_green(m) = building_pars(46,building_type)
3588	surf_usm_h%c_surface_window(m) = building_pars(45,building_type)
3589	surf_usm_h%lambda_surf_window(m) = building_pars(46,building_type)
3590
3591	ENDDO
3592
3593	DO l = 0, 3
3594	DO m = 1, surf_usm_v(l)%ns
3595
3596	surf_usm_v(l)%surface_types(m) = wall_category !< default category for root surface
3597	!
3598	!-- In order to distinguish between ground floor level and
3599	!-- above-ground-floor level surfaces, set input indices.
3600	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3601	surf_usm_v(l)%ground_level(m) )
3602	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3603	surf_usm_v(l)%ground_level(m) )
3604	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3605	surf_usm_v(l)%ground_level(m) )
3606	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3607	surf_usm_v(l)%ground_level(m) )
3608	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3609	surf_usm_v(l)%ground_level(m) )
3610	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
3611	surf_usm_v(l)%ground_level(m) )
3612	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
3613	surf_usm_v(l)%ground_level(m) )
3614	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
3615	surf_usm_v(l)%ground_level(m) )
3616	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
3617	surf_usm_v(l)%ground_level(m) )
3618	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
3619	surf_usm_v(l)%ground_level(m) )
3620	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
3621	surf_usm_v(l)%ground_level(m) )
3622	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
3623	surf_usm_v(l)%ground_level(m) )
3624	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
3625	surf_usm_v(l)%ground_level(m) )
3626	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
3627	surf_usm_v(l)%ground_level(m) )
3628	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
3629	surf_usm_v(l)%ground_level(m) )
3630	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
3631	surf_usm_v(l)%ground_level(m) )
3632	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3633	surf_usm_v(l)%ground_level(m) )
3634	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3635	surf_usm_v(l)%ground_level(m) )
3636
3637	!
3638	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3639	surf_usm_v(l)%frac(0,m) = building_pars(ind_wall_frac,building_type)
3640	surf_usm_v(l)%frac(1,m) = building_pars(ind_green_frac_w,building_type)
3641	surf_usm_v(l)%frac(2,m) = building_pars(ind_win_frac,building_type)
3642	surf_usm_v(l)%lai(m) = building_pars(ind_lai_w,building_type)
3643
3644	surf_usm_v(l)%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,building_type)
3645	surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3646	surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3647	surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
3648
3649	surf_usm_v(l)%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,building_type)
3650	surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3651	surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3652	surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
3653
3654	surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,building_type)
3655	surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
3656	surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
3657	surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
3658
3659	surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars(ind_tc1,building_type)
3660	surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
3661	surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3662	surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3663
3664	surf_usm_v(l)%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,building_type)
3665	surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
3666	surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3667	surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3668
3669	surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,building_type)
3670	surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
3671	surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
3672	surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
3673
3674	surf_usm_v(l)%target_temp_summer(m) = building_pars(12,building_type)
3675	surf_usm_v(l)%target_temp_winter(m) = building_pars(13,building_type)
3676	!
3677	!-- emissivity of wall-, green- and window fraction
3678	surf_usm_v(l)%emissivity(0,m) = building_pars(ind_emis_wall,building_type)
3679	surf_usm_v(l)%emissivity(1,m) = building_pars(ind_emis_green,building_type)
3680	surf_usm_v(l)%emissivity(2,m) = building_pars(ind_emis_win,building_type)
3681
3682	surf_usm_v(l)%transmissivity(m) = building_pars(ind_trans,building_type)
3683
3684	surf_usm_v(l)%z0(m) = building_pars(ind_z0,building_type)
3685	surf_usm_v(l)%z0h(m) = building_pars(ind_z0qh,building_type)
3686	surf_usm_v(l)%z0q(m) = building_pars(ind_z0qh,building_type)
3687
3688	surf_usm_v(l)%albedo_type(0,m) = INT( building_pars(ind_alb_wall,building_type) )
3689	surf_usm_v(l)%albedo_type(1,m) = INT( building_pars(ind_alb_green,building_type) )
3690	surf_usm_v(l)%albedo_type(2,m) = INT( building_pars(ind_alb_win,building_type) )
3691
3692	surf_usm_v(l)%zw(nzb_wall,m) = building_pars(ind_thick_1,building_type)
3693	surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
3694	surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
3695	surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
3696
3697	surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars(ind_thick_1,building_type)
3698	surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
3699	surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
3700	surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
3701
3702	surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars(ind_thick_1,building_type)
3703	surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
3704	surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
3705	surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
3706
3707	surf_usm_v(l)%c_surface(m) = building_pars(45,building_type)
3708	surf_usm_v(l)%lambda_surf(m) = building_pars(46,building_type)
3709	surf_usm_v(l)%c_surface_green(m) = building_pars(45,building_type)
3710	surf_usm_v(l)%lambda_surf_green(m) = building_pars(46,building_type)
3711	surf_usm_v(l)%c_surface_window(m) = building_pars(45,building_type)
3712	surf_usm_v(l)%lambda_surf_window(m) = building_pars(46,building_type)
3713
3714	ENDDO
3715	ENDDO
3716	!
3717	!-- Level 2 - initialization via building type read from file
3718	IF ( building_type_f%from_file ) THEN
3719	DO m = 1, surf_usm_h%ns
3720	i = surf_usm_h%i(m)
3721	j = surf_usm_h%j(m)
3722	!
3723	!-- For the moment, limit building type to 6 (to overcome errors in input file).
3724	st = building_type_f%var(j,i)
3725	IF ( st /= building_type_f%fill ) THEN
3726
3727	!
3728	!-- In order to distinguish between ground floor level and
3729	!-- above-ground-floor level surfaces, set input indices.
3730	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3731	surf_usm_h%ground_level(m) )
3732	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3733	surf_usm_h%ground_level(m) )
3734	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3735	surf_usm_h%ground_level(m) )
3736	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3737	surf_usm_h%ground_level(m) )
3738	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3739	surf_usm_h%ground_level(m) )
3740	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
3741	surf_usm_h%ground_level(m) )
3742	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
3743	surf_usm_h%ground_level(m) )
3744	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
3745	surf_usm_h%ground_level(m) )
3746	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
3747	surf_usm_h%ground_level(m) )
3748	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
3749	surf_usm_h%ground_level(m) )
3750	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
3751	surf_usm_h%ground_level(m) )
3752	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
3753	surf_usm_h%ground_level(m) )
3754	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
3755	surf_usm_h%ground_level(m) )
3756	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
3757	surf_usm_h%ground_level(m) )
3758	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
3759	surf_usm_h%ground_level(m) )
3760	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
3761	surf_usm_h%ground_level(m) )
3762	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3763	surf_usm_h%ground_level(m) )
3764	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3765	surf_usm_h%ground_level(m) )
3766
3767	!
3768	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3769	surf_usm_h%frac(0,m) = building_pars(ind_wall_frac,st)
3770	surf_usm_h%frac(1,m) = building_pars(ind_green_frac_r,st)
3771	surf_usm_h%frac(2,m) = building_pars(ind_win_frac,st)
3772	surf_usm_h%lai(m) = building_pars(ind_green_frac_r,st)
3773
3774	surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,st)
3775	surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st)
3776	surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st)
3777	surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st)
3778	surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1,st)
3779	surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st)
3780	surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st)
3781	surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st)
3782
3783	surf_usm_h%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,st)
3784	surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,st)
3785	surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,st)
3786	surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,st)
3787	surf_usm_h%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,st)
3788	surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,st)
3789	surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,st)
3790	surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,st)
3791
3792	surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,st)
3793	surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,st)
3794	surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,st)
3795	surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,st)
3796	surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,st)
3797	surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,st)
3798	surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,st)
3799	surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,st)
3800
3801	surf_usm_h%target_temp_summer(m) = building_pars(12,st)
3802	surf_usm_h%target_temp_winter(m) = building_pars(13,st)
3803	!
3804	!-- emissivity of wall-, green- and window fraction
3805	surf_usm_h%emissivity(0,m) = building_pars(ind_emis_wall,st)
3806	surf_usm_h%emissivity(1,m) = building_pars(ind_emis_green,st)
3807	surf_usm_h%emissivity(2,m) = building_pars(ind_emis_win,st)
3808
3809	surf_usm_h%transmissivity(m) = building_pars(ind_trans,st)
3810
3811	surf_usm_h%z0(m) = building_pars(ind_z0,st)
3812	surf_usm_h%z0h(m) = building_pars(ind_z0qh,st)
3813	surf_usm_h%z0q(m) = building_pars(ind_z0qh,st)
3814	!
3815	!-- albedo type for wall fraction, green fraction, window fraction
3816	surf_usm_h%albedo_type(0,m) = INT( building_pars(ind_alb_wall,st) )
3817	surf_usm_h%albedo_type(1,m) = INT( building_pars(ind_alb_green,st) )
3818	surf_usm_h%albedo_type(2,m) = INT( building_pars(ind_alb_win,st) )
3819
3820	surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1,st)
3821	surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2,st)
3822	surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3,st)
3823	surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4,st)
3824
3825	surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1,st)
3826	surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,st)
3827	surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,st)
3828	surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,st)
3829
3830	surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1,st)
3831	surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,st)
3832	surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,st)
3833	surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,st)
3834
3835	surf_usm_h%c_surface(m) = building_pars(45,st)
3836	surf_usm_h%lambda_surf(m) = building_pars(46,st)
3837	surf_usm_h%c_surface_green(m) = building_pars(45,st)
3838	surf_usm_h%lambda_surf_green(m) = building_pars(46,st)
3839	surf_usm_h%c_surface_window(m) = building_pars(45,st)
3840	surf_usm_h%lambda_surf_window(m) = building_pars(46,st)
3841
3842	ENDIF
3843	ENDDO
3844
3845	DO l = 0, 3
3846	DO m = 1, surf_usm_v(l)%ns
3847	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
3848	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
3849	!
3850	!-- For the moment, limit building type to 6 (to overcome errors in input file).
3851
3852	st = building_type_f%var(j,i)
3853	IF ( st /= building_type_f%fill ) THEN
3854
3855	!
3856	!-- In order to distinguish between ground floor level and
3857	!-- above-ground-floor level surfaces, set input indices.
3858	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3859	surf_usm_v(l)%ground_level(m) )
3860	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3861	surf_usm_v(l)%ground_level(m) )
3862	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3863	surf_usm_v(l)%ground_level(m) )
3864	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3865	surf_usm_v(l)%ground_level(m) )
3866	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3867	surf_usm_v(l)%ground_level(m) )
3868	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
3869	surf_usm_v(l)%ground_level(m) )
3870	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
3871	surf_usm_v(l)%ground_level(m) )
3872	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
3873	surf_usm_v(l)%ground_level(m) )
3874	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
3875	surf_usm_v(l)%ground_level(m) )
3876	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
3877	surf_usm_v(l)%ground_level(m) )
3878	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
3879	surf_usm_v(l)%ground_level(m) )
3880	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
3881	surf_usm_v(l)%ground_level(m) )
3882	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
3883	surf_usm_v(l)%ground_level(m) )
3884	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
3885	surf_usm_v(l)%ground_level(m) )
3886	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
3887	surf_usm_v(l)%ground_level(m) )
3888	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
3889	surf_usm_v(l)%ground_level(m) )
3890	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3891	surf_usm_v(l)%ground_level(m) )
3892	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3893	surf_usm_v(l)%ground_level(m) )
3894
3895	!
3896	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3897	surf_usm_v(l)%frac(0,m) = building_pars(ind_wall_frac,st)
3898	surf_usm_v(l)%frac(1,m) = building_pars(ind_green_frac_w,st)
3899	surf_usm_v(l)%frac(2,m) = building_pars(ind_win_frac,st)
3900	surf_usm_v(l)%lai(m) = building_pars(ind_lai_w,st)
3901
3902	surf_usm_v(l)%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,st)
3903	surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st)
3904	surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st)
3905	surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st)
3906
3907	surf_usm_v(l)%rho_c_green(nzb_wall,m) = building_pars(ind_hc1,st)
3908	surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = building_pars(ind_hc1,st)
3909	surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = building_pars(ind_hc2,st)
3910	surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = building_pars(ind_hc3,st)
3911
3912	surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars(ind_hc1,st)
3913	surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1,st)
3914	surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2,st)
3915	surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3,st)
3916
3917	surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars(ind_tc1,st)
3918	surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st)
3919	surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st)
3920	surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st)
3921
3922	surf_usm_v(l)%lambda_h_green(nzb_wall,m) = building_pars(ind_tc1,st)
3923	surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = building_pars(ind_tc1,st)
3924	surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = building_pars(ind_tc2,st)
3925	surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = building_pars(ind_tc3,st)
3926
3927	surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1,st)
3928	surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1,st)
3929	surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2,st)
3930	surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3,st)
3931
3932	surf_usm_v(l)%target_temp_summer(m) = building_pars(12,st)
3933	surf_usm_v(l)%target_temp_winter(m) = building_pars(13,st)
3934	!
3935	!-- emissivity of wall-, green- and window fraction
3936	surf_usm_v(l)%emissivity(0,m) = building_pars(ind_emis_wall,st)
3937	surf_usm_v(l)%emissivity(1,m) = building_pars(ind_emis_green,st)
3938	surf_usm_v(l)%emissivity(2,m) = building_pars(ind_emis_win,st)
3939
3940	surf_usm_v(l)%transmissivity(m) = building_pars(ind_trans,st)
3941
3942	surf_usm_v(l)%z0(m) = building_pars(ind_z0,st)
3943	surf_usm_v(l)%z0h(m) = building_pars(ind_z0qh,st)
3944	surf_usm_v(l)%z0q(m) = building_pars(ind_z0qh,st)
3945
3946	surf_usm_v(l)%albedo_type(0,m) = INT( building_pars(ind_alb_wall,st) )
3947	surf_usm_v(l)%albedo_type(1,m) = INT( building_pars(ind_alb_green,st) )
3948	surf_usm_v(l)%albedo_type(2,m) = INT( building_pars(ind_alb_win,st) )
3949
3950	surf_usm_v(l)%zw(nzb_wall,m) = building_pars(ind_thick_1,st)
3951	surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars(ind_thick_2,st)
3952	surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars(ind_thick_3,st)
3953	surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars(ind_thick_4,st)
3954
3955	surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars(ind_thick_1,st)
3956	surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,st)
3957	surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,st)
3958	surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,st)
3959
3960	surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars(ind_thick_1,st)
3961	surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2,st)
3962	surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3,st)
3963	surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4,st)
3964
3965	surf_usm_v(l)%c_surface(m) = building_pars(45,st)
3966	surf_usm_v(l)%lambda_surf(m) = building_pars(46,st)
3967	surf_usm_v(l)%c_surface_green(m) = building_pars(45,st)
3968	surf_usm_v(l)%lambda_surf_green(m) = building_pars(46,st)
3969	surf_usm_v(l)%c_surface_window(m) = building_pars(45,st)
3970	surf_usm_v(l)%lambda_surf_window(m) = building_pars(46,st)
3971
3972
3973	ENDIF
3974	ENDDO
3975	ENDDO
3976	ENDIF
3977
3978	!
3979	!-- Level 3 - initialization via building_pars read from file
3980	IF ( building_pars_f%from_file ) THEN
3981	DO m = 1, surf_usm_h%ns
3982	i = surf_usm_h%i(m)
3983	j = surf_usm_h%j(m)
3984
3985	!
3986	!-- In order to distinguish between ground floor level and
3987	!-- above-ground-floor level surfaces, set input indices.
3988	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3989	surf_usm_h%ground_level(m) )
3990	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3991	surf_usm_h%ground_level(m) )
3992	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3993	surf_usm_h%ground_level(m) )
3994	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3995	surf_usm_h%ground_level(m) )
3996	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3997	surf_usm_h%ground_level(m) )
3998	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
3999	surf_usm_h%ground_level(m) )
4000	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
4001	surf_usm_h%ground_level(m) )
4002	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
4003	surf_usm_h%ground_level(m) )
4004	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
4005	surf_usm_h%ground_level(m) )
4006	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
4007	surf_usm_h%ground_level(m) )
4008	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
4009	surf_usm_h%ground_level(m) )
4010	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
4011	surf_usm_h%ground_level(m) )
4012	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
4013	surf_usm_h%ground_level(m) )
4014	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
4015	surf_usm_h%ground_level(m) )
4016	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
4017	surf_usm_h%ground_level(m) )
4018	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
4019	surf_usm_h%ground_level(m) )
4020	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
4021	surf_usm_h%ground_level(m) )
4022	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
4023	surf_usm_h%ground_level(m) )
4024
4025	!
4026	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
4027	IF ( building_pars_f%pars_xy(ind_wall_frac,j,i) /= building_pars_f%fill ) &
4028	surf_usm_h%frac(0,m) = building_pars_f%pars_xy(ind_wall_frac,j,i)
4029	IF ( building_pars_f%pars_xy(ind_green_frac_r,j,i) /= building_pars_f%fill ) &
4030	surf_usm_h%frac(1,m) = building_pars_f%pars_xy(ind_green_frac_r,j,i)
4031	IF ( building_pars_f%pars_xy(ind_win_frac,j,i) /= building_pars_f%fill ) &
4032	surf_usm_h%frac(2,m) = building_pars_f%pars_xy(ind_win_frac,j,i)
4033
4034
4035	IF ( building_pars_f%pars_xy(ind_lai_r,j,i) /= building_pars_f%fill ) &
4036	surf_usm_h%lai(m) = building_pars_f%pars_xy(ind_lai_r,j,i)
4037
4038	IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN
4039	surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4040	surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4041	ENDIF
4042	IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) &
4043	surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i)
4044	IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) &
4045	surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i)
4046	IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN
4047	surf_usm_h%rho_c_green(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4048	surf_usm_h%rho_c_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4049	ENDIF
4050	IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) &
4051	surf_usm_h%rho_c_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i)
4052	IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) &
4053	surf_usm_h%rho_c_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i)
4054	IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= building_pars_f%fill ) THEN
4055	surf_usm_h%rho_c_window(nzb_wall,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4056	surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_hc1,j,i)
4057	ENDIF
4058	IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= building_pars_f%fill ) &
4059	surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_hc2,j,i)
4060	IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= building_pars_f%fill ) &
4061	surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_hc3,j,i)
4062
4063	IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN
4064	surf_usm_h%lambda_h(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i)
4065	surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i)
4066	ENDIF
4067	IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) &
4068	surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i)
4069	IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) &
4070	surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i)
4071	IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN
4072	surf_usm_h%lambda_h_green(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i)
4073	surf_usm_h%lambda_h_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i)
4074	ENDIF
4075	IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) &
4076	surf_usm_h%lambda_h_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i)
4077	IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) &
4078	surf_usm_h%lambda_h_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i)
4079	IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= building_pars_f%fill ) THEN
4080	surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars_f%pars_xy(ind_tc1,j,i)
4081	surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_tc1,j,i)
4082	ENDIF
4083	IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= building_pars_f%fill ) &
4084	surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_tc2,j,i)
4085	IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= building_pars_f%fill ) &
4086	surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_tc3,j,i)
4087
4088	IF ( building_pars_f%pars_xy(12,j,i) /= building_pars_f%fill ) &
4089	surf_usm_h%target_temp_summer(m) = building_pars_f%pars_xy(12,j,i)
4090	IF ( building_pars_f%pars_xy(13,j,i) /= building_pars_f%fill ) &
4091	surf_usm_h%target_temp_winter(m) = building_pars_f%pars_xy(13,j,i)
4092
4093	IF ( building_pars_f%pars_xy(ind_emis_wall,j,i) /= building_pars_f%fill ) &
4094	surf_usm_h%emissivity(0,m) = building_pars_f%pars_xy(ind_emis_wall,j,i)
4095	IF ( building_pars_f%pars_xy(ind_emis_green,j,i) /= building_pars_f%fill )&
4096	surf_usm_h%emissivity(1,m) = building_pars_f%pars_xy(ind_emis_green,j,i)
4097	IF ( building_pars_f%pars_xy(ind_emis_win,j,i) /= building_pars_f%fill ) &
4098	surf_usm_h%emissivity(2,m) = building_pars_f%pars_xy(ind_emis_win,j,i)
4099
4100	IF ( building_pars_f%pars_xy(ind_trans,j,i) /= building_pars_f%fill ) &
4101	surf_usm_h%transmissivity(m) = building_pars_f%pars_xy(ind_trans,j,i)
4102
4103	IF ( building_pars_f%pars_xy(ind_z0,j,i) /= building_pars_f%fill ) &
4104	surf_usm_h%z0(m) = building_pars_f%pars_xy(ind_z0,j,i)
4105	IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill ) &
4106	surf_usm_h%z0h(m) = building_pars_f%pars_xy(ind_z0qh,j,i)
4107	IF ( building_pars_f%pars_xy(ind_z0qh,j,i) /= building_pars_f%fill ) &
4108	surf_usm_h%z0q(m) = building_pars_f%pars_xy(ind_z0qh,j,i)
4109
4110	IF ( building_pars_f%pars_xy(ind_alb_wall,j,i) /= building_pars_f%fill ) &
4111	surf_usm_h%albedo_type(0,m) = building_pars_f%pars_xy(ind_alb_wall,j,i)
4112	IF ( building_pars_f%pars_xy(ind_alb_green,j,i) /= building_pars_f%fill ) &
4113	surf_usm_h%albedo_type(1,m) = building_pars_f%pars_xy(ind_alb_green,j,i)
4114	IF ( building_pars_f%pars_xy(ind_alb_win,j,i) /= building_pars_f%fill ) &
4115	surf_usm_h%albedo_type(2,m) = building_pars_f%pars_xy(ind_alb_win,j,i)
4116
4117	IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) &
4118	surf_usm_h%zw(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i)
4119	IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) &
4120	surf_usm_h%zw(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i)
4121	IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) &
4122	surf_usm_h%zw(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i)
4123	IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) &
4124	surf_usm_h%zw(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i)
4125	IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) &
4126	surf_usm_h%zw_green(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i)
4127	IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) &
4128	surf_usm_h%zw_green(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i)
4129	IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) &
4130	surf_usm_h%zw_green(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i)
4131	IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) &
4132	surf_usm_h%zw_green(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i)
4133	IF ( building_pars_f%pars_xy(ind_thick_1,j,i) /= building_pars_f%fill ) &
4134	surf_usm_h%zw_window(nzb_wall,m) = building_pars_f%pars_xy(ind_thick_1,j,i)
4135	IF ( building_pars_f%pars_xy(ind_thick_2,j,i) /= building_pars_f%fill ) &
4136	surf_usm_h%zw_window(nzb_wall+1,m) = building_pars_f%pars_xy(ind_thick_2,j,i)
4137	IF ( building_pars_f%pars_xy(ind_thick_3,j,i) /= building_pars_f%fill ) &
4138	surf_usm_h%zw_window(nzb_wall+2,m) = building_pars_f%pars_xy(ind_thick_3,j,i)
4139	IF ( building_pars_f%pars_xy(ind_thick_4,j,i) /= building_pars_f%fill ) &
4140	surf_usm_h%zw_window(nzb_wall+3,m) = building_pars_f%pars_xy(ind_thick_4,j,i)
4141
4142	IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) &
4143	surf_usm_h%c_surface(m) = building_pars_f%pars_xy(45,j,i)
4144	IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) &
4145	surf_usm_h%lambda_surf(m) = building_pars_f%pars_xy(46,j,i)
4146	IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) &
4147	surf_usm_h%c_surface_green(m) = building_pars_f%pars_xy(45,j,i)
4148	IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) &
4149	surf_usm_h%lambda_surf_green(m) = building_pars_f%pars_xy(46,j,i)
4150	IF ( building_pars_f%pars_xy(45,j,i) /= building_pars_f%fill ) &
4151	surf_usm_h%c_surface_window(m) = building_pars_f%pars_xy(45,j,i)
4152	IF ( building_pars_f%pars_xy(46,j,i) /= building_pars_f%fill ) &
4153	surf_usm_h%lambda_surf_window(m) = building_pars_f%pars_xy(46,j,i)
4154	ENDDO
4155
4156
4157
4158	DO l = 0, 3
4159	DO m = 1, surf_usm_v(l)%ns
4160	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
4161	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
4162
4163	!
4164	!-- In order to distinguish between ground floor level and
4165	!-- above-ground-floor level surfaces, set input indices.
4166	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
4167	surf_usm_v(l)%ground_level(m) )
4168	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
4169	surf_usm_v(l)%ground_level(m) )
4170	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
4171	surf_usm_v(l)%ground_level(m) )
4172	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
4173	surf_usm_v(l)%ground_level(m) )
4174	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
4175	surf_usm_v(l)%ground_level(m) )
4176	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
4177	surf_usm_v(l)%ground_level(m) )
4178	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
4179	surf_usm_v(l)%ground_level(m) )
4180	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
4181	surf_usm_v(l)%ground_level(m) )
4182	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
4183	surf_usm_v(l)%ground_level(m) )
4184	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
4185	surf_usm_v(l)%ground_level(m) )
4186	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
4187	surf_usm_v(l)%ground_level(m) )
4188	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
4189	surf_usm_v(l)%ground_level(m) )
4190	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
4191	surf_usm_v(l)%ground_level(m) )
4192	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
4193	surf_usm_v(l)%ground_level(m) )
4194	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
4195	surf_usm_v(l)%ground_level(m) )
4196	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
4197	surf_usm_v(l)%ground_level(m) )
4198	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
4199	surf_usm_v(l)%ground_level(m) )
4200	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
4201	surf_usm_v(l)%ground_level(m) )
4202
4203	!
4204	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
4205	IF ( building_pars_f%pars_xy(ind_wall_frac,j,i) /= &
4206	building_pars_f%fill ) &
4207	surf_usm_v(l)%frac(0,m) = building_pars_f%pars_xy(ind_wall_frac,j,i)
4208	IF ( building_pars_f%pars_xy(ind_green_frac_w,j,i) /= &
4209	building_pars_f%fill ) &
4210	surf_usm_v(l)%frac(1,m) = building_pars_f%pars_xy(ind_green_frac_w,j,i)
4211	IF ( building_pars_f%pars_xy(ind_win_frac,j,i) /= &
4212	building_pars_f%fill ) &
4213	surf_usm_v(l)%frac(2,m) = building_pars_f%pars_xy(ind_win_frac,j,i)
4214
4215	IF ( building_pars_f%pars_xy(ind_lai_w,j,i) /= building_pars_f%fill ) &
4216	surf_usm_v(l)%lai(m) = building_pars_f%pars_xy(ind_lai_w,j,i)
4217
4218