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

Last change on this file since 2765 was 2765, checked in by maronga, 3 years ago

major bugfix in calculation of aerodynamic resistance for vertical surface elements

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