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

Last change on this file since 4329 was 4329, checked in by motisi, 20 months ago

Renamed wall_flags_0 to wall_flags_static_0

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