Introduce a better treatment of the shocktube boundaries.

* The 2d diagonal boundary conditions are correctly applied
when using multiprocessors.

* The ymin (ymax) surfaces have their points in the region
indsum < minsum (indsum > maxsum) fixed after a first sync
operation. A second sync guarantees the neighbour grid
receives the correct value of the boundary on those regions.

* Mirror symmetries applied at the faces intercepting the upper
right corner assumed cubic local domain. However, rectangular
local domains are common after domain decomposition. This commit
fixes the indexes applying this symmetry for these cases.

* Upper right corner/edge/surface condition should use sten instead
of stenp1. Thanks for catching this Josh!

* Attempt at a 3d diagonal boundary treatment (Josh's commit)

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@135 ac85fae7-cede-4708-beff-ae01c7fa1c26

3 files changed, 2370 insertions, 113 deletions

diff --git a/par/balsara1d.par b/par/balsara1d.par
index 3d1980f..f688b16 100644
--- a/par/balsara1d.par
+++ b/par/balsara1d.par
@@ -81,7 +81,7 @@ Cactus::terminate="iteration"
 #Cactus::cctk_final_time = 0.02
 #Cactus::cctk_final_time = 0.4
 #cactus::cctk_itlast = 40
-cactus::cctk_itlast = 400
+cactus::cctk_itlast = 800
 #cactus::cctk_itlast = 2
 
 IO::out_dir = $parfile
@@ -97,12 +97,12 @@ CarpetIOBasic::outInfo_every=1
 ##CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel grhydro::psidc grhydro::dens grhydro::tau grhydro::scon HydroBase::Bvec"
 ##CarpetIOASCII::out1D_vars = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel  grhydro::dens grhydro::tau grhydro::scon "
 
-#CarpetIOHDF5::out_every = 1
-#CarpetIOHDF5::out_vars  = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel HydroBase::w_lorentz  grhydro::dens grhydro::divB grhydro::tau grhydro::scon HydroBase::Bvec"
+CarpetIOHDF5::out_every = 10
+CarpetIOHDF5::out_vars  = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel HydroBase::w_lorentz  grhydro::dens grhydro::divB grhydro::tau grhydro::scon HydroBase::Bvec"
 ##CarpetIOHDF5::out_vars  = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel HydroBase::w_lorentz grhydro::psidc grhydro::dens grhydro::divB grhydro::tau grhydro::scon HydroBase::Bvec"
 ##CarpetIOHDF5::out_vars  = "HydroBase::rho HydroBase::press HydroBase::eps  HydroBase::vel HydroBase::w_lorentz  grhydro::dens grhydro::tau grhydro::scon "
 
-CarpetIOHDF5::out2D_every                    = 1
+CarpetIOHDF5::out2D_every                    = 10
 CarpetIOHDF5::out2D_xy                       = "yes"
 CarpetIOHDF5::out2D_xz                       = "no"
 CarpetIOHDF5::out2D_yz                       = "no"
diff --git a/schedule.ccl b/schedule.ccl
index 2e32506..2b72ab4 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -43,20 +43,46 @@ if (CCTK_Equals(initial_hydro,"shocktube")) {
 
     if(CCTK_Equals(shocktube_type,"diagshock")) 
     {
-      schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
-      {
-        LANG: Fortran	
-      } "Diagonal shock boundary conditions"
+      if(clean_divergence){
+        schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+        {
+          LANG: Fortran	
+          SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons, GRHydro::psidc
+        } "Diagonal shock boundary conditions"
+      } else {
+        schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+        {
+          LANG: Fortran	
+          SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons
+        } "Diagonal shock boundary conditions"
+      }
     }
 
     if(CCTK_Equals(shocktube_type,"diagshock2d")) 
     {
-      schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
-      {
-        LANG: Fortran	
-      } "2-D Diagonal shock boundary conditions"
+        schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+        {
+          LANG: Fortran	
+        } "2-D Diagonal shock boundary conditions"
     }
 
+#    if(CCTK_Equals(shocktube_type,"diagshock2d")) 
+#    {
+#      if(clean_divergence){
+#        schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+#        {
+#          LANG: Fortran	
+#          SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons, GRHydro::psidc
+#        } "2-D Diagonal shock boundary conditions"
+#      } else {
+#        schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection
+#        {
+#          LANG: Fortran	
+#          SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons
+#        } "2-D Diagonal shock boundary conditions"
+#      }
+#    }
+
   } else {
     schedule GRHydro_shocktube in HydroBase_Initial
     {
diff --git a/src/GRHydro_ShockTubeM.F90 b/src/GRHydro_ShockTubeM.F90
index 6418baf..fc6a45b 100644
--- a/src/GRHydro_ShockTubeM.F90
+++ b/src/GRHydro_ShockTubeM.F90
@@ -57,7 +57,7 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
   DECLARE_CCTK_FUNCTIONS
   
   CCTK_INT :: i, j, k, nx, ny, nz
-  CCTK_REAL :: direction, det
+  CCTK_REAL :: direction, det, DIRECTION_TINY
   CCTK_REAL :: rhol, rhor, velxl, velxr, velyl, velyr, &
        velzl, velzr, epsl, epsr
   CCTK_REAL :: bvcxl,bvcyl,bvczl,bvcxr,bvcyr,bvczr
@@ -531,8 +531,10 @@ subroutine GRHydro_shocktubeM(CCTK_ARGUMENTS)
           call CCTK_WARN(0,"Shock case not recognized")
         end if
 
-        if ( ((change_shock_direction==0).and.(direction .lt. 0.0d0)).or.& 
-             ((change_shock_direction==1).and.(direction .gt. 0.0d0)) ) then
+        DIRECTION_TINY = 0.0001*CCTK_DELTA_SPACE(1)
+
+        if ( ((change_shock_direction==0).and.(direction .lt. DIRECTION_TINY)).or.& 
+             ((change_shock_direction==1).and.(direction .gt. DIRECTION_TINY)) ) then
 
 !!$ Left state
 
@@ -651,10 +653,18 @@ subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
   
-  CCTK_INT :: i, j, k, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew
+  CCTK_INT :: i, j, k, nx, ny, nz, sten, stenp1, minsum, maxsum
+  CCTK_INT :: inew, jnew, knew
   CCTK_INT :: xoff,yoff,zoff,indsum
   CCTK_REAL :: det
+  CCTK_INT :: status
+  character(len=200) :: error_message
+  integer :: gindex_tau, gindex_dens, gindex_scon, gindex_Bcons
+  CCTK_INT, parameter :: num_groups = 4
+  CCTK_INT :: groups(num_groups)
+
 
+  sten=GRHydro_stencil
   stenp1=GRHydro_stencil + 1
 
   nx = cctk_lsh(1)
@@ -668,42 +678,879 @@ subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS)
   yoff=0
   zoff=0
         
-  do k=1,nz
-     if(k.lt.stenp1)then
-              zoff=k-stenp1
-           else if(k.gt.nz-stenp1+1) then
-              zoff=k-(nz-stenp1+1)
-           else
-              zoff=0
-           endif
+!!$ Loop over the cubical domain 6 faces:
+
+  !!$ 1) xmin:
+  do k = stenp1, nz-sten
+    do j = stenp1, ny-sten
+      do i = 1, sten
+
+        indsum = i+j+k
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          knew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - stenp1 + 1
+          jnew = ny - stenp1 + 1
+          knew = nz - stenp1 + 1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+           call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+          dens(i,j,k) = dens(inew,jnew,knew)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,knew)
+        sy(i,j,k) = sy(inew,jnew,knew)
+        sz(i,j,k) = sz(inew,jnew,knew)
+        tau(i,j,k) = tau(inew,jnew,knew)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,knew)
+        endif
 
-     do j=1,ny
-        if(j.lt.stenp1)then
-           yoff=j-stenp1
-        else if(j.gt.ny-stenp1+1) then
-           yoff=j-(ny-stenp1+1)
+      end do 
+    end do
+  end do 
+
+  !!$ 2) xmax:
+  do k = stenp1, nz-sten
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx
+
+        indsum = i+j+k
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          knew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - stenp1 + 1
+          jnew = ny - stenp1 + 1
+          knew = nz - stenp1 + 1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
         else
-           yoff=0
+           call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+           dens(i,j,k) = dens(inew,jnew,knew)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,knew)
+        sy(i,j,k) = sy(inew,jnew,knew)
+        sz(i,j,k) = sz(inew,jnew,knew)
+        tau(i,j,k) = tau(inew,jnew,knew)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,knew)
         endif
 
-        do i=1,nx
-           if(i.lt.stenp1) then
-              xoff=i-stenp1
-           else if(i.gt.nx-stenp1+1) then
-              xoff=i-(nx-stenp1+1)
-           else
-              xoff=0
-           endif
+      end do 
+    end do
+  end do 
+
+  !!$ 3) ymin:
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i =  stenp1, nx-sten
+
+        indsum = i+j+k
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          knew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - stenp1 + 1
+          jnew = ny - stenp1 + 1
+          knew = nz - stenp1 + 1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+           call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+           dens(i,j,k) = dens(inew,jnew,knew)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,knew)
+        sy(i,j,k) = sy(inew,jnew,knew)
+        sz(i,j,k) = sz(inew,jnew,knew)
+        tau(i,j,k) = tau(inew,jnew,knew)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,knew)
+        endif
 
-           indsum = i+j+k
+      end do 
+    end do
+  end do 
+
+  !!$ 4) ymax:
+  do k = stenp1, nz-sten
+    do j = ny-sten+1, ny
+      do i = stenp1, nx-sten
+
+        indsum = i+j+k
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          knew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - stenp1 + 1
+          jnew = ny - stenp1 + 1
+          knew = nz - stenp1 + 1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+           call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+           dens(i,j,k) = dens(inew,jnew,knew)
+        end if
            
-           if( (xoff.ne.0.or.yoff.ne.0.or.zoff.ne.0) .and. &
-                indsum.ge.minsum.and.indsum.le.maxsum) then
-!!$ We can map the point to the interior diagonal, orthogonal to the shock.
+        sx(i,j,k) = sx(inew,jnew,knew)
+        sy(i,j,k) = sy(inew,jnew,knew)
+        sz(i,j,k) = sz(inew,jnew,knew)
+        tau(i,j,k) = tau(inew,jnew,knew)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,knew)
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$  the following two cases are different for a 3d diagonal than a 2d diagonal
+
+  !!$ 5) zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten
+      do i =  stenp1, nx-sten
+ 
+        indsum = i+j+k
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          knew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - stenp1 + 1
+          jnew = ny - stenp1 + 1
+          knew = nz - stenp1 + 1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+           call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+           dens(i,j,k) = dens(inew,jnew,knew)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,knew)
+        sy(i,j,k) = sy(inew,jnew,knew)
+        sz(i,j,k) = sz(inew,jnew,knew)
+        tau(i,j,k) = tau(inew,jnew,knew)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,knew)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 6) zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = stenp1, nx-sten
+ 
+        indsum = i+j+k
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          knew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - stenp1 + 1
+          jnew = ny - stenp1 + 1
+          knew = nz - stenp1 + 1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+           call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+           dens(i,j,k) = dens(inew,jnew,knew)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,knew)
+        sy(i,j,k) = sy(inew,jnew,knew)
+        sz(i,j,k) = sz(inew,jnew,knew)
+        tau(i,j,k) = tau(inew,jnew,knew)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,knew)
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Loop over the cubical domain 12 edges:
+
+  !!$ 1) xmin,ymin:
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i = 1, sten
+
+         indsum = i+j+k
+!!$ indsum > maxsum impossible!
+         if(indsum.lt.minsum) then
+            inew = stenp1
+            jnew = stenp1
+            knew = stenp1
+         else
+            call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         end if
+
+!!$  Don't poison edges!
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 2) xmin,ymax:
+  do k = stenp1, nz-sten
+    do j =  ny-sten+1, ny
+      do i = 1, sten
+
+         indsum = i+j+k
+!!$ indsum > maxsum and indsum < minsum both impossible!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 3) xmax,ymin:
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i = nx-sten+1, nx 
+
+         indsum = i+j+k
+!!$ indsum > maxsum and indsum < minsum both impossible!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 4) xmax,ymax:
+  do k = stenp1, nz-sten
+    do j = ny-sten+1, ny 
+      do i = nx-sten+1, nx 
+
+         indsum = i+j+k
+!!$ indsum < minsum impossible!
+         if(indsum.gt.maxsum) then
+            inew = nx - stenp1 + 1
+            jnew = ny - stenp1 + 1
+            knew = nz - stenp1 + 1
+         else
+            call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         end if
+         
+!!$  Don't poison edges!
+         dens(i,j,k) = dens(inew,jnew,knew)
+
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 5) ymin,zmin:
+  do k = 1, sten 
+    do j = 1, sten 
+      do i =  stenp1, nx-sten
+
+         indsum = i+j+k
+!!$ indsum > maxsum impossible!
+         if(indsum.lt.minsum) then
+            inew = stenp1
+            jnew = stenp1
+            knew = stenp1
+         else
+            call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         end if
+         
+!!$  Don't poison edges!
+         dens(i,j,k) = dens(inew,jnew,knew)
+
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+ 
+      end do 
+    end do
+  end do 
+
+  !!$ 6) ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten 
+      do i =  stenp1, nx-sten
+ 
+         indsum = i+j+k
+!!$ indsum > maxsum and indsum < minsum both impossible!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 7) ymax,zmin:
+  do k = 1, sten 
+    do j = ny-sten+1, ny
+      do i =  stenp1, nx-sten
+ 
+         indsum = i+j+k
+!!$ indsum > maxsum and indsum < minsum both impossible!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 8) ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny
+      do i =  stenp1, nx-sten
+ 
+         indsum = i+j+k
+!!$ indsum < minsum impossible!
+         if(indsum.gt.maxsum) then
+            inew = nx - stenp1 + 1
+            jnew = ny - stenp1 + 1
+            knew = nz - stenp1 + 1
+         else
+            call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         end if
+
+!!$  Don't poison edges!
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 9) xmin,zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten 
+      do i = 1, sten 
+ 
+         indsum = i+j+k
+!!$ indsum > maxsum impossible!
+         if(indsum.lt.minsum) then
+            inew = stenp1
+            jnew = stenp1
+            knew = stenp1
+         else
+            call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         end if
+         
+!!$  Don't poison edges!
+         dens(i,j,k) = dens(inew,jnew,knew)
+
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 10) xmin,zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = 1, sten 
+
+         indsum = i+j+k
+!!$ indsum > maxsum and indsum < minsum both impossible!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 11) xmax,zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx 
+ 
+         indsum = i+j+k
+!!$ indsum > maxsum and indsum < minsum both impossible!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 12) xmax,zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx 
+
+         indsum = i+j+k
+!!$ indsum < minsum impossible!
+         if(indsum.gt.maxsum) then
+            inew = nx - stenp1 + 1
+            jnew = ny - stenp1 + 1
+            knew = nz - stenp1 + 1
+         else
+            call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         end if
+                  
+!!$  Don't poison edges!
+         dens(i,j,k) = dens(inew,jnew,knew)
+
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Loop over the cubical domain 8 corners:
+
+  !!$ 1) xmin,ymin,zmin:
+  do k = 1, sten
+    do j = 1, sten
+      do i = 1, sten
+ 
+!!$ indsum < minsum guaranteed!
+
+         inew = stenp1
+         jnew = stenp1
+         knew = stenp1
+
+!!$  Don't poison corners!
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 2) xmin,ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten
+      do i = 1, sten
+ 
+         indsum = i+j+k
+!!$ indsum < maxsum and indsum > minsum guaranteed!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 3) xmin,ymax,zmin:
+  do k = 1, sten
+    do j = ny-sten+1, ny 
+      do i = 1, sten
+ 
+         indsum = i+j+k
+!!$ indsum < maxsum and indsum > minsum guaranteed!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 4) xmin,ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny
+      do i = 1, sten
+ 
+         indsum = i+j+k
+!!$ indsum < maxsum and indsum > minsum guaranteed!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 5) xmax,ymin,zmin:
+  do k = 1, sten
+    do j = 1, sten
+      do i = nx-sten+1, nx
+ 
+         indsum = i+j+k
+!!$ indsum < maxsum and indsum > minsum guaranteed!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 6) xmax,ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten
+      do i = nx-sten+1, nx
+ 
+         indsum = i+j+k
+!!$ indsum < maxsum and indsum > minsum guaranteed!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 7) xmax,ymax,zmin:
+  do k = 1, sten
+    do j = ny-sten+1, ny 
+      do i = nx-sten+1, nx
+ 
+         indsum = i+j+k
+!!$ indsum < maxsum and indsum > minsum guaranteed!
+         call find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+         dens(i,j,k) = dens(inew,jnew,knew)         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 8) xmax,ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny 
+      do i = nx-sten+1, nx
+ 
+!!$ indsum > maxsum guaranteed!
+
+         inew = nx - stenp1 + 1
+         jnew = ny - stenp1 + 1
+         knew = nz - stenp1 + 1
+
+!!$  Don't poison corners!
+         dens(i,j,k) = dens(inew,jnew,knew)
+         
+         sx(i,j,k) = sx(inew,jnew,knew)
+         sy(i,j,k) = sy(inew,jnew,knew)
+         sz(i,j,k) = sz(inew,jnew,knew)
+         tau(i,j,k) = tau(inew,jnew,knew)
+         Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+         Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+         Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+         if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,knew)
+         endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Synchronize ghost zones (first time):
+
+  call CCTK_GroupIndex(gindex_tau,"GRHydro::tau")
+  call CCTK_GroupIndex(gindex_dens,"GRHydro::dens")
+  call CCTK_GroupIndex(gindex_scon,"GRHydro::scon")
+  call CCTK_GroupIndex(gindex_Bcons,"GRHydro::Bcons")
+
+  groups(1) = gindex_tau
+  groups(2) = gindex_dens
+  groups(3) = gindex_scon
+  groups(4) = gindex_Bcons
+
+  call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups)
+  if(status.lt.0) then
+    write(error_message,'(a33,i3,a5,i3,a7)' ) &
+    'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.'
+    call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+  endif
+  
+  if(clean_divergence.ne.0) then
+    call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc")
+    if(status.lt.0) then
+      write(error_message,'(a31,i3,a19)' ) &
+      'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc'
+      call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+    endif
+  endif 
+
+
+!!$ Fix the indsum < minsum area in one go?
+  do k=1,minsum-3
+     do j=1,minsum-k-2
+        do i=1,minsum-k-j-1
+
+           indsum=i+j+k
+           if(dens(i,j,k).lt.0.0d0) then
+              if(dens(j,i,k).gt.0.0d0) then
+                 inew=j
+                 jnew=i
+                 knew=k
+              else if(dens(k,j,i).gt.0.0d0) then
+                 inew=k
+                 jnew=j
+                 knew=i
+              else if(dens(i,k,j).gt.0.0d0) then
+                 inew=i
+                 jnew=k
+                 knew=j
+              else
+                 inew=stenp1
+                 jnew=stenp1
+                 knew=stenp1
+              endif
 
-              inew=indsum/3
-              jnew=(indsum-inew)/2
-              knew=indsum-inew-jnew
+              dens(i,j,k) = dens(inew,jnew,knew)
+              sx(i,j,k) = sx(inew,jnew,knew)
+              sy(i,j,k) = sy(inew,jnew,knew)
+              sz(i,j,k) = sz(inew,jnew,knew)
+              tau(i,j,k) = tau(inew,jnew,knew)
+              Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+              Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+              Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+              if(clean_divergence.ne.0) then 
+                 psidc(i,j,k)=psidc(inew,jnew,knew)
+              endif
+           end if
+
+        end do
+     end do
+  end do
+
+!!$ Fix the indsum > maxsum area in one go?
+  do k=nz-3*stenp1+4,nz
+     do j=ny-3*stenp1+4+(nz-k),ny
+        do i=nx-3*stenp1+4+(nz-k)+(ny-j),nx
+
+           indsum=i+j+k
+           xoff=nx-i
+           yoff=ny-j
+           zoff=nz-k
+           if(dens(i,j,k).lt.0.0d0) then
+
+              if(dens(nx-yoff,ny-xoff,k).gt.0.0d0) then
+                 inew=nx-yoff
+                 jnew=ny-xoff
+                 knew=k
+              else if(dens(nx-zoff,j,nz-xoff).gt.0.0d0) then
+                 inew=nx-zoff
+                 jnew=j
+                 knew=nz-xoff
+              else if(dens(i,ny-zoff,nz-yoff).gt.0.0d0) then
+                 inew=i
+                 jnew=ny-zoff
+                 knew=nz-xoff
+              else
+                 inew=nx - stenp1 + 1
+                 jnew=ny - stenp1 + 1
+                 knew=nz - stenp1 + 1
+              endif
               
               dens(i,j,k) = dens(inew,jnew,knew)
               sx(i,j,k) = sx(inew,jnew,knew)
@@ -713,16 +1560,214 @@ subroutine GRHydro_Diagshock_BoundaryM(CCTK_ARGUMENTS)
               Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
               Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
               Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
-             if(clean_divergence.ne.0) then 
-               psidc(i,j,k)=psidc(inew,jnew,knew)
-             endif
-
-           endif
+              if(clean_divergence.ne.0) then 
+                 psidc(i,j,k)=psidc(inew,jnew,knew)
+              endif
+           end if
            
-        enddo
-     enddo
-  enddo
+        end do
+     end do
+  end do
   
+
+!!$ Fix the ymin face in the region indsum < minsum
+!!$  do k = stenp1, nz-sten
+!!$    do j = 1, sten
+!!$      do i =  stenp1, nx-sten
+!!$
+!!$        indsum = i+j
+!!$        if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then
+!!$          inew = j
+!!$          jnew = i
+!!$          if(dens(inew,jnew,k).lt.0.0d0) then
+!!$            inew = stenp1
+!!$            jnew = stenp1
+!!$          endif 
+!!$           
+!!$          dens(i,j,k) = dens(inew,jnew,k)
+!!$          sx(i,j,k) = sx(inew,jnew,k)
+!!$          sy(i,j,k) = sy(inew,jnew,k)
+!!$          sz(i,j,k) = sz(inew,jnew,k)
+!!$          tau(i,j,k) = tau(inew,jnew,k)
+!!$          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+!!$          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+!!$          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+!!$          if(clean_divergence.ne.0) then 
+!!$            psidc(i,j,k)=psidc(inew,jnew,k)
+!!$          endif
+!!$        end if
+!!$
+!!$      end do 
+!!$    end do
+!!$  end do 
+
+!!$ Fix the xmin face in the region indsum < minsum
+!!$  do k = stenp1, nz-sten
+!!$    do j = stenp1, ny-sten
+!!$      do i = 1, sten
+!!$
+!!$        indsum = i+j
+!!$        if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then
+!!$          inew = j
+!!$          jnew = i
+!!$          if(dens(inew,jnew,k).lt.0.0d0) then
+!!$            inew = stenp1
+!!$            jnew = stenp1
+!!$          endif 
+!!$           
+!!$          dens(i,j,k) = dens(inew,jnew,k)
+!!$          sx(i,j,k) = sx(inew,jnew,k)
+!!$          sy(i,j,k) = sy(inew,jnew,k)
+!!$          sz(i,j,k) = sz(inew,jnew,k)
+!!$          tau(i,j,k) = tau(inew,jnew,k)
+!!$          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+!!$          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+!!$          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+!!$          if(clean_divergence.ne.0) then 
+!!$            psidc(i,j,k)=psidc(inew,jnew,k)
+!!$          endif
+!!$        end if
+!!$
+!!$
+!!$      end do 
+!!$    end do
+!!$  end do 
+
+!!$ Fix the ymax face in the region indsum > maxsum
+!!$  do k = stenp1, nz-sten
+!!$    do j = ny-sten+1, ny
+!!$      do i = stenp1, nx-sten
+!!$
+!!$        indsum = i+j
+!!$        if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then
+!!$          inew = j-ny+nx
+!!$          jnew = i-nx+ny
+!!$          if(dens(inew,jnew,k).lt.0.0d0) then
+!!$            inew = nx - stenp1
+!!$            jnew = ny - stenp1
+!!$          endif 
+!!$           
+!!$          dens(i,j,k) = dens(inew,jnew,k)
+!!$          sx(i,j,k) = sx(inew,jnew,k)
+!!$          sy(i,j,k) = sy(inew,jnew,k)
+!!$          sz(i,j,k) = sz(inew,jnew,k)
+!!$          tau(i,j,k) = tau(inew,jnew,k)
+!!$          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+!!$          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+!!$          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+!!$          if(clean_divergence.ne.0) then 
+!!$            psidc(i,j,k)=psidc(inew,jnew,k)
+!!$          endif
+!!$        endif
+!!$
+!!$      end do 
+!!$    end do
+!!$  end do 
+
+!!$ Fix the xmax face in the region indsum > maxsum
+!!$  do k = stenp1, nz-sten
+!!$    do j = stenp1, ny-sten
+!!$      do i = nx-sten+1, nx
+!!$
+!!$        indsum = i+j
+!!$        if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then
+!!$          inew = j-ny+nx
+!!$          jnew = i-nx+ny
+!!$          if(dens(inew,jnew,k).lt.0.0d0) then
+!!$            inew = nx - stenp1
+!!$            jnew = ny - stenp1
+!!$          endif 
+!!$           
+!!$          dens(i,j,k) = dens(inew,jnew,k)
+!!$          sx(i,j,k) = sx(inew,jnew,k)
+!!$          sy(i,j,k) = sy(inew,jnew,k)
+!!$          sz(i,j,k) = sz(inew,jnew,k)
+!!$          tau(i,j,k) = tau(inew,jnew,k)
+!!$          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+!!$          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+!!$          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+!!$          if(clean_divergence.ne.0) then 
+!!$            psidc(i,j,k)=psidc(inew,jnew,k)
+!!$          endif
+!!$        endif
+!!$
+!!$      end do 
+!!$    end do
+!!$  end do 
+
+!!$ Synchronize ghost zones (second time):
+
+  call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups)
+  if(status.lt.0) then
+    write(error_message,'(a33,i3,a5,i3,a7)' ) &
+    'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.'
+    call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+  endif
+  
+  if(clean_divergence.ne.0) then
+    call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc")
+    if(status.lt.0) then
+      write(error_message,'(a31,i3,a19)' ) &
+      'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc'
+      call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+    endif
+  endif 
+
+!!$  do k=1,nz
+!!$     if(k.lt.stenp1)then
+!!$              zoff=k-stenp1
+!!$           else if(k.gt.nz-stenp1+1) then
+!!$              zoff=k-(nz-stenp1+1)
+!!$           else
+!!$              zoff=0
+!!$           endif
+!!$
+!!$     do j=1,ny
+!!$        if(j.lt.stenp1)then
+!!$           yoff=j-stenp1
+!!$        else if(j.gt.ny-stenp1+1) then
+!!$           yoff=j-(ny-stenp1+1)
+!!$        else
+!!$           yoff=0
+!!$        endif
+!!$
+!!$        do i=1,nx
+!!$           if(i.lt.stenp1) then
+!!$              xoff=i-stenp1
+!!$           else if(i.gt.nx-stenp1+1) then
+!!$              xoff=i-(nx-stenp1+1)
+!!$           else
+!!$              xoff=0
+!!$           endif
+!!$
+!!$           indsum = i+j+k
+!!$           
+!!$           if( (xoff.ne.0.or.yoff.ne.0.or.zoff.ne.0) .and. &
+!!$                indsum.ge.minsum.and.indsum.le.maxsum) then
+!!$    We can map the point to the interior diagonal, orthogonal to the shock.
+!!$
+!!$              inew=indsum/3
+!!$              jnew=(indsum-inew)/2
+!!$              knew=indsum-inew-jnew
+!!$              
+!!$              dens(i,j,k) = dens(inew,jnew,knew)
+!!$              sx(i,j,k) = sx(inew,jnew,knew)
+!!$              sy(i,j,k) = sy(inew,jnew,knew)
+!!$              sz(i,j,k) = sz(inew,jnew,knew)
+!!$              tau(i,j,k) = tau(inew,jnew,knew)
+!!$              Bconsx(i,j,k)=Bconsx(inew,jnew,knew)
+!!$              Bconsy(i,j,k)=Bconsy(inew,jnew,knew)
+!!$              Bconsz(i,j,k)=Bconsz(inew,jnew,knew)
+!!$             if(clean_divergence.ne.0) then 
+!!$               psidc(i,j,k)=psidc(inew,jnew,knew)
+!!$             endif
+!!$
+!!$           endif
+!!$           
+!!$        enddo
+!!$     enddo
+!!$  enddo
+ 
 end subroutine GRHydro_Diagshock_BoundaryM
 
 subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS)
@@ -733,10 +1778,16 @@ subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS
   DECLARE_CCTK_FUNCTIONS
   
-  CCTK_INT :: i, j, k, kc, nx, ny, nz, stenp1, minsum, maxsum, inew, jnew, knew
-  CCTK_INT :: xoff,yoff,zoff,indsum
+  CCTK_INT :: i, j, k, kc, nx, ny, nz, sten, stenp1, minsum, maxsum, inew, jnew, knew
+  CCTK_INT :: xoff,yoff,zoff,indsum,numoff
   CCTK_REAL :: det
+  CCTK_INT :: status
+  character(len=200) :: error_message
+  integer :: gindex_tau, gindex_dens, gindex_scon, gindex_Bcons
+  CCTK_INT, parameter :: num_groups = 4
+  CCTK_INT :: groups(num_groups)
   
+  sten=GRHydro_stencil
   stenp1=GRHydro_stencil + 1
   
   nx = cctk_lsh(1)
@@ -744,80 +1795,1260 @@ subroutine GRHydro_Diagshock2D_BoundaryM(CCTK_ARGUMENTS)
   nz = cctk_lsh(3)
   
   minsum = 2*stenp1
-  maxsum = nx+ny-2*stenp1 + 2
+  maxsum = nx+ny-2*sten
   
   xoff=0
   yoff=0
   zoff=0
-        
-  do k=1,nz
-     if(k.lt.stenp1)then
-              zoff=k-stenp1
-           else if(k.gt.nz-stenp1+1) then
-              zoff=k-(nz-stenp1+1)
-           else
-              zoff=0
-           endif
 
-     do j=1,ny
-        if(j.lt.stenp1)then
-           yoff=j-stenp1
-        else if(j.gt.ny-stenp1+1) then
-           yoff=j-(ny-stenp1+1)
+!!$ Loop over the cubical domain 6 faces:
+
+  !!$ 1) xmin:
+  do k = stenp1, nz-sten
+    do j = stenp1, ny-sten
+      do i = 1, sten
+
+        xoff=i-stenp1
+        yoff=0
+        zoff=0
+        indsum = i+j
+        inew = i - xoff
+        jnew = j + xoff
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - sten
+          jnew = ny - sten
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
         else
-           yoff=0
+          dens(i,j,k) = dens(inew,jnew,k)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
         endif
 
-        do i=1,nx
-           if(i.lt.stenp1) then
-              xoff=i-stenp1
-           else if(i.gt.nx-stenp1+1) then
-              xoff=i-(nx-stenp1+1)
-           else
-              xoff=0
-           endif
-  
+      end do 
+    end do
+  end do 
+
+  !!$ 2) xmax:
+  do k = stenp1, nz-sten
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx
+
+        xoff=i-(nx-stenp1+1)
+        yoff=0
+        zoff=0
         indsum = i+j
-        
-        if( (xoff.ne.0.or.yoff.ne.0) .and. &
-             indsum.ge.minsum.and.indsum.le.maxsum) then
-!!$ We can map the point to the interior diagonal, orthogonal to the shock.
+        inew = i - xoff
+        jnew = j + xoff
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - sten
+          jnew = ny - sten
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+          dens(i,j,k) = dens(inew,jnew,k)
+        end if
            
-           inew=indsum/2
-           jnew=indsum-inew
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 3) ymin:
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i =  stenp1, nx-sten
+
+        xoff=0
+        yoff=j-stenp1
+        zoff=0
+        indsum = i+j
+        inew = i + yoff
+        jnew = j - yoff
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - sten
+          jnew = ny - sten
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+          dens(i,j,k) = dens(inew,jnew,k)
+        end if
            
-           dens(i,j,k) = dens(inew,jnew,k)
-           sx(i,j,k) = sx(inew,jnew,k)
-           sy(i,j,k) = sy(inew,jnew,k)
-           sz(i,j,k) = sz(inew,jnew,k)
-           tau(i,j,k) = tau(inew,jnew,k)
-           Bconsx(i,j,k)=Bconsx(inew,jnew,k)
-           Bconsy(i,j,k)=Bconsy(inew,jnew,k)
-           Bconsz(i,j,k)=Bconsz(inew,jnew,k)
-           if(clean_divergence.ne.0) then 
-              psidc(i,j,k)=psidc(inew,jnew,k)
-           endif
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
 
-        else if( zoff.ne.0) then
+  !!$ 4) ymax:
+  do k = stenp1, nz-sten
+    do j = ny-sten+1, ny
+      do i = stenp1, nx-sten
+
+        xoff=0
+        yoff=j-(ny-stenp1+1)
+        zoff=0
+        indsum = i+j
+        inew = i + yoff
+        jnew = j - yoff
+        if(indsum.lt.minsum) then
+          inew = stenp1
+          jnew = stenp1
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else if(indsum.gt.maxsum) then
+          inew = nx - sten
+          jnew = ny - sten
+          dens(i,j,k) = -1.0d0 !!$ poison at shock propagation direction corners
+        else
+          dens(i,j,k) = dens(inew,jnew,k)
+        end if
+           
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 5) zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten
+      do i =  stenp1, nx-sten
  
-           kc = nz/2
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
 
-           dens(i,j,k) = dens(i,j,kc)
-           sx(i,j,k) = sx(i,j,kc)
-           sy(i,j,k) = sy(i,j,kc)
-           sz(i,j,k) = sz(i,j,kc)
-           tau(i,j,k) = tau(i,j,kc)
-           Bconsx(i,j,k)=Bconsx(i,j,kc)
-           Bconsy(i,j,k)=Bconsy(i,j,kc)
-           Bconsz(i,j,k)=Bconsz(i,j,kc)
-           if(clean_divergence.ne.0) then
-              psidc(i,j,k)=psidc(i,j,kc)
-           endif
+      end do 
+    end do
+  end do 
+
+  !!$ 6) zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Loop over the cubical domain 12 edges:
+
+  !!$ 1) xmin,ymin:
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i = 1, sten
+
+        !!$ We already know that for this edge indsum.lt.minsum.
+        inew = stenp1
+        jnew = stenp1
+           
+        dens(i,j,k) = dens(inew,jnew,k)
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 2) xmin,ymax:
+  do k = stenp1, nz-sten
+    do j =  ny-sten+1, ny
+      do i = 1, sten
+
+        xoff=i-stenp1
+        yoff=j-(ny-stenp1+1)
+        zoff=0
+ 
+        numoff = max(abs(xoff),abs(yoff))
+        if( abs(xoff).eq.numoff) then
+          inew = i + numoff
+          jnew = j - numoff
+        else
+          jnew = j - numoff
+          inew = i + numoff
+        endif  
+
+        dens(i,j,k) = dens(inew,jnew,k)
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 3) xmax,ymin:
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i = nx-sten+1, nx 
+
+        xoff=i-(nx-stenp1+1)
+        yoff=j-stenp1
+        zoff=0
+        numoff = max(abs(xoff),abs(yoff))
+        if( abs(xoff).eq.numoff) then
+          inew = i - numoff
+          jnew = j + numoff
+        else
+          jnew = j + numoff
+          inew = i - numoff
+        endif  
+           
+        dens(i,j,k) = dens(inew,jnew,k)
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 4) xmax,ymax:
+  do k = stenp1, nz-sten
+    do j = ny-sten+1, ny 
+      do i = nx-sten+1, nx 
+
+        !!$ We already know that for this edge indsum.gt.maxsum
+        inew = nx - sten
+        jnew = ny - sten
+          
+        dens(i,j,k) = dens(inew,jnew,k)
+        sx(i,j,k) = sx(inew,jnew,k)
+        sy(i,j,k) = sy(inew,jnew,k)
+        sz(i,j,k) = sz(inew,jnew,k)
+        tau(i,j,k) = tau(inew,jnew,k)
+        Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+        Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+        Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+        if(clean_divergence.ne.0) then 
+          psidc(i,j,k)=psidc(inew,jnew,k)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 5) ymin,zmin:
+  do k = 1, sten 
+    do j = 1, sten 
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 6) ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten 
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 7) ymax,zmin:
+  do k = 1, sten 
+    do j = ny-sten+1, ny
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 8) ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 9) xmin,zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten 
+      do i = 1, sten 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 10) xmin,zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = 1, sten 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 11) xmax,zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 12) xmax,zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Loop over the cubical domain 8 corners:
+
+  !!$ 1) xmin,ymin,zmin:
+  do k = 1, sten
+    do j = 1, sten
+      do i = 1, sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 2) xmin,ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten
+      do i = 1, sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 3) xmin,ymax,zmin:
+  do k = 1, sten
+    do j = ny-sten+1, ny 
+      do i = 1, sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 4) xmin,ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny
+      do i = 1, sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 5) xmax,ymin,zmin:
+  do k = 1, sten
+    do j = 1, sten
+      do i = nx-sten+1, nx
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 6) xmax,ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten
+      do i = nx-sten+1, nx
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 7) xmax,ymax,zmin:
+  do k = 1, sten
+    do j = ny-sten+1, ny 
+      do i = nx-sten+1, nx
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 8) xmax,ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny 
+      do i = nx-sten+1, nx
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Synchronize ghost zones (first time):
+
+  call CCTK_GroupIndex(gindex_tau,"GRHydro::tau")
+  call CCTK_GroupIndex(gindex_dens,"GRHydro::dens")
+  call CCTK_GroupIndex(gindex_scon,"GRHydro::scon")
+  call CCTK_GroupIndex(gindex_Bcons,"GRHydro::Bcons")
+
+  groups(1) = gindex_tau
+  groups(2) = gindex_dens
+  groups(3) = gindex_scon
+  groups(4) = gindex_Bcons
+
+  call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups)
+  if(status.lt.0) then
+    write(error_message,'(a33,i3,a5,i3,a7)' ) &
+    'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.'
+    call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+  endif
+  
+  if(clean_divergence.ne.0) then
+    call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc")
+    if(status.lt.0) then
+      write(error_message,'(a31,i3,a19)' ) &
+      'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc'
+      call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+    endif
+  endif 
+
+!!$ Fix the ymin face in the region indsum < minsum
+  do k = stenp1, nz-sten
+    do j = 1, sten
+      do i =  stenp1, nx-sten
+
+        indsum = i+j
+        if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then
+          inew = j
+          jnew = i
+          if(dens(inew,jnew,k).lt.0.0d0) then
+            inew = stenp1
+            jnew = stenp1
+          endif 
+           
+          dens(i,j,k) = dens(inew,jnew,k)
+          sx(i,j,k) = sx(inew,jnew,k)
+          sy(i,j,k) = sy(inew,jnew,k)
+          sz(i,j,k) = sz(inew,jnew,k)
+          tau(i,j,k) = tau(inew,jnew,k)
+          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+          if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,k)
+          endif
+        end if
+
+      end do 
+    end do
+  end do 
+
+!!$ Fix the xmin face in the region indsum < minsum
+  do k = stenp1, nz-sten
+    do j = stenp1, ny-sten
+      do i = 1, sten
+
+        indsum = i+j
+        if(indsum.lt.minsum.and.dens(i,j,k).lt.0.0d0) then
+          inew = j
+          jnew = i
+          if(dens(inew,jnew,k).lt.0.0d0) then
+            inew = stenp1
+            jnew = stenp1
+          endif 
+           
+          dens(i,j,k) = dens(inew,jnew,k)
+          sx(i,j,k) = sx(inew,jnew,k)
+          sy(i,j,k) = sy(inew,jnew,k)
+          sz(i,j,k) = sz(inew,jnew,k)
+          tau(i,j,k) = tau(inew,jnew,k)
+          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+          if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,k)
+          endif
+        end if
+
+
+      end do 
+    end do
+  end do 
+
+!!$ Fix the ymax face in the region indsum > maxsum
+  do k = stenp1, nz-sten
+    do j = ny-sten+1, ny
+      do i = stenp1, nx-sten
+
+        indsum = i+j
+        if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then
+          inew = j-ny+nx
+          jnew = i-nx+ny
+          if(dens(inew,jnew,k).lt.0.0d0) then
+            inew = nx - sten
+            jnew = ny - sten
+          endif 
+           
+          dens(i,j,k) = dens(inew,jnew,k)
+          sx(i,j,k) = sx(inew,jnew,k)
+          sy(i,j,k) = sy(inew,jnew,k)
+          sz(i,j,k) = sz(inew,jnew,k)
+          tau(i,j,k) = tau(inew,jnew,k)
+          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+          if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,k)
+          endif
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Fix the xmax face in the region indsum > maxsum
+  do k = stenp1, nz-sten
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx
+
+        indsum = i+j
+        if(indsum.gt.maxsum.and.dens(i,j,k).lt.0.0d0) then
+          inew = j-ny+nx
+          jnew = i-nx+ny
+          if(dens(inew,jnew,k).lt.0.0d0) then
+            inew = nx - sten
+            jnew = ny - sten
+          endif 
            
+          dens(i,j,k) = dens(inew,jnew,k)
+          sx(i,j,k) = sx(inew,jnew,k)
+          sy(i,j,k) = sy(inew,jnew,k)
+          sz(i,j,k) = sz(inew,jnew,k)
+          tau(i,j,k) = tau(inew,jnew,k)
+          Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+          Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+          Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+          if(clean_divergence.ne.0) then 
+            psidc(i,j,k)=psidc(inew,jnew,k)
+          endif
+        endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Fix the upper and lower (5-12 below) domain edges since 
+!!$ they still contain poison got from the faces in the initial
+!!$ step: 
+
+  !!$ 5) ymin,zmin:
+  do k = 1, sten 
+    do j = 1, sten 
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 6) ymin,zmax:
+  do k = nz-sten+1, nz
+    do j = 1, sten 
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 7) ymax,zmin:
+  do k = 1, sten 
+    do j = ny-sten+1, ny
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 8) ymax,zmax:
+  do k = nz-sten+1, nz
+    do j = ny-sten+1, ny
+      do i =  stenp1, nx-sten
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 9) xmin,zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten 
+      do i = 1, sten 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 10) xmin,zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = 1, sten 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 11) xmax,zmin:
+  do k = 1, sten 
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
+        endif
+
+      end do 
+    end do
+  end do 
+
+  !!$ 12) xmax,zmax:
+  do k = nz-sten+1, nz
+    do j = stenp1, ny-sten
+      do i = nx-sten+1, nx 
+ 
+        kc = nz/2
+
+        dens(i,j,k) = dens(i,j,kc)
+        sx(i,j,k) = sx(i,j,kc)
+        sy(i,j,k) = sy(i,j,kc)
+        sz(i,j,k) = sz(i,j,kc)
+        tau(i,j,k) = tau(i,j,kc)
+        Bconsx(i,j,k)=Bconsx(i,j,kc)
+        Bconsy(i,j,k)=Bconsy(i,j,kc)
+        Bconsz(i,j,k)=Bconsz(i,j,kc)
+        if(clean_divergence.ne.0) then
+          psidc(i,j,k)=psidc(i,j,kc)
         endif
+
+      end do 
+    end do
+  end do 
+
+!!$ Synchronize ghost zones (second time):
+
+  call CCTK_SyncGroupsI(status,cctkGH,num_groups,groups)
+  if(status.lt.0) then
+    write(error_message,'(a33,i3,a5,i3,a7)' ) &
+    'CCTK_SyncGroupsI returned status ',status,' for ',num_groups,'groups.'
+    call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+  endif
+  
+  if(clean_divergence.ne.0) then
+    call CCTK_SyncGroup(status,cctkGH,"GRHydro::psidc")
+    if(status.lt.0) then
+      write(error_message,'(a31,i3,a19)' ) &
+      'CCTK_SyncGroup returned status ',status,' for GRHydro::psidc'
+      call CCTK_WARN(CCTK_WARN_ABORT,error_message)
+    endif
+  endif 
         
-        enddo
-     enddo
-  enddo
+!!$  do k=1,nz
+!!$     if(k.lt.stenp1)then
+!!$        zoff=k-stenp1
+!!$     else if(k.gt.nz-stenp1+1) then
+!!$        zoff=k-(nz-stenp1+1)
+!!$     else
+!!$        zoff=0
+!!$     endif
+!!$
+!!$     do j=1,ny
+!!$        if(j.lt.stenp1)then
+!!$           yoff=j-stenp1
+!!$        else if(j.gt.ny-stenp1+1) then
+!!$           yoff=j-(ny-stenp1+1)
+!!$        else
+!!$           yoff=0
+!!$        endif
+!!$
+!!$        do i=1,nx
+!!$           if(i.lt.stenp1) then
+!!$              xoff=i-stenp1
+!!$           else if(i.gt.nx-stenp1+1) then
+!!$              xoff=i-(nx-stenp1+1)
+!!$           else
+!!$              xoff=0
+!!$           endif
+!!$  
+!!$        indsum = i+j
+!!$        
+!!$        if(xoff.ne.0.or.yoff.ne.0) then
+!!$
+!!$          if(indsum.ge.minsum.and.indsum.le.maxsum) then
+!!$
+!!$           numoff = max(abs(xoff),abs(yoff))
+!!$           if( abs(xoff).eq.numoff) then
+!!$             if(xoff.gt.0) then
+!!$               inew = i - numoff
+!!$               jnew = j + numoff
+!!$             else
+!!$               inew = i + numoff
+!!$               jnew = j - numoff
+!!$             endif
+!!$           else
+!!$             if(yoff.gt.0) then
+!!$               jnew = j - numoff
+!!$               inew = i + numoff
+!!$             else
+!!$               jnew = j + numoff
+!!$               inew = i - numoff
+!!$             endif
+!!$           endif  
+!!$!!$ Note that the following two conditions helps but
+!!$!!$ don't solve the problem when running in several
+!!$!!$ processors
+!!$         else if(indsum.lt.minsum) then
+!!$           inew = stenp1
+!!$           jnew = stenp1
+!!$
+!!$         else if(indsum.gt.maxsum) then
+!!$           inew = nx - stenp1
+!!$           jnew = ny - stenp1
+!!$         end if
+!!$           
+!!$           dens(i,j,k) = dens(inew,jnew,k)
+!!$           sx(i,j,k) = sx(inew,jnew,k)
+!!$           sy(i,j,k) = sy(inew,jnew,k)
+!!$           sz(i,j,k) = sz(inew,jnew,k)
+!!$           tau(i,j,k) = tau(inew,jnew,k)
+!!$           Bconsx(i,j,k)=Bconsx(inew,jnew,k)
+!!$           Bconsy(i,j,k)=Bconsy(inew,jnew,k)
+!!$           Bconsz(i,j,k)=Bconsz(inew,jnew,k)
+!!$           if(clean_divergence.ne.0) then 
+!!$              psidc(i,j,k)=psidc(inew,jnew,k)
+!!$           endif
+!!$
+!!$        else if( zoff.ne.0) then
+!!$ 
+!!$           kc = nz/2
+!!$
+!!$           dens(i,j,k) = dens(i,j,kc)
+!!$           sx(i,j,k) = sx(i,j,kc)
+!!$           sy(i,j,k) = sy(i,j,kc)
+!!$           sz(i,j,k) = sz(i,j,kc)
+!!$           tau(i,j,k) = tau(i,j,kc)
+!!$           Bconsx(i,j,k)=Bconsx(i,j,kc)
+!!$           Bconsy(i,j,k)=Bconsy(i,j,kc)
+!!$           Bconsz(i,j,k)=Bconsz(i,j,kc)
+!!$           if(clean_divergence.ne.0) then
+!!$              psidc(i,j,k)=psidc(i,j,kc)
+!!$           endif
+!!$           
+!!$        endif
+!!$        
+!!$        enddo
+!!$     enddo
+!!$  enddo
   
 end subroutine GRHydro_Diagshock2D_BoundaryM
+
+subroutine find_mapped_point(i,j,k,nx,ny,nz,stenp1,inew,jnew,knew)
+
+!!$ Find the nearest mapped point in the interior, assuming it exists
+
+  implicit none
+  CCTK_INT :: i,j,k,n,nx,ny,nz,stenp1
+  CCTK_INT :: xoff,yoff,zoff,indsum,minsum,maxsum
+  CCTK_INT :: inew,jnew,knew,dspace,newsum
+  
+  minsum = 3*stenp1
+  maxsum = nx+ny+nz - 3*(stenp1-1)
+
+  indsum = i+j+k
+
+  if(indsum.lt.minsum.or.indsum.gt.maxsum)write(6,*)'HUGE ERROR - ILLEGAL INDEX SUM!!!!!!'
+
+  if(i.lt.stenp1)then
+     xoff = i - stenp1 
+  else if(i.gt.nx-stenp1+1) then
+     xoff=i - (nx-stenp1+1)
+  else
+     xoff=0
+  endif
+  
+  if(j.lt.stenp1)then
+     yoff = j - stenp1 
+  else if(j.gt.ny-stenp1+1) then
+     yoff=j - (ny-stenp1+1)
+  else
+     yoff=0
+  endif
+     
+  if(k.lt.stenp1)then
+     zoff = k - stenp1 
+  else if(k.gt.nz-stenp1+1) then
+     zoff=k-(nz-stenp1+1)
+  else
+     zoff=0
+  endif
+  
+  if(xoff.eq.0.and.yoff.eq.0.and.zoff.eq.0) then
+     write(6,*)'Why did you try to map an interior point?'
+     inew=i
+     jnew=j
+     knew=k
+     return
+  else
+     
+!!$ We have something to actually do!
+     
+     inew = i-xoff
+     jnew = j-yoff
+     knew = k-zoff
+     dspace = xoff + yoff + zoff
+     
+     do n=1,abs(dspace)
+        
+        if(dspace.lt.0) then
+           
+!!$     dspace negative; point sum corrected upward
+!!$     need to move other coords down
+           
+           if(inew.gt.stenp1) then
+              inew = inew - 1
+           else if(jnew.gt.stenp1) then
+              jnew = jnew - 1
+           else if (knew.gt.stenp1) then
+              knew = knew - 1
+           else
+              write(6,*)'Big error!!!!'
+           endif
+           
+        else if (dspace.gt.0) then
+           
+!!$     dspace positive; point sum corrected downward
+!!$     need to move other coords up
+           
+           if(inew.lt.nx-stenp1+1) then
+              inew=inew+1
+           else if(jnew.lt.ny-stenp1+1) then
+              jnew=jnew+1
+           else if(knew.lt.nz-stenp1+1) then
+              knew=knew+1
+           else
+              write(6,*)'Big error2!!!!'
+           endif
+
+        end if
+
+     end do
+     
+     newsum = inew+jnew+knew
+     
+     if(inew.lt.stenp1.or.inew.gt.(nx-(stenp1-1)).or. &
+          jnew.lt.stenp1.or.jnew.gt.(ny-(stenp1-1)).or. &
+          knew.lt.stenp1.or.knew.gt.(nz-(stenp1-1)).or. &
+          newsum.ne.indsum)write(6,*)'Big error 3!!!!'
+     
+  end if
+
+  return
+
+end subroutine find_mapped_point
+