source: palm/trunk/SOURCE/indoor_model_mod.f90

Last change on this file was 4850, checked in by suehring, 7 months ago

Bugfix in indoor model: consider previous indoor temperature during restarts; Further bugfix in mpi-io restart mechanism for the waste-heat flux from buildings

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