Home

Context Navigation

urban_surface_mod.f90 @ 2723

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