Home

Context Navigation

urban_surface_mod.f90 @ 4051

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