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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

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