path: root/src/AHFinder.F

/*@@
  @file      AHFinder.F
  @date      April 1998
  @author    Miguel Alcubierre
  @desc
             Find an apparent horizon using a spherical harmonic
             expansion and either a minimization or Gundlach's
             fast flow algorithm.

             Cactus4 by Lars Nerger / Gerd Lanfermann
  @enddesc
  @version   $Header$
@@*/

#include "cctk.h"

#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"

      subroutine AHFinder(CCTK_ARGUMENTS)

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      logical status,horizon,gauss,ahfconformal
      logical giveup,report

      integer handle,x_index,y_index,z_index,ierror
      integer ahfso,ahfafter
      integer i,j,l,m
      integer NN
      integer mtype
      integer getoutpfx
      integer istat
      integer nhorizon

      CCTK_REAL atime
      CCTK_REAL rmx,r0
      CCTK_REAL intexp_h,intexp2_h,intarea_h
      CCTK_REAL aux,expin,expout,small1,small2
      CCTK_REAL c0_temp(0:18)
      CCTK_REAL maxval(3),minval(3)
      CCTK_REAL rhor

      CCTK_REAL, parameter :: zero = 0.0D0
      CCTK_REAL, parameter :: half = 0.5D0
      CCTK_REAL, parameter :: one  = 1.0D0

      CCTK_REAL, dimension (1:3) :: ahfmass
      CCTK_INT :: save_conformal_state

      character(200) logf

      save ahfafter,ahfso,atime

!     Description of variables:
!
!     status    Surface found?
!
!     horizon   Horizon found?
!
!     gauss     Output Gaussian curvature?
!
!     ahfso     How often we want to find a horizon?
!     ahfafter  After how many iterations start looking for horizons?
!
!     i,l,m     Counters.
!
!     NN        Total number of terms in expansion.
!
!     rmx       Radius of largest sphere that fits in the grid.
!
!     mtype     Type of surface:
!                  0.  No surface found.
!                  1.  Integral of H^2 is positive out and negative in.
!                  2.  Integral of H^2 is negative out and positive in.
!                  3.  Integral of H^2 is negative out and in.
!                  4.  Integral of H^2 is positive out and in.
!
!     nhorizon  Number of horizons to look for.


!     ***************************
!     ***   INITIALIZE MASK   ***
!     ***************************

!     Initialize mask to one, but only on first iteration.
!     Notice that I do this even if we do not want to find
!     a horizon from the very beginning.  This is to make
!     sure that the mask is 1 and not 0 in case some other
!     thorn wants to use it.

      if (cctk_iteration.eq.0 .or. CCTK_EQUALS(ahf_mask, "off")) then
         ahmask = one
      end if


!     ***********************************************************
!     ***   CHECK IF WE WANT TO FIND A HORIZON AT THIS TIME   ***
!     ***********************************************************

      atime = ahf_findaftertime

      if (atime.le.zero) then
         ahfafter = ahf_findafter
         if (cctk_iteration.lt.ahfafter) return
      else
         if (cctk_time.lt.atime) return
         if (ahf_ncall.eq.0) ahfafter = cctk_iteration
      end if

      ahfso = ahf_findevery
      if (mod(int(cctk_iteration-ahfafter),ahfso).ne.0) return


!     *******************************************
!     ***   ADD TO COUNT OF CALLS TO FINDER   ***
!     *******************************************

!     Add to count of calls to finder.  Notice that this counter was
!     initialized to zero, so it must be increased on first call to 1.

      ahf_ncall = ahf_ncall + 1


!     *********************************************************
!     ***   FIND OUT IF WE WANT TO FIND A TRAPPED SURFACE   ***
!     *********************************************************

      if (ahf_trapped_surface.ne.0) then
         find_trapped_surface = .true.
      else
         find_trapped_surface = .false.
      end if


!     ***************************
!     ***   GRID PARAMETERS   ***
!     ***************************

      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)

      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)


!     *********************************************
!     ***   SET SYMMETRIES FOR GRID FUNCTIONS   ***
!     *********************************************

      call AHFinder_SetSym(CCTK_ARGUMENTS)


!     **********************************************
!     ***   FIND OUT ON WHICH PROCESSOR WE ARE   ***
!     **********************************************

      myproc = CCTK_MyProc(cctkGH)


!     *******************************************
!     ***   FIND  {xmn,xmx,ymn,ymx,zmn,zmx}   ***
!     *******************************************

      call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d")
      call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d")
      call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d")


!     **************************************************
!     ***   TRANSFORM CONFORMAL TO PHYSICAL METRIC   ***
!     **************************************************

      ahfconformal = .false.

!     The horizon finding routines know nothing about conformal
!     metrics, so we need to make sure that we are working with
!     the physical metric.  I change not only the metric, but
!     also the "conformal state" flag to make sure that any
!     macros I use have always the physical metric.  However,
!     I need to remember the original value of the flag to change
!     it back at the end if necesary.

      if ( CCTK_Equals(metric_type,"static conformal") == 1 ) then
         ahfconformal    = .true.
         save_conformal_state = conformal_state
         conformal_state = 0

         call ConfToPhysInPlace (nx,ny,nz,psi,
     &        gxx,gxy,gxz,gyy,gyz,gzz)

      end if


!     ************************
!     ***   FIND  {lmax}   ***
!     ************************

      lmax = ahf_lmax

      if (lmax.ge.20) then
         write(*,*)
         write(*,*) 'ahf_lmax must be smaller than 20.'
         write(*,*) 'STOPPED IN AHFinder.F'
         write(*,*)
         call CCTK_WARN(0,"ahf_lmax must be smaller than 20.")
      end if


!     ***********************************************
!     ***   CHECK IF WE WANT TO FIND 3 HORIZONS   ***
!     ***********************************************

      if (ahf_find3.eq.0) then
         find3 = .false.
         nhorizon = 1
      else
         find3 = .true.
         nhorizon = 3
      end if


!     ***************************************************
!     ***   ALLOCATE STORAGE FOR COEFFICIENT ARRAYS   ***
!     ***************************************************

      if (ahf_ncall.eq.1) then

         allocate(c0(0:lmax),c0_old(0:lmax))
         allocate(cc(lmax,lmax),cc_old(lmax,lmax))
         allocate(cs(lmax,lmax),cs_old(lmax,lmax))

         c0 = zero
         cc = zero
         cs = zero

         c0_old = zero
         cc_old = zero
         cs_old = zero

         allocate(c0_0(0:lmax),c0_1(0:lmax),c0_2(0:lmax))
         allocate(cc_0(lmax,lmax),cc_1(lmax,lmax),cc_2(lmax,lmax))
         allocate(cs_0(lmax,lmax),cs_1(lmax,lmax),cs_2(lmax,lmax))

         c0_0 = zero
         cc_0 = zero
         cs_0 = zero

         c0_1 = zero
         cc_1 = zero
         cs_1 = zero

         c0_2 = zero
         cc_2 = zero
         cs_2 = zero

      end if

!     ***************************************************
!     ***   START LOOP TO LOOK FOR SEVERAL HORIZONS   ***
!     ***************************************************

!     IMPORTANT:  This is the point where we start the loop
!     to look for more than one horizon. Because of this,
!     anything that is likely to be different for the different
!     horizons should NOT be initialized until after this point!

      ahfmass = 0.0D0

      do mfind=1,nhorizon

!     Say what we are doing.

      if (nhorizon.gt.1) then
         if (mfind.eq.1) then
            write(*,*)
            write(*,*)
            write(*,*) 'AHFinder: MULTIPLE HORIZON MODE!'
         end if
         write(*,*)
         write(*,*) 'Searching for horizon ',mfind-1
      else
         write(*,*)
         write(*,*)
         write(*,*) 'AHFinder: SINGLE HORIZON MODE!'
         write(*,*)
         write(*,*) 'Searching for horizon'
      end if


!     ********************************
!     ***   INITIALIZE  {status}   ***
!     ********************************

      status = .false.


!     *********************************
!     ***   INITIALIZE FIRSTCALLS   ***
!     *********************************

      firstleg = .true.
      firstint = .true.


!     ***************************************
!     ***   FIND OUT THE TYPE OF OUTPUT   ***
!     ***************************************

      call CCTK_FortranString(nfile,out_dir,filestr)
      if (ahf_logfile.ne.0) then
         logfile = .true.
         logf = filestr(1:nfile)//"/ahf_logfile"
         if (find3) then
            logf = trim(logf) // "_" // char(ichar('0') + mfind - 1)
         end if
      else
         logfile = .false.
      end if

      if (ahf_verbose.ne.0) then
         verbose = .true.
      else
         verbose = .false.
      end if

      if (ahf_veryverbose.ne.0) then
         veryver = .true.
      else
         veryver = .false.
      end if

      if (ahf_guessverbose.ne.0) then
         guessverbose = .true.
      else
         guessverbose = .false.
      end if

      if (logfile.and.(myproc.eq.0)) then
         if (logfile) then
            if (myproc.eq.0) then
               open(11,file=logf,form='formatted',status='replace')
               write(11,*)
               write(11,*) 'LOG FILE FOR APPARENT HORIZON FINDER'
               write(11,*)
               write(11,*) 'Note:  During an evolution, this file will only'
               write(11,*) 'contain information about the last time that the'
               write(11,*) 'finder was called.'
               write(11,*)
               close(11)
            end if
         end if
      end if


!     **********************************
!     ***   INITIALIZE  {xc,yc,zc}   ***
!     **********************************

      if (find3) then
         if (mfind.eq.1) then
            xc = ahf_xc_0
            yc = ahf_yc_0
            zc = ahf_zc_0
         else if (mfind.eq.2) then
            xc = ahf_xc_1
            yc = ahf_yc_1
            zc = ahf_zc_1
         else
            xc = ahf_xc_2
            yc = ahf_yc_2
            zc = ahf_zc_2
         end if
      else
         xc = ahf_xc
         yc = ahf_yc
         zc = ahf_zc
      end if

!     Check if the center of the expansion is
!     out of the grid.  If it is, give up.

      giveup = .false.

      if ((xc.gt.xmx).or.(xc.lt.xmn)) then
         write(*,*)
         write(*,*) 'xc is out of the grid, giving up.'
         giveup = .true.
      end if

      if ((yc.gt.ymx).or.(yc.lt.ymn)) then
         write(*,*)
         write(*,*) 'yc is out of the grid, giving up.'
         giveup = .true.
      end if

      if ((zc.gt.zmx).or.(zc.lt.zmn)) then
         write(*,*)
         write(*,*) 'zc is out of the grid, giving up.'
         giveup = .true.
      end if

!     The end of this if statement is very far down,
!     it basically jumps all of the code.  To find
!     it look for the word "giveup" on a comment.

      if (.not.giveup) then


!     **************************
!     ***  FIND SYMMETRIES   ***
!     **************************

!     Find out if we want to move the center.

      if (ahf_offset.ne.0) then
         offset = .true.
      else
         offset = .false.
      end if

      if (ahf_wander.ne.0) then
         wander = .true.
      else
         wander = .false.
      end if

      if ((xc.ne.zero).or.(yc.ne.zero).or.(zc.ne.zero).or.
     .    (wander)) offset = .true.

!     Find out if we are in non-axisymmetric mode.

      if (ahf_phi.ne.0) then
         nonaxi = .true.
      else
         nonaxi = .false.
      end if


      call AHFinder_SetReflections(CCTK_ARGUMENTS)

!     Find the values of {stepx,stepy,stepz} depending
!     on the symmetries.

      stepx = 0
      stepy = 0
      stepz = 0

      if (refx) stepx=1
      if (refy) stepy=1
      if (refz) stepz=1

      if (CCTK_Equals(ahf_octant,"high").eq.1) stepx=3


!     *******************************
!     ***   FIND  {ntheta,nphi}   ***
!     *******************************

      ntheta = ahf_ntheta
      nphi   = ahf_nphi

!     For cartoon, I use only 2 points in the phi direction.
!     I could use one, but using 2 minimizes the changes in
!     the integration routine.

      if (cartoon) then
         nphi = 1
      end if


!     ********************************************************
!     ***   FIND TOTAL NUMBER OF TERMS IN EXPANSION {NN}   ***
!     ********************************************************

      NN = 0

!     Count {xc,yc,zc}.

      if (wander) then
         if (.not.refx) NN = NN + 1
         if (.not.refy) NN = NN + 1
         if (.not.refz) NN = NN + 1
      end if

!     Count c0(l).

      do l=0,lmax,1+stepz
         NN = NN + 1
      end do

      if (nonaxi) then

!        Count cc(l,m).

         do l=1,lmax
            do m=1+stepx,l,1+stepx
               if (stepz*mod(l-m,2).eq.0) NN = NN +1
            end do
         end do

!        Count cs(l,m).

         if (.not.refy) then
            do l=1,lmax
               do m=1,l,1+stepx
                  if (stepz*mod(l-m,2).eq.0) NN = NN + 1
               end do
            end do
         end if

      end if


!     *****************************************
!     ***   FIND RADIUS OF INITIAL SPHERE   ***
!     *****************************************

!     Find largest sphere that fits into the grid, taking
!     into account the position of the center.  Then multiply
!     this number with 0.9 to leave a safety margin at the edge
!     of the grid.

!     Cartoon.

      if (cartoon) then

         xc = zero
         yc = zero

         rmx = min(xmx,zmx-zc,zc-zmn)

!     Octant.

      else if (CCTK_Equals(domain,"octant").eq.1) then

         if ((xc.eq.zero).and.(yc.eq.zero).and.(zc.eq.zero)) then
            rmx = min(xmx,ymx,zmx)
         else if ((xc.eq.zero).and.(yc.eq.zero)) then
            rmx = min(xmx,ymx,zmx-zc,zc)
         else if ((xc.eq.zero).and.(zc.eq.zero)) then
            rmx = min(xmx,ymx-yc,zmx,yc)
         else if ((yc.eq.zero).and.(zc.eq.zero)) then
            rmx = min(xmx-xc,ymx,zmx,xc)
         else if (xc.eq.zero) then
            rmx = min(xmx,ymx-yc,zmx-zc,yc,zc)
         else if (yc.eq.zero) then
            rmx = min(xmx-xc,ymx,zmx-zc,xc,zc)
         else if (zc.eq.zero) then
            rmx = min(xmx-xc,ymx-yc,zmx,xc,yc)
         else
            rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc)
         end if

!     Quadrant.

      else if (CCTK_Equals(domain,"quadrant").eq.1) then
         if ((xc.eq.zero).and.(yc.eq.zero)) then
            rmx = min(xmx,ymx,zmx-zc,zc-zmn)
         else if (xc.eq.zero) then
            rmx = min(xmx,ymx-yc,zmx-zc,yc,zc-zmn)
         else if (yc.eq.zero) then
            rmx = min(xmx-xc,ymx,zmx-zc,xc,zc-zmn)
         else
            rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc-zmn)
         end if

!     Bitant.

      else if (CCTK_Equals(domain,"bitant").eq.1) then

         if (CCTK_Equals(bitant_plane,"xy").eq.1) then

            if (zc.eq.zero) then
               rmx = min(xmx-xc,ymx-yc,xc-xmn,yc-ymn,zmx)
            else
               rmx = min(xmx-xc,ymx-yc,xc-xmn,yc-ymn,zmx-zc,zc)
            end if

         else if (CCTK_Equals(bitant_plane,"xz").eq.1) then

            if (yc.eq.zero) then
               rmx = min(xmx-xc,zmx-zc,xc-xmn,zc-zmn,ymx)
            else
               rmx = min(xmx-xc,zmx-zc,xc-xmn,zc-zmn,ymx-yc,yc)
            end if

         else

            if (xc.eq.zero) then
               rmx = min(ymx-yc,zmx-zc,yc-ymn,zc-zmn,xmx)
            else
               rmx = min(ymx-yc,zmx-zc,yc-ymn,zc-zmn,xmx-xc,xc)
            end if

         end if

!     Full.

      else

         rmx = min(xmx-xc,ymx-yc,zmx-zc,xc-xmn,yc-ymn,zc-zmn)

      end if

      rmx = 0.9*rmx

!     Check if a special radius for the initial sphere
!     has been forced by the parameter file.

      if (find3) then
         if (mfind.eq.1) then
            r0 = ahf_r0_0
         else if (mfind.eq.2) then
            r0 = ahf_r0_1
         else
            r0 = ahf_r0_2
         end if
      else
         r0 = ahf_r0
      end if

      if (r0.ne.zero) then
         if (r0.le.rmx) then
            rmx = r0
         else
            write (*,*)
            write (*,*) 'ahf_r0 is too large, rescaling ...'
         end if
      end if


!     *****************************
!     ***   GET INITIAL GUESS   ***
!     *****************************

!     Initializing flags

      sloppy = (ahf_sloppyguess.eq.1)

      if (cartoon) sloppy = .true.

      flow = (ahf_flow.eq.1)

      inner = (ahf_inner.eq.1)

      if (ahf_guessold.ne.0) then
         guessold = .true.
      else
         guessold = .false.
      end if

      if (ahf_guess_absmin.ne.0) then
         guess_absmin = .true.
      else
         guess_absmin = .false.
      end if

      if (ahf_manual_guess.ne.0) then
         manual_guess = .true.
      else
         manual_guess = .false.
      end if

      if (ahf_minarea.ne.0) then
         minarea = .true.
      else
         minarea = .false.
      end if

!     Get initial guess.

      if (find3) then
         if (mfind.eq.1) then
            status_old = status_old_0
         else if (mfind.eq.2) then
            status_old = status_old_1
         else
            status_old = status_old_2
         end if
      end if

      if (guessold.and.status_old) then

!        Old horizon as initial guess.

         if ((myproc.eq.0).and.verbose) then
            write(*,*)
            write(*,*)  'Using old horizon as initial guess'
         end if

!        If we are only looking for 1 horizon, then we are OK
!        since the arrays (c0,cc,cs) where saved from the last
!        call.  If we are looking for 3 horizons, we saved all
!        3 sets of coefficients, but now we need to copy them
!        back to the right place.

         if (find3) then
            if (mfind.eq.1) then
               c0 = c0_0
               cc = cc_0
               cs = cs_0
            else if (mfind.eq.2) then
               c0 = c0_1
               cc = cc_1
               cs = cs_1
            else
               c0 = c0_2
               cc = cc_2
               cs = cs_2
            end if
         end if

      else if (manual_guess) then

!        Manual guess.

         if ((myproc.eq.0).and.verbose) then
            write(*,*)
            write(*,*)  'Using manual initial guess'
         end if

         c0_temp(0) = ahf_l0_guess
         c0_temp(2) = ahf_l2_guess
         c0_temp(4) = ahf_l4_guess
         c0_temp(6) = ahf_l6_guess
         c0_temp(8) = ahf_l8_guess
         c0_temp(10) = ahf_l10_guess
         c0_temp(12) = ahf_l12_guess
         c0_temp(14) = ahf_l14_guess
         c0_temp(16) = ahf_l16_guess
         c0_temp(18) = ahf_l18_guess

         if (lmax.le.18) then
            do i=0,lmax
               c0(i) = c0_temp(i)
            end do
         else
            do i=0,18
               c0(i) = c0_temp(i)
            end do
         end if

      else

!        Standard initial guess.

         if (.not.flow) then

!           For minimization, call initial guess routine.

            call AHFinder_initguess(CCTK_ARGUMENTS,rmx)

         else

!           For flow, take as initial guess a large sphere.

            c0 = zero
            cc = zero
            cs = zero

            c0_old = zero
            cc_old = zero
            cs_old = zero

            if ((myproc.eq.0).and.verbose) then
               write(*,*)
               write(*,"(A35,ES11.3)")  ' Initial guess sphere with radius: ',rmx
               write(*,"(A13,3ES11.3)") ' centered on:',real(xc),real(yc),real(zc)
            end if

            c0(0) = rmx

         end if

      end if


!     ***************************
!     ***   LOOK FOR HORIZON  ***
!     ***************************

      if (.not.flow) then

!        Minimization algorithm.

         call AHFinder_min(CCTK_ARGUMENTS,NN,status,logf)

      else

!        Flow algorithm.

         if ((ahf_ncall.eq.1).and.(mfind.eq.1)) then
            allocate(hflow0(0:lmax),cflow0(0:lmax),nflow0(0:lmax))
            allocate(hflowc(lmax,lmax),cflowc(lmax,lmax),nflowc(lmax,lmax))
            allocate(hflows(lmax,lmax),cflows(lmax,lmax),nflows(lmax,lmax))
         end if

         call AHFinder_flow(CCTK_ARGUMENTS,rmx,status,logf)
         flow = .false.

      end if


!     ***************************************
!     ***   WRITE MESSAGES AND LOG FILE   ***
!     ***************************************

      if (find3) then
         if (mfind.eq.1) then
            status_old_0 = status
         else if (mfind.eq.2) then
            status_old_1 = status
         else
            status_old_2 = status
         end if
      else
         status_old = status
      end if

      if (status) then

!        Check if the expansion on the surface is small enough
!        for a horizon.  How small is small enough is of course
!        pretty arbitrary, so to reduce the subjectivity I use
!        three separate criteria:
!
!        1) The integral of the square of the expansion at the surface
!           must be smaller than a given threshold (small1).
!
!        2) The mean value of the expansion on the surface must be
!           smaller in absolute value than its standard deviation.
!           This means that the mean value of the expansion is
!           consistent with zero.
!
!        3) The integral of the square of the expansion at the surface
!           must be smaller by a factor of "small2" than the same
!           integral slightly outside and slightly inside.

         horizon = .true.

         small1 = 1.0D-2
         small2 = 0.5D0

!        Surface found.

         call AHFinder_fun(CCTK_ARGUMENTS)
         call AHFinder_exp(CCTK_ARGUMENTS)
         call AHFinder_int(CCTK_ARGUMENTS)

         if (ahf_gaussout.ne.0) then
            call AHFinder_gau(CCTK_ARGUMENTS)
         else
            circ_eq = 0.0D0
            meri_p1 = 0.0D0
            meri_p2 = 0.0D0
         end if

         intexp_h  = intexp
         intexp2_h = intexp2
         intarea_h = intarea

         if (intexp2_h.gt.small1) horizon = .false.

!        Find standard deviation:  sqrt(<exp^2> - <exp>^2)
!        The term inside the parenthesis is supposed to be positive,
!        but may fail to be because of numerical errors in pathological
!        situations, so be careful.

         aux = intexp2_h/intarea_h - (intexp_h/intarea_h)**2

         if (aux.ge.zero) then
            aux = sqrt(aux)
            if (abs(intexp_h/intarea_h)-aux.gt.zero) horizon = .false.
         else
            horizon = .false.
         end if

!        Surface slightly inside

         aux = c0(0)

         c0(0) = aux - half*min(dx,dy,dz)
         call AHFinder_int(CCTK_ARGUMENTS)

         expin = intexp

         if (intexp2_h.gt.small2*intexp2) horizon = .false.

         if (veryver) then
            write(*,*)
            write(*,*) 'Surface slightly inside the minimum surface: '
            write(*,"(A21,ES14.6)") ' Integral of H      =',intexp
            write(*,"(A21,ES14.6)") ' Integral of H^2    =',intexp2
            write(*,"(A21,ES14.6)") ' Mean value of H    =',intexp/intarea
            write(*,"(A21,ES14.6)") ' Mean value of H^2  =',intexp2/intarea
            write(*,*) 'Number of interpolated points:     ',
     .         inside_min_count
            write(*,*) 'Number of those that are negative: ',
     .         inside_min_neg_count
         endif

!        Surface slightly outside.

         c0(0) = aux + half*min(dx,dy,dz)
         call AHFinder_int(CCTK_ARGUMENTS)

         expout = intexp

         if (intexp2_h.gt.small2*intexp2) horizon = .false.

         if (veryver) then
            write(*,*)
            write(*,*) 'Surface slightly outside the minimum surface: '
            write(*,"(A21,ES14.6)") ' Integral of H      =',intexp
            write(*,"(A21,ES14.6)") ' Integral of H^2    =',intexp2
            write(*,"(A21,ES14.6)") ' Mean value of H    =',intexp/intarea
            write(*,"(A21,ES14.6)") ' Mean value of H^2  =',intexp2/intarea
            write(*,*) 'Number of interpolated points:     ',
     .         inside_min_count
            write(*,*) 'Number of those that are negative: ',
     .         inside_min_neg_count
         endif

         c0(0) = aux

!        What kind of horizon?

         if ((expin.le.zero).and.(expout.ge.zero)) then
            mtype = 1
         else if ((expin.ge.zero).and.(expout.le.zero)) then
            mtype = 2
         else if ((expin.le.zero).and.(expout.le.zero)) then
            mtype = 3
         else
            mtype = 4
         end if

!        Messages to screen and logfile

         if (ahf_ReportAlways.eq.1) then
            report = .true.
         else
            report = mtype.eq.1
         end if

!        Report on the details of the surface.

         if (report.and.(myproc.eq.0)) then

!           Open logfile.

            if (logfile) then
               open(11,file=logf,form='formatted',
     .         status='old',position='append')
            end if

!           Find inverse of area.  Be carefull not to divide by zero!
!           Also, a value of 1.0D10 indicates an error in the integration
!           routine, and all integrals end up with the same value.
!           So do not divide by it or I will be confused later.

            if ((intarea_h.ne.zero).and.(intarea.lt.1.0D10)) then
               aux = one/intarea_h
            else
               aux = one
            end if

!           Surface information.  Notice that here I define the
!           surface "mass" as:
!
!           M = sqrt( A / (16 pi) ) = 0.141047396 sqrt(A)

            ahfmass(mfind) = 0.141047396D0*sqrt(intarea_h)

            write(*,*)
            write(*,*) 'Surface found, details below.'
            write(*,*)
            write(*,"(A21,ES14.6)") ' Surface area       =',intarea_h
            write(*,"(A21,ES14.6)") ' Surface mass       =',
     .      0.141047396D0*sqrt(intarea_h)
            write(*,"(A21,ES14.6)") ' Mean value of H    =',aux*intexp_h
            write(*,"(A21,ES14.6)") ' Mean value of H^2  =',aux*intexp2_h
            write(*,"(A21,ES14.6)") ' Standard deviation =',
     .      sqrt(abs(aux*intexp2_h - (aux*intexp_h)**2))
            write(*,*)
            if (ahf_gaussout.ne.0) then
               write(*,"(A11,ES14.6)") ' circ_eq  =',circ_eq
               write(*,"(A11,ES14.6)") ' meri_p1  =',meri_p1
               write(*,"(A11,ES14.6)") ' meri_p2  =',meri_p2
               write(*,*)
            end if

            if (logfile) then
               write(11,*)
               write(11,*) 'Surface found, details below'
               write(11,*)
               write(11,"(A21,ES14.6)") ' Surface area       =',intarea_h
               write(11,"(A21,ES14.6)") ' Surface mass       =',
     .         0.141047396D0*sqrt(intarea_h)
               write(11,"(A21,ES14.6)") ' Mean value of H    =',aux*intexp_h
               write(11,"(A21,ES14.6)") ' Mean value of H^2  =',aux*intexp2_h
               write(11,"(A21,ES14.6)") ' Standard deviation =',
     .         sqrt(abs(aux*intexp2_h - (aux*intexp_h)**2))
               write(11,*)
               if (ahf_gaussout.ne.0) then
                  write(11,"(A11,ES14.6)") ' circ_eq  =',circ_eq
                  write(11,"(A11,ES14.6)") ' meri_p1  =',meri_p1
                  write(11,"(A11,ES14.6)") ' meri_p2  =',meri_p2
                  write(11,*)
               end if
            end if

            if (mtype.eq.1) then
               if (horizon) then
                 write(*,*) 'The surface seems to be an outer horizon.'
                 if (logfile) then
                    write(11,*) 'The surface seems to be an outer horizon.'
                 end if
               else
                 write(*,*) 'The surface is probably an outer horizon, but'
                 write(*,*) 'the mean of the square of the expansion is'
                 write(*,*) 'not small enough. Try a larger lmax, or a'
                 write(*,*) 'smaller grid size.'
                 if (logfile) then
                    write(11,*) 'The surface is probably an outer horizon, but'
                    write(11,*) 'the mean of the square of the expansion is'
                    write(11,*) 'not small enough. Try a larger lmax, or a'
                    write(11,*) 'smaller grid size.'
                 end if
               end if
            else if (mtype.eq.2) then
               if (horizon) then
                 write(*,*) 'The surface seems to be an inner horizon.'
                 if (logfile) then
                    write(11,*) 'The surface seems to be an inner horizon.'
                 end if
               else
                 write(*,*) 'The surface is probably an inner horizon, but'
                 write(*,*) 'the mean of the square of the expansion is'
                 write(*,*) 'not small enough. Try a larger lmax, or a'
                 write(*,*) 'smaller grid size.'
                 if (logfile) then
                    write(11,*) 'The surface is probably an inner horizon, but'
                    write(11,*) 'the mean of the square of the expansion is'
                    write(11,*) 'not small enough. Try a larger lmax, or a'
                    write(11,*) 'smaller grid size.'
                 end if
               end if
            else if (mtype.eq.3)then
               if (horizon) then
                 write(*,*) 'The surface seems to be a marginal horizon'
                 write(*,*) '(the integral of the expansion is negative'
                 write(*,*) 'on both sides.)'
                 if (logfile) then
                    write(11,*) 'The surface seems to be a marginal horizon'
                    write(11,*) '(the integral of the expansion is negative'
                    write(11,*) 'on both sides.)'
                 end if
               else
                 write(*,*) 'The surface seems to be a trapped surface.'
                 if (logfile) then
                    write(11,*) 'The surface seems to be a trapped surface.'
                 end if
               end if
            else if (mtype.eq.4)then
               if (horizon) then
                 write(*,*) 'The surface seems to be a marginal horizon'
                 write(*,*) '(the integral of the expansion is positive'
                 write(*,*) 'on both sides.)'
                 if (logfile) then
                    write(11,*) 'The surface seems to be a marginal horizon'
                    write(11,*) '(the integral of the expansion is positive'
                    write(11,*) 'on both sides.)'
                 end if
               else
                 write(*,*) 'The surface does not seem to be a horizon.'
                 if (logfile) then
                    write(11,*) 'The surface does not seem to be a horizon.'
                 end if
               end if
            end if

!           Shape.

            if (verbose) then

               write(*,*)
               write(*,*) 'Shape coefficients:'

               if (offset) then
                  write(*,*)
                  write(*,"(A6,ES14.6)") ' xc  =',xc
                  write(*,"(A6,ES14.6)") ' yc  =',yc
                  write(*,"(A6,ES14.6)") ' zc  =',zc
               end if

               write(*,*)

               write(*,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0)

               do l=1+stepz,lmax,1+stepz
                  write(*,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l)
               end do

            end if

            if (logfile) then

               write(11,*)
               write(11,*) 'Shape coefficients:'

               if (offset) then
                  write(11,*)
                  write(11,"(A6,ES14.6)") ' xc  =',xc
                  write(11,"(A6,ES14.6)") ' yc  =',yc
                  write(11,"(A6,ES14.6)") ' zc  =',zc
               end if

               write(11,*)

               write(11,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0)

               do l=1+stepz,lmax,1+stepz
                  write(11,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l)
               end do

            end if

            if (nonaxi) then

               if (.not.refy) then

                  if (verbose) then

                     write(*,*)
                     do l=1,lmax
                        do m=1,l
                           if (stepz*mod(l-m,2).eq.0) then
                              write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(',
     .                        l,',',m,')  =',cc(l,m)
                           end if
                        end do
                     end do
                     write(*,*)
                     do l=1,lmax
                        do m=1,l
                           if (stepz*mod(l-m,2).eq.0) then
                              write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(',
     .                        l,',',m,')  =',cs(l,m)
                           end if
                        end do
                     end do

                  end if

                  if (logfile) then
                     write(11,*)

                     do l=1,lmax
                        do m=1,l
                           if (stepz*mod(l-m,2).eq.0) then
                              write(11,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(',
     .                        l,',',m,')  =',cc(l,m)
                           end if
                        end do
                     end do
                     write(11,*)
                     do l=1,lmax
                        do m=1,l
                           if (stepz*mod(l-m,2).eq.0) then
                              write(11,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(',
     .                        l,',',m,')  =',cs(l,m)
                           end if
                        end do
                     end do

                  end if

               else

                  if (verbose) then

                     write(*,*)
                     do l=1,lmax
                        do m=1,l
                           if (stepz*mod(l-m,2).eq.0) then
                              write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(',
     .                        l,',',m,')  =',cc(l,m)
                           end if
                        end do
                     end do

                  end if

                  if (logfile) then

                     write(11,*)
                     do l=1,lmax
                        do m=1,l
                           if (stepz*mod(l-m,2).eq.0) then
                              write(11,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(',
     .                        l,',',m,')  =',cc(l,m)
                           end if
                        end do
                     end do

                  end if

               end if
            end if

!           Close logfile

            write(*,*)

            if (logfile) then
               write(11,*)
               write(11,*) 'END OF LOG FILE'
               write(11,*)
               close(11)
            end if

!        If we found a surface that is not a horizon and we
!        do not want to report all its details, just say that
!        we did not find a horizon.

         else if (myproc.eq.0) then

!           Open logfile.

            if (logfile) then
               open(11,file=logf,form='formatted',
     .         status='old',position='append')
            end if

!           Write messages.

            write(*,*)
            write(*,*) 'No horizon found.'

            if (logfile) then
               write(11,*)
               write(11,*) 'No horizon found.'
               write(11,*)
            end if

!           Close logfile

            write(*,*)

            if (logfile) then
               write(11,*)
               write(11,*) 'END OF LOG FILE'
               write(11,*)
               close(11)
            end if

         end if

!     No surface found.

      else

         mtype = 0
         report = .false.

         if (myproc.eq.0) then

!           Open logfile.

            if (logfile) then
               open(11,file=logf,form='formatted',
     .         status='old',position='append')
            end if

!           Write messages.

            write(*,*)
            write(*,*) 'No horizon found.'

            if (logfile) then
               write(11,*)
               write(11,*) 'No horizon found.'
               write(11,*)
            end if

!           Close logfile

            write(*,*)

            if (logfile) then
               write(11,*)
               write(11,*) 'END OF LOG FILE'
               write(11,*)
               close(11)
            end if

         end if
      end if


!     ******************************************************
!     ***   PREPARING  {ahfgrid3,ahf_exp3}  FOR OUTPUT   ***
!     ******************************************************

      if (find3) then
         call AHFinder_find3(CCTK_ARGUMENTS,mtype)
      end if


!     ****************************
!     ***   WRITE DATA FILES   ***
!     ****************************

      call AHFinder_Output(CCTK_ARGUMENTS,report,status,horizon,mtype,intarea_h)


!     ****************
!     ***   MASK   ***
!     ****************

!     Initialize the sides of the excised region to the center
!     of the expansion.

      if (ahf_ncall.eq.1) then

         if (find3) then

            if (mfind.eq.1) then

               dhole1_xmin = xc
               dhole1_ymin = yc
               dhole1_zmin = zc

               dhole1_xmax = xc
               dhole1_ymax = yc
               dhole1_zmax = zc

            else if (mfind.eq.2) then

               dhole2_xmin = xc
               dhole2_ymin = yc
               dhole2_zmin = zc

               dhole2_xmax = xc
               dhole2_ymax = yc
               dhole2_zmax = zc

            else

               dhole3_xmin = xc
               dhole3_ymin = yc
               dhole3_zmin = zc

               dhole3_xmax = xc
               dhole3_ymax = yc
               dhole3_zmax = zc

            end if

         else

            dhole1_xmin = xc
            dhole1_ymin = yc
            dhole1_zmin = zc

            dhole1_xmax = xc
            dhole1_ymax = yc
            dhole1_zmax = zc

         end if

      end if

!     The mask should only be different from 1 inside a horizon.
!     I therefore need to check first if I do have a horizon, and
!     also if the mask is desired.

      if ((horizon.and.(CCTK_EQUALS(ahf_mask,'strong'))).or.
     .   (status.and.(mtype.eq.1).and.(CCTK_EQUALS(ahf_mask,'weak')))) then

         rhor = 0.0D0

         if (find3) then
            if ((mfind.eq.1).and.(ahf_mask_0.eq.1)) then
               call AHFinder_mask(CCTK_ARGUMENTS,rhor)
            else if ((mfind.eq.2).and.(ahf_mask_1.eq.1)) then
               call AHFinder_mask(CCTK_ARGUMENTS,rhor)
            else if ((mfind.eq.3).and.(ahf_mask_2.eq.1)) then
               call AHFinder_mask(CCTK_ARGUMENTS,rhor)
            end if
         else
            call AHFinder_mask(CCTK_ARGUMENTS,rhor)
         end if

!        Here I alter the parameters for SimpleExcision.  Notice that at
!        the moment, I only allow the excised cube to grow, and never
!        to shrink.  Also, notice that the sides of the cube have been
!        set to 0.57*rhor.  This is because if we want the corners of
!        the cube to touch a sphere of radius rhor, we need to make the
!        sides equal to 1/sqrt(3), which is roughly 0.57.  If you worry
!        that this is too close to the horizon, remember that rhor in fact
!        has already been shrunk by a safety margin.

         aux = 0.57D0*rhor

         if (find3) then

            if (mfind.eq.1) then

!              First hole.

               if (dhole1_xmin.gt.(xc-aux)) dhole1_xmin = xc-aux
               if (dhole1_ymin.gt.(yc-aux)) dhole1_ymin = yc-aux
               if (dhole1_zmin.gt.(zc-aux)) dhole1_zmin = zc-aux

               if (dhole1_xmax.lt.(xc+aux)) dhole1_xmax = xc+aux
               if (dhole1_ymax.lt.(yc+aux)) dhole1_ymax = yc+aux
               if (dhole1_zmax.lt.(zc+aux)) dhole1_zmax = zc+aux

            else if (mfind.eq.2) then

!              Second hole.

               if (dhole2_xmin.gt.(xc-aux)) dhole2_xmin = xc-aux
               if (dhole2_ymin.gt.(yc-aux)) dhole2_ymin = yc-aux
               if (dhole2_zmin.gt.(zc-aux)) dhole2_zmin = zc-aux

               if (dhole2_xmax.lt.(xc+aux)) dhole2_xmax = xc+aux
               if (dhole2_ymax.lt.(yc+aux)) dhole2_ymax = yc+aux
               if (dhole2_zmax.lt.(zc+aux)) dhole2_zmax = zc+aux

            else

!              Third hole.

               if (dhole3_xmin.gt.(xc-aux)) dhole3_xmin = xc-aux
               if (dhole3_ymin.gt.(yc-aux)) dhole3_ymin = yc-aux
               if (dhole3_zmin.gt.(zc-aux)) dhole3_zmin = zc-aux

               if (dhole3_xmax.lt.(xc+aux)) dhole3_xmax = xc+aux
               if (dhole3_ymax.lt.(yc+aux)) dhole3_ymax = yc+aux
               if (dhole3_zmax.lt.(zc+aux)) dhole3_zmax = zc+aux

            end if

         else

            if (dhole1_xmin.gt.(xc-aux)) dhole1_xmin = xc-aux
            if (dhole1_ymin.gt.(yc-aux)) dhole1_ymin = yc-aux
            if (dhole1_zmin.gt.(zc-aux)) dhole1_zmin = zc-aux

            if (dhole1_xmax.lt.(xc+aux)) dhole1_xmax = xc+aux
            if (dhole1_ymax.lt.(yc+aux)) dhole1_ymax = yc+aux
            if (dhole1_zmax.lt.(zc+aux)) dhole1_zmax = zc+aux

         end if

      end if


!     ********************************************************
!     ***   END OF LOOP FOR LOOKING FOR SEVERAL HORIZONS   ***
!     ********************************************************

!     End of "giveup" if statement.

      end if

!     If we are looking for 3 horizons, save values of coefficients.

      if (find3) then
         if (mfind.eq.1) then
            c0_0 = c0
            cc_0 = cc
            cs_0 = cs
         else if (mfind.eq.2) then
            c0_1 = c0
            cc_1 = cc
            cs_1 = cs
         else
            c0_2 = c0
            cc_2 = cc
            cs_2 = cs
         end if
      end if

      write (*,*)

      end do


!     *******************************************************
!     ***   CALCULATE PROPER DISTANCES BETWEEN HORIZONS   ***
!     *******************************************************

      if (find3) then
         call AHFinder_dis(CCTK_ARGUMENTS)
      end if


!     **************************************************
!     ***   TRANSFORM PHYSICAL TO CONFORMAL METRIC   ***
!     **************************************************

!     If necesary, change back the metric to conformal
!     and change back the "conformal_state" flag.

      if (ahfconformal) then

         conformal_state = save_conformal_state

         call PhysToConfInPlace (nx,ny,nz,psi,
     &        gxx,gxy,gxz,gyy,gyz,gzz)

      end if


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

      end subroutine AHFinder