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

Last change on this file since 4493 was 4493, checked in by pavelkrc, 5 years ago

Merge brach resler

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