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

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