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

Last change on this file since 3987 was 3987, checked in by kanani, 2 years ago

clean up location, debug and error messages

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