Home

Context Navigation

urban_surface_mod.f90 @ 2696

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