SEISRISK Ill: 

A Computer Program for 

Seismic Hazard Estimation 

U.S. GEOLOGICAL SURVEY BULLETIN 1772 

SEISRISK Ill: 

A Computer Program for 

Seismic Hazard Estimation 

By BERNICE BENDER and DAVID M. PERKINS 

U.S. GEOLOGICAL SURVEY BULLETIN 1772 
DEPARTMENT OF THE INTERIOR 

DONALD PAUL HODEL, Secretary 

U.S. GEOLOGICAL SURVEY 

Dallas l. Peck, Director 

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON: 1987 

For sale by the 

Books and Open-File Reports Section 

U.S. Geological Survey 

Federal Center 

Box 25425 

Denver, CO 80225 

The computer program listed in this report has been used successfully for the purposes and on the system 

described. The U.S. Geological Survey is not responsible for its further success under other circumstances. 
CONTENTS 

Summary 1 

Introduction 1 

Modeling earthquake location uncertainty 2 

Seismic hazard model for homogeneous sources 5 

The fault-rupture model 6 

Smoothing of accelerations calculated using fault ruptures 6 

Magnitude smoothing of closest distance ruptures 7 

Distance smoothing for ruptures along artificial parallel faults 8 

Modeling acceleration variability 9 

Calculating ground-motion exceedance probabilities 10 

Model implementation--details and requirements 10 

Geometry 10 

Acceleration table levels 11 

Source zone computations 11 

Area computations 11 

Grid spacing requirements for modeling uncertainty 12 

Using sources with different cr values 13 

Fault computations 14 

Rupture length computations 14 

Magnitude smoothing 15 

Attenuation function 15 

Halving the magnitude intervals 15 

Output options 15 

Accumulating rates in successive runs 16 

References cited 16 

Appendix A 16 

Appendix B 20 

Appendix C 24 

FIGURES 

1. Effect on source-zone seismicity of four different degrees of smoothing 3 

2. Probabilistic accelerations for three different degrees of smoothing 4 

3. Rupture spanning three fault segments 6 

4. Acceleration exceedance curves calculated using three different magnitude 

increments 7 

5. Effect of distance smoothing for traverse across three faults 8 

6. Differences in acceleration levels calculated with and without attenuation 

variability 9 

7. Four points defining new equator and seismic-felt area 11 

8. Inputs for identifying quadrilaterals in a source-zone 12 

9. Arcs used to approximate fraction of source zone at a given distance 

from site 12 

10. Extension of local grid for calculating effect of future earthquake location 

uncertainty 13 

11. Flowchart of source and fault computations 14 

A 1 . Directions of increasing distance along fault from site 17 

A2. Relocating fault segment to x-axis 18 

A3. Effective local coordinates when rupture overlaps three fault segments 19 

Contents Ill 

SEISRISK Ill: 

A Computer Program for Seismic Hazard Estimation 

By Bernice Bender and David M. Perkins 

SUMMARY 

SEISRISK Ill is a revision of SEISRISK II. a computer program for seismic hazard analysis. described in 

Open File Report 82-293. Specifically. SEISRISK II and 

SEISRISK Ill were designed to compute maximum groundmotion levels that have a specified probability of not being 

exceeded during fixed time periods at each of a set of sites 

uniformly spaced on a two-dimensional grid. The principal 

differences between SEISRISK II and SEISRISK Ill are: 

1. SEISRISK II assumes that seismicity within a seismic source zone is uniform; that is. each point within a 

source zone has the same probability of being the epicenter of a future earthquake. This assumption means that 

the projected rate of earthquakes changes abruptly at a 

source zone boundary; an effect of such a change in seismicity is that calculated probabilistic ground-motion levels may differ substantially at sites a few kilometers apart 

near a boundary. SEISRISK Ill retains the concept of seismic source zones but allows earthquakes within a zone to 

be normally rather than uniformly distributed. This allows some of the earthquakes that previously would have 

occurred within the zone to occur outside the zone. permitting seismicity to vary smoothly across the boundaries 

of the zone. The result is that calculated acceleration levels also vary more smoothly at sites near a boundary. 

2. SEISRISK Ill does a partial "magnitude smoothing" 

that treats the "'closest distance" ruptures as if they occurred over a range of magnitudes. In the fault-rupture 

model employed. ground motion (for example. acceleration) at a site resulting from a rupture along a fault is 

regarded as a function of closest site-to-rupture distance. 

This means that for the longer ruptures associated with 

higher magnitude earthquakes. a substantial fraction of 

the possible ruptures of a given magnitude along the fault 

may be at the same (closest) distance from the site. and 

the same ground-motion level will be calculated for all of 

these ruptures. Because the computations are performed 

for a discrete set of magnitudes. the ground-motion levels produced by the ruptures at closest distance for each 

of these magnitudes will be calculated as occurring more 

frequently than other ground-motion levels. This has the 

effect that the calculated probability of exceeding a specified ground-motion level does not vary smoothly as a 

function of ground-motion level. Increasing the number 

of magnitudes at which the calculations are performed reduces these irregularities but increases the computation 

time. By treating the closest distance ruptures for each 

magnitude as if they occurred over a range of magnitudes. 

SEISRISK Ill smooths the acceleration densities resulting 

from these ruptures. 

3. The fault pattern may be quite complex within an active fault zone. and faults may be spread over a wide area. 

In both SEISRISK II and SEISRISK Ill. ··artificial" parallel faults may be modeled when the presence of faulting 

is inferred from geologic considerations. and the orientation of the faults is assumed to be known. If a series 

of equidistant parallel faults are input into the program. 

SEISRISK Ill can do a partial "'distance smoothing" to 

simulate a finer spacing between faults in order to better 

approximate a uniform distribution. (SEISRISK II did not 

permit this smoothing.) 

4. Both SEISRISK II and SEISRISK Ill calculate probabilistic ground motions at all sites resulting from earthquakes in a single source. before proceeding to the next 

source. This requires retaining intermediate calculations 

for each site to accumulate ground-motions from successive sources. SEISRISK II wrote these intermediate results onto a disk; SEISRISK Ill saves the results in twodimensional arrays in memory. Currently. the program requires two arrays dimensioned (56.N). where N =number 

of sites on grid of sites for which calculations are performed. 

INTRODUCTION 

In most seismic hazard analyses, earthquakes are 

modeled as 

1. points located randomly within seismically homogeneous areal source zones; and (or) 

2. finite length ruptures that occur randomly along linear fault segments. 

Probabilistic models are used because sizes and locations of future earthquakes cannot be predicted precisely. Source zones are typically regarded as seismically 

homogeneous, that is, each point within a source zone is 

assumed to have the same probability of being the epicenter of a future earthquake, because the analyst has 

Introduction 
no reason to assume that an earthquake is more likely 

to occur at one point than at another within the zone. 

The analyst uses catalogs of recorded earthquakes, 

and whatever geologic, geophysical or other information is available to help him to determine seismic som·ce 

zones. However, as Thenhaus (1983, p.31) noted, "The 

procedures used in delineating seismic source zones are 

ill defined. No single standard exists by which source 

zones across the nation can be drawn, primarily because 

of the nonuniform level of pertinent seismological, geological and geophysical information available for areas 

of vastly differing tectonic and geologic settings. The 

equivocal association of seismicity with geologic structure throughout most of the United States compounds 

the problem." Because the evaluation and interpretation of the various types of information available depends strongly on individual judgment or opinion, different experts may define source zones differently. 

Assuming that seismicity is uniform within a source 

zone means that seismicity changes abruptly at each 

source-zone boundary. An effect of these abrupt changes 

in seismicity is that probabilistic ground-motion levels calculated at, sites a few kilometers apart near a 

boundary may differ significantly. Bender ( 1986) noted 

changes of 50-80 percent or more in acceleration levels calculated for longer exposure times at sites 20 km 

apart near a boundary of a source zone. Such changes 

seem unreasonable if, in fact, source zones are not well 

determined. 

To avoid sudden changes in seismicity at sourcezone boundaries, SEISRISK III, unlike SEISRISK II, 

permits an option of "earthquake location uncertainty". 

Source-zone boundaries are specified in SEISRISK III as 

in SEISRISK II, but now each point within the zone is 

regarded as the mean or most likely location of a future 

earthquake; the locations of actual earthquakes are normally distributed with standard deviation a about their 

mean locations (a is an input parameter). This permits 

earthquake rates to vary smoothly near a boundary, 

without being significantly different in the center of the 

region from the rates assumed in the uniform-seismicity 

model (fig.1). As a increases, a higher percentage of 

the earthquakes that were expected to occur within the 

source zone when no location uncertainty (a=O) is assumed, occur instead beyond the boundaries. 

In the following sections, we discuss in more detail earthquake location uncertainty and how this uncertainty is implemented in SElSRISK III. We also discuss 

"magnitude smoothing" and "distance smoothing" for 

ground-motions resulting from ruptures along a fault, 

and the parameters that control this smoothing in the 

computer program. A listing of SEISRISK III (Appendix C) and information regarding data input and 

formats (Appendix B) are also included. 

2 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

TERMINOLOGY 

In this report, we use acceleration as ?., generic measure of ground shaking; we could equally well have used 

velocity or any othfr measure. By acceleration occurrences in an interval a3

·_ 1 <a~ ai we mean the average annual occurrence rate of accelerations in the range 

ai-1 < a ~ ai experienced at a specific site. Acceleration levels or probabilistic acceleration levels refer to 

the accelerations that are calculated as having a specified probability of not being exceeded during a given 

time period at particular sites. 

MODELING EARTHQUAKE 

LOCATION UNCERTAINTY 

In the model implemented in SEISRISK III, the 

seismicity associated with each point in a source zone 

is regarded as normally distributed about that point. 

The possible earthquakes, rather than being uniformly 

distributed throughout a source region, have mean locations that are uniformly distributed, with scatter about 

each mean location given by a circular-normal distribution. Stated more precisely, in this model, if a is 

the standard deviation in location, and if the expected 

(mean) location of an earthquake is at (Xe, Ye), the 

probability that the earthquake will occur instead in 

a small area A surrounding (Xe+Llx, Ye+Lly) is given 

by 

A Llx2 + .Lly2 

PA(.Llx,.Lly)=-2 exp( 2 ). (1) 21ra 2a 

Figure 2 (A-C) shows acceleration levels calculated 

as having a 90-percent probability of not being exceeded 

during exposure times of 10, 50, 250, and 1000 years 

when no variability in earthquake location (a = 0) is 

assumed, and when a= 10 and 20 km (the attenuation 

function of Joyner and Boore (1981) with aa = 0.5 in 

loge acceleration was used). The sites are located 5 km 

apart on a line through the centers of two contiguous 

source zones. Each zone is a square 60 km wide; the 

second zone has an earthquake rate 10 times higher than 

the first. Unless a is a substantial fraction of the length 

or width of the source zone (for example, a = 20 km 

in fig.2), the acceleration levels calculated at sites near 

the center and at sites on the boundaries of the source 

zone are nearly the same for various values of a. For 

a site outside the source, if a =F 0 is assumed, some 

of the earthquakes will occur closer to the site than 

when a=O is assumed, and these closer earthquakes will 

produce higher accelerations at the site than will any 

earthquakes of the same magnitude occurring within 

the fixed sourde zone. 

Figure 2D shows the difference in acceleration level 

calculated for the 1000-year exposure time at sites 5 
A 

Figure 1. Earthquake rates within gridded areas ( 10 

km x 10 km). In A earthquakes occur only within two 

well-defined seismically homogeneous source zones. In B. 

earthquake locations are normally distributed with standard deviation a= 10 km; in C, a= 25 km: in D. a= 50 

km. 

km apart on the line when a= 0, 10, and 20 km; for 

example, a change in acceleration of +0.05 g at the site 

at X km means that the acceleration level calculated 

at the site at X+5 km is 0.05 g higher than the level 

calculated at X. 

To explain how acceleration levels are calculated in 

SEISRISK III when earthquake-location variability is 

taken into account, we first observe that, if earthquakes 

are normally distributed, those that occur at. distance 

(Ax, Ay) from their mean locations, in effect, create a 

new source zone translated by (Ax, Ay) from the original zone. This means that the possible locations of the 

source zone (rather than of the earthquakes) may be 

B 

regarded as having a normal distribution. ~ore explicitly, the source. zone is said to be located at (Ax, Ay) 

if an arbitrary point in the zone that would have been 

at (Xb, Yb) if the zone were at its mean or most likely 

location, lies instead at (X,,+ Ax, Yl,+Ay). The probability that the zone will be located at (Ax, Ay), is 

Ax2 +Ay2 

proportional to exp(- 2 ), where a is now the 

2a 

standard deviation in source location. Furthermore, the 

acceleration occurrences in an interval ai-I < a~ ai that 

are calculated for a site at (X, Y) when the source location is moved to (Ax, Ay) are the same as those calculated for a site a.t (X- Ax, Y-Ay) when the source 

is a.t its original location. This means that the effect 

of source-zone location uncertainty a.t a site at (X, Y) 

can be produced by computing a weighted average of 

the acceleration occurrences at various sites, calculated 

using the original seismically homogeneous source zone. 

This last. fact was used to implement earthquake location (or, equivalently, source-zone location) uncertainty 

in SEISRISK III as follows. 

Modeling Earthquake Location Uncertainty 3 
0.50 

A 

u=O 

> 0.40 

~ 

> 4: 

a: 0.30 ~ 

z 

0 

i= 

4: 0.20 a: 

w 

....J 

w 

u 

u 0.10 4: 

0.00 

0 60 120 180 240 

SITES (km) 

0.50 

~ 0.40 

> 

-~ 

~ 0.30 

z 

0 

i= 

4: 

a: 0.20 w 

....J 

w 

u 

u 

4: 0.10 

SITES (km) 

Figure 2. (A)-( C): Accelerations calculated to have a 

90 percent probability of not being exceeded during various time periods. at sites on a line through the centers 

of two seismic sources. Source II has an earthquake rate 

10 times higher than Source I. Earthquake locations are 

normally distributed with A, a=O km; B, a=10 km; C, 

a=20 km. D, Differences in acceleration levels calculated 

above for 1000-yr period at sites 5 km apart on the line 

of sites. 

SEISRISK III begins, as does SEISRISK II, by calculating expected annual accelerations at each site on a 

grid of sites, resulting from earthquakes in well-defined, 

seismically homogeneous zones. Sites are spaced at increments dx, dy = dx. The acceleration occurrences 

at a site at (XI, Y1) that result when earthquake location variability is taken into account is a. weighted 

sum of the acceleration occurrences that were calculated at each point on the grid, assuming homogeneous 

4 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

0.50 

B 

u=lO 

0.40 

0.30 

0.20 

0.10 

0.00 

0 60 120 180 240 

SITES (km) 

0.10 

1000 

0.08 Year D 

E 

.::.!. 

It) 0.06 -

z 

0 0.04 i= 

4: 

a: 0.02 w 

....J 

w 0.00 u 

u 

4: -0.02 z 

w -0.04 (.!) 

z -0.06 4: 

I 

u -0.08 

-0.10 

0 60 120 180 240 

SITES (km) 

source zones. The weight assigned to acceleration occurrences in the interval ai-l <a ~ ai at the point at 

(X1 +k1dx, Y1 +k2dy) is proportional to 

(2) 

The occcurrences in the interval a,-_1 < a ~ a,- that 

result when earthquake location variability is taken into 

account is given by 

Lk Lk., w(kb k2)PI(XI +k1dx, Y1 +k2dy), j] D = 1 -

l:kl Lk2 w(kl, k2) 

(3) 

where pi( X 1 + k 1 dx, Y1 + k2dy), j]=the acceleration occurrences calculated in the interval ai-l < a ~ ai at 

the site at xl +kl dx, yl +k2dy assuming no earthquake 

location variability (i.e., homogeneous source zones). 
SEISMIC HAZARD MODEL 

FOR HOMOGENEOUS SOURCES 

To determine acceleration rates in SEISRISK III 

that result when earthquake location variability is taken 

into account, we first compute acceleration rates, as 

in SEISRISK II, assuming that earthquakes occur randomly within seismically homogeneous source zones or 

along linear fault segments. In this section, we briefly 

describe the model before incorporating earthquake location uncertainty. 

Earthquake occurrences are assumed to have a Poisson distribution, and rates that remain constant during 

the time periods of interest. Mean or median groundmotion from an earthquake is an increasing function 

of magnitude and decreasing function of site-to-source 

distance. Earthquakes occur as points within quadrilateral source zones or as finite-length ruptures along 

linear fault segments. Ruptures of a given length may, 

with equal likelihood, be centered at any point on the 

fault for which neither end of the rupture would extend 

beyond the end of the fault. 

Earthquakes in a given zone or fault are restricted 

to occur within a specified magnitude range; the range 

may be different for different zones and faults. More 

specifically, for each zone (and fault) a maximum magnitude mmax, a minimum magnitude mo, and a magnitude interval .6.m are assumed. All magnitudes occurring within a magnitude interval are grouped at the 

center of the interval, that is, at a set of n distinct magnitudes mi: 

mi=mo+(i+l/2).6.m, 0 ~ j ~ n-1, 

Llm = mmax -mo. (

4) 

n 

The rate of earthquakes per unit time in the zone (or on 

the fault) in each magnitude interval is specified by the 

program user; it frequently follows (but is not required 

to follow) a relationship of the form 

logNm=a-bm (5) 

where a and b are constants. If Nrnli) is the total number of earthquakes in the ph magnitude interval, the 

fractional occurrences expected in any small subarea 

AA is 

AA 

ANm(il 

where A =total area of zone. 

(6) 

Seismicity is assumed to remain constant during 

the time periods being considered; that is, the av~rage 

rate of earthquakes per unit time for each magmtude 

interval does not change with time. 

A table of median (or mean) ground-motion values (for instance, acceleration) for the desired attenuation function is entered into the program for a set 

of earthquake magnitudes and distances. The program 

assumes that ground-motion increases with magnitude 

and decreases with site-to-source distance; ground motions for various magnitudes and distances are obtained 

as needed by interpolating values in the table. 

Given the source zones, earthquake rates and attenuation function, an average yearly rate of accelerations within each of 55 acceleration intervals is constructed for each site. (In SEISRISK II 100 intervals 

were used; the number was decreased to save time and 

space, but can be adjusted within the program.) Acceleration occurrences in the range ai-l < a ~ ai are 

h ·th . accumulated into reg(i), the accumulator fort e t mterval, as described below. 

For the given attenuation function, let dm (i) be 

the distance at which a magnitude m earthquake produces median acceleration ai at a site. The probability that an acceleration in the range ai-l < a ~ ai occurs at the site is the probability that a magnitude m 

earthquake occurs at a distance drn(i) ::; d < dm(i-1) 

integrated (summed) over all permissible magnitudes. 

(Note dm(i)::; dm(i-1) because shorter di_stances produce higher accelerations.) The probability that an 

earthquake of magnitude m occurring within the zone 

(or on the fault) occurs within the distance interval 

drn ( i) ~ d < drn ( i-1) is equivalent to the fraction of the 

area of the zone (or length of the fault) that lies within 

this interval. The area (or length) within the interval 

A[m(i)] = f[dm(i), Adm(i)], (7) 

is a function of interval width 

D.d(i)=dm(i-1)-dm(i) (B) 

and interval midpoint 

- (") dm(i-1}+dm(i) d t = . m 2 (9) 

The rate of accelerations in the range ai-l < a ::; ai at 

the site resulting from magnitude m earthquakes withi~. 

the zone then 1s equal to the rate of magnitude m earthquakes within A[m(i)]: 

. ( .) _ A[m(i)] N 

I'm ,z - A m (10) 

where Nm=rate of magnitude m earthquakes for the 

entire zone or fault 

A '=total area of zone (or total fault length). 

Seismic Hazard Model for Homogeneous Sources 5 
The program determines A[m(i)] for each magnitude and adds the corresponding rate Pm (i) into the ith 

accumulator, reg(i). The total rate accumulated for accelerations in the range ai-l <a ~ ai from one source 

zone then' is the sum over magnitude of 

n 

p(i) = L Pm1 (i) ( 11) 

j=l 

where n=total number of magnitude intervals. 

The process is repeated for all 1 ~ i ~55, for all sources. 

These calculations assume a single value of ground 

motion (acceleration) for earthquakes of a given magnitude and distance. Calculations that take into account 

ground motion variability will be discussed later. First, 

we shall consider the calculation of (median) accelerations resulting from ruptures along faults and some 

ways of smoothing the calculated acceleration densities. 

THE FAULT-RUPTURE MODEL 

As previously noted, both SEISRISK II and SEISRISK III allow earthquakes to be modeled as finitelength ruptures along linear fault segments. In both 

programs, acceleration occurrences resulting from fault. 

ruptures are added to those resulting from earthquakes 

in source areas. However, unlike SEISRISK II, SElSRISK III permits a partial "smoothing" of acceleration occurrences at a site caused by ruptures along a 

fault. We shall discuss the fault-rupture model and the 

changes that have been made to smootl1 those occurrences that result from fault ruptures. 

In the fault-rupture model, sets of one or more 

faults are defined as fault zones; if a fault zo11e contains 

more than one fault, earthquakes in the zone are distributed among the faults in proportion to the lengths of 

the faults. An individual fault of length L may consist of 

as many as 24 connected straight-line segments. (The 

maximum number of segments may. be altered within 

the program.) A rupture may span several segments 

(fig.3). Ruptures occur with equal likelihood anywhere 

along the fault, so long as they are wholly contained 

within the fault. Rupture lengths are a function of 

magnitude and are lognormally distributed with median length l( m) and standard deviation a,. 

As with source areas, the ground-motion at a site is 

a function of earthquake magnitude and distance. However, in the case of faults, distance from the site to a 

rupture is always taken to be the distance to the point 

on the rupture closest to the site. (The fault-rupture 

model is described in detail in Bender, 1984b.) All magnitude m ruptures along the fault whose closest distance 

6 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

......... ......._ _______ _ RUPTURE 

Figure 3. Fault consisting of three connected straightline segments. A rupture may span several segments. 

Closest distance from rupture to site is used to determine 

ground motion at the site. 

d to the site is in the range dm (J') ~ d < dm (J'- 1) will 

produce acceleration a, ai-l <a ~ ai, at the site. The 

calculations to determine site-to-rupture distance and 

the fraction of ruptures of each magnitude within each 

distance range are described in Appendix A. 

Smoothing of Accelerations Calculated 

Using Fault Ruptures 

In the fault.-rupture model, a substantial fraction of 

the ruptures of earthquakes of a given magnitude may 

occur at the same closest distance from the site and may 

produce the' same median acceleration at the site. For 

example, assume the fault extends on the x-axis from 

0 ~ x ~ L. A rupture of length L/4 may be centered at 

any point x along the fault such that 1/BL < x ~ 7 /BL. 

If the site is located on a line perpendicular to the center 

of the fault at. (L/2, P), all ruptures that have their 

centers in the range 1/4L~ x~ 3/4L along the fault will 

be at the same (closest) distance P from the site. This 

means that two-thirds of the possible ruptures of length 

L/4 will be at distance P from the site; for rupture 

length L/2, all possible ruptures will be at distance P 

from the site. Because longer ruptures are associated 

with higher magnitude earthquakes, the fraction of the 

ruptures that will be at closest distance P from the site 

increases with magnitude. 

Because calculations are normally done only at a 

set of magnitudes mi spaced ~m apart, ruptures of 

magnitude mi earthquakes at distance P from the site. 

will cause a "spike" of acceleration occurrences to be 

accumulated in some reg(j). (Recall that occurrences 

of accelerations in an interval ai-l < a ~ ai are accumulated in a corresponding accumulator reg(j).) The 

number of spikes in the acceleration occurrences at a 

site depends on the number of magnitudes used; the 

locations and sizes of the spikes depend upon the actual magnitude values. If no acceleration variability 

is assumed (for example, mean or median acceleration 

values only are used), or if aa, the standard deviation 

in log acceleration, is small, spikes in acceleration occurrences will result in sharp changes in the accelera-
tion exceedance curves at the accelerations at which the 

spikes are present. Furthermore, these spikes may occur 

at different accelerations at different sites, causing the 

probabilistic accelerations calculated at different sites to 

vary irregularly. SEISRISK III does a partial smoothing 

of acceleration occurrences resulting from the ruptures 

at distance P, as described below. 

Magnitude Smoothing for Closest 

Distance Ruptures 

Earthquakes at a fixed distance P from the site, 

with magnitudes in the range m(j -1) ~ m < m(j) will 

produce median accelerations in the corresponding interval aj-1 ~ a < a3. If earthquakes within the magnitude interval mi -l:::.m/2 ~ mi ~ mi+l:::.m/2 are assumed 

to follow a Gutenberg-Richter magnitude-frequency relationship, the number h(j) of those earthquakes with 

magnitudes that lie between m(i- 1) ~ m ~ m(j) can 

be determined (h(j) contains magnitudes in the range 

max!mi-l:::.m/2, m(i-1)] ~ m~ min!m(j.), nli+l:::.m/2]). 

If f(i) is the fraction of ruptures of length li{m) of magnitude mi earthquakes that occur at distance P from 

the site, approximately the same fraction of earthquakes 

will be at distance P for all magnitudes in the range 

mi -l:::.m/2 ~ nli ~ nli +l:::.m/2, for small values of l:::.m. 

This means that a rate of accelerations r(i, j) = f (i)h(j) 

resulting from ruptures at distance P for magnitudes 

mi-l:::.m/2 ~ mi ~ mi+l:::.m/2 can be assigned to the 

acceleration interval a3-1 ~a< a3. 

The solid lines in figure 4 show the yearly rate of 

exceedances of various accelerations calculated at a site 

at a perpendicular distance of 1 km from the center of 

a 300 km long fault, when the magnitude range is 5.0 ~ 

m ~ 8.0 and the magnitude increments are l:::.m = 0.6, 

0.3, and 0.15, respectively. The dashed lines in figure 

4 were obtained in each case by magnitude smoothing 

of the spikes of acceleration occurrences. Very little 

difference appears in the smoothed results for .tlm= 0.6, 

0.3 and 0.15, suggesting that for sites near the center 

of the fault a magnitude increment of Am= 0.6 might 

suffice if the spikes are smoothed in this manner. 

Figure 4. Yearly acceleration exceedance rates calculated at a site at a perpendicular distance of 1 km from 

the center of a 300-km-long fault. using magnitude increments A. 6.m=0.6: B. 6.m =0.3: C. 6.m=0.15. The 

magnitude range is 5.0 S m ::::_: 8.0. Attenuation curves 

are those of Schnabel and Seed. (1973). The dashed line 

in each case show the exceedances calculated when the 

acceleration density spikes resulting from ruptures that 

intersect the point on the fault closest to the site are 

smoothed for each magnitude. 

10-4 

A 

0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9 1.0 

ACCELERATION-PERCENT GRAVITY 

0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9 1.0 

ACCELERATION-PERCENT GRAVITY 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

ACCELERATION-PERCENT GRAVITY 

The Fault-Rupture Model 7 
In smoothing acceleration occurrences, occurrences 

normally assigned to magnitude mma.x -l:l.m/2 are distributed throughout the interval mmax -l:l.m/2 ~ m ~ 

mma.x + l:l.m/2; in particular, some of the occurrences 

are assigned to the interval containing the acceleration 

produced at distance P by the maximum magnitude, 

mma:r. However, when magnitude smoothing is not 

done, the highest. magnitude used in the calculations 

continues to be m 111 a:1 -l:l.m/2, (equation 4). For a site 

beyond the end of the fault at (X, P) {either X < 0 

or X > L), there are no ruptures at distance P, and, 

hence, no acceleration occurrence spikes from ruptures 

at c1osest distance; therefore, no magnitude smoothing 

is done for acceleration occurrences at such a site. This 

implies, in effect, a lower maximum magnitude for sites 

beyond the end of the fault. 

Distance Smoothing for Ruptures Along 

Artificial Parallel Faults 

The fault pattern may be quite complex within an 

active fault zone, and faults may be spread over a wide 

area. SEISRISK III allows one to model artificial parallel faults when the presence of faulting is inferred from 

geologic considerations and the orientation of the faults 

is assumed to be known. In this model, faults with a 

specified orientation are assumed to be equally likely 

anywhere within a given area; for example, faults parallel to the x-axis are equally likely to have their ycoordinate anywhere in the range 0 ~ y ~ w. (The 

y-coordinate of the fault is uniformly distributed in 

this model.) To approximate a uniform distribution of 

faults, coordinates of evenly spaced parallel faults must 

be entered into the computer program. However, unless 

the spacing between faults is quite small (for example, 1 

or 2 km), or unless some sort of smoothing is done, the 

acceleration levels calculated at sites a short distance 

apart can be quite different. Figure 5 illustrates how 

calculated acceleration levels may change as a function 

Figure 5. Accelerations calculated to have a 90 percent 

probability of not being exceeded during various time periods at sites on a line perpendicular to the centers of three 

parallel faults. A, 20 km apart. B. the scale is expanded 

and the dashed lines show values obtained when faults are 

placed instead at increments of 1 km over a distance of 

60 km and the original seismicity is divided evenly among 

the faults. C. the solid lines show the results when 20-km 

fault spacing is used in the calculations. but ruptures on 

each fault that occur "as close as possible" to the site 

are placed instead on a line perpendicular to the original 

faults; the dashed lines are the same as in B. The Schnabel and Seed (1973) attenuation curves are used. 

8 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

s:-

I-

> <{ 

a: 

~ 

z 

0 

i= 

<{ 

a: 

w 

...J 

w 

u 

u 

<{ 

~ 

> <{ 

a: 

~ 

z 

0 

i= 

<( 

a: 

w 

...J 

w 

u 

u 

<( 

s:-

I-

> <{ 

a: 

~ 

z 

0 

i= 

<{ 

a: 

w 

...J 

w 

u 

u 

<{ 

0.90 

0.80 

0.70 

0.60 

0.50 

0.40 

0.30 

0.20 

0.10 

0.00 

0 

0.90 

0.80 

0.70 

0.60 

0.50 

0.40 

0.30 , 

, 

0.20 "' 

0.10 

0.00 

250 

0.90 

0.80 

0.70 

0.60 

0.50 

0.40 

0.30 

0.20 

0.10 

0.00 

250 

------1l~ 1000 Year 

Faults 

250 Year 

50 

Year 

100 200 300 400 500 600 

SITES (km) 

1000 Year 

,-

, - -\ 

, \ 

I I 

I I 

I \ 

I I 

I I 

I 

I ' \ 

I 

I ' \ 

I \ 

I \ 

\ I , \ 

\ I \ I ' , 

, ' ' ' 

10 Year 

260 270 280 290 300 310 320 330 340 350 

SITES (km) 

1000 Year 

260 270 280 290 300 310 320 330 340 350 

SITES (km) 
of distance from a fault. Figure SA shows levels calculated at sites on a line perpendicular to the centers 

of three parallel faults, 20 km apart. In figure SB the 

scale was expanded; the solid lines again show acceleration levels for the three fixed fault. locations; the dashed 

lines give results obtained when faults were placed 1 km 

apart and the seismicity adjusted accordingly. 

A distance smoothing is desirable when a single 

fault is used to represent. a uniform distribution of faults. 

For example, assume a fault that extends from (0, 0) to 

(L, 0) was entered into the program t.o represent faults 

with y coordinate anywhere in the range -y0 S y S y0 . 

Accelerations are calculated only for ruptures aloug the 

fault at y=O. If the site is located "above" the fault at 

(X, P) (0 ::; X ::; L), an acceleration density spike will 

result from ruptures at distance P. To smooth the accelerations resulting from these ruptures, assume that 

any rupture with end points (x 1 ,0), (x2 ,0) along the 

original fault, at dista11ce P from the site, may instead 

have end points ( x 1 , y), ( x2 , y) on any fault such that 

-yo :::; y ::; y,,. The seismicity associated with such a 

rupture along a fault at y concentrat.ed at the point on 

the fault closest to the site at (X, y). This means that if 

faults [with end points at (0, y), (L, y), -y0 ::; yS. y0 ], are 

uniformly distributed parallel to the original fault, ruptures at closest distance on each fault will produce the 

same accelerations as points uniformly distributed on a 

line through the site perpendicular to, and extending 

±yo km from, the original fault. Acceleration occurrences in each interval ai-l ::; a:::; ai may be determined 

from the frequency of the "point-source" earthquakes of 

each magnitude at the appropriate distances along this 

perpendicular line. In figure 6C, the ruptures at closest 

distance to each site were treated as point-source earthquakes along a line through the site and perpendicular 

to the faults; the calculation was done using Am= 0.6. 

MODELING ACCELERATION VARIABILITY 

The ground-motion computations that have been 

described assume a single value of acceleration results 

from earthquakes of each magnitude and distance. However, in fact, a range of accelerations results; this variability in acceleration is commonly modeled by assuming that accelerations from earthquakes of a given magnitude and distance are lognormally distributed with 

standard deviation aa in loge acceleration. Assuming 

a lognormal distribution of accelerations tends to cause 

the calculated probabilistic ground-motion levels at all 

sites to be higher than the levels calculated when only 

median accelerations are used (for example, Bender, 

1984a). 

E 

~ 

10 

........ 

z 

0 

i= 

<( 

a: 

UJ 

..J 

UJ 

u 

u 

<( 

~ 

UJ 

(!) 

z 

<( 

J: 

u 

0.10 I 

1000 Year : 

0.08 I 

I 

I 

I 

I 

0.06 I 

I 

I 

I 

0.04 I 

I 

I 

0.02 

0.00 

-0.02 

-0.04 

-o.06 

-o.o8 

-Q.1 0 +---.---+----.,....---+---.---1-----.r-----. 

0 60 120 1.80 240 

SITES (km) 

Figure 6. Differences in acceleration levels calculated 

for 1000-yr period at sites 5 km. apart. assuming a0 =0 

and assuming a 0 =0.5 in loge-acceleration in the attenuation function of Joyner and Boore (1981): no uncertainty 

(cr=O) in earthquake location. 

Figure 6 shows the differences in probabilistic acceleration levels calculated for a 1000-year exposure time 

at sites 5 km apart, when the sites and sources are the 

same as in figure 2. The calculation was done assuming no uncertainty in earthquake location using median values in the Joyner and Boore attenuation and 

using Oa = 0.5 in logr, acceleration. (The calculation 

for Oa = 0.5 is the same as that shown in figure 2D for 

the 1000-year exposure period for a= 0.) Figure 3 was 

included to illustrate that taking variability in the attenuation function into account, rather than reducing 

the differences in levels at sites near a boundary1 may 

increase the site-to-site acceleration differences above 

those calculated when only median values are used in 

the calculations. This is discussed in det.ail in Bender 

(1986). 

Assuming that a11 , the standard deviation iu log 

acceleration: i~ independent of magnitude and distauce 

permits SEISRISK III to take acceleration variability 

into account only after doiug calculations using median 

acceleration values for all earthquakes. Median rather 

than mean acceleration values are used initially because 

log acceleration (ratl1er than acceleration} is assumed to 

be 11onnally distributed; the median value of acceleration corresponds to the mean value of log-acceleration. 

Let us assume occurrences of median accelerations 

m each interval ar1 < a ::; a.i have been accumulated 

m the corresponding reg(J'). Let regs(J') be a second 

Modeling Acceleration Variability 9 
array m which t.o accumulate occurrences of accelerations aJ_1 < a :::; aJ after acceleration variability has 

been taken into account. SEISRISK III redistributes a 

fraction of the acceleration occurrences accumulated in 

reg(j) to regs(k) (for each j and k) as follows. 

The entry in reg(j) is assumed to be concentrated 

at acceleration a where 

1 (

-) loge[aj-l]+loge[aiJ 

oge a = 2 . (12) 

In a lognormal distribution, given loge a, the probability 

of an acceleration in the range a 1 ~ a ~ a2 is the area 

under the normal probability curve: 

[ (x-logea) 2 ]d exp - 2 x. 2cra 

(13) 

From equation 13, a. fraction f(k) of the accelerations 

in reg(i) is placed in regs(k), 

(14) 

where p(i)=occurrences of (median) accelerations in 

the range ai _ 1 < a ~ ai accumulated in reg{j). Accounting for acceleration variability in this manner may 

be considerably faster computationally than methods 

which include acceleration variability for earthquakes 

of each magnitude and distance. Accounting for acceleration variability after calculations for median accelerations have been completed is possible because of 

the assumption that era, the standard deviation in log 

acceleration is independent of distance and magnitude. 

CALCULATING GROUND-MOTION 

EXCEEDANCE PROBABILITIES 

Given that average annual acceleration occurrence 

rates in the interval aj-l < ~ a~ ai have been accumulated in regs{j), 1 ~ j ~ 55, Ex[ai], the yearly occurrence rate of accelerations a> ai, is simply 

Ex[aiJ = L p* (i) (15) 

i=j+l 

where p* (i) =accelerations in the range ai-l< a~ ai accumulated in regs(i). The yearly exceedance rate Ex( a) 

of acceleration a can be obtained by interpolating the 

Ex[ai J in equation 15. 

We are now in a position to determine the probability associated with the exceedance of any acceleration. 

Earthquake occurrences are assumed to have a Poisson 

distribution; that is, the probability t.ha.t an earthquake 

will occur is the same at all times regardless of when the 

1 0 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

last earthquake occurred. For a Poisson distribution 

with mean rate JJ-, the probability of exactly k events 

during time t is 

(16) 

The probability of no event during the same time interval is, therefore, 

P(O, t) = exp( -JJ-t). (17) 

Setting J.L =Ex( a), the probability of not exceeding a in 

t years is 

P(O, t) =exp[-Ex(a)t]. (18) 

Taking the natural, logarithm (loge) of both sides of 

equation 18 gives 

Ex( a)= loge[P(O, t)] (19) 

That is, the value of a for which equation 19 is satisfied 

is the acceleration which has probability p = P(O, t) of 

not being exceeded during the time interval t. (Equivalently, a has probability q=1-P(O,t) of being exceeded 

one or more times during time t.) By specifying any 

two of the three parameters (p, a, t), one can determine 

the remaining parameter. However, in SEISRISK III a 

single probability level p and several values (up to 20) 

of time are specified and the corresponding values of a 

are obtained. 

MODEL IMPLEMENTATION-DETAILS 

AND REQUIREMENTS 

In this section, we further describe the model and 

how it. is implemented in SEISRISK III. Input parameters required to run SEISRISK III are identified and described briefly at the beginning of the program listing in 

Appendix A; here we give a more detailed explanation 

of the meaning and roles of some of these parameters. 

Geometry: Transformation to Local Equator 

In both SEISRISK II and SEISRISK III, a pair of 

latitude and longitude values are selected (input) to be 

two points on a great circle that becomes the equator 

in a transformed coordinate system. The great circle is 

chosen to pass through the region of interest, and all 

subsequent input coordinates (for example, end points 

of source areas and fault segments) are transformed relative to the new equator. Furthermore, the subsequent 

calculations are done assuming a rectangular, rather 

than a spherical coordinate system; in the new system 

the (new) equator becomes the x-axis. 
AJ, f/J3 r- - - - - - I ------

(Xz, f/Jz) 

I 

I 

1 --, 

I 1 

I I 

I 

I 

I 1 

I I 

L_ I 

------- I ------.1 A4, r/J4 

Figure 7. Four points define transformation and seismicfelt area. Great circle through () 1 , ¢1 ), () 2 , ¢2 ) becomes 

new equator. {A1, ¢I) -+ (d, 0}, {A2, ¢2) -+ {0, 0). 

()3, ¢3) and (-\4, ¢4) become end points of rectangu!u 

area containing seismic felt points. (-\3, ¢3) -+ (x3, ~· J, 

(.A4, ¢4)-+ (x4, Y4). 

For points in the vicinity of the equator, calculations done using a rectangular coordinate system give 

a good approximation to those obtained using spherical coordinates. Because of convergence of parallels 

of latitude, the error could be substantial if t.he original (long, lat) coordinates are treated as being on a 

flat surface. However, this error is minimized when the 

coordinates are transformed to lie in the vicinity of a 

new equator, so that, for example, within a square region 12 degrees wide and centered on this equator, the 

maximum error introduced by using the straight line 

dist.ance between two points (xi, yi), (xj, Yi) instead of 

the great circJe distance is less than 0.6 percent. 

Sites at which accelerations are calculated lie on a 

rectangular m by n grid in which the site (xi,Yi) is 

Xi = ck (longitude of row i), 1::; i::; m 

Yi = ck (latitude of col j), 1::; j ~ n 

where ck=conversion factor for degrees to kilometers at 

the equator. 

Coordinates (long, lat) of two opposite corners (upper left, lower right) of the seismic felt area in the initial system and grid spaciug are input. to the program 

(fig. 7). Sites in the new coordinate system are located 

at even increments in D.x, D.y (D.y = D.x) at intersections of grid lines. (In SEISRISK II, D.y = D.x was 

not required.) At the end of the computation, sites are 

transformed to the coordinates they would have in the 

original system. Sites may also be points evenly spaced 

on a numbe,~ of lines; in this case, (long, lat) of two end 

points of each line, number of sites per line, and number 

of lines are specified as input. However, calculations for 

sites on a line are done assuming only seismically homogenous source zones; earthquake location variability 

cannot be implemented for sites that are not on a grid 

of sites. 

I I I I I I 

-+-+-+-+-+-+-+-+- 1 I I I I I 

r-- + -+ -+ -+ -+ -+ -+ -+ -

(0, 0) -+---1-+--+-----11---t--+--+-+-~ (d, 0) 

I I 

-+-+-+-+-+-+-+-+-

-+-+-+-+-+-+-+-+- t I a. I I I I I 

-+-+-+-+-+-+-+.-+-

Acceleration Table levels 

A table of 55 acceleration levels ak, such that ak < 

aJM-1 for 1 ~ k::; 55, is constructed within the program. 

A scale factor (input) allows the original range of acceleration values in the table to be expanded or contracted. For a scale factor of one, 50 of the acceleration 

levels are in the range 0.02-1.0 g at increments of 0.02 

g. Four acceleration levels are less than 0.02 g; accelerations at these levels are saved for calculations that 

are done when acceleration variability is taken into account; these accelerations are not printed in the output 

and, hence, are invisible to the analyst. A scale factor 

of 0.5 would permit values from 0.01g-0.5g. (The table 

is created in the subroutine "box" and other values may 

be substituted by replacing this subroutine.) 

Given the table of acceleration levels ak, a corresponding array reg(k) is set aside to accumulate fractional acceleration occurrences, with occurrences in the 

range ai-l < a ::; a3· being placed in reg(j). Accelerations less than a(1) are placed in reg(1) and those exceeding the highest entry a(max) (where max =54) are 

placed into reg( max+ 1). Ear-thquakes as point events 

within quadrilateral source areas and as finite length 

ruptures along linear fault segments may contribute to 

accelerations at a site. We shall discuss these two cases 

separately. 

Source Zone Computations 

Area Computations 

A seismic source zone is defined as a seismically 

homogeneou~ area enc1osed by one or more arbitrary 

quadrilat.erals, connected or disjoint (fig.8). 

Recall that an acceleration "box" represents occurrences of accelerations in some range aj-1 <a::; ai. 

For a specified magnitude, acceleration box boundaries 

[ aj-l, ai] correspond to earthquakes at distances in the 

range r(j)::; r < r(j -1) from the site. The contribution 

to accelerations a3·_1 < a ::; ai produced by magnitude 

m earthquakes in a seismically homogeneous quadrilateral is equal to the rate per unit area within the 

Model Implementation-Details and Requirements 11 
Figure 8. Seismic source zone consisting of two sets of 

quadrilaterals 

Input: For each quadrilateral set 

( 1) jseg ifr tot 

jseg=number+1 of quadrilaterals in this set 

ifr identifies current set, 1 ~ ifr ~ tot 

tot=total number of sets of quadrilaterals 

(2) pairs of quadrilateral corner points. 

In this example. input is 

412 

.A1 f/>1 .A2 f/>2 

.A3 ¢>3 .A4 ¢>4 

.As 4>s .A6 ¢>6 

A7 ¢>7 .As 4>s 

222 

.A9 ¢>9 .Aw ¢>w 

.Au 4>u .A12 f/>12 

quadrilateral, multiplied by the area within the quadrilateral at distances r(j) ~ r < r(j- 1) from the site. 

The area of the quadrilateral within the distance range 

r(j) ~ r < r(j -1) is approximately the area contained 

within a ring of radius r and width dr = !r(J' -1)-r(j) J 

centered at the site (fig.9). 

The distance range over which earthquakes produce accelerations in the interval a3_1 < a ~ a3 depends on magnitude. The accelerations from earthquakes of each magnitude in the appropriate distance 

range within a source area are summed to obtain the 

total contribution from these earthquakes added into 

reg(j). The procedure is repeated for all acceleration 

intervals, and subsequently for all sources. 

To save computation time, for a particular site and 

source, both SEISRISK II and SEISRISK III create a. 

table of distance versus arc-length (for a predetermined 

set of radii) and then later interpolate in this table to 

obtain areas at desired distances. In this case 

Area= ¢>rdr (20) 

where ¢=arc length in radians of that portion of the 

circle of radius r = r3 that lies wholly within 

the quadrilateral 

dr=interval width. 

12 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

The quantity ¢r is tabulated as a. function of r for seven 

evenly spaced radii r per quadrila.tera.l and for two to 

eight additional radii to vertices and edges, depending 

upon the location of the site relative to the source. The 

accuracy of the calculated areas depends on the geometry of the source, location of the site, and on the number 

of radii used in the ca1culat.ions. 

To test the accuracy of th~ arc-area computations 

for a fixed value of itst, areas for a set of 585G quadrilateral-site pairs (used in a California study) were evaluated using itst=5, 7, 9, 11, and 13. For itst=7, the 

{"arc estimate" area - "true" area ) / "true" area 

exceeded 0.03 only five times and never exceeded 0.04; 

for itst=5, the error exceeded 0.04 seven times; for 

itst=9, 11 and 13, the largest error was 0.02. We concluded that itst = 7 gave sufficient accuracy. The number of evenly spaced radii is set by parameter itst = 7 in 

subroutine RRISK; itst may be increased be recompiling the program, but this should probably be done only 

under unusual circumstances (for example, for testing 

purposes), because using larger values of itst increases 

the computation time . 

Grid Spacing Requirements 

For Modeling Uncertainty 

SEISRISK III initially calculates expected acceleration occurrences at each site on a grid assuming that 

seismicity is uniform throughout a source zone. Then to 

determine acceleration occurrences at a site at (X11 YI), 

(a) (b) (c) 

·Site 

·Site 

(d) (e) 

Site Site Site 

Figure 9. Arcs of circles with centers at the site are used 

to approximate the area of a quadrilateral. Area = ¢>rdr 

in D. E. and F. Shaded area in E and F is the difference 

between the quadrilateral area and arc area. In E errors 

nearly cancel: in F the arc area is an overestimate. 
with earthquake location uncertainty included, SEISRISK III computes a weighted average of acceleration 

occurrences previously calculated at various sites assuming uniform seismicity. The fractional contribution 

or weight assigned to the accelerations calculated at a 

site at (Xi, Yk) depends on a and on the distance between (X1, Yt) and (Xi, Yk) (equation 2). 

The grid spacing or distance between sites is of major importance in these calculations. Accelerations at 

sites at distances greater than 2.0a- 2.5a from (XIJ YI) 

have low weights and little effect on the weighted accelerations at (X1 , YI). This means that if the spacing 

between sites is large relative to a, the only accelerations that will have much weight at (X1 , YI) are those 

originally calculated at (Xb YI) (for uniform seismicity), and, in this case, the weighted results will give a 

poor approximation to the true effects of earthquake 

location uncertainty. In tests to determine the maximum acceptable grid spacing, we found that we could 

use a spacing dx as large as a and still obtain reasonable results; in these tests, for dx =a results were quite 

consistent with those obtained using dx=a /2. 

In order to take into account accelerations at sites 

at distances 2a from every site of interest, SEISRISK 

III extends the gridded region to include additional sites 

(fig.lO). If the sites of interest are at row and column 

intersections on an ( n) ( n) grid, and the distance between rows and columns is dx, then the grid must be 

extended in each direction by at least nrc = 2.0a / dx 

rows and columns; that is, the calculations must be 

done at (n + 2nrc) 2 sites in order to model effects of 

earthquake location variability at each of the original n 

sites. The number of additional sites can be quite large 

if a is small relative to dx. For example, if dx = 10 and 

a=40, the number of sites increases from n2 to (n+16) 2 . 

Increasing the number of sites at which the calculations are performed increases the storage requirements 

as well as the computation time. Because acceleration 

occurrences calculated for uniform seismicity at various sites are used in calculating the weighted averages, 

we retain these occurrences in a two-dimensional array. 

This requires an array of dimension of at least (I I, J J), 

where I /=number of acceleration levels in tl1e reg array, 

and J }=the number of sites (including the sites added 

for earthquake location variability) at which the cakulation is performed (J J 2: (n+2nrc) 2 ). Furthermore, a. 

second array of the same size is set aside in which to 

accumulate the weighted acceleration occurrences. In 

SEISRISK III, the two arrays RAWBIN and SMOBIN 

for saving the accelerations are dimensioned (55, 1600). 

A total of 1600 sites (including the sites added for earthquake location variability) can be processed. ( Dimensions of arrays SMOBIN and RAWBIN can be adjusted 

to change the number of sites permitted; the parameter 

maxar must be set equal to the site dimension.) SavFigure 10. Sites of interest located at intersections of 

grid lines within darkened area. If the spacing between 

grid lines is dx and the standard deviation in earthquake 

location is u. the gridded region is extended in each direction by Rrc = 2.0u/dx (rounded to the next highest 

integer) rows and columns. 

ing results in two-dimensional arrays contrasts with the 

procedure in SEISRISK II, in which acceleration densities calculated at a single site were written onto disc 

{logical unit 3 or logical unit 4) rather than being saved 

in a large array. 

Using Sources with Different a Values 

SEISRISK III permits earthquakes associated with 

different source-zones to have different values of a. ·If 

a is the same for earthquakes in a number of zones, 

weighted averages of accelerations can be calculated 

after the calculations assuming uniform seismicity are 

completed for earthquakes in all of the zones. On the 

other hand, if a different a is associated with earthquakes in each source zone, the weighting (equation 2) 

of accelerations resulting from earthquakes in a zone 

must be performed after the calculations for uniform 

seismicity have been done for that zone. 

To instruct the program when to cakulate weighted 

averages, a value of an input. parameter num is entered 

for each source zone along with other data for that zone; 

num=98 indicates averaging should be performed for all 

accelerations accumulated since the averaging was last 

performed. The current value of a is used in this averaging, and one need enter a value of a only for source 

zones for which num = 98. Figure llA gives an abbreviated flowchart for the source zone computations. 

Model Implementation-Details and Requirements 13 
Next site: 

Calculate accelerations 

from current source &: 

earthquakes with uniform 

spatial distribution. 

Add accelerations for 

this source to those 

accumulated In RAWBIN 

for this site. 

each site into 

SMOBIN array. 

At each site: 

Averaae accelerations accumulated in RAWBIN 

for various sites to 

adjust for earthquake 

location uncertainty. 

Save in SMOBIN array. 

Compute statistics 

at this llme? 

yes 

Calculate probabilistic 

acceleration levels from 

loE-------------------1 accumulated accelerations. 

~rint histogram and (or) 

save levels tor plotting. 

no 

Figure 11. Abbreviated flowcharts for (a) source-zone 

calculations and (b) fault-rupture calculations. 

When considering using different values of u for 

earthquakes in different zones, we should keep the following considerations in mind. 

1. The number of rows and columns and the spacing 

between sites remains the same throughtout the analysis. The value of ~rc, the number of additional rows 

and columns by which the initial gridded area must be 

extended in each direction to account for earthquake 

location uncertainty, is determined by the largest value 

of u to be used for any source-zone. However, nrc is 

calculated before the source zone inputs are processed, 

and therefore, the largest value of u must be identified 

in the program initially. 

2. H considerably different values of u are assigned to 

earthquakes in different areas, the grid spacing must be 

sufficiently small to adequately represent the intended 

variability for those earthquakes with the smallest u. 

At the same time, a sufficient number of marginal sites 

(determined by nrc) must be provided to permit accurate calculations when u is large. Accomodating both 

these requirements can lead to evaluating accelerations 

at considerably more sites than necessary, both when u 

is small (because of an unnecessarily large value of nrc 

14 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

Calculate accelerations 

from ruptures along 

current fault. Subdivide 

magnitudes for smoothing 

accelerations from closestdistance ruptures 

Add accelerations for 

ruptures on this fault 

to those accumulated in 

RAWBIN for this site 

Calculate accelerations 

from ruptures along 

current fault. Spread 

closest-distance 

ruptures over a distance 

of TOTL/2 km parallel to 

fault in each direction 

in this case) and when u is large (because of the small 

increment dx relative to u). 

3. Allowing different values of u can lead to an anomalous ridge of seismicity along one side of a joint boundary of contiguous source zones. In general, we suggest 

that two source zones with different u values be separated by at least two to three times the difference in the 

u values. 

Fault Computations 

Rupture Length Computations 

A rupture of a magnitude m earthquake is defined 

to have length li(m) with probability p(i), where 

where a and b are constants. 

p(i) =probability of rupture of length li(m) 

/(i+l) 

p(i) = vk j exp(- x; )dz 

f(i) 

(21) 

(22) 
where f(i) is selected so that 71(i) is 

f(i) < rJ(i) ~ f(i+ 1); /(1) = -oo; f(nr + 1) = oo; 

nr 

LP(i) = 1. 

i=l 

nr=number of rupture lengths per magnitude. 

SEISRISK III allows either one rupture length or 

or five rupture lengths per magnitude depending on 

whether median rupture lengths only are assumed or 

rupture length variability is to be modeled in the length 

versus magnitude relationship. Thus 

{ 

0, if only median rupture length is used 

( .) = (one rupture length per magnitude), 

'1 t five values in the range - 2 ~ '1 ( i) ~ 2 

(five lengths per magnitude). 

(23) 

The values of a, b and Ut may be provided as input 

for the rupture-length magnitude relationship (equation 

21). H a, b and Ut are all zero or blank, previous values 

of these quantities are used. H no values are provided, 

default values of Bonilla and Buchanan (1970) are used: 

a=-1.085 b=0.389 Ut =0.52 

H a f:. 0 , b f:. 0, and Ut = 0 are input, a single rupturelength value is used. 

Magnitude Smoothing: 

Calculation of Gutenberg-Richter b-Value 

The Gutenberg-Richter b-value used to determine 

the number of earthquakes in an arbitrary magnitude 

interval is calculated from the number N( ml) of earthquakes at ·magnitude m1 and N(m2 ), the number at 

magnitude m 2 ; that is, 

b = __ 1_ 1 [ N ( m2) J ~m og N(ml) · (24) 

where m1 is the lowest magnitude used in the calculations and m2 = m1 + ~m (equation 4). However, if 

only a single magnitude is specified, neither the magnitude interval nor the b-value can be calculated, and 

hence default values of b and ~m are provided. (The 

default values are b = 1.0 and ~m = 0.6.) 

Distance Smoothing Option 

To instruct the program to do distance (rather than 

magnitude) smoothing, the distanced between the parallel faults must be input. 

{ 

d distance between faults if this is a set of 

totl = 

0 

"dummy" faults used to approximate 

a uniform field of faults. 

if individual fault. 

If totl = 0 (or is blank), the program assumes that 

the faults are individual faults, and the previously discussed magnitude (rather than distance) smoothing is 

performed. Distance smoothing over a range of ±d/2 

km for a single fault can be forced by setting totl = d 

for a single fault. The program cannot do both distance and magnitude smoothing. Figure 11B shows an 

abbreviated flowchart for the fault rupture calculations. 

Attenuation Function: Input Table, 

Interpolation, Scale Factor 

The program computes ground-motion ( acceleration) by interpolating in a table of acceleration as a 

function of magnitude and distance. Interpolations are 

linear in distance and linear in log acceleration and magnitude; in SEISRISK II, interpolations were linear in log 

magnitude rather than in magnitude. 

The table of ground-motions is read in for up to 

eight magnitudes and 20 distances. Velocity or some 

other ground-motion parameter may be substituted for 

acceleration. H the ground-motion values read in are 

in a range different from 0.02 to 1.0, the appropriate 

scale factor must be input to make the "box" levels 

compatible with the desired values. 

Halving the Magnitude Intervals 

The program requires that earthquake rates be entered for each source at equally spaced magnitude intervals. The magnitudes that are entered are assumed 

to be values at the centers of the intervals. H input 

paramater i nos = 1, the program will divide each magnitude interval into two parts and perform the calculations for magnitudes at the midpoints of the new intervals. Earthquake rates are calculated for each interval 

from the total earthquakes assigned to the original intervals using a Gutenberg-Richter b-value based on the 

number of earthquakes in the first two (original) magnitude intervals, as in equation 24. 

Output Options 

As stated earlier, we seek the acceleration that, 

with probability q, will be exceeded (or probability p = 

1- q will not be exceeded) during t years. We have seen 

that this is the value a of acceleration for which Ex( a) 

the yearly exceedance rate of a at the site is 

Ex( a) = _loge(P). 

t (25) 

Assuming yearly exceedance rates (as in equation 15) 

have been calculated for each of 55 acceleration levels, 

the solution a that corresponds to the specified (p, t) 

may be found by interpolating between the calculated 

values. 

Model Implementation-Details and Requirements 1 S 
Solutions are computed for up to 20 time periods 

t using rates derived for median accelerations only and 

also for the same set of times using rates that include 

acceleration variability. 

The results, including the acceleration histogram, 

may optionally be written to file 16 for later printing. 

Summary results may be saved on file 2 for future plotting. 

Accelerations from the earthquakes within a source 

area or along a set of faults are calculated for all sites 

before proceeding to the next source area or fault, with 

acceleration occurrences being added to those obtained 

previously for each site. Therefore, partial results may 

be examined at any stage. Output options are given for 

each source set by the input parameter iprint. 

iprint= 

-1: continue to next set; no print, no statistics 

+1: compute solutions a; write histogram and solutions to file 16, then continue to next set 

+2: same as +1 but in addition save the solutions 

on file 2 

+3: no write to file 16; save solutions on file 2. 

Accumulating Rates In 

Sucessive Runs 

Rates that have been accumulated in each acceleration interval for each site are saved in file 3 at the 

en~ of a run. These values may be read back into 

the 2-dimensional array SMOBIN for a subsequent run. 

(The saved rates are rates that have been adjusted for 

earthquake location uncertainty.) Doing successive runs 

makes it possible, for example, to accumulate accelerations from sources having different attenuation functions. The equator, maximum a, gridded region, and 

rows and columns containing sites of interest should not 

be changed from their original values in a run continuation. 

. { 1 

'tSW = O 

if this is a run continuation 

if this is a new run. 

REFERENCES CITED 

Bender, B., and Perkins, D. M., 1982, SEISRISK II, A computer program for seismic hazard estimation: U.S. Geological Survey Open-File Report 82-293, 103 p. 

Bender, B., 1984a, Incorporating acceleration variability into 

seismic hazard analysis: Bulletin of Seismological Society 

of America v. 74, p. 1451-1462. 

--1984b. Seismic hazard analysis using a finite fault rupture model: Bulletin of Seismological Society of America 

v. 74, p. 1899-1923. 

Bender, B., 1986. Modeling source zone boundary uncertainty 

in Seismic hazard analysis: Bulletin of Seismological Society of America v. 76, p. 329-341. 

16 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

Bonilla, M. G., and Buchanan, J. M., 1970, Interim report on 

world wide historic and surface faulting: U. S. Geological 

Survey Open-File Report, 32 p. 

Joyner, W. B., and Boore, D. M., 1981. Peak horizontal acceleration and velocity from strong motion records including 

records from the 1979 Imperial Valley, California earthquake. Bulletin of Seismological Society of America v. 

71, p. 2011-2038. 

Schnabel, P. B., and Seed, H. B., 1973, Acceleration in rock 

for earthquakes in the Western United States: Bulletin 

of Seismological Society of America v. 63, p. 501-516. 

Thenhaus, P. C., 1983. Summary of workshops concerning regional seismic source zones of parts of the conterminous 

United States, convened by the U.S. Geological Survey 

1979-1980, Golden, Colo. U.S. Geological Survey Circular 898, 36 p. 

APPENDIX A 

Details of Fault-Rupture Computations 

A fault consists of one or more connected line segments. A rupture may be wholly contained within one 

segment or overlap any number of fault segments. It 

may be as long as the fault, but cannot extend beyond 

the end of the fault. The distance from a rupture to 

a site is the distance from the site to the point on the 

rupture closest to the site. 

A conceptually simple approach for calculating accelerations at a site resulting from ruptures along the 

fault pretends that ruptures of length li(m) of magnitude m earthquakes occur only at a set of n distinct 

locations, with rupture centers evenly spaced at increments ~li(m) on the fault. In this case, the kth rupture 

center is at a distance di ( k) from one end of the fault 

where 

~ (k) = li(m) +(k-1)~li(m) (A1) 2 

and 

~li(m) = L-li(m), 1 ~ k~ n. (A2) n-1 

The fractional seismicity C ( i, m) associated with one 

such rupture is 

C( . ) _ p(m)p(i) f 

t,m- , n (A3) 

p(m)=rate of magnitude m earthquakes for the set 

of faults, 

/=length of this fault/total length of faults in set, 

p(i)=probability of rupture of length li(m), for 

either 1 or 5 rupture lengths per magnitude. 

For each rupture, the closest distance to the site is 

determined (fig.A1). For one such rupture of magnitude 

m and length li ( m), the acceleration interval is found 

that contains the acceleration produced by a rupture of 

this magnitude at this distance. The fractional seismicity G(i, m) is added to the accumulator reg(i), for the 

appropriate j. 
TSite 

I 

I 

I 

- - - ..l.F-=--=======:::;::-::--

A 

TSite 

I 

I 

I 

I 

~Site 

1\ 

I \ 

I \ 

I \ 1\/ I \ 

I \ 

:;:=====::t:~\ 

c 

Figure AI. Arrows indicate direction of increasing distance along the fault from the site. 

A Closest distance from site to fault is at left end of fault. 

Distance increases monotonically along the fault as one 

moves in the direction of the arrows. 

B Distance to site decreases as one moves along the fault 

from left end of fault to point A. then increases monotonically as one continues to other end. 

C Distance to site alternately decreases and increases as 

one moves along the fault. 

The computer time required may be substantial if 

the computation of nearest distance from the site to 

a rupture is done repeatedly for different rupture locations along the fault. Some transformations are made in 

SEISRISK III that simplify the site-to-rupture computations. The following describes these transformations. 

For an arbitrary fault segment 1 ~ i ~ n, let 

L(i) = J[x!(i) -x2(i)j2+ [y!(i) -y2(i)]2 {A4) 

= length of segment i, 

where [xl(i), yl(i)], [x2(i), y2 (i)] =coordinates of endpoints of segment i. 

The coordinates are transformed so that fault end points 

[xt(i), Yl(i)j--+ (0, 0) 

[x2(i), Y2(i)]--+ (L(i), OJ (A5) 

and Li1e site at (Fx, Fy) 

(Fx, Fy) --+ [X(i), P(i)] 

where 

X(i) = [Fx-x1 (i)] cos a+[Fy-yl(i)] sin a, 

P(i) = -[Fx-xt(i)] sin a+(Fy-y1(i)] cos a. (A6) 

=perpendicular distance from site to line i. 

. Y2(i) -yl(i) 

sma = L(i). 

After L(i), X(i) and P(i) have been determined, 

the fault segment may be regarded as located on the xaxis extending from {0,0) to [L(i), OJ and the site located 

at [X(i), P(i)] (fig.A2). 

If the fault consists of N connected line segments, 

SEISRISK III computes and saves L(i), X(i), P(i) for 

each segment, 1 ~ i ~ N. Assume an arbitrary rupture 

has one end point on segment i and the second on segment k, where 1 ~ i ~ k ~ N. The coordinates of the end 

points of the part of the rupture that lies on segment i 

(i ~i ~ k) are defined as [Rl(i), 0], [R2(i), 0] (fig.A3). 

0 ~ R!(i) ~ R2(i) ~ L(i) fori= i or i = k, (A7) 

Ri(i)=O,~(i)=L(i) for i<i<k. 

If a rupture extends over three or more fault segments, 

the rupture overlap with an interior segment is the entire segment. 

In the new system, the relative end and overlap 

points Rl(i), R2(i) are determined from the total distance 8 along the fault of the initial end point (start) 

of the given rupture and from the length of the rupture 

and fault segments, by calculating Rem(i) for each segment i, where Rem(J') is defined as 

i-1 

Rem(i)=8- E L(i), where L(O)=O. 

i=O 

The first segment containing the rupture beginning at 

8 is the segment i for which Rem(J') >0 and Rem(J' + 

1) ~ 0. (For this i, R1 (J') = Rem(J').) The segment 

containing the final end point of a rupture of length l is 

determined similarly by substituting ( 8-tl) for ( 8) above. 

If the rupture overlaps segment i, let D2(i) =square 

of the closest distance from the site to the rupture on 

segment i. If the rupture extends over several segments, 

the value of D2 ( i), that is the minimum for all such 

segments, is required. D2 (i) is defined as follows: 

{ 

[Rl(i)-X(i)]2+P(i)2 if X(i)<R1(i), 

D2 (i) = [X(i)-R2(i)]2+P(i)2 if X(i)>L(i), 

P(i)2 if Rl(i)~X(i)~~(i). 

(A8) 

Appendix A 17 
(Fx,Fy) 

I 

Figure A2. Fault segment and site. A original coordinates, B relocated so that segment lies on X-axis. A Original coordinates. Segment of fault extends from(xb y!) 

to (x2,Y2); site is at (Fx,F11 ). Perpendicular from site 

intersects extended fault line at (x3 , y3 ). B Relocated. 

Fault and site moved so that fault extends from (0.0) to 

(L, 0). Site moved to (X, P). 

Thus, having calculated and saved L(i), P(i), X(i), 

determining the (square of the) closest distance requires 

only a simple computation for each fault aegment that 

the rupture overlaps. 

Assigning the fractional seismicity (equation A3) 

associated with a rupture to a single acceleration interval results in chunks of seismicity in a set of intervals. This means that the calculated acceleration densities do not vary smoothly as a function of acceleration 

level. In some cases, SEISRISK III can avoid doing calculations for a distinct set of ruptures and instead directly calculate the fraction of ruptures of magnitude m 

earthquakes that produce accelerations in each interval, 

a.i-l <a ~,a,.. When the more direct approach cannot 

be employed, SEISRISK III calculates site-to-rupture 

distance for n ruptures of magnitude m earthquakes, 

evenly spaced along the fault (equation A1). H rupture 

i is at distance di from the site, and rupture i + 1 is 

at distance di-+1 from the site, SEISRISK III assumes a 

fraction 1/n of the possible ruptures are in the distance 

range di ~ d ~ di+l (assuming di ~ di+I). It then finds 

the accelerations a(i) and a(i+ 1) produced by magnitude m ruptures at distances ~ and di+I· Assuming 

a(i + 1) ~ a; ~ aft-1 ~ a(i), the fractional seismicity assigned to the acceleration interval a;~ a~ ai+I is given 

by 

. aj-t-1 -a; C(~, m) (") (" 1), a~ -a~+ 

(A9) 

where G(i, m) is defined as in equation A3. 

We next consider two cases in which a more direct 

calculation of fractional accelerations can be done for 

each acceleration interval. We use the notation and 

definitions already developed. 

Case 1. X(i) < 0 for all 1 ~ i ~ n (or L(i) < X(i) 

for all 1 ~ i ~ n). The point on the rupture nearest the 

18 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

(X,P) 

- _1_(0,0-)-------:-:-(L,O) 

B 

site is at one end of the rupture with the distance to 

the site increasing (or decreasing) monotonically along 

the rupture. 

If X(i) < 0 for all i, for any rupture overlapping 

segments j through k, where j ~ k, then [ R1 (j), 0] is 

the point closest to the site. Similarly if X(i) > L(i) for 

all i, then [R2(k), 0] is the point of the rupture closest 

to the site. 

Assume X(i) < O, for all i. (L(i) < X(i) could be 

considered in a symmetrical argument.) In this case the 

rupture point nearest to the site is RI(i) for the lowest 

i containing the rupture. Recall that, for a rupture of 

length l, rupture centers may be located at any distance 

de along the fault for which 

l l -<d <L-- 2- C- 2 (A10) 

or, equivalently, the lowest end point of a rupture may 

be at any distance 6 for which 

O~S~L-l. 

The acceleration produced by a rupture of magnitude m and length li(m) decreases as the rupture 

moves along the fault from a maximum acceleration 

amax [li{m)] for a rupture with closest end point at 6 =0 

to a minimum acceleration amin[li(m)] for a rupture 

with closest end point at 6=L-l(m). Therefore, for an 

acceleration interval such that 

a corresponding interval exists along the fault 

hm(i) ~S~hm(j-1) 

for magnitude m earthquakes in which (closest) end 

points of ruptures producing these accelerations may 

be located. Because these ruptures are at a distance 

dm (i) ~ d ~ dm (i -1) from the site, 

hm(j -1) = V dm(i -1)2 -P(i)2+X(i). (A12) 
Ruptures at distance along the fault hm(j) ~ 5 < 

hm (j - 1) therefore contribute to accelerations in the 

range ai-l <a~ ai at the rate 

[ ( )] ( ) hm(j-1)- hm(i) (") 

rat ai, li m =p m L-li(m) f P z (A13) 

where /, p(i) and p(m) are defined as in equation A3. 

For n rupture lengths, given that 

l1 (m) > l2{m) ... > ln(m), (A14) 

a corresponding section of fault 0 ~ 5 ~ L-li(m) may 

contain closest end points of ruptures and 

L-ln(m) >L-ln-dm) > ... >L-h(m). (A15) 

Now define 

h ( ") { hm(j) m,i J = L-li(m) 

if hm {i) < L -li ( m) 

otherwise, 1 ~ i ~ n. 

The total contribution to accelerations ai-l < a ~ ai is 

the sum over rupture lengths for magnitude m ruptures: 

rat(a) = p(m) t hm,i(i- 1)-hm,i(i) f p{i) {A16) 

i=l L-li(m) 

where hm,i (0) = 0. 

Thus, one need only compute one set of distances 

hm (i) per magnitude along the fault rather than repeating the process for each rupture length, so that 

including rupture length variability in the calculation 

does not substantially increase the computation time. 

(0,0)----tt ILuO) 

(0,0)---------(L2 ,0) 

R1(1) 

(0,0)----+-] --(Ll,O) 

rupture 

----4t R2(1) R1(2)---------R2 (2) 

R1(3)-- -+-] --

R2(3) 

Figure AS. Rupture overlaps three fault segments. Each 

segment is relocated to lie on the X-axis. End points of 

the part of the rupture on segment i are [Rl(i), 0] and 

[R2(i), 0]. 

Case 2. The distance from the fault to the site decreases monotonically as one moves a distance de along 

the fault to some point on segment k [(Ik,O) in the new 

system] and then increases steadily as one continues to 

the other end. This corresponds to: 

X(i)<O for 1~i~k-1, 

k-1 

X(k)=lk=de- L L{i), 

i+l 

X(i)>L(i) for k+1~i~n. 

In this case, all ruptures of length li(m) with midpoints at a distance 5 along the fault 

d li(m) £ d li(m) ---<o< +-- e 2 e 2 (A17) 

overlap the point de and hence are regarded as being at 

the same closest distance from the site. The acceleration at the site produced by these magnitude m ruptures 

is the acceleration at distance P( k), the perpendicular 

distance to the kth segment. Because ruptures cannot 

extend past the end of the fault, the location of possible 

rupture centers is restricted to a part of the fault and 

the contribution at P( k) becomes 

rat[a(j), l;(m)J = p(m)f t Du, ~~;:{~j (m) p(i), 

i=l ' 

(AlB) 

where 

D ( ) _ . [d + li { m) L _ ~ { m) J u; m - mm e 2 , 2 , {A19) 

D ( ) _ [d _ li(m) li(m)] l; m - max e 2 , 2 . 

Now the bookkeeping becomes more complicated, but 

we can look at all ruptures of length li(m) with centers 

at distance 5 on the fault 

li(m) < 5 < L _ li(m) 

2 - - 2 (A20) 

and remove those already accounted for, those with centers 

d _ li(m) < £ d li(m) e 2 _u < e+ 2 . (A21) 

The remaining ruptures may be separated into two sets: 

those (if any) with centers 

li(m) < 5 <d _ ~{m) 2 - e 2 (A23) 

and those for which 

d + li ( m) < 5 < L _ li( m) . e 2 - 2 

Appendix A 19 
The last two intervals correspond to ruptures whose 

closest end points to the site are at 6, where 

(A24) 

for one set of ruptures, and at 

(A25) 

for the other set. These two sets of ruptures do not 

overlap and the fault may be regarded as consisting of 

two distinct segments in which the distance from the 

rupture to the site decreases (or increases) monotonically as one moves along the fault. Hence, each of these 

segments may be treated separately as in Case 1, and 

the accelerations resulting from ruptures along these 

segments may be added to those already accounted for. 

APPENDIX B 

Sample Inputs and Outputs 

Sample inputs and outputs are provided in this appendix. An explanation of the input variables (file 15) 

and formats is provided at the beginning of the program 

listing in Appendix C. 

The output file (file 16) gives ground-motion statistics for each site. The site coordinates (longitude, latitude) in decimal degrees appear first, followed by the 

shortest distance from the site to the closest fault. Then 

two arrays are given, one for zero attenuation variability and one for attenuation variability Ua. The columns 

20 Seisrisk Ill: A Computer Program for Seismic Hazard Estimation 

in each array give for the ground-motion (g.m.) in the 

ith row: 

occlyr: occurrence rate per year of earthquakes 

producing ground motions at the site in the 

range g.m.(i-1) < g.m.~ g.m.(i); 

exclyr: exceedance rate per year of g.m.(i); 

r( events): the average number of events required 

to produce an exceedance of g.m.(i), 

( ) 

total yearly events r events I (") ; exc yr of g.m. 1 

r(yrs): return period in years, 

1 

r(yrs) = I ("). exc yr of g.m. 1 

Following the arrays is "total yearly events," the 

total number of earthquakes from all sources in this 

run. This total should be nearly the same at all sites; 

slight differences between sites are due to inaccuracies in 

calculating fractional areas and assigning proportional 

seismicity to these areas. Large differences probably 

indicate an error in the computation. 

The concluding lines for each site read 

"xxx ext prob=yyy for zzz years", 

that is, xxx is the extreme probability that groundmotion yyy will not be exceeded in zzz years. 

Included also in this appendix are a listing of a program to read (binary) output file 02 and the output of 

that program for the sample problem. File 02 contains 

only the summary information needed for mapping purposes: site coordinates (long, lat) and ground-motions 

calculated at the specified probability level for the times 

of interest (for both zero-attenuation variability and one 

value of ua)· 
f 

Q. 

~· 

= 

~ 

File 15--Inputs 

SEISRISK III Sample Problem 

0 20. 

. 90 3 10 50 250 

1. 0.50 

123.0 38.0 120.0 38.0 

123.00 39.00 120.00 

5 6 5 6 

1 

3 

123. 39. 123. 38 

6 12 

'watashh79' 

3.22 .74 

6. 43 . 64 

16. 09 . 49 

32. 18 . 36 

64. 3 . 21 

96. 5 . 12 

160.9 .045 

321. 8 . 013 

643.0 . 0034 

1288.0 .00085 

2570.0 .00021 

5140.0 .0001 

00 100. -1 

3 1 1 

8. 5 

. 73 

. 62 

. 43 

. 28 

. 14 

. 07 

. 025 

. 0076 

. 0019 

. 00047 

. 00012 

. 00006 

zn01 

122.20 37.22 122.00 37.32 

123.85 38.97 123.68 39.05 

124.12 39.32 123.89 39.40 

. 1 06 . 281 . 7 4 1. 97 

6.100 5. 500 4.900 4.300 

00 65.2 -1 zn02 

2 1 1 

120.05 34.09 119.26 34.20 

120.40 33.50 119.77 33.39 

37.00 

7.6 

. 67 

. 53 

. 32 

. 17 

. 06 

. 03 

. 015 

. 0026 

. 00065 

. 00016 

. 00004 

. 00002 

.031 .073 . 172 .405 .955 2.248 

7.300 6.700 6.100 5. 500 4.900 4.300 

98 1. 0 -1 zn03 20. 

7 1 2 

120.40 34.68 

120. 56 34. 93 

120.78 35.20 

121. 00 35. 48 

121.51 36.00 

121.82 36.45 

122. 00 36. 73 

2 2 2 

120.00 

120.28 

120.38 

121. 00 

121.47 

121. 80 

122.21 

34. 78 

35. 00 

35.39 

36.00 

36. 50 

36.80 

37. 11 

121.80 36.80 121.77 36.88 

122.21 37.11 122.20 37.22 

.0198 .0641 .2074 .6712 

6.100 5. 500 4.900 4.300 

. 100 

6.6 

. 45 

. 36 

. 19 

. 09 

. 035 

. 0138 

. 005 

. 0012 

. 0003 

. 00007 

. 00002 

. 00001 

99 1.0 -1 ft01-1.085 .389. 50 

12 1 1 

117.37 34.17 

119.26 34.90 

122.10 37.27 

. 0021 . 0057 . 015 

8. 500 7.900 7.300 

118. 09 

119. 54 

123. 77 

. 040 

6. 700 

34. 59 

35. 06 

39.01 

118.67 

121. 27 

124.01 

. 100 

5. 6 

. 195 

. 135 

. 052 

. 02 

. 0051 

. 0023 

. 00083 

. 00021 

. 00005 

. 00001 

. 00001 

. 00001 

34. 79 

36.66 

39.36 

5.2 

. 072 

. 047 

. 02 

. 0052 

. 0013 

. 00042 

. 0001 

. 00003 

. 00001 

. 00001 

. 00001 

. 00001 

119. 03 

121. 69 

124. 16 

4.2 

34.85 

36.94 

39.93 

00 10. +2 

4 1 4 

122. 40 38. 10 

3 2 4 

122.60 38. 68 

5 3 4 

122. 16 38. 16 

123. 40 39.98 

3 4 4 

122. 50 39. 10 

.0072. 0200 .057 

7.900 7. 300 6.700 

99 

99 

ft02-1.085 .389. 50 

122.90 38.66 123.38 39. 10 123.94 39. 79 

123. 10 39.20 123. 74 39.95 

122. 46 38. 80 122. 75 38.98 123. 18 39. 56 

122. 80 39.25 122.95 39.43 
~ 

~ 

VI 

~ 

~· 

~ 

> 

~ i 

li 

.. 

., 

i 

3 

Q' 

.. 

~ 

~· 

;::;· 

:z::: 

=: 

Ill 

a. 

~ 

[ 

~· 

File 16--0utputs 

SEISRISK III Sample Problem 

isw=O: new run--no previous results included 

extreme probability 0.900 

for exposure times <years> 10 50 250 

scale factor for ground motion "box" levels= 1.00 

coordinates input in decimal degrees 

coordinates are printed in decimal degrees 

variability in attenuation, sigma= 0. 50 

grid oriented parallel to great circle thru < 123.00, 38. 00), < 120.00, 38.00) 

corners of gridded area-upper left= 123.00, 39.00 

lower right= 120.00, 37.00 

longitude increment= 0. 1000 (decimal degrees) 

latitude increment = 0. 1000 (decimal degrees> 

gridded region contains 27 rows. 30 cols' including border 3 rows and cols 

for this run begin at row 5 end row 6, begin col 5 end col 6 

new coordinates (km) gridded area 

upper left= 294.73 0.00; lower right= 144.70 -111.30 

sites are also located on 1 line
 

3 sites per line 

line 1 end points at 123.000, 39.000 and 123.000, 38.000 

attenuation function watashh79 

magnitude 

dist<km> 

3.22 

6.43 

16.09 

32. 18 

64.30 

96. 50 

160.90 

321. 80 

643.00 

1288.00 

2570.00 

5140.00 

8.50 

0. 74000 

0.64000 

0.49000 

0.36000 

0.21000 

0. 12000 

0.04500 

0.01300 

0.00340 

0.00085 

0.00021 

0.00010 

7.60 

0.73000 

0.62000 

0.43000 

0.28000 

0. 14000 

0.07000 

0.02500 

0.00760 

0.00190 

0.00047 

0.00012 

0.00006 

yrnoc= 100. iprint=-1 for area 

122.200 37.220 122.000 

123.850 38.970 123.680 

124.120 39.320 123.890 

nr of levels of seismicity = 4 

zn01 beta= -1.6319 

before normalizing to rate/year 

6.60 

0.67000 

0. 53000 

0.32000 

0. 17000 

0.06000 

0.03000 

0.01500 

0.00260 

0.00065 

0.00016 

0.00004 

0.00002 

zn01 

37. 320 

39.050 

39. 400 

5.60 

0.45000 

0. 36000 

0. 19000 

0.09000 

0.03500 

0. 01380 

0.00500 

0.00120 

0.00030 

0.00007 

0.00002 

0.00001 

occurrences• 1.970000 0.740000 0.281000 0.106000 

earthquake rate I year 

occurrences= 0.019700 0.007400 0.002810 0.001060 

magnitudes= 4.30 4.90 5. 50 6.10 

5.20 

0. 19500 

0. 13500 

0.05200 

0.02000 

0.00510 

0.00230 

0.00083 

0.00021 

0.00005 

0.00001 

0.00001 

0.00001 

4. 20 

0.07200 

0.04700 

0.02000 

0.00520 

0.00130 

0.00042 

0.00010 

0.00003 

0.00001 

0.00001 

0.00001 

0.00001 

zn01 area= 5457. sq km, rate/sq km= 0.36101E-05 for mags 4.00- 4.60 

yrnoc= 6~ iprint=-1 for area zn02 

120.050 34.090 119.260 34.200 

120.400 33.500 119.770 33.390 

nr of levels of seismicity = 6 

zn02 beta= -1.4268 

before normalizing to rate/year 

occurrences= 2.248000 0.955000 0.405000 0.172000 0.073000 0.031000 

earthquake rate I year 

occurrences= 0.034479 0.014647 0.006212 0.002638 0.001120 0.000475 

magnitudes= 4.30 4.90 5. 50 6.10 6. 70 7.30 

zn02 area= 5212. sq km, rate/sq km= 0.66157E-05 for mags 4.00- 4.60 

yrnoc= 1. iprint=-1 for area zn03 

smooth for earthquake location uncertainty sigma= 20.0 

120.400 34.680 120.000 34.780 

120.560 34.930 120.280 35.000 

120.780 35.200 120.380 35.390 

121.000 35.480 121.000 36.000 

121.510 36.000 121.470 36.500 

121.820 36.450 121.800 36.800 

122.000 36.730 122.210 37.110 

121.800 36.800 121.770 36.880 

122.210 37.110 122.200 37.220 nr of levels of seismicity = 4 

zn03 beta= -1.9574 

earthquake rate I year 

occurrences= 0.671200 0.207400 0.064100 0.019800 

magnitudes= 4.30 4.90 5. 50 6.10 

zn03 area= 10075. sq km, rate/sq km= 0.66621E-04 for mags 4.00- 4.60 

yrnoc= 1. iprint=-1 for area ft01 

Distance between dummy faults= 0.0 

fault 1 of 1 

117.37 34. 17. 118.09 

119.26 34.90, 119.54 

34. 59. 

35. 06. 

122. 10 37. 27, 123. 77 39. 01o 

118. 67 

121.27 

124.01 

34. 79. 

36. 66. 

39. 36, nr of levels of seismicity = 4 

ft01 beta= -1.6347 

earthquake rate I year 

occurrences= 0.040000 0.015000 0.005700 

magnitudes= 6.70 7.30 7.90 

fault rupture length parameters al= -1.085 

yrnoc= 10. iprint= 2 for area ft02 

Distance between dummy faults= 0.0 

fault 1 of 4 

122. 40 38. 10. 122. 90 38. 66. 

fault 2 of 4 

122. 60 38. 68. 123. 10 39. 20. 

fault 3 of 4 

122. 16 38. 16. 122. 46 38. so. 123. 40 39. 98, 

fault 4 of 4 

122. 50 39. 10, 122. 80 39. 25. 

nr of levels of seismicity = 3 

ft02 beta= -1.7455 

before normalizing to rate/year 

occurrences= 0.057000 0.020000 0.007200 

earthquake rate I year 

occurrences= 0.005700 0.002000 0.000720 

magnitudes= 6.70 7.30 7.90 

fault rupture length parameters al= -1.085 

0.002100 

8. 50 

bl= 0.389 sigl= 

123.38 39. 10. 

123. 74 39. 95, 

122.75 38. 98. 

122. 95 39. 43. 

bl= 0.389 sigl= 

119. 03 

121. 69 

124. 16 

0. 50 

123.94 

123. 18 

0. 50 

34. 85, 

36. 94, 

39. 93, 

39. 79, 

39. 56. 
t 

>e• 

= 

h) 

~ 

SEISRISK III Sample Problem 

site at long 122.480, lat 38.605 

shortest dist to fault• ~ 128 km 

zero attenuation variability 

g.m. occ/yr exc/yr r<eventsl r<yrsl 

0. 02 1. 08967 0. 03459 32. 5 28. 9 

0.04 0.01039 0.02421 46.4 41.3 

0.06 0.00559 0.01862 60.4 53.7 

0.08 0.00627 0.01235 91.0 81.0 

0.10 0.00241 0.00994 113.2 100.7 

0.12 0.00220 0.00774 145.3 129.3 

0.14 0.00092 0.00682 164.8 146.6 

0.16 0.00112 0.00570 197.2 175.4 

0.18 0.00101 0.00469 239.9 213.4 

0.20 0.00077 0.00392 286.8 255.1 

0.22 0.00083 0.00309 363.4 323.2 

0.24 0.00048 0.00261 430.7 383.1 

0.26 0.00029 0.00232 484.3 430.7 

0.28 0.00021 0.00211 533.3 474.3 

0.30 0.00018 0.00193 583.9 519.3 

0.32 0.00013 0.00180 625.8 556.6 

0.34 0.00015 0.00164 684.4 608.7 

0.36 0.00018 0.00146 769.8 684.7 

0.38 0.00013 0.00134 842.0 748.9 

0.40 0.00008 0.00126 893.0 794.3 

0.42 0.00007 0.00119 943.4 839.1 

0.44 0.00024 0.00095 1180.6 1050.1 

0.46 0.00021 0.00074 1518.3 1350.4 

0.48 0.00018 0.00056 2002.5 1781.2 

0.50 0.00014 0.00042 2672.0 2376.7 

0. 52 0.00014 0.00028 4015.7 3571.9 

0. 54 0.00009 0.00019 5936.6 5280.4 

0. 56 0.00006 0.00013 8550.0 7604.9 

0. 58 0.00012 0.00001 99999.9 99999.9 

0.60 0.00001 0. 00000 99999.9 99999.9 

0.62 0.00000 0.00000 99999.9 99999.9 

0.64 0.00000 0. 00000 99999.9 99999.9 

0.66 0.00000 0.00000 99999.9 99999.9 

0.68 0.00000 0.00000 99999.9 99999.9 

0.70 0.00000 0.00000 99999.9 99999.9 

0.72 0.00000 0.00000 99999.9 99999.9 

0.74 0.00000 0.00000 99999.9 99999.9 

0.76 0.00000 0.00000 99999.9 99999.9 

0.78 0.00000 0.00000 99999.9 99999.9 

0.80 0.00000 0.00000 99999.9 99999.9 

0.82 0.00000 0.00000 99999.9 99999.9 

0.84 0.00000 0.00000 99999.9 99999.9 

0.86 0.00000 0.00000 99999.9 99999.9 

0.88 0.00000 0.00000 99999.9 99999.9 

0.90 0.00000 0.00000 99999.9 99999.9 

0.92 0.00000 0.00000 99999.9 99999.9 

0.94 0.00000 0.00000 99999.9 99999.9 

0.96 0.00000 0.00000 99999.9 99999.9 

0.98 0.00000 0.00000 99999.9 99999.9 

1.00 0.00000 0.00000 99999.9 99999.9 

1.07 0.00000 0.00000 99999.9 99999.9 

1.14 0.00000 0.00000 99999.9 99999.9 

1.21 0.00000 0.00000 99999.9 99999.9 

1.27 0.00000 0.00000 99999.9 99999.9 

1. 34 0.00000 0.00000 99999.9 99999.9 

1.41 0.00000 0.00000 99999.9 99999.9 

1.48 0.00000 0.00000 99999.9 99999.9 

1. 55 0.00000 0.00000 99999.9 99999.9 

1.62 0.00000 0.00000 99999.9 99999.9 

1.69 0.00000 0.00000 99999.9 99999.9 

total yearly events 1. 12427 

zero attenuation variability 

0.900 ext prob = 0.095 for 10 years 

0.900 ext prob = 0.280 for 50 years 

0.900 ext prob = 0. 500 for 250 years 

variability in atten. sigma= 0. 50 

g.m. occ/yr exc/yr r<eventsl r(yrsl 

0. 02 1.08964 0. 03462 32.5 28.9 

0.04 0.01057 0.02405 46.7 41.6 

0.06 0.00611 0.01794 62.7 55.8 

0.08 0.00403 0.01390 80.9 71.9 

0.10 ~00281 0.01109 101.4 90.2 

0.12 0.00203 0.00906 124.1 110.4 

0.14 0.00152 0.00754 149.0 132.6 

0.16 0.00116 0.00638 176.2 156.7 

0.18 0.00092 0.00546 205.8 183.0 

0.20 0.00074 0.00473 237.9 211.6 

0.22 0.00060 0.00412 272.6 242.5 

0.24 0.00050 0.00362 310.3 276.0 

0.26 0.00042 0.00320 351.2 312.4 

0.28 0.00036 0.00284 395.4 351.7 

0.30 0.00031 0.00254 443.3 394.3 

0.32 0.00027 0.00227 495.2 440.5 

0.34 0.00023 0.00204 551.5 490.5 

0.36 0.00020 0.00184 612.4 544.7 

0.38 0.00018 0.00166 678.6 603.6 

0.40 0.00016 0.00150 750.4 667.4 

0.42 0.00014 0.00136 828.3 736.7 

0. 44 0.00013 0.00123 912.9 812.0 

0.46 0.00011 0.00112 1004.7 893.7 

0.48 0.00010 0.00102 1104.5 982.4 

0. 50 0.00009 0.00093 1212.8 1078.8 

0. 52 0.00008 0.00085 1330.5 1183.4 

0. 54 0.00007 0.00077 1458.2 1297.1 

0. 56 0.00007 0.00070 1597.0 1420.5 

0. 58 0.00006 0.00064 1747.6 1554.4 

0.60 0.00006 0.00059 1911.0 1699.8 

0.62 0.00005 0.00054 2088.4 1857.6 

0.64 0.00005 0. 00049 2280.8 2028.7 

0.66 0.00004 0.00045 2489.4 2214.2 

0.68 0.00004 0.00041 2715.5 2415.4 

0. 70 0.00003 0.00038 2960.5 2633.3 

0.72 0.00003 0.00035 3225.9 2869.4 

0. 74 0.00003 0.00032 3513.3 3125.0 

0. 76 0.00003 0.00029 3824.4 3401.7 

0.78 0.00002 0.00027 4160.9 3701.0 

0.80 0.00002 0.00025 4524.8 4024.7 

0.82 0.00002 0.00023 4918.3 4374.7 

0.84 0.00002 0.00021 5343.4 4752.8 

0.86 0.00002 0.00019 5802.5 5161.2 

0.88 0.00002 0.00018 6298.2 5602.1 

0.90 0.00001 0.00016 6833.2 6077.9 

0.92 0.00001 0.00015 7410.3 6591.3 

0.94 0.00001 0.00014 8032.6 7144.8 

0.96 0.00001 0.00013 8703.4 7741.4 

0.98 0.00001 0.00012 9426.1 8384.3 

1.00 0.00001 0.00011 10204.4 9076.5 

1.07 0.00003 0.00008 13359.0 11882.4 

1.14 0.00002 0.00006 17403.6 15480.0 

1.21 0.00001 0.00005 22565.4 20071.3 

1.27 0.00001 0.00004 29122.6 25903.7 

1.34 0.00001 0.00003 37417.7 33282.0 

1.41 0.00001 0.00002 47868.3 42577.4 

1.48 0.00001 0.00002 60977.2 54237.5 

1. 55 0.00000 0.00001 77362.9 68812.0 

1.62 0.00000 0.00001 97755.8 86951.0 

1.69 0.00000 0.00001 99999.9 99999.9 

total yearly events 1. 12426 

variability in atten. sigma=O. 50 

0. 105 for 10 years 

0.334 for 50 years 

0.676 for 250 years 

ratio 250 yr 0.900 extreme value to 10 yr val= 5.28 6.43 

·program to read summary file02 created by seisrisk iit 

file 02 contains for each site the levei of ground motion 

that has probability "perct" of not being exceeded during 

"Jtim" vears for each of "ntim" values of "Jtim". 

ground motion values are given first for "zero variability in 

attenuation" for all times "Jtim"• then for "variability 

in attenuation sigma = sd" for all values of "Jtim". 

first record contains: irow1, irow2. icoll. icol2 = first and last rows. first and last 

columns in gridded area containing sites for which 

ground motions were calculated in seisrisk ii. 

indvpt = total number of sites on lines 

ntims a number of exposure times 

Jtims = exposure times in years 

sd =standard deviation in attehuatiow variability. 

perct: ground motion calculateJ has p~obability "perct• 

of not being exceeded duri~g "Jtim" years 

all records after first: xlong.ylat =longitude, latitude of site in decimal degrees 

sol(il =ground motion at <xlat.ylong) corresponding to 

"per~t" probability and Jtim(k) years where 

i~1.2 •..... ntim. k=i for zero attenuation variability 

i = ntim+1, ntim+2 ..... 2*ntim. k=i-ntim 

for variability in attenuation sigma = sd 

dimension sol(48),Jtim(24l 

·read (2) irow1, irow2, icol1, icol2, indvpt. ntim. < Jtim< Jvl. JV=1· 

&ntiml.perct.sd 

ntim2=2•ntim 

write(44. 10), ( Jtim( JV), JV=1, ntim)', ( Jtim< Jvl. JV=1, ntiml 

10 format(30x'exposure times'/' long lat '6i9) 

format for printout will depend on scale used for solution 

20 read (2, end=40l xlong. ylat, (sol ( Jtl. Jt=1, ntim2l 

write (44. 30), x long. y lat. <sol (Jtl. Jt=t. ntim2) 

30 format(1h 2f7.3.6f9.3) 

go to 20 

40 call exit 

end 

The following are the summarized results that were written 

to file 02 for the sample problem 

no acceleration variability sigma•. 5 in log acceleration 

long 

122.48 

122.35 

122.48 

122.35 

123.00 

123.00 

123.00 

lat 

38.61 

38.61 

38.51 

38.p1 

39.00 

38.50 

38.00 

10yr 50yr 250\jr 10yr 50yr 250yr 

0.095 

0. 079 

! 0. 104 

0.086 

0. 129 

0. 239 

0.219 

0.280 

0. 224 

0.277 

0.224 

0.351 

0. 398 

0.443 

0. 500 

o. ·n2 0. 418 

0. 714 

0. 558 

0.454 

0.499 

0. 105 

0. 087 

0. 115 

0. 096 

0. 138 

0. 210 

0. 203 

0.334 

0.295 

0.333 

0.301 

0. 422 

0. 520 

0. 570 

0. 676 

0.871 

0.614 

0. 802 

0.784 

0.834 

0.925 
N 

~ 

fl) 

!!. 

Ill 

:::!. 

fQ. c 

c .. 

> 

n 

= 3 

~ ;-

.. 

., 

a 

CJCI 

ill 

3 

;-

.. 

fl) 

!!. 

Ill 

it. 

1"1 

:z: 

~ Dol 

a. 

m 

Ill 

~ 

3 

~. 

= :a 

SEISRISK III 

computer program for determining ground motions 

for use in seismic hazard mapping 

******Allows locations of future earthquakes in seismic 

source zones to have a normal distribution******* 

SEISRISK III is a rev1s1on of SEISRISK II, documented in Open 

File Report 82-293. 

See model description of SEISRISK III and additional 

documentation in USGS Bulletin 1772 

files: 

file 15 

file 16 

file 02 

file 03: 

inputs 

outputs 

summary for plots--see text 

accelerations saved on file 03 for use in next 

run if next run is a continuation 

inputs 

1. title-up to 80 characters <free field> 

2. isw, sigmax <free field) 

isw=O new run <usual case> isw=l restart or continue 

from previous run 

sigmax=maximum standard deviation (km> for earthquake 

location variability in seismic source zones for 

this run. 

3. prob, ntims, Jtim<l>, Jtim<2> .... Jtim<ntims> 

acceleration is sought for which there is probability prob of 

not being exceeded in Jtim< 1 >, Jtim<2>, ... Jtim<ntim) 

years 

prob = extreme probability in decimal 

ntims = number of times for which calculation is done 

Jtim<i> = durations <years> for which extreme motions are to be 

calculated at the prob probability level, 1 <= i <= ntims 

4. scale, dsw, sd, inos (free field) 

scale=scaling factor for ground motion boxes 

scale=1 for motions .02 to 1.0--scale accordingly 

dsw: =1 if inputs are degrees and minutes 

=0 if inputs are decimal degrees 

sd=standard deviation in log acceleration 

for acceleration variability around mean value 

inos=1: divide magnitude interval in half and do calculations 

at twice as many magnitudes <Assumes Gutenberg-Richter 

b-value is same for all intervals.> 

inos= any number other than 1: do calculations for original 

magnitude intervals 

5. x1, yl, x2, y2 <long, lat> in decimal degrees (free field> 

transform great circle thru <x1.y1>, <x2.y2>, <O,O> 

to equator for new coordiate system 

6. fl1,phl,fl2,ph2,phinc 

(fl1.ph1> upper left <fl2,ph2) lower right <long, lat> 

corners of seismic felt region for risk computation 

rectangular region with sides parallel to arc 

of great circle thru <x1,y1), (x2,y2> (0,0) --defined by 5. 

other two sides perpendicular and thru given end points 

c 

c 

c 

c 

c 

c 

c 

c 

c 

+ (fl1, ph1) 

1 

<x1.y1> 

1 

(fl2.ph2>+ < x2, y2> 

area within dots represents felt region 

'1' s represents line JOining <x1.y1> <x2.y2) 

phinc <long, lat> increment in degrees (in new coordinate 

system) for which risk is to be computed 

7. irow1, irow2, icol1, icol2 (free field) 

starting and ending rows and columns for this run 

accelerations computed for sites in these rows 

and columns--may be subset of seismic felt region defined 

in input 6. Also irow1=irow2=icol1=icol2=0 permitted 

if only individual sites on lines (input B> are selected. B. indv <free field> 

number of line segments containing individual sites at 

which acceleration is to be computed; zero if 

only fixed grid is used. 

If indv > 0 read next inputs, otherwise skip to input 9. nvs (free field) 

number of sites per line segment 

xel,ye1, xe2.ye2 <free field--indv lines• one pair per line) 

end points (long, lat) of line segment: sites will be 

evenly spaced on line in new coordinate system 

NOTE: NO SMOOTHING OF ACCELERATIONS AT SITES ON LINE FOR 

ACCELERATIONS FROM EARTHQUAKES WITHIN AREAL SOURCE ZONES. 

CALCULATIONS FOR UNIFORM SEISMICITY ONLY AT THESE SITES. 

9. Jent.mdis (free field> 

Jent=no of magnitudes for which acceleration is 

tabulated as a function of distance <maximum 8) 

mdis=number of distances in attenuation table <maximum 20> 

10. nam, tm < Jent values of tm> 

nam=identifier up to 10 characters-for example 'schn-seed' 

identifies attenuation curves used 

tm=magnitude for which table of distance versus acceleration 

values follows (free field) 

*******magnitudes must be in descending order*** 

11. rtab(i), <atab(i. J), J=1, Jent> <mdis lines, i=1,mdis) 

r-

;· 

!:!'. 

:::s 

OQ 

0 

-ott 

n 

0 

3 

~ 

c 

.. 

~ 

.. 

"'C 

.. 

0 

OQ 

.. 

r:» 

3 

(I) 

m 

(I) 

:::tJ 

(I) 

" 

rtab<i>=ith tabular distance in kilometers from earthquake source 

atab(i, J)=mean peak acceleration <or other ground motion parameter 

at i th distance for J th magnitude. 

source area inputs: 

each source deck - -

1s. num, yrnoc, iprint, totl, dumid, als, bls, sigls 

format< i2, flO. O, i2, flO. O, a4, 2f6. 2, f5. 2> 

num=O for seismic source areas if no averaging for earthquake 

location uncertainty is to be done at end of ground 

motion computations for this source zone 

num=98 for seismic source zones if averaging for earthquake 

location uncertainty is to be done at end of computations 

)> 

"'C 

"'C 

m 

2 

c 

>< 

n 
)oo 

"CC 

i 

:I 

~ 

>l!" 

n 

N 

t.n 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

for this source zone 

=99 for 1st fault set 

=0 for all fault sets after 1st 

=99 at end of computation 

yrnoc=number of years over which the earthquake occurrences 

take place 

iprint = -1 no statistics for this <intermediate> set of 

occurrences 

=1 statistical calculations and printout Con unit 16) 

=2 same as=1 plus summary file for plot, etc <unit 2) 

=3 omit printout; do summary file for plot <unit 2) 

totl ignored here but used for faults 

dumid = four character identifier for source zone or fault 

als=earthquake location uncertainty standard deviation in km1 

current value of als is ignored when num=O 

bls, sigls ignored here but used for faults 

2s. Jseg, ifr, i tot 

JSeg=num of pairs of quadrilatersl corner points in this source 

set (Jseg-1 quads). 

seismicity to be apportioned by fractional area among itot sources. 

ifr =set number (ifr =1.2 ... itot in sequence> 

itot=10 max; total quad pairs for itot sources maximum 50 

3s. Jseg quadrilateral endpoint lines (two points per line>: 

quad corner points Cleft long, left lath <right long, right lat), 

in decimal degrees (if dsw•O> degrees, minutes (if dsw=1) 

X 1' y 1 

1 

1 

1 

1 

x3,y3 

x2,y2 

x4,y4 

1 

1 1 

1 1 

1 1 

1 1 

x5, y5 x6,y6 

JSeg=3 in this example: quad endpoint pairs are 

< x 1, y 1) --- < x2, y2) 

Cx3.y3>---<x4,y4) 

<x5,y5>---<x6,y6> 

subregions of a set are defined as shown 

subregions of one set are separate from those of other sets 

repeat 2s, 3s for itot sources 

4s. noc<l> 1=1, lev <lev=12f6. 2> <lev=12 maximum> 

numbers of events expected in yrnoc years in each 

magnitude interval for earthquake occurrences 

5s. fm( 1), 1==1, lev < 12f6. 2> 

fm<l>=center of magnitude interval for which noc(l) occur 

repeat 1s. thru 5s. for remaining sources 

fault inputs: 

1f. num, yrnoc, iprint. totl. dumid, als. bls, sigls 

parameters defined as in 1s. above 

**note num=99 for first fault, =0 thereafter for successive faults 

als.bls,sigls rupture length parameters 

if als, bls non zero: 

length=10**<als +bls*m +fr*sigl> 

where m=magnitude; sigl=standard deviat on <in log length> 

fr=<normally) 5 values in range <-2,+2) if s gl non zero 

fr=O (1 value>--mean rupture length only if s gl=O 

if als.bls.sigls=O Cor blank> previous values are used 

if no previous values input, default values are used: 

als=-1.085 bls=.389 sigls=. 52 

totl=distance between faults if this set of faults is a set 

c of "dummy" faults used to approximate a uniform field c of faults. 

c totl=O if individual fault or a number of well defined faults are used. 

totl determines whether ground motions for a site are to 

be smoothed in distance or smoothed in magnitude 

<smoothed in magnitude if totl=O> 

<smoothed in distance <totl nonzero>. 

c 2f. Jseg, ifr, itot (free field) 

Jseg=num of fault segment end points to be connected 

into single fault (Jseg-1 segments> 

ifr, itot as defined in 2s 

3f. <xl<i, i(fr), yl<i. i(fr), i=t. Jseg) (free field> 

xl=long yl=lat <degrees) Cifr=fault number) 

Jseg=24 max itot=26 max 4f. noc < 1) same as in 4s. 

5f. fm<l> same as Ss. 

repeat 2f, 3f for itot faults 

repeat 1f. thru Sf. for remaining faults 

end with num=99 <omit other inputs on line> 

c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

character idsw*16, nam*10,tytle*BO,dumid*4 

data imaxsq/2500/ 

common/cc/dumfac,fac(2500> 

dimension idsw<2>, ibrC10), lsub<26), dtotC26) 

dimension xpt2C500>.ypt2C500), xpt1<500>.ypt1(500) 

dimension prlsC10>,rllsC10.25> 

dimension nrlsvC25),prlrlC10,25),rll<10.25) 

dimension disl<SO>,anglCSO> 

common/ata/ataC20,25> 

common/inout/mdum<4>.stsize 

common/slinds/sld1<4>.perpaC4>.sld2<4>.sllqC4>, 

& minaC4),perpaqC4),sldum2(8) 

common/extra/isid.shdis.distl<24), inear 

common /xlyl/xll(4), yl(4), xrC4), yr(4), xcC4), yc<4> 

common/itabs/dtabC105,25>.maxa<25> 

common/magdis/Jent.mdis 

dimension tmdif(10>.atabC20,25),tmC10),atadif(20,25) 

common/rtb/rtabC20> 

common/xyarea/xsav(2,50),ysav<2,50),sareaC50) 

t*******need maxdim=dimension of reg & of rawbin(maxdim, 1600) 

data maxdim.maxar/55, 1600/ 

dimension reggC1600> 

dimension regC55),rawbinC55, 1600), smobin<55, 1600) 

common/rdis/rdis 

common/wds/max.maxp.maxsp 

common/ext/prob,ex,exl 

common/tims/ftim(2Q), JtimC20>.ntims.nwt. iprint 
~ 

<" 

VI 

tD 

;;· 

.. 

;;· 

~ 

> 

~ 3 

'i 

;-

.. 

-.::1 

~ 

~ Q' 

.. 

~ ;;· 

3 

rr 

:I: 

E 

a. 

~ "" 

~ ~· 

data iscoml5001 

c*** iscom =dimension of following sin.cos arrays 

commonlsincoslsinx<500>,cosx<500),siny(500),cosy<500> 

commonl~uadlmaxq,ac(105),dumlog,aclog(105) 

dimension x<4>,y<4>,noc(24),fm(25) 

commonllinpsldist<24,2b>.disq<24,2b), Jseg, Jsegm, 

& xl<24,2b>.yp<24,2b>.coa<24.2b),sia<24.2b> 

commonllindisldls<24>.perp<24>,min<24),persq<24> 

dimension Jm(26), disrem<24> 

commonlaccellgm(25>,exdif(101),atah<20,25>,aclim<2> 

c arrays on previous line inserted for spreading out accelerations 

real noc 

integer dsw 

data !inset, line,nrlsiO,O,OI 

data idswl 'decimal degrees ', 'degrees, minutes' I 

data r.pi2.pi.pihl6378 .• 6.2831854,3.1415927,1.57079631 

rad=p i 1180. 

input on unit 15 

rewind 15 

saved on unit 3 for later use or continuation 

rewind 3 

allerr=O. 

dlevsv=O. 

flevst=O. 

levno=O. 

tot5=0. 

read(15,20> tytle 

20 format<a80) 

write<lb.30> tytle 

30 format<'1'a80) 

isw tells whether to include data from a previous run, zero if no. 

sigmax•maximum value of sigma <km. in earthquake location> 

read<15,*>isw.sigmax 

if<isw. e~.O> write<16.40) 

40 format<' isw=O: new run--no previous results included') 

if(isw. e~. 1> write <16, 50) 

50 format<' isw=1: run continuation; add to previous results') 

if<isw.e~.O> go to 70 

idrc=l 

bO read(3,end=70> <smobin(iqrp. idrc), iqrp=l.maxdim> 

idrc=idrc+l 

go to bO 

70 rewind 3 

prob=extreme probability in decimal 

ntims are number of 'Jtims' fpr which calculation is done 

Jtims are the durations for which the extreme motions are 

to be calculated at the prob extreme probability level 

read < 15, *> prob, ntims. < Jtim( JV), JV=l• ntims > 

if <prob. lt. 1. > go to 90 

write<lb,80> prob 

80 format<' prob ='f8.3' must be decimal less than one'> 

call exit 

90 continue 

write <lb. 100) prob, < Jtim< JV), JV=L ntims > 

100 format<' extreme probability 'f6.3 

&I' for exposure times <years> '2<10i51)) 

ex=-alog<prob> 

exl=alog<ex) 

do 110 JV=l,ntims 

110 ftim<Jv>=Jtim(JV) 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

read< 15, *> scale, dsw, sd, inos 

write<16, 120> scale, idsw<dsw+1), idsw<dsw+1), sd 

120 format<' scale factor for ground motion "box" levels=' 

&f7.21' coordinates input in 'a16, 

&I' coordinates are printed in 'alb, 

&I' variability in attenuation, sigma='fb.2> 

scale=multiplicative factor for acceleration levels 

dsw=O if input is decimal degrees; =1 for degrees,minutes 

sd=standard deviation for log acceleration 

read <x1, yl), <x2, y2) =long, lat in dec. degrees of 2 points to transform to equator 

<xl,yl) to <O.,O><x2,y2> to (dist. between points,O> 

130read<15,*) xl.y1.x2.y2 

write <16, 140) xl. yl, x2, y2 

140 format<' grid oriented parallel to great circle thru (' 

& f7. 2', 'f6. 2'>, < 'f7. 2', 'fb. 2'> '> 

convert to decimal degrees unless dsw=O. 

if<dsw. ne. O> call cond<xl. x2, yl. y2> 

x1=xl*rad 

x2=x2*rad 

y1=y1*rad 

y2=y2*rad 

transform to equator 

call inieqr(x2, y2, x!, yl> 

call toeqr<xl, y1, x1t, y1t> 

call toeqr<x2.y2,x2t.y2t> 

read in upper left and lower right corners of affected rectangle 

--degrees, minutes ie 20.30=20 degrees.30 minutes <dsw=1> 

<decimal degrees 20. 50 = 20 112 degrees (dsw=O> ) 

opposite corners of gridded area <in latitude and longitude) 

gridded area becomes rectangle surrounding new equator. 

sites are at uniform increments in latitude and longitude 

within this gridded area in new coordinate system. 

in order to limit affected area you can specify 

beginning and ending rows and columns as irowl and 2 

and icoll and 2 <next input line> 

phl=phi=latitude 

fl1=lambda=longitude 

read(15,*> fl1.ph1,fl2.ph2,phinc 

write < lb, 160) fll, ph 1, f12, ph2 

lbO format(' corners of gridded area-upper left='f7. 2', 'f7. 2, 

&123x' lower right='f7. 2', 'f7. 2> 

convert to decimal degrees <unless dsw=O> 

if(dsw. e~.O) go to 170 

call cond<ph1, fll, ph2, fl2> 

call condac<phinc> 

170 continue 

180 

determine number of subregions in map area and set up subscripting irow=no of rows 

icol=no of cols 

write<t6, 180> phinc.phinc 

format(' longitude increment='f7.4' <decimal degrees)' 

&I' latitude increment ='f7.4' (decimal degrees>'> 

phinc=phinc*rad 

flinc=phinc 
t 

>e· 

~ 

N 

....... 

dphinc=phinc*r 

dphinch•dphinc/2. 

dflinc=flinc*r 

fl1=f11*rad 

ph1-ph1*rad 

fl2=fl2*rad 

ph2=ph2*rad 

transform to e~uator and x-y rectangular coordinate system 

convert to kilometers 

call toe~r(fl1.ph1,flam1.phi1) 

call to e~r(fl2.ph2,flam2.phi2> 

if(phil. gt. O> go to 210 

if(phi2. gt.O> go to 220 

apparent error in inputs or xformataion 

190 write<16.200> 

200 format(' probable input error--felt region upper left,' 

&' lower right'/' must be on opposite sides of new e~uator') 

call exit 

210 if< ph i2. le. O> go to 230 

go to 190 

220 phs=phil 

phil=phi2 

phi2=phs 

fls=flaml 

flaml=flam2 

flam2=fls 

have proper orientation--continue 

230 ph1=phi1 

if(flam1. lt. O> go to 240 

fl1=flam1 

find total number of cols 

go to 260 

240 write<16,250)flam1.flam2 

250 format<' apparent error in xforming lat and long limits' 

& 'to e~r--'1' flam1•'f9. 3' fl2='f9.3> 

call exit 

260 continue 

ph1r=ph1*r 

phi2r=phi2*r 

fllr=fll*r 

flam2r=flam2*r 

c adJust to include area beyond edges for later smoothing 

determine nbord=number of additional rows and columns at 

c which ground motions must be calculated to account for 

c earth~uake location variability. 

nbord=2.*sigmax/dphinc+.01 

if(sigmax. ·~· 0. > go to 280 

if(nbord.ge.2> go to 280 

write (16,270) dphinc.sigmax.nbord 

270 format<' increments between sites ='e12. 5' too large for' 

& 'sigmax='e12. 5' nbord='i5) 

call exit 

280 ph1=ph1+nbord*phinc 

if <flam2.gt.fl1> flinc=-flinc 

fl1=fll+nbord*flinc 

irow=<ph1-phi2)/phinc+1.01+nbord 

icol=<fl1-flam2)/flinc+1.01+nbord 

write<16,290)irow. icol.nbord 

290 format<' gridded region contains 'i4.' rows, 'i4' cols' 

&' including border 'i4' rows and cols') 

read<15,*)irow1. irow2. icoll, icol2 

limits rows and columns of felt points for which calculation is c to be made. 

write<16,300)irow1. irow2, icoll, icol2 

300 format<' for this run begin at row 'i3' end row 'i3 & ', begin col 'i3' end col 'i3) 

for earth~uake location variability, compute at more sites 

c calculate at additional rows and cols 

irow2=irow2+2*nbord 

icol2=icol2+2*nbord 

icols=icol2-icol1+1 

irows=irow2-irow1+1 

if<irow1. e~.O> irows=O 

if< icoll. e~. O> icols=O 

if<irowl. le. irow2> go to 320 

write<16,310) irow1 

310 format(' irow1 ='i4' too large--error stop') 

call exit 

320 if ( i c 0 11. 1e. i c 0 12 ) g 0 t 0 340 

write(16,330> icol1 

330 format(' icol1 ='i4' too large--error stop'> 

call exit 

340 if< icols*irows. le. maxar) go to 370 

350 write<16,360)irows. icols, indvpt.maxdim 

360 format(' irows='i4' icols='i4' pts on line='i4 

&' exceeds maxdim of 'i4) 

call exit 

370 continue 

c read in number of line segments containing sites at which 

c calculation is to be done (in addition to or instead of 

c computing for sites on a fixed grid) 

c zero if no line segments 

uply=phl*r 

uplx=f11*r 

boty=phi2*r 

botx=fllam2*r 

write<16,380) uplx.botx.uply.boty 

380 format(' new coordinates <km) gridded area' 

&I' upper left= '2f9. 2 '; lower right= '2f9. 2> 

read(15,*)indv 

write(16,390) indv 

390 format<' sites are also located on 'i3' line(s)') 

if <indv. e~.O) go to 460 

calculate accelerations at nvs sites on each of indv lines 

read in end point pairs <long. lat> xe1.ye1, xe2.ye2 

for each of indv lines 

nvs sites are evenly spaced on line* 

c nvs=2 gives accelerations at end points of line; 

nvs=3 acceleration at 2 end points plus center of line.etc. 

read<15,*> nvs 

write<16,400) nvs 

400 format<i4' sites per line') 

indvpt•O 

do 450 i=1, indv 

read ( 15, *> xe1, ye1. xe2, ye2 

write<16.410> i, xel. ye1, xe2, ye2 

410 format<' line 'i2' end points at 'f7. 3', 'f7. 3' and 'f7. 3 

&', 'f7. 3) 
N 

= 

VI 

t'D 

iii' 

.. 

~ 

> 

g 

3 

"CC 

c 

;' 

.. 

~ 

~ ~ 

3 

~ 

~ iii' 

3 

;:;· 

:I 

~ Dol 

a. 

m 

=. 

~ cs· 

= 

if(dsw. ne. O> call cond<xeL xe2, ye1, ye2> 

xe1=xe1*rad 

xe2=xe2*rad 

ye1=ye1*rad 

ye2=ye2*rad 

call toeq_r( xeL ye1, xo1, yol > 

call toeq_r<xe2.ye2, xo2,yo2) 

fr=O. 

if <nvs. gt. 1 > frinc=l. I (nvs-1) 

do 440 iq_=1,nvs 

indvpt=indvpt+1 

if( indvpt. le. 500) go to 430 

write ( 16, 420> 

420 format(' more than 500 single points req_uested as sites' 

8c 1 X I , d 0 500 , ) 

go to 460 

430 xl1=xo1+fr*<xo2-xo1> 

yp1=yo1+fr*<yo2-yo1> 

cosy<1>=cos<yp1> 

siny(1)=sin<yp1) 

sinx<1>=sin<xl1> 

cosx<1>=cos<xl1> 

call backw(1, L xptl<indvpt>, yptl<indvpt)) 

xpt2<indvpt>=xl1*r 

ypt2<indvpt>=yp1*r 

fr=fr+frinc 

440 continue 

450 continue 

460 if(icols*irows+indvpt. gt.maxar) go to 350 

nrc=icols*irows 

Jent•number of magnitude entries in table of acceleration as a 

function of magnitude and distance 

mdis=number of distance entries 

read(1~•*> Jent.mdis 

if < Jent. le. 10) go to 480 

write (16.470) Jent 

470 format(' Jent='i3' too large--max 10 allowed') 

call exit 

480 if <mdis. le.20) go to 500 

write<16,490) mdis 

490 format<' mdis='i4' too large--max 20 distances allowed') 

call exit 

500 continue 

read identifier and magnitudes for which acceleration versus 

distance table follows. magnitudes must be in decreasing order 

read(1~•*>nam, <tm(i), i=L Jent> 

print magnitude headings for attenuation function. 

write <16.~10> nam 

510 format<' attenuation function 'a10/30x'magnitude') 

write<16,520) <tm(i), i=1, Jent> 

520 format<' dist<km> 'f7.2.10f10.2> 

read in table of acceleration vs distance values for set 

of magnitudes <distance increasing in table> 

do ~80 m=1,mdis 

read<1~•*> rtab(m), <atab(m, J), J=1, Jent> 

i f < at a b < 1, 1 >. 1 t. 1 0 > g o to 540 

write <16, 530) rtab(m), <atab<m, J), J==L Jent) 

530 format<1h f10.2, 10f10.2> 

go to 560 

540 write<16, 550) rtab(m), <atab(m, J), J=L Jent) 

550 format<lh f10.2,8f10. 5> 

take logs of tabular values since interp. is in log ace. etc. 

560 do 570 J•1, Jent 

570 atab(m, J)=alog(atab(m, J)) 

580 continue 

define limits for acceleration boxes: accelerations in the 

range ac(J-1) < ac < ac(J) will be collected the Jth box. 

max=number of boxes; acceleration values too large for 

any of the box levels correspond to maxp=max+1 

as here constituted boxes have limits ranging from 

.02 to 1.0 

c For motions exceeding 1<g>., provide 

c appropriate distance vs ground table and select 'scale' 

for multiplying boxes to get consistent range 

or rewrite and recompile subroutine BOX 

600 call box(scale.sd> 

basic step size for distances in source areas 

do 610 i=1,maxsp 

610 aclog(i)=alog(ac(i)) 

number of digits after decimal for print of 

accelerations (subroutine out> determined by scale 

dentify rows and cols for file02 

620 r2=irow2-2*nbord 

cl2=icol2-2*nbord 

write<2> irow1, ir2, icolL icl2, indvpt,ntims, (Jtim<Jv), JV=1,ntims), 

8cprob.sd 

c create table of sines and cosines for later back transformation 

c of felt points 

if<icol1. eq_.O) go to 680 

yp 1=<ph 1-< irow1-1. >*ph inc) 

xl1=<fl1-<icol1-1. )*flinc> 

irow2p•irow2+1 

icol2p=icol2+1 

if< irow2p. le. iscom> go to 650 

630 write(16,640> irow2p, icol2p. iscom 

640 format < i~' rows, ' i5' cols needed in sincos arrays, ' 

8c' dimension='i4/' need to increase dimension') 

call exit 

650 if(icol2p. gt. iscom> go to 630 

do 660 ira=irow1, irow2p 

cosy<ira>=cos(yp1> 

siny<ira>•sin<yp1> 

660 yp1=yp1-phinc 

do 670 i i=icoll. icol2p 

cosx<ii>=cos<xl1) 

sinx<ii)=sin<xll) 

670 xl1=xl1-flinc 

680 read(15,700>num,yrnoc,iprint,totl,dumid•als,bls,sigls 

numsav=num 

690 format ( 1 0yrnoc= 1 f6. o, 1 iprint= 1 i2, 1 for area 'a4) 

700 format<i2,f10.0, i2 ,f10.2,a4.2f6.2,f5.2) 

if<num.ne.99) go to 710 

if<line.eq_.l> go to 2950 

yrnoc=O here if only source areas. no faults included 
> 

l 

><. 

t"') 

N 

I.C 

c 

c 

c 

c 

c 

c 

if<11rnoc. eq. 0. > go to 2950 

line=l 

710 write<16,690)\lrnoc. iprint.dumid 

720 if< 1 ine. eq. 1> go to 830 

if<num. eq.98) write <16,705> als 

705 format(' smooth for earthquake location uncertaint\1 sigma=' 1 f5. 1) 

nbr=O 

ist=O 

iend=O 

read identifiers for source area inputs 

730 read<15,*>Jseg, ifr, itot 

write<16,740> ifr, itot 

740 format(' source 'i2' of 'i2) 

JSeg•number of pairs of quadrilateral end points in this 

single source 

itot=number of sources in this source area ifr=identifies which source. ifr=1,2 •.... itot 

if ( i t 0 t. le. 1 0 ) g 0 t 0 7 60 

write(16, 750> itot. JSeg, ifr 

750 format<' itot ='i4' to large--max 10 JSeg='i4 &c ' ifr='i4) 

call exit 

760 ist=iend+l 

iend=ist+Jseg-1 

if < iend. le. 50> go to 780 

write<16,770) 

770 format<' too man11 quadrilaterals in source--max 50') 

call exit 

read in boundaries for quadrilaterals for all subareas 

in this source 

read limits of seismic area x<i>=lambda(i), \l<i>=phi(i) 

upper left(x,\1> then upper right<x,\1> --1st quadrilateral 

2nd line lower left, lower right 1st quadrilateral= 

2nd card upper left,upper right 2nd quadrilateral, etc. 

780 do 810 ii•ist. iend 

read<15,*> <xsav<i.ii),\lsav(i,ii>• i=1•2> 

790 format(4f6.0) 

BOO format<lh 4f10.3) 

write<16,BOO><xsav(i, ii>• \ISav(i, ii), i=1·2> 

810 continue 

if<ifr.eq_.itot> go to 820 

nbr=nbr+1 

ibr<nbr>=iend 

go to 730 

820 num•iend 

go to 960 

fault line read 

830 J•1 

write<16,840> totl 

840 format<' Distance between dumm\1 faults='f5. 1> 

dsum=O. 

fault inputs 

850 read(15,*> Jseg, ifr,itot 

write<16,860) ifr, itot 

860 format<' fault 'i2' of 'i2) 

Jseg=number of end points of connected segments 

for the current fault 

JSeg=one plus number of segments 

c 

itot=number of distinct faults in this zone 

c ifr identifies current fault, ifr=L 2, ... itot 

if (Jseg. le. 24> go to 890 

870 write <16.880> Jseg, itot 

880 format<' JSeg='i4' (max 24> itot='i4' <max 26>-stop') 890 if<itot. gt.26) go to 870 

read in end points long, lat of each segment 

read(15,*><xl(i, J), \IP<i. J), 1=1. JSeg) 

wr i t e < 16, 900 > < x 1 ( i , J > , \1 p ( i , J ) , i = 1 , J s e g > 900 format<1h 4(f10. 2, fB. 2, '• ')) 

Jm(J)=Jseg 

do 920 i=1, Jseg 

convert to decimal degrees unless dsw=O 

if<dsw. eq.O> go to 910 

call condec<\lp(i, J» 

call condec(xl<i, J» 

910 continue 

transform to equator--rectangular coordinates 

\lin=\lp<i, J)*rad 

xin=xl<i, J)*rad 

call to eqr < x in. \1 in. xout, \lOUt> 

xl(i. J)=xout*r 

\IP(i, J)""\IOUt*r 

920 continue 

compute 1 ine parameters, lengths for fault segments Jsegm=Jseg-1 

dtot<J>=O. 

do 930 i=l. Jsegm 

xdelta=xl<i+l. J)-xl(i, J) 

\ldelta•\lp(i+l. J)-\lp(i. J) 

disq(i, J)=xdelta*xdelta+\ldelta*\ldelta 

dist<i, J)=sqrt(disq(i, J» 

c oa ( i, J) = ( x 1 < i + 1. J >- x 1 ( i, J > ) I dis t ( i, J > 

s i a< i. J > = ( \1 p ( i + 1. J >-\I p < i. J > >I dis t ( i, J) 930 dtot(J) 111 dtot<J>+dist(i, J) 

dsum=dsum+dtot(J) 

if< ifr. eq. itot > go to 950 

940 J""J+l 

go to 850 

read in num of occurrences for each magnitude <12f6.2) both area and fault input 950 Jtot=J 

960 read < 15, 1170> <noc ( i), i=l, 12> 

lev=12 

970 if(noc(lav>.ne.O.O> go to 990 

lev=lev-1 

if<lev. gt.O) go to 970 

write(16.9BO> 

980 format<' apparent input err for no of occurrences at each level') go to 2990 

990 continue 

write(16, lOOO>lev 

1000 format<' nr of levels of seismicit\1 = 'i2> 

read in corresponding mag. interval center points. 

1010 format(' before normalizing to ratel11ear'> 

1020 format( ' earthquake rate I 11ear'> 

<lev of them) 

read(15, 1170)(fm(i), i=1, lev ) 

c magnitudes must be evenl\1 spaced and in increasing order. 
Wol 

= 

VI 

tD 

;;· 

.. 

~ 

> 

~ 3 

~ 

if 

.. 

"'1:1 

ci 

iil 

3 

;-

.. 

~ ;;· 

3 

;:;-

:::1: 

$ 

a. 

~ 3" 

9t. 

~· 

if not entered in increasing order, program reverses order 

of magnitudes and corresponding earthquake occurrences 

at each magnitude. 

Don't confuse this with the original table of ground 

motions as a function of magnitude and distance. 

The magnitudes in that table must b~ read in decreasing 

<rather than increasing) order. 

ifCfmCU. lt. fmC2> > go to 1040 

if( lev. eq. 1) go to 1040 

levh=lev/2 

ll=lev 

do 1030 lz1,levh 

sav=fm(l) 

fmC 1 >=fmC 11 > 

fmCll>=sav 

sav=noc(l) 

noc<l>•nocCll> 

noc(ll)zsav 

1030 11=11-1 

1040 continue 

if( lev. gt. 1 > go to 10:50 

we cannot calculate beta from numbers of earthquakes 

if only one level--so give a default beta so we can later spread out 

beta=-2.0 

fm<2>=fm<1>+.6 

delm•.6 

go to 1060 

1050 continue 

dlev=absCfmC1>-fm<2>> 

beta=alogC nocC2)/nocC1) )/dlev 

1060 writeC16, 1070)dumid, beta 

1070 format< 1 x, a4, ' beta= 'f8. 4 > 

ifCyrnoc. eq. 1. ) go to 1090 

writeC16,1010> 

writeC16, 1100> CnocCll), 11=1, lev) 

this is a normalization to rate per year 

do 1080 lm=1, lev 

noc<lm>•nocClm>lyrnoc 

1080 continue 

1090 writeC16, 1020) 

writeC16, 1100)(noc(l), 1=1, lev> 

1100 format<' occurrences='12f10.6> 

writeC16, 1110)(fm(l), 1=1, lev> 

1110 format<' magnitudes='12f10.2) 

summoc=O. 

do 1120 lm•1, lev 

1120 summoc=summoc+nocClm> 

evaluate at twice as many magnitudes as initially read 

in if inos•1 

ifCinos. ne. 1) go to 1140 

dlev=<fmC2>-fm(1))/2. 

zac•expCabsCbeta>*dlev> 

noc<l>=noc<1>*zac/Cl. +zac> 

summoc2=noc<l> 

lev=2*lev 

fm(l)=fm(1)-dlev/2. 

do 1130 11=2, lev 

fmCll>•fm<ll-1>+dlev 

nocCll>•nocCll-1)/zac 

1130 summoc2=summoc2+noc(ll> 

writeC16, 1150> summoc,summoc2 

writeC16, 1100> Cnoc<l>,l=1,lev> 

writeC16, 1110> CfmCl), 1=1. lev> 

1140 continue 

1150 format(' total eqs='2e12. 5> 

1160 continue 

1170 format<12f6.2> 

if Cline.eq.O> 

determine which 

ifCals.eq_.O) 

al=als 

bl=bls 

sig l=sig h 

go to 1250 

rupture length magnitude relationship to use go to 1240 

c compute break length 

1180 wr i t e C 16, 1190 > a 1, b 1, s i g 1 

1190 format<' fault rupture length parameters al='f7.3' bl='f7.3 

c 

1200 

1210 

1220 

ll<' sigl='f6. 2) 

linset•1 

if Csigl. gt.O> go to 1230 

which single break default do we want 

if bl•O, al non zero use old Algermissen-Perkins default values 

use parameter values supplied for one break only 

do 1210 i=L lev 

rllsC1, i)=10. **Cal+bl*fm(i)) 

prlsC1>=1. 

nrls=1 

go to 1350 

set up rupture lengths and probabilities as fct of magnitude 

nrls lengths per magnitude 

c al, bl, sigl non zero 

1230 call pbreak<fm,al,bl,sigl,rlls,prls, lev.nrls> go to 1350 

1240 if Cnrls. eq_. 1 > go to 1200 

if <1 i n set . e q_. 1> g o to 1230 

compute default values 

al=-1. 085 

bl=.389 

sigl=. 52 

linset=1 

go to 1180 

convert quad corner points to decimal degrees Cif dsw =1> c 

1250 do 1280 ii=l, num 

do 1280 i=1,2 

1260 

c 

1270 

1280 

c 

if<dsw. eq_.O> go to 1260 

call condec < xsav< i, i i)) 

call condec<ysav(i, ii)) 

continue 

x in=xsav< i, i i >*rad 

yin=ysavCi, ii>*rad 

call to eq_rCxin,yin, xout,yout> 

xsav<i, ii)=xout*r 

ysav(i, ii)=yout*r 

writeC16, 1270) xsav<i, ii), ysav(i, ii), ii 

format(' xsav,ysav='2f9.3' ii='i3) 

continue 

1290 nm=num-1 
f 

a. = 

jji" 

n 

~ 

... 

stot=O. 

do 1330 ii=l.nm 

if<nbr.e~.O> go to 1310 

do 1300 i~=1,nbr 

if(ii. ·~· ibr<i~)) go to 1330 

1300 continue 

c set line parameters for subsource 

1310 call setreg(ii) 

do 1320 i-1.2 

x<i>=xsav(i, ii) 

x ( i+2>=xsav( i, i i+1 > 

y < i > =y sav < i, i i > 

1320 y(i+2>=ysav(i, ii+1> 

c find area of subregion ii 

s1=abs<x<1>*<y<2>-y<3>>+x<2>*<y<3>-y<1>>+x<3>*<y<1>-y<2>>> 

s2=abs<x<4>*<y<2>-y<3>>+x<2>*<y<3>-y<4>>+x<3>*<y<4>-y<2>>> 

sarea(ii)= (s1+s2>12. 

stot=stot+sarea<ii> 

1330 continue 

eua=noc<1> 

c find earth~uake rate per unit area for lowest magnitude 

rate will be printed but not used in program 

delm=<fm<2>-fm<1>>12. 

fm1=fm(1)-delm 

fm2=fm(1)+delm 

eua=eua/stot 

write<16, 1340) dumid.stot.eua,fm1,fm2 

1340 format(1h a4' area='f9.0' sq km, rate/sq km='e12. 5, 

& 'for mags 'f5.2'-'f5.2) 

1350 dlev=fm(1)-fm<2> 

create table of log accelerations for set of magnitudes 

c in <tm> for the original set of mdis distances as a function 

c of distance 

next three lines test whether table has already 

been computed--if so skip over 

if< lev. gt. levno > go to 1360 

if(fm(1).ne. flevst> go to 1360 

if<dlev.e~. dlevsv> go to 1520 

1360 do 1370 il1=2· Jent 

1370 tmdif(ill>=tm<ill>-tm<ill-1> 

do 1420 11=1, lev 

do 1380 il•2• Jent 

if(fm<l1>. gt. tm(il)) go to 1390 

1380 continue 

il=Jent 

1390 i lm=il-1 

do 1400 J•1, mdis 

1400 ata<J• ll>=atab<J• ilm>+<fm<ll>-tm(ilm))/tmdif(il) 

& *<atab < J• i 1>-atab < J• ilm)) 

do 1410 J=2,mdis 

1410 atadif<J• ll>=ata<J• 11>-ata<J-1• 11> 

1420 lsub<ll>=il 

find distance corresponding to each acceleration in 

table ac <set up in subroutine box> for each magnitude 

do 1310 1=1, lev 

do 1500 J•l.maxp 

do 1430 i=2,mdis 

if(aclog<JL gt. ata(i, 1)) go to 1450 

1430 continue 

c beyond table 

1440 d tab< J• 1 >=rtab <md is> 

go to 1500 

1450 if<atadif(i, l>.ne.O> go to 1490 

c no change in acceleration at this distance 

1460 if(i. e~.2> go to 1480 

write(16, 1470) i, 1. J 

1470 format<' apparent error indist-table i. 1, J='3i4) 

call exit 

1480 dtab < J• 1 >=0. 

go to 1500 

1490 dtab(J, l>=rtab<i>-<rtab<i>-rtab<i-1»*<ata(i, 1>-aclog<J» 

&/atadif(i, 1> 

if(dtab<J• 1>. lt. O> dtab<J• 1>=0. 

1500 continue 

1510 continue 

dlevsv111 dlev 

levno=lev 

flevst=fm< 1) 

1520 continue 

c add half values of mag for later spreading out 

dmag=fm<2>-fm<1> 

bmag=fm<1>-dmag/2. 

levp=lev+1 

do 1530 i=1,levp 

gm<i>=bmag 

1530 bmag=bmag+dmag 

do 1570 11=1,1evp 

do 1540 i1=2, Jent 

if<gm<ll>. gt. tm<il)) go to 1550 

1540 continue 

il=Jent 

1550 do 1560 J=1.mdis 

atah(J, ll)=atab<J• il-1)+(gm<ll>-tm<il-1))/tmdif<il> 

& *<atab(J, il>-atab<J• il-1» 

1560 continue 

1570 continue 

take exp differences for later accel spread out 

the "exdif" values are used in subroutine "spread" 

if <lev. e~. 1) betad=2. 0 

if(lev .ne.l> betad=alog<noc(1>/noc(2)) 

beta ... betad/dmag 

fmu1=1. /( 1. -exp <-beta*dmag >) 

frak=exp <-beta*dmag/100. > 

frm•frak 

exd if< 1>=0. 

do 1580 i•2.101 

exdif<i>=<1.-frm>*fmul 

1580 frm=frm*frak 

dtmax•dtab < 1. 1 > 

if<dlev. lt. O> dtmax=dtab(1, lev> 

find entry for each magnitude in distance-acceleration table 

for zero distance <highest acceleration for that magnitude> 

do 1600 1•1, lev 

J=maxp 

1590 if(dtab<J• 1>. ne. 0. > go to 1600 

J""J-1 

go to 1590 

1600 maxa<l>•J 

irow1s=irow1 
w 

N 

VI 

II> 

;;;· 

:::!. 

~ 

>-

~ 3 

"CC 

c 

If 

.. 

., 

~ Dl 

3 

0' 

.. 

f 

~· 

;:;· 

::1: 

E 

a. 

!r 

f 

~· 

irow2s=irow2 

icol1s=icol1 

icol2s=icol2 

idrc=O 

do 2770 ispt=1,2 

if< ispt. eq_. 1 > go to 1610 

case of calculation done for single points on line 

if<indvpt. eq_.O) go to 2770 

irow1s=1 

irow2s=1 

ico11s=1 

icol2s=indvpt 

go to 1620 

1610 if<irow1. eq_.O) go to 2770 

1620 do 2760 ira=irow1s,irow2s 

ymid=Cph1-Cira-1. >*phinc>*r 

do 2750 ii=icol1s, icol2s 

idrc=idrc+1 

if<ispt.eq.2> go to 1630 

fi i=i i 

xmi d=Cfl1-< fi i-1. >*flinc >*r 

go to 1640 

1630 xmid=xpt2Cii) 

ymid•ypt2< ii > 

1640 if(isw.eq_.O) go to 1660 

c after first time thru, read accelerations calculated 

c for previous sources at this site <saved in 2 dimensional 

array rawbin) into reg (1 dimensional> array. 

do 1650 iq_rp=1,maxdim 

1650 regCiqrp>•rawbin<iqrp, idrc> 

go to 1680 

c first time thru clear reg array if not continuation run 

c 

c 

1660 do 1670 11•1,maxp 

1670 regCll>=O. 

reggCidrc>=9999.999 

set regg to min dist from fault 

for fixed affected area in row ira, col ii 

c source area computation 

1680 if< 1 in e. eq_. 1> go to 2120 

do 2110 k•1,nm 

if(nbr. eq_.O) go to 1700 

do 1690 iq_=1,nbr 

if( k. eq_. ibr< iq_)) go to 2110 

1690 continue 

1700 approx=O. 

do 1710 1•1,2 

xll<i>=xsavCi, k> 

x 11 < i+2>•xsavC i, k+1) 

yl<i>•ysavCi, k) 

y 1 < i+2>•ysavC 1, k+1 > 

1710 continue 

xr<3>•xll<1> 

yr<3>=y1Cl) 

xr<1>=xll(2) 

yrCl>=ylC2) 

xr<2>=xll<4> 

yr<2>=yl<4> 

xr<4>=xll<3> 

yrC4>=ylC3) 

aread=sarea(k)/stot 

determine whether site is inside or outside k th quadrilateral 

subregion. find distance rc to closest endpoint 

and rf to farthest end point. Set up table of distance versus 

arc length for distances in the range rc to rf using 

subroutines outsid or inside. 

itst=no of radii at fixed increments between rc and rf 

at which arc lengths are to be evaluated in table 

itst=7 set in subroutine rrisk--recompile to change 

call rrisk < xmid. ymid, in, k, dis, rc, rf, itst> 

if< in. eq. 1> go to 1740 

if(rc. lt. dtmax> go to 1730 

do 1720 11=1, lev 

1720 reg<1>=reg<1>+aread*noc(ll> 

go to 2110 

1730 rcs=rc 

dislC1>=rc 

anglC1>=0. 

itst=itst+1 

ibeg=2 

go to 1750 

1740 disl<l>=O. 

ang 1< 1>=0 

d isl <2>=rc 

anglC2)=pi2*rc 

ibeg=3 

if<rc. eq.O> ibeg=2 

rcs=O. 

itst=itst+2 

compute distance versus arc length table for even increments 

c in distance 

1750 do 1810 iii=ibeg, itst 

if< in. eq. 0. > go to 1760 

call inside<xmid, ymid, k, dis, angle) 

go to 1780 

1760 call outsid ( xmid, ymid, k, dis, angle, kerr) 

1770 formatC'err outside for ira='i4' ii='i4> 

1780 continue 

dislCiii>=dis 

ang 1 <iii >=angle*dis 

ifCkerr.eq.O> go to 1810 

write <16, 1770) ira, ii 

do 1790 imkl•1, iii 

1790 writeC16, 1800> dislCimkl>,anglCimkl> 

1800 format<' dist, ang='2e12. 5> 

1810 dis=dis+stsize 

itst=itst+1 

dislCitst>•rf 

anglCitst>=O. 

compute distance versus arc length table for distances 

to selected vertices and perpendiculars to edges 

rcsq_=rc*rc 

rfsq_=rf*r# 
i 

~· 

~ 

~ 

~ 

do 1940 i=1,4 

ire=O 

if(in. e~.O> ire=2 

i f < s 11 ~ < i >. le. r c s ~ + 1. > g o to 1 930 

if < s 11 q < i > . g e. r fs q -1. > g o to 1 930 

dis=sqrt<sllq<i>>-.001 

i f < d i s. le. r c ) g o to 1 930 

add angle, distance corresponding to vertex 

if <in. eq.O> go to 1830 

call inside<xmid, vmid, k, dis. angle> 

go to 1860 

1820 if(mina<i>.ne.3) go to 1940 

c closest side next 

dis•abs(perpa<i>> 

if<dis. le. rc+. 5) go to 1940 

if<dis. ge.rf-.5> go to 1940 

dis•dis-.001 

call inside< xmid, ymid. k. dis. angle) 

ire=1 

go to 1860 

1830 continue 

call outsid < xmid. ymid. k, dis. angle. kerr> 

go to1860 

1840 if(mina(i).ne.3> go to 1940 

dis=abs(perpa(i)) 

ire=1 

do 1850 iql=1,4 

if(abs<s11q(iql>-perpaq(i)). le. 2. > go to 1940 

1850 continue 

if(dis. le. rc+1.) go to 1940 

if( dis. ge. rf-1. > go to 1940 

dis=dis-.001 

call outsid<xmid, vmid. k, dis, angle, kerr> 

1860 continue 

if<kerr.eq.O> go to 1880 

write <16. 1770> ira. ii 

do 1870 imkl=1,iii 

1870 write(16, 1800) disl<imkl>.angl<imkl> 

1880 continue 

do 1890 11=1, itst 

if(dis-disl<ll>> 1900,1920,1890 

1890 continue 

c insert new distance. angle 

1900 Jk=itst 

do 1910 mm=ll, itst 

angl(Jk+1>•angl(Jk) 

disl(Jk+1)•disl(Jk) 

1910 Jk=Jk-1 

disl<ll)==dis 

angl(ll>=angle*dis 

itst•itst+1 

1920 if(ire.eq.O> go to 1820 

if<ire.e~.2> go to 1840 

go to1940 

1930 if<in>1840, 1840,1820 

1940 continue 

itst=itst+1 

disl<itst>•rf 

ang 1< i tst >=0. 

c for each magnitude, source. determine fractional source area giving 

acceleration a: ac(i-1) <a < ac(i) for each entry in acceleration 

table ac. equivalent to finding fractional source area at distance 

from site between successive entries in acceleration-distance 

table: dtab(i, 1> < d < dtab<i-1.1> for magnitude 1. 

interpolates for distance is disl vs angl table 

area=integral<angle*dist*de1ta(dist)) from dtab(i-1.1> to dtab(i.l) 

do 2090 11~~:1, lev 

stnoc•noc<ll)/stot 

iswd=O 

maxm=maxa(ll) 

if<maxm. le. 0.) go to 1960 

approx•O. 

if<in. eq.O> go to 1950 

J1o=maxm 

go to 1980 

1950 if<rcs. lt. dtab(1, 11» go to 1970 

1960 reg<1>=reg(l)+aread*noc(ll> 

go to 2080 

1970 J1o=indpt(dtab<1, 11>.rcs.maxm) 

1980 J•J10 

iqs=2 

do 2060 Jk=l,Jlo 

i f ( J k . g t. 1 ) g 0 t 0 1 990 

dbot=rcs 

ddis=dtab<J• 11> 

i f < d tab < J .11> . 1 t. r f > g o t o 20 1 0 

reg(J+1>=reg(J+1)+sarea(k)*stnoc 

go to 2080 

1990 if(dtab(J, 11>. 1t. rf> go to 2000 

dbot=dtabCJ+1, 11> 

ddis=rf 

iswd=1 

go to 2010 

2000 ddis=dtab<J• 11> 

dbot=dtab < J+1, 11 > 

2010 anarea=O. 

do 2020 iq=iqs.itst 

i f < d b o t. 1 e. d i s 1 < i q > > g o to 2030 

2020 continue 

iq=itst 

c distance exceeds last value-error since table set up to handle this 

2030 dtop=disl<iq> 

if(disl<iq>. gt. ddis) dtop=ddis 

dmid=<dtop+dbot)/2. 

delr=dtop-dbot 

fr=<dmid-disl<iq-1))/(disl<iq)-disl<iq-1)) 

ang•angl<iq-1)+(angl<iq>-angl<iq-1>>*fr 

anarea=anarea+ang*delr 

if<dtop. eq. ddis> go to 2050 

dbot=dtop 

iq=iq+1 

if(iq. le. itst> go to 2030 

wr i t e < 16, 2040 > 

2040 format<' err in arc area computation') 

go to 2030 

c interpolate in angle , distance table 

2050 approx=approx+anarea 

reg(J+l>•reg(J+l)+anarea*stnoc 

J=J-1 

if(iswd. eq.O> go to 2060 
w 

"" 

f.l) 

II) 

!;;" 

.. 

!;;" 

~ 

> 

~ ~ 

ii 

.. 

"'1:1 

~ 3 

Q 

.. 

~ 

~· 

;::;· 

:t 

AI 

tl 

a. 

m 

~ 

~ g· 

arerr=Capprox-sarea(k))/sarea<k> 

allerr=allerr+1 . 

if Cabs Car err>. 1 t. . 05) go to 2090 

tot5=tot5+1. 

if'Cabs(arerr>. le .. 10> go to 2090 

writeC16,2100) arerr, ira, ii,rc.rf',k, 11, in,dumid 

go to 2090 

2060 iqs=iq_ 

need to add rema1n1ng area if any to lowest acceleration 

dif•sareaCk>-approx 

if ( d i f. 1 t. 0) g 0 t 0 2080 

2070 regC1)=regC1)+dif*stnoc 

2080 approx=sareaCk> 

2090 continue 

2100 format(' int err='f'6. 3' row 'i2' col 'i2' rc=' 

8c f7. 2' rf='f7. 2' k='i2' ll='i2' in='i1, 1x, a4) 

2110 continue 

go to 2730 

c fault computation 

2120 do 2720 J=1, Jtot 

JSeg=Jm(J) 

JSegm=Jseg-1 

call cldisCxmid, ymid, J) 

do 2160 11•1, lev 

fill in rupture lengths--make certain that rupture length 

c does not exceed total fault length 

nrl=nrls 

do 2130 irl=l,nrl 

r 1 1 < i r 1 , 11 > =r 11 s < i r 1 , 11 > 

2130 prlrl<irl, ll>=prls<irl> 

if < r 11 < 1 .11 >. 1 t. d tot < J » g o to 2160 

r11C1. ll>=dtotCJ) 

if' <nrl. eq. 1 > go to 2160 

2140 ifCrllC2, 11>. lt. dtot<J» go to 2160 

prlrl<1, ll>=prlrl<1, ll>+prlrl<2, 11> 

nrl•nrl-1 

if<nrl. eq. 1 > go to 2160 

do 2150 i•2,nrl 

prlrl< 1, 11 >=prlrl < 1+1, 11 > 

2150 rll(i, ll>=rl1Ci+1, 11> 

go to 2140 

2160 nrlsv(ll)•nrl 

ptl•sqrtCshdis> 

ifCptl.lt.reggCidrc>> reggCidrc>•pt1 

if Cisid. eq.O> go to 2490 

distance is monotonically increasing for all segments 

<the decreasing case has been converted to increasing> 

or single shortest distance is interior to fault 

do 2480 11•1, lev 

nrl•nrlsvCll> 

frcon•dtot(J)/dsum*noc<ll> 

maxm=maxaCll> 

J lo•indptCdtab C 1, 11 ), pt1, maxm> 

if (Jlo.gt.O> go to 2170 

c distance beyond table--lump acceleration in first box 

reg<1>•reg<1>+frcon 

go to 2480 

2170 continue 

if < isid. eq. 1 > go to 2190 

2180 ibeg=1 

ifin=Jsegm 

ixt=1 

dnear=dtotCJ) 

inear=Jsegm 

go to 2310 

2190 ibeg=1 

ifin=inear 

ixt=2 

dnear=distance along fault from one end to point on fault 

closest to the site 

dnear•total length of' segments 1, 2, ... inear 

inear defined sub cldis--in common/extra/ 

dnear=O. 

do 2200 iq=1, inear 

2200 dnear=dnear+distlCiq) 

dfar=dtot(J)-dnear 

totac=O. 

do 2220 irl=1,nrl 

smd=r 11 < i r 1, 11 > 

if<smd.lt. dtot(J)) go to 2210 

break is entire segment 

totac=totac+frcon*prlrlCirl, 11> 

sumr=totac 

go to 2220 

2210 smd2=smd/2. 

c find fraction of fault at closest distance to site 

2220 

for interior site 

arcn=smd2 

arcf=smd2 

ifCdnear. lt. smd> arcn=dnear-smd2 

ifCdfar. lt. smd) arcf=dfar-smd2 

arct•arcn+arcf 

if(arct. lt .. 0001> go to 2220 

totac•totac+frcon/Cdtot<J>-smd>*arct*prlrlCirl, 11> 

continue 

sumr=totac 

totacs=totacs+totac 

Jlop=Jlo+l 

ifCtotl. gt.O> go to 2230 

spread dummy faults over "totl" km distance if "totl" > 0 

if "totl" = o, spread over magnitude instead 

call spread<reg, Jlop, 11, totac, pt1> 

go to 2310 

2230 width=totl/2. 

sumr==sumr/totl 

add=O 

ifCabsCptl>. gt.width> go to 2280 

q_tot=width+absCpt1> 

dsvs=O. 

Jlo2=maxm 

do 2270 iqk=1,2 

Jrp•Jlo2 

do 2250 Jk•1, Jlo2 

dl3 =dtabCJrp, 11> 

if(dl3.gt. qtot> go to 2260 

cary=dl3-dsvs 

reg(Jrp+1)=reg(Jrp+1)+sumr*cary 

add=add+sumr*cary 

Jrp=Jrp-1 

2240 format<' add•'e12. S> 
f 

a. 

i(" 

n 

~ 

"" 

2250 dsvs•dl3 

Jrp=O 

2260 car~=~tot-dsvs 

add=add+sumr*car~ 

Jrpp•Jrp+1 

reg(Jrp+1)arag<Jrp+1>+sumr*cary 

dsvs=O. 

~tot=width-abs<pt1> 

2270 Jrp•Jlo2 

Jrpp==Jrp+1 

go to 2310 

2280 dstart=absCpt1>-width 

~tot•totl+dstart 

Jlo2=indpt<dtab C 1, 11 ), dstart, maxm> 

dsvs=dstart 

Jrp=Jlo2 

do 2290 Jk•l, Jlo2 

dl3=dtab<Jrp, 11> 

if(dl3.gt. ~tot) go to 2300 

cary•dl3-dsvs 

reg(Jrp+l>=reg(Jrp+1)+sumr*cary 

add=add+sumr*cary 

Jrp=Jrp-1 

2290 dsvs•dl3 

Jrp=O 

2300 cary•~tot-dsvs 

reg(Jrp+l)=reg<Jrp+1)+cary*sumr 

c find frac of fault in distance range dtabCJr,ll> < d <dtab(Jr-1, 11> 

c for each magnitude 11 

2310 do 2470 ithru=1,ixt 

Jr=JlO 

ilo•1 

iii=ibeg 

casum=O. 

irl•1 

dsav•-dlsCiii) 

2320 frac•O. 

Jk•l 

smd=rllCirl, 11> 

if C smd. lt. dtot < J » go to 2330 

if<ixt.e~.2>go to 2450 

c add in entire fault length 

totac=frcon*prlrl<irl.ll> 

Jlop=Jlo+l 

call spread Crag, Jlop. 11, totac, ptl> 

go to 2450 

2330 cuml•dnear-smd 

if<cuml. le. O> go to 2450 

do 2340 i~=irl,nrl 

2340 frac=frac+prlrl<i~, 11>/CdtotCJ>-rll<i~, 11)) 

frac•frac*frcon 

dlen•O. 

do 2350 i~•ibeg.ifin 

disrem<i~>=distl<i~) 

dlen•dlen+disremCi~) 

if (dlen.ge.cuml> go to 2370 

2350 continue 

write<16,2360)dlen,ibeg, ifin.cuml 

2360 format<' error for fault line dlen•'e12.5' ibeg='i3 

&' ifin='i3' cuml='e12.5> 

call exit 

2370 disrem(i~>=cuml-dlen+disrem<i~) 

2380 ifins=i~ 

2390 dend=disrem<iii)-dlsCiii) 

c dsav=starting distance along segment iii 

2400 continue 

do 2410 Jk=ilo, Jlo 

dl3=s~rt(dtab(Jr• 11>**2-pers~Ciii)) 

ifCdl3. gt. dend> go to 2440 

carea=dl3-dsav 

casum=casum+carea 

reg(Jr+1>•reg(Jr+1>+carea*frac 

Jr=Jr-1 

2410 dsav=dl3 

c at top of table--remaining dist beyond table 

2420 do 2430 i~=irl.nrl 

cuml-dnear-rll < i~, 11 > 

2430 reg<1>=reg<1>+<cuml-casum)*prlr1Ci~, ll>*frcon 8c I C d tot C J ) -r 11< i ~, 11 )) 

go to 2460 

2440 carea=dend-dsav 

2450 

2460 

2470 

2480 

casum=casum+carea 

regCJr+1>•reg(Jr+1)+carea*frac 

dsav=dend 

at end of segment iii 

if< iii. e~. ifins> go to 2450 i ii•i ii+1 

dsav=-dls(iii) 

ilo•Jk 

go to 2390 

irl=irl+1 

ilo=Jk 

if(irl. le. nrl> go to 2320 

if( ithru. e~. ixt> go to2480 

irl=1 

ibeg•inear+1 

ifin=Jsegm 

dnear=dtotCJ)-dnear 

continue 

continue 

go to 2720 

more than one turning point 

brute force. find contribution for single break 

c 

2490 

increment break centers by distance deltal along fault from one end to other. 

continue 

c 

del tal•10. 

ifCpt1. gt. 100) deltal=20. 

ifCptl. lt. 15. > deltal=5. 

do 2710 11=1. lev 

nrl=nrlsvCll> 

do 2710 irl=l,nrl 

smd •r 11 C i r 1, 11 > 

prnoc=noc C 11 >*prlrl C irl, 11 > 

dJ=dtot<J>-smd 

nd=dJ/deltal+1. 5 

ifCnd. gt. 1> go to 2500 

nd=2 

ifCdJ. gt. 1. > go to 2500 

break is entire fault length 
(N 

~ 

(II 

tD 

iii" 

... 

iii" 

Ill:" 

> 

~ 3 

i 

;-

... 

"'1:1 

a 

OQ 

Dl 

3 

§' 

~ iii" 

3 

;:;· 

:z: 

$ 

a. 

1%' 

§· 

~ 

~· 

2500 

2510 

2520 

2530 

2540 

2550 

c 

2560 

c 

2570 

2580 

maxm=maxa(11) 

J1o=indptCdtabC1, 11), pt1, maxm> 

add=prnoc*dtotCJ)/dsum 

Jlop=Jlo+1 

call spread <reg, Jlop, 11, add, pt1 > 

go to 2710 

continue 

dmvp=prnoc*dtot(J)/Cdsum*<nd-1)) 

del=dJ/(nd-1> 

do 2700 m=1,nd 

dst=<m-1>*del 

do 2510 i=1, Jsegm 

if<dst. lt. dist(i, J) > go to 2530 

dst=dst-dist<L J) 

continue 

writeC16,2520)dst, J 

format(' apparent error fault dist--dst='e12. 5' J='i4> 

go to 2990 

ddsav=1. e8 

dfin=dst+smd 

dover=O. 

if(dfin. le. dist(i, J» go to 2550 

dover=dfin-distCi, J) 

dfin=dist<L J) 

ifCdls(i). 1t. dst> go to 2560 

ifCdls(i>. gt. dfin) go to 2570 

shortest distance to segment at point within segment 

dd=persq_Ci> 

go to 2580 

shortest dist at start of segment 

dd=Cd1s(i)-dst>**2+persq_(i) 

go to 2580 

shortest distance at end of segment 

dd=Cd1s(i)-dfin>**2+persq_(i) 

continue 

if(dd. 1t. ddsav> ddsav=dd 

if(abs<dover>. 1t. 0. 001> go to 2600 

dst•O 

i=i+1 

dfin=dover 

if(i.le. Jsegm> go to 2540 

2590 format(' i too large i='i2) 

write<16,2590) i 

call exit 

2600 dis•sq_rt<ddsav> 

have distance to subsegment of arc of fault 

if <m. gt. 1) go to 2610 

J en 2• in d p t < d tab < 1 , 11 > , d i s , ma x a < 1 1 ) ) 

disav•dis 

go to 2700 

2610 disdif=dis-disav 

ifCdisdif. eq_.O) go to 2620 

d fac=dmvp/d isd if 

spread out contributions into accel boxes 

ifCdisdif) 2630,2620,2660 

2620 J1op=Jen2+1 

call spreadCreg, J1op, ll,dmvp,dis) 

go to 2700 

distance decreasing 

:2630 Jen2=Jen2+1 

2640 distab=dis-dtabCJen2, 11> 

if (distab. 1t. O> go to 2650 

regCJen2)=regCJen2)+(dis-disav)*dfac 

add=Cdis-disav)*dfac 

if(add. lt.O> print *•2740 

Jen2=Jen2-1 

disav=dis 

go to 2700 

2650 regCJen2>=reg(Jen2>+Cdtab(Jen2, 11>-disav>*dfac 

add=Cdtab(Jen2, 11)-disav>*dfac 

ifCadd. lt. O> print*, 2750 

d isav=dtab < Jen2, 11 > 

J en2=J Em2+1 

go to 2640 

2660 if(Jen2. eq_.O> go to 2680 

2670 distab=dis-dtabCJen2, 11) 

ifCdistab. gt.O> goto 2690 

2680 regCJen2+1)=regCJen2+1)+Cdis-disav)*dfac 

add=Cdis-disav)*dfac 

if Cad d. 1 t. 0) print *, 2780 

disav=dis 

go to 2700 

2690 regCJen2+1)=regCJen2+1>+Cdtab(Jen2, 11>-disav>*dfac 

add=Cdtab(Jen2,11>-disav>*dfac 

if(add. lt. 0) print*, 2790 

disav=dtab(Jen2, 11> 

Jen2=Jen2-1 

if CJen2. ge. 1) go to 2670 

go to 2680 

2700 continue 

2710 continue 

c above ends loop for more than one turning point 

2720 continue 

2730 continue 

do 2740 iq_r=1,maxdim 

2740 rawbin(iq_r, idrc>=regCiq_r) 

save accelerations for this source & site in rawbin array 

ends loop for one site, one area or fault 

repeat for all sites for this source before going on to next c source 

2750 continue 

2760 continue 

c one source completed for all sites 

2770 continue 

c----smooth here if sigma for earthq_uake locations has changed 

c----from previous value 

if( line. eq_. 1> go to 2850 

c-----no smoothing for faults 

ifCnumsav.ne.98) go to 2850 

if sigma=O, no smoothing 

ifCals. eq_.O) go to 2850 

ifCnrc. eq_.O> go to 2850 

c-----

ismsw=l 

sigd=Cdphinc/als>**2 

do 2780 im=1, imaxsq_ 

2780 fac(im)=expCC-im*sigd)/2. 

kmax=10 
> i 

i 

:;;c· 

r"' 

~ 

...... 

c 

eq,sum=O. 

dumfac=l. 

compute icon based on sigma and dphinc 

probably enough to go out 2 sigma on either side of point 

find how many rows or columns that corresponds to 

icon=2*als/dphinc+1.01 

if (icon. 1 t. 5) icon=5 

if' icon<5. set icon=5 

irm=irow2-nbord 

icm=icol2-nbord 

isubl=irowl+nbord 

isub2=icol1+nbord 

do 2830 iia=isub2, icm 

J l=i ia-icon 

J2=iia+icon 

if< J2. gt. icol2> J2= col2 

if(Jl. lt. icoll> J1= colt 

do 2830 iraa=isub1. rm 

i1=iraa-icon 

i2=iraa+icon 

if( il. lt. irowl> i1=irow1 

i2=iraa+icon 

if( i2. gt. irow2> i2=irow2 

do 2800 iq=l.maxp 

add=O. 

wtt=O. 

d 0 2790 i = i 1. i 2 

d 0 2790 J = J 1. J 2 

ksq,=<iraa-i)**2+(iia-J)**2 

if( ksct. gt. imaxsct> go to 2790 

we need to convert to proper subscript in rawbin 

in rawbin. first all cols for one row 

<i-1>*icol(s)+J=row and col with single subscript in rawbin 

isub=<i-irowl>*icols+J-icol1+1 

add=add+f'ac<ksct>*rawbin(iq, isub> 

wtt=wtt+f•c<ksct> 

2790 continue 

reg<ict>=O 

if'<wtt.ne.O> reg<ict>=add/wtt 

ectsum=ectsum+reg<ict> 

2800 continue 

2810 continue 

isub=(ir•a-irow1>*icols+iia-icol1+1 

do 2820 ict=1.maxdim 

2820 smobin<ict• isub>=reg<ict>+smobin<ict• isub) 

2830 continue 

c now delete entries from rawbin array 

do 2840 isub=1.nrc 

do 2840 ict=1.maxdim 

2840 rawbin<ict• isub>=O. 

2850 continue 

if!(iprint. lt. O> go to 2940 

c output all results to this point 

do 2930 ispt=1.2 

if' < ispt. ect. 1) go to 2860 

if'(indvpt. ect.O> go to 2930 

irow1s=1 

irow2s=1 

icolls=l 

icol2s=indvpt 

go to 2870 

2860 if'(irowl. ect.O> go to 2930 

irow1s=irow1+nbord 

irow2s=irow2-nbord 

icolls=icoll+nbord 

icol2s=icol2-nbord 

2870 do 2920 ira=irow1s. irow2s 

do 2920 ii=icol1s, icol2s 

iF(ispt. •ll· 1) go to 2890 

2880 xouts=xpt1<ii) 

youts=yptl<ii) 

isub=irows*icols+ii 

go to 2900 

c back transform point from ectuator to dec degrees 

2890 call backw<ii. ira. xouts.youts> 

results calculated and printed in subroutine out 

isub=<ira-irowl>*icols+ii-icol1+1 

2900 do 2910 ictr=l.maxp 

2910 reg<ictr>=rawbin<ictr• isub>+smobin<iqr. isub) 

rdis=regg(isub) 

i F < i p r in t. 1e. 2 > wr i t e < 16, 30 > t y t 1 e 

call out(reg. naf, xouts. youts, sd) 

2920 continue 

2930 continue 

2940 continue 

isw=l 

go to 680 

c Finished 

2950 nvals=icols*irows+indvpt 

write to unit 3 for possible run continuation 

may delete to save time and space if' no run 

continuations are desired. 

do 2970 i=l.nvals 

do 2960 ictr=l.maxdim 

2960 reg<ictr>=rawbin<ictr• i)+smobin<ictr• i) 

write<3> reg 

2970 continue 

c transfer to unit 3 for possible continuation or restart 

2980 continue 

2990 continue 

c write(16.2530) tot5.allerr 

3000 Format(/' sources with error in area> .05='f7.0 

~ 'out of 'f7.0' sources ') 

3010 format<' number of brute force points='i4) 

stop 

end 

subroutine cond<x1.x2.yl.y2> 

call condec<xl> 

call condec<x2> 

call condec<yl> 

call condec <y2> 

return 

end 

subroutine condec<phi> 

i=phi 

c Fr=phi-i 

phi=float<i>+<phi-float(i))/.60 

return 

end 

c *************************************************************** 
~ 

= 

~ ;;· 

.. 

;;· 

~ 

> 

~ Q 

3 

i 

tD 

.. 

""C' 

a 

OQ 

iii: 

3 

§' 

~ ;;· 

3 

;:;· 

:::t: 

~ a. 

; 3. 

~. 

~ 

subroutine lin< x1, 1,11. x2, 1,12. a, b, c. x3. 1,13> 

c coefficients a. b. c are determined so that points which are in the 

same half plane as x3,\13 will be a positive distance from the line 

c JOining x 1. 1,11 and x2. 1,12 

c a x + b \1 + c=O 

if(xl. eq,. x2) go to 10 

a=<~,~1-1,12)/(x1-x2> 

b=-1 

c=~,~1-a*x1 

go to 20 

10 a=l. 

b=O 

c=-xl 

20 d=a*x3+b*1,13+c 

if<d.gt.O. >return 

a=-a 

b=-b 

c=-c 

return 

end 

c *************************************************************** subroutine linint<reg, t, sol> 

calculates acceleration "sol" that has probabilit\1 

"prob" of not being exceeded during "t" 1,1ears at 

site. 

expected 1,1earl1,1 accelerations in the range 

ac(J-1> <a < ac(J) are in reg(J). 

linint is called from subroutine out. 

common/tims/dumm\1<42), iprint 

common/quad/maxp.ac(105),dumlog.aclog(105) 

common/ext/prob.ex.exl 

dimension reg(106) 

ii=maxp-1 

sum=O. 

10 if( ii. le. 2> go to 20 

if(reg<ii>. gt. 0.) go to 40 

ii=ii-1 

go to 10 

20 sol=O. 

if<iprint.le.2> write<16,30) 

30 format<' sol not obtained for time='f7.0) 

return 

40 sum=sum+reg<ii) 

if<t*sum. gt.ex) go to 50 

i i=ii-1 

if(ii. ge.2> go to 40 

go to 20 

50 if<sum. eq.reg(ii>> go to 70 

1,12=alog<t*<sum-reg<ii>>> 

l,ll=alog<t*sum> 

exlt=exl 

60 a=<l,ll-1,12)/(ac<ii-1>-ac::<ii))' 

b=\11-a*ac(ii-1> 

sol=<exlt-b)/a 

return 

70 1,1l=t*sum 

1,12=0. 

exlt=ex 

if(iprint.le.2> write(16.80> 

80 format<' log interpolation fails--use linear'> 

go to 60 

end 

c *************************************************************** subroutine box<scale,sd) 

determine boundaries ac(i) for acceleration levels 

accelerations in the range ac<i-1> <a <ac<i> will be 

accumulated in reg(i) 

scale multiplies the basic levels computed 

scale=! for accelerations in the range .02 to 1. 

common /q,uad/maxq,,ac(105>.dumlog.aclog<105> 

common/wds/max.maxp.maxsp 

maxsp=104 

max=54 

maxp=55 

ac<5>=.02*scale 

dif=ac<5> 

do 10 i=6,55 

ac<i>=ac<i-l)+dif 

10 continue 

c -----------------note important--------------

the following inserted to allow acceleration boxes for small 

values below range of interest for sd computation 

this should be changed later as desired 

J=4 

efac=exp <-sd/2. ) 

do 20 i=1,4 

ac<J>=ac(J+l>*efac 

20 J=J-1 

for later attenuation variabilit\1, expand acceleration range 

used onl\1 in subroutine "out" 

l,l=alog(ac(max>> 

w=exp(1,1+4.4*sd) 

wl=exp(l,l+2.0*sd) 

dimension of ac, aclog is 105; we cannot exceed this dimension 

n=maxsp-max 

n2=n/2 

i2s=max+n2 

dif=<w1-ac(max>>ln2 

do 30 i=maxp, i2s 

ac<i>=ac<i-l)+dif 

30 continue 

i2sp=i2s+1 

dif=<w-ac<i2s))/n2 

do 40 i=i2sp.maxsp 

ac<i>=ac(i-l)+dif 

40 continue 

return 

end 

c ************************************************************** subroutine inieqr<xl,\11• x2.1,12) 

read xlong, 1,1lat in radians for 1st point 

read xlong, 1,1lat in radians for 2nd point 

points 1 and 2 to be transformed to equator 

point 1 will be transformed to <O,O> 

point 2 will be transformed to (d,O> 

where d•great circle distance between <x1.\11) and <x2.1,12) 

output is transformation matrices f (forward) fi (backward) 

generated and stored in common 

common/eq,r/ r1(3,3),r2(3,3),r3<3,3),d(3,3),f(3,3),vec<3> 

common/inv/fi<3,3),res(3) 
~ 

l >C. 

~ 

~ 

~ 

cp1=cos(y1) 

cp2=cos(y2) 

cl1=cos<x1) 

sp1=sin<y1> 

sp2=sin(y2) 

sll=sin< x 1) 

initialize matrices 

do 10 i=1.3 

do 10 J=1.3 

r 1 ( i, J) =0. 

r2( i. J >=0. 

10 r3<i•J>•O. 

rl< 1. 1>=1. 

r2<2• 2>=1. 

r3(3,3)=1. 

r2<1• 1>=cp1 

r2(3,3>=cp1 

r2<1.3)•sp1 

r2(3, 1>=-sp1 

r3< 1. 1 >=c 11 

r3(2,2>=cl1 

r3( 1. 2>=sl1 

r3(2, 1>=-sl1 

s1=cp1*sp2-sp1*cp2*cos<x2-x1> 

c1=cp2*sin<x2-x1> 

den=s~rt<c1*c1+sl*s1> 

sinth=s1/den 

costh=c1/den 

r1(2,2)=costh 

r1<3.3)=costh 

r1<2.3>=sinth 

r1(3.2>=-sinth 

call matmpyCr1,r2,d) 

call matmpy(d,r3,f) 

transformation matrix in f! inverse matrix in fi 

compute inverse 

do 20 i=1.3 

do 20 J=1,3 

r 1 ( i. J) =0. 

r2< i· J >=0. 

r3< i· J >=0. 

20 continue 

r1 (3, 3>=1. 

r2<2• 2>=1. 

r3( 1, 1 >=1. 

rl< 1. l>=c 11 

r1 (2, 2)=c 11 

r1 <2. 1 >=sl1 

r1(1,2)=-sl1 

r2< 1. 1 >=cp1 

r2(3,3)=cp1 

r2C1.3)=-sp1 

r2 < 3, 1 ) =s p 1 

r3C2.2>=costh 

r3C3.3)=costh 

r3(3,2>=sinth 

r3(2,3>=-sinth 

call matmpy(rl.r2.d) 

call matmpyCd.r3.fi) 

return 

end 

c *************************************************************** 

subroutine toe~r<x3.y3,flam6,phi6> 

transform point at <x3.y3) <long. lat in radians> to new 

equator using transformation defined in subroutine ineqr 

common/eqr/ r1(3,3>.r2C3.3),r3(3,3),d(3,3),f(3,3),vec<3> 

common/inv/fi(3,3),res(3) 

c xmid=cos ( x3> 

cymid=cosCy3> 

sxmid=sin(x3) 

symid=sin(y3) 

vec(1)=cxmid*cymid 

vec<2>=sxmid*cymid 

vec<3>=symid 

call vecmpy(f,vec.res) 

phi6=asin<res(3)) 

cp6=cos(phi6) 

cl6=res(1)/cp6 

sl6=res(2)/cp6 

if(cl6. gt. 1. > cl6=1. 

flam6=acos(cl6) 

if<sl6. lt. 0.) flam6=-flam6 

return 

end 

c *************************************************************** 

subroutine backw<ix. iy. flii. phii> 

for printout and identification: point defined by 

row ix and col iy in new coordinate system is transformed 

to location <long, lat. in dec degrees> in original 

coordinate system. 

common/inv/fi(3,3),res<3> 

dimension vec(3) 

c need cos. sins in common array 

common/sincos/sinxC500>.cosxC500>.siny(500),cosy<500) 

data rad/.0174533/ 

sl6=sinxCix) 

cl6=cosxCix> 

cp6=cosyCiy> 

sp6=sinyCiy) 

vec<1>=cl6*cp6 

vec<2>=sl6*cp6 

vec(3)=sp6 

call vecmpyCfi,vec.res> 

if < a b s ( r e s ( 3 ) ) . 1 t. 1. 0 ) g o to 30 

if<res(3). ge. 1.) go to 10 

phi i=-90. 

go to 20 

10 phi i=90. 

20 flam=O. 

return 

30 phii=asin(resC3)) 

cpi=cosCphii) 

phi i=ph i ilrad 

cli=res(l)/cpi 

sli=res(2)/cpi 

flii=acos(cli) 

if(sli. lt. 0. > flii=-flii 

flii=fliilrad 

return 

end 
... 

= 

~ ij;• 

.. i;;. 

Ill:" 

> 

Q 

3 

"C 

c 

;-

.. 

"1::1 

i 

0' 

.. 

~ ij;• 

3 

;:.;· 

::z: 

E 

a. 

~ 

~r 

!:. §. 

c *************************************************************** 

subroutine matmpy<a.b.c) 

multiply two 3 x 3 matrices A x B = C 

dimension a(3,3),b(3,3),c(3,3) 

do 10 i=1. 3 

do 10 J=1.3 

c ( i. J >=0. 

do 10 k=l. 3 

10 c(i, J)=c<i. J)+a(i, k)*b(k, J) 

return 

end 

c *************************************************************** 

subroutine vecmpy<f.vec.r> 

vector multiplication used by subroutines backw, toeqr 

dimension f(3,3),vec(3),r(3) 

do 10 i=1.3 

r< i >=0. 

do 10 k=1.3 

10 r<U=r<i>+f(i, k)*vec<k> 

return 

end 

c **************************************************************** 

subroutine csqdis<x, Y• xl. yl, xr. yr. isub) 

for a quadrilateral source area. calculates 

nearest distance from site at <x.y) to each side 

and distance to each vertex. 

called from subroutine rrisk 

dimension x 1<4), y 1<4), xr <4>, yr <4> 

common/slinp/fsa(4,50),fsb(4,50),fsc<4,50),rtaC4.50).dist<4,50) 

8c .disq(4.50) 

common /slinds/dls(4).perpC4),sl1(4),slsq(4),minC4).persq<4> 

s., srsq<4>. dlp<4> 

do 10 i=1,4 

10 perp<i>=<fsa(i, isub>*x+fsb(i, isub>*y+fsc(i, isub> )/rta(i, isub> 

do 20 i=1.4 

slsq<i>=<x-xl(i))**2+<y-yl(i))**2 

srsq<i>=<x-xr(i))**2+<y-yr(i))**2 

20 continue 

do 50 i•1.4 

persq<i>=perp<i>*perp(i) 

dl1=slsq(i)-persq(i) 

dl2=srsq(i)-persq(i) 

if(dll. lt. 0.) dll=O. 

if(dl2. lt. 0.) dl2=0. 

ii'Cdl1+dl2. le. disq(i, isub)) go to 40 

shortest distance external to segment 

if<dll. lt. dl2) go to 30 

min(i)•1 

dlp(i)=sq,rt<dl2) 

d 1 s ( i ) =dis t ( i. is u b ) +d 1 p ( i ) 

go to 50 

30 min(i)=2 

dls(i)=-sqrt<dll) 

go to 50 

40 min< i >=3 

dls<i>=sqrt<dl1> 

dlp(i)=O. 

50 continue 

return 

end 
f 

!. 

>e· 

n 

!: 

c *************************************************************** 

c 

c 

c 

subroutine rrisk < xnot, ynot, in, isub, r. rc. rf, itst > 

commonlinoutlic, ics, ifar.azimf,stsize 

commonlxyarealxsav(2,50),ysav<2,50),sarea(50) 

commonlslindsldls(4),perp<4>,sl1<4>.slsq(4),min<4>,persq<4> 

&c. srsq(4), dlp(4) 

commonlslinplaa(4,50),bb(4,50>.cc(4,50),rta(4,50),dist(4,50), 

lie disq<4,50> 

commonlxlyllxl(4), yl(4), xr<4>, yr<4>, xc<4>, yc<4> 

subroutine which loads temporary arrays with this 

subsource's corners and determines if this site is 

inside or outside the subrource. 

data pi216.28318531 

fixed number of radii at which distance-vs arc length 

table is computed is set here itst=7 *************** 

itst=7 

rc2=10000000. 

ics=O 

call csqdis<xnot, ynot, xl, yl, xr. yr. isub) 

rf2=0. 

do 20 ii=1,4 

if(slsq(ii). gt. rc2) go to 10 

rc2=slsq( ii) 

10 if(slsq(ii). lt. rf2> go to 20 

rf2=slsq<ii) 

ifar=ii 

20 continue 

30 

40 

50 

60 

70 

80 

do 30 ii=1,4 

if<min(ii>.ne.3> go to 30 

if<persq<ii>. gt. rc2> go to 30 

rc2=persq(ii) 

ics=ii 

continue 

rc=sqrt<rc2> 

rf=sqrt<rf2> 

stsize=<rf-rc>litst 

step thru source area. 

r=rc+(O. 5)*stsize 

determine if site is inside source area. 

do 50 ii=1,4 

d=aa(ii, isub>*xnot+bb(ii. isub)*ynot+cc(ii. isub) 

if(d. lt. O>go to 90 

continue 

in=l 

determine azimuth of farthest point with respect to site 

azimf=acos((xl(ifar)-xnot)lrf) 

if <yl<ifar>-ynot) 60,70,70 

azimf=pi2 - azimf 

continue 

rc is now distance from site to closest side. 

rf is now distance from site to farthest corner. 

pick step size based on fraction of area left 

return 

if do loop finished, point lies within area. 

90 in=O 

100 return 

end 

************************************************************** subroutine setreg<isub) 

find equation of line thru each side of quadrilateral 

c 

c 

find length of each side, etc. 

commonloutslxouts.youts 

dimension xl (4), yl<4>, xr<4>, yr<4> 

commonlxyarealxsav(2,50),ysav<2,50),sarea(50> 

commonlslinplaa(4,50),bb(4,50),cc<4,50>.rta(4,50),dist<4,50), 

lie disq(4,50) 

data rl6378. I 

do 10 i=l. 2 

xl<i>=xsav(i. isub> 

xl<i+2>=xsav(i, isub+l) 

yl<i>=ysav(i. isub> 

yl< i+2>=ysav< i, isub+l> 

10 continue 

xr<3>=xl<l> 

yr<3>=yl<1> 

xr<l>=xl<2> 

yr<l>=1Jl<2> 

xr<2>=xl<4> 

yr<2>=yl<4> 

xr<4>=xl<3> 

yr<4>=yl<3> 

determine if any sides are vertical lines 

do 70 ii=1.4 

H<xl<ii>. eq. xr(ii)) go tc 20 

a a ( ii, is u b > = < y 1 < i i > -y r < i i > >I< x 1 < i i >- x r < ii > > b b ( i i. i sub ) =-1 . 

cc<ii, isub)=yl<ii>-aa(ii. isub>*xl<ii) go to 30 

20 aa<ii, isub)=l. 

bb<ii, isub>=O. 

C C ( i i • iS U b ) =-X 1 ( i i ) 

30 if ( i i. g t. 2) g 0 t 0 40 

d=aa<ii, isub>*xl<3)+bb<ii. isub>*yl<3>+cc<ii. isub) go to 50 

40 d=aa<ii, isub>*xl<2>+bb(ii. isub>*yl<2)+cc<iL isub> 

50 if ( d . g t. 0. ) g 0 t 0 60 

aa<ii, isub>=-aa(ii, isub> 

bb(ii. isub>=-bb(ii. isub> 

cc(ii. isub>=-cc(ii, isub> 

60 disq<ii, isub>=<xr(ii>-xl<ii>>**2+<yr<ii)-yl(ii))**2 

dist(ii. isub>=sqrt<disq(ii. isub)) 

70 rta<ii, isub)=sqrt<aa<ii, isub)*aa(ii, isub)+bb(ii, isub)* 

llcbb<iL isub)) 

return 

end 

function indpt<y.yval.n) 

points in descending order--non zero values from points 

1 to n: find first value in table exceeding yval 

dimension y<102) 

H<yval. lt. y<1» go to 10 

indpt=O 

return 

10 if<yval.gt.y<n» go to 20 

indpt=n 

return 

20 num•n 

idiv=2 

30 nh=numlidiv 

if(nh. le. 1> go to 50 

if<yval. lt. y<nh)) go to 40 
~ 

~ 

~ iii' 

~. 

!Q-

> 

Q 

3 

~ ;-

., 

., 

a 

OQ 

D'S 

3 

S' ., 

~ iii' 

3 

;:;· 

:I: 

~ a. 

~ 3' 

~ j' 

idiv=2*idiv 

go to 30 

yval lies approx between nh and 2*nh or 2*nh 

40 idiv=2*idiv 

nh2=nh+num/idiv 

if<nh. eq. nh2> goto 80 

if< yval. gt. y (nh2> > go to 40 

yval further down on table than y(nh2> 

nh=nh2 

go to 40 

50 do 60 i=l. n 

if<11<i>. lt. 11va1> go to 70 

60 continue 

70 indpt=i-1 

return 

80 continue 

do 90 i=nh.n 

if(y(i). lt. yval> go to70 

90 continue 

end 

c ************************************************************** subroutine cldis(x, y, Jl) 

subroutine to compute shortest distance. etc to line 

containing fault segment i of Jl th fault <for each i) 

determines whether closest point to fault is at 

low end <min(i)=2> 

interior <min(i)=3> 

high end <min(i)=1> 

set isid=O if multiple turning points 

isid=1 if single closest distance from site to fault 

interior to fault 

isid=2 if closest distance from site to fault at one end 

of the fault 

dimension dlp<20) 

common/linps/distl(24,26),disq(24.26), JSeg, JSegm, 

& xl<24,26>•11P<24,26),coa<24,26),gia(24,26l 

common/extra/isid.shdis,dist(24), inear 

common /lindis/dls<24l.perp<24l,min(24),persq<24) 

shdis=l. e9 

do 10 i i=1, Jsegm 

10 dist<iil=distl<ii. Jl> 

i1=0 

i2=0 

i3=0 

isid=O 

do 20 i=l, JSIIg 

slsq=<x-xl(i, Jl>>**2+<y-yp<i. Jl)l**2 

if(slsq. It. shdis> shdis=slsq 

20 continue 

do 50 i=l,Jsegm 

xd=x-xl(i, Jl> 

\jd=y-yp(i, Jl> 

perp<il=-xd*sia(i, Jl>+yd*coa<i. Jl> 

dls<il=xd*coa(i, Jl>+lld*sia(i, Jl> 

persq(i)=perp(i)*perp(i) 

if(dls(i). lt. O> go to 30 

if(dls(i). le. dist(i)) go to 40 

shortest distance external to segment 

min< il=1 

i1=i1+1 

dlp<i>=dls<i>-dist(i) 

go to 50 

30 min<i>=2 

i2=i2+1 

go to 50 

40 min< i >=3 

i3=i3+1 

if<persq<i>. lt. shdis) shdis=persq(i) 

d 1 p ( i >=O. 

50 continue 

if< i 1. ne. JSegm> go to 90 

iup=Jsegm 

isid=2 

60 iuh=iup/2 

if(2*iuh.ne. iup> dls(iuh+1>=-dlp(iuh+1> 

if<iup. eq. 1> go to 80 

J=iup 

do 70 i=1, iuh 

t=persq<i> 

persq<i>=persq<J> 

persq<Jl=t 

d ls ( i) =-d 1 p ( J) 

dls(J)=-dlp<i> 

t=dist<i> 

dis t ( i ) =dis t ( J) 

min<i>=2 

min(Jl=2 

dist<J>=t 

70 J=rl 

80 return 

90 if<i2.ne. Jsegm> go to 100 

isid=2 

return 

100 if< i3. gt. 1) return 

if(min<1l. eq. 2) return 

iturn=O 

iq=l 

do 110 k=l. JSegm 

if < m i n < k ) . n e. 1l g o to 120 

iturn=iturn+l 

110 continue 

120 if<iturn.ne. ill return 

if(i3.eq.1 .and. min(k>.eq.2> return 

we have single turning point 

isid=1 

inear=k-1 

if<min< k l. ne. 3) go to 140 

inear=k 

adJUSt to 2 sep segments 

J=Jseg 

do 130 i=k, Jsegm 

dls < J >=dls< J-1) 

persq(J)=persq(J-1) 

dist<J>=dist<J-1> 

130 J=J-1 

JSegm=Jseg 

Jseg=Jseg+l 

dist(k+ll=dist(k+l)-dls<k> 

dist<k>=dls(k) 

dls(k+l>=O 
f 

8. 

>C' 

n 

.... 

~ 

140 iup=inear 

go to 60 

end 

c ************************************************************** subroutine pbreak<fm, aL bl, sigL rlls. prls, lev. nrls> 

dimension fmC12>.prlsC12>.rllsC10.12> 

dimension pmid<5>.prang<5> 

data pmid/1. 575 •. 74,0 .• -.74.-1. 575/ 

data prang/. 1151,. 2295,. 3108,. 2295,. 1151/ 

data nrl/5/ 

c evaluate al+bl*fm(l)+pmid(i)*sigl 

c where 1 <1 <lev and 1 < i < nrl 

prang probability associated with Jth break 

c for each magnitude; compute nrl rupture lengths 

c :::::::: must be ordered in descending order--largest break 

::::::::: length first for each magnitude 

do 10 J=1, lev 

con=al + bl*fmCJ) 

do 10 i=1.nrl 

c=con+pmid<i>*sigl 

10 rlls(i, J>=10.**c 

do 20 i=1,nrl 

20 prls(i)=prang<i> 

nrls=nrl 

return 

end 

c ************************************************************* subroutine bbgaus(x,gg> 

c computes normal probability integral gg from -<infinity> to xx 

sets gg=1 if XX < -6; gg=O if XX ) +6 

uses approximation Handbook of Mathematical FunctionsNBS Applied Math Series 55--p299 7. 1. 26 for error function 

ax=absCx)/1.4142136 

if(ax-4.2426408> 20.20,10 

10 gg=O. 

if ( X • 1 t. 0 ) g g =+ 1. 

return 

20 continue 

d=1.+<<<<<4.30638e-5*ax+2.765672e-4>*ax+1. 520143e-4 

& >*ax+9.2705272e-3)*ax+4.2282012e-2)*ax+7.0523078e-2)*ax 

d=d*d 

d=d*d 

d•d*d 

gg=. 5/(d*d) 

if( X) 30, 30,40 

30 gg=1.-gg 

40 return 

end 

c ************************************************************ subroutine outsid < xnot, ynot. isub, r. angle) 

common/inout/ic, ics. ifar.azimf,stsize 

common/xysav/xsav(2,50),ysav(2,50),sarea<50> 

common/xlyl/xl<4), yl<4), xr(4), yrC4h xc(4), yc<4> 

common/slinp I aa(4,SO>.bbC4,50>.ccC4,50),rtaC4,50),distC4,50> 

&.disq_C4,50) 

common/t;linds/dls(4), perp (4), sll <4>, slsq_(4), min<4>, persq<4> 

&,srsq_C4),dlp(4) 

dimension JSUb(3) 

subroutine for calculating angle when site is outside 

c <quadrilateral> source area. 

c 

c 

c 

c 

r is the radial distance for which angle is sought 

<xnot.ynot)=point outside quadrilateral source region isub 

r=radius of circle (or radius of annular ring of width dr) 

centered at <xnot.ynot> 

isub = index identifying quadrilateral source region 

the area contained in the intersection of the region and annular 

ring is determined by the fraction of the ring within the region. 

area = angle x r x dr 

where angle = the angle Cin radians> subtended by the 

part of the annular ring contained within the region. 

Crc is closest distance between site and source.> 

<rf is furthest distance between site and source> 

subroutine csq_dis has computed distance to vertices of subregion 

and shortest distance to each side 

adapted from subroutine by mcguire 

data Jsub/2,3, 1/ 

npt.,O 

angle=O. 

signal=l. 

loop on each side 

rsq_=r*r 

do 140 ii=1,4 

if<min<ii>-2> 10.40.70 

closest point is xr<ii),yr(ii> 

at most one intersection 

10 if <rsq_. gt. slsq_( i i » go to 140 

if(rsq_. lt. srsq_( i i)) go to 140 

d=<sqrt<rsq_-persq_Cii>>-dlp(ii>>ldist(ii· isub) 

20 y1=<yl<ii>-yr<ii>>*d+yr<ii> 

x1=<xl(ii)-xr<ii))*d+xrCii> 

intersects in range 

does circle <radius r> intersect it? 

store first point 

30 npt=npt+1 

xc<npt>=x1 

yc<npt>=y1 

go to 140 

closest point is at xl<ii).yl<ii> 

40 if(rsq_. lt. slsq_(ii)) go to 140 

if C r s 11.. g t. s r s 11. C i i )) g o to 140 

50 d=Csq_rtCrsq_-persq_Cii>>+dls(ii))/dist(ii• isub> 

60 x1=xlCii)+(xrCii>-xlCii>>*d 

y1=ylCii)+(yrCii)-ylCii))*d 

go to 30 

70 if(rsq_. lt. persq_Cii)) go to 140 

i f < r s q_. le. s ls 11. C ii )) go to 80 

i f < r s 11.. g t. s r s 11. < i i » g o t o 140 

single intersection 

on side nearer xr.yr 

go to 50 

80 if(rsq_. gt. srsq_Cii>> go to 90 

two intersections on this side 

can compute angle directly without determining coordinates 

arg=absCperpCii))/r 

as=acos(arg> 

go to 100 
ollo 

ollo 

~ 

tD 

(;;" 

::!. 

Ill 

Jll:"' 

.. 

> 

n 

Q 

3 

"CI 

c 

;-

., 

"'CI 

a 

OQ 

Dl 

3 c 

S' c ., 

~ 

tD 

Iii" 

3 

;:;· 

::t 

AI 

N 

AI 

a. 

m 

~ 3" 

~ 

s· 

:I 

90 d=<dls(ii)-sqrt<rsq,-persq,(ii)) )/dist(ii. isub> 

go to 60 

single intersection 

see if' this side is closest to point. if' so. treat specially. 

100 if! (ii-i c s) 110. 120. 110 

both points are on boundary. calculate angle betwen them . 

110 sign=-1. 

go to 130 

120 sign=l. 

signal=-1. 

130 angle=sign*2.*as +angle 

see if' second point only is on boundary 

140 continue 

if'<npt.gt.O> go to 150 

if'(signal> 310.320.320 

150 go to (320.240, 160.250>.npt 

this is an error unless radius is on a vertex and 

was counted twice. 

160 do 170 i=1.3 

J=JSUb(i) 

if(abs<xc<i>-xc<J» .ge .. 01> go to 170 

if'<abs(yc(i)-yc<J» . le .. 01) go to 220 

170 continue 

check for one of' three points being a vertex 

do 210 i=l. 4 

do 180 J=1.3 

if'<abs<xc(J)-xl<i». gt .. 001> go to 180 

if'(abs(yc(J)-yl(i)) . le .. 001> go to 190 

180 continue 

go to 210 

190 if(J-2) 200.200,230 

200 xc < J >=xc (3) 

yc < J >=yc <3> 

go to 230 

210 continue 

go to 320 

220 if'( i. ne. 1> go to 230 

xc<2>=xc<3> 

yc<2>=yc<3> 

230 npt=2 

go to 240 

240 ad=sqrt<<xc<1>-xc<2>>*<xc<1>-xc<2>>+<yc<1>-yc<2>>*<yc<1>-yc<2>>> 

angle=angle + signal*2.*asin<ad1<2.*r>> 

go to 310 

four intersection points <each on a different side). 

determine angle by finding closest 2 intersections to 

farthest corner. calculate angle between. and add angle 

between other two intersections. 

250 dist1=10000000. 

il=O 

i2=0 

i3=0 

i4=0 

do 300 JJ=1.4 

distn=<xl<if'ar>-xc(JJ))*(xl<if'ar)-xc(JJ)) 

& +<yl<if'ar)-yc(JJ))*(yl(ifar)-yc(JJ)) 

if' (distn-distl) 260.260.270 

260 d i st2=d istl 

dist1=distn 

i4=i3 

i3=i2 

i2=i 1 

il=JJ 

go to 300 

270 if' (distn-dist2) 280,280.290 

280 dist2=distn 

i4=i3 

i3=i2 

i2=JJ 

go to 300 

290 i4=i3 

i3=JJ 

300 continue 

calculate angle between 2 closest points to farthest corner. 

ad=sqrt((xc(il)-xc(i2>>*<xc<i1>-xc(i2)) 

& + <yc(i1)-yc<i2>>*<yc<il>-yc<i2>>> 

angle=2.*asin(ad/(2.*r)) 

calculate angle between 2 points farthest from 

farthest corner and add to previous angle . 

ad=sqrt<<xc(i3>-xc(i4>>*<xc(i3)-xc(i4)) 

& +<yc<i3>-yc(i4))*(yc(i3)-yc<i4>>> 

angle=2.*asin<ad/(2.*r)) +angle 

310 continue 

compute rate of earthquakes in this annular source 

return 

error printout 

320 wr i t e < 16, 330 ) i sub , i c , n p t, x not, y no t , < x 1 < i > , y 1 < i) , i = 1 , 4 > , r , an g 1 e, 

&<xc(i), yc(i), i=1.4> 

330 format<' ***** error in subroutine outsid. source no. '• i3. 

&' debug values follow ..... '• /lOx, 2i5, 5(/lOx, 4f14. 6) > 

340 call exit 

end 

c ************************************************************** subroutine inside(xnot. ynot, isub, r, pangle) 

common/xyarea/xsav<2,5Q),ysav(2,50),sarea<50) 

common/inout/ic, ics, ifar. azimf, stsize, stepo 

common/slinds/dls(4), perp(4), sl1(4), slsq(4), min(4), persq(4) 

c 

&, srsq(4), dlp<4> 

common/slinp/aa(4,50),bb(4,50),cc(4,50),rta<4,50),dist<4•50>, 

& disq(4,50) 

common/x 1111/x 1 <4>, y 1 <4>, xr <4>, yr (4), xc <4 >, yc (4) 

dimension Jsub<3> 

<xmid.ymid)=point inside quadrilateral source region ii 

r=radius of circle <or radius of annular ring of width dr> 

centered at <xmid,ymid) 

ii = index identifying quadrilateral source region 

c the area contained in the intersection of' the region and annular 

c ring is determined by the fraction of' the ring within the region. 

area = pangle x r x dr 

where pangle = the angle (in radians> subtended by the 

part of the annular ring contained within the region. 

subroutine insid determines pangle. 

adapted from subroutine by mcguire 

data pi, p i2/3. 1415927,6. 2831853/ 

data Jsub/2,3, 1/ 

rsq=r*r 

npt=O 
f 

a. 

;;c· 

n 

.,. 

VI 

ang le=O. 

loop on each side 

do 110 ii=1, 4 

if<min(ii>-2> 10,40.70 

closest point is xr<ii>.yr<ii> 

c at most one intersection 

10 H<rsq_.gt.slsq_(ii)) go to 110 

if(rsq_. lt. srsq_(ii)) go to 110 

d=<sq_rt(rsq_-persq_(ii))-dlp(ii))/dist(ii• isub> 

20 yl=(yl(ii)-yr(ii))*d+yr(ii) 

xl=<xl(ii>-xr(ii))*d+xr(ii) 

intersects in range 

c does circle <radius r> intersect it? 

c store first point 

30 npt=npt+1 

xc<npt)=x1 

yc<npt>=y1 

go to 110 

closest point is at xl<ii),yl(ii) 

40 if(rsq_.lt. slsq_(ii)) go to 110 

H<rsq_. gt. srsq_(ii)) go to 110 

50 d=(sq_rt<rsq_-persq_(ii))+dls(ii))/dist(ii• isub) 

60 x1=xl(ii)+(xr(ii)-xl(ii))*d 

yl=yl(ii)+(yr(ii)-yl(ii>>*d 

go to 30 

70 if(rsq_.lt.persq_(ii)) go to 110 

i f ( r sq_. 1 e. s 1 s 11. ( i i )) g o to 80 

H<rsq_. gt. srsq_<ii)) go to 110 

single intersection 

on side nearer xr,yr 

go to 50 

80 if(rsq_. gt. srsq_(ii)) go to 90 

c two intersections on this side 

can compute angle directly without determining coordinates 

arg=abs(perp(ii))/r 

as=acos<arg> 

go to 100 

90 d=<dls(ii)-sq_rt<rsq_-persq_(ii)))/dist<ii• isub) 

go to 60 

single intersection 

both points are on boundary. calculate angle betwen them. 

100 angle=2.*as +angle 

see if second point only is on boundary 

110 continue 

if <npt> 570, 120, 140 

120 if (angle-0.001) 570,130,130 

following is for case of no single intersection points; 

angle is 2 * pi - angle calculated so far. 

130 pangle=pi2-angle 

go to 560 

140 go to <570,230. 150.500),npt 

c this is an error unless radius is on a vertex and 

c was counted twice. 

150 do 160 i=1.3 

J'"'JSUb(i) 

if ( a b s < x c < i > - x c < J )) . g e. . 0 1) g o t o 160 

if(abs<yc(i)-yc(J)) .le .. 01> go to 210 

160 continue 

c check for one of three points being a vertex 

do 200 i=1.4 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

c 

c 

c 

c 

500 

do 170 J=1, 3 

if <a b s < x c < J ) - x 1 < i > >. 1 e. . 001 > go to 170 

if(abs(yc(J)-yl(i)) . gt .. 001> go to 180 

continue 

go to 200 

if(J-2) 190.190.220 

xc(J)=xc(3) 

yc(J)=yc(3) 

go to 220 

continue 

go to 570 

i f < i . n e. 1 > g o to 220 

xc<2>=xc(3) 

yc<2>=yc<3> 

npt=2 

go to 230 

2 intersection points; determine azimuths. 

if <xc(1)-xnot-r> 260,250.240 

if <xc(1)-xnot-r-0.001) 250.250,570 

azim1=0.0 

go to 300 

if <xc<1>-xnot+r) 270.280.290 

if <xc<1>-xnot+r+0.001) 570.280.280 

a z im1=p i 

go to 320 

azim1=acos<<xc(1)-xnot)/r) 

if (yc<1>-ynot> 310.320,320 

aziml=pi2 - azim1 

if (xc<2>-xnot-r> 350,340.330 

if <xc<2>-xnot-r-O. 001> 340,340.570 

azim2=0.0 

go to 390 

if <xc<2>-xnot+r) 360,370.380 

if <xc<2>-xnot+r+0.001) 570,370,370 

azim2=pi 

go to 410 

azim2=acos<<xc<2>-xnot)/r) 

if <yc(2)-ynot> 400,410,410 

azim2=pi2 -azim2 

pangle=azim2-azim1 

if (pangle> 420,570.460 

if (azim1-azimf) 430.570,440 

pangle=pi2 +pangle -angle 

go'to 560 

if (azimf-azim2) 430.570,450 

pangle=-pangle-angle 

go to 560 

if (azim2-azimf) 470,570.480 

pangle=pi2 -pangle -angle 

go to 560 

if (azimf-azim1) 470.570.490 

pangle•pangle-angle 

go to 560 

four intersection points <each on a different side). 

determine angle by finding closest 2 intersections to 

farthest corner. calculate angle betwen, and add angle 

between other two intersections. 

d ist1=10000000. 

i 1=0 

i2=0 
""' a'l 

VI 

~ 

;;· 

.. 

;;· 

~ 

> 

~ 

I .. 

., 

ci 

Dl 

3 

0' 

.. 

~ ;;· 

3 

;:;· 

:z: 

~ a. 

~ 3 

~ 8" 

i3=0 

i4=0 

do 550 J J=l. 4 

distn=<xl<ifar>-xc(JJ))*(xlCifar>-xcCJJ)) 

& +(ylCifar>-yc(JJ))*CylCifar>-ycCJJ)) 

if Cdistn-dist1> 510.510.520 

510 d ist2=d istl 

dist1=distn 

i4=i3 

i3=i2 

i2=i 1 

i1=JJ 

go to 550 

520 if(distn-dist2) 530.530.540 

530 dist2=distn 

i4=i3 

i3=i2 

i2=JJ 

go to 550 

540 i4=i3 

i3=JJ 

550 continue 

c calculate angle between 2 closest points to farthest corner 

ad=s~rtCCxc(i1>-xc(i2>>*<xcCi1>-xc(i2)) 

& + <yc<i1>-yc<i2>>*<yc<i1l-yc<i2>>> 

pangle=2.*asin(ad/(2.*r>> 

c calculate angle between 2 farthest points from 

c farthest corner and add to previous angle. 

ad=s~rt<<xc<i3>-xc<i4>>*<xcCi3)-xc(i4>> 

& + CycCi3>-yc<i4))*(ycCi3)-ycCi4>>> 

pangle=2.*asin<ad/C2.*r))+pangle 

angle for this radius is now known. calculate risk 

560 return 

c anarea=pangle*r*stsize 

c error printout 

570 write < 16. 580) isub. ifar. npt. xnot. ynot. < x 1 < i), y 1 < i), i•l. 4), r, 

&pangle, <xc<i>. yc<i>• i=1, 4) 

580 format (' ***** error in subroutine inside. source no. '• 

& i3,' debug values follow ...... ',/10x.2i10,10C/10x.2f12.6>> 

590 c a 1 1 e x i t 

end 

C************************************************************ 

subroutine spread(reg, Jlo. ll,totac.dmid) 

magnitude smoother for fault ruptures: 

totac acceleration occurrences calculated for ruptures of 

magnitude fm(ll> at distance dmid from site 

are "spread" to acceleration levels corresponding 

to magnitudes in the range 

fm(lll-delta/2 < fm<ll> < fm<lll+delta/2 

[where fm<ll-ll=fm(lll-deltaJ. 

earth~uake occurrences follow a Gutenberg-Richter 

c magnitude-frequency relationship. calculated from number 

c of earthquakes input at 2 successive magnitudes. 

common/magdis/Jent.mdis 

common/rtb/rtab(20> 

dimension reg<106) 

common/accel/gm<25l.exdif<101),atah(20,25>,aclim<2> 

common/wds/max.maxp 

common/~uad/maxq.ac(105),dumlog.aclog(105) 

10 continue 

if< dmid. e~. -1. > go to 60 

do 20 i=2.mdis 

if(dmid.le.rtabCi)) go to 40 

20 continue 

distance exceeds table--we should not even be here 

write(16,30) dmid 

30 format<' in spread. dmid='e12. 5' too large-should not be here') 

call exit 

40 fr=<rtab(i)-dmidl/CrtabCil-rtab<i-1>> 

ilim=1 

lp=ll+1 

do 50 lm=ll, lp 

aclimCi1im)=atah(i, lm>-fr*<atah(i, lm>-atah<i-1, lm)) 

50 ilim=2 

60 acdif=aclim<2>-aclimC1) 

if(acdif. ne. 0.) go to 90 

Jbot=Jlo+1 

70 if < a c 1 o g < J bot ) . 1 e. a c 1 i m < 1 > > go to 80 

Jbot=Jbot-1 

if< J bot. g e. 1 > go to 70 

80 JUp=Jbot+1 

regCJup>=reg(Jup>+totac 

return 

90 add=O. 

Jbot=Jlo+1 

1 00 if < a c 1 o g < J b o t ) . le . a c 1 i m ( 1 ) ) g o to 11 0 

Jbot=Jbot-1 

i f < J b o t . g e. 1 ) g o to 1 00 

we are alower accel level--for lowest mag Of range 110 aclo=exdif(1l 

JUp=JbOt+1 

120 upl=aclog(Jup> 

if< up 1. g t. a c 1 i m < 2 > > up 1 =a c 1 i m < 2 > 

fr=<upl-aclim(1))/acdif 

iachi=100.*fr+1.00001 

quan=<+exdif(iachil-aclo>*totac 

add=add+quan 

reg(Jup>=reg(JUp)+quan 

if(aclim<2>. e~. upl> go to 130 

aclo=exdif<iachi) 

JUP=Jup+1 

if < J up . 1 e. ma x p ) g o to 120 

write<16, 125> 

125 format<' spread routine needs higher acceleration level;' 

1/' rewrite BOX or increase scale factor') 

call exit 

130 continue 

return 

end 

C*********************************************************** 

subroutine out <reg, naf, xouts. youts. sd) 

common/tims/ftimC20), Jtim<20>.ntims.nwt. iprint 

common/rdis/rdis 

dimension ~p<140),solC20),sols(20> 

common/quad/maxq,ac<105>.dumlog.aclog(105) 

dimension regC106) 

common/ext/prob.ex.exl 

common/wds/max.maxp.maxsp 

dimension regs<106),cdfs(105),cdf(105) 

if( iprint. le. 2> write<16. 10) xouts. youts. rdis. sd 
t 

>(" 

n 

~ 

'I 

c 

c 

c 

10 format(' site at long 'f!8.3', lat 'f!8.3, 

&/'shortest dist to fault='f!9.3' km'/lOx 

&'zero attenuation variability'18x'variability in atten, sigma=' 

& f5.2/2(' g.m. occ/yr exc/yr r<events> r<yrs> ')) 

do 20 1=1,maxsp 

20 regs<l>=O 

30 sd6=4.4*sd 

kls=l 

do 70 1•1. max 

if <reg< 1+1 >. eq. O> go to 70 

asav=<aclog(l)+aclog(l+1))/2. 

kl=kls 

40 if<aclog<kl>-asav.ge.-sd6) go to 50 

kl=kl+1 

if(kl. le. max> go to 40 

go to 70 

:'.50 k ls=kl 

w=<aclag<kl)-asav>lsd 

call bbgaus(w,g1) 

if!(kl. eq. 1> regs<1>=regs(1>+(1. -gl>*reg<l+l) 

60 kl=kl+1 

wh=<aclog<kl>-asav)/sd 

call bbgaus<wh.g2> 

regs<kl)•regs<kl>+<g1-g2)*reg<l+1) 

if<wh. ge.4.4) go to 70 

if<kl.ge.maxsp) go to 70 

g1•g2 

go to 60 

70 continue 

we are really not interested in 1st 4 acceleration boxes but used 

them for computation of variability in acceleration--correct so 

box 5 accumulates all accelerations below ac<S> 

ac(5l is the first value printed 

80 cQrltinue 

regs<1>=regs<l)+reg<1> 

regsav=reg<S> 

regssv=regs(5) 

do 90 i=1,4 

regs<5>=regs(5)+regs(i) 

90 reg<5>=reg(5)+reg(i) 

cdf!s(maxsp >=0. 

cdf <maxp >=0. 

maxq=ma-xp 

k=max 

do 100 i•S.max 

cdf<k>=cdf(k+1)+reg(k+1) 

100 k=k-1 

mm•maxsp-1 

k=mm 

do 110 i=S,mm 

cdfs(k)•cdfs(k+1)+regs<k+l) 

110 k=k-1 

if( <cdfs<S>. ne. 0. 0>. and. (cdf(5). ne. 0. 0)) go to 140 

if<iprint. ge.2> write(16, 120) 

120 format<' no acceleration events accumulated. statistics calc. 

&bypassed') 

do 130 ier=1.ntims 

sol<ier>=O. 0 

130 sols(ier>=O. 

go to 380 

c 

c 

c 

c 

140 continue 

if accelerations do not exceed ac<5), the first acceleration 

of interest, skip statistics calculation to avoid dividing by 

zero or upsetting interpolation routines. 

cnums=cdfs(5)+regs<5> 

cnum=cdf<5>+reg(5) 

yrncnum=l. /cnum 

nrep=O 

idone•O 

if(iprint. eq.3) go to 330 

do 300 i=5,maxsp 

if < i. g t . ma x > g o to 1 50 

if<cdf(i). eq.~O> go to 170 

if<reg(i).ne.O> go to 160 

1 50 i f < r e g s < i ) . e q. 0. > g o to 1 90 

compute return period in number of events and number of years 

return period in events = average number of earthquakes needed to produce an acceleration exceeding a 

r<events>=<yearly earthquakes)/(yearly exceedances of a> 

return period in years = average time between earthquakes 

which cause acceleration a to be exceeded 

r<yrs>•l/(yearly exceedances of a) 

160 if(i. le. max> retev=cnum/cdf(i) 

retyr=retev*yrn cnum 

compute acceleration which has probability prob of not being 

exceeded during t years. this is given by the value of a for which 

prob • exp<-t *number of exceedances of a per year) 

go to 180 

170 retev=99999.9 

retyr=99999.9 

180 if<cdfs(i). eq. 0. > go to 190 

retevs=cnums/cdfs(i) 

retyrs=retevs/cnums 

go to 200 

190 retevs=99999.9 

retyrs=99999.9 

if'(cdf(i). ne. 0.) go to 200 

200 

210 

220 

idone=1 

continue 

qp<nrep+l>=ac(i) 

if(i. le. max) go to 

qp<nrep+2>=0 

qp<nrep+3)=0 

retev=99999.9 

retyr=99999.9 

go to 220 

qp<nrep+2>=reg(i) 

q_p(nrep+3>=cdf(i) 

if(retev. ge. 100000. 

if(retyr. ge. 100000. 

qp<nrep+4>=retev 

q_p(nrep+5)=retyr 

qp<nrep+6>=ac<i> 

210 

retev=99999.9 

retyr=99999.9 
.. 

co 

~ 

tD 

{ii" 

.. 

~ 

> 

f 

;-

.. 

"a 

ci 

~ 0' 

.. 

~ {ii" 

3 

;:;· 

:c 

E 

a. 

~ 

l i" 

l} 

c 

0 

G> 

~ m 

JJ 

z 

3: 

m 

~ 

-o 

JJ 

~ z 

G> 

0 

,.. 

,.. 

0 

!1:1 

~ 

I 

~ ~ 

~ JJ 

!B 

0 

z 

~ CXI 

230 

240 

250 

260 

270 

280 

300 

310 

320 

330 

340 

qp<nrep+7>=regs<i> 

qp<nrep+B>=cdfs(i) 

if<retevs. ge. 100000. 

if< ret11rs. ge. 100000. 

qp<nrep+9l=retevs 

qp<nrep+10l=ret\jrs 

retevs=99999.9 

ret\jrs=99999.9 

if< nrep. ne. 130) go to 260 

if<naf. le. 3> go to 230 

write<16.280) qp 

go to 240 

write <16.270> qp 

continue 

nrep=O 

format< 1 h 1) 

if(idone. eq. 1) go to 330 

go to 300 

nrep=nrep+10 

format ( 1 h f8. 2, 2f9. 5, 2f8. 1. f12. 2, 2f9. 5, 2f8. 1) 

format< 1h fB. 4, 2f9. 5, 2f8. 1. f12. 4, 2f9. 5, 2f8. 1 > 

if<retev. eq.99999.9 .and. retevs. eq.99999.9> go 

continue 

if<nrep. eq.O) go to 330 

if(naf. le. 3) go to 320 

write<16,280) <qp(ll), 11=1,nrep> 

go to 330 

write(16,270) (qp(ll), ll=Lnrep> 

if <iprint.gt.2> go to 350 

write<16,340> cnum.cnums, sd 

format(2(' total \jearl\1 events 'flO. 5, 12x> 

&//10x'zero attenuation variabilit\1'18x, 

& 'variabilit\1 in atten, sigma='f4.2> 

to 310 

compute acceleration which has probabilit\1 prob of not being 

exceeded during t 11ears. this is given b\1 the value of a 

for which 

prob = exp<-t * number of exceedances of a per \jear> 

350 reg<5>=regsav 

regs<5>=regssv 

do 360 JV=l.ntims 

maxq=max 

call 1 inint <reg, ftim( JV), sol< JV) > 

maxq=maxsp 

call linint(regs,ftim(Jv),sol~<Jv>> 

if(iprint. le. 2> write<16,370>prob, sol<Jv), Jtim<Jv), sols(Jv), 

&Jtim< JV) 

360 continue 

370 format<f10.3' ext prob ='f7.3' for 'i5' \jears'19x, 

& f7.3' for 'i4' 11ears'> 

380 rat=O.O 

if<sol(l>. ne. 0. O> rat=sol<ntimsl/sol<l> 

ratsd=O 

if<sols<l>.ne.O> ratsd=sols<ntimsl/sols(l) 

if< iprint. le. 2>wri te < 16, 390) Jtim(ntims l, prob, Jtim< 1 >,rat, ratsd 

if< iprint. ge. 2> write <2> xouts. \lOUts, <sol< Jt >, Jt=l. ntims l, 

&<sols<Jt), Jt=l.ntims> 

390 format<' ratio'i5' \lr'f6.3' extreme value to'i5' 11r val 

&='f7.2,7x,f7.2> 

return 

end