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

Last change on this file since 3115 was 3115, checked in by suehring, 3 years ago

Separate bridges as 3D building objects from normal surface-mounted buildings in terms of correct referencing onto the terrain

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