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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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