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

Last change on this file since 4214 was 4214, checked in by suehring, 22 months ago

Bugfix, missing initialization and clearing of soil-moisture tendency (J.Resler)

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