Index: palm/trunk/SOURCE/indoor_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/indoor_model_mod.f90	(revision 4681)
+++ palm/trunk/SOURCE/indoor_model_mod.f90	(revision 4687)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Bugfix: values of theta_m_t_prev were not stored for individual surfaces and thus re-used by all
+! surfaces and buildings, which led to excessive indoor temperatures
+!
+! 4681 2020-09-16 10:23:06Z pavelkrc
 ! Bugfix for implementation of downward surfaces
 !
@@ -126,4 +130,5 @@
 ! @author Matthias Suehring
 ! @author Sascha RiÃmann
+! @author BjÃ¶rn Maronga
 !
 !
@@ -144,5 +149,5 @@
 !> @note How to write indoor temperature output to pt array?
 !>
-!> @bug  <Enter known bugs here>
+!> @bug  Calculation of iwghf_eb and iwghf_eb_window is faulty
 !--------------------------------------------------------------------------------------------------!
  MODULE indoor_model_mod
@@ -234,4 +239,6 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in       !< mean building indoor temperature, height dependent
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in_l     !< mean building indoor temperature on local subdomain, height dependent
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  theta_m_t_prev_h !< [degree_C] value of theta_m_t from previous time step (horizontal)
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  theta_m_t_prev_v !< [degree_C] value of theta_m_t from previous time step (vertical)
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  volume     !< total building volume, height dependent
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  vol_frac   !< fraction of local on total building volume, height dependent
@@ -319,6 +326,5 @@
     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_m_t                           !< [degree_C] (Fictive) component temperature during timestep
     REAL(wp) ::  theta_op                            !< [degree_C] operative temperature
     REAL(wp) ::  theta_s                             !< [degree_C] surface temperature of the RC-node
@@ -449,5 +455,5 @@
 !--------------------------------------------------------------------------------------------------!
  SUBROUTINE im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,                        &
-                                   near_facade_temperature, phi_hc_nd_dummy )
+                                   near_facade_temperature, phi_hc_nd_dummy, theta_m_t_prev )
 
     INTEGER(iwp) ::  i
@@ -458,4 +464,5 @@
     REAL(wp) ::  near_facade_temperature
     REAL(wp) ::  phi_hc_nd_dummy
+    REAL(wp), INTENT(IN) :: theta_m_t_prev
 !
 !-- Calculation of total mass specific thermal load (internal and external)
@@ -468,5 +475,5 @@
                )                                                                !< [degree_C] Eq. (C.5)
 !
-!-- Calculation of component temperature at factual timestep
+!-- Calculation of component temperature at current timestep
     theta_m_t = ( ( theta_m_t_prev                                                                 &
                     * ( ( c_m / 3600.0_wp ) - 0.5_wp * ( h_t_3 + h_t_wm ) )                        &
@@ -476,9 +483,9 @@
                 )                                                               !< [degree_C] Eq. (C.4)
 !
-!-- Calculation of mean inner temperature for the RC-node in actual timestep
+!-- Calculation of mean inner temperature for the RC-node in current 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
+!-- Calculation of mean surface temperature of the RC-node in current timestep
     theta_s = ( (   h_t_ms * theta_m + phi_st + h_t_es * pt(k,j,i)                                 &
                   + h_t_1  * ( near_facade_temperature                                             &
@@ -490,6 +497,9 @@
 !
 !-- 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
@@ -925,4 +935,7 @@
           ALLOCATE( buildings(nb)%l_v(1:buildings(nb)%num_facades_per_building_v_l) )
           ALLOCATE( buildings(nb)%m_v(1:buildings(nb)%num_facades_per_building_v_l) )
+
+          ALLOCATE( buildings(nb)%theta_m_t_prev_h(1:buildings(nb)%num_facades_per_building_h_l) )
+          ALLOCATE( buildings(nb)%theta_m_t_prev_v(1:buildings(nb)%num_facades_per_building_v_l) )
        ENDIF
 
@@ -1125,11 +1138,15 @@
     ENDDO
 !
-!-- Initial room temperature [K]
-!-- (after first loop, use theta_m_t as theta_m_t_prev)
-    theta_m_t_prev = initial_indoor_temperature
-!
 !-- Initialize indoor temperature. Actually only for output at initial state.
     DO  nb = 1, num_build
-       IF ( buildings(nb)%on_pe )  buildings(nb)%t_in(:) = initial_indoor_temperature
+       IF ( buildings(nb)%on_pe )  THEN
+          buildings(nb)%t_in(:) = initial_indoor_temperature
+
+!
+!--       (after first loop, use theta_m_t as theta_m_t_prev)
+          buildings(nb)%theta_m_t_prev_h(:) = initial_indoor_temperature
+          buildings(nb)%theta_m_t_prev_v(:) = initial_indoor_temperature
+
+       ENDIF
     ENDDO
 
@@ -1337,5 +1354,5 @@
 
              CALL  im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,                 &
-                                          near_facade_temperature, phi_hc_nd )
+                                          near_facade_temperature, phi_hc_nd, buildings(nb)%theta_m_t_prev_h(fa) )
 !
 !--          If air temperature between border temperatures of heating and cooling, assign output
@@ -1366,5 +1383,5 @@
 
                 CALL  im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,              &
-                                             near_facade_temperature, phi_hc_nd )
+                                             near_facade_temperature, phi_hc_nd, buildings(nb)%theta_m_t_prev_h(fa) )
                 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)                        &
@@ -1395,10 +1412,11 @@
 !--             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 )
+                                             near_facade_temperature, phi_hc_nd, buildings(nb)%theta_m_t_prev_h(fa) )
                 theta_air_ac = theta_air
              ENDIF
 !
 !--          Update theta_m_t_prev
-             theta_m_t_prev = theta_m_t
+             buildings(nb)%theta_m_t_prev_h(fa) = theta_m_t
+
 
              q_vent = h_v * ( theta_air - near_facade_temperature )
@@ -1408,4 +1426,5 @@
 !--          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(l)%t_indoor(m) = theta_op                               !< not integrated now
 !
@@ -1440,5 +1459,5 @@
                             ) / facade_element_area                                             !< [W/m2] , observe the directional
                                                                                                 !< convention in PALM!
-             surf_usm_h(l)%waste_heat(m) = q_waste_heat
+             surf_usm_h(l)%waste_heat(m) = 0.0_wp !q_waste_heat
           ENDDO !< Horizontal surfaces loop
 !
@@ -1524,7 +1543,8 @@
              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(m,ind_veg_wall) * t_wall_v(l)%val(nzt_wall,m) &
+             indoor_wall_window_temperature =                                                          &
+                                    surf_usm_v(l)%frac(m,ind_veg_wall) * t_wall_v(l)%val(nzt_wall,m)   &
                                   + surf_usm_v(l)%frac(m,ind_wat_win)  * t_window_v(l)%val(nzt_wall,m)
+
 !
 !--          Solar thermal gains. If net_sw_in larger than sun-protection
@@ -1563,5 +1583,5 @@
              phi_hc_nd = 0.0_wp
              CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,                  &
-                                         near_facade_temperature, phi_hc_nd )
+                                         near_facade_temperature, phi_hc_nd, buildings(nb)%theta_m_t_prev_v(fa) )
 !
 !--          If air temperature between border temperatures of heating and cooling, assign output
@@ -1592,5 +1612,5 @@
 
                 CALL  im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,              &
-                                             near_facade_temperature, phi_hc_nd )
+                                             near_facade_temperature, phi_hc_nd, buildings(nb)%theta_m_t_prev_v(fa) )
 
                 theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
@@ -1623,10 +1643,11 @@
 !--             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 )
+                                             near_facade_temperature, phi_hc_nd, buildings(nb)%theta_m_t_prev_v(fa) )
                 theta_air_ac = theta_air
              ENDIF
 !
 !--          Update theta_m_t_prev
-             theta_m_t_prev = theta_m_t
+             buildings(nb)%theta_m_t_prev_v(fa) = theta_m_t
+
 
              q_vent = h_v * ( theta_air - near_facade_temperature )
@@ -1635,4 +1656,5 @@
 !--          Will be used for thermal comfort calculations.
              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
 !
@@ -1665,5 +1687,5 @@
                                                     ) / facade_element_area  !< [W/m2] , observe the directional convention in
                                                                              !< PALM!
-             surf_usm_v(l)%waste_heat(m) = q_waste_heat
+             surf_usm_v(l)%waste_heat(m) = 0.0_wp !q_waste_heat
           ENDDO !< Vertical surfaces loop
        ENDIF !< buildings(nb)%on_pe
