program hicon
C PROGRAM TO MODEL DYNAMICS OF FLOW OF HAWAIIAN MAGMA/GAS
C MIXTURES IN VERTICAL ERUPTIVE CONDUITS. THE MODEL ASSUMES
C STEADY STATE, HOMOGENEOUS FLOW OF LIQUID MAGMA AND IDEAL
C GAS. THE COMPLETE DOCUMENTATION TO THIS CODE IS EXPLAINED
C IN THE U.S.G.S. OPEN-FILE REPORT "A NUMERICAL PROGRAM FOR
C STEADY-STATE FLOW OF HAWAIIAN MAGMA-GAS MIXTURES THROUGH
C VERTICAL ERUPTIVE CONDUITS", OPEN FILE REPORT NUMBER 95-756
C
C WRITTEN BY LARRY G. MASTIN
C U.S. GEOLOGICAL SURVEY
C CASCADES VOLCANO OBSERVATORY
C 1300 SE Cardinal Court,
C Bldg.10, Suite 100
C VANCOUVER, WA 98683
C TEL. (360) 696-7518
C E-MAIL lgmastin@usgs.gov
C
C DECLARE VARIABLES
implicit double precision (a-h,o-z)
real*8 area(2000),output(37,2000),pres(2000),re(2000),rho(2000),
* time(2000),vel(2000),vesic(2000),visc(2000),z(2000)
real*8 mach(2000)
real*8 blkmag,cm,co2,cp,cph2o,cpco2,cps,cv,cvh2o,
* cvco2,cvs,dadz,diam,dpdz,dz,dznext,
* dzout,eps,eta,exco2,exh2o,exmol,exsulfur,f,fo,g,gamma,
* h2o,m,mco2,mdot,mf,mh2o,mm,ms,p,pi,pout,pgrad,r,rey,rhof,rhom,
* rhomix,sulfur,sv,t1,temp,v,vfgas,wtdepth,wtmin,
* xco2,xh2o,xsarea,xsulfur,zin
integer i,icalc,iend,ifinal,im,in,ip,io,ivar(12),ivars,ives,
* ivt,ix
character title(37)*10, outfile*40
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
common/ints/i,icalc,iend,im,ip,io,ives
data title/' area',' mach #',' p (MPa)',
* ' log Re',' density',' time (s)',' vel (m/s)',
* ' vfgas',' log visc',' z (m)',' da/dz',
* ' log(p)',' dpdz',' dz (m)',' wt% exco2',
* ' wt% exh2o',' wt% exs',' f',' gamma',
* ' mgas/mtot',' mmag/mtot',' R (J/m K)',' gas dens',
* ' sonic vel',' xco2',' xh2o',' xsulfur',
* ' cph2o',' cpco2',' cpsulfur',' cvh2o',
* ' cvco2',' cvsulfur',' temp(C)',' h (kJ/kg)',
* ' cp',' radius'/
C OPEN INPUT AND OUTPUT FILES
open (unit=8,file='hicin')
C QUERY USER FOR NAME OF OUTPUT FILE
write (6,1342)
1342 format (/,5x,'enter name of output file:'$)
read (5,1343) outfile
1343 format (a40)
open (unit=9,file=outfile)
C READ INPUT DATA
read (8,*)
read (8,*) icalc
read (8,*) pres(1)
read (8,*) pgrad
read (8,*) ivt
read (8,*) vel(1)
read (8,*) rhom
read (8,*) t1
30 read (8,*) h2o, co2, sulfur
if ((h2o.eq.0.).and.(sulfur.ne.0.).and.(co2.ne.0.)) then
write (6,265)
265 format (/,' Sorry, I can"t handle zero-values of H2O and',/,
* ' non-zero values of CO2 and sulfur. Please try',/,
* ' another combination.')
stop
end if
read (8,*) ives
read (8,*) z(1)
if (z(1).gt.0.) z(1)=-z(1)
read (8,*) diam
read (8,*) fo
C READ PARAMETERS TO WRITE TO OUTPUT FILE
do 1039 iparam=1,18
1039 read (8,*)
ivars=0
do 1041 iin=1,37
read (8,31,err=1041) ix
31 format (i1)
if (ix.le.0) go to 1041
ivar(ix)=iin
ivars = ivars+1
1041 continue
do 1044 ix=1,ivars
if (ivar(ix).eq.0) then
write (6,1043)
1043 format (/, ' output parameters incorrectly listed.',
* /,' program stopped.')
stop
end if
1044 continue
C SET INITIAL VALUES
g = 9.81d+00
if (icalc.eq.2) then
C CONVERT PRESSURE GRADIENT FROM Pa/m TO MPa/km
dpdz = -pgrad*1000.
pres(1) = 1.013d+05 + dpdz*z(1)
ivt = 1
eps = 1.0d-03
imax = 2000
else
dadz = 0.
pres(1) = pres(1)*1.0d+06
eps = 3.0d-07
imax = 1000
end if
blkmag = 1.0d+11
cm = 1.0d+03
i = 1
iend = 0
im = 1
ip = 1
iruns = 0
mco2 = 0.0440
mh2o = 0.01801
ms = 0.032064
pi = 3.14159d+00
r = 8.314
area(1) = pi * (diam/2.)**2
xsarea = area(1)
C WRITE OUT INITIAL CONDITIONS
write (6,134) vel(1), rhom, t1, fo, h2o, co2, sulfur
134 format (/,10x,'INPUT VALUES:',//,
* 10x,' input velocity = ',f8.4,' m/s',/,
* 10x,' magma density=',f7.0,' kg/m3',/,
* 10x,' input temperature =',f6.0,' degrees Celsius',/,
* 10x,' fo (wall roughness) =',f7.4,/,
* 10x,' initial dissolved h2o=',f6.3,' wt %',/,
* 10x,' initial dissolved co2=',f6.3,' wt %',/,
* 10x,' initial dissolved S=',f6.3,' wt %')
if (icalc.eq.1) then
write (6,1538) diam, pres(1)/1.0d+06
1538 format (/,10x,'assume constant conduit diameter:',/,
* 10x,'diameter = ',f6.3,' meters',/,
* 10x,'input pressure = ',f6.2,' MPa')
else
write (6,1539) pgrad, diam
1539 format (/,10x,'assume constant pressure gradient',/,
* 10x,'pressure gradient = ',f7.2,' MPa/km',/,
* 10x,'initial conduit diameter = ',f6.3,' meters')
end if
if (ivt.eq.1) then
write (6,1534)
1534 format (/,10x,'no velocity adjustment',/)
else
write (6,1535)
1535 format (/,10x,'automatic velocity adjustment',/)
end if
if (ives.eq.1) then
write (6,1536)
1536 format (/,10x,
* 'continued exsolution after fragmentation',/)
else
write (6,1537)
1537 format (10x,'no exsolution after fragmentation',/)
end if
C CALCULATE PARAMETERS AT BOTTOM OF CONDUIT
50 continue
iruns = iruns+1
write (6,1433) iruns
1433 format (' STARTING RUN NUMBER ',i3,':',$)
p = pres(1)
dz = -z(1)/100.
wtmin = 5000.
temp = t1 + 2.7315d+02
C CALCULATE AND WRITE OUT MASS FLUX
call exsolv(p)
call dens(p)
C STOP PROGRAM IF MAGMA PRESSURE GRADIENT > CONDUIT PRESSURE GRADIENT
if ((icalc.eq.2).and.(rhomix.gt.(pgrad*100.))) then
write (6,1552) rhomix, rhomix/100.
1552 format (/,10x,
* 'Density of magma/gas mixture = ',f6.0,' kg/m3.',/,
* 10x,'Thus its pressure gradient is ',f6.0,' MPa/km.',/,
* 10x,'This is greater than that specified for the conduit',
* /,10x,
* 'It must be LESS THAN that of the conduit',
* /,10x,
* 'or else magma will not erupt.',//,
* 10x,'program stopped.',/)
stop
end if
mdot = rhomix*xsarea*vel(1)
write (6,135) mdot
135 format(10x,' mass flux=',e12.4,' kg/s')
C CALCULATE AND WRITE INITIAL VALUES OF VARIABLES. IF ICALC=2,
C ADJUST DZ VALUE SO THAT XSAREA DOESN'T CHANGE BY MORE THAN
C 10% BETWEEN STEPS.
if (icalc.eq.1) then
call derivs(z(1),pres(1),dpdz)
else
call derivs(z(1),xsarea,dadz)
dz = abs(0.1*xsarea/dadz)
end if
ho = temp*(mf*cp+mm*cm)+pres(1)/rhom + vel(1)**2/2.
area(1) = xsarea
mach(1) = m
re(1) = rey
rho(1) = rhomix
time(1) = 0.
vesic(1) = vfgas
visc(1) = eta
C ASSIGN VARIABLES TO OUTPUT ARRAY
output(1,1) = area(1)
output(2,1) = mach(1)
output(3,1) = pres(1)/1.e+06
output(4,1) = log10(re(1))
output(5,1) = rho(1)
output(6,1) = time(1)
output(7,1) = vel(1)
output(8,1) = vesic(1)
output(9,1) = log10(visc(1))
output(10,1) = z(1)
output(11,1) = dadz
output(12,1) = dlog10(pres(1))-6.0
output(13,1) = dpdz
output(14,1) = dz
output(15,1) = exco2
output(16,1) = exh2o
output(17,1) = exsulfur
output(18,1) = f
output(19,1) = gamma
output(20,1) = mf
output(21,1) = mm
output(22,1) = r
output(23,1) = rhof
output(24,1) = sv
output(25,1) = xco2
output(26,1) = xh2o
output(27,1) = xsulfur
output(28,1) = cph2o
output(29,1) = cpco2
output(30,1) = cps
output(31,1) = cvh2o
output(32,1) = cvco2
output(33,1) = cvs
output(34,1) = temp-273.15
output(35,1) = ho/1.0d+03
output(36,1) = cp
output(37,1) = diam/2.
write (6,13) ' i', (title(ivar(ix)), ix=1,ivars)
write (6,1402) i, (output(ivar(ix),1), ix=1,ivars)
C**********************************************************************
C BEGIN LOOP
i = 1
100 continue
zin = z(i)
p = pres(i)
C CALL SELF-ADJUSTING RUNGE-KUTTA PROGRAM TO CALCULATE NEXT Z STEP
C AND PRSSURE ETC AT THAT STEP.
if (icalc.eq.1) then
call rkqc(p,dpdz,zin,dz,eps,p,dzout,dznext,derivs)
i = i+1
z(i) = zin
pres(i) = p
else
call rkqc(xsarea,dadz,zin,dz,eps,xsarea,dzout,dznext,derivs)
i = i+1
z(i) = zin
pres(i) = 1.013d+05 + dpdz*zin
end if
if ((i.eq.1).and.(dzout.ne.dz)) output(14,1) = dzout
dz = dznext
C IF Z>0., ADJUST DZ AND CALL RK4 TO CALCULATE P AT Z=0.
if (z(i).gt.0.) then
dz=-z(i-1)
if (icalc.eq.1) then
call rk4(pres(i-1),dpdz,z(i-1),dz,pout,derivs)
pres(i) = pout
else
call rk4(area(i-1),dadz,z(i-1),dz,xsarea,derivs)
pres(i) = 1.013d+05
area(i) = xsarea
end if
z(i) = 0.
end if
C ASSIGN NEW VALUES OF PARAMETERS TO ARRAYS
if (icalc.eq.1) then
call derivs(z(i),pres(i),dpdz)
else
call derivs(z(i),xsarea,dadz)
end if
area(i) = xsarea
mach(i) = m
re(i) = rey
rho(i) = rhomix
vel(i) = v
vavg = (vel(i)+vel(i-1))/2.
time(i) = time(i-1)+(z(i)-z(i-1))/vavg
vesic(i) = vfgas
visc(i) = eta
C CALCULATE NEW TEMPERATURE AND ENTHALPY BY ASSUMING THAT THE SUM
C (ENTHALPY + SP. KINETIC ENERGY OF MIXTURE) IS CONSTANT
h = ho - vel(i)**2/2. - g*(z(i)-z(1))
temp = (h - pres(i)/rhom)/(mf*cp + mm*cm)
C WRITE VALUES TO ARRAY "OUTPUT"
output(1,i) = area(i)
output(2,i) = mach(i)
output(3,i) = pres(i)/1.e+06
output(4,i) = log10(re(i))
output(5,i) = rho(i)
output(6,i) = time(i)
output(7,i) = vel(i)
output(8,i) = vesic(i)
output(9,i) = log10(visc(i))
output(10,i) = z(i)
output(11,i) = dadz
output(12,i) = dlog10(pres(i))-6.0
output(13,i) = dpdz
output(14,i) = dz
output(15,i) = exco2
output(16,i) = exh2o
output(17,i) = exsulfur
output(18,i) = f
output(19,i) = gamma
output(20,i) = mf
output(21,i) = mm
output(22,i) = r
output(23,i) = rhof
output(24,i) = sv
output(25,i) = xco2
output(26,i) = xh2o
output(27,i) = xsulfur
output(28,i) = cph2o
output(29,i) = cpco2
output(30,i) = cps
output(31,i) = cvh2o
output(32,i) = cvco2
output(33,i) = cvs
output(34,i) = temp-273.15
output(35,i) = h/1.0d+03
output(36,i) = cp
output(37,i) = diam/2.
C CALCULATE WATER TABLE DEPTH REQUIRED FOR HYDROSTATIC PRESSURE
C TO EXCEED PRES(I) (ASSUME WATER DENSITY OF 0.9, OR 9000 PA/M)
wtdepth = z(i) + (pres(i)-1.013d+05)/9000.
if (wtdepth.lt.wtmin) wtmin=wtdepth
if (icalc.eq.2) go to 1500
C IF ICALC=1 AND P<1 ATM, OR M>1, ADJUST THE INPUT VELOCITY
if (((pres(i).lt.1.013e+05).or.(mach(i).gt.1.00)).
* and.(z(i).lt.-0.05)) then
if (ivt.eq.1) then
1491 if (pres(i).lt.1.013e+05) then
write (6,533)
533 format (/,5x,'p < 1 atm. program stopped')
else
write (6,1434)
1434 format (/,5x,'mach number > 1. program stopped')
end if
go to 1501
else
1492 write (6,1402) i, (output(ivar(ix),i), ix=1,ivars)
if ((vel(1).eq.0.001).and.(pres(i).lt.1.013e+05)) then
write (6,1445)
1445 format(/,10x,'pressure insufficient to produce eruption',/)
iend = 1
go to 1501
end if
call adjust
if (iend.eq.1) go to 1501
go to 50
end if
end if
C IF M>0.99999 AND Z>-0.05, CALL IT GOOD
if ((z(i).gt.-0.05).and.(mach(i).gt.0.99999)) go to 1501
C RETURN TO BEGINNING OF LOOP IF Z < 0. AND P > 1 ATM
1500 if ((z(i).lt.0.).and.(i.lt.imax)) go to 100
C IF imax ITERATIONS IS REACHED, OR IVT=1, STOP PROGRAM
if ((i.eq.imax).or.(ivt.eq.1)) go to 1501
C IF Z=0. AND P>1 ATM, ADJUST THE INPUT VELOCITY AND START OVER
if (pres(i).gt.1.014e+05) then
write (6,1402) i, (output(ivar(ix),i), ix=1,ivars)
call adjust
if (iend.eq.1) then
write (6,5109)
5109 format (/,10x,'limit of resolution reached',/)
go to 1501
end if
go to 50
C IF PRESSURE IS BETWEEN 0.1012 AND 0.1014 MPA, CALL IT GOOD
else if ((pres(i).le.1.014e+05).and.(pres(i).gt.1.012e+05))
* then
go to 1501
C IF PRESSURE IS LESS THAN .1012 MPA, ADJUST THE INPUT VELOCITY
else
write (6,1402) i, (output(ivar(ix),i), ix=1,ivars)
write (6,1424)
1424 format (/,10x,'pressure < 1 atm',/)
call adjust
if (iend.eq.1) go to 1501
go to 50
end if
C***********************************************************************
C ONCE THE A SOLUTION HAS BEEN REACHED, WRITE OUT FINAL RESULTS
1501 continue
ifinal = i
C WRITE OUT HEADERS
write (9,1301) ifinal,' i', (title(ivar(ix)), ix=1,ivars)
13 format (/,1x,a4,7a10)
1301 format (i4,a3,7a10)
C WRITE FINAL VALUES TO STANDARD OUTPUT
if (iend.ne.1) then
write (6,1402) ifinal, (output(ivar(ix),ifinal), ix=1,ivars)
if ((ifinal.eq.imax).and.(icalc.eq.2)) then
write (6,1902)
1902 format (/,10x,'Number of iterations has exceeded the limit.',/,
* 10x,'You should be able to perform this calculation in',/,
* 10x,'fewer iterations. Check and see if the radius is',/,
* 10x,'changing rapidly at the bottom of the conduit.',/,
* 10x,'(if icalc=2). If so, try adjusting these',/,
* 10x,'parameters and running it over again.')
go to 1034
end if
write (6,1401)
1401 format (/,' successful completion')
C CALCULATE MAXIMUM THEORETICAL VELOCITY, FINAL TEMPERATURES
call velmax
C WRITE OUT WATER TABLE DEPTH REQUIRED FOR WATER INFLUX
write (6,1042) wtmin
1042 format (/,' maximum water table depth that will',
* ' allow g.w. influx = ',f8.2,' meters',/,
* ' (negative values are below ground surface,',
* ' positive values are above)',/)
end if
C WRITE FLOW PROPERTIES TO OUTPUT FILE
1034 do 1100 in=1,ifinal
write (9,1403) in, (output(ivar(ix),in), ix=1,ivars)
1402 format (1x,i4,7f10.3)
1403 format (3x,i4,7f10.3)
1100 continue
close(8)
close(9)
end
C END OF PROGRAM
C************************************************************************
C************************************************************************
C SUBROUTINE RKQC
C SUBROUTINE TAKEN FROM NUMERICAL RECIPES, P. 558, THAT
C CALCULATES THE PRESSURE AT THE NEW VALUE OF Z AND ADJUSTS
C THE NEXT DZ VALUE FOR THE CHANGING CONDITIONS
subroutine rkqc(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
implicit double precision (a-h,o-z)
parameter (pgrow=-0.20,pshrink=-0.25,fcor=1./15.,
* one=1.,safety=0.9,errcon=6.d-04)
real*8 dydx,eps,errmax,h,hdid,hh,hnext,
* htry,x,xsav,y,ysav,
* yscal,ytemp
external derivs
xsav = x
ysav = y
dysav = dydx
h = htry
1 hh = 0.5*h
call rk4(ysav,dysav,xsav,hh,ytemp,derivs)
x = xsav+hh
call derivs(x,ytemp,dydx)
call rk4(ytemp,dydx,x,hh,y,derivs)
x = xsav+h
if (x.eq.xsav) then
write (6,1400)
1400 format (/,10x,'Stepsize not significant in rkqc.',/,
* 10x,'If you are using a constant pressure gradient',/,
* 10x,'chances are that your conduit radius is too large',/,
* 10x,'(or too small),or your entrance velocity is too low',/,
* 10x,'(or too high). Try again.',//,
* 10x,'program stopped',/)
stop
end if
call rk4(ysav,dysav,xsav,h,ytemp,derivs)
ytemp = y-ytemp
errmax = (abs(ytemp/yscal))/eps
if (errmax.gt.one) then
h = safety*h*(errmax**pshrink)
go to 1
else
hdid=h
if (errmax.gt.errcon) then
hnext = safety*h*(errmax**pgrow)
else
hnext = 1.7*h
end if
end if
y = y + ytemp*fcor
return
end
C************************************************************************
C SUBROUTINE RK4
C SUBROUTINE TO CALCULATE NEW PRESSURE AT THE NEW ELEVATION
C USING THE FOURTH-ORDER RUNGE-KUTTA METHOD. THIS SUBROUTINE
C WAS MODIFIED FROM THAT IN NUMERICAL RECIPES, P. 553.
subroutine rk4(y,dydx,x,h,yout,derivs)
implicit double precision (a-h,o-z)
real*8 dydx,dydx2,dym,dyt,h,h6,hh,x,xh,y,yout,yt
hh = h*0.5
h6 = h/6.
xh = x + hh
dydx2 = dydx
yt = y + hh*dydx
call derivs(xh,yt,dyt)
yt = y + hh*dyt
call derivs(xh,yt,dym)
yt = y + h*dym
dym = dyt + dym
call derivs(x+h,yt,dyt)
dydx = dydx2
yout = y + h6*(dydx + dyt + 2.*dym)
return
end
C************************************************************************
subroutine derivs(znow,yp,grad)
C SUBROUTINE TO CALCULATE PRESSURE GRADIENT BETWEEN THE PREVIOUS
C AND THE CURRENT DEPTH.
C DECLARE VARIABLES USED IN THIS SUBROUTINE
implicit double precision (a-h,o-z)
real*8 area(2000),dadz,diam,eta,f,fo,g,grad,m,
* mdot,p,rey,rhomix,sv,v,vfgas,xsarea,yp,z(2000),znow
integer i, icalc
C DECLARE VARIABLES LISTED IN THE COMMON BLOCK BUT NOT USED
C IN THIS SUBROUTINE
real*8 blkmag,cm,co2,cp,cph2o,cpco2,cps,cvh2o,cvco2,cv,cvs,eps,
* exco2,exh2o,exmol,exsulfur,gamma,h2o,
* mach(2000),mco2,mf,mh2o,mm,ms,pres(2000),r,re(2000),
* rho(2000),rhof,rhom,sulfur,temp,vel(2000),vesic(2000),
* visc(2000),xco2,xh2o,xsulfur
integer iend,im,ip,io,ives
C COMMON BLOCKS
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
common/ints/i,icalc,iend,im,ip,io,ives
if (icalc.eq.1) then
p = yp
else
xsarea = yp
if (xsarea.lt.0.) pause ' xsarea < 0'
diam = 2.*sqrt(xsarea/3.14159)
p = 1.013d+05 + dpdz*znow
end if
C CALCULATE MIXTURE DENSITY
call dens(p)
C CALCULATE VISCOSITY
call vis
C CALCULATE VELOCITY, REYNOLDS, NUMBER, FRICTION FACTOR
v = mdot/(rhomix*xsarea)
rey = rhomix*v*diam/eta
f = 16./rey + fo
C CALCULATE SONIC VELOCITY, MACH NUMBER
call svel
m = v/sv
C CALCULATE GRADIENT IN PRESSURE OR XSAREA
if (icalc.eq.1) then
term1 = 2.*f*v**2/diam
grad = (rhomix*(-g - 2.*f*v**2/diam + v**2/xsarea * dadz))/
* (1.-m**2)
dpdz = grad
else
grad = xsarea/v**2 * ( dpdz*(1.-m**2)/rhomix + g +
* (2.*f*v**2)/diam )
dadz = grad
end if
return
end
C*******************************************************************
C SUBROUTINE DENS
C SUBROUTINE TO CALCULATE MIXTURE DENSITY. IT TAKES THE
C PRESSURE P FROM THE CALLING PROGRAM AND RETURNS THE VALUE
C RHOMIX
subroutine dens(p)
C DECLARE VARIABLES USED IN THIS SUBROUTINE
implicit double precision (a-h,o-z)
real*8 exmol, mf, mm, p, r, rhof, rhom, rhomix, temp
integer ives
C DECLARE VARIABLES LISTED IN THE COMMON BLOCK BUT NOT USED
C IN THIS SUBROUTINE
real*8 area(2000),blkmag,cm,co2,cp,cpco2,cph2o,cps,
* cv,cvco2,cvh2o,cvs,dadz,diam,eps,eta,exco2,exh2o,
* exsulfur,f,fo,g,gamma,h2o,m,mach(2000),mco2,mdot,
* mh2o,ms,pres(2000),re(2000),rey,rho(2000),sulfur,
* sv,v,vel(2000),vesic(2000),
* vfgas,visc(2000),xsarea,z(2000)
integer i,iend,im,io,ip
C COMMON BLOCKS
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
common/ints/i,icalc,iend,im,ip,io,ives
C DETERMINE NEW MASS FRACTION GAS ONLY IF MAGMA IS UNFRAGMENTED.
C IF MAGMA IS FRAGMENTED, ASSUME THAT NO NEW GAS EXSOLVES.
if ((vfgas.lt.0.75).or.(ives.eq.1)) then
call exsolv(p)
end if
C CALCULATE GAS DENSITY, MIXTURE DENSITY, VOLUME FRACTION GAS
if (mf.gt.0.) then
rhof = mf*p/(exmol*r*temp)
rhomix = 1./(mm/rhom + mf/rhof)
vfgas = mf*rhomix/rhof
else
rhof = 1.
rhomix = rhom
vfgas = 0.
end if
return
end
C****************************************************************
C SUBROUTINE EXSOLV
C SUBROUTINE TO CALCULATE MASS FRACTION OF EXSOLVED GAS
C USING EQUATIONS TAKEN FROM GERLACH (1986)
subroutine exsolv(p)
implicit double precision (a-h,o-z)
C DECLARE VARIABLES LISTED IN THE COMMON BLOCK BUT NOT USED
C IN THIS SUBROUTINE
real*8 area(2000),blkmag,cpco2,dadz,diam,eta,eps,
* f,fo,g,m,mach(2000),mdot,pres(2000),
* re(2000),rey,
* rho(2000),rhof,rhom,rhomix,sv,
* temp,v,vel(2000),vesic(2000),
* visc(2000),xsarea,z(2000)
C DECLARE VARIABLES USED IN THIS SUBROUTINE
real*8 a1,a2,b1,b2,cm,co2,cp,cph2o,cps,cvco2,cvh2o,cv,cvs,dco2,
* dh2o, dph2o1, dph2o2, dsulfur,
* exco2, exdif, exh2o, exh2o1, exh2o2, exmol, exsulfur,
* extotal,gamma,h2o,mco2,mf,mh2o,mm,ms,p,ph2o,ph2onew,
* r,sulfur,vfgas, vgvm, xco2, xh2o, xsulfur
C COMMON BLOCKS
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
if (p.gt.0.) then
C CALCULATE DISSOLVED CO2 USING EQUATION 7 FROM GERLACH (1986)
dco2 = 0.0005 + 5.9e-04*(p/1.0d+06)
if (dco2.ge.co2) then
mf = 0.
vfgas = 0.
gamma = 1.167
exco2 = 0.
exh2o = 0.
exsulfur = 0.
xco2 = 0.333
xh2o = 0.333
xsulfur = 0.333
extotal = 0.
exmol = 0.
go to 300
end if
exco2 = co2 - dco2
C BEGIN LOOP THAT ITERATES TO FIND THE PARTIAL PRESSURE OF H2O
C (PH2O, WHICH IS NEEDED TO CALCULATE DISSOLVED H2O). START
C BY GUESSING THAT PH2O = P/2.
ph2onew = p/2.
1000 continue
ph2o = ph2onew
C CALCULATE EXSOLVED WATER USING EQUATION 6 FROM GERLACH (1986)
dh2o = 1801./(8394.19*sqrt(p/1.0d+06)*
* ((ph2o/1.0d+06)**(-.9917)) - 356.98)
exh2o1 = h2o - dh2o
C CALCULATE EXSOLVED WATER USING A REARRANGED VERSION OF EQUATION
C 13 FROM GERLACH (1986)
exh2o2 = (21.07*ph2o*exco2)/(44.*p - 51.48*ph2o)
C COMPARE EXSOLVED WATER FROM THE TWO METHODS
exdif = exh2o2-exh2o1
if (abs(exdif).gt.0.00001) then
C IF THEY DON'T AGREE, USE THE FIRST DERIVATIVES OF EXH2O1 AND EXH2O2
C (I.E. D(EXH2O1)/D(PH2O)) AND D(EXH2O2)/D(PH2O) TO CALCULATE A NEW
C PH2O
dph2o1 = (-1.5001e+07*sqrt(p/1.0d+06)*((ph2o/1.0d+06)**
* (-1.9917))) /
* (8394.19*sqrt(p/1.0d+06)*((ph2o/1.0d+06)**(-0.9917))
* - 356.98)**2
dph2o1 = dph2o1/1.0d+06
dph2o2 = ((21.06*exco2)/(44.*p-51.48*ph2o)) *
* ( 1. + (51.48*ph2o)/(44.*p - 51.48*ph2o))
C THESE DERIVATIVES MAKE UP THE SLOPES OF LINES THAT GO THROUGH PH2O AND
C THE RESPECTIVE EXH2O VALUES. THE LINES ARE DEFINED BY THE EQUATIONS
C Y1 = A1 * X + B1 FOR EXH2O1 VS PH2O, AND
C Y2 = A2 * X + B2 FOR EXH2O2 VS PH2O.
C THE OBJECTIVE IS TO FIND WHERE THESE LINES CROSS. THAT WILL BE OUR
C NEW ESTIMATED VALUE OF PH2O.
a1 = dph2o1
a2 = dph2o2
b1 = exh2o1 - ph2o * a1
b2 = exh2o2 - ph2o * a2
C THE LINES CROSS AT THE VALUE PH2O = (B1-B2)/(A2-A1).
ph2onew = (b1-b2)/(a2-a1)
C MAKE SURE THE NEW PH2O DOES NOT EXCEED 0.847*P (IN WHICH CASE THE
C SOLUTION BLOWS UP).
if (ph2onew.gt.(0.847*p)) then
ph2onew = (ph2o + 0.847*p)/2.
end if
go to 1000
end if
C ONCE WE HAVE FOUND THE VALUE OF PH2O, CALCULATE TOTAL EXSOLVED
C FLUIDS AND VESICULARITY
exh2o = (exh2o1+exh2o2)/2.
if (exh2o.lt.0.) exh2o = 0.
dh2o = h2o-exh2o
exsulfur = 0.3025*exh2o + 0.1238*exco2
if (exsulfur.lt.0.) exsulfur = 0.
dsulfur = sulfur-exsulfur
extotal = exh2o + exco2 + exsulfur
exmol = .01*(exh2o/mh2o + exco2/mco2 + exsulfur/ms)
mf = extotal/100.
mm = 1.-mf
xh2o = ph2o/p
xco2 = 0.855 - xh2o
xsulfur = 0.155
vgvm = exmol*r*temp*rhom/(p*mm)
vfgas = vgvm/(vgvm+1.)
C CALCULATE CP FOR GAS SPECIES USING
C EMPIRICAL FORMULAS TAKEN FROM TABLE A-15 IN MORAN IN SHAPIRO
C AT VARIOUS TEMPERATUES (P ASSUMED =1 ATM)
cph2o = (4.070 - 1.108d-3*temp + 4.152e-6*temp**2 -
* 2.964d-9*temp**3 + 0.807d-12*temp**4) * r/mh2o
cpco2 = (2.401 + 8.735d-3*temp - 6.607d-6*temp**2 +
* 2.002d-9*temp**3) * r/mco2
cps = (3.267 + 5.324d-3*temp + 0.684d-06*temp**2 -
* 5.281d-9*temp**3 + 2.559d-12*temp**4) * r/ms
cp = (exh2o*cph2o+exco2*cpco2+exsulfur*cps)/extotal
cvh2o = cph2o - 8.314/mh2o
cvco2 = cpco2 - 8.314/mco2
cvs = cps - 8.314/ms
C CALCULATE CV, GAMMA
cv = (exh2o*cvh2o+exco2*cvco2+exsulfur*cvs)/extotal
gamma = cp/cv
C IF PRESSURE < 0., MAKE UP NUMBERS FOR GAMMA, MF, VFGAS
else
mf = (h2o+co2+sulfur)/100.
gamma = 1.25
vfgas = 0.99
end if
300 continue
mm = 1. - mf
return
end
C****************************************************************
C SUBROUTINE SVEL
C SUBROUTINE TO CALCULATE SONIC VELOCITY
subroutine svel
C DECLARE VARIABLES USED IN THIS SUBROUTINE
implicit double precision (a-h,o-z)
real*8 blkgas, blkmag, blkmix, gamma, mf, mm, r,
* rhof, rhom, rhomix, sv, temp, vfgas, vfmag
C DECLARE VARIABLES LISTED IN THE COMMON BLOCK BUT NOT USED
C IN THIS SUBROUTINE
real*8 area(2000),cm,co2,cp,cpco2,cph2o,cps,cvco2,cvh2o,
* cv,cvs,dadz, diam,eps,eta,exco2,exh2o,exmol,exsulfur,
* f,fo,g,h2o,m,mach(2000),mco2,mdot,mh2o,ms,pres(2000),
* re(2000),rey,rho(2000),sulfur,v,vel(2000),
* vesic(2000),visc(2000),xsarea,z(2000)
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
C CALCULATE VOLUME FRACTION OF ROCK
vfmag = 1.-vfgas
C CALCULATE BULK MODULI AND SONIC VELOCITY
if (mf.gt.0.) then
blkgas = rhof*gamma*exmol*r*temp/mf
blkmix = 1./(vfgas/blkgas + vfmag/blkmag)
sv = sqrt(blkmix/rhomix)
else
sv = sqrt(blkmag/rhom)
end if
return
end
C*******************************************************************
C SUBROUTINE VIS
C SUBROUTINE THAT CALCULATES VISCOSITY AS A FUNCTION OF T,
C GAS CONTENT, AND VESICULARITY. THE SECTION THAT GIVES
C VISCOSITY CONTENT VERSUS T AND WATER CONTENT WAS TAKEN
C FROM SHAW (1972) FOR KILAUEAN BASALTS.
C THE SECTION THAT GIVES VISCOSITY AS A
C FUNCTION OF VESICULARITY IS MODIFIED FROM EQUATIONS IN
C DOBRAN (1992, EQ. 24)
subroutine vis
C DECLARE VARIABLES USED IN THIS SUBROUTINE
implicit double precision (a-h,o-z)
real*8 eta, etagas, etagaslg, etamag, etamaglg, etamixlg,
* exp, temp, vfgas
C DECLARE VARIABLES THAT ARE IN THE COMMON BLOCK BUT NOT USED
C IN THIS SUBROUTINE
real*8 area(2000),blkmag,cm,co2,cp,cpco2,cph2o,cps,cv,
* cvco2,cvh2o,cvs,dadz,diam,eps,exco2,exh2o,exmol,
* exsulfur,f,fo,g,gamma,h2o,m,mach(2000),
* mco2,mdot,mf,mh2o,mm,ms,pres(2000),
* r,re(2000),rey,rho(2000),rhof,rhom,rhomix,sulfur,sv,v,
* vel(2000),vesic(2000),visc(2000),xsarea,z(2000)
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
C CALCULATE MAGMA VISCOSITY AS A FUNCTION OF TEMPERATURE
exp = -10.737 + 1.8183 * (1.0d+04/temp)
eta = 10.**exp
C CALCULATE VISCOSITY AS A FUNCTION OF VFGAS, USING RELATIONSHIPS
C FROM DOBRAN (1992). IF VESICULARITY IS > 50%, USE A WEIGHTING
C FUNCTION TO CALCULATE THE VISCOSITY OF MAGMA AND GAS NEAR THE
C FRAGMENTATION POINT.
if (vfgas.lt.0.5) then
eta = eta / (1.-vfgas)
else if (vfgas.lt.1.0) then
etamag = eta / (1.-vfgas)
etagas = (5.3d-05 * (1.-(1.-vfgas)/0.62)**(-1.5))
etamaglg = dlog10(etamag)
etagaslg = dlog10(etagas)
etamixlg = etamaglg * 2.**(-((vfgas/0.75)**40))
* + etagaslg * 2.**(-((0.75/vfgas)**40))
eta = 10.**(etamixlg)
else
eta = 5.3d-05
end if
return
end
C******************************************************************************
C SUBROUTINE ADJUST
C SUBROUTINE THAT ADJUSTS INPUT VELOCITY
subroutine adjust
C DECLARE VARIABLES USED IN THIS SUBROUTINE
implicit double precision (a-h,o-z)
real*8 eps,int, mach(2000),pres(2000),slope,vel(2000),v1m(100),
* v1o(100),v1p(100),rmaxo,rminv,z(2000),zlm(100),zlp(100)
integer i,im,ip
C DECLARE VARIABLES IN THE COMMON BLOCK BUT NOT USED IN THIS
C SUBROUTINE
real*8 area(2000),blkmag,cm,co2,cp,cpco2,cph2o,cps,cv,
* cvco2,cvh2o,cvs,dadz,diam,eta,exco2,exh2o,exmol,
* exsulfur,f,fo,g,gamma,h2o,m,mco2,mdot,mf,mh2o,
* mm,ms,r,re(2000),
* rey,rho(2000),rhof,rhom,rhomix,sulfur,sv,temp,v,vfgas,
* vesic(2000),visc(2000),xsarea
integer icalc,iend,io,ives
C COMMON BLOCKS
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
common/ints/i,icalc,iend,im,ip,io,ives
C*********************************************************************
C SOLUTIONS FOR Z=0.
C*********************************************************************
if (z(i).eq.0.) then
C IF Z(I)=0 AND P>1 ATM, SEE IF THERE ARE PREVIOUS CALCULATIONS
C RESULTING IN M>1 OR P<1 ATM. IF SO, INCREASE VEL(1) BY 0.2
C TIMES THE DIFFERENCE BETWEEN THE MINIMUM V1M OR V1P AND THE
C MAXIMUM V1O
if (pres(i).gt.1.014e+05) then
write (6,*) 'exit pressure > 1 atm and M < 1'
v1o(io) = vel(1)
io = io+1
C FIND MINIMUM OF V1P AND V1M
rminv = 5.d+05
do 150 j=1,im-1
150 rminv=min(rminv,v1m(j))
do 160 j=1,ip-1
160 rminv=min(rminv,v1p(j))
C FIND MAXIMUM OF V1O
rmaxo = 0.
do 170 j=1,io-1
170 rmaxo=max(rmaxo,v1o(j))
C IF THE MINIMUM OF V1P OR V1M AND THE MAXIMUM OF V1O DON'T DIFFER BY MORE
C THAN ABOUT 100*EPS, THEN CALL THE SOLUTION GOOD.
if ((rmaxo.gt.0.).and.((rminv-rmaxo)/rmaxo).lt.(10.*eps)) then
iend=1
return
C CALCULATE NEW VALUE OF VEL(1)
else if (rminv.lt.5.d+05) then
vel(1) = rmaxo + 0.25*(rminv-rmaxo)
else
vel(1) = 1.5*vel(1)
end if
go to 500
C IF Z=0. AND THE FINAL PRESSURE IS <0.1012 MPA, ADJUST THE
C INITIAL VELOCITY SO THAT IT'S HALFWAY BETWEEN THE CURRENT
C VEL(1) AND THE MAXIMUM V1O
else
rmaxo = 0.
do 171 j=1,io-1
171 rmaxo=max(rmaxo,v1o(j))
vel(1) = rmaxo + 0.5*(vel(1)-rmaxo)
if (vel(1).lt.0.001) then
write (6,143)
143 format (/,10x,'pressure is too low to erupt')
iend = 1
return
end if
go to 500
end if
end if
C*********************************************************************
C SOLUTIONS FOR Z<0.
C*********************************************************************
C IF M>1, THE RELATIONSHIP BETWEEN Z(I) AND VEL(1) SHOULD BE
C ALMOST LINEAR. THEREFORE, IF WE HAVE RUN AT LEAST
C TWO ITERATIONS WHERE VEL(1) AND Z(I) HAVE BEEN OBTAINED,
C WE CAN USE THOSE VALUES TO DEFINE THE LINE RELATING Z(I) AND
C VEL(1), AND CALCULATE THE VEL(1) NEEDED TO OBTAIN M=1 AT
C Z(I) = 0. THIS VALUE IS VEL(1)=-INT/SLOPE, WHERE INT IS THE
C Z(I) INTERCEPT OF THE Z(I) VS VEL(1) LINE, AND SLOPE IS THE
C SLOPE OF THE LINE.
C IT TURNS OUT THAT, IF ONE ASSUMES THE RELATIONSHIP
C BETWEEN Z(I) AND VEL(1) TO BE PERFECTLY LINEAR, ONE OVERSHOOTS
C AND ENDS UP BACK IN THE M<1 FIELD. THEREFORE, I ADD 10% OF
C THE DIFFERENCE BETWEEN THE OLD VEL(1) AND THAT OBTAINED FROM
C VEL(1)=-INT/SLOPE, TO MAKE SURE THAT VEL(1) STAYS WITHIN
C THE M>=1 FIELD.
if (mach(i).gt.1.000) then
write (6,*) ' mach number > 1. adjusting vi'
zlm(im) = z(i)
v1m(im) = vel(1)
C CALCULATE SLOPE AND INTERCEPT OF V(1) VS Z(I) LINE,
C AND VALUE OF V(1) ALONG THIS LINE WHEN Z(I)=0.
if (im.ge.2) then
slope = (zlm(im)-zlm(im-1))/(v1m(im)-v1m(im-1))
int = zlm(im)-slope*v1m(im)
vel(1) = -int/slope + 0.1*(min(v1m(im),v1m(im-1))+
* int/slope)
else
vel(1) = -vel(1)*(z(i)-z(1))/z(1)
end if
im = im+1
C FOR P<1 ATM, ADJUST V(1). ASSUME, AS WHEN M>1, THAT THE
C RELATIONSHIP BETWEEN V(1) AND Z(I) IS LINEAR. CALCULATE V(1)
C ALONG THIS LINE FOR Z(I)=0.
else if (pres(i).lt.1.013d+05) then
write (6,*) 'pressure is less than 1 atm. adjusting vi'
zlp(ip) = z(i)
v1p(ip) = vel(1)
if (ip.ge.2) then
slope = (zlp(ip)-zlp(ip-1))/(v1p(ip)-v1p(ip-1))
int = zlp(ip)-slope*v1p(ip)
if (((-int/slope).le.0.).or.((-int/slope).gt.
* (vel(1)))) then
vel(1) = vel(1)/2.
else
vel(1) = -int/slope
end if
else
vel(1) = -vel(1)*(z(i)-z(1))/z(1)
end if
if (vel(1).le.0.001) vel(1)=0.001
ip = ip+1
end if
C RE-INITIALIZE VARIABLES AND WRITE OUT NEW INITIAL VELOCITY
500 continue
do 100 j=1,i
re(j) = 0.
rho(j) = 0.
mach(j) = 0.
vesic(j) = 0.
visc(j) = 0.
100 continue
do 200 j=2,i
vel(j) = 0.
pres(j) = 0.
z(j) = 0.
200 continue
i = 1
write (6,1044) vel(1)
1044 format (/,' trying new input velocity ',f8.5,' m/s',/)
return
end
C***********************************************************************
C SUBROUTINE VELMAX
subroutine velmax
C CALCULATES MAXIMUM THEORETICAL VELOCITY. THIS IS DONE BY FIRST
C CALCULATING THE AMOUNT OF COOLING THAT ACCOMPANIES DECOMPRESSION
C OF THE MIXTURE TO ATMOSPHERIC PRESSURE. USING THIS AMOUNT OF
C COOLING, A NEW ENTHALPY IS CALCULATED. THEN, FROM THE NEW
C ENTHALPY, A MAXIMUM THEORETICAL VELOCITY IS CALCULATED FROM THE
C EQUATION:
C VMAX = SQRT(2.*(H2-H1))
C
C WHERE H1 IS THE ENTHALPY OF THE MIXTURE AT THE VENT EXIT, AND
C H2 IS THE ENTHALPY OF THE MIXTURE ONCE ITS PRESSURE EQUILIBRATES
C WITH ATMOSPHERIC.
C DECLARE VARIABLES USED IN THIS SUBROUTINE
implicit double precision (a-h,o-z)
real*8 cm, cp, deltah, deltat, gmix, mf, mm, pres(2000), r,
* rhom, temp, temp2, temp2c, vel(2000), vmax
integer i
C DECLARE VARIABLES LISTED IN COMMON BLOCK BUT NOT USED IN
C THIS SUBROUTINE
real*8 area(2000),blkmag,co2,cph2o,cpco2,cps,cv,cvh2o,
* cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,exmol,
* exsulfur,f,fo,g,gamma,h2o,m,mach(2000),mco2,
* mdot,mh2o,ms,re(2000),
* rey,rho(2000),rhof,rhomix,sulfur,sv,v,vfgas,
* vesic(2000),visc(2000),xco2,xh2o,xsarea,xsulfur,
* z(2000)
integer icalc,iend,im,ip,io,ives
C COMMON BLOCK
Common/block1/area,blkmag,cm,co2,cp,cph2o,cpco2,cps,
* cv,cvh2o,cvco2,cvs,dadz,diam,dpdz,eps,eta,exco2,exh2o,
* exmol,exsulfur,f,fo,g,gamma,h2o,m,mach,mdot,mf,mm,
* mco2,mh2o,ms,pres,r,re, rey,rho,rhof,rhom,rhomix,
* sulfur,sv,temp,v,vel,vfgas,vesic,visc,xco2,
* xh2o,xsarea,xsulfur,z
common/ints/i,icalc,iend,im,ip,io,ives
C 1. CALCULATE GAMMA FOR MIXTURE:
gmix = (mf*cp + mm*cm)/(mf*cv + mm*cm)
C 2. CALCULATE NEW TEMPERATURE:
temp2 = temp*(0.1013e+06/pres(i))**((gmix-1.)/gmix)
temp2c = temp2-273.15
deltat = temp-temp2
C 3. CALCULATE CHANGE IN ENTHALPY
deltah = mf*cp*deltat + mm*(cm*deltat +
* (pres(i)-0.1013e+06)/rhom)
C 4. CALCULATE MAXIMUM THEORETICAL VELOCITY
if (deltah.ge.0.) then
vmax = vel(i) + sqrt(2.*deltah)
else
vmax = vel(i) - sqrt(-2.*deltah)
end if
write (6,1432) temp2c, deltat, deltah, vmax
1432 format(/,5x,'AFTER ISENTROPIC EQUILIBRATION TO 1 ATM PRESSURE:',/,
* 10x,' final temperature = ',f7.2,' deg. C',/,
* 10x,' temperature change = ',f7.3,' deg. K',/,
* 10x,' enthalpy change = ',e12.4,' J/kg',/,
* 10x,' max. theoretical velocity = ',f9.2,' m/s',/)
return
end