subroutine cylind(xq,zq,a,rho,xp,zp,gx,gz)
c
c  Subroutine CYLINDer calculates x and z components of gravitational
c  attraction due to cylinder lying parallel to y axis.
c
c  Input parameters:
c    Point of observation is (xp,zp).  Axis of cylinder penetrates 
c    x,z plane at (xq,zq).  Radius of cylinder is a and density is 
c    rho.  Density in kg/(m**3).  All distance parameters in km. 
c
c  Output parameters:
c    Components of gravitational attraction (gx,gz) in mGal.
c
      real km2m
      data gamma/6.67e-11/,si2mg/1.e5/,pi/3.14159265/,km2m/1.e3/
      rx=xp-xq
      rz=zp-zq
      r2=rx**2+rz**2
      if(r2.eq.0.)pause 'CYLIND:  Bad argument detected.'
      tmass=pi*(a**2)*rho
      gx=-2.*gamma*tmass*rx/r2
      gz=-2.*gamma*tmass*rz/r2
      gx=gx*si2mg*km2m
      gz=gz*si2mg*km2m
      return
      end