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

Last change on this file since 2943 was 2943, checked in by suehring, 3 years ago

Calculate exner function at all height levels in USM In order to avoid segmentation faults in case radiation-interactions are switched-off, calculate exner function at all height levels in USM .

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