Finite Element Glacier Dynamics Model 

Applied to Columbia Glacier, Alaska 

By W. G. SIKONIA 

STUDIES OF COLUMBIA GLACIER, ALASKA 

GEOLOGICAL SURVEY PROFESSIONAL PAPER 1258-B 

Description of a computer model indicating 

that Columbia Glacier will begin 

a rapid, catastrophic retreat in 1983 

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON:1982 
UNITED STATES DEPARTMENT OF THE INTERIOR 

JAMES G. WATT, Secretary 

GEOLOGICAL SURVEY 

Dallas L. Peck, Director 

Library of Congress catalog-card No. 81-607563 

For sale by the Superintendent of Documents, U.S. Government Printing Office 

Washington, D.C. 20402 
CONTENTS 

Page 

Symbols and abbreviations------------------------------------------------------- BIV 

Abstract --------------------------------------------------------------------- 1 

Introduction ------------------------------------------------------------------ 1 

The finite element model -------------------------------------------------------- 1 

Calving law------------------------------------------------------------------- 3 

Modeling Columbia Glacier ------------------------------------------------------ 5 

References cited --------------------------------------------------------------- 28 

ILLUSTRATIONS 

Page 

FIGURE 1. Index map of Columbia Glacier, Alaska----------------------------------------------------------------------- B2 

2. Graph showing calving law --------------------------------------------------------------------------------- 3 3. Contour map showing flow line used in modeling ___________________ .:____________________________________________ 7 

4. Finite element grid --------------------------------------------------------------------------------------- 8 

5. Flow-line profiles showing projected Columbia Glacier behavior--------------------------------------------------- 9 

6-9. Graphs showing: 

6. Terminus position in time ---------------------------------------------------------------------------- 23 

7. Unsupported height at the terminus, in time------------------------------------------------------------- 24 

8. Surface horizontal speed at flow line ------------------------------------------------------------------- 26 

9. Rate of iceberg production--------------------------------------------------------------------------- 27 

TABLES 

Page 

TABLE 1. Calving law data----------------------------------------------------------------------------------------- B4 

2. Flow-line data------------------------------------------------------------------------------------------- 6 

APPENDIXES 

Page 

APPENDIX A. Mathematics ------------------------------------------------------------------------------------------ B29 

Tangent stiffness matrix ----------------------------------------------------------------------------- 29 

Upper free-surface dynamics------------------------------------------------------------------------- 30 

Surface tractions and body forces --------------------------------------------------------------------- 31 

Updating the mesh --------------------------------------------------------------------------------- 31 

Future modification--------------------------------------------------------------------------------- 31 

B. User instructions for the computer program ---------------------------------------------------------------- 32 

Coordinate data ------------------------------------------------------------------------------------ 32 

Material property data ------------------------------------------------------------------------------ 32 

Boundary restraint data ----------------------------------------------------------------------------- 32 

Nodal forced boundary data -------------------------------------------------------------------------- 32 

Nodal temperature data----------------------------------------------------------------------------- 32 

TANF macroinstruction----------------------------------------------------------------------------- 33 

REST and SAVE macroinstructions ------------------------------------------------------------------- 33 

UPNO macroinstruction -----------------------------------------------------'------------------------- 33 

PLOT macroinstruction ----------------------------------------------------------------------------- 33 

C. New subroutines -------------------------------------------------------------------------------------- 34 

D. Modifications of the general purpose finite element program--------------------------------------------------- 47 

E. Sarp.ple run ------------------------------------------------------------------------------------------- 54 

III 
IV CONTENTS 

SYMBOLS AND ABBREVIATIONS 

A ----- Ice flow Jaw parameter a ___ Year 

a ___ Calving law paramet:er 

a" __ Rate of ice-thickness-change vector 

aP ____ Nodal pressure vector 

a" ___ Nodal velocity vector 

B ----- Constant in Budd sliding law 

B ----- Deformation rate shape function 

b ----- Calving law parameter 

b ----- Annual ice balance 

c ----- Calving law parameter 

D ----- Discharge 

D' ---- Constitutive matrix 

DOF __ Degree of freedom 

d ___ Day 

d ---- Effective shear deformation rate 

d41 -- Deformation rate tensor 

d ----- Deformation rate vector 

6 ----- Indicates variation of following variable 

a ----- used in partial derivative 

E ----- •(11~&) D' 

F(a) -- Nonlinear ice-tlow function of nodal variables a 

f ----- Force vector 

G ----- Auxiliary matrix used in tangent matrix derivation 

'Y ----- Angle of surface 

. r. ---- Ice surface 

h ---- Total ice thickness; related quantities have this as 

superscript 

h, ---- Altitude of bed 

h, ----- Altitude of surface 

h.. -- Height of the ice column unsupported by water pressure 

h., __ Water depth 

i __ Coefficiellt matrix related •to spatial variation of flow field 

• 1 .It"' -- Coefficient matrix related to ;rate of ice thickness change 

1[, K• _ Coefficient matrix describing effect of ice viscosity 

KP ___ .Coefficient amtrix related to pressure distribution 

km __ Kilometer 

Mg --- Megagram M ____ Mass matrix 

m ___ Meter 

" ----- Viscosity N ___ Shape function matrix 

n ____ Ice flow law parameter 

0 ----- Volume 

p ----- Pressure; related quantities have this as superscript 

r ____ Coefficient of determination 

p, ---- Density of ice 

p. --- Density of water s ____ Second 

s ---- Matrix used in determining emergence velocity a41' ____ Deviatoric stress tensor 

t1 ----- Deviatoric stress vector 

T ----- Denotes transpose matrix as superscript 

.,, ----- Shear stress at the bed 

UTM __ Universal Transverse Mercator 

u ___ Denotes velocity; related quantities have this as 

superscript 

u ----- Velocity vector 

V, ___ Sliding speed at the bed 

Vc ---- Calving speed 

V. ____ Emergence velocity 

V. __ Surface ice speed 

X(t) ___ Position of calving terminus at timet 
STUDIES OF COLUMBIA GLACIER, ALASKA 

FINITE ELEMENT GLACIER DYNAMICS MODEL 

APPLIED TO COLUMBIA GLACIER, ALASKA 

By w. G. SIKONIA 

ABSTRACT 

A two-dimensional finite element computer model describes the 

dynamics of glacier flow and includes the flow of ice, sliding, calving, 

and surface balance. Surface profile and terminus position are followed 

in time. A calving law that relates calving speed to water discharge at 

the bed of the glacier and to the height of the ice column unsupported 

by water pressure is shown to fit observed calving at Columbia Glacier. 

Application of the flow model to Columbia Glacier indicates that catastrophic, rapid retreat from the terminal moraine will begin in 1983. 

This will produce an increase in the production of icebergs that can 

drift into oil tanker lanes to and from Valdez, Alaska. 

INTRODUCTION 

Columbia Glacier, Alaska (fig. 1), is in a state precariously close to catastrophic retreat (Post, 1975; Meier 

and others, 1978; Sikonia and Post, 1980; Meier and 

others, 1980a, b). A dynamic model of glacier behavior 

was sought that could adequately represent at least the 

lowest 14 km of the glacier; the model must be capable 

of dealing with a large amount of sliding, with flow up a. 

terminal moraine, and with a moving calving face. Bindschadler applied a one-dimensional model (Bindschadler, 

1978; Meier and others, 1980b), but it seems questionable whether the approximations of the model adequately represent reality near the calving face, especially with regard to surface-slope averaging schemes and 

longitudinal stresses. 

THE FINITE ELEMENT MODEL 

The two-dimensional finite element model described in 

this paper is an attempt to relate more adequately the 

motion of the ice mass to the physical processes involved. However, the two-dimensional nature of the 

model dictates that it also is an idealization of reality. 

The model represents a vertical cross section through a 

glacier along a central flow line, but no attempt has been 

made yet to include effects of lateral shear stresses or 

lateral deformation rates. 

The ice flow law describing the deformational contribution to ice motion is 

(1) 

where du is the deformation rate tensor; o' u. the 

deviatoric stress tensor; r, the effective shear stress 

defined by 

(2) 

(summation convention on subscripts i,J); and A and n 

are constants, with approximate values of 

A=0.140bar..,.a·•andn=3. (3) 

Sliding is modeled as a layer of very soft ice (large 

valued A) along the bed. In the Columbia Glacier application, the layer is assumed to be a constant 10 m thick, 

with n=3 and A ranging from 7 to 40 bar-"a-•, fitted 

along the bed so that the correct sliding is obtained to 

match measured surface velocities. This is approximately equivalent to use of a sliding law suggested by Budd 

and others (1979) in the form 

(4) 

Here Vb is the sliding speed on the bed; B, a constant 

perhaps approximately 1,800 m1

a-•bar-" for n=3; rb, the 

shear stress at the bed; and h.., the ice-column height unsupported by water pressure 

h..=h-h.,pjp, (5) 

where his total ice depth; h.,, water depth; and p., p1, the 

densities of water and ice, respectively. Bindschadler 

(written commun., 1978) found that this sliding law fits 

Variegated Glacier, Alaska, data. 

The solution proceeds in a sequence of time steps. At 

each step, equilibrium iterations using Newton's method 

are done first to produce consistent stress-velocity fields 

described by the nonlinear equation 1. By using the 

velocity field, the surface and calving front then are projected to new positions, and the process stepped in time. 

Bl 
B2 STUDIES OF COLUMBIA GLACIER, ALASKA 

0 5 10 KILOMETERS I 1 e I I I I I I I I 

FIGURE 1.-Index map of Columbia Glacier, Alaska. Arrows show direction of flow. Main ice stream is indicated by longer arrows and dots at 

2-km intervals along the longitudinal coordinate system. 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B3 

The finite element framework for the computer program is a general purpose program (Taylor, 1980; 

Zienkiewicz, 1977, p. 733). Velocity and pressure are 

chosen as the variables, and the treatment follows that 

of Zienkiewicz (1977, p. 612) except that, for ice, the 

viscosity, p., is not constant but is rather a function of the 

effective shear stress, 

p. = 1/(2A T"-•). (6) 

to water depth has been suggested (Meier and others, 

1980b). The suggested water-depth law was restricted 

to yearly averaged data that was averaged also across 

the width of the glacier. It definitely does not apply for 

intervals of less than a year for calving at a single location on the terminus. When applied to the calving of Columbia Glacier at the head of the embayment, by using 

photogrammetrically determined terminus positions 

taken at approximately 2-month intervals from July 24, 

1976, to September 2, 1980, the water-depth law yields a 

coefficient of determination, r, of only 0.08. However, 

for Columbia Glacier, calving speed appears to be 

related to water discharge at the bed (Sikonia and Post, 

1980). A law that does appear to fit the observed intrayearly calving at Columbia Glacier is (see also fig. 2 and 

table 1) 

(8) 

Sr-------------------------------------------------------1 

Also, kinetic effects can be ignored because of the slow 

creep of ice. Because the viscosity is nonconstant, the 

equations are nonlinear, and the tangent stiffness 

matrix for use in Newton-Raphson iteration takes the 

appropriate derivatives of viscosity into account. Other 

finite element programs with a similar approach are attributable to Schmidt (1977), Raymond (1978), and 

Hooke and others (1980), but these treat only geometrically static, purely creeping flow. The program 

described in this paper treats deformational motion, 

sliding, balance, and calving effects in a geometrically 

dynamic setting. It also starts with a good, readily 

available general purpose finite element framework; the 

advantage of this approach is that only the specific program modules directly concerned with the glacier 

dynamics need be considered, greatly simplifying pro- a: graming, documentation, and orientation of new users. ~ 

To take into account the dynamics, the movement of > 

the upper free surface must be followed. This surface is ~ 

followed essentially as in Zienkiewicz (1977, p. 623); the 

free surface must be projected along the surface velocity 

vectors to its new position. The program makes this projection by defining an auxiliary variable-the vertical 

surface speed- decomposing the surface velocity into :.:: 

general nonorthogonal components along the surface z 

and in the vertical direction. In addition to the move- g 

ment of the upper free surface, the movement of the ~ 

calving face is projected from one time step to the next en 

by using the calving law described in the section on 

"Calving Law." Because of large changes in the 

configuration of the ice mass, the entire mesh must be 

reconstructed at each time step. 

CALVING LAW 

The position of the calving terminus, X(t) (assumed 

vertical), in time is determined by the imbalance between the horizontal surface ice speed, V., at the terminus (positive downglacier) and the calving speed, Vc 

(positive upglacier), 

(7) 

Alternatively, equation 7 defines calving speed, Vc, as 

the rate of retreat of the calving terminus in the moving 

ice stream. A calving law relating calving speed linearly 

Q 

w 

> 

a: 

w 

en 

In 

0 

COMPUTED CALVING SPEED, IN KILOMETERS PER YEAR 

FIGURE 2.-Calving law V.=aD"h,;, where a=1.092xl03 (km/a) 

(m3/s)-"m-c, b=0.5689, and c= -2.175. V. is speed of calving at the 

head of the embayment, D is the discharge of the Knik River, and h.. 

is the height of the unsupported ice column. The data points, as given 

in table 1, were determined from analysis of aerial photographs. 

Point 8 corresponds to a time interval during which Kadin Lake suddenly released a large volume of water into the subglacial flow; in 

this interval, the Knik River discharge would not correlate well with 

this flow. Points 8, 23, and 24 were not used in the regression that 

fitted the constants a, b, and c. 
B4 STUDIES OF COLUMBIA GLACIER, ALASKA 

TABLE 1.-Calving law data 

[N, data interval number. t1 and ~. beginning and ending dates of interval. flt, number of days in the interval. h.~ and h.a, ice-cliff heights at head of embayment at 

times t1 and~. respectively. hw1 and hwa, water depths at head of embayment at times t1 and~. respectively. D, average discharge ofKnik River during interval. h.., 

average unsupported height during interval. (V.). and (V.)., observed and calculated calving speeds at head of embayment. fl V., difference between calculated and 

observed calving speed. Calculated speed is based on the calving law V.= aD'h.:] 

N tl ~ flt h,l h,a hwt 

(d) (m) (m) (m) 

1 ---- 07/24176 10/01176 69 98 86 127 

2 ---- 10/01176 11/17176 47 86 86 180 

3--- 11/17176 01/19177 63 86 88 163 

4_ 01119/77 03/07177 47 85 89 155 

5 ---- 03/07177 04123177 47 88 94 137 

6 ---- 04/23/77 06/02177 40 94 100 100 

7 ---- 06/02177 07/07177 35 99 95 80 

8 --- 07/07/77 08/29177 53 93 99 102 

9 --- 08/29177 11/08177 71 100 78 200 

10 ---- 11/08/77 02/28178 112 81 97 170 

11--- 02128177 04/19178 50 96 84 125 

12--- 04/19178 06/11178 53 68 81 110 

13 ---- 06/11178 07/30/78 49 80 94 112 

14 ---- 07/30/78 08/26178 27 90 90 160 

15 ---- 08/26/78 11/08178 74 92 77 185 

16 ---- 11108178 01/06/79 59 72 76 215 

17 ---- 01/06179 04/12/79 96 72 77 182 

18 ---- 04/12/79 08/18179 128 93 111 85 

19 ---- 08/18179 10/20/79 63 70 76 82 

,20 ---- 10/20179 02/29/80 132 77 75 195 

21--- 02/29/80 06112180 73 74 90 167 

22 ---- 05/12/80 07/22/80 71 80 101 110 

23 ---- 07/22/80 09/02/80 42 92 89 90 

24 ---- 09/02/80 10/30/80 58 91 85 110 

Vc is the calving speed at the head of the embayment, 

and hv is the ice column height unsupported by water 

pressure, as defined in equation 5. The variable D would 

most appropriately be the actual discharge at the bed of 

the Columbia Glacier, but because actual discharge 

would be almost impossible to measure, a related 

discharge, that of nearby Knik River, was used instead 

(see Sikonia and Post, 1980). Constants a, b, and c have 

fitted values as follows: 

a= 1.092 x 106 (m/a)(m3/s)-"m-c, 

b=0.5689, 

C= -2.175, 

p ... = 1.018 Mg/m3 , 

p~=0.900 Mg/m3 • 

(9) 

(10) 

hwa D h.. (V.). (V.). flV. 

(m) (m3

/s) (m) (m/a) (m/a) (m/a) 

180 464.65 71.87 3,964 3,293 -671 

175 119.82 62.73 2,008 2,047 39 

146 59.96 66.74 1,190 1,207 17 

142 63.45 67.53 1,697 1,215 -482 

100 30.68 75.46 307 631 324 

80 110~07 85.20 1,259 1,002 -257 

109 486.96 84.61 2,256 2,371 115 

278 765.94 71.09 7,059 4,481 -2,578 

200 231.72 62.78 2,738 2,974 236 

125 21.12 69.66 1,175 607 -568 

112 30.90 74.46 712 652 -60 

113 126.76 59.88 1,832 2,339 509 

177 460.50 68.05 4,362 3,689 -673 

182 684.48 67.58 4,588 4,692 104 

210 230.27 58.61 3,150 3,442 292 

182 22.46 47.97 1,657 1,415 -242 

140 17.36 53.39 1,004 969 -35 

135 412.75 87.58 2,318 2,002 -316 

195 350.34 54.84 6,669 5,049 -1,620 

167 44.71 52.27 1,478 1,737 259 

130 35.61 62.53 683 1,034 351 

137 390.94 74.31 1,748 2,775 1,027 

118 587.54 76.86 2,880 3,252 372 

133 207.33 72.07 2,205 2,068 -137 

The standard error for the Columbia data is 529 m/a (in 

calving speeds ranging from 307 to 4,588 m/a), with a 

coefficient of determination of the logarithmically derived corresponding linear regression problem of 

r = 0.83. (Although 1.018 Mg/m3 is the measured density 

of seawater in front of Columbia Glacier, the modeling 

computer runs (see section entitled "Modeling Columbia 

Glacier') actually were made by using p ... =1.000 Mg/m3 

and p~=0.900 Mg/m3 • In this case, a=1.798x106 

(m/a)(m3/s)-"m-c, b=0.5510, C= -2.239, rl=0.83, and the 

standard error was 369 m/a, based on fewer data points 

than the constants in equation 9 because of limited data 

availability at the time. It was not considered worth rerunning the program for these minor changes in constants.) 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B5 

It is uncertain exactly what the relation between the 

discharge-unsupported height calving law and the 

water-depth law is, but there appear to be at least three 

possibilities: 

1. One or the other law is wrong. 

2. Some, as yet, not understood process averages the 

above law into the water-depth law for time periods 

of a year or longer. 

3. The two laws are describing two, as yet, undistinguished distinct calving processes. 

The discharge-unsupported height law is used in this 

study for reasons discussed in this section. This law also 

accommodates the observation that embayments form 

at the location of the main water discharge from the terminus and supports the physical idea that the closer the 

terminus is to floating, the more easily this discharge 

can bring about calving and ice removal from the face. 

MODELING COLUMBIA GLACIER 

The finite element model was applied to the lowest 14 

km of Columbia Glacier to gain a better understanding 

of its possible future behavior (see also figs. 3 and 4 and 

table 2). The constant a in the calving law is too large for 

an idealized two-dimensional glacier; it produces an immediate catastrophic retreat, from the starting modeling date 1978.150, that has not occurred on the real 

glacier. This is understandable because embayments 

and related surface lowering are filled partly by flow 

from the sides, a factor which the two-dimensional 

model does not take into account. Also there is an "archbridging'' stress effect within the ice around an embayment that could partly support the ice cliff at the head of 

the embayment, thus providing a larger unsupported 

height, h.,, and a smaller calving speed, Vc. This effect 

also is not taken into account in the two-dimensional 

model. An a corresponding to width-averaged calving is 

needed. Rather than trying to derive the averaging 

effect, the exponents b and c were fixed, and the a 

parameter adjusted for the two-dimensional glacier to 

produce future behavior that seems reasonable in view 

of the actual configuration of the glacier from early 1978 

to the present. 

Future behavior is extremely sensitive to changes in 

the a parameter. For example, a= 1168750 produced a 

collapse that could not be followed with a time step of 

0.025 year by 1982.075; the terminus had retreated 

2,285 m by 1982.025. Decreasing a by only 0.27 percent 

to a= 1165625 extended extremely rapid collapse to 

1984.600, the terminus retreating 3,031 m by 1984.525. 

Decreasing a by another 0.27 percent to a= 1162500 

produced a gradual advance to the crest of the terminal 

moraine by 1982.400 (see also figs. 5, 6, 7, 8, and 9). The 

modeling also showed that a collapse is preceded by thinning of the glacier to the point that the unsupported 

height factor, h.,, of the calving law is small enough to 

produce a large calving speed when applied in equation 

8; for example, h .. = 7 m at the terminus, time 1982.250 

and a= 1168750. 

The sensitivity of the glacier's behavior to the a 

parameter may be overstated or at least compressed in 

time scale because of the two-dimensional idealization, 

but, nevertheless, the behavior of the real glacier is 

probably still quite dependent on relatively small 

changes in the calving. As indicated by the calving law, 

for the real glacier, decreased calving that would stall 

collapse could be brought about by low runoff or by 

outlet of the major discharge at a location having a 

relatively high ice cliff and shallow water depth. Reference to figure 7 shows that h,. was relatively large in 

1980, and, as indicated on figure 6, the retreat was, indeed, slowed that year. However, heavy rains in July 

and August 1981, together with a partial outburst of 

Kadin Lake, produced a large discharge that caused a 

large embayment to form. The data presented in figure 

6 show that the observed terminus position is following 

again the retreat predicted by the model using 

a= 1165625 and that the indicated collapse of Columbia 

Glacier is imminent. 
B6 STUDIES OF COLUMBIA GLACIER, ALASKA 

TABLE 2. -Flow-line data 

[Col., column number. Easting and Northing, UTM coordinates. h,, ice surface altitude. h., bed altitude. b, annual balance. V., surface speed. h, annual surface altitude change] 

Col. Easting Northing h. h. 6 v. 1i 

(m) (m) (m) (m) (m/a) (m/a) (m/a) 

1.00 ----------- 496047.77 6776061.73 465 -278 -2.67 950 -2.6 

2.00 ----------- 496190.00 6775312.92 430 -247 -3.42 1,030 -3.8 

3.00 ----------- 496212.11 6774551.05 411 -262 -3.98 960 -4.0 

4.00 ----------- 496069.20 6773802.37 401 -238 -4.45 922 -4.7 

5.00 ----------- 495989.61 6773044.34 370 -204 -4.90 935 -3.2 

6.00 ----------- 496058.28 6772285.25 340 -191 -5.41 930 -2.8 

7.00 ----------- 496296.14 6771561.12 314 -224 -5.98 900 -1.9 

8.00 ----------- 496606.30 6770864.89 290 -320 -6.30 839 -6.3 

9.00 ----------- 496867.26 6770148.76 275 -382 -6.32 808 -8.6 

10.00 ----------- 497094.55 6769421.24 264 -395 -6.29 812 -6.0 

11.00 ----------- 497283.81 6768682.92 248 -384 -6.26 866 -6.9 

12.00 ----------- 497449.85 6767939.03 226 -383 -6.27 904 -7.3 

13.00 ----------- 497519.90 6767180.06 210 -402 -6.36 933 -6.9 

14.00 ----------- 497522.32 6766417.87 198 -401 -6.53 1,012 -8.5 

15.00 ----------- 497459.13 6765657.99 178 -402 -6.68 1,102 -9.6 

15.50 ----------- 497402.49 6765281.07 167 -402 -6.76 1,150 -9.5 

16.00 ----------- 497340.59 6764905.06 159 -382 -6.91 1,210 -9.0 

16.50 ----------- 497273.44 6764529.95 152 -355 -7.09 1,320 -8.4 

17.00 ----------- 497201.03 6764155.75 140 -320 -7.21 1,460 -8.2 

17.50 ----------- 497115.89 6763784.20 129 -286 - 7.38 ------------------ -7.3 

18.00 ----------- 497024.00 6763414.39 118 -260 - 7.63 ------------------ -3.0 

18.25 ----------- 496975.53 6763230.14 113 -242 - 7. 77 ------------------ -0.7 

18.50 ----------- 496925.37 6763046.32 108 -221 - 7.87 ------------------ 1.0 

18.75 ----------- 496873.53 6762862.94 100 -196 - 7. 93 ------------------ 1.1 

19.00 ----------- 496820.00 6762680.00 96 -162 -7.94 ------------------ 0.1 

19.25 ----------- 496759.80 6762499.20 96 -135 - 7. 94 ------------------ -1.1 

19.50 ----------- 496699.05 6762318.60 96 -106 - 7. 94 ------------------ -2.5 

19.75 ----------- 496637.76 6762138.18 96 - 60 - 7.94 -------------------------- 20.00 ----------- 496575.92 6761957.94 96 - 10 - 7. 94 -------------------------- 20.25 ----------- 496513.54 6761777.91 ---------------- - 68 ---------------------------------------- 20.50 ----------- 496450.61 6761598.05 ---------------- -132 ---------------------------------------- 20.75 ----------- 496387.14 6761418.37 ---------------- -162 ---------------------------------------- 21.00 ----------- 496323.12 6761238.89 ---------------- -168 ---------------------------------------- 21.25 ----------- 496257.12 6761060.15 ---------------- --178 ---------------------------------------- 21.50 ----------- 496190.01 6760881.82 ---------------- -182 ---------------------------------------- 21.75 ----------- 496i21.77 6760703.90 ---------------- -193 ---------------------------------------- 22.00 ----------- 496052.42 6760526.38 ---------------- -202 ----------------------------------------
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B7 

+ 

0 

I 

~---­

( 

I 

......... ............. '\ 

,..J 

1 

I 

2 

I 

( 

) 

3 

I 

4 

I 

5 Kl LOMETERS 

I 

+ 

+ 11" 

147° 

• 67IOOoom.N. 

50000011\f. 

FIGURE 3.-Flow line used in the modeling. The 

numbers along the flow line are the column 

numbers, and UTM locations and data correspond to those given in table 2. Bed contours 

are solid, surface contours are dashed, and 

the terminal moraine is dotted (column locations 21 and 22 are in Columbia Bay). Kadin 

Lake is at UTM coordinates 490000 m.E., 

677 4000 m.N. 
B8 

ELEMENT 1 

ELEMENT BOUNDARY 

.. • 

- ....... 

~ 

I ----- ... • 7 I I 8 

9 

13 14 

STUDIES OF COLUMBIA GLACIER, ALASKA 

ELEMENT BOUNDARY 

ELEMENT 4 

ELEMENT 6 

(10m THICK) 

• 

I 

10 

5 

ELEMENT 7 

• 

• 

6 

DISTANCE SCALE 0 0.5 1 km 

I I I I I I I I I I I 0 10 20 kmla 

VElOCITY VECTOR SCALE 

- ... ... I I 

11 12 

~ ... 

... .... !" 

I I 

13 

FIGURE 4.-Finite element grid. There are three elements through the thickness, with the bottom element a uniform 10m in thickness. The ice 

velocity vectors are scaled to represent 0.05 year of displacement (if the ice continued beyond the velocity sampling points in straight lines, 

which, in general, it does not do). 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B9 

DISTANCE SCALE 

0 0.5 1 km 

11111111111 

0 10 20km/a 

VELOCITY VECTOR SCALE 

; ; ; ; ; ; ; ~ ~ ~ t~~Fo 11 1 9 7 8 . 1 50 14 15 16 17 

I 1 

1978.275 

I 

~ I ~ ; I ~ 1978.525 ' 

i ; ; 

I 

; 

I 

1978.650 

; 

I 

; ; 

I I 

1978.775 

I 1 I 1 I I 

FIGURE 5.-FuL~ ll~bia~lgig behavior, projected by using modeling constant a= 1165625. The individual illustrations show 

the flow-line profile from column locations 11 through 20 as a function of time. 
BlO STUDIES OF COLUMBIA GLACIER, ALASKA 

:~ 11 

I 

; ~ ·; 

I I 

1979.025 14 

I 

; ; ; 

I I 

1979.150 

: ~ I 

; 

I ~ ~ . I I I 

1979.275 

I 

; 

I 

; ; I I I 

1979.400 

I 

; 

I 

; 

I 

~ I I I I 

1979.525 

I I I 

1979.649 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA Bll 

; 

I 

11 

: ; : !' r , 

1980.024 

I I 1 I 

1980.149 

I I 

1980.274 

~ ~ . . I I I 

1980.399 

I 

~ I ~ I 

~ I 

1980.524 

; : ; I I I 

1980.649 

; 

I 

; ; I I 

1980.774 

; ; ; 

I I I 

1980.899 

: : : ; ~; ~ ;_,_.,... ; ; ~~ r 1 

r r I 1 17 15 16 

I • 

. . . 

; 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using modeling constant a=1165625. The individual 

illustrations show the flow-line profile from column locations 11 through 20 as a function of time. 
B12 STUDIES OF COLUMBIA GLACIER, ALASKA 

; ; ; ; ; ; ~ ~ ~ ~ ~~~~0 I 17 18 19 

11 1981 .. 024 14 15 16 

I 

~ I 

; 

I 

; I 

; I 

; I 

; 

r 

; t 

; I 

; I 

; 111 ==~ 1981.150 

; ; ~ ; ; ; ; ; ; ; ;:~~ I I I I 

I I I I I I I I I 

1981.275 

~ ; ~ ; ~ ; ; ; ; ; ;;~ I I I I 

I I I I I I I I I 

1981.400 

~ ~ ; ; ; ; ; ; ; ; ; ; ~q "(1"1 

I I I I . I I I I I I I I I 

1 981 . 525 

; 

I 

; ; 

I I 

1981.650 

; 

I : ; 

I I 

I~ 

1981.775 

; 

I 

; ; I I I I I I 

1981.900 

FIGURE 5 (Continued). -Future Columbia glacier behavior, projected by using modeling constant a= 1165625. The individual 

illustrations show the flow-line profile from column locations 11 through 20 as a function of time. 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA 

; ; ; ; ; ; ~ ; ~ ~ 18 ~~~0 19 

',11 1982.025 14 15 16 17 

; ~ ; ; ; ; ; ; I 

; I 

; I =~ I 

I I I I I I I I 

1 98 2. 1 5.0 

I ~ I 

; 

I ~ I 

; ~ I' 

; I 

; I 

; I 

; I 

; I ~~ 1982.275 

~ ; ·~ ~ .. ~ ; ; ; ; I 

>; ,· .• ~··~~ I I I I I I I I 

1 9-82. 400 

~ ~ ~ ~ ~ ; ; ; ; I 

; ,:,;,~ I I I I I I I I I 

1982.525 

; 

• ·• 

: ; 

• ; • ; I ; I ;;:~ I I I 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using,modeling constant a=1165625. The individual 

illustrations show the flow-line profile from column locations 11 through 20 as a function of time. 

B13 
B14 STUDIES OF COLUMBIA GLACIER, ALASKA 

; • • . • • T r 

1983.149 

I 

; ~ ; I I 

1983.274 

I ~ I 

~ I 

~ 1983.399 

r---~~----~1 --~~~----~--~--ri-JI 

.. 

• • • 

1983.524 

I I 

1983.649 

I I 

1983.774 

- : :: I I • 

1983.899 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using modeling constant a=1165625. The individual 

illustrations show the flow-line profile from column locations 11 through 20 as a function of time. 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA Bl5 

I 1 I I I 

1984.150 

I I 

1984.275 

• 

-

- 1984.650 

r 

I I 

1984.775 

I I 

1984.900 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using modeling constant a= 1165625. The individual 

illustrations show the flow-line profile from column locations 10 through 18 as a function of time. 
B16 STUDIES OF COLUMBIA GLACIER, ALASKA 

- - ------

9 10 

1985.775 ' ----~--~------, 

FIGURE 5 (Continued).- Future Columbia glacier behavior, projected by using modeling constant a= 1165625 . .TI:ie indiVidual 

illustrations show the flow-line profile from column locations 1 through· 10 as a function of time.' 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA Bl7 

9 10 

1 986 · 1 50 II --r------. 

1 986 · 27 5 Tl -----r-__, 

~;~~;;;;;~~:~;;;;~; qw--

1 986 · 400 II -~-----. 

----,, ---.----~ 

I 1986.900 -----11 ------.------. 

FIGURE 5 (Continued).- Future Columbia glacier behavior, projected by using modeling constant a= 1165625. The individual 

illustrations show the flow-line profile from column locations 1 through 10 as a function of time. 
B18 STUDIES OF COLUMBIA GLACIER, ALASKA 

I I I 

1 2 3 4 1 987. 025 7 -----,1 -..._,.._ _ ___, 

8 91 I 10 

~ ~ ; ~ ; ; ; ; ; ; 1; ; ; ~ :; ; ; =i ; ; :1: 1 

I 1987.150 ---- 11 ---.---____, 

I 1 9 8 7 · 2 7 5 ---II -..._,.._-____, 

~ ; ~ ~ ; ; ; ; ; ; ~;; ~ :;; ; ; :rF: 

I 1987.400 

1987.525 ---II ---,....------. 

~;; ;1;;;~;;;~;;:;1;;;; 1 q; 

1987.650 

I---.. II ----r------. 

;:....i.J; ;~;~;~: :~~ ;~; ;~~ ;~;; ;~;; ; ;...,; ;li~------.-1 ----- 1 9 8 7 • 7 7 5 II ----.--.,.1 ·---. 

;;;~I;;=~;;;~;;;~;:;;;:;= I1 987 · 900 I I I 

FIGURE 5 (Continued).- Future Columbia glacier behavior, projected by using modeling constant a= 1165625. The individual illustrations show the flow~line profile from column locations 1 through 10 as a function of time. 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B19 

1 2 3 4 1 988.025 ·7 k-~9,------., 10 

LL~WLL~jU~1J~,~~~--~--

1 9 8 8 · 1 50 II --,..--_ ____, 

~ ; ; ;, ; ; ; ~ ; ; ~ ; ; ; i;;,; : : ; ;m= 

1988.275 ---·~-~--

;L-L i i~~. ;l_ij; :~; ;E-JU; ;~,; ;~; ;~~=:: ;~qe..__. -r----

1988.650 --- II ----r---

; ; ~ ;;;;;~;;;~:;,;;: ;;); I 

1988.775 ---I 

:__u~ ; ;~; ~~= ; ;~~ ;:..JLi£.; ; -i£-il-i~ ; ,~;;; i i~l i&---,--1 --r--_ 

1 9 8 8 • 9 0 0 ~-._,._ _ __, 

FIGURE 5 (Continued).- Future Columbia glacier behavior, projected by using modeling constant a= 1165625. The individual 

illustrations show the flow-line profile from column locations 1 through 10 as a function of time. 
B20 STUDIES OF COLUMBIA GLACIER, ALASKA 

; ~ ; ; ~ ; ; ; ; ; ; ; ; ; ;1;1;1;;; :~~ 1 

I I I 

------ I 1 2 3 4 1 989. 025 7 I I 8 9 10 

~ ~ ; ; : ; ; ; ; ; ; ; ; ; ; ; ; ~ ; ; ;;;; 

I I I 

1

1 9 8 9 · 

1

1 5 0 ------ I 

~ ; ; ; : ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; :;;; 

r I I ;989.~75 ----

; ; ; ; ~;;; ~;;;~;,;;;;;;;a; I 

------ I 1989.400 

;;;;~;:;~;;;:;~;;;;;Iff 1 

------ I 1989.525 

;;;;~;;;:;:; ;, ;;; : :r; I 

----- I 1989.650 

1 9 8 9 • 7 7 5 --r----r-----""T 

!_U~;;Ul.if,;~~;;;~;~.i£$;;;~;,;_.;:;;:~;:1; ~I ----r---1 ---... 

1 9 8 9 . 9 0 0 II-----,----. 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using modeling constant a=1165625. The individual 

illustrations show the flow-line profile from column locations 1 through 10 as a function of time; 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B21 

1 2 3 4 1 9 9 0 • 0 2 5 7 8

11 -___,.._

9

, -------., 

10 

1990.150 

1990.400 

r 1990.525 

~~:-~:;~:~~;;:~;~~;;~:;;~::;:~~=~,~~~~~~ 1990.650 r---~--~ 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using modeling constant a= 1165625. The individual 

illustrations show the flow-line profile from column locations 1 through 10 as a function of time. 
B22 STUDIES OF COLUMBIA GLACIER, ALASKA 

1 2 3 4 1991.025 7 'k----~~----~1 9 10 

1991.150 

1991.275 ~~--~--~ 

;; :; ; ; ; ~ :: : : : : ; .. 1i--i: .,.--....-I --r-1 ---.----... -r--r--------. 

1991.525 

FIGURE 5 (Continued).-Future Columbia glacier behavior, projected by using modeling constant a=1165625. The individual 

illustrations show the flow-line profile from column locations 1 through 10 as a function of time. 
en 

a: 

w 

1-

w 

:E 

0 

..J 

~ 

z 

w 

z 

..J 

~ 

g 

LL. 

CJ 

z 

0 

..J 

<( 

w 

0 

z 

<C 

1-

en 

0 

14.5 

14.0 

13.5 

13.0 

12.5 

12.0 

FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B23 

OBSERVED 

I 

\ -, I I 

\ ,_./ \ .A:: 

\ I \ I "' '-.// 

/ 

START OF MODELING, 

ALL THREE CASES 

14 

C/) CONTINUATION a:1165625 a: 

w 

~12 

:::11 

0 

_, 

:;;: 10 

~ 

w· 

~ 8 _, 

3:: 

0 

it 6 

C) 

z 

0 

~ 4 

w 

(.) 

z 

<( 

1-

C/) 

0 

19 

15 a: 

w 

al 

:::11 

13 ~ 

z 

11 ~ _, 

0 9 (.) 

w 

z 7 :::i 

I 

3:: 

0 5 _, 

LL 

3 

\_.a:1168750 

. 

. 

\ . 

• 

19.0~ 

al 

:E 

::J 

z 

18.5 z 

:E 

::J 

..J 

0 

18.0° 

w 

z 

:J 

I 

17.5 ~ ...1 

LL. 

11.5~----~------~----~------~----~------------~~----~----~ 

1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 

FIGURE 6.-Terminus position in time, for three values of the modeling constant a of the calving law. The observed position represents an average 

across the width of the terminus. 
B24 STUDIES OF COLUMBIA GLACIER, ALASKA 

160r------------------------------------------------------------------

140 

120 a=1162500 

tn 

a: 100 w 

1-

w 

:t 

z 

....... 

J: 

(!J 

w 

J: 

c 

w 

1-

a: 

0 

a. 

0.. 

;::) 

tn 

z 

;::) 

80 

60 

40 

20 

OBSERVED 

\I' I \ 

I \ " \ I \ I \ 

\ rJ \ l ' r, 'J ' ,, \ 

'J \ 

START OF MODELING, 

ALL THREE CASES 

~ 

" I ' I ,,,...., 

I 

I 

I 

a:1165625 

a=1168750 

0~------~------~------~-------L------~------~------~~------~----~ 

1976 1977 1978 1979 1980 1981 1982 1983 1984 

FIGURE 7.-Unsupported height, h,., at the terminus, in time, for three values of the modeling constant a of the calving 

law. The observed is at the head of the embayment. The anomolous peak in this curve in early 1979 is due to a sudden, 

temporary switch of the location of the subglacial river discharge to a relatively undisturbed location on the terminus 

1985 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B25 

r----------------------------------------------------------160 

'~984 

140 

120 

100 

80 

60 

40 

20 

0 

1985 1986 1987 1988 1989 1990 1991 1992 

with a higher ice cliff. The dotted curve represents approximately what would have resulted had this temporary shift 

in embayment location not occurred. The greater variability of the observed curve compared to that of the calculated 

curve probably is due also to the shifting of the subglacial discharge across the terminus. 
B26 STUDIES OF COLUMBIA GLACIER, ALASKA 

a: 

<( 5 w 

> 

a: 

w 

c. 

fJ) 4 a:1178650 

a: 

w 

1-

w 

:E 

0 3 

-I 

~ START OF MODELING, 

z (ALL THREE CASES 

2 OBSERVED\ r, Q '\ I \ 

w ,., I \ I \ 

w '\ I \ I \ 

0. I "'"'"" I \ ' fJ) I \I\ \ I \ 

I " I \ ' I ' -I 

<( ' I ' J 1 I \J 

1-

z a=116;;;:7 0 

N 

a: 

0 0

%: 1976 1984 1985 

FIGURE 8.-Surface horizontal speed at flow-line-column number 16.9, in time. A phase shift of unknown origin (a three-dimensional effect?) between computed speed and observed speed indicates computed speed leading the observed speed, (V.)0 , by about 1

/s year. Note the low speeds 

in 1980, when the retreat was interrupted. 
FINITE ELEMENT GLACIER DYNAMICS MODEL APPLIED TO COLUMBIA GLACIER, ALASKA B27 

30 

a: 28 

<( 

w 

> 26 

a: 

w 

a. 24 

en 

a: 22 

w 

.... 

w 20 ~ 

0 

..J 18 ~ 

0 16 

m 

::l 

0 14 

z 

12 .. X 

::l 10 ..J 

u. 

~ 8 z 

> 

..J 

<( 

0 

0;---~--~--~--~--~~--~------------------~--~--~----~--T-~ 

76 77 78 79 1980 81 82 83 84 85 86 87 88 89 1990 91 92 

FIGURE 9.-Rate of iceberg production, in volume, projected by using modeling constant a= 1165625. This graph was derived by converting the 

two-dimensiol).al flux of the model to three dimensions, using bed configuration and spatial velocity distributions for the real glacier. 
B28 STUDIES OF COLUMBIA GLACIER, ALASKA 

REFERENCES CITED . 

Bindschadler, R. A., 1978, A time-dependent model of temperate 

glacier flow and its application to predict changes in the surge-type 

Variegated Glacier during its quiescent phase: Seattle, University 

of Washington Geophysics Program, unpublished thesis, ·244 p. 

Budd, W. F., Keage, P. L., and Blundy, N. A., 1979, Empirical studies 

of ice sliding: Journal of Glaciology, v. 23, no. 89, p. 157-170. 

Hooke, R. L., Gustafson, R. J., Raymond, C. F., and Hotchkiss, R. L., 

1980, Computer programs for calculating temperatures and· 

velocities in polar glaciers using the finite element method: Minneapolis, University of Minnesota, Report UM-GMPH-2, 95 p. 

Meier, M. F., Post, Austin, Brown, C. S., Frank, David, Hodge, S.M., 

Mayo, L. R., Rasmussen, L.A., Senear, E. A., Sikonia, W. G., Trabant, D. C., and Watts, R. D., 1978, Columbia Glacier progress 

report-December 1977: U.S. Geological Survey Open-File Report 

78-264, 56 p. 

Meier, M. F., Post, Austin, Rasmussen, L.A., Sikonia, W. G., and 

Mayo, L. R., 1980a, Retreat of Columbia Glacier-A preliminary 

prediction: U.S. Geological Survey Open-File Report 80-10, 12 p. 

Meier, M. F., Rasmussen, L.A., Post, Austin, Brown, C. S., Sikonia, 

W. G., Bindschadler, R. A., Mayo, L. R., and Trabant, D. C., 

1980b, Predicted timing of the disintegration of the lower reach of 

Columbia Glacier, Alaska: U.S. Geological Survey Open-File 

Report 80-582, 58 p. 

Post, Austin, 1975 [1976], Preliminary hydrography and historic terminal changes of Columbia Glacier, Alaska: U.S. Geological Survey 

Hydrologic Investigations Atlas 559, 3 sheets. 

Raymond, C. J., 1978, Numerical calculation of glacier flow by finite 

element methods: Seattle, University of Washington, final technical report for National Science Foundation grant no. 

DPP74-19075, geophysics program AK50, 201 p. 

Schmidt, W. F., 1977, Finite. element analysis of planar creeping 

viscous ice flow: Orono, University of Maine, Department of 

Mechanical Engineering, Report No. 77-1, 36 p. 

Sikonia, W. G., and Post, Austin, 1980, Columbia Glacier, Alaska: 

Recent ice loss and its relationship to seasonal terminal embayment, thinning, and glacier flow: U.S. Geological Survey Hydrologic Investigations Atlas 619, 3 sheets. 

Taylor, R. L., 1980, Finite element analysis program, in Zienkiewicz, 

0. C., The finite element method (3d ed.): Berkeley, University of 

California, Department of Civil Engineering •. [Available on cards 

(for $100 as of 1980) from R. L. Taylor at the univ.ersity.] 

Zienkiewicz, .. O: C., 1977, The finite element method .(3d ed~): 

McGraw-Hill, 787 p; 
APPENDIX A: MATHEMATICS 

The formulation used is that of Zienkiewicz (1977, p. 

612, sec. 22.3.1). The K and M matrices can be omitted 

as they represent contributions from kinetic effects that 

can be ignored in the slow movement of glaciers. 

Quadratic Lagrange shape functions were used for 

velocities, and so each element is three nodes by three 

nodes. Pressure shape functions were chosen to be 

linear in the vertical direction and quadratic in the horizontal direction, and so pressure unknowns are 

associated with only the top and bottom nodes of each 

element; that is, with nodes number 4, 7, 3, 1, 5, and 2 in 

figure 24.25 of Zienkiewicz (1977, p. 748). Three-bythree Gaussian intregration is used. The starting point 

for the element subroutine is Taylor's fluid element 

(Zienkiewicz, 1977, p. 755). The viscosity, p., for ice is 

nonconstant, and so the resultant system of equations is 

solved iteratively and has a contribution to the tangent 

stiffness matrix (Zienkiewicz, p. 699, equation 24.2) arising from appropriate partial derivatives of p.. 

TANGENT STIFFNESS MATRIX 

Assume for simplicity that a", aP are partitioned within 

an element as in section 22.3.1 of Zienkiewicz (1977, p. 

612). (Actually, in the· program, the velocities and 

pressures alternate, grouped according to node 

number.) In Zienkiewicz' notation, the equation for ice is 

(see equation 22.17) 

F(a) E Ka" + KPaP + fu = 0, (Al) 

where 

(A2) 

This equation is nonlinear because p. (in K) is a function 

of the effective shear stress, and so of a". To solve 

iteratively by Newton's method, expand to the linear 

term in Taylor's series: 

F(a)+F'(a)Oa::O (A3) 

or 

(A4) 

Here a is the trial, or "current," approximation of the 

solution of equation Al, and oa represents a correction 

to be applied to improve this approximation. Considering the first term in which 

Ka" = [! BTO'BdO]a", (A5) 0 

define 

G= _!_ BTO'B. (A6) 

p. 

Differentiating, 

~Ga")] oa"=[ap. oa"]Ga"+p.Goa". (A7) aau aau 

Note that~ is a row matrix, while aa" and Ga" are colaa" 

umn matrices. Because of the special nature of these 

matrices, the order of operations in the first term on the 

right can be rearranged to give 

[~Ga·>] aa·=[ Ga·::. +,.G ]aa•. (A8) 

Thus equation A4 becomes 

[!Ga·~O+Ka·]aa•+K•6a•e -F(a). (A9) 

The integral represents the contribution due to the nonconstant viscosity, p.. To express it in workable form, 

note that 

ap. ap. ad" -=--' 

aau ad" aau 

(AlO) 

where 

(All) 

(summation convention on subscripts) is the effective ~ 

shear deformation rate. Using the notation of Zienkiewicz, section 22.3.1, and writing column matrix din 

place of d1" 

we have 

where 

d=Ba", 

d" = lf4(Ba")TE(Ba"), 

E•ln. p. 

(Al2) 

(Al3) 

(Al4) 

B29 
B30 STUDIES OF COLUMBIA GLACIER, ALASKA 

Taking the derivative, 

ad2 

aa"fJa" = 1/4[(fJa")T8TEBa" + a"TBTEBfJa"], 

or rearranging, 

that is, 

(Al5) 

(Al6) 

From equations 1, 6, and All, the ice viscosity, p., can be 

written in terms of the effective deformation rate, d, as 

(dl)((l-n)/ln) • 

p. 2A 1/11 (Al8) 

This yields 

ap. (1- n)d< (1/n)-l) 

ad'}.= 4nA11n (Al9) 

Thus, the integral term in equation A9-the contribution 

in Newton's method due to nonconstant p.-is 

[! 'h~Ga•)(Ga•)•dll ]oa•. (A20) 

The column matrix 

(A21) 

(defining Q as found in the element subroutine, not the Q 

of equation 24.10 of Zienkiewicz (1977, p. 710)) is evaluated by using the "current" velocities, a", of the iterative 

Newton's method process. Q must be calculated as part 

of accounting for "current" internal forces (Zienkiewicz, 

equation 24.2, p. 699), so the dependence of the tangential matrix on p. can be obtained with little additional 

effort. 

Some distribution of p. must be chosen initially to start 

the iterative process. In the program, this initial p. is 

assumed to be given by equation 6 with effective shear 

stress, T, identically 1 bar throughout the ice mass. 

UPPER FREE-SURFACE DYNAMICS 

To follow the motion of the glacier surface in time, the 

ice velocity along the surface, r ., is resolved into 

generally nonorthogonal components along the surface 

and in the vertical direction; this vertical component is 

called the emergence velocity, V •. If 'Y is the angle between horizontal and an upward normal to the surface, 

define 

(A22) 

then 

(A23) 

For V. evaluated on the surface, the only nonzero elements in N" are those which correspond to surface 

nodes. Let h represent the (unknown) rate of thickness 

change and suppose that 

(A24) 

where the scalar shape function elements in N11 correspond to scalar velocity shape function elements for the 

surface nodes, for evaluation on the surface only; that is, 

they corresP.ond to the nonzero elements of N• in equation A23. If b is the balance rate, then at any point on the 

surface, 

h- v.-b=O. (A25) 

To obtain a suitable numerical equation, use the 

Galerkin process and premultiply equation A25 by (fJh)r, 

then integrate over the surface: 

(fJa11)rJ N11T (h- v.- b) dr.=O. (A26) 

r. 

Define 

K</=aJ (N,A)IN/dr., r. 

KIJ"= -aJ (N~11 )rsl"N/dr., (A27) r. 

f111 = -aJ (N~11 )r6dr., r. 

where b is actually stored for the surface nodes and interpolated along the surface as h is in equation A24 and 

where a is a constant to be explained in the next paragraph. Thus, equation A26 is 

K11a11 +K"a"+f"=0. (A28) 

The equation which replaces equation 22.22 of Zienkiewicz (1977, p. 613) is as follows: 

[

K

11 K" 0] {a11

KuT K KP a" } + {f"} f = 0. (A29) 

0 KpT 0 aP 0 

K"r was added to the stiffness matrix to preserve symmetry, and the value assigned to constant a of equations 

A27 was chosen to be a very small number (l0-16) so that 

the K"r would not significantly alter the second (set of) 

rows of equation A29. 

Here the unknowns all, a", and ap have been partitioned into groups for simplicity, but, in the program, 
APPENDIX A: MATHEMATICS B31 

they alternate, grouped by node number. The unknowns 

ah occupy some of the third nodal degree of freedom 

(DOF) not needed for pressure unknowns ap (since 

pressure is only linear in the vertical direction). The 

third DOF of the second node from the surface in a column, thus, is assigned the unknown a;h, which then represents the rate of surface altitude change at the current 

column location. 

The tangent stiffness matrix has Kh, K", and K"T added 

in the same locations as in equation A29. The right side 

of the Newton's method equation must include the contribution arising from "current" v .. ,b, and h, 

P;h=a! (N;h}T(V .. +b-h)dl\. (A30) r. 

(Caution: Pas defined by Taylor in the element routines 

is the entire right side of the Newton's method equation; 

that is, it is - f-P of equation 24.2 or equation 18.10 

(Zienkiewicz, 1977, p. 699, 453).) 

SURFACE TRACTIONS AND BODY FORCES 

Water and atmospheric pressure (surface tractions) 

and gravity (a body force) are added to the equations, as 

in equation 22.16 of Zienkiewicz (1977, p. 612). 

UPDATING THE MESH 

If h6 denotes the ice surface altitude and hb is the bed 

altitude and 

h = h6 - hb meters, (A31) 

the eight nodes within a vertical column have been 

chosen to be at locations h6 , h6 -0.4 (h-10), h6 -0.8 

(h-10), h6 -0.9 (h-10), hb+ 10, hb+ 10, hb+5 and hb, with 

a greater concentration of nodes near the bed where 

deformation rates are higher. (In an infinitely wide slab 

glacier with n= 3, deformational velocity varies through 

the thickness as a quartic; the program would approximate this quartic with two quadratics.) 

The upper free surface and calving face are updated as 

in section 22.6 of Zienkiewicz (1977, p. 623), using ah 

calculated as described in the section of Appendix A entitled "Upper Free-Surface Dynamics" and the calving 

. speed from the calving law of equation 8. Because of 

large deformations, the entire mesh is updated at each 

step. The updating is done by calculating glacier length 

at each time step and then by requiring that each column 

~ of nodes be at the same fractional location within this 

new length as it was initially. Nodes are then distributed 

within each column as a fraction of the new thickness, as 

described in the preceding paragraph. 

FUTURE MODIFICATION 

The program was written to allow easy future 

modification to three space dimensions; a large part of 

this change can be done by simply replacing twodimensional shape functions by three-dimensional ones. 
APPENDIX B: USER INSTRUCTIONS FOR THE COMPUTER PROGRAM 

The general purpose finite element program 

resembles closely that described by Taylor in chapter 24 

of Zienkiewicz (1977); only relevant changes and 

specialized modules are described in this appendix. The 

program now (1980) available from Taylor contains a 

data prerocessor that converts free format input to 

fields always of length 10-if an integer is expected, the 

field is 110; if a FORTRAN real is expected, the field is 

F10.0. 

COORDINATE DATA 

The program assumes that nodes are specified in vertical columns. Each column is, to some extent, considered a unit. The sample run of Appendix C (p. 

B34-B46) uses columns of 8 nodes (the user can select 

more elements per column if he chooses). The third and 

fourth nodes frQm the bottom in a column are at the 

same location on the column and have the velocity DOF 

tied together, allowing only a pressure discontinuity to 

accommodate the discontinuity in the A flow law parameter at that point. Extra "dummy" columns can be specified beyond the initial terminus position to accommodate 

an advancing glacier. The NUMNP on the control card 

should be the total number of nodes including these extras. The true terminus position is as specified on the 

material property data card. 

MATERIAL PROPERTY DATA 

Columns Description Variable 

01 to 10 _____ A of ice-flow law------------------- D(1,1) 

11 to 20 _____ n of ice-flow law------------------- D(2,1) 

21 to 30 _____ Uppermost element in terminal face___ D(3,1) 

31 to 40 _____ Number of nodes in each column______ D(4,1) 

BOUNDARY RESTRAINT DATA 

Boundary Condition Cards-FORMAT(B/10).- Use 

only 0, + 1, and -1 for the constraint conditions that 

Taylor describes. A boundary code describing tieing is 

the node number of a preceding like-numbered DOF. 

For example, to tie DOF 1 and DOF 2 of node 14 to DOF 

1 and DOF 2 of node 13, specify the following boundary 

condition codes: 

B32 

14 

198 

8 

0 

13 

197 

13 

197 

The tieing uses the generation increment also to tie the 

first two degrees of freedom for node 22 to 21, node 30 

to 29, ... , node 198 to 197. When tieing is specified, only the DOF of the node tied to is retained in the system 

of equations. Do not tie degrees of freedom that have 

specified displacements. 

Unused degrees of freedom should have displacements 

specified. Thus, the middle row of nodes in each element, except for the surface elements, should have 

specified displacements for the third DOF, because 

pressure is only modeled linearly in the vertical direction. (The middle row of nodes of the surface elements 

have the h unknown as the third DOF.) Also, the extra 

"dummy" nodes beyond the initial terminus position 

should have all degrees of freedom with displacements 

specified. 

NODAL FORCED BOUNDARY DATA 

Boundary forces F(3, ·) corresponding to unused 

degrees of fre~dom in a column are used to store the 

balance rate, b, and sliding flow law parameter, Ab. 

These are considered to be at the location of the column 

when dynamic modeling changes the column position. 

Within a column, storage is as follows: 

"Current" parameters F(3, ·) Location in column 

------------------------- 1 (surface) 

------------------------- 2 ------------------------- 3 

0-------------------------- 4 ------------------------- 5 

------------------------- NDEEP-2 

Ab ------------------------ NDEEP-1 ------------------------- NDEEP (bed) 

For example, if the second node from the bottom in a 

column is node number 55, the sliding flow law 

parameter, Ab, for the current position of the column is 

stored in F(3,55). 

NODAL TEMPERATURE DATA 

In a manner similar to column storage of parameters 

in F(3, • ), nodal "temperatures," T( • ), are used to store 

information corresponding to the initial configuration of 

the ice mass. Note the distinction that the T( ·) refer to 

fixed initial column coordinates, while the F(3, ·) refer to 
APPENDIX B: USER INSTRUCTIONS FOR THE COMPUTER PROGRAM B33 

current dynamic column coordinates. Within a column, 

storage is as follows: 

"Initial" parameters, T(") Location in column 

;: =============================== ~ (surface) ~b (future use)--------------------- 3 

~-------------------------------- 4 

h (future use)---------------------- 5 

V. ------------------------------1 NDEEP-2 

Ab ------------------------------- NDEEP-1 

------------------------- NDEEP (bed) 

Here Xb, Yb, .and Zb are the bed coordinates; b, the 

balanc~ rate; h, the rate of surface altitude change; Vs, 

the horizontal surface speed; and Ab, the sliding flow law 

parameter- all representing measured or estimated 

values for the initial column position. For example, if the 

second node from the surface in a column is node 

number 50, the y coordinate of the bed at the initial column position is stored in T(50). 

TANF MACROINSTRUCTION 

A few additional macroinstructions have been added 

to Taylor's set. TANF simply combines the functions of 

his TANG and FORM. Proceeding this way was more 

efficient for the case at hand and is, in fact, now required 

· by the element subroutine ELMT04 for ice flow rather 

than the TANG and FORM routines. 

REST AND SAVE MACROINSTRUCTIONS 

Program restart capability is provided through REST 

and SAVE macroinstructions. SAVE writes the node 

number, nodal coordinates, and nodal degrees of freedom to FORTRAN file 15 with 

FORMAT(I16,4(1PE16.9)) 

or 

FORMAT(I16,4(1PE16.9)/5(1PE16.9)) 

depending on the number of coordinates and degrees of 

freedom. Macro REST reads this file as FORTRAN file 

16. If 

REST 1. 

is specified, where the 1. is in columns 21 to 30 the coor- . ' dinates, b, Ab, and the degrees of freedom are updated 

from the rest~rt file. If 0. is given in place of the 1., the 

coordinates, b and Ab are left as on the initial input data, 

and only the degrees of freedom are updated. 

UPNO MACROINSTRUCTION 

The mesh is updated as in section 22.6 of Zienkiewicz 

(1977, p. 623) by using the macroinstruction (UPNO 

(mnemonically, "update nodes"). The upper surface is 

adjusted according to ice velocity vectors and balance 

rate, and the calving face is adjusted according to the 

calving law of equation 7. 

PLOT MACROINSTRUCTION 

The PLOT macroinstruction produces a plot giving 

surface, bed, and calving face outlines and velocity vectors. The plot routine expects scaling information immediately following a SCAL macro card that is placed 

after the END macro that terminates MACR processing. The scaling parameters are as follows: 

Columns Description Variable 

01 to 10 ___ Minimum x-coordinate in plot----------- XMIN 

11 to 20 ___ Maximum x-coordinate in plot----------- XMAX 

21 to 30 ___ Minimum y-coordinate in plot----------- YMIN 

31 to 40 ___ Maximum y-coordinate in plot----------- YMAX 

41 to 50 ___ Scale factor in units per inch------------ SCALE 

51 to 60 --- Vector scale factor in speed units per inch _ VSCALE 

61 to 70--- Arrowhead length in inches------------ AHLEN 

71 to 80 ___ Plot rotation, 0 or 90 degrees (only)------ ROT 

This scaling parameter card is followed by an END 

macro; that is, the plot routine reads the SCAL card, the 

parameter card, and an END card. 
APPENDIX C: NEW SUBROUTINES 

The following subroutines were written by the author specifically to treat the modeling of glaciers. 

SUBROUTINE ELMT04(D,UL,XL,IX,TL,S,P,NDF,NDM,NST,ISW) 

c 

C •••• 'IWO DIMENSIONAL ELEMENT FOR ICE FLOW 

c 

COMMON /CDATA/ O,HEAD(20),NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR 

COMMON /ELDATA/ DM,N,MA,MCT,IEL,NEL 

COMMON /RSTART/ !START 

COMMON /JACOBI/ XS(2,2),SX(2,2) 

COMMON M(1) 

COMMON /POINTR/ NO,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13, 

* N14,NE 

LOGICAL SFLAG 

DIMENSION XL(NDM,1),UL(NDF,1),P(NDF,1),S(NST,1),D(1),V(3),DV(2,2) 

1,SHP(3,9),IX(1),SG(9),TG(9),WG(9),SHPP(3,9),TL(1) 

DIMENSION SG1D(3),TG1D(3),WG1D(3),IANODE(9),IXP(9),INTERP(9) 

DIMENSION ITOP(3),MID(3),Q(27),ANODE(3) 

C •••• !ANODE GIVES VISCOSITY "A" INTERPOLATION (REDEFINITION OF NODES). 

C •••• INTERP YIELDS PRESSURE IN,TERPOLATION (1=NODE USED, O=NOT). 

DATA IANODE/1,2,2,1,3,2,3,1,3/,INTERP/1,1,1,1,1,0,1,0,0/ 

DATA ITOP/3,4,7/,MID/6,8,9/,Q/27*0./ 

GO TO (1,2,3,4,5,3),ISW 

C •••• INPUT/OUTPUT FLUID PROPERTIES 

C •••• LINT1D IS THE NUMBER OF INTEGRATION POINTS IN EACH DIRECTION 

C •••• D(1) =A (ICE), D(2) = N 

C •••• D(3) IS THE LOWEST NUMBERED NODE IN THE CALVING FACE -

C •••• ALL NODES FROM D(3) TO NUMNP ARE ASSUMED TO BE IN THE FACE. 

C •••• D(4) IS THE NUMBER OF NODES IN A VERTICAL COLUMN THROUGH DEPTH. 

B34 

1 READ(5,1000) (D(I),I=1,4) 

WRITE ( 6 , 2 0 0 0) ( D (I) , I= 1, 4) 

LINT1D = 3 

LINT = 0 

!START = 0 

F1EPS = 1./D(2) 

F2EPS = 1./D(2) - 1. 

RHOG = 0.0882 

RHOGW = 0. 0980 

A'!MOSP = 1.013 

ALPHA = 1.E-16 

NFACE = D(3) + 0.1 

NDEEP = D(4) + 0.1 

ELM 001 

ELM 002 

ELM 003 

ELM 004 

ELM 005 

ELM 006 

ELM 007 

ELM 008 

ELM 009 

ELM 010 

ELM 011 

ELM 012 

ELM 013 

ELM 014 

ELM 015 

ELM 016 

ELM 017 

ELM 018 

ELM 019 

ELM 020 

ELM 021 

ELM 022 

ELM 023 

ELM 024 

ELM 025 

ELM 026 

ELM 027 

ELM 028 

ELM 029 

ELM 030 

ELM 031 

ELM 032 

ELM 033 

ELM 034 

ELM 035 

ELM 036 

ELM 037 

ELM 038 

ELM 039 

ELM 040 
APPENDIX C: NEW SUBROUTINES 

DO 20 !=1,9 

IXP(I) = 0 

IF(INTERP(I).EQ.1)NELPC=I 

20 SHPP(3,I) = 0. 

RETURN 

2 RETURN 

C •••• COMPUTE SYMMETRIC TANGENT STIFFNESS AND OUT OF BALANCE FORCES. 

C. • • • USE TANF, NOT TANG AND FORM -- ONLY ONE PASS PER ELEMENT. 

3 !START = !START + 1 

C •••• REDEFINE NUMNP TO EXCLUDE "EXTRA" NODES BEYOND CALVING FACE. 

NUMNP = NFACE + NDEEP - 1 

L = LINT1D 

IF(L*L.NE.LINT) CALL PGAUSS(L,LINT,SG,TG,WG) 

DO 65 L=1,LINT 

CALL SHAPE(SG(L),TG(L),XL,SHP,XSJ,NDM,NEL,IX,.FALSE.) 

XSJWG = XSJ*NG(L) 

C. • • • CALL PRESSURE SHAPE FUNCTIONS -- • TRUE. BECAUSE PRESSURE 

C. • • • GRADIENT NOT NEEDED, AND WOULD NOT BE CORRECT ANYWAY 

C. • • • BECAUSE NOT ISOPARAMETRIC. 

NELP = MINO(NELPC,NEL) 

DO 25 I=1,NELPC 

25 IXP(I) = IX(I)*INTERP(I) 

CALL SHAPE(SG(L),TG(L),XL,SHPP,XSJP,NDM,NELP,IXP,.TRUE.) 

C •••• COMPUTE VELOCITIES, VELOCITY GRADIENTS, PRESSURES = V(3), 

C •••• AND FLOW LAW PARAMETER A (ENTERED AS NODAL "TEMPERATURE"). 

C •••• USE ONLY LOWER NODES TO ALLOW VISCOSITY DISCONTINUITY IN DEPTH, 

C •••• BUT STILL RETAIN CONTINUITY ALONG FLOW DIRECTION. 

DO 32 I = 1,2 

V(I) = 0. 

DO 31 K = 1,NEL 

31 V(I) = V(I) + SHP(3,K)*UL(I,K) 

DO 32 J = 1,2 

DV (I , J) = 0 • 0 

DO 32 K = 1,NEL 

32 DV(I#J) = DV(I,J) + SHP(J,K)*UL(I,K) 

V(3) = 0. 

DO 280 K=1,NELP 

280 V(3) = V(3) + SHPP(3,K)*UL(3,K) 

A = 0. 

C •••• CHECK IF AT BED. 

IF(MOD(IX(1),NDEEP).NE.O)GO TO 283 

C •••• COMPUTE "A" FLOW LAW COEFFICIENT ALONG BED TO SIMULATE SLIDING. 

DO 281 K=1,3 

281 CALL VALUES(ANODE(K),M(N10),3,IX(MID(K)),NDF) 

DO 282 K=1,NEL 

282 A= A+ SHP(3,K)*ANODE(IANODE(K)) 

283 A= AMAX1(D(1),A) 

C •••• COMPUTE EFFECTIVE STRAIN RATE EPSEFF 

EPSEFF = 0. 

DO 290 !=1,2 

B35 

ELM 041 

ELM 042 

ELM 043 

ELM 044 

ELM 045 

ELM 046 

ELM 047 

ELM 048 

ELM 049 

ELM 050 

ELM 051 

ELM 052 

ELM 053 

ELM 054 

ELM 055 

ELM 056 

ELM 057 

ELM 058 

ELM 059 

ELM 060 

ELM 061 

ELM 062 

ELM 063 

ELM 064 

ELM 065 

ELM 066 

ELM 067 

ELM 068 

ELM 069 

ELM 070 

ELM 071 

ELM 072 

ELM 073 

ELM 074 

ELM 075 

ELM 076 

ELM 077 

ELM 078 

ELM 079 

ELM 080 

ELM 081 

ELM 082 

ELM 083 

ELM 084 

ELM 085 

ELM 086 

ELM 087 

ELM 088 

ELM 089 

ELM 090 
B36 STUDIES OF COLUMBIA GLACIER, ALASKA 

DO 290 J=1,2 

290 EPSEFF = EPSEFF + (DV(I,J) + DV(J,I))**2 

EPSEFF = SQRT(EPSEFF/8.) 

C •••• NOW COMPUTE VISCOSITY AMU --USE AMUCON IF 

C •••• FIRST PASS (ALL ELEMENTS) THROUGH TANF -- INDICATED 

C •••• BY FLAG ISTART.LE.NUMEL. 

C •••• FOR FIRST PASS USE VISCOSITY AT BED WITH ~AUBED = 1. 

C •••• DMU IS DERIVATIVE OF VISCOSITY MU W.R.T. EPSEFF SQUARED--

C •••• USED IN NEWTON'S METHOD -- DMU = 0. FOR FIRST PASS. 

AMU = 0.5/A 

DMU = 0. 

IF(ISTART.LE.NUMEL)GO '.00 291 

AMU = (0.5/A**F1EPS)*EPSEFF**F2EPS 

DMU = (0.25*F2EPS/A**F1EPS)*EPSEFF**(F2EPS-2.) 

291 XMU = AMU*XSJWG 

XDMU2 = DMU/2. *XSJWG 

C •••• COMPUTE DIVERGENCE TERM 

DIV = DV(1,1) + DV(2,2) 

GAMXY = DV(1,2) + DV(2,1) 

C. • • • COMPUTE INTERNAL FORCES 

C. • • • Q RELATES 'ID NEW'IDN' S METHOD ACCOUNTING OF THE DERIVATIVE OF 

C •••• VISCOSITY. 

J1 = 1 

DO 64 J=1,NEL 

Q(J1) = 2.*DV(1,1)*SHP(1,J) + GAMXY*SHP(2,J) 

P(1,J) = P(1,J) + SHP(1,J)*V(3)*XSJWG- Q(J1)*XMU 

Q(J1+1) = 2.*DV(2,2)*SHP(2,J) + GAMXY*SHP(1,J) 

P(2,J) = P(2,J) + SHP(2,J)*V(3)*XSJWG- Q(J1+1)*XMU 

* - XSJWG*RHOG*SHP(3,J) 

C. • • • THE LAST TERM ABOVE IS GRAVITY LOADING 

P(3,J) = P(3,J) + SHPP(3,J)*DIV*XSJWG 

J1 = J1 + NDF 

64 CONTINUE 

C •••• COMPUTE TANGENT, LOOP OVER COLUMNS OF S 

K1 = 1 

DO 34 K = 1,NEL 

A1 = XMU*SHP(1,K) 

A2 = XMU*SHP(2,K) 

QK1 = Q(K1)*XDMU2 

QK1P1 = Q(K1+1)*XDMU2 

C •••• LOOP OVER ROWS OF S 

J1 = 1 

DO 33 J = 1,K 

S(J1 ,K1 ) = S(J1 ,K1 + SHP(1,J)*(A1+A1)+SHP(2,J)*A2 

* + Q(J1)*QK1 

S(J1 ,K1+1) = S(J1 ,K1+1) + SHP(2,J)*Al + Q(J1)*QK1P1 

S(J1 ,K1+2) = S(J1 ,K1+2) - SHP(1,J)*SHPP(3,K)*XSJWG 

S(J1+1,K1 ) = S(J1+1,K1 + SHP(1,J)*A2 + Q(J1+1)*QK1 

S(J1+1,K1+1) = S(J1+1,K1+1) + SHP(1,J)*Al+SHP(2,J)*(A2+A2) 

* + Q(J1+1)*QK1P1 

ELM 091 

ELM 092 

ELM 093 

ELM 094 

ELM 095 

ELM 096 

ELM 097 

ELM 098 

ELM 099 

ELM 100 

ELM 101 

ELM 102 

ELM 103 

ELM 104 

ELM 105 

ELM 106 

ELM 107 

ELM 108 

ELM 109 

ELM 110 

ELM 111 

ELM 112 

ELM 113 

ELM 114 

ELM 115 

ELM 116 

ELM 117 

ELM 118 

ELM 119 

ELM 120 

ELM 121 

ELM 122 

ELM 123 

ELM 124 

ELM 125 

ELM 126 

ELM 127 

ELM 128 

ELM 129 

ELM 130 

ELM 131 

ELM 132 

ELM 133 

ELM 134 

ELM 135 

ELM 136 

ELM 137 

ELM 138 

ELM 139 

ELM 140 
S(J1+1,K1+2) 

S(J1+2,K1 ) 

S(J1+2,K1+1) 

33 J1 = J1 + NDF 

34 K1 = K1 + NDF 

65 CONTINUE 

= 

= 

= 

APPENDIX C: NEW SUBROUTINES 

S(J1+1,K1+2) - SHP(2,J)*SHPP(3,K)*XSJWG 

S(J1+2,K1 ) - SHPP(3,J)*SHP(1,K)*XSJWG 

S(J1+2,K1+1) - SHPP(3,J)*SHP(2,K)*XSJWG 

C •••• CHECK IF AT CALVING FACE WHERE ICOLUM=2. 

IF(IX(3).LT.NFACE)GO TO 302 

C •••• COMPUTE 1-D GAUSS POINTS AND WEIGHTS FROM 2-D FOR CALVING FACE. 

DO 301 L1D=1,LINT1D 

L2D = L1D 

IF(L1D.EQ.3)L2D=9 

SG1D(L1D) = 1. 

TG1D(L1D) = SG(L2D) 

301 WG1D(L1D) = SQRT(WG(L2D)) 

ICOLUM = 2 

DELTA= -2.*ATAN(1.) 

GO TO 304 

302 CONTINUE 

C •••• CHECK IF AT FREE SURFACE WHERE ICOLUM=1. 

IF(MOD(IX(3),NDEEP).NE.1)GO TO 66 

C. • • • COMPUTE 1-D GAUSS POINTS AND WEIGHTS FROM 2-D FOR FREE SURFACE. 

DO 303 L1D=1,LINT1D 

L2D = L1D 

IF(L1D.EQ.3)L2D = 9 

SG1D(L1D) = SG(L2D) 

TG1D(L1D) = 1. 

303 WG1D(L1D) = SQRT(WG(L2D)) 

ICOLUM = 1 

DELTA= 2.*ATAN(1.) 

C •••• COMPUTE WATER + A'.IMOSPHERIC FORCE AT INTEGRATION POINT ON FACE. 

304 DO 330 L1D=1,LINT1D 

C •••• CALL 2-D SHAPE FUNCTION WITH LAST ARGUMENT .TRUE.~ THIS RETURNS 

C •••• NATURAL (NOT GLOBAL) DERIVATIVES. 

CALL SHAPE(SG1D(L1D),TG1D(L1D),XL,SHP,XSJ,NDM,NEL,IX,.TRUE.) 

C •••• COMPUTE GLOBAL COORDINATES CORRESPONDING TO INTEGRATION POINTS. 

C •••• NATURAL DERIVATIVES FROM COMMON /JACOBI/ ARE 

C •••• XS(2,ICOLUM)=DY, XS(1,ICOLUM)=DX. 

C. • • • NOTE THAT XS IS THE TRANSPOSE OF THE JACOBIAN. 

YCOORD = 0. 

DO 305 K=1,NEL 

305 YCOORD = YCOORD + XL(2,K)*SHP(3,K) 

C •••• COMPUTE JACOBIAN= LINE LENGTH TRANSFORMATION 

XSJ = SQRT(XS(1,ICOLUM)**2 + XS(2,ICOLUM)**2) 

C. • • • COMPUTE ANGLE OF FACE, +X TO +NORMAL COORD. ALONG FACE (DY /DX) • 

BETA= ATAN2(XS(2,ICOLUM),XS(1,ICOLUM)) 

C •••• COMPUTE ANGLE TO OU'IWARD NORMAL TO FACE. 

GAMMA = BETA + DELTA 

C. • • • ADD TO INTERNAL FORCES 

DO 310 J=1,NEL 

ELM 141 

ELM 142 

ELM 143 

ELM 144 

ELM 145 

ELM 146 

ELM 147 

ELM 148 

ELM 149 

ELM 150 

ELM 151 

ELM 152 

ELM 153 

ELM 154 

ELM 155 

ELM 156 

ELM 157 

ELM 158 

ELM 159 

ELM 160 

ELM 161 

ELM 162 

ELM 163 

ELM 164 

ELM 165 

ELM 166 

ELM 167 

ELM 168 

ELM 169 

ELM 170 

ELM 171 

ELM 172 

ELM 173 

ELM 174 

ELM 175 

ELM 176 

ELM 177 

ELM 178 

ELM 179 

ELM 180 

ELM 181 

ELM 182 

ELM 183 

ELM 184 

ELM 185 

ELM 186 

ELM 187 

ELM 188 

ELM 189 

ELM 190 

B37 
B38 STUDIES OF COLUMBIA GLACIER, ALASKA 

WFORCE = (AMAX1(-YCOORD,O.)*RHOGW+ATMOSP)*XSJ*WG1D{L1D)*SHP(3,J) 

P{l,J) = P{1,J) - WFORCE*COS(GAMMA) 

310 P(2,J) = P{2,J) - WFORCE*SIN(GAMMA) 

C •••• CHECK IF AT UPPER FREE SURFACE, WHERE ICOLUM=1. 

IF(ICOLUM.NE.1)GO TO 330 

C. • • • ADD HDOT SUBMATRIX TO TANGENT STIFFNESS. 

C. • • • TERMS INVOLVING J2, K2, INCLUDE EMERGENCE CHANGE DUE TO VELOCITY 

C. • • • CHANGE CALCULATED IN THIS ITERATION. HDOT TERMS ARE MULTIPLIED 

C. • • • BY VERY SMALL ALPHA SO SYMMETRY CAN BE PRESERVED WITHOUT CHANGING 

C •••• VELOCITY CALCULATION. 

CTNGAM = COS(GAMMA)/SIN(GAMMA) 

DO 315 K=1,3 

K1 = MID(K)*NDF 

K2 = ITOP(K)*NDF 

SHK = ALPHA*SHP(3,ITOP(K))*XSJ*WG1D(L1D) 

DO 315 J=1,3 

J1 = MID(J)*NDF 

J2 = ITOP(J)*NDF 

SHJSHK = SHP(3,ITOP(J))*SHK 

S(J2-2,K1) = S(J2-2,K1) - CTNGAM*SHJSHK 

S(J2-1,K1) = S(J2-l,K1) - SHJSHK 

S(J1,K2-1) = S(Jl,K2-1) - SHJSHK 

S(J1,K2-2) = S(J1,K2-2) - CTNGAM*SHJSHK 

315 S(J1,Kl) = S{J1,K1) + SHJSHK 

C •••• COMPUTE BALANCE AND OLD HOOT (SO ITERATION YIELDS DELTA HDOT) 

C •••• AT INTEGRATION POINT. 

BALANC = 0. 

HDOT = 0. 

DO 316 J=l,3 

HDOT =HOOT+ UL(3,MID(J))*SHP(3,ITOP(J)) 

CALL VALUES(BALNOD,M(N10),3,IX(MID(J))+2,NDF) 

316 BALANC = BALANC + BALNOD*SHP(3,ITOP(J)) 

C. • • • COMPUTE SURFACE VELOCITY AT INTEGRATION POINT 

DO 320 I=1,2 

V(I) = 0. 

DO 320 J=1,3 

320 V(I) = V(I) + SHP(3,ITOP(J))*UL(I,ITOP(J)) 

C. • • • COMPUTE EMERGENCE 

EMERG = (V(1)*COS(GAMMA) + V(2)*SIN(GAMMA))/SIN(GAMMA) 

C •••• ADD TO RIGHT HAND SIDE 

DO 325 J=1,3 

325 P(3,MID(J)) = P(3,MID(J)) 

* + ALPHA*(EMERG + BALANC- HDOT)*SHP(3,ITOP(J))*XSJ*WG1D(L1D) 

330 CONTINUE 

IF(ICOLUM.EQ.2)GO TO 302 

C •••• FORM LOWER PART OF TANGENT STIFFNESS MATRIX BY SYMMETRY. 

66 NL = NDF*NEL 

DO 68 J=2,NL 

DO 68 K=1,J 

68 S(J,K) = S(K,J) 

ELM 191 

ELM 192 

ELM 193 

ELM 194 

ELM 195 

ELM 196 

ELM 197 

ELM 198 

ELM 199 

ELM 200 

ELM 201 

ELM 202 

ELM 203 

ELM 204 

ELM 205 

ELM 206 

ELM 207 

ELM 208 

ELM 209 

ELM 210 

ELM 211 

ELM 212 

ELM 213 

ELM 214 

ELM 215 

ELM 216 

ELM 217 

ELM 218 

ELM 219 

ELM 220 

ELM 221 

ELM 222 

ELM 223 

ELM 224 

ELM 225 

ELM 226 

ELM 227 

ELM 228 

ELM 229 

ELM 230 

ELM 231 

ELM 232 

ELM 233 

ELM 234 

ELM 235 

ELM 236 

ELM 237 

ELM 238 

ELM 239 

ELM 240 
APPENDIX C: NEW SUBROUTINES 

RETURN 

C. • • • COMPUTE STRESSES AND VELOCITY GRADIENTS 

4 RETURN 

C. • • • COMPUTE MASS MATRIX 

5 RETURN 

C •••• FORMATS 

1000 FORMAT(10F10.0) 

2000 FORMAT(5X,32HTWO DIMENSIONAL ICE FLOW ELEMENT// 

1 10X,12HVISCOS. A =,E12.5/ 

2 10X,12HVISCOS. N =,E12.5/ 

3 10X,12HCALVE NODE =,E12.5/ 

4 10X,12HNODES DEEP =,E12.5/) 

END 

SUBROUTINE PLT(ID,X,B,F,T,D,TIME,NDM,NDF) 

c 

C. • • • PLOT DISPLACEMENT VECTORS B AT NODAL COORDINATES X 

c 

LOGICAL PCOMP 

COMMON /PRLOD/PROP 

COMMON /CDATA/ O,HEAD(20),NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR 

COMMON /NORIG/ NUMEXT 

DIMENSION X(NDM,1),B(1),UL(6),ID(NDF,1),F(NDF,1),WD(2),D(1),T(1) 

DATA BL/4HBLAN/,DOTFAC/0.3/,TICFAC/2./,NWD/2/ 

DATA WD/4HSCAL,4HEND /,FIRST/0./ 

C •••• FUNCTIONS TO APPLY ORIGIN SHIFT, SCALE, AND ROTATION 

X2(XP,YP) = (SHIFT+ COSINE*(XP-XMIN) - SINE*(YP-YMIN))/SCALE 

Y2(XP,YP) = (SINE*(XP-XMIN) + COSINE*(YP-YMIN))/SCALE 

C. • • • FIRST=O MEANS READ IN NEW SCALES AS THIS IS THE FIRST 

C •••• TIME THROUGH THIS PLOT MACRO 

IF(ABS(FIRST).GT •• 5)GO TO 2 

FIRST = 1. 

10 READ(5,5000)CC 

DO 20 J=1,NWD 

20 IF(PCOMP(CC,WD(J)))GO TO 30 

GO TO 10 

30 GO TO (1,2),J 

C •••• READ SCALE INFORMATION 

1 READ(5,5001)XMIN,XMAX,YMIN,YMAX,SCALE,VSCALE,AHLEN,ROT 

DOTLEN = DOTFAC*AHLEN 

BLOWUP = SCALE/VSCALE 

ANGLE = 0. 

COSINE = 1. 

SINE = 0. 

SHIFT = 0. 

IF(ABS(ROT-90.).GT •• 5)GO TO 15 

COSINE = 0. 

SINE = 1. 

ANGLE = 90. 

SHIFT = YMAX - YMIN 

15 CONTINUE 

ELM 241 

ELM 242 

ELM 243 

ELM 244 

ELM 245 

ELM 246 

ELM 247 

ELM 248 

ELM 249 

ELM 250 

ELM 251 

ELM 252 

ELM 253 

PLT 001 

PLT 002 

PLT 003 

PLT 004 

PLT 005 

PLT 006 

PLT 007 

PLT 008 

PLT 009 

PLT 010 

PLT 011 

PLT 012 

PLT 013 

PLT 014 

PLT 015 

PLT 016 

PLT 017 

PLT 018 

PLT 019 

PLT 020 

PLT 021 

PLT 022 

PLT 023 

PLT 024 

PLT 025 

PLT 026 

PLT 027 

PLT 028 

PLT 029 

PLT 030 

PLT 031 

PLT 032 

PLT 033 

PLT 034 

PLT 035 

PLT 036 

PLT 037 

B39 
B40 

CALL PLOTS(0,0,99) 

GO TO 10 

STUDIES OF COLUMBIA GLACIER, ALASKA 

C.... PLOT DISPLACEMENT VECTORS 

2 CALL NEWPEN ( 2) 

DO 202 N=1,NUMNP 

IF(PCOMP(X(1,N),BL))GO TO 202 

DO 200 I=1,NDF 

UL(I) = F(I,N)*PROP*BLOWUP 

K = IABS(ID(I,N)) 

200 IF(K.GT.O)UL(I)=B(K)*BLOWUP 

XTAIL = X(1,N) 

YTAIL = X(2,N) 

IF(XTAIL.LT.XMIN.OR.XTAIL.GT.XMAX.OR.YTAIL.LT.YMIN.OR. 

* YTAIL.GT.YMAX)GO TO 202 

CALL SYMBOL(X2(XTAIL,YTAIL),Y2(XTAIL,YTAIL),DOTLEN,1,ANGLE,-1) 

DV = SQRT(UL(1)**2 + UL(2)**2)/SCALE 

IF(DV.LT.DOTLEN)GO TO 202 

AHLEN2 = AMINI (AHLEN,DV) 

XTIP = XTAIL + UL(1) 

YTIP = YTAIL + UL(2) 

IF(XTIP.LT.XMIN.OR.XTIP.GT.XMAX.OR.YTIP.LT.YMIN.OR. 

* YTIP.GT.YMAX)GO TO 202 

CALL AROHD(X2(XTAIL,YTAIL),Y2(XTAIL,YTAIL),X2(XTIP,YTIP), 

1 Y2(XTIP,YTIP),AHLEN2,0.,14) 

202 CONTINUE 

NDEEP = D(4) + 0.1 

CALL NEWPEN ( 1) 

C •••• PLOT BED, SMOOTHED. 

TICK = TICFAC*AHLEN*SCALE 

XR = X2(T(1),T(2)) 

YR = Y2(T(1),T(2)) 

CALL PLOT(XR,YR,3) 

CALL SMOOT(XR,YR,O) 

NDEEP1 = NDEEP + 1 

NUMEXM = NUMEXT - NDEEP 

NUMNPM = NUMNP - NDEEP 

DO 330 I=NDEEP1,NUMEXM,NDEEP 

XR = X2(T(I),T(I+1)) 

YR = Y2(T(I),T(I+1)) 

CALL SMOOT ( XR, YR,-2) 

XP = T (I-NDEEP) 

YP = T(I-NDEEP+1) - TICK 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL SMOOT(XR,YR,2) 

YP = YP + TICK 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

330.CALL SMOOT(XR,YR,2) 

XP = T(NUMEXM+1) 

PLT 038 

PLT 039 

PLT 040 

PLT 041 

PLT 042 

PLT 043 

PLT 044 

PLT 045 

PLT 046 

PLT 047 

PLT 048 

PLT 049 

PLT 050 

PLT 051 

PLT 052 

PLT 053 

PLT 054 

PLT 055 

PLT 056 

PLT 057 

PLT 058 

PLT 059 

PLT 060 

PLT 061 

PLT 062 

PLT 063 

PLT 064 

PLT 065 

PLT 066 

PLT 067 

PLT 068 

PLT 069 

PLT 070 

PLT 071 

PLT 072 

PLT 073 

PLT 074 

PLT 075 

PLT 076 

PLT 077 

PLT 078 

PLT 079 

PLT 080 

PLT 081 

PLT 082 

PLT 083 

PLT 084 

PLT 085 

PLT 086 

PLT 087 
YP = T(NUMEXM+2) 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL SMOOT(XR,YR,-24) 

YP = YP - TICK 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,2) 

YP = YP + TICK 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,2) 

K = NUMEXM 

- NDEEP + 1 

XP = T(K) 

YP = T(K+1) 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,3) 

YP = YP - TICK 

XR := X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,2) 

YP = YP + TICK 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,2) 

APPENDIX C: NEW SUBROUTINES 

C. • • • PLOT SURFACE, SMOOTHED. 

XR = X2(X(1,1),X(2,1)) 

YR = Y2(X(1,1),X(2,1)) 

CALL PLOT(XR,YR,3) 

CALL SMOOT(XR,YR,O) 

DO 340 I=NDEEP1,NUMNPM,NDEEP 

XR = X2(X(1,I),X(2,I)) 

YR = Y2(X(1,I),X(2,1)) 

340 CALL SMOOT (XR,YR,-2) 

XR = X2(X(1,NUMNPM+1),X(2,NUMNPM+1)) 

YR = Y2(X(1,NUMNPM+1),X(2,NUMNPM+1)) 

CALL SMOOT(XR,YR,-24) 

C.... PLOT CALVING FACE 

XP = X(1,NUMNPM+1) 

YP = X(2,NUMNPM+1) 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,3) 

XP = X(1,NUMNP) 

YP = X(2,NUMNP) 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,2) 

C. • • • PLOT SEALEVEL. 

PLT 088 

PLT 089 

PLT 090 

PLT 091 

PLT 092 

PLT 093 

PLT 094 

PLT 095 

PLT 096 

PLT 097 

PLT 098 

PLT 099 

PLT 100 

PLT 101 

PLT 102 

PLT 103 

PLT 104 

PLT 105 

PLT 106 

PLT 107 

PLT 108 

PLT 109 

PLT 110 

PLT 111 

PLT 112 

PLT 113 

PLT 114 

PLT 115 

PLT 116 

PLT 117 

PLT 118 

PLT 119 

PLT 120 

PLT 121 

PLT 122 

PLT 123 

PLT 124 

PLT 125 

PLT 126 

PLT 127 

PLT 128 

PLT 129 

PLT 130 

PLT 131 

PLT 132 

PLT 133 

PLT 134 

PLT 135 

PLT 136 

PLT 137 

B41 
B42 

XP = X (1,NUMEXT) 

yp = o. 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,3) 

XP = X(1,1) 

yp = o. 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL PLOT(XR,YR,2) 

STUDIES OF COLUMBIA GLACIER, ALASKA 

C •••• PLOT TIME. 

c 

K = NUMNPM/2 - NDEEP + 1 

XP = T(K) 

YP = T(K+1) - TICK*3. 

XR = X2(XP,YP) 

YR = Y2(XP,YP) 

CALL NUMBER(XR,YR,3.*AHLEN,TIME+1900.,ANGLE,3) 

CALL PLOT(X2(XMAX,YMIN),0.,-3) 

RETURN 

5000 FORMAT(A4) 

5001 FORMAT(8F10.0) 

END 

SUBROUTINE REST(ID,X,B,F,T,D,NDM,NDF,V1) 

C •••• READ NODE NUMBER N, COORDINATES X, AND DISPLACEMENTS B 

C •••• FROM RESTART TAPE15 

c 

COMMON /CDATA/ O,HEAD(20),NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR 

COMMON /RSTART/ !START 

COMMON /NORIG/ NUMEXT 

DIMENSION X(NDM,1),B(1),UL(6),ID(NDF,1),XR(3),T(1),F(NDF,1),D(1) 

REWIND 15 

C •••• SHIFT "FIXED" BED COORDINATES, BOOT, AND "A" TO TEMPERATURE ARRAY 

NDEEP = D(4) + 0.1 

DO 5 N=NDEEP I NUMNP I NDEEP 

DO 4 I=1,NDM 

4 T(N-NDEEP+I) = X(I,N) 

T(N-NDEEP+4) = F(3,N-NDEEP+4) 

5 T(N-1) = F(3,N-1) 

C. • • • SET TOTAL NUMBER OF NODES IN EXTENDED GRID, BEFORE 

C •••• NUMNP IS REDUCED IN SUBROUTINE EIMT04. 

NUMEXT = NUMNP 

NVBLS = NDM + NDF 

10 IF(NVBLS.LE.4)READ(15,15000)N,(XR(I),I=1,NDM),(UL(I),I=1,NDF) 

IF(NVBLS.GT.4)READ(15,15001)N,(XR(I),I=1,NDM),(UL(I),I=1,NDF) 

C •••• CHECK FOR BLANK EOF INDICATOR 

IF(N.LE.O.OR.N.GT.NUMNP)GO TO 103 

C •••• SKIP REDEFINITION OF CURRENT COORDINATES, BOOT, ABED UNLESS 

C •••• V1 IS POSITIVE. 

IF(V1.LE.O.)GO TO 95 

PLT 138 

PLT 139 

PLT 140 

PLT 141 

PLT 142 

PLT 143 

PLT 144 

PLT 145 

PLT 146 

PLT 147 

PLT 148 

PLT 149 

PLT 150 

PLT 151 

PLT 152 

PLT 153 

PLT 154 

PLT 155 

PLT 156 

PLT 157 

PLT 158 

PLT 159 

RES 01 

RES 02 

RES 03 

RES 04 

RES 05 

RES 06 

RES 07 

RES 08 

RES 09 

RES 10 

RES 11 

RES 12 

RES 13 

RES 14 

RES 15 

RES 16 

RES 17 

RES 18 

RES 19 

RES 20 

RES 21 

RES 22 

RES 23 

RES 24 

RES 25 

RES 26 

RES 27 

RES 28 
APPENDIX C: NEW SUBROUTINES 

DO 90 I=1,NDM 

90 X(I,N) = XR(I) 

IF(MOD(N,NDEEP).EQ.4.0R.MOD(N,NDEEP).EQ.NDEEP-1)F(3,N) = UL(3) 

95 DO 100 I=1,NDF 

K = IABS(ID(I,N)) 

100 IF(K.GT.O)B(K)=UL(I) 

GO TO 10 

103 REWIND 15 

C •••• SET FLAG SO STRAIN-DEPENDENT VISCOSITY USED 

ISTART = NUMEL + 1 

RETURN 

15000 FORMAT(I16,4(1PE16.9)) 

15001 FORMAT(I16,4(1PE16.9)/5(1PE16.9)) 

END 

SUBROUTINE SAVE(ID,X,B,F,NDM,NDF) 

c 

C •••• SAVE NODE NUMBER N, COORDINATES X, AND DIS:PLACEMENTS B 

C •••• ON TAPE16. 

c 

LOGICAL PCOMP 

COMMON /PRLOD/PROP 

COMMON /CDATA/ O,HEAD(20),NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR 

DIMENSION X(NDM,1),B(I),UL(6),ID(NDF,1),F(NDF,1) 

DATA BL/4HBLAN/ 

REWIND 16 

NVBLS = NDM + NDF 

DO 102 N=1,NUMNP 

IF(PCOMP(X(1,N),BL))GO TO 102 

DO 100 I=1,NDF 

UL(I) = F(I,N)*PROP 

K = IABS(ID(I,N)) 

100 IF(K.GT.O)UL(I) = B(K) 

IF(NVBLS.LE.4)WRITE(16,16000)N,(X(I,N),I=1,NDM),(UL(I),I=1,NDF) 

IF(NVBLS.GT.4)WRITE(16,16002)N,(X(I,N),I=1,NDM),(UL(I),I=1,NDF) 

102 CONTINUE 

C.... WRITE END OF DATA INDICATOR -- BLANK CARD 

WRITE(16,16001) 

IF(NVBLS.GT.4)WRITE(16,16001) 

16001 FORMAT(1H ) 

END FILE 16 

REWIND 16 

RETURN 

16000 FORMAT(I16,4(1PE16.9)) 

16002 FORMAT(I16,4(1PE16.9)/5(1PE16.9)) 

END 

SUBROUTINE UPNODE ( D I X I B ~ DT IF IT I TIME I ID., NOM, NDF) 

C •••• UPDATE MESH FOR 2-D DYNAMIC GLACIER PROBLEM 

COMMON /CDATA/O,HEAD(20),NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR 

COMMON /NORIG/ NUMEXT 

DIMENSION D(1),X(NDM,1),B(1),ID(NDF,1),FACT(3),T(1),F(NDF,1) 

RES 29 

RES 30 

RES 31 

RES 32 

RES 33 

RES 34 

RES 35 

RES 36 

RES 37 

RES 38 

RES 39 

RES 40 

RES 41 

RES 42 

SAV 01 

SAV 02 

SAV 03 

SAV 04 

SAV 05 

SAV 06 

SAV 07 

SAV 08 

SAV 09 

SAV 10 

SAV 11 

SAV 12 

SAV 13 

SAV 14 

SAV 15 

SAV 16 

SAV 17 

SAV 18 

SAV 19 

SAV 20 

SAV 21 

SAV 22 

SAV 23 

SAV 24 

SAV 25 

SAV 26 

SAV 27 

SAV 28 

SAV 29 

SAV 30 

SAV 31 

UPN 001 

UPN 002 

UPN 003 

UPN 004 

UPN 005 

B43 
B44 STUDIES OF COLUMBIA GLACIER, ALASKA 

DIMENSION DKNIK(40),IXS(3),IXB(3),IX(9),XDUM(1,5),SHP(3,9) 

DIMENSION ISURF(3) 

DATA FACT/0.4,0.8,0.9/ 

DATA RHOW/1./,RHOI/0.9/,ACONST/1165625./,BCONST/0.550989/ 

DATA CCONST/-2.23920/,IX/1,1,0,0,1,0,0,0,0/,TT/-1./ 

DATA XDUM/5*0./,ISURF/1,2,5/,SHP/27*0./ 

DATA DKNIK/24.27, 23.69, 25.38, 28.15, 29.78, 29.36, 27.44, 

* 21.51, 22.95, 23.86, 25.43, 26.12, 44.59, 84.95, 99.69, 143.10, 

* 220.03, 317.97, 426.72, 479.86, 542.15, 607.59, 656.56, 720.51, 

* 805.18, 705.21, 476.12, 356.48, 372.13, 257.24, 157.36, 135.86, 

* 83.29, 61.79, 62.53, 44.41, 32.71, 30.47, 27.13, 25.63/ 

NDEEP = D(4) + 0.1 

NDEEP1 = NDEEP + 1 

NDPM4 = NDEEP - 4 

NDPM5 = NDEEP - 5 

C. • • • FIRST DETERMINE UNSUPPORTED HEIGHT HU. 

C. • • • DO BEFORE PROJECTING SURFACE TO (TIME+DT) • 

HG = X(2,NUMNP-NDEEP+1) 

HW = AMAX1 (0. ,-X (2,NUMNP)) 

HU = HG - ( RHOW/RHOI - 1. ) *HW 

C.... LOCATE ADJUSTED SURFACE ALTITUDE. 

DO 10 I=NDEEP1 ,NUMNP ,NDEEP 

K = IABS(ID(3,I+1)) 

10 X(2,I) = X(2,I) + B(K)*DT 

C •••• DETERMINE DELX. 

C •••• FIRST DETERMINE V. 

K = IABS(ID(1,NUMNP-NDEEP+1)) 

V = B(K) 

C •••• NEXT DETERMINE DISCHARGE "DISCH" FOR THE TIME INTERVAL "TIME" TO 

C •••• "TIME+ DT". 

C •••• NEXT DETERMINE DISCHARGE FROM THE SAME INTERVAL. 

I = AMOD(TIME,1.0)/DT + 1.1 

DISCH = QKNIK(I) 

C •••• DETERMINE CALVING SPEED USING CALVING LAW. 

VC = ACONST*DISCH**BCONST*HU**CCONST 

C •••• DETERMINE DELX 

DELX = (V- VC)*DT 

C •••• FIND ORIGINAL TOTAL GLACIER LENGTH. 

XTORIG = T(NUMNP-NDEEP+1) 

XOORIG = T(1) 

ORIGL = XTORIG - XOORIG 

C. • • • FIND THE BED LENGTH. 

XBED = T(NUMEXT-NDEEP+1) 

BEDLEN = XBED - XOORIG 

C. • • • FIND NEW CURRENT TOTAL LENGTH. 

CURRL = X (1,NUMNP-NDEEP+1) + DELX 

C •••• INSURE THAT TERMINUS REMAINS WITHIN THE BED LENGTH. 

CURRL = AMIN1(BEDLEN,CURRL) 

CURRL = AMAX1(CURRL,O.) 

C •••• FIND NEW COLUMN X-COORDINATES -- PROPORTIONED 

UPN 006 

UPN 007 

UPN 008 

UPN 009 

UPN 010 

UPN 011 

UPN 012 

UPN 013 

UPN 014 

UPN 015 

UPN 016 

UPN 017 

UPN 018 

UPN 019 

UPN 020 

UPN 021 

UPN 022 

UPN 023 

UPN 024 

UPN 025 

UPN 026 

UPN 027 

UPN 028 

UPN 029 

UPN 030 

UPN 031 

UPN 032 

UPN 033 

UPN 034 

UPN 035 

UPN 036 

UPN 037 

UPN 038 

UPN 039 

UPN 040 

UPN 041 

UPN 042 

UPN 043 

UPN 044 

UPN 045 

UPN 046 

UPN 047 

UPN 048 

UPN 049 

UPN 050 

UPN 051 

UPN 052 

UPN 053 

UPN 054 

UPN 055 
APPENDIX C: NEW SUBROUTINES 

C. • • • AS IN ORIGINAL MESH. 

C. • • • THE FIRST COLUMN REMAINS FIXED -- START WITH THE SOCOND. 

IXS(1) = 1 

IXS(3) = 1 + NDEEP 

IXS(2) = 1 + 2*NDEEP 

IXB(1) =· 1. 

IXB(3) = 1 +·NDEEP 

IXB(2) = 1 + 2*NDEEP 

DO 300 I=NDEEP1 ,NUMNP ,NDEEP 

C •••• FIND ORIGINAL X-COORDINATE OF COLUMN. 

XIORIG = T(I) 

C •••• FIND PROPORTIONATE NEW X-GOORDINATE. 

XSURF = (XIORIG - XOORIG)IORIGL*CURRL + XOORIG 

C •••• REPLACE X-COORDINATE OF ALL NODES IN THIS COLUMN EXCEPT 

C •••• THE SURFACE NODE. 

DO 40 J=2,NDEEP 

40 X(1,I+J-1) = XSURF 

C~ .... DETERMINE CORRECT ELEMENT FOR SURFACE ALTITUDE INTERPOLATION. 

50 IF(IXS(2).GE.NUMNP-NDEEP+l)GO TO 60 

IF(XSURF.LE.X(1,IXS(2)))GO 'ID 60 

DO 55 K=1,3 

55 IXS'{K) = IXS (K) + NDEEP 

GO TO. 50 

60 CONTINUE · 

C. • • • INTERPOLATE SURFACE. 

C •••• FIRST FIND CORRESPONDING NATURAL COORDINATES. 

IF(XSURF.GE.X(1,IXS(3))) 

* ss = (XSURF - X(1,IXS(3·))) I (X(1,IXS(2)) - X(1,IXS(3))) 

IF(XSURF.LT.X(1,IXS(3))) 

* ss = (XSURF- X(1,IXS(3))J I (X(1,IXS(3)) - X(1,IXS(1))) 

CALL SHAPE(SS,TT,XDUM,SHP;XSJ,1,5,IX,.TRUE.) 

YSURF = 0. 

DO 70 K=1,3 

70 YSURF = YSURF + SHP(3,ISURF(K))*X(2,IXS(K))' 

C •••• STORE INTERPOLATED NEW CURRENT Y-COORDINATE OF THE SURFACE NODE 

C •••• IN PLACE OF THEY-COORDINATE OF THE SECOND NODE OF THE COLUMN, 

C. • • • TEMPORARILY. 

X(2,I+1) = YSURF 

C. • • • DETERMINE CORRECT ELEMENT FOR BED INTERPOLATION. 

250 IF(IXB(2).GE.NUMEXT-NDEEP+1)GO TO 260 

IF(XSURF.LE.T(IXB(2)))GO TO 260 

DO 255 K=1,3 

255 IXB(K) = IXB(K) + NDEEP 

GO TO 250 

260 CONTINUE 

C. • • • INTERPOLATE BED. 

C ..... FIRST FIND. CORRESPONDING NATURAL COORDINATES. 

IF(XSURF.GE.T(IXB(3))) 

* ss = (XSURF - T(IXB(3))) I (T'(IXB(2)) - T(IXB(3))) 

IF(XSURF.LT.T(IXB(3))) 

UPN 056 

UPN 057 

UPN 058 

UPN 059: 

UPN 060 

UPN 061 

UPN 062 

UPN 063 

UPN 064 

UPN 065 

UPN 066 

UPN 067 

UPN 068 

UPN 069 

UPN 070 

UPN 071 

UPN 072 

UPN 073 

UPN 074 

UPN 075 

UPN 076 

UPN 077 

UPN 078 

UPN 079 

UPN 080 

UPN 081: 

UPN 082 

UPN 083 

UPN 084 

UPN 085 

UPN 086 

UPN 0.87 

UPN 088 

UPN 089 

UPN 090 

UPN 091 

UPN 092 

UPN 093 

UPN 094 

UPN 095 

UPN 096 

UPN 097 

UPN 098 

UPN 099 

UPN' 100 

UPN' 101 

UPN' 102 

UPNI 103 

UPNI 104 

UPN 105 

B45 
B46 STUDIES OF COLUMBIA GLACIER, ALASKA 

* ss = (XSURF- T(IXB(3))) I (T(IXB(3)) - T(IXB(1))) 

CALL SHAPE(SS,TT,XDUM,SHP,XSJ,1,5,IX,.TRUE.) 

YBED = 0. 

BCURR = 0. 

ACURR = 0. 

DO 270 K=1,3 

YBED = YBED + SHP(3,ISURF(K))*T(IXB(K)+1) 

BCURR = BCURR + SHP(3,ISURF(K))*T(IXB(K)+3) 

270 ACURR = ACURR + SHP(3,ISURF(K))*T(IXB(K)+NDEEP-2) 

C •••• REPLACE OLD "CURRENT" BED ALTITUDE, BOOT, AND ABED WITH 

C •••• NEW ONES (CORRESPONDING TO XSURF NOW). 

X(2,I+NDEEP-1) = YBED 

X(2,I+NDEEP-2) = YBED + 5. 

X(2,I+NDEEP-3) = YBED + 10. 

X(2,I+NDEEP-4) = YBED + 10. 

F(3,I+3) = BCURR 

F(3,I+NDEEP-2) = ACURR 

300 CONTINUE 

C •••• REPLACE COLUMN COORDINATES WITH NEW INTERPOLATED ONES. 

DO 310 I=NDEEP1 ,NUMNP ,NDEEP 

X(1,I) = X(1,I+1) 

X(2,I) = X(2,I+1) 

DELTA = X(2,I) - X(2,I+NDPM4) 

DO 310 K=1,NDPM5 

310 X(2,I+K) = X(2,I) - DELTA*FACT(K) 

RETURN 

END 

SUBROUTINE VALUES(VALUE,F,I,J,NDF) 

C ••••• LOCATE BOUNDARY VALUES FOR SUBROUTINE ELMT04 IN COMMON M. 

DIMENSION F(NDF,1) 

VALUE = F(I,J) 

RETURN 

END 

UPN 106 

UPN 107 

UPN 108 

UPN 109 

UPN 110 

UPN 111 

UPN 112 

UPN 113 

UPN 114 

UPN 115 

UPN 116 

UPN 117 

UPN 118 

UPN 119 

UPN 120 

UPN 121 

UPN 122 

UPN 123 

UPN 124 

UPN 125 

UPN 126 

UPN 127 

UPN 128 

UPN 129 

UPN 130 

UPN 131 

UPN 132 

VAL 1 

VAL 2 

VAL 3 

VAL 4 

VAL 5 

VAL 6 
APPENDIX D: MODIFICATIONS OF THE 

GENERAL PURPOSE FINITE ELEMENT PROGRAM 

Included below are modifications to the general pur- through 80 of the card. Changes with blanks in these colpose finite element program as it is listed in Zienkiewicz umns are attributable to the author. Enough lines of the 

(1977). The subroutine FREEFM is from R. L. Taylor general purpose subroutines of Zienkiewicz (1977) are 

(written commun., 1980), as received from him in the given so that the reader can locate the changes. The incard deck of the general purpose program. It is included elusion of modifications by Taylor is meant only to 

here for completeness, so that the reader may refer to specify the exact form of the program used in the modelthe program listed in Zienkiewicz. Some single line ing. It is strongly recommended that the interested 

modifications listed below are also from Taylor; these reader obtain the latest version of the program directly 

are indicated either by explicit statement to that effect from Taylor. 

or by identification of the form CR 4/80C in columns 73 

C**** COMMENT OUT CDC PROGRAM CARD NOT REQUIRED ON IBM. 

C PROGRAM FEAPB(INPUT=101,0UTPUT=101,TAPE4=INPUT,TAPE5=101, 

C 1 TAPE6=0UTPUT) 

C MASTER MINIFEM 

C. • • • SET PROGRAM CAPACITY * MAX MUST AGREE WITH DIMENSION OF M 

C**** CHANGE PROGRAM CAPACITY FOR PROBLEM AT HAND. 

COMMON M(20000) 

COMMON/PSIZE/ MAX 

C**** CHANGE MAX TO AGREE WITH DIMENSION OF M. 

MAX ='20000 

CALL FREEFM(M,M(11),M(21),M(31),M(111)) 

CALL PCONTR 

STOP 

END 

SUBROUTINE FREEFM(IS,IE,IT,IN,IWD) 

C. • • • FREE FORMA TED READ * * CONVERTED TO FIELD WIDTHS OF TEN ( 10) 

C •••• CARDS WHICH BEGIN WITH A STAR (*) ARE ALPHAI,IDMERIC 

C •••• A CARD WHICH BEGINS WITH 'IWO STARS (**) WILL TERMINATE INPUT 

LOGICAL BL,FL 

DIMENSION IS(10),IE(10),IT(10),IN(80),IWD(80) 

C**** CHANGE 10H TO 01H TO ACCOMMODATE IBM (STILL OK ON CDC). 

C DATA NDIG,IBL/10,10H /,ISL/1H// 

DATA NDIG,IBL/10,01H /,ISL/1H// 

DATA IA,IC,IR,IZ/1HA,1H,,1HR,1HZ/,IST/1H*/ 

REWIND 5 . 

WRITE(6,2000) 

FL = .TRUE. 

10 READ(4,1000) IN 

ILM = 1 

DO 20 I = 1,80 

IF(IN(I).NE.IBL) ILM =I +1 

20 CONTINUE 

IL = 1 

IFEM 1C 

IFEM 2t:= 

IFEM 4C 

CR 4/80C 

IFEM 6C 

IFEM 7C 

IFEM 8C 

FRE 01 

FRE 02 

FRE 03 

FRE 04 

FRE 05 

FRE 06 

FRE 08 

FRE 09 

FRE 10 

FRE 11 

FRE 12 

FRE 13 

FRE 14 

FRE 15 

FRE 16 

FRE 17 

B47 
B48 

30 

35. 

40 

45 

50 

100 

105 

'110 

• 

120 

125 

140 

144 

145 

150 

STUDIES OF COLUMBIA GLACIER, ALASKA 

JL = IL + 1 

DO 35 I = 1,80 

IWD(I) = IBL 

IF(IN(IL).NE.IST) GO TO 50 

IF(IN·(JL) .;EQ. IST) GO TO 160 

DO 40 ·1 :: JL,ILM 

IFJI'N(I).EQ.ISL) GO TO 45 

IWD(I-IL) = IN(I) 

lL = ILM + 1 

GO TO 150 

IL = I + 1 

GO TO 150 

DO 100 I = 1,1.0 

IS(I) = 0 

IT(I) = IR 

BL ·-= • TRUE. 

K = 1 

J = IN(IL) 

IF(J.E:Q.IBL.AND.BL) GO TO 120 

IF(J .EQ. IBL.OR.J .EQ. IC.OR.J .EQ •. ISL) GO TO 110 

IF(~NOT.BL) GO TO 120 

BL = .FALSE. 

IS(K) = IL 

IF (J .GE. IA.AND.J..,LE. IZ) .I'll(K) = IA 

GO TO 120· 

BL = .TRUE. 

IE(K) = IL - 1 

K=.K+l 

IL·= IL + 1 

IF(K.GT.S.OR.J.EQ.ISL) GO TO 125 

IF(IL.LE.ILM) GO TO 105 

K = K - 1 

IF(K.EQ.O) GO TO 150 

KL = 1 

KR = NDIG 

DO 145 I = 1,K 

JL = IS(I) 

IF(JL.EQ.O) GO TO 144 

JR = IE(I) 

IF(JR-JL.GT.9) GO TO 200 

II = KR- JR 

IF(IT(I).EQ.IA) II= KL- JL 

DO 140 J = JL,JR 

IWD(J+II) = IN(J) 

KL = KR + 1 

KR = KR + NDIG 

WRITE(6,2001) IWD 

WRITE(5,1000) IWD 

IF(IL.LE.ILM.OR.IN(IL-l).EQ.ISL) GO TO 30 

GO TO 10 

REWIND 5 

.WRITE( 6, 2002) 

IF (FL) RETURN 

STOP 

FRE 18 

FRE 19 

FRE 20 

FRE 21 

FRE 22 

FRE 23 

FRE 24 

FRE 25 

:.FRE 26 

FRE 27 

:FRE 28 

FRE 29 

FRE 30 

FRE 31 

FRE 32 

FRE 33 

FRE 34 

FRE 35 

FRE 36 

FRE 37 

FRE 38 

FRE 39 

FRE 40 

FRE ,41 

FRE '42 

FRE 43 

FRE 44 

FRE 45 

FRE 46 

FRE 47 

FRE 48 

FRE 49 

FRE 50 

FRE 51 

FRE 52 

FRE 53 

FRE 54 

FRE 55 

FRE 56 

FRE 57 

FRE 58 

FRE 59 

FRE 60 

FRE 61 

FRE 62 

FRE 63 

FRE 64 

FRE 65 

FRE 66 

FRE 67 

FRE 68 

FRE 69 

FRE 70 

FRE 71 
APPENDIX D: MODIFICATIONS, GENERAL PURPOSE FINITE ELEMENT PROGRAM B49 

200 WRITE(6,3000) IN 

FL = .FALSE. 

GO TO 10 

1000 FORMAT(80A1) 

2000 FORMAT(1H1,54HFIXED FORMAT FORM OF DATA CARDS FROM FREE FORMAT 

1UT ///1X) 

2001 FORMAT(5X,80A1) 

2002 FORMAT(1H0,17HEND OF DATA INPUT ///1X) 

3000 FORMAT(5X,48HMORE THAN TEN CHARACTERS/DIGITS ON THE DATA CARD// 

1 10X,80A1) 

END 

SUBROUTINE PCONTR 

COMMON M(1) 

C**** ADD THE FOLLOWING COMMON SO VARIABLES IN M CAN BE FOUND. 

COMMON /POINTR/ NO,N1,N2,N3,N4,NS,N6,N7,N8,N9,N10,N11,N12,N13, 

* N14,NE 

DDKENSION TITL(20),WD(3) 

CALL SE'IMEM (NE) 

CALL PZERO(M(N13),NEQ) 

GO TO 1 

C**** FROM R.L. TAYLOR -- NEW INPUT FORMAT. 

1001 FORMAT(8I10) 

END 

BLOCK DATA 

C  SINGLE PRECISION IS SUFFICIENT. 

DATA 0/1HO/,IPR/1/ 

END 

SUBROUTINE GENVEC(NDM,X,CD,PRT,ERR) 

C** ** FROM R.L. TAYLOR -- NEW INPUT FORMAT. 

1000 FORMAT(2I10,6F10.0) 

END 

FRE 72 

FRE 73 

FRE 74 

FRE 75 

INP FRE 76 

FRE 77 

FRE 78 

FRE 79 

FRE 80 

FRE 81 

FRE 82 

PCO 1C 

PC014C 

PCO 15C 

PCO 59C 

CR10/77C 

PCO 60C 

PCO 68C 

PCO 74C 

lC 

9C 

GEN 1C 

GEN 38C 

GEN 44C 
B50 STUDIES OF COLUMBIA GLACIER, ALASKA 

SUBROUTINE PMESH(IDL,IE,D,ID,X,IX,F,T,NDF,NDM,NEN1,III) PME 1C 

C  REMOVE DATA SPECIFICATION FOR UNUSED VARIABLE VA. 

DATA WD/4HCOOR, 4HELEM, 4HMATE, 4HBOUN, 4HFORC, 4HTEMP, 4HEND , 4HPRIN, PME 11C 

C 1 4HNOPR,4HPAGE/,BL/4HBLAN/,VA/4H VAL,2HUE/,LIST/10/,PRT/.TRUE./ 

1 4HNOPR,4HPAGE/,BL/4HBLAN/, LIST/10/,PRT/.TRUE./ 

C •••• INITIALIZE ARRAYS PME 13C 

ID(I,N) = 0 

X(1,N) = BL 

100 F(I,N) = 0.0 

C** ** FROM R. L. TAYLOR -- OMIT XHED FROM VARIABLES READ. 

READ(5,1002) MA,IEL 

C  FROM R. L. TAYLOR -- OMIT XHED FROM VARIABLES WRITTEN. 

WRITE(6,2003)MA,IEL 

IE(MA) = IEL 

300 CALL ELMLIB(D(1,MA),IDL,X,IX,IDL,IDL,IDL,NDF,NDM,NDF,1) 

GO TO 10 

DO 53 I = 1,NDF 

C** ** GENERATION FOR TIED DEGREES OF FREEDOM. 

IF(ID(I,L-LG).GT.1)ID(I,L)=ID(I,L-LG) + LG 

53 IF(ID(I,L-LG).LT.O) ID(I,L) = -1 

C  FROM R.L. TAYLOR -- NEW INPUT FORMAT. 

1001 FORMAT(8I10) 

C** ** FROM R.L. TAYLOR -- NEW INPUT FORMAT. 

1002 FORMAT(2I10) 

END 

SUBROUTINE PMACR (UL,XL,TL,LD,P,S,IE,D,ID,X,IX,F,T,JDIAG,B,DR,CT 

1,NDF,NDM,NEN1,NST,NEND) 

COMMON /TDATA/ TIME,DT,C1,C2,C3,C4,C5 

C* *** INCREASE THE NUMBER OF MACROS. 

DIMENSION WD(27),CT(4,1),CTL(4),LVS(9),LVE(9),JDIAG(1), 

1 UL(1),XL(1),TL(1),LD(1),P(1),S(1),IE(1),D(1),ID(1),X(1), 

2 IX(1),F(1),T(1),B(1),DR(1) 

PME 18C 

CR 5/79C 

PME 19C 

PME 60C 

PME 61C 

PME 62C 

CR 8/77C 

PME 64C . 

PME SOC 

PME 81C 

PME105C 

PME106C 

PME117C 

PMA 1C 

PMA 2C 

PMA 14C 

PMA 16C 

PMA 17C 
APPENDIX D: MODIFICATIONS, GENERAL PURPOSE FINITE ELEMENT PROGRAM 

C* *** INCREASE THE NUMBER OF MACROS. 

DATA WD/4HTOL , 4HDT , 4HSTRE, 4HDISP, 4HTANG, 4HFORM, 4HLOOP, 4HNEXT, 

1 4HPROP,4HDATA,4HTIME,4HCONV,4HSOLV,4HLMAS,4HCMAS,4HMESH, 

2 4HEIGE,4HEXCD,4HUTAN,4HREAC,4HCHEC,4HSAVE,4HREST,4HPLOT, 

3 4HUTAF,4HTANF,4HUPNO/ 

DATA NWD/27/,ENDM/4HEND /,NV,NC/1,1/ 

C •••• SET INITIAL VALUES OF PARAMETERS 

IF(L.NE.1.AND.L.NE.LL) 

1WRITE(6,2010) I,(CT(K,L),K = 1,4) 

C**** INCREASE THE NUMBER OF MACROS. 

GO TO (1,2,3,4,5,6,7,S,9,10,11,12,13,14,15,16,17,1S,19,20,21, 

* 22,23,24,25,26,27),J 

C •••• SET SOLUTION TOLERANCE 

CALL PFORM(UL,XL,TL,LD,P,S,IE,D,ID,X,IX,F,T,JDIAG,DR,DR,DR, 

2 NDF,NDM,NEN1,NST,6,B,M(NV),.FALSE.,.TRUE.,.FALSE.,.FALSE.) 

C**** ADD STATEMENT NUMBER 60. 

60 BFR = .TRUE. 

RN = 0. 

21 CALL PFORM(UL,XL,TL,LD,P,S,IE,D,ID,X,IX,F,T,JDIAG,DR,DR,DR, 

1 NDF,NDM,NEN1,NST,2,B,F,.FALSE.,.FALSE.,.FALSE~,.FALSE.) 

C**** ADD MACROS. 

GO TO 330 

C •••• SAVE NODE NUMBER N, COORDINATES X, AND DISPLACEMENTS B 

C •••• ON TAPE16. 

22 CALL SAVE(ID,X,B,F,NDM,NDF) 

GO TO 330 

C •••• READ NODE NUMBER N, COORDINATES X (IF CT(3,L) = 1.), 

C •••• AND DISPLACEMENTS B FROM RESTART TAPE15. 

C •••• IF CT(3,L) = 0., READ ONLY DISPLACEMENTS. 

23 CALL REST(ID,X,B,F,T,D,NDM,NDF,CT(3,L)) 

GO TO 330 

C •••• PLOT 

24 CALL PLT(ID,X,B,F,T,D,TIME,NDM,NDF) 

GO TO 330 

C. • • • FORM TANGENT STIFFNESS AND OUT OF BALANCE FORCES, ONE CALL 

C. • • • PER ELEMENT. 

25 IF(CFL)CALL PSETM(NC,NE,JDIAG(NEQ)*IPR,CFL) 

CALL PZERO(M(NC),JDIAG(NEQ)) 

CFR = .TRUE. 

26 IF(J.EQ.26)CFR = .FALSE. 

IF(GFL)CALL PSETM(NA,NE,JDIAG(NEQ)*IPR,GFL) 

IF(NPLD.GT.O)PROP = PROPLD(TIME,O) 

CALL PZERO(M(NA),JDIAG(NEQ)) 

CALL PLOAD ( ID, F, DR, NNEQ, PROP) 

B51 

PMA 1SC 

PMA 19C 

PMA 22C 

PMA sse 

PMA S6C 

PMA sse 

PMA121C 

PMA122C 

PMA124C 

PMA229C 

PMA230C 
B52 STUDIES OF COLUMBIA GLACIER, ALASKA 

CALL PFORM(UL,XL,TL,LD,P,S,IE,D,ID,X,IX,F,T,JDIAG,DR,M(NA),M(NC), 

2 NDF,NDM,NEN1,NST,3,B,M(NV) ,.TRUE., .TRUE. ,CFR, .• FALSE.) 

AFR = .TRUE. 

GO TO 60 

C •••• UPDATE MESH 

27 CALL UPNODE(D,X,B,DT,F,T,TIME,ID,NDM,NDF) 

330 L = L + 1 

C**** FROM R.L. TAYLOR -- NEW INPUT FORMAT. 

1000 FORMAT(A4,6X,A4,6X,2F10.0) 

END 

SUBROUTINE ACTCOL(A,B,JDIAG,NEQ,AFAC,BACK) 

ID = JDIAG(K+I) 

C  GIVE ERROR MESSAGE FOR DIVISION BY 0. 

IF(A(ID).NE.O.O) GO TO 310 

INDEX = K + I 

GO TO 1100 

310 D = A(I) 

A(I) = A(I)/A(ID) 

ID = JDIAG(I) 

C** ** GIVE ERROR MESSAGE FOR DIVISION BY 0. 

IF(A(ID).NE.O.O)GO TO 610 

INDEX = I 

GO TO 1100 

610 B(I) = B(I)/A(ID) 

700 AENGY = AENGY + B(I)*B(I)*A(ID) 

GO TO 800 

C**** ERROR MESSAGE AND STOP FOR SINGULAR SYSTEM. 

110·0 WRITE(6,6000) INDEX 

6000 FORMAT(46HOZERO DIAGONAL IN PARTIALLY REDUCED SYSTEM AT , 

* 16HI.>IAGONAL NUMBER , 'I 6) 

STOP 

END 

SUBROUTINE PROFIL ( JD IAG, ID, IX, NDF, NEN1, NAD) 

GO TO 40· 

C**** RESTRAINED. (J .LE.1) OR TIED (J:.GT.1) NODES 

PMA231C 

PMA24lC 

PMA264C 

ACT 1C 

ACT 33C 

ACT 36C 

ltCT45C 

ACT 47C 

AC'l' 61C 

ACT 62C 

PRO lC 

PRO 16C 
30 

40 

APPENDIX D: MODIFICATIONS, GENERAL PURPOSE FINITE ELEMENT PROGRAM 

IF(J.LE.1)ID(I,N) = 0 

IF(J.GT.l)ID(I,N)=ID(I,MOD(J-1,NUMNP)+1) 

CONTINUE 

END 

SUBROUTINE PRTDIS(ID1 X,B,F,NDM,NDF) 

RETURN 

C  CHANGE TIME FORMAT FROM E13.5 TO E16.8. 

2000 FORMAT(A1,20A4//5X,19HNODAL D.ISPLACEMENTS,SX,4HTIME,E16.8// 

1 6X,4HNODE,9(I7,A4,A2)) 

C** ** 

END 

SUBROUTINE SHAP2(S,T,SHP,IX,NEL) 

T2 = (1.-T*T)/2. 

DO 100 I=5,9 

DO 100 J = 1,3 

IF(IX(9).EQ.O) GO To 107 

FROM R.L. TAYLOR -- MULTIPLY BY 

SHP(1,9) = -4.*S*T2 

SHP(2,9) = -4.*T*S2 

SHP(3,9) = 4.*S2*T2 

END 

4. 

B53 

PRO 18C 

PRO 49C 

PRT 1C 

PRT 25C 

PRT 27C 

PRT 28C 

SHA 1C 

SHA 7C 

CR 3/78C 

SHA 9C 

SHA 33C 

SHA 34C 

SHA 35C 

SHA 36C 

SHA 51C 
APPENDIX E: SAMPLE RUN 

The following sample run demonstrates input and output associated with going from time 1981.150 to time 

1981.275 by using calving law constant a= 1165625. The 

corresponding plots can be found in figure 4. During the 

actual modeling runs, 3 years were run at a time; the only modification necessary for such a multiple run is to 

change the macroinstruction LOOP*1 to LOOP*24. The 

asterisks appearing in some of the macros are merely 

spacers so that the subroutine FREEFM will place the 

constants that follow them in the proper card position. 

The input data contained in data set $$$$$$$.COLDATA is simply the freeform version of the data as 

listed at the beginning of the output, with single spaces 

between data items and an appended card with** in columns 1 and 2. In addition, the title card has an asterisk 

in column 1 to indicate alphanumeric data. The general 

purpose program is contained in subroutine library 

$$$$$$$.LOADLIB, with the element subroutine 

ELMT04 for glacier flow in element ICE2 of that 

library. The data sets $$$$$$$.0LDSTART and 

$$$$$$$.RESTART are restart files read by macro 

REST and written by macro SAVE, respectively. Data 

set $$$$$$$.COLPLOT is the plot file. 

SYSLRS232.DATA100.CALCOMP is a library of 

Calcomp1 software designed to drive a local Calcomp 

drum plotter. Step COPYPLOT copies the plot file to the 

appropriate output class K for local plotting. 

1The use of the brand name in this report is for identification purposes only and 

does not imply endorsement by the U.S. Geological Survey. 

H A S P J 0 B L 0 G 

$21.50.22 JOB 8770 -- $$$$$$$$ -- BEGINNING EXEC - INIT 7 - CLASS F 

$21.51.25 JOB 8770 ENDED 001105L 

---- HASP-II RE1 STATISTICS ----

4 0 CARDS READ 

1,116 SYSOUT PRINT RECORDS 

335 SYSOUT PUNCH RECORDS 

1.06MINUTES ELAPSED TIME 

//$$$$$$$$JOB ($$$$$$$$$,$$$$,$$,10,1000,,,N),'SIKONIA*R*', J8770 

II CLASS=F,MSGLEVEL=(2,0) 

***ROUTE PRINT REMOTE241 RE1 

***ROUTE PUNCH REMOTE21 REI 

//FEAPRUN EXEC FORTXCG,TIME=5,COND=(O,LT), 

II PARM.FORT='NOMAP,NAME(FEAP),OPTIMIZE(2),', 

II ULIB='$$$$$$$.LOADLIB',PARM.GO='NOMAP,EP=FEAP,SIZE=210K', 

II REGION.G0=236K 

//FORT.SYSPRINT DD SYSOUT=A 

//FORT.SYSIN DD * 

IEF-142I - STEP WAS EXECUTED - COND CODE 0000 

IEF373I STEP /FORT / START 81037.2150 

IEF374I STEP /FORT / STOP 81037.2150 CPU OMIN 00.17SEC MAIN 248K LCS OK 

CCD006I STEP /FORT / I/0 IN DD SEQ: 4A4:2 563:1 155:0 

CCD006I STEP /FORT / I/0 IN DD SEQ: 563:0 4Dl:1 

CCD001I STEP /FORT / MEMORY RESERVED 250K: ACCESSES: DISK 4/ TAPE 0 

//GO.SYSLIB DD 

II DD DSN=SYS1.RS232.DATA100.CALCOMP,DISP=SHR 

II DD DSN=SYS1.FORTGLIB,DISP=SHR 

II DD DSN=XTENT.LIB,DISP=SHR 

//GO.SYSLINB DD DSN=$$$$$$$.LOADLIB(ICE2),DISP=SHR 

//GO.FT06F001 DD SYSOUT=A 

//GO.SYSIN DD UNIT=SYSDK,SPACE=(TRK,(34,34)), 

II DCB=(RECFM=FB,LRECL=80,BLKSIZE=640) 

B54 
APPENDIX E: SAMPLE RUN 

//GO.FT04F001 DD DSN=$$$$$$$.COLDATA,DISP=SHR 

//GO.FT99F001 DD DSN=$$$$$$$.COLPLOT,DISP=OLD 

//GO.FT15F001 DD DSN=$$$$$$$.0LDSTART,DISP=SHR 

//GO.FT16F001 DD DSN=$$$$$$$.RESTART,DISP=(OLD,KEEP,KEEP) 

IEF142I - STEP WAS EXECUTED - COND CODE 0000 

B55 

IEF373I STEP /GO / START 81037.2150 

IEF374I STEP /GO / STOP 81037.2151 CPU OMIN 12.14SEC MAIN 234K LCS OK 

CCD006I STEP /GO / I/0 IN DD SEQ~ 4A3:1 634:260 314:232 

CCD006I STEP /GO / I/0 IN DD SEQ~ 314:0 314:0 563:3 

CCD006I STEP /GO / I/0 IN DD SEQ~ 634:10 155:89 4A4:46 

CCD006I STEP /GO / I/0 IN DD SEQ~ 4BO:O 4A5:0 746:42 

CCD006I STEP /GO / I/0 IN DD SEQ~ 740:42 741:4 746:102 

CCD001I STEP /GO / MEMORY RESERVED 236K~ ACCESSES: DISK 831/ TAPE 0 

//COPYPLOT EXEC PGM=IEBGENER,TIME=1,COND=(O,LT),REGION=104K 

//SYSPRINT DO SYSOUT=A 

//SYSIN DD DUMMY 

//SYSUT1 DD DSN=$$$$$$$.COLPLOT,DISP=SHR 

//SYSUT2 DD SYSOUT=(K,,WS81),DCB=(RECFM=U,BLKSIZE=80) 

IEF142I - STEP WAS EXECUTED - COND CODE 0000 

IEF373I STEP /COPYPLOT/ START 81037.2151 

IEF374I STEP /COPYPLOT/ STOP 81037.2151 CPU OMIN 00.22SEC MAIN 34K LCS OK 

CCD006I STEP /COPYPLOT/ I/0 IN DD SEQ~ 4A2:1 000:0 740:43 

CCD006I STEP /COPYPLOT/ I/0 IN DD SEQ~ 4B0:19 

CCD001I STEP /COPYPLOT/ MEMORY RESERVED 104K~ ACCESSES: DISK 63/ TAPE 0 

IEF375I JOB /$$$$$$$$/ START 81037.2150 

IEF376I JOB /$$$$$$$$/ STOP 81037.2151 CPU OMIN 12.53SEC 

CCD007A JOB /$$$$$$$$/ BASIC CHRG: FIXED CPU DSK-IO TPE-IO KBYTE-HR TOTAL 

.41 2.79 .so .54 4.54 

CCD007D JOB /$$$$$$$$/ LOCAL CHRG: FORMS PRINTING CRDS-RD CRDS-PU TOTAL 

.oo .00 

CCD002I JOB /$$$$$$$$/ TOTAL CHARGE RECORDED = $4.54 

CCD003I JOB /$$$$$$$$/ US GEOLOGICAL SURVEY RE1 OS REL 218 TIME = 21.51.25 DATE = 

02/06/81 

LEVEL 2.3.0 (JUNE 78) OS/360 FORTRAN H EXTENDED DATE 81.037/21.50.27 

REQUESTED OPTIONS: NOMAP,NAME(FEAP),OPTIMIZE(2), 

OPTIONS IN EFFECT: NAME(FEAP) OPTIMIZE(2) LINECOUNT(54) SIZE(MAX) AUTODBL(NONE) 

SOURCE EBCDIC NOLIST NODECK OBJECT NOMAP NOFORMAT GOSTMT NOXREF 

NOALC NOANSF NOTERM IBM 

C PROGRAM FEAPB(INPUT=101,0UTPUT=101,TAPE4=INPUT,TAPE5=101, 

C 1 TAPE6=0UTPUT) 

C MASTER MINIFEM IFEM 1C 

C. • • • SET PROGRAM CAPACITY * MAX MUST AGREE WITH 

DIMENSION OF M IFEM 2C 

ISN 0002 COMMON M(20000) 

ISN 0003 COMMON/PSIZE/ MAX IFEM 4C 

ISN 0004 MAX = 20000 

ISN 0005 CALL FREEFM(M,M(11),M(21),M(31),M(111)) 

ISN 0006 CALL PCONTR IFEM 6C 

ISN 0007 STOP IFEM 7C 

ISN 0008 END IFEM sc 
B56 STUDIES OF COLUMBIA GLACIER, ALASKA 

*OPTIONS IN EFFECT*NAME(FEAP) OPTIMIZE(2) LINECOUNT(54) SIZE(MAX) AUTODBL(NONE) 

*OPTIONS IN EFFECT*SOURCE EBCDIC NOLIST NODECK OBJECT NOMAP NOFORMAT GOST.MT NOXREF 

NOALC NOANSF NOTERM IBM 

*STATISTICS* SOURCE STATEMENTS = 7, PROGRAM SIZE = 318, 

SUBPROGRAM NAME = FEAP 

*STATISTICS* NO DIAGNOSTICS GENERATED 

****** END OF COMPILATION ****** 58K BYTES OF CORE NOT USED 

OS/360 LOADER 

OPTIONS USED - PRINT,NOMAP,NOLET,CALL,RES,NOTERM,SIZE=215040,NAME=**GO 

EP=FEAP 

TOTAL LENGTH 2EC28 

ENTRY ADDRESS 4BA810 

FIXED FORMAT FORM OF DATA CARDS FROM FREE FORMAT INPUT 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

232 36 1 2 3 9 0 

NOPR 

COOR 

1 1 o.o 465.00 

3 1 o.o -121.40 

5 0 o.o -268.00 

6 1 0.0 -268.00 

8 0 o.o -278.00 

9 1 762.50 430.00 

11 1 762.50 -103.60 

13 0 762.50 -237.00 

14 1 762.50 -237.00 

16 0 762.50 -247.00 

17 1 1525.00 411.00 

19 1 1525.00 -119.40 

21 0 1525.00 -252.00 

22 1 1525.00 -252.00 

24 0 1525.00 -262.00 

25 1 2287.50 401.00 

27 1 2287.50 -102.20 

29 0 2287.50 -228.00 

30 1 2287.50 -228.00 

32 0 2287.50 -238.00 

33 1 3050.00 370.00 

35 1 3050.00 -81.20 

37 0 3050.00 -194.00 

38 1 3050.00 -194.00 

40 0 3050.00 -204.00 

41 1 3812.50 340.00 

43 1 3812.50 -76.80 

45 0 3812.50 -181.00 

46 1 3812.50 -181.00 

48 0 3812.50 -191.00 

49 1 4575.00 314.00 

51 1 4575.00 -108.40 
APPENDIX E: SAMPLE RUN B57 

53 0 4575.00 -214.00 

54 1 4575.00 -214.00 

56 0 4575.00 -224.00 

57 1 5337.50 290.00 

59- 1 5337.5a -190.00 

61 0 5337.50 -310~~ 00 

62 1 5337.50 -310.00 

64 a 5337.50 -320.00 

65 1 6100.0.0 275.00 

67 1 6100.00 -242.60 

69 0 6100.00 -372.00 

70' 1 6100.00 -372.00· 

72' 0 6100.00 -382.00 

73 1 6862.50 264.00 

75 1 6862.50 -255.20 

77 0 6862.50 -385.00 

78 1 6862.50 -385.00 

80 0 6862.5'0 -395.00 

81 1 7625·. 00 248.00· 

83 1 762'5. 00 -249.60 

85 0 7625.00 -3.74.00 

86 1 7625.00 -374.00 

88 0 7625.00 -384.00 

89 1 8387.50 226.00 

91 1 8387.50 -253.20 

93 o· 8387.50 -373.00 

94 1 8387.50 -373.00 

96 0 8387.50 -3'83. 00 

97 1 9150.00 210.00 gg·~ 1 9150.00 -271.60 

101 0 9150.00 -392.00 

102 1 9150~00 -392.00' 

104 0 9150.00 -402.00 

105 1 9912.50 198.00 

107 1 9912.50 -273.20 

109 0 9912.50 -391.00 

110 1 9.912.50 ·-391. 00 

112 0 9912.50 -401~00 

113 1 10675.00 178.00 

115 1 10675.00 -278.00 

117 0 10675 .. 00 1 -392.00 

118· 1 10675.?0 -392.00 

120 0 10675.00 -402.00 

121 1 11056.25 167.00 

123 1 11056.25 -280.20. 

125 0 11056.25 -392.00 

126 1 11056.25 -392.00 

128 0 11056.25 -402:·. 00 

129 1 11437~50 159.00. 

131 1'· 11437.50 -265.80' 

13.3· 0 11437.50 -372.00 

134 1 11437.50 -372.00 

136 0 11437.50 -382.00 

137 1 11818.]5.· 152.00 
B58 STUDIES OF COLUMBIA GLACIER, ALASKA 

139 1 11818.75 -245.60 

141 0 11818.75 -345.00 

142 1 11818.75 -345.00 

144 0 11818.75 -355.00 

145 1 12200.00 140.00 

147 1 12200.00 -220.00 

149 0 12200.00 -310.00 

150 1 12200.00 -310.00 

152 0 12200.00 -320.00 

153 1 12581.25 129.00 

155 1 12581.25 -195.00 

157 0 12581.25 -276.00 

158 1 12581.25 -276.00 

160 0 12581.25 -286.00 

161 1 12962.50 118.00 

163 1 12962.50 -176.40 

165 0 12962.50 -250.00 

166 1 12962.50 -250.00 

168 0 12962.50 -260.00 

169 1 13153.12 113.00 

171 1 13153.12 -163.00 

173 0 13153.12 -232.00 

174 1 13153.12 -232.00 

176 0 13153.12 -242.00 

177 1 13343.75 108.00 

179 1 13343.75 -147.20 

181 0 13343.75 -211.00 

182 1 13343.75 -211.00 

184 0 13343.75 -221.00 

185 1 13534.37 100.00 

187 1 13534.37 -128.80 

189 0 13534.37 -186.00 

190 1 13534.37 -186.00 

192 0 13534.37 -196.00 

193 1 13725.00 96.00 

195 1 13725.00 -102.40 

197 0 13725.00 -152.00 

198 1 13725.00 

--152.00 

200 0 13725.00 -162.00 

201 1 13915.62 96.00 

203 1 13915.62 -80.80 

205 0 13915.62 -125.00 

206 1 13915.62 -125.00 

208 0 13915.62 -135.00 

209 1 14106.25 96.00 

211 1 14106.25 -57.60 

213 0 14106.25 -96.00 

214 1 14106.25 -96.00 

216 0 14106~25 -106.00 

217 1 14296.87 96.00 

219 1 14296.87 -20.80 

221 0 14296.87 -50.00 

222 1 14296.87 -50.00 

224 0 14296.87 -60.00 
APPENDIX E: SAMPLE RUN B59 

225 1 14487.50 96.00 

227 1 14487.50 19.20 

229 0 14487.50 0.0 

230 1 14487.50 o.o 

232 0 14487.50 -10.00 

ELEM 

01 1 3 19 17 1 11 18 

9 2 10 2 

03 1 8 24 22 6 16 23 

14 7 15 0 

04 1 19 35 33 17 27 34 

25 18 26 2 

06 1 24 40 38 22 32 39 

30 23 31 0 

07 1 35 51 49 33 43 so 

41 34 42 2 

09 1 40 56 54 38 48 55 

46 39 47 0 

10 1 51 67 65 49 59 66 

57 so 58 2 

12 1 56 72 70 54 64 71 

62 55 63 0 

13 1 67 83 81 65 75 82 

73 66 74 2 

15 1 72 88 86 70 80 87 

78 71 79 0 

16 1 83 99 97 81 91 98 

89 82 90 2 

18 1 88 104 102 86 96 103 

94 87 95 0 

19 1 99 115 113 97 107 114 

105 98 106 2 

21 1 104 120 118 102 112 119 

110 103 111 0 

22 1 115 131 129 113 123 130 

121 114 122 2 

24 1 120 136 134 118 128 135 

126 119 127 0 

25 1 131 147 145 129 139 146 

137 130 138 2 

27 1 136 152 150 134 144 151 

142 135 143 0 

28 1 147 163 161 145 155 162 

153 146 154 2 

30 1 152 168 166 150 160 167 

158 151 159 0 

31 1 163 179 177 161 171 178 

169 162 170 2 

33 1 168 184 182 166 176 183 

174 167 175 0 

34 1 179 195 193 177 187 194 

185 178 186 2 

36 1 184 200 198 182 192 199 

190 183 191 0 
B60 STUDIES OF COLUMBIA GLACIER, ALASKA 

MATE 

1 4 

.140 3.00 193. 8. 

BOUN 

8 8 -1 -1 

200 0 1 ·.1 

.J. 0 1 

2 0 1 

3 .o 1 

5 0 1 

6 0 1 5 

4 8 1 0 -1 

196 0 0 0 1 

7 8 1 0 -1 

199 0 0 0 1 

14 8 13 13 

198 0 197 197 

201 1 -1 -1 -1 

232 0 1 1 1 

FORC 

1 0 950.00 

2 0 943.56 

3 0 846.88 

4 0 784.83 ·o -2.67 

5 0 698.25 

6 0 698.25 

7 0 349.13 0 40. 

12 0 0 0 -3.42 

20 0 0 0 -3.98 

28 0 0 0 -4.45 

36 0 0 0 -4.90 

44 0 0 0 -5.41 

52 0 0 0 -5.98 

60 0 0 0 -6.30 

68 0 0 0 -6.32 

76 0 0 0 -6.29 

84 0 0 0 -6.26 

92 0 0 0 -6.27 

100 .0 0 0 -6.36 

108 0 0 0 -6.53 

116 0 0 0 -6.68 

124 0 0 0 -6.76 

132 0 0 0 -6.91 

140 0 0 0 -7.09 

148 0 0 0 -7.21 

156 0 0 0 ·-7.38 

.164 0 0 ·.0 -7.63 

17-2 0 ·o 0 -7.77 

.180 0 0 0 -7.87 

1;88 0 0 0 -7.93 

196 0 0 0 -7.. 9.4 

204 0 0 ~o -7 •. 94 

212 0 0 0 -7.94 

.220 0 0 0 -7.94 
APPENDIX E: SAMPLE RUN B61 

228 0 0 0 -7.94 

15 0 0 0 37. 

23 0 0 0 25.14 

31 0 0 0 8.99 

39 0 0 0 8.38 

47 0 0 0 6.91 

55 0 0 0 8.63 

63 0 0 0 11.24 

71 0 0 0 7.51 

79 0 0 0 7.36 

87 0 0 0 15.91 

95 0 0 0 21.35 

103 0 (} 0 22.58 

111 0 0 0 18.01 

119 0 0 0 20.84 

127 0 0 0 28. 

135 0 0 0 33. 

143 0 0 0 37.76 

151 0 0 0 40.25 

159 0 0 0 40.25 

167 0 0 0 40.25 

175 0 0 0 40.25 

183 0 0 0 40.25 

191 0 0 0 40.25 

199 0 0 0 40.25 

207 0 0 0 40.25 

215 0 0 0 40.25 

223 0 0 0 40.25 

231 0 0 0 40.25 

TPMP 

6 0 950 

14 0 1030 

22 0 960 

30 0 922 

38 0 935 

46 0 930 

54 0 900 

62 0 839 

70 0 808 

78 0 812 

86 0 866 

94 0 904 

102 0 933 

110 0 1012 

118 0 1102 

126 0 1150 

134 0 1210 

142 0 1320 

150 0 1460 

END 

MACR 

TOL * 1.E-2 

DT * 81.15 

TIME 
. B62 STUDIES OF COLUMBIA GLACIER, ALASKA 

REST * 1. 

LOOP * 10. 

TANF 

SOLV 

CONV 

DT * o. 

TIME 

NEXT 

PLOT 

DISP 

SAVE 

LOOP * 1 

LOOP * 5 

DT * .025 

UPNO 

TIME 

LOOP * 10 

TANF 

SOLV 

CONV 

DT * o. 

TIME 

NEXT 

NEXT 

PLOT 

DISP 

SAVE 

NEXT 

END 

SCAL 

-5oo. 16000. -1000. 1000. 635. 12700. 

END 

STOP 

END OF DATA INPUT 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NUMBER OF NODAL POINTS = 232 

NUMBER OF ELEMENTS = 36 

NUMBER OF MATERIAL SETS = 1 

DIMENSION OF COORDINATE SPACE= 2 

DEGREE OF FREEDOMS/NODE = 3 

NODES PER ELEMENT (MAXIMUM) = 9 

EXTRA D.O.F. TO ELEMENT = 0 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

MATERIAL PROPERTIES 

MATERIAL SET 1 FOR ELEMENT TYPE 4 

'IWO DIMENSIONAL ICE FLOW ELEMENT 

VISCOS. A = 0.14000E+00 

VISCOS. N = 0.30000E+01 

CALVE NODE = 0.19300E+03 

NODES DEEP = 0.80000E+01 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

MACRO INSTRUCTIONS 

• 08 90 • 
APPENDIX E: SAMPLE RUN 

MACRO STATEMENT 

TOL 

VARIABLE 1 VARIABLE 

.10000E-01 .o 

DT 81.150 .o 

TIME .o .o 

REST 1.0000 .o 

LOOP 10.000 .o 

TANF .o .o 

SOLV .o .o 

CONV .o .o 

DT .o .o 

TIME .o .o 

NEXT .o .o 

PLOT .o .o 

DISP .o .o 

SAVE .o .o 

LOOP 1.0000 .o 

LOOP 5. 0000 . .o 

DT .25000E-01 .o 

UPNO .o .o 

TIME .o .o 

LOOP 10.000 .o 

TANF .o .o 

SOLV .o .o 

CONV .o .o 

DT .o .o 

TIME .o .o 

NEXT .o .o 

NEXT .o .o 

PLOT .o .o 

DISP .o .o 

SAVE .o .o 

NEXT .o .o 

END .o .o 

**MACRO INSTRUCTION 1 EXECUTED** TOL 

**MACRO INSTRUCTION 2 EXECUTED** DT 

**MACRO INSTRUCTION 3 EXECUTED** TIME 

**MACRO INSTRUCTION 4 EXECUTED** REST 

**MACRO INSTRUCTION 5 EXECUTED** LOOP 

**MACRO INSTRUCTION 6 EXECUTED** TANF 

FORCE CONVERGENCE TEST 

RNMAX = .39112 RN 

**MACRO INSTRUCTION 7 EXECUTED** SOLV 

**MACRO INSTRUCTION 8 EXECUTED** CONV 

DISPLACEMENT CONVERGENCE TEST 

= 

UNMAX = 13151. UN = 

**MACRO INSTRUCTION 9 EXECUTED** DT 

**MACRO INSTRUCTION 10 EXECUTED** TIME 

**MACRO INSTRUCTION 11 EXECUTED** NEXT 

**MACRO INSTRUCTION 12 EXECUTED** PLOT 

**MACRO INSTRUCTION 13 EXECUTED** DISP 

2 

V1 = 

V1 = 

V1 = 

V1 = 

V1 = 

V1 = 

.39112 

V1 = 

V1 = 

.26771E-01 

V1 = 

V1 = 

V1 = 

V1 = 

V1 = 

B63 

.1000E-01 V2 = .o 

81.15 , V2 = .o 

.o , V2 = .o 

1.000 , V2 = .o 

10.00 , V2 = 12.00 

.o V2 = .o 

TOL = .10000E-01 

.o , V2 = .o 

.o , V2 = .o 

'l'OL = .10000E-01 

.o V2 = .o 

.o V2 = .o 

10.00 V2 = 6.000 

.o V2 = .o 

.o , V2 = .0 
B64 STUDIES OF COLUMBIA GLACIER, ALASKA 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** TIME 81.150 

PROPORTIONAL LOAD 1.0000 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81149994E+02 

NODE !COORD 2 C()()RD 1 DISPL 2 DISPL 3 DISPL 3 Explanation 

1 o.o 465.0000 0.9500E+03 0.2658E+02 0.1454E+Ol pressure 

2 o.o 171.8000 0.9436E+03 0.3344E+02 0.1315E+02 hdot 

3 o.o -121.4000 0.8469E+03 0.4342E+02 0.5157E+02 pressure 

4 o.o -194.7000 0.7848E+03 0.4757E+02 -0.2670E+Ol bdot 

5 o.o -268.0000 0.6982E+03 0.6042E+02 0.6282E+02 pressure 

6 o.o -268.0000 0.6982E+03 0.6042E+02 0.6363E+02 pressure 

7 o.o -273.0000 0.3491E+03 0.3069E+02 0.4000E+02 Abed 

8 o.o -278.0000 o.o 0.0 0.6445E+02 pressure 

9 768.9395 456.6458 0.9785E+03 -0.2287E+02 0.3865E+00 

10 768.9395 179.2141 0.9582E+03 -0.1728E+02 0.8934E+00 

11 768.9395 -98.2178 0.9041E+03 -0.7195E+Ol 0.4951E+02 

12 768.9395 -167.5757 0.8807E+03 -0.2426E+Ol -0.3426E+Ol 

13 768.9395 -236.9341 0.8561E+03 0.3056E+Ol 0.6186E+02 

14 768.9395 -236.9341 0.8561E+03 0.3056E+Ol 0.6235E+02 

15 768.9395 -241.9341 0.4316E+03 0.1645E+Ol 0.3694E+02 

16 768.9395 -246.9341 0.0 o.o 0.6325E+02 

17 1537.8792 429.7080 0.9810E+03 -0.1349E+02 0.6445E+00 

18 1537.8792 157.0576 0.9806E+03 -0.6178E+Ol 0.4564E+00 

19 1537.8792 -115.5928 0.9235E+03 -0.7839E+00 0.4921E+02 

20 1537.8792 -183.7554 0.8882E+03 -0.3930E+Ol -0.3989E+01 

21 1537.8792 -251.9182 0.8384E+03 -0.1038E+02 0.6127E+02 

22 1537.8792 -251.9182 0.8384E+03 -0.1038E+02 0.6072E+02 

23 1537.8792 -256.9182 0.4251E+03 -0.5466E+Ol 0.2490E+02 

24 1537.8792 -261.9182 0.0 o.o 0.6156E+02 

25 2306.8184 422.9641 0.1045E+04 -0.1377E+02 -0.6725E-01 

26 2306.8184 162.8738 0.1024E+04 0.1247E+02 0.2698E+01 

27 2306.8184 -97.2166 0.9461E+03 0.3053E+02 0.4666E+02 

28 2306.8184 -162.2393 0.8948E+03 0.3436E+02 -0.4462E+Ol 

29 2306.8184 -227.2621 0.8223E+03 0.3676E+02 0.5836E+02 

30 2306.8184 -227.2621 0.8223E+03 0.3676E+02 0.5858E+02 

31 2306.8184 -232.2621 0.4134E+03 0.1864E+02 0.8783E+01 

32 2306.8184 -237.2621 o.o 0.0 0.5948E+02 

33 3075.7583 397.1182 0.1138E+04 -0.3268E+02 0.9516E-01 

34 3075.7583 160.9837 0.1114E+04 -0.3918E+Ol 0.2955E+Ol 

35 3075.7583 -75.1506 0.1025E+04 0.2363E+02 0.4208E+02 

36 3075.7583 -134.1843 0.9871E+03 0.2973E+02 -0.4916E+Ol 

37 3075.7583 -193.2181 0.9433E+03 0.3318E+02 0.5281E+02 

38 3075.7583 -193.2181 0.9433E+03 0.3318E+02 0.5299E+02 

39 3075.7583 -198.2181 0.4804E+03 0.1690E+02 0.8344E+01 

40 3075.7583 -203.2181 0.0 0.0 0.5386E+02 

41 3844.6975 369.8826 0.1206E+04 -0.5064E+02 0.2668E+OO 

42 3844.6975 149.3443 0.1182E+04 -0.3620E+02 0.2964E+Ol 

43 3844.6975 -71.1938 0.1100E+04 -0.1966E+02 0.3966E+02 

44 3844.6975 -126.3284 0.1039E+04 -0.1402E+02 -0.5433E+01 

45 3844.6975 -181.4632 0.9477E+03 -0.9526E+01 0.4966E+02 

46 3844.6975 -181.4632 0.9477E+03 -0.9526E+Ol 0.4974E+02 

47 3844.6975 -186.4632 0.4826E+03 -0.4756E+01 0.6918E+01 

48 3844.6975 -191.4632 o.o o.o 0.5063E+02 

49 4613.6367 319.5710 0.1219E+04 -0.8927E+02 0.2130E+01 

so 4613.6367 104.8030 0.1196E+04 -0.1060E+03 -0.3566E+00 
APPENDIX E: SAMPLE RUN 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81149994E+02 

NODE 1 COORD 2 COORD 1 DISPL 

51 4613.6367 -109.9651 0.1124E+04 

52 4613.6367 -163.6570 O.l078E+04 

53 4613.6367 -217.3492 0.1009E+04 

54 4613.6367 -217.3492 0.1009E+04 

55 4613.6367 -222.3492 0.5209E+03 

56 4613.6367 -227.3492 o.o 

57 5382.5742 281.2563 0.1089E+04 

58 5382.5742 42.9097 O.l076E+04 

59 5382.5742 -195.4370 0.9849E+03 

60 5382.5742 -255.0234 0.9346E+03 

61 5382.5742 -314.6104 0.8710E+03 

62 5382.5742 -314.6104 0.8710E+03 

63 5382.5742 -319.6104 0.4392E+03 

64 5382.5742 -324.6104 o.o 

65 6151.5156 270.0581 0.1001E+04 

66 6151.5156 12.2664 0.9778E+03 

67 6151.5156 -245.5256 0.8962E+03 

68 6151.5156 -309.9736 0.8555E+03 

69 6151.5156 . -374.4216 O.B090E+03 

70 6151.5156 -374.4216 0.8090E+03 

71 6151.5156 -379.4216 0.4065E+03 

72 6151.5156 -384.4216 o.o 

73 6920.4570 256.9724 0.9923E+03 

74 6920.4570 0.1812 0.9711E+03 

75 6920.4570 -256.6106 0.8839E+03 

76 6920.4570 -320.8083 0.8274E+03 

77 6920.4570 -385.0063 0.7512E+03 

78 6920.4570 -385.0063 0.7512E+03 

79 6920.4570 -390.0063 0.3774E+03 

80 6920.4570 -395.0063 o.o 

81 7689.3945 239.8603 O.l001E+04 

82 7689.3945 -5.4953 0.9865E+03 

83 7689.3945 -250.8506 0.9322E+03 

84 7689.3945 -312.1895 0.8920E+03 

85 7689.3945 -373.5288 0.8343E+03 

86 7689.3945 -373.5288 0.8343E+03 

87 7689.3945 -378.5288 0.4222E+03 

88 7689.3945 -383.5288 0~0 

89 8458.3359 216.2448 O.l006E+04 

90 8458.3359 -19.8220 0.9992E+03 

91 8458.3359 -255.8888 0.9684E+03 

92 8458.3359 -314.9053 0.9475E+03 

93 8458.3359 -373.9224 0.9187E+03 

94 8458.3359 -373.9224 0.9187E+03 

95 8458.3359 -378.9224 0.4685E+03 

96 8458.3359 -383.9224 o.o 

97 9227.2734 193.8667 0.1005E+04 

98 9227.2734 -40.8036 0.1001E+04 

99 9227.2734 -275.4739 0.9639E+03 

100 9227.2734 -334.1414 0.9376E+03 

2 DISPL 

-0.1178E+03 

-0.1176E+03 

-0.1131E+03 

-0.1131E+03 

-0.5836E+02 

o.o 

-0.1921E+02 

-0.5648E+02 

-0.8647E+02 

-0.9266E+02 

-0.9498E+02 

-0.9498E+02 

-0.4815E+02 

o.o 

0.3958E+01 

-0.8922E+01 

-0.2281E+02 

-0.2435E+02 

-0.2707E+02 

-0.2707E+02 

-0.1364E+02 

o.o 

-0.7751E+Ol 

-0.5736E+Ol 

-0.1563E+OO 

0.8642E+OO 

0.1675E+01 

0.1675E+Ol 

0.8385E+OO 

o.o 

-0.1462E+02 

-0.1085E+02 

0.1088E+Ol 

0.8432E+Ol 

0.1788E+02 

0.1788E+02 

0.9346E+01 

0.0 

-0.2367E+02 

-0.2227E+02 

-0.1869E+02 

-0.1722E+02 

-0.1509E+02 

-0.1509E+02 

-0.7634E+Ol 

o.o 

-0.1446E+02 

-0.1389E+02 

-0.1330E+02 

-0.1393E+02 

3 DISPL 

0.3881E+02 

-0.6002E+01 

0.4791E+02 

0.4783E+02 

0.8741E+Ol 

0.4870E+02 

0.2112E+01 

0.6542E+Ol 

0.4367E+02 

-0.6310E+01 

0.5376E+02 

0.5342E+02 

0.1120E+02 

0.5429E+02 

0.2001E+01 

0.5873E+01 

0.4647E+02 

-0.6320E+01 

0.5796E+02 

0.5769E+02 

0.7387E+01 

0.5858E+02 

0.6951E+OO 

0.6039E+Ol 

0.4632E+02 

-0.6288E+01 

0.5767E+02 

0.5769E+02 

0.7704E+01 

0.5857E+02 

0.7923E+00 

0.5455E+01 

0.4402E+02 

-0.6259E+01 

0.5492E+02 

0.5547E+02 

0.1649E+02 

0.5638E+02 

0.7183E+OO 

0.1431E+01 

0.4247E+02 

-0.6275E+01 

0.5284E+02 

0.5299E+02 

0.2164E+02 

0.5387E+02 

0.8011E+OO 

0.2497E+01 

0.4242E+02 

-0.6374E+01 

B65 
B66 STUDIES OF COLUMBIA GLACIER, ALASKA 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81149994E+02 

NODE 1 COORD 2 COORD 1 DISPL 

101 9227.2734 -392.8093 0.9023E+03 

102 9227.2734 -392.8093 0.9023E+03 

103 9227.2734 -397.8093 0.4529E+03 

104 9227.2734 -402.8093 o.o 

105 9996.2148 177.0835 0.1014E+04 

106 9996.2148 -50.1546 0.1008E+04 

107 9996.2148 -277.3923 0.9791E+03 

108 9996.2148 -334.2019 0.9559E+03 

109 9996.2148 -391.0117 0.9225E+03 

110 9996.2148 -391.0117 0.9225E+03 

111 9996.2148 -396.0117 0.4677E+03 

112 9996.2148 -401.0117 o.o 

113 10765.1523 156.6192 0.1020E+04 

114 10765.1523 -62.8644 0.1015E+04 

115 10765.1523 -282.3479 0.9908E+03 

116 10765.1523 -337.2187 0.9737E+03 

117 10765.1523 -392.0901 0.9501E+03 

118 10765.1523 -392.0901 0.9501E+03 

119 10765.1523 -397.0901 0.4841E+03 

120 10765.1523 -402.0901 o.o 

121 11149.6211 149.5192 0.1024E+04 

122 11149.6211 -65.8687 0.1026E+04 

123 11149.6211 -281.2566 0.9940E+03 

124 11149.6211 -335.1035 0.9772E+03 

125 11149.6211 -388.9509 0.9582E+03 

126 11149.6211 -388.9509 0.9582E+03 

127 11149.6211 -393.9509 0.4736E+03 

128 11149.6211 -398.9509 o.o 

129 11534.0937 148.0525 0.1061E+04 

130 11534.0937 -57.4969 0.1056E+04 

131 11534.0937 -263.0461 0.1029E+04 

132 11534.0937 -314.4336 0.1013E+04 

133 11534.0937 -365.8213 0.9896E+03 

134 11534.0937 -365.8213 0.9896E+03 

135 11534.0937 -370.8213 0.5074E+03 

136 11534.0937 -375.8213 o.o 

137 11918.5625 140.1805 0.1124E+04 

138 11918.5625 -50.5356 0.1108E+04 

139 11918.5625 -241.2516 0.1066E+04 

140 11918.5625 -288.9304 0.1050E+04 

141 11918.5625 -336.6099 0.1030E+04 

142 11918.5625 -336.6099 0.1030E+04 

143 11918.5625 -341.6099 0.5193E+03 

144 11918.5625 -346.6099 o.o 

145 12303.0312 131.3288 0.1210E+04 

146 12303.0312 -41.4878 0.1192E+04 

147 12303.0312 -214.3042 0.1136E+04 

148 12303.0312 -257.5083 0.1114E+04 

149 12303.03~2 -300.7129 0.1088E+04 

150 12303.0312 -300.7129 0.1088E+04 

2 DISPL 

-0.1475E+02 

-0.1475E+02 

-0.7455E+01 

o.o 

-0.6122E+01 

-0.3719E+01 

0.1995E+00 

0.1658E+01 

0.3707E+01 

0.3707E+01 

0.1906E+01 

0.0 

-0.1305E+02 

-0.1277E+02 

-0.1292E+02 

-0.1291E+02 

-0.1263E+02 

-0.1263E+02 

-0.6267E+01 

o.o 

0.6732E+01 

0.1787E+02 

0.2880E+02 

0.3101E+02 

0.3309E+02 

0.3309E+02 

0.1643E+02 

o.o 

0.2465E+02 

0.4870E+02 

0.6691E+02 

0.7144E+02 

0.7523E+02 

0.7523E+02 

0.3856E+02 

0.0 

0.1301E+02 

0.4891E+02 

0.7625E+02 

0.8179E+02 

0.8593E+02 

0.8593E+02 

0.4347E+02 

o.o 

0.6849E+01 

0.5015E+02 

0.9043E+02 

0.9956E+02 

0.1066E+03 

0.1066E+03 

3 DISPL 

0.5275E+02 

0.5262E+02 

0.2238E+02 

0.5350E+02 

0.6530E+OO 

0.1126E+02 

0.4115E+02 

-0.6547E+01 

0.5117E+02 

0.5135E+02 

0.1796E+02 

0.5223E+02 

0.7924E+00 

0.1112E+02 

0.3919E+02 

-0.6705E+01 

0.4880E+02 

0.4882E+02 

0.2214E+02 

0.4971E+02 

0.3306E+OO 

0.9019E+01 

0.3899E+02 

-0.6790E+01 

0.4855E+02 

0.4881E+02 

0.2942E+02 

0.4968E+02 

-0.7999E-01 

0.2809E+02 

0.3699E+02 

-0.6953E+01 

0.4636E+02 

0. 4682E+0'2 

0.3423E+02 

0.4772E+02 

-0.1882E+OO 

0.3065E+02 

0.3374E+02 

-0.7127E+01 

0.4247E+02 

0.4300E+02 

0.3863E+02 

0.4388E+02 

-0.2639E+OO 

0.3638E+02 

0.3090E+02 

-0.7251E+01 

0.3895E+02 

0.3971E+02 
APPENDIX E: SAMPLE RUN B67 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81149994E+02 

NODE 1 COORD 2COORD 1 DISPL 2 DISPL 3 DISPL 

151 12303.0312 -305.7129 0.5562E+03 0.5478E+02 0.4050E+02 

152 12303.0312 -310.7129 o.o o.o 0.4061E+02 

153 12687.5039 116.6096 0.1298E+04 -0.1562E+01 -0.2111E+OO 

154 12687.5039 -37.2140 0.1287E+04 0.3523E+02 0.4183E+02 

155 12687.5039 -191.0374 0.1237E+04 0.7301E+02 0.2695E+02 

156 12687.5039 -229.4932 0.1215E+04 0.8263E+02 -0.7442E+01 

157 12687.5039 -267.9495 0.1189E+04 0.9161E+02 0.3394E+02 

158 12687.5039 -267.9495 0.1189E+04 0.9161E+02 0.3488E+02 

159 12687.5039 -272.9495 0.6078E+03 0.4715E+02 0.4025E+02 

160 12687.5039 -277.9495 o.o 0.0 0.3576E+02 

161 13071.9727 101.4698 0.1391E+04 0.1581E+02 -0.3078E+OO 

162 13071.9727 -34.5920 0.1386E+04 0.5316E+02 0.6919E+02 

163 13071.9727 -170.6537 0.1342E+04 0.8741E+02 0.2428E+02 

164 13071.9727 -204.6691 0.1325E+04 0.9566E+02 -0.7724E+01 

165 13071.9727 -238.6848 0.1304E+04 0.1033E+03 0.3045E+02 

166 13071.9727 -238.6848 0.1304E+04 0.1033E+03 0.3135E+02 

167 13071.9727 -243.6848 0.6799E+03 0.5462E+02 0.4025E+02 

168 13071.9727 -248.6848 o.o 0.0 0.3225E+02 

169 13264.2031 94.3128 0.1473E+04 0.9351E+01 -0.5730E+OO 

170 13264.2031 -31.4634 0.1454E+04 0.7239E+02 0.3975E+02 

171 13264.2031 -157.2396 0.1403E+04 0.1234E+03 0.2189E+02 

172 13264.2031 -188.6836 0.1385E+04 0.1339E+03 -0.7833E+01 

173 13264.2031 -220.1277 0.1363E+04 0.1437E+03 0.2756E+02 

174 13264.2031 -220.1277 0.1363E+04 0.1437E+03 0.2865E+02 

175 13264.2031 -225.1277 0.6702E+03 0.7095E+02 0.4025E+02 

176 13264 •, 2031 -230.1277 o.o 0.0 0.2952E+02 

177 13456.4414 89.8647 0.1602E+04 -0.7009E+01 -0.1291E+01 

178 13456.4414 -24.7626 0.1576E+04 0.1070E+03 0.6118E+02 

179 13456.4414 -139.3900 0.1492E+04 0.1856E+03 0.2025E+02 

180 13456.4414 -168.0467 0.1469E+04 0.2002E+03 -0.7910E+01 

181 13456.4414 -196.7039 0.1450E+04 0.2093E+03 0.2610E+02 

182 13456.4414 -196.7039 0.1450E+04 0.2093E+03 0.2713E+02 

183 13456.4414 -201.7039 0.8329E+03 0.1207E+03 0.4025E+02 

184 13456.4414 -206.7039 o.o 0.0 0.2811E+02 

185 13648.6719 76.2207 0.1824E+04 -0.4276E+02 -0.8349E+OO 

186 13648.6719 -20.9452 0.1795E+04 0.6911E+02 0.6717E+02 

187 13648.6719 -118.1111 0.1689E+04 0.1934E+03 0.1660E+02 

188 13648.6719 -142.4025 0.1643E+04 0.2221E+03 -0.7942E+01 

189 13648.6719 -166.6940 0.1591E+04 0.2464E+03 0.2117E+02 

190 13648.6719 -166.6940 0.1591E+04 0.2464E+03 0.2295E+02 

191 13648.6719 -171.6940 0.7263E+03 0.1141E+03 0.4025E+02 

192 13648.6719 -176.6940 o.o o.o 0.2379E+02 

193 13840.9141 66.9960 0.1955E+04 -0.2757E+02 -0.1211E+01 

194 13840.9141 -13.7014 0.2063E+04 0.4808E+02 0.1659E+02 

195 13840.9141 -94.3988 0.2034E+04 0.1827E+03 0.1231E+02 

196 13840.9141 -114.5732 0.2000E+04 0.2339E+03 -0.7941E+01 

197 13840.9141 -134.7476 0.1959E+04 0.2749E+03 0.1806E+02 

198 13840.9141 -134.7476 0.1959E+04 0.2749E+03 0.2018E+02 

199 13840.9141 -139.7476 0.1570E+04 0.2394E+03 0.4025E+02 

200 13840.9141 -144.7476 o.o o.o 0.2173E+02 
B68 STUDIES OF COLUMBIA GLACIER, ALASKA 

**MACRO INSTRUCTION 14 EXECUTED** SAVE V1 = • 0' I V2 = .o 

**MACRO INSTRUCTION 15 EXECUTED** LOOP V1 = 1.000 I V2 = 32.00 

**MACRO INSTRUCTION 16 EXECUTED** LOOP V1 = 5.000 I V2 = 28.00 

**MACRO INSTRUCTION 17 EXECUTED** DT V1 = .2500E-01 I V2 = .o 

**MACRO INSTRUCTION 18 EXECUTED** UPNO V1 = .o I V2 = .o 

**MACRO INSTRUCTION 19 EXECUTED** TIME V1 = .o V2 = .o 

**MACRO INSTRUCTION 20 EXECUTED** LOOP V1 = 10.00 I V2 = 27.00 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 1985.6 RN = 1985.6 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o I V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o I V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 13008. UN = 166.76 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o I V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o I V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 1.000 I V2 = 21.00 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 15.457 RN = 15.457 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o I V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o I V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 13010. UN = 2.4093 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o I V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o I V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 10.00 V2 = 21.00 

**MACRO INSTRUCTION 27 EXECUTED** NEXT V1 = 1.000 V2 = 17.00 

**MACRO INSTRUCTION 17 EXECUTED** DT V1 = .2500E-01 I V2 = .o 

**MACRO INSTRUCTION 18 EXECUTED** UPNO V1 = .o V2 = .o 

**MACRO INSTRUCTION 19 EXECUTED** TIME V1 = .o I V2 = .o 

**MACRO INSTRUCTION 20 EXECUTED** LOOP V1 = 10.00 I V2 = 27.00 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o I V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 2440.2 RN = 2440.2 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o I V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o I V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 12846. UN = 195.34 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o I V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 1.000 I V2 = 21.00 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 19.483 RN = 19.483 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o I V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o I V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 12848. UN = 3.0509 TOL = .10000E-01 

. **MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o I V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o I V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 10.00 V2 = 21.00 
APPENDIX E: SAMPLE RUN B69 

**MACRO INSTRUCTION 27 EXECUTED** NEXT V1 = 2.000 V2 = 17.00 

**MACRO INSTRUCTION 17 EXECUTED** DT V1 = .2500E-01 ' V2 = .o 

**MACRO INSTRUCTION 18 EXECUTED** UPNO V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 19 EXECUTED** TIME V1 = .o V2 = .o 

**MACRO INSTRUCTION 20 EXECUTED** LOOP V1 = 10.00 ' V2 = 27.00 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o ' V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 2394.2 RN = 2394.2 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o , V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 12704. UN = 175.48 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 1.000 V2 = 21.00· 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 27.012 RN = 27.012 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o , V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 12706. UN = 2.7596 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 10.00 ' V2 = 21.00 

**MACRO INSTRUCTION 27 EXECUTED** NEXT V1 = 3.000 ' V2 = 17.00 

**MACRO INSTRUCTION 17 EXECUTED** DT V1 = .2500E-01 V2 = .o 

**MACRO INSTRUCTION 18 EXECUTED** UPNO V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 19 EXECUTED** TIME V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 20 EXECUTED** LOOP V1 = 10.00 ' V2 = 27.00 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 1750.1 RN = 1750.1 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o , V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 12572. UN = 152.09 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 1.000 ' V2 = 21.0Q 

**MACRO INSTRUCTION 21 EXECUTED** TANF V1 = .o ' V2 = .o 

FORCE CONVERGENCE TEST 

RNMAX = 127.53 RN = 127.53 TOL = .10000E-01 

**MACRO INSTRUCTION 22 EXECUTED** SOLV V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 23 EXECUTED** CONV V1 = .o , V2 = .o 

DISPLACEMENT CONVERGENCE TEST 

UNMAX = 12574. UN = 2.6256 TOL = .10000E-01 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .o ' V2 = .o 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 10.00 ' V2 = 21.00 

**MACRO INSTRUCTION 27 EXECUTED** NEXT V1 = 4.000 ' V2 = 17.00 

**MACRO INSTRUCTION 17 EXECUTED** DT V1 = .2500E-01 ' V2 = .o 
B70 STUDIES OF COLUMBIA GLACIER, ALASKA 

**MACRO INSTRUCTION 18 EXECUTED** UPNO .o 

**MACRO INSTRUCTION 19 EXECUTED** TIME 

**MACRO INSTRUCTION 20 EXECUTED** LOOP 

**MACRO INSTRUCTION 21 EXECUTED** TANF 

V1 = 

V1 ·= 

V1 = 

.o 

10.00 

FORCE CONVERGENCE TEST 

RNMAX = 1291.0 RN 

**MACRO INSTRUCTION 22 EXECUTED** SOLV 

**MACRO INSTRUCTION 23 EXECUTED** CONV 

DISPLACEMENT CONVERGENCE TEST 

= 

V1 = 

1291.0 

V1 = 

V1 = 

UNMAX = 1246~. UN = 123.86 

.o 

.o 

.o 

**MACRO INSTRUCTION 24 EXECUTED** DT V1 = .0 

**MACRO INSTRUCTION 25 EXECUTED** TIME V1 = .0 

TOL 

TOL 

**MACRO INSTRUCTION 26 EXECUTED** NEXT V1 = 10. 00 

**MACRO INSTRUCTION 27 EXECUTED** NEXT V1 = 5.000 

**MACRO INSTRUCTION 28 EXECUTED** PLOT V1 = • 0 

**MACRO INSTRUCTION 29 EXECUTED** DISP V1 = • 0 

V2 = 

V2 = 

V2 = , V2 = 

= 

I V2 = 

, V2 = 

= 

, V2 = , V2 = , V2 = 

V2 = , V2 = 

V2 = 

.o 

.o 

27.00 

.o 

.10000E-01 

.o 

.o 

.10000E-01 

.o 

.o 

21.00 

17.00 

.o 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

PROPORTIONAL LOAD 1.0000 

TIME 

.o 

81.275 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81274963E+02 

NODE 1 COORD 2 COORD 1 DISPL 

1 0.0 465.0000 0.9500E+03 

2 0.0 171.8000 0.9436E+03 

3 0.0 -121.4000 0.8469E+03 

4 0.0 -194.7000 0.7848E+03 

5 0.0 -268.0000 0.6982E+03 

6 0.0 -268.0000 0.6982E+03 

7 0.0 -273.0000 0.3491E+03 

8 o.o -278.0000 o.o 

9 774.9919 456.7146 0.9753E+03 

10 774.9919 179.2788 0.9557E+03 

11 774.9919 -98.1570 0.9024E+03 

12 774.9919 -167.5159 0.8791E+03 

13 774.9919 -236.8751 0.8548E+03 

14 774.9919 -236.8751 0.8548E+03 

15 774.9919 -241.8751 0.4309E+03 

16 774.9919 -246.8751 o.o 

17 1549.9841 429.6428 0.9762E+03 

18 1549.9841 157.0532 0.9760E+03 

19 1549.9841 -115.5364 0.9192E+03 

20 1549.9841 -183.6838 0.8838E+03 

21 1549.9841 -251.8313 0.8338E+03 

22 1549.9841 -251.8313 0.8338E+03 

23 1549.9841 -256.8313 0.4227E+03 

24 1549.9841 -261.8313 o.o 

25 2324.9758 423.0266 0.1037E+04 

26 2324.9758 163.1912 0.1017E+04 

27 2324.9758 -96.6445 0.9395E+03 

28 2324.9758 -161.6035 0.8878E+03 

29 2324.9758 -226.5626 0.8145E+03 

30 2324.9758 -226.5626 0.8145E+03 

31 2324.9758 -231.5626 0.4095E+03 

32 2324.9758 -236.5626 o.o 

33 3099.9685 396.8281 0.1126E+04 

2 DISPL 

0.2882E+02 

0.3435E+02 

0.4353E+02 

0.4749E+02 

0.6017E+02 

0.6017E+02 

0.3056E+02 

o.o 

-0.2072E+02 

-0.1599E+02 

-0.6765E+01 

-0.2233E+01 

0.3018E+01 

0.3018E+01 

0.1622E+01 

o.o 

-0.1030E+02 

-0.3824E+01 

0.5947E+OO 

-0.2891E+01 

-0.9678E+01 

-0.9678E+01 

-0.5113E+01 

0.0 

-0.1124E+02 

0.1365E+02 

0.3074E+02 

0.3431E+02 

0.3651E+02 

0.3651E+02 

0.1851E+02 

o.o 

-0.2900E+02 

3 DISPL 

0.1535E+01 

0.1517E+02 

0.5159E+02 

-0.2670E+01 

0.6284E+02 

0.6365E+02 

0.4000E+02 

0.6447E+02 

0.4242E+OO 

0.2853E+01 

0.4952E+02 

-0.3431E+01 

0.6187E+02 

0.6234E+02 

0.3688E+02 

0.6324E+02 

0.6830E+00 

0.3248E+Ol 

0.4924E+02 

-0.3997E+01 

0.6130E+02 

0.6073E+02 

0.2468E+02 

0.6157E+02 

-0.5320E-01 

0.5053E+01 

0.4662E+02 

-0.4473E+01 

0.5830E+02 

0.5852E+02 

0.8597E+01 

0.5942E+02 

0.1372E+00 
APPENDIX E: SAMPLE RUN B71 

34 3099.9685 161.0949 0.1103E+04 -0.2955E+Ol 0.6358E+Ol 

35 3099.9685 -74.6382 0.1017E+04 0.2284E+02 0.4203E+02 

36 3099.9685 -133.5715 0.9795E+03 0.2885E+02 -0.4932E+Ol 

37 3099.9685 -192.5051 0.9362E+03 0.3245E+02 0.5274E+02 

38 3099.9685 -192.5051 0.9362E+03 0.3245E+02 0.5293E+02 

39 3099.9685 -197.5051 0.4769E+03 0.1656E+02 0.8310E+Ol 

40 3099.9685 -202.5051 o.o o.o 0.5380E+02 

41 3874.9602 368.8999 0.1185E+04 -0.4704E+02 0.3037E+OO 

42 3874.9602 148.5506 0.1162E+04 -0.3464E+02 0.6679E+Ol 

43 3874.9602 -71.7986 0.1084E+04 -0.2039E+02 0.3967E+02 

44 3874.9602 -126.8860 0.1026E+04 -0.1568E+02 -0.5454E+Ol 

45 3874.9602 -181.9735 0.9358E+03 -0.1186E+02 0.4961E+02 

46 3874.9602 -181.9735 0.9358E+03 -0.1186E+02 0.4968E+02 

47 3874.9602 -186.9735 0.4768E+03 -0.5977E+Ol 0.6931E+Ol 

48 3874.9602 -191.9735 o.o o.o 0.5056E+02 

49 4649.9492 317.6731 0.1192E+04 -0.8114E+02 0.2195E+Ol 

50 4649.9492 102.3462 0.1170E+04 -0.1003E+03 0.3837E+Ol 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81274963E+02 

NODE !COORD 2COORD 1 DISPL 2 DISPL 3 DISPL 

51 4649.9492 -112.9807 0.1094E+04 -0.1151E+03 0.3902E+02 

52 4649.9492 -166.8125 0.1048E+04 -0.1162E+03 -0.6023E+Ol 

53 4649.9492 -220.6443 0.9791E+03 -0.1126E+03 0.4809E+02 

54 4649.9492 -220.6443 0.9791E+03 -0.1126E+03 0.4797E+02 

55 4649.9492 -225.6443 0.5050E+03 -0.5810E+02 0.8847E+Ol 

56 4649.9492 -230.6443 o.o o.o 0.4884E+02 

57 5424.9414 281.0791 0.1054E+04 -0.9754E+Ol 0.2136E+Ol 

58 5424.9414 41.1132 0.1039E+04 -0.4953E+02 0.1298E+02 

59 5424.9414 -198.8525 0.9476E+03 -0.8097E+02 0.4398E+02 

60 5424.9414 -258.8440 0.8986E+03 -0.8734E+02 -0.6318E+Ol 

61 5424.9414 -318.8357 0.8368E+03 -0.8974E+02 0.5414E+02 

62 5424.9414 -318.8357 0.8368E+03 -0.8974E+02 0.5380E+02 

63 5424.9414 -323.8357 0.4217E+03 -0.4547E+02 0.1113E+02 

64 5424.9414 -328.8357 0.0 o.o 0.5467E+02 

65 6199.9336 270.1453 0.9592E+03 0.5920E+Ol 0.2010E+Ol 

66 6199.9336 11.4897 0.9382E+03 -0.6883E+Ol 0.8312E+Ol 

67 6199.9336 -247.1658 0.8640E+03 -0.1952E+02 0.4661E+02 

68 6199.9336 -311.8296 0.8260E+03 -0.2080E+02 -0.6319E+Ol 

69 6199.9336 -376.4937 0.7824E+03 -0.2302E+02 0.5812E+02 

70 6199.9336 -376.4937 0.7824E+03 -0.2302E+02 0.5789E+02 

71 6199.9336 -381.4937 0.3934E+03 -o·.ll61E+02 0.7286E+Ol 

72 6199.9336 -386.4937 o.o o.o 0.5878E+02 

73 6974.9297 256.8132 0.9498E+03 -0.3999E+Ol 0.8032E+00 

74 6974.9297 0.1335 0.9307E+03 -0.2694E+Ol 0.9165E+Ol 

75 6974.9297 -256.5461 0.8516E+03 0.1483E+Ol 0.4629E+02 

76 6974.9297 -320.7161 0.7993E+03 0.2469E+Ol -0.6286E+Ol 

77 6974.9297 -384.8862 0.7284E+03 0.2891E+Ol 0.5766E+02 

78 6974.9297 -384.8862 0.7284E+03 0.2891E+Ol 0.5767E+02 

79 6974.9297 -389.8862 0.3661E+03 0.1450E+Ol 0.8074E+Ol 

80 6974.9297 -394.8862 o.o 0.0 0.5855E+02 

81 7749.9180 239.0933 0.9555E+03 -0.1150E+02 0.1017E+Ol 

82 7749.9180 -5.8043 0.9420E+03 -0.9521E+Ol 0.8511E+Ol 

83 7749.9180 -250.7019 0.8951E+03 0.1285E+00 0.4395E+02 

84 7749.9180 -311.9260 0.8576E+03 0.7086E+Ol -0.6259E+Ol 
B72 STUDIES OF COLUMBIA GLACIER, ALASKA 

85 7749.9180 -373.1509 0_. 8030E+03 0.1561E+02 0.5488E+02 

86 7749.9180 -373.1509 0.8030E+03 0.1561E+02 0.5538E+02 

87 7749.9180 -378.1509 0.4069E+03 0.8168E+01 0.1701E+02 

88 7749.9180 -383.1509 0.0 o.o 0.5629E+02 

89 8524.9102 214.6174 0.9552E+03 -0.1915E+02 0.9557E+OO 

90 8524.9102 -21.2081 0.9497E+03 -0.1881E+02 0.4678E+01 

91 8524.9102 -257.0334 0.9224E+03 -0.1738E+02 0.4242E+02 

92 8524.9102 -315.9897 0.9030E+03 -0.1636E+02 -0.6280E+01 

93 8524.9102 -374.9465 0.8762E+03 -0.1495E+02 0.5280E+02 

94 8524.9102 -374.9465 0.8762E+03 -0.1495E+02 0.5290E+02 

95 8524.9102 -379.9465 0.4469E+03 -0.7574E+01 0.2188E+02 

96 8524.9102 -384.9465 0.0 o.o 0.5378E+02 

97 9299.9023 192.5049 0.9525E+03 -0.8651E+01 0.1247E+01 

98 9299.9023 -41.8501 0.9479E+03 -0.9233E+O~ 0.4868E+01 

99 9299.9023 -276.2048 0.9116E+03 -0.1030E+02 0.4240E+02 

100 9299.9023 -334.7937 0.8873E+03 -0.1063E+02 -0.6387E+01 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81274963E+02 

NODE 1 COORD 2 COORD 1 DISPL 

101 9299.9023 -393.3828 0.8549E+03 

102 9299.9023 -393.3828 0.8549E+03 

103 9299.9023 -398.3828 0.4285E+03 

104 9299.9023 -403.3828 o.o 

105 10074.8945 176.7471 0.9544E+03 

106 10074.8945 -50.3697 0.9492E+03 

107 10074.8945 -277.4863 0.9231E+03 

108 10074.8945 -334.2656 0.9017E+03 

109 10074.8945 -391.0452 0.8706E+03 

110 10074.8945 -391.0452 0.8706E+03 

111 10074.8945 -396.0452 0.4419E+03 

112 10074.8945 -401.0452 0.0 

113 10849.8867 155.4968 0.9552E+03 

114 10849.8867 -63.5514 0.9507E+03 

115 10849.8867 -282.5994 0.9269E+03 

116 10849.8867 -337.3613 0.9104E+03 

117 10849.8867 -392.1238 0.8885E+03 

118 10849.8867 -392.1238 0.8885E+03 

119 10849.8867 -397.1238 0.4510E+03 

120 10849.8867 -402.1238 o.o 

121 11237.3828 150.2279 0.9566E+03 

122 11237.3828 -63.8597 0.9571E+03 

123 11237.3828 -277.9473 0.9271E+03 

124 11237.3828 -331.4692 0.9118E+03 

125 11237.3828 -384.9915 0.8948E+03 

126 11237.3828 -384.9915 0.8948E+03 

127 11237.3828 -389.9915 0.4419E+03 

128 11237.3828 -394.9915 o.o 

129 11624.8789 150.7312 0.9885E+03 

130 11624.8789 -53.4030 0.9822E+03 

131 11624.8789 -257.5371 0.9575E+03 

132 11624.8789 -308.5706 0.9423E+03 

133 11624.8789 -359.6042 0.9215E+03 

134 11624.8789 -359.6042 0.9215E+03 

135 11624.8789 -364.6042 0.4721E+03 

2 DISPL 

-0.1176E+02 

-0.1176E+02 

-0.5936E+01 

o.o 

-0.2155E+01 

-0.1550E+01 

0.6115E+00 

0.1667E+01 

0.3158E+01 

0.3158E+01 

0.1616E+01 

o.o 

-0.3793E+01 

-0.4707E+01 

-0.6922E+01 

-0.7465E+01 

-0.7910E+01 

-0.7910E+01 

-0.3895E+01 

o.o 

0.1928E+02 

0.2765E+02 

0.3572E+02 

0.3718E+02 

0.3876E+02 

0.3876E+02 

0.1922E+02 

o.o 

0.2897E+02 

0.5029E+02 

0.6617E+02 

0.7013E+02 

0.7316E+02 

0.7316E+02 

0.3746E+02 

3 DISPL 

0.5282E+02 

0.5264E+02 

0.2214E+02 

0.5353E+02 

0.8901E+OO 

0.1362E+02 

0.4112E+02 

-0.6564E+01 

0.5113E+02 

0.5127E+02 

0.1799E+02 

0.5215E+02 

0.1113E+01 

0.1871E+02 

0.3919E+02 

-0.6725E+01 

0.4878E+02 

0.4873E+02 

0.2359E+02 

0.4963E+02 

0.4064E+OO 

0.1498E+02 

0.3892E+02 

-0.6823E+01 

0.4837E+02 

0.4863E+02 

0.3064E+02 

0.4951E+02 

-0.4567E-01 

0.3111E+02 

0.3661E+02 

-0.6995E+01 

0.4593E+02 

0.4634E+02 

0.3537E+02 
APPENDIX E: SAMPLE RUN B73 

136 11624.8789 -369.6042 o.o o.o 0.4724E+02 

137 12012.3750 142.2204 0.1045E+04 0.1843E+02 -0.1379E+00 

138 12012.3750 -45.9574 0.1031E+04 0.4992E+02 0.3542E+02 

139 12012.3750 -234.1351 0.9915E+03 0.7445E+02 0.3338E+02 

140 12012.3750 -281.1792 0.9765E+03 0.7946E+02 -0.7158E+01 

141 12012.3750 -328.2241 0.9587E+03 0.8321E+02 0.4200E+02 

142 12012.3750 -328.2241 0.9587E+03 0.8321E+02 0.4252E+02 

143 12012.3750 -333.2241 0.4834E+03 0.4212E+02 0.3931E+02 

144 12012.3750 -338.2241 o.o o.o 0.4341E+02 

145 12399.8711 133.1550 0.1119E+04 0.1276E+02 -0.1695E+OO 

146 12399.8711 -36.9271 0.1104E+04 0.4695E+02 0.4065E+02 

147 12399.8711 -207.0090 0.1056E+04 0.7955E+02 0.3036E+02 

148 12399.8711 -249.5296 0.1036E+04 0.8724E+02 -0.7293E+01 

149 12399.8711 -292.0505 0.1012E+04 0.9337E+02 0.3825E+02 

lSD 12199.8711 -292.0505 0.1012E+04 0.9337E+02 0.3897E+02 

FEAP COLUMBIA GLACIER VERTICAL CENTERLINE PROFILE ******************** 

NODAL DISPLACEMENTS TIME 0.81274963E+02 

NODE !COORD 2 COORD 1 DISPL 2 DISPL 3 DISPL 

151 12399.8711 -297.0505 0.5201E+03 0.4817E+02 0.4056E+02 

152 12399.8711 -302.0505 o.o o.o 0.3987E+02 

153 12787.3672 118.2029 0.1189E+04 0.6555E+Ol -0.1346E+OO 

154 12787.3672 -33.4582 0.1181E+04 0.3665E+02 0.4326E+02 

155 12787.3672 -185.1193 0.1139E+04 0.6797E+02 0.2667E+02 

156 12787.3672 -223.0344 0.1121E+04 0.7596E+02 -0.7505E+Ol 

157 12787.3672 -260.9500 0.1099E+04 0.8338E+02 0.3354E+02 

158 12787.3672 -260.9500 0.1099E+04 0.8338E+02 0.3440E+02 

159 12787.3672 -265.9500 0.5602E+03 0.4281E+02 0.4025E+02 

160 12787.3672 -270.9500 o.o o.o 0.3529E+02 

161 13174.8633 104.6301 0.1271E+04 0.2728E+02 -0.3111E+OO 

162 13174.8633 -29.1245 0.1270E+04 0.6326E+02 0.6827E+02 

163 13174.8633 -162.8789 0.1233E+04 0.9604E+02 0.2384E+02 

164 13174.8633 -196.3176 0.1218E+04 0.1039E+03 -0.7783E+01 

165 13174.8633 -229.7564 0.1199E+04 0.1111E+03 0.2995E+02 

166 13174.8633 -229.7564 0.1199E+04 0.1111E+03 0.3087E+02 

167 13174.8633 -234.7564 0.6223E+03 0.5801E+02 0.4025E+02 

168 13174.8633 -239.7564 o.o o.o 0.3177E+02 

169 13368.6094 97.0733 0.1357E+04 0.2381E+02 -0.5945E+00 

170 13368.6094 -24.9426 0.1339E+04 0.8648E+02 0.5714E+02 

171 13368.6094 -146.9585 0.1291E+04 0.1334E+03 0.2130E+02 

172 13368.6094 -177.4623 0.1274E+04 0.1423E+03 -0.7880E+01 

173 13368.6094 -207.9665 0.1257E+04 0.1504E+03 0.2677E+02 

174 13368.6094 -207.9665 0.1257E+04 0.1504E+03 0.2781E+02 

175 13368.6094 -212.9665 0.6225E+03 0.7506E+02 0.4025E+02 

176 13368.6094 -217.9665 o.o o.o 0.2869E+02 

177 13562.3594 92.0420 0.1496E+04 -0.5806E-01 -0.1345E+01 

178 13562.3594 -17.4035 0.1464E+04 0.1197E+03 0.5321E+02 

179 13562.3594 -126.8490 0.1375E+04 0.2031E+03 0~1920E+02 

180 13562.3594 -154.2104 0.1351E+04 0.2184E+03 -0.7935E+01 

181 13562.3594 -181.5718 0.1331E+04 0.2276E+03 0.2491E+02 

182 13562.3594 -181.5718 0.1331E+04 0.2276E+03 0.2601E+02 

183 13562.3594 -186.5718 0.7400E+03 0.1243E+03 0.4025E+02 

184 13562.3594 -191.5718 o.o o.o 0.2697E+02 

185 13756.1055 78.2083 O.l743E+04 -0.4757E+02 -0.9306E+00 

186 13756.1055 -11.9215 0.1707E+04 0.6358E+02 0.8268E+02 
B74 STUDIES OF COLUMBIA GLACIER, ALASKA 

187 13756.1055 -102.0513 0.1595E+04 0.1812E+03 0.1544E+02 

188 13756.1055 -124.5837 0.1548E+04 0.2081E+03 -0.7941E+01 

189 13756.1055 -147.1162 0.1495E+04 0.2310E+03 0.1974E+02 

190 13756.1055 -147.1162 0.1495E+04 0.2310E+03 0.2155E+02 

191 13756.1055 -152.1162 0.7126E+03 0.1115E+03 0.4025E+02 

192 13756.1055 -157.1162 o.o o.o 0.2240E+02 

193 13949.8594 62.7661 0.1895E+04 -0.7402E+02 -0.1078E+01 

194 13949.8594 -10.3158 0.1965E+04 -0.1467E+02 0.6399E+02 

195 13949.8594 -83.3978 0.1922E+04 0.9215E+02 O.l077E+02 

196 13949.8594 -101.6682 0.1889E+04 0.1333E+03 -0.7940E+01 

197 13949.8594 -119.9387 0.1850E+04 0.1682E+03 0.1576E+02 

198 13949.8594 -119.9387 0.1850E+04 0.1682E+03 0.1780E+02 

199 13949.8594 -124.9387 0.1454E+04 0.1565E+03 0.4025E+02 

200 13949.8594 -129.9387 o.o 0.0 o.1928E+02 

**MACRO INSTRUCTION 30 EXECUTED** SAVE V1 = .o I V2 = .o 

**MACRO INSTRUCTION 31 EXECUTED** NEXT V1 = 1.000 I V2 = 16.00 

DATA SET UTILITY - GENERATE PAGE 0001 

PROCESSING ENDED AT EOD 

'k U.S. GOVERNMENT PRINTING OFFICE: 1982-361-614/11