Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |