NBSIR 88-3814
Progress Toward A General Analytical
Method for Predicting Indoor Air
Pollution in Buildings
Indoor Air Quality Modeling
Phase III Report
James Axley
U.S. DEPARTMENT OF COMMERCE
National Bureau of Standards
National Engineering Laboratory
Center for Building Technology
Building Environment Division
Gaithersburg, MD 20899
July 1988
75 Years Stimulating America's Progress
1913-1988
U.S. Environmental Protection Agency
U.S. Department of Energy
U.S. Consumer Products Safety Commission
-------
NBSIR 88-3814
PROGRESS TOWARD A GENERAL ANALYTICAL
METHOD FOR PREDICTING INDOOR AIR
POLLUTION IN BUILDINGS
INDOOR AIR QUALITY MODELING
PHASE III REPORT
James Axley
U.S. DEPARTMENT OF COMMERCE
National Bureau of Standards
National Engineering Laboratory
Center for Building Technology
Building Environment Division
Gaithersburg, MD 20899
July 1988
Prepared for:
U.S. Environmental Protection Agency
U.S. Department of Energy
U.S. Consumer Products Safety Commission
U.S. DEPARTMENT OF COMMERCE, C. William Verity, Secretary
NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
-------
Indoor Air Quality Model ABSTRACT
Phase III Report
ABSTRACT
This interim report presents the results of Phase III of the NBS General Indoor Air
Pollution Concentration Model Project. It describes;
a) a general element-assembly formulation of multi-zone contaminant dispersal
analysis theory that provides a general framework for the development of detailed
(element) models of mass transport phenomena that may affect contaminant
dispersal in buildings,
b) an approach to modeling the dispersal of interactive contaminants involving
contaminant mass transport phenomena governed by basic principals of kinetics
and introduces a linear first order kinetics element to achieve this end,
c) an approach to modeling the details of contaminant dispersal driven by
convection-diffusion processes in one-dimensional flow situations (e.g., HVAC
ductwork) and introduces a convection-diffusion flow element to achieve this
end,
and
d) the features and use of CONTAM87, a program that provides a computational
implementation of the theory and methods discussed.
The theory and methods presented are based upon a generalization of the building
idealization employed earlier [Axley, 1987]. Here, building air flow systems are
idealized as assemblages of mass transport elements, rather than simply flow elements
as used previously, connected to discrete system nodes corresponding to well-mixed
air zones within the building and its HVAC system. Equations governing contaminant
dispersal in the whole building air flow system due to air flow and reaction or sorption
mass transport phenomena are formulated by assembling element equations, that
characterize a specific instance of mass transport in the building air flow system, in such
a manner that the fundamental requirement of conservation of mass is satisfied in each
zone.
KEY WORDS: building simulation, indoor air quality, contaminant dispersal analysis,
element assembly, discrete modeling techniques
-------
Indoor Air Quality Model ACKNOWLEDGEMENTS
Phase III Report
ACKNOWLEDGEMENT
This work was supported by the following Interagency Agreements:
IAG: DE AI01 -86-CD 21013 Amendment Number 4 with the Department of Energy,
IAG: DW 1391103-01-2 with the Environmental Protection Agency, and
IAG: 74-25 Task Number 87-4 with the Consumer Products Safety Commission.
iii
-------
Indoor Air Quality Model CONTENTS
Phase III Report
CONTENTS
ABSTRACT ii
ACKNOWLEDGEMENTS iii
CONTENTS iv
PREFACE vi
NOTATION viii
PART I - THEORY 1-1
1. Introduction 1-1
1.1 Indoor Air Quality Analysis 1-1
1.2 The Well-Mixed Microscopic Model 1-4
2. Contaminant Dispersal Analysis 2-1
2.1 Element Equations 2-1
2.2 System Equations 2-5
2.3 Solution of System Equations 2-9
3. Interactive Contaminant Dispersal Analysis 3-1
3.1 Multiple, Noninteractive, Contaminant Dispersal 3-1
3.2 Basic Concepts of Reaction Kinetics 3-4
3.3 Kinetics Element Equations 3-10
4. One Dimensional Convection-Diffusion Flow 4-1
4.1 Convection-Diffusion Equation 4-1
4.2 Convection-Diffusion Element Equations 4-6
4.3 Use of the Convection-Diffusion Flow Element 4-8
4.4 Analytical Properties of the Convection-Diffusion ... 4-14
Element Equations
4.5 Comparison to Tanks-in-Series Idealizations 4-15
PART II - CONTAM87 USERS MANUAL 5-1
5. General Instructions 5-1
6. Command Conventions 6-1
7. Introductory Example 7-1
iv
-------
Indoor Air Quality Model CONTENTS
Phase III Report
8. Command Reference 8-1
8.1 Intrinsic Commands 8-1
8.1.1 HELP 8-1
8.1.2 ECHO 8-1
8.1.3 LIST 8-1
8.1.4 PRINT A= 8-1
8.1.5 DIAGRAM A= 8-2
8.1.6 SUBMIT F= 8-2
8.1.7 PAUSE 8-2
8.1.8 RETURN 8-2
8.1.9 QUIT 8-2
8.2 CONTAM87 Commands 8-3
8.2.1 FLOWSYS 8-3
8.2.2 FLOWELEM 8-5
8.2.3 KINELEM 8-7
8.2.4 FORM-[W] 8-9
8.2.5 STEADY 8-10
8.2.6 TIMECONS 8-10
8.2.7 Dynamic Analysis 8-11
8.2.7.1 FLOWDAT 8-11
8.2.7.2 EXCITDAT 8-13
8.2.7.3 DYNAMIC 8-13
8.2.7.4 Dynamic Analysis Example 8-15
8.2.8 RESET 8-15
9. Example Applications 9-1
9.1 IBR Test House Study 9-1
9.2 Carnegie-Mellon Townhouse Study 9-5
9.3 NBS Office Building Study 9-10
10. Summary and Directions of Future Work 10-1
REFERENCES Ref-1
APPENDIX - CONTAM87 FORTRAN 77 Source Code . . . Append-1
-------
Indoor Air Quality Model PREFACE
Phase III Report
PREFACE
The work reported here is a product of the General Indoor Air Pollution
Concentration Model Project initiated in 1985 at the National Bureau of Standards and
supported by the U. S. Environmental Protection Agency, the U.S. Department of
Energy, and the Consumer Products Safety Commission. The fundamental objective of
this project is to develop a comprehensive validated computer model to simulate
dynamic pollutant movement and concentration variation in buildings. The scope of the
project is ambitious; a full-scale, multi-zone building contaminant dispersal model that
simulates flow processes (e.g., infiltration, dilution, & exfiltration) and contaminant
generation, reaction, and removal processes is being developed.
During the planning stage of this project it was decided to organize efforts into three
distinct phases:
Phase I: formulation of a general framework for the development of general indoor
air quality analysis models (see [McNall et.al., 1986] for report of Phase I
work),
Phase II: development of a residential-scale model, based on the simplifying
assumption that air is well-mixed within each building zone, providing
simple simulation of HVAC system interaction, and
Phase III: extension of modeling capabilities to allow more complete simulation of
HVAC system interaction and consideration of rooms that are not well-
mixed.
This report presents analytical methods that, together with those methods developed
during Phase II of the project, satisfies the scope and objectives set for Phase III of the
"General Indoor Air Pollution Concentration Model" Project and, as such, completes
Phase III efforts. The report is organized in two parts. In the first part of the report the
underlying theory is presented;
Section 1: outlines the general aspects of indoor air quality analysis - making the
distinction between contaminant dispersal analysis, inverse contaminant
dispersal analysis, and air flow analysis - that the project has attempted
to address, and defines the approach taken to modeling,
Section 2: presents a general formulation of multi-zone contaminant dispersal
theory, using an element assembly approach,
Section 3: applies the theory from Section 2 to develop an interactive, multiple-
contaminant dispersal analysis method based upon the formulation of a
VI
-------
Indoor Air Quality Model PREFACE
Phase III Report
kinetics element designed to model mass transport phenomena governed
by the principals of kinetics,
Section 4: applies the theory from Section 2 to model the details of one-dimensional
contaminant dispersal driven by combined convection-diffusion mass
transport processes,
The second part of the report presents the practical implementation of the
contaminant dispersal analysis theory in the program CONTAM87;
Sections 5 -8: provide a users manual for the program CONTAM87, and
Section 9: gives examples of application of CONTAM87 to representative problems
of contaminant dispersal analysis.
The last section, Section 10, provides a summary of the work reported here and
outlines possible directions of future study.
The complete source code for CONTAM87 is listed in the appendix.
VII
-------
Indoor Air Quality Mode! PREFACE
Phase III Report
NOTATION
Scalars
Scalar variables will be designated by lower and upper case, plain, english and
greek characters. An equals sign enclosed in square brackets, [=], is used to designate
typical units of a variable. The more commonly used scalars are listed below:
C concentration in terms of mass fraction (mass-species/mass-air)
[=] Ib-species/lb-air or kg-species/kg-air
G species generation rate (mass/unit time)
[=] Ib-species/hr or kg-species/s
M mass of the volume of air within a given zone (mass-air)
[=] Ib-air or kg-air
P pressure (force/unit area)
[=] Ib/ft2 or Pascals (Pa)
t time
T temperature
[=] °F or °C
v velocity (e.g., of a fluid particle)
[=] ft/hr or m/s
w mass transport rate (e.g., due to flow, chemical reaction, etc.)
[=] Ib/hr or kg/s
x.y.z spacial coordinates
K reaction rate coefficient
[=] s-1
p density (mass/unit volume)
[=] Ib/ft3 or kg/m3
Indices
The indicial notation used in this report is modeled after the conventions that are
commonly used in structural analysis and Finite Element analysis literature. A variety of
indices may be associated with any single variable including pre-subscripts, pre-
superscripts, post-subscripts, and post-superscripts. Although the meaning of any
index should be clear from the context of the discussion, the conventions diagrammed
below will be followed to help maintain clarity.
viii
-------
Indoor Air Quality Model PREFACE
Phase III Report
species index -> a a <- element index
X <- variable
descriptive index -» con j <- node index
(In some contexts the post-subscript is used to indicate an element of an array as well.)
In addition to these indices which serve to specifically identify a variable, it will be
necessary to use additional indices to indicate the value (or approximation to the value)
of a variable at a discrete point in time or for a discrete step in an iterative scheme. In
both cases, additional superscripts and subscripts may have to be introduced. To
distinguish these time-step or iterative-step indices they shall be enclosed in
parenthesis as indicated below.
g(k) «_ the "kth" iterate
X.
con '(n + 1)<- the "n+1th" time step
Generally, the following conventions will be used for superscripts and subscripts:
a,b,c,... specific element indices
e general element index
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 1. Introduction
PART I - THEORY
The first section of this part of the report gives definition to the meaning of indoor air
quality analysis and describes the modeling approach taken to develop practical indoor
analysis tools. It will be seen that indoor air quality analysis involves (forward)
contaminant dispersal analysis, inverse contaminant dispersal analysis, air flow
analysis, and thermal analysis. The following three sections extend the contaminant
dispersal analysis theory developed during Phase II of the present project [Axley, 1987]
by first presenting a more general formulation of multi-zone contaminant dispersal
analysis theory and then applying this more general theory to a) dispersal problems
involving interactions between contaminant species and/or the building fabric and b)
dispersal problems where the details of convection-diffusion flow processes are
important (e.g., HVAC ductwork). Current research efforts focused on the inverse
contaminant dispersal analysis problem have led to promising new multi-zone tracer gas
techniques, the Pulse Tracer Techniques and have provided a better understanding of
existing tracer gas techniques. Formulations of building air flow and thermal analysis
theories that are compatible with the formulations of the forward and inverse contaminant
dispersal analysis theories have been presented earlier [Axley 1987; Axley 1985] and
will become the focus of future work.
1. Introduction
During the past decade, indoor air pollution emerged as an. international health
issue and, as a result, a new field of simulation, indoor air quality analysis, is emerging
to provide the means to predict concentration variation of indoor air contaminants in
existing and proposed buildings and, thereby, to assess the nature and severity of
potential indoor air quality problems. It may be expected that this new field will come to
play a key role in the development of strategies to mitigate indoor air quality problems
and, eventually, become central to the design of high quality indoor air environments.
1.1 Indoor Air Quality Analysis
The central concern of indoor air quality analysis is the prediction of airborne
contaminant dispersal in buildings. Airborne contaminants disperse throughout
buildings in a complex manner that depends on the nature of air movement in-to, out-of,
and within the building system; the influence of the heating, ventilating, and air
conditioning (HVAC) systems; the possibility of removal,, by filtration, or contribution, by
generation, of contaminants; and the possibility of chemical reaction, radio-chemical
decay, settling, or sorption of contaminants. In indoor air quality analysis we seek to
comprehensively model all of these phenomena.
More precisely, in indoor air quality analysis we consider building air flow systems
to be three dimensional fields within which we seek to completely describe the state of
infinitesimal air parcels. The state of such an air parcel is defined by its temperature,
1-1
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
1. Introduction
pressure, velocity, and contaminant concentration(s) - the state variables of indoor air
quality analysis.
Air Parcel State Variables
Temperature: T(x,y,Z,t)
Pressure: P(x,y,z,t)
Flow Velocity: v( X ,y,Z ,t)
Concentration: aC(x,y,Z ,t), PC(x,y,Z ,t v
ttj1". :•• • >'
' .-" ' for species a , (5,...
Figure 1-1 Indoor Air Quality State Variables
The central problem of indoor air quality analysis is, then, the determination of the
spacial (x.y.z) and temporal (t) variation of contaminant species concentrations or
contaminant dispersal analysis.
For a single noninteractive^ species, a, contaminant dispersal is driven by the air
velocity field and its variation with time and thus the contaminant dispersal analysis
problem, for this case, may be represented as:
Noninteractive Contaminant Dispersal Analysis
«C(x,y,z,t) = «C( v(x,y,z,t),...)
where the ellipses are used to indicate initial and boundary conditions required to
complete the definition of the analytical problem. To solve the contaminant dispersal
problem, then, the flow field must be either specified or determined.
Two approaches to flow determination exist. In the first approach a nonlinear flow
analysis problem and, in general, a coupled thermal analysis problem is formulated and
solved, given the environmental excitation (e.g., wind, solar, and thermal excitation)
acting on the building system. Alternatively, for existing buildings it may be possible to
"measure" building air flows using tracer gas techniques. These techniques are based
on the formulation and solution of the inverse contaminant dispersal analysis problem.
Functionally, these related problems take the following forms:
1 Noninteractive Contaminant: a contaminant whose dispersal is not affected by kinetics of
reaction, sorption, settling, or other similar or related mass transport phenomena.
1-2
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 1. Introduction
Coupled Flow/Thermal Analysis Inverse Contaminant Dispersal A nalysis
v(x,y,z,t) = v( P(x,y,z,t),...) '
P(x,y,z,t) = P(T(x,y,z,t),...)
T(x,y,z,t) = T( v(x,y,z,t),...)
Flow Analysis V(X,y,Z,t) = V( aC(X,y,Z,t), ...)
Buoyancy Effects (the basis of tracer gas techniqu es)
Thermal Analysis
When contaminant reaction, settling, sorption, etc. kinetics is important, the
contaminant dispersal analysis problem becomes a coupled (and, generally, nonlinear)
analysis problem as (the rate of change of) each species' concentration will depend
upon both species' concentrations and the air flow velocity field:
Interactive Contaminant Dispersal Analysis
«C(x,y,z,t) = «C( v(x,y,z,t), «C(x,y,z,t), PC(x,y,z,t), ... )
For such cases we say the contaminant is an interactive contaminant and describe the
analytical problem as a problem of interactive contaminant dispersal analysis.
A complete indoor air quality analysis package should provide the analyst with tools
to consider this relatively complex set of analytical problems related to the central task
of contaminant dispersal analysis. As indicated in Figure 1-2 one may anticipate three
basic indoor air quality analysis scenarios;
1 ) in some instances the analyst may choose to simply specify the flow field (e.g., in
design situations or in those cases where the HVAC system substantially
determines air flow in the building system) and directly consider the contaminant
dispersal analysis problem for specific indoor air pollutant sources or sinks,
2) for existing buildings, tracer gas techniques, based upon inverse contaminant
dispersal analysis methods, may be used to determine airflows that may then be
used to complete the required contaminant dispersal analysis for any number of
specific indoor air pollutant sources or sinks, or
3) in some instances the analyst may choose to complete an airflow analysis of the
building system, given building and wind characteristics, to determine the
airflows needed to complete the contaminant dispersal analysis task.
Many specific pollutant source or sink models will involve chemical or mass transport
governed by the kinetics of the mass transport phenomena; thus analytical tools are
needed to properly account for this. Finally, when airflow analysis is elected the
analyst will either have to specify the temperature field or determine it by solving the
coupled flow-thermal analysis problem to properly account for buoyancy effects; a
complete indoor air quality analysis package should provide this capability.
1-3
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
1. Introduction
or
Flow
Field
Specified
Kinetics
Inverse
Contaminant
Dispersal Anal.
Contaminant
Disperal
Analysis
or
or
Flow Analysis
Pressure B.C.s
Specified
Buoyancy Effects
Temperature
Field
Specified
or
Thermal
Analysis
Figure 1-2 The Central and Related Problems of Indoor Air Quality Analysis
1.2 The Well-Mixed Macroscopic Model
To develop this needed indoor air quality analysis capability we follow the tradition
established by others in the field of indoor air quality analysis [Sinden, 1979, Sandberg,
1984, Walton, 1985] and model building air flow systems using a well-mixed zone
simplification of the macroscopic equations of motion (i.e., mass, momentum, and
energy balances for flow systems) that, in essence, transforms the indicated field
problems discussed above into spatially discrete, but temporally continuous, ordinary
differential equations. The present approach breaks from this tradition, however, in that
an element assembly approach is taken to formulate the respective analytical problems.
That is to say:
building air flow systems (fields) will be idealized as assemblages of discrete
flow elements, that are used to model specific flow transport processes between
well-mixed zones, and kinetics elements, that are used to model specific
transport processes that occur within the well-mixed zones that may be
described using the principals of kinetics.
Such idealizations of building air flow systems may be represented graphically in a
1-4
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
1. Introduction
direct and intuitive way as illustrated in Figure 1-3 for a hypothetical building system.
0 —Zone Node
2 — Node Number
Actual Building
-I
Air-Flow System Idealization
Flow Element—
Element Flow-
Fiure 1-3 Idealization of A Hothetical Building Air Flow Sstem
With a knowledge of the air flow paths in the building system the analyst selects from
the library of available air flow elements to assemble graphically, and, hence
mathematically, the building air-flow idealization. Kinetics elements may then be added
to this assemblage to account for the nonflow transport processes that may occur within
a given zone. Thus, for example, the idealization presented above would, conceivably,
be appropriate for the analysis of carbon monoxide dispersal. The indicated flow
elements would model HVAC, infiltration, and exfiltration flow paths and the single
kinetics element (labeled Rx ) would model the kinetics of carbon monoxide generation
within the furnace system. Note that in this case a well-mixed zone is associated with
the furnace, a junction of the HVAC ducts, the exterior environment and each of the
rooms of the building system.
Presently, the library of flow elements contains those indicated in Figure 1-4. The
kinetics element and the convection-diffusion element are presented in the next section
of this report; the other elements were presented earlier [Axley 1987].
Library of Contaminant Dispersal Elements
<•>• 1 I-l
Element
— !]""=> Flow w/Filter
" Element
i!. . . .- Convect. -Diftus.
Flow Element
(p^ Kinetics
^-^ Element
Library of Flow Elements
Flow Resistance
Element
Fan/Pump
Element
Figure 1-4 Current Library of Indoor Air Quality Analysis Elements
1-5
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 1. Introduction
With each well-mixed zone we associate a set of discrete state variables with a
distinct, but arbitrary point in the zone, the zone node. These discrete variables are
meant to approximate the corresponding field variables in the zone at that point. For a
system idealized as n well-mixed zones, then, the key discrete state variables would
include:
{P} s {PL Pa... Pn} . the vector of system pressure variables (1.1)
("0 = fl"i. T2. - Tn} : the vector of system temperature variables (1.2)
and the vector of system concentration variables defined as:
• for the dispersal of a single species, a:
/r*\ = / <-f* &f* Qf* \
l / \ i * t> ""* n/ i
• for the dispersal of two species, a and (i:
^
r^i r G/^ B/^ O/^ B/"* Ot^ IV"* i
*'^t»'^tl"*"ll''lirMJ I
• etc.
where the subscripts are zone/node indices. These variables will be referred to as the
system (state) variables.
With each element "e" in the system assembly we associate one or more element
nodes. With each node we associate variables that define the state of the element - the
element (state) variables, which will normally be subsets of the system variables2,
and note their association with the system variables. Thus, for example, a contaminant
dispersal element having three nodes, i, j, and k, would have the element state
variables;
• for the dispersal of a single species, a:
ft a A T
- / cy*. *V"». ot/^ \
for the dispersal of two species, a and (3:
.cue p«e cue p_e cue p_e T
= ( Uj, Oj, Oj, Uj, Ok, U|
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 1. Introduction
general.
With these element variables in hand, element equations are formulated that
describe the specific mass and/or energy transport phenomena that the element is
meant to represent and, by demanding conservation of mass or energy transport at each
of the system nodes, these element equations are then assembled to form system
equations governing the behavior of the building air flow system as a whole.
From a practical point of view, the element assembly approach is intuitively
satisfying and allows consideration of systems of arbitrary complexity. From a research
and development point of view this approach separates the general problem of indoor
air quality analysis into two primary subproblems; element development and
development of solution method. Research efforts can, thus, focus on the modeling of
specific transport phenomena to develop improved or new elements or, alternatively,
focus on developing improved methods of solving the resulting equations while
accounting for the complex coupling that exists between the thermal, dispersal, and flow
analysis problems.
The approach has been formulated to be completely analogous and compatible with
approaches based upon Generalized Finite Element Method [Zienkiewicz, 1983]
solutions of the microscopic equation of motion for fluids and makes use of the numerical
methods and computational strategies that have been developed to support the Finite
Element and associated methods. It is expected that this compatibility will, eventually,
allow the analyst to employ mixed idealizations of building air flow systems wherein a
portion of the building air flow system is modeled in detail using microscopic elements
(e.g., elements based upon Finite Element approximations of the governing microscopic
continuum equations) while the rest of the air flow system is modeled using discrete
elements. In this way the analyst may study the details of dispersal in one area of the
system, accounting for whole-system interaction, without the overhead of modeling the
entire system microscopically. The one-dimensional convection-diffusion element
presented in the next section represents the first step in this direction.
1-7
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 2. Contaminant Dispersal Analysis
2. Contaminant Dispersal Analysis
Multi-zone building contaminant dispersal analysis theory has placed a singular
emphasis on contaminant dispersal driven by flow mass transport processes [ Sinden,
1979, Sandberg 1984, Walton 1985] even though it has long been recognized that the
dispersal of many important indoor air contaminants are affected by other mass transport
processes as well, most notably, processes associated with reaction, sorption, and
settling phenomena. The flow-element-assembly formulation of multi-zone building
contaminant dispersal analysis theory developed during Phase II of the current project
provides a conceptual framework to extend existing dispersal analysis theory to account
for these other mass transport processes. Extending the flow-element-assembly
approach, this section presents a general formulation of a multi-zone contaminant
dispersal analysis theory that provides a basis for the development of more complete
models of contaminant dispersal in buildings.
The general formulation of multi-zone contaminant dispersal theory is straightforward.
We first establish a restricted, but very general, form for equations that will be used to
describe mass transport at the element level1. Then, by establishing the
correspondence between the element concentration variables, {Ce}, and the system
concentration variables, {C}, and demanding species mass conservation at each of the
system nodes we show that these element equations may be assembled to form
equations governing the system as a whole. Consideration of boundary conditions, the
qualitative character of these equations, and the solution of these equations was
presented earlier [Axley 1987] and, therefore, will not be emphasized here.
2.1 Element Equations
As indicated above, it will be useful to distinguish those elements that model the
transport of mass from zone to zone by flow processes from those elements that model
the transport mass within a zone from species to species (e.g., by chemical or radio-
chemical reaction) or, possibly, from species to the environment of the zone, itself (e.g.,
by chemical or radio-chemical decay to a "noncontaminant" product, absorption,
adsorption, or settling processes). In either case, we shall attempt to describe the
behavior of an element by equations of the form:
(2.1)
where;
{w ) is a vector of species mass transport rates into the element
({C }) is a transformation of {C8} that has the form of a linear transformation
that is specific to a given class of elements
1 That is to say, to model specific instances of mass transport phenomena in a building's air flow
system.
2-1
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
2. Contaminant Dispersal Analysis
to!
is a vector of element derived species generation rates.
Element Mass Transport Rates The vector of species mass transport rates, {
for the dispersal of a single species a, may be represented diagrammatically as shown
below for a hypothetical three-node flow element and a single-node kinetics element,
where the arrows indicate positive mass transport rates. In the case of the flow element,
mass is transported physically
node k
flow element Q ; \
zone k
kinetics element e
J_
zonei rt|)
a e
CT
a e
Wj
•node i ""
Figure 2-1 Element Mass Transport Rate Variables
by the air flow moving from the zone into the element and, thus, the arrows used indicate
the sense or direction of the averaged or bulk fluid velocity - a common convention.
For the kinetics element mass transport is somewhat more subtle; mass may not literally
be transported out of the zone, rather mass of species a is, typically, converted from a
form that is considered to be a contaminant to another form (e.g., another compound or
phase) that is not of any special interest. Thus, for the kinetics element the arrow
indicating mass transport is directed into the "element" from the zone node to indicate
only that mass of species a is being removed from the zone by the element. It should be
noted that:
for each of the element concentration variables, cf, there exists a
corresponding element mass transport rate variables, wf.
This will also be true for the dispersal of more than one species.
For contaminant dispersal involving multiple species, then, a single physical flow
element might be thought to transport each individual species from zone-to-zone while a
kinetics element might be thought to transport mass, by conversion, from each of the
species to any or all of the other species and/or from any of the species to a form that
has no special interest, within the single zone associated with the kinetics element.
Extending this notion of an element, a combined flow and kinetics element, that not only
transports mass of one species from one zone to another, but transports, by conversion,
that species to another species or form during the flow passage, is also not only
2-2
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 2. Contaminant Dispersal Analysis
conceivable but reasonable for many reactive contaminants such as NC>2 and the radon
chain of decay products. (Inasmuch as it is difficult to represent these possible multi-
species mass transport/conversion phenomena diagrammatically we shall not attempt to
do so, here.)
Element Transformation Operator The element transformation operator L () is
restricted to the form of a linear transformation:
dt (2.2)
where;
[x8], [ye], [ze] are transformation coefficient matrices
[xe] is the element transport matrix
[y8] is the element mass matrix
but we admit transformation coefficient-matrices that may, in fact, vary with time and/or
depend, nonlinearly, on
[x6] = [x. ? jn general (2.3a)
[y°I=-[ye(t,{C^)] , in general (2.3b)
[z6] = [ ze( t, {C6}) ] > in general (2.3c)
etc.
thus a practically endless variety of element equations may be formulated that have this
form and, as such, the restriction to this form should not lead to any serious limitation.
Simple Flow Elements By assuming flow through a two-node flow element is
practically instantaneous and well-mixed, the mass transport of a single species, say a,
from element node i to j, due to an air mass flow rate we(t) from i to j, may be described
by the following element equations [Axley 1987]:
; we(t) > 0
(2.4a)
or, in this case we have:
(2.4b)
[x6] -rt = we(t) \ J] ; [y6] =[0] ; [z6) =[0] ; {ge} = {0}
2-3
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
2. Contaminant Dispersal Analysis
= {
a e,T
W,
(2.4C)
where we identify the element transport matrix for this case as the element mass flow
rate matrix, [f6]. It should be noted that the transformation matrix [xe] is seen to vary with
time to account for the time variation of flow through the element. (Figure 2-2, below,
should help to clarify the meaning of the element variables in this case.)
element "e"
Figure 2-2 Simple Two-Node Contaminant Dispersal Flow Element Variables
Simple Flow Element with Filter The simple flow element equations, above, may
be modified to account for the action of a filter that removes a fraction, rj, of the
contaminant as it passes through the element [Axley 1987], to yield the following
element equations;
a e
a e
Wj
- we(t)
0
0
(2.53)
or, in this case we have:
[x6! = [f6] = we(t)
1
(n-D
(2.5b)
where, again, we identify the element transport matrix as the element mass flow rate
matrix, [f6]. In this case the time variation of the first transformation matrix, [xe], could be
due to both the time variation of flow through the element and the time variation of the
filter efficiency, rj = t|(t). Furthermore, it should be recognized that the filter efficiency
will, in general, vary with each contaminant so that this first transformation matrix may be
expected to be species dependent. Following the notational convention established to
distinguish species types (i.e., the use of a leading superscript) we shall indicate this
species dependency as [axe], [Pxe] ,... for species a, p,... where:
2-4
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 2. Contaminant Dispersal Analysis
(Vi) o
0
(2.5C)
2.2 System Equations
Equations governing the dispersal of contaminants in the system as a whole may be
assembled from element equations of the form of equation (2.1) by first transforming the
element equations so that they are expressed in terms of system variables. To this end
we recognize that there exists a one-to-one correspondence between an element's
concentration variables, {C6}, and the system's concentration variables, {C}, that may
be described by a simple Boolean transformation as:
(2.6)
where;
B-
IB1 is an m x n Boolean Transformation matrix (i.e., consisting of only
ones and zeros) for an m-node element within an n-node system
idealization
The Boolean transformation is simply a means to express the equality of each of the
element concentrations variables with its associated system concentration variable
within the framework of concise vector notation; it defines the relation between the
(larger) vector of system concentration variables and the (smaller) vector of a specific
elements concentration variables, a subset of the system variables.
This same Boolean transformation matrix may be used to transform the vector
element mass transport rates, {w6}, into a "system-sized" vector of mass transport rates
for element "e", {W6}, as:
T
{W6} = [B6] {w8} (2.7)
This vector {We} will have the same number of elements as the system concentration
vector {C} , providing a correspondence between each system concentration variable
and a "system-sized" mass transport rate for the element "e". It represents the net
species mass transport rate from each of the system nodes into a specific element "e"
and, therefore, the sum of these mass transport vectors for all elements in the system
assemblage will equal a vector describing the total mass transport from the system
nodes into all elements combined.
2-5
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
2. Contaminant Dispersal Analysis
-ah
/ total species
e, I mass transport
{W } 55 / fromeachnode
1 into connected elements
I ateachnode
(2.8)
Demanding the conservation of species mass at each of the system nodes, the sum
of the quantity above plus the rate of change of species mass within each zone must be
equal to any species mass generated within the zone:
where;
(2.9)
the (diagonal) zone air volume mass matrix defined as (for n zones):
• for a single species (i.e., with {C} defined by equation (1.3a)):
M! 0 ... 0
0 M2 ... 0
0 0 ... Mn
for two species2 (i.e., with {C} defined by equation (1.3b)):
(2.10a)
"M!
0
0
0
0
0
0
M!
0
0
0
0
0
0
M2
0
0
0
0 ...
0 ...
0 ...
M2...
0 ...
0 ...
0
0
0
0
Mn
0
0
0
0
0
0
Mn_
(2.1 Ob)
• etc.
Mj = the mass of the air in the volume of zone i
[G] is the zone species generation rate vector, defined as
• for a single species, a:
,... "Gn}T
2 One could conceivably associate a different "active" zone air volume with each species and have, in
this case, a diagonal mass matrix of the form: diag ( aM-| , ^M-| , atJ\2 , ^2 , — aMn , ^Mn ,).
2-6
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 2. Contaminant Dispersal Analysis
• for two species, a and (J:
/r. . a_ p a- p_ a-. p_ T
{G} s { G! , G!, G2. Q2,... Gn> Gn/ (2.
• etc.
ot
Gj = the mass generation rate of species a in zone i
General Expression for Multi-Zone Dispersal Analysis Substituting the
transformation relations, equations (2.6) and (2.7), along with the general form of the
element equations, equation (2.1), into the species mass conservation relation,
equation (2.9), we obtain the final result — a general expression for multi-zone
contaminant dispersal analysis:
(2.12a)
where;
e =a< b> ;•• the system (mass) transport matrix (2.1 2b)
8 = a> b> - the system mass matrix (2. 1 2c)
[Z]=
e
etc.
{G} = {G}
the system generation vector (2.1 2e)
8 =a- b> -
It should be noted that in this general formulation the system mass matrix and system
generation vector have element contributions that add to the more familar zonal values.
The kinetics element, that will be presented in Section 2.2, will be seen to provide
element contributions to the system generation vector. The convection-diffusion
element, that will be presented in Section 2.3, will be seen to provide element
contributions to both the system mass matrix and the system generation vector.
2-7
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 2. Contaminant Dispersal Analysis
The Assembly Operator The summation and Boolean transformation of element
matrices, contained in the expressions above, is an operation that recurs frequently in
the Finite Element and related discrete modeling literature and, therefore, has come to
be defined as a standard operation - the assembly operation - designated by the
symbol A, the so-called assembly operator where;
A {v6} - I pflvl
e=a.b,... e=a.b.... for element vector assembly (2.1 3a)
and
A
e=a.b,... e=a.b.... for element matrix assembly (2.1 3b)
The assembly operator is therefore simply a generalization of the conventional
summation operator, S, and equal to the summation operator when the Boolean
transformation matrices equal the identity matrix.
The assembly operation is important, theoretically, in that it provides the necessary
formal definition of the assembly process required for subsequent mathematical
analysis. It does not, however, define an efficient numerical procedure for assembling
the element arrays needed to form the system equations for practical contaminant
dispersal analysis — the indicated Boolean transformations involve multiplications by
either zero or one and, therefore, need not be actually implemented. Practically, then,
the assembly operation is carried out using relatively simple algorithms that accumulate
element array terms within system array memory locations according to the physical
connectivity of each element. The "LM Algorithm" presented by Bathe [ Bathe 1982]
provides an example of one possible algorithm.
Relation to the Phase II Formulation Previously [Axley 1987], we considered
contaminant dispersal due only to flow mass transport via simple elements, with and
without filtration, and at that time the element transport matrices, [xe] , defined above by
equations (2.4b) and (2.5b) were identified as element mass flow rate matrices and
given the symbol [f8]. These element mass flow matrices were then assembled to form
the system transport matrix identified, then, as the system flow matrix and given the
symbol [F]. We shall see that, in the general case, the system transport matrix [W],
defined above, will be equal to the sum of the system flow matrix and what will be called
the system kinetics matrix, [K], assembled from element kinetics transport matrices, [k8],
as:
= [F] + [K] (2.i4a)
[F]= A
flow elements (2.14b)
2-8
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 2. Contaminant Dispersal Analysis
[K] = A
kinetics elem ents /2
2.3 Solution of System Equations
The contaminant dispersal elements developed to date are all described in terms of
first order linear transformations (i.e., involve only [xe] and [ye] transformation matrices)
having, in some cases, time varying element transport matrices (i.e., having [xe] =
[xe(t)], consequently system equations assembled from these element equations will be
limited to the first order form:
= (G}
where the system transport matrix will, in general, vary with time: [W] = [W(t)].
System equations of this form are identical, in form, to those system equations that
result from idealizations restricted to assemblages of simple flow elements. Therefore,
the procedures used to account for boundary conditions and possibilities of massless
nodes, the solution options that may be considered (i.e., steady state analysis,
eigenanalysis, and dynamic analysis), and the numerical methods that may be
employed to affect these solutions are identical to those discussed earlier [Axley 1987]
and will not be considered here. The qualitative character of equations (2.15) depends
critically upon the qualitative character of the element equations from which they are
assembled, therefore, the qualitative analysis presented earlier [Axley 1987] will have
to be reconsidered with the introduction of each new element.
2-9
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
3. Interactive Contaminant Dispersal Analysis
Often, indoor air quality analysis will involve consideration of several contaminants
and their dispersal in a building. Some of these contaminants may:
a) be absorbed or adsorbed by building materials or other contaminant particles,
b) settle from suspension or precipitate from (gaseous) solution, or
c) decay radiochemically, decompose chemically, or react with other
contaminants to produce product contaminants (or other substances that are of
no particular interest).
That is to say, contaminant dispersal processes may be complicated by the kinetics of:
a) sorption processes,
b) settling or precipitation processes, or
c) chemical or radio-chemical reaction processes
that must be accounted for.
In this section we shall introduce a general approach to extend noninteractive
contaminant dispersal theory to account for mass transport phenomena governed by the
kinetics of these processes, based upon the principals of reaction kinetics, to develop
an interactive contaminant dispersal theory. We shall set the stage by first considering
possible forms of system equations for multiple, noninteractive contaminant dispersal
analysis then, after a review of some basic concepts of reaction kinetics, go on to
develop the so-called kinetics element equations that will become the basis of the
interactive contaminant dispersal theory.
3.1 Multiple, Noninteractive, Contaminant Dispersal
The dispersal of each contaminant of a given set of noninteractive contaminants will
be governed by the single-species contaminant dispersal equation, thus the dispersal of
the noninteractive contaminant set may be represented by a set of equations of the form:
{aG}
{PG}
(3.1)
3-1
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
3. Interactive Contaminant Dispersal Analysis
where;
a, (5,... are species indices
I n is defined by equations (2.5c) and (2.13b)
[M] is by equations (2.10)
(Note that, in general, the flow matrices for each species may not be identical because
individual flow elements may act to filter contaminant species differently.)
Contaminant dispersal analysis for the set could, then, be computed by simply
completing a separate analysis for each species of interest. If, however, the system
characteristics change with time (e.g., airflow within the building is nonsteady and thus
the flow matrix, [F], changes with time) and the flow matrices for each species are
identical then it would be computationally more efficient to simultaneously compute the
response of each species while stepping through time as suggested by rewriting
equation (3.1) in the form:
[*F]KaC},{PC},...]
for;
d{aC} d{°C}
dt ' dt ""
-...•[>]
(3.2)
As an alternative approach, that will help set the stage for multiple interactive
contaminants, we may write the uncoupled set of equations given by equation (3.1) as
an expanded system of equations of the form:
={G}
(3.3a)
where system variables are organized by species for each node of the system
idealization as:
mi r
{«} = (
V n \
,-.. Gn, Gn,...}
(3.3b)
(3.3c)
The system flow transport matrix, [F], in this case, may be assembled by species and,
then, by element as:
[F] = A
e = a,b,...
A [V]
o=a, p,...
(3.3d)
where;
a is a general species index ( a, (3,... are specific species indices)
3-2
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
3. Interactive Contaminant Dispersal Analysis
e is a general element index (a, b,... are specific element indices)
The inner assembly sum, by species, may be thought to generate element equations
for a noninteractive, multi-species flow element. This multi-species flow element could,
then, be assembled in the usual manner to form the system equations. For three
contaminant species, a, p, y, the noninteractive, multi-species, flow element transport
matrix (with filtration of each species) would have the form:
A
o=a,P,
the noninteractive, multi-species, flow element transport matrix
(
1
0
0
VD
o (
0
0
1
0
0
V
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
(3.43)
for element flow from node i to node j and element concentration variables organized as:
e 8 8
YeT
This noninteractive. multi-species, element flow matrix is seen to be very sparse.
Consequently, assemblies of such elements would result in extremely sparse system
equations. It would, therefore, be computationally impractical to employ this approach
for noninteractive, multi-species, contaminant dispersal analysis — one should use the
strategies indicated by equations (3.1) or (3.2) instead. It will be seen, however, that
kinetic interactions among contaminant species will act to couple the species variables
at the system nodes (zones) with the result that system matrices will tend not only to be
much less sparse, but reasonably well-banded. The use of noninteractive, multi-
species, flow elements will, therefore, become attractive when combined with the
kinetics elements presented subsequently for interactive contaminant dispersal
analysis.
The program CONTAM87, documented in the second part of this report, organizes
system and element variables following equations (3.3b), (3.3c), and (3.4b) and makes
use of the double assembly process for simple flow elements defined by equations
(3.3d) and (3.4a).
3-3
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
3.2 Basic Concepts of Reaction Kinetics
Reaction kinetics involves the study of the rate of change of chemical components
m a single or related series of chemical reactions. Some basic concepts of (isothermal,
constant volume) gas reaction kinetics will be reviewed here; greater detail may be
found in one of several texts on the subject [Moore 1981, Nicholas 1976, Walas 1959].
Much of this material may also be applied to sorption, settling or precipitation, radio-
chemical, and other chemical phenomena that may be important for some interactive
contaminants in indoor air quality analysis.
A general form of a chemical reaction involving reactants, a, (J that react to form
products, p, a may be represented as;
catalyst
a + p + ...—> p +a +...
catalyst
where _> indicates the possible affect of catalysts on the reaction.
(3.5)
Given the rate of change of a selected component's concentration, say a, is defined
as;
d°C
aR = :^-r^- ; rate of reaction (\nterms of reactanta)
dt
(3.6)
where:
aC is the concentration of species a measured in terms of mass of a per
unit mass of air (i.e., strictly speaking the mass fraction of a)
and the stoichiometry of the reaction, expressed in terms of relative masses. am, Pm, ....
of reactants and products as;
catalyst
a • B „ p a
ma + mp + ... —> mp + ma + ..
(3.7)
the rate of change of the other components' concentrations may be related to that of the
selected component's as;
P.
m
m
m
m
°R
(3.8)
3-4
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
thus, the rate of a given chemical reaction may be described in terms of the rate of
change of concentration of any one of the reactants or products.
In general, the rate of a given chemical reaction may depend upon a variety of
factors including reactant and catalyst concentrations, temperature, T, pressure, P, and
the detailed mechanisms of the chemical reaction, therefore, rate expressions take the
general functional form of;
C, PC,... PC, "C,...T,P,...) (3.9)
Constant Rate Expressions In some instances the rate of reaction may remain more
or less constant:
= o (3.10a)
or depend solely on temperature and pressure:
R = R0(T, P) (3.1 Ob)
Examples include the catalytic decomposition of some gases, such as ammonia, the
radioactive decay of isotopes with very long half lives, such as Ra226, the controlled
burning of fossil fuels, and other relatively slow reactions driven by reactant and
product concentrations that remain, more or less, constant over the time of interest.
Power Law Rate Expressions In many cases the explicit form of a rate expression
will prove to be rather complex. In some cases, however, (empirical or semi-empirical)
rate expressions may be employed that take the form of so-called power expressions:
°R = aK(T,P)(aC)a(PC)b..(PC)rC)S... (3.11)
where:
«K = the rate constant [=] 1/time
a, b, ...r, s,... = constant exponents
Reactions governed by such power expression are classified by their overall order- the
sum of the constant exponents - or by their order with respect to each kinetically active
component - the constant exponent of that component. For the reaction described by
equation (3.11), then, the overall order will be (a + b + ... + r + s + ...) and the order with
respect to component a will be simply (a).
3-5
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
3. Interactive Contaminant Dispersal Analysis
First Order Rate Expressions Rate expressions for certain general classes of
reactions, including single-reactant, consecutive, opposing, and concurrent first order
reactions, often take the form of linear combinations of contaminant concentrations:
{R}=-[K]{C}
or
(3.12a)
act aB ao
K - K ... - K
PP
K
oa <
- K -
ao
(3.12b)
where we have included the constant component, {R0}, for completeness and
recognize that, again, the rate coefficient matrix, [K], and the constant component
vector, {R0}, will, in general, vary with temperature and pressure. It should be noted
that equation (3.12b) has been written so that all rate coefficients, V^K ; \|/, co = oc.p,... o ,
will be positive for realistic reactions.
In fact, it is possible, in principal, to linearize any given rate expression about some
(current) state of concentration, say (aC0 , PC0 , ...), by employing a Taylor's expansion
about that state, to obtain an approximate rate expression expressed as the sum of a
series of first order rate expressions, as:
aR(aC,PC,...)
3aR(cb0,lfc0l...)I
apc
(3.13)
that, together with equation (3.8), may be used to form a linearized system of first order
rate expressions of the form of equations (3.12). One could, conceivably, employ this
linearization strategy, within an appropriate nonlinear solution method, to account for
arbitrarily complex kinetics. The first .order kinetics element, that will be presented
subsequently, provides a first step in this direction.
Linear systems of first order reaction expressions are defined by the characteristics
of their reaction rate coefficient matrices, [K]. To gain a better understanding of the types
and characteristics of such reactions several specific classes of reactions are
described below.
3-6
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
3. Interactive Contaminant Dispersal Analysis
Single-Reactant First Order Reactions For reactions involving single
contaminant reactants that decompose or decay to form products (that are of little
particular interest):
a -» products
P -> products
products
the rate coefficient matrix takes the following form:
(3.14)
[K] =
•
act
K
0
0
0 ...
PPK ...
0 ...
0
0
oa
K
(3.15)
Consecutive First Order Reactions The radioactive decay chain of Radon gas is
an especially important example of a consecutive first order reaction series. The
reaction rate expression for a simple two-step consecutive reaction will be discussed
first then the general case will be considered.
For a two-step consecutive reaction series involving a single reactant at each step:
a -» p
P -» products
with reactions governed by rate expressions of the following form:
(3.16)
PaKaC- PPKPC
the matrix of rate coefficients becomes:
[K] =
act
K 0
(3.17)
(3.18)
Here, the generalization to a multi-step reaction involving single reactants at each
step is straightforward. For a general multi-step consecutive reaction series:
3-7
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
a -» p
p -* ...
p -» a
a -> products
governed by rate expressions of the form:
(3.19)
"R = -aaKaC
KK °C - PPK PC
°PKPC-
the matrix of rate coefficients becomes:
[K] =
aa «
K 0
0 0
0 0
... 0
... 0
pp
... K
0
0
0
oo
K
(3.20)
(3.21)
Opposing First Order Reactions For simple reversible reactions involving a single
reactant and single product:
a <-» p
governed by rate expressions of the form:
°R = -aaKaC+ aPKPC
(3.22)
(3.23)
the matrix of rate coefficients becomes:
3-8
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
3. Interactive Contaminant Dispersal Analysis
[K] =
act aB
K - K
PP
(3.24)
For the more general case of a series of reversible reactions involving single
reactants and products:
a <->P<-»y<-»8<-» ...
governed by rate expressions of the form:
(3.25)
°R = -
apKPC
paKac- PPKPC+ Vc
*KPC - VC + Y5K5C
the matrix of rate coefficients takes the tridiagonal form:
(3.26)
[K] =
•
""K -aPK 0
-P\ PPK -\
0 A Tc
0 0 -5YK
•
0 ...
0 ...
68K ...
_
(3.27)
Concurrent Linear Reactions Due to their linearity, rate coefficient matrices for
concurrent linear reactions may simply be added to obtain an effective rate coefficient
matrix for all reactions combined. For example, consider a set of concurrent reactions
involving a single contaminant reactant, a, that decays or decomposes to produce
products (3, y, a :
a -» p ; reaction "a"
a -» Y ; reaction "b"
a -» 8 ; reaction "c"
governed by rate expressions of the form:
(3.28)
3-9
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
o_ . aa aa aa . a_
H = -(a K+ b K + c K) C
'R =
8R =
The matrix of rate coefficients is:
[K] =
aa aa aa 1
a K + b K + c K)
Pa
"I*
"b ^
5a
C
0
0
0
0
0
0
0
0
0
0
0
0
(3.29)
(3.30a)
which is seen to be the sum of the individual reaction rate coefficient matrices:
aUK 0 0 0
-J°K 0 0 0
0000
0 000
• at
b°K 0 0 0
0 000
-b°K 000
0 000
c°K 0 0 0
0 000
0 000
Sa
-c K 000
(3.30b)
3.3 Kinetics Element Equations
The development of a general kinetics element is straightforward. Limiting our
considerations to mass transport phenomena occurring within a specific zone "i",
containing a set of contaminant species, a, (3, y,.... we first identify the relevant element
variables as:
- /
= I
ae
Ye
; the element state variables
(3.31)
and
a raepeve .T
} = t W| , w, , -W| ,...)
; the element mass flow rate vector (3.32)
Then, assuming that the mass transport phenomena to be modeled is governed by the
3-10
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
kinetics discussed above, a general kinetics element equation follows directly from the
definition of rate of reaction, equation (3.6) and the general form of rate expressions,
equation (3.9), as:
{w } = - [Mn {Rj ({Cl, T, P)} . tne generai kinetics element equation (3.33a)
where;
"" diO 0 ...
0 0 Mj...
(3.33b)
(3.33C)
(The superscript e has been added to identify the specific kinetics element and the
subscript i has been added to identify the specific node/zone being considered. The
negative sign is needed as species mass transport "into" the element (i.e., removed from
the zone) is defined to be positive.)
Using the notation introduced earlier (equations (2..1) and (2.2)) it is seen that the
general kinetics element is one defined by the following element arrays:
= 0 ; [y8] = 0 ; [z6] = 0 ; ... ; = .. (333d)
an element that is defined in terms of only element derived species generation rates.
If the rate expressions are constant rate expressions, then the element derived
generation rate terms will simply add to any direct species generation rates specified
within a zone, after equation (2.12e) as:
kinetics dements
(3.34)
where {Gj} and {Gj} are the subsets of the system vectors {G} and {G} corresponding to
node/zone i. In these cases there will be no practical difference between the physical
generation of species mass (e.g., by physical release of a contaminant) and the kinetics
generation of species mass (e.g., by chemical or physical-chemical processes) and,
3-11
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
therefore, the analyst may model either using simple noninteractive contaminant
dispersal analysis techniques.
The form of equations (3.33) is deceptively simple. The rate expressions defining
these element derived species generation rates depend on species concentration, in
general, so that the general kinetics element introduces a nonlinear species generation
contribution (i.e., a species generation rate that depends nonlinearly on the solution
vector {C}), which is distinctly different from the (constant or time dependent) nodal
direct generation contribution. The solution of the contaminant dispersal problem
involving general kinetics elements will, therefore, require the application of a nonlinear
solution strategy in the solution process. While this adds complexity to the analysis
process it should not be difficult to develop an appropriate nonlinear solution strategy
(e.g., using one or more of the strategies considered earlier to solve the nonlinear air
flow analysis problem [Axley 1987]) for certain classes of kinetics.
Few interactive indoor contaminants have been studied in sufficient detail to
completely define their kinetics, therefore, the development of nonlinear solution
techniques for arbitrary nonlinear kinetics would be premature at this time. Instead we
have limited our attention to kinetics described by linear systems of first order rate
expressions (that lead directly to linear systems of equations for interactive contaminant
dispersal analysis) to develop a practically useful interactive contaminant dispersal
analysis method. It seems likely, however, that more complex kinetics may be modeled
in the future by employing the combination of the linear method described below and the
Taylor's expansion presented earlier (equation (3.13)) to linearize rate expressions,
within an appropriate iterative solution scheme.
First Order Kinetics Element Equations For reaction kinetics described by
systems of first order equations, equations (3.12):
, (335)
the kinetics element equations (3.33) become:
{w8} =
oi (336a)
or:
[x6] = [M*HK°] ; [y6] ,= 0 ; [z6] = 0 ; ... ; {g6} = [M^R*.} (336b)
and, again, one must keep in mind that the rate coefficient matrix and constant rate
component will, in general, be temperature and pressure dependent:
3-12
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 3. Interactive Contaminant Dispersal Analysis
e e
[KJ] = [K^Lr/j (3.36c)
{R0j} = {R0i(T,P)} (336d)
It will be convenient to introduce a new variable for the linear first order element
kinetics transport matrix, as:
Q
[kj = [Mjj[Kj] : the element kinetics transport matrix (3.37)
and a corresponding variable for the system kinetics transport matrix, that is assembled
from the element kinetics transport matrices in the usual manner, as:
M • A [k6!
kinetics elements • tne system kinetics transport matrix (3.38)
so that the system transport matrix, [W], (equation (2.12b)) may be thought to equal the
sum of the familiar system flow matrix, [F], and the system kinetics transport matrix as
noted in equations (2.13) and repeated here:
= [F] + [K] = A [fl + A
flow elements kinetics elements /3 3g\
The program CONTAM87, presented in the second part of this report implements the
interactive contaminant dispersal theory, presented above, providing the linear first
order kinetics element in its library of elements and ordering system concentration
variables as discussed above (equation (3.3b)) to enable the proper assembly of flow
elements (i.e., by equation (3.3d)) for multi-contaminant dispersal analysis. Examples of
the application of these techniques are presented in Section 9.
3-13
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
4. One Dimensional Convection-Diffusion Flow
The flow element presented earlier [Axley 1987] provides the simplest modeling of
species mass flow from one zone to another. This simple element is based on the implicit
assumption that the volume and the length of the flow passage is negligible and,
therefore does not account for any dynamic dispersal phenomena occurring within the
flow passage.
For problems where zonal dynamics is of primary interest and it is suspected that
flow passage dynamics need not be considered (i.e., for systems with zonal volumetric
masses much greater than flow passage volumetric masses and for which flow through
these passages is practically instantaneous) the simple flow element should suffice. For
those problems where some interest is focused on the detail of flow passage dynamics,
or where it is believed that dynamic phenomena in the flow passages can not be
ignored, an alternative flow element is required. In this section a convection-diffusion
flow element will be developed that will answer this need.
4.1 Convection-Diffusion Equation
Consider the flow through a flow passage (e.g., a section of duct work) that connects
zone i to zone j;
a e
Wi
node i
X
J / 'V
fel In
AX ^(X,
t)+A°C(
»
x,t);
node j
Figure 4-1 One Dimensional Flow Passage
Isolating a segment AX of the flow passage and demanding the conservation of
species mass flowing through this segment we may write;
we(aC-(aC+AaC))
31 (4.1)
or, in the limit as Ax -»0;
4-1
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
edU d W a . d U
-w + g = pA—
dX dX 3t /42)
where;
we is the total fluid (air) mass flow rate through the flow passage
aC is the a species mass concentration
aw' is the a species mass flow rate relative to total fluid mass flow rate
p is the density of the fluid (air)
A is the cross-sectional area of the flow passage
ag is the a species mass generation rate per unit length of passage
x, t are distance and time, respectively
The first term above accounts for species mass flow due to convection and the
second term accounts for species mass flow due to species diffusion that is
superimposed on the bulk flow.
The generation term, «g, may be replaced by an appropriate generation (kinetics)
rate expression. For example, if the generation involves the single species a one could
replace the generation term with a linear generation rate expression of the form;
a a , aa a_
g = g0 + pA KG
(4.3)
a owx
where g0 is a constant rate component and K is a generation rate constant. This
form of generation rate expression would be appropriate, for example, for formaldehyde
emission from the flow passage walls [Mathews et. al. 1984, Grot et. al. 1985].
The diffusion of species mass relative to the (bulk) total mass fluid flow, aw', will, in
general, depend upon the details of the fluid velocity profile, and its turbulence
characteristics (e.g., the eddy diffusivity), and molecular diffusion along and
perpendicular to the flow passage length. In many practical situations, however, this
diffusion component may be modeled using an expression based upon Pick's law of
diffusion which may be written as;
O . A «r-vd C
w' = - pA D
d* : By Analogy to Pick's Law (4.4)
where aD is the axial dispersion coefficient of a in the flow fluid (air).
Substituting equation (4.4), equation (4.2) may be rewritten as;
4-2
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
,o_9C a . 3 C e 3 C
pAT) - + g = pA - + w -
(4.5a)
an equation that is commonly called the one dimensional convection-diffusion equation
that also appears in thermal convection-diffusion problems. The convection-diffusion
equation is often expresses in dimensionless form as;
1 3 °b a 9aC 3aC
_ + y = —- + —
Pe
3* "* (4.5b)
where;
\ uL
Pe.
pA D D tne dimensionless Peclet Number (4.5c)
L is the (characteristic) length of the flow passage
X is the dimensionless length = x/L
i is the dimensionless time = t/t
°y is the dimensionless generation rate = °gLAve
t is the nominal transit time = L/ u
U is the bulk fluid velocity = we/pA
The Peclet number alone, then, characterizes the convection-diffusion process in a flow
passage not involving a kinetic rate' expression. It provides a measure of the importance
of convection mass transport relative to diffusion mass transport.
The convection-diffusion equation presented above, equation (4.5), is referred to as
the axial dispersion model or axial-dispersed plug-flow model in the chemical
engineering literature where it has played an important role in the simulation of flow
systems found in the chemical process industries since 1908 when Langmuir introduced
it. As one might suspect, the utility of this equation depends critically upon the
determination of the dispersion coefficient to be used for a given set of flow
circumstances. A complete discussion this problem is well beyond the scope of this
report and the reader is, therefore, directed to the excellent general discussion of this
approach by Nauman and Buffham [1983] and the more practically useful, reference
work by Wen and Fan [1975]. Suffice it to say that for turbulent, isothermal flow
conditions in relatively long flow passages, the dispersal coefficient is reasonably well
correlated to the characteristic Reynolds number, Re, of the flow and is practically
independent of species molecular diffusivity as indicated by the Taylor expression
(reported by Wen and Fan [1975 p. 47]):
4-3
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
3 D- 2uR
3.0x10
1.35
Re
Re
0.125
Re> 2000
(4.6)
where;
R is the flow passage radius
Re s 2puR4i
|i is the flow fluid's viscosity
Under turbulent flow conditions, the fluid velocity profile is relatively flat and diffusion is
dominated by mass transport by flow eddies rather than by molecular diffusion, thus the
dispersion coefficient becomes primarily dependent upon the turbulence
characteristics, as measured by the Reynolds number, and is, therefore, often identified
as the eddy diffusivity.
Under laminar flow conditions, on the other hand, the velocity profile tends to be
parabolic, turbulence subsides and as a result both radial and axial molecular diffusion
come to play an important role and the diffusion becomes two dimensional in nature.
Nevertheless, if the fluid can be assumed to be homogenous, so that the radial and axial
molecular diffusivities may be assumed to be identical, then a solution of the complete
two dimensional convection diffusion problem reveals that the asymptotic behavior is
equivalent to that described by the one dimensional convection diffusion equation using
the Taylor-Aris dispersal coefficient [Nauman & Buffham 1983 pp.112-113]:
«D
48
; Re<2000 ;^<0.16 —
n a
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
In practice, building air ducts are normally designed for bulk air flow velocities greater
than or equal to 2 m/s [ASHRAE 1985]. For this minimum operational flow rate, then, the
Reynolds number for the flow in this duct would be (pajr = 1.2 kg/m3 & p.air « 1.8 x 10'5
kg/m-s at 20°C):
Re =
= 2(1.2kg/s)(2.0m/s)(Q.5m) =
-5.
(1.8x10 kg/m-s)
and, by the Taylor correlation, equation (4.6), the dispersal coefficient would be:
_ 2.1 _ 0.125]
Re Re
2(2.0 m/s )(0.5 m)
3.0x10
1.35
0.125
(1.33x1 05) f
= 6.2x10"^^
Typical diffusivities of gas pairs are on the order of 1 to 2 x 10"5 m2/s. Thus, even under
these relatively low flow conditions the dispersal coefficient (eddy diffusivity, in this
case) is seen to be over four orders of magnitude greater than molecular diffusion.
Continuing, the Peclet number, for this case would be:
Pe = Hk = (2-0m/sX10m) = 32.3
D -12
6.2x10 m/s
An examination of the turbulent correlation expression reveals that the dispersal
coefficient for turbulent conditions is dependent upon the average flow velocity and the
flow passage radius. This dependency is plotted below for a range of flow velocities
and radii:
4-5
-------
Indoor Air Quality Model
Phase III Report
4. One Dimensional
PART I - THEORY
Convection-Diffusion Flow
(m2/s)
0.00
0.00 0.50 1.00 1.50 2.00 2.50
R (m)
Figure 4-3 Dispersal Coefficient for Turbulent Flow in Ducts (20°C. 1 atrn^
4.2 Convection-Diffusion Element Equations
Finite element solutions of connective-diffusion equations of the form of equation
(4.5) have received considerable attention in recent years. Following the one-
dimensional example discussed by Huebner and Thornton [1982] element equations for
a two-node flow element may be developed from equation (4.5) using linear shape
functions (i.e., assuming species concentrations vary in a piece-wise linear manner
along the flow passage) and applying either the (conventional) Galerkin method or the
(upwind) Petrov-Galerkin method in the formulation of these element equations. The
resulting element equations are:
(4.8a)
where;
. a e a e,T
= { W, , Wj }
W
1 1
-1-1
W [1 -1
2-11
(4.8b)
the convection component of the element flow transport matrix
4-6
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
the so-called upwind parameter, 0 < <|> < 1
L (4.8C)
the diffusion component of the element flow transport matrix
Le the length of the element (i.e., a portion of the length of the flow path)
211 ApAJ/
12 4
the element volume mass matrix
(4.8d)
(4.8e)
the internal generation rate vector
for total fluid mass flow rate. we, through the flow passage from node i to node j as
indicated in Figure 4-1.
The Upwind Parameter If is set equal to 0, the element equations become identical
to those that would be obtained using the (conventional) Galerkin approach of element
formulation. Unfortunately, these conventional element equations lead to solution
approximations that exhibit spurious spacial variations when convective transport is
large relative to transport by diffusion. The upwind parameter, , has been introduced,
using the Petrov-Galerkin approach, to control this form of numerical instability, but at
the cost of introducing artificial diffusion (vis a vis the second term on the right of
equation (4.8b)) that introduces inaccuracies.
Generation Kinetics If the generation term of equation (4.5) is replaced by the
generation (kinetics) rate expression of equation (4.3) the element equations will have
the slightly modified form of:
(4.9a)
where;
. eaa
pAL K
(DpAL9aaK
-1-1
1 1
(4.9b)
4-7
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
the kinetics component of the element flow transport matrix
(4.9c)
the constant component species generation rate vector
Lumped Element Mass Matrix The convection diffusion element equations defined
by either equations (4.8) or (4.9) may be assembled, in the usual manner, along with the
simple flow element equations (equations (2.4) or (2.5)) and kinetics element equations
(equations (3.36)) to form the system equations. The convection-diffusion element
introduces, however, nondiagonal contributions to the system mass matrix, [M], that
adds some complexity to the assembly and solution algorithms used in the
computational implementation of the contaminant dispersal theory. To avoid this
complexity one may replace the so-called consistent element volume mass matrix,
equation (4.6d), with a diagonal lumped mass approximation to it, given by:
[m6]
pALe
the element lumped mass matrix (4.10)
This approximation may be expected, however, to introduce some additional error. The
program CONTAM87 provides convection-diffusion elements having this lumped mass
approximation.
4.3 Use of The Convection-Diffusion Flow Element
In this application, we are using the Finite Element method to approximate the spatial
variation of contaminant concentration along a flow passage in a piece-wise linear
manner where each linear segment of the approximation will correspond to an individual
convection-diffusion flow element. Clearly, if the form of the actual concentration
variation along the flow path is linear we could model the flow passage with a single
element. In general, however, we will not have a priori knowledge about the form of the
concentration variation and, therefore, should employ a series of flow elements to obtain
a piece-wise approximation to the actual concentration variation.
Without some experience, the analyst may not know how many convection-diffusion
elements to use in a given situation. In such a situation, however, a first analysis may
be completed with a trial subdivision of the flow path then the analysis may be repeated
with a finer subdivision. The finer subdivision may be expected to provide a better
approximation to the solution (providing numerically stability has been achieved) and,
therefore, may be used to access the accuracy of the solution. This process of
subdivision could, then, be repeated until the solution converges to within an
acceptable accuracy.
4-8
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
Steady State Analysis When considering steady-state flow without internal
generation, Huebner and Thornton show that instability may be avoided if an upwind
parameter is selected satisfying the conditions;
>1 -4 ; Pee>2
Pe
= 0 ; Pe<2
6 (4.11)
where;
e. e —,e
w L uL
PA°b °b (4.-12)
is the element Peclet number
(Note: Pe = (Pe/n) for a flow passage idealized by an assembly of n equal-length
convection-diffusion elements.)
To fix these ideas, consider the problem presented by Huebner and Thornton [1982];
the dispersal of a contaminant along a straight flow passage, under steady flow
conditions, without generation, and with inlet contaminant concentration maintained at
C0 and outlet concentration maintained at zero, as diagrammed in Figure 4-4.
C(0) 1 C(L) _ Q
Co "
f
t
t
we-
->
i
i
L i L
T
1
«
i
Steadv
Node Number
2^ 3
1\t 2
• 3
State
4
* 4
i
Convection-Diffusion Problem
5
+ 5
6
f 6
7 8 9 10 11
^7 »8 f9 •lOf
Element Number
Ten Convection-Diffusion Element Idealization
Figure 4-4 A Steady State Convection-Diffusion Problem and Corresponding Finite
Element Idealization
For this problem the convection-diffusion equation simplifies to:
4-9
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
1 d2C
Pe
_ 0
2 " dx
dx
which may be solved for the boundary conditions:
C(x=0) =C0
C(x=L) = 0
to obtain an exact solution:
(4.133)
(4.13b)
C(x/L) _ e L - e L .
Pe
1 - 6L
; 0 = 0.0
and $ = 1.0. The exact with the approximate solutions are compared below in Figure 4-
5.
C(x)/Co 0.8
1.0 B
00 .
J
*
*
*
4T
— Pe=0.2, Exact • Pe=0.2, 0=0.0 ° Pe=0.2, 0=1.0
•"-> Pe=20, Exact A pe=20. 0=0.0 A pe=20, 0=1.0
t
%
^
^
^
*
-*^!
J
^
*
m
*
m
«
%
«
%
«
«
'— —•D
A
m
*
*
m
m
f
m
a
m
g
. f i
9
g
k
•
%
*
%
*
«
«
•
*
%
•
i
%
\
y^^:\
3
g
t
»
m
»
*
•
*
•
•
*
•
*
Bl-»_
^*"^i
5^^^
^^i
1 >
3^^^
\
\
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
X/L
Figure 4-5 Comparison of Exact and Finite Element Solutions for a-Steady Convection-
Diffusion Problem
The results clearly demonstrate the numerical instability that may result when
4-10
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
upwinding is not used for high element Peclet numbers. For convection-dominated flow,
which should be expected to be typical in building HVAC ductwork under operating
conditions, the analyst may, then, choose to either use a fine subdivision of a given
duct or employ upwinding to maintain numerical stability, keeping in mind that the
upwinding will introduce artificial diffusion that may add error. (A close examination of
the results above reveals that full upwinding underestimated the concentration variation
slightly.)
Dynamic Analysis The convection-diffusion flow element may also be employed for
dynamic analysis, but the analyst must take special care to assure an accurate solution
has been obtained. In dynamic analysis, accuracy is affected not only by element size
(i.e., the subdivision of the flow path), and the degree of upwinding chosen, but also by
the integration time step selected to complete the dynamic solution. Furthermore, the
use of the lumped mass approximation, while avoiding the complexity demanded by
nondiagonal mass contributions, tends to introduce spurious anomolies in the computed
solution in some cases [Huebner & Thornton 1982].
Partly because of the challenge of these difficulties and partly because of the
importance of the convection-diffusion equation in the area of fluid mechanics, finite
element solutions of the convection-diffusion equation have become the focus of
considerable research in recent years. Strategies have been put forward to improve the
accuracy of the finite element approximation presented above that are, regrettably,
beyond the scope of this presentation and the reader is, therefore, advised to review the
current and emerging literature. The papers by Hughes and Brooks [1982], Tezduyar
and Ganjoo [1986], and Yu and Heinrich [1986] are particularly useful in this regard.
In spite of the numerical pitfalls that await the use of the convection-diffusion flow
element presented above we shall proceed and employ these elements (with the lumped
mass approximation) to compute the transport of a contaminant pulse through a length of
ductwork. The conditions of this problem are diagrammed in Figure 4-6: fluid flows
through a duct of length L and radius R at a mass flow rate we; a contaminant is injected
into the inlet stream at a rate G(t) for a short time interval introducing a pulse of
contaminant of mass I into the inlet stream; the pulse is convected and dispersed as it
moves along the duct. We seek to determine the concentration time history of the
contaminant as it emerges from the outlet of the duct.
The exact solution to this problem is available for an impulse (i.e., a pulse defined by
the dirac delta function), for "closed" inlet and outlet conditions, but it is expressed as
an infinite sum that is practically difficult to use [Wen & Fan 1975 pp. 133-137]. For
Pe=0 the duct becomes a well-mixed system, the initial concentration throughout the
duct becomes, simply, (l/pAL), and the outlet concentration decays exponentially:
C(L.t) = e-t/f
(l/pAL) ;Pe=0 (4.15)
4-11
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
For relatively large Peclet numbers the outlet concentration is well approximated by the
the following expression reported by Nauman and Buffham [1983 pp. 101-103]:
C(L.t)
(1/pAL)
Pe
e\
-Pe(1 -t/F)
4t/f
_ 3
;Pe > 16
(4.16)
and for very large Peclet numbers the outlet concentration approaches a Gaussian
distribution [Wen & Fan 1975 p. 133]:
C(L.t) _,/PT
(l/pAL) "
-2
-Pe(1 -t/t)
4
;Pe » 16
(4.17)
Approximate solutions to this problem were computed using a 10-element
subdivision, as shown in Figure 4-6, and a twenty-element subdivision. The "closed"
boundary condition was modeled using the simple flow element as this element models
(instantaneous) plug flow conditions as required. The impulse was approximated by a
pulse of finite but small duration. In all studies the upwind parameter, , was chosen to
satisfy the lower bound (i.e., equality) of the stability requirement of equation (4.11). The
results are compared below, Figure 4-7, to the solutions discussed above, equations
(4.15) to (4.17).
G(t)J,
_l
~l
t
,/GWdt-l °11""
,t «
1 W
Rjl >
->*
/\^
1 \
V. Plug-Flow
Impulse Transport Problem
"Closed Boundary"
1 •
"Exterior
Zone"
Node Number
2 3'
-*
^^•2
\.
4
+ 3
567
•4*5*
8
6 f 7
Element Number
9 10 1
+ 8 • 9 4
Simple Flow
1 12
> 10 <>7>
/
Element
Ten Convection-Diffusion Element Idealization
Figure 4-6 The Transport of a Pulse in a Duct and the Corresponding Finite Element
Idealization
4-12
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
It is seen that in this case the approximate, finite element solution for the low Peclet
number, Pe=1, approaches the exact well-mixed solution, as expected. The
approximate solution for the higher Peclet numbers has some difficulty in capturing the
amplitude of the exit pulse, although, the timing and the form of the pulse appear to be
well-approximated. Some part of this error may be attributed to approximating the
impulse of the analytic solutions by a pulse of finite duration in the computed solutions.
In these studies the pulse duration was set at 0.001 units of dimensionless time,
increasing the pulse duration by a factor of 4 resulted in an additional underestimation of
the pulse amplitude at Pe=20 of approximately 5%.
Some part of the error may be attributed to the coarseness of the finite element
subdivision. A comparison of the results of the 10-element and 20-element
approximations for Pe=10 indicate that a convergent solution was obtained (i.e., further
subdivision would not alter the solution), yet when these results are compared to the
exact results reported by Wen and Fan [Wen & Fan 1975 Fig. 5-8 p. 136] the amplitude
appears to be underestimated by approximately 10%. This same comparison for Pe=20
indicates that a convergent solution was almost but not quite achieved. An additional
subdivision would presumably reveal convergence, and the error in amplitude
estimation was approximately 20%. It is interesting to note that the element Peclet
numbers for these two (nearly) convergent solutions - the 10-element solution at Pe=10
and the 20-element solution at Pe=20 - are both equal to 1.0, a condition that demands
no upwinding (i.e., for which $ may be set to 0) to maintain numerical stability. Results
were also computed for cases violating the stability requirement of equation (4.11) and,
as expected, spurious variations in concentration responses - "wiggles" - were
observed.
»" Exact: Pe=0
— Nauman: Pe=20
Gauss: Pe=20
••• 20 Elem: Pe=20
a loElem:Pe=20
-A- 20 Elem: Pe=10
A 10Elem:Pe=10
••• 20 Elem: Pe=1
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Dimensionless Time t/(LAl)
Figure 4-7 Comparison of Analytic Solutions with Finite Element Solutions for the Pulse
Transport Problem
It may be useful to relate these nondimensional studies to more conventional units.
4-13
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
The study for Pe=20 corresponds to studying the transport of a pulse through a circular
duct of 1 m radius having a length of 10 m with a bulk flow velocity of 2 m/s (the practical
minimum operational flow rate in HVAC ducts). For these conditions, by Figure 4-3, the
dispersal coefficient may be expected to be about 1.0 m2/s. The results reported in
Figure 4-7 were computed using a pulse duration of 0.005 sec (i.e., 0.001 times the
nominal transit time, F = L/u = 10 m 12 m/s = 5 s). The dynamic solution was computed
using a time step of 0.001 second, in part to capture the short-time pulse accurately and
partly to achieve a practically convergent solution.
In practical situations the inaccuracies revealed in these studies are likely to be
considered very small and, thus, the convection-diffusion flow element should provide a
practically useful analytical tool. Nevertheless, to minimize error the analyst is well
advised to seek a convergent solution through both mesh refinement (i.e., repeated
subdivision of the flow path), starting, perhaps, with a subdivision that results in an
element Peclet number of 1.0, and time step refinement, starting with a time step
sufficiently small to capture the dynamic variation of any excitation with reasonable
accuracy, being careful to select an upwind factor so that the stability requirement of
equation (4.11) is always satisfied. When employing convection-diffusion elements in
an idealization of a building airflow system it is very likely that extremely small time steps
will be required to obtain a convergent solution.
4.4 Analytical Properties of the Convection-Diffusion Element
Equations
The numerical properties of the convection-diffusion flow element have been seen to
be dependent on the element Peclet number. To investigate this dependency in greater
detail we may rewrite combined convection and diffusion components of the element
flow transport matrices, equations (4.8b) and (4.8c), in terms of the element Peclet
number, as:
-l
°DA 1 -1
w
~2
1 1
-1-1
(4.18)
The stability requirement of equation (4.11), which may be rewritten as,:
; Pe>2
-------
Indoor Air Quality Model
Phase III Report
PART I - THEORY
4. One Dimensional Convection-Diffusion Flow
assures that the flow transport matrix will be an M-matrix, which is to say it is a real
square matrix with positive diagonal elements and nonpositive off-diagonal elements
such that 11 u + £['] ] is strictly diagonally dominant for all scalars £ > 0.
It was shown earlier [Axley 1987] that element flow transport matrices satisfying this
condition (coupled with mass matrices that are positive diagonal matrices) lead to
system transport matrices that are not only nonsingular, but may be decomposed to
[L][U] form by a variant of Gauss elimination without the need for pivoting in an efficient
and numerically stable manner and will have stable homogeneous forms.
4.5 Comparison to Tanks-in-Series Idealizations
In the chemical engineering literature the so-called tanks-in-series idealization is
frequently employed to model the behavior of one dimensional convection-diffusion
transport processes or other processes whose inlet-outlet transformation characteristics
appear to match those described by one dimensional convection-diffusion processes.
Below, in Figure 4-8, we compare a 5-element/6-node finite element idealization to a
corresponding 6-node tanks-in-series idealization where the the fluid mass of volume of
the flow path has been subdivided into four "unit" tanks containing one fifth of the total
fluid mass each and two "half-unit" tanks containing one tenth of the total fluid mass
each.
Node
Conv-Diff Flow Element
w
Simple Flow Element
" l\v\A
Convection-Diffusion Element Idealization
w f
Tanks-ln-Series Idealization
wef
Figure 4-8 Comparison of Tanks-in-Series Idealization with Finite Element Idealization
In the tanks-in-series idealization a portion of the flow, f, assumed to recirculate between
adjacent tanks is used to model the nature of turbulent and molecular diffusion.
The subassemblage of this tanks-in-series idealization consisting of half of two
adjacent "unit" tanks and the connecting simple flow elements, which we shall refer to
as a tanks-in-series element, may be compared directly to the convection-diffusion flow
element, as indicated in Figure 4-9, below:
4-15
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
PAL
Conv-Diff Flow ElemQnt
"Tanks-ln-Series" Element
Figure 4-9 The Equivalence of the Convection-Diffusion Flow Element and a Tanks-in-
Series Element
Element equations for the tanks-in-series flow element may be assembled directly from
the simple flow element equations, equation (2.4a), producing:
= w
1 0
-1 0
+ fw
d°cr
pAL'hOl7 dt
(4.20)
Comparing these equations with the convection-diffusion element equations, equations
(4.8), we see that they become equivalent when:
e. e
W L
Pee
(4.21)
and full upwinding, $ = 1.0, is used.
It is interesting to note that the extreme of pure plug flow in the convection-diffusion
case corresponds to conditions having a dispersal coefficient equal to zero. A fine
subdivision of the flow path into a large number of finite elements would be required to
provide a good approximation of the plug flow behavior. In comparison, by equation
(4.21), plug flow would correspond to a tanks-in-series element with f = 0.0 - that is to
say an element without the recirculating backflow. Nauman and Buffham [1983 pp. 58-
59] show that as the number of tanks-in-series without backflow becomes large, the
behavior of the assemblage of tanks approaches plug flow. The other extreme of well-
mixed conditions may be modeled with an infinite dispersal coefficient in the
convection-diffusion case. This corresponds to infinte recirculation in the tank-in-series
element and we obtain the behavior of a simple well-mixed zone with either a single
convection-diffusion element or a single tanks-in-series element.
The Imperfectly Mixed "Zone Element" The comparison of the convection-
diffusion element and the tanks-in-series idealization, supports the conclusion drawn
above that, in general, modeling high Peclet number flows will demand a fine
subdivision of elements and modeling low Peclet number flows will not. It also points out
4-16
-------
Indoor Air Quality Model PART I - THEORY
Phase III Report 4. One Dimensional Convection-Diffusion Flow
the fact that the convection-diffusion element may be used to model a zone that is not
perfectly mixed. In fact, although we developed the convection-diffusion element to
model flow transport situations, it should now become apparent that this element
provides one means to model imperfectly mixed zones. If the exit flow response of a flow
zone to a supply flow pulse takes the form of the solutions presented above, equations
(4.15) to (4.17), then one may employ an assemblage of convection-diffusion flow
elements to model the global characteristics of that zone, even though the internal
mechanisms governing the imperfect mixing are not apparent.
It may be shown that the variance of the nondimensional response, a2, is related
directly to the Peclet number of the flow, for a "closed" system, as:
Pe (4.22)
For large Peclet flows the form of the (nondimensional) exit response is well
approximated by the form of a Gaussian distribution (e.g., see the results of Figure 4-7)
which has a variance of:
2 2
aG = •=—
Pe ; for the Gaussian approximation equation (4.17) (4.23)
Either of these two expression provides a means to determine an effective Peclet
number for a zone, from a rather straightforward statistical reduction of actual pulse
response measurements, that may, then, be used for modeling purposes.
Chapter 2 presented the general formulation of a multi-zone contaminant dispersal
analysis theory, based on element assembly techniques, and briefly presented the
element equations developed earlier: the simple flow element with and without filtration.
Chapter 3 outlined a means to organize the multi-zone contaminant dispersal analysis
equations for the consideration of multiple contaminants and introduced a kinetics
element to account for chemical and physical interactions between contaminants and
the materials of the building construction and furnishings. The present chapter
introduced a fourth contaminant dispersal element, the one-dimensional convection
diffusion element that may be employed to either study the details of dispersal in one-
dimensional flow regimes (e.g., as found in HVAC ducts) or to simulate the general
behavior of certain imperfectly mixed zones.
The program CONTAM87 implements this theory. Chapters 5 through 8 provide a
users manual to this program and Chapter 9 provides some examples of its use.
4-17
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
5. General Instructions
PART II - CONTAM87 USERS MANUAL
5. General Instructions
The program CONTAM87 is a command processor; it responds to commands in the
order that they are presented and processes data associated with each command.
Commands may be presented to the program interactively, using keyboard and monitor,
or through the use of command/data input files; that is to say, it offers two modes of
operation - interactive and batch modes.
For most practical problems of contaminant dispersal analysis the batch mode of
operation will be preferred. For these problems, analysis involves three basic steps;
Step 1:
Idealization of the Building System and Excitation
Actual Building
•—Zone Node Flow Element-
2 - Node Number Element Flow-
Air-Flow System Idealization
Idealization of the building flow system involves a) discretization of the system as an
assemblage of appropriate flow elements connected at system nodes, b) identification
of boundary conditions, and c) numbering of system nodes optimally (i.e., to minimize
the bandwidth or, equivalently, node number difference of the system equations).
The excitation (i.e., specified contaminant concentrations and generation rates) may
be modeled to be steady or defined in terms of arbitrary time histories. For the latter case
initial conditions of nodal contaminant concentration will have to also be specified.
5-1
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
5. General Instructions
Step 2:
Preparation of Command/Data Input File
Edit
In the batch mode, the program reads ASCII text files of commands and associated
data, collected together in distinct data groups, that define the building flow idealization
and excitation. The command/data input file may be prepared with any available ASCII
text editing program and given a file name, , specified by the user. The
must, however, consist of 8 or less alphanumeric characters and can not
include an extension (i.e., characters separated from the filename by a period,".").
Step 3:
Execution of CONTAM87
Plot Output File
.PLT
Disk Storage Data File?
Command/data
Input File
0010ft
01100
00011
0010^
01100
00011
filename >.FEL filename >.WDT
COIMTAM
4-
0010^1
01100
00011
0010'
01100
00011
.RXN filename >.EDT
Results Output File f^T^i
.OUT
CONTAM87 is then executed. Initially CONTAM87 will be in the interactive mode.
To enter the batch mode the command "SUBMIT F=" may be used to submit
the command/data input file to the program. The program will then proceed to form
element and system arrays and compute the solution to the posed problem.
CONTAM87 reads the ASCII command/data input file and creates an ASCII (i.e.,
printable) output file .OUT. The results of an analysis, .OUT, may
be conveniently reviewed using an ASCII editor and, from the editor, portions or all of
the results may be printed out. Key response results are also written to the ASCII file
.PLT in a format that may easily be transferred to spreadsheet and plotting
programs (data values within each line are separated by the tab character and lines of
data are separated by a carriage return) for plotting or subsequent processing.
5-2
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 5. General Instructions
Depending upon the commands processed, CONTAM87 will also create a variety of
binary files for disk storage needed for subsequent processing. A summary of files read
and created includes;
Files Read
an ASCII input file specified by the user that contains commands and
associated data
Files Created
.OUT a printable ASCII output file that contains analysis results
.PLT an ASCII output file that contains key analysis results in a form
that may be transferred to spreadsheet and/or plotting programs
.FEL a binary file used for disk storage of flow element data
.KIN a binary file used for disk storage of kinetics element data
.WDT a binary file used for disk storage of element flow time
history data
.EDT a binary file used for disk storage of excitation time history
data
In the interactive mode is set to the default value of "CONTAM87" and
commands are read from the keyboard. A help command, "HELP" or "H", will produce a
screen listing of intrinsic commands.
5-3
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 6. Command Conventions
6. Command Conventions
Commands and their associated data (if any) may be single-line or multiple-line
command/data groups.
Single-Line Commands Single line command/data groups begin with the command
keyword and may have any number of associated data items identified by data identifies
of the typical form;
COMMAND A=n1,n2,n3 B=n4 C=n5,n6 D=c1c2c3
where n1,n2,n3,... is numeric data and c1c2c3 is character data. In this example the
keyword COMMAND is the command keyword and the data identifiers are A=, B=,
C=, and D=.
Multiple-Line Commands Multiple-line command/data groups are delimited by the
command keyword and the keyword END and may have any number of data subgroups
terminated by the less-than character "<" within. They have the typical form of;
COMMAND A=n1,n2
n1 I=n2,n3,n4 B=n5 C=c1c2c3c4
n1 I=n2,n3,n4 B=n5 C=c1c2c3c4
n1 I=n2,n3,n4 B=n5 C=c1c2c3c4
<
n1,n2,n3 D=n4,n5,n6 E=n7 F=c1c2c3
n1,n2,n3 D=n4,n5,n6 E=n7 F=c1c2c3
n1,n2,n3 D=n4,n5,n6 E=n7 F=c1c2c3
<
Clc2c3c4c5c6
END
Classes of Commands Two general groups of commands are available, the
Intrinsic Commands and the CONTAM87 Commands. The intrinsic commands are
used to control the operation of the command processor CONTAM87 and to examine
arrays generated by the CONTAM87 commands. The CONTAM87 commands provide
contaminant dispersal analysis operations.
Command/data Lines Normally the line length (i.e., the number of character and
spaces on a line) is limited to 80. A backslash "\" at the end of information on any line
will, however, allow the next line to be interpreted as a continuation of the first line
providing an effective line length of 160.
A less-than character "<" indicates the end of information on any line. Information
entered to the right of the less-than character is ignored by the program and may,
therefore, be used to annotate a command/data input file.
6-1
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 6. Command Conventions
An asterisk "*" at the beginning of any line will cause the line to be echoed as a
comment on the console and to the results output file. Lines marked in this way may,
then, be used to annotate the results output file. Comment lines may also help indicate
the progress of computation when using the batch mode of operation.
Data Identifiers Data identifiers and their associated data may be placed in any order
within each line of the command/data group (with the exception that the first line of a
command/data group must begin with the command keyword). In some instances data
may not be associated with a data identifier, such data must be placed first in a line.
Data Decimal points are not required for real numeric data. Scientific notation of the
form nnE+nn or nn.nnE+nn (e.g., 5.79E-13) may be used. Simple arithmetic expressions
employing the conventional operators +, -, *, and / may be used. The order of
evaluation is sequential from left to right - unlike FORTRAN or other programing
languages where other "precedence" rules are used.
If fewer data values are supplied than required, the missing data will assumed to be
zero, blank, or set to default values as appropriate.
6-2
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
7. Introductory Example
7. Introductory Example
For purposes of contaminant dispersal analysis the specific command/data groups
that need to be included in a command/data input file will depend upon the details of the
flow system idealization, the nature of the excitation, and the type of analysis to be
computed. A specific introductory example, should however, provide some useful
insight into the more general aspects of contaminant dispersal analysis using
CONTAM87
Consider the two-story residence with basement shown, in section, below. In this
residence interior air is circulated by a forced-air furnace and exterior air infiltrates the
house through leaks around the two first-floor windows. The flow system may be
idealized using flow elements to model the ductwork, room-to-room, and infiltration flow
paths as shown below.
N02
Generation
0.11
g/hr
0.1 m3
20 m3/hr
70 m3/hr
• — Zone Node
2 — Node Number
70 m3/hr
Flow Element —
Element Flow _
Actual Building
Air-Flow System Idealization
Figure 7-1 Hypothetical Residential Example
For this building idealization we shall consider the hypothetical problem of
determining the steady state distribution of N02 generated by a kerosene heater placed
in room "2", distributed by the furnace flow system operated at constant conditions, and
diluted by infiltration at a constant rate. The N02 generation rate is assumed to be 0.11
g/hr, exterior N02 concentration is assumed to be negligible, and the assumed air
volumetric flow rates are indicated on the drawings above. Inasmuch as N02 is a
reactive gas it will also be assumed that N02 is constantly transformed into other
products that, here, are of no particular interest, as;
N02 -> products
7-1
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 7. Introductory Example
This reaction is often assumed to be governed by rate expressions of the form:
thus the matrix of rate coefficients for this case is a 1 x 1 matrix:
[«] = I
where 2K, the reaction rate constant for this reaction, will be assumed to have a value
of 0.40 hr- 1. (These values of N02 generation and reaction rate are based on values
reported by Traynor [Traynor et. al. 1983] and Nitschke [Nitschke et. al. 1985]. The
generation rate is representative of that produced by portable kerosene heaters. The
reaction rate constant is thought to be representative of that to be expected indoors, but
the kinetics of N02 chemical or physical-chemical behavior indoors is not yet well
understood.)
The CONTAM87 command/data file to complete 'this steady state analysis is listed
below. Command/data groups needed to complete a time constant analysis and
dynamic analysis for this building idealization are presented as examples in the
reference section of this manual.
7-2
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
7. Introductory Example
Command/data File for Residential Example
Note: CONTAM87 keywords and identifiers are displayed in boldface below.
Description
Column
Comments:
Comments
Comments
Comments
Comments
Comments
System Definition:
No.Nodes &.Species, Species IDs
Boundary Conditions
Nodal Volumetric Mass
Flow Element Data:
Element Number & Connectivity
Kinetics Element Data:
Rate Coef. (Matrix): Rxn 1
Kinetics Elem. Location & Type
Steady State Solution:
Flow Element Mass Flow Rates
Contaminant Excitation
Return to Interactive Mode
Command/data File
1
* Six-Zone (7-Node) Example
* Units: g, m, hr
* Concentration [=] g-N02/g-air
* Generation rate [=] g-N02/hr
*
FLOWSYS
N=7 S=l ID=N02
7 BC=C < Ext."Zone" Cone.Spec.
< (Air Dens. 1.2E+3 g/m3)
Node 1 Vol. Mass
Nodes 2 & 3 Vol. Mass
Nodes 4 & 5 Vol. Mass
Node 6 Vol. Mass
Node 7 Ext. Vol. Mass
Flow Element 1
Flow Element 2
Flow Element 3
Flow Element 4
Flow Element 5
Flow Element 6
Flow Element 7
Flow Element 8
Flow Element 9
Flow Element 10
Flow Element 11
1 V=1.2E+3*1.0
2,3 V=1.2E+3*40.0
4,5 V=1.2E+3*30.0
6 V=1.2E+3*0.1
7 V=1.2E+3*l.E+6
END
FLOWELEM
1
2
3
4
5
6
7
8
9
10
11
=1,2
=1,3
=7,2
=2,7
=7,3
=3,7
=2,4
=3,5
=4,6
=5,6
-6,1
END
KINELEM
K=l
0.4
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
<
1 1=1 K=l
6 1=6 GEN=1 K=l
END
STEADY
1,2 W=70*1.2E+3
3,6 W=20*1.2E+3
7,10 W=70*1.2E+3
11 W=140*1.2E+3
2 CG=0.11
7 CG=0.0
END
RETURN
Rxn 1:N02 -> products
< Node 1: Rxn 1
< Nodes 2 to 6: Rxn 1
<(Air Dens.l.2E+3 g/m3)
< Supply Ducts
< Infiltration
< Return Loop
< Main Return Duct
< Node 2: N02 Gen. Rate
< Node 7:Ext. N02 Cone.
Details are given on the following pages for each CONTAM87 command.
7-3
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
8. Command Reference
CONTAM87 provides two general classes of commands; Intrinsic Commands and
CONTAM Commands.
8.1 Intrinsic Commands
The intrinsic commands are used to control the operation of the command processor
CONTAM87 and to examine arrays generated by the CONTAM87 commands.
These intrinsic commands have been developed to provide general command
processor operations that, together with the general command conventions outlined
earlier, define a standard user-machine interface that may be used in the development
of other simulation software.
8.1.1 HELP
The command HELP, or simply H, will produce a list of available intrinsic
commands, in abbreviated form.
8.1.2 ECHO
The command ECHO-ON acts to cause computed results normally directed to the
results output file to be echoed to the screen. The command ECHO-OFF turns this
feature off. At start-up CONTAM87 is set to ECHO-ON. Selective use of ECHO-ON
and ECHO-OFF can act to speed computation as writing results to the screen
consumes a significant amount of time.
8.1.3 LIST
The command LIST, or simply L, will produce a list of all arrays currently in the
array database.
8.1.4 PRINT A=
The command PRINT A= or simply P A= will print
array named , a one-to four character name, to the screen. Arrays
currently in memory are listed by name using the LIST command.
8-1
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
8.1.5 DIAGRAM A=
The command DIAGRAM A= or simply D A= will print
a diagram of array named , a one-to four character name, to the screen
indicating position of zero and nonzero terms. (Character arrays can not be diagramed.)
8.1.6 SUBMIT F=
The command SUBMIT F= or simply S F= will cause the
program to switch to batch mode and read all subsequent commands from the batch file
.
8.1.7 PAUSE
The command PAUSE will cause the execution of CONTAM87 to pause until a
carriage return is entered from the keyboard. Selective use of PAUSE in a batch
command/data input file will allow the user time to view results of interim calculations.
(Note: PAUSE is a single line command and, therefore, cannot be placed within other
multiline command/data groups.)
8.1.8 RETURN
The command RETURN returns the operation of the program from batch mode to
interactive mode. RETURN or QUIT will normally be the last command of batch
command/data input files.
8.1.9 QUIT
The command QUIT or simply Q terminates execution of the program and returns the
user to the control of the operating system.
8-2
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
8.2 CONTAM87 Commands
The GONTAM87 Commands implement basic contaminant dispersal analysis
operations. These operations are based upon the dimensionally homogeneous theory
presented in the first part of this report, thus, the analyst may use any dimensional units
that are convenient so long as a consistent set of units are employed. Following the
underlying theory one may elect to express all quantities in terms of units of mass and
time;
species concentration [=] mass-species/mass-air
species generation rate [=] mass-species/time
zone ("volumetric") mass [=] mass-air
air flow rates [=] mass-air/time
kinetics rate constants [=] 1/time
or, if consideration is limited to isothermal cases, one may elect to use volumetric
quantities;
species concentration [=] volume-species/volume-air
species generation rate [=] volume-species/time
zone volume [=] volume-air
air flow rates [=] volume-air/time
kinetics rate constants [=] 1/time
Using the simple in-line arithmetic expressions allowed for numeric data, one may
easily convert from one quantity to another while maintaining a record of the conversion
in the input command/data file for future reference. For example, the zone "volumetric"
mass (i.e., the mass of the air within the volume of each zone) required by the
FLOWSYS command could be expressed in terms of zone dimensions and air density
as v=s .0*10.0*2.5*1.2 for a room that is 5 m x 10 m x 2.5 m containing air with a
density of 1.2 kg/m3.
The following conventions will be used for the command definitions presented in this
section;
ellipses,'...', indicate unlimited repetition of similar data items or data lines
within a data subgroup
- square brackets, [...], indicate optional data,
numeric data is indicated by lower case n, as n1,n2,..., and
- character data by lower case c, as c1.
8.2.1 FLOWSYS
The number of the flow system nodes and species, boundary conditions, and
volumetric masses of system nodes are defined with the following command/data group;
8-3
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
FLOWSYS
N=n1 S=n2 [ID=c1,c2, ...]
n3,n4,n5 BC=c3,c4, ...
n3,n4,n5 V=n6
END
where; n 1 = the number of flow nodes
n2 = the number of contaminant species
c1,c2, ... = species ID's; a four character (or less) identification for each
species used for labeling results; species are identified by
species-number and, optionally, by species ID given in species-
number order; omitted species IDs will be set to species-number
n3,n4,n5 = first node, last node, node increment of a series of nodes with
identical boundary conditions
c3,c4, ... = boundary condition codes for each species by species number
order; a single character code of C for concentration prescribed
nodes or G , for generation prescribed nodes; (default = G),
n6 = nodal volumetric mass; (default = 0.0)
The direct species mass generation rate 01 the species concentration, but not both.
may be specified for each species at each node to establish the discrete boundary
conditions of the analysis problem being posed.
Omitted boundary condition data will be assumed to be generation-prescribed.
Typically, nodes associated with the outdoor environment will be assigned specific
contaminant concentrations and nodes associated with indoor air zones will be
assigned specific species generation rates (a zero generation rate will often be
appropriate for the interior species/node combinations).
Volumetric mass data omitted will be assumed to be zero. The present version of
CONTAM does not eliminate system variables associated with zero mass nodes. For
time constant analysis, and in some instances dynamic analysis, a zero nodal mass
value will result in numerical difficulties. From a practical point of view, all nodes of a
flow system idealization will have some volume of air associated with them, although
some may seem insignificantly small, and, therefore, to avoid numerical difficulties all of
these volumes should be modeled with nonzero volumetric mass values.
At the other extreme, some nodes, such as those corresponding to the outdoor
environment , may have practically infinite volumes associated with them. The analyst
should realize practically accurate analysis results for these "infinite" nodes if their
volumes are modeled with volumetric masses several orders of magnitude larger than
that of the largest "non-infinite" node.
8-4
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
8.2.2 FLOWELEM
Presently two types of flow elements are available for contaminant dispersal
analysis;
- a simple flow element that models fluid flow from one node to another ignoring
the dynamic effects of diffusion and convection that result in species flow delay
along the flow path (i.e., flow of a fluid parcel in simple flow elements is
instantaneous) and,
- a convection-diffusion flow element that models fluid flow from one node to
another accounting for these dynamic effects (presently limited to constant
cross-section flow passage idealizations and lumped mass idealizations).
To use these elements effectively and reliably the analyst should be familiar with
their underlying theoretical basis and numerical characteristics. This is especially
important when using the convection-diffusion element. An inexperienced analyst is
well-advised to avoid the use of convection-diffusion elements altogether.
Both simple flow elements and convection-diffusion flow elements may be added to
the flow system assemblage with a command/data group having unique formats of data
lines for each flow element type of the form;
FLOWELEM
n1 I=n2,n3 [GEN=n4] [T=SIMP] [E=n5,n6,...]
or
n1 I=n2,n3 [GEN=n4] T=CNDF M=n7 L=n8 [D==n9,n10,...] [F=n11]
END
where; n1 = the element number
n2, n3 = the system node numbers to which the element is connected
n4 = generation increment (default = 1)
For Simple Flow Elements: [T=SIMP]
n5,n6,... = the element filter efficiency for each species being considered,
in species-number order, (must be > 0.0; default = 0.0),
For Convection-Diffusion Elements: T=CNDF
n7 = the fluid mass per unit length of the (equivalent) constant cross-
section element (must be > 0.0; default = 0.0),
n8 = the flow passage length (must be > 0.0; default = 0.0),
n9,n10,... = species dispersal coefficient for each species considered, in
8-5
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
species-number order, (must be > 0.0; default = 0.0)
n11 = upwind factor, ; where 0 < < 1; (default = the lower bound of
the stability criteria , equation (4.11))
For assemblages consisting of only Simple flow elements the command/data group
would consist of data lines of the form with the type identifier T=SIMP. For
assemblages consisting of only Convection-Diffusion flow elements the command/data
group would consist of data lines of the form with the type identifier T=CNDF. For
mixed assemblages the appropriate mix of the two forms of data lines would be used;
there are no special restrictions on the use of mixed assemblages.
If the element type identifier is omitted the element will be assumed to be of type
SIMP.
Normally, the analyst should accept the default upwind factor for the Convection -
Diffusion element. This default will ensure that numerical solutions to the posed problem
may be determined in an efficient and stable manner. (The option to specify the upwind
parameter is provided to allow one to study the numerical characteristics of the
upwinding strategy rather than the practical behavior of flow systems.)
Element data must be supplied in numerical order. Omitted data is automatically
generated by incrementing the preceding node numbers by the current generation
increment. Generated elements will have the properties of the current element. If, for
example, an HVAC duct, included as part of a air flow system, was to modeled by a
series of, say, ten convection diffusion elements, as illustrated below in Figure 8-1, then
one could conveniently use the generation option to "generate" the intermediate
elements by specifying only the first and last Convection-Diffusion flow element in the
series. The portion of the input command/data file needed to implement this example is
listed below.
FLOWELEM
21 1=12,15 T=CNDP M=1.2E+03 L=l. 0
30 1=39,42 T=CNDF M=1.2E+03 L=1.0 GEN=3
END
Node Number
12 15^ 18 21
24 27 30 33 36 39 42
f 21V f 22 • 23 + 24 » 25 • 26 f 27 • 28 + 29 f 30
Element Number
Figure 8-1 Hypothetical Conduction-Diffusion Element Subassemblv
The command FLOWELEM may be invoked more than once to incrementally add
flow elements to the assemblage. Using this feature an analyst may consider a series of
successively more complex flow system assemblages and their response to specified
8-6
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
excitations.
8.2.3 KINELEM
Interactive species behavior governed by first order kinetics may be accounted for
in the model through the assembly of kinetics elements. A kinetics element will,
typically, model chemical, radio-chemical, or sorption kinetics between a contaminant
species and the immediate environment or other species within a well-mixed zone. As
such they may be associated only with those system nodes that correspond to well-
mixed zones. These elements may be added to the assembly with the following
command/data group that first defines pertinent rate coefficient matrices and then
assigns specific rate coefficient matrices to specific nodes of the system;
KINELEM
K=n1
n2, n3, ...
n4, n5, ...
K=n1
n2, n3, ...
n4, n5, ...
n6 I=n7 K= n8 [GEN=n9]
• • •
END
where; n1 = the kinetics ID of the following rate coefficient matrix,
n2, n3, ...
n4, n5, .. = the rate coefficient matrix for kinetics ID n1
n6 = the kinetics element number
n7 = the well-mixed zone system node number at which the kinetics is
to be applied,
n8 = kinetics ID
n9 = generation increment (default = 1),
The rate coefficient matrices are entered by rows, in species-number order, and
must be defined in terms of all species considered whether or not all of the species are
involved in a given rate expression. Therefore, for a system involving N species all rate
coefficient matrices will be square matrices with N x N terms. If a particular species is
not involved in a given rate expression the terms in the columns and rows
corresponding to this species will simply be zero values.
Element data must be supplied in numerical order. Omitted data is automatically
generated by incrementing the preceding node number by the current generation
8-7
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
8. Command Reference
increment. Generated elements will have the properties of the current element.
Example
Given a system involving three species, say, A, B, and C involved in the following
chemical reactions;
1) a simple reversible reaction
A <-> B
governed by the rate expression;
d[A]
dt
d[B]
dt
d[C]
dt
= -0.45[A] + 0.45[B] + 0.0[C]
= 0.45[B] - 0.45[A] + 0.0[C]
= 0.0[A] + 0.0[B] + 0.0[C]
or, more concisely by the first order rate coefficient matrix;
0.45-0.45 0
-0.45 0.45 0
000
and,
2) a single-step reaction
B -» products (that are of no particular interest)
governed by the rate expression;
d[A]
dt
d[B]
dt
d[C]
dt
= 0.0[A] + 0.0[B] + 0.0[C]
= 0.0[A] - 0.35[B] + 0.0[C]
= 0.0[A] + 0.0[B] + 0.0[C]
or, more concisely by the first order rate coefficient matrix;
8-8
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
8. Command Reference
000
0 0.35 0
000
where the first reaction occurs at system nodes 3,5,7,9 while the second reaction only
occurs at nodes 5 and 7 (e.g., due to the action of a specific catalyst at these nodes)
kinetics elements could be added to the contaminant dispersal system using the
following command/data group;
KINELEM
K=l
0
-0
0
K=
0
0
0
1
4
5
6
.45 -0.
.45
• Q
2
.0
.0
.0
1=3
1=9
1=5
1=7
0.
0.
0.
0.
0.
45 0.0
45 0.0
0 0.0
0 0.
35 0.
0 0.
K=l
K=l
K=2
K=2
0
0
0
G
END
< Kinetics ID 1: A <=> 3
< Kinetics ID 2: B => products
< Kin Elem 1: Node 3: Kinetics ID 1:
< Kin Elems 2,3, & 4: Nodes 5,7,& 9: Kinetics ID 1:
< Kin Elem 5: Node 5: Kinetics ID 2:
< Kin Elem 6: Node 7: Kinetics ID 2:
The rate expressions, above, have been written in terms of all contaminant species,
including the nonreactive species C, to emphasize the manner in which rate
expressions are defined through the use of rate coefficient matrices.
8.2.4 FORM-[W]
In some instances an analyst may wish to examine the details of the mass transport
matrix, [W]. The command FORM-[W] answers this special (and unusual) need. This
command is not required as an interim step to complete any of the analyses options
offered by subsequently defined commands.
The system mass transport matrix, [W], assembled from flow element and reaction
element matrices may be formed with the following command/data group;
FORM-[W] [F=c1c2c3c4]
n1,n2,n3 W=n4
END
where; c1c2c3c4 = FULL or BAND; (default = FULL)
n1,n2,n3 = first element number, last element number, element number
increment of a series of elements with identical mass flow rates
n4 = element total mass flow rate
8-9
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
The matrix may be formed in its full form or compacted form (i.e., only the nonzero
band of the [W] matrix). The system mass transport matrix may be printed or diagrammed
using the intrinsic commands PRINT and DIAGRAM.
8.2.5 STEADY
The response of the system to steady contaminant generation with steady element
mass flow may be computed with the following command/data group;
STEADY
n1,n2,n3 W=n4
n5,n6,n7 CG=n8,n9,...
END
where; n1 ,n2,n3 = first element number, last element number, element number
increment of a series of elements with identical mass flow rates
n4 = element total mass flow rate
n5,n6,n7 = first node, last node, node increment of a series of nodes with
identical excitation
n8,n9,... = contaminant concentration or contaminant generation rate for
each species considered, as appropriate to the boundary
condition of the node/species combination specified with the
FLOWSYS command; (default = 0.0)
Net total mass flow rate at each system node will be reported, but computation will
not be aborted if net mass flow is nonzero. The analyst must assume the responsibility to
check continuity of mass flow from these reported values.
8.2.6 TIMECONS
System time constants, nominal and actual, may be computed with the following
command/data group;
TIMECONS [E=n1]
n2,n3,n4 W=n5
END
where; n1 = optional convergence parameter, epsilon ; (default = machine
precision)
n2,n3,n4 = first element number, last element number, element number
8-10
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
increment of a series of elements with identical mass flow rates
n5 = element total mass flow rate
The nominal time constants are computed for each node as the quotient of the nodal
volumetric mass divided by the total air flow out of a zone. The actual time constants
are computed using an eigenanalysis routine that is a variant of Jacob! iteration
adapted for nonsymmetric matrices [Eberlein et. al. 1971]. It should be noted that the
actual time constants are likely to be very different from the nominal time constants for
systems having well-coupled zones and the eigenanalysis of the flow system matrices is
a time consuming task.
Example
To determine the time constants associated with the building idealization presented
earlier, in the introductory example, the following command/data group would have to
be added to the command/data file.
TIMECONS < (Air Density 1.2E+3 g/m3)
1,2 W=70*1.2E+3 < Supply Ducts
3,6 W=20*1.2E+3 < Infiltration
7,10 W=70*1.2E+3 < Return Loop
11 W=140*1.2E+3 < Main Return Duct
END
8.2.7 Dynamic Analysis
The response of the system, including transients, to general dynamic excitation,
may be computed using the command DYNAMIC. The dynamic solution procedure
used is driven by discrete time histories of excitation and element mass flow rate data
that must first be generated with the commands FLOWDAT and EXCITDAT.
8.2.7.1 FLOWDAT
Discrete time histories of element mass flow rate may be defined, in step-wise
manner, from given element mass flow data, as illustrated below;
Data
Value
' Given Data Value
• Time
TIME(I) TlhE(2) TlhE(3) TIME(4)
Figure 8-1 Arbitrarily Defined Time History Data
8-11
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
8. Command Reference
or, alternatively, discrete time histories of element mass flow data, defined in a step-wise
manner at equal time-step intervals along piece-wise linear segments, may be
generated from given element mass flow data over a time range defined by an initial
time, TJ, a final time, Tf, and a generation time increment, AT, as illustrated below;
Data
Value
—Given Data Value
Generated Value
AT
Time
TIME(1)
TIME(2) TIME(3)
TIME(4)
Figure 8-2 Eaual-Time-Steo-Generated Time History Data
using the following command/data group;
FLOWDAT [T=n1,n2,n3]
TIME=n4
n5,n6,n7 W=n8
TIME=n4
n5,n6,n7 W=n8
[additional TIME data to define the complete excitation time history]
END
where; n1 ,n2,n3 = initial time, final time, time step increment used for the generation
option
n4 = time value for subsequent data subgroups
n5,n6,n7 = first element number, last element number, element number
increment of a series of elements with identical mass flow data
n8 = prescribed element mass flow: (default = 0.0)
If data values n1,n2,n3 are specified, step-wise time histories will be generated from
the given data, along piece-wise linear segments as illustrated in Fig. 7.2 above,
otherwise the given data will be used directly, as illustrated in Fig. 7.1 above.
At least two "TIME" data subgroups must be provided. FLOWDAT writes the
generated time history to the file .WDT so that this data may subsequently be
accessed by the command DYNAMIC.
8-12
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 8. Command Reference
8.2.7.2 EXCITDAT
Discrete time histories of excitation data may be defined in the two ways discussed
above for the FLOWDAT command using the following command/data group;
EXCITDAT [T=n1,n2,n3]
TIME=n4
n5,n6,n7 CG=n8,n9,...
TIME=n4 [additional TIME data to define the complete excitation time history]
n5,n6,n7 CG=n8,n9,...
END
where; n1,n2,n3 = initial time, final time, time step increment used for generation
option
n4 = time value for subsequent data subgroups
n5,n6,n7 = first node, last node, node number increment of a series of
nodes with identical excitation data
n8,n9,... = prescribed contaminant concentration o_r prescribed
contaminant generation rate (as appropriate to node boundary
condition) for each species considered: (default = 0.0)
If data values n1 ,n2,n3 are specified, step-wise time histories will be generated, from
the given data, along piece-wise linear segments as illustrated in Fig. 7.2 above,
otherwise the given data will be used directly, as illustrated in Fig. 7.1 above.
At least two "TIME" data subgroups must be provided. EXCITDAT writes the
generated time history to the file .EDT so that it may subsequently be
accessed by the command DYNAMIC.
8.2.7.3 DYNAMIC
The response of the system to excitation defined by the EXCITDAT command,
using the prescribed element flow data defined by the FLOWDAT command, may be
computed using the following command/data group;
DYNAMIC
T=n1,n2,n3 [A=n4] [Rl=n5] [PS=n6]
n7,n8,n9 IC=n10,n11,...
END
8-13
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
8. Command Reference
where; n1,n2,n3 = initial time, final time, time step increment
n4 = integration parameter, a, where 0.PLT, of concentration response results will be created
with values scaled by the factor n6
n7,n8,n9 = first node, last node, node increment of a series of nodes with
identical data
n10,n11,...= initial nodal concentration for each species in species order;
(default = 0.0)
The response is computed using the. predictor-corrector method presented earlier
[Axley 1987]. With this method the system flow matrix is updated at the discrete times
used to define element flow rate time histories and the system excitation is updated at the
discrete times used to define excitation time histories, as illustrated below;
Row
Pteta
Vail 10
I
,
<
Gene
c\ r\
rl_U
' Flow Time Step
Time
Excitation
Data
Value
4
'
1
...
: Generated by
> : EXCITDAT
Excitation Time Step
Response
Value
Computed by
DYNAMIC
"Time
Integration Time Step
Figure 8-3 Flow and Excitation Driven Dynamic Solution Procedure
The accuracy of the computed response is, therefore, dependent upon the choice of
the flow data time step, the excitation data time step, and the integration time step
chosen by the analyst. Furthermore, the flow data and excitation data time steps may be
nonconstant. The analyst should, therefore, consider investigating the effects of the
choice of these time constants to gain a sense of the error they induce.
8-14
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
8. Command Reference
8.2.7.4 Dynamic Analysis Example
To provide an example of a command/data sequence needed for dynamic analysis
we may consider an extension to the introductory example presented earlier; the
analysis of the dynamic response of the given building system, under conditions of
constant air flows, to a step change in NC>2 generation. Specifically, to consider the
case where the kerosene heater is turned on and then turned off 133 minutes later, the
following command/data group would have to be added to the command/data file used
in the introductory example;
FLOWDAT
TIME=0.0
1,2 W=-70*1.2E+3
3,6 W=20*1.2E+3
7,10 W=70*1.2E+3
11 W=140*1.2E+3
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
9. Example Applications
Examples of the application of CONTAM87 to practical problems of building
contaminant dispersal are presented in this section. The reader will also find a
discussion of the use of the convection-diffusion flow element for both steady state and
dynamic analysis of contaminant dispersal in one-dimerisional flow paths presented in
section 4.3.
9.1 IBR Test House Study
While working at the National Swedish Institute for Building Research (IBR) Kai Siren
developed a program for multi-zone contaminant dispersal analysis, MULTIC, and
applied it to the analysis of the dynamic behavior of the five-room test house maintained
by the IBR [Siren 1986], Figure 9-1. Using data reported by Sir6n the building
idealization shown below, Figure 9-2, was formulated and the (dynamic) decay
response of the idealization to an initial concentration in the bed room, node 2, was
computed (for steady flow conditions) and compared to the results reported by Sir6n. Air
flow rates, zonal volumes, and initial conditions are reported below, Figure 9-2.
Figure 9-1 The IBR Five-Room Test House
9-1
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
13.28
"exterior zone"
9.45
V = 55.6
C(0) = 350
35.67
22.39
V = 13.0
C(0) = 350
V = 36.2
C(0) = 35C
9.45
V . 36.1
C(0) = 2500
42.89^
27.40
V = 34.8
C(0) = 350
0.33
7.80
0.23
1.21
16.70
Figure 9-2 Idealization of IBR Test House
(all flows [=] l/s; volumes [=] m3; initial concentrations [=] ppm)
CONTAM87 Command/Data File The CONTAM87 command/data file used to
complete the analysis is listed below.
* IBR5ZONE
*
*
*
FLOWSYS
N=6 S=l ID=C02
6 BC=C
Swedish IBR Test House: 5 Zone Model
Units m, hr
Analysis by volume rather than mass (i.e.
contaminant density set to unity).
air and
1
2
3
4
5
6
END
V=55
V=36
V=36
V=34
V=13
V=l.
.
.
f
.
B
6
1
2
8
0
OE+06
<
<
<
<
<
<
Node
Node
Node-
Node
Node
Node
1
2
3
4
5
6
FLOWELEM
1
2
3
1=6,
1=6,
1=1,
1
2
3
<"Ext. Zone" Cone. Specified
< Air density set to unity.
Vol. Mass
Vol. Mass
Vol. Mass
Vol. Mass
Vol. Mass
Exterior Vol. Mass
9-2
-------
Indoor Air Quality Model PART II - CONTAM87 USERS MANUAL
Phase III Report 9. Example Applications
4
5
6
7
8
9
10
11
12
13
14
END
*
1=3,1
1=2,3
1=3,5
1=5,3
1=3,4
1=4,3
1=6,5
1=5,6
1=6,3
1-6,4
1=4,6
* STEADY STATE ANALYSIS: C02 generation assumed constant at 80 I/sec.
*
STEADY < Flow rates - m3/sec x 3600 sec/hr
1 W=13.28E-03*3600
2 W=9.45E-03*3600
3 W=35.67E-03*3600
4 W=22.39E-03*3600
5 W=9.45E-03*3600
6 W=19.52E-03*3600
7 W=12.05E-03*3600
8 W=42.89E-03*3600
9 W=27.40E-03*3600
10 W=0.33E-03*3600
11 W=7.80E-03*3600
12 W=0.23E-03*3600
13 W=1.21E-03*3600
14 W=16.70E-03*3600
< < Generation Rates & Specified Concentration
4 CG=80E-03 < Steady Generation [=] m3-C02/hr
6 CG=350E-06 < Exterior Concentration [=] m3-CO2/m3-air
END
*
* DYNAMIC ANALYSIS: Flow steady, C02 generation 80 I/sec for 1.0 hr.
*
FLOWDAT
TIME=0.0
1 W=13.28E-03*3600
2 W=9.45E-03*3600
3 W=35.67E-03*3600
4 W=22.39E-03*3600
5 W=9.45E-03*3600
6 W=19.52E-03*3600
7 W=12.05E-03*3600
8 W=42.89E-03*3600
9 W=27.40E-03*3600
10 W=0.33E-03*3600
11 W=7.80E-03*3600
12 W=0.23E-03*3600
13 W=1.21E-03*3600
14 W=16.70E-03*3600
<
TIME=10.1
1 W=13.28E-03*3600
2 W=9.45E-03*3600
3 W=35.67E-03*3600
4 W=22.39E-03*3600
5 W=9.45E-03*3600
6 W=19.52E-03*3600
7 W=12.05E-03*3600
8 W=42.89E-03*3600
9 W=27.40E-03*3600
9-3
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
W=0.33E-03*3600
W=7.80E-03*3600
W=0.23E-03*3600
W=1.21E-03*3600
W=16.70E-03*3600
10
11
12
13
14
END
EXCITDAT
TIME=0.0
4 CG=80E-03
6 CG-350E-06
TIME=1.0
4 CG=0.0
6 CG=350E-06
TIME=10.1
4 CG=0.0
6 CG-350E-06
END
DYNAMIC
T=0,10,0.1 PS=1.0E+6
1 IC-350E-06
2 IC=2500E-06
3,6,1 IC=350E-06
END
RETURN
< Steady Generation [=] m3-C02/hr
< Exterior Concentration [-] m3-C02/m3-air
< Steady Generation [=] m3-CO2/hr
< Exterior Concentration [=] m3-C02/m3-air
< Steady Generation [=] m3-C02/hr
< Exterior Concentration [=] m3-C02/m3-air
Results The results obtained using CONTAM87 are compared to those using Siren's
program MULTIC below, Figure 9-3. These results are practically identical, as they
should be, as this study provides, in effect, a comparison of two numerical solutions of
the identical system of equations. Nevertheless, the results do indicate that both
programs have numerical procedures that have been correctly coded.
2500 0
2000
Cone.
(ppm)
1500
1000
500
• Node 1
O Node 2
• Node 3
D Node 4
A Node 5
012345678910
Time (hr)
Figure 9-3 Comparison of MULTIC Results (markers) to CONTAM87 Results (lines)
9-4
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
9.2 Carnegie-Mellon Townhouse Study
Borrazzo and his colleagues at Carnegie-Mellon University, Pittsburg,
Pennsylvania, have conducted detailed field investigations of a two-story townhouse
measuring CO, NO, and N02 emissions characteristics of the gas appliances within the
townhouse and the dispersal of these contaminants throughout the townhouse under a
variety of different weather conditions [Borazzo et. al. 1987]. Illustrated in Figure 9-4 is
an idealization of the townhouse and the measured dynamic emission characteristics of
the principal pollutant source, the gas range. The instantaneous emission rate, G(t), is
plotted relative to the steady state value, Gss. The N02 emission characteristics were
more or less constant and are, therefore, not illustrated. N02 is a reactive contaminant
and was modeled as so using the measured reactivity of K=2.4 hr1-
O.U]
G(t)/G
2.0
1.0
0.0(
SS~
•^^
^^H
^MM
(
\
2<
>
) 0.5 Time(hr) 1-
Building Idealization
G(t)/C
1.0
0.5
0.0(
Son
^^^
(55!
1
VHHWt
^<
>
) 0.5 Time(hr) 1.
Emission Characterics
Figure 9-4 Townhouse Building Idealization and Range Emission Characteristics
CONTAM87 Command/Data File The CONTAM87 command/data file used to
complete the analysis for the N02 response is listed below.
* Borrazzo et al's. Townhouse 4-Node, 3-Zone Example
* Units: g, hr, m
*
FLOWSYS
N=4 S=l ID=N02
4 BC=C < Node 4 is exterior node
9-5
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
1,3V=126.5*1.2E+03
4 V=126.5E+09*1.2E+03
END
FLOWELEM
1 1=1,2
2 1=2,1
3 1=2,3
4 1=3,2
5 1=4,2
6 1=2,4
7 1=4,3
8 1=3,4
9 1=4,1
10 1=1,4
END
FLOWDAT
TIME-10.0
1,2,1 W=0.4*126.5*1.2E+03
3,4 W=7.5*126.5*1.2E+3
5,8 W-0.21*126.5*1.2E+03
9,10 W=0.21*126.5*1.2E+03
TIME=20.0
1,2,1
3,4
5,8
9,10
END
KINELEM
K=l
2.4
W-0,
W=7,
4*126.5*1.2E+03
5*126.5*1.2E+3
W=0.21*126.5*1.2E+03
W=0.21*126.5*1.2E+03
< Air Density = 1.2E+03 g/m3
< Exterior Zone Set At Large Value
< ACH bsmnt-to-first
< ACH first-to-second
< 0.21 ACH ext-to-first &
< ACH ext-to-bsmnt
second
< ACH bsmnt-to-first
< ACH first-to-second
< 0.21 ACH ext-to-first & second
< ACH ext-to-bsmnt
< Rxn 1: N02 => products
< Kinelem 1: Node 1: Rxn 1
< Kinelem 2: Node 2: Rxn 1
< Kinelem 2: Node 2: Rxn 1
1 1=1 K=l
2 1=2 K=l
3 1=3 K=l
END
* Transient N02 Emission Model (basis :Gss = 12 (J.g/kJ
EXCITDAT
TIME=10.0
2 CG=0.0068
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
4
END
RETURN
IC=0.0206E-06
<0.013 ppm * (46/28.98)
The CONTAM87 command /data files for the CO and NO analysis would be
identical to the file listed above with the species IDs changed from N02 to CO and NO,
respectively, the KINELEM command removed, and the EXCITDAT and DYNAM
commands replaced with those listed below. Note that in both cases a constant pilot
light generation contribution plus a dynamically varying burner generation contribution
is accounted for.
For CO:
* Transient CO Emission Model (basis:Gss = 98 |J.g/kJ ; Igas •» 150 kJ/min)
EXCITDAT
TIME=10.0
2 CG-0.0415 0.0415
4 CG=0.389E-06
END
DYNAMIC
T=10, 16. 0,0.1 RI=1
IC=0.389E-06
IC=0.389E-06
IC=0.389E-06
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
For NO:
* Transient NO Emission Model (basis:Gss
EXCITDAT
TIME=10.0
2 CG-0.0038
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
exchange rate was assumed to be 0.4 ACH, and all flows were assumed to be constant.
As may be seen, the CO response was under-predicted and the NO response was
over-predicted, but both are within the reported uncertainty of the emission
characteristics (CO: 18% & NO: 6.5%).
[CO]
(ppm)
10
Upstairs
Bsmnl
O Lvng Rm
'~ Model Bsmnt
Kitch n Outside
Model 1st Fir — Model 2nd Fir
11
12
13
Time (hour)
0.55
Q-a--n-n-n—o-n-o-n-a -a-rj—a-n-a~n-a-a-n-
13
Time (hour)
Figure 9-5 Comparison of Computed and Measured CO & NO Response
Although, the measured N02 data is quite suspect, because of its scatter and
negative values, there appears to be some agreement between this data and the
computed response. Importantly, it is noted that the NO2 concentration can fall below
ambient levels out-of-doors due to the reactivity of this contaminant.
data, reporting an interzonal exchange rate of 1.35 ACH. Their method was considered to be very
poorly conditioned, thus, their results unreliable, and, therefore, the interzonal air change rate was
assumed based upon the computed behavior of the townhouse arid past experience.
9-9
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
Inasmuch as the measured data was used to determine the reactivity constant the
agreement here may be an artifice. The basis of determination of the reactivity and the
basis of the computed response are more or less the same as the system behaves,
practically, as a single zone system, therefore the agreement may reflect no more than
this.
[NO2]
(ppm)
11
12
13
Time (hour)
14
15
16
Figure 9-6 Comparison of Computed and Measured NOo Response
^"^**^^""^"^^~ - - * ,__-_-._- , -,« -,„._-,j-r-n^ ._-_n_ - £^^.^__M»*B*_«i^H_
(NOa Reactivity = 2.4 hr -1)
9.3 NBS Office Building Study
Infiltration and ventilation studies of a fifteen story office building are presently being
conducted by members of the Indoor Air Quality and Ventilation Group at NBS. Some of
these studies involve periodic injections of a commonly used tracer gas, SF6, into the
fresh air supply ports of the building HVAC system. Flows in the supply ducts were
measured (with significant uncertainty) by hot-wire anemometer traverse, SF6
concentration time histories were recorded, and outdoor air change rates were
estimated by tracer decay. Using the air flow measurements the upper two floors of this
building were idealized as shown in Figure 9-7.
As indicated by this idealization, fresh air was supplied to each floor through a
ceiling plenum space and exhausted via an exhaust duct to the outside. In Figure 9-8
we compare measured SF6 concentration time histories (measured centrally within the
"space" and at the "exhaust" ports) to computed values of the 15th floor for two supply
flow rates: 100% and 75% of the measured flow. In this case, the agreement between
measured and computed time histories is within the uncertainty of the measured flows
and validation is therefore indicated.
9-10
-------
Indoor Air Quality Model
Phase III Report
PART II - CONTAM87 USERS MANUAL
9. Example Applications
8
rs|
E
7
*
_
3
3
? J
x L
II
T
4j
7
4d
4-
_ . „ . 100% Supply
Exterior "Zone" f1 g^
(1.68 A
15th Ceiling Plenum «3 I
03 4550 c
15th Floor 9 4
• J M .
I...J I...I
0.20 ACH
14th Ceiling Plenum • 5 f
06 4550 Cfi
14th Floor « 6
..... .....
• y y -
0.20 ACH
Flow
fm
OH,
fm
5
Tl
Ra
M
L
2
o
Q
§;
CO
te
Figure 9-7 Idealization of the 14th and 15th Floors of an Office Building
1000
15th "Space"-12/1/86
15th "Exhaust"- 12/1/86
— Model 3 Node 4- 10095
— Model 3 Node 4 - 75S5
100
200
600
700
300 400 500
Time (min)
Hours 9-8 Comparison of Computed and Measured Response for an Office Building
800
9-11
-------
Indoor Air Quality Model
Phase III Report 10. Summary and Directions of Future Work
10. Summary and Directions of Future Work
In the first section of this report we have attempted to give clearer definition to the
emerging field of indoor air quality analysis. It has been argued that the central problem
of indoor air quality analysis is contaminant dispersal analysis and that the related
problems of inverse contaminant dispersal analysis, flow analysis, and thermal
analysis may be thought to serve the needs of contaminant dispersal analysis.
Furthermore, we have suggested that the central problem and these related problems
may be addressed with an integrated set of computational tools based on an element
assembly formulation of the familiar well-mixed zone simplification of the macroscopic
equations of motion for multi-zone building systems of arbitrary complexity. The
CONTAM family of programs is presently under development to provide one
demonstration of this integrated approach; the first two members of the family
CONTAM86 and CONTAM87 are presently available and provide contaminant
dispersal analysis tools.
The noninteractive contaminant dispersal theory presented in the Phase II report of
this project [Axley 1987] has been extended in this report through:
a) a discussion of strategies of forming contaminant dispersal analysis equations for
multi-zone systems involving multiple contaminant species,
b) the introduction of element equations that may be used to model mass transport
phenomena governed by first order kinetics, and
c) through the introduction of element equations that may be used to model the
details of mass transport driven by conduction and diffusion processes in one
dimensional flow paths.
CONTAM87 provides a complete computational implementation of the contaminant
dispersal theory presented earlier and that introduced here. As such, CONTAM87
provides a set of indoor air quality analysis commands1 that are a superset of those
made available in CONTAM86. Future members of the CONTAM family of programs will
provide additional indoor air quality analysis commands superseding or complimenting
those made available by earlier members of the family.
Although it is well recognized that kinetics plays an important role in chemical,
radio-chemical, sorption, and settling processes that affect contaminant dispersal
processes in buildings, the detailed knowledge needed to apply the kinetics analysis
techniques presented here is often not available and actual field or experimental
measured data needed to validate any modeling effort is scarce. The application of the
1 A command, here, is a set of computational procedures that completes a basic indoor air
quality analysis task.
10-1
-------
Indoor Air Quality Model
Phase III Report 10. Summary and Directions of Future Work
kinetics techniques presented here, nominally known as source or sink modeling, has
become an area of emphasis in the direction of our future work.
Although "good" source and sink models are essential to interactive contaminant
dispersal analysis, they introduce a source of uncertainty, as they are inevitably based
upon empirical correlations, that is not a problem in noninteractive contaminant
dispersal analysis. Therefore, while validation of the contaminant dispersal analysis
techniques developed for noninteractive contaminant dispersal involved, primarily, the
verification program logic (i.e., the primary assumption involved was the assumption of
conservation of mass), the validation of techniques developed for interactive
contaminant dispersal analysis will, necessarily, focus on the validity of the specific
source or sink models being employed.
At this time the kinetics of radon decay is well understood and simple models of the
kinetics of formaldehyde emission and N02 reaction are available, yet multi-zone field
measurements needed to validate the use of these models in the multi-zone context are
wanting.
In the development and application of the one dimensional convection-diffusion
element presented in this report, it was recognized that this element provided one means
to model certain classes of imperfectly mixed zones, those zones that behave as-/Ythey
were one dimensional flow systems. Thus this mass transport element could be
considered, also, to be a imperfectly-mixed zone element. Following this train of
thought, a well-mixed zone, whose mass transport behavior is defined by its volumetric
mass, may be thought to be modeled by a well-mixed zone element, rather than being
considered a basic assumption of the underlying theory, and, therefore, the contaminant
dispersal theory presented here may be generalized to remove the restricting
assumption of perfect mixing. Presently, an attempt is being made to recast the element
assembly approach to contaminant dispersal analysis in such a way as to avoid the
limiting assumption of perfect mixing. In this new formulation of the theory the well-mixed
model becomes one special case and a framework is provided for the development of
other imperfectly mixed zone elements.
Two parallel research efforts in the areas of inverse contaminant dispersal analysis
and flow analysis, respectively, are also presently underway, complimenting the
contaminant dispersal analysis work reported here. In the former area integral
formulations of the multi-zone inverse analysis problem have been formulated and used
to develop a new multi-zone tracer gas technique. Field applications of this technique
have proven the technique to be promising. In the flow analysis area, the flow elements
developed during Phase II of this project [Axley 1987] have undergone further
refinement and some additional elements have been formulated. The results of first
applications of these new flow elements have been encouraging.
10-2
-------
Indoor Air Quality Model REFERENCES
Phase III Report
REFERENCES
Axley, J.W., DTAM1: A Discrete Thermal Analysis Method for Building Energy Simulation: Part I Linear
Thermal Systems with DTAM1 Users' Manual. 1985 (presently under review for publication by the
National Bureau of Standards, Building Environment Division, Center for Building Technology)
Axley, James, Indoor Air Quality Modeling: Phase II Report. NBSIR 87-3661, CBT, National Bureau of
Standards, Gaithersburg, MO, Oct., 87
ASHRAE. ASHRAE Handbook 1985 Fundamentals. SI Edition. ASHRAE, Atlanta, GA. , 1985, pp.
33.11-12
Bathe, K.J., Finite Element Procedures in Engineering Analysis. Prentice-Hall, Inc., 1982, pp. 702-706
Borrazzo, J.E., Osborn, J.F., Fortmann, R.C., & Davidson, C.L., " Modeling and Monitoring of CO, NO
and NO2 in a modern Townhouse", Atmospheric Environment, Vol. 21, No. 2, Pergamon Press,
1987, pp. 299-311
Eberlein, P.J. & Boothroyd, J., "Contribution 11/12: Solution to the Eigenproblem by a Norm Reducing
Jacob! Type Method," Handbook for Automatic Computation: Volume II: Linear Algebra, Wilkinson,
J.H. & Reinsch,& Reinsch, C. - editors, Springer-Verlag, 1971
Grot, R.A., Silberstein, S., & Ishiguro, K., Validation of Models for Predicting Formaldehyde
Concentrations in Residences Due to Pressed Wood Products: Phase I. NBSIR 85-3255, CBT,
NBS, Gaithersburg, MD, Sept., 1985
Huebner, K.H. & Thornton, E.A., The Finite Element Method for Engineers. 2nd Edition. John Wiley &
Sons, New York, 1982 . pp. 444-451
Hughes, T.J.R. & Brooks, A., "A Theoretical Framework for Petrov-Galerkin Methods with Discontinuous
Weighting Functions: Application to the Streamline-Upwind Procedure," Finite Elements in Fluids,
Volume 4, Edited by Gallagher et.al.. John Wiley & Sons, 1982, pp. 47-65
Matthews, T.G., Reed, T.J., Daffron, C.R., & Hawthorne, A.R., "Environmental Dependence of
Formaldehyde Emissions from Pressed-Wood Products: Experimental Studies and Modeling,"
Proceedings. 18th International Washington State University Particleboard/Composite Materials
Symposium, pp. 10-23,1984
McNall, P., Walton. G., Silberstein, S., Axley. J., Ishiguro, K., Grot, R., & Kusuda, T., Indoor Air Quality
Modeling Phase I Report: Framework for Development of General Models. NBSIR 85-3265, U. S.
Dept. of Commerce, National Bureau of Standards, October 1986
Moore, J.W. & Pearson, R.G., Kinetics and Mechanism Third Edition. John Wiley & Sons, New York,
1981
Nauman, E.B. & Buff ham, B.A., Mixing in Continuous Flow Systems. John Wiley & Sons, NY, 1983
Nicholas, J.E.. Chemical Kinetics: A Modern Survey of Gas Reactions. John Wiley & Sons, New York,
1976
Ref-1
-------
Indoor Air Quality Model REFERENCES
Phase III Report
Nitschke, I. A., Traynor, G.W., Wadach, J.B., Clarkin, M.E., & Clarke, W.A., Indoor Air Quality.
Infiltration and Ventilation in Residential Buildings. NYSERDA Report 85-10, N.Y. State ERDA,
Albany. N.Y.. Mar. 1985
Sandberg, Mats, "The Multi-Chamber Theory Reconsidered from the Viewpoint of Air Quality Studies,"
Building and Environment, Vol. 19, No. 4, Pergamon Press, Great Britain, 1984 pp. 221-233
Sinden, F.W., "Mufti-Chamber Theory of Air Infiltration," Building and Environment, Vol. 13, Pergamon
Press, Great Britain, 1978 pp. 21-28
Siren, Kai, A Computer Program to Calculate the Concentration Histories and Soma Air Quality Related-
Quantities in a Multi-Chamber System. Helsinki University of Technology, Institute of Energy
Engineering, Otaniemi, 1986
Tezduyar, T.E. & Ganjoo, O.K., "Petrov-Galerkin Formulations with Weighting Functions Dependent
Upon Spatial and Temporal Discretizations: Applications to Transient Convection-Diffusion
Problems," Computer Methods in Applied Mechanics and Engineering, Vol. 59, Elsevier Science
Publishers, North Holland, 1986, pp. 49-71
Traynor, G.W., Allen, J.R., Apte, M.G., Girman, J.R., & Holowell, C.D., "Pollution Emissions from
Portable Kerosene-Fired Space Heaters", Environmental Science & Technology, Vol. 17. June 1983,
pp. 369-371
Walas, S.M., Reaction Kinetics for Chemical Engineers. McGraw-Hill, New York, 1959
Walton, G.N., Estimating Interroom Contaminant Movements. NBSIR 85-3229,U.S. DOC, NBS,
Gaithersburg, MD, August, 1985
Wen, C.Y., & Fan, L.T., Models for Flow Systems and Chemical Reactors. Marcel Dekker, Inc., NY.NY,
1975
Yu, C.C. & Heinrich, J.C., "Petrov-Galerkin Methods for the Time-Dependent Convective Transport
Equation,", International Journal for Numerical Methods in Engineering, Vol. 23, John Wiley & Sons,
1986, pp. 883-901
Zienkiewicz, O.C. & Morgan, K., Finite Elements and Approximation. John Wiley & Sons, NY, 1983
Ref-2
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
APPENDIX
CONTAM87 FORTRAN SOURCE CODE
The FORTRAN?? source code for CONTAM87 is
listed below. In this listing you will note the use of
the compiler directives "INCLUDE". These directives
are commonly available in most compilers but are not
part of the FORTRAN language. They are used to
"include" code stored in separate "include files" that,
here, contain common block data specifications
shared by many subroutines. The contents of these
include files are listed on the last page of this
appendix.
PROGRAM CONTAM
•PRO:CONTAM - BUILDING CONTAMINANT DISPERSAL ANALYSIS PROGRAM
VERSION FY87
Developed by James Axley
Building Environment Division, NBS
Spring 1917
Using;
A) CAL-SAP Library of subroutines developed by Ed Wilson,
u.c. Berkeley
B) Microsoft FORTRAN V2.2 Compiler for Apple Macintosh
For Mac
1. Set logical unit numbers, In SUBROUTINE INITIO, as;
NTR - 9 ; NTW " 9 ; NCMD - 9
2. INCLUDE statements use .INC (I.e.. without ')
3. In SUBROUTINE PROMPT use: WRITE(NTW,'(A,M') STRING
4. In SUBROUTINE EIGEN2 use: WRITE!... A.\) at Section 2.0
C) IBM PC Professional FORTRAN (Ryan-McFarland)
1. Set logical unit numbers. In SUBROUTINE INITIO. as;
NTR = 5 ; NTH <• 6 ; NCMD - 5
2. INCLUDE statements use '.INC' (I.e., with ')
3. In SUBROUTINE PROMPT use: WRITE(NTW.'(A)') STRING
4. In SUBROUTINE EIGEN2 don't use: WRITE!... A) at Section 2.0
Memory for dynamically allocated/defined arrays is located in
vector IA(MTOT) in blank common. To increase or decrease this
area alter the dimension of IA. in the section 0.0 below, set
MTOT, in section 1.0 below, equal to this new dimension, and
recompile the code. As Integers are 4 bytes wide, memory
dedicated to IA(MTOT) Is equal to MTOT*4 bytes.
The number of species is presently limited to MAXSP£°25.
This may be changed by altering MAXSPE In the CNTCOM.INC file.
IMPLICIT REAL'S (A-H.O-Z)
C—0.0 DATA SPECIFICATIONS < COMMON STORAGE
COMMON MTOT,NP,IA<50000)
INCLUDE 'ARYCOM.INC1
INCLUDE 'IOCOM.INC1
INCLUDE 'CMDCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
IF(MODE.EQ.'INTER') CALL PROMPT!' CMND>'//CHAR(7))
CALL FREE
IF(MODE.EQ.'BATCH') CALL FREEHR(NTH)
C
C—3.3 INTERPRET COMMAND LINE
C
C GET COMMAND ( ARRAY NAMES. IF ANY
C
CALL FREECC ' .NCMND, 8.1)
CALL FREEC('A',M1U),4.7)
C
c INTRINSIC COMMANDS
C
IFKNNCMND.EQ.'H'I .OR. (NNCMND.EO.'HELP')) THEN
IF(MODE.ZQ.'BATCH') THEN
CALL ALERT('ERROR: Command not defined In BATCH mode.',
+ 'S'.'S')
CALL RETRN
ELSE
CALL HELP
ENDIF
ELSEIFI(NNCMND.EO.'ECHO-ON').OR.(NNCMND.EO.'ON')) THEN
ECHO - .TRUE.
ELSEIF((NNCMND.EO.'ECHO-OFF').OR.(NNCMND.EO.'OFF')) THEN
ECHO = .FALSE.
ELSEIF((NNCKND.EQ.*L').OR.(NNCHND.EQ.'LIST')) THEN
IF(MODE.EQ.'BATCH11 THEN
CALL ALERT('ERROR: Command not defined In BATCH mode.',
+ 'S'.'S')
CALL RETRN
ELSE
CALL LIST
ENDIF
ELSEIF((NKCMND.EQ.'P').OR.(NNCMND.EQ.'PRINT')) THEN
CALL PRINT
ELSEIF((NNCKND.EO.'D').OR.(NNCMND.EQ.'DIAGRAM')I THEN
CALL DII.GRM
ELSEIF(NNCMND.EQ.'PAUSE') THEN
PAUSE ' "" PAUSE: Enter to continue.'
ELSEIF((NHCMNO.EQ.'S').OR.(NNCMND.EO.'SUBMIT')) THEN
IF(MODE.EQ.'BATCH') THEN
CALL ALERT('ERROR: Command not defined In BATCH mode.',
+ 'S'.'S')
CALL 11ETRN
ELSE
CALL .'SUBMIT
ENDIF
ELSEIFKNHCMND.EO.'R'I.OR. (NNCMND.EQ.'RETURN') I THEN
IF(MODE.EQ.'INTER') THEN
WRITE(NTW,2320)
ELSE
CALL KETRN
ENDIF
2320 FORMAT!' **•* ERROR: Command not defined in INTERACTIVE mode.1)
ELSEIF((NNCMND.EQ.'0').OR.(NNCMND.EQ.'QUIT1)) THEN
CLOSE(NOT)
STOP
CONTAM COMMANDS
ELSEIF(NNCMND.EQ.'FLOWSYS'I THEN
CALL FLOSYS
ELSEIF(NNCMND.EQ.'FLOWELEM') THEN
CALL FLOELM
C--1.0 INITIALIZE INTERNAL VARIABLES
c
MTOT •> 50000
CALL INITARIMTOT)
CALL INITIO
CALL INITCN
ERR=.FALSE.
C—2.0 WRITE BANNER
CALL BANNER(NTW)
CALL BANNER(NOT)
WRITE(NOT.2200) (FNAME(1:LFNAME)//'.OUT1)
2200 FORMAT!/1 •"•= RESULTS OUTPUT FILE: ', (A) I
ELSEIF(NNCMND.EO.'KINELEM') THEN
CALL XINELM
ELSEIFWCMND.EQ. 'FORM-[H| ') THEN
CALL FORMFO
ELSEIF(NNCMND.EQ.'STEADY') THEN
CALL STEADY
ELSEIF(NMCMND.EO.'TIMECONS') THEN
CALL T1IMCON
ELSEIF (N1ICMND.EQ. 'FLOWDAT'I THEN
CALL F1.0DAT
ELSEIF(NMCMND.EQ.'EXCITDAT') THEN
CALL EXCDAT
C--3.-0 COMMAND PROCESSOR LOOP
C
C 3.1 CHECK BLANK COMMON STORAGE
C
30 NSTOR -
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
SUBROUTINE INITMUMTOT)
C—SUB:INITAR - INITIALIZES DYNAMIC ARRAY MANAGER VARIABLES
C IN BLANK COMMON AND LABELED COMMON /ARYCOM/
INCLUDE 'ARYCOM.INC'
DUMA - 0
NEXT = 1
IDIR " MTOT
IP<1) - 4
IP (2) = 8 '
IPO) - 1
RETURN
END
SUBROUTINE INITIO
C—SUBUNITIO - INITIALIZES LABELED COMMON /ICCOM/
C OPENS DEFAULT RESULTS OUTPUT FILE
INCLUDE 'IOCOM.INC1
NTR = 5
NTH = «
NCMD « 5
NIN - 10
NOT = 11
NPLT - 12
ND1 = 13
ND2 - 14
ND3 - 15
ND4 - 16
FNAME = 'CONTAM'
LFNAME = 6
EXT = ' '
CALL NOPENINOT,(FNAME(1:LFNAME)//'.OUT'),'FORMATTED')
MODE - 'INTER'
ECHO - .TRUE.
RETURN
END
SUBROUTINE INITCN
C—SUBlINITCN - INITIALIZES CONTAM LABELED COMMON /CNTCOM/
INCLUDE 'CNTCOM.INC'
C INITIALIZE CONTAM CONTROL VARIABLES
IHITAR C
C-
(H)ELP
List available Intrinsic commands.',/,
NSNOD
NSSPE
NSEQ
MS BAN
NFELM
NKINEL
C INITIALIZE POINTERS
MPV =0
MPVM = 0
MPF =0
MFC =0
MPE -0
MPKSEQ - 0
MPNE - 0
MPEFF = 0
KPDIFF= 0
MPCENR- 0
DO 10 N-1,9
10 MPKIK(N) = 0
MPTEMP - 0
C INITIALIZE OTHERS
EP - l.OD-16
RETURN
END
INCLUDE 'IOCOM.INC'
HRITE(NTN,2000)
HELP LIST ~
2000 FORMAT!/, ' -
(H)ELP
ECHO-(ON)
ECHO-(OFF)
(L)IST
(P)RINT A-
(D)IAGRAM A-
INTRINSIC COMMANDS',//,
List available intrinsic commands.', /,
Echo results to screen.',/,
Do not echo results to screen.'./.
List the directory of all arrays.',/,
Print array named .'./.
Diagram array named .',/,
(S)UBMIT F*> Read connands from batch .',
(R)ETURN Return to Interactive mode.',/,
(Q)UIT Quit program.'/)
SUBROUTINE LIST
C—SUB:LIST - LIST DIRECTORY OF ALL ARRAYS IN BLANK COMMON
C
C--HELP LIST
C
C .' (L)IST List the directory of all arrays.',/,
COMMON HTOT.NP.IAU)
INCLUDE 'AXYCOM.INC'
INCLUDE 'ICCOM.INC'
CHARACTER*! NAM(4),LOC(4,2).TYPE(9,3),STOR(13,2)
CHARACTER*! CHK
DATA TYPE/'I'.'N'
1 'R'.'E'
2 'C'.'H'
'T','E','G'. '£', 'R'.
•A'.'L',' '.' '.' ',
'A','R'.'A','C','T',
E'.'R'/
DATA LOC/'C','0'.'R'.'E'
SUBROUTINE BANNER(LUN)
C—SUB: BANNER - WRITES PROGRAM BANNER TO LOGICAL UNIT LUN
COMMON MTOT,NP,IA(1)
HRITE(LUN,2000) MTOT
2000 FORMAT!//, IX, 78UH-)./.
.'I C 0 N T A M 9 T.T79. 'I'./.
.' I Contaminant Dispersal Analysis for Building Systems'
.,T79,'|'./,
.' I Version 4/87 - Jim Axley - NBS',
.179. 'I './.IX.TBdH-) ,/,«5X. 'MTOT:',19)
RETURN
END
C
C INTRINSIC COMMANDS
C
C
C
C
SUBROUTINE HELP
C—SUB: HELP - PROVIDES ON-SCREEN HELP
C
C—HELP LIST .
C
DATA STOR/'S','E'.'Q'.'U'.'E'.'N'.'T','I1.'A'.'L'
1 'D','I','R'.'E','C1.'T'.' '.'A1,'C','C'
C
C LIST DIRECTORY OF ALL ARRAYS IN DATA BASE
IF(NUMA.EQ.O) GO TO 900
C
C WRITE HEADER FOR SCREEN LISTING OF FILE DATA
WRITE(NTH,1000)
C
C START COUNT OF LINES TO SCREEN
IL - 5
C
1C = IDIR
DO 100 I-l.NUMA
IL - IL t 1
1LOC = 1
1ST » 0
IA6 • IAUC+6)
IA7 •• IA(IC*7)
IA9 * IAdCn-91
C CHECK FOR LOCATION AND STORAGE TYPE
IFUA9.GT.O) ILOC-2
IFIIA7.LT.O) ILOC-2
IFIIA7.EQ.-1) 1ST-1
IFdA7.EQ.-2l IST-2
IFIIA9.GT.O) IST-3
IPN = 1C - 1
DO 10 J=1.4
IPN = IPN + 1
10 NAM(J) - CHAR(IAIIPN) I
C WRITE DATA TO TERMINAL
IF(IST.EQ.O) WRITE(NTH,1100) (NAM(J),J-l.4),
* IA(IC+4),IA(IC*5). (TYPE(K.IA6),K=1.9). •
* (LOC(L.ILOC).L*1,4)
C
IF(IST.EO.l) WRITE(NTW.1100) (NAM(J) . J-l, 4) ,
* IA(ICt4).IA(IC+5),(TYPE(K,IA6),K=1.9),
• (LOC(L.ILOC).L-1,4). , (LOC (L, ILOC), L-1,4),
• (STOR(M,2),M-l,13)
C
1C - 1C + 10
C CHECK FOR NUMBER OF LINES PRINTED
IFdL.LT.20l GO TO 100
IFII.EO.NUMA) GO TO 100
CALL PROMPT!' ** Do you want more ? (Y/N) ')
READ(NTR,2200)
IFUCHX.EO. 'n'l .OR. (CHK.EQ. 'N')) GO TO 900
IL - 0
WRITE(NTW.2000]
100 CONTINUE
C
900 RETURN
C
1000 FORMAT!' —— LIST: ARRAY LIST',//,
* • Name',2X.'Number1.2X.'Number',5X,'Data'.5X.
" 'Location',SX,'Storage',/,8X,'Rows',2X,
* 'Columns',5X.'Type*,19X,'Type',/)
1100 FORMAT(IX,4A1.2X,I4.4X.I4.5X.9A1,4X.4A1,4X,13A1)
Append -2
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
1200 FORMAT (IX, 4A1, ' NI=',I4,' NR=',I4.' NC=',14,5X.4A1,4X.13A1)
1300 FORMAT(IX.4A1.3X,-RECORD LENGTH = '. 16,7X,4A1,4X,13A1)
2000 FORMAT!)
2200 FORMAT (1 AD
2000 FORMAT!/' COLi ='.6112)
2001 FORMAT(' ROM',14,6E12.S)
2002 FORMAT!' ROH',14.6F12.S)
END
• SUBROUTINE PRINT
C—SUB:PRINT - COMMAND TO "PRINT" ARRAY TO RESULTS OUTPUT FILE
C
C—HELP LIST
C
C .' (P)RINT A- Print array named .'./.
' PRINT C
SUBROUTINE CPRT(C,NR,NO
C—SUB: CPRT - PRINTS CHARACTER*! ARRAY TO RESULTS OUTPUT FILE
COMMON MTOT.NP.IA(l)
INCLUDE 'ARYCOM.INC'
INCLUDE 'IOCOM.INC'
INCLUDE 'CMDCOM.INC'
CHARACTER MA*4
EQUIVALENCE (MA.Ml(D)
C PRINT OF REAL OR INTEGER ARRAY
CALL PROMH(l)
C LOCATE MATRIX TO BE PRINTED
IF(ECHO) WRITE(NTH,2000) Ml
WRITE(NOT, 2000) Ml
CALL LOCATE (Ml. NA.NR.NO
IF(NA.EO.O) THEN
CALL ALERTCERROR: Array '//MA//' does not exist.'.'S'
CALL ABORT
RETURN
ELSEIF(NA.LT.O) THEN
CALL ALERTCERROR: Array '//MA//' Is out of core.','!'
CALL ABORT
RETURN
ELSE
IFWP.EQ.l) CALL IPRT(IA(NA),NR,NC)
IF(NP.EQ.2) CALL RPRT (IA (NA) , NR.NC)
IFOJP.EQ.3) CALL CPRT (lA(NA).NR. NCI
ENDIF
'$')
'$')
RETURN
2000 FORMAT!/'
END
PRINT OF ARRAY "',4A1,'"')
SUBROUTINE IPRT(N,NR.NC)
C—SUB: IPRT - PRINTS INTEGER ARRAY TO RESULTS OUTPUT FILE
DIMENSION N(NR.NC)
INCLUDE 'IOCOM.INC'
NUMC = 14
DO 100 1=1.NC,NUMC
IN » I + NUMC - 1
IF(IN.GT.NC) IN = NC
WRITE(NOT.2000) Diagram array named .'./,
COMMON MTOT.NP.IA(l)
INCLUDE 'IOCOM.INC'
INCLUDE 'CMDCOM.INC'
CHARACTER MA*4
EQUIVALENCE (MA. Ml(1))
IPRT C PRINT OF REAL OR INTEGER ARRAY
CALL PROMH(l)
c LOCATE MATRIX TO BE PRINTED
IFIECHO) HRITE (NTH. 2000) Ml
HRITE(NOT,2000) Ml
CALL LOCATEIMl.NA.NR.NC)
IF(NA.EQ.O) THEN
CALL ALERTCERBOR: Array '//MA//1 does not exist.','S','S')
CALL ABORT
RETURN
ELSEIF(NA.LT.O) THEN
CALL ALERTCERROR: Array '//MA//' Is out of core.', ' S'. ' S'l
CALL ABORT
RETURN
ELSEIFINP.EQ.3) THEN
CALL ALERTCERROR: Array '//MA//' Is a character array.'.
+ 'S'.'S'I
CALL ABORT
RETURN
ELSE
IF(NP.EQ.l) CALL IDIAGRIIAINA).NR.NC)
IFINP.F.Q.2) CALL RDIAGR (IAINA) ,NR,NC)
ENDIF
RETURN
2000 FORMAT!/'
END
> DIAGRAM OF ARRAY
SUBROUTINE IDIAGR(N.NR.NC)
C--SUB: IDIAGR - "DIAGRAMS" INTEGER ARRAY TO RESULTS OUTPUT FILE
INTEGER N(NR.NC)
CHARACTER*! ICON (36)
INCLUDE 'IOCOM.INC1
C DIAGRAM INTEGER ARRAY
NUMC = 36
DO 200 1=1.NC,NUMC
IN = I » NUMC - 1
IFIIN.GT.NC) IN = NC
WRITE(NOT.2000) (INT(K/10).K-I,IN)
WRITE(NOT,2010) ((K-INT(X/10)*10),K=I,IN)
IFIECHO) HRITE(NTH.2000) (INT(K/10).K=I.IN)
IFIECHO) HRITE(NTH.2010) (
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
SUBROUTINE RDIAGR(A,NR,NC>
C—SUB: RDIAGR - 'DIAGRAMS- REAL ARRAY TO RESULTS OUTPUT FILE
REAL'S A(NR.NC)
CHARACTER*! ICON(36)
INCLUDE 'IOCOH.INC'
C DIAGRAM INTEGER ARRAY
NUMC - 36
DO 200 1=1.NC. NUMC
IN - I + NUMC - 1
IF (IN.GT.NO IN - NC
WRITE(NOT,2000) (INT(K/10) ,K-I. IN)
WRITE(NOT.2010) ((K-INT(K/10)«10),K-I.IN)
IF(ECHO) WRITE(NTH,2000) (INT(K/10),K=I,IN)
IF(ECHO) WRITE(NTS,2010) ((K-INT(K/10)*10) , K-I, IN)
DO 200 J-l.NR
DO 100 K-I.IN
ICON(K) - '*'
IF(AIJ.K) .EO. O.ODO) ICON(K) - ' '
100 CONTINUE
WRITE(NOT,2020) J, (ICON (X),K-I.IN)
IF(ECHO) WRITE(NTW.2020) J.(ICON(K),K=I,IN)
200 CONTINUE
C
RETURN
2000 FORMAT(/• COLI ••,36(IX,II))
2010 FORMAT (7X. 36 (IX. ID)
2020 FORMAT!' ROW'. 13.36(IX.Al))
END
SUBROUTINE SUBMIT
C—SUB: SUBMIT - SWITCHES TO BATCH MODE AND OPENS BATCH COMMAND FILE
C
C—HELP LIST
C
C .' (S)UBMIT F- Read commands from batch .',/,
C — SUB
C
C
C — BEL
C
C
C
C
C
C
C
C
C
SUBROUTINE FLOSYS
FLOSYS - COMMAND TO READ < PROCESS FLOW SYSTEM CONTROL VARIABLES
ESTABLISHES FLOW SYSTEM EQUATION NUMBERS t B.C.
FLONSYS N=nl
N-nl S=n2 ID-cl,c2, .
n3,n4,n5 BC-C3
n3,n4,n5 V=n6
END'!//,
Flowsystem control variables.'./.
. nl ** no. flow nodes; n2- no. species
cl,c2,... species IDs (4 chars)',/.
c3 » boundary condition; G or C'./,
n6 •* nodal volumetric mass',/,
COMMON MTOT.NP.IA(l)
INCLUDE 'IOCOH.INC1
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
INTEGER IJK<3)
EXTERNAL BCDATO.VDATO
ERR * .FALSE.
C
C—1.0 GET NUMBER OF FLOW SYSTEM NODES, NUMBER OF SPECIES, f SPECIES IDS
C
IF (ECHO.OR.(MODE.EO.'INTER')) WRITE(NTW.2100)
WRITE(NOT,2100)
2100 FORMAT (/.' — FLOWSYS: FLOW SYSTEM CHARACTERISTICS')
IF(NSNOO.NE.O) CALL RESET
IF(MODE.EQ.'INTER') CALL PROMPT|'DATA>')
CALL FREE
IF (MODE. EQ.'BATCH') CALL FREEWRtNTW)
INCLUDE 'IOCOM.INC'
LOGICAL FOUND
CALL FREEC('F',FNAME.12.1)
INQUIRE(FILE-FNAME(1:LENTRM(FNAME)),EXIST-FOUND)
IF(FOUND) THEN
MODE = 'BATCH'
NCMD » NIN
LFNAME - LENTRM(FNAME)
WRITE(NTW,2010) FNAME
WRITE(NOT,*)
WRITE(NOT,2010) FNAME
2010 FORMAT(• •«•• CONTAM set to BATCH mode uslno file: ',A)
OPEN(NCMD,FILE=FNAME(1:LFNAHE).STATUS-'OLD')
REWIND NCMD
CLOSE(NOT)
CALL NOPENINOT,(FNAME(1:LFNAME)//'.OUT1),'FORMATTED')
CALL BANNER(NOT)
WRITE(NOT.2020) (FNAME(1:LFNAME)//'.OUT')
2020 FORMAT(/' —— RESULTS OUTPUT FILE: '.(A))
ELSE
WRITE(NTW,2030)
2030 FORMAT!' — NOTE: Submit file not found.')
CALL ABORT
ENDIF
RETURN
END
SUBROUTINE RETRN
C—SUB:RETRN - RETURNS TO INTERACTIVE MODE
C
C—HELP LIST
C
C .' (R)ETURN Return to Interactive mode.'./,
C—1.1 NUMBER OF FLOW NODES
CALL FREEI('N',NSNOD,1)
CALL CKIZERCthe number of flow system nodes'.NSNOD.l,'GT',ERR)
IF(ERR) GO TO 999
IF(ECHO) WRITE(NTW,2120) NSNOD
WRITE(NOT,2120) NSNOD
2120 FORMAT!/,1 Number of flow system nodes ',15)
C 1.2 NUMBER OF SPECIES
CALL FREEH'S',NSSPE, 1)
CALL CXIRNGCthe nunfcer of contaminant species'.NSSPE, 1,
•fO, 'LTLE' .MAXSPE.ERR)
IF(ERR) GO TO 999
IF(ECHO) WRITE(NTW,2140) NSSPE
WRITE(NOT.2140) NSSPE
2140 FORMAT(' Number of contaminant species ....'.15)
C 1.3 SPECIES IDS
CALL ZEROC(SID,4,NSSPE)
CALL GETIDS(SID,NSSPE,'Contaminant species IDs:')
C
C—2.0 DEFINE XSEQ ARRAY AND NUMBER EQUATIONS IN (NODE.SPECIES) ORDER
C
IF (ECHO.OR.(MODE.EO.'INTER')) WRITE(NTW,2150)
WRITE(NOT,2150)
2150 FORMAT!/,' =• Boundary Conditions and Equation Nunfcers')
NSEQ - NSNOD'NSSPE
CALL DELETE('XSEQ')
CALL DEF1NI('KSEQ',MPKSEQ,NSNOD,NSSPE)
CALL EQNUMdA (MPKSEQ),NSNOD.NSSPE)
INCLUDE 'IOCOM.INC'
WRITE(NOT, *)
WRITE(NOT. 2010)
CLOSE(NCMD)
CLOSE (NOT)
FNAME - 'CONTAM'
LFNAME - 6
OPEN (NOT. FILE" (FNAME (1: LFNAME) //'.OUT') , STATUS-'OLD',
+FORM-'FORMATTED')
CALL APPEND(NOT)
NCMD - NTR
MODE = 'INTER1
WRITE(NTW, 2010)
WRITE(NOT, •)
WRITE(NOT, 2010)
2010 FORMAT!' •••* CONTAM returned to INTERACTIVE node.')
RETURN
END
CONTAM COMMAND
C—3.0 PROCESS BOUNDARY CONDITION DATA < REPORT EQUATION NUMBERS
C
CALL DATGEN(BCDATO,NSNOD,ERR)
IF(ERR) GO TO 999
CALL RPRTX(IA(MPKSEQ),SID,NSNOD.NSSPE)
C
C--4.0 GET NODAL VOLUMETRIC MASSES
C
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTW,2400)
WRITE(NOT,2400)
2400 FORMAT)/,' = Nodal Volumetric Mass')
CALL DELETE('V ')
CALL DEFINRCV ', MPV.NSNOD, 1)
CALL ZEROR(IA(MPV),NSNOD.l)
CALL DATGEN(VDATO,NSNOD,ERR)
IF(ERR) GO TO 999
CALL RPRTNOIIA(MPV),NSNOD,'Node')
C
C--5.0 ORDERLY COMPLETION: SXIP TO "END"
C
-C 500 IF(EOC) RETURN
C IFIMODE.EQ.'INTER') CALL PROMPTCEND?>')
C CALL FREE
C GO TO 500
Append -4
-------
Indoor Air Quality Model
Phase III Report
CONTAIM87 FORTRAN
APPENDIX
77 Source Code
C—9.0 ABORT IF ERR
C
999 IF(ERR) THEN
CALL DELETE CKSEQ')
CALL DELETECV ')
NSNOD " 0
NSSPE = 0
NPKSEO = 0
HPV - 0
ERR - .FALSE.
CALL ABORT
RETURN
ENDIF
END
SUBROUTINE EONUM(KSEQ.NSNOD,NSSPE)
C--SUB:EQNUM - ESTABLISHES EQUATION NUMBERS
INTEGER KSEQ (NSNCO,NSSPE)
NN - 0
DO 10 N=l,NSNOD
DO 10 M»l,NSSPE
NN = NN+1
10 KSEQ(N,H) » NN
RETURN
END
SUBROUTINE VDATO(N.ERR)
C~SUB:VDATO - CALLS VDAT1 PASSING ARRAYS
C
COMMON MTCT,NP,IA(1)
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
CALL VDAT1(IA(MPV),NSNOD,N,ERR)
RETURN
END
SUBROUTINE BCDATO(N,ERR)
C—SUB:BCDATO - CALLS BCDAT1 PASSING TEMPORARY ARRAY
C
C POINTER—ARRAY
C MPBCBC (NSSPE) *1 : TEMPORARY STORAGE OF BC CODES
COMMON MTOT.NP.IAIl)
INCLUDE 'lOCOH.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
CALL BCDAT1 (LMMPKSEQ) .CDATA.N, NSNOD,NSSPE. ERR)
RETURN
END
SUBROUTINE BCDAT1(KSEQ,BC,N.NSNOD.NSSPE,ERR)
C—SUB: BCDAT1 - PROCESSES BC DATA
INCLUDE 'IOCOM.INC1
INTEGER KSEQ (NSNOD, NSSPE)
CHARACTER BC(NSSPE)*!, BCH'l
LOGICAL ERR
CALL FREECCC',BC(1) .1. NSSPE)
DO 10 H-l,NSSPE
BCM - BC(M)
IFI(BCM.NE.'C').AND.(BCM.NE.'G1).AND.(BCM.NE.' ')) THEN
WRITE(NTH.2000) BCM.N.M
WRITE(NOT,*)
WRITE(NOT,2000) BCM,N,M
ERR « .TRUE.
RETURN
2000 FORMAT(' *•** ERROR: Boundary condition "',A1,'- not available.',
+/, • Node:',14,' Speclea:',I4)
ELSEIF(BCM.EO.'C') THEN
XSEQ(N,M) = -KSEQ(N.M)
10 ENDIF
RETURN
END
c RPRTK
SUBROUTINE RPRTK(KSEQ,SID,NSNOD,NSSPE)
C—SUB:RPRTK - REPORTS SYSTEM EQUATION NUMBERS STORED IN ARRAY
C KSEQ(NSNOD,NSSPE)
INCLUDE 'IOCOM.INC'
INTEGER KSEQ(NSNOD,NSSPE)
CHARACTER SID(NSSPE)•<
IF(ECHO) WRITE(NTH.2000) (SID(M),M=1,NSSPE)
WRITE(NOT,2000) (SID(M),H-l,NSSPE)
DO 10 H-l.NSNOD
IF(ECHO) WRITE(NTH,2010) N, (KSEQ(N, M),M=1,NSSPE)
HRITEINOT.2010) N,(KSEQIN.M),M-1.NSSPE)
10 CONTINUE
RETURN
2000 FORMAT!/,
.6X,'Neg. Eqtn-f •* concentration-prescribed (Independent DOF).',/.
.fix,'Pos. Eqtn-l *> generation-prescribed (dependent DOF).',//,
.13X,'Species ID:',/.
.EX,'Node1,10I3X.A4))
2010 FORMAT16X,14,10<3X,14)>
END
C VDATO
SUBROUTINE VDAT1(V,NSNOD,N.ERR)
C—SUB:VDAT1 - READS NODE VOLUMETRIC MASS DATA
C
• INCLUDE 'IOCOM.INC'
REAL'S V(NSNOD),VDAT
LOGICAL ERR
CALL FREEIM'V'.VDAT, 1)
CALL CKRZURCnodal volumetric mass',VDAT. 1. 'GE' .ERR)
IF(ERR) RETURN
V(N) • VDAT
RETURN
END
C .................. i. mm. .... HI.. ............. ............................ ...... . ..... i-j.i i. n FLOELM
SUBROUTINE FLOELM
C— SUB: FLOELM - COMMAND TO READ t WRITE FLOW ELEMENT DATA TO FILE *.FEL
C
C--HELP LIST ------------------------------------------------------------
C
FLOWELEM Flow elements: Simple or Conv-dlf f . ' , /.
nl T-SIMP I=n2,n3 E=n4'./.
or',/.
L=i>6 D=n7,n8,... G=n9,nlO, . . . F-nll',/,
nl = elera. number; n2,n3 a end nodes ',/,
n4 «" filter eff.; nS = mass/length;',/,
n6 » length; n7,n8,...= disp. coef.',/,
n9.nlO = generation; nil » upwind fact.')
nl T-CNDF I=n2,n3
END
COMMON MTOT.NP.IA(l)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL F.RR, FOUND
EXTERNAL FLOELO
C VARIABLE DESCRIPTION
C ERR ERROR FLAG
C FOUND FILE FOUND FLAG
C FLOELO SUB. TO READ I WRITE FLOW ELEM DATA
C NENOD NUMBER OF ELEMENT NODES
C NESTRT FIRST ELEMENT NUMBER
C—0.0 INITIALIZATION
C
ERR = .FALSE.
NENOD - 2
IFIECHO.OR.(MODE.EQ.'INTER')) WRITE(NTW,2000)
WRITE(NOr.2000)
2000 FORMAT!/,' — FLOWELEM: FLOW ELEMENTS')
C
C--1.0 CHECK TO SEE IF SYSTEM NODES I EQUATION NUMBERS ARE DEFINED
C
CALL CKSYSd.ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C—2.0 OPEN .FEL
C
INQUIRE (FILE-FNAME(1:LFNAME)//'.FEL',EXIST-FOUND)
IF((.NOV.FOUND).OR.(NFELM.EQ.0)) THEN
CALL MOPEN(ND1,(FNAME(1:LFNAME)//'.FEL'),'UNFORMATTED1)
ELSEIF(FOUND) THEN
WRITE(NTW,2200)
WRITE(NOT,•)
WRITE(NOT, 2200)
2200 FORMAT(
+ ' *" NOTE: Additional flow elements being added to system.')
OPEN (ND1, FILE=(FNAMEU: LFNAME)//1. FEL1 ),STATUS=' OLD',
» FORM"'UNFORMATTED')
CALL APPEND(ND1) ~
ENDIF
C
C—3
C
0 DEFINE TEMPORARY ARRAYS; GENERATE ELEMENT DATA; REPORT BANDWIDTH
CALL DELETE('EFF ')
CALL DELETECDIFF')
CALL DELETECGENR')
CALL DEFINRCGENR' .MPCENR.NSSPE, 1)
CALL DEFINRCDIFF'.MPDIFF,NSSPE,1)
CALL OEFINRCEFF '. MPEFF. NSSPE, 1)
NESTRT = NFELM+1
CALL ELGEN(FLOELO,NENOD.NESTRT,NSNOD,ERR)
Append -5
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
CALL DELETE('EFF ')
CALL DELETE('DIFF') C- READ ELEMENT DATA
CALL DELETE('GENR') CALL ZEROR(EFF,NSSPE.1)
CALL FREERCE'.EFFU),NSSPE)
IF(ERR) THEN CALL CXRZERCfilter offIclency'.EFF.NSSPE,'GE'.ERR)
CALL ABORT IF (ERR) GO TO 999
RETURN
ENDIF C UPDATE SYSTEM BANDWIDTH
DO 30 N=l,NSSPE
IF(ECHO) WRITE(NTH.2300) MSBAN LS(l) = N
WRITE(NOT,") 30 CALL ELBAN(KSEQ.NSNOD,NSSPE,LN.NENOD.LS.NESPE.MSBAN)
WRITE(NOT.2300) MSBAN
2300 FORMAT!' •« NOTE: Current system bandwidth la:', IS) C WRITE ELEM. DATA TO ND1
C WRITE(ND1) TYPE
C—4.0 ORDERLY COMPLETION: CLOSE FILE ND1; SKIP TO 'END- WRITE(NDl) LN(1), LN(2), (EFF (N) ,11-1, NSSPE)
C
CLOSE(ND1) C UPDATE ELEMENT COUNT
NFELM - NEL
IF(HCOE.EQ.'INTER') RETURN
500 IF(EOC) RETURN C REPORT ELEMENT DATA
CALL FREE WRITE(NOT,2030) NEL,TYPE, LN (1), LN(2), SID(1),EFF(1)
GO TO 500 IF(NSSPE.GE.2) WRITE(NOT,2032) (SID(N),EFF(N),N=2.NSSPE)
IF(ECHO) THEN
END WRITE(NTW,2030) NEL,TYPE, LN(1), LN (2), SID(l), EFF(l)
IF(NSSPE.GE.2) WRITE(NTW. 2032) (SID (N) , EFF (N) .N-2,NSSPE)
SUBROUTINE FLOELO(NEL,LN,ERR) 2030 FORMAT(2X,14,1X.A4,214,2X.A4,1X.G11.4)
C--SUB:FLOELO - CALLS FLOEL1 PASSING ARRAYS 2032 FORMAT!(21X.A4,1X.G11.4) I
COMMON MTOT,NP,IA(1) C
C—4.0 CONVECTION-DIFFUSION ELEMENTS
INCLUDE 'CNTCOM.INC' C
ELSEIF(TYPE.EO.'CNDF') THEN
LOGICAL ERR
INTEGER LN(2) C READ ELEMENT DATA
CALL FLOEL1(IA(MPKSEQ),IA(MP£FF).IA(MPDIFF).IA(MPGENR),NEL,LN,ERR) MASSL - O.ODO
CALL FREEH CM', MASSL, 1)
RETURN CALL CXRZER('mass/lengtn',MASSL,1,'GE'. ERR)
END IF(ERR) GO TO 999
c FLOEL1 LENGTH - O.ODO
SUBROUTINE FLOEL1(KSEO.EFF.DIFF,GENR,NEL,LN.ERR) CALL FREER)'L',LENGTH,1)
C—SUB:FLOEL1 - READS FLOW ELEMENT PROPERTY DATA, CALL CXRZER('flow passage length'.LENGTH,1,'GT',ERR)
C UPDATES SYSTEM BANDWIDTH MSBAN, IF(ERR) GO TO 999
C WRITES FLOW ELEMENT DATA TO LOGICAL UNIT ND1, AND
C REPORTS ELEMENT DATA TO RESULTS OUTPUT FILE CALL ZEROR(DIFF,NSSPE,1)
CALL FREERCD'.DIFFd),NSSPE)
INCLUDE 'IOCOM.INC' CALL CXRZER('dispersal coe£.',DIFF,NSSPE,'GE'.ERR)
INCLUDE 'CNTCOM.INC' IF(ERR) GO TO 999
REAL'S EFF(NSSPE), DIFF(NSSPE), GENR(NSSPE).MASSL.LENGTH,FACTOR
INTEGER KSEO (NSNOD. NSSPE), LN (2) , LS (1), NEL.
CHARACTER TYPE*4,TYPEON«4
LOGICAL ERR
SAVE TYPEON
C VARIABLE DESCRIPTION-
TYPE*4 : ELEMENT TYPE: 'SIMP' OR 'CNDF'
TYPEON*4 : CURRENT ELEMENT TYPE "ON" OR ACTIVE
MASSL : MASS PER UNIT LENGTH - 'CNDF' ELEMENTS.
LENGTH : FLOW PASSAGE LENGTH - 'CNDF' ELEMENTS.
FACTOR : UPWIND FACTOR - 'CNDF' ELEMENTS.
ERR : ERROR FLAG
C CALL ZEROR (GENR, NSSPE, 1)
C CALL FREER('G',GENR(1), NSSPE)
FACTOR - l.ODO
CALL FREERCF'.FACTOR.l)
CALL CKRRNGC upwind factor' , FACTOR, 1. 0. 000. 'LELE' . 1 . ODO, ERR)
IF (ERR) GO TO 999
IF(NOOATA) FACTOR - -l.ODO
UPDATE SYSTEM BANDWIDTH
DO 40 N-l. NSSPE
LS(1) - N
CALL ELBAN (KSEO, NSNOD, NSSPE, LN.NENOD, LS.NESPE, MSBAN)
C—0.0 INITIALIZATION
C
NESPE - 1
NENOD - 2
C
C—1.0 GET ELEMENT TYPE
C
TYPE = 'SIMP1
CALL FREECCT'.TYPE,4,1)
WRITE ELEM. DATA TO ND1
WRITE(ND1) TYPE
WRITE(ND1) LN(1),LN(2),MASSL.LENGTH,FACTOR.
(DIFF(N),N=1,NSSPE)
UPDATE ELEMENT COUNT
NFELM = NEL
-UNDEFINED ELEMENTS
IF((TYPE.NE.'SIMP').AND.(TYPE.NE.'CNDF')) THEN
ERR • .TRUE.
CALL ALERT(
+ 'ERROR: Flow element type '//TYPE//' 15 not available',
+ 'S'.'S')
GO TO 999
ENDIF
C
C—2.0 REPORT TABLE HEADER IF NECESSARY
C
IF((NEL.EQ.NESTRT).OR.(TYPE.NE.TYPEON)) THEN
TYPEON = TYPE
SIMPLE ELEMENTS
IF(TYPEON.EQ.'SIMP') THEN
IF(ECHO) WRITE(NTW,2020)
WRITE(NOT.2020)
ELSEIF(TYPEON.EQ.'CNDF') THEN
IF(ECHO) WRITE(NTH,2022)
WRITE(NOT,2022)
ENDIF
ENDIF
CONV-DIFF ELEMENTS
2020 FORMAT(/,3X,'Num Type I
2022 FORMAT!/,3«,'Num Type I
+Fact Spec Dlsp.Coef.',/)
C—3.0 SIMPLE ELEMENTS
C
IF(TYPE.EQ.'SIMP') THEN
J Spec Fllt.Eff',/)
J M/Length Length
Upw.
c REPORT ELEMENT DATA
IF(FACTOR.NE.-l.ODO) THEN
WRITE(NOT, 2040) NEL,TYPE, LN(1),LN(2),MASSL,LENGTH,
» FACTOR,SID(l),DIFF(l)
ELSE
WRITE(NOT,2041) NEL,TYPE, LN(1),LN(2),MASSL,LENGTH,
+ ' default '.SID(1),DIFF(1)
ENDIF
IF (NSSPE.GE.2)
* WRITE(NOT,2042) (SID(N).DIFF(N).N-2.NSSPE)
IF(ECHO) THEN
IF(FACTOR.NE.-l.ODO) THEN
WRITE(NTW.2040) NEL.TYPE.LN(1).LN(2).MASSL,LENGTH.
+ FACTOR,SID(1),DIFF(1)
ELSE
WRITE (NTW,2041) NEL,TYPE,LN(1),LN(2),MASSL,LENGTH,
+ ' default '.SID(1),DIFF(1)
ENDIF
IF (NSSPE.GE.2)
* WRITE(NTW,2042) (SID(N),DIFF(N).N=2.NSSPE)
ENDIF
2040 FORMAT(2X,14.IX,A4,214.3(Gil.4),IX,A4.1(Oil.4))
2041 FORMAT(2X,14.IX,A4.214.2(Gil.4),A11,IX,A4,1(Gil.4))
2042 FORMAT((53X,A4,1(G11.4)))
ENDIF
RETURN
999 CALL ALERT(
^'WARNING: All flow element data has been deleted.'.'S','S')
NFELM - 0
CLOSE(ND1. STATUS-'DELETE')
RETURN
END
Append -6
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
SUBROUTINE KINELM
C— SUB: KINELM - COMMAND
C
C
C
C .
C
C
C
C .
C
C
C .
C
KINELEM
K=nl
n2,n3,...
M,n5....
K=nl* , /,
••«'»/»
<
n6 I=n7 K«n8
C . END' )
TO READ ft WRITE KIN ELEMENT DATA TO FILE *.KIN
Kinetics elements:',/.
nl ** rate coefficient matrix ID number',/,
n2,n3,...B lat row rate coef. matrix',/,
n4,n5,...s 2nd row rate coef. matrix'./.
additional rows as necessary',/.
end of rate coef. matrices subgroup',/.
n6 = elera, number; n7 = node number',/.
n8 = rate coefficient matrix ID number',/.
COMMON MTOT.N?,IAU)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOH.INC'
INTEGER NK
LOGICAL ERR,FOUND
EXTERNAL KINELO
C VARIABLE-
-DESCRIPTION
C ERR ERROR FLAG
C FOUND FILE FOUND FLAG
C KINELO SUB. TO READ t WRITE KIN ELEM DATA
C NK KATE COEF. MATRIX ID NUMBER
C NENOD NUMBER OF ELEMENT NODES
C NESPE NUMBER OF SPECIES PER ELEMENT (CURRENTLY=NSSPE)
C NESTRT FIRST ELEMENT NUMBER
C—0.0 INITIALIZATION
C
ERR = .FALSE.
NENOD - 1
IF(ECHO.OR.(MODE.EQ.'INTER')) HXITE(NTH,2000)
WRITE(NOT,2000)
2000 FORMAT!/,' ===== KINELEM: KINETICS ELEMENTS')
C
C~1.0 CHECK TO SEE IF SYSTEM NODES C EQUATION NUMBERS ARE DEFINED
C
CALL CXSYSU.ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C—2.0 OPEN .KIN
C
INQUIRE (FILE=FNAME(1:LFNAME)//'.KIN',EXIST=FOUND)
IF((.NOT.FOUND).OR.(NKINEL.EQ.0)) THEN
CALL NOPEN(ND1,(FNAME(1:LFNAME)//'.KIN'),'UNFORMATTED')
ELSEIF(FOUND) THEN
WRITE(NTH.2200)
HRITEINOT,")
WRITE(NOT,2200)
2200 FORMAT I
+ ' "" NOTE: Additional kin. elements being added to system.')
OPEN(ND1,FILE-(FNAME(1:LFNAME)//'.KIN').STATUS-'OLD',
+ FORM-'UNFORMATTED')
CALL APPEND(ND1)
ENDIF
C
C--3.0 GET RATE COEFFICIENT ARRAYS
C
CALL DELETE('TEMP')
CALL OEFINR('TEMP'.MPTEMP.NSSPE.l)
30 CALL FREE
IF (ECO) GO TO 40
NK - 0
CALL FREEICK'.NK, 1)
CALL CKIRNGCrate coef. matrix ID', NK. 1, 0, ' LTLE', 9. ERR)
IF(ERR) GO TO 999
CALL DELETE('KIK'//CHAR(NK+48))
CALL DEFINR!'KIK'//CHAR(NK+48), MPKIK(NK).NSSPE.NSSPE)
CALL GETKIK(IA(MPKIK(NK)),IA(MPTEMP),NK,ERR)
IF(ERR) GO TO 999
GO TO 30
C
C—4.0 GENERATE ELEMENT DATA; REPORT BANDWIDTH
C
40 WRITE(NTW,2400)
HRITEINOT.2400)
2400 FORMAT!/,' — Kinetics Elements'.//.6X.'Elem
IF(ECHO) WRITE(NTW,2410) MSBAN
HRITE(NOT.«)
HRITEINOT,2410) MSBAN
2410 FORMAT!' "" NOTE: Current system bandwidth Is:1, IS)
C
C--5.0 ORDERLY COMPLETION: CLOSE FILE ND1; SKIP TO "END"
C
CLOSE (HDD
IFIMOOE.EC.'INTER') RETURN
500 IF(EOC) RETURN
CALL FREE
GO TO 500
C
C— ABORT COMMAND
C
999 CALL ALERT<
(•'WARNING: All kinetics element data has been deleted.','$','$')
NKINEL=0
CLOSE (ND1, STATUS*>'DELETE')
DO 900 NK-1,9
900 CALL DELETE('KIK'//CHAR(NK+48) )
CALL DELETECTEMP')
CALL ABORT
RETURN
END
SUBROUTINE GETKIK(XIK,TEMP,NK, ERR)
C—SUB:GETKIK - READS AND REPORTS KINETICS RATE COEF. ARRAYS
INCLUDE 'CNTCOM.INC1
INCLUDE 'IOCOM.INC'
REAL'S KIK(NSSPE.NSSPE), TEMP(NSSPE)
INTEGER NK
LOGICAL ERR
IF (ECHO) HRITE(NTW,2000) NK
WRITE(NOT,2000) NK
2000 FORMAT!/,' = Kinetics Rate Coef. Matrix: KlnlD =',13)
C
C—1.0 READ IK)
C
DO 110 1-1,NSSPE
IF(MODE.EQ.'INTER') WRITE(NTW,2100) I
2100 FORMAT(' •• Enter terms In row number: '.14)
CALL FREE
IF(EOD) THEN
CALL ALERT!
+ 'ERROR: Data expected. Data subgroup terminator found.',
+ 'S'.'S'I
ERR •' .TRUE.
RETURN
ENDIF
CALL FKEERC ', TEMP,NSSPE)
DO 100 J=l. NSSPE
100 KIK(I.J) = TEMP(J)
110 CONTINUE
C
C—2.0 REPORT FIVE COLUMNS AT A TIME
C
DO 200 J1<=1,NSSPE,5
J2 = MIN(NSSPE,Jl+4)
IF(ECHO) WRITE(NTH, 2200) (SID(J),J-J1,J2)
«RITE(NOT,2200) (SID(J) , J=J1,J2)
2200 FORMAT!/, 12X, 5 (:3X,A4.6X»
DO 200 1=1,NSSPE
IF(ECHO) WRITE(NTW,2210) SID (I) , (KIK (I. J) . J-J1. J2)
HRITEINOT, 2210) SID(I), (KIK (I, J) , J-J1, J2)
200 CONTINUE
2210 FORMAT(6X,A4,2X,S(:G11.3,2X)|
RETURN
END
SUBROUTINE KINELO(NEL,LN,ERR)
C—SUB:KINELO - READS ADDITIONAL KINETICS ELEMENT DATA,
C WRITE KINETICS ELEMENT DATA TO FILE ND1,
C UPDATES SYSTEM BANDWIDTH, AND REPORTS ELEMENT DATA
COMMON MTOT,NP,IA(1)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC1
LOGICAL ERR
INTEGER LN(1), NK
NENOD - 1
NESPE = NSSPE
Node KlnlD')
C
C 'TEMP' STORES SPECIES CONNECTIVITY ARRAY, LS (NSSPE), USED BY ELBAN C—1.0 GET RATE COEFFICIENT MATRIX ID
CALL DELETE('TEMP')
CALL DEFINICTEMP',MPTEMP.NSSPE,1)
DO 42 N-l, NSSPE
42 IAIHPTEMP+N-1) = N
NESTRT • NKINEL+1
CALL ELGEN(KINELO,NENOD.NESTRT.NSNOD,ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
CALL DELETE I'TEMP'I
CALL FREEI CK'.NK, 1)
IFIMPKIXINK).EQ.0) THEN
.CALL ALERTI'ERROR: Rate coefficient matrix not defined',
* 'S'.'S')
ERR = .TRUE.
RETURN
ENDIF
-2.0 WRITE DATA TO ND1
HRITEIND1) LN(l), NX
-3.0 UPDATE SYSTEM BANDWIDTH (NOTE: IA(MPTEMP)=LS(NSSPE) )
Append -7
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
CALL ELBAN(IA(MPKSEQ),NSNOD.NSSPE,LN.NENOD,IA(HPTEMP).NESPE.MSBAN)
C
C—4.0 REPORT DATA
C
IF(ECHO) WRITE(NTH,2000) NEL, LN(1), NK
KRITE (NOT. 2000) NEL. LN(1). NX
CALL ABORT
RETURN
ENDIF
CLOSE(N01)
IF(NKINEL.GT.O) CLOSE(ND2)
2000 FORMAT(6X,I4,4X,I4,4X,I4)
NKINEL = NEL
RETURN
END
SUBROUTINE FORMFO
C — SUB:FORHFO - COMMAND TO FORM CONTAM. DISPERSAL MASS TRANS ARRAY [N]
C
C .' FORM- [Ml F«ceee Form (F| , cccc - FULL or BAND',/,
C .' n2,n3,M W«nS n2,n3,n4 = elem: first, last, Incr.1,/,
C .' ... n5 » element flow rate',/.
C .' END')
C
IMPLICIT REALM (A-H.O-Z)
COMMON MTOT. NP.IA(l)
INCLUDE 'IOCOM.INC1
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
C
C— 5.0 DELETE ARRAYS
C
CALL DELETE (' TEMP '
CALL DELETE CCONT'
CALL DELETE ('VCD '
CALL DELETE ('EFF '
CALL DELETE ('DIFF'
CALL DELETE COENR'
CALL DELETE ('HE '
RETURN
END
SUBROUTINE STEADY
C SOLUTION (C) IS WRITTEN OVER (E)
C
C
C . STEADY Steady state solution.'./.
C . ... n4 H element flow rate',/,
C . ... • n8,n9, . . . *> spec. cone, or gen. rate* , /,
C . END',//,
ERR - .FALSE.
C
C—0.0 WRITE HEADER AND READ ARRAY FORM
C
WRITE(NOT. 2000)
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTW,2000)
2000 FORMAT)/,
t, a== FORM-IW]: FORM MASS TRANSPORT MATRIX [W] ')
FORM = 'FULL'
CALL FREECCF'.FORM,4,1)
IF(FORM.EQ.'FULL') THEN
IF (ECHO) WRITE (NT*. 2002)
WRITE(NOT,*)
WRITE(NOT,2002)
2002 FORMATC *• NOTE: [W] being formed In FULL form.')
ELSEIF(FORM.EQ.'BAND') THEN
IF (ECHO) WRITE(NTW,2004)
WRITE(NOT,«)
WRITE(NOT,2004)
2004 FORMAT!' ** NOTE: (W| being formed In BAND form.')
ELSE
CALL ALERT('ERROR: '//FORM//' not defined.','$'.•S')
CALL ABORT
RETURN
ENDIF
C
C--1.0 CHECK IF FLOW SYSTEM AND ELEMENT DATA ARE DEFINED
C
CALL CKSYS(2,ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C—2.0 DEFINE AND INITIALIZE ARRAYS
CALL DELETE('TEMP'
CALL DELETE CCONT'
CALL DELETE!'VCD '
CALL DELETE!'EFF '
CALL DELETE('DIFF1
CALL DELETE I'GENR'
CALL DELETE CHE '
CALL DELETE('W '
IF (FORM.EQ.'FULL') CALL DEFINRCW
IF (FORM.EQ.'BAND') CALL DEFINRCW
CALL DEFINRCWE ' , MPWE.NFELM, 1)
CALL DEFINR COENR', MPGENR. NSSPE. 1)
CALL DEFINR CDIFF'.MPDIFF, NSSPE, 1)
CALL DEFINRCEFF ',MPEFF.NSSPE,1)
CALL DEFINR('VCD ',MPVCD,NSNCO,1)
CALL DEFINRCCONT',MPCONT,NSNOD,1)
IMPLICIT REAL*8(A-H, 0-Z)
COMMON MTOT.NP.IAd)
INCLUDE 'IOCOH.INC'
INCLUDE 'CMDCOM.INC'
INCLUDE 'CNTCOM.INC'
COMMON /STDCOM/ MPEDAT
LOGICAL ERR
ERR - .FALSE.
WRITE(NOT,2000)
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTW, 2000)
2000 FORMAT!/,' -= STEADY: STEADY STATE SOLUTION')
C
C—1.0 CHECK IF FLOW SYSTEM AND ELEMENT DATA ARE DEFINED I AVAIL
C
CALL CKSYSI2.ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C—2.0 DEFINE AND INITIALIZEARRAYS
C
',MPF,NSEQ,NSEQ>
1,MPF.NSEQ. 2'MSBAN-l)
CALL DELETE CEDAT'
CALL DELETE CCONT •
CALL DELETE!' VCD '
CALL DELETE! 'EFF '
CALL DELETE! 'DIFF'
CALL DELETE!' GENR'
CALL DELETE CE
CALL DELETE CHE '
CALL DELETE CF '
CALL DEFINR ('F '
CALL DEFINRCWE '
CALL DEFINRCE '
CALL DEFINR ('GENR'
CALL OEFINRCDIFF'
CALL DEFINR I ' EFF '
CALL DEFINR C VCD '
CALL DEFINR CCONT'
CALL DEFINR CEDAT'
)
)
)
)
)
)
)
)
>
,MPF,NSEQ,2*MSBAN-1)
.KPWE.NFELM, 1)
, MPE.NSEQ, 1)
.MPGENR, NSSPE, 1)
.MPOIFF, NSSPE, 1)
, HP EFF, NSSPE. 1)
.HPVCD, NSNOD. 1)
, HP CONT, NSNOD, 1)
.MPEDAT, NSSPE, 1)
C
C--3,
C
C
C—4.
C
0 GET ELEMENT FLOW RATES (WE)
CALL ZEROR(IA(MPWE),NFELM.1)
CALL READWE(ERR)
IF(ERR) THEN.
CALL ABORT
RETURN
ENOIF
0 FORM (W|
OPEN (FNAME (1: LFNAME) //'.KIN') ,STATUS='OLD',
+FORH-'UNFORMATTED')
ENDIF
CALL FORMFIIA(MPKSEQ),IA(MPF).IA(MPHE),FORM,ERR)
IF(ERR) THEN
C
C—3.0 GET ELEMENT FLOW RATES (WE)
C
CALL ZEROR(IA(MPWE),NFELM,1)
CALL README(ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C--4.0 FORM (F)
C
OPEN(N01,FILE=(FNAME(1:LFNAME)//'.FEL'),STATUS='OLD'
+FORM-'UNFORMATTED')
IF (NKINEL.GT.O) THEN
OPEN(ND2,FILE"(FNAME(1:LFNAME)//1.KIN'),STATUS-'OLD'
+FORM-'UNFORMATTED'I
ENDIF
CALL FORMF(IA(MPKSEO),IA(MPF),IA(MPWE),
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
CLOSE(ND1)
IFINXINEL.GT.O) CLOSE (ND2)
'BAND',ERR)
Append -8
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
c
C—6.
C
c
C--7.
C
c
C--5.0 FORM (E)
C
CALL ZEROR(IA(HPE),NSEQ, 1)
CALL FORMEX(ERR)
IF(EKR) THEN
CALL ABORT
RETURN
ENDIF
0 MODIFY (E) AND [F] FOR PRESCRIBED CONCENTRATIONS
CALL MODIF(IA(MPKSEQ),IA(MPF),IA(MPE))
0 SOLVE
CALL FACTCA(LMM?F),NSEQ,MSBAN,ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
CALL SOLVCA (IA (MPF) . IA (MPE) . NSEQ. MSBAN, ERR)
IF (ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C—8.0 REPORT SOLUTION
C
IF(ECHO) WRITE(NTH.2800)
WRITE(NOT.2800)
2800 FORMAT!/,' == Response: Node Concentrations')
CALL RPRTEN(IA(MPE),IA(MPKSEO))
C
C—9.0 DELETE ARRAYS
C
CALL DELETE CEDAT1)
CALL DELETE('CONT')
CALL DELETE('VCD ')
CALL DELETE('EFF ')
CALL DELETE('DIFF')
CALL DELETE CGENR')
CALL DELETE CE ')
CALL DELETE ('HE ')
CALL DELETECF ')
RETURN
END
c FORMEX
SUBROUTINE FORMEX(ERR)
C—SUB:FORMEX - READS I REPORTS NODAL CONTAMINANT EXCITATION DATA
COMMON MTOT.NP.IAU)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
EXTERNAL EXDATO
WRITE(NOT,2100)
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTH.2100)
2100 FORMAT!/,
+' «• Excitation: Contaminant Concentration or Generation')
CALL DATGEN(EXDATO,NSNOD,ERR)
CALL RPRTEN (IA (MPE) , IA (MPKSEQ) )
RETURN
END
SUBROUTINE EXDATO(N.ERR)
C—SUB:EXDATO - CALLS EXDAT1 PASSING ARRAYS
C
COMMON MTOT.NP.IAd)
INCLUDE 'CNTCOM.INC'
COMMON' /STDCOM/ MPEDAT
LOGICAL ERR
CALL EXDAT1(IA (MPE).IA(MPKSEQ),IA(MPEDAT) ,N. ERR)
RETURN
END
ELSE
WRITE(NTH,2000) N
NRITE(NOT,«)
WRITE(NOT,2000) N
000 FORMAT(' "*•• ERROR: Node ',15,'
ERR - .TRUE.
RETURN
ENDIF
10 CONTINUE
RETURN
END
Is not a defined flow node.')
SUBROUTINE MODIF(XSEQ,F. E)
C—SUB:MODIF - MODIFIES [F| AND (E) FOR C-PRESCRIBED DOFS
INCLUDE 'CNTCOM.INC'
REAL'S F(NSEQ,2*MSBAN-1),E(NSEQ)
INTEGER KSEQ(NSNOD,NSSPE)
DO 10 N-l,NSNOD
DO 10 M=l, NSSPE
NEQ - KSEOIN.M)
NNEQ - ABS(NEQ)
IF(NEO.LT.O) THEN
FINKEQ.MSBAN) - F(NNEQ,MSBAN)•!.OD15
E(NKEO) » E(NNEQ)*F(NNEQ,MSBAN)
ENDIF
10 CONTINUE
RETURN
END
SUBROUTINE TIMCON
SUB:TIMCON - COMMAND TO FORM CCNTAM. DISPERSAL EIGENVALUE PROBLEM
C
C
C
C
c
c
c
c
c
c
C—HELP LIST—•
C
C
c
c
c
c-
c
c
c
[[V|-1[F| - (1/THIHE) = (0)
IIHERE: [V| - FLOW VOLUMETRIC MASS MATRIX (DIAGONAL)
[F| • FLOW SYSTEM FLOW MATRIX
(E) - (RIGHT) EIGENVECTORS
T = CONTAM. DISPERSAL TIME CONSTANTS
TO EVALUATE TIME CONSTANTS. EIGENVECTORS ARE NOT FOUND.
TIMECOHS E=nl
n2,n3,n4 W=n5
END')
Time constant solution, nl ° epsllon',/,
n2.n3.n4 » elem: first, last, Incr.'./,
nS - element flow rate',/.
MPTC TC (NSEQ) TEMPORARY ARRAY FOR STORAGE OF TIME CONS
IMPLICIT REAL'S(A-H.O-Z)
COMMON MTOT.NP.IAd)
INCLUDE 'IOCOM.INC1
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
ERR = .FALSE.
C
C—0.0 WRITE HEADER AND READ PRECISION
C
NRITE(NOT,2000)
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTH,2000)
2000 FORMAT!/.
+' ran riMECONS: TIME CONSTANTS - CONTAMINANT DISPERSAL SYSTEM ')
EP1 = El'
CALL FREERCE'.EPl.l)
HRITE(NOT.2010) EP1
IF (ECHO) HRITE (NTH, 2010) EP1
2010 FORMAT(/' Ba Convergence parameter, epsllon: ', G10.3)
C
c—i.o CHECK :;F FLOW SYSTEM AND ELEMENT DATA ARE DEFINED
c
CALL CK!iYS(2.ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
SUBROUTINE EXDAT1(E,KSEQ.EDAT.N, ERR)
C—SUB:GDAT1 - READS CONTAMINANT EXCITATION DATA
C
COMMON MTOT.NP.IAd)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
REAL'S E(NSEQ). EDAT (NSSPE)
INTEGER KSEQ(NSNOD,NSSPE)
LOGICAL ERR
CALL ZEROR(EDAT.NSSPE,1)
CALL FREER I'G',EDAT, NSSPE)
DO 10 M=l,NSSPE
NEQ = ABS (KSEQ (N, M))
IFfNEQ.NE.O) THEN
E(NEQ) - EDAT(H)
C—2.0 DEFINE ARRAYS
C
CALL DELETE('TEMP
CALL DELETE ('CONT
CALL DELETE)'VCD
CALL DELETE I'EFF
CALL DELETE CDIFF
CALL DELETE CGENR
CALL DELETE('TC
CALL DELETE CVM
CALL DELETE CHE
CALL DELETE)'F
CALL DEFINRCF
CALL DEFINRCWE
CALL DEFINRCVM
CALL DEFINRCTC
CALL DEFINR CGENR
CALL DEFINR CDIFF
CALL DEFINRCEFF
CALL DEFINRCVCD
. MPF, NSEQ.NSEQ)
,MPWE.NFELM, 1)
.MPVM.NSEQ. 1)
,MPTC.NSEQ, 1)
,MPGENR,NSSPE. 1)
,MPDIFF.NSSPE,1)
,HPEFF.NSSPE,1)
,HPVCD.NSNOD,1)
Append -9
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
CALL DEFINRCCONT'.MPCONT.NSNOO.l)
C
C—3.0 GET ELEMENT FLOW RATES {ME)
C
CALL ZERORUA(MPNE).NFELM,1)
CALL README(ERR)
IF(ERR) THEM
CALL ABORT
RETURN
ENOIF
C
C—4.0 FORM [F|
C
OPEN(NDl,FILE-(FNAME(l:LFNAME)//'.FEL'),STATUSa'OLD',
+FORM-'UNFORMATTEO')
IF(NKINEL.GT.O) THEN
OPEN(ND2,FILE-(FNAHE(1:LFNAHE)//'.KIN'),STATUS-'OLD',
+FORM-'UN FORMATTED')
ENDIF
CALL FORMF (IA (HPKSEQ), IA (MPF), IA (HFWE) , ' FULL'. ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENOIF
CLOSE (HDD
IF(NKINEL.GT.O) CLOSE(ND2)
C
C—5.0 FORM VOLUMETRIC MASS MATRIX
C
CALL ZEROR(IA(MPVM),NSEQ.l)
CALL FORHVM (IA (MPKSEQ), IA (M?V). IA (HPVCD). IA (MPVM))
C
C—6.0 COMPUTE i REPORT NOMINAL TIME CONSTANTS
C
IF(ECHO) WRITE(NTH,2600)
MRITE(NOT.2600)
2600 FORMAT!/.1 — Nominal Tims Constanta')
CALL ZEROR(IA(MPTC),NSEQ,1)
CALL NOMTC (IA (MPKSEO), IA (MPVM) , IA (MPF) . IA (MPTC))
CALL RPRTEN(IA(MPTC).IA.HDT
OPTIONALLY EQUAL STEP TIME HISTORIES MAY BE GENERATED
FLONDAT [T-nl,n2,n3l Generate element Clow time histories.',/,
TIME=nl nl = time'./.
nl,n2,n3 N=n4 nl,n2,n3 » node: first, last, Incr.',/,
... n4 « element mass flow rate.',/.
:'./,
END',//,
IMPLICIT REAL*8(A-H, 0-Z)
CAL-SAP: DATA I COMMON STORAGE
COMMON MTOT.NP.IAd)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
c
C— FLODAT: DATA t COMMON STORAGE
C
C— — — — D I C T
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
I 0 N A 1
POINTER VARIABLE
TIME (3)
MPNE ME (NFELM)
TIME
DAT(l)
DAT (2)
HIST
DESCRIPTION
! START TIME, ENDTIME, TIMESTEP
: CURRENT ELEMENT MASS FLOH VALUES
O R ¥ DATA
* - - - Time histories of excitation data
I defined as step-wise functions of
are
time
*- using arbitrary values or, optionally.
1
*_
1
- - *-
TM(2)
MPTDAT TDATI2)
MPHDT1 HDT1 (NFELM)
MPNDT2 NDT2 (NFELM)
generated Intermediate values of
equal step size.
1 _
TM(1)
: CURRENT ARBITRARY TIME VALUES
: ELEM. FLOH DATA AT TDAT(l)
: ELEM. FLOH DATA AT TDAT(l)
Append -10
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
COMMON /FLODT/ MPTDAT.MPHDT1.MPHDT2
REAL'S TIME(3)
LOGICAL ERR
ERR « .FALSE.
WRITE(NOT.2000)
IF(ECHO.OR.(MODE.Ed.'INTER')) WRITE(NTH,2000)
2000 FORMAT!/,' — FLOMDAT: ELEMENT FLOW TIME HISTORY DATA')
C
C—1.0 CHECK TO SEE IF PERTINENT DATA HAS BEEN DEFINED
C
CALL CKSYS(2,ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C—2.0 GET DATA GENERATION CONTROL DATA
C
TIME(l) - O.ODO
TIME(2) - O.ODO
TIME(3) - O.ODO
CALL FREER CT', TIME (1) , 3)
CALL CKIZERCtlme step', TIME (3) . 1,'GE', ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ELSEIF(TIME(3).GT.O.ODO) THEN
IF (TIME (2) .LT.TIME(D) THEN
CALL ALERT(
+ 'ERROR: Final time must be greater than Initial time.',
+ '$','$')
CALL ABORT
RETURN
ENDIF
IF (ECHO) WRITE (NTH, 2220)
WRITE(NOT,2220)
2220 FORMAT!/,' «= Generation Control Variables')
IF (ECHO) WRITE (NTH, 2230) (TIME (I) , 1=1, 3)
WRITE(NOT,2230) (TIME(I),1=1,3)
2230 FORMAT(/,
.' Initial time '.G10.3,/,
.' Final time '.G10.3,/,
.' Time step increment *,G10.3)
ENDIF
C
C—3.0 OPEN .WDT
C
CALL NOPENIND1, (FNAMEU :LFNAME) //' .WDT'), 'UNFORMATTED')
C
C—4.0 READ < GENERATE FLOW DATA
C
WRITE(NOT,2400)
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTH,2400)
2400 FORMAT!/,' — Element Mass Flow Time History Data')
C
C 4.1 DEFINE t INITIALIZE ARRAYS
C
CALL DELETE CTDAT')
CALL DELETE CWDT1')
CALL DEFINR('HDT1',MPHDT1,NFELM,1)
CALL DEFINR('TDAT',MPTDAT,1,2)
CALL ZEROR(IA(MPHDT1),NFELM,1)
CALL ZEROR(IA(MPTDAT),1,2)
IF(TIME(3).GT.O.ODO) THEN
CALL DELETE('HDT2')
CALL DELETE CHE •)
CALL DEFINRCHE ' .MPHE.NFELM, 1)
CALL DEFINR('WDT2',MPWDT2,NFELM,1)
CALL ZEROR(IA(KPNDT2),NFELM, 1)
CALL ZEROR(lA(MPNE),NFELM,1)
ENDIF
C
C—4.2 GENERATE VALUES ( WRITE TO .HDT
C
IF(TIME(3).GT.O.ODO) THEN
CALL GENWD1(IA(MPHE),IA(MPTDAT),IA(MPHDT1),IA(MPHDT2),TIME,ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
ELSE
CALL GENHD2 (IA (MPTDAT). IA (MPWDT1) , ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
ENDIF
C--5.0 DELETE ARRAYS, CLOSE ELEMENT FLOB DATA FILE, SKIP TO -END"
C
CALL DELETE CHE ')
CALL DELETE CHDT2')
CALL DELETE CHDT1'I
CALL DELETE CTDAT')
CLOSE (ND1)
IF(MODE.EQ.'BATCH') THEN
500 IF(EOC) RETURN
CALL FREE
GO TO SOO
ENDIF
RETURN
END
c GENWD1
SUBROUTINE GENWD1(HE,TDAT.HDTl,WDT2.TIME, ERR)
C—SUB: GENWD1 - GENERATES ELEMENT MASS FLOW DATA, AT EQUAL TIME STEP
C INTERVALS, FROM GIVEN ARBITRARY DISCRETE TIME DATA
IMPLICIT REAL*8(A-H,O-Z)
INCLUDE 'ICCOM.INC'
INCLUDE 'CNTCOM.INC'
C
C FLOHDAT: DATA t COMMON STORAGE
C
COMMON /FLODT/ MPTDAT,HPWDT1,HPKDT2
LOGICAL ERR
-— GENHD1: DATA I COMMON STORAGE
REAL'8 NF. (NFELM), TDAT (2), WDT1 (NFELM) ,HDT2 (NFELM) , TIME (3)
•1.0 GET FIRST TWO TIME HISTORY RECORDS
CALL GEl'TDT(TDAT)
IF(EOC) THEN
CALL ALERT('ERROR: Insufficient data.','$','S')
ERR - .TRUE.
RETURN
ENDIF
CALL GETWDT(WDT1,HDT2.ERR)
IF(ERR) RETURN
CALL GETTDT(TDAT)
IF(EOC) THEN
CALL ALERT('ERROR: Insufficient data.'. 'S', 'S')
ERR - .TRUE.
RETURN
ELSEIF(TDATIl).LT.TDAT(21) THEN
CALL ALERT('ERROR: Time data out of sequence.','$','$')
ERR - .TRUE.
RETURN
ENDIF
CALL GETWDT(WDT1,WDT2.ERR)
IF(ERR) RETURN
•2.0 GENERATION TIME LOOP
DO 200 T-TIME(1),TIME(2).TIME<3)
—2.1 UPDATE FLOW TIME HISTORY DATA IF NEEDED
20 IF(T.GT.TDAT(D) THEN
CALL GETTDT(TDAT)
IF(EOC) THEN
CALL ALERT('ERROR: Insufficient data.','$','$')
ERR = .TRUE.
RETURN
ELSEIF(TDATd) .LT.TDAT(2)| THEN
CALL ALERTCERROR: Time data out of sequence. '.'S'.'S'l
ERR = .TRUE.
RETURN
ENDIF
CALL GETWDT(WDT1,WDT2.ERR(
IF(ERR) RETURN
GO TO 20
ENDIF
-2.2 COMPUTE INTERPOLATION FRACTION
XT = (T--TDAT(2))/(TDAT(1)-TDAT(2)(
—2.3 COMPUTE (WE(T) )
DO 23 N'-l,NFELM
HE(N) - HDT2(N) -f XT* (HDT1 (N)-HDT2 (N) )
23 CONTINUE
C
C 2.4 WRITE TIME, I WE (TM TO ND1
C
WRITE (HDD T
WRITE(ND1) (WE(I),1=1,NFELM)
200 CONTINUE
—3.0 WRITE ONE ADDITIONAL TIME VALUE TO DISK
WRITE(ND1) T
' RETURN
END
SUBROUTINE OETWDT(WDT1,HDT2,ERR)
—SUB: GETHDT - UPDATES ELEMENT FLOW DATA VALUES
INCLUDI! 'CNTCOM.INC'
LOGICAL ERR
REAL'S WDT1(NFELM), HDT2(NFELM)
EXTERNAL WDATO
C
C—1.0 UPDAT1! 'OLD' DATA VALUES; INITIALIZE 'NEW' DATA VALUES
C
DO 10 N=l.NFELM
WDT2 (N) - HDT1 (N)
WDT1(N) a 0.ODO
10 CONTINUE
C
C—2.0 READ NEW VALUES
C
CALL DATGENIHDATO,NFELM,ERR)
IF(ERR) RETURN
Append -1 1
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
CALL RPRTNO(HDTKl), NFELM, 'Elem')
RETURN
END
SUBROUTINE NOATO (N , ERR)
C — SUB:NDATO - CALLS WDAT11 PASSING ARRAYS
C
COMMON MTOT.NP.IAU)
INCLUDE 'CNTCOM.INC'
COMMON /FLODT/ MPTDAT, MPNDT1 , MPHDT2
LOGICAL ERR
CALL WDAT1 (IA (MPWDT1 ) , NFELM, N)
RETURN
END
SUBROUTINE WDAT1 (WDT1, NFELM, N)
C— SUB:»DAT1 - READS ELEMENT MASS FLOW RATE TIME HISTORY DATA
C
REAL'S NDTI (NFELM)
CALL FREER('W',WDTl(N).l)
RETURN
END
SUBROUTINE GENWD2 (TDAT, WDT1, ERR)
C— SUB: GENWD2 - GENERATES ELEMENT MASS FLOW DATA. AT GIVEN TIME STEP
C INTERVALS, FROM GIVEN DISCRETE TIME DATA
IMPLICIT REAL'S (A-H.O-Z)
INCLUDE 'ICCOM.INC'
INCLUDE 'CNTCOM.INC'
I 0 0 0 0 1
.' :'./.
.' END'.//.
IMPLICIT REAL'S (A-H.O-Z)
COMMON MTOT.NP.IA(l)
INCLUDE 'ICCOM.INC'
INCLUDE 'CNTCOM.INC'
POINTER VARIABLE DESCRIPTION
TIME (3) : START TIME, ENDTIME, TIMESTEP
MPE E (NSEQ, NSSPE) : CURRENT EXCITATION VALUES
TIME HISTORY DATA
1 defined as atop-wise functions of time
*- using arbitrary values or. optionally,
t generated Intermediate values of
*- equal step size.
1
DAT (2) - - •-
TM(2) TM(1)
MPTDAT TDAT(2) : CURRENT ARBITRARY TIME VALUES
MPEDT1 EDT1 (NSSPE, NSNOD) : EXCITATION DATA AT TDAT(l)
MPEDT2 EDT2 (NSSPE. NSNOD) : EXCITATION DATA AT TDATI2)
COMMON /EXCDT/ MPTDAT, MPEOT1, MPEDT2
REAL'S TIME (3)
LOGICAL ERR
ERR - .FALSE.
WRITE (NOT. 2000)
C FLOHDAT: DATA t COMMON STORAGE
C
COMMON /FLODT/ MPTDAT,HPWDT1.MPWDT2
LOGICAL ERR
EXTERNAL HDATO
C •
C GENHD2: DATA I COMMON STORAGE
C
REAL'S TDAT(2),MDT1(NFELM)
C—1.0 GET FIRST TIME HISTORY RECORD ( TDAT(l), WDT1 (NFELM) )
C
CALL GETTDT(TDAT)
IF(EOC) RETURN
TDAT(2) - TDAT(l)
CALL ZEROR(WDT1,NFELM,1)
CALL DATGENIMDATO,NFELM,ERR)
IF(ERR) RETURN
CALL RPRTNO(WDT1(1),NFELM,'Elem')
WRITE(ND1) TDAT(l)
WRITE(NDl) (WDTKIl.I'l.NFELM)
C
C—2.0 GET ADDITIONAL TIME HISTORY RECORDS
C
20 CALL GETTDT(TDAT)
IF(EOC) GO TO 300
IFITDAT(l).LT.TDAT(2)) THEN
CALL ALERTCERROR: Time data out of sequence.','$','S')
ERR - .TRUE.
RETURN
ENDIF
TDAT(2) - TDAT(l)
CALL ZEROR(WDT1,NFELM. 1)
CALL DATGEN(WDATO.NFELM,ERR)
IF(ERR) RETURN
CALL RPRTNO(HDT1(1),NFELM,'Elem')
WRITE(ND1) TDAT(l)
WRITE (NOD (WDTKI) ,1=1. NFELM)
GO TO 20
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTW.2000)
2000 FORMAT!/, ' — EXCITDAT: EXCITATION TIME HISTORY DATA')
C
C—1,
c
c
C—2.
C
0 CHECK TO SEE IF FLOW SYSTEM HAS BEEN DEFINED
CALL CKSYS(l.ERR)
IF (ERR) THEN
CALL ABORT
RETURN
ENOIF
0 GET DATA GENERATION CONTROL DATA
TIME(l) - O.ODO
TIME(2) - O.ODO
TIME(3) = O.ODO
CALL FREERCT'.TIME(l) ,3)
IF(TIME(3).LT.O.ODO) THEN
CALL ALERTCERROR: Tim step may not be negative.',1?1
CALL ABORT
RETURN
ELSEIF(TIME(3).GT.O.ODO) THEN
IF(TIME(2).LT.TIME(1)) THEN
CALL ALERT!
+ 'ERROR: Final time must be greater than initial time.
+ 'S'.'S'I
CALL ABORT
RETURN
ENDIF
IF(ECHO) WRITE(NTW,2220)
WRITE(NOT,2220)
I FORMAT!/,' = Generation Control Variables')
IF(ECHO) WRITE(NTW,2230) (TIME(I),1-1,3)
WRITE(NOT.2230) (TIME(I),1=1,3)
I FORMAT)/,
.' Initial time '.G10.3./.
.' Final tlrao '.G10.3,/,
.' Time step increment '.G10.3)
ENDIF
C--3.0 WRITE ONE ADDITIONAL TIKE VALUE TO DISK
C
300 WRITE(ND1) TDAT(l)
RETURN
END
SUBROUTINE EXCDAT
C—SUB:EXCDAT - COMMAND TO READ EXCITATION DATA I GENERATE STEPWISE
C TIME HISTORIES OF EXCITATION VALUES, EINSEQ).AND
C WRITES TIME HISTORIES IN FORMAT;
C
C • TIME
C .EDT
C—HELP LIST
C
C .' EXCITDAT [Tanl.n2.n3) Generate excitation time histories.',/.
C .' TIME-nl nl • time'./.
C .' nl.n2,n3 CG*n4,n5.... nl,n2,n3 = node: first, last, incr.',/.
C—3.0 OPEN .EDT
C
CALL NOPEN(ND1.(FNAME(1:LFNAME)//'.EOT'),'UNFORMATTED')
C
C—4.0 READ < GENERATE EXCITATION DATA
C
WRITE(NOT,2400)
IF(ECHO.OR.(MODE.EQ.'INTER')) WRITE(NTH,2400)
2400 FORMAT!/,• = Nodal Excitation Tim- History Data')
c
C 4.1 DEFINE ( INITIALIZE ARRAYS
C
CALL DELETE!TDAT')
CALL DELETE('EDT1')
CALL DELETE!'E ')
CALL DEFINRCE ' , HPE.NSEQ, 1)
CALL DEFINR('EDT1',HPEDT1,NSSPE,NSNOD)
CALL DEFINRCTDAT',MPTDAT.1,2)
CALL ZERORIIA(MPE).NSEQ.1)
CALL ZEROR(lA(MPEDTl).NSSPE.NSNOD)
CALL ZEROR(IA(MPTDAT).1.2)
IF(TIME(3).GT.O.ODO) THEN
CALL DELETECEDT2')
CALL DEFINR('EDT2'.MPEDT2,NSSPE.NSNOD)
Append -12
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
CALL ZEROR(IA(MPEDT2),NSSPE,NSNOD)
ENDIF
C
C 4.2 GENERATE VALUES t WRITE TO .EDT
C
IF(TIME(3) .GT.O.ODO) THEN
CALL GENED1(IA(MPKSEO),IA(MCE),IA(MPTDAT),IA(HPEDT1).IA(HCEDT2),
* TIME.ERK)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
ELSE
CALL GENED2 (IA (MPKSEQ) . IA (MCE) . IA (MPTDAT) . IA (MPEDT1) , ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
ENDIF
C
C—5.0 DELETE ARRAYS, CLOSE ELEMENT FLOW DATA FILE, SKIP TO "END"
C
CALL DELETECEDT2')
CALL DELETE('TDAT')
CALL DELETE CEDT1')
CALL DELETE CE ')
CLOSE(ND1)
II- (MODE. EQ. 'BATCH') THEN
500 IF(EOC) RETURN
CALL FREE
GO TO SOO
ENDIF
NEQ - ABS(KSEO(N,M)>
IF(NEQ.NE.O) E(NEQ) - EDT2(M,N) + XT"(EDT1(M.N)-EDT2(M,N))
23 CONTINUE
C
C 2.4 WRITE TIME, (E(TH TO ND1
C
HRITE(NDl) T
WRITE(ND1)
200 CONTINUE
-3.0 WRITE ONE ADDITIONAL TIME VALUE TO DISK
WRITE(ND1) T
RETURN
END
SUBROUTINE GETEDT(KSEO,E.EDT1.EDT2,ERR)
-SUB: GETEDT - UPDATES EXCITATION DATA VALUES
COMMON MTOT.NP.IAd)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
C
C— GETEDT: DATA I COMMON STORAGE
C
LOGICAL ERR
REAL'S EINSEOI,EDT1(NSSPE.NSNOD), EDT2(NSSPE.NSNOD)
INTEGER KSEQ(NSNOD,NSSCE)
EXTERNAL EDATO
RETURN
END
SUBROUTINE GENED1 (KSEO. E, TDAT, EDT1 . EDT2 , TIME , ERR)
C
C— 1.
c
C — SUB: GENED1 - GENERATES EXCITATION DATA, AT EQUAL TIKE STEP
C
C
C
c
C
c
C--1.
c
INTERVALS, FROM GIVEN ARBITRARY TIME DATA
IMPLICIT REAL'S (A-H.O-Z)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
GENEDI: DATA < COMMON STORAGE
REAL'S E (NSEQ) , TDAT (2) , EDT1 (NSSPE.NSNOD) , EDT2 (NSSPE. NSNOD) , TIKE (3)
INTEGER KSEQ (NSNOD, NSSPE)
0 GET FIRST TWO TIME HISTORY RECORDS
CALL GETTDT(TDAT)
IF(EOC) THEN
CALL ALERTCERROR: Insufficient data. ' , ' S ' , ' S ' )
ERR a .TRUE.
RETURN
ENDIF
CALL GETEDT (KSEQ, E.EDT1.EDT2, ERR)
IF (ERR) RETURN
CALL GETTDT(TDAT)
10
C--2.
C
c
C—3.
C
30
0 UPDATE 'OLD' DATA VALUES; INITIALIZE 'NEW' DATA VALUES
DO 10 N-l.NSNOO
EDT2(H,N) - EDTl(M.N)
EDTKM.N) • O.ODO
CONTINUE
0 READ NEW VALUES
CALL DATSEN (EDATO, NSNOD, ERR)
IF (ERR) RETURN
0 REPORT. VALUES
DO 30 N-l, NSNOD
DO 30 K=l. NSSPE
NEQ = ABS(KSEQ(N,M)I
IF(NEQ.NE.O) E(NEQ) = EDTKM.N)
CONTINUE
CALL RPRTENIEd) . IA (MPKSEQ))
RETURN
END
SUBROUTINE EDATO (N, ERR)
C— SUB: EDATO - CALLS EDAT1 CASSING ARRAYS
C
COMMON MTOT.NP.IAd)
INCLUDE 'CNTCOM.INC'
IF(EOC) THEN
CALL ALERTCERROR: Insufficient data.'.'$'.'$•)
ERR = .TRUE.
RETURN
ELSEIF(TDATd) .LT.TDATI2)) THEN
CALL ALERTCERROR: Tine data out of sequence. '.'S'.'S')
ERR = .TRUE.
RETURN
ENDIF
CALL GETEDT(KSEQ. E.EDT1.EDT2, ERR)
IF(ERR) RETURN
—2.0 GENERATION TIME LOOP
DO 200 T=TIHE(1),TIME(2),TIME<3)
2.1 UPDATE EXCITATION TIME HISTORY DATA IF NEEDED
20 IF(T.GT.TDATd)) THEN
CALL GETTDT(TDAT)
IF(ECC) THEN
CALL ALERTCERROR: Insufficient data.'.'$','S')
ERR - .TRUE.
RETURN
ELSEIF(TDATd) .LT.TDAT(2» THEN
CALL ALERTCERROR: Time data out of sequence. ','S'.'S')
ERR = .TRUE.
RETURN
ENDIF
CALL GETEDT(KSEQ.E.EDT1.EDT2,ERR)
IF(ERR) RETURN
GO TO 20
ENDIF
2.2 COMPUTE INTERPOLATION FRACTION
XT - (T-TDAT(2))/(TDAT<1)-TDAT(2))
2.3 COMPUTE (EIT))
DO 23 N-l,NSNOD
DO 23 M-l,NSSPE
COMMON /EXCDT/ MPTDAT.MPEDT1, MPEDT2
LOGICAL ERR
CALL EDST1(lA(MPEDTl),N)
RETURN
END
c EDAT1
SUBROUTINE EDAT1(EDT1,N)
-SUB:EDATO - READS EXCITATION TIME HISTORY DATA
INCLUDE 'CNTCOM.INC'
REAL"8 S;DT1 (NSSPE. NSNOD)
CALL FRKERCG'.EDTKI.N),NSSPE)
RETURN
END
SUBROUTINE GENED2(KSEQ.E.TDAT.EDT1.ERR)
-SUB: GENED3 - GENERATES EXCITATION DATA FROM GIVEN TIME DATA
IMPLICIT REAL*8(A-H.O-Z)
COMMON MTOT.NP.IAd)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
EXTERNAL EDATO
C
C GENE02: DATA ( COMMON STORAGE
C
REAL'S TDAT (2), EOT1 (NSSPE.NSNOD), E (NSEQ)
INTEGER KSEQ(NSNOD.NSSPE)
Append -13
-------
Indoor Air Quality Model
Phase II! Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
C—1.0 GET FIRST TIME HISTORY RECORD ( TDAT(l), EDT1 (N3SPE. NSNOD) )
C
CALL GZTTDT(TDAT)
IF(EOC) RETURN
TDAT(2) - TDAT(l)
CALL ZERORIEDT1.NSSPE.NSNOO)
CALL OATGENIEDATO. NSNOD. ERR)
IF(ERR) RETURN
DO 10 N-l,NSNOD
DO 10 M-l.NSSPE
NEQ - ABS(KSEO(N,M))
IFINEQ.NE.O) E(NEQ) - EDTKM.N)
10 CONTINUE
CALL RPRTEN(E(1),IA(KPKSEQ))
WRITE')
CALL FREE
IF (MODE.EQ.'BATCH') CALL FREEHR (NTW)
TIME(l) - O.ODO
TIKE(2) - O.ODO
TIKE(3) - O.ODO
CALL FREERCT',TIME(1),3)
IF(TIME(3).LE.O.ODO) THEN
CALL ALERTCERROR: Tims step must be greater than O.'.'S'.'S'I
CALL ABORT
RETURN
ELSEIF(TIME(2).LT.TIME(1)) THEN
CALL ALERT!
+ 'ERROR: Final time must be greater than Initial time.',
+ '$','$')
CALL ABORT
RETURN
ENDIF
ALPHA - 0.7SDO
CALL FREER CA', ALPHA, 1)
CALL CKRRNGC alpha', ALPHA, 1, O.ODO, 'LELE*. 1. ODO, ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
PINT - 1
CALL FREEK'I'.PINT, 1)
CALL CKIZERCresults print Interval',PINT, 1,'GT',ERR)
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
PSCALE = O.ODO
CALL FREERCS',PSCALE, 1)
IF(ECHO) WRITE(NTH,2250) (TIME(I) , 1=1, 3).ALPHA,PINT
HRITE(NOT,2250) (TIME(I), 1=1, 3).ALPHA.PINT
2250 FORMAT!/.
.' Initial time '.G10.3,/,
.' Final time '.G10.3,/.
.' Time step Increment ',G10.3,/,
.' Integration parameter: alpha .... ',G10.3,/,
.' Reaulta print Interval ',16)
IF(PSCALE.NE.O.ODO) THEN
IF(ECHO) WRITE(NTW. 2260) PSCALE
WRITE(NOT,2260) PSCALE
ENDIF
2260 FORMAT!' Reaulta plot-file scale factor .. \G10.3)
C
C--3.
C
0 DEFINE AND INITIALIZE SYSTEM ARRAYS
IMPLICIT REAL'S(A-H.O-Z)
COMMON MTOT.NP.IAU)
INCLUDE 'IOCCH.INC'
INCLUDE 'CNTCOM.INC'
COMMON /DYNM/ THDAT,TEDAT,MPCDAT
LOGICAL ERR, FOUND
REAL'S TIME(3), PSCALE
INTEGER PINT
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
— D I C T
TIME (3) START TIME, END TIME, TIME INCREMENT
TWDAT
TEDAT
PINT
PSCALE
P 0 I N
MPFS
MPC
MPCD
MPCDD
MPCDAT
KPE
MPID
TIME OF NEXT ELEMENT FLOW RATE RECORD
TIME OF NEXT EXCITATION RECORD
RESPONSE RESULTS PRINT INTERVAL
RESULTS PLOT FILE SCALE FACTOR
TERS TO BLANK COMMON LOC
A T I 0 N S
FS (NSEQ. 2'MSBAN-l) : [F*l DYNAM ALG. MATRIX (ASYM-COMPACT)
C(NSEO) CURRENT (C)
CD(NSEQ) CURRENT d(C)/dt
CDDINSEQ) CURRENT d/dt (d|C)/dt)
CDAT(NSSPE) TEMP. STORAGE OF INITIAL
E(NSEQ) CURRENT (E)
IDCNID) LIST OF INDEPENDENT DOF
CONDS. DATA
EQUATION NOS.
CALL DELETE C TEMP ')
CALL COUNTKIA(MPKSEQ) ,NID)
CALL DELETE ('ID
CALL DELETE CFS
CALL DELETE ('CD
CALL DELETE CCDD
CALL DELETE! 'COAT
CALL DELETE I 'CONT
CALL DELETE ('VCD
CALL DELETE CEFF
CALL DELETE CDIFF
CALL DELETE CGENR
CALL DELETE ('WE
CALL DELETE CC
CALL DELETE) 'E
CALL DELETE! 'F
CALL DELETE ('VM
CALL DEFINRCVM
CALL DEFINRCF
CALL DEFINRCE
CALL DEFINRCC
CALL DEFINRCWE
CALL DEFINR CGENR
CALL DEFINR CDIFF
CALL DEFINR CEFF
CALL DEFINR! 'VCD
CALL DEFINR! 'CONT
CALL DEFINR CCDAT
CALL DEFINR ( 'COD
CALL DEFINR! 'CD
CALL DEFINR CFS
CALL DEFINICID
, HPVM.NSEQ, 1)
, KPF. NSEO. 2*MSBAN-1)
.MPE.NSEO, 1)
.MPC.NSEQ.l)
.MPHE.NFCLM.l)
.MPGENR.NSSPE.l)
. MPDIFF.NSSPE.il
.MPEFF.NSSPE.l)
.MPVCD, NSNOD. 1)
.MPCONT, NSNOD. 1)
, MPCDAT, NSSPE.l)
, MPCDD.NSEQ. 1)
, MPCD.NSEQ. 1)
, MPFS.NSEO. 2'MSBAN-l)
.MPID.NID.l)
ERR - .FALSE.
HRITE(NOT.2000)
C
C— 4
C
0 GET NODAL INITIAL CONCENTRATIONS
CALL ZEROR(IA(KPC).NSEQ. 1)
CALL CETIC(ERR)
Append -14
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
IF(ERR) THEN
CALL ABORT
RETURN
ENDIF
C
C--5.0 OPEN ELEMENT, FLOW AND EXCITATION DATA FILES, ( PLOT FILE
C
OPEN(ND1,FILE"(FNAME(1:LFNAME)//'.FEL'),STATUS-'OLD',
+FORM-'UNFORMATTED•)
IF(NKINEL.GT.O) THEN
OPEN(ND2,FILE-(FNAME(1:LFNAME)//'.KIN'),STATUS-'OLD'.
+ FORM-'UNFORMATTED')
ENDIF
OPEN
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
C (I.E.. TDOF EQUATION NUMBERS)
C X(NEQN.2*MBAN-1) : [K] MATRIX: ASYM-BANDED COMPACT-STORED
C KSINEQN, 2'MBAN-l) : IK'I = [C| + aDT[K| MATRIX (SCALED FOR NEC ID)
C C(NEQN) : CURRENT (C] (ORDERED BY EQTN I)
C E(NEQN) : CURRENT (E) (ORDERED BY EQTN I)
C TINEQN) : CURRENT (T) (ORDERED BY EQTN I)
C TD(NEQN) : CURRENT ERROR: Can not compute for step change In',
+ ' dependent variable number:'.IS)
ERR - .TRUE.
RETURN
ELSE
roil) - O.ODO
ENDIF
c 'E'-DOF: FORM [E|-[K1(T) WHERE [K] IS IN COMPACT STORAGE
ELSE
TEMP - Ed)
Kl - MAXd.MBAN-I+l)
K2 - MIN(2'MBAN-l,MBAN+NEQN-I)
DO 20 KK-K1.K2
J - I * KK - MBAN
TEMP - TEMP - K(I,KK)«T(J)
20 CONTINUE
C SOLVE
TO (I) - TEMP/C (I)
22 ENDIF
C
C—3.0 COMPUTE TAYLOR'S TIMESTEP CHECK
C
IF(ECHO) WRITE(NTW. 2300)
WRITE(NOT.2300)
2300 FORHAT
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
COMMON /DYNM/ TWDAT, TEDAT, MPCDAT
REAL'S KINSEQ,2'MSBAN-l), C(NSEQ).TM, TWDAT, TEDAT
LOGICAL ERR. KUPDAT
C
C—1.0 UPDATE ELEMENT FLOH RATES IFITH.GE.TNDAT)
C
CALL UPDAT .' Time: ',010.3)
CALL RPRTEN(E,IA(MPKSEQ)>
ENDIF
RETURN
END
------------------------------------------------------------------ UPDAT
SUBROUTINE UPDAT (LUN. T. TD, D. ND, UPDATE. ERR)
—SUB: UPDAT
SEARCHES A SEQUENTIAL DATA RECORD, ON UNIT LUN, OF THE FORM:
TD
(D (I) ,1=1, NO)
TD
TO UPDATE DATA VALUES TO CURRENT TIME. "T" . IF DATA VALUES ARE
UPDATED LOGICAL -UPDATE" IS SET TO TRUE.
: DISCRETE TIME VALUE
: UPDATED TO NEXT VALUE
: CORRESPONDING DISCRETE DATA VALUES
TD
DID
UPDAT MUST BE "PRIMED" BY READING FIRST TD VALUE TO MEMORY
INCLUDE 'IOCOM.INC'
REAL'S D(ND).T,TD
LOGICAL ERR, UPDATE
UPDATE = .FALSE.
10 IF(T.GE.TD) THEN
UPDATE DISCRETE DATA VALUES
READ(LUN, ERR-800, END=900I (D(I).1=1,NDI
IF(ERR) RETURN
UPDATE o .TRUE.
GET NEXT DISCRETE TIME
READ (LUN. ERR=800. END=900) TD
IF (ERR) RETURN
GO TO 10
ELSE
RETURN
ENDIF
800 ERR = .TRUE.
CALL ALERTCERROR: Time historydata file read error. ','$'.' $')
RETURN
900 ERR = .TRUE.
CALL ALERT (
+'ERROR: EOF encountered on time history data file.',
+'Insufficient tlras history data.'.'$')
RETURN
END
IMPLICIT REAL'S (A-H.O-Z)
REAL'S T (NEON) ,TD (NEQN) ,E (NEQN).K(NEON,2-MBAN-lI.
+KS(NEON,2«MBAN-1)
INTEGER ID (NID)
LOGICAL TDOF
DO 20 I..I.NEQN
C SCALE BV DIAGONAL FOR TEMP PRESCRIBED NODES
IF(TDOF(I,ID, NID)) THEN
TO (I) » TD(I)*KS(I,MBAN)
C FORM (EI-tKUT) WHERE IK| IS IN COMPACT STORAGE
ELSE
TEMP •• E(I)
IO. - HAXU.MBAN-I+1)
K2 - I«N(2«MBAN-1.MBAN+NEQN-I)
DO 10 KX°K1,K2
J - I + KK - MBAN
TEMP '- TEMP - K(I,KK)*T(J)
10 CONTINUE
TD(I) - TEMP
ENDIF
20 CONTINUE
RETURN
END
FUNCTION TDOF(N,ID,NID)
C-FUN:TDOF - DETERMINES IF EQUATION NUMBER N IS A TEMPERATURE DOF
LOGICAL TDOF
INTEGER ID(NID)
TDOF - .FALSE.
DO 10 NN-l.NID
IF«ID(NN).EQ.N)) THEN
TDOF - .TRUE.
RETURN
ENDIF
10 CONTINUE
RETURN
END
RESULT
SUBROUTINE FORMKS(ID,K.KS.C.ALPHA. DT,NID,NEQN,MBAN)
-SUB:FORMXS - FORMS;
[K*| = [Cl + aDT(K|
SCALES [K«] = [K*1*1.0D15 FOR 'T'-DOF
IMPLICIT REAL'S(A-H.O-Z)
REAL'S K(NEON.2*KBAN-1),KS(NEQN,2'MBAN-l).C(NEQN)
INTEGER ID(NID)
SUBROUTINE RESULT(T)
C—SUB:RESULT - REPORTS RESPONSE RESULT VECTOR (T)
COMMON MTOT.NP.IA(l)
INCLUDE 'CNTCOM.INC'
REAL'S T(NSEQ)
CALL RPRTENIT,IAIMPKSEQ))
RETURN
END
C • •*-<—'-< ' ' ^=—==0=. RESET
SUBROUTINE RESET
C—SUB:RESET - COMMAND TO RESET CONTAM BY RE-INITIALIZING POINTERS AND
C COUNTERS AND DELETES ARRAYS LEFT BY CONTAM IN BLANK COMMON
C--HELP LIST
C .' RESET Reset CONTAM for new problem.'
c
INCLUDE 'IOCOM.INC'
LOGICAL FOUND
CHARACTER BINEXT(4)*3
INTEGER NBIN
DATA BINEXT/'FEL','KIN'.'WDT','EOT'/, NBIN/4/
C
C—1.0 RE-INITIALIZE CONTAM CONTROL VARIABLES t DELETE CONTAM ARRAYS
C
CALL ::NITCN
CALL DELETE! '0 ')
CALL DELETE('VM ')
CALL DELETECC ')
CALL DELETE!'F ')
CALL DELETE('WE ')
CALL DELETE!'V ')
CALL DELETE CKSEO'I
Append -17
-------
Indoor Air Quality Model
Phase III Report
CONTAM87 FORTRAN
APPENDIX
77 Source Code
DO 100 NK=1.»
100 CALL DELETE('KIK'//CBAR(NK+48))
CALL DELETE ('TEMP')
CALL DELETE CCONT')
CALL DELETE ('VCD ')
C
C — 2.0 DELETE CONTAM BINARY FILES
C
DO 20 N-l.NBIN
INQUIRE (FILE-FNAME (1 :LFNAME) //BINEXT
RETURN
END
c READWE
SUBROUTINE READWE(ERR)
C—SUB:README - READS I REPORTS ELEMENT TOTAL MASS FLOW RATE DATA
C
COMMON MTOT.NP.IAU)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
EXTERNAL NEDATO
IF(ECHO.OR.(MODE.EQ.'INTER')) HRITE(NTH,2000)
HRITE(NOT.2000)
2000 FORMAT!/,' -= Element Mass Flow Rates')
CALL DATGENIHEDATO.NFELM.ERR)
IF(ERR) RETURN
CALL RPRTNO(IA(MPHE),NFELM,'Elem')
RETURN
END
SUBROUTINE HEDATO(N,ERR)
C—SUB:HEDATO - CALLS HEDAT1 PASSING ARRAYS
C
COMMON MTOT.NP.IAd)
INCLUDE 'CNTCOM.INC'
Append -18
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
DO 100 NX-1,9
100 CALL DELETE('XIX'//CHAR 0.',
+ 'FLOHSYS command must be executed.','$')
ERR - .TRUE.
RETURN
ENDIF
CALL LOCATE('KSEQ',MPXSEQ.NR,NO
IF(MPKSEQ.EQ.O) THEN
CALL ALERT!
+ 'ERROR: System equation number army "KSEQ" not found.',
+ 'FLOHSYS command must be executed.','$')
ERR - .TRUE.
RETURN
ENDIF
CALL LOCATE CV '.MPV.NR.NO
IF(MPV.EQ.O) THEN
CALL ALERT('ERROR: Nodal volumetric mass array ~V" not found.
+ 'FLOHSYS command must be executed.','$')
ERR - .TRUE.
RETURN
ENDIF
C
C— 2
C
IF(NOPT.EQ.l) RETURN
0 FLOHELEM I XINELEM DATA VERIFICATION
IF((NFELM.EQ.O).AND.(NKINEL.EQ.0)) THEN
CALL ALERT!
+ 'ERROR: Number of flow t kinetics elements both = 0.',
+ 'FLOHELEM C/or KINELEM must be executed.'.'$')
ERR = .TRUE.
RETURN
ELSEIF(NFELM.EQ.O) THEN
CALL ALERTCERROR: Number of flow elements - 0.',
+ 'FLOHELEM must be executed.'.'$')
ERR - .TRUE.
RETURN
ELSEIF(NKINEL.EO.O) THEN
IF (ECHO) WRITE (NTH, 2210)
WRITE(NOT.•)
WRITE(NOT,2210)
0 FORMAT(
+ ' •* NOTE: Number of kinetics elements - 0.')
ENDIF
INQUIRE(FILE-(FNAME(1:LFNAME)//'.FEL'),EXIST-FOUND)
IF(.NOT.FOUND) THEN
CALL ALERT(
+ 'ERROR: Data file '//FNAME(1:LFNAME)//'.FEL not found.1,
+ 'FLOHELEM command must be executed.'.'$')
ERR - .TRUE.
RETURN
ENDIF
IF(NKINEL.ST.0) THEN
INQUIRE(FILE-(FNAME(1:LFNAME)//'.KIN'),EXIST-FOUND)
IF(.NOT.FOUND) THEN
CALL ALERT!
+ 'ERROIl: Data file '//FNAME(1:LFNAME)//'.KIN not found.'
+ 'KINELEM command must be executed.','$')
ERR - .TRUE.
RETURN
ENDIF
ERR - .TRUE.
DO 200 IIK-1,9
CALL LOCATE('XIX'//CHAR(48+NK), MPKIK(NK), NR.NC)
IF(MPKIIUNK).NE.O) ERR- .FALSE.
IF(ERR) CALL ALERT(
+ 'ERROR: Kinetics rate coefflcent arrays not found.',
+ 'KINELEM command must be executed.','$')
RETURN
ENDIF
IF(NOPT.EQ.2) RETURN
C— 3.0 FLOHDAT DATA VERIFICATION
C
INQUIRE(FILE-(FNAME(1:LFNAME)//'.TOT'),EXIST-FOUND)
IF(.NOT.FOUND) THEN
CALL AlERTI
+ 'ERROR: Data file '//FNAME(1:LFNAME)//'.HOT not found.',
+ 'FLONDAT command must be executed.','$')
ERR - .TRUE.
RETURN
ENDIF
C
C— 4.0 EXCITDAT DATA VERIFICATION
C
INQUIRE(FILE-(FNAME(1:LFNAME)//'.EOT').EXIST-FOUND)
IF!.NOT.FOUND) THEN
CALL ALERT!
+ 'ERROR: Data file '//FNAME(1:LFNAME)//'.EOT not found.',
+ 'EXCITDAT command must be executed.','$')
ERR - .TRUE.
RETURN
ENDIF
IF(NOPT.EO.3) RETURN
RETURN
END
SUBROUTINE GETIDS(IDS,NIDS,LABEL)
C—SUB: GETIDS - GETS CHARM IDS FROM COMMAND/DATA LINE
C SETS ID" NUMBER FOR BLANK IDS
INCLUDE 'IOCOM.INC'
CHARACTER IDS (NIDS) M, LABEL*!*), HEADER* (*)
PARAMETER (HEADER-'Nun. ID')
CALL FRLEC('D',IDS<1),4,NIDS)
DO 10 N-l.NIDS
IF (IDS («) .EQ.' ') THEN
NCENT - N/100
NTENS - (N - NCENT*100)/10
NONES - N - NCENTMOO - NTENSMO
IDS(N) - ' '//CHAR(NCENT+48)//CHAR(NTENS+48)//CHAR(NONES+4e)
10 ENDIF
NCOLS - MIN(NIDS.S)
IF(ECHO) THEN
WRITE(NTH,2000) LABEL, (HEADER,N-l,NCOLS)
HRITE(NTH,2010) (I, IDS(I),I-l.NIDS)
ENDIF
WRITE(NOT,2000) LABEL, (HEADER,N-l,NCOLS)
WRITE(NOT,2010) (I,IDS(I),1-1,NIDS)
2000 FORMAT!/,' — ', (A),//.8X.5(: (A),3X))
2010 FORMAT((8X,5(:I3,2X,A4,2X)))
RETURN
END
SUBROUTINE README(ERR)
C—SUB:README - READS < REPORTS ELEMENT TOTAL MASS FLOW RATE DATA
C
COMMON MTOT.NP.IAU)
INCLUDE 'IOCOM.INC'
INCLUDE 'CNTCOM.INC'
LOGICAL ERR
EXTERNAL HEDATO
IF(ECHO.OR.(MODE.EO.'INTER')) WRITE(NTH.2000)
WRITE(KOT,2000)
2000 FORMAT(/,' — Element Mass Flow Rates')
CALL DJiTGEN(HEDATO,NFELM, ERR)
IF(ERR) RETURN
CALL RI'RTNO(IAIMPHE),NFELM,'Elem')
RETURN
END
SUBROUTINE HEDATO(N, ERR)
C—SUB:HEDATO - CALLS HEDAT1 PASSING ARRAYS
C
COMMON MTOT.NP.IA(l)
INCLUDE 'CNTCOM.INC'
Append -18
-------
Indoor Air Quality Model APPENDIX
Phase III Report CONTAM87 FORTRAN 77 Source Code
HIN - MAX C NSP : (CURRENT) SPECIES NUMBER
DO 10 I-l.NENOD C FORM : FORM OF SYSTEM ARRAY 'FULL1 OR 'BAND-
CD 10 J-l.NESPE C CONT(NSNOD) : NODAL MASS CONTINUITY ACCUMULATOR
NN - ABS(KSEQ(LN(I),LS(J») C LN : (LOCATION) NODE Of KINETICS ELEMENT
IF(NN.GT.MAX) MAX-NN C
IFINN.LT.MIN) MIN-NN
10 CONTINUE C
C C—0.0 INITIALIZE SYSTEM ARRAYS
C—2.0 COMPUTE ELEM. BADWIDTH AND COMPARE TO CURRENT MAX SYST. BANDWIDTH C
C IF(FORM.EO.'BAND') CALL ZEROR(F.NSEQ,2*MSBAN-1)
MEBAN - MAX-MINtl IF (FORM.EQ.'FULL') CALL ZEROR(F.NSEQ. NSEQ)
IF(MEBAN.GT.MSBAN) MSBAN-MEBAN C
RETURN C—1.0 PROCESS FLOH ELEMENTS
END C
C 1.1 FORM AND ASSEMBLE ELEMENT ARRAYS
SUBROUTINE RPRTNO(X,NX,LABEL) REWIND (ND1)
C~SUB:RPRTNO - REPORTS REAL VECTOR(X) BY INDEX NUMBER CALL ZEROR(IA (MPVCD) ,NSNOD. 1)
CALL ZEROR(IA(MPCONT),NSNOD.l)
C VARIABLE DESCRIPTION DO 10 NEL-1.NFELM
C X(NX) VECTOR OF REAL VALUES ORDERED BY INDEX NUMBER READ(ND1, ERR-900, END-900) TYPE
C LABEL TABLE LABEL CHARACTER**
C IF(TYPE.EQ.'SIMP') THEN
CALL SIMP(NEL.WE,lA(MPEFF).KSEQ.F.IA(MPCONT).FORM.ERR)
IMPLICIT REAL'S(A-H.O-Z) IF(ERR) RETURN
ELSEIF(TYPE.EO.'CNDF') THEN
INCLUDE 'IOCOM.INC' CALL CNDFINEL.HE, lA(MPDIFF) . lA(MPGENR) .KSEO.F.IA(MPVCO) ,
+ IA(MPCONT).FORM.ERR)
REAL'S X(NX) IF(ERR) RETURN
CHARACTER LABEL'4 ELSE
GO TO 900
WRITE (NOT. 2000) (LABEL.N-1.4) 10 ENDIF
IF(ECHO) WRITE(NTH,2000) (LABEL.N-1.4) C
WRITE (NOT. 2010) (N, X(N), N-l.NX) C 1.2 REPORT NET TOTAL MASS FLOW
IF (ECHO) WRITE (NTW, 2010) (N, X(N), N-l.NX) C
WRITE(NOT,2200)
2000 FORMAT!/,6X.4I2X.A4,' Value'. 3X1) IF (ECHO) WRITE (NTW, 2200)
2010 FORMAT((6X,1(I6.1X,G11.3)H 2200 FORMAT(/.' — Nee Total Mass Flow')
CALL RPRTNO(IA(MPCONT),NSNOD,'Node')
RETURN C
END C—2.0 PROCESS KINETICS ELEMENTS
C '
SUBROUTINE RPRTEN(X.KSEQ) C 'TEMP' STORES SPECIES CONNECTIVITY ARRAY, LM(NSSPE)
C—SUB:RPRTEN - REPORTS(X)IN NODE ORDER SEQUENCE FOR EACH SPECIES
CALL DELETE('TEMP')
C VARIABLE DESCRIPTION CALL DEFINR (• TEMP', MPTEMP, NSSPE, 11
C X(NSEQ) VECTOR OF VALUES ORDERED BY EQUATION NUMBER
IMPLICIT REAL'8(A-H.O-Z) DO 200 NEL=1.NKINEL
READ(ND2. ERR-950. END-950) LN, NK
INCLUDE 'IOCOM.INC' 200 CALL KINELX (LN, KSEQ, IA(MPKZK(NK)), F. IA (MPTEMP). IA (MPV). FORM)
INCLUDE 'CNTCOM.INC'
CALL DELETECTEMP')
REAL'S X (NSEQ), XX M)
INTEGER KSEQ (NSNOD, NSSPE) RETURN
CHARACTER FLG(4)*1 C
C—3.0 READ ERROR TERMINATION
WRITE(NOT,2000) C
IF(ECHO) WRITE(NTW,2000) 900 ERR = .TRUE.
2000 FORMAT!/, CALL ALERT(
.13X,'V - independent DOFa "U" - undefined DOFa.'l +'ERROR: Read error or EOF In file '//FNAME//'.FEL','$','$')
RETURN
DO 100 M-l,NSSPE
WRITE(NOT.2010) SID(M) 950 ERR - .TRUE.
IF (ECHO) WRITE (NTW, 2010) SID (M) CALL ALERT (
2010 FORMAT!/,8X,'Species: '.A4,/. *'ERROR: Read error or EOF In file '//FNAME//'.KIN','$','$')
* 6X,4(2X,'Node Value',3X1) RETURN
DO 100 N=l,NSNOD,« END
NN = MIN(N«3, NSNOD)
DO 10 I-N.NN.l C SIMP
NEQ - KSEQII.M) SUBROUTINE SIMP (NEL,WE, EFF, KSEO.F.CCNT, FORM.ERR)
NNEO * ABS(NEO) C—SUB:SIMP - FORMS AND ASSEMBLES SIMPLE FLOH ELEMENT EQUATIONS
IF(NEQ.LT.O) THEN C FOR ALL SPECIES CONSIDERED
XXd-N+1) - X(NNEQ)
FLG(I-N*1) - '" INCLUDE 'IOCOM.INC'
ELSEIF(NEQ.EQ.O) THEN INCLUDE 'CNTCOM.INC'
XXU-N+1) • O.ODO
FLGd-N+1) - 'U' REAL'S WE (NFELM), EFF (NSSPE) ,F (NSEQ, 1) ,CCNT(NSNOD) , ELFI2, 2) ,W
ELSE INTEGER KSEO(NSNOD.NSSPE).LN(2>, LM(2)
XX(I-N-l-l) - X(NNEQ) CHARACTER FORM'4
FLG(I-N»1) - ' ' LOGICAL ERR
ENDIF
10 CONTINUE C VARIABLE DESCRIPTION
IF(ECHO) WRITE(NTH,2020) (I,FLG(I-Ntl),XX(I-N+1). I=N, NN| C NEL : (CURRENT) ELEMENT NUMBER
WRITE(NOT,2020) (I.FLO(I-N+1).XX(I-N+1),I=N,NN) C FORM : FORM OF SYSTEM ARRAY 'FULL' OR 'BAND'
2020 FORMAT((6X,4(K,1A1,G11.3))) C CONT(NSNOD) : NODAL MASS CONTINUITY ACCUMULATOR
100 CONTINUE C EF (NENOO, NENOD) : ELEMENT (F| ARRAY
C W : ELEMENT TOTAL MASS FLOW RATE
RETURN C LN(2) : ELEMENT NODE LOCATION/CONNECTIVITY
END C NSP : (CURRENT) SPECIES NUMBER
C LM(2) : SYSTEM DOF CORRESPONDING TO EACH ELEMENT DOF
SUBROUTINE FORMF(KSEQ,F.WE.FORM,ERR)
C~SUB:FORMF - FORMS SYSTEM FLOW MATRIX AND CONDF CONTRIBUTION TO [V] C
C ARRAY CONT(NSEQ) USED TO CHECK NODAL MASS FLOH CONTINUITY C—1.0 GET ELEMENT DATA
C
COMMON MTOT.NP.IA(l) READ(ND1. END=900, ERR=900) LNU) .LN (2), (EFF(I) , 1-1,NSSPE)
W = HE (NEL)
INCLUDE 'IOCOM.INC' C
INCLUDE 'CNTCOM.INC' C—2.0 FORM ELEMENT ARRAYS
C
REAL'S F(NSEQ,1), HE (NFELM) IF (W. GT. 0 . ODO) THEN
INTEGER KSEQ(NSNOD.NSSPE). MPCONT ELF(1.1) °W
LOGICAL ERR, LN . ELF(1,2) - O.ODO
CHARACTER FORM«4, TYPE*4 , ELFI2.2) • O.ODO
CONT(LNU)) = CONT(LNd)) + H
C VARIABLE DESCRIPTION CONT(LN(2» -CONT(LN(2» -W
C NEL : (CURRENT) ELEMENT NUMBER ELSEIF(H.LT.0.ODO) THEN
Append -20
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
ELF(1,1) - O.ODO
ELF(2,1) - O.ODO
ELF(2,2) - -W
CONT(LNU)) - CONT(LNU)) + W
CONT(LN(2)I - CONT(LN(2)) - W
ELSE
ELFU.l) - O.ODO
ELF(1,2) - O.ODO
O.ODO
O.ODO
ELF(2.1)
ELFI2.2)
ENDIF
10 CONTINUE
RETURN
900 ERR - .TRUE.
CALL ALERT!
•(•'ERROR: Read error or EOF In file •//FNAHE//'.FEL'.'S'
RETURN
END
'$•>
C 2.1 LOOP OVER SPECIES FOX NONZERO OFF-DIAGONAL TERM
C
DO 10 NSP-1,NSSPE
LH(1) - ABS(XSEQ(LN(1),NSP))
LM(2) - ABS(KSEQ(LN(2).NSP))
IF(W.GT.O.ODO) THEN
ELF(2,1) - -«•(!.ODO-EFFfliSP))
ELSEIFW.LT. O.ODO) THEN
ELF (1.2) - NMl.ODO-EFF(NSP))
ENDIF
SUBROUTINE KJNELK(LN,XSEO,XIX,F,LM,V,FORM)
C~SUB:XINELX - FORMS KINETICS ELEMENT ARRAY FROM KIN RATE COEF. MATRIX
C
LN : (LOCATION) NODE OF KINETICS
XIK(NSSPE.NSSPE): KINETICS RATE COEF. MATRIX
FINSEQ,1) : SYSTEM FLOW MATRIX
LM(NSSPE) : SYSTEM DOF CORRESPONDING TO EACH ELEMENT DOF
V(NSNOO) : NODAL VOLUMETRIC MASSES
FORM : TORM OF [F]: 'FULL' OR 'BAND'
C
C 2.
C
2 ASSEMBLE ELEMENT ARRAYS
CALL ADDA(ELF,2,F.NSEQ.HSBAN.LM,1.ODD,FORM)
10 CONTINUE
RETURN
900 ERR - .TRUE.
CALL ALERT!
+ 'ERROR: Read error or EOF In file '//FNAME//'.FEL', '$','$•)
RETURN
END
SUBROUTINE CNOF(NEL,HE,DIFF,GENR,KSEQ,F,VCD,CCNT,FORM,ERR)
C—SUB:CNDF - FORMS AND ASSEMBLES CONV-DIFF FLCH ELEMENT EQUATIONS
C FOR ALL SPECIES CONSIDERED
INCLUDE
INCLUDE
'IOCOM.INC'
'CNTCOM.INC'
INCLUDE 'CNTCOM.INC'
REAL'S laK (NSSPE. NSSPE), F(NSEQ.l), V(NSNOD), SCALE
INTEGER 1.M (NSSPE) , LN, KSEQ (NSNOO. NSSPE)
CHARACTER FORUM
DO 100 N»l,NSSPE
100 LM(N) - ABS(KSEQ(LN,N))
SCALE - V(LN)
CALL ADDA(KIK,NSSPE,F,NSEQ,MSBAN,LM,SCALE,FORM)
RETURN
END
CNDF SUBROUTINE ADDA (AE.NEDOF, AS,NSDOF.MSBAN.LM, SCALE, FORM)
C—SUB:ADDA - ADDS SCALED ELEMENT ARRAY, SCALE*[AEI, TO SYSTEM ARRAY, AS
C
REAL'S AE(NEDOF.NEDOF), AS(NSDOF,1). SCALE
INTEGER LM(NEDOF)
CHARACTER FORM'4
IMPLICIT REAL*8(A-H.O-Z)
REAL*! HE(NFELM),DIFF(NSSPE),GENR(NSSPE),F(NSEQ,1),VCD(NSNOD),
+CONT(NSNOD),ELF(2,2),H,MASSL,LENGTH,FACTOR
INTEGER KSEQ(NSNOO.NSSPE).LN(2),LM(2)
CHARACTER FORM*4
LOGICAL ERR
C VARIABLE DESCRIPTION
C NEL : (CURRENT) ELEMENT NUMBER
C VARIABLE DESCRIPTION
C AE(NEDOF.NEDOF) : ELEMENT ARRAY
C NEDOF : NUMBER OF ELEMENT DOFS
C AS(NSDOF,2'MSBAN-l) : (COMPACTED) BANDED ASYM. SYSTEM ARRAY
— OR --
AS (NSDOF, NSDOF) : FULL ASYM. SYSTEM ARRAY
LM (NEDOF) : SYSTEM DOF CORRESPONDING TO EACH ELEMENT DOF
SCALE : SCALAR FACTOR
FORM : FORM OF SYSTEM ARRAY 'FULL' OR 'BAND'
C
c
c
c
c
c
c
c
c
c
FORM : FOR
CCNT (NSNOD)
ELF(NENOD,NE>
N : ELE
LN(2)
NSP : (CU
LM(2)
MASSL
LENGTH
FACTOR
M OF SYSTEM ARRAY 'FULL' OR 'BAND'
NODAL MASS CONTINUITY ACCUMULATOR
CO) : ELEMENT [F| ARRAY
MENT TOTAL MASS FLOW RATE
ELEMENT NODE LOCATION/CONNECTIVITY
RXENT) SPECIES NUMBER
SYSTEM DOF CORRESPONDING TO EACH ELEMENT DOF
AIR MASS PER UNIT LENGTH - CNDF ELEMENTS
FLOH PASSAGE LENGTH - CNDF ELEMENTS
UPWIND FACTOR - CNDF ELEMENTS
c
C— 1.0 GET ELEMENT DATA
C
DO 20 I«l.NEDOF
II - LM(I)
DO 10 J-l. NEDOF
IF(FORM.EQ.'BAND') JJ = MSBAN - II + LM(J)
IF(TORM.EQ.'FULL') JJ-LM(J)
ASCII,JJ) = AS(II.JJ) + SCALE'AEd, J)
10 CONTINUE
20 CONTINUE
RETURN
END
REAO(ND1.END-900,ERR=900) LN (1) ,LN (2) .MASSL. LENGTH. FACTOR.
+ (DIFF (N).N-l, NSSPE)
W - WE (NEL)
C
C— 2.0 FORM ELEMENT ARRAYS
C
IF (W.LT. O.ODO) THEN
LNTEMP - LNI2)
LN(2) - LNU)
LN(1) » LNTEMP
ENDIF
C
C ---- 2.1 FORM ELEMENT LUMPED VOLUMETRIC MASS TERMS
C
VCD(LN(1)) - VCD(LNU))
VCD (LN (2)1 = VCD(LN(2)I
MASSL*LENGTH*0.50DO
MASSL'LENGTH'O . 50DO
C
C ---- 2.2 ACCUMULATE MASS FLOW RATES FOR CONTINUITY CHECK C-
C
CONT(LNIl)) -CONT(LNU)) + W C-
CONT(LN(2») -CONT(LN(2)I -W C
C
C ---- 2.3 LOOP OVER SPECIES TO FORM ELEM MASS TRANSPORT MATRIX 4 ASSEMBLE
C
DO 10 NSP-1, NSSPE C
LM(l) = ABS(XSEQ(LN(1),NSP)I C-
LM(2) = ABS(KSEQ(LN(2) ,NSP)) C
COEF - MASSL'DIFF (NSP) /LENGTH C
IF(FACTOR.EQ.-l.ODO) FACTOR>MAX (0. ODO. 1 .ODO-COEF/W) c-
ELF(l.l) = (W/2)*< 1 + FACTOR) + COEF c
ELF (1,2) = (W/2)*( 1 - FACTOR) - COEF c.
ELFI2.1) » (W/2)*(-l - FACTOR) - COEF
ELF (2, 2) - (W/2JM-1 * FACTOR) + COEF
—-2.2 ASSEMBLE ELEMENT FLOW ARRAYS
CALL ADDA(ELF,2,F.NSEQ,MSBAN,LM,1.ODO,FORM)
SUBROUTINE FORMVM(KSEO,V,VCD,VM)
-SUB:FORMVM - FORMS VOLUMETRIC MASS MATRIX (A DIAGONAL ARRAY)
INCLUDE 'CNTCOM.INC'
REAL'S V(NSNOD) , VM (NSEQ), VCD (NSNOD) ,VN
INTEGER KSEQ(NSNOD.NSSPE)
CALL ZEROR(VM,NSEQ,1)
00 10 N-l,NSNOD
VN - »(N) + VCD(N)
DO 10 M-l,NSSPE
NEQ - ABSIKSEQ(N.M))
IF(NEQ.NE.O) VM(NEQ) - VN
10 CONTINUE
RETURN
END
SUBROUTINE GETTDT(TDAT)
-SUB: GETTCO - UPDATES TIME DATA VALUES
INCLUDE 'IOCOM.INC'
REAL'S TDATI2)
-1.0 UPDATE OLD VALUES
TDATI2) = TDAT(l)
-2.0 READ NEW VALUE
— CHECK TOR END-OF-COMMAND "END"
IF(EOC) THEN
EOD •• .TRUE.
RETURN
ENDIF
IFIHODK.EQ.'INTER') CALL PROMPT('TIME>')
CALL FREE
IFIMOOE.EQ.'BATCH') CALL FREEWR(NTW)
Append -21
-------
Indoor Air Quality Model APPENDIX
Phase III Report CONTAM87 FORTRAN 77 Source Code
C ---- CHECK FOR END-OF-COMMAND "END" CALL FREEIC '.NR,1)
IF(EOC) THEM GO TO 100
EDO - .TRUE. C
RETURN 200 IF(NC.GT.O) GO TO 900
ZNDIF CALL PROMPT C " Enter number of columns: ')
CALL FREER('E',TDAT(1),1) CALL FREE
C ---- REPORT CALL FREEIC '.NC,D
IF (ECHO) WRITE (NTH, 2020) TDAT(l) GO TO 200
WRITE (NOT, 2020) TDAT(l) C
2020 FORMAT!/,' —Tina: '.G10.3) 900 RETURN
END
RETURN
END C ----------------------------------------------------------------- ABORT
SUBROUTINE ABORT
C — SUB: ABORT - ABORTS COMMAND AND RETURNS TO INTERACTIVE MODE
C C
C COMMAND PROCESSOR UTILITIES C INCLUDE 'IOCOM.INC'
C C
C+-M-++++++++++-M-M * * ' * M+-f+4-»-f-M-»-f-f+-f-f-M"f-f-M"f+-M"f-f+-H"»-»-*-f4"f-f+-f+-M-f-f-M-4"fC WRITE (NTW, 2000)
WRITE (NOT, ")
C -------------------------------------------------------------- NOPEN WRITE (NOT, 2000)
SUBROUTINE NOPEN (LUN.FNAME. FRM) 2000 FORMAT!' **•* COMMAND ABORTED')
C — SUB: NOPEN - OPENS A FILE AS A NEW FILE WHETHER IT EXISTS OR NOT IF (MODE. EQ. 'BATCH' ) CALL RETRN
C LUN •> LOGICAL UNIT NUMBER
C FNAME - FILENAME RETURN
C FRH - FORM; 'UNFORMATTED' OR 'FORMATTED' END
INTEGER LUN C+++-f++++4-+++4^++-
CHARACTER FNAME* (•), FRM«(«> C C
LOGICAL FOUND C CALSAPX LIBRARY C
C C
INQUIRE (FILE-FNAME.EXIST-FOUND) C AN EXTENSION OF "CAL-SAP" LIBRARY OF SUBROUTINES C
IF (FOUND) THEN C DEVELOPED BY ED WILSON, U.C. BERKELEY C
OPEN (LUN, FILE=FNAME.STATUS='OLD',FORM=FRM) C*++«++«+++«++«++«++«++-M++++++-M.+»-M-++++++t+++-M.++++++<-f++-H-++++C
WRITE (LUN, 2000) LUN C 1.0 FREE-FIELD INPUT SUBROUTINES
2000 FORMAT (1 6) C+++^++^+++^++4+++4>++^++^++^++4^++^++^++4+++4+++^++^+++-*-++-t~K:
ELSEIFIFRM.EQ. 'UNFORMATTED') THEN C ---------------------------- - ------------------------------------- FREE
WRITE (LUN) LUN SUBROUTINE FREE
ENDIF C~SUB:FREE - READ LINE OF FREE FIELD DATA
CLOSE(LUN,STATUS-'DELETE') C COMMENTS LINES ECHOED TO SCREEN
OPEN (LUN, FILE-FNAME. STATUS-'NEW' , FORM-FRM)
ELSE INCLUDE 'IOCOM.INC'
OPEN (LUN. FILE-FNAME, STATUS-'NEW', FORM=FRM) INCLUDE 'FRECOM.INC'
ENDIF C
C-O.O-INITIALIZE VARIABLES
RETURN C
END EOD - .FALSE.
EOC - .FALSE.
SUBROUTINE APPEND(LUN) 3 LINE(I)-' '
C--SUB: APPEND - POSITIONS 'OLD' FILE AT LAST RECORD SO ADDITIONAL C
C RECORDS MAY BE APPENDED C-1.0 GET LINE OF DATA
C LUN - LOGICAL UNIT NUMBER C
INTEGER LUN II- 80
READ(NC«D.1000.ERR=100> (LINE (KX) .XK-I.II)
REWIND LUN
10 READ(LUN,*,END-20) C ----- CHECK FOR ADDITIONAL LINE
GO TO 10
20 BACKSPACE (LUN) JJ - LENTRM(LLINE)
RETURN DO 12 K-I.JJ
END IF(LINE(K).EQ. '\'( THEN
I - K
c --------------------------------------------------------------- PROMPT 11= K+79
SUBROUTINE PROMPT (STRING) READINCMD, 1000, ERR-100) (LINE (KK) , KK=I, II)
C— SUB:PROMPT - INLINE PROMPT 1000 FORMAT (80A1)
GO TO 14
INCLUDE 'IOCOM.INC' ENDIF
12 CONTINUE
CHARACTER STRING*!*)
C ----- CHECK FOR COMMENT
WRITE (NTW, ' (A) ' ) STRING
RETURN 14 IF(LINEU).EO.'*') THEN
END IF (MODE. EQ. 'BATCH') CALL FREEHR (NTW)
CALL FREEWR(NOT)
c --------------------------------------------------------------- PROMH GO TO 10
SUBROUTINE PROMH (N) ENDIF
C— SUB: PROMH - "HOLLERITH PROMPT" C
C-2.0 DETERMINE LENGTH-OF-INFORMATION
COMMON MTOT.NP.IAU) C
JJ - LENTRM(LLINE)
INCLUDE 'IOCOM.INC' C
CHARACTER*! NCMND, M C-3.0 DETERMINE LENGTH-OF-DATA AND CONVERT DATA TO UPPER CASE
COMMON /CMND/ NCMND(6) , M(4, 7) C
ISP - ICHARC ')
C ----- PROMPT FOR ARRAY NAMES IA - ICHAR('a')
IF (MODE. EQ.' BATCH') GO TO 900 DO 30 I-l.JJ
DO 200 1=1, N IF(LINE(I).EQ.'
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
SUBROUTINE FREEWR(LUN)
C—SUBlFREEWR - WRITE COMMAND/DATA LINE TO FILE LUN
C LUN - LOGICAL UNIT NUMBER TO WRITE TO
INCLUDE
INCLUDE
'lOCOM.INC'
FRECOM.INC'
WRITE(LUN.2000) (LINE(I),1=1,JJ)
2000 FORMAT (IX,BOA!)
RETURN
END
SUBROUTINE FREEFN(SEP,NC,FOUND)
C—SUB:FREEFN - FINDS NEXT NC-CHARACTER SEPARATOR IN INPUT FILE
C SEP (NO M - CHARACTER STRING
INCLUDE
INCLUDE
'lOCOM.INC'
'FRECOM.INC'
CHARACTER*! SEP(NC)
LOGICAL FOUND
FOUND = .FALSE.
50 CALL FREE
IF(NC.LE.II) THEN
DO 60 N-l.NC
CO IF(SEP(N).NE.LINE(N» GO TO 50
FOUND = .TRUE.
RETURN
ELSE
GO TO SO
ENDIF
Y=0
IS-1
xx=o.o
IF(LINE(I+1).EQ.'-') THEN
IS—1
1=1+1
ELSEIF 'LELE' : (RMIN <• RVALUE (N) <= RMAX) IS O.K.
OPT - 'LTLE1 : (VMIN < RVALUE(N) <= RMAX) IS O.K.
OPT - 'LELT' : (RMIN <= RVALUE(N) < RMAX) IS O.K.
OPT = 'LTLT' : (RMIN < RVALUE(N) < RMAX) IS O.K.
STRING - SINGULAR NOUN DESCRIBING RVALUE
INCLUDE 'IOCOM.INC1
CHARACTER STRING*(•), OPT«4
REAL'S RVALUE(NUH).RMAX.RMIN.RVAL
LOGICAL ERR
DO 500 N-l.NUM
RVAL - RVALUE(N)
IF(OPT.EQ.'LELE') THEN
IF(.NOT.((RMIN.LE.RVAL).AND.(RVAL.LE.RMAX))) THEN
ERR = .TRUE.
WRITE(NTH.2000) STRING,RMIN,' <= value <= ',RMAX, RVAL
WRITE(NOT.*)
WRITE(NOT. 2000) STRING. RMIN.' <= value <- '.RMAX,RVAL
RETURN
ENDIF
FORMAT!
»' •*** ERROR: The value of •, IX, A, IX,'la limited Co the range:',/,
•H2X.G11.4.A.G11.4./.
•f' The given or generated value ls:f,G11.4)
ELSEIF(OPT.EQ.'LELT'I THEN
IF(.NOT.((RMIN.LE.RVAL).AND.(RVAL.LT.RMAX))) THEN
ERR - .TRUE.
WRITE(NTW,2000) STRING,RMIN,' <= value < '.RMAX,RVAL
WRITE(NOT.*)
WRITE(NOT,2000) STRING,RMIN,' <= value < ',RMAX,RVAL
RETURN
ENDIF
ELSEIF(OPT.EQ.'LTLE') THEN
IF(.NOT.((RMIN.LT.RVAL).AND.(RVAL.LE.RMAX))) THEN
ERR = .TRUE.
WRITE (NTH. 2000) STRING. RMIN, ' < value <= ', RMAX. RVAL
WRITE (NOT. •)
WRITE (NOT. 2000) STRING, RMIN, ' < value <- ', RMAX. RVAL
RETURN
ENDIF
ELSEIFIOPT.EQ.'LTLT') THEN
IFI.NOT.((RMIN.LT.RVAL).AND.(RVAL.LT.RMAX))) THEN
ERR = .TRUE.
WRITE(NTW,2000) STRING,RMIN,' < value < '.RMAX,RVAL
WRITE(NOT,*)
WRITE(NOT,2000) STRING,RMIN,' < value < '.RHAX,RVAL
RETURN
ENDIF
ELSE
ERR - .TRUE.
WRITE(NTW,2040)
WRITE(NOT. •)
WRITE(NOT,2040)
RETURN
ENDIF
2040 FORMAT!' *•** ERROR: Call to CKRRNG passed an undefined option.')
500 CONTINUE
RETURN
END
C————————————— ———————————————————————————————— ——-CKIRNG
SUBROUTINE CKIRNG(STRING,IVALUE,NUM,IMIN.OPT.IMAX,ERR)
C—SUB:CKIRNG - CHECKS INTEGER VALUE RANGE
C RETURN ERR-.TRUE. IF NOT O.K.
C VALUE IS A VECTOR OF DIMENSION IVALUE (NUM)
C
C OPT - 'LELE' : (IMIN <- IVALUE(N) <- IMAX) IS O.K.
C OPT - 'LTLE' : (IMIN < IVALUE(N) <- IMAX) IS O.K.
C OPT - 'LELT' : (IMIN <- IVALUE(N) < IMAX) IS O.K.
C OPT " 'LTLT' : (IMIN < IVALUE(N) < IMAX) IS O.K.
C
C STRING - SINGULAR NOUN DESCRIBING IVALUE
INCLUDE 'IOCOM.INC'
CHARACTER STRING*!*), OPT*4
INTEGER IVALUE(NUM),IMAX,IMIN, IVAL
LOGICAL ERR
DO SOO N-l.NUM
IVAL - IVALUE (N)
IF(OPT.EQ.'LELE1) THEN
IF(.NOT.((IMIN.LE.IVAL).AND.(IVAL.LE.IMAX))) THEN
ERR - .TRUE.
WRITE(NTW,2000) STRING,IMIN,' <- value «• ',IMAX,IVAL
WRITE(NOT.*)
WRITE (NOT, 2000) STRING. IMIN,' «• value <= MHAX.IVAL
RETURN
ENDIF
2000 FORMAT!
+' •••* ERROR: The value of'.IX.A,IX,'Is limited to the range:',/,
+12X,I6,A.16,/,
+ ' The given or generated value is:',16)
ELSEIFIOPT.EQ.'LELT') THEN
IF(.NOT.((IMIN.LE.IVAL).AND.(IVAL.LT.IMAX))) THEN
ERR " .TRUE.
WRITE(NTW,2000) STRING,IMIN,' <- value < ',IMAX,IVAL
WRITE(NOT,*)
WRITE(NOT, 2000) STRING,IMIN,' <= value < ',IMAX,IVAL
RETURN
ENDIF
ELSEIF(OPT.EQ.'LTLE') THEN
IF(.NOT.((IMIN.LT.IVAL).AND.(IVAL.LE.IMAX))) THEN
ERR « .TRUE.
WRITE(NTW,2000) STRING,IMIN,' < value <= '.IMAX,IVAL
WRITE(NOT,*)
WRITE(NOT,2000) STRING,IMIN,' < value <= ',IMAX,IVAL
RETURN
ENDIF
ELSEIF(OPT.EQ.'LTLT') THEN
IFI.NOT.((IMIN.LT.IVAL).AND.(IVAL.LT.IMAX))) THEN
ERR = .TRUE.
WRITE(NTW,2000) STRING,IMIN,' < value < '.IMAX.IVAL
WRITE(NOT.*)
WRITE(NOT,2000) STRING,IMIN,' < value < ',IMAX,IVAL
RETURN
ENDIF
ELSE
ERR « .TRUE.
WRITE(NTH,2010)
WRITE(NOT,*)
WRITE(NOT,2010)
RETURN
ENDIF
2010 FORMATC **** ERROR: Call to CKIRNC passed an undefined option.')
SOO CONTINUE
RETURN
END
C CKRZER
SUBROUTINE CKRZER(STRING, RVALUE.NUM.OPT.ERR)
C—SUB:CKRZER - CHECKS REAL VALUE RELATIVE TO ZERO
C RETURNS ERR-.TRUE.IF NOT O.K.
C WHERE VALUE IS A VECTOR OF DIMENSION RVALUE(NUM)
C
OPT - 'LT' : RVALUE(N) .LT. 0.ODD IS O.K.
OPT = 'LE' : RVALUE(N) .LE. O.ODO IS O.K.
OPT - 'GE' : RVALUE (N) .GE. O.ODO IS O.K..
OPT = 'GT' : RVALUE(N) .GT. O.ODO IS O.K.
OPT - 'NE' : RVALUE(N) .NE. O.ODO IS O.K.
Append -24
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
STRING - SINGULAR NOUN DESCRIBING RVALUE
INCLUDE 'IOCOM.INC'
CHARACTER STRING*!*). OPI-2
REAL*8 RVALUE(NUM).RVAL
LOGICAL ERR
DO SOO N-l.NUH
RVAL - RVALUE(N)
IFIOPT.EQ.'LT') THEN
IF(.NOT.(RVAL.LT.O.ODO)) THEN
ERR - .TRUE.
WRITE(NTH,2000) STRING,'must be < O.'.RVAL
WRITE (NOT, •)
WRITE (NOT, 2000) STRING, 'must be < O.'.RVAL
RETURN
ENDIF
2000 FORMAT!' ""•* ERROR: the value of',1X,A.1X.A,/.
+ ' The Qlven or generated value la:',Gil.4)
ELSEIF(OPT.EQ.'LE') THEN
IF(.NOT.(RVAL.LE.O.ODO)) THEN
ERR * .TRUE.
HRITE(NTH,2000) STRING,'muat be <° O.'.RVAL
WRITE (NOT. *)
WRITE(NOT.2000) STRING,'must be <- O.'.RVAL
RETURN
ENDIF
ELSEIFIOPT.EQ.'GE') THEN
IF(.NOT.(RVAL.GE.O.ODO)) THEN
ERR = .TRUE.
WRITE(NTW,2000) STRING,'must be >- O.'.RVAL
WRITE(NOT,•)
NRITE(NOT,2000) STRING, 'mat be >- O.'.RVAL
RETURN
ENDIF
ELSEIF(OPT.EQ.'GT') THEN
IF(.NOT.(RVAL.GT.0.000)1 THEN
ERR = .TRUE.
WRITE (NTW, 2000) STRING,'muat be > O.'.RVAL
WRITE (NOT, ")
WRITE(NOT.2000) STRING,'muat be > O.'.RVAL
RETURN
ENDIF
ELSEIFIOPT.EQ.'NE') THEN
IF(RVAL.EQ.O.ODO) THEN
ERR - .TRUE.
WRITE(NTW,2000) STRING,'must not - O.'.RVAL
WRITE(NOT. •)
WRITE(NOT.2000) STRING.'must not " O.'.RVAL
RETURN
ENDIF
ELSE
ERR - .TRUE.
WRITE(NTH.2010)
WRITE(NOT.«)
WRITE(NOT.2010)
RETURN
ENDIF
ELSEIF(OPT.EQ.'LE') THEN
IFI.NOT. (IVAL.LE.O)) THEN
ERR •' .TRUE.
WRITE, (NTW. 2000) STRING,'muat be <= O.'.IVAL
WRITE. (NOT, •)
WRITE. (NOT, 2000) STRING.'muat be <= O.'.IVAL
RETURN
ENDIF
ELSEIF(OCT.EO.'GE') THEN
IFI.NOT.(IVAL.GE.O)) THEN
ERR - .TRUE.
WRITE; (NTW, 2000) STRING,'must be >•> O.'.IVAL
WRITE(NOT,*)
WRITE(NOT,2000) STRING,'muat be >= O.'.TVAL
RETURN
ENDIF
ELSEIF(OFT.EQ.'GT') THEN
IF(.NOT.(IVAL.GT.0)) THEN
ERR •< .TRUE.
WRITE:(NTW,2000) STRING,'must be > O.'.IVAL
WRITE; (NOT, «»
WRITi: (NOT. 2000) STRING,'muat be > O.'.IVAL
RETURN
ENDIF
ELSEIF(OFT.EQ.'NE') THEN
IF(IVAL.EQ.O) THEN
ERR " .TRUE.
WRITE;(NTW,2000) STRING,'muat not = O.'.IVAL
WRITE(NOT.*)
WRITE;(NOT,2000) STRING,'muat not - O.'.IVAL
RETUP.N
ENDIF
ELSE
ERR - .TRUE.
WRITE(NTH, 2010)
WRITE(MOT,«)
WRITE(MOT,2010)
RETURN
ENDIF
2010 FORMAT!' «**« ERROR: Call to CKIZER paaaed an undefined option.')
500 CONTINUE
RETURN
END
2010 FORMAT(' •*** ERROR: Call to CKRZER paaaed an undefined option.') 2001 FORMAT!1 ••*• '. (A))
SUBROUTINE ALERT(LINE1,LINE2,LINE3)
C—SUB: ERRMSG - WRITES ALERT MESSAGE TO TERMINAL AND OUTPUT FILE
C LINE1 IS ALWAYS WRITTEN'
C LINE2 IS WRITTEN IF LINE2(1:1).NE.'S'
C LINE3 IS WRITTEN IF LINE3(1:1).NE.'S'
INCLUDE 'IOCOM.INC'
CHARACTER LINE1*!*), LINE2*<•). LINE3* («)
WRITE(NTW,2001) LINE1
WRITE(NOT,*)
WRITE(NOT,2001) LINE1
500 CONTINUE
RETURN
END
SUBROUTINE CXIZER(STRING,IVALUE.HUM,OPT,ERR)
C—SUB:CKIZER - CHECKS INTEGER VALUE RELATIVE TO ZERO
C RETURNS ERR".TRUE.IF NOT O.K.
C WHERE VALUE IS A VECTOR OF DIMENSION IVALUE(HUM)
C
IF(LINE2(1:1).NE.'S') THEN
WRITE(NTW.2002) LINE2
WRITE(HOT.2002) LINE2
ENDIF
2002 FORMAT (UX. (A))
IFILINE3(1:1).NE.'S') THEN
WRITE(UTW,2003) LINE3
WRITE(NOT. 2003) LINE3
ENDIF
2003 FORMAT(13X,(A))
C
C
C
C
C
OPT
OPT
OPT
OPT
S
B
a
-
•LT'
•LE'
'GE'
•GT'
: IVALUE (N)
: IVALUE (N)
: IVALUE (N)
: IVALUE (N)
SINGULAR NOUN
.LT.
.LE.
.GE.
.GT.
0
0
0
0
IS 0
IS O
IS 0
IS O
.X.
.K.
.X.
.K.
C 3
RETURN
END
. 0 DYNAMIC &RRAY MANAGEMENT ROUTINES
INCLUDE 'IOCOM.INC'
CHARACTER STRING*!*), OPT*2
INTEGER IVALUE(NUM),IVAL
LOGICAL ERR
DO 500 N=1.NUM
IVAL - IVALUE(N)
IFIOPT.EQ.'LT') THEN
IFI.NOT.(IVAL.LT.O)) THEN
ERR = .TRUE.
WRITEINTW.2000) STRING,'muat be < O.'.IVAL
WRITE(NOT,«)
WRITE(NOT.2000) STRING,'muat be < O.'.IVAL
RETURN
ENDIF
2000 FORMAT)' **** ERROR: The value of',IX,A,IX,A,/,
+' The given or generated value Is:',16)
SUBROUTINE DEFINR(NAME, NA.NR.NC)
C—SUB:DEFINR - DEFINE DIRECTORY AND RESERVE STORAGE
C FOR REAL ARRAY IN DATABASE
C NAME - NAME OF ARRAY
C NA - BLANK COMMON POINTER TO ARRAY (RETURNED)
C NR NUMBER OF ROWS
C NC NUMBER OF COLUMNS
COMMON MTOT,NP.IA(1)
CHARACTER*! NAME(4)
NP - 2
CALL DEFIN(NAME,NA.NR.NC)
RETURN
END
SUBROUTINE DEFINI(NAME.NA,NR.NC)
C—SUB:DEFINI - DEFINE DIRECTORY AND RESERVE STORAGE
C FOR INTEGER ARRAY IN DATABASE
C NAME = NAME OF ARRAY
C NA o BLANK COMMON POINTER TO ARRAY (RETURNED)
Append -25
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
NR - NUMBER OF ROWS
NC - NUMBER OF COLUMNS
COMMON MTOT.NP.IAd)
CHARACTER*! NAME(4)
HP •> 1
CALL DEFIN(NAME,NA,NR,NC)
RETURN
END
SUBROUTINE DEFIN(NAME, NA. NR, NC)
c DEFINE AND RESERVE STORAGE FOR ARRAY
100 CALL ICON(NAME.IA(I))
IAU + 4) NR
IA (1+51 NC
IA (1+61 NP
IAd+7) ISTR
IA(I + 8I 0
IA(I+9) 0
900 RETURN
2000 FORMAT!
•• .... ERROR: Insufficient blank COMMON storage.',/,
•• Storage required MTOT -Ml./,
*' Storage available MTOT =',I7)
END
COMMON MTOT.NP.IAd)
INCLUDE 'ARYCOM.INC'
INCLUDE 'IOCOM.INC'
CHARACTER*1 NAME(4)
c DEFIN VARIABLES
C
C
C
c
c
c
c
c
c
c
c
c
c
c
c-
c
c
c
c
c
c
c
c
c
NC '
MTOT
NUMA
NEXT
IDIR
IP
LENR
NAME - NAME OF ARRAY - 4 LOGICALS MAXIMUM
NA - LOCATION OF ARRAY IF IN BLANK COMMON
NR - NUMBER OF RONS
' NUMBER OF COLUMNS
END OF DIRECTORY
NUMBER OF ARRAYS IN DATA BASE
NEXT AVAILABLE STORAGE LOCATION
START OF DIRECTORY IN BLANK COMMON
NUMBER OF LOGICALS CONTAINED IN DATA TYPE
NUMBER OF LOGICALS IN PHYSICAL RECORD
NP - TYPE OF DATA
- 1 INTEGER DATA
- 2 REAL DATA
- 3 LOGICAL DATA
-DIRECTORY DEFINITION FOR CORE OR SEQUENTIAL FILES
IDIR(l.N) - NAME OF ARRAY - INAME (4 CHAR.)
IDIR(S.N) - NUMBER OF RONS - NR
IDIRI6.N) - NUMBER OF COLUMNS - NC
IDIR(7,N) - TYPE OF DATA - NP
IDIRI8.N) - INCORE ADDRESS - NA
- -1 IF SEQUENTIAL FILE ON DISK
- -2 IF DIRECT ACCESS ON DISK
SIZE OF ARRAY
- 0 IF IN CORE STORAGE
SUBROUTINE LOCATE (NAME, NA.NR.NC)
C—SUB:LOCATE - LOCATE ARRAY "NAME" AND RETURN
C NA - POINTER TO LOCATION IN BLANK COMMON
C NR = NUMBER OF RONS
C NC - NUMBER OF COLUMNS
COMMON MTOT.NP.IAd)
IDIRI9.N) •
IDIR(IO.N)
CHARACTER*! NAME
DIMENSION NAME(4),INAME<4)
c LOCATE AND RETURN PROPERTIES ON ARRAY
NA - 0
CALL ICON(NAME.INAME)
I - IFIND(INAME,0)
IF(I.EQ.O) GO TO 900
C RETURN ARRAY PROPERTIES
NA - IA(I+7)
NR - IAd + 4)
NC - IAd+5)
NP - IAd+6)
900 RETURN
END
C DIRECTORY DEFINITION FOR DIRECT ACCESS FILES
C IDIR(S.N) - NUMBER OF INTEGERS
C IDIR(6,N) = NUMBER OF REAL WORDS
C IDIR(7,N) - NUMBER OF LOGICALS
C IDIRie.N) - NUMBER OF LOGICAL RECORDS
C IDIR(9,N) - LOGICAL RECORD NUMBER
C IDIR(IO.N) - LUN IF ON LOGICAL UNIT LUN
c EVALUATE STORAGE REQUIREMENTS
NSIZE = (NR*NC*IP(NP) -1)/(IP(1)*2)
NSIZE - NSI2E*2 + 2
NA - NEXT
NEXT - NEXT + NSIZE
C SET UP NEW DIRECTORY
NUMA = NUMA + 1
IDIR - IDIR - 10
I » IDIR
C CHECK STORAGE LIMITS
IF(I.GE.NEXT) GO TO 100
I - NEXT - I + MTOT - 1
WRITE(NTW.2000) I,MTOT
WRITE(NOT,*)
WRITE(NOT,2000) I,MTOT
PAUSE
STOP
100 CALL ICON (NAME. IA(I) )
LAd+O NR
IAd+5) NC
IAd+6) NP
IAd+1) NA
IAd+8) NSIZE
IAd+9) 0
900 RETURN
2000 FORMAT<
•' •*•* ERROR: Insufficient blank COMMON storage.'./,
«' Storage required MTOT -',I7,/,
•' Storage available MTOT -M7)
END
SUBROUTINE DEFDIR(NAME,NR.NC,ISTR)
C—SUB:OEFDIR - DEFINE DIRECTORY FOR OUT-OF-CORE FILE
C NAME - NAME OF ARRAY
C NR ' NUMBER OF RONS
C NC = NUMBER OF COLUMNS
C ISTR - OUT OF CORE FLAG (=-1)
SUBROUTINE DELETE(NAME)
C—SUB:DELETE - DELETE ARRAY "NAME- FROM DATABASE
COMMON MTOT.NP.IAd)
INCLUDE 'ARYCOM.INC'
INCLUDE 'IOCOM.INC'
CHARACTER*! NAME
DIMENSION NAME (4), INAME (4)
C DELETE ARRAY FROM STORAGE
100 CALL ICON(NAME,INAME)
I - IFIND(INAME,0)
IF(I.EQ.O) GO TO 900
C CHECK ON STORAGE LOCATION
200 NSIZE - IAII+8)
C SET SIZE OF ARRAY
' NEXT = NEXT - NSIZE
NUMA - NUMA - 1
NA • IAd+7)
C CHECK IF OUT OF CORE OR DIRECT ACCESS
IF(NA.GT.O) GO TO 500
WRITE(NTH,1000) NAME
WRITE(NOT,•)
WRITE(NOT.1000) NAME
GO TO SOO
500 IF(NA.EQ.NEXT) GO TO 800
C COMPACT STORAGE
II - NA + NSIZE
NNXT = NEXT - 1
DO 700 J=NA,NNXT
IA(J) ° IAIII)
700 II - II + 1
C COMPACT AND UPDATE DIRECTORY
600 NA - I - IDIR
IDIR - IDIR + 10
IF(NA.EQ.O) GO TO 900
NA - NA/10
DO 860 K-l.NA
II - I + 9
DO 850 J=l,10
IA(II) = IAdI-10)
850 II - II - 1
IF(IA(I»7).LE.O) GO TO 860
IF(IA(I+9).EQ.O) IA(I+7) - IA(I+7) - NSIZE
860 I » I - 10
C
900 RETURN
1000 FORMAT!' — Name '.4A1.' la being used for
* ' OUT-OF-CORE file.1,/)
END
COMMON MTOT.NP.IAd)
INCLUDE 'ARYCOM.INC'
INCLUDE 'IOCOM.INC'
CHARACTER*1 NAME(4)
C EVALUATE STORAGE REQUIREMENTS
IF(NP.EQ.O) NP - 2
C SET UP NEW DIRECTORY
NUMA - NUMA + 1
IDIR - IDIR - 10
I - IDIR
C CHECK STORAGE LIMITS
IF(I.GE.NEXT) GO TO 100
I - NEXT - I » MTOT - 1
WRITE(NTW.2000) I,MTOT
WRITE (NOT. • I
WRITE (NOT. 2000) I, MTOT
PAUSE
STOP
SUBROUTINE ICON (NAME. INAME)
CHARACTER*! NAME(4)
DIMENSION INAME(4)
C CONVERT LOGICALS TO INTEGER DATA
DO 100 1=1.4
100 INAME (II - ICHAJU NAME(I) )
C
RETURN
END
FUNCTION IFIND(INAME.LUN)
C—FUN: IFIND - FIND
COMMON MTOT.NP.IAd)
INCLUDE 'ARYCOM.INC'
DIMENSION INAME(4)
c FIND ARRAY LOCATION
I m IDIR
Append -26
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
DO 100 N-l.NUMA
IF(LUN.NE.IA(I+9)) GO TO 100
IF (INAMEU) .NE.IAII )) GO TO 100
IF (INAHEI2).NE.IAII*!)) GO TO 100
IF (INAMEO).NE.IA(H-2)) GO TO 100
IF 0 : generate right, [Tl. transformations
Eigenvectors of real eigenvalues occurr as rows (cols) of [T|-l
((T)). Eigenvectors for a complex eigenvalue pair aj.l ± laj,j+1
may be formed by tj ± ltj+1 where t}, tj+1 are the corresponding
rows (cols) of [T|-l (IT|)
Iterations are limited to 50 maximum. On exit from the procedure
TMX records the number of iterations performed. Failure to
converge is indicated by TMX«50 or. if all transformations in
one Iteration are the identity matrix, by TMX<0.
The machine dependent variable EP is set to IE-OS and should be
reset for machine precision available.
: DICTIONARY OF VARIABLES
VARIABLI:-
--INPUT
A(N.N)
N
TMS
—OUTPUT
T(N,N)
TMX
—LOCAL
DESCRIPTION
Array to be analyzed.
System size
Control parameter
Array to receive eigenvectors.
Iteration count/iteration flag
EP
Precision
IMPLICIT REAL*8(A-H,O-Z)
REAL-8 MN.N) ,T(N,N) ,EP
INTEGER N.TMX
LOGICAL MARK, LEFT, RIGHT
C
C—0.0 INITIALIZE CONTROL VARIABLES
C
IFIEP.LI1.O.ODO) EP - l.OD-9
EPS - SQRT(EP)
LEFT - .FALSE.
Append -27
-------
Indoor Air Quality Model
Phase III Report
APPENDIX
CONTAM87 FORTRAN 77 Source Code
RIGHT - .FALSE.
IF(TMX.LT.O) THEN
LEFT - .TRUE.
ELSEIF(TMX.OT.O) THEN
RIGHT - .TRUE.
ENDIF
MARK - .FALSE.
C
C--1.0 INITIALIZE [T| AS IDENTITY MATRIX
C
IF(THX.NE.O) THEN
DO 10 I-l.N
T(I,I) - l.ODO
DO 10 J-I+l.N
TII.J) • O.ODO
TIJ.I) - O.ODO
10 CONTINUE
ENDIF
C
C—2.0 MAIN LOOP
C
C-MAC WRITE!*,MSX.A.M') ' '
DO 26 IT-1,50
•MAC WRITE!*,'(A,\)') '+'
•-2.1 IF MARK IS SET
TRANSFORMATIONS OF PREVIOUS ITERATION WERE OMITTED
PROCEDURE WILL NOT CONVERGE
IF(MARK) THEN
TMX - 1-IT
RETURN
ENDIF
C
C—2.2 COMPUTE CONVERGENCE CRITERIA
C
DO 20 I-l.N-1
All - All,I)
DO 20 J-l+l.N
AIJ - A(I.J)
AJI - AIJ,I)
IF((ABS(AIJ+AJI).GT.EPS) .OR.
+ ((ABS(AIJ-AJI). GT.EPS). AND. (ABS (AII-AU, J)) . GT.EPS))) THEN
GOTO 21
ENDIF
20 CONTINUE
TMX - IT -1
RETURN
C
C—2.3 BEGIN NEXT TRANSFORMATION
C
21 MARK = .TRUE.
• DO 25 K-l.N-1
DO 25 M-K-H.N
B > O.ODO
G = O.ODO
HJ - O.ODO
YH - 0.000
DO 22 I-l.N
AIK - A(I,K)
AIM = Ad,HI
TE - AIK*AIK
TEE - AIM-AIM
YH = YH + TE - TEE
IFIII.NE.K) .AND. (I.NE.MM THEN
CHY - l.ODO
SHY • O.ODO
ELSE
CHY » 1.0DO/SQRTI1.0DO - TANHY*TANHY)
SHY - CHY*TANHY
ENDIF
C
C
C
C
PUT
Cl - CHY*CX - SHY*SX
C2 » CHY*CX + SHY*SX
SI • CHY*SX + SHY*CX
32 - -CHY«SX + SHY*CX
— — APPLY TRANSFORMATION IF WARRANTED
IF((ABS(SI).GT.EP).OR.(ABSIS2).GT.EP)) THEN
MARK - .FALSE.
TRANSFORMATION ON THE LEFT
DO 23 I-l.N
AKI - AIK, I)
AMI - AIM. I)
A(K.I) - Cl'AKI + SI-AMI
AIM,I) » 32•AKI + C2*AMI
IF(LEFT) THEN
TKI - T(K,I)
TMI - KM, I)
T(K,I) • Cl'TKI + S1*TMI
T(M,I) - S2*TKI + C2*TMI
ENDIF
23 CONTINUE
TRANSFORMATION ON THE RIGHT
DO 24 I-l.N
AIK - A(I,K)
AIM - A(I.M)
A(I,K) • C2*AIK - S2»AIM
Ad.MI i -Sl'AIK * Cl'AIM
IF(RIGHT) THEN
TIK - TII.K)
TIM - T(I,M)
TII.K) - C2«TIK - S2*TIM
TII.M) - -Sl'TIK + C1«TIM
ENDIF
24 CONTINUE
ENDIF
25 CONTINUE
26 CONTINUE
TMX - 50
RETURN
END
AKI
AMI
H
TEP
TEN
G
HJ
ENDIF
CONTINUE
H
D
AXM
AMK
C
E
AIK, I)
AIM. I)
H + AKI*AMI - AIK'AIM
TE + AMI"AMI
TEE * AKI'AKI
G » TEP + TEM
BJ - TEP + TEM
H
A(K.K) - A(M.M)
A(K.M)
A(M.K)
AKM + AMK
AKM - AMK
COMPUTE ELEMENTS OF (Ri|
IF(ABSIC).LE.EP) THEN
CX = l.ODO
SX = O.ODO
ELSE
COT2X • D/C
SIG - SIGN(1.ODO.COT2X)
COTX - COT2X + (SIG'SQRT(l.ODO + COT2X«COT2X)I
SX = SIG/SQRT(l.ODO + COTX'COTX)
cx = sx-corx
ENDIF
IFIYH.LT.O.ODO) THEN
TEM m CX
CX = SX
SX - -TEM
ENDIF
COS2X - CX*CX - SX*SX
SIN2X « 2.0DO*SX«CX
D » D*COS2X + C*SIN2X
H - H*COS2X - HJ'SIN2X
DEN " G * 2.0DO*IE*E + D'DI
TANHY = (E*D - H/2.0DOI/DEN
COMPUTE ELEMENTS OF [SI]
IF IABS(TANHY).LE.EP) THEN
Append -28
-------
Indoor Air Quality Model APPENDIX
Phase III Report CONTAM87 FORTRAN 77 Source Code
Include Files
C CALSAPX FR HE-FIELD INPUT SFRECOM.INC
C CALSAPX ARRAY MANAGEMENT $ARYCOM.INC CHARACTER LINE'l. LLINE*1«0
c
COMMON /CLINE1/ LINE (160)
COMMON /ARYCOM/ NUMA.NEXT. IDIR. IP <3) COMMON /CLINE2/ II. JJ
C MTOT SIZE OF BLANK COMMON VECTOR IA
C NP CURRENT DATA TYPE: 1-INTEGER: 2-REAL; 3-CHAR. c ____ VARIABLE _________ DESCRIPTION
C IA
-------
NBS.114A (REV. 2-aci
U.S. DEPT. OF COMM.
BIBLIOGRAPHIC DATA
SHEET (See instructions)
I. PUBLICATION OR
REPORT NO.
NBSIR 88-3814
2. Performing Organ. Report No.
3. Publication Date
JULY 1988
4. TITLE AND SUBTITLE
Progress Toward a General Analytical Method for Predicting Indoor Air Pollution in
Buildings - Indoor Air Quality Modeling Phase III Report
5. AUTHOR(S)
James Axley
6. PERFORMING ORGANIZATION (If joint or other than NBS. see instructions;
NATIONAL BUREAU OF STANDARDS
U.S. DEPARTMENT OF COMMERCE
GATTHERSBURG, MD 20899
7. Contract/Grant No.
8. Type of Report & Period Covered
9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State, ZIP)
National Bureau of Standards
Building Environment Division 747
Center for Building Technology
ftaithersbure. MD 20899
10= SUPPLEMENTARY NOTES
[ | Document describes a computer program; SF-185, FIPS Software Summary, is attached.
11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant
bibliography or literature survey, mention it here) This interim report presents the results of Phaselll
of the NBS General Indoor Air Pollution Concentration Model Project. It describes;
a)a general element-assembly formulation of multi-zone contaminant dispersal analysis
theory that provides a general framework for the development of detailed (element)
models of mass transport phenomena that may affect contaminant dispersal in buildings.
b)an approach to modeling the dispersal of interactive contaminants involving con-
taminant mass transport phenomena governed by basic principals of kinetics and intro-
duces a linear first order kinetics element to achieve this end, c)an approach to
modeling the details of contaminant dispersal driven by convection-diffusion processes
in one-dimensional flow situations (e.g., HVAC ductwork) and introduces a convection-
diffusion flow element to achieve this end, and d) the features and use of CONTAM87, a
program that provides a computational implementation of the theory and methods discussed
The theory and methods presented are based upon a slight generalization of the building
idealization employed earlier (Axley, 1987). Here, building air flow systems are
idealized as assemblages of mass transport elements, rather than simply flow elements a
used previously, connected to discrete system nodes corresponding to well-mixed air
zones within the building and its HVAC system. Equations governing contaminant dis-
persal in the whole building air flow system due to air flow and reaction or sorption
mass transport phenomena are formulated by assembling element equations, that character-
ize a specific instance of mass transport in the building air flow system, in such a
manner that the fundamental requirement of conservation of mass is satisfied in each
ZOTl
12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)
contaminant dispersal analysis, inverse contaminant dispersal analysis, tracer gas
techniques, building simulation
13. AVAILABILITY
[JQ Unlimited
| | For Official Distribution. Do Not Release to NTIS
| | Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.
20402.
[y~] Order From National Technical Information Service (NTIS), Springfield, VA. 22161
14. NO. OF
PRINTED PAGES
124
15. Price
$18.95
USCOMM-DC 0043-P80
------- |