source: palm/trunk/SOURCE/urban_surface_mod.f90 @ 3347

Last change on this file since 3347 was 3347, checked in by suehring, 3 years ago

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

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