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

Last change on this file since 4230 was 4230, checked in by suehring, 5 years ago

Several bugfixes: profile initialization of chemical species in restart runs; Runge-Kutta tendency array not initialized in chemistry model in restart runs; fixed determination of time indices for chemical emissions (introduced with commit -4227); Update chemistry profiles in offline nesting; initialize canopy resistances for greened building walls (even if green fraction is zero)

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