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

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

last commit documented

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