Appendix 5. Listing of FORTRAN programs used in depth estimation.
c Program bedrk_depth
c-------------------------------------------------------------------------
c---------> Function: calculate grid of depth to bedrock beneath 2 density
c functions from gravity anomaly and depth-to-contact grids
c-------------------------------------------------------------------------
c
c-------------------------------------------------------------------------
c---------> Parameters and notes:
c Model is pile of vertical axis cylinders with coincident axes
c for prism beneath grid point, pile of vertical lines of mass
c for prisms away from grid point. Requires a gravity anomaly
c grid, a depth to bottom of density fn.1 grid, and
c a mask grid for no computation over bedrock or control (drill
c hole, etc) points. Mask val -100. or less means compute; >-100.
c and <0.0 will use absolute value of mask if calculated estimate
c is less than mask value
c Requires a 2 files of {ztop,zbot,drho,radius}, max 101 layers
c describing the density contrast and radius of the pile of cyls
c for the shallow and deep density function
c-------------------------------------------------------------------------
c date programmer event
c-------------------------------------------------------------------------
c 7 Jul 97 MEG created
c 21 Jun 99 MEG minimum depth constraint from mask added
c*-----------------------------------------------------------------------*/
c*Disclaimer: Although this program has been used by the USGS, no */
c*warranty, expressed or implied, is made by the USGS or the United */
c*States Government as to the accuracy and functioning of the program */
c*and related program material nor shall the fact of distribution */
c*constitute any such warranty, and no responsibility is assumed by the */
c*USGS in connection therewith. */
c*-----------------------------------------------------------------------*/
c
character filnam1*50, dflt50*50, qu*80, fmt*50
real rhoeff, gxsl, gxcy, zz(501),dgz(501),raddep
real rho1d(101,4), rho2d(101,4)
real rho1(101,4), rho2(101,4), zbar1(101), zbar2(101)
character dflt1*1,qr*1
character id*56,dflt56*56,pgm*8,id2*56,pgm2*8
character id3*56,pgm3*8
real z1(110), z2(110), z3(110)
real grv(100,100), dpth(100,100), dmask(100,100)
real zc(100,100), res(100,100), dgm(100,100)
real*8 rms
logical verify
verify=.false.
itin=5
itout=6
dval=1.0E+38
c
c Values .ge. 1.0E+38 flag grid points for no compute
c
c
c First get density functions and compute dg(mgals) vs depth to
c contact function to use to estimate new contact depth changes
c
io2 = 11
odpth = 0.0
odz = 0.1
rhoeff = 0.0
gxcy = 0.0
c
c get filename & read file of depths to layers, drho, radius
c
qu='Enter shallow density contrast & radius filename [dshlow.log]-->'
dflt50 = 'dshlow.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
nlayrs1 = 1
write(itout,202)' Shallow density function'
202 format(a)
write(itout,208)
208 format(6x,'ztop',6x,'zbot',6x,'drho',7x,'rad',6x,'zbar',6x,'dgcum')
500 read(io2,fmt,end=510) (rho1(nlayrs1,i), i=1,4)
zbar1(nlayrs1) = (rho1(nlayrs1,1)+rho1(nlayrs1,2))*0.5
rhoeff = rhoeff + rho1(nlayrs1,3)*(rho1(nlayrs1,2)-rho1(nlayrs1,1))
zt = rho1(nlayrs1,1)
zb = rho1(nlayrs1,2)
dr = rho1(nlayrs1,3)
rad = rho1(nlayrs1,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad)
write(itout,207) (rho1(nlayrs1,i), i=1,4), zbar1(nlayrs1), gxcy
207 format(6f10.3)
nlayrs1 = nlayrs1 + 1
go to 500
510 nlayrs1= nlayrs1 - 1
rhoeff = rhoeff / (rho1(nlayrs1,2) - rho1(1,1))
gxsl = 41.908846 * rhoeff * (rho1(nlayrs1,2) - rho1(1,1))
write(itout,204) rhoeff, gxsl, gxcy
204 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
qu='Enter deeper density contrast & radius filename [ddeep.log]-->'
dflt50 = 'ddeep.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
rhoeff1 = 0.0
gxcy = 0.0
nlayrs2 = 1
write(itout,202)' Deep density function'
write(itout,208)
502 read(io2,fmt,end=512) (rho2(nlayrs2,i), i=1,4)
zbar2(nlayrs2) = (rho2(nlayrs2,1)+rho2(nlayrs2,2))*0.5
rhoeff1 = rhoeff1 + rho2(nlayrs2,3)*(rho2(nlayrs2,2)-rho2(nlayrs2,1))
zt = rho2(nlayrs2,1)
zb = rho2(nlayrs2,2)
dr = rho2(nlayrs2,3)
rad = rho2(nlayrs2,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad)
write(itout,207) (rho2(nlayrs2,i), i=1,4), zbar2(nlayrs2), gxcy
nlayrs2 = nlayrs2 + 1
go to 502
512 nlayrs2 = nlayrs2 - 1
rhoeff1 = rhoeff1 / (rho2(nlayrs2,2) - rho2(1,1))
gxsl = 41.908846 * rhoeff1 * (rho2(nlayrs2,2) - rho2(1,1))
write(itout,203) rhoeff1, gxsl, gxcy
203 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
qu = 'Enter radius for depth estimates (km)[10.0]-->'
call ttinr4(qu, raddep, 10.0, itin, itout, verify)
c
c now get grids & estimate depth to contact
c grv is residual (basin) anomaly, mgals
c dpth is depth to density fn 1- i.e. contact, km;
c dmask is mask telling whether
c on bedrock (dval), point of computation(-1.0), or control point
c (value is depth to contact)
c Use USGS BGP standard grid format grids
c
call getfile(21,'Grav grid filename:','in','unformatted')
call getfile(22,'Contact dpth grid filename:','in','unformatted')
call getfile(24,'Mask grid filename:','in','unformatted')
read(21) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
read(22) id2,pgm2,ncol2,nrow2,nz2,xo2,dx2,yo2,dy2
read(24) id3,pgm3,ncol3,nrow3,nz3,xo3,dx3,yo3,dy3
if(nz.ne.1.or.nz2.ne.1.or.nz3.ne.1) then
print *, ' nz not = 1 for input grid...'
stop
endif
call compspecs(ncol,nrow,xo,dx,yo,dy,
& ncol2,nrow2,xo2,dx2,yo2,dy2)
call compspecs(ncol,nrow,xo,dx,yo,dy,
& ncol3,nrow3,xo3,dx3,yo3,dy3)
c
c store sqrt(dx * dy)/sqrt(pi) as radius in density fn
c (radius of circle same area as dx*dy
c
eqrad = sqrt(dx*dy/3.141527)
do 40 i=1, nlayrs1
rho1d(i,1) = rho1(i,1)
rho1d(i,2) = rho1(i,2)
rho1d(i,3) = rho1(i,3)
rho1d(i,4) = raddep
40 rho1(i,4) = eqrad
do 42 i=1, nlayrs2
rho2d(i,1) = rho2(i,1)
rho2d(i,2) = rho2(i,2)
rho2d(i,3) = rho2(i,3)
rho2d(i,4) = raddep
42 rho2(i,4) = eqrad
do 60 j=1,nrow
call getrow(21,ncol,z1)
call getrow(22,ncol,z2)
call getrow(24,ncol,z3)
do 62 i=1,ncol
grv(j,i) = z1(i)
dpth(j,i) = z2(i)
dmask(j,i) = z3(i)
c
c load dval if any input grid is dval or gravity anomaly > 0.0
c
if(grv(j,i).ge.0.0.or.dpth(j,i).ge.dval.or.dmask(j,i).ge.dval)
& zc(j,i) = dval
c
c estimate depth to interface; if control point, enter mask value
c
if(dmask(j,i).ge.0.0) zc(j,i) = dmask(j,i)
if(dmask(j,i).lt.0.0) then
zmin = dpth(j,i)
320 zmax = 5.0
nzval = 101
dz = (zmax - zmin)/(nzval-1)
c
c calculate gravity as fn of interface depth
c
do 50 k=1,nzval
zct = zmin + (k-1)*dz
call dgtot(0.0, zmin, zct, nlayrs1, rho1d, nlayrs2, rho2d, dg)
zz(nzval-k+1) = zct
dgz(nzval-k+1) = dg
50 continue
c
c interpolate depth to bottom from g vs z fn
c
call piksr2(nzval,dgz,zz)
grvji = grv(j,i)
call linterp(dgz,zz,nzval,grvji,zcji)
if(zcji.lt.zmin) zcji = zmin
if(zcji.gt.zmax) zcji = zmax
if(zcji.lt.0.0) zcji = 0.0
c
c Check if this is a deeper than constraint from a well
c
if(dmask(j,i).gt.-100.and.zcji.lt.abs(dmask(j,i)))
& zcji = abs(dmask(j,i))
zc(j,i) = zcji
endif
62 continue
60 continue
c
c grids loaded; calculate model gravity and residual
c
300 rms = 0.0
do 70 j=1,nrow
do 71 i=1,ncol
if(grv(j,i).ge.0.0.or.dpth(j,i).ge.dval.or.dmask(j,i).ge.dval) then
dgm(j,i) = dval
go to 71
endif
dgm(j,i) = 0.0
do 72 jj=1,nrow
do 73 ii=1,ncol
if(zc(jj,ii).ge.dval) go to 73
yy = (jj-j)*dy
xx = (ii-i)*dx
radius = sqrt(xx*xx+yy*yy)
if(radius.gt.10.0) go to 73
call dgtot(radius, zc(jj,ii), dpth(jj,ii), nlayrs1, rho1,
& nlayrs2, rho2, dg)
dgm(j,i) = dgm(j,i) + dg
73 continue
72 continue
resid = grv(j,i) - dgm(j,i)
rms = rms + resid*resid
res(j,i) = resid
71 continue
70 continue
rmserr = sqrt(rms/(ncol*nrow))
write(itout,'(a,f10.2)')' rms error = ',rmserr
qu = 'another iteration (a) or stop (q)[a]'
dflt1 = 'a'
call ttinaa(qu,qr,1,dflt1,itin,itout,verify)
if(qr.eq.'q'.or.qr.eq.'Q') go to 310
c
c iterate by adjusting contact depths
c
do 80 j=1,nrow
do 80 i=1,ncol
if(zc(j,i).ge.dval) go to 80
if(dmask(j,i).ge.0.0) zc(j,i) = dmask(j,i)
if(dmask(j,i).lt.0.0) then
zmin = dpth(j,i)
330 zmax = 5.0
nzval = 51
dz = (zmax - zmin)/(nzval-1)
c
c calculate gravity as fn of interface depth
c
call dgtot(0.0, zmin, zc(j,i), nlayrs1, rho1d, nlayrs2, rho2d, dgcur)
do 90 k=1,nzval
zct = zmin + (k-1)*dz
call dgtot(0.0, zmin, zct, nlayrs1, rho1d, nlayrs2, rho2d, dg)
zz(nzval-k+1) = zct
dgz(nzval-k+1) = dg - dgcur
90 continue
c
c interpolate interface depth from g vs z fn
c
call piksr2(nzval,dgz,zz)
resji = res(j,i)
call linterp(dgz,zz,nzval,resji,zcji)
if(zcji.lt.zmin) zcji = zmin
if(zcji.gt.zmax) zcji = zmax
if(zcji.lt.0.0) zcji = 0.0
c
c Check if this is a deeper than constraint from a well
c
if(dmask(j,i).gt.-100.and.zcji.lt.abs(dmask(j,i)))
& zcji = abs(dmask(j,i))
zc(j,i) = zcji
endif
80 continue
c
c recalc model, residual
c
go to 300
c
c output grids
c
310 call getfile(23,'Bedrock depth grid filename:','out','unformatted')
qu='Enter id for output grid (<57c)[Bedrk_depth bdrk depth grid]'
dflt56='Bedrk_depth bdrk depth grid'
call ttinaa(qu,id,56,dflt56,itin,itout,verify)
pgm='bedrk_de'
write(23) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
do 92 j=1,nrow
do 93 i=1,ncol
93 z1(i) = zc(j,i)
call putrow(23,ncol,z1)
92 continue
close(23)
call getfile(23,'Output resid grid filename:','out','unformatted')
qu='Enter id for output grid (<57c)[Bedrk_depth output residual grid]'
dflt56='grvdpth output residual grid'
call ttinaa(qu,id,56,dflt56,itin,itout,verify)
pgm='bedrk_de'
write(23) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
do 94 j=1,nrow
do 95 i=1,ncol
95 z1(i) = res(j,i)
call putrow(23,ncol,z1)
94 continue
close(23)
call getfile(23,'Gravity model grid filename:','out','unformatted')
qu='Enter id for output grid (<57c)[Bedrk_depth output grav model grid]'
dflt56='grvdpth output grav model grid'
call ttinaa(qu,id,56,dflt56,itin,itout,verify)
pgm='bedrk_de'
write(23) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
do 96 j=1,nrow
do 97 i=1,ncol
97 z1(i) = dgm(j,i)
call putrow(23,ncol,z1)
96 continue
900 close(21)
close(22)
close(23)
close(24)
stop
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getfile(unit,question,inout,form)
c Opens files...filenames not recognized by VAX will generate
c error message and request a second try.
c Sample call:
c call getfile(21,'*Give input file:','in','formatted')
integer unit
character filename*50, status*8
character*(*) question, inout, form
print *, ' '
10 write(*,101) question
101 format(a,' ',$)
read '(a50)', filename
if(inout.eq.'in') then
status='old'
else
status='unknown'
endif
if(form.eq.'formatted'.and.inout.eq.'out') then
open(unit,file=filename,form=form,status=status,
& err=900)
else
open(unit,file=filename,form=form,status=status,err=900)
endif
return
900 print *,'ERROR IN OPENING ',filename
print *,' try again......'
close(unit)
go to 10
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine compspecs(ncol,nrow,xo,dx,yo,dy,
& ncol2,nrow2,xo2,dx2,yo2,dy2)
c Compares the parameters for two grids...
c If parameters are different, stops execution.
if(ncol.ne.ncol2.or.nrow.ne.nrow2
& .or.xo.ne.xo2.or.dx.ne.dx2
& .or.yo.ne.yo2.or.dy.ne.dy2) then
print *,'***ERROR...GRIDSPECS DO NOT MATCH'
stop
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getrow(iunit,ncol,zrow)
real zrow(ncol)
read(iunit) dum,zrow
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine putrow(iunit,ncol,zrow)
real zrow(ncol)
dum=0.0
write(iunit) dum,zrow
return
end
function stkcrcyl(h1,h2,nlayrs,rho)
c
c calc dg of stack of coincident axis vert cylinders from depth h1
c down to depth h2 using array rho for ztop,zbot,drho,rad
c NOTE rho(nlayrs,2) must be greater than H
c
real h1, h2, rho(101,4)
real*8 dg
integer nlayrs
dg = 0.0d0
do 50 i=1,nlayrs
if(h1.ge.rho(i,2)) go to 50
zt = rho(i,1)
zb = rho(i,2)
dr = rho(i,3)
rad = rho(i,4)
if(h1.ge.rho(i,1).and.h1.lt.rho(i,2)) then
if(h2.lt.rho(i,2)) go to 57
dg = dg + cylfn(h1,zb,dr,rad)
go to 50
endif
if(h2.lt.rho(i,2)) go to 55
dg = dg + cylfn(zt,zb,dr,rad)
50 continue
go to 60
55 dg = dg + cylfn(zt,h2,dr,rad)
go to 60
57 dg = dg + cylfn(h1,h2,dr,rad)
60 stkcrcyl = dg
return
end
function cylfn(ztop,zbot,drho,rad)
c gravity effect of vert. axis circular cylinder at pt on axis
real*8 t
t = zbot - ztop + sqrt(rad*rad + ztop*ztop)
1 - sqrt(rad*rad + zbot*zbot)
cylfn = 41.908846 * drho * t
return
end
subroutine dgtot(r, zct, ztot, nlayrs1, rho1, nlayrs2, rho2, dg)
c
c subroutine to calculate gravity effect of a stack of prisms
c with density fn rho1 to depth zct, rho2 from zct to zbot
c uses right circular cylinder if r=0, line of mass formula
c if r > 0
c MEG/USGS/Dec96
c
dimension rho1(101,4),rho2(101,4)
c
c decide whether to use cylinder or line mass formula
c
if(r.lt.0.0001) then
c
c r (horizontal distance) is 0 - use circular cyl on axis formula
c
dg = stkcrcyl(0.0,zct,nlayrs1,rho1)
dg = dg + stkcrcyl(zct,ztot,nlayrs2,rho2)
return
else
c
c r (horizonatal distance > 0 - use line of mass formula
c
dg = stklinms(r,0.0,zct,nlayrs1,rho1)
dg = dg + stklinms(r,zct,ztot,nlayrs2,rho2)
return
endif
end
function stklinms(r,h1,h2,nlayrs,rho)
c
c calc dg of stack of coincident axis vert mass lines from depth h1
c down to depth h2 using array rho for ztop,zbot,drho,rad
c NOTE rho(nlayrs,2) must be greater than H; rad is area prism
c so rad*drho is mass/unit length
c
real h1, h2, rho(101,4), r, lmsfn, area
real*8 dg
integer nlayrs
dg = 0.0d0
do 50 i=1,nlayrs
if(h1.ge.rho(i,2)) go to 50
zt = rho(i,1)
zb = rho(i,2)
dr = rho(i,3)
rad = rho(i,4)*1.772454
area = rad * rad
c note sqrt pi factor for program grvdpth only where rad stored in rho(i,4)
c is sqrt((dx*dy)/pi) remove factor for general use
if(h1.ge.rho(i,1).and.h1.lt.rho(i,2)) then
if(h2.lt.rho(i,2)) go to 57
dg = dg + lmsfn(r,h1,zb,dr,area)
go to 50
endif
if(h2.lt.rho(i,2)) go to 55
dg = dg + lmsfn(r,zt,zb,dr,area)
50 continue
go to 60
55 dg = dg + lmsfn(r,zt,h2,dr,area)
go to 60
57 dg = dg + lmsfn(r,h1,h2,dr,area)
60 stklinms = dg
return
end
real function lmsfn(r,zt,zb,dr,area)
c gravity effect of vertical line of mass from zt to zb,
c horiz. distance r, area*dr mass/unit z,units km,g/cc,mgal
real*8 t1, t2
real area, dr, r, zt, zb
t1 = 1.0/sqrt(r*r+zt*zt)
t2 = 1.0/sqrt(r*r+zb*zb)
lmsfn = area*dr*6.67*(t1-t2)
return
end
SUBROUTINE PIKSR2(N,ARR,BRR)
DIMENSION ARR(N),BRR(N)
DO 12 J=2,N
A=ARR(J)
B=BRR(J)
DO 11 I=J-1,1,-1
IF(ARR(I).LE.A)GO TO 10
ARR(I+1)=ARR(I)
BRR(I+1)=BRR(I)
11 CONTINUE
I=0
10 ARR(I+1)=A
BRR(I+1)=B
12 CONTINUE
RETURN
END
subroutine linterp(x,y,n,xx,yy)
c
c linear interpolation on x, gives endpoint value if x outside
c interval (x(1),x(n)). assumes x monotonic increasing.
c assumes number smaller in abs value 1.0e-15 is zero. meg/jul97
c
dimension x(n),y(n)
tol = 1.0e-15
tolinv = 1.0e15
if(xx.lt.x(1)) then
yy = y(1)
return
endif
if(xx.gt.x(n)) then
yy = y(n)
return
endif
do 50 i=1,n-1
if(xx.ge.x(i).and.xx.lt.x(i+1)) then
dy = y(i+1)-y(i)
dx = x(i+1)-x(i)
if(abs(dx).le.tol) then
slop = dy*tolinv
else
slop = dy/dx
endif
yy = slop*(xx-x(i))+y(i)
return
endif
50 continue
end
c Program contact_depth
c-------------------------------------------------------------------------
c---------> Function: calculate grid of depth to contact between 2 density
c functions from gravity anomaly and total depth grids
c-------------------------------------------------------------------------
c
c-------------------------------------------------------------------------
c---------> Parameters and notes:
c Model is pile of vertical axis cylinders with coincident axes
c for prism beneath grid point, pile of vertical lines of mass
c for prisms away from grid point. Requires a gravity anomaly
c grid, a depth to bottom of density fn.2 grid, and
c a mask grid for no computation over bedrock or control (drill
c hole, etc) points. Mask val -100. or less means compute; >-100.
c and <0.0 will use absolute value of mask if calculated estimate
c is less than mask value
c Requires a 2 files of {ztop,zbot,drho,radius}, max 101 layers
c describing the density contrast and radius of the pile of cyls
c for the shallow and deep density function
c-------------------------------------------------------------------------
c date programmer event
c-------------------------------------------------------------------------
c 6 Dec 96 MEG created
c 21 Jun 99 MEG added mininimum depth constraint for mask
c*-----------------------------------------------------------------------*/
c*Disclaimer: Although this program has been used by the USGS, no */
c*warranty, expressed or implied, is made by the USGS or the United */
c*States Government as to the accuracy and functioning of the program */
c*and related program material nor shall the fact of distribution */
c*constitute any such warranty, and no responsibility is assumed by the */
c*USGS in connection therewith. */
c*-----------------------------------------------------------------------*/
c
character filnam1*50, dflt50*50, qu*80, fmt*50
real rhoeff, gxsl, gxcy, zz(501),dgz(501),raddep
real rho1d(101,4), rho2d(101,4)
real rho1(101,4), rho2(101,4), zbar1(101), zbar2(101)
character dflt1*1,qr*1
character id*56,dflt56*56,pgm*8,id2*56,pgm2*8
character id3*56,pgm3*8
real z1(110), z2(110), z3(110)
real grv(100,100), dpth(100,100), dmask(100,100)
real zc(100,100), res(100,100), dgm(100,100)
real*8 rms
logical verify
verify=.false.
itin=5
itout=6
dval=1.0E+38
c
c Values .ge. 1.0E+38 flag grid points for no compute
c
c
c First get density functions and compute dg(mgals) vs depth to
c contact function to use to estimate new contact depth changes
c
io2 = 11
odpth = 0.0
odz = 0.1
rhoeff = 0.0
gxcy = 0.0
c
c get filename & read file of depths to layers, drho, radius
c
qu='Enter shallow density contrast & radius filename [dshlow.log]-->'
dflt50 = 'dshlow.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
nlayrs1 = 1
write(itout,202)' Shallow density function'
202 format(a)
write(itout,208)
208 format(6x,'ztop',6x,'zbot',6x,'drho',7x,'rad',6x,'zbar',6x,'dgcum')
500 read(io2,fmt,end=510) (rho1(nlayrs1,i), i=1,4)
zbar1(nlayrs1) = (rho1(nlayrs1,1)+rho1(nlayrs1,2))*0.5
rhoeff = rhoeff + rho1(nlayrs1,3)*(rho1(nlayrs1,2)-rho1(nlayrs1,1))
zt = rho1(nlayrs1,1)
zb = rho1(nlayrs1,2)
dr = rho1(nlayrs1,3)
rad = rho1(nlayrs1,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad)
write(itout,207) (rho1(nlayrs1,i), i=1,4), zbar1(nlayrs1), gxcy
207 format(6f10.3)
nlayrs1 = nlayrs1 + 1
go to 500
510 nlayrs1= nlayrs1 - 1
rhoeff = rhoeff / (rho1(nlayrs1,2) - rho1(1,1))
gxsl = 41.908846 * rhoeff * (rho1(nlayrs1,2) - rho1(1,1))
write(itout,204) rhoeff, gxsl, gxcy
204 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
qu='Enter deeper density contrast & radius filename [ddeep.log]-->'
dflt50 = 'ddeep.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
rhoeff1 = 0.0
gxcy = 0.0
nlayrs2 = 1
write(itout,202)' Deep density function'
write(itout,208)
502 read(io2,fmt,end=512) (rho2(nlayrs2,i), i=1,4)
zbar2(nlayrs2) = (rho2(nlayrs2,1)+rho2(nlayrs2,2))*0.5
rhoeff1 = rhoeff1 + rho2(nlayrs2,3)*(rho2(nlayrs2,2)-rho2(nlayrs2,1))
zt = rho2(nlayrs2,1)
zb = rho2(nlayrs2,2)
dr = rho2(nlayrs2,3)
rad = rho2(nlayrs2,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad)
write(itout,207) (rho2(nlayrs2,i), i=1,4), zbar2(nlayrs2), gxcy
nlayrs2 = nlayrs2 + 1
go to 502
512 nlayrs2 = nlayrs2 - 1
rhoeff1 = rhoeff1 / (rho2(nlayrs2,2) - rho2(1,1))
gxsl = 41.908846 * rhoeff1 * (rho2(nlayrs2,2) - rho2(1,1))
write(itout,203) rhoeff1, gxsl, gxcy
203 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
qu = 'Enter radius for depth estimates (km)[10.0]-->'
call ttinr4(qu, raddep, 10.0, itin, itout, verify)
c
c now get grids & estimate depth to contact
c grv is residual (basin) anomaly, mgals
c dpth is depth to bottom of basin, km; dmask is mask telling whether
c on bedrock (dval), point of computation(-1.0), or control point
c (value is depth to contact)
c Use USGS BGP standard grid format grids
c
call getfile(21,'Grav grid filename:','in','unformatted')
call getfile(22,'Depth grid filename:','in','unformatted')
call getfile(24,'Mask grid filename:','in','unformatted')
read(21) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
read(22) id2,pgm2,ncol2,nrow2,nz2,xo2,dx2,yo2,dy2
read(24) id3,pgm3,ncol3,nrow3,nz3,xo3,dx3,yo3,dy3
if(nz.ne.1.or.nz2.ne.1.or.nz3.ne.1) then
print *, ' nz not = 1 for input grid...'
stop
endif
call compspecs(ncol,nrow,xo,dx,yo,dy,
& ncol2,nrow2,xo2,dx2,yo2,dy2)
call compspecs(ncol,nrow,xo,dx,yo,dy,
& ncol3,nrow3,xo3,dx3,yo3,dy3)
c
c store sqrt(dx * dy)/sqrt(pi) as radius in density fn
c (radius of circle same area as dx*dy
c
eqrad = sqrt(dx*dy/3.141527)
do 40 i=1, nlayrs1
rho1d(i,1) = rho1(i,1)
rho1d(i,2) = rho1(i,2)
rho1d(i,3) = rho1(i,3)
rho1d(i,4) = raddep
40 rho1(i,4) = eqrad
do 42 i=1, nlayrs2
rho2d(i,1) = rho2(i,1)
rho2d(i,2) = rho2(i,2)
rho2d(i,3) = rho2(i,3)
rho2d(i,4) = raddep
42 rho2(i,4) = eqrad
do 60 j=1,nrow
call getrow(21,ncol,z1)
call getrow(22,ncol,z2)
call getrow(24,ncol,z3)
do 62 i=1,ncol
grv(j,i) = z1(i)
dpth(j,i) = z2(i)
dmask(j,i) = z3(i)
c
c load dval if any input grid is dval or gravity anomaly > 0.0
c
if(grv(j,i).ge.0.0.or.dpth(j,i).ge.dval.or.dmask(j,i).ge.dval)
& zc(j,i) = dval
c
c estimate depth to interface; if control point, enter mask value
c
if(dmask(j,i).ge.0.0) zc(j,i) = dmask(j,i)
if(dmask(j,i).lt.0.0) then
zmin = 0.0
zmax = dpth(j,i)
nzval = 101
dz = (zmax - zmin)/(nzval-1)
c
c calculate gravity as fn of interface depth
c
do 50 k=1,nzval
zct = zmin + (k-1)*dz
call dgtot(0.0, zct, zmax, nlayrs1, rho1d, nlayrs2, rho2d, dg)
zz(nzval-k+1) = zct
dgz(nzval-k+1) = dg
50 continue
c
c interpolate interface depth from g vs z fn
c
call piksr2(nzval,dgz,zz)
grvji = grv(j,i)
call linterp(dgz,zz,nzval,grvji,zcji)
if(zcji.lt.zmin) zcji = zmin
if(zcji.gt.zmax) zcji = zmax
if(zcji.lt.0.0) zcji = 0.0
c
c Check if this is a deeper than constraint from a well
c
if(dmask(j,i).gt.-100.and.zcji.lt.abs(dmask(j,i)))
& zcji = abs(dmask(j,i))
zc(j,i) = zcji
endif
62 continue
60 continue
c
c grids loaded; calculate model gravity and residual
c
300 rms = 0.0
do 70 j=1,nrow
do 71 i=1,ncol
if(grv(j,i).ge.0.0.or.dpth(j,i).ge.dval.or.dmask(j,i).ge.dval) then
dgm(j,i) = dval
go to 71
endif
dgm(j,i) = 0.0
do 72 jj=1,nrow
do 73 ii=1,ncol
if(zc(jj,ii).ge.dval) go to 73
yy = (jj-j)*dy
xx = (ii-i)*dx
radius = sqrt(xx*xx+yy*yy)
if(radius.gt.10.0) go to 73
call dgtot(radius, zc(jj,ii), dpth(jj,ii), nlayrs1, rho1,
& nlayrs2, rho2, dg)
dgm(j,i) = dgm(j,i) + dg
73 continue
72 continue
resid = grv(j,i) - dgm(j,i)
rms = rms + resid*resid
res(j,i) = resid
71 continue
70 continue
rmserr = sqrt(rms/(ncol*nrow))
write(itout,'(a,f10.2)')' rms error = ',rmserr
qu = 'another iteration (a) or stop (q)[a]'
dflt1 = 'a'
call ttinaa(qu,qr,1,dflt1,itin,itout,verify)
if(qr.eq.'q'.or.qr.eq.'Q') go to 310
c
c iterate by adjusting contact depths
c
do 80 j=1,nrow
do 80 i=1,ncol
if(zc(j,i).ge.dval) go to 80
if(dmask(j,i).ge.0.0) zc(j,i) = dmask(j,i)
if(dmask(j,i).lt.0.0) then
zmin = 0.0
zmax = dpth(j,i)
nzval = 51
dz = (zmax - zmin)/(nzval-1)
c
c calculate gravity as fn of interface depth
c
call dgtot(0.0, zc(j,i), zmax, nlayrs1, rho1d, nlayrs2, rho2d, dgcur)
do 90 k=1,nzval
zct = zmin + (k-1)*dz
call dgtot(0.0, zct, zmax, nlayrs1, rho1d, nlayrs2, rho2d, dg)
zz(nzval-k+1) = zct
dgz(nzval-k+1) = dg - dgcur
90 continue
c
c interpolate interface depth from g vs z fn
c
call piksr2(nzval,dgz,zz)
resji = res(j,i)
call linterp(dgz,zz,nzval,resji,zcji)
if(zcji.lt.zmin) zcji = zmin
if(zcji.gt.zmax) zcji = zmax
if(zcji.lt.0.0) zcji = 0.0
c
c Check if this is a deeper than constraint from a well
c
if(dmask(j,i).gt.-100.and.zcji.lt.abs(dmask(j,i)))
& zcji = abs(dmask(j,i))
zc(j,i) = zcji
endif
80 continue
c
c recalc model, residual
c
go to 300
c
c output grids
c
310 call getfile(23,'Contact depth grid filename:','out','unformatted')
qu='Enter id for output grid (<57c)[grvdpth contact depth grid]'
dflt56='contact_depth contact depth grid'
call ttinaa(qu,id,56,dflt56,itin,itout,verify)
pgm='contact_'
write(23) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
do 92 j=1,nrow
do 93 i=1,ncol
93 z1(i) = zc(j,i)
call putrow(23,ncol,z1)
92 continue
close(23)
call getfile(23,'Output resid grid filename:','out','unformatted')
qu='Enter id for output grid (<57c)[grvdpth output residual grid]'
dflt56='contact_depth output residual grid'
call ttinaa(qu,id,56,dflt56,itin,itout,verify)
pgm='contact_'
write(23) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
do 94 j=1,nrow
do 95 i=1,ncol
95 z1(i) = res(j,i)
call putrow(23,ncol,z1)
94 continue
close(23)
call getfile(23,'Gravity model grid filename:','out','unformatted')
qu='Enter id for output grid (<57c)[grvdpth output grav model grid]'
dflt56='contact_depth output grav model grid'
call ttinaa(qu,id,56,dflt56,itin,itout,verify)
pgm='contact_'
write(23) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
do 96 j=1,nrow
do 97 i=1,ncol
97 z1(i) = dgm(j,i)
call putrow(23,ncol,z1)
96 continue
900 close(21)
close(22)
close(23)
close(24)
stop
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getfile(unit,question,inout,form)
c Opens files...filenames not recognized by VAX will generate
c error message and request a second try.
c Sample call:
c call getfile(21,'*Give input file:','in','formatted')
integer unit
character filename*50, status*8
character*(*) question, inout, form
print *, ' '
10 write(*,101) question
101 format(a,' ',$)
read '(a50)', filename
if(inout.eq.'in') then
status='old'
else
status='unknown'
endif
if(form.eq.'formatted'.and.inout.eq.'out') then
open(unit,file=filename,form=form,status=status,
& err=900)
else
open(unit,file=filename,form=form,status=status,err=900)
endif
return
900 print *,'ERROR IN OPENING ',filename
print *,' try again......'
close(unit)
go to 10
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine compspecs(ncol,nrow,xo,dx,yo,dy,
& ncol2,nrow2,xo2,dx2,yo2,dy2)
c Compares the parameters for two grids...
c If parameters are different, stops execution.
if(ncol.ne.ncol2.or.nrow.ne.nrow2
& .or.xo.ne.xo2.or.dx.ne.dx2
& .or.yo.ne.yo2.or.dy.ne.dy2) then
print *,'***ERROR...GRIDSPECS DO NOT MATCH'
stop
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getrow(iunit,ncol,zrow)
real zrow(ncol)
read(iunit) dum,zrow
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine putrow(iunit,ncol,zrow)
real zrow(ncol)
dum=0.0
write(iunit) dum,zrow
return
end
function stkcrcyl(h1,h2,nlayrs,rho)
c
c calc dg of stack of coincident axis vert cylinders from depth h1
c down to depth h2 using array rho for ztop,zbot,drho,rad
c NOTE rho(nlayrs,2) must be greater than H
c
real h1, h2, rho(101,4)
real*8 dg
integer nlayrs
dg = 0.0d0
do 50 i=1,nlayrs
if(h1.ge.rho(i,2)) go to 50
zt = rho(i,1)
zb = rho(i,2)
dr = rho(i,3)
rad = rho(i,4)
if(h1.ge.rho(i,1).and.h1.lt.rho(i,2)) then
if(h2.lt.rho(i,2)) go to 57
dg = dg + cylfn(h1,zb,dr,rad)
go to 50
endif
if(h2.lt.rho(i,2)) go to 55
dg = dg + cylfn(zt,zb,dr,rad)
50 continue
go to 60
55 dg = dg + cylfn(zt,h2,dr,rad)
go to 60
57 dg = dg + cylfn(h1,h2,dr,rad)
60 stkcrcyl = dg
return
end
function cylfn(ztop,zbot,drho,rad)
c gravity effect of vert. axis circular cylinder at pt on axis
real*8 t
t = zbot - ztop + sqrt(rad*rad + ztop*ztop)
1 - sqrt(rad*rad + zbot*zbot)
cylfn = 41.908846 * drho * t
return
end
subroutine dgtot(r, zct, ztot, nlayrs1, rho1, nlayrs2, rho2, dg)
c
c subroutine to calculate gravity effect of a stack of prisms
c with density fn rho1 to depth zct, rho2 from zct to zbot
c uses right circular cylinder if r=0, line of mass formula
c if r > 0
c MEG/USGS/Dec96
c
dimension rho1(101,4),rho2(101,4)
c
c decide whether to use cylinder or line mass formula
c
if(r.lt.0.0001) then
c
c r (horizontal distance) is 0 - use circular cyl on axis formula
c
dg = stkcrcyl(0.0,zct,nlayrs1,rho1)
dg = dg + stkcrcyl(zct,ztot,nlayrs2,rho2)
return
else
c
c r (horizonatal distance > 0 - use line of mass formula
c
dg = stklinms(r,0.0,zct,nlayrs1,rho1)
dg = dg + stklinms(r,zct,ztot,nlayrs2,rho2)
return
endif
end
function stklinms(r,h1,h2,nlayrs,rho)
c
c calc dg of stack of coincident axis vert mass lines from depth h1
c down to depth h2 using array rho for ztop,zbot,drho,rad
c NOTE rho(nlayrs,2) must be greater than H; rad is area prism
c so rad*drho is mass/unit length
c
real h1, h2, rho(101,4), r, lmsfn, area
real*8 dg
integer nlayrs
dg = 0.0d0
do 50 i=1,nlayrs
if(h1.ge.rho(i,2)) go to 50
zt = rho(i,1)
zb = rho(i,2)
dr = rho(i,3)
rad = rho(i,4)*1.772454
area = rad * rad
c note sqrt pi factor for program grvdpth only where rad stored in rho(i,4)
c is sqrt((dx*dy)/pi) remove factor for general use
if(h1.ge.rho(i,1).and.h1.lt.rho(i,2)) then
if(h2.lt.rho(i,2)) go to 57
dg = dg + lmsfn(r,h1,zb,dr,area)
go to 50
endif
if(h2.lt.rho(i,2)) go to 55
dg = dg + lmsfn(r,zt,zb,dr,area)
50 continue
go to 60
55 dg = dg + lmsfn(r,zt,h2,dr,area)
go to 60
57 dg = dg + lmsfn(r,h1,h2,dr,area)
60 stklinms = dg
return
end
real function lmsfn(r,zt,zb,dr,area)
c gravity effect of vertical line of mass from zt to zb,
c horiz. distance r, area*dr mass/unit z,units km,g/cc,mgal
real*8 t1, t2
real area, dr, r, zt, zb
t1 = 1.0/sqrt(r*r+zt*zt)
t2 = 1.0/sqrt(r*r+zb*zb)
lmsfn = area*dr*6.67*(t1-t2)
return
end
SUBROUTINE PIKSR2(N,ARR,BRR)
DIMENSION ARR(N),BRR(N)
DO 12 J=2,N
A=ARR(J)
B=BRR(J)
DO 11 I=J-1,1,-1
IF(ARR(I).LE.A)GO TO 10
ARR(I+1)=ARR(I)
BRR(I+1)=BRR(I)
11 CONTINUE
I=0
10 ARR(I+1)=A
BRR(I+1)=B
12 CONTINUE
RETURN
END
subroutine linterp(x,y,n,xx,yy)
c
c linear interpolation on x, gives endpoint value if x outside
c interval (x(1),x(n)). assumes x monotonic increasing.
c assumes number smaller in abs value 1.0e-15 is zero. meg/jul97
c
dimension x(n),y(n)
tol = 1.0e-15
tolinv = 1.0e15
if(xx.lt.x(1)) then
yy = y(1)
return
endif
if(xx.gt.x(n)) then
yy = y(n)
return
endif
do 50 i=1,n-1
if(xx.ge.x(i).and.xx.lt.x(i+1)) then
dy = y(i+1)-y(i)
dx = x(i+1)-x(i)
if(abs(dx).le.tol) then
slop = dy*tolinv
else
slop = dy/dx
endif
yy = slop*(xx-x(i))+y(i)
return
endif
50 continue
end
subroutine ttinr4(qu, r4, defalt, itin,itout, verify)
c Terminal entry of real*4 variable - F77 version
c qu - Character string prompt que,80 C long.
c r4 - variable to be input
c defalt - default value if <cr> only given
c itin,itout - read, write LUN's of terminal
c verify - if true, asks for verification of input
c if false no verification
c typical call:
c call ttinr4(que,a,1.67,itin,itout,verify)
c MEG/ Feb 86; SUN conversion Jul 89/ MEG
c
logical verify
character ans*10, blank*10, qu*80
character que(80)*1, qr*1
real r4, defalt
blank = ' '
read(unit=qu, fmt=200) que
200 format(80a1)
call lnnobl(que, nq, 80)
300 write(unit=itout, fmt=204)
204 format(1h ,$)
do 40 i = 1, nq
40 write(unit=itout, fmt=201) que(i)
201 format(a,$)
read(unit=itin, fmt=101, err=300) ans
101 format(a)
if (ans .eq. blank) goto 910
read(unit=ans, fmt=100) r4
100 format(bn,f10.0)
if (.not. verify) goto 900
write(unit=itout, fmt=202) r4
202 format(1x,6hValue ,1pe14.7,5h ok? ,$)
read(unit=itin, fmt=101) qr
if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300
900 return
910 r4 = defalt
goto 900
end
subroutine lnnobl(a, na, lngth)
c determine length of string of max length lngth with no blanks on left
c changed to specify string max length from 80 by MEG/10Jan91
character a(lngth), bl
bl = ' '
do 50 i = lngth, 1, -1
na = i
if (a(i) .ne. bl) goto 60
50 continue
na = 0
60 return
end
subroutine ttinr8(qu, r8, defalt, itin,itout, verify)
c Terminal entry of real*8 variable - F77 version
c qu - Character string prompt que,80 C long.
c r8 - variable to be input
c def - default value if <cr> only given
c itin,itout - read, write LUN's of terminal
c verify - if true, asks for verification of input
c if false no verification
c typical call:
c call ttinr8(que,a,1.67d-03,itin,itout,verify)
c MEG/ Feb 86; SUN conversion Jul 89/ MEG
c
logical verify
character ans*16, blank*16, qu*80
character que(80)*1, qr*1
double precision r8, defalt
blank = ' '
read(unit=qu, fmt=200) que
200 format(80a1)
call lnnobl(que, nq, 80)
300 write(unit=itout, fmt=204)
204 format(1h ,$)
do 40 i = 1, nq
40 write(unit=itout, fmt=201) que(i)
201 format(a,$)
read(unit=itin, fmt=101, err=300) ans
101 format(a)
if (ans .eq. blank) goto 910
read(unit=ans, fmt=100) r8
100 format(bn,f16.0)
if (.not. verify) goto 900
write(unit=itout, fmt=202) r8
202 format(1x,6hValue ,1pd24.17,5h ok? ,$)
read(unit=itin, fmt=101) qr
if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300
900 return
910 r8 = defalt
goto 900
end
subroutine ttini4(qu, i4, defalt, itin,itout, verify)
c Terminal entry of integer*4 variable - F77 version
c qu - Character string prompt que,80 C long.
c i4 - variable to be input
c defalt - default value if <cr> only given
c itin,itout - read, write LUN's of terminal
c verify - if true, asks for verification of input
c if false no verification
c typical call:
c call ttini4(que,k,16,itin,itout,verify)
c MEG/ Feb 86; SUN conversion Jul 89/ MEG
c
logical verify
character ans*16, blank*16, qu*80
character que(80)*1, qr*1
integer*4 i4, defalt
blank = ' '
read(unit=qu, fmt=200) que
200 format(80a1)
call lnnobl(que, nq, 80)
300 write(unit=itout, fmt=204)
204 format(1h ,$)
do 40 i = 1, nq
40 write(unit=itout, fmt=201) que(i)
201 format(a,$)
read(unit=itin, fmt=101, err=300) ans
101 format(a)
if (ans .eq. blank) goto 910
read(unit=ans, fmt=100) i4
100 format(bn,i16)
if (.not. verify) goto 900
write(unit=itout, fmt=202) i4
202 format(1x,6hValue ,i16,5h ok? ,$)
read(unit=itin, fmt=101) qr
if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300
900 return
910 i4 = defalt
goto 900
end
subroutine ttini2(qu, i2, defalt, itin,itout, verify)
c Terminal entry of integer*2 variable - F77 version
c qu - Character string prompt que,80 C long.
c i2 - variable to be input
c defalt - default value if <cr> only given
c itin,itout - read, write LUN's of terminal
c verify - if true, asks for verification of input
c if false no verification
c typical call:
c call ttini2(que,k,16,itin,itout,verify)
c MEG/ Feb 86; SUN conversion Jul 89/ MEG
c
logical verify
character ans*6, blank*6, qu*80
character que(80)*1, qr*1
integer*2 i2, defalt
blank = ' '
read(unit=qu, fmt=200) que
200 format(80a1)
call lnnobl(que, nq, 80)
300 write(unit=itout, fmt=204)
204 format(1h ,$)
do 40 i = 1, nq
40 write(unit=itout, fmt=201) que(i)
201 format(a,$)
read(unit=itin, fmt=101, err=300) ans
101 format(a)
if (ans .eq. blank) goto 910
read(unit=ans, fmt=100) i2
100 format(bn,i6)
if (.not. verify) goto 900
write(unit=itout, fmt=202) i2
202 format(1x,6hValue ,i6,5h ok? ,$)
read(unit=itin, fmt=101) qr
if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300
900 return
910 i2 = defalt
goto 900
end
subroutine ttinaa(qu, cs, ncs, defalt, itin,itout, verify)
c Terminal entry of character string - F77 version
c qu - 80C string of prompt que
c cs - character string to be input-array of c*1 but can
c be equivalenced to char string in calling pgm
c ncs - # characters in cs
c defalt - default string for cs
c itin,itout - read, write LUN's of terminal
c verify - if true, asks for verification of input
c if false no verification
c typical call:
c call ttinaa(que,fnam,16,defalt,itin,itout,verify)
c MEG/ Dec 85; SUN conversion Jul 89/ MEG
c
logical verify
character css(80)
character qu*80
character que(80)*1, qr*1
c character defalt(*), cs(*) replaced with following:
character*(*) cs,defalt
read(unit=qu, fmt=200) que
200 format(80a1)
call lnnobl(que, nq, 80)
300 write(unit=itout, fmt=204)
204 format(' ',$)
do 40 i = 1, nq
40 write(unit=itout, fmt=201) que(i)
201 format(a1,$)
read(unit=itin, fmt=101, err=300) css
101 format(80a1)
call lnnobl(css, kcss, 80)
c <cr> go to default
if (kcss .eq. 0) goto 910
write(cs,'(80a)')(css(i),i=1,ncs)
if (.not. verify) goto 900
write(unit=itout, fmt=202) (css(i),i = 1, ncs)
202 format(1x,17hCharacter string /1x,80a:)
write(unit=itout, fmt=203)
203 format(5h ok? ,$)
read(unit=itin, fmt=101) qr
if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300
900 return
910 cs = defalt
goto 900
end
c Program dg4ztcandzbr
c-------------------------------------------------------------------------
c---------> Function: calculate gravity anomaly from
c 2 density functions with known depth to interface between
c density functions and known depth to bottom.
c-------------------------------------------------------------------------
c
c-------------------------------------------------------------------------
c---------> Parameters and notes:
c Model is pile of vertical axis cylinders on coincident axis
c Requires a 2 files of {ztop,zbot,drho,radius}, max 101 layers each
c describing the density contrast and radius of the pile of cyls
c for the shallow and deep density function
c
c and a third file of gravity anomaly and and r/ro to use if total
c model gravity is less than specified amount. I.e., if gravity anom
c aly was 9 mgal, table might be 3.0, 1.0, 9.0,0.5,
c meaning to use whole cylinder(s) gravity effect for first 3 mgal
c and use factor from r / r0 =.5 for each cylinder gravity effect
c from 3mgal to 9 mgal Note factor nonlinear but r/r0=1.0 is whole
c cylinder and r/r0 = 0.0 is half the cylinder; r is radius to plane
c cutting cylinder, ro is radius of cylinder
c This would be appropriate for a well over the range front fault
c where 3mgal due to fill over pediment, rest due to fill against
c hanging wall
c-------------------------------------------------------------------------
c date programmer event
c-------------------------------------------------------------------------
c 26 Dec 96 MEG created
c 27 Apr 99 MEG cylinder factor as fn of grav anomaly added
c
c*-----------------------------------------------------------------------*/
c*Disclaimer: Although this program has been used by the USGS, no */
c*warranty, expressed or implied, is made by the USGS or the United */
c*States Government as to the accuracy and functioning of the program */
c*and related program material nor shall the fact of distribution */
c*constitute any such warranty, and no responsibility is assumed by the */
c*USGS in connection therewith. */
c*-----------------------------------------------------------------------*/
c
character filnam1*50, dflt50*50, qu*80, fmt*50, ids*20
real rhoeff, gxsl, gxcy
real rho1(101,4), rho2(101,4), zbar1(101), zbar2(101), fvsdg(21,2)
real fdg(6), fro(6)
logical verify
data fdg/0.500,0.803,0.915,0.966,0.991,1.000/
data fro/0.0,0.2,0.4,0.6,0.8,1.0/
c
c
c
itin = 5
itout = 6
io2 = 11
verify = .false.
odpth = 0.0
odz = 0.1
rhoeff = 0.0
gxcy = 0.0
gxslc = 0.0
c
c enter identification string
c
qu='Enter identification string (<21c) [dpth2fun estimate]-->'
dflt50 = 'dpth2fun estimate'
call ttinaa(qu, ids, 20, dflt50, itin, itout, verify)
c
c get filename & read file of depths to layers, drho, radius
c
qu='Enter shallow density contrast filename [dshlow.log]-->'
dflt50 = 'dshlow.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
nlayrs1 = 1
write(itout,202)' Shallow density function'
202 format(a)
write(itout,208)
208 format(6x,'ztop',6x,'zbot',6x,'drho',7x,'rad',6x,'zbar',4x,'cum dg',
1 2x,'cum slab')
500 read(io2,fmt,end=510) (rho1(nlayrs1,i), i=1,4)
zbar1(nlayrs1) = (rho1(nlayrs1,1)+rho1(nlayrs1,2))*0.5
rhoeff = rhoeff + rho1(nlayrs1,3)*(rho1(nlayrs1,2)-rho1(nlayrs1,1))
zt = rho1(nlayrs1,1)
zb = rho1(nlayrs1,2)
dr = rho1(nlayrs1,3)
rad = rho1(nlayrs1,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad, 1.0)
gxslc = gxslc + 41.908846 * dr * (zb-zt)
write(itout,207) (rho1(nlayrs1,i), i=1,4), zbar1(nlayrs1), gxcy, gxslc
207 format(7f10.3)
nlayrs1 = nlayrs1 + 1
go to 500
510 nlayrs1= nlayrs1 - 1
rhoeff = rhoeff / (rho1(nlayrs1,2) - rho1(1,1))
gxsl = 41.908846 * rhoeff * (rho1(nlayrs1,2) - rho1(1,1))
write(itout,204) rhoeff, gxsl, gxcy
204 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
qu='Enter deeper density contrast & radius filename [ddeep.log]-->'
dflt50 = 'ddeep.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
rhoeff = 0.0
gxcy = 0.0
gxslc = 0.0
nlayrs2 = 1
write(itout,202)' Deep density function'
write(itout,209)
209 format(6x,'ztop',6x,'zbot',6x,'drho',7x,'rad',6x,'zbar',4x,'cum dg',
1 2x,'cum slab')
502 read(io2,fmt,end=512) (rho2(nlayrs2,i), i=1,4)
zbar2(nlayrs2) = (rho2(nlayrs2,1)+rho2(nlayrs2,2))*0.5
rhoeff = rhoeff + rho2(nlayrs2,3)*(rho2(nlayrs2,2)-rho2(nlayrs2,1))
zt = rho2(nlayrs2,1)
zb = rho2(nlayrs2,2)
dr = rho2(nlayrs2,3)
rad = rho2(nlayrs2,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad, 1.0)
gxslc = gxslc + 41.908846 * dr * (zb-zt)
write(itout,211) (rho2(nlayrs2,i), i=1,4), zbar2(nlayrs2), gxcy, gxslc
211 format(7f10.3)
nlayrs2 = nlayrs2 + 1
go to 502
512 nlayrs2 = nlayrs2 - 1
rhoeff = rhoeff / (rho2(nlayrs2,2) - rho2(1,1))
gxsl = 41.908846 * rhoeff * (rho2(nlayrs2,2) - rho2(1,1))
write(itout,203) rhoeff, gxsl, gxcy
203 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
c
c get cylinder factor vs gravity effect function
c
qu='Enter gravity anomaly vs cylinder factor filename
1 [fctrvsdg.dat]-->'
dflt50 = 'fctrvsdg.dat'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(2f10.3)]-->'
dflt50 = '(2f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
nftr = 1
write(itout,202)' Gravity anomaly vs cylinder factor function'
write(itout,212)
212 format(4x,'deltag',6x,'R/Ro',4x,'factor')
503 read(io2,fmt,end=513) fvsdg(nftr,1), roverro
call interp(fro,fdg,6,4,roverro,fvsdg(nftr,2))
write(itout,211) fvsdg(nftr,1), roverro, fvsdg(nftr,2)
nftr = nftr + 1
go to 503
513 nftr = nftr - 1
c
c get parameters
c
ozct = 0.0
ozmax = 0.0
300 write(qu,200) ozct
200 format('Enter depth to interface, km [',f8.3,']-->')
call ttinr4(qu, zct, ozct, itin, itout, verify)
write(qu,205) ozmax
205 format('Enter depth to bottom, km [',f8.3,']-->')
call ttinr4(qu, zmax, ozmax, itin, itout, verify)
c
c calculate gravity for zTc and ztot
c
c write(itout,202)' ztot deltag'
c call dgtot(0.0, zct, zmax, nlayrs1, rho1, nlayrs2, rho2, dg)
call dgtot1(zct, zmax, nlayrs1, rho1, nlayrs2, rho2,
1 nftr, fvsdg, dg)
write(itout,'(a,a,2f9.3)') ids, ' ztot, deltag: ', zmax, dg
50 continue
write(itout,'(1x)')
ozct = zct
ozmax = zmax
qu = 'Enter 1 to quit, 0 to recalculate [0]'
call ttini4(qu, iquit, 0, itin, itout, verify)
if(iquit.eq.0) go to 300
stop
end
function stkcrcyl(h1,h2,nlayrs,rho,c,nftr,fvsdg)
c
c calc dg of stack of coincident axis vert cylinders from depth h1
c down to depth h2 using array rho for ztop,zbot,drho,rad
c NOTE rho(nlayrs,2) must be greater than H;
c c is offset for model grav dg, removed before return
c nftr is # entries in factor vs dg table
c fvsdg is factor vs dg table dg first entry, f is 2nd
c MEG mod for factor Mar99
c
real h1, h2, rho(101,4), fvsdg(21,4), c
real*8 dg
integer nlayrs, nftr
dg = c * 1.0d0
do 50 i=1,nlayrs
if(h1.ge.rho(i,2)) go to 50
zt = rho(i,1)
zb = rho(i,2)
dr = rho(i,3)
rad = rho(i,4)
fctr = 1.0
do 70 j=1,nftr
if(dg.ge.fvsdg(j,1)) then
fctr = fvsdg(j,2)
go to 72
endif
70 continue
72 continue
c write(*,'(a,2f10.3)')'factor ',fctr,dg
if(h1.ge.rho(i,1).and.h1.lt.rho(i,2)) then
if(h2.lt.rho(i,2)) go to 57
dg = dg + cylfn(h1,zb,dr,rad,fctr)
go to 50
endif
if(h2.lt.rho(i,2)) go to 55
dg = dg + cylfn(zt,zb,dr,rad,fctr)
50 continue
go to 60
55 dg = dg + cylfn(zt,h2,dr,rad,fctr)
go to 60
57 dg = dg + cylfn(h1,h2,dr,rad,fctr)
60 stkcrcyl = dg - c
return
end
function cylfn(ztop,zbot,drho,rad,f)
c gravity effect of vert. axis circular cylinder at pt on axis
c multiplied by factor f
real*8 t
t = zbot - ztop + sqrt(rad*rad + ztop*ztop)
1 - sqrt(rad*rad + zbot*zbot)
cylfn = 41.908846 * drho * t * f
return
end
subroutine dgtot1(zct, ztot, nlayrs1, rho1, nlayrs2, rho2,
1 nftr, fvsdg, dg)
c
c subroutine to calculate gravity effect of a stack coaxial cylinders
c with density fn rho1 to depth zct, rho2 from zct to zbot
c array fvsdg gives proportion of cylinder (<=1.0) to use if grav
c effect is >= to fvsdg(i,1)which is a negative #
c MEG/USGS/Dec96; Mar 99 fvsdg added
c
dimension rho1(101,4),rho2(101,4), fvsdg(21,2)
dg = stkcrcyl(0.0,zct,nlayrs1,rho1,0.0,nftr,fvsdg)
dg = dg + stkcrcyl(zct,ztot,nlayrs2,rho2,dg,nftr,fvsdg)
return
end
0.00 0.500
0.20 0.803
0.40 0.915
0.60 0.966
0.80 0.991
1.00 1.000
c Program dpth2fun
c-------------------------------------------------------------------------
c---------> Function: calculate depth to bottom from gravity anomaly from
c 2 density functions with known depth to interface between
c density functions, and a given table of proportion of each
c cylinder to use as function of gravity anomaly. This is to
c obtain better estimates near the edge of basin
c-------------------------------------------------------------------------
c
c-------------------------------------------------------------------------
c---------> Parameters and notes:
c Model is pile of vertical axis cylinders on coincident axis
c Requires a 2 files of {ztop,zbot,drho,radius}, max 101 layers each
c describing the density contrast and radius of the pile of cyls
c for the shallow and deep density function
c and a third file of gravity anomaly and r/ro to use if total
c model gravity is less than specified amount. I.e., if gravity anom
c aly was 9 mgal, table might be 3.0, 1.0, 9.0,0.5,
c meaning to use whole cylinder(s) gravity effect for first 3 mgal
c and use factor from r / r0 =.5 for each cylinder gravity effect
c from 3mgal to 9 mgal Note factor nonlinear but r/r0=1.0 is whole
c cylinder and r/r0 = 0.0 is half the cylinder; r is radius to plane
c cutting cylinder, ro is radius of cylinder
c This would be appropriate for a well over the range front fault
c where 3mgal due to fill over pediment, rest due to fill against
c hanging wall
c-------------------------------------------------------------------------
c date programmer event
c-------------------------------------------------------------------------
c 26 Dec 96 MEG created
c 11 Dec 96 MEG dpth 2 bottom by bisection added
c 23 Mar 99 MEG cylinder factor as fn of grav anomaly added
c*-----------------------------------------------------------------------*/
c*Disclaimer: Although this program has been used by the USGS, no */
c*warranty, expressed or implied, is made by the USGS or the United */
c*States Government as to the accuracy and functioning of the program */
c*and related program material nor shall the fact of distribution */
c*constitute any such warranty, and no responsibility is assumed by the */
c*USGS in connection therewith. */
c*-----------------------------------------------------------------------*/
c
character filnam1*50, dflt50*50, qu*80, fmt*50, ids*20
real rhoeff, gxsl, gxcy
real rho1(101,4), rho2(101,4), zbar1(101), zbar2(101), fvsdg(21,2)
real fdg(6), fro(6)
logical verify
data fdg/0.500,0.803,0.915,0.966,0.991,1.000/
data fro/0.0,0.2,0.4,0.6,0.8,1.0/
c
c
c
itin = 5
itout = 6
io2 = 11
verify = .false.
odpth = 0.0
odz = 0.1
rhoeff = 0.0
gxcy = 0.0
gxslc = 0.0
c
c enter identification string
c
qu='Enter identification string (<21c) [dpth2fun estimate]-->'
dflt50 = 'dpth2fun estimate'
call ttinaa(qu, ids, 20, dflt50, itin, itout, verify)
c
c get filename & read file of depths to layers, drho, radius
c
qu=
1'Enter shallow density contrast & radius filename [dshlow.log]-->'
dflt50 = 'dshlow.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
nlayrs1 = 1
write(itout,202)' Shallow density function'
202 format(a)
write(itout,208)
208 format(6x,'ztop',6x,'zbot',6x,'drho',7x,'rad',6x,'zbar',4x,'cum dg',
1 2x,'cum slab')
500 read(io2,fmt,end=510) (rho1(nlayrs1,i), i=1,4)
zbar1(nlayrs1) = (rho1(nlayrs1,1)+rho1(nlayrs1,2))*0.5
rhoeff = rhoeff + rho1(nlayrs1,3)*(rho1(nlayrs1,2)-rho1(nlayrs1,1))
zt = rho1(nlayrs1,1)
zb = rho1(nlayrs1,2)
dr = rho1(nlayrs1,3)
rad = rho1(nlayrs1,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad, 1.0)
gxslc = gxslc + 41.908846 * dr * (zb-zt)
write(itout,207) (rho1(nlayrs1,i), i=1,4), zbar1(nlayrs1), gxcy, gxslc
207 format(7f10.3)
nlayrs1 = nlayrs1 + 1
go to 500
510 nlayrs1= nlayrs1 - 1
rhoeff = rhoeff / (rho1(nlayrs1,2) - rho1(1,1))
gxsl = 41.908846 * rhoeff * (rho1(nlayrs1,2) - rho1(1,1))
write(itout,204) rhoeff, gxsl, gxcy
204 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
qu='Enter deeper density contrast & radius filename [ddeep.log]-->'
dflt50 = 'ddeep.log'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(4f10.3)]-->'
dflt50 = '(4f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
rhoeff = 0.0
gxcy = 0.0
gxslc = 0.0
nlayrs2 = 1
write(itout,202)' Deep density function'
write(itout,209)
209 format(6x,'ztop',6x,'zbot',6x,'drho',7x,'rad',6x,'zbar',4x,'cum dg',
1 2x,'cum slab')
502 read(io2,fmt,end=512) (rho2(nlayrs2,i), i=1,4)
zbar2(nlayrs2) = (rho2(nlayrs2,1)+rho2(nlayrs2,2))*0.5
rhoeff = rhoeff + rho2(nlayrs2,3)*(rho2(nlayrs2,2)-rho2(nlayrs2,1))
zt = rho2(nlayrs2,1)
zb = rho2(nlayrs2,2)
dr = rho2(nlayrs2,3)
rad = rho2(nlayrs2,4)
gxcy = gxcy + cylfn(zt,zb,dr,rad, 1.0)
gxslc = gxslc + 41.908846 * dr * (zb-zt)
write(itout,211) (rho2(nlayrs2,i), i=1,4), zbar2(nlayrs2), gxcy, gxslc
211 format(7f10.3)
nlayrs2 = nlayrs2 + 1
go to 502
512 nlayrs2 = nlayrs2 - 1
rhoeff = rhoeff / (rho2(nlayrs2,2) - rho2(1,1))
gxsl = 41.908846 * rhoeff * (rho2(nlayrs2,2) - rho2(1,1))
write(itout,203) rhoeff, gxsl, gxcy
203 format(/'Thickness weighted mean density=',f8.3/
1 'Max.grav.effect infinite slab with this density=',f9.3/
2 'Max.grav.effect cylinders of density model=',f9.3/)
c
c get cylinder factor vs gravity effect function
c
qu='Enter gravity anomaly vs cylinder factor filename
1 [fctrvsdg.dat]-->'
dflt50 = 'fctrvsdg.dat'
call ttinaa(qu, filnam1, 50, dflt50, itin, itout, verify)
open(io2, file=filnam1, status = 'old')
qu = 'Enter data format (<51c)[(2f10.3)]-->'
dflt50 = '(2f10.3)'
call ttinaa(qu, fmt, 50, dflt50, itin, itout, verify)
nftr = 1
write(itout,202)' Gravity anomaly vs cylinder factor function'
write(itout,212)
212 format(4x,'deltag',6x,'R/Ro',4x,'factor')
503 read(io2,fmt,end=513) fvsdg(nftr,1), roverro
call interp(fro,fdg,6,4,roverro,fvsdg(nftr,2))
write(itout,211) fvsdg(nftr,1), roverro, fvsdg(nftr,2)
nftr = nftr + 1
go to 503
513 nftr = nftr - 1
c
c get parameters
c
ozct = 0.0
ozmax = 0.0
odgobs = 0.0
300 write(qu,200) ozct
200 format('Enter depth to interface, km [',f8.3,']-->')
call ttinr4(qu, zct, ozct, itin, itout, verify)
write(qu,205) ozmax
205 format('Enter max depth, km [',f8.3,']-->')
call ttinr4(qu, zmax, ozmax, itin, itout, verify)
write(qu,210) odgobs
210 format('Enter grav.anomaly [',f8.3,']-->')
call ttinr4(qu, dgobs, odgobs, itin, itout, verify)
zacc = 0.001
njmax = 40
c
c calculate depth to bottom by bisection
c
call dgtot1(zct, zmax, nlayrs1, rho1, nlayrs2, rho2,
1 nftr, fvsdg, dg)
zmid = zct
call dgtot1(zct, zmid, nlayrs1, rho1, nlayrs2, rho2,
1 nftr, fvsdg, dgmid)
f = dgobs - dg
fmid = dgobs - dgmid
c write(*,'(a,2f9.3)')'f,fmid:',f,fmid
if(f*fmid.ge.0.0) then
write(itout,'(a,f9.3,a,f9.3)')'no solution between dg(zct)',
1 dgmid,'and dg(zmax)',dg
go to 302
endif
if (f.lt.0.0) then
rtbis = zmax
dz = zct - zmax
else
rtbis = zct
dz = zmax - zct
endif
do 50 i=1,njmax
dz = dz * 0.5
zmid = rtbis + dz
call dgtot1(zct, zmid, nlayrs1, rho1, nlayrs2, rho2,
1 nftr, fvsdg, dgmid)
fmid = dgobs - dgmid
c write(*,'(a,f9.3)')'fmid:',fmid
if(fmid.le.0.0) rtbis = zmid
if(abs(dz).lt.zacc.or.abs(fmid).lt.0.001) go to 304
50 continue
304 write(itout,'(a,a,f9.3,a,f9.3)') ids,'depth to bottom:',zmid,
1 ' misfit:', fmid
write(itout,'(1x)')
ozct = zct
ozmax = zmax
odgobs = dgobs
302 qu = 'Enter 1 to quit, 0 to recalculate [0]'
call ttini4(qu, iquit, 0, itin, itout, verify)
if(iquit.eq.0) go to 300
stop
end
function stkcrcyl(h1,h2,nlayrs,rho,c,nftr,fvsdg)
c
c calc dg of stack of coincident axis vert cylinders from depth h1
c down to depth h2 using array rho for ztop,zbot,drho,rad
c NOTE rho(nlayrs,2) must be greater than H;
c c is offset for model grav dg, removed before return
c nftr is # entries in factor vs dg table
c fvsdg is factor vs dg table dg first entry, f is 2nd
c MEG mod for factor Mar99
c
real h1, h2, rho(101,4), fvsdg(21,4), c
real*8 dg
integer nlayrs, nftr
dg = c * 1.0d0
do 50 i=1,nlayrs
if(h1.ge.rho(i,2)) go to 50
zt = rho(i,1)
zb = rho(i,2)
dr = rho(i,3)
rad = rho(i,4)
fctr = 1.0
do 70 j=1,nftr
if(dg.ge.fvsdg(j,1)) then
fctr = fvsdg(j,2)
go to 72
endif
70 continue
72 continue
c write(*,'(a,2f10.3)')'factor ',fctr,dg
if(h1.ge.rho(i,1).and.h1.lt.rho(i,2)) then
if(h2.lt.rho(i,2)) go to 57
dg = dg + cylfn(h1,zb,dr,rad,fctr)
go to 50
endif
if(h2.lt.rho(i,2)) go to 55
dg = dg + cylfn(zt,zb,dr,rad,fctr)
50 continue
go to 60
55 dg = dg + cylfn(zt,h2,dr,rad,fctr)
go to 60
57 dg = dg + cylfn(h1,h2,dr,rad,fctr)
60 stkcrcyl = dg - c
return
end
function cylfn(ztop,zbot,drho,rad,f)
c gravity effect of vert. axis circular cylinder at pt on axis
c multiplied by factor f
real*8 t
t = zbot - ztop + sqrt(rad*rad + ztop*ztop)
1 - sqrt(rad*rad + zbot*zbot)
cylfn = 41.908846 * drho * t * f
return
end
subroutine dgtot1(zct, ztot, nlayrs1, rho1, nlayrs2, rho2,
1 nftr, fvsdg, dg)
c
c subroutine to calculate gravity effect of a stack coaxial cylinders
c with density fn rho1 to depth zct, rho2 from zct to zbot
c array fvsdg gives proportion of cylinder (<=1.0) to use if grav
c effect is >= to fvsdg(i,1)which is a negative #
c MEG/USGS/Dec96; Mar 99 fvsdg added
c
dimension rho1(101,4),rho2(101,4), fvsdg(21,2)
dg = stkcrcyl(0.0,zct,nlayrs1,rho1,0.0,nftr,fvsdg)
dg = dg + stkcrcyl(zct,ztot,nlayrs2,rho2,dg,nftr,fvsdg)
return
end