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

Last change on this file since 3151 was 3151, checked in by raasch, 3 years ago

output of PALM code revision number in palmrun header, some preprocessor define strings removed, small revisions in configuration files

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