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

Last change on this file since 3246 was 3246, checked in by sward, 6 years ago

Added error handling for wrong input parameters

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