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

Last change on this file since 4401 was 4346, checked in by motisi, 5 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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