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

Last change on this file since 2737 was 2737, checked in by kanani, 3 years ago

Removed unused variables t_surf_whole...

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