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

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

Bugix: made array declaration in urban_surface_mod Fortran Standard conform

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