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

Last change on this file since 4259 was 4259, checked in by suehring, 21 months ago

Normalize relative wall-, green- and window fractions so that they sum-up to one. Give an informative message in case normalization is required.

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