Home

Context Navigation

urban_surface_mod.f90 @ 2921

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