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

Last change on this file since 4214 was 4214, checked in by suehring, 23 months ago

Bugfix, missing initialization and clearing of soil-moisture tendency (J.Resler)

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