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

Last change on this file since 3896 was 3896, checked in by suehring, 2 years ago

Bugfix, wrong index used for accessing building_pars from PIDS

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