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

Last change on this file since 4209 was 4209, checked in by suehring, 2 years ago

Bugfix in initialization of indoor temperature; Prescibe default indoor temperature in case it is not given in the namelist input

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