path: root/src/AHFinder_int.F

/*@@
  @file      AHFinder_int.F
  @date      April 1998
  @author    Miguel Alcubierre
  @desc
             Find surface integrals. The integrals are done
             in parallel, but the number of processors is
             assumed to be either the square of an integer
             or a power of two (if this is not the case,
             I just use less processors).
  @enddesc
  @version   $Header$
@@*/

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

      subroutine AHFinder_int(CCTK_ARGUMENTS)

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer i,j,l,m
      integer npt,npp
      integer l_ntheta,l_nphi,theta0,phi0
      integer npoints
      integer auxi,vindex
      integer ierror

      integer interp_handle,coord_system_handle
      integer reduction_handle,param_table_handle
      character(30) options_string
      character(128) operator
      CCTK_INT red_tmp, nchars

      CCTK_POINTER, dimension(3) :: interp_coords
      CCTK_INT, dimension(9) :: in_array_indices
      CCTK_POINTER, dimension(9) :: out_arrays
      CCTK_INT, dimension(9) :: out_array_type_codes


      CCTK_REAL LEGEN
      CCTK_REAL xp,yp,zp,rp
      CCTK_REAL xw,yw,zw
      CCTK_REAL cost,sint,cosp,sinp
      CCTK_REAL trr,ttt,tpp,trt,trp,ttp,ft,fp
      CCTK_REAL dtheta,dphi,dtp,idtheta,idphi,phistart
      CCTK_REAL det,idet
      CCTK_REAL zero,quarter,half,one,two,three,four,pi
      CCTK_REAL sina,cosa,intw,grad,sigma,h0
      CCTK_REAL theta,phi
      CCTK_REAL aux,aux1,aux2

      CCTK_REAL tempv(0:lmax),tempm(lmax,lmax)

      CCTK_REAL, dimension(3,3) :: gupij

      CCTK_REAL, allocatable, dimension(:) :: costheta,sintheta
      CCTK_REAL, allocatable, dimension(:) :: cosphi,sinphi

      CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za,da,exp,gradn
      CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz
      CCTK_REAL, allocatable, dimension(:,:) :: intmask

!     Reduction variables
      CCTK_POINTER_TO_CONST, dimension(1) :: input_array
      CCTK_POINTER, dimension(1) :: reduction_value
      CCTK_INT, dimension(1) ::input_array_type
      
!     Variables to be saved for next call.

      save dtheta,dphi,dtp,idtheta,idphi,phistart
      save npt,npp,l_ntheta,l_nphi,theta0,phi0

!     Description of variables:
!
!     i,j,l,m    Counters.
!
!     theta      Latitude.
!     phi        Longitude.
!
!     npt        number of processors in theta direction.
!     npp        number of processors in phi direction.
!
!     l_ntheta   Local number of grid points in theta direction.
!     l_nphi     Local number of grid points in phi direction.
!
!     theta0     Local origin in theta direction.
!     phi0       Local origin in phi direction.
!
!     npoints    Number of points to interpolate .
!
!     xp         x coordinate from surface centre.
!     yp         y coordinate from surface centre.
!     zp         z coordinate from surface centre.
!
!     rp         Radius.
!     rr         Radius array.
!
!     xa         Array with x coordinates from grid centre.
!     ya         Array with y coordinates from grid centre.
!     za         Array with z coordinates from grid centre.
!
!     exp        Interpolated expansion.
!
!     txx        Array with interpolated gxx metric component.
!     tyy        Array with interpolated gyy metric component.
!     tzz        Array with interpolated gzz metric component.
!     txy        Array with interpolated gxy metric component.
!     txz        Array with interpolated gxz metric component.
!     tyz        Array with interpolated gyz metric component.
!
!     gradn      Array with interpolated norm of gradient of horizon function.
!
!     intmask    Array with interpolated mask.
!
!     cost       cos(theta)
!     sint       sin(theta)
!     cosp       cos(phi)
!     sinp       sin(phi)
!
!     trr        Metric component {r,r}.
!     ttt        Metric component {theta,theta}.
!     tpp        Metric component {phi,phi}.
!     trt        Metric component {r,theta}.
!     trp        Metric component {r,phi}.
!     ttp        Metric component {theta,phi}.
!
!     ft         dr/dtheta
!     fp         dr/dphi
!
!     dtheta     Grid spacing in theta.
!     dphi       Grid spacing in phi.
!     dtp        dtheta*dphi.
!
!     idtheta    1/(2 dtheta)
!     idphi      1/(2 dphi)
!
!     phistart   Origin for phi (normally 0, but for some symmetries -pi/2).
!
!     da         Area element.


!     *******************
!     ***   NUMBERS   ***
!     *******************

      zero    = 0.0D0
      quarter = 0.25D0
      half    = 0.5D0
      one     = 1.0D0
      two     = 2.0D0
      three   = 3.0D0
      four    = 4.0D0

      pi = acos(-one)


!     *****************************************************************
!     ***   TOTAL NUMBER OF PROCESSORS AND LOCAL PROCESSOR NUMBER   ***
!     *****************************************************************

      myproc = CCTK_MyProc(cctkGH)
      nprocs = CCTK_nProcs(cctkGH)


!     ****************************************
!     ***   GET REDUCTION HANDLE FOR SUM   ***
!     ****************************************

      call CCTK_LocalArrayReductionHandle(reduction_handle,"sum")

      if (reduction_handle.lt.0) then
        call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??")
      end if


!     ******************************************************
!     ***   THINGS TO DO ON FIRST CALL TO THIS ROUTINE   ***
!     ******************************************************

      if (firstint) then

         firstint = .false.


!        ******************************************
!        ***   FIND  phistart,dtheta,dphi,dtp   ***
!        ******************************************

         phistart = zero

         dtheta = pi/dble(ntheta)

         if (cartoon) then

            dphi    = zero
            dtp     = dtheta
            idphi   = one
            idtheta = half/dtheta

         else

            dphi = two*pi/dble(nphi)

            if (refz) dtheta = half*dtheta

            if (refx.and.refy) then
               dphi = quarter*dphi
            else if (refx) then
               dphi = half*dphi
               phistart = - half*pi
            else if (refy) then
               dphi = half*dphi
            end if

            dtp = dtheta*dphi

            idtheta = half/dtheta
            idphi   = half/dphi

         end if


!        *********************************************
!        ***   DIVIDE 2D DOMAIN AMONG PROCESSORS   ***
!        *********************************************

!        For cartoon this is trivial.

         if (cartoon) then

            npt = nprocs
            npp = 1

!        Otherwise it is more complicated.  At the moment I
!        only consider numbers of processors that are either
!        the square of an integer, or a power of two.
!        If neither of this is the case, then I use fewer
!        processors to make sure that I always deal with
!        either a square or a power of two.  This can
!        certainly be improved, but it might not be that
!        urgent since in most cases the number of processors
!        will be a power of two.

         else

            aux = sqrt(dble(nprocs))

!           One processor.

            if (nprocs.eq.1) then

               npt = 1
               npp = 1

!           Square of an integer.

            else if (aux.eq.int(aux)) then

               npt = int(aux)
               npp = int(aux)

!           The number of processors is not the square of an integer.
!           We will now check if it is a power of two.

            else

               l = 1
               m = 1

               i = 1

               do while (l*m.lt.nprocs)

                  npt = l
                  npp = m

                  if (mod(i,2).ne.0) then
                     m = 2*m
                  else
                     l = 2*l
                  end if

                  i = i + 1

               end do

!              "nprocs" is a power of two.

               if (l*m.eq.nprocs) then

                  npt = l
                  npp = m

!              "nprocs" is not a power of two.  This is where I
!              would need to do something clever.  At the moment,
!              I will just use fewer processors.  The number of
!              processors I use is the power of two or the square
!              of an integer that is closest to "nprocs".

               else

                  if (npt*npp.lt.int(aux)**2) then
                     npt = int(aux)
                     npp = int(aux)
                  end if

               end if

            end if

         end if

!        Now be careful not to use too many procs and end up
!        with very few points per proc.

         if (ntheta/npt.lt.3) npt = ntheta/3

         if (.not.cartoon) then
            if (nphi/npp.lt.3) npp = nphi/3
         end if

!        Figure out the number of grid points per processor,
!        and the "origin" for a given processor.

         if (myproc.ge.npt*npp) then

            l_ntheta = 0
            l_nphi   = 0

            theta0 = 0
            phi0   = 0

         else

!           First take the number of grid points per processor
!           in a given direction to be equal to the total number
!           of grid points divided by the number of processors
!           in that direction.

            l_ntheta = ntheta/npt
            l_nphi   = nphi/npp

!           And take the "origin" as the corresponding displacement
!           from the real origin.  Notice that mod(myproc,npt)
!           tells my on which "theta" bin I am, while myproc/npt
!           tells me on which "phi" bin I am.

            i = mod(myproc,npt)
            j = myproc/npt

            theta0 = i*l_ntheta
            phi0   = j*l_nphi

!           Now correct all this if the total number of grid points
!           in a given direction was not exactly divisible by the
!           number of processors on that direction.

            if (l_ntheta*npt.ne.ntheta) then

!              Find residue on theta direction.

               l = ntheta - l_ntheta*npt

!              Distribute residue in first few processors, and
!              displace theta0 accordingly.

               if (i.lt.l) then
                  l_ntheta = l_ntheta + 1
                  theta0 = theta0 + i
               else
                  theta0 = theta0 + l
               end if

            end if

            if (l_nphi*npp.ne.nphi) then

!              Find residue on phi direction.

               l = nphi - l_nphi*npp

!              Distribute residue in first few processors, and
!              displace phi0 accordingly.

               if (j.lt.l) then
                  l_nphi = l_nphi + 1
                  phi0 = phi0 + j
               else
                  phi0 = phi0 + l
               end if

            end if

         end if

      end if


!     **************************************
!     ***   ALLOCATE MEMORY FOR ARRAYS   ***
!     **************************************

      allocate(costheta(1:l_ntheta+1))
      allocate(sintheta(1:l_ntheta+1))

      allocate(cosphi(1:l_nphi+1),sinphi(1:l_nphi+1))

      allocate(rr(1:l_ntheta+1,1:l_nphi+1))

      allocate(xa(1:l_ntheta+1,1:l_nphi+1))
      allocate(ya(1:l_ntheta+1,1:l_nphi+1))
      allocate(za(1:l_ntheta+1,1:l_nphi+1))

      allocate(da(1:l_ntheta+1,1:l_nphi+1))
      allocate(exp(1:l_ntheta+1,1:l_nphi+1))
      allocate(gradn(1:l_ntheta+1,1:l_nphi+1))

      allocate(txx(1:l_ntheta+1,1:l_nphi+1))
      allocate(tyy(1:l_ntheta+1,1:l_nphi+1))
      allocate(tzz(1:l_ntheta+1,1:l_nphi+1))
      allocate(txy(1:l_ntheta+1,1:l_nphi+1))
      allocate(txz(1:l_ntheta+1,1:l_nphi+1))
      allocate(tyz(1:l_ntheta+1,1:l_nphi+1))

      allocate(intmask(1:l_ntheta+1,1:l_nphi+1))


!     ********************************
!     ***   INITIALIZE VARIABLES   ***
!     ********************************

!     Initialize error flags.

      interror1 = 0
      interror2 = 0
      interror3 = 0

!     Initialize surface integrals.

      intexp  = zero
      intexp2 = zero
      intarea = zero
      intexpdel2 = zero

!     For flow algorithm, initialize spectral components.

      if (flow) then

         hflow0 = zero
         hflowc = zero
         hflows = zero

         cflow0 = zero
         cflowc = zero
         cflows = zero

         nflow0 = zero
         nflowc = zero
         nflows = zero

      end if


!     *************************************
!     ***   FIND INTERPOLATING POINTS   ***
!     *************************************

      if (myproc.ge.npt*npp) then

         xa = half*(xmx+xmn)
         ya = half*(ymx+ymn)
         za = half*(zmx+zmn)

      else

         do j=1,l_nphi+1

!           Find phi.

            phi = dphi*dble(j-1+phi0) + phistart

!           Find sines and cosines of phi.

            cosp = cos(phi)
            sinp = sin(phi)

            cosphi(j) = cosp
            sinphi(j) = sinp

            do i=1,l_ntheta+1

!              Find theta.

               theta = dtheta*dble(i-1+theta0)

!              Find sines and cosines of theta.

               cost = cos(theta)
               sint = sin(theta)

               costheta(i) = cost
               sintheta(i) = sint

!              Find radius rp.

               rp = c0(0)

               do l=1+stepz,lmax,1+stepz
                  rp = rp + c0(l)*LEGEN(l,0,cost)
               end do

!              Notice how the sum over m is first.  This will allow
!              me to use the recursion relations to avoid having to
!              start from scratch every time.  Also, I sum over all
!              l even if I do not want some terms.  This is because
!              in order to use the recursion relations I need all
!              polynomials.

               if (nonaxi) then
                  do m=1,lmax
                     aux = dble(m)*phi
                     sina = sin(aux)
                     cosa = cos(aux)
                     do l=m,lmax
                        aux = LEGEN(l,m,cost)
                        rp = rp + aux*cc(l,m)*cosa
                        if (.not.refy) then
                           rp = rp + aux*cs(l,m)*sina
                        end if
                     end do
                  end do
               end if

!              Check for negative radius.

               if (rp.le.zero) then
                  interror1 = 1
               end if

!              Find cartesian coordinates.

               aux = rp*sint

               xp = aux*cosp + xc
               yp = aux*sinp + yc
               zp = rp*cost + zc

!              Save arrays.

               rr(i,j) = rp

               xa(i,j) = xp
               ya(i,j) = yp
               za(i,j) = zp

!              Check if we are within bounds.

               if ((xp.gt.xmx).or.(xp.lt.xmn).or.
     .             (yp.gt.ymx).or.(yp.lt.ymn).or.
     .             (zp.gt.zmx).or.(zp.lt.zmn)) then
                  interror2 = 1
               end if

            end do
         end do

      end if


!     *******************************************
!     ***   REDUCE ERRORS ACROSS PROCESSORS   ***
!     *******************************************

      input_array_type(1) = CCTK_VARIABLE_INT
      reduction_value(1)= CCTK_PointerTo(red_tmp)
      input_array(1)    = CCTK_PointerTo(interror1)
      call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .   1, input_array,0, 0, input_array_type, 1,
     .   input_array_type, reduction_value)
     
      if (ierror.ne.0) then
         call CCTK_WARN(1,"Reduction of norm failed!")
      end if
      interror1 = red_tmp
      
      
      input_array_type(1) = CCTK_VARIABLE_INT
      reduction_value(1)= CCTK_PointerTo(red_tmp)
      input_array(1)    = CCTK_PointerTo(interror2)
      call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .   1, input_array,0, 0, input_array_type, 1,
     .   input_array_type, reduction_value)
     
      if (ierror.ne.0) then
         call CCTK_WARN(1,"Reduction of norm failed!")
      end if
      interror2 = red_tmp

!     If there was an error on any processor then assign
!     large values to the integrals and return from the
!     subroutine (but remember to deallocate the arrays
!     first!).

      interror = .false.

      if ((interror1.ne.0).or.(interror2.ne.0)) then
         intexp  = 1.0D10
         intexp2 = 1.0D10
         intarea = 1.0D10
         intexpdel2 = 1.0D10
         interror = .true.
         goto 10
      end if


!     ***********************
!     ***   INTERPOLATE   ***
!     ***********************

!     Number of points to interpolate.

      npoints = (l_ntheta+1)*(l_nphi+1)

!     parameter, local interpolator, and coordinate system handle.
      param_table_handle = -1
      interp_handle = -1
      coord_system_handle = -1

      options_string = "order = " // char(ichar('0') + interpolation_order)
      call Util_TableCreateFromString (param_table_handle, options_string)
      if (param_table_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
      endif

      call CCTK_FortranString (nchars, interpolation_operator, operator)
      call CCTK_InterpHandle (interp_handle, operator)
      if (interp_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
      endif

      call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
      if (coord_system_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
      endif

!     fill in the input/output arrays for the interpolator
      interp_coords(1) = CCTK_PointerTo(xa)
      interp_coords(2) = CCTK_PointerTo(ya)
      interp_coords(3) = CCTK_PointerTo(za)

      call CCTK_VarIndex (vindex, "admbase::gxx")
      in_array_indices(1) = vindex
      call CCTK_VarIndex (vindex, "admbase::gyy")
      in_array_indices(2) = vindex
      call CCTK_VarIndex (vindex, "admbase::gzz")
      in_array_indices(3) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxy")
      in_array_indices(4) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxz")
      in_array_indices(5) = vindex
      call CCTK_VarIndex (vindex, "admbase::gyz")
      in_array_indices(6) = vindex
      call CCTK_VarIndex (vindex, "ahfinder::ahf_exp")
      in_array_indices(7) = vindex
      call CCTK_VarIndex (vindex, "ahfinder::ahmask")
      in_array_indices(8) = vindex
      call CCTK_VarIndex (vindex, "ahfinder::ahfgradn")
      in_array_indices(9) = vindex

      out_arrays(1) = CCTK_PointerTo(txx)
      out_arrays(2) = CCTK_PointerTo(tyy)
      out_arrays(3) = CCTK_PointerTo(tzz)
      out_arrays(4) = CCTK_PointerTo(txy)
      out_arrays(5) = CCTK_PointerTo(txz)
      out_arrays(6) = CCTK_PointerTo(tyz)
      out_arrays(7) = CCTK_PointerTo(exp)
      out_arrays(8) = CCTK_PointerTo(intmask)
      out_arrays(9) = CCTK_PointerTo(gradn)

      out_array_type_codes = CCTK_VARIABLE_REAL


!     For minimization we need the interpolated metric,
!     the expansion and the mask.

      if (.not.flow) then

         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               8, in_array_indices,
     .                               8, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif


!     For N flow, we need the interpolated metric, the expansion,
!     the norm of the gradient of the horizon function, and the mask.

      else if (nw.ne.zero) then

         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               9, in_array_indices,
     .                               9, out_array_type_codes, out_arrays)
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

!     For H or C flows, we need the interpolated expansion, the
!     norm of the gradient of the horizon function, and the mask.

      else

         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices(7),
     .                               3, out_array_type_codes(7), out_arrays(7))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

      end if

!     release parameter table
      call Util_TableDestroy (ierror, param_table_handle)

!     ***************************************
!     ***   CHECK IF WE ARE INSIDE MASK   ***
!     ***************************************

!     Check if we are either inside mask, or to close for comfort.
!     The way I do this is very simple, I just check the interpolated
!     values of the mask.  If any of these values is not 1.0, it means
!     that a point inside the mask is contaminating the interpolation,
!     and this means in turn that we are too close to the mask.

!     Notice that in practice I do not ask for the interpolated values
!     to be exactly 1.0, since round off error can modify the value
!     slightly without really meaning that we are close to the mask.

      if (.not.CCTK_EQUALS(ahf_mask,'off')) then

         if (myproc.lt.npt*npp) then
            do j=1,l_nphi+1
               do i=1,l_ntheta+1
                  if (intmask(i,j).lt.0.99D0) then
                     interror3 = 1
                  end if
               end do
            end do
         end if

         input_array_type(1) = CCTK_VARIABLE_INT
         reduction_value(1)= CCTK_PointerTo(red_tmp)
         input_array(1)    = CCTK_PointerTo(interror3)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)
     
         interror3 = red_tmp

         if (interror3.ne.0) then
            intexp  = 1.0D10
            intexp2 = 1.0D10
            intarea = 1.0D10
            intexpdel2 = 1.0D10
            interror = .true.
            goto 10
         end if

      end if


!     *************************
!     ***   AREA ELEMENTS   ***
!     *************************

!     Find 2D array for area elements.  Notice that
!     this is not needed for the flow algorithm!

      if (.not.flow) then

         if (myproc.ge.npt*npp) then

            da = zero

         else

            do j=1,l_nphi+1

!              Find sines and cosines of phi.

               cosp = cosphi(j)
               sinp = sinphi(j)

               do i=1,l_ntheta+1

!                 Find {rp}.

                  rp = rr(i,j)

!                 Find sines and cosines of theta.

                  cost = costheta(i)
                  sint = sintheta(i)

!                 Find spherical metric components.

                  trr = sint**2*(txx(i,j)*cosp**2
     .                + tyy(i,j)*sinp**2) + tzz(i,j)*cost**2
     .                + two*sint*(txy(i,j)*sint*cosp*sinp
     .                + cost*(txz(i,j)*cosp + tyz(i,j)*sinp))

                  ttt = rp**2*(cost**2*(txx(i,j)*cosp**2
     .                + tyy(i,j)*sinp**2) + tzz(i,j)*sint**2
     .                + two*cost*(txy(i,j)*cost*cosp*sinp
     .                - sint*(txz(i,j)*cosp + tyz(i,j)*sinp)))

                  tpp = (rp*sint)**2*(txx(i,j)*sinp**2 + tyy(i,j)*cosp**2
     .                - two*txy(i,j)*cosp*sinp)

                  trt = rp*(sint*cost*(txx(i,j)*cosp**2 + tyy(i,j)*sinp**2
     .                - tzz(i,j) + two*txy(i,j)*sinp*cosp)
     .                + (cost**2-sint**2)*(txz(i,j)*cosp + tyz(i,j)*sinp))

                  trp = rp*sint*(sint*(cosp*sinp*(tyy(i,j) - txx(i,j))
     .                + txy(i,j)*(cosp**2 - sinp**2))
     .                + cost*(cosp*tyz(i,j) - sinp*txz(i,j)))

                  ttp = rp**2*sint*(cost*(sinp*cosp*(tyy(i,j) - txx(i,j))
     .                + (cosp**2 - sinp**2)*txy(i,j))
     .                + sint*(sinp*txz(i,j) - cosp*tyz(i,j)))

!                 Find derivatives  {ft,fp}.  For interior points I
!                 use centered differences, and for the boundary points
!                 I use second order one sided differences. Notice that
!                 since this is done even in inter-processor boundaries,
!                 the final result will depend on the number of processors,
!                 but the differences will converge away at high order.

                  if ((i.ne.1).and.(i.ne.l_ntheta+1)) then
                     ft = idtheta*(rr(i+1,j) - rr(i-1,j))
                  else if (i.eq.1) then
                     ft = - idtheta*(three*rr(1,j)
     .                  - four*rr(2,j) + rr(3,j))
                  else
                     ft = + idtheta*(three*rr(l_ntheta+1,j)
     .                  - four*rr(l_ntheta,j) + rr(l_ntheta-1,j))
                  end if

                  if (nonaxi) then
                     if ((j.ne.1).and.(j.ne.l_nphi+1)) then
                        fp = idphi*(rr(i,j+1) - rr(i,j-1))
                     else if (j.eq.1) then
                        fp = - idphi*(three*rr(i,1)
     .                     - four*rr(i,2) + rr(i,3))
                     else
                        fp = + idphi*(three*rr(i,l_nphi+1)
     .                     - four*rr(i,l_nphi) + rr(i,l_nphi-1))
                     end if
                  else
                     fp = zero
                  end if

!                 Find area element.  For this I use the fact that if
!                 we have  r = f(theta,phi), then the area element is:
!
!                        /                  2
!                 da  =  | [ gtt + grr (f,t)  +  2 grt f,t ]
!                        \
!                                      2
!                     [ gpp + grr (f,p)  +  2 grp f,p ]  -  [ gtp
!
!                                                          2 \ 1/2
!                     + grr f,t f,p  +  grt f,p + grp f,t ]  |     dtheta dphi
!                                                            /
!
!                 where  f,t = df/dtheta  and  f,p = df/dphi.

                  aux = (ttt + trr*ft**2 + two*trt*ft)
     .                * (tpp + trr*fp**2 + two*trp*fp)
     .                - (ttp + ft*(trr*fp + trp) + trt*fp)**2

                  if (aux.gt.zero) then
                     aux = sqrt(aux)
                  else
                     aux = zero
                  end if

!                 Multiply with dtheta*dphi.

                  da(i,j) = aux*dtp

               end do
            end do

         end if

      end if


!     *******************************************
!     ***   AREA AND INTEGRALS OF EXPANSION   ***
!     *******************************************

!     To find the integrals I use a second order method:
!     I approximate the value of the integrand at each
!     cell centre by an average of the values at the
!     four cell corners, and then add up all cells.

!     Notice that for the flow algorithm I do not need to
!     calculate the area nor the integrals of the expansion
!     except at the very end, but by then the flag 'flow' has
!     been switched off.

      inside_min_count     = 0
      inside_min_neg_count = 0

      if (.not.flow) then

         if (myproc.lt.npt*npp) then

            do j=1,l_nphi
               do i=1,l_ntheta

                  intarea = intarea + quarter
     .                    *(da(i,j  ) + da(i+1,j  )
     .                    + da(i,j+1) + da(i+1,j+1))

                  intexp  = intexp  + quarter
     .                    *(da(i  ,j  )*exp(i  ,j  )
     .                    + da(i+1,j  )*exp(i+1,j  )
     .                    + da(i  ,j+1)*exp(i  ,j+1)
     .                    + da(i+1,j+1)*exp(i+1,j+1))
                  intexp2 = intexp2 + quarter
     .                    *(da(i  ,j  )*exp(i  ,j  )**2
     .                    + da(i+1,j  )*exp(i+1,j  )**2
     .                    + da(i  ,j+1)*exp(i  ,j+1)**2
     .                    + da(i+1,j+1)*exp(i+1,j+1)**2)

                  intexpdel2 = intexpdel2 + quarter
     .                 *(da(i  ,j  )*(exp(i  ,j  )+trapped_surface_delta)**2
     .                 + da(i+1,j  )*(exp(i+1,j  )+trapped_surface_delta)**2
     .                 + da(i  ,j+1)*(exp(i  ,j+1)+trapped_surface_delta)**2
     .                 + da(i+1,j+1)*(exp(i+1,j+1)+trapped_surface_delta)**2)

                  inside_min_count = inside_min_count + 1

                  if ((exp(i  ,j  ).lt.0.0D0).and.
     .                (exp(i+1,j  ).lt.0.0D0).and.
     .                (exp(i  ,j+1).lt.0.0D0).and.
     .                (exp(i+1,j+1).lt.0.0D0)) then
                     inside_min_neg_count = inside_min_neg_count + 1
                  end if

               end do
            end do

         end if

!        Reduce integrals.

         input_array_type(1) = CCTK_VARIABLE_REAL
         reduction_value(1)= CCTK_PointerTo(aux)
         input_array(1)    = CCTK_PointerTo(intarea)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)
     
	 intarea=aux

         input_array_type(1) = CCTK_VARIABLE_REAL
         reduction_value(1)= CCTK_PointerTo(aux)
         input_array(1)    = CCTK_PointerTo(intexp)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)
     
         intexp=aux

         input_array_type(1) = CCTK_VARIABLE_REAL
         reduction_value(1)= CCTK_PointerTo(aux)
         input_array(1)    = CCTK_PointerTo(intexp2)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)
     
         intexp2=aux
         
         input_array_type(1) = CCTK_VARIABLE_REAL
         reduction_value(1)= CCTK_PointerTo(aux)
         input_array(1)    = CCTK_PointerTo(intexpdel2)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)

         intexpdel2=aux

         

         input_array_type(1) = CCTK_VARIABLE_INT
         reduction_value(1)= CCTK_PointerTo(red_tmp)
         input_array(1)    = CCTK_PointerTo(inside_min_neg_count)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)
     
         inside_min_neg_count=red_tmp
         	 
         input_array_type(1) = CCTK_VARIABLE_INT
         reduction_value(1)= CCTK_PointerTo(red_tmp)
         input_array(1)    = CCTK_PointerTo(inside_min_count)
         call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .      1, input_array,0, 0, input_array_type, 1,
     .      input_array_type, reduction_value)
     
         inside_min_count=red_tmp

      end if


!     *******************************
!     ***   SPECTRAL COMPONENTS   ***
!     *******************************

!     For flow algorithm, find spectral components.
!     For this I need to recalculate the Legendre polynomials,
!     but I see no way around it.  Also, to avoid having to
!     define lots of 2D arrays I do the integrals inside
!     the loop.  This means that I have to take care to see
!     if I am on an edge or corner point. Notice also that
!     these integrals do not involve the area element.

      if (flow) then

         if (myproc.lt.npt*npp) then

!           Find h0.

            if (find_trapped_surface) then
               h0 = - trapped_surface_delta
            else
               h0 = zero
            end if

!           Find integrals.

            do j=1,l_nphi+1

!              Find phi.

               phi = dphi*dble(j+phi0) + phistart

               do i=1,l_ntheta+1

!                 Find sines and cosines of theta.

                  cost = costheta(i)
                  sint = sintheta(i)

!                 Find norm of gradient of horizon function.

                  grad = gradn(i,j)

!                 Find weight factor.

                  intw = dtp*(exp(i,j)-h0)*sint**2

!                 Modify weight factor at corners and edges.

                  if (((j.eq.1).or.(j.eq.l_nphi+1)).and.
     .               ((i.eq.1).or.(i.eq.l_ntheta+1))) then
                     intw = quarter*intw
                  else if (((j.eq.1).or.(j.eq.l_nphi+1)).or.
     .               ((i.eq.1).or.(i.eq.l_ntheta+1))) then
                     intw = half*intw
                  end if

!                 Find "sigma" for Nakamura flow (see Carstens paper).

                  if (nw.eq.0.0) then

                     sigma = zero

                  else

                     xw = xa(i,j) - xc
                     yw = ya(i,j) - yc
                     zw = za(i,j) - zc

!                    Calculate metric on surface.

                     det = txx(i,j)*tyy(i,j)*tzz(i,j)
     .                   + two*txy(i,j)*txz(i,j)*tyz(i,j)
     .                   - txx(i,j)*tyz(i,j)**2 - tyy(i,j)*txz(i,j)**2
     .                   - tzz(i,j)*txy(i,j)**2

                     idet = one/det

                     gupij(1,1) = idet*(tyy(i,j)*tzz(i,j)-tyz(i,j)**2)
                     gupij(2,2) = idet*(txx(i,j)*tzz(i,j)-txz(i,j)**2)
                     gupij(3,3) = idet*(txx(i,j)*tyy(i,j)-txy(i,j)**2)

                     gupij(1,2) = idet*(txz(i,j)*tyz(i,j)-tzz(i,j)*txy(i,j))
                     gupij(2,1) = gupij(1,2)
                     gupij(1,3) = idet*(txy(i,j)*tyz(i,j)-tyy(i,j)*txz(i,j))
                     gupij(3,1) = gupij(1,3)
                     gupij(2,3) = idet*(txy(i,j)*txz(i,j)-txx(i,j)*tyz(i,j))
                     gupij(3,2) = gupij(2,3)

                     call AHFinder_calcsigma(CCTK_ARGUMENTS,
     .               xw,yw,zw,gupij,sigma)

                  end if

!                 Monopole components.

                  aux = intw*grad

                  hflow0(0) = hflow0(0) + intw
                  cflow0(0) = cflow0(0) + aux
                  nflow0(0) = nflow0(0) + aux*sigma

!                 Axisymmetric components.

                  do l=1+stepz,lmax,1+stepz

                     aux  = intw*LEGEN(l,0,cost)
                     aux1 = aux*grad

                     hflow0(l) = hflow0(l) + aux
                     cflow0(l) = cflow0(l) + aux1
                     nflow0(l) = nflow0(l) + aux1*sigma

                  end do

!                 Non-axisymmetric components.

                  if (nonaxi) then

                     do m=1,lmax

                        aux = dble(m)*phi
                        sina = sin(aux)
                        cosa = cos(aux)

                        do l=m,lmax

                           aux  = intw*LEGEN(l,m,cost)

                           aux1 = aux*cosa
                           aux2 = aux1*grad
                           hflowc(l,m) = hflowc(l,m) + aux1
                           cflowc(l,m) = cflowc(l,m) + aux2
                           nflowc(l,m) = nflowc(l,m) + aux2*sigma

                           if (.not.refy) then
                              aux1 = aux*sina
                              aux2 = aux1*grad
                              hflows(l,m) = hflows(l,m) + aux1
                              cflows(l,m) = cflows(l,m) + aux2
                              nflows(l,m) = nflows(l,m) + aux2*sigma
                           end if

                        end do
                     end do

                  end if

               end do
            end do

         end if

!        Reduce integrals.

!        Axisymmetric components.

         auxi = lmax+1

         if (hw.ne.zero) then
            input_array_type(1) = CCTK_VARIABLE_REAL
            reduction_value(1)= CCTK_PointerTo(tempv)
            input_array(1)    = CCTK_PointerTo(hflow0)
            call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .         1, input_array,1, auxi, input_array_type, 1,
     .         input_array_type, reduction_value)

            hflow0 = tempv
	 
         end if

         if (cw.ne.zero) then
            input_array_type(1) = CCTK_VARIABLE_REAL
            reduction_value(1)= CCTK_PointerTo(tempv)
            input_array(1)    = CCTK_PointerTo(cflow0)
            call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .         1, input_array,1, auxi, input_array_type, 1,
     .         input_array_type, reduction_value)
     
            cflow0 = tempv

         end if

         if (nw.ne.zero) then
            input_array_type(1) = CCTK_VARIABLE_REAL
            reduction_value(1)= CCTK_PointerTo(tempv)
            input_array(1)    = CCTK_PointerTo(nflow0)
            call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .         1, input_array,1, auxi, input_array_type, 1,
     .         input_array_type, reduction_value)
            nflow0 = tempv

         end if

!        Non-axisymmetric components.

         if (nonaxi) then

            auxi = lmax*lmax

            if (hw.ne.zero) then
               input_array_type(1) = CCTK_VARIABLE_REAL
               reduction_value(1)= CCTK_PointerTo(tempm)
               input_array(1)    = CCTK_PointerTo(hflowc)
               call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .            1, input_array,1, auxi, input_array_type, 1,
     .            input_array_type, reduction_value)

               hflowc = tempm
            end if

            if (cw.ne.zero) then
               input_array_type(1) = CCTK_VARIABLE_REAL
               reduction_value(1)= CCTK_PointerTo(tempm)
               input_array(1)    = CCTK_PointerTo(cflowc)
               call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .            1, input_array,1, auxi, input_array_type, 1,
     .            input_array_type, reduction_value)

               cflowc = tempm
            end if

            if (nw.ne.zero) then
               input_array_type(1) = CCTK_VARIABLE_REAL
               reduction_value(1)= CCTK_PointerTo(tempm)
               input_array(1)    = CCTK_PointerTo(nflowc)
               call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .            1, input_array,1, auxi, input_array_type, 1,
     .            input_array_type, reduction_value)

               nflowc = tempm
            end if

            if (.not.refy) then

               if (hw.ne.zero) then
                  input_array_type(1) = CCTK_VARIABLE_REAL
                  reduction_value(1)= CCTK_PointerTo(tempm)
                  input_array(1)    = CCTK_PointerTo(hflows)
                  call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .               1, input_array,1, auxi, input_array_type, 1,
     .               input_array_type, reduction_value)

                  hflows = tempm
               end if

               if (cw.ne.zero) then
                  input_array_type(1) = CCTK_VARIABLE_REAL
                  reduction_value(1)= CCTK_PointerTo(tempm)
                  input_array(1)    = CCTK_PointerTo(cflows)
                  call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .               1, input_array,1, auxi, input_array_type, 1,
     .               input_array_type, reduction_value)

                  cflows = tempm
               end if

               if (nw.ne.zero) then
                  input_array_type(1) = CCTK_VARIABLE_REAL
                  reduction_value(1)= CCTK_PointerTo(tempm)
                  input_array(1)    = CCTK_PointerTo(nflows)
                  call CCTK_ReduceArraysGlobally(ierror, cctkGH, -1,reduction_handle, -1,
     .               1, input_array,1, auxi, input_array_type, 1,
     .               input_array_type, reduction_value)

                  nflows = tempm
               end if

            end if

         end if

      end if


!     *****************************************************
!     ***   RESCALE INTEGRALS ACCORDING TO SYMMETRIES   ***
!     *****************************************************

!     For cartoon multiply the integrals with 2*pi.

      if (cartoon) then

         intw = 6.283185307D0

         if (.not.flow) then
            intexp  = intw*intexp
            intexp2 = intw*intexp2
            intarea = intw*intarea
            intexpdel2 = intw*intexpdel2
         end if

         if (flow) then
            hflow0 = intw*hflow0
            cflow0 = intw*cflow0
            nflow0 = intw*nflow0
         end if

!     Other symmetries.

      else if (refx.or.refy.or.refz) then

         intw = one

         if (refx) intw = two*intw
         if (refy) intw = two*intw
         if (refz) intw = two*intw

         if (.not.flow) then
            intexp  = intw*intexp
            intexp2 = intw*intexp2
            intarea = intw*intarea
            intexpdel2 = intw*intexpdel2
         end if

!        Spectral components for flow.

         if (flow) then

!           Axisymmetric.

            hflow0 = intw*hflow0
            cflow0 = intw*cflow0
            nflow0 = intw*nflow0

!           Non-axisymmetric.

            if (nonaxi) then

               hflowc = intw*hflowc
               cflowc = intw*cflowc
               nflowc = intw*nflowc

               if (.not.refy) then
                  hflows = intw*hflows
                  cflows = intw*cflows
                  nflows = intw*nflows
               end if

            end if

         end if

      end if


!     *****************************
!     ***   DEALLOCATE MEMORY   ***
!     *****************************

   10 continue

      deallocate(costheta,sintheta)
      deallocate(cosphi,sinphi)
      deallocate(rr,xa,ya,za)
      deallocate(da,exp,gradn)
      deallocate(txx,tyy,tzz,txy,txz,tyz)
      deallocate(intmask)


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

      end subroutine AHFinder_int