minor cleaning up, plus a stupid fix to a problem with the absoft compiler.

It turns out I needed to change the order of two definitions to stop the
damm compiler from inventing a non-existent NaN.  &%$#*!


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@316 89daf98e-ef62-4674-b946-b8ff9de2216c

3 files changed, 202 insertions, 142 deletions

diff --git a/src/AHFinder_dat.F b/src/AHFinder_dat.F
index 4e31ed8..a7fc589 100644
--- a/src/AHFinder_dat.F
+++ b/src/AHFinder_dat.F
@@ -49,6 +49,8 @@
       CCTK_REAL hw,cw,nw
       CCTK_REAL dx,dy,dz
 
+      CCTK_REAL avgx,avgy,avgz
+
       CCTK_REAL, allocatable, dimension (:)   :: c0,c0_0,c0_1,c0_2,c0_old
       CCTK_REAL, allocatable, dimension (:,:) :: cc,cc_0,cc_1,cc_2,cc_old
       CCTK_REAL, allocatable, dimension (:,:) :: cs,cs_0,cs_1,cs_2,cs_old
@@ -59,8 +61,6 @@
 
       character*200 filestr
 
-      CCTK_REAL avgx, avgy
-
       data ahf_ncall / 0 /
       data status_old   / .false. /
       data status_old_0 / .false. /
@@ -163,9 +163,11 @@
 !
 !     inside_min_neg_count    Number of elements in integral that have
 !                             negative expansion.
+!
+!     avgx        X-position of centroid of apparent horizon
+!     avgy        Y-position of centroid of apparent horizon
+!     avgz        Z-position of centroid of apparent horizon
 
-!     avgx                 X-position of centroid of apparent horizon
-!     avgy                 Y-position of centroid of apparent horizon
 
 !     ***************
 !     ***   END   ***
diff --git a/src/AHFinder_exp.F b/src/AHFinder_exp.F
index e990501..388b30f 100644
--- a/src/AHFinder_exp.F
+++ b/src/AHFinder_exp.F
@@ -48,7 +48,8 @@
       CCTK_REAL c2fxx,c2fyy,c2fzz,c2fxy,c2fxz,c2fyz
       CCTK_REAL dfx,dfy,dfz,dfux,dfuy,dfuz
       CCTK_REAL aux,T0,T1,T2,T3,T4
-      CCTK_REAL idx,idy,idz,idx2,idy2,idz2,idxy,idxz,idyz
+      CCTK_REAL idx,idy,idz
+      CCTK_REAL idx2,idy2,idz2,idxy,idxz,idyz
       CCTK_REAL zero,half,one,two
 
 !     Description of variables:
@@ -105,14 +106,14 @@
       idy = half/dy
       idz = half/dz
 
-      idxy = idx*idy
-      idxz = idx*idz
-      idyz = idy*idz
-
       idx2 = one/dx**2
       idy2 = one/dy**2
       idz2 = one/dz**2
 
+      idxy = idx*idy
+      idxz = idx*idz
+      idyz = idy*idz
+
       do k=2,nz-1
          do j=2,ny-1
             do i=2,nx-1
@@ -163,11 +164,13 @@
                det = gdxx*gdyy*gdzz + two*gdxy*gdxz*gdyz
      .             - gdxx*gdyz**2 - gdyy*gdxz**2 - gdzz*gdxy**2
 
+!              If determinant is not zero proceed.
+
                if (det.gt.zero) then
 
                   idet = one/det
 
-!              Find inverse spatial metric.
+!                 Find inverse spatial metric.
 
                   guxx = idet*(gdyy*gdzz - gdyz**2)
                   guyy = idet*(gdxx*gdzz - gdxz**2)
@@ -177,7 +180,7 @@
                   guxz = idet*(gdxy*gdyz - gdyy*gdxz)
                   guyz = idet*(gdxy*gdxz - gdxx*gdyz)
 
-!              Find spatial derivatives of f.
+!                 Find spatial derivatives of f.
                   
                   T0 = two*ahfgrid(i,j,k)
               
@@ -199,13 +202,13 @@
                   dfz   = (T1 - T2)*idz
                   d2fzz = (T1 - T0 + T2)*idz2
 
-!              Save gradient of horizon function and its norm
-!              (they will be needed later in the flow algorithm
-!              and in the gaussian curvature).  Here I also try
-!              to avoid possible division by zero later by not
-!              allowing   ahfgradn(i,j,k) = 0.  This should only
-!              ever happen far from the horizon, so resetting this
-!              to 1 should have no important effects.
+!                 Save gradient of horizon function and its norm
+!                 (they will be needed later in the flow algorithm
+!                 and in the gaussian curvature).  Here I also try
+!                 to avoid possible division by zero later by not
+!                 allowing   ahfgradn(i,j,k) = 0.  This should only
+!                 ever happen far from the horizon, so resetting this
+!                 to 1 should have no important effects.
 
                   ahfgradx(i,j,k) = dfx
                   ahfgrady(i,j,k) = dfy
@@ -217,22 +220,22 @@
                   ahfgradn(i,j,k) = sqrt(aux)
                   if (ahfgradn(i,j,k).eq.zero) ahfgradn(i,j,k) = one
 
-!              Find crossed derivatives.
+!                 Find crossed derivatives.
 
                   d2fxy = (ahfgrid(i+1,j+1,k) + ahfgrid(i-1,j-1,k)
-     .                 -  ahfgrid(i+1,j-1,k) - ahfgrid(i-1,j+1,k))*idxy
+     .                  -  ahfgrid(i+1,j-1,k) - ahfgrid(i-1,j+1,k))*idxy
                   d2fxz = (ahfgrid(i+1,j,k+1) + ahfgrid(i-1,j,k-1)
-     .                 -  ahfgrid(i+1,j,k-1) - ahfgrid(i-1,j,k+1))*idxz
+     .                  -  ahfgrid(i+1,j,k-1) - ahfgrid(i-1,j,k+1))*idxz
                   d2fyz = (ahfgrid(i,j+1,k+1) + ahfgrid(i,j-1,k-1)
-     .                 -  ahfgrid(i,j+1,k-1) - ahfgrid(i,j-1,k+1))*idyz
-                  
-!              Raise indices in first derivatives.
+     .                  -  ahfgrid(i,j+1,k-1) - ahfgrid(i,j-1,k+1))*idyz
+
+!                 Raise indices in first derivatives.
 
                   dfux = guxx*dfx + guxy*dfy + guxz*dfz
                   dfuy = guxy*dfx + guyy*dfy + guyz*dfz
                   dfuz = guxz*dfx + guyz*dfy + guzz*dfz
                   
-!              Find second covariant derivatives of f.
+!                 Find second covariant derivatives of f.
 
                   c2fxx = d2fxx - half*(dfux*ddxxx
      .                 + dfuy*(two*ddxxy - ddyxx)
@@ -255,35 +258,38 @@
                   c2fyz = d2fyz - half*(dfuy*ddzyy + dfuz*ddyzz
      .                 + dfux*(ddyxz + ddzxy - ddxyz))
 
-!                           /  m      \ 1/2
-!              Find:  u  =  | d f d f |
-!                           \      m  /
+!                              /  m      \ 1/2
+!                 Find:  u  =  | d f d f |
+!                              \      m  /
 
                   T0 = dfx*dfux + dfy*dfuy + dfz*dfuz
 
+!                 if T0 is positive proceed.
+
                   if (T0.gt.zero) then
+
                      T0 = one/sqrt(T0)
 
-!                     __2
-!              Find:  \/ f / u
+!                           __2
+!                    Find:  \/ f / u
 
                      T1 = guxx*c2fxx + guyy*c2fyy + guzz*c2fzz
      .                    + two*(guxy*c2fxy + guxz*c2fxz + guyz*c2fyz)
 
                      T1 = T1*T0
 
-!                          a   b       2
-!              Find:  K   d f d f  /  u
-!                      ab
+!                                a   b       2
+!                    Find:  K   d f d f  /  u
+!                            ab
                      
                      T2 = kdxx*dfux**2 + kdyy*dfuy**2 + kdzz*dfuz**2
      .                    + two*(dfux*(kdxy*dfuy + kdxz*dfuz) + kdyz*dfuy*dfuz)
 
                      T2 = T2*T0**2
 
-!                      __a__b                3
-!              Find:   \/ \/ f  d f d f  /  u
-!                                a   b
+!                            __a__b                3
+!                    Find:   \/ \/ f  d f d f  /  u
+!                                      a   b
 
                      T3 = c2fxx*dfux**2 + c2fyy*dfuy**2 + c2fzz*dfuz**2
      .                    + two*(dfux*(c2fxy*dfuy + c2fxz*dfuz)
@@ -291,132 +297,169 @@
 
                      T3 = T3*T0**3
 
-!              Find:  trK
+!                    Find:  trK
 
                      T4 = guxx*kdxx + guyy*kdyy + guzz*kdzz
      .                    + two*(guxy*kdxy + guxz*kdxz + guyz*kdyz)
 
-!              Find the expansion.
+!                    Find the expansion.
 
                      ahf_exp(i,j,k) = T1 + T2 - T3 - T4
 
+!                 T0 was not positive.
+
                   else
                      ahf_exp(i,j,k) = zero
-                  endif
+                  end if
+
+!              Determinant was zero.
+
                else
                   ahf_exp(i,j,k) = zero
-               endif
+               end if
+
             end do
          end do
       end do
 
-!     Zero out the edges before doing the boundaries
-      ahf_exp( 1, 1, :) = 0
-      ahf_exp( 1,ny, :) = 0
-      ahf_exp( 1, :, 1) = 0
-      ahf_exp( 1, :,nz) = 0
-      ahf_exp(nx, 1, :) = 0
-      ahf_exp(nx,ny, :) = 0
-      ahf_exp(nx, :, 1) = 0
-      ahf_exp(nx, :,nz) = 0
-      ahf_exp( :, 1, 1) = 0
-      ahf_exp( :,ny, 1) = 0
-      ahf_exp( :, 1,nz) = 0
-      ahf_exp( :,ny,nz) = 0
-
-      ahfgradx( 1, 1, :) = 0
-      ahfgradx( 1,ny, :) = 0
-      ahfgradx( 1, :, 1) = 0
-      ahfgradx( 1, :,nz) = 0
-      ahfgradx(nx, 1, :) = 0
-      ahfgradx(nx,ny, :) = 0
-      ahfgradx(nx, :, 1) = 0
-      ahfgradx(nx, :,nz) = 0
-      ahfgradx( :, 1, 1) = 0
-      ahfgradx( :,ny, 1) = 0
-      ahfgradx( :, 1,nz) = 0
-      ahfgradx( :,ny,nz) = 0
-
-      ahfgrady( 1, 1, :) = 0
-      ahfgrady( 1,ny, :) = 0
-      ahfgrady( 1, :, 1) = 0
-      ahfgrady( 1, :,nz) = 0
-      ahfgrady(nx, 1, :) = 0
-      ahfgrady(nx,ny, :) = 0
-      ahfgrady(nx, :, 1) = 0
-      ahfgrady(nx, :,nz) = 0
-      ahfgrady( :, 1, 1) = 0
-      ahfgrady( :,ny, 1) = 0
-      ahfgrady( :, 1,nz) = 0
-      ahfgrady( :,ny,nz) = 0
-
-      ahfgradz( 1, 1, :) = 0
-      ahfgradz( 1,ny, :) = 0
-      ahfgradz( 1, :, 1) = 0
-      ahfgradz( 1, :,nz) = 0
-      ahfgradz(nx, 1, :) = 0
-      ahfgradz(nx,ny, :) = 0
-      ahfgradz(nx, :, 1) = 0
-      ahfgradz(nx, :,nz) = 0
-      ahfgradz( :, 1, 1) = 0
-      ahfgradz( :,ny, 1) = 0
-      ahfgradz( :, 1,nz) = 0
-      ahfgradz( :,ny,nz) = 0
-
-      ahfgradn( 1, 1, :) = 0
-      ahfgradn( 1,ny, :) = 0
-      ahfgradn( 1, :, 1) = 0
-      ahfgradn( 1, :,nz) = 0
-      ahfgradn(nx, 1, :) = 0
-      ahfgradn(nx,ny, :) = 0
-      ahfgradn(nx, :, 1) = 0
-      ahfgradn(nx, :,nz) = 0
-      ahfgradn( :, 1, 1) = 0
-      ahfgradn( :,ny, 1) = 0
-      ahfgradn( :, 1,nz) = 0
-      ahfgradn( :,ny,nz) = 0
+!     Zero out the edges before doing the boundaries.
+
+      ahf_exp( 1, 1, :) = 0.0d0
+      ahf_exp( 1,ny, :) = 0.0d0
+      ahf_exp( 1, :, 1) = 0.0d0
+      ahf_exp( 1, :,nz) = 0.0d0
+      ahf_exp(nx, 1, :) = 0.0d0
+      ahf_exp(nx,ny, :) = 0.0d0
+      ahf_exp(nx, :, 1) = 0.0d0
+      ahf_exp(nx, :,nz) = 0.0d0
+      ahf_exp( :, 1, 1) = 0.0d0
+      ahf_exp( :,ny, 1) = 0.0d0
+      ahf_exp( :, 1,nz) = 0.0d0
+      ahf_exp( :,ny,nz) = 0.0d0
+
+      ahfgradx( 1, 1, :) = 0.0d0
+      ahfgradx( 1,ny, :) = 0.0d0
+      ahfgradx( 1, :, 1) = 0.0d0
+      ahfgradx( 1, :,nz) = 0.0d0
+      ahfgradx(nx, 1, :) = 0.0d0
+      ahfgradx(nx,ny, :) = 0.0d0
+      ahfgradx(nx, :, 1) = 0.0d0
+      ahfgradx(nx, :,nz) = 0.0d0
+      ahfgradx( :, 1, 1) = 0.0d0
+      ahfgradx( :,ny, 1) = 0.0d0
+      ahfgradx( :, 1,nz) = 0.0d0
+      ahfgradx( :,ny,nz) = 0.0d0
+
+      ahfgrady( 1, 1, :) = 0.0d0
+      ahfgrady( 1,ny, :) = 0.0d0
+      ahfgrady( 1, :, 1) = 0.0d0
+      ahfgrady( 1, :,nz) = 0.0d0
+      ahfgrady(nx, 1, :) = 0.0d0
+      ahfgrady(nx,ny, :) = 0.0d0
+      ahfgrady(nx, :, 1) = 0.0d0
+      ahfgrady(nx, :,nz) = 0.0d0
+      ahfgrady( :, 1, 1) = 0.0d0
+      ahfgrady( :,ny, 1) = 0.0d0
+      ahfgrady( :, 1,nz) = 0.0d0
+      ahfgrady( :,ny,nz) = 0.0d0
+
+      ahfgradz( 1, 1, :) = 0.0d0
+      ahfgradz( 1,ny, :) = 0.0d0
+      ahfgradz( 1, :, 1) = 0.0d0
+      ahfgradz( 1, :,nz) = 0.0d0
+      ahfgradz(nx, 1, :) = 0.0d0
+      ahfgradz(nx,ny, :) = 0.0d0
+      ahfgradz(nx, :, 1) = 0.0d0
+      ahfgradz(nx, :,nz) = 0.0d0
+      ahfgradz( :, 1, 1) = 0.0d0
+      ahfgradz( :,ny, 1) = 0.0d0
+      ahfgradz( :, 1,nz) = 0.0d0
+      ahfgradz( :,ny,nz) = 0.0d0
+
+      ahfgradn( 1, 1, :) = 0.0d0
+      ahfgradn( 1,ny, :) = 0.0d0
+      ahfgradn( 1, :, 1) = 0.0d0
+      ahfgradn( 1, :,nz) = 0.0d0
+      ahfgradn(nx, 1, :) = 0.0d0
+      ahfgradn(nx,ny, :) = 0.0d0
+      ahfgradn(nx, :, 1) = 0.0d0
+      ahfgradn(nx, :,nz) = 0.0d0
+      ahfgradn( :, 1, 1) = 0.0d0
+      ahfgradn( :,ny, 1) = 0.0d0
+      ahfgradn( :, 1,nz) = 0.0d0
+      ahfgradn( :,ny,nz) = 0.0d0
 
 !     Boundaries on x direction.
 
-      ahf_exp(1,2:ny-1,2:nz-1)  = 2.0D0*ahf_exp(2,2:ny-1,2:nz-1) - ahf_exp(3,2:ny-1,2:nz-1)
-      ahfgradx(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(2,2:ny-1,2:nz-1) - ahfgradx(3,2:ny-1,2:nz-1)
-      ahfgrady(1,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(2,2:ny-1,2:nz-1) - ahfgrady(3,2:ny-1,2:nz-1)
-      ahfgradz(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(2,2:ny-1,2:nz-1) - ahfgradz(3,2:ny-1,2:nz-1)
-      ahfgradn(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(2,2:ny-1,2:nz-1) - ahfgradn(3,2:ny-1,2:nz-1)
-
-      ahf_exp(nx,2:ny-1,2:nz-1)  = 2.0D0*ahf_exp(nx-1,2:ny-1,2:nz-1) - ahf_exp(nx-2,2:ny-1,2:nz-1)
-      ahfgradx(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(nx-1,2:ny-1,2:nz-1) - ahfgradx(nx-2,2:ny-1,2:nz-1)
-      ahfgrady(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(nx-1,2:ny-1,2:nz-1) - ahfgrady(nx-2,2:ny-1,2:nz-1)
-      ahfgradz(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(nx-1,2:ny-1,2:nz-1) - ahfgradz(nx-2,2:ny-1,2:nz-1)
-      ahfgradn(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(nx-1,2:ny-1,2:nz-1) - ahfgradn(nx-2,2:ny-1,2:nz-1)
+      ahf_exp(1,2:ny-1,2:nz-1)  = 2.0D0*ahf_exp(2,2:ny-1,2:nz-1)
+     .   - ahf_exp(3,2:ny-1,2:nz-1)
+      ahfgradx(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(2,2:ny-1,2:nz-1)
+     .   - ahfgradx(3,2:ny-1,2:nz-1)
+      ahfgrady(1,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(2,2:ny-1,2:nz-1)
+     .   - ahfgrady(3,2:ny-1,2:nz-1)
+      ahfgradz(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(2,2:ny-1,2:nz-1)
+     .   - ahfgradz(3,2:ny-1,2:nz-1)
+      ahfgradn(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(2,2:ny-1,2:nz-1)
+     .   - ahfgradn(3,2:ny-1,2:nz-1)
+
+      ahf_exp(nx,2:ny-1,2:nz-1)  = 2.0D0*ahf_exp(nx-1,2:ny-1,2:nz-1)
+     .   - ahf_exp(nx-2,2:ny-1,2:nz-1)
+      ahfgradx(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(nx-1,2:ny-1,2:nz-1)
+     .   - ahfgradx(nx-2,2:ny-1,2:nz-1)
+      ahfgrady(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(nx-1,2:ny-1,2:nz-1)
+     .   - ahfgrady(nx-2,2:ny-1,2:nz-1)
+      ahfgradz(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(nx-1,2:ny-1,2:nz-1)
+     .   - ahfgradz(nx-2,2:ny-1,2:nz-1)
+      ahfgradn(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(nx-1,2:ny-1,2:nz-1)
+     .   - ahfgradn(nx-2,2:ny-1,2:nz-1)
 
 !     Boundaries on y direction.
 
-      ahf_exp(2:nx-1,1,2:nz-1)  = 2.0D0*ahf_exp(2:nx-1,2,2:nz-1) - ahf_exp(2:nx-1,3,2:nz-1)
-      ahfgradx(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,2,2:nz-1) - ahfgradx(2:nx-1,3,2:nz-1)
-      ahfgrady(2:nx-1,1,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,2,2:nz-1) - ahfgrady(2:nx-1,3,2:nz-1)
-      ahfgradz(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,2,2:nz-1) - ahfgradz(2:nx-1,3,2:nz-1)
-      ahfgradn(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,2,2:nz-1) - ahfgradn(2:nx-1,3,2:nz-1)
-
-      ahf_exp(2:nx-1,ny,2:nz-1)  = 2.0D0*ahf_exp(2:nx-1,ny-1,2:nz-1) - ahf_exp(2:nx-1,ny-2,2:nz-1)
-      ahfgradx(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,ny-1,2:nz-1) - ahfgradx(2:nx-1,ny-2,2:nz-1)
-      ahfgrady(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,ny-1,2:nz-1) - ahfgrady(2:nx-1,ny-2,2:nz-1)
-      ahfgradz(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,ny-1,2:nz-1) - ahfgradz(2:nx-1,ny-2,2:nz-1)
-      ahfgradn(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,ny-1,2:nz-1) - ahfgradn(2:nx-1,ny-2,2:nz-1)
+      ahf_exp(2:nx-1,1,2:nz-1)  = 2.0D0*ahf_exp(2:nx-1,2,2:nz-1)
+     .   - ahf_exp(2:nx-1,3,2:nz-1)
+      ahfgradx(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,2,2:nz-1)
+     .   - ahfgradx(2:nx-1,3,2:nz-1)
+      ahfgrady(2:nx-1,1,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,2,2:nz-1)
+     .   - ahfgrady(2:nx-1,3,2:nz-1)
+      ahfgradz(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,2,2:nz-1)
+     .   - ahfgradz(2:nx-1,3,2:nz-1)
+      ahfgradn(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,2,2:nz-1)
+     .   - ahfgradn(2:nx-1,3,2:nz-1)
+
+      ahf_exp(2:nx-1,ny,2:nz-1)  = 2.0D0*ahf_exp(2:nx-1,ny-1,2:nz-1)
+     .   - ahf_exp(2:nx-1,ny-2,2:nz-1)
+      ahfgradx(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,ny-1,2:nz-1)
+     .   - ahfgradx(2:nx-1,ny-2,2:nz-1)
+      ahfgrady(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,ny-1,2:nz-1)
+     .   - ahfgrady(2:nx-1,ny-2,2:nz-1)
+      ahfgradz(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,ny-1,2:nz-1)
+     .   - ahfgradz(2:nx-1,ny-2,2:nz-1)
+      ahfgradn(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,ny-1,2:nz-1)
+     .   - ahfgradn(2:nx-1,ny-2,2:nz-1)
 
 !     Boundaries on z direction.
 
-      ahf_exp(2:nx-1,2:ny-1,1)  = 2.0D0*ahf_exp(2:nx-1,2:ny-1,2) - ahf_exp(2:nx-1,2:ny-1,3)
-      ahfgradx(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,2) - ahfgradx(2:nx-1,2:ny-1,3)
-      ahfgrady(2:nx-1,2:ny-1,1) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,2) - ahfgrady(2:nx-1,2:ny-1,3)
-      ahfgradz(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,2) - ahfgradz(2:nx-1,2:ny-1,3)
-      ahfgradn(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,2) - ahfgradn(2:nx-1,2:ny-1,3)
-
-      ahf_exp(2:nx-1,2:ny-1,nz)  = 2.0D0*ahf_exp(2:nx-1,2:ny-1,nz-1) - ahf_exp(2:nx-1,2:ny-1,nz-2)
-      ahfgradx(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,nz-1) - ahfgradx(2:nx-1,2:ny-1,nz-2)
-      ahfgrady(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,nz-1) - ahfgrady(2:nx-1,2:ny-1,nz-2)
-      ahfgradz(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,nz-1) - ahfgradz(2:nx-1,2:ny-1,nz-2)
-      ahfgradn(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,nz-1) - ahfgradn(2:nx-1,2:ny-1,nz-2)
+      ahf_exp(2:nx-1,2:ny-1,1)  = 2.0D0*ahf_exp(2:nx-1,2:ny-1,2)
+     .   - ahf_exp(2:nx-1,2:ny-1,3)
+      ahfgradx(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,2)
+     .   - ahfgradx(2:nx-1,2:ny-1,3)
+      ahfgrady(2:nx-1,2:ny-1,1) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,2)
+     .   - ahfgrady(2:nx-1,2:ny-1,3)
+      ahfgradz(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,2)
+     .   - ahfgradz(2:nx-1,2:ny-1,3)
+      ahfgradn(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,2)
+     .   - ahfgradn(2:nx-1,2:ny-1,3)
+
+      ahf_exp(2:nx-1,2:ny-1,nz)  = 2.0D0*ahf_exp(2:nx-1,2:ny-1,nz-1)
+     .   - ahf_exp(2:nx-1,2:ny-1,nz-2)
+      ahfgradx(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,nz-1)
+     .   - ahfgradx(2:nx-1,2:ny-1,nz-2)
+      ahfgrady(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,nz-1)
+     .   - ahfgrady(2:nx-1,2:ny-1,nz-2)
+      ahfgradz(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,nz-1)
+     .   - ahfgradz(2:nx-1,2:ny-1,nz-2)
+      ahfgradn(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,nz-1)
+     .   - ahfgradn(2:nx-1,2:ny-1,nz-2)
 
 !     Synchronize.
 
diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index 05160c3..28ab392 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -505,6 +505,7 @@
 
          avgx = 0.0D0
          avgy = 0.0D0
+         avgz = 0.0D0
 
          do j=1,nphi+1
             do i=1,ntheta+1
@@ -573,9 +574,6 @@
                xp = xa(i,j)
                yp = ya(i,j)
                zp = za(i,j)
- 
-               avgx = avgx + xp
-               avgy = avgy + yp
 
                if ((xp.gt.xmx).or.(xp.lt.xmn).or.
      .             (yp.gt.ymx).or.(yp.lt.ymn).or.
@@ -583,22 +581,39 @@
                   error2 = 1
                end if
 
+!              Add to sums for average position.
+
+               avgx = avgx + xp
+               avgy = avgy + yp
+               avgz = avgz + zp
+
             end do
          end do
 
+!        Find position of horizon centroid.
+
          avgx = avgx/dble(npoints)
          avgy = avgy/dble(npoints)
+         avgz = avgz/dble(npoints)
 
          if ((.not.find3).or.
      .       (find3.and.(mfind-1.eq.horizon_centroid))) then
+
             ahf_centroid_x = avgx
             ahf_centroid_y = avgy
 
+!           Take symmetries into account.
+
+            if (refx) ahf_centroid_x = xc
+            if (refy) ahf_centroid_y = yc
+
             if (verbose) then
+              write(*,*)
               write(6,*) 'AHFinder_gau: ahf_centroid_x = ', ahf_centroid_x
               write(6,*) 'AHFinder_gau: ahf_centroid_y = ', ahf_centroid_y
             end if
-         endif
+
+         end if
 
 !     Other processors.