Home

Context Navigation

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

Last change on this file since 4845 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 112.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |