Home

Context Navigation

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

Last change on this file since 4083 was 3885, checked in by kanani, 6 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 85.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |