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

Last change on this file since 4442 was 4442, checked in by suehring, 4 years ago

last commit documented

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