Complete rewrite to properly handle removal and addition of points to the

excised region in various symmetry modes and on multi processors. All
symmetry boundary condisions and syncs are scheduled in the schedule.ccl
file.


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

1 files changed, 315 insertions, 135 deletions

diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90
index 8b979cb..995b3a5 100644
--- a/src/EHFinder_SetMask.F90
+++ b/src/EHFinder_SetMask.F90
@@ -19,49 +19,51 @@ subroutine EHFinder_MaskInit(CCTK_ARGUMENTS)
 
   eh_mask = 0
 
-  if ( cctk_bbox(1) .eq. 1 ) then
+  # include "include/physical_part.h"
+
+  if ( ( cctk_bbox(1) .eq. 1 ) .and. ( ixl .eq. 1 ) ) then
     eh_mask(1,:,:) = eh_mask(1,:,:) + ll(0)
-  else
-    eh_mask(1:ngx,:,:) = -2
+!  else
+!    eh_mask(1:ngx,:,:) = -2
   end if
-  if ( cctk_bbox(2) .eq. 1 ) then
+  if ( ( cctk_bbox(2) .eq. 1 ) .and. ( ixr .eq. nx ) ) then
     eh_mask(nx,:,:) = eh_mask(nx,:,:) + ll(1)
-  else
-   eh_mask(nx-ngx+1:nx,:,:) = -2
+!  else
+!   eh_mask(nx-ngx+1:nx,:,:) = -2
   end if
-  if ( cctk_bbox(3) .eq. 1 ) then
+  if ( ( cctk_bbox(3) .eq. 1 ) .and. ( jyl .eq. 1 ) ) then
     eh_mask(:,1,:) = eh_mask(:,1,:) + ll(2)
-  else
-    eh_mask(:,1:ngy,:) = -2
+!  else
+!    eh_mask(:,1:ngy,:) = -2
   end if
-  if ( cctk_bbox(4) .eq. 1 ) then
+  if ( ( cctk_bbox(4) .eq. 1 ) .and. ( jyr .eq. ny ) ) then
     eh_mask(:,ny,:) = eh_mask(:,ny,:) + ll(3)
-  else
-    eh_mask(:,ny-ngy+1:ny,:) = -2
+!  else
+!    eh_mask(:,ny-ngy+1:ny,:) = -2
   end if
-  if ( cctk_bbox(5) .eq. 1 ) then
+  if ( ( cctk_bbox(5) .eq. 1 ) .and. ( kzl .eq. 1 ) ) then
     eh_mask(:,:,1) = eh_mask(:,:,1) + ll(4)
-  else
-    eh_mask(:,:,1:ngz) = -2
+!  else
+!    eh_mask(:,:,1:ngz) = -2
   end if
-  if ( cctk_bbox(6) .eq. 1 ) then
+  if ( ( cctk_bbox(6) .eq. 1 ) .and. ( kzr .eq. nz ) ) then
     eh_mask(:,:,nz) = eh_mask(:,:,nz) + ll(5)
-  else
-    eh_mask(:,:,nz-ngz+1:nz) = -2
+!  else
+!    eh_mask(:,:,nz-ngz+1:nz) = -2
   end if
 
-  if ( CCTK_Equals(domain,'octant') ) then
-    if ( cctk_lbnd(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
-    if ( cctk_lbnd(2) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
-    if ( cctk_lbnd(3) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
-  end if
+!  if ( CCTK_Equals(domain,'octant') ) then
+!    if ( cctk_lbnd(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
+!    if ( cctk_lbnd(2) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
+!    if ( cctk_lbnd(3) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
+!  end if
     
 !   call CCTK_WARN(0,'bbox')
   return
 end subroutine EHFinder_MaskInit
 
 
-subroutine EHFinder_SetMask(CCTK_ARGUMENTS)
+subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
 
   use EHFinder_mod
 
@@ -71,6 +73,7 @@ subroutine EHFinder_SetMask(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
 
+  CCTK_INT :: i, j, k
   logical :: active
 
   active = .false.
@@ -81,127 +84,304 @@ subroutine EHFinder_SetMask(CCTK_ARGUMENTS)
   end if
 
   if ( active ) then
-    mask_first = .false.
 
-    where ( ( eh_mask .ge. 0 ) .and. ( f .lt. -shell_width*delta ) ) 
-!    where ( ( f .lt. -shell_width*delta ) ) 
-     eh_mask = -1
-     f = 0
-    end where
+    # include "include/physical_part.h"
+
+    tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr)
+    ftmp(ixl:ixr,jyl:jyr,kzl:kzr) = f(ixl:ixr,jyl:jyr,kzl:kzr)
+
+!    print*,'Mask 1'
+!    print*,ixl,ixr,jyl,jyr,kzl,kzr
+!    print*,eh_mask(:,jyl,kzl)
+!    print*,f(:,jyl,kzl)
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
+
+          if ( eh_mask(i,j,k) .eq. -1 ) then
+
+            if ( eh_mask(i-1,j,k) .ge. 0 ) then
+              if ( f(i-1,j,k) - delta .gt. -shell_width * delta ) then
+                tm_mask(i,j,k) = 0
+                ftmp(i,j,k) = f(i-1,j,k) - delta
+              end if
+            end if
+
+            if ( eh_mask(i+1,j,k) .ge. 0 ) then
+              if ( f(i+1,j,k) - delta .gt. -shell_width * delta ) then
+                tm_mask(i,j,k) = 0
+                ftmp(i,j,k) = f(i+1,j,k) - delta
+              end if
+            end if
+
+            if ( eh_mask(i,j-1,k) .ge. 0 ) then
+              if ( f(i,j-1,k) - delta .gt. -shell_width * delta ) then
+                tm_mask(i,j,k) = 0
+                ftmp(i,j,k) = f(i,j-1,k) - delta
+              end if
+            end if
+
+            if ( eh_mask(i,j+1,k) .ge. 0 ) then
+              if ( f(i,j+1,k) - delta .gt. -shell_width * delta ) then
+                tm_mask(i,j,k) = 0
+                ftmp(i,j,k) = f(i,j+1,k) - delta
+              end if
+            end if
+
+            if ( eh_mask(i,j,k-1) .ge. 0 ) then
+              if ( f(i,j,k-1) - delta .gt. -shell_width * delta ) then
+                tm_mask(i,j,k) = 0
+                ftmp(i,j,k) = f(i,j,k-1) - delta
+              end if
+            end if
+
+            if ( eh_mask(i,j,k+1) .ge. 0 ) then
+              if ( f(i,j,k+1) - delta .gt. -shell_width * delta ) then
+                tm_mask(i,j,k) = 0
+                ftmp(i,j,k) = f(i,j,k+1) - delta
+              end if
+            end if
+
+          end if
+        end do
+      end do
+    end do
 
-    tm_mask = eh_mask
-
-    where ( eh_mask(2:nx-1,2:ny-1,2:nz-1) .eq. -1 )
-      where ( eh_mask(1:nx-2,2:ny-1,2:nz-1) .ge. 0 )
-        where ( f(1:nx-2,2:ny-1,2:nz-1) - delta .gt. -shell_width * delta )
-          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
-          f(2:nx-1,2:ny-1,2:nz-1) = f(1:nx-2,2:ny-1,2:nz-1) - delta
-        end where
-      end where
-      where ( eh_mask(3:nx,2:ny-1,2:nz-1) .ge. 0 )
-        where ( f(3:nx,2:ny-1,2:nz-1) - delta .gt. -shell_width * delta )
-          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
-          f(2:nx-1,2:ny-1,2:nz-1) = f(3:nx,2:ny-1,2:nz-1) - delta
-        end where
-      end where
-      where ( eh_mask(2:nx-1,1:ny-2,2:nz-1) .ge. 0 )
-        where ( f(2:nx-1,1:ny-2,2:nz-1) - delta .gt. -shell_width * delta )
-          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
-          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,1:ny-2,2:nz-1) - delta
-        end where
-      end where
-      where ( eh_mask(2:nx-1,3:ny,2:nz-1) .ge. 0 )
-        where ( f(2:nx-1,3:ny,2:nz-1) - delta .gt. -shell_width * delta )
-          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
-          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,3:ny,2:nz-1) - delta
-        end where
-      end where
-      where ( eh_mask(2:nx-1,2:ny-1,1:nz-2) .ge. 0 )
-        where ( f(2:nx-1,2:ny-1,1:nz-2) - delta .gt. -shell_width * delta )
-          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
-          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,2:ny-1,1:nz-2) - delta
-        end where
-      end where
-      where ( eh_mask(2:nx-1,2:ny-1,3:nz) .ge. 0 )
-        where ( f(2:nx-1,2:ny-1,3:nz) - delta .gt. -shell_width * delta )
-          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
-          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,2:ny-1,3:nz) - delta
-        end where
-      end where
+    eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = tm_mask(ixl:ixr,jyl:jyr,kzl:kzr)
+    f(ixl:ixr,jyl:jyr,kzl:kzr) = ftmp(ixl:ixr,jyl:jyr,kzl:kzr)
+
+    where ( ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) .and. &
+            ( f(ixl:ixr,jyl:jyr,kzl:kzr) .lt. -shell_width*delta ) ) 
+      eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = -1
+      f(ixl:ixr,jyl:jyr,kzl:kzr) = zero
     end where
 
-    eh_mask = tm_mask
+!    where ( eh_mask(2:nx-1,2:ny-1,2:nz-1) .eq. -1 )
+!      where ( eh_mask(1:nx-2,2:ny-1,2:nz-1) .ge. 0 )
+!        where ( f(1:nx-2,2:ny-1,2:nz-1) - delta .gt. -shell_width * delta )
+!          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
+!          f(2:nx-1,2:ny-1,2:nz-1) = f(1:nx-2,2:ny-1,2:nz-1) - delta
+!        end where
+!      end where
+!      where ( eh_mask(3:nx,2:ny-1,2:nz-1) .ge. 0 )
+!        where ( f(3:nx,2:ny-1,2:nz-1) - delta .gt. -shell_width * delta )
+!          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
+!          f(2:nx-1,2:ny-1,2:nz-1) = f(3:nx,2:ny-1,2:nz-1) - delta
+!        end where
+!      end where
+!      where ( eh_mask(2:nx-1,1:ny-2,2:nz-1) .ge. 0 )
+!        where ( f(2:nx-1,1:ny-2,2:nz-1) - delta .gt. -shell_width * delta )
+!          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
+!          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,1:ny-2,2:nz-1) - delta
+!        end where
+!      end where
+!      where ( eh_mask(2:nx-1,3:ny,2:nz-1) .ge. 0 )
+!        where ( f(2:nx-1,3:ny,2:nz-1) - delta .gt. -shell_width * delta )
+!          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
+!          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,3:ny,2:nz-1) - delta
+!        end where
+!      end where
+!      where ( eh_mask(2:nx-1,2:ny-1,1:nz-2) .ge. 0 )
+!        where ( f(2:nx-1,2:ny-1,1:nz-2) - delta .gt. -shell_width * delta )
+!          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
+!          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,2:ny-1,1:nz-2) - delta
+!        end where
+!      end where
+!      where ( eh_mask(2:nx-1,2:ny-1,3:nz) .ge. 0 )
+!        where ( f(2:nx-1,2:ny-1,3:nz) - delta .gt. -shell_width * delta )
+!          tm_mask(2:nx-1,2:ny-1,2:nz-1) = 0
+!          f(2:nx-1,2:ny-1,2:nz-1) = f(2:nx-1,2:ny-1,3:nz) - delta
+!        end where
+!      end where
+!    end where
+!
+!    eh_mask = tm_mask
+
+!    print*,'After'
+!    print*,eh_mask(:,jyl,kzl)
+!    print*,f(:,jyl,kzl)
+!    print*
+  end if
+end subroutine EHFinder_SetMask1
+
+
+subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
+
+  use EHFinder_mod
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k
+  logical :: active
+
+  active = .false.
+
+  if ( mask_first ) active = .true.
+  if ( reparametrize_every .gt. 0 ) then
+    if ( mod(cctk_iteration,reparametrize_every) .eq. 0 ) active = .true.
+  end if
+
+  if ( active ) then
+
+    # include "include/physical_part.h"
     
-    tm_mask = 0
-
-    where ( eh_mask(2:nx-1,2:ny-1,2:nz-1) .ge. 0 )
-      where ( eh_mask(1:nx-2,2:ny-1,2:nz-1) .eq. -1 )
-        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(0)
-      end where
-      where ( eh_mask(3:nx,2:ny-1,2:nz-1) .eq. -1 )
-        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(1)
-      end where
-      where ( eh_mask(2:nx-1,1:ny-2,2:nz-1) .eq. -1 )
-        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(2)
-      end where
-      where ( eh_mask(2:nx-1,3:ny,2:nz-1) .eq. -1 )
-        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(3)
-      end where
-      where ( eh_mask(2:nx-1,2:ny-1,1:nz-2) .eq. -1 )
-        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(4)
-      end where
-      where ( eh_mask(2:nx-1,2:ny-1,3:nz) .eq. -1 )
-        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(5)
-      end where
-    end where
+    if ( ixl .eq. 1 ) ixl = 2
+    if ( ixr .eq. nx ) ixr = nx - 1
+    if ( jyl .eq. 1 ) jyl = 2
+    if ( jyr .eq. ny ) jyr = ny - 1
+    if ( kzl .eq. 1 ) kzl = 2
+    if ( kzr .eq. nz ) kzr = nz - 1
+!    if ( symx ) then
+!      im = 1
+!    else
+!      im = 0
+!    end if
+!    if ( symy ) then
+!      jm = 1
+!    else
+!      jm = 0
+!    end if
+!    if ( symz ) then
+!      km = 1
+!    else
+!      km = 0
+!    end if
 
-    where ( eh_mask(2:nx-1,2:ny-1,2:nz-1) .ge. 0 )
-      eh_mask(2:nx-1,2:ny-1,2:nz-1) = tm_mask(2:nx-1,2:ny-1,2:nz-1)
-    end where
+    tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) = 0
+
+!    print*,'Iteration number ',cctk_iteration
+!    print*,count(eh_mask(ixl+1-im:ixr-1,jyl+1-jm:jyr-1,kzl+1-km:kzr-1) .eq. 0)
+!    print*
+!    print*,'Mask 2'
+!    print*,ixl,ixr,jyl,jyr,kzl,kzr
+!    print*,eh_mask(:,jyl,kzl)
+!    print*,f(:,jyl,kzl)
+    do k = kzl, kzr
+      do j = jyl, jyr
+        do i = ixl, ixr
 
-    tm_mask = 0
-    if ( cctk_bbox(1) .eq. 1 ) then
-      tm_mask(1,:,:) = tm_mask(1,:,:) + ll(0)
-    else
-      tm_mask(1:ngx,:,:) = -2
-    end if
-    if ( cctk_bbox(2) .eq. 1 ) then
-      tm_mask(nx,:,:) = tm_mask(nx,:,:) + ll(1)
-    else
-     tm_mask(nx-ngx+1:nx,:,:) = -2
-    end if
-    if ( cctk_bbox(3) .eq. 1 ) then
-      tm_mask(:,1,:) = tm_mask(:,1,:) + ll(2)
-    else
-      tm_mask(:,1:ngy,:) = -2
-    end if
-    if ( cctk_bbox(4) .eq. 1 ) then
-      tm_mask(:,ny,:) = tm_mask(:,ny,:) + ll(3)
-    else
-      tm_mask(:,ny-ngy+1:ny,:) = -2
-    end if
-    if ( cctk_bbox(5) .eq. 1 ) then
-      tm_mask(:,:,1) = tm_mask(:,:,1) + ll(4)
-    else
-      tm_mask(:,:,1:ngz) = -2
-    end if
-    if ( cctk_bbox(6) .eq. 1 ) then
-      tm_mask(:,:,nz) = tm_mask(:,:,nz) + ll(5)
-    else
-      tm_mask(:,:,nz-ngz+1:nz) = -2
-    end if
-
-    where ( tm_mask .ne. 0 )
-      eh_mask = tm_mask
+          if ( eh_mask(i,j,k) .ge. 0 ) then
+
+            if ( eh_mask(i-1,j,k) .eq. -1 ) then
+              tm_mask(i,j,k) = tm_mask(i,j,k) + ll(0)
+            end if
+
+            if ( eh_mask(i+1,j,k) .eq. -1 ) then
+              tm_mask(i,j,k) = tm_mask(i,j,k) + ll(1)
+            end if
+
+            if ( eh_mask(i,j-1,k) .eq. -1 ) then
+              tm_mask(i,j,k) = tm_mask(i,j,k) + ll(2)
+            end if
+
+            if ( eh_mask(i,j+1,k) .eq. -1 ) then
+              tm_mask(i,j,k) = tm_mask(i,j,k) + ll(3)
+            end if
+
+            if ( eh_mask(i,j,k-1) .eq. -1 ) then
+              tm_mask(i,j,k) = tm_mask(i,j,k) + ll(4)
+            end if
+
+            if ( eh_mask(i,j,k+1) .eq. -1 ) then
+              tm_mask(i,j,k) = tm_mask(i,j,k) + ll(5)
+            end if
+
+          end if
+
+        end do
+      end do
+    end do
+
+    where ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 )
+      eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = &
+                    tm_mask(ixl:ixr,jyl:jyr,kzl:kzr)
     end where
-    if ( CCTK_Equals(domain,'octant') ) then
-      if ( cctk_lbnd(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
-      if ( cctk_lbnd(2) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
-      if ( cctk_lbnd(3) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
-    end if
+
+!    print*,'Iteration number ',cctk_iteration
+!    print*,count(eh_mask(ixl+1-im:ixr-1,jyl+1-jm:jyr-1,kzl+1-km:kzr-1) .eq. 0)
+!    print*
+
+!    eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = &
+!                    tm_mask(ixl:ixr,jyl:jyr,kzl:kzr)
+
+!    where ( eh_mask(2:nx-1,2:ny-1,2:nz-1) .ge. 0 )
+!      where ( eh_mask(1:nx-2,2:ny-1,2:nz-1) .eq. -1 )
+!        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(0)
+!      end where
+!      where ( eh_mask(3:nx,2:ny-1,2:nz-1) .eq. -1 )
+!        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(1)
+!      end where
+!      where ( eh_mask(2:nx-1,1:ny-2,2:nz-1) .eq. -1 )
+!        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(2)
+!      end where
+!      where ( eh_mask(2:nx-1,3:ny,2:nz-1) .eq. -1 )
+!        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(3)
+!      end where
+!      where ( eh_mask(2:nx-1,2:ny-1,1:nz-2) .eq. -1 )
+!        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(4)
+!      end where
+!      where ( eh_mask(2:nx-1,2:ny-1,3:nz) .eq. -1 )
+!        tm_mask(2:nx-1,2:ny-1,2:nz-1)=tm_mask(2:nx-1,2:ny-1,2:nz-1) + ll(5)
+!      end where
+!    end where
+!
+!    where ( eh_mask(2:nx-1,2:ny-1,2:nz-1) .ge. 0 )
+!      eh_mask(2:nx-1,2:ny-1,2:nz-1) = tm_mask(2:nx-1,2:ny-1,2:nz-1)
+!    end where
+!
+!    tm_mask = 0
+!    if ( cctk_bbox(1) .eq. 1 ) then
+!      tm_mask(1,:,:) = tm_mask(1,:,:) + ll(0)
+!    else
+!      tm_mask(1:ngx,:,:) = -2
+!    end if
+!    if ( cctk_bbox(2) .eq. 1 ) then
+!      tm_mask(nx,:,:) = tm_mask(nx,:,:) + ll(1)
+!    else
+!     tm_mask(nx-ngx+1:nx,:,:) = -2
+!    end if
+!    if ( cctk_bbox(3) .eq. 1 ) then
+!      tm_mask(:,1,:) = tm_mask(:,1,:) + ll(2)
+!    else
+!      tm_mask(:,1:ngy,:) = -2
+!    end if
+!    if ( cctk_bbox(4) .eq. 1 ) then
+!      tm_mask(:,ny,:) = tm_mask(:,ny,:) + ll(3)
+!    else
+!      tm_mask(:,ny-ngy+1:ny,:) = -2
+!    end if
+!    if ( cctk_bbox(5) .eq. 1 ) then
+!      tm_mask(:,:,1) = tm_mask(:,:,1) + ll(4)
+!    else
+!      tm_mask(:,:,1:ngz) = -2
+!    end if
+!    if ( cctk_bbox(6) .eq. 1 ) then
+!      tm_mask(:,:,nz) = tm_mask(:,:,nz) + ll(5)
+!    else
+!      tm_mask(:,:,nz-ngz+1:nz) = -2
+!    end if
+!
+!    where ( tm_mask .ne. 0 )
+!      eh_mask = tm_mask
+!    end where
+!    if ( CCTK_Equals(domain,'octant') ) then
+!      if ( cctk_lbnd(1) .eq. 0 ) eh_mask(1:ngx,:,:) = -2
+!      if ( cctk_lbnd(2) .eq. 0 ) eh_mask(:,1:ngy,:) = -2
+!      if ( cctk_lbnd(3) .eq. 0 ) eh_mask(:,:,1:ngz) = -2
+!    end if
 
     call CCTK_INFO('Mask modified')
 
+    mask_first = .false.
+
+!    print*,'After'
+!    print*,eh_mask(:,jyl,kzl)
+!    print*,f(:,jyl,kzl)
+!    print*
   end if
 
-end subroutine EHFinder_SetMask
+end subroutine EHFinder_SetMask2