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

Last change on this file since 2920 was 2920, checked in by kanani, 3 years ago

Optimize SVF calculation, clean-up, bugfixes

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