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

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

Bugfixes and clean-up for output quantity theta_2m*

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