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

Last change on this file since 3152 was 3152, checked in by suehring, 3 years ago

Further adjustments for surface structure

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