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

Last change on this file since 4305 was 4305, checked in by suehring, 17 months ago

urban-surface model: bugfix, consider m_liq in restart runs; remove unused variables

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