Home

Context Navigation

urban_surface_mod.f90 @ 2943

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