/*@@
  @file      AHFinder_mask.F
  @date      November 1998
  @author    Miguel Alcubierre
  @desc 
             Set up horizon mask for black hole excision.
  @enddesc 
@@*/

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

      subroutine AHFinder_mask(CCTK_FARGUMENTS,rhor)

!     Set up the horizon mask.  All this routine does is
!     set up the value of the grid function "ahmask" to zero
!     for points inside the horizon.
!
!     It also finds out the radius of the largest sphere that
!     fits inside the horizon and it centered in (xc,yc,zc).
!     This is useful for simple excision.

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

      integer i,j,k,reduce_handle

      CCTK_REAL zero,one
      CCTK_REAL aux,rhor,rhortemp


!     *****************************
!     ***   DEFINE {zero,one}   ***
!     *****************************

      zero = 0.0D0
      one  = 1.0D0


!     *************************
!     ***   START ROUTINE   ***
!     *************************

!     If inside of horizon (minus a bufferzone) set the mask to 0.
!     To find out how big is the buffer zone  I check the value
!     of the parameter ahf_maskshrink.  This parameter says
!     by which factor of the monopole term to shrink the mask
!     with respect to the real horizon.  Here I use the fact that
!     the function ahfgrid is zero at the horizon, and is linear
!     on the monopole term (which measures the overall scale).
!
!     Notice that if I want a cubical mask, I first find out the radius
!     of the largest sphere that fits inside the horizon and only after
!     that I set the mask.

      rhor = 1.0D10

      call AHFinder_fun(CCTK_FARGUMENTS)

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

               if (ahfgrid(i,j,k).ge.(ahf_maskshrink-one)*c0(0)) then
                  aux = dsqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2
     .                + (z(i,j,k)-zc)**2)
                  if (rhor.gt.aux) rhor = aux
               else if (ahf_maskcube.eq.0) then
                  if (ahfgrid(i,j,k).lt.(0.8D0*ahf_maskshrink-one)*c0(0)) then
                     ahmask(i,j,k) = zero
                  else
                     ahmask(i,j,k) = (ahfgrid(i,j,k) - (0.8D0*ahf_maskshrink-one)*c0(0))
     .                             / (0.2D0*ahf_maskshrink*c0(0))
                  end if
               end if

            end do
         end do
      end do

!     Now find the minumum of rhor across processors.

      call CCTK_ReductionArrayHandle(reduce_handle,"minimum")
      rhortemp = rhor
      call CCTK_ReduceLocalScalar(ierr ,cctkGH,-1,reduce_handle,
     .        rhortemp,rhor,CCTK_VARIABLE_REAL)

      if (ierr.ne.0) then
         call CCTK_WARN(1,"AHFinder_mask:Reduction of rhor failed!");
      end if


!     ************************
!     ***   CUBICAL MASK   ***
!     ************************

!     If we want a cubical mask, set the mask based
!     on the value of rhor.

      if (ahf_maskcube.ne.0) then

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

                  if ((x(i,j,k).gt.(xc-0.5D0*rhor)).and.
     .                (y(i,j,k).gt.(yc-0.5D0*rhor)).and.
     .                (z(i,j,k).gt.(zc-0.5D0*rhor)).and.
     .                (x(i,j,k).lt.(xc+0.5D0*rhor)).and.
     .                (y(i,j,k).lt.(yc+0.5D0*rhor)).and.
     .                (z(i,j,k).lt.(zc+0.5D0*rhor))) then
                     ahmask(i,j,k) = zero
                  end if

               end do
            end do
         end do

      end if


!     ****************************************************
!     ***  APPLY SYMMETRY BOUNDARIES AND SYNCHRONIZE   ***
!     ****************************************************

      call CCTK_SyncGroup(cctkGH,"ahfinder::ahfmask")
      call CartSymGN(ierr,cctkGH,"ahfinder::ahfmask")

!     If desired, copy ahf_mask to einstein mask (emask).

      if (use_mask.eq.1) then
         emask = ahmask
      end if


!     ***************
!     ***   END   ***
!     ***************

      end subroutine AHFinder_mask