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

!> @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-2020 Czech Technical University in Prague
! Copyright 2015-2020 Institute of Computer Science of the
!                     Czech Academy of Sciences, Prague
! Copyright 1997-2020 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: urban_surface_mod.f90 4495 2020-04-13 20:11:20Z raasch $
! restart data handling with MPI-IO added
! 
! 4493 2020-04-10 09:49:43Z pavelkrc
! J.Resler, 2020/03/19
! - remove reading of deprecated input parameters c_surface and lambda_surf
! - and calculate them from parameters of the outer wall/roof layer
! 
! 4481 2020-03-31 18:55:54Z maronga
! use statement for exchange horiz added
! 
! 4442 2020-03-04 19:21:13Z suehring
! Change order of dimension in surface arrays %frac, %emissivity and %albedo 
! to allow for better vectorization in the radiation interactions.
! 
! 4441 2020-03-04 19:20:35Z suehring
! Removed wall_flags_static_0 from USE statements as it's not used within
! the module
! 
! 4329 2019-12-10 15:46:36Z motisi
! Renamed wall_flags_0 to wall_flags_static_0
! 
! 4309 2019-11-26 18:49:59Z suehring
! - Bugfix, include m_liq into restarts
! - Remove unused arrays for liquid water and saturation moisture at vertical
!   walls
!
! 4305 2019-11-25 11:15:40Z suehring
! Revision of some indoor-model parameters
! 
! 4259 2019-10-09 10:05:22Z suehring
! Instead of terminate the job in case the relative wall fractions do not 
! sum-up to one, give only an informative message and normalize the fractions.
! 
! 4258 2019-10-07 13:29:08Z suehring
! - Add checks to ensure that relative fractions of walls, windowns and green
!   surfaces sum-up to one.
! - Revise message calls dealing with local checks.
! 
! 4245 2019-09-30 08:40:37Z pavelkrc
! Initialize explicit per-surface parameters from building_surface_pars
! 
! 4238 2019-09-25 16:06:01Z suehring
! Indoor-model parameters for some building types adjusted in order to avoid
! unrealistically high indoor temperatures (S. Rissmann) 
! 
! 4230 2019-09-11 13:58:14Z suehring
! Bugfix, initialize canopy resistance. Even if no green fraction is set, 
! r_canopy must be initialized for output purposes. 
! 
! 4227 2019-09-10 18:04:34Z gronemeier
! implement new palm_date_time_mod
! 
! 4214 2019-09-02 15:57:02Z suehring
! Bugfix, missing initialization and clearing of soil-moisture tendency 
! (J.Resler)
! 
! 4182 2019-08-22 15:20:23Z scharf
! Corrected "Former revisions" section
! 
! 4168 2019-08-16 13:50:17Z suehring
! Replace function get_topography_top_index by topo_top_ind
! 
! 4148 2019-08-08 11:26:00Z suehring
! - Add anthropogenic heat output factors for heating and cooling to building
!   data base
! - Move definition of building_pars to usm_init_arrays since it is already
!   required in the indoor model
! 
! 4127 2019-07-30 14:47:10Z suehring
! Do not add anthopogenic energy during wall/soil spin-up 
! (merge from branch resler)
! 
! 4077 2019-07-09 13:27:11Z gronemeier
! Set roughness length z0 and z0h/q at ground-floor level to same value as
! those above ground-floor level
! 
! 4051 2019-06-24 13:58:30Z suehring
! Remove work-around for green surface fraction on buildings 
! (do not set it zero) 
! 
! 4050 2019-06-24 13:57:27Z suehring
! In order to avoid confusion with global control parameter, rename the 
! USM-internal flag spinup into during_spinup.
! 
! 3987 2019-05-22 09:52:13Z kanani
! Introduce alternative switch for debug output during timestepping
! 
! 3943 2019-05-02 09:50:41Z maronga
! Removed qsws_eb. Bugfix in calculation of qsws.
! 
! 3933 2019-04-25 12:33:20Z kanani
! Remove allocation of pt_2m, this is done in surface_mod now (surfaces%pt_2m)
! 
! 3921 2019-04-18 14:21:10Z suehring
! Undo accidentally commented initialization  
! 
! 3918 2019-04-18 13:33:11Z suehring
! Set green fraction to zero also at vertical surfaces
! 
! 3914 2019-04-17 16:02:02Z suehring
! In order to obtain correct surface temperature during spinup set window 
! fraction to zero (only during spinup) instead of just disabling
! time-integration of window-surface temperature. 
! 
! 3901 2019-04-16 16:17:02Z suehring
! Workaround - set green fraction to zero ( green-heat model crashes ).
! 
! 3896 2019-04-15 10:10:17Z suehring
! 
! 
! 3896 2019-04-15 10:10:17Z suehring
! Bugfix, wrong index used for accessing building_pars from PIDS 
! 
! 3885 2019-04-11 11:29:34Z kanani
! Changes related to global restructuring of location messages and introduction 
! of additional debug messages
! 
! 3882 2019-04-10 11:08:06Z suehring
! Avoid different type kinds
! Move definition of building-surface properties from declaration block
! to an extra routine
! 
! 3881 2019-04-10 09:31:22Z suehring
! Revise determination of local ground-floor level height. 
! Make level 3 initalization conform with Palm-input-data standard
! Move output of albedo and emissivity to radiation module
! 
! 3832 2019-03-28 13:16:58Z raasch
! instrumented with openmp directives
! 
! 3824 2019-03-27 15:56:16Z pavelkrc
! Remove unused imports
!
! 
! 3814 2019-03-26 08:40:31Z pavelkrc
! unused subroutine commented out
! 
! 3769 2019-02-28 10:16:49Z moh.hefny
! removed unused variables
!
! 3767 2019-02-27 08:18:02Z raasch
! unused variables removed from rrd-subroutines parameter list
! 
! 3748 2019-02-18 10:38:31Z suehring
! Revise conversion of waste-heat flux (do not divide by air density, will
! be done in diffusion_s)
! 
! 3745 2019-02-15 18:57:56Z suehring
! - Remove internal flag indoor_model (is a global control parameter)
! - add waste heat from buildings to the kinmatic heat flux
! - consider waste heat in restart data
! - remove unused USE statements
! 
! 3744 2019-02-15 18:38:58Z suehring
! fixed surface heat capacity in the building parameters
! convert the file back to unix format
!
! 3730 2019-02-11 11:26:47Z moh.hefny
! Formatting and clean-up (rvtils)
! 
! 3710 2019-01-30 18:11:19Z suehring
! Check if building type is set within a valid range.
! 
! 3705 2019-01-29 19:56:39Z suehring
! make nzb_wall public, required for virtual-measurements
! 
! 3704 2019-01-29 19:51:41Z suehring
! Some interface calls moved to module_interface + cleanup
! 
! 3655 2019-01-07 16:51:22Z knoop
! Implementation of the PALM module interface
! 
! 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:
!> -------------
!> @todo Output of _av variables in case of restarts
!> @todo Revise flux conversion in energy-balance solver
!> @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 remove reading of old csv inputs
!------------------------------------------------------------------------------!
 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, sigma_sb

    USE control_parameters,                                                    &
        ONLY:  coupling_start_time, topography,                                &
               debug_output, debug_output_timestep, debug_string,              &
               dt_3d, humidity, indoor_model,                                  &
               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_pt_mean, spinup_time, time_do3d, dt_do3d,                &
               average_count_3d, varnamelength, urban_surface, dz, restart_data_format_output

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

    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, topo_top_ind

    USE, INTRINSIC :: iso_c_binding 

    USE kinds
              
    USE palm_date_time_mod,                                                    &
        ONLY:  get_date_time, seconds_per_hour

    USE pegrid
        
    USE radiation_model_mod,                                                   &
        ONLY:  albedo_type, radiation_interaction,                             &
               radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
               force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u,       &
               iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id,       &
               nz_urban_b, nz_urban_t, unscheduled_radiation_calls

    USE restart_data_mpi_io_mod,                                                                   &
        ONLY:  rd_mpi_io_surface_filetypes, wrd_mpi_io, wrd_mpi_io_surface

    USE statistics,                                                            &
        ONLY:  hom, statistic_regions

    USE surface_mod,                                                           &
        ONLY:  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, &  !< soil 1
                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, &  !< soil 2
                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, &  !< soil 3
                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, &  !< soil 4
                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, &  !< soil 5
                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, &  !< soil 6
                                 0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  &  !< soil 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, &  !< soil 1
                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, &  !< soil 2
                                 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, &  !< soil 3
                                 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, &  !< soil 4
                                 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, &  !< soil 5
                                 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, &  !< soil 6
                                 0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  &  !< soil 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
                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 ::  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 in building_pars for (above) ground floor level 
    INTEGER(iwp) ::  ind_alb_wall_agfl     = 38   !< index in input list for albedo_type of wall above ground floor level
    INTEGER(iwp) ::  ind_alb_wall_gfl      = 66   !< index in input list for albedo_type of wall ground floor level
    INTEGER(iwp) ::  ind_alb_wall_r        = 101  !< index in input list for albedo_type of wall roof
    INTEGER(iwp) ::  ind_alb_green_agfl    = 39   !< index in input list for albedo_type of green above ground floor level
    INTEGER(iwp) ::  ind_alb_green_gfl     = 78   !< index in input list for albedo_type of green ground floor level
    INTEGER(iwp) ::  ind_alb_green_r       = 117  !< index in input list for albedo_type of green roof
    INTEGER(iwp) ::  ind_alb_win_agfl      = 40   !< index in input list for albedo_type of window fraction above ground floor level
    INTEGER(iwp) ::  ind_alb_win_gfl       = 77   !< index in input list for albedo_type of window fraction ground floor level
    INTEGER(iwp) ::  ind_alb_win_r         = 115  !< index in input list for albedo_type of window fraction roof
    INTEGER(iwp) ::  ind_emis_wall_agfl    = 14   !< index in input list for wall emissivity, above ground floor level
    INTEGER(iwp) ::  ind_emis_wall_gfl     = 32   !< index in input list for wall emissivity, ground floor level
    INTEGER(iwp) ::  ind_emis_wall_r       = 100  !< index in input list for wall emissivity, roof
    INTEGER(iwp) ::  ind_emis_green_agfl   = 15   !< index in input list for green emissivity, above ground floor level
    INTEGER(iwp) ::  ind_emis_green_gfl    = 34   !< index in input list for green emissivity, ground floor level
    INTEGER(iwp) ::  ind_emis_green_r      = 116  !< index in input list for green emissivity, roof
    INTEGER(iwp) ::  ind_emis_win_agfl     = 16   !< index in input list for window emissivity, above ground floor level
    INTEGER(iwp) ::  ind_emis_win_gfl      = 33   !< index in input list for window emissivity, ground floor level
    INTEGER(iwp) ::  ind_emis_win_r        = 113  !< index in input list for window emissivity, roof
    INTEGER(iwp) ::  ind_gflh              = 20   !< index in input list for ground floor level height
    INTEGER(iwp) ::  ind_green_frac_w_agfl = 2    !< index in input list for green fraction on wall, above ground floor level
    INTEGER(iwp) ::  ind_green_frac_w_gfl  = 23   !< index in input list for green fraction on wall, ground floor level
    INTEGER(iwp) ::  ind_green_frac_r_agfl = 3    !< index in input list for green fraction on roof, above ground floor level
    INTEGER(iwp) ::  ind_green_frac_r_gfl  = 24   !< index in input list for green fraction on roof, ground floor level
    INTEGER(iwp) ::  ind_hc1_agfl          = 6    !< index in input list for heat capacity at first wall layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_hc1_gfl           = 26   !< index in input list for heat capacity at first wall layer, ground floor level
    INTEGER(iwp) ::  ind_hc1_wall_r        = 94   !< index in input list for heat capacity at first wall layer, roof
    INTEGER(iwp) ::  ind_hc1_win_agfl      = 83   !< index in input list for heat capacity at first window layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_hc1_win_gfl       = 71   !< index in input list for heat capacity at first window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_hc1_win_r         = 107  !< index in input list for heat capacity at first window layer, roof
    INTEGER(iwp) ::  ind_hc2_agfl          = 7    !< index in input list for heat capacity at second wall layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_hc2_gfl           = 27   !< index in input list for heat capacity at second wall layer, ground floor level
    INTEGER(iwp) ::  ind_hc2_wall_r        = 95   !< index in input list for heat capacity at second wall layer, roof
    INTEGER(iwp) ::  ind_hc2_win_agfl      = 84   !< index in input list for heat capacity at second window layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_hc2_win_gfl       = 72   !< index in input list for heat capacity at second window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_hc2_win_r         = 108  !< index in input list for heat capacity at second window layer, roof
    INTEGER(iwp) ::  ind_hc3_agfl          = 8    !< index in input list for heat capacity at third wall layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_hc3_gfl           = 28   !< index in input list for heat capacity at third wall layer, ground floor level
    INTEGER(iwp) ::  ind_hc3_wall_r        = 96   !< index in input list for heat capacity at third wall layer, roof
    INTEGER(iwp) ::  ind_hc3_win_agfl      = 85   !< index in input list for heat capacity at third window layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_hc3_win_gfl       = 73   !< index in input list for heat capacity at third window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_hc3_win_r         = 109  !< index in input list for heat capacity at third window layer, roof
    INTEGER(iwp) ::  ind_indoor_target_temp_summer = 12
    INTEGER(iwp) ::  ind_indoor_target_temp_winter = 13
    INTEGER(iwp) ::  ind_lai_r_agfl        = 4    !< index in input list for LAI on roof, above ground floor level
    INTEGER(iwp) ::  ind_lai_r_gfl         = 4  !< index in input list for LAI on roof, ground floor level
    INTEGER(iwp) ::  ind_lai_w_agfl        = 5    !< index in input list for LAI on wall, above ground floor level
    INTEGER(iwp) ::  ind_lai_w_gfl         = 25   !< index in input list for LAI on wall, ground floor level
    INTEGER(iwp) ::  ind_tc1_agfl          = 9    !< index in input list for thermal conductivity at first wall layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_tc1_gfl           = 29   !< index in input list for thermal conductivity at first wall layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc1_wall_r        = 97   !< index in input list for thermal conductivity at first wall layer, roof
    INTEGER(iwp) ::  ind_tc1_win_agfl      = 86   !< index in input list for thermal conductivity at first window layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_tc1_win_gfl       = 74   !< index in input list for thermal conductivity at first window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc1_win_r         = 110  !< index in input list for thermal conductivity at first window layer, roof
    INTEGER(iwp) ::  ind_tc2_agfl          = 10   !< index in input list for thermal conductivity at second wall layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_tc2_gfl           = 30   !< index in input list for thermal conductivity at second wall layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc2_wall_r        = 98   !< index in input list for thermal conductivity at second wall layer, roof
    INTEGER(iwp) ::  ind_tc2_win_agfl      = 87   !< index in input list for thermal conductivity at second window layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_tc2_win_gfl       = 75   !< index in input list for thermal conductivity at second window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc2_win_r         = 111  !< index in input list for thermal conductivity at second window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc3_agfl          = 11   !< index in input list for thermal conductivity at third wall layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_tc3_gfl           = 31   !< index in input list for thermal conductivity at third wall layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc3_wall_r        = 99   !< index in input list for thermal conductivity at third wall layer, roof
    INTEGER(iwp) ::  ind_tc3_win_agfl      = 88   !< index in input list for thermal conductivity at third window layer,
                                                  !< above ground floor level
    INTEGER(iwp) ::  ind_tc3_win_gfl       = 76   !< index in input list for thermal conductivity at third window layer,
                                                  !< ground floor level
    INTEGER(iwp) ::  ind_tc3_win_r         = 112  !< index in input list for thermal conductivity at third window layer, roof
    INTEGER(iwp) ::  ind_thick_1_agfl      = 41   !< index for wall layer thickness - 1st layer above ground floor level
    INTEGER(iwp) ::  ind_thick_1_gfl       = 62   !< index for wall layer thickness - 1st layer ground floor level
    INTEGER(iwp) ::  ind_thick_1_wall_r    = 90   !< index for wall layer thickness - 1st layer roof
    INTEGER(iwp) ::  ind_thick_1_win_agfl  = 79   !< index for window layer thickness - 1st layer above ground floor level
    INTEGER(iwp) ::  ind_thick_1_win_gfl   = 67   !< index for window layer thickness - 1st layer ground floor level
    INTEGER(iwp) ::  ind_thick_1_win_r     = 103  !< index for window layer thickness - 1st layer roof
    INTEGER(iwp) ::  ind_thick_2_agfl      = 42   !< index for wall layer thickness - 2nd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_2_gfl       = 63   !< index for wall layer thickness - 2nd layer ground floor level
    INTEGER(iwp) ::  ind_thick_2_wall_r    = 91   !< index for wall layer thickness - 2nd layer roof
    INTEGER(iwp) ::  ind_thick_2_win_agfl  = 80   !< index for window layer thickness - 2nd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_2_win_gfl   = 68   !< index for window layer thickness - 2nd layer ground floor level
    INTEGER(iwp) ::  ind_thick_2_win_r     = 104  !< index for window layer thickness - 2nd layer roof
    INTEGER(iwp) ::  ind_thick_3_agfl      = 43   !< index for wall layer thickness - 3rd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_3_gfl       = 64   !< index for wall layer thickness - 3rd layer ground floor level
    INTEGER(iwp) ::  ind_thick_3_wall_r    = 92   !< index for wall layer thickness - 3rd layer roof
    INTEGER(iwp) ::  ind_thick_3_win_agfl  = 81   !< index for window layer thickness - 3rd layer above ground floor level
    INTEGER(iwp) ::  ind_thick_3_win_gfl   = 69   !< index for window layer thickness - 3rd layer ground floor level  
    INTEGER(iwp) ::  ind_thick_3_win_r     = 105  !< index for window layer thickness - 3rd layer roof
    INTEGER(iwp) ::  ind_thick_4_agfl      = 44   !< index for wall layer thickness - 4th layer above ground floor level
    INTEGER(iwp) ::  ind_thick_4_gfl       = 65   !< index for wall layer thickness - 4th layer ground floor level
    INTEGER(iwp) ::  ind_thick_4_wall_r    = 93   !< index for wall layer thickness - 4st layer roof
    INTEGER(iwp) ::  ind_thick_4_win_agfl  = 82   !< index for window layer thickness - 4th layer above ground floor level
    INTEGER(iwp) ::  ind_thick_4_win_gfl   = 70   !< index for window layer thickness - 4th layer ground floor level
    INTEGER(iwp) ::  ind_thick_4_win_r     = 106  !< index for window layer thickness - 4th layer roof
    INTEGER(iwp) ::  ind_trans_agfl        = 17   !< index in input list for window transmissivity, above ground floor level
    INTEGER(iwp) ::  ind_trans_gfl         = 35   !< index in input list for window transmissivity, ground floor level
    INTEGER(iwp) ::  ind_trans_r           = 114  !< index in input list for window transmissivity, roof
    INTEGER(iwp) ::  ind_wall_frac_agfl    = 0    !< index in input list for wall fraction, above ground floor level
    INTEGER(iwp) ::  ind_wall_frac_gfl     = 21   !< index in input list for wall fraction, ground floor level
    INTEGER(iwp) ::  ind_wall_frac_r       = 89   !< index in input list for wall fraction, roof
    INTEGER(iwp) ::  ind_win_frac_agfl     = 1    !< index in input list for window fraction, above ground floor level
    INTEGER(iwp) ::  ind_win_frac_gfl      = 22   !< index in input list for window fraction, ground floor level
    INTEGER(iwp) ::  ind_win_frac_r        = 102  !< index in input list for window fraction, roof
    INTEGER(iwp) ::  ind_z0_agfl           = 18   !< index in input list for z0, above ground floor level
    INTEGER(iwp) ::  ind_z0_gfl            = 36   !< index in input list for z0, ground floor level
    INTEGER(iwp) ::  ind_z0qh_agfl         = 19   !< index in input list for z0h / z0q, above ground floor level
    INTEGER(iwp) ::  ind_z0qh_gfl          = 37   !< index in input list for z0h / z0q, ground floor level
    INTEGER(iwp) ::  ind_green_type_roof   = 118  !< index in input list for type of green roof
!
!-- Indices of input attributes in building_surface_pars (except for
!-- radiation-related, which are in radiation_model_mod)
    INTEGER(iwp) ::  ind_s_wall_frac                 = 0  !< index for wall fraction (0-1)
    INTEGER(iwp) ::  ind_s_win_frac                  = 1  !< index for window fraction (0-1)
    INTEGER(iwp) ::  ind_s_green_frac_w              = 2  !< index for green fraction on wall (0-1)
    INTEGER(iwp) ::  ind_s_green_frac_r              = 3  !< index for green fraction on roof (0-1)
    INTEGER(iwp) ::  ind_s_lai_r                     = 4  !< index for leaf area index of green fraction
    INTEGER(iwp) ::  ind_s_hc1                       = 5  !< index for heat capacity of wall layer 1
    INTEGER(iwp) ::  ind_s_hc2                       = 6  !< index for heat capacity of wall layer 2
    INTEGER(iwp) ::  ind_s_hc3                       = 7  !< index for heat capacity of wall layer 3
    INTEGER(iwp) ::  ind_s_tc1                       = 8  !< index for thermal conducivity of wall layer 1
    INTEGER(iwp) ::  ind_s_tc2                       = 9  !< index for thermal conducivity of wall layer 2
    INTEGER(iwp) ::  ind_s_tc3                       = 10 !< index for thermal conducivity of wall layer 3
    INTEGER(iwp) ::  ind_s_indoor_target_temp_summer = 11 !< index for indoor target summer temperature
    INTEGER(iwp) ::  ind_s_indoor_target_temp_winter = 12 !< index for indoor target winter temperature
    INTEGER(iwp) ::  ind_s_emis_wall                 = 13 !< index for emissivity of wall fraction (0-1)
    INTEGER(iwp) ::  ind_s_emis_green                = 14 !< index for emissivity of green fraction (0-1)
    INTEGER(iwp) ::  ind_s_emis_win                  = 15 !< index for emissivity o f window fraction (0-1)
    INTEGER(iwp) ::  ind_s_trans                     = 16 !< index for transmissivity of window fraction (0-1)
    INTEGER(iwp) ::  ind_s_z0                        = 17 !< index for roughness length for momentum (m)
    INTEGER(iwp) ::  ind_s_z0qh                      = 18 !< index for roughness length for heat (m)

    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                         ', &  !< type 0 
                                   'residential - 1950                   ', &  !< type  1 
                                   'residential 1951 - 2000              ', &  !< type  2
                                   'residential 2001 -                   ', &  !< type  3
                                   'office - 1950                        ', &  !< type  4 
                                   'office 1951 - 2000                   ', &  !< type  5
                                   'office 2001 -                        ', &  !< type  6
                                   'bridges                              '  &  !< type  7
                                                                     /)


!
!-- Building facade/wall/green/window properties (partly according to PIDS). 
!-- Initialization of building_pars is outsourced to usm_init_pars. This is
!-- needed because of the huge number of attributes given in building_pars 
!-- (>700), while intel and gfortran compiler have hard limit of continuation 
!-- lines of 511.
    REAL(wp), DIMENSION(0:135,1:7) ::  building_pars
!
!-- 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), TARGET ::  tm_liq_usm_h_m      !< liquid water reservoir tendency (m), horizontal 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 /)
    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_green  = (/0.25_wp,   0.5_wp,    0.75_wp,  1.0_wp /)
                                                                         !< normalized soil, wall and roof, window and
                                                                         !<green layer depths (m/m)

    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
!
!-- 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_3d_data_averaging
       MODULE PROCEDURE usm_3d_data_averaging
    END INTERFACE usm_3d_data_averaging

    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
       MODULE PROCEDURE usm_init
    END INTERFACE usm_init

    INTERFACE usm_init_arrays
       MODULE PROCEDURE usm_init_arrays
    END INTERFACE usm_init_arrays

    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

    
    SAVE

    PRIVATE 

!
!-- Public functions
    PUBLIC usm_boundary_condition, usm_check_parameters, usm_init,               &
           usm_rrd_local,                                                        & 
           usm_surface_energy_balance, usm_material_heat_model,                  &
           usm_swap_timelevel, usm_check_data_output, usm_3d_data_averaging,     &
           usm_data_output_3d, usm_define_netcdf_grid, usm_parin,                &
           usm_wrd_local, usm_init_arrays

!
!-- Public parameters, constants and initial values
    PUBLIC usm_anthropogenic_heat, usm_material_model, usm_wall_mod, &
           usm_green_heat_model, building_pars,                      &
           nzb_wall, 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_init_arrays
    
        IMPLICIT NONE
       
        INTEGER(iwp) ::  l

        IF ( debug_output )  CALL debug_message( 'usm_init_arrays', 'start' )

!
!--Â  Â  Â 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(1:surf_usm_h%ns,0:2) )
        ALLOCATE ( surf_usm_h%albedo(1:surf_usm_h%ns,0:2)      )
        surf_usm_h%albedo_type = albedo_type
        DO  l = 0, 3
           ALLOCATE ( surf_usm_v(l)%albedo_type(1:surf_usm_v(l)%ns,0:2) )
           ALLOCATE ( surf_usm_v(l)%albedo(1:surf_usm_v(l)%ns,0:2)      )
           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
!
!--Â  Â  Â In case the indoor model is applied, allocate memory for waste heat 
!--Â  Â  Â and indoor temperature.
        IF ( indoor_model )  THEN
           ALLOCATE ( surf_usm_h%waste_heat(1:surf_usm_h%ns) )
           surf_usm_h%waste_heat = 0.0_wp
           DO  l = 0, 3
              ALLOCATE ( surf_usm_v(l)%waste_heat(1:surf_usm_v(l)%ns) )
              surf_usm_v(l)%waste_heat = 0.0_wp
           ENDDO
        ENDIF
!
!--Â  Â  Â 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(1:surf_usm_h%ns,0:2) )
         surf_usm_h%frac = 0.0_wp
         DO  l = 0, 3
            ALLOCATE ( surf_usm_v(l)%frac(1:surf_usm_v(l)%ns,0:2) )
            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(1:surf_usm_h%ns,0:2)     )
        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(1:surf_usm_v(l)%ns,0:2)     )
           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%pt_10cm(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)%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) )  
        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

!
!--Â  Â  Â 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)                   )
       tm_liq_usm_h_m%var_usm_1d = 0.0_wp
!
!--Â  Â  Â 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
!
!--Â  Â  Â Initialize building-surface properties, which are also required by other modules, 
!--Â  Â  Â e.g. the indoor model.
        CALL usm_define_pars
        
        IF ( debug_output )  CALL debug_message( 'usm_init_arrays', 'end' )
        
    END SUBROUTINE usm_init_arrays


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Sum up and time-average urban surface output quantities as well as allocate
!> the array necessary for storing the average.
!------------------------------------------------------------------------------!
    SUBROUTINE usm_3d_data_averaging( 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  !< runnin indices
        CHARACTER(LEN=varnamelength)                       :: var                                     !< trimmed variable
        INTEGER(iwp), PARAMETER                            :: nd = 5                                  !< number of directions
        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 /)

        IF ( variable(1:4) == 'usm_' )  THEN  ! is such a check really rquired?

!
!--Â  Â  Â 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_av) )  THEN
                        ALLOCATE ( surf_usm_h%qsws_av(1:surf_usm_h%ns) )
                        surf_usm_h%qsws_av = 0.0_wp
                    ELSE
                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%qsws_av) )  THEN
                           ALLOCATE ( surf_usm_v(l)%qsws_av(1:surf_usm_v(l)%ns) )
                           surf_usm_v(l)%qsws_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_av(m) =                              &
                                          surf_usm_h%qsws_av(m) +           &
                                          surf_usm_h%qsws(m) * l_v
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_av(m) =                        &
                                          surf_usm_v(l)%qsws_av(m) +        &
                                          surf_usm_v(l)%qsws(m) * l_v
                       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
                    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_av(m) =                              &
                                          surf_usm_h%qsws_av(m) /           &
                                          REAL( average_count_3d, kind=wp )
                    ENDDO
                    ELSE
                       DO  m = 1, surf_usm_v(l)%ns
                          surf_usm_v(l)%qsws_av(m) =                        &
                                          surf_usm_v(l)%qsws_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

        ENDIF

    END SUBROUTINE usm_3d_data_averaging



!------------------------------------------------------------------------------!
! 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 = 15      !< 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_t_surf_wall               ', &
                    'usm_t_surf_window             ', &
                    'usm_t_surf_green              ', &
                    'usm_t_green                   ', &
                    'usm_qsws                      ', &
                    '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' )  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

       USE netcdf_data_input_mod,                                             &
            ONLY:  building_type_f

       IMPLICIT NONE

       INTEGER(iwp) ::  i        !< running index, x-dimension
       INTEGER(iwp) ::  j        !< running index, y-dimension

!
!--Â  Â  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( 'usm_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( 'usm_check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
       ENDIF
!
!--Â  Â  Check if building types are set within a valid range.
       IF ( building_type < LBOUND( building_pars, 2 )  .AND.                  &
            building_type > UBOUND( building_pars, 2 ) )  THEN
          WRITE( message_string, * ) 'building_type = ', building_type,        &
                                     ' is out of the valid range'
          CALL message( 'usm_check_parameters', 'PA0529', 2, 2, 0, 6, 0 )
       ENDIF
       IF ( building_type_f%from_file )  THEN
          DO  i = nxl, nxr
             DO  j = nys, nyn
                IF ( building_type_f%var(j,i) /= building_type_f%fill  .AND.   &
              ( building_type_f%var(j,i) < LBOUND( building_pars, 2 )  .OR.    &
                building_type_f%var(j,i) > UBOUND( building_pars, 2 ) ) )      &
                THEN
                   WRITE( message_string, * ) 'building_type = is out of ' //  &
                                              'the valid range at (j,i) = ', j, i
                   CALL message( 'usm_check_parameters', 'PA0529', 2, 2, myid, 6, 0 )
                ENDIF
             ENDDO
          ENDDO
       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        !< flag if averaged
        CHARACTER (len=*), INTENT(IN)  ::  variable  !< variable name
        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(sp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)     ::  temp_pf    !< temp array for urban surface output procedure
        
        CHARACTER (len=varnamelength)                      :: var     !< trimmed variable name
        INTEGER(iwp), PARAMETER                            :: nd = 5  !< number of directions
        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
        INTEGER(iwp)                   :: i,j,k,iwl,istat, l, m  !< running indices

        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 = sp) )
                 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 = sp) + 1.0_sp )
                 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_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(m) * l_v
                    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(m) * l_v
                    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_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_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

                 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
                DO  k = nzb_do, nzt_do
                    local_pf(i,j,k) = temp_pf(k,j,i)
                ENDDO
            ENDDO
        ENDDO
        
    END SUBROUTINE usm_data_output_3d
    

!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Soubroutine defines appropriate grid for netcdf variables.
!> It is called out from subroutine netcdf.
!------------------------------------------------------------------------------!
    SUBROUTINE usm_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
    
        IMPLICIT NONE

        CHARACTER (len=*), INTENT(IN)  ::  variable    !< 
        LOGICAL, INTENT(OUT)           ::  found       !< 
        CHARACTER (len=*), INTENT(OUT) ::  grid_x      !< 
        CHARACTER (len=*), INTENT(OUT) ::  grid_y      !< 
        CHARACTER (len=*), INTENT(OUT) ::  grid_z      !< 

        CHARACTER (len=varnamelength)  :: var

        var = TRIM(variable)
        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:9) == 'usm_qsws_'  .OR.  var(1:13) == 'usm_qsws_veg_'  .OR.              &
             var(1:13) == 'usm_qsws_liq_' .OR.                                              &
             var(1:15) == 'usm_t_surf_wall'  .OR.  var(1:10) == 'usm_t_wall'  .OR.          &
             var(1:17) == 'usm_t_surf_window'  .OR.  var(1:12) == 'usm_t_window'  .OR.      &
             var(1:16) == 'usm_t_surf_green'  .OR. var(1:11) == 'usm_t_green' .OR.          &
             var(1:15) == 'usm_theta_10cm' .OR.                                             &
             var(1:9) == 'usm_surfz'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                &
             var(1:16) == 'usm_surfwintrans'  .OR. var(1:7) == 'usm_swc' ) THEN

            found = .TRUE.
            grid_x = 'x'
            grid_y = 'y'
            grid_z = 'zu'
        ELSE
            found  = .FALSE.
            grid_x = 'none'
            grid_y = 'none'
            grid_z = 'none'
        ENDIF

    END SUBROUTINE usm_define_netcdf_grid
    

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization of the wall surface model
!------------------------------------------------------------------------------!
    SUBROUTINE usm_init_material_model

        IMPLICIT NONE

        INTEGER(iwp) ::  k, l, m            !< running indices
        
        IF ( debug_output )  CALL debug_message( 'usm_init_material_model', 'start' )

!
!--Â  Â  Â Calculate wall grid spacings. 
!--Â  Â  Â Temperature is defined at the center of the wall layers,
!--Â  Â  Â whereas gradients/fluxes are defined at the edges (_stag)Â  Â  Â  
!--Â  Â  Â apply for all particular surface grids. First for horizontal surfaces
        DO  m = 1, surf_usm_h%ns

           surf_usm_h%dz_wall(nzb_wall,m) = surf_usm_h%zw(nzb_wall,m)
           DO k = nzb_wall+1, nzt_wall
               surf_usm_h%dz_wall(k,m) = surf_usm_h%zw(k,m) -                  &
                                         surf_usm_h%zw(k-1,m)
           ENDDO
           surf_usm_h%dz_window(nzb_wall,m) = surf_usm_h%zw_window(nzb_wall,m)
           DO k = nzb_wall+1, nzt_wall
               surf_usm_h%dz_window(k,m) = surf_usm_h%zw_window(k,m) -         &
                                         surf_usm_h%zw_window(k-1,m)
           ENDDO
           
           surf_usm_h%dz_wall(nzt_wall+1,m) = surf_usm_h%dz_wall(nzt_wall,m)

           DO k = nzb_wall, nzt_wall-1
               surf_usm_h%dz_wall_stag(k,m) = 0.5 * (                          &
                           surf_usm_h%dz_wall(k+1,m) + surf_usm_h%dz_wall(k,m) )
           ENDDO
           surf_usm_h%dz_wall_stag(nzt_wall,m) = surf_usm_h%dz_wall(nzt_wall,m)
           
           surf_usm_h%dz_window(nzt_wall+1,m) = surf_usm_h%dz_window(nzt_wall,m)

           DO k = nzb_wall, nzt_wall-1
               surf_usm_h%dz_window_stag(k,m) = 0.5 * (                        &
                           surf_usm_h%dz_window(k+1,m) + surf_usm_h%dz_window(k,m) )
           ENDDO
           surf_usm_h%dz_window_stag(nzt_wall,m) = surf_usm_h%dz_window(nzt_wall,m)

           IF (surf_usm_h%green_type_roof(m) == 2.0_wp ) THEN
!
!-- extensive green roof
!-- set ratio of substrate layer thickness, soil-type and LAI
              soil_type = 3
              surf_usm_h%lai(m) = 2.0_wp
             
              surf_usm_h%zw_green(nzb_wall,m)   = 0.05_wp
              surf_usm_h%zw_green(nzb_wall+1,m) = 0.10_wp
              surf_usm_h%zw_green(nzb_wall+2,m) = 0.15_wp
              surf_usm_h%zw_green(nzb_wall+3,m) = 0.20_wp
           ELSE
!
!-- intensiv green roof
!-- set ratio of substrate layer thickness, soil-type and LAI
              soil_type = 6
              surf_usm_h%lai(m) = 4.0_wp
             
              surf_usm_h%zw_green(nzb_wall,m)   = 0.05_wp
              surf_usm_h%zw_green(nzb_wall+1,m) = 0.10_wp
              surf_usm_h%zw_green(nzb_wall+2,m) = 0.40_wp
              surf_usm_h%zw_green(nzb_wall+3,m) = 0.80_wp
           ENDIF
           
           surf_usm_h%dz_green(nzb_wall,m) = surf_usm_h%zw_green(nzb_wall,m)
           DO k = nzb_wall+1, nzt_wall
               surf_usm_h%dz_green(k,m) = surf_usm_h%zw_green(k,m) -           &
                                         surf_usm_h%zw_green(k-1,m)
           ENDDO
           surf_usm_h%dz_green(nzt_wall+1,m) = surf_usm_h%dz_green(nzt_wall,m)

           DO k = nzb_wall, nzt_wall-1
               surf_usm_h%dz_green_stag(k,m) = 0.5 * (                         &
                           surf_usm_h%dz_green(k+1,m) + surf_usm_h%dz_green(k,m) )
           ENDDO
           surf_usm_h%dz_green_stag(nzt_wall,m) = surf_usm_h%dz_green(nzt_wall,m)
           
          IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
             alpha_vangenuchten = soil_pars(0,soil_type)
          ENDIF

          IF ( l_vangenuchten == 9999999.9_wp )  THEN
             l_vangenuchten = soil_pars(1,soil_type)
          ENDIF

          IF ( n_vangenuchten == 9999999.9_wp )  THEN
             n_vangenuchten = soil_pars(2,soil_type)            
          ENDIF

          IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
             hydraulic_conductivity = soil_pars(3,soil_type)            
          ENDIF

          IF ( saturation_moisture == 9999999.9_wp )  THEN
             saturation_moisture = m_soil_pars(0,soil_type)           
          ENDIF

          IF ( field_capacity == 9999999.9_wp )  THEN
             field_capacity = m_soil_pars(1,soil_type)           
          ENDIF

          IF ( wilting_point == 9999999.9_wp )  THEN
             wilting_point = m_soil_pars(2,soil_type)            
          ENDIF

          IF ( residual_moisture == 9999999.9_wp )  THEN
             residual_moisture = m_soil_pars(3,soil_type)       
          ENDIF
          
          DO k = nzb_wall, nzt_wall+1
             swc_h(k,m) = field_capacity
             rootfr_h(k,m) = 0.5_wp
             surf_usm_h%alpha_vg_green(m)      = alpha_vangenuchten
             surf_usm_h%l_vg_green(m)          = l_vangenuchten
             surf_usm_h%n_vg_green(m)          = n_vangenuchten 
             surf_usm_h%gamma_w_green_sat(k,m) = hydraulic_conductivity
             swc_sat_h(k,m)                    = saturation_moisture
             fc_h(k,m)                         = field_capacity
             wilt_h(k,m)                       = wilting_point
             swc_res_h(k,m)                    = residual_moisture
          ENDDO

        ENDDO

        surf_usm_h%ddz_wall        = 1.0_wp / surf_usm_h%dz_wall
        surf_usm_h%ddz_wall_stag   = 1.0_wp / surf_usm_h%dz_wall_stag
        surf_usm_h%ddz_window      = 1.0_wp / surf_usm_h%dz_window
        surf_usm_h%ddz_window_stag = 1.0_wp / surf_usm_h%dz_window_stag
        surf_usm_h%ddz_green       = 1.0_wp / surf_usm_h%dz_green
        surf_usm_h%ddz_green_stag  = 1.0_wp / surf_usm_h%dz_green_stag
!        
!--Â  Â  Â For vertical surfaces
        DO  l = 0, 3
           DO  m = 1, surf_usm_v(l)%ns
              surf_usm_v(l)%dz_wall(nzb_wall,m) = surf_usm_v(l)%zw(nzb_wall,m)
              DO k = nzb_wall+1, nzt_wall
                  surf_usm_v(l)%dz_wall(k,m) = surf_usm_v(l)%zw(k,m) -         &
                                               surf_usm_v(l)%zw(k-1,m)
              ENDDO
              surf_usm_v(l)%dz_window(nzb_wall,m) = surf_usm_v(l)%zw_window(nzb_wall,m)
              DO k = nzb_wall+1, nzt_wall
                  surf_usm_v(l)%dz_window(k,m) = surf_usm_v(l)%zw_window(k,m) - &
                                               surf_usm_v(l)%zw_window(k-1,m)
              ENDDO
              surf_usm_v(l)%dz_green(nzb_wall,m) = surf_usm_v(l)%zw_green(nzb_wall,m)
              DO k = nzb_wall+1, nzt_wall
                  surf_usm_v(l)%dz_green(k,m) = surf_usm_v(l)%zw_green(k,m) - &
                                               surf_usm_v(l)%zw_green(k-1,m)
              ENDDO
           
              surf_usm_v(l)%dz_wall(nzt_wall+1,m) =                            &
                                              surf_usm_v(l)%dz_wall(nzt_wall,m)

              DO k = nzb_wall, nzt_wall-1
                  surf_usm_v(l)%dz_wall_stag(k,m) = 0.5 * (                    &
                                                surf_usm_v(l)%dz_wall(k+1,m) + &
                                                surf_usm_v(l)%dz_wall(k,m) )
              ENDDO
              surf_usm_v(l)%dz_wall_stag(nzt_wall,m) =                         &
                                              surf_usm_v(l)%dz_wall(nzt_wall,m)
              surf_usm_v(l)%dz_window(nzt_wall+1,m) =                          &
                                              surf_usm_v(l)%dz_window(nzt_wall,m)

              DO k = nzb_wall, nzt_wall-1
                  surf_usm_v(l)%dz_window_stag(k,m) = 0.5 * (                    &
                                                surf_usm_v(l)%dz_window(k+1,m) + &
                                                surf_usm_v(l)%dz_window(k,m) )
              ENDDO
              surf_usm_v(l)%dz_window_stag(nzt_wall,m) =                         &
                                              surf_usm_v(l)%dz_window(nzt_wall,m)
              surf_usm_v(l)%dz_green(nzt_wall+1,m) =                             &
                                              surf_usm_v(l)%dz_green(nzt_wall,m)

              DO k = nzb_wall, nzt_wall-1
                  surf_usm_v(l)%dz_green_stag(k,m) = 0.5 * (                    &
                                                surf_usm_v(l)%dz_green(k+1,m) + &
                                                surf_usm_v(l)%dz_green(k,m) )
              ENDDO
              surf_usm_v(l)%dz_green_stag(nzt_wall,m) =                         &
                                              surf_usm_v(l)%dz_green(nzt_wall,m)
           ENDDO
           surf_usm_v(l)%ddz_wall        = 1.0_wp / surf_usm_v(l)%dz_wall
           surf_usm_v(l)%ddz_wall_stag   = 1.0_wp / surf_usm_v(l)%dz_wall_stag
           surf_usm_v(l)%ddz_window      = 1.0_wp / surf_usm_v(l)%dz_window
           surf_usm_v(l)%ddz_window_stag = 1.0_wp / surf_usm_v(l)%dz_window_stag
           surf_usm_v(l)%ddz_green       = 1.0_wp / surf_usm_v(l)%dz_green
           surf_usm_v(l)%ddz_green_stag  = 1.0_wp / surf_usm_v(l)%dz_green_stag
        ENDDO      

        
        IF ( debug_output )  CALL debug_message( 'usm_init_material_model', 'end' )

    END SUBROUTINE usm_init_material_model


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization of the urban surface model
!------------------------------------------------------------------------------!
    SUBROUTINE usm_init

        USE arrays_3d,                                                         &
            ONLY:  zw

        USE netcdf_data_input_mod,                                             &
            ONLY:  building_pars_f, building_surface_pars_f, building_type_f,  &
                   terrain_height_f
    
        IMPLICIT NONE

        INTEGER(iwp) ::  i                   !< loop index x-dirction
        INTEGER(iwp) ::  ind_alb_green       !< index in input list for green albedo
        INTEGER(iwp) ::  ind_alb_wall        !< index in input list for wall albedo
        INTEGER(iwp) ::  ind_alb_win         !< index in input list for window albedo
        INTEGER(iwp) ::  ind_emis_wall       !< index in input list for wall emissivity
        INTEGER(iwp) ::  ind_emis_green      !< index in input list for green emissivity
        INTEGER(iwp) ::  ind_emis_win        !< index in input list for window emissivity
        INTEGER(iwp) ::  ind_green_frac_w    !< index in input list for green fraction on wall
        INTEGER(iwp) ::  ind_green_frac_r    !< index in input list for green fraction on roof
        INTEGER(iwp) ::  ind_hc1             !< index in input list for heat capacity at first wall layer
        INTEGER(iwp) ::  ind_hc1_win         !< index in input list for heat capacity at first window layer
        INTEGER(iwp) ::  ind_hc2             !< index in input list for heat capacity at second wall layer
        INTEGER(iwp) ::  ind_hc2_win         !< index in input list for heat capacity at second window layer
        INTEGER(iwp) ::  ind_hc3             !< index in input list for heat capacity at third wall layer
        INTEGER(iwp) ::  ind_hc3_win         !< index in input list for heat capacity at third window layer
        INTEGER(iwp) ::  ind_lai_r           !< index in input list for LAI on roof
        INTEGER(iwp) ::  ind_lai_w           !< index in input list for LAI on wall
        INTEGER(iwp) ::  ind_tc1             !< index in input list for thermal conductivity at first wall layer
        INTEGER(iwp) ::  ind_tc1_win         !< index in input list for thermal conductivity at first window layer
        INTEGER(iwp) ::  ind_tc2             !< index in input list for thermal conductivity at second wall layer
        INTEGER(iwp) ::  ind_tc2_win         !< index in input list for thermal conductivity at second window layer
        INTEGER(iwp) ::  ind_tc3             !< index in input list for thermal conductivity at third wall layer
        INTEGER(iwp) ::  ind_tc3_win         !< index in input list for thermal conductivity at third window layer
        INTEGER(iwp) ::  ind_thick_1         !< index in input list for thickness of first wall layer 
        INTEGER(iwp) ::  ind_thick_1_win     !< index in input list for thickness of first window layer 
        INTEGER(iwp) ::  ind_thick_2         !< index in input list for thickness of second wall layer 
        INTEGER(iwp) ::  ind_thick_2_win     !< index in input list for thickness of second window layer
        INTEGER(iwp) ::  ind_thick_3         !< index in input list for thickness of third wall layer
        INTEGER(iwp) ::  ind_thick_3_win     !< index in input list for thickness of third window layer
        INTEGER(iwp) ::  ind_thick_4         !< index in input list for thickness of fourth wall layer
        INTEGER(iwp) ::  ind_thick_4_win     !< index in input list for thickness of fourth window layer
        INTEGER(iwp) ::  ind_trans           !< index in input list for window transmissivity
        INTEGER(iwp) ::  ind_wall_frac       !< index in input list for wall fraction
        INTEGER(iwp) ::  ind_win_frac        !< index in input list for window fraction
        INTEGER(iwp) ::  ind_z0              !< index in input list for z0
        INTEGER(iwp) ::  ind_z0qh            !< index in input list for z0h / z0q
        INTEGER(iwp) ::  is                  !< loop index input surface element
        INTEGER(iwp) ::  j                   !< loop index y-dirction
        INTEGER(iwp) ::  k                   !< loop index z-dirction
        INTEGER(iwp) ::  l                   !< loop index surface orientation
        INTEGER(iwp) ::  m                   !< loop index surface element
        INTEGER(iwp) ::  st                  !< dummy

        LOGICAL      ::  relative_fractions_corrected !< flag indicating if relative surface fractions require normalization

        REAL(wp)     ::  c, tin, twin
        REAL(wp)     ::  ground_floor_level_l         !< local height of ground floor level
        REAL(wp)     ::  sum_frac                     !< sum of the relative material fractions at a surface element
        REAL(wp)     ::  z_agl                        !< height of the surface element above terrain

        IF ( debug_output )  CALL debug_message( 'usm_init', 'start' )

        CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
!
!--Â  Â  Â surface forcing have to be disabled for LSF 
!--Â  Â  Â in case of enabled urban surface module
        IF ( large_scale_forcing )  THEN
            lsf_surf = .FALSE.
        ENDIF
!
!--Â  Â  Â Flag surface elements belonging to the ground floor level. Therefore, 
!--Â  Â  Â use terrain height array from file, if available. This flag is later used
!--Â  Â  Â to control initialization of surface attributes.
!--Â  Â  Â Todo: for the moment disable initialization of building roofs with
!--Â  Â  Â ground-floor-level properties. 
        surf_usm_h%ground_level = .FALSE. 

        DO  l = 0, 3
           surf_usm_v(l)%ground_level = .FALSE.
           DO  m = 1, surf_usm_v(l)%ns
              i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
              j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
              k = surf_usm_v(l)%k(m)
!
!--Â  Â  Â  Â  Â  Â Determine local ground level. Level 1 - default value,
!--Â  Â  Â  Â  Â  Â level 2 - initialization according to building type, 
!--Â  Â  Â  Â  Â  Â level 3 - initialization from value read from file.
              ground_floor_level_l = ground_floor_level
              
              IF ( building_type_f%from_file )  THEN
                  ground_floor_level_l =                                       &
                              building_pars(ind_gflh,building_type_f%var(j,i))
              ENDIF
              
              IF ( building_pars_f%from_file )  THEN
                 IF ( building_pars_f%pars_xy(ind_gflh,j,i) /=                 &
                      building_pars_f%fill )                                   &
                    ground_floor_level_l = building_pars_f%pars_xy(ind_gflh,j,i)
              ENDIF
!
!--Â  Â  Â  Â  Â  Â Determine height of surface element above ground level. Please 
!--Â  Â  Â  Â  Â  Â note, height of surface element is determined with respect to
!--Â  Â  Â  Â  Â  Â its height above ground of the reference grid point in atmosphere,
!--Â  Â  Â  Â  Â  Â Therefore, substract the offset values when assessing the terrain
!--Â  Â  Â  Â  Â  Â height.
              IF ( terrain_height_f%from_file )  THEN
                 z_agl = zw(k) - terrain_height_f%var(j-surf_usm_v(l)%joff,    &
                                                      i-surf_usm_v(l)%ioff)
              ELSE
                 z_agl = zw(k)
              ENDIF
!
!--Â  Â  Â  Â  Â  Â Set flag for ground level
              IF ( z_agl <= ground_floor_level_l )                             &
                 surf_usm_v(l)%ground_level(m) = .TRUE.

           ENDDO
        ENDDO
!
!--Â  Â  Â Initialization of resistances. 
        DO  m = 1, surf_usm_h%ns
           surf_usm_h%r_a(m)        = 50.0_wp
           surf_usm_h%r_a_green(m)  = 50.0_wp
           surf_usm_h%r_a_window(m) = 50.0_wp
        ENDDO
        DO  l = 0, 3
           DO  m = 1, surf_usm_v(l)%ns
              surf_usm_v(l)%r_a(m)        = 50.0_wp
              surf_usm_v(l)%r_a_green(m)  = 50.0_wp
              surf_usm_v(l)%r_a_window(m) = 50.0_wp
           ENDDO
        ENDDO
        
!
!--Â  Â  Map values onto horizontal elemements
       DO  m = 1, surf_usm_h%ns
          surf_usm_h%r_canopy(m)     = 200.0_wp !< canopy_resistance
          surf_usm_h%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
          surf_usm_h%g_d(m)          = 0.0_wp   !< canopy_resistance_coefficient
       ENDDO
!
!--Â  Â  Map values onto vertical elements, even though this does not make
!--Â  Â  much sense.
       DO  l = 0, 3
          DO  m = 1, surf_usm_v(l)%ns
             surf_usm_v(l)%r_canopy(m)     = 200.0_wp !< canopy_resistance
             surf_usm_v(l)%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
             surf_usm_v(l)%g_d(m)          = 0.0_wp   !< canopy_resistance_coefficient
          ENDDO
       ENDDO

!
!--Â  Â  Â Initialize urban-type surface attribute. According to initialization in 
!--Â  Â  Â land-surface model, follow a 3-level approach. 
!--Â  Â  Â Level 1 - initialization via default attributes 
        DO  m = 1, surf_usm_h%ns
!
!--Â  Â  Â  Â  Now, all horizontal surfaces are roof surfaces (?)
           surf_usm_h%isroof_surf(m)   = .TRUE.
           surf_usm_h%surface_types(m) = roof_category         !< default category for root surface
!
!--Â  Â  Â  Â  In order to distinguish between ground floor level and 
!--Â  Â  Â  Â  above-ground-floor level surfaces, set input indices.

           ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
                                     surf_usm_h%ground_level(m) )
           ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
                                     surf_usm_h%ground_level(m) )
           ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
                                     surf_usm_h%ground_level(m) )
           ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
                                     surf_usm_h%ground_level(m) )
!
!--Â  Â  Â  Â  Store building type and its name on each surface element
           surf_usm_h%building_type(m)      = building_type
           surf_usm_h%building_type_name(m) = building_type_name(building_type)
!
!--Â  Â  Â  Â  Initialize relatvie wall- (0), green- (1) and window (2) fractions
           surf_usm_h%frac(m,ind_veg_wall)  = building_pars(ind_wall_frac_r,building_type)   
           surf_usm_h%frac(m,ind_pav_green) = building_pars(ind_green_frac_r,building_type)  
           surf_usm_h%frac(m,ind_wat_win)   = building_pars(ind_win_frac_r,building_type)  
           surf_usm_h%lai(m)                = building_pars(ind_lai_r,building_type)  

           surf_usm_h%rho_c_wall(nzb_wall,m)   = building_pars(ind_hc1_wall_r,building_type)  
           surf_usm_h%rho_c_wall(nzb_wall+1,m) = building_pars(ind_hc1_wall_r,building_type)
           surf_usm_h%rho_c_wall(nzb_wall+2,m) = building_pars(ind_hc2_wall_r,building_type)
           surf_usm_h%rho_c_wall(nzb_wall+3,m) = building_pars(ind_hc3_wall_r,building_type)    
           surf_usm_h%lambda_h(nzb_wall,m)   = building_pars(ind_tc1_wall_r,building_type)  
           surf_usm_h%lambda_h(nzb_wall+1,m) = building_pars(ind_tc1_wall_r,building_type) 
           surf_usm_h%lambda_h(nzb_wall+2,m) = building_pars(ind_tc2_wall_r,building_type)
           surf_usm_h%lambda_h(nzb_wall+3,m) = building_pars(ind_tc3_wall_r,building_type)    
           surf_usm_h%rho_c_green(nzb_wall,m)   = rho_c_soil !building_pars(ind_hc1_wall_r,building_type)  
           surf_usm_h%rho_c_green(nzb_wall+1,m) = rho_c_soil !building_pars(ind_hc1_wall_r,building_type)
           surf_usm_h%rho_c_green(nzb_wall+2,m) = rho_c_soil !building_pars(ind_hc2_wall_r,building_type)
           surf_usm_h%rho_c_green(nzb_wall+3,m) = rho_c_soil !building_pars(ind_hc3_wall_r,building_type)    
           surf_usm_h%lambda_h_green(nzb_wall,m)   = lambda_h_green_sm !building_pars(ind_tc1_wall_r,building_type)  
           surf_usm_h%lambda_h_green(nzb_wall+1,m) = lambda_h_green_sm !building_pars(ind_tc1_wall_r,building_type) 
           surf_usm_h%lambda_h_green(nzb_wall+2,m) = lambda_h_green_sm !building_pars(ind_tc2_wall_r,building_type)
           surf_usm_h%lambda_h_green(nzb_wall+3,m) = lambda_h_green_sm !building_pars(ind_tc3_wall_r,building_type)
           surf_usm_h%rho_c_window(nzb_wall,m)   = building_pars(ind_hc1_win_r,building_type)  
           surf_usm_h%rho_c_window(nzb_wall+1,m) = building_pars(ind_hc1_win_r,building_type)
           surf_usm_h%rho_c_window(nzb_wall+2,m) = building_pars(ind_hc2_win_r,building_type)
           surf_usm_h%rho_c_window(nzb_wall+3,m) = building_pars(ind_hc3_win_r,building_type)    
           surf_usm_h%lambda_h_window(nzb_wall,m)   = building_pars(ind_tc1_win_r,building_type)  
           surf_usm_h%lambda_h_window(nzb_wall+1,m) = building_pars(ind_tc1_win_r,building_type) 
           surf_usm_h%lambda_h_window(nzb_wall+2,m) = building_pars(ind_tc2_win_r,building_type)
           surf_usm_h%lambda_h_window(nzb_wall+3,m) = building_pars(ind_tc3_win_r,building_type)    

           surf_usm_h%target_temp_summer(m)  = building_pars(ind_indoor_target_temp_summer,building_type)    
           surf_usm_h%target_temp_winter(m)  = building_pars(ind_indoor_target_temp_winter,building_type)    
!
!--Â  Â  Â  Â  emissivity of wall-, green- and window fraction 
           surf_usm_h%emissivity(m,ind_veg_wall)  = building_pars(ind_emis_wall_r,building_type)
           surf_usm_h%emissivity(m,ind_pav_green) = building_pars(ind_emis_green_r,building_type)
           surf_usm_h%emissivity(m,ind_wat_win)   = building_pars(ind_emis_win_r,building_type)

           surf_usm_h%transmissivity(m)      = building_pars(ind_trans_r,building_type)

           surf_usm_h%z0(m)                  = building_pars(ind_z0,building_type)
           surf_usm_h%z0h(m)                 = building_pars(ind_z0qh,building_type)
           surf_usm_h%z0q(m)                 = building_pars(ind_z0qh,building_type)
!
!--Â  Â  Â  Â  albedo type for wall fraction, green fraction, window fraction
           surf_usm_h%albedo_type(m,ind_veg_wall)  = INT( building_pars(ind_alb_wall_r,building_type)  )
           surf_usm_h%albedo_type(m,ind_pav_green) = INT( building_pars(ind_alb_green_r,building_type) )
           surf_usm_h%albedo_type(m,ind_wat_win)   = INT( building_pars(ind_alb_win_r,building_type)   )

           surf_usm_h%zw(nzb_wall,m)         = building_pars(ind_thick_1_wall_r,building_type)
           surf_usm_h%zw(nzb_wall+1,m)       = building_pars(ind_thick_2_wall_r,building_type)
           surf_usm_h%zw(nzb_wall+2,m)       = building_pars(ind_thick_3_wall_r,building_type)
           surf_usm_h%zw(nzb_wall+3,m)       = building_pars(ind_thick_4_wall_r,building_type)
           
           surf_usm_h%zw_green(nzb_wall,m)         = building_pars(ind_thick_1_wall_r,building_type)
           surf_usm_h%zw_green(nzb_wall+1,m)       = building_pars(ind_thick_2_wall_r,building_type)
           surf_usm_h%zw_green(nzb_wall+2,m)       = building_pars(ind_thick_3_wall_r,building_type)
           surf_usm_h%zw_green(nzb_wall+3,m)       = building_pars(ind_thick_4_wall_r,building_type)
           
           surf_usm_h%zw_window(nzb_wall,m)         = building_pars(ind_thick_1_win_r,building_type)
           surf_usm_h%zw_window(nzb_wall+1,m)       = building_pars(ind_thick_2_win_r,building_type)
           surf_usm_h%zw_window(nzb_wall+2,m)       = building_pars(ind_thick_3_win_r,building_type)
           surf_usm_h%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4_win_r,building_type)

           surf_usm_h%green_type_roof(m)     = building_pars(ind_green_type_roof,building_type)

        ENDDO

        DO  l = 0, 3
           DO  m = 1, surf_usm_v(l)%ns

              surf_usm_v(l)%surface_types(m) = wall_category         !< default category for root surface
!
!--Â  Â  Â  Â  Â  Â In order to distinguish between ground floor level and 
!--Â  Â  Â  Â  Â  Â above-ground-floor level surfaces, set input indices.
              ind_alb_green    = MERGE( ind_alb_green_gfl,    ind_alb_green_agfl,    &
                                        surf_usm_v(l)%ground_level(m) )
              ind_alb_wall     = MERGE( ind_alb_wall_gfl,     ind_alb_wall_agfl,     &
                                        surf_usm_v(l)%ground_level(m) )
              ind_alb_win      = MERGE( ind_alb_win_gfl,      ind_alb_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_wall_frac    = MERGE( ind_wall_frac_gfl,    ind_wall_frac_agfl,    &
                                        surf_usm_v(l)%ground_level(m) )
              ind_win_frac     = MERGE( ind_win_frac_gfl,     ind_win_frac_agfl,     &
                                        surf_usm_v(l)%ground_level(m) )
              ind_green_frac_w = MERGE( ind_green_frac_w_gfl, ind_green_frac_w_agfl, &
                                        surf_usm_v(l)%ground_level(m) )
              ind_green_frac_r = MERGE( ind_green_frac_r_gfl, ind_green_frac_r_agfl, &
                                        surf_usm_v(l)%ground_level(m) )
              ind_lai_r        = MERGE( ind_lai_r_gfl,        ind_lai_r_agfl,        &
                                        surf_usm_v(l)%ground_level(m) )
              ind_lai_w        = MERGE( ind_lai_w_gfl,        ind_lai_w_agfl,        &
                                        surf_usm_v(l)%ground_level(m) )
              ind_hc1          = MERGE( ind_hc1_gfl,          ind_hc1_agfl,          &
                                        surf_usm_v(l)%ground_level(m) )
              ind_hc1_win      = MERGE( ind_hc1_win_gfl,      ind_hc1_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_hc2          = MERGE( ind_hc2_gfl,          ind_hc2_agfl,          &
                                        surf_usm_v(l)%ground_level(m) )
              ind_hc2_win      = MERGE( ind_hc2_win_gfl,      ind_hc2_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_hc3          = MERGE( ind_hc3_gfl,          ind_hc3_agfl,          &
                                        surf_usm_v(l)%ground_level(m) )
              ind_hc3_win      = MERGE( ind_hc3_win_gfl,      ind_hc3_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_tc1          = MERGE( ind_tc1_gfl,          ind_tc1_agfl,          &
                                        surf_usm_v(l)%ground_level(m) )
              ind_tc1_win      = MERGE( ind_tc1_win_gfl,      ind_tc1_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_tc2          = MERGE( ind_tc2_gfl,          ind_tc2_agfl,          &
                                        surf_usm_v(l)%ground_level(m) )
              ind_tc2_win      = MERGE( ind_tc2_win_gfl,      ind_tc2_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_tc3          = MERGE( ind_tc3_gfl,          ind_tc3_agfl,          &
                                        surf_usm_v(l)%ground_level(m) )
              ind_tc3_win      = MERGE( ind_tc3_win_gfl,      ind_tc3_win_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_1      = MERGE( ind_thick_1_gfl,      ind_thick_1_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_1_win  = MERGE( ind_thick_1_win_gfl,  ind_thick_1_win_agfl,  &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_2      = MERGE( ind_thick_2_gfl,      ind_thick_2_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_2_win  = MERGE( ind_thick_2_win_gfl,  ind_thick_2_win_agfl,  &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_3      = MERGE( ind_thick_3_gfl,      ind_thick_3_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_3_win  = MERGE( ind_thick_3_win_gfl,  ind_thick_3_win_agfl,  &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_4      = MERGE( ind_thick_4_gfl,      ind_thick_4_agfl,      &
                                        surf_usm_v(l)%ground_level(m) )
              ind_thick_4_win  = MERGE( ind_thick_4_win_gfl,  ind_thick_4_win_agfl,  &
                                        surf_usm_v(l)%ground_level(m) )
              ind_emis_wall    = MERGE( ind_emis_wall_gfl,    ind_emis_wall_agfl,    &
                                        surf_usm_v(l)%ground_level(m) )
              ind_emis_green   = MERGE( ind_emis_green_gfl,   ind_emis_green_agfl,   &
                                        surf_usm_v(l)%ground_level(m) )
              ind_emis_win     = MERGE( ind_emis_win_gfl,     ind_emis_win_agfl,     &
                                        surf_usm_v(l)%ground_level(m) )
              ind_trans        = MERGE( ind_trans_gfl,       ind_trans_agfl,         &
                                        surf_usm_v(l)%ground_level(m) )
              ind_z0           = MERGE( ind_z0_gfl,           ind_z0_agfl,           &
                                        surf_usm_v(l)%ground_level(m) )
              ind_z0qh         = MERGE( ind_z0qh_gfl,         ind_z0qh_agfl,         &
                                        surf_usm_v(l)%ground_level(m) )
!
!--Â  Â  Â  Â  Â  Â Store building type and its name on each surface element
              surf_usm_v(l)%building_type(m)      = building_type
              surf_usm_v(l)%building_type_name(m) = building_type_name(building_type)
!
!--Â  Â  Â  Â  Â  Â Initialize relatvie wall- (0), green- (1) and window (2) fractions
            Â