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

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

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

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