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

Last change on this file since 3993 was 3993, checked in by suehring, 4 years ago

Rename the USM-internal flag spinup into during_spinup, in order to avoid confusion with global control parameter

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