Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |