Delete this file since the functionality is now done in

EHFinder_ReInitialize.F90


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

1 files changed, 0 insertions, 645 deletions

diff --git a/src/EHFinder_ReParametrize.F90 b/src/EHFinder_ReParametrize.F90
deleted file mode 100644
index 7559395..0000000
--- a/src/EHFinder_ReParametrize.F90
+++ /dev/null
@@ -1,645 +0,0 @@
-! Re-Parametrization of the level set function
-! $Header$
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-! Control routine for reparametrization of the level set function.
-
-subroutine EHFinder_ReParametrizeControl(CCTK_ARGUMENTS)
-
-  use EHFinder_mod
-
-  implicit none
-
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_FUNCTIONS
- 
-  CCTK_REAL, dimension(eh_number_level_sets) :: fmin, fmax, fminloc, fmaxloc
-  CCTK_INT :: i, j, k, l
-
-!  print*,'EHFinder_ReParametrize1 Entered'
-! Initialize the control variables. The default means no re-parametrization.
-  re_param_control_pde = 0
-  re_param_control_approx = 0
-
-! If the method is 'approx' or 'mixed' check if it is time for an
-! approximation re-parametrization and set the control variable to 1.
-  if ( ( CCTK_EQUALS ( re_param_method, 'approx' ) ) .or. &
-       ( CCTK_EQUALS ( re_param_method, 'mixed' ) ) ) then
-    if ( reparametrize_every_approx .gt. 0 ) then
-      if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then
-        re_param_control_approx = 1
-      end if
-    end if
-  end if
-
-! If the method is 'pde' or 'mixed' check if it is time for a pde
-! re-parametrization and set the control variable to 1. Note it overrides
-! the approximation re-parametrization.
-  if ( ( CCTK_EQUALS ( re_param_method, 'pde' ) ) .or. &
-       ( CCTK_EQUALS ( re_param_method, 'mixed' ) ) ) then
-    if ( reparametrize_every_pde .gt. 0 ) then
-      if ( mod(cctk_iteration,reparametrize_every_pde) .eq. 0 ) then
-        re_param_control_approx = 0
-        re_param_control_pde = 1
-      end if 
-    end if 
-  end if 
-    
-! If it is time for a re-parametrization get reduction handles for
-! minimum and maximum reductions.
-  if ( ( re_param_control_approx .gt. 0 ) .or. &
-       ( re_param_control_pde .gt. 0 ) ) then
-    call CCTK_ReductionArrayHandle ( max_handle, 'maximum' )
-    if ( max_handle .lt. 0 ) then
-      call CCTK_WARN(0,'Could not obtain a handle for maximum reduction')
-    end if
-    call CCTK_ReductionArrayHandle ( min_handle, 'minimum' )
-    if ( min_handle .lt. 0 ) then
-      call CCTK_WARN(0,'Could not obtain a handle for minimum reduction')
-    end if
-
-!      if ( cctk_bbox(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
-!      if ( cctk_bbox(2) .eq. 0 ) eh_mask(nx-ngx+1:nx,:,:) = -2
-!      if ( cctk_bbox(3) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
-!      if ( cctk_bbox(4) .eq. 0 ) eh_mask(:,ny-ngy+1:ny,:) = -2
-!      if ( cctk_bbox(5) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
-!      if ( cctk_bbox(6) .eq. 0 ) eh_mask(:,:,nz-ngz+1:nz) = -2
-  
-! Set the courant factor for the 'Euler' re-parametrization scheme.
-    if ( CCTK_EQUALS ( re_param_int_method, 'Euler' ) ) then
-!      hfac = one / eight
-      hfac = one / four
-    end if
-
-! Set the courant factor for the 'rk2' re-parametrization scheme.
-    if ( CCTK_EQUALS ( re_param_int_method, 'rk2' ) ) then
-!      hfac = half
-      hfac = one / four
-    end if
-
-! Set up counter for the number of iterations. Copy f into ftmp and
-! fbak and eh_mask into eh_mask_bak (in case the re-parametrization has 
-! to be undone) and find the min and max values of f.
-    ncalls = 0
-    ftmp = f
-    fbak = f
-    eh_mask_bak = eh_mask
-!    rep_mask = 0
-    do l = 1, eh_number_level_sets
-      fminloc(l) = minval(f(:,:,:,l))
-      fmaxloc(l) = maxval(f(:,:,:,l))
-    end do
-
-    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, &
-                                  fmaxloc, fmax, eh_number_level_sets, &
-                                  CCTK_VARIABLE_REAL )
-    if ( ierr .ne. 0 ) then
-      call CCTK_WARN(0,'Reduction of fmax failed')
-    end if
-    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
-                                      fminloc, fmin, eh_number_level_sets, &
-                                      CCTK_VARIABLE_REAL )
-    if ( ierr .ne. 0 ) then
-      call CCTK_WARN(0,'Reduction of fmin failed')
-    end if
-
-    reparam_this_level_set = .true.
-
-! If fmin and fmax have the same sign, there is no surface present and
-! re-parametrization will not be performed for the given surface.
-    do l = 1, eh_number_level_sets
-      if ( fmin(l) * fmax(l) .gt. zero ) then
-        reparam_this_level_set(l) = .false.
-      end if
-    end do
-
-! If none of the surfaces should be reparametrized set the control variables
-! to zero.
-    if ( all ( .not. reparam_this_level_set ) ) then
-      re_param_control_approx = 0
-      re_param_control_pde = 0
-      call CCTK_INFO ( 'No zero-points of the level set functions. No re-parametrization performed.' )
-    end if
-  end if
-!  print*,'EHFinder_ReParametrizeControl Exited'
-
-end subroutine EHFinder_ReParametrizeControl
-
-
-!subroutine EHFinder_ReParametrizeRK2_1(CCTK_ARGUMENTS)
-!
-!  use EHFinder_mod
-!
-!  implicit none
-!
-!  DECLARE_CCTK_PARAMETERS
-!  DECLARE_CCTK_ARGUMENTS
-!  DECLARE_CCTK_FUNCTIONS
-!
-!  CCTK_INT :: i, j, k
-!  CCTK_REAL :: idx, idy, idz
-!  CCTK_REAL :: al, ar, bl, br, cl, cr
-!  CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
-!  CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
-!
-!!  print*,'EHFinder_ReParametrizeRK2_1 Entered'
-!  h = hfac*minval(cctk_delta_space)
-!
-!  if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then
-!
-!# include "include/centered_second.h"
-!
-!  end if
-!
-!  if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then
-!
-!# include "include/upwind_first.h"
-!
-!  end if
-!
-!  if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then
-!
-!# include "include/upwind_second.h"
-!
-!  end if
-!
-!  where ( eh_mask .ge. 0 )
-!    sftmp =  sqrt( dfx**2 + dfy**2 + dfz**2 )
-!    sftmp = - half * h * f / sqrt(f**2+one) * ( sftmp - one )
-!    f = f + sftmp
-!  elsewhere
-!    sftmp = one
-!  end where
-!!  print*,'EHFinder_ReParametrizeRK2_1 Exited'
-!
-!end subroutine EHFinder_ReParametrizeRK2_1
-
-
-!subroutine EHFinder_ReParametrizeRK2_2(CCTK_ARGUMENTS)
-!
-!
-!  use EHFinder_mod
-!
-!  implicit none
-!
-!  DECLARE_CCTK_PARAMETERS
-!  DECLARE_CCTK_ARGUMENTS
-!  DECLARE_CCTK_FUNCTIONS
-!
-!  CCTK_INT :: i, j, k
-!  CCTK_REAL :: idx, idy, idz, maxf, maxdfloc, maxdf, mindfloc, mindf
-!  CCTK_REAL :: al, ar, bl, br, cl, cr
-!  CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
-!  CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
-!  character(len=128) :: info_message
-!
-!!  print*,'EHFinder_ReParametrizeRK2_2 Entered'
-!
-!  if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then
-!
-!# include "include/centered_second.h"
-!
-!  end if
-!
-!  if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then
-!
-!# include "include/upwind_first.h"
-!
-!  end if
-!
-!  if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then
-!
-!# include "include/upwind_second.h"
-!
-!  end if
-!
-!  where ( eh_mask .ge. 0 )
-!
-!    sftmp =  sqrt( dfx**2 + dfy**2 + dfz**2 )
-!    sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one )
-!    f = ftmp + sftmp
-!  elsewhere
-!    sftmp = one
-!  end where
-!
-!  maxdfloc = zero
-!  mindfloc = 1.0d23
-!  do k = kzl, kzr
-!    do j = jyl, jyr
-!      do i =ixl, ixr
-!        if ( eh_mask(i,j,k) .ge. 0 ) then
-!          maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k)) )
-!          mindfloc = min ( mindfloc, abs(sftmp(i,j,k)) )
-!        end if
-!      end do
-!    end do
-!  end do
-!
-!  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, &
-!                                    maxdfloc, maxdf, CCTK_VARIABLE_REAL )
-!  if ( ierr .ne. 0 ) then
-!    call CCTK_WARN(0,'Reduction of maxdf failed')
-!  end if
-!  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, &
-!                                    mindfloc, mindf, CCTK_VARIABLE_REAL )
-!  if ( ierr .ne. 0 ) then
-!    call CCTK_WARN(0,'Reduction of mindf failed')
-!  end if
-!
-!  if ( maxdf .lt. h*minval(cctk_delta_space)**2 ) re_param_control_pde = 0
-!
-!  ncalls = ncalls + 1
-!
-!  if ( re_param_control_pde .eq. 0 ) then
-!    write(info_message,'(a35,i5,a12)') &
-!        'PDE re-parametrization complete in ',ncalls,' iterations.'
-!    call CCTK_INFO(info_message)
-!  end if
-!
-!  if ( ncalls .gt. re_param_max_iter ) then
-!    call CCTK_WARN(1,'Re-parametrization failed to converge')
-!    re_param_control_pde = 0
-!  end if
-!
-!  ftmp = f
-!!  print*,ncalls,maxdf,mindf
-!!  print*,mx1,my1,mz2,mx2,my2,mz2
-!!  print*,f(mx1,my1,mz1),f(mx2,my2,mz2),eh_mask(mx1,my1,mz1),eh_mask(mx2,my2,mz2)
-!!   print*,f(27,53,27),dfsq(27,53,27),ftmp(27,53,27)
-!!  
-!!  call CCTK_INFO('I was called')
-!!  print*,'EHFinder_ReParametrizeRK2_2 Exited'
-!
-!end subroutine EHFinder_ReParametrizeRK2_2
-
-
-subroutine EHFinder_ReParametrizeEuler(CCTK_ARGUMENTS)
-
-  use EHFinder_mod
-
-  implicit none
-
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_FUNCTIONS
-
-  CCTK_INT :: i, j, k, l
-  CCTK_REAL :: idx, idy, idz, maxf
-  CCTK_REAL, dimension(eh_number_level_sets) :: maxdfloc, maxdf, &
-                                                mindfloc, mindf
-  CCTK_REAL :: al, ar, bl, br, cl, cr
-  CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
-  CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
-  character(len=128) :: info_message
-
-!  print*,'EHFinder_ReParametrizeEuler Entered'
-  h = hfac*minval(cctk_delta_space)
-
-  if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then
-
-# include "include/centered_second.h"
-
-  end if
-
-  if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then
-
-# include "include/upwind_first.h"
-
-  end if
-
-  if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then
-
-# include "include/upwind_second.h"
-
-  end if
-
-  where ( eh_mask .ge. 0 )
-    sftmp =  sqrt( dfx**2 + dfy**2 + dfz**2 )
-    sftmp =  - h * f / sqrt(f**2+one) * ( sftmp - one )
-    f = f + sftmp
-  elsewhere
-    sftmp = one
-  end where
-!  print*,f(2:4,2:4,2:4,5)
-!  print*,sftmp(3,3,3,5)
-!  stop
-  
-!  do l = 1, eh_number_level_sets
-!    maxdfloc = zero
-!    mindfloc = 1.0d23
-!    do k = kzl, kzr
-!      do j = jyl, jyr
-!        do i = ixl, ixr
-!          if ( eh_mask(i,j,k,l) .ge. 0 ) then
-!            maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k,l)) ) 
-!            mindfloc = min ( mindfloc, abs(sftmp(i,j,k,l)) ) 
-!          end if
-!        end do
-!      end do
-!    end do
-!  end do
-  
-  do l = 1, eh_number_level_sets
-    maxdfloc(l) = maxval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), &
-                             mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0)
-    mindfloc(l) = minval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), &
-                             mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0)
-  end do
-
-  call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, &
-                                    maxdfloc, maxdf, eh_number_level_sets, &
-                                    CCTK_VARIABLE_REAL )
-  if ( ierr .ne. 0 ) then
-    call CCTK_WARN(0,'Reduction of maxdf failed')
-  end if
-  call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
-                                    mindfloc, mindf, eh_number_level_sets, &
-                                    CCTK_VARIABLE_REAL )
-  if ( ierr .ne. 0 ) then
-    call CCTK_WARN(0,'Reduction of mindf failed')
-  end if
-
-  if ( all ( maxdf .lt. h*minval(cctk_delta_space)**2 ) ) then
-    re_param_control_pde = 0
-  endif
-
-  ncalls = ncalls + 1
-
-  if ( re_param_control_pde .eq. 0 ) then
-    write(info_message,'(a35,i5,a12)') &
-        'PDE re-parametrization complete in ',ncalls,' iterations.'
-    call CCTK_INFO(info_message)
-  end if
-
-  if ( ncalls .gt. re_param_max_iter ) then
-    call CCTK_WARN(1,'Re-parametrization failed to converge')
-    re_param_control_pde = 0
-  end if
-
-!  print*,ncalls
-!  print*,maxdf
-!  print*,mindf
-!  print*,'EHFinder_ReParametrizeEuler Exited'
-
-end subroutine EHFinder_ReParametrizeEuler
-
-!subroutine EHFinder_ReParametrize5(CCTK_ARGUMENTS)
-!
-!  use EHFinder_mod
-!
-!  implicit none
-!
-!  DECLARE_CCTK_PARAMETERS
-!  DECLARE_CCTK_ARGUMENTS
-!  DECLARE_CCTK_FUNCTIONS
-!
-!  CCTK_INT :: i, j, k
-!  CCTK_INT :: rep_countloc, rep_totalloc
-!  CCTK_REAL :: idx2, idy2, idz2, idx, idy, idz
-!  CCTK_REAL :: d2fdx2, d2fdy2, d2fdz2, d2fdxdy, d2fdxdz, d2fdydz
-!  CCTK_REAL :: n_x, n_y, n_z
-!  CCTK_REAL :: a, b, c, lambda
-!
-!  interface
-!    CCTK_REAL function solve_quadratic ( a, b, c )
-!      CCTK_REAL, intent(in) :: a, b, c
-!    end function solve_quadratic
-!  end interface
-!
-!  if ( re_param_control_approx .eq. 1 ) then
-!    if ( reparametrize_every_approx .gt. 0 ) then
-!      if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then
-!        call CCTK_ReductionArrayHandle ( sum_handle, 'sum' )
-!        if ( sum_handle .lt. 0 ) then
-!          call CCTK_WARN(0,'Could not obtain a handle for sum reduction')
-!        end if
-!
-!# include "include/centered_second.h"
-!
-!        ncalls = 0
-!
-!        rep_totalloc = count ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 )
-!
-!        rep_countloc = 0
-!
-!        ftmp = f
-!
-!        ixr = min ( ixr, nx-1 )
-!        jyr = min ( jyr, ny-1 )
-!        kzr = min ( kzr, nz-1 )
-!
-!        do k = kzl, kzr
-!          do j = jyl, jyr
-!            do i = ixl, ixr
-!              if ( eh_mask(i,j,k) .eq. 0 ) then
-!                if ( ( f(i,j,k)*f(i-1,j,k) .le. zero ) .or. &
-!                     ( f(i,j,k)*f(i+1,j,k) .le. zero ) .or. &
-!                     ( f(i,j,k)*f(i,j-1,k) .le. zero ) .or. &
-!                     ( f(i,j,k)*f(i,j+1,k) .le. zero ) .or. &
-!                     ( f(i,j,k)*f(i,j,k-1) .le. zero ) .or. &
-!                     ( f(i,j,k)*f(i,j,k+1) .le. zero ) ) then
-!
-!                  rep_mask(i,j,k) = 1
-!                  rep_countloc = rep_countloc + 1
-!
-!                  d2fdx2 = idx2 * ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) )
-!                  d2fdy2 = idy2 * ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) )
-!                  d2fdz2 = idz2 * ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) )
-!
-!                  d2fdxdy = idx * idy * ( f(i-1,j-1,k) - f(i+1,j-1,k) - &
-!                                          f(i-1,j+1,k) + f(i+1,j+1,k) )
-!                  d2fdxdz = idx * idz * ( f(i-1,j,k-1) - f(i+1,j,k-1) - &
-!                                          f(i-1,j,k+1) + f(i+1,j,k+1) )
-!                  d2fdydz = idy * idz * ( f(i,j-1,k-1) - f(i,j+1,k-1) - &
-!                                          f(i,j-1,k+1) + f(i,j+1,k+1) )
-!
-!                  c = f(i,j,k)
-!                  b = sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
-!                  n_x = dfx(i,j,k) / b
-!                  n_y = dfy(i,j,k) / b
-!                  n_z = dfz(i,j,k) / b
-!                  a = half*( n_x**2*d2fdx2 + n_y**2*d2fdy2 + n_z**2*d2fdz2 ) +&
-!                      n_x*n_y*d2fdxdy + n_x*n_z*d2fdxdz + n_y*n_z*d2fdydz
-!
-!                  if ( a .ne. 0 ) then
-!                    lambda = solve_quadratic ( a, b, c )
-!                  else
-!                    lambda = abs ( c / b )
-!                  endif
-!                  if ( lambda .eq. huge ) then
-!                    lambda = abs ( c / b )
-!                  end if
-!
-!                  if ( lambda .gt. two*delta ) then
-!                    print*,i,j,k
-!                    print*,f(i-1:i+1,j-1:j+1,k-1:k+1)
-!                    print*,eh_mask(i-1:i+1,j-1:j+1,k-1:k+1)
-!                    print*,a,b,c,lambda
-!                    call CCTK_WARN(0,'Re-parametrization failed')
-!                  end if
-!
-!                  ftmp(i,j,k) = sign ( lambda, f(i,j,k) )
-!!                  print*,i,j,k
-!!                  print*,lambda,f(i,j,k), ftmp(i,j,k)
-!!                  print*,f(i-1,j,k),f(i+1,j,k)
-!!                  print*,f(i,j-1,k),f(i,j+1,k)
-!!                  print*,f(i,j,k-1),f(i,j,k+1)
-!                end if
-!              end if
-!            end do
-!          end do
-!        end do
-!
-!        where ( eh_mask .ge. 0 ) f = ftmp
-!
-!!        print*,rep_countloc,rep_totalloc
-!        call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
-!                                    rep_countloc, rep_count, CCTK_VARIABLE_INT )
-!        if ( ierr .ne. 0 ) then
-!          call CCTK_WARN(0,'Sum of rep_count failed')
-!        end if
-!      
-!        call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
-!                                    rep_totalloc, rep_total, CCTK_VARIABLE_INT )
-!        if ( ierr .ne. 0 ) then
-!          call CCTK_WARN(0,'Sum of rep_total failed')
-!        end if
-!!        print*,'First step'
-!!        print*,rep_count,rep_total
-!!        call CCTK_WARN(0,'Finished')
-! 
-!!        call CartSymGN(ierr,cctkGH,'ehfinder::level_set')
-!!        call CartSymGN(ierr,cctkGH,'ehfinder::rep_mask')
-! 
-!      end if
-!    end if
-!  
-!  end if            
-!
-!end subroutine EHFinder_ReParametrize5
-
-
-!subroutine EHFinder_ReParametrize6(CCTK_ARGUMENTS)
-!
-!  use EHFinder_mod
-!
-!  implicit none
-!
-!  DECLARE_CCTK_PARAMETERS
-!  DECLARE_CCTK_ARGUMENTS
-!  DECLARE_CCTK_FUNCTIONS
-!
-!  CCTK_INT :: i, j, k, l, nneigh, rep_loc1, rep_loc2
-!  CCTK_INT, dimension(1) :: lmax
-!  CCTK_REAL, dimension(3) :: n0, n1, vec
-!  CCTK_REAL, dimension(0:5) :: dp1, dp2, fnew
-!  CCTK_REAL :: dfs1, dfs2
-!
-!# include "include/physical_part.h"
-!
-!  ncalls = ncalls + 1
-!
-!  rep_loc1 = 0
-!
-!  do k = kzl, kzr
-!    do j = jyl, jyr
-!      do i = ixl, ixr
-!        if (  ( eh_mask(i,j,k) .ge. 0 ) .and. ( rep_mask(i,j,k) .eq. 0 ) ) then
-!          nneigh = 0
-!          fnew = zero
-!          dfs1 = one / sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
-!          n0(1) = dfx(i,j,k) * dfs1
-!          n0(2) = dfy(i,j,k) * dfs1
-!          n0(3) = dfz(i,j,k) * dfs1
-!          do l = 0, 5
-!            if ( .not. btest(eh_mask(i,j,k),l) ) then
-!              if ( rep_mask(i+ix(l),j+jy(l),k+kz(l)) .eq. 1 ) then
-!                vec(1) = -ix(l) * cctk_delta_space(1)
-!                vec(2) = -jy(l) * cctk_delta_space(2)
-!                vec(3) = -kz(l) * cctk_delta_space(3)
-!                dfs2 = one / sqrt ( dfx(i+ix(l),j+jy(l),k+kz(l))**2 + &
-!                                    dfy(i+ix(l),j+jy(l),k+kz(l))**2 + &
-!                                    dfz(i+ix(l),j+jy(l),k+kz(l))**2 )
-!                n1(1) = dfx(i+ix(l),j+jy(l),k+kz(l)) * dfs2
-!                n1(2) = dfy(i+ix(l),j+jy(l),k+kz(l)) * dfs2
-!                n1(3) = dfz(i+ix(l),j+jy(l),k+kz(l)) * dfs2
-!                dp1(nneigh) = vec(1) * n0(1) + vec(2) * n0(2) + vec(3) * n0(3)
-!                dp2(nneigh) = vec(1) * n1(1) + vec(2) * n1(2) + vec(3) * n1(3)
-!                fnew(nneigh) = f(i+ix(l),j+jy(l),k+kz(l)) + &
-!                               half * ( dp1(nneigh) + dp2(nneigh) )
-!                nneigh = nneigh + 1
-!           
-!              end if
-!            end if
-!          end do
-!
-!          if ( nneigh .gt. 0 ) then
-!            f(i,j,k) = sum(fnew(0:nneigh-1)) / nneigh
-!            rep_mask(i,j,k) = 2
-!            rep_loc1 = rep_loc1 + 1 
-!          end if
-!        end if
-!      end do
-!    end do
-!  end do
-!
-!  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
-!                              rep_loc1, rep_loc2, CCTK_VARIABLE_INT )
-!  if ( ierr .ne. 0 ) then
-!    call CCTK_WARN(0,'Sum of rep_count failed')
-!  end if
-!!  print*,rep_loc1,rep_loc2
-!      
-!  rep_count = rep_count + rep_loc2
-!!  print*,ncalls,rep_count,rep_total
-!
-!  if ( rep_count .ge. rep_total ) then
-!    re_param_control_approx = 0
-!    call CCTK_INFO('Approximation re-parametrization complete')
-!  end if
-!
-!!  call CCTK_WARN(0,'Testing')
-!!  print*,'after'
-!!  print*,rep_mask(24,21,1:4),eh_mask(24,21,1:4)
-!!  print*,f(24,21,1:4)
-!!  print*
-!end subroutine EHFinder_ReParametrize6
-
-
-!subroutine EHFinder_ReParametrize7(CCTK_ARGUMENTS)
-!
-!  use EHFinder_mod
-!
-!  implicit none
-!
-!  DECLARE_CCTK_PARAMETERS
-!  DECLARE_CCTK_ARGUMENTS
-!  DECLARE_CCTK_FUNCTIONS
-!
-!  where ( rep_mask .eq. 2 ) rep_mask = 1
-!
-!end subroutine EHFinder_ReParametrize7
-
-
-!CCTK_REAL function solve_quadratic ( a, b, c )
-!
-!  use EHFinder_Constants
-!
-!  CCTK_REAL, intent(in) :: a, b, c
-!  CCTK_REAL :: d2, q
-!
-!  d2 = b**2 - four * a * c
-!
-!  if ( d2 .ge. zero ) then
-!    q = - half * ( b + sign(one,b) * sqrt(d2) )
-!
-!    solve_quadratic = min ( abs(q/a), abs(c/q) )
-!  else
-!    solve_quadratic = huge
-!  end if
-!
-!end function solve_quadratic