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

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

New vertical stretching procedure has been introduced

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