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

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

Undo accidentally commented initialization

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