Home

Context Navigation

urban_surface_mod.f90 @ 2746

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