function gfilt(kx,ky,z1,z2)
c
c  Function GFILT calculates the value of the gravitational
c  earth filter at a single (kx,ky) location.
c
c  Input parameters:
c    kx - the wavenumber coordinate in the kx direction, in 
c         units of 1/km.
c    ky - the wavenumber coordinate in the ky direction, in 
c         units of 1/km.
c    z1 - the depth to the top of the layer, in km.
c    z2 - the depth to the bottom of the layer, in km.
c
c  Output parameters:
c    gfilt - the value of the earth filter.
c
      real kx,ky,k
      data pi/3.14159265/,gamma/6.67e-11/
      k=sqrt(kx**2+ky**2)
      if(k.eq.0.)then
         gfilt=2.*pi*gamma*(z2-z1)
         else
            gfilt=2.*pi*gamma*(exp(-k*z1)-exp(-k*z2))/k
            end if
      return
      end