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

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

Output of ground-heat flux at natural- and urban-type surfaces in one output variable; enable restart data of _av variables that belong to both land- and urban-surface model

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