source: palm/trunk/SOURCE/urban_surface_mod.f90

Last change on this file was 4893, checked in by raasch, 4 years ago

revised output of surface data via MPI-IO for better performance

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