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

Last change on this file since 2720 was 2720, checked in by kanani, 4 years ago

update of palm version

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