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

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

Removal of chem directive, plus minor changes

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