Forgot these files in the last commit.

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

1 files changed, 252 insertions, 239 deletions

diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90
index 48e946a..746fd3f 100644
--- a/src/EHFinder_Sources.F90
+++ b/src/EHFinder_Sources.F90
@@ -7,205 +7,217 @@
 
 subroutine EHFinder_Sources(CCTK_ARGUMENTS)
 
+  use EHFinder_mod
+
   implicit none
 
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
-  CCTK_INT :: i, j, k, nx, ny, nz, ngx, ngy, ngz, ierr
-  CCTK_INT :: ix1, ix2, iy1, iy2, iz1, iz2
-  CCTK_INT :: jx1, jx2, jy1, jy2, jz1, jz2
-!  CCTK_REAL, allocatable, dimension(:,:,:) :: dfx, dfy, dfz, dfsq
-  CCTK_REAL, allocatable, dimension(:,:,:) :: dfsq
+  CCTK_INT :: i, j, k
   CCTK_REAL :: idx, idy, idz
   CCTK_REAL :: a, b, c
   CCTK_REAL :: g4tt, g4tx, g4ty, g4tz, g4xx, g4xy, g4xz, g4yy, g4yz, g4zz
   CCTK_REAL :: idetg, alp2, tmp1, tmp2, tmp3
-  CCTK_REAL :: dfsq, ratio
-
-  nx = cctk_lsh(1)
-  ny = cctk_lsh(2)
-  nz = cctk_lsh(3)
-  
-  ngx = cctk_nghostzones(1)
-  ngy = cctk_nghostzones(2)
-  ngz = cctk_nghostzones(3)
-
-!  print*,nx,ny,nz
-!  print*,cctk_delta_space
-
-! allocate storage for derivatives of the level set function
-
-!  allocate ( dfx(nx,ny,nz), dfy(nx,ny,nz), dfz(nx,ny,nz), dfsq(nx,ny,nz), stat=ierr )
-  allocate ( dfsq(nx,ny,nz), stat=ierr )
-  if ( ierr .gt. 0 ) then
-    call CCTK_WARN(0,'Unable to allocate work arrays')
-  end if
+  CCTK_REAL :: ratio
+  CCTK_REAL, dimension(3) :: maxpos
 
 ! calculate 1/(2*delta) in each direction
 
-  idx = 1.0 / ( 2.0 * cctk_delta_space(1) )
-  idy = 1.0 / ( 2.0 * cctk_delta_space(2) )
-  idz = 1.0 / ( 2.0 * cctk_delta_space(3) )
-
-  if ( cctk_bbox(1) .eq. 1 ) then
-    ix1 = 0
-    jx1 = 1
-  else
-    ix1 = ngx
-    jx1 = ngx
-  end if
-  if ( cctk_bbox(2) .eq. 1 ) then
-    ix2 = 0
-    jx2 = 1
-  else
-    ix2 = ngx
-    jx2 = ngx
-  end if
-
-  if ( cctk_bbox(3) .eq. 1 ) then
-    iy1 = 0
-    jy1 = 1
-  else 
-    iy1 = ngy
-    jy1 = ngy
-  endif
-  if ( cctk_bbox(4) .eq. 1 ) then
-    iy2 = 0
-    jy2 = 1
-  else
-    iy2 = ngy
-    jy2 = ngy
-  endif
-
-  if ( cctk_bbox(5) .eq. 1 ) then
-    iz1 = 0
-    jz1 = 1
-  else 
-    iz1 = ngz
-    jz1 = ngz
-  endif
-  if ( cctk_bbox(6) .eq. 1 ) then
-    iz2 = 0
-    jz2 = 1
-  else
-    iz2 = ngz
-    jz2 = ngz
-  endif
-
-! calculate one sided derivatives for the outer boundary points
-! for x
-
-  if ( cctk_bbox(1) .eq. 1 ) then
-    dfx(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
-                                  -3.0*f(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
-                                   4.0*f(2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
-                                       f(3, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
-  endif
-  if ( cctk_bbox(2) .eq. 1 ) then
-    dfx(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
-                                   3.0*f(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
-                                   4.0*f(nx-1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
-                                       f(nx-2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
-  endif
-
-! for y
-
-  if ( cctk_bbox(3) .eq. 1 ) then
-    dfy(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) = idy * ( &
-                                  -3.0*f(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) + &
-                                   4.0*f(1+ix1:nx-ix2, 2, 1+iz1:nz-iz2) - &
-                                       f(1+ix1:nx-ix2, 3, 1+iz1:nz-iz2) )
-  endif
-  if ( cctk_bbox(4) .eq. 1 ) then
-    dfy(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) = idy * ( &
-                                   3.0*f(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) - &
-                                   4.0*f(1+ix1:nx-ix2, ny-1, 1+iz1:nz-iz2) + &
-                                       f(1+ix1:nx-ix2, ny-2, 1+iz1:nz-iz2) )
-  endif
-
-! for z
-
-  if ( cctk_bbox(5) .eq. 1 ) then
-    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) = idz * ( &
-                                  -3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) + &
-                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2) - &
-                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 3) )
-  endif
-  if ( cctk_bbox(6) .eq. 1 ) then
-    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) = idz * ( &
-                                   3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) - &
-                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-1) + &
-                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-2) )
-  endif
-
-! calculate normal centered derivatives for interior points
-
-  if ( one_sided ) then
-    do k = 1+iz1, nz-iz2
-      do j = 1+iy1, ny-iy2
-        do i = 1+jx1, nx-jx2
-          if ( f(i,j,k) - f(i-1,j,k) .ne. 0 ) then
-            ratio = abs ( f(i+1,j,k) - f(i,j,k) ) / ( f(i,j,k) - f(i-1,j,k) )
-            if ( ratio .lt. 1 ) then
-              dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
-            else
-              dfx(i,j,k) = idx * ( 3.0*f(i,j,k) - 4.0*f(i-1,j,k) + f(i-2,j,k) )
-            end if
+  idx = half / cctk_delta_space(1)
+  idy = half / cctk_delta_space(2)
+  idz = half / cctk_delta_space(3)
+
+  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
+
+! calculate one sided derivatives where necessary for x, y and z
+
+!  where ( eh_mask .ge. 0 )
+!    where ( btest(eh_mask,0) )
+!      dfx = idx * ( -three * f + four * cshift(f,1,1) - cshift(f,2,1) )
+!    end where
+!    where ( btest(eh_mask,1) )
+!      dfx = idx * ( three * f - four * cshift(f,-1,1) + cshift(f,-2,1) )
+!    end where
+!    where ( btest(eh_mask,2) )
+!      dfy = idy * ( -three * f + four * cshift(f,1,2) - cshift(f,2,2) )
+!    end where
+!    where ( btest(eh_mask,3) )
+!      dfy = idy * ( three * f - four * cshift(f,-1,2) + cshift(f,-2,2) )
+!    end where
+!    where ( btest(eh_mask,4) )
+!      dfz = idz * ( -three * f + four * cshift(f,1,3) - cshift(f,2,3) )
+!    end where
+!    where ( btest(eh_mask,5) )
+!      dfz = idz * ( three * f - four * cshift(f,-1,3) + cshift(f,-2,3) )
+!    end where
+!  
+!    where ( ( .not. btest(eh_mask,0) ) .and. ( .not. btest(eh_mask,1) ) )
+!      dfx =  idx * ( cshift(f,1,1) - cshift(f,-1,1) )
+!    end where
+!    where ( ( .not. btest(eh_mask,2) ) .and. ( .not. btest(eh_mask,3) ) )
+!      dfy =  idy * ( cshift(f,1,2) - cshift(f,-1,2) )
+!    end where
+!    where ( ( .not. btest(eh_mask,4) ) .and. ( .not. btest(eh_mask,5) ) )
+!      dfz =  idz * ( cshift(f,1,3) - cshift(f,-1,3) )
+!    end where
+!  end where
+
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        if ( eh_mask(i,j,k) .ge. 0 ) then
+          if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
+               ( .not. btest(eh_mask(i,j,k),1) ) ) then
+            dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
+          else if ( btest(eh_mask(i,j,k),0) ) then
+            dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
+                                  four * f(i+1,j,k) - f(i+2,j,k) )
           else
-            dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
+            dfx(i,j,k) = idx * ( three * f(i,j,k) - &
+                                 four * f(i-1,j,k) + f(i-2,j,k) )
           end if
-        end do
-      end do
-    end do
-          
-    do k = 1+iz1, nz-iz2
-      do j = 1+jy1, ny-jy2
-        do i = 1+ix1, nx-ix2
-          if ( f(i,j,k) - f(i,j-1,k) .ne. 0 ) then
-            ratio = abs ( f(i,j+1,k) - f(i,j,k) ) / ( f(i,j,k) - f(i,j-1,k) )
-            if ( ratio .lt. 1 ) then
-              dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
-            else
-              dfy(i,j,k) = idy * ( 3.0*f(i,j,k) - 4.0*f(i,j-1,k) + f(i,j-2,k) )
-            end if
+          if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
+               ( .not. btest(eh_mask(i,j,k),3) ) ) then
+            dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
+          else if ( btest(eh_mask(i,j,k),2) ) then
+             dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
+                                   four * f(i,j+1,k) - f(i,j+2,k) )
           else
-            dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
+             dfy(i,j,k) = idy * ( three * f(i,j,k) - &
+                                  four * f(i,j-1,k) + f(i,j-2,k) )
           end if
-        end do
-      end do
-    end do
-            
-    do k = 1+jz1, nz-jz2
-      do j = 1+iy1, ny-iy2
-        do i = 1+ix1, nx-ix2
-          if ( f(i,j,k) - f(i,j,k-1) .ne. 0 ) then
-            ratio = abs ( f(i,j,k+1) - f(i,j,k) ) / ( f(i,j,k) - f(i,j,k-1) )
-            if ( ratio .lt. 1 ) then
-              dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
-            else
-              dfz(i,j,k) = idz * ( 3.0*f(i,j,k) - 4.0*f(i,j,k-1) + f(i,j,k-2) )
-            end if
+          if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
+               ( .not. btest(eh_mask(i,j,k),5) ) ) then
+            dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
+          else if ( btest(eh_mask(i,j,k),4) ) then
+             dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
+                                   four * f(i,j,k+1) - f(i,j,k+2) )
           else
-            dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
+             dfz(i,j,k) = idz * ( three * f(i,j,k) - &
+                                  four * f(i,j,k-1) + f(i,j,k-2) )
           end if
-        end do
+        end if
       end do
     end do
-
-  else
-
-    dfx(1+jx1:nx-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = &
-                  idx * ( f(2+jx1:nx+1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
-                          f(jx1:nx-1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
-    dfy(1+ix1:nx-ix2, 1+jy1:ny-jy2, 1+iz1:nz-iz2) = &
-                  idy * ( f(1+ix1:nx-ix2, 2+jy1:ny+1-jy2, 1+iz1:nz-iz2) - &
-                          f(1+ix1:nx-ix2, jy1:ny-1-jy2, 1+iz1:nz-iz2) )
-    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1+jz1:nz-jz2) = &
-                  idz * ( f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2+jz1:nz+1-jz2) - &
-                          f(1+ix1:nx-ix2, 1+iy1:ny-iy2, jz1:nz-1-jz2) )
-  end if
+  end do
+
+    
+!  if ( cctk_bbox(1) .eq. 1 ) then
+!    dfx(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
+!                                  -3.0*f(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
+!                                   4.0*f(2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
+!                                       f(3, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
+!  endif
+!  if ( cctk_bbox(2) .eq. 1 ) then
+!    dfx(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
+!                                   3.0*f(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
+!                                   4.0*f(nx-1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
+!                                       f(nx-2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
+!  endif
+!
+!! for y
+!
+!  if ( cctk_bbox(3) .eq. 1 ) then
+!    dfy(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) = idy * ( &
+!                                  -3.0*f(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) + &
+!                                   4.0*f(1+ix1:nx-ix2, 2, 1+iz1:nz-iz2) - &
+!                                       f(1+ix1:nx-ix2, 3, 1+iz1:nz-iz2) )
+!  endif
+!  if ( cctk_bbox(4) .eq. 1 ) then
+!    dfy(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) = idy * ( &
+!                                   3.0*f(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) - &
+!                                   4.0*f(1+ix1:nx-ix2, ny-1, 1+iz1:nz-iz2) + &
+!                                       f(1+ix1:nx-ix2, ny-2, 1+iz1:nz-iz2) )
+!  endif
+!
+!! for z
+!
+!  if ( cctk_bbox(5) .eq. 1 ) then
+!    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) = idz * ( &
+!                                  -3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) + &
+!                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2) - &
+!                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 3) )
+!  endif
+!  if ( cctk_bbox(6) .eq. 1 ) then
+!    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) = idz * ( &
+!                                   3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) - &
+!                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-1) + &
+!                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-2) )
+!  endif
+!
+!! calculate normal centered derivatives for interior points
+!
+!  if ( one_sided ) then
+!    do k = 1+iz1, nz-iz2
+!      do j = 1+iy1, ny-iy2
+!        do i = 1+jx1, nx-jx2
+!          if ( f(i,j,k) - f(i-1,j,k) .ne. 0 ) then
+!            ratio = abs ( f(i+1,j,k) - f(i,j,k) ) / ( f(i,j,k) - f(i-1,j,k) )
+!            if ( ratio .lt. 1 ) then
+!              dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
+!            else
+!              dfx(i,j,k) = idx * ( 3.0*f(i,j,k) - 4.0*f(i-1,j,k) + f(i-2,j,k) )
+!            end if
+!          else
+!            dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
+!          end if
+!        end do
+!      end do
+!    end do
+!          
+!    do k = 1+iz1, nz-iz2
+!      do j = 1+jy1, ny-jy2
+!        do i = 1+ix1, nx-ix2
+!          if ( f(i,j,k) - f(i,j-1,k) .ne. 0 ) then
+!            ratio = abs ( f(i,j+1,k) - f(i,j,k) ) / ( f(i,j,k) - f(i,j-1,k) )
+!            if ( ratio .lt. 1 ) then
+!              dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
+!            else
+!              dfy(i,j,k) = idy * ( 3.0*f(i,j,k) - 4.0*f(i,j-1,k) + f(i,j-2,k) )
+!            end if
+!          else
+!            dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
+!          end if
+!        end do
+!      end do
+!    end do
+!            
+!    do k = 1+jz1, nz-jz2
+!      do j = 1+iy1, ny-iy2
+!        do i = 1+ix1, nx-ix2
+!          if ( f(i,j,k) - f(i,j,k-1) .ne. 0 ) then
+!            ratio = abs ( f(i,j,k+1) - f(i,j,k) ) / ( f(i,j,k) - f(i,j,k-1) )
+!            if ( ratio .lt. 1 ) then
+!              dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
+!            else
+!              dfz(i,j,k) = idz * ( 3.0*f(i,j,k) - 4.0*f(i,j,k-1) + f(i,j,k-2) )
+!            end if
+!          else
+!            dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
+!          end if
+!        end do
+!      end do
+!    end do
+!
+!  else
+!
+!    dfx(1+jx1:nx-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = &
+!                  idx * ( f(2+jx1:nx+1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
+!                          f(jx1:nx-1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
+!    dfy(1+ix1:nx-ix2, 1+jy1:ny-jy2, 1+iz1:nz-iz2) = &
+!                  idy * ( f(1+ix1:nx-ix2, 2+jy1:ny+1-jy2, 1+iz1:nz-iz2) - &
+!                          f(1+ix1:nx-ix2, jy1:ny-1-jy2, 1+iz1:nz-iz2) )
+!    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1+jz1:nz-jz2) = &
+!                  idz * ( f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2+jz1:nz+1-jz2) - &
+!                          f(1+ix1:nx-ix2, 1+iy1:ny-iy2, jz1:nz-1-jz2) )
+!  end if
 
 !  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
 !    print*,dfx(1,1+iy1:ny-iy2, 1+iz1:nz-iz2)
@@ -222,59 +234,66 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
 ! the 4 metric. This does not effect the solution to the quadratic
 ! equation, but avoids potential problems of dividing with the square 
 ! of the lapse in regions where it is close to zero.
+!
+!        
+!  if ( normalize ) then
+!    dfsq = 1.0 / sqrt(dfx**2+dfy**2+dfz**2)
+!    dfx = dfx*dfsq
+!    dfy = dfy*dfsq
+!    dfz = dfz*dfsq
+!  end if
 
-        
-  if ( normalize ) then
-    dfsq = 1.0 / sqrt(dfx**2+dfy**2+dfz**2)
-    dfx = dfx*dfsq
-    dfy = dfy*dfsq
-    dfz = dfz*dfsq
-  end if
-
-  do k = 1+iz1, nz-iz2
-    do j = 1+iy1, ny-iy2
-      do i = 1+ix1, nx-ix2
-        alp2 = alp(i,j,k)**2
-        g4tt = -1.0
-        g4tx = betax(i,j,k)
-        g4ty = betay(i,j,k)
-        g4tz = betaz(i,j,k)
-
-        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 = 1.0 / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 )
-
-        g4xx = tmp1 * idetg
-        g4xy = tmp2 * idetg
-        g4xz = tmp3 * idetg
-        g4yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
-        g4yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg
-        g4zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
-        
-        g4xx = alp2*g4xx - betax(i,j,k)**2
-        g4xy = alp2*g4xy - betax(i,j,k)*betay(i,j,k)
-        g4xz = alp2*g4xz - betax(i,j,k)*betaz(i,j,k)
-        g4yy = alp2*g4yy - betay(i,j,k)**2
-        g4yz = alp2*g4yz - betay(i,j,k)*betaz(i,j,k)
-        g4zz = alp2*g4zz - betaz(i,j,k)**2
-      
-        a = g4tt
-        b = 2.0 * ( g4tx * dfx(i,j,k) + &
-                    g4ty * dfy(i,j,k) + &
-                    g4tz * dfz(i,j,k) )
-        c = g4xx * dfx(i,j,k)**2 + &
-            g4yy * dfy(i,j,k)**2 + &
-            g4zz * dfz(i,j,k)**2 + &
-            2.0 * ( g4xy * dfx(i,j,k) * dfy(i,j,k) + &
-                    g4xz * dfx(i,j,k) * dfz(i,j,k) + &
-                    g4yz * dfy(i,j,k) * dfz(i,j,k) )
-
-        if ( b .lt. 0. ) then
-          sf(i,j,k) = - ( -b + sqrt(b**2 - 4*a*c) ) / ( 2.0*a )
+!  do k = 1+iz1, nz-iz2
+!    do j = 1+iy1, ny-iy2
+!      do i = 1+ix1, nx-ix2
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        if ( eh_mask(i,j,k) .ge. 0 ) then
+          alp2 = alp(i,j,k)**2
+          g4tt = -one
+          g4tx = betax(i,j,k)
+          g4ty = betay(i,j,k)
+          g4tz = betaz(i,j,k)
+
+          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 )
+
+          g4xx = tmp1 * idetg
+          g4xy = tmp2 * idetg
+          g4xz = tmp3 * idetg
+          g4yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
+          g4yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg
+          g4zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
+
+          g4xx = alp2*g4xx - betax(i,j,k)**2
+          g4xy = alp2*g4xy - betax(i,j,k)*betay(i,j,k)
+          g4xz = alp2*g4xz - betax(i,j,k)*betaz(i,j,k)
+          g4yy = alp2*g4yy - betay(i,j,k)**2
+          g4yz = alp2*g4yz - betay(i,j,k)*betaz(i,j,k)
+          g4zz = alp2*g4zz - betaz(i,j,k)**2
+
+          a = g4tt
+          b = two * ( g4tx * dfx(i,j,k) + &
+                      g4ty * dfy(i,j,k) + &
+                      g4tz * dfz(i,j,k) )
+          c = g4xx * dfx(i,j,k)**2 + &
+              g4yy * dfy(i,j,k)**2 + &
+              g4zz * dfz(i,j,k)**2 + &
+              two * ( g4xy * dfx(i,j,k) * dfy(i,j,k) + &
+                      g4xz * dfx(i,j,k) * dfz(i,j,k) + &
+                      g4yz * dfy(i,j,k) * dfz(i,j,k) )
+
+          if ( b .lt. zero ) then
+            sf(i,j,k) = - ( -b + sqrt(b**2 - four*a*c) ) / ( two*a )
+          else
+            sf(i,j,k) = - two*c / ( -b - sqrt(b**2 - 4*a*c) )
+          end if
         else
-          sf(i,j,k) = - 2.0*c / ( -b - sqrt(b**2 - 4*a*c) )
+          sf(i,j,k) = zero
         end if
 !        if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
 !          if ( abs(b) .gt. 0.00000001 ) then
@@ -292,20 +311,14 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS)
       end do
     end do
   end do 
-!  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
-!    print*,sf
-!    print*
-!  end if
-!  print*,cctk_bbox
-!  print*,cctk_nghostzones
+ 
+!  print*,maxval(sf)
+!  print*
+!  sf = zero
 
-!  deallocate ( dfx, dfy, dfz, dfsq )
-  deallocate ( dfsq )
 !  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
 !    call CCTK_WARN(0,'stopping')
 !  end if
   
   return
 end subroutine EHFinder_Sources
-
-