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

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

Limit roughness for heat and moisture in case it exceeds the model-assumed surface-layer height; Limit surface resistance to positive values.

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