Home

Context Navigation

urban_surface_mod.f90 @ 4050

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