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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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