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

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

Move definition of building-surface properties from declaration block to an extra routine; avoid different type kinds

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