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

Last change on this file since 4258 was 4258, checked in by suehring, 21 months ago

Land-surface model: Revise limitation for soil moisture in case it exceeds its saturation value; Revise initialization of soil moisture and temperature in a nested run in case dynamic input information is available. This case, the soil within the child domains can be initialized separately; As part of this revision, migrate the netcdf input of soil temperature / moisture to this module, as well as the routine to inter/extrapolate soil profiles between different grids.; Plant-canopy: Check if any LAD is prescribed when plant-canopy model is applied, in order to avoid crashes in the radiation transfer model; Surface-layer fluxes: Initialization of Obukhov length also at vertical surfaces (if allocated); Urban-surface model: Add checks to ensure that relative fractions of walls, windowns and green surfaces sum-u to one; Revise message calls dealing with local checks

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