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

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

Enable 3D data output also with 64-bit precision

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