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

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

Fixes for radiative transfer model

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