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

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

merge fixes of radiation branch (r3362) to trunk

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