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

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

Bugfix in re-mapping surface-element data in case of restarts; bugfix in initialization of water-surface roughness; bugfix - uninitialized resistance at urban-type surfaces

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