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

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

Bugfix: error in placements of CPP statement nopointer fixed

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