Home

Context Navigation

urban_surface_mod.f90 @ 3382

Last change on this file since 3382 was 3382, checked in by knoop, 6 years ago
Bugix: made array declaration in urban_surface_mod Fortran Standard conform
Property svn:keywords set to `Id`
File size: 446.4 KB

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