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

Last change on this file since 2735 was 2735, checked in by suehring, 4 years ago

Output of resistance also urban-type surfaces

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