Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |