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

Last change on this file since 4360 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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