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

Last change on this file since 4843 was 4843, checked in by raasch, 7 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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