Home

Context Navigation

urban_surface_mod.f90 @ 4493

Last change on this file since 4493 was 4493, checked in by pavelkrc, 5 years ago
Merge brach resler
Property svn:keywords set to `Id`
File size: 560.6 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-2020 Czech Technical University in Prague
18	! Copyright 2015-2020 Institute of Computer Science of the
19	! Czech Academy of Sciences, Prague
20	! Copyright 1997-2020 Leibniz Universitaet Hannover
21	!------------------------------------------------------------------------------!
22	!
23	! Current revisions:
24	! ------------------
25	!
26	!
27	! Former revisions:
28	! -----------------
29	! $Id: urban_surface_mod.f90 4493 2020-04-10 09:49:43Z pavelkrc $
30	! J.Resler, 2020/03/19
31	! - remove reading of deprecated input parameters c_surface and lambda_surf
32	! - and calculate them from parameters of the outer wall/roof layer
33	!
34	! 4481 2020-03-31 18:55:54Z maronga
35	! use statement for exchange horiz added
36	!
37	! 4442 2020-03-04 19:21:13Z suehring
38	! Change order of dimension in surface arrays %frac, %emissivity and %albedo
39	! to allow for better vectorization in the radiation interactions.
40	!
41	! 4441 2020-03-04 19:20:35Z suehring
42	! Removed wall_flags_static_0 from USE statements as it's not used within
43	! the module
44	!
45	! 4329 2019-12-10 15:46:36Z motisi
46	! Renamed wall_flags_0 to wall_flags_static_0
47	!
48	! 4309 2019-11-26 18:49:59Z suehring
49	! - Bugfix, include m_liq into restarts
50	! - Remove unused arrays for liquid water and saturation moisture at vertical
51	! walls
52	!
53	! 4305 2019-11-25 11:15:40Z suehring
54	! Revision of some indoor-model parameters
55	!
56	! 4259 2019-10-09 10:05:22Z suehring
57	! Instead of terminate the job in case the relative wall fractions do not
58	! sum-up to one, give only an informative message and normalize the fractions.
59	!
60	! 4258 2019-10-07 13:29:08Z suehring
61	! - Add checks to ensure that relative fractions of walls, windowns and green
62	! surfaces sum-up to one.
63	! - Revise message calls dealing with local checks.
64	!
65	! 4245 2019-09-30 08:40:37Z pavelkrc
66	! Initialize explicit per-surface parameters from building_surface_pars
67	!
68	! 4238 2019-09-25 16:06:01Z suehring
69	! Indoor-model parameters for some building types adjusted in order to avoid
70	! unrealistically high indoor temperatures (S. Rissmann)
71	!
72	! 4230 2019-09-11 13:58:14Z suehring
73	! Bugfix, initialize canopy resistance. Even if no green fraction is set,
74	! r_canopy must be initialized for output purposes.
75	!
76	! 4227 2019-09-10 18:04:34Z gronemeier
77	! implement new palm_date_time_mod
78	!
79	! 4214 2019-09-02 15:57:02Z suehring
80	! Bugfix, missing initialization and clearing of soil-moisture tendency
81	! (J.Resler)
82	!
83	! 4182 2019-08-22 15:20:23Z scharf
84	! Corrected "Former revisions" section
85	!
86	! 4168 2019-08-16 13:50:17Z suehring
87	! Replace function get_topography_top_index by topo_top_ind
88	!
89	! 4148 2019-08-08 11:26:00Z suehring
90	! - Add anthropogenic heat output factors for heating and cooling to building
91	! data base
92	! - Move definition of building_pars to usm_init_arrays since it is already
93	! required in the indoor model
94	!
95	! 4127 2019-07-30 14:47:10Z suehring
96	! Do not add anthopogenic energy during wall/soil spin-up
97	! (merge from branch resler)
98	!
99	! 4077 2019-07-09 13:27:11Z gronemeier
100	! Set roughness length z0 and z0h/q at ground-floor level to same value as
101	! those above ground-floor level
102	!
103	! 4051 2019-06-24 13:58:30Z suehring
104	! Remove work-around for green surface fraction on buildings
105	! (do not set it zero)
106	!
107	! 4050 2019-06-24 13:57:27Z suehring
108	! In order to avoid confusion with global control parameter, rename the
109	! USM-internal flag spinup into during_spinup.
110	!
111	! 3987 2019-05-22 09:52:13Z kanani
112	! Introduce alternative switch for debug output during timestepping
113	!
114	! 3943 2019-05-02 09:50:41Z maronga
115	! Removed qsws_eb. Bugfix in calculation of qsws.
116	!
117	! 3933 2019-04-25 12:33:20Z kanani
118	! Remove allocation of pt_2m, this is done in surface_mod now (surfaces%pt_2m)
119	!
120	! 3921 2019-04-18 14:21:10Z suehring
121	! Undo accidentally commented initialization
122	!
123	! 3918 2019-04-18 13:33:11Z suehring
124	! Set green fraction to zero also at vertical surfaces
125	!
126	! 3914 2019-04-17 16:02:02Z suehring
127	! In order to obtain correct surface temperature during spinup set window
128	! fraction to zero (only during spinup) instead of just disabling
129	! time-integration of window-surface temperature.
130	!
131	! 3901 2019-04-16 16:17:02Z suehring
132	! Workaround - set green fraction to zero ( green-heat model crashes ).
133	!
134	! 3896 2019-04-15 10:10:17Z suehring
135	!
136	!
137	! 3896 2019-04-15 10:10:17Z suehring
138	! Bugfix, wrong index used for accessing building_pars from PIDS
139	!
140	! 3885 2019-04-11 11:29:34Z kanani
141	! Changes related to global restructuring of location messages and introduction
142	! of additional debug messages
143	!
144	! 3882 2019-04-10 11:08:06Z suehring
145	! Avoid different type kinds
146	! Move definition of building-surface properties from declaration block
147	! to an extra routine
148	!
149	! 3881 2019-04-10 09:31:22Z suehring
150	! Revise determination of local ground-floor level height.
151	! Make level 3 initalization conform with Palm-input-data standard
152	! Move output of albedo and emissivity to radiation module
153	!
154	! 3832 2019-03-28 13:16:58Z raasch
155	! instrumented with openmp directives
156	!
157	! 3824 2019-03-27 15:56:16Z pavelkrc
158	! Remove unused imports
159	!
160	!
161	! 3814 2019-03-26 08:40:31Z pavelkrc
162	! unused subroutine commented out
163	!
164	! 3769 2019-02-28 10:16:49Z moh.hefny
165	! removed unused variables
166	!
167	! 3767 2019-02-27 08:18:02Z raasch
168	! unused variables removed from rrd-subroutines parameter list
169	!
170	! 3748 2019-02-18 10:38:31Z suehring
171	! Revise conversion of waste-heat flux (do not divide by air density, will
172	! be done in diffusion_s)
173	!
174	! 3745 2019-02-15 18:57:56Z suehring
175	! - Remove internal flag indoor_model (is a global control parameter)
176	! - add waste heat from buildings to the kinmatic heat flux
177	! - consider waste heat in restart data
178	! - remove unused USE statements
179	!
180	! 3744 2019-02-15 18:38:58Z suehring
181	! fixed surface heat capacity in the building parameters
182	! convert the file back to unix format
183	!
184	! 3730 2019-02-11 11:26:47Z moh.hefny
185	! Formatting and clean-up (rvtils)
186	!
187	! 3710 2019-01-30 18:11:19Z suehring
188	! Check if building type is set within a valid range.
189	!
190	! 3705 2019-01-29 19:56:39Z suehring
191	! make nzb_wall public, required for virtual-measurements
192	!
193	! 3704 2019-01-29 19:51:41Z suehring
194	! Some interface calls moved to module_interface + cleanup
195	!
196	! 3655 2019-01-07 16:51:22Z knoop
197	! Implementation of the PALM module interface
198	!
199	! 2007 2016-08-24 15:47:17Z kanani
200	! Initial revision
201	!
202	!
203	! Description:
204	! ------------
205	! 2016/6/9 - Initial version of the USM (Urban Surface Model)
206	! authors: Jaroslav Resler, Pavel Krc
207	! (Czech Technical University in Prague and Institute of
208	! Computer Science of the Czech Academy of Sciences, Prague)
209	! with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
210	! partly inspired by PALM LSM (B. Maronga)
211	! parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
212	!> Module for Urban Surface Model (USM)
213	!> The module includes:
214	!> 1. radiation model with direct/diffuse radiation, shading, reflections
215	!> and integration with plant canopy
216	!> 2. wall and wall surface model
217	!> 3. surface layer energy balance
218	!> 4. anthropogenic heat (only from transportation so far)
219	!> 5. necessary auxiliary subroutines (reading inputs, writing outputs,
220	!> restart simulations, ...)
221	!> It also make use of standard radiation and integrates it into
222	!> urban surface model.
223	!>
224	!> Further work:
225	!> -------------
226	!> @todo Output of _av variables in case of restarts
227	!> @todo Revise flux conversion in energy-balance solver
228	!> @todo Check divisions in wtend (etc.) calculations for possible division
229	!> by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
230	!> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
231	!> @todo remove reading of old csv inputs
232	!------------------------------------------------------------------------------!
233	MODULE urban_surface_mod
234
235	USE arrays_3d, &
236	ONLY: hyp, zu, pt, p, u, v, w, tend, exner, hyrho, prr, q, ql, vpt
237
238	USE calc_mean_profile_mod, &
239	ONLY: calc_mean_profile
240
241	USE basic_constants_and_equations_mod, &
242	ONLY: c_p, g, kappa, pi, r_d, rho_l, l_v, sigma_sb
243
244	USE control_parameters, &
245	ONLY: coupling_start_time, topography, &
246	debug_output, debug_output_timestep, debug_string, &
247	dt_3d, humidity, indoor_model, &
248	intermediate_timestep_count, initializing_actions, &
249	intermediate_timestep_count_max, simulated_time, end_time, &
250	timestep_scheme, tsc, coupling_char, io_blocks, io_group, &
251	message_string, time_since_reference_point, surface_pressure, &
252	pt_surface, large_scale_forcing, lsf_surf, &
253	spinup_pt_mean, spinup_time, time_do3d, dt_do3d, &
254	average_count_3d, varnamelength, urban_surface, dz
255
256	USE bulk_cloud_model_mod, &
257	ONLY: bulk_cloud_model, precipitation
258
259	USE cpulog, &
260	ONLY: cpu_log, log_point, log_point_s
261
262	USE grid_variables, &
263	ONLY: dx, dy, ddx, ddy, ddx2, ddy2
264
265	USE indices, &
266	ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, &
267	nysg, nzb, nzt, nbgp, topo_top_ind
268
269	USE, INTRINSIC :: iso_c_binding
270
271	USE kinds
272
273	USE palm_date_time_mod, &
274	ONLY: get_date_time, seconds_per_hour
275
276	USE pegrid
277
278	USE radiation_model_mod, &
279	ONLY: albedo_type, radiation_interaction, &
280	radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, &
281	force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u, &
282	iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id, &
283	nz_urban_b, nz_urban_t, unscheduled_radiation_calls
284
285	USE statistics, &
286	ONLY: hom, statistic_regions
287
288	USE surface_mod, &
289	ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_usm_h, &
290	surf_usm_v, surface_restore_elements
291
292
293	IMPLICIT NONE
294
295	!
296	!-- USM model constants
297
298	REAL(wp), PARAMETER :: &
299	b_ch = 6.04_wp, & !< Clapp & Hornberger exponent
300	lambda_h_green_dry = 0.19_wp, & !< heat conductivity for dry soil
301	lambda_h_green_sm = 3.44_wp, & !< heat conductivity of the soil matrix
302	lambda_h_water = 0.57_wp, & !< heat conductivity of water
303	psi_sat = -0.388_wp, & !< soil matrix potential at saturation
304	rho_c_soil = 2.19E6_wp, & !< volumetric heat capacity of soil
305	rho_c_water = 4.20E6_wp !< volumetric heat capacity of water
306	! m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
307
308	!
309	!-- Soil parameters I alpha_vg, l_vg_green, n_vg, gamma_w_green_sat
310	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/ &
311	3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & !< soil 1
312	3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & !< soil 2
313	0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, & !< soil 3
314	3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & !< soil 4
315	2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & !< soil 5
316	1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, & !< soil 6
317	0.00_wp, 0.00_wp, 0.00_wp, 0.57E-6_wp & !< soil 7
318	/), (/ 4, 7 /) )
319
320	!
321	!-- Soil parameters II swc_sat, fc, wilt, swc_res
322	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
323	0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & !< soil 1
324	0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & !< soil 2
325	0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & !< soil 3
326	0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & !< soil 4
327	0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & !< soil 5
328	0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & !< soil 6
329	0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp & !< soil 7
330	/), (/ 4, 7 /) )
331	!
332	!-- value 9999999.9_wp -> generic available or user-defined value must be set
333	!-- otherwise -> no generic variable and user setting is optional
334	REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !< NAMELIST alpha_vg
335	field_capacity = 9999999.9_wp, & !< NAMELIST fc
336	hydraulic_conductivity = 9999999.9_wp, & !< NAMELIST gamma_w_green_sat
337	l_vangenuchten = 9999999.9_wp, & !< NAMELIST l_vg
338	n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg
339	residual_moisture = 9999999.9_wp, & !< NAMELIST m_res
340	saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat
341	wilting_point = 9999999.9_wp !< NAMELIST m_wilt
342
343	!
344	!-- configuration parameters (they can be setup in PALM config)
345	LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used
346	LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources
347	!< (e.g.transportation) are used
348	LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls
349	LOGICAL :: read_wall_temp_3d = .FALSE.
350	LOGICAL :: usm_wall_mod = .FALSE. !< reduces conductivity of the first 2 wall layers by factor 0.1
351
352
353	INTEGER(iwp) :: building_type = 1 !< default building type (preleminary setting)
354	INTEGER(iwp) :: land_category = 2 !< default category for land surface
355	INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone
356	INTEGER(iwp) :: pedestrian_category = 2 !< default category for wall surface in pedestrian zone
357	INTEGER(iwp) :: roof_category = 2 !< default category for root surface
358	REAL(wp) :: roughness_concrete = 0.001_wp !< roughness length of average concrete surface
359	!
360	!-- Indices of input attributes in building_pars for (above) ground floor level
361	INTEGER(iwp) :: ind_alb_wall_agfl = 38 !< index in input list for albedo_type of wall above ground floor level
362	INTEGER(iwp) :: ind_alb_wall_gfl = 66 !< index in input list for albedo_type of wall ground floor level
363	INTEGER(iwp) :: ind_alb_wall_r = 101 !< index in input list for albedo_type of wall roof
364	INTEGER(iwp) :: ind_alb_green_agfl = 39 !< index in input list for albedo_type of green above ground floor level
365	INTEGER(iwp) :: ind_alb_green_gfl = 78 !< index in input list for albedo_type of green ground floor level
366	INTEGER(iwp) :: ind_alb_green_r = 117 !< index in input list for albedo_type of green roof
367	INTEGER(iwp) :: ind_alb_win_agfl = 40 !< index in input list for albedo_type of window fraction above ground floor level
368	INTEGER(iwp) :: ind_alb_win_gfl = 77 !< index in input list for albedo_type of window fraction ground floor level
369	INTEGER(iwp) :: ind_alb_win_r = 115 !< index in input list for albedo_type of window fraction roof
370	INTEGER(iwp) :: ind_emis_wall_agfl = 14 !< index in input list for wall emissivity, above ground floor level
371	INTEGER(iwp) :: ind_emis_wall_gfl = 32 !< index in input list for wall emissivity, ground floor level
372	INTEGER(iwp) :: ind_emis_wall_r = 100 !< index in input list for wall emissivity, roof
373	INTEGER(iwp) :: ind_emis_green_agfl = 15 !< index in input list for green emissivity, above ground floor level
374	INTEGER(iwp) :: ind_emis_green_gfl = 34 !< index in input list for green emissivity, ground floor level
375	INTEGER(iwp) :: ind_emis_green_r = 116 !< index in input list for green emissivity, roof
376	INTEGER(iwp) :: ind_emis_win_agfl = 16 !< index in input list for window emissivity, above ground floor level
377	INTEGER(iwp) :: ind_emis_win_gfl = 33 !< index in input list for window emissivity, ground floor level
378	INTEGER(iwp) :: ind_emis_win_r = 113 !< index in input list for window emissivity, roof
379	INTEGER(iwp) :: ind_gflh = 20 !< index in input list for ground floor level height
380	INTEGER(iwp) :: ind_green_frac_w_agfl = 2 !< index in input list for green fraction on wall, above ground floor level
381	INTEGER(iwp) :: ind_green_frac_w_gfl = 23 !< index in input list for green fraction on wall, ground floor level
382	INTEGER(iwp) :: ind_green_frac_r_agfl = 3 !< index in input list for green fraction on roof, above ground floor level
383	INTEGER(iwp) :: ind_green_frac_r_gfl = 24 !< index in input list for green fraction on roof, ground floor level
384	INTEGER(iwp) :: ind_hc1_agfl = 6 !< index in input list for heat capacity at first wall layer,
385	!< above ground floor level
386	INTEGER(iwp) :: ind_hc1_gfl = 26 !< index in input list for heat capacity at first wall layer, ground floor level
387	INTEGER(iwp) :: ind_hc1_wall_r = 94 !< index in input list for heat capacity at first wall layer, roof
388	INTEGER(iwp) :: ind_hc1_win_agfl = 83 !< index in input list for heat capacity at first window layer,
389	!< above ground floor level
390	INTEGER(iwp) :: ind_hc1_win_gfl = 71 !< index in input list for heat capacity at first window layer,
391	!< ground floor level
392	INTEGER(iwp) :: ind_hc1_win_r = 107 !< index in input list for heat capacity at first window layer, roof
393	INTEGER(iwp) :: ind_hc2_agfl = 7 !< index in input list for heat capacity at second wall layer,
394	!< above ground floor level
395	INTEGER(iwp) :: ind_hc2_gfl = 27 !< index in input list for heat capacity at second wall layer, ground floor level
396	INTEGER(iwp) :: ind_hc2_wall_r = 95 !< index in input list for heat capacity at second wall layer, roof
397	INTEGER(iwp) :: ind_hc2_win_agfl = 84 !< index in input list for heat capacity at second window layer,
398	!< above ground floor level
399	INTEGER(iwp) :: ind_hc2_win_gfl = 72 !< index in input list for heat capacity at second window layer,
400	!< ground floor level
401	INTEGER(iwp) :: ind_hc2_win_r = 108 !< index in input list for heat capacity at second window layer, roof
402	INTEGER(iwp) :: ind_hc3_agfl = 8 !< index in input list for heat capacity at third wall layer,
403	!< above ground floor level
404	INTEGER(iwp) :: ind_hc3_gfl = 28 !< index in input list for heat capacity at third wall layer, ground floor level
405	INTEGER(iwp) :: ind_hc3_wall_r = 96 !< index in input list for heat capacity at third wall layer, roof
406	INTEGER(iwp) :: ind_hc3_win_agfl = 85 !< index in input list for heat capacity at third window layer,
407	!< above ground floor level
408	INTEGER(iwp) :: ind_hc3_win_gfl = 73 !< index in input list for heat capacity at third window layer,
409	!< ground floor level
410	INTEGER(iwp) :: ind_hc3_win_r = 109 !< index in input list for heat capacity at third window layer, roof
411	INTEGER(iwp) :: ind_indoor_target_temp_summer = 12
412	INTEGER(iwp) :: ind_indoor_target_temp_winter = 13
413	INTEGER(iwp) :: ind_lai_r_agfl = 4 !< index in input list for LAI on roof, above ground floor level
414	INTEGER(iwp) :: ind_lai_r_gfl = 4 !< index in input list for LAI on roof, ground floor level
415	INTEGER(iwp) :: ind_lai_w_agfl = 5 !< index in input list for LAI on wall, above ground floor level
416	INTEGER(iwp) :: ind_lai_w_gfl = 25 !< index in input list for LAI on wall, ground floor level
417	INTEGER(iwp) :: ind_tc1_agfl = 9 !< index in input list for thermal conductivity at first wall layer,
418	!< above ground floor level
419	INTEGER(iwp) :: ind_tc1_gfl = 29 !< index in input list for thermal conductivity at first wall layer,
420	!< ground floor level
421	INTEGER(iwp) :: ind_tc1_wall_r = 97 !< index in input list for thermal conductivity at first wall layer, roof
422	INTEGER(iwp) :: ind_tc1_win_agfl = 86 !< index in input list for thermal conductivity at first window layer,
423	!< above ground floor level
424	INTEGER(iwp) :: ind_tc1_win_gfl = 74 !< index in input list for thermal conductivity at first window layer,
425	!< ground floor level
426	INTEGER(iwp) :: ind_tc1_win_r = 110 !< index in input list for thermal conductivity at first window layer, roof
427	INTEGER(iwp) :: ind_tc2_agfl = 10 !< index in input list for thermal conductivity at second wall layer,
428	!< above ground floor level
429	INTEGER(iwp) :: ind_tc2_gfl = 30 !< index in input list for thermal conductivity at second wall layer,
430	!< ground floor level
431	INTEGER(iwp) :: ind_tc2_wall_r = 98 !< index in input list for thermal conductivity at second wall layer, roof
432	INTEGER(iwp) :: ind_tc2_win_agfl = 87 !< index in input list for thermal conductivity at second window layer,
433	!< above ground floor level
434	INTEGER(iwp) :: ind_tc2_win_gfl = 75 !< index in input list for thermal conductivity at second window layer,
435	!< ground floor level
436	INTEGER(iwp) :: ind_tc2_win_r = 111 !< index in input list for thermal conductivity at second window layer,
437	!< ground floor level
438	INTEGER(iwp) :: ind_tc3_agfl = 11 !< index in input list for thermal conductivity at third wall layer,
439	!< above ground floor level
440	INTEGER(iwp) :: ind_tc3_gfl = 31 !< index in input list for thermal conductivity at third wall layer,
441	!< ground floor level
442	INTEGER(iwp) :: ind_tc3_wall_r = 99 !< index in input list for thermal conductivity at third wall layer, roof
443	INTEGER(iwp) :: ind_tc3_win_agfl = 88 !< index in input list for thermal conductivity at third window layer,
444	!< above ground floor level
445	INTEGER(iwp) :: ind_tc3_win_gfl = 76 !< index in input list for thermal conductivity at third window layer,
446	!< ground floor level
447	INTEGER(iwp) :: ind_tc3_win_r = 112 !< index in input list for thermal conductivity at third window layer, roof
448	INTEGER(iwp) :: ind_thick_1_agfl = 41 !< index for wall layer thickness - 1st layer above ground floor level
449	INTEGER(iwp) :: ind_thick_1_gfl = 62 !< index for wall layer thickness - 1st layer ground floor level
450	INTEGER(iwp) :: ind_thick_1_wall_r = 90 !< index for wall layer thickness - 1st layer roof
451	INTEGER(iwp) :: ind_thick_1_win_agfl = 79 !< index for window layer thickness - 1st layer above ground floor level
452	INTEGER(iwp) :: ind_thick_1_win_gfl = 67 !< index for window layer thickness - 1st layer ground floor level
453	INTEGER(iwp) :: ind_thick_1_win_r = 103 !< index for window layer thickness - 1st layer roof
454	INTEGER(iwp) :: ind_thick_2_agfl = 42 !< index for wall layer thickness - 2nd layer above ground floor level
455	INTEGER(iwp) :: ind_thick_2_gfl = 63 !< index for wall layer thickness - 2nd layer ground floor level
456	INTEGER(iwp) :: ind_thick_2_wall_r = 91 !< index for wall layer thickness - 2nd layer roof
457	INTEGER(iwp) :: ind_thick_2_win_agfl = 80 !< index for window layer thickness - 2nd layer above ground floor level
458	INTEGER(iwp) :: ind_thick_2_win_gfl = 68 !< index for window layer thickness - 2nd layer ground floor level
459	INTEGER(iwp) :: ind_thick_2_win_r = 104 !< index for window layer thickness - 2nd layer roof
460	INTEGER(iwp) :: ind_thick_3_agfl = 43 !< index for wall layer thickness - 3rd layer above ground floor level
461	INTEGER(iwp) :: ind_thick_3_gfl = 64 !< index for wall layer thickness - 3rd layer ground floor level
462	INTEGER(iwp) :: ind_thick_3_wall_r = 92 !< index for wall layer thickness - 3rd layer roof
463	INTEGER(iwp) :: ind_thick_3_win_agfl = 81 !< index for window layer thickness - 3rd layer above ground floor level
464	INTEGER(iwp) :: ind_thick_3_win_gfl = 69 !< index for window layer thickness - 3rd layer ground floor level
465	INTEGER(iwp) :: ind_thick_3_win_r = 105 !< index for window layer thickness - 3rd layer roof
466	INTEGER(iwp) :: ind_thick_4_agfl = 44 !< index for wall layer thickness - 4th layer above ground floor level
467	INTEGER(iwp) :: ind_thick_4_gfl = 65 !< index for wall layer thickness - 4th layer ground floor level
468	INTEGER(iwp) :: ind_thick_4_wall_r = 93 !< index for wall layer thickness - 4st layer roof
469	INTEGER(iwp) :: ind_thick_4_win_agfl = 82 !< index for window layer thickness - 4th layer above ground floor level
470	INTEGER(iwp) :: ind_thick_4_win_gfl = 70 !< index for window layer thickness - 4th layer ground floor level
471	INTEGER(iwp) :: ind_thick_4_win_r = 106 !< index for window layer thickness - 4th layer roof
472	INTEGER(iwp) :: ind_trans_agfl = 17 !< index in input list for window transmissivity, above ground floor level
473	INTEGER(iwp) :: ind_trans_gfl = 35 !< index in input list for window transmissivity, ground floor level
474	INTEGER(iwp) :: ind_trans_r = 114 !< index in input list for window transmissivity, roof
475	INTEGER(iwp) :: ind_wall_frac_agfl = 0 !< index in input list for wall fraction, above ground floor level
476	INTEGER(iwp) :: ind_wall_frac_gfl = 21 !< index in input list for wall fraction, ground floor level
477	INTEGER(iwp) :: ind_wall_frac_r = 89 !< index in input list for wall fraction, roof
478	INTEGER(iwp) :: ind_win_frac_agfl = 1 !< index in input list for window fraction, above ground floor level
479	INTEGER(iwp) :: ind_win_frac_gfl = 22 !< index in input list for window fraction, ground floor level
480	INTEGER(iwp) :: ind_win_frac_r = 102 !< index in input list for window fraction, roof
481	INTEGER(iwp) :: ind_z0_agfl = 18 !< index in input list for z0, above ground floor level
482	INTEGER(iwp) :: ind_z0_gfl = 36 !< index in input list for z0, ground floor level
483	INTEGER(iwp) :: ind_z0qh_agfl = 19 !< index in input list for z0h / z0q, above ground floor level
484	INTEGER(iwp) :: ind_z0qh_gfl = 37 !< index in input list for z0h / z0q, ground floor level
485	INTEGER(iwp) :: ind_green_type_roof = 118 !< index in input list for type of green roof
486	!
487	!-- Indices of input attributes in building_surface_pars (except for
488	!-- radiation-related, which are in radiation_model_mod)
489	INTEGER(iwp) :: ind_s_wall_frac = 0 !< index for wall fraction (0-1)
490	INTEGER(iwp) :: ind_s_win_frac = 1 !< index for window fraction (0-1)
491	INTEGER(iwp) :: ind_s_green_frac_w = 2 !< index for green fraction on wall (0-1)
492	INTEGER(iwp) :: ind_s_green_frac_r = 3 !< index for green fraction on roof (0-1)
493	INTEGER(iwp) :: ind_s_lai_r = 4 !< index for leaf area index of green fraction
494	INTEGER(iwp) :: ind_s_hc1 = 5 !< index for heat capacity of wall layer 1
495	INTEGER(iwp) :: ind_s_hc2 = 6 !< index for heat capacity of wall layer 2
496	INTEGER(iwp) :: ind_s_hc3 = 7 !< index for heat capacity of wall layer 3
497	INTEGER(iwp) :: ind_s_tc1 = 8 !< index for thermal conducivity of wall layer 1
498	INTEGER(iwp) :: ind_s_tc2 = 9 !< index for thermal conducivity of wall layer 2
499	INTEGER(iwp) :: ind_s_tc3 = 10 !< index for thermal conducivity of wall layer 3
500	INTEGER(iwp) :: ind_s_indoor_target_temp_summer = 11 !< index for indoor target summer temperature
501	INTEGER(iwp) :: ind_s_indoor_target_temp_winter = 12 !< index for indoor target winter temperature
502	INTEGER(iwp) :: ind_s_emis_wall = 13 !< index for emissivity of wall fraction (0-1)
503	INTEGER(iwp) :: ind_s_emis_green = 14 !< index for emissivity of green fraction (0-1)
504	INTEGER(iwp) :: ind_s_emis_win = 15 !< index for emissivity o f window fraction (0-1)
505	INTEGER(iwp) :: ind_s_trans = 16 !< index for transmissivity of window fraction (0-1)
506	INTEGER(iwp) :: ind_s_z0 = 17 !< index for roughness length for momentum (m)
507	INTEGER(iwp) :: ind_s_z0qh = 18 !< index for roughness length for heat (m)
508
509	REAL(wp) :: roof_height_limit = 4.0_wp !< height for distinguish between land surfaces and roofs
510	REAL(wp) :: ground_floor_level = 4.0_wp !< default ground floor level
511
512
513	CHARACTER(37), DIMENSION(0:7), PARAMETER :: building_type_name = (/ &
514	'user-defined ', & !< type 0
515	'residential - 1950 ', & !< type 1
516	'residential 1951 - 2000 ', & !< type 2
517	'residential 2001 - ', & !< type 3
518	'office - 1950 ', & !< type 4
519	'office 1951 - 2000 ', & !< type 5
520	'office 2001 - ', & !< type 6
521	'bridges ' & !< type 7
522	/)
523
524
525	!
526	!-- Building facade/wall/green/window properties (partly according to PIDS).
527	!-- Initialization of building_pars is outsourced to usm_init_pars. This is
528	!-- needed because of the huge number of attributes given in building_pars
529	!-- (>700), while intel and gfortran compiler have hard limit of continuation
530	!-- lines of 511.
531	REAL(wp), DIMENSION(0:135,1:7) :: building_pars
532	!
533	!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
534	TYPE t_surf_vertical
535	REAL(wp), DIMENSION(:), ALLOCATABLE :: t
536	END TYPE t_surf_vertical
537	!
538	!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
539	TYPE t_wall_vertical
540	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t
541	END TYPE t_wall_vertical
542
543	TYPE surf_type_usm
544	REAL(wp), DIMENSION(:), ALLOCATABLE :: var_usm_1d !< 1D prognostic variable
545	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_usm_2d !< 2D prognostic variable
546	END TYPE surf_type_usm
547
548	TYPE(surf_type_usm), POINTER :: m_liq_usm_h, & !< liquid water reservoir (m), horizontal surface elements
549	m_liq_usm_h_p !< progn. liquid water reservoir (m), horizontal surface elements
550
551	TYPE(surf_type_usm), TARGET :: m_liq_usm_h_1, & !<
552	m_liq_usm_h_2 !<
553
554	TYPE(surf_type_usm), TARGET :: tm_liq_usm_h_m !< liquid water reservoir tendency (m), horizontal surface elements
555	!
556	!-- anthropogenic heat sources
557	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2)
558	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheatprof !< diurnal profiles of anthropogenic heat
559	!< for particular layers
560	INTEGER(iwp) :: naheatlayers = 1 !< number of layers of anthropogenic heat
561
562	!
563	!-- wall surface model
564	!-- wall surface model constants
565	INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched)
566	INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched)
567	INTEGER(iwp), PARAMETER :: nzw = 4 !< number of wall layers (fixed for now)
568
569	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
570	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
571	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
572	!< normalized soil, wall and roof, window and
573	!<green layer depths (m/m)
574
575	REAL(wp) :: wall_inner_temperature = 295.0_wp !< temperature of the inner wall
576	!< surface (~22 degrees C) (K)
577	REAL(wp) :: roof_inner_temperature = 295.0_wp !< temperature of the inner roof
578	!< surface (~22 degrees C) (K)
579	REAL(wp) :: soil_inner_temperature = 288.0_wp !< temperature of the deep soil
580	!< (~15 degrees C) (K)
581	REAL(wp) :: window_inner_temperature = 295.0_wp !< temperature of the inner window
582	!< surface (~22 degrees C) (K)
583
584	REAL(wp) :: m_total = 0.0_wp !< weighted total water content of the soil (m3/m3)
585	INTEGER(iwp) :: soil_type
586
587	!
588	!-- surface and material model variables for walls, ground, roofs
589	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m)
590	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_window !< normalized window layer depths (m)
591	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_green !< normalized green layer depths (m)
592
593	REAL(wp), DIMENSION(:), POINTER :: t_surf_wall_h
594	REAL(wp), DIMENSION(:), POINTER :: t_surf_wall_h_p
595	REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h
596	REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h_p
597	REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h
598	REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h_p
599
600	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_wall_h_1
601	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_wall_h_2
602	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_1
603	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_2
604	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_1
605	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_2
606
607	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_wall_v
608	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_wall_v_p
609	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v
610	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v_p
611	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v
612	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v_p
613
614	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_wall_v_1
615	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_wall_v_2
616	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_1
617	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_2
618	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_1
619	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_2
620
621	!
622	!-- Energy balance variables
623	!-- parameters of the land, roof and wall surfaces
624
625	REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p
626	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_1, t_wall_h_2
627	REAL(wp), DIMENSION(:,:), POINTER :: t_window_h, t_window_h_p
628	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_1, t_window_h_2
629	REAL(wp), DIMENSION(:,:), POINTER :: t_green_h, t_green_h_p
630	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_1, t_green_h_2
631	REAL(wp), DIMENSION(:,:), POINTER :: swc_h, rootfr_h, wilt_h, fc_h, swc_sat_h, swc_h_p, swc_res_h
632	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: swc_h_1, rootfr_h_1, &
633	wilt_h_1, fc_h_1, swc_sat_h_1, swc_h_2, swc_res_h_1
634
635
636	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p
637	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_1, t_wall_v_2
638	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_window_v, t_window_v_p
639	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_1, t_window_v_2
640	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_green_v, t_green_v_p
641	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_1, t_green_v_2
642	!
643	!-- Surface and material parameters classes (surface_type)
644	!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
645	INTEGER(iwp) :: n_surface_types !< number of the wall type categories
646	INTEGER(iwp), PARAMETER :: n_surface_params = 9 !< number of parameters for each type of the wall
647	INTEGER(iwp), PARAMETER :: ialbedo = 1 !< albedo of the surface
648	INTEGER(iwp), PARAMETER :: iemiss = 2 !< emissivity of the surface
649	INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity lambda S between surface
650	!< and material ( W m-2 K-1 )
651	INTEGER(iwp), PARAMETER :: irough = 4 !< roughness length z0 for movements
652	INTEGER(iwp), PARAMETER :: iroughh = 5 !< roughness length z0h for scalars
653	!< (heat, humidity,...)
654	INTEGER(iwp), PARAMETER :: icsurf = 6 !< Surface skin layer heat capacity (J m-2 K-1 )
655	INTEGER(iwp), PARAMETER :: ithick = 7 !< thickness of the surface (wall, roof, land) ( m )
656	INTEGER(iwp), PARAMETER :: irhoC = 8 !< volumetric heat capacity rho*C of
657	!< the material ( J m-3 K-1 )
658	INTEGER(iwp), PARAMETER :: ilambdah = 9 !< thermal conductivity lambda H
659	!< of the wall (W m-1 K-1 )
660	CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports)
661	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types
662	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: surface_params !< parameters of wall types
663
664	!
665	!-- interfaces of subroutines accessed from outside of this module
666	INTERFACE usm_3d_data_averaging
667	MODULE PROCEDURE usm_3d_data_averaging
668	END INTERFACE usm_3d_data_averaging
669
670	INTERFACE usm_boundary_condition
671	MODULE PROCEDURE usm_boundary_condition
672	END INTERFACE usm_boundary_condition
673
674	INTERFACE usm_check_data_output
675	MODULE PROCEDURE usm_check_data_output
676	END INTERFACE usm_check_data_output
677
678	INTERFACE usm_check_parameters
679	MODULE PROCEDURE usm_check_parameters
680	END INTERFACE usm_check_parameters
681
682	INTERFACE usm_data_output_3d
683	MODULE PROCEDURE usm_data_output_3d
684	END INTERFACE usm_data_output_3d
685
686	INTERFACE usm_define_netcdf_grid
687	MODULE PROCEDURE usm_define_netcdf_grid
688	END INTERFACE usm_define_netcdf_grid
689
690	INTERFACE usm_init
691	MODULE PROCEDURE usm_init
692	END INTERFACE usm_init
693
694	INTERFACE usm_init_arrays
695	MODULE PROCEDURE usm_init_arrays
696	END INTERFACE usm_init_arrays
697
698	INTERFACE usm_material_heat_model
699	MODULE PROCEDURE usm_material_heat_model
700	END INTERFACE usm_material_heat_model
701
702	INTERFACE usm_green_heat_model
703	MODULE PROCEDURE usm_green_heat_model
704	END INTERFACE usm_green_heat_model
705
706	INTERFACE usm_parin
707	MODULE PROCEDURE usm_parin
708	END INTERFACE usm_parin
709
710	INTERFACE usm_rrd_local
711	MODULE PROCEDURE usm_rrd_local
712	END INTERFACE usm_rrd_local
713
714	INTERFACE usm_surface_energy_balance
715	MODULE PROCEDURE usm_surface_energy_balance
716	END INTERFACE usm_surface_energy_balance
717
718	INTERFACE usm_swap_timelevel
719	MODULE PROCEDURE usm_swap_timelevel
720	END INTERFACE usm_swap_timelevel
721
722	INTERFACE usm_wrd_local
723	MODULE PROCEDURE usm_wrd_local
724	END INTERFACE usm_wrd_local
725
726
727	SAVE
728
729	PRIVATE
730
731	!
732	!-- Public functions
733	PUBLIC usm_boundary_condition, usm_check_parameters, usm_init, &
734	usm_rrd_local, &
735	usm_surface_energy_balance, usm_material_heat_model, &
736	usm_swap_timelevel, usm_check_data_output, usm_3d_data_averaging, &
737	usm_data_output_3d, usm_define_netcdf_grid, usm_parin, &
738	usm_wrd_local, usm_init_arrays
739
740	!
741	!-- Public parameters, constants and initial values
742	PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod, &
743	usm_green_heat_model, building_pars, &
744	nzb_wall, nzt_wall, t_wall_h, t_wall_v, &
745	t_window_h, t_window_v, building_type
746
747
748
749	CONTAINS
750
751	!------------------------------------------------------------------------------!
752	! Description:
753	! ------------
754	!> This subroutine creates the necessary indices of the urban surfaces
755	!> and plant canopy and it allocates the needed arrays for USM
756	!------------------------------------------------------------------------------!
757	SUBROUTINE usm_init_arrays
758
759	IMPLICIT NONE
760
761	INTEGER(iwp) :: l
762
763	IF ( debug_output ) CALL debug_message( 'usm_init_arrays', 'start' )
764
765	!
766	!-- Allocate radiation arrays which are part of the new data type.
767	!-- For horizontal surfaces.
768	ALLOCATE ( surf_usm_h%surfhf(1:surf_usm_h%ns) )
769	ALLOCATE ( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
770	!
771	!-- For vertical surfaces
772	DO l = 0, 3
773	ALLOCATE ( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) )
774	ALLOCATE ( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
775	ENDDO
776
777	!
778	!-- Wall surface model
779	!-- allocate arrays for wall surface model and define pointers
780	!-- allocate array of wall types and wall parameters
781	ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
782	ALLOCATE ( surf_usm_h%building_type(1:surf_usm_h%ns) )
783	ALLOCATE ( surf_usm_h%building_type_name(1:surf_usm_h%ns) )
784	surf_usm_h%building_type = 0
785	surf_usm_h%building_type_name = 'none'
786	DO l = 0, 3
787	ALLOCATE ( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
788	ALLOCATE ( surf_usm_v(l)%building_type(1:surf_usm_v(l)%ns) )
789	ALLOCATE ( surf_usm_v(l)%building_type_name(1:surf_usm_v(l)%ns) )
790	surf_usm_v(l)%building_type = 0
791	surf_usm_v(l)%building_type_name = 'none'
792	ENDDO
793	!
794	!-- Allocate albedo_type and albedo. Each surface element
795	!-- has 3 values, 0: wall fraction, 1: green fraction, 2: window fraction.
796	ALLOCATE ( surf_usm_h%albedo_type(1:surf_usm_h%ns,0:2) )
797	ALLOCATE ( surf_usm_h%albedo(1:surf_usm_h%ns,0:2) )
798	surf_usm_h%albedo_type = albedo_type
799	DO l = 0, 3
800	ALLOCATE ( surf_usm_v(l)%albedo_type(1:surf_usm_v(l)%ns,0:2) )
801	ALLOCATE ( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2) )
802	surf_usm_v(l)%albedo_type = albedo_type
803	ENDDO
804
805	!
806	!-- Allocate indoor target temperature for summer and winter
807	ALLOCATE ( surf_usm_h%target_temp_summer(1:surf_usm_h%ns) )
808	ALLOCATE ( surf_usm_h%target_temp_winter(1:surf_usm_h%ns) )
809	DO l = 0, 3
810	ALLOCATE ( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) )
811	ALLOCATE ( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) )
812	ENDDO
813	!
814	!-- In case the indoor model is applied, allocate memory for waste heat
815	!-- and indoor temperature.
816	IF ( indoor_model ) THEN
817	ALLOCATE ( surf_usm_h%waste_heat(1:surf_usm_h%ns) )
818	surf_usm_h%waste_heat = 0.0_wp
819	DO l = 0, 3
820	ALLOCATE ( surf_usm_v(l)%waste_heat(1:surf_usm_v(l)%ns) )
821	surf_usm_v(l)%waste_heat = 0.0_wp
822	ENDDO
823	ENDIF
824	!
825	!-- Allocate flag indicating ground floor level surface elements
826	ALLOCATE ( surf_usm_h%ground_level(1:surf_usm_h%ns) )
827	DO l = 0, 3
828	ALLOCATE ( surf_usm_v(l)%ground_level(1:surf_usm_v(l)%ns) )
829	ENDDO
830	!
831	!-- Allocate arrays for relative surface fraction.
832	!-- 0 - wall fraction, 1 - green fraction, 2 - window fraction
833	ALLOCATE ( surf_usm_h%frac(1:surf_usm_h%ns,0:2) )
834	surf_usm_h%frac = 0.0_wp
835	DO l = 0, 3
836	ALLOCATE ( surf_usm_v(l)%frac(1:surf_usm_v(l)%ns,0:2) )
837	surf_usm_v(l)%frac = 0.0_wp
838	ENDDO
839
840	!
841	!-- wall and roof surface parameters. First for horizontal surfaces
842	ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) )
843	ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) )
844	ALLOCATE ( surf_usm_h%lambda_surf_window(1:surf_usm_h%ns) )
845	ALLOCATE ( surf_usm_h%lambda_surf_green(1:surf_usm_h%ns) )
846	ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) )
847	ALLOCATE ( surf_usm_h%c_surface_window(1:surf_usm_h%ns) )
848	ALLOCATE ( surf_usm_h%c_surface_green(1:surf_usm_h%ns) )
849	ALLOCATE ( surf_usm_h%transmissivity(1:surf_usm_h%ns) )
850	ALLOCATE ( surf_usm_h%lai(1:surf_usm_h%ns) )
851	ALLOCATE ( surf_usm_h%emissivity(1:surf_usm_h%ns,0:2) )
852	ALLOCATE ( surf_usm_h%r_a(1:surf_usm_h%ns) )
853	ALLOCATE ( surf_usm_h%r_a_green(1:surf_usm_h%ns) )
854	ALLOCATE ( surf_usm_h%r_a_window(1:surf_usm_h%ns) )
855	ALLOCATE ( surf_usm_h%green_type_roof(1:surf_usm_h%ns) )
856	ALLOCATE ( surf_usm_h%r_s(1:surf_usm_h%ns) )
857
858	!
859	!-- For vertical surfaces.
860	DO l = 0, 3
861	ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) )
862	ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) )
863	ALLOCATE ( surf_usm_v(l)%lambda_surf_window(1:surf_usm_v(l)%ns) )
864	ALLOCATE ( surf_usm_v(l)%c_surface_window(1:surf_usm_v(l)%ns) )
865	ALLOCATE ( surf_usm_v(l)%lambda_surf_green(1:surf_usm_v(l)%ns) )
866	ALLOCATE ( surf_usm_v(l)%c_surface_green(1:surf_usm_v(l)%ns) )
867	ALLOCATE ( surf_usm_v(l)%transmissivity(1:surf_usm_v(l)%ns) )
868	ALLOCATE ( surf_usm_v(l)%lai(1:surf_usm_v(l)%ns) )
869	ALLOCATE ( surf_usm_v(l)%emissivity(1:surf_usm_v(l)%ns,0:2) )
870	ALLOCATE ( surf_usm_v(l)%r_a(1:surf_usm_v(l)%ns) )
871	ALLOCATE ( surf_usm_v(l)%r_a_green(1:surf_usm_v(l)%ns) )
872	ALLOCATE ( surf_usm_v(l)%r_a_window(1:surf_usm_v(l)%ns) )
873	ALLOCATE ( surf_usm_v(l)%r_s(1:surf_usm_v(l)%ns) )
874	ENDDO
875
876	!
877	!-- allocate wall and roof material parameters. First for horizontal surfaces
878	ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) )
879	ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns) )
880	ALLOCATE ( surf_usm_h%thickness_green(1:surf_usm_h%ns) )
881	ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
882	ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
883	ALLOCATE ( surf_usm_h%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
884	ALLOCATE ( surf_usm_h%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
885	ALLOCATE ( surf_usm_h%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
886	ALLOCATE ( surf_usm_h%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
887
888	ALLOCATE ( surf_usm_h%rho_c_total_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
889	ALLOCATE ( surf_usm_h%n_vg_green(1:surf_usm_h%ns) )
890	ALLOCATE ( surf_usm_h%alpha_vg_green(1:surf_usm_h%ns) )
891	ALLOCATE ( surf_usm_h%l_vg_green(1:surf_usm_h%ns) )
892	ALLOCATE ( surf_usm_h%gamma_w_green_sat(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
893	ALLOCATE ( surf_usm_h%lambda_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
894	ALLOCATE ( surf_usm_h%gamma_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
895	ALLOCATE ( surf_usm_h%tswc_h_m(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
896
897	!
898	!-- For vertical surfaces.
899	DO l = 0, 3
900	ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) )
901	ALLOCATE ( surf_usm_v(l)%thickness_window(1:surf_usm_v(l)%ns) )
902	ALLOCATE ( surf_usm_v(l)%thickness_green(1:surf_usm_v(l)%ns) )
903	ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
904	ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
905	ALLOCATE ( surf_usm_v(l)%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
906	ALLOCATE ( surf_usm_v(l)%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
907	ALLOCATE ( surf_usm_v(l)%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
908	ALLOCATE ( surf_usm_v(l)%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
909	ENDDO
910
911	!
912	!-- allocate green wall and roof vegetation and soil parameters. First horizontal surfaces
913	ALLOCATE ( surf_usm_h%g_d(1:surf_usm_h%ns) )
914	ALLOCATE ( surf_usm_h%c_liq(1:surf_usm_h%ns) )
915	ALLOCATE ( surf_usm_h%qsws_liq(1:surf_usm_h%ns) )
916	ALLOCATE ( surf_usm_h%qsws_veg(1:surf_usm_h%ns) )
917	ALLOCATE ( surf_usm_h%r_canopy(1:surf_usm_h%ns) )
918	ALLOCATE ( surf_usm_h%r_canopy_min(1:surf_usm_h%ns) )
919	ALLOCATE ( surf_usm_h%pt_10cm(1:surf_usm_h%ns) )
920
921	!
922	!-- For vertical surfaces.
923	DO l = 0, 3
924	ALLOCATE ( surf_usm_v(l)%g_d(1:surf_usm_v(l)%ns) )
925	ALLOCATE ( surf_usm_v(l)%c_liq(1:surf_usm_v(l)%ns) )
926	ALLOCATE ( surf_usm_v(l)%qsws_liq(1:surf_usm_v(l)%ns) )
927	ALLOCATE ( surf_usm_v(l)%qsws_veg(1:surf_usm_v(l)%ns) )
928	ALLOCATE ( surf_usm_v(l)%r_canopy(1:surf_usm_v(l)%ns) )
929	ALLOCATE ( surf_usm_v(l)%r_canopy_min(1:surf_usm_v(l)%ns) )
930	ALLOCATE ( surf_usm_v(l)%pt_10cm(1:surf_usm_v(l)%ns) )
931	ENDDO
932
933	!
934	!-- allocate wall and roof layers sizes. For horizontal surfaces.
935	ALLOCATE ( zwn(nzb_wall:nzt_wall) )
936	ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
937	ALLOCATE ( zwn_window(nzb_wall:nzt_wall) )
938	ALLOCATE ( surf_usm_h%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
939	ALLOCATE ( zwn_green(nzb_wall:nzt_wall) )
940	ALLOCATE ( surf_usm_h%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
941	ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
942	ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
943	ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
944	ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
945	ALLOCATE ( surf_usm_h%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
946	ALLOCATE ( surf_usm_h%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
947	ALLOCATE ( surf_usm_h%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
948	ALLOCATE ( surf_usm_h%zw_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
949	ALLOCATE ( surf_usm_h%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
950	ALLOCATE ( surf_usm_h%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
951	ALLOCATE ( surf_usm_h%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
952	ALLOCATE ( surf_usm_h%zw_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
953
954	!
955	!-- For vertical surfaces.
956	DO l = 0, 3
957	ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
958	ALLOCATE ( surf_usm_v(l)%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
959	ALLOCATE ( surf_usm_v(l)%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
960	ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
961	ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
962	ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
963	ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
964	ALLOCATE ( surf_usm_v(l)%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
965	ALLOCATE ( surf_usm_v(l)%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
966	ALLOCATE ( surf_usm_v(l)%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
967	ALLOCATE ( surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
968	ALLOCATE ( surf_usm_v(l)%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
969	ALLOCATE ( surf_usm_v(l)%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
970	ALLOCATE ( surf_usm_v(l)%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
971	ALLOCATE ( surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
972	ENDDO
973
974	!
975	!-- allocate wall and roof temperature arrays, for horizontal walls
976	!
977	!-- Allocate if required. Note, in case of restarts, some of these arrays
978	!-- might be already allocated.
979	IF ( .NOT. ALLOCATED( t_surf_wall_h_1 ) ) &
980	ALLOCATE ( t_surf_wall_h_1(1:surf_usm_h%ns) )
981	IF ( .NOT. ALLOCATED( t_surf_wall_h_2 ) ) &
982	ALLOCATE ( t_surf_wall_h_2(1:surf_usm_h%ns) )
983	IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) &
984	ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
985	IF ( .NOT. ALLOCATED( t_wall_h_2 ) ) &
986	ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
987	IF ( .NOT. ALLOCATED( t_surf_window_h_1 ) ) &
988	ALLOCATE ( t_surf_window_h_1(1:surf_usm_h%ns) )
989	IF ( .NOT. ALLOCATED( t_surf_window_h_2 ) ) &
990	ALLOCATE ( t_surf_window_h_2(1:surf_usm_h%ns) )
991	IF ( .NOT. ALLOCATED( t_window_h_1 ) ) &
992	ALLOCATE ( t_window_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
993	IF ( .NOT. ALLOCATED( t_window_h_2 ) ) &
994	ALLOCATE ( t_window_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
995	IF ( .NOT. ALLOCATED( t_surf_green_h_1 ) ) &
996	ALLOCATE ( t_surf_green_h_1(1:surf_usm_h%ns) )
997	IF ( .NOT. ALLOCATED( t_surf_green_h_2 ) ) &
998	ALLOCATE ( t_surf_green_h_2(1:surf_usm_h%ns) )
999	IF ( .NOT. ALLOCATED( t_green_h_1 ) ) &
1000	ALLOCATE ( t_green_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1001	IF ( .NOT. ALLOCATED( t_green_h_2 ) ) &
1002	ALLOCATE ( t_green_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1003	IF ( .NOT. ALLOCATED( swc_h_1 ) ) &
1004	ALLOCATE ( swc_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1005	IF ( .NOT. ALLOCATED( swc_sat_h_1 ) ) &
1006	ALLOCATE ( swc_sat_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1007	IF ( .NOT. ALLOCATED( swc_res_h_1 ) ) &
1008	ALLOCATE ( swc_res_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1009	IF ( .NOT. ALLOCATED( swc_h_2 ) ) &
1010	ALLOCATE ( swc_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1011	IF ( .NOT. ALLOCATED( rootfr_h_1 ) ) &
1012	ALLOCATE ( rootfr_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1013	IF ( .NOT. ALLOCATED( wilt_h_1 ) ) &
1014	ALLOCATE ( wilt_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1015	IF ( .NOT. ALLOCATED( fc_h_1 ) ) &
1016	ALLOCATE ( fc_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1017
1018	IF ( .NOT. ALLOCATED( m_liq_usm_h_1%var_usm_1d ) ) &
1019	ALLOCATE ( m_liq_usm_h_1%var_usm_1d(1:surf_usm_h%ns) )
1020	IF ( .NOT. ALLOCATED( m_liq_usm_h_2%var_usm_1d ) ) &
1021	ALLOCATE ( m_liq_usm_h_2%var_usm_1d(1:surf_usm_h%ns) )
1022
1023	!
1024	!-- initial assignment of the pointers
1025	t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2
1026	t_window_h => t_window_h_1; t_window_h_p => t_window_h_2
1027	t_green_h => t_green_h_1; t_green_h_p => t_green_h_2
1028	t_surf_wall_h => t_surf_wall_h_1; t_surf_wall_h_p => t_surf_wall_h_2
1029	t_surf_window_h => t_surf_window_h_1; t_surf_window_h_p => t_surf_window_h_2
1030	t_surf_green_h => t_surf_green_h_1; t_surf_green_h_p => t_surf_green_h_2
1031	m_liq_usm_h => m_liq_usm_h_1; m_liq_usm_h_p => m_liq_usm_h_2
1032	swc_h => swc_h_1; swc_h_p => swc_h_2
1033	swc_sat_h => swc_sat_h_1
1034	swc_res_h => swc_res_h_1
1035	rootfr_h => rootfr_h_1
1036	wilt_h => wilt_h_1
1037	fc_h => fc_h_1
1038
1039	!
1040	!-- allocate wall and roof temperature arrays, for vertical walls if required
1041	!
1042	!-- Allocate if required. Note, in case of restarts, some of these arrays
1043	!-- might be already allocated.
1044	DO l = 0, 3
1045	IF ( .NOT. ALLOCATED( t_surf_wall_v_1(l)%t ) ) &
1046	ALLOCATE ( t_surf_wall_v_1(l)%t(1:surf_usm_v(l)%ns) )
1047	IF ( .NOT. ALLOCATED( t_surf_wall_v_2(l)%t ) ) &
1048	ALLOCATE ( t_surf_wall_v_2(l)%t(1:surf_usm_v(l)%ns) )
1049	IF ( .NOT. ALLOCATED( t_wall_v_1(l)%t ) ) &
1050	ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1051	IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) ) &
1052	ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1053	IF ( .NOT. ALLOCATED( t_surf_window_v_1(l)%t ) ) &
1054	ALLOCATE ( t_surf_window_v_1(l)%t(1:surf_usm_v(l)%ns) )
1055	IF ( .NOT. ALLOCATED( t_surf_window_v_2(l)%t ) ) &
1056	ALLOCATE ( t_surf_window_v_2(l)%t(1:surf_usm_v(l)%ns) )
1057	IF ( .NOT. ALLOCATED( t_window_v_1(l)%t ) ) &
1058	ALLOCATE ( t_window_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1059	IF ( .NOT. ALLOCATED( t_window_v_2(l)%t ) ) &
1060	ALLOCATE ( t_window_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1061	IF ( .NOT. ALLOCATED( t_surf_green_v_1(l)%t ) ) &
1062	ALLOCATE ( t_surf_green_v_1(l)%t(1:surf_usm_v(l)%ns) )
1063	IF ( .NOT. ALLOCATED( t_surf_green_v_2(l)%t ) ) &
1064	ALLOCATE ( t_surf_green_v_2(l)%t(1:surf_usm_v(l)%ns) )
1065	IF ( .NOT. ALLOCATED( t_green_v_1(l)%t ) ) &
1066	ALLOCATE ( t_green_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1067	IF ( .NOT. ALLOCATED( t_green_v_2(l)%t ) ) &
1068	ALLOCATE ( t_green_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1069	ENDDO
1070	!
1071	!-- initial assignment of the pointers
1072	t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2
1073	t_surf_wall_v => t_surf_wall_v_1; t_surf_wall_v_p => t_surf_wall_v_2
1074	t_window_v => t_window_v_1; t_window_v_p => t_window_v_2
1075	t_green_v => t_green_v_1; t_green_v_p => t_green_v_2
1076	t_surf_window_v => t_surf_window_v_1; t_surf_window_v_p => t_surf_window_v_2
1077	t_surf_green_v => t_surf_green_v_1; t_surf_green_v_p => t_surf_green_v_2
1078
1079	!
1080	!-- Allocate intermediate timestep arrays. For horizontal surfaces.
1081	ALLOCATE ( surf_usm_h%tt_surface_wall_m(1:surf_usm_h%ns) )
1082	ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1083	ALLOCATE ( surf_usm_h%tt_surface_window_m(1:surf_usm_h%ns) )
1084	ALLOCATE ( surf_usm_h%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1085	ALLOCATE ( surf_usm_h%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1086	ALLOCATE ( surf_usm_h%tt_surface_green_m(1:surf_usm_h%ns) )
1087
1088	!
1089	!-- Allocate intermediate timestep arrays
1090	!-- Horizontal surfaces
1091	ALLOCATE ( tm_liq_usm_h_m%var_usm_1d(1:surf_usm_h%ns) )
1092	tm_liq_usm_h_m%var_usm_1d = 0.0_wp
1093	!
1094	!-- Set inital values for prognostic quantities
1095	IF ( ALLOCATED( surf_usm_h%tt_surface_wall_m ) ) surf_usm_h%tt_surface_wall_m = 0.0_wp
1096	IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp
1097	IF ( ALLOCATED( surf_usm_h%tt_surface_window_m ) ) surf_usm_h%tt_surface_window_m = 0.0_wp
1098	IF ( ALLOCATED( surf_usm_h%tt_window_m ) ) surf_usm_h%tt_window_m = 0.0_wp
1099	IF ( ALLOCATED( surf_usm_h%tt_green_m ) ) surf_usm_h%tt_green_m = 0.0_wp
1100	IF ( ALLOCATED( surf_usm_h%tt_surface_green_m ) ) surf_usm_h%tt_surface_green_m = 0.0_wp
1101	!
1102	!-- Now, for vertical surfaces
1103	DO l = 0, 3
1104	ALLOCATE ( surf_usm_v(l)%tt_surface_wall_m(1:surf_usm_v(l)%ns) )
1105	ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1106	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_wall_m ) ) surf_usm_v(l)%tt_surface_wall_m = 0.0_wp
1107	IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp
1108	ALLOCATE ( surf_usm_v(l)%tt_surface_window_m(1:surf_usm_v(l)%ns) )
1109	ALLOCATE ( surf_usm_v(l)%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1110	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_window_m ) ) surf_usm_v(l)%tt_surface_window_m = 0.0_wp
1111	IF ( ALLOCATED( surf_usm_v(l)%tt_window_m ) ) surf_usm_v(l)%tt_window_m = 0.0_wp
1112	ALLOCATE ( surf_usm_v(l)%tt_surface_green_m(1:surf_usm_v(l)%ns) )
1113	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_green_m ) ) surf_usm_v(l)%tt_surface_green_m = 0.0_wp
1114	ALLOCATE ( surf_usm_v(l)%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1115	IF ( ALLOCATED( surf_usm_v(l)%tt_green_m ) ) surf_usm_v(l)%tt_green_m = 0.0_wp
1116	ENDDO
1117	!
1118	!-- allocate wall heat flux output array and set initial values. For horizontal surfaces
1119	! ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns) ) !can be removed
1120	ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
1121	ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
1122	ALLOCATE ( surf_usm_h%wghf_eb_window(1:surf_usm_h%ns) )
1123	ALLOCATE ( surf_usm_h%wghf_eb_green(1:surf_usm_h%ns) )
1124	ALLOCATE ( surf_usm_h%iwghf_eb(1:surf_usm_h%ns) )
1125	ALLOCATE ( surf_usm_h%iwghf_eb_window(1:surf_usm_h%ns) )
1126	IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp
1127	IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp
1128	IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp
1129	IF ( ALLOCATED( surf_usm_h%wghf_eb_window ) ) surf_usm_h%wghf_eb_window = 0.0_wp
1130	IF ( ALLOCATED( surf_usm_h%wghf_eb_green ) ) surf_usm_h%wghf_eb_green = 0.0_wp
1131	IF ( ALLOCATED( surf_usm_h%iwghf_eb ) ) surf_usm_h%iwghf_eb = 0.0_wp
1132	IF ( ALLOCATED( surf_usm_h%iwghf_eb_window ) ) surf_usm_h%iwghf_eb_window = 0.0_wp
1133	!
1134	!-- Now, for vertical surfaces
1135	DO l = 0, 3
1136	! ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns) ) ! can be removed
1137	ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
1138	ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
1139	ALLOCATE ( surf_usm_v(l)%wghf_eb_window(1:surf_usm_v(l)%ns) )
1140	ALLOCATE ( surf_usm_v(l)%wghf_eb_green(1:surf_usm_v(l)%ns) )
1141	ALLOCATE ( surf_usm_v(l)%iwghf_eb(1:surf_usm_v(l)%ns) )
1142	ALLOCATE ( surf_usm_v(l)%iwghf_eb_window(1:surf_usm_v(l)%ns) )
1143	IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp
1144	IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp
1145	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp
1146	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_window ) ) surf_usm_v(l)%wghf_eb_window = 0.0_wp
1147	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_green ) ) surf_usm_v(l)%wghf_eb_green = 0.0_wp
1148	IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb ) ) surf_usm_v(l)%iwghf_eb = 0.0_wp
1149	IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) ) surf_usm_v(l)%iwghf_eb_window = 0.0_wp
1150	ENDDO
1151	!
1152	!-- Initialize building-surface properties, which are also required by other modules,
1153	!-- e.g. the indoor model.
1154	CALL usm_define_pars
1155
1156	IF ( debug_output ) CALL debug_message( 'usm_init_arrays', 'end' )
1157
1158	END SUBROUTINE usm_init_arrays
1159
1160
1161	!------------------------------------------------------------------------------!
1162	! Description:
1163	! ------------
1164	!> Sum up and time-average urban surface output quantities as well as allocate
1165	!> the array necessary for storing the average.
1166	!------------------------------------------------------------------------------!
1167	SUBROUTINE usm_3d_data_averaging( mode, variable )
1168
1169	IMPLICIT NONE
1170
1171	CHARACTER(LEN=*), INTENT(IN) :: mode
1172	CHARACTER(LEN=*), INTENT(IN) :: variable
1173
1174	INTEGER(iwp) :: i, j, k, l, m, ids, idsint, iwl, istat !< runnin indices
1175	CHARACTER(LEN=varnamelength) :: var !< trimmed variable
1176	INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions
1177	CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
1178	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
1179
1180	IF ( variable(1:4) == 'usm_' ) THEN ! is such a check really rquired?
1181
1182	!
1183	!-- find the real name of the variable
1184	ids = -1
1185	l = -1
1186	var = TRIM(variable)
1187	DO i = 0, nd-1
1188	k = len(TRIM(var))
1189	j = len(TRIM(dirname(i)))
1190	IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN
1191	ids = i
1192	idsint = dirint(ids)
1193	var = var(:k-j)
1194	EXIT
1195	ENDIF
1196	ENDDO
1197	l = idsint - 2 ! horisontal direction index - terible hack !
1198	IF ( l < 0 .OR. l > 3 ) THEN
1199	l = -1
1200	END IF
1201	IF ( ids == -1 ) THEN
1202	var = TRIM(variable)
1203	ENDIF
1204	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
1205	!
1206	!-- wall layers
1207	READ(var(12:12), '(I1)', iostat=istat ) iwl
1208	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1209	var = var(1:10)
1210	ELSE
1211	!
1212	!-- wrong wall layer index
1213	RETURN
1214	ENDIF
1215	ENDIF
1216	IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN
1217	!
1218	!-- wall layers
1219	READ(var(14:14), '(I1)', iostat=istat ) iwl
1220	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1221	var = var(1:12)
1222	ELSE
1223	!
1224	!-- wrong window layer index
1225	RETURN
1226	ENDIF
1227	ENDIF
1228	IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN
1229	!
1230	!-- wall layers
1231	READ(var(13:13), '(I1)', iostat=istat ) iwl
1232	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1233	var = var(1:11)
1234	ELSE
1235	!
1236	!-- wrong green layer index
1237	RETURN
1238	ENDIF
1239	ENDIF
1240	IF ( var(1:8) == 'usm_swc_' .AND. len(TRIM(var)) >= 9 ) THEN
1241	!
1242	!-- swc layers
1243	READ(var(9:9), '(I1)', iostat=istat ) iwl
1244	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1245	var = var(1:7)
1246	ELSE
1247	!
1248	!-- wrong swc layer index
1249	RETURN
1250	ENDIF
1251	ENDIF
1252
1253	IF ( mode == 'allocate' ) THEN
1254
1255	SELECT CASE ( TRIM( var ) )
1256
1257	CASE ( 'usm_wshf' )
1258	!
1259	!-- array of sensible heat flux from surfaces
1260	!-- land surfaces
1261	IF ( l == -1 ) THEN
1262	IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN
1263	ALLOCATE ( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1264	surf_usm_h%wshf_eb_av = 0.0_wp
1265	ENDIF
1266	ELSE
1267	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN
1268	ALLOCATE ( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1269	surf_usm_v(l)%wshf_eb_av = 0.0_wp
1270	ENDIF
1271	ENDIF
1272
1273	CASE ( 'usm_qsws' )
1274	!
1275	!-- array of latent heat flux from surfaces
1276	!-- land surfaces
1277	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%qsws_av) ) THEN
1278	ALLOCATE ( surf_usm_h%qsws_av(1:surf_usm_h%ns) )
1279	surf_usm_h%qsws_av = 0.0_wp
1280	ELSE
1281	IF ( .NOT. ALLOCATED(surf_usm_v(l)%qsws_av) ) THEN
1282	ALLOCATE ( surf_usm_v(l)%qsws_av(1:surf_usm_v(l)%ns) )
1283	surf_usm_v(l)%qsws_av = 0.0_wp
1284	ENDIF
1285	ENDIF
1286
1287	CASE ( 'usm_qsws_veg' )
1288	!
1289	!-- array of latent heat flux from vegetation surfaces
1290	!-- land surfaces
1291	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%qsws_veg_av) ) THEN
1292	ALLOCATE ( surf_usm_h%qsws_veg_av(1:surf_usm_h%ns) )
1293	surf_usm_h%qsws_veg_av = 0.0_wp
1294	ELSE
1295	IF ( .NOT. ALLOCATED(surf_usm_v(l)%qsws_veg_av) ) THEN
1296	ALLOCATE ( surf_usm_v(l)%qsws_veg_av(1:surf_usm_v(l)%ns) )
1297	surf_usm_v(l)%qsws_veg_av = 0.0_wp
1298	ENDIF
1299	ENDIF
1300
1301	CASE ( 'usm_qsws_liq' )
1302	!
1303	!-- array of latent heat flux from surfaces with liquid
1304	!-- land surfaces
1305	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%qsws_liq_av) ) THEN
1306	ALLOCATE ( surf_usm_h%qsws_liq_av(1:surf_usm_h%ns) )
1307	surf_usm_h%qsws_liq_av = 0.0_wp
1308	ELSE
1309	IF ( .NOT. ALLOCATED(surf_usm_v(l)%qsws_liq_av) ) THEN
1310	ALLOCATE ( surf_usm_v(l)%qsws_liq_av(1:surf_usm_v(l)%ns) )
1311	surf_usm_v(l)%qsws_liq_av = 0.0_wp
1312	ENDIF
1313	ENDIF
1314	!
1315	!-- Please note, the following output quantities belongs to the
1316	!-- individual tile fractions - ground heat flux at wall-, window-,
1317	!-- and green fraction. Aggregated ground-heat flux is treated
1318	!-- accordingly in average_3d_data, sum_up_3d_data, etc..
1319	CASE ( 'usm_wghf' )
1320	!
1321	!-- array of heat flux from ground (wall, roof, land)
1322	IF ( l == -1 ) THEN
1323	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN
1324	ALLOCATE ( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1325	surf_usm_h%wghf_eb_av = 0.0_wp
1326	ENDIF
1327	ELSE
1328	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN
1329	ALLOCATE ( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1330	surf_usm_v(l)%wghf_eb_av = 0.0_wp
1331	ENDIF
1332	ENDIF
1333
1334	CASE ( 'usm_wghf_window' )
1335	!
1336	!-- array of heat flux from window ground (wall, roof, land)
1337	IF ( l == -1 ) THEN
1338	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN
1339	ALLOCATE ( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
1340	surf_usm_h%wghf_eb_window_av = 0.0_wp
1341	ENDIF
1342	ELSE
1343	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) ) THEN
1344	ALLOCATE ( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
1345	surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
1346	ENDIF
1347	ENDIF
1348
1349	CASE ( 'usm_wghf_green' )
1350	!
1351	!-- array of heat flux from green ground (wall, roof, land)
1352	IF ( l == -1 ) THEN
1353	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN
1354	ALLOCATE ( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
1355	surf_usm_h%wghf_eb_green_av = 0.0_wp
1356	ENDIF
1357	ELSE
1358	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) ) THEN
1359	ALLOCATE ( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
1360	surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
1361	ENDIF
1362	ENDIF
1363
1364	CASE ( 'usm_iwghf' )
1365	!
1366	!-- array of heat flux from indoor ground (wall, roof, land)
1367	IF ( l == -1 ) THEN
1368	IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN
1369	ALLOCATE ( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
1370	surf_usm_h%iwghf_eb_av = 0.0_wp
1371	ENDIF
1372	ELSE
1373	IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_av) ) THEN
1374	ALLOCATE ( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
1375	surf_usm_v(l)%iwghf_eb_av = 0.0_wp
1376	ENDIF
1377	ENDIF
1378
1379	CASE ( 'usm_iwghf_window' )
1380	!
1381	!-- array of heat flux from indoor window ground (wall, roof, land)
1382	IF ( l == -1 ) THEN
1383	IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN
1384	ALLOCATE ( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
1385	surf_usm_h%iwghf_eb_window_av = 0.0_wp
1386	ENDIF
1387	ELSE
1388	IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) ) THEN
1389	ALLOCATE ( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
1390	surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
1391	ENDIF
1392	ENDIF
1393
1394	CASE ( 'usm_t_surf_wall' )
1395	!
1396	!-- surface temperature for surfaces
1397	IF ( l == -1 ) THEN
1398	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_wall_av) ) THEN
1399	ALLOCATE ( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) )
1400	surf_usm_h%t_surf_wall_av = 0.0_wp
1401	ENDIF
1402	ELSE
1403	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_wall_av) ) THEN
1404	ALLOCATE ( surf_usm_v(l)%t_surf_wall_av(1:surf_usm_v(l)%ns) )
1405	surf_usm_v(l)%t_surf_wall_av = 0.0_wp
1406	ENDIF
1407	ENDIF
1408
1409	CASE ( 'usm_t_surf_window' )
1410	!
1411	!-- surface temperature for window surfaces
1412	IF ( l == -1 ) THEN
1413	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN
1414	ALLOCATE ( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
1415	surf_usm_h%t_surf_window_av = 0.0_wp
1416	ENDIF
1417	ELSE
1418	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_window_av) ) THEN
1419	ALLOCATE ( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
1420	surf_usm_v(l)%t_surf_window_av = 0.0_wp
1421	ENDIF
1422	ENDIF
1423
1424	CASE ( 'usm_t_surf_green' )
1425	!
1426	!-- surface temperature for green surfaces
1427	IF ( l == -1 ) THEN
1428	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN
1429	ALLOCATE ( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
1430	surf_usm_h%t_surf_green_av = 0.0_wp
1431	ENDIF
1432	ELSE
1433	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_green_av) ) THEN
1434	ALLOCATE ( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
1435	surf_usm_v(l)%t_surf_green_av = 0.0_wp
1436	ENDIF
1437	ENDIF
1438
1439	CASE ( 'usm_theta_10cm' )
1440	!
1441	!-- near surface (10cm) temperature for whole surfaces
1442	IF ( l == -1 ) THEN
1443	IF ( .NOT. ALLOCATED(surf_usm_h%pt_10cm_av) ) THEN
1444	ALLOCATE ( surf_usm_h%pt_10cm_av(1:surf_usm_h%ns) )
1445	surf_usm_h%pt_10cm_av = 0.0_wp
1446	ENDIF
1447	ELSE
1448	IF ( .NOT. ALLOCATED(surf_usm_v(l)%pt_10cm_av) ) THEN
1449	ALLOCATE ( surf_usm_v(l)%pt_10cm_av(1:surf_usm_v(l)%ns) )
1450	surf_usm_v(l)%pt_10cm_av = 0.0_wp
1451	ENDIF
1452	ENDIF
1453
1454	CASE ( 'usm_t_wall' )
1455	!
1456	!-- wall temperature for iwl layer of walls and land
1457	IF ( l == -1 ) THEN
1458	IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN
1459	ALLOCATE ( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1460	surf_usm_h%t_wall_av = 0.0_wp
1461	ENDIF
1462	ELSE
1463	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN
1464	ALLOCATE ( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1465	surf_usm_v(l)%t_wall_av = 0.0_wp
1466	ENDIF
1467	ENDIF
1468
1469	CASE ( 'usm_t_window' )
1470	!
1471	!-- window temperature for iwl layer of walls and land
1472	IF ( l == -1 ) THEN
1473	IF ( .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN
1474	ALLOCATE ( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1475	surf_usm_h%t_window_av = 0.0_wp
1476	ENDIF
1477	ELSE
1478	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_window_av) ) THEN
1479	ALLOCATE ( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1480	surf_usm_v(l)%t_window_av = 0.0_wp
1481	ENDIF
1482	ENDIF
1483
1484	CASE ( 'usm_t_green' )
1485	!
1486	!-- green temperature for iwl layer of walls and land
1487	IF ( l == -1 ) THEN
1488	IF ( .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN
1489	ALLOCATE ( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1490	surf_usm_h%t_green_av = 0.0_wp
1491	ENDIF
1492	ELSE
1493	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_green_av) ) THEN
1494	ALLOCATE ( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1495	surf_usm_v(l)%t_green_av = 0.0_wp
1496	ENDIF
1497	ENDIF
1498	CASE ( 'usm_swc' )
1499	!
1500	!-- soil water content for iwl layer of walls and land
1501	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%swc_av) ) THEN
1502	ALLOCATE ( surf_usm_h%swc_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1503	surf_usm_h%swc_av = 0.0_wp
1504	ELSE
1505	IF ( .NOT. ALLOCATED(surf_usm_v(l)%swc_av) ) THEN
1506	ALLOCATE ( surf_usm_v(l)%swc_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1507	surf_usm_v(l)%swc_av = 0.0_wp
1508	ENDIF
1509	ENDIF
1510
1511	CASE DEFAULT
1512	CONTINUE
1513
1514	END SELECT
1515
1516	ELSEIF ( mode == 'sum' ) THEN
1517
1518	SELECT CASE ( TRIM( var ) )
1519
1520	CASE ( 'usm_wshf' )
1521	!
1522	!-- array of sensible heat flux from surfaces (land, roof, wall)
1523	IF ( l == -1 ) THEN
1524	DO m = 1, surf_usm_h%ns
1525	surf_usm_h%wshf_eb_av(m) = &
1526	surf_usm_h%wshf_eb_av(m) + &
1527	surf_usm_h%wshf_eb(m)
1528	ENDDO
1529	ELSE
1530	DO m = 1, surf_usm_v(l)%ns
1531	surf_usm_v(l)%wshf_eb_av(m) = &
1532	surf_usm_v(l)%wshf_eb_av(m) + &
1533	surf_usm_v(l)%wshf_eb(m)
1534	ENDDO
1535	ENDIF
1536
1537	CASE ( 'usm_qsws' )
1538	!
1539	!-- array of latent heat flux from surfaces (land, roof, wall)
1540	IF ( l == -1 ) THEN
1541	DO m = 1, surf_usm_h%ns
1542	surf_usm_h%qsws_av(m) = &
1543	surf_usm_h%qsws_av(m) + &
1544	surf_usm_h%qsws(m) * l_v
1545	ENDDO
1546	ELSE
1547	DO m = 1, surf_usm_v(l)%ns
1548	surf_usm_v(l)%qsws_av(m) = &
1549	surf_usm_v(l)%qsws_av(m) + &
1550	surf_usm_v(l)%qsws(m) * l_v
1551	ENDDO
1552	ENDIF
1553
1554	CASE ( 'usm_qsws_veg' )
1555	!
1556	!-- array of latent heat flux from vegetation surfaces (land, roof, wall)
1557	IF ( l == -1 ) THEN
1558	DO m = 1, surf_usm_h%ns
1559	surf_usm_h%qsws_veg_av(m) = &
1560	surf_usm_h%qsws_veg_av(m) + &
1561	surf_usm_h%qsws_veg(m)
1562	ENDDO
1563	ELSE
1564	DO m = 1, surf_usm_v(l)%ns
1565	surf_usm_v(l)%qsws_veg_av(m) = &
1566	surf_usm_v(l)%qsws_veg_av(m) + &
1567	surf_usm_v(l)%qsws_veg(m)
1568	ENDDO
1569	ENDIF
1570
1571	CASE ( 'usm_qsws_liq' )
1572	!
1573	!-- array of latent heat flux from surfaces with liquid (land, roof, wall)
1574	IF ( l == -1 ) THEN
1575	DO m = 1, surf_usm_h%ns
1576	surf_usm_h%qsws_liq_av(m) = &
1577	surf_usm_h%qsws_liq_av(m) + &
1578	surf_usm_h%qsws_liq(m)
1579	ENDDO
1580	ELSE
1581	DO m = 1, surf_usm_v(l)%ns
1582	surf_usm_v(l)%qsws_liq_av(m) = &
1583	surf_usm_v(l)%qsws_liq_av(m) + &
1584	surf_usm_v(l)%qsws_liq(m)
1585	ENDDO
1586	ENDIF
1587
1588	CASE ( 'usm_wghf' )
1589	!
1590	!-- array of heat flux from ground (wall, roof, land)
1591	IF ( l == -1 ) THEN
1592	DO m = 1, surf_usm_h%ns
1593	surf_usm_h%wghf_eb_av(m) = &
1594	surf_usm_h%wghf_eb_av(m) + &
1595	surf_usm_h%wghf_eb(m)
1596	ENDDO
1597	ELSE
1598	DO m = 1, surf_usm_v(l)%ns
1599	surf_usm_v(l)%wghf_eb_av(m) = &
1600	surf_usm_v(l)%wghf_eb_av(m) + &
1601	surf_usm_v(l)%wghf_eb(m)
1602	ENDDO
1603	ENDIF
1604
1605	CASE ( 'usm_wghf_window' )
1606	!
1607	!-- array of heat flux from window ground (wall, roof, land)
1608	IF ( l == -1 ) THEN
1609	DO m = 1, surf_usm_h%ns
1610	surf_usm_h%wghf_eb_window_av(m) = &
1611	surf_usm_h%wghf_eb_window_av(m) + &
1612	surf_usm_h%wghf_eb_window(m)
1613	ENDDO
1614	ELSE
1615	DO m = 1, surf_usm_v(l)%ns
1616	surf_usm_v(l)%wghf_eb_window_av(m) = &
1617	surf_usm_v(l)%wghf_eb_window_av(m) + &
1618	surf_usm_v(l)%wghf_eb_window(m)
1619	ENDDO
1620	ENDIF
1621
1622	CASE ( 'usm_wghf_green' )
1623	!
1624	!-- array of heat flux from green ground (wall, roof, land)
1625	IF ( l == -1 ) THEN
1626	DO m = 1, surf_usm_h%ns
1627	surf_usm_h%wghf_eb_green_av(m) = &
1628	surf_usm_h%wghf_eb_green_av(m) + &
1629	surf_usm_h%wghf_eb_green(m)
1630	ENDDO
1631	ELSE
1632	DO m = 1, surf_usm_v(l)%ns
1633	surf_usm_v(l)%wghf_eb_green_av(m) = &
1634	surf_usm_v(l)%wghf_eb_green_av(m) + &
1635	surf_usm_v(l)%wghf_eb_green(m)
1636	ENDDO
1637	ENDIF
1638
1639	CASE ( 'usm_iwghf' )
1640	!
1641	!-- array of heat flux from indoor ground (wall, roof, land)
1642	IF ( l == -1 ) THEN
1643	DO m = 1, surf_usm_h%ns
1644	surf_usm_h%iwghf_eb_av(m) = &
1645	surf_usm_h%iwghf_eb_av(m) + &
1646	surf_usm_h%iwghf_eb(m)
1647	ENDDO
1648	ELSE
1649	DO m = 1, surf_usm_v(l)%ns
1650	surf_usm_v(l)%iwghf_eb_av(m) = &
1651	surf_usm_v(l)%iwghf_eb_av(m) + &
1652	surf_usm_v(l)%iwghf_eb(m)
1653	ENDDO
1654	ENDIF
1655
1656	CASE ( 'usm_iwghf_window' )
1657	!
1658	!-- array of heat flux from indoor window ground (wall, roof, land)
1659	IF ( l == -1 ) THEN
1660	DO m = 1, surf_usm_h%ns
1661	surf_usm_h%iwghf_eb_window_av(m) = &
1662	surf_usm_h%iwghf_eb_window_av(m) + &
1663	surf_usm_h%iwghf_eb_window(m)
1664	ENDDO
1665	ELSE
1666	DO m = 1, surf_usm_v(l)%ns
1667	surf_usm_v(l)%iwghf_eb_window_av(m) = &
1668	surf_usm_v(l)%iwghf_eb_window_av(m) + &
1669	surf_usm_v(l)%iwghf_eb_window(m)
1670	ENDDO
1671	ENDIF
1672
1673	CASE ( 'usm_t_surf_wall' )
1674	!
1675	!-- surface temperature for surfaces
1676	IF ( l == -1 ) THEN
1677	DO m = 1, surf_usm_h%ns
1678	surf_usm_h%t_surf_wall_av(m) = &
1679	surf_usm_h%t_surf_wall_av(m) + &
1680	t_surf_wall_h(m)
1681	ENDDO
1682	ELSE
1683	DO m = 1, surf_usm_v(l)%ns
1684	surf_usm_v(l)%t_surf_wall_av(m) = &
1685	surf_usm_v(l)%t_surf_wall_av(m) + &
1686	t_surf_wall_v(l)%t(m)
1687	ENDDO
1688	ENDIF
1689
1690	CASE ( 'usm_t_surf_window' )
1691	!
1692	!-- surface temperature for window surfaces
1693	IF ( l == -1 ) THEN
1694	DO m = 1, surf_usm_h%ns
1695	surf_usm_h%t_surf_window_av(m) = &
1696	surf_usm_h%t_surf_window_av(m) + &
1697	t_surf_window_h(m)
1698	ENDDO
1699	ELSE
1700	DO m = 1, surf_usm_v(l)%ns
1701	surf_usm_v(l)%t_surf_window_av(m) = &
1702	surf_usm_v(l)%t_surf_window_av(m) + &
1703	t_surf_window_v(l)%t(m)
1704	ENDDO
1705	ENDIF
1706
1707	CASE ( 'usm_t_surf_green' )
1708	!
1709	!-- surface temperature for green surfaces
1710	IF ( l == -1 ) THEN
1711	DO m = 1, surf_usm_h%ns
1712	surf_usm_h%t_surf_green_av(m) = &
1713	surf_usm_h%t_surf_green_av(m) + &
1714	t_surf_green_h(m)
1715	ENDDO
1716	ELSE
1717	DO m = 1, surf_usm_v(l)%ns
1718	surf_usm_v(l)%t_surf_green_av(m) = &
1719	surf_usm_v(l)%t_surf_green_av(m) + &
1720	t_surf_green_v(l)%t(m)
1721	ENDDO
1722	ENDIF
1723
1724	CASE ( 'usm_theta_10cm' )
1725	!
1726	!-- near surface temperature for whole surfaces
1727	IF ( l == -1 ) THEN
1728	DO m = 1, surf_usm_h%ns
1729	surf_usm_h%pt_10cm_av(m) = &
1730	surf_usm_h%pt_10cm_av(m) + &
1731	surf_usm_h%pt_10cm(m)
1732	ENDDO
1733	ELSE
1734	DO m = 1, surf_usm_v(l)%ns
1735	surf_usm_v(l)%pt_10cm_av(m) = &
1736	surf_usm_v(l)%pt_10cm_av(m) + &
1737	surf_usm_v(l)%pt_10cm(m)
1738	ENDDO
1739	ENDIF
1740
1741	CASE ( 'usm_t_wall' )
1742	!
1743	!-- wall temperature for iwl layer of walls and land
1744	IF ( l == -1 ) THEN
1745	DO m = 1, surf_usm_h%ns
1746	surf_usm_h%t_wall_av(iwl,m) = &
1747	surf_usm_h%t_wall_av(iwl,m) + &
1748	t_wall_h(iwl,m)
1749	ENDDO
1750	ELSE
1751	DO m = 1, surf_usm_v(l)%ns
1752	surf_usm_v(l)%t_wall_av(iwl,m) = &
1753	surf_usm_v(l)%t_wall_av(iwl,m) + &
1754	t_wall_v(l)%t(iwl,m)
1755	ENDDO
1756	ENDIF
1757
1758	CASE ( 'usm_t_window' )
1759	!
1760	!-- window temperature for iwl layer of walls and land
1761	IF ( l == -1 ) THEN
1762	DO m = 1, surf_usm_h%ns
1763	surf_usm_h%t_window_av(iwl,m) = &
1764	surf_usm_h%t_window_av(iwl,m) + &
1765	t_window_h(iwl,m)
1766	ENDDO
1767	ELSE
1768	DO m = 1, surf_usm_v(l)%ns
1769	surf_usm_v(l)%t_window_av(iwl,m) = &
1770	surf_usm_v(l)%t_window_av(iwl,m) + &
1771	t_window_v(l)%t(iwl,m)
1772	ENDDO
1773	ENDIF
1774
1775	CASE ( 'usm_t_green' )
1776	!
1777	!-- green temperature for iwl layer of walls and land
1778	IF ( l == -1 ) THEN
1779	DO m = 1, surf_usm_h%ns
1780	surf_usm_h%t_green_av(iwl,m) = &
1781	surf_usm_h%t_green_av(iwl,m) + &
1782	t_green_h(iwl,m)
1783	ENDDO
1784	ELSE
1785	DO m = 1, surf_usm_v(l)%ns
1786	surf_usm_v(l)%t_green_av(iwl,m) = &
1787	surf_usm_v(l)%t_green_av(iwl,m) + &
1788	t_green_v(l)%t(iwl,m)
1789	ENDDO
1790	ENDIF
1791
1792	CASE ( 'usm_swc' )
1793	!
1794	!-- soil water content for iwl layer of walls and land
1795	IF ( l == -1 ) THEN
1796	DO m = 1, surf_usm_h%ns
1797	surf_usm_h%swc_av(iwl,m) = &
1798	surf_usm_h%swc_av(iwl,m) + &
1799	swc_h(iwl,m)
1800	ENDDO
1801	ELSE
1802	ENDIF
1803
1804	CASE DEFAULT
1805	CONTINUE
1806
1807	END SELECT
1808
1809	ELSEIF ( mode == 'average' ) THEN
1810
1811	SELECT CASE ( TRIM( var ) )
1812
1813	CASE ( 'usm_wshf' )
1814	!
1815	!-- array of sensible heat flux from surfaces (land, roof, wall)
1816	IF ( l == -1 ) THEN
1817	DO m = 1, surf_usm_h%ns
1818	surf_usm_h%wshf_eb_av(m) = &
1819	surf_usm_h%wshf_eb_av(m) / &
1820	REAL( average_count_3d, kind=wp )
1821	ENDDO
1822	ELSE
1823	DO m = 1, surf_usm_v(l)%ns
1824	surf_usm_v(l)%wshf_eb_av(m) = &
1825	surf_usm_v(l)%wshf_eb_av(m) / &
1826	REAL( average_count_3d, kind=wp )
1827	ENDDO
1828	ENDIF
1829
1830	CASE ( 'usm_qsws' )
1831	!
1832	!-- array of latent heat flux from surfaces (land, roof, wall)
1833	IF ( l == -1 ) THEN
1834	DO m = 1, surf_usm_h%ns
1835	surf_usm_h%qsws_av(m) = &
1836	surf_usm_h%qsws_av(m) / &
1837	REAL( average_count_3d, kind=wp )
1838	ENDDO
1839	ELSE
1840	DO m = 1, surf_usm_v(l)%ns
1841	surf_usm_v(l)%qsws_av(m) = &
1842	surf_usm_v(l)%qsws_av(m) / &
1843	REAL( average_count_3d, kind=wp )
1844	ENDDO
1845	ENDIF
1846
1847	CASE ( 'usm_qsws_veg' )
1848	!
1849	!-- array of latent heat flux from vegetation surfaces (land, roof, wall)
1850	IF ( l == -1 ) THEN
1851	DO m = 1, surf_usm_h%ns
1852	surf_usm_h%qsws_veg_av(m) = &
1853	surf_usm_h%qsws_veg_av(m) / &
1854	REAL( average_count_3d, kind=wp )
1855	ENDDO
1856	ELSE
1857	DO m = 1, surf_usm_v(l)%ns
1858	surf_usm_v(l)%qsws_veg_av(m) = &
1859	surf_usm_v(l)%qsws_veg_av(m) / &
1860	REAL( average_count_3d, kind=wp )
1861	ENDDO
1862	ENDIF
1863
1864	CASE ( 'usm_qsws_liq' )
1865	!
1866	!-- array of latent heat flux from surfaces with liquid (land, roof, wall)
1867	IF ( l == -1 ) THEN
1868	DO m = 1, surf_usm_h%ns
1869	surf_usm_h%qsws_liq_av(m) = &
1870	surf_usm_h%qsws_liq_av(m) / &
1871	REAL( average_count_3d, kind=wp )
1872	ENDDO
1873	ELSE
1874	DO m = 1, surf_usm_v(l)%ns
1875	surf_usm_v(l)%qsws_liq_av(m) = &
1876	surf_usm_v(l)%qsws_liq_av(m) / &
1877	REAL( average_count_3d, kind=wp )
1878	ENDDO
1879	ENDIF
1880
1881	CASE ( 'usm_wghf' )
1882	!
1883	!-- array of heat flux from ground (wall, roof, land)
1884	IF ( l == -1 ) THEN
1885	DO m = 1, surf_usm_h%ns
1886	surf_usm_h%wghf_eb_av(m) = &
1887	surf_usm_h%wghf_eb_av(m) / &
1888	REAL( average_count_3d, kind=wp )
1889	ENDDO
1890	ELSE
1891	DO m = 1, surf_usm_v(l)%ns
1892	surf_usm_v(l)%wghf_eb_av(m) = &
1893	surf_usm_v(l)%wghf_eb_av(m) / &
1894	REAL( average_count_3d, kind=wp )
1895	ENDDO
1896	ENDIF
1897
1898	CASE ( 'usm_wghf_window' )
1899	!
1900	!-- array of heat flux from window ground (wall, roof, land)
1901	IF ( l == -1 ) THEN
1902	DO m = 1, surf_usm_h%ns
1903	surf_usm_h%wghf_eb_window_av(m) = &
1904	surf_usm_h%wghf_eb_window_av(m) / &
1905	REAL( average_count_3d, kind=wp )
1906	ENDDO
1907	ELSE
1908	DO m = 1, surf_usm_v(l)%ns
1909	surf_usm_v(l)%wghf_eb_window_av(m) = &
1910	surf_usm_v(l)%wghf_eb_window_av(m) / &
1911	REAL( average_count_3d, kind=wp )
1912	ENDDO
1913	ENDIF
1914
1915	CASE ( 'usm_wghf_green' )
1916	!
1917	!-- array of heat flux from green ground (wall, roof, land)
1918	IF ( l == -1 ) THEN
1919	DO m = 1, surf_usm_h%ns
1920	surf_usm_h%wghf_eb_green_av(m) = &
1921	surf_usm_h%wghf_eb_green_av(m) / &
1922	REAL( average_count_3d, kind=wp )
1923	ENDDO
1924	ELSE
1925	DO m = 1, surf_usm_v(l)%ns
1926	surf_usm_v(l)%wghf_eb_green_av(m) = &
1927	surf_usm_v(l)%wghf_eb_green_av(m) / &
1928	REAL( average_count_3d, kind=wp )
1929	ENDDO
1930	ENDIF
1931
1932	CASE ( 'usm_iwghf' )
1933	!
1934	!-- array of heat flux from indoor ground (wall, roof, land)
1935	IF ( l == -1 ) THEN
1936	DO m = 1, surf_usm_h%ns
1937	surf_usm_h%iwghf_eb_av(m) = &
1938	surf_usm_h%iwghf_eb_av(m) / &
1939	REAL( average_count_3d, kind=wp )
1940	ENDDO
1941	ELSE
1942	DO m = 1, surf_usm_v(l)%ns
1943	surf_usm_v(l)%iwghf_eb_av(m) = &
1944	surf_usm_v(l)%iwghf_eb_av(m) / &
1945	REAL( average_count_3d, kind=wp )
1946	ENDDO
1947	ENDIF
1948
1949	CASE ( 'usm_iwghf_window' )
1950	!
1951	!-- array of heat flux from indoor window ground (wall, roof, land)
1952	IF ( l == -1 ) THEN
1953	DO m = 1, surf_usm_h%ns
1954	surf_usm_h%iwghf_eb_window_av(m) = &
1955	surf_usm_h%iwghf_eb_window_av(m) / &
1956	REAL( average_count_3d, kind=wp )
1957	ENDDO
1958	ELSE
1959	DO m = 1, surf_usm_v(l)%ns
1960	surf_usm_v(l)%iwghf_eb_window_av(m) = &
1961	surf_usm_v(l)%iwghf_eb_window_av(m) / &
1962	REAL( average_count_3d, kind=wp )
1963	ENDDO
1964	ENDIF
1965
1966	CASE ( 'usm_t_surf_wall' )
1967	!
1968	!-- surface temperature for surfaces
1969	IF ( l == -1 ) THEN
1970	DO m = 1, surf_usm_h%ns
1971	surf_usm_h%t_surf_wall_av(m) = &
1972	surf_usm_h%t_surf_wall_av(m) / &
1973	REAL( average_count_3d, kind=wp )
1974	ENDDO
1975	ELSE
1976	DO m = 1, surf_usm_v(l)%ns
1977	surf_usm_v(l)%t_surf_wall_av(m) = &
1978	surf_usm_v(l)%t_surf_wall_av(m) / &
1979	REAL( average_count_3d, kind=wp )
1980	ENDDO
1981	ENDIF
1982
1983	CASE ( 'usm_t_surf_window' )
1984	!
1985	!-- surface temperature for window surfaces
1986	IF ( l == -1 ) THEN
1987	DO m = 1, surf_usm_h%ns
1988	surf_usm_h%t_surf_window_av(m) = &
1989	surf_usm_h%t_surf_window_av(m) / &
1990	REAL( average_count_3d, kind=wp )
1991	ENDDO
1992	ELSE
1993	DO m = 1, surf_usm_v(l)%ns
1994	surf_usm_v(l)%t_surf_window_av(m) = &
1995	surf_usm_v(l)%t_surf_window_av(m) / &
1996	REAL( average_count_3d, kind=wp )
1997	ENDDO
1998	ENDIF
1999
2000	CASE ( 'usm_t_surf_green' )
2001	!
2002	!-- surface temperature for green surfaces
2003	IF ( l == -1 ) THEN
2004	DO m = 1, surf_usm_h%ns
2005	surf_usm_h%t_surf_green_av(m) = &
2006	surf_usm_h%t_surf_green_av(m) / &
2007	REAL( average_count_3d, kind=wp )
2008	ENDDO
2009	ELSE
2010	DO m = 1, surf_usm_v(l)%ns
2011	surf_usm_v(l)%t_surf_green_av(m) = &
2012	surf_usm_v(l)%t_surf_green_av(m) / &
2013	REAL( average_count_3d, kind=wp )
2014	ENDDO
2015	ENDIF
2016
2017	CASE ( 'usm_theta_10cm' )
2018	!
2019	!-- near surface temperature for whole surfaces
2020	IF ( l == -1 ) THEN
2021	DO m = 1, surf_usm_h%ns
2022	surf_usm_h%pt_10cm_av(m) = &
2023	surf_usm_h%pt_10cm_av(m) / &
2024	REAL( average_count_3d, kind=wp )
2025	ENDDO
2026	ELSE
2027	DO m = 1, surf_usm_v(l)%ns
2028	surf_usm_v(l)%pt_10cm_av(m) = &
2029	surf_usm_v(l)%pt_10cm_av(m) / &
2030	REAL( average_count_3d, kind=wp )
2031	ENDDO
2032	ENDIF
2033
2034
2035	CASE ( 'usm_t_wall' )
2036	!
2037	!-- wall temperature for iwl layer of walls and land
2038	IF ( l == -1 ) THEN
2039	DO m = 1, surf_usm_h%ns
2040	surf_usm_h%t_wall_av(iwl,m) = &
2041	surf_usm_h%t_wall_av(iwl,m) / &
2042	REAL( average_count_3d, kind=wp )
2043	ENDDO
2044	ELSE
2045	DO m = 1, surf_usm_v(l)%ns
2046	surf_usm_v(l)%t_wall_av(iwl,m) = &
2047	surf_usm_v(l)%t_wall_av(iwl,m) / &
2048	REAL( average_count_3d, kind=wp )
2049	ENDDO
2050	ENDIF
2051
2052	CASE ( 'usm_t_window' )
2053	!
2054	!-- window temperature for iwl layer of walls and land
2055	IF ( l == -1 ) THEN
2056	DO m = 1, surf_usm_h%ns
2057	surf_usm_h%t_window_av(iwl,m) = &
2058	surf_usm_h%t_window_av(iwl,m) / &
2059	REAL( average_count_3d, kind=wp )
2060	ENDDO
2061	ELSE
2062	DO m = 1, surf_usm_v(l)%ns
2063	surf_usm_v(l)%t_window_av(iwl,m) = &
2064	surf_usm_v(l)%t_window_av(iwl,m) / &
2065	REAL( average_count_3d, kind=wp )
2066	ENDDO
2067	ENDIF
2068
2069	CASE ( 'usm_t_green' )
2070	!
2071	!-- green temperature for iwl layer of walls and land
2072	IF ( l == -1 ) THEN
2073	DO m = 1, surf_usm_h%ns
2074	surf_usm_h%t_green_av(iwl,m) = &
2075	surf_usm_h%t_green_av(iwl,m) / &
2076	REAL( average_count_3d, kind=wp )
2077	ENDDO
2078	ELSE
2079	DO m = 1, surf_usm_v(l)%ns
2080	surf_usm_v(l)%t_green_av(iwl,m) = &
2081	surf_usm_v(l)%t_green_av(iwl,m) / &
2082	REAL( average_count_3d, kind=wp )
2083	ENDDO
2084	ENDIF
2085
2086	CASE ( 'usm_swc' )
2087	!
2088	!-- soil water content for iwl layer of walls and land
2089	IF ( l == -1 ) THEN
2090	DO m = 1, surf_usm_h%ns
2091	surf_usm_h%swc_av(iwl,m) = &
2092	surf_usm_h%swc_av(iwl,m) / &
2093	REAL( average_count_3d, kind=wp )
2094	ENDDO
2095	ELSE
2096	DO m = 1, surf_usm_v(l)%ns
2097	surf_usm_v(l)%swc_av(iwl,m) = &
2098	surf_usm_v(l)%swc_av(iwl,m) / &
2099	REAL( average_count_3d, kind=wp )
2100	ENDDO
2101	ENDIF
2102
2103
2104	END SELECT
2105
2106	ENDIF
2107
2108	ENDIF
2109
2110	END SUBROUTINE usm_3d_data_averaging
2111
2112
2113
2114	!------------------------------------------------------------------------------!
2115	! Description:
2116	! ------------
2117	!> Set internal Neumann boundary condition at outer soil grid points
2118	!> for temperature and humidity.
2119	!------------------------------------------------------------------------------!
2120	SUBROUTINE usm_boundary_condition
2121
2122	IMPLICIT NONE
2123
2124	INTEGER(iwp) :: i !< grid index x-direction
2125	INTEGER(iwp) :: ioff !< offset index x-direction indicating location of soil grid point
2126	INTEGER(iwp) :: j !< grid index y-direction
2127	INTEGER(iwp) :: joff !< offset index x-direction indicating location of soil grid point
2128	INTEGER(iwp) :: k !< grid index z-direction
2129	INTEGER(iwp) :: koff !< offset index x-direction indicating location of soil grid point
2130	INTEGER(iwp) :: l !< running index surface-orientation
2131	INTEGER(iwp) :: m !< running index surface elements
2132
2133	koff = surf_usm_h%koff
2134	DO m = 1, surf_usm_h%ns
2135	i = surf_usm_h%i(m)
2136	j = surf_usm_h%j(m)
2137	k = surf_usm_h%k(m)
2138	pt(k+koff,j,i) = pt(k,j,i)
2139	ENDDO
2140
2141	DO l = 0, 3
2142	ioff = surf_usm_v(l)%ioff
2143	joff = surf_usm_v(l)%joff
2144	DO m = 1, surf_usm_v(l)%ns
2145	i = surf_usm_v(l)%i(m)
2146	j = surf_usm_v(l)%j(m)
2147	k = surf_usm_v(l)%k(m)
2148	pt(k,j+joff,i+ioff) = pt(k,j,i)
2149	ENDDO
2150	ENDDO
2151
2152	END SUBROUTINE usm_boundary_condition
2153
2154
2155	!------------------------------------------------------------------------------!
2156	!
2157	! Description:
2158	! ------------
2159	!> Subroutine checks variables and assigns units.
2160	!> It is called out from subroutine check_parameters.
2161	!------------------------------------------------------------------------------!
2162	SUBROUTINE usm_check_data_output( variable, unit )
2163
2164	IMPLICIT NONE
2165
2166	CHARACTER(LEN=*),INTENT(IN) :: variable !<
2167	CHARACTER(LEN=*),INTENT(OUT) :: unit !<
2168
2169	INTEGER(iwp) :: i,j,l !< index
2170	CHARACTER(LEN=2) :: ls
2171	CHARACTER(LEN=varnamelength) :: var !< TRIM(variable)
2172	INTEGER(iwp), PARAMETER :: nl1 = 15 !< number of directional usm variables
2173	CHARACTER(LEN=varnamelength), DIMENSION(nl1) :: varlist1 = & !< list of directional usm variables
2174	(/'usm_wshf ', &
2175	'usm_wghf ', &
2176	'usm_wghf_window ', &
2177	'usm_wghf_green ', &
2178	'usm_iwghf ', &
2179	'usm_iwghf_window ', &
2180	'usm_surfz ', &
2181	'usm_surfwintrans ', &
2182	'usm_surfcat ', &
2183	'usm_t_surf_wall ', &
2184	'usm_t_surf_window ', &
2185	'usm_t_surf_green ', &
2186	'usm_t_green ', &
2187	'usm_qsws ', &
2188	'usm_theta_10cm '/)
2189
2190	INTEGER(iwp), PARAMETER :: nl2 = 3 !< number of directional layer usm variables
2191	CHARACTER(LEN=varnamelength), DIMENSION(nl2) :: varlist2 = & !< list of directional layer usm variables
2192	(/'usm_t_wall ', &
2193	'usm_t_window ', &
2194	'usm_t_green '/)
2195
2196	INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions
2197	CHARACTER(LEN=6), DIMENSION(nd), PARAMETER :: dirname = & !< direction names
2198	(/'_roof ','_south','_north','_west ','_east '/)
2199	LOGICAL :: lfound !< flag if the variable is found
2200
2201
2202	lfound = .FALSE.
2203
2204	var = TRIM(variable)
2205
2206	!
2207	!-- check if variable exists
2208	!-- directional variables
2209	DO i = 1, nl1
2210	DO j = 1, nd
2211	IF ( TRIM(var) == TRIM(varlist1(i))//TRIM(dirname(j)) ) THEN
2212	lfound = .TRUE.
2213	EXIT
2214	ENDIF
2215	IF ( lfound ) EXIT
2216	ENDDO
2217	ENDDO
2218	IF ( lfound ) GOTO 10
2219	!
2220	!-- directional layer variables
2221	DO i = 1, nl2
2222	DO j = 1, nd
2223	DO l = nzb_wall, nzt_wall
2224	WRITE(ls,'(A1,I1)') '_',l
2225	IF ( TRIM(var) == TRIM(varlist2(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
2226	lfound = .TRUE.
2227	EXIT
2228	ENDIF
2229	ENDDO
2230	IF ( lfound ) EXIT
2231	ENDDO
2232	ENDDO
2233	IF ( .NOT. lfound ) THEN
2234	unit = 'illegal'
2235	RETURN
2236	ENDIF
2237	10 CONTINUE
2238
2239	IF ( var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
2240	var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. &
2241	var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. &
2242	var(1:17) == 'usm_surfwintrans_' .OR. &
2243	var(1:9) == 'usm_qsws_' .OR. var(1:13) == 'usm_qsws_veg_' .OR. &
2244	var(1:13) == 'usm_qsws_liq_' ) THEN
2245	unit = 'W/m2'
2246	ELSE IF ( var(1:15) == 'usm_t_surf_wall' .OR. var(1:10) == 'usm_t_wall' .OR. &
2247	var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR. &
2248	var(1:16) == 'usm_t_surf_green' .OR. &
2249	var(1:11) == 'usm_t_green' .OR. var(1:7) == 'usm_swc' .OR. &
2250	var(1:14) == 'usm_theta_10cm' ) THEN
2251	unit = 'K'
2252	ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:11) == 'usm_surfcat' ) THEN
2253	unit = '1'
2254	ELSE
2255	unit = 'illegal'
2256	ENDIF
2257
2258	END SUBROUTINE usm_check_data_output
2259
2260
2261	!------------------------------------------------------------------------------!
2262	! Description:
2263	! ------------
2264	!> Check parameters routine for urban surface model
2265	!------------------------------------------------------------------------------!
2266	SUBROUTINE usm_check_parameters
2267
2268	USE control_parameters, &
2269	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing, &
2270	lsf_surf, topography
2271
2272	USE netcdf_data_input_mod, &
2273	ONLY: building_type_f
2274
2275	IMPLICIT NONE
2276
2277	INTEGER(iwp) :: i !< running index, x-dimension
2278	INTEGER(iwp) :: j !< running index, y-dimension
2279
2280	!
2281	!-- Dirichlet boundary conditions are required as the surface fluxes are
2282	!-- calculated from the temperature/humidity gradients in the urban surface
2283	!-- model
2284	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
2285	message_string = 'urban surface model requires setting of '// &
2286	'bc_pt_b = "dirichlet" and '// &
2287	'bc_q_b = "dirichlet"'
2288	CALL message( 'usm_check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
2289	ENDIF
2290
2291	IF ( .NOT. constant_flux_layer ) THEN
2292	message_string = 'urban surface model requires '// &
2293	'constant_flux_layer = .T.'
2294	CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
2295	ENDIF
2296
2297	IF ( .NOT. radiation ) THEN
2298	message_string = 'urban surface model requires '// &
2299	'the radiation model to be switched on'
2300	CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
2301	ENDIF
2302	!
2303	!-- Surface forcing has to be disabled for LSF in case of enabled
2304	!-- urban surface module
2305	IF ( large_scale_forcing ) THEN
2306	lsf_surf = .FALSE.
2307	ENDIF
2308	!
2309	!-- Topography
2310	IF ( topography == 'flat' ) THEN
2311	message_string = 'topography /= "flat" is required '// &
2312	'when using the urban surface model'
2313	CALL message( 'usm_check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2314	ENDIF
2315	!
2316	!-- naheatlayers
2317	IF ( naheatlayers > nzt ) THEN
2318	message_string = 'number of anthropogenic heat layers '// &
2319	'"naheatlayers" can not be larger than'// &
2320	' number of domain layers "nzt"'
2321	CALL message( 'usm_check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
2322	ENDIF
2323	!
2324	!-- Check if building types are set within a valid range.
2325	IF ( building_type < LBOUND( building_pars, 2 ) .AND. &
2326	building_type > UBOUND( building_pars, 2 ) ) THEN
2327	WRITE( message_string, * ) 'building_type = ', building_type, &
2328	' is out of the valid range'
2329	CALL message( 'usm_check_parameters', 'PA0529', 2, 2, 0, 6, 0 )
2330	ENDIF
2331	IF ( building_type_f%from_file ) THEN
2332	DO i = nxl, nxr
2333	DO j = nys, nyn
2334	IF ( building_type_f%var(j,i) /= building_type_f%fill .AND. &
2335	( building_type_f%var(j,i) < LBOUND( building_pars, 2 ) .OR. &
2336	building_type_f%var(j,i) > UBOUND( building_pars, 2 ) ) ) &
2337	THEN
2338	WRITE( message_string, * ) 'building_type = is out of ' // &
2339	'the valid range at (j,i) = ', j, i
2340	CALL message( 'usm_check_parameters', 'PA0529', 2, 2, myid, 6, 0 )
2341	ENDIF
2342	ENDDO
2343	ENDDO
2344	ENDIF
2345	END SUBROUTINE usm_check_parameters
2346
2347
2348	!------------------------------------------------------------------------------!
2349	!
2350	! Description:
2351	! ------------
2352	!> Output of the 3D-arrays in netCDF and/or AVS format
2353	!> for variables of urban_surface model.
2354	!> It resorts the urban surface module output quantities from surf style
2355	!> indexing into temporary 3D array with indices (i,j,k).
2356	!> It is called from subroutine data_output_3d.
2357	!------------------------------------------------------------------------------!
2358	SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2359
2360	IMPLICIT NONE
2361
2362	INTEGER(iwp), INTENT(IN) :: av !< flag if averaged
2363	CHARACTER (len=*), INTENT(IN) :: variable !< variable name
2364	INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0)
2365	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
2366	LOGICAL, INTENT(OUT) :: found !<
2367	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< sp - it has to correspond to module data_output_3d
2368	REAL(sp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: temp_pf !< temp array for urban surface output procedure
2369
2370	CHARACTER (len=varnamelength) :: var !< trimmed variable name
2371	INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions
2372	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
2373	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
2374	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: diridx = (/ -1, 1, 0, 3, 2 /)
2375	!< index for surf_*_v: 0:3 = (North, South, East, West)
2376	INTEGER(iwp) :: ids,idsint,idsidx
2377	INTEGER(iwp) :: i,j,k,iwl,istat, l, m !< running indices
2378
2379	found = .TRUE.
2380	temp_pf = -1._wp
2381
2382	ids = -1
2383	var = TRIM(variable)
2384	DO i = 0, nd-1
2385	k = len(TRIM(var))
2386	j = len(TRIM(dirname(i)))
2387	IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN
2388	ids = i
2389	idsint = dirint(ids)
2390	idsidx = diridx(ids)
2391	var = var(:k-j)
2392	EXIT
2393	ENDIF
2394	ENDDO
2395	IF ( ids == -1 ) THEN
2396	var = TRIM(variable)
2397	ENDIF
2398	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
2399	!
2400	!-- wall layers
2401	READ(var(12:12), '(I1)', iostat=istat ) iwl
2402	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2403	var = var(1:10)
2404	ENDIF
2405	ENDIF
2406	IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN
2407	!
2408	!-- window layers
2409	READ(var(14:14), '(I1)', iostat=istat ) iwl
2410	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2411	var = var(1:12)
2412	ENDIF
2413	ENDIF
2414	IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN
2415	!
2416	!-- green layers
2417	READ(var(13:13), '(I1)', iostat=istat ) iwl
2418	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2419	var = var(1:11)
2420	ENDIF
2421	ENDIF
2422	IF ( var(1:8) == 'usm_swc_' .AND. len(TRIM(var)) >= 9 ) THEN
2423	!
2424	!-- green layers soil water content
2425	READ(var(9:9), '(I1)', iostat=istat ) iwl
2426	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2427	var = var(1:7)
2428	ENDIF
2429	ENDIF
2430
2431	SELECT CASE ( TRIM(var) )
2432
2433	CASE ( 'usm_surfz' )
2434	!
2435	!-- array of surface height (z)
2436	IF ( idsint == iup_u ) THEN
2437	DO m = 1, surf_usm_h%ns
2438	i = surf_usm_h%i(m)
2439	j = surf_usm_h%j(m)
2440	k = surf_usm_h%k(m)
2441	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, KIND = sp) )
2442	ENDDO
2443	ELSE
2444	l = idsidx
2445	DO m = 1, surf_usm_v(l)%ns
2446	i = surf_usm_v(l)%i(m)
2447	j = surf_usm_v(l)%j(m)
2448	k = surf_usm_v(l)%k(m)
2449	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, KIND = sp) + 1.0_sp )
2450	ENDDO
2451	ENDIF
2452
2453	CASE ( 'usm_surfcat' )
2454	!
2455	!-- surface category
2456	IF ( idsint == iup_u ) THEN
2457	DO m = 1, surf_usm_h%ns
2458	i = surf_usm_h%i(m)
2459	j = surf_usm_h%j(m)
2460	k = surf_usm_h%k(m)
2461	temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2462	ENDDO
2463	ELSE
2464	l = idsidx
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)%surface_types(m)
2470	ENDDO
2471	ENDIF
2472
2473	CASE ( 'usm_surfwintrans' )
2474	!
2475	!-- transmissivity window tiles
2476	IF ( idsint == iup_u ) THEN
2477	DO m = 1, surf_usm_h%ns
2478	i = surf_usm_h%i(m)
2479	j = surf_usm_h%j(m)
2480	k = surf_usm_h%k(m)
2481	temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
2482	ENDDO
2483	ELSE
2484	l = idsidx
2485	DO m = 1, surf_usm_v(l)%ns
2486	i = surf_usm_v(l)%i(m)
2487	j = surf_usm_v(l)%j(m)
2488	k = surf_usm_v(l)%k(m)
2489	temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
2490	ENDDO
2491	ENDIF
2492
2493	CASE ( 'usm_wshf' )
2494	!
2495	!-- array of sensible heat flux from surfaces
2496	IF ( av == 0 ) THEN
2497	IF ( idsint == iup_u ) THEN
2498	DO m = 1, surf_usm_h%ns
2499	i = surf_usm_h%i(m)
2500	j = surf_usm_h%j(m)
2501	k = surf_usm_h%k(m)
2502	temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2503	ENDDO
2504	ELSE
2505	l = idsidx
2506	DO m = 1, surf_usm_v(l)%ns
2507	i = surf_usm_v(l)%i(m)
2508	j = surf_usm_v(l)%j(m)
2509	k = surf_usm_v(l)%k(m)
2510	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2511	ENDDO
2512	ENDIF
2513	ELSE
2514	IF ( idsint == iup_u ) THEN
2515	DO m = 1, surf_usm_h%ns
2516	i = surf_usm_h%i(m)
2517	j = surf_usm_h%j(m)
2518	k = surf_usm_h%k(m)
2519	temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2520	ENDDO
2521	ELSE
2522	l = idsidx
2523	DO m = 1, surf_usm_v(l)%ns
2524	i = surf_usm_v(l)%i(m)
2525	j = surf_usm_v(l)%j(m)
2526	k = surf_usm_v(l)%k(m)
2527	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2528	ENDDO
2529	ENDIF
2530	ENDIF
2531
2532
2533	CASE ( 'usm_qsws' )
2534	!
2535	!-- array of latent heat flux from surfaces
2536	IF ( av == 0 ) THEN
2537	IF ( idsint == iup_u ) THEN
2538	DO m = 1, surf_usm_h%ns
2539	i = surf_usm_h%i(m)
2540	j = surf_usm_h%j(m)
2541	k = surf_usm_h%k(m)
2542	temp_pf(k,j,i) = surf_usm_h%qsws(m) * l_v
2543	ENDDO
2544	ELSE
2545	l = idsidx
2546	DO m = 1, surf_usm_v(l)%ns
2547	i = surf_usm_v(l)%i(m)
2548	j = surf_usm_v(l)%j(m)
2549	k = surf_usm_v(l)%k(m)
2550	temp_pf(k,j,i) = surf_usm_v(l)%qsws(m) * l_v
2551	ENDDO
2552	ENDIF
2553	ELSE
2554	IF ( idsint == iup_u ) THEN
2555	DO m = 1, surf_usm_h%ns
2556	i = surf_usm_h%i(m)
2557	j = surf_usm_h%j(m)
2558	k = surf_usm_h%k(m)
2559	temp_pf(k,j,i) = surf_usm_h%qsws_av(m)
2560	ENDDO
2561	ELSE
2562	l = idsidx
2563	DO m = 1, surf_usm_v(l)%ns
2564	i = surf_usm_v(l)%i(m)
2565	j = surf_usm_v(l)%j(m)
2566	k = surf_usm_v(l)%k(m)
2567	temp_pf(k,j,i) = surf_usm_v(l)%qsws_av(m)
2568	ENDDO
2569	ENDIF
2570	ENDIF
2571
2572	CASE ( 'usm_qsws_veg' )
2573	!
2574	!-- array of latent heat flux from vegetation surfaces
2575	IF ( av == 0 ) THEN
2576	IF ( idsint == iup_u ) THEN
2577	DO m = 1, surf_usm_h%ns
2578	i = surf_usm_h%i(m)
2579	j = surf_usm_h%j(m)
2580	k = surf_usm_h%k(m)
2581	temp_pf(k,j,i) = surf_usm_h%qsws_veg(m)
2582	ENDDO
2583	ELSE
2584	l = idsidx
2585	DO m = 1, surf_usm_v(l)%ns
2586	i = surf_usm_v(l)%i(m)
2587	j = surf_usm_v(l)%j(m)
2588	k = surf_usm_v(l)%k(m)
2589	temp_pf(k,j,i) = surf_usm_v(l)%qsws_veg(m)
2590	ENDDO
2591	ENDIF
2592	ELSE
2593	IF ( idsint == iup_u ) THEN
2594	DO m = 1, surf_usm_h%ns
2595	i = surf_usm_h%i(m)
2596	j = surf_usm_h%j(m)
2597	k = surf_usm_h%k(m)
2598	temp_pf(k,j,i) = surf_usm_h%qsws_veg_av(m)
2599	ENDDO
2600	ELSE
2601	l = idsidx
2602	DO m = 1, surf_usm_v(l)%ns
2603	i = surf_usm_v(l)%i(m)
2604	j = surf_usm_v(l)%j(m)
2605	k = surf_usm_v(l)%k(m)
2606	temp_pf(k,j,i) = surf_usm_v(l)%qsws_veg_av(m)
2607	ENDDO
2608	ENDIF
2609	ENDIF
2610
2611	CASE ( 'usm_qsws_liq' )
2612	!
2613	!-- array of latent heat flux from surfaces with liquid
2614	IF ( av == 0 ) THEN
2615	IF ( idsint == iup_u ) THEN
2616	DO m = 1, surf_usm_h%ns
2617	i = surf_usm_h%i(m)
2618	j = surf_usm_h%j(m)
2619	k = surf_usm_h%k(m)
2620	temp_pf(k,j,i) = surf_usm_h%qsws_liq(m)
2621	ENDDO
2622	ELSE
2623	l = idsidx
2624	DO m = 1, surf_usm_v(l)%ns
2625	i = surf_usm_v(l)%i(m)
2626	j = surf_usm_v(l)%j(m)
2627	k = surf_usm_v(l)%k(m)
2628	temp_pf(k,j,i) = surf_usm_v(l)%qsws_liq(m)
2629	ENDDO
2630	ENDIF
2631	ELSE
2632	IF ( idsint == iup_u ) THEN
2633	DO m = 1, surf_usm_h%ns
2634	i = surf_usm_h%i(m)
2635	j = surf_usm_h%j(m)
2636	k = surf_usm_h%k(m)
2637	temp_pf(k,j,i) = surf_usm_h%qsws_liq_av(m)
2638	ENDDO
2639	ELSE
2640	l = idsidx
2641	DO m = 1, surf_usm_v(l)%ns
2642	i = surf_usm_v(l)%i(m)
2643	j = surf_usm_v(l)%j(m)
2644	k = surf_usm_v(l)%k(m)
2645	temp_pf(k,j,i) = surf_usm_v(l)%qsws_liq_av(m)
2646	ENDDO
2647	ENDIF
2648	ENDIF
2649
2650	CASE ( 'usm_wghf' )
2651	!
2652	!-- array of heat flux from ground (land, wall, roof)
2653	IF ( av == 0 ) THEN
2654	IF ( idsint == iup_u ) THEN
2655	DO m = 1, surf_usm_h%ns
2656	i = surf_usm_h%i(m)
2657	j = surf_usm_h%j(m)
2658	k = surf_usm_h%k(m)
2659	temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2660	ENDDO
2661	ELSE
2662	l = idsidx
2663	DO m = 1, surf_usm_v(l)%ns
2664	i = surf_usm_v(l)%i(m)
2665	j = surf_usm_v(l)%j(m)
2666	k = surf_usm_v(l)%k(m)
2667	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2668	ENDDO
2669	ENDIF
2670	ELSE
2671	IF ( idsint == iup_u ) THEN
2672	DO m = 1, surf_usm_h%ns
2673	i = surf_usm_h%i(m)
2674	j = surf_usm_h%j(m)
2675	k = surf_usm_h%k(m)
2676	temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2677	ENDDO
2678	ELSE
2679	l = idsidx
2680	DO m = 1, surf_usm_v(l)%ns
2681	i = surf_usm_v(l)%i(m)
2682	j = surf_usm_v(l)%j(m)
2683	k = surf_usm_v(l)%k(m)
2684	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2685	ENDDO
2686	ENDIF
2687	ENDIF
2688
2689	CASE ( 'usm_wghf_window' )
2690	!
2691	!-- array of heat flux from window ground (land, wall, roof)
2692	IF ( av == 0 ) THEN
2693	IF ( idsint == iup_u ) THEN
2694	DO m = 1, surf_usm_h%ns
2695	i = surf_usm_h%i(m)
2696	j = surf_usm_h%j(m)
2697	k = surf_usm_h%k(m)
2698	temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
2699	ENDDO
2700	ELSE
2701	l = idsidx
2702	DO m = 1, surf_usm_v(l)%ns
2703	i = surf_usm_v(l)%i(m)
2704	j = surf_usm_v(l)%j(m)
2705	k = surf_usm_v(l)%k(m)
2706	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m)
2707	ENDDO
2708	ENDIF
2709	ELSE
2710	IF ( idsint == iup_u ) THEN
2711	DO m = 1, surf_usm_h%ns
2712	i = surf_usm_h%i(m)
2713	j = surf_usm_h%j(m)
2714	k = surf_usm_h%k(m)
2715	temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
2716	ENDDO
2717	ELSE
2718	l = idsidx
2719	DO m = 1, surf_usm_v(l)%ns
2720	i = surf_usm_v(l)%i(m)
2721	j = surf_usm_v(l)%j(m)
2722	k = surf_usm_v(l)%k(m)
2723	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
2724	ENDDO
2725	ENDIF
2726	ENDIF
2727
2728	CASE ( 'usm_wghf_green' )
2729	!
2730	!-- array of heat flux from green ground (land, wall, roof)
2731	IF ( av == 0 ) THEN
2732	IF ( idsint == iup_u ) THEN
2733	DO m = 1, surf_usm_h%ns
2734	i = surf_usm_h%i(m)
2735	j = surf_usm_h%j(m)
2736	k = surf_usm_h%k(m)
2737	temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
2738	ENDDO
2739	ELSE
2740	l = idsidx
2741	DO m = 1, surf_usm_v(l)%ns
2742	i = surf_usm_v(l)%i(m)
2743	j = surf_usm_v(l)%j(m)
2744	k = surf_usm_v(l)%k(m)
2745	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
2746	ENDDO
2747	ENDIF
2748	ELSE
2749	IF ( idsint == iup_u ) THEN
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%wghf_eb_green_av(m)
2755	ENDDO
2756	ELSE
2757	l = idsidx
2758	DO m = 1, surf_usm_v(l)%ns
2759	i = surf_usm_v(l)%i(m)
2760	j = surf_usm_v(l)%j(m)
2761	k = surf_usm_v(l)%k(m)
2762	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m)
2763	ENDDO
2764	ENDIF
2765	ENDIF
2766
2767	CASE ( 'usm_iwghf' )
2768	!
2769	!-- array of heat flux from indoor ground (land, wall, roof)
2770	IF ( av == 0 ) THEN
2771	IF ( idsint == iup_u ) THEN
2772	DO m = 1, surf_usm_h%ns
2773	i = surf_usm_h%i(m)
2774	j = surf_usm_h%j(m)
2775	k = surf_usm_h%k(m)
2776	temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
2777	ENDDO
2778	ELSE
2779	l = idsidx
2780	DO m = 1, surf_usm_v(l)%ns
2781	i = surf_usm_v(l)%i(m)
2782	j = surf_usm_v(l)%j(m)
2783	k = surf_usm_v(l)%k(m)
2784	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m)
2785	ENDDO
2786	ENDIF
2787	ELSE
2788	IF ( idsint == iup_u ) THEN
2789	DO m = 1, surf_usm_h%ns
2790	i = surf_usm_h%i(m)
2791	j = surf_usm_h%j(m)
2792	k = surf_usm_h%k(m)
2793	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
2794	ENDDO
2795	ELSE
2796	l = idsidx
2797	DO m = 1, surf_usm_v(l)%ns
2798	i = surf_usm_v(l)%i(m)
2799	j = surf_usm_v(l)%j(m)
2800	k = surf_usm_v(l)%k(m)
2801	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m)
2802	ENDDO
2803	ENDIF
2804	ENDIF
2805
2806	CASE ( 'usm_iwghf_window' )
2807	!
2808	!-- array of heat flux from indoor window ground (land, wall, roof)
2809	IF ( av == 0 ) THEN
2810	IF ( idsint == iup_u ) THEN
2811	DO m = 1, surf_usm_h%ns
2812	i = surf_usm_h%i(m)
2813	j = surf_usm_h%j(m)
2814	k = surf_usm_h%k(m)
2815	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
2816	ENDDO
2817	ELSE
2818	l = idsidx
2819	DO m = 1, surf_usm_v(l)%ns
2820	i = surf_usm_v(l)%i(m)
2821	j = surf_usm_v(l)%j(m)
2822	k = surf_usm_v(l)%k(m)
2823	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m)
2824	ENDDO
2825	ENDIF
2826	ELSE
2827	IF ( idsint == iup_u ) THEN
2828	DO m = 1, surf_usm_h%ns
2829	i = surf_usm_h%i(m)
2830	j = surf_usm_h%j(m)
2831	k = surf_usm_h%k(m)
2832	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
2833	ENDDO
2834	ELSE
2835	l = idsidx
2836	DO m = 1, surf_usm_v(l)%ns
2837	i = surf_usm_v(l)%i(m)
2838	j = surf_usm_v(l)%j(m)
2839	k = surf_usm_v(l)%k(m)
2840	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m)
2841	ENDDO
2842	ENDIF
2843	ENDIF
2844
2845	CASE ( 'usm_t_surf_wall' )
2846	!
2847	!-- surface temperature for surfaces
2848	IF ( av == 0 ) THEN
2849	IF ( idsint == iup_u ) THEN
2850	DO m = 1, surf_usm_h%ns
2851	i = surf_usm_h%i(m)
2852	j = surf_usm_h%j(m)
2853	k = surf_usm_h%k(m)
2854	temp_pf(k,j,i) = t_surf_wall_h(m)
2855	ENDDO
2856	ELSE
2857	l = idsidx
2858	DO m = 1, surf_usm_v(l)%ns
2859	i = surf_usm_v(l)%i(m)
2860	j = surf_usm_v(l)%j(m)
2861	k = surf_usm_v(l)%k(m)
2862	temp_pf(k,j,i) = t_surf_wall_v(l)%t(m)
2863	ENDDO
2864	ENDIF
2865	ELSE
2866	IF ( idsint == iup_u ) THEN
2867	DO m = 1, surf_usm_h%ns
2868	i = surf_usm_h%i(m)
2869	j = surf_usm_h%j(m)
2870	k = surf_usm_h%k(m)
2871	temp_pf(k,j,i) = surf_usm_h%t_surf_wall_av(m)
2872	ENDDO
2873	ELSE
2874	l = idsidx
2875	DO m = 1, surf_usm_v(l)%ns
2876	i = surf_usm_v(l)%i(m)
2877	j = surf_usm_v(l)%j(m)
2878	k = surf_usm_v(l)%k(m)
2879	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_wall_av(m)
2880	ENDDO
2881	ENDIF
2882	ENDIF
2883
2884	CASE ( 'usm_t_surf_window' )
2885	!
2886	!-- surface temperature for window surfaces
2887	IF ( av == 0 ) THEN
2888	IF ( idsint == iup_u ) THEN
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) = t_surf_window_h(m)
2894	ENDDO
2895	ELSE
2896	l = idsidx
2897	DO m = 1, surf_usm_v(l)%ns
2898	i = surf_usm_v(l)%i(m)
2899	j = surf_usm_v(l)%j(m)
2900	k = surf_usm_v(l)%k(m)
2901	temp_pf(k,j,i) = t_surf_window_v(l)%t(m)
2902	ENDDO
2903	ENDIF
2904
2905	ELSE
2906	IF ( idsint == iup_u ) THEN
2907	DO m = 1, surf_usm_h%ns
2908	i = surf_usm_h%i(m)
2909	j = surf_usm_h%j(m)
2910	k = surf_usm_h%k(m)
2911	temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
2912	ENDDO
2913	ELSE
2914	l = idsidx
2915	DO m = 1, surf_usm_v(l)%ns
2916	i = surf_usm_v(l)%i(m)
2917	j = surf_usm_v(l)%j(m)
2918	k = surf_usm_v(l)%k(m)
2919	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_window_av(m)
2920	ENDDO
2921
2922	ENDIF
2923
2924	ENDIF
2925
2926	CASE ( 'usm_t_surf_green' )
2927	!
2928	!-- surface temperature for green surfaces
2929	IF ( av == 0 ) THEN
2930	IF ( idsint == iup_u ) THEN
2931	DO m = 1, surf_usm_h%ns
2932	i = surf_usm_h%i(m)
2933	j = surf_usm_h%j(m)
2934	k = surf_usm_h%k(m)
2935	temp_pf(k,j,i) = t_surf_green_h(m)
2936	ENDDO
2937	ELSE
2938	l = idsidx
2939	DO m = 1, surf_usm_v(l)%ns
2940	i = surf_usm_v(l)%i(m)
2941	j = surf_usm_v(l)%j(m)
2942	k = surf_usm_v(l)%k(m)
2943	temp_pf(k,j,i) = t_surf_green_v(l)%t(m)
2944	ENDDO
2945	ENDIF
2946
2947	ELSE
2948	IF ( idsint == iup_u ) THEN
2949	DO m = 1, surf_usm_h%ns
2950	i = surf_usm_h%i(m)
2951	j = surf_usm_h%j(m)
2952	k = surf_usm_h%k(m)
2953	temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
2954	ENDDO
2955	ELSE
2956	l = idsidx
2957	DO m = 1, surf_usm_v(l)%ns
2958	i = surf_usm_v(l)%i(m)
2959	j = surf_usm_v(l)%j(m)
2960	k = surf_usm_v(l)%k(m)
2961	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_green_av(m)
2962	ENDDO
2963
2964	ENDIF
2965
2966	ENDIF
2967
2968	CASE ( 'usm_theta_10cm' )
2969	!
2970	!-- near surface temperature for whole surfaces
2971	IF ( av == 0 ) THEN
2972	IF ( idsint == iup_u ) THEN
2973	DO m = 1, surf_usm_h%ns
2974	i = surf_usm_h%i(m)
2975	j = surf_usm_h%j(m)
2976	k = surf_usm_h%k(m)
2977	temp_pf(k,j,i) = surf_usm_h%pt_10cm(m)
2978	ENDDO
2979	ELSE
2980	l = idsidx
2981	DO m = 1, surf_usm_v(l)%ns
2982	i = surf_usm_v(l)%i(m)
2983	j = surf_usm_v(l)%j(m)
2984	k = surf_usm_v(l)%k(m)
2985	temp_pf(k,j,i) = surf_usm_v(l)%pt_10cm(m)
2986	ENDDO
2987	ENDIF
2988
2989
2990	ELSE
2991	IF ( idsint == iup_u ) THEN
2992	DO m = 1, surf_usm_h%ns
2993	i = surf_usm_h%i(m)
2994	j = surf_usm_h%j(m)
2995	k = surf_usm_h%k(m)
2996	temp_pf(k,j,i) = surf_usm_h%pt_10cm_av(m)
2997	ENDDO
2998	ELSE
2999	l = idsidx
3000	DO m = 1, surf_usm_v(l)%ns
3001	i = surf_usm_v(l)%i(m)
3002	j = surf_usm_v(l)%j(m)
3003	k = surf_usm_v(l)%k(m)
3004	temp_pf(k,j,i) = surf_usm_v(l)%pt_10cm_av(m)
3005	ENDDO
3006
3007	ENDIF
3008	ENDIF
3009
3010	CASE ( 'usm_t_wall' )
3011	!
3012	!-- wall temperature for iwl layer of walls and land
3013	IF ( av == 0 ) THEN
3014	IF ( idsint == iup_u ) THEN
3015	DO m = 1, surf_usm_h%ns
3016	i = surf_usm_h%i(m)
3017	j = surf_usm_h%j(m)
3018	k = surf_usm_h%k(m)
3019	temp_pf(k,j,i) = t_wall_h(iwl,m)
3020	ENDDO
3021	ELSE
3022	l = idsidx
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_wall_v(l)%t(iwl,m)
3028	ENDDO
3029	ENDIF
3030	ELSE
3031	IF ( idsint == iup_u ) THEN
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_wall_av(iwl,m)
3037	ENDDO
3038	ELSE
3039	l = idsidx
3040	DO m = 1, surf_usm_v(l)%ns
3041	i = surf_usm_v(l)%i(m)
3042	j = surf_usm_v(l)%j(m)
3043	k = surf_usm_v(l)%k(m)
3044	temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
3045	ENDDO
3046	ENDIF
3047	ENDIF
3048
3049	CASE ( 'usm_t_window' )
3050	!
3051	!-- window temperature for iwl layer of walls and land
3052	IF ( av == 0 ) THEN
3053	IF ( idsint == iup_u ) 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_window_h(iwl,m)
3059	ENDDO
3060	ELSE
3061	l = idsidx
3062	DO m = 1, surf_usm_v(l)%ns
3063	i = surf_usm_v(l)%i(m)
3064	j = surf_usm_v(l)%j(m)
3065	k = surf_usm_v(l)%k(m)
3066	temp_pf(k,j,i) = t_window_v(l)%t(iwl,m)
3067	ENDDO
3068	ENDIF
3069	ELSE
3070	IF ( idsint == iup_u ) THEN
3071	DO m = 1, surf_usm_h%ns
3072	i = surf_usm_h%i(m)
3073	j = surf_usm_h%j(m)
3074	k = surf_usm_h%k(m)
3075	temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
3076	ENDDO
3077	ELSE
3078	l = idsidx
3079	DO m = 1, surf_usm_v(l)%ns
3080	i = surf_usm_v(l)%i(m)
3081	j = surf_usm_v(l)%j(m)
3082	k = surf_usm_v(l)%k(m)
3083	temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m)
3084	ENDDO
3085	ENDIF
3086	ENDIF
3087
3088	CASE ( 'usm_t_green' )
3089	!
3090	!-- green temperature for iwl layer of walls and land
3091	IF ( av == 0 ) THEN
3092	IF ( idsint == iup_u ) THEN
3093	DO m = 1, surf_usm_h%ns
3094	i = surf_usm_h%i(m)
3095	j = surf_usm_h%j(m)
3096	k = surf_usm_h%k(m)
3097	temp_pf(k,j,i) = t_green_h(iwl,m)
3098	ENDDO
3099	ELSE
3100	l = idsidx
3101	DO m = 1, surf_usm_v(l)%ns
3102	i = surf_usm_v(l)%i(m)
3103	j = surf_usm_v(l)%j(m)
3104	k = surf_usm_v(l)%k(m)
3105	temp_pf(k,j,i) = t_green_v(l)%t(iwl,m)
3106	ENDDO
3107	ENDIF
3108	ELSE
3109	IF ( idsint == iup_u ) THEN
3110	DO m = 1, surf_usm_h%ns
3111	i = surf_usm_h%i(m)
3112	j = surf_usm_h%j(m)
3113	k = surf_usm_h%k(m)
3114	temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
3115	ENDDO
3116	ELSE
3117	l = idsidx
3118	DO m = 1, surf_usm_v(l)%ns
3119	i = surf_usm_v(l)%i(m)
3120	j = surf_usm_v(l)%j(m)
3121	k = surf_usm_v(l)%k(m)
3122	temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m)
3123	ENDDO
3124	ENDIF
3125	ENDIF
3126
3127	CASE ( 'usm_swc' )
3128	!
3129	!-- soil water content for iwl layer of walls and land
3130	IF ( av == 0 ) THEN
3131	IF ( idsint == iup_u ) THEN
3132	DO m = 1, surf_usm_h%ns
3133	i = surf_usm_h%i(m)
3134	j = surf_usm_h%j(m)
3135	k = surf_usm_h%k(m)
3136	temp_pf(k,j,i) = swc_h(iwl,m)
3137	ENDDO
3138	ELSE
3139
3140	ENDIF
3141	ELSE
3142	IF ( idsint == iup_u ) THEN
3143	DO m = 1, surf_usm_h%ns
3144	i = surf_usm_h%i(m)
3145	j = surf_usm_h%j(m)
3146	k = surf_usm_h%k(m)
3147	temp_pf(k,j,i) = surf_usm_h%swc_av(iwl,m)
3148	ENDDO
3149	ELSE
3150	l = idsidx
3151	DO m = 1, surf_usm_v(l)%ns
3152	i = surf_usm_v(l)%i(m)
3153	j = surf_usm_v(l)%j(m)
3154	k = surf_usm_v(l)%k(m)
3155	temp_pf(k,j,i) = surf_usm_v(l)%swc_av(iwl,m)
3156	ENDDO
3157	ENDIF
3158	ENDIF
3159
3160
3161	CASE DEFAULT
3162	found = .FALSE.
3163	RETURN
3164	END SELECT
3165
3166	!
3167	!-- Rearrange dimensions for NetCDF output
3168	!-- FIXME: this may generate FPE overflow upon conversion from DP to SP
3169	DO j = nys, nyn
3170	DO i = nxl, nxr
3171	DO k = nzb_do, nzt_do
3172	local_pf(i,j,k) = temp_pf(k,j,i)
3173	ENDDO
3174	ENDDO
3175	ENDDO
3176
3177	END SUBROUTINE usm_data_output_3d
3178
3179
3180	!------------------------------------------------------------------------------!
3181	!
3182	! Description:
3183	! ------------
3184	!> Soubroutine defines appropriate grid for netcdf variables.
3185	!> It is called out from subroutine netcdf.
3186	!------------------------------------------------------------------------------!
3187	SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
3188
3189	IMPLICIT NONE
3190
3191	CHARACTER (len=*), INTENT(IN) :: variable !<
3192	LOGICAL, INTENT(OUT) :: found !<
3193	CHARACTER (len=*), INTENT(OUT) :: grid_x !<
3194	CHARACTER (len=*), INTENT(OUT) :: grid_y !<
3195	CHARACTER (len=*), INTENT(OUT) :: grid_z !<
3196
3197	CHARACTER (len=varnamelength) :: var
3198
3199	var = TRIM(variable)
3200	IF ( var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
3201	var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. &
3202	var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. &
3203	var(1:9) == 'usm_qsws_' .OR. var(1:13) == 'usm_qsws_veg_' .OR. &
3204	var(1:13) == 'usm_qsws_liq_' .OR. &
3205	var(1:15) == 'usm_t_surf_wall' .OR. var(1:10) == 'usm_t_wall' .OR. &
3206	var(1:17) == 'usm_t_surf_window' .OR. var(1:12) == 'usm_t_window' .OR. &
3207	var(1:16) == 'usm_t_surf_green' .OR. var(1:11) == 'usm_t_green' .OR. &
3208	var(1:15) == 'usm_theta_10cm' .OR. &
3209	var(1:9) == 'usm_surfz' .OR. var(1:11) == 'usm_surfcat' .OR. &
3210	var(1:16) == 'usm_surfwintrans' .OR. var(1:7) == 'usm_swc' ) THEN
3211
3212	found = .TRUE.
3213	grid_x = 'x'
3214	grid_y = 'y'
3215	grid_z = 'zu'
3216	ELSE
3217	found = .FALSE.
3218	grid_x = 'none'
3219	grid_y = 'none'
3220	grid_z = 'none'
3221	ENDIF
3222
3223	END SUBROUTINE usm_define_netcdf_grid
3224
3225
3226	!------------------------------------------------------------------------------!
3227	! Description:
3228	! ------------
3229	!> Initialization of the wall surface model
3230	!------------------------------------------------------------------------------!
3231	SUBROUTINE usm_init_material_model
3232
3233	IMPLICIT NONE
3234
3235	INTEGER(iwp) :: k, l, m !< running indices
3236
3237	IF ( debug_output ) CALL debug_message( 'usm_init_material_model', 'start' )
3238
3239	!
3240	!-- Calculate wall grid spacings.
3241	!-- Temperature is defined at the center of the wall layers,
3242	!-- whereas gradients/fluxes are defined at the edges (_stag)
3243	!-- apply for all particular surface grids. First for horizontal surfaces
3244	DO m = 1, surf_usm_h%ns
3245
3246	surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
3247	DO k = nzb_wall+1, nzt_wall
3248	surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) - &
3249	surf_usm_h%zw(k-1,m)
3250	ENDDO
3251	surf_usm_h%dz_window(nzb_wall,m) = surf_usm_h%zw_window(nzb_wall,m)
3252	DO k = nzb_wall+1, nzt_wall
3253	surf_usm_h%dz_window(k,m) = surf_usm_h%zw_window(k,m) - &
3254	surf_usm_h%zw_window(k-1,m)
3255	ENDDO
3256
3257	surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
3258
3259	DO k = nzb_wall, nzt_wall-1
3260	surf_usm_h%dz_wall_stag(k,m) = 0.5 * ( &
3261	surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
3262	ENDDO
3263	surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
3264
3265	surf_usm_h%dz_window(nzt_wall+1,m) = surf_usm_h%dz_window(nzt_wall,m)
3266
3267	DO k = nzb_wall, nzt_wall-1
3268	surf_usm_h%dz_window_stag(k,m) = 0.5 * ( &
3269	surf_usm_h%dz_window(k+1,m) + surf_usm_h%dz_window(k,m) )
3270	ENDDO
3271	surf_usm_h%dz_window_stag(nzt_wall,m) = surf_usm_h%dz_window(nzt_wall,m)
3272
3273	IF (surf_usm_h%green_type_roof(m) == 2.0_wp ) THEN
3274	!
3275	!-- extensive green roof
3276	!-- set ratio of substrate layer thickness, soil-type and LAI
3277	soil_type = 3
3278	surf_usm_h%lai(m) = 2.0_wp
3279
3280	surf_usm_h%zw_green(nzb_wall,m) = 0.05_wp
3281	surf_usm_h%zw_green(nzb_wall+1,m) = 0.10_wp
3282	surf_usm_h%zw_green(nzb_wall+2,m) = 0.15_wp
3283	surf_usm_h%zw_green(nzb_wall+3,m) = 0.20_wp
3284	ELSE
3285	!
3286	!-- intensiv green roof
3287	!-- set ratio of substrate layer thickness, soil-type and LAI
3288	soil_type = 6
3289	surf_usm_h%lai(m) = 4.0_wp
3290
3291	surf_usm_h%zw_green(nzb_wall,m) = 0.05_wp
3292	surf_usm_h%zw_green(nzb_wall+1,m) = 0.10_wp
3293	surf_usm_h%zw_green(nzb_wall+2,m) = 0.40_wp
3294	surf_usm_h%zw_green(nzb_wall+3,m) = 0.80_wp
3295	ENDIF
3296
3297	surf_usm_h%dz_green(nzb_wall,m) = surf_usm_h%zw_green(nzb_wall,m)
3298	DO k = nzb_wall+1, nzt_wall
3299	surf_usm_h%dz_green(k,m) = surf_usm_h%zw_green(k,m) - &
3300	surf_usm_h%zw_green(k-1,m)
3301	ENDDO
3302	surf_usm_h%dz_green(nzt_wall+1,m) = surf_usm_h%dz_green(nzt_wall,m)
3303
3304	DO k = nzb_wall, nzt_wall-1
3305	surf_usm_h%dz_green_stag(k,m) = 0.5 * ( &
3306	surf_usm_h%dz_green(k+1,m) + surf_usm_h%dz_green(k,m) )
3307	ENDDO
3308	surf_usm_h%dz_green_stag(nzt_wall,m) = surf_usm_h%dz_green(nzt_wall,m)
3309
3310	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
3311	alpha_vangenuchten = soil_pars(0,soil_type)
3312	ENDIF
3313
3314	IF ( l_vangenuchten == 9999999.9_wp ) THEN
3315	l_vangenuchten = soil_pars(1,soil_type)
3316	ENDIF
3317
3318	IF ( n_vangenuchten == 9999999.9_wp ) THEN
3319	n_vangenuchten = soil_pars(2,soil_type)
3320	ENDIF
3321
3322	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
3323	hydraulic_conductivity = soil_pars(3,soil_type)
3324	ENDIF
3325
3326	IF ( saturation_moisture == 9999999.9_wp ) THEN
3327	saturation_moisture = m_soil_pars(0,soil_type)
3328	ENDIF
3329
3330	IF ( field_capacity == 9999999.9_wp ) THEN
3331	field_capacity = m_soil_pars(1,soil_type)
3332	ENDIF
3333
3334	IF ( wilting_point == 9999999.9_wp ) THEN
3335	wilting_point = m_soil_pars(2,soil_type)
3336	ENDIF
3337
3338	IF ( residual_moisture == 9999999.9_wp ) THEN
3339	residual_moisture = m_soil_pars(3,soil_type)
3340	ENDIF
3341
3342	DO k = nzb_wall, nzt_wall+1
3343	swc_h(k,m) = field_capacity
3344	rootfr_h(k,m) = 0.5_wp
3345	surf_usm_h%alpha_vg_green(m) = alpha_vangenuchten
3346	surf_usm_h%l_vg_green(m) = l_vangenuchten
3347	surf_usm_h%n_vg_green(m) = n_vangenuchten
3348	surf_usm_h%gamma_w_green_sat(k,m) = hydraulic_conductivity
3349	swc_sat_h(k,m) = saturation_moisture
3350	fc_h(k,m) = field_capacity
3351	wilt_h(k,m) = wilting_point
3352	swc_res_h(k,m) = residual_moisture
3353	ENDDO
3354
3355	ENDDO
3356
3357	surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall
3358	surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag
3359	surf_usm_h%ddz_window = 1.0_wp / surf_usm_h%dz_window
3360	surf_usm_h%ddz_window_stag = 1.0_wp / surf_usm_h%dz_window_stag
3361	surf_usm_h%ddz_green = 1.0_wp / surf_usm_h%dz_green
3362	surf_usm_h%ddz_green_stag = 1.0_wp / surf_usm_h%dz_green_stag
3363	!
3364	!-- For vertical surfaces
3365	DO l = 0, 3
3366	DO m = 1, surf_usm_v(l)%ns
3367	surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
3368	DO k = nzb_wall+1, nzt_wall
3369	surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - &
3370	surf_usm_v(l)%zw(k-1,m)
3371	ENDDO
3372	surf_usm_v(l)%dz_window(nzb_wall,m) = surf_usm_v(l)%zw_window(nzb_wall,m)
3373	DO k = nzb_wall+1, nzt_wall
3374	surf_usm_v(l)%dz_window(k,m) = surf_usm_v(l)%zw_window(k,m) - &
3375	surf_usm_v(l)%zw_window(k-1,m)
3376	ENDDO
3377	surf_usm_v(l)%dz_green(nzb_wall,m) = surf_usm_v(l)%zw_green(nzb_wall,m)
3378	DO k = nzb_wall+1, nzt_wall
3379	surf_usm_v(l)%dz_green(k,m) = surf_usm_v(l)%zw_green(k,m) - &
3380	surf_usm_v(l)%zw_green(k-1,m)
3381	ENDDO
3382
3383	surf_usm_v(l)%dz_wall(nzt_wall+1,m) = &
3384	surf_usm_v(l)%dz_wall(nzt_wall,m)
3385
3386	DO k = nzb_wall, nzt_wall-1
3387	surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * ( &
3388	surf_usm_v(l)%dz_wall(k+1,m) + &
3389	surf_usm_v(l)%dz_wall(k,m) )
3390	ENDDO
3391	surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = &
3392	surf_usm_v(l)%dz_wall(nzt_wall,m)
3393	surf_usm_v(l)%dz_window(nzt_wall+1,m) = &
3394	surf_usm_v(l)%dz_window(nzt_wall,m)
3395
3396	DO k = nzb_wall, nzt_wall-1
3397	surf_usm_v(l)%dz_window_stag(k,m) = 0.5 * ( &
3398	surf_usm_v(l)%dz_window(k+1,m) + &
3399	surf_usm_v(l)%dz_window(k,m) )
3400	ENDDO
3401	surf_usm_v(l)%dz_window_stag(nzt_wall,m) = &
3402	surf_usm_v(l)%dz_window(nzt_wall,m)
3403	surf_usm_v(l)%dz_green(nzt_wall+1,m) = &
3404	surf_usm_v(l)%dz_green(nzt_wall,m)
3405
3406	DO k = nzb_wall, nzt_wall-1
3407	surf_usm_v(l)%dz_green_stag(k,m) = 0.5 * ( &
3408	surf_usm_v(l)%dz_green(k+1,m) + &
3409	surf_usm_v(l)%dz_green(k,m) )
3410	ENDDO
3411	surf_usm_v(l)%dz_green_stag(nzt_wall,m) = &
3412	surf_usm_v(l)%dz_green(nzt_wall,m)
3413	ENDDO
3414	surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall
3415	surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag
3416	surf_usm_v(l)%ddz_window = 1.0_wp / surf_usm_v(l)%dz_window
3417	surf_usm_v(l)%ddz_window_stag = 1.0_wp / surf_usm_v(l)%dz_window_stag
3418	surf_usm_v(l)%ddz_green = 1.0_wp / surf_usm_v(l)%dz_green
3419	surf_usm_v(l)%ddz_green_stag = 1.0_wp / surf_usm_v(l)%dz_green_stag
3420	ENDDO
3421
3422
3423	IF ( debug_output ) CALL debug_message( 'usm_init_material_model', 'end' )
3424
3425	END SUBROUTINE usm_init_material_model
3426
3427
3428	!------------------------------------------------------------------------------!
3429	! Description:
3430	! ------------
3431	!> Initialization of the urban surface model
3432	!------------------------------------------------------------------------------!
3433	SUBROUTINE usm_init
3434
3435	USE arrays_3d, &
3436	ONLY: zw
3437
3438	USE netcdf_data_input_mod, &
3439	ONLY: building_pars_f, building_surface_pars_f, building_type_f, &
3440	terrain_height_f
3441
3442	IMPLICIT NONE
3443
3444	INTEGER(iwp) :: i !< loop index x-dirction
3445	INTEGER(iwp) :: ind_alb_green !< index in input list for green albedo
3446	INTEGER(iwp) :: ind_alb_wall !< index in input list for wall albedo
3447	INTEGER(iwp) :: ind_alb_win !< index in input list for window albedo
3448	INTEGER(iwp) :: ind_emis_wall !< index in input list for wall emissivity
3449	INTEGER(iwp) :: ind_emis_green !< index in input list for green emissivity
3450	INTEGER(iwp) :: ind_emis_win !< index in input list for window emissivity
3451	INTEGER(iwp) :: ind_green_frac_w !< index in input list for green fraction on wall
3452	INTEGER(iwp) :: ind_green_frac_r !< index in input list for green fraction on roof
3453	INTEGER(iwp) :: ind_hc1 !< index in input list for heat capacity at first wall layer
3454	INTEGER(iwp) :: ind_hc1_win !< index in input list for heat capacity at first window layer
3455	INTEGER(iwp) :: ind_hc2 !< index in input list for heat capacity at second wall layer
3456	INTEGER(iwp) :: ind_hc2_win !< index in input list for heat capacity at second window layer
3457	INTEGER(iwp) :: ind_hc3 !< index in input list for heat capacity at third wall layer
3458	INTEGER(iwp) :: ind_hc3_win !< index in input list for heat capacity at third window layer
3459	INTEGER(iwp) :: ind_lai_r !< index in input list for LAI on roof
3460	INTEGER(iwp) :: ind_lai_w !< index in input list for LAI on wall
3461	INTEGER(iwp) :: ind_tc1 !< index in input list for thermal conductivity at first wall layer
3462	INTEGER(iwp) :: ind_tc1_win !< index in input list for thermal conductivity at first window layer
3463	INTEGER(iwp) :: ind_tc2 !< index in input list for thermal conductivity at second wall layer
3464	INTEGER(iwp) :: ind_tc2_win !< index in input list for thermal conductivity at second window layer
3465	INTEGER(iwp) :: ind_tc3 !< index in input list for thermal conductivity at third wall layer
3466	INTEGER(iwp) :: ind_tc3_win !< index in input list for thermal conductivity at third window layer
3467	INTEGER(iwp) :: ind_thick_1 !< index in input list for thickness of first wall layer
3468	INTEGER(iwp) :: ind_thick_1_win !< index in input list for thickness of first window layer
3469	INTEGER(iwp) :: ind_thick_2 !< index in input list for thickness of second wall layer
3470	INTEGER(iwp) :: ind_thick_2_win !< index in input list for thickness of second window layer
3471	INTEGER(iwp) :: ind_thick_3 !< index in input list for thickness of third wall layer
3472	INTEGER(iwp) :: ind_thick_3_win !< index in input list for thickness of third window layer
3473	INTEGER(iwp) :: ind_thick_4 !< index in input list for thickness of fourth wall layer
3474	INTEGER(iwp) :: ind_thick_4_win !< index in input list for thickness of fourth window layer
3475	INTEGER(iwp) :: ind_trans !< index in input list for window transmissivity
3476	INTEGER(iwp) :: ind_wall_frac !< index in input list for wall fraction
3477	INTEGER(iwp) :: ind_win_frac !< index in input list for window fraction
3478	INTEGER(iwp) :: ind_z0 !< index in input list for z0
3479	INTEGER(iwp) :: ind_z0qh !< index in input list for z0h / z0q
3480	INTEGER(iwp) :: is !< loop index input surface element
3481	INTEGER(iwp) :: j !< loop index y-dirction
3482	INTEGER(iwp) :: k !< loop index z-dirction
3483	INTEGER(iwp) :: l !< loop index surface orientation
3484	INTEGER(iwp) :: m !< loop index surface element
3485	INTEGER(iwp) :: st !< dummy
3486
3487	LOGICAL :: relative_fractions_corrected !< flag indicating if relative surface fractions require normalization
3488
3489	REAL(wp) :: c, tin, twin
3490	REAL(wp) :: ground_floor_level_l !< local height of ground floor level
3491	REAL(wp) :: sum_frac !< sum of the relative material fractions at a surface element
3492	REAL(wp) :: z_agl !< height of the surface element above terrain
3493
3494	IF ( debug_output ) CALL debug_message( 'usm_init', 'start' )
3495
3496	CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
3497	!
3498	!-- surface forcing have to be disabled for LSF
3499	!-- in case of enabled urban surface module
3500	IF ( large_scale_forcing ) THEN
3501	lsf_surf = .FALSE.
3502	ENDIF
3503	!
3504	!-- Flag surface elements belonging to the ground floor level. Therefore,
3505	!-- use terrain height array from file, if available. This flag is later used
3506	!-- to control initialization of surface attributes.
3507	!-- Todo: for the moment disable initialization of building roofs with
3508	!-- ground-floor-level properties.
3509	surf_usm_h%ground_level = .FALSE.
3510
3511	DO l = 0, 3
3512	surf_usm_v(l)%ground_level = .FALSE.
3513	DO m = 1, surf_usm_v(l)%ns
3514	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
3515	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
3516	k = surf_usm_v(l)%k(m)
3517	!
3518	!-- Determine local ground level. Level 1 - default value,
3519	!-- level 2 - initialization according to building type,
3520	!-- level 3 - initialization from value read from file.
3521	ground_floor_level_l = ground_floor_level
3522
3523	IF ( building_type_f%from_file ) THEN
3524	ground_floor_level_l = &
3525	building_pars(ind_gflh,building_type_f%var(j,i))
3526	ENDIF
3527
3528	IF ( building_pars_f%from_file ) THEN
3529	IF ( building_pars_f%pars_xy(ind_gflh,j,i) /= &
3530	building_pars_f%fill ) &
3531	ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)
3532	ENDIF
3533	!
3534	!-- Determine height of surface element above ground level. Please
3535	!-- note, height of surface element is determined with respect to
3536	!-- its height above ground of the reference grid point in atmosphere,
3537	!-- Therefore, substract the offset values when assessing the terrain
3538	!-- height.
3539	IF ( terrain_height_f%from_file ) THEN
3540	z_agl = zw(k) - terrain_height_f%var(j-surf_usm_v(l)%joff, &
3541	i-surf_usm_v(l)%ioff)
3542	ELSE
3543	z_agl = zw(k)
3544	ENDIF
3545	!
3546	!-- Set flag for ground level
3547	IF ( z_agl <= ground_floor_level_l ) &
3548	surf_usm_v(l)%ground_level(m) = .TRUE.
3549
3550	ENDDO
3551	ENDDO
3552	!
3553	!-- Initialization of resistances.
3554	DO m = 1, surf_usm_h%ns
3555	surf_usm_h%r_a(m) = 50.0_wp
3556	surf_usm_h%r_a_green(m) = 50.0_wp
3557	surf_usm_h%r_a_window(m) = 50.0_wp
3558	ENDDO
3559	DO l = 0, 3
3560	DO m = 1, surf_usm_v(l)%ns
3561	surf_usm_v(l)%r_a(m) = 50.0_wp
3562	surf_usm_v(l)%r_a_green(m) = 50.0_wp
3563	surf_usm_v(l)%r_a_window(m) = 50.0_wp
3564	ENDDO
3565	ENDDO
3566
3567	!
3568	!-- Map values onto horizontal elemements
3569	DO m = 1, surf_usm_h%ns
3570	surf_usm_h%r_canopy(m) = 200.0_wp !< canopy_resistance
3571	surf_usm_h%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
3572	surf_usm_h%g_d(m) = 0.0_wp !< canopy_resistance_coefficient
3573	ENDDO
3574	!
3575	!-- Map values onto vertical elements, even though this does not make
3576	!-- much sense.
3577	DO l = 0, 3
3578	DO m = 1, surf_usm_v(l)%ns
3579	surf_usm_v(l)%r_canopy(m) = 200.0_wp !< canopy_resistance
3580	surf_usm_v(l)%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
3581	surf_usm_v(l)%g_d(m) = 0.0_wp !< canopy_resistance_coefficient
3582	ENDDO
3583	ENDDO
3584
3585	!
3586	!-- Initialize urban-type surface attribute. According to initialization in
3587	!-- land-surface model, follow a 3-level approach.
3588	!-- Level 1 - initialization via default attributes
3589	DO m = 1, surf_usm_h%ns
3590	!
3591	!-- Now, all horizontal surfaces are roof surfaces (?)
3592	surf_usm_h%isroof_surf(m) = .TRUE.
3593	surf_usm_h%surface_types(m) = roof_category !< default category for root surface
3594	!
3595	!-- In order to distinguish between ground floor level and
3596	!-- above-ground-floor level surfaces, set input indices.
3597
3598	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3599	surf_usm_h%ground_level(m) )
3600	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3601	surf_usm_h%ground_level(m) )
3602	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3603	surf_usm_h%ground_level(m) )
3604	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3605	surf_usm_h%ground_level(m) )
3606	!
3607	!-- Store building type and its name on each surface element
3608	surf_usm_h%building_type(m) = building_type
3609	surf_usm_h%building_type_name(m) = building_type_name(building_type)
3610	!
3611	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3612	surf_usm_h%frac(m,ind_veg_wall) = building_pars(ind_wall_frac_r,building_type)
3613	surf_usm_h%frac(m,ind_pav_green) = building_pars(ind_green_frac_r,building_type)
3614	surf_usm_h%frac(m,ind_wat_win) = building_pars(ind_win_frac_r,building_type)
3615	surf_usm_h%lai(m) = building_pars(ind_lai_r,building_type)
3616
3617	surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1_wall_r,building_type)
3618	surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1_wall_r,building_type)
3619	surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2_wall_r,building_type)
3620	surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3_wall_r,building_type)
3621	surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1_wall_r,building_type)
3622	surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1_wall_r,building_type)
3623	surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2_wall_r,building_type)
3624	surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3_wall_r,building_type)
3625	surf_usm_h%rho_c_green(nzb_wall,m) = rho_c_soil !building_pars(ind_hc1_wall_r,building_type)
3626	surf_usm_h%rho_c_green(nzb_wall+1,m) = rho_c_soil !building_pars(ind_hc1_wall_r,building_type)
3627	surf_usm_h%rho_c_green(nzb_wall+2,m) = rho_c_soil !building_pars(ind_hc2_wall_r,building_type)
3628	surf_usm_h%rho_c_green(nzb_wall+3,m) = rho_c_soil !building_pars(ind_hc3_wall_r,building_type)
3629	surf_usm_h%lambda_h_green(nzb_wall,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,building_type)
3630	surf_usm_h%lambda_h_green(nzb_wall+1,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,building_type)
3631	surf_usm_h%lambda_h_green(nzb_wall+2,m) = lambda_h_green_sm !building_pars(ind_tc2_wall_r,building_type)
3632	surf_usm_h%lambda_h_green(nzb_wall+3,m) = lambda_h_green_sm !building_pars(ind_tc3_wall_r,building_type)
3633	surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1_win_r,building_type)
3634	surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1_win_r,building_type)
3635	surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2_win_r,building_type)
3636	surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3_win_r,building_type)
3637	surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1_win_r,building_type)
3638	surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1_win_r,building_type)
3639	surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2_win_r,building_type)
3640	surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3_win_r,building_type)
3641
3642	surf_usm_h%target_temp_summer(m) = building_pars(ind_indoor_target_temp_summer,building_type)
3643	surf_usm_h%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,building_type)
3644	!
3645	!-- emissivity of wall-, green- and window fraction
3646	surf_usm_h%emissivity(m,ind_veg_wall) = building_pars(ind_emis_wall_r,building_type)
3647	surf_usm_h%emissivity(m,ind_pav_green) = building_pars(ind_emis_green_r,building_type)
3648	surf_usm_h%emissivity(m,ind_wat_win) = building_pars(ind_emis_win_r,building_type)
3649
3650	surf_usm_h%transmissivity(m) = building_pars(ind_trans_r,building_type)
3651
3652	surf_usm_h%z0(m) = building_pars(ind_z0,building_type)
3653	surf_usm_h%z0h(m) = building_pars(ind_z0qh,building_type)
3654	surf_usm_h%z0q(m) = building_pars(ind_z0qh,building_type)
3655	!
3656	!-- albedo type for wall fraction, green fraction, window fraction
3657	surf_usm_h%albedo_type(m,ind_veg_wall) = INT( building_pars(ind_alb_wall_r,building_type) )
3658	surf_usm_h%albedo_type(m,ind_pav_green) = INT( building_pars(ind_alb_green_r,building_type) )
3659	surf_usm_h%albedo_type(m,ind_wat_win) = INT( building_pars(ind_alb_win_r,building_type) )
3660
3661	surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1_wall_r,building_type)
3662	surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2_wall_r,building_type)
3663	surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3_wall_r,building_type)
3664	surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4_wall_r,building_type)
3665
3666	surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1_wall_r,building_type)
3667	surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2_wall_r,building_type)
3668	surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3_wall_r,building_type)
3669	surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4_wall_r,building_type)
3670
3671	surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1_win_r,building_type)
3672	surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2_win_r,building_type)
3673	surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3_win_r,building_type)
3674	surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4_win_r,building_type)
3675
3676	surf_usm_h%green_type_roof(m) = building_pars(ind_green_type_roof,building_type)
3677
3678	ENDDO
3679
3680	DO l = 0, 3
3681	DO m = 1, surf_usm_v(l)%ns
3682
3683	surf_usm_v(l)%surface_types(m) = wall_category !< default category for root surface
3684	!
3685	!-- In order to distinguish between ground floor level and
3686	!-- above-ground-floor level surfaces, set input indices.
3687	ind_alb_green = MERGE( ind_alb_green_gfl, ind_alb_green_agfl, &
3688	surf_usm_v(l)%ground_level(m) )
3689	ind_alb_wall = MERGE( ind_alb_wall_gfl, ind_alb_wall_agfl, &
3690	surf_usm_v(l)%ground_level(m) )
3691	ind_alb_win = MERGE( ind_alb_win_gfl, ind_alb_win_agfl, &
3692	surf_usm_v(l)%ground_level(m) )
3693	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3694	surf_usm_v(l)%ground_level(m) )
3695	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3696	surf_usm_v(l)%ground_level(m) )
3697	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3698	surf_usm_v(l)%ground_level(m) )
3699	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3700	surf_usm_v(l)%ground_level(m) )
3701	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3702	surf_usm_v(l)%ground_level(m) )
3703	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
3704	surf_usm_v(l)%ground_level(m) )
3705	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
3706	surf_usm_v(l)%ground_level(m) )
3707	ind_hc1_win = MERGE( ind_hc1_win_gfl, ind_hc1_win_agfl, &
3708	surf_usm_v(l)%ground_level(m) )
3709	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
3710	surf_usm_v(l)%ground_level(m) )
3711	ind_hc2_win = MERGE( ind_hc2_win_gfl, ind_hc2_win_agfl, &
3712	surf_usm_v(l)%ground_level(m) )
3713	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
3714	surf_usm_v(l)%ground_level(m) )
3715	ind_hc3_win = MERGE( ind_hc3_win_gfl, ind_hc3_win_agfl, &
3716	surf_usm_v(l)%ground_level(m) )
3717	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
3718	surf_usm_v(l)%ground_level(m) )
3719	ind_tc1_win = MERGE( ind_tc1_win_gfl, ind_tc1_win_agfl, &
3720	surf_usm_v(l)%ground_level(m) )
3721	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
3722	surf_usm_v(l)%ground_level(m) )
3723	ind_tc2_win = MERGE( ind_tc2_win_gfl, ind_tc2_win_agfl, &
3724	surf_usm_v(l)%ground_level(m) )
3725	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
3726	surf_usm_v(l)%ground_level(m) )
3727	ind_tc3_win = MERGE( ind_tc3_win_gfl, ind_tc3_win_agfl, &
3728	surf_usm_v(l)%ground_level(m) )
3729	ind_thick_1 = MERGE( ind_thick_1_gfl, ind_thick_1_agfl, &
3730	surf_usm_v(l)%ground_level(m) )
3731	ind_thick_1_win = MERGE( ind_thick_1_win_gfl, ind_thick_1_win_agfl, &
3732	surf_usm_v(l)%ground_level(m) )
3733	ind_thick_2 = MERGE( ind_thick_2_gfl, ind_thick_2_agfl, &
3734	surf_usm_v(l)%ground_level(m) )
3735	ind_thick_2_win = MERGE( ind_thick_2_win_gfl, ind_thick_2_win_agfl, &
3736	surf_usm_v(l)%ground_level(m) )
3737	ind_thick_3 = MERGE( ind_thick_3_gfl, ind_thick_3_agfl, &
3738	surf_usm_v(l)%ground_level(m) )
3739	ind_thick_3_win = MERGE( ind_thick_3_win_gfl, ind_thick_3_win_agfl, &
3740	surf_usm_v(l)%ground_level(m) )
3741	ind_thick_4 = MERGE( ind_thick_4_gfl, ind_thick_4_agfl, &
3742	surf_usm_v(l)%ground_level(m) )
3743	ind_thick_4_win = MERGE( ind_thick_4_win_gfl, ind_thick_4_win_agfl, &
3744	surf_usm_v(l)%ground_level(m) )
3745	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
3746	surf_usm_v(l)%ground_level(m) )
3747	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
3748	surf_usm_v(l)%ground_level(m) )
3749	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
3750	surf_usm_v(l)%ground_level(m) )
3751	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
3752	surf_usm_v(l)%ground_level(m) )
3753	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3754	surf_usm_v(l)%ground_level(m) )
3755	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3756	surf_usm_v(l)%ground_level(m) )
3757	!
3758	!-- Store building type and its name on each surface element
3759	surf_usm_v(l)%building_type(m) = building_type
3760	surf_usm_v(l)%building_type_name(m) = building_type_name(building_type)
3761	!
3762	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3763	surf_usm_v(l)%frac(m,ind_veg_wall) = building_pars(ind_wall_frac,building_type)
3764	surf_usm_v(l)%frac(m,ind_pav_green) = building_pars(ind_green_frac_w,building_type)
376