Home

Context Navigation

urban_surface_mod.f90 @ 3378

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