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

Last change on this file since 3704 was 3704, checked in by suehring, 2 years ago

Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.

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