Index: palm/trunk/SOURCE/indoor_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/indoor_model_mod.f90	(revision 4346)
+++ palm/trunk/SOURCE/indoor_model_mod.f90	(revision 4402)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! Major bugfix in calculation of energy demand - floor-area-per-facade was wrongly
+! calculated leading to unrealistically high energy demands and thus to
+! unreallistically high waste-heat fluxes.
+! 
+! 4346 2019-12-18 11:55:56Z motisi
 ! Introduction of wall_flags_total_0, which currently sets bits based on static
 ! topography information used in wall_flags_static_0
@@ -133,4 +138,9 @@
  MODULE indoor_model_mod 
 
+    USE arrays_3d,                                                             &
+        ONLY:  ddzw,                                                           &
+               dzw,                                                            &
+               pt
+
     USE control_parameters,                                                    &
         ONLY:  initializing_actions
@@ -188,4 +198,5 @@
                                           !< (very light), 110000 (light), 165000 (medium), 260000 (heavy), 370000 (very heavy)
        REAL(wp) ::  f_c_win               !< [-] shading factor
+       REAL(wp) ::  fapf                  !< [m2/m2] floor area per facade
        REAL(wp) ::  g_value_win           !< [-] SHGC factor
        REAL(wp) ::  h_es                  !< [W/(m2 K)] surface-related heat transfer coefficient between extern and surface 
@@ -408,12 +419,5 @@
  SUBROUTINE im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,    &
                                    near_facade_temperature, phi_hc_nd_dummy )
- 
-    USE arrays_3d,                                                             &
-        ONLY:  pt
-    
-    
-    IMPLICIT NONE
-    
-    
+
     INTEGER(iwp) ::  i
     INTEGER(iwp) ::  j
@@ -477,7 +481,4 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE im_init
-  
-    USE arrays_3d,                                                             &
-        ONLY:  dzw
 
     USE control_parameters,                                                    &
@@ -497,6 +498,4 @@
     USE urban_surface_mod,                                                     &
         ONLY:  building_pars, building_type
-
-    IMPLICIT NONE
 
     INTEGER(iwp) ::  bt          !< local building type
@@ -908,32 +907,9 @@
           ALLOCATE( buildings(nb)%m_v(1:buildings(nb)%num_facades_per_building_v_l) )
        ENDIF
-!
-!--    Determine volume per facade element (vpf)
+
        IF ( buildings(nb)%on_pe )  THEN
           ALLOCATE( buildings(nb)%vpf(buildings(nb)%kb_min:buildings(nb)%kb_max) )
           buildings(nb)%vpf = 0.0_wp
-          
-          DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
-!
-!--          In order to avoid division by zero, check if the number of facade
-!--          elements is /= 0. This can e.g. happen if a building is embedded
-!--          in higher terrain and at a given k-level neither horizontal nor
-!--          vertical facade elements are located.
-             IF ( buildings(nb)%num_facade_h(k)                                &
-                + buildings(nb)%num_facade_v(k) > 0 )  THEN  
-                buildings(nb)%vpf(k) = buildings(nb)%volume(k) /               &
-                                REAL( buildings(nb)%num_facade_h(k) +          &
-                                      buildings(nb)%num_facade_v(k), KIND = wp )
-             ENDIF
-          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
-         
+
           facade_area_v = 0.0_wp          
           DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
@@ -941,8 +917,18 @@
                              * dzw(k+1) * dx
           ENDDO
-         
+!
+!--       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!
           buildings(nb)%vpf = buildings(nb)%vol_tot /                          &
                         ( buildings(nb)%num_facades_per_building_h * dx * dy + &
                           facade_area_v )
+!
+!--       Determine floor-area-per-facade.
+          buildings(nb)%fapf = buildings(nb)%num_facades_per_building_h        &
+                             * dx * dy                                         &
+                             / ( buildings(nb)%num_facades_per_building_h      &
+                               * dx * dy + facade_area_v )
        ENDIF
     ENDDO
@@ -1146,7 +1132,4 @@
  SUBROUTINE im_main_heatcool
 
-    USE arrays_3d,                                                             &
-        ONLY:  ddzw, dzw
-
 !     USE basic_constants_and_equations_mod,                                     &
 !         ONLY:  c_p
@@ -1166,7 +1149,4 @@
         ONLY:  nzt_wall, t_wall_h, t_wall_v, t_window_h, t_window_v,           &
                building_type
-
-
-    IMPLICIT NONE
 
     INTEGER(iwp) ::  i    !< index of facade-adjacent atmosphere grid point in x-direction
@@ -1237,9 +1217,9 @@
 !--          Determine building height level index. 
              kk = surf_usm_h%k(m) + surf_usm_h%koff
-!            
+!
 !--          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+1)                        !< [m2/m2] floor area per facade area
-             indoor_volume_per_facade     = buildings(nb)%vpf(kk)                                     !< [m3/m2] indoor air volume per facade area 
+             facade_element_area          = dx * dy                               !< [m2] surface area per facade element    
+             floor_area_per_facade        = buildings(nb)%fapf                    !< [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 +              &
@@ -1447,12 +1427,13 @@
 !--          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+1)                          !< [m2] surface area per facade element
-             IF ( l == 2  .OR. l == 3 ) facade_element_area = dy * dzw(kk+1)                          !< [m2] surface area per facade element
-             floor_area_per_facade        = buildings(nb)%vpf(kk) * ddzw(kk+1)                        !< [m2/m2] floor area per facade area
-             indoor_volume_per_facade     = buildings(nb)%vpf(kk)                                     !< [m3/m2] indoor air volume per facade area 
+             IF ( l == 0  .OR. l == 1 ) facade_element_area = dx * dzw(kk+1)    !< [m2] surface area per facade element
+             IF ( l == 2  .OR. l == 3 ) facade_element_area = dy * dzw(kk+1)    !< [m2] surface area per facade element
+
+             floor_area_per_facade        = buildings(nb)%fapf                  !< [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
+             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 )
@@ -1571,5 +1552,4 @@
                 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)
@@ -1642,7 +1622,5 @@
                             ) / facade_element_area !< [W/m2] , observe the directional convention in PALM!
              surf_usm_v(l)%waste_heat(m) = q_waste_heat
-             
           ENDDO !< Vertical surfaces loop
-
        ENDIF !< buildings(nb)%on_pe
     ENDDO !< buildings loop
@@ -1724,7 +1702,5 @@
 !-----------------------------------------------------------------------------! 
  SUBROUTINE im_check_data_output( var, unit )
-        
-    IMPLICIT NONE
-    
+
     CHARACTER (LEN=*) ::  unit   !<
     CHARACTER (LEN=*) ::  var    !<
@@ -1771,6 +1747,4 @@
 !   USE control_parameters,
 !       ONLY: message_string
-       
-   IMPLICIT NONE
    
  END SUBROUTINE im_check_parameters
@@ -1784,5 +1758,4 @@
  SUBROUTINE im_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
 
-   IMPLICIT NONE
    
    CHARACTER (LEN=*), INTENT(IN)  ::  var
@@ -1831,10 +1804,8 @@
                                nzb_do, nzt_do )
 
-   USE indices
-   
-   USE kinds
-          
-   IMPLICIT NONE
-   
+    USE indices
+
+    USE kinds
+
     CHARACTER (LEN=*) ::  variable !< 
 
@@ -1971,6 +1942,5 @@
     USE control_parameters,                                                    &
         ONLY:  indoor_model
-    
-    IMPLICIT NONE
+
 
     CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
