Home

Context Navigation

urban_surface_mod.f90 @ 2977

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