Home

Context Navigation

urban_surface_mod.f90 @ 2718

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