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

Last change on this file since 4281 was 4267, checked in by suehring, 19 months ago

Indoor model: revision of some parameters and implementation of seasonal-dependent parameters

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