************************************************************************
*  Computes apparent horizon edge flag and differencing stencil flag
*  from the mask function.
*     Randall R. Correll
*     Mijan F. Huq
*     Scott A. Klasky
*     Richard A. Matzner
*
*     COPYRIGHT 1995-1996
************************************************************************
      subroutine flags(nx,ny,nz,nlmax,xi,yj,zk,fx,fy,fz,mask,m,d,nl)

      implicit none

      integer   nx,ny,nz,nlmax,nl
      real*8    mask(nx,ny,nz)
      integer   xi(nlmax),yj(nlmax),zk(nlmax)
      real*8    fx(nlmax),fy(nlmax),fz(nlmax)
      integer   m(3*nlmax)
      real*8    d(3*nlmax)

c  local variables

      integer   i, j, k, l
      integer   nlxmax,nlymax,nlzmax
		real*8    zero, one
		parameter ( zero = 0.0e0, one = 1.0e0 )

c  Computes lists of grid points for one-sided stencils.
c  xi, xj, and zk are the grid indices for each point and
c  fx, fy, and fz tells if the differencing should be
c  forward (=1.0), centered (=0.0), or backward (=-1.0).

c  Initialize grid lists and difference types

      do l = 1, nlmax

         xi(l) = 0
         yj(l) = 0
         zk(l) = 0

         fx(l) = zero 
         fy(l) = zero
         fz(l) = zero

      end do

c  Compute lists and difference types at interior grid points

      l = 0

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

      if (mask(i,j,k).eq.one .and. 
     $   ((mask(i+1,j,k).eq.zero).or.(mask(i-1,j,k).eq.zero).or.
     $    (mask(i,j+1,k).eq.zero).or.(mask(i,j-1,k).eq.zero).or.
     $    (mask(i,j,k+1).eq.zero).or.(mask(i,j,k-1).eq.zero)))
     $                                 then

         l = l + 1
         xi(l) = i
         yj(l) = j
         zk(l) = k

         if ( mask(i+1,j,k).eq.0.0e0 ) fx(l) = -1.0e0
         if ( mask(i-1,j,k).eq.0.0e0 ) fx(l) =  1.0e0
         if ( mask(i,j+1,k).eq.0.0e0 ) fy(l) = -1.0e0
         if ( mask(i,j-1,k).eq.0.0e0 ) fy(l) =  1.0e0
         if ( mask(i,j,k+1).eq.0.0e0 ) fz(l) = -1.0e0
         if ( mask(i,j,k-1).eq.0.0e0 ) fz(l) =  1.0e0

      end if


      end do
      end do
      end do

      nl = l

c  Builds masks m and d

      do l = 1, nl

         m(l)      = xi(l)
         m(l+nl)   = yj(l)
         m(l+2*nl) = zk(l)
         d(l)      = fx(l)
         d(l+nl)   = fy(l)
         d(l+2*nl) = fz(l)

      end do

      return
      end