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

Last change on this file since 4245 was 4245, checked in by pavelkrc, 21 months ago

Merge branch resler into trunk

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