subroutine gbox(x0,y0,z0,x1,y1,z1,x2,y2,z2,rho,g)
c
c  Subroutine GBOX computes the vertical attraction of a 
c  rectangular prism.  Sides of prism are parallel to x,y,z axes,
c  and z axis is vertical down.  
c
c  Input parameters:
c    Observation point is (x0,y0,z0).  The prism extends from x1
c    to x2, from y1 to y2, and from z1 to z2 in the x, y, and z 
c    directions, respectively.  Density of prism is rho.  All 
c    distance parameters in units of km; rho in units of 
c    kg/(m**3). 
c
c  Output parameters:
c    Vertical attraction of gravity, g, in mGal. 
c
      real km2m
      dimension x(2),y(2),z(2),isign(2)
      data isign/-1,1/,gamma/6.670e-11/,twopi/6.2831853/,
     &     si2mg/1.e5/,km2m/1.e3/
      x(1)=x0-x1
      y(1)=y0-y1
      z(1)=z0-z1
      x(2)=x0-x2
      y(2)=y0-y2
      z(2)=z0-z2
      sum=0.
      do 1 i=1,2
         do 1 j=1,2
            do 1 k=1,2
               rijk=sqrt(x(i)**2+y(j)**2+z(k)**2)
               ijk=isign(i)*isign(j)*isign(k)
               arg1=atan2((x(i)*y(j)),(z(k)*rijk))
               if(arg1.lt.0.)arg1=arg1+twopi
               arg2=rijk+y(j)
               arg3=rijk+x(i)
               if(arg2.le.0.)pause 'GBOX:  Bad field point'
               if(arg3.le.0.)pause 'GBOX:  Bad field point'
               arg2=alog(arg2)
               arg3=alog(arg3)
               sum=sum+ijk*(z(k)*arg1-x(i)*arg2-y(j)*arg3)
    1          continue
      g=rho*gamma*sum*si2mg*km2m
      return
      end