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

Last change on this file since 3029 was 3029, checked in by raasch, 3 years ago

bugfix: close unit 151 instead of 90

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