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

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

renamed all Fortran NAMELISTS

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