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

Last change on this file since 4227 was 4227, checked in by gronemeier, 22 months ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

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