Added support for finding out how many surfaces are present and then do the

integrals over all surfaces. Still only supports full mode. Only weakly tested.
There is a problem with the outputs in multiprocessor mode. The results of
the integration are stored in 1D grid arrays. However if I use DISTRIB=DEFAULT
(taking care only to store results on the correct processor), the
output fails if the number of processors are bigger than the number of
elements in the grid array (essentially some processors contains a chunk
of length zero). If I on the other hand uses DISTRIB=CONSTANT, it works
but on multiple processors i get multiple copies of the data in the output.


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

5 files changed, 232 insertions, 70 deletions

diff --git a/interface.ccl b/interface.ccl
index f4fe515..9af0506 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -53,7 +53,8 @@ CCTK_INT surface_counter TYPE=SCALAR
 CCTK_INT points_counter TYPE=SCALAR
 CCTK_INT more_surfaces TYPE=SCALAR
 CCTK_INT more_points TYPE=SCALAR
- 
+CCTK_INT integrate_counter TYPE=SCALAR 
+
 CCTK_REAL surface_arrays TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=ntheta,nphi GHOSTSIZE=n_array_ghosts,n_array_ghosts DISTRIB=DEFAULT
 {
   ctheta, cphi, rsurf, sintheta, costheta, sinphi, cosphi,
@@ -90,9 +91,31 @@ CCTK_INT rep_mask TYPE=GF TIMELEVELS=1
 
 CCTK_INT re_param_control_pde TYPE=SCALAR
 CCTK_INT re_param_control_approx TYPE=SCALAR
-CCTK_REAL eh_area TYPE=SCALAR
-CCTK_REAL eh_area2 TYPE=SCALAR
-CCTK_REAL eh_centroid TYPE=SCALAR
+
+#CCTK_REAL eh_area TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT
+
+CCTK_REAL eh_area TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
+
+#CCTK_REAL eh_area[maximum_surface_number] TYPE=SCALAR
+
+#CCTK_REAL eh_area2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number GHOSTSIZE=0 DISTRIB=DEFAULT #"The area of the event horizon(s)"
+CCTK_REAL eh_area2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT #"The area of the event horizon(s)"
+
+#CCTK_REAL eh_area2[maximum_surface_number] TYPE=SCALAR
+#CCTK_REAL eh_area3 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
+
+CCTK_REAL eh_centroid TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
 {
   eh_centroid_x, eh_centroid_y, eh_centroid_z
-} "The centroid of the event horizon"
+} "The centroid of the event horizon(s)"
+
+CCTK_REAL eh_centroid2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
+{
+  eh_centroid2_x, eh_centroid2_y, eh_centroid2_z
+} "The centroid of the event horizon(s)"
+
+#CCTK_REAL eh_centroid3 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_number DISTRIB=CONSTANT
+#{
+#  eh_centroid3_x, eh_centroid3_y, eh_centroid3_z
+#} "The centroid of the event horizon(s)"
+
diff --git a/param.ccl b/param.ccl
index ca7335f..b48dd45 100644
--- a/param.ccl
+++ b/param.ccl
@@ -178,7 +178,8 @@ INT n_array_ghosts "Number of ghost points in the surface grid array"
 1:               :: "Positive please"
 } 1
 
-INT maximum_surface_number "The maximum number of surfaces expected in the data"{
+INT maximum_surface_number "The maximum number of surfaces expected in the data"
+{
   1:*            :: "Positive please"
 } 1
 
diff --git a/schedule.ccl b/schedule.ccl
index 572ea93..7652251 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -17,7 +17,8 @@ else
   STORAGE: mask_bak
   STORAGE: surface_arrays
   STORAGE: center_arrays
-  STORAGE: eh_area2
+#  STORAGE: eh_area, eh_centroid
+  STORAGE: eh_area2, eh_centroid2
 #  STORAGE: surface_index
 }
 
@@ -88,77 +89,99 @@ schedule EHFinder_InitWeights at CCTK_POSTINITIAL
 schedule GROUP EHFinder_Surfaces at CCTK_ANALYSIS
 {
   STORAGE: more_surfaces, surface_counter, more_points, points_counter
+  STORAGE: integrate_counter
   STORAGE: surface_index
-  TRIGGER: eh_area2, eh_centroid
-} "Schedule group for counting the number of surfaces and integrating over them"
+#  STORAGE: eh_area2, eh_centroid2
+  TRIGGERS: eh_area, eh_centroid
+} "Count the number of surfaces and integrate over them"
 
 
 schedule EHFinder_CountSurfacesInit in EHFinder_Surfaces
 {
   LANG: Fortran
-  TRIGGER: eh_area2, eh_centroid
-} "Initialising while loop control for the group EHFinder_CountMarkSurfaces"
+  TRIGGERS: eh_area, eh_centroid
+} "Initialise while loop control"
 
 schedule GROUP EHFinder_CountMarkSurfaces in EHFinder_Surfaces after EHFinder_CountSurfacesInit WHILE ehfinder::more_surfaces
 {
-  TRIGGER: eh_area2, eh_centroid
-} "Scedule group for counting and marking surfaces"
+  TRIGGERS: eh_area, eh_centroid
+} "Counting and mark surfaces"
 
 schedule EHFinder_CountSurfaces in EHFinder_CountMarkSurfaces
 {
   LANG: Fortran
-  TRIGGER: eh_area2, eh_centroid
+  TRIGGERS: eh_area, eh_centroid
   SYNC: surface_index
-} "Check if there are more surfaces and initialise while loop control for EHFinder_MarkSurfaces"
+} "Check if there are more surfaces"
 
 schedule EHFinder_MarkSurfaces in EHFinder_CountMarkSurfaces after EHFinder_CountSurfaces WHILE ehfinder::more_points
 {
   LANG: Fortran
-  TRIGGER: eh_area2, eh_centroid
+  TRIGGERS: eh_area, eh_centroid
   SYNC: surface_index
-} "Mark points inside the current surface until they are all marked"
+} "Mark points inside the current surface"
 
 schedule EHFinder_InfoSurfaces in EHFinder_Surfaces after EHFinder_CountMarkSurfaces
 {
   LANG: Fortran
-  TRIGGER: eh_area2, eh_centroid
-} "Output input about found surfaces"
+  TRIGGERS: eh_area, eh_centroid
+} "Output info about found surfaces"
 
-schedule EHFinder_FindSurface at CCTK_ANALYSIS
+schedule group EHFinder_Integration in EHFinder_Surfaces after EHFinder_InfoSurfaces while ehfinder::integrate_counter
 {
-  LANG: Fortran
-  STORAGE: surface_tmp_arrays, surface_int_array
-  TRIGGER: eh_area2, eh_centroid
-} "Find Surface"
+  TRIGGERS: eh_area, eh_centroid
+} "Find and integrate over surfaces"
 
-schedule EHFinder_Integrate at CCTK_ANALYSIS after EHFinder_FindSurface
+schedule EHFinder_FindSurface in EHFinder_Integration
 {
   LANG: Fortran
-  STORAGE: eh_area
-  TRIGGER: eh_area
-} "Calculate integrals"
+  STORAGE: surface_tmp_arrays, surface_int_array
+  TRIGGERS: eh_area, eh_centroid
+} "Find Surface"
 
-schedule EHFinder_FindSurfaceElement at CCTK_ANALYSIS after EHFinder_FindSurface
+schedule EHFinder_FindSurfaceElement in EHFinder_Integration after EHFinder_FindSurface
 {
   LANG: Fortran
   STORAGE: surface_tmp_arrays, interp_metric_arrays 
-  TRIGGER: eh_area2, eh_centroid
+  TRIGGERS: eh_area, eh_centroid
 } "Find Surface Area Element"
 
-schedule EHFinder_IntegrateArea at CCTK_ANALYSIS after EHFinder_FindSurfaceElement
+schedule EHFinder_IntegrateArea in EHFinder_Integration after EHFinder_FindSurfaceElement
 {
   LANG: Fortran
-  STORAGE: integrate_tmp_array
-  TRIGGER: eh_area2, eh_centroid
+  STORAGE: integrate_tmp_array 
+  TRIGGERS: eh_area, eh_centroid
 } "Calculate area integrals"
 
-schedule EHFinder_IntegrateCentroid at CCTK_ANALYSIS after EHFinder_IntegrateArea
+schedule EHFinder_IntegrateCentroid in EHFinder_Integration after EHFinder_IntegrateArea
 {
   LANG: Fortran
-  STORAGE: eh_centroid, surface_tmp_arrays, integrate_tmp_array
-  TRIGGER: eh_centroid
+  STORAGE: surface_tmp_arrays, integrate_tmp_array
+  TRIGGERS: eh_centroid
 } "Calculate centroid integrals"
 
+schedule EHFinder_CopyArea at CCTK_ANALYSIS after EHFinder_Surfaces
+{
+  LANG: Fortran
+  STORAGE: eh_area
+  TRIGGERS: eh_area
+} "Copy areas to output variable"
+
+schedule EHFinder_CopyCentroid at CCTK_ANALYSIS after EHFinder_Surfaces
+{
+  LANG: Fortran
+  STORAGE: eh_centroid
+  TRIGGERS: eh_centroid
+} "Copy areas to output variable"
+
+#schedule EHFinder_Integrate at CCTK_ANALYSIS
+#schedule EHFinder_Integrate at CCTK_ANALYSIS
+#{
+#  LANG: Fortran
+#  STORAGE: eh_area
+#  TRIGGERS: eh_area
+#} "Calculate integrals"
+
 # Read in the data used in reconstructing the 4-metric if necessary
 
 if (CCTK_Equals(eh_metric_type,"numerical"))
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index fe41fd1..41bbac4 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -15,11 +15,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i, j, k, im, jm
+  CCTK_INT :: i, j, k, im, jm, sn
   CCTK_REAL, parameter :: eps = 1.0d-10
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
   CCTK_INT :: interp_handle2
 
+  character(len=30) :: info_message
   CCTK_INT, dimension(4) :: bbox
   CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl
 
@@ -42,7 +43,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 ! the 3D grid functions.
 #include "include/physical_part.h"
 
-  call CCTK_INFO ( 'Finding points on surface' )
+  sn = surface_counter - integrate_counter + 1
+
+  write(info_message,'(a26,i3)') 'Finding points on surface ',sn
+  call CCTK_INFO ( info_message )
 
   ! Find information about the local and global properties of the 2D
   ! Grid arrays.
@@ -77,26 +81,33 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   end if
   
   ! If the user defined center point should be used...
-  if ( use_user_center ) then
+  if ( use_user_center .gt. 0 ) then
     center(1) = center_x
     center(2) = center_y
     center(3) = center_z
   else
 
     ! IF not find an approximation for it by averaging the coodinates that
-    ! are within one grid cell distance from the surface.
+    ! are within one grid cell distance from the surface in question.
     nc_loc = 0; center_loc = zero
     delta = minval ( cctk_delta_space )
     do k = kzl, kzr
       do j = jyl, jyr
         do i = ixl, ixr
           if ( eh_mask(i,j,k) .eq. 0 ) then
-            if ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. &
-                 ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) then
+            if ( ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. &
+                   ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) .and. &
+                 ( ( sc(i,j,k) .eq. sn ) .or. &
+                   ( sc(i-1,j,k) .eq. sn ) .or. & 
+                   ( sc(i+1,j,k) .eq. sn ) .or. &
+                   ( sc(i,j-1,k) .eq. sn ) .or. &
+                   ( sc(i,j+1,k) .eq. sn ) .or. &
+                   ( sc(i,j,k-1) .eq. sn ) .or. &
+                   ( sc(i,j,k+1) .eq. sn ) ) ) then
               nc_loc = nc_loc + 1
               center_loc(1) = center_loc(1) + x(i,j,k)
               center_loc(2) = center_loc(2) + y(i,j,k)
@@ -115,7 +126,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                         CCTK_VARIABLE_REAL )
     ncinv = one / nc
     center = ncinv * center
-!    print*,center
+!    print*,'EHFinder_FindSurfaces : ', sn, center
 
     ! For now symmetries are not handled properly
     if ( symx ) center(1) = zero
@@ -183,8 +194,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   ! Get the 3D coordinate system handle.
   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 ??" )
+    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system.  Forgot to activate an implementation providing coordinates ??" )
   endif
 
   ! Get pointers to the grid arrays containing the interpolation points.
@@ -427,6 +437,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
   call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
 
+  integrate_counter = integrate_counter - 1
 !  print*,'Calculated derivatives'
 !  print*,rsurf
   call Util_TableDestroy ( status, table_handle )
@@ -634,4 +645,7 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS)
   end if
   call CCTK_INFO(info_message)
 
+  integrate_counter = surface_counter
+
+!  print*,'EHFinder_InfoSurfaces : ',integrate_counter
 end subroutine EHFinder_InfoSurfaces
diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90
index f82cd54..e463a43 100644
--- a/src/EHFinder_Integrate2.F90
+++ b/src/EHFinder_Integrate2.F90
@@ -72,8 +72,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   
   call CCTK_InterpHandle ( interp_handle, "Lagrange polynomial interpolation" )
   if ( interp_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot get handle for interpolation. &
- &  Forgot to activate an implementation providing interpolation operators??" )
+    call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
   end if
 
   call Util_TableCreateFromString ( table_handle, "order=3" )
@@ -83,8 +82,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS)
   
   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 ??" )
+    call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" )
   endif
 
   interp_coords(1) = CCTK_PointerTo(interp_x)
@@ -201,10 +199,14 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: int_var
+  CCTK_INT :: int_var, sn
+  CCTK_REAL :: eh_area_tmp
   character(len=30) :: info_message
+!  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
 
   call CCTK_INFO ( 'Integrating area' )
+
+  sn = surface_counter - integrate_counter
   int_tmp = weights * da
 
   call CCTK_VarIndex ( int_var, "ehfinder::int_tmp" )
@@ -212,10 +214,33 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS)
     Call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' )
   end if
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_area2, 1, int_var )
-  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2
+                     eh_area2(sn), 1, int_var )
+  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2(sn)
+
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     eh_area_tmp, 1, int_var )
+!  write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp
+
   call CCTK_INFO(info_message)
 
+!  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 ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then
+!    eh_area2(sn-lbnd(1)) = eh_area_tmp
+!  end if
+!  print*,lbnd,ubnd,lsh
+
 end subroutine EHFinder_IntegrateArea
 
 
@@ -229,10 +254,15 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: int_var
+  CCTK_INT :: int_var, sn, k
   character(len=50) :: info_message
+!  CCTK_INT, dimension(1) :: lbnd, ubnd, lsh
+  CCTK_REAL :: centroid_x, centroid_y, centroid_z
 
   call CCTK_INFO ( 'Integrating centroid' )
+
+  sn = surface_counter - integrate_counter
+
   interp_x = center(1) + rsurf * sintheta * cosphi
   interp_y = center(2) + rsurf * sintheta * sinphi
   interp_z = center(3) + rsurf * costheta
@@ -244,31 +274,66 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS)
 
   int_tmp = interp_x * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid_x, 1, int_var )
+                     eh_centroid2_x(sn), 1, int_var )
 
   int_tmp = interp_y * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid_y, 1, int_var )
+                     eh_centroid2_y(sn), 1, int_var )
 
   int_tmp = interp_z * weights * da
   call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
-                     eh_centroid_z, 1, int_var )
-
-  eh_centroid_x = eh_centroid_x / eh_area2
-  eh_centroid_y = eh_centroid_y / eh_area2
-  eh_centroid_z = eh_centroid_z / eh_area2
-
-  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ',eh_centroid_x, &
-                                                             eh_centroid_y, &
-                                                             eh_centroid_z
+                     eh_centroid2_z(sn), 1, int_var )
+
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     centroid_x, 1, int_var )
+!
+!  int_tmp = interp_y * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     centroid_y, 1, int_var )
+!
+!  int_tmp = interp_z * weights * da
+!  call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, &
+!                     centroid_z, 1, int_var )
+!
+!  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
+
+    eh_centroid2_x(sn) = eh_centroid2_x(sn) / eh_area2(sn)
+    eh_centroid2_y(sn) = eh_centroid2_y(sn) / eh_area2(sn)
+    eh_centroid2_z(sn) = eh_centroid2_z(sn) / eh_area2(sn)
+
+!  if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then
+!    k = sn - lbnd(1)
+!    eh_centroid2_x(k) = eh_centroid2_x(k) / eh_area2(k)
+!    eh_centroid2_y(k) = eh_centroid2_y(k) / eh_area2(k)
+!    eh_centroid2_z(k) = eh_centroid2_z(k) / eh_area2(k)
+!  end if
+
+  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
+                                       eh_centroid2_x(sn), &
+                                       eh_centroid2_y(sn), &
+                                       eh_centroid2_z(sn)
+
+!  write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', &
+!                                       centroid_x, &
+!                                       centroid_y, &
+!                                       centroid_z
   call CCTK_INFO(info_message)
 
 end subroutine EHFinder_IntegrateCentroid
 
 
 ! This routine sets up the weights for the Simpsons rule integration
-
-! This routine sets up the weights for the Simpsons rule integration
 ! over the surface.
 subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
 
@@ -340,3 +405,39 @@ subroutine EHFinder_InitWeights(CCTK_ARGUMENTS)
   !  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
+
+
+subroutine EHFinder_CopyArea(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  call CCTK_INFO('Copying area data')
+
+  eh_area = eh_area2
+
+end subroutine EHFinder_CopyArea
+
+
+subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  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 subroutine EHFinder_CopyCentroid