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

Last change on this file since 3943 was 3943, checked in by maronga, 2 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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