Changeset 4148 for palm/trunk

Timestamp:

Aug 8, 2019 11:26:00 AM (5 years ago)

Author:

suehring

Message:

indoor_model_mod: Revision of indoor model; urban_surface_mod: parameters for waste heat from cooling and heating are introduced to building data base; initialization of building data base moved to an earlier point of time before indoor model will be initialized; radiation_model_mod: minor improvement in some comments; synthetic_turbulence_generator: unused variable removed

Location:

palm/trunk/SOURCE

Files:

• indoor_model_mod.f90 (modified) (47 diffs)
• radiation_model_mod.f90 (modified) (5 diffs)
• synthetic_turbulence_generator_mod.f90 (modified) (2 diffs)
• urban_surface_mod.f90 (modified) (11 diffs)

palm/trunk/SOURCE/indoor_model_mod.f90

@@ -25 +25 @@
 ! Former revisions:
 ! -----------------
-! $Id$
-! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
+! $Id: indoor_model_mod.f90 
+! - change calculation of a_m and c_m
+! - change calculation of u-values (use h_es in building array)
+! - rename h_tr_... to  h_t_...
+!          h_tr_em  to  h_t_wm
+!          h_tr_op  to  h_t_wall
+!          h_tr_w   to  h_t_es
+! - rename h_ve     to  h_v
+! - rename h_is     to  h_ms
+! - inserted net_floor_area
+! - inserted params_waste_heat_h, params_waste_heat_c from building database
+!   in building array
+! - change calculation of q_waste_heat
+! - bugfix in averaging mean indoor temperature 
 ! 
-! 3885 2019-04-11 11:29:34Z kanani
-! Changes related to global restructuring of location messages and introduction 
-! of additional debug messages
 ! 
-! 3786 2019-03-06 16:58:03Z raasch
-! unused variables removed
-! 
-! 3759 2019-02-21 15:53:45Z suehring
+! 3759 2019-02-21 15:53:45Z suehring $
 ! - Calculation of total building volume
 ! - Several bugfixes
@@ -85 +91 @@
 !> @todo Replace window_area_per_facade by %frac(1,m) for window
 !> @todo emissivity change for window blinds if solar_protection_on=1 
-!> @todo write datas in netcdf file as output data
-!> @todo reduce the building volume with netto ground surface to take respect costruction areas like walls and ceilings. Have effect on factor_a, factor_c, airchange and lambda_at
-!>
+
 !> @note Do we allow use of integer flags, or only logical flags? (concerns e.g. cooling_on, heating_on)
 !> @note How to write indoor temperature output to pt array?
@@ -96 +100 @@
 
     USE control_parameters,                                                    &
-        ONLY:  debug_output, initializing_actions
+        ONLY:  initializing_actions
 
     USE kinds
@@ -120 +124 @@
        INTEGER(iwp) ::  num_facades_per_building_v = 0    !< total number of vertical facades elements
        INTEGER(iwp) ::  num_facades_per_building_v_l = 0  !< number of vertical facade elements on local subdomain
+       INTEGER(iwp) ::  ventilation_int_loads             !< [-] allocation of activity in the building
 
        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  l_v            !< index array linking surface-element orientation index
@@ -131 +136 @@
        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) ::  building_height   !< building height
-       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) ::  vol_tot           !< total building volume
+       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) ::  area_facade           !< [m2] area of total facade
+       REAL(wp) ::  building_height       !< building height
+       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) ::  f_c_win               !< [-] shading factor
+       REAL(wp) ::  g_value_win           !< [-] SHGC factor
+       REAL(wp) ::  h_es                  !< [W/(m2 K)] surface-related heat transfer coefficient between extern and surface 
+       REAL(wp) ::  height_cei_con        !< [m] ceiling construction heigth
+       REAL(wp) ::  height_storey         !< [m] storey heigth
+       REAL(wp) ::  params_waste_heat_c   !< [-] anthropogenic heat outputs for cooling e.g. 1.33 for KKM with COP = 3
+       REAL(wp) ::  params_waste_heat_h   !< [-] anthropogenic heat outputs for heating e.g. 1 - 0.9 = 0.1 for combustion with eta = 0.9 or -2 for WP with COP = 3
+       REAL(wp) ::  phi_c_max             !< [W] Max. Cooling capacity (negative)
+       REAL(wp) ::  phi_h_max             !< [W] Max. Heating capacity (positive)
+       REAL(wp) ::  q_c_max               !< [W/m2] Max. Cooling heat flux per netto floor area (negative)
+       REAL(wp) ::  q_h_max               !< [W/m2] Max. Heating heat flux per netto floor area (positive)
+       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) ::  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_floor_area        !< [m2] netto ground area
+       REAL(wp) ::  s_layer3              !< [m] half thickness of the inner layer (layer_3)
+       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) ::  u_value_win           !< [W/(m2*K)] transmittance
+       REAL(wp) ::  vol_tot               !< [m3] total building volume
 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in       !< mean building indoor temperature, height dependent
@@ -171 +182 @@
 
     INTEGER(iwp) ::  num_build   !< total number of buildings in domain
-
 !
 !-- Declare all global variables within the module
@@ -179 +189 @@
     INTEGER(iwp) ::  solar_protection_on     !< Solar protection on
 
-    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) ::  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) ::  h_tr_1                      !<! [W/K] Heat transfer coefficient auxiliary variable 1
-    REAL(wp) ::  h_tr_2                      !<! [W/K] Heat transfer coefficient auxiliary variable 2
-    REAL(wp) ::  h_tr_3                      !<! [W/K] Heat transfer coefficient auxiliary variable 3 
-    REAL(wp) ::  h_tr_em                     !<! [W/K] Heat transfer coefficient of the emmision (got with h_tr_ms the thermal mass)
-    REAL(wp) ::  h_tr_is                     !<! [W/K] thermal coupling conductance (Thermischer Kopplungsleitwert)
-    REAL(wp) ::  h_tr_ms                     !<! [W/K] Heat transfer conductance term (got with h_tr_em the thermal mass)
-    REAL(wp) ::  h_tr_op                     !<! [W/K] heat transfer coefficient of opaque components (assumption: got all thermal mass) contains of h_tr_em and h_tr_ms
-    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) ::  initial_indoor_temperature  !< namelist parameter
-    REAL(wp) ::  net_sw_in                   !< net short-wave radiation (in - out; was i_global --> CORRECT?)
-    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 
-    REAL(wp) ::  phi_hc_nd_ac                !<! [W] actual heating demand and/or cooling demand
-    REAL(wp) ::  phi_hc_nd_un                !<! [W] unlimited heating demand and/or cooling demand which is necessary to reach the demanded required temperature (heating is positive, cooling is negative) 
-    REAL(wp) ::  phi_ia                      !< [W] internal air load = internal loads * 0.5, Eq. (C.1)
-    REAL(wp) ::  phi_m                       !<! [W] mass specific thermal load (internal and external)
-    REAL(wp) ::  phi_mtot                    !<! [W] total mass specific thermal load (internal and external)
-    REAL(wp) ::  phi_sol                     !< [W] solar loads
-    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_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) ::  schedule_d                  !< activation for internal loads (low or high + low)
+    REAL(wp) ::  a_m                           !< [m2] the effective mass-related area
+    REAL(wp) ::  air_change                    !< [1/h] Airflow
+    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) ::  facade_element_area           !< [m2_facade] building surface facade
+    REAL(wp) ::  floor_area_per_facade         !< [m2/m2] floor area per facade area
+    REAL(wp) ::  h_t_1                         !< [W/K] Heat transfer coefficient auxiliary variable 1
+    REAL(wp) ::  h_t_2                         !< [W/K] Heat transfer coefficient auxiliary variable 2
+    REAL(wp) ::  h_t_3                         !< [W/K] Heat transfer coefficient auxiliary variable 3 
+    REAL(wp) ::  h_t_wm                        !< [W/K] Heat transfer coefficient of the emmision (got with h_t_ms the thermal mass)
+    REAL(wp) ::  h_t_is                        !< [W/K] thermal coupling conductance (Thermischer Kopplungsleitwert)
+    REAL(wp) ::  h_t_ms                        !< [W/K] Heat transfer conductance term (got with h_t_wm the thermal mass)
+    REAL(wp) ::  h_t_wall                      !< [W/K] heat transfer coefficient of opaque components (assumption: got all thermal mass) contains of h_t_wm and h_t_ms
+    REAL(wp) ::  h_t_es                        !< [W/K] heat transfer coefficient of doors, windows, curtain walls and glazed walls (assumption: thermal mass=0)
+    REAL(wp) ::  h_v                           !< [W/K] heat transfer of ventilation
+    REAL(wp) ::  indoor_volume_per_facade      !< [m3] indoor air volume per facade element
+    REAL(wp) ::  initial_indoor_temperature    !< [K] initial indoor temperature (namelist parameter)
+    REAL(wp) ::  net_sw_in                     !< [W/m2] net short-wave radiation
+    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 
+    REAL(wp) ::  phi_hc_nd_ac                  !< [W] actual heating demand and/or cooling demand
+    REAL(wp) ::  phi_hc_nd_un                  !< [W] unlimited heating demand and/or cooling demand which is necessary to reach the demanded required temperature (heating is positive, cooling is negative) 
+    REAL(wp) ::  phi_ia                        !< [W] internal air load = internal loads * 0.5, Eq. (C.1)
+    REAL(wp) ::  phi_m                         !< [W] mass specific thermal load (internal and external)
+    REAL(wp) ::  phi_mtot                      !< [W] total mass specific thermal load (internal and external)
+    REAL(wp) ::  phi_sol                       !< [W] solar loads
+    REAL(wp) ::  phi_st                        !< [W] mass specific thermal load implied non thermal mass 
+    REAL(wp) ::  q_wall_win                    !< [W/m2]heat flux from indoor into wall/window
+    REAL(wp) ::  q_waste_heat                  !< [W/m2]waste heat, sum of waste heat over the roof to Palm
+    
+    REAL(wp) ::  q_c_m                         !< [W] Energy of thermal storage mass specific thermal load for internal and external heatsources (for energy bilanz)
+    REAL(wp) ::  q_c_st                        !< [W] Energy of thermal storage mass specific thermal load implied non thermal mass (for energy bilanz)
+    REAL(wp) ::  q_int                         !< [W] Energy of internal air load (for energy bilanz)
+    REAL(wp) ::  q_sol                         !< [W] Energy of solar (for energy bilanz)
+    REAL(wp) ::  q_trans                       !< [W] Energy of transmission (for energy bilanz)
+    REAL(wp) ::  q_vent                        !< [W] Energy of ventilation (for energy bilanz)
+    
+    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/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_m                     !<! [degree_C} inner temperature of the RC-node
-    REAL(wp) ::  theta_m_t                   !<! [degree_C] (Fictive) component temperature timestep
-    REAL(wp) ::  theta_m_t_prev              !< [degree_C] (Fictive) component temperature previous timestep (do not change)
-    REAL(wp) ::  theta_op                    !< [degree_C] operative temperature
-    REAL(wp) ::  theta_s                     !<! [degree_C] surface temperature of the RC-node
-    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) ::  f_sr                        !< [-] factor surface reduction
-    REAL(wp) ::  f_cei                       !< [-] ceiling reduction factor
-    REAL(wp) ::  ngs                         !< [m2] netto ground surface
-    
+    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/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_m                       !< [degree_C} inner temperature of the RC-node
+    REAL(wp) ::  theta_m_t                     !< [degree_C] (Fictive) component temperature timestep
+    REAL(wp) ::  theta_m_t_prev                !< [degree_C] (Fictive) component temperature previous timestep (do not change)
+    REAL(wp) ::  theta_op                      !< [degree_C] operative temperature
+    REAL(wp) ::  theta_s                       !< [degree_C] surface temperature of the RC-node
+    REAL(wp) ::  time_indoor = 0.0_wp          !< [s] time since last call of indoor model
+    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), PARAMETER ::  h_is                     = 3.45_wp     !< [W/(m2 K)] surface-related heat transfer coefficient between surface and air (chap. 7.2.2.2)
+    REAL(wp), PARAMETER ::  h_ms                     = 9.1_wp      !< [W/(m2 K)] surface-related heat transfer coefficient between component and surface (chap. 12.2.2)
     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) 
     REAL(wp), PARAMETER ::  params_f_win             = 0.5_wp      !< [-] proportion of window area, Database A_win aus Datenbank 27 window_area_per_facade_percent
     REAL(wp), PARAMETER ::  params_solar_protection  = 300.0_wp    !< [W/m2] chap. G.5.3.1 sun protection closed, if the radiation on facade exceeds this value
-    REAL(wp), PARAMETER ::  params_waste_heat_c      = 4.0_wp      !< [-] anthropogenic heat outputs for cooling e.g. 4 for KKM with COP = 3
-    REAL(wp), PARAMETER ::  params_waste_heat_h      = 1.111_wp    !< [-] anthropogenic heat outputs for heating e.g. 1 / 0.9 = 1.111111 for combustion with eta = 0.9 or -3 for WP with COP = 4
-    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
@@ -332 +340 @@
     REAL(wp) ::  near_facade_temperature
     REAL(wp) ::  phi_hc_nd_dummy
-   
-
-    !< Calculation of total mass specific thermal load (internal and external)
-    phi_mtot = ( phi_m + h_tr_em * indoor_wall_window_temperature              &
-                       + h_tr_3  * ( phi_st + h_tr_w * pt(k,j,i)               &
-                                            + h_tr_1 *                         &
-                                               ( ( ( phi_ia + phi_hc_nd_dummy ) / h_ve )  &
+!
+!-- Calculation of total mass specific thermal load (internal and external)
+    phi_mtot = ( phi_m + h_t_wm * indoor_wall_window_temperature               &
+                       + h_t_3  * ( phi_st + h_t_es * pt(k,j,i)                &
+                                            + h_t_1 *                          &
+                                    ( ( ( phi_ia + phi_hc_nd_dummy ) / h_v )   &
                                                  + near_facade_temperature )   &
-                                   ) / h_tr_2                                  &
+                                   ) / h_t_2                                   &
                )                                                                !< [degree_C] Eq. (C.5)
-    
-    !< Calculation of component temperature at factual timestep
-    theta_m_t = ( ( theta_m_t_prev                                               &
-                    * ( ( c_m / 3600.0_wp ) - 0.5_wp * ( h_tr_3 + h_tr_em ) ) + phi_mtot &
-                  )                                                              &
-                  /   ( ( c_m / 3600.0_wp ) + 0.5_wp * ( h_tr_3 + h_tr_em ) )            &
+! 
+!-- Calculation of component temperature at factual timestep
+    theta_m_t = ( ( theta_m_t_prev                                             &
+                    * ( ( c_m / 3600.0_wp ) - 0.5_wp * ( h_t_3 + h_t_wm ) )    &
+                     + phi_mtot                                                &
+                  )                                                            &
+                  /   ( ( c_m / 3600.0_wp ) + 0.5_wp * ( h_t_3 + h_t_wm ) )    &
                 )                                                               !< [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_wp                              !< [degree_C] Eq. (C.9)
-    
-    !< Calculation of mean surface temperature of the RC-node in actual timestep
-    theta_s = ( (   h_tr_ms * theta_m + phi_st + h_tr_w * pt(k,j,i)                         &
-                  + h_tr_1  * ( near_facade_temperature + ( phi_ia + phi_hc_nd_dummy ) / h_ve ) &
-                )                                                                           &
-                / ( h_tr_ms + h_tr_w + h_tr_1 )                                             &
+!
+!-- Calculation of mean inner temperature for the RC-node in actual timestep
+    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
+    theta_s = ( (   h_t_ms * theta_m + phi_st + h_t_es * pt(k,j,i)             &
+                  + h_t_1  * ( near_facade_temperature                         &
+                           + ( phi_ia + phi_hc_nd_dummy ) / h_v )              &
+                )                                                              &
+                / ( h_t_ms + h_t_es + h_t_1 )                                  &
               )                                                                 !< [degree_C] Eq. (C.10)
     
-    !< Calculation of the air temperature of the RC-node
-    theta_air = ( h_tr_is * theta_s + h_ve * near_facade_temperature         &
-                                    + phi_ia + phi_hc_nd_dummy ) / ( h_tr_is + h_ve ) !< [degree_C] Eq. (C.11)
+!
+!-- Calculation of the air temperature of the RC-node
+    theta_air = ( h_t_is * theta_s + h_v * near_facade_temperature             &
+                + phi_ia + phi_hc_nd_dummy ) / ( h_t_is + h_v )                 !< [degree_C] Eq. (C.11)
 
  END SUBROUTINE im_calc_temperatures
@@ -375 +386 @@
 !-- Input values
 !-- Input datas from Palm, M4
-!     i_global             -->  net_sw_in                         !global radiation [W/m2]
-!     theta_e              -->  pt(k,j,i)                         !undisturbed outside temperature, 1. PALM volume, for windows
+!     i_global             -->  net_sw_in                         !< global radiation [W/m2]
+!     theta_e              -->  pt(k,j,i)                         !< undisturbed outside temperature, 1. PALM volume, for windows
 !     theta_sup = theta_f  -->  surf_usm_h%pt_10cm(m)
-!                               surf_usm_v(l)%pt_10cm(m)   !Air temperature, facade near (10cm) air temperature from 1. Palm volume
+!                               surf_usm_v(l)%pt_10cm(m)          !< Air temperature, facade near (10cm) air temperature from 1. Palm volume
 !     theta_node           -->  t_wall_h(nzt_wall,m)
-!                               t_wall_v(l)%t(nzt_wall,m)         !Temperature of innermost wall layer, for opaque wall
+!                               t_wall_v(l)%t(nzt_wall,m)         !< Temperature of innermost wall layer, for opaque wall
 !------------------------------------------------------------------------------!
  SUBROUTINE im_init
@@ -431 +442 @@
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  receive_dum_h       !< dummy array used for MPI_ALLREDUCE  
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  receive_dum_v       !< dummy array used for MPI_ALLREDUCE  
-
+    
     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) ::  u_tmp                                     !< dummy for temporary calculation of u-value without h_is
+    REAL(wp) ::  du_tmp                                    !< 1/u_tmp
+    REAL(wp) ::  du_win_tmp                                !< 1/building(nb)%u_value_win
+
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  volume         !< total building volume at each discrete height level
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  volume_l       !< total building volume at each discrete height level,
                                                            !< on local subdomain
 
-    IF ( debug_output )  CALL debug_message( 'im_init', 'start' )
-
+    CALL location_message( 'initializing indoor model', .FALSE. )
 !
 !-- Initializing of indoor model is only possible if buildings can be 
@@ -457 +471 @@
           IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
              IF ( num_buildings_l(myid) > 0 )  THEN
-                IF ( ANY( building_id_f%var(j,i) ==  build_ids_l ) )  THEN
+                IF ( ANY( building_id_f%var(j,i) .EQ.  build_ids_l ) )  THEN 
                    CYCLE
                 ELSE
@@ -516 +530 @@
     build_ids = build_ids_l
 #endif
+
 !
 !-- Note: in parallel mode, building IDs can occur mutliple times, as 
@@ -524 +539 @@
 
        IF ( ALLOCATED(build_ids_final) )  THEN
-          IF ( ANY( build_ids(n) == build_ids_final ) )  THEN    !FK: Warum ANY?, Warum .EQ.? --> s.o
+          IF ( ANY( build_ids(n) == build_ids_final ) )  THEN
              CYCLE
           ELSE
@@ -623 +638 @@
 !--    these building is present. 
        IF ( buildings(nb)%on_pe )  THEN
+
           ALLOCATE( buildings(nb)%volume(buildings(nb)%kb_min:buildings(nb)%kb_max)   )
           ALLOCATE( buildings(nb)%vol_frac(buildings(nb)%kb_min:buildings(nb)%kb_max) )
@@ -659 +675 @@
        IF ( ALLOCATED( buildings(nb)%volume ) )                                &
           buildings(nb)%vol_tot = SUM( buildings(nb)%volume )
-       
+
        DEALLOCATE( volume   )
        DEALLOCATE( volume_l )
 
     ENDDO
-
 !
 !-- Allocate arrays for indoor temperature.  
@@ -716 +731 @@
     ENDDO
 !
-!-- Vertical facades
+!-- Vertical facades!
     buildings(:)%num_facades_per_building_v_l = 0
     DO  l = 0, 3
@@ -777 +792 @@
                            comm2d,                                             &
                            ierr )
-       IF ( ALLOCATED( buildings(nb)%num_facade_h ) )                          &  !FK: Was wenn not allocated? --> s.o.
+       IF ( ALLOCATED( buildings(nb)%num_facade_h ) )                          &
            buildings(nb)%num_facade_h = receive_dum_h
        IF ( ALLOCATED( buildings(nb)%num_facade_v ) )                          &
@@ -785 +800 @@
        buildings(nb)%num_facade_v = num_facades_v
 #endif
+
 !
 !--    Deallocate dummy arrays
@@ -812 +828 @@
        ENDIF
 !
-!--    Determine volume per facade element (vpf)
+! --    Determine volume per facade element (vpf)
        IF ( buildings(nb)%on_pe )  THEN
           ALLOCATE( buildings(nb)%vpf(buildings(nb)%kb_min:buildings(nb)%kb_max) )
@@ -818 +834 @@
           DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
              buildings(nb)%vpf(k) = buildings(nb)%volume(k) /                  &
-                                    ( buildings(nb)%num_facade_h(k) +          &
-                                      buildings(nb)%num_facade_v(k) )
+                                REAL( buildings(nb)%num_facade_h(k) +          &
+                                      buildings(nb)%num_facade_v(k), KIND = wp )
           ENDDO
+       ENDIF
+    
+!
+!--    Determine volume per total facade area (vpf). For the horizontal facade
+!--    area num_facades_per_building_h can be taken, multiplied with dx*dy. 
+!--    However, due to grid stretching, vertical facade elements must be
+!--    summed-up vertically. Please note, if dx /= dy, an error is made!
+       IF ( buildings(nb)%on_pe )  THEN
+         
+          buildings(nb)%vpf = buildings(nb)%vol_tot /                          &
+               ( buildings(nb)%num_facades_per_building_h * dx * dy +          &
+                 SUM( buildings(nb)%num_facade_v(buildings(nb)%kb_min:         &
+                                                 buildings(nb)%kb_max)         &
+                    * dzw(buildings(nb)%kb_min:buildings(nb)%kb_max) ) * dx )
        ENDIF
     ENDDO
@@ -863 +893 @@
 !-- 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)
+    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(:)%q_h_max             = building_pars(128,building_type)    
+    buildings(:)%q_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)
+    buildings(:)%params_waste_heat_h = building_pars(134,building_type)
+    buildings(:)%params_waste_heat_c = building_pars(135,building_type)
 !
 !-- Initialize ventilaation load. Please note, building types > 7 are actually
@@ -909 +941 @@
           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) ),    &
+                 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)                 
+                 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)%q_h_max             = building_pars(128,bt)    
+                 buildings(nb)%q_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) 
+                 buildings(nb)%params_waste_heat_h = building_pars(134,bt)
+                 buildings(nb)%params_waste_heat_c = building_pars(135,bt)
 !
 !--              Initialize ventilaation load. Please note, building types > 7 
@@ -958 +992 @@
         ENDDO
     ENDIF
+
+!
+!-- Calculation of surface-related heat transfer coeffiecient
+!-- out of standard u-values from building database 
+!-- only amount of extern and surface is used
+!-- otherwise amount between air and surface taken account twice
+    DO nb = 1, num_build
+       IF ( buildings(nb)%on_pe ) THEN   
+          du_win_tmp = 1.0_wp / buildings(nb)%u_value_win
+          u_tmp = buildings(nb)%u_value_win * ( du_win_tmp / ( du_win_tmp -    &
+                  0.125_wp + ( 1.0_wp / h_is ) ) )
+                  
+          du_tmp = 1.0_wp / u_tmp
+          
+          buildings(nb)%h_es = ( du_tmp / ( du_tmp - ( 1.0_wp / h_is ) ) ) *   &
+                                 u_tmp   
+       ENDIF
+    ENDDO
 !
 !-- Initial room temperature [K]
@@ -968 +1020 @@
     ENDDO
 
-    IF ( debug_output )  CALL debug_message( 'im_init', 'end' )
+    CALL location_message( 'finished', .TRUE. )
 
  END SUBROUTINE im_init
@@ -983 +1035 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzw, dzw
+
+!     USE basic_constants_and_equations_mod,                                     &
+!         ONLY:  c_p
+
+!     USE control_parameters,                                                    &
+!         ONLY:  rho_surface
 
     USE date_and_time_mod,                                                     &
@@ -1071 +1129 @@
 !            
 !--          Building geometries --> not time-dependent
-             facade_element_area      = dx * dy                             !< [m2] surface area per facade element    
-             floor_area_per_facade    = buildings(nb)%vpf(kk) * ddzw(kk)    !< [m2] net floor area per facade element       
-             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
-             
-! 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                    = buildings(nb)%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
-
+             facade_element_area          = dx * dy                                                   !< [m2] surface area per facade element    
+             floor_area_per_facade        = buildings(nb)%vpf(kk) * ddzw(kk)                          !< [m2/m2] floor area per facade area
+             indoor_volume_per_facade     = buildings(nb)%vpf(kk)                                     !< [m3/m2] indoor air volume per facade area 
+             buildings(nb)%area_facade    = facade_element_area *                                   &
+                                          ( buildings(nb)%num_facades_per_building_h +              &
+                                            buildings(nb)%num_facades_per_building_v )                !< [m2] area of total facade
+             window_area_per_facade       = surf_usm_h%frac(ind_wat_win,m)  * facade_element_area     !< [m2] window area per facade element
+
+             buildings(nb)%net_floor_area = buildings(nb)%vol_tot / ( buildings(nb)%height_storey )
+             total_area                   = buildings(nb)%net_floor_area                              !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
+             a_m                          = buildings(nb)%factor_a * total_area *                   &
+                                          ( facade_element_area / buildings(nb)%area_facade ) *     &
+                                            buildings(nb)%lambda_at                                   !< [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 * total_area *                   &
+                                          ( facade_element_area / buildings(nb)%area_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)             
+!
 !--          Calculation of heat transfer coefficient for transmission --> not time-dependent
-             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.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
+             h_t_es   = window_area_per_facade * buildings(nb)%h_es                                   !< [W/K] only for windows
+
+             h_t_is  = buildings(nb)%area_facade  * h_is                                                             !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
+             h_t_ms  = a_m * h_ms                                                                     !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
+             h_t_wall  = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade )     & !< [W/K] 
+                                    * buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp &
+                                             ) + 1.0_wp / h_t_ms )                                    !< [W/K] opaque components
+             h_t_wm  = 1.0_wp / ( 1.0_wp / h_t_wall - 1.0_wp / h_t_ms )                               !< [W/K] emmision 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_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )            &
-                              * ngs )         !< [W] Eq. (C.1)
+             phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )  &
+                              * floor_area_per_facade ) 
+             q_int = phi_ia / total_area
 !
 !--          Airflow dependent on the occupacy of the room
@@ -1110 +1168 @@
 !--          not less than 0.01 W/K to provide division by 0 in further calculations
 !--          with heat capacity of air 0.33 Wh/m2K
-             h_ve   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *      &
+             h_v   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *      &
                                     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.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.0_wp / ( ( 1.0_wp / h_tr_2 ) + ( 1.0_wp / h_tr_ms ) )  !< [W/K] Eq. (C.8)
+             h_t_1 = 1.0_wp / ( ( 1.0_wp / h_v )   + ( 1.0_wp / h_t_is ) )  !< [W/K] Eq. (C.6)
+             h_t_2 = h_t_1 + h_t_es                                         !< [W/K] Eq. (C.7)
+             h_t_3 = 1.0_wp / ( ( 1.0_wp / h_t_2 ) + ( 1.0_wp / h_t_ms ) )  !< [W/K] Eq. (C.8)
 !
 !--          Net short-wave radiation through window area (was i_global)
@@ -1129 +1187 @@
                   surf_usm_h%frac(ind_veg_wall,m) * t_wall_h(nzt_wall,m)       &
                 + surf_usm_h%frac(ind_wat_win,m)  * t_window_h(nzt_wall,m)
+!
+!--          Solar thermal gains. If net_sw_in larger than sun-protection 
+!--          threshold parameter (params_solar_protection), sun protection will 
+!--          be activated
+             IF ( net_sw_in <= params_solar_protection )  THEN 
+                solar_protection_off = 1
+                solar_protection_on  = 0
+             ELSE 
+                solar_protection_off = 0
+                solar_protection_on  = 1
+             ENDIF
+!
+!--          Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection 
+!--          DIN 4108 - 2 chap.8
+             phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off                           &
+                         + 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
+             q_sol = phi_sol           
+!
+!--          Calculation of the mass specific thermal load for internal and external heatsources of the inner node
+             phi_m   = (a_m / total_area) * ( phi_ia + phi_sol )                                    !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
+             q_c_m = phi_m
+!
+!--          Calculation mass specific thermal load implied non thermal mass 
+             phi_st  =   ( 1.0_wp - ( a_m / total_area ) - ( h_t_es / ( 9.1_wp * total_area ) ) ) &
+                       * ( phi_ia + phi_sol )                                                       !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
+             q_c_st = phi_st            
+!
+!--          Calculations for deriving indoor temperature and heat flux into the wall
+!--          Step 1: Indoor temperature without heating and cooling
+!--          section C.4.1 Picture C.2 zone 3)
+             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 ( 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/m2 heating/cooling power and calculate indoor temperature
+!--          again.
+             ELSE
+!
+!--             Temperature not correct, calculation method according to section C4.2
+                theta_air_0 = theta_air                                                  !< temperature without heating/cooling  
+!
+!--             Heating or cooling?
+                IF ( theta_air_0 > buildings(nb)%theta_int_c_set )  THEN
+                   theta_air_set = buildings(nb)%theta_int_c_set
+                ELSE 
+                   theta_air_set = buildings(nb)%theta_int_h_set 
+                ENDIF
+!
+!--             Calculate the temperature with phi_hc_nd_10 
+                phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
+                phi_hc_nd    = phi_hc_nd_10
+                
+                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
+                                            near_facade_temperature, phi_hc_nd )
+                theta_air_10 = theta_air                                                !< temperature with 10 W/m2 of heating
+                phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0)          &
+                                            / (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)
+                buildings(nb)%phi_c_max = buildings(nb)%q_c_max * floor_area_per_facade              
+                buildings(nb)%phi_h_max = buildings(nb)%q_h_max * floor_area_per_facade
+                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  
+                ELSE
+!
+!--             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.0_wp )  THEN
+                      phi_hc_nd_ac = buildings(nb)%phi_h_max                            !< Limit heating
+                   ELSE 
+                      phi_hc_nd_ac = buildings(nb)%phi_c_max                            !< Limit cooling
+                   ENDIF
+                ENDIF
+                phi_hc_nd = phi_hc_nd_ac   
+!
+!--             Calculate the temperature with phi_hc_nd_ac (new)
+                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
+                                            near_facade_temperature, phi_hc_nd )
+                theta_air_ac = theta_air
+             ENDIF
+!
+!--          Update theta_m_t_prev
+             theta_m_t_prev = theta_m_t
+             
+             q_vent = h_v * ( theta_air - near_facade_temperature )
+!
+!--          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_wp * theta_air_ac + 0.7_wp * theta_s          !< [degree_C] operative Temperature Eq. (C.12)
+!              surf_usm_h%t_indoor(m) = theta_op                               !< not integrated now
+!
+!--          Heat flux into the wall. Value needed in urban_surface_mod to 
+!--          calculate heat transfer through wall layers towards the facade
+!--          (use c_p * rho_surface to convert [W/m2] into [K m/s])
+             q_wall_win = h_t_ms * ( theta_s - theta_m )                       &
+                                    / (   facade_element_area                  &
+                                        - window_area_per_facade )
+             q_trans = q_wall_win * facade_element_area                                        
+!
+!--          Transfer q_wall_win back to USM (innermost wall/window layer)
+             surf_usm_h%iwghf_eb(m)        = q_wall_win
+             surf_usm_h%iwghf_eb_window(m) = q_wall_win
+!
+!--          Sum up operational indoor temperature per kk-level. Further below,
+!--          this temperature is reduced by MPI to one temperature per kk-level
+!--          and building (processor overlapping)
+             buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
+!
+!--          Calculation of waste heat
+!--          Anthropogenic heat output
+             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 * (                                    &
+                              buildings(nb)%params_waste_heat_h * heating_on + &
+                              buildings(nb)%params_waste_heat_c * cooling_on ) &
+                            ) / facade_element_area                                               !< [W/m2] , observe the directional convention in PALM!
+             surf_usm_h%waste_heat(m) = q_waste_heat
+          ENDDO !< Horizontal surfaces loop
+!
+!--       Vertical surfaces
+          DO  fa = 1, buildings(nb)%num_facades_per_building_v_l
+!
+!--          Determine indices where corresponding surface-type information 
+!--          is stored. 
+             l = buildings(nb)%l_v(fa)
+             m = buildings(nb)%m_v(fa)
+!
+!--          Determine building height level index. 
+             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
+             IF ( l == 2  .OR. l == 3 ) facade_element_area = dy * dzw(kk)                            !< [m2] surface area per facade element
+             floor_area_per_facade        = buildings(nb)%vpf(kk) * ddzw(kk)                          !< [m2/m2] floor area per facade area
+             indoor_volume_per_facade     = buildings(nb)%vpf(kk)                                     !< [m3/m2] indoor air volume per facade area 
+             buildings(nb)%area_facade    = facade_element_area *                                   &
+                                          ( buildings(nb)%num_facades_per_building_h +              &
+                                            buildings(nb)%num_facades_per_building_v )                !< [m2] area of total facade
+             window_area_per_facade       = surf_usm_v(l)%frac(ind_wat_win,m)  * facade_element_area     !< [m2] window area per facade element
+
+             buildings(nb)%net_floor_area = buildings(nb)%vol_tot / ( buildings(nb)%height_storey )
+             total_area                   = buildings(nb)%net_floor_area                              !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
+             a_m                          = buildings(nb)%factor_a * total_area *                   &
+                                          ( facade_element_area / buildings(nb)%area_facade ) *     &
+                                            buildings(nb)%lambda_at                                   !< [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 * total_area *                   &
+                                          ( facade_element_area / buildings(nb)%area_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)
+!
+!--          Calculation of heat transfer coefficient for transmission --> not time-dependent
+             h_t_es   = window_area_per_facade * buildings(nb)%h_es                                   !< [W/K] only for windows
+
+             h_t_is  = buildings(nb)%area_facade  * h_is                                                             !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
+             h_t_ms  = a_m * h_ms                                                                     !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
+             h_t_wall  = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade )     & !< [W/K] 
+                                    * buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp &
+                                             ) + 1.0_wp / h_t_ms )                                    !< [W/K] opaque components
+             h_t_wm  = 1.0_wp / ( 1.0_wp / h_t_wall - 1.0_wp / h_t_ms )                               !< [W/K] emmision 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_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )  &
+                              * floor_area_per_facade )
+             q_int = phi_ia
+
+!
+!--          Airflow dependent on the occupacy of the room
+!--          basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
+             air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low )  
+!
+!--          Heat transfer of ventilation 
+!--          not less than 0.01 W/K to provide division by 0 in further calculations
+!--          with heat capacity of air 0.33 Wh/m2K
+             h_v   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *                       &
+                                    0.33_wp * (1.0_wp - buildings(nb)%eta_ve ) ) )                    !< [W/K] from ISO 13789 Eq.(10)
+                                    
+!--          Heat transfer coefficient auxiliary variables
+             h_t_1 = 1.0_wp / ( ( 1.0_wp / h_v )   + ( 1.0_wp / h_t_is ) )                            !< [W/K] Eq. (C.6)
+             h_t_2 = h_t_1 + h_t_es                                                                   !< [W/K] Eq. (C.7)
+             h_t_3 = 1.0_wp / ( ( 1.0_wp / h_t_2 ) + ( 1.0_wp / h_t_ms ) )                            !< [W/K] Eq. (C.8)
+!
+!--          Net short-wave radiation through window area (was i_global)
+             net_sw_in = surf_usm_v(l)%rad_sw_in(m) - surf_usm_v(l)%rad_sw_out(m)
+!
+!--          Quantities needed for im_calc_temperatures
+             i = surf_usm_v(l)%i(m)
+             j = surf_usm_v(l)%j(m)   
+             k = surf_usm_v(l)%k(m)
+             near_facade_temperature = surf_usm_v(l)%pt_10cm(m)
+             indoor_wall_window_temperature =                                                       &
+                  surf_usm_v(l)%frac(ind_veg_wall,m) * t_wall_v(l)%t(nzt_wall,m)                    &
+                + surf_usm_v(l)%frac(ind_wat_win,m)  * t_window_v(l)%t(nzt_wall,m)
 !
 !--          Solar thermal gains. If net_sw_in larger than sun-protection 
@@ -1141 +1409 @@
              ENDIF
 !
-!--          Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection 
+!--          Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
 !--          DIN 4108 - 2 chap.8
-             phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off               &
+             phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off                             &
                          + 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
-             phi_m   = (eff_mass_area / total_area) * ( phi_ia + phi_sol )          !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
+                       * buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w
+             q_sol = phi_sol
+!
+!--          Calculation of the mass specific thermal load for internal and external heatsources
+             phi_m   = (a_m / total_area) * ( phi_ia + phi_sol )          !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
+             q_c_m = phi_m
 !
 !--          Calculation mass specific thermal load implied non thermal mass 
-             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
+             phi_st  =   ( 1.0_wp - ( a_m / total_area ) - ( h_t_es / ( 9.1_wp * total_area ) ) )                &
+                       * ( phi_ia + phi_sol )                                                                       !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
+             q_c_st = phi_st 
 !
 !--          Calculations for deriving indoor temperature and heat flux into the wall
@@ -1158 +1429 @@
 !--          section C.4.1 Picture C.2 zone 3)
              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 air temperature between border temperatures of heating and cooling, assign output variable, then ready  
              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
@@ -1174 +1444 @@
 !--             Temperature not correct, calculation method according to section C4.2
                 theta_air_0 = theta_air !< Note temperature without heating/cooling
-
+!
 !--             Heating or cooling?
-                IF ( theta_air > buildings(nb)%theta_int_c_set )  THEN
+                IF ( theta_air_0 > buildings(nb)%theta_int_c_set )  THEN
                    theta_air_set = buildings(nb)%theta_int_c_set
                 ELSE 
@@ -1182 +1452 @@
                 ENDIF
 
-!--             Calculate the temperature with phi_hc_nd_10 
+!--             Calculate the temperature with phi_hc_nd_10
                 phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
                 phi_hc_nd    = phi_hc_nd_10
@@ -1190 +1460 @@
 
                 theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
-!
-   
-                phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0)                        &
-                                            / (theta_air_10  - theta_air_0)            !< Eq. (C.13)
-                                            
+
+                
+                phi_hc_nd_un = phi_hc_nd_10 * ( theta_air_set - theta_air_0 )  &
+                                            / ( 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)
+                buildings(nb)%phi_c_max = buildings(nb)%q_c_max * floor_area_per_facade
+                buildings(nb)%phi_h_max = buildings(nb)%q_h_max * floor_area_per_facade
                 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    
+                   phi_hc_nd = phi_hc_nd_un
                 ELSE
+!
 !--             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)
@@ -1210 +1483 @@
                    ENDIF
                 ENDIF
-       
-                phi_hc_nd = phi_hc_nd_ac   
+                phi_hc_nd = phi_hc_nd_ac  
 !
 !--             Calculate the temperature with phi_hc_nd_ac (new)
                 CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
                                             near_facade_temperature, phi_hc_nd )
-
                 theta_air_ac = theta_air
-       
              ENDIF
 !
 !--          Update theta_m_t_prev
              theta_m_t_prev = theta_m_t
-!
-!--          Calculate the operating temperature with weighted mean temperature of air and mean solar temperature
+             
+             q_vent = h_v * ( theta_air - near_facade_temperature )
+!
+!--          Calculate the operating temperature with weighted mean of temperature of air and mean
 !--          Will be used for thermal comfort calculations 
-             theta_op     = 0.3_wp * theta_air_ac + 0.7_wp * theta_s          !< [degree_C] operative Temperature Eq. (C.12)
+             theta_op     = 0.3_wp * theta_air_ac + 0.7_wp * theta_s
+!              surf_usm_v(l)%t_indoor(m) = theta_op                  !< not integrated yet
 !
 !--          Heat flux into the wall. Value needed in urban_surface_mod to 
 !--          calculate heat transfer through wall layers towards the facade
-!--          (use c_p * rho_surface to convert [W/m2] into [K m/s])
-             q_wall_win = h_tr_ms * ( theta_s - theta_m )                    &
+             q_wall_win = h_t_ms * ( theta_s - theta_m )                       &
                                     / (   facade_element_area                  &
                                         - window_area_per_facade )
+             q_trans = q_wall_win * facade_element_area
 !
 !--          Transfer q_wall_win back to USM (innermost wall/window layer)
-             surf_usm_h%iwghf_eb(m)        = q_wall_win
-             surf_usm_h%iwghf_eb_window(m) = q_wall_win
+             surf_usm_v(l)%iwghf_eb(m)        = q_wall_win
+             surf_usm_v(l)%iwghf_eb_window(m) = q_wall_win
 !
 !--          Sum up operational indoor temperature per kk-level. Further below,
@@ -1251 +1524 @@
              ELSE 
                 heating_on = 0
-                cooling_on = 1
+                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
+             q_waste_heat = ( phi_hc_nd * (                                    &
+                    buildings(nb)%params_waste_heat_h * heating_on +           &
+                    buildings(nb)%params_waste_heat_c * cooling_on )           &
+                            ) / facade_element_area !< [W/m2] , observe the directional convention in PALM!
+             surf_usm_v(l)%waste_heat(m) = q_waste_heat
              
-          ENDDO !< Horizontal surfaces loop
-!
-!--       Vertical surfaces
-          DO  fa = 1, buildings(nb)%num_facades_per_building_v_l
-!
-!--          Determine indices where corresponding surface-type information 
-!--          is stored. 
-             l = buildings(nb)%l_v(fa)
-             m = buildings(nb)%m_v(fa)
-!
-!--          Determine building height level index. 
-             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
-             IF ( l == 2  .OR. l == 3 ) facade_element_area = dy * dzw(kk)  !< [m2] surface area per facade element
-             floor_area_per_facade    = buildings(nb)%vpf(kk) * ddzw(kk)    !< [m2] net floor area per facade element       
-             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
-
-             f_cei                    = buildings(nb)%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 * 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.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_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 = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low )  
-!
-!--          Heat transfer of ventilation 
-!--          not less than 0.01 W/K to provide division by 0 in further calculations
-!--          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 - buildings(nb)%eta_ve ) ) )                    !< [W/K] from ISO 13789 Eq.(10)
-                                    
-!--          Heat transfer coefficient auxiliary variables
-             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.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)
-             net_sw_in = surf_usm_v(l)%rad_sw_in(m) - surf_usm_v(l)%rad_sw_out(m)
-!
-!--          Quantities needed for im_calc_temperatures
-             i = surf_usm_v(l)%i(m)
-             j = surf_usm_v(l)%j(m)
-             k = surf_usm_v(l)%k(m)
-             near_facade_temperature = surf_usm_v(l)%pt_10cm(m)
-             indoor_wall_window_temperature =                                    &
-                  surf_usm_v(l)%frac(ind_veg_wall,m) * t_wall_v(l)%t(nzt_wall,m) &
-                + surf_usm_v(l)%frac(ind_wat_win,m)  * t_window_v(l)%t(nzt_wall,m)
-!
-!--          Solar thermal gains. If net_sw_in larger than sun-protection 
-!--          threshold parameter (params_solar_protection), sun protection will 
-!--          be activated
-             IF ( net_sw_in <= params_solar_protection )  THEN 
-                solar_protection_off = 1
-                solar_protection_on  = 0 
-             ELSE 
-                solar_protection_off = 0
-                solar_protection_on  = 1 
-             ENDIF
-!
-!--          Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
-!--          DIN 4108 - 2 chap.8
-             phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off               &
-                         + 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
-             phi_m   = (eff_mass_area / total_area) * ( phi_ia + phi_sol )          !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
-!
-!--          Calculation mass specific thermal load implied non thermal mass 
-             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
-!
-!--          Calculations for deriving indoor temperature and heat flux into the wall
-!--          Step 1: Indoor temperature without heating and cooling
-!--          section C.4.1 Picture C.2 zone 3)
-             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 ( 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/m2 heating/cooling power and calculate indoor temperature
-!--          again.
-             ELSE
-!
-!--             Temperature not correct, calculation method according to section C4.2
-                theta_air_0 = theta_air !< Note temperature without heating/cooling
-
-!--             Heating or cooling?
-                IF ( theta_air > buildings(nb)%theta_int_c_set )  THEN
-                   theta_air_set = buildings(nb)%theta_int_c_set
-                ELSE 
-                   theta_air_set = buildings(nb)%theta_int_h_set 
-                ENDIF
-
-!--             Calculate the temperature with phi_hc_nd_10
-                phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
-                phi_hc_nd    = phi_hc_nd_10
-       
-                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
-                                            near_facade_temperature, phi_hc_nd )
-
-                theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
-
-                
-                phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0)                        &
-                                            / (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 ( 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
-                ELSE
-!--             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.0_wp )  THEN
-                      phi_hc_nd_ac = buildings(nb)%phi_h_max                                         !< Limit heating
-                   ELSE 
-                      phi_hc_nd_ac = buildings(nb)%phi_c_max                                         !< Limit cooling
-                   ENDIF
-                ENDIF
-       
-                phi_hc_nd = phi_hc_nd_ac   
-!
-!--             Calculate the temperature with phi_hc_nd_ac (new)
-                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
-                                            near_facade_temperature, phi_hc_nd )
-
-                theta_air_ac = theta_air
-       
-             ENDIF
-!
-!--          Update theta_m_t_prev
-             theta_m_t_prev = theta_m_t
-!
-!--          Calculate the operating temperature with weighted mean of temperature of air and mean
-!--          Will be used for thermal comfort calculations 
-             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 
-!--          calculate heat transfer through wall layers towards the facade
-             q_wall_win = h_tr_ms * ( theta_s - theta_m )                    &
-                                    / (   facade_element_area                  &
-                                        - window_area_per_facade )
-!
-!--          Transfer q_wall_win back to USM (innermost wall/window layer)
-             surf_usm_v(l)%iwghf_eb(m)        = q_wall_win
-             surf_usm_v(l)%iwghf_eb_window(m) = q_wall_win
-!
-!--          Sum up operational indoor temperature per kk-level. Further below,
-!--          this temperature is reduced by MPI to one temperature per kk-level
-!--          and building (processor overlapping)
-             buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
-
-!
-!--          Calculation of waste heat
-!--          Anthropogenic heat output
-             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
 
@@ -1473 +1553 @@
           t_in_l_send = buildings(nb)%t_in_l
        ENDIF
+
 
 #if defined( __parallel ) 
@@ -1489 +1570 @@
           buildings(nb)%t_in = buildings(nb)%t_in_l
 #endif
-       IF ( ALLOCATED( buildings(nb)%t_in ) )                                  &
-          buildings(nb)%t_in =   buildings(nb)%t_in /                          &
-                               ( buildings(nb)%num_facade_h +                  &
-                                 buildings(nb)%num_facade_v )
+
+       IF ( ALLOCATED( buildings(nb)%t_in ) )  THEN
+!
+!--       Average indoor temperature. Note, in case a building is completely
+!--       surrounded by higher buildings, it may have no facade elements
+!--       at some height levels, will will lead to a divide by zero. If this
+!--       is the case, indoor temperature will be set to -1.0. 
+          DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
+             IF ( buildings(nb)%num_facade_h(k) +                              &
+                  buildings(nb)%num_facade_v(k) > 0 )  THEN
+                buildings(nb)%t_in(k) = buildings(nb)%t_in(k) /                &
+                               REAL( buildings(nb)%num_facade_h(k) +           &
+                                     buildings(nb)%num_facade_v(k), KIND = wp )
+             ELSE
+                buildings(nb)%t_in(k) = -1.0_wp
+             ENDIF
+          ENDDO
+       ENDIF
+       
+
 !
 !--    Deallocate dummy arrays
@@ -1499 +1596 @@
 
     ENDDO
-
-
+    
  END SUBROUTINE im_main_heatcool
 
@@ -1530 +1626 @@
            unit = 'W m-2'
         
-        CASE ( 'im_t_indoor' )
+        CASE ( 'im_t_indoor_mean' )
+           unit = 'K'
+           
+        CASE ( 'im_t_indoor_roof' )
+           unit = 'K'
+           
+        CASE ( 'im_t_indoor_wall_win' )
            unit = 'K'
         
@@ -1548 +1650 @@
  SUBROUTINE im_check_parameters
 
-!!!!   USE control_parameters,
-!!!!       ONLY: message_string
+!   USE control_parameters,
+!       ONLY: message_string
        
    IMPLICIT NONE
@@ -1587 +1689 @@
           grid_y = 'y'
           grid_z = 'zu'
-          
-       CASE ( 'im_t_indoor' )
+
+       CASE ( 'im_t_indoor_mean', 'im_t_indoor_roof', 'im_t_indoor_wall_win')
           grid_x = 'x'
           grid_y = 'y'
@@ -1627 +1729 @@
     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 !< 
@@ -1643 +1744 @@
 !--     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' )
+        CASE ( 'im_t_indoor_mean' )
            IF ( av == 0 ) THEN
               DO  i = nxl, nxr
@@ -1676 +1777 @@
               ENDDO
            ENDIF 
-    
+
         CASE ( 'im_hf_roof_waste' )
            IF ( av == 0 ) THEN
@@ -1685 +1786 @@
                  local_pf(i,j,k) = surf_usm_h%waste_heat(m)
               ENDDO
-           ENDIF                     
-       
+           ENDIF
+
        CASE ( 'im_hf_wall_win' )
            IF ( av == 0 ) THEN
@@ -1698 +1799 @@
               ENDDO
            ENDIF
-           
+
         CASE ( 'im_hf_wall_win_waste' )
            IF ( av == 0 ) THEN
@@ -1709 +1810 @@
                  ENDDO
               ENDDO
-           ENDIF      
-    
+           ENDIF
+
+!
+!< NOTE im_t_indoor_roof and im_t_indoor_wall_win not work yet
+
+!         CASE ( 'im_t_indoor_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%t_indoor(m)
+!               ENDDO
+!            ENDIF
+! 
+!         CASE ( 'im_t_indoor_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)%t_indoor(m)
+!                  ENDDO
+!               ENDDO
+!            ENDIF
+
         CASE DEFAULT
            found = .FALSE.

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Comments added
+! 
+! 4134 2019-08-02 18:39:57Z suehring
 ! Bugfix in formatted write statement
 ! 
@@ -3569 +3572 @@
 
        IF ( average_radiation ) THEN
-
+!
+!--       Determine minimum topography top index. 
           k_topo_l = MINVAL( get_topography_top_index( 's' ) )
 #if defined( __parallel )
@@ -3705 +3709 @@
                              rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1)
           
-!              CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
-!              rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
-!              rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
-!              rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
-!              rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
-!              rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
-!              rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
-!              rrtm_reliq      , rrtm_lw_tauaer,                               &
-!              rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
-!              rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
-!              rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
-
              CALL rrtmg_lw( 1,                                                 &                                        
                             nzt_rad-k_topo,                                    &
@@ -3809 +3801 @@
              rrtm_sw_asmaer_dum = rrtm_sw_asmaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
              rrtm_sw_ecaer_dum  = rrtm_sw_ecaer(0:0,k_topo+1:nzt_rad+1,1:naerec+1)
-!              CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld     , rrtm_iaer      ,&
-!              rrtm_play      , rrtm_plev     , rrtm_tlay     , rrtm_tlev      ,&
-!              rrtm_tsfc      , rrtm_h2ovmr   , rrtm_o3vmr    , rrtm_co2vmr    ,&
-!              rrtm_ch4vmr    , rrtm_n2ovmr   , rrtm_o2vmr    , rrtm_asdir     ,&
-!              rrtm_asdif     , rrtm_aldir    , rrtm_aldif    , zenith         ,&
-!              0.0_wp         , day_of_year   , solar_constant, rrtm_inflgsw   ,&
-!              rrtm_iceflgsw  , rrtm_liqflgsw , rrtm_cldfr    , rrtm_sw_taucld ,&
-!              rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp    ,&
-!              rrtm_cliqwp    , rrtm_reice    , rrtm_reliq    , rrtm_sw_tauaer ,&
-!              rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx    ,&
-!              rrtm_swdflx    , rrtm_swhr     , rrtm_swuflxc  , rrtm_swdflxc   ,&
-!              rrtm_swhrc     , rrtm_dirdflux , rrtm_difdflux )
+
              CALL rrtmg_sw( 1,                                                 &
                             nzt_rad-k_topo,                                    &
@@ -3879 +3860 @@
   
 !
-!--          Save fluxes:
-!--          - whole domain
+!--          Save radiation fluxes for the entire depth of the model domain
              DO k = nzb, nzt+1
                 rad_sw_in(k,:,:)  = rrtm_swdflx(0,k)
                 rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
              ENDDO
-!--          - direct and diffuse SW at urban-surface-layer (required by RTM)
+!--          Save direct and diffuse SW radiation at the surface (required by RTM)
              rad_sw_in_dir(:,:) = rrtm_dirdflux(0,k_topo)
              rad_sw_in_diff(:,:) = rrtm_difdflux(0,k_topo)

palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Remove unused variable
+! 
+! 4144 2019-08-06 09:11:47Z raasch
 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
 ! 
@@ -2269 +2272 @@
         
     USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
 
     IMPLICIT NONE

palm/trunk/SOURCE/urban_surface_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! - Add anthropogenic heat output factors for heating and cooling to building
+!   data base
+! - Move definition of building_pars to usm_init_arrays since it is already
+!   required in the indoor model
+! 
+! 4127 2019-07-30 14:47:10Z suehring
 ! Do not add anthopogenic energy during wall/soil spin-up 
 ! (merge from branch resler)
@@ -781 +787 @@
 !-- (>700), while intel and gfortran compiler have hard limit of continuation 
 !-- lines of 511.
-    REAL(wp), DIMENSION(0:133,1:7) ::  building_pars
+    REAL(wp), DIMENSION(0:135,1:7) ::  building_pars
 !
 !-- Type for surface temperatures at vertical walls. Is not necessary for horizontal walls. 
@@ -1430 +1436 @@
            IF ( ALLOCATED( surf_usm_v(l)%iwghf_eb_window ) )  surf_usm_v(l)%iwghf_eb_window = 0.0_wp
         ENDDO
-
+!
+!--     Initialize building-surface properties, which are also required by other modules, 
+!--     e.g. the indoor model.
+        CALL usm_define_pars
+        
         IF ( debug_output )  CALL debug_message( 'usm_init_arrays', 'end' )
         
@@ -3778 +3788 @@
 
         CALL cpu_log( log_point_s(78), 'usm_init', 'start' )
-!
-!--     Initialize building-surface properties
-        CALL usm_define_pars
 !
 !--     surface forcing have to be disabled for LSF 
@@ -9068 +9075 @@
            10.0_wp,        &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           0.1_wp,         &  !< parameter 134 - anthropogenic heat output for heating
+           1.333_wp        &  !< parameter 135 - anthropogenic heat output for cooling
                             /)
                             
@@ -9205 +9214 @@
            8.0_wp,         &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           0.1_wp,         &  !< parameter 134 - anthropogenic heat output for heating
+           1.333_wp        &  !< parameter 135 - anthropogenic heat output for cooling
                             /)
                             
@@ -9342 +9353 @@
            8.0_wp,         &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           0.1_wp,         &  !< parameter 134 - anthropogenic heat output for heating
+           1.333_wp        &  !< parameter 135 - anthropogenic heat output for cooling
                             /)   
                             
@@ -9479 +9492 @@
            10.0_wp,        &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           0.1_wp,         &  !< parameter 134 - anthropogenic heat output for heating
+           1.333_wp        &  !< parameter 135 - anthropogenic heat output for cooling
                             /)   
                             
@@ -9616 +9631 @@
            20.0_wp,        &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           0.0_wp,         &  !< parameter 134 - anthropogenic heat output for heating
+           2.54_wp         &  !< parameter 135 - anthropogenic heat output for cooling
                             /)
                             
@@ -9753 +9770 @@
            15.0_wp,        &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           -2.0_wp,        &  !< parameter 134 - anthropogenic heat output for heating
+           1.25_wp         &  !< parameter 135 - anthropogenic heat output for cooling
                             /)
                             
@@ -9890 +9909 @@
            0.0_wp,         &  !< parameter 131 - basic internal heat gains without occupancy of the room
            3.0_wp,         &  !< parameter 132 - storey height
-           0.2_wp          &  !< parameter 133 - ceiling construction height
+           0.2_wp,         &  !< parameter 133 - ceiling construction height
+           0.0_wp,         &  !< parameter 134 - anthropogenic heat output for heating
+           0.0_wp          &  !< parameter 135 - anthropogenic heat output for cooling
                         /)