Index: palm/trunk/SOURCE/indoor_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/indoor_model_mod.f90	(revision 3685)
+++ 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 )
