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

Last change on this file since 3014 was 3014, checked in by maronga, 3 years ago

series of bugfixes

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