Preparing for release. Added comments, deleted obsolete pieces of code, fixed

a few bugs.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@143 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 704 insertions, 0 deletions

diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90
new file mode 100644
index 0000000..db2a315
--- /dev/null
+++ b/src/EHFinder_Integrate.F90
@@ -0,0 +1,704 @@
+! Integrating over the surface(s). Right now only in full mode.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k, l, im, jm
+  CCTK_INT :: interp_handle, table_handle, coord_system_handle
+
+  character(len=200) :: area_interp
+  character(len=128) :: warn_message
+  CCTK_INT :: area_interp_len
+  character(len=7) :: area_order
+
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
+
+  CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv
+  CCTK_REAL :: dxdth, dxdph, dydth, dydph, dzdth, dzdph
+  CCTK_POINTER, dimension(3) :: interp_coords
+  CCTK_POINTER, dimension(7) :: out_array
+  CCTK_INT, dimension(7) :: in_array
+  CCTK_INT, dimension(7), parameter :: out_types = (/ CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL, &
+                                                      CCTK_VARIABLE_REAL /)
+  
+! If finding of surface failed do not try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
+! Store the levelset_counter in a shorter named variable for convenience.
+  l = levelset_counter
+
+! Write an Info message about what we are doing.
+  call CCTK_INFO ( 'Finding surface element' )
+
+! Obtain the dimensions and characteristics of the grid arrays.
+  call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get bounding box for surface arrays' )
+  end if
+  call CCTK_GroupgshGN ( ierr, cctkGH, 2, gsh, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get global size for surface arrays' )
+  end if
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 2, ubnd, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get upper bounds for surface arrays' )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
+  end if
+  call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, &
+                                           'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
+  end if
+
+! Setup the theta and phi spacings and their inverse.
+  dtheta = ctheta(2,1) - ctheta(1,1)
+  dphi = cphi(1,2) - cphi(1,1)
+  dthetainv = one / dtheta
+  dphiinv = one / dphi
+
+! Convert the interpolator parameter into a fortran string.
+  call CCTK_FortranString ( area_interp_len, area_interpolator, &
+                                             area_interp )
+
+! Get an interpolation handle.
+  call CCTK_InterpHandle ( interp_handle, area_interp(1:area_interp_len) )
+
+  if ( interp_handle .lt. 0 ) then
+    warn_message = 'Cannot get handle for interpolation. '
+    warn_message = warn_message//'Forgot to activate an implementation '
+    warn_message = warn_message//'providing interpolation operators??'
+    call CCTK_WARN( 0, warn_message )
+  end if
+
+! Write the interpolation order parameter into a string
+  write(area_order,'(a6,i1)') 'order=',area_interpolation_order
+
+! Create a table from the interpolation order string.
+  call Util_TableCreateFromString ( table_handle, area_order )
+  if ( table_handle .lt. 0 ) then
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
+  end if
+
+! Get a handle for the coordinate system.
+  call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' )
+  if ( coord_system_handle .lt. 0) then
+    warn_message = 'Cannot get handle for cart3d coordinate system. '
+    warn_message = warn_message//'Forgot to activate an implementation '
+    warn_message = warn_message//'providing coordinates??'
+    call CCTK_WARN( 0, warn_message )
+  endif
+
+! Get the pointers to the interpolation points.
+  interp_coords(1) = CCTK_PointerTo(interp_x)
+  interp_coords(2) = CCTK_PointerTo(interp_y)
+  interp_coords(3) = CCTK_PointerTo(interp_z)
+
+! Get the pointers to the interpolation return arrays.
+  out_array(1) = CCTK_PointerTo(gxxi)
+  out_array(2) = CCTK_PointerTo(gxyi)
+  out_array(3) = CCTK_PointerTo(gxzi)
+  out_array(4) = CCTK_PointerTo(gyyi)
+  out_array(5) = CCTK_PointerTo(gyzi)
+  out_array(6) = CCTK_PointerTo(gzzi)
+  out_array(7) = CCTK_PointerTo(psii)
+
+! find the cartesian coordinates for the interpolation points
+  do j = 1, lsh(2)
+    do i = 1, lsh(1)
+      interp_x(i,j) = center(1) + rsurf(i,j) * sintheta(i,j) * cosphi(i,j)
+      interp_y(i,j) = center(2) + rsurf(i,j) * sintheta(i,j) * sinphi(i,j)
+      interp_z(i,j) = center(3) + rsurf(i,j) * costheta(i,j)
+    end do
+  end do
+
+! Get the indices for the grid functions to be interpolated.
+  call CCTK_VarIndex ( in_array(1), 'admbase::gxx' )
+  call CCTK_VarIndex ( in_array(2), 'admbase::gxy' )
+  call CCTK_VarIndex ( in_array(3), 'admbase::gxz' )
+  call CCTK_VarIndex ( in_array(4), 'admbase::gyy' )
+  call CCTK_VarIndex ( in_array(5), 'admbase::gyz' )
+  call CCTK_VarIndex ( in_array(6), 'admbase::gzz' )
+
+! If the static conformal factor is used, we also need to interpolate it.
+  if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
+    call CCTK_VarIndex ( in_array(7), 'staticconformal::psi' )
+
+!   Call the interpolator for the metric and the conformal factor.
+    call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 7, in_array, &
+                                 7, out_types, out_array )
+
+!   Get the physical metric.
+    gxxi = psii**4 * gxxi
+    gxyi = psii**4 * gxyi
+    gxzi = psii**4 * gxzi
+    gyyi = psii**4 * gyyi
+    gyzi = psii**4 * gyzi
+    gzzi = psii**4 * gzzi
+  else
+!   Else call the interpolater only for the metric.
+    call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, &
+                                 table_handle, coord_system_handle, &
+                                 lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
+                                 interp_coords, 6, in_array(1:6), &
+                                 6, out_types(1:6), out_array(1:6) )
+  end if
+
+! For all points on the surface...
+  do j = 1, lsh(2)
+    do i = 1, lsh(1)
+
+!     Calculate the coordinate transformation between the 3-d cartesian
+!     coordinates and the 2-d (theta,phi) coordinates on the surface.
+      dxdth = ( drdtheta(i,j) * sintheta(i,j) + &
+                rsurf(i,j) * costheta(i,j) ) * cosphi(i,j)
+      dydth = ( drdtheta(i,j) * sintheta(i,j) + &
+                rsurf(i,j) * costheta(i,j) ) * sinphi(i,j)
+      dzdth = drdtheta(i,j) * costheta(i,j) - rsurf(i,j) * sintheta(i,j)
+      dxdph = ( drdphi(i,j) * cosphi(i,j) - &
+                rsurf(i,j) * sinphi(i,j) ) * sintheta(i,j)
+      dydph = ( drdphi(i,j) * sinphi(i,j) + &
+                rsurf(i,j) * cosphi(i,j) ) * sintheta(i,j)
+      dzdph = drdphi(i,j) * costheta(i,j)
+
+!     Using this coordinate transformation, calculate the induced 2-metric
+!     on the surface.
+      gtt(i,j) = dxdth**2 * gxxi(i,j) + dydth**2 * gyyi(i,j) + &
+            dzdth**2 * gzzi(i,j) + &
+            two * ( dxdth * dydth * gxyi(i,j) + &
+                    dxdth * dzdth * gxzi(i,j) + &
+                    dydth * dzdth * gyzi(i,j) )
+      gtp(i,j) = dxdth * dxdph * gxxi(i,j) + &
+            ( dxdth * dydph + dydth * dxdph ) * gxyi(i,j) + &
+            ( dxdth * dzdph + dzdth * dxdph ) * gxzi(i,j) + &
+            dydth * dydph * gyyi(i,j) + &
+            ( dydth * dzdph + dzdth * dydph ) * gyzi(i,j) + &
+            dzdth * dzdph * gzzi(i,j)
+      gpp(i,j) = dxdph**2 * gxxi(i,j) + dydph**2 * gyyi(i,j) + &
+            dzdph**2 * gzzi(i,j) + &
+            two * ( dxdph * dydph * gxyi(i,j) + &
+                    dxdph * dzdph * gxzi(i,j) + &
+                    dydph * dzdph * gyzi(i,j) )
+
+!     Calculate the area element as the square-root of the determinant
+!     of the two metric multiplied by dtheta*dphi.
+      da(i,j) = dtheta * dphi * sqrt ( gtt(i,j) * gpp(i,j) - gtp(i,j)**2 )
+
+!     Calculate the arc length element for the polar circumference integration.
+      dltheta(i,j) = dtheta * sqrt ( gtt(i,j) )
+
+!     Calculate the arc length element for the equatorial circumference
+!     integration.
+      dlphi(i,j) = dphi * sqrt ( gpp(i,j) )
+
+    end do
+  end do
+ 
+end subroutine EHFinder_FindSurfaceElement
+
+
+subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: int_var, sn, l
+  CCTK_REAL :: eh_area_tmp
+  character(len=30) :: info_message
+  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
+
+! If finding of surface failed do not try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
+! Store the levelset_counter in a shorter named variable for convenience.
+  l = levelset_counter
+
+! Write an info message, to indicate what we are doing.
+  call CCTK_INFO ( 'Integrating area' )
+
+  sn = surface_counter - integrate_counter
+
+! Taking into account the sym_factor and the weights, the area is simply
+! The sum of the elements of this array.
+  int_tmp = sym_factor * weights * da
+
+! Get the index to the area integration array.
+  call CCTK_VarIndex ( int_var, 'ehfinder::int_tmp' )
+  if ( int_var .lt. 0 ) then
+    call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
+  end if
+
+! Do a sum reduce over all processors. The result is the area.
+  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+                     eh_area_tmp, 1, int_var )
+
+! Write an info message reporting the area found.
+  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp
+
+  call CCTK_INFO(info_message)
+
+! Figure out the bounds of the area storage grid array.
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, 'ehfinder::eh_area2' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get lower bounds for area array' )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, 'ehfinder::eh_area2' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get upper bounds for area array' )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, 'ehfinder::eh_area2' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for area array' )
+  end if
+
+! If this processor contains the right part of the area storage grid array
+! save the area, otherwise do nothing.
+  if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then
+    eh_area2(sn-lbnd(1),l) = eh_area_tmp
+  end if
+
+end subroutine EHFinder_IntegrateArea
+
+
+subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: int_var, sn, k, l
+  character(len=50) :: info_message
+  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
+  CCTK_REAL :: centroid_x, centroid_y, centroid_z
+
+! If finding of surface failed do not try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
+! Store the levelset_counter in a shorter named variable for convenience.
+  l = levelset_counter
+
+! Write an info message, to indicate what we are doing.
+  call CCTK_INFO ( 'Integrating centroid' )
+
+  sn = surface_counter - integrate_counter
+
+! Calculate the (x,y,z) coordinates of the points on the surface.
+  interp_x = center(1) + rsurf * sintheta * cosphi
+  interp_y = center(2) + rsurf * sintheta * sinphi
+  interp_z = center(3) + rsurf * costheta
+
+! Get the index to the temporary integration array.
+  call CCTK_VarIndex ( int_var, 'ehfinder::int_tmp' )
+  if ( int_var .lt. 0 ) then
+    call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
+  end if
+
+! Take into account the symmetry, weights, x-coordinate and area element
+! to calculate the integral of x*da over the surface. This gives 
+! centroid_x * Area, so later we divide by the area.
+  int_tmp = sym_factor * interp_x * weights * da
+
+! Sum up the all elements of int_tmp.
+  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+                     centroid_x, 1, int_var )
+
+! Take into account the symmetry, weights, y-coordinate and area element
+! to calculate the integral of y*da over the surface. This gives
+! centroid_y * Area, so later we divide by the area.
+  int_tmp = sym_factor * interp_y * weights * da
+
+! Sum up the all elements of int_tmp.
+  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+                     centroid_y, 1, int_var )
+
+! Take into account the symmetry, weights, z-coordinate and area element
+! to calculate the integral of z*da over the surface. This gives
+! centroid_z * Area, so later we divide by the area.
+  int_tmp = sym_factor * interp_z * weights * da
+
+! Sum up the all elements of int_tmp.
+  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+                     centroid_z, 1, int_var )
+! If this surface has symmetries set the corresponding centroids to zero.
+  if ( s_symx ) centroid_x = zero
+  if ( s_symy ) centroid_y = zero
+  if ( s_symz ) centroid_z = zero
+
+! Get the bounds for the centroid arrays. It should be enough to do
+! it just for the x-array, since they are all defined in the same way.
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, 'ehfinder::eh_centroid2_x' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get lower bounds for centroid array' )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, 'ehfinder::eh_centroid2_x' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get upper bounds for centroid array' )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, 'ehfinder::eh_centroid2_x' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for centroid array' )
+  end if
+
+! If we are on the right processor, store the centroid, otherwise do nothing.
+  if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then
+    k = sn - lbnd(1)
+    eh_centroid2_x(k,l) = centroid_x / eh_area2(k,l)
+    eh_centroid2_y(k,l) = centroid_y / eh_area2(k,l)
+    eh_centroid2_z(k,l) = centroid_z / eh_area2(k,l)
+  end if
+
+! Write an info message with the calculated centroids.
+  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
+                                       eh_centroid2_x(sn,l), &
+                                       eh_centroid2_y(sn,l), &
+                                       eh_centroid2_z(sn,l)
+  call CCTK_INFO(info_message)
+
+end subroutine EHFinder_IntegrateCentroid
+
+
+subroutine EHFinder_IntegrateCircumference(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: int_var, sn, k, l
+  character(len=50) :: info_message
+  CCTK_INT, dimension(1) :: lbnd1, ubnd1, lsh1
+  CCTK_INT, dimension(2) :: lbnd, lsh, nghost
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha
+  CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ
+
+! If finding of surface failed do not try to integrate but exit.
+  if ( find_surface_status .lt. 0 ) then
+    return
+  endif
+
+! Store the levelset_counter in a shorter named variable for convenience.
+  l = levelset_counter
+
+  sn = surface_counter - integrate_counter
+
+  ! Find out the lower bounds of the distributed integration grid arrays.
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' )
+  end if
+
+  ! Find out the local size of the distributed integration grid arrays
+  call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
+  end if
+
+  ! Find out the bounding box of the distributed integration grid arrays
+  call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get bounding box for surface arrays' )
+  end if
+
+  ! Find out the ghost size of the distributed integration grid arrays
+  call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, &
+                                       'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get ghost zones for surface arrays' )
+  end if
+ 
+! Initialize the lower and upper local index for the theta and phi arrays.
+  ithl = 1; ithr = lsh(1)
+  jphl = 1; jphr = lsh(2)
+  
+! If the boundaries are internal, adjust the lower and upper local index
+! to take into account the ghost cells.
+  if ( bbox(1) .eq. 0 ) ithl = ithl + nghost(1)
+  if ( bbox(2) .eq. 0 ) ithr = ithr - nghost(1)
+  if ( bbox(3) .eq. 0 ) jphl = jphl + nghost(2)
+  if ( bbox(4) .eq. 0 ) jphr = jphr - nghost(2)
+
+! Write and info message about what we are doing.
+  call CCTK_INFO ( 'Integrating equatorial circumference' )
+
+! Setup the integration array for the equatorial circumference
+! calculation. Note this is still a 2d-array, so we have to
+! do the sum only for the correct 1d-slice of the array.
+  int_tmp = phi_sym_factor * phiweights * dlphi
+
+! If we are on a processor that contains part of the correct 1-d slice
+! sum up the appropriate elemets, otherwise this processor will
+! contribute zero.
+  if ( ( ithl+lbnd(1) .le. githeta ) .and. ( githeta .le. ithr+lbnd(1) ) ) then
+    itha = githeta - lbnd(1)
+    eq_circ_loc = sum ( int_tmp(itha,jphl:jphr) )
+  else
+    eq_circ_loc = zero
+  end if
+
+! Then reduce the results over all processors and send the result to
+! all processors.
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
+                              eq_circ_loc, eq_circ, CCTK_VARIABLE_REAL )
+
+! Get the upper and lower bounds for the circumference arrays.
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd1, 'ehfinder::eh_circ_eq2' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get lower bounds for area array' )
+  end if
+  call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd1, 'ehfinder::eh_circ_eq2' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get upper bounds for area array' )
+  end if
+  call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh1, 'ehfinder::eh_circ_eq2' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for area array' )
+  end if
+
+! If we are on the correct processor store the result.
+  if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then
+    k = sn - lbnd1(1)
+    eh_circ_eq2(k,l) = eq_circ
+  end if
+
+! Write an info message with the calculated circunference.
+  write(info_message,'(a35,f10.5)') 'Horizon equatorial circumference = ', &
+                                    eq_circ
+  call CCTK_INFO(info_message)
+
+! Now do the same for the polar circomference.
+  call CCTK_INFO ( 'Integrating polar circumference' )
+
+! Setup the integration array for the polar circumference
+! calculation. Note this is still a 2d-array, so we have to
+! do the sum only for the correct 1d-slice of the array.
+  int_tmp = theta_sym_factor * thetaweights * dltheta
+
+! If we are on a processor that contains part of the correct 1-d slice
+! sum up the appropriate elemets, otherwise this processor will
+! contribute zero.
+  if ( ( jphl+lbnd(2) .le. gjphi ) .and. ( gjphi .le. jphr+lbnd(2) ) ) then
+    jpha = gjphi - lbnd(2)
+    pol_circ_loc = sum ( int_tmp(ithl:ithr,jpha) )
+  else
+    pol_circ_loc = zero
+  end if
+
+! Then reduce the results over all processors and send the result to
+! all processors.
+  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
+                              pol_circ_loc, pol_circ, CCTK_VARIABLE_REAL )
+
+
+! If we are on the correct processor store the result.
+  if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then
+    k = sn - lbnd1(1)
+    eh_circ_pol2(k,l) = pol_circ
+  end if
+
+! Write an info message with the calculated circumference.
+  write(info_message,'(a30,f10.5)') 'Horizon polar circumference = ',pol_circ
+  call CCTK_INFO(info_message)
+  
+end subroutine EHFinder_IntegrateCircumference
+
+
+! This routine sets up the weights for the Simpsons rule integration
+! over the surface.
+subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k, im, jm
+
+  CCTK_INT, dimension(4) :: bbox
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
+
+! Find out the lower bounds of the distributed integration grid arrays.
+  call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' )
+  end if
+
+! Find out the local size of the distributed integration grid arrays
+  call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, 'ehfinder::surface_arrays' )
+  if ( ierr .lt. 0 ) then
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
+  end if
+
+  weights = one
+
+! Initialise the weight grid array for the 2D Simpsons rule integration.
+! To do this I need to figure out the global location of the given point.
+! There are 3 cases in the one dimensional case. If the point is on the
+!  boundary the weight is 1/3. If it is at an even position the weight
+! is 4/3 and if it is at an odd position the weight is 2/3. 
+
+  do j = 1, lsh(2)
+
+!   This is first done in the phi direction. Meaning that all points with
+!   the same theta coordinate are set to the same weight.
+    if ( ( lbnd(2)+j .eq. 1 ) .or. ( lbnd(2)+j .eq. nphi ) ) then
+      weights(:,j) = onethird
+      phiweights(:,j) = onethird
+    else if ( mod(lbnd(2)+j,2) .eq. 0 ) then
+      weights(:,j) = fourthirds
+      phiweights(:,j) = fourthirds
+    else
+      weights(:,j) = twothirds
+      phiweights(:,j) = twothirds
+    end if
+
+!   Then it is done in the theta direction with the one-dimensional
+!   weights beeing multiplied.
+    do i = 1, lsh(1)
+      if ( ( lbnd(1)+i .eq. 1 ) .or. ( lbnd(1)+i .eq. ntheta ) ) then
+        weights(i,j) = onethird * weights(i,j)
+        thetaweights(i,j) = onethird
+      else if ( mod(lbnd(1)+i,2) .eq. 0 ) then
+        weights(i,j) = fourthirds * weights(i,j)
+        thetaweights(i,j) = fourthirds
+      else
+        weights(i,j) = twothirds * weights(i,j)
+        thetaweights(i,j) = twothirds
+      end if
+    end do
+  end do
+
+! The end result is a 2D array with the weights in the following pattern.
+
+!  1/9   4/9  2/9   4/9  2/9   4/9  1/9
+!  4/9  16/9  8/9  16/9  8/9  16/9  4/9
+!  2/9   8/9  4/9   8/9  4/9   8/9  2/9
+!  4/9  16/9  8/9  16/9  8/9  16/9  4/9
+!  2/9   8/9  4/9   8/9  4/9   8/9  2/9
+!  4/9  16/9  8/9  16/9  8/9  16/9  4/9
+!  1/9   4/9  2/9   4/9  2/9   4/9  1/9
+end subroutine EHFinder_InitWeights
+
+
+! This is the routine where the stored areas (in eh_area2) are copied
+! into the trigger variable eh_area.
+subroutine EHFinder_CopyArea(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+! If finding of surface failed set area to zero.
+  if ( find_surface_status .lt. 0 ) then
+    eh_area = zero
+    return
+  else
+    call CCTK_INFO('Copying area data')
+    eh_area = eh_area2
+  endif
+end subroutine EHFinder_CopyArea
+
+
+! This is the routine where the stored centroids (in eh_centroid2_[xyz])
+! are copied into the trigger variable eh_centroid_[xyz].
+subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+! If finding of surface failed set centroids to zero.
+  if ( find_surface_status .lt. 0 ) then
+    eh_centroid_x = zero
+    eh_centroid_y = zero
+    eh_centroid_z = zero
+    return
+  else
+    call CCTK_INFO('Copying centroid data')
+    eh_centroid_x = eh_centroid2_x
+    eh_centroid_y = eh_centroid2_y
+    eh_centroid_z = eh_centroid2_z
+  end if
+
+end subroutine EHFinder_CopyCentroid
+
+
+! This is the routine where the stored circumferences (in eh_circ_[eq,pol]2)
+! are copied into the trigger variable eh_circ_[eq,pol].
+subroutine EHFinder_CopyCircumference(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  if ( find_surface_status .lt. 0 ) then
+    eh_circ_eq = zero
+    eh_circ_pol = zero
+    return
+  else
+    call CCTK_INFO('Copying circumference data')
+    eh_circ_eq = eh_circ_eq2
+    eh_circ_pol = eh_circ_pol2
+  end if
+
+end subroutine EHFinder_CopyCircumference