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

Last change on this file since 3046 was 3046, checked in by Giersch, 3 years ago

Remaining error messages revised, comments extended

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