Home

Context Navigation

urban_surface_mod.f90 @ 2894

Last change on this file since 2894 was 2894, checked in by Giersch, 6 years ago
Reading/Writing? data in case of restart runs revised
Property svn:keywords set to `Id`
File size: 413.5 KB

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