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

14 files changed, 1269 insertions, 710 deletions

diff --git a/src/EHFinder_Check.F90 b/src/EHFinder_Check.F90
index 22b3703..1fd07d4 100644
--- a/src/EHFinder_Check.F90
+++ b/src/EHFinder_Check.F90
@@ -17,7 +17,7 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
 
   CCTK_INT :: i, j, k, l
-  CCTK_REAL :: fmin_extremum, fmin_extremum_loc, fmax_loc, fmin_bound_loc
+  CCTK_REAL :: fmin_extremum, fmin_extremum_loc, fmax_loc
   character(len=128) :: info_message
   logical :: check_re_init
   CCTK_INT, dimension(3) :: ex_loc
@@ -45,15 +45,20 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS)
 !     values in the extremum.
       fmax_loc = zero
   
+!     On the local processor find the maximum value of f and the maximum value
+!     of f in any local extrema where f is negative.
       fmin_extremum_loc = ex_value
       do k = kzl, kzr
         do j = jyl, jyr
           do i = ixl, ixr
             fmax_loc = max ( f(i,j,k,l), fmax_loc )
             if ( f(i,j,k,l) .lt. zero ) then
-              if ( ( (f(i,j,k,l)-f(i-1,j,k,l))*(f(i+1,j,k,l)-f(i,j,k,l)).le.zero ).and. &
-                   ( (f(i,j,k,l)-f(i,j-1,k,l))*(f(i,j+1,k,l)-f(i,j,k,l)).le.zero ).and. &
-                   ( (f(i,j,k,l)-f(i,j,k-1,l))*(f(i,j,k+1,l)-f(i,j,k,l)).le.zero ) ) then
+              if ( ( (f(i,j,k,l) - f(i-1,j,k,l)) * &
+                     (f(i+1,j,k,l) - f(i,j,k,l)) .le. zero ) .and. &
+                   ( (f(i,j,k,l) - f(i,j-1,k,l)) * &
+                     (f(i,j+1,k,l) - f(i,j,k,l)) .le. zero ) .and. &
+                   ( (f(i,j,k,l) - f(i,j,k-1,l)) * &
+                     (f(i,j,k+1,l) - f(i,j,k,l)) .le. zero ) ) then
                 fmin_extremum_loc = max ( f(i,j,k,l), fmin_extremum_loc )
               end if
             end if
@@ -61,6 +66,7 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS)
         end do
       end do
   
+!     Maximum reduction over all processors.
       call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, &
                          fmin_extremum_loc, fmin_extremum, CCTK_VARIABLE_REAL )
   
@@ -68,15 +74,18 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS)
         call CCTK_WARN(0,'Reduction of fmax_extremum failed')
       end if
   
+!     Write an info message about the local extremum.
       write(info_message,'(a39,i3,a6,f12.5)') &
           'Minimum f at the extrema for level set ', l, ' is : ',fmin_extremum
       call CCTK_INFO(info_message)
   
+!     If the local extremum indicate and imminent pinch off undo the
+!     re-initialization if this is requested by the user.
       if ( re_init_undo .gt. 0 ) then
         if ( fmin_extremum .gt. -two * minval(cctk_delta_space) ) then
-          write(info_message,'(a45,i3,a6,f12.5)') &
+          write(info_message,'(a45,i3,a33)') &
                'Re-initialization of level set ', l, &
-               'undone due to imminent pinch-off'
+               ' undone due to imminent pinch-off'
           call CCTK_INFO(info_message)
           f(:,:,:,l) = fbak(:,:,:,l)
           eh_mask(:,:,:,l) = eh_mask_bak(:,:,:,l)
diff --git a/src/EHFinder_Constants.F90 b/src/EHFinder_Constants.F90
index fdd5b6f..2162631 100644
--- a/src/EHFinder_Constants.F90
+++ b/src/EHFinder_Constants.F90
@@ -10,13 +10,10 @@ module EHFinder_Constants
   CCTK_REAL, parameter :: three = 3.0d0
   CCTK_REAL, parameter :: four = 4.0d0
   CCTK_REAL, parameter :: eight = 8.0d0
-  CCTK_REAL, parameter :: ten = 10.0d0
   CCTK_REAL, parameter :: half = 0.5d0
   CCTK_REAL, parameter :: onethird = one / three
   CCTK_REAL, parameter :: twothirds = two / three
   CCTK_REAL, parameter :: fourthirds = four / three
   CCTK_REAL, parameter :: quarter = 0.25d0
-  CCTK_REAL, parameter :: tenth = 0.1d0
-  CCTK_REAL, parameter :: huge = 1d23
   CCTK_REAL, parameter :: pi = 3.1415926535897932385d0
 end module EHFinder_Constants
diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90
index 2537408..d53db6e 100644
--- a/src/EHFinder_FindSurface.F90
+++ b/src/EHFinder_FindSurface.F90
@@ -21,13 +21,14 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   CCTK_INT :: interp_handle2
 
   character(len=80) :: info_message
+  character(len=128) :: warn_message
   character(len=200) :: surface_interp
   CCTK_INT :: surface_interp_len
   character(len=7) :: surface_order
   character(len=15) :: vname
 
   CCTK_INT, dimension(4) :: bbox
-  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl
+  CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost
 
   CCTK_REAL, dimension(3) :: center_loc
   CCTK_INT :: nc_loc, nc
@@ -55,52 +56,52 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 ! set to -1.
   find_surface_status = 0
   sn = surface_counter - integrate_counter + 1
+
+! Set the index to the current level set function.
   l = levelset_counter
+
+! Write an info message about what is going on.
   write(info_message,'(a26,i3,a14,i3)') 'Finding points on surface ', sn, &
                                  ' in level set ',l
   call CCTK_INFO ( info_message )
 
-  ! Find information about the local and global properties of the 2D
-  ! Grid arrays.
-  call CCTK_GroupbboxGN ( status, cctkGH, 4, bbox, "ehfinder::surface_arrays" )
+! Find information about the local and global properties of the 2D
+! Grid arrays.
+  call CCTK_GroupbboxGN ( status, cctkGH, 4, bbox, 'ehfinder::surface_arrays' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get bounding box for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get bounding box for surface arrays' )
   end if
-!  print*,bbox
-  call CCTK_GroupgshGN ( status, cctkGH, 2, gsh, "ehfinder::surface_arrays" )
+  call CCTK_GroupgshGN ( status, cctkGH, 2, gsh, 'ehfinder::surface_arrays' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get global size for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get global size for surface arrays' )
   end if
-!  print*,gsh
-  call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd, "ehfinder::surface_arrays" )
+  call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' )
   end if
-!  print*,lbnd
-  call CCTK_GroupubndGN ( status, cctkGH, 2, ubnd, "ehfinder::surface_arrays" )
+  call CCTK_GroupubndGN ( status, cctkGH, 2, ubnd, 'ehfinder::surface_arrays' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get upper bounds for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get upper bounds for surface arrays' )
   end if
-!  print*,ubnd
-  call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, "ehfinder::surface_arrays" )
+  call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, 'ehfinder::surface_arrays' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
   end if
-!  print*,lsh
-  call CCTK_GroupnghostzonesGN ( status, cctkGH, 2, nghost, "ehfinder::surface_arrays" )
+  call CCTK_GroupnghostzonesGN ( status, cctkGH, 2, nghost, &
+                                         'ehfinder::surface_arrays' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
   end if
   
-  ! If the user defined center point should be used...
+! If the user defined center point should be used...
   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 in question.
+!   If not find an approximation for it by averaging the coodinates that
+!   are within one grid cell distance from the surface in question.
     nc_loc = 0; center_loc = zero
     min_delta = minval ( cctk_delta_space )
     do k = kzl, kzr
@@ -130,7 +131,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end do
     end do
 
-    ! Reduce over the processors.
+!   Reduce over the processors.
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, sum_handle, &
                                     nc_loc, nc, CCTK_VARIABLE_INT )
     call CCTK_ReduceLocArrayToArray1D ( status, cctkGH, -1, sum_handle, &
@@ -138,14 +139,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                         CCTK_VARIABLE_REAL )
     ncinv = one / nc
     center = ncinv * center
-!    print*,'EHFinder_FindSurfaces : ', sn, center
-!    print*
   end if
-!  print*,'Center : ',center
-!  print*,'use_user_center : ',use_user_center
 
-  ! Try to figure out the symmetries. First find the range in coordinates
-  ! for the points inside the current surface.
+! Try to figure out the symmetries. First find the range in coordinates
+! for the points inside the current surface.
   crange_min_loc(1) = minval ( x, mask = sc .eq. sn )
   crange_min_loc(2) = minval ( y, mask = sc .eq. sn )
   crange_min_loc(3) = minval ( z, mask = sc .eq. sn )
@@ -153,18 +150,15 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
   crange_max_loc(2) = maxval ( y, mask = sc .eq. sn )
   crange_max_loc(3) = maxval ( z, mask = sc .eq. sn )
 
+! Reduce over all processors.
   call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
                                       crange_min_loc, crange_min_glob, &
                                       3, CCTK_VARIABLE_REAL )
   call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, &
                                       crange_max_loc, crange_max_glob, &
                                       3, CCTK_VARIABLE_REAL )
-!  print*,crange_min_loc,crange_max_loc
-!  print*,crange_min_glob,crange_max_glob
-!  print*
-!  print*,'Symmetries : ', symx, symy, symz
 
-  ! Full mode. Modify if necessary.
+! Full mode. Modify if necessary.
   ltheta = zero; utheta = pi
   lphi = zero; uphi = two * pi
   sym_factor = one
@@ -172,9 +166,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
   if ( CCTK_EQUALS ( domain, 'bitant' ) ) then
 
-    ! Bitant with symmetry in the z-direction
+!   Bitant with symmetry in the z-direction
 
     if ( CCTK_EQUALS ( bitant_plane, 'xy' ) ) then
+!     If the minimum z-coordinate of the points inside the surface is less
+!     than zero, then the surface must share the grid symmetry across the
+!     xy-plane.
       if ( crange_min_glob(3) .lt. zero ) then
         ltheta = zero; utheta = half * pi
         lphi = zero; uphi = two * pi
@@ -184,7 +181,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end if 
     end if
 
-    ! Bitant with symmetry in the y-direction
+!   Bitant with symmetry in the y-direction
   
     if ( CCTK_EQUALS ( bitant_plane, 'xz' ) ) then
       if ( crange_min_glob(2) .lt. zero ) then
@@ -196,7 +193,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end if 
     end if
 
-    ! Bitant with symmetry in the x-direction
+!   Bitant with symmetry in the x-direction
 
     if ( CCTK_EQUALS ( bitant_plane, 'yz' ) ) then
       if ( crange_min_glob(1) .lt. zero ) then
@@ -211,7 +208,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
 
   if ( CCTK_EQUALS ( domain, 'quadrant' ) ) then
 
-    ! Quadrant in x-direction
+!   Quadrant in x-direction
 
     if ( CCTK_EQUALS ( quadrant_direction, 'x' ) ) then
       if ( ( crange_min_glob(3) .lt. zero ) .and. &
@@ -236,7 +233,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end if
     end if
     
-    ! Quadrant in y-direction
+!   Quadrant in y-direction
 
     if ( CCTK_EQUALS ( quadrant_direction, 'y' ) ) then
       if ( ( crange_min_glob(3) .lt. zero ) .and. &
@@ -261,7 +258,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end if
     end if
 
-    ! Quadrant in z-direction
+!   Quadrant in z-direction
 
     if ( CCTK_EQUALS ( quadrant_direction,'z' ) ) then
       if ( ( crange_min_glob(2) .lt. zero ) .and. &
@@ -287,7 +284,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     end if
   end if
     
-  ! Octant
+! Octant
 
   if ( CCTK_EQUALS ( domain, 'octant' ) ) then
     if ( ( crange_min_glob(3) .lt. zero ) .and. &
@@ -340,24 +337,30 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     end if
   end if
 
+! Set the theta index to be used for the equatorial circumference integration
+! based on the symmetries.
   if ( .not. s_symz ) then
     githeta = ( ntheta - 1 ) / 2 + 1
   else
     githeta = ntheta
   end if
   
+! Set the phi index to be used for the polar circumference integration based
+! on the symmetries.
   if ( ( .not. s_symy ) .and. ( s_symx ) ) then
     gjphi = ( nphi - 1 ) / 2 + 1
   else
     gjphi = 1
   end if
 
+! Set up the symmetry muliplication factors based on the range in the
+! angular coordinates.
   theta_sym_factor = two * pi / ( utheta - ltheta )
   phi_sym_factor = two * pi / ( uphi - lphi )
 
-  ! Find dtheta and dphi and initialise the theta and phi grid arrays.
-  ! Here i + lbnd(1) - 1 is the global index for theta and
-  !      j + lbnd(2) - 1 is the global index for phi.
+! Find dtheta and dphi and initialise the theta and phi grid arrays.
+! Here i + lbnd(1) - 1 is the global index for theta and
+!      j + lbnd(2) - 1 is the global index for phi.
 
   dtheta = ( utheta - ltheta ) / ( ntheta - 1 )
   dphi = ( uphi - lphi ) / ( nphi - 1 )
@@ -370,64 +373,71 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     cphi(:,j) = lphi + dphi * ( j +lbnd(2) - 1 )
   end do
 
-!  print*,'Delta : ',dtheta,dphi
-  ! Calculate sines and cosines and store them in grid arrays since they
-  ! are expensive.
+! Calculate sines and cosines and store them in grid arrays since they
+! are computationally expensive and memory for 2D grid arrays is cheap.
   sintheta = sin(ctheta)
   costheta = cos(ctheta)
   sinphi = sin(cphi)
   cosphi = cos(cphi)
   
-  ! rsurf is the radial grid array (as function of theta and phi) and is
-  ! initialised to zero. The grid array drsurf contains the changes to rsurf
-  ! and is initialised to 2*Delta. The grid array n_since_last_reduction is
-  ! a help in keeping track on how long ago drsurf for a given point has been
-  ! reduced and is initialised to -1 to indicate that no reduction has taken
-  ! place.
+! rsurf is the radial grid array (as function of theta and phi) and is
+! initialised to zero. The grid array drsurf contains the changes to rsurf
+! and is initialised to 2*Delta. The grid array n_since_last_reduction is
+! a help in keeping track on how long ago drsurf for a given point has been
+! reduced and is initialised to -1 to indicate that no reduction has taken
+! place.
   drsurf = two * min_delta
   rsurf = zero
   n_since_last_reduction = -1
 
-  ! Get and interpolation handle for required polynomial interpolation.
-  ! If Lagrange polynomial interpolation does not work try Hermite
-  ! polynomial interpolation, which avoids problems with non-continuous
-  !  derivatives in the Newton iteration later on.
+! Get and interpolation handle for required polynomial interpolation.
+! If Lagrange polynomial interpolation does not work try Hermite
+! polynomial interpolation, which avoids problems with non-continuous
+! derivatives in the Newton iteration later on.
 
-  ! First convert the string parameter to a Fortran string.
+! First convert the string parameter to a Fortran string.
   call CCTK_FortranString ( surface_interp_len, surface_interpolator, &
                                                 surface_interp )
 
-  ! Then get the hande ...
+! Then get the handle ...
   call CCTK_InterpHandle ( interp_handle, surface_interp(1:surface_interp_len) )
 
+! The following is done because some C-preprocessors seem to have problems
+! handling character strings that are continued on several lines.
   if ( interp_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
+    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 required interpolation order parameter into a fortran string.
+! Write the required interpolation order parameter into a fortran string.
   write(surface_order,'(a6,i1)') 'order=',surface_interpolation_order
 
-  ! Then create a table from the string.
+! Then create a table from the string.
   call Util_TableCreateFromString ( table_handle, surface_order )
   if ( table_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
   end if
   
-  ! Get the 3D coordinate system handle.
-  call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
+! 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 ??" )
+    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 pointers to the grid arrays containing the interpolation points.
+! Get pointers to the grid arrays containing 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 pointer to the interpolation output array
+! Get pointer to the interpolation output array
   out_array(1) = CCTK_PointerTo(f_interp)
 
-  ! Get the variable index for f
+! Get the variable index for f
   write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1, ']'
   call CCTK_VarIndex ( in_array, vname )
 
@@ -442,22 +452,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     do j = 1, lsh(2)
       do i = 1, lsh(1)
 
-        ! The radius that are tested in this direction.
+!       The radius that are tested in this direction.
         rloc = rsurf(i,j) + drsurf(i,j)
 
-        ! Calculate the interpolation points.
+!       Calculate the interpolation points.
         interp_x(i,j) = center(1) + rloc * sintheta(i,j) * cosphi(i,j)
         interp_y(i,j) = center(2) + rloc * sintheta(i,j) * sinphi(i,j)
         interp_z(i,j) = center(3) + rloc * costheta(i,j)
 
-        ! If the point has been reduced once, increment n_since_last_reduction.
-        ! 
+!       If the point has been reduced once, increment n_since_last_reduction.
         if ( n_since_last_reduction(i,j) .ge. 0 ) then 
           n_since_last_reduction(i,j) = n_since_last_reduction(i,j) + 1
         end if
       end do
     end do
 
+!   Call the interpolator.
     call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
                                  table_handle, coord_system_handle, &
                                  lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
@@ -470,18 +480,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
        return
     end if
                                
+!   Find the maximal change in rsurf.
     maxdr_loc = maxval ( drsurf )
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
                                     maxdr_loc, maxdr, CCTK_VARIABLE_REAL )
 
+!   Find the maximal value of interpolated f values.
     maxf_loc = maxval ( abs ( f_interp ) )
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
                                     maxf_loc, maxf, CCTK_VARIABLE_REAL )
     
-    maxl = maxloc ( abs ( f_interp ) )
-
+!   If the maximal change in rsurf is small enough, we are satisfied that
+!   the Newton iteration will be able to find the root.
     if ( maxdr .lt. min_delta ) exit find_starting_points
 
+!   Check if we have crossed the zero point and if we have reduce the
+!   change in rsurf we are going to try next. 
     do j = 1, lsh(2)
       do i = 1, lsh(1)
         if ( f_interp(i,j) .gt. zero ) then
@@ -497,28 +511,28 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
     end do
     
   end do find_starting_points
-!  print*,'Max values: ',maxval(rsurf),maxval(f_interp)
-!  print*,'Min values: ',minval(rsurf),minval(f_interp)
-!  print*,'found starting points' 
+
+! Now for the Newton iteration we also need the interpolates derivatives.
   out_array(2) = CCTK_PointerTo(dfdx_interp)
   out_array(3) = CCTK_PointerTo(dfdy_interp)
   out_array(4) = CCTK_PointerTo(dfdz_interp)
 
   call Util_TableSetIntArray ( status, table_handle, 4, &
-                               op_indices, "operand_indices" )
+                               op_indices, 'operand_indices' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+    call CCTK_WARN ( 0, 'Cannot set operand indices array in parameter table' )
   endif
   call Util_TableSetIntArray ( status, table_handle, 4, &
-                               op_codes, "operation_codes" )
+                               op_codes, 'operation_codes' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+    call CCTK_WARN ( 0, 'Cannot set operation codes array in parameter table' )
   endif
 
 ! find the surface points in all directions simultaneously.
 
   find_surface_points: do k = 1, 15
 
+!   Calculate 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)
@@ -527,6 +541,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
       end do
     end do
 
+!   Call the interpolator.
     call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
                                  table_handle, coord_system_handle, &
                                  lsh(1) * lsh(2), CCTK_VARIABLE_REAL, &
@@ -539,12 +554,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
        return
     end if
  
+!   Find the maximum value of interpolated f values.
     maxf_loc = maxval ( abs ( f_interp ) )
     call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, &
                                     maxf_loc, maxf, CCTK_VARIABLE_REAL )
     
-    maxl = maxloc ( abs ( f_interp ) )
-
+!   If the maximum interpolated value of f is to big, perform a Newton
+!   iteration
     if ( maxf .gt. eps ) then
       do j = 1, lsh(2)
         do i = 1, lsh(1)
@@ -555,16 +571,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
         end do
       end do
     else
+!     Else we are done
       exit find_surface_points
     end if
-!    call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
 
   end do find_surface_points
-  if ( k .gt. 15 ) call CCTK_WARN(1,'Newton iteration did not converge! Area may be inaccurate')
-!  print*,'found_surface_points'
+
+! If we used to many iterations, the Newton iteration did not converge, so
+! print a warning.  
+  if ( k .gt. 15 ) then
+    warn_message = 'Newton iteration did not converge! Area may be inaccurate'
+    call CCTK_WARN(1, warn_message )
+  end if
 
   drdtheta = zero
   drdphi = zero
+
 ! Calculate the derivatives of r.
 
 ! First for interior points
@@ -661,13 +683,11 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS)
                                        rsurf(im,jm-2) )
   end if
 
-  call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" )
-
+! I do not remember why this is done. I will have to investigate.
   integrate_counter = integrate_counter - 1
-!  print*,'Calculated derivatives'
-!  print*,rsurf
+
+! Get rid of the interpolation table.
   call Util_TableDestroy ( status, table_handle )
-!  print*,'gjphi = ',gjphi
   
 end subroutine EHFinder_FindSurface
 
@@ -718,7 +738,6 @@ subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-!  print*,'EHFinder_CountSurfacesInit called'
   more_surfaces = 1
   surface_counter = 0
   sc = 0
@@ -744,27 +763,27 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
 ! the 3D grid functions.
 #include "include/physical_part.h"
 
-  ! Find the location of the minimum among active points not already found
-  ! to be inside a surface.
+! Find the location of the minimum among active points not already found
+! to be inside a surface.
   minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter), &
                ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. &
                ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) )
 
-  ! Find the value of f at that location.
+! Find the value of f at that location.
   minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1,levelset_counter)
 
-  ! Find the minimum over all processors.
+! Find the minimum over all processors.
   call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, &
                               minf_loc, minf_glob, CCTK_VARIABLE_REAL )
   if ( ierr .ne. 0 ) then
     call CCTK_WARN(0,'Reduction of minf_glob failed')
   end if
 
-  ! To figure out on which processor the minimum is located I set procpos_loc
-  ! to the local processor number if the local minimum is equal to the global
-  ! minimum otherwise I set it to the total number of processors. When I then
-  ! find the minimum of procpos_loc over all processors I get the processor
-  ! containing the minimum.
+! To figure out on which processor the minimum is located I set procpos_loc
+! to the local processor number if the local minimum is equal to the global
+! minimum otherwise I set it to the total number of processors. When I then
+! find the minimum of procpos_loc over all processors I get the processor
+! containing the minimum.
   if ( minf_loc .eq. minf_glob ) then
     procpos_loc = CCTK_MyProc ( cctkGH )
   else
@@ -777,8 +796,8 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS)
     call CCTK_WARN(0,'Reduction of minf_glob failed')
   end if
 
-  ! If minf_glob<0 there must exist another surface so lets set up the variables
-  ! to start marking and counting these points.
+! If minf_glob<0 there must exist another surface so lets set up the variables
+! to start marking and counting these points.
   if ( minf_glob .lt. zero ) then
     surface_counter = surface_counter + 1
     points_counter = 0
@@ -813,13 +832,13 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
 ! the 3D grid functions.
 #include "include/physical_part.h"
   
-  ! Initialise the local point counter.
+! Initialise the local point counter.
   n_loc = 0
 
-  ! Check if any points have a marked point as a neighbour. Only look at
-  ! active points with f<0 that are not already marked. Take care at the
-  ! boundaries (physical and excision) not to look at inactive or not
-  ! exisiting points.
+! Check if any points have a marked point as a neighbour. Only look at
+! active points with f<0 that are not already marked. Take care at the
+! boundaries (physical and excision) not to look at inactive or not
+! exisiting points.
   do k = kzl, kzr
     do j = jyl, jyr
       do i = ixl, ixr
@@ -828,7 +847,7 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
              ( f(i,j,k,levelset_counter) .lt. 0 ) .and. &
              ( sc(i,j,k) .ne. surface_counter ) ) then
 
-          ! Find out if we are on the boundary
+!         Find out if we are on the boundary
           if ( i + cctk_lbnd(1) .eq. 1 ) then
             i1 = 0
           else
@@ -881,31 +900,25 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS)
     end do
   end do
 
-  ! Reduce the number of newly marked points over all processors
+! Reduce the number of newly marked points over all processors
   call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
                               n_loc, n_glob, CCTK_VARIABLE_INT )
   if ( ierr .ne. 0 ) then
     call CCTK_WARN(0,'Reduction of n_glob failed')
   end if
 
-  ! If the total number of newly marked points is 0, set more_points to 0
-  ! to indicate that we are finished. Otherwise continue and add the number
-  ! of points to the points_counter.
+! If the total number of newly marked points is 0, set more_points to 0
+! to indicate that we are finished. Otherwise continue and add the number
+! of points to the points_counter.
   if ( n_glob .eq. 0 ) then
     more_points = 0
   end if
   points_counter = points_counter + n_glob
 
-  ! This should not be necessary, since I have a SYNC in the schedule.ccl,
-  ! but for some reason that does not work.
-  call CCTK_SyncGroup ( ierr, cctkGH, "ehfinder::surface_index" )
-
-  call CartSymVN(ierr,cctkGH,'ehfinder::sc')
-
 end subroutine EHFinder_MarkSurfaces
 
 
-!This routine prints out information about the found surfaces.
+! This routine prints out information about the found surfaces.
 subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -918,7 +931,7 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS)
 
   character(len=80) :: info_message
 
-  ! Write out how many surfaces were found.
+! Write out how many surfaces were found.
   if ( surface_counter .eq. 1 ) then
     write(info_message,'(a9,i3,a22,i3)') 'There is ',surface_counter, &
                                       ' surface in level set ',levelset_counter
diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90
index bb8e350..28ce33a 100644
--- a/src/EHFinder_Generator_Sources.F90
+++ b/src/EHFinder_Generator_Sources.F90
@@ -19,6 +19,7 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
 
   character(len=200) :: gen_interp
+  character(len=128) :: warn_message
   CCTK_INT :: gen_interp_len
   character(len=7) :: gen_order
   character(len=15) :: vname
@@ -37,46 +38,53 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign
   CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
 
+! Set the sign correctly depending on the surface direction.
   if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
   if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
 
   out_types = CCTK_VARIABLE_REAL
 
-  ! Convert the generator_interpolator string parameter to a Fortran string.
+! Convert the generator_interpolator string parameter to a Fortran string.
   call CCTK_FortranString ( gen_interp_len, generator_interpolator, &
                                             gen_interp )
 
-  ! Get the corresponding interpolator handle.
+! Get the corresponding interpolator handle.
   call CCTK_InterpHandle ( interp_handle, gen_interp )
 
   if ( interp_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
+    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
 
-  ! Convert the interpolation order parameter to a Fortran string to be placed
-  ! in the interpolator table. Note that the order is assumed to contain only
-  ! 1 digit.
+! Convert the interpolation order parameter to a Fortran string to be placed
+! in the interpolator table. Note that the order is assumed to contain only
+! 1 digit.
   write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order
 
-  ! Create the table directly from the string.
+! Create the table directly from the string.
   call Util_TableCreateFromString ( table_handle, gen_order )
   if ( table_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
   end if
 
-  ! Get the 3D coordinate system handle.
-  call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
+! 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 ??" )
+    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
 
-  ! Find out how many interpolation points are located on this processor.
-  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" )
+! Find out how many interpolation points are located on this processor.
+  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
   end if
 
-  ! Set the pointers to the output arrays.
+! Set the pointers to the output arrays.
   out_arrays(1) = CCTK_PointerTo(alpg)
   out_arrays(2) = CCTK_PointerTo(betaxg)
   out_arrays(3) = CCTK_PointerTo(betayg)
@@ -92,51 +100,53 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
   out_arrays(13) = CCTK_PointerTo(dfzg)
   out_arrays(14) = CCTK_PointerTo(psig)
 
-  ! Loop over the level sets
+! Loop over the level sets
   do l = 1, eh_number_level_sets
 
-    ! Set the pointers to the points to be interpolated to.
+!   Set the pointers to the points to be interpolated to.
     interp_coords(1) = CCTK_PointerTo(xg(:,l))
     interp_coords(2) = CCTK_PointerTo(yg(:,l))
     interp_coords(3) = CCTK_PointerTo(zg(:,l))
 
-    ! Set the indices to the input grid functions.
-    call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
-    call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
-    call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
-    call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
-    call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
-    call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
-    call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
-    call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
-    call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
-    call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
+!   Set the indices to the input grid functions.
+    call CCTK_VarIndex ( in_arrays(1), 'admbase::alp' )
+    call CCTK_VarIndex ( in_arrays(2), 'admbase::betax' )
+    call CCTK_VarIndex ( in_arrays(3), 'admbase::betay' )
+    call CCTK_VarIndex ( in_arrays(4), 'admbase::betaz' )
+    call CCTK_VarIndex ( in_arrays(5), 'admbase::gxx' )
+    call CCTK_VarIndex ( in_arrays(6), 'admbase::gxy' )
+    call CCTK_VarIndex ( in_arrays(7), 'admbase::gxz' )
+    call CCTK_VarIndex ( in_arrays(8), 'admbase::gyy' )
+    call CCTK_VarIndex ( in_arrays(9), 'admbase::gyz' )
+    call CCTK_VarIndex ( in_arrays(10), 'admbase::gzz' )
     write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1,']'
     call CCTK_VarIndex ( in_arrays(11), vname )
-    call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" )
+    call CCTK_VarIndex ( in_arrays(12), 'staticconformal::psi' )
 
-    ! Check the metric type. At present physical and static_conformal are
-    ! supported.
+!   Check the metric type. At present physical and static_conformal are
+!   supported.
     if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
 
-      ! Set the operand indices table entry, corresponding
-      ! to interpolation of admbase::alp (1), admbase::shift (3),
-      ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for
-      ! a total of 13 output arrays.
+!     Set the operand indices table entry, corresponding
+!     to interpolation of admbase::alp (1), admbase::shift (3),
+!     admbase::metric(6) and the first derivatives of ehfinder::f (3) for
+!     a total of 13 output arrays.
       call Util_TableSetIntArray ( status, table_handle, 13, &
-                                 op_indices(1:13), "operand_indices" )
+                                 op_indices(1:13), 'operand_indices' )
       if ( status .lt. 0 ) then
-        call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+        warn_message = 'Cannot set operand indices array in parameter table' 
+        call CCTK_WARN ( 0, warn_message )
       endif
 
-      ! Set the corresponding table entry for the operation codes.
+!     Set the corresponding table entry for the operation codes.
       call Util_TableSetIntArray ( status, table_handle, 13, &
-                                 op_codes(1:13), "operation_codes" )
+                                 op_codes(1:13), 'operation_codes' )
       if ( status .lt. 0 ) then
-        call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+        warn_message = 'Cannot set operation codes array in parameter table'
+        call CCTK_WARN ( 0, warn_message )
       endif
 
-      ! Call the interpolator.
+!     Call the interpolator.
       call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
                                    table_handle, coord_system_handle, &
                                    lsh(1), CCTK_VARIABLE_REAL, &
@@ -147,14 +157,14 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
         call CCTK_INFO ( 'Interpolation failed.' )
       end if
 
-      ! For each point on this processor calculate the right hand side of the
-      ! characteristic evolution equation.
+!     For each point on this processor calculate the right hand side of the
+!     characteristic evolution equation.
       do i = 1, lsh(1)
 
-        ! calculate the square of the lapse.
+!       calculate the square of the lapse.
         alp2 = alpg(i)**2
 
-        ! Calculate the inverse of the 3-metric.
+!       Calculate the inverse of the 3-metric.
         guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
         guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
         guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
@@ -169,32 +179,34 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
         guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
         guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
       
-        ! Raise the index of the partial derivatives of f.
+!       Raise the index of the partial derivatives of f.
         dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
         dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
         dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
 
-        ! Calculate the overall multiplication factor.
+!       Calculate the overall multiplication factor.
         factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
                                         dfuy * dfyg(i) + &
                                         dfuz * dfzg(i) ) )
 
-        ! Finally obtain dx^i/dt.
+!       Finally obtain dx^i/dt.
         dxg(i,l) = - betaxg(i) + ssign * factor * dfux
         dyg(i,l) = - betayg(i) + ssign * factor * dfuy
         dzg(i,l) = - betazg(i) + ssign * factor * dfuz
       end do
     else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
       call Util_TableSetIntArray ( status, table_handle, 14, &
-                                 op_indices, "operand_indices" )
+                                 op_indices, 'operand_indices' )
       if ( status .lt. 0 ) then
-        call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+        warn_message = 'Cannot set operand indices array in parameter table'
+        call CCTK_WARN ( 0, warn_message )
       endif
 
       call Util_TableSetIntArray ( status, table_handle, 14, &
-                                 op_codes, "operation_codes" )
+                                 op_codes, 'operation_codes' )
       if ( status .lt. 0 ) then
-        call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+        warn_message = 'Cannot set operation codes array in parameter table'
+        call CCTK_WARN ( 0, warn_message )
       endif
 
       call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
@@ -217,7 +229,7 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
         guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
         guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
 
-!       The determinant divided by psi^4.
+!       The inverse of the determinant divided by psi^4.
         idetg = psi4 / ( gxxg(i) * guxx + &
                          gxyg(i) * guxy + &
                          gxzg(i) * guxz )
diff --git a/src/EHFinder_Generator_Sources2.F90 b/src/EHFinder_Generator_Sources2.F90
index b557381..27a9f90 100644
--- a/src/EHFinder_Generator_Sources2.F90
+++ b/src/EHFinder_Generator_Sources2.F90
@@ -19,6 +19,7 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
   CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
 
   character(len=200) :: gen_interp
+  character(len=128) :: warn_message
   CCTK_INT :: gen_interp_len
   character(len=7) :: gen_order
 
@@ -32,82 +33,89 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
   CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign
   CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
 
+! Set the sign depending on the surface direction.
   if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
   if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
 
   out_types = CCTK_VARIABLE_REAL
 
-  ! Convert the generator_interpolator string parameter to a Fortran string.
+! Convert the generator_interpolator string parameter to a Fortran string.
   call CCTK_FortranString ( gen_interp_len, generator_interpolator, &
                                             gen_interp )
 
-  ! Get the corresponding interpolator handle.
+! Get the corresponding interpolator handle.
   call CCTK_InterpHandle ( interp_handle, gen_interp )
 
   if ( interp_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
+    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
 
-  ! Convert the interpolation order parameter to a Fortran string to be placed
-  ! in the interpolator table. Note that the order is assumed to contain only
-  ! 1 digit.
+! Convert the interpolation order parameter to a Fortran string to be placed
+! in the interpolator table. Note that the order is assumed to contain only
+! 1 digit.
   write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order
 
-  ! Create the table directly from the string.
+! Create the table directly from the string.
   call Util_TableCreateFromString ( table_handle, gen_order )
   if ( table_handle .lt. 0 ) then
-    call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
+    call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' )
   end if
 
-  ! Get the 3D coordinate system handle.
-  call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
+! 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 ??" )
+    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
 
 #include "include/physical_part.h"
 
-  ! Find out how many interpolation points are located on this processor.
-  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" )
+! Find out how many interpolation points are located on this processor.
+  call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
+    call CCTK_WARN ( 0, 'cannot get local size for surface arrays' )
   end if
 
-  ! Set the indices to the input grid functions.
-  call CCTK_VarIndex ( in_arrays(1), "ehfinder::xgf" )
-  call CCTK_VarIndex ( in_arrays(2), "ehfinder::ygf" )
-  call CCTK_VarIndex ( in_arrays(3), "ehfinder::zgf" )
+! Set the indices to the input grid functions.
+  call CCTK_VarIndex ( in_arrays(1), 'ehfinder::xgf' )
+  call CCTK_VarIndex ( in_arrays(2), 'ehfinder::ygf' )
+  call CCTK_VarIndex ( in_arrays(3), 'ehfinder::zgf' )
 
-  ! Set the operand indices table entry, corresponding
-  ! to interpolation of ehfinder::generator_gf (3)
+! Set the operand indices table entry, corresponding
+! to interpolation of ehfinder::generator_gf (3)
   call Util_TableSetIntArray ( status, table_handle, 3, &
-                               op_indices, "operand_indices" )
+                               op_indices, 'operand_indices' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
+    call CCTK_WARN ( 0, 'Cannot set operand indices array in parameter table' )
   endif
 
-  ! Set the corresponding table entry for the operation codes.
+! Set the corresponding table entry for the operation codes.
   call Util_TableSetIntArray ( status, table_handle, 3, &
-                               op_codes, "operation_codes" )
+                               op_codes, 'operation_codes' )
   if ( status .lt. 0 ) then
-    call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
+    call CCTK_WARN ( 0, 'Cannot set operation codes array in parameter table' )
   endif
 
-  ! Loop over the level sets
+! Loop over the level sets
   do l = 1, eh_number_level_sets
 
-    ! Set the pointers to the points to be interpolated to.
+!   Set the pointers to the points to be interpolated to.
     interp_coords(1) = CCTK_PointerTo(xg(:,l))
     interp_coords(2) = CCTK_PointerTo(yg(:,l))
     interp_coords(3) = CCTK_PointerTo(zg(:,l))
 
-    ! Set the pointers to the output arrays.
+!   Set the pointers to the output arrays.
     out_arrays(1) = CCTK_PointerTo(dxg(:,l))
     out_arrays(2) = CCTK_PointerTo(dyg(:,l))
     out_arrays(3) = CCTK_PointerTo(dzg(:,l))
 
-    ! Check the metric type. At present physical and static_conformal are
-    ! supported.
+!   Check the metric type. At present physical and static_conformal are
+!   supported.
     if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
 
       do k = kzl, kzr
@@ -115,10 +123,10 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
           do i = ixl, ixr
             if ( eh_mask(i,j,k,l) .ge. 0 ) then
 
-              ! calculate the square of the lapse.
+!             calculate the square of the lapse.
               alp2 = alp(i,j,k)**2
 
-              ! Calculate the inverse of the 3-metric.
+!             Calculate the inverse of the 3-metric.
               guxx = gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k)**2
               guxy = gxz(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gzz(i,j,k)
               guxz = gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k)
@@ -132,20 +140,25 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
               guxz = idetg * guxz
 
               guyy = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
-              guyz = ( gxy(i,j,k) * gxz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) ) * idetg
-              guzz = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
+              guyz = ( gxy(i,j,k) * gxz(i,j,k) - &
+                       gxx(i,j,k) * gyz(i,j,k) ) * idetg
+              guzz = ( gxx(i,j,k) * gyy(i,j,k) - &
+                       gxy(i,j,k)**2 ) * idetg
         
-              ! Raise the index of the partial derivatives of f.
-              dfux = guxx * dfx(i,j,k,l) + guxy * dfy(i,j,k,l) + guxz * dfz(i,j,k,l)
-              dfuy = guxy * dfx(i,j,k,l) + guyy * dfy(i,j,k,l) + guyz * dfz(i,j,k,l)
-              dfuz = guxz * dfx(i,j,k,l) + guyz * dfy(i,j,k,l) + guzz * dfz(i,j,k,l)
-
-              ! Calculate the overall multiplication factor.
+!             Raise the index of the partial derivatives of f.
+              dfux = guxx * dfx(i,j,k,l) + guxy * dfy(i,j,k,l) + &
+                                           guxz * dfz(i,j,k,l)
+              dfuy = guxy * dfx(i,j,k,l) + guyy * dfy(i,j,k,l) + &
+                                           guyz * dfz(i,j,k,l)
+              dfuz = guxz * dfx(i,j,k,l) + guyz * dfy(i,j,k,l) + &
+                                           guzz * dfz(i,j,k,l)
+
+!             Calculate the overall multiplication factor.
               factor = alp2 / sqrt ( alp2 * ( dfux * dfx(i,j,k,l) + &
                                               dfuy * dfy(i,j,k,l) + &
                                               dfuz * dfz(i,j,k,l) ) )
 
-              ! Finally obtain dx^i/dt.
+!             Finally obtain dx^i/dt.
               xgf(i,j,k) = - betax(i,j,k) + ssign * factor * dfux
               ygf(i,j,k) = - betay(i,j,k) + ssign * factor * dfuy
               zgf(i,j,k) = - betaz(i,j,k) + ssign * factor * dfuz
@@ -166,14 +179,14 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
             if ( eh_mask(i,j,k,l) .ge. 0 ) then
               alp2 = alp(i,j,k)**2
 
-              ! The inverse of psi^4
+!             The inverse of psi^4
               psi4 = one / psi(i,j,k)**4
 
               guxx = gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k)**2
               guxy = gxz(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gzz(i,j,k)
               guxz = gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k)
 
-!             The determinant divided by psi^4.
+!             The inverse of the determinant divided by psi^4.
               idetg = psi4 / ( gxx(i,j,k) * guxx + &
                                gxy(i,j,k) * guxy + &
                                gxz(i,j,k) * guxz )
@@ -185,12 +198,16 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
               guxz = idetg * guxz
 
               guyy = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
-              guyz = ( gxy(i,j,k) * gxz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) ) * idetg
+              guyz = ( gxy(i,j,k) * gxz(i,j,k) - &
+                       gxx(i,j,k) * gyz(i,j,k) ) * idetg
               guzz = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
       
-              dfux = guxx * dfx(i,j,k,l) + guxy * dfy(i,j,k,l) + guxz * dfz(i,j,k,l)
-              dfuy = guxy * dfx(i,j,k,l) + guyy * dfy(i,j,k,l) + guyz * dfz(i,j,k,l)
-              dfuz = guxz * dfx(i,j,k,l) + guyz * dfy(i,j,k,l) + guzz * dfz(i,j,k,l)
+              dfux = guxx * dfx(i,j,k,l) + guxy * dfy(i,j,k,l) + &
+                                           guxz * dfz(i,j,k,l)
+              dfuy = guxy * dfx(i,j,k,l) + guyy * dfy(i,j,k,l) + &
+                                           guyz * dfz(i,j,k,l)
+              dfuz = guxz * dfx(i,j,k,l) + guyz * dfy(i,j,k,l) + &
+                                           guzz * dfz(i,j,k,l)
               factor = alp2 / sqrt ( alp2 * ( dfux * dfx(i,j,k,l) + &
                                               dfuy * dfy(i,j,k,l) + &
                                               dfuz * dfz(i,j,k,l) ) )
@@ -207,7 +224,7 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS)
       end do
     end if
 
-    ! Call the interpolator.
+!   Call the interpolator.
     call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
                                  table_handle, coord_system_handle, &
                                  lsh(1), CCTK_VARIABLE_REAL, &
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
index c291376..32c50a4 100644
--- a/src/EHFinder_Init.F90
+++ b/src/EHFinder_Init.F90
@@ -23,25 +23,22 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
   CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el
   CCTK_INT, dimension(1) :: lsh, lbnd
 
-  ! Initialize the current_iteration variable.
-  current_iteration = last_iteration_number
-
-  ! Initialize the last_time varible. Note the parameters should be
-  ! chosen to be consistent with the run producing the numerical data.
-  ! F.ex. if dt was 0.1 but the data was only stored every 4 iterations, the
-  ! parameters should be chosen so that dt now is 0.4.
+! Initialize the last_time varible. Note the parameters should be
+! chosen to be consistent with the run producing the numerical data.
+! F.ex. if dt was 0.1 but the data was only stored every 4 iterations, the
+! parameters should be chosen so that dt now is 0.4.
   last_time = abs(cctk_delta_time) * last_iteration_number / &
                                      saved_iteration_every
 
   cctk_time = last_time
 
-  ! Allocate the logical array containing the flag determining if the
-  ! corresponding level set should be re-initialized.
+! Allocate the logical array containing the flag determining if the
+! corresponding level set should be re-initialized.
 
-  if ( allocated(reparam_this_level_set) ) then
-    deallocate ( reparam_this_level_set )
+  if ( allocated(re_init_this_level_set) ) then
+    deallocate ( re_init_this_level_set )
   end if
-  allocate ( reparam_this_level_set(eh_number_level_sets) )
+  allocate ( re_init_this_level_set(eh_number_level_sets) )
   if ( allocated(re_initialize_undone) ) then
     deallocate ( re_initialize_undone )
   end if
@@ -49,13 +46,13 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
   if ( evolve_generators .gt. 0 ) then
 
-    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::xg" )
+    call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, 'ehfinder::xg' )
     if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" )
+      call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' )
     end if
-    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" )
+    call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' )
     if ( status .lt. 0 ) then
-      call CCTK_WARN ( 0, "cannot get local size for generator arrays" )
+      call CCTK_WARN ( 0, 'cannot get local size for generator arrays' )
     end if
 
     if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then
@@ -89,11 +86,11 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
 
   do l = 1, eh_number_level_sets
 
-    ! If a sphere is requested...
+!   If a sphere is requested...
     if ( CCTK_EQUALS( initial_f(l), 'sphere' ) ) then
 
-      ! Set up a sphere of radius initial_rad and translated 
-      ! (translate_x,translate_y,translate_z) away from the origin.
+!     Set up a sphere of radius initial_rad and translated 
+!     (translate_x,translate_y,translate_z) away from the origin.
       f(:,:,:,l) = sqrt( ( x - translate_x(l) )**2 + &
                          ( y - translate_y(l) )**2 + &
                          ( z - translate_z(l) )**2 ) - initial_rad(l)
@@ -109,10 +106,10 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       end if
     end if
 
-    ! If an ellipsoid is requested...
+!   If an ellipsoid is requested...
     if ( CCTK_EQUALS( initial_f(l), 'ellipsoid' ) ) then
 
-      ! Calculate sines and cosines of the rotation parameters.
+!     Calculate sines and cosines of the rotation parameters.
       cosa = cos(rotation_alpha(l))
       sina = sin(rotation_alpha(l))
       cosb = cos(rotation_beta(l))
@@ -120,8 +117,8 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       cosc = cos(rotation_gamma(l))
       sinc = sin(rotation_gamma(l))
 
-      ! Set up the rotation matrix. The order is alpha around the z-axis,
-      ! beta around the y-axis and finally gamma around the x-axis.
+!     Set up the rotation matrix. The order is alpha around the z-axis,
+!     beta around the y-axis and finally gamma around the x-axis.
       txyz(1,1) = cosa * cosb
       txyz(1,2) = sina * cosb
       txyz(1,3) = -sinb
@@ -131,11 +128,10 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       txyz(3,1) = cosa * sinb * cosc + sina * sinc
       txyz(3,2) = sina * sinb * cosc - cosa * sinc
       txyz(3,3) = cosb * cosc
-!      print*,txyz
 
-    ! Apply the rotations and translation for all points on the grid.
-    ! Even though at first glance it looks like the translation is done
-    ! first, the opposite is actually true.
+!     Apply the rotations and translation for all points on the grid.
+!     Even though at first glance it looks like the translation is done
+!     first, the opposite is actually true.
       do k = 1, nz
         do j = 1, ny
           do i = 1, nx
@@ -168,14 +164,14 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS)
       end if
     end if
 
-    ! if an ovaloid of Cassini is requested...
+!   if an ovaloid of Cassini is requested...
     if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then
       f(:,:,:,l) = (x**2+y**2+z**2)**2 + cas_a(l)**4 - &
             2*cas_a(l)**2*(x**2 - (y**2+z**2)) - cas_b(l)**4
     end if
   end do
 
-  ! Initialise the internal mask.
+! Initialise the internal mask.
   eh_mask = 0
 
   return
@@ -191,26 +187,13 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  ! Set up the value used in interiour inactive cells.
+! Set up the value used in interiour inactive cells.
   ex_value = - ( one + shell_width ) * maxval(cctk_delta_space)
 
-  ! Initialize nx, ny and nz based on the local gris size. They are defined
-  ! in the module EHFinder_mod and will therefore be known from now on by all
-  ! routine that uses this module.
-  nx = cctk_lsh(1)
-  ny = cctk_lsh(2)
-  nz = cctk_lsh(3)
-
-  ! Find the maximal grid spacing.
+! Find the maximal grid spacing.
   delta = maxval ( cctk_delta_space )
 
-  ! Initialise the ghostzone information. These variables are defined in
-  ! the module EHFinder_mod.
-  ngx = cctk_nghostzones(1)
-  ngy = cctk_nghostzones(2)
-  ngz = cctk_nghostzones(3)
-
-  ! Get handles for various reduction operations.
+! Get handles for various reduction operations.
   call CCTK_ReductionArrayHandle ( max_handle, 'maximum' )
   if ( max_handle .lt. 0 ) then
     call CCTK_WARN(0,'Could not obtain a handle for maximum reduction')
@@ -224,24 +207,4 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS)
     call CCTK_WARN(0,'Could not obtain a handle for sum reduction')
   end if
 
-!  ! Register a coodinatsystem for the 2D surface grid arrays. This has no
-!  ! effect yet, but should do something useful, when TAGS can be used to
-!  ! assign coordinates to grid arrays.
-!  call CCTK_CoordRegisterSystem ( ierr, 2, 'cart2d' )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN(1,'Could not register a 2D coordinate system as "cart2d"')
-!  end if
-!  call CCTK_CoordRegisterData ( ierr, 1, 'ehfinder::ctheta', &
-!                                                 'x', 'cart2d' )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN(1,'Could not register ctheta as the 1st coordinate for "cart2d"')
-!  end if
-!  call CCTK_CoordRegisterData ( ierr, 2, 'ehfinder::cphi', &
-!                                                 'y', 'cart2d' )
-!  if ( ierr .lt. 0 ) then
-!    call CCTK_WARN(1,'Could not register cphi as the 2nd coordinate for "cart2d"')
-!  end if
-!
-!  call CCTK_INFO('2d coordinate system registered')
-
 end subroutine EHFinder_Init
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
diff --git a/src/EHFinder_ParamCheck.F90 b/src/EHFinder_ParamCheck.F90
index ae9ec1e..9640d45 100644
--- a/src/EHFinder_ParamCheck.F90
+++ b/src/EHFinder_ParamCheck.F90
@@ -12,40 +12,21 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-!  @@
-!  @routine    EHFinder_ParamCheck
-!  @date       Tue May 21 21:26:45 2002
-!  @author     Peter Diener
-!  @desc 
-!  Scheduled routine to detect invalid parameter settings.
-!  @enddesc 
-!  @calls     
-!  @calledby   
-!  @history 
-!
-!  @endhistory 
-!
-! @@
-
 subroutine EHFinder_ParamCheck(CCTK_ARGUMENTS)
  
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
+  character(len=128) :: param_warn
+
   ! Check if the metric is one of the known types: physical or
   ! static conformal
   if ( ( .not. CCTK_EQUALS(metric_type,'physical') ) .and. &
        ( .not. CCTK_EQUALS(metric_type,'static conformal') ) ) then
-    call CCTK_PARAMWARN('Unknown ADMBase::metric_type - known types are "physical" and "static conformal"')
+    param_warn = 'Unknown ADMBase::metric_type - known types are "physical" '
+    param_warn = param_warn//'and "static conformal"'
+    call CCTK_PARAMWARN ( param_warn )
   end if
 
-  ! Check if ntheta and nphi are odd as required for the Simpsons rulei
-  ! integration.
-!  if ( mod(ntheta,2) .eq. 0 ) then
-!    call CCTK_PARAMWARN('EHFinder::ntheta has to be an odd number')
-!  end if
-!  if ( mod(nphi,2) .eq. 0 ) then
-!    call CCTK_PARAMWARN('EHFinder::nphi has to be an odd number')
-!  end if
 end subroutine EHFinder_ParamCheck
diff --git a/src/EHFinder_ReadData.F90 b/src/EHFinder_ReadData.F90
index 08e3ce8..0185f5b 100644
--- a/src/EHFinder_ReadData.F90
+++ b/src/EHFinder_ReadData.F90
@@ -21,25 +21,24 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
   character(len=10) :: iteration_string
   CCTK_INT :: i, nc, ntot, res
 
-  ! Figure out which iteration number to read, based on the parameters
-  ! last_iteration_number (the last iteration in the numerical evolution
-  ! producing the metric) and saved_iteration_every (how often was the metric
-  ! saved) and the current iteration and save it in a string variable.
-
+! Figure out which iteration number to read, based on the parameters
+! last_iteration_number (the last iteration in the numerical evolution
+! producing the metric) and saved_iteration_every (how often was the metric
+! saved) and the current iteration and save it in a string variable.
   write(iteration_string,'(i10)') last_iteration_number - &
                                   saved_iteration_every * cctk_iteration
 
-  ! Trim the string variable.
+! Trim the string variable.
   iteration_string = adjustl(iteration_string)
   nc = len_trim(iteration_string)
 
-  ! Generate a string with the filenames. Note at present the requirement
-  ! therefore is that all iterations of one variable is written in the
-  ! same file.
+! Generate a string with the filenames. Note at present the requirement
+! therefore is that all iterations of one variable is written in the
+! same file.
   in_files = 'gxx gxy gxz gyy gyz gzz'
 
-  ! Fill in the string array, used to requesting the metric components at
-  ! the specified iteration number.
+! Fill in the string array, used to requesting the metric components at
+! the specified iteration number.
   var_names(1) = 'admbase::gxx[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(2) = 'admbase::gxy[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(3) = 'admbase::gxz[cctk_iteration='//iteration_string(1:nc)//']'
@@ -47,7 +46,7 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
   var_names(5) = 'admbase::gyz[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(6) = 'admbase::gzz[cctk_iteration='//iteration_string(1:nc)//']'
 
-  ! merge all the varible names into a single string.
+! merge all the variable names into a single string.
   in_vars = ' '
   ntot = 0
   do i = 1, 6
@@ -56,8 +55,8 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS)
     ntot = ntot + nc + 1
   end do
    
-  ! Call the routine that actualle does the file access. Note failures are
-  ! just silently ignored.
+! Call the routine that actualle does the file access. Note failures are
+! just silently ignored.
   call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
 
 end subroutine EHFinder_Read_Metric
@@ -78,14 +77,28 @@ subroutine EHFinder_Read_Lapse(CCTK_ARGUMENTS)
   character(len=10) :: iteration_string
   CCTK_INT :: nc, res
 
+! Figure out which iteration number to read, based on the parameters
+! last_iteration_number (the last iteration in the numerical evolution
+! producing the metric) and saved_iteration_every (how often was the lapse
+! saved) and the current iteration and save it in a string variable.
   write(iteration_string,'(i10)') last_iteration_number - &
                                   saved_iteration_every * cctk_iteration
+
+! Trim the string variable.
   iteration_string = adjustl(iteration_string)
   nc = len_trim(iteration_string)
 
+! Generate a string with the filename. Note at present the requirement
+! therefore is that all iterations of one variable is written in the
+! same file.
   in_files = 'alp'
+
+! Fill in the string used to requesting the lapse at the specified
+! iteration number.
   in_vars = 'admbase::alp[cctk_iteration='//iteration_string(1:nc)//']'
 
+! Call the routine that actualle does the file access. Note failures are
+! just silently ignored.
   call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
 
 end subroutine EHFinder_Read_Lapse
@@ -107,15 +120,28 @@ subroutine EHFinder_Read_Shift(CCTK_ARGUMENTS)
   character(len=10) :: iteration_string
   CCTK_INT :: i, nc, ntot, res
 
+! Figure out which iteration number to read, based on the parameters
+! last_iteration_number (the last iteration in the numerical evolution
+! producing the metric) and saved_iteration_every (how often was the metric
+! saved) and the current iteration and save it in a string variable.
   write(iteration_string,'(i10)') last_iteration_number - &
                                   saved_iteration_every * cctk_iteration
+! Trim the string variable.
   iteration_string = adjustl(iteration_string)
   nc = len_trim(iteration_string)
 
+! Generate a string with the filenames. Note at present the requirement
+! therefore is that all iterations of one variable is written in the
+! same file.
   in_files = 'betax betay betaz'
+
+! Fill in the string array, used to requesting the shift components at
+! the specified iteration number.
   var_names(1) = 'admbase::betax[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(2) = 'admbase::betay[cctk_iteration='//iteration_string(1:nc)//']'
   var_names(3) = 'admbase::betaz[cctk_iteration='//iteration_string(1:nc)//']'
+
+! merge all the variable names into a single string.
   in_vars = ' '
   ntot = 0
   do i = 1, 3
@@ -124,6 +150,8 @@ subroutine EHFinder_Read_Shift(CCTK_ARGUMENTS)
     ntot = ntot + nc + 1
   end do
  
+! Call the routine that actualle does the file access. Note failures are
+! just silently ignored.
   call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
 
 end subroutine EHFinder_Read_Shift
@@ -144,18 +172,34 @@ subroutine EHFinder_Read_Conformal(CCTK_ARGUMENTS)
   character(len=10) :: iteration_string
   CCTK_INT :: nc, res
 
+! Figure out which iteration number to read, based on the parameters
+! last_iteration_number (the last iteration in the numerical evolution
+! producing the metric) and saved_iteration_every (how often was the metric
+! saved) and the current iteration and save it in a string variable. Note
+! that for this variable the default option is to only to read it in once at
+! the first timestep.
   if ( read_conformal_factor_once .gt. 0 ) then
     write(iteration_string,'(i10)') 0
   else
     write(iteration_string,'(i10)') last_iteration_number - &
                                     saved_iteration_every * cctk_iteration
   end if
+
+! Trim the string variable.
   iteration_string = adjustl(iteration_string)
   nc = len_trim(iteration_string)
 
+! Generate a string with the filenames. Note at present the requirement
+! therefore is that all iterations of one variable is written in the
+! same file.
   in_files = 'psi'
+
+! Fill in the string array, used to requesting the conformal factor at
+! the specified iteration number.
   in_vars = 'staticconformal::psi[cctk_iteration='//iteration_string(1:nc)//']'
 
+! Call the routine that actualle does the file access. Note failures are
+! just silently ignored.
   call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
 
 end subroutine EHFinder_Read_Conformal
@@ -176,66 +220,28 @@ subroutine EHFinder_Read_Mask(CCTK_ARGUMENTS)
   character(len=10) :: iteration_string
   CCTK_INT :: nc, res
 
+! Figure out which iteration number to read, based on the parameters
+! last_iteration_number (the last iteration in the numerical evolution
+! producing the metric) and saved_iteration_every (how often was the metric
+! saved) and the current iteration and save it in a string variable.
   write(iteration_string,'(i10)') last_iteration_number - &
                                   saved_iteration_every * cctk_iteration
+
+! Trim the string variable.
   iteration_string = adjustl(iteration_string)
   nc = len_trim(iteration_string)
 
+! Generate a string with the filename. Note at present the requirement
+! therefore is that all iterations of one variable is written in the
+! same file.
   in_files = 'emask'
+
+! Fill in the string array, used to requesting the old style mask at
+! the specified iteration number.
   in_vars = 'spacemask::emask[cctk_iteration='//iteration_string(1:nc)//']'
 
+! Call the routine that actualle does the file access. Note failures are
+! just silently ignored.
   call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
 
 end subroutine EHFinder_Read_Mask
-
-
-! This routine reads in everything.
-subroutine EHFinder_ReadData(CCTK_ARGUMENTS)
-
-  use EHFinder_mod
-
-  implicit none
-
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_FUNCTIONS
-
-  character(len=1024) :: in_files, in_vars
-  character(len=64), dimension(11)  :: var_names
-  character(len=10) :: iteration_string
-  CCTK_INT :: i, nc, ntot, res
-
-  write(iteration_string,'(i10)') last_iteration_number - &
-                                  saved_iteration_every * cctk_iteration
-  iteration_string = adjustl(iteration_string)
-  nc = len_trim(iteration_string)
-  
-  in_files = 'alp betax betay betaz gxx gxy gxz gyy gyz gzz psi'
-  var_names(1) = 'admbase::alp[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(2) = 'admbase::betax[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(3) = 'admbase::betay[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(4) = 'admbase::betaz[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(5) = 'admbase::gxx[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(6) = 'admbase::gxy[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(7) = 'admbase::gxz[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(8) = 'admbase::gyy[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(9) = 'admbase::gyz[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(10) = 'admbase::gzz[cctk_iteration='//iteration_string(1:nc)//']'
-  var_names(11) = 'staticconformal::psi[cctk_iteration='//iteration_string(1:nc)//']'
- 
-  in_vars = ' '
-  ntot = 0
-  do i = 1, 11
-    nc = len_trim(var_names(i))
-    in_vars(ntot+1:ntot+1+nc+1) = var_names(i)(1:nc+1)
-    ntot = ntot + nc + 1
-  end do
-
-  print*,in_vars(1:len_trim(in_vars)),len_trim(in_vars)
-
-  call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars )
-  print*,res
-  current_iteration = current_iteration - 1
-  
-  return
-end subroutine EHFinder_ReadData
diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90
index 566f654..eb13858 100644
--- a/src/EHFinder_SetMask.F90
+++ b/src/EHFinder_SetMask.F90
@@ -65,9 +65,8 @@ end subroutine EHFinder_MaskInit
 
 
 ! This routine will excise points and re-activate excised points
-! after need. EHFinder_Setmask2 will then find the boundary of the
+! as needed. EHFinder_Setmask2 will then find the boundary of the
 ! excsion region and make sure that the mask is correct there.
-
 subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -218,9 +217,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
             call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' )
           end if
 
-!          print*, imin_glob(1), imax_glob(1), fimin_glob(1), fimax_glob(1)
-!          print*, imin_glob(2), imax_glob(2), fimin_glob(2), fimax_glob(2)
-!          print*, imin_glob(3), imax_glob(3), fimin_glob(3), fimax_glob(3)
         end if
            
 !       Now check and see if any interior excised points need to be activated.
@@ -307,11 +303,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
             end if
           end do
 
-!          print*,'New range'
-!          print*,imin_n(1),imax_n(1)
-!          print*,imin_n(2),imax_n(2)
-!          print*,imin_n(3),imax_n(3)
-
 !         Use the new indices to actually activate points, taking care to
 !         activate points that are actually on the local grid.  First do the
 !         faces of the rectangular box.
@@ -322,11 +313,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
               j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr )
               k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl )
               k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr )
-!              print*,'x1 ',i1,j1,j2,k1,k2
               if ( ( j1 .le. j2 ) .and. ( k1 .le. k2 ) ) then
                 tm_mask(i1,j1:j2,k1:k2,l) = 0
                 ftmp(i1,j1:j2,k1:k2,l) = ftmp(i1+1,j1:j2,k1:k2,l) + delta
-!                print*,'done'
               end if
             end if
           end if
@@ -337,11 +326,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
               j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr )
               k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl )
               k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr )
-!              print*,'x2 ',i2,j1,j2,k1,k2
               if ( ( j1 .le. j2 ) .and. ( k1 .le. k2 ) ) then
                 tm_mask(i2,j1:j2,k1:k2,l) = 0
                 ftmp(i2,j1:j2,k1:k2,l) = ftmp(i2-1,j1:j2,k1:k2,l) + delta
-!                print*,'done'
               end if
             end if
           end if
@@ -352,11 +339,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
               i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr )
               k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl )
               k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr )
-!              print*,'y1 ',j1,i1,i2,k1,k2
               if ( ( i1 .le. i2 ) .and. ( k1 .le. k2 ) ) then
                 tm_mask(i1:i2,j1,k1:k2,l) = 0
                 ftmp(i1:i2,j1,k1:k2,l) = ftmp(i1:i2,j1+1,k1:k2,l) + delta
-!                print*,'done'
               end if
             end if
           end if
@@ -367,11 +352,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
               i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr )
               k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl )
               k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr )
-!              print*,'y2 ',j2,i1,i2,k1,k2
               if ( ( i1 .le. i2 ) .and. ( k1 .le. k2 ) ) then
                 tm_mask(i1:i2,j2,k1:k2,l) = 0
                 ftmp(i1:i2,j2,k1:k2,l) = f(i1:i2,j2-1,k1:k2,l) + delta
-!                print*,'done'
               end if
             end if
           end if
@@ -382,11 +365,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
               i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr )
               j1 = max ( max ( imin_n(2), imin_glob(2) ) - cctk_lbnd(2), jyl )
               j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr )
-!              print*,'z1 ',k1,i1,i2,j1,j2
               if ( ( i1 .le. i2 ) .and. ( j1 .le. j2 ) ) then
                 tm_mask(i1:i2,j1:j2,k1,l) = 0
                 ftmp(i1:i2,j1:j2,k1,l) = ftmp(i1:i2,j1:j2,k1+1,l) + delta
-!                print*,'done'
               end if
             end if
           end if
@@ -397,11 +378,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
               i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr )
               j1 = max ( max ( imin_n(2), imin_glob(2) ) - cctk_lbnd(2), jyl )
               j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr )
-!              print*,'z2 ',k2,i1,i2,j1,j2
               if ( ( i1 .le. i2 ) .and. ( j1 .le. j2 ) ) then
                 tm_mask(i1:i2,j1:j2,k2,l) = 0
                 ftmp(i1:i2,j1:j2,k2,l) = ftmp(i1:i2,j1:j2,k2-1,l) + delta
-!                print*,'done'
               end if
             end if
           end if
@@ -674,7 +653,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
             i2 = imax_n(1) - cctk_lbnd(1)
             j2 = imax_n(2) - cctk_lbnd(2)
             k2 = imax_n(3) - cctk_lbnd(3)
-!            print*,'Debug 3 : ',i2,j2,k2
             if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. &
                  ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. &
                  ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then
@@ -829,7 +807,6 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
                         tm_mask(ixl:ixr,jyl:jyr,kzl:kzr,l)
         end where
       end if
-!      call CCTK_INFO('Mask modified')
     end do
 
 !   Indicate that it is not anymore the first time the mask is set.
@@ -839,6 +816,10 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
 end subroutine EHFinder_SetMask2
 
 
+! This routine checks if any of the excision regions are to close together
+! so that some points do not have enough active neigbouring points to be
+! able to calculate derivatives. If this is the case these points are 
+! excised as well.
 subroutine EHFinder_SetMask3(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -884,81 +865,76 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS)
           do j = jyl, jyr
             do i = ixl, ixr
               if ( eh_mask(i,j,k,l) .ge. 0 ) then
-                ! If we are on an excision boundary in the -x-direction...
+
+!               If we are on an excision boundary in the -x-direction...
                 if ( btest(eh_mask(i,j,k,l),0) ) then
-                  ! If any of the two closest neighbours in the +x-direction is
-                  ! excised we have to excise this point
+
+!                 If any of the two closest neighbours in the +x-direction is
+!                 excised we have to excise this point
                   if ( ( eh_mask(i+1,j,k,l) .eq. -1 ) .or. &
                        ( eh_mask(i+2,j,k,l) .eq. -1 ) ) then
                     tm_mask(i,j,k,l) = -1
                     mask_modified = .true.
-                    print*, i, j, k, 'x-'
-                    print*, eh_mask(i:i+2,j,k,l)
                   end if
                 end if
 
-                ! If we are on an excision boundary in the +x-direction...
+!               If we are on an excision boundary in the +x-direction...
                 if ( btest(eh_mask(i,j,k,l),1) ) then
-                  ! If any of the two closest neighbours in the -x-direction is
-                  ! excised we have to excise this point
+
+!                 If any of the two closest neighbours in the -x-direction is
+!                 excised we have to excise this point
                   if ( ( eh_mask(i-1,j,k,l) .eq. -1 ) .or. &
                        ( eh_mask(i-2,j,k,l) .eq. -1 ) ) then
                     tm_mask(i,j,k,l) = -1
                     mask_modified = .true.
-                    print*, i, j, k, 'x+'
-                    print*, eh_mask(i:i-2,j,k,l)
                   end if
                 end if
 
-                ! If we are on an excision boundary in the -y-direction...
+!               If we are on an excision boundary in the -y-direction...
                 if ( btest(eh_mask(i,j,k,l),2) ) then
-                  ! If any of the two closest neighbours in the +y-direction is
-                  ! excised we have to excise this point
+
+!                 If any of the two closest neighbours in the +y-direction is
+!                 excised we have to excise this point
                   if ( ( eh_mask(i,j+1,k,l) .eq. -1 ) .or. &
                        ( eh_mask(i,j+2,k,l) .eq. -1 ) ) then
                     tm_mask(i,j,k,l) = -1
                     mask_modified = .true.
-                    print*, i, j, k, 'y-'
-                    print*, eh_mask(i,j:j+2,k,l)
                   end if
                 end if
 
-                ! If we are on an excision boundary in the +y-direction...
+!               If we are on an excision boundary in the +y-direction...
                 if ( btest(eh_mask(i,j,k,l),3) ) then
-                  ! If any of the two closest neighbours in the -y-direction is
-                  ! excised we have to excise this point
+
+!                 If any of the two closest neighbours in the -y-direction is
+!                 excised we have to excise this point
                   if ( ( eh_mask(i,j-1,k,l) .eq. -1 ) .or. &
                        ( eh_mask(i,j-2,k,l) .eq. -1 ) ) then
                     tm_mask(i,j,k,l) = -1
                     mask_modified = .true.
-                    print*, i, j, k, 'y+'
-                    print*, eh_mask(i,j:j-2,k,l)
                   end if
                 end if
 
-                ! If we are on an excision boundary in the -z-direction...
+!               If we are on an excision boundary in the -z-direction...
                 if ( btest(eh_mask(i,j,k,l),4) ) then
-                  ! If any of the two closest neighbours in the +z-direction is
-                  ! excised we have to excise this point
+
+!                 If any of the two closest neighbours in the +z-direction is
+!                 excised we have to excise this point
                   if ( ( eh_mask(i,j,k+1,l) .eq. -1 ) .or. &
                        ( eh_mask(i,j,k+2,l) .eq. -1 ) ) then
                     tm_mask(i,j,k,l) = -1
                     mask_modified = .true.
-                    print*, i, j, k, 'z-'
-                    print*, eh_mask(i,j,k:k+2,l)
                   end if
                 end if
 
-                ! If we are on an excision boundary in the +z-direction...
+!               If we are on an excision boundary in the +z-direction...
                 if ( btest(eh_mask(i,j,k,l),5) ) then
-                  ! If any of the two closest neighbours in the -z-direction is
-                  ! excised we have to excise this point
+
+!                 If any of the two closest neighbours in the -z-direction is
+!                 excised we have to excise this point
                   if ( ( eh_mask(i,j,k-1,l) .eq. -1 ) .or. &
                        ( eh_mask(i,j,k-2,l) .eq. -1 ) ) then
                     tm_mask(i,j,k,l) = -1
                     mask_modified = .true.
-                    print*, i, j, k, 'z+'
-                    print*, eh_mask(i,j,k:k-2,l)
                   end if
                 end if
               end if
@@ -968,7 +944,7 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS)
         end do
 
         if ( mask_modified ) then
-          call CCTK_WARN ( 1, "Mask modified because it was not convex" )
+          call CCTK_WARN ( 1, 'Mask modified because it was not convex' )
         end if
       end if
     end do
diff --git a/src/EHFinder_SetSym.F90 b/src/EHFinder_SetSym.F90
index b04672d..7811406 100644
--- a/src/EHFinder_SetSym.F90
+++ b/src/EHFinder_SetSym.F90
@@ -1,4 +1,4 @@
-! Initialisation of the level set function
+! Registration of symmetries for the necessary grid functions.
 ! $Header$
 
 #include "cctk.h"
@@ -15,33 +15,28 @@ subroutine EHFinder_SetSym(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
 
   CCTK_INT, dimension(3) :: sym
-  CCTK_INT :: l
-  character(len=14) :: fname
 
+! All grid functions have even symmetries in all directions.
   sym = 1
 
-  do l = 0, eh_number_level_sets - 1
-    write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']'
-    call SetCartSymVN(ierr,cctkGH,sym,fname)
-    if ( ierr .gt. 0 ) then
-      call CCTK_WARN(1,'Failed to register symmetry for the level set function')
-    end if
-!    call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::f[2]')
-!    if ( ierr .gt. 0 ) then
-!      call CCTK_WARN(1,'Failed to register symmetry for the level set function')
-!    end if
-  end do
+! Set up symmetries for the level set function.
+  call SetCartSymGN ( ierr, cctkGH, sym, 'ehfinder::f' )
+  if ( ierr .gt. 0 ) then
+    call CCTK_WARN(1,'Failed to register symmetry for the level set function')
+  end if
+
+! Set up symmetries for the mask function.
+  call SetCartSymGN( ierr, cctkGH, sym, 'ehfinder::eh_mask' )
+  if ( ierr .gt. 0 ) then
+    call CCTK_WARN(1,'Failed to register symmetry for eh_mask')
+  end if
 
+! Set up symmetries for the surface index function.
   call SetCartSymGN(ierr,cctkGH,sym,'ehfinder::surface_index')
   if ( ierr .gt. 0 ) then
     call CCTK_WARN(1,'Failed to register symmetry for sc')
   end if
 
-!  call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::rep_mask')
-!  if ( ierr .gt. 0 ) then
-!    call CCTK_WARN(1,'Failed to register symmetry for the level reparametrization mask')
-!  end if
-
   return
 end subroutine EHFinder_SetSym
 
@@ -56,28 +51,22 @@ subroutine EHFinder_ApplySymAll(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: l
-  character(len=14) :: fname
-
-  do l = 0, eh_number_level_sets - 1
-    write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']'
-    call CartSymVN(ierr,cctkGH,fname)
-  end do
+  character(len=80) :: warn_message
 
-# include "include/physical_part.h"
-
-  if ( symx ) then
-    eh_mask(ngx:1:-1,:,:,:) = eh_mask(ngx+1:ngx+ngx,:,:,:)
-!    rep_mask(ngx:1:-1,:,:,:) = rep_mask(ngx+1:ngx+ngx,:,:,:)
-  end if
-  if ( symy ) then
-    eh_mask(:,ngy:1:-1,:,:) = eh_mask(:,ngy+1:ngy+ngy,:,:)
-!    rep_mask(:,ngy:1:-1,:,:) = rep_mask(:,ngy+1:ngy+ngy,:,:)
+! Apply symmetry for the level set function.
+  call CartSymGN ( ierr, cctkGH, 'ehfinder::f' )
+  if ( ierr .gt. 0 ) then
+    warn_message = 'Failed to perform symmetry operation on the level '
+    warn_message = warn_message//'set function'
+    call CCTK_WARN( 1, warn_message )
   end if
-  if ( symz ) then
-    eh_mask(:,:,ngz:1:-1,:) = eh_mask(:,:,ngz+1:ngz+ngz,:)
-!    rep_mask(:,:,ngz:1:-1,:) = rep_mask(:,:,ngz+1:ngz+ngz,:)
+
+! Apply symmetry for the mask function.
+  call CartSymGN ( ierr, cctkGH, 'ehfinder::eh_mask' )
+  if ( ierr .gt. 0 ) then
+    call CCTK_WARN( 1, 'Failed to perform symmetry operation on the mask')
   end if
+
   return
 end subroutine EHFinder_ApplySymAll
 
@@ -92,19 +81,21 @@ subroutine EHFinder_ApplySymF(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: l
-  character(len=14) :: fname
+  character(len=80) :: warn_message
 
-  do l = 0, eh_number_level_sets - 1
-    write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']'
-    call CartSymVN(ierr,cctkGH,fname)
-  end do
+! Apply symmetry for the level set function.
+  call CartSymGN ( ierr, cctkGH, 'ehfinder::f' )
+  if ( ierr .gt. 0 ) then
+    warn_message = 'Failed to perform symmetry operation on the level '
+    warn_message = warn_message//'set function'
+    call CCTK_WARN( 1, warn_message )
+  end if
 
   return
 end subroutine EHFinder_ApplySymF
 
 
-subroutine EHFinder_ApplySymRep(CCTK_ARGUMENTS)
+subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS)
 
   use EHFinder_mod
 
@@ -114,55 +105,17 @@ subroutine EHFinder_ApplySymRep(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
-# include "include/physical_part.h"
-
-!  if ( symx ) then
-!    rep_mask(ngx:1:-1,:,:,:) = rep_mask(ngx+1:ngx+ngx,:,:,:)
-!  end if
-!  if ( symy ) then
-!    rep_mask(:,ngy:1:-1,:,:) = rep_mask(:,ngy+1:ngy+ngy,:,:)
-!  end if
-!  if ( symz ) then
-!    rep_mask(:,:,ngz:1:-1,:) = rep_mask(:,:,ngz+1:ngz+ngz,:)
-!  end if
-  return
-end subroutine EHFinder_ApplySymRep
-
-
-subroutine EHFinder_ApplySymFRep(CCTK_ARGUMENTS)
-
-  use EHFinder_mod
-
-  implicit none
-
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_FUNCTIONS
+! Apply symmetry for the mask function.
+  call CartSymGN ( ierr, cctkGH, 'ehfinder::eh_mask' )
+  if ( ierr .gt. 0 ) then
+    call CCTK_WARN( 1, 'Failed to perform symmetry operation on the mask')
+  end if
 
-  CCTK_INT :: l
-  character(len=14) :: fname
-
-# include "include/physical_part.h"
-
-  do l = 0, eh_number_level_sets - 1
-    write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']'
-    call CartSymVN(ierr,cctkGH,fname)
-  end do
-
-!  if ( symx ) then
-!    rep_mask(ngx:1:-1,:,:,:) = rep_mask(ngx+1:ngx+ngx,:,:,:)
-!  end if
-!  if ( symy ) then
-!    rep_mask(:,ngy:1:-1,:,:) = rep_mask(:,ngy+1:ngy+ngy,:,:)
-!  end if
-!  if ( symz ) then
-!    rep_mask(:,:,ngz:1:-1,:) = rep_mask(:,:,ngz+1:ngz+ngz,:)
-!  end if
   return
-end subroutine EHFinder_ApplySymFRep
+end subroutine EHFinder_ApplySymMask
 
 
-subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS)
+subroutine EHFinder_ApplySymSC(CCTK_ARGUMENTS)
 
   use EHFinder_mod
 
@@ -172,16 +125,14 @@ subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
 
-# include "include/physical_part.h"
+  character(len=80) :: warn_message
 
-  if ( symx ) then
-    eh_mask(ngx:1:-1,:,:,:) = eh_mask(ngx+1:ngx+ngx,:,:,:)
-  end if
-  if ( symy ) then
-    eh_mask(:,ngy:1:-1,:,:) = eh_mask(:,ngy+1:ngy+ngy,:,:)
-  end if
-  if ( symz ) then
-    eh_mask(:,:,ngz:1:-1,:) = eh_mask(:,:,ngz+1:ngz+ngz,:)
+! Apply the symmetry for the surface index function.
+  call CartSymGN ( ierr, cctkGH, 'ehfinder::surface_index' )
+  if ( ierr .gt. 0 ) then
+    warn_message = 'Failed to perform symmetry operation on the surface index'
+    call CCTK_WARN( 1, warn_message )
   end if
+
   return
-end subroutine EHFinder_ApplySymMask
+end subroutine EHFinder_ApplySymSC
diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90
index 65e0488..69a5464 100644
--- a/src/EHFinder_Sources.F90
+++ b/src/EHFinder_Sources.F90
@@ -27,37 +27,44 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
   CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
   CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
 
-! calculate 1/(2*delta) in each direction
-
 #include "include/physical_part.h"
 
+! calculate 1/(2*delta) in each direction
   idx = half / cctk_delta_space(1)
   idy = half / cctk_delta_space(2)
   idz = half / cctk_delta_space(3)
 
   mdelta = maxval ( cctk_delta_space )
 
+! Set the sign depending on the surface direction.
   if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one
   if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one
 
-  if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then
-    do l = 1, eh_number_level_sets
+  do l = 1, eh_number_level_sets
+
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+!         Calculate the inverse of the 3-metric
+# include "include/metric.h"
+        end do
+      end do
+    end do
+
+!   Calculate the derivatives of the level set function using the intrinsic
+!   scheme. Note, this should never be used and may disappear in later versions.
+    if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then
       do k = kzl, kzr
         do j = jyl, jyr
           do i = ixl, ixr
-            if ( f(i,j,k,l) .gt. fmin_bound - three*mdelta ) then
-#include "include/centered_second2.h"
-            else
 #include "include/upwind_second2.h"
-            end if
           end do
         end do
       end do
-    end do
-  end if
+    end if
     
-  if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then
-    do l = 1, eh_number_level_sets
+!   Calculate the derivatives of the level set function using shift upwinding.
+    if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then
       do k = kzl, kzr
         do j = jyl, jyr
           do i = ixl, ixr
@@ -65,119 +72,89 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
           end do
         end do
       end do
-    end do
-  end if
+    end if
 
-  if ( CCTK_EQUALS ( upwind_type, 'characteristic' ) ) then
-    do l = 1, eh_number_level_sets
+!   If the three metric is the static conformal metric we convert the inverse
+!   three metric to the physical inverse three metric by multiplying with
+!   psi^(-4).
+    if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
       do k = kzl, kzr
-        do j = jyl, jyr
+        do j = jyl,jyr
           do i = ixl, ixr
-#include "include/metric.h"
-#include "include/upwind_characteristic_second2.h"
+            if ( eh_mask(i,j,k,l) .ge. 0 ) then
+              psito4 = psi(i,j,k)**(-4)
+              g3xx(i,j,k) = g3xx(i,j,k) * psito4
+              g3xy(i,j,k) = g3xy(i,j,k) * psito4
+              g3xz(i,j,k) = g3xz(i,j,k) * psito4
+              g3yy(i,j,k) = g3yy(i,j,k) * psito4
+              g3yz(i,j,k) = g3yz(i,j,k) * psito4
+              g3zz(i,j,k) = g3zz(i,j,k) * psito4
+            end if
           end do
         end do
       end do
-    end do
-  end if
+    end if
 
-  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
-    do l = 1, eh_number_level_sets
+!   Calculate the derivatives of the level set using characteristic upwinding.
+    if ( CCTK_EQUALS ( upwind_type, 'characteristic' ) ) then
       do k = kzl, kzr
         do j = jyl, jyr
           do i = ixl, ixr
-            if ( eh_mask(i,j,k,l) .ge. 0 ) then
-              alp2 = alp(i,j,k)**2
-  
-              tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2
-              tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k)
-              tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k)
-   
-              idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 )
-  
-              g3xx = tmp1 * idetg
-              g3xy = tmp2 * idetg
-              g3xz = tmp3 * idetg
-              g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
-              g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg
-              g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
-
-              tmp1 = betax(i,j,k) * dfx(i,j,k,l) + &
-                     betay(i,j,k) * dfy(i,j,k,l) + &
-                     betaz(i,j,k) * dfz(i,j,k,l)
-
-              tmp2 = g3xx * dfx(i,j,k,l)**2 + &
-                     g3yy * dfy(i,j,k,l)**2 + &
-                     g3zz * dfz(i,j,k,l)**2 + &
-                     two * ( g3xy * dfx(i,j,k,l) * dfy(i,j,k,l) + &
-                             g3xz * dfx(i,j,k,l) * dfz(i,j,k,l) + &
-                             g3yz * dfy(i,j,k,l) * dfz(i,j,k,l) )
-              if ( tmp2 .ge. zero ) then
-                sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 )
-              else
-                call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" )
-              end if
-            else
-              sf(i,j,k,l) = zero
-            end if
+
+!           We use centered derivatives to figure out which direction
+!           to upwind in.
+#include "include/upwind_characteristic_second2.h"
           end do
         end do
-      end do 
-    end do
-  end if
+      end do
+    end if
 
-  if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then
-    do l = 1, eh_number_level_sets
-      do k = kzl, kzr
-        do j = jyl, jyr
-          do i = ixl, ixr
-            if ( eh_mask(i,j,k,l) .ge. 0 ) then
-              alp2 = alp(i,j,k)**2
- 
-              psito4 = psi(i,j,k)**4
-              gxxc = gxx(i,j,k) * psito4
-              gxyc = gxy(i,j,k) * psito4
-              gxzc = gxz(i,j,k) * psito4
-              gyyc = gyy(i,j,k) * psito4
-              gyzc = gyz(i,j,k) * psito4
-              gzzc = gzz(i,j,k) * psito4
-            
-              tmp1 = gyyc*gzzc - gyzc**2
-              tmp2 = gxzc*gyzc - gxyc*gzzc
-              tmp3 = gxyc*gyzc - gxzc*gyyc
- 
-              idetg = one / ( gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 )
-    
-              g3xx = tmp1 * idetg
-              g3xy = tmp2 * idetg
-              g3xz = tmp3 * idetg
-              g3yy = ( gxxc*gzzc - gxzc**2 ) * idetg
-              g3yz = ( gxyc*gxzc - gxxc*gyzc ) * idetg
-              g3zz = ( gxxc*gyyc - gxyc**2 ) * idetg
-
-              tmp1 = betax(i,j,k) * dfx(i,j,k,l) + &
-                     betay(i,j,k) * dfy(i,j,k,l) + &
-                     betaz(i,j,k) * dfz(i,j,k,l)
-
-              tmp2 = g3xx * dfx(i,j,k,l)**2 + &
-                     g3yy * dfy(i,j,k,l)**2 + &
-                     g3zz * dfz(i,j,k,l)**2 + &
-                     two * ( g3xy * dfx(i,j,k,l) * dfy(i,j,k,l) + &
-                             g3xz * dfx(i,j,k,l) * dfz(i,j,k,l) + &
-                             g3yz * dfy(i,j,k,l) * dfz(i,j,k,l) )
-              if ( tmp2 .ge. zero ) then
-                sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 )
-              else
-                call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" )
-              end if
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+
+!         If the current point is active ...
+          if ( eh_mask(i,j,k,l) .ge. 0 ) then
+
+!           Square the lapse.
+            alp2 = alp(i,j,k)**2
+  
+!           Calculate beta^i df_i.
+            tmp1 = betax(i,j,k) * dfx(i,j,k,l) + &
+                   betay(i,j,k) * dfy(i,j,k,l) + &
+                   betaz(i,j,k) * dfz(i,j,k,l)
+
+!           Calculate gamma^ij df_i df_j.
+            tmp2 = g3xx(i,j,k) * dfx(i,j,k,l)**2 + &
+                   g3yy(i,j,k) * dfy(i,j,k,l)**2 + &
+                   g3zz(i,j,k) * dfz(i,j,k,l)**2 + &
+                   two * ( g3xy(i,j,k) * dfx(i,j,k,l) * dfy(i,j,k,l) + &
+                           g3xz(i,j,k) * dfx(i,j,k,l) * dfz(i,j,k,l) + &
+                           g3yz(i,j,k) * dfy(i,j,k,l) * dfz(i,j,k,l) )
+
+!           If the metric is positive definite ...
+            if ( tmp2 .ge. zero ) then
+
+!             Calculate the right hand side.
+              sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 )
+
+!             If the lapse is negative we change the sign of the right hand
+!             side function. This is done to be able to handle for example
+!             Schwarzschild in isotropic coordinates with the isotropic
+!             lapse.          
+              sf(i,j,k,l) = sf(i,j,k,l) * sign ( one, alp(i,j,k) )
             else
-              sf(i,j,k,l) = zero
+
+!             Otherwise print a level 0 warning.
+              call CCTK_WARN ( 0, '3-metric not positive definite: Stopping' )
             end if
-          end do
+          else
+            sf(i,j,k,l) = zero
+          end if
         end do
-      end do 
-    end do
-  end if
- 
+      end do
+    end do 
+  end do
+
   return
 end subroutine EHFinder_Sources
diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90
index 9c731a0..3fe10a5 100644
--- a/src/EHFinder_mod.F90
+++ b/src/EHFinder_mod.F90
@@ -6,33 +6,21 @@ module EHFinder_mod
 
   use EHFinder_Constants
 
-  CCTK_INT, dimension(0:5), parameter :: ix = (/ -1, 1, 0, 0, 0, 0 /)
-  CCTK_INT, dimension(0:5), parameter :: jy = (/ 0, 0, -1, 1, 0, 0 /)
-  CCTK_INT, dimension(0:5), parameter :: kz = (/ 0, 0, 0, 0, -1, 1 /)
   CCTK_INT, dimension(0:5), parameter :: ll = (/ 1, 2, 4, 8, 16, 32 /)
   CCTK_INT :: ixl, ixr, jyl, jyr, kzl, kzr
-!  CCTK_INT, parameter :: ixl = 1
-!  CCTK_INT, parameter :: ixr = 2
-!  CCTK_INT, parameter :: iyl = 4
-!  CCTK_INT, parameter :: iyr = 8
-!  CCTK_INT, parameter :: izl = 16
-!  CCTK_INT, parameter :: izr = 32
   CCTK_REAL :: hfac
-  CCTK_INT :: nx, ny, nz, ngx, ngy, ngz, ierr
+  CCTK_INT :: nx, ny, nz, ngx, ngy, ngz
   CCTK_INT :: max_handle, min_handle, sum_handle
-  CCTK_INT :: ncalls, rep_count, rep_total
+  CCTK_INT :: ncalls, ierr
   CCTK_REAL :: delta
-  CCTK_REAL :: h, hlocal
+  CCTK_REAL :: h
   CCTK_REAL :: ex_value
-  CCTK_REAL :: rdel = 0.8d0
-  CCTK_REAL :: fmin_bound
   CCTK_INT, dimension(3) :: imin_glob, imax_glob
   CCTK_REAL, dimension(3) :: fimin_glob, fimax_glob
-  CCTK_INT :: current_iteration
   CCTK_INT :: githeta, gjphi
   logical :: mask_first = .true.
   logical :: symx, symy, symz
   logical :: s_symx, s_symy, s_symz
-  logical, dimension(:), allocatable :: reparam_this_level_set, &
+  logical, dimension(:), allocatable :: re_init_this_level_set, &
                                         re_initialize_undone
 end module EHFinder_mod
diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90
index 4683c09..58a5611 100644
--- a/src/MoL_Init.F90
+++ b/src/MoL_Init.F90
@@ -10,23 +10,19 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
   implicit none
 
   CCTK_INT :: ierr, ierr_cum, varindex, rhsindex
-  CCTK_INT :: MoLRegisterEvolved, MoLRegisterEvolvedGroup
+  CCTK_INT :: MoLRegisterEvolvedGroup!, MoLRegisterEvolved
   CCTK_INT :: l
 
-  character(len=15) :: vname
-
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
   ierr_cum = 0
-  do l = 0, eh_number_level_sets - 1
-    write(vname,'(a12,i1,a1)') 'ehfinder::f[', l, ']'
-    call CCTK_VarIndex(varindex, vname )
-    write(vname,'(a13,i1,a1)') 'ehfinder::sf[', l, ']'
-    call CCTK_VarIndex(rhsindex, vname)
-    ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex)
-  end do
+
+  call CCTK_GroupIndex ( varindex, 'ehfinder::f' )
+  call CCTK_GroupIndex ( rhsindex, 'ehfinder::sf' )
+
+  ierr_cum = ierr_cum + MolRegisterEvolvedGroup ( varindex, rhsindex )
 
   if ( evolve_generators .gt. 0 ) then
 
@@ -47,37 +43,6 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
 
   end if
  
-!  call CCTK_VarIndex(varindex, 'admbase::gxx')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::gxy')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::gxz')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::gyy')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::gyz')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::gzz')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::alp')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::betax')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::betay')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-!  call CCTK_VarIndex(varindex, 'admbase::betaz')
-!  call MoL_RegisterPrimitive(ierr, varindex)
-!  ierr_cum = ierr_cum + ierr
-  print*,'MoLRegister ', ierr_cum
   if ( ierr_cum .gt. 0 ) then
     call CCTK_WARN(0,'Problems registering variables with MoL')
   end if