Home

Context Navigation

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

Last change on this file since 4842 was 4842, checked in by raasch, 3 years ago
reading of namelist file and actions in case of namelist errors revised so that statement labels and goto statements are not required any more, deprecated namelists removed
Property svn:keywords set to `Id`
File size: 111.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |