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

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

Limit aerodynamic resistance at vertial building walls; check for roughness not exceeding surface-layer height

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