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

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

Minor revision of static input file checks, bugfix in initialization of surface-fractions in LSM; minor bugfix in initialization of albedo at window-surfaces; for clearer access of albedo and emissivity introduce index for vegetation/wall, pavement/green-wall and water/window surfaces

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