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

!> @file urban_surface_mod.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 2015-2018 Czech Technical University in Prague
! Copyright 2015-2018 Institute of Computer Science of the 
!                     Czech Academy of Sciences, Prague
! Copyright 1997-2018 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: urban_surface_mod.f90 3636 2018-12-19 13:48:34Z raasch $
! nopointer option removed
! 
! 3614 2018-12-10 07:05:46Z raasch
! unused variables removed
! 
! 3607 2018-12-07 11:56:58Z suehring
! Output of radiation-related quantities migrated to radiation_model_mod.
! 
! 3597 2018-12-04 08:40:18Z maronga
! Fixed calculation method of near surface air potential temperature at 10 cm
! and moved to surface_layer_fluxes. Removed unnecessary _eb strings.
! 
! 3524 2018-11-14 13:36:44Z raasch
! bugfix concerning allocation of t_surf_wall_v
! 
! 3502 2018-11-07 14:45:23Z suehring
! Disable initialization of building roofs with ground-floor-level properties, 
! since this causes strong oscillations of surface temperature during the 
! spinup.
! 
! 3469 2018-10-30 20:05:07Z kanani
! Add missing PUBLIC variables for new indoor model
! 
! 3449 2018-10-29 19:36:56Z suehring
! Bugfix: Fix average arrays allocations in usm_average_3d_data (J.Resler)
! Bugfix: Fix reading wall temperatures (J.Resler)
! Bugfix: Fix treating of outputs for wall temperature and sky view factors (J.Resler)
!
! 
! 3435 2018-10-26 18:25:44Z gronemeier
! Bugfix: allocate gamma_w_green_sat until nzt_wall+1
! 
! 3418 2018-10-24 16:07:39Z kanani
! (rvtils, srissman)
! -Updated building databse, two green roof types (ind_green_type_roof)
! -Latent heat flux for green walls and roofs, new output of latent heatflux
!  and soil water content of green roof substrate
! -t_surf changed to t_surf_wall
! -Added namelist parameter usm_wall_mod for lower wall tendency
!  of first two wall layers during spinup
! -Window calculations deactivated during spinup
! 
! 3382 2018-10-19 13:10:32Z knoop
! Bugix: made array declaration Fortran Standard conform
! 
! 3378 2018-10-19 12:34:59Z kanani
! merge from radiation branch (r3362) into trunk
! (moh.hefny):
! - check the requested output variables if they are correct
! - added unscheduled_radiation_calls switch to control force_radiation_call
! - minor formate changes
! 
! 3371 2018-10-18 13:40:12Z knoop
! Set flag indicating that albedo at urban surfaces is already initialized
! 
! 3347 2018-10-15 14:21:08Z suehring
! Enable USM initialization with default building parameters in case no static
! input file exist. 
! 
! 3343 2018-10-15 10:38:52Z suehring
! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
! 
! 3274 2018-09-24 15:42:55Z knoop
! Modularization of all bulk cloud physics code components
! 
! 3248 2018-09-14 09:42:06Z sward
! Minor formating changes
! 
! 3246 2018-09-13 15:14:50Z sward
! Added error handling for input namelist via parin_fail_message
! 
! 3241 2018-09-12 15:02:00Z raasch
! unused variables removed
! 
! 3223 2018-08-30 13:48:17Z suehring
! Bugfix for commit 3222
! 
! 3222 2018-08-30 13:35:35Z suehring
! Introduction of surface array for type and its name
! 
! 3203 2018-08-23 10:48:36Z suehring
! Revise bulk parameter for emissivity at ground-floor level
! 
! 3196 2018-08-13 12:26:14Z maronga
! Added maximum aerodynamic resistance of 300 for horiztonal surfaces.
! 
! 3176 2018-07-26 17:12:48Z suehring
! Bugfix, update virtual potential surface temparture, else heat fluxes on 
! roofs might become unphysical
! 
! 3152 2018-07-19 13:26:52Z suehring
! Initialize q_surface, which might be used in surface_layer_fluxes
! 
! 3151 2018-07-19 08:45:38Z raasch
! remaining preprocessor define strings __check removed
! 
! 3136 2018-07-16 14:48:21Z suehring
! Limit also roughness length for heat and moisture where necessary
! 
! 3123 2018-07-12 16:21:53Z suehring
! Correct working precision for INTEGER number
! 
! 3115 2018-07-10 12:49:26Z suehring
! Additional building type to represent bridges 
! 
! 3091 2018-06-28 16:20:35Z suehring
! - Limit aerodynamic resistance at vertical walls.
! - Add check for local roughness length not exceeding surface-layer height and
!   limit roughness length where necessary.
! 
! 3065 2018-06-12 07:03:02Z Giersch
! Unused array dxdir was removed, dz was replaced by dzu to consider vertical
! grid stretching
! 
! 3049 2018-05-29 13:52:36Z Giersch
! Error messages revised
! 
! 3045 2018-05-28 07:55:41Z Giersch
! Error message added
! 
! 3029 2018-05-23 12:19:17Z raasch
! bugfix: close unit 151 instead of 90
! 
! 3014 2018-05-09 08:42:38Z maronga
! Added pc_transpiration_rate
! 
! 2977 2018-04-17 10:27:57Z kanani
! Implement changes from branch radiation (r2948-2971) with minor modifications.
! (moh.hefny):
! Extended exn for all model domain height to avoid the need to get nzut.
! 
! 2963 2018-04-12 14:47:44Z suehring
! Introduce index for vegetation/wall, pavement/green-wall and water/window 
! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
! 
! 2943 2018-04-03 16:17:10Z suehring
! Calculate exner function at all height levels and remove some un-used
! variables.
! 
! 2932 2018-03-26 09:39:22Z maronga
! renamed urban_surface_par to urban_surface_parameters
! 
! 2921 2018-03-22 15:05:23Z Giersch
! The activation of spinup has been moved to parin
! 
! 2920 2018-03-22 11:22:01Z kanani
! Remove unused pcbl, npcbl from ONLY list
! moh.hefny:
! Fixed bugs introduced by new structures and by moving radiation interaction
! into radiation_model_mod.f90.
! Bugfix: usm data output 3D didn't respect directions
! 
! 2906 2018-03-19 08:56:40Z Giersch
! Local variable ids has to be initialized with a value of -1 in 
! usm_average_3d_data
! 
! 2894 2018-03-15 09:17:58Z Giersch
! Calculations of the index range of the subdomain on file which overlaps with
! the current subdomain are already done in read_restart_data_mod,
! usm_read/write_restart_data have been renamed to usm_r/wrd_local, variable 
! named found has been introduced for checking if restart data was found, 
! reading of restart strings has been moved completely to 
! read_restart_data_mod, usm_rrd_local is already inside the overlap loop 
! programmed in read_restart_data_mod, SAVE attribute added where necessary, 
! deallocation and allocation of some arrays have been changed to take care of 
! different restart files that can be opened (index i), the marker *** end usm 
! *** is not necessary anymore, strings and their respective lengths are 
! written out and read now in case of restart runs to get rid of prescribed 
! character lengths
! 
! 2805 2018-02-14 17:00:09Z suehring
! Initialization of resistances.
! 
! 2797 2018-02-08 13:24:35Z suehring
! Comment concerning output of ground-heat flux added.
! 
! 2766 2018-01-22 17:17:47Z kanani
! Removed redundant commas, added some blanks
! 
! 2765 2018-01-22 11:34:58Z maronga
! Major bugfix in calculation of f_shf. Adjustment of roughness lengths in
! building_pars
! 
! 2750 2018-01-15 16:26:51Z knoop
! Move flag plant canopy to modules
! 
! 2737 2018-01-11 14:58:11Z kanani
! Removed unused variables t_surf_whole...
! 
! 2735 2018-01-11 12:01:27Z suehring
! resistances are saved in surface attributes
! 
! 2723 2018-01-05 09:27:03Z maronga
! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
! 
! 2720 2018-01-02 16:27:15Z kanani
! Correction of comment
! 
! 2718 2018-01-02 08:49:38Z maronga
! Corrected "Former revisions" section
! 
! 2705 2017-12-18 11:26:23Z maronga
! Changes from last commit documented
! 
! 2703 2017-12-15 20:12:38Z maronga
! Workaround for calculation of r_a
!
! 2696 2017-12-14 17:12:51Z kanani
! - Change in file header (GPL part)
! - Bugfix in calculation of pt_surface and related fluxes. (BM)
! - Do not write surface temperatures onto pt array as this might cause 
!   problems with nesting. (MS)
! - Revised calculation of pt1 (now done in surface_layer_fluxes).
!   Bugfix, f_shf_window and f_shf_green were not set at vertical surface 
!   elements. (MS)
! - merged with branch ebsolver
!   green building surfaces do not evaporate yet
!   properties of green wall layers and window layers are taken from wall layers
!   this input data is missing. (RvT)
! - Merged with branch radiation (developed by Mohamed Salim)
! - Revised initialization. (MS)
! - Rename emiss_surf into emissivity, roughness_wall into z0, albedo_surf into 
!   albedo. (MS)
! - Move first call of usm_radiatin from usm_init to init_3d_model
! - fixed problem with near surface temperature
! - added near surface temperature pt_10cm_h(m), pt_10cm_v(l)%t(m)
! - does not work with temp profile including stability, ol
!   pt_10cm = pt1 now
! - merged with 2357 bugfix, error message for nopointer version
! - added indoor model coupling with wall heat flux
! - added green substrate/ dry vegetation layer for buildings
! - merged with 2232 new surface-type structure
! - added transmissivity of window tiles
! - added MOSAIK tile approach for 3 different surfaces (RvT)
! 
! 2583 2017-10-26 13:58:38Z knoop
! Bugfix: reverted MPI_Win_allocate_cptr introduction in last commit 
!
! 2582 2017-10-26 13:19:46Z hellstea
! Workaround for gnufortran compiler added in usm_calc_svf. CALL MPI_Win_allocate is
! replaced by CALL MPI_Win_allocate_cptr if defined ( __gnufortran ).
!
! 2544 2017-10-13 18:09:32Z maronga
! Date and time quantities are now read from date_and_time_mod. Solar constant is
! read from radiation_model_mod
! 
! 2516 2017-10-04 11:03:04Z suehring
! Remove tabs
! 
! 2514 2017-10-04 09:52:37Z suehring
! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
! no output of ghost layer data
! 
! 2350 2017-08-15 11:48:26Z kanani
! Bugfix and error message for nopointer version. 
! Additional "! defined(__nopointer)" as workaround to enable compilation of 
! nopointer version.
! 
! 2318 2017-07-20 17:27:44Z suehring
! Get topography top index via Function call 
! 
! 2317 2017-07-20 17:27:19Z suehring
! Bugfix: adjust output of shf. Added support for spinups
! 
! 2287 2017-06-15 16:46:30Z suehring
! Bugfix in determination topography-top index
! 
! 2269 2017-06-09 11:57:32Z suehring
! Enable restart runs with different number of PEs
! Bugfixes nopointer branch
! 
! 2258 2017-06-08 07:55:13Z suehring
! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
! 
! 2233 2017-05-30 18:08:54Z suehring
!
! 2232 2017-05-30 17:47:52Z suehring
! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux; 
! insteat, heat fluxes are directly applied in diffusion_s.
! 
! 2213 2017-04-24 15:10:35Z kanani
! Removal of output quantities usm_lad and usm_canopy_hr
! 
! 2209 2017-04-19 09:34:46Z kanani
! cpp switch __mpi3 removed,
! minor formatting,
! small bugfix for division by zero (Krc)
! 
! 2113 2017-01-12 13:40:46Z kanani
! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
! 
! 2071 2016-11-17 11:22:14Z maronga
! Small bugfix (Resler)
! 
! 2031 2016-10-21 15:11:58Z knoop
! renamed variable rho to rho_ocean
! 
! 2024 2016-10-12 16:42:37Z kanani
! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
! optimization of MPI-RMA operations,
! declaration of pcbl as integer,
! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
! use of new control parameter varnamelength,
! added output variables usm_rad_ressw, usm_rad_reslw,
! minor formatting changes,
! minor optimizations.
! 
! 2011 2016-09-19 17:29:57Z kanani
! Major reformatting according to PALM coding standard (comments, blanks,
! alphabetical ordering, etc.),
! removed debug_prints, 
! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
! defined in MODULE control_parameters (modules.f90) to avoid circular 
! dependencies,
! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
! 
! 2007 2016-08-24 15:47:17Z kanani
! Initial revision
!
!
! Description:
! ------------
! 2016/6/9 - Initial version of the USM (Urban Surface Model)
!            authors: Jaroslav Resler, Pavel Krc
!                     (Czech Technical University in Prague and Institute of
!                      Computer Science of the Czech Academy of Sciences, Prague)
!            with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
!            partly inspired by PALM LSM (B. Maronga)
!            parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
!> Module for Urban Surface Model (USM)
!> The module includes:
!>    1. radiation model with direct/diffuse radiation, shading, reflections
!>       and integration with plant canopy
!>    2. wall and wall surface model
!>    3. surface layer energy balance
!>    4. anthropogenic heat (only from transportation so far)
!>    5. necessary auxiliary subroutines (reading inputs, writing outputs,
!>       restart simulations, ...)
!> It also make use of standard radiation and integrates it into
!> urban surface model.
!>
!> Further work:
!> -------------
!> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity
!>    from surfaces that are visible from local surfaces (i.e. there is a SVF
!>    where target is local). To do that, radiosity will be exchanged after each
!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
!>
!> 2. Temporarily large values of surface heat flux can be observed, up to
!>    1.2 Km/s, which seem to be not realistic.
!>
!> @todo Output of _av variables in case of restarts
!> @todo Revise flux conversion in energy-balance solver
!> @todo Bugfixing in nopointer branch 
!> @todo Check optimizations for RMA operations
!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
!>       factor 3 between min and max time
!> @todo Move setting of flag indoor_model to indoor_model_mod once available
!> @todo Check divisions in wtend (etc.) calculations for possible division
!>       by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
!> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
!> @todo Move plant canopy stuff into plant canopy code
!------------------------------------------------------------------------------!
 MODULE urban_surface_mod

    USE arrays_3d,                                                             &
        ONLY:  hyp, zu, pt, p, u, v, w, tend, exner, hyrho, prr, q, ql, vpt

    USE calc_mean_profile_mod,                                                 &
        ONLY:  calc_mean_profile

    USE basic_constants_and_equations_mod,                                     &
        ONLY:  c_p, g, kappa, pi, r_d, rho_l, l_v

    USE control_parameters,                                                    &
        ONLY:  coupling_start_time, topography, dt_3d, humidity,               &
               intermediate_timestep_count, initializing_actions,              &
               intermediate_timestep_count_max, simulated_time, end_time,      &
               timestep_scheme, tsc, coupling_char, io_blocks, io_group,       &
               message_string, time_since_reference_point, surface_pressure,   &
               pt_surface, large_scale_forcing, lsf_surf, spinup,              &
               spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
               average_count_3d, varnamelength, urban_surface,                 &
               plant_canopy, dz

    USE bulk_cloud_model_mod,                                                  &
        ONLY: bulk_cloud_model, precipitation
               
    USE cpulog,                                                                &
        ONLY:  cpu_log, log_point, log_point_s

    USE date_and_time_mod,                                                     &
        ONLY:  time_utc_init

    USE grid_variables,                                                        &
        ONLY:  dx, dy, ddx, ddy, ddx2, ddy2

    USE indices,                                                               &
        ONLY:  nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys,    &
               nysg, nzb, nzt, nbgp, wall_flags_0

    USE, INTRINSIC :: iso_c_binding 

    USE kinds
              
    USE pegrid
    
    USE plant_canopy_model_mod,                                                &
        ONLY:  pc_heating_rate, pc_transpiration_rate, pc_latent_rate
    
    USE radiation_model_mod,                                                   &
        ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,        &
               radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
               sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon,              &
               force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u,       &
               iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id,       &
               iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
               idcsf, ndcsf, kdcsf, pct,                                       &
               nzub, nzut, unscheduled_radiation_calls

    USE statistics,                                                            &
        ONLY:  hom, statistic_regions

    USE surface_mod,                                                           &
        ONLY:  get_topography_top_index_ji, get_topography_top_index,          &
               ind_pav_green, ind_veg_wall, ind_wat_win, surf_usm_h,           &
               surf_usm_v, surface_restore_elements


    IMPLICIT NONE

!
!-- USM model constants

    REAL(wp), PARAMETER ::                     &
              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
              lambda_h_green_dry = 0.19_wp,    & ! heat conductivity for dry soil   
              lambda_h_green_sm  = 3.44_wp,    & ! heat conductivity of the soil matrix
              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
              rho_c_water        = 4.20E6_wp !,  & ! volumetric heat capacity of water
!               m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)

!
!-- Soil parameters I           alpha_vg,      l_vg_green,    n_vg, gamma_w_green_sat
    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
                                 /), (/ 4, 7 /) )

!
!-- Soil parameters II              swc_sat,     fc,   wilt,    swc_res  
    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
                                 /), (/ 4, 7 /) )
                                
    !   value 9999999.9_wp -> generic available or user-defined value must be set
!   otherwise -> no generic variable and user setting is optional
    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
                field_capacity = 9999999.9_wp,          & !< NAMELIST fc
                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_green_sat
                lambda_h_green_sat = 0.0_wp,            & !< heat conductivity for saturated soil
                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
                saturation_moisture = 9999999.9_wp,     & !< NAMELIST m_sat
                wilting_point = 9999999.9_wp!,           & !< NAMELIST m_wilt
    

!-- configuration parameters (they can be setup in PALM config)
    LOGICAL ::  usm_material_model = .TRUE.        !< flag parameter indicating wheather the  model of heat in materials is used
    LOGICAL ::  usm_anthropogenic_heat = .FALSE.   !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used
    LOGICAL ::  force_radiation_call_l = .FALSE.   !< flag parameter for unscheduled radiation model calls
    LOGICAL ::  indoor_model = .FALSE.             !< whether to use the indoor model
    LOGICAL ::  read_wall_temp_3d = .FALSE.
    LOGICAL ::  usm_wall_mod = .FALSE.             !< reduces conductivity of the first 2 wall layers by factor 0.1


    INTEGER(iwp) ::  building_type = 1               !< default building type (preleminary setting)
    INTEGER(iwp) ::  land_category = 2               !< default category for land surface
    INTEGER(iwp) ::  wall_category = 2               !< default category for wall surface over pedestrian zone
    INTEGER(iwp) ::  pedestrian_category = 2         !< default category for wall surface in pedestrian zone
    INTEGER(iwp) ::  roof_category = 2               !< default category for root surface
    REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
!
!-- Indices of input attributes for (above) ground floor level 
    INTEGER(iwp) ::  ind_alb_wall_agfl     = 65  !< index in input list for albedo_type of wall above ground floor level
    INTEGER(iwp) ::  ind_alb_wall_gfl      = 32  !< index in input list for albedo_type of wall ground floor level
    INTEGER(iwp) ::  ind_alb_wall_r        = 96  !< index in input list for albedo_type of wall roof
    INTEGER(iwp) ::  ind_alb_green_agfl    = 83  !< index in input list for albedo_type of green above ground floor level
    INTEGER(iwp) ::  ind_alb_green_gfl     = 50  !< index in input list for albedo_type of green ground floor level
    INTEGER(iwp) ::  ind_alb_green_r       = 115 !< index in input list for albedo_type of green roof
    INTEGER(iwp) ::  ind_alb_win_agfl      = 79  !< index in input list for albedo_type of window fraction above ground floor level
    INTEGER(iwp) ::  ind_alb_win_gfl       = 46  !< index in input list for albedo_type of window fraction ground floor level 
    INTEGER(iwp) ::  ind_alb_win_r         = 110 !< index in input list for albedo_type of window fraction roof
    INTEGER(iwp) ::  ind_emis_wall_agfl    = 64  !< index in input list for wall emissivity, above ground floor level
    INTEGER(iwp) ::  ind_emis_wall_gfl     = 31  !< index in input list for wall emissivity, ground floor level
    INTEGER(iwp) ::  ind_emis_wall_r       = 95  !< index in input list for wall emissivity, roof
    INTEGER(iwp) ::  ind_emis_green_agfl   = 82  !< index in input list for green emissivity, above ground floor level
    INTEGER(iwp) ::  ind_emis_green_gfl    = 49  !< index in input list for green emissivity, ground floor level
    INTEGER(iwp) ::  ind_emis_green_r      = 114 !< index in input list for green emissivity, roof
    INTEGER(iwp) ::  ind_emis_win_agfl     = 77  !< index in input list for window emissivity, above ground floor level
    INTEGER(iwp) ::  ind_emis_win_gfl      = 44  !< index in input list for window emissivity, ground floor level
    INTEGER(iwp) ::  ind_emis_win_r        = 108 !< index in input list for window emissivity, roof
    INTEGER(iwp) ::  ind_green_frac_w_agfl = 80  !< index in input list for green fraction on wall, above ground floor level
    INTEGER(iwp) ::  ind_green_frac_w_gfl  = 47  !< index in input list for green fraction on wall, ground floor level
    INTEGER(iwp) ::  ind_green_frac_r_agfl = 112 !< index in input list for green fraction on roof, above ground floor level
    INTEGER(iwp) ::  ind_green_frac_r_gfl  = 111 !< index in input list for green fraction on roof, ground floor level
    INTEGER(iwp) ::  ind_hc1_agfl          = 58  !< index in input list for heat capacity at first wall layer, above ground floor level
    INTEGER(iwp) ::  ind_hc1_gfl           = 25  !< index in input list for heat capacity at first wall layer, ground floor level
    INTEGER(iwp) ::  ind_hc1_wall_r        = 89  !< index in input list for heat capacity at first wall layer, roof
    INTEGER(iwp) ::  ind_hc1_win_agfl      = 71  !< index in input list for heat capacity at first window layer, above ground floor level
    INTEGER(iwp) ::  ind_hc1_win_gfl       = 38  !< index in input list for heat capacity at first window layer, ground floor level
    INTEGER(iwp) ::  ind_hc1_win_r         = 102 !< index in input list for heat capacity at first window layer, roof
    INTEGER(iwp) ::  ind_hc2_agfl          = 59  !< index in input list for heat capacity at second wall layer, above ground floor level
    INTEGER(iwp) ::  ind_hc2_gfl           = 26  !< index in input list for heat capacity at second wall layer, ground floor level
    INTEGER(iwp) ::  ind_hc2_wall_r        = 90  !< index in input list for heat capacity at second wall layer, roof
    INTEGER(iwp) ::  ind_hc2_win_agfl      = 72  !< index in input list for heat capacity at second window layer, above ground floor level
    INTEGER(iwp) ::  ind_hc2_win_gfl       = 39  !< index in input list for heat capacity at second window layer, ground floor level
    INTEGER(iwp) ::  ind_hc2_win_r         = 103 !< index in input list for heat capacity at second window layer, roof
    INTEGER(iwp) ::  ind_hc3_agfl          = 60  !< index in input list for heat capacity at third wall layer, above ground floor level
    INTEGER(iwp) ::  ind_hc3_gfl           = 27  !< index in input list for heat capacity at third wall layer, ground floor level
    INTEGER(iwp) ::  ind_hc3_wall_r        = 91  !< index in input list for heat capacity at third wall layer, roof
    INTEGER(iwp) ::  ind_hc3_win_agfl      = 73  !< index in input list for heat capacity at third window layer, above ground floor level
    INTEGER(iwp) ::  ind_hc3_win_gfl       = 40  !< index in input list for heat capacity at third window layer, ground floor level
    INTEGER(iwp) ::  ind_hc3_win_r         = 104 !< index in input list for heat capacity at third window layer, roof
    INTEGER(iwp) ::  ind_gflh              = 17  !< index in input list for ground floor level height
    INTEGER(iwp) ::  ind_lai_r_agfl        = 113 !< index in input list for LAI on roof, above ground floor level
    INTEGER(iwp) ::  ind_lai_r_gfl         = 113  !< index in input list for LAI on roof, ground floor level
    INTEGER(iwp) ::  ind_lai_w_agfl        = 81  !< index in input list for LAI on wall, above ground floor level
    INTEGER(iwp) ::  ind_lai_w_gfl         = 48  !< index in input list for LAI on wall, ground floor level
    INTEGER(iwp) ::  ind_tc1_agfl          = 61  !< index in input list for thermal conductivity at first wall layer, above ground floor level
    INTEGER(iwp) ::  ind_tc1_gfl           = 28  !< index in input list for thermal conductivity at first wall layer, ground floor level
    INTEGER(iwp) ::  ind_tc1_wall_r        = 92  !< index in input list for thermal conductivity at first wall layer, roof
    INTEGER(iwp) ::  ind_tc1_win_agfl      = 74  !< index in input list for thermal conductivity at first window layer, above ground floor level
    INTEGER(iwp) ::  ind_tc1_win_gfl       = 41  !< index in input list for thermal conductivity at first window layer, ground floor level
    INTEGER(iwp) ::  ind_tc1_win_r         = 105 !< index in input list for thermal conductivity at first window layer, roof
    INTEGER(iwp) ::  ind_tc2_agfl          = 62  !< index in input list for thermal conductivity at second wall layer, above ground floor level
    INTEGER(iwp) ::  ind_tc2_gfl           = 29  !< index in input list for thermal conductivity at second wall layer, ground floor level
    INTEGER(iwp) ::  ind_tc2_wall_r        = 93  !< index in input list for thermal conductivity at second wall layer, roof
    INTEGER(iwp) ::  ind_tc2_win_agfl      = 75  !< index in input list for thermal conductivity at second window layer, above ground floor level
    INTEGER(iwp) ::  ind_tc2_win_gfl       = 42  !< index in input list for thermal conductivity at second window layer, ground floor level
    INTEGER(iwp) ::  ind_tc2_win_r         = 106 !< index in input list for thermal conductivity at second window layer, ground floor level
    INTEGER(iwp) ::  ind_tc3_agfl          = 63  !< index in input list for thermal conductivity at third wall layer, above ground floor level
    INTEGER(iwp) ::  ind_tc3_gfl           = 30  !< index in input list for thermal conductivity at third wall layer, ground floor level
    INTEGER(iwp) ::  ind_tc3_wall_r        = 94  !< index in input list for thermal conductivity at third wall layer, roof
    INTEGER(iwp) ::  ind_tc3_win_agfl      = 76  !< index in input list for thermal conductivity at third window layer, above ground floor level
    INTEGER(iwp) ::  ind_tc3_win_gfl       = 43  !< index in input list for thermal conductivity at third window layer, ground floor level
    INTEGER(iwp) ::  ind_tc3_win_r         = 107 !< index in input list for thermal conductivity at third window layer, roof
    INTEGER(iwp) ::  ind_thick_1_agfl      = 54  !< index for wall layer thickness - 1st layer above ground floor level
    INTEGER(iwp) ::  ind_thick_1_gfl       = 21  !< index for wall layer thickness - 1st layer ground floor level
    INTEGER(iwp) ::  ind_thick_1_wall_r    = 85  !< index for wall layer thickness - 1st layer roof
    INTEGER(iwp) ::  ind_thick_1_win_agfl  = 67  !< index for window layer thickness - 1st layer above ground floor level
    INTEGER(iwp) ::  ind_thick_1_win_gfl   = 34  !< index for window layer thickness - 1st layer ground floor level
    INTEGER(iwp) ::  ind_thick_1_win_r     = 98  !< index for window layer thickness - 1st layer roof
    INTEGER(iwp) ::  ind_thick_2_agfl      = 55  !< index for wall layer thickness - 2nd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_2_gfl       = 22  !< index for wall layer thickness - 2nd layer ground floor level
    INTEGER(iwp) ::  ind_thick_2_wall_r    = 86  !< index for wall layer thickness - 2nd layer roof
    INTEGER(iwp) ::  ind_thick_2_win_agfl  = 68  !< index for window layer thickness - 2nd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_2_win_gfl   = 35  !< index for window layer thickness - 2nd layer ground floor level
    INTEGER(iwp) ::  ind_thick_2_win_r     = 99  !< index for window layer thickness - 2nd layer roof
    INTEGER(iwp) ::  ind_thick_3_agfl      = 56  !< index for wall layer thickness - 3rd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_3_gfl       = 23  !< index for wall layer thickness - 3rd layer ground floor level
    INTEGER(iwp) ::  ind_thick_3_wall_r    = 87  !< index for wall layer thickness - 3rd layer roof
    INTEGER(iwp) ::  ind_thick_3_win_agfl  = 69  !< index for window layer thickness - 3rd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_3_win_gfl   = 36  !< index for window layer thickness - 3rd layer ground floor level  
    INTEGER(iwp) ::  ind_thick_3_win_r     = 100 !< index for window layer thickness - 3rd layer roof
    INTEGER(iwp) ::  ind_thick_4_agfl      = 57  !< index for wall layer thickness - 4th layer above ground floor level
    INTEGER(iwp) ::  ind_thick_4_gfl       = 24  !< index for wall layer thickness - 4th layer ground floor level
    INTEGER(iwp) ::  ind_thick_4_wall_r    = 88  !< index for wall layer thickness - 4st layer roof
    INTEGER(iwp) ::  ind_thick_4_win_agfl  = 70  !< index for window layer thickness - 4th layer above ground floor level
    INTEGER(iwp) ::  ind_thick_4_win_gfl   = 37  !< index for window layer thickness - 4th layer ground floor level
    INTEGER(iwp) ::  ind_thick_4_win_r     = 101 !< index for window layer thickness - 4th layer roof
    INTEGER(iwp) ::  ind_trans_agfl        = 78  !< index in input list for window transmissivity, above ground floor level
    INTEGER(iwp) ::  ind_trans_gfl         = 45  !< index in input list for window transmissivity, ground floor level
    INTEGER(iwp) ::  ind_trans_r           = 109 !< index in input list for window transmissivity, roof
    INTEGER(iwp) ::  ind_wall_frac_agfl    = 53  !< index in input list for wall fraction, above ground floor level
    INTEGER(iwp) ::  ind_wall_frac_gfl     = 20  !< index in input list for wall fraction, ground floor level
    INTEGER(iwp) ::  ind_wall_frac_r       = 84  !< index in input list for wall fraction, roof
    INTEGER(iwp) ::  ind_win_frac_agfl     = 66  !< index in input list for window fraction, above ground floor level
    INTEGER(iwp) ::  ind_win_frac_gfl      = 33  !< index in input list for window fraction, ground floor level
    INTEGER(iwp) ::  ind_win_frac_r        = 97  !< index in input list for window fraction, roof
    INTEGER(iwp) ::  ind_z0_agfl           = 51  !< index in input list for z0, above ground floor level
    INTEGER(iwp) ::  ind_z0_gfl            = 18  !< index in input list for z0, ground floor level
    INTEGER(iwp) ::  ind_z0qh_agfl         = 52  !< index in input list for z0h / z0q, above ground floor level
    INTEGER(iwp) ::  ind_z0qh_gfl          = 19  !< index in input list for z0h / z0q, ground floor level
    INTEGER(iwp) ::  ind_green_type_roof   = 116 !< index in input list for type of green roof


    REAL(wp)  ::  roof_height_limit = 4.0_wp          !< height for distinguish between land surfaces and roofs
    REAL(wp)  ::  ground_floor_level = 4.0_wp        !< default ground floor level


    CHARACTER(37), DIMENSION(0:7), PARAMETER :: building_type_name = (/     &
                                   'user-defined                         ', & !  0 
                                   'residential - 1950                   ', & !  1 
                                   'residential 1951 - 2000              ', & !  2
                                   'residential 2001 -                   ', & !  3
                                   'office - 1950                        ', & !  4 
                                   'office 1951 - 2000                   ', & !  5
                                   'office 2001 -                        ', & !  6
                                   'bridges                              '  & !  7
                                                                     /)
!
!-- building parameters, 6 different types
!-- Parameter for urban surface model
!-- 0 - heat capacity wall surface, 1 - heat capacity of window surface, 2 - heat capacity of green surface
!-- 3 - thermal conductivity of wall surface, 4 - thermal conductivity of window surface, 
!-- 5 - thermal conductivty of green surface, 6 - wall fraction ground plate, 
!-- 7 - 1st wall layer thickness ground plate, 8 - 2nd wall layer thickness ground plate
!-- 9 - 3rd wall layer thickness ground plate, 10 - 4th wall layer thickness ground plate,
!-- 11 - heat capacity 1st/2nd wall layer ground plate, 12 - heat capacity 3rd wall layer ground plate
!-- 13 - heat capacity 4th wall layer ground plate, 14 - thermal conductivity 1st/2nd wall layer ground plate,
!-- 15 - thermal conductivity 3rd wall layer ground plate, 16 - thermal conductivity 4th wall layer ground plate
!-- 17 - ground floor level height, 18 - z0 roughness ground floor level, 19 - z0h/z0g roughness heaat/humidity,
!-- 20 - wall fraction ground floor level, 21 - 1st wall layer thickness ground floor level, 
!-- 22 - 2nd wall layer thickness ground floor level, 23 - 3rd wall layer thickness ground floor level,
!-- 24 - 4th wall layer thickness ground floor level, 25 - heat capacity 1st/2nd wall layer ground floor level,
!-- 26 - heat capacity 3rd wall layer ground floor level, 27 - heat capacity 4th wall layer ground floor level,
!-- 28 - thermal conductivity 1st/2nd wall layer ground floor level, 
!-- 29 - thermal conductivity 3rd wall layer ground floor level, 30 - thermal conductivity 4th wall layer ground floor level
!-- 31 - wall emissivity ground floor level, 32 - wall albedo ground floor level, 33 - window fraction ground floor level,
!-- 34 - 1st window layer thickness ground floor level, 35 - 2nd window layer thickness ground floor level,
!-- 36 - 3rd window layer thickness ground floor level, 37 - 4th window layer thickness ground floor level,
!-- 38 - heat capacity 1st/2nd window layer ground floor level, 39 - heat capacity 3rd window layer ground floor level,
!-- 40 - heat capacity 4th window layer ground floor level, 
!-- 41 - thermal conductivity 1st/2nd window layer ground floor level, 
!-- 42 - thermal conductivity 3rd window layer ground floor level, 
!-- 43 - thermal conductivity 4th window layer ground floor level, 44 - window emissivity ground floor level,
!-- 45 - window transmissivity ground floor level, 46 - window albedo ground floor level,
!-- 47 - green fraction ground floor level, 48 - LAI on wall ground floor level, 49 - green emissivity ground floor level,
!-- 50 - green albedo ground floor level, 51 - z0 roughness above ground floor level, 
!-- 52 - z0h/z0g roughness heat/humidity above ground floor level, 53 - wall fraction above ground floor level
!-- 54 - 1st wall layer thickness above ground floor level, 55 - 2nd wall layer thickness above ground floor level
!-- 56 - 3rd wall layer thickness above ground floor level, 57 - 4th wall layer thickness above ground floor level 
!-- 58 - heat capacity 1st/2nd wall layer above ground floor level,
!-- 59 - heat capacity 3rd wall layer above ground floor level, 
!-- 60 - heat capacity 4th wall layer above ground floor level,
!-- 61 - thermal conductivity 1st/2nd wall layer above ground floor level,
!-- 62 - thermal conductivity 3rd wall layer above ground floor level,
!-- 63 - thermal conductivity 4th wall layer above ground floor level,
!-- 64 - wall emissivity above ground floor level, 65 - wall albedo above ground floor level,
!-- 66 - window fraction above ground floor level, 67 - 1st window layer thickness above ground floor level,
!-- 68 - 2nd thickness window layer above ground floor level, 69 - 3rd window layer thickness above ground floor level,
!-- 70 - 4th window layer thickness above ground floor level, 
!-- 71 - heat capacity 1st/2nd window layer above ground floor level,
!-- 72 - heat capacity 3rd window layer above ground floor level,
!-- 73 - heat capacity 4th window layer above ground floor level, 
!-- 74 - conductivity 1st/2nd window layer above ground floor level,
!-- 75 - thermal conductivity 3rd window layer above ground floor level, 
!-- 76 - thermal conductivity 4th window layer above ground floor level, 77 - window emissivity above ground floor level,
!-- 78 - window transmissivity above ground floor level, 79 - window albedo above ground floor level,
!-- 80 - green fraction above ground floor level, 81 - LAI on wall above ground floor level,
!-- 82 - green emissivity above ground floor level, 83 - green albedo above ground floor level,
!-- 84 - wall fraction roof, 85 - 1st wall layer thickness roof, 86 - 2nd wall layer thickness roof,
!-- 87 - 3rd wall layer thickness roof, 88 - 4th wall layer thickness roof,
!-- 89 - heat capacity 1st/2nd wall layer roof, 90 - heat capacity 3rd wall layer roof,
!-- 91 - heat capacity 4th wall layer roof, 92 - thermal conductivity 1st/2nd wall layer roof,
!-- 93 - thermal conductivity 3rd wall layer roof, 94 - thermal conductivity 4th wall layer roof,
!-- 95 - wall emissivity roof, 96 - wall albedo roof, 97 - window fraction roof,
!-- 98 - window 1st layer thickness roof, 99 - window 2nd layer thickness roof, 100 - window 3rd layer thickness roof,
!-- 101 - window 4th layer thickness, 102 - heat capacity 1st/2nd window layer roof, 
!-- 103 - heat capacity 3rd window layer roof, 104 - heat capacity 4th window layer roof,
!-- 105 - thermal conductivity 1st/2nd window layer roof, 106 - thermal conductivity 3rd window layer roof,
!-- 107 - thermal conductivity 4th window layer roof, 108 - window emissivity roof, 109 - window transmissivity roof,
!-- 110 - window albedo roof, 111 - green fraction roof ground floor level, 
!-- 112 - green fraction roof above ground floor level, 113 - LAI roof, 114 - green emissivity roof,
!-- 115 - green albedo roof, 116 - green type roof,
!-- Parameter for indoor model
!-- 117 - indoor target summer temperature, 118 - indoor target winter temperature, 
!-- 119 - shading factor, 120 - g-value windows, 121 - u-value windows, 122 - basical airflow without occupancy of the room,
!-- 123 - additional airflow depend of occupancy of the room, 124 - heat recovery efficiency,
!-- 125 - dynamic parameter specific effective surface, 126 - dynamic parameter innner heatstorage,
!-- 127 - ratio internal surface/floor area, 128 - maximal heating capacity, 129 - maximal cooling capacity,
!-- 130 - additional internal heat gains dependent on occupancy of the room, 
!-- 131 - basic internal heat gains without occupancy of the room, 132 - storey height, 133 - ceiling construction height


    REAL(wp), DIMENSION(0:133,1:7), PARAMETER :: building_pars = RESHAPE( (/   &
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.005_wp, 0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp,             & !parameter 6-11
        1400000.0_wp, 1300000.0_wp, 0.35_wp, 0.8_wp, 2.1_wp, 4.0_wp,           & !parameter 12-17
        0.01_wp, 0.001_wp, 0.75_wp,                                            & !parameter 18-20
        0.005_wp, 0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp,                     & !parameter 21-25
        1400000.0_wp, 1300000.0_wp, 0.35_wp,                                   & !parameter 26-28                     
        0.8_wp, 2.1_wp, 0.93_wp,                                               & !parameter 29-31       
        27.0_wp, 0.25_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,              & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.57_wp, 0.57_wp, 0.57_wp, 0.91_wp,                                    & !parameter 41-44
        0.75_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                             & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 0.7_wp, 0.005_wp,                        & !parameter 50-54
        0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp,                               & !parameter 55-58
        1400000.0_wp, 1300000.0_wp, 0.35_wp, 0.8_wp,                           & !parameter 59-62
        2.1_wp, 0.93_wp, 27.0_wp, 0.3_wp,                                      & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.57_wp, 0.57_wp, 0.57_wp, 0.91_wp, 0.75_wp,                           & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.005_wp, 0.01_wp, 0.31_wp, 0.63_wp, 2200000.0_wp, 1400000.0_wp,       & !parameter 85-90
        1300000.0_wp, 0.35_wp, 0.8_wp, 2.1_wp, 0.93_wp, 27.0_wp, 0.0_wp,       & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.57_wp, 0.57_wp, 0.57_wp,                 & !parameter 103-107
        0.91_wp, 0.75_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                     & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 0.76_wp, 5.0_wp,                         & !parameter 117-121
        0.1_wp, 0.5_wp, 0.0_wp, 3.5_wp, 370000.0_wp, 4.5_wp,                   & !parameter 122-127
        100000.0_wp, 0.0_wp, 3.0_wp, 10.0_wp, 3.0_wp, 0.2_wp,                  & !parameter 128-133- end of type 1
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.005_wp, 0.01_wp, 0.31_wp, 0.42_wp, 2000000.0_wp,             & !parameter 6-11
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.38_wp, 0.04_wp, 4.0_wp,           & !parameter 12-17
        0.01_wp, 0.001_wp, 0.78_wp,                                            & !parameter 18-20
        0.005_wp, 0.01_wp, 0.31_wp, 0.43_wp, 2000000.0_wp,                     & !parameter 21-25
        103000.0_wp, 900000.0_wp, 0.35_wp,                                     & !parameter 26-28                     
        0.38_wp, 0.04_wp, 0.92_wp,                                             & !parameter 29-31       
        27.0_wp, 0.22_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,              & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.11_wp, 0.11_wp, 0.11_wp, 0.11_wp,                                    & !parameter 41-44
        0.7_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                              & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 0.73_wp, 0.005_wp,                       & !parameter 50-54
        0.01_wp, 0.31_wp, 0.43_wp, 2000000.0_wp,                               & !parameter 55-58
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.38_wp,                            & !parameter 59-62
        0.04_wp, 0.92_wp, 27.0_wp, 0.27_wp,                                    & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.11_wp, 0.11_wp, 0.11_wp, 0.87_wp, 0.7_wp,                            & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.005_wp, 0.01_wp, 0.5_wp, 0.79_wp, 2000000.0_wp, 103000.0_wp,         & !parameter 85-90
        900000.0_wp, 0.35_wp, 0.38_wp, 0.04_wp, 0.93_wp, 27.0_wp, 0.0_wp,      & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.11_wp, 0.11_wp, 0.11_wp,                 & !parameter 103-107
        0.87_wp, 0.7_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                      & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 0.6_wp, 3.0_wp,                          & !parameter 117-121
        0.1_wp, 0.5_wp, 0.0_wp, 2.5_wp, 165000.0_wp, 4.5_wp,                   & !parameter 122-127
        100000.0_wp, 0.0_wp, 4.0_wp, 8.0_wp, 3.0_wp, 0.2_wp,                   & !parameter 128-133- end of type 2
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.005_wp, 0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp,              & !parameter 6-11
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.14_wp, 0.035_wp, 4.0_wp,          & !parameter 12-17
        0.01_wp, 0.001_wp, 0.75_wp,                                            & !parameter 18-20
        0.005_wp, 0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp,                      & !parameter 21-25
        103000.0_wp, 900000.0_wp, 0.35_wp,                                     & !parameter 26-28                     
        0.14_wp, 0.035_wp, 0.92_wp,                                            & !parameter 29-31       
        27.0_wp, 0.25_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,              & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.037_wp, 0.037_wp, 0.037_wp, 0.8_wp,                                     & !parameter 41-44
        0.6_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                              & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 0.7_wp, 0.005_wp,                        & !parameter 50-54
        0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp,                                & !parameter 55-58
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.14_wp,                            & !parameter 59-62
        0.035_wp, 0.92_wp, 27.0_wp, 0.3_wp,                                    & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.037_wp, 0.037_wp, 0.037_wp, 0.8_wp, 0.6_wp,                             & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.005_wp, 0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp, 103000.0_wp,         & !parameter 85-90
        900000.0_wp, 0.35_wp, 0.14_wp, 0.035_wp, 0.93_wp, 27.0_wp, 0.0_wp,     & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.037_wp, 0.037_wp, 0.037_wp,                 & !parameter 103-107
        0.8_wp, 0.6_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                       & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 0.5_wp, 0.6_wp,                          & !parameter 117-121
        0.1_wp, 0.5_wp, 0.8_wp, 2.5_wp, 80000.0_wp, 4.5_wp,                    & !parameter 122-127
        100000.0_wp, 0.0_wp, 3.0_wp, 8.0_wp, 3.0_wp, 0.2_wp,                   & !parameter 128-133- end of type 3
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.005_wp, 0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp,             & !parameter 6-11
        1400000.0_wp, 1300000.0_wp, 0.35_wp, 0.8_wp, 2.1_wp, 4.0_wp,           & !parameter 12-17
        0.01_wp, 0.001_wp, 0.55_wp,                                            & !parameter 18-20
        0.005_wp, 0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp,                     & !parameter 21-25
        1400000.0_wp, 1300000.0_wp, 0.35_wp,                                   & !parameter 26-28                     
        0.8_wp, 2.1_wp, 0.93_wp,                                               & !parameter 29-31       
        27.0_wp, 0.45_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,              & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.57_wp, 0.57_wp, 0.57_wp, 0.91_wp,                                    & !parameter 41-44
        0.75_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                             & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 0.5_wp, 0.005_wp,                        & !parameter 50-54
        0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp,                               & !parameter 55-58
        1400000.0_wp, 1300000.0_wp, 0.35_wp, 0.8_wp,                           & !parameter 59-62
        2.1_wp, 0.93_wp, 27.0_wp, 0.5_wp,                                      & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.57_wp, 0.57_wp, 0.57_wp, 0.91_wp, 0.75_wp,                           & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.005_wp, 0.01_wp, 0.39_wp, 0.63_wp, 2200000.0_wp, 1400000.0_wp,       & !parameter 85-90
        1300000.0_wp, 0.35_wp, 0.8_wp, 2.1_wp, 0.93_wp, 27.0_wp, 0.0_wp,       & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.57_wp, 0.57_wp, 0.57_wp,                 & !parameter 103-107
        0.91_wp, 0.75_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                     & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 0.76_wp, 5.0_wp,                         & !parameter 117-121
        0.1_wp, 1.5_wp, 0.0_wp, 3.5_wp, 370000.0_wp, 4.5_wp,                   & !parameter 122-127
        100000.0_wp, 0.0_wp, 3.0_wp, 10.0_wp, 3.0_wp, 0.2_wp,                  & !parameter 128-133- end of type 4
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.005_wp, 0.01_wp, 0.31_wp, 0.43_wp, 2000000.0_wp,             & !parameter 6-11
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.38_wp, 0.04_wp, 4.0_wp,           & !parameter 12-17
        0.01_wp, 0.001_wp, 0.55_wp,                                            & !parameter 18-20
        0.005_wp, 0.01_wp, 0.31_wp, 0.43_wp, 2000000.0_wp,                     & !parameter 21-25
        103000.0_wp, 900000.0_wp, 0.35_wp,                                     & !parameter 26-28                     
        0.38_wp, 0.04_wp, 0.92_wp,                                             & !parameter 29-31       
        27.0_wp, 0.45_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,              & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.11_wp, 0.11_wp, 0.11_wp, 0.87_wp,                                    & !parameter 41-44
        0.7_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                              & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 0.5_wp, 0.005_wp,                        & !parameter 50-54
        0.01_wp, 0.31_wp, 0.43_wp, 2000000.0_wp,                               & !parameter 55-58
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.38_wp,                            & !parameter 59-62
        0.04_wp, 0.92_wp, 27.0_wp, 0.5_wp,                                     & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.11_wp, 0.11_wp, 0.11_wp, 0.87_wp, 0.7_wp,                            & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.005_wp, 0.01_wp, 0.31_wp, 0.43_wp, 2000000.0_wp, 103000.0_wp,        & !parameter 85-90
        900000.0_wp, 0.35_wp, 0.38_wp, 0.04_wp, 0.91_wp, 27.0_wp, 0.0_wp,      & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.11_wp, 0.11_wp, 0.11_wp,                 & !parameter 103-107
        0.87_wp, 0.7_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                      & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 0.6_wp, 3.0_wp,                          & !parameter 117-121
        0.1_wp, 1.5_wp, 0.65_wp, 2.5_wp, 165000.0_wp, 4.5_wp,                  & !parameter 122-127
        100000.0_wp, 0.0_wp, 7.0_wp, 20.0_wp, 3.0_wp, 0.2_wp,                  & !parameter 128-133- end of type 5
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.005_wp, 0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp,              & !parameter 6-11
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.14_wp, 0.035_wp, 4.0_wp,          & !parameter 12-17
        0.01_wp, 0.001_wp, 0.475_wp,                                           & !parameter 18-20
        0.005_wp, 0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp,                      & !parameter 21-25
        103000.0_wp, 900000.0_wp, 0.35_wp,                                     & !parameter 26-28                     
        0.14_wp, 0.035_wp, 0.92_wp,                                            & !parameter 29-31       
        27.0_wp, 0.525_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,             & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.037_wp, 0.037_wp, 0.037_wp, 0.8_wp,                                     & !parameter 41-44
        0.6_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                              & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 0.425_wp, 0.005_wp,                      & !parameter 50-54
        0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp,                                & !parameter 55-58
        103000.0_wp, 900000.0_wp, 0.35_wp, 0.14_wp,                            & !parameter 59-62
        0.035_wp, 0.92_wp, 27.0_wp, 0.575_wp,                                  & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.037_wp, 0.037_wp, 0.037_wp, 0.8_wp, 0.6_wp,                             & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.005_wp, 0.01_wp, 0.41_wp, 0.7_wp, 2000000.0_wp, 103000.0_wp,         & !parameter 85-90
        900000.0_wp, 0.35_wp, 0.14_wp, 0.035_wp, 0.91_wp, 27.0_wp, 0.0_wp,     & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.037_wp, 0.037_wp, 0.037_wp,                 & !parameter 103-107
        0.8_wp, 0.6_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                       & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 0.5_wp, 0.6_wp,                          & !parameter 117-121
        0.1_wp, 1.5_wp, 0.9_wp, 2.5_wp, 80000.0_wp, 4.5_wp,                    & !parameter 122-127
        100000.0_wp, 0.0_wp, 5.0_wp, 15.0_wp, 3.0_wp, 0.2_wp,                  & !parameter 128-133- end of type 6    
        10.0_wp, 10.0_wp, 20000.0_wp, 23.0_wp, 23.0_wp, 10.0_wp,               & !parameter 0-5
        1.0_wp, 0.29_wp, 0.295_wp, 0.695_wp, 0.985_wp, 1950400.0_wp,           & !parameter 6-11
        1848000.0_wp, 1848000.0_wp, 0.7_wp, 1.0_wp, 1.0_wp, 4.0_wp,            & !parameter 12-17
        0.01_wp, 0.001_wp, 1.0_wp,                                             & !parameter 18-20
        0.29_wp, 0.295_wp, 0.695_wp, 0.985_wp, 1950400.0_wp,                   & !parameter 21-25
        1848000.0_wp, 1848000.0_wp, 0.7_wp,                                    & !parameter 26-28                     
        1.0_wp, 1.0_wp, 0.9_wp,                                                & !parameter 29-31       
        27.0_wp, 0.0_wp, 0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,               & !parameter 32-37
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 38-40
        0.57_wp, 0.57_wp, 0.57_wp, 0.8_wp,                                     & !parameter 41-44
        0.6_wp, 27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp,                              & !parameter 45-49
        5.0_wp, 0.001_wp, 0.0001_wp, 1.0_wp, 0.29_wp,                         & !parameter 50-54
        0.295_wp, 0.695_wp, 0.985_wp, 1950400.0_wp,                            & !parameter 55-58
        1848000.0_wp, 1848000.0_wp, 0.7_wp, 1.0_wp,                            & !parameter 59-62
        1.0_wp, 0.9_wp, 27.0_wp, 0.0_wp,                                       & !parameter 63-66
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp,                                & !parameter 67-70
        1736000.0_wp, 1736000.0_wp, 1736000.0_wp,                              & !parameter 71-73
        0.57_wp, 0.57_wp, 0.57_wp, 0.8_wp, 0.6_wp,                             & !parameter 74-78
        27.0_wp, 0.0_wp, 1.5_wp, 0.86_wp, 5.0_wp, 1.0_wp,                      & !parameter 79-84
        0.29_wp, 0.295_wp, 0.695_wp, 0.985_wp, 1950400.0_wp, 1848000.0_wp,     & !parameter 85-90
        1848000.0_wp, 0.7_wp, 1.0_wp, 1.0_wp, 0.9_wp, 27.0_wp, 0.0_wp,         & !parameter 91-97
        0.003_wp, 0.006_wp, 0.012_wp, 0.018_wp, 1736000.0_wp,                  & !parameter 98-102 
        1736000.0_wp, 1736000.0_wp, 0.57_wp, 0.57_wp, 0.57_wp,                 & !parameter 103-107
        0.8_wp, 0.6_wp, 27.0_wp, 0.0_wp, 0.0_wp, 1.5_wp,                       & !parameter 108-113
        0.86_wp, 5.0_wp, 0.0_wp,                                              & !parameter 114-116
        299.15_wp, 293.15_wp, 0.8_wp, 100.0_wp, 100.0_wp,                      & !parameter 117-121
        20.0_wp, 20.0_wp, 0.0_wp, 1.0_wp, 1.0_wp, 4.5_wp,                      & !parameter 122-127
        100000.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp, 0.2_wp                    & !parameter 128-133- end of type 7 (bridge)
                                                                       /),     &
                                                               (/134, 7/) )

!
!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls. 
    TYPE t_surf_vertical
       REAL(wp), DIMENSION(:), ALLOCATABLE         :: t
    END TYPE t_surf_vertical
!
!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls. 
    TYPE t_wall_vertical
       REAL(wp), DIMENSION(:,:), ALLOCATABLE       :: t
    END TYPE t_wall_vertical

    TYPE surf_type_usm
       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_usm_1d !< 1D prognostic variable
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_usm_2d !< 2D prognostic variable
    END TYPE surf_type_usm
    
    TYPE(surf_type_usm), POINTER  ::  m_liq_usm_h,        & !< liquid water reservoir (m), horizontal surface elements 
                                      m_liq_usm_h_p         !< progn. liquid water reservoir (m), horizontal surface elements 

    TYPE(surf_type_usm), TARGET   ::  m_liq_usm_h_1,      & !<
                                      m_liq_usm_h_2         !<

    TYPE(surf_type_usm), DIMENSION(:), POINTER  ::    &
                                      m_liq_usm_v,        & !< liquid water reservoir (m), vertical surface elements 
                                      m_liq_usm_v_p         !< progn. liquid water reservoir (m), vertical surface elements 

    TYPE(surf_type_usm), DIMENSION(0:3), TARGET   ::  &
                                      m_liq_usm_v_1,      & !<
                                      m_liq_usm_v_2         !<

    TYPE(surf_type_usm), TARGET ::  tm_liq_usm_h_m      !< liquid water reservoir tendency (m), horizontal surface elements 
    TYPE(surf_type_usm), DIMENSION(0:3), TARGET ::  tm_liq_usm_v_m      !< liquid water reservoir tendency (m), vertical surface elements 


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- anthropogenic heat sources
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE        ::  aheat             !< daily average of anthropogenic heat (W/m2)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheatprof         !< diurnal profiles of anthropogenic heat for particular layers
    INTEGER(iwp)                                   ::  naheatlayers = 1  !< number of layers of anthropogenic heat

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- wall surface model
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- wall surface model constants
    INTEGER(iwp), PARAMETER                        :: nzb_wall = 0       !< inner side of the wall model (to be switched)
    INTEGER(iwp), PARAMETER                        :: nzt_wall = 3       !< outer side of the wall model (to be switched)
    INTEGER(iwp), PARAMETER                        :: nzw = 4            !< number of wall layers (fixed for now)

    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
                                                                         !< normalized soil, wall and roof layer depths (m/m)
!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default = (/0.33_wp, 0.66_wp, 1.0_wp /)
    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.33_wp, 0.66_wp, 1.0_wp /)
!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
                                                                         !< normalized window layer depths (m/m)
!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
                                                                         !< normalized green layer depths (m/m)
    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /)


    REAL(wp)                                       :: wall_inner_temperature = 295.0_wp    !< temperature of the inner wall surface (~22 degrees C) (K)
    REAL(wp)                                       :: roof_inner_temperature = 295.0_wp    !< temperature of the inner roof surface (~22 degrees C) (K)
    REAL(wp)                                       :: soil_inner_temperature = 288.0_wp    !< temperature of the deep soil (~15 degrees C) (K)
    REAL(wp)                                       :: window_inner_temperature = 295.0_wp  !< temperature of the inner window surface (~22 degrees C) (K)

    REAL(wp)                                       ::   m_total = 0.0_wp !< weighted total water content of the soil (m3/m3)
    INTEGER(iwp)                                   ::   soil_type
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- surface and material model variables for walls, ground, roofs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn                !< normalized wall layer depths (m)
    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn_window         !< normalized window layer depths (m)
    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn_green          !< normalized green layer depths (m)

    REAL(wp), DIMENSION(:), POINTER                :: t_surf_wall_h
    REAL(wp), DIMENSION(:), POINTER                :: t_surf_wall_h_p 
    REAL(wp), DIMENSION(:), POINTER                :: t_surf_window_h
    REAL(wp), DIMENSION(:), POINTER                :: t_surf_window_h_p 
    REAL(wp), DIMENSION(:), POINTER                :: t_surf_green_h
    REAL(wp), DIMENSION(:), POINTER                :: t_surf_green_h_p 

    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_wall_h_1
    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_wall_h_2
    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_1
    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_2
    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_1
    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_2

    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_wall_v
    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_wall_v_p
    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_window_v
    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_window_v_p
    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_green_v
    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_green_v_p

    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_wall_v_1
    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_wall_v_2
    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_window_v_1
    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_window_v_2
    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_green_v_1
    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_green_v_2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- Energy balance variables
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-- parameters of the land, roof and wall surfaces

    REAL(wp), DIMENSION(:,:), POINTER                :: t_wall_h, t_wall_h_p
    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_1, t_wall_h_2
    REAL(wp), DIMENSION(:,:), POINTER                :: t_window_h, t_window_h_p
    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_1, t_window_h_2
    REAL(wp), DIMENSION(:,:), POINTER                :: t_green_h, t_green_h_p
    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_1, t_green_h_2
    REAL(wp), DIMENSION(:,:), POINTER                :: swc_h, rootfr_h, wilt_h, fc_h, swc_sat_h, swc_h_p, swc_res_h
    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: swc_h_1, rootfr_h_1, &
                                                        wilt_h_1, fc_h_1, swc_sat_h_1, swc_h_2, swc_res_h_1
    

    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_wall_v, t_wall_v_p
    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_1, t_wall_v_2
    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_window_v, t_window_v_p
    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_1, t_window_v_2
    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_green_v, t_green_v_p
    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_1, t_green_v_2
    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: swc_v, swc_v_p
    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: swc_v_1, swc_v_2

!-- Surface and material parameters classes (surface_type)
!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
    INTEGER(iwp)                                   :: n_surface_types      !< number of the wall type categories
    INTEGER(iwp), PARAMETER                        :: n_surface_params = 9 !< number of parameters for each type of the wall
    INTEGER(iwp), PARAMETER                        :: ialbedo  = 1         !< albedo of the surface
    INTEGER(iwp), PARAMETER                        :: iemiss   = 2         !< emissivity of the surface
    INTEGER(iwp), PARAMETER                        :: ilambdas = 3         !< heat conductivity lambda S between surface and material ( W m-2 K-1 )
    INTEGER(iwp), PARAMETER                        :: irough   = 4         !< roughness length z0 for movements
    INTEGER(iwp), PARAMETER                        :: iroughh  = 5         !< roughness length z0h for scalars (heat, humidity,...)
    INTEGER(iwp), PARAMETER                        :: icsurf   = 6         !< Surface skin layer heat capacity (J m-2 K-1 )
    INTEGER(iwp), PARAMETER                        :: ithick   = 7         !< thickness of the surface (wall, roof, land)  ( m )
    INTEGER(iwp), PARAMETER                        :: irhoC    = 8         !< volumetric heat capacity rho*C of the material ( J m-3 K-1 )
    INTEGER(iwp), PARAMETER                        :: ilambdah = 9         !< thermal conductivity lambda H of the wall (W m-1 K-1 )
    CHARACTER(12), DIMENSION(:), ALLOCATABLE       :: surface_type_names   !< names of wall types (used only for reports)
    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        :: surface_type_codes   !< codes of wall types
    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: surface_params       !< parameters of wall types

    
!-- interfaces of subroutines accessed from outside of this module
    INTERFACE usm_boundary_condition
       MODULE PROCEDURE usm_boundary_condition
    END INTERFACE usm_boundary_condition

    INTERFACE usm_check_data_output
       MODULE PROCEDURE usm_check_data_output
    END INTERFACE usm_check_data_output
    
    INTERFACE usm_check_parameters
       MODULE PROCEDURE usm_check_parameters
    END INTERFACE usm_check_parameters
    
    INTERFACE usm_data_output_3d
       MODULE PROCEDURE usm_data_output_3d
    END INTERFACE usm_data_output_3d
    
    INTERFACE usm_define_netcdf_grid
       MODULE PROCEDURE usm_define_netcdf_grid
    END INTERFACE usm_define_netcdf_grid

    INTERFACE usm_init_urban_surface
       MODULE PROCEDURE usm_init_urban_surface
    END INTERFACE usm_init_urban_surface

    INTERFACE usm_material_heat_model
       MODULE PROCEDURE usm_material_heat_model
    END INTERFACE usm_material_heat_model
    
    INTERFACE usm_green_heat_model
       MODULE PROCEDURE usm_green_heat_model
    END INTERFACE usm_green_heat_model
    
    INTERFACE usm_parin
       MODULE PROCEDURE usm_parin
    END INTERFACE usm_parin

    INTERFACE usm_rrd_local 
       MODULE PROCEDURE usm_rrd_local
    END INTERFACE usm_rrd_local

    INTERFACE usm_surface_energy_balance
       MODULE PROCEDURE usm_surface_energy_balance
    END INTERFACE usm_surface_energy_balance
    
    INTERFACE usm_swap_timelevel
       MODULE PROCEDURE usm_swap_timelevel
    END INTERFACE usm_swap_timelevel
        
    INTERFACE usm_wrd_local
       MODULE PROCEDURE usm_wrd_local
    END INTERFACE usm_wrd_local

    INTERFACE usm_allocate_surface
       MODULE PROCEDURE usm_allocate_surface
    END INTERFACE usm_allocate_surface

    INTERFACE usm_average_3d_data
       MODULE PROCEDURE usm_average_3d_data
    END INTERFACE usm_average_3d_data

    
    SAVE

    PRIVATE 
    
!-- Public functions
    PUBLIC usm_boundary_condition, usm_check_parameters, usm_init_urban_surface,&
           usm_rrd_local,                                                      & 
           usm_surface_energy_balance, usm_material_heat_model,                &
           usm_swap_timelevel, usm_check_data_output, usm_average_3d_data,     &
           usm_data_output_3d, usm_define_netcdf_grid, usm_parin,              &
           usm_wrd_local, usm_allocate_surface

!-- Public parameters, constants and initial values
    PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod,           &
           usm_green_heat_model, building_pars,  &
           nzt_wall, t_wall_h, t_wall_v,         &
           t_window_h, t_window_v, building_type



 CONTAINS

!------------------------------------------------------------------------------!
! Description:
! ------------
!> This subroutine creates the necessary indices of the urban surfaces
!> and plant canopy and it allocates the needed arrays for USM
!------------------------------------------------------------------------------!
    SUBROUTINE usm_allocate_surface
    
        IMPLICIT NONE
       
        INTEGER(iwp) ::  l

!
!--     Allocate radiation arrays which are part of the new data type. 
!--     For horizontal surfaces.
        ALLOCATE( surf_usm_h%surfhf(1:surf_usm_h%ns)    )
        ALLOCATE( surf_usm_h%rad_net_l(1:surf_usm_h%ns) )
!
!--     For vertical surfaces
        DO  l = 0, 3
           ALLOCATE( surf_usm_v(l)%surfhf(1:surf_usm_v(l)%ns)    )
           ALLOCATE( surf_usm_v(l)%rad_net_l(1:surf_usm_v(l)%ns) )
        ENDDO

!--     Wall surface model
!--     allocate arrays for wall surface model and define pointers
       
!--     allocate array of wall types and wall parameters
        ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns)      )
        ALLOCATE ( surf_usm_h%building_type(1:surf_usm_h%ns)      )
        ALLOCATE ( surf_usm_h%building_type_name(1:surf_usm_h%ns) )
        surf_usm_h%building_type      = 0
        surf_usm_h%building_type_name = 'none'
        DO  l = 0, 3
           ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%building_type(1:surf_usm_v(l)%ns)      )
           ALLOCATE ( surf_usm_v(l)%building_type_name(1:surf_usm_v(l)%ns) )
           surf_usm_v(l)%building_type      = 0
           surf_usm_v(l)%building_type_name = 'none'
        ENDDO
!
!--     Allocate albedo_type and albedo. Each surface element
!--     has 3 values, 0: wall fraction, 1: green fraction, 2: window fraction.
        ALLOCATE( surf_usm_h%albedo_type(0:2,1:surf_usm_h%ns) )
        ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)      )
        surf_usm_h%albedo_type = albedo_type
        DO  l = 0, 3
           ALLOCATE( surf_usm_v(l)%albedo_type(0:2,1:surf_usm_v(l)%ns) )
           ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns)      )
           surf_usm_v(l)%albedo_type = albedo_type
        ENDDO       


!
!--     Allocate indoor target temperature for summer and winter
        ALLOCATE( surf_usm_h%target_temp_summer(1:surf_usm_h%ns) )
        ALLOCATE( surf_usm_h%target_temp_winter(1:surf_usm_h%ns) )
        DO  l = 0, 3
           ALLOCATE( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) )
           ALLOCATE( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) )
        ENDDO   
!
!--     Allocate flag indicating ground floor level surface elements
        ALLOCATE ( surf_usm_h%ground_level(1:surf_usm_h%ns) ) 
        DO  l = 0, 3
           ALLOCATE( surf_usm_v(l)%ground_level(1:surf_usm_v(l)%ns) )
        ENDDO   
!
!--      Allocate arrays for relative surface fraction. 
!--      0 - wall fraction, 1 - green fraction, 2 - window fraction
         ALLOCATE( surf_usm_h%frac(0:2,1:surf_usm_h%ns) )
         surf_usm_h%frac = 0.0_wp
         DO  l = 0, 3
            ALLOCATE( surf_usm_v(l)%frac(0:2,1:surf_usm_v(l)%ns) )
            surf_usm_v(l)%frac = 0.0_wp
         ENDDO
       
!--     wall and roof surface parameters. First for horizontal surfaces
        ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns)        )
        ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns)        )
        ALLOCATE ( surf_usm_h%lambda_surf_window(1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%lambda_surf_green(1:surf_usm_h%ns)  )
        ALLOCATE ( surf_usm_h%c_surface(1:surf_usm_h%ns)          )
        ALLOCATE ( surf_usm_h%c_surface_window(1:surf_usm_h%ns)   )
        ALLOCATE ( surf_usm_h%c_surface_green(1:surf_usm_h%ns)    )
        ALLOCATE ( surf_usm_h%transmissivity(1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%lai(1:surf_usm_h%ns)                )
        ALLOCATE ( surf_usm_h%emissivity(0:2,1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%r_a(1:surf_usm_h%ns)                )
        ALLOCATE ( surf_usm_h%r_a_green(1:surf_usm_h%ns)          )
        ALLOCATE ( surf_usm_h%r_a_window(1:surf_usm_h%ns)         )
        ALLOCATE ( surf_usm_h%green_type_roof(1:surf_usm_h%ns)    )
        ALLOCATE ( surf_usm_h%r_s(1:surf_usm_h%ns)                )
        
!
!--     For vertical surfaces.
        DO  l = 0, 3
           ALLOCATE ( surf_usm_v(l)%lambda_surf(1:surf_usm_v(l)%ns)     )
           ALLOCATE ( surf_usm_v(l)%c_surface(1:surf_usm_v(l)%ns)       )
           ALLOCATE ( surf_usm_v(l)%lambda_surf_window(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%c_surface_window(1:surf_usm_v(l)%ns)   )
           ALLOCATE ( surf_usm_v(l)%lambda_surf_green(1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%c_surface_green(1:surf_usm_v(l)%ns)    )
           ALLOCATE ( surf_usm_v(l)%transmissivity(1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%lai(1:surf_usm_v(l)%ns)             )
           ALLOCATE ( surf_usm_v(l)%emissivity(0:2,1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%r_a(1:surf_usm_v(l)%ns)             )
           ALLOCATE ( surf_usm_v(l)%r_a_green(1:surf_usm_v(l)%ns)       )
           ALLOCATE ( surf_usm_v(l)%r_a_window(1:surf_usm_v(l)%ns)      )           
           ALLOCATE ( surf_usm_v(l)%r_s(1:surf_usm_v(l)%ns)                )
        ENDDO

!       
!--     allocate wall and roof material parameters. First for horizontal surfaces
        ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns)               )
        ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns)                  )
        ALLOCATE ( surf_usm_h%thickness_green(1:surf_usm_h%ns)                   )
        ALLOCATE ( surf_usm_h%lambda_h(nzb_wall:nzt_wall,1:surf_usm_h%ns)   )
        ALLOCATE ( surf_usm_h%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_h%ns)    )
        ALLOCATE ( surf_usm_h%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
        ALLOCATE ( surf_usm_h%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )

        ALLOCATE ( surf_usm_h%rho_c_total_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%n_vg_green(1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%alpha_vg_green(1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%l_vg_green(1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%gamma_w_green_sat(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%lambda_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%gamma_w_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%tswc_h_m(nzb_wall:nzt_wall,1:surf_usm_h%ns)     )

!
!--     For vertical surfaces.
        DO  l = 0, 3
           ALLOCATE ( surf_usm_v(l)%thickness_wall(1:surf_usm_v(l)%ns)               )
           ALLOCATE ( surf_usm_v(l)%thickness_window(1:surf_usm_v(l)%ns)                  )
           ALLOCATE ( surf_usm_v(l)%thickness_green(1:surf_usm_v(l)%ns)                   )
           ALLOCATE ( surf_usm_v(l)%lambda_h(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)   )
           ALLOCATE ( surf_usm_v(l)%rho_c_wall(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%lambda_h_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%rho_c_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)    )
           ALLOCATE ( surf_usm_v(l)%lambda_h_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%rho_c_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)     )
        ENDDO

!
!--     allocate green wall and roof vegetation and soil parameters. First horizontal surfaces
        ALLOCATE ( surf_usm_h%g_d(1:surf_usm_h%ns)              )
        ALLOCATE ( surf_usm_h%c_liq(1:surf_usm_h%ns)            )
        ALLOCATE ( surf_usm_h%qsws_liq(1:surf_usm_h%ns)         )
        ALLOCATE ( surf_usm_h%qsws_veg(1:surf_usm_h%ns)         )
        ALLOCATE ( surf_usm_h%r_canopy(1:surf_usm_h%ns)         )
        ALLOCATE ( surf_usm_h%r_canopy_min(1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%qsws_eb(1:surf_usm_h%ns)          )
        ALLOCATE ( surf_usm_h%pt_10cm(1:surf_usm_h%ns)          ) 
        ALLOCATE ( surf_usm_h%pt_2m(1:surf_usm_h%ns)            ) 

!
!--     For vertical surfaces.
        DO  l = 0, 3
          ALLOCATE ( surf_usm_v(l)%g_d(1:surf_usm_v(l)%ns)              )
          ALLOCATE ( surf_usm_v(l)%c_liq(1:surf_usm_v(l)%ns)            )
          ALLOCATE ( surf_usm_v(l)%qsws_liq(1:surf_usm_v(l)%ns)         )
          ALLOCATE ( surf_usm_v(l)%qsws_veg(1:surf_usm_v(l)%ns)         )
          ALLOCATE ( surf_usm_v(l)%qsws_eb(1:surf_usm_v(l)%ns)          )
          ALLOCATE ( surf_usm_v(l)%r_canopy(1:surf_usm_v(l)%ns)         )
          ALLOCATE ( surf_usm_v(l)%r_canopy_min(1:surf_usm_v(l)%ns)     )
          ALLOCATE ( surf_usm_v(l)%pt_10cm(1:surf_usm_v(l)%ns)     )
        ENDDO
        
!--     allocate wall and roof layers sizes. For horizontal surfaces.
        ALLOCATE ( zwn(nzb_wall:nzt_wall) )
        ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
        ALLOCATE ( zwn_window(nzb_wall:nzt_wall) )
        ALLOCATE ( surf_usm_h%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
        ALLOCATE ( zwn_green(nzb_wall:nzt_wall) )
        ALLOCATE ( surf_usm_h%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)      )
        ALLOCATE ( surf_usm_h%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)    )
        ALLOCATE ( surf_usm_h%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
        ALLOCATE ( surf_usm_h%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%zw(nzb_wall:nzt_wall,1:surf_usm_h%ns)            )
        ALLOCATE ( surf_usm_h%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)    )
        ALLOCATE ( surf_usm_h%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
        ALLOCATE ( surf_usm_h%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%zw_window(nzb_wall:nzt_wall,1:surf_usm_h%ns)       )
        ALLOCATE ( surf_usm_h%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)     )
        ALLOCATE ( surf_usm_h%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)   )
        ALLOCATE ( surf_usm_h%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_h%ns)  )
        ALLOCATE ( surf_usm_h%zw_green(nzb_wall:nzt_wall,1:surf_usm_h%ns)        )
!
!--     For vertical surfaces.
        DO  l = 0, 3
           ALLOCATE ( surf_usm_v(l)%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
           ALLOCATE ( surf_usm_v(l)%dz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
           ALLOCATE ( surf_usm_v(l)%dz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)      )
           ALLOCATE ( surf_usm_v(l)%ddz_wall(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)    )
           ALLOCATE ( surf_usm_v(l)%dz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%ddz_wall_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%zw(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)            )
           ALLOCATE ( surf_usm_v(l)%ddz_window(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)    )
           ALLOCATE ( surf_usm_v(l)%dz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%ddz_window_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)       )
           ALLOCATE ( surf_usm_v(l)%ddz_green(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)     )
           ALLOCATE ( surf_usm_v(l)%dz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)   )
           ALLOCATE ( surf_usm_v(l)%ddz_green_stag(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)  )
           ALLOCATE ( surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns)        )
        ENDDO

!--     allocate wall and roof temperature arrays, for horizontal walls
!
!--     Allocate if required. Note, in case of restarts, some of these arrays 
!--     might be already allocated.
        IF ( .NOT. ALLOCATED( t_surf_wall_h_1 ) )                                   &
           ALLOCATE ( t_surf_wall_h_1(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( t_surf_wall_h_2 ) )                                   &
           ALLOCATE ( t_surf_wall_h_2(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( t_wall_h_1 ) )                                   &           
           ALLOCATE ( t_wall_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( t_wall_h_2 ) )                                   &           
           ALLOCATE ( t_wall_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )          
        IF ( .NOT. ALLOCATED( t_surf_window_h_1 ) )                            &
           ALLOCATE ( t_surf_window_h_1(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( t_surf_window_h_2 ) )                            &
           ALLOCATE ( t_surf_window_h_2(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( t_window_h_1 ) )                                 &           
           ALLOCATE ( t_window_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( t_window_h_2 ) )                                 &           
           ALLOCATE ( t_window_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )          
        IF ( .NOT. ALLOCATED( t_surf_green_h_1 ) )                             &
           ALLOCATE ( t_surf_green_h_1(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( t_surf_green_h_2 ) )                             &
           ALLOCATE ( t_surf_green_h_2(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( t_green_h_1 ) )                                  &           
           ALLOCATE ( t_green_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( t_green_h_2 ) )                                  &           
           ALLOCATE ( t_green_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )          
        IF ( .NOT. ALLOCATED( swc_h_1 ) )                                  &           
           ALLOCATE ( swc_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( swc_sat_h_1 ) )                                  &           
           ALLOCATE ( swc_sat_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( swc_res_h_1 ) )                                  &           
           ALLOCATE ( swc_res_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( swc_h_2 ) )                                  &           
           ALLOCATE ( swc_h_2(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( rootfr_h_1 ) )                                  &           
           ALLOCATE ( rootfr_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( wilt_h_1 ) )                                  &           
           ALLOCATE ( wilt_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 
        IF ( .NOT. ALLOCATED( fc_h_1 ) )                                  &           
           ALLOCATE ( fc_h_1(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) 

        IF ( .NOT. ALLOCATED( m_liq_usm_h_1%var_usm_1d ) )                         &
           ALLOCATE ( m_liq_usm_h_1%var_usm_1d(1:surf_usm_h%ns) )
        IF ( .NOT. ALLOCATED( m_liq_usm_h_2%var_usm_1d ) )                         &
           ALLOCATE ( m_liq_usm_h_2%var_usm_1d(1:surf_usm_h%ns) )
           
!           
!--     initial assignment of the pointers
        t_wall_h    => t_wall_h_1;    t_wall_h_p    => t_wall_h_2
        t_window_h    => t_window_h_1;    t_window_h_p    => t_window_h_2
        t_green_h    => t_green_h_1;    t_green_h_p    => t_green_h_2
        t_surf_wall_h => t_surf_wall_h_1; t_surf_wall_h_p => t_surf_wall_h_2           
        t_surf_window_h => t_surf_window_h_1; t_surf_window_h_p => t_surf_window_h_2  
        t_surf_green_h => t_surf_green_h_1; t_surf_green_h_p => t_surf_green_h_2           
        m_liq_usm_h     => m_liq_usm_h_1;     m_liq_usm_h_p     => m_liq_usm_h_2
        swc_h       => swc_h_1; swc_h_p     => swc_h_2
        swc_sat_h    => swc_sat_h_1
        swc_res_h    => swc_res_h_1
        rootfr_h       => rootfr_h_1
        wilt_h       => wilt_h_1
        fc_h       => fc_h_1

!--     allocate wall and roof temperature arrays, for vertical walls if required
!
!--     Allocate if required. Note, in case of restarts, some of these arrays 
!--     might be already allocated.
        DO  l = 0, 3
           IF ( .NOT. ALLOCATED( t_surf_wall_v_1(l)%t ) )                           &
              ALLOCATE ( t_surf_wall_v_1(l)%t(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( t_surf_wall_v_2(l)%t ) )                           &
              ALLOCATE ( t_surf_wall_v_2(l)%t(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( t_wall_v_1(l)%t ) )                           &           
              ALLOCATE ( t_wall_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
           IF ( .NOT. ALLOCATED( t_wall_v_2(l)%t ) )                           &           
              ALLOCATE ( t_wall_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )  
           IF ( .NOT. ALLOCATED( t_surf_window_v_1(l)%t ) )                    &
              ALLOCATE ( t_surf_window_v_1(l)%t(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( t_surf_window_v_2(l)%t ) )                    &
              ALLOCATE ( t_surf_window_v_2(l)%t(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( t_window_v_1(l)%t ) )                         &           
              ALLOCATE ( t_window_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
           IF ( .NOT. ALLOCATED( t_window_v_2(l)%t ) )                         &           
              ALLOCATE ( t_window_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )  
           IF ( .NOT. ALLOCATED( t_surf_green_v_1(l)%t ) )                     &
              ALLOCATE ( t_surf_green_v_1(l)%t(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( t_surf_green_v_2(l)%t ) )                     &
              ALLOCATE ( t_surf_green_v_2(l)%t(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( t_green_v_1(l)%t ) )                          &           
              ALLOCATE ( t_green_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
           IF ( .NOT. ALLOCATED( t_green_v_2(l)%t ) )                          &           
              ALLOCATE ( t_green_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )  
           IF ( .NOT. ALLOCATED( m_liq_usm_v_1(l)%var_usm_1d ) )              &
              ALLOCATE ( m_liq_usm_v_1(l)%var_usm_1d(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( m_liq_usm_v_2(l)%var_usm_1d ) )              &
              ALLOCATE ( m_liq_usm_v_2(l)%var_usm_1d(1:surf_usm_v(l)%ns) )
           IF ( .NOT. ALLOCATED( swc_v_1(l)%t ) )                           &           
              ALLOCATE ( swc_v_1(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) ) 
           IF ( .NOT. ALLOCATED( swc_v_2(l)%t ) )                           &           
              ALLOCATE ( swc_v_2(l)%t(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )  
        ENDDO
!
!--     initial assignment of the pointers
        t_wall_v    => t_wall_v_1;    t_wall_v_p    => t_wall_v_2
        t_surf_wall_v => t_surf_wall_v_1; t_surf_wall_v_p => t_surf_wall_v_2
        t_window_v    => t_window_v_1;    t_window_v_p    => t_window_v_2
        t_green_v    => t_green_v_1;    t_green_v_p    => t_green_v_2
        t_surf_window_v => t_surf_window_v_1; t_surf_window_v_p => t_surf_window_v_2
        t_surf_green_v => t_surf_green_v_1; t_surf_green_v_p => t_surf_green_v_2
        m_liq_usm_v     => m_liq_usm_v_1;     m_liq_usm_v_p     => m_liq_usm_v_2
        swc_v    => swc_v_1;    swc_v_p    => swc_v_2

!
!--     Allocate intermediate timestep arrays. For horizontal surfaces.
        ALLOCATE ( surf_usm_h%tt_surface_wall_m(1:surf_usm_h%ns)                  )
        ALLOCATE ( surf_usm_h%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%tt_surface_window_m(1:surf_usm_h%ns)             )
        ALLOCATE ( surf_usm_h%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_h%ns)  )
        ALLOCATE ( surf_usm_h%tt_surface_green_m(1:surf_usm_h%ns)              )

!
!--    Allocate intermediate timestep arrays
!--    Horizontal surfaces
       ALLOCATE ( tm_liq_usm_h_m%var_usm_1d(1:surf_usm_h%ns)                   )
!
!--    Horizontal surfaces
       DO  l = 0, 3
          ALLOCATE ( tm_liq_usm_v_m(l)%var_usm_1d(1:surf_usm_v(l)%ns)          )
       ENDDO 
        
!
!--     Set inital values for prognostic quantities
        IF ( ALLOCATED( surf_usm_h%tt_surface_wall_m ) )  surf_usm_h%tt_surface_wall_m = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%tt_wall_m    ) )  surf_usm_h%tt_wall_m    = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%tt_surface_window_m ) )  surf_usm_h%tt_surface_window_m = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%tt_window_m    )      )  surf_usm_h%tt_window_m         = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%tt_green_m    )       )  surf_usm_h%tt_green_m          = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%tt_surface_green_m )  )  surf_usm_h%tt_surface_green_m  = 0.0_wp
!
!--     Now, for vertical surfaces
        DO  l = 0, 3
           ALLOCATE ( surf_usm_v(l)%tt_surface_wall_m(1:surf_usm_v(l)%ns)                  )
           ALLOCATE ( surf_usm_v(l)%tt_wall_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_wall_m ) )  surf_usm_v(l)%tt_surface_wall_m = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%tt_wall_m    ) )  surf_usm_v(l)%tt_wall_m    = 0.0_wp
           ALLOCATE ( surf_usm_v(l)%tt_surface_window_m(1:surf_usm_v(l)%ns)             )
           ALLOCATE ( surf_usm_v(l)%tt_window_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns) )
           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_window_m ) )  surf_usm_v(l)%tt_surface_window_m = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%tt_window_m  ) )  surf_usm_v(l)%tt_window_m    = 0.0_wp
           ALLOCATE ( surf_usm_v(l)%tt_surface_green_m(1:surf_usm_v(l)%ns)              )
           IF ( ALLOCATED( surf_usm_v(l)%tt_surface_green_m ) )  surf_usm_v(l)%tt_surface_green_m = 0.0_wp
           ALLOCATE ( surf_usm_v(l)%tt_green_m(nzb_wall:nzt_wall+1,1:surf_usm_v(l)%ns)  )
           IF ( ALLOCATED( surf_usm_v(l)%tt_green_m   ) )  surf_usm_v(l)%tt_green_m    = 0.0_wp
        ENDDO

!--     allocate wall heat flux output array and set initial values. For horizontal surfaces
!         ALLOCATE ( surf_usm_h%wshf(1:surf_usm_h%ns)    )  !can be removed 
        ALLOCATE ( surf_usm_h%wshf_eb(1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%wghf_eb(1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%wghf_eb_window(1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%wghf_eb_green(1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%iwghf_eb(1:surf_usm_h%ns) )
        ALLOCATE ( surf_usm_h%iwghf_eb_window(1:surf_usm_h%ns) )
        IF ( ALLOCATED( surf_usm_h%wshf    ) )  surf_usm_h%wshf    = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%wshf_eb ) )  surf_usm_h%wshf_eb = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%wghf_eb ) )  surf_usm_h%wghf_eb = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%wghf_eb_window ) )  surf_usm_h%wghf_eb_window = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%wghf_eb_green ) )  surf_usm_h%wghf_eb_green = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%iwghf_eb ) )  surf_usm_h%iwghf_eb = 0.0_wp
        IF ( ALLOCATED( surf_usm_h%iwghf_eb_window ) )  surf_usm_h%iwghf_eb_window = 0.0_wp
!
!--     Now, for vertical surfaces
        DO  l = 0, 3
!            ALLOCATE ( surf_usm_v(l)%wshf(1:surf_usm_v(l)%ns)    )    ! can be removed
           ALLOCATE ( surf_usm_v(l)%wshf_eb(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%wghf_eb(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%wghf_eb_window(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%wghf_eb_green(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%iwghf_eb(1:surf_usm_v(l)%ns) )
           ALLOCATE ( surf_usm_v(l)%iwghf_eb_window(1:surf_usm_v(l)%ns) )
           IF ( ALLOCATED( surf_usm_v(l)%wshf    ) )  surf_usm_v(l)%wshf    = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%wshf_eb ) )  surf_usm_v(l)%wshf_eb = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb ) )  surf_usm_v(l)%wghf_eb = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_window ) )  surf_usm_v(l)%wghf_eb_window = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%wghf_eb_green ) )  surf_usm_v(l)%wghf_eb_green = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb ) )  surf_usm_v(l)%iwghf_eb = 0.0_wp
           IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) )  surf_usm_v(l)%iwghf_eb_window = 0.0_wp
        ENDDO
        
    END SUBROUTINE usm_allocate_surface


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Sum up and time-average urban surface output quantities as well as allocate
!> the array necessary for storing the average.
!------------------------------------------------------------------------------!
    SUBROUTINE usm_average_3d_data( mode, variable )

        IMPLICIT NONE

        CHARACTER(LEN=*), INTENT(IN) ::  mode
        CHARACTER(LEN=*), INTENT(IN) :: variable
  
        INTEGER(iwp)                                       :: i, j, k, l, m, ids, idsint, iwl, istat
        CHARACTER(LEN=varnamelength)                       :: var
        INTEGER(iwp), PARAMETER                            :: nd = 5
        CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER     :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)

!--     find the real name of the variable
        ids = -1
        l = -1
        var = TRIM(variable)
        DO i = 0, nd-1
            k = len(TRIM(var))
            j = len(TRIM(dirname(i)))
            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                ids = i
                idsint = dirint(ids)
                var = var(:k-j)
                EXIT
            ENDIF
        ENDDO
        l = idsint - 2  ! horisontal direction index - terible hack !
        IF ( l < 0 .OR. l > 3 ) THEN
           l = -1
        END IF
        IF ( ids == -1 )  THEN
            var = TRIM(variable)
        ENDIF
        IF ( var(1:11) == 'usm_t_wall_'  .AND.  len(TRIM(var)) >= 12 )  THEN
!--          wall layers
            READ(var(12:12), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:10)
            ELSE
!--             wrong wall layer index
                RETURN
            ENDIF
        ENDIF
        IF ( var(1:13) == 'usm_t_window_'  .AND.  len(TRIM(var)) >= 14 )  THEN
!--          wall layers
            READ(var(14:14), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:12)
            ELSE
!--             wrong window layer index
                RETURN
            ENDIF
        ENDIF
        IF ( var(1:12) == 'usm_t_green_'  .AND.  len(TRIM(var)) >= 13 )  THEN
!--          wall layers
            READ(var(13:13), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:11)
            ELSE
!--             wrong green layer index
                RETURN
            ENDIF
        ENDIF
        IF ( var(1:8) == 'usm_swc_'  .AND.  len(TRIM(var)) >= 9 )  THEN
!--          swc layers
            READ(var(9:9), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:7)
            ELSE
!--             wrong swc layer index
                RETURN
            ENDIF
        ENDIF

        IF ( mode == 'allocate' )  THEN
           
           SELECT CASE ( TRIM( var ) )

                CASE ( 'usm_wshf' )
!--                 array of sensible heat flux from surfaces
!--                 land surfaces
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
                          ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
                          surf_usm_h%wshf_eb_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%wshf_eb_av = 0.0_wp
                       ENDIF
                    ENDIF
                    
                CASE ( 'usm_qsws' )
!--                 array of latent heat flux from surfaces
!--                 land surfaces
                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%qsws_eb_av) )  THEN
                        ALLOCATE( surf_usm_h%qsws_eb_av(1:surf_usm_h%ns) )
                        surf_usm_h%qsws_eb_av = 0.0_wp
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%qsws_eb_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%qsws_eb_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%qsws_eb_av = 0.0_wp
                       ENDIF
                    ENDIF
                    
                CASE ( 'usm_qsws_veg' )
!--                 array of latent heat flux from vegetation surfaces
!--                 land surfaces
                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%qsws_veg_av) )  THEN
                        ALLOCATE( surf_usm_h%qsws_veg_av(1:surf_usm_h%ns) )
                        surf_usm_h%qsws_veg_av = 0.0_wp
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%qsws_veg_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%qsws_veg_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%qsws_veg_av = 0.0_wp
                       ENDIF
                    ENDIF
                    
                CASE ( 'usm_qsws_liq' )
!--                 array of latent heat flux from surfaces with liquid
!--                 land surfaces
                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%qsws_liq_av) )  THEN
                        ALLOCATE( surf_usm_h%qsws_liq_av(1:surf_usm_h%ns) )
                        surf_usm_h%qsws_liq_av = 0.0_wp
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%qsws_liq_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%qsws_liq_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%qsws_liq_av = 0.0_wp
                       ENDIF
                    ENDIF
!
!--             Please note, the following output quantities belongs to the 
!--             individual tile fractions - ground heat flux at wall-, window-, 
!--             and green fraction. Aggregated ground-heat flux is treated
!--             accordingly in average_3d_data, sum_up_3d_data, etc..
                CASE ( 'usm_wghf' )
!--                 array of heat flux from ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
                           ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
                           surf_usm_h%wghf_eb_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%wghf_eb_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_wghf_window' )
!--                 array of heat flux from window ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
                           ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
                           surf_usm_h%wghf_eb_window_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_wghf_green' )
!--                 array of heat flux from green ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
                           ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
                           surf_usm_h%wghf_eb_green_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_iwghf' )
!--                 array of heat flux from indoor ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
                           ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
                           surf_usm_h%iwghf_eb_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%iwghf_eb_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_iwghf_window' )
!--                 array of heat flux from indoor window ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
                           ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
                           surf_usm_h%iwghf_eb_window_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_t_surf_wall' )
!--                 surface temperature for surfaces
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_wall_av) )  THEN
                           ALLOCATE( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) )
                           surf_usm_h%t_surf_wall_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_wall_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%t_surf_wall_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%t_surf_wall_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_t_surf_window' )
!--                 surface temperature for window surfaces
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
                           ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
                           surf_usm_h%t_surf_window_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%t_surf_window_av = 0.0_wp
                       ENDIF
                    ENDIF
                    
                CASE ( 'usm_t_surf_green' )
!--                 surface temperature for green surfaces
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
                           ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
                           surf_usm_h%t_surf_green_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%t_surf_green_av = 0.0_wp
                       ENDIF
                    ENDIF
                
                CASE ( 'usm_theta_10cm' )
!--                 near surface (10cm) temperature for whole surfaces
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%pt_10cm_av) )  THEN
                           ALLOCATE( surf_usm_h%pt_10cm_av(1:surf_usm_h%ns) )
                           surf_usm_h%pt_10cm_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%pt_10cm_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%pt_10cm_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%pt_10cm_av = 0.0_wp
                       ENDIF
                    ENDIF
                 
                CASE ( 'usm_t_wall' )
!--                 wall temperature for iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
                           ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
                           surf_usm_h%t_wall_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%t_wall_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_t_window' )
!--                 window temperature for iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
                           ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
                           surf_usm_h%t_window_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%t_window_av = 0.0_wp
                       ENDIF
                    ENDIF

                CASE ( 'usm_t_green' )
!--                 green temperature for iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
                           ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
                           surf_usm_h%t_green_av = 0.0_wp
                       ENDIF
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%t_green_av = 0.0_wp
                       ENDIF
                    ENDIF
                CASE ( 'usm_swc' )
!--                 soil water content for iwl layer of walls and land
                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%swc_av) )  THEN
                        ALLOCATE( surf_usm_h%swc_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
                        surf_usm_h%swc_av = 0.0_wp
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%swc_av) )  THEN
                           ALLOCATE( surf_usm_v(l)%swc_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%swc_av = 0.0_wp
                       ENDIF
                    ENDIF

               CASE DEFAULT
                   CONTINUE

           END SELECT

        ELSEIF ( mode == 'sum' )  THEN
           
           SELECT CASE ( TRIM( var ) )

                CASE ( 'usm_wshf' )
!--                 array of sensible heat flux from surfaces (land, roof, wall)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wshf_eb_av(m) =                              &
                                             surf_usm_h%wshf_eb_av(m) +           &
                                             surf_usm_h%wshf_eb(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wshf_eb_av(m) =                        &
                                          surf_usm_v(l)%wshf_eb_av(m) +        &
                                          surf_usm_v(l)%wshf_eb(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_qsws' )
!--                 array of latent heat flux from surfaces (land, roof, wall)
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%qsws_eb_av(m) =                              &
                                          surf_usm_h%qsws_eb_av(m) +           &
                                          surf_usm_h%qsws_eb(m)
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_eb_av(m) =                        &
                                          surf_usm_v(l)%qsws_eb_av(m) +        &
                                          surf_usm_v(l)%qsws_eb(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_qsws_veg' )
!--                 array of latent heat flux from vegetation surfaces (land, roof, wall)
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%qsws_veg_av(m) =                              &
                                          surf_usm_h%qsws_veg_av(m) +           &
                                          surf_usm_h%qsws_veg(m)
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_veg_av(m) =                        &
                                          surf_usm_v(l)%qsws_veg_av(m) +        &
                                          surf_usm_v(l)%qsws_veg(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_qsws_liq' )
!--                 array of latent heat flux from surfaces with liquid (land, roof, wall)
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%qsws_liq_av(m) =                              &
                                          surf_usm_h%qsws_liq_av(m) +           &
                                          surf_usm_h%qsws_liq(m)
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_liq_av(m) =                        &
                                          surf_usm_v(l)%qsws_liq_av(m) +        &
                                          surf_usm_v(l)%qsws_liq(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_wghf' )
!--                 array of heat flux from ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wghf_eb_av(m) =                              &
                                             surf_usm_h%wghf_eb_av(m) +           &
                                             surf_usm_h%wghf_eb(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wghf_eb_av(m) =                        &
                                          surf_usm_v(l)%wghf_eb_av(m) +        &
                                          surf_usm_v(l)%wghf_eb(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_wghf_window' )
!--                 array of heat flux from window ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wghf_eb_window_av(m) =                              &
                                             surf_usm_h%wghf_eb_window_av(m) +           &
                                             surf_usm_h%wghf_eb_window(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
                                          surf_usm_v(l)%wghf_eb_window_av(m) +        &
                                          surf_usm_v(l)%wghf_eb_window(m)
                       ENDDO
                    ENDIF

                CASE ( 'usm_wghf_green' )
!--                 array of heat flux from green ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wghf_eb_green_av(m) =                              &
                                             surf_usm_h%wghf_eb_green_av(m) +           &
                                             surf_usm_h%wghf_eb_green(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
                                          surf_usm_v(l)%wghf_eb_green_av(m) +        &
                                          surf_usm_v(l)%wghf_eb_green(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_iwghf' )
!--                 array of heat flux from indoor ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%iwghf_eb_av(m) =                              &
                                             surf_usm_h%iwghf_eb_av(m) +           &
                                             surf_usm_h%iwghf_eb(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
                                          surf_usm_v(l)%iwghf_eb_av(m) +        &
                                          surf_usm_v(l)%iwghf_eb(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_iwghf_window' )
!--                 array of heat flux from indoor window ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%iwghf_eb_window_av(m) =                              &
                                             surf_usm_h%iwghf_eb_window_av(m) +           &
                                             surf_usm_h%iwghf_eb_window(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
                                          surf_usm_v(l)%iwghf_eb_window_av(m) +        &
                                          surf_usm_v(l)%iwghf_eb_window(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_surf_wall' )
!--                 surface temperature for surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                       surf_usm_h%t_surf_wall_av(m) =                               & 
                                          surf_usm_h%t_surf_wall_av(m) +            &
                                          t_surf_wall_h(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_surf_wall_av(m) =                         &
                                          surf_usm_v(l)%t_surf_wall_av(m) +         &
                                          t_surf_wall_v(l)%t(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_surf_window' )
!--                 surface temperature for window surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_surf_window_av(m) =                               &
                                             surf_usm_h%t_surf_window_av(m) +            &
                                             t_surf_window_h(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_surf_window_av(m) =                         &
                                          surf_usm_v(l)%t_surf_window_av(m) +         &
                                          t_surf_window_v(l)%t(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_surf_green' )
!--                 surface temperature for green surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_surf_green_av(m) =                               &
                                             surf_usm_h%t_surf_green_av(m) +            &
                                             t_surf_green_h(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_surf_green_av(m) =                         &
                                          surf_usm_v(l)%t_surf_green_av(m) +         &
                                          t_surf_green_v(l)%t(m)
                       ENDDO
                    ENDIF
                
                CASE ( 'usm_theta_10cm' )
!--                 near surface temperature for whole surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%pt_10cm_av(m) =                               &
                                             surf_usm_h%pt_10cm_av(m) +            &
                                             surf_usm_h%pt_10cm(m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%pt_10cm_av(m) =                         &
                                          surf_usm_v(l)%pt_10cm_av(m) +         &
                                          surf_usm_v(l)%pt_10cm(m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_wall' )
!--                 wall temperature for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_wall_av(iwl,m) =                           &
                                             surf_usm_h%t_wall_av(iwl,m) +        &
                                             t_wall_h(iwl,m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
                                          surf_usm_v(l)%t_wall_av(iwl,m) +     &
                                          t_wall_v(l)%t(iwl,m)
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_window' )
!--                 window temperature for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_window_av(iwl,m) =                           &
                                             surf_usm_h%t_window_av(iwl,m) +        &
                                             t_window_h(iwl,m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_window_av(iwl,m) =                     &
                                          surf_usm_v(l)%t_window_av(iwl,m) +     &
                                          t_window_v(l)%t(iwl,m)
                       ENDDO
                    ENDIF

                CASE ( 'usm_t_green' )
!--                 green temperature for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_green_av(iwl,m) =                           &
                                             surf_usm_h%t_green_av(iwl,m) +        &
                                             t_green_h(iwl,m)
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_green_av(iwl,m) =                     &
                                          surf_usm_v(l)%t_green_av(iwl,m) +     &
                                          t_green_v(l)%t(iwl,m)
                       ENDDO
                    ENDIF

                CASE ( 'usm_swc' )
!--                 soil water content for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%swc_av(iwl,m) =                           &
                                          surf_usm_h%swc_av(iwl,m) +        &
                                          swc_h(iwl,m)
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%swc_av(iwl,m) =                     &
                                          surf_usm_v(l)%swc_av(iwl,m) +     &
                                          swc_v(l)%t(iwl,m)
                       ENDDO
                    ENDIF

                CASE DEFAULT
                    CONTINUE

           END SELECT

        ELSEIF ( mode == 'average' )  THEN
           
           SELECT CASE ( TRIM( var ) )

                CASE ( 'usm_wshf' )
!--                 array of sensible heat flux from surfaces (land, roof, wall)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wshf_eb_av(m) =                              &
                                             surf_usm_h%wshf_eb_av(m) /           &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wshf_eb_av(m) =                        &
                                          surf_usm_v(l)%wshf_eb_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_qsws' )
!--                 array of latent heat flux from surfaces (land, roof, wall)
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%qsws_eb_av(m) =                              &
                                          surf_usm_h%qsws_eb_av(m) /           &
                                          REAL( average_count_3d, kind=wp )
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_eb_av(m) =                        &
                                          surf_usm_v(l)%qsws_eb_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF

                CASE ( 'usm_qsws_veg' )
!--                 array of latent heat flux from vegetation surfaces (land, roof, wall)
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%qsws_veg_av(m) =                              &
                                          surf_usm_h%qsws_veg_av(m) /           &
                                          REAL( average_count_3d, kind=wp )
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_veg_av(m) =                        &
                                          surf_usm_v(l)%qsws_veg_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_qsws_liq' )
!--                 array of latent heat flux from surfaces with liquid (land, roof, wall)
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%qsws_liq_av(m) =                              &
                                          surf_usm_h%qsws_liq_av(m) /           &
                                          REAL( average_count_3d, kind=wp )
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_liq_av(m) =                        &
                                          surf_usm_v(l)%qsws_liq_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_wghf' )
!--                 array of heat flux from ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wghf_eb_av(m) =                              &
                                             surf_usm_h%wghf_eb_av(m) /           &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wghf_eb_av(m) =                        &
                                          surf_usm_v(l)%wghf_eb_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_wghf_window' )
!--                 array of heat flux from window ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wghf_eb_window_av(m) =                              &
                                             surf_usm_h%wghf_eb_window_av(m) /           &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
                                          surf_usm_v(l)%wghf_eb_window_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF

                CASE ( 'usm_wghf_green' )
!--                 array of heat flux from green ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%wghf_eb_green_av(m) =                              &
                                             surf_usm_h%wghf_eb_green_av(m) /           &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
                                          surf_usm_v(l)%wghf_eb_green_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF

                CASE ( 'usm_iwghf' )
!--                 array of heat flux from indoor ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%iwghf_eb_av(m) =                              &
                                             surf_usm_h%iwghf_eb_av(m) /           &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
                                          surf_usm_v(l)%iwghf_eb_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_iwghf_window' )
!--                 array of heat flux from indoor window ground (wall, roof, land)
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%iwghf_eb_window_av(m) =                              &
                                             surf_usm_h%iwghf_eb_window_av(m) /           &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
                                          surf_usm_v(l)%iwghf_eb_window_av(m) /        &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_surf_wall' )
!--                 surface temperature for surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                       surf_usm_h%t_surf_wall_av(m) =                               & 
                                          surf_usm_h%t_surf_wall_av(m) /            &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_surf_wall_av(m) =                         &
                                          surf_usm_v(l)%t_surf_wall_av(m) /         &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_surf_window' )
!--                 surface temperature for window surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_surf_window_av(m) =                               &
                                             surf_usm_h%t_surf_window_av(m) /            &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_surf_window_av(m) =                         &
                                          surf_usm_v(l)%t_surf_window_av(m) /         &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_t_surf_green' )
!--                 surface temperature for green surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_surf_green_av(m) =                               &
                                             surf_usm_h%t_surf_green_av(m) /            &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_surf_green_av(m) =                         &
                                          surf_usm_v(l)%t_surf_green_av(m) /         &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_theta_10cm' )
!--                 near surface temperature for whole surfaces
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%pt_10cm_av(m) =                               &
                                             surf_usm_h%pt_10cm_av(m) /            &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%pt_10cm_av(m) =                         &
                                          surf_usm_v(l)%pt_10cm_av(m) /         &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF

                    
                CASE ( 'usm_t_wall' )
!--                 wall temperature for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_wall_av(iwl,m) =                           &
                                             surf_usm_h%t_wall_av(iwl,m) /        &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
                                          surf_usm_v(l)%t_wall_av(iwl,m) /     &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF

                CASE ( 'usm_t_window' )
!--                 window temperature for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_window_av(iwl,m) =                           &
                                             surf_usm_h%t_window_av(iwl,m) /        &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_window_av(iwl,m) =                     &
                                          surf_usm_v(l)%t_window_av(iwl,m) /     &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF

                CASE ( 'usm_t_green' )
!--                 green temperature for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                       DO  m = 1, surf_usm_h%ns
                          surf_usm_h%t_green_av(iwl,m) =                           &
                                             surf_usm_h%t_green_av(iwl,m) /        &
                                             REAL( average_count_3d, kind=wp )
                       ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%t_green_av(iwl,m) =                     &
                                          surf_usm_v(l)%t_green_av(iwl,m) /     &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF
                    
                CASE ( 'usm_swc' )
!--                 soil water content for  iwl layer of walls and land
                    IF ( l == -1 ) THEN
                    DO  m = 1, surf_usm_h%ns
                       surf_usm_h%swc_av(iwl,m) =                           &
                                          surf_usm_h%swc_av(iwl,m) /        &
                                          REAL( average_count_3d, kind=wp )
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%swc_av(iwl,m) =                     &
                                          surf_usm_v(l)%swc_av(iwl,m) /     &
                                          REAL( average_count_3d, kind=wp )
                       ENDDO
                    ENDIF


           END SELECT

        ENDIF

    END SUBROUTINE usm_average_3d_data



!------------------------------------------------------------------------------!
! Description:
! ------------
!> Set internal Neumann boundary condition at outer soil grid points 
!> for temperature and humidity. 
!------------------------------------------------------------------------------!
 SUBROUTINE usm_boundary_condition
 
    IMPLICIT NONE

    INTEGER(iwp) :: i      !< grid index x-direction
    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
    INTEGER(iwp) :: j      !< grid index y-direction
    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
    INTEGER(iwp) :: k      !< grid index z-direction
    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
    INTEGER(iwp) :: l      !< running index surface-orientation
    INTEGER(iwp) :: m      !< running index surface elements

    koff = surf_usm_h%koff
    DO  m = 1, surf_usm_h%ns
       i = surf_usm_h%i(m)
       j = surf_usm_h%j(m)
       k = surf_usm_h%k(m)
       pt(k+koff,j,i) = pt(k,j,i)
    ENDDO

    DO  l = 0, 3
       ioff = surf_usm_v(l)%ioff
       joff = surf_usm_v(l)%joff
       DO  m = 1, surf_usm_v(l)%ns
          i = surf_usm_v(l)%i(m)
          j = surf_usm_v(l)%j(m)
          k = surf_usm_v(l)%k(m)
          pt(k,j+joff,i+ioff) = pt(k,j,i)
       ENDDO
    ENDDO

 END SUBROUTINE usm_boundary_condition


!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Subroutine checks variables and assigns units.
!> It is called out from subroutine check_parameters.
!------------------------------------------------------------------------------!
    SUBROUTINE usm_check_data_output( variable, unit )

        IMPLICIT NONE

        CHARACTER(LEN=*),INTENT(IN)    ::  variable   !< 
        CHARACTER(LEN=*),INTENT(OUT)   ::  unit       !<

        INTEGER(iwp)                                  :: i,j,l        !< index
        CHARACTER(LEN=2)                              :: ls
        CHARACTER(LEN=varnamelength)                  :: var          !< TRIM(variable)
        INTEGER(iwp), PARAMETER                       :: nl1 = 16     !< number of directional usm variables
        CHARACTER(LEN=varnamelength), DIMENSION(nl1)  :: varlist1 = & !< list of directional usm variables
                  (/'usm_wshf                      ', &
                    'usm_wghf                      ', &
                    'usm_wghf_window               ', &
                    'usm_wghf_green                ', &
                    'usm_iwghf                     ', &
                    'usm_iwghf_window              ', &
                    'usm_surfz                     ', &
                    'usm_surfwintrans              ', &
                    'usm_surfcat                   ', &
                    'usm_surfalb                   ', &
                    'usm_surfemis                  ', &
                    'usm_t_surf_wall               ', &
                    'usm_t_surf_window             ', &
                    'usm_t_surf_green              ', &
                    'usm_t_green                   ', &
                    'usm_theta_10cm                '/)

        INTEGER(iwp), PARAMETER                       :: nl2 = 3      !< number of directional layer usm variables
        CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of directional layer usm variables
                  (/'usm_t_wall                    ', &
                    'usm_t_window                  ', &
                    'usm_t_green                   '/)

        INTEGER(iwp), PARAMETER                       :: nd = 5     !< number of directions
        CHARACTER(LEN=6), DIMENSION(nd), PARAMETER  :: dirname = &  !< direction names
                  (/'_roof ','_south','_north','_west ','_east '/)
        LOGICAL                                       :: lfound     !< flag if the variable is found


        lfound = .FALSE.

        var = TRIM(variable)

!--     check if variable exists
!       directional variables
        DO i = 1, nl1
           DO j = 1, nd
              IF ( TRIM(var) == TRIM(varlist1(i))//TRIM(dirname(j)) ) THEN
                 lfound = .TRUE.
                 EXIT
              ENDIF
              IF ( lfound ) EXIT
           ENDDO
        ENDDO
        IF ( lfound ) GOTO 10
!       directional layer variables
        DO i = 1, nl2
           DO j = 1, nd
              DO l = nzb_wall, nzt_wall
                 WRITE(ls,'(A1,I1)') '_',l
                 IF ( TRIM(var) == TRIM(varlist2(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
                    lfound = .TRUE.
                    EXIT
                 ENDIF
              ENDDO
              IF ( lfound ) EXIT
           ENDDO
        ENDDO
        IF ( .NOT.  lfound ) THEN
           unit = 'illegal'
           RETURN
        ENDIF
10      CONTINUE

        IF ( var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
             var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_'    .OR.      &
             var(1:17) == 'usm_surfwintrans_' .OR.                                        &
             var(1:9)  == 'usm_qsws_'  .OR.  var(1:13)  == 'usm_qsws_veg_'  .OR.          &
             var(1:13) == 'usm_qsws_liq_' ) THEN
            unit = 'W/m2'
        ELSE IF ( var(1:15) == 'usm_t_surf_wall'   .OR.  var(1:10) == 'usm_t_wall' .OR.   &
                  var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR.  &
                  var(1:16) == 'usm_t_surf_green'  .OR.                                   &
                  var(1:11) == 'usm_t_green' .OR.  var(1:7) == 'usm_swc' .OR.             &
                  var(1:14) == 'usm_theta_10cm' )  THEN
            unit = 'K'
        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:11) == 'usm_surfcat'  .OR.         &
                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  )  THEN
            unit = '1'
        ELSE
            unit = 'illegal'
        ENDIF

    END SUBROUTINE usm_check_data_output


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check parameters routine for urban surface model
!------------------------------------------------------------------------------!
    SUBROUTINE usm_check_parameters
    
       USE control_parameters,                                                 &
           ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing,   &
                  lsf_surf, topography

!
!--    Dirichlet boundary conditions are required as the surface fluxes are
!--    calculated from the temperature/humidity gradients in the urban surface
!--    model
       IF ( bc_pt_b == 'neumann'   .OR.   bc_q_b == 'neumann' )  THEN
          message_string = 'urban surface model requires setting of '//        &
                           'bc_pt_b = "dirichlet" and '//                      &
                           'bc_q_b  = "dirichlet"'
          CALL message( 'usm_check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
       ENDIF

       IF ( .NOT.  constant_flux_layer )  THEN
          message_string = 'urban surface model requires '//                   &
                           'constant_flux_layer = .T.'
          CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
       ENDIF

       IF (  .NOT.  radiation )  THEN
          message_string = 'urban surface model requires '//                   &
                           'the radiation model to be switched on'
          CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
       ENDIF
!        
!--    Surface forcing has to be disabled for LSF in case of enabled 
!--    urban surface module
       IF ( large_scale_forcing )  THEN
          lsf_surf = .FALSE.
       ENDIF
!
!--    Topography
       IF ( topography == 'flat' )  THEN
          message_string = 'topography /= "flat" is required '//               &
                           'when using the urban surface model'
          CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
       ENDIF
!
!--    naheatlayers
       IF ( naheatlayers > nzt )  THEN
          message_string = 'number of anthropogenic heat layers '//            &
                           '"naheatlayers" can not be larger than'//           &
                           ' number of domain layers "nzt"'
          CALL message( 'check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
       ENDIF

    END SUBROUTINE usm_check_parameters


!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Output of the 3D-arrays in netCDF and/or AVS format 
!> for variables of urban_surface model.
!> It resorts the urban surface module output quantities from surf style
!> indexing into temporary 3D array with indices (i,j,k).
!> It is called from subroutine data_output_3d.
!------------------------------------------------------------------------------!
    SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
        
        IMPLICIT NONE

        INTEGER(iwp), INTENT(IN)       ::  av        !< 
        CHARACTER (len=*), INTENT(IN)  ::  variable  !< 
        INTEGER(iwp), INTENT(IN)       ::  nzb_do    !< lower limit of the data output (usually 0)
        INTEGER(iwp), INTENT(IN)       ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
        LOGICAL, INTENT(OUT)           ::  found     !< 
        REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< sp - it has to correspond to module data_output_3d
        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)     ::  temp_pf    !< temp array for urban surface output procedure
        
        CHARACTER (len=varnamelength)                          :: var
        INTEGER(iwp), PARAMETER                                :: nd = 5
        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER         :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: dirint =  (/    iup_u, isouth_u, inorth_u,  iwest_u,  ieast_u /)
        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: diridx =  (/       -1,        1,        0,        3,        2 /)
                                                                     !< index for surf_*_v: 0:3 = (North, South, East, West)
        INTEGER(iwp)                                           :: ids,idsint,idsidx,isvf
        INTEGER(iwp)                                           :: i,j,k,iwl,istat, l, m

        found = .TRUE.
        temp_pf = -1._wp
        
        ids = -1
        var = TRIM(variable)
        DO i = 0, nd-1
            k = len(TRIM(var))
            j = len(TRIM(dirname(i)))
            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                ids = i
                idsint = dirint(ids)
                idsidx = diridx(ids)
                var = var(:k-j)
                EXIT
            ENDIF
        ENDDO
        IF ( ids == -1 )  THEN
            var = TRIM(variable)
        ENDIF
        IF ( var(1:11) == 'usm_t_wall_'  .AND.  len(TRIM(var)) >= 12 )  THEN
!--         wall layers
            READ(var(12:12), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:10)
            ENDIF
        ENDIF
        IF ( var(1:13) == 'usm_t_window_'  .AND.  len(TRIM(var)) >= 14 )  THEN
!--         window layers
            READ(var(14:14), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:12)
            ENDIF
        ENDIF
        IF ( var(1:12) == 'usm_t_green_'  .AND.  len(TRIM(var)) >= 13 )  THEN
!--         green layers
            READ(var(13:13), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:11)
            ENDIF
        ENDIF
        IF ( var(1:8) == 'usm_swc_'  .AND.  len(TRIM(var)) >= 9 )  THEN
!--         green layers soil water content
            READ(var(9:9), '(I1)', iostat=istat ) iwl
            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
                var = var(1:7)
            ENDIF
        ENDIF
        
        SELECT CASE ( TRIM(var) )

          CASE ( 'usm_surfz' )
!--           array of surface height (z)
              IF ( idsint == iup_u )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
                    j = surf_usm_h%j(m)
                    k = surf_usm_h%k(m)
                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
                 ENDDO
              ELSE
                 l = idsidx
                 DO  m = 1, surf_usm_v(l)%ns
                    i = surf_usm_v(l)%i(m)
                    j = surf_usm_v(l)%j(m)
                    k = surf_usm_v(l)%k(m)
                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
                 ENDDO
              ENDIF

          CASE ( 'usm_surfcat' )
!--           surface category
              IF ( idsint == iup_u )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
                    j = surf_usm_h%j(m)
                    k = surf_usm_h%k(m)
                    temp_pf(k,j,i) = surf_usm_h%surface_types(m)
                 ENDDO
              ELSE
                 l = idsidx
                 DO  m = 1, surf_usm_v(l)%ns
                    i = surf_usm_v(l)%i(m)
                    j = surf_usm_v(l)%j(m)
                    k = surf_usm_v(l)%k(m)
                    temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
                 ENDDO
              ENDIF
              
          CASE ( 'usm_surfalb' )
!--           surface albedo, weighted average
              IF ( idsint == iup_u )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
                    j = surf_usm_h%j(m)
                    k = surf_usm_h%k(m)
                    temp_pf(k,j,i) = surf_usm_h%frac(ind_veg_wall,m)     *     &
                                     surf_usm_h%albedo(ind_veg_wall,m)  +      &
                                     surf_usm_h%frac(ind_pav_green,m)    *     &
                                     surf_usm_h%albedo(ind_pav_green,m) +      &
                                     surf_usm_h%frac(ind_wat_win,m)      *     &
                                     surf_usm_h%albedo(ind_wat_win,m)
                 ENDDO
              ELSE
                 l = idsidx
                 DO  m = 1, surf_usm_v(l)%ns
                    i = surf_usm_v(l)%i(m)
                    j = surf_usm_v(l)%j(m)
                    k = surf_usm_v(l)%k(m)
                    temp_pf(k,j,i) = surf_usm_v(l)%frac(ind_veg_wall,m)     *  &
                                     surf_usm_v(l)%albedo(ind_veg_wall,m)  +   &
                                     surf_usm_v(l)%frac(ind_pav_green,m)    *  &
                                     surf_usm_v(l)%albedo(ind_pav_green,m) +   &
                                     surf_usm_v(l)%frac(ind_wat_win,m)      *  &
                                     surf_usm_v(l)%albedo(ind_wat_win,m)
                 ENDDO
              ENDIF
              
          CASE ( 'usm_surfemis' )
!--           surface emissivity, weighted average
              IF ( idsint == iup_u )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
                    j = surf_usm_h%j(m)
                    k = surf_usm_h%k(m)
                    temp_pf(k,j,i) =  surf_usm_h%frac(ind_veg_wall,m)      *   &
                                      surf_usm_h%emissivity(ind_veg_wall,m)  + &
                                      surf_usm_h%frac(ind_pav_green,m)     *   &
                                      surf_usm_h%emissivity(ind_pav_green,m) + &
                                      surf_usm_h%frac(ind_wat_win,m)       *   &
                                      surf_usm_h%emissivity(ind_wat_win,m)
                 ENDDO
              ELSE
                 l = idsidx
                 DO  m = 1, surf_usm_v(l)%ns
                    i = surf_usm_v(l)%i(m)
                    j = surf_usm_v(l)%j(m)
                    k = surf_usm_v(l)%k(m)
                    temp_pf(k,j,i) = surf_usm_v(l)%frac(ind_veg_wall,m)       *&
                                     surf_usm_v(l)%emissivity(ind_veg_wall,m) +&
                                     surf_usm_v(l)%frac(ind_pav_green,m)      *&
                                     surf_usm_v(l)%emissivity(ind_pav_green,m)+&
                                     surf_usm_v(l)%frac(ind_wat_win,m)        *&
                                     surf_usm_v(l)%emissivity(ind_wat_win,m)
                 ENDDO
              ENDIF

          CASE ( 'usm_surfwintrans' )
!--           transmissivity window tiles
              IF ( idsint == iup_u )  THEN
                 DO  m = 1, surf_usm_h%ns
                    i = surf_usm_h%i(m)
                    j = surf_usm_h%j(m)
                    k = surf_usm_h%k(m)
                    temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
                 ENDDO
              ELSE
                 l = idsidx
                 DO  m = 1, surf_usm_v(l)%ns
                    i = surf_usm_v(l)%i(m)
                    j = surf_usm_v(l)%j(m)
                    k = surf_usm_v(l)%k(m)
                    temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
                 ENDDO
              ENDIF

          CASE ( 'usm_wshf' )
!--           array of sensible heat flux from surfaces
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
                    ENDDO
                 ENDIF
              ENDIF
              
              
          CASE ( 'usm_qsws' )
!--           array of latent heat flux from surfaces
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%qsws_eb(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%qsws_eb(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%qsws_eb_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%qsws_eb_av(m)
                    ENDDO
                 ENDIF
              ENDIF
              
          CASE ( 'usm_qsws_veg' )
!--           array of latent heat flux from vegetation surfaces
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%qsws_veg(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%qsws_veg(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%qsws_veg_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%qsws_veg_av(m)
                    ENDDO
                 ENDIF
              ENDIF
              
          CASE ( 'usm_qsws_liq' )
!--           array of latent heat flux from surfaces with liquid
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%qsws_liq(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%qsws_liq(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%qsws_liq_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%qsws_liq_av(m)
                    ENDDO
                 ENDIF
              ENDIF

          CASE ( 'usm_wghf' )
!--           array of heat flux from ground (land, wall, roof)
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
                    ENDDO
                 ENDIF
              ENDIF

          CASE ( 'usm_wghf_window' )
!--           array of heat flux from window ground (land, wall, roof)

              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
                    ENDDO
                 ENDIF
              ENDIF

          CASE ( 'usm_wghf_green' )
!--           array of heat flux from green ground (land, wall, roof)

              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green_av(m)
                    ENDDO
                 ENDIF
              ENDIF

          CASE ( 'usm_iwghf' )
!--           array of heat flux from indoor ground (land, wall, roof)
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_av(m)
                    ENDDO
                 ENDIF
              ENDIF

          CASE ( 'usm_iwghf_window' )
!--           array of heat flux from indoor window ground (land, wall, roof)

              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%iwghf_eb_window_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%iwghf_eb_window_av(m)
                    ENDDO
                 ENDIF
              ENDIF
              
          CASE ( 'usm_t_surf_wall' )
!--           surface temperature for surfaces
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = t_surf_wall_h(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = t_surf_wall_v(l)%t(m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%t_surf_wall_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_wall_av(m)
                    ENDDO
                 ENDIF
              ENDIF
              
          CASE ( 'usm_t_surf_window' )
!--           surface temperature for window surfaces

              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = t_surf_window_h(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = t_surf_window_v(l)%t(m)
                    ENDDO
                 ENDIF

              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%t_surf_window_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_window_av(m)
                    ENDDO

                 ENDIF

              ENDIF

          CASE ( 'usm_t_surf_green' )
!--           surface temperature for green surfaces

              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = t_surf_green_h(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = t_surf_green_v(l)%t(m)
                    ENDDO
                 ENDIF

              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%t_surf_green_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%t_surf_green_av(m)
                    ENDDO

                 ENDIF

              ENDIF

          CASE ( 'usm_theta_10cm' )
!--           near surface temperature for whole surfaces

              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%pt_10cm(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%pt_10cm(m)
                    ENDDO
                 ENDIF
              
              
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%pt_10cm_av(m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%pt_10cm_av(m)
                    ENDDO

                  ENDIF
              ENDIF
             
          CASE ( 'usm_t_wall' )
!--           wall temperature for  iwl layer of walls and land
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = t_wall_h(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = t_wall_v(l)%t(iwl,m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%t_wall_av(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%t_wall_av(iwl,m)
                    ENDDO
                 ENDIF
              ENDIF
             
          CASE ( 'usm_t_window' )
!--           window temperature for iwl layer of walls and land
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = t_window_h(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = t_window_v(l)%t(iwl,m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%t_window_av(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%t_window_av(iwl,m)
                    ENDDO
                 ENDIF
              ENDIF

          CASE ( 'usm_t_green' )
!--           green temperature for  iwl layer of walls and land
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = t_green_h(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = t_green_v(l)%t(iwl,m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%t_green_av(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%t_green_av(iwl,m)
                    ENDDO
                 ENDIF
              ENDIF
              
              CASE ( 'usm_swc' )
!--           soil water content for  iwl layer of walls and land
              IF ( av == 0 )  THEN
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = swc_h(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = swc_v(l)%t(iwl,m)
                    ENDDO
                 ENDIF
              ELSE
                 IF ( idsint == iup_u )  THEN
                    DO  m = 1, surf_usm_h%ns
                       i = surf_usm_h%i(m)
                       j = surf_usm_h%j(m)
                       k = surf_usm_h%k(m)
                       temp_pf(k,j,i) = surf_usm_h%swc_av(iwl,m)
                    ENDDO
                 ELSE
                    l = idsidx
                    DO  m = 1, surf_usm_v(l)%ns
                       i = surf_usm_v(l)%i(m)
                       j = surf_usm_v(l)%j(m)
                       k = surf_usm_v(l)%k(m)
                       temp_pf(k,j,i) = surf_usm_v(l)%swc_av(iwl,m)
                    ENDDO
                 ENDIF
              ENDIF

             
          CASE DEFAULT
              found = .FALSE.
              RETURN
        END SELECT

!
!--     Rearrange dimensions for NetCDF output
!--     FIXME: this may generate FPE overflow upon conversion from DP to SP
        DO  j = nys, nyn
            DO  i = nxl, nxr