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

Last change on this file since 4329 was 4329, checked in by motisi, 2 years ago

Renamed wall_flags_0 to wall_flags_static_0

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