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

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

In order to obtain correct surface temperature during spinup set window fraction to zero (only during spinup) instead of just disabling time-integration of window-surface temperature.

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