Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 3743)
+++ /palm/trunk/SOURCE/Makefile	(revision 3744)
@@ -20,5 +20,6 @@
 # Current revisions:
 # ------------------
-# 
+# Add Netcdf data output for indoor model
+# Remove some dependencies
 # 
 # Former revisions:
@@ -33,14 +34,14 @@
 # 3700 2019-01-26 17:03:42Z knoop
 # Rename surface_output_mod into surface_data_output_mod
-#
+# 
 # 3637 2018-12-20 01:51:36Z knoop
 # Implementation of the PALM module interface
-#
+# 
 # 3634 2018-12-18 12:31:28Z knoop
 # OpenACC port for SPEC
-#
+# 
 # 3579 2018-11-29 15:32:39Z suehring
 # Dependency for check_parameters on nesting_offl_mod added
-#
+# 
 # 3569 2018-11-27 17:03:40Z kanani
 # dom_dwd_user, Schrempf:
@@ -50,22 +51,22 @@
 # 3525 2018-11-14 16:06:14Z kanani
 # Changes related to clean-up of biometeorology (dom_dwd_user)
-#
+# 
 # 3522 2018-11-13 12:14:36Z suehring
 # Dependencies for virtual measurement module added
-#
+# 
 # 3494 2018-11-06 14:51:27Z suehring
 # Surface output revised
-#
+# 
 # 3474 2018-10-30 21:07:39Z kanani
 # Add virtual measurement module
-#
+# 
 # 3472 2018-10-30 20:43:50Z suehring
 # Add indoor model (kanani, srissman, tlang),
 # minor formatting
-#
+# 
 # 3467 2018-10-30 19:05:21Z suehring
 # Implementation of a new aerosol module salsa.
-#
-#
+# 
+# 
 # 3458 2018-10-30 14:51:23Z kanani
 # from chemistry branch r3443, banzhafs, Russo, forkel, basit:
@@ -74,23 +75,23 @@
 # Added chemistry emission module
 # chemistry_model_mod added to flow_statistics
-#
+# 
 # 3448 2018-10-29 18:14:31Z kanani
 # Adjustment of biometeorology dependencies
-#
+# 
 # 3436 2018-10-26 18:35:15Z gronemeier
 # Add surface_mod to user_data_output_mask
-#
+# 
 # 3435 2018-10-26 18:25:44Z gronemeier
 # - Add surface_mod to data_output_mask
 # - Add chemistry_model_mod and surface_mod to init_masks
-#
+# 
 # 3421 2018-10-24 18:39:32Z gronemeier
 # Add netcdf_data_input_mod to netcdf_interface_mod
 # bugfix: add dependencies to chemistry_model_mod
 # Add module for surface data output
-#
+# 
 # 3381 2018-10-19 13:09:06Z raasch
 # dependencies for ocean_mod fixed
-#
+# 
 # 3355 2018-10-16 14:03:34Z knoop
 # Add module for offline nesting;
@@ -98,13 +99,13 @@
 # Bugfix, missing dependency for turbulence generator in init_3d_model;
 # Some formatting ajdustments
-#
+# 
 # 3343 2018-10-15 10:38:52Z suehring
 # (from branch resler)
 # Add biometeorology
-#
-#
+# 
+# 
 # 3322 2018-10-09 10:02:39Z kanani
 # Formatting and cleanup
-#
+# 
 # 3298 2018-10-02 12:21:11Z kanani
 # Added missing dependencies and replaced blanks with tabs (forkel)
@@ -115,68 +116,68 @@
 # changes related to modularization of the ocean mode,
 # bugfix: dependency to advec_ws was missed in chemistry_model_mod
-#
+# 
 # 3274 2018-09-24 15:42:55Z knoop
 # Added palm dependency of multi_agent_system_mod, because of mas_last_actions
 # call at the end of palm run
-#
+# 
 # 3167 2018-07-24 18:17:30Z suehring
 # Bugfix, add missing dependencies for multi-agent system
-#
+# 
 # 3159 2018-07-20 11:20:01Z sward
 # Added multi agent system
-#
+# 
 # 3130 2018-07-16 11:08:55Z gronemeier
 # add surface_layer_fluxes_mod to turbulence_closure_mod
-#
+# 
 # 3129 2018-07-16 07:45:13Z gronemeier
 # add turbulence_closure_mod to parin
-#
+# 
 # 2963 2018-04-12 14:47:44Z suehring
-# Introduce index for vegetation/wall, pavement/green-wall and water/window
+# Introduce index for vegetation/wall, pavement/green-wall and water/window 
 # surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
-#
+# 
 # 2955 2018-04-09 15:14:01Z suehring
 # Add log-points to measure CPU time of NetCDF data input.
-#
+# 
 # 2938 2018-03-27 15:52:42Z suehring
 # No initialization of child domains via dynamic input file, except for soil
 # moisture and temperature
 # Apply turbulence generator at non-cyclic lateral boundary in nesting case
-#
+# 
 # 2936 2018-03-27 14:49:27Z suehring
 # Added dependencies for parent and child synchronization
-#
+# 
 # 2921 2018-03-22 15:05:23Z Giersch
 # date_and_time_mod dependency has been added to read/write_restart_data_mod
-#
+# 
 # 2918 2018-03-21 15:52:14Z gronemeier
-# read/write_3d_binary and read/write_var_list has been removed,
+# read/write_3d_binary and read/write_var_list has been removed, 
 # read/write_restart_data_mod, wrd_write_string and
-# user_read/write_restart_data_mod has been added, dependencies with respect to
+# user_read/write_restart_data_mod has been added, dependencies with respect to 
 # the aforementioned routines have been added/removed
 #
 # 2847 2018-03-02 21:45:58Z suehring
 # Changed format and enforced sorting
-#
+# 
 # 2817 2018-02-19 16:32:21Z knoop
 # Preliminary gust module interface implemented
-#
+# 
 # 2802 2018-02-14 16:21:39Z thiele
 # Changed lpm from subroutine to module.
 # Introduce particle transfer in nested models.
-#
+# 
 # 2773 2018-01-30 14:12:54Z suehring
 # Nesting of chemical species
-#
+# 
 # 2718 2018-01-02 08:49:38Z maronga
 # Corrected "Former revisions" section
-#
+# 
 # 2697 2017-12-14 17:57:20Z kanani
 # Bugfix, missing dependencies
-#
+# 
 # 2696 2017-12-14 17:12:51Z kanani
 # Change in file header (GPL part)
 # Implementation of uv exposure model (FK)
-# Bugfix, removed loop dependcy for vertical_nesting_mod and
+# Bugfix, removed loop dependcy for vertical_nesting_mod and 
 # turbulence_closure_mod, added depencies for vertical_nesting_mod (TG)
 # implemented turbulence_closure_mod (TG)
@@ -186,13 +187,13 @@
 # For LSM, add dependency on calc_mean_profile (??)
 # poismg_noopt modularized and renamed into poismg_noopt_mod
-# add dependencies for netcdf_data_input_mod, calc_mean_profile,
+# add dependencies for netcdf_data_input_mod, calc_mean_profile, 
 # radiation_module_mod, land_surface_model_mod (MS)
-#
+# 
 # 2608 2017-11-13 14:04:26Z schwenkel
-# Added diagnostic_quantities_mod
-#
+# Added diagnostic_quantities_mod 
+# 
 # 2600 2017-11-01 14:11:20Z raasch
 # comment line concerning bound checks removed
-#
+# 
 # 2599 2017-11-01 13:18:45Z hellstea
 # virtual_flight_mod, synthetic_turbulence_generator_mod and
@@ -201,40 +202,40 @@
 # 2563 2017-10-19 15:36:10Z Giersch
 # wind_turbine_model_mod and synthetic_turbulence_generator_mod were added to
-# write_var_list and virtual_flight_mod was deleted from read_var_list
+# write_var_list and virtual_flight_mod was deleted from read_var_list 
 #
 # 2544 2017-10-13 18:09:32Z maronga
-# Added date_and_time_mod
-#
+# Added date_and_time_mod 
+# 
 # 2371 2017-08-24 13:01:17Z kanani
 # Corrected dependencies for vertical_nesting_mod
-#
+# 
 # 2370 2017-08-23 06:11:43Z raasch
 # dependency bugfix for synthetic_turbulence_generator
-#
+# 
 # 2365 2017-08-21 14:59:59Z kanani
 # Added dependencies for vertical_nesting_mod
-#
+# 
 # 2339 2017-08-07 13:55:26Z gronemeier
 # corrected timestamp in header
-#
+# 
 # 2338 2017-08-07 12:15:38Z gronemeier
 # Modularize 1D model
-#
+# 
 # 2320 2017-07-21 12:47:43Z suehring
 # -ls_forcing nudging
 # +large_scale_forcing_nudging
-#
+# 
 # 2318 2017-07-20 17:27:44Z suehring
 # Add further dependencies on surface_mod
-#
+# 
 # 2317 2017-07-20 17:27:19Z suehring
 # Added time_integration_spinup
-#
+# 
 # 2269 2017-06-09 11:57:32Z suehring
 # Add dependency in read_3d_binary
-#
+# 
 # 2263 2017-06-08 14:59:01Z schwenkel
 # Implemented splitting and merging algorithm
-#
+# 
 # 2259 2017-06-08 09:09:11Z gronemeier
 # Implemented synthetic turbulence generator
@@ -242,12 +243,12 @@
 # 2256 2017-06-07 13:58:08Z suehring
 # Remove ring dependency in init_pegrid
-#
+# 
 # 2238 2017-05-31 16:49:16Z suehring
 # Bugfix, further missing dependency on surface_mod
-#
+# 
 # 2237 2017-05-31 10:34:53Z suehring
-# Bugfix, add dependencies on surface_mod for surface_coupler,
-# plant_canopy_model_mod and ls_forcing_mod
-#
+# Bugfix, add dependencies on surface_mod for surface_coupler, 
+# plant_canopy_model_mod and ls_forcing_mod 
+# 
 # 2233 2017-05-30 18:08:54Z suehring
 #
@@ -255,5 +256,5 @@
 # +dependencies for surface_mod
 # -wall_fluxes
-#
+# 
 # 2130 2017-01-24 16:25:39Z raasch
 # dependency for timestep updated
@@ -261,13 +262,13 @@
 # 2118 2017-01-17 16:38:49Z raasch
 # -cuda_fft_interfaces_mod
-#
+# 
 # 2050 2016-11-08 15:00:55Z gronemeier
 # Implement turbulent outflow method
-#
+# 
 # 2007 2016-08-24 15:47:17Z kanani
 # urban surface module added,
 # cleaned up some lines (compiler flags/options), which were accidentally
 # added in rev1938
-#
+# 
 # 1998 2016-08-20 18:45:34Z knoop
 # Bugfix: added netcdf_interface to dependency list for user_init_land_surface
@@ -275,8 +276,8 @@
 # 1986 2016-08-10 14:07:17Z gronemeier
 # POSIX-calls module added
-#
+# 
 # 1972 2016-07-26 07:52:02Z maronga
 # Removed some dependencies due to further modularization of land surface model
-#
+# 
 # 1957 2016-07-07 10:43:48Z suehring
 # flight module added
@@ -286,9 +287,9 @@
 #
 # 1938 2016-06-13 15:26:05Z hellstea
-# Some dependency errors corrected
-#
+# Some dependency errors corrected 
+# 
 # 1934 2016-06-13 09:46:57Z hellstea
 # poismg renamed poismg_noopt, poismg_fast_mod renamed poismg_mod
-#
+# 
 # 1914 2016-05-26 14:44:07Z witha
 # Added wind_turbine_model_mod
@@ -301,5 +302,5 @@
 #
 # 1850 2016-04-08 13:29:27Z maronga
-# Adapted for modularization of microphysics
+# Adapted for modularization of microphysics 
 # Several files renamed --> _mod
 # Bugfix for previous commit
@@ -316,16 +317,16 @@
 # 1833 2016-04-07 14:23:03Z raasch
 # spectrum renamed spectra_mod, depencies for spectra changed
-#
+# 
 # 1826 2016-04-07 12:01:39Z maronga
 # Renamed radiation_model to radiation_model_mod.
 # Renamed plant_canopy_model to plant_canopy_model_mod.
-#
+# 
 # 1822 2016-04-07 07:49:42Z hoffmann
 # Tails removed. lpm_release_set removed. calc_precipitation, impact_of_latent_heat
-# removed.
+# removed. 
 #
 # 1817 2016-04-06 15:44:20Z maronga
 # Renamed land_surface_model to land_surface_model_mod.
-#
+# 
 # 1808 2016-04-05 19:44:00Z raasch
 # -local_flush, -local_getenv
@@ -353,5 +354,5 @@
 #
 # 1762 2016-02-25 12:31:13Z hellstea
-# +pmc_interface, +pmc routines
+# +pmc_interface, +pmc routines 
 #
 # 1747 2016-02-08 12:25:53Z raasch
@@ -359,10 +360,10 @@
 #
 # 1691 2015-10-26 16:17:44Z maronga
-# Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to
+# Replaced prandtl_fluxes with surface_layer_fluxes. Added radiation_model to 
 # prognostic_equations
-#
+# 
 # 1585 2015-04-30 07:05:52Z maronga
 # Added user_init_radiation.f90
-#
+# 
 # 1575 2015-03-27 09:56:27Z raasch
 # +poismg_fast
@@ -370,5 +371,5 @@
 # 1551 2015-03-03 14:18:16Z maronga
 # Bugfix: further adjustments for the land surface model and radiation model
-#
+# 
 # 1517 2015-01-07 19:12:25Z hoffmann
 # advec_s_bc added to prognostic_equations
@@ -376,16 +377,16 @@
 # 1500 2014-12-03 17:42:41Z maronga
 # Bugfix: missing adjustments for land surface model and radiation model
-#
+# 
 # 1496 2014-12-02 17:25:50Z maronga
-# Added land surface model and radiation model files: land_surface_model,
+# Added land surface model and radiation model files: land_surface_model, 
 # radiation_model, user_init_land_surface
-#
+# 
 # 1484 2014-10-21 10:53:05Z kanani
 # plant_canopy_model-dependency added for check_parameters, header, init_3d_model,
 # package_parin, read_var_list, user_init_plant_canopy, write_var_list
-#
+# 
 # 1444 2014-08-02 20:10:32Z letzel
 # bugfix: cpulog added to lpm_advec
-#
+# 
 # 1404 2014-05-14 09:01:39Z keck
 # bugfix: dependencies added for progress_bar
@@ -393,11 +394,11 @@
 # 1402 2014-05-09 14:25:13Z raasch
 # progress_bar added
-#
+# 
 # 1400 2014-05-09 14:03:54Z knoop
 # Added new module random_generator_parallel
-#
+# 
 # 1380 2014-04-28 12:40:45Z heinze
 # bugfix: mod_particle_attributes  added to check_open
-# nudging added to time_integration
+# nudging added to time_integration 
 #
 # 1374 2014-04-25 12:55:07Z raasch
@@ -407,5 +408,5 @@
 # Added new module calc_mean_profile, previously in module buoyancy,
 # removed buoyancy dependency from nudging
-#
+# 
 # 1363 2014-04-17 12:28:49Z keck
 # bugfix: cpulog added to lpm_pack_arrays
@@ -413,7 +414,7 @@
 # 1361 2014-04-16 15:17:48Z hoffmann
 # cpulog added to microphysics
-#
+# 
 # 1359 2014-04-11 17:15:14Z hoffmann
-# mod_particle_attributes added, lpm_sort_arrays removed,
+# mod_particle_attributes added, lpm_sort_arrays removed, 
 # lpm_extend_particle_array removed
 #
@@ -722,5 +723,5 @@
 
 CC = cc
-CFLAGS = -O
+CFLAGS = -O 
 
 F90 =
@@ -1349,4 +1350,5 @@
 	chemistry_model_mod.o \
 	gust_mod.o \
+	indoor_model_mod.o \
 	land_surface_model_mod.o \
 	mod_kinds.o \
@@ -1750,5 +1752,4 @@
 	modules.o \
 	netcdf_data_input_mod.o \
-	plant_canopy_model_mod.o \
 	radiation_model_mod.o \
 	surface_mod.o
Index: /palm/trunk/SOURCE/indoor_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/indoor_model_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/indoor_model_mod.f90	(revision 3744)
@@ -21,5 +21,10 @@
 ! Current revisions:
 ! -----------------
-! 
+! - remove building_type from module
+! - initialize parameters for each building individually instead of a bulk 
+!   initializaion with  identical building type for all
+! - output revised
+! - add missing _wp
+! - some restructuring of variables in building data structure
 ! 
 ! Former revisions:
@@ -76,4 +81,7 @@
 
     USE kinds
+    
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  building_id_f, building_type_f
 
     USE surface_mod,                                                           &
@@ -105,6 +113,31 @@
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  num_facade_v   !< number of vertical facades elements per buidling
                                                                   !< and height level
+                                                                  
+       INTEGER(iwp) ::  ventilation_int_loads
 
        LOGICAL ::  on_pe = .FALSE.   !< flag indicating whether a building with certain ID is on local subdomain
+       
+       
+       REAL(wp) ::  lambda_layer3     !< [W/(m*K)] Thermal conductivity of the inner layer  
+       REAL(wp) ::  s_layer3          !< [m] half thickness of the inner layer (layer_3)
+       REAL(wp) ::  f_c_win           !< [-] shading factor
+       REAL(wp) ::  g_value_win       !< [-] SHGC factor
+       REAL(wp) ::  u_value_win       !< [W/(m2*K)] transmittance
+       REAL(wp) ::  air_change_low    !< [1/h] air changes per time_utc_hour
+       REAL(wp) ::  air_change_high   !< [1/h] air changes per time_utc_hour
+       REAL(wp) ::  eta_ve            !< [-] heat recovery efficiency
+       REAL(wp) ::  factor_a          !< [-] Dynamic parameters specific effective surface according to Table 12; 2.5 
+                                      !< (very light, light and medium), 3.0 (heavy), 3.5 (very heavy)
+       REAL(wp) ::  factor_c          !< [J/(m2 K)] Dynamic parameters inner heatstorage according to Table 12; 80000 
+                                      !< (very light), 110000 (light), 165000 (medium), 260000 (heavy), 370000 (very heavy)
+       REAL(wp) ::  lambda_at         !< [-] ratio internal surface/floor area chap. 7.2.2.2.
+       REAL(wp) ::  theta_int_h_set   !< [degree_C] Max. Setpoint temperature (winter)
+       REAL(wp) ::  theta_int_c_set   !< [degree_C] Max. Setpoint temperature (summer)
+       REAL(wp) ::  phi_h_max         !< [W] Max. Heating capacity (negative)
+       REAL(wp) ::  phi_c_max         !< [W] Max. Cooling capacity (negative)
+       REAL(wp) ::  qint_high         !< [W/m2] internal heat gains, option Database qint_0-23
+       REAL(wp) ::  qint_low          !< [W/m2] internal heat gains, option Database qint_0-23
+       REAL(wp) ::  height_storey     !< [m] storey heigth
+       REAL(wp) ::  height_cei_con    !< [m] ceiling construction heigth
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in       !< mean building indoor temperature, height dependent
@@ -113,5 +146,5 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  vol_frac   !< fraction of local on total building volume, height dependent
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  vpf        !< building volume volume per facade element, height dependent
-
+       
     END TYPE build
 
@@ -130,5 +163,5 @@
 !-- Declare all global variables within the module
 
-    INTEGER(iwp) ::  building_type = 1       !< namelist parameter with 
+!     INTEGER(iwp) ::  building_type = 1       !< namelist parameter with 
                                              !< X1=construction year (cy) 1950, X2=cy 2000, X3=cy 2050
                                              !< R=Residental building, O=Office, RW=Enlarged Windows, P=Panel type (Plattenbau) WBS 70, H=Hospital (in progress), I=Industrial halls (in progress), S=Special Building (in progress)
@@ -139,21 +172,13 @@
     INTEGER(iwp) ::  solar_protection_on     !< Solar protection on
 
-    REAL(wp) ::  air_change_high             !< [1/h] air changes per time_utc_hour
-    REAL(wp) ::  air_change_low              !< [1/h] air changes per time_utc_hour
-    REAL(wp) ::  eff_mass_area               !< [mÂ²] the effective mass-related area
-    REAL(wp) ::  floor_area_per_facade       !< [mÂ²] net floor area (Sum of all floors)
-    REAL(wp) ::  total_area                  !<! [mÂ²] area of all surfaces pointing to zone
+    REAL(wp) ::  eff_mass_area               !< [m2] the effective mass-related area
+    REAL(wp) ::  floor_area_per_facade       !< [m2] net floor area (Sum of all floors)
+    REAL(wp) ::  total_area                  !<! [m2] area of all surfaces pointing to zone
     REAL(wp) ::  window_area_per_facade      !< [m2] window area per facade element
     REAL(wp) ::  air_change                  !< [1/h] Airflow
-    REAL(wp) ::  bldg_part_surf_i    = 4     !< [mÂ²_surf,i] part building surface, from Palm, das mÃŒsste mittlerweile "facade_element_area" sein!
-    REAL(wp) ::  facade_element_area         !< [mÂ²_facade] building surface facade
-    REAL(wp) ::  indoor_volume_per_facade    !< [mÂ³] indoor air volume per facade element
+    REAL(wp) ::  facade_element_area         !< [m2_facade] building surface facade
+    REAL(wp) ::  indoor_volume_per_facade    !< [m3] indoor air volume per facade element
     REAL(wp) ::  c_m                         !< [J/K] internal heat storage capacity
     REAL(wp) ::  dt_indoor = 3600.0_wp       !< [s] namelist parameter: time interval for indoor-model application
-    REAL(wp) ::  eta_ve                      !< [-] heat recovery efficiency
-    REAL(wp) ::  f_c_win                     !< [-] shading factor
-    REAL(wp) ::  factor_a                    !< [-] Dynamic parameters specific effective surface according to Table 12; 2.5 (very light, light and medium), 3.0 (heavy), 3.5 (very heavy)
-    REAL(wp) ::  factor_c                    !< [J/(m2 K)] Dynamic parameters inner heatstorage according to Table 12; 80000 (very light), 110000 (light), 165000 (medium), 260000 (heavy), 370000 (very heavy)
-    REAL(wp) ::  g_value_win                 !< [-] SHGC factor
     REAL(wp) ::  h_tr_1                      !<! [W/K] Heat transfer coefficient auxiliary variable 1
     REAL(wp) ::  h_tr_2                      !<! [W/K] Heat transfer coefficient auxiliary variable 2
@@ -165,14 +190,6 @@
     REAL(wp) ::  h_tr_w                      !<! [W/K] heat transfer coefficient of doors, windows, curtain walls and glazed walls (assumption: thermal mass=0)
     REAL(wp) ::  h_ve                        !<! [W/K] heat transfer of ventilation
-    REAL(wp) ::  height_storey               !< [m] storey heigth
-    REAL(wp) ::  height_cei_con              !< [m] ceiling construction heigth    
     REAL(wp) ::  initial_indoor_temperature  !< namelist parameter
-    REAL(wp) ::  lambda_at                   !< [-] ratio internal surface/floor area chap. 7.2.2.2.
-    REAL(wp) ::  lambda_layer3               !< [W/(m*K)] Thermal conductivity of the inner layer  
     REAL(wp) ::  net_sw_in                   !< net short-wave radiation (in - out; was i_global --> CORRECT?)
-    REAL(wp) ::  qint_high                   !< [W/m2] internal heat gains, option Database qint_0-23
-    REAL(wp) ::  qint_low                    !< [W/m2] internal heat gains, option Database qint_0-23
-    REAL(wp) ::  phi_c_max                   !< [W] Max. Cooling capacity (negative)
-    REAL(wp) ::  phi_h_max                   !< [W] Max. Heating capacity (negative)
     REAL(wp) ::  phi_hc_nd                   !<! [W] heating demand and/or cooling demand
     REAL(wp) ::  phi_hc_nd_10                !<! [W] heating demand and/or cooling demand for heating or cooling 
@@ -185,17 +202,14 @@
     REAL(wp) ::  phi_st                      !<! [W] mass specific thermal load implied non thermal mass 
     REAL(wp) ::  q_emission                  !< emissions, in first version = 0, option for second part of the project
-    REAL(wp) ::  q_wall_win               !< heat flux from indoor into wall/window
+    REAL(wp) ::  q_wall_win                  !< heat flux from indoor into wall/window
     REAL(wp) ::  q_waste_heat                !< waste heat, sum of waste heat over the roof to Palm
     REAL(wp) ::  q_waste_heat_bldg           !< [W/building] waste heat of the complete building, in Palm sum of all indoor_model-calculations
-    REAL(wp) ::  s_layer3                    !< [m] half thickness of the inner layer (layer_3)
     REAL(wp) ::  schedule_d                  !< activation for internal loads (low or high + low)
     REAL(wp) ::  skip_time_do_indoor = 0.0_wp  !< [s] Indoor model is not called before this time
     REAL(wp) ::  theta_air                   !<! [degree_C] air temperature of the RC-node
     REAL(wp) ::  theta_air_0                 !<! [degree_C] air temperature of the RC-node in equilibrium
-    REAL(wp) ::  theta_air_10                !<! [degree_C] air temperature of the RC-node from a heating capacity of 10 W/mÂ²
+    REAL(wp) ::  theta_air_10                !<! [degree_C] air temperature of the RC-node from a heating capacity of 10 W/m2
     REAL(wp) ::  theta_air_ac                !< [degree_C] actual room temperature after heating/cooling
     REAL(wp) ::  theta_air_set               !< [degree_C] Setpoint_temperature for the room
-    REAL(wp) ::  theta_int_c_set             !< [degree_C] Max. Setpoint temperature (summer)
-    REAL(wp) ::  theta_int_h_set             !< [degree_C] Max. Setpoint temperature (winter)
     REAL(wp) ::  theta_m                     !<! [degree_C} inner temperature of the RC-node
     REAL(wp) ::  theta_m_t                   !<! [degree_C] (Fictive) component temperature timestep
@@ -205,9 +219,11 @@
     REAL(wp) ::  time_indoor = 0.0_wp        !< [s] time since last call of indoor model
     REAL(wp) ::  time_utc_hour               !< Time in hours per day (UTC)
-    REAL(wp) ::  u_value_win                 !< [W/(m2*K)] transmittance
     REAL(wp) ::  ventilation_int_loads       !< Zuteilung der GebÃ€ude fÃŒr Verlauf/AktivitÃ€t der LÃŒftung und internen Lasten
-
-!
-!-- Declare all global parameters within the module    
+    
+    REAL(wp) ::  f_sr                        !< [-] factor surface reduction
+    REAL(wp) ::  f_cei                       !< [-] ceiling reduction factor
+    REAL(wp) ::  ngs                         !< [m2] netto ground surface
+    REAL(wp) ::  building_height
+    
     REAL(wp), PARAMETER ::  params_f_f               = 0.3_wp      !< [-] frame ratio chap. 8.3.2.1.1 for buildings with mostly cooling 2.0_wp
     REAL(wp), PARAMETER ::  params_f_w               = 0.9_wp      !< [-] correction factor (fuer nicht senkrechten Stahlungseinfall DIN 4108-2 chap.8, (hier konstant, keine WinkelabhÃ€ngigkeit) 
@@ -218,5 +234,7 @@
     REAL(wp), PARAMETER ::  h_is                     = 3.45_wp     !< [W/(m^2 K)]  h_is = 3.45 between surface and air (chap. 7.2.2.2)
     REAL(wp), PARAMETER ::  h_ms                     = 9.1_wp      !< [W/K] h_ms = 9.10 W / (m2 K) between component and surface (chap. 12.2.2)
-  
+
+    
+    
     SAVE
 
@@ -226,5 +244,7 @@
 !
 !-- Add INTERFACES that must be available to other modules
-    PUBLIC im_init, im_main_heatcool, im_parin 
+    PUBLIC im_init, im_main_heatcool, im_parin, im_define_netcdf_grid,          &
+           im_check_data_output, im_data_output_3d, im_check_parameters
+    
 
 !
@@ -232,4 +252,37 @@
     PUBLIC dt_indoor, skip_time_do_indoor, time_indoor
 
+!
+!-- PALM interfaces:
+!-- Data output checks for 2D/3D data to be done in check_parameters
+     INTERFACE im_check_data_output
+        MODULE PROCEDURE im_check_data_output
+     END INTERFACE im_check_data_output
+!
+!-- Input parameter checks to be done in check_parameters
+     INTERFACE im_check_parameters
+        MODULE PROCEDURE im_check_parameters
+     END INTERFACE im_check_parameters
+!
+!-- Data output of 3D data
+     INTERFACE im_data_output_3d
+        MODULE PROCEDURE im_data_output_3d
+     END INTERFACE im_data_output_3d
+
+!
+!-- Definition of data output quantities
+     INTERFACE im_define_netcdf_grid
+        MODULE PROCEDURE im_define_netcdf_grid
+     END INTERFACE im_define_netcdf_grid
+! 
+! !
+! !-- Output of information to the header file
+!     INTERFACE im_header
+!        MODULE PROCEDURE im_header
+!     END INTERFACE im_header
+
+!-- Data Output
+!    INTERFACE im_data_output
+!       MODULE PROCEDURE im_data_output
+!    END INTERFACE im_data_output
 !
 !-- Calculations for indoor temperatures  
@@ -242,5 +295,4 @@
        MODULE PROCEDURE im_init
     END INTERFACE im_init
- 
 !
 !-- Main part of indoor model  
@@ -268,6 +320,6 @@
     USE arrays_3d,                                                             &
         ONLY:  pt
-
-
+    
+    
     IMPLICIT NONE
     
@@ -280,4 +332,5 @@
     REAL(wp) ::  near_facade_temperature
     REAL(wp) ::  phi_hc_nd_dummy
+   
 
     !< Calculation of total mass specific thermal load (internal and external)
@@ -292,11 +345,11 @@
     !< Calculation of component temperature at factual timestep
     theta_m_t = ( ( theta_m_t_prev                                               &
-                    * ( ( c_m / 3600 ) - 0.5 * ( h_tr_3 + h_tr_em ) ) + phi_mtot &
+                    * ( ( c_m / 3600.0_wp ) - 0.5_wp * ( h_tr_3 + h_tr_em ) ) + phi_mtot &
                   )                                                              &
-                  /   ( ( c_m / 3600 ) + 0.5 * ( h_tr_3 + h_tr_em ) )            &
+                  /   ( ( c_m / 3600.0_wp ) + 0.5_wp * ( h_tr_3 + h_tr_em ) )            &
                 )                                                               !< [degree_C] Eq. (C.4)
 
     !< Calculation of mean inner temperature for the RC-node in actual timestep
-    theta_m = ( theta_m_t + theta_m_t_prev ) * 0.5                              !< [degree_C] Eq. (C.9)
+    theta_m = ( theta_m_t + theta_m_t_prev ) * 0.5_wp                              !< [degree_C] Eq. (C.9)
     
     !< Calculation of mean surface temperature of the RC-node in actual timestep
@@ -343,7 +396,4 @@
         ONLY:  dx, dy
 
-    USE netcdf_data_input_mod,                                                 &
-        ONLY:  building_id_f
-
     USE pegrid
 
@@ -356,4 +406,5 @@
     IMPLICIT NONE
 
+    INTEGER(iwp) ::  bt   !< local building type
     INTEGER(iwp) ::  fa   !< running index for facade elements of each building
     INTEGER(iwp) ::  i    !< running index along x-direction
@@ -384,5 +435,5 @@
     INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings         !< number of buildings with different ID on entire model domain
     INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings_l       !< number of buildings with different ID on local subdomain
-
+    
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  local_weight   !< dummy representing fraction of local on total building volume, 
                                                            !< height dependent
@@ -503,4 +554,5 @@
 !-- the respective building is present (in order to reduce memory demands).
     ALLOCATE( buildings(1:num_build) )
+
 !
 !-- Store building IDs and check if building with certain ID is present on
@@ -792,30 +844,112 @@
 
 !
-!-- Building parameters by type of building. Assigned in urban_surface_mod.f90
-
-    lambda_layer3            = building_pars(63, building_type)
-    s_layer3                 = building_pars(57, building_type)
-    f_c_win                  = building_pars(119, building_type)
-    g_value_win              = building_pars(120, building_type)    
-    u_value_win              = building_pars(121, building_type)    
-    air_change_low           = building_pars(122, building_type)    
-    air_change_high          = building_pars(123, building_type)    
-    eta_ve                   = building_pars(124, building_type)    
-    factor_a                 = building_pars(125, building_type)    
-    factor_c                 = building_pars(126, building_type)
-    lambda_at                = building_pars(127, building_type)    
-    theta_int_h_set          = building_pars(118, building_type)    
-    theta_int_c_set          = building_pars(117, building_type)
-    phi_h_max                = building_pars(128, building_type)    
-    phi_c_max                = building_pars(129, building_type)          
-    qint_high                = building_pars(130, building_type)
-    qint_low                 = building_pars(131, building_type)
-    height_storey            = building_pars(132, building_type)
-    height_cei_con           = building_pars(133, building_type)
-
-!
-!-- Setting of initial room temperature [K]
+!-- Initialize building parameters, first by mean building type. Note,
+!-- in this case all buildings have the same type. 
+!-- In a second step initialize with building tpyes from static input file, 
+!-- where building types can be individual for each building.
+    buildings(:)%lambda_layer3    = building_pars(63,building_type)
+    buildings(:)%s_layer3         = building_pars(57,building_type)
+    buildings(:)%f_c_win          = building_pars(119,building_type)
+    buildings(:)%g_value_win      = building_pars(120,building_type)    
+    buildings(:)%u_value_win      = building_pars(121,building_type)    
+    buildings(:)%air_change_low   = building_pars(122,building_type)    
+    buildings(:)%air_change_high  = building_pars(123,building_type)    
+    buildings(:)%eta_ve           = building_pars(124,building_type)    
+    buildings(:)%factor_a         = building_pars(125,building_type)    
+    buildings(:)%factor_c         = building_pars(126,building_type)
+    buildings(:)%lambda_at        = building_pars(127,building_type)    
+    buildings(:)%theta_int_h_set  = building_pars(118,building_type)    
+    buildings(:)%theta_int_c_set  = building_pars(117,building_type)
+    buildings(:)%phi_h_max        = building_pars(128,building_type)    
+    buildings(:)%phi_c_max        = building_pars(129,building_type)          
+    buildings(:)%qint_high        = building_pars(130,building_type)
+    buildings(:)%qint_low         = building_pars(131,building_type)
+    buildings(:)%height_storey    = building_pars(132,building_type)
+    buildings(:)%height_cei_con   = building_pars(133,building_type)
+!
+!-- Initialize ventilaation load. Please note, building types > 7 are actually
+!-- not allowed (check already in urban_surface_mod and netcdf_data_input_mod. 
+!-- However, the building data base may be later extended. 
+    IF ( building_type ==  1  .OR.  building_type ==  2  .OR.                  &
+         building_type ==  3  .OR.  building_type == 10  .OR.                  &
+         building_type == 11  .OR.  building_type == 12 )  THEN
+       buildings(nb)%ventilation_int_loads = 1
+!
+!-- Office, building with large windows
+    ELSEIF ( building_type ==  4  .OR.  building_type ==  5  .OR.              &
+             building_type ==  6  .OR.  building_type ==  7  .OR.              &
+             building_type ==  8  .OR.  building_type ==  9)  THEN
+       buildings(nb)%ventilation_int_loads = 2
+!
+!-- Industry, hospitals
+    ELSEIF ( building_type == 13  .OR.  building_type == 14  .OR.              &
+             building_type == 15  .OR.  building_type == 16  .OR.              &
+             building_type == 17  .OR.  building_type == 18 )  THEN
+       buildings(nb)%ventilation_int_loads = 3
+    ENDIF
+!
+!-- Initialization of building parameters - level 2
+    IF ( building_type_f%from_file )  THEN
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+              IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
+                 nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ),    &
+                              DIM = 1 )
+                 bt = building_type_f%var(j,i)
+                 
+                 buildings(nb)%lambda_layer3    = building_pars(63,bt)
+                 buildings(nb)%s_layer3         = building_pars(57,bt)
+                 buildings(nb)%f_c_win          = building_pars(119,bt)
+                 buildings(nb)%g_value_win      = building_pars(120,bt)    
+                 buildings(nb)%u_value_win      = building_pars(121,bt)    
+                 buildings(nb)%air_change_low   = building_pars(122,bt)    
+                 buildings(nb)%air_change_high  = building_pars(123,bt)    
+                 buildings(nb)%eta_ve           = building_pars(124,bt)    
+                 buildings(nb)%factor_a         = building_pars(125,bt)    
+                 buildings(nb)%factor_c         = building_pars(126,bt)
+                 buildings(nb)%lambda_at        = building_pars(127,bt)    
+                 buildings(nb)%theta_int_h_set  = building_pars(118,bt)    
+                 buildings(nb)%theta_int_c_set  = building_pars(117,bt)
+                 buildings(nb)%phi_h_max        = building_pars(128,bt)    
+                 buildings(nb)%phi_c_max        = building_pars(129,bt)          
+                 buildings(nb)%qint_high        = building_pars(130,bt)
+                 buildings(nb)%qint_low         = building_pars(131,bt)
+                 buildings(nb)%height_storey    = building_pars(132,bt)
+                 buildings(nb)%height_cei_con   = building_pars(133,bt)                 
+!
+!--              Initialize ventilaation load. Please note, building types > 7 
+!--              are actually not allowed (check already in urban_surface_mod  
+!--              and netcdf_data_input_mod. However, the building data base may 
+!--              be later extended. 
+                 IF ( bt ==  1  .OR.  bt ==  2  .OR.                           &
+                      bt ==  3  .OR.  bt == 10  .OR.                           &
+                      bt == 11  .OR.  bt == 12 )  THEN
+                    buildings(nb)%ventilation_int_loads = 1
+!                    
+!--              Office, building with large windows
+                 ELSEIF ( bt ==  4  .OR.  bt ==  5  .OR.                       &
+                          bt ==  6  .OR.  bt ==  7  .OR.                       &
+                          bt ==  8  .OR.  bt ==  9)  THEN
+                    buildings(nb)%ventilation_int_loads = 2
+!
+!--              Industry, hospitals
+                 ELSEIF ( bt == 13  .OR.  bt == 14  .OR.                       &
+                          bt == 15  .OR.  bt == 16  .OR.                       &
+                          bt == 17  .OR.  bt == 18 )  THEN
+                    buildings(nb)%ventilation_int_loads = 3
+                 ENDIF
+              ENDIF
+           ENDDO
+        ENDDO
+    ENDIF
+!
+!-- Initial room temperature [K]
 !-- (after first loop, use theta_m_t as theta_m_t_prev)
     theta_m_t_prev = initial_indoor_temperature
+!
+!-- Initialize indoor temperature. Actually only for output at initial state.
+    DO  nb = 1, num_build
+       buildings(nb)%t_in(:) = initial_indoor_temperature
+    ENDDO
 
     CALL location_message( 'finished', .TRUE. )
@@ -874,74 +1008,7 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in_l_send   !< dummy send buffer used for summing-up indoor temperature per kk-level
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in_recv     !< dummy recv buffer used for summing-up indoor temperature per kk-level
-
-!
-!-- Daily schedule, here for 08:00-18:00 = 1, other hours = 0.
-!-- time_utc_hour is calculated here based on time_utc [s] from 
-!-- date_and_time_mod.
-!-- (kanani: Does this schedule not depend on if it's an office or resident 
-!-- building?)
+!
+!-- Determine time of day in hours. 
     time_utc_hour = time_utc / 3600.0_wp
-
-!
-!-- Allocation of the load profiles to the building types
-!-- Residental Building, panel WBS 70
-    if      (building_type ==  1 .OR. &
-             building_type ==  2 .OR. &
-             building_type ==  3 .OR. &
-             building_type == 10 .OR. &
-             building_type == 11 .OR. &
-             building_type == 12) then
-                                        ventilation_int_loads = 1
-!-- Office, building with large windows
-    else if (building_type ==  4 .OR. &
-             building_type ==  5 .OR. &
-             building_type ==  6 .OR. &
-             building_type ==  7 .OR. &
-             building_type ==  8 .OR. &
-             building_type ==  9) then
-                                        ventilation_int_loads = 2
-!-- Industry, hospitals
-    else if (building_type == 13 .OR. &
-             building_type == 14 .OR. &
-             building_type == 15 .OR. &
-             building_type == 16 .OR. &
-             building_type == 17 .OR. &
-             building_type == 18) then
-                                        ventilation_int_loads = 3
-
-    end if
-
-!-- Residental Building, panel WBS 70
-    
-    if (ventilation_int_loads == 1) THEN
-       if ( time_utc_hour >= 6.0_wp .AND. time_utc_hour <= 8.0_wp )  THEN
-          schedule_d = 1
-       else if ( time_utc_hour >= 18.0_wp .AND. time_utc_hour <= 23.0_wp )  THEN
-          schedule_d = 1
-       else
-          schedule_d = 0
-       end if
-    end if
-
-!-- Office, building with large windows
-
-    if (ventilation_int_loads == 2) THEN
-       if ( time_utc_hour >= 8.0_wp  .AND.  time_utc_hour <= 18.0_wp )  THEN
-          schedule_d = 1
-       else
-          schedule_d = 0
-       end if
-    end if
-       
-!-- Industry, hospitals
-    if (ventilation_int_loads == 3) THEN
-       if ( time_utc_hour >= 6.0_wp  .AND.  time_utc_hour <= 22.0_wp )  THEN
-          schedule_d = 1
-       else
-          schedule_d = 0
-       end if
-    end if
-
-
 !
 !-- Following calculations must be done for each facade element.
@@ -951,6 +1018,35 @@
        IF ( buildings(nb)%on_pe )  THEN
 !
+!--       Determine daily schedule. 08:00-18:00 = 1, other hours = 0.
+!--       Residental Building, panel WBS 70    
+          IF ( buildings(nb)%ventilation_int_loads == 1 )  THEN
+             IF ( time_utc_hour >= 6.0_wp  .AND.  time_utc_hour <= 8.0_wp )  THEN
+                schedule_d = 1
+             ELSEIF ( time_utc_hour >= 18.0_wp  .AND.  time_utc_hour <= 23.0_wp )  THEN
+                schedule_d = 1
+             ELSE
+                schedule_d = 0
+             ENDIF
+          ENDIF
+!
+!--       Office, building with large windows
+          IF ( buildings(nb)%ventilation_int_loads == 2 )  THEN
+             IF ( time_utc_hour >= 8.0_wp  .AND.  time_utc_hour <= 18.0_wp )  THEN
+                schedule_d = 1
+             ELSE
+                schedule_d = 0
+             ENDIF
+          ENDIF
+!       
+!--       Industry, hospitals
+          IF ( buildings(nb)%ventilation_int_loads == 3 )  THEN
+             IF ( time_utc_hour >= 6.0_wp  .AND.  time_utc_hour <= 22.0_wp )  THEN
+                schedule_d = 1
+             ELSE
+                schedule_d = 0
+             ENDIF
+          ENDIF
+!
 !--       Initialize/reset indoor temperature
-          buildings(nb)%t_in   = 0.0_wp
           buildings(nb)%t_in_l = 0.0_wp 
 !
@@ -970,24 +1066,37 @@
              indoor_volume_per_facade = buildings(nb)%vpf(kk)               !< [m3] indoor air volume per facade element           
              window_area_per_facade   = surf_usm_h%frac(ind_wat_win,m)  * facade_element_area  !< [m2] window area per facade element
-             eff_mass_area            = factor_a * floor_area_per_facade    !< [m2] standard values according to Table 12 section 12.3.1.2  (calculate over Eq. (65) according to section 12.3.1.2)
-             c_m                      = factor_c * floor_area_per_facade    !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
-             total_area               = lambda_at * floor_area_per_facade   !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
+             
+!              building_height          = buildings(nb)%num_facades_per_building_v_l * 0.1 * dzw(kk)
+             building_height          = buildings(nb)%kb_max * dzw(kk)
+             
+! print*, "building_height", building_height
+! print*, "num_facades_v_l", buildings(nb)%num_facades_per_building_v_l
+! print*, "num_facades_v", buildings(nb)%num_facades_per_building_v
+! print*, "kb_min_max", buildings(nb)%kb_min, buildings(nb)%kb_max
+! print*, "dzw kk", dzw(kk), kk
+
+             f_cei                    = building_height/(buildings(nb)%height_storey-buildings(nb)%height_cei_con) !< [-] factor for ceiling redcution
+             ngs                      = buildings(nb)%vpf(kk)/f_cei                    !< [m2] calculation of netto ground surface
+             f_sr                     = ngs/floor_area_per_facade                      !< [-] factor for surface reduction
+             eff_mass_area            = buildings(nb)%factor_a * ngs    !< [m2] standard values according to Table 12 section 12.3.1.2  (calculate over Eq. (65) according to section 12.3.1.2)
+             c_m                      = buildings(nb)%factor_c * ngs    !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
+             total_area               = buildings(nb)%lambda_at * floor_area_per_facade   !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
 
 !--          Calculation of heat transfer coefficient for transmission --> not time-dependent
-             h_tr_w   = window_area_per_facade * u_value_win   !< [W/K] only for windows
+             h_tr_w   = window_area_per_facade * buildings(nb)%u_value_win   !< [W/K] only for windows
              h_tr_is  = total_area * h_is                      !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
              h_tr_ms  = eff_mass_area * h_ms                    !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
-             h_tr_op  = 1 / ( 1 / ( ( facade_element_area - window_area_per_facade ) &
-                                    * lambda_layer3 / s_layer3 * 0.5 ) + 1 / h_tr_ms )
-             h_tr_em  = 1 / ( 1 / h_tr_op - 1 / h_tr_ms )      !< [W/K] Eq. (63), Section 12.2.2
+             h_tr_op  = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade ) &
+                                    * buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp ) + 1.0_wp / h_tr_ms )
+             h_tr_em  = 1.0_wp / ( 1.0_wp / h_tr_op - 1.0_wp / h_tr_ms )      !< [W/K] Eq. (63), Section 12.2.2
 !
 !--          internal air loads dependent on the occupacy of the room
 !--          basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
-             phi_ia = 0.5 * ( ( qint_high * schedule_d + qint_low )            &
-                              * floor_area_per_facade )         !< [W] Eq. (C.1)
+             phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )            &
+                              * ngs )         !< [W] Eq. (C.1)
 !
 !--          Airflow dependent on the occupacy of the room
 !--          basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
-             air_change = ( air_change_high * schedule_d + air_change_low )  !< [1/h]? 
+             air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low )  !< [1/h]? 
 !
 !--          Heat transfer of ventilation 
@@ -995,10 +1104,10 @@
 !--          with heat capacity of air 0.33 Wh/m2K
              h_ve   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *      &
-                                    0.33_wp * (1 - eta_ve ) ) )    !< [W/K] from ISO 13789 Eq.(10)
+                                    0.33_wp * (1.0_wp - buildings(nb)%eta_ve ) ) )    !< [W/K] from ISO 13789 Eq.(10)
 
 !--          Heat transfer coefficient auxiliary variables
-             h_tr_1 = 1 / ( ( 1 / h_ve )   + ( 1 / h_tr_is ) )  !< [W/K] Eq. (C.6)
+             h_tr_1 = 1.0_wp / ( ( 1.0_wp / h_ve )   + ( 1.0_wp / h_tr_is ) )  !< [W/K] Eq. (C.6)
              h_tr_2 = h_tr_1 + h_tr_w                           !< [W/K] Eq. (C.7)
-             h_tr_3 = 1 / ( ( 1 / h_tr_2 ) + ( 1 / h_tr_ms ) )  !< [W/K] Eq. (C.8)
+             h_tr_3 = 1.0_wp / ( ( 1.0_wp / h_tr_2 ) + ( 1.0_wp / h_tr_ms ) )  !< [W/K] Eq. (C.8)
 !
 !--          Net short-wave radiation through window area (was i_global)
@@ -1028,6 +1137,6 @@
 !--          DIN 4108 - 2 chap.8
              phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off               &
-                         + window_area_per_facade * net_sw_in * f_c_win * solar_protection_on )    &
-                       * g_value_win * ( 1 - params_f_f ) * params_f_w              !< [W]
+                         + window_area_per_facade * net_sw_in * buildings(nb)%f_c_win * solar_protection_on )    &
+                       * buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w              !< [W]
 !
 !--          Calculation of the mass specific thermal load for internal and external heatsources of the inner node
@@ -1035,5 +1144,5 @@
 !
 !--          Calculation mass specific thermal load implied non thermal mass 
-             phi_st  =   ( 1 - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1 * total_area ) ) ) &
+             phi_st  =   ( 1.0_wp - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1_wp * total_area ) ) ) &
                        * ( phi_ia + phi_sol )                                       !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
 !
@@ -1041,5 +1150,5 @@
 !--          Step 1: Indoor temperature without heating and cooling
 !--          section C.4.1 Picture C.2 zone 3)
-             phi_hc_nd = 0
+             phi_hc_nd = 0.0_wp
              
              CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
@@ -1047,10 +1156,10 @@
 !
 !--          If air temperature between border temperatures of heating and cooling, assign output variable, then ready   
-             IF ( theta_int_h_set <= theta_air  .AND.  theta_air <= theta_int_c_set )  THEN
-                phi_hc_nd_ac = 0
+             IF ( buildings(nb)%theta_int_h_set <= theta_air  .AND.  theta_air <= buildings(nb)%theta_int_c_set )  THEN
+                phi_hc_nd_ac = 0.0_wp
                 phi_hc_nd    = phi_hc_nd_ac
                 theta_air_ac = theta_air
 !
-!--          Step 2: Else, apply 10 W/mÂ² heating/cooling power and calculate indoor temperature
+!--          Step 2: Else, apply 10 W/m2 heating/cooling power and calculate indoor temperature
 !--          again.
              ELSE
@@ -1060,8 +1169,8 @@
 
 !--             Heating or cooling?
-                IF ( theta_air > theta_int_c_set )  THEN
-                   theta_air_set = theta_int_c_set
+                IF ( theta_air > buildings(nb)%theta_int_c_set )  THEN
+                   theta_air_set = buildings(nb)%theta_int_c_set
                 ELSE 
-                   theta_air_set = theta_int_h_set 
+                   theta_air_set = buildings(nb)%theta_int_h_set 
                 ENDIF
 
@@ -1079,10 +1188,8 @@
                                             / (theta_air_10  - theta_air_0)            !< Eq. (C.13)
                                             
-                                            
-                                        
 !--             Step 3: With temperature ratio to determine the heating or cooling capacity    
 !--             If necessary, limit the power to maximum power
 !--             section C.4.1 Picture C.2 zone 2) and 4)
-                IF ( phi_c_max < phi_hc_nd_un  .AND.  phi_hc_nd_un < phi_h_max )  THEN
+                IF ( buildings(nb)%phi_c_max < phi_hc_nd_un  .AND.  phi_hc_nd_un < buildings(nb)%phi_h_max )  THEN
                    phi_hc_nd_ac = phi_hc_nd_un
                    phi_hc_nd = phi_hc_nd_un    
@@ -1090,8 +1197,8 @@
 !--             Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
 !--             section C.4.1 Picture C.2 zone 1) and 5)
-                   IF ( phi_hc_nd_un > 0 )  THEN
-                      phi_hc_nd_ac = phi_h_max                                         !< Limit heating
+                   IF ( phi_hc_nd_un > 0.0_wp )  THEN
+                      phi_hc_nd_ac = buildings(nb)%phi_h_max                                         !< Limit heating
                    ELSE 
-                      phi_hc_nd_ac = phi_c_max                                         !< Limit cooling
+                      phi_hc_nd_ac = buildings(nb)%phi_c_max                                         !< Limit cooling
                    ENDIF
                 ENDIF
@@ -1112,5 +1219,5 @@
 !--          Calculate the operating temperature with weighted mean temperature of air and mean solar temperature
 !--          Will be used for thermal comfort calculations 
-             theta_op     = 0.3 * theta_air_ac + 0.7 * theta_s          !< [degree_C] operative Temperature Eq. (C.12)
+             theta_op     = 0.3_wp * theta_air_ac + 0.7_wp * theta_s          !< [degree_C] operative Temperature Eq. (C.12)
 !
 !--          Heat flux into the wall. Value needed in urban_surface_mod to 
@@ -1132,14 +1239,14 @@
 !--          Calculation of waste heat
 !--          Anthropogenic heat output
-              IF ( phi_hc_nd_ac > 0 )  THEN 
-                   heating_on = 1
-                   cooling_on = 0
-              ELSE 
-                   heating_on = 0
-                   cooling_on = 1
-              ENDIF
-
-             q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on))  !< [W/m2]  anthropogenic heat output
-!              surf_usm_h%shf(m)=q_waste_heat
+             IF ( phi_hc_nd_ac > 0.0_wp )  THEN 
+                heating_on = 1
+                cooling_on = 0
+             ELSE 
+                heating_on = 0
+                cooling_on = 1
+             ENDIF
+
+             q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on))!< [W/GebÃ€udemodell] , observe the directional convention in PALM!
+             surf_usm_h%waste_heat(m) = q_waste_heat
              
           ENDDO !< Horizontal surfaces loop
@@ -1156,4 +1263,6 @@
              kk = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
 !
+!--          (SOME OF THE FOLLOWING (not time-dependent COULD PROBABLY GO INTO A FUNCTION
+!--          EXCEPT facade_element_area, EVERYTHING IS CALCULATED EQUALLY)
 !--          Building geometries  --> not time-dependent
              IF ( l == 0  .OR. l == 1 ) facade_element_area = dx * dzw(kk)  !< [m2] surface area per facade element
@@ -1162,24 +1271,29 @@
              indoor_volume_per_facade = buildings(nb)%vpf(kk)               !< [m3] indoor air volume per facade element           
              window_area_per_facade   = surf_usm_v(l)%frac(ind_wat_win,m)  * facade_element_area  !< [m2] window area per facade element
-             eff_mass_area            = factor_a * floor_area_per_facade    !< [m2] standard values according to Table 12 section 12.3.1.2  (calculate over Eq. (65) according to section 12.3.1.2)
-             c_m                      = factor_c * floor_area_per_facade    !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
-             total_area               = lambda_at * floor_area_per_facade   !< [m2] area of all surfaces pointing to zone Eq. (9) according to section 7.2.2.2
+             
+             building_height          = buildings(nb)%kb_max * dzw(kk)
+             f_cei                    = building_height/(buildings(nb)%height_storey-buildings(nb)%height_cei_con) !< [-] factor for ceiling redcution
+             ngs                      = buildings(nb)%vpf(kk)/f_cei                    !< [m2] calculation of netto ground surface
+             f_sr                     = ngs/floor_area_per_facade                      !< [-] factor for surface reduction
+             eff_mass_area            = buildings(nb)%factor_a * ngs    !< [m2] standard values according to Table 12 section 12.3.1.2  (calculate over Eq. (65) according to section 12.3.1.2)
+             c_m                      = buildings(nb)%factor_c * ngs    !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
+             total_area               = buildings(nb)%lambda_at * floor_area_per_facade   !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
 !
 !--          Calculation of heat transfer coefficient for transmission --> not time-dependent
-             h_tr_w   = window_area_per_facade * u_value_win   !< [W/K] only for windows
+             h_tr_w   = window_area_per_facade * buildings(nb)%u_value_win   !< [W/K] only for windows
              h_tr_is  = total_area * h_is                      !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
              h_tr_ms  = eff_mass_area * h_ms                   !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
-             h_tr_op  = 1 / ( 1 / ( ( facade_element_area - window_area_per_facade ) &
-                                    * lambda_layer3 / s_layer3 * 0.5 ) + 1 / h_tr_ms )
-             h_tr_em  = 1 / ( 1 / h_tr_op - 1 / h_tr_ms )      !< [W/K] Eq. (63), Section 12.2.2
+             h_tr_op  = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade ) &
+                                    * buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp ) + 1.0_wp / h_tr_ms )
+             h_tr_em  = 1.0_wp / ( 1.0_wp / h_tr_op - 1.0_wp / h_tr_ms )      !< [W/K] Eq. (63), Section 12.2.2
 !
 !--          internal air loads dependent on the occupacy of the room
 !--          basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
-             phi_ia = 0.5 * ( ( qint_high * schedule_d + qint_low )            &
-                              * floor_area_per_facade )                      !< [W] Eq. (C.1)
+             phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )            &
+                              * ngs )                      !< [W] Eq. (C.1)
 !
 !--          Airflow dependent on the occupacy of the room
 !--          basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
-             air_change = ( air_change_high * schedule_d + air_change_low )  
+             air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low )  
 !
 !--          Heat transfer of ventilation 
@@ -1187,10 +1301,10 @@
 !--          with heat capacity of air 0.33 Wh/m2K
              h_ve   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *      &
-                                    0.33_wp * (1 - eta_ve ) ) )                    !< [W/K] from ISO 13789 Eq.(10)
+                                    0.33_wp * (1 - buildings(nb)%eta_ve ) ) )                    !< [W/K] from ISO 13789 Eq.(10)
                                     
 !--          Heat transfer coefficient auxiliary variables
-             h_tr_1 = 1 / ( ( 1 / h_ve )   + ( 1 / h_tr_is ) )                  !< [W/K] Eq. (C.6)
+             h_tr_1 = 1.0_wp / ( ( 1.0_wp / h_ve )   + ( 1.0_wp / h_tr_is ) )                  !< [W/K] Eq. (C.6)
              h_tr_2 = h_tr_1 + h_tr_w                                           !< [W/K] Eq. (C.7)
-             h_tr_3 = 1 / ( ( 1 / h_tr_2 ) + ( 1 / h_tr_ms ) )                  !< [W/K] Eq. (C.8)
+             h_tr_3 = 1.0_wp / ( ( 1.0_wp / h_tr_2 ) + ( 1.0_wp / h_tr_ms ) )                  !< [W/K] Eq. (C.8)
 !
 !--          Net short-wave radiation through window area (was i_global)
@@ -1220,6 +1334,6 @@
 !--          DIN 4108 - 2 chap.8
              phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off               &
-                         + window_area_per_facade * net_sw_in * f_c_win * solar_protection_on )    &
-                       * g_value_win * ( 1 - params_f_f ) * params_f_w
+                         + window_area_per_facade * net_sw_in * buildings(nb)%f_c_win * solar_protection_on )    &
+                       * buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w
 !
 !--          Calculation of the mass specific thermal load for internal and external heatsources
@@ -1227,5 +1341,5 @@
 !
 !--          Calculation mass specific thermal load implied non thermal mass 
-             phi_st  =   ( 1 - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1 * total_area ) ) ) &
+             phi_st  =   ( 1.0_wp - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1_wp * total_area ) ) ) &
                        * ( phi_ia + phi_sol )                                       !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
 !
@@ -1233,16 +1347,15 @@
 !--          Step 1: Indoor temperature without heating and cooling
 !--          section C.4.1 Picture C.2 zone 3)
-             phi_hc_nd = 0
-             
+             phi_hc_nd = 0.0_wp
              CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
                                          near_facade_temperature, phi_hc_nd )
 !
 !--          If air temperature between border temperatures of heating and cooling, assign output variable, then ready  
-             IF ( theta_int_h_set <= theta_air  .AND.  theta_air <= theta_int_c_set )  THEN
-                phi_hc_nd_ac = 0
+             IF ( buildings(nb)%theta_int_h_set <= theta_air  .AND.  theta_air <= buildings(nb)%theta_int_c_set )  THEN
+                phi_hc_nd_ac = 0.0_wp
                 phi_hc_nd    = phi_hc_nd_ac
                 theta_air_ac = theta_air
 !
-!--          Step 2: Else, apply 10 W/mÂ² heating/cooling power and calculate indoor temperature
+!--          Step 2: Else, apply 10 W/m2 heating/cooling power and calculate indoor temperature
 !--          again.
              ELSE
@@ -1252,8 +1365,8 @@
 
 !--             Heating or cooling?
-                IF ( theta_air > theta_int_c_set )  THEN
-                   theta_air_set = theta_int_c_set
+                IF ( theta_air > buildings(nb)%theta_int_c_set )  THEN
+                   theta_air_set = buildings(nb)%theta_int_c_set
                 ELSE 
-                   theta_air_set = theta_int_h_set 
+                   theta_air_set = buildings(nb)%theta_int_h_set 
                 ENDIF
 
@@ -1274,5 +1387,5 @@
 !--             If necessary, limit the power to maximum power
 !--             section C.4.1 Picture C.2 zone 2) and 4)
-                IF ( phi_c_max < phi_hc_nd_un  .AND.  phi_hc_nd_un < phi_h_max )  THEN
+                IF ( buildings(nb)%phi_c_max < phi_hc_nd_un  .AND.  phi_hc_nd_un < buildings(nb)%phi_h_max )  THEN
                    phi_hc_nd_ac = phi_hc_nd_un
                    phi_hc_nd = phi_hc_nd_un
@@ -1280,8 +1393,8 @@
 !--             Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
 !--             section C.4.1 Picture C.2 zone 1) and 5)
-                   IF ( phi_hc_nd_un > 0 )  THEN
-                      phi_hc_nd_ac = phi_h_max                                         !< Limit heating
+                   IF ( phi_hc_nd_un > 0.0_wp )  THEN
+                      phi_hc_nd_ac = buildings(nb)%phi_h_max                                         !< Limit heating
                    ELSE 
-                      phi_hc_nd_ac = phi_c_max                                         !< Limit cooling
+                      phi_hc_nd_ac = buildings(nb)%phi_c_max                                         !< Limit cooling
                    ENDIF
                 ENDIF
@@ -1302,5 +1415,5 @@
 !--          Calculate the operating temperature with weighted mean of temperature of air and mean
 !--          Will be used for thermal comfort calculations 
-             theta_op     = 0.3 * theta_air_ac + 0.7 * theta_s
+             theta_op     = 0.3_wp * theta_air_ac + 0.7_wp * theta_s
 !
 !--          Heat flux into the wall. Value needed in urban_surface_mod to 
@@ -1322,15 +1435,15 @@
 !--          Calculation of waste heat
 !--          Anthropogenic heat output
-              IF ( phi_hc_nd_ac > 0 )  THEN 
-                   heating_on = 1
-                   cooling_on = 0
-              ELSE 
-                   heating_on = 0
-                   cooling_on = 1
-              ENDIF
-
-             q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on))   !< [W/m2] , anthropogenic heat output
-!              surf_usm_v(l)%waste_heat(m)=q_waste_heat
-             
+             IF ( phi_hc_nd_ac > 0.0_wp )  THEN 
+                heating_on = 1
+                cooling_on = 0
+             ELSE 
+                heating_on = 0
+                cooling_on = 1
+             ENDIF
+
+             q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on))!< [W/GebÃ€udemodell] , observe the directional convention in PALM!
+             surf_usm_v(l)%waste_heat(m) = q_waste_heat
+
           ENDDO !< Vertical surfaces loop
 
@@ -1339,6 +1452,5 @@
 
 !
-!-- Determine total number of facade elements per building and assign number to 
-!-- building data type.
+!-- Determine the mean building temperature. 
     DO  nb = 1, num_build
 !
@@ -1383,4 +1495,216 @@
  END SUBROUTINE im_main_heatcool
 
+!-----------------------------------------------------------------------------!
+! Description:
+!-------------
+!> Check data output for plant canopy model
+!-----------------------------------------------------------------------------! 
+ SUBROUTINE im_check_data_output( var, unit )
+        
+    IMPLICIT NONE
+    
+    CHARACTER (LEN=*) ::  unit   !<
+    CHARACTER (LEN=*) ::  var    !<
+        
+    SELECT CASE ( TRIM( var ) )
+    
+    
+        CASE ( 'im_hf_roof')
+           unit = 'W m-2'
+        
+        CASE ( 'im_hf_wall_win' )
+           unit = 'W m-2'
+           
+        CASE ( 'im_hf_wall_win_waste' )
+           unit = 'W m-2'
+           
+        CASE ( 'im_hf_roof_waste' )
+           unit = 'W m-2'
+        
+        CASE ( 'im_t_indoor' )
+           unit = 'K'
+        
+        CASE DEFAULT
+           unit = 'illegal'
+           
+    END SELECT
+    
+ END SUBROUTINE
+
+
+!-----------------------------------------------------------------------------!
+! Description:
+!-------------
+!> Check parameters routine for plant canopy model
+!-----------------------------------------------------------------------------! 
+ SUBROUTINE im_check_parameters
+
+!!!!   USE control_parameters,
+!!!!       ONLY: message_string
+       
+   IMPLICIT NONE
+   
+ END SUBROUTINE im_check_parameters
+
+!-----------------------------------------------------------------------------!
+! Description:
+!-------------
+!> Subroutine defining appropriate grid for netcdf variables.
+!> It is called from subroutine netcdf.
+!-----------------------------------------------------------------------------! 
+ SUBROUTINE im_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+
+   IMPLICIT NONE
+   
+   CHARACTER (LEN=*), INTENT(IN)  ::  var
+   LOGICAL, INTENT(OUT)           ::  found
+   CHARACTER (LEN=*), INTENT(OUT) ::  grid_x
+   CHARACTER (LEN=*), INTENT(OUT) ::  grid_y
+   CHARACTER (LEN=*), INTENT(OUT) ::  grid_z
+   
+   found   = .TRUE.
+   
+!
+!-- Check for the grid
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'im_hf_roof', 'im_hf_roof_waste' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zw'
+!
+!--    Heat fluxes at vertical walls are actually defined on stagged grid, i.e. xu, yv.
+       CASE ( 'im_hf_wall_win', 'im_hf_wall_win_waste' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zu'
+          
+       CASE ( 'im_t_indoor' )
+          grid_x = 'x'
+          grid_y = 'y'
+          grid_z = 'zw'
+          
+       CASE DEFAULT
+          found  = .FALSE.
+          grid_x = 'none'
+          grid_y = 'none'
+          grid_z = 'none'
+    END SELECT
+    
+ END SUBROUTINE im_define_netcdf_grid
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE im_data_output_3d( av, variable, found, local_pf, fill_value,      &
+                               nzb_do, nzt_do )
+
+   USE indices
+   
+   USE kinds
+          
+   IMPLICIT NONE
+   
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av    !< 
+    INTEGER(iwp) ::  i     !< 
+    INTEGER(iwp) ::  j     !< 
+    INTEGER(iwp) ::  k     !< 
+    INTEGER(iwp) ::  l     !<
+    INTEGER(iwp) ::  m     !< 
+    INTEGER(iwp) ::  nb    !< index of the building in the building data structure 
+    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
+    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
+
+    
+    LOGICAL      ::  found !< 
+
+    REAL(wp), INTENT(IN) ::  fill_value !< value for the _FillValue attribute
+
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
+    
+    local_pf = fill_value
+    
+    found = .TRUE.
+    
+    SELECT CASE ( TRIM( variable ) )
+!
+!--     Output of indoor temperature. All grid points within the building are
+!--     filled with values, while atmospheric grid points are set to _FillValues.
+        CASE ( 'im_t_indoor' )
+           IF ( av == 0 ) THEN
+              DO  i = nxl, nxr
+                 DO  j = nys, nyn
+                    IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
+!
+!--                    Determine index of the building within the building data structure. 
+                       nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ),   &
+                                    DIM = 1 )
+!
+!--                    Write mean building temperature onto output array. Please note,
+!--                    in contrast to many other loops in the output, the vertical 
+!--                    bounds are determined by the lowest and hightest vertical index
+!--                    occupied by the building. 
+                       DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
+                          local_pf(i,j,k) = buildings(nb)%t_in(k)
+                       ENDDO
+                    ENDIF
+                 ENDDO
+              ENDDO
+           ENDIF  
+
+        CASE ( 'im_hf_roof' )
+           IF ( av == 0 ) THEN
+              DO  m = 1, surf_usm_h%ns
+                  i = surf_usm_h%i(m) !+ surf_usm_h%ioff
+                  j = surf_usm_h%j(m) !+ surf_usm_h%joff
+                  k = surf_usm_h%k(m) !+ surf_usm_h%koff
+                  local_pf(i,j,k) = surf_usm_h%iwghf_eb(m)
+              ENDDO
+           ENDIF 
+    
+        CASE ( 'im_hf_roof_waste' )
+           IF ( av == 0 ) THEN
+              DO m = 1, surf_usm_h%ns 
+                 i = surf_usm_h%i(m) !+ surf_usm_h%ioff
+                 j = surf_usm_h%j(m) !+ surf_usm_h%joff
+                 k = surf_usm_h%k(m) !+ surf_usm_h%koff
+                 local_pf(i,j,k) = surf_usm_h%waste_heat(m)
+              ENDDO
+           ENDIF                     
+       
+       CASE ( 'im_hf_wall_win' )
+           IF ( av == 0 ) THEN
+              DO l = 0, 3
+                 DO m = 1, surf_usm_v(l)%ns
+                    i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
+                    j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
+                    k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
+                    local_pf(i,j,k) = surf_usm_v(l)%iwghf_eb(m)
+                 ENDDO
+              ENDDO
+           ENDIF
+           
+        CASE ( 'im_hf_wall_win_waste' )
+           IF ( av == 0 ) THEN
+              DO l = 0, 3
+                 DO m = 1, surf_usm_v(l)%ns 
+                    i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
+                    j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
+                    k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
+                    local_pf(i,j,k) =  surf_usm_v(l)%waste_heat(m) 
+                 ENDDO
+              ENDDO
+           ENDIF      
+    
+        CASE DEFAULT
+           found = .FALSE.
+           
+    END SELECT    
+
+ END SUBROUTINE im_data_output_3d          
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1397,10 +1721,5 @@
     CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
 
-    
-    
-    NAMELIST /indoor_parameters/  building_type, dt_indoor,                           &
-                                  initial_indoor_temperature
-
-!    line = ' '
+    NAMELIST /indoor_parameters/  dt_indoor, initial_indoor_temperature
 
 !
@@ -1410,5 +1729,4 @@
     DO   WHILE ( INDEX( line, '&indoor_parameters' ) == 0 )
        READ ( 11, '(A)', END=10 )  line
-!    PRINT*, 'line: ', line
     ENDDO
     BACKSPACE ( 11 )
Index: /palm/trunk/SOURCE/module_interface.f90
===================================================================
--- /palm/trunk/SOURCE/module_interface.f90	(revision 3743)
+++ /palm/trunk/SOURCE/module_interface.f90	(revision 3744)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Add indoor model
 ! 
 ! Former revisions:
@@ -174,4 +174,7 @@
    USE indoor_model_mod,                                                       &
        ONLY:  im_parin,                                                        &
+              im_check_data_output,                                            &
+              im_check_parameters,                                             &
+              im_data_output_3d,                                               &
               im_init
 
@@ -497,4 +500,5 @@
    IF ( air_chemistry )        CALL chem_check_parameters
    IF ( gust_module_enabled )  CALL gust_check_parameters
+   IF ( indoor_model )         CALL im_check_parameters
    IF ( land_surface )         CALL lsm_check_parameters
    IF ( large_scale_forcing  .OR.  nudging )  CALL lsf_nudging_check_parameters ! ToDo: create single module switch
@@ -647,4 +651,8 @@
    ENDIF
 
+   IF ( unit == 'illegal' .AND. indoor_model ) THEN
+       CALL im_check_data_output( variable, unit )
+   ENDIF
+   
    IF ( unit == 'illegal'  .AND.  urban_surface                      &
         .AND.  variable(1:4) == 'usm_' )  THEN  ! ToDo: remove aditional conditions
@@ -1022,4 +1030,9 @@
    IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
       CALL gust_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
+      resorted = .TRUE.
+   ENDIF
+   
+   IF ( .NOT. found  .AND.  indoor_model )  THEN
+      CALL im_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
       resorted = .TRUE.
    ENDIF
Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3744)
@@ -20,5 +20,4 @@
 ! Current revisions:
 ! -----------------
-! 
 ! 
 ! Former revisions:
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3744)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Add indoor model
 ! 
 ! Former revisions:
@@ -660,6 +660,6 @@
                do2d_yz_time_count, dt_data_output_av, dt_do2d_xy, dt_do2d_xz,  &
                dt_do2d_yz, dt_do3d, dt_write_agent_data, mask_size,            &
-               do2d_xy_time_count,                                             &
-               do3d_time_count, domask_time_count, end_time, land_surface,     &
+               do2d_xy_time_count, do3d_time_count, domask_time_count,         &
+               end_time, indoor_model, land_surface,                           &
                mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global,      &
                mask_k_global, mask_surface,                                    &
@@ -683,4 +683,7 @@
     USE kinds
 
+    USE indoor_model_mod,                                                      &
+        ONLY: im_define_netcdf_grid
+    
     USE land_surface_model_mod,                                                &
         ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
@@ -1911,4 +1914,10 @@
                       CALL gust_define_netcdf_grid( do3d(av,i), found, grid_x, &
                                                     grid_y, grid_z )
+                   ENDIF
+!
+!--                Check for indoor model quantities                   
+                   IF ( .NOT. found .AND. indoor_model ) THEN 
+                      CALL im_define_netcdf_grid( do3d(av,i), found,           &
+                                                  grid_x, grid_x, grid_z )
                    ENDIF
 
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3744)
@@ -22,5 +22,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix in transpiration, floating invalid when temperature 
+! becomes > 40 degrees
 ! 
 ! Former revisions:
@@ -236,4 +237,6 @@
 !> 2) Calculation of sinks and sources of momentum, heat and scalar concentration 
 !> due to canopy elements (subroutine pcm_tendency).
+!
+! @todo - precalculate constant terms in pcm_calc_transpiration_rate
 !------------------------------------------------------------------------------!
  MODULE plant_canopy_model_mod
@@ -422,7 +425,7 @@
      sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
 !--  Wind speed
-     wind_speed = SQRT( ((u(k,j,i) + u(k,j,i+1))/2.0_wp )**2 +           &
-                        ((v(k,j,i) + v(k,j,i+1))/2.0_wp )**2 +           &
-                        ((w(k,j,i) + w(k,j,i+1))/2.0_wp )**2 )
+     wind_speed = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 +            &
+                        ( 0.5_wp * ( v(k,j,i) + v(k,j,i+1) ) )**2 +            &
+                        ( 0.5_wp * ( w(k,j,i) + w(k,j,i+1) ) )**2 )
 !--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
      g_b = 0.01_wp * wind_speed + 0.0071_wp
@@ -434,5 +437,5 @@
 !--  Second function for calculation of stomatal conductivity (temperature dependency)
 !--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
-     f2 = MAX(t2_min, (temp-t2_min) * (t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
+     f2 = MAX(t2_min, (temp-t2_min) * MAX(0.0_wp,t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
               ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
 !--  Water pressure deficit
Index: /palm/trunk/SOURCE/surface_data_output_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/surface_data_output_mod.f90	(revision 3744)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Output of waste_heat and innermost wall flux from indoor model
 ! 
 ! Former revisions:
@@ -28,5 +28,5 @@
 ! clean-up coding layout
 ! 
-! 3735 2019-02-12 09:52:40Z dom_dwd_user
+! 3735 2019-02-12 09:52:40Z suehring
 ! - Split initialization into initialization of arrays and further initialization
 !   in order to enable reading of restart data.
@@ -76,5 +76,5 @@
 ! Authors:
 ! --------
-! @author Klaus Ketelsen, Matthias Suehring
+! @author Klaus Ketelsen, Matthias Suehring, Tobias Gronemeier
 !
 ! Description:
@@ -2778,5 +2778,5 @@
 
                ENDIF
-
+               
             CASE ( 'uvw1' )
 !
@@ -2807,5 +2807,64 @@
 
                ENDIF
-
+!
+!--         Waste heat from indoor model
+            CASE ( 'waste_heat' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_data_output_collect( surf_def_h(0)%waste_heat,  &
+                                               surf_def_h(1)%waste_heat,       &
+                                               surf_lsm_h%waste_heat,          &
+                                               surf_usm_h%waste_heat,          &
+                                               surf_def_v(0)%waste_heat,       &
+                                               surf_lsm_v(0)%waste_heat,       &
+                                               surf_usm_v(0)%waste_heat,       &
+                                               surf_def_v(1)%waste_heat,       &
+                                               surf_lsm_v(1)%waste_heat,       &
+                                               surf_usm_v(1)%waste_heat,       &
+                                               surf_def_v(2)%waste_heat,       &
+                                               surf_lsm_v(2)%waste_heat,       &
+                                               surf_usm_v(2)%waste_heat,       &
+                                               surf_def_v(3)%waste_heat,       &
+                                               surf_lsm_v(3)%waste_heat,       &
+                                               surf_usm_v(3)%waste_heat )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
+!
+!--         Innermost building wall flux from indoor model
+            CASE ( 'im_hf' )
+!
+!--            Output of instantaneous data
+               IF ( av == 0 )  THEN
+                  CALL surface_data_output_collect( surf_def_h(0)%iwghf_eb,    &
+                                               surf_def_h(1)%iwghf_eb,         &
+                                               surf_lsm_h%iwghf_eb,            &
+                                               surf_usm_h%iwghf_eb,            &
+                                               surf_def_v(0)%iwghf_eb,         &
+                                               surf_lsm_v(0)%iwghf_eb,         &
+                                               surf_usm_v(0)%iwghf_eb,         &
+                                               surf_def_v(1)%iwghf_eb,         &
+                                               surf_lsm_v(1)%iwghf_eb,         &
+                                               surf_usm_v(1)%iwghf_eb,         &
+                                               surf_def_v(2)%iwghf_eb,         &
+                                               surf_lsm_v(2)%iwghf_eb,         &
+                                               surf_usm_v(2)%iwghf_eb,         &
+                                               surf_def_v(3)%iwghf_eb,         &
+                                               surf_lsm_v(3)%iwghf_eb,         &
+                                               surf_usm_v(3)%iwghf_eb )
+               ELSE
+!
+!--               Output of averaged data
+                  surfaces%var_out(:) = surfaces%var_av(:,n_out) /             &
+                                        REAL( average_count_surf, KIND=wp )
+                  surfaces%var_av(:,n_out) = 0.0_wp
+
+               ENDIF
 !
 !--            Add further variables:
@@ -3706,4 +3765,40 @@
                                            surf_lsm_v(3)%uvw_abs,              &
                                            surf_usm_v(3)%uvw_abs, n_out )
+                                           
+            CASE ( 'waste_heat' )
+               CALL surface_data_output_sum_up( surf_def_h(0)%waste_heat,      &
+                                           surf_def_h(1)%waste_heat,           &
+                                           surf_lsm_h%waste_heat,              &
+                                           surf_usm_h%waste_heat,              &
+                                           surf_def_v(0)%waste_heat,           &
+                                           surf_lsm_v(0)%waste_heat,           &
+                                           surf_usm_v(0)%waste_heat,           &
+                                           surf_def_v(1)%waste_heat,           &
+                                           surf_lsm_v(1)%waste_heat,           &
+                                           surf_usm_v(1)%waste_heat,           &
+                                           surf_def_v(2)%waste_heat,           &
+                                           surf_lsm_v(2)%waste_heat,           &
+                                           surf_usm_v(2)%waste_heat,           &
+                                           surf_def_v(3)%waste_heat,           &
+                                           surf_lsm_v(3)%waste_heat,           &
+                                           surf_usm_v(3)%waste_heat, n_out )
+                                           
+            CASE ( 'im_hf' )
+               CALL surface_data_output_sum_up( surf_def_h(0)%iwghf_eb,        &
+                                           surf_def_h(1)%iwghf_eb,             &
+                                           surf_lsm_h%iwghf_eb,                &
+                                           surf_usm_h%iwghf_eb,                &
+                                           surf_def_v(0)%iwghf_eb,             &
+                                           surf_lsm_v(0)%iwghf_eb,             &
+                                           surf_usm_v(0)%iwghf_eb,             &
+                                           surf_def_v(1)%iwghf_eb,             &
+                                           surf_lsm_v(1)%iwghf_eb,             &
+                                           surf_usm_v(1)%iwghf_eb,             &
+                                           surf_def_v(2)%iwghf_eb,             &
+                                           surf_lsm_v(2)%iwghf_eb,             &
+                                           surf_usm_v(2)%iwghf_eb,             &
+                                           surf_def_v(3)%iwghf_eb,             &
+                                           surf_lsm_v(3)%iwghf_eb,             &
+                                           surf_usm_v(3)%iwghf_eb, n_out )
 
          END SELECT
@@ -4161,6 +4256,6 @@
 
        USE control_parameters,                                                 &
-           ONLY:  averaging_interval, dt_data_output, initializing_actions,    &
-                  message_string
+           ONLY:  averaging_interval, dt_data_output, indoor_model,            &
+                  initializing_actions, message_string
 
        USE pegrid,                                                             &
@@ -4293,4 +4388,14 @@
             CASE ( 'r_a', 'r_canopy', 'r_soil', 'r_s' )
                unit = 's/m'
+               
+            CASE ( 'waste_heat', 'im_hf' )
+               IF ( .NOT. indoor_model )  THEN
+                  message_string = TRIM( trimvar ) //                          &
+                             ' requires the indoor model'
+               CALL message( 'surface_data_output_check_parameters',           &
+                             'PA0588', 1, 2, 0, 6, 0 )
+               ENDIF
+            
+               unit = 'W/m2'
 
             CASE DEFAULT
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 3744)
@@ -21,5 +21,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! Bugfix, missing calculation of 10cm temperature at vertical building walls, 
+! required for indoor model
 ! 
 ! Former revisions:
@@ -443,7 +444,7 @@
              CALL calc_pt_near_surface ( '2m' )
           ENDIF
-          IF ( indoor_model )  THEN
-             CALL calc_pt_near_surface ( '10cm' )
-          ENDIF
+!
+!--       Calculate 10cm temperature, required in indoor model
+          IF ( indoor_model )  CALL calc_pt_near_surface ( '10cm' )
        ENDIF
 
@@ -548,4 +549,7 @@
 !--          Compute respective surface fluxes for momentum and TKE
              CALL calc_surface_fluxes
+!
+!--          Calculate 10cm temperature, required in indoor model
+             IF ( indoor_model )  CALL calc_pt_near_surface ( '10cm' )
           ENDIF
        ENDDO
@@ -566,4 +570,5 @@
 !--          Compute respective surface fluxes for momentum and TKE
              CALL calc_surface_fluxes
+             
           ENDIF
        ENDDO
@@ -625,4 +630,7 @@
 !--          Compute respective surface fluxes for momentum and TKE
              CALL calc_surface_fluxes
+!
+!--          Calculate 10cm temperature, required in indoor model             
+             IF ( indoor_model )  CALL calc_pt_near_surface ( '10cm' )
           ENDIF
        ENDDO
@@ -2020,9 +2028,9 @@
                 k = surf%k(m)
 
-                surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa  &
-                                   * ( log( 0.1_wp /  surf%z0h(m) )                    &
-                                  - psi_h( 0.1_wp / surf%ol(m) )                    &
+                surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa      &
+                                   * ( LOG( 0.1_wp /  surf%z0h(m) )            &
+                                  - psi_h( 0.1_wp / surf%ol(m) )               &
                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
-
+                                    
              ENDDO
 
@@ -2036,7 +2044,7 @@
                 k = surf%k(m)
 
-                surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa  &
-                                   * ( log( 2.0_wp /  surf%z0h(m) )                    &
-                                     - psi_h( 2.0_wp / surf%ol(m) )                    &
+                surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa        &
+                                   * ( LOG( 2.0_wp /  surf%z0h(m) )            &
+                                     - psi_h( 2.0_wp / surf%ol(m) )            &
                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
 
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 3744)
@@ -21,5 +21,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +waste_heat
 ! 
 ! Former revisions:
@@ -455,4 +455,7 @@
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rho_c_total       !< volumetric heat capacity of the actual soil matrix (J/m3/K)
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  root_fr           !< root fraction within the soil layers
+       
+!--    Indoor model variables
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  waste_heat          !< waste heat 
 !
 !--    Urban surface variables
@@ -1243,5 +1246,5 @@
 !--    Characteristic temperature and surface flux of sensible heat
        DEALLOCATE ( surfaces%ts )    
-       DEALLOCATE ( surfaces%shf )    
+       DEALLOCATE ( surfaces%shf )
 !
 !--    surface temperature
@@ -1371,5 +1374,5 @@
 !--    Characteristic temperature and surface flux of sensible heat
        ALLOCATE ( surfaces%ts(1:surfaces%ns)  )    
-       ALLOCATE ( surfaces%shf(1:surfaces%ns) )    
+       ALLOCATE ( surfaces%shf(1:surfaces%ns) )
 !
 !--    Surface temperature
@@ -1537,5 +1540,5 @@
 !
 !--    Sensible heat flux
-       DEALLOCATE ( surfaces%shf )    
+       DEALLOCATE ( surfaces%shf )
 !
 !--    Latent heat flux
@@ -1609,5 +1612,5 @@
 !
 !--    Sensible heat flux
-       ALLOCATE ( surfaces%shf(1:surfaces%ns) )    
+       ALLOCATE ( surfaces%shf(1:surfaces%ns) )
 !
 !--    Latent heat flux
@@ -3637,5 +3640,5 @@
              CALL wrd_write_string( 'surf_h(' // dum // ')%sasws' ) 
              WRITE ( 14 )  surf_h(l)%sasws
-          ENDIF     
+          ENDIF
   
        ENDDO
@@ -3781,5 +3784,5 @@
              WRITE ( 14 )  surf_v(l)%mom_flux_tke
           ENDIF 
-
+         
        ENDDO
 
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3743)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3744)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Call indoor model after first timestep
 ! 
 ! Former revisions:
@@ -1413,10 +1413,13 @@
 !--    If required, calculate indoor temperature, waste heat, heat flux
 !--    through wall, etc.
-!--    dt_indoor steers the frequency of the indoor model calculations
+!--    dt_indoor steers the frequency of the indoor model calculations.
+!--    Note, at first timestep indoor model is called, in order to provide
+!--    a waste heat flux. 
        IF ( indoor_model )  THEN
 
           time_indoor = time_indoor + dt_3d
 
-          IF ( time_indoor >= dt_indoor )  THEN
+          IF ( time_indoor >= dt_indoor  .OR.  current_timestep_number == 0 )  &
+          THEN
 
              time_indoor = time_indoor - dt_indoor
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3743)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3744)
@@ -23,5 +23,8 @@
 ! Current revisions:
 ! ------------------
-! 
+! - Remove internal flag indoor_model (is a global control parameter)
+! - add waste heat from buildings to the kinmatic heat flux
+! - consider waste heat in restart data
+! - remove unused USE statements
 ! 
 ! Former revisions:
@@ -401,10 +404,8 @@
 !> @todo Output of _av variables in case of restarts
 !> @todo Revise flux conversion in energy-balance solver
-!> @todo Bugfixing in nopointer branch 
 !> @todo Check optimizations for RMA operations
 !> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
 !> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
 !>       factor 3 between min and max time
-!> @todo Move setting of flag indoor_model to indoor_model_mod once available
 !> @todo Check divisions in wtend (etc.) calculations for possible division
 !>       by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
@@ -424,5 +425,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  coupling_start_time, topography, dt_3d, humidity,               &
+        ONLY:  coupling_start_time, topography, dt_3d, humidity, indoor_model, &
                intermediate_timestep_count, initializing_actions,              &
                intermediate_timestep_count_max, simulated_time, end_time,      &
@@ -431,6 +432,5 @@
                pt_surface, large_scale_forcing, lsf_surf, spinup,              &
                spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
-               average_count_3d, varnamelength, urban_surface,                 &
-               plant_canopy, dz
+               average_count_3d, varnamelength, urban_surface, dz
 
     USE bulk_cloud_model_mod,                                                  &
@@ -455,8 +455,5 @@
               
     USE pegrid
-    
-    USE plant_canopy_model_mod,                                                &
-        ONLY:  pc_heating_rate, pc_transpiration_rate, pc_latent_rate
-    
+        
     USE radiation_model_mod,                                                   &
         ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,        &
@@ -535,5 +532,4 @@
                                                    !< (e.g.transportation) are used
     LOGICAL ::  force_radiation_call_l = .FALSE.   !< flag parameter for unscheduled radiation model calls
-    LOGICAL ::  indoor_model = .FALSE.             !< whether to use the indoor model
     LOGICAL ::  read_wall_temp_3d = .FALSE.
     LOGICAL ::  usm_wall_mod = .FALSE.             !< reduces conductivity of the first 2 wall layers by factor 0.1
@@ -1257,5 +1253,16 @@
            ALLOCATE ( surf_usm_v(l)%target_temp_summer(1:surf_usm_v(l)%ns) )
            ALLOCATE ( surf_usm_v(l)%target_temp_winter(1:surf_usm_v(l)%ns) )
-        ENDDO   
+        ENDDO
+!
+!--     In case the indoor model is applied, allocate memory for waste heat 
+!--     and indoor temperature.
+        IF ( indoor_model )  THEN
+           ALLOCATE ( surf_usm_h%waste_heat(1:surf_usm_h%ns) )
+           surf_usm_h%waste_heat = 0.0_wp
+           DO  l = 0, 3
+              ALLOCATE ( surf_usm_v(l)%waste_heat(1:surf_usm_v(l)%ns) )
+              surf_usm_v(l)%waste_heat = 0.0_wp
+           ENDDO
+        ENDIF
 !
 !--     Allocate flag indicating ground floor level surface elements
@@ -6054,5 +6061,4 @@
                            usm_material_model,                                 &
                            wall_category,                                      &
-                           indoor_model,                                       &
                            wall_inner_temperature,                             &
                            roof_inner_temperature,                             &
@@ -6073,5 +6079,4 @@
                            usm_material_model,                                 &
                            wall_category,                                      &
-                           indoor_model,                                       &
                            wall_inner_temperature,                             &
                            roof_inner_temperature,                             &
@@ -6264,10 +6269,22 @@
 
        LOGICAL, INTENT(OUT)  ::  found 
+!!!    suehring: Why the SAVE attribute?       
+       REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE   ::  tmp_surf_wall_h
+       REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE   ::  tmp_surf_window_h
+       REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE   ::  tmp_surf_green_h
+       REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE   ::  tmp_surf_waste_h
        
-       REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE   ::  tmp_surf_wall_h, tmp_surf_window_h, tmp_surf_green_h
-       REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  tmp_wall_h, tmp_window_h, tmp_green_h
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  tmp_wall_h
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  tmp_window_h
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  tmp_green_h
        
-       TYPE( t_surf_vertical ), DIMENSION(0:3), SAVE ::  tmp_surf_wall_v, tmp_surf_window_v, tmp_surf_green_v
-       TYPE( t_wall_vertical ), DIMENSION(0:3), SAVE ::  tmp_wall_v, tmp_window_v, tmp_green_v
+       TYPE( t_surf_vertical ), DIMENSION(0:3), SAVE ::  tmp_surf_wall_v
+       TYPE( t_surf_vertical ), DIMENSION(0:3), SAVE ::  tmp_surf_window_v
+       TYPE( t_surf_vertical ), DIMENSION(0:3), SAVE ::  tmp_surf_green_v
+       TYPE( t_surf_vertical ), DIMENSION(0:3), SAVE ::  tmp_surf_waste_v
+       
+       TYPE( t_wall_vertical ), DIMENSION(0:3), SAVE ::  tmp_wall_v
+       TYPE( t_wall_vertical ), DIMENSION(0:3), SAVE ::  tmp_window_v
+       TYPE( t_wall_vertical ), DIMENSION(0:3), SAVE ::  tmp_green_v
 
 
@@ -6289,4 +6306,6 @@
                       DEALLOCATE( tmp_surf_green_h ) 
                    IF ( ALLOCATED( tmp_green_h) ) DEALLOCATE( tmp_green_h )
+                   IF ( ALLOCATED( tmp_surf_waste_h) )                         &
+                      DEALLOCATE( tmp_surf_waste_h )
   
 !
@@ -6305,4 +6324,5 @@
                    ALLOCATE( tmp_green_h(nzb_wall:nzt_wall+1,                  &
                                          1:ns_h_on_file_usm) )
+                   ALLOCATE( tmp_surf_waste_h(1:ns_h_on_file_usm) )
 
                 ENDIF
@@ -6325,4 +6345,6 @@
                       IF ( ALLOCATED( tmp_green_v(l)%t ) )                     &
                          DEALLOCATE( tmp_green_v(l)%t )
+                      IF ( ALLOCATED( tmp_surf_waste_v(l)%t ) )                &
+                         DEALLOCATE( tmp_surf_waste_v(l)%t )
                    ENDDO 
 
@@ -6343,4 +6365,5 @@
                       ALLOCATE( tmp_green_v(l)%t(nzb_wall:nzt_wall+1,          &
                                                  1:ns_v_on_file_usm(l) ) )
+                      ALLOCATE( tmp_surf_waste_v(l)%t(1:ns_v_on_file_usm(l)) )
                    ENDDO
 
@@ -6614,4 +6637,89 @@
                                         t_surf_window_v_1(3)%t,                &  
                                         tmp_surf_window_v(3)%t,                &
+                                        surf_usm_v(3)%start_index,             & 
+                                        start_index_on_file,                   &
+                                        end_index_on_file,                     &
+                                        nxlc, nysc,                            &
+                                        nxlf, nxrf, nysf, nynf,                &
+                                        nys_on_file, nyn_on_file,              &
+                                        nxl_on_file,nxr_on_file )
+
+             CASE ( 'waste_heat_h' )
+                IF ( k == 1 )  THEN
+                   IF ( .NOT.  ALLOCATED( surf_usm_h%waste_heat ) )            &
+                      ALLOCATE( surf_usm_h%waste_heat(1:surf_usm_h%ns) )
+                   READ ( 13 )  tmp_surf_waste_h
+                ENDIF              
+                CALL surface_restore_elements(                                 &
+                                        surf_usm_h%waste_heat,                 &
+                                        tmp_surf_waste_h,                      &
+                                        surf_usm_h%start_index,                &
+                                        start_index_on_file,                   &
+                                        end_index_on_file,                     &
+                                        nxlc, nysc,                            &
+                                        nxlf, nxrf, nysf, nynf,                &
+                                        nys_on_file, nyn_on_file,              &
+                                        nxl_on_file,nxr_on_file )                 
+                                        
+             CASE ( 'waste_heat_v(0)' )
+                IF ( k == 1 )  THEN
+                   IF ( .NOT.  ALLOCATED( surf_usm_v(0)%waste_heat ) )         &
+                      ALLOCATE( surf_usm_v(0)%waste_heat(1:surf_usm_v(0)%ns) )
+                   READ ( 13 )  tmp_surf_waste_v(0)%t
+                ENDIF
+                CALL surface_restore_elements(                                 &
+                                        surf_usm_v(0)%waste_heat,              &
+                                        tmp_surf_waste_v(0)%t,                 &
+                                        surf_usm_v(0)%start_index,             & 
+                                        start_index_on_file,                   &
+                                        end_index_on_file,                     &
+                                        nxlc, nysc,                            &
+                                        nxlf, nxrf, nysf, nynf,                &
+                                        nys_on_file, nyn_on_file,              &
+                                        nxl_on_file,nxr_on_file )
+                      
+             CASE ( 'waste_heat_v(1)' )
+                IF ( k == 1 )  THEN
+                   IF ( .NOT.  ALLOCATED( surf_usm_v(1)%waste_heat ) )         &
+                      ALLOCATE( surf_usm_v(1)%waste_heat(1:surf_usm_v(1)%ns) )
+                   READ ( 13 )  tmp_surf_waste_v(1)%t
+                ENDIF
+                CALL surface_restore_elements(                                 &
+                                        surf_usm_v(1)%waste_heat,              &
+                                        tmp_surf_waste_v(1)%t,                 &
+                                        surf_usm_v(1)%start_index,             & 
+                                        start_index_on_file,                   &
+                                        end_index_on_file,                     &
+                                        nxlc, nysc,                            &
+                                        nxlf, nxrf, nysf, nynf,                &
+                                        nys_on_file, nyn_on_file,              &
+                                        nxl_on_file,nxr_on_file )
+
+             CASE ( 'waste_heat_v(2)' )
+                IF ( k == 1 )  THEN
+                   IF ( .NOT.  ALLOCATED( surf_usm_v(2)%waste_heat ) )         &
+                      ALLOCATE( surf_usm_v(2)%waste_heat(1:surf_usm_v(2)%ns) )
+                   READ ( 13 )  tmp_surf_waste_v(2)%t
+                ENDIF
+                CALL surface_restore_elements(                                 &
+                                        surf_usm_v(2)%waste_heat,              &
+                                        tmp_surf_waste_v(2)%t,                 &
+                                        surf_usm_v(2)%start_index,             & 
+                                        start_index_on_file,                   &
+                                        end_index_on_file,                     &
+                                        nxlc, nysc,                            &
+                                        nxlf, nxrf, nysf, nynf,                &
+                                        nys_on_file, nyn_on_file,              &
+                                        nxl_on_file,nxr_on_file )
+                      
+             CASE ( 'waste_heat_v(3)' )
+                IF ( k == 1 )  THEN
+                   IF ( .NOT.  ALLOCATED( surf_usm_v(3)%waste_heat ) )         &
+                      ALLOCATE( surf_usm_v(3)%waste_heat(1:surf_usm_v(3)%ns) )
+                   READ ( 13 )  tmp_surf_waste_v(3)%t
+                ENDIF
+                CALL surface_restore_elements(                                 &
+                                        surf_usm_v(3)%waste_heat,              &
+                                        tmp_surf_waste_v(3)%t,                 &
                                         surf_usm_v(3)%start_index,             & 
                                         start_index_on_file,                   &
@@ -7627,5 +7735,5 @@
         j_off = surf_usm_h%joff
         i_off = surf_usm_h%ioff
-
+        
 !        
 !--     First, treat horizontal surface elements
@@ -7957,4 +8065,10 @@
 !--        diffusion_s, surface_layer_fluxes,...
            surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / c_p
+!
+!--        If the indoor model is applied, further add waste heat from buildings to the 
+!--        kinematic flux.
+           IF ( indoor_model )  THEN
+              surf_usm_h%shf(m) = surf_usm_h%shf(m) + surf_usm_h%waste_heat(m) / rho_cp
+           ENDIF
       
 
@@ -8382,4 +8496,11 @@
 !--           diffusion_s, surface_layer_fluxes,...
               surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / c_p
+!
+!--           If the indoor model is applied, further add waste heat from buildings to the 
+!--           kinematic flux.
+              IF ( indoor_model )  THEN
+                 surf_usm_v(l)%shf(m) = surf_usm_v(l)%shf(m) +                       &
+                                        surf_usm_v(l)%waste_heat(m) / rho_cp
+              ENDIF              
 
               IF ( surf_usm_v(l)%frac(ind_pav_green,m) > 0.0_wp ) THEN
@@ -8685,5 +8806,17 @@
         CALL wrd_write_string( 't_surf_green_h' )
         WRITE ( 14 )  t_surf_green_h
- 
+!
+!--     Write restart data which is especially needed for the urban-surface 
+!--     model. In order to do not fill up the restart routines in 
+!--     surface_mod. 
+!--     Output of waste heat from indoor model. Restart data is required in 
+!--     this special case, because the indoor model where waste heat is 
+!--     computed is call each hour (current default), so that waste heat would
+!--     have zero value until next call of indoor model. 
+        IF ( indoor_model )  THEN
+           CALL wrd_write_string( 'waste_heat_h' )
+           WRITE ( 14 )  surf_usm_h%waste_heat
+        ENDIF   
+           
         DO  l = 0, 3
  
@@ -8703,5 +8836,10 @@
  
            CALL wrd_write_string( 't_surf_green_v(' // dum // ')' )
-           WRITE ( 14 ) t_surf_green_v(l)%t    
+           WRITE ( 14 ) t_surf_green_v(l)%t  
+           
+           IF ( indoor_model )  THEN
+              CALL wrd_write_string( 'waste_heat_v(' // dum // ')' )
+              WRITE ( 14 )  surf_usm_v(l)%waste_heat
+           ENDIF
            
         ENDDO
@@ -8742,5 +8880,4 @@
         
         ENDDO
- 
         
      END SUBROUTINE usm_wrd_local
