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

Last change on this file since 3337 was 3337, checked in by kanani, 3 years ago

reintegrate branch resler to trunk

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