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

Last change on this file since 4313 was 4309, checked in by suehring, 18 months ago

Synthetic turbulence generator: Computation of velocity seeds optimized. This implies that random numbers are computed now using the parallel random number generator. Random number are now only computed and normalized locally, while distributed over all mpi ranks afterwards, instead of computing random numbers on a global array. urther, the number of calls for the time-consuming velocity-seed generation is reduced - now the left and right, as well as the north and south boundary share the same velocity-seed matrices.

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