Home

Context Navigation

urban_surface_mod.f90 @ 4148

Last change on this file since 4148 was 4148, checked in by suehring, 22 months ago
indoor_model_mod: Revision of indoor model; urban_surface_mod: parameters for waste heat from cooling and heating are introduced to building data base; initialization of building data base moved to an earlier point of time before indoor model will be initialized; radiation_model_mod: minor improvement in some comments; synthetic_turbulence_generator: unused variable removed
Property svn:keywords set to `Id`
File size: 549.6 KB

Line
1	!> @file urban_surface_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2015-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 4148 2019-08-08 11:26:00Z suehring $
30	! - Add anthropogenic heat output factors for heating and cooling to building
31	! data base
32	! - Move definition of building_pars to usm_init_arrays since it is already
33	! required in the indoor model
34	!
35	! 4127 2019-07-30 14:47:10Z suehring
36	! Do not add anthopogenic energy during wall/soil spin-up
37	! (merge from branch resler)
38	!
39	! 4077 2019-07-09 13:27:11Z gronemeier
40	! Set roughness length z0 and z0h/q at ground-floor level to same value as
41	! those above ground-floor level
42	!
43	! 4051 2019-06-24 13:58:30Z suehring
44	! Remove work-around for green surface fraction on buildings
45	! (do not set it zero)
46	!
47	! 4050 2019-06-24 13:57:27Z suehring
48	! In order to avoid confusion with global control parameter, rename the
49	! USM-internal flag spinup into during_spinup.
50	!
51	! 3987 2019-05-22 09:52:13Z kanani
52	! Introduce alternative switch for debug output during timestepping
53	!
54	! 3943 2019-05-02 09:50:41Z maronga
55	! Removed qsws_eb. Bugfix in calculation of qsws.
56	!
57	! 3933 2019-04-25 12:33:20Z kanani
58	! Remove allocation of pt_2m, this is done in surface_mod now (surfaces%pt_2m)
59	!
60	! 3921 2019-04-18 14:21:10Z suehring
61	! Undo accidentally commented initialization
62	!
63	! 3918 2019-04-18 13:33:11Z suehring
64	! Set green fraction to zero also at vertical surfaces
65	!
66	! 3914 2019-04-17 16:02:02Z suehring
67	! In order to obtain correct surface temperature during spinup set window
68	! fraction to zero (only during spinup) instead of just disabling
69	! time-integration of window-surface temperature.
70	!
71	! 3901 2019-04-16 16:17:02Z suehring
72	! Workaround - set green fraction to zero ( green-heat model crashes ).
73	!
74	! 3896 2019-04-15 10:10:17Z suehring
75	!
76	!
77	! 3896 2019-04-15 10:10:17Z suehring
78	! Bugfix, wrong index used for accessing building_pars from PIDS
79	!
80	! 3885 2019-04-11 11:29:34Z kanani
81	! Changes related to global restructuring of location messages and introduction
82	! of additional debug messages
83	!
84	! 3882 2019-04-10 11:08:06Z suehring
85	! Avoid different type kinds
86	! Move definition of building-surface properties from declaration block
87	! to an extra routine
88	!
89	! 3881 2019-04-10 09:31:22Z suehring
90	! Revise determination of local ground-floor level height.
91	! Make level 3 initalization conform with Palm-input-data standard
92	! Move output of albedo and emissivity to radiation module
93	!
94	! 3832 2019-03-28 13:16:58Z raasch
95	! instrumented with openmp directives
96	!
97	! 3824 2019-03-27 15:56:16Z pavelkrc
98	! Remove unused imports
99	!
100	!
101	! 3814 2019-03-26 08:40:31Z pavelkrc
102	! unused subroutine commented out
103	!
104	! 3769 2019-02-28 10:16:49Z moh.hefny
105	! removed unused variables
106	!
107	! 3767 2019-02-27 08:18:02Z raasch
108	! unused variables removed from rrd-subroutines parameter list
109	!
110	! 3748 2019-02-18 10:38:31Z suehring
111	! Revise conversion of waste-heat flux (do not divide by air density, will
112	! be done in diffusion_s)
113	!
114	! 3745 2019-02-15 18:57:56Z suehring
115	! - Remove internal flag indoor_model (is a global control parameter)
116	! - add waste heat from buildings to the kinmatic heat flux
117	! - consider waste heat in restart data
118	! - remove unused USE statements
119	!
120	! 3744 2019-02-15 18:38:58Z suehring
121	! fixed surface heat capacity in the building parameters
122	! convert the file back to unix format
123	!
124	! 3730 2019-02-11 11:26:47Z moh.hefny
125	! Formatting and clean-up (rvtils)
126	!
127	! 3710 2019-01-30 18:11:19Z suehring
128	! Check if building type is set within a valid range.
129	!
130	! 3705 2019-01-29 19:56:39Z suehring
131	! make nzb_wall public, required for virtual-measurements
132	!
133	! 3704 2019-01-29 19:51:41Z suehring
134	! Some interface calls moved to module_interface + cleanup
135	!
136	! 3655 2019-01-07 16:51:22Z knoop
137	! Implementation of the PALM module interface
138	!
139	! 3636 2018-12-19 13:48:34Z raasch
140	! nopointer option removed
141	!
142	! 3614 2018-12-10 07:05:46Z raasch
143	! unused variables removed
144	!
145	! 3607 2018-12-07 11:56:58Z suehring
146	! Output of radiation-related quantities migrated to radiation_model_mod.
147	!
148	! 3597 2018-12-04 08:40:18Z maronga
149	! Fixed calculation method of near surface air potential temperature at 10 cm
150	! and moved to surface_layer_fluxes. Removed unnecessary _eb strings.
151	!
152	! 3524 2018-11-14 13:36:44Z raasch
153	! bugfix concerning allocation of t_surf_wall_v
154	!
155	! 3502 2018-11-07 14:45:23Z suehring
156	! Disable initialization of building roofs with ground-floor-level properties,
157	! since this causes strong oscillations of surface temperature during the
158	! spinup.
159	!
160	! 3469 2018-10-30 20:05:07Z kanani
161	! Add missing PUBLIC variables for new indoor model
162	!
163	! 3449 2018-10-29 19:36:56Z suehring
164	! Bugfix: Fix average arrays allocations in usm_3d_data_averaging (J.Resler)
165	! Bugfix: Fix reading wall temperatures (J.Resler)
166	! Bugfix: Fix treating of outputs for wall temperature and sky view factors (J.Resler)
167	!
168	!
169	! 3435 2018-10-26 18:25:44Z gronemeier
170	! Bugfix: allocate gamma_w_green_sat until nzt_wall+1
171	!
172	! 3418 2018-10-24 16:07:39Z kanani
173	! (rvtils, srissman)
174	! -Updated building databse, two green roof types (ind_green_type_roof)
175	! -Latent heat flux for green walls and roofs, new output of latent heatflux
176	! and soil water content of green roof substrate
177	! -t_surf changed to t_surf_wall
178	! -Added namelist parameter usm_wall_mod for lower wall tendency
179	! of first two wall layers during spinup
180	! -Window calculations deactivated during spinup
181	!
182	! 3382 2018-10-19 13:10:32Z knoop
183	! Bugix: made array declaration Fortran Standard conform
184	!
185	! 3378 2018-10-19 12:34:59Z kanani
186	! merge from radiation branch (r3362) into trunk
187	! (moh.hefny):
188	! - check the requested output variables if they are correct
189	! - added unscheduled_radiation_calls switch to control force_radiation_call
190	! - minor formate changes
191	!
192	! 3371 2018-10-18 13:40:12Z knoop
193	! Set flag indicating that albedo at urban surfaces is already initialized
194	!
195	! 3347 2018-10-15 14:21:08Z suehring
196	! Enable USM initialization with default building parameters in case no static
197	! input file exist.
198	!
199	! 3343 2018-10-15 10:38:52Z suehring
200	! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
201	!
202	! 3274 2018-09-24 15:42:55Z knoop
203	! Modularization of all bulk cloud physics code components
204	!
205	! 3248 2018-09-14 09:42:06Z sward
206	! Minor formating changes
207	!
208	! 3246 2018-09-13 15:14:50Z sward
209	! Added error handling for input namelist via parin_fail_message
210	!
211	! 3241 2018-09-12 15:02:00Z raasch
212	! unused variables removed
213	!
214	! 3223 2018-08-30 13:48:17Z suehring
215	! Bugfix for commit 3222
216	!
217	! 3222 2018-08-30 13:35:35Z suehring
218	! Introduction of surface array for type and its name
219	!
220	! 3203 2018-08-23 10:48:36Z suehring
221	! Revise bulk parameter for emissivity at ground-floor level
222	!
223	! 3196 2018-08-13 12:26:14Z maronga
224	! Added maximum aerodynamic resistance of 300 for horiztonal surfaces.
225	!
226	! 3176 2018-07-26 17:12:48Z suehring
227	! Bugfix, update virtual potential surface temparture, else heat fluxes on
228	! roofs might become unphysical
229	!
230	! 3152 2018-07-19 13:26:52Z suehring
231	! Initialize q_surface, which might be used in surface_layer_fluxes
232	!
233	! 3151 2018-07-19 08:45:38Z raasch
234	! remaining preprocessor define strings __check removed
235	!
236	! 3136 2018-07-16 14:48:21Z suehring
237	! Limit also roughness length for heat and moisture where necessary
238	!
239	! 3123 2018-07-12 16:21:53Z suehring
240	! Correct working precision for INTEGER number
241	!
242	! 3115 2018-07-10 12:49:26Z suehring
243	! Additional building type to represent bridges
244	!
245	! 3091 2018-06-28 16:20:35Z suehring
246	! - Limit aerodynamic resistance at vertical walls.
247	! - Add check for local roughness length not exceeding surface-layer height and
248	! limit roughness length where necessary.
249	!
250	! 3065 2018-06-12 07:03:02Z Giersch
251	! Unused array dxdir was removed, dz was replaced by dzu to consider vertical
252	! grid stretching
253	!
254	! 3049 2018-05-29 13:52:36Z Giersch
255	! Error messages revised
256	!
257	! 3045 2018-05-28 07:55:41Z Giersch
258	! Error message added
259	!
260	! 3029 2018-05-23 12:19:17Z raasch
261	! bugfix: close unit 151 instead of 90
262	!
263	! 3014 2018-05-09 08:42:38Z maronga
264	! Added pc_transpiration_rate
265	!
266	! 2977 2018-04-17 10:27:57Z kanani
267	! Implement changes from branch radiation (r2948-2971) with minor modifications.
268	! (moh.hefny):
269	! Extended exn for all model domain height to avoid the need to get nzut.
270	!
271	! 2963 2018-04-12 14:47:44Z suehring
272	! Introduce index for vegetation/wall, pavement/green-wall and water/window
273	! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
274	!
275	! 2943 2018-04-03 16:17:10Z suehring
276	! Calculate exner function at all height levels and remove some un-used
277	! variables.
278	!
279	! 2932 2018-03-26 09:39:22Z maronga
280	! renamed urban_surface_par to urban_surface_parameters
281	!
282	! 2921 2018-03-22 15:05:23Z Giersch
283	! The activation of spinup has been moved to parin
284	!
285	! 2920 2018-03-22 11:22:01Z kanani
286	! Remove unused pcbl, npcbl from ONLY list
287	! moh.hefny:
288	! Fixed bugs introduced by new structures and by moving radiation interaction
289	! into radiation_model_mod.f90.
290	! Bugfix: usm data output 3D didn't respect directions
291	!
292	! 2906 2018-03-19 08:56:40Z Giersch
293	! Local variable ids has to be initialized with a value of -1 in
294	! usm_3d_data_averaging
295	!
296	! 2894 2018-03-15 09:17:58Z Giersch
297	! Calculations of the index range of the subdomain on file which overlaps with
298	! the current subdomain are already done in read_restart_data_mod,
299	! usm_read/write_restart_data have been renamed to usm_r/wrd_local, variable
300	! named found has been introduced for checking if restart data was found,
301	! reading of restart strings has been moved completely to
302	! read_restart_data_mod, usm_rrd_local is already inside the overlap loop
303	! programmed in read_restart_data_mod, SAVE attribute added where necessary,
304	! deallocation and allocation of some arrays have been changed to take care of
305	! different restart files that can be opened (index i), the marker *** end usm
306	! *** is not necessary anymore, strings and their respective lengths are
307	! written out and read now in case of restart runs to get rid of prescribed
308	! character lengths
309	!
310	! 2805 2018-02-14 17:00:09Z suehring
311	! Initialization of resistances.
312	!
313	! 2797 2018-02-08 13:24:35Z suehring
314	! Comment concerning output of ground-heat flux added.
315	!
316	! 2766 2018-01-22 17:17:47Z kanani
317	! Removed redundant commas, added some blanks
318	!
319	! 2765 2018-01-22 11:34:58Z maronga
320	! Major bugfix in calculation of f_shf. Adjustment of roughness lengths in
321	! building_pars
322	!
323	! 2750 2018-01-15 16:26:51Z knoop
324	! Move flag plant canopy to modules
325	!
326	! 2737 2018-01-11 14:58:11Z kanani
327	! Removed unused variables t_surf_whole...
328	!
329	! 2735 2018-01-11 12:01:27Z suehring
330	! resistances are saved in surface attributes
331	!
332	! 2723 2018-01-05 09:27:03Z maronga
333	! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
334	!
335	! 2720 2018-01-02 16:27:15Z kanani
336	! Correction of comment
337	!
338	! 2718 2018-01-02 08:49:38Z maronga
339	! Corrected "Former revisions" section
340	!
341	! 2705 2017-12-18 11:26:23Z maronga
342	! Changes from last commit documented
343	!
344	! 2703 2017-12-15 20:12:38Z maronga
345	! Workaround for calculation of r_a
346	!
347	! 2696 2017-12-14 17:12:51Z kanani
348	! - Change in file header (GPL part)
349	! - Bugfix in calculation of pt_surface and related fluxes. (BM)
350	! - Do not write surface temperatures onto pt array as this might cause
351	! problems with nesting. (MS)
352	! - Revised calculation of pt1 (now done in surface_layer_fluxes).
353	! Bugfix, f_shf_window and f_shf_green were not set at vertical surface
354	! elements. (MS)
355	! - merged with branch ebsolver
356	! green building surfaces do not evaporate yet
357	! properties of green wall layers and window layers are taken from wall layers
358	! this input data is missing. (RvT)
359	! - Merged with branch radiation (developed by Mohamed Salim)
360	! - Revised initialization. (MS)
361	! - Rename emiss_surf into emissivity, roughness_wall into z0, albedo_surf into
362	! albedo. (MS)
363	! - Move first call of usm_radiatin from usm_init to init_3d_model
364	! - fixed problem with near surface temperature
365	! - added near surface temperature pt_10cm_h(m), pt_10cm_v(l)%t(m)
366	! - does not work with temp profile including stability, ol
367	! pt_10cm = pt1 now
368	! - merged with 2357 bugfix, error message for nopointer version
369	! - added indoor model coupling with wall heat flux
370	! - added green substrate/ dry vegetation layer for buildings
371	! - merged with 2232 new surface-type structure
372	! - added transmissivity of window tiles
373	! - added MOSAIK tile approach for 3 different surfaces (RvT)
374	!
375	! 2583 2017-10-26 13:58:38Z knoop
376	! Bugfix: reverted MPI_Win_allocate_cptr introduction in last commit
377	!
378	! 2582 2017-10-26 13:19:46Z hellstea
379	! Workaround for gnufortran compiler added in usm_calc_svf. CALL MPI_Win_allocate is
380	! replaced by CALL MPI_Win_allocate_cptr if defined ( __gnufortran ).
381	!
382	! 2544 2017-10-13 18:09:32Z maronga
383	! Date and time quantities are now read from date_and_time_mod. Solar constant is
384	! read from radiation_model_mod
385	!
386	! 2516 2017-10-04 11:03:04Z suehring
387	! Remove tabs
388	!
389	! 2514 2017-10-04 09:52:37Z suehring
390	! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
391	! no output of ghost layer data
392	!
393	! 2350 2017-08-15 11:48:26Z kanani
394	! Bugfix and error message for nopointer version.
395	! Additional "! defined(__nopointer)" as workaround to enable compilation of
396	! nopointer version.
397	!
398	! 2318 2017-07-20 17:27:44Z suehring
399	! Get topography top index via Function call
400	!
401	! 2317 2017-07-20 17:27:19Z suehring
402	! Bugfix: adjust output of shf. Added support for spinups
403	!
404	! 2287 2017-06-15 16:46:30Z suehring
405	! Bugfix in determination topography-top index
406	!
407	! 2269 2017-06-09 11:57:32Z suehring
408	! Enable restart runs with different number of PEs
409	! Bugfixes nopointer branch
410	!
411	! 2258 2017-06-08 07:55:13Z suehring
412	! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
413	!
414	! 2233 2017-05-30 18:08:54Z suehring
415	!
416	! 2232 2017-05-30 17:47:52Z suehring
417	! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
418	! insteat, heat fluxes are directly applied in diffusion_s.
419	!
420	! 2213 2017-04-24 15:10:35Z kanani
421	! Removal of output quantities usm_lad and usm_canopy_hr
422	!
423	! 2209 2017-04-19 09:34:46Z kanani
424	! cpp switch __mpi3 removed,
425	! minor formatting,
426	! small bugfix for division by zero (Krc)
427	!
428	! 2113 2017-01-12 13:40:46Z kanani
429	! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
430	!
431	! 2071 2016-11-17 11:22:14Z maronga
432	! Small bugfix (Resler)
433	!
434	! 2031 2016-10-21 15:11:58Z knoop
435	! renamed variable rho to rho_ocean
436	!
437	! 2024 2016-10-12 16:42:37Z kanani
438	! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
439	! optimization of MPI-RMA operations,
440	! declaration of pcbl as integer,
441	! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
442	! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
443	! use of new control parameter varnamelength,
444	! added output variables usm_rad_ressw, usm_rad_reslw,
445	! minor formatting changes,
446	! minor optimizations.
447	!
448	! 2011 2016-09-19 17:29:57Z kanani
449	! Major reformatting according to PALM coding standard (comments, blanks,
450	! alphabetical ordering, etc.),
451	! removed debug_prints,
452	! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
453	! defined in MODULE control_parameters (modules.f90) to avoid circular
454	! dependencies,
455	! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
456	!
457	! 2007 2016-08-24 15:47:17Z kanani
458	! Initial revision
459	!
460	!
461	! Description:
462	! ------------
463	! 2016/6/9 - Initial version of the USM (Urban Surface Model)
464	! authors: Jaroslav Resler, Pavel Krc
465	! (Czech Technical University in Prague and Institute of
466	! Computer Science of the Czech Academy of Sciences, Prague)
467	! with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
468	! partly inspired by PALM LSM (B. Maronga)
469	! parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
470	!> Module for Urban Surface Model (USM)
471	!> The module includes:
472	!> 1. radiation model with direct/diffuse radiation, shading, reflections
473	!> and integration with plant canopy
474	!> 2. wall and wall surface model
475	!> 3. surface layer energy balance
476	!> 4. anthropogenic heat (only from transportation so far)
477	!> 5. necessary auxiliary subroutines (reading inputs, writing outputs,
478	!> restart simulations, ...)
479	!> It also make use of standard radiation and integrates it into
480	!> urban surface model.
481	!>
482	!> Further work:
483	!> -------------
484	!> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity
485	!> from surfaces that are visible from local surfaces (i.e. there is a SVF
486	!> where target is local). To do that, radiosity will be exchanged after each
487	!> reflection step using MPI_Alltoall instead of current MPI_Allgather.
488	!>
489	!> 2. Temporarily large values of surface heat flux can be observed, up to
490	!> 1.2 Km/s, which seem to be not realistic.
491	!>
492	!> @todo Output of _av variables in case of restarts
493	!> @todo Revise flux conversion in energy-balance solver
494	!> @todo Check optimizations for RMA operations
495	!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
496	!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
497	!> factor 3 between min and max time
498	!> @todo Check divisions in wtend (etc.) calculations for possible division
499	!> by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
500	!> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
501	!> @todo Move plant canopy stuff into plant canopy code
502	!------------------------------------------------------------------------------!
503	MODULE urban_surface_mod
504
505	USE arrays_3d, &
506	ONLY: hyp, zu, pt, p, u, v, w, tend, exner, hyrho, prr, q, ql, vpt
507
508	USE calc_mean_profile_mod, &
509	ONLY: calc_mean_profile
510
511	USE basic_constants_and_equations_mod, &
512	ONLY: c_p, g, kappa, pi, r_d, rho_l, l_v, sigma_sb
513
514	USE control_parameters, &
515	ONLY: coupling_start_time, topography, &
516	debug_output, debug_output_timestep, debug_string, &
517	dt_3d, humidity, indoor_model, &
518	intermediate_timestep_count, initializing_actions, &
519	intermediate_timestep_count_max, simulated_time, end_time, &
520	timestep_scheme, tsc, coupling_char, io_blocks, io_group, &
521	message_string, time_since_reference_point, surface_pressure, &
522	pt_surface, large_scale_forcing, lsf_surf, &
523	spinup_pt_mean, spinup_time, time_do3d, dt_do3d, &
524	average_count_3d, varnamelength, urban_surface, dz
525
526	USE bulk_cloud_model_mod, &
527	ONLY: bulk_cloud_model, precipitation
528
529	USE cpulog, &
530	ONLY: cpu_log, log_point, log_point_s
531
532	USE date_and_time_mod, &
533	ONLY: time_utc_init
534
535	USE grid_variables, &
536	ONLY: dx, dy, ddx, ddy, ddx2, ddy2
537
538	USE indices, &
539	ONLY: nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, &
540	nysg, nzb, nzt, nbgp, wall_flags_0
541
542	USE, INTRINSIC :: iso_c_binding
543
544	USE kinds
545
546	USE pegrid
547
548	USE radiation_model_mod, &
549	ONLY: albedo_type, radiation_interaction, &
550	radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out, &
551	force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u, &
552	iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id, &
553	iz, iy, ix, nsurf, idsvf, ndsvf, &
554	idcsf, ndcsf, kdcsf, pct, &
555	nz_urban_b, nz_urban_t, unscheduled_radiation_calls
556
557	USE statistics, &
558	ONLY: hom, statistic_regions
559
560	USE surface_mod, &
561	ONLY: get_topography_top_index_ji, get_topography_top_index, &
562	ind_pav_green, ind_veg_wall, ind_wat_win, surf_usm_h, &
563	surf_usm_v, surface_restore_elements
564
565
566	IMPLICIT NONE
567
568	!
569	!-- USM model constants
570
571	REAL(wp), PARAMETER :: &
572	b_ch = 6.04_wp, & !< Clapp & Hornberger exponent
573	lambda_h_green_dry = 0.19_wp, & !< heat conductivity for dry soil
574	lambda_h_green_sm = 3.44_wp, & !< heat conductivity of the soil matrix
575	lambda_h_water = 0.57_wp, & !< heat conductivity of water
576	psi_sat = -0.388_wp, & !< soil matrix potential at saturation
577	rho_c_soil = 2.19E6_wp, & !< volumetric heat capacity of soil
578	rho_c_water = 4.20E6_wp !< volumetric heat capacity of water
579	! m_max_depth = 0.0002_wp ! Maximum capacity of the water reservoir (m)
580
581	!
582	!-- Soil parameters I alpha_vg, l_vg_green, n_vg, gamma_w_green_sat
583	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/ &
584	3.83_wp, 1.250_wp, 1.38_wp, 6.94E-6_wp, & !< soil 1
585	3.14_wp, -2.342_wp, 1.28_wp, 1.16E-6_wp, & !< soil 2
586	0.83_wp, -0.588_wp, 1.25_wp, 0.26E-6_wp, & !< soil 3
587	3.67_wp, -1.977_wp, 1.10_wp, 2.87E-6_wp, & !< soil 4
588	2.65_wp, 2.500_wp, 1.10_wp, 1.74E-6_wp, & !< soil 5
589	1.30_wp, 0.400_wp, 1.20_wp, 0.93E-6_wp, & !< soil 6
590	0.00_wp, 0.00_wp, 0.00_wp, 0.57E-6_wp & !< soil 7
591	/), (/ 4, 7 /) )
592
593	!
594	!-- Soil parameters II swc_sat, fc, wilt, swc_res
595	REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
596	0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & !< soil 1
597	0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & !< soil 2
598	0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & !< soil 3
599	0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & !< soil 4
600	0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & !< soil 5
601	0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & !< soil 6
602	0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp & !< soil 7
603	/), (/ 4, 7 /) )
604	!
605	!-- value 9999999.9_wp -> generic available or user-defined value must be set
606	!-- otherwise -> no generic variable and user setting is optional
607	REAL(wp) :: alpha_vangenuchten = 9999999.9_wp, & !< NAMELIST alpha_vg
608	field_capacity = 9999999.9_wp, & !< NAMELIST fc
609	hydraulic_conductivity = 9999999.9_wp, & !< NAMELIST gamma_w_green_sat
610	l_vangenuchten = 9999999.9_wp, & !< NAMELIST l_vg
611	n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg
612	residual_moisture = 9999999.9_wp, & !< NAMELIST m_res
613	saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat
614	wilting_point = 9999999.9_wp !< NAMELIST m_wilt
615
616	!
617	!-- configuration parameters (they can be setup in PALM config)
618	LOGICAL :: usm_material_model = .TRUE. !< flag parameter indicating wheather the model of heat in materials is used
619	LOGICAL :: usm_anthropogenic_heat = .FALSE. !< flag parameter indicating wheather the anthropogenic heat sources
620	!< (e.g.transportation) are used
621	LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls
622	LOGICAL :: read_wall_temp_3d = .FALSE.
623	LOGICAL :: usm_wall_mod = .FALSE. !< reduces conductivity of the first 2 wall layers by factor 0.1
624
625
626	INTEGER(iwp) :: building_type = 1 !< default building type (preleminary setting)
627	INTEGER(iwp) :: land_category = 2 !< default category for land surface
628	INTEGER(iwp) :: wall_category = 2 !< default category for wall surface over pedestrian zone
629	INTEGER(iwp) :: pedestrian_category = 2 !< default category for wall surface in pedestrian zone
630	INTEGER(iwp) :: roof_category = 2 !< default category for root surface
631	REAL(wp) :: roughness_concrete = 0.001_wp !< roughness length of average concrete surface
632	!
633	!-- Indices of input attributes in building_pars for (above) ground floor level
634	INTEGER(iwp) :: ind_alb_wall_agfl = 38 !< index in input list for albedo_type of wall above ground floor level
635	INTEGER(iwp) :: ind_alb_wall_gfl = 66 !< index in input list for albedo_type of wall ground floor level
636	INTEGER(iwp) :: ind_alb_wall_r = 101 !< index in input list for albedo_type of wall roof
637	INTEGER(iwp) :: ind_alb_green_agfl = 39 !< index in input list for albedo_type of green above ground floor level
638	INTEGER(iwp) :: ind_alb_green_gfl = 78 !< index in input list for albedo_type of green ground floor level
639	INTEGER(iwp) :: ind_alb_green_r = 117 !< index in input list for albedo_type of green roof
640	INTEGER(iwp) :: ind_alb_win_agfl = 40 !< index in input list for albedo_type of window fraction above ground floor level
641	INTEGER(iwp) :: ind_alb_win_gfl = 77 !< index in input list for albedo_type of window fraction ground floor level
642	INTEGER(iwp) :: ind_alb_win_r = 115 !< index in input list for albedo_type of window fraction roof
643	INTEGER(iwp) :: ind_c_surface = 45 !< index in input list for heat capacity wall surface
644	INTEGER(iwp) :: ind_c_surface_green = 48 !< index in input list for heat capacity green surface
645	INTEGER(iwp) :: ind_c_surface_win = 47 !< index in input list for heat capacity window surface
646	INTEGER(iwp) :: ind_emis_wall_agfl = 14 !< index in input list for wall emissivity, above ground floor level
647	INTEGER(iwp) :: ind_emis_wall_gfl = 32 !< index in input list for wall emissivity, ground floor level
648	INTEGER(iwp) :: ind_emis_wall_r = 100 !< index in input list for wall emissivity, roof
649	INTEGER(iwp) :: ind_emis_green_agfl = 15 !< index in input list for green emissivity, above ground floor level
650	INTEGER(iwp) :: ind_emis_green_gfl = 34 !< index in input list for green emissivity, ground floor level
651	INTEGER(iwp) :: ind_emis_green_r = 116 !< index in input list for green emissivity, roof
652	INTEGER(iwp) :: ind_emis_win_agfl = 16 !< index in input list for window emissivity, above ground floor level
653	INTEGER(iwp) :: ind_emis_win_gfl = 33 !< index in input list for window emissivity, ground floor level
654	INTEGER(iwp) :: ind_emis_win_r = 113 !< index in input list for window emissivity, roof
655	INTEGER(iwp) :: ind_gflh = 20 !< index in input list for ground floor level height
656	INTEGER(iwp) :: ind_green_frac_w_agfl = 2 !< index in input list for green fraction on wall, above ground floor level
657	INTEGER(iwp) :: ind_green_frac_w_gfl = 23 !< index in input list for green fraction on wall, ground floor level
658	INTEGER(iwp) :: ind_green_frac_r_agfl = 3 !< index in input list for green fraction on roof, above ground floor level
659	INTEGER(iwp) :: ind_green_frac_r_gfl = 24 !< index in input list for green fraction on roof, ground floor level
660	INTEGER(iwp) :: ind_hc1_agfl = 6 !< index in input list for heat capacity at first wall layer,
661	!< above ground floor level
662	INTEGER(iwp) :: ind_hc1_gfl = 26 !< index in input list for heat capacity at first wall layer, ground floor level
663	INTEGER(iwp) :: ind_hc1_wall_r = 94 !< index in input list for heat capacity at first wall layer, roof
664	INTEGER(iwp) :: ind_hc1_win_agfl = 83 !< index in input list for heat capacity at first window layer,
665	!< above ground floor level
666	INTEGER(iwp) :: ind_hc1_win_gfl = 71 !< index in input list for heat capacity at first window layer,
667	!< ground floor level
668	INTEGER(iwp) :: ind_hc1_win_r = 107 !< index in input list for heat capacity at first window layer, roof
669	INTEGER(iwp) :: ind_hc2_agfl = 7 !< index in input list for heat capacity at second wall layer,
670	!< above ground floor level
671	INTEGER(iwp) :: ind_hc2_gfl = 27 !< index in input list for heat capacity at second wall layer, ground floor level
672	INTEGER(iwp) :: ind_hc2_wall_r = 95 !< index in input list for heat capacity at second wall layer, roof
673	INTEGER(iwp) :: ind_hc2_win_agfl = 84 !< index in input list for heat capacity at second window layer,
674	!< above ground floor level
675	INTEGER(iwp) :: ind_hc2_win_gfl = 72 !< index in input list for heat capacity at second window layer,
676	!< ground floor level
677	INTEGER(iwp) :: ind_hc2_win_r = 108 !< index in input list for heat capacity at second window layer, roof
678	INTEGER(iwp) :: ind_hc3_agfl = 8 !< index in input list for heat capacity at third wall layer,
679	!< above ground floor level
680	INTEGER(iwp) :: ind_hc3_gfl = 28 !< index in input list for heat capacity at third wall layer, ground floor level
681	INTEGER(iwp) :: ind_hc3_wall_r = 96 !< index in input list for heat capacity at third wall layer, roof
682	INTEGER(iwp) :: ind_hc3_win_agfl = 85 !< index in input list for heat capacity at third window layer,
683	!< above ground floor level
684	INTEGER(iwp) :: ind_hc3_win_gfl = 73 !< index in input list for heat capacity at third window layer,
685	!< ground floor level
686	INTEGER(iwp) :: ind_hc3_win_r = 109 !< index in input list for heat capacity at third window layer, roof
687	INTEGER(iwp) :: ind_indoor_target_temp_summer = 12
688	INTEGER(iwp) :: ind_indoor_target_temp_winter = 13
689	INTEGER(iwp) :: ind_lai_r_agfl = 4 !< index in input list for LAI on roof, above ground floor level
690	INTEGER(iwp) :: ind_lai_r_gfl = 4 !< index in input list for LAI on roof, ground floor level
691	INTEGER(iwp) :: ind_lai_w_agfl = 5 !< index in input list for LAI on wall, above ground floor level
692	INTEGER(iwp) :: ind_lai_w_gfl = 25 !< index in input list for LAI on wall, ground floor level
693	INTEGER(iwp) :: ind_lambda_surf = 46 !< index in input list for thermal conductivity of wall surface
694	INTEGER(iwp) :: ind_lambda_surf_green = 50 !< index in input list for thermal conductivity of green surface
695	INTEGER(iwp) :: ind_lambda_surf_win = 49 !< index in input list for thermal conductivity of window surface
696	INTEGER(iwp) :: ind_tc1_agfl = 9 !< index in input list for thermal conductivity at first wall layer,
697	!< above ground floor level
698	INTEGER(iwp) :: ind_tc1_gfl = 29 !< index in input list for thermal conductivity at first wall layer,
699	!< ground floor level
700	INTEGER(iwp) :: ind_tc1_wall_r = 97 !< index in input list for thermal conductivity at first wall layer, roof
701	INTEGER(iwp) :: ind_tc1_win_agfl = 86 !< index in input list for thermal conductivity at first window layer,
702	!< above ground floor level
703	INTEGER(iwp) :: ind_tc1_win_gfl = 74 !< index in input list for thermal conductivity at first window layer,
704	!< ground floor level
705	INTEGER(iwp) :: ind_tc1_win_r = 110 !< index in input list for thermal conductivity at first window layer, roof
706	INTEGER(iwp) :: ind_tc2_agfl = 10 !< index in input list for thermal conductivity at second wall layer,
707	!< above ground floor level
708	INTEGER(iwp) :: ind_tc2_gfl = 30 !< index in input list for thermal conductivity at second wall layer,
709	!< ground floor level
710	INTEGER(iwp) :: ind_tc2_wall_r = 98 !< index in input list for thermal conductivity at second wall layer, roof
711	INTEGER(iwp) :: ind_tc2_win_agfl = 87 !< index in input list for thermal conductivity at second window layer,
712	!< above ground floor level
713	INTEGER(iwp) :: ind_tc2_win_gfl = 75 !< index in input list for thermal conductivity at second window layer,
714	!< ground floor level
715	INTEGER(iwp) :: ind_tc2_win_r = 111 !< index in input list for thermal conductivity at second window layer,
716	!< ground floor level
717	INTEGER(iwp) :: ind_tc3_agfl = 11 !< index in input list for thermal conductivity at third wall layer,
718	!< above ground floor level
719	INTEGER(iwp) :: ind_tc3_gfl = 31 !< index in input list for thermal conductivity at third wall layer,
720	!< ground floor level
721	INTEGER(iwp) :: ind_tc3_wall_r = 99 !< index in input list for thermal conductivity at third wall layer, roof
722	INTEGER(iwp) :: ind_tc3_win_agfl = 88 !< index in input list for thermal conductivity at third window layer,
723	!< above ground floor level
724	INTEGER(iwp) :: ind_tc3_win_gfl = 76 !< index in input list for thermal conductivity at third window layer,
725	!< ground floor level
726	INTEGER(iwp) :: ind_tc3_win_r = 112 !< index in input list for thermal conductivity at third window layer, roof
727	INTEGER(iwp) :: ind_thick_1_agfl = 41 !< index for wall layer thickness - 1st layer above ground floor level
728	INTEGER(iwp) :: ind_thick_1_gfl = 62 !< index for wall layer thickness - 1st layer ground floor level
729	INTEGER(iwp) :: ind_thick_1_wall_r = 90 !< index for wall layer thickness - 1st layer roof
730	INTEGER(iwp) :: ind_thick_1_win_agfl = 79 !< index for window layer thickness - 1st layer above ground floor level
731	INTEGER(iwp) :: ind_thick_1_win_gfl = 67 !< index for window layer thickness - 1st layer ground floor level
732	INTEGER(iwp) :: ind_thick_1_win_r = 103 !< index for window layer thickness - 1st layer roof
733	INTEGER(iwp) :: ind_thick_2_agfl = 42 !< index for wall layer thickness - 2nd layer above ground floor level
734	INTEGER(iwp) :: ind_thick_2_gfl = 63 !< index for wall layer thickness - 2nd layer ground floor level
735	INTEGER(iwp) :: ind_thick_2_wall_r = 91 !< index for wall layer thickness - 2nd layer roof
736	INTEGER(iwp) :: ind_thick_2_win_agfl = 80 !< index for window layer thickness - 2nd layer above ground floor level
737	INTEGER(iwp) :: ind_thick_2_win_gfl = 68 !< index for window layer thickness - 2nd layer ground floor level
738	INTEGER(iwp) :: ind_thick_2_win_r = 104 !< index for window layer thickness - 2nd layer roof
739	INTEGER(iwp) :: ind_thick_3_agfl = 43 !< index for wall layer thickness - 3rd layer above ground floor level
740	INTEGER(iwp) :: ind_thick_3_gfl = 64 !< index for wall layer thickness - 3rd layer ground floor level
741	INTEGER(iwp) :: ind_thick_3_wall_r = 92 !< index for wall layer thickness - 3rd layer roof
742	INTEGER(iwp) :: ind_thick_3_win_agfl = 81 !< index for window layer thickness - 3rd layer above ground floor level
743	INTEGER(iwp) :: ind_thick_3_win_gfl = 69 !< index for window layer thickness - 3rd layer ground floor level
744	INTEGER(iwp) :: ind_thick_3_win_r = 105 !< index for window layer thickness - 3rd layer roof
745	INTEGER(iwp) :: ind_thick_4_agfl = 44 !< index for wall layer thickness - 4th layer above ground floor level
746	INTEGER(iwp) :: ind_thick_4_gfl = 65 !< index for wall layer thickness - 4th layer ground floor level
747	INTEGER(iwp) :: ind_thick_4_wall_r = 93 !< index for wall layer thickness - 4st layer roof
748	INTEGER(iwp) :: ind_thick_4_win_agfl = 82 !< index for window layer thickness - 4th layer above ground floor level
749	INTEGER(iwp) :: ind_thick_4_win_gfl = 70 !< index for window layer thickness - 4th layer ground floor level
750	INTEGER(iwp) :: ind_thick_4_win_r = 106 !< index for window layer thickness - 4th layer roof
751	INTEGER(iwp) :: ind_trans_agfl = 17 !< index in input list for window transmissivity, above ground floor level
752	INTEGER(iwp) :: ind_trans_gfl = 35 !< index in input list for window transmissivity, ground floor level
753	INTEGER(iwp) :: ind_trans_r = 114 !< index in input list for window transmissivity, roof
754	INTEGER(iwp) :: ind_wall_frac_agfl = 0 !< index in input list for wall fraction, above ground floor level
755	INTEGER(iwp) :: ind_wall_frac_gfl = 21 !< index in input list for wall fraction, ground floor level
756	INTEGER(iwp) :: ind_wall_frac_r = 89 !< index in input list for wall fraction, roof
757	INTEGER(iwp) :: ind_win_frac_agfl = 1 !< index in input list for window fraction, above ground floor level
758	INTEGER(iwp) :: ind_win_frac_gfl = 22 !< index in input list for window fraction, ground floor level
759	INTEGER(iwp) :: ind_win_frac_r = 102 !< index in input list for window fraction, roof
760	INTEGER(iwp) :: ind_z0_agfl = 18 !< index in input list for z0, above ground floor level
761	INTEGER(iwp) :: ind_z0_gfl = 36 !< index in input list for z0, ground floor level
762	INTEGER(iwp) :: ind_z0qh_agfl = 19 !< index in input list for z0h / z0q, above ground floor level
763	INTEGER(iwp) :: ind_z0qh_gfl = 37 !< index in input list for z0h / z0q, ground floor level
764	INTEGER(iwp) :: ind_green_type_roof = 118 !< index in input list for type of green roof
765
766
767	REAL(wp) :: roof_height_limit = 4.0_wp !< height for distinguish between land surfaces and roofs
768	REAL(wp) :: ground_floor_level = 4.0_wp !< default ground floor level
769
770
771	CHARACTER(37), DIMENSION(0:7), PARAMETER :: building_type_name = (/ &
772	'user-defined ', & !< type 0
773	'residential - 1950 ', & !< type 1
774	'residential 1951 - 2000 ', & !< type 2
775	'residential 2001 - ', & !< type 3
776	'office - 1950 ', & !< type 4
777	'office 1951 - 2000 ', & !< type 5
778	'office 2001 - ', & !< type 6
779	'bridges ' & !< type 7
780	/)
781
782
783	!
784	!-- Building facade/wall/green/window properties (partly according to PIDS).
785	!-- Initialization of building_pars is outsourced to usm_init_pars. This is
786	!-- needed because of the huge number of attributes given in building_pars
787	!-- (>700), while intel and gfortran compiler have hard limit of continuation
788	!-- lines of 511.
789	REAL(wp), DIMENSION(0:135,1:7) :: building_pars
790	!
791	!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
792	TYPE t_surf_vertical
793	REAL(wp), DIMENSION(:), ALLOCATABLE :: t
794	END TYPE t_surf_vertical
795	!
796	!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
797	TYPE t_wall_vertical
798	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: t
799	END TYPE t_wall_vertical
800
801	TYPE surf_type_usm
802	REAL(wp), DIMENSION(:), ALLOCATABLE :: var_usm_1d !< 1D prognostic variable
803	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_usm_2d !< 2D prognostic variable
804	END TYPE surf_type_usm
805
806	TYPE(surf_type_usm), POINTER :: m_liq_usm_h, & !< liquid water reservoir (m), horizontal surface elements
807	m_liq_usm_h_p !< progn. liquid water reservoir (m), horizontal surface elements
808
809	TYPE(surf_type_usm), TARGET :: m_liq_usm_h_1, & !<
810	m_liq_usm_h_2 !<
811
812	TYPE(surf_type_usm), DIMENSION(:), POINTER :: &
813	m_liq_usm_v, & !< liquid water reservoir (m), vertical surface elements
814	m_liq_usm_v_p !< progn. liquid water reservoir (m), vertical surface elements
815
816	TYPE(surf_type_usm), DIMENSION(0:3), TARGET :: &
817	m_liq_usm_v_1, & !<
818	m_liq_usm_v_2 !<
819
820	TYPE(surf_type_usm), TARGET :: tm_liq_usm_h_m !< liquid water reservoir tendency (m), horizontal surface elements
821	TYPE(surf_type_usm), DIMENSION(0:3), TARGET :: tm_liq_usm_v_m !< liquid water reservoir tendency (m),
822	!< vertical surface elements
823
824	!
825	!-- anthropogenic heat sources
826	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: aheat !< daily average of anthropogenic heat (W/m2)
827	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aheatprof !< diurnal profiles of anthropogenic heat
828	!< for particular layers
829	INTEGER(iwp) :: naheatlayers = 1 !< number of layers of anthropogenic heat
830
831	!
832	!-- wall surface model
833	!-- wall surface model constants
834	INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched)
835	INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched)
836	INTEGER(iwp), PARAMETER :: nzw = 4 !< number of wall layers (fixed for now)
837
838	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
839	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
840	REAL(wp), DIMENSION(nzb_wall:nzt_wall) :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
841	!< normalized soil, wall and roof, window and
842	!<green layer depths (m/m)
843
844	REAL(wp) :: wall_inner_temperature = 295.0_wp !< temperature of the inner wall
845	!< surface (~22 degrees C) (K)
846	REAL(wp) :: roof_inner_temperature = 295.0_wp !< temperature of the inner roof
847	!< surface (~22 degrees C) (K)
848	REAL(wp) :: soil_inner_temperature = 288.0_wp !< temperature of the deep soil
849	!< (~15 degrees C) (K)
850	REAL(wp) :: window_inner_temperature = 295.0_wp !< temperature of the inner window
851	!< surface (~22 degrees C) (K)
852
853	REAL(wp) :: m_total = 0.0_wp !< weighted total water content of the soil (m3/m3)
854	INTEGER(iwp) :: soil_type
855
856	!
857	!-- surface and material model variables for walls, ground, roofs
858	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m)
859	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_window !< normalized window layer depths (m)
860	REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_green !< normalized green layer depths (m)
861
862	REAL(wp), DIMENSION(:), POINTER :: t_surf_wall_h
863	REAL(wp), DIMENSION(:), POINTER :: t_surf_wall_h_p
864	REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h
865	REAL(wp), DIMENSION(:), POINTER :: t_surf_window_h_p
866	REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h
867	REAL(wp), DIMENSION(:), POINTER :: t_surf_green_h_p
868
869	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_wall_h_1
870	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_wall_h_2
871	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_1
872	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_window_h_2
873	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_1
874	REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET :: t_surf_green_h_2
875
876	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_wall_v
877	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_wall_v_p
878	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v
879	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_window_v_p
880	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v
881	TYPE(t_surf_vertical), DIMENSION(:), POINTER :: t_surf_green_v_p
882
883	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_wall_v_1
884	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_wall_v_2
885	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_1
886	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_window_v_2
887	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_1
888	TYPE(t_surf_vertical), DIMENSION(0:3), TARGET :: t_surf_green_v_2
889
890	!
891	!-- Energy balance variables
892	!-- parameters of the land, roof and wall surfaces
893
894	REAL(wp), DIMENSION(:,:), POINTER :: t_wall_h, t_wall_h_p
895	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_wall_h_1, t_wall_h_2
896	REAL(wp), DIMENSION(:,:), POINTER :: t_window_h, t_window_h_p
897	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_window_h_1, t_window_h_2
898	REAL(wp), DIMENSION(:,:), POINTER :: t_green_h, t_green_h_p
899	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_green_h_1, t_green_h_2
900	REAL(wp), DIMENSION(:,:), POINTER :: swc_h, rootfr_h, wilt_h, fc_h, swc_sat_h, swc_h_p, swc_res_h
901	REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: swc_h_1, rootfr_h_1, &
902	wilt_h_1, fc_h_1, swc_sat_h_1, swc_h_2, swc_res_h_1
903
904
905	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_wall_v, t_wall_v_p
906	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_wall_v_1, t_wall_v_2
907	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_window_v, t_window_v_p
908	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_window_v_1, t_window_v_2
909	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: t_green_v, t_green_v_p
910	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: t_green_v_1, t_green_v_2
911	TYPE(t_wall_vertical), DIMENSION(:), POINTER :: swc_v, swc_v_p
912	TYPE(t_wall_vertical), DIMENSION(0:3), TARGET :: swc_v_1, swc_v_2
913
914	!
915	!-- Surface and material parameters classes (surface_type)
916	!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
917	INTEGER(iwp) :: n_surface_types !< number of the wall type categories
918	INTEGER(iwp), PARAMETER :: n_surface_params = 9 !< number of parameters for each type of the wall
919	INTEGER(iwp), PARAMETER :: ialbedo = 1 !< albedo of the surface
920	INTEGER(iwp), PARAMETER :: iemiss = 2 !< emissivity of the surface
921	INTEGER(iwp), PARAMETER :: ilambdas = 3 !< heat conductivity lambda S between surface
922	!< and material ( W m-2 K-1 )
923	INTEGER(iwp), PARAMETER :: irough = 4 !< roughness length z0 for movements
924	INTEGER(iwp), PARAMETER :: iroughh = 5 !< roughness length z0h for scalars
925	!< (heat, humidity,...)
926	INTEGER(iwp), PARAMETER :: icsurf = 6 !< Surface skin layer heat capacity (J m-2 K-1 )
927	INTEGER(iwp), PARAMETER :: ithick = 7 !< thickness of the surface (wall, roof, land) ( m )
928	INTEGER(iwp), PARAMETER :: irhoC = 8 !< volumetric heat capacity rho*C of
929	!< the material ( J m-3 K-1 )
930	INTEGER(iwp), PARAMETER :: ilambdah = 9 !< thermal conductivity lambda H
931	!< of the wall (W m-1 K-1 )
932	CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports)
933	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: surface_type_codes !< codes of wall types
934	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: surface_params !< parameters of wall types
935
936	!
937	!-- interfaces of subroutines accessed from outside of this module
938	INTERFACE usm_3d_data_averaging
939	MODULE PROCEDURE usm_3d_data_averaging
940	END INTERFACE usm_3d_data_averaging
941
942	INTERFACE usm_boundary_condition
943	MODULE PROCEDURE usm_boundary_condition
944	END INTERFACE usm_boundary_condition
945
946	INTERFACE usm_check_data_output
947	MODULE PROCEDURE usm_check_data_output
948	END INTERFACE usm_check_data_output
949
950	INTERFACE usm_check_parameters
951	MODULE PROCEDURE usm_check_parameters
952	END INTERFACE usm_check_parameters
953
954	INTERFACE usm_data_output_3d
955	MODULE PROCEDURE usm_data_output_3d
956	END INTERFACE usm_data_output_3d
957
958	INTERFACE usm_define_netcdf_grid
959	MODULE PROCEDURE usm_define_netcdf_grid
960	END INTERFACE usm_define_netcdf_grid
961
962	INTERFACE usm_init
963	MODULE PROCEDURE usm_init
964	END INTERFACE usm_init
965
966	INTERFACE usm_init_arrays
967	MODULE PROCEDURE usm_init_arrays
968	END INTERFACE usm_init_arrays
969
970	INTERFACE usm_material_heat_model
971	MODULE PROCEDURE usm_material_heat_model
972	END INTERFACE usm_material_heat_model
973
974	INTERFACE usm_green_heat_model
975	MODULE PROCEDURE usm_green_heat_model
976	END INTERFACE usm_green_heat_model
977
978	INTERFACE usm_parin
979	MODULE PROCEDURE usm_parin
980	END INTERFACE usm_parin
981
982	INTERFACE usm_rrd_local
983	MODULE PROCEDURE usm_rrd_local
984	END INTERFACE usm_rrd_local
985
986	INTERFACE usm_surface_energy_balance
987	MODULE PROCEDURE usm_surface_energy_balance
988	END INTERFACE usm_surface_energy_balance
989
990	INTERFACE usm_swap_timelevel
991	MODULE PROCEDURE usm_swap_timelevel
992	END INTERFACE usm_swap_timelevel
993
994	INTERFACE usm_wrd_local
995	MODULE PROCEDURE usm_wrd_local
996	END INTERFACE usm_wrd_local
997
998
999	SAVE
1000
1001	PRIVATE
1002
1003	!
1004	!-- Public functions
1005	PUBLIC usm_boundary_condition, usm_check_parameters, usm_init, &
1006	usm_rrd_local, &
1007	usm_surface_energy_balance, usm_material_heat_model, &
1008	usm_swap_timelevel, usm_check_data_output, usm_3d_data_averaging, &
1009	usm_data_output_3d, usm_define_netcdf_grid, usm_parin, &
1010	usm_wrd_local, usm_init_arrays
1011
1012	!
1013	!-- Public parameters, constants and initial values
1014	PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod, &
1015	usm_green_heat_model, building_pars, &
1016	nzb_wall, nzt_wall, t_wall_h, t_wall_v, &
1017	t_window_h, t_window_v, building_type
1018
1019
1020
1021	CONTAINS
1022
1023	!------------------------------------------------------------------------------!
1024	! Description:
1025	! ------------
1026	!> This subroutine creates the necessary indices of the urban surfaces
1027	!> and plant canopy and it allocates the needed arrays for USM
1028	!------------------------------------------------------------------------------!
1029	SUBROUTINE usm_init_arrays
1030
1031	IMPLICIT NONE
1032
1033	INTEGER(iwp) :: l
1034
1035	IF ( debug_output ) CALL debug_message( 'usm_init_arrays', 'start' )
1036
1037	!
1038	!-- Allocate radiation arrays which are part of the new data type.
1039	!-- For horizontal surfaces.
1040	ALLOCATE ( surf_usm_h%surfhf(1:surf_usm_h%ns) )
1041	ALLOCATE ( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
1042	!
1043	!-- For vertical surfaces
1044	DO l = 0, 3
1045	ALLOCATE ( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns) )
1046	ALLOCATE ( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
1047	ENDDO
1048
1049	!
1050	!-- Wall surface model
1051	!-- allocate arrays for wall surface model and define pointers
1052	!-- allocate array of wall types and wall parameters
1053	ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) )
1054	ALLOCATE ( surf_usm_h%building_type(1:surf_usm_h%ns) )
1055	ALLOCATE ( surf_usm_h%building_type_name(1:surf_usm_h%ns) )
1056	surf_usm_h%building_type = 0
1057	surf_usm_h%building_type_name = 'none'
1058	DO l = 0, 3
1059	ALLOCATE ( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
1060	ALLOCATE ( surf_usm_v(l)%building_type(1:surf_usm_v(l)%ns) )
1061	ALLOCATE ( surf_usm_v(l)%building_type_name(1:surf_usm_v(l)%ns) )
1062	surf_usm_v(l)%building_type = 0
1063	surf_usm_v(l)%building_type_name = 'none'
1064	ENDDO
1065	!
1066	!-- Allocate albedo_type and albedo. Each surface element
1067	!-- has 3 values, 0: wall fraction, 1: green fraction, 2: window fraction.
1068	ALLOCATE ( surf_usm_h%albedo_type(0:2,1:surf_usm_h%ns) )
1069	ALLOCATE ( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) )
1070	surf_usm_h%albedo_type = albedo_type
1071	DO l = 0, 3
1072	ALLOCATE ( surf_usm_v(l)%albedo_type(0:2,1:surf_usm_v(l)%ns) )
1073	ALLOCATE ( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
1074	surf_usm_v(l)%albedo_type = albedo_type
1075	ENDDO
1076
1077	!
1078	!-- Allocate indoor target temperature for summer and winter
1079	ALLOCATE ( surf_usm_h%target_temp_summer(1:surf_usm_h%ns) )
1080	ALLOCATE ( surf_usm_h%target_temp_winter(1:surf_usm_h%ns) )
1081	DO l = 0, 3
1082	ALLOCATE ( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) )
1083	ALLOCATE ( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) )
1084	ENDDO
1085	!
1086	!-- In case the indoor model is applied, allocate memory for waste heat
1087	!-- and indoor temperature.
1088	IF ( indoor_model ) THEN
1089	ALLOCATE ( surf_usm_h%waste_heat(1:surf_usm_h%ns) )
1090	surf_usm_h%waste_heat = 0.0_wp
1091	DO l = 0, 3
1092	ALLOCATE ( surf_usm_v(l)%waste_heat(1:surf_usm_v(l)%ns) )
1093	surf_usm_v(l)%waste_heat = 0.0_wp
1094	ENDDO
1095	ENDIF
1096	!
1097	!-- Allocate flag indicating ground floor level surface elements
1098	ALLOCATE ( surf_usm_h%ground_level(1:surf_usm_h%ns) )
1099	DO l = 0, 3
1100	ALLOCATE ( surf_usm_v(l)%ground_level(1:surf_usm_v(l)%ns) )
1101	ENDDO
1102	!
1103	!-- Allocate arrays for relative surface fraction.
1104	!-- 0 - wall fraction, 1 - green fraction, 2 - window fraction
1105	ALLOCATE ( surf_usm_h%frac(0:2,1:surf_usm_h%ns) )
1106	surf_usm_h%frac = 0.0_wp
1107	DO l = 0, 3
1108	ALLOCATE ( surf_usm_v(l)%frac(0:2,1:surf_usm_v(l)%ns) )
1109	surf_usm_v(l)%frac = 0.0_wp
1110	ENDDO
1111
1112	!
1113	!-- wall and roof surface parameters. First for horizontal surfaces
1114	ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) )
1115	ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) )
1116	ALLOCATE ( surf_usm_h%lambda_surf_window(1:surf_usm_h%ns) )
1117	ALLOCATE ( surf_usm_h%lambda_surf_green(1:surf_usm_h%ns) )
1118	ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns) )
1119	ALLOCATE ( surf_usm_h%c_surface_window(1:surf_usm_h%ns) )
1120	ALLOCATE ( surf_usm_h%c_surface_green(1:surf_usm_h%ns) )
1121	ALLOCATE ( surf_usm_h%transmissivity(1:surf_usm_h%ns) )
1122	ALLOCATE ( surf_usm_h%lai(1:surf_usm_h%ns) )
1123	ALLOCATE ( surf_usm_h%emissivity(0:2,1:surf_usm_h%ns) )
1124	ALLOCATE ( surf_usm_h%r_a(1:surf_usm_h%ns) )
1125	ALLOCATE ( surf_usm_h%r_a_green(1:surf_usm_h%ns) )
1126	ALLOCATE ( surf_usm_h%r_a_window(1:surf_usm_h%ns) )
1127	ALLOCATE ( surf_usm_h%green_type_roof(1:surf_usm_h%ns) )
1128	ALLOCATE ( surf_usm_h%r_s(1:surf_usm_h%ns) )
1129
1130	!
1131	!-- For vertical surfaces.
1132	DO l = 0, 3
1133	ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns) )
1134	ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns) )
1135	ALLOCATE ( surf_usm_v(l)%lambda_surf_window(1:surf_usm_v(l)%ns) )
1136	ALLOCATE ( surf_usm_v(l)%c_surface_window(1:surf_usm_v(l)%ns) )
1137	ALLOCATE ( surf_usm_v(l)%lambda_surf_green(1:surf_usm_v(l)%ns) )
1138	ALLOCATE ( surf_usm_v(l)%c_surface_green(1:surf_usm_v(l)%ns) )
1139	ALLOCATE ( surf_usm_v(l)%transmissivity(1:surf_usm_v(l)%ns) )
1140	ALLOCATE ( surf_usm_v(l)%lai(1:surf_usm_v(l)%ns) )
1141	ALLOCATE ( surf_usm_v(l)%emissivity(0:2,1:surf_usm_v(l)%ns) )
1142	ALLOCATE ( surf_usm_v(l)%r_a(1:surf_usm_v(l)%ns) )
1143	ALLOCATE ( surf_usm_v(l)%r_a_green(1:surf_usm_v(l)%ns) )
1144	ALLOCATE ( surf_usm_v(l)%r_a_window(1:surf_usm_v(l)%ns) )
1145	ALLOCATE ( surf_usm_v(l)%r_s(1:surf_usm_v(l)%ns) )
1146	ENDDO
1147
1148	!
1149	!-- allocate wall and roof material parameters. First for horizontal surfaces
1150	ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) )
1151	ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns) )
1152	ALLOCATE ( surf_usm_h%thickness_green(1:surf_usm_h%ns) )
1153	ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1154	ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1155	ALLOCATE ( surf_usm_h%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1156	ALLOCATE ( surf_usm_h%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1157	ALLOCATE ( surf_usm_h%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1158	ALLOCATE ( surf_usm_h%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1159
1160	ALLOCATE ( surf_usm_h%rho_c_total_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1161	ALLOCATE ( surf_usm_h%n_vg_green(1:surf_usm_h%ns) )
1162	ALLOCATE ( surf_usm_h%alpha_vg_green(1:surf_usm_h%ns) )
1163	ALLOCATE ( surf_usm_h%l_vg_green(1:surf_usm_h%ns) )
1164	ALLOCATE ( surf_usm_h%gamma_w_green_sat(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1165	ALLOCATE ( surf_usm_h%lambda_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1166	ALLOCATE ( surf_usm_h%gamma_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1167	ALLOCATE ( surf_usm_h%tswc_h_m(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1168
1169	!
1170	!-- For vertical surfaces.
1171	DO l = 0, 3
1172	ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns) )
1173	ALLOCATE ( surf_usm_v(l)%thickness_window(1:surf_usm_v(l)%ns) )
1174	ALLOCATE ( surf_usm_v(l)%thickness_green(1:surf_usm_v(l)%ns) )
1175	ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1176	ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1177	ALLOCATE ( surf_usm_v(l)%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1178	ALLOCATE ( surf_usm_v(l)%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1179	ALLOCATE ( surf_usm_v(l)%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1180	ALLOCATE ( surf_usm_v(l)%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1181	ENDDO
1182
1183	!
1184	!-- allocate green wall and roof vegetation and soil parameters. First horizontal surfaces
1185	ALLOCATE ( surf_usm_h%g_d(1:surf_usm_h%ns) )
1186	ALLOCATE ( surf_usm_h%c_liq(1:surf_usm_h%ns) )
1187	ALLOCATE ( surf_usm_h%qsws_liq(1:surf_usm_h%ns) )
1188	ALLOCATE ( surf_usm_h%qsws_veg(1:surf_usm_h%ns) )
1189	ALLOCATE ( surf_usm_h%r_canopy(1:surf_usm_h%ns) )
1190	ALLOCATE ( surf_usm_h%r_canopy_min(1:surf_usm_h%ns) )
1191	ALLOCATE ( surf_usm_h%pt_10cm(1:surf_usm_h%ns) )
1192
1193	!
1194	!-- For vertical surfaces.
1195	DO l = 0, 3
1196	ALLOCATE ( surf_usm_v(l)%g_d(1:surf_usm_v(l)%ns) )
1197	ALLOCATE ( surf_usm_v(l)%c_liq(1:surf_usm_v(l)%ns) )
1198	ALLOCATE ( surf_usm_v(l)%qsws_liq(1:surf_usm_v(l)%ns) )
1199	ALLOCATE ( surf_usm_v(l)%qsws_veg(1:surf_usm_v(l)%ns) )
1200	ALLOCATE ( surf_usm_v(l)%r_canopy(1:surf_usm_v(l)%ns) )
1201	ALLOCATE ( surf_usm_v(l)%r_canopy_min(1:surf_usm_v(l)%ns) )
1202	ALLOCATE ( surf_usm_v(l)%pt_10cm(1:surf_usm_v(l)%ns) )
1203	ENDDO
1204
1205	!
1206	!-- allocate wall and roof layers sizes. For horizontal surfaces.
1207	ALLOCATE ( zwn(nzb_wall:nzt_wall) )
1208	ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1209	ALLOCATE ( zwn_window(nzb_wall:nzt_wall) )
1210	ALLOCATE ( surf_usm_h%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1211	ALLOCATE ( zwn_green(nzb_wall:nzt_wall) )
1212	ALLOCATE ( surf_usm_h%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1213	ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1214	ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1215	ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1216	ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1217	ALLOCATE ( surf_usm_h%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1218	ALLOCATE ( surf_usm_h%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1219	ALLOCATE ( surf_usm_h%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1220	ALLOCATE ( surf_usm_h%zw_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1221	ALLOCATE ( surf_usm_h%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1222	ALLOCATE ( surf_usm_h%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1223	ALLOCATE ( surf_usm_h%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1224	ALLOCATE ( surf_usm_h%zw_green(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1225
1226	!
1227	!-- For vertical surfaces.
1228	DO l = 0, 3
1229	ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1230	ALLOCATE ( surf_usm_v(l)%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1231	ALLOCATE ( surf_usm_v(l)%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1232	ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1233	ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1234	ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1235	ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1236	ALLOCATE ( surf_usm_v(l)%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1237	ALLOCATE ( surf_usm_v(l)%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1238	ALLOCATE ( surf_usm_v(l)%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1239	ALLOCATE ( surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1240	ALLOCATE ( surf_usm_v(l)%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1241	ALLOCATE ( surf_usm_v(l)%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1242	ALLOCATE ( surf_usm_v(l)%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1243	ALLOCATE ( surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1244	ENDDO
1245
1246	!
1247	!-- allocate wall and roof temperature arrays, for horizontal walls
1248	!
1249	!-- Allocate if required. Note, in case of restarts, some of these arrays
1250	!-- might be already allocated.
1251	IF ( .NOT. ALLOCATED( t_surf_wall_h_1 ) ) &
1252	ALLOCATE ( t_surf_wall_h_1(1:surf_usm_h%ns) )
1253	IF ( .NOT. ALLOCATED( t_surf_wall_h_2 ) ) &
1254	ALLOCATE ( t_surf_wall_h_2(1:surf_usm_h%ns) )
1255	IF ( .NOT. ALLOCATED( t_wall_h_1 ) ) &
1256	ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1257	IF ( .NOT. ALLOCATED( t_wall_h_2 ) ) &
1258	ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1259	IF ( .NOT. ALLOCATED( t_surf_window_h_1 ) ) &
1260	ALLOCATE ( t_surf_window_h_1(1:surf_usm_h%ns) )
1261	IF ( .NOT. ALLOCATED( t_surf_window_h_2 ) ) &
1262	ALLOCATE ( t_surf_window_h_2(1:surf_usm_h%ns) )
1263	IF ( .NOT. ALLOCATED( t_window_h_1 ) ) &
1264	ALLOCATE ( t_window_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1265	IF ( .NOT. ALLOCATED( t_window_h_2 ) ) &
1266	ALLOCATE ( t_window_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1267	IF ( .NOT. ALLOCATED( t_surf_green_h_1 ) ) &
1268	ALLOCATE ( t_surf_green_h_1(1:surf_usm_h%ns) )
1269	IF ( .NOT. ALLOCATED( t_surf_green_h_2 ) ) &
1270	ALLOCATE ( t_surf_green_h_2(1:surf_usm_h%ns) )
1271	IF ( .NOT. ALLOCATED( t_green_h_1 ) ) &
1272	ALLOCATE ( t_green_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1273	IF ( .NOT. ALLOCATED( t_green_h_2 ) ) &
1274	ALLOCATE ( t_green_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1275	IF ( .NOT. ALLOCATED( swc_h_1 ) ) &
1276	ALLOCATE ( swc_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1277	IF ( .NOT. ALLOCATED( swc_sat_h_1 ) ) &
1278	ALLOCATE ( swc_sat_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1279	IF ( .NOT. ALLOCATED( swc_res_h_1 ) ) &
1280	ALLOCATE ( swc_res_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1281	IF ( .NOT. ALLOCATED( swc_h_2 ) ) &
1282	ALLOCATE ( swc_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1283	IF ( .NOT. ALLOCATED( rootfr_h_1 ) ) &
1284	ALLOCATE ( rootfr_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1285	IF ( .NOT. ALLOCATED( wilt_h_1 ) ) &
1286	ALLOCATE ( wilt_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1287	IF ( .NOT. ALLOCATED( fc_h_1 ) ) &
1288	ALLOCATE ( fc_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1289
1290	IF ( .NOT. ALLOCATED( m_liq_usm_h_1%var_usm_1d ) ) &
1291	ALLOCATE ( m_liq_usm_h_1%var_usm_1d(1:surf_usm_h%ns) )
1292	IF ( .NOT. ALLOCATED( m_liq_usm_h_2%var_usm_1d ) ) &
1293	ALLOCATE ( m_liq_usm_h_2%var_usm_1d(1:surf_usm_h%ns) )
1294
1295	!
1296	!-- initial assignment of the pointers
1297	t_wall_h => t_wall_h_1; t_wall_h_p => t_wall_h_2
1298	t_window_h => t_window_h_1; t_window_h_p => t_window_h_2
1299	t_green_h => t_green_h_1; t_green_h_p => t_green_h_2
1300	t_surf_wall_h => t_surf_wall_h_1; t_surf_wall_h_p => t_surf_wall_h_2
1301	t_surf_window_h => t_surf_window_h_1; t_surf_window_h_p => t_surf_window_h_2
1302	t_surf_green_h => t_surf_green_h_1; t_surf_green_h_p => t_surf_green_h_2
1303	m_liq_usm_h => m_liq_usm_h_1; m_liq_usm_h_p => m_liq_usm_h_2
1304	swc_h => swc_h_1; swc_h_p => swc_h_2
1305	swc_sat_h => swc_sat_h_1
1306	swc_res_h => swc_res_h_1
1307	rootfr_h => rootfr_h_1
1308	wilt_h => wilt_h_1
1309	fc_h => fc_h_1
1310
1311	!
1312	!-- allocate wall and roof temperature arrays, for vertical walls if required
1313	!
1314	!-- Allocate if required. Note, in case of restarts, some of these arrays
1315	!-- might be already allocated.
1316	DO l = 0, 3
1317	IF ( .NOT. ALLOCATED( t_surf_wall_v_1(l)%t ) ) &
1318	ALLOCATE ( t_surf_wall_v_1(l)%t(1:surf_usm_v(l)%ns) )
1319	IF ( .NOT. ALLOCATED( t_surf_wall_v_2(l)%t ) ) &
1320	ALLOCATE ( t_surf_wall_v_2(l)%t(1:surf_usm_v(l)%ns) )
1321	IF ( .NOT. ALLOCATED( t_wall_v_1(l)%t ) ) &
1322	ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1323	IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) ) &
1324	ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1325	IF ( .NOT. ALLOCATED( t_surf_window_v_1(l)%t ) ) &
1326	ALLOCATE ( t_surf_window_v_1(l)%t(1:surf_usm_v(l)%ns) )
1327	IF ( .NOT. ALLOCATED( t_surf_window_v_2(l)%t ) ) &
1328	ALLOCATE ( t_surf_window_v_2(l)%t(1:surf_usm_v(l)%ns) )
1329	IF ( .NOT. ALLOCATED( t_window_v_1(l)%t ) ) &
1330	ALLOCATE ( t_window_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1331	IF ( .NOT. ALLOCATED( t_window_v_2(l)%t ) ) &
1332	ALLOCATE ( t_window_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1333	IF ( .NOT. ALLOCATED( t_surf_green_v_1(l)%t ) ) &
1334	ALLOCATE ( t_surf_green_v_1(l)%t(1:surf_usm_v(l)%ns) )
1335	IF ( .NOT. ALLOCATED( t_surf_green_v_2(l)%t ) ) &
1336	ALLOCATE ( t_surf_green_v_2(l)%t(1:surf_usm_v(l)%ns) )
1337	IF ( .NOT. ALLOCATED( t_green_v_1(l)%t ) ) &
1338	ALLOCATE ( t_green_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1339	IF ( .NOT. ALLOCATED( t_green_v_2(l)%t ) ) &
1340	ALLOCATE ( t_green_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1341	IF ( .NOT. ALLOCATED( m_liq_usm_v_1(l)%var_usm_1d ) ) &
1342	ALLOCATE ( m_liq_usm_v_1(l)%var_usm_1d(1:surf_usm_v(l)%ns) )
1343	IF ( .NOT. ALLOCATED( m_liq_usm_v_2(l)%var_usm_1d ) ) &
1344	ALLOCATE ( m_liq_usm_v_2(l)%var_usm_1d(1:surf_usm_v(l)%ns) )
1345	IF ( .NOT. ALLOCATED( swc_v_1(l)%t ) ) &
1346	ALLOCATE ( swc_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1347	IF ( .NOT. ALLOCATED( swc_v_2(l)%t ) ) &
1348	ALLOCATE ( swc_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1349	ENDDO
1350	!
1351	!-- initial assignment of the pointers
1352	t_wall_v => t_wall_v_1; t_wall_v_p => t_wall_v_2
1353	t_surf_wall_v => t_surf_wall_v_1; t_surf_wall_v_p => t_surf_wall_v_2
1354	t_window_v => t_window_v_1; t_window_v_p => t_window_v_2
1355	t_green_v => t_green_v_1; t_green_v_p => t_green_v_2
1356	t_surf_window_v => t_surf_window_v_1; t_surf_window_v_p => t_surf_window_v_2
1357	t_surf_green_v => t_surf_green_v_1; t_surf_green_v_p => t_surf_green_v_2
1358	m_liq_usm_v => m_liq_usm_v_1; m_liq_usm_v_p => m_liq_usm_v_2
1359	swc_v => swc_v_1; swc_v_p => swc_v_2
1360
1361	!
1362	!-- Allocate intermediate timestep arrays. For horizontal surfaces.
1363	ALLOCATE ( surf_usm_h%tt_surface_wall_m(1:surf_usm_h%ns) )
1364	ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1365	ALLOCATE ( surf_usm_h%tt_surface_window_m(1:surf_usm_h%ns) )
1366	ALLOCATE ( surf_usm_h%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1367	ALLOCATE ( surf_usm_h%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
1368	ALLOCATE ( surf_usm_h%tt_surface_green_m(1:surf_usm_h%ns) )
1369
1370	!
1371	!-- Allocate intermediate timestep arrays
1372	!-- Horizontal surfaces
1373	ALLOCATE ( tm_liq_usm_h_m%var_usm_1d(1:surf_usm_h%ns) )
1374	!
1375	!-- Horizontal surfaces
1376	DO l = 0, 3
1377	ALLOCATE ( tm_liq_usm_v_m(l)%var_usm_1d(1:surf_usm_v(l)%ns) )
1378	ENDDO
1379
1380	!
1381	!-- Set inital values for prognostic quantities
1382	IF ( ALLOCATED( surf_usm_h%tt_surface_wall_m ) ) surf_usm_h%tt_surface_wall_m = 0.0_wp
1383	IF ( ALLOCATED( surf_usm_h%tt_wall_m ) ) surf_usm_h%tt_wall_m = 0.0_wp
1384	IF ( ALLOCATED( surf_usm_h%tt_surface_window_m ) ) surf_usm_h%tt_surface_window_m = 0.0_wp
1385	IF ( ALLOCATED( surf_usm_h%tt_window_m ) ) surf_usm_h%tt_window_m = 0.0_wp
1386	IF ( ALLOCATED( surf_usm_h%tt_green_m ) ) surf_usm_h%tt_green_m = 0.0_wp
1387	IF ( ALLOCATED( surf_usm_h%tt_surface_green_m ) ) surf_usm_h%tt_surface_green_m = 0.0_wp
1388	!
1389	!-- Now, for vertical surfaces
1390	DO l = 0, 3
1391	ALLOCATE ( surf_usm_v(l)%tt_surface_wall_m(1:surf_usm_v(l)%ns) )
1392	ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1393	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_wall_m ) ) surf_usm_v(l)%tt_surface_wall_m = 0.0_wp
1394	IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m ) ) surf_usm_v(l)%tt_wall_m = 0.0_wp
1395	ALLOCATE ( surf_usm_v(l)%tt_surface_window_m(1:surf_usm_v(l)%ns) )
1396	ALLOCATE ( surf_usm_v(l)%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1397	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_window_m ) ) surf_usm_v(l)%tt_surface_window_m = 0.0_wp
1398	IF ( ALLOCATED( surf_usm_v(l)%tt_window_m ) ) surf_usm_v(l)%tt_window_m = 0.0_wp
1399	ALLOCATE ( surf_usm_v(l)%tt_surface_green_m(1:surf_usm_v(l)%ns) )
1400	IF ( ALLOCATED( surf_usm_v(l)%tt_surface_green_m ) ) surf_usm_v(l)%tt_surface_green_m = 0.0_wp
1401	ALLOCATE ( surf_usm_v(l)%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
1402	IF ( ALLOCATED( surf_usm_v(l)%tt_green_m ) ) surf_usm_v(l)%tt_green_m = 0.0_wp
1403	ENDDO
1404	!
1405	!-- allocate wall heat flux output array and set initial values. For horizontal surfaces
1406	! ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns) ) !can be removed
1407	ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
1408	ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
1409	ALLOCATE ( surf_usm_h%wghf_eb_window(1:surf_usm_h%ns) )
1410	ALLOCATE ( surf_usm_h%wghf_eb_green(1:surf_usm_h%ns) )
1411	ALLOCATE ( surf_usm_h%iwghf_eb(1:surf_usm_h%ns) )
1412	ALLOCATE ( surf_usm_h%iwghf_eb_window(1:surf_usm_h%ns) )
1413	IF ( ALLOCATED( surf_usm_h%wshf ) ) surf_usm_h%wshf = 0.0_wp
1414	IF ( ALLOCATED( surf_usm_h%wshf_eb ) ) surf_usm_h%wshf_eb = 0.0_wp
1415	IF ( ALLOCATED( surf_usm_h%wghf_eb ) ) surf_usm_h%wghf_eb = 0.0_wp
1416	IF ( ALLOCATED( surf_usm_h%wghf_eb_window ) ) surf_usm_h%wghf_eb_window = 0.0_wp
1417	IF ( ALLOCATED( surf_usm_h%wghf_eb_green ) ) surf_usm_h%wghf_eb_green = 0.0_wp
1418	IF ( ALLOCATED( surf_usm_h%iwghf_eb ) ) surf_usm_h%iwghf_eb = 0.0_wp
1419	IF ( ALLOCATED( surf_usm_h%iwghf_eb_window ) ) surf_usm_h%iwghf_eb_window = 0.0_wp
1420	!
1421	!-- Now, for vertical surfaces
1422	DO l = 0, 3
1423	! ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns) ) ! can be removed
1424	ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
1425	ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
1426	ALLOCATE ( surf_usm_v(l)%wghf_eb_window(1:surf_usm_v(l)%ns) )
1427	ALLOCATE ( surf_usm_v(l)%wghf_eb_green(1:surf_usm_v(l)%ns) )
1428	ALLOCATE ( surf_usm_v(l)%iwghf_eb(1:surf_usm_v(l)%ns) )
1429	ALLOCATE ( surf_usm_v(l)%iwghf_eb_window(1:surf_usm_v(l)%ns) )
1430	IF ( ALLOCATED( surf_usm_v(l)%wshf ) ) surf_usm_v(l)%wshf = 0.0_wp
1431	IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) ) surf_usm_v(l)%wshf_eb = 0.0_wp
1432	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) ) surf_usm_v(l)%wghf_eb = 0.0_wp
1433	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_window ) ) surf_usm_v(l)%wghf_eb_window = 0.0_wp
1434	IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_green ) ) surf_usm_v(l)%wghf_eb_green = 0.0_wp
1435	IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb ) ) surf_usm_v(l)%iwghf_eb = 0.0_wp
1436	IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) ) surf_usm_v(l)%iwghf_eb_window = 0.0_wp
1437	ENDDO
1438	!
1439	!-- Initialize building-surface properties, which are also required by other modules,
1440	!-- e.g. the indoor model.
1441	CALL usm_define_pars
1442
1443	IF ( debug_output ) CALL debug_message( 'usm_init_arrays', 'end' )
1444
1445	END SUBROUTINE usm_init_arrays
1446
1447
1448	!------------------------------------------------------------------------------!
1449	! Description:
1450	! ------------
1451	!> Sum up and time-average urban surface output quantities as well as allocate
1452	!> the array necessary for storing the average.
1453	!------------------------------------------------------------------------------!
1454	SUBROUTINE usm_3d_data_averaging( mode, variable )
1455
1456	IMPLICIT NONE
1457
1458	CHARACTER(LEN=*), INTENT(IN) :: mode
1459	CHARACTER(LEN=*), INTENT(IN) :: variable
1460
1461	INTEGER(iwp) :: i, j, k, l, m, ids, idsint, iwl, istat !< runnin indices
1462	CHARACTER(LEN=varnamelength) :: var !< trimmed variable
1463	INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions
1464	CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
1465	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
1466
1467	IF ( variable(1:4) == 'usm_' ) THEN ! is such a check really rquired?
1468
1469	!
1470	!-- find the real name of the variable
1471	ids = -1
1472	l = -1
1473	var = TRIM(variable)
1474	DO i = 0, nd-1
1475	k = len(TRIM(var))
1476	j = len(TRIM(dirname(i)))
1477	IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN
1478	ids = i
1479	idsint = dirint(ids)
1480	var = var(:k-j)
1481	EXIT
1482	ENDIF
1483	ENDDO
1484	l = idsint - 2 ! horisontal direction index - terible hack !
1485	IF ( l < 0 .OR. l > 3 ) THEN
1486	l = -1
1487	END IF
1488	IF ( ids == -1 ) THEN
1489	var = TRIM(variable)
1490	ENDIF
1491	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
1492	!
1493	!-- wall layers
1494	READ(var(12:12), '(I1)', iostat=istat ) iwl
1495	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1496	var = var(1:10)
1497	ELSE
1498	!
1499	!-- wrong wall layer index
1500	RETURN
1501	ENDIF
1502	ENDIF
1503	IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN
1504	!
1505	!-- wall layers
1506	READ(var(14:14), '(I1)', iostat=istat ) iwl
1507	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1508	var = var(1:12)
1509	ELSE
1510	!
1511	!-- wrong window layer index
1512	RETURN
1513	ENDIF
1514	ENDIF
1515	IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN
1516	!
1517	!-- wall layers
1518	READ(var(13:13), '(I1)', iostat=istat ) iwl
1519	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1520	var = var(1:11)
1521	ELSE
1522	!
1523	!-- wrong green layer index
1524	RETURN
1525	ENDIF
1526	ENDIF
1527	IF ( var(1:8) == 'usm_swc_' .AND. len(TRIM(var)) >= 9 ) THEN
1528	!
1529	!-- swc layers
1530	READ(var(9:9), '(I1)', iostat=istat ) iwl
1531	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
1532	var = var(1:7)
1533	ELSE
1534	!
1535	!-- wrong swc layer index
1536	RETURN
1537	ENDIF
1538	ENDIF
1539
1540	IF ( mode == 'allocate' ) THEN
1541
1542	SELECT CASE ( TRIM( var ) )
1543
1544	CASE ( 'usm_wshf' )
1545	!
1546	!-- array of sensible heat flux from surfaces
1547	!-- land surfaces
1548	IF ( l == -1 ) THEN
1549	IF ( .NOT. ALLOCATED(surf_usm_h%wshf_eb_av) ) THEN
1550	ALLOCATE ( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1551	surf_usm_h%wshf_eb_av = 0.0_wp
1552	ENDIF
1553	ELSE
1554	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wshf_eb_av) ) THEN
1555	ALLOCATE ( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1556	surf_usm_v(l)%wshf_eb_av = 0.0_wp
1557	ENDIF
1558	ENDIF
1559
1560	CASE ( 'usm_qsws' )
1561	!
1562	!-- array of latent heat flux from surfaces
1563	!-- land surfaces
1564	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%qsws_av) ) THEN
1565	ALLOCATE ( surf_usm_h%qsws_av(1:surf_usm_h%ns) )
1566	surf_usm_h%qsws_av = 0.0_wp
1567	ELSE
1568	IF ( .NOT. ALLOCATED(surf_usm_v(l)%qsws_av) ) THEN
1569	ALLOCATE ( surf_usm_v(l)%qsws_av(1:surf_usm_v(l)%ns) )
1570	surf_usm_v(l)%qsws_av = 0.0_wp
1571	ENDIF
1572	ENDIF
1573
1574	CASE ( 'usm_qsws_veg' )
1575	!
1576	!-- array of latent heat flux from vegetation surfaces
1577	!-- land surfaces
1578	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%qsws_veg_av) ) THEN
1579	ALLOCATE ( surf_usm_h%qsws_veg_av(1:surf_usm_h%ns) )
1580	surf_usm_h%qsws_veg_av = 0.0_wp
1581	ELSE
1582	IF ( .NOT. ALLOCATED(surf_usm_v(l)%qsws_veg_av) ) THEN
1583	ALLOCATE ( surf_usm_v(l)%qsws_veg_av(1:surf_usm_v(l)%ns) )
1584	surf_usm_v(l)%qsws_veg_av = 0.0_wp
1585	ENDIF
1586	ENDIF
1587
1588	CASE ( 'usm_qsws_liq' )
1589	!
1590	!-- array of latent heat flux from surfaces with liquid
1591	!-- land surfaces
1592	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%qsws_liq_av) ) THEN
1593	ALLOCATE ( surf_usm_h%qsws_liq_av(1:surf_usm_h%ns) )
1594	surf_usm_h%qsws_liq_av = 0.0_wp
1595	ELSE
1596	IF ( .NOT. ALLOCATED(surf_usm_v(l)%qsws_liq_av) ) THEN
1597	ALLOCATE ( surf_usm_v(l)%qsws_liq_av(1:surf_usm_v(l)%ns) )
1598	surf_usm_v(l)%qsws_liq_av = 0.0_wp
1599	ENDIF
1600	ENDIF
1601	!
1602	!-- Please note, the following output quantities belongs to the
1603	!-- individual tile fractions - ground heat flux at wall-, window-,
1604	!-- and green fraction. Aggregated ground-heat flux is treated
1605	!-- accordingly in average_3d_data, sum_up_3d_data, etc..
1606	CASE ( 'usm_wghf' )
1607	!
1608	!-- array of heat flux from ground (wall, roof, land)
1609	IF ( l == -1 ) THEN
1610	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_av) ) THEN
1611	ALLOCATE ( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1612	surf_usm_h%wghf_eb_av = 0.0_wp
1613	ENDIF
1614	ELSE
1615	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_av) ) THEN
1616	ALLOCATE ( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1617	surf_usm_v(l)%wghf_eb_av = 0.0_wp
1618	ENDIF
1619	ENDIF
1620
1621	CASE ( 'usm_wghf_window' )
1622	!
1623	!-- array of heat flux from window ground (wall, roof, land)
1624	IF ( l == -1 ) THEN
1625	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_window_av) ) THEN
1626	ALLOCATE ( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
1627	surf_usm_h%wghf_eb_window_av = 0.0_wp
1628	ENDIF
1629	ELSE
1630	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) ) THEN
1631	ALLOCATE ( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
1632	surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
1633	ENDIF
1634	ENDIF
1635
1636	CASE ( 'usm_wghf_green' )
1637	!
1638	!-- array of heat flux from green ground (wall, roof, land)
1639	IF ( l == -1 ) THEN
1640	IF ( .NOT. ALLOCATED(surf_usm_h%wghf_eb_green_av) ) THEN
1641	ALLOCATE ( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
1642	surf_usm_h%wghf_eb_green_av = 0.0_wp
1643	ENDIF
1644	ELSE
1645	IF ( .NOT. ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) ) THEN
1646	ALLOCATE ( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
1647	surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
1648	ENDIF
1649	ENDIF
1650
1651	CASE ( 'usm_iwghf' )
1652	!
1653	!-- array of heat flux from indoor ground (wall, roof, land)
1654	IF ( l == -1 ) THEN
1655	IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_av) ) THEN
1656	ALLOCATE ( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
1657	surf_usm_h%iwghf_eb_av = 0.0_wp
1658	ENDIF
1659	ELSE
1660	IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_av) ) THEN
1661	ALLOCATE ( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
1662	surf_usm_v(l)%iwghf_eb_av = 0.0_wp
1663	ENDIF
1664	ENDIF
1665
1666	CASE ( 'usm_iwghf_window' )
1667	!
1668	!-- array of heat flux from indoor window ground (wall, roof, land)
1669	IF ( l == -1 ) THEN
1670	IF ( .NOT. ALLOCATED(surf_usm_h%iwghf_eb_window_av) ) THEN
1671	ALLOCATE ( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
1672	surf_usm_h%iwghf_eb_window_av = 0.0_wp
1673	ENDIF
1674	ELSE
1675	IF ( .NOT. ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) ) THEN
1676	ALLOCATE ( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
1677	surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
1678	ENDIF
1679	ENDIF
1680
1681	CASE ( 'usm_t_surf_wall' )
1682	!
1683	!-- surface temperature for surfaces
1684	IF ( l == -1 ) THEN
1685	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_wall_av) ) THEN
1686	ALLOCATE ( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) )
1687	surf_usm_h%t_surf_wall_av = 0.0_wp
1688	ENDIF
1689	ELSE
1690	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_wall_av) ) THEN
1691	ALLOCATE ( surf_usm_v(l)%t_surf_wall_av(1:surf_usm_v(l)%ns) )
1692	surf_usm_v(l)%t_surf_wall_av = 0.0_wp
1693	ENDIF
1694	ENDIF
1695
1696	CASE ( 'usm_t_surf_window' )
1697	!
1698	!-- surface temperature for window surfaces
1699	IF ( l == -1 ) THEN
1700	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_window_av) ) THEN
1701	ALLOCATE ( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
1702	surf_usm_h%t_surf_window_av = 0.0_wp
1703	ENDIF
1704	ELSE
1705	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_window_av) ) THEN
1706	ALLOCATE ( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
1707	surf_usm_v(l)%t_surf_window_av = 0.0_wp
1708	ENDIF
1709	ENDIF
1710
1711	CASE ( 'usm_t_surf_green' )
1712	!
1713	!-- surface temperature for green surfaces
1714	IF ( l == -1 ) THEN
1715	IF ( .NOT. ALLOCATED(surf_usm_h%t_surf_green_av) ) THEN
1716	ALLOCATE ( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
1717	surf_usm_h%t_surf_green_av = 0.0_wp
1718	ENDIF
1719	ELSE
1720	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_surf_green_av) ) THEN
1721	ALLOCATE ( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
1722	surf_usm_v(l)%t_surf_green_av = 0.0_wp
1723	ENDIF
1724	ENDIF
1725
1726	CASE ( 'usm_theta_10cm' )
1727	!
1728	!-- near surface (10cm) temperature for whole surfaces
1729	IF ( l == -1 ) THEN
1730	IF ( .NOT. ALLOCATED(surf_usm_h%pt_10cm_av) ) THEN
1731	ALLOCATE ( surf_usm_h%pt_10cm_av(1:surf_usm_h%ns) )
1732	surf_usm_h%pt_10cm_av = 0.0_wp
1733	ENDIF
1734	ELSE
1735	IF ( .NOT. ALLOCATED(surf_usm_v(l)%pt_10cm_av) ) THEN
1736	ALLOCATE ( surf_usm_v(l)%pt_10cm_av(1:surf_usm_v(l)%ns) )
1737	surf_usm_v(l)%pt_10cm_av = 0.0_wp
1738	ENDIF
1739	ENDIF
1740
1741	CASE ( 'usm_t_wall' )
1742	!
1743	!-- wall temperature for iwl layer of walls and land
1744	IF ( l == -1 ) THEN
1745	IF ( .NOT. ALLOCATED(surf_usm_h%t_wall_av) ) THEN
1746	ALLOCATE ( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1747	surf_usm_h%t_wall_av = 0.0_wp
1748	ENDIF
1749	ELSE
1750	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_wall_av) ) THEN
1751	ALLOCATE ( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1752	surf_usm_v(l)%t_wall_av = 0.0_wp
1753	ENDIF
1754	ENDIF
1755
1756	CASE ( 'usm_t_window' )
1757	!
1758	!-- window temperature for iwl layer of walls and land
1759	IF ( l == -1 ) THEN
1760	IF ( .NOT. ALLOCATED(surf_usm_h%t_window_av) ) THEN
1761	ALLOCATE ( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1762	surf_usm_h%t_window_av = 0.0_wp
1763	ENDIF
1764	ELSE
1765	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_window_av) ) THEN
1766	ALLOCATE ( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1767	surf_usm_v(l)%t_window_av = 0.0_wp
1768	ENDIF
1769	ENDIF
1770
1771	CASE ( 'usm_t_green' )
1772	!
1773	!-- green temperature for iwl layer of walls and land
1774	IF ( l == -1 ) THEN
1775	IF ( .NOT. ALLOCATED(surf_usm_h%t_green_av) ) THEN
1776	ALLOCATE ( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1777	surf_usm_h%t_green_av = 0.0_wp
1778	ENDIF
1779	ELSE
1780	IF ( .NOT. ALLOCATED(surf_usm_v(l)%t_green_av) ) THEN
1781	ALLOCATE ( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1782	surf_usm_v(l)%t_green_av = 0.0_wp
1783	ENDIF
1784	ENDIF
1785	CASE ( 'usm_swc' )
1786	!
1787	!-- soil water content for iwl layer of walls and land
1788	IF ( l == -1 .AND. .NOT. ALLOCATED(surf_usm_h%swc_av) ) THEN
1789	ALLOCATE ( surf_usm_h%swc_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1790	surf_usm_h%swc_av = 0.0_wp
1791	ELSE
1792	IF ( .NOT. ALLOCATED(surf_usm_v(l)%swc_av) ) THEN
1793	ALLOCATE ( surf_usm_v(l)%swc_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1794	surf_usm_v(l)%swc_av = 0.0_wp
1795	ENDIF
1796	ENDIF
1797
1798	CASE DEFAULT
1799	CONTINUE
1800
1801	END SELECT
1802
1803	ELSEIF ( mode == 'sum' ) 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	surf_usm_h%wshf_eb(m)
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	surf_usm_v(l)%wshf_eb(m)
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	surf_usm_h%qsws(m) * l_v
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	surf_usm_v(l)%qsws(m) * l_v
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	surf_usm_h%qsws_veg(m)
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	surf_usm_v(l)%qsws_veg(m)
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	surf_usm_h%qsws_liq(m)
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	surf_usm_v(l)%qsws_liq(m)
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	surf_usm_h%wghf_eb(m)
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	surf_usm_v(l)%wghf_eb(m)
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	surf_usm_h%wghf_eb_window(m)
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	surf_usm_v(l)%wghf_eb_window(m)
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	surf_usm_h%wghf_eb_green(m)
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	surf_usm_v(l)%wghf_eb_green(m)
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	surf_usm_h%iwghf_eb(m)
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	surf_usm_v(l)%iwghf_eb(m)
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	surf_usm_h%iwghf_eb_window(m)
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	surf_usm_v(l)%iwghf_eb_window(m)
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	t_surf_wall_h(m)
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	t_surf_wall_v(l)%t(m)
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	t_surf_window_h(m)
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	t_surf_window_v(l)%t(m)
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	t_surf_green_h(m)
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	t_surf_green_v(l)%t(m)
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	surf_usm_h%pt_10cm(m)
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	surf_usm_v(l)%pt_10cm(m)
2025	ENDDO
2026	ENDIF
2027
2028	CASE ( 'usm_t_wall' )
2029	!
2030	!-- wall temperature for iwl layer of walls and land
2031	IF ( l == -1 ) THEN
2032	DO m = 1, surf_usm_h%ns
2033	surf_usm_h%t_wall_av(iwl,m) = &
2034	surf_usm_h%t_wall_av(iwl,m) + &
2035	t_wall_h(iwl,m)
2036	ENDDO
2037	ELSE
2038	DO m = 1, surf_usm_v(l)%ns
2039	surf_usm_v(l)%t_wall_av(iwl,m) = &
2040	surf_usm_v(l)%t_wall_av(iwl,m) + &
2041	t_wall_v(l)%t(iwl,m)
2042	ENDDO
2043	ENDIF
2044
2045	CASE ( 'usm_t_window' )
2046	!
2047	!-- window temperature for iwl layer of walls and land
2048	IF ( l == -1 ) THEN
2049	DO m = 1, surf_usm_h%ns
2050	surf_usm_h%t_window_av(iwl,m) = &
2051	surf_usm_h%t_window_av(iwl,m) + &
2052	t_window_h(iwl,m)
2053	ENDDO
2054	ELSE
2055	DO m = 1, surf_usm_v(l)%ns
2056	surf_usm_v(l)%t_window_av(iwl,m) = &
2057	surf_usm_v(l)%t_window_av(iwl,m) + &
2058	t_window_v(l)%t(iwl,m)
2059	ENDDO
2060	ENDIF
2061
2062	CASE ( 'usm_t_green' )
2063	!
2064	!-- green temperature for iwl layer of walls and land
2065	IF ( l == -1 ) THEN
2066	DO m = 1, surf_usm_h%ns
2067	surf_usm_h%t_green_av(iwl,m) = &
2068	surf_usm_h%t_green_av(iwl,m) + &
2069	t_green_h(iwl,m)
2070	ENDDO
2071	ELSE
2072	DO m = 1, surf_usm_v(l)%ns
2073	surf_usm_v(l)%t_green_av(iwl,m) = &
2074	surf_usm_v(l)%t_green_av(iwl,m) + &
2075	t_green_v(l)%t(iwl,m)
2076	ENDDO
2077	ENDIF
2078
2079	CASE ( 'usm_swc' )
2080	!
2081	!-- soil water content for iwl layer of walls and land
2082	IF ( l == -1 ) THEN
2083	DO m = 1, surf_usm_h%ns
2084	surf_usm_h%swc_av(iwl,m) = &
2085	surf_usm_h%swc_av(iwl,m) + &
2086	swc_h(iwl,m)
2087	ENDDO
2088	ELSE
2089	DO m = 1, surf_usm_v(l)%ns
2090	surf_usm_v(l)%swc_av(iwl,m) = &
2091	surf_usm_v(l)%swc_av(iwl,m) + &
2092	swc_v(l)%t(iwl,m)
2093	ENDDO
2094	ENDIF
2095
2096	CASE DEFAULT
2097	CONTINUE
2098
2099	END SELECT
2100
2101	ELSEIF ( mode == 'average' ) THEN
2102
2103	SELECT CASE ( TRIM( var ) )
2104
2105	CASE ( 'usm_wshf' )
2106	!
2107	!-- array of sensible heat flux from surfaces (land, roof, wall)
2108	IF ( l == -1 ) THEN
2109	DO m = 1, surf_usm_h%ns
2110	surf_usm_h%wshf_eb_av(m) = &
2111	surf_usm_h%wshf_eb_av(m) / &
2112	REAL( average_count_3d, kind=wp )
2113	ENDDO
2114	ELSE
2115	DO m = 1, surf_usm_v(l)%ns
2116	surf_usm_v(l)%wshf_eb_av(m) = &
2117	surf_usm_v(l)%wshf_eb_av(m) / &
2118	REAL( average_count_3d, kind=wp )
2119	ENDDO
2120	ENDIF
2121
2122	CASE ( 'usm_qsws' )
2123	!
2124	!-- array of latent heat flux from surfaces (land, roof, wall)
2125	IF ( l == -1 ) THEN
2126	DO m = 1, surf_usm_h%ns
2127	surf_usm_h%qsws_av(m) = &
2128	surf_usm_h%qsws_av(m) / &
2129	REAL( average_count_3d, kind=wp )
2130	ENDDO
2131	ELSE
2132	DO m = 1, surf_usm_v(l)%ns
2133	surf_usm_v(l)%qsws_av(m) = &
2134	surf_usm_v(l)%qsws_av(m) / &
2135	REAL( average_count_3d, kind=wp )
2136	ENDDO
2137	ENDIF
2138
2139	CASE ( 'usm_qsws_veg' )
2140	!
2141	!-- array of latent heat flux from vegetation surfaces (land, roof, wall)
2142	IF ( l == -1 ) THEN
2143	DO m = 1, surf_usm_h%ns
2144	surf_usm_h%qsws_veg_av(m) = &
2145	surf_usm_h%qsws_veg_av(m) / &
2146	REAL( average_count_3d, kind=wp )
2147	ENDDO
2148	ELSE
2149	DO m = 1, surf_usm_v(l)%ns
2150	surf_usm_v(l)%qsws_veg_av(m) = &
2151	surf_usm_v(l)%qsws_veg_av(m) / &
2152	REAL( average_count_3d, kind=wp )
2153	ENDDO
2154	ENDIF
2155
2156	CASE ( 'usm_qsws_liq' )
2157	!
2158	!-- array of latent heat flux from surfaces with liquid (land, roof, wall)
2159	IF ( l == -1 ) THEN
2160	DO m = 1, surf_usm_h%ns
2161	surf_usm_h%qsws_liq_av(m) = &
2162	surf_usm_h%qsws_liq_av(m) / &
2163	REAL( average_count_3d, kind=wp )
2164	ENDDO
2165	ELSE
2166	DO m = 1, surf_usm_v(l)%ns
2167	surf_usm_v(l)%qsws_liq_av(m) = &
2168	surf_usm_v(l)%qsws_liq_av(m) / &
2169	REAL( average_count_3d, kind=wp )
2170	ENDDO
2171	ENDIF
2172
2173	CASE ( 'usm_wghf' )
2174	!
2175	!-- array of heat flux from ground (wall, roof, land)
2176	IF ( l == -1 ) THEN
2177	DO m = 1, surf_usm_h%ns
2178	surf_usm_h%wghf_eb_av(m) = &
2179	surf_usm_h%wghf_eb_av(m) / &
2180	REAL( average_count_3d, kind=wp )
2181	ENDDO
2182	ELSE
2183	DO m = 1, surf_usm_v(l)%ns
2184	surf_usm_v(l)%wghf_eb_av(m) = &
2185	surf_usm_v(l)%wghf_eb_av(m) / &
2186	REAL( average_count_3d, kind=wp )
2187	ENDDO
2188	ENDIF
2189
2190	CASE ( 'usm_wghf_window' )
2191	!
2192	!-- array of heat flux from window ground (wall, roof, land)
2193	IF ( l == -1 ) THEN
2194	DO m = 1, surf_usm_h%ns
2195	surf_usm_h%wghf_eb_window_av(m) = &
2196	surf_usm_h%wghf_eb_window_av(m) / &
2197	REAL( average_count_3d, kind=wp )
2198	ENDDO
2199	ELSE
2200	DO m = 1, surf_usm_v(l)%ns
2201	surf_usm_v(l)%wghf_eb_window_av(m) = &
2202	surf_usm_v(l)%wghf_eb_window_av(m) / &
2203	REAL( average_count_3d, kind=wp )
2204	ENDDO
2205	ENDIF
2206
2207	CASE ( 'usm_wghf_green' )
2208	!
2209	!-- array of heat flux from green ground (wall, roof, land)
2210	IF ( l == -1 ) THEN
2211	DO m = 1, surf_usm_h%ns
2212	surf_usm_h%wghf_eb_green_av(m) = &
2213	surf_usm_h%wghf_eb_green_av(m) / &
2214	REAL( average_count_3d, kind=wp )
2215	ENDDO
2216	ELSE
2217	DO m = 1, surf_usm_v(l)%ns
2218	surf_usm_v(l)%wghf_eb_green_av(m) = &
2219	surf_usm_v(l)%wghf_eb_green_av(m) / &
2220	REAL( average_count_3d, kind=wp )
2221	ENDDO
2222	ENDIF
2223
2224	CASE ( 'usm_iwghf' )
2225	!
2226	!-- array of heat flux from indoor ground (wall, roof, land)
2227	IF ( l == -1 ) THEN
2228	DO m = 1, surf_usm_h%ns
2229	surf_usm_h%iwghf_eb_av(m) = &
2230	surf_usm_h%iwghf_eb_av(m) / &
2231	REAL( average_count_3d, kind=wp )
2232	ENDDO
2233	ELSE
2234	DO m = 1, surf_usm_v(l)%ns
2235	surf_usm_v(l)%iwghf_eb_av(m) = &
2236	surf_usm_v(l)%iwghf_eb_av(m) / &
2237	REAL( average_count_3d, kind=wp )
2238	ENDDO
2239	ENDIF
2240
2241	CASE ( 'usm_iwghf_window' )
2242	!
2243	!-- array of heat flux from indoor window ground (wall, roof, land)
2244	IF ( l == -1 ) THEN
2245	DO m = 1, surf_usm_h%ns
2246	surf_usm_h%iwghf_eb_window_av(m) = &
2247	surf_usm_h%iwghf_eb_window_av(m) / &
2248	REAL( average_count_3d, kind=wp )
2249	ENDDO
2250	ELSE
2251	DO m = 1, surf_usm_v(l)%ns
2252	surf_usm_v(l)%iwghf_eb_window_av(m) = &
2253	surf_usm_v(l)%iwghf_eb_window_av(m) / &
2254	REAL( average_count_3d, kind=wp )
2255	ENDDO
2256	ENDIF
2257
2258	CASE ( 'usm_t_surf_wall' )
2259	!
2260	!-- surface temperature for surfaces
2261	IF ( l == -1 ) THEN
2262	DO m = 1, surf_usm_h%ns
2263	surf_usm_h%t_surf_wall_av(m) = &
2264	surf_usm_h%t_surf_wall_av(m) / &
2265	REAL( average_count_3d, kind=wp )
2266	ENDDO
2267	ELSE
2268	DO m = 1, surf_usm_v(l)%ns
2269	surf_usm_v(l)%t_surf_wall_av(m) = &
2270	surf_usm_v(l)%t_surf_wall_av(m) / &
2271	REAL( average_count_3d, kind=wp )
2272	ENDDO
2273	ENDIF
2274
2275	CASE ( 'usm_t_surf_window' )
2276	!
2277	!-- surface temperature for window surfaces
2278	IF ( l == -1 ) THEN
2279	DO m = 1, surf_usm_h%ns
2280	surf_usm_h%t_surf_window_av(m) = &
2281	surf_usm_h%t_surf_window_av(m) / &
2282	REAL( average_count_3d, kind=wp )
2283	ENDDO
2284	ELSE
2285	DO m = 1, surf_usm_v(l)%ns
2286	surf_usm_v(l)%t_surf_window_av(m) = &
2287	surf_usm_v(l)%t_surf_window_av(m) / &
2288	REAL( average_count_3d, kind=wp )
2289	ENDDO
2290	ENDIF
2291
2292	CASE ( 'usm_t_surf_green' )
2293	!
2294	!-- surface temperature for green surfaces
2295	IF ( l == -1 ) THEN
2296	DO m = 1, surf_usm_h%ns
2297	surf_usm_h%t_surf_green_av(m) = &
2298	surf_usm_h%t_surf_green_av(m) / &
2299	REAL( average_count_3d, kind=wp )
2300	ENDDO
2301	ELSE
2302	DO m = 1, surf_usm_v(l)%ns
2303	surf_usm_v(l)%t_surf_green_av(m) = &
2304	surf_usm_v(l)%t_surf_green_av(m) / &
2305	REAL( average_count_3d, kind=wp )
2306	ENDDO
2307	ENDIF
2308
2309	CASE ( 'usm_theta_10cm' )
2310	!
2311	!-- near surface temperature for whole surfaces
2312	IF ( l == -1 ) THEN
2313	DO m = 1, surf_usm_h%ns
2314	surf_usm_h%pt_10cm_av(m) = &
2315	surf_usm_h%pt_10cm_av(m) / &
2316	REAL( average_count_3d, kind=wp )
2317	ENDDO
2318	ELSE
2319	DO m = 1, surf_usm_v(l)%ns
2320	surf_usm_v(l)%pt_10cm_av(m) = &
2321	surf_usm_v(l)%pt_10cm_av(m) / &
2322	REAL( average_count_3d, kind=wp )
2323	ENDDO
2324	ENDIF
2325
2326
2327	CASE ( 'usm_t_wall' )
2328	!
2329	!-- wall temperature for iwl layer of walls and land
2330	IF ( l == -1 ) THEN
2331	DO m = 1, surf_usm_h%ns
2332	surf_usm_h%t_wall_av(iwl,m) = &
2333	surf_usm_h%t_wall_av(iwl,m) / &
2334	REAL( average_count_3d, kind=wp )
2335	ENDDO
2336	ELSE
2337	DO m = 1, surf_usm_v(l)%ns
2338	surf_usm_v(l)%t_wall_av(iwl,m) = &
2339	surf_usm_v(l)%t_wall_av(iwl,m) / &
2340	REAL( average_count_3d, kind=wp )
2341	ENDDO
2342	ENDIF
2343
2344	CASE ( 'usm_t_window' )
2345	!
2346	!-- window temperature for iwl layer of walls and land
2347	IF ( l == -1 ) THEN
2348	DO m = 1, surf_usm_h%ns
2349	surf_usm_h%t_window_av(iwl,m) = &
2350	surf_usm_h%t_window_av(iwl,m) / &
2351	REAL( average_count_3d, kind=wp )
2352	ENDDO
2353	ELSE
2354	DO m = 1, surf_usm_v(l)%ns
2355	surf_usm_v(l)%t_window_av(iwl,m) = &
2356	surf_usm_v(l)%t_window_av(iwl,m) / &
2357	REAL( average_count_3d, kind=wp )
2358	ENDDO
2359	ENDIF
2360
2361	CASE ( 'usm_t_green' )
2362	!
2363	!-- green temperature for iwl layer of walls and land
2364	IF ( l == -1 ) THEN
2365	DO m = 1, surf_usm_h%ns
2366	surf_usm_h%t_green_av(iwl,m) = &
2367	surf_usm_h%t_green_av(iwl,m) / &
2368	REAL( average_count_3d, kind=wp )
2369	ENDDO
2370	ELSE
2371	DO m = 1, surf_usm_v(l)%ns
2372	surf_usm_v(l)%t_green_av(iwl,m) = &
2373	surf_usm_v(l)%t_green_av(iwl,m) / &
2374	REAL( average_count_3d, kind=wp )
2375	ENDDO
2376	ENDIF
2377
2378	CASE ( 'usm_swc' )
2379	!
2380	!-- soil water content for iwl layer of walls and land
2381	IF ( l == -1 ) THEN
2382	DO m = 1, surf_usm_h%ns
2383	surf_usm_h%swc_av(iwl,m) = &
2384	surf_usm_h%swc_av(iwl,m) / &
2385	REAL( average_count_3d, kind=wp )
2386	ENDDO
2387	ELSE
2388	DO m = 1, surf_usm_v(l)%ns
2389	surf_usm_v(l)%swc_av(iwl,m) = &
2390	surf_usm_v(l)%swc_av(iwl,m) / &
2391	REAL( average_count_3d, kind=wp )
2392	ENDDO
2393	ENDIF
2394
2395
2396	END SELECT
2397
2398	ENDIF
2399
2400	ENDIF
2401
2402	END SUBROUTINE usm_3d_data_averaging
2403
2404
2405
2406	!------------------------------------------------------------------------------!
2407	! Description:
2408	! ------------
2409	!> Set internal Neumann boundary condition at outer soil grid points
2410	!> for temperature and humidity.
2411	!------------------------------------------------------------------------------!
2412	SUBROUTINE usm_boundary_condition
2413
2414	IMPLICIT NONE
2415
2416	INTEGER(iwp) :: i !< grid index x-direction
2417	INTEGER(iwp) :: ioff !< offset index x-direction indicating location of soil grid point
2418	INTEGER(iwp) :: j !< grid index y-direction
2419	INTEGER(iwp) :: joff !< offset index x-direction indicating location of soil grid point
2420	INTEGER(iwp) :: k !< grid index z-direction
2421	INTEGER(iwp) :: koff !< offset index x-direction indicating location of soil grid point
2422	INTEGER(iwp) :: l !< running index surface-orientation
2423	INTEGER(iwp) :: m !< running index surface elements
2424
2425	koff = surf_usm_h%koff
2426	DO m = 1, surf_usm_h%ns
2427	i = surf_usm_h%i(m)
2428	j = surf_usm_h%j(m)
2429	k = surf_usm_h%k(m)
2430	pt(k+koff,j,i) = pt(k,j,i)
2431	ENDDO
2432
2433	DO l = 0, 3
2434	ioff = surf_usm_v(l)%ioff
2435	joff = surf_usm_v(l)%joff
2436	DO m = 1, surf_usm_v(l)%ns
2437	i = surf_usm_v(l)%i(m)
2438	j = surf_usm_v(l)%j(m)
2439	k = surf_usm_v(l)%k(m)
2440	pt(k,j+joff,i+ioff) = pt(k,j,i)
2441	ENDDO
2442	ENDDO
2443
2444	END SUBROUTINE usm_boundary_condition
2445
2446
2447	!------------------------------------------------------------------------------!
2448	!
2449	! Description:
2450	! ------------
2451	!> Subroutine checks variables and assigns units.
2452	!> It is called out from subroutine check_parameters.
2453	!------------------------------------------------------------------------------!
2454	SUBROUTINE usm_check_data_output( variable, unit )
2455
2456	IMPLICIT NONE
2457
2458	CHARACTER(LEN=*),INTENT(IN) :: variable !<
2459	CHARACTER(LEN=*),INTENT(OUT) :: unit !<
2460
2461	INTEGER(iwp) :: i,j,l !< index
2462	CHARACTER(LEN=2) :: ls
2463	CHARACTER(LEN=varnamelength) :: var !< TRIM(variable)
2464	INTEGER(iwp), PARAMETER :: nl1 = 15 !< number of directional usm variables
2465	CHARACTER(LEN=varnamelength), DIMENSION(nl1) :: varlist1 = & !< list of directional usm variables
2466	(/'usm_wshf ', &
2467	'usm_wghf ', &
2468	'usm_wghf_window ', &
2469	'usm_wghf_green ', &
2470	'usm_iwghf ', &
2471	'usm_iwghf_window ', &
2472	'usm_surfz ', &
2473	'usm_surfwintrans ', &
2474	'usm_surfcat ', &
2475	'usm_t_surf_wall ', &
2476	'usm_t_surf_window ', &
2477	'usm_t_surf_green ', &
2478	'usm_t_green ', &
2479	'usm_qsws ', &
2480	'usm_theta_10cm '/)
2481
2482	INTEGER(iwp), PARAMETER :: nl2 = 3 !< number of directional layer usm variables
2483	CHARACTER(LEN=varnamelength), DIMENSION(nl2) :: varlist2 = & !< list of directional layer usm variables
2484	(/'usm_t_wall ', &
2485	'usm_t_window ', &
2486	'usm_t_green '/)
2487
2488	INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions
2489	CHARACTER(LEN=6), DIMENSION(nd), PARAMETER :: dirname = & !< direction names
2490	(/'_roof ','_south','_north','_west ','_east '/)
2491	LOGICAL :: lfound !< flag if the variable is found
2492
2493
2494	lfound = .FALSE.
2495
2496	var = TRIM(variable)
2497
2498	!
2499	!-- check if variable exists
2500	!-- directional variables
2501	DO i = 1, nl1
2502	DO j = 1, nd
2503	IF ( TRIM(var) == TRIM(varlist1(i))//TRIM(dirname(j)) ) THEN
2504	lfound = .TRUE.
2505	EXIT
2506	ENDIF
2507	IF ( lfound ) EXIT
2508	ENDDO
2509	ENDDO
2510	IF ( lfound ) GOTO 10
2511	!
2512	!-- directional layer variables
2513	DO i = 1, nl2
2514	DO j = 1, nd
2515	DO l = nzb_wall, nzt_wall
2516	WRITE(ls,'(A1,I1)') '_',l
2517	IF ( TRIM(var) == TRIM(varlist2(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
2518	lfound = .TRUE.
2519	EXIT
2520	ENDIF
2521	ENDDO
2522	IF ( lfound ) EXIT
2523	ENDDO
2524	ENDDO
2525	IF ( .NOT. lfound ) THEN
2526	unit = 'illegal'
2527	RETURN
2528	ENDIF
2529	10 CONTINUE
2530
2531	IF ( var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
2532	var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. &
2533	var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. &
2534	var(1:17) == 'usm_surfwintrans_' .OR. &
2535	var(1:9) == 'usm_qsws_' .OR. var(1:13) == 'usm_qsws_veg_' .OR. &
2536	var(1:13) == 'usm_qsws_liq_' ) THEN
2537	unit = 'W/m2'
2538	ELSE IF ( var(1:15) == 'usm_t_surf_wall' .OR. var(1:10) == 'usm_t_wall' .OR. &
2539	var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR. &
2540	var(1:16) == 'usm_t_surf_green' .OR. &
2541	var(1:11) == 'usm_t_green' .OR. var(1:7) == 'usm_swc' .OR. &
2542	var(1:14) == 'usm_theta_10cm' ) THEN
2543	unit = 'K'
2544	ELSE IF ( var(1:9) == 'usm_surfz' .OR. var(1:11) == 'usm_surfcat' ) THEN
2545	unit = '1'
2546	ELSE
2547	unit = 'illegal'
2548	ENDIF
2549
2550	END SUBROUTINE usm_check_data_output
2551
2552
2553	!------------------------------------------------------------------------------!
2554	! Description:
2555	! ------------
2556	!> Check parameters routine for urban surface model
2557	!------------------------------------------------------------------------------!
2558	SUBROUTINE usm_check_parameters
2559
2560	USE control_parameters, &
2561	ONLY: bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing, &
2562	lsf_surf, topography
2563
2564	USE netcdf_data_input_mod, &
2565	ONLY: building_type_f
2566
2567	IMPLICIT NONE
2568
2569	INTEGER(iwp) :: i !< running index, x-dimension
2570	INTEGER(iwp) :: j !< running index, y-dimension
2571
2572	!
2573	!-- Dirichlet boundary conditions are required as the surface fluxes are
2574	!-- calculated from the temperature/humidity gradients in the urban surface
2575	!-- model
2576	IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN
2577	message_string = 'urban surface model requires setting of '// &
2578	'bc_pt_b = "dirichlet" and '// &
2579	'bc_q_b = "dirichlet"'
2580	CALL message( 'usm_check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
2581	ENDIF
2582
2583	IF ( .NOT. constant_flux_layer ) THEN
2584	message_string = 'urban surface model requires '// &
2585	'constant_flux_layer = .T.'
2586	CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
2587	ENDIF
2588
2589	IF ( .NOT. radiation ) THEN
2590	message_string = 'urban surface model requires '// &
2591	'the radiation model to be switched on'
2592	CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
2593	ENDIF
2594	!
2595	!-- Surface forcing has to be disabled for LSF in case of enabled
2596	!-- urban surface module
2597	IF ( large_scale_forcing ) THEN
2598	lsf_surf = .FALSE.
2599	ENDIF
2600	!
2601	!-- Topography
2602	IF ( topography == 'flat' ) THEN
2603	message_string = 'topography /= "flat" is required '// &
2604	'when using the urban surface model'
2605	CALL message( 'usm_check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2606	ENDIF
2607	!
2608	!-- naheatlayers
2609	IF ( naheatlayers > nzt ) THEN
2610	message_string = 'number of anthropogenic heat layers '// &
2611	'"naheatlayers" can not be larger than'// &
2612	' number of domain layers "nzt"'
2613	CALL message( 'usm_check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
2614	ENDIF
2615	!
2616	!-- Check if building types are set within a valid range.
2617	IF ( building_type < LBOUND( building_pars, 2 ) .AND. &
2618	building_type > UBOUND( building_pars, 2 ) ) THEN
2619	WRITE( message_string, * ) 'building_type = ', building_type, &
2620	' is out of the valid range'
2621	CALL message( 'usm_check_parameters', 'PA0529', 2, 2, 0, 6, 0 )
2622	ENDIF
2623	IF ( building_type_f%from_file ) THEN
2624	DO i = nxl, nxr
2625	DO j = nys, nyn
2626	IF ( building_type_f%var(j,i) /= building_type_f%fill .AND. &
2627	( building_type_f%var(j,i) < LBOUND( building_pars, 2 ) .OR. &
2628	building_type_f%var(j,i) > UBOUND( building_pars, 2 ) ) ) &
2629	THEN
2630	WRITE( message_string, * ) 'building_type = is out of ' // &
2631	'the valid range at (j,i) = ', j, i
2632	CALL message( 'usm_check_parameters', 'PA0529', 2, 2, 0, 6, 0 )
2633	ENDIF
2634	ENDDO
2635	ENDDO
2636	ENDIF
2637	END SUBROUTINE usm_check_parameters
2638
2639
2640	!------------------------------------------------------------------------------!
2641	!
2642	! Description:
2643	! ------------
2644	!> Output of the 3D-arrays in netCDF and/or AVS format
2645	!> for variables of urban_surface model.
2646	!> It resorts the urban surface module output quantities from surf style
2647	!> indexing into temporary 3D array with indices (i,j,k).
2648	!> It is called from subroutine data_output_3d.
2649	!------------------------------------------------------------------------------!
2650	SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2651
2652	IMPLICIT NONE
2653
2654	INTEGER(iwp), INTENT(IN) :: av !< flag if averaged
2655	CHARACTER (len=*), INTENT(IN) :: variable !< variable name
2656	INTEGER(iwp), INTENT(IN) :: nzb_do !< lower limit of the data output (usually 0)
2657	INTEGER(iwp), INTENT(IN) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
2658	LOGICAL, INTENT(OUT) :: found !<
2659	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< sp - it has to correspond to module data_output_3d
2660	REAL(sp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: temp_pf !< temp array for urban surface output procedure
2661
2662	CHARACTER (len=varnamelength) :: var !< trimmed variable name
2663	INTEGER(iwp), PARAMETER :: nd = 5 !< number of directions
2664	CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
2665	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
2666	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: diridx = (/ -1, 1, 0, 3, 2 /)
2667	!< index for surf_*_v: 0:3 = (North, South, East, West)
2668	INTEGER(iwp) :: ids,idsint,idsidx
2669	INTEGER(iwp) :: i,j,k,iwl,istat, l, m !< running indices
2670
2671	found = .TRUE.
2672	temp_pf = -1._wp
2673
2674	ids = -1
2675	var = TRIM(variable)
2676	DO i = 0, nd-1
2677	k = len(TRIM(var))
2678	j = len(TRIM(dirname(i)))
2679	IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) ) THEN
2680	ids = i
2681	idsint = dirint(ids)
2682	idsidx = diridx(ids)
2683	var = var(:k-j)
2684	EXIT
2685	ENDIF
2686	ENDDO
2687	IF ( ids == -1 ) THEN
2688	var = TRIM(variable)
2689	ENDIF
2690	IF ( var(1:11) == 'usm_t_wall_' .AND. len(TRIM(var)) >= 12 ) THEN
2691	!
2692	!-- wall layers
2693	READ(var(12:12), '(I1)', iostat=istat ) iwl
2694	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2695	var = var(1:10)
2696	ENDIF
2697	ENDIF
2698	IF ( var(1:13) == 'usm_t_window_' .AND. len(TRIM(var)) >= 14 ) THEN
2699	!
2700	!-- window layers
2701	READ(var(14:14), '(I1)', iostat=istat ) iwl
2702	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2703	var = var(1:12)
2704	ENDIF
2705	ENDIF
2706	IF ( var(1:12) == 'usm_t_green_' .AND. len(TRIM(var)) >= 13 ) THEN
2707	!
2708	!-- green layers
2709	READ(var(13:13), '(I1)', iostat=istat ) iwl
2710	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2711	var = var(1:11)
2712	ENDIF
2713	ENDIF
2714	IF ( var(1:8) == 'usm_swc_' .AND. len(TRIM(var)) >= 9 ) THEN
2715	!
2716	!-- green layers soil water content
2717	READ(var(9:9), '(I1)', iostat=istat ) iwl
2718	IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN
2719	var = var(1:7)
2720	ENDIF
2721	ENDIF
2722
2723	SELECT CASE ( TRIM(var) )
2724
2725	CASE ( 'usm_surfz' )
2726	!
2727	!-- array of surface height (z)
2728	IF ( idsint == iup_u ) THEN
2729	DO m = 1, surf_usm_h%ns
2730	i = surf_usm_h%i(m)
2731	j = surf_usm_h%j(m)
2732	k = surf_usm_h%k(m)
2733	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, KIND = sp) )
2734	ENDDO
2735	ELSE
2736	l = idsidx
2737	DO m = 1, surf_usm_v(l)%ns
2738	i = surf_usm_v(l)%i(m)
2739	j = surf_usm_v(l)%j(m)
2740	k = surf_usm_v(l)%k(m)
2741	temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, KIND = sp) + 1.0_sp )
2742	ENDDO
2743	ENDIF
2744
2745	CASE ( 'usm_surfcat' )
2746	!
2747	!-- surface category
2748	IF ( idsint == iup_u ) THEN
2749	DO m = 1, surf_usm_h%ns
2750	i = surf_usm_h%i(m)
2751	j = surf_usm_h%j(m)
2752	k = surf_usm_h%k(m)
2753	temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2754	ENDDO
2755	ELSE
2756	l = idsidx
2757	DO m = 1, surf_usm_v(l)%ns
2758	i = surf_usm_v(l)%i(m)
2759	j = surf_usm_v(l)%j(m)
2760	k = surf_usm_v(l)%k(m)
2761	temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
2762	ENDDO
2763	ENDIF
2764
2765	CASE ( 'usm_surfwintrans' )
2766	!
2767	!-- transmissivity window tiles
2768	IF ( idsint == iup_u ) THEN
2769	DO m = 1, surf_usm_h%ns
2770	i = surf_usm_h%i(m)
2771	j = surf_usm_h%j(m)
2772	k = surf_usm_h%k(m)
2773	temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
2774	ENDDO
2775	ELSE
2776	l = idsidx
2777	DO m = 1, surf_usm_v(l)%ns
2778	i = surf_usm_v(l)%i(m)
2779	j = surf_usm_v(l)%j(m)
2780	k = surf_usm_v(l)%k(m)
2781	temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
2782	ENDDO
2783	ENDIF
2784
2785	CASE ( 'usm_wshf' )
2786	!
2787	!-- array of sensible heat flux from surfaces
2788	IF ( av == 0 ) THEN
2789	IF ( idsint == iup_u ) THEN
2790	DO m = 1, surf_usm_h%ns
2791	i = surf_usm_h%i(m)
2792	j = surf_usm_h%j(m)
2793	k = surf_usm_h%k(m)
2794	temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2795	ENDDO
2796	ELSE
2797	l = idsidx
2798	DO m = 1, surf_usm_v(l)%ns
2799	i = surf_usm_v(l)%i(m)
2800	j = surf_usm_v(l)%j(m)
2801	k = surf_usm_v(l)%k(m)
2802	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2803	ENDDO
2804	ENDIF
2805	ELSE
2806	IF ( idsint == iup_u ) THEN
2807	DO m = 1, surf_usm_h%ns
2808	i = surf_usm_h%i(m)
2809	j = surf_usm_h%j(m)
2810	k = surf_usm_h%k(m)
2811	temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2812	ENDDO
2813	ELSE
2814	l = idsidx
2815	DO m = 1, surf_usm_v(l)%ns
2816	i = surf_usm_v(l)%i(m)
2817	j = surf_usm_v(l)%j(m)
2818	k = surf_usm_v(l)%k(m)
2819	temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2820	ENDDO
2821	ENDIF
2822	ENDIF
2823
2824
2825	CASE ( 'usm_qsws' )
2826	!
2827	!-- array of latent heat flux from surfaces
2828	IF ( av == 0 ) THEN
2829	IF ( idsint == iup_u ) THEN
2830	DO m = 1, surf_usm_h%ns
2831	i = surf_usm_h%i(m)
2832	j = surf_usm_h%j(m)
2833	k = surf_usm_h%k(m)
2834	temp_pf(k,j,i) = surf_usm_h%qsws(m) * l_v
2835	ENDDO
2836	ELSE
2837	l = idsidx
2838	DO m = 1, surf_usm_v(l)%ns
2839	i = surf_usm_v(l)%i(m)
2840	j = surf_usm_v(l)%j(m)
2841	k = surf_usm_v(l)%k(m)
2842	temp_pf(k,j,i) = surf_usm_v(l)%qsws(m) * l_v
2843	ENDDO
2844	ENDIF
2845	ELSE
2846	IF ( idsint == iup_u ) THEN
2847	DO m = 1, surf_usm_h%ns
2848	i = surf_usm_h%i(m)
2849	j = surf_usm_h%j(m)
2850	k = surf_usm_h%k(m)
2851	temp_pf(k,j,i) = surf_usm_h%qsws_av(m)
2852	ENDDO
2853	ELSE
2854	l = idsidx
2855	DO m = 1, surf_usm_v(l)%ns
2856	i = surf_usm_v(l)%i(m)
2857	j = surf_usm_v(l)%j(m)
2858	k = surf_usm_v(l)%k(m)
2859	temp_pf(k,j,i) = surf_usm_v(l)%qsws_av(m)
2860	ENDDO
2861	ENDIF
2862	ENDIF
2863
2864	CASE ( 'usm_qsws_veg' )
2865	!
2866	!-- array of latent heat flux from vegetation surfaces
2867	IF ( av == 0 ) THEN
2868	IF ( idsint == iup_u ) THEN
2869	DO m = 1, surf_usm_h%ns
2870	i = surf_usm_h%i(m)
2871	j = surf_usm_h%j(m)
2872	k = surf_usm_h%k(m)
2873	temp_pf(k,j,i) = surf_usm_h%qsws_veg(m)
2874	ENDDO
2875	ELSE
2876	l = idsidx
2877	DO m = 1, surf_usm_v(l)%ns
2878	i = surf_usm_v(l)%i(m)
2879	j = surf_usm_v(l)%j(m)
2880	k = surf_usm_v(l)%k(m)
2881	temp_pf(k,j,i) = surf_usm_v(l)%qsws_veg(m)
2882	ENDDO
2883	ENDIF
2884	ELSE
2885	IF ( idsint == iup_u ) THEN
2886	DO m = 1, surf_usm_h%ns
2887	i = surf_usm_h%i(m)
2888	j = surf_usm_h%j(m)
2889	k = surf_usm_h%k(m)
2890	temp_pf(k,j,i) = surf_usm_h%qsws_veg_av(m)
2891	ENDDO
2892	ELSE
2893	l = idsidx
2894	DO m = 1, surf_usm_v(l)%ns
2895	i = surf_usm_v(l)%i(m)
2896	j = surf_usm_v(l)%j(m)
2897	k = surf_usm_v(l)%k(m)
2898	temp_pf(k,j,i) = surf_usm_v(l)%qsws_veg_av(m)
2899	ENDDO
2900	ENDIF
2901	ENDIF
2902
2903	CASE ( 'usm_qsws_liq' )
2904	!
2905	!-- array of latent heat flux from surfaces with liquid
2906	IF ( av == 0 ) THEN
2907	IF ( idsint == iup_u ) THEN
2908	DO m = 1, surf_usm_h%ns
2909	i = surf_usm_h%i(m)
2910	j = surf_usm_h%j(m)
2911	k = surf_usm_h%k(m)
2912	temp_pf(k,j,i) = surf_usm_h%qsws_liq(m)
2913	ENDDO
2914	ELSE
2915	l = idsidx
2916	DO m = 1, surf_usm_v(l)%ns
2917	i = surf_usm_v(l)%i(m)
2918	j = surf_usm_v(l)%j(m)
2919	k = surf_usm_v(l)%k(m)
2920	temp_pf(k,j,i) = surf_usm_v(l)%qsws_liq(m)
2921	ENDDO
2922	ENDIF
2923	ELSE
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) = surf_usm_h%qsws_liq_av(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) = surf_usm_v(l)%qsws_liq_av(m)
2938	ENDDO
2939	ENDIF
2940	ENDIF
2941
2942	CASE ( 'usm_wghf' )
2943	!
2944	!-- array of heat flux from ground (land, wall, roof)
2945	IF ( av == 0 ) THEN
2946	IF ( idsint == iup_u ) THEN
2947	DO m = 1, surf_usm_h%ns
2948	i = surf_usm_h%i(m)
2949	j = surf_usm_h%j(m)
2950	k = surf_usm_h%k(m)
2951	temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2952	ENDDO
2953	ELSE
2954	l = idsidx
2955	DO m = 1, surf_usm_v(l)%ns
2956	i = surf_usm_v(l)%i(m)
2957	j = surf_usm_v(l)%j(m)
2958	k = surf_usm_v(l)%k(m)
2959	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2960	ENDDO
2961	ENDIF
2962	ELSE
2963	IF ( idsint == iup_u ) THEN
2964	DO m = 1, surf_usm_h%ns
2965	i = surf_usm_h%i(m)
2966	j = surf_usm_h%j(m)
2967	k = surf_usm_h%k(m)
2968	temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2969	ENDDO
2970	ELSE
2971	l = idsidx
2972	DO m = 1, surf_usm_v(l)%ns
2973	i = surf_usm_v(l)%i(m)
2974	j = surf_usm_v(l)%j(m)
2975	k = surf_usm_v(l)%k(m)
2976	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2977	ENDDO
2978	ENDIF
2979	ENDIF
2980
2981	CASE ( 'usm_wghf_window' )
2982	!
2983	!-- array of heat flux from window ground (land, wall, roof)
2984	IF ( av == 0 ) THEN
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%wghf_eb_window(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)%wghf_eb_window(m)
2999	ENDDO
3000	ENDIF
3001	ELSE
3002	IF ( idsint == iup_u ) THEN
3003	DO m = 1, surf_usm_h%ns
3004	i = surf_usm_h%i(m)
3005	j = surf_usm_h%j(m)
3006	k = surf_usm_h%k(m)
3007	temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
3008	ENDDO
3009	ELSE
3010	l = idsidx
3011	DO m = 1, surf_usm_v(l)%ns
3012	i = surf_usm_v(l)%i(m)
3013	j = surf_usm_v(l)%j(m)
3014	k = surf_usm_v(l)%k(m)
3015	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
3016	ENDDO
3017	ENDIF
3018	ENDIF
3019
3020	CASE ( 'usm_wghf_green' )
3021	!
3022	!-- array of heat flux from green ground (land, wall, roof)
3023	IF ( av == 0 ) THEN
3024	IF ( idsint == iup_u ) THEN
3025	DO m = 1, surf_usm_h%ns
3026	i = surf_usm_h%i(m)
3027	j = surf_usm_h%j(m)
3028	k = surf_usm_h%k(m)
3029	temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
3030	ENDDO
3031	ELSE
3032	l = idsidx
3033	DO m = 1, surf_usm_v(l)%ns
3034	i = surf_usm_v(l)%i(m)
3035	j = surf_usm_v(l)%j(m)
3036	k = surf_usm_v(l)%k(m)
3037	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
3038	ENDDO
3039	ENDIF
3040	ELSE
3041	IF ( idsint == iup_u ) THEN
3042	DO m = 1, surf_usm_h%ns
3043	i = surf_usm_h%i(m)
3044	j = surf_usm_h%j(m)
3045	k = surf_usm_h%k(m)
3046	temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
3047	ENDDO
3048	ELSE
3049	l = idsidx
3050	DO m = 1, surf_usm_v(l)%ns
3051	i = surf_usm_v(l)%i(m)
3052	j = surf_usm_v(l)%j(m)
3053	k = surf_usm_v(l)%k(m)
3054	temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m)
3055	ENDDO
3056	ENDIF
3057	ENDIF
3058
3059	CASE ( 'usm_iwghf' )
3060	!
3061	!-- array of heat flux from indoor ground (land, wall, roof)
3062	IF ( av == 0 ) THEN
3063	IF ( idsint == iup_u ) THEN
3064	DO m = 1, surf_usm_h%ns
3065	i = surf_usm_h%i(m)
3066	j = surf_usm_h%j(m)
3067	k = surf_usm_h%k(m)
3068	temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
3069	ENDDO
3070	ELSE
3071	l = idsidx
3072	DO m = 1, surf_usm_v(l)%ns
3073	i = surf_usm_v(l)%i(m)
3074	j = surf_usm_v(l)%j(m)
3075	k = surf_usm_v(l)%k(m)
3076	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m)
3077	ENDDO
3078	ENDIF
3079	ELSE
3080	IF ( idsint == iup_u ) THEN
3081	DO m = 1, surf_usm_h%ns
3082	i = surf_usm_h%i(m)
3083	j = surf_usm_h%j(m)
3084	k = surf_usm_h%k(m)
3085	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
3086	ENDDO
3087	ELSE
3088	l = idsidx
3089	DO m = 1, surf_usm_v(l)%ns
3090	i = surf_usm_v(l)%i(m)
3091	j = surf_usm_v(l)%j(m)
3092	k = surf_usm_v(l)%k(m)
3093	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m)
3094	ENDDO
3095	ENDIF
3096	ENDIF
3097
3098	CASE ( 'usm_iwghf_window' )
3099	!
3100	!-- array of heat flux from indoor window ground (land, wall, roof)
3101	IF ( av == 0 ) THEN
3102	IF ( idsint == iup_u ) THEN
3103	DO m = 1, surf_usm_h%ns
3104	i = surf_usm_h%i(m)
3105	j = surf_usm_h%j(m)
3106	k = surf_usm_h%k(m)
3107	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
3108	ENDDO
3109	ELSE
3110	l = idsidx
3111	DO m = 1, surf_usm_v(l)%ns
3112	i = surf_usm_v(l)%i(m)
3113	j = surf_usm_v(l)%j(m)
3114	k = surf_usm_v(l)%k(m)
3115	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m)
3116	ENDDO
3117	ENDIF
3118	ELSE
3119	IF ( idsint == iup_u ) THEN
3120	DO m = 1, surf_usm_h%ns
3121	i = surf_usm_h%i(m)
3122	j = surf_usm_h%j(m)
3123	k = surf_usm_h%k(m)
3124	temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
3125	ENDDO
3126	ELSE
3127	l = idsidx
3128	DO m = 1, surf_usm_v(l)%ns
3129	i = surf_usm_v(l)%i(m)
3130	j = surf_usm_v(l)%j(m)
3131	k = surf_usm_v(l)%k(m)
3132	temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m)
3133	ENDDO
3134	ENDIF
3135	ENDIF
3136
3137	CASE ( 'usm_t_surf_wall' )
3138	!
3139	!-- surface temperature for surfaces
3140	IF ( av == 0 ) THEN
3141	IF ( idsint == iup_u ) THEN
3142	DO m = 1, surf_usm_h%ns
3143	i = surf_usm_h%i(m)
3144	j = surf_usm_h%j(m)
3145	k = surf_usm_h%k(m)
3146	temp_pf(k,j,i) = t_surf_wall_h(m)
3147	ENDDO
3148	ELSE
3149	l = idsidx
3150	DO m = 1, surf_usm_v(l)%ns
3151	i = surf_usm_v(l)%i(m)
3152	j = surf_usm_v(l)%j(m)
3153	k = surf_usm_v(l)%k(m)
3154	temp_pf(k,j,i) = t_surf_wall_v(l)%t(m)
3155	ENDDO
3156	ENDIF
3157	ELSE
3158	IF ( idsint == iup_u ) THEN
3159	DO m = 1, surf_usm_h%ns
3160	i = surf_usm_h%i(m)
3161	j = surf_usm_h%j(m)
3162	k = surf_usm_h%k(m)
3163	temp_pf(k,j,i) = surf_usm_h%t_surf_wall_av(m)
3164	ENDDO
3165	ELSE
3166	l = idsidx
3167	DO m = 1, surf_usm_v(l)%ns
3168	i = surf_usm_v(l)%i(m)
3169	j = surf_usm_v(l)%j(m)
3170	k = surf_usm_v(l)%k(m)
3171	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_wall_av(m)
3172	ENDDO
3173	ENDIF
3174	ENDIF
3175
3176	CASE ( 'usm_t_surf_window' )
3177	!
3178	!-- surface temperature for window surfaces
3179	IF ( av == 0 ) THEN
3180	IF ( idsint == iup_u ) THEN
3181	DO m = 1, surf_usm_h%ns
3182	i = surf_usm_h%i(m)
3183	j = surf_usm_h%j(m)
3184	k = surf_usm_h%k(m)
3185	temp_pf(k,j,i) = t_surf_window_h(m)
3186	ENDDO
3187	ELSE
3188	l = idsidx
3189	DO m = 1, surf_usm_v(l)%ns
3190	i = surf_usm_v(l)%i(m)
3191	j = surf_usm_v(l)%j(m)
3192	k = surf_usm_v(l)%k(m)
3193	temp_pf(k,j,i) = t_surf_window_v(l)%t(m)
3194	ENDDO
3195	ENDIF
3196
3197	ELSE
3198	IF ( idsint == iup_u ) THEN
3199	DO m = 1, surf_usm_h%ns
3200	i = surf_usm_h%i(m)
3201	j = surf_usm_h%j(m)
3202	k = surf_usm_h%k(m)
3203	temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
3204	ENDDO
3205	ELSE
3206	l = idsidx
3207	DO m = 1, surf_usm_v(l)%ns
3208	i = surf_usm_v(l)%i(m)
3209	j = surf_usm_v(l)%j(m)
3210	k = surf_usm_v(l)%k(m)
3211	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_window_av(m)
3212	ENDDO
3213
3214	ENDIF
3215
3216	ENDIF
3217
3218	CASE ( 'usm_t_surf_green' )
3219	!
3220	!-- surface temperature for green surfaces
3221	IF ( av == 0 ) THEN
3222	IF ( idsint == iup_u ) THEN
3223	DO m = 1, surf_usm_h%ns
3224	i = surf_usm_h%i(m)
3225	j = surf_usm_h%j(m)
3226	k = surf_usm_h%k(m)
3227	temp_pf(k,j,i) = t_surf_green_h(m)
3228	ENDDO
3229	ELSE
3230	l = idsidx
3231	DO m = 1, surf_usm_v(l)%ns
3232	i = surf_usm_v(l)%i(m)
3233	j = surf_usm_v(l)%j(m)
3234	k = surf_usm_v(l)%k(m)
3235	temp_pf(k,j,i) = t_surf_green_v(l)%t(m)
3236	ENDDO
3237	ENDIF
3238
3239	ELSE
3240	IF ( idsint == iup_u ) THEN
3241	DO m = 1, surf_usm_h%ns
3242	i = surf_usm_h%i(m)
3243	j = surf_usm_h%j(m)
3244	k = surf_usm_h%k(m)
3245	temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
3246	ENDDO
3247	ELSE
3248	l = idsidx
3249	DO m = 1, surf_usm_v(l)%ns
3250	i = surf_usm_v(l)%i(m)
3251	j = surf_usm_v(l)%j(m)
3252	k = surf_usm_v(l)%k(m)
3253	temp_pf(k,j,i) = surf_usm_v(l)%t_surf_green_av(m)
3254	ENDDO
3255
3256	ENDIF
3257
3258	ENDIF
3259
3260	CASE ( 'usm_theta_10cm' )
3261	!
3262	!-- near surface temperature for whole surfaces
3263	IF ( av == 0 ) THEN
3264	IF ( idsint == iup_u ) THEN
3265	DO m = 1, surf_usm_h%ns
3266	i = surf_usm_h%i(m)
3267	j = surf_usm_h%j(m)
3268	k = surf_usm_h%k(m)
3269	temp_pf(k,j,i) = surf_usm_h%pt_10cm(m)
3270	ENDDO
3271	ELSE
3272	l = idsidx
3273	DO m = 1, surf_usm_v(l)%ns
3274	i = surf_usm_v(l)%i(m)
3275	j = surf_usm_v(l)%j(m)
3276	k = surf_usm_v(l)%k(m)
3277	temp_pf(k,j,i) = surf_usm_v(l)%pt_10cm(m)
3278	ENDDO
3279	ENDIF
3280
3281
3282	ELSE
3283	IF ( idsint == iup_u ) THEN
3284	DO m = 1, surf_usm_h%ns
3285	i = surf_usm_h%i(m)
3286	j = surf_usm_h%j(m)
3287	k = surf_usm_h%k(m)
3288	temp_pf(k,j,i) = surf_usm_h%pt_10cm_av(m)
3289	ENDDO
3290	ELSE
3291	l = idsidx
3292	DO m = 1, surf_usm_v(l)%ns
3293	i = surf_usm_v(l)%i(m)
3294	j = surf_usm_v(l)%j(m)
3295	k = surf_usm_v(l)%k(m)
3296	temp_pf(k,j,i) = surf_usm_v(l)%pt_10cm_av(m)
3297	ENDDO
3298
3299	ENDIF
3300	ENDIF
3301
3302	CASE ( 'usm_t_wall' )
3303	!
3304	!-- wall temperature for iwl layer of walls and land
3305	IF ( av == 0 ) THEN
3306	IF ( idsint == iup_u ) THEN
3307	DO m = 1, surf_usm_h%ns
3308	i = surf_usm_h%i(m)
3309	j = surf_usm_h%j(m)
3310	k = surf_usm_h%k(m)
3311	temp_pf(k,j,i) = t_wall_h(iwl,m)
3312	ENDDO
3313	ELSE
3314	l = idsidx
3315	DO m = 1, surf_usm_v(l)%ns
3316	i = surf_usm_v(l)%i(m)
3317	j = surf_usm_v(l)%j(m)
3318	k = surf_usm_v(l)%k(m)
3319	temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
3320	ENDDO
3321	ENDIF
3322	ELSE
3323	IF ( idsint == iup_u ) THEN
3324	DO m = 1, surf_usm_h%ns
3325	i = surf_usm_h%i(m)
3326	j = surf_usm_h%j(m)
3327	k = surf_usm_h%k(m)
3328	temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
3329	ENDDO
3330	ELSE
3331	l = idsidx
3332	DO m = 1, surf_usm_v(l)%ns
3333	i = surf_usm_v(l)%i(m)
3334	j = surf_usm_v(l)%j(m)
3335	k = surf_usm_v(l)%k(m)
3336	temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
3337	ENDDO
3338	ENDIF
3339	ENDIF
3340
3341	CASE ( 'usm_t_window' )
3342	!
3343	!-- window temperature for iwl layer of walls and land
3344	IF ( av == 0 ) THEN
3345	IF ( idsint == iup_u ) THEN
3346	DO m = 1, surf_usm_h%ns
3347	i = surf_usm_h%i(m)
3348	j = surf_usm_h%j(m)
3349	k = surf_usm_h%k(m)
3350	temp_pf(k,j,i) = t_window_h(iwl,m)
3351	ENDDO
3352	ELSE
3353	l = idsidx
3354	DO m = 1, surf_usm_v(l)%ns
3355	i = surf_usm_v(l)%i(m)
3356	j = surf_usm_v(l)%j(m)
3357	k = surf_usm_v(l)%k(m)
3358	temp_pf(k,j,i) = t_window_v(l)%t(iwl,m)
3359	ENDDO
3360	ENDIF
3361	ELSE
3362	IF ( idsint == iup_u ) THEN
3363	DO m = 1, surf_usm_h%ns
3364	i = surf_usm_h%i(m)
3365	j = surf_usm_h%j(m)
3366	k = surf_usm_h%k(m)
3367	temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
3368	ENDDO
3369	ELSE
3370	l = idsidx
3371	DO m = 1, surf_usm_v(l)%ns
3372	i = surf_usm_v(l)%i(m)
3373	j = surf_usm_v(l)%j(m)
3374	k = surf_usm_v(l)%k(m)
3375	temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m)
3376	ENDDO
3377	ENDIF
3378	ENDIF
3379
3380	CASE ( 'usm_t_green' )
3381	!
3382	!-- green temperature for iwl layer of walls and land
3383	IF ( av == 0 ) THEN
3384	IF ( idsint == iup_u ) THEN
3385	DO m = 1, surf_usm_h%ns
3386	i = surf_usm_h%i(m)
3387	j = surf_usm_h%j(m)
3388	k = surf_usm_h%k(m)
3389	temp_pf(k,j,i) = t_green_h(iwl,m)
3390	ENDDO
3391	ELSE
3392	l = idsidx
3393	DO m = 1, surf_usm_v(l)%ns
3394	i = surf_usm_v(l)%i(m)
3395	j = surf_usm_v(l)%j(m)
3396	k = surf_usm_v(l)%k(m)
3397	temp_pf(k,j,i) = t_green_v(l)%t(iwl,m)
3398	ENDDO
3399	ENDIF
3400	ELSE
3401	IF ( idsint == iup_u ) THEN
3402	DO m = 1, surf_usm_h%ns
3403	i = surf_usm_h%i(m)
3404	j = surf_usm_h%j(m)
3405	k = surf_usm_h%k(m)
3406	temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
3407	ENDDO
3408	ELSE
3409	l = idsidx
3410	DO m = 1, surf_usm_v(l)%ns
3411	i = surf_usm_v(l)%i(m)
3412	j = surf_usm_v(l)%j(m)
3413	k = surf_usm_v(l)%k(m)
3414	temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m)
3415	ENDDO
3416	ENDIF
3417	ENDIF
3418
3419	CASE ( 'usm_swc' )
3420	!
3421	!-- soil water content for iwl layer of walls and land
3422	IF ( av == 0 ) THEN
3423	IF ( idsint == iup_u ) THEN
3424	DO m = 1, surf_usm_h%ns
3425	i = surf_usm_h%i(m)
3426	j = surf_usm_h%j(m)
3427	k = surf_usm_h%k(m)
3428	temp_pf(k,j,i) = swc_h(iwl,m)
3429	ENDDO
3430	ELSE
3431	l = idsidx
3432	DO m = 1, surf_usm_v(l)%ns
3433	i = surf_usm_v(l)%i(m)
3434	j = surf_usm_v(l)%j(m)
3435	k = surf_usm_v(l)%k(m)
3436	temp_pf(k,j,i) = swc_v(l)%t(iwl,m)
3437	ENDDO
3438	ENDIF
3439	ELSE
3440	IF ( idsint == iup_u ) THEN
3441	DO m = 1, surf_usm_h%ns
3442	i = surf_usm_h%i(m)
3443	j = surf_usm_h%j(m)
3444	k = surf_usm_h%k(m)
3445	temp_pf(k,j,i) = surf_usm_h%swc_av(iwl,m)
3446	ENDDO
3447	ELSE
3448	l = idsidx
3449	DO m = 1, surf_usm_v(l)%ns
3450	i = surf_usm_v(l)%i(m)
3451	j = surf_usm_v(l)%j(m)
3452	k = surf_usm_v(l)%k(m)
3453	temp_pf(k,j,i) = surf_usm_v(l)%swc_av(iwl,m)
3454	ENDDO
3455	ENDIF
3456	ENDIF
3457
3458
3459	CASE DEFAULT
3460	found = .FALSE.
3461	RETURN
3462	END SELECT
3463
3464	!
3465	!-- Rearrange dimensions for NetCDF output
3466	!-- FIXME: this may generate FPE overflow upon conversion from DP to SP
3467	DO j = nys, nyn
3468	DO i = nxl, nxr
3469	DO k = nzb_do, nzt_do
3470	local_pf(i,j,k) = temp_pf(k,j,i)
3471	ENDDO
3472	ENDDO
3473	ENDDO
3474
3475	END SUBROUTINE usm_data_output_3d
3476
3477
3478	!------------------------------------------------------------------------------!
3479	!
3480	! Description:
3481	! ------------
3482	!> Soubroutine defines appropriate grid for netcdf variables.
3483	!> It is called out from subroutine netcdf.
3484	!------------------------------------------------------------------------------!
3485	SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
3486
3487	IMPLICIT NONE
3488
3489	CHARACTER (len=*), INTENT(IN) :: variable !<
3490	LOGICAL, INTENT(OUT) :: found !<
3491	CHARACTER (len=*), INTENT(OUT) :: grid_x !<
3492	CHARACTER (len=*), INTENT(OUT) :: grid_y !<
3493	CHARACTER (len=*), INTENT(OUT) :: grid_z !<
3494
3495	CHARACTER (len=varnamelength) :: var
3496
3497	var = TRIM(variable)
3498	IF ( var(1:9) == 'usm_wshf_' .OR. var(1:9) == 'usm_wghf_' .OR. &
3499	var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR. &
3500	var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR. &
3501	var(1:9) == 'usm_qsws_' .OR. var(1:13) == 'usm_qsws_veg_' .OR. &
3502	var(1:13) == 'usm_qsws_liq_' .OR. &
3503	var(1:15) == 'usm_t_surf_wall' .OR. var(1:10) == 'usm_t_wall' .OR. &
3504	var(1:17) == 'usm_t_surf_window' .OR. var(1:12) == 'usm_t_window' .OR. &
3505	var(1:16) == 'usm_t_surf_green' .OR. var(1:11) == 'usm_t_green' .OR. &
3506	var(1:15) == 'usm_theta_10cm' .OR. &
3507	var(1:9) == 'usm_surfz' .OR. var(1:11) == 'usm_surfcat' .OR. &
3508	var(1:16) == 'usm_surfwintrans' .OR. var(1:7) == 'usm_swc' ) THEN
3509
3510	found = .TRUE.
3511	grid_x = 'x'
3512	grid_y = 'y'
3513	grid_z = 'zu'
3514	ELSE
3515	found = .FALSE.
3516	grid_x = 'none'
3517	grid_y = 'none'
3518	grid_z = 'none'
3519	ENDIF
3520
3521	END SUBROUTINE usm_define_netcdf_grid
3522
3523
3524	!------------------------------------------------------------------------------!
3525	! Description:
3526	! ------------
3527	!> Initialization of the wall surface model
3528	!------------------------------------------------------------------------------!
3529	SUBROUTINE usm_init_material_model
3530
3531	IMPLICIT NONE
3532
3533	INTEGER(iwp) :: k, l, m !< running indices
3534
3535	IF ( debug_output ) CALL debug_message( 'usm_init_material_model', 'start' )
3536
3537	!
3538	!-- Calculate wall grid spacings.
3539	!-- Temperature is defined at the center of the wall layers,
3540	!-- whereas gradients/fluxes are defined at the edges (_stag)
3541	!-- apply for all particular surface grids. First for horizontal surfaces
3542	DO m = 1, surf_usm_h%ns
3543
3544	surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
3545	DO k = nzb_wall+1, nzt_wall
3546	surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) - &
3547	surf_usm_h%zw(k-1,m)
3548	ENDDO
3549	surf_usm_h%dz_window(nzb_wall,m) = surf_usm_h%zw_window(nzb_wall,m)
3550	DO k = nzb_wall+1, nzt_wall
3551	surf_usm_h%dz_window(k,m) = surf_usm_h%zw_window(k,m) - &
3552	surf_usm_h%zw_window(k-1,m)
3553	ENDDO
3554
3555	surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)
3556
3557	DO k = nzb_wall, nzt_wall-1
3558	surf_usm_h%dz_wall_stag(k,m) = 0.5 * ( &
3559	surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
3560	ENDDO
3561	surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
3562
3563	surf_usm_h%dz_window(nzt_wall+1,m) = surf_usm_h%dz_window(nzt_wall,m)
3564
3565	DO k = nzb_wall, nzt_wall-1
3566	surf_usm_h%dz_window_stag(k,m) = 0.5 * ( &
3567	surf_usm_h%dz_window(k+1,m) + surf_usm_h%dz_window(k,m) )
3568	ENDDO
3569	surf_usm_h%dz_window_stag(nzt_wall,m) = surf_usm_h%dz_window(nzt_wall,m)
3570
3571	IF (surf_usm_h%green_type_roof(m) == 2.0_wp ) THEN
3572	!
3573	!-- extensive green roof
3574	!-- set ratio of substrate layer thickness, soil-type and LAI
3575	soil_type = 3
3576	surf_usm_h%lai(m) = 2.0_wp
3577
3578	surf_usm_h%zw_green(nzb_wall,m) = 0.05_wp
3579	surf_usm_h%zw_green(nzb_wall+1,m) = 0.10_wp
3580	surf_usm_h%zw_green(nzb_wall+2,m) = 0.15_wp
3581	surf_usm_h%zw_green(nzb_wall+3,m) = 0.20_wp
3582	ELSE
3583	!
3584	!-- intensiv green roof
3585	!-- set ratio of substrate layer thickness, soil-type and LAI
3586	soil_type = 6
3587	surf_usm_h%lai(m) = 4.0_wp
3588
3589	surf_usm_h%zw_green(nzb_wall,m) = 0.05_wp
3590	surf_usm_h%zw_green(nzb_wall+1,m) = 0.10_wp
3591	surf_usm_h%zw_green(nzb_wall+2,m) = 0.40_wp
3592	surf_usm_h%zw_green(nzb_wall+3,m) = 0.80_wp
3593	ENDIF
3594
3595	surf_usm_h%dz_green(nzb_wall,m) = surf_usm_h%zw_green(nzb_wall,m)
3596	DO k = nzb_wall+1, nzt_wall
3597	surf_usm_h%dz_green(k,m) = surf_usm_h%zw_green(k,m) - &
3598	surf_usm_h%zw_green(k-1,m)
3599	ENDDO
3600	surf_usm_h%dz_green(nzt_wall+1,m) = surf_usm_h%dz_green(nzt_wall,m)
3601
3602	DO k = nzb_wall, nzt_wall-1
3603	surf_usm_h%dz_green_stag(k,m) = 0.5 * ( &
3604	surf_usm_h%dz_green(k+1,m) + surf_usm_h%dz_green(k,m) )
3605	ENDDO
3606	surf_usm_h%dz_green_stag(nzt_wall,m) = surf_usm_h%dz_green(nzt_wall,m)
3607
3608	IF ( alpha_vangenuchten == 9999999.9_wp ) THEN
3609	alpha_vangenuchten = soil_pars(0,soil_type)
3610	ENDIF
3611
3612	IF ( l_vangenuchten == 9999999.9_wp ) THEN
3613	l_vangenuchten = soil_pars(1,soil_type)
3614	ENDIF
3615
3616	IF ( n_vangenuchten == 9999999.9_wp ) THEN
3617	n_vangenuchten = soil_pars(2,soil_type)
3618	ENDIF
3619
3620	IF ( hydraulic_conductivity == 9999999.9_wp ) THEN
3621	hydraulic_conductivity = soil_pars(3,soil_type)
3622	ENDIF
3623
3624	IF ( saturation_moisture == 9999999.9_wp ) THEN
3625	saturation_moisture = m_soil_pars(0,soil_type)
3626	ENDIF
3627
3628	IF ( field_capacity == 9999999.9_wp ) THEN
3629	field_capacity = m_soil_pars(1,soil_type)
3630	ENDIF
3631
3632	IF ( wilting_point == 9999999.9_wp ) THEN
3633	wilting_point = m_soil_pars(2,soil_type)
3634	ENDIF
3635
3636	IF ( residual_moisture == 9999999.9_wp ) THEN
3637	residual_moisture = m_soil_pars(3,soil_type)
3638	ENDIF
3639
3640	DO k = nzb_wall, nzt_wall+1
3641	swc_h(k,m) = field_capacity
3642	rootfr_h(k,m) = 0.5_wp
3643	surf_usm_h%alpha_vg_green(m) = alpha_vangenuchten
3644	surf_usm_h%l_vg_green(m) = l_vangenuchten
3645	surf_usm_h%n_vg_green(m) = n_vangenuchten
3646	surf_usm_h%gamma_w_green_sat(k,m) = hydraulic_conductivity
3647	swc_sat_h(k,m) = saturation_moisture
3648	fc_h(k,m) = field_capacity
3649	wilt_h(k,m) = wilting_point
3650	swc_res_h(k,m) = residual_moisture
3651	ENDDO
3652
3653	ENDDO
3654
3655	surf_usm_h%ddz_wall = 1.0_wp / surf_usm_h%dz_wall
3656	surf_usm_h%ddz_wall_stag = 1.0_wp / surf_usm_h%dz_wall_stag
3657	surf_usm_h%ddz_window = 1.0_wp / surf_usm_h%dz_window
3658	surf_usm_h%ddz_window_stag = 1.0_wp / surf_usm_h%dz_window_stag
3659	surf_usm_h%ddz_green = 1.0_wp / surf_usm_h%dz_green
3660	surf_usm_h%ddz_green_stag = 1.0_wp / surf_usm_h%dz_green_stag
3661	!
3662	!-- For vertical surfaces
3663	DO l = 0, 3
3664	DO m = 1, surf_usm_v(l)%ns
3665	surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
3666	DO k = nzb_wall+1, nzt_wall
3667	surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) - &
3668	surf_usm_v(l)%zw(k-1,m)
3669	ENDDO
3670	surf_usm_v(l)%dz_window(nzb_wall,m) = surf_usm_v(l)%zw_window(nzb_wall,m)
3671	DO k = nzb_wall+1, nzt_wall
3672	surf_usm_v(l)%dz_window(k,m) = surf_usm_v(l)%zw_window(k,m) - &
3673	surf_usm_v(l)%zw_window(k-1,m)
3674	ENDDO
3675	surf_usm_v(l)%dz_green(nzb_wall,m) = surf_usm_v(l)%zw_green(nzb_wall,m)
3676	DO k = nzb_wall+1, nzt_wall
3677	surf_usm_v(l)%dz_green(k,m) = surf_usm_v(l)%zw_green(k,m) - &
3678	surf_usm_v(l)%zw_green(k-1,m)
3679	ENDDO
3680
3681	surf_usm_v(l)%dz_wall(nzt_wall+1,m) = &
3682	surf_usm_v(l)%dz_wall(nzt_wall,m)
3683
3684	DO k = nzb_wall, nzt_wall-1
3685	surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * ( &
3686	surf_usm_v(l)%dz_wall(k+1,m) + &
3687	surf_usm_v(l)%dz_wall(k,m) )
3688	ENDDO
3689	surf_usm_v(l)%dz_wall_stag(nzt_wall,m) = &
3690	surf_usm_v(l)%dz_wall(nzt_wall,m)
3691	surf_usm_v(l)%dz_window(nzt_wall+1,m) = &
3692	surf_usm_v(l)%dz_window(nzt_wall,m)
3693
3694	DO k = nzb_wall, nzt_wall-1
3695	surf_usm_v(l)%dz_window_stag(k,m) = 0.5 * ( &
3696	surf_usm_v(l)%dz_window(k+1,m) + &
3697	surf_usm_v(l)%dz_window(k,m) )
3698	ENDDO
3699	surf_usm_v(l)%dz_window_stag(nzt_wall,m) = &
3700	surf_usm_v(l)%dz_window(nzt_wall,m)
3701	surf_usm_v(l)%dz_green(nzt_wall+1,m) = &
3702	surf_usm_v(l)%dz_green(nzt_wall,m)
3703
3704	DO k = nzb_wall, nzt_wall-1
3705	surf_usm_v(l)%dz_green_stag(k,m) = 0.5 * ( &
3706	surf_usm_v(l)%dz_green(k+1,m) + &
3707	surf_usm_v(l)%dz_green(k,m) )
3708	ENDDO
3709	surf_usm_v(l)%dz_green_stag(nzt_wall,m) = &
3710	surf_usm_v(l)%dz_green(nzt_wall,m)
3711	ENDDO
3712	surf_usm_v(l)%ddz_wall = 1.0_wp / surf_usm_v(l)%dz_wall
3713	surf_usm_v(l)%ddz_wall_stag = 1.0_wp / surf_usm_v(l)%dz_wall_stag
3714	surf_usm_v(l)%ddz_window = 1.0_wp / surf_usm_v(l)%dz_window
3715	surf_usm_v(l)%ddz_window_stag = 1.0_wp / surf_usm_v(l)%dz_window_stag
3716	surf_usm_v(l)%ddz_green = 1.0_wp / surf_usm_v(l)%dz_green
3717	surf_usm_v(l)%ddz_green_stag = 1.0_wp / surf_usm_v(l)%dz_green_stag
3718	ENDDO
3719
3720
3721	IF ( debug_output ) CALL debug_message( 'usm_init_material_model', 'end' )
3722
3723	END SUBROUTINE usm_init_material_model
3724
3725
3726	!------------------------------------------------------------------------------!
3727	! Description:
3728	! ------------
3729	!> Initialization of the urban surface model
3730	!------------------------------------------------------------------------------!
3731	SUBROUTINE usm_init
3732
3733	USE arrays_3d, &
3734	ONLY: zw
3735
3736	USE netcdf_data_input_mod, &
3737	ONLY: building_pars_f, building_type_f, terrain_height_f
3738
3739	IMPLICIT NONE
3740
3741	INTEGER(iwp) :: i !< loop index x-dirction
3742	INTEGER(iwp) :: ind_alb_green !< index in input list for green albedo
3743	INTEGER(iwp) :: ind_alb_wall !< index in input list for wall albedo
3744	INTEGER(iwp) :: ind_alb_win !< index in input list for window albedo
3745	INTEGER(iwp) :: ind_emis_wall !< index in input list for wall emissivity
3746	INTEGER(iwp) :: ind_emis_green !< index in input list for green emissivity
3747	INTEGER(iwp) :: ind_emis_win !< index in input list for window emissivity
3748	INTEGER(iwp) :: ind_green_frac_w !< index in input list for green fraction on wall
3749	INTEGER(iwp) :: ind_green_frac_r !< index in input list for green fraction on roof
3750	INTEGER(iwp) :: ind_hc1 !< index in input list for heat capacity at first wall layer
3751	INTEGER(iwp) :: ind_hc1_win !< index in input list for heat capacity at first window layer
3752	INTEGER(iwp) :: ind_hc2 !< index in input list for heat capacity at second wall layer
3753	INTEGER(iwp) :: ind_hc2_win !< index in input list for heat capacity at second window layer
3754	INTEGER(iwp) :: ind_hc3 !< index in input list for heat capacity at third wall layer
3755	INTEGER(iwp) :: ind_hc3_win !< index in input list for heat capacity at third window layer
3756	INTEGER(iwp) :: ind_lai_r !< index in input list for LAI on roof
3757	INTEGER(iwp) :: ind_lai_w !< index in input list for LAI on wall
3758	INTEGER(iwp) :: ind_tc1 !< index in input list for thermal conductivity at first wall layer
3759	INTEGER(iwp) :: ind_tc1_win !< index in input list for thermal conductivity at first window layer
3760	INTEGER(iwp) :: ind_tc2 !< index in input list for thermal conductivity at second wall layer
3761	INTEGER(iwp) :: ind_tc2_win !< index in input list for thermal conductivity at second window layer
3762	INTEGER(iwp) :: ind_tc3 !< index in input list for thermal conductivity at third wall layer
3763	INTEGER(iwp) :: ind_tc3_win !< index in input list for thermal conductivity at third window layer
3764	INTEGER(iwp) :: ind_thick_1 !< index in input list for thickness of first wall layer
3765	INTEGER(iwp) :: ind_thick_1_win !< index in input list for thickness of first window layer
3766	INTEGER(iwp) :: ind_thick_2 !< index in input list for thickness of second wall layer
3767	INTEGER(iwp) :: ind_thick_2_win !< index in input list for thickness of second window layer
3768	INTEGER(iwp) :: ind_thick_3 !< index in input list for thickness of third wall layer
3769	INTEGER(iwp) :: ind_thick_3_win !< index in input list for thickness of third window layer
3770	INTEGER(iwp) :: ind_thick_4 !< index in input list for thickness of fourth wall layer
3771	INTEGER(iwp) :: ind_thick_4_win !< index in input list for thickness of fourth window layer
3772	INTEGER(iwp) :: ind_trans !< index in input list for window transmissivity
3773	INTEGER(iwp) :: ind_wall_frac !< index in input list for wall fraction
3774	INTEGER(iwp) :: ind_win_frac !< index in input list for window fraction
3775	INTEGER(iwp) :: ind_z0 !< index in input list for z0
3776	INTEGER(iwp) :: ind_z0qh !< index in input list for z0h / z0q
3777	INTEGER(iwp) :: j !< loop index y-dirction
3778	INTEGER(iwp) :: k !< loop index z-dirction
3779	INTEGER(iwp) :: l !< loop index surface orientation
3780	INTEGER(iwp) :: m !< loop index surface element
3781	INTEGER(iwp) :: st !< dummy
3782
3783	REAL(wp) :: c, tin, twin
3784	REAL(wp) :: ground_floor_level_l !< local height of ground floor level
3785	REAL(wp) :: z_agl !< height above ground
3786
3787	IF ( debug_output ) CALL debug_message( 'usm_init', 'start' )
3788
3789	CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
3790	!
3791	!-- surface forcing have to be disabled for LSF
3792	!-- in case of enabled urban surface module
3793	IF ( large_scale_forcing ) THEN
3794	lsf_surf = .FALSE.
3795	ENDIF
3796	!
3797	!-- Flag surface elements belonging to the ground floor level. Therefore,
3798	!-- use terrain height array from file, if available. This flag is later used
3799	!-- to control initialization of surface attributes.
3800	!-- Todo: for the moment disable initialization of building roofs with
3801	!-- ground-floor-level properties.
3802	surf_usm_h%ground_level = .FALSE.
3803
3804	DO l = 0, 3
3805	surf_usm_v(l)%ground_level = .FALSE.
3806	DO m = 1, surf_usm_v(l)%ns
3807	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
3808	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
3809	k = surf_usm_v(l)%k(m)
3810	!
3811	!-- Determine local ground level. Level 1 - default value,
3812	!-- level 2 - initialization according to building type,
3813	!-- level 3 - initialization from value read from file.
3814	ground_floor_level_l = ground_floor_level
3815
3816	IF ( building_type_f%from_file ) THEN
3817	ground_floor_level_l = &
3818	building_pars(ind_gflh,building_type_f%var(j,i))
3819	ENDIF
3820
3821	IF ( building_pars_f%from_file ) THEN
3822	IF ( building_pars_f%pars_xy(ind_gflh,j,i) /= &
3823	building_pars_f%fill ) &
3824	ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)
3825	ENDIF
3826	!
3827	!-- Determine height of surface element above ground level. Please
3828	!-- note, height of surface element is determined with respect to
3829	!-- its height above ground of the reference grid point in atmosphere,
3830	!-- Therefore, substract the offset values when assessing the terrain
3831	!-- height.
3832	IF ( terrain_height_f%from_file ) THEN
3833	z_agl = zw(k) - terrain_height_f%var(j-surf_usm_v(l)%joff, &
3834	i-surf_usm_v(l)%ioff)
3835	ELSE
3836	z_agl = zw(k)
3837	ENDIF
3838	!
3839	!-- Set flag for ground level
3840	IF ( z_agl <= ground_floor_level_l ) &
3841	surf_usm_v(l)%ground_level(m) = .TRUE.
3842
3843	ENDDO
3844	ENDDO
3845	!
3846	!-- Initialization of resistances.
3847	DO m = 1, surf_usm_h%ns
3848	surf_usm_h%r_a(m) = 50.0_wp
3849	surf_usm_h%r_a_green(m) = 50.0_wp
3850	surf_usm_h%r_a_window(m) = 50.0_wp
3851	ENDDO
3852	DO l = 0, 3
3853	DO m = 1, surf_usm_v(l)%ns
3854	surf_usm_v(l)%r_a(m) = 50.0_wp
3855	surf_usm_v(l)%r_a_green(m) = 50.0_wp
3856	surf_usm_v(l)%r_a_window(m) = 50.0_wp
3857	ENDDO
3858	ENDDO
3859
3860	!
3861	!-- Map values onto horizontal elemements
3862	DO m = 1, surf_usm_h%ns
3863	surf_usm_h%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
3864	surf_usm_h%g_d(m) = 0.0_wp !< canopy_resistance_coefficient
3865	ENDDO
3866	!
3867	!-- Map values onto vertical elements, even though this does not make
3868	!-- much sense.
3869	DO l = 0, 3
3870	DO m = 1, surf_usm_v(l)%ns
3871	surf_usm_v(l)%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
3872	surf_usm_v(l)%g_d(m) = 0.0_wp !< canopy_resistance_coefficient
3873	ENDDO
3874	ENDDO
3875
3876	!
3877	!-- Initialize urban-type surface attribute. According to initialization in
3878	!-- land-surface model, follow a 3-level approach.
3879	!-- Level 1 - initialization via default attributes
3880	DO m = 1, surf_usm_h%ns
3881	!
3882	!-- Now, all horizontal surfaces are roof surfaces (?)
3883	surf_usm_h%isroof_surf(m) = .TRUE.
3884	surf_usm_h%surface_types(m) = roof_category !< default category for root surface
3885	!
3886	!-- In order to distinguish between ground floor level and
3887	!-- above-ground-floor level surfaces, set input indices.
3888
3889	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3890	surf_usm_h%ground_level(m) )
3891	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
3892	surf_usm_h%ground_level(m) )
3893	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
3894	surf_usm_h%ground_level(m) )
3895	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
3896	surf_usm_h%ground_level(m) )
3897	!
3898	!-- Store building type and its name on each surface element
3899	surf_usm_h%building_type(m) = building_type
3900	surf_usm_h%building_type_name(m) = building_type_name(building_type)
3901	!
3902	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
3903	surf_usm_h%frac(ind_veg_wall,m) = building_pars(ind_wall_frac_r,building_type)
3904	surf_usm_h%frac(ind_pav_green,m) = building_pars(ind_green_frac_r,building_type)
3905	surf_usm_h%frac(ind_wat_win,m) = building_pars(ind_win_frac_r,building_type)
3906	surf_usm_h%lai(m) = building_pars(ind_lai_r,building_type)
3907
3908	surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1_wall_r,building_type)
3909	surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1_wall_r,building_type)
3910	surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2_wall_r,building_type)
3911	surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3_wall_r,building_type)
3912	surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1_wall_r,building_type)
3913	surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1_wall_r,building_type)
3914	surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2_wall_r,building_type)
3915	surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3_wall_r,building_type)
3916	surf_usm_h%rho_c_green(nzb_wall,m) = rho_c_soil !building_pars(ind_hc1_wall_r,building_type)
3917	surf_usm_h%rho_c_green(nzb_wall+1,m) = rho_c_soil !building_pars(ind_hc1_wall_r,building_type)
3918	surf_usm_h%rho_c_green(nzb_wall+2,m) = rho_c_soil !building_pars(ind_hc2_wall_r,building_type)
3919	surf_usm_h%rho_c_green(nzb_wall+3,m) = rho_c_soil !building_pars(ind_hc3_wall_r,building_type)
3920	surf_usm_h%lambda_h_green(nzb_wall,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,building_type)
3921	surf_usm_h%lambda_h_green(nzb_wall+1,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,building_type)
3922	surf_usm_h%lambda_h_green(nzb_wall+2,m) = lambda_h_green_sm !building_pars(ind_tc2_wall_r,building_type)
3923	surf_usm_h%lambda_h_green(nzb_wall+3,m) = lambda_h_green_sm !building_pars(ind_tc3_wall_r,building_type)
3924	surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1_win_r,building_type)
3925	surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1_win_r,building_type)
3926	surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2_win_r,building_type)
3927	surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3_win_r,building_type)
3928	surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1_win_r,building_type)
3929	surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1_win_r,building_type)
3930	surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2_win_r,building_type)
3931	surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3_win_r,building_type)
3932
3933	surf_usm_h%target_temp_summer(m) = building_pars(ind_indoor_target_temp_summer,building_type)
3934	surf_usm_h%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,building_type)
3935	!
3936	!-- emissivity of wall-, green- and window fraction
3937	surf_usm_h%emissivity(ind_veg_wall,m) = building_pars(ind_emis_wall_r,building_type)
3938	surf_usm_h%emissivity(ind_pav_green,m) = building_pars(ind_emis_green_r,building_type)
3939	surf_usm_h%emissivity(ind_wat_win,m) = building_pars(ind_emis_win_r,building_type)
3940
3941	surf_usm_h%transmissivity(m) = building_pars(ind_trans_r,building_type)
3942
3943	surf_usm_h%z0(m) = building_pars(ind_z0,building_type)
3944	surf_usm_h%z0h(m) = building_pars(ind_z0qh,building_type)
3945	surf_usm_h%z0q(m) = building_pars(ind_z0qh,building_type)
3946	!
3947	!-- albedo type for wall fraction, green fraction, window fraction
3948	surf_usm_h%albedo_type(ind_veg_wall,m) = INT( building_pars(ind_alb_wall_r,building_type) )
3949	surf_usm_h%albedo_type(ind_pav_green,m) = INT( building_pars(ind_alb_green_r,building_type) )
3950	surf_usm_h%albedo_type(ind_wat_win,m) = INT( building_pars(ind_alb_win_r,building_type) )
3951
3952	surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1_wall_r,building_type)
3953	surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2_wall_r,building_type)
3954	surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3_wall_r,building_type)
3955	surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4_wall_r,building_type)
3956
3957	surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1_wall_r,building_type)
3958	surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2_wall_r,building_type)
3959	surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3_wall_r,building_type)
3960	surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4_wall_r,building_type)
3961
3962	surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1_win_r,building_type)
3963	surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2_win_r,building_type)
3964	surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3_win_r,building_type)
3965	surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4_win_r,building_type)
3966
3967	surf_usm_h%c_surface(m) = building_pars(ind_c_surface,building_type)
3968	surf_usm_h%lambda_surf(m) = building_pars(ind_lambda_surf,building_type)
3969	surf_usm_h%c_surface_green(m) = building_pars(ind_c_surface_green,building_type)
3970	surf_usm_h%lambda_surf_green(m) = building_pars(ind_lambda_surf_green,building_type)
3971	surf_usm_h%c_surface_window(m) = building_pars(ind_c_surface_win,building_type)
3972	surf_usm_h%lambda_surf_window(m) = building_pars(ind_lambda_surf_win,building_type)
3973
3974	surf_usm_h%green_type_roof(m) = building_pars(ind_green_type_roof,building_type)
3975
3976	ENDDO
3977
3978	DO l = 0, 3
3979	DO m = 1, surf_usm_v(l)%ns
3980
3981	surf_usm_v(l)%surface_types(m) = wall_category !< default category for root surface
3982	!
3983	!-- In order to distinguish between ground floor level and
3984	!-- above-ground-floor level surfaces, set input indices.
3985	ind_alb_green = MERGE( ind_alb_green_gfl, ind_alb_green_agfl, &
3986	surf_usm_v(l)%ground_level(m) )
3987	ind_alb_wall = MERGE( ind_alb_wall_gfl, ind_alb_wall_agfl, &
3988	surf_usm_v(l)%ground_level(m) )
3989	ind_alb_win = MERGE( ind_alb_win_gfl, ind_alb_win_agfl, &
3990	surf_usm_v(l)%ground_level(m) )
3991	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
3992	surf_usm_v(l)%ground_level(m) )
3993	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
3994	surf_usm_v(l)%ground_level(m) )
3995	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
3996	surf_usm_v(l)%ground_level(m) )
3997	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
3998	surf_usm_v(l)%ground_level(m) )
3999	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
4000	surf_usm_v(l)%ground_level(m) )
4001	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
4002	surf_usm_v(l)%ground_level(m) )
4003	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
4004	surf_usm_v(l)%ground_level(m) )
4005	ind_hc1_win = MERGE( ind_hc1_win_gfl, ind_hc1_win_agfl, &
4006	surf_usm_v(l)%ground_level(m) )
4007	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
4008	surf_usm_v(l)%ground_level(m) )
4009	ind_hc2_win = MERGE( ind_hc2_win_gfl, ind_hc2_win_agfl, &
4010	surf_usm_v(l)%ground_level(m) )
4011	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
4012	surf_usm_v(l)%ground_level(m) )
4013	ind_hc3_win = MERGE( ind_hc3_win_gfl, ind_hc3_win_agfl, &
4014	surf_usm_v(l)%ground_level(m) )
4015	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
4016	surf_usm_v(l)%ground_level(m) )
4017	ind_tc1_win = MERGE( ind_tc1_win_gfl, ind_tc1_win_agfl, &
4018	surf_usm_v(l)%ground_level(m) )
4019	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
4020	surf_usm_v(l)%ground_level(m) )
4021	ind_tc2_win = MERGE( ind_tc2_win_gfl, ind_tc2_win_agfl, &
4022	surf_usm_v(l)%ground_level(m) )
4023	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
4024	surf_usm_v(l)%ground_level(m) )
4025	ind_tc3_win = MERGE( ind_tc3_win_gfl, ind_tc3_win_agfl, &
4026	surf_usm_v(l)%ground_level(m) )
4027	ind_thick_1 = MERGE( ind_thick_1_gfl, ind_thick_1_agfl, &
4028	surf_usm_v(l)%ground_level(m) )
4029	ind_thick_1_win = MERGE( ind_thick_1_win_gfl, ind_thick_1_win_agfl, &
4030	surf_usm_v(l)%ground_level(m) )
4031	ind_thick_2 = MERGE( ind_thick_2_gfl, ind_thick_2_agfl, &
4032	surf_usm_v(l)%ground_level(m) )
4033	ind_thick_2_win = MERGE( ind_thick_2_win_gfl, ind_thick_2_win_agfl, &
4034	surf_usm_v(l)%ground_level(m) )
4035	ind_thick_3 = MERGE( ind_thick_3_gfl, ind_thick_3_agfl, &
4036	surf_usm_v(l)%ground_level(m) )
4037	ind_thick_3_win = MERGE( ind_thick_3_win_gfl, ind_thick_3_win_agfl, &
4038	surf_usm_v(l)%ground_level(m) )
4039	ind_thick_4 = MERGE( ind_thick_4_gfl, ind_thick_4_agfl, &
4040	surf_usm_v(l)%ground_level(m) )
4041	ind_thick_4_win = MERGE( ind_thick_4_win_gfl, ind_thick_4_win_agfl, &
4042	surf_usm_v(l)%ground_level(m) )
4043	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
4044	surf_usm_v(l)%ground_level(m) )
4045	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
4046	surf_usm_v(l)%ground_level(m) )
4047	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
4048	surf_usm_v(l)%ground_level(m) )
4049	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
4050	surf_usm_v(l)%ground_level(m) )
4051	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
4052	surf_usm_v(l)%ground_level(m) )
4053	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
4054	surf_usm_v(l)%ground_level(m) )
4055	!
4056	!-- Store building type and its name on each surface element
4057	surf_usm_v(l)%building_type(m) = building_type
4058	surf_usm_v(l)%building_type_name(m) = building_type_name(building_type)
4059	!
4060	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
4061	surf_usm_v(l)%frac(ind_veg_wall,m) = building_pars(ind_wall_frac,building_type)
4062	surf_usm_v(l)%frac(ind_pav_green,m) = building_pars(ind_green_frac_w,building_type)
4063	surf_usm_v(l)%frac(ind_wat_win,m) = building_pars(ind_win_frac,building_type)
4064	surf_usm_v(l)%lai(m) = building_pars(ind_lai_w,building_type)
4065
4066	surf_usm_v(l)%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,building_type)
4067	surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,building_type)
4068	surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,building_type)
4069	surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,building_type)
4070
4071	surf_usm_v(l)%rho_c_green(nzb_wall,m) = rho_c_soil !building_pars(ind_hc1,building_type)
4072	surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = rho_c_soil !building_pars(ind_hc1,building_type)
4073	surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = rho_c_soil !building_pars(ind_hc2,building_type)
4074	surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = rho_c_soil !building_pars(ind_hc3,building_type)
4075
4076	surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars(ind_hc1_win,building_type)
4077	surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1_win,building_type)
4078	surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2_win,building_type)
4079	surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3_win,building_type)
4080
4081	surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars(ind_tc1,building_type)
4082	surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,building_type)
4083	surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,building_type)
4084	surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,building_type)
4085
4086	surf_usm_v(l)%lambda_h_green(nzb_wall,m) = lambda_h_green_sm !building_pars(ind_tc1,building_type)
4087	surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = lambda_h_green_sm !building_pars(ind_tc1,building_type)
4088	surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = lambda_h_green_sm !building_pars(ind_tc2,building_type)
4089	surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = lambda_h_green_sm !building_pars(ind_tc3,building_type)
4090
4091	surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1_win,building_type)
4092	surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1_win,building_type)
4093	surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2_win,building_type)
4094	surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3_win,building_type)
4095
4096	surf_usm_v(l)%target_temp_summer(m) = building_pars(ind_indoor_target_temp_summer,building_type)
4097	surf_usm_v(l)%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,building_type)
4098	!
4099	!-- emissivity of wall-, green- and window fraction
4100	surf_usm_v(l)%emissivity(ind_veg_wall,m) = building_pars(ind_emis_wall,building_type)
4101	surf_usm_v(l)%emissivity(ind_pav_green,m) = building_pars(ind_emis_green,building_type)
4102	surf_usm_v(l)%emissivity(ind_wat_win,m) = building_pars(ind_emis_win,building_type)
4103
4104	surf_usm_v(l)%transmissivity(m) = building_pars(ind_trans,building_type)
4105
4106	surf_usm_v(l)%z0(m) = building_pars(ind_z0,building_type)
4107	surf_usm_v(l)%z0h(m) = building_pars(ind_z0qh,building_type)
4108	surf_usm_v(l)%z0q(m) = building_pars(ind_z0qh,building_type)
4109
4110	surf_usm_v(l)%albedo_type(ind_veg_wall,m) = INT( building_pars(ind_alb_wall,building_type) )
4111	surf_usm_v(l)%albedo_type(ind_pav_green,m) = INT( building_pars(ind_alb_green,building_type) )
4112	surf_usm_v(l)%albedo_type(ind_wat_win,m) = INT( building_pars(ind_alb_win,building_type) )
4113
4114	surf_usm_v(l)%zw(nzb_wall,m) = building_pars(ind_thick_1,building_type)
4115	surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
4116	surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
4117	surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
4118
4119	surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars(ind_thick_1,building_type)
4120	surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,building_type)
4121	surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,building_type)
4122	surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,building_type)
4123
4124	surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars(ind_thick_1_win,building_type)
4125	surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2_win,building_type)
4126	surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3_win,building_type)
4127	surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4_win,building_type)
4128
4129	surf_usm_v(l)%c_surface(m) = building_pars(ind_c_surface,building_type)
4130	surf_usm_v(l)%lambda_surf(m) = building_pars(ind_lambda_surf,building_type)
4131	surf_usm_v(l)%c_surface_green(m) = building_pars(ind_c_surface_green,building_type)
4132	surf_usm_v(l)%lambda_surf_green(m) = building_pars(ind_lambda_surf_green,building_type)
4133	surf_usm_v(l)%c_surface_window(m) = building_pars(ind_c_surface_win,building_type)
4134	surf_usm_v(l)%lambda_surf_window(m) = building_pars(ind_lambda_surf_win,building_type)
4135
4136	ENDDO
4137	ENDDO
4138	!
4139	!-- Level 2 - initialization via building type read from file
4140	IF ( building_type_f%from_file ) THEN
4141	DO m = 1, surf_usm_h%ns
4142	i = surf_usm_h%i(m)
4143	j = surf_usm_h%j(m)
4144	!
4145	!-- For the moment, limit building type to 6 (to overcome errors in input file).
4146	st = building_type_f%var(j,i)
4147	IF ( st /= building_type_f%fill ) THEN
4148
4149	!
4150	!-- In order to distinguish between ground floor level and
4151	!-- above-ground-floor level surfaces, set input indices.
4152
4153	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
4154	surf_usm_h%ground_level(m) )
4155	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
4156	surf_usm_h%ground_level(m) )
4157	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
4158	surf_usm_h%ground_level(m) )
4159	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
4160	surf_usm_h%ground_level(m) )
4161	!
4162	!-- Store building type and its name on each surface element
4163	surf_usm_h%building_type(m) = st
4164	surf_usm_h%building_type_name(m) = building_type_name(st)
4165	!
4166	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
4167	surf_usm_h%frac(ind_veg_wall,m) = building_pars(ind_wall_frac_r,st)
4168	surf_usm_h%frac(ind_pav_green,m) = building_pars(ind_green_frac_r,st)
4169	surf_usm_h%frac(ind_wat_win,m) = building_pars(ind_win_frac_r,st)
4170	surf_usm_h%lai(m) = building_pars(ind_lai_r,st)
4171
4172	surf_usm_h%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1_wall_r,st)
4173	surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1_wall_r,st)
4174	surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2_wall_r,st)
4175	surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3_wall_r,st)
4176	surf_usm_h%lambda_h(nzb_wall,m) = building_pars(ind_tc1_wall_r,st)
4177	surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1_wall_r,st)
4178	surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2_wall_r,st)
4179	surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3_wall_r,st)
4180
4181	surf_usm_h%rho_c_green(nzb_wall,m) = rho_c_soil !building_pars(ind_hc1_wall_r,st)
4182	surf_usm_h%rho_c_green(nzb_wall+1,m) = rho_c_soil !building_pars(ind_hc1_wall_r,st)
4183	surf_usm_h%rho_c_green(nzb_wall+2,m) = rho_c_soil !building_pars(ind_hc2_wall_r,st)
4184	surf_usm_h%rho_c_green(nzb_wall+3,m) = rho_c_soil !building_pars(ind_hc3_wall_r,st)
4185	surf_usm_h%lambda_h_green(nzb_wall,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,st)
4186	surf_usm_h%lambda_h_green(nzb_wall+1,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,st)
4187	surf_usm_h%lambda_h_green(nzb_wall+2,m) = lambda_h_green_sm !building_pars(ind_tc2_wall_r,st)
4188	surf_usm_h%lambda_h_green(nzb_wall+3,m) = lambda_h_green_sm !building_pars(ind_tc3_wall_r,st)
4189
4190	surf_usm_h%rho_c_window(nzb_wall,m) = building_pars(ind_hc1_win_r,st)
4191	surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1_win_r,st)
4192	surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2_win_r,st)
4193	surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3_win_r,st)
4194	surf_usm_h%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1_win_r,st)
4195	surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1_win_r,st)
4196	surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2_win_r,st)
4197	surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3_win_r,st)
4198
4199	surf_usm_h%target_temp_summer(m) = building_pars(ind_indoor_target_temp_summer,st)
4200	surf_usm_h%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,st)
4201	!
4202	!-- emissivity of wall-, green- and window fraction
4203	surf_usm_h%emissivity(ind_veg_wall,m) = building_pars(ind_emis_wall_r,st)
4204	surf_usm_h%emissivity(ind_pav_green,m) = building_pars(ind_emis_green_r,st)
4205	surf_usm_h%emissivity(ind_wat_win,m) = building_pars(ind_emis_win_r,st)
4206
4207	surf_usm_h%transmissivity(m) = building_pars(ind_trans_r,st)
4208
4209	surf_usm_h%z0(m) = building_pars(ind_z0,st)
4210	surf_usm_h%z0h(m) = building_pars(ind_z0qh,st)
4211	surf_usm_h%z0q(m) = building_pars(ind_z0qh,st)
4212	!
4213	!-- albedo type for wall fraction, green fraction, window fraction
4214	surf_usm_h%albedo_type(ind_veg_wall,m) = INT( building_pars(ind_alb_wall_r,st) )
4215	surf_usm_h%albedo_type(ind_pav_green,m) = INT( building_pars(ind_alb_green_r,st) )
4216	surf_usm_h%albedo_type(ind_wat_win,m) = INT( building_pars(ind_alb_win_r,st) )
4217
4218	surf_usm_h%zw(nzb_wall,m) = building_pars(ind_thick_1_wall_r,st)
4219	surf_usm_h%zw(nzb_wall+1,m) = building_pars(ind_thick_2_wall_r,st)
4220	surf_usm_h%zw(nzb_wall+2,m) = building_pars(ind_thick_3_wall_r,st)
4221	surf_usm_h%zw(nzb_wall+3,m) = building_pars(ind_thick_4_wall_r,st)
4222
4223	surf_usm_h%zw_green(nzb_wall,m) = building_pars(ind_thick_1_wall_r,st)
4224	surf_usm_h%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2_wall_r,st)
4225	surf_usm_h%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3_wall_r,st)
4226	surf_usm_h%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4_wall_r,st)
4227
4228	surf_usm_h%zw_window(nzb_wall,m) = building_pars(ind_thick_1_win_r,st)
4229	surf_usm_h%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2_win_r,st)
4230	surf_usm_h%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3_win_r,st)
4231	surf_usm_h%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4_win_r,st)
4232
4233	surf_usm_h%c_surface(m) = building_pars(ind_c_surface,st)
4234	surf_usm_h%lambda_surf(m) = building_pars(ind_lambda_surf,st)
4235	surf_usm_h%c_surface_green(m) = building_pars(ind_c_surface_green,st)
4236	surf_usm_h%lambda_surf_green(m) = building_pars(ind_lambda_surf_green,st)
4237	surf_usm_h%c_surface_window(m) = building_pars(ind_c_surface_win,st)
4238	surf_usm_h%lambda_surf_window(m) = building_pars(ind_lambda_surf_win,st)
4239
4240	surf_usm_h%green_type_roof(m) = building_pars(ind_green_type_roof,st)
4241
4242	ENDIF
4243	ENDDO
4244
4245	DO l = 0, 3
4246	DO m = 1, surf_usm_v(l)%ns
4247	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
4248	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
4249	!
4250	!-- For the moment, limit building type to 6 (to overcome errors in input file).
4251
4252	st = building_type_f%var(j,i)
4253	IF ( st /= building_type_f%fill ) THEN
4254
4255	!
4256	!-- In order to distinguish between ground floor level and
4257	!-- above-ground-floor level surfaces, set input indices.
4258	ind_alb_green = MERGE( ind_alb_green_gfl, ind_alb_green_agfl, &
4259	surf_usm_v(l)%ground_level(m) )
4260	ind_alb_wall = MERGE( ind_alb_wall_gfl, ind_alb_wall_agfl, &
4261	surf_usm_v(l)%ground_level(m) )
4262	ind_alb_win = MERGE( ind_alb_win_gfl, ind_alb_win_agfl, &
4263	surf_usm_v(l)%ground_level(m) )
4264	ind_wall_frac = MERGE( ind_wall_frac_gfl, ind_wall_frac_agfl, &
4265	surf_usm_v(l)%ground_level(m) )
4266	ind_win_frac = MERGE( ind_win_frac_gfl, ind_win_frac_agfl, &
4267	surf_usm_v(l)%ground_level(m) )
4268	ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
4269	surf_usm_v(l)%ground_level(m) )
4270	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
4271	surf_usm_v(l)%ground_level(m) )
4272	ind_lai_r = MERGE( ind_lai_r_gfl, ind_lai_r_agfl, &
4273	surf_usm_v(l)%ground_level(m) )
4274	ind_lai_w = MERGE( ind_lai_w_gfl, ind_lai_w_agfl, &
4275	surf_usm_v(l)%ground_level(m) )
4276	ind_hc1 = MERGE( ind_hc1_gfl, ind_hc1_agfl, &
4277	surf_usm_v(l)%ground_level(m) )
4278	ind_hc1_win = MERGE( ind_hc1_win_gfl, ind_hc1_win_agfl, &
4279	surf_usm_v(l)%ground_level(m) )
4280	ind_hc2 = MERGE( ind_hc2_gfl, ind_hc2_agfl, &
4281	surf_usm_v(l)%ground_level(m) )
4282	ind_hc2_win = MERGE( ind_hc2_win_gfl, ind_hc2_win_agfl, &
4283	surf_usm_v(l)%ground_level(m) )
4284	ind_hc3 = MERGE( ind_hc3_gfl, ind_hc3_agfl, &
4285	surf_usm_v(l)%ground_level(m) )
4286	ind_hc3_win = MERGE( ind_hc3_win_gfl, ind_hc3_win_agfl, &
4287	surf_usm_v(l)%ground_level(m) )
4288	ind_tc1 = MERGE( ind_tc1_gfl, ind_tc1_agfl, &
4289	surf_usm_v(l)%ground_level(m) )
4290	ind_tc1_win = MERGE( ind_tc1_win_gfl, ind_tc1_win_agfl, &
4291	surf_usm_v(l)%ground_level(m) )
4292	ind_tc2 = MERGE( ind_tc2_gfl, ind_tc2_agfl, &
4293	surf_usm_v(l)%ground_level(m) )
4294	ind_tc2_win = MERGE( ind_tc2_win_gfl, ind_tc2_win_agfl, &
4295	surf_usm_v(l)%ground_level(m) )
4296	ind_tc3 = MERGE( ind_tc3_gfl, ind_tc3_agfl, &
4297	surf_usm_v(l)%ground_level(m) )
4298	ind_tc3_win = MERGE( ind_tc3_win_gfl, ind_tc3_win_agfl, &
4299	surf_usm_v(l)%ground_level(m) )
4300	ind_thick_1 = MERGE( ind_thick_1_gfl, ind_thick_1_agfl, &
4301	surf_usm_v(l)%ground_level(m) )
4302	ind_thick_1_win = MERGE( ind_thick_1_win_gfl, ind_thick_1_win_agfl, &
4303	surf_usm_v(l)%ground_level(m) )
4304	ind_thick_2 = MERGE( ind_thick_2_gfl, ind_thick_2_agfl, &
4305	surf_usm_v(l)%ground_level(m) )
4306	ind_thick_2_win = MERGE( ind_thick_2_win_gfl, ind_thick_2_win_agfl, &
4307	surf_usm_v(l)%ground_level(m) )
4308	ind_thick_3 = MERGE( ind_thick_3_gfl, ind_thick_3_agfl, &
4309	surf_usm_v(l)%ground_level(m) )
4310	ind_thick_3_win = MERGE( ind_thick_3_win_gfl, ind_thick_3_win_agfl, &
4311	surf_usm_v(l)%ground_level(m) )
4312	ind_thick_4 = MERGE( ind_thick_4_gfl, ind_thick_4_agfl, &
4313	surf_usm_v(l)%ground_level(m) )
4314	ind_thick_4_win = MERGE( ind_thick_4_win_gfl, ind_thick_4_win_agfl, &
4315	surf_usm_v(l)%ground_level(m) )
4316	ind_emis_wall = MERGE( ind_emis_wall_gfl, ind_emis_wall_agfl, &
4317	surf_usm_v(l)%ground_level(m) )
4318	ind_emis_green = MERGE( ind_emis_green_gfl, ind_emis_green_agfl, &
4319	surf_usm_v(l)%ground_level(m) )
4320	ind_emis_win = MERGE( ind_emis_win_gfl, ind_emis_win_agfl, &
4321	surf_usm_v(l)%ground_level(m) )
4322	ind_trans = MERGE( ind_trans_gfl, ind_trans_agfl, &
4323	surf_usm_v(l)%ground_level(m) )
4324	ind_z0 = MERGE( ind_z0_gfl, ind_z0_agfl, &
4325	surf_usm_v(l)%ground_level(m) )
4326	ind_z0qh = MERGE( ind_z0qh_gfl, ind_z0qh_agfl, &
4327	surf_usm_v(l)%ground_level(m) )
4328	!
4329	!-- Store building type and its name on each surface element
4330	surf_usm_v(l)%building_type(m) = st
4331	surf_usm_v(l)%building_type_name(m) = building_type_name(st)
4332	!
4333	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
4334	surf_usm_v(l)%frac(ind_veg_wall,m) = building_pars(ind_wall_frac,st)
4335	surf_usm_v(l)%frac(ind_pav_green,m) = building_pars(ind_green_frac_w,st)
4336	surf_usm_v(l)%frac(ind_wat_win,m) = building_pars(ind_win_frac,st)
4337	surf_usm_v(l)%lai(m) = building_pars(ind_lai_w,st)
4338
4339	surf_usm_v(l)%rho_c_wall(nzb_wall,m) = building_pars(ind_hc1,st)
4340	surf_usm_v(l)%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1,st)
4341	surf_usm_v(l)%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2,st)
4342	surf_usm_v(l)%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3,st)
4343
4344	surf_usm_v(l)%rho_c_green(nzb_wall,m) = rho_c_soil !building_pars(ind_hc1,st)
4345	surf_usm_v(l)%rho_c_green(nzb_wall+1,m) = rho_c_soil !building_pars(ind_hc1,st)
4346	surf_usm_v(l)%rho_c_green(nzb_wall+2,m) = rho_c_soil !building_pars(ind_hc2,st)
4347	surf_usm_v(l)%rho_c_green(nzb_wall+3,m) = rho_c_soil !building_pars(ind_hc3,st)
4348
4349	surf_usm_v(l)%rho_c_window(nzb_wall,m) = building_pars(ind_hc1_win,st)
4350	surf_usm_v(l)%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1_win,st)
4351	surf_usm_v(l)%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2_win,st)
4352	surf_usm_v(l)%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3_win,st)
4353
4354	surf_usm_v(l)%lambda_h(nzb_wall,m) = building_pars(ind_tc1,st)
4355	surf_usm_v(l)%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1,st)
4356	surf_usm_v(l)%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2,st)
4357	surf_usm_v(l)%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3,st)
4358
4359	surf_usm_v(l)%lambda_h_green(nzb_wall,m) = lambda_h_green_sm !building_pars(ind_tc1,st)
4360	surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) = lambda_h_green_sm !building_pars(ind_tc1,st)
4361	surf_usm_v(l)%lambda_h_green(nzb_wall+2,m) = lambda_h_green_sm !building_pars(ind_tc2,st)
4362	surf_usm_v(l)%lambda_h_green(nzb_wall+3,m) = lambda_h_green_sm !building_pars(ind_tc3,st)
4363
4364	surf_usm_v(l)%lambda_h_window(nzb_wall,m) = building_pars(ind_tc1_win,st)
4365	surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1_win,st)
4366	surf_usm_v(l)%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2_win,st)
4367	surf_usm_v(l)%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3_win,st)
4368
4369	surf_usm_v(l)%target_temp_summer(m) = building_pars(ind_indoor_target_temp_summer,st)
4370	surf_usm_v(l)%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,st)
4371	!
4372	!-- emissivity of wall-, green- and window fraction
4373	surf_usm_v(l)%emissivity(ind_veg_wall,m) = building_pars(ind_emis_wall,st)
4374	surf_usm_v(l)%emissivity(ind_pav_green,m) = building_pars(ind_emis_green,st)
4375	surf_usm_v(l)%emissivity(ind_wat_win,m) = building_pars(ind_emis_win,st)
4376
4377	surf_usm_v(l)%transmissivity(m) = building_pars(ind_trans,st)
4378
4379	surf_usm_v(l)%z0(m) = building_pars(ind_z0,st)
4380	surf_usm_v(l)%z0h(m) = building_pars(ind_z0qh,st)
4381	surf_usm_v(l)%z0q(m) = building_pars(ind_z0qh,st)
4382
4383	surf_usm_v(l)%albedo_type(ind_veg_wall,m) = INT( building_pars(ind_alb_wall,st) )
4384	surf_usm_v(l)%albedo_type(ind_pav_green,m) = INT( building_pars(ind_alb_green,st) )
4385	surf_usm_v(l)%albedo_type(ind_wat_win,m) = INT( building_pars(ind_alb_win,st) )
4386
4387	surf_usm_v(l)%zw(nzb_wall,m) = building_pars(ind_thick_1,st)
4388	surf_usm_v(l)%zw(nzb_wall+1,m) = building_pars(ind_thick_2,st)
4389	surf_usm_v(l)%zw(nzb_wall+2,m) = building_pars(ind_thick_3,st)
4390	surf_usm_v(l)%zw(nzb_wall+3,m) = building_pars(ind_thick_4,st)
4391
4392	surf_usm_v(l)%zw_green(nzb_wall,m) = building_pars(ind_thick_1,st)
4393	surf_usm_v(l)%zw_green(nzb_wall+1,m) = building_pars(ind_thick_2,st)
4394	surf_usm_v(l)%zw_green(nzb_wall+2,m) = building_pars(ind_thick_3,st)
4395	surf_usm_v(l)%zw_green(nzb_wall+3,m) = building_pars(ind_thick_4,st)
4396
4397	surf_usm_v(l)%zw_window(nzb_wall,m) = building_pars(ind_thick_1_win,st)
4398	surf_usm_v(l)%zw_window(nzb_wall+1,m) = building_pars(ind_thick_2_win,st)
4399	surf_usm_v(l)%zw_window(nzb_wall+2,m) = building_pars(ind_thick_3_win,st)
4400	surf_usm_v(l)%zw_window(nzb_wall+3,m) = building_pars(ind_thick_4_win,st)
4401
4402	surf_usm_v(l)%c_surface(m) = building_pars(ind_c_surface,st)
4403	surf_usm_v(l)%lambda_surf(m) = building_pars(ind_lambda_surf,st)
4404	surf_usm_v(l)%c_surface_green(m) = building_pars(ind_c_surface_green,st)
4405	surf_usm_v(l)%lambda_surf_green(m) = building_pars(ind_lambda_surf_green,st)
4406	surf_usm_v(l)%c_surface_window(m) = building_pars(ind_c_surface_win,st)
4407	surf_usm_v(l)%lambda_surf_window(m) = building_pars(ind_lambda_surf_win,st)
4408
4409
4410	ENDIF
4411	ENDDO
4412	ENDDO
4413	ENDIF
4414
4415	!
4416	!-- Level 3 - initialization via building_pars read from file. Note, only
4417	!-- variables that are also defined in the input-standard can be initialized
4418	!-- via file. Other variables will be initialized on level 1 or 2.
4419	IF ( building_pars_f%from_file ) THEN
4420	DO m = 1, surf_usm_h%ns
4421	i = surf_usm_h%i(m)
4422	j = surf_usm_h%j(m)
4423
4424	!
4425	!-- In order to distinguish between ground floor level and
4426	!-- above-ground-floor level surfaces, set input indices.
4427	ind_wall_frac = MERGE( ind_wall_frac_gfl, &
4428	ind_wall_frac_agfl, &
4429	surf_usm_h%ground_level(m) )
4430	ind_green_frac_r = MERGE( ind_green_frac_r_gfl, &
4431	ind_green_frac_r_agfl, &
4432	surf_usm_h%ground_level(m) )
4433	ind_win_frac = MERGE( ind_win_frac_gfl, &
4434	ind_win_frac_agfl, &
4435	surf_usm_h%ground_level(m) )
4436	ind_lai_r = MERGE( ind_lai_r_gfl, &
4437	ind_lai_r_agfl, &
4438	surf_usm_h%ground_level(m) )
4439	ind_z0 = MERGE( ind_z0_gfl, &
4440	ind_z0_agfl, &
4441	surf_usm_h%ground_level(m) )
4442	ind_z0qh = MERGE( ind_z0qh_gfl, &
4443	ind_z0qh_agfl, &
4444	surf_usm_h%ground_level(m) )
4445	ind_hc1 = MERGE( ind_hc1_gfl, &
4446	ind_hc1_agfl, &
4447	surf_usm_h%ground_level(m) )
4448	ind_hc2 = MERGE( ind_hc2_gfl, &
4449	ind_hc2_agfl, &
4450	surf_usm_h%ground_level(m) )
4451	ind_hc3 = MERGE( ind_hc3_gfl, &
4452	ind_hc3_agfl, &
4453	surf_usm_h%ground_level(m) )
4454	ind_tc1 = MERGE( ind_tc1_gfl, &
4455	ind_tc1_agfl, &
4456	surf_usm_h%ground_level(m) )
4457	ind_tc2 = MERGE( ind_tc2_gfl, &
4458	ind_tc2_agfl, &
4459	surf_usm_h%ground_level(m) )
4460	ind_tc3 = MERGE( ind_tc3_gfl, &
4461	ind_tc3_agfl, &
4462	surf_usm_h%ground_level(m) )
4463	ind_emis_wall = MERGE( ind_emis_wall_gfl, &
4464	ind_emis_wall_agfl, &
4465	surf_usm_h%ground_level(m) )
4466	ind_emis_green = MERGE( ind_emis_green_gfl, &
4467	ind_emis_green_agfl, &
4468	surf_usm_h%ground_level(m) )
4469	ind_emis_win = MERGE( ind_emis_win_gfl, &
4470	ind_emis_win_agfl, &
4471	surf_usm_h%ground_level(m) )
4472	ind_trans = MERGE( ind_trans_gfl, &
4473	ind_trans_agfl, &
4474	surf_usm_h%ground_level(m) )
4475
4476	!
4477	!-- Initialize relatvie wall- (0), green- (1) and window (2) fractions
4478	IF ( building_pars_f%pars_xy(ind_wall_frac,j,i) /= &
4479	building_pars_f%fill ) &
4480	surf_usm_h%frac(ind_veg_wall,m) = &
4481	building_pars_f%pars_xy(ind_wall_frac,j,i)
4482
4483	IF ( building_pars_f%pars_xy(ind_green_frac_r,j,i) /= &
4484	building_pars_f%fill ) &
4485	surf_usm_h%frac(ind_pav_green,m) = &
4486	building_pars_f%pars_xy(ind_green_frac_r,j,i)
4487
4488	IF ( building_pars_f%pars_xy(ind_win_frac,j,i) /= &
4489	building_pars_f%fill ) &
4490	surf_usm_h%frac(ind_wat_win,m) = &
4491	building_pars_f%pars_xy(ind_win_frac,j,i)
4492
4493	IF ( building_pars_f%pars_xy(ind_lai_r,j,i) /= &
4494	building_pars_f%fill ) &
4495	surf_usm_h%lai(m) = building_pars_f%pars_xy(ind_lai_r,j,i)
4496
4497	IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= &
4498	building_pars_f%fill ) THEN
4499	surf_usm_h%rho_c_wall(nzb_wall,m) = &
4500	building_pars_f%pars_xy(ind_hc1,j,i)
4501	surf_usm_h%rho_c_wall(nzb_wall+1,m) = &
4502	building_pars_f%pars_xy(ind_hc1,j,i)
4503	ENDIF
4504
4505
4506	IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= &
4507	building_pars_f%fill ) &
4508	surf_usm_h%rho_c_wall(nzb_wall+2,m) = &
4509	building_pars_f%pars_xy(ind_hc2,j,i)
4510
4511	IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= &
4512	building_pars_f%fill ) &
4513	surf_usm_h%rho_c_wall(nzb_wall+3,m) = &
4514	building_pars_f%pars_xy(ind_hc3,j,i)
4515
4516	IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= &
4517	building_pars_f%fill ) THEN
4518	surf_usm_h%rho_c_green(nzb_wall,m) = &
4519	building_pars_f%pars_xy(ind_hc1,j,i)
4520	surf_usm_h%rho_c_green(nzb_wall+1,m) = &
4521	building_pars_f%pars_xy(ind_hc1,j,i)
4522	ENDIF
4523	IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= &
4524	building_pars_f%fill ) &
4525	surf_usm_h%rho_c_green(nzb_wall+2,m) = &
4526	building_pars_f%pars_xy(ind_hc2,j,i)
4527
4528	IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= &
4529	building_pars_f%fill ) &
4530	surf_usm_h%rho_c_green(nzb_wall+3,m) = &
4531	building_pars_f%pars_xy(ind_hc3,j,i)
4532
4533	IF ( building_pars_f%pars_xy(ind_hc1,j,i) /= &
4534	building_pars_f%fill ) THEN
4535	surf_usm_h%rho_c_window(nzb_wall,m) = &
4536	building_pars_f%pars_xy(ind_hc1,j,i)
4537	surf_usm_h%rho_c_window(nzb_wall+1,m) = &
4538	building_pars_f%pars_xy(ind_hc1,j,i)
4539	ENDIF
4540	IF ( building_pars_f%pars_xy(ind_hc2,j,i) /= &
4541	building_pars_f%fill ) &
4542	surf_usm_h%rho_c_window(nzb_wall+2,m) = &
4543	building_pars_f%pars_xy(ind_hc2,j,i)
4544
4545	IF ( building_pars_f%pars_xy(ind_hc3,j,i) /= &
4546	building_pars_f%fill ) &
4547	surf_usm_h%rho_c_window(nzb_wall+3,m) = &
4548	building_pars_f%pars_xy(ind_hc3,j,i)
4549
4550	IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= &
4551	building_pars_f%fill ) THEN
4552	surf_usm_h%lambda_h(nzb_wall,m) = &
4553	building_pars_f%pars_xy(ind_tc1,j,i)
4554	surf_usm_h%lambda_h(nzb_wall+1,m) = &
4555	building_pars_f%pars_xy(ind_tc1,j,i)
4556	ENDIF
4557	IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= &
4558	building_pars_f%fill ) &
4559	surf_usm_h%lambda_h(nzb_wall+2,m) = &
4560	building_pars_f%pars_xy(ind_tc2,j,i)
4561
4562	IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= &
4563	building_pars_f%fill ) &
4564	surf_usm_h%lambda_h(nzb_wall+3,m) = &
4565	building_pars_f%pars_xy(ind_tc3,j,i)
4566
4567	IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= &
4568	building_pars_f%fill ) THEN
4569	surf_usm_h%lambda_h_green(nzb_wall,m) = &
4570	building_pars_f%pars_xy(ind_tc1,j,i)
4571	surf_usm_h%lambda_h_green(nzb_wall+1,m) = &
4572	building_pars_f%pars_xy(ind_tc1,j,i)
4573	ENDIF
4574	IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= &
4575	building_pars_f%fill ) &
4576	surf_usm_h%lambda_h_green(nzb_wall+2,m) = &
4577	building_pars_f%pars_xy(ind_tc2,j,i)
4578
4579	IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= &
4580	building_pars_f%fill ) &
4581	surf_usm_h%lambda_h_green(nzb_wall+3,m) = &
4582	building_pars_f%pars_xy(ind_tc3,j,i)
4583
4584	IF ( building_pars_f%pars_xy(ind_tc1,j,i) /= &
4585	building_pars_f%fill ) THEN
4586	surf_usm_h%lambda_h_window(nzb_wall,m) = &
4587	building_pars_f%pars_xy(ind_tc1,j,i)
4588	surf_usm_h%lambda_h_window(nzb_wall+1,m) = &
4589	building_pars_f%pars_xy(ind_tc1,j,i)
4590	ENDIF
4591	IF ( building_pars_f%pars_xy(ind_tc2,j,i) /= &
4592	building_pars_f%fill ) &
4593	surf_usm_h%lambda_h_window(nzb_wall+2,m) = &
4594	building_pars_f%pars_xy(ind_tc2,j,i)
4595
4596	IF ( building_pars_f%pars_xy(ind_tc3,j,i) /= &
4597	building_pars_f%fill ) &
4598	surf_usm_h%lambda_h_window(nzb_wall+3,m) = &
4599	building_pars_f%pars_xy(ind_tc3,j,i)
4600
4601	IF ( building_pars_f%pars_xy(ind_indoor_target_temp_summer,j,i) /=&
4602	building_pars_f%fill ) &
4603	surf_usm_h%target_temp_summer(m) = &
4604	building_pars_f%pars_xy(ind_indoor_target_temp_summer,j,i)
4605	IF ( building_pars_f%pars_xy(ind_indoor_target_temp_winter,j,i) /=&
4606	building_pars_f%fill ) &
4607	surf_usm_h%target_temp_winter(m) = &
4608	building_pars_f%pars_xy(ind_indoor_target_temp_winter,j,i)
4609
4610	�