! $Header$

! Take a mask where the boundary has been been marked.
! Return information about the normals in these locations.

#include "cctk.h"
#include "cctk_Parameters.h"

#include "maskvalues.h"

subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk)

  implicit none

  DECLARE_CCTK_FUNCTIONS
  DECLARE_CCTK_PARAMETERS
  
  ! Arguments
  
  ! out: zero for success, nonzero for error
  integer      :: ierr
  
  ! in: array sizes for grid functions
  !     (you can pass in cctk_lsh(:) for these)
  integer      :: ni,nj,nk

  ! in: mask
  CCTK_REAL    :: mask(ni,nj,nk)
  
  ! out: normal directions to use for interpolation
  CCTK_REAL    :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk)

! Internal variables.

  integer i,j,k,ii,jj,kk

  CCTK_REAL sx,sy,sz,smax

! Initialise direction arrays to zero.

  dirx = 0
  diry = 0
  dirz = 0

! Loop over all grid points.

  do k=2,nk-1
     do j=2,nj-1
        do i=2,ni-1

!          Check if the point is in the boundary, if it
!          is not forget about ir.

           if (mask(i,j,k)==MASK_BOUNDARY) then

!             Initialise the vector (sx,sy,sz) to zero.

              sx = 0; sy = 0; sz = 0

!             Loop around all nearest neighbours.

              do kk=-1,1
                 do jj=-1,1
                    do ii=-1,1

!                      If a given neighbour is active, then add its
!                      relative coordinates to (sx,sy,sz).

                       if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then
                          sx = sx + dble(ii)
                          sy = sy + dble(jj)
                          sz = sz + dble(kk)
                       end if

                    end do
                 end do
              end do

!             Find the maximum component of (sx,sy,sz).

              smax = max(dabs(sx),dabs(sy),dabs(sz))

!             Divide (sx,sy,sz) with smax.

              if (smax /= 0.0D0) then
                 sx = sx/smax
                 sy = sy/smax
                 sz = sz/smax
              end if

!             Now we have a real vector where the largest component
!             has absolute value 1.  We want to change this into
!             the closest vector with components given by 0 or +-1.
!             To do this, I add 1/2 to all components (with the
!             correct sign), ang get the integer part.

              if (sx /= 0.0D0) sx = dble(int(dabs(sx)+0.5D0))*sx/dabs(sx)
              if (sy /= 0.0D0) sy = dble(int(dabs(sy)+0.5D0))*sy/dabs(sy)
              if (sz /= 0.0D0) sz = dble(int(dabs(sz)+0.5D0))*sz/dabs(sz)

!             Copy the vector (sx,sy,sz) into the normal arrays.

              dirx(i,j,k) = sx
              diry(i,j,k) = sy
              dirz(i,j,k) = sz

!             Sanity check.  If the check fails then the mask really
!             has a crazy shape.

              if (mask(i+int(sx),j+int(sy),k+int(sz)) /= MASK_ACTIVE) then
                 call CCTK_WARN (0, "Mask boundary layer is too thick")
              end if

           end if
           
        end do
     end do
  end do
  
  ! Success

  ierr = 0
  
end subroutine excision_findnormals