!*********************************************************************** 
!  Computes apparent horizon edge flag and differencing stencil flag 
!  from the mask function. 
!*********************************************************************** 
subroutine flags(params,xi,yj,zk,fx,fy,fz,mask,m,d,nl) 
 
implicit none 
integer                       :: nx,ny,nz,nl
real*8,dimension(:)           :: fx,fy,fz 
real*8,dimension(:)           :: params 
integer,dimension(:)          :: xi,yj,zk 
real*8,dimension(:)           :: d
integer,dimension(:)          :: m
real*8,dimension(:,:,:)       :: mask

!  local variables 
 
integer                      :: i,j,k,l,nlxmax,nlymax,nlzmax
real*8 ,parameter            :: zero = 0.0d0
real*8 ,parameter            :: one  = 1.0d0
 
!  Computes lists of grid points for one-sided stencils. 
!  xi, xj, and zk are the grid indices for each point and 
!  fx, fy, and fz tells if the differencing should be 
!  forward (=1.0), centered (=0.0), or backward (=-1.0). 
 
!  Initialize grid lists and difference types 

   nx = int(params(1))
   ny = int(params(2))
   nz = int(params(3))
   xi = 0 
   yj = 0 
   zk = 0 
 
   fx = zero  
   fy = zero 
   fz = zero 
 
!  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.0d0 ) fx(l) = -1.0d0 
         if ( mask(i-1,j,k).eq.0.0d0 ) fx(l) =  1.0d0 
         if ( mask(i,j+1,k).eq.0.0d0 ) fy(l) = -1.0d0 
         if ( mask(i,j-1,k).eq.0.0d0 ) fy(l) =  1.0d0 
         if ( mask(i,j,k+1).eq.0.0d0 ) fz(l) = -1.0d0 
         if ( mask(i,j,k-1).eq.0.0d0 ) fz(l) =  1.0d0 
 
      end if 
 
 
      end do 
      end do 
      end do 
 
      nl = l 
 
!  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 
