Home

Context Navigation

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

Last change on this file since 3775 was 3759, checked in by suehring, 6 years ago
Indoor model: provide total building volume; several bugfixes; calculation of building height in IM revised
Property svn:keywords set to `Id`
File size: 85.9 KB

Line
1	!> @file indoor_model_mod.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
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 3759 2019-02-21 15:53:45Z gronemeier $
28	! - Calculation of total building volume
29	! - Several bugfixes
30	! - Calculation of building height revised
31	!
32	! 3745 2019-02-15 18:57:56Z suehring
33	! - remove building_type from module
34	! - initialize parameters for each building individually instead of a bulk
35	! initializaion with identical building type for all
36	! - output revised
37	! - add missing _wp
38	! - some restructuring of variables in building data structure
39	!
40	! 3744 2019-02-15 18:38:58Z suehring
41	! Some interface calls moved to module_interface + cleanup
42	!
43	! 3597 2018-12-04 08:40:18Z maronga
44	! Renamed t_surf_10cm to pt_10cm
45	!
46	! 3593 2018-12-03 13:51:13Z kanani
47	! Replace degree symbol by degree_C
48	!
49	! 3524 2018-11-14 13:36:44Z raasch
50	! working precision added to make code Fortran 2008 conform
51	!
52	! 3469 2018-10-30 20:05:07Z kanani
53	! Initial revision (tlang, suehring, kanani, srissman)
54	!
55	!
56	!
57	! Authors:
58	! --------
59	! @author Tobias Lang
60	! @author Jens Pfafferott
61	! @author Farah Kanani-Suehring
62	! @author Matthias Suehring
63	! @author Sascha RiÃmann
64	!
65	!
66	! Description:
67	! ------------
68	!> <Description of the new module>
69	!> Module for Indoor Climate Model (ICM)
70	!> The module is based on the DIN EN ISO 13790 with simplified hour-based procedure.
71	!> This model is a equivalent circuit diagram of a three-point RC-model (5R1C).
72	!> This module differ between indoor-air temperature an average temperature of indoor surfaces which make it prossible to determine thermal comfort
73	!> the heat transfer between indoor and outdoor is simplified
74
75	!> @todo Replace window_area_per_facade by %frac(1,m) for window
76	!> @todo emissivity change for window blinds if solar_protection_on=1
77	!> @todo write datas in netcdf file as output data
78	!> @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
79	!>
80	!> @note Do we allow use of integer flags, or only logical flags? (concerns e.g. cooling_on, heating_on)
81	!> @note How to write indoor temperature output to pt array?
82	!>
83	!> @bug <Enter known bugs here>
84	!------------------------------------------------------------------------------!
85	MODULE indoor_model_mod
86
87	USE control_parameters, &
88	ONLY: initializing_actions
89
90	USE kinds
91
92	USE netcdf_data_input_mod, &
93	ONLY: building_id_f, building_type_f
94
95	USE surface_mod, &
96	ONLY: surf_usm_h, surf_usm_v
97
98
99	IMPLICIT NONE
100
101	!
102	!-- Define data structure for buidlings.
103	TYPE build
104
105	INTEGER(iwp) :: id !< building ID
106	INTEGER(iwp) :: kb_min !< lowest vertical index of a building
107	INTEGER(iwp) :: kb_max !< highest vertical index of a building
108	INTEGER(iwp) :: num_facades_per_building_h = 0 !< total number of horizontal facades elements
109	INTEGER(iwp) :: num_facades_per_building_h_l = 0 !< number of horizontal facade elements on local subdomain
110	INTEGER(iwp) :: num_facades_per_building_v = 0 !< total number of vertical facades elements
111	INTEGER(iwp) :: num_facades_per_building_v_l = 0 !< number of vertical facade elements on local subdomain
112
113	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: l_v !< index array linking surface-element orientation index
114	!< for vertical surfaces with building
115	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: m_h !< index array linking surface-element index for
116	!< horizontal surfaces with building
117	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: m_v !< index array linking surface-element index for
118	!< vertical surfaces with building
119	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facade_h !< number of horizontal facade elements per buidling
120	!< and height level
121	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facade_v !< number of vertical facades elements per buidling
122	!< and height level
123
124	INTEGER(iwp) :: ventilation_int_loads
125
126	LOGICAL :: on_pe = .FALSE. !< flag indicating whether a building with certain ID is on local subdomain
127
128	REAL(wp) :: building_height !< building height
129	REAL(wp) :: lambda_layer3 !< [W/(m*K)] Thermal conductivity of the inner layer
130	REAL(wp) :: s_layer3 !< [m] half thickness of the inner layer (layer_3)
131	REAL(wp) :: f_c_win !< [-] shading factor
132	REAL(wp) :: g_value_win !< [-] SHGC factor
133	REAL(wp) :: u_value_win !< [W/(m2*K)] transmittance
134	REAL(wp) :: air_change_low !< [1/h] air changes per time_utc_hour
135	REAL(wp) :: air_change_high !< [1/h] air changes per time_utc_hour
136	REAL(wp) :: eta_ve !< [-] heat recovery efficiency
137	REAL(wp) :: factor_a !< [-] Dynamic parameters specific effective surface according to Table 12; 2.5
138	!< (very light, light and medium), 3.0 (heavy), 3.5 (very heavy)
139	REAL(wp) :: factor_c !< [J/(m2 K)] Dynamic parameters inner heatstorage according to Table 12; 80000
140	!< (very light), 110000 (light), 165000 (medium), 260000 (heavy), 370000 (very heavy)
141	REAL(wp) :: lambda_at !< [-] ratio internal surface/floor area chap. 7.2.2.2.
142	REAL(wp) :: theta_int_h_set !< [degree_C] Max. Setpoint temperature (winter)
143	REAL(wp) :: theta_int_c_set !< [degree_C] Max. Setpoint temperature (summer)
144	REAL(wp) :: phi_h_max !< [W] Max. Heating capacity (negative)
145	REAL(wp) :: phi_c_max !< [W] Max. Cooling capacity (negative)
146	REAL(wp) :: qint_high !< [W/m2] internal heat gains, option Database qint_0-23
147	REAL(wp) :: qint_low !< [W/m2] internal heat gains, option Database qint_0-23
148	REAL(wp) :: height_storey !< [m] storey heigth
149	REAL(wp) :: height_cei_con !< [m] ceiling construction heigth
150	REAL(wp) :: vol_tot !< total building volume
151
152	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in !< mean building indoor temperature, height dependent
153	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_l !< mean building indoor temperature on local subdomain, height dependent
154	REAL(wp), DIMENSION(:), ALLOCATABLE :: volume !< total building volume, height dependent
155	REAL(wp), DIMENSION(:), ALLOCATABLE :: vol_frac !< fraction of local on total building volume, height dependent
156	REAL(wp), DIMENSION(:), ALLOCATABLE :: vpf !< building volume volume per facade element, height dependent
157
158	END TYPE build
159
160	TYPE(build), DIMENSION(:), ALLOCATABLE :: buildings !< building array
161
162	INTEGER(iwp) :: num_build !< total number of buildings in domain
163
164	REAL(wp) :: volume_fraction
165
166	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in !< dummy array for indoor temperature for the
167	!< total building volume at each discrete height level
168	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_l !< dummy array for indoor temperature for the
169	!< local building volume fraction at each discrete height level
170
171	!
172	!-- Declare all global variables within the module
173	INTEGER(iwp) :: cooling_on !< Indoor cooling flag (0=off, 1=on)
174	INTEGER(iwp) :: heating_on !< Indoor heating flag (0=off, 1=on)
175	INTEGER(iwp) :: solar_protection_off !< Solar protection off
176	INTEGER(iwp) :: solar_protection_on !< Solar protection on
177
178	REAL(wp) :: eff_mass_area !< [m2] the effective mass-related area
179	REAL(wp) :: floor_area_per_facade !< [m2] net floor area (Sum of all floors)
180	REAL(wp) :: total_area !<! [m2] area of all surfaces pointing to zone
181	REAL(wp) :: window_area_per_facade !< [m2] window area per facade element
182	REAL(wp) :: air_change !< [1/h] Airflow
183	REAL(wp) :: facade_element_area !< [m2_facade] building surface facade
184	REAL(wp) :: indoor_volume_per_facade !< [m3] indoor air volume per facade element
185	REAL(wp) :: c_m !< [J/K] internal heat storage capacity
186	REAL(wp) :: dt_indoor = 3600.0_wp !< [s] namelist parameter: time interval for indoor-model application
187	REAL(wp) :: h_tr_1 !<! [W/K] Heat transfer coefficient auxiliary variable 1
188	REAL(wp) :: h_tr_2 !<! [W/K] Heat transfer coefficient auxiliary variable 2
189	REAL(wp) :: h_tr_3 !<! [W/K] Heat transfer coefficient auxiliary variable 3
190	REAL(wp) :: h_tr_em !<! [W/K] Heat transfer coefficient of the emmision (got with h_tr_ms the thermal mass)
191	REAL(wp) :: h_tr_is !<! [W/K] thermal coupling conductance (Thermischer Kopplungsleitwert)
192	REAL(wp) :: h_tr_ms !<! [W/K] Heat transfer conductance term (got with h_tr_em the thermal mass)
193	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
194	REAL(wp) :: h_tr_w !<! [W/K] heat transfer coefficient of doors, windows, curtain walls and glazed walls (assumption: thermal mass=0)
195	REAL(wp) :: h_ve !<! [W/K] heat transfer of ventilation
196	REAL(wp) :: initial_indoor_temperature !< namelist parameter
197	REAL(wp) :: net_sw_in !< net short-wave radiation (in - out; was i_global --> CORRECT?)
198	REAL(wp) :: phi_hc_nd !<! [W] heating demand and/or cooling demand
199	REAL(wp) :: phi_hc_nd_10 !<! [W] heating demand and/or cooling demand for heating or cooling
200	REAL(wp) :: phi_hc_nd_ac !<! [W] actual heating demand and/or cooling demand
201	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)
202	REAL(wp) :: phi_ia !< [W] internal air load = internal loads * 0.5, Eq. (C.1)
203	REAL(wp) :: phi_m !<! [W] mass specific thermal load (internal and external)
204	REAL(wp) :: phi_mtot !<! [W] total mass specific thermal load (internal and external)
205	REAL(wp) :: phi_sol !< [W] solar loads
206	REAL(wp) :: phi_st !<! [W] mass specific thermal load implied non thermal mass
207	REAL(wp) :: q_emission !< emissions, in first version = 0, option for second part of the project
208	REAL(wp) :: q_wall_win !< heat flux from indoor into wall/window
209	REAL(wp) :: q_waste_heat !< waste heat, sum of waste heat over the roof to Palm
210	REAL(wp) :: q_waste_heat_bldg !< [W/building] waste heat of the complete building, in Palm sum of all indoor_model-calculations
211	REAL(wp) :: schedule_d !< activation for internal loads (low or high + low)
212	REAL(wp) :: skip_time_do_indoor = 0.0_wp !< [s] Indoor model is not called before this time
213	REAL(wp) :: theta_air !<! [degree_C] air temperature of the RC-node
214	REAL(wp) :: theta_air_0 !<! [degree_C] air temperature of the RC-node in equilibrium
215	REAL(wp) :: theta_air_10 !<! [degree_C] air temperature of the RC-node from a heating capacity of 10 W/m2
216	REAL(wp) :: theta_air_ac !< [degree_C] actual room temperature after heating/cooling
217	REAL(wp) :: theta_air_set !< [degree_C] Setpoint_temperature for the room
218	REAL(wp) :: theta_m !<! [degree_C} inner temperature of the RC-node
219	REAL(wp) :: theta_m_t !<! [degree_C] (Fictive) component temperature timestep
220	REAL(wp) :: theta_m_t_prev !< [degree_C] (Fictive) component temperature previous timestep (do not change)
221	REAL(wp) :: theta_op !< [degree_C] operative temperature
222	REAL(wp) :: theta_s !<! [degree_C] surface temperature of the RC-node
223	REAL(wp) :: time_indoor = 0.0_wp !< [s] time since last call of indoor model
224	REAL(wp) :: time_utc_hour !< Time in hours per day (UTC)
225	REAL(wp) :: ventilation_int_loads !< Zuteilung der GebÃ€ude fÃŒr Verlauf/AktivitÃ€t der LÃŒftung und internen Lasten
226
227	REAL(wp) :: f_sr !< [-] factor surface reduction
228	REAL(wp) :: f_cei !< [-] ceiling reduction factor
229	REAL(wp) :: ngs !< [m2] netto ground surface
230
231	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
232	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)
233	REAL(wp), PARAMETER :: params_f_win = 0.5_wp !< [-] proportion of window area, Database A_win aus Datenbank 27 window_area_per_facade_percent
234	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
235	REAL(wp), PARAMETER :: params_waste_heat_c = 4.0_wp !< [-] anthropogenic heat outputs for cooling e.g. 4 for KKM with COP = 3
236	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
237	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)
238	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)
239
240
241
242	SAVE
243
244
245	PRIVATE
246
247	!
248	!-- Add INTERFACES that must be available to other modules
249	PUBLIC im_init, im_main_heatcool, im_parin, im_define_netcdf_grid, &
250	im_check_data_output, im_data_output_3d, im_check_parameters
251
252
253	!
254	!-- Add VARIABLES that must be available to other modules
255	PUBLIC dt_indoor, skip_time_do_indoor, time_indoor
256
257	!
258	!-- PALM interfaces:
259	!-- Data output checks for 2D/3D data to be done in check_parameters
260	INTERFACE im_check_data_output
261	MODULE PROCEDURE im_check_data_output
262	END INTERFACE im_check_data_output
263	!
264	!-- Input parameter checks to be done in check_parameters
265	INTERFACE im_check_parameters
266	MODULE PROCEDURE im_check_parameters
267	END INTERFACE im_check_parameters
268	!
269	!-- Data output of 3D data
270	INTERFACE im_data_output_3d
271	MODULE PROCEDURE im_data_output_3d
272	END INTERFACE im_data_output_3d
273
274	!
275	!-- Definition of data output quantities
276	INTERFACE im_define_netcdf_grid
277	MODULE PROCEDURE im_define_netcdf_grid
278	END INTERFACE im_define_netcdf_grid
279	!
280	! !
281	! !-- Output of information to the header file
282	! INTERFACE im_header
283	! MODULE PROCEDURE im_header
284	! END INTERFACE im_header
285	!
286	!-- Calculations for indoor temperatures
287	INTERFACE im_calc_temperatures
288	MODULE PROCEDURE im_calc_temperatures
289	END INTERFACE im_calc_temperatures
290	!
291	!-- Initialization actions
292	INTERFACE im_init
293	MODULE PROCEDURE im_init
294	END INTERFACE im_init
295	!
296	!-- Main part of indoor model
297	INTERFACE im_main_heatcool
298	MODULE PROCEDURE im_main_heatcool
299	END INTERFACE im_main_heatcool
300	!
301	!-- Reading of NAMELIST parameters
302	INTERFACE im_parin
303	MODULE PROCEDURE im_parin
304	END INTERFACE im_parin
305
306	CONTAINS
307
308	!------------------------------------------------------------------------------!
309	! Description:
310	! ------------
311	!< Calculation of the air temperatures and mean radiation temperature
312	!< This is basis for the operative temperature
313	!< Based on a Crank-Nicholson scheme with a timestep of a hour
314	!------------------------------------------------------------------------------!
315	SUBROUTINE im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
316	near_facade_temperature, phi_hc_nd_dummy )
317
318	USE arrays_3d, &
319	ONLY: pt
320
321
322	IMPLICIT NONE
323
324
325	INTEGER(iwp) :: i
326	INTEGER(iwp) :: j
327	INTEGER(iwp) :: k
328
329	REAL(wp) :: indoor_wall_window_temperature !< weighted temperature of innermost wall/window layer
330	REAL(wp) :: near_facade_temperature
331	REAL(wp) :: phi_hc_nd_dummy
332
333
334	!< Calculation of total mass specific thermal load (internal and external)
335	phi_mtot = ( phi_m + h_tr_em * indoor_wall_window_temperature &
336	+ h_tr_3 * ( phi_st + h_tr_w * pt(k,j,i) &
337	+ h_tr_1 * &
338	( ( ( phi_ia + phi_hc_nd_dummy ) / h_ve ) &
339	+ near_facade_temperature ) &
340	) / h_tr_2 &
341	) !< [degree_C] Eq. (C.5)
342
343	!< Calculation of component temperature at factual timestep
344	theta_m_t = ( ( theta_m_t_prev &
345	* ( ( c_m / 3600.0_wp ) - 0.5_wp * ( h_tr_3 + h_tr_em ) ) + phi_mtot &
346	) &
347	/ ( ( c_m / 3600.0_wp ) + 0.5_wp * ( h_tr_3 + h_tr_em ) ) &
348	) !< [degree_C] Eq. (C.4)
349
350	!< Calculation of mean inner temperature for the RC-node in actual timestep
351	theta_m = ( theta_m_t + theta_m_t_prev ) * 0.5_wp !< [degree_C] Eq. (C.9)
352
353	!< Calculation of mean surface temperature of the RC-node in actual timestep
354	theta_s = ( ( h_tr_ms * theta_m + phi_st + h_tr_w * pt(k,j,i) &
355	+ h_tr_1 * ( near_facade_temperature + ( phi_ia + phi_hc_nd_dummy ) / h_ve ) &
356	) &
357	/ ( h_tr_ms + h_tr_w + h_tr_1 ) &
358	) !< [degree_C] Eq. (C.10)
359
360	!< Calculation of the air temperature of the RC-node
361	theta_air = ( h_tr_is * theta_s + h_ve * near_facade_temperature &
362	+ phi_ia + phi_hc_nd_dummy ) / ( h_tr_is + h_ve ) !< [degree_C] Eq. (C.11)
363
364	END SUBROUTINE im_calc_temperatures
365
366	!------------------------------------------------------------------------------!
367	! Description:
368	! ------------
369	!> Initialization of the indoor model.
370	!> Static information are calculated here, e.g. building parameters and
371	!> geometrical information, everything that doesn't change in time.
372	!
373	!-- Input values
374	!-- Input datas from Palm, M4
375	! i_global --> net_sw_in !global radiation [W/m2]
376	! theta_e --> pt(k,j,i) !undisturbed outside temperature, 1. PALM volume, for windows
377	! theta_sup = theta_f --> surf_usm_h%pt_10cm(m)
378	! surf_usm_v(l)%pt_10cm(m) !Air temperature, facade near (10cm) air temperature from 1. Palm volume
379	! theta_node --> t_wall_h(nzt_wall,m)
380	! t_wall_v(l)%t(nzt_wall,m) !Temperature of innermost wall layer, for opaque wall
381	!------------------------------------------------------------------------------!
382	SUBROUTINE im_init
383
384	USE arrays_3d, &
385	ONLY: dzw
386
387	USE control_parameters, &
388	ONLY: message_string
389
390	USE indices, &
391	ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
392
393	USE grid_variables, &
394	ONLY: dx, dy
395
396	USE pegrid
397
398	USE surface_mod, &
399	ONLY: surf_usm_h, surf_usm_v
400
401	USE urban_surface_mod, &
402	ONLY: building_pars, building_type
403
404	IMPLICIT NONE
405
406	INTEGER(iwp) :: bt !< local building type
407	INTEGER(iwp) :: fa !< running index for facade elements of each building
408	INTEGER(iwp) :: i !< running index along x-direction
409	INTEGER(iwp) :: j !< running index along y-direction
410	INTEGER(iwp) :: k !< running index along z-direction
411	INTEGER(iwp) :: l !< running index for surface-element orientation
412	INTEGER(iwp) :: m !< running index surface elements
413	INTEGER(iwp) :: n !< building index
414	INTEGER(iwp) :: nb !< building index
415
416	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids !< building IDs on entire model domain
417	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final !< building IDs on entire model domain,
418	!< multiple occurences are sorted out
419	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_final_tmp !< temporary array used for resizing
420	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l !< building IDs on local subdomain
421	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: build_ids_l_tmp !< temporary array used to resize array of building IDs
422	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displace_dum !< displacements of start addresses, used for MPI_ALLGATHERV
423	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_max_l !< highest vertical index of a building on subdomain
424	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: k_min_l !< lowest vertical index of a building on subdomain
425	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: n_fa !< counting array
426	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facades_h !< dummy array used for summing-up total number of
427	!< horizontal facade elements
428	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: num_facades_v !< dummy array used for summing-up total number of
429	!< vertical facade elements
430	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: receive_dum_h !< dummy array used for MPI_ALLREDUCE
431	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: receive_dum_v !< dummy array used for MPI_ALLREDUCE
432
433	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings !< number of buildings with different ID on entire model domain
434	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_buildings_l !< number of buildings with different ID on local subdomain
435
436	REAL(wp), DIMENSION(:), ALLOCATABLE :: local_weight !< dummy representing fraction of local on total building volume,
437	!< height dependent
438	REAL(wp), DIMENSION(:), ALLOCATABLE :: volume !< total building volume at each discrete height level
439	REAL(wp), DIMENSION(:), ALLOCATABLE :: volume_l !< total building volume at each discrete height level,
440	!< on local subdomain
441
442	CALL location_message( 'initializing indoor model', .FALSE. )
443
444	!
445	!-- Initializing of indoor model is only possible if buildings can be
446	!-- distinguished by their IDs.
447	IF ( .NOT. building_id_f%from_file ) THEN
448	message_string = 'Indoor model requires information about building_id'
449	CALL message( 'im_init', 'PA0999', 1, 2, 0, 6, 0 )
450	ENDIF
451	!
452	!-- Determine number of different building IDs on local subdomain.
453	num_buildings_l = 0
454	num_buildings = 0
455	ALLOCATE( build_ids_l(1) )
456	DO i = nxl, nxr
457	DO j = nys, nyn
458	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
459	IF ( num_buildings_l(myid) > 0 ) THEN
460	IF ( ANY( building_id_f%var(j,i) .EQ. build_ids_l ) ) THEN
461	CYCLE
462	ELSE
463	num_buildings_l(myid) = num_buildings_l(myid) + 1
464	!
465	!-- Resize array with different local building ids
466	ALLOCATE( build_ids_l_tmp(1:SIZE(build_ids_l)) )
467	build_ids_l_tmp = build_ids_l
468	DEALLOCATE( build_ids_l )
469	ALLOCATE( build_ids_l(1:num_buildings_l(myid)) )
470	build_ids_l(1:num_buildings_l(myid)-1) = &
471	build_ids_l_tmp(1:num_buildings_l(myid)-1)
472	build_ids_l(num_buildings_l(myid)) = building_id_f%var(j,i)
473	DEALLOCATE( build_ids_l_tmp )
474	ENDIF
475	!
476	!-- First occuring building id on PE
477	ELSE
478	num_buildings_l(myid) = num_buildings_l(myid) + 1
479	build_ids_l(1) = building_id_f%var(j,i)
480	ENDIF
481	ENDIF
482	ENDDO
483	ENDDO
484	!
485	!-- Determine number of building IDs for the entire domain. (Note, building IDs
486	!-- can appear multiple times as buildings might be distributed over several
487	!-- PEs.)
488	#if defined( __parallel )
489	CALL MPI_ALLREDUCE( num_buildings_l, num_buildings, numprocs, &
490	MPI_INTEGER, MPI_SUM, comm2d, ierr )
491	#else
492	num_buildings = num_buildings_l
493	#endif
494	ALLOCATE( build_ids(1:SUM(num_buildings)) )
495	!
496	!-- Gather building IDs. Therefore, first, determine displacements used
497	!-- required for MPI_GATHERV call.
498	ALLOCATE( displace_dum(0:numprocs-1) )
499	displace_dum(0) = 0
500	DO i = 1, numprocs-1
501	displace_dum(i) = displace_dum(i-1) + num_buildings(i-1)
502	ENDDO
503
504	#if defined( __parallel )
505	CALL MPI_ALLGATHERV( build_ids_l(1:num_buildings_l(myid)), &
506	num_buildings(myid), &
507	MPI_INTEGER, &
508	build_ids, &
509	num_buildings, &
510	displace_dum, &
511	MPI_INTEGER, &
512	comm2d, ierr )
513
514	DEALLOCATE( displace_dum )
515
516	#else
517	build_ids = build_ids_l
518	#endif
519	!
520	!-- Note: in parallel mode, building IDs can occur mutliple times, as
521	!-- each PE has send its own ids. Therefore, sort out building IDs which
522	!-- appear multiple times.
523	num_build = 0
524	DO n = 1, SIZE(build_ids)
525
526	IF ( ALLOCATED(build_ids_final) ) THEN
527	IF ( ANY( build_ids(n) .EQ. build_ids_final ) ) THEN !FK: Warum ANY?, Warum .EQ.? --> s.o
528	CYCLE
529	ELSE
530	num_build = num_build + 1
531	!
532	!-- Resize
533	ALLOCATE( build_ids_final_tmp(1:num_build) )
534	build_ids_final_tmp(1:num_build-1) = build_ids_final(1:num_build-1)
535	DEALLOCATE( build_ids_final )
536	ALLOCATE( build_ids_final(1:num_build) )
537	build_ids_final(1:num_build-1) = build_ids_final_tmp(1:num_build-1)
538	build_ids_final(num_build) = build_ids(n)
539	DEALLOCATE( build_ids_final_tmp )
540	ENDIF
541	ELSE
542	num_build = num_build + 1
543	ALLOCATE( build_ids_final(1:num_build) )
544	build_ids_final(num_build) = build_ids(n)
545	ENDIF
546	ENDDO
547
548	!
549	!-- Allocate building-data structure array. Note, this is a global array
550	!-- and all building IDs on domain are known by each PE. Further attributes,
551	!-- e.g. height-dependent arrays, however, are only allocated on PEs where
552	!-- the respective building is present (in order to reduce memory demands).
553	ALLOCATE( buildings(1:num_build) )
554
555	!
556	!-- Store building IDs and check if building with certain ID is present on
557	!-- subdomain.
558	DO nb = 1, num_build
559	buildings(nb)%id = build_ids_final(nb)
560
561	IF ( ANY( building_id_f%var(nys:nyn,nxl:nxr) == buildings(nb)%id ) ) &
562	buildings(nb)%on_pe = .TRUE.
563	ENDDO
564	!
565	!-- Determine the maximum vertical dimension occupied by each building.
566	ALLOCATE( k_min_l(1:num_build) )
567	ALLOCATE( k_max_l(1:num_build) )
568	k_min_l = nzt + 1
569	k_max_l = 0
570
571	DO i = nxl, nxr
572	DO j = nys, nyn
573	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
574	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
575	DIM = 1 )
576	DO k = nzb+1, nzt+1
577	!
578	!-- Check if grid point belongs to a building.
579	IF ( BTEST( wall_flags_0(k,j,i), 6 ) ) THEN
580	k_min_l(nb) = MIN( k_min_l(nb), k )
581	k_max_l(nb) = MAX( k_max_l(nb), k )
582	ENDIF
583
584	ENDDO
585	ENDIF
586	ENDDO
587	ENDDO
588
589	DO nb = 1, num_build
590	#if defined( __parallel )
591	CALL MPI_ALLREDUCE( k_min_l(nb), buildings(nb)%kb_min, 1, MPI_INTEGER, &
592	MPI_MIN, comm2d, ierr )
593	CALL MPI_ALLREDUCE( k_max_l(nb), buildings(nb)%kb_max, 1, MPI_INTEGER, &
594	MPI_MAX, comm2d, ierr )
595	#else
596	buildings(nb)%kb_min = k_min_l(nb)
597	buildings(nb)%kb_max = k_max_l(nb)
598	#endif
599
600	ENDDO
601
602	DEALLOCATE( k_min_l )
603	DEALLOCATE( k_max_l )
604	!
605	!-- Calculate building height.
606	DO nb = 1, num_build
607	buildings(nb)%building_height = 0.0_wp
608	DO k = buildings(nb)%kb_min, buildings(nb)%kb_max
609	buildings(nb)%building_height = buildings(nb)%building_height &
610	+ dzw(k)
611	ENDDO
612	ENDDO
613	!
614	!-- Calculate building volume
615	DO nb = 1, num_build
616	!
617	!-- Allocate temporary array for summing-up building volume
618	ALLOCATE( volume(buildings(nb)%kb_min:buildings(nb)%kb_max) )
619	ALLOCATE( volume_l(buildings(nb)%kb_min:buildings(nb)%kb_max) )
620	volume = 0.0_wp
621	volume_l = 0.0_wp
622	!
623	!-- Calculate building volume per height level on each PE where
624	!-- these building is present.
625	IF ( buildings(nb)%on_pe ) THEN
626	ALLOCATE( buildings(nb)%volume(buildings(nb)%kb_min:buildings(nb)%kb_max) )
627	ALLOCATE( buildings(nb)%vol_frac(buildings(nb)%kb_min:buildings(nb)%kb_max) )
628	buildings(nb)%volume = 0.0_wp
629	buildings(nb)%vol_frac = 0.0_wp
630
631	IF ( ANY( building_id_f%var(nys:nyn,nxl:nxr) == buildings(nb)%id ) ) &
632	THEN
633	DO i = nxl, nxr
634	DO j = nys, nyn
635	DO k = buildings(nb)%kb_min, buildings(nb)%kb_max
636	IF ( building_id_f%var(j,i) /= building_id_f%fill ) &
637	volume_l(k) = volume_l(k) + dx * dy * dzw(k)
638	ENDDO
639	ENDDO
640	ENDDO
641	ENDIF
642	ENDIF
643	!
644	!-- Sum-up building volume from all subdomains
645	#if defined( __parallel )
646	CALL MPI_ALLREDUCE( volume_l, volume, SIZE(volume), MPI_REAL, MPI_SUM, &
647	comm2d, ierr )
648	#else
649	volume = volume_l
650	#endif
651	!
652	!-- Save total building volume as well as local fraction on volume on
653	!-- building data structure.
654	IF ( ALLOCATED( buildings(nb)%volume ) ) buildings(nb)%volume = volume
655	!
656	!-- Determine fraction of local on total building volume
657	IF ( buildings(nb)%on_pe ) buildings(nb)%vol_frac = volume_l / volume
658	!
659	!-- Calculate total building volume
660	IF ( ALLOCATED( buildings(nb)%volume ) ) &
661	buildings(nb)%vol_tot = SUM( buildings(nb)%volume )
662
663	DEALLOCATE( volume )
664	DEALLOCATE( volume_l )
665
666	ENDDO
667
668	!
669	!-- Allocate arrays for indoor temperature.
670	DO nb = 1, num_build
671	IF ( buildings(nb)%on_pe ) THEN
672	ALLOCATE( buildings(nb)%t_in(buildings(nb)%kb_min:buildings(nb)%kb_max) )
673	ALLOCATE( buildings(nb)%t_in_l(buildings(nb)%kb_min:buildings(nb)%kb_max) )
674	buildings(nb)%t_in = 0.0_wp
675	buildings(nb)%t_in_l = 0.0_wp
676	ENDIF
677	ENDDO
678	!
679	!-- Allocate arrays for number of facades per height level. Distinguish between
680	!-- horizontal and vertical facades.
681	DO nb = 1, num_build
682	IF ( buildings(nb)%on_pe ) THEN
683	ALLOCATE( buildings(nb)%num_facade_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
684	ALLOCATE( buildings(nb)%num_facade_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
685
686	buildings(nb)%num_facade_h = 0
687	buildings(nb)%num_facade_v = 0
688	ENDIF
689	ENDDO
690	!
691	!-- Determine number of facade elements per building on local subdomain.
692	!-- Distinguish between horizontal and vertical facade elements.
693	!
694	!-- Horizontal facades
695	buildings(:)%num_facades_per_building_h_l = 0
696	DO m = 1, surf_usm_h%ns
697	!
698	!-- For the current facade element determine corresponding building index.
699	!-- First, obtain j,j,k indices of the building. Please note the
700	!-- offset between facade/surface element and building location (for
701	!-- horizontal surface elements the horizontal offsets are zero).
702	i = surf_usm_h%i(m) + surf_usm_h%ioff
703	j = surf_usm_h%j(m) + surf_usm_h%joff
704	k = surf_usm_h%k(m) + surf_usm_h%koff
705	!
706	!-- Determine building index and check whether building is on PE
707	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
708	IF ( buildings(nb)%on_pe ) THEN
709	!
710	!-- Count number of facade elements at each height level.
711	buildings(nb)%num_facade_h(k) = buildings(nb)%num_facade_h(k) + 1
712	!
713	!-- Moreover, sum up number of local facade elements per building.
714	buildings(nb)%num_facades_per_building_h_l = &
715	buildings(nb)%num_facades_per_building_h_l + 1
716	ENDIF
717	ENDDO
718	!
719	!-- Vertical facades
720	buildings(:)%num_facades_per_building_v_l = 0
721	DO l = 0, 3
722	DO m = 1, surf_usm_v(l)%ns
723	!
724	!-- For the current facade element determine corresponding building index.
725	!-- First, obtain j,j,k indices of the building. Please note the
726	!-- offset between facade/surface element and building location (for
727	!-- vertical surface elements the vertical offsets are zero).
728	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
729	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
730	k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
731
732	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
733	DIM = 1 )
734	IF ( buildings(nb)%on_pe ) THEN
735	buildings(nb)%num_facade_v(k) = buildings(nb)%num_facade_v(k) + 1
736	buildings(nb)%num_facades_per_building_v_l = &
737	buildings(nb)%num_facades_per_building_v_l + 1
738	ENDIF
739	ENDDO
740	ENDDO
741
742	!
743	!-- Determine total number of facade elements per building and assign number to
744	!-- building data type.
745	DO nb = 1, num_build
746	!
747	!-- Allocate dummy array used for summing-up facade elements.
748	!-- Please note, dummy arguments are necessary as building-date type
749	!-- arrays are not necessarily allocated on all PEs.
750	ALLOCATE( num_facades_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
751	ALLOCATE( num_facades_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
752	ALLOCATE( receive_dum_h(buildings(nb)%kb_min:buildings(nb)%kb_max) )
753	ALLOCATE( receive_dum_v(buildings(nb)%kb_min:buildings(nb)%kb_max) )
754	num_facades_h = 0
755	num_facades_v = 0
756	receive_dum_h = 0
757	receive_dum_v = 0
758
759	IF ( buildings(nb)%on_pe ) THEN
760	num_facades_h = buildings(nb)%num_facade_h
761	num_facades_v = buildings(nb)%num_facade_v
762	ENDIF
763
764	#if defined( __parallel )
765	CALL MPI_ALLREDUCE( num_facades_h, &
766	receive_dum_h, &
767	buildings(nb)%kb_max - buildings(nb)%kb_min + 1, &
768	MPI_INTEGER, &
769	MPI_SUM, &
770	comm2d, &
771	ierr )
772
773	CALL MPI_ALLREDUCE( num_facades_v, &
774	receive_dum_v, &
775	buildings(nb)%kb_max - buildings(nb)%kb_min + 1, &
776	MPI_INTEGER, &
777	MPI_SUM, &
778	comm2d, &
779	ierr )
780	IF ( ALLOCATED( buildings(nb)%num_facade_h ) ) & !FK: Was wenn not allocated? --> s.o.
781	buildings(nb)%num_facade_h = receive_dum_h
782	IF ( ALLOCATED( buildings(nb)%num_facade_v ) ) &
783	buildings(nb)%num_facade_v = receive_dum_v
784	#else
785	buildings(nb)%num_facade_h = num_facades_h
786	buildings(nb)%num_facade_v = num_facades_v
787	#endif
788	!
789	!-- Deallocate dummy arrays
790	DEALLOCATE( num_facades_h )
791	DEALLOCATE( num_facades_v )
792	DEALLOCATE( receive_dum_h )
793	DEALLOCATE( receive_dum_v )
794	!
795	!-- Allocate index arrays which link facade elements with surface-data type.
796	!-- Please note, no height levels are considered here (information is stored
797	!-- in surface-data type itself).
798	IF ( buildings(nb)%on_pe ) THEN
799	!
800	!-- Determine number of facade elements per building.
801	buildings(nb)%num_facades_per_building_h = SUM( buildings(nb)%num_facade_h )
802	buildings(nb)%num_facades_per_building_v = SUM( buildings(nb)%num_facade_v )
803	!
804	!-- Allocate arrays which link the building with the horizontal and vertical
805	!-- urban-type surfaces. Please note, linking arrays are allocated over all
806	!-- facade elements, which is required in case a building is located at the
807	!-- subdomain boundaries, where the building and the corresponding surface
808	!-- elements are located on different subdomains.
809	ALLOCATE( buildings(nb)%m_h(1:buildings(nb)%num_facades_per_building_h_l) )
810
811	ALLOCATE( buildings(nb)%l_v(1:buildings(nb)%num_facades_per_building_v_l) )
812	ALLOCATE( buildings(nb)%m_v(1:buildings(nb)%num_facades_per_building_v_l) )
813	ENDIF
814	!
815	!-- Determine volume per facade element (vpf)
816	IF ( buildings(nb)%on_pe ) THEN
817	ALLOCATE( buildings(nb)%vpf(buildings(nb)%kb_min:buildings(nb)%kb_max) )
818
819	DO k = buildings(nb)%kb_min, buildings(nb)%kb_max
820	buildings(nb)%vpf(k) = buildings(nb)%volume(k) / &
821	( buildings(nb)%num_facade_h(k) + &
822	buildings(nb)%num_facade_v(k) )
823	ENDDO
824	ENDIF
825	ENDDO
826	!
827	!-- Link facade elements with surface data type.
828	!-- Allocate array for counting.
829	ALLOCATE( n_fa(1:num_build) )
830	n_fa = 1
831
832	DO m = 1, surf_usm_h%ns
833	i = surf_usm_h%i(m) + surf_usm_h%ioff
834	j = surf_usm_h%j(m) + surf_usm_h%joff
835
836	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
837
838	IF ( buildings(nb)%on_pe ) THEN
839	buildings(nb)%m_h(n_fa(nb)) = m
840	n_fa(nb) = n_fa(nb) + 1
841	ENDIF
842	ENDDO
843
844	n_fa = 1
845	DO l = 0, 3
846	DO m = 1, surf_usm_v(l)%ns
847	i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
848	j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
849
850	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), DIM = 1 )
851
852	IF ( buildings(nb)%on_pe ) THEN
853	buildings(nb)%l_v(n_fa(nb)) = l
854	buildings(nb)%m_v(n_fa(nb)) = m
855	n_fa(nb) = n_fa(nb) + 1
856	ENDIF
857	ENDDO
858	ENDDO
859	DEALLOCATE( n_fa )
860
861	!
862	!-- Initialize building parameters, first by mean building type. Note,
863	!-- in this case all buildings have the same type.
864	!-- In a second step initialize with building tpyes from static input file,
865	!-- where building types can be individual for each building.
866	buildings(:)%lambda_layer3 = building_pars(63,building_type)
867	buildings(:)%s_layer3 = building_pars(57,building_type)
868	buildings(:)%f_c_win = building_pars(119,building_type)
869	buildings(:)%g_value_win = building_pars(120,building_type)
870	buildings(:)%u_value_win = building_pars(121,building_type)
871	buildings(:)%air_change_low = building_pars(122,building_type)
872	buildings(:)%air_change_high = building_pars(123,building_type)
873	buildings(:)%eta_ve = building_pars(124,building_type)
874	buildings(:)%factor_a = building_pars(125,building_type)
875	buildings(:)%factor_c = building_pars(126,building_type)
876	buildings(:)%lambda_at = building_pars(127,building_type)
877	buildings(:)%theta_int_h_set = building_pars(118,building_type)
878	buildings(:)%theta_int_c_set = building_pars(117,building_type)
879	buildings(:)%phi_h_max = building_pars(128,building_type)
880	buildings(:)%phi_c_max = building_pars(129,building_type)
881	buildings(:)%qint_high = building_pars(130,building_type)
882	buildings(:)%qint_low = building_pars(131,building_type)
883	buildings(:)%height_storey = building_pars(132,building_type)
884	buildings(:)%height_cei_con = building_pars(133,building_type)
885	!
886	!-- Initialize ventilaation load. Please note, building types > 7 are actually
887	!-- not allowed (check already in urban_surface_mod and netcdf_data_input_mod.
888	!-- However, the building data base may be later extended.
889	IF ( building_type == 1 .OR. building_type == 2 .OR. &
890	building_type == 3 .OR. building_type == 10 .OR. &
891	building_type == 11 .OR. building_type == 12 ) THEN
892	buildings(nb)%ventilation_int_loads = 1
893	!
894	!-- Office, building with large windows
895	ELSEIF ( building_type == 4 .OR. building_type == 5 .OR. &
896	building_type == 6 .OR. building_type == 7 .OR. &
897	building_type == 8 .OR. building_type == 9) THEN
898	buildings(nb)%ventilation_int_loads = 2
899	!
900	!-- Industry, hospitals
901	ELSEIF ( building_type == 13 .OR. building_type == 14 .OR. &
902	building_type == 15 .OR. building_type == 16 .OR. &
903	building_type == 17 .OR. building_type == 18 ) THEN
904	buildings(nb)%ventilation_int_loads = 3
905	ENDIF
906	!
907	!-- Initialization of building parameters - level 2
908	IF ( building_type_f%from_file ) THEN
909	DO i = nxl, nxr
910	DO j = nys, nyn
911	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
912	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
913	DIM = 1 )
914	bt = building_type_f%var(j,i)
915
916	buildings(nb)%lambda_layer3 = building_pars(63,bt)
917	buildings(nb)%s_layer3 = building_pars(57,bt)
918	buildings(nb)%f_c_win = building_pars(119,bt)
919	buildings(nb)%g_value_win = building_pars(120,bt)
920	buildings(nb)%u_value_win = building_pars(121,bt)
921	buildings(nb)%air_change_low = building_pars(122,bt)
922	buildings(nb)%air_change_high = building_pars(123,bt)
923	buildings(nb)%eta_ve = building_pars(124,bt)
924	buildings(nb)%factor_a = building_pars(125,bt)
925	buildings(nb)%factor_c = building_pars(126,bt)
926	buildings(nb)%lambda_at = building_pars(127,bt)
927	buildings(nb)%theta_int_h_set = building_pars(118,bt)
928	buildings(nb)%theta_int_c_set = building_pars(117,bt)
929	buildings(nb)%phi_h_max = building_pars(128,bt)
930	buildings(nb)%phi_c_max = building_pars(129,bt)
931	buildings(nb)%qint_high = building_pars(130,bt)
932	buildings(nb)%qint_low = building_pars(131,bt)
933	buildings(nb)%height_storey = building_pars(132,bt)
934	buildings(nb)%height_cei_con = building_pars(133,bt)
935	!
936	!-- Initialize ventilaation load. Please note, building types > 7
937	!-- are actually not allowed (check already in urban_surface_mod
938	!-- and netcdf_data_input_mod. However, the building data base may
939	!-- be later extended.
940	IF ( bt == 1 .OR. bt == 2 .OR. &
941	bt == 3 .OR. bt == 10 .OR. &
942	bt == 11 .OR. bt == 12 ) THEN
943	buildings(nb)%ventilation_int_loads = 1
944	!
945	!-- Office, building with large windows
946	ELSEIF ( bt == 4 .OR. bt == 5 .OR. &
947	bt == 6 .OR. bt == 7 .OR. &
948	bt == 8 .OR. bt == 9) THEN
949	buildings(nb)%ventilation_int_loads = 2
950	!
951	!-- Industry, hospitals
952	ELSEIF ( bt == 13 .OR. bt == 14 .OR. &
953	bt == 15 .OR. bt == 16 .OR. &
954	bt == 17 .OR. bt == 18 ) THEN
955	buildings(nb)%ventilation_int_loads = 3
956	ENDIF
957	ENDIF
958	ENDDO
959	ENDDO
960	ENDIF
961	!
962	!-- Initial room temperature [K]
963	!-- (after first loop, use theta_m_t as theta_m_t_prev)
964	theta_m_t_prev = initial_indoor_temperature
965	!
966	!-- Initialize indoor temperature. Actually only for output at initial state.
967	DO nb = 1, num_build
968	buildings(nb)%t_in(:) = initial_indoor_temperature
969	ENDDO
970
971	CALL location_message( 'finished', .TRUE. )
972
973	END SUBROUTINE im_init
974
975
976	!------------------------------------------------------------------------------!
977	! Description:
978	! ------------
979	!> Main part of the indoor model.
980	!> Calculation of .... (kanani: Please describe)
981	!------------------------------------------------------------------------------!
982	SUBROUTINE im_main_heatcool
983
984	USE arrays_3d, &
985	ONLY: ddzw, dzw
986
987	USE basic_constants_and_equations_mod, &
988	ONLY: c_p
989
990	USE control_parameters, &
991	ONLY: rho_surface
992
993	USE date_and_time_mod, &
994	ONLY: time_utc
995
996	USE grid_variables, &
997	ONLY: dx, dy
998
999	USE pegrid
1000
1001	USE surface_mod, &
1002	ONLY: ind_veg_wall, ind_wat_win, surf_usm_h, surf_usm_v
1003
1004	USE urban_surface_mod, &
1005	ONLY: nzt_wall, t_wall_h, t_wall_v, t_window_h, t_window_v, &
1006	building_type
1007
1008
1009	IMPLICIT NONE
1010
1011	INTEGER(iwp) :: i !< index of facade-adjacent atmosphere grid point in x-direction
1012	INTEGER(iwp) :: j !< index of facade-adjacent atmosphere grid point in y-direction
1013	INTEGER(iwp) :: k !< index of facade-adjacent atmosphere grid point in z-direction
1014	INTEGER(iwp) :: kk !< vertical index of indoor grid point adjacent to facade
1015	INTEGER(iwp) :: l !< running index for surface-element orientation
1016	INTEGER(iwp) :: m !< running index surface elements
1017	INTEGER(iwp) :: nb !< running index for buildings
1018	INTEGER(iwp) :: fa !< running index for facade elements of each building
1019
1020	REAL(wp) :: indoor_wall_window_temperature !< weighted temperature of innermost wall/window layer
1021	REAL(wp) :: near_facade_temperature !< outside air temperature 10cm away from facade
1022	REAL(wp) :: time_utc_hour !< time of day (hour UTC)
1023
1024	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_l_send !< dummy send buffer used for summing-up indoor temperature per kk-level
1025	REAL(wp), DIMENSION(:), ALLOCATABLE :: t_in_recv !< dummy recv buffer used for summing-up indoor temperature per kk-level
1026	!
1027	!-- Determine time of day in hours.
1028	time_utc_hour = time_utc / 3600.0_wp
1029	!
1030	!-- Following calculations must be done for each facade element.
1031	DO nb = 1, num_build
1032	!
1033	!-- First, check whether building is present on local subdomain.
1034	IF ( buildings(nb)%on_pe ) THEN
1035	!
1036	!-- Determine daily schedule. 08:00-18:00 = 1, other hours = 0.
1037	!-- Residental Building, panel WBS 70
1038	IF ( buildings(nb)%ventilation_int_loads == 1 ) THEN
1039	IF ( time_utc_hour >= 6.0_wp .AND. time_utc_hour <= 8.0_wp ) THEN
1040	schedule_d = 1
1041	ELSEIF ( time_utc_hour >= 18.0_wp .AND. time_utc_hour <= 23.0_wp ) THEN
1042	schedule_d = 1
1043	ELSE
1044	schedule_d = 0
1045	ENDIF
1046	ENDIF
1047	!
1048	!-- Office, building with large windows
1049	IF ( buildings(nb)%ventilation_int_loads == 2 ) THEN
1050	IF ( time_utc_hour >= 8.0_wp .AND. time_utc_hour <= 18.0_wp ) THEN
1051	schedule_d = 1
1052	ELSE
1053	schedule_d = 0
1054	ENDIF
1055	ENDIF
1056	!
1057	!-- Industry, hospitals
1058	IF ( buildings(nb)%ventilation_int_loads == 3 ) THEN
1059	IF ( time_utc_hour >= 6.0_wp .AND. time_utc_hour <= 22.0_wp ) THEN
1060	schedule_d = 1
1061	ELSE
1062	schedule_d = 0
1063	ENDIF
1064	ENDIF
1065	!
1066	!-- Initialize/reset indoor temperature
1067	buildings(nb)%t_in_l = 0.0_wp
1068	!
1069	!-- Horizontal surfaces
1070	DO fa = 1, buildings(nb)%num_facades_per_building_h_l
1071	!
1072	!-- Determine index where corresponding surface-type information
1073	!-- is stored.
1074	m = buildings(nb)%m_h(fa)
1075	!
1076	!-- Determine building height level index.
1077	kk = surf_usm_h%k(m) + surf_usm_h%koff
1078	!
1079	!-- Building geometries --> not time-dependent
1080	facade_element_area = dx * dy !< [m2] surface area per facade element
1081	floor_area_per_facade = buildings(nb)%vpf(kk) * ddzw(kk) !< [m2] net floor area per facade element
1082	indoor_volume_per_facade = buildings(nb)%vpf(kk) !< [m3] indoor air volume per facade element
1083	window_area_per_facade = surf_usm_h%frac(ind_wat_win,m) * facade_element_area !< [m2] window area per facade element
1084
1085	! print*, "building_height", building_height
1086	! print*, "num_facades_v_l", buildings(nb)%num_facades_per_building_v_l
1087	! print*, "num_facades_v", buildings(nb)%num_facades_per_building_v
1088	! print*, "kb_min_max", buildings(nb)%kb_min, buildings(nb)%kb_max
1089	! print*, "dzw kk", dzw(kk), kk
1090
1091	f_cei = buildings(nb)%building_height / &
1092	(buildings(nb)%height_storey-buildings(nb)%height_cei_con) !< [-] factor for ceiling redcution
1093	ngs = buildings(nb)%vpf(kk)/f_cei !< [m2] calculation of netto ground surface
1094	f_sr = ngs/floor_area_per_facade !< [-] factor for surface reduction
1095	eff_mass_area = buildings(nb)%factor_a * ngs !< [m2] standard values according to Table 12 section 12.3.1.2 (calculate over Eq. (65) according to section 12.3.1.2)
1096	c_m = buildings(nb)%factor_c * ngs !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
1097	total_area = buildings(nb)%lambda_at * floor_area_per_facade !< [m2] area of all surfaces pointing to zone Eq. (9) according to section 7.2.2.2
1098
1099	!-- Calculation of heat transfer coefficient for transmission --> not time-dependent
1100	h_tr_w = window_area_per_facade * buildings(nb)%u_value_win !< [W/K] only for windows
1101	h_tr_is = total_area * h_is !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
1102	h_tr_ms = eff_mass_area * h_ms !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
1103	h_tr_op = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade ) &
1104	* buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp ) + 1.0_wp / h_tr_ms )
1105	h_tr_em = 1.0_wp / ( 1.0_wp / h_tr_op - 1.0_wp / h_tr_ms ) !< [W/K] Eq. (63), Section 12.2.2
1106	!
1107	!-- internal air loads dependent on the occupacy of the room
1108	!-- basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
1109	phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low ) &
1110	* ngs ) !< [W] Eq. (C.1)
1111	!
1112	!-- Airflow dependent on the occupacy of the room
1113	!-- basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
1114	air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low ) !< [1/h]?
1115	!
1116	!-- Heat transfer of ventilation
1117	!-- not less than 0.01 W/K to provide division by 0 in further calculations
1118	!-- with heat capacity of air 0.33 Wh/m2K
1119	h_ve = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade * &
1120	0.33_wp * (1.0_wp - buildings(nb)%eta_ve ) ) ) !< [W/K] from ISO 13789 Eq.(10)
1121
1122	!-- Heat transfer coefficient auxiliary variables
1123	h_tr_1 = 1.0_wp / ( ( 1.0_wp / h_ve ) + ( 1.0_wp / h_tr_is ) ) !< [W/K] Eq. (C.6)
1124	h_tr_2 = h_tr_1 + h_tr_w !< [W/K] Eq. (C.7)
1125	h_tr_3 = 1.0_wp / ( ( 1.0_wp / h_tr_2 ) + ( 1.0_wp / h_tr_ms ) ) !< [W/K] Eq. (C.8)
1126	!
1127	!-- Net short-wave radiation through window area (was i_global)
1128	net_sw_in = surf_usm_h%rad_sw_in(m) - surf_usm_h%rad_sw_out(m)
1129	!
1130	!-- Quantities needed for im_calc_temperatures
1131	i = surf_usm_h%i(m)
1132	j = surf_usm_h%j(m)
1133	k = surf_usm_h%k(m)
1134	near_facade_temperature = surf_usm_h%pt_10cm(m)
1135	indoor_wall_window_temperature = &
1136	surf_usm_h%frac(ind_veg_wall,m) * t_wall_h(nzt_wall,m) &
1137	+ surf_usm_h%frac(ind_wat_win,m) * t_window_h(nzt_wall,m)
1138	!
1139	!-- Solar thermal gains. If net_sw_in larger than sun-protection
1140	!-- threshold parameter (params_solar_protection), sun protection will
1141	!-- be activated
1142	IF ( net_sw_in <= params_solar_protection ) THEN
1143	solar_protection_off = 1
1144	solar_protection_on = 0
1145	ELSE
1146	solar_protection_off = 0
1147	solar_protection_on = 1
1148	ENDIF
1149	!
1150	!-- Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
1151	!-- DIN 4108 - 2 chap.8
1152	phi_sol = ( window_area_per_facade * net_sw_in * solar_protection_off &
1153	+ window_area_per_facade * net_sw_in * buildings(nb)%f_c_win * solar_protection_on ) &
1154	* buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w !< [W]
1155	!
1156	!-- Calculation of the mass specific thermal load for internal and external heatsources of the inner node
1157	phi_m = (eff_mass_area / total_area) * ( phi_ia + phi_sol ) !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
1158	!
1159	!-- Calculation mass specific thermal load implied non thermal mass
1160	phi_st = ( 1.0_wp - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1_wp * total_area ) ) ) &
1161	* ( phi_ia + phi_sol ) !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
1162	!
1163	!-- Calculations for deriving indoor temperature and heat flux into the wall
1164	!-- Step 1: Indoor temperature without heating and cooling
1165	!-- section C.4.1 Picture C.2 zone 3)
1166	phi_hc_nd = 0.0_wp
1167
1168	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1169	near_facade_temperature, phi_hc_nd )
1170	!
1171	!-- If air temperature between border temperatures of heating and cooling, assign output variable, then ready
1172	IF ( buildings(nb)%theta_int_h_set <= theta_air .AND. theta_air <= buildings(nb)%theta_int_c_set ) THEN
1173	phi_hc_nd_ac = 0.0_wp
1174	phi_hc_nd = phi_hc_nd_ac
1175	theta_air_ac = theta_air
1176	!
1177	!-- Step 2: Else, apply 10 W/m2 heating/cooling power and calculate indoor temperature
1178	!-- again.
1179	ELSE
1180	!
1181	!-- Temperature not correct, calculation method according to section C4.2
1182	theta_air_0 = theta_air !< Note temperature without heating/cooling
1183
1184	!-- Heating or cooling?
1185	IF ( theta_air > buildings(nb)%theta_int_c_set ) THEN
1186	theta_air_set = buildings(nb)%theta_int_c_set
1187	ELSE
1188	theta_air_set = buildings(nb)%theta_int_h_set
1189	ENDIF
1190
1191	!-- Calculate the temperature with phi_hc_nd_10
1192	phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
1193	phi_hc_nd = phi_hc_nd_10
1194
1195	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1196	near_facade_temperature, phi_hc_nd )
1197
1198	theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
1199	!
1200
1201	phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0) &
1202	/ (theta_air_10 - theta_air_0) !< Eq. (C.13)
1203
1204	!-- Step 3: With temperature ratio to determine the heating or cooling capacity
1205	!-- If necessary, limit the power to maximum power
1206	!-- section C.4.1 Picture C.2 zone 2) and 4)
1207	IF ( buildings(nb)%phi_c_max < phi_hc_nd_un .AND. phi_hc_nd_un < buildings(nb)%phi_h_max ) THEN
1208	phi_hc_nd_ac = phi_hc_nd_un
1209	phi_hc_nd = phi_hc_nd_un
1210	ELSE
1211	!-- Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
1212	!-- section C.4.1 Picture C.2 zone 1) and 5)
1213	IF ( phi_hc_nd_un > 0.0_wp ) THEN
1214	phi_hc_nd_ac = buildings(nb)%phi_h_max !< Limit heating
1215	ELSE
1216	phi_hc_nd_ac = buildings(nb)%phi_c_max !< Limit cooling
1217	ENDIF
1218	ENDIF
1219
1220	phi_hc_nd = phi_hc_nd_ac
1221	!
1222	!-- Calculate the temperature with phi_hc_nd_ac (new)
1223	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1224	near_facade_temperature, phi_hc_nd )
1225
1226	theta_air_ac = theta_air
1227
1228	ENDIF
1229	!
1230	!-- Update theta_m_t_prev
1231	theta_m_t_prev = theta_m_t
1232	!
1233	!-- Calculate the operating temperature with weighted mean temperature of air and mean solar temperature
1234	!-- Will be used for thermal comfort calculations
1235	theta_op = 0.3_wp * theta_air_ac + 0.7_wp * theta_s !< [degree_C] operative Temperature Eq. (C.12)
1236	!
1237	!-- Heat flux into the wall. Value needed in urban_surface_mod to
1238	!-- calculate heat transfer through wall layers towards the facade
1239	!-- (use c_p * rho_surface to convert [W/m2] into [K m/s])
1240	q_wall_win = h_tr_ms * ( theta_s - theta_m ) &
1241	/ ( facade_element_area &
1242	- window_area_per_facade )
1243	!
1244	!-- Transfer q_wall_win back to USM (innermost wall/window layer)
1245	surf_usm_h%iwghf_eb(m) = q_wall_win
1246	surf_usm_h%iwghf_eb_window(m) = q_wall_win
1247	!
1248	!-- Sum up operational indoor temperature per kk-level. Further below,
1249	!-- this temperature is reduced by MPI to one temperature per kk-level
1250	!-- and building (processor overlapping)
1251	buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
1252	!
1253	!-- Calculation of waste heat
1254	!-- Anthropogenic heat output
1255	IF ( phi_hc_nd_ac > 0.0_wp ) THEN
1256	heating_on = 1
1257	cooling_on = 0
1258	ELSE
1259	heating_on = 0
1260	cooling_on = 1
1261	ENDIF
1262
1263	q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on))!< [W/GebÃ€udemodell] , observe the directional convention in PALM!
1264	surf_usm_h%waste_heat(m) = q_waste_heat
1265
1266	ENDDO !< Horizontal surfaces loop
1267	!
1268	!-- Vertical surfaces
1269	DO fa = 1, buildings(nb)%num_facades_per_building_v_l
1270	!
1271	!-- Determine indices where corresponding surface-type information
1272	!-- is stored.
1273	l = buildings(nb)%l_v(fa)
1274	m = buildings(nb)%m_v(fa)
1275	!
1276	!-- Determine building height level index.
1277	kk = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
1278	!
1279	!-- (SOME OF THE FOLLOWING (not time-dependent COULD PROBABLY GO INTO A FUNCTION
1280	!-- EXCEPT facade_element_area, EVERYTHING IS CALCULATED EQUALLY)
1281	!-- Building geometries --> not time-dependent
1282	IF ( l == 0 .OR. l == 1 ) facade_element_area = dx * dzw(kk) !< [m2] surface area per facade element
1283	IF ( l == 2 .OR. l == 3 ) facade_element_area = dy * dzw(kk) !< [m2] surface area per facade element
1284	floor_area_per_facade = buildings(nb)%vpf(kk) * ddzw(kk) !< [m2] net floor area per facade element
1285	indoor_volume_per_facade = buildings(nb)%vpf(kk) !< [m3] indoor air volume per facade element
1286	window_area_per_facade = surf_usm_v(l)%frac(ind_wat_win,m) * facade_element_area !< [m2] window area per facade element
1287
1288	f_cei = buildings(nb)%building_height / &
1289	(buildings(nb)%height_storey-buildings(nb)%height_cei_con) !< [-] factor for ceiling redcution
1290	ngs = buildings(nb)%vpf(kk)/f_cei !< [m2] calculation of netto ground surface
1291	f_sr = ngs/floor_area_per_facade !< [-] factor for surface reduction
1292	eff_mass_area = buildings(nb)%factor_a * ngs !< [m2] standard values according to Table 12 section 12.3.1.2 (calculate over Eq. (65) according to section 12.3.1.2)
1293	c_m = buildings(nb)%factor_c * ngs !< [J/K] standard values according to table 12 section 12.3.1.2 (calculate over Eq. (66) according to section 12.3.1.2)
1294	total_area = buildings(nb)%lambda_at * floor_area_per_facade !< [m2] area of all surfaces pointing to zone Eq. (9) according to section 7.2.2.2
1295	!
1296	!-- Calculation of heat transfer coefficient for transmission --> not time-dependent
1297	h_tr_w = window_area_per_facade * buildings(nb)%u_value_win !< [W/K] only for windows
1298	h_tr_is = total_area * h_is !< [W/K] with h_is = 3.45 W / (m2 K) between surface and air, Eq. (9)
1299	h_tr_ms = eff_mass_area * h_ms !< [W/K] with h_ms = 9.10 W / (m2 K) between component and surface, Eq. (64)
1300	h_tr_op = 1.0_wp / ( 1.0_wp / ( ( facade_element_area - window_area_per_facade ) &
1301	* buildings(nb)%lambda_layer3 / buildings(nb)%s_layer3 * 0.5_wp ) + 1.0_wp / h_tr_ms )
1302	h_tr_em = 1.0_wp / ( 1.0_wp / h_tr_op - 1.0_wp / h_tr_ms ) !< [W/K] Eq. (63), Section 12.2.2
1303	!
1304	!-- internal air loads dependent on the occupacy of the room
1305	!-- basical internal heat gains (qint_low) with additional internal heat gains by occupancy (qint_high) (0,5*phi_int)
1306	phi_ia = 0.5_wp * ( ( buildings(nb)%qint_high * schedule_d + buildings(nb)%qint_low ) &
1307	* ngs ) !< [W] Eq. (C.1)
1308	!
1309	!-- Airflow dependent on the occupacy of the room
1310	!-- basical airflow (air_change_low) with additional airflow gains by occupancy (air_change_high)
1311	air_change = ( buildings(nb)%air_change_high * schedule_d + buildings(nb)%air_change_low )
1312	!
1313	!-- Heat transfer of ventilation
1314	!-- not less than 0.01 W/K to provide division by 0 in further calculations
1315	!-- with heat capacity of air 0.33 Wh/m2K
1316	h_ve = MAX( 0.01_wp , ( air_change * indoor_volume_per_facade * &
1317	0.33_wp * (1 - buildings(nb)%eta_ve ) ) ) !< [W/K] from ISO 13789 Eq.(10)
1318
1319	!-- Heat transfer coefficient auxiliary variables
1320	h_tr_1 = 1.0_wp / ( ( 1.0_wp / h_ve ) + ( 1.0_wp / h_tr_is ) ) !< [W/K] Eq. (C.6)
1321	h_tr_2 = h_tr_1 + h_tr_w !< [W/K] Eq. (C.7)
1322	h_tr_3 = 1.0_wp / ( ( 1.0_wp / h_tr_2 ) + ( 1.0_wp / h_tr_ms ) ) !< [W/K] Eq. (C.8)
1323	!
1324	!-- Net short-wave radiation through window area (was i_global)
1325	net_sw_in = surf_usm_v(l)%rad_sw_in(m) - surf_usm_v(l)%rad_sw_out(m)
1326	!
1327	!-- Quantities needed for im_calc_temperatures
1328	i = surf_usm_v(l)%i(m)
1329	j = surf_usm_v(l)%j(m)
1330	k = surf_usm_v(l)%k(m)
1331	near_facade_temperature = surf_usm_v(l)%pt_10cm(m)
1332	indoor_wall_window_temperature = &
1333	surf_usm_v(l)%frac(ind_veg_wall,m) * t_wall_v(l)%t(nzt_wall,m) &
1334	+ surf_usm_v(l)%frac(ind_wat_win,m) * t_window_v(l)%t(nzt_wall,m)
1335	!
1336	!-- Solar thermal gains. If net_sw_in larger than sun-protection
1337	!-- threshold parameter (params_solar_protection), sun protection will
1338	!-- be activated
1339	IF ( net_sw_in <= params_solar_protection ) THEN
1340	solar_protection_off = 1
1341	solar_protection_on = 0
1342	ELSE
1343	solar_protection_off = 0
1344	solar_protection_on = 1
1345	ENDIF
1346	!
1347	!-- Calculation of total heat gains from net_sw_in through windows [W] in respect on automatic sun protection
1348	!-- DIN 4108 - 2 chap.8
1349	phi_sol = ( window_area_per_facade * net_sw_in * solar_protection_off &
1350	+ window_area_per_facade * net_sw_in * buildings(nb)%f_c_win * solar_protection_on ) &
1351	* buildings(nb)%g_value_win * ( 1.0_wp - params_f_f ) * params_f_w
1352	!
1353	!-- Calculation of the mass specific thermal load for internal and external heatsources
1354	phi_m = (eff_mass_area / total_area) * ( phi_ia + phi_sol ) !< [W] Eq. (C.2) with phi_ia=0,5*phi_int
1355	!
1356	!-- Calculation mass specific thermal load implied non thermal mass
1357	phi_st = ( 1.0_wp - ( eff_mass_area / total_area ) - ( h_tr_w / ( 9.1_wp * total_area ) ) ) &
1358	* ( phi_ia + phi_sol ) !< [W] Eq. (C.3) with phi_ia=0,5*phi_int
1359	!
1360	!-- Calculations for deriving indoor temperature and heat flux into the wall
1361	!-- Step 1: Indoor temperature without heating and cooling
1362	!-- section C.4.1 Picture C.2 zone 3)
1363	phi_hc_nd = 0.0_wp
1364	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1365	near_facade_temperature, phi_hc_nd )
1366	!
1367	!-- If air temperature between border temperatures of heating and cooling, assign output variable, then ready
1368	IF ( buildings(nb)%theta_int_h_set <= theta_air .AND. theta_air <= buildings(nb)%theta_int_c_set ) THEN
1369	phi_hc_nd_ac = 0.0_wp
1370	phi_hc_nd = phi_hc_nd_ac
1371	theta_air_ac = theta_air
1372	!
1373	!-- Step 2: Else, apply 10 W/m2 heating/cooling power and calculate indoor temperature
1374	!-- again.
1375	ELSE
1376	!
1377	!-- Temperature not correct, calculation method according to section C4.2
1378	theta_air_0 = theta_air !< Note temperature without heating/cooling
1379
1380	!-- Heating or cooling?
1381	IF ( theta_air > buildings(nb)%theta_int_c_set ) THEN
1382	theta_air_set = buildings(nb)%theta_int_c_set
1383	ELSE
1384	theta_air_set = buildings(nb)%theta_int_h_set
1385	ENDIF
1386
1387	!-- Calculate the temperature with phi_hc_nd_10
1388	phi_hc_nd_10 = 10.0_wp * floor_area_per_facade
1389	phi_hc_nd = phi_hc_nd_10
1390
1391	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1392	near_facade_temperature, phi_hc_nd )
1393
1394	theta_air_10 = theta_air !< Note the temperature with 10 W/m2 of heating
1395
1396
1397	phi_hc_nd_un = phi_hc_nd_10 * (theta_air_set - theta_air_0) &
1398	/ (theta_air_10 - theta_air_0) !< Eq. (C.13)
1399	!
1400	!-- Step 3: With temperature ratio to determine the heating or cooling capacity
1401	!-- If necessary, limit the power to maximum power
1402	!-- section C.4.1 Picture C.2 zone 2) and 4)
1403	IF ( buildings(nb)%phi_c_max < phi_hc_nd_un .AND. phi_hc_nd_un < buildings(nb)%phi_h_max ) THEN
1404	phi_hc_nd_ac = phi_hc_nd_un
1405	phi_hc_nd = phi_hc_nd_un
1406	ELSE
1407	!-- Step 4: Inner temperature with maximum heating (phi_hc_nd_un positive) or cooling (phi_hc_nd_un negative)
1408	!-- section C.4.1 Picture C.2 zone 1) and 5)
1409	IF ( phi_hc_nd_un > 0.0_wp ) THEN
1410	phi_hc_nd_ac = buildings(nb)%phi_h_max !< Limit heating
1411	ELSE
1412	phi_hc_nd_ac = buildings(nb)%phi_c_max !< Limit cooling
1413	ENDIF
1414	ENDIF
1415
1416	phi_hc_nd = phi_hc_nd_ac
1417	!
1418	!-- Calculate the temperature with phi_hc_nd_ac (new)
1419	CALL im_calc_temperatures ( i, j, k, indoor_wall_window_temperature, &
1420	near_facade_temperature, phi_hc_nd )
1421
1422	theta_air_ac = theta_air
1423
1424	ENDIF
1425	!
1426	!-- Update theta_m_t_prev
1427	theta_m_t_prev = theta_m_t
1428	!
1429	!-- Calculate the operating temperature with weighted mean of temperature of air and mean
1430	!-- Will be used for thermal comfort calculations
1431	theta_op = 0.3_wp * theta_air_ac + 0.7_wp * theta_s
1432	!
1433	!-- Heat flux into the wall. Value needed in urban_surface_mod to
1434	!-- calculate heat transfer through wall layers towards the facade
1435	q_wall_win = h_tr_ms * ( theta_s - theta_m ) &
1436	/ ( facade_element_area &
1437	- window_area_per_facade )
1438	!
1439	!-- Transfer q_wall_win back to USM (innermost wall/window layer)
1440	surf_usm_v(l)%iwghf_eb(m) = q_wall_win
1441	surf_usm_v(l)%iwghf_eb_window(m) = q_wall_win
1442	!
1443	!-- Sum up operational indoor temperature per kk-level. Further below,
1444	!-- this temperature is reduced by MPI to one temperature per kk-level
1445	!-- and building (processor overlapping)
1446	buildings(nb)%t_in_l(kk) = buildings(nb)%t_in_l(kk) + theta_op
1447
1448	!
1449	!-- Calculation of waste heat
1450	!-- Anthropogenic heat output
1451	IF ( phi_hc_nd_ac > 0.0_wp ) THEN
1452	heating_on = 1
1453	cooling_on = 0
1454	ELSE
1455	heating_on = 0
1456	cooling_on = 1
1457	ENDIF
1458
1459	q_waste_heat = (phi_hc_nd * (params_waste_heat_h * heating_on + params_waste_heat_c * cooling_on))!< [W/GebÃ€udemodell] , observe the directional convention in PALM!
1460	surf_usm_v(l)%waste_heat(m) = q_waste_heat
1461
1462	ENDDO !< Vertical surfaces loop
1463
1464	ENDIF !< buildings(nb)%on_pe
1465	ENDDO !< buildings loop
1466
1467	!
1468	!-- Determine the mean building temperature.
1469	DO nb = 1, num_build
1470	!
1471	!-- Allocate dummy array used for summing-up facade elements.
1472	!-- Please note, dummy arguments are necessary as building-date type
1473	!-- arrays are not necessarily allocated on all PEs.
1474	ALLOCATE( t_in_l_send(buildings(nb)%kb_min:buildings(nb)%kb_max) )
1475	ALLOCATE( t_in_recv(buildings(nb)%kb_min:buildings(nb)%kb_max) )
1476	t_in_l_send = 0.0_wp
1477	t_in_recv = 0.0_wp
1478
1479	IF ( buildings(nb)%on_pe ) THEN
1480	t_in_l_send = buildings(nb)%t_in_l
1481	ENDIF
1482
1483	#if defined( __parallel )
1484	CALL MPI_ALLREDUCE( t_in_l_send, &
1485	t_in_recv, &
1486	buildings(nb)%kb_max - buildings(nb)%kb_min + 1, &
1487	MPI_REAL, &
1488	MPI_SUM, &
1489	comm2d, &
1490	ierr )
1491
1492	IF ( ALLOCATED( buildings(nb)%t_in ) ) &
1493	buildings(nb)%t_in = t_in_recv
1494	#else
1495	IF ( ALLOCATED( buildings(nb)%t_in ) ) &
1496	buildings(nb)%t_in = buildings(nb)%t_in_l
1497	#endif
1498	IF ( ALLOCATED( buildings(nb)%t_in ) ) &
1499	buildings(nb)%t_in = buildings(nb)%t_in / &
1500	( buildings(nb)%num_facade_h + &
1501	buildings(nb)%num_facade_v )
1502	!
1503	!-- Deallocate dummy arrays
1504	DEALLOCATE( t_in_l_send )
1505	DEALLOCATE( t_in_recv )
1506
1507	ENDDO
1508
1509
1510	END SUBROUTINE im_main_heatcool
1511
1512	!-----------------------------------------------------------------------------!
1513	! Description:
1514	!-------------
1515	!> Check data output for plant canopy model
1516	!-----------------------------------------------------------------------------!
1517	SUBROUTINE im_check_data_output( var, unit )
1518
1519	IMPLICIT NONE
1520
1521	CHARACTER (LEN=*) :: unit !<
1522	CHARACTER (LEN=*) :: var !<
1523
1524	SELECT CASE ( TRIM( var ) )
1525
1526
1527	CASE ( 'im_hf_roof')
1528	unit = 'W m-2'
1529
1530	CASE ( 'im_hf_wall_win' )
1531	unit = 'W m-2'
1532
1533	CASE ( 'im_hf_wall_win_waste' )
1534	unit = 'W m-2'
1535
1536	CASE ( 'im_hf_roof_waste' )
1537	unit = 'W m-2'
1538
1539	CASE ( 'im_t_indoor' )
1540	unit = 'K'
1541
1542	CASE DEFAULT
1543	unit = 'illegal'
1544
1545	END SELECT
1546
1547	END SUBROUTINE
1548
1549
1550	!-----------------------------------------------------------------------------!
1551	! Description:
1552	!-------------
1553	!> Check parameters routine for plant canopy model
1554	!-----------------------------------------------------------------------------!
1555	SUBROUTINE im_check_parameters
1556
1557	!!!! USE control_parameters,
1558	!!!! ONLY: message_string
1559
1560	IMPLICIT NONE
1561
1562	END SUBROUTINE im_check_parameters
1563
1564	!-----------------------------------------------------------------------------!
1565	! Description:
1566	!-------------
1567	!> Subroutine defining appropriate grid for netcdf variables.
1568	!> It is called from subroutine netcdf.
1569	!-----------------------------------------------------------------------------!
1570	SUBROUTINE im_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1571
1572	IMPLICIT NONE
1573
1574	CHARACTER (LEN=*), INTENT(IN) :: var
1575	LOGICAL, INTENT(OUT) :: found
1576	CHARACTER (LEN=*), INTENT(OUT) :: grid_x
1577	CHARACTER (LEN=*), INTENT(OUT) :: grid_y
1578	CHARACTER (LEN=*), INTENT(OUT) :: grid_z
1579
1580	found = .TRUE.
1581
1582	!
1583	!-- Check for the grid
1584	SELECT CASE ( TRIM( var ) )
1585
1586	CASE ( 'im_hf_roof', 'im_hf_roof_waste' )
1587	grid_x = 'x'
1588	grid_y = 'y'
1589	grid_z = 'zw'
1590	!
1591	!-- Heat fluxes at vertical walls are actually defined on stagged grid, i.e. xu, yv.
1592	CASE ( 'im_hf_wall_win', 'im_hf_wall_win_waste' )
1593	grid_x = 'x'
1594	grid_y = 'y'
1595	grid_z = 'zu'
1596
1597	CASE ( 'im_t_indoor' )
1598	grid_x = 'x'
1599	grid_y = 'y'
1600	grid_z = 'zw'
1601
1602	CASE DEFAULT
1603	found = .FALSE.
1604	grid_x = 'none'
1605	grid_y = 'none'
1606	grid_z = 'none'
1607	END SELECT
1608
1609	END SUBROUTINE im_define_netcdf_grid
1610
1611	!------------------------------------------------------------------------------!
1612	! Description:
1613	! ------------
1614	!> Subroutine defining 3D output variables
1615	!------------------------------------------------------------------------------!
1616	SUBROUTINE im_data_output_3d( av, variable, found, local_pf, fill_value, &
1617	nzb_do, nzt_do )
1618
1619	USE indices
1620
1621	USE kinds
1622
1623	IMPLICIT NONE
1624
1625	CHARACTER (LEN=*) :: variable !<
1626
1627	INTEGER(iwp) :: av !<
1628	INTEGER(iwp) :: i !<
1629	INTEGER(iwp) :: j !<
1630	INTEGER(iwp) :: k !<
1631	INTEGER(iwp) :: l !<
1632	INTEGER(iwp) :: m !<
1633	INTEGER(iwp) :: nb !< index of the building in the building data structure
1634	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1635	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1636
1637
1638	LOGICAL :: found !<
1639
1640	REAL(wp), INTENT(IN) :: fill_value !< value for the _FillValue attribute
1641
1642	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
1643
1644	local_pf = fill_value
1645
1646	found = .TRUE.
1647
1648	SELECT CASE ( TRIM( variable ) )
1649	!
1650	!-- Output of indoor temperature. All grid points within the building are
1651	!-- filled with values, while atmospheric grid points are set to _FillValues.
1652	CASE ( 'im_t_indoor' )
1653	IF ( av == 0 ) THEN
1654	DO i = nxl, nxr
1655	DO j = nys, nyn
1656	IF ( building_id_f%var(j,i) /= building_id_f%fill ) THEN
1657	!
1658	!-- Determine index of the building within the building data structure.
1659	nb = MINLOC( ABS( buildings(:)%id - building_id_f%var(j,i) ), &
1660	DIM = 1 )
1661	IF ( buildings(nb)%on_pe ) THEN
1662	!
1663	!-- Write mean building temperature onto output array. Please note,
1664	!-- in contrast to many other loops in the output, the vertical
1665	!-- bounds are determined by the lowest and hightest vertical index
1666	!-- occupied by the building.
1667	DO k = buildings(nb)%kb_min, buildings(nb)%kb_max
1668	local_pf(i,j,k) = buildings(nb)%t_in(k)
1669	ENDDO
1670	ENDIF
1671	ENDIF
1672	ENDDO
1673	ENDDO
1674	ENDIF
1675
1676	CASE ( 'im_hf_roof' )
1677	IF ( av == 0 ) THEN
1678	DO m = 1, surf_usm_h%ns
1679	i = surf_usm_h%i(m) !+ surf_usm_h%ioff
1680	j = surf_usm_h%j(m) !+ surf_usm_h%joff
1681	k = surf_usm_h%k(m) !+ surf_usm_h%koff
1682	local_pf(i,j,k) = surf_usm_h%iwghf_eb(m)
1683	ENDDO
1684	ENDIF
1685
1686	CASE ( 'im_hf_roof_waste' )
1687	IF ( av == 0 ) THEN
1688	DO m = 1, surf_usm_h%ns
1689	i = surf_usm_h%i(m) !+ surf_usm_h%ioff
1690	j = surf_usm_h%j(m) !+ surf_usm_h%joff
1691	k = surf_usm_h%k(m) !+ surf_usm_h%koff
1692	local_pf(i,j,k) = surf_usm_h%waste_heat(m)
1693	ENDDO
1694	ENDIF
1695
1696	CASE ( 'im_hf_wall_win' )
1697	IF ( av == 0 ) THEN
1698	DO l = 0, 3
1699	DO m = 1, surf_usm_v(l)%ns
1700	i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
1701	j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
1702	k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
1703	local_pf(i,j,k) = surf_usm_v(l)%iwghf_eb(m)
1704	ENDDO
1705	ENDDO
1706	ENDIF
1707
1708	CASE ( 'im_hf_wall_win_waste' )
1709	IF ( av == 0 ) THEN
1710	DO l = 0, 3
1711	DO m = 1, surf_usm_v(l)%ns
1712	i = surf_usm_v(l)%i(m) !+ surf_usm_v(l)%ioff
1713	j = surf_usm_v(l)%j(m) !+ surf_usm_v(l)%joff
1714	k = surf_usm_v(l)%k(m) !+ surf_usm_v(l)%koff
1715	local_pf(i,j,k) = surf_usm_v(l)%waste_heat(m)
1716	ENDDO
1717	ENDDO
1718	ENDIF
1719
1720	CASE DEFAULT
1721	found = .FALSE.
1722
1723	END SELECT
1724
1725	END SUBROUTINE im_data_output_3d
1726	!------------------------------------------------------------------------------!
1727	! Description:
1728	! ------------
1729	!> Parin for &indoor_parameters for indoor model
1730	!------------------------------------------------------------------------------!
1731	SUBROUTINE im_parin
1732
1733	USE control_parameters, &
1734	ONLY: indoor_model
1735
1736	IMPLICIT NONE
1737
1738	CHARACTER (LEN=80) :: line !< string containing current line of file PARIN
1739
1740	NAMELIST /indoor_parameters/ dt_indoor, initial_indoor_temperature
1741
1742	!
1743	!-- Try to find indoor model package
1744	REWIND ( 11 )
1745	line = ' '
1746	DO WHILE ( INDEX( line, '&indoor_parameters' ) == 0 )
1747	READ ( 11, '(A)', END=10 ) line
1748	ENDDO
1749	BACKSPACE ( 11 )
1750
1751	!
1752	!-- Read user-defined namelist
1753	READ ( 11, indoor_parameters )
1754	!
1755	!-- Set flag that indicates that the indoor model is switched on
1756	indoor_model = .TRUE.
1757
1758	!
1759	!-- Activate spinup (maybe later
1760	! IF ( spinup_time > 0.0_wp ) THEN
1761	! coupling_start_time = spinup_time
1762	! end_time = end_time + spinup_time
1763	! IF ( spinup_pt_mean == 9999999.9_wp ) THEN
1764	! spinup_pt_mean = pt_surface
1765	! ENDIF
1766	! spinup = .TRUE.
1767	! ENDIF
1768
1769	10 CONTINUE
1770
1771	END SUBROUTINE im_parin
1772
1773
1774	END MODULE indoor_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |