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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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