************************************************************************
*  Locates apparent horizon and sets inout mask
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996

************************************************************************
      subroutine horizon(nx,ny,nz,g11,g12,g13,g22,g23,g33,
     &           k11,k12,k13,k21,k22,k23,k31,k32,k33,
     &           x,y,z,bha,bhr,inout)

      implicit none

      integer  nx,ny,nz
c     metric:  g_ij
      real*8   g11(nx,ny,nz), g12(nx,ny,nz), g13(nx,ny,nz),
     &                        g22(nx,ny,nz), g23(nx,ny,nz),
     &                                       g33(nx,ny,nz)
c     extrinsic curvature:  K^i_j
      real*8   k11(nx,ny,nz), k12(nx,ny,nz), k13(nx,ny,nz),
     &         k21(nx,ny,nz), k22(nx,ny,nz), k23(nx,ny,nz),
     &         k31(nx,ny,nz), k32(nx,ny,nz), k33(nx,ny,nz)
      real*8   x(nx), y(ny), z(nz), bha(2), bhr(2,3)
      real*8  inout(nx,ny,nz)

c     local variables
      integer  i, j, k
      real*8   test, r

c     compute coordinate radius from singularity

      test = bha(1)

      do k = 1, nz
      do j = 1, ny
      do i = 1, nx

         r = sqrt( (x(i)-bhr(1,1))**2
     &            +(y(j)-bhr(1,2))**2
     &            +(z(k)-bhr(1,3))**2 )
c        use masking funtion over singularity
         if ( r .lt. test ) then
            inout(i,j,k) = 0.0e0
         else
            inout(i,j,k) = 1.0e0
         end if

      end do
      end do
      end do

      return
      end
