path: root/src/AHFinder_int.F

/*@@
  @file      AHFinder_int.F
  @date      April 1998
  @author    Miguel Alcubierre
  @desc 
             Find area of surface and integral of expansion.
             This is 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 use less processors).
  @enddesc 
@@*/

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

      subroutine AHFinder_int(CCTK_FARGUMENTS)

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

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

      CCTK_REAL LEGEN
      CCTK_REAL theta,phi,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 xmn,ymn,zmn,xmx,ymx,zmx
      CCTK_REAL dtheta,dphi,dtp,idtheta,idphi
      CCTK_REAL det,idet
      CCTK_REAL zero,quarter,half,one,two,three,four,pi
      CCTK_REAL sina,cosa,intw,grad,sigma,h0
      CCTK_REAL aux,aux1,aux2

      CCTK_REAL gupij(3,3)
      CCTK_REAL tempv(0:lmax),tempm(lmax,lmax)
      CCTK_REAL, dimension(3) :: origin,maxval,minval
      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 related things.

      CCTK_INT reduction_handle,red_tmp

!     Variables to be saved for next call.

      save xmn,ymn,zmn,xmx,ymx,zmx
      save dtheta,dphi,dtp,idtheta,idphi
      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
!
!     xmn        Minimum value of x over the grid.
!     xmx        Maximum value of x over the grid.
!
!     ymn        Minimum value of y over the grid.
!     ymx        Maximum value of y over the grid.
!
!     zmn        Minimum value of z over the grid.
!     zmx        Maximum value of z over the grid.
!
!     dtheta     Grid spacing in theta.
!     dphi       Grid spacing in phi.
!     dtp        dtheta*dphi.
!
!     idtheta    1/(2 dtheta)
!     idphi      1/(2 dphi)
!
!     da         Area element.

!     Get the reduction handle for the sum operation.

      call CCTK_ReductionArrayHandle(reduction_handle,"sum")

      if (reduction_handle.lt.0) then 
        call CCTK_WARN(1,"Cannot get reduction handle for SUM operation.")
      endif


!     **************************
!     ***   DEFINE NUMBERS   ***
!     **************************

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

      pi = 3.141592654D0

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


!     *********************************
!     ***   INITIALIZE PARAMETERS   ***
!     *********************************

      if (firstint) then

         firstint = .false.

!        Find grid boundaries.
         
         call CCTK_CoordRange(ierror,cctkGH,xmn,xmx,"x")
         call CCTK_CoordRange(ierror,cctkGH,ymn,ymx,"y")
         call CCTK_CoordRange(ierror,cctkGH,zmn,zmx,"z")

!        Find {dtheta,dphi,dtp}.

         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
               call CCTK_INFO('This combination of symmetries has not been properly implemented yet!')
            else if (refy) then
               dphi = half*dphi
            end if

            dtp = dtheta*dphi

            idtheta = half/dtheta
            idphi   = half/dphi

         end if

!        Check the number of processors to figure out
!        how to divide the domain.

!        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))

            if (nprocs.eq.1) then

!              One processor.

               npt = 1
               npp = 1

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

!              "nprocs" is the square of an integer.

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

            else

!              Check if "nprocs" is a power of two.

               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

!        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

!           Take first 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(int(myproc),int(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(rr(0:l_ntheta,0:l_nphi))

      allocate(xa(0:l_ntheta,0:l_nphi))
      allocate(ya(0:l_ntheta,0:l_nphi))
      allocate(za(0:l_ntheta,0:l_nphi))

      allocate(da(0:l_ntheta,0:l_nphi))
      allocate(exp(0:l_ntheta,0:l_nphi))
      allocate(gradn(0:l_ntheta,0:l_nphi))

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

      allocate(intmask(0:l_ntheta,0:l_nphi))

!     Initialize.

      rr = zero

      xa = zero; ya = zero; za = zero

      da = zero; exp = zero; gradn = zero

      txx = one;  tyy = one;  tzz = one
      txy = zero; txz = zero; tyz = zero


!     ***************************
!     ***   START MAIN LOOP   ***
!     ***************************

!     Initialize {interror1,interror2}.

      interror1 = 0
      interror2 = 0

!     Initialize {intexp,intexp2,intarea}.

      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=0,l_nphi
            do i=0,l_ntheta

!              Find {theta,phi}.

               theta = dtheta*dble(i+theta0)
               phi   = dphi*dble(j+phi0)

!              Find sines and cosines.

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

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

!              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 s 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

!              Save radius array.

               rr(i,j) = rp

!              Find {xa,ya,za}

               xa(i,j) = rp*sint*cosp + xc
               ya(i,j) = rp*sint*sinp + yc
               za(i,j) = rp*cost + zc

!              Check if we are within bounds.

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

               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 the errors across processors.

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .   interror1,red_tmp,CCTK_VARIABLE_INT)

      if (ierror.ne.0) then
         call CCTK_WARN(1,"Reduction of norm failed!")
      end if

      interror1 = red_tmp

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .   interror2,red_tmp,CCTK_VARIABLE_INT) 

      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 expansion and metric.

      npoints = (l_ntheta+1)*(l_nphi+1)
      origin  = cctk_origin_space + cctk_lbnd*cctk_delta_space

      call CCTK_GetInterpHandle(handle,"simple")

      call CCTK_Interp(ierror,cctkGH,handle,npoints,3,9,9,
     .   nx,ny,nz,xa,ya,za,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     .   CCTK_VARIABLE_REAL,origin(1),origin(2),origin(3),
     .   dx,dy,dz,gxx,gyy,gzz,gxy,gxz,gyz,ahf_exp,ahfgradn,ahmask,
     .   CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     .   CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     .   CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     .   txx,tyy,tzz,txy,txz,tyz,exp,gradn,intmask,
     .   CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     .   CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     .   CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL)

!     Check if we are either inside mask, or to close for comfort.

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

         interror3 = 0

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

         call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .      interror3,red_tmp,CCTK_VARIABLE_INT)
         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

!     Find 2D array for area element.

      if (myproc.ge.npt*npp) then

         da = zero

      else

         do j=0,l_nphi
            do i=0,l_ntheta

!              Find {theta,phi}.

               theta = dtheta*dble(i+theta0)
               phi   = dphi*dble(j+phi0)

!              Find {rp}.

               rp = rr(i,j)

!              Find sines and cosines.

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

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

!              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.0).and.(i.ne.l_ntheta)) then
                  ft = idtheta*(rr(i+1,j) - rr(i-1,j))
               else if (i.eq.0) then
                  ft = - idtheta*(three*rr(0,j)
     .               - four*rr(1,j) + rr(2,j))
               else
                  ft = + idtheta*(three*rr(l_ntheta,j)
     .               - four*rr(l_ntheta-1,j) + rr(l_ntheta-2,j))
               end if

               if (nonaxi) then
                  if ((j.ne.0).and.(j.ne.l_nphi)) then
                     fp = idphi*(rr(i,j+1) - rr(i,j-1))
                  else if (j.eq.0) then
                     fp = - idphi*(three*rr(i,0)
     .                  - four*rr(i,1) + rr(i,2))
                  else
                     fp = + idphi*(three*rr(i,l_nphi)
     .                  - four*rr(i,l_nphi-1) + rr(i,l_nphi-2))
                  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

!        Find integrals.  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.

         inside_min_count     = 0
         inside_min_neg_count = 0

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

               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  ) .le. 0.0d0) .and.
     .             (exp(i+1,j  ) .le. 0.0d0) .and.
     .             (exp(i  ,j+1) .le. 0.0d0) .and.
     .             (exp(i+1,j+1) .le. 0.0d0)) then
                  inside_min_neg_count = inside_min_neg_count + 1
               endif

            end do
         end do

!        If necessary, find spectral components for flow.
!        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

!           Find h0.

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

!           Find integrals.

            do j=0,l_nphi
               do i=0,l_ntheta

!                 Find {theta,phi}.

                  theta = dtheta*dble(i+theta0)
                  phi   = dphi*dble(j+phi0)

!                 Find cost,sint.

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

!                 Find norm of gradient of horizon function.

                  grad = gradn(i,j)

!                 Find weight factor.

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

                  if (((j.eq.0).or.(j.eq.l_nphi)).and.
     .               ((i.eq.0).or.(i.eq.l_ntheta))) then
                     intw = quarter*intw
                  else if (((j.eq.0).or.(j.eq.l_nphi)).or.
     .               ((i.eq.0).or.(i.eq.l_ntheta))) 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_FARGUMENTS,
     .               xw,yw,zw,gupij,sigma)

                  end if

!                 Monopole integrals.

                  aux = intw*grad

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

!                 Axisymmetric integrals.

                  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 integrals.

                  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

      end if

!     Reduce the integrals across processors.

!     Area and expansion.

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .     intarea,aux,CCTK_VARIABLE_REAL) 
      intarea=aux

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .     intexp,aux,CCTK_VARIABLE_REAL) 
      intexp=aux

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .     intexp2,aux,CCTK_VARIABLE_REAL) 
      intexp2=aux

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .     intexpdel2,aux,CCTK_VARIABLE_REAL) 
      intexpdel2=aux

!     negative expansion elements on surface

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .     inside_min_neg_count,auxi,CCTK_VARIABLE_INT) 
      inside_min_neg_count=auxi

      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
     .     inside_min_count,auxi,CCTK_VARIABLE_INT) 
      inside_min_count=auxi

!     Spectral components for flow.

      if (flow) then

!        Axisymmetric.

         auxi = lmax+1

	 if (hw.ne.zero) then
            call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .           hflow0,tempv,auxi,CCTK_VARIABLE_REAL) 
            hflow0 = tempv
	 end if

         if (cw.ne.zero) then      
            call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .          cflow0,tempv,auxi,CCTK_VARIABLE_REAL) 
            cflow0 = tempv
         end if

         if (nw.ne.zero) then
            call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .           nflow0,tempv,auxi,CCTK_VARIABLE_REAL) 
             nflow0 = tempv
         end if

!        Non-axisymmetric.

         if (nonaxi) then

            auxi = lmax*lmax

	    if (nw.ne.zero) then
               call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .              hflowc,tempm,auxi,CCTK_VARIABLE_REAL) 
               hflowc = tempm
	    end if

            if (cw.ne.zero) then
               call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .              cflowc,tempm,auxi,CCTK_VARIABLE_REAL) 
               cflowc = tempm
            end if

            if (nw.ne.zero) then
               call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .              nflowc,tempm,auxi,CCTK_VARIABLE_REAL) 
               nflowc = tempm
            end if

            if (.not.refy) then

	       if (hw.ne.zero) then
                  call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .                 hflows,tempm,auxi,CCTK_VARIABLE_REAL) 
                  hflows = tempm
	       end if

               if (cw.ne.zero) then
                  call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .                 cflows,tempm,auxi,CCTK_VARIABLE_REAL) 
                  cflows = tempm
               end if

               if (nw.ne.zero) then
                  call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
     .                 nflows,tempm,auxi,CCTK_VARIABLE_REAL) 
                  nflows = tempm
               end if

            end if

         end if

      end if

!     For cartoon multiply the integrals with 2*pi.

      if (cartoon) then

         intw = 6.283185307D0

         intexp  = intw*intexp
         intexp2 = intw*intexp2
         intarea = intw*intarea
         intexpdel2 = intw*intexpdel2

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

!     Rescale the integrals according to 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

         intexp  = intw*intexp
         intexp2 = intw*intexp2
         intarea = intw*intarea

!        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(rr,xa,ya,za)
      deallocate(da,exp,gradn)
      deallocate(txx,tyy,tzz,txy,txz,tyz)
      deallocate(intmask)


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

      end subroutine AHFinder_int