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

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

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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