Home

Context Navigation

urban_surface_mod.f90 @ 3196

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