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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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