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

Last change on this file since 3248 was 3248, checked in by sward, 3 years ago

Minor format changes

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