Added parameters to separately turn off inner and outer excision and added

some comments.


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

1 files changed, 318 insertions, 307 deletions

diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90
index ea8fa57..f521ba2 100644
--- a/src/EHFinder_SetMask.F90
+++ b/src/EHFinder_SetMask.F90
@@ -31,16 +31,16 @@ subroutine EHFinder_MaskInit(CCTK_ARGUMENTS)
 
   CCTK_INT :: i, j, k
 
-  ! Initialise the mask to 0.
+! Initialise the mask to 0.
   eh_mask = 0
 
 ! Figure out the local dimensions of the grid excluding ghostzones and
 ! symmetry zones.
 # include "include/physical_part.h"
 
-  ! Check if the point is located at a physical outer boundary and set the
-  ! eh_mask appropriately. The array ll is defined in EHFinder_mod.F90
-  ! and contains the values ( 1, 2, 4, 8, 16, 32 ).
+! Check if the point is located at a physical outer boundary and set the
+! eh_mask appropriately. The array ll is defined in EHFinder_mod.F90
+! and contains the values ( 1, 2, 4, 8, 16, 32 ).
   if ( ( cctk_bbox(1) .eq. 1 ) .and. ( ixl .eq. 1 ) ) then
     eh_mask(1,:,:) = eh_mask(1,:,:) + ll(0)
   end if
@@ -64,6 +64,10 @@ subroutine EHFinder_MaskInit(CCTK_ARGUMENTS)
 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
+! excsion region and make sure that the mask is correct there.
+
 subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -82,8 +86,13 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
 
   active = .false.
 
+! If the mask has not been set before, we should do it now so 
+! we set active=.true.
+! mask_first is initialized to .true. in EHFinder_mod.F90
   if ( mask_first ) active = .true.
 
+! If re-parametrization has just been done, we should reset the mask so we
+! set active=.true.
   if ( CCTK_EQUALS(re_param_method,'approx') .or. &
        CCTK_EQUALS(re_param_method,'mixed') ) then
     if ( reparametrize_every_approx .gt. 0 ) then
@@ -99,365 +108,350 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS)
     end if
   end if
 
+! If the reparametrization was not undone...
   if ( active .and. .not.reparam_undone ) then
 
+! Get the minimum and maximum index excluding ghost and symmetry cells.
 # include "include/physical_part.h"
 
+!   Store the current mask and level set function
     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)
 
-    fimin_loc = -ex_value; fimax_loc = -ex_value
-    imin_loc = cctk_gsh; imax_loc = 0
-
-    do k = kzl, kzr
-      do j = jyl, jyr
-        do i = ixl, ixr
-          if ( eh_mask(i,j,k) .ge. 0 ) then
-            if ( f(i,j,k) .lt. shell_width * delta ) then
-              if ( i+cctk_lbnd(1) .le. imin_loc(1) ) then
-                imin_loc(1) = i+cctk_lbnd(1)
-              end if
-              if ( i+cctk_lbnd(1) .ge. imax_loc(1) ) then
-                imax_loc(1) = i+cctk_lbnd(1)
-              end if
-              if ( j+cctk_lbnd(2) .le. imin_loc(2) ) then
-                imin_loc(2) = j+cctk_lbnd(2)
-              end if
-              if ( j+cctk_lbnd(2) .ge. imax_loc(2) ) then
-                imax_loc(2) = j+cctk_lbnd(2)
-              end if
-              if ( k+cctk_lbnd(3) .le. imin_loc(3) ) then
-                imin_loc(3) = k+cctk_lbnd(3)
-              end if
-              if ( k+cctk_lbnd(3) .ge. imax_loc(3) ) then
-                imax_loc(3) = k+cctk_lbnd(3)
+!   Next we try to locate the minimum and maximum of the global indeces of
+!   the currently active cells giving us the smallest rectangular box
+!   containing all active cells.
+
+    if ( use_outer_excision .gt. 0 ) then
+!     First initialize some variables.
+!     fimin_loc, and fimax_loc are 3 element arrays that will contain the 
+!     minimum value of f on the boundaries of this rectangular box. They
+!     will be used to decide if the box should be changed. They are
+!     initialized to the negative of the value used in excised cells.
+
+!     The 3 element arrays imin_loc and imax_loc will contain the min and
+!     max global indices of the ractangular box. They are initialised to
+!     the maximum global index and 0, respectively.
+
+      fimin_loc = -ex_value; fimax_loc = -ex_value
+      imin_loc = cctk_gsh; imax_loc = 0
+
+!     Find all cells with f<shell_width*delta and find their minimum
+!     and maximum global cell index.
+      do k = kzl, kzr
+        do j = jyl, jyr
+          do i = ixl, ixr
+            if ( eh_mask(i,j,k) .ge. 0 ) then
+              if ( f(i,j,k) .lt. shell_width * delta ) then
+                if ( i+cctk_lbnd(1) .le. imin_loc(1) ) then
+                  imin_loc(1) = i+cctk_lbnd(1)
+                end if
+                if ( i+cctk_lbnd(1) .ge. imax_loc(1) ) then
+                  imax_loc(1) = i+cctk_lbnd(1)
+                end if
+                if ( j+cctk_lbnd(2) .le. imin_loc(2) ) then
+                  imin_loc(2) = j+cctk_lbnd(2)
+                end if
+                if ( j+cctk_lbnd(2) .ge. imax_loc(2) ) then
+                  imax_loc(2) = j+cctk_lbnd(2)
+                end if
+                if ( k+cctk_lbnd(3) .le. imin_loc(3) ) then
+                  imin_loc(3) = k+cctk_lbnd(3)
+                end if
+                if ( k+cctk_lbnd(3) .ge. imax_loc(3) ) then
+                  imax_loc(3) = k+cctk_lbnd(3)
+                end if
               end if
             end if
-          end if
+          end do
         end do
       end do
-    end do
-    print*, 'Local values'
-    print*, imin_loc
-    print*, imax_loc
-
-!    do k = kzl, kzr
-!      do j = jyl, jyr
-!        do i = ixl, ixr
-!          if ( eh_mask(i,j,k) .ge. 0 ) then
-!            if ( f(i,j,k) .lt. shell_width * delta ) then
-!              if ( x(i,j,k) .le. cmin_loc(1) ) then
-!                cmin_loc(1) = x(i,j,k)
-!                imin_loc(1) = i
-!              end if
-!              if ( x(i,j,k) .ge. cmax_loc(1) ) then
-!                cmax_loc(1) = x(i,j,k)
-!                imax_loc(1) = i
-!              end if
-!              if ( y(i,j,k) .le. cmin_loc(2) ) then
-!                cmin_loc(2) = y(i,j,k)
-!                imin_loc(2) = j
-!              end if
-!              if ( y(i,j,k) .ge. cmax_loc(2) ) then
-!                cmax_loc(2) = y(i,j,k)
-!                imax_loc(2) = j
-!              end if
-!              if ( z(i,j,k) .le. cmin_loc(3) ) then
-!                cmin_loc(3) = z(i,j,k)
-!                imin_loc(3) = k
-!              end if
-!              if ( z(i,j,k) .ge. cmax_loc(3) ) then
-!                cmax_loc(3) = z(i,j,k)
-!                imax_loc(3) = k
-!              end if
-!            end if
-!          end if
-!        end do
-!      end do
-!    end do
-
-    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
-                                   imin_loc, imin_glob, 3, CCTK_VARIABLE_INT )
-    if ( ierr .ne. zero ) then
-      call CCTK_WARN ( 0, 'Reduction of array imin_loc failed' )
-    end if
-
-    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, &
-                                   imax_loc, imax_glob, 3, CCTK_VARIABLE_INT )
-    if ( ierr .ne. zero ) then
-      call CCTK_WARN ( 0, 'Reduction of array imax_loc failed' )
-    end if
 
-    i1 = imin_glob(1)-cctk_lbnd(1)
-    i2 = imax_glob(1)-cctk_lbnd(1)
-    j1 = imin_glob(2)-cctk_lbnd(2)
-    j2 = imax_glob(2)-cctk_lbnd(2)
-    k1 = imin_glob(3)-cctk_lbnd(3)
-    k2 = imax_glob(3)-cctk_lbnd(3)
+!     Reduce over all processors to get the global indeces of the
+!     rectangular box containing all cells with f<shell_width*delta
+      call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
+                                     imin_loc, imin_glob, 3, CCTK_VARIABLE_INT )
+      if ( ierr .ne. zero ) then
+        call CCTK_WARN ( 0, 'Reduction of array imin_loc failed' )
+      end if
 
-    if ( ( 1 .le. i1 ) .and. ( i1 .le. nx ) ) then
-      fimin_loc(1) = minval ( f(i1,jyl:jyr,kzl:kzr) )
-    end if
-    if ( ( 1 .le. i2 ) .and. ( i2 .le. nx ) ) then
-      fimax_loc(1) = minval ( f(i2,jyl:jyr,kzl:kzr) )
-    end if
-    if ( ( 1 .le. j1 ) .and. ( j1 .le. ny ) ) then
-      fimin_loc(2) = minval ( f(ixl:ixr,j1,kzl:kzr) )
-    end if
-    if ( ( 1 .le. j2 ) .and. ( j2 .le. ny ) ) then
-      fimax_loc(2) = minval ( f(ixl:ixr,j2,kzl:kzr) )
-    end if
-    if ( ( 1 .le. k1 ) .and. ( k1 .le. nz ) ) then
-      fimin_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k1) )
-    end if
-    if ( ( 1 .le. k2 ) .and. ( k2 .le. nz ) ) then
-      fimax_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k2) )
-    end if
+      call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, &
+                                     imax_loc, imax_glob, 3, CCTK_VARIABLE_INT )
+      if ( ierr .ne. zero ) then
+        call CCTK_WARN ( 0, 'Reduction of array imax_loc failed' )
+      end if
 
-    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
-                                fimin_loc, fimin_glob, 3, CCTK_VARIABLE_REAL )
-    if ( ierr .ne. zero ) then
-      call CCTK_WARN ( 0, 'Reduction of array fimin_loc failed' )
-    end if
+!     Convert into local indeces. Note that these might be less than 1 or
+!     larger than the local grid size if the box does not overlap with the
+!     local grid.
+      i1 = imin_glob(1)-cctk_lbnd(1)
+      i2 = imax_glob(1)-cctk_lbnd(1)
+      j1 = imin_glob(2)-cctk_lbnd(2)
+      j2 = imax_glob(2)-cctk_lbnd(2)
+      k1 = imin_glob(3)-cctk_lbnd(3)
+      k2 = imax_glob(3)-cctk_lbnd(3)
+
+!     Find the minimum value of f on the various faces of the rectangular box
+!     if part of the face is present on the current grid.
+      if ( ( 1 .le. i1 ) .and. ( i1 .le. nx ) ) then
+        fimin_loc(1) = minval ( f(i1,jyl:jyr,kzl:kzr) )
+      end if
+      if ( ( 1 .le. i2 ) .and. ( i2 .le. nx ) ) then
+        fimax_loc(1) = minval ( f(i2,jyl:jyr,kzl:kzr) )
+      end if
+      if ( ( 1 .le. j1 ) .and. ( j1 .le. ny ) ) then
+        fimin_loc(2) = minval ( f(ixl:ixr,j1,kzl:kzr) )
+      end if
+      if ( ( 1 .le. j2 ) .and. ( j2 .le. ny ) ) then
+        fimax_loc(2) = minval ( f(ixl:ixr,j2,kzl:kzr) )
+      end if
+      if ( ( 1 .le. k1 ) .and. ( k1 .le. nz ) ) then
+        fimin_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k1) )
+      end if
+      if ( ( 1 .le. k2 ) .and. ( k2 .le. nz ) ) then
+        fimax_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k2) )
+      end if
 
-    call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
-                                fimax_loc, fimax_glob, 3, CCTK_VARIABLE_REAL )
-    if ( ierr .ne. zero ) then
-      call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' )
-    end if
+!     Reduce to get the minimum values of f on the faces of the box
+      call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
+                                  fimin_loc, fimin_glob, 3, CCTK_VARIABLE_REAL )
+      if ( ierr .ne. zero ) then
+        call CCTK_WARN ( 0, 'Reduction of array fimin_loc failed' )
+      end if
 
-!    imin_loc = imin_loc
-!    imax_loc = imax_loc
+      call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, &
+                                  fimax_loc, fimax_glob, 3, CCTK_VARIABLE_REAL )
+      if ( ierr .ne. zero ) then
+        call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' )
+      end if
 
-    print*, 'Active range'
-    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)
+!      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
            
-    do k = kzl, kzr
-      do j = jyl, jyr
-        do i = ixl, ixr
-
-          if ( ( eh_mask(i,j,k) .eq. -1 ) .and. ( f(i,j,k) .lt. 0 ) ) 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
+!   Now check and see if any interior excised points need to be activated.
+    if ( use_inner_excision .gt. 0 ) then
+      do k = kzl, kzr
+        do j = jyl, jyr
+          do i = ixl, ixr
+
+!           If an interior excised point has a non-excised neighbour where
+!           f-delta>-shell_width*delta then activate the point by setting
+!           the temporary mask to zero. The new value of f will be the value
+!           if f in its neighbour point - delta. Do this for all directions.
+
+            if ( ( eh_mask(i,j,k) .eq. -1 ) .and. ( f(i,j,k) .lt. 0 ) ) 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
-            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
+              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
-            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
+              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
-            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
+              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
-            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
+              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
 
-            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
+              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 if
+            end if
+          end do
         end do
       end do
-    end do
+    end if
 
-    do i = 1, 3
-      imin_n(i) = imin_glob(i)
-      imax_n(i) = imax_glob(i)
-      if ( fimin_glob(i) + delta .lt. shell_width * delta ) then
-        if ( imin_glob(i) .gt. 1 ) then
-          imin_n(i) = imin_glob(i) - 1
-        endif
-      end if
-      if ( fimax_glob(i) + delta .lt. shell_width * delta ) then
-        if ( imax_glob(i) .lt. cctk_gsh(i) ) then
-          imax_n(i) = imax_glob(i) + 1
-        endif
-      end if
-    end do
+!   Check and see if the boundary of the exterior excision region should
+!   be changed and if so find the indices describing the new excision region.
+    if ( use_outer_excision .gt. 0 ) then
+      do i = 1, 3
+        imin_n(i) = imin_glob(i)
+        imax_n(i) = imax_glob(i)
+!       If the minimum value of f on a face on the box plus delta is less
+!       than shell_width * delta then the active region has to be increased
+!       in size in the corresponding region, i.e. inactive cells has to be
+!       activated.
+        if ( fimin_glob(i) + delta .lt. shell_width * delta ) then
+          if ( imin_glob(i) .gt. 1 ) then
+            imin_n(i) = imin_glob(i) - 1
+          endif
+        end if
+        if ( fimax_glob(i) + delta .lt. shell_width * delta ) then
+          if ( imax_glob(i) .lt. cctk_gsh(i) ) then
+            imax_n(i) = imax_glob(i) + 1
+          endif
+        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)
-
-    if ( imin_n(1) .ne. imin_glob(1) ) then
-      i1 = imin_n(1) - cctk_lbnd(1)
-      if ( ( ixl .le. i1 ) .and. ( i1 .lt. nx-1 ) ) then
-        j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
-        j2 = min ( imax_n(2) - cctk_lbnd(2), jyr )
-        k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
-        k2 = min ( imax_n(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) = 0
-          ftmp(i1,j1:j2,k1:k2) = ftmp(i1+1,j1:j2,k1:k2) + delta
-          print*,'done'
+!      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.
+      if ( imin_n(1) .ne. imin_glob(1) ) then
+        i1 = imin_n(1) - cctk_lbnd(1)
+        if ( ( ixl .le. i1 ) .and. ( i1 .lt. nx-1 ) ) then
+          j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
+          j2 = min ( imax_n(2) - cctk_lbnd(2), jyr )
+          k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
+          k2 = min ( imax_n(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) = 0
+            ftmp(i1,j1:j2,k1:k2) = ftmp(i1+1,j1:j2,k1:k2) + delta
+!            print*,'done'
+          end if
         end if
       end if
-    end if
-    if ( imax_n(1) .ne. imax_glob(1) ) then
-      i2 = imax_n(1) - cctk_lbnd(1)
-      if ( ( ixl .lt. i2 ) .and. ( i2 .le. nx-1 ) ) then
-        j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
-        j2 = min ( imax_n(2) - cctk_lbnd(2), jyr )
-        k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
-        k2 = min ( imax_n(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) = 0
-          ftmp(i2,j1:j2,k1:k2) = ftmp(i2-1,j1:j2,k1:k2) + delta
-          print*,'done'
+      if ( imax_n(1) .ne. imax_glob(1) ) then
+        i2 = imax_n(1) - cctk_lbnd(1)
+        if ( ( ixl .lt. i2 ) .and. ( i2 .le. nx-1 ) ) then
+          j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
+          j2 = min ( imax_n(2) - cctk_lbnd(2), jyr )
+          k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
+          k2 = min ( imax_n(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) = 0
+            ftmp(i2,j1:j2,k1:k2) = ftmp(i2-1,j1:j2,k1:k2) + delta
+!            print*,'done'
+          end if
         end if
       end if
-    end if
-    if ( imin_n(2) .ne. imin_glob(2) ) then
-      j1 = imin_n(2) - cctk_lbnd(2)
-      if ( ( jyl .le. j1 ) .and. ( j1 .lt. ny-1 ) ) then
-        i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
-        i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
-        k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
-        k2 = min ( imax_n(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) = 0
-          ftmp(i1:i2,j1,k1:k2) = ftmp(i1:i2,j1+1,k1:k2) + delta
-          print*,'done'
+      if ( imin_n(2) .ne. imin_glob(2) ) then
+        j1 = imin_n(2) - cctk_lbnd(2)
+        if ( ( jyl .le. j1 ) .and. ( j1 .lt. ny-1 ) ) then
+          i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
+          i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
+          k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
+          k2 = min ( imax_n(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) = 0
+            ftmp(i1:i2,j1,k1:k2) = ftmp(i1:i2,j1+1,k1:k2) + delta
+!            print*,'done'
+          end if
         end if
       end if
-    end if
-    if ( imax_n(2) .ne. imax_glob(2) ) then
-      j2 = imax_n(2) - cctk_lbnd(2)
-      if ( ( jyl .lt. j2 ) .and. ( j2 .le. ny-1 ) ) then
-        i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
-        i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
-        k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
-        k2 = min ( imax_n(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) = 0
-          ftmp(i1:i2,j2,k1:k2) = ftmp(i1:i2,j2-1,k1:k2) + delta
-          print*,'done'
+      if ( imax_n(2) .ne. imax_glob(2) ) then
+        j2 = imax_n(2) - cctk_lbnd(2)
+        if ( ( jyl .lt. j2 ) .and. ( j2 .le. ny-1 ) ) then
+          i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
+          i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
+          k1 = max ( imin_n(3) - cctk_lbnd(3), kzl )
+          k2 = min ( imax_n(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) = 0
+            ftmp(i1:i2,j2,k1:k2) = ftmp(i1:i2,j2-1,k1:k2) + delta
+!            print*,'done'
+          end if
         end if
       end if
-    end if
-    if ( imin_n(3) .ne. imin_glob(3) ) then
-      k1 = imin_n(3) - cctk_lbnd(3)
-      if ( ( kzl .le. k1 ) .and. ( k1 .lt. nz-1 ) ) then
-        i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
-        i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
-        j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
-        j2 = min ( imax_n(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) = 0
-          ftmp(i1:i2,j1:j2,k1) = ftmp(i1:i2,j1:j2,k1+1) + delta
-          print*,'done'
+      if ( imin_n(3) .ne. imin_glob(3) ) then
+        k1 = imin_n(3) - cctk_lbnd(3)
+        if ( ( kzl .le. k1 ) .and. ( k1 .lt. nz-1 ) ) then
+          i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
+          i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
+          j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
+          j2 = min ( imax_n(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) = 0
+            ftmp(i1:i2,j1:j2,k1) = ftmp(i1:i2,j1:j2,k1+1) + delta
+!            print*,'done'
+          end if
         end if
       end if
-    end if
-    if ( imax_n(3) .ne. imax_glob(3) ) then
-      k2 = imax_n(3) - cctk_lbnd(3)
-      if ( ( kzl .lt. k2 ) .and. ( k2 .le. nz-1 ) ) then
-        i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
-        i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
-        j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
-        j2 = min ( imax_n(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) = 0
-          ftmp(i1:i2,j1:j2,k2) = ftmp(i1:i2,j1:j2,k2-1) + delta
-          print*,'done'
+      if ( imax_n(3) .ne. imax_glob(3) ) then
+        k2 = imax_n(3) - cctk_lbnd(3)
+        if ( ( kzl .lt. k2 ) .and. ( k2 .le. nz-1 ) ) then
+          i1 = max ( imin_n(1) - cctk_lbnd(1), ixl )
+          i2 = min ( imax_n(1) - cctk_lbnd(1), ixr )
+          j1 = max ( imin_n(2) - cctk_lbnd(2), jyl )
+          j2 = min ( imax_n(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) = 0
+            ftmp(i1:i2,j1:j2,k2) = ftmp(i1:i2,j1:j2,k2-1) + delta
+!            print*,'done'
+          end if
         end if
       end if
-    end if
-
-!    if ( imax_n(1) .ne. imax_loc(1) ) then
-!        tm_mask(imax_n(1),imin_n(2):imax_n(2),imin_n(3):imax_n(3)) = 0
-!        ftmp(imax_n(1),imin_n(2):imax_n(2),imin_n(3):imax_n(3)) = &
-!          ftmp(imax_n(1)-1,imin_n(2):imax_n(2),imin_n(3):imax_n(3)) + delta
-!    end if
-!    if ( imin_n(2) .ne. imin_loc(2) ) then
-!        tm_mask(imin_n(1):imax_n(1),imin_n(2),imin_n(3):imax_n(3)) = 0
-!        ftmp(imin_n(1):imax_n(1),imin_n(2),imin_n(3):imax_n(3)) = &
-!          ftmp(imin_n(1):imax_n(1),imin_n(2)+1,imin_n(3):imax_n(3)) + delta
-!    end if
-!    if ( imax_n(2) .ne. imax_loc(2) ) then
-!        tm_mask(imin_n(1):imax_n(1),imax_n(2),imin_n(3):imax_n(3)) = 0
-!        ftmp(imin_n(1):imax_n(1),imax_n(2),imin_n(3):imax_n(3)) = &
-!          ftmp(imin_n(1):imax_n(1),imax_n(2)-1,imin_n(3):imax_n(3)) + delta
-!    end if
-!    if ( imin_n(3) .ne. imin_loc(3) ) then
-!        tm_mask(imin_n(1):imax_n(1),imin_n(2):imax_n(2),imin_n(3)) = 0
-!        ftmp(imin_n(1):imax_n(1),imin_n(2):imax_n(2),imin_n(3)) = &
-!          ftmp(imin_n(1):imax_n(1),imin_n(2):imax_n(2),imin_n(3)+1) + delta
-!    end if
-!    if ( imax_n(3) .ne. imax_loc(3) ) then
-!        tm_mask(imin_n(1):imax_n(1),imin_n(2):imax_n(2),imax_n(3)) = 0
-!        ftmp(imin_n(1):imax_n(1),imin_n(2):imax_n(2),imax_n(3)) = &
-!          ftmp(imin_n(1):imax_n(1),imin_n(2):imax_n(2),imax_n(3)-1) + delta
-!    end if
 
+    end if
 
+!   Copy the modified mask and level set function into the proper place.
     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) = ex_value
-    end where
+!   Now check if any interior points are far enough away from the f=0
+!   surface and if so excise them.
+    if ( use_inner_excision .gt. 0 ) then
+      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) = ex_value
+      end where
+    end if
 
-    do k = kzl, kzr
-      do j = jyl, jyr
-        do i = ixl, ixr
-          if ( ( ( i + cctk_lbnd(1) .lt. imin_n(1) ) .or. &
-                 ( i + cctk_lbnd(1) .gt. imax_n(1) ) .or. &
-                 ( j + cctk_lbnd(2) .lt. imin_n(2) ) .or. &
-                 ( j + cctk_lbnd(2) .gt. imax_n(2) ) .or. &
-                 ( k + cctk_lbnd(3) .lt. imin_n(3) ) .or. &
-                 ( k + cctk_lbnd(3) .gt. imax_n(3) ) ) .and. &
-                 ( eh_mask(i,j,k) .ge. 0 ) ) then
-             eh_mask(i,j,k) = -1
-             f(i,j,k) = -ex_value
-           end if
-         end do
-       end do
-     end do
-!    where ( ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) .and. &
-!            ( f(ixl:ixr,jyl:jyr,kzl:kzr) .gt. shell_width*delta ) ) 
-!      eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = -1
-!      f(ixl:ixr,jyl:jyr,kzl:kzr) = -ex_value
-!    end where
+!   Make sure to mark all points outside of the rectangular box as excised
+!   points and set the value of f to -ex_value.
+    if ( use_outer_excision .gt. 0 ) then
+      do k = kzl, kzr
+        do j = jyl, jyr
+          do i = ixl, ixr
+            if ( ( ( i + cctk_lbnd(1) .lt. imin_n(1) ) .or. &
+                   ( i + cctk_lbnd(1) .gt. imax_n(1) ) .or. &
+                   ( j + cctk_lbnd(2) .lt. imin_n(2) ) .or. &
+                   ( j + cctk_lbnd(2) .gt. imax_n(2) ) .or. &
+                     ( k + cctk_lbnd(3) .lt. imin_n(3) ) .or. &
+                   ( k + cctk_lbnd(3) .gt. imax_n(3) ) ) .and. &
+                   ( eh_mask(i,j,k) .ge. 0 ) ) then
+              eh_mask(i,j,k) = -1
+              f(i,j,k) = -ex_value
+            end if
+          end do
+        end do
+      end do
+    end if
   end if
 end subroutine EHFinder_SetMask1
 
 
+! This routine finds the excision boundary after it has been changed and
+! makes sure that it has the right mask value.
 subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
 
   use EHFinder_mod
@@ -473,6 +467,8 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
 
   active = .false.
 
+! If re-parametrization has just been done, we are in the process of
+! resetting the mask so we set active=.true.
   if ( mask_first ) active = .true.
   if ( CCTK_EQUALS(re_param_method,'approx') .or. &
        CCTK_EQUALS(re_param_method,'mixed') ) then
@@ -489,10 +485,13 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
     end if
   end if
 
+! If the reparametrization was not undone...
   if ( active .and. .not.reparam_undone ) then
 
+! Get the minimum and maximum index excluding ghost and symmetry cells.
 # include "include/physical_part.h"
     
+!   Make sure we are not on the physical outer boundary.
     if ( ixl .eq. 1 ) ixl = 2
     if ( ixr .eq. nx ) ixr = nx - 1
     if ( jyl .eq. 1 ) jyl = 2
@@ -500,34 +499,44 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
     if ( kzl .eq. 1 ) kzl = 2
     if ( kzr .eq. nz ) kzr = nz - 1
 
+!   Initialize the temporary mask to zero.
     tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) = 0
 
+!   We loop over all points and adjust the mask if the point is
+!   on the excision boundary.
     do k = kzl, kzr
       do j = jyl, jyr
         do i = ixl, ixr
 
+!         If the point is active.....
           if ( eh_mask(i,j,k) .ge. 0 ) then
 
+!           If the neighbour in the -x-directions is excised....
             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 the neighbour in the +x-directions is excised....
             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 the neighbour in the -y-directions is excised....
             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 the neighbour in the +y-directions is excised....
             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 the neighbour in the -z-directions is excised....
             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 the neighbour in the +z-directions is excised....
             if ( eh_mask(i,j,k+1) .eq. -1 ) then
               tm_mask(i,j,k) = tm_mask(i,j,k) + ll(5)
             end if
@@ -538,6 +547,7 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
       end do
     end do
 
+!   Copy the changes back into eh_mask
     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)
@@ -545,6 +555,7 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS)
 
     call CCTK_INFO('Mask modified')
 
+!   Indicate that it is not anymore the first time the mask is set.
     mask_first = .false.
   end if