Home

Context Navigation

urban_surface_mod.f90 @ 4238

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