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

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

Set green-fraction on buildings to zero for now, in order to avoid crashes in the green-heat model

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