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

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

Introduction of addtional surface variables indicating type and name of the surface elements; Bugfix in LSM

  • Property svn:keywords set to Id
File size: 434.0 KB
RevLine 
[2737]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
[2920]18! Copyright 2015-2018 Institute of Computer Science of the
19!                     Czech Academy of Sciences, Prague
[2737]20! Copyright 1997-2018 Leibniz Universitaet Hannover
21!------------------------------------------------------------------------------!
22!
23! Current revisions:
24! ------------------
[2705]25!
[3049]26!
[2737]27! Former revisions:
28! -----------------
29! $Id: urban_surface_mod.f90 3222 2018-08-30 13:35:35Z suehring $
[3222]30! Introduction of surface array for type and its name
31!
32! 3203 2018-08-23 10:48:36Z suehring
[3203]33! Revise bulk parameter for emissivity at ground-floor level
34!
35! 3196 2018-08-13 12:26:14Z maronga
[3196]36! Added maximum aerodynamic resistance of 300 for horiztonal surfaces.
37!
38! 3176 2018-07-26 17:12:48Z suehring
[3176]39! Bugfix, update virtual potential surface temparture, else heat fluxes on
40! roofs might become unphysical
41!
42! 3152 2018-07-19 13:26:52Z suehring
[3152]43! Initialize q_surface, which might be used in surface_layer_fluxes
44!
45! 3151 2018-07-19 08:45:38Z raasch
[3151]46! remaining preprocessor define strings __check removed
47!
48! 3136 2018-07-16 14:48:21Z suehring
[3136]49! Limit also roughness length for heat and moisture where necessary
50!
51! 3123 2018-07-12 16:21:53Z suehring
[3123]52! Correct working precision for INTEGER number
53!
54! 3115 2018-07-10 12:49:26Z suehring
[3115]55! Additional building type to represent bridges
56!
57! 3091 2018-06-28 16:20:35Z suehring
[3091]58! - Limit aerodynamic resistance at vertical walls.
59! - Add check for local roughness length not exceeding surface-layer height and
60!   limit roughness length where necessary.
61!
62! 3065 2018-06-12 07:03:02Z Giersch
[3065]63! Unused array dxdir was removed, dz was replaced by dzu to consider vertical
64! grid stretching
65!
66! 3049 2018-05-29 13:52:36Z Giersch
[3049]67! Error messages revised
68!
69! 3045 2018-05-28 07:55:41Z Giersch
[3045]70! Error message added
71!
72! 3029 2018-05-23 12:19:17Z raasch
[3029]73! bugfix: close unit 151 instead of 90
74!
75! 3014 2018-05-09 08:42:38Z maronga
[3014]76! Added pc_transpiration_rate
77!
78! 2977 2018-04-17 10:27:57Z kanani
[2977]79! Implement changes from branch radiation (r2948-2971) with minor modifications.
80! (moh.hefny):
81! Extended exn for all model domain height to avoid the need to get nzut.
82!
83! 2963 2018-04-12 14:47:44Z suehring
[2963]84! Introduce index for vegetation/wall, pavement/green-wall and water/window
85! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
86!
87! 2943 2018-04-03 16:17:10Z suehring
[2943]88! Calculate exner function at all height levels and remove some un-used
89! variables.
90!
91! 2932 2018-03-26 09:39:22Z maronga
[2932]92! renamed urban_surface_par to urban_surface_parameters
93!
94! 2921 2018-03-22 15:05:23Z Giersch
[2921]95! The activation of spinup has been moved to parin
96!
97! 2920 2018-03-22 11:22:01Z kanani
[2920]98! Remove unused pcbl, npcbl from ONLY list
99! moh.hefny:
100! Fixed bugs introduced by new structures and by moving radiation interaction
101! into radiation_model_mod.f90.
102! Bugfix: usm data output 3D didn't respect directions
103!
104! 2906 2018-03-19 08:56:40Z Giersch
[2906]105! Local variable ids has to be initialized with a value of -1 in
106! usm_average_3d_data
107!
108! 2894 2018-03-15 09:17:58Z Giersch
[2894]109! Calculations of the index range of the subdomain on file which overlaps with
110! the current subdomain are already done in read_restart_data_mod,
111! usm_read/write_restart_data have been renamed to usm_r/wrd_local, variable
112! named found has been introduced for checking if restart data was found,
113! reading of restart strings has been moved completely to
114! read_restart_data_mod, usm_rrd_local is already inside the overlap loop
115! programmed in read_restart_data_mod, SAVE attribute added where necessary,
116! deallocation and allocation of some arrays have been changed to take care of
117! different restart files that can be opened (index i), the marker *** end usm
118! *** is not necessary anymore, strings and their respective lengths are
119! written out and read now in case of restart runs to get rid of prescribed
120! character lengths
121!
122! 2805 2018-02-14 17:00:09Z suehring
[2805]123! Initialization of resistances.
124!
125! 2797 2018-02-08 13:24:35Z suehring
[2797]126! Comment concerning output of ground-heat flux added.
127!
128! 2766 2018-01-22 17:17:47Z kanani
[2766]129! Removed redundant commas, added some blanks
130!
131! 2765 2018-01-22 11:34:58Z maronga
[2765]132! Major bugfix in calculation of f_shf. Adjustment of roughness lengths in
133! building_pars
134!
135! 2750 2018-01-15 16:26:51Z knoop
[2746]136! Move flag plant canopy to modules
137!
138! 2737 2018-01-11 14:58:11Z kanani
[2737]139! Removed unused variables t_surf_whole...
140!
141! 2735 2018-01-11 12:01:27Z suehring
142! resistances are saved in surface attributes
143!
[2735]144! 2723 2018-01-05 09:27:03Z maronga
[2737]145! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
146!
147! 2720 2018-01-02 16:27:15Z kanani
148! Correction of comment
149!
150! 2718 2018-01-02 08:49:38Z maronga
151! Corrected "Former revisions" section
152!
153! 2705 2017-12-18 11:26:23Z maronga
154! Changes from last commit documented
155!
156! 2703 2017-12-15 20:12:38Z maronga
157! Workaround for calculation of r_a
158!
159! 2696 2017-12-14 17:12:51Z kanani
160! - Change in file header (GPL part)
161! - Bugfix in calculation of pt_surface and related fluxes. (BM)
162! - Do not write surface temperatures onto pt array as this might cause
163!   problems with nesting. (MS)
164! - Revised calculation of pt1 (now done in surface_layer_fluxes).
165!   Bugfix, f_shf_window and f_shf_green were not set at vertical surface
166!   elements. (MS)
167! - merged with branch ebsolver
168!   green building surfaces do not evaporate yet
169!   properties of green wall layers and window layers are taken from wall layers
170!   this input data is missing. (RvT)
171! - Merged with branch radiation (developed by Mohamed Salim)
172! - Revised initialization. (MS)
173! - Rename emiss_surf into emissivity, roughness_wall into z0, albedo_surf into
174!   albedo. (MS)
175! - Move first call of usm_radiatin from usm_init to init_3d_model
176! - fixed problem with near surface temperature
177! - added near surface temperature t_surf_10cm_h(m), t_surf_10cm_v(l)%t(m)
178! - does not work with temp profile including stability, ol
179!   t_surf_10cm = pt1 now
180! - merged with 2357 bugfix, error message for nopointer version
181! - added indoor model coupling with wall heat flux
182! - added green substrate/ dry vegetation layer for buildings
183! - merged with 2232 new surface-type structure
184! - added transmissivity of window tiles
185! - added MOSAIK tile approach for 3 different surfaces (RvT)
186!
187! 2583 2017-10-26 13:58:38Z knoop
188! Bugfix: reverted MPI_Win_allocate_cptr introduction in last commit
189!
190! 2582 2017-10-26 13:19:46Z hellstea
191! Workaround for gnufortran compiler added in usm_calc_svf. CALL MPI_Win_allocate is
192! replaced by CALL MPI_Win_allocate_cptr if defined ( __gnufortran ).
193!
194! 2544 2017-10-13 18:09:32Z maronga
195! Date and time quantities are now read from date_and_time_mod. Solar constant is
196! read from radiation_model_mod
197!
198! 2516 2017-10-04 11:03:04Z suehring
199! Remove tabs
200!
201! 2514 2017-10-04 09:52:37Z suehring
202! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
203! no output of ghost layer data
204!
205! 2350 2017-08-15 11:48:26Z kanani
206! Bugfix and error message for nopointer version.
207! Additional "! defined(__nopointer)" as workaround to enable compilation of
208! nopointer version.
209!
210! 2318 2017-07-20 17:27:44Z suehring
211! Get topography top index via Function call
212!
213! 2317 2017-07-20 17:27:19Z suehring
214! Bugfix: adjust output of shf. Added support for spinups
215!
216! 2287 2017-06-15 16:46:30Z suehring
217! Bugfix in determination topography-top index
218!
219! 2269 2017-06-09 11:57:32Z suehring
220! Enable restart runs with different number of PEs
221! Bugfixes nopointer branch
222!
223! 2258 2017-06-08 07:55:13Z suehring
224! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
225!
226! 2233 2017-05-30 18:08:54Z suehring
227!
228! 2232 2017-05-30 17:47:52Z suehring
229! Adjustments according to new surface-type structure. Remove usm_wall_heat_flux;
230! insteat, heat fluxes are directly applied in diffusion_s.
231!
232! 2213 2017-04-24 15:10:35Z kanani
233! Removal of output quantities usm_lad and usm_canopy_hr
234!
235! 2209 2017-04-19 09:34:46Z kanani
236! cpp switch __mpi3 removed,
237! minor formatting,
238! small bugfix for division by zero (Krc)
239!
240! 2113 2017-01-12 13:40:46Z kanani
241! cpp switch __mpi3 added for MPI-3 standard code (Ketelsen)
242!
243! 2071 2016-11-17 11:22:14Z maronga
244! Small bugfix (Resler)
245!
246! 2031 2016-10-21 15:11:58Z knoop
247! renamed variable rho to rho_ocean
248!
249! 2024 2016-10-12 16:42:37Z kanani
250! Bugfixes in deallocation of array plantt and reading of csf/csfsurf,
251! optimization of MPI-RMA operations,
252! declaration of pcbl as integer,
253! renamed usm_radnet -> usm_rad_net, usm_canopy_khf -> usm_canopy_hr,
254! splitted arrays svf -> svf & csf, svfsurf -> svfsurf & csfsurf,
255! use of new control parameter varnamelength,
256! added output variables usm_rad_ressw, usm_rad_reslw,
257! minor formatting changes,
258! minor optimizations.
259!
260! 2011 2016-09-19 17:29:57Z kanani
261! Major reformatting according to PALM coding standard (comments, blanks,
262! alphabetical ordering, etc.),
263! removed debug_prints,
264! removed auxiliary SUBROUTINE get_usm_info, instead, USM flag urban_surface is
265! defined in MODULE control_parameters (modules.f90) to avoid circular
266! dependencies,
267! renamed canopy_heat_flux to pc_heating_rate, as meaning of quantity changed.
268!
269! 2007 2016-08-24 15:47:17Z kanani
270! Initial revision
271!
272!
273! Description:
274! ------------
275! 2016/6/9 - Initial version of the USM (Urban Surface Model)
276!            authors: Jaroslav Resler, Pavel Krc
277!                     (Czech Technical University in Prague and Institute of
278!                      Computer Science of the Czech Academy of Sciences, Prague)
279!            with contributions: Michal Belda, Nina Benesova, Ondrej Vlcek
280!            partly inspired by PALM LSM (B. Maronga)
281!            parameterizations of Ra checked with TUF3D (E. S. Krayenhoff)
282!> Module for Urban Surface Model (USM)
283!> The module includes:
284!>    1. radiation model with direct/diffuse radiation, shading, reflections
285!>       and integration with plant canopy
286!>    2. wall and wall surface model
287!>    3. surface layer energy balance
288!>    4. anthropogenic heat (only from transportation so far)
289!>    5. necessary auxiliary subroutines (reading inputs, writing outputs,
290!>       restart simulations, ...)
291!> It also make use of standard radiation and integrates it into
292!> urban surface model.
293!>
294!> Further work:
295!> -------------
[2920]296!> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity
[2737]297!>    from surfaces that are visible from local surfaces (i.e. there is a SVF
298!>    where target is local). To do that, radiosity will be exchanged after each
299!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
300!>
[2920]301!> 2. Temporarily large values of surface heat flux can be observed, up to
[2737]302!>    1.2 Km/s, which seem to be not realistic.
303!>
304!> @todo Output of _av variables in case of restarts
305!> @todo Revise flux conversion in energy-balance solver
306!> @todo Bugfixing in nopointer branch
307!> @todo Check optimizations for RMA operations
308!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
309!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
310!>       factor 3 between min and max time
311!> @todo Move setting of flag indoor_model to indoor_model_mod once available
312!> @todo Check divisions in wtend (etc.) calculations for possible division
313!>       by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
314!> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
[3014]315!> @todo Move plant canopy stuff into plant canopy code
[2737]316!------------------------------------------------------------------------------!
317 MODULE urban_surface_mod
318
319#if ! defined( __nopointer )
320    USE arrays_3d,                                                             &
[3065]321        ONLY:  dzu, hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
[2737]322#endif
323
324    USE cloud_parameters,                                                      &
325        ONLY:  cp, r_d
326
327    USE constants,                                                             &
328        ONLY:  pi
329   
330    USE control_parameters,                                                    &
[3152]331        ONLY:  coupling_start_time, topography, dt_3d, humidity,               &
[2737]332               intermediate_timestep_count, initializing_actions,              &
333               intermediate_timestep_count_max, simulated_time, end_time,      &
334               timestep_scheme, tsc, coupling_char, io_blocks, io_group,       &
335               message_string, time_since_reference_point, surface_pressure,   &
336               g, pt_surface, large_scale_forcing, lsf_surf, spinup,           &
337               spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
[2746]338               average_count_3d, varnamelength, urban_surface, kappa,          &
339               plant_canopy
[2737]340
341    USE cpulog,                                                                &
342        ONLY:  cpu_log, log_point, log_point_s
343
344    USE date_and_time_mod,                                                     &
[2920]345        ONLY:  time_utc_init
[2737]346
347    USE grid_variables,                                                        &
348        ONLY:  dx, dy, ddx, ddy, ddx2, ddy2
349
350    USE indices,                                                               &
351        ONLY:  nx, ny, nnx, nny, nnz, nxl, nxlg, nxr, nxrg, nyn, nyng, nys,    &
352               nysg, nzb, nzt, nbgp, wall_flags_0
353
354    USE, INTRINSIC :: iso_c_binding
355
356    USE kinds
357             
358    USE pegrid
359   
360    USE plant_canopy_model_mod,                                                &
[3014]361        ONLY:  pc_heating_rate, pc_transpiration_rate
[2737]362   
363    USE radiation_model_mod,                                                   &
[2943]364        ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,        &
[3045]365               radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
[2737]366               sigma_sb, solar_constant, sun_direction, sun_dir_lat,           &
367               sun_dir_lon,                                                    &
368               force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
369               surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
[2920]370               surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
[2737]371               rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
[2920]372               surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,             &
373               iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,             &
374               inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
[2943]375               iz, iy, ix, idir, jdir, kdir,  nsurf_type, nsurf, idsvf, ndsvf, &
[2920]376               iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,           &
377               idcsf, ndcsf, kdcsf, pct,                                       &
378               startland, endland, startwall, endwall, skyvf, skyvft
[2737]379
380    USE statistics,                                                            &
381        ONLY:  hom, statistic_regions
382
[2963]383    USE surface_mod,                                                           &
384        ONLY:  get_topography_top_index_ji, get_topography_top_index,          &
385               ind_pav_green, ind_veg_wall, ind_wat_win, surf_usm_h,           &
386               surf_usm_v, surface_restore_elements
[2737]387
388
389    IMPLICIT NONE
390
391
392!-- configuration parameters (they can be setup in PALM config)
[3151]393    LOGICAL ::  usm_material_model = .TRUE.        !< flag parameter indicating wheather the  model of heat in materials is used
394    LOGICAL ::  usm_anthropogenic_heat = .FALSE.   !< flag parameter indicating wheather the anthropogenic heat sources (e.g.transportation) are used
395    LOGICAL ::  force_radiation_call_l = .FALSE.   !< flag parameter for unscheduled radiation model calls
396    LOGICAL ::  indoor_model = .FALSE.             !< whether to use the indoor model
397    LOGICAL ::  read_wall_temp_3d = .FALSE.
[2737]398
[2920]399
[3151]400    INTEGER(iwp) ::  building_type = 1               !< default building type (preleminary setting)
401    INTEGER(iwp) ::  land_category = 2               !< default category for land surface
402    INTEGER(iwp) ::  wall_category = 2               !< default category for wall surface over pedestrian zone
403    INTEGER(iwp) ::  pedestrian_category = 2         !< default category for wall surface in pedestrian zone
404    INTEGER(iwp) ::  roof_category = 2               !< default category for root surface
405    REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
[2737]406!
407!-- Indices of input attributes for (above) ground floor level
408    INTEGER(iwp) ::  ind_alb_wall          = 38 !< index in input list for albedo_type of wall fraction
409    INTEGER(iwp) ::  ind_alb_green         = 39 !< index in input list for albedo_type of green fraction
410    INTEGER(iwp) ::  ind_alb_win           = 40 !< index in input list for albedo_type of window fraction
411    INTEGER(iwp) ::  ind_emis_wall_agfl    = 14 !< index in input list for wall emissivity, above ground floor level
412    INTEGER(iwp) ::  ind_emis_wall_gfl     = 32 !< index in input list for wall emissivity, ground floor level
413    INTEGER(iwp) ::  ind_emis_green_agfl   = 15 !< index in input list for green emissivity, above ground floor level
414    INTEGER(iwp) ::  ind_emis_green_gfl    = 33 !< index in input list for green emissivity, ground floor level
415    INTEGER(iwp) ::  ind_emis_win_agfl     = 16 !< index in input list for window emissivity, above ground floor level
416    INTEGER(iwp) ::  ind_emis_win_gfl      = 34 !< index in input list for window emissivity, ground floor level
417    INTEGER(iwp) ::  ind_green_frac_w_agfl = 2  !< index in input list for green fraction on wall, above ground floor level
418    INTEGER(iwp) ::  ind_green_frac_w_gfl  = 23 !< index in input list for green fraction on wall, ground floor level
419    INTEGER(iwp) ::  ind_green_frac_r_agfl = 3  !< index in input list for green fraction on roof, above ground floor level
420    INTEGER(iwp) ::  ind_green_frac_r_gfl  = 24 !< index in input list for green fraction on roof, ground floor level
421    INTEGER(iwp) ::  ind_hc1_agfl          =  6 !< index in input list for heat capacity at first wall layer, above ground floor level
422    INTEGER(iwp) ::  ind_hc1_gfl           = 26 !< index in input list for heat capacity at first wall layer, ground floor level
423    INTEGER(iwp) ::  ind_hc2_agfl          = 7  !< index in input list for heat capacity at second wall layer, above ground floor level
424    INTEGER(iwp) ::  ind_hc2_gfl           = 27 !< index in input list for heat capacity at second wall layer, ground floor level
425    INTEGER(iwp) ::  ind_hc3_agfl          = 8  !< index in input list for heat capacity at third wall layer, above ground floor level
426    INTEGER(iwp) ::  ind_hc3_gfl           = 28 !< index in input list for heat capacity at third wall layer, ground floor level
427    INTEGER(iwp) ::  ind_gflh              = 20 !< index in input list for ground floor level height
428    INTEGER(iwp) ::  ind_lai_r_agfl        = 4  !< index in input list for LAI on roof, above ground floor level
429    INTEGER(iwp) ::  ind_lai_r_gfl         = 4  !< index in input list for LAI on roof, ground floor level
430    INTEGER(iwp) ::  ind_lai_w_agfl        = 5  !< index in input list for LAI on wall, above ground floor level
431    INTEGER(iwp) ::  ind_lai_w_gfl         = 25 !< index in input list for LAI on wall, ground floor level
432    INTEGER(iwp) ::  ind_tc1_agfl          = 9  !< index in input list for thermal conductivity at first wall layer, above ground floor level
433    INTEGER(iwp) ::  ind_tc1_gfl           = 29 !< index in input list for thermal conductivity at first wall layer, ground floor level
434    INTEGER(iwp) ::  ind_tc2_agfl          = 10 !< index in input list for thermal conductivity at second wall layer, above ground floor level
435    INTEGER(iwp) ::  ind_tc2_gfl           = 30 !< index in input list for thermal conductivity at second wall layer, ground floor level
436    INTEGER(iwp) ::  ind_tc3_agfl          = 11 !< index in input list for thermal conductivity at third wall layer, above ground floor level
437    INTEGER(iwp) ::  ind_tc3_gfl           = 31 !< index in input list for thermal conductivity at third wall layer, ground floor level
438    INTEGER(iwp) ::  ind_thick_1           = 41 !< index for wall layer thickness - 1st layer
439    INTEGER(iwp) ::  ind_thick_2           = 42 !< index for wall layer thickness - 2nd layer
440    INTEGER(iwp) ::  ind_thick_3           = 43 !< index for wall layer thickness - 3rd layer
441    INTEGER(iwp) ::  ind_thick_4           = 44 !< index for wall layer thickness - 4th layer
442    INTEGER(iwp) ::  ind_trans_agfl        = 17 !< index in input list for window transmissivity, above ground floor level
443    INTEGER(iwp) ::  ind_trans_gfl         = 35 !< index in input list for window transmissivity, ground floor level
444    INTEGER(iwp) ::  ind_wall_frac_agfl    = 0  !< index in input list for wall fraction, above ground floor level
445    INTEGER(iwp) ::  ind_wall_frac_gfl     = 21 !< index in input list for wall fraction, ground floor level
446    INTEGER(iwp) ::  ind_win_frac_agfl     = 1  !< index in input list for window fraction, above ground floor level
447    INTEGER(iwp) ::  ind_win_frac_gfl      = 22 !< index in input list for window fraction, ground floor level
448    INTEGER(iwp) ::  ind_z0_agfl           = 18 !< index in input list for z0, above ground floor level
449    INTEGER(iwp) ::  ind_z0_gfl            = 36 !< index in input list for z0, ground floor level
450    INTEGER(iwp) ::  ind_z0qh_agfl         = 19 !< index in input list for z0h / z0q, above ground floor level
451    INTEGER(iwp) ::  ind_z0qh_gfl          = 37 !< index in input list for z0h / z0q, ground floor level
452
453
454    REAL(wp)  ::  roof_height_limit = 4._wp          !< height for distinguish between land surfaces and roofs
455    REAL(wp)  ::  ground_floor_level = 4.0_wp        !< default ground floor level
456
[2920]457
[3115]458    CHARACTER(37), DIMENSION(0:7), PARAMETER :: building_type_name = (/     &
[2737]459                                   'user-defined                         ', & !  0
460                                   'residential - 1950                   ', & !  1
461                                   'residential 1951 - 2000              ', & !  2
462                                   'residential 2001 -                   ', & !  3
463                                   'office - 1950                        ', & !  4
464                                   'office 1951 - 2000                   ', & !  5
[3115]465                                   'office 2001 -                        ', & !  6
466                                   'bridges                              '  & !  7
[2737]467                                                                     /)
468!
469!-- building parameters, 4 different types
470!-- 0 - wall fraction, 1- window fraction, 2 - green fraction on wall, 3- green fraction
471!-- at roof, 4 - lai of green fraction at roof,  5 - lai of green fraction at wall,
472!-- 6 - heat capacity of wall layer 1, 7 - heat capacity of wall layer 2,
473!-- 8 - heat capacity of wall layer 3, 9 - thermal conductivity of wall layer 1,
474!-- 10 - thermal conductivity of wall layer 2, 11 - thermal conductivity of wall layer 3, 
475!-- 12 - indoor target summer temperature ( K ), 13 - indoor target winter temperature (K),
476!-- 14 - emissivity of wall fraction, 15 - emissivity of green fraction, 16 - emissivity of window fraction,
477!-- 17 - transmissivity of window fraction, 18 - z0, 19 - z0h/z0q, 20 - ground floor height,
478!-- 21 - ground floor wall fraction, 22 - ground floor window fraction, 23 ground floor green fraction,
479!-- 24 - ground floor green fraction on roof, 25 - ground floor lai of green fraction,
480!-- 26 - ground floor heat capacity of wall layer 1, 27 - ground floor heat capacity of wall layer 1,
481!-- 28 - ground floor heat capacity of wall layer 3, 29 - ground floor thermal conductivity of wall layer 1,
482!-- 30 - ground floor thermal conductivity of wall layer 2, 31 - ground floor thermal conductivity of wall layer 3,
483!-- 32 - ground floor emissivity of wall fraction, 33 - ground floor emissivity of green fraction,
484!-- 34 - ground floor emissivity of window fraction, 35 - ground floor transmissivity of window fraction,
485!-- 36 - ground floor z0, 37 - ground floor z0h/z0q, 38 - albedo type wall fraction
486!-- 39 - albedo type green fraction, 40 - albedo type window fraction
487!-- 41 - wall layer thickness - 1st layer, 42 - wall layer thickness - 2nd layer,
488!-- 43 - wall layer thickness - 3rd layer, 44 - wall layer thickness - 4th layer,
489!-- 45 - heat capacity of the wall surface, 46 - heat conductivity
490!-- Please note, only preleminary dummy values so far!
[3115]491    REAL(wp), DIMENSION(0:46,1:7), PARAMETER :: building_pars = RESHAPE( (/    &
[2737]492        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
493        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
494        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
[3115]495        0.001_wp, 0.0001_wp, 4.0_wp,                                           & !parameter 18-20
[2737]496        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
497        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
498        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]499        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[2737]500        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
501        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
502        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 1
503        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
504        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
505        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
[3115]506        0.001_wp, 0.0001_wp, 4.0_wp,                                           & !parameter 18-20
[2737]507        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
508        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
509        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]510        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[2737]511        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
512        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
513        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 2
514        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
515        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
516        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
[3115]517        0.001_wp, 0.0001_wp, 4.0_wp,                                           & !parameter 18-20
[2737]518        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
519        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
520        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]521        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[2737]522        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
523        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
524        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 3
525        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
526        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
527        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
528        0.01_wp, 0.001_wp, 4.0_wp,                                             & !parameter 18-20
529        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
530        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
531        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]532        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[2737]533        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
534        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
535        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 4
536        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
537        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
538        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
[3115]539        0.001_wp, 0.0001_wp, 4.0_wp,                                           & !parameter 18-20
[2737]540        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
541        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
542        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]543        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[2737]544        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
545        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
546        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 5
547        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
548        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
549        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
[3115]550        0.001_wp, 0.0001_wp, 4.0_wp,                                           & !parameter 18-20
[2737]551        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
552        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
553        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]554        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[2737]555        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
556        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
[3115]557        20000.0_wp, 10.0_wp,                                                   & !parameter 45-46 - end of type 6
558        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 1.0_wp,                        & !parameter 0-5
559        1000000.0_wp, 1000000.0_wp, 1000000.0_wp, 0.3_wp, 0.3_wp, 0.3_wp,      & !parameter 6-11
560        296.15_wp, 293.15_wp, 0.9_wp, 0.9_wp, 0.01_wp, 0.99_wp,                & !parameter 12-17
561        0.001_wp, 0.0001_wp, 0.0_wp,                                           & !parameter 18-20
562        1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 3.0_wp,                                & !parameter 21-25
563        1000000.0_wp, 1000000.0_wp, 1000000.0_wp,                              & !parameter 26-28                     
564        0.3_wp, 0.3_wp, 0.3_wp,                                                & !parameter 29-31       
[3203]565        0.9_wp, 0.4_wp, 0.4_wp, 0.4_wp, 0.01_wp, 0.001_wp,                     & !parameter 32-37
[3115]566        24.0_wp, 24.0_wp, 24.0_wp,                                             & !parameter 38-40
567        0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp,                                & !parameter 41-44
568        20000.0_wp, 10.0_wp                                                    & !parameter 45-46 - end of type 7 (bridges)
[2737]569                                                                          /),  &
[3115]570                                                               (/47, 7/) )
[2737]571
572!
573!-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls.
574    TYPE t_surf_vertical
575       REAL(wp), DIMENSION(:), ALLOCATABLE         :: t
576    END TYPE t_surf_vertical
577!
578!-- Type for wall temperatures at vertical walls. Is not necessary for horizontal walls.
579    TYPE t_wall_vertical
580       REAL(wp), DIMENSION(:,:), ALLOCATABLE       :: t
581    END TYPE t_wall_vertical
582
583
584!-- arrays for time averages
585!-- 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
586    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rad_net_av       !< average of rad_net_l
587    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw_av      !< average of sw radiation falling to local surface including radiation from reflections
588    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw_av      !< average of lw radiation falling to local surface including radiation from reflections
589    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdir_av   !< average of direct sw radiation falling to local surface
590    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswdif_av   !< average of diffuse sw radiation from sky and model boundary falling to local surface
591    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwdif_av   !< average of diffuse lw radiation from sky and model boundary falling to local surface
592    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinswref_av   !< average of sw radiation falling to surface from reflections
593    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlwref_av   !< average of lw radiation falling to surface from reflections
594    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw_av     !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
595    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw_av     !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
596    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
597    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
598    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf_av        !< average of total radiation flux incoming to minus outgoing from local surface 
599    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_av       !< average of wghf_eb
600    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wshf_eb_av       !< average of wshf_eb
601    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  t_wall_av        !< Average of t_wall
602    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_green_av !< average of wghf_eb_green
603    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  t_green_av       !< Average of t_green
604    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  wghf_eb_window_av !< average of wghf_eb_window
605    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  t_window_av      !< Average of t_window   
606   
607
608!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
609!-- anthropogenic heat sources
610!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2920]611    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE        ::  aheat             !< daily average of anthropogenic heat (W/m2)
612    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheatprof         !< diurnal profiles of anthropogenic heat for particular layers
[3123]613    INTEGER(iwp)                                   ::  naheatlayers = 1  !< number of layers of anthropogenic heat
[2737]614
615!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
616!-- wall surface model
617!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
618!-- wall surface model constants
619    INTEGER(iwp), PARAMETER                        :: nzb_wall = 0       !< inner side of the wall model (to be switched)
620    INTEGER(iwp), PARAMETER                        :: nzt_wall = 3       !< outer side of the wall model (to be switched)
621    INTEGER(iwp), PARAMETER                        :: nzw = 4            !< number of wall layers (fixed for now)
622
623    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
624                                                                         !< normalized soil, wall and roof layer depths (m/m)
625!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default = (/0.33_wp, 0.66_wp, 1.0_wp /)
626    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
627!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.33_wp, 0.66_wp, 1.0_wp /)
628!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_window = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
629                                                                         !< normalized window layer depths (m/m)
630!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.0242_wp, 0.0969_wp, 0.346_wp, 1.0_wp /)
631                                                                         !< normalized green layer depths (m/m)
632    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.25_wp, 0.5_wp, 0.75_wp, 1.0_wp /)
633!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /)
634
635
[2920]636    REAL(wp)                                       :: wall_inner_temperature = 295.0_wp    !< temperature of the inner wall surface (~22 degrees C) (K)
637    REAL(wp)                                       :: roof_inner_temperature = 295.0_wp    !< temperature of the inner roof surface (~22 degrees C) (K)
638    REAL(wp)                                       :: soil_inner_temperature = 288.0_wp    !< temperature of the deep soil (~15 degrees C) (K)
639    REAL(wp)                                       :: window_inner_temperature = 295.0_wp  !< temperature of the inner window surface (~22 degrees C) (K)
640
[2737]641!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
642!-- surface and material model variables for walls, ground, roofs
643!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
644    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn                !< normalized wall layer depths (m)
645    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn_window         !< normalized window layer depths (m)
646    REAL(wp), DIMENSION(:), ALLOCATABLE            :: zwn_green          !< normalized green layer depths (m)
647
648#if defined( __nopointer )
649    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h           !< wall surface temperature (K) at horizontal walls
650    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_p         !< progn. wall surface temperature (K) at horizontal walls
651    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h    !< window surface temperature (K) at horizontal walls
652    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_p  !< progn. window surface temperature (K) at horizontal walls
653    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h     !< green surface temperature (K) at horizontal walls
654    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_p   !< progn. green surface temperature (K) at horizontal walls
655    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h      !< near surface temperature (10cm) (K) at horizontal walls
656    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h_p    !< progn. near surface temperature (10cm) (K) at horizontal walls
657    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v
658    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_v_p
659    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_window_v
660    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_window_v_p
661    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_green_v
662    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_green_v_p
663    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_10cm_v
664    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  ::  t_surf_10cm_v_p
665#else
666    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h
667    REAL(wp), DIMENSION(:), POINTER                :: t_surf_h_p
668    REAL(wp), DIMENSION(:), POINTER                :: t_surf_window_h
669    REAL(wp), DIMENSION(:), POINTER                :: t_surf_window_h_p
670    REAL(wp), DIMENSION(:), POINTER                :: t_surf_green_h
671    REAL(wp), DIMENSION(:), POINTER                :: t_surf_green_h_p
672    REAL(wp), DIMENSION(:), POINTER                :: t_surf_10cm_h
673    REAL(wp), DIMENSION(:), POINTER                :: t_surf_10cm_h_p
674
675    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_1
676    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_h_2
677    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_1
678    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_h_2
679    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_1
680    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_h_2
681    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h_1
682    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_h_2
683
684    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v
685    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_v_p
686    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_window_v
687    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_window_v_p
688    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_green_v
689    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_green_v_p
690    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_10cm_v
691    TYPE(t_surf_vertical), DIMENSION(:), POINTER ::  t_surf_10cm_v_p
692
693    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_1
694    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_v_2
695    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_window_v_1
696    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_window_v_2
697    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_green_v_1
698    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_green_v_2
699    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_10cm_v_1
700    TYPE(t_surf_vertical), DIMENSION(0:3), TARGET  :: t_surf_10cm_v_2
701   
702#endif
703    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_av          !< average of wall surface temperature (K)
704    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_window_av   !< average of window surface temperature (K)
705    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_green_av    !< average of green wall surface temperature (K)
706    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET    :: t_surf_10cm_av    !< average of whole wall surface temperature (K)
707
708!-- Temporal tendencies for time stepping           
709    REAL(wp), DIMENSION(:), ALLOCATABLE            :: tt_surface_m       !< surface temperature tendency of wall (K)
710    REAL(wp), DIMENSION(:), ALLOCATABLE            :: tt_surface_window_m !< surface temperature tendency of window (K)
711    REAL(wp), DIMENSION(:), ALLOCATABLE            :: tt_surface_green_m !< surface temperature tendency of green wall (K)
712
713!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
714!-- Energy balance variables
715!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
716!-- parameters of the land, roof and wall surfaces
717
718#if defined( __nopointer )
719    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h             !< Wall temperature (K)
720    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_av          !< Average of t_wall
721    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_p           !< Prog. wall temperature (K)
722    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h           !< Window temperature (K)
723    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_av        !< Average of t_window
724    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_p         !< Prog. window temperature (K)
725    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h            !< Green temperature (K)
726    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_av         !< Average of t_green
727    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_p          !< Prog. green temperature (K)
728
729    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v             !< Wall temperature (K)
730    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av          !< Average of t_wall
731    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_p           !< Prog. wall temperature (K)
732    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v           !< Window temperature (K)
733    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_av        !< Average of t_window
734    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_p         !< Prog. window temperature (K)
735    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v            !< Green temperature (K)
736    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_av         !< Average of t_green
737    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_p          !< Prog. green temperature (K)
738#else
739    REAL(wp), DIMENSION(:,:), POINTER                :: t_wall_h, t_wall_h_p
740    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_wall_h_av, t_wall_h_1, t_wall_h_2
741    REAL(wp), DIMENSION(:,:), POINTER                :: t_window_h, t_window_h_p
742    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_window_h_av, t_window_h_1, t_window_h_2
743    REAL(wp), DIMENSION(:,:), POINTER                :: t_green_h, t_green_h_p
744    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET    :: t_green_h_av, t_green_h_1, t_green_h_2
745
746    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_wall_v, t_wall_v_p
747    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_wall_v_av, t_wall_v_1, t_wall_v_2
748    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_window_v, t_window_v_p
749    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_window_v_av, t_window_v_1, t_window_v_2
750    TYPE(t_wall_vertical), DIMENSION(:), POINTER   :: t_green_v, t_green_v_p
751    TYPE(t_wall_vertical), DIMENSION(0:3), TARGET  :: t_green_v_av, t_green_v_1, t_green_v_2
752#endif
753
754!-- Wall temporal tendencies for time stepping
755    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: tt_wall_m          !< t_wall prognostic array
756    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: tt_window_m        !< t_window prognostic array
757    REAL(wp), DIMENSION(:,:), ALLOCATABLE          :: tt_green_m         !< t_green prognostic array
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
[2920]762    INTEGER(iwp), PARAMETER                        :: n_surface_params = 9 !< number of parameters for each type of the wall
[2737]763    INTEGER(iwp), PARAMETER                        :: ialbedo  = 1         !< albedo of the surface
764    INTEGER(iwp), PARAMETER                        :: iemiss   = 2         !< emissivity of the surface
[2920]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 )
[2737]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
[2894]818    INTERFACE usm_rrd_local
819       MODULE PROCEDURE usm_rrd_local
820    END INTERFACE usm_rrd_local
[2737]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       
[2894]830    INTERFACE usm_wrd_local
831       MODULE PROCEDURE usm_wrd_local
832    END INTERFACE usm_wrd_local
[2737]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,&
[2920]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,              &
[2894]853           usm_wrd_local, usm_allocate_surface
[2737]854
855!-- Public parameters, constants and initial values
[2920]856    PUBLIC usm_anthropogenic_heat, usm_material_model,                          &
[2737]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
[3222]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 = 'none'
[2737]896        DO  l = 0, 3
897           ALLOCATE( surf_usm_v(l)%surface_types(1:surf_usm_v(l)%ns) )
[3222]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 = 'none'
[2737]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 
[2920]1272        INTEGER(iwp)                                       :: i, j, k, l, m, ids, idsint, iwl, istat
[2737]1273        CHARACTER (len=varnamelength)                      :: var, surfid
1274        INTEGER(iwp), PARAMETER                            :: nd = 5
1275        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER     :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
[2920]1276        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         :: dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
[2737]1277
1278!--     find the real name of the variable
[2906]1279        ids = -1
[2737]1280        var = TRIM(variable)
1281        DO i = 0, nd-1
1282            k = len(TRIM(var))
1283            j = len(TRIM(dirname(i)))
1284            IF ( var(k-j+1:k) == dirname(i) )  THEN
1285                ids = i
[2920]1286                idsint = dirint(ids)
[2737]1287                var = var(:k-j)
1288                EXIT
1289            ENDIF
1290        ENDDO
1291        IF ( ids == -1 )  THEN
1292            var = TRIM(variable)
1293        ENDIF
1294        IF ( var(1:11) == 'usm_t_wall_'  .AND.  len(TRIM(var)) >= 12 )  THEN
1295!--          wall layers
1296            READ(var(12:12), '(I1)', iostat=istat ) iwl
1297            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
1298                var = var(1:10)
1299            ELSE
1300!--             wrong wall layer index
1301                RETURN
1302            ENDIF
1303        ENDIF
1304        IF ( var(1:13) == 'usm_t_window_'  .AND.  len(TRIM(var)) >= 14 )  THEN
1305!--          wall layers
1306            READ(var(14:14), '(I1)', iostat=istat ) iwl
1307            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
1308                var = var(1:12)
1309            ELSE
1310!--             wrong window layer index
1311                RETURN
1312            ENDIF
1313        ENDIF
1314        IF ( var(1:12) == 'usm_t_green_'  .AND.  len(TRIM(var)) >= 13 )  THEN
1315!--          wall layers
1316            READ(var(13:13), '(I1)', iostat=istat ) iwl
1317            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
1318                var = var(1:11)
1319            ELSE
1320!--             wrong green layer index
1321                RETURN
1322            ENDIF
1323        ENDIF
1324
1325        IF ( mode == 'allocate' )  THEN
1326           
1327           SELECT CASE ( TRIM( var ) )
1328               
1329                CASE ( 'usm_rad_net' )
1330!--                 array of complete radiation balance
1331                    IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
1332                        ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
1333                        surf_usm_h%rad_net_av = 0.0_wp
1334                    ENDIF
1335                    DO  l = 0, 3
1336                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
1337                           ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
1338                           surf_usm_v(l)%rad_net_av = 0.0_wp
1339                       ENDIF
1340                    ENDDO
1341                   
1342                CASE ( 'usm_rad_insw' )
1343!--                 array of sw radiation falling to surface after i-th reflection
[2920]1344                    IF ( .NOT.  ALLOCATED(surfinsw_av) )  THEN
1345                        ALLOCATE( surfinsw_av(nsurfl) )
1346                        surfinsw_av = 0.0_wp
[2737]1347                    ENDIF
[2920]1348
[2737]1349                CASE ( 'usm_rad_inlw' )
1350!--                 array of lw radiation falling to surface after i-th reflection
[2920]1351                    IF ( .NOT.  ALLOCATED(surfinlw_av) )  THEN
1352                        ALLOCATE( surfinlw_av(nsurfl) )
1353                        surfinlw_av = 0.0_wp
[2737]1354                    ENDIF
1355
1356                CASE ( 'usm_rad_inswdir' )
1357!--                 array of direct sw radiation falling to surface from sun
1358                    IF ( .NOT.  ALLOCATED(surfinswdir_av) )  THEN
[2920]1359                        ALLOCATE( surfinswdir_av(nsurfl) )
[2737]1360                        surfinswdir_av = 0.0_wp
1361                    ENDIF
1362
1363                CASE ( 'usm_rad_inswdif' )
1364!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
1365                    IF ( .NOT.  ALLOCATED(surfinswdif_av) )  THEN
[2920]1366                        ALLOCATE( surfinswdif_av(nsurfl) )
[2737]1367                        surfinswdif_av = 0.0_wp
1368                    ENDIF
1369
1370                CASE ( 'usm_rad_inswref' )
1371!--                 array of sw radiation falling to surface from reflections
1372                    IF ( .NOT.  ALLOCATED(surfinswref_av) )  THEN
[2920]1373                        ALLOCATE( surfinswref_av(nsurfl) )
[2737]1374                        surfinswref_av = 0.0_wp
1375                    ENDIF
1376
1377                CASE ( 'usm_rad_inlwdif' )
1378!--                 array of sw radiation falling to surface after i-th reflection
1379                   IF ( .NOT.  ALLOCATED(surfinlwdif_av) )  THEN
[2920]1380                        ALLOCATE( surfinlwdif_av(nsurfl) )
[2737]1381                        surfinlwdif_av = 0.0_wp
1382                    ENDIF
1383
1384                CASE ( 'usm_rad_inlwref' )
1385!--                 array of lw radiation falling to surface from reflections
1386                    IF ( .NOT.  ALLOCATED(surfinlwref_av) )  THEN
[2920]1387                        ALLOCATE( surfinlwref_av(nsurfl) )
[2737]1388                        surfinlwref_av = 0.0_wp
1389                    ENDIF
1390
1391                CASE ( 'usm_rad_outsw' )
1392!--                 array of sw radiation emitted from surface after i-th reflection
1393                    IF ( .NOT.  ALLOCATED(surfoutsw_av) )  THEN
[2920]1394                        ALLOCATE( surfoutsw_av(nsurfl) )
[2737]1395                        surfoutsw_av = 0.0_wp
1396                    ENDIF
1397
1398                CASE ( 'usm_rad_outlw' )
1399!--                 array of lw radiation emitted from surface after i-th reflection
1400                    IF ( .NOT.  ALLOCATED(surfoutlw_av) )  THEN
[2920]1401                        ALLOCATE( surfoutlw_av(nsurfl) )
[2737]1402                        surfoutlw_av = 0.0_wp
1403                    ENDIF
1404                CASE ( 'usm_rad_ressw' )
1405!--                 array of residua of sw radiation absorbed in surface after last reflection
1406                    IF ( .NOT.  ALLOCATED(surfins_av) )  THEN
[2920]1407                        ALLOCATE( surfins_av(nsurfl) )
[2737]1408                        surfins_av = 0.0_wp
1409                    ENDIF
1410                                   
1411                CASE ( 'usm_rad_reslw' )
1412!--                 array of residua of lw radiation absorbed in surface after last reflection
1413                    IF ( .NOT.  ALLOCATED(surfinl_av) )  THEN
[2920]1414                        ALLOCATE( surfinl_av(nsurfl) )
[2737]1415                        surfinl_av = 0.0_wp
1416                    ENDIF
1417                                   
1418                CASE ( 'usm_rad_hf' )
1419!--                 array of heat flux from radiation for surfaces after i-th reflection
1420                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
1421                        ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
1422                        surf_usm_h%surfhf_av = 0.0_wp
1423                    ENDIF
1424                    DO  l = 0, 3
1425                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
1426                           ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
1427                           surf_usm_v(l)%surfhf_av = 0.0_wp
1428                       ENDIF
1429                    ENDDO
1430
1431                CASE ( 'usm_wshf' )
1432!--                 array of sensible heat flux from surfaces
1433!--                 land surfaces
1434                    IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
1435                        ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
1436                        surf_usm_h%wshf_eb_av = 0.0_wp
1437                    ENDIF
1438                    DO  l = 0, 3
1439                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
1440                           ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
1441                           surf_usm_v(l)%wshf_eb_av = 0.0_wp
1442                       ENDIF
1443                    ENDDO
[2797]1444!
1445!--             Please note, the following output quantities belongs to the
1446!--             individual tile fractions - ground heat flux at wall-, window-,
1447!--             and green fraction. Aggregated ground-heat flux is treated
1448!--             accordingly in average_3d_data, sum_up_3d_data, etc..
[2737]1449                CASE ( 'usm_wghf' )
1450!--                 array of heat flux from ground (wall, roof, land)
1451                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
1452                        ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
1453                        surf_usm_h%wghf_eb_av = 0.0_wp
1454                    ENDIF
1455                    DO  l = 0, 3
1456                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
1457                           ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
1458                           surf_usm_v(l)%wghf_eb_av = 0.0_wp
1459                       ENDIF
1460                    ENDDO
1461
1462                CASE ( 'usm_wghf_window' )
1463!--                 array of heat flux from window ground (wall, roof, land)
1464                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
1465                        ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
1466                        surf_usm_h%wghf_eb_window_av = 0.0_wp
1467                    ENDIF
1468                    DO  l = 0, 3
1469                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
1470                           ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
1471                           surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
1472                       ENDIF
1473                    ENDDO
1474
1475                CASE ( 'usm_wghf_green' )
1476!--                 array of heat flux from green ground (wall, roof, land)
1477                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
1478                        ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
1479                        surf_usm_h%wghf_eb_green_av = 0.0_wp
1480                    ENDIF
1481                    DO  l = 0, 3
1482                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
1483                           ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
1484                           surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
1485                       ENDIF
1486                    ENDDO
1487
1488                CASE ( 'usm_iwghf' )
1489!--                 array of heat flux from indoor ground (wall, roof, land)
1490                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
1491                        ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
1492                        surf_usm_h%iwghf_eb_av = 0.0_wp
1493                    ENDIF
1494                    DO  l = 0, 3
1495                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
1496                           ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
1497                           surf_usm_v(l)%iwghf_eb_av = 0.0_wp
1498                       ENDIF
1499                    ENDDO
1500
1501                CASE ( 'usm_iwghf_window' )
1502!--                 array of heat flux from indoor window ground (wall, roof, land)
1503                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
1504                        ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
1505                        surf_usm_h%iwghf_eb_window_av = 0.0_wp
1506                    ENDIF
1507                    DO  l = 0, 3
1508                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
1509                           ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
1510                           surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
1511                       ENDIF
1512                    ENDDO
1513                   
1514                CASE ( 'usm_t_surf' )
1515!--                 surface temperature for surfaces
1516                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
1517                        ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
1518                        surf_usm_h%t_surf_av = 0.0_wp
1519                    ENDIF
1520                    DO  l = 0, 3
1521                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_av) )  THEN
1522                           ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
1523                           surf_usm_v(l)%t_surf_av = 0.0_wp
1524                       ENDIF
1525                    ENDDO
1526
1527                CASE ( 'usm_t_surf_window' )
1528!--                 surface temperature for window surfaces
1529                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
1530                        ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
1531                        surf_usm_h%t_surf_window_av = 0.0_wp
1532                    ENDIF
1533                    DO  l = 0, 3
1534                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
1535                           ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
1536                           surf_usm_v(l)%t_surf_window_av = 0.0_wp
1537                       ENDIF
1538                    ENDDO
1539                   
1540                CASE ( 'usm_t_surf_green' )
1541!--                 surface temperature for green surfaces
1542                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
1543                        ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
1544                        surf_usm_h%t_surf_green_av = 0.0_wp
1545                    ENDIF
1546                    DO  l = 0, 3
1547                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
1548                           ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
1549                           surf_usm_v(l)%t_surf_green_av = 0.0_wp
1550                       ENDIF
1551                    ENDDO
1552               
1553                CASE ( 'usm_t_surf_10cm' )
1554!--                 near surface temperature for whole surfaces
1555                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
1556                        ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
1557                        surf_usm_h%t_surf_10cm_av = 0.0_wp
1558                    ENDIF
1559                    DO  l = 0, 3
1560                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) )  THEN
1561                           ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) )
1562                           surf_usm_v(l)%t_surf_10cm_av = 0.0_wp
1563                       ENDIF
1564                    ENDDO
1565
1566                CASE ( 'usm_t_wall' )
1567!--                 wall temperature for iwl layer of walls and land
1568                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
1569                        ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1570                        surf_usm_h%t_wall_av = 0.0_wp
1571                    ENDIF
1572                    DO  l = 0, 3
1573                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
1574                           ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1575                           surf_usm_v(l)%t_wall_av = 0.0_wp
1576                       ENDIF
1577                    ENDDO
1578
1579                CASE ( 'usm_t_window' )
1580!--                 window temperature for iwl layer of walls and land
1581                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
1582                        ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1583                        surf_usm_h%t_window_av = 0.0_wp
1584                    ENDIF
1585                    DO  l = 0, 3
1586                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
1587                           ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1588                           surf_usm_v(l)%t_window_av = 0.0_wp
1589                       ENDIF
1590                    ENDDO
1591
1592                CASE ( 'usm_t_green' )
1593!--                 green temperature for iwl layer of walls and land
1594                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
1595                        ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
1596                        surf_usm_h%t_green_av = 0.0_wp
1597                    ENDIF
1598                    DO  l = 0, 3
1599                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
1600                           ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
1601                           surf_usm_v(l)%t_green_av = 0.0_wp
1602                       ENDIF
1603                    ENDDO
1604
1605               CASE DEFAULT
1606                   CONTINUE
1607
1608           END SELECT
1609
1610        ELSEIF ( mode == 'sum' )  THEN
1611           
1612           SELECT CASE ( TRIM( var ) )
1613               
1614                CASE ( 'usm_rad_net' )
1615!--                 array of complete radiation balance
1616                    DO  m = 1, surf_usm_h%ns
1617                       surf_usm_h%rad_net_av(m) =                              &
1618                                          surf_usm_h%rad_net_av(m) +           &
1619                                          surf_usm_h%rad_net_l(m)
1620                    ENDDO
1621                    DO  l = 0, 3
1622                       DO  m = 1, surf_usm_v(l)%ns
1623                          surf_usm_v(l)%rad_net_av(m) =                        &
1624                                          surf_usm_v(l)%rad_net_av(m) +        &
1625                                          surf_usm_v(l)%rad_net_l(m)
1626                       ENDDO
1627                    ENDDO
1628                   
1629                CASE ( 'usm_rad_insw' )
1630!--                 array of sw radiation falling to surface after i-th reflection
[2920]1631                    DO l = 1, nsurfl
1632                        IF ( surfl(id,l) == idsint )  THEN
[2737]1633                            surfinsw_av(l) = surfinsw_av(l) + surfinsw(l)
1634                        ENDIF
1635                    ENDDO
1636                             
1637                CASE ( 'usm_rad_inlw' )
1638!--                 array of lw radiation falling to surface after i-th reflection
[2920]1639                    DO l = 1, nsurfl
1640                        IF ( surfl(id,l) == idsint )  THEN
[2737]1641                            surfinlw_av(l) = surfinlw_av(l) + surfinlw(l)
1642                        ENDIF
1643                    ENDDO
1644                   
1645                CASE ( 'usm_rad_inswdir' )
1646!--                 array of direct sw radiation falling to surface from sun
[2920]1647                    DO l = 1, nsurfl
1648                        IF ( surfl(id,l) == idsint )  THEN
[2737]1649                            surfinswdir_av(l) = surfinswdir_av(l) + surfinswdir(l)
1650                        ENDIF
1651                    ENDDO
1652                   
1653                CASE ( 'usm_rad_inswdif' )
1654!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
[2920]1655                    DO l = 1, nsurfl
1656                        IF ( surfl(id,l) == idsint )  THEN
[2737]1657                            surfinswdif_av(l) = surfinswdif_av(l) + surfinswdif(l)
1658                        ENDIF
1659                    ENDDO
1660                   
1661                CASE ( 'usm_rad_inswref' )
1662!--                 array of sw radiation falling to surface from reflections
[2920]1663                    DO l = 1, nsurfl
1664                        IF ( surfl(id,l) == idsint )  THEN
[2737]1665                            surfinswref_av(l) = surfinswref_av(l) + surfinsw(l) - &
1666                                                surfinswdir(l) - surfinswdif(l)
1667                        ENDIF
1668                    ENDDO
1669
1670                   
1671                CASE ( 'usm_rad_inlwdif' )
1672!--                 array of sw radiation falling to surface after i-th reflection
[2920]1673                    DO l = 1, nsurfl
1674                        IF ( surfl(id,l) == idsint )  THEN
1675                            surfinlwdif_av(l) = surfinlwdif_av(l) + surfinlwdif(l)
[2737]1676                        ENDIF
1677                    ENDDO
1678!                     
1679                CASE ( 'usm_rad_inlwref' )
1680!--                 array of lw radiation falling to surface from reflections
[2920]1681                    DO l = 1, nsurfl
1682                        IF ( surfl(id,l) == idsint )  THEN
1683                            surfinlwref_av(l) = surfinlwref_av(l) + &
1684                                                surfinlw(l) - surfinlwdif(l)
[2737]1685                        ENDIF
1686                    ENDDO
1687                   
1688                CASE ( 'usm_rad_outsw' )
1689!--                 array of sw radiation emitted from surface after i-th reflection
[2920]1690                    DO l = 1, nsurfl
1691                        IF ( surfl(id,l) == idsint )  THEN
1692                            surfoutsw_av(l) = surfoutsw_av(l) + surfoutsw(l)
[2737]1693                        ENDIF
1694                    ENDDO
1695                   
1696                CASE ( 'usm_rad_outlw' )
1697!--                 array of lw radiation emitted from surface after i-th reflection
[2920]1698                    DO l = 1, nsurfl
1699                        IF ( surfl(id,l) == idsint )  THEN
1700                            surfoutlw_av(l) = surfoutlw_av(l) + surfoutlw(l)
[2737]1701                        ENDIF
1702                    ENDDO
1703                   
1704                CASE ( 'usm_rad_ressw' )
1705!--                 array of residua of sw radiation absorbed in surface after last reflection
[2920]1706                    DO l = 1, nsurfl
1707                        IF ( surfl(id,l) == idsint )  THEN
1708                            surfins_av(l) = surfins_av(l) + surfins(l)
[2737]1709                        ENDIF
1710                    ENDDO
1711                                   
1712                CASE ( 'usm_rad_reslw' )
1713!--                 array of residua of lw radiation absorbed in surface after last reflection
[2920]1714                    DO l = 1, nsurfl
1715                        IF ( surfl(id,l) == idsint )  THEN
1716                            surfinl_av(l) = surfinl_av(l) + surfinl(l)
[2737]1717                        ENDIF
1718                    ENDDO
1719                   
1720                CASE ( 'usm_rad_hf' )
1721!--                 array of heat flux from radiation for surfaces after i-th reflection
1722                    DO  m = 1, surf_usm_h%ns
1723                       surf_usm_h%surfhf_av(m) =                               &
1724                                          surf_usm_h%surfhf_av(m) +            &
1725                                          surf_usm_h%surfhf(m)
1726                    ENDDO
1727                    DO  l = 0, 3
1728                       DO  m = 1, surf_usm_v(l)%ns
1729                          surf_usm_v(l)%surfhf_av(m) =                         &
1730                                          surf_usm_v(l)%surfhf_av(m) +         &
1731                                          surf_usm_v(l)%surfhf(m)
1732                       ENDDO
1733                    ENDDO
1734                   
1735                CASE ( 'usm_wshf' )
1736!--                 array of sensible heat flux from surfaces (land, roof, wall)
1737                    DO  m = 1, surf_usm_h%ns
1738                       surf_usm_h%wshf_eb_av(m) =                              &
1739                                          surf_usm_h%wshf_eb_av(m) +           &
1740                                          surf_usm_h%wshf_eb(m)
1741                    ENDDO
1742                    DO  l = 0, 3
1743                       DO  m = 1, surf_usm_v(l)%ns
1744                          surf_usm_v(l)%wshf_eb_av(m) =                        &
1745                                          surf_usm_v(l)%wshf_eb_av(m) +        &
1746                                          surf_usm_v(l)%wshf_eb(m)
1747                       ENDDO
1748                    ENDDO
1749                   
1750                CASE ( 'usm_wghf' )
1751!--                 array of heat flux from ground (wall, roof, land)
1752                    DO  m = 1, surf_usm_h%ns
1753                       surf_usm_h%wghf_eb_av(m) =                              &
1754                                          surf_usm_h%wghf_eb_av(m) +           &
1755                                          surf_usm_h%wghf_eb(m)
1756                    ENDDO
1757                    DO  l = 0, 3
1758                       DO  m = 1, surf_usm_v(l)%ns
1759                          surf_usm_v(l)%wghf_eb_av(m) =                        &
1760                                          surf_usm_v(l)%wghf_eb_av(m) +        &
1761                                          surf_usm_v(l)%wghf_eb(m)
1762                       ENDDO
1763                    ENDDO
1764                   
1765                CASE ( 'usm_wghf_window' )
1766!--                 array of heat flux from window ground (wall, roof, land)
1767                    DO  m = 1, surf_usm_h%ns
1768                       surf_usm_h%wghf_eb_window_av(m) =                              &
1769                                          surf_usm_h%wghf_eb_window_av(m) +           &
1770                                          surf_usm_h%wghf_eb_window(m)
1771                    ENDDO
1772                    DO  l = 0, 3
1773                       DO  m = 1, surf_usm_v(l)%ns
1774                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
1775                                          surf_usm_v(l)%wghf_eb_window_av(m) +        &
1776                                          surf_usm_v(l)%wghf_eb_window(m)
1777                       ENDDO
1778                    ENDDO
1779
1780                CASE ( 'usm_wghf_green' )
1781!--                 array of heat flux from green ground (wall, roof, land)
1782                    DO  m = 1, surf_usm_h%ns
1783                       surf_usm_h%wghf_eb_green_av(m) =                              &
1784                                          surf_usm_h%wghf_eb_green_av(m) +           &
1785                                          surf_usm_h%wghf_eb_green(m)
1786                    ENDDO
1787                    DO  l = 0, 3
1788                       DO  m = 1, surf_usm_v(l)%ns
1789                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
1790                                          surf_usm_v(l)%wghf_eb_green_av(m) +        &
1791                                          surf_usm_v(l)%wghf_eb_green(m)
1792                       ENDDO
1793                    ENDDO
1794                   
1795                CASE ( 'usm_iwghf' )
1796!--                 array of heat flux from indoor ground (wall, roof, land)
1797                    DO  m = 1, surf_usm_h%ns
1798                       surf_usm_h%iwghf_eb_av(m) =                              &
1799                                          surf_usm_h%iwghf_eb_av(m) +           &
1800                                          surf_usm_h%iwghf_eb(m)
1801                    ENDDO
1802                    DO  l = 0, 3
1803                       DO  m = 1, surf_usm_v(l)%ns
1804                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
1805                                          surf_usm_v(l)%iwghf_eb_av(m) +        &
1806                                          surf_usm_v(l)%iwghf_eb(m)
1807                       ENDDO
1808                    ENDDO
1809                   
1810                CASE ( 'usm_iwghf_window' )
1811!--                 array of heat flux from indoor window ground (wall, roof, land)
1812                    DO  m = 1, surf_usm_h%ns
1813                       surf_usm_h%iwghf_eb_window_av(m) =                              &
1814                                          surf_usm_h%iwghf_eb_window_av(m) +           &
1815                                          surf_usm_h%iwghf_eb_window(m)
1816                    ENDDO
1817                    DO  l = 0, 3
1818                       DO  m = 1, surf_usm_v(l)%ns
1819                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
1820                                          surf_usm_v(l)%iwghf_eb_window_av(m) +        &
1821                                          surf_usm_v(l)%iwghf_eb_window(m)
1822                       ENDDO
1823                    ENDDO
1824                   
1825                CASE ( 'usm_t_surf' )
1826!--                 surface temperature for surfaces
1827                    DO  m = 1, surf_usm_h%ns
1828                       surf_usm_h%t_surf_av(m) =                               & 
1829                                          surf_usm_h%t_surf_av(m) +            &
1830                                          t_surf_h(m)
1831                    ENDDO
1832                    DO  l = 0, 3
1833                       DO  m = 1, surf_usm_v(l)%ns
1834                          surf_usm_v(l)%t_surf_av(m) =                         &
1835                                          surf_usm_v(l)%t_surf_av(m) +         &
1836                                          t_surf_v(l)%t(m)
1837                       ENDDO
1838                    ENDDO
1839                   
1840                CASE ( 'usm_t_surf_window' )
1841!--                 surface temperature for window surfaces
1842                    DO  m = 1, surf_usm_h%ns
1843                       surf_usm_h%t_surf_window_av(m) =                               & 
1844                                          surf_usm_h%t_surf_window_av(m) +            &
1845                                          t_surf_window_h(m)
1846                    ENDDO
1847                    DO  l = 0, 3
1848                       DO  m = 1, surf_usm_v(l)%ns
1849                          surf_usm_v(l)%t_surf_window_av(m) =                         &
1850                                          surf_usm_v(l)%t_surf_window_av(m) +         &
1851                                          t_surf_window_v(l)%t(m)
1852                       ENDDO
1853                    ENDDO
1854                   
1855                CASE ( 'usm_t_surf_green' )
1856!--                 surface temperature for green surfaces
1857                    DO  m = 1, surf_usm_h%ns
1858                       surf_usm_h%t_surf_green_av(m) =                               & 
1859                                          surf_usm_h%t_surf_green_av(m) +            &
1860                                          t_surf_green_h(m)
1861                    ENDDO
1862                    DO  l = 0, 3
1863                       DO  m = 1, surf_usm_v(l)%ns
1864                          surf_usm_v(l)%t_surf_green_av(m) =                         &
1865                                          surf_usm_v(l)%t_surf_green_av(m) +         &
1866                                          t_surf_green_v(l)%t(m)
1867                       ENDDO
1868                    ENDDO
1869               
1870                CASE ( 'usm_t_surf_10cm' )
1871!--                 near surface temperature for whole surfaces
1872                    DO  m = 1, surf_usm_h%ns
1873                       surf_usm_h%t_surf_10cm_av(m) =                               & 
1874                                          surf_usm_h%t_surf_10cm_av(m) +            &
1875                                          t_surf_10cm_h(m)
1876                    ENDDO
1877                    DO  l = 0, 3
1878                       DO  m = 1, surf_usm_v(l)%ns
1879                          surf_usm_v(l)%t_surf_10cm_av(m) =                         &
1880                                          surf_usm_v(l)%t_surf_10cm_av(m) +         &
1881                                          t_surf_10cm_v(l)%t(m)
1882                       ENDDO
1883                    ENDDO
1884
1885                   
1886                CASE ( 'usm_t_wall' )
1887!--                 wall temperature for  iwl layer of walls and land
1888                    DO  m = 1, surf_usm_h%ns
1889                       surf_usm_h%t_wall_av(iwl,m) =                           &
1890                                          surf_usm_h%t_wall_av(iwl,m) +        &
1891                                          t_wall_h(iwl,m)
1892                    ENDDO
1893                    DO  l = 0, 3
1894                       DO  m = 1, surf_usm_v(l)%ns
1895                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
1896                                          surf_usm_v(l)%t_wall_av(iwl,m) +     &
1897                                          t_wall_v(l)%t(iwl,m)
1898                       ENDDO
1899                    ENDDO
1900                   
1901                CASE ( 'usm_t_window' )
1902!--                 window temperature for  iwl layer of walls and land
1903                    DO  m = 1, surf_usm_h%ns
1904                       surf_usm_h%t_window_av(iwl,m) =                           &
1905                                          surf_usm_h%t_window_av(iwl,m) +        &
1906                                          t_window_h(iwl,m)
1907                    ENDDO
1908                    DO  l = 0, 3
1909                       DO  m = 1, surf_usm_v(l)%ns
1910                          surf_usm_v(l)%t_window_av(iwl,m) =                     &
1911                                          surf_usm_v(l)%t_window_av(iwl,m) +     &
1912                                          t_window_v(l)%t(iwl,m)
1913                       ENDDO
1914                    ENDDO
1915
1916                CASE ( 'usm_t_green' )
1917!--                 green temperature for  iwl layer of walls and land
1918                    DO  m = 1, surf_usm_h%ns
1919                       surf_usm_h%t_green_av(iwl,m) =                           &
1920                                          surf_usm_h%t_green_av(iwl,m) +        &
1921                                          t_green_h(iwl,m)
1922                    ENDDO
1923                    DO  l = 0, 3
1924                       DO  m = 1, surf_usm_v(l)%ns
1925                          surf_usm_v(l)%t_green_av(iwl,m) =                     &
1926                                          surf_usm_v(l)%t_green_av(iwl,m) +     &
1927                                          t_green_v(l)%t(iwl,m)
1928                       ENDDO
1929                    ENDDO
1930
1931                CASE DEFAULT
1932                    CONTINUE
1933
1934           END SELECT
1935
1936        ELSEIF ( mode == 'average' )  THEN
1937           
1938           SELECT CASE ( TRIM( var ) )
1939               
1940                CASE ( 'usm_rad_net' )
1941!--                 array of complete radiation balance
1942                    DO  m = 1, surf_usm_h%ns
1943                       surf_usm_h%rad_net_av(m) =                              &
1944                                          surf_usm_h%rad_net_av(m) /           &
1945                                          REAL( average_count_3d, kind=wp )
1946                    ENDDO
1947                    DO  l = 0, 3
1948                       DO  m = 1, surf_usm_v(l)%ns
1949                          surf_usm_v(l)%rad_net_av(m) =                        &
1950                                          surf_usm_v(l)%rad_net_av(m) /        &
1951                                          REAL( average_count_3d, kind=wp )
1952                       ENDDO
1953                    ENDDO
1954                   
1955                CASE ( 'usm_rad_insw' )
1956!--                 array of sw radiation falling to surface after i-th reflection
[2920]1957                    DO l = 1, nsurfl
1958                        IF ( surfl(id,l) == idsint )  THEN
[2737]1959                            surfinsw_av(l) = surfinsw_av(l) / REAL( average_count_3d, kind=wp )
1960                        ENDIF
1961                    ENDDO
1962                             
1963                CASE ( 'usm_rad_inlw' )
1964!--                 array of lw radiation falling to surface after i-th reflection
[2920]1965                    DO l = 1, nsurfl
1966                        IF ( surfl(id,l) == idsint )  THEN
[2737]1967                            surfinlw_av(l) = surfinlw_av(l) / REAL( average_count_3d, kind=wp )
1968                        ENDIF
1969                    ENDDO
1970                   
1971                CASE ( 'usm_rad_inswdir' )
1972!--                 array of direct sw radiation falling to surface from sun
[2920]1973                    DO l = 1, nsurfl
1974                        IF ( surfl(id,l) == idsint )  THEN
[2737]1975                            surfinswdir_av(l) = surfinswdir_av(l) / REAL( average_count_3d, kind=wp )
1976                        ENDIF
1977                    ENDDO
1978                   
1979                CASE ( 'usm_rad_inswdif' )
1980!--                 array of difusion sw radiation falling to surface from sky and borders of the domain
[2920]1981                    DO l = 1, nsurfl
1982                        IF ( surfl(id,l) == idsint )  THEN
[2737]1983                            surfinswdif_av(l) = surfinswdif_av(l) / REAL( average_count_3d, kind=wp )
1984                        ENDIF
1985                    ENDDO
1986                   
1987                CASE ( 'usm_rad_inswref' )
1988!--                 array of sw radiation falling to surface from reflections
[2920]1989                    DO l = 1, nsurfl
1990                        IF ( surfl(id,l) == idsint )  THEN
[2737]1991                            surfinswref_av(l) = surfinswref_av(l) / REAL( average_count_3d, kind=wp )
1992                        ENDIF
1993                    ENDDO
1994                   
1995                CASE ( 'usm_rad_inlwdif' )
1996!--                 array of sw radiation falling to surface after i-th reflection
[2920]1997                    DO l = 1, nsurfl
1998                        IF ( surfl(id,l) == idsint )  THEN
[2737]1999                            surfinlwdif_av(l) = surfinlwdif_av(l) / REAL( average_count_3d, kind=wp )
2000                        ENDIF
2001                    ENDDO
2002                   
2003                CASE ( 'usm_rad_inlwref' )
2004!--                 array of lw radiation falling to surface from reflections
[2920]2005                    DO l = 1, nsurfl
2006                        IF ( surfl(id,l) == idsint )  THEN
[2737]2007                            surfinlwref_av(l) = surfinlwref_av(l) / REAL( average_count_3d, kind=wp )
2008                        ENDIF
2009                    ENDDO
2010                   
2011                CASE ( 'usm_rad_outsw' )
2012!--                 array of sw radiation emitted from surface after i-th reflection
[2920]2013                    DO l = 1, nsurfl
2014                        IF ( surfl(id,l) == idsint )  THEN
[2737]2015                            surfoutsw_av(l) = surfoutsw_av(l) / REAL( average_count_3d, kind=wp )
2016                        ENDIF
2017                    ENDDO
2018                   
2019                CASE ( 'usm_rad_outlw' )
2020!--                 array of lw radiation emitted from surface after i-th reflection
[2920]2021                    DO l = 1, nsurfl
2022                        IF ( surfl(id,l) == idsint )  THEN
[2737]2023                            surfoutlw_av(l) = surfoutlw_av(l) / REAL( average_count_3d, kind=wp )
2024                        ENDIF
2025                    ENDDO
2026                   
2027                CASE ( 'usm_rad_ressw' )
2028!--                 array of residua of sw radiation absorbed in surface after last reflection
[2920]2029                    DO l = 1, nsurfl
2030                        IF ( surfl(id,l) == idsint )  THEN
[2737]2031                            surfins_av(l) = surfins_av(l) / REAL( average_count_3d, kind=wp )
2032                        ENDIF
2033                    ENDDO
2034                                   
2035                CASE ( 'usm_rad_reslw' )
2036!--                 array of residua of lw radiation absorbed in surface after last reflection
[2920]2037                    DO l = 1, nsurfl
2038                        IF ( surfl(id,l) == idsint )  THEN
[2737]2039                            surfinl_av(l) = surfinl_av(l) / REAL( average_count_3d, kind=wp )
2040                        ENDIF
2041                    ENDDO
2042                   
2043                CASE ( 'usm_rad_hf' )
2044!--                 array of heat flux from radiation for surfaces after i-th reflection
2045                    DO  m = 1, surf_usm_h%ns
2046                       surf_usm_h%surfhf_av(m) =                               &
2047                                          surf_usm_h%surfhf_av(m) /            &
2048                                          REAL( average_count_3d, kind=wp )
2049                    ENDDO
2050                    DO  l = 0, 3
2051                       DO  m = 1, surf_usm_v(l)%ns
2052                          surf_usm_v(l)%surfhf_av(m) =                         &
2053                                          surf_usm_v(l)%surfhf_av(m) /         &
2054                                          REAL( average_count_3d, kind=wp )
2055                       ENDDO
2056                    ENDDO
2057                   
2058                CASE ( 'usm_wshf' )
2059!--                 array of sensible heat flux from surfaces (land, roof, wall)
2060                    DO  m = 1, surf_usm_h%ns
2061                       surf_usm_h%wshf_eb_av(m) =                              &
2062                                          surf_usm_h%wshf_eb_av(m) /           &
2063                                          REAL( average_count_3d, kind=wp )
2064                    ENDDO
2065                    DO  l = 0, 3
2066                       DO  m = 1, surf_usm_v(l)%ns
2067                          surf_usm_v(l)%wshf_eb_av(m) =                        &
2068                                          surf_usm_v(l)%wshf_eb_av(m) /        &
2069                                          REAL( average_count_3d, kind=wp )
2070                       ENDDO
2071                    ENDDO
2072                   
2073                CASE ( 'usm_wghf' )
2074!--                 array of heat flux from ground (wall, roof, land)
2075                    DO  m = 1, surf_usm_h%ns
2076                       surf_usm_h%wghf_eb_av(m) =                              &
2077                                          surf_usm_h%wghf_eb_av(m) /           &
2078                                          REAL( average_count_3d, kind=wp )
2079                    ENDDO
2080                    DO  l = 0, 3
2081                       DO  m = 1, surf_usm_v(l)%ns
2082                          surf_usm_v(l)%wghf_eb_av(m) =                        &
2083                                          surf_usm_v(l)%wghf_eb_av(m) /        &
2084                                          REAL( average_count_3d, kind=wp )
2085                       ENDDO
2086                    ENDDO
2087                   
2088                CASE ( 'usm_wghf_window' )
2089!--                 array of heat flux from window ground (wall, roof, land)
2090                    DO  m = 1, surf_usm_h%ns
2091                       surf_usm_h%wghf_eb_window_av(m) =                              &
2092                                          surf_usm_h%wghf_eb_window_av(m) /           &
2093                                          REAL( average_count_3d, kind=wp )
2094                    ENDDO
2095                    DO  l = 0, 3
2096                       DO  m = 1, surf_usm_v(l)%ns
2097                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
2098                                          surf_usm_v(l)%wghf_eb_window_av(m) /        &
2099                                          REAL( average_count_3d, kind=wp )
2100                       ENDDO
2101                    ENDDO
2102
2103                CASE ( 'usm_wghf_green' )
2104!--                 array of heat flux from green ground (wall, roof, land)
2105                    DO  m = 1, surf_usm_h%ns
2106                       surf_usm_h%wghf_eb_green_av(m) =                              &
2107                                          surf_usm_h%wghf_eb_green_av(m) /           &
2108                                          REAL( average_count_3d, kind=wp )
2109                    ENDDO
2110                    DO  l = 0, 3
2111                       DO  m = 1, surf_usm_v(l)%ns
2112                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
2113                                          surf_usm_v(l)%wghf_eb_green_av(m) /        &
2114                                          REAL( average_count_3d, kind=wp )
2115                       ENDDO
2116                    ENDDO
2117
2118                CASE ( 'usm_iwghf' )
2119!--                 array of heat flux from indoor ground (wall, roof, land)
2120                    DO  m = 1, surf_usm_h%ns
2121                       surf_usm_h%iwghf_eb_av(m) =                              &
2122                                          surf_usm_h%iwghf_eb_av(m) /           &
2123                                          REAL( average_count_3d, kind=wp )
2124                    ENDDO
2125                    DO  l = 0, 3
2126                       DO  m = 1, surf_usm_v(l)%ns
2127                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
2128                                          surf_usm_v(l)%iwghf_eb_av(m) /        &
2129                                          REAL( average_count_3d, kind=wp )
2130                       ENDDO
2131                    ENDDO
2132                   
2133                CASE ( 'usm_iwghf_window' )
2134!--                 array of heat flux from indoor window ground (wall, roof, land)
2135                    DO  m = 1, surf_usm_h%ns
2136                       surf_usm_h%iwghf_eb_window_av(m) =                              &
2137                                          surf_usm_h%iwghf_eb_window_av(m) /           &
2138                                          REAL( average_count_3d, kind=wp )
2139                    ENDDO
2140                    DO  l = 0, 3
2141                       DO  m = 1, surf_usm_v(l)%ns
2142                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
2143                                          surf_usm_v(l)%iwghf_eb_window_av(m) /        &
2144                                          REAL( average_count_3d, kind=wp )
2145                       ENDDO
2146                    ENDDO
2147                   
2148                CASE ( 'usm_t_surf' )
2149!--                 surface temperature for surfaces
2150                    DO  m = 1, surf_usm_h%ns
2151                       surf_usm_h%t_surf_av(m) =                               & 
2152                                          surf_usm_h%t_surf_av(m) /            &
2153                                          REAL( average_count_3d, kind=wp )
2154                    ENDDO
2155                    DO  l = 0, 3
2156                       DO  m = 1, surf_usm_v(l)%ns
2157                          surf_usm_v(l)%t_surf_av(m) =                         &
2158                                          surf_usm_v(l)%t_surf_av(m) /         &
2159                                          REAL( average_count_3d, kind=wp )
2160                       ENDDO
2161                    ENDDO
2162                   
2163                CASE ( 'usm_t_surf_window' )
2164!--                 surface temperature for window surfaces
2165                    DO  m = 1, surf_usm_h%ns
2166                       surf_usm_h%t_surf_window_av(m) =                               & 
2167                                          surf_usm_h%t_surf_window_av(m) /            &
2168                                          REAL( average_count_3d, kind=wp )
2169                    ENDDO
2170                    DO  l = 0, 3
2171                       DO  m = 1, surf_usm_v(l)%ns
2172                          surf_usm_v(l)%t_surf_window_av(m) =                         &
2173                                          surf_usm_v(l)%t_surf_window_av(m) /         &
2174                                          REAL( average_count_3d, kind=wp )
2175                       ENDDO
2176                    ENDDO
2177                   
2178                CASE ( 'usm_t_surf_green' )
2179!--                 surface temperature for green surfaces
2180                    DO  m = 1, surf_usm_h%ns
2181                       surf_usm_h%t_surf_green_av(m) =                               & 
2182                                          surf_usm_h%t_surf_green_av(m) /            &
2183                                          REAL( average_count_3d, kind=wp )
2184                    ENDDO
2185                    DO  l = 0, 3
2186                       DO  m = 1, surf_usm_v(l)%ns
2187                          surf_usm_v(l)%t_surf_green_av(m) =                         &
2188                                          surf_usm_v(l)%t_surf_green_av(m) /         &
2189                                          REAL( average_count_3d, kind=wp )
2190                       ENDDO
2191                    ENDDO
2192                   
2193                CASE ( 'usm_t_surf_10cm' )
2194!--                 near surface temperature for whole surfaces
2195                    DO  m = 1, surf_usm_h%ns
2196                       surf_usm_h%t_surf_10cm_av(m) =                               & 
2197                                          surf_usm_h%t_surf_10cm_av(m) /            &
2198                                          REAL( average_count_3d, kind=wp )
2199                    ENDDO
2200                    DO  l = 0, 3
2201                       DO  m = 1, surf_usm_v(l)%ns
2202                          surf_usm_v(l)%t_surf_10cm_av(m) =                         &
2203                                          surf_usm_v(l)%t_surf_10cm_av(m) /         &
2204                                          REAL( average_count_3d, kind=wp )
2205                       ENDDO
2206                    ENDDO
2207                   
2208                CASE ( 'usm_t_wall' )
2209!--                 wall temperature for  iwl layer of walls and land
2210                    DO  m = 1, surf_usm_h%ns
2211                       surf_usm_h%t_wall_av(iwl,m) =                           &
2212                                          surf_usm_h%t_wall_av(iwl,m) /        &
2213                                          REAL( average_count_3d, kind=wp )
2214                    ENDDO
2215                    DO  l = 0, 3
2216                       DO  m = 1, surf_usm_v(l)%ns
2217                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
2218                                          surf_usm_v(l)%t_wall_av(iwl,m) /     &
2219                                          REAL( average_count_3d, kind=wp )
2220                       ENDDO
2221                    ENDDO
2222
2223                CASE ( 'usm_t_window' )
2224!--                 window temperature for  iwl layer of walls and land
2225                    DO  m = 1, surf_usm_h%ns
2226                       surf_usm_h%t_window_av(iwl,m) =                           &
2227                                          surf_usm_h%t_window_av(iwl,m) /        &
2228                                          REAL( average_count_3d, kind=wp )
2229                    ENDDO
2230                    DO  l = 0, 3
2231                       DO  m = 1, surf_usm_v(l)%ns
2232                          surf_usm_v(l)%t_window_av(iwl,m) =                     &
2233                                          surf_usm_v(l)%t_window_av(iwl,m) /     &
2234                                          REAL( average_count_3d, kind=wp )
2235                       ENDDO
2236                    ENDDO
2237
2238                CASE ( 'usm_t_green' )
2239!--                 green temperature for  iwl layer of walls and land
2240                    DO  m = 1, surf_usm_h%ns
2241                       surf_usm_h%t_green_av(iwl,m) =                           &
2242                                          surf_usm_h%t_green_av(iwl,m) /        &
2243                                          REAL( average_count_3d, kind=wp )
2244                    ENDDO
2245                    DO  l = 0, 3
2246                       DO  m = 1, surf_usm_v(l)%ns
2247                          surf_usm_v(l)%t_green_av(iwl,m) =                     &
2248                                          surf_usm_v(l)%t_green_av(iwl,m) /     &
2249                                          REAL( average_count_3d, kind=wp )
2250                       ENDDO
2251                    ENDDO
2252
2253
2254           END SELECT
2255
2256        ENDIF
2257
2258    END SUBROUTINE usm_average_3d_data
2259
2260
2261
2262!------------------------------------------------------------------------------!
2263! Description:
2264! ------------
2265!> Set internal Neumann boundary condition at outer soil grid points
2266!> for temperature and humidity.
2267!------------------------------------------------------------------------------!
2268 SUBROUTINE usm_boundary_condition
2269 
2270    IMPLICIT NONE
2271
2272    INTEGER(iwp) :: i      !< grid index x-direction
2273    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
2274    INTEGER(iwp) :: j      !< grid index y-direction
2275    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
2276    INTEGER(iwp) :: k      !< grid index z-direction
2277    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
2278    INTEGER(iwp) :: l      !< running index surface-orientation
2279    INTEGER(iwp) :: m      !< running index surface elements
2280
2281    koff = surf_usm_h%koff
2282    DO  m = 1, surf_usm_h%ns
2283       i = surf_usm_h%i(m)
2284       j = surf_usm_h%j(m)
2285       k = surf_usm_h%k(m)
2286       pt(k+koff,j,i) = pt(k,j,i)
2287    ENDDO
2288
2289    DO  l = 0, 3
2290       ioff = surf_usm_v(l)%ioff
2291       joff = surf_usm_v(l)%joff
2292       DO  m = 1, surf_usm_v(l)%ns
2293          i = surf_usm_v(l)%i(m)
2294          j = surf_usm_v(l)%j(m)
2295          k = surf_usm_v(l)%k(m)
2296          pt(k,j+joff,i+ioff) = pt(k,j,i)
2297       ENDDO
2298    ENDDO
2299
2300 END SUBROUTINE usm_boundary_condition
2301
2302
2303!------------------------------------------------------------------------------!
2304!
2305! Description:
2306! ------------
2307!> Subroutine checks variables and assigns units.
2308!> It is called out from subroutine check_parameters.
2309!------------------------------------------------------------------------------!
2310    SUBROUTINE usm_check_data_output( variable, unit )
2311       
2312        IMPLICIT NONE
2313 
2314        CHARACTER (len=*),INTENT(IN)    ::  variable !:
2315        CHARACTER (len=*),INTENT(OUT)   ::  unit     !:
2316       
2317        CHARACTER (len=varnamelength)   :: var
2318
2319        var = TRIM(variable)
2320        IF ( var(1:12) == 'usm_rad_net_'  .OR.  var(1:13) == 'usm_rad_insw_'  .OR.        &
2321             var(1:13) == 'usm_rad_inlw_'  .OR.  var(1:16) == 'usm_rad_inswdir_'  .OR.    &
2322             var(1:16) == 'usm_rad_inswdif_'  .OR.  var(1:16) == 'usm_rad_inswref_'  .OR. &
2323             var(1:16) == 'usm_rad_inlwdif_'  .OR.  var(1:16) == 'usm_rad_inlwref_'  .OR. &
2324             var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.     &
2325             var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.     &
2326             var(1:11) == 'usm_rad_hf_'  .OR.                                             &
2327             var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
2328             var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
2329             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' )  THEN
2330            unit = 'W/m2'
[2920]2331        ELSE IF ( var(1:10) == 'usm_t_surf'   .OR.  var(1:10) == 'usm_t_wall' .OR.         &
[2737]2332                  var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR.  &
2333                  var(1:16) == 'usm_t_surf_green'  .OR.                                   &
2334                  var(1:11) == 'usm_t_green' .OR.                                         &
2335                  var(1:15) == 'usm_t_surf_10cm')  THEN
2336            unit = 'K'
2337        ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              & 
2338                  var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
[2920]2339                  var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.     &
2340                  var(1:9) == 'usm_skyvf' .OR. var(1:9) == 'usm_skyvft' )  THEN
[2737]2341            unit = '1'
2342        ELSE
2343            unit = 'illegal'
2344        ENDIF
2345
2346    END SUBROUTINE usm_check_data_output
2347
2348
2349!------------------------------------------------------------------------------!
2350! Description:
2351! ------------
2352!> Check parameters routine for urban surface model
2353!------------------------------------------------------------------------------!
2354    SUBROUTINE usm_check_parameters
2355   
2356       USE control_parameters,                                                 &
2357           ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, large_scale_forcing,   &
2358                  lsf_surf, topography
2359
2360!
2361!--    Dirichlet boundary conditions are required as the surface fluxes are
2362!--    calculated from the temperature/humidity gradients in the urban surface
2363!--    model
2364       IF ( bc_pt_b == 'neumann'   .OR.   bc_q_b == 'neumann' )  THEN
2365          message_string = 'urban surface model requires setting of '//        &
2366                           'bc_pt_b = "dirichlet" and '//                      &
2367                           'bc_q_b  = "dirichlet"'
[3045]2368          CALL message( 'usm_check_parameters', 'PA0590', 1, 2, 0, 6, 0 )
[2737]2369       ENDIF
2370
2371       IF ( .NOT.  constant_flux_layer )  THEN
2372          message_string = 'urban surface model requires '//                   &
2373                           'constant_flux_layer = .T.'
[3045]2374          CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
[2737]2375       ENDIF
[3045]2376
2377       IF (  .NOT.  radiation )  THEN
2378          message_string = 'urban surface model requires '//                   &
2379                           'the radiation model to be switched on'
2380          CALL message( 'usm_check_parameters', 'PA0084', 1, 2, 0, 6, 0 )
2381       ENDIF
[2737]2382!       
2383!--    Surface forcing has to be disabled for LSF in case of enabled
2384!--    urban surface module
2385       IF ( large_scale_forcing )  THEN
2386          lsf_surf = .FALSE.
2387       ENDIF
2388!
2389!--    Topography
2390       IF ( topography == 'flat' )  THEN
2391          message_string = 'topography /= "flat" is required '//               &
2392                           'when using the urban surface model'
2393          CALL message( 'check_parameters', 'PA0592', 1, 2, 0, 6, 0 )
2394       ENDIF
[2920]2395!
2396!--    naheatlayers
2397       IF ( naheatlayers > nzt )  THEN
2398          message_string = 'number of anthropogenic heat layers '//            &
2399                           '"naheatlayers" can not be larger than'//           &
2400                           ' number of domain layers "nzt"'
2401          CALL message( 'check_parameters', 'PA0593', 1, 2, 0, 6, 0 )
2402       ENDIF
[2737]2403
2404    END SUBROUTINE usm_check_parameters
2405
2406
2407!------------------------------------------------------------------------------!
2408!
2409! Description:
2410! ------------
2411!> Output of the 3D-arrays in netCDF and/or AVS format
2412!> for variables of urban_surface model.
2413!> It resorts the urban surface module output quantities from surf style
2414!> indexing into temporary 3D array with indices (i,j,k).
2415!> It is called from subroutine data_output_3d.
2416!------------------------------------------------------------------------------!
2417    SUBROUTINE usm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
2418       
2419        IMPLICIT NONE
2420
2421        INTEGER(iwp), INTENT(IN)       ::  av        !<
2422        CHARACTER (len=*), INTENT(IN)  ::  variable  !<
2423        INTEGER(iwp), INTENT(IN)       ::  nzb_do    !< lower limit of the data output (usually 0)
2424        INTEGER(iwp), INTENT(IN)       ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
2425        LOGICAL, INTENT(OUT)           ::  found     !<
2426        REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< sp - it has to correspond to module data_output_3d
2427        REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr)     ::  temp_pf    !< temp array for urban surface output procedure
2428       
2429        CHARACTER (len=varnamelength)                          :: var, surfid
2430        INTEGER(iwp), PARAMETER                                :: nd = 5
2431        CHARACTER(len=6), DIMENSION(0:nd-1), PARAMETER         :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
[2920]2432        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: dirint =  (/    iup_u, isouth_u, inorth_u,  iwest_u,  ieast_u /)
2433        INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER             :: diridx =  (/       -1,        1,        0,        3,        2 /)
2434                                                                     !< index for surf_*_v: 0:3 = (North, South, East, West)
[2737]2435        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
2436        INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
[2920]2437        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt
[2737]2438        INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
[3014]2439        INTEGER(iwp)                                           :: k_topo    !< topography top index
[2737]2440
2441        dirstart = (/ startland, startwall, startwall, startwall, startwall /)
2442        dirend = (/ endland, endwall, endwall, endwall, endwall /)
2443
2444        found = .TRUE.
2445        temp_pf = -1._wp
2446       
2447        ids = -1
2448        var = TRIM(variable)
2449        DO i = 0, nd-1
2450            k = len(TRIM(var))
2451            j = len(TRIM(dirname(i)))
2452            IF ( var(k-j+1:k) == dirname(i) )  THEN
2453                ids = i
[2920]2454                idsint = dirint(ids)
2455                idsidx = diridx(ids)
[2737]2456                var = var(:k-j)
2457                EXIT
2458            ENDIF
2459        ENDDO
2460        IF ( ids == -1 )  THEN
2461            var = TRIM(variable)
2462        ENDIF
2463        IF ( var(1:11) == 'usm_t_wall_'  .AND.  len(TRIM(var)) >= 12 )  THEN
2464!--         wall layers
2465            READ(var(12:12), '(I1)', iostat=istat ) iwl
2466            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
2467                var = var(1:10)
2468            ENDIF
2469        ENDIF
2470        IF ( var(1:13) == 'usm_t_window_'  .AND.  len(TRIM(var)) >= 14 )  THEN
2471!--         window layers
2472            READ(var(14:14), '(I1)', iostat=istat ) iwl
2473            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
2474                var = var(1:12)
2475            ENDIF
2476        ENDIF
2477        IF ( var(1:12) == 'usm_t_green_'  .AND.  len(TRIM(var)) >= 13 )  THEN
2478!--         green layers
2479            READ(var(13:13), '(I1)', iostat=istat ) iwl
2480            IF ( istat == 0  .AND.  iwl >= nzb_wall  .AND.  iwl <= nzt_wall )  THEN
2481                var = var(1:11)
2482            ENDIF
2483        ENDIF
2484        IF ( (var(1:8) == 'usm_svf_'  .OR.  var(1:8) == 'usm_dif_')  .AND.  len(TRIM(var)) >= 13 )  THEN
2485!--         svf values to particular surface
2486            surfid = var(9:)
2487            i = index(surfid,'_')
2488            j = index(surfid(i+1:),'_')
2489            READ(surfid(1:i-1),*, iostat=istat ) is
2490            IF ( istat == 0 )  THEN
2491                READ(surfid(i+1:i+j-1),*, iostat=istat ) js
2492            ENDIF
2493            IF ( istat == 0 )  THEN
2494                READ(surfid(i+j+1:),*, iostat=istat ) ks
2495            ENDIF
2496            IF ( istat == 0 )  THEN
2497                var = var(1:7)
2498            ENDIF
2499        ENDIF
2500       
2501        SELECT CASE ( TRIM(var) )
2502
2503          CASE ( 'usm_surfz' )
2504!--           array of lw radiation falling to local surface after i-th reflection
[2920]2505              IF ( idsint == iup_u )  THEN
2506                 DO  m = 1, surf_usm_h%ns
2507                    i = surf_usm_h%i(m)
2508                    j = surf_usm_h%j(m)
2509                    k = surf_usm_h%k(m)
2510                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) )
2511                 ENDDO
2512              ELSE
2513                 l = idsidx
[2737]2514                 DO  m = 1, surf_usm_v(l)%ns
2515                    i = surf_usm_v(l)%i(m)
2516                    j = surf_usm_v(l)%j(m)
2517                    k = surf_usm_v(l)%k(m)
2518                    temp_pf(0,j,i) = MAX( temp_pf(0,j,i), REAL( k, kind=wp) + 1.0_wp )
2519                 ENDDO
[2920]2520              ENDIF
[2737]2521
2522          CASE ( 'usm_surfcat' )
2523!--           surface category
[2920]2524              IF ( idsint == iup_u )  THEN
2525                 DO  m = 1, surf_usm_h%ns
2526                    i = surf_usm_h%i(m)
2527                    j = surf_usm_h%j(m)
2528                    k = surf_usm_h%k(m)
2529                    temp_pf(k,j,i) = surf_usm_h%surface_types(m)
2530                 ENDDO
2531              ELSE
2532                 l = idsidx
[2737]2533                 DO  m = 1, surf_usm_v(l)%ns
2534                    i = surf_usm_v(l)%i(m)
2535                    j = surf_usm_v(l)%j(m)
2536                    k = surf_usm_v(l)%k(m)
2537                    temp_pf(k,j,i) = surf_usm_v(l)%surface_types(m)
2538                 ENDDO
[2920]2539              ENDIF
[2737]2540             
2541          CASE ( 'usm_surfalb' )
2542!--           surface albedo, weighted average
[2920]2543              IF ( idsint == iup_u )  THEN
2544                 DO  m = 1, surf_usm_h%ns
2545                    i = surf_usm_h%i(m)
2546                    j = surf_usm_h%j(m)
2547                    k = surf_usm_h%k(m)
[2963]2548                    temp_pf(k,j,i) = surf_usm_h%frac(ind_veg_wall,m)     *     &
2549                                     surf_usm_h%albedo(ind_veg_wall,m)  +      &
2550                                     surf_usm_h%frac(ind_pav_green,m)    *     &
2551                                     surf_usm_h%albedo(ind_pav_green,m) +      &
2552                                     surf_usm_h%frac(ind_wat_win,m)      *     &
2553                                     surf_usm_h%albedo(ind_wat_win,m)
[2920]2554                 ENDDO
2555              ELSE
2556                 l = idsidx
[2737]2557                 DO  m = 1, surf_usm_v(l)%ns
2558                    i = surf_usm_v(l)%i(m)
2559                    j = surf_usm_v(l)%j(m)
2560                    k = surf_usm_v(l)%k(m)
[2963]2561                    temp_pf(k,j,i) = surf_usm_v(l)%frac(ind_veg_wall,m)     *  &
2562                                     surf_usm_v(l)%albedo(ind_veg_wall,m)  +   &
2563                                     surf_usm_v(l)%frac(ind_pav_green,m)    *  &
2564                                     surf_usm_v(l)%albedo(ind_pav_green,m) +   &
2565                                     surf_usm_v(l)%frac(ind_wat_win,m)      *  &
2566                                     surf_usm_v(l)%albedo(ind_wat_win,m)
[2737]2567                 ENDDO
[2920]2568              ENDIF
[2737]2569             
2570          CASE ( 'usm_surfemis' )
2571!--           surface emissivity, weighted average
[2920]2572              IF ( idsint == iup_u )  THEN
2573                 DO  m = 1, surf_usm_h%ns
2574                    i = surf_usm_h%i(m)
2575                    j = surf_usm_h%j(m)
2576                    k = surf_usm_h%k(m)
[2963]2577                    temp_pf(k,j,i) =  surf_usm_h%frac(ind_veg_wall,m)      *   &
2578                                      surf_usm_h%emissivity(ind_veg_wall,m)  + &
2579                                      surf_usm_h%frac(ind_pav_green,m)     *   &
2580                                      surf_usm_h%emissivity(ind_pav_green,m) + &
2581                                      surf_usm_h%frac(ind_wat_win,m)       *   &
2582                                      surf_usm_h%emissivity(ind_wat_win,m)
[2920]2583                 ENDDO
2584              ELSE
2585                 l = idsidx
[2737]2586                 DO  m = 1, surf_usm_v(l)%ns
2587                    i = surf_usm_v(l)%i(m)
2588                    j = surf_usm_v(l)%j(m)
2589                    k = surf_usm_v(l)%k(m)
[2963]2590                    temp_pf(k,j,i) = surf_usm_v(l)%frac(ind_veg_wall,m)       *&
2591                                     surf_usm_v(l)%emissivity(ind_veg_wall,m) +&
2592                                     surf_usm_v(l)%frac(ind_pav_green,m)      *&
2593                                     surf_usm_v(l)%emissivity(ind_pav_green,m)+&
2594                                     surf_usm_v(l)%frac(ind_wat_win,m)        *&
2595                                     surf_usm_v(l)%emissivity(ind_wat_win,m)
[2737]2596                 ENDDO
[2920]2597              ENDIF
[2737]2598
2599          CASE ( 'usm_surfwintrans' )
2600!--           transmissivity window tiles
[2920]2601              IF ( idsint == iup_u )  THEN
2602                 DO  m = 1, surf_usm_h%ns
2603                    i = surf_usm_h%i(m)
2604                    j = surf_usm_h%j(m)
2605                    k = surf_usm_h%k(m)
2606                    temp_pf(k,j,i) = surf_usm_h%transmissivity(m)
2607                 ENDDO
2608              ELSE
2609                 l = idsidx
[2737]2610                 DO  m = 1, surf_usm_v(l)%ns
2611                    i = surf_usm_v(l)%i(m)
2612                    j = surf_usm_v(l)%j(m)
2613                    k = surf_usm_v(l)%k(m)
2614                    temp_pf(k,j,i) = surf_usm_v(l)%transmissivity(m)
2615                 ENDDO
[2920]2616              ENDIF
[2737]2617
[2920]2618          CASE ( 'usm_skyvf' )
2619!--           sky view factor
2620              DO isurf = dirstart(ids), dirend(ids)
2621                 IF ( surfl(id,isurf) == idsint )  THEN
2622                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvf(isurf)
2623                 ENDIF
[2737]2624              ENDDO
[2920]2625             
2626          CASE ( 'usm_skyvft' )
2627!--           sky view factor
2628              DO isurf = dirstart(ids), dirend(ids)
2629                 IF ( surfl(id,isurf) == ids )  THEN
2630                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = skyvft(isurf)
2631                 ENDIF
2632              ENDDO
[2737]2633
2634!
2635!-- Not adjusted so far             
2636          CASE ( 'usm_svf', 'usm_dif' )
2637!--           shape view factors or iradiance factors to selected surface
2638              IF ( TRIM(var)=='usm_svf' )  THEN
2639                  k = 1
2640              ELSE
2641                  k = 2
2642              ENDIF
2643              DO isvf = 1, nsvfl
2644                  isurflt = svfsurf(1, isvf)
2645                  isurfs = svfsurf(2, isvf)
2646                             
2647                  IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND.       &
[2920]2648                       surf(iz,isurfs) == ks  .AND.  surf(id,isurfs) == idsint )  THEN
[2737]2649  !--                 correct source surface
2650                      temp_pf(surfl(iz,isurflt),surfl(iy,isurflt),surfl(ix,isurflt)) = svf(k,isvf)
2651                  ENDIF
2652              ENDDO
2653
2654          CASE ( 'usm_rad_net' )
2655!--           array of complete radiation balance
2656              IF ( av == 0 )  THEN
[2920]2657                 IF ( idsint == iup_u )  THEN
2658                    DO  m = 1, surf_usm_h%ns
2659                       i = surf_usm_h%i(m)
2660                       j = surf_usm_h%j(m)
2661                       k = surf_usm_h%k(m)
2662                       temp_pf(k,j,i) = surf_usm_h%rad_net_l(m)
2663                    ENDDO
2664                 ELSE
2665                    l = idsidx
[2737]2666                    DO  m = 1, surf_usm_v(l)%ns
2667                       i = surf_usm_v(l)%i(m)
2668                       j = surf_usm_v(l)%j(m)
2669                       k = surf_usm_v(l)%k(m)
2670                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_l(m)
2671                    ENDDO
[2920]2672                 ENDIF
[2737]2673              ELSE
[2920]2674                 IF ( idsint == iup_u )  THEN
2675                    DO  m = 1, surf_usm_h%ns
2676                       i = surf_usm_h%i(m)
2677                       j = surf_usm_h%j(m)
2678                       k = surf_usm_h%k(m)
2679                       temp_pf(k,j,i) = surf_usm_h%rad_net_av(m)
2680                    ENDDO
2681                 ELSE
2682                    l = idsidx
[2737]2683                    DO  m = 1, surf_usm_v(l)%ns
2684                       i = surf_usm_v(l)%i(m)
2685                       j = surf_usm_v(l)%j(m)
2686                       k = surf_usm_v(l)%k(m)
2687                       temp_pf(k,j,i) = surf_usm_v(l)%rad_net_av(m)
2688                    ENDDO
[2920]2689                 ENDIF
[2737]2690              ENDIF
2691
2692          CASE ( 'usm_rad_insw' )
2693!--           array of sw radiation falling to surface after i-th reflection
2694              DO isurf = dirstart(ids), dirend(ids)
[2920]2695                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2696                   IF ( av == 0 )  THEN
2697                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw(isurf)
2698                   ELSE
2699                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinsw_av(isurf)
2700                   ENDIF
2701                 ENDIF
2702              ENDDO
2703
2704          CASE ( 'usm_rad_inlw' )
2705!--           array of lw radiation falling to surface after i-th reflection
2706              DO isurf = dirstart(ids), dirend(ids)
[2920]2707                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2708                   IF ( av == 0 )  THEN
2709                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf)
2710                   ELSE
2711                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw_av(isurf)
2712                   ENDIF
2713                 ENDIF
2714              ENDDO
2715
2716          CASE ( 'usm_rad_inswdir' )
2717!--           array of direct sw radiation falling to surface from sun
2718              DO isurf = dirstart(ids), dirend(ids)
[2920]2719                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2720                   IF ( av == 0 )  THEN
2721                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir(isurf)
2722                   ELSE
2723                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdir_av(isurf)
2724                   ENDIF
2725                 ENDIF
2726              ENDDO
2727
2728          CASE ( 'usm_rad_inswdif' )
2729!--           array of difusion sw radiation falling to surface from sky and borders of the domain
2730              DO isurf = dirstart(ids), dirend(ids)
[2920]2731                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2732                   IF ( av == 0 )  THEN
2733                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif(isurf)
2734                   ELSE
2735                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswdif_av(isurf)
2736                   ENDIF
2737                 ENDIF
2738              ENDDO
2739
2740          CASE ( 'usm_rad_inswref' )
2741!--           array of sw radiation falling to surface from reflections
2742              DO isurf = dirstart(ids), dirend(ids)
[2920]2743                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2744                   IF ( av == 0 )  THEN
2745                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = &
2746                       surfinsw(isurf) - surfinswdir(isurf) - surfinswdif(isurf)
2747                   ELSE
2748                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinswref_av(isurf)
2749                   ENDIF
2750                 ENDIF
2751              ENDDO
2752
[2920]2753          CASE ( 'usm_rad_inlwdif' )
2754!--           array of difusion lw radiation falling to surface from sky and borders of the domain
2755              DO isurf = dirstart(ids), dirend(ids)
2756                 IF ( surfl(id,isurf) == idsint )  THEN
2757                   IF ( av == 0 )  THEN
2758                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif(isurf)
2759                   ELSE
2760                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwdif_av(isurf)
2761                   ENDIF
2762                 ENDIF
2763              ENDDO
2764
[2737]2765          CASE ( 'usm_rad_inlwref' )
2766!--           array of lw radiation falling to surface from reflections
2767              DO isurf = dirstart(ids), dirend(ids)
[2920]2768                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2769                   IF ( av == 0 )  THEN
2770                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
2771                   ELSE
2772                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinlwref_av(isurf)
2773                   ENDIF
2774                 ENDIF
2775              ENDDO
2776
2777          CASE ( 'usm_rad_outsw' )
2778!--           array of sw radiation emitted from surface after i-th reflection
2779              DO isurf = dirstart(ids), dirend(ids)
[2920]2780                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2781                   IF ( av == 0 )  THEN
2782                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw(isurf)
2783                   ELSE
2784                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutsw_av(isurf)
2785                   ENDIF
2786                 ENDIF
2787              ENDDO
2788
2789          CASE ( 'usm_rad_outlw' )
2790!--           array of lw radiation emitted from surface after i-th reflection
2791              DO isurf = dirstart(ids), dirend(ids)
[2920]2792                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2793                   IF ( av == 0 )  THEN
2794                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw(isurf)
2795                   ELSE
2796                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfoutlw_av(isurf)
2797                   ENDIF
2798                 ENDIF
2799              ENDDO
2800
2801          CASE ( 'usm_rad_ressw' )
2802!--           average of array of residua of sw radiation absorbed in surface after last reflection
2803              DO isurf = dirstart(ids), dirend(ids)
[2920]2804                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2805                   IF ( av == 0 )  THEN
2806                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins(isurf)
2807                   ELSE
2808                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfins_av(isurf)
2809                   ENDIF
2810                 ENDIF
2811              ENDDO
2812
2813          CASE ( 'usm_rad_reslw' )
2814!--           average of array of residua of lw radiation absorbed in surface after last reflection
2815              DO isurf = dirstart(ids), dirend(ids)
[2920]2816                 IF ( surfl(id,isurf) == idsint )  THEN
[2737]2817                   IF ( av == 0 )  THEN
2818                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl(isurf)
2819                   ELSE
2820                     temp_pf(surfl(iz,isurf),surfl(iy,isurf),surfl(ix,isurf)) = surfinl_av(isurf)
2821                   ENDIF
2822                 ENDIF
2823              ENDDO
2824 
2825          CASE ( 'usm_rad_hf' )
2826!--           array of heat flux from radiation for surfaces after all reflections
2827              IF ( av == 0 )  THEN
[2920]2828                 IF ( idsint == iup_u )  THEN
2829                    DO  m = 1, surf_usm_h%ns
2830                       i = surf_usm_h%i(m)
2831                       j = surf_usm_h%j(m)
2832                       k = surf_usm_h%k(m)
2833                       temp_pf(k,j,i) = surf_usm_h%surfhf(m)
2834                    ENDDO
2835                 ELSE
2836                    l = idsidx
[2737]2837                    DO  m = 1, surf_usm_v(l)%ns
2838                       i = surf_usm_v(l)%i(m)
2839                       j = surf_usm_v(l)%j(m)
2840                       k = surf_usm_v(l)%k(m)
2841                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf(m)
2842                    ENDDO
[2920]2843                 ENDIF
[2737]2844              ELSE
[2920]2845                 IF ( idsint == iup_u )  THEN
2846                    DO  m = 1, surf_usm_h%ns
2847                       i = surf_usm_h%i(m)
2848                       j = surf_usm_h%j(m)
2849                       k = surf_usm_h%k(m)
2850                       temp_pf(k,j,i) = surf_usm_h%surfhf_av(m)
2851                    ENDDO
2852                 ELSE
2853                    l = idsidx
[2737]2854                    DO  m = 1, surf_usm_v(l)%ns
2855                       i = surf_usm_v(l)%i(m)
2856                       j = surf_usm_v(l)%j(m)
2857                       k = surf_usm_v(l)%k(m)
2858                       temp_pf(k,j,i) = surf_usm_v(l)%surfhf_av(m)
2859                    ENDDO
[2920]2860                 ENDIF
[2737]2861              ENDIF
2862 
2863          CASE ( 'usm_wshf' )
2864!--           array of sensible heat flux from surfaces
2865              IF ( av == 0 )  THEN
[2920]2866                 IF ( idsint == iup_u )  THEN
2867                    DO  m = 1, surf_usm_h%ns
2868                       i = surf_usm_h%i(m)
2869                       j = surf_usm_h%j(m)
2870                       k = surf_usm_h%k(m)
2871                       temp_pf(k,j,i) = surf_usm_h%wshf_eb(m)
2872                    ENDDO
2873                 ELSE
2874                    l = idsidx
[2737]2875                    DO  m = 1, surf_usm_v(l)%ns
2876                       i = surf_usm_v(l)%i(m)
2877                       j = surf_usm_v(l)%j(m)
2878                       k = surf_usm_v(l)%k(m)
2879                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb(m)
2880                    ENDDO
[2920]2881                 ENDIF
[2737]2882              ELSE
[2920]2883                 IF ( idsint == iup_u )  THEN
2884                    DO  m = 1, surf_usm_h%ns
2885                       i = surf_usm_h%i(m)
2886                       j = surf_usm_h%j(m)
2887                       k = surf_usm_h%k(m)
2888                       temp_pf(k,j,i) = surf_usm_h%wshf_eb_av(m)
2889                    ENDDO
2890                 ELSE
2891                    l = idsidx
[2737]2892                    DO  m = 1, surf_usm_v(l)%ns
2893                       i = surf_usm_v(l)%i(m)
2894                       j = surf_usm_v(l)%j(m)
2895                       k = surf_usm_v(l)%k(m)
2896                       temp_pf(k,j,i) = surf_usm_v(l)%wshf_eb_av(m)
2897                    ENDDO
[2920]2898                 ENDIF
[2737]2899              ENDIF
2900
2901
2902          CASE ( 'usm_wghf' )
2903!--           array of heat flux from ground (land, wall, roof)
2904              IF ( av == 0 )  THEN
[2920]2905                 IF ( idsint == iup_u )  THEN
2906                    DO  m = 1, surf_usm_h%ns
2907                       i = surf_usm_h%i(m)
2908                       j = surf_usm_h%j(m)
2909                       k = surf_usm_h%k(m)
2910                       temp_pf(k,j,i) = surf_usm_h%wghf_eb(m)
2911                    ENDDO
2912                 ELSE
2913                    l = idsidx
[2737]2914                    DO  m = 1, surf_usm_v(l)%ns
2915                       i = surf_usm_v(l)%i(m)
2916                       j = surf_usm_v(l)%j(m)
2917                       k = surf_usm_v(l)%k(m)
2918                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb(m)
2919                    ENDDO
[2920]2920                 ENDIF
[2737]2921              ELSE
[2920]2922                 IF ( idsint == iup_u )  THEN
2923                    DO  m = 1, surf_usm_h%ns
2924                       i = surf_usm_h%i(m)
2925                       j = surf_usm_h%j(m)
2926                       k = surf_usm_h%k(m)
2927                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_av(m)
2928                    ENDDO
2929                 ELSE
2930                    l = idsidx
[2737]2931                    DO  m = 1, surf_usm_v(l)%ns
2932                       i = surf_usm_v(l)%i(m)
2933                       j = surf_usm_v(l)%j(m)
2934                       k = surf_usm_v(l)%k(m)
2935                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_av(m)
2936                    ENDDO
[2920]2937                 ENDIF
[2737]2938              ENDIF
2939
2940          CASE ( 'usm_wghf_window' )
2941!--           array of heat flux from window ground (land, wall, roof)
2942
2943              IF ( av == 0 )  THEN
[2920]2944                 IF ( idsint == iup_u )  THEN
2945                    DO  m = 1, surf_usm_h%ns
2946                       i = surf_usm_h%i(m)
2947                       j = surf_usm_h%j(m)
2948                       k = surf_usm_h%k(m)
2949                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_window(m)
2950                    ENDDO
2951                 ELSE
2952                    l = idsidx
[2737]2953                    DO  m = 1, surf_usm_v(l)%ns
2954                       i = surf_usm_v(l)%i(m)
2955                       j = surf_usm_v(l)%j(m)
2956                       k = surf_usm_v(l)%k(m)
2957                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window(m)
2958                    ENDDO
[2920]2959                 ENDIF
[2737]2960              ELSE
[2920]2961                 IF ( idsint == iup_u )  THEN
2962                    DO  m = 1, surf_usm_h%ns
2963                       i = surf_usm_h%i(m)
2964                       j = surf_usm_h%j(m)
2965                       k = surf_usm_h%k(m)
2966                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_window_av(m)
2967                    ENDDO
2968                 ELSE
2969                    l = idsidx
[2737]2970                    DO  m = 1, surf_usm_v(l)%ns
2971                       i = surf_usm_v(l)%i(m)
2972                       j = surf_usm_v(l)%j(m)
2973                       k = surf_usm_v(l)%k(m)
2974                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_window_av(m)
2975                    ENDDO
[2920]2976                 ENDIF
[2737]2977              ENDIF
2978
2979          CASE ( 'usm_wghf_green' )
2980!--           array of heat flux from green ground (land, wall, roof)
2981
2982              IF ( av == 0 )  THEN
[2920]2983                 IF ( idsint == iup_u )  THEN
2984                    DO  m = 1, surf_usm_h%ns
2985                       i = surf_usm_h%i(m)
2986                       j = surf_usm_h%j(m)
2987                       k = surf_usm_h%k(m)
2988                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_green(m)
2989                    ENDDO
2990                 ELSE
2991                    l = idsidx
[2737]2992                    DO  m = 1, surf_usm_v(l)%ns
2993                       i = surf_usm_v(l)%i(m)
2994                       j = surf_usm_v(l)%j(m)
2995                       k = surf_usm_v(l)%k(m)
2996                       temp_pf(k,j,i) = surf_usm_v(l)%wghf_eb_green(m)
2997                    ENDDO
[2920]2998                 ENDIF
[2737]2999              ELSE
[2920]3000                 IF ( idsint == iup_u )  THEN
3001                    DO  m = 1, surf_usm_h%ns
3002                       i = surf_usm_h%i(m)
3003                       j = surf_usm_h%j(m)
3004                       k = surf_usm_h%k(m)
3005                       temp_pf(k,j,i) = surf_usm_h%wghf_eb_green_av(m)
3006                    ENDDO
3007                 ELSE
3008                    l = idsidx
[2737]3009                    DO  m = 1, surf_usm_v(l)%ns
3010                       i = surf_usm_v(l)%i(m)
3011                       j = surf_usm_v(l)%j(m)