United States

Department of

Agriculture

Forest Service

Forest

Products

Laboratory

Research

Paper

FPL–RP–561

Anton TenWolde

FPL Roof Temperature 

and Moisture Model

Description and Verification

solar radiation infrared radiation

indoor air

ventilation

air

wind
Abstract

This paper describes a mathematical model developed by the

Forest Products Laboratory to predict attic temperatures,

relative humidities, and roof sheathing moisture content.

Comparison of data from model simulation and measured

data provided limited validation of the model and led to the

following conclusions: (1) the model can provide reasonably

accurate estimates for temperatures of roof sheathing and attic

air, although heat storage effects often cause delay of 1 to 2 h

in attic air temperatures; (2) the model can accurately predict

the frequency of occurrence of high roof sheathing temperatures (> 120°F (49°C)) during summer, but accuracy is

highly dependent on solar absorptance and emissivity values

of the roof shingles; (3) the model consistently overpredicts

the extent of night-time cooling from sky radiation losses,

leading to predicted temperatures that are too low; (4) treatment of the effect of snow cover is too simplistic, but no

better alternatives are apparent for simulating this very complex behavior; (5) the model apparently can predict average

moisture conditions in the sheathing with reasonable accuracy, generally within 1% moisture content, when moisture

content is not excessively high or low; and (6) hourly moisture behavior is not represented as well as is daily or seasonal behavior, especially for north-facing sheathing. The

model would benefit from verification with data that include

measured emissivity and solar absorptance of the shingles,

addition of thermal mass in attic and roof, better algorithms

to calculate direct and diffuse solar radiation, and verification

for roof with east–west orientation.

Keywords: attic, fire retardant, humidity, moisture, plywood,

roof, temperature, ventilation.

April 1997

TenWolde, Anton. 1997. FPL roof temperature and moisture model:

description and verification. Res. Pap. FPL–RP–561. Madison, WI: U.S.

Department of Agriculture, Forest Service, Forest Products Laboratory.

48 p.

A limited number of free copies of this publication are available to the

public from the Forest Products Laboratory, One Gifford Pinchot Drive,

Madison, WI 53705–2398. Laboratory publications are sent to more than

1,000 libraries in the United States and elsewhere.

The Forest Products Laboratory is maintained in cooperation with the

University of Wisconsin.

The United States Department of Agriculture (USDA) prohibits discrimination in its programs on the basis of race, color, national origin, sex,

religion, age, disability, political beliefs, and marital or familial status.

Persons with disabilities who require alternative means of communication

of program information (braille, large print, audiotape, etc.) should contact

the USDA Office of Communications at (202) 720–2791. To file a complaint, write the Secretary of Agriculture, U.S. Department of Agriculture,

Washington, DC 20250, or call 1–800–245–6340, or (202) 720–1127

(TTD). USDA is an equal employment opportunity employer.

Acknowledgments

Financial support for model revisions, documentation, and

validation was provided by a grant from the New Jersey

Department of Community Affairs. The author thanks Tom

Gorman of the University of Idaho for his work on initial

development of the model, and William Rose and Jeffrey

Gordon of the Building Research Center, University of Illinois–Champaign, for their assistance with model verification

and transfer of measured data. The author thanks Clifford

King of Wood Construction Technologies, Inc., for his

suggestions for improvements of the model and other helpful

comments. The author also acknowledges Hoover Treated

Wood Products and Wood Construction Technologies, Inc.,

for financial assistance. The initial development of the FPL

roof temperature and moisture model was made possible

through a U.S. Department of Agriculture Competitive

Grant.

Contents

Page

Introduction................................................................1

Description of FPL Model.............................................1

Origin....................................................................1

Capabilities.............................................................2

Structure.................................................................2

Home Configuration .................................................3

Ventilation..............................................................4

Attic Heat Balance....................................................4

Attic Moisture Balance..............................................6

Model Validation ........................................................7

Introduction ............................................................7

Description of Test House..........................................7

Data Collection........................................................7

Model Runs............................................................7

Results...................................................................8

Conclusions .............................................................14

References................................................................15

Nomenclature............................................................16

Appendix A—Convective Heat Transfer Correlations.......17

Appendix B—Program Listing for FPL Model ...............18

Appendix C—SUNDAT Program ................................43
FPL Roof Temperature

and Moisture Model

Description and Verification

Anton TenWolde, Research Physicist

Forest Products Laboratory, Madison, Wisconsin

Introduction

Attic ventilation has historically been considered an effective

strategy for preventing moisture problems in attics. Ventilation usually relies on fixed openings in the end walls, eaves,

and ridge, although an increasing number of attics are

equipped with mechanical fans. Recommendations for minimum attic ventilation rates necessary to prevent condensation

during winter are primarily based on the results from laboratory measurements of small test houses at the University of

Minnesota (Rowley and others 1939) that were insulated

well below current insulation standards.

Recently, concerns have been raised about the effectiveness of

attic ventilation in reducing moisture levels in the attic or

roof (TenWolde and Carll 1992). Depending on climate,

ventilation air may occasionally carry moisture into the attic

rather than remove it. During the night, the roof may be as

much as 10°F to 15°F (5°C to 8°C) colder than the ambient

air as a result of radiation to the sky. This can cause moisture from the ventilation air to be absorbed by, or condense

on, the sheathing, increasing its moisture content. Ventilation also limits the opportunities to insulate the roof, especially a flat roof or cathedral ceiling. When converting an

attic space into heated living space, the need for ventilation

between the sheathing and insulation has been called into

question.

The Forest Products Laboratory (FPL) roof temperature and

moisture model was created to investigate the effect of construction, ventilation strategies, and indoor and outdoor

climate on moisture conditions in the attic (Gorman 1987).

Although the model was conceived to analyze winter conditions, it was soon used to predict temperatures and moisture

content in fire-retardant-treated (FRT) plywood roof sheathing during the summer. This work was precipitated by the

failure of several FRT plywood roof decks, presumably as a

result, in part, of thermal degradation (LeVan and Collet

1989).

The purposes of the model were (1) to provide an estimate of

temperature and moisture content of roof sheathing and

temperature and humidity of attic air at any geographical

location for which hourly weather data are available and

(2) to evaluate the effect of design and indoor climate. The

intent was to achieve this goal in the simplest way possible

with the greatest practical benefits. More detail in heat and

moisture transfer equations does not necessarily provide

better information because of uncertainties in input data and

parameters, such as actual attic ventilation rates and other air

flows.

Since the conception of the FPL attic model, Burch (1992)

has adapted the MOIST computer model, originally developed for moisture analysis of walls, to examine moisture in

roofs. The MOIST program contains much more detail on

moisture transfer in building materials than does the FPL

attic model, and it is capable of calculating liquid moisture

flows. However, the MOIST program does not currently

include interaction between different sheathing surfaces, and it

contains less detailed convective heat transfer coefficient

calculations for the exterior roof surface. Both models make

many simplifying assumptions that may affect their accuracy.

Description of FPL Model

Origin

Burch and Luna (1980) developed the first mathematical

model that simulates moisture transfer to and from the attic.

Burch and others (1984) modified this model to account for

moisture sorption in the sheathing. Cleary (1984) was able

to predict daily variations in attic humidity by assuming

moisture equilibrium between the attic air and the sheathing.

This means that the moisture content and temperature of the

sheathing determine the absolute humidity of the attic air.

Cleary observed that the absolute humidity of the attic air

peaks at mid-day because solar radiation drives moisture

from the sheathing.

The FPL roof temperature and moisture model is based on

the models developed by Burch and Luna (1980) and Cleary

(1984). By adding three heat and moisture balance equations

to Burch and Luna’s model, we are able to differentiate
2

between two separate surfaces of a gable roof and account for

differences in incident solar radiation (Fig. 1). Where Cleary

treated the roof sheathing as one single element, we distinguish between moisture transfer between the attic air and a

thin surface layer, and moisture transfer between this surface

layer and the interior of the sheathing. The interior of the

sheathing is subdivided into two layers. As in the previous

models, we ignored moisture storage in other wood framing

members and other hygroscopic materials in the attic.

An improved subroutine to calculate surface heat transfer

coefficients was also added to the heat balance. We incorporated a subroutine developed by Wilkes (1989) in the existing model. The heat transfer coefficients are a function of

temperature, wind speed, and direction of heat flow. The

formulas used in this subroutine were taken from Holman

(1981) (App. A). We also added the ability to evaluate the

effect of mechanical attic ventilation, running continuously or

controlled by timer, humidistat, or photosensor.

Capabilities

The data required for the FPL model are described in

Table 1. The FPL model (App. B) runs on a personal

computer, programmed in FORTRAN, and it calculates the

following:

• Hourly temperatures in sheathing (two sides of gable roof)

• Hourly average moisture content of sheathing

(two sides of gable roof)

• Hourly moisture content of inside sheathing surface

(two sides of gable roof)

• Hourly attic air temperature

• Hourly attic air relative humidity

• Fraction of time attic fan is turned on

Structure

The program sequentially solves a heat balance and moisture

balance for the attic space (Fig. 1), which are calculated every

10 min, using interpolated hourly weather data (snow cover

data are not interpolated) and interpolated solar radiation

data. Temperatures and attic air relative humidity (RH)

during the first 10-min period and moisture content results

from the last 10-min period of the previous hour are reported

as hourly results.

Figure 2 shows the interaction and function of the main

program components. SUNDAT calculates hourly solar

Figure 1—Schematic diagram of attic.

Table 1—Data required for FPL roof temperature and

moisture modela

Hourly

 weather data

Outdoor temperature (°F)

Outdoor dew point (°F)

Wind speed (knots)

Solar radiation (Btu/h⋅ft2

)

b

Cloud cover (0–10)

Snow cover (0 or 1)

Building data Length of building (ft)

Width of building (ft)

Distance from floor of first floor to ceiling

 of top floor (ft)

Roof overhang (ft)

Thickness of roof sheathing (in.)

Building orientation (input for

 SUNDAT, degrees from north)

Pitch of roof (rise/run)

R-value of attic floor (h⋅ft2

/Btu)

Combined R-value of roof sheathing

 and roofing materials (h⋅ft2

/Btu)

R-value of roof sheathing (h⋅ft2

/Btu)

R-value of end walls and eaves (h⋅ft2

/Btu)

Solar absorptance of roof shingles

Fan capacity (ft3

/min), type of attic fan

 control,c

 and control set points

Other data Wind shading (0–1)

Ventilation rate of living space (ach)d

Fraction of living space ventilation air

 escaping to attic

Natural attic ventilation rate (ach)

Living space air temperature (°F)

Living space air relative humidity (0–1)

Initial moisture content of sheathing

 (two roof orientations, 0–1)

a

TF = TC (1.8) +32; 1 J = 9.5 × 10−4 Btu; 1 m = 3.28 ft;

 1 mm = 0.039 in.

b

Solar radiation on the two roof sides must first be calculated

 with a separate program (SUNDAT) from solar radiation

 weather data (in Langley units; 1 Langley = 1 cal/cm2

 =

 3.687 Btu/ft2

) and data on home location and roof

 orientation.

c

Humidistat, timer, or photo-sensor.

d

ach = air change per hour.
3

radiation on each roof surface from measured solar radiation

(on horizontal surface, in Langley units), geographical location, roof orientation, and roof pitch, and converts the results

into British thermal units per hour per square foot

(Btu/h ft2

).1 SUNDAT creates two input files for the main

program, one for each roof surface. SUNDAT is described in

detail in Appendix C. Figure 3 shows all programs, subroutines, and functions with their major purpose.

The main program uses hourly solar radiation and weather

data in a steady-state heat balance to calculate temperatures in

the attic and roof sheathing. Because the surface heat transfer

coefficients depend on the temperatures, the heat balance is

solved iteratively. The moisture balance is solved on the

basis of these temperatures and previous moisture content.

The heat balance is solved again, adjusting for calculated

latent heat flows. This process is repeated until latent heat

 

1

1 Btu/h⋅ft2 = 3.155 W/m2.

flows converge within 5%. This iterative procedure is

repeated six times each hour, representing 10-min intervals,

using interpolated weather and solar radiation data.

Home Configuration

Building

The building is assumed to be rectangular with a simple

gable roof. Building length, width, and height (distance from

floor of first floor to ceiling of top floor) are used to determine

the approximate volume of the living space. This volume, in

turn, is used to calculate air flow into the attic from below,

using the air exchange rate and the fraction escaping through

the ceiling to the attic.

Location,

roof orientation

Weather data SUNDAT

Program

Attic moisture

balance

Attic heat

balance

Air exchange

living space

conditions

Attic

ventilation

Temperatures

Moisture

contents

Attic air RH

Attic

temperatures

Latent

heat effects

Solar

radiation

Solar radiation

on roof surfaces

Figure 2—Structure of FPL roof temperature and

moisture model.

SUNDAT

•

AHM

• reads input files

• reads run parameters

• calls subprograms

PARAM

• reads/changes parameters

HEAT

• solves attic heat balance

INFRED

• calculates sky radiation loss

HCON

• calculates heat transfer coefficients

MOIST

• solves attic moisture balance

• calculates latent heat

AHR

• calculates attic air humidity ratio

PSAT

• calculates saturation vapor pressure

WHR

• calculates equilibrium humidity ratio

prepares solar

radiation input files

Figure 3—FORTRAN programs, subroutines, and

functions in FPL model.
4

Attic

The dimensions of the attic are calculated from building

length, building width, roof overhang, and roof pitch. The

model assumes a simple isosceles triangular shape (see

Fig. 1). The R-values of the attic floor, roof (including all

roofing materials), roof sheathing, and end walls are also

required. Solar absorptance of the exterior roof surface is used

in the heat balance to calculate the amount of solar radiation

absorbed.

Ventilation

Attic Ventilation

Air enters the attic directly from the outside and from the

living space below. The natural attic ventilation rate specified by the user sets the amount of outside air entering the

attic. The net amount leaving the attic is the sum of the two

air flows. Thus, even when attic ventilation is specified as

zero, an equal amount of air escapes from the attic as enters

from below. All flows are held constant throughout the

simulation run.

When the user specifies an attic fan, the fan is assumed to

provide all the outside air to the attic while the fan is running. When the fan is off, the specified natural attic ventilation rate is assumed. Air leakage through the ceiling is

assumed to be unaffected by the fan. The fan can be controlled with a humidistat, a timer, or a photosensor. If the

user specifies a timer control, the user is asked to specify the

hours at which the fan is turned on and off. The fan can also

be operated continuously by providing the appropriate timer

control settings. With a humidistat, the user specifies the

minimum attic RH setting and with a photosensor, the

minimum solar radiation level. When snow is on the roof,

the photosensor will not receive any solar radiation and the

attic fan will remain off.

Living Space Ventilation

The user specifies the overall ventilation rate of the living

space in air changes per hour and the fraction of the exhaust

air escaping to the attic. The user needs to judge the relationship between air flow through the ceiling and total ventilation rate by choosing appropriate values. These air flows

are used in the attic heat balance and moisture balance.

Attic Heat Balance

The attic heat balance is solved in the HEAT subroutine.

Six simultaneous equations are used to calculate the temperature at the nodes shown in Figure 1. Heat balance equations

for the attic air (node 1), attic floor (node 2), and underside of

the sheathing (nodes 3a and 3b) are the same as those used

by Burch and Luna (1980). We added heat balance equations

for the exterior roof surface (nodes 4a and 4b) to allow for the

effects of solar radiation and night radiant heat loss. The attic

heat balance includes heat conduction through the roof,

ceiling, end walls, and soffits; air flows from the living

space and attic ventilation; radiation exchange between attic

floor and roof sheathing; solar radiation; radiation loss to the

sky; and latent heat effects.

Solar Radiation

Solar radiation on each roof surface is calculated with the

SUNDAT program. The program calculates the ratio of

radiation on a surface with a specified tilt and orientation and

on a horizontal surface at a given location and hour. The

calculation is performed for a cloudless sky, apportioning

direct and diffuse radiation based on chapter 27 of the

ASHRAE Handbook of Fundamentals (ASHRAE 1993).

This ratio is then applied to measured solar radiation data

(horizontal surface) to obtain solar radiation incident on the

specified tilted surface. SUNDAT produces solar radiation

input files for the main program. In the main program, total

radiation absorbed by the roof surface is obtained by multiplying incident radiation by the roof’s solar absorptance

(see App. C).

Surface Heat Transfer Coefficients

Calculation of surface heat transfer coefficients takes place in

the HCON subroutine (Fig. 3). The formulas for natural and

forced convective heat transfer coefficients are shown in Appendix A. Natural convection is assumed for the floor surface

and underside of the roof, and a combination of natural and

forced convection applies to the exterior roof surface. The

formulas used to calculate natural and forced convection

(App. A) are taken from Holman (1981) and were also used

by Wilkes (1989) in his roof thermal performance model. For

the characteristic length of the plate, which is needed to

calculate the Reynolds, Raleigh, and Nusselt numbers, we

use the average of length and width of the surface (e.g., roof

surface), primarily because the direction of the convective or

forced air flow is unknown. We use wind speed directly for

calculating the forced convective transfer coefficient. Wind

direction is not considered. The user can specify a wind

shading factor from 0 to 1 to modify the effect of wind, where

0 represents complete shelter and 1 no shelter.

Long-Wave Radiation

The heat balance includes long-wave radiation exchange

between the underside of the sheathing and the attic floor,

and between the top surface of the roof and the sky. The

emissivity of all surfaces has been assumed at 0.9, a typical

value for many nonmetallic building materials (ASHRAE

1993). To facilitate the calculation, the nonlinear formulas

for radiation exchange were converted to linear equations in

the following way:

 

Q TT

T T T T T T FT T

x y, x y

x y x yx y x y

→ = −

= + + −≈ −

net εσ

εσ

( )

( )( )( ) ( )

4 4

2 2 (1)

where

Qx y, → net = net radiation from surface x to surface y

Tx, Ty = temperatures of surfaces x and y,

respectively (°R)

ε = surface emissivity

σ = Stefan–Boltzman constant (Btu/h⋅ft2⋅°R4)

The equation is linearized by treating the factor F as a constant and calculating F from initial values for Tx and Ty.
5

The error introduced by this linearization is relatively small

if the absolute temperatures are relatively close.

Snow

Measurements in an attic in Madison, Wisconsin, suggest

that the upper roof surface temperatures were held steady at

approximately 32°F (0°C) with significant snow accumulation. Winandy and Beaumont (1995) also reported this

phenomenon in tests on a series of small test chambers near

Madison, Wisconsin. Therefore, we included an alternate

heat balance for the sheathing that fixes the temperature of the

exterior roof surface at 32°F (0°C) when snow is present,

regardless of outdoor temperature or solar radiation.

Equations

Attic Air—The heat balance for the attic air (node 1) is given

by

 

A h T T fI V C T T

A

hTT

A

hTT

A

R

T T

I CV T T

cf f a hh p i a

r r,1 a r,1 r r,2 a r,2 es

es

a o

a pa a o

() ()

2 ( 2 ( )( )

 ( )

−+ −

= −+ −+ −

+ −

ρ

ρ

 ) (2)

where

Ac, Ar,

Aes = surface areas of ceiling, total roof area, and

 combined area of soffit and end walls (ft2)

Cp = specific heat of air (Btu/lb⋅°F)

f = fraction of house exfiltration that transfers

 into attic

hf, hr,n = convective heat transfer coefficients at

attic floor and underside of roof (surface 

n), respectively (Btu/h⋅ft2⋅°F)

Ih = house exfiltration rate (h−1)

Ia = attic ventilation rate, i.e., outdoor air

entering attic (h−1)

Re,s = average thermal resistance of end walls

and eaves (h⋅ft2⋅°F/Btu)

Ta = attic air temperature (°R)

Tf = temperature of attic floor surface (°R)

Ti = indoor air temperature (°R)

To = outdoor air temperature (°R)

Tr,1 = temperature of roof sheathing underside a (°R)

Tr,2 = temperature of roof sheathing underside b (°R)

Va = volume of attic space (ft3)

Vh = volume of house (ft3)

ρ = density of air (lb/ft3)

The convective heat transfer coefficients are calculated iteratively with the HCON subroutine. If the iterations do not

converge, the value of 1 Btu/h⋅ft2 is used and an error message is printed in the MESSAGE output file.

Attic Floor—The heat balance for the attic floor (node 2) is

 T T

R

hT T F T T F T T i f

c ff a

1

f r,1

2

f r,2 ( ) 2 ( ) 2 ( ) − = −+ − + − (3)

where 

Rc = thermal resistance of attic floor (ceiling)

(h⋅ft2⋅°F/Btu)

F1, F2 = radiative heat transfer coefficients between

attic floor and undersides of roof

(Btu/h⋅ft2⋅°F), with the following

definition:

 

 

F T T TT n n nn = + −++=

σ

(/ ) (/ ) ε ε

( )( ) , 1 111 2

f r

f

2 r,2

f r, (4)

where

εf = emissivity of attic floor surface

(assumed to be 0.9)

εr = emissivity of roof sheathing surface

(underside, assumed to be 0.9)

σ = Stefan–Boltzman constant (Btu/h⋅ft2⋅°R4)

To facilitate the calculation, radiation exchange between the

floor and the sheathing is approximated in a linear equation

and included in Equation (3). Radiation exchange between

the two sheathing surfaces is neglected because the aspect

ratio is usually low. Radiation exchange with end walls is

also neglected.

Sheathing—The heat balance at the underside of sheathing

(nodes 3a and 3b) is

T T

R

hT T

A

A

FT T

A

LW n

n n ,n

nn n

r, s,

r

ra r

c

r

f r,

r

h r,

( )

( ) 2 1,2

− = −

+ −+ =

 (5)

where

Lh = latent heat of vaporization (1,050 Btu/lb)

Rr = thermal resistance of roof (h⋅ft2⋅°F/Btu)

Ts,n = temperature of outside roof surface (°R)

Wr,n = rate of moisture adsorption into sheathing

(lb/h)

The approximate latent heat of sorption/desorption or condensation/evaporation is calculated from the moisture deposition rate on the sheathing surface, using a value of

1,050 Btu/lb. The adjusted heat balance and associated

moisture balance are solved iteratively until the latent heat

adjustment is within 5% of the previous value. If convergence is not reached, the latent heat is ignored and an error

message is printed in the MESSAGE file.

The heat balance at the top surface of sheathing (nodes 4a

and 4b) is

T T

R

I h h TT L n n n

nn n

s, r,

r

o, IR o s, IR ( )( ) 1, 2 − =+ + − − = α (6)

where

hIR,LIR = adjustments for infrared radiation

exchange with sky (see Eqs. (7), (8), and (9))

ho,n = convective heat transfer coefficient at

exterior roof surface n (Btu/h⋅ft2⋅°F)

In = total solar radiation incident on roof surface 

(Btu/h⋅ft2)

α = solar absorptance
6

The adjustments for infrared sky radiation were adapted from

Martin and Berdahl (1984) and take the effect of outdoor

temperature, dew point temperature, and cloud cover into

account. The base of the clouds is assumed at 1.24 mi

(2 km). The radiation exchange is approximated with corrections to the surface heat transfer coefficients and incident solar

radiation. In Equation (5), the contribution of sky radiation

was linearized to facilitate calculation. The terms LIR and hIR

are computed in the INFRED subroutine:

 

h T IR s o = 4 3 ε σ (7)

 

L T IR s o IR = − εσ ε 4( ) 1 (8)

where

εs = emissivity of roof shingles (assumed to be 0.9)

εIR = sky emissivity with clouds

The sky emissivity is calculated with equations from Martin

and Berdahl (1984). The emissivity for a clear sky is

 

 

ε

π

0

d d =+ + 









 + 



 



 0 711 0 56

100

0 73

100

0 013 2

24

2

. . . . cos T T t (9)

where

ε0 = emissivity of clear sky

Td = outdoor dew point temperature (°C)

t = time (h)

and the emissivity of the sky with clouds is

εε ε IR =+ − 0 0 0 784 1 .( ) C (10)

where C = total cloud cover as recorded by National

Climatic Center (Asheville, NC). Values of C range from

0 to 10.

Attic Moisture Balance

The attic moisture balance requires that the water vapor

transfer by convective air flow through the ceiling equal the

combined net loss of moisture from the attic air by attic

ventilation and transfer to the two sheathing surfaces:

 

fI V I V W W h h a a a a o r,1 r,2 ρω ω ρω ω ()( ) i −= −++ (11)

where

ωi, ωa, ωo = humidity ratios of indoor, attic, and

outside air, respectively (dry air basis)

Wr,n = rate of moisture adsorption into sheathing

surfaces (lb/h)

Water vapor transfer by diffusion through the ceiling is

ignored because it is generally much less than that by air

leakage through the ceiling. Moisture exchange with the

outside through the roof is also ignored because it is likely

to be much smaller than losses through ventilation.

Moisture transfer from the attic air to the two sheathing

surfaces is given by

W A r,n n h r D,r a r, 2 = − ( ) ω ω (12)

where

hD,r = air-surface moisture transfer coefficient

(set at 1.1 lb/h⋅ft2)

ωr = equilibrium humidity ratio of surface air in

equilibrium with wood sheathing surface

(dry air basis)

For the equilibrium humidity ratio ωr at the sheathing surface, we used the correlation developed by Cleary (1984):

 

ωr T = ++ + e B CM DM EM / ( ) A 23 (13)

where

T = wood surface temperature (°F)

M = moisture content of wood, expressed as a ratio

A = 28.6°F

B =−0.00049

C = 0.0172

D =−0.060

E = 0.076

Cleary based this correlation on isothermal sorption data for

untreated clear wood. It is not adjusted for any fire-retardant

or other treatment.

The sheathing is divided into three layers: a 1/16-in.

(1.6-mm) surface layer facing the attic, an adjacent 1/16-in.

(1.6-mm) layer, and a core layer representing the rest of the

sheathing. With plywood sheathing, the first two layers

combined represent the surface veneer layer that faces the

attic. Moisture diffusion from layer m to layer m + 1 (surface

layer is layer 1) is approximated with

 

 

w K

d d m,m m m m

m m + ++

= −

+ 1 = 1

1 2

1 2 w( )

( ) , ω ω/ (14)

where

wm,m+1 = moisture flux from layer m to

next layer (lb/h⋅ft2)

Kw = diffusion coefficient (lb/h⋅ft)

ωm, m+1 = equilibrium humidity ratio, layer m,

layer m + 1, respectively

dm,dm+1 = thickness of layer m, layer m + 1,

respectively (ft)

The equilibrium humidity ratio is calculated with Equation (13). We obtained approximate values for Kw from data

for isothermal moisture diffusion coefficients for clear wood

published by Choong (1965). Choong’s diffusion coefficient

applies to moisture content (MC) gradients rather than

humidity ratios. An approximate conversion yielded

 

K TT

 MM

w =− − +

+ −+

11 48 0 000083 0 0222

168 27 102 4 37 9

2

3 2

. . . 

. .. M

 (15)
7

where

T = average temperature of the two layers (°F)

M = average MC of the two layers, expressed as a ratio

Choong’s diffusion coefficient applies to clear wood under

isothermal conditions. The effect of glue lines, temperature

gradients, and fire retardants on diffusion coefficients is unknown. In addition, the coefficients only apply below fiber

saturation. We assumed that the coefficient remains constant

at MC above 30%.

Model Validation

Introduction

To assess the accuracy of temperature and moisture predictions with the FPL roof temperature and moisture model, we

obtained data collected by William Rose and Jeffrey Gordon

in a test building at the Building Research Council (BRC),

University of Illinois, Champaign–Urbana. Temperature and

moisture were measured in locations similar to the nodes in

the model, and most of the necessary input data were also

available, such as weather and solar radiation. Missing were

data for attic ventilation and absorptance of the shingles.

We compared predictions of sheathing temperatures, attic air

temperature, attic air RH, and sheathing MC with measured

data taken during the winter and summer of 1993.

Description of Test House

The test building was designed, configured, and instrumented to measure the performance of various roof/attic

assemblies under natural conditions. The 80-ft- (24-m-) long,

20-ft- (6-m-) wide building had one story and a gable roof.

There were 10 separate roof bays, of which 8 were instrumented; the roof bays at the gable ends were not instrumented. Of the instrumented bays, only data from bay 2 were

used for validation of the FPL roof. Bay 2 had black shingles, solid vinyl soffit panels, and no ridge vent. It had a flat

ceiling and therefore represented traditional attic construction.

There were no penetrations in the ceiling during the winter

and summer of 1993.

All the flat ceiling bays (attics) were constructed with standard roof trusses at a 5:12 pitch, 24 in. (61 cm) on center

(o.c.), with 2-ft (0.6-m) overhangs. The orientation of the

ridge was east to west. The roof was covered with single

thickness, triple-tab asphalt shingles over 15-lb felt underlayment. The roof sheathing material was 7/16-in. (11-mm)

oriented strandboard (OSB). The ceiling was 5/8-in.

(16-mm) gypsum drywall, overlaid with R-30 unfaced glass

fiber batt insulation. There was no vapor retarder. The flat

ceiling bays were thermally separated by 1-in. (25-mm) foilfaced polyisocyanurate foam panels; cracks were filled with

urethane foam.

Data Collection

Temperature

Roof surface temperatures were taken with type-T thermocouples. A resistance temperature detector (RTD) was used

as a reference junction. The temperatures were sampled every

10 min, averaged, and stored as hourly values.

In each bay, thermocouples were located at the following

sites:

• at underside of sheathing and between sheathing and felt

underlayment,

• on north- and south-facing slopes, and

• at low and high locations on the roof: 3 ft (0.9 m) away

from wall plate and ridge, respectively.

Attic air temperature was measured with one platinum RTD

in each bay.

Moisture Content and Humidity

Attic air humidity was measured with one polyamide humidity sensor in each bay. Although the error of measurement

should not have been more than 4% RH, there is reason to

believe that these sensors were not capable of handling the

wide temperature swings in the attic bays and did not

produce reliable data.

Humidity conditions in the sheathing were measured with

small, wood electrical resistance sensors modified from the

type described by Duff (1966). The modifications included

the use of conductive epoxy and a polyethylene “sleeve”

around part of the sensor. Sensors were inserted in the middle of the sheathing as well as installed on its surface.

Although the accuracy of these sensors is not very high

(TenWolde and Courville 1985), they are nevertheless the

most practical available equipment for research on building

moisture. The sensors were not individually calibrated, and

the data were recorded as electrical resistance. We converted

these data to MC, using generic calibration data for electrical

resistance wood moisture meters. We used measured MC

data for the north-facing sheathing of bay 2 only.

Other Instrumentation

On-site measurements were taken of wind speed and direction, air temperature, RH, total horizontal solar radiation,

and total infrared radiation. Barometric pressure, snow cover,

and rainfall were measured at a meteorology station 1/4 mi

(0.4 km) from the site.

Model Runs

We ran the FPL roof temperature and moisture model from

February 3 to March 28, 1993, and May 1 to June 30, 1993,

using measured data for outdoor temperature and humidity.

Solar radiation on the roof surfaces was determined by running the SUNDAT program with measured solar radiation

(horizontal surface) and the location of the test house and roof

slope and orientation.
8

We selected input parameters to match the conditions in bay

2 (black shingles and unvented). Because we did not know

the solar absorptance of the shingles, we did several runs

with varying values for absorptance. We then compared the

temperatures at the top surface of the sheathing with the

measured values to determine which value of absorptance

produced the closest match.

Selection of Input Values

Because the model could not handle the exact configuration

of the test house with separate attic/roof bays and shingles of

different colors, we selected a configuration that would most

closely approximate the actual situation, keeping in mind the

assumptions and equations used in the model. The heat

transfer coefficient at the exterior roof surface as calculated in

the model depends on the building dimensions. We therefore

used the full dimensions of the building as input, instead of

the dimensions of the individual bays. Heat loss through the

side walls separating each bay from its neighbors was assumed to be minimal because the bays are separated by 1-in.

(25-mm) foil-faced foam insulation and the measured air

temperature differences between bays were relatively small.2

Because the model assumes outdoor air temperature on the

exterior of the end walls, we specified a very high R-value

(R-100, RSI-18) for the end walls to simulate minimal heat

flow between bays. Table 2 lists the input parameters for bay

2, with the exception of solar absorptance.

Selection of Shingle Absorptance

We did several runs for bay 2 with varying absorptance

values and compared the temperatures at the top surface of the

sheathing with the measured values to determine which

absorptance value produced the closest match. We used

weather data for the first 15 days of May 1993. Reagan and

Acklam (1979) recommend the following values for solar

absorptance of shingles:

Shingle color Absorptance

Very light 0.25

Light 0.35

Medium 0.55

Dark 0.75

Very dark 0.90

Since bay 2 had dark shingles, it would have an expected

absorptance value of approximately 0.75.

Figure 4 compares the measured temperature on the top of

the south- and north-facing sheathing surfaces, respectively,

of bay 2 with the model results, using absorptance values of

0.6 and 0.65. Solar absorptance plays a role only during

 

2

Measured temperature differences between bays 1 and 2

averaged less than 6°F (3°C) during May and June 1993,

and even less during February and March 1993. Maximum

temperature difference was around 27°F (15°C). Temperature

differences between bays 2 and 3 were much less than that

between bays 2 and 1.

daylight hours and the comparison should therefore be limited to that period. For the purpose of selecting solar absorptance values, we focused on the peak temperatures during

mid-day.

Both absorptance values of 0.6 and 0.65 produced reasonably

accurate temperature profiles during mid-day on most days

and would be acceptable choices. Some differences between

model results and measurements were likely the result of

errors in the calculation of solar radiation on the south- and

north-facing roof surfaces, because the amount of diffuse

radiation was not measured, but estimated from the total

amount of horizontal radiation. In addition, wind is likely to

provide a discrepancy between model results and measurements. The 0.65 value produced slightly better results and

was therefore selected for the “dark” shingles. This compares

with the value of 0.75 recommended by Reagan and Acklam

(1979) for dark shingles.

Results

Attic Air Temperature—Winter

Figure 5 compares measured attic air temperatures in bay 2

from February 3 to 16, 1993, and air temperatures from the

FPL model. Most discrepancies were due to a time delay of

about 1 to 2 h in the measured temperatures, most likely the

result of heat storage in the building materials in the roof and

attic. If the time delay is ignored, the model predictions were

usually within 5°F (3°C), but greater discrepancies of

approximately 10°F (6°C) occurred, especially during the

night. As will be discussed later, the model consistently

overpredicted the extent of night-time cooling from infrared

Table 2—Input parameters used for model validation,

bay 2

Building length 80 ft (24 m)

Building width 20 ft (6 m)

Ceiling height 8 ft (2.4 m)

Roof overhang 2 ft (0.6 m)

Roof pitch 0.42

Sheathing thickness 0.44 in. (11 mm)

R-value of ceiling 11 h⋅ft2

°F/Btu (1.94 m2⋅K/W)

R-value of roof, total 1.04 h⋅ft2

⋅°F/Btu (0.183 m2⋅K/W)

R-value of roof sheathing 0.6 h⋅ft2

⋅°F/Btu (0.1 m2⋅K/W)

R-value of end walls 100 h⋅ft2

⋅°F/Btu (17.6 m2⋅K/W)

Wind shade factor 1

House ventilation rate 0.1 ach

Fraction of exhaust air to

 attic

0.01

Attic ventilation rate 0.01 ach

Indoor temperature, winter 70°F (21°C)

Indoor temperature, summer 75°F (24°C)

Indoor relative humidity 50%

Initial sheathing moisture

 content, February 1

9.3%

Initial sheathing moisture

 content, May 1

7.0%
9

Figure 4—Exterior sheathing surface temperature, bay 2, May 1–15, 1993. (a) south facing; (b) north-facing.
10

radiation to the sky. This, led to predicted attic air temperatures that were too low during the night, especially when

there were no clouds.

The weather data indicate that day 42 (February 11) was very

cloudy (little solar radiation) and windy, with temperatures

dropping throughout the day. This explains the steady low

attic temperatures during that day. Snow could have also

been a contributing factor, although snow cover data do not

indicate snow on the ground that day.

Snow cover on the ground on day 47 (February 16) caused

the model to hold the exterior surface temperature of the roof

at 32°F (0°C), which then held the attic air temperature

steady at about 35°F (2°C). Although this is a crude approach to the effect of snow cover, the model did provide a

fair approximation of measured temperatures. However, snow

cover data did not necessarily match actual snow cover on

the roof. The model was also unable to deal with partial

snow cover on the roof, leading to additional errors.

Attic Air Relative Humidity—Winter

Figure 6 shows a consistently large discrepancy between

predicted and measured RH of the attic air during the same

period (February 3 to 16, 1993). The measured humidity

data are inconsistent with the measured MC in the sheathing;

in private conversations, William Rose of the University of

Illinois expressed skepticism about the accuracy of the sensor

when used in the highly variable attic environment. We

therefore defer judgment on the accuracy of the model’s RH

predictions.

Sheathing Moisture Content—Winter

Figure 7 shows the measured MC on the underside of the

north-facing sheathing near the eaves and the ridge compared

with the results from the model. Both the MC at the under

surface and the average MC are shown. The results prompt

several observations:

• measured MC varies significantly with location,

• swings in measured MC suggest that values represent

MC of sheathing surface, rather than average MC,

• measured MC maxima and minima coincide with model

minima and maxima, respectively, and

• the difference between calculated average MC and measured MC is usually within 2%, fairly close when considering that estimated measurement error was 2% to

3% MC.

The discrepancy in maxima and minima of the calculated and

measured MC of the north-facing sheathing is explored further in Figure 8. Although the variations in the average

measured MC of the north-facing sheathing are different than

the simulated surface MC of the north-facing sheathing, there

is much better agreement between measured MC of the northfacing sheathing and the model results for the south-facing

Figure 5—Attic air temperature, bay 2, February 3–16, 1993; absorptance = 0.65.
11

Figure 6—Attic air relative humidity, bay 2, February 3–16, 1993; absorptance = 0.65.

Figure 7—Sheathing moisture content, bay 2, north-facing, February 3–16, 1993; absorptance = 0.65.
12

sheathing (i.e., the minima and maxima coincide). The

model results for the north- and south-facing sheathing

surfaces behave in this fashion because the model assumes

that most or all of the moisture driven from the south-facing

sheathing is deposited on the north-facing sheathing when

sunlight hits the south-facing roof. The measurements show

that the north-facing sheathing in fact behaves very similar to

the south-facing sheathing. This means that when sunlight

hits the roof, moisture is driven from both sides and is absorbed elsewhere in the attic (for example, rafters, floors, end

walls, insulation). The model only accounts for moisture in

the sheathing and does not provide for moisture storage

elsewhere. Figure 8 demonstrates the importance and role of

the temporary moisture storage in materials other than the

sheathing. The lack of such a mechanism in the FPL model

seriously hampers its ability to predict hourly changes in

sheathing MC and tends to exaggerate the difference in MC

between north- and south-facing sheathing.

Sheathing Temperature—Summer

Figure 9 compares measured temperatures on the top surface

of the plywood sheathing (immediately under the roof felt)

during May 1 to 15, 1993, with results from the model. In

contrast to the attic air temperature, there is no discernible

time shift between measured top surface temperatures and

model results. This supports the view that the time delay in

attic air temperatures was caused by heat storage in the

sheathing and other materials in the attic, which the model

does not take into account. On most days, measurements and

simulation results were reasonably close (within 5°F (3°C)),

except during clear nights and during mid-day on some days,

especially May 13 to 15.

The data show discrepancies during the night of up to 15°F

(8°C). This error is primarily due to an overestimate of

cooling by sky radiation. Research into alternative equations

for sky radiation did not reveal substantial differences between results from these equations. We therefore surmise that

the fault did not lie in the equations or calculations. Rather,

radiation loss to the sky was perhaps partially compensated

by radiation received from the surroundings, or the choice of

emissivity (0.9) for the roof surface was too high.

The discrepancies during mid-day could be the result of

errors in the estimated effect of wind cooling. Inspection of

the measured temperature and wind data did show an effect of

wind direction on roof surface temperatures, but wind direction is not considered in the model. The discrepancies could

also be caused by errors in the calculated radiation input data

for each roof surface. This calculation was prone to errors

partly because of inaccurate cloud cover data, which influence

the ratio of direct to diffuse radiation.

To determine the potential for thermal degradation of the

sheathing, the model can be used to generate temperature–

duration distributions. Figures 10 and 11 show temperature–

time distributions from May 1 through June 30, 1993, for

the top surface of north- and south-facing sheathing, as compiled from the measured data and the model results. For

the north side, the model underpredicted the maximum

Figure 8—Sheathing moisture content, bay 2, north- and south-facing, February 3–16, 1993; absorptance = 0.65.
13

Figure 9—Sheathing temperatures, bay 2, top surface, May 1–15, 1993; absorptance = 0.65.

............................................

Number of hours (h)

0 50 100 150 200 250 300

Temperature (°F)

150-160

140-150

130-140

120-130

110-120

100-110

 90-100

 80-90

 70-80

 60-70

 50-60

 40-50

 <40

........................................................

.........................................................

.....................................................

.......................

.................

..............

....

.........................

.............................

......................

..................

....................

...................

.................

..................

...................

...................

...........

..............

.....

.........

.

....

North sheathing

Measured

Model

Figure 10—Temperature–time distribution, bay 2, top of

north sheathing, May 1–15, 1993; absorptance = 0.65.

 Tc= [TF − 32]/1.8.

...........................................

Number of hours (h)

0 50 100 150 200 250 300

Temperature (°F)

150-160

140-150

130-140

120-130

110-120

100-110

 90-100

 80-90

 70-80

 60-70

 50-60

 40-50

 <40

.......................................................

...........................................................

.....................................................

.........................

.................

..............

....

........................

.............................

..................

.................

................

...................

...............

..................

...............

...................

..............

.............

.........

.........

...

....

South sheathing

Measured

Model

Figure 11—Temperature–time distribution, bay 2, top of

south sheathing, May 1–15, 1993; absorptance = 0.65.

Tc= [TF − 32]/1.8.
14

temperature and frequency of temperatures above 120°F

(49°C) (Fig. 10). The model overpredicted the number of

hours below 60°F (16°C), primarily because of its overprediction of night sky cooling. On the south side (Fig. 11), the

model appears quite accurate for temperatures above 120°F

(49°C), generally underpredicted the frequency of mid-range

temperatures between 60°F and 120°F (15°C and 49°C), and

overpredicted the occurrence of temperatures below 60°F

(49°C). Usually, for thermal degradation concerns, the frequency of occurrence of high temperatures is of the most

interest, and for that purpose the model results appear to have

acceptable accuracy.

Sheathing Moisture Content—Summer

Figure 12 compares the MC of the north-facing sheathing as

calculated from the model (average and surface MC) with the

measured MC (average of two locations) during the last

2 weeks of June 1993. The conclusions are similar to those

made earlier. On average, the measured MC was close to the

model results (about 1% lower), especially when considering

an estimated error in the measured MC of at least 2%. However, the data show that moisture storage in the attic needs to

be added to the model to accurately simulate the moisture

behavior of the north-facing sheathing.

Conclusions

The model validation described in this report is not a complete test of the model’s accuracy over the full range of conditions and can therefore provide only an indication of its

performance. However, we feel that the results justify the

following conclusions:

• Comparison of measured data and data from the FPL roof

temperature and moisture model shows that the model is

capable of providing reasonably accurate estimates for

temperatures of the roof sheathing and attic air, with the

difference usually within 5°F (3°C). However, heat storage

effects, which are not accounted for in the model, often

cause a time shift (delay) of 1 to 2 h in attic air temperatures.

• The model can accurately predict the frequency of occurrence of high roof sheathing temperatures (above 120°F

(49°C)) during summer, but its accuracy is highly dependent on the values chosen for solar absorptance and

emissivity of the roof shingles. However, data on emissivity and solar absorptance are usually not available for specific shingles.

• The model consistently overpredicts the extent of nighttime cooling from sky radiation losses, leading to

Figure 12—Sheathing moisture content, bay 2, north-facing, May 1–15, 1993; absorptance = 0.65.
15

predicted temperatures that are too low. A possible

explanation is that the roof surfaces may receive night-time

radiation from the surrounding buildings, lights, and other

objects, which may not be included in the measured radiation input data, or that the choice of 0.9 for the emissivity

of the shingles is too high.

• The model’s treatment of the effect of snow cover is too

simplistic, but no better alternatives are apparent for simulating this very complex behavior.

• We were unable to determine the model’s accuracy of attic

air relative humidity predictions because of the suspect

quality of measured relative humidity data.

• The model appears capable of predicting average moisture

conditions in the sheathing with reasonable accuracy, generally within 1%, when moisture content is not excessively high or low. Moisture behavior at high and very

low moisture content and humidity conditions was not

tested.

• Hourly moisture behavior is not well represented by the

model as daily or seasonal behavior, especially when considering the north-facing sheathing. To improve the

model’s performance, a mass with moisture storage needs

to be added to the model to represent moisture sorption in

the attic floor, insulation, structural members, and other

hygroscopic materials.

• The model would further benefit from

– verification with data that include measured

emissivity and solar absorptance of the shingles,

– addition of thermal mass in the attic and roof,

– better algorithms to calculate direct and diffuse

solar radiation, and

– verification for roof with east–west orientation.

References

ASHRAE. 1993. ASHRAE Handbook of Fundamentals.

Atlanta, Ga: American Society of Heating, Refrigerating, and

Air-conditioning Engineers.

Burch, D.M. 1992. Controlling moisture in the roof cavities

of manufactured housing. NISTIR 4916. Gaithersburg, MD:

National Institute of Standards and Technology,

Burch, D.M.; Luna, D.E. 1980. A mathematical model for

predicting attic ventilation rate required for preventing condensation on roof sheathing. ASHRAE Transactions 86(1).

Atlanta, Ga: American Society of Heating, Refrigerating, and

Air-conditioning Engineers.

Burch, D.M.; Lemay, M.R.; Rian, B.R.; Parker, E.J.

1984. Experimental validation of an attic condensation

model. ASHRAE Transactions 90(2). Atlanta, Ga: American

Society of Heating, Refrigerating, and Air-conditioning

Engineers.

Choong, E.T. 1965. Diffusion coefficients of softwoods by

steady-state and theoretical methods. Forest Products

Journal. 15(1): 21–27.

Cleary, P. 1984. Moisture control by attic ventilation—an

in situ study. Report LBL–18062. Berkeley, CA: Lawrence

Berkeley Laboratory.

Duff, J.E. 1966. A probe for accurate delamination of moisture content of wood products in use. Res. Note FPL–RN–

0142. Madison, WI: U.S. Department of Agriculture, Forest

Service, Forest Products Laboratory.

Gorman, T.M. 1987. Modeling attic humidity as a function

of weather, building construction, and ventilation rates.

Syracuse, NY: State University of New York, College of

Environmental Science and Forestry. Ph.D. thesis.

Holman, J.P. 1981. Heat transfer. 5th ed. New York, NY:

McGraw–Hill.

LeVan, S.; Collet, M. 1989. Choosing and applying fireretardant-treated plywood and lumber for roof designs. Gen.

Tech. Rep. FPL–GTR–62. Madison, WI: U.S. Department

of Agriculture, Forest Service, Forest Products Laboratory.

Martin, M.; Berdahl, P. 1984. Characteristics of infrared

sky radiation in the United States. Solar Energy. 33:

321–336.

Reagan, J.A.; Acklam, D.M. 1979. Solar reflectivity of

common roofing materials and its influence on the roof heat

gain of typical southwestern residences. In: Summer attic and

whole-house ventilation. NBS Special Publication 548.

Gaithersburg, MD: U.S. Department of Commerce, National

Bureau of Standards: 1–13.

Rowley, F.B.; Algren, A.B.; Lund, C.E. 1939. Condensation of moisture and its relation to building construction and

operation. American Society of Heating and Ventilating

Engineers Transactions. 45: 231–252.

TenWolde, A.; Carll, C. 1992. Effect of cavity ventilation

on moisture in walls and roofs. In: Thermal performance of

the exterior envelope V. ASHRAE Special Publication.

Atlanta, Ga: American Society of Heating, Refrigerating, and

Air-conditioning Engineers: 555–562.

TenWolde, A.; Courville, G.E. 1985. Instrumentation for

measuring moisture in building envelopes. Atlanta, Ga:

American Society of Heating, Refrigerating, and Airconditioning Engineers Transactions: 91, part 2.

Wilkes, K.E. 1989. Model for roof thermal performance.

ORNL/CON–274. Oak Ridge, TN: Oak Ridge National

Laboratory.

Winandy, J.E.; Beaumont, R. 1995. Roof temperatures in

simulated attics. Res. Pap. FPL–RP–543. Madison, WI:

U.S. Department of Agriculture, Forest Service, Forest

Products Laboratory.
16

Nomenclature

Ac = surface areas of ceiling (ft2)

Aes = combined surface area of soffit and end walls (ft2)

Ar = total roof surface area (ft2)

C = total cloud cover (0 to 10)

Cp = specific heat of air, Btu/lb⋅°F

dm = thickness of sheathing layer m (ft)

f = fraction of house exfiltration that transfers into attic

F1, F2 = radiative heat transfer coefficients between attic floor and undersides of roof (Btu/h⋅ft2⋅°F)

hD,r = air-surface moisture transfer coefficient (set at 1.1 lb/h⋅ft2)

hf = convective heat transfer coefficient at attic floor (Btu/h⋅ft2⋅°F)

hIR = adjustment in heat transfer coefficient for infrared radiation exchange with sky (Btu/h⋅ft2⋅°F)

ho,n = convective heat transfer coefficient at exterior roof surface n (Btu/h⋅ft2⋅°F)

hr,n = convective heat transfer coefficient at underside of roof (surface n) (Btu/h⋅ft2⋅°F)

Ia = attic ventilation rate, i.e., outdoor air entering attic (h−1)

Ih = house exfiltration rate (h−1)

In = total solar radiation incident on roof surface (Btu/h⋅ft2)

Kw = diffusion coefficient (lb/h⋅ft)

LIR = partial adjustment for infrared radiation exchange with sky (Btu/h⋅ft2)

Lh = latent heat of vaporization (1,050 Btu/lb)

M = moisture content of wood, unitless

Rc = thermal resistance of attic floor (ceiling) (h⋅ft2⋅°F/Btu)

Re,s = average thermal resistance of end walls and eaves (h⋅ft2⋅°F/Btu)

Rr = thermal resistance of roof (h⋅ft2⋅°F/Btu)

t = time (h)

T = temperature (°F)

Ta = temperature of attic air (°R)

Td = temperature of outdoor dew point (°C)

Tf = temperature of attic floor surface (°R)

Ti = temperature of indoor air (°R)

To = temperature of outdoor air (°R)

Tr,1 = temperature of underside of roof sheathing a (°R)

Tr,2 = temperature of underside of roof sheathing b (°R)

Ts = temperature of outside roof surface (°R)

Va = volume of attic space (ft3)

Vh = volume of house (ft3)

wm,m+1 = moisture flux from layer m to next layer (lb/h⋅ft2)

Wr,n = rate of moisture adsorption into sheathing surfaces (lb/h)

α = solar absorptance

εf = emissivity of attic floor surface

εIR = emissivity of sky with clouds

εr = emissivity of roof sheathing surface (underside)

εs = emissivity of roof shingles (assumed to be 0.9)

ε0 = emissivity of clear sky

ρ = density of air (lb/ft3)

σ = Stefan–Boltzman constant (Btu/h⋅ft2⋅°R4)

ωa = humidity ratio of attic air, dry air basis

ωi = humidity ratio of indoor air, dry air basis

ωm = equilibrium humidity ratio, sheathing layer m, dry air basis

ωo = humidity ratio of outside air, dry air basis

ωr = equilibrium humidity ratio of surface air in equilibrium with wood sheathing surface, dry air basis
17

Appendix A—Convective Heat Transfer Correlations

The following convective heat transfer correlations were adapted from Wilkes (1989).

A. Natural convection, horizontal or nearly horizontal surface (tilt angle φ < 2°):

a. Heat flow up:

Nu = 0.54 Ra1/4 for Ra < 8 × 106

Nu = 0.15 Ra1/3 for Ra > 8 × 106

b. Heat flow down:

Nu = 0.58 Ra1/5

B. Natural convection, tilted surface (tilt angle φ > 2°):

 a. Heat flow up:

Nu = 0.56 (Ra cos(φ))1/4 for Ra/Pr < Grc

Nu = 0.14 (Ra1/3 − (GrcPr)1/3)

 + 0.56 (Ra cos(φ))1/4 for Ra/Pr > Grc

with Grc = 1 × 106 for φ < 15°

 Grc = 10(φ/(1.1870 + 0.0870φ)) for 15° < φ < 75°

 Grc = 5 × 109 for φ > 75°

 b.Heat flow down:

Nu = 0.56 (Ra cos(φ))1/4

C. Forced convection:

Nu = 0.664 Pr1/3 Re1/2 for Re < 5 × 105

Nu = Pr1/3 (0.037 Re4/5 − 850) for Re > 5 × 105

Combined forced and natural convection:

hcombined = (hforced

3

 + hnatural

3

)

1/3

Nu = Nusselt number = hL/k

Ra = Raleigh number = gβρCp∆TL3/νk

Gr = Grasshof number = Ra/Pr

Pr = Prandtl number = ν/α

Re = Reynolds number = vL/ν

h = convective heat transfer coefficient

L = characteristic length of “plate” (average of length and width)

k = thermal conductivity of air

g = acceleration of gravity

β = volume coefficient of expansion of air

ρ = density of air

Cp = specific heat of air

∆T = temperature difference between surface and air

ν = kinematic viscosity of air

α = thermal diffusivity of air

v = velocity of air (wind speed)
18

Appendix B—Program Listing for FPL Model

PROGRAM AHM

$DEBUG

C VERSION 7A

C Date: April 15, 1992 revision of COMMON block

C Works with MOIST7A, PARAM7A, HEAT7A, HCON6

C Requires input of dewpoint temperatures and cloudcover data

C in range of 0-10

C New convective heat transfer and infrared sky radiation

C Added wind speeds from weather data (in knots)

C Added option to bypass latent heat calculations

C Changed output file format to single column (version 5C)

C Version 6: WEXTRA introduced

C Updated August 13, 1992: changed name of message file to MESSAGE.PRN

C Updated August 15, 1992: streamlining

C Major revision August 17, 1992: All hourly input parameters (except

C snow cover) are linearly interpolated between hours in 10-minute

C increments.

C Major revision February 28, 1993: change to 3-layer sheathing (from 5)

C Revision April 15, 1993: made sheathing thickness and R-value variable

C input parameters.

C

 DOUBLE PRECISION WMC(6),WEXTRA

 COMMON A,AC,AES,AF,AR,ALR,ALF,F,I1,I2,IA,IAF,IFLAG,IH,ITMOIST,

 1JJ,KK,L,NFAN,RC,RES,RHMAX,RR,RS,VA,VH,RHA,RHI,PITCH,VENT,PHI,

 2VLHA,VLHB,T(10),TARR(184,24),RAARR(184,24),WIND(184,24),WSHADE,

 3RBARR(184,24),SKYCOV(184,24),SNOCOV(184,24),TATT(184,24), 

 4TRA(185,24),TRB(185,24),OUTDP(184,24),WMCA(185,24), 

 5WMCB(185,24),RHAT(184,24),FAN(184,24),TEMPA(184,24),TEMPB(184,24), 

 6TEMPOA(184,24),TEMPOB(184,24),AVEMC(1,24),D,IIILAT,SLOPE(6),JJJ

 REAL I1,I2,IA,IAF,IH,L

 CHARACTER*12 FLTEMP,FLSOLA,FLSOLB,FLCLD,FLSNOW,FLDEWP,FLWND

C A ROOF SURFACE ABSORPTANCE 

C AC AREA OF CEILING, SQ. FT. 

C AES COMBINED AREA OF ENDWALLS AND SOFFITS, SQ. FT. 

C AR AREA OF ROOF. SQ. FT. 

C AF ATTIC FAN VENTILATION RATE, CFM 

C ALF CHARACTERISTIC LENGTH OF FLOOR (FOR HCON SUBROUTINE)

C ALR CHARACTERISTIC LENGTH OF ROOF (FOR HCON SUBROUTINE)

C DC DISTANCE FROM FLOOR TO CEILING, FT. 

C DL LENGTH OF HOUSE, FT. 

C DOH AMOUNT OF ROOF OVERHANG, FT. 

C DW WIDTH OF HOUSE, FT. 

C D THICKNESS OF ROOF SHEATHING (set at 0.47" in PARAM)

C DWA WOOD-TO-WOOD MOISTURE TRANSFER COEFFICIENT FOR SHEATHING 

C "A", LB/FT.*HR 

C DWB WOOD-TO-WOOD MOISTURE TRANSFER COEFFICIENT FOR SHEATHING 

C "B", LB/FT.*HR 

C F FRACTION OF HOUSE EXFILTRATION THAT TRANSFERS INTO ATTIC 

C HDR AIR-TO-WOOD SURFACE MOISTURE TRANSFER COEFFICIENT, 

C LB/FT.SQ.*HR 

C HF CONVECTIVE HEAT TRANSFER COEFFICIENT AT ATTIC FLOOR, 

C BTU/H*FT.SQ.*F. 

C HRA,
19

C HRB CONVECTIVE HEAT TRANSFER COEFFICIENT AT UNDERSIDE OF ROOF, 

C BTU/H*FT.SQ.*F. 

C HRI HUMIDITY RATIO OF INDOOR AIR, (UNITLESS) 

C HRO HUMIDITY RATIO OF OUTSIDE AIR, (UNITLESS) 

C HRRA HUMIDITY RATIO OF INTERIOR SURFACE OF SHEATHING "A" 

C HRRB HUMIDITY RATIO OF INTERIOR SURFACE OF SHEATHING "B" 

C HOA,

C HOB CONVECTIVE HEAT TRANSFER COEFFICIENT AT ROOF SURFACE, 

C BTU/H*FT.SQ.*F. 

C I1 SOLAR RADIATION ON ROOF SURFACE "A", BTU/H*FT.SQ. 

C I2 SOLAR RADIATION ON ROOF SURFACE "B", BTU/H*FT.SQ. 

C IA NATURAL ATTIC VENTILATION RATE, VOLUME CHANGES/H. 

C IAF ATTIC FAN VENTILATION CAPACITY, VOLUME CHANGES/H. 

C IH HOUSE EXFILTRATION RATE, VOLUME CHANGES.H. 

C L RADIATION LOSS FROM ROOF SURFACE TO SKY, BTU/H*FT.SQ. 

C NFAN VARIABLE USED TO INDICATE NUMBER OF QUARTER HOURS FAN IS 

C TO BE RUN DURING AN HOUR'S CALCULATIONS (INTEGER) 

C PITCH SLOPE OF ROOF, RISE/RUN 

C PHI ROOF ANGLE, DEGREES

C RC THERMAL RESISTANCE OF ATTIC FLOOR, H*FT.SQ.*F./BTU 

C RES AVERAGE THERMAL RESISTANCE OF ENDWALLS AND EAVES, 

C H*FT.SQ.*F./BTU

C RS THERMAL RESISTANCE OF ROOF SHEATHING, H*FT.SQ.*F/BTU

C RHA RELATIVE HUMIDITY OF THE ATTIC 

C RHI RELATIVE HUMIDITY OF THE INDOOR LIVING SPACE 

C DPO DEWPOINT TEMPERATURE OF THE OUTSIDE AIR, R.

C RR THERMAL RESISTANCE OF ROOF, INCL. SHINGLES, H*FT.SQ.*F./BTU

C T(1) ATTIC AIR TEMPERATURE, R. 

C T(2) ATTIC FLOOR TEMPERATURE, R. 

C T(3) TEMPERATURE OF UNDERSIDE OF ROOF "A" SHEATHING, R. 

C T(4) TEMPERATURE OF UNDERSIDE OF ROOF "B" SHEATHING, R. 

C T(5) TEMPERATURE OF OUTSIDE SURFACE OF ROOF "A", R. 

C T(6) TEMPERATURE OF OUTSIDE SURFACE OF ROOF "B", R. 

C T(7) TEMPERATURE OF THE INDOOR LIVING SPACE, R. 

C T(8) TEMPERATURE OF OUTSIDE AIR, R. 

C VA VOLUME OF THE ATTIC SPACE, CU.FT. 

C VENRAT RATIO OF ATTIC VENT AREA/CEILING AREA 

C VENT TOTAL ATTIC VENTILATION RATE, VOLUME CHANGES/H. 

C VH VOLUME OF HOUSE LIVING SPACE, CU.FT.

C WEXTRA SURPLUS MOISTURE, LB/HR

C WIND WINDSPEED (WEATHER DATA), KNOTS

C WMC(1) MOISTURE CONTENT OF UNDERSIDE SURFACE SHEATHING "A" 

C WMC(2) MOISTURE CONTENT OF UNDERSIDE SURFACE SHEATHING "B" 

C WMC(3) MOISTURE CONTENT OF SHEATHING "A" LAYER 2

C WMC(4) MOISTURE CONTENT OF SHEATHING "B" LAYER 2

C WMC(5) MOISTURE CONTENT OF SHEATHING "A" LAYER 3

C WMC(6) MOISTURE CONTENT OF SHEATHING "B" LAYER 3

C 

C The first step in the program is to read in the house parameters: 

 CALL PARAM(WMC)

C 

C 

C All sheathing layers are assigned a moisture content 

C equal to that of the initial moisture contents: 

 DO 20 K=3,5,2
20

 WMC(K)=WMC(1) 

 WMC(K+1)=WMC(2) 

 20 CONTINUE

C 

C The capacity of the attic fan is translated into volume changes 

C per hour: 

 IAF=AF*60./VA 

 IRHMAX=0

C 

 WRITE (*,99) 

 99 FORMAT (' ENTER NUMBER OF DAYS ') 

 READ (*,*) IDAYS 

C 

C The fan control variables are set to artificially high values so 

C as not to interfere with operation of the program. If a humidity 

C or photoperiod fan control is desired, these variables will be 

C reassigned in the next step: 

 RHMAX=10.0 

 BTU=1000.0 

C 

C Menu-driven selection of fan control devices is displayed on the 

C screen and control variables are assigned: 

 100 WRITE (*,101) 

 101 FORMAT (' HUMIDISTAT-CONTROLLED ATTIC FAN? (0=NO, 1=YES)') 

 READ (*,*) HCONT 

 IF (HCONT) 120,120,102 

 102 WRITE (*,103) 

 103 FORMAT (' ENTER MAXIMUM HUMIDITY DESIRED (DECIMAL)') 

 READ (*,*) RHMAX

 PCONT=0

 TCONT=0

 GOTO 138

C 

 120 WRITE (*,121) 

 121 FORMAT (' PHOTOPERIOD-CONTROLLED ATTIC FAN? (0=NO, 1=YES)') 

 READ (*,*) PCONT 

 IF (PCONT) 130,130,122 

 122 WRITE (*,123) 

 123 FORMAT (' ','ON WHICH ROOF SURFACE IS PHOTOSENSOR LOCATED?') 

 WRITE (*,124) 

 124 FORMAT (' ','0) SURFACE "A" OR 1) SURFACE "B" (ENTER NUMBER)') 

 READ (*,*) ISURF 

 WRITE (*,125) 

 125 FORMAT (' ENTER MINIMUM BTU/SQFT/H NEEDED TO TRIGGER FAN') 

 READ (*,*) BTU

 TCONT=0

 GOTO 138

C 

 130 WRITE (*,131) 

 131 FORMAT (' TIMER-CONTROLLED ATTIC FAN? (0=NO, 1=YES)') 

 READ (*,*) TCONT 

 IF (TCONT) 138,138,132 

 132 WRITE (*,133) 

 133 FORMAT (' ENTER HOUR TO BEGIN FAN (INTEGER)') 

 READ (*,*) IHOUR1 
21

 134 WRITE (*,135) 

 135 FORMAT (' ENTER LAST HOUR FAN IS TO RUN (INTEGER)') 

 READ (*,*) IHOUR2 

C 

 138 WRITE (*,136)

 136 FORMAT (' INCLUDE LATENT HEAT EFFECTS? (0=NO, 1=YES)'/

 A' Note: inclusion greatly increases run time!')

 READ (*,*) IIILAT

C 

C Specify input file names

 WRITE (*,400)

400 FORMAT (' ENTER 7 FILENAMES FOR',/,' DRYBULB TEMP',/,' SOLAR',

 1' RADIATION, A AND B',/,' CLOUDCOVER',/,' SNOWCOVER',/,

 2' DEWPOINT TEMP',/,' WINDSPEED')

 READ (*,410) FLTEMP

410 FORMAT (A12)

 READ (*,410) FLSOLA

 READ (*,410) FLSOLB

 READ (*,410) FLCLD

 READ (*,410) FLSNOW

 READ (*,410) FLDEWP

 READ (*,410) FLWND

C 

C The data files are read from disk. To facilitate operation of the 

C program in its present format, the names of the files are given 

C here rather than asked for as input: 

C 

C Outdoor temperature: 

 OPEN (4,FILE=FLTEMP) 

 READ (4,*) ((TARR(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (4,STATUS='KEEP') 

C 

C Solar radiation on roof surface "A": 

 OPEN (4,FILE=FLSOLA) 

 READ (4,*) ((RAARR(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (4,STATUS='KEEP') 

C 

C Solar radiation roof surface "B": 

 OPEN (4,FILE=FLSOLB) 

 READ (4,*) ((RBARR(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (4,STATUS='KEEP') 

C 

C Cloud cover data: 

 OPEN (4,FILE=FLCLD) 

 READ (4,*) ((SKYCOV(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (4,STATUS='KEEP') 

C 

C Depth of snow on the roof: 

 OPEN (4,FILE=FLSNOW) 

 READ (4,*) ((SNOCOV(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (4,STATUS='KEEP') 

C 

C Outdoor dewpoint temperature:

 OPEN (4,FILE=FLDEWP) 

 READ (4,*) ((OUTDP(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 
22

 CLOSE (4,STATUS='KEEP')

C

C Windspeed

 OPEN (4,FILE=FLWND) 

 READ (4,*) ((WIND(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (4,STATUS='KEEP')

C 

C Open the message file

 OPEN (6,FILE='MESSAGE.PRN',STATUS='UNKNOWN')

C 

C Initialize surface and attic air temperatures

 T(2)=TARR(1,1)+459.67

 T(1)=T(2)

 T(3)=T(2)

 T(4)=T(2)

 T(5)=T(2)

 T(6)=T(2)

C Initialize surplus moisture

 WEXTRA=0.0D0

C

C The beginning of the DO loop to perform each hour's calculations: 

 DO 210 JJ=1,IDAYS

 WRITE (*,139) JJ 

 139 FORMAT (' ','DAY ',I3) 

 DO 200 KK=1,24 

 WRITE (6,140) JJ,KK

 140 FORMAT (' ','DAY ',I3,' HOUR ',I2) 

C

C Calculate 10-minute increment for hourly inputs:

 NXTDAY=JJ

 NXTHR=KK+1

 IF (KK.EQ.24) THEN

NXTDAY=JJ+1

NXTHR=1

IF (JJ.EQ.IDAYS) THEN

 NXTDAY=JJ

 NXTHR=24

ENDIF

 ENDIF

 SLOPE(1)=(TARR(NXTDAY,NXTHR)-TARR(JJ,KK))/6

 SLOPE(2)=(RAARR(NXTDAY,NXTHR)-RAARR(JJ,KK))/6

 SLOPE(3)=(RBARR(NXTDAY,NXTHR)-RBARR(JJ,KK))/6

 SLOPE(4)=(SKYCOV(NXTDAY,NXTHR)-SKYCOV(JJ,KK))/6

 SLOPE(5)=(OUTDP(NXTDAY,NXTHR)-OUTDP(JJ,KK))/6

 SLOPE(6)=(WIND(NXTDAY,NXTHR)-WIND(JJ,KK))/6

C

C Initializing ventilation rates and flags:

 IF (IRHMAX.EQ.1) THEN

 NFAN=1

 VENT=IAF

 ELSE

 NFAN=0

 VENT=IA

 ENDIF

C
23

C The next section, to line 150, determines if an attic fan control 

C device was selected, and if so determines the attic ventilation 

C rate appropriate for the type of control and the attic condition. 

C 

 IF (PCONT) 144,144,141 

 141 IF (ISURF) 142,142,143 

 142 IF ((RAARR(JJ,KK) .GT. BTU) .AND. (SNOCOV(JJ,KK) .LT. 0.5)) GO TO 

 A 147 

 143 IF ((RBARR(JJ,KK) .GT. BTU) .AND. (SNOCOV(JJ,KK) .LT. 0.5)) GO TO 

 A 147 

 144 IF ((KK .GE. IHOUR1) .AND. (KK .LE. IHOUR2)) GO TO 147 

C 

 145 GO TO 148 

 147 NFAN=6

 VENT=IAF

C 

 148 CONTINUE

C

C 

C The beginning of the DO loop to perform the 6 iterations within

C each hour:

 DO 180 JJJ=1,6

C 

C Initialize latent heat:

 VLHA=0.0 

 VLHB=0.0 

C The heat balance subroutine is called: 

 155 CALL HEAT

C 

C

C More initializing:

 IFLAG=0 

 ITMOIST=0 

C The moisture balance subroutine is called: 

 160 CALL MOIST(WMC,WEXTRA)

C

 IF (NFAN.EQ.6) GO TO 180

C Check if humidistat-controlled fan should be turned on next period

 IF (RHA.GT.RHMAX) THEN

 VENT=IAF

 IRHMAX=1

 IF (JJJ.LT.6) THEN

 NFAN=NFAN+1

 VLHA=0.0

 VLHB=0.0

 CALL HEAT

 ENDIF

 ELSE

 VENT=IA

 IF (IRHMAX.EQ.1.AND.JJJ.LT.6) THEN

 VLHA=0.0

 VLHB=0.0

 CALL HEAT

 ENDIF 

 IRHMAX=0
24

 ENDIF

C 

 180 CONTINUE

C

C The fraction of the hour that the fan was operated is calculated 

C and assigned to the output matrix: 

 190 FAN(JJ,KK)=NFAN/6

C 

C 

 200 CONTINUE 

 210 CONTINUE 

C 

C Close message file

 CLOSE (6)

C 

C Output files are opened and written to disk: 

 OPEN (7,FILE='XTAXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TEMPA(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XTBXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TEMPB(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XTCXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TEMPOA(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XTDXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TEMPOB(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XTTXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TATT(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XRHXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((RHAT(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='WMCA.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((WMCA(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XMAXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TRA(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='WMCB.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((WMCB(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='XMBXXXXX.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((TRB(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

 OPEN (7,FILE='FAN.PRN',STATUS='UNKNOWN') 

 WRITE (7,300) ((FAN(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS) 

 CLOSE (7,STATUS='KEEP') 

C 

 300 FORMAT (' ',4416(F7.3,/,' ')) 

C 

 STOP 

 END 
25

 SUBROUTINE PARAM(WMC)

$DEBUG

C VERSION 7B

C May 7, 1994, corrected June 3, 1994

C As 7A but with correction of error in reading of windshade factor WSHADE 

C Revised April 16, 1993: made sheathing R and D variable, eliminated VENRAT

C Revised February 28, 1993, Anton TenWolde

C As version 6A but with 3-layered sheathing.

C WARNING: SHEATHING THICKNESS IS SET AT 0.47 INCH

C History:

C Added: windspeed, windshade factor, char. lengths (1/13/90)

C Deleted: heat transfer coefficients (1/13/90) [handled by HCON]

C Corrected attic volume calculations (8/21/91)

C Revised August 17, 1992

C

 DOUBLE PRECISION WMC(6)

 COMMON A,AC,AES,AF,AR,ALR,ALF,F,I1,I2,IA,IAF,IFLAG,IH,ITMOIST,

 1JJ,KK,L,NFAN,RC,RES,RHMAX,RR,RS,VA,VH,RHA,RHI,PITCH,VENT,PHI,

 2VLHA,VLHB,T(10),TARR(184,24),RAARR(184,24),WIND(184,24),WSHADE,

 3RBARR(184,24),SKYCOV(184,24),SNOCOV(184,24),TATT(184,24),

 4TRA(185,24),TRB(185,24),OUTDP(184,24),WMCA(185,24),

 5WMCB(185,24),RHAT(184,24),FAN(184,24),TEMPA(184,24),TEMPB(184,24),

 6TEMPOA(184,24),TEMPOB(184,24),AVEMC(1,24),D,IIILAT,SLOPE(6),JJJ

 REAL I1,I2,IA,IAF,IH,L

 CHARACTER*12 FNAME

C

C

 WRITE (*,20)

 20 FORMAT (' ENTER NAME OF HOUSE PARAMETER FILE')

 READ (*,21) FNAME

 21 FORMAT (A12)

 OPEN (4,FILE=FNAME)

 READ (4,*) DL,DW,DC,DOH,PITCH,D,F,WSHADE,A,IH,RC,RES,RR,RS,

 1T(7),RHI,WMC(1),WMC(2),AF,IA

 CLOSE (4)

 26 DOW=DW+2*DOH

 AC=DL*DW

 ALF=(DW+DL)/2

 EDGE=SQRT(DOW**2+(PITCH*DOW)**2)

 AR=DL*EDGE

 ALR=(DL+EDGE/2)/2

 AES=(2*DOH*DL)+2*PITCH*(DOW/2)**2

 VA=DL*PITCH*(DOW/2)**2

 VH=AC*DC

 PHI=ATAN(PITCH)*180/3.14159265

 TRA(1,1)=WMC(1)

 TRB(1,1)=WMC(2)

 WMCA(1,1)=WMC(1)

 WMCB(1,1)=WMC(2)

C

 WRITE (*,27)

 27 FORMAT (' HOUSE PARAMETERS ARE AS FOLLOWS:')

 WRITE (*,28)
26

 28 FORMAT ('0DIMENSIONS:')

 WRITE (*,31)

 31 FORMAT (' ',15X,'CEILING',3X,'ROOF',5X,'ROOF',3X,'SHEATHING')

 WRITE (*,32)

 32 FORMAT (' LENGTH WIDTH HEIGHT OVERHANG PITCH THICKNESS')

 WRITE (*,33) DL,DW,DC,DOH,PITCH,D

 33 FORMAT (F6.2,1X,F6.2,2X,F6.2,4X,F6.2,2X,F6.2,3X,F6.2)

 WRITE (*,34)

 34 FORMAT (' MATERIAL PROPERTIES:')

 WRITE (*,35)

 35 FORMAT (1X,12X,'R-VALUES')

 WRITE (*,36)

 36 FORMAT (' CEILING ROOF SHEATHING ENDWALLS/SOFFITS WINDSHADE FA

 1CTOR SHINGLE ABSORPT.')

 WRITE (*,37) RC,RR,RS,RES,WSHADE,A

 37 FORMAT (' ',3(F6.2,' '),6X,F6.2,13X,F6.2,8X,F6.2)

 WRITE (*,38)

 38 FORMAT (' VENTILATION DATA:')

 WRITE (*,39)

 39 FORMAT (' HOUSE EXFILTRATION ATTIC VENTILATION')

 WRITE (*,40)

 40 FORMAT (' VOL. CHG. FRACTION NATURAL FAN')

 WRITE (*,41)

 41 FORMAT (' PER HOUR TO ATTIC (VOL./HR) (CFM)')

 WRITE (*,42) IH,F,IA,AF

 42 FORMAT (' ',2(F6.2,' '),' ',2(F6.2,' '))

 WRITE (*,43)

 43 FORMAT (' ENVIRONMENT:')

 WRITE (*,44)

 44 FORMAT (' INDOOR INDOOR MOISTURE CONTENT MOISTURE CON

 1TENT')

 WRITE (*,45)

 45 FORMAT (' TEMP. (F.) RH SHEATHING "A" SHEATHING "

 1B"')

 WRITE (*,46) T(7)-459.67,RHI,WMC(1),WMC(2)

 46 FORMAT (' ',2(' ',F6.2,' '),2(' ',F6.3,' '))

 WRITE (*,50)

 50 FORMAT ('01 CHANGE DIMENSIONS')

 WRITE (*,51)

 51 FORMAT (' 2 CHANGE MATERIAL PROPERTIES')

 WRITE (*,52)

 52 FORMAT (' 3 CHANGE VENTILATION DATA')

 WRITE (*,53)

 53 FORMAT (' 4 CHANGE ENVIRONMENT')

 WRITE (*,54)

 54 FORMAT (' 5 PROCEED')

 WRITE (*,55)

 55 FORMAT (' (ENTER CHOICE)')

 READ (*,*) K

 IF (K .NE. 1) GO TO 70

C

 WRITE (*,60) DL

 60 FORMAT (' LENGTH OF HOUSE (',F6.2,')')

 READ (*,*) ADL

 DL=ADL
27

 61 WRITE (*,62) DW

 62 FORMAT (' WIDTH OF HOUSE (',F6.2,')')

 READ (*,*) ADW

 DW=ADW

 WRITE (*,63) DC

 63 FORMAT (' HEIGHT OF CEILING (',F6.2,')')

 READ (*,*) ADC

 DC=ADC

 WRITE (*,64) DOH

 64 FORMAT (' WIDTH OF ROOF OVERHANG (',F6.2,')')

 READ (*,*) ADOH

 DOH=ADOH

 WRITE (*,65) PITCH

 65 FORMAT (' ROOF PITCH (RISE/RUN) (',F6.3,')')

 READ (*,*) APITCH

 PITCH=APITCH

 WRITE (*,66) D

 66 FORMAT (' SHEATHING THICKNESS (',F6.3,')')

 READ (*,*) AD

 D=AD

 GO TO 26

 70 IF (K .NE. 2) GO TO 80

 WRITE (*,71) RC

 71 FORMAT (' R-VALUE OF CEILING INSULATION (',F6.2,')')

 READ (*,*) ARC

 RC=ARC

 WRITE (*,72) RR

 72 FORMAT (' TOTAL R-VALUE OF ROOF (',F6.2,')')

 READ (*,*) ARR

 RR=ARR

 WRITE (*,75) RS

 75 FORMAT (' R-VALUE OF WOOD SHEATHING (',F6.2,')')

 READ (*,*) ARS

 RS=ARS

 WRITE (*,73) RES

 73 FORMAT (' R-VALUE OF ENDWALLS AND SOFFITS (',F6.2,')')

 READ (*,*) ARES

 RES=ARES

 WRITE (*,74) WSHADE

 74 FORMAT (' WINDSHADE FACTOR: 1=NO SHADE, 0=FULL SHADING (',

 1F6.2,')')

 READ (*,*) ASHADE

 WSHADE=ASHADE

 WRITE (*,76) A

 76 FORMAT (' ABSORBTIVITY EXTERIOR ROOF SURFACE (',F6.2,')')

 READ (*,*) A

 GO TO 26

 80 IF (K .NE. 3) GO TO 90

 WRITE (*,81) IH

 81 FORMAT (' HOUSE AIR EXCHANGE RATE (VOLUME CHANGES/HR) (',F6.2,')')

 READ (*,*) AIH

 IH=AIH

 WRITE (*,82) F

 82 FORMAT (' FRACTION OF HOUSE EXFILTRATION ENTERING ATTIC (',F6.2,')

 1')
28

 READ (*,*) AFA

 F=AFA

 85 WRITE (*,86) IA

 86 FORMAT (' NATURAL ATTIC VENTILATION (AIR CHANGES/HR) (',F6.2,')')

 READ (*,*) AIA

 IA=AIA

 WRITE (*,87) AF

 87 FORMAT (' ATTIC FAN CAPACITY (CFM) (',F6.2,')')

 READ (*,*) AAF

 AF=AAF

 GO TO 26

 90 IF (K .NE. 4) GO TO 100

 WRITE (*,91) T(7)-459.67

 91 FORMAT (' INDOOR TEMP., F. (',F6.2,')')

 READ (*,*) ATEM

 T(7)=ATEM+459.67

 WRITE (*,92) RHI

 92 FORMAT (' INDOOR RELATIVE HUMIDITY (',F6.2,')')

 READ (*,*) ARHI

 RHI=ARHI

 WRITE (*,93) WMC(1)

 93 FORMAT (' MOISTURE CONTENT OF SHEATHING "A" (',F6.3,')')

 READ (*,*) AWMCA

 WMC(1)=AWMCA

 WRITE (*,94) WMC(2)

 94 FORMAT (' MOISTURE CONTENT OF SHEATHING "B" (',F6.3,')')

 READ (*,*) AWMCB

 WMC(2)=AWMCB

 GO TO 26

 100 CONTINUE

C

 RETURN

 END

 SUBROUTINE HEAT

$DEBUG

C VERSION 7A

C Correction: 5/7/94. Corrected handling of windshade factor

C Date April 15, 1993: updated COMMON block

C Correction: March 2, 1993

C Cloud cover data in range 1 - 10

C Convective heat transfer routine for top and bottom of roof

C from Wilkes (1989). (for other sources see HCON subroutine)

C Infrared radiation routine on basis of Martin and Berdahl,

C Solar Energy vol 33, no 3/4, pp 321-336, 1984

C Revised August 4, 1992, Anton TenWolde

C Major revision August 17, 1992: incorporated SLOPE

 COMMON A,AC,AES,AF,AR,ALR,ALF,F,I1,I2,IA,IAF,IFLAG,IH,ITMOIST,

 1JJ,KK,L,NFAN,RC,RES,RHMAX,RR,RS,VA,VH,RHA,RHI,PITCH,VENT,PHI,

 2VLHA,VLHB,T(10),TARR(184,24),RAARR(184,24),WIND(184,24),WSHADE,

 3RBARR(184,24),SKYCOV(184,24),SNOCOV(184,24),TATT(184,24),

 4TRA(185,24),TRB(185,24),OUTDP(184,24),WMCA(185,24),

 5WMCB(185,24),RHAT(184,24),FAN(184,24),TEMPA(184,24),TEMPB(184,24),
29

 6TEMPOA(184,24),TEMPOB(184,24),AVEMC(1,24),D,IIILAT,SLOPE(6),JJJ

 REAL I1,I2,IA,IAF,IH,L

 DOUBLE PRECISION AK(6,6),Q(6),X(6)

 DIMENSION TOLD(6),TDEFLT(6)

C

C

C Initialize

 T(8)=TARR(JJ,KK)+(JJJ-1)*SLOPE(1)+459.67

C

 ITERAT=0

C

C If the depth of the snow layer is 1/2 inch or greater, determine

C the heat balance with nodes 4a and 4b set at 32 degrees F.

C Otherwise, go to line 1000 and perform the heat balance for

C all 6 nodes:

 IF (SNOCOV(JJ,KK) .LT. 0.5) GO TO 1000

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

C Snow-covered roof:

C

C Start iterative loop

 200 ITERAT=ITERAT+1

 IF (ITERAT.GE.10) THEN

 WRITE (6,10)

 10 FORMAT(' ERROR: TEMPERATURE CALCULATION DID NOT CONVERGE')

C Assign default temperatures (temperatures at end of first iteration)

 DO 201 MI=1,4

 T(MI)=TDEFLT(MI)

 201 CONTINUE

 GO TO 2000

 ENDIF

 FA=0.15426E-8*(T(2)*T(2)+T(3)*T(3))*(T(2)+T(3))

 FB=0.15426E-8*(T(2)*T(2)+T(4)*T(4))*(T(2)+T(4))

 DO 210 MI=1,4

 TOLD(MI)=T(MI)

 210 CONTINUE

C

C For first calculation set all heat transfer coeff. at 1.

 IF (ITERAT.EQ.1) THEN

 HRA=1.

 HRB=1.

 HF=1.

 ELSE

C Calculate heat transfer coefficients HRA, HRB, HF

 TS=T(3)-459.67

 TA=T(1)-459.67

 PPHI=PHI

 AL=ALR

 CALL HCON(TS,TA,PPHI,AL,2,0.0,HRA)

 TS=T(4)-459.67

 CALL HCON(TS,TA,PPHI,AL,2,0.0,HRB)

 TS=T(2)-459.67

 AL=ALF

 CALL HCON(TS,TA,0.,AL,1,0.0,HF)

 ENDIF

C
30

 AK(1,1)=-(AC*HF)-(F*IH*VH*.018)-(AR/2)*(HRA+HRB)-(AES/RES)

 A -(VENT*.018*VA)

 AK(1,2)=AC*HF

 AK(1,3)=.5*AR*HRA

 AK(1,4)=.5*AR*HRB

 AK(2,1)=-HF

 AK(2,2)=HF+FA/2+FB/2+1/RC

 AK(2,3)=-FA/2

 AK(2,4)=-FB/2

 AK(3,1)=HRA

 AK(3,2)=(AC/AR)*FA

 AK(3,3)=-HRA-(AC/AR)*FA-1/RR

 AK(3,4)=0.0D0

 AK(4,1)=HRB

 AK(4,2)=(AC/AR)*FB

 AK(4,3)=0.0D0

 AK(4,4)=-HRB-(AC/AR)*FB-1/RR

C

 X(1)=-0.018*F*IH*VH*T(7)-(T(8)*AES/RES)-(0.018*VENT*VA*T(8))

 X(2)=T(7)/RC

 X(3)=-VLHA-491.67/RR

 X(4)=-VLHB-491.67/RR

C

 DO 300 I=1,4

 300 Q(I)=X(I)

C Solve the simultaneous equations using the Gauss-Jordan method:

 500 DO 530 K=1,3

 DO 515 I=K+1,4

 IF (K .EQ. 1) GO TO 510

 DO 505 J=1,K-1

 505 AK(I,K)=AK(I,K)-AK(I,J)*AK(J,K)

 510 AK(I,K)=AK(I,K)/AK(K,K)

 515 CONTINUE

 DO 525 I=K+1,4

 DO 520 J=1,K

 520 AK(K+1,I)=AK(K+1,I)-AK(K+1,J)*AK(J,I)

 525 CONTINUE

 530 CONTINUE

 DO 540 I=2,4

 DO 535 K=1,I-1

 535 Q(I)=Q(I)-AK(I,K)*Q(K)

 540 CONTINUE

 DO 555 I=4,1,-1

 IF(I .EQ. 4) GO TO 550

 DO 545 K=I+1,4

 545 Q(I)=Q(I)-AK(I,K)*Q(K)

 550 Q(I)=Q(I)/AK(I,I)

 555 CONTINUE

 DO 560 I=1,4

 560 T(I)=Q(I)

C

C Set temperatures at nodes 4a and 4b to 32 degrees F (491.67 R):

 T(5)=491.67

 T(6)=491.67

C
31

C Set default temperatures

 IF (ITERAT.EQ.1) THEN

 DO 570 MI=1,4

 570 TDEFLT(MI)=T(MI)

 ENDIF

C

C Check for convergence in surface temperatures (change<0.1F)

 NCNVRG=0

 DO 600 MI=2,4

 IF (ABS(TOLD(MI)-T(MI)).GT.0.1) NCNVRG=1

 600 CONTINUE

 IF (NCNVRG.EQ.1) GOTO 200

 GO TO 2000

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

C Without snow on the roof:

C

C Calculate contribution of infrared radiation to sky

 1000 TD=OUTDP(JJ,KK)+(JJJ-1)*SLOPE(5)

 CLDAMT=SKYCOV(JJ,KK)+(JJJ-1)*SLOPE(4)

 TO=T(8)

 CALL INFRED(TO,TD,KK,CLDAMT,QINFRD,HINFRD)

C

C Convert windspeed from knots to feet per hour,

C and correct for windshading

 V=6077.28*WSHADE*(WIND(JJ,KK)+(JJJ-1)*SLOPE(6))

C

C Start iterative loop

 1200 ITERAT=ITERAT+1

 IF (ITERAT.GE.10) THEN

C Write error message

 WRITE (6,10)

C Assign default temperatures (temperatures at end of first iteration) 

 DO 1201 MI=1,6

 T(MI)=TDEFLT(MI)

 1201 CONTINUE

 GO TO 2000

 ENDIF

C

 FA=0.15426E-8*(T(2)*T(2)+T(3)*T(3))*(T(2)+T(3))

 FB=0.15426E-8*(T(2)*T(2)+T(4)*T(4))*(T(2)+T(4))

 DO 1010 MI=1,6

 TOLD(MI)=T(MI)

 1010 CONTINUE

C

C For first calculation set all heat transfer coeff. at 1.

 IF (ITERAT.EQ.1) THEN

 HRA=1.

 HRB=1.

 HF=1.

 HOA=1.

 HOB=1.

 ELSE

C Calculate heat transfer coefficients HRA, HRB, HF, HOA, HOB

 TS=T(3)-459.67

 TA=T(1)-459.67
32

 PPHI=PHI

 AL=ALR

 CALL HCON(TS,TA,PPHI,AL,2,0.0,HRA)

 TS=T(4)-459.67

 CALL HCON(TS,TA,PPHI,AL,2,0.0,HRB)

 TS=T(2)-459.67

 AL=ALF

 CALL HCON(TS,TA,0.,AL,1,0.0,HF)

 TA=T(8)-459.67

 TS=T(5)-459.67

 AL=ALR

 VV=V

 CALL HCON(TS,TA,PPHI,AL,1,VV,HOA)

 TS=T(6)-459.67

 CALL HCON(TS,TA,PPHI,AL,1,VV,HOB)

 ENDIF

C

 AK(1,1)=-(AC*HF)-(F*IH*VH*0.018)-(AR/2.0D0)*(HRA+HRB)-(AES/RES)

 A -(VENT*0.018D0*VA)

 AK(1,2)=AC*HF

 AK(1,3)=0.5D0*AR*HRA

 AK(1,4)=0.5D0*AR*HRB

 AK(1,5)=0.0D0

 AK(1,6)=0.0D0

 AK(2,1)=-HF

 AK(2,2)=HF+FA/2.0D0+FB/2.0D0+1.0D0/RC

 AK(2,3)=-FA/2.0D0

 AK(2,4)=-FB/2.0D0

 AK(2,5)=0.0D0

 AK(2,6)=0.0D0

 AK(3,1)=HRA

 AK(3,2)=(AC/AR)*FA

 AK(3,3)=-HRA-(AC/AR)*FA-1.0D0/RR

 AK(3,4)=0.0D0

 AK(3,5)=1.0D0/RR

 AK(3,6)=0.0D0

 AK(4,1)=HRB

 AK(4,2)=(AC/AR)*FB

 AK(4,3)=0.0D0

 AK(4,4)=-HRB-(AC/AR)*FB-1.0D0/RR

 AK(4,5)=0.0D0

 AK(4,6)=1.0D0/RR

 AK(5,1)=0.0D0

 AK(5,2)=0.0D0

 AK(5,3)=1.0D0/RR

 AK(5,4)=0.0D0

 AK(5,5)=-1.0D0/RR-HOA-HINFRD

 AK(5,6)=0.0D0

 AK(6,1)=0.0D0

 AK(6,2)=0.0D0

 AK(6,3)=0.0D0

 AK(6,4)=1.0D0/RR

 AK(6,5)=0.0D0

 AK(6,6)=-1.0D0/RR-HOB-HINFRD

 X(1)=-0.018D0*F*IH*VH*T(7)-(T(8)*AES/RES)-(0.018D0*VENT*VA*T(8))
33

 X(2)=T(7)/RC

 X(3)=-VLHA

 X(4)=-VLHB

 X(5)=-A*(RAARR(JJ,KK)+(JJJ-1)*SLOPE(2))+QINFRD-T(8)*(HOA+HINFRD)

 X(6)=-A*(RBARR(JJ,KK)+(JJJ-1)*SLOPE(3))+QINFRD-T(8)*(HOB+HINFRD)

C

 DO 1300 I=1,6

 1300 Q(I)=X(I)

C

C Solve the simultaneous equations using the Gauss-Jordan method:

 1500 DO 1530 K=1,5

 DO 1515 I=K+1,6

 IF (K .EQ. 1) GO TO 1510

 DO 1505 J=1,K-1

 1505 AK(I,K)=AK(I,K)-AK(I,J)*AK(J,K)

 1510 AK(I,K)=AK(I,K)/AK(K,K)

 1515 CONTINUE

 DO 1525 I=K+1,6

 DO 1520 J=1,K

 1520 AK(K+1,I)=AK(K+1,I)-AK(K+1,J)*AK(J,I)

 1525 CONTINUE

 1530 CONTINUE

 DO 1540 I=2,6

 DO 1535 K=1,I-1

 1535 Q(I)=Q(I)-AK(I,K)*Q(K)

 1540 CONTINUE

 DO 1555 I=6,1,-1

 IF (I .EQ. 6) GO TO 1550

 DO 1545 K=I+1,6

 1545 Q(I)=Q(I)-AK(I,K)*Q(K)

 1550 Q(I)=Q(I)/AK(I,I)

 1555 CONTINUE

 DO 1560 I=1,6

 1560 T(I)=Q(I)

C

C

C Set default temperatures

 IF (ITERAT.EQ.1) THEN

 DO 1570 MI=1,6

 1570 TDEFLT(MI)=T(MI)

 ENDIF

C

C Check for convergence in surface temperatures (change<0.1F)

 NCNVRG=0

 DO 1600 MI=2,6

 IF (ABS(TOLD(MI)-T(MI)).GT.0.1) NCNVRG=1

 1600 CONTINUE

 IF (NCNVRG.EQ.1) GOTO 1200

 GO TO 2000

C

C End of iterative loops

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

 2000 TATT(JJ,KK)=T(1)-459.67

 RETURN

 END
34

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

 SUBROUTINE INFRED(TO,TD,IHR,CLDAMT,QINFRD,HINFRD)

C

C Calculates infrared radiation to and from sky. Part (QINFRD) is

C subtracted from solar radiation; the rest is handled with a

C modification of heat transfer coefficient HO (HINFRD)

C Base of clouds is assumed at 2 km.

C

 EPSRF=0.9

 TDD=(TD-32)/1.8

C Clear sky emissivity

 EPS0=0.711+0.0056*TDD+0.000073*TDD*TDD+0.013*COS(3.14159*IHR/12)

C With clouds:

 EPS=EPS0+0.784*(1-EPS0)*CLDAMT/10

 HINFRD=4*0.17117E-8*EPSRF*TO*TO*TO

 QINFRD=HINFRD*(1-EPS)*TO*0.25

 RETURN

 END

 SUBROUTINE HCON(TS,TA,PHI,AL,IFLAG,V,HC)

$DEBUG

C Version 7

C Corrected May 7, 1994, Anton TenWolde

C Latest update August 7, 1992

C

C Adapted from HCON routine in "Model for roof thermal performance",

C K.E. Wilkes, ORNL/CON-274, Oak Ridge Nat.Lab., 1989.

C Limited to flat surface and surfaces with tilt over 2 degrees

C TS = surface temperature, F

C TA = air temperature, F

C PHI = tilt angle, degrees (horizontal=0)

C AL = characteristic length ,(length + width)/2, ft

C IFLAG = 1 for surface facing up

C IFLAG = 2 for surface facing down

C V = air speed, feet per hour

C HCF = forced convection coefficient

C HCN = natural convection coefficient

C HC = total convection coefficient

C

 REAL NUS,K,MU,NU

 NERROR=0

 DT=TS-TA

 ABSPHI=ABS(PHI)

 IF (IFLAG.EQ.2) DT=-DT

C Calculate film temperature

 IF (ABSPHI.LT.0.001) THEN

 TF=(TS+TA)/2

 TF1=TF

 ELSE

 TF=TS-0.25*(TS-TA)

 TF1=TA+0.25*(TS-TA)

 ENDIF
35

 TK=(TF+459.67)/1.8

C Check for negative temperatures

 IF (TK.LE.0) NERROR=1

 IF (TF1.LT.-459.) NERROR=2

 IF (NERROR.NE.0) GOTO 100

C

 PI=3.14159265

C Thermal conductivity of air, Btu/hr-ft-F

 K=0.6325E-5*SQRT(TK)/(1+(245.4*10**(-12/TK))/TK)*241.77

C Dynamic viscosity of air, lb/hr-ft

 MU=(145.8*TK*SQRT(TK)/(TK+110.4))*241.9E-7

C Prandtl number

 PR=0.788-2.631E-4*TK

 IF (PR.LE.0.) NERROR=3

C In case of errors:

 IF (NERROR.NE.0) THEN

100 WRITE (6,10) NERROR

10 FORMAT(' ERROR ',I1,' IN CONVECTIVE HEAT TRANSFER COEFFICIENT')

 HC=1.

 RETURN

 ENDIF

C

C Volume expansion coefficient of air, perfect gas

 BETA=1/(TF1+459.67)

C Density of air, lbs/ft3

 RHO=22.0493/TK

C Kinematic viscosity

 NU=MU/RHO

C Specific heat of air, Btu/lb-ft

 CP=(3.4763+1.066E-4*TK)*0.068559

C Raleigh number

 RA=4.16975E8*BETA*RHO*CP*ABS(DT)*(AL**3)/(NU*K)

C

C Calculate convection coefficient (w/o forced convection)

C Branch to different correlations depending on surface orientation

C Roof is considered flat if tilt is less than 2 degrees

 IF (ABSPHI.LT.2.0) THEN

C Horizontal surface (0-2 degrees)

 IF (DT.GE.0.0) THEN

C Heat flow upwards

 NUS=0.15*RA**0.33333

 IF (RA.LT.8.0E6) NUS=0.54*RA**0.25

 ELSE

C Heat flow down

 NUS=0.85*RA**0.2

 ENDIF

 ELSE

C Tilted surface (over 2 degrees)

 IF (DT.GT.0.0) THEN

C Heat flow upwards

C Calculate critical Grasshof number

 GRC=10**(PHI/(1.1870+0.0870*PHI))

 IF (ABSPHI.LT.15.) GRC=1.0E6

 IF (ABSPHI.GT.75.) GRC=5.0E9

 GR=RA/PR
36

 NUS=0.56*(RA*COS((90.-PHI)*PI/180))**0.25

 IF (GR.GT.GRC) NUS=0.14*(RA**0.33333-(GRC*PR)**0.33333)

 A +NUS

 ELSE

C Heat flow down

 NUS=0.56*(RA*COS((90.-PHI)*PI/180))**0.25

 ENDIF

 ENDIF

 HCN=NUS*K/AL

C Set a minimum of 0.1 for HCN to avoid swings in

C iterative calculations

 HCN=MAX(0.1,HCN)

 IF (V.LE.1.0E-12) THEN

 HC=HCN

 RETURN

 ENDIF

C

C Calculate forced convection coefficient

 RE=V*AL/NU

 IF (RE.LT.5.0E5) THEN

 NUS=0.664*(PR**0.33333)*SQRT(RE)

 ELSE

 NUS=(PR**0.33333)*(0.037*(RE**0.8)-850)

 ENDIF

 HCF=NUS*K/AL

C Combine natural and forced convection coefficients

 HC=(HCF**3+HCN**3)**0.3333333

 RETURN

 END

 SUBROUTINE MOIST(WMC,WEXTRA)

$DEBUG

C VERSION 7A

C April, 1993

C Works with outdoor dewpoint temperature and cloud cover data 0-10

C

C This version of the moisture balance subroutine utilizes a three-layer

C system to represent the plywood sheathing. The first and second layers

C make up the first ply (1/8") facing the interior of the attic.

C The third layer makes up the remaining thickness of the plywood.

C

C History:

C Revised August 4, 1992 by Anton TenWolde

C Revised August 8, 1992 by Anton TenWolde: introduced surplus moisture

C storage and carry over to next time iteration (WEXTRA).

C Revised August 11, 1992: simplified diffusion coefficient;

C [reintroduced liquid transfer (proportional to delta MC)](removed 8/18)

C Revised August 12, 1992: Limited moisture loss from sheathing to

C equilibrium MC with attic air when RH attic=100%;

C Adjusted definition of WEXTRA

C Revised August 13, 1992: added messages when equilibrium limits are

C invoked

C Revised August 14, 1992: added temperature protection for PSAT
37

C August 15, 1992: added MC protection for PSAT

C Major revision August 17, 1992: incorporated SLOPE, synchronized timing

C of output reports

C August 18, 1992: tightened latent heat convergance criterion to 5%;

C removed liquid transfer; fixed error in WEXTRA calculation

C

C February 28, 1993: Created version 7: 3-layered sheathing, set TLAYER(1)

C and TLAYER(3) at surface temperature, changed definition of MCAVG(5) and

C MCAVG(6)

C April 15, 1993: incorporated variable sheathing thickness and R-value

C

 DOUBLE PRECISION WMC(6),WRA,WRB,WW(4),ULAYER(6),WDA,WDB,SATMC,

 1WHR,WEXTRA,WRAVVL,WRBVVL,WRMIN

 COMMON A,AC,AES,AF,AR,ALR,ALF,F,I1,I2,IA,IAF,IFLAG,IH,ITMOIST,

 1JJ,KK,L,NFAN,RC,RES,RHMAX,RR,RS,VA,VH,RHA,RHI,PITCH,VENT,PHI,

 2VLHA,VLHB,T(10),TARR(184,24),RAARR(184,24),WIND(184,24),WSHADE,

 3RBARR(184,24),SKYCOV(184,24),SNOCOV(184,24),TATT(184,24), 

 4TRA(185,24),TRB(185,24),OUTDP(184,24),WMCA(185,24),

 5WMCB(185,24),RHAT(184,24),FAN(184,24),TEMPA(184,24),TEMPB(184,24), 

 6TEMPOA(184,24),TEMPOB(184,24),AVEMC(1,24),D,IIILAT,SLOPE(6),JJJ

 REAL I1,I2,IA,IAF,IH,L,MCAV(4)

 DIMENSION TLAYER(6),DW(4)

C 

 LERROR=0

C Set the surface diffusion coefficient: 

 HDR=1.1

C 

C Read the outdoor dewpoint temperature (F) and convert to R: 

 DPO=OUTDP(JJ,KK)+(JJJ-1)*SLOPE(5)+459.67

C 

C Calculate humidity ratios for indoor air, outdoor air, surface of 

C sheathing "A", and surface of sheathing "B", respectively: 

 HRI=AHR(T(7),RHI) 

 HRO=AHR(DPO,1.0) 

10 WDA=WMC(1)

 WDB=WMC(2)

 HRRA=WHR(T(3),WDA) 

 HRRB=WHR(T(4),WDB) 

C

C

C Determine the temperature of the interface between the 1/8" bottom ply

C of the sheathing and the remaining mass of each sheathing:

C NOTE: The R-value of the shingles = R roof - R sheathing

C E.g.: R roof = 1.04, R sheathing = 0.6, then R shingles = 0.44

C

 TINTA=(0.125/D)*(RS/RR)*(T(5)-T(3))+T(3)

 TINTB=(0.125/D)*(RS/RR)*(T(6)-T(4))+T(4)

C

C Determine temperature at each node: 1st layer is at surface temperature,

C 2nd and 3rd layers are at temperature in middle of layer

 TLAYER(1)=T(3)

 TLAYER(2)=T(4)

 TLAYER(3)=0.75*(TINTA-T(3))+T(3)

 TLAYER(4)=0.75*(TINTB-T(4))+T(4)

 TOUTA=T(3)+RS*(T(5)-T(3))/RR
38

 TOUTB=T(4)+RS*(T(6)-T(4))/RR

 TLAYER(5)=(TOUTA+TINTA)/2

 TLAYER(6)=(TOUTB+TINTB)/2

C

C 

C Calculate the humidity ratio at dewpoint for sheathing surfaces: 

 DPTSA=PSAT(T(3))/48.1 

 DPTSB=PSAT(T(4))/48.1 

C 

C Set humidity ratios for each sheathing surface to the minimum of 

C the humidity ratios based on current MC and humidity ratios based 

C on dewpoint temperatures: 

C 

 HRRAAIR=MIN(HRRA,DPTSA) 

 HRRBAIR=MIN(HRRB,DPTSB) 

C 

C 

C Calculate the humidity ratio of the attic air (revised 8/7/92): 

 HRA=(WEXTRA+0.075*F*IH*VH*HRI+0.075*VENT*VA*HRO+AR*HDR*HRRAAIR*

 a0.5+AR*HDR*HRRBAIR*0.5)/(0.075*F*IH*VH+0.075*VENT*VA+AR*HDR)

C 

C Calculate the saturation humidity ratio of the attic air:

 HRASAT=PSAT(T(1))/48.1

C Convert the humidity ratio to relative humidity: 

 RHA=HRA/HRASAT 

C

C If attic RH is less than 100%:

 IF (RHA .LE. 1.0) THEN

WRAVVL=HDR*(HRA-HRRAAIR)

WRBVVL=HDR*(HRA-HRRBAIR)

MRH100=0

 ELSE

C Attic RH=100%:

HRA=HRASAT

RHA=1

MRH100=1

WRAVVL=HDR*(HRA-HRRAAIR)

WRBVVL=HDR*(HRA-HRRBAIR)

C The next section only applies if moisture is driven from sheathing

C surface A.

IF (WRAVVL.LT.0.0D0) THEN

C Calculate equilibrium MC at surface. This is minimum MC of

C surface and moisture loss is allowed down to this minimum MC

C Note: WRAVVL is negative: WRMIN is a maximum negative flow

 Z=HRA/(PSAT(T(3))*15.8514)

 WMCEQL=0.029+169.6*Z-269200*Z**2+2.2E08*Z**3

 WRMIN=0.08125*6*(WMCEQL-WMC(1))

 IF (WRMIN.GE.0.0D0) WRMIN=0.0D0

 IF (WRAVVL.LT.WRMIN) THEN

 WRAVVL=WRMIN

 WRITE(6,28)

28 FORMAT(' WRAVVL has been limited')

 ENDIF

ENDIF

C Same for B side:
39

IF (WRBVVL.LT.0.0D0) THEN

 Z=HRA/(PSAT(T(4))*15.8514)

 WMCEQL=0.029+169.6*Z-269200*Z**2+2.2E08*Z**3

 WRMIN=0.08125*6*(WMCEQL-WMC(1))

 IF (WRMIN.GE.0.0D0) WRMIN=0.0D0

 IF (WRBVVL.LT.WRMIN) THEN

 WRBVVL=WRMIN

 WRITE(6,29)

29 FORMAT(' WRBVVL has been limited')

 ENDIF

ENDIF

 ENDIF

C

 WRA=0.5*AR*WRAVVL/6

 WRB=0.5*AR*WRBVVL/6

C

C Bypass latent heat section if IIILAT=0 or LERROR=1

 IF (IIILAT.EQ.0.OR.LERROR.EQ.1) GOTO 53

C 

C Calculate the heat flux attributed to the latent heat of 

C vaporization for each sheathing surface: 

 VVLHA=WRAVVL*1050

 VVLHB=WRBVVL*1050

C 

C If this is the first iteration of the moisture subprogram,

C prepare for next iteration (latent heat).

 IF (ITMOIST.EQ. 0) GO TO 200 

C 

C Compare the new latent heat values with the old. If the new value 

C is very small or within 5% of the old or the difference is very small,

C set the flag and continue: 

 IF (ABS(VVLHA-VLHA).LT.1.0E-20.OR.ABS(VVLHA).LT.1.0E-20) GOTO 50

 IF (ABS(VVLHA-VLHA).GT.ABS(VLHA*.05)) GOTO 200

 50 IF (ABS(VVLHB-VLHB).LT.1.0E-20.OR.ABS(VVLHB).LT.1.0E-20) GOTO 53

 IF (ABS(VVLHB-VLHB).GT.ABS(VLHB*.05)) GOTO 200

 53 IFLAG=1

 ITMOIST=0

C

C If RH attic =100%, calculate amount of "surplus" moisture, lbs/hr:

C Note: surplus is added to moisture balance in next time iteration.

 WEXTRA=0.0D0

 IF (MRH100.EQ.1) WEXTRA=0.075*(F*IH*VH*(HRI-HRASAT)-

 a VENT*VA*(HRASAT-HRO))-6*WRA-6*WRB

C

 IF (WEXTRA.GT.0.0D0) THEN

 WRITE (6,31) WEXTRA

 ELSE

 WEXTRA=0.0D0

 ENDIF

 31 FORMAT (' ','EXCESS MOISTURE (LBS/HR): ',G9.2) 

C

C 

C These next steps determine wood humidity ratios for the sheathing

C layers, based on current MC and temps: 

 SATMC=.3 
40

 DO 100 J=1,5,2

C 

C If the wood moisture content is above saturation, humidity ratio is

C the same as at saturation (liquid transfer ignored).

 IF (WMC(J)-SATMC) 55,55,60 

 55 WDA=WMC(J)

 ULAYER(J)=WHR(TLAYER(J),WDA) 

 GO TO 62

 60 ULAYER(J)=WHR(TLAYER(J),SATMC)

 62 CONTINUE

 IF (WMC(J+1)-SATMC) 65,65,70

 65 WDB=WMC(J+1)

 ULAYER(J+1)=WHR(TLAYER(J+1),WDB) 

 GO TO 72

 70 ULAYER(J+1)=WHR(TLAYER(J+1),SATMC) 

 72 CONTINUE

C

 100 CONTINUE 

C 

C Determine diffusion coefficients for moisture transfer between

C layers, and calculate amount of moisture transferred:

 DO 110 III=1,4

C

C First, calculate average MC between adjacent sheathing layers:

 IF (WMC(III).LT.1000.0D0.AND.WMC(III+2).LT.1000.0D0) THEN

MCAV(III)=WMC(III)/2+WMC(III+2)/2

 ELSE

MCAV(III)=SATMC

 ENDIF

C 

C If the average MC is above saturation, calculate the interior

C diffusion coefficient at saturation; also check for negative MC.

 IF (MCAV(III).GT.SATMC) MCAV(III)=SATMC

 IF (MCAV(III).LT.0) MCAV(III)=1.0E-10

C Calculate transport coefficients

 DW(III)=EXP(-9.37455+19.97645*MCAV(III))

 110 CONTINUE

C Calculate the amount of moisture transferred between first two

C sheathing layers (every 1/6 hour):

 DO 130 III=1,2

 WW(III)=(DW(III)*AR/12)*(ULAYER(III)-ULAYER(III+2))/((1.0/16)/12)

 130 CONTINUE

C Same for layers 5 and 6

 DO 170 III=3,4

 WW(III)=(DW(III)*AR/12)*(ULAYER(III)-ULAYER(III+2))

 a/(((D-.0625)/2)/12)

 170 CONTINUE

C

C Given the following mass for sheathing surface and interior wood 

C layers, calculate the new MC for each portion: 

 WMASS=(AR/2)*((.125/2)/12)*31.2 

 WMC(1)=WMC(1)+(WRA-WW(1))/WMASS

 WMC(2)=WMC(2)+(WRB-WW(2))/WMASS

 WMC(3)=WMC(3)+(WW(1)-WW(3))/WMASS

 WMC(4)=WMC(4)+(WW(2)-WW(4))/WMASS
41

 WMC(5)=WMC(5)+(WW(3))/((AR/2)*((D-.125)/12)*31.2)

 WMC(6)=WMC(6)+(WW(4))/((AR/2)*((D-.125)/12)*31.2)

C

C Assign RH and temperatures at start of hour to output array:

 IF (JJJ.EQ.1) THEN

RHAT(JJ,KK)=RHA

TEMPA(JJ,KK)=T(3)-459.67

TEMPB(JJ,KK)=T(4)-459.67

TEMPOA(JJ,KK)=TOUTA-459.67

TEMPOB(JJ,KK)=TOUTB-459.67

 ENDIF

C

C Calculate average MC for the sheathing and at end of hour.

C Assign MC values at end of hour to next hour out put array:

 IF (JJJ.EQ.6) THEN

NXTHR=KK+1

NXTDAY=JJ

IF (KK.EQ.24) THEN

 NXTHR=1

 NXTDAY=JJ+1

ENDIF

TRA(NXTDAY,NXTHR)=((WMC(1)+WMC(3))/2)*(.125/D)

 A +WMC(5)*((D-.125)/D)

TRB(NXTDAY,NXTHR)=((WMC(2)+WMC(4))/2)*(.125/D)

 A +WMC(6)*((D-.125)/D)

C Assign current calculated MC values to output arrays:

WMCA(NXTDAY,NXTHR)=WMC(1)

WMCB(NXTDAY,NXTHR)=WMC(2)

 ENDIF

C

 GO TO 210

C

C Preparation for next latent heat iteration:

 200 ITMOIST=ITMOIST+1

 IF (ITMOIST.GE.10) THEN

 WRITE (6,202)

 202 FORMAT (' ERROR: LATENT HEAT ITERATIONS DID NOT CONVERGE')

C Substitute temperatures and mc's without latent heat

 LERROR=1

 VVLHA=0.

 VVLHB=0.

 ENDIF

 VLHA=VVLHA

 VLHB=VVLHB

 CALL HEAT

 GO TO 10

C

C

 210 CONTINUE

C

 RETURN

 END

C

C

C THE FOLLOWING FUNCTION RETURNS AN AIR HUMIDITY RATIO WHEN GIVEN
42

C AIR TEMP. AND RELATIVE HUMIDITY

C

 FUNCTION AHR(TEM,RH)

 AHR=RH*PSAT(TEM)/48.1

 RETURN

 END

C

C

C THE FOLLOWING FUNCTION RETURNS THE SATURATED VAPOR PRESSURE WHEN

C GIVEN AIR TEMP (vapor pressure in inch Hg., temperature in degree R) 

C

 FUNCTION PSAT(TEM)

C Set limits to prevent overflow

 IF (TEM.LT.70) THEN

EW=0.

 ELSE IF (TEM.GT.671.69) THEN

EW=1013.246

C Note: = saturation pressure at steam point, in mb

 ELSE

RATIO=671.69/TEM

EWEXP=-7.90298*(RATIO-1)+5.02808*LOG10(RATIO)-

 A 0.13816E-06*(10**(11.344*(1-1/RATIO))-1)+

 B .81328E-02*(10**(-3.49149*(RATIO-1))-1)+3.005715

C Check to prevent overflow

IF (EWEXP.LT.-30) EWEXP=-30

IF (EWEXP.GT.30) EWEXP=30

EW=10**EWEXP

 END IF

 PSAT=EW*.02953

 RETURN

 END

C

C

C THE FOLLOWING FUNCTION RETURNS A WOOD HUMIDITY RATIO WHEN GIVEN

C WOOD TEMPERATURE (degree R) AND MOISTURE CONTENT (ratio)

C

 FUNCTION WHR(TEM,XMC)

 DOUBLE PRECISION WHR,XMC,YMC,OMEGA50

 IF (XMC.GT.0.0D0) THEN

 YMC=XMC

 IF (YMC.GT.1.0D0) YMC=1.0D0

C Calculate humidity ratio according to Cleary's formula at 50 F

C Then adjust for temperature

C (EXP(50F/28.6)=5.744551, saturation pressure @ 50F = 0.3624 in Hg)

 OMEGA50=5.744551*(-.49E-03+.0172*YMC-.06*YMC**2+.076*YMC**3)

 WHR=OMEGA50*PSAT(TEM)/0.3624

 IF (WHR.LT.0.0D0) WHR=0.0D0

 ELSE

 WHR=0.0D0

 ENDIF

 RETURN

 END
43

Appendix C—SUNDAT Program

The SUNDAT program can be used to calculate solar radiation on each roof surface. The program calculates the ratio of radiation on a surface with a specified tilt and orientation and on a horizontal surface at a given location and hour. The calculation is

performed for a cloudless sky, apportioning direct and diffuse radiation based on chapter 27 of the ASHRAE Handbook of Fundamentals (1989). This ratio is then applied to measured solar radiation data (horizontal surface) to obtain solar radiation incident on the specified tilted surface. SUNDAT produces solar radiation input files for the main program. In the main program,

radiation absorption is obtained by multiplying incident radiation by the roof's solar absorptance. Inputs to the program are

measured solar radiation on horizontal surface (Langley units), latitude, longitude, time zone, roof tilt angle (degrees from horizontal), and roof orientation. To generate the solar data for one building, SUNDAT has to be run for each roof surface (two

orientations), using the same solar radiation input data.

A listing of the SUNDAT program follows below:

C Program SUNDAT

C Version 2, Date May 12, 1994;

C revised July 19, 1994: added zero divide protection if S(13)=0

C revised July 15, 1994: added zeroing of output array (initialization)

C revised June 7, 1994: fixed old mistake in program

C revised June 6, 1994 (changed diffuse radiation on surface calculations)

C Based on SUN4 program (does not use cloud cover data)

C Calculates solar radiation on roof surface from solar radiation data

C for horizontal surface (in Btu/h.ft2)

C Use only for tilted roof surfaces

 REAL LATD,LONG,MERID,LOND

 INTEGER FDAY,LDAY,LCOUNT

 CHARACTER*12 FILENM

 CHARACTER*12 FNAME

 DIMENSION A0(5),A1(5),A2(5),A3(5),B1(5),B2(5),B3(5),SRAD(184,24),

 1RADDAT(184,24)

 DOUBLE PRECISION S(35)

 DATA A0/0.302,-0.0002,368.44,0.1717,0.0905/,

 2 A1/-22.93,0.4197,24.52,-0.0344,-0.0410/,

 3 A2/-0.229,-3.2265,-1.14,0.0032,0.0073/,

 4 A3/-0.243,-0.0903,-1.09,0.0024,0.0015/,

 5 B1/3.851,-7.351,0.58,-0.0043,-0.0034/,

 6 B2/0.002,-9.3912,-0.18,0.,0.0004/,

 7 B3/-0.055,-0.3361,0.28,-0.0008,-0.0006/

C S(1) = latitude, degrees (+ North, - South)

C S(2) = longitude, degrees (+ West, - East)

C S(3) = time zone number, standard time

C Atlantic 4

C Eastern 5

C Central 6

C Mountain 7

C Pacific 8

C S(4) = days (from start of year)

C S(5) = hour, after midnight

C S(6) = daylight saving time indicator (set at 0)

C S(7) = ground reflectivity

C S(8) = clearness number

C S(9) = azimuth angle, degrees from south

C S(10) = tilt angle, degrees from horizon

C S(11) = sun rise time

C S(12) = sunset time

C S(13) = cosine of angle between ray and normal (horiz. surface)
44

C S(14) = cosine of angle between ray and West

C S(15) = cosine of angle between ray and South

C S(16) = direction cosines normal to surface

C S(17),S(18) = angles

C S(19) = cosine of incidence angle

C S(20) = solar altitude angle

C S(21) = solar azimuth angle

C S(22) = diffuse sky radiation on horizontal surface

C S(23) = diffuse ground reflected radiation

C S(24) = direct normal radiation

C S(25) = total solar radiation intensity on surface

C S(26) = diffuse sky radiation intensity, for roof tilts >45 degree

C S(27) = ground reflected diffuse radiation intensity

C S(28) = sun declination angle, degrees

C S(29) = equation of time, hours

C S(30),S(31),S(32) = solar factors

C S(33) = cloud cover modifier

C S(34) = intensity of direct solar radiation on surface

C S(35) = hour angle, degrees

 S(6)=0

 S(7)=0.2

 S(8)=1.

 S(33)=1.

 PI=3.1415927

C

 WRITE (*,2)

 2 FORMAT (' ENTER DAYS OF YEAR TO BE INCLUDED IN FILE')

 WRITE (*,3)

 3 FORMAT (' FIRST DAY '\)

 READ (*,*) FDAY

 WRITE (*,4)

 4 FORMAT (' LAST DAY '\)

 READ (*,*) LDAY

 WRITE (*,15)

 15 FORMAT (' LATITUDE (DEGREES, NORTH = +, SOUTH = -)'\)

 READ (*,*) S(1)

 WRITE (*,16)

 16 FORMAT (' LONGITUDE (DEGREES, + = WEST, - = EAST)'\)

 READ (*,*) S(2)

 WRITE (*,17)

 17 FORMAT (' TIMEZONE (ATL=4, EAST=5, CENTR=6, MOUNT=7,PAC=7)'\)

 READ (*,*) S(3)

C

 IDAYS=LDAY-FDAY+1

C

C Read in solar radiation data

 WRITE (*,8)

 8 FORMAT (' ENTER NAME OF FILE CONTAINING SOLAR RADIATION DATA '\)

 READ (*,7) FNAME

 OPEN (7,FILE=FNAME,STATUS='OLD')

 READ (7,*) ((RADDAT(IROW,ICOL), ICOL=1,24), IROW=1,IDAYS)

 CLOSE (7)

C

 1000 WRITE (*,1)

 1 FORMAT (' ENTER NAME OF FILE TO BE CREATED '\)
45

 READ (*,7) FILENM

C

 WRITE (*,5)

 5 FORMAT (' ENTER ROOF TILT ANGLE (DEGREES FROM HORIZONTAL) '\)

 READ (*,*) S(10)

 WRITE (*,6)

 6 FORMAT (' ENTER ROOF SURFACE AZIMUTH ANGLE (DEGREES FROM SOUTH '\

 1' + = WEST, - = EAST)'\)

 READ (*,*) S(9)

C

 7 FORMAT (A12)

C

C Initialize output array (added on 7/15/94)

 DO 100 MM=1,184

 DO 100 MMM=1,24

 SRAD(MM,MMM)=0.

100 CONTINUE 

C

 LCOUNT=0

C

 DO 1300 JJ=FDAY,LDAY

 S(4)=JJ

 WRITE (*,200) JJ

200 FORMAT(' DAY: ',I3)

 LCOUNT=LCOUNT+1

 DO 1200 KK=1,24

 S(5)=KK-1.

C

 X=2*PI/366.*S(4)

 C1=COS(X)

 C2=COS(2*X)

 C3=COS(3*X)

 S1=SIN(X)

 S2=SIN(2*X)

 S3=SIN(3*X)

 DO 10 K=1,5

 KS=(K-1)+28

 10 S(KS)=A0(K)+A1(K)*C1+A2(K)*C2+A3(K)*C3+B1(K)*S1+B2(K)*S2+B3(K)*S3

 S(29)=S(29)/60.

 LATD=S(1)

 LONG=S(2)

 MERID=15*S(3)

 LOND=LONG-MERID

 Y=S(28)*PI/180.

 YY=LATD*PI/180.

 HP=-TAN(Y)*TAN(YY)

 TR=12/PI*ACOS(HP)

 S(11)=(12-TR)-S(29)+LOND/15.

 S(12)=24.-S(11)

 H=15*(S(5)-12+S(3)+S(29)-S(6))-S(2)

 S(35)=H

 S13=SIN(YY)*SIN(Y)+COS(YY)*COS(Y)*COS(H*PI/180.)

 S(13)=S13

 HP1=180.*ACOS(HP)/PI

 X1=ABS(HP1)
46

 X2=ABS(H)

 IF (X1.LT.X2) THEN

 DO 140 J=14,27

 140 S(J)=0.

 S(34)=0.

 GO TO 150

 ENDIF

C

 20 S(14)=COS(Y)*SIN(H*PI/180.)

 COSSSQ=1.-S(13)*S(13)-S(14)*S(14)

 IF (COSSSQ.LE.0) THEN

 S(15)=0

 ELSE

 S(15)=SQRT(COSSSQ)

 ENDIF

 STEST=S(15)

 STEST1=COS(H*PI/180.)-TAN(Y)/TAN(YY)

 IF (STEST1) 40,30,30

 30 S(15)=STEST

 GO TO 50

 40 S(15)=-STEST

 50 IF (S(13).GE.1) THEN

 S(20)=PI/2

 ELSEIF (S(13).LE.-1) THEN

 S(20)=-PI/2

 ELSE

 S(20)=ASIN(S(13))

 ENDIF

C fixed error on 6/7/94

 RATIO=S(14)/COS(S(20))

 IF (S(15)) 70,60,60

 60 IF (RATIO.GE.1) THEN

 S(21)=PI/2

 ELSEIF (RATIO.LE.-1) THEN

 S(21)=-PI/2

 ELSE

 S(21)=ASIN(RATIO)

 ENDIF

 GO TO 80

 70 IF (RATIO.GE.1) THEN

 S(21)=PI/2

 ELSEIF (RATIO.LE.-1) THEN

 S(21)=3*PI/2

 ELSE

 S(21)=PI-ASIN(RATIO)

 ENDIF

 80 S(20)=180.*S(20)/PI

 S(21)=180.*S(21)/PI

C

C Calculate direct normal, diffuse and ground reflected radiation

C (no clouds)

C Revision 7/19/94 

 IF (S(13).LT.0.000001) THEN

 S(24)=0.

 ELSE
47

 S(24)=S(30)*S(8)*S(33)*EXP(-S(31)/S(13))

 ENDIF

C end of revision

 S(22)=S(32)*S(24)/S(8)/S(8)

 S(23)=S(7)*(S(22)+S(24)*S(13))

C CALCULATE SOLAR RADIATION ON SURFACE (CLOUDLESS)

 WT=S(10)*PI/180.

 WA=S(9)*PI/180.

 S(16)=COS(WT)

 S(17)=SIN(WA)*SIN(WT)

 S(18)=COS(WA)*SIN(WT)

 S(19)=S(16)*S(13)+S(17)*S(14)+S(18)*S(15)

 S(34)=S(24)*S(19)

C If roof tilt <45, diffuse radiation same as on horizontal surface

C Tilt >45: use equations for vertical surface

 IF (S(10).LT.45.) THEN

 Y=1.

 ELSE

 Y=0.55+0.437*S(19)+0.313*S(19)**2

 IF (S(19).LE.-0.2) Y=0.45 

 ENDIF

 IF (S(19)) 110,110,120

 110 S(19)=0.

 S(34)=0.

 120 CONTINUE

 S(26)=S(22)*Y

 S(27)=S(23)*(1-S(16))/2. 

 S(25)=S(34)+S(26)+S(27)

 IF (S(25).LT.0.000001) S(25)=0.00

C

C Calculate radiation on horizontal surface

 150 SHORIZ=S(24)*S(13)+S(22)

C

C Convert measured solar data from Langleys to Btu/h.ft2

 SDATA=RADDAT(LCOUNT,KK)

 IF (SDATA.GT.9000) SDATA=0.

 IF (SDATA.LT.0.) SDATA=0.

 SDATA=0.3687*SDATA

C

 IF (SHORIZ.LT.0.000001) THEN

C If SDATA not zero, radiation is assumed to be diffuse:

 SRAD(LCOUNT,KK)=SDATA

 ELSE

C Calculate adjustment factor from measured data (factor includes

C effect of cloudcover)

 FACTOR=SDATA/SHORIZ

 SRAD(LCOUNT,KK)=FACTOR*S(25)

 ENDIF

C

 1200 CONTINUE

 1300 CONTINUE

 OPEN(7,FILE=FILENM,STATUS='NEW')

 WRITE (7,1500) ((SRAD(L,I), I=1,24), L=1,LCOUNT)

 CLOSE (7,STATUS='KEEP')

 1500 FORMAT (' ',184(24(F7.2,' '),/))
48

C

 WRITE (*,300)

 300 FORMAT(' MORE RUNS WITH SAME MEASURED SOLAR DATA? (0=NO 1=YES)'\)

 READ (*,*) NRUN

 IF (NRUN.EQ.1) GO TO 1000

 END