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

Last change on this file since 4238 was 4238, checked in by suehring, 20 months ago

Building data base: Indoor-model parameters for some building types adjusted in order to avoid unrealistically high indoor temperatures (S. Rissmann); Indoor model: Bugfix in determination of minimum facade height and in location message, avoid divisions by zero, minor optimizations; radiation model: Modify check in order to avoid equality comparisons of floating points

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