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

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

Remove work-around for green surface fraction on buildings - do not set it zero

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