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

Last change on this file since 3241 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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