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

Last change on this file since 2705 was 2705, checked in by maronga, 4 years ago

minor bugfix for restarts

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