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

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

Reading/Writing? data in case of restart runs revised

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