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

Last change on this file since 3614 was 3614, checked in by raasch, 3 years ago

unused variables removed, abort renamed inifor_abort to avoid intrinsic problem in Fortran

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