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

Last change on this file since 4077 was 4077, checked in by gronemeier, 2 years ago

Set roughness length z0 and z0h/q at ground-floor level to same value as those above ground-floor level (urban_surface_mod.f90)

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