source: palm/trunk/SOURCE/indoor_model_mod.f90 @ 4180

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

  • Property svn:keywords set to Id
File size: 92.6 KB
Line 
1!> @file indoor_model_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2018-2019 Leibniz Universitaet Hannover
18! Copyright 2018-2019 Hochschule Offenburg
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: indoor_model_mod.f90
28! Bugfix in case of non grid-resolved buildings. Further, vertical grid spacing
29! is now considered at the correct level. 
30!
31!
32! - change calculation of a_m and c_m
33! - change calculation of u-values (use h_es in building array)
34! - rename h_tr_... to  h_t_...
35!          h_tr_em  to  h_t_wm
36!          h_tr_op  to  h_t_wall
37!          h_tr_w   to  h_t_es
38! - rename h_ve     to  h_v
39! - rename h_is     to  h_ms
40! - inserted net_floor_area
41! - inserted params_waste_heat_h, params_waste_heat_c from building database
42!   in building array
43! - change calculation of q_waste_heat
44! - bugfix in averaging mean indoor temperature
45!
46!
47! 3759 2019-02-21 15:53:45Z suehring $
48! - Calculation of total building volume
49! - Several bugfixes
50! - Calculation of building height revised
51!
52! 3745 2019-02-15 18:57:56Z suehring
53! - remove building_type from module
54! - initialize parameters for each building individually instead of a bulk
55!   initializaion with  identical building type for all
56! - output revised
57! - add missing _wp
58! - some restructuring of variables in building data structure
59!
60! 3744 2019-02-15 18:38:58Z suehring
61! Some interface calls moved to module_interface + cleanup
62!
63!
64! Description:
65! ------------
66!> <Description of the new module>
67!> Module for Indoor Climate Model (ICM)
68!> The module is based on the DIN EN ISO 13790 with simplified hour-based procedure.
69!> This model is a equivalent circuit diagram of a three-point RC-model (5R1C).
70!> This module differ between indoor-air temperature an average temperature of indoor surfaces which make it prossible to determine thermal comfort
71!> the heat transfer between indoor and outdoor is simplified
72
73!> @todo Replace window_area_per_facade by %frac(1,m) for window
74!> @todo emissivity change for window blinds if solar_protection_on=1
75
76!> @note Do we allow use of integer flags, or only logical flags? (concerns e.g. cooling_on, heating_on)
77!> @note How to write indoor temperature output to pt array?
78!>
79!> @bug  <Enter known bugs here>
80!------------------------------------------------------------------------------!
81 MODULE indoor_model_mod
82
83    USE control_parameters,                                                    &
84        ONLY:  initializing_actions
85
86    USE kinds
87   
88    USE netcdf_data_input_mod,                                                 &
89        ONLY:  building_id_f, building_type_f
90
91    USE surface_mod,                                                           &
92        ONLY:  surf_usm_h, surf_usm_v
93
94
95    IMPLICIT NONE
96
97!
98!-- Define data structure for buidlings.
99    TYPE build
100
101       INTEGER(iwp) ::  id                                !< building ID
102       INTEGER(iwp) ::  kb_min                            !< lowest vertical index of a building
103       INTEGER(iwp) ::  kb_max                            !< highest vertical index of a building
104       INTEGER(iwp) ::  num_facades_per_building_h = 0    !< total number of horizontal facades elements
105       INTEGER(iwp) ::  num_facades_per_building_h_l = 0  !< number of horizontal facade elements on local subdomain
106       INTEGER(iwp) ::  num_facades_per_building_v = 0    !< total number of vertical facades elements
107       INTEGER(iwp) ::  num_facades_per_building_v_l = 0  !< number of vertical facade elements on local subdomain
108       INTEGER(iwp) ::  ventilation_int_loads             !< [-] allocation of activity in the building
109
110       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  l_v            !< index array linking surface-element orientation index
111                                                                  !< for vertical surfaces with building
112       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  m_h            !< index array linking surface-element index for
113                                                                  !< horizontal surfaces with building
114       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  m_v            !< index array linking surface-element index for
115                                                                  !< vertical surfaces with building
116       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  num_facade_h   !< number of horizontal facade elements per buidling
117                                                                  !< and height level
118       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  num_facade_v   !< number of vertical facades elements per buidling
119                                                                  !< and height level
120                                                                 
121
122       LOGICAL ::  on_pe = .FALSE.   !< flag indicating whether a building with certain ID is on local subdomain
123       
124       REAL(wp) ::  air_change_high       !< [1/h] air changes per time_utc_hour
125       REAL(wp) ::  air_change_low        !< [1/h] air changes per time_utc_hour
126       REAL(wp) ::  area_facade           !< [m2] area of total facade
127       REAL(wp) ::  building_height       !< building height
128       REAL(wp) ::  eta_ve                !< [-] heat recovery efficiency
129       REAL(wp) ::  factor_a              !< [-] Dynamic parameters specific effective surface according to Table 12; 2.5
130                                          !< (very light, light and medium), 3.0 (heavy), 3.5 (very heavy)
131       REAL(wp) ::  factor_c              !< [J/(m2 K)] Dynamic parameters inner heatstorage according to Table 12; 80000
132                                          !< (very light), 110000 (light), 165000 (medium), 260000 (heavy), 370000 (very heavy)
133       REAL(wp) ::  f_c_win               !< [-] shading factor
134       REAL(wp) ::  g_value_win           !< [-] SHGC factor
135       REAL(wp) ::  h_es                  !< [W/(m2 K)] surface-related heat transfer coefficient between extern and surface
136       REAL(wp) ::  height_cei_con        !< [m] ceiling construction heigth
137       REAL(wp) ::  height_storey         !< [m] storey heigth
138       REAL(wp) ::  params_waste_heat_c   !< [-] anthropogenic heat outputs for cooling e.g. 1.33 for KKM with COP = 3
139       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
140       REAL(wp) ::  phi_c_max             !< [W] Max. Cooling capacity (negative)
141       REAL(wp) ::  phi_h_max             !< [W] Max. Heating capacity (positive)
142       REAL(wp) ::  q_c_max               !< [W/m2] Max. Cooling heat flux per netto floor area (negative)
143       REAL(wp) ::  q_h_max               !< [W/m2] Max. Heating heat flux per netto floor area (positive)
144       REAL(wp) ::  qint_high             !< [W/m2] internal heat gains, option Database qint_0-23
145       REAL(wp) ::  qint_low              !< [W/m2] internal heat gains, option Database qint_0-23
146       REAL(wp) ::  lambda_at             !< [-] ratio internal surface/floor area chap. 7.2.2.2.
147       REAL(wp) ::  lambda_layer3         !< [W/(m*K)] Thermal conductivity of the inner layer
148       REAL(wp) ::  net_floor_area        !< [m2] netto ground area
149       REAL(wp) ::  s_layer3              !< [m] half thickness of the inner layer (layer_3)
150       REAL(wp) ::  theta_int_c_set       !< [degree_C] Max. Setpoint temperature (summer)
151       REAL(wp) ::  theta_int_h_set       !< [degree_C] Max. Setpoint temperature (winter)
152       REAL(wp) ::  u_value_win           !< [W/(m2*K)] transmittance
153       REAL(wp) ::  vol_tot               !< [m3] total building volume
154
155       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in       !< mean building indoor temperature, height dependent
156       REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in_l     !< mean building indoor temperature on local subdomain, height dependent
157       REAL(wp), DIMENSION(:), ALLOCATABLE ::  volume     !< total building volume, height dependent
158       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vol_frac   !< fraction of local on total building volume, height dependent
159       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vpf        !< building volume volume per facade element, height dependent
160       
161    END TYPE build
162
163    TYPE(build), DIMENSION(:), ALLOCATABLE ::  buildings   !< building array
164
165    INTEGER(iwp) ::  num_build   !< total number of buildings in domain
166!
167!-- Declare all global variables within the module
168    INTEGER(iwp) ::  cooling_on              !< Indoor cooling flag (0=off, 1=on)
169    INTEGER(iwp) ::  heating_on              !< Indoor heating flag (0=off, 1=on)
170    INTEGER(iwp) ::  solar_protection_off    !< Solar protection off
171    INTEGER(iwp) ::  solar_protection_on     !< Solar protection on
172
173    REAL(wp) ::  a_m                           !< [m2] the effective mass-related area
174    REAL(wp) ::  air_change                    !< [1/h] Airflow
175    REAL(wp) ::  c_m                           !< [J/K] internal heat storage capacity
176    REAL(wp) ::  dt_indoor = 3600.0_wp         !< [s] namelist parameter: time interval for indoor-model application
177    REAL(wp) ::  facade_element_area           !< [m2_facade] building surface facade
178    REAL(wp) ::  floor_area_per_facade         !< [m2/m2] floor area per facade area
179    REAL(wp) ::  h_t_1                         !< [W/K] Heat transfer coefficient auxiliary variable 1
180    REAL(wp) ::  h_t_2                         !< [W/K] Heat transfer coefficient auxiliary variable 2
181    REAL(wp) ::  h_t_3                         !< [W/K] Heat transfer coefficient auxiliary variable 3
182    REAL(wp) ::  h_t_wm                        !< [W/K] Heat transfer coefficient of the emmision (got with h_t_ms the thermal mass)
183    REAL(wp) ::  h_t_is                        !< [W/K] thermal coupling conductance (Thermischer Kopplungsleitwert)
184    REAL(wp) ::  h_t_ms                        !< [W/K] Heat transfer conductance term (got with h_t_wm the thermal mass)
185    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
186    REAL(wp) ::  h_t_es                        !< [W/K] heat transfer coefficient of doors, windows, curtain walls and glazed walls (assumption: thermal mass=0)
187    REAL(wp) ::  h_v                           !< [W/K] heat transfer of ventilation
188    REAL(wp) ::  indoor_volume_per_facade      !< [m3] indoor air volume per facade element
189    REAL(wp) ::  initial_indoor_temperature    !< [K] initial indoor temperature (namelist parameter)
190    REAL(wp) ::  net_sw_in                     !< [W/m2] net short-wave radiation
191    REAL(wp) ::  phi_hc_nd                     !< [W] heating demand and/or cooling demand
192    REAL(wp) ::  phi_hc_nd_10                  !< [W] heating demand and/or cooling demand for heating or cooling
193    REAL(wp) ::  phi_hc_nd_ac                  !< [W] actual heating demand and/or cooling demand
194    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)
195    REAL(wp) ::  phi_ia                        !< [W] internal air load = internal loads * 0.5, Eq. (C.1)
196    REAL(wp) ::  phi_m                         !< [W] mass specific thermal load (internal and external)
197    REAL(wp) ::  phi_mtot                      !< [W] total mass specific thermal load (internal and external)
198    REAL(wp) ::  phi_sol                       !< [W] solar loads
199    REAL(wp) ::  phi_st                        !< [W] mass specific thermal load implied non thermal mass
200    REAL(wp) ::  q_wall_win                    !< [W/m2]heat flux from indoor into wall/window
201    REAL(wp) ::  q_waste_heat                  !< [W/m2]waste heat, sum of waste heat over the roof to Palm
202   
203    REAL(wp) ::  q_c_m                         !< [W] Energy of thermal storage mass specific thermal load for internal and external heatsources (for energy bilanz)
204    REAL(wp) ::  q_c_st                        !< [W] Energy of thermal storage mass specific thermal load implied non thermal mass (for energy bilanz)
205    REAL(wp) ::  q_int                         !< [W] Energy of internal air load (for energy bilanz)
206    REAL(wp) ::  q_sol                         !< [W] Energy of solar (for energy bilanz)
207    REAL(wp) ::  q_trans                       !< [W] Energy of transmission (for energy bilanz)
208    REAL(wp) ::  q_vent                        !< [W] Energy of ventilation (for energy bilanz)
209   
210    REAL(wp) ::  schedule_d                    !< [-] activation for internal loads (low or high + low)
211    REAL(wp) ::  skip_time_do_indoor = 0.0_wp  !< [s] Indoor model is not called before this time
212    REAL(wp) ::  theta_air                     !< [degree_C] air temperature of the RC-node
213    REAL(wp) ::  theta_air_0                   !< [degree_C] air temperature of the RC-node in equilibrium
214    REAL(wp) ::  theta_air_10                  !< [degree_C] air temperature of the RC-node from a heating capacity of 10 W/m2
215    REAL(wp) ::  theta_air_ac                  !< [degree_C] actual room temperature after heating/cooling
216    REAL(wp) ::  theta_air_set                 !< [degree_C] Setpoint_temperature for the room
217    REAL(wp) ::  theta_m                       !< [degree_C} inner temperature of the RC-node
218    REAL(wp) ::  theta_m_t                     !< [degree_C] (Fictive) component temperature timestep
219    REAL(wp) ::  theta_m_t_prev                !< [degree_C] (Fictive) component temperature previous timestep (do not change)
220    REAL(wp) ::  theta_op                      !< [degree_C] operative temperature
221    REAL(wp) ::  theta_s                       !< [degree_C] surface temperature of the RC-node
222    REAL(wp) ::  time_indoor = 0.0_wp          !< [s] time since last call of indoor model
223    REAL(wp) ::  total_area                    !< [m2] area of all surfaces pointing to zone
224    REAL(wp) ::  window_area_per_facade        !< [m2] window area per facade element
225   
226    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)
227    REAL(wp), PARAMETER ::  h_ms                     = 9.1_wp      !< [W/(m2 K)] surface-related heat transfer coefficient between component and surface (chap. 12.2.2)
228    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
229    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)
230    REAL(wp), PARAMETER ::  params_f_win             = 0.5_wp      !< [-] proportion of window area, Database A_win aus Datenbank 27 window_area_per_facade_percent
231    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
232
233   
234    SAVE
235
236
237    PRIVATE
238   
239!
240!-- Add INTERFACES that must be available to other modules
241    PUBLIC im_init, im_main_heatcool, im_parin, im_define_netcdf_grid,          &
242           im_check_data_output, im_data_output_3d, im_check_parameters
243   
244
245!
246!-- Add VARIABLES that must be available to other modules
247    PUBLIC dt_indoor, skip_time_do_indoor, time_indoor
248
249!
250!-- PALM interfaces:
251!-- Data output checks for 2D/3D data to be done in check_parameters
252     INTERFACE im_check_data_output
253        MODULE PROCEDURE im_check_data_output
254     END INTERFACE im_check_data_output
255!
256!-- Input parameter checks to be done in check_parameters
257     INTERFACE im_check_parameters
258        MODULE PROCEDURE im_check_parameters
259     END INTERFACE im_check_parameters
260!
261!-- Data output of 3D data
262     INTERFACE im_data_output_3d
263        MODULE PROCEDURE im_data_output_3d
264     END INTERFACE im_data_output_3d
265
266!
267!-- Definition of data output quantities
268     INTERFACE im_define_netcdf_grid
269        MODULE PROCEDURE im_define_netcdf_grid
270     END INTERFACE im_define_netcdf_grid
271!
272! !
273! !-- Output of information to the header file
274!     INTERFACE im_header
275!        MODULE PROCEDURE im_header
276!     END INTERFACE im_header
277!
278!-- Calculations for indoor temperatures 
279    INTERFACE im_calc_temperatures
280       MODULE PROCEDURE im_calc_temperatures
281    END INTERFACE im_calc_temperatures
282!
283!-- Initialization actions 
284    INTERFACE im_init
285       MODULE PROCEDURE im_init
286    END INTERFACE im_init
287!
288!-- Main part of indoor model 
289    INTERFACE im_main_heatcool
290       MODULE PROCEDURE im_main_heatcool
291    END INTERFACE im_main_heatcool
292!
293!-- Reading of NAMELIST parameters
294    INTERFACE im_parin
295       MODULE PROCEDURE im_parin
296    END INTERFACE im_parin
297
298 CONTAINS
299
300!------------------------------------------------------------------------------!
301! Description:
302! ------------
303!< Calculation of the air temperatures and mean radiation temperature
304!< This is basis for the operative temperature
305!< Based on a Crank-Nicholson scheme with a timestep of a hour
306!------------------------------------------------------------------------------!
307 SUBROUTINE im_calc_temperatures ( i, j, k, indoor_wall_window_temperature,    &
308                                   near_facade_temperature, phi_hc_nd_dummy )
309 
310    USE arrays_3d,                                                             &
311        ONLY:  pt
312   
313   
314    IMPLICIT NONE
315   
316   
317    INTEGER(iwp) ::  i
318    INTEGER(iwp) ::  j
319    INTEGER(iwp) ::  k
320   
321    REAL(wp) ::  indoor_wall_window_temperature  !< weighted temperature of innermost wall/window layer
322    REAL(wp) ::  near_facade_temperature
323    REAL(wp) ::  phi_hc_nd_dummy
324!
325!-- Calculation of total mass specific thermal load (internal and external)
326    phi_mtot = ( phi_m + h_t_wm * indoor_wall_window_temperature               &
327                       + h_t_3  * ( phi_st + h_t_es * pt(k,j,i)                &
328                                            + h_t_1 *                          &
329                                    ( ( ( phi_ia + phi_hc_nd_dummy ) / h_v )   &
330                                                 + near_facade_temperature )   &
331                                   ) / h_t_2                                   &
332               )                                                                !< [degree_C] Eq. (C.5)
333!
334!-- Calculation of component temperature at factual timestep
335    theta_m_t = ( ( theta_m_t_prev                                             &
336                    * ( ( c_m / 3600.0_wp ) - 0.5_wp * ( h_t_3 + h_t_wm ) )    &
337                     + phi_mtot                                                &
338                  )                                                            &
339                  /   ( ( c_m / 3600.0_wp ) + 0.5_wp * ( h_t_3 + h_t_wm ) )    &
340                )                                                               !< [degree_C] Eq. (C.4)
341!
342!-- Calculation of mean inner temperature for the RC-node in actual timestep
343    theta_m = ( theta_m_t + theta_m_t_prev ) * 0.5_wp                           !< [degree_C] Eq. (C.9)
344   
345!
346!-- Calculation of mean surface temperature of the RC-node in actual timestep
347    theta_s = ( (   h_t_ms * theta_m + phi_st + h_t_es * pt(k,j,i)             &
348                  + h_t_1  * ( near_facade_temperature                         &
349                           + ( phi_ia + phi_hc_nd_dummy ) / h_v )              &
350                )                                                              &
351                / ( h_t_ms + h_t_es + h_t_1 )                                  &
352              )                                                                 !< [degree_C] Eq. (C.10)
353   
354!
355!-- Calculation of the air temperature of the RC-node
356    theta_air = ( h_t_is * theta_s + h_v * near_facade_temperature             &
357                + phi_ia + phi_hc_nd_dummy ) / ( h_t_is + h_v )                 !< [degree_C] Eq. (C.11)
358
359 END SUBROUTINE im_calc_temperatures
360
361!------------------------------------------------------------------------------!
362! Description:
363! ------------
364!> Initialization of the indoor model.
365!> Static information are calculated here, e.g. building parameters and
366!> geometrical information, everything that doesn't change in time.
367!
368!-- Input values
369!-- Input datas from Palm, M4
370!     i_global             -->  net_sw_in                         !< global radiation [W/m2]
371!     theta_e              -->  pt(k,j,i)                         !< undisturbed outside temperature, 1. PALM volume, for windows
372!     theta_sup = theta_f  -->  surf_usm_h%pt_10cm(m)
373!                               surf_usm_v(l)%pt_10cm(m)          !< Air temperature, facade near (10cm) air temperature from 1. Palm volume
374!     theta_node           -->  t_wall_h(nzt_wall,m)
375!                               t_wall_v(l)%t(nzt_wall,m)         !< Temperature of innermost wall layer, for opaque wall
376!------------------------------------------------------------------------------!
377 SUBROUTINE im_init
378 
379    USE arrays_3d,                                                             &
380        ONLY:  dzw
381
382    USE control_parameters,                                                    &
383        ONLY:  message_string
384
385    USE indices,                                                               &
386        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
387
388    USE grid_variables,                                                        &
389        ONLY:  dx, dy
390
391    USE pegrid
392
393    USE surface_mod,                                                           &
394        ONLY:  surf_usm_h, surf_usm_v
395       
396    USE urban_surface_mod,                                                     &
397        ONLY:  building_pars, building_type
398
399    IMPLICIT NONE
400
401    INTEGER(iwp) ::  bt   !< local building type
402    INTEGER(iwp) ::  i    !< running index along x-direction
403    INTEGER(iwp) ::  j    !< running index along y-direction
404    INTEGER(iwp) ::  k    !< running index along z-direction
405    INTEGER(iwp) ::  l    !< running index for surface-element orientation
406    INTEGER(iwp) ::  m    !< running index surface elements
407    INTEGER(iwp) ::  n    !< building index
408    INTEGER(iwp) ::  nb   !< building index
409
410    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids           !< building IDs on entire model domain
411    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final     !< building IDs on entire model domain,
412                                                                    !< multiple occurences are sorted out
413    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_final_tmp !< temporary array used for resizing
414    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l         !< building IDs on local subdomain
415    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  build_ids_l_tmp     !< temporary array used to resize array of building IDs
416    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displace_dum        !< displacements of start addresses, used for MPI_ALLGATHERV
417    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_max_l             !< highest vertical index of a building on subdomain
418    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k_min_l             !< lowest vertical index of a building on subdomain
419    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  n_fa                !< counting array
420    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  num_facades_h       !< dummy array used for summing-up total number of
421                                                                    !< horizontal facade elements
422    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  num_facades_v       !< dummy array used for summing-up total number of
423                                                                    !< vertical facade elements
424    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  receive_dum_h       !< dummy array used for MPI_ALLREDUCE 
425    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  receive_dum_v       !< dummy array used for MPI_ALLREDUCE 
426   
427    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings         !< number of buildings with different ID on entire model domain
428    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  num_buildings_l       !< number of buildings with different ID on local subdomain
429                                                             
430    REAL(wp) ::  u_tmp                                     !< dummy for temporary calculation of u-value without h_is
431    REAL(wp) ::  du_tmp                                    !< 1/u_tmp
432    REAL(wp) ::  du_win_tmp                                !< 1/building(nb)%u_value_win
433    REAL(wp) ::  facade_area_v                             !< dummy to compute the total facade area from vertical walls
434
435    REAL(wp), DIMENSION(:), ALLOCATABLE ::  volume         !< total building volume at each discrete height level
436    REAL(wp), DIMENSION(:), ALLOCATABLE ::  volume_l       !< total building volume at each discrete height level,
437                                                           !< on local subdomain
438
439    CALL location_message( 'initializing indoor model', .FALSE. )
440!
441!-- Initializing of indoor model is only possible if buildings can be
442!-- distinguished by their IDs.
443    IF ( .NOT. building_id_f%from_file )  THEN
444       message_string = 'Indoor model requires information about building_id'
445       CALL message( 'im_init', 'PA0999', 1, 2, 0, 6, 0  )
446    ENDIF
447!
448!-- Determine number of different building IDs on local subdomain.
449    num_buildings_l = 0
450    num_buildings   = 0
451    ALLOCATE( build_ids_l(1) )
452    DO  i = nxl, nxr
453       DO  j = nys, nyn
454          IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
455             IF ( num_buildings_l(myid) > 0 )  THEN
456                IF ( ANY( building_id_f%var(j,i) .EQ.  build_ids_l ) )  THEN
457                   CYCLE
458                ELSE
459                   num_buildings_l(myid) = num_buildings_l(myid) + 1
460!
461!--                Resize array with different local building ids
462                   ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
463                   build_ids_l_tmp = build_ids_l
464                   DEALLOCATE( build_ids_l )
465                   ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
466                   build_ids_l(1:num_buildings_l(myid)-1) =                 &
467                               build_ids_l_tmp(1:num_buildings_l(myid)-1)
468                   build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
469                   DEALLOCATE( build_ids_l_tmp )
470                ENDIF
471!
472!--          First occuring building id on PE
473             ELSE
474                num_buildings_l(myid) = num_buildings_l(myid) + 1
475                build_ids_l(1) = building_id_f%var(j,i)
476             ENDIF
477          ENDIF
478       ENDDO
479    ENDDO
480!
481!-- Determine number of building IDs for the entire domain. (Note, building IDs
482!-- can appear multiple times as buildings might be distributed over several
483!-- PEs.)
484#if defined( __parallel )
485    CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs,              &
486                        MPI_INTEGER, MPI_SUM, comm2d, ierr ) 
487#else
488    num_buildings = num_buildings_l
489#endif
490    ALLOCATE( build_ids(1:SUM(num_buildings)) )
491!
492!-- Gather building IDs. Therefore, first, determine displacements used
493!-- required for MPI_GATHERV call.
494    ALLOCATE( displace_dum(0:numprocs-1) )
495    displace_dum(0) = 0
496    DO i = 1, numprocs-1
497       displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
498    ENDDO
499
500#if defined( __parallel )
501    CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)),                 &
502                         num_buildings(myid),                                  &
503                         MPI_INTEGER,                                          &
504                         build_ids,                                            &
505                         num_buildings,                                        &
506                         displace_dum,                                         & 
507                         MPI_INTEGER,                                          &
508                         comm2d, ierr )   
509
510    DEALLOCATE( displace_dum )
511
512#else
513    build_ids = build_ids_l
514#endif
515
516!
517!-- Note: in parallel mode, building IDs can occur mutliple times, as
518!-- each PE has send its own ids. Therefore, sort out building IDs which
519!-- appear multiple times.
520    num_build = 0
521    DO  n = 1, SIZE(build_ids)
522
523       IF ( ALLOCATED(build_ids_final) )  THEN
524          IF ( ANY( build_ids(n) == build_ids_final ) )  THEN
525             CYCLE
526          ELSE
527             num_build = num_build + 1
528!
529!--          Resize
530             ALLOCATE( build_ids_final_tmp(1:num_build) )
531             build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
532             DEALLOCATE( build_ids_final )
533             ALLOCATE( build_ids_final(1:num_build) )
534             build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
535             build_ids_final(num_build) = build_ids(n)
536             DEALLOCATE( build_ids_final_tmp )
537          ENDIF             
538       ELSE
539          num_build = num_build + 1
540          ALLOCATE( build_ids_final(1:num_build) )
541          build_ids_final(num_build) = build_ids(n)
542       ENDIF
543    ENDDO
544
545!
546!-- Allocate building-data structure array. Note, this is a global array
547!-- and all building IDs on domain are known by each PE. Further attributes,
548!-- e.g. height-dependent arrays, however, are only allocated on PEs where
549!-- the respective building is present (in order to reduce memory demands).
550    ALLOCATE( buildings(1:num_build) )
551
552!
553!-- Store building IDs and check if building with certain ID is present on
554!-- subdomain.
555    DO  nb = 1, num_build
556       buildings(nb)%id = build_ids_final(nb)
557
558       IF ( ANY( building_id_f%var(nys:nyn,nxl:nxr) == buildings(nb)%id ) )    &
559          buildings(nb)%on_pe = .TRUE.
560    ENDDO
561!
562!-- Determine the maximum vertical dimension occupied by each building.
563    ALLOCATE( k_min_l(1:num_build) )
564    ALLOCATE( k_max_l(1:num_build) )
565    k_min_l = nzt + 1
566    k_max_l = 0   
567
568    DO  i = nxl, nxr
569       DO  j = nys, nyn
570          IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
571             nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ),    &
572                         DIM = 1 )
573             DO  k = nzb, nzt+1
574!
575!--             Check if grid point belongs to a building.
576                IF ( BTEST( wall_flags_0(k,j,i), 6 ) )  THEN
577                   k_min_l(nb) = MIN( k_min_l(nb), k )
578                   k_max_l(nb) = MAX( k_max_l(nb), k )
579                ENDIF
580
581             ENDDO
582          ENDIF
583       ENDDO
584    ENDDO
585
586    DO nb = 1, num_build
587#if defined( __parallel )
588       CALL MPI_ALLREDUCE( k_min_l(nb), buildings(nb)%kb_min, 1, MPI_INTEGER,  &
589                           MPI_MIN, comm2d, ierr )
590       CALL MPI_ALLREDUCE( k_max_l(nb), buildings(nb)%kb_max, 1, MPI_INTEGER,  &
591                           MPI_MAX, comm2d, ierr )
592#else
593       buildings(nb)%kb_min = k_min_l(nb)
594       buildings(nb)%kb_max = k_max_l(nb)
595#endif
596
597    ENDDO
598
599    DEALLOCATE( k_min_l )
600    DEALLOCATE( k_max_l )
601!
602!-- Calculate building height.
603    DO  nb = 1, num_build
604       buildings(nb)%building_height = 0.0_wp
605       DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
606          buildings(nb)%building_height = buildings(nb)%building_height        &
607                                        + dzw(k+1)
608       ENDDO
609    ENDDO
610!
611!-- Calculate building volume
612    DO  nb = 1, num_build
613!
614!--    Allocate temporary array for summing-up building volume
615       ALLOCATE( volume(buildings(nb)%kb_min:buildings(nb)%kb_max)   )
616       ALLOCATE( volume_l(buildings(nb)%kb_min:buildings(nb)%kb_max) )
617       volume   = 0.0_wp
618       volume_l = 0.0_wp
619!
620!--    Calculate building volume per height level on each PE where
621!--    these building is present.
622       IF ( buildings(nb)%on_pe )  THEN
623
624          ALLOCATE( buildings(nb)%volume(buildings(nb)%kb_min:buildings(nb)%kb_max)   )
625          ALLOCATE( buildings(nb)%vol_frac(buildings(nb)%kb_min:buildings(nb)%kb_max) )
626          buildings(nb)%volume   = 0.0_wp
627          buildings(nb)%vol_frac = 0.0_wp
628         
629          IF ( ANY( building_id_f%var(nys:nyn,nxl:nxr) == buildings(nb)%id ) ) &
630          THEN
631             DO  i = nxl, nxr
632                DO  j = nys, nyn
633                   DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
634                      IF ( building_id_f%var(j,i) /= building_id_f%fill )      &
635                         volume_l(k) = volume_l(k) + dx * dy * dzw(k+1)
636                   ENDDO
637                ENDDO
638             ENDDO
639          ENDIF
640       ENDIF
641!
642!--    Sum-up building volume from all subdomains
643#if defined( __parallel )
644       CALL MPI_ALLREDUCE( volume_l, volume, SIZE(volume), MPI_REAL, MPI_SUM,  &
645                           comm2d, ierr )
646#else
647       volume = volume_l
648#endif
649!
650!--    Save total building volume as well as local fraction on volume on
651!--    building data structure.
652       IF ( ALLOCATED( buildings(nb)%volume ) )  buildings(nb)%volume = volume
653!
654!--    Determine fraction of local on total building volume
655       IF ( buildings(nb)%on_pe )  buildings(nb)%vol_frac = volume_l / volume
656!
657!--    Calculate total building volume
658       IF ( ALLOCATED( buildings(nb)%volume ) )                                &
659          buildings(nb)%vol_tot = SUM( buildings(nb)%volume )
660
661       DEALLOCATE( volume   )
662       DEALLOCATE( volume_l )
663
664    ENDDO
665!
666!-- Allocate arrays for indoor temperature. 
667    DO  nb = 1, num_build
668       IF ( buildings(nb)%on_pe )  THEN
669          ALLOCATE( buildings(nb)%t_in(buildings(nb)%kb_min:buildings(nb)%kb_max)   )
670          ALLOCATE( buildings(nb)%t_in_l(buildings(nb)%kb_min:buildings(nb)%kb_max) )
671          buildings(nb)%t_in   = 0.0_wp
672          buildings(nb)%t_in_l = 0.0_wp
673       ENDIF
674    ENDDO
675!
676!-- Allocate arrays for number of facades per height level. Distinguish between
677!-- horizontal and vertical facades.
678    DO  nb = 1, num_build
679       IF ( buildings(nb)%on_pe )  THEN
680          ALLOCATE( buildings(nb)%num_facade_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
681          ALLOCATE( buildings(nb)%num_facade_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
682
683          buildings(nb)%num_facade_h = 0
684          buildings(nb)%num_facade_v = 0
685       ENDIF
686    ENDDO
687!
688!-- Determine number of facade elements per building on local subdomain.
689!-- Distinguish between horizontal and vertical facade elements.
690!
691!-- Horizontal facades
692    buildings(:)%num_facades_per_building_h_l = 0
693    DO  m = 1, surf_usm_h%ns
694!
695!--    For the current facade element determine corresponding building index.
696!--    First, obtain j,j,k indices of the building. Please note the
697!--    offset between facade/surface element and building location (for
698!--    horizontal surface elements the horizontal offsets are zero).
699       i = surf_usm_h%i(m) + surf_usm_h%ioff
700       j = surf_usm_h%j(m) + surf_usm_h%joff
701       k = surf_usm_h%k(m) + surf_usm_h%koff
702!
703!--    Determine building index and check whether building is on PE
704       nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
705
706       IF ( buildings(nb)%on_pe )  THEN
707!
708!--       Count number of facade elements at each height level.
709          buildings(nb)%num_facade_h(k) = buildings(nb)%num_facade_h(k) + 1 
710!
711!--       Moreover, sum up number of local facade elements per building.
712          buildings(nb)%num_facades_per_building_h_l =                         &
713                                buildings(nb)%num_facades_per_building_h_l + 1
714       ENDIF
715    ENDDO
716!
717!-- Vertical facades!
718    buildings(:)%num_facades_per_building_v_l = 0
719    DO  l = 0, 3
720       DO  m = 1, surf_usm_v(l)%ns
721!
722!--       For the current facade element determine corresponding building index.
723!--       First, obtain j,j,k indices of the building. Please note the
724!--       offset between facade/surface element and building location (for
725!--       vertical surface elements the vertical offsets are zero).
726          i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
727          j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
728          k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
729
730          nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ),        &
731                       DIM = 1 )
732          IF ( buildings(nb)%on_pe )  THEN
733             buildings(nb)%num_facade_v(k) = buildings(nb)%num_facade_v(k) + 1 
734             buildings(nb)%num_facades_per_building_v_l =                      &
735                                buildings(nb)%num_facades_per_building_v_l + 1
736          ENDIF
737       ENDDO
738    ENDDO
739
740!
741!-- Determine total number of facade elements per building and assign number to
742!-- building data type.
743    DO  nb = 1, num_build
744!
745!--    Allocate dummy array used for summing-up facade elements.
746!--    Please note, dummy arguments are necessary as building-date type
747!--    arrays are not necessarily allocated on all PEs.
748       ALLOCATE( num_facades_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
749       ALLOCATE( num_facades_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
750       ALLOCATE( receive_dum_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
751       ALLOCATE( receive_dum_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
752       num_facades_h = 0
753       num_facades_v = 0
754       receive_dum_h = 0
755       receive_dum_v = 0
756
757       IF ( buildings(nb)%on_pe )  THEN
758          num_facades_h = buildings(nb)%num_facade_h
759          num_facades_v = buildings(nb)%num_facade_v
760       ENDIF
761
762#if defined( __parallel )
763       CALL MPI_ALLREDUCE( num_facades_h,                                      &
764                           receive_dum_h,                                      &
765                           buildings(nb)%kb_max - buildings(nb)%kb_min + 1,    &
766                           MPI_INTEGER,                                        &
767                           MPI_SUM,                                            &
768                           comm2d,                                             &
769                           ierr )
770
771       CALL MPI_ALLREDUCE( num_facades_v,                                      &
772                           receive_dum_v,                                      &
773                           buildings(nb)%kb_max - buildings(nb)%kb_min + 1,    &
774                           MPI_INTEGER,                                        &
775                           MPI_SUM,                                            &
776                           comm2d,                                             &
777                           ierr )
778       IF ( ALLOCATED( buildings(nb)%num_facade_h ) )                          &
779           buildings(nb)%num_facade_h = receive_dum_h
780       IF ( ALLOCATED( buildings(nb)%num_facade_v ) )                          &
781           buildings(nb)%num_facade_v = receive_dum_v
782#else
783       buildings(nb)%num_facade_h = num_facades_h
784       buildings(nb)%num_facade_v = num_facades_v
785#endif
786
787!
788!--    Deallocate dummy arrays
789       DEALLOCATE( num_facades_h )
790       DEALLOCATE( num_facades_v )
791       DEALLOCATE( receive_dum_h )
792       DEALLOCATE( receive_dum_v )
793!
794!--    Allocate index arrays which link facade elements with surface-data type.
795!--    Please note, no height levels are considered here (information is stored
796!--    in surface-data type itself).
797       IF ( buildings(nb)%on_pe )  THEN
798!
799!--       Determine number of facade elements per building.
800          buildings(nb)%num_facades_per_building_h = SUM( buildings(nb)%num_facade_h )
801          buildings(nb)%num_facades_per_building_v = SUM( buildings(nb)%num_facade_v )
802!
803!--       Allocate arrays which link the building with the horizontal and vertical
804!--       urban-type surfaces. Please note, linking arrays are allocated over all
805!--       facade elements, which is required in case a building is located at the
806!--       subdomain boundaries, where the building and the corresponding surface
807!--       elements are located on different subdomains.
808          ALLOCATE( buildings(nb)%m_h(1:buildings(nb)%num_facades_per_building_h_l) )
809
810          ALLOCATE( buildings(nb)%l_v(1:buildings(nb)%num_facades_per_building_v_l) )
811          ALLOCATE( buildings(nb)%m_v(1:buildings(nb)%num_facades_per_building_v_l) )
812       ENDIF
813!
814! --    Determine volume per facade element (vpf)
815       IF ( buildings(nb)%on_pe )  THEN
816          ALLOCATE( buildings(nb)%vpf(buildings(nb)%kb_min:buildings(nb)%kb_max) )
817         
818          DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
819             buildings(nb)%vpf(k) = buildings(nb)%volume(k) /                  &
820                                REAL( buildings(nb)%num_facade_h(k) +          &
821                                      buildings(nb)%num_facade_v(k), KIND = wp )
822          ENDDO
823       ENDIF
824   
825!
826!--    Determine volume per total facade area (vpf). For the horizontal facade
827!--    area num_facades_per_building_h can be taken, multiplied with dx*dy.
828!--    However, due to grid stretching, vertical facade elements must be
829!--    summed-up vertically. Please note, if dx /= dy, an error is made!
830       IF ( buildings(nb)%on_pe )  THEN
831         
832          facade_area_v = 0.0_wp         
833          DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
834             facade_area_v = facade_area_v + buildings(nb)%num_facade_v(k)     &
835                             * dzw(k+1) * dx
836          ENDDO
837         
838          buildings(nb)%vpf = buildings(nb)%vol_tot /                          &
839                        ( buildings(nb)%num_facades_per_building_h * dx * dy + &
840                          facade_area_v )
841       ENDIF
842    ENDDO
843!
844!-- Link facade elements with surface data type.
845!-- Allocate array for counting.
846    ALLOCATE( n_fa(1:num_build) )
847    n_fa = 1
848
849    DO  m = 1, surf_usm_h%ns
850       i = surf_usm_h%i(m) + surf_usm_h%ioff
851       j = surf_usm_h%j(m) + surf_usm_h%joff
852
853       nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
854
855       IF ( buildings(nb)%on_pe )  THEN
856          buildings(nb)%m_h(n_fa(nb)) = m
857          n_fa(nb) = n_fa(nb) + 1 
858       ENDIF
859    ENDDO
860
861    n_fa = 1
862    DO  l = 0, 3
863       DO  m = 1, surf_usm_v(l)%ns
864          i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
865          j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
866
867          nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
868
869          IF ( buildings(nb)%on_pe )  THEN
870             buildings(nb)%l_v(n_fa(nb)) = l
871             buildings(nb)%m_v(n_fa(nb)) = m
872             n_fa(nb) = n_fa(nb) + 1   
873          ENDIF
874       ENDDO
875    ENDDO
876    DEALLOCATE( n_fa )
877
878!
879!-- Initialize building parameters, first by mean building type. Note,
880!-- in this case all buildings have the same type.
881!-- In a second step initialize with building tpyes from static input file,
882!-- where building types can be individual for each building.
883    buildings(:)%lambda_layer3       = building_pars(63,building_type)
884    buildings(:)%s_layer3            = building_pars(57,building_type)
885    buildings(:)%f_c_win             = building_pars(119,building_type)
886    buildings(:)%g_value_win         = building_pars(120,building_type)   
887    buildings(:)%u_value_win         = building_pars(121,building_type)   
888    buildings(:)%air_change_low      = building_pars(122,building_type)   
889    buildings(:)%air_change_high     = building_pars(123,building_type)   
890    buildings(:)%eta_ve              = building_pars(124,building_type)   
891    buildings(:)%factor_a            = building_pars(125,building_type)   
892    buildings(:)%factor_c            = building_pars(126,building_type)
893    buildings(:)%lambda_at           = building_pars(127,building_type)   
894    buildings(:)%theta_int_h_set     = building_pars(118,building_type)   
895    buildings(:)%theta_int_c_set     = building_pars(117,building_type)
896    buildings(:)%q_h_max             = building_pars(128,building_type)   
897    buildings(:)%q_c_max             = building_pars(129,building_type)         
898    buildings(:)%qint_high           = building_pars(130,building_type)
899    buildings(:)%qint_low            = building_pars(131,building_type)
900    buildings(:)%height_storey       = building_pars(132,building_type)
901    buildings(:)%height_cei_con      = building_pars(133,building_type)
902    buildings(:)%params_waste_heat_h = building_pars(134,building_type)
903    buildings(:)%params_waste_heat_c = building_pars(135,building_type)
904!
905!-- Initialize ventilaation load. Please note, building types > 7 are actually
906!-- not allowed (check already in urban_surface_mod and netcdf_data_input_mod.
907!-- However, the building data base may be later extended.
908    IF ( building_type ==  1  .OR.  building_type ==  2  .OR.                  &
909         building_type ==  3  .OR.  building_type == 10  .OR.                  &
910         building_type == 11  .OR.  building_type == 12 )  THEN
911       buildings(nb)%ventilation_int_loads = 1
912!
913!-- Office, building with large windows
914    ELSEIF ( building_type ==  4  .OR.  building_type ==  5  .OR.              &
915             building_type ==  6  .OR.  building_type ==  7  .OR.              &
916             building_type ==  8  .OR.  building_type ==  9)  THEN
917       buildings(nb)%ventilation_int_loads = 2
918!
919!-- Industry, hospitals
920    ELSEIF ( building_type == 13  .OR.  building_type == 14  .OR.              &
921             building_type == 15  .OR.  building_type == 16  .OR.              &
922             building_type == 17  .OR.  building_type == 18 )  THEN
923       buildings(nb)%ventilation_int_loads = 3
924    ENDIF
925!
926!-- Initialization of building parameters - level 2
927    IF ( building_type_f%from_file )  THEN
928       DO  i = nxl, nxr
929          DO  j = nys, nyn
930              IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
931                 nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
932                              DIM = 1 )
933                 bt = building_type_f%var(j,i)
934                 
935                 buildings(nb)%lambda_layer3       = building_pars(63,bt)
936                 buildings(nb)%s_layer3            = building_pars(57,bt)
937                 buildings(nb)%f_c_win             = building_pars(119,bt)
938                 buildings(nb)%g_value_win         = building_pars(120,bt)   
939                 buildings(nb)%u_value_win         = building_pars(121,bt)   
940                 buildings(nb)%air_change_low      = building_pars(122,bt)   
941                 buildings(nb)%air_change_high     = building_pars(123,bt)   
942                 buildings(nb)%eta_ve              = building_pars(124,bt)   
943                 buildings(nb)%factor_a            = building_pars(125,bt)   
944                 buildings(nb)%factor_c            = building_pars(126,bt)
945                 buildings(nb)%lambda_at           = building_pars(127,bt)   
946                 buildings(nb)%theta_int_h_set     = building_pars(118,bt)   
947                 buildings(nb)%theta_int_c_set     = building_pars(117,bt)
948                 buildings(nb)%q_h_max             = building_pars(128,bt)   
949                 buildings(nb)%q_c_max             = building_pars(129,bt)         
950                 buildings(nb)%qint_high           = building_pars(130,bt)
951                 buildings(nb)%qint_low            = building_pars(131,bt)
952                 buildings(nb)%height_storey       = building_pars(132,bt)
953                 buildings(nb)%height_cei_con      = building_pars(133,bt) 
954                 buildings(nb)%params_waste_heat_h = building_pars(134,bt)
955                 buildings(nb)%params_waste_heat_c = building_pars(135,bt)
956!
957!--              Initialize ventilaation load. Please note, building types > 7
958!--              are actually not allowed (check already in urban_surface_mod 
959!--              and netcdf_data_input_mod. However, the building data base may
960!--              be later extended.
961                 IF ( bt ==  1  .OR.  bt ==  2  .OR.                           &
962                      bt ==  3  .OR.  bt == 10  .OR.                           &
963                      bt == 11  .OR.  bt == 12 )  THEN
964                    buildings(nb)%ventilation_int_loads = 1
965!                   
966!--              Office, building with large windows
967                 ELSEIF ( bt ==  4  .OR.  bt ==  5  .OR.                       &
968                          bt ==  6  .OR.  bt ==  7  .OR.                       &
969                          bt ==  8  .OR.  bt ==  9)  THEN
970                    buildings(nb)%ventilation_int_loads = 2
971!
972!--              Industry, hospitals
973                 ELSEIF ( bt == 13  .OR.  bt == 14  .OR.                       &
974                          bt == 15  .OR.  bt == 16  .OR.                       &
975                          bt == 17  .OR.  bt == 18 )  THEN
976                    buildings(nb)%ventilation_int_loads = 3
977                 ENDIF
978              ENDIF
979           ENDDO
980        ENDDO
981    ENDIF
982
983!
984!-- Calculation of surface-related heat transfer coeffiecient
985!-- out of standard u-values from building database
986!-- only amount of extern and surface is used
987!-- otherwise amount between air and surface taken account twice
988    DO nb = 1, num_build
989       IF ( buildings(nb)%on_pe ) THEN   
990          du_win_tmp = 1.0_wp / buildings(nb)%u_value_win
991          u_tmp = buildings(nb)%u_value_win * ( du_win_tmp / ( du_win_tmp -    &
992                  0.125_wp + ( 1.0_wp / h_is ) ) )
993                 
994          du_tmp = 1.0_wp / u_tmp
995         
996          buildings(nb)%h_es = ( du_tmp / ( du_tmp - ( 1.0_wp / h_is ) ) ) *   &
997                                 u_tmp   
998       ENDIF
999    ENDDO
1000!
1001!-- Initial room temperature [K]
1002!-- (after first loop, use theta_m_t as theta_m_t_prev)
1003    theta_m_t_prev = initial_indoor_temperature
1004!
1005!-- Initialize indoor temperature. Actually only for output at initial state.
1006    DO  nb = 1, num_build
1007       buildings(nb)%t_in(:) = initial_indoor_temperature
1008    ENDDO
1009
1010    CALL location_message( 'finished', .TRUE. )
1011
1012 END SUBROUTINE im_init
1013
1014
1015!------------------------------------------------------------------------------!
1016! Description:
1017! ------------
1018!> Main part of the indoor model.
1019!> Calculation of .... (kanani: Please describe)
1020!------------------------------------------------------------------------------!
1021 SUBROUTINE im_main_heatcool
1022
1023    USE arrays_3d,                                                             &
1024        ONLY:  ddzw, dzw
1025
1026!     USE basic_constants_and_equations_mod,                                     &
1027!         ONLY:  c_p
1028
1029!     USE control_parameters,                                                    &
1030!         ONLY:  rho_surface
1031
1032    USE date_and_time_mod,                                                     &
1033        ONLY:  time_utc
1034
1035    USE grid_variables,                                                        &
1036        ONLY:  dx, dy
1037
1038    USE pegrid
1039   
1040    USE surface_mod,                                                           &
1041        ONLY:  ind_veg_wall, ind_wat_win, surf_usm_h, surf_usm_v
1042
1043    USE urban_surface_mod,                                                     &
1044        ONLY:  nzt_wall, t_wall_h, t_wall_v, t_window_h, t_window_v,           &
1045               building_type
1046
1047
1048    IMPLICIT NONE
1049
1050    INTEGER(iwp) ::  i    !< index of facade-adjacent atmosphere grid point in x-direction
1051    INTEGER(iwp) ::  j    !< index of facade-adjacent atmosphere grid point in y-direction
1052    INTEGER(iwp) ::  k    !< index of facade-adjacent atmosphere grid point in z-direction
1053    INTEGER(iwp) ::  kk   !< vertical index of indoor grid point adjacent to facade
1054    INTEGER(iwp) ::  l    !< running index for surface-element orientation
1055    INTEGER(iwp) ::  m    !< running index surface elements
1056    INTEGER(iwp) ::  nb   !< running index for buildings
1057    INTEGER(iwp) ::  fa   !< running index for facade elements of each building
1058
1059    REAL(wp) ::  indoor_wall_window_temperature   !< weighted temperature of innermost wall/window layer
1060    REAL(wp) ::  near_facade_temperature          !< outside air temperature 10cm away from facade
1061    REAL(wp) ::  time_utc_hour                    !< time of day (hour UTC)
1062
1063    REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in_l_send   !< dummy send buffer used for summing-up indoor temperature per kk-level
1064    REAL(wp), DIMENSION(:), ALLOCATABLE ::  t_in_recv     !< dummy recv buffer used for summing-up indoor temperature per kk-level
1065!
1066!-- Determine time of day in hours.
1067    time_utc_hour = time_utc / 3600.0_wp
1068!
1069!-- Following calculations must be done for each facade element.
1070    DO  nb = 1, num_build
1071!
1072!--    First, check whether building is present on local subdomain.
1073       IF ( buildings(nb)%on_pe )  THEN
1074!
1075!--       Determine daily schedule. 08:00-18:00 = 1, other hours = 0.
1076!--       Residental Building, panel WBS 70   
1077          IF ( buildings(nb)%ventilation_int_loads == 1 )  THEN
1078             IF ( time_utc_hour >= 6.0_wp  .AND.  time_utc_hour <= 8.0_wp )  THEN
1079                schedule_d = 1
1080             ELSEIF ( time_utc_hour >= 18.0_wp  .AND.  time_utc_hour <= 23.0_wp )  THEN
1081                schedule_d = 1
1082             ELSE
1083                schedule_d = 0
1084             ENDIF
1085          ENDIF
1086!
1087!--       Office, building with large windows
1088          IF ( buildings(nb)%ventilation_int_loads == 2 )  THEN
1089             IF ( time_utc_hour >= 8.0_wp  .AND.  time_utc_hour <= 18.0_wp )  THEN
1090                schedule_d = 1
1091             ELSE
1092                schedule_d = 0
1093             ENDIF
1094          ENDIF
1095!       
1096!--       Industry, hospitals
1097          IF ( buildings(nb)%ventilation_int_loads == 3 )  THEN
1098             IF ( time_utc_hour >= 6.0_wp  .AND.  time_utc_hour <= 22.0_wp )  THEN
1099                schedule_d = 1
1100             ELSE
1101                schedule_d = 0
1102             ENDIF
1103          ENDIF
1104!
1105!--       Initialize/reset indoor temperature
1106          buildings(nb)%t_in_l = 0.0_wp 
1107!
1108!--       Horizontal surfaces
1109          DO  fa = 1, buildings(nb)%num_facades_per_building_h_l
1110!
1111!--          Determine index where corresponding surface-type information
1112!--          is stored.
1113             m = buildings(nb)%m_h(fa)
1114!
1115!--          Determine building height level index.
1116             kk = surf_usm_h%k(m) + surf_usm_h%koff
1117!           
1118!--          Building geometries --> not time-dependent
1119             facade_element_area          = dx * dy                                                   !< [m2] surface area per facade element   
1120             floor_area_per_facade        = buildings(nb)%vpf(kk) * ddzw(kk+1)                        !< [m2/m2] floor area per facade area
1121             indoor_volume_per_facade     = buildings(nb)%vpf(kk)                                     !< [m3/m2] indoor air volume per facade area
1122             buildings(nb)%area_facade    = facade_element_area *                                   &
1123                                          ( buildings(nb)%num_facades_per_building_h +              &
1124                                            buildings(nb)%num_facades_per_building_v )                !< [m2] area of total facade
1125             window_area_per_facade       = surf_usm_h%frac(ind_wat_win,m)  * facade_element_area     !< [m2] window area per facade element
1126
1127             buildings(nb)%net_floor_area = buildings(nb)%vol_tot / ( buildings(nb)%height_storey )
1128             total_area                   = buildings(nb)%net_floor_area                              !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
1129             a_m                          = buildings(nb)%factor_a * total_area *                   &
1130                                          ( facade_element_area / buildings(nb)%area_facade ) *     &
1131                                            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)
1132             c_m                          = buildings(nb)%factor_c * total_area *                   &
1133                                          ( 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)             
1134!
1135!--          Calculation of heat transfer coefficient for transmission --> not time-dependent
1136             h_t_es   = window_area_per_facade * buildings(nb)%h_es                                   !< [W/K] only for windows
1137
1138             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)
1139             h_t_ms  = a_m * h_ms                                                                     !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
1140             h_t_wall  = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade )     & !< [W/K]
1141                                    * buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp &
1142                                             ) + 1.0_wp / h_t_ms )                                    !< [W/K] opaque components
1143             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
1144!
1145!--          internal air loads dependent on the occupacy of the room
1146!--          basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
1147             phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )  &
1148                              * floor_area_per_facade ) 
1149             q_int = phi_ia / total_area
1150!
1151!--          Airflow dependent on the occupacy of the room
1152!--          basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
1153             air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low )  !< [1/h]?
1154!
1155!--          Heat transfer of ventilation
1156!--          not less than 0.01 W/K to provide division by 0 in further calculations
1157!--          with heat capacity of air 0.33 Wh/m2K
1158             h_v   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *      &
1159                                    0.33_wp * (1.0_wp - buildings(nb)%eta_ve ) ) )    !< [W/K] from ISO 13789 Eq.(10)
1160
1161!--          Heat transfer coefficient auxiliary variables
1162             h_t_1 = 1.0_wp / ( ( 1.0_wp / h_v )   + ( 1.0_wp / h_t_is ) )  !< [W/K] Eq. (C.6)
1163             h_t_2 = h_t_1 + h_t_es                                         !< [W/K] Eq. (C.7)
1164             h_t_3 = 1.0_wp / ( ( 1.0_wp / h_t_2 ) + ( 1.0_wp / h_t_ms ) )  !< [W/K] Eq. (C.8)
1165!
1166!--          Net short-wave radiation through window area (was i_global)
1167             net_sw_in = surf_usm_h%rad_sw_in(m) - surf_usm_h%rad_sw_out(m)
1168!
1169!--          Quantities needed for im_calc_temperatures
1170             i = surf_usm_h%i(m)
1171             j = surf_usm_h%j(m)
1172             k = surf_usm_h%k(m)
1173             near_facade_temperature = surf_usm_h%pt_10cm(m)
1174             indoor_wall_window_temperature =                                  &
1175                  surf_usm_h%frac(ind_veg_wall,m) * t_wall_h(nzt_wall,m)       &
1176                + surf_usm_h%frac(ind_wat_win,m)  * t_window_h(nzt_wall,m)
1177!
1178!--          Solar thermal gains. If net_sw_in larger than sun-protection
1179!--          threshold parameter (params_solar_protection), sun protection will
1180!--          be activated
1181             IF ( net_sw_in <= params_solar_protection )  THEN
1182                solar_protection_off = 1
1183                solar_protection_on  = 0
1184             ELSE
1185                solar_protection_off = 0
1186                solar_protection_on  = 1
1187             ENDIF
1188!
1189!--          Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
1190!--          DIN 4108 - 2 chap.8
1191             phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off                           &
1192                         + window_area_per_facade * net_sw_in * buildings(nb)%f_c_win * solar_protection_on )  &
1193                       * buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w
1194             q_sol = phi_sol           
1195!
1196!--          Calculation of the mass specific thermal load for internal and external heatsources of the inner node
1197             phi_m   = (a_m / total_area) * ( phi_ia + phi_sol )                                    !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
1198             q_c_m = phi_m
1199!
1200!--          Calculation mass specific thermal load implied non thermal mass
1201             phi_st  =   ( 1.0_wp - ( a_m / total_area ) - ( h_t_es / ( 9.1_wp * total_area ) ) ) &
1202                       * ( phi_ia + phi_sol )                                                       !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
1203             q_c_st = phi_st           
1204!
1205!--          Calculations for deriving indoor temperature and heat flux into the wall
1206!--          Step 1: Indoor temperature without heating and cooling
1207!--          section C.4.1 Picture C.2 zone 3)
1208             phi_hc_nd = 0.0_wp
1209             
1210             CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1211                                         near_facade_temperature, phi_hc_nd )
1212!
1213!--          If air temperature between border temperatures of heating and cooling, assign output variable, then ready   
1214             IF ( buildings(nb)%theta_int_h_set <= theta_air  .AND.  theta_air <= buildings(nb)%theta_int_c_set )  THEN
1215                phi_hc_nd_ac = 0.0_wp
1216                phi_hc_nd    = phi_hc_nd_ac           
1217                theta_air_ac = theta_air
1218!
1219!--          Step 2: Else, apply 10 W/m2 heating/cooling power and calculate indoor temperature
1220!--          again.
1221             ELSE
1222!
1223!--             Temperature not correct, calculation method according to section C4.2
1224                theta_air_0 = theta_air                                                  !< temperature without heating/cooling 
1225!
1226!--             Heating or cooling?
1227                IF ( theta_air_0 > buildings(nb)%theta_int_c_set )  THEN
1228                   theta_air_set = buildings(nb)%theta_int_c_set
1229                ELSE
1230                   theta_air_set = buildings(nb)%theta_int_h_set
1231                ENDIF
1232!
1233!--             Calculate the temperature with phi_hc_nd_10
1234                phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
1235                phi_hc_nd    = phi_hc_nd_10
1236               
1237                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1238                                            near_facade_temperature, phi_hc_nd )
1239                theta_air_10 = theta_air                                                !< temperature with 10 W/m2 of heating
1240                phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0)          &
1241                                            / (theta_air_10  - theta_air_0)             !< Eq. (C.13)
1242!
1243!--             Step 3: With temperature ratio to determine the heating or cooling capacity   
1244!--             If necessary, limit the power to maximum power
1245!--             section C.4.1 Picture C.2 zone 2) and 4)
1246                buildings(nb)%phi_c_max = buildings(nb)%q_c_max * floor_area_per_facade             
1247                buildings(nb)%phi_h_max = buildings(nb)%q_h_max * floor_area_per_facade
1248                IF ( buildings(nb)%phi_c_max < phi_hc_nd_un  .AND.  phi_hc_nd_un < buildings(nb)%phi_h_max )  THEN
1249                   phi_hc_nd_ac = phi_hc_nd_un
1250                   phi_hc_nd = phi_hc_nd_un 
1251                ELSE
1252!
1253!--             Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
1254!--             section C.4.1 Picture C.2 zone 1) and 5)
1255                   IF ( phi_hc_nd_un > 0.0_wp )  THEN
1256                      phi_hc_nd_ac = buildings(nb)%phi_h_max                            !< Limit heating
1257                   ELSE
1258                      phi_hc_nd_ac = buildings(nb)%phi_c_max                            !< Limit cooling
1259                   ENDIF
1260                ENDIF
1261                phi_hc_nd = phi_hc_nd_ac   
1262!
1263!--             Calculate the temperature with phi_hc_nd_ac (new)
1264                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1265                                            near_facade_temperature, phi_hc_nd )
1266                theta_air_ac = theta_air
1267             ENDIF
1268!
1269!--          Update theta_m_t_prev
1270             theta_m_t_prev = theta_m_t
1271             
1272             q_vent = h_v * ( theta_air - near_facade_temperature )
1273!
1274!--          Calculate the operating temperature with weighted mean temperature of air and mean solar temperature
1275!--          Will be used for thermal comfort calculations
1276             theta_op     = 0.3_wp * theta_air_ac + 0.7_wp * theta_s          !< [degree_C] operative Temperature Eq. (C.12)
1277!              surf_usm_h%t_indoor(m) = theta_op                               !< not integrated now
1278!
1279!--          Heat flux into the wall. Value needed in urban_surface_mod to
1280!--          calculate heat transfer through wall layers towards the facade
1281!--          (use c_p * rho_surface to convert [W/m2] into [K m/s])
1282             q_wall_win = h_t_ms * ( theta_s - theta_m )                       &
1283                                    / (   facade_element_area                  &
1284                                        - window_area_per_facade )
1285             q_trans = q_wall_win * facade_element_area                                       
1286!
1287!--          Transfer q_wall_win back to USM (innermost wall/window layer)
1288             surf_usm_h%iwghf_eb(m)        = q_wall_win
1289             surf_usm_h%iwghf_eb_window(m) = q_wall_win
1290!
1291!--          Sum up operational indoor temperature per kk-level. Further below,
1292!--          this temperature is reduced by MPI to one temperature per kk-level
1293!--          and building (processor overlapping)
1294             buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
1295!
1296!--          Calculation of waste heat
1297!--          Anthropogenic heat output
1298             IF ( phi_hc_nd_ac > 0.0_wp )  THEN
1299                heating_on = 1
1300                cooling_on = 0
1301             ELSE
1302                heating_on = 0
1303                cooling_on = -1
1304             ENDIF
1305
1306             q_waste_heat = ( phi_hc_nd * (                                    &
1307                              buildings(nb)%params_waste_heat_h * heating_on + &
1308                              buildings(nb)%params_waste_heat_c * cooling_on ) &
1309                            ) / facade_element_area                                               !< [W/m2] , observe the directional convention in PALM!
1310             surf_usm_h%waste_heat(m) = q_waste_heat
1311          ENDDO !< Horizontal surfaces loop
1312!
1313!--       Vertical surfaces
1314          DO  fa = 1, buildings(nb)%num_facades_per_building_v_l
1315!
1316!--          Determine indices where corresponding surface-type information
1317!--          is stored.
1318             l = buildings(nb)%l_v(fa)
1319             m = buildings(nb)%m_v(fa)
1320!
1321!--          Determine building height level index.
1322             kk = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
1323!
1324!--          (SOME OF THE FOLLOWING (not time-dependent COULD PROBABLY GO INTO A FUNCTION
1325!--          EXCEPT facade_element_area, EVERYTHING IS CALCULATED EQUALLY)
1326!--          Building geometries  --> not time-dependent
1327             IF ( l == 0  .OR. l == 1 ) facade_element_area = dx * dzw(kk+1)                          !< [m2] surface area per facade element
1328             IF ( l == 2  .OR. l == 3 ) facade_element_area = dy * dzw(kk+1)                          !< [m2] surface area per facade element
1329             floor_area_per_facade        = buildings(nb)%vpf(kk) * ddzw(kk+1)                        !< [m2/m2] floor area per facade area
1330             indoor_volume_per_facade     = buildings(nb)%vpf(kk)                                     !< [m3/m2] indoor air volume per facade area
1331             buildings(nb)%area_facade    = facade_element_area *                                   &
1332                                          ( buildings(nb)%num_facades_per_building_h +              &
1333                                            buildings(nb)%num_facades_per_building_v )                !< [m2] area of total facade
1334             window_area_per_facade       = surf_usm_v(l)%frac(ind_wat_win,m)  * facade_element_area     !< [m2] window area per facade element
1335
1336             buildings(nb)%net_floor_area = buildings(nb)%vol_tot / ( buildings(nb)%height_storey )
1337             total_area                   = buildings(nb)%net_floor_area                              !< [m2] area of all surfaces pointing to zone  Eq. (9) according to section 7.2.2.2
1338             a_m                          = buildings(nb)%factor_a * total_area *                   &
1339                                          ( facade_element_area / buildings(nb)%area_facade ) *     &
1340                                            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)
1341             c_m                          = buildings(nb)%factor_c * total_area *                   &
1342                                          ( 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)
1343!
1344!--          Calculation of heat transfer coefficient for transmission --> not time-dependent
1345             h_t_es   = window_area_per_facade * buildings(nb)%h_es                                   !< [W/K] only for windows
1346
1347             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)
1348             h_t_ms  = a_m * h_ms                                                                     !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
1349             h_t_wall  = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade )     & !< [W/K]
1350                                    * buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp &
1351                                             ) + 1.0_wp / h_t_ms )                                    !< [W/K] opaque components
1352             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
1353!
1354!--          internal air loads dependent on the occupacy of the room
1355!--          basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
1356             phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low )  &
1357                              * floor_area_per_facade )
1358             q_int = phi_ia
1359
1360!
1361!--          Airflow dependent on the occupacy of the room
1362!--          basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
1363             air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low ) 
1364!
1365!--          Heat transfer of ventilation
1366!--          not less than 0.01 W/K to provide division by 0 in further calculations
1367!--          with heat capacity of air 0.33 Wh/m2K
1368             h_v   = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade *                       &
1369                                    0.33_wp * (1.0_wp - buildings(nb)%eta_ve ) ) )                    !< [W/K] from ISO 13789 Eq.(10)
1370                                   
1371!--          Heat transfer coefficient auxiliary variables
1372             h_t_1 = 1.0_wp / ( ( 1.0_wp / h_v )   + ( 1.0_wp / h_t_is ) )                            !< [W/K] Eq. (C.6)
1373             h_t_2 = h_t_1 + h_t_es                                                                   !< [W/K] Eq. (C.7)
1374             h_t_3 = 1.0_wp / ( ( 1.0_wp / h_t_2 ) + ( 1.0_wp / h_t_ms ) )                            !< [W/K] Eq. (C.8)
1375!
1376!--          Net short-wave radiation through window area (was i_global)
1377             net_sw_in = surf_usm_v(l)%rad_sw_in(m) - surf_usm_v(l)%rad_sw_out(m)
1378!
1379!--          Quantities needed for im_calc_temperatures
1380             i = surf_usm_v(l)%i(m)
1381             j = surf_usm_v(l)%j(m)   
1382             k = surf_usm_v(l)%k(m)
1383             near_facade_temperature = surf_usm_v(l)%pt_10cm(m)
1384             indoor_wall_window_temperature =                                                       &
1385                  surf_usm_v(l)%frac(ind_veg_wall,m) * t_wall_v(l)%t(nzt_wall,m)                    &
1386                + surf_usm_v(l)%frac(ind_wat_win,m)  * t_window_v(l)%t(nzt_wall,m)
1387!
1388!--          Solar thermal gains. If net_sw_in larger than sun-protection
1389!--          threshold parameter (params_solar_protection), sun protection will
1390!--          be activated
1391             IF ( net_sw_in <= params_solar_protection )  THEN
1392                solar_protection_off = 1
1393                solar_protection_on  = 0 
1394             ELSE
1395                solar_protection_off = 0
1396                solar_protection_on  = 1 
1397             ENDIF
1398!
1399!--          Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
1400!--          DIN 4108 - 2 chap.8
1401             phi_sol = (   window_area_per_facade * net_sw_in * solar_protection_off                             &
1402                         + window_area_per_facade * net_sw_in * buildings(nb)%f_c_win * solar_protection_on )    &
1403                       * buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w
1404             q_sol = phi_sol
1405!
1406!--          Calculation of the mass specific thermal load for internal and external heatsources
1407             phi_m   = (a_m / total_area) * ( phi_ia + phi_sol )          !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
1408             q_c_m = phi_m
1409!
1410!--          Calculation mass specific thermal load implied non thermal mass
1411             phi_st  =   ( 1.0_wp - ( a_m / total_area ) - ( h_t_es / ( 9.1_wp * total_area ) ) )                &
1412                       * ( phi_ia + phi_sol )                                                                       !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
1413             q_c_st = phi_st
1414!
1415!--          Calculations for deriving indoor temperature and heat flux into the wall
1416!--          Step 1: Indoor temperature without heating and cooling
1417!--          section C.4.1 Picture C.2 zone 3)
1418             phi_hc_nd = 0.0_wp
1419             CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1420                                         near_facade_temperature, phi_hc_nd )
1421!
1422!--          If air temperature between border temperatures of heating and cooling, assign output variable, then ready 
1423             IF ( buildings(nb)%theta_int_h_set <= theta_air  .AND.  theta_air <= buildings(nb)%theta_int_c_set )  THEN
1424                phi_hc_nd_ac = 0.0_wp
1425                phi_hc_nd    = phi_hc_nd_ac
1426                theta_air_ac = theta_air
1427!
1428!--          Step 2: Else, apply 10 W/m2 heating/cooling power and calculate indoor temperature
1429!--          again.
1430             ELSE
1431!
1432!--             Temperature not correct, calculation method according to section C4.2
1433                theta_air_0 = theta_air !< Note temperature without heating/cooling
1434!
1435!--             Heating or cooling?
1436                IF ( theta_air_0 > buildings(nb)%theta_int_c_set )  THEN
1437                   theta_air_set = buildings(nb)%theta_int_c_set
1438                ELSE
1439                   theta_air_set = buildings(nb)%theta_int_h_set
1440                ENDIF
1441
1442!--             Calculate the temperature with phi_hc_nd_10
1443                phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
1444                phi_hc_nd    = phi_hc_nd_10
1445       
1446                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1447                                            near_facade_temperature, phi_hc_nd )
1448
1449                theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
1450
1451               
1452                phi_hc_nd_un = phi_hc_nd_10 * ( theta_air_set - theta_air_0 )  &
1453                                            / ( theta_air_10  - theta_air_0 )            !< Eq. (C.13)
1454!
1455!--             Step 3: With temperature ratio to determine the heating or cooling capacity   
1456!--             If necessary, limit the power to maximum power
1457!--             section C.4.1 Picture C.2 zone 2) and 4)
1458                buildings(nb)%phi_c_max = buildings(nb)%q_c_max * floor_area_per_facade
1459                buildings(nb)%phi_h_max = buildings(nb)%q_h_max * floor_area_per_facade
1460                IF ( buildings(nb)%phi_c_max < phi_hc_nd_un  .AND.  phi_hc_nd_un < buildings(nb)%phi_h_max )  THEN
1461                   phi_hc_nd_ac = phi_hc_nd_un
1462                   phi_hc_nd = phi_hc_nd_un
1463                ELSE
1464!
1465!--             Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
1466!--             section C.4.1 Picture C.2 zone 1) and 5)
1467                   IF ( phi_hc_nd_un > 0.0_wp )  THEN
1468                      phi_hc_nd_ac = buildings(nb)%phi_h_max                                         !< Limit heating
1469                   ELSE
1470                      phi_hc_nd_ac = buildings(nb)%phi_c_max                                         !< Limit cooling
1471                   ENDIF
1472                ENDIF
1473                phi_hc_nd = phi_hc_nd_ac 
1474!
1475!--             Calculate the temperature with phi_hc_nd_ac (new)
1476                CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1477                                            near_facade_temperature, phi_hc_nd )
1478                theta_air_ac = theta_air
1479             ENDIF
1480!
1481!--          Update theta_m_t_prev
1482             theta_m_t_prev = theta_m_t
1483             
1484             q_vent = h_v * ( theta_air - near_facade_temperature )
1485!
1486!--          Calculate the operating temperature with weighted mean of temperature of air and mean
1487!--          Will be used for thermal comfort calculations
1488             theta_op     = 0.3_wp * theta_air_ac + 0.7_wp * theta_s
1489!              surf_usm_v(l)%t_indoor(m) = theta_op                  !< not integrated yet
1490!
1491!--          Heat flux into the wall. Value needed in urban_surface_mod to
1492!--          calculate heat transfer through wall layers towards the facade
1493             q_wall_win = h_t_ms * ( theta_s - theta_m )                       &
1494                                    / (   facade_element_area                  &
1495                                        - window_area_per_facade )
1496             q_trans = q_wall_win * facade_element_area
1497!
1498!--          Transfer q_wall_win back to USM (innermost wall/window layer)
1499             surf_usm_v(l)%iwghf_eb(m)        = q_wall_win
1500             surf_usm_v(l)%iwghf_eb_window(m) = q_wall_win
1501!
1502!--          Sum up operational indoor temperature per kk-level. Further below,
1503!--          this temperature is reduced by MPI to one temperature per kk-level
1504!--          and building (processor overlapping)
1505             buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
1506!
1507!--          Calculation of waste heat
1508!--          Anthropogenic heat output
1509             IF ( phi_hc_nd_ac > 0.0_wp )  THEN
1510                heating_on = 1
1511                cooling_on = 0
1512             ELSE
1513                heating_on = 0
1514                cooling_on = -1
1515             ENDIF
1516
1517             q_waste_heat = ( phi_hc_nd * (                                    &
1518                    buildings(nb)%params_waste_heat_h * heating_on +           &
1519                    buildings(nb)%params_waste_heat_c * cooling_on )           &
1520                            ) / facade_element_area !< [W/m2] , observe the directional convention in PALM!
1521             surf_usm_v(l)%waste_heat(m) = q_waste_heat
1522             
1523          ENDDO !< Vertical surfaces loop
1524
1525       ENDIF !< buildings(nb)%on_pe
1526    ENDDO !< buildings loop
1527
1528!
1529!-- Determine the mean building temperature.
1530    DO  nb = 1, num_build
1531!
1532!--    Allocate dummy array used for summing-up facade elements.
1533!--    Please note, dummy arguments are necessary as building-date type
1534!--    arrays are not necessarily allocated on all PEs.
1535       ALLOCATE( t_in_l_send(buildings(nb)%kb_min:buildings(nb)%kb_max) )
1536       ALLOCATE( t_in_recv(buildings(nb)%kb_min:buildings(nb)%kb_max) )
1537       t_in_l_send = 0.0_wp
1538       t_in_recv   = 0.0_wp
1539
1540       IF ( buildings(nb)%on_pe )  THEN
1541          t_in_l_send = buildings(nb)%t_in_l
1542       ENDIF
1543
1544
1545#if defined( __parallel )
1546       CALL MPI_ALLREDUCE( t_in_l_send,                                        &
1547                           t_in_recv,                                          &
1548                           buildings(nb)%kb_max - buildings(nb)%kb_min + 1,    &
1549                           MPI_REAL,                                           &
1550                           MPI_SUM,                                            &
1551                           comm2d,                                             &
1552                           ierr )
1553
1554       IF ( ALLOCATED( buildings(nb)%t_in ) )                                  &
1555           buildings(nb)%t_in = t_in_recv
1556#else
1557       IF ( ALLOCATED( buildings(nb)%t_in ) )                                  &
1558          buildings(nb)%t_in = buildings(nb)%t_in_l
1559#endif
1560
1561       IF ( ALLOCATED( buildings(nb)%t_in ) )  THEN
1562!
1563!--       Average indoor temperature. Note, in case a building is completely
1564!--       surrounded by higher buildings, it may have no facade elements
1565!--       at some height levels, will will lead to a divide by zero. If this
1566!--       is the case, indoor temperature will be set to -1.0.
1567          DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
1568             IF ( buildings(nb)%num_facade_h(k) +                              &
1569                  buildings(nb)%num_facade_v(k) > 0 )  THEN
1570                buildings(nb)%t_in(k) = buildings(nb)%t_in(k) /                &
1571                               REAL( buildings(nb)%num_facade_h(k) +           &
1572                                     buildings(nb)%num_facade_v(k), KIND = wp )
1573             ELSE
1574                buildings(nb)%t_in(k) = -1.0_wp
1575             ENDIF
1576          ENDDO
1577       ENDIF
1578       
1579
1580!
1581!--    Deallocate dummy arrays
1582       DEALLOCATE( t_in_l_send )
1583       DEALLOCATE( t_in_recv )
1584
1585    ENDDO
1586   
1587 END SUBROUTINE im_main_heatcool
1588
1589!-----------------------------------------------------------------------------!
1590! Description:
1591!-------------
1592!> Check data output for plant canopy model
1593!-----------------------------------------------------------------------------!
1594 SUBROUTINE im_check_data_output( var, unit )
1595       
1596    IMPLICIT NONE
1597   
1598    CHARACTER (LEN=*) ::  unit   !<
1599    CHARACTER (LEN=*) ::  var    !<
1600       
1601    SELECT CASE ( TRIM( var ) )
1602   
1603   
1604        CASE ( 'im_hf_roof')
1605           unit = 'W m-2'
1606       
1607        CASE ( 'im_hf_wall_win' )
1608           unit = 'W m-2'
1609           
1610        CASE ( 'im_hf_wall_win_waste' )
1611           unit = 'W m-2'
1612           
1613        CASE ( 'im_hf_roof_waste' )
1614           unit = 'W m-2'
1615       
1616        CASE ( 'im_t_indoor_mean' )
1617           unit = 'K'
1618           
1619        CASE ( 'im_t_indoor_roof' )
1620           unit = 'K'
1621           
1622        CASE ( 'im_t_indoor_wall_win' )
1623           unit = 'K'
1624       
1625        CASE DEFAULT
1626           unit = 'illegal'
1627           
1628    END SELECT
1629   
1630 END SUBROUTINE
1631
1632
1633!-----------------------------------------------------------------------------!
1634! Description:
1635!-------------
1636!> Check parameters routine for plant canopy model
1637!-----------------------------------------------------------------------------!
1638 SUBROUTINE im_check_parameters
1639
1640!   USE control_parameters,
1641!       ONLY: message_string
1642       
1643   IMPLICIT NONE
1644   
1645 END SUBROUTINE im_check_parameters
1646
1647!-----------------------------------------------------------------------------!
1648! Description:
1649!-------------
1650!> Subroutine defining appropriate grid for netcdf variables.
1651!> It is called from subroutine netcdf.
1652!-----------------------------------------------------------------------------!
1653 SUBROUTINE im_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1654
1655   IMPLICIT NONE
1656   
1657   CHARACTER (LEN=*), INTENT(IN)  ::  var
1658   LOGICAL, INTENT(OUT)           ::  found
1659   CHARACTER (LEN=*), INTENT(OUT) ::  grid_x
1660   CHARACTER (LEN=*), INTENT(OUT) ::  grid_y
1661   CHARACTER (LEN=*), INTENT(OUT) ::  grid_z
1662   
1663   found   = .TRUE.
1664   
1665!
1666!-- Check for the grid
1667    SELECT CASE ( TRIM( var ) )
1668
1669       CASE ( 'im_hf_roof', 'im_hf_roof_waste' )
1670          grid_x = 'x'
1671          grid_y = 'y'
1672          grid_z = 'zw'
1673!
1674!--    Heat fluxes at vertical walls are actually defined on stagged grid, i.e. xu, yv.
1675       CASE ( 'im_hf_wall_win', 'im_hf_wall_win_waste' )
1676          grid_x = 'x'
1677          grid_y = 'y'
1678          grid_z = 'zu'
1679
1680       CASE ( 'im_t_indoor_mean', 'im_t_indoor_roof', 'im_t_indoor_wall_win')
1681          grid_x = 'x'
1682          grid_y = 'y'
1683          grid_z = 'zw'
1684         
1685       CASE DEFAULT
1686          found  = .FALSE.
1687          grid_x = 'none'
1688          grid_y = 'none'
1689          grid_z = 'none'
1690    END SELECT
1691   
1692 END SUBROUTINE im_define_netcdf_grid
1693
1694!------------------------------------------------------------------------------!
1695! Description:
1696! ------------
1697!> Subroutine defining 3D output variables
1698!------------------------------------------------------------------------------!
1699 SUBROUTINE im_data_output_3d( av, variable, found, local_pf, fill_value,      &
1700                               nzb_do, nzt_do )
1701
1702   USE indices
1703   
1704   USE kinds
1705         
1706   IMPLICIT NONE
1707   
1708    CHARACTER (LEN=*) ::  variable !<
1709
1710    INTEGER(iwp) ::  av    !<
1711    INTEGER(iwp) ::  i     !<
1712    INTEGER(iwp) ::  j     !<
1713    INTEGER(iwp) ::  k     !<
1714    INTEGER(iwp) ::  l     !<
1715    INTEGER(iwp) ::  m     !<
1716    INTEGER(iwp) ::  nb    !< index of the building in the building data structure
1717    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
1718    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1719   
1720    LOGICAL      ::  found !<
1721
1722    REAL(wp), INTENT(IN) ::  fill_value !< value for the _FillValue attribute
1723
1724    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
1725   
1726    local_pf = fill_value
1727   
1728    found = .TRUE.
1729   
1730    SELECT CASE ( TRIM( variable ) )
1731!
1732!--     Output of indoor temperature. All grid points within the building are
1733!--     filled with values, while atmospheric grid points are set to _FillValues.
1734        CASE ( 'im_t_indoor_mean' )
1735           IF ( av == 0 ) THEN
1736              DO  i = nxl, nxr
1737                 DO  j = nys, nyn
1738                    IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
1739!
1740!--                    Determine index of the building within the building data structure.
1741                       nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ),   &
1742                                    DIM = 1 )
1743                       IF ( buildings(nb)%on_pe )  THEN
1744!
1745!--                       Write mean building temperature onto output array. Please note,
1746!--                       in contrast to many other loops in the output, the vertical
1747!--                       bounds are determined by the lowest and hightest vertical index
1748!--                       occupied by the building.
1749                          DO  k = buildings(nb)%kb_min, buildings(nb)%kb_max
1750                             local_pf(i,j,k) = buildings(nb)%t_in(k)
1751                          ENDDO
1752                       ENDIF
1753                    ENDIF
1754                 ENDDO
1755              ENDDO
1756           ENDIF 
1757
1758        CASE ( 'im_hf_roof' )
1759           IF ( av == 0 ) THEN
1760              DO  m = 1, surf_usm_h%ns
1761                  i = surf_usm_h%i(m) !+ surf_usm_h%ioff
1762                  j = surf_usm_h%j(m) !+ surf_usm_h%joff
1763                  k = surf_usm_h%k(m) !+ surf_usm_h%koff
1764                  local_pf(i,j,k) = surf_usm_h%iwghf_eb(m)
1765              ENDDO
1766           ENDIF
1767
1768        CASE ( 'im_hf_roof_waste' )
1769           IF ( av == 0 ) THEN
1770              DO m = 1, surf_usm_h%ns
1771                 i = surf_usm_h%i(m) !+ surf_usm_h%ioff
1772                 j = surf_usm_h%j(m) !+ surf_usm_h%joff
1773                 k = surf_usm_h%k(m) !+ surf_usm_h%koff
1774                 local_pf(i,j,k) = surf_usm_h%waste_heat(m)
1775              ENDDO
1776           ENDIF
1777
1778       CASE ( 'im_hf_wall_win' )
1779           IF ( av == 0 ) THEN
1780              DO l = 0, 3
1781                 DO m = 1, surf_usm_v(l)%ns
1782                    i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
1783                    j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
1784                    k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
1785                    local_pf(i,j,k) = surf_usm_v(l)%iwghf_eb(m)
1786                 ENDDO
1787              ENDDO
1788           ENDIF
1789
1790        CASE ( 'im_hf_wall_win_waste' )
1791           IF ( av == 0 ) THEN
1792              DO l = 0, 3
1793                 DO m = 1, surf_usm_v(l)%ns
1794                    i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
1795                    j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
1796                    k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
1797                    local_pf(i,j,k) =  surf_usm_v(l)%waste_heat(m) 
1798                 ENDDO
1799              ENDDO
1800           ENDIF
1801
1802!
1803!< NOTE im_t_indoor_roof and im_t_indoor_wall_win not work yet
1804
1805!         CASE ( 'im_t_indoor_roof' )
1806!            IF ( av == 0 ) THEN
1807!               DO  m = 1, surf_usm_h%ns
1808!                   i = surf_usm_h%i(m) !+ surf_usm_h%ioff
1809!                   j = surf_usm_h%j(m) !+ surf_usm_h%joff
1810!                   k = surf_usm_h%k(m) !+ surf_usm_h%koff
1811!                   local_pf(i,j,k) = surf_usm_h%t_indoor(m)
1812!               ENDDO
1813!            ENDIF
1814!
1815!         CASE ( 'im_t_indoor_wall_win' )
1816!            IF ( av == 0 ) THEN           
1817!               DO l = 0, 3
1818!                  DO m = 1, surf_usm_v(l)%ns
1819!                     i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
1820!                     j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
1821!                     k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
1822!                     local_pf(i,j,k) = surf_usm_v(l)%t_indoor(m)
1823!                  ENDDO
1824!               ENDDO
1825!            ENDIF
1826
1827        CASE DEFAULT
1828           found = .FALSE.
1829           
1830    END SELECT   
1831
1832 END SUBROUTINE im_data_output_3d         
1833!------------------------------------------------------------------------------!
1834! Description:
1835! ------------
1836!> Parin for &indoor_parameters for indoor model
1837!------------------------------------------------------------------------------!
1838 SUBROUTINE im_parin
1839   
1840    USE control_parameters,                                                    &
1841        ONLY:  indoor_model
1842   
1843    IMPLICIT NONE
1844
1845    CHARACTER (LEN=80) ::  line  !< string containing current line of file PARIN
1846
1847    NAMELIST /indoor_parameters/  dt_indoor, initial_indoor_temperature
1848
1849!
1850!-- Try to find indoor model package
1851    REWIND ( 11 )
1852    line = ' '
1853    DO   WHILE ( INDEX( line, '&indoor_parameters' ) == 0 )
1854       READ ( 11, '(A)', END=10 )  line
1855    ENDDO
1856    BACKSPACE ( 11 )
1857
1858!
1859!-- Read user-defined namelist
1860    READ ( 11, indoor_parameters )
1861!
1862!-- Set flag that indicates that the indoor model is switched on
1863    indoor_model = .TRUE.
1864
1865!
1866!--    Activate spinup (maybe later
1867!        IF ( spinup_time > 0.0_wp )  THEN
1868!           coupling_start_time = spinup_time
1869!           end_time = end_time + spinup_time
1870!           IF ( spinup_pt_mean == 9999999.9_wp )  THEN
1871!              spinup_pt_mean = pt_surface
1872!           ENDIF
1873!           spinup = .TRUE.
1874!        ENDIF
1875
1876 10 CONTINUE
1877   
1878 END SUBROUTINE im_parin
1879
1880
1881END MODULE indoor_model_mod
Note: See TracBrowser for help on using the repository browser.