Home

Context Navigation

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

Last change on this file since 3497 was 3469, checked in by kanani, 6 years ago
Implement indoor climate and energy demand model
Property svn:keywords set to `Id`
File size: 69.5 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
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2018-2018 Leibniz Universitaet Hannover
18	! Copyright 2018-2018 Hochschule Offenburg
19	!--------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! -----------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: indoor_model_mod.f90 3469 2018-10-30 20:05:07Z suehring $
28	! Initial revision (tlang, suehring, kanani, srissman)
29	!
30	!
31	!
32	! Authors:
33	! --------
34	! @author Tobias Lang
35	! @author Jens Pfafferott
36	! @author Farah Kanani-Suehring
37	! @author Matthias Suehring
38	! @author Sascha RiÃmann
39	!
40	!
41	! Description:
42	! ------------
43	!> <Description of the new module>
44	!> Module for Indoor Climate Model (ICM)
45	!> The module is based on the DIN EN ISO 13790 with simplified hour-based procedure.
46	!> This model is a equivalent circuit diagram of a three-point RC-model (5R1C).
47	!> This module differ between indoor-air temperature an average temperature of indoor surfaces which make it prossible to determine thermal comfort
48	!> the heat transfer between indoor and outdoor is simplified
49
50	!> @todo Replace window_area_per_facade by %frac(1,m) for window
51	!> @todo emissivity change for window blinds if solar_protection_on=1
52	!> @todo write datas in netcdf file as output data
53	!> @todo reduce the building volume with netto ground surface to take respect costruction areas like walls and ceilings. Have effect on factor_a, factor_c, airchange and lambda_at
54	!>
55	!> @note Do we allow use of integer flags, or only logical flags? (concerns e.g. cooling_on, heating_on)
56	!> @note How to write indoor temperature output to pt array?
57	!>
58	!> @bug <Enter known bugs here>
59	!------------------------------------------------------------------------------!
60	MODULE indoor_model_mod
61
62	USE control_parameters, &
63	ONLY: initializing_actions
64
65	USE kinds
66
67	USE surface_mod, &
68	ONLY: surf_usm_h, surf_usm_v
69
70
71	IMPLICIT NONE
72
73	!
74	!-- Define data structure for buidlings.
75	TYPE build
76
77	INTEGER(iwp) :: id !< building ID
78	INTEGER(iwp) :: kb_min !< lowest vertical index of a building
79	INTEGER(iwp) :: kb_max !< highest vertical index of a building
80	INTEGER(iwp) :: num_facades_per_building_h !< total number of horizontal facades elements
81	INTEGER(iwp) :: num_facades_per_building_h_l !< number of horizontal facade elements on local subdomain
82	INTEGER(iwp) :: num_facades_per_building_v !< total number of vertical facades elements
83	INTEGER(iwp) :: num_facades_per_building_v_l !< number of vertical facade elements on local subdomain
84
85	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: l_v !< index array linking surface-element orientation index
86	!< for vertical surfaces with building
87	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: m_h !< index array linking surface-element index for
88	!< horizontal surfaces with building
89	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: m_v !< index array linking surface-element index for
90	!< vertical surfaces with building
91	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facade_h !< number of horizontal facade elements per buidling
92	!< and height level
93	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facade_v !< number of vertical facades elements per buidling
94	!< and height level
95
96	LOGICAL :: on_pe = .FALSE. !< flag indicating whether a building with certain ID is on local subdomain
97
98	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in !< mean building indoor temperature, height dependent
99	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_l !< mean building indoor temperature on local subdomain, height dependent
100	REAL(wp), DIMENSION(:), ALLOCATABLE :: volume !< total building volume, height dependent
101	REAL(wp), DIMENSION(:), ALLOCATABLE :: vol_frac !< fraction of local on total building volume, height dependent
102	REAL(wp), DIMENSION(:), ALLOCATABLE :: vpf !< building volume volume per facade element, height dependent
103
104	END TYPE build
105
106	TYPE(build), DIMENSION(:), ALLOCATABLE :: buildings !< building array
107
108	INTEGER(iwp) :: num_build !< total number of buildings in domain
109
110	REAL(wp) :: volume_fraction
111
112	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in !< dummy array for indoor temperature for the
113	!< total building volume at each discrete height level
114	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_l !< dummy array for indoor temperature for the
115	!< local building volume fraction at each discrete height level
116
117	!
118	!-- Declare all global variables within the module
119
120	INTEGER(iwp) :: building_type = 1 !< namelist parameter with
121	!< X1=construction year (cy) 1950, X2=cy 2000, X3=cy 2050
122	!< R=Residental building, O=Office, RW=Enlarged Windows, P=Panel type (Plattenbau) WBS 70, H=Hospital (in progress), I=Industrial halls (in progress), S=Special Building (in progress)
123	!< (0=R1, 1=R2, 2=R3, 3=O1, 4=O2, 5=O3,...)
124	INTEGER(iwp) :: cooling_on !< Indoor cooling flag (0=off, 1=on)
125	INTEGER(iwp) :: heating_on !< Indoor heating flag (0=off, 1=on)
126	INTEGER(iwp) :: solar_protection_off !< Solar protection off
127	INTEGER(iwp) :: solar_protection_on !< Solar protection on
128
129	REAL(wp) :: air_change_high !< [1/h] air changes per time_utc_hour
130	REAL(wp) :: air_change_low !< [1/h] air changes per time_utc_hour
131	REAL(wp) :: eff_mass_area !< [mÂ²] the effective mass-related area
132	REAL(wp) :: floor_area_per_facade !< [mÂ²] net floor area (Sum of all floors)
133	REAL(wp) :: total_area !<! [mÂ²] area of all surfaces pointing to zone
134	REAL(wp) :: window_area_per_facade !< [m2] window area per facade element
135	REAL(wp) :: air_change !< [1/h] Airflow
136	REAL(wp) :: bldg_part_surf_i = 4 !< [mÂ²_surf,i] part building surface, from Palm, das mÃŒsste mittlerweile "facade_element_area" sein!
137	REAL(wp) :: facade_element_area !< [mÂ²_facade] building surface facade
138	REAL(wp) :: indoor_volume_per_facade !< [mÂ³] indoor air volume per facade element
139	REAL(wp) :: c_m !< [J/K] internal heat storage capacity
140	REAL(wp) :: dt_indoor = 3600.0_wp !< [s] namelist parameter: time interval for indoor-model application
141	REAL(wp) :: eta_ve !< [-] heat recovery efficiency
142	REAL(wp) :: f_c_win !< [-] shading factor
143	REAL(wp) :: factor_a !< [-] Dynamic parameters specific effective surface according to Table 12; 2.5 (very light, light and medium), 3.0 (heavy), 3.5 (very heavy)
144	REAL(wp) :: factor_c !< [J/(m2 K)] Dynamic parameters inner heatstorage according to Table 12; 80000 (very light), 110000 (light), 165000 (medium), 260000 (heavy), 370000 (very heavy)
145	REAL(wp) :: g_value_win !< [-] SHGC factor
146	REAL(wp) :: h_tr_1 !<! [W/K] Heat transfer coefficient auxiliary variable 1
147	REAL(wp) :: h_tr_2 !<! [W/K] Heat transfer coefficient auxiliary variable 2
148	REAL(wp) :: h_tr_3 !<! [W/K] Heat transfer coefficient auxiliary variable 3
149	REAL(wp) :: h_tr_em !<! [W/K] Heat transfer coefficient of the emmision (got with h_tr_ms the thermal mass)
150	REAL(wp) :: h_tr_is !<! [W/K] thermal coupling conductance (Thermischer Kopplungsleitwert)
151	REAL(wp) :: h_tr_ms !<! [W/K] Heat transfer conductance term (got with h_tr_em the thermal mass)
152	REAL(wp) :: h_tr_op !<! [W/K] heat transfer coefficient of opaque components (assumption: got all thermal mass) contains of h_tr_em and h_tr_ms
153	REAL(wp) :: h_tr_w !<! [W/K] heat transfer coefficient of doors, windows, curtain walls and glazed walls (assumption: thermal mass=0)
154	REAL(wp) :: h_ve !<! [W/K] heat transfer of ventilation
155	REAL(wp) :: height_storey !< [m] storey heigth
156	REAL(wp) :: height_cei_con !< [m] ceiling construction heigth
157	REAL(wp) :: initial_indoor_temperature !< namelist parameter
158	REAL(wp) :: lambda_at !< [-] ratio internal surface/floor area chap. 7.2.2.2.
159	REAL(wp) :: lambda_layer3 !< [W/(m*K)] Thermal conductivity of the inner layer
160	REAL(wp) :: net_sw_in !< net short-wave radiation (in - out; was i_global --> CORRECT?)
161	REAL(wp) :: qint_high !< [W/m2] internal heat gains, option Database qint_0-23
162	REAL(wp) :: qint_low !< [W/m2] internal heat gains, option Database qint_0-23
163	REAL(wp) :: phi_c_max !< [W] Max. Cooling capacity (negative)
164	REAL(wp) :: phi_h_max !< [W] Max. Heating capacity (negative)
165	REAL(wp) :: phi_hc_nd !<! [W] heating demand and/or cooling demand
166	REAL(wp) :: phi_hc_nd_10 !<! [W] heating demand and/or cooling demand for heating or cooling
167	REAL(wp) :: phi_hc_nd_ac !<! [W] actual heating demand and/or cooling demand
168	REAL(wp) :: phi_hc_nd_un !<! [W] unlimited heating demand and/or cooling demand which is necessary to reach the demanded required temperature (heating is positive, cooling is negative)
169	REAL(wp) :: phi_ia !< [W] internal air load = internal loads * 0.5, Eq. (C.1)
170	REAL(wp) :: phi_m !<! [W] mass specific thermal load (internal and external)
171	REAL(wp) :: phi_mtot !<! [W] total mass specific thermal load (internal and external)
172	REAL(wp) :: phi_sol !< [W] solar loads
173	REAL(wp) :: phi_st !<! [W] mass specific thermal load implied non thermal mass
174	REAL(wp) :: q_emission !< emissions, in first version = 0, option for second part of the project
175	REAL(wp) :: q_wall_win !< heat flux from indoor into wall/window
176	REAL(wp) :: q_waste_heat !< waste heat, sum of waste heat over the roof to Palm
177	REAL(wp) :: q_waste_heat_bldg !< [W/building] waste heat of the complete building, in Palm sum of all indoor_model-calculations
178	REAL(wp) :: s_layer3 !< [m] half thickness of the inner layer (layer_3)
179	REAL(wp) :: schedule_d !< activation for internal loads (low or high + low)
180	REAL(wp) :: skip_time_do_indoor = 0.0_wp !< [s] Indoor model is not called before this time
181	REAL(wp) :: theta_air !<! [Â°C] air temperature of the RC-node
182	REAL(wp) :: theta_air_0 !<! [oC] air temperature of the RC-node in equilibrium
183	REAL(wp) :: theta_air_10 !<! [oC] air temperature of the RC-node from a heating capacity of 10 W/mÂ²
184	REAL(wp) :: theta_air_ac !< [oC] actual room temperature after heating/cooling
185	REAL(wp) :: theta_air_set !< [oC] Setpoint_temperature for the room
186	REAL(wp) :: theta_int_c_set !< [oC] Max. Setpoint temperature (summer)
187	REAL(wp) :: theta_int_h_set !< [oC] Max. Setpoint temperature (winter)
188	REAL(wp) :: theta_m !<! [oC} inner temperature of the RC-node
189	REAL(wp) :: theta_m_t !<! [oC] (Fictive) component temperature timestep
190	REAL(wp) :: theta_m_t_prev !< [oC] (Fictive) component temperature previous timestep (do not change)
191	REAL(wp) :: theta_op !< [oC] operative temperature
192	REAL(wp) :: theta_s !<! [oC] surface temperature of the RC-node
193	REAL(wp) :: time_indoor = 0.0_wp !< [s] time since last call of indoor model
194	REAL(wp) :: time_utc_hour !< Time in hours per day (UTC)
195	REAL(wp) :: u_value_win !< [W/(m2*K)] transmittance
196	REAL(wp) :: ventilation_int_loads !< Zuteilung der GebÃ€ude fÃŒr Verlauf/AktivitÃ€t der LÃŒftung und internen Lasten
197
198	!
199	!-- Declare all global parameters within the module
200	REAL(wp), PARAMETER :: params_f_f = 0.3_wp !< [-] frame ratio chap. 8.3.2.1.1 for buildings with mostly cooling 2.0_wp
201	REAL(wp), PARAMETER :: params_f_w = 0.9_wp !< [-] correction factor (fuer nicht senkrechten Stahlungseinfall DIN 4108-2 chap.8, (hier konstant, keine WinkelabhÃ€ngigkeit)
202	REAL(wp), PARAMETER :: params_f_win = 0.5_wp !< [-] proportion of window area, Database A_win aus Datenbank 27 window_area_per_facade_percent
203	REAL(wp), PARAMETER :: params_solar_protection = 300.0_wp !< [W/m2] chap. G.5.3.1 sun protection closed, if the radiation on facade exceeds this value
204	REAL(wp), PARAMETER :: params_waste_heat_c = 4.0_wp !< [-] anthropogenic heat outputs for cooling e.g. 4 for KKM with COP = 3
205	REAL(wp), PARAMETER :: params_waste_heat_h = 1.111_wp !< [-] anthropogenic heat outputs for heating e.g. 1 / 0.9 = 1.111111 for combustion with eta = 0.9 or -3 for WP with COP = 4
206	REAL(wp), PARAMETER :: h_is = 3.45_wp !< [W/(m^2 K)] h_is = 3.45 between surface and air (chap. 7.2.2.2)
207	REAL(wp), PARAMETER :: h_ms = 9.1_wp !< [W/K] h_ms = 9.10 W / (m2 K) between component and surface (chap. 12.2.2)
208
209	SAVE
210
211
212	PRIVATE
213
214	!
215	!-- Add INTERFACES that must be available to other modules
216	PUBLIC im_init, im_main_heatcool, im_parin
217
218	!
219	!-- Add VARIABLES that must be available to other modules
220	PUBLIC dt_indoor, skip_time_do_indoor, time_indoor
221
222	!
223	!-- Calculations for indoor temperatures
224	INTERFACE im_calc_temperatures
225	MODULE PROCEDURE im_calc_temperatures
226	END INTERFACE im_calc_temperatures
227	!
228	!-- Initialization actions
229	INTERFACE im_init
230	MODULE PROCEDURE im_init
231	END INTERFACE im_init
232
233	!
234	!-- Main part of indoor model
235	INTERFACE im_main_heatcool
236	MODULE PROCEDURE im_main_heatcool
237	END INTERFACE im_main_heatcool
238	!
239	!-- Reading of NAMELIST parameters
240	INTERFACE im_parin
241	MODULE PROCEDURE im_parin
242	END INTERFACE im_parin
243
244	CONTAINS
245
246	!------------------------------------------------------------------------------!
247	! Description:
248	! ------------
249	!< Calculation of the air temperatures and mean radiation temperature
250	!< This is basis for the operative temperature
251	!< Based on a Crank-Nicholson scheme with a timestep of a hour
252	!------------------------------------------------------------------------------!
253	SUBROUTINE im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
254	near_facade_temperature, phi_hc_nd_dummy )
255
256	USE arrays_3d, &
257	ONLY: pt
258
259
260	IMPLICIT NONE
261
262
263	INTEGER(iwp) :: i
264	INTEGER(iwp) :: j
265	INTEGER(iwp) :: k
266
267	REAL(wp) :: indoor_wall_window_temperature !< weighted temperature of innermost wall/window layer
268	REAL(wp) :: near_facade_temperature
269	REAL(wp) :: phi_hc_nd_dummy
270
271	!< Calculation of total mass specific thermal load (internal and external)
272	phi_mtot = ( phi_m + h_tr_em * indoor_wall_window_temperature &
273	+ h_tr_3 * ( phi_st + h_tr_w * pt(k,j,i) &
274	+ h_tr_1 * &
275	( ( ( phi_ia + phi_hc_nd_dummy ) / h_ve ) &
276	+ near_facade_temperature ) &
277	) / h_tr_2 &
278	) !< [oC] Eq. (C.5)
279
280	!< Calculation of component temperature at factual timestep
281	theta_m_t = ( ( theta_m_t_prev &
282	* ( ( c_m / 3600 ) - 0.5 * ( h_tr_3 + h_tr_em ) ) + phi_mtot &
283	) &
284	/ ( ( c_m / 3600 ) + 0.5 * ( h_tr_3 + h_tr_em ) ) &
285	) !< [oC] Eq. (C.4)
286
287	!< Calculation of mean inner temperature for the RC-node in actual timestep
288	theta_m = ( theta_m_t + theta_m_t_prev ) * 0.5 !< [oC] Eq. (C.9)
289
290	!< Calculation of mean surface temperature of the RC-node in actual timestep
291	theta_s = ( ( h_tr_ms * theta_m + phi_st + h_tr_w * pt(k,j,i) &
292	+ h_tr_1 * ( near_facade_temperature + ( phi_ia + phi_hc_nd_dummy ) / h_ve ) &
293	) &
294	/ ( h_tr_ms + h_tr_w + h_tr_1 ) &
295	) !< [oC] Eq. (C.10)
296
297	!< Calculation of the air temperature of the RC-node
298	theta_air = ( h_tr_is * theta_s + h_ve * near_facade_temperature &
299	+ phi_ia + phi_hc_nd_dummy ) / ( h_tr_is + h_ve ) !< [oC] Eq. (C.11)
300
301	END SUBROUTINE im_calc_temperatures
302
303	!------------------------------------------------------------------------------!
304	! Description:
305	! ------------
306	!> Initialization of the indoor model.
307	!> Static information are calculated here, e.g. building parameters and
308	!> geometrical information, everything that doesn't change in time.
309	!
310	!-- Input values
311	!-- Input datas from Palm, M4
312	! i_global --> net_sw_in !global radiation [W/m2]
313	! theta_e --> pt(k,j,i) !undisturbed outside temperature, 1. PALM volume, for windows
314	! theta_sup = theta_f --> surf_usm_h%t_surf_10cm(m)
315	! surf_usm_v(l)%t_surf_10cm(m) !Air temperature, facade near (10cm) air temperature from 1. Palm volume
316	! theta_node --> t_wall_h(nzt_wall,m)
317	! t_wall_v(l)%t(nzt_wall,m) !Temperature of innermost wall layer, for opaque wall
318	!------------------------------------------------------------------------------!
319	SUBROUTINE im_init
320
321	USE arrays_3d, &
322	ONLY: dzw
323
324	USE control_parameters, &
325	ONLY: message_string
326
327	USE indices, &
328	ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
329
330	USE grid_variables, &
331	ONLY: dx, dy
332
333	USE netcdf_data_input_mod, &
334	ONLY: building_id_f
335
336	USE pegrid
337
338	USE surface_mod, &
339	ONLY: surf_usm_h, surf_usm_v
340
341	USE urban_surface_mod, &
342	ONLY: building_pars, building_type
343
344	IMPLICIT NONE
345
346	INTEGER(iwp) :: fa !< running index for facade elements of each building
347	INTEGER(iwp) :: i !< running index along x-direction
348	INTEGER(iwp) :: j !< running index along y-direction
349	INTEGER(iwp) :: k !< running index along z-direction
350	INTEGER(iwp) :: l !< running index for surface-element orientation
351	INTEGER(iwp) :: m !< running index surface elements
352	INTEGER(iwp) :: n !< building index
353	INTEGER(iwp) :: nb !< building index
354
355	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids !< building IDs on entire model domain
356	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final !< building IDs on entire model domain,
357	!< multiple occurences are sorted out
358	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final_tmp !< temporary array used for resizing
359	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l !< building IDs on local subdomain
360	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l_tmp !< temporary array used to resize array of building IDs
361	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displace_dum !< displacements of start addresses, used for MPI_ALLGATHERV
362	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_max_l !< highest vertical index of a building on subdomain
363	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_min_l !< lowest vertical index of a building on subdomain
364	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: n_fa !< counting array
365	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facades_h !< dummy array used for summing-up total number of
366	!< horizontal facade elements
367	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facades_v !< dummy array used for summing-up total number of
368	!< vertical facade elements
369	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: receive_dum_h !< dummy array used for MPI_ALLREDUCE
370	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: receive_dum_v !< dummy array used for MPI_ALLREDUCE
371
372	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings !< number of buildings with different ID on entire model domain
373	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings_l !< number of buildings with different ID on local subdomain
374
375	REAL(wp), DIMENSION(:), ALLOCATABLE :: local_weight !< dummy representing fraction of local on total building volume,
376	!< height dependent
377	REAL(wp), DIMENSION(:), ALLOCATABLE :: volume !< total building volume at each discrete height level
378	REAL(wp), DIMENSION(:), ALLOCATABLE :: volume_l !< total building volume at each discrete height level,
379	!< on local subdomain
380
381	!
382	!-- Initializing of indoor model is only possible if buildings can be
383	!-- distinguished by their IDs.
384	IF ( .NOT. building_id_f%from_file ) THEN
385	message_string = 'Indoor model requires information about building_id'
386	CALL message( 'im_init', 'PA0999', 1, 2, 0, 6, 0 )
387	ENDIF
388	!
389	!-- Determine number of different building IDs on local subdomain.
390	num_buildings_l = 0
391	num_buildings = 0
392	ALLOCATE( build_ids_l(1) )
393	DO i = nxl, nxr
394	DO j = nys, nyn
395	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
396	IF ( num_buildings_l(myid) > 0 ) THEN
397	IF ( ANY( building_id_f%var(j,i) .EQ. build_ids_l ) ) THEN
398	CYCLE
399	ELSE
400	num_buildings_l(myid) = num_buildings_l(myid) + 1
401	!
402	!-- Resize array with different local building ids
403	ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
404	build_ids_l_tmp = build_ids_l
405	DEALLOCATE( build_ids_l )
406	ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
407	build_ids_l(1:num_buildings_l(myid)-1) = &
408	build_ids_l_tmp(1:num_buildings_l(myid)-1)
409	build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
410	DEALLOCATE( build_ids_l_tmp )
411	ENDIF
412	!
413	!-- First occuring building id on PE
414	ELSE
415	num_buildings_l(myid) = num_buildings_l(myid) + 1
416	build_ids_l(1) = building_id_f%var(j,i)
417	ENDIF
418	ENDIF
419	ENDDO
420	ENDDO
421	!
422	!-- Determine number of building IDs for the entire domain. (Note, building IDs
423	!-- can appear multiple times as buildings might be distributed over several
424	!-- PEs.)
425	#if defined( __parallel )
426	CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs, &
427	MPI_INTEGER, MPI_SUM, comm2d, ierr )
428	#else
429	num_buildings = num_buildings_l
430	#endif
431	ALLOCATE( build_ids(1:SUM(num_buildings)) )
432	!
433	!-- Gather building IDs. Therefore, first, determine displacements used
434	!-- required for MPI_GATHERV call.
435	ALLOCATE( displace_dum(0:numprocs-1) )
436	displace_dum(0) = 0
437	DO i = 1, numprocs-1
438	displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
439	ENDDO
440
441	#if defined( __parallel )
442	CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)), &
443	num_buildings(myid), &
444	MPI_INTEGER, &
445	build_ids, &
446	num_buildings, &
447	displace_dum, &
448	MPI_INTEGER, &
449	comm2d, ierr )
450
451	DEALLOCATE( displace_dum )
452
453	#else
454	build_ids = build_ids_l
455	#endif
456	!
457	!-- Note: in parallel mode, building IDs can occur mutliple times, as
458	!-- each PE has send its own ids. Therefore, sort out building IDs which
459	!-- appear multiple times.
460	num_build = 0
461	DO n = 1, SIZE(build_ids)
462
463	IF ( ALLOCATED(build_ids_final) ) THEN
464	IF ( ANY( build_ids(n) .EQ. build_ids_final ) ) THEN !FK: Warum ANY?, Warum .EQ.? --> s.o
465	CYCLE
466	ELSE
467	num_build = num_build + 1
468	!
469	!-- Resize
470	ALLOCATE( build_ids_final_tmp(1:num_build) )
471	build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
472	DEALLOCATE( build_ids_final )
473	ALLOCATE( build_ids_final(1:num_build) )
474	build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
475	build_ids_final(num_build) = build_ids(n)
476	DEALLOCATE( build_ids_final_tmp )
477	ENDIF
478	ELSE
479	num_build = num_build + 1
480	ALLOCATE( build_ids_final(1:num_build) )
481	build_ids_final(num_build) = build_ids(n)
482	ENDIF
483	ENDDO
484
485	!
486	!-- Allocate building-data structure array. Note, this is a global array
487	!-- and all building IDs on domain are known by each PE. Further attributes,
488	!-- e.g. height-dependent arrays, however, are only allocated on PEs where
489	!-- the respective building is present (in order to reduce memory demands).
490	ALLOCATE( buildings(1:num_build) )
491	!
492	!-- Store building IDs and check if building with certain ID is present on
493	!-- subdomain.
494	DO nb = 1, num_build
495	buildings(nb)%id = build_ids_final(nb)
496
497	IF ( ANY( building_id_f%var == buildings(nb)%id ) ) &
498	buildings(nb)%on_pe = .TRUE.
499	ENDDO
500	!
501	!-- Determine the maximum vertical dimension occupied by each building.
502	ALLOCATE( k_min_l(1:num_build) )
503	ALLOCATE( k_max_l(1:num_build) )
504	k_min_l = nzt + 1
505	k_max_l = 0
506
507	DO i = nxl, nxr
508	DO j = nys, nyn
509	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
510	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
511	DIM = 1 )
512	DO k = nzb+1, nzt+1
513	!
514	!-- Check if grid point belongs to a building.
515	IF ( BTEST( wall_flags_0(k,j,i), 6 ) ) THEN
516	k_min_l(nb) = MIN( k_min_l(nb), k )
517	k_max_l(nb) = MAX( k_max_l(nb), k )
518	ENDIF
519
520	ENDDO
521	ENDIF
522	ENDDO
523	ENDDO
524
525	DO nb = 1, num_build
526	#if defined( __parallel )
527	CALL MPI_ALLREDUCE( k_min_l(nb), buildings(nb)%kb_min, 1, MPI_INTEGER, &
528	MPI_MIN, comm2d, ierr )
529	CALL MPI_ALLREDUCE( k_max_l(nb), buildings(nb)%kb_max, 1, MPI_INTEGER, &
530	MPI_MAX, comm2d, ierr )
531	#else
532	buildings(nb)%kb_min = k_min_l(nb)
533	buildings(nb)%kb_max = k_max_l(nb)
534	#endif
535
536	ENDDO
537
538	DEALLOCATE( k_min_l )
539	DEALLOCATE( k_max_l )
540	!
541	!-- Calculate building volume
542	DO nb = 1, num_build
543	!
544	!-- Allocate temporary array for summing-up building volume
545	ALLOCATE( volume(buildings(nb)%kb_min:buildings(nb)%kb_max) )
546	ALLOCATE( volume_l(buildings(nb)%kb_min:buildings(nb)%kb_max) )
547	volume = 0.0_wp
548	volume_l = 0.0_wp
549	!
550	!-- Calculate building volume per height level on each PE where
551	!-- these building is present.
552	IF ( buildings(nb)%on_pe ) THEN
553	ALLOCATE( buildings(nb)%volume(buildings(nb)%kb_min:buildings(nb)%kb_max) )
554	ALLOCATE( buildings(nb)%vol_frac(buildings(nb)%kb_min:buildings(nb)%kb_max) )
555	buildings(nb)%volume = 0.0_wp
556	buildings(nb)%vol_frac = 0.0_wp
557
558	IF ( ANY( building_id_f%var == buildings(nb)%id ) ) THEN
559	DO i = nxl, nxr
560	DO j = nys, nyn
561	DO k = buildings(nb)%kb_min, buildings(nb)%kb_max
562	IF ( building_id_f%var(j,i) /= building_id_f%fill ) &
563	volume_l(k) = dx * dy * dzw(k)
564	ENDDO
565	ENDDO
566	ENDDO
567	ENDIF
568	ENDIF
569	!
570	!-- Sum-up building volume from all subdomains
571	#if defined( __parallel )
572	CALL MPI_ALLREDUCE( volume_l, volume, SIZE(volume), MPI_REAL, MPI_SUM, &
573	comm2d, ierr )
574	#else
575	volume = volume_l
576	#endif
577	!
578	!-- Save total building volume as well as local fraction on volume on
579	!-- building data structure.
580	IF ( ALLOCATED( buildings(nb)%volume ) ) buildings(nb)%volume = volume
581	!
582	!-- Determine fraction of local on total building volume
583	IF ( buildings(nb)%on_pe ) buildings(nb)%vol_frac = volume_l / volume
584
585	DEALLOCATE( volume )
586	DEALLOCATE( volume_l )
587
588	ENDDO
589
590	!
591	!-- Allocate arrays for indoor temperature.
592	DO nb = 1, num_build
593	IF ( buildings(nb)%on_pe ) THEN
594	ALLOCATE( buildings(nb)%t_in(buildings(nb)%kb_min:buildings(nb)%kb_max) )
595	ALLOCATE( buildings(nb)%t_in_l(buildings(nb)%kb_min:buildings(nb)%kb_max) )
596	buildings(nb)%t_in = 0.0_wp
597	buildings(nb)%t_in_l = 0.0_wp
598	ENDIF
599	ENDDO
600	!
601	!-- Allocate arrays for number of facades per height level. Distinguish between
602	!-- horizontal and vertical facades.
603	DO nb = 1, num_build
604	IF ( buildings(nb)%on_pe ) THEN
605	ALLOCATE( buildings(nb)%num_facade_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
606	ALLOCATE( buildings(nb)%num_facade_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
607
608	buildings(nb)%num_facade_h = 0
609	buildings(nb)%num_facade_v = 0
610	ENDIF
611	ENDDO
612	!
613	!-- Determine number of facade elements per building on local subdomain.
614	!-- Distinguish between horizontal and vertical facade elements.
615	!
616	!-- Horizontal facades
617	buildings(:)%num_facades_per_building_h_l = 0
618	DO m = 1, surf_usm_h%ns
619	!
620	!-- For the current facade element determine corresponding building index.
621	!-- First, obtain j,j,k indices of the building. Please note the
622	!-- offset between facade/surface element and building location (for
623	!-- horizontal surface elements the horizontal offsets are zero).
624	i = surf_usm_h%i(m) + surf_usm_h%ioff
625	j = surf_usm_h%j(m) + surf_usm_h%joff
626	k = surf_usm_h%k(m) + surf_usm_h%koff
627	!
628	!-- Determine building index and check whether building is on PE
629	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
630	IF ( buildings(nb)%on_pe ) THEN
631	!
632	!-- Count number of facade elements at each height level.
633	buildings(nb)%num_facade_h(k) = buildings(nb)%num_facade_h(k) + 1
634	!
635	!-- Moreover, sum up number of local facade elements per building.
636	buildings(nb)%num_facades_per_building_h_l = &
637	buildings(nb)%num_facades_per_building_h_l + 1
638	ENDIF
639	ENDDO
640	!
641	!-- Vertical facades
642	buildings(:)%num_facades_per_building_v_l = 0
643	DO l = 0, 3
644	DO m = 1, surf_usm_v(l)%ns
645	!
646	!-- For the current facade element determine corresponding building index.
647	!-- First, obtain j,j,k indices of the building. Please note the
648	!-- offset between facade/surface element and building location (for
649	!-- vertical surface elements the vertical offsets are zero).
650	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
651	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
652	k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
653
654	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
655	DIM = 1 )
656	IF ( buildings(nb)%on_pe ) THEN
657	buildings(nb)%num_facade_v(k) = buildings(nb)%num_facade_v(k) + 1
658	buildings(nb)%num_facades_per_building_v_l = &
659	buildings(nb)%num_facades_per_building_v_l + 1
660	ENDIF
661	ENDDO
662	ENDDO
663
664	!
665	!-- Determine total number of facade elements per building and assign number to
666	!-- building data type.
667	DO nb = 1, num_build
668	!
669	!-- Allocate dummy array used for summing-up facade elements.
670	!-- Please note, dummy arguments are necessary as building-date type
671	!-- arrays are not necessarily allocated on all PEs.
672	ALLOCATE( num_facades_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
673	ALLOCATE( num_facades_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
674	ALLOCATE( receive_dum_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
675	ALLOCATE( receive_dum_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
676	num_facades_h = 0
677	num_facades_v = 0
678	receive_dum_h = 0
679	receive_dum_v = 0
680
681	IF ( buildings(nb)%on_pe ) THEN
682	num_facades_h = buildings(nb)%num_facade_h
683	num_facades_v = buildings(nb)%num_facade_v
684	ENDIF
685
686	#if defined( __parallel )
687	CALL MPI_ALLREDUCE( num_facades_h, &
688	receive_dum_h, &
689	buildings(nb)%kb_max - buildings(nb)%kb_min + 1, &
690	MPI_INTEGER, &
691	MPI_SUM, &
692	comm2d, &
693	ierr )
694
695	CALL MPI_ALLREDUCE( num_facades_v, &
696	receive_dum_v, &
697	buildings(nb)%kb_max - buildings(nb)%kb_min + 1, &
698	MPI_INTEGER, &
699	MPI_SUM, &
700	comm2d, &
701	ierr )
702	IF ( ALLOCATED( buildings(nb)%num_facade_h ) ) & !FK: Was wenn not allocated? --> s.o.
703	buildings(nb)%num_facade_h = receive_dum_h
704	IF ( ALLOCATED( buildings(nb)%num_facade_v ) ) &
705	buildings(nb)%num_facade_v = receive_dum_v
706	#else
707	buildings(nb)%num_facade_h = num_facades_h
708	buildings(nb)%num_facade_v = num_facades_v
709	#endif
710	!
711	!-- Deallocate dummy arrays
712	DEALLOCATE( num_facades_h )
713	DEALLOCATE( num_facades_v )
714	DEALLOCATE( receive_dum_h )
715	DEALLOCATE( receive_dum_v )
716	!
717	!-- Allocate index arrays which link facade elements with surface-data type.
718	!-- Please note, no height levels are considered here (information is stored
719	!-- in surface-data type itself).
720	IF ( buildings(nb)%on_pe ) THEN
721	!
722	!-- Determine number of facade elements per building.
723	buildings(nb)%num_facades_per_building_h = SUM( buildings(nb)%num_facade_h )
724	buildings(nb)%num_facades_per_building_v = SUM( buildings(nb)%num_facade_v )
725	!
726	!-- Allocate arrays which link the building with the horizontal and vertical
727	!-- urban-type surfaces. Please note, linking arrays are allocated over all
728	!-- facade elements, which is required in case a building is located at the
729	!-- subdomain boundaries, where the building and the corresponding surface
730	!-- elements are located on different subdomains.
731	ALLOCATE( buildings(nb)%m_h(1:buildings(nb)%num_facades_per_building_h_l) )
732
733	ALLOCATE( buildings(nb)%l_v(1:buildings(nb)%num_facades_per_building_v_l) )
734	ALLOCATE( buildings(nb)%m_v(1:buildings(nb)%num_facades_per_building_v_l) )
735	ENDIF
736	!
737	!-- Determine volume per facade element (vpf)
738	IF ( buildings(nb)%on_pe ) THEN
739	ALLOCATE( buildings(nb)%vpf(buildings(nb)%kb_min:buildings(nb)%kb_max) )
740
741	DO k = buildings(nb)%kb_min, buildings(nb)%kb_max
742	buildings(nb)%vpf(k) = buildings(nb)%volume(k) / &
743	( buildings(nb)%num_facade_h(k) + &
744	buildings(nb)%num_facade_v(k) )
745	ENDDO
746	ENDIF
747	ENDDO
748	!
749	!-- Link facade elements with surface data type.
750	!-- Allocate array for counting.
751	ALLOCATE( n_fa(1:num_build) )
752	n_fa = 1
753
754	DO m = 1, surf_usm_h%ns
755	i = surf_usm_h%i(m) + surf_usm_h%ioff
756	j = surf_usm_h%j(m) + surf_usm_h%joff
757
758	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
759
760	buildings(nb)%m_h(n_fa(nb)) = m
761	n_fa(nb) = n_fa(nb) + 1
762	ENDDO
763
764	n_fa = 1
765	DO l = 0, 3
766	DO m = 1, surf_usm_v(l)%ns
767	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
768	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
769
770	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
771
772	buildings(nb)%l_v(n_fa(nb)) = l
773	buildings(nb)%m_v(n_fa(nb)) = m
774	n_fa(nb) = n_fa(nb) + 1
775	ENDDO
776	ENDDO
777	DEALLOCATE( n_fa )
778
779	!
780	!-- Building parameters by type of building. Assigned in urban_surface_mod.f90
781
782	lambda_layer3 = building_pars(63, building_type)
783	s_layer3 = building_pars(57, building_type)
784	f_c_win = building_pars(119, building_type)
785	g_value_win = building_pars(120, building_type)
786	u_value_win = building_pars(121, building_type)
787	air_change_low = building_pars(122, building_type)
788	air_change_high = building_pars(123, building_type)
789	eta_ve = building_pars(124, building_type)
790	factor_a = building_pars(125, building_type)
791	factor_c = building_pars(126, building_type)
792	lambda_at = building_pars(127, building_type)
793	theta_int_h_set = building_pars(118, building_type)
794	theta_int_c_set = building_pars(117, building_type)
795	phi_h_max = building_pars(128, building_type)
796	phi_c_max = building_pars(129, building_type)
797	qint_high = building_pars(130, building_type)
798	qint_low = building_pars(131, building_type)
799	height_storey = building_pars(132, building_type)
800	height_cei_con = building_pars(133, building_type)
801
802	!
803	!-- Setting of initial room temperature [K]
804	!-- (after first loop, use theta_m_t as theta_m_t_prev)
805	theta_m_t_prev = initial_indoor_temperature
806
807
808	END SUBROUTINE im_init
809
810
811	!------------------------------------------------------------------------------!
812	! Description:
813	! ------------
814	!> Main part of the indoor model.
815	!> Calculation of .... (kanani: Please describe)
816	!------------------------------------------------------------------------------!
817	SUBROUTINE im_main_heatcool
818
819	USE arrays_3d, &
820	ONLY: ddzw, dzw
821
822	USE basic_constants_and_equations_mod, &
823	ONLY: c_p
824
825	USE control_parameters, &
826	ONLY: rho_surface
827
828	USE date_and_time_mod, &
829	ONLY: time_utc
830
831	USE grid_variables, &
832	ONLY: dx, dy
833
834	USE pegrid
835
836	USE surface_mod, &
837	ONLY: ind_veg_wall, ind_wat_win, surf_usm_h, surf_usm_v
838
839	USE urban_surface_mod, &
840	ONLY: nzt_wall, t_surf_10cm_h, t_surf_10cm_v, t_wall_h, t_wall_v, &
841	t_window_h, t_window_v, building_type
842
843
844	IMPLICIT NONE
845
846	INTEGER(iwp) :: i !< index of facade-adjacent atmosphere grid point in x-direction
847	INTEGER(iwp) :: j !< index of facade-adjacent atmosphere grid point in y-direction
848	INTEGER(iwp) :: k !< index of facade-adjacent atmosphere grid point in z-direction
849	INTEGER(iwp) :: kk !< vertical index of indoor grid point adjacent to facade
850	INTEGER(iwp) :: l !< running index for surface-element orientation
851	INTEGER(iwp) :: m !< running index surface elements
852	INTEGER(iwp) :: nb !< running index for buildings
853	INTEGER(iwp) :: fa !< running index for facade elements of each building
854
855	REAL(wp) :: indoor_wall_window_temperature !< weighted temperature of innermost wall/window layer
856	REAL(wp) :: near_facade_temperature !< outside air temperature 10cm away from facade
857	REAL(wp) :: time_utc_hour !< time of day (hour UTC)
858
859	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_l_send !< dummy send buffer used for summing-up indoor temperature per kk-level
860	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_recv !< dummy recv buffer used for summing-up indoor temperature per kk-level
861
862	!
863	!-- Daily schedule, here for 08:00-18:00 = 1, other hours = 0.
864	!-- time_utc_hour is calculated here based on time_utc [s] from
865	!-- date_and_time_mod.
866	!-- (kanani: Does this schedule not depend on if it's an office or resident
867	!-- building?)
868	time_utc_hour = time_utc / 3600.0_wp
869
870	!
871	!-- Allocation of the load profiles to the building types
872	!-- Residental Building, panel WBS 70
873	if (building_type == 1 .OR. &
874	building_type == 2 .OR. &
875	building_type == 3 .OR. &
876	building_type == 10 .OR. &
877	building_type == 11 .OR. &
878	building_type == 12) then
879	ventilation_int_loads = 1
880	!-- Office, building with large windows
881	else if (building_type == 4 .OR. &
882	building_type == 5 .OR. &
883	building_type == 6 .OR. &
884	building_type == 7 .OR. &
885	building_type == 8 .OR. &
886	building_type == 9) then
887	ventilation_int_loads = 2
888	!-- Industry, hospitals
889	else if (building_type == 13 .OR. &
890	building_type == 14 .OR. &
891	building_type == 15 .OR. &
892	building_type == 16 .OR. &
893	building_type == 17 .OR. &
894	building_type == 18) then
895	ventilation_int_loads = 3
896
897	end if
898
899	!-- Residental Building, panel WBS 70
900
901	if (ventilation_int_loads == 1) THEN
902	if ( time_utc_hour >= 6.0_wp .AND. time_utc_hour <= 8.0_wp ) THEN
903	schedule_d = 1
904	else if ( time_utc_hour >= 18.0_wp .AND. time_utc_hour <= 23.0_wp ) THEN
905	schedule_d = 1
906	else
907	schedule_d = 0
908	end if
909	end if
910
911	!-- Office, building with large windows
912
913	if (ventilation_int_loads == 2) THEN
914	if ( time_utc_hour >= 8.0_wp .AND. time_utc_hour <= 18.0_wp ) THEN
915	schedule_d = 1
916	else
917	schedule_d = 0
918	end if
919	end if
920
921	!-- Industry, hospitals
922	if (ventilation_int_loads == 3) THEN
923	if ( time_utc_hour >= 6.0_wp .AND. time_utc_hour <= 22.0_wp ) THEN
924	schedule_d = 1
925	else
926	schedule_d = 0
927	end if
928	end if
929
930
931	!
932	!-- Following calculations must be done for each facade element.
933	DO nb = 1, num_build
934	!
935	!-- First, check whether building is present on local subdomain.
936	IF ( buildings(nb)%on_pe ) THEN
937	!
938	!-- Initialize/reset indoor temperature
939	buildings(nb)%t_in = 0.0_wp
940	buildings(nb)%t_in_l = 0.0_wp
941	!
942	!-- Horizontal surfaces
943	DO fa = 1, buildings(nb)%num_facades_per_building_h_l
944	!
945	!-- Determine index where corresponding surface-type information
946	!-- is stored.
947	m = buildings(nb)%m_h(fa)
948	!
949	!-- Determine building height level index.
950	kk = surf_usm_h%k(m) + surf_usm_h%koff
951	!
952	!-- Building geometries --> not time-dependent
953	facade_element_area = dx * dy !< [m2] surface area per facade element
954	floor_area_per_facade = buildings(nb)%vpf(kk) * ddzw(kk) !< [m2] net floor area per facade element
955	indoor_volume_per_facade = buildings(nb)%vpf(kk) !< [m3] indoor air volume per facade element
956	window_area_per_facade = surf_usm_h%frac(ind_wat_win,m) * facade_element_area !< [m2] window area per facade element
957	eff_mass_area = factor_a * floor_area_per_facade !< [m2] standard values according to Table 12 section 12.3.1.2 (calculate over Eq. (65) according to section 12.3.1.2)
958	c_m = factor_c * floor_area_per_facade !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
959	total_area = lambda_at * floor_area_per_facade !< [m2] area of all surfaces pointing to zone Eq. (9) according to section 7.2.2.2
960
961	!-- Calculation of heat transfer coefficient for transmission --> not time-dependent
962	h_tr_w = window_area_per_facade * u_value_win !< [W/K] only for windows
963	h_tr_is = total_area * h_is !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
964	h_tr_ms = eff_mass_area * h_ms !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
965	h_tr_op = 1 / ( 1 / ( ( facade_element_area - window_area_per_facade ) &
966	* lambda_layer3 / s_layer3 * 0.5 ) + 1 / h_tr_ms )
967	h_tr_em = 1 / ( 1 / h_tr_op - 1 / h_tr_ms ) !< [W/K] Eq. (63), Section 12.2.2
968	!
969	!-- internal air loads dependent on the occupacy of the room
970	!-- basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
971	phi_ia = 0.5 * ( ( qint_high * schedule_d + qint_low ) &
972	* floor_area_per_facade ) !< [W] Eq. (C.1)
973	!
974	!-- Airflow dependent on the occupacy of the room
975	!-- basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
976	air_change = ( air_change_high * schedule_d + air_change_low ) !< [1/h]?
977	!
978	!-- Heat transfer of ventilation
979	!-- not less than 0.01 W/K to provide division by 0 in further calculations
980	!-- with heat capacity of air 0.33 Wh/m2K
981	h_ve = MAX( 0.01 , ( air_change * indoor_volume_per_facade * &
982	0.33 * (1 - eta_ve ) ) ) !< [W/K] from ISO 13789 Eq.(10)
983
984	!-- Heat transfer coefficient auxiliary variables
985	h_tr_1 = 1 / ( ( 1 / h_ve ) + ( 1 / h_tr_is ) ) !< [W/K] Eq. (C.6)
986	h_tr_2 = h_tr_1 + h_tr_w !< [W/K] Eq. (C.7)
987	h_tr_3 = 1 / ( ( 1 / h_tr_2 ) + ( 1 / h_tr_ms ) ) !< [W/K] Eq. (C.8)
988	!
989	!-- Net short-wave radiation through window area (was i_global)
990	net_sw_in = surf_usm_h%rad_sw_in(m) - surf_usm_h%rad_sw_out(m)
991	!
992	!-- Quantities needed for im_calc_temperatures
993	i = surf_usm_h%i(m)
994	j = surf_usm_h%j(m)
995	k = surf_usm_h%k(m)
996	near_facade_temperature = t_surf_10cm_h(m)
997	indoor_wall_window_temperature = &
998	surf_usm_h%frac(ind_veg_wall,m) * t_wall_h(nzt_wall,m) &
999	+ surf_usm_h%frac(ind_wat_win,m) * t_window_h(nzt_wall,m)
1000	!
1001	!-- Solar thermal gains. If net_sw_in larger than sun-protection
1002	!-- threshold parameter (params_solar_protection), sun protection will
1003	!-- be activated
1004	IF ( net_sw_in <= params_solar_protection ) THEN
1005	solar_protection_off = 1
1006	solar_protection_on = 0
1007	ELSE
1008	solar_protection_off = 0
1009	solar_protection_on = 1
1010	ENDIF
1011	!
1012	!-- Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
1013	!-- DIN 4108 - 2 chap.8
1014	phi_sol = ( window_area_per_facade * net_sw_in * solar_protection_off &
1015	+ window_area_per_facade * net_sw_in * f_c_win * solar_protection_on ) &
1016	* g_value_win * ( 1 - params_f_f ) * params_f_w !< [W]
1017	!
1018	!-- Calculation of the mass specific thermal load for internal and external heatsources of the inner node
1019	phi_m = (eff_mass_area / total_area) * ( phi_ia + phi_sol ) !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
1020	!
1021	!-- Calculation mass specific thermal load implied non thermal mass
1022	phi_st = ( 1 - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1 * total_area ) ) ) &
1023	* ( phi_ia + phi_sol ) !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
1024	!
1025	!-- Calculations for deriving indoor temperature and heat flux into the wall
1026	!-- Step 1: Indoor temperature without heating and cooling
1027	!-- section C.4.1 Picture C.2 zone 3)
1028	phi_hc_nd = 0
1029
1030	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1031	near_facade_temperature, phi_hc_nd )
1032	!
1033	!-- If air temperature between border temperatures of heating and cooling, assign output variable, then ready
1034	IF ( theta_int_h_set <= theta_air .AND. theta_air <= theta_int_c_set ) THEN
1035	phi_hc_nd_ac = 0
1036	phi_hc_nd = phi_hc_nd_ac
1037	theta_air_ac = theta_air
1038	!
1039	!-- Step 2: Else, apply 10 W/mÂ² heating/cooling power and calculate indoor temperature
1040	!-- again.
1041	ELSE
1042	!
1043	!-- Temperature not correct, calculation method according to section C4.2
1044	theta_air_0 = theta_air !< Note temperature without heating/cooling
1045
1046	!-- Heating or cooling?
1047	IF ( theta_air > theta_int_c_set ) THEN
1048	theta_air_set = theta_int_c_set
1049	ELSE
1050	theta_air_set = theta_int_h_set
1051	ENDIF
1052
1053	!-- Calculate the temperature with phi_hc_nd_10
1054	phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
1055	phi_hc_nd = phi_hc_nd_10
1056
1057	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1058	near_facade_temperature, phi_hc_nd )
1059
1060	theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
1061	!
1062
1063	phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0) &
1064	/ (theta_air_10 - theta_air_0) !< Eq. (C.13)
1065
1066
1067
1068	!-- Step 3: With temperature ratio to determine the heating or cooling capacity
1069	!-- If necessary, limit the power to maximum power
1070	!-- section C.4.1 Picture C.2 zone 2) and 4)
1071	IF ( phi_c_max < phi_hc_nd_un .AND. phi_hc_nd_un < phi_h_max ) THEN
1072	phi_hc_nd_ac = phi_hc_nd_un
1073	phi_hc_nd = phi_hc_nd_un
1074	ELSE
1075	!-- Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
1076	!-- section C.4.1 Picture C.2 zone 1) and 5)
1077	IF ( phi_hc_nd_un > 0 ) THEN
1078	phi_hc_nd_ac = phi_h_max !< Limit heating
1079	ELSE
1080	phi_hc_nd_ac = phi_c_max !< Limit cooling
1081	ENDIF
1082	ENDIF
1083
1084	phi_hc_nd = phi_hc_nd_ac
1085	!
1086	!-- Calculate the temperature with phi_hc_nd_ac (new)
1087	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1088	near_facade_temperature, phi_hc_nd )
1089
1090	theta_air_ac = theta_air
1091
1092	ENDIF
1093	!
1094	!-- Update theta_m_t_prev
1095	theta_m_t_prev = theta_m_t
1096	!
1097	!-- Calculate the operating temperature with weighted mean temperature of air and mean solar temperature
1098	!-- Will be used for thermal comfort calculations
1099	theta_op = 0.3 * theta_air_ac + 0.7 * theta_s !< [Â°C] operative Temperature Eq. (C.12)
1100	!
1101	!-- Heat flux into the wall. Value needed in urban_surface_mod to
1102	!-- calculate heat transfer through wall layers towards the facade
1103	!-- (use c_p * rho_surface to convert [W/m2] into [K m/s])
1104	q_wall_win = h_tr_ms * ( theta_s - theta_m ) &
1105	/ ( facade_element_area &
1106	- window_area_per_facade )
1107	!
1108	!-- Transfer q_wall_win back to USM (innermost wall/window layer)
1109	surf_usm_h%iwghf_eb(m) = q_wall_win
1110	surf_usm_h%iwghf_eb_window(m) = q_wall_win
1111	!
1112	!-- Sum up operational indoor temperature per kk-level. Further below,
1113	!-- this temperature is reduced by MPI to one temperature per kk-level
1114	!-- and building (processor overlapping)
1115	buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
1116	!
1117	!-- Calculation of waste heat
1118	!-- Anthropogenic heat output
1119	IF ( phi_hc_nd_ac > 0 ) THEN
1120	heating_on = 1
1121	cooling_on = 0
1122	ELSE
1123	heating_on = 0
1124	cooling_on = 1
1125	ENDIF
1126
1127	q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on)) !< [W/m2] anthropogenic heat output
1128	! surf_usm_h%shf(m)=q_waste_heat
1129
1130	ENDDO !< Horizontal surfaces loop
1131	!
1132	!-- Vertical surfaces
1133	DO fa = 1, buildings(nb)%num_facades_per_building_v_l
1134	!
1135	!-- Determine indices where corresponding surface-type information
1136	!-- is stored.
1137	l = buildings(nb)%l_v(fa)
1138	m = buildings(nb)%m_v(fa)
1139	!
1140	!-- Determine building height level index.
1141	kk = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
1142	!
1143	!-- Building geometries --> not time-dependent
1144	IF ( l == 0 .OR. l == 1 ) facade_element_area = dx * dzw(kk) !< [m2] surface area per facade element
1145	IF ( l == 2 .OR. l == 3 ) facade_element_area = dy * dzw(kk) !< [m2] surface area per facade element
1146	floor_area_per_facade = buildings(nb)%vpf(kk) * ddzw(kk) !< [m2] net floor area per facade element
1147	indoor_volume_per_facade = buildings(nb)%vpf(kk) !< [m3] indoor air volume per facade element
1148	window_area_per_facade = surf_usm_v(l)%frac(ind_wat_win,m) * facade_element_area !< [m2] window area per facade element
1149	eff_mass_area = factor_a * floor_area_per_facade !< [m2] standard values according to Table 12 section 12.3.1.2 (calculate over Eq. (65) according to section 12.3.1.2)
1150	c_m = factor_c * floor_area_per_facade !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
1151	total_area = lambda_at * floor_area_per_facade !< [m2] area of all surfaces pointing to zone Eq. (9) according to section 7.2.2.2
1152	!
1153	!-- Calculation of heat transfer coefficient for transmission --> not time-dependent
1154	h_tr_w = window_area_per_facade * u_value_win !< [W/K] only for windows
1155	h_tr_is = total_area * h_is !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
1156	h_tr_ms = eff_mass_area * h_ms !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
1157	h_tr_op = 1 / ( 1 / ( ( facade_element_area - window_area_per_facade ) &
1158	* lambda_layer3 / s_layer3 * 0.5 ) + 1 / h_tr_ms )
1159	h_tr_em = 1 / ( 1 / h_tr_op - 1 / h_tr_ms ) !< [W/K] Eq. (63), Section 12.2.2
1160	!
1161	!-- internal air loads dependent on the occupacy of the room
1162	!-- basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
1163	phi_ia = 0.5 * ( ( qint_high * schedule_d + qint_low ) &
1164	* floor_area_per_facade ) !< [W] Eq. (C.1)
1165	!
1166	!-- Airflow dependent on the occupacy of the room
1167	!-- basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
1168	air_change = ( air_change_high * schedule_d + air_change_low )
1169	!
1170	!-- Heat transfer of ventilation
1171	!-- not less than 0.01 W/K to provide division by 0 in further calculations
1172	!-- with heat capacity of air 0.33 Wh/m2K
1173	h_ve = MAX( 0.01 , ( air_change * indoor_volume_per_facade * &
1174	0.33 * (1 - eta_ve ) ) ) !< [W/K] from ISO 13789 Eq.(10)
1175
1176	!-- Heat transfer coefficient auxiliary variables
1177	h_tr_1 = 1 / ( ( 1 / h_ve ) + ( 1 / h_tr_is ) ) !< [W/K] Eq. (C.6)
1178	h_tr_2 = h_tr_1 + h_tr_w !< [W/K] Eq. (C.7)
1179	h_tr_3 = 1 / ( ( 1 / h_tr_2 ) + ( 1 / h_tr_ms ) ) !< [W/K] Eq. (C.8)
1180	!
1181	!-- Net short-wave radiation through window area (was i_global)
1182	net_sw_in = surf_usm_v(l)%rad_sw_in(m) - surf_usm_v(l)%rad_sw_out(m)
1183	!
1184	!-- Quantities needed for im_calc_temperatures
1185	i = surf_usm_v(l)%i(m)
1186	j = surf_usm_v(l)%j(m)
1187	k = surf_usm_v(l)%k(m)
1188	near_facade_temperature = t_surf_10cm_v(l)%t(m)
1189	indoor_wall_window_temperature = &
1190	surf_usm_v(l)%frac(ind_veg_wall,m) * t_wall_v(l)%t(nzt_wall,m) &
1191	+ surf_usm_v(l)%frac(ind_wat_win,m) * t_window_v(l)%t(nzt_wall,m)
1192	!
1193	!-- Solar thermal gains. If net_sw_in larger than sun-protection
1194	!-- threshold parameter (params_solar_protection), sun protection will
1195	!-- be activated
1196	IF ( net_sw_in <= params_solar_protection ) THEN
1197	solar_protection_off = 1
1198	solar_protection_on = 0
1199	ELSE
1200	solar_protection_off = 0
1201	solar_protection_on = 1
1202	ENDIF
1203	!
1204	!-- Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
1205	!-- DIN 4108 - 2 chap.8
1206	phi_sol = ( window_area_per_facade * net_sw_in * solar_protection_off &
1207	+ window_area_per_facade * net_sw_in * f_c_win * solar_protection_on ) &
1208	* g_value_win * ( 1 - params_f_f ) * params_f_w
1209	!
1210	!-- Calculation of the mass specific thermal load for internal and external heatsources
1211	phi_m = (eff_mass_area / total_area) * ( phi_ia + phi_sol ) !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
1212	!
1213	!-- Calculation mass specific thermal load implied non thermal mass
1214	phi_st = ( 1 - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1 * total_area ) ) ) &
1215	* ( phi_ia + phi_sol ) !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
1216	!
1217	!-- Calculations for deriving indoor temperature and heat flux into the wall
1218	!-- Step 1: Indoor temperature without heating and cooling
1219	!-- section C.4.1 Picture C.2 zone 3)
1220	phi_hc_nd = 0
1221
1222	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1223	near_facade_temperature, phi_hc_nd )
1224	!
1225	!-- If air temperature between border temperatures of heating and cooling, assign output variable, then ready
1226	IF ( theta_int_h_set <= theta_air .AND. theta_air <= theta_int_c_set ) THEN
1227	phi_hc_nd_ac = 0
1228	phi_hc_nd = phi_hc_nd_ac
1229	theta_air_ac = theta_air
1230	!
1231	!-- Step 2: Else, apply 10 W/mÂ² heating/cooling power and calculate indoor temperature
1232	!-- again.
1233	ELSE
1234	!
1235	!-- Temperature not correct, calculation method according to section C4.2
1236	theta_air_0 = theta_air !< Note temperature without heating/cooling
1237
1238	!-- Heating or cooling?
1239	IF ( theta_air > theta_int_c_set ) THEN
1240	theta_air_set = theta_int_c_set
1241	ELSE
1242	theta_air_set = theta_int_h_set
1243	ENDIF
1244
1245	!-- Calculate the temperature with phi_hc_nd_10
1246	phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
1247	phi_hc_nd = phi_hc_nd_10
1248
1249	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1250	near_facade_temperature, phi_hc_nd )
1251
1252	theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
1253
1254
1255	phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0) &
1256	/ (theta_air_10 - theta_air_0) !< Eq. (C.13)
1257	!
1258	!-- Step 3: With temperature ratio to determine the heating or cooling capacity
1259	!-- If necessary, limit the power to maximum power
1260	!-- section C.4.1 Picture C.2 zone 2) and 4)
1261	IF ( phi_c_max < phi_hc_nd_un .AND. phi_hc_nd_un < phi_h_max ) THEN
1262	phi_hc_nd_ac = phi_hc_nd_un
1263	phi_hc_nd = phi_hc_nd_un
1264	ELSE
1265	!-- Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
1266	!-- section C.4.1 Picture C.2 zone 1) and 5)
1267	IF ( phi_hc_nd_un > 0 ) THEN
1268	phi_hc_nd_ac = phi_h_max !< Limit heating
1269	ELSE
1270	phi_hc_nd_ac = phi_c_max !< Limit cooling
1271	ENDIF
1272	ENDIF
1273
1274	phi_hc_nd = phi_hc_nd_ac
1275	!
1276	!-- Calculate the temperature with phi_hc_nd_ac (new)
1277	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1278	near_facade_temperature, phi_hc_nd )
1279
1280	theta_air_ac = theta_air
1281
1282	ENDIF
1283	!
1284	!-- Update theta_m_t_prev
1285	theta_m_t_prev = theta_m_t
1286	!
1287	!-- Calculate the operating temperature with weighted mean of temperature of air and mean
1288	!-- Will be used for thermal comfort calculations
1289	theta_op = 0.3 * theta_air_ac + 0.7 * theta_s
1290	!
1291	!-- Heat flux into the wall. Value needed in urban_surface_mod to
1292	!-- calculate heat transfer through wall layers towards the facade
1293	q_wall_win = h_tr_ms * ( theta_s - theta_m ) &
1294	/ ( facade_element_area &
1295	- window_area_per_facade )
1296	!
1297	!-- Transfer q_wall_win back to USM (innermost wall/window layer)
1298	surf_usm_v(l)%iwghf_eb(m) = q_wall_win
1299	surf_usm_v(l)%iwghf_eb_window(m) = q_wall_win
1300	!
1301	!-- Sum up operational indoor temperature per kk-level. Further below,
1302	!-- this temperature is reduced by MPI to one temperature per kk-level
1303	!-- and building (processor overlapping)
1304	buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
1305
1306	!
1307	!-- Calculation of waste heat
1308	!-- Anthropogenic heat output
1309	IF ( phi_hc_nd_ac > 0 ) THEN
1310	heating_on = 1
1311	cooling_on = 0
1312	ELSE
1313	heating_on = 0
1314	cooling_on = 1
1315	ENDIF
1316
1317	q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on)) !< [W/m2] , anthropogenic heat output
1318	! surf_usm_v(l)%waste_heat(m)=q_waste_heat
1319
1320	ENDDO !< Vertical surfaces loop
1321
1322	ENDIF !< buildings(nb)%on_pe
1323	ENDDO !< buildings loop
1324
1325	!
1326	!-- Determine total number of facade elements per building and assign number to
1327	!-- building data type.
1328	DO nb = 1, num_build
1329	!
1330	!-- Allocate dummy array used for summing-up facade elements.
1331	!-- Please note, dummy arguments are necessary as building-date type
1332	!-- arrays are not necessarily allocated on all PEs.
1333	ALLOCATE( t_in_l_send(buildings(nb)%kb_min:buildings(nb)%kb_max) )
1334	ALLOCATE( t_in_recv(buildings(nb)%kb_min:buildings(nb)%kb_max) )
1335	t_in_l_send = 0.0_wp
1336	t_in_recv = 0.0_wp
1337
1338	IF ( buildings(nb)%on_pe ) THEN
1339	t_in_l_send = buildings(nb)%t_in_l
1340	ENDIF
1341
1342	#if defined( __parallel )
1343	CALL MPI_ALLREDUCE( t_in_l_send, &
1344	t_in_recv, &
1345	buildings(nb)%kb_max - buildings(nb)%kb_min + 1, &
1346	MPI_REAL, &
1347	MPI_SUM, &
1348	comm2d, &
1349	ierr )
1350
1351	IF ( ALLOCATED( buildings(nb)%t_in ) ) &
1352	buildings(nb)%t_in = t_in_recv
1353	#else
1354	buildings(nb)%t_in = buildings(nb)%t_in_l
1355	#endif
1356
1357	buildings(nb)%t_in = buildings(nb)%t_in / &
1358	( buildings(nb)%num_facade_h + &
1359	buildings(nb)%num_facade_v )
1360	!
1361	!-- Deallocate dummy arrays
1362	DEALLOCATE( t_in_l_send )
1363	DEALLOCATE( t_in_recv )
1364
1365	ENDDO
1366
1367
1368	END SUBROUTINE im_main_heatcool
1369
1370	!------------------------------------------------------------------------------!
1371	! Description:
1372	! ------------
1373	!> Parin for &indoor_parameters for indoor model
1374	!------------------------------------------------------------------------------!
1375	SUBROUTINE im_parin
1376
1377	USE control_parameters, &
1378	ONLY: indoor_model
1379
1380	IMPLICIT NONE
1381
1382	CHARACTER (LEN=80) :: line !< string containing current line of file PARIN
1383
1384
1385
1386	NAMELIST /indoor_parameters/ building_type, dt_indoor, &
1387	initial_indoor_temperature
1388
1389	! line = ' '
1390
1391	!
1392	!-- Try to find indoor model package
1393	REWIND ( 11 )
1394	line = ' '
1395	DO WHILE ( INDEX( line, '&indoor_parameters' ) == 0 )
1396	READ ( 11, '(A)', END=10 ) line
1397	! PRINT*, 'line: ', line
1398	ENDDO
1399	BACKSPACE ( 11 )
1400
1401	!
1402	!-- Read user-defined namelist
1403	READ ( 11, indoor_parameters )
1404	!
1405	!-- Set flag that indicates that the indoor model is switched on
1406	indoor_model = .TRUE.
1407
1408	!
1409	!-- Activate spinup (maybe later
1410	! IF ( spinup_time > 0.0_wp ) THEN
1411	! coupling_start_time = spinup_time
1412	! end_time = end_time + spinup_time
1413	! IF ( spinup_pt_mean == 9999999.9_wp ) THEN
1414	! spinup_pt_mean = pt_surface
1415	! ENDIF
1416	! spinup = .TRUE.
1417	! ENDIF
1418
1419	10 CONTINUE
1420
1421	END SUBROUTINE im_parin
1422
1423
1424	END MODULE indoor_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |