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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

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