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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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