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

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

Disable also green fraction at vertical building walls; update test configuration for urban_environment test case

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