Appendix E—
Heat Flow (Two-Dimensional Diffusion Code)
Time-dependent heat diffusion expressed by Fick's second law of diffusion [Eq. (2-9); Chapter 2] is an important tool for evaluating thermal resource [Eq. (2-6); Chapter 2]. Because the temporal and spatial variation of temperature in a sequence of rocks of varying character is the desired solution, Eq. (2-9) can be written in a two-dimensional, nonlinear Cartesian form using thermal conductivity (kt ) and temperature (T):
where ktx and kty are the thermal conductivities in the x and y directions, respectively. The following computer program uses a two-dimensional, finite-difference numerical solution based on an algorithm given by Harbaugh and Bonham-Carter (1970). It calculates time-dependent heat diffusion for an array of variable rock properties, using an averaging technique to obtain the spatially varying thermal conductivity coefficients.
Important initial conditions include regional thermal gradient, variation of vertical and horizontal conductivities as a function of user-specified stratigraphy, and time and spatial step size. This code is written for interactive use; results are printed out to an external file as well as to the screen.
The program requests
(1) the number of rows and columns in the array that represents a vertical cross section of the desired model,
(2) the spatial and temporal step for the array, and
(3) the number of rock types.
After the numerical representation of the simulated cross-section is entered by rock type, the array is saved to an external file for
access in future runs. By specifying conductivity coefficients for each rock type, the program calculates effective convective coefficients for locations in which vertical conductivities are greater than horizontal ones. Finally, the program calculates whether the specified time step is small enough to result in stable solutions for the desired mesh.
Source Code Listing
The source code, listed below, is written in VAXÔ FORTRAN, which makes use of virtual memory. Implementation of this program on a personal computer requires a FORTRAN compiler, and the source code may require some changes of array sizes, variable declarations, and assignment of external files. Standard code format requires that each statement is indented at least 7 spaces and line continuation characters (denoted here as '&') are typed in the sixth space of each line.
c............program 2DHF.....................
c........... Two dimensional heat diffusion in heterogeneous
c............materials by solution of Fick's second law
c...........of diffusion, utilizing finite difference technique
c...........of Harbaugh and Bonham-Carter,
c..........."Computer Simulation in Geololgy",
c...........Wiley Interscience, 1970, p. 225
common krock(50,50),temp(50,50),tempi(50,50),
&ntempi(50,50),coefx(50,50),coefy(50,50),ntgradi(50),
&ntemp(1000,50,50),l(50),tempnew(50,50),
&ntgrad(1000,50)
write(5,10)
1 format('1','* * * * * * * 2-D HETEROGENEOUS',
&'HEAT FLOW* * * * * * * * * * * * *')
c....Define grid size and geology..................
write(5,10)
10 format(///,' Enter ncol, nrow, dxy (km), and nrock',/,
&'(number of rock types up to 4, excluding the magma,'/,
&'for the computational mesh:',/)
read(5,*) ncol, nrow, ddxy, nrock
dxy = ddxy* 1.0e+05
nj = ncol/2
write(5,15)
15 format(/,' Enter magma temperature and',
&'regional geothermal gradient (deg C/km):',/)
read(5,*) tm, tg
cv = 0.6
c....Read in geologic stucture and rock characteristics..
50 call grid(nrow,ncol,nrock,tm,tg,dxy,ddxy,cv,coef0,
&coef1x,coef1y,coef2x,coef2y,coef3x,coef3y,coef4x,
&coef4y,conv0,conv1,conv2,conv3,conv4)
c.....Stabilize initial heat flow in grid...........
hfb = 1.0
hft = 1.0
stab = 0.0
call stabil (stab,hfb,hft,nrow,ncol,nj,tm,tg,dxy,cv,
&coef1x,coef1y,conv0,conv1,conv2,conv3,conv4)
c.....Set dt for time step...................
90 write(5,92) hfb,hft
92 format(/,' hfb = ',f7.4,'; hft = ',f7.4,//,
&'What is the time-step (yrs)?',/,
&'If time-step is too large, program will go unstable',//)
read(5,*) dt
dt = dt*3.1536e+7
c.....Reset grid ..................
do 94 i=1,nrow
do 93 j=1,ncol
temp(i,j) = tempi(i,j)
if(krock(i,j).eq.0) then
coefx(i,j) = coef0/cv
coefy(i,j) = coef0/cv
end if
93 continue
94 continue
nt = 1
ntt = 0
factor = dt/dxy**2.0
tmax = 0.0
tmin = tm
itt = 1
write(5,95)
95 format(//,' Grid all set up for calculation,',/,
&'do you wish to continue?<y>')
read(5,'(a)') ano 1
if(ano1.eq.'n'.or.ano1.eq.'N') go to 200
write(5,96)
96 format(",//,10x,' nt',5x,' tmax ',5x,' tmin',/)
c....Begin diffusion loop.........
100 call diffus (stab,hfb,hft,ntt,nt,tmax,tmin,nrow,nj,
& ncol,factor,conv0,conv1,conv2,conv3,conv4)
c....Record every 1000*dt calculation of thermal state.....
id = 1.0/ddxy + 0.5
do 170 i=1,nrow
do 160 j=1,ncol
ntemp(nt,i,j) = temp(i,j) + 0.5
if(i.eq.id) ntgrad(nt,j) = (temp(i,j)-20.0) + 0.5
160 continue
170 continue
write(5,175)nt,tmax,tmin
175 format(10x,i3,5x,f6.2,5x,f6.2)
if(nt.eq.1.or.nt.eq.5.or.nt.eq.(itt*10)) then
write(5,180)
read(5,'(a)') ano
if(ano.eq.'y'.or.ano.eq.'Y') go to 90
write(5,181)
read(5,'(a)') ayes1
if(ayes1.eq.'n'.or.ayes1.eq.'N') go to 190
end if
180 format(//,' Do you wish to change time step? <n>')
181 format(' Do you wish to continue? <y>')
190 if(tm.ge.300.0) then
if(tmax.le.300.0) go to 200
if(nt.ge.199.or.ntt.lt.1000) go to 200
end if
if(ntt.eq.1000) ntt = 0
if(ayes1.eq.'n'.or.ayes1.eq.'N') go to 200
nt = nt + 1
if(nt.gt.(itt*10)) itt = itt + 1
go to 100
c....Write out results, timef (ka).................
200 ncy = nt
nccy = ntt
timef = ((ncy*1000.0+nccy)*dt)/3.1536e+07
tntt = (dt*1000.0)/3.1536e+07
tint = 0.0
space = dxy/1.0e+05
time = 1.0
write(5,250) timef,tmax,nt,tntt
250 format('1',///,' Calculation Complete',/,' ',
&'___________________________________________',//,
&' Cooling time = ',8x,f10.1,' years',/,
&' Maximum magma temperature = ',f5.1,' deg C'/,
&' Number of plots = ',i3,'at',f10.2,' year intervals',//,
&' Do you wish printout of initial geometry? <y>',//)
read(5,'(a)') ano
if(ano.eq.'n'.or.ano.eq.'N') goto 320
c....Initial time plots...............
write(5,309) space, tint, (l(j),j=1,ncol)
do 300 i=1,nrow
write(5,311) i, (krock(i,j), j=1,ncol)
300 continue
309 format('1',' 2-D HEAT FLOW PLOT',//,
&' Grid spacing = ',f4.2,' km',
&' Time = ',f12.1,' yrs ',//,5x,<ncol>i4,//)
310 format('1',' 2-D HEAT FLOW PLOT',//,
&' Grid spacing = ',f4.2,' km',
&' Time = ',f12.1,' yrs',//,
&5x,<ncol>i4,/,' d/km',<ncol>i4,//,' 0',2x,<ncol>('20'))
311 format(",i2,2x,<ncol>i4)
320 write(5,330)
330 format(/,' Do you wish printout of init. temperatures?<y>')
read(5,'(a)') ano
if(ano.eq.'n'.or.ano.eq.'N') goto 350
write(5,310) space, tint, (l(j),j=1,ncol),(ntgradi(j),j=1,ncol)
do 340 i=1,nrow
write(5,311) i, (ntempi(i,j),j=1,ncol)
340 continue
c....Calculated-time plots............
350 if(ano1.eq.'n'.or.ano1.eq.'N') go to 402
360 write(5,370)
370 format('For which time interval do you wish a plot?',/,
&' enter 0 for initial time plot, (-1) for none.')
read(5,*)time
if(time.eq.-1.0) go to 400
if(time.eq.0.) go to 320
ncy = time
nccy = 0
if(ncy.eq.nt) nccy = ntt
tint = (ncy*1000.0+nccy)*dt/3.1536e+07
write(5,310) space,tint,(l(j),j=1,ncol),
&(ntgrad(ncy,j),j=1,ncol)
do 380 i=1,nrow
write(5,311)i, (ntemp(ncy,i,j),j=1,ncol)
380 continue
write(5,390)
390 format(",'Do you wish another time plot?<y>')
read(5,'(a)')ano
if(ano.eq.'n'.or.ano.eq.'N') goto 400
go to 360
400 if(ayes1.eq.'n'.or.ayes1.eq.'N') then
402 write(5,405)
read(5,'(a)')ano
if(ano.eq.'n') go to 409
ayes 1 = 'y'
go to 190
end if
405 format(/,'Do you wish to continue this calculation? <y>')
409 write(5,410)
410 format(/,'Do you wish to reset the conductivities and',
&'time step? <n>')
read(5,'(a)')ayes
if(ayes.eq.'y') go to 50
500 stop
end
c.........SUBROUTINE GRID.......................
subroutine grid (nrow,ncol,nrock,tm,tg,dxy,ddxy,cv,
&coef0,coef1x,coef1y,coef2x,coef2y,coef3x,coef3y,
&coef4x,coef4y,conv0,conv1,conv2,conv3,conv4)
common krock(50,50),temp(50,50),tempi(50,50),
&ntempi(50,50),coefx(50,50),coefy(50,50),ntgradi(50),
&ntemp(1000,50,50),l(50),tempnew(50,50),
&ntgrad(1000,50)
open(unit=8,status='unknown',file='2dhdif.dat')
c....Set up initial rock geometry in grid......
write(5,10)
10 format(//' Read rock-grid file? <y>')
read(5,'(a)') ano
if(ano.eq.'n'.or.ano.eq.'N') go to 30
do 20 i=1,nrow
read(8,45) (krock(i,j), j=1,ncol)
20 continue
go to 50
30 write(5,35)
35 format(/,'Enter rock types in grid by row,',/,
&' 0 = magma, 1 for deepest basement rock,',
&'2,3,4 for other rocks:',/)
do 40 i=1,nrow
read(5,*) (krock(i,j), j=1,ncol)
write(8,45) (krock(i,j), j=1,ncol)
40 continue
45 format(<ncol>i3)
c......Enter x and y conductivity constants.......
c......set constant volume heat capacity .........
c......for cp = 0.24 cal/g-deg, rho = 2.5 g/cm3..
c......to cv = 0.60..........................................
50 coef0 = 0.005
coef = coef0*1000.0
conv0 = 3.0
convf = 1000.0
write(5,51) coef
51 format(/,' Enter x and y conductivities (mcal/cm-s-deg)',
&'by rock type.',/,' For veritcal convection y > x,',
&' for horizontal convection (to the right) x > y):'/,
&' Magma conductivity = ',f3.1,' mcal/cm-s-deg'//,
&' Rock 1 conductivities = ')
read(5,*) coef1x,coef1y
conv1 = 1.0
if(coef1y.gt.coef1x) then
write (5,59)
read(5,'(a)') ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv1 = convf*coef1y/coef1x
end if
if(coef1x.gt.coef1y) then
write (5,59)
read(5,'(a)') ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv1 = -convf* coef1x/coef1y
end if
if(nrock.eq.1) go to 60
write(5,53)
53 format(/,' Rock 2 conductivities = ')
read(5,*) coef2x,coef2y
conv2 = 1.0
if(coef2y.gt.coef2x) then
write (5,59)
read(5,'(a)') ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv2 = convf* coef2y/coef2x
end if
if(coef2x.gt.coef2y) then
write (5,59)
read(5,'(a)') ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv2 = -convf* coef2x/coef2y
end if
if(nrock.eq.2) go to 60
write(5,55)
55 format(/,' Rock 3 conductivities = ')
read(5,*) coef3x,coef3y
conv3 = 1.0
if(coef3y.gt.coef3x) then
write (5,59)
read(5,'(a)') ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv3 = convf* coef3y/coef3x
end if
if(coef3x.gt.coef3y) then
write (5,59)
read(5,'(a)')ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv3 =-convf* coef3x/coef3y
end if
if(nrock.eq.3) go to 60
write(5,57)
57 format(/,' Rock 4 conductivities =')
read(5,*) coef4x,coef4y
conv4 = 1.0
if(coef4y.gt.coef4x) then
write (5,59)
read(5,'(a)') ayes
f(ayes.eq.'y'.or.ayes.eq.'Y')
& conv4 = convf*coef4y/coef4x
end if
if(coef4x.gt.coef4y) then
write (5,59)
read(5,'(a)') ayes
if(ayes.eq.'y'.or.ayes.eq.'Y')
& conv4 = -convf*coef4x/coef4y
end if
59 format(/'Is this rock unit convective?<n>')
c....Place initial temperatures in grid.........
60 idd = 1.0/ddxy + 0.5
do 80 i=1, nrow
do 70 j=1,ncol
l(j) = j
if(krock(i,j).eq.0) then
tempi(i,j) = tm
ntempi(i,j) = tempi(i,j) + 0.5
end if
if(krock(i,j).gt.0) then
tempi(i,j) = 20.0 + (tg*i*dxy/1.0e+05)
ntempi(i,j) = tempi(i,j) + 0.5
end if
if(i.eq.idd) ntgradi(j) = (tempi(idd,j)-20.0) + 0.5
70 continue
80 continue
c....Place diffusivity coefficients in grid.........
do 100 i=1, nrow
do 90 j=1,ncol
if(krock(i,j).eq.0) then
coefx(i,j) = coef0/cv
coefy(i,j) = coef0/cv
end if
if(krock(i,j).eq.1) then
coefx(i,j) = 0.001*coef1x/cv
coefy(i,j) = 0.001*coef1y/cv
end if
if(nrock.eq.1) go to 90
if(krock(i,j).eq.2) then
coefx(i,j) = 0.001*coef2x/cv
coefy(i,j) = 0.001*coef2y/cv
end if
if(nrock.eq.2) go to 90
if(krock(i,j).eq.3) then
coefx(i,j) = 0.001*coef3x/cv
coefy(i,j) = 0.001*coef3y/cv
end if
if(nrock.eq.3) go to 90
if(krock(i,j).eq.4) then
coefx(i,j) = 0.001*coef4x/cv
coefy(i,j) = 0.001*coef4y/cv
end if
90 continue
100 continue
close(unit=8,status='keep')
return
end
c..........SUBROUTINE STABIL.......................
Subroutine stabil (stab,hfb,hft,nrow,ncol,nj,tm,tg,dxy,cv,
& coef1x,coef1y,conv0,conv1,conv2,conv3,conv4)
common krock(50,50),temp(50,50),tempi(50,50),
&ntempi(50,50),coefx(50,50),coefy(50,50),ntgradi(50),
&ntemp(1000,50,50),l(50),tempnew(50,50),
&ntgrad(1000,50)
if(stab.eq.1.0) go to 100
10 write(5,10)
10 format(///,' Stabilizing ambient heat flow')
do 30 i=1,nrow
do 20 j=1,ncol
temp(i,j) = 20.0 + (tg*i*dxy/1.0e+05)
if(krock(i,j).eq.0)coefx(i,j) = 0.001*coef1x/cv
coefy(i,j) = coefx(i,j)
20 continue
30 continue
nt = 0
ntt = 0
tob = tempi(nrow,nj)
tot = tempi(1,nj)
tmax = 0.0
tmin = tm
dt = 10.0 * 3.1536e+7 * ((dxy/1.0e+05)**2.0)
factor = dt/dxy**2.0
40 call diffus (stab,hfb,hft,ntt,nt,tmax,tmin,nrow,nj,
& ncol,factor,conv0,conv1,conv2,conv3,conv4)
if(ntt.eq.2000) go to 100
tnb = temp(nrow,nj)
tnt = temp(1,nj)
hf = hfb
tn = tnb
to = tob
i2 = 1
50 if(i2.eq.2) then
hf = hft
tn = tnt
to = tot
end if
if(tn.gt.(to+0.01)) hf=hf+0.1
if(tn.gt.(to+0.005)) hf=hf+0.01
if(tn.gt.(to+0.001)) hf=hf+0.002
if(tn.gt.(to+0.0005)) hf=hf+0.0005
if(tn.gt.(to+0.0001)) hf=hf+0.00002
if(tn.gt.(to+0.00005))hf=hf+0.000005
if(tn.lt.(to-0.01)) hf=hf-0.1
if(tn.lt.(to-0.005)) hf=hf-0.01
if(tn.lt.(to-0.001)) hf=hf-0.002
if(tn.lt.(to-0.0005)) hf=hf-0.0005
if(tn.lt.(to-0.0001)) hf=hf-0.00002
if(tn.lt.(to-0.0001)) hf=hf-0.000005
if(i2.eq.1) then
hfb = hf
tob = tn
end if
if(i2.eq.2) then
hft = hf
tot = tn
end if
if(i2.eq.2) go to 40
i2 = i2 + 1
go to 50
100 stab = 1.0
return
end
c.........SUBROUTINE DIFFUS...................
subroutine diffus (stab,hfb,hft,ntt,nt,tmax,tmin,nrow,nj,
& ncol,factor,conv0,conv1,conv2,conv3,conv4)
common krock(50,50),temp(50,50),tempi(50,50)
&ntempi(50,50),coefx(50,50),coefy(50,50),ntgradi(50),
&ntemp(1000,50,50),|(50),tempnew(50,50),
&ntgrad(1000,50)
character large*28/' Time-step too large!'/
c....Begin diffusional loop..........
c....First determine heat flow to stabilize grid...
105 ntt = ntt + 1
tmaxnew = 0.0
do 140 i=1,nrow
do 130 j=1,ncol
c......Calculate new temperatures from diffusion
equation.
c......with appropriate boundary conditions and
diffusion.
c......coefficients......................................
if(i.eq.1.and.j.eq.1) go to 110
if(i.eq.1.and.j.gt.1.and.j.lt.ncol) go to 112
if(i.eq.1.and.j.eq.ncol) go to 114
if(i.gt.1.and.i.lt.nrow.and.j.eq.ncol) go to 116
if(i.eq.nrow.and.j.eq.ncol) go to 118
if(i.eq.nrow.and.j.gt.1.and.j.lt.ncol) go to 120
if(i.eq.nrow.and.j.eq.1) go to 122
if(i.gt.1.and.i.lt.nrow.and.j.eq.1) go to 124
c....Non-boundary cells with convection.................c
convy = 1.0
convx = 1.0
if(stab.eq.0.0) go to 109
if(krock(i,j).eq.0) then
if(conv0.lt.0.0) convx = -conv0
if(conv0.gt.0.0) convy = conv0
end if
if(krock(i,j).eq.1) then
if(conv1.lt.0.0) convx = -conv1
if(conv1.gt.0.0) convy = conv1
end if
if(krock(i,j).eq.2) then
if(conv2.lt.0.0) convx = -conv2
if(conv2.gt.0.0) convy = conv2
end if
if(krock(i,j).eq.3) then
if(conv3.lt.0.0) convx = -conv3
if(conv3.gt.0.0) convy = conv3
end if
if(krock(i,j).eq.4) then
if(conv4.lt.0.0) convx = -conv4
if(conv4.gt.0.0) convy = conv4
end if
109 cx1 = (convy * coefx(i,j) + coefx(i+1,j))/2.0
cx2 = (coefx(i-1,j) + coefx(i,j))/2.0
cx = (cx1 + cx2)/2.0
cy1 = (convx * coefy(i,j) + coefy(i,j+1))/2.0
cy2 = (coefy(i,j-1) + coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *-
& ( (temp(i-1,j)*cx2) + (temp(i+1,j)*cx1) -
& (2.0*temp(i,j)*cx) +
& (temp(i,j-1)*cy2) + (temp(i,j+1)*cy1) -
& (2.0*temp(i,j)*cy))
go to 125
c....Top left corner cell.......................
110 cx1 = (coefx(i,j) + coefx(i+1,j))/2.0
cx2 = (2.0 * coefx(i,j))/2.0
cx = hft * (cx1 + cx2)/2.0
cy1 = (coefy(i,j) + coefy(i,j+1))/2.0
cy2 = (2.0 * coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ((temp(i,j+1)*cy1) - (temp(i,j)*cy) +
& (temp(i+1,j)*cx1) - (temp(i,j)*cx))
go to 125
c....Top margin of grid........................c
112 cx1 = (coefx(i,j) + coefx(i+1,j))/2.0
cx2 = (2.0 * coefx(i,j))/2.0
cx = hft * (cx1 + cx2)/2.0
cy1 = (coefy(i,j) + coefy(i,j+1))/2.0
cy2 = (coefy(i,j-1) + coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ((temp(i,j-1)*cy2) + (temp(i,j+1)*cy1)-
& (2.0*temp(i,j)*cy) +
& (temp(i+1,j)*cx1) - (temp(i,j)*cx))
go to 125
c....Top right corner cell.......................
114 cx1 = (coefx(i,j) + coefx(i+1,j))/2.0
cx2 = 2.0 * coefx(i,j))/2.0
cx = hft * (cx1 + cx2)/2.0
cy1 = (2.0 * coefy(i,j))/2.0
cy2 = (coefy(i,j-1) + coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor*
& ((temp(i,j-1)*cy2) - (temp(i,j)*cy) +
& (temp(i+1,j)*cx1) - (temp(i,j)*cx))
go to 125
c....Right margin of grid.......................
116 cx1 = (coefx(i,j) + coefx(i+1,j))/2.0
cx2 = (coefx(i-1,j) + coefx(i,j))/2.0
cx = (cx1 + cx2)/2.0
cy1 = (2.0 * coefy(i,j))/2.0
cy2 = (coefy(i,j-1) + coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ( (temp(i,j-1)*cy2) - (temp(i,j)*cy) +
& (temp(i-1,j)*cx2) + (temp(i+1,j)*cx1)-
& (2.0*temp(i,j)*cx))
go to 125
c....Bottom right corner cell, maintain heat flow.........
118 cx1 = (2.0 * coefx(i,j))/2.0
cx2 = (coefx(i-1,j) + coefx(i,j))/2.0
cx = hfb * (cx1 + cx2)/2.0
cy1 = (2.0 * coefy(i,j))/2.0
cy2 = (coefy(i,j-1) + coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ( (temp(i,j-1)*cy2) - (temp(i,j)*cy) +
& (temp(i-1,j)*cx2) - (temp(i,j)*cx))
go to 125
c....Bottom margin of grid, maintain heat flow.........
120 cx1 = (2.0 * coefx(i,j))/2.0
cx2 = (coefx(i-1,j) + coefx(i,j))/2.0
cx = hfb * (cx1 + cx2)/2.0
cy1 = (coefy(i,j) + coefy(i,j+1))/2.0
cy2 = (coefy(i,j-1) + coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ((temp(i,j-1)* cy2) + (temp(i,j+1)*cy1)-
& (2.0*temp(i,j)*cy) +
& (temp(i-1,j)*cx2) - (temp(i,j)*cx))
go to 125
c....Bottom left corner cell, maintain heat
flow..............
122 cx1 = (2.0 * coefx(i,j))/2.0
cx2 = (coefx(i-1,j) + coefx(i,j))/2.0
cx = hfb * (cx1 + cx2)/2.0
cy1 = (coefy(i,j) + coefy(i,j+1))/2.0
cy2 = (2.0 * coefy(i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ((temp(i,j+1)*cy1) - (temp(i,j)*cy) +
& (temp(i-1,j)*cx2) - (temp(i,j)*cx))
go to 125
c....Left margin of grid........................
124 cx1 = (coefx(i,j) + coefx(i+1,j))/2.0
cx2 = (coefx(i-1,j) + coefx(i,j))/2.0
cx = (cx1 + cx2)/2.0
cy1 = (coefy(i,j) + coefy (i,j+1))/2.0
cy2 = (2.0 * coefy (i,j))/2.0
cy = (cy1 + cy2)/2.0
tempnew(i,j) = temp(i,j) + factor *
& ((temp(i,j+1)*cy1) - (temp(i,j)*cy) +
& (temp(i-1,j)*cx2) + (temp(i+1,j)*cx1)-
& (2.0*temp(i,j)*cx))
go to 125
c....Set new max and min magma
temperature......................
125 if((temp(i,j)-tempnew(i,j)).gt.500) then
write(5,'(//,28a)') large
go to 500
end if
if(krock(i,j).eq.0) then
if(tempnew(i,j).ge.tmaxnew)
& tmaxnew = tempnew(i,j)
if(tempnew(i,j).le. tmin) tmin = tempnew(i,j)
end if
130 continue
140 continue
c....Reset temp and test for stability.............
do 144 i=1, nrow
do 142 j=1, ncol
temp(i,j) = tempnew (i,j)
142 continue
144 continue
if(stab.eq.0.0) go to 500
if(tmaxnew.le.600.0) conv0 = 1.0
if(tmaxnew.gt.tmax.and.nt.gt.5) then
write(5,200)
go to 500
end if
200 format(///,'MELT DOWN!!! (Magma is heating up)')
tmax = tmaxnew
if(ntt.lt.1000) go to 105
500 return
end