Home

Context Navigation

urban_surface_mod.f90 @ 4441

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