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

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

further inipar parameter has been added to restart data, bugfix in spinup mechanism

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