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

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

new procedure for reading/writing svf data, initialization of local variable ids

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