Changes in methodology for filling in the interior of Cook-Pfeiffer data as

used in the "turducken" paper. This includes some spline blending from
Erik and the most recent method of filling in the interior using an elliptic
equation with the "good" points as boundary values. This method uses a
conjugate gradient solver (built in), that should be robust (i.e. it should
always work) but may not be the fastest. However, since this is only done on
initial data this shouldn't be an issue. The elliptic method currently support
a second order, fourth order and sixth order equation, that will give C0, C1
and C2 solutions, respectively, across the boundary of the fill in region.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/NoExcision/trunk@7 4ec1db94-0e4f-0410-ada3-8bed251432c9

10 files changed, 1603 insertions, 80 deletions

diff --git a/src/NoExcision_mod.F90 b/src/NoExcision_mod.F90
new file mode 100644
index 0000000..646d335
--- /dev/null
+++ b/src/NoExcision_mod.F90
@@ -0,0 +1,266 @@
+! $Header$
+
+#include "cctk.h"
+
+module NoExcision_mod
+
+  implicit none
+
+  CCTK_REAL, parameter :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  CCTK_REAL, parameter :: one = 1.0_wp
+  CCTK_REAL, parameter :: two = 2.0_wp
+  CCTK_INT :: infnorm_handle, sum_handle
+  CCTK_INT :: sym_selector
+  CCTK_INT :: loop_counter
+  CCTK_REAL, dimension(16) :: infnormresid
+  CCTK_REAL, dimension(16) :: delta_new, delta_old, alpha, beta, absres
+  character(len=18), dimension(16) :: red_names = (/ 'NoExcision::redgxx', &
+                                                     'NoExcision::redgxy', &
+                                                     'NoExcision::redgxz', &
+                                                     'NoExcision::redgyy', &
+                                                     'NoExcision::redgyz', &
+                                                     'NoExcision::redgzz', &
+                                                     'NoExcision::redkxx', &
+                                                     'NoExcision::redkxy', &
+                                                     'NoExcision::redkxz', &
+                                                     'NoExcision::redkyy', &
+                                                     'NoExcision::redkyz', &
+                                                     'NoExcision::redkzz', &
+                                                     'NoExcision::red', &
+                                                     'NoExcision::redx', &
+                                                     'NoExcision::redy', &
+                                                     'NoExcision::redz' /)
+  character(len=5), dimension(16) :: var_names = (/ '  gxx', '  gxy', &
+                                                    '  gxz', '  gyy', &
+                                                    '  gyz', '  gzz', &
+                                                    '  kxx', '  kxy', &
+                                                    '  kxz', '  kyy', &
+                                                    '  kyz', '  kzz', &
+                                                    '  alp', 'betax', &
+                                                    'betay', 'betaz' /)
+  integer :: nx, ny, nz
+  integer :: ii
+  logical, dimension(16) :: cont = (/ (.true., ii=1, 16) /)
+
+  contains
+
+    ! This is the basic residual calculation routine...
+
+    subroutine residual (val, mask, res, si, order )
+      CCTK_REAL, dimension(:,:,:), intent(in) :: val
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+      CCTK_REAL, dimension(:,:,:), intent(out) :: res
+      CCTK_REAL, intent(in) :: si
+      CCTK_INT, intent(in) :: order
+      CCTK_REAL, parameter :: six = 6.0_wp
+      CCTK_REAL, parameter :: sixth = 1.0_wp / six
+      CCTK_REAL, parameter :: four = 4.0_wp
+      CCTK_REAL, parameter :: eighteen = 18.0_wp
+      CCTK_REAL, parameter :: eighteenth = 1.0_wp / eighteen
+      CCTK_REAL, parameter :: fifteen = 15.0_wp
+      CCTK_REAL, parameter :: sixty = 60.0_wp
+      CCTK_REAL, parameter :: sixtieth = 1.0_wp / sixty
+      integer :: i, j, k
+  
+      select case (order)
+      case (2) 
+
+        do k = 2, nz-1
+          do j = 2, ny-1
+            do i = 2, nx-1
+    
+              if (mask(i,j,k) > 0) then
+  
+                res(i,j,k) = val(i+1,j,k) + val(i-1,j,k) + val(i,j+1,k) + &
+                             val(i,j-1,k) + val(i,j,k+1) + val(i,j,k-1) - &
+                             six * val(i,j,k)
+                res(i,j,k) = sign(one,si) * res(i,j,k)
+  
+              end if
+    
+            end do
+          end do
+        end do
+
+      case (4)
+
+        do k = 3, nz-2
+          do j = 3, ny-2
+            do i = 3, nx-2
+
+              if (mask(i,j,k) > 0) then
+  
+                res(i,j,k) = - ( val(i+2,j,k) + val(i-2,j,k) + &
+                                 val(i,j+2,k) + val(i,j-2,k) + &
+                                 val(i,j,k+2) + val(i,j,k-2) ) + &
+                        four * ( val(i+1,j,k) + val(i-1,j,k) + &
+                                 val(i,j+1,k) + val(i,j-1,k) + &
+                                 val(i,j,k+1) + val(i,j,k-1) ) - &
+                      eighteen * val(i,j,k)
+                res(i,j,k) = sign(one,si) * res(i,j,k)
+
+              end if
+
+            end do
+          end do
+        end do
+
+      case (6)
+
+        do k = 4, nz-3
+          do j = 4, ny-3
+            do i = 4, nx-3
+
+              if (mask(i,j,k) > 0) then
+
+                res(i,j,k) =   ( val(i+3,j,k) + val(i-3,j,k) + &
+                                 val(i,j+3,k) + val(i,j-3,k) + &
+                                 val(i,j,k+3) + val(i,j,k-3) ) - &
+                         six * ( val(i+2,j,k) + val(i-2,j,k) + &
+                                 val(i,j+2,k) + val(i,j-2,k) + &
+                                 val(i,j,k+2) + val(i,j,k-2) ) + &
+                     fifteen * ( val(i+1,j,k) + val(i-1,j,k) + &
+                                 val(i,j+1,k) + val(i,j-1,k) + &
+                                 val(i,j,k+1) + val(i,j,k-1) ) - &
+                         sixty * val(i,j,k)
+                res(i,j,k) = sign(one,si) * res(i,j,k)
+
+              end if
+
+            end do
+          end do
+        end do
+
+      case default
+
+        call CCTK_WARN ( 0, 'Internal error: Order out of range' )
+  
+      end select
+
+    end subroutine residual
+
+
+    ! This is the routine calculating residuals for all variables, that
+    ! have not yet converged.
+
+    subroutine residual_all ( v1, v2, v3, v4, v5, v6, v7, v8, &
+                               v9, v10, v11, v12, v13, v14, v15, v16, &
+                               r1, r2, r3, r4, r5, r6, r7, r8, &
+                               r9, r10, r11, r12, r13, r14, r15, r16, &
+                               mask, si, order )
+
+      CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+                                                 v7, v8, v9, v10, v11, v12, &
+                                                 v13, v14, v15, v16
+      CCTK_REAL, dimension(:,:,:), intent(out) :: r1, r2, r3, r4, r5, r6, &
+                                                  r7, r8, r9, r10, r11, r12, &
+                                                  r13, r14, r15, r16
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+      CCTK_REAL, intent(in) :: si
+      CCTK_INT, intent(in) :: order
+
+      if ( cont(1) ) call residual ( v1, mask, r1, si, order )
+      if ( cont(2) ) call residual ( v2, mask, r2, si, order )
+      if ( cont(3) ) call residual ( v3, mask, r3, si, order )
+      if ( cont(4) ) call residual ( v4, mask, r4, si, order )
+      if ( cont(5) ) call residual ( v5, mask, r5, si, order )
+      if ( cont(6) ) call residual ( v6, mask, r6, si, order )
+      if ( cont(7) ) call residual ( v7, mask, r7, si, order )
+      if ( cont(8) ) call residual ( v8, mask, r8, si, order )
+      if ( cont(9) ) call residual ( v9, mask, r9, si, order )
+      if ( cont(10) ) call residual ( v10, mask, r10, si, order )
+      if ( cont(11) ) call residual ( v11, mask, r11, si, order )
+      if ( cont(12) ) call residual ( v12, mask, r12, si, order )
+      if ( cont(13) ) call residual ( v13, mask, r13, si, order )
+      if ( cont(14) ) call residual ( v14, mask, r14, si, order )
+      if ( cont(15) ) call residual ( v15, mask, r15, si, order )
+      if ( cont(16) ) call residual ( v16, mask, r16, si, order )
+
+    end subroutine residual_all
+
+
+    ! This routine multiplies two variables. Only does it according
+    ! to the mask and only for variables that have not converged yet.
+
+    subroutine multiply ( u1, u2, u3, u4, u5, u6, u7, u8, &
+                          u9, u10, u11, u12, u13, u14, u15, u16, &
+                          v1, v2, v3, v4, v5, v6, v7, v8, &
+                          v9, v10, v11, v12, v13, v14, v15, v16, &
+                          r1, r2, r3, r4, r5, r6, r7, r8, &
+                          r9, r10, r11, r12, r13, r14, r15, r16, &
+                          mask )
+
+      CCTK_REAL, dimension(:,:,:), intent(in) :: u1, u2, u3, u4, u5, u6, &
+                                                 u7, u8, u9, u10, u11, u12, &
+                                                 u13, u14, u15, u16
+      CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+                                                 v7, v8, v9, v10, v11, v12, &
+                                                 v13, v14, v15, v16
+      CCTK_REAL, dimension(:,:,:), intent(out) :: r1, r2, r3, r4, r5, r6, &
+                                                  r7, r8, r9, r10, r11, r12, &
+                                                  r13, r14, r15, r16
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+
+      if ( cont(1) ) where ( mask > 0 ) r1 = u1*v1
+      if ( cont(2) ) where ( mask > 0 ) r2 = u2*v2
+      if ( cont(3) ) where ( mask > 0 ) r3 = u3*v3
+      if ( cont(4) ) where ( mask > 0 ) r4 = u4*v4
+      if ( cont(5) ) where ( mask > 0 ) r5 = u5*v5
+      if ( cont(6) ) where ( mask > 0 ) r6 = u6*v6
+      if ( cont(7) ) where ( mask > 0 ) r7 = u7*v7
+      if ( cont(8) ) where ( mask > 0 ) r8 = u8*v8
+      if ( cont(9) ) where ( mask > 0 ) r9 = u9*v9
+      if ( cont(10) ) where ( mask > 0 ) r10 = u10*v10
+      if ( cont(11) ) where ( mask > 0 ) r11 = u11*v11
+      if ( cont(12) ) where ( mask > 0 ) r12 = u12*v12
+      if ( cont(13) ) where ( mask > 0 ) r13 = u13*v13
+      if ( cont(14) ) where ( mask > 0 ) r14 = u14*v14
+      if ( cont(15) ) where ( mask > 0 ) r15 = u15*v15
+      if ( cont(16) ) where ( mask > 0 ) r16 = u16*v16
+
+    end subroutine multiply
+                  
+
+    ! This routine takes the product of a scalar and a variable and
+    ! adds it to the product of another scalar and variable. Does it
+    ! according to the mask and only for variables that have not yet
+    ! converged.
+
+    subroutine multiply_sum ( u1, u2, u3, u4, u5, u6, u7, u8, &
+                              u9, u10, u11, u12, u13, u14, u15, u16, &
+                              v1, v2, v3, v4, v5, v6, v7, v8, &
+                              v9, v10, v11, v12, v13, v14, v15, v16, &
+                              c1, c2, mask )
+
+      CCTK_REAL, dimension(:,:,:), intent(inout) :: u1, u2, u3, u4, u5, u6, &
+                                                    u7, u8, u9, u10, u11, u12, &
+                                                    u13, u14, u15, u16
+      CCTK_REAL, dimension(:,:,:), intent(in) :: v1, v2, v3, v4, v5, v6, &
+                                                 v7, v8, v9, v10, v11, v12, &
+                                                 v13, v14, v15, v16
+      CCTK_REAL, dimension(16), intent(in) :: c1, c2
+      CCTK_INT, dimension(:,:,:), intent(in) :: mask
+
+      
+      if ( cont(1) ) where ( mask > 0 ) u1 = c1(1)*u1 + c2(1)*v1
+      if ( cont(2) ) where ( mask > 0 ) u2 = c1(2)*u2 + c2(2)*v2
+      if ( cont(3) ) where ( mask > 0 ) u3 = c1(3)*u3 + c2(3)*v3
+      if ( cont(4) ) where ( mask > 0 ) u4 = c1(4)*u4 + c2(4)*v4
+      if ( cont(5) ) where ( mask > 0 ) u5 = c1(5)*u5 + c2(5)*v5
+      if ( cont(6) ) where ( mask > 0 ) u6 = c1(6)*u6 + c2(6)*v6
+      if ( cont(7) ) where ( mask > 0 ) u7 = c1(7)*u7 + c2(7)*v7
+      if ( cont(8) ) where ( mask > 0 ) u8 = c1(8)*u8 + c2(8)*v8
+      if ( cont(9) ) where ( mask > 0 ) u9 = c1(9)*u9 + c2(9)*v9
+      if ( cont(10) ) where ( mask > 0 ) u10 = c1(10)*u10 + c2(10)*v10
+      if ( cont(11) ) where ( mask > 0 ) u11 = c1(11)*u11 + c2(11)*v11
+      if ( cont(12) ) where ( mask > 0 ) u12 = c1(12)*u12 + c2(12)*v12
+      if ( cont(13) ) where ( mask > 0 ) u13 = c1(13)*u13 + c2(13)*v13
+      if ( cont(14) ) where ( mask > 0 ) u14 = c1(14)*u14 + c2(14)*v14
+      if ( cont(15) ) where ( mask > 0 ) u15 = c1(15)*u15 + c2(15)*v15
+      if ( cont(16) ) where ( mask > 0 ) u16 = c1(16)*u16 + c2(16)*v16
+
+    end subroutine multiply_sum
+
+
+end module NoExcision_mod
diff --git a/src/cg.F90 b/src/cg.F90
new file mode 100644
index 0000000..71f56c0
--- /dev/null
+++ b/src/cg.F90
@@ -0,0 +1,634 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_CGInit_1 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: my_level, n_levels
+  CCTK_INT :: ierr
+  CCTK_INT, dimension(3) :: sym
+
+  nx = cctk_lsh(1); ny = cctk_lsh(2); nz = cctk_lsh(3)
+
+  loop_control = 0
+
+  call NoExcision_levelinfo ( cctkGH, my_level, n_levels )
+
+  if ( my_level == n_levels-1 ) then
+
+    where (abs(gxx) < 1d-8 .and. abs(alp) < 1d-8)
+      nes_mask = 1
+    elsewhere
+      nes_mask = 0
+    end where
+
+    call CCTK_INFO ( 'Starting smoothing procedure' )
+
+    resgxx = zero; resgxy = zero; resgxz = zero
+    resgyy = zero; resgyz = zero; resgzz = zero
+    reskxx = zero; reskxy = zero; reskxz = zero
+    reskyy = zero; reskyz = zero; reskzz = zero
+    res = zero; resx = zero; resy = zero; resz = zero
+
+    dgxx = zero; dgxy = zero; dgxz = zero
+    dgyy = zero; dgyz = zero; dgzz = zero
+    dkxx = zero; dkxy = zero; dkxz = zero
+    dkyy = zero; dkyz = zero; dkzz = zero
+    d = zero; dx = zero; dy = zero; dz = zero
+
+    qgxx = zero; qgxy = zero; qgxz = zero
+    qgyy = zero; qgyz = zero; qgzz = zero
+    qkxx = zero; qkxy = zero; qkxz = zero
+    qkyy = zero; qkyz = zero; qkzz = zero
+    q = zero; qx = zero; qy = zero; qz = zero
+
+    redgxx = zero; redgxy = zero; redgxz = zero
+    redgyy = zero; redgyz = zero; redgzz = zero
+    redkxx = zero; redkxy = zero; redkxz = zero
+    redkyy = zero; redkyz = zero; redkzz = zero
+    red = zero; redx = zero; redy = zero; redz = zero
+
+    ! r = b - A x.
+    ! Since x=0 and we actually use A':   b = -A' 0 and  r = b = -A' 0.
+
+    call residual_all ( gxx, gxy, gxz, gyy, gyz, gzz, &
+                        kxx, kxy, kxz, kyy, kyz, kzz, &
+                        alp, betax, betay, betaz, &
+                        resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                        reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                        res, resx, resy, resz, &
+                        nes_mask, -one, smoothing_order )
+
+    ! d = r = b.
+
+    call residual_all ( gxx, gxy, gxz, gyy, gyz, gzz, &
+                        kxx, kxy, kxz, kyy, kyz, kzz, &
+                        alp, betax, betay, betaz, &
+                        dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+                        dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+                        d, dx, dy, dz, &
+                        nes_mask, -one, smoothing_order )
+
+    ! red = r*r.
+
+    call multiply ( resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                    reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                    res, resx, resy, resz, &
+                    resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                    reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                    res, resx, resy, resz, &
+                    redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
+                    redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
+                    red, redx, redy, redz, nes_mask )
+                     
+    call CCTK_ReductionArrayHandle ( sum_handle, 'sum' )
+    if ( sum_handle .lt. 0 ) then
+      call CCTK_WARN(0,'Could not obtain a handle for sum reduction')
+    end if
+
+    call CCTK_ReductionArrayHandle ( infnorm_handle, 'norm_inf' )
+    if ( infnorm_handle .lt. 0 ) then
+      call CCTK_WARN(0,'Could not obtain a handle for norm_inf reduction')
+    end if
+
+    loop_counter = 0
+    loop_control = 1
+
+    sym_selector = 1
+
+  end if
+
+end subroutine NoExcision_CGInit_1
+
+
+subroutine NoExcision_CGInit_2 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT :: res_var, ierr
+  integer :: i
+  character(len=56) :: conv_message
+
+  if ( loop_control == 1 ) then
+    do i = 1, 16
+      if ( cont(i) ) then
+        call CCTK_VarIndex ( res_var, red_names(i) )
+        if ( res_var < 0 ) then
+          call CCTK_WARN ( 0, 'Could not get index to grid function red' )
+        end if
+
+        ! delta_new = r^T r.
+
+        call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
+                           CCTK_VARIABLE_REAL, delta_new(i), 1, res_var )
+
+        call CCTK_Reduce ( ierr, cctkGH, -1, infnorm_handle, 1, &
+                           CCTK_VARIABLE_REAL, infnormresid(i), 1, res_var )
+
+        ! Since we are reducing r*r, we have to take the square root.
+
+        infnormresid(i) = sqrt(infnormresid(i))
+
+      end if
+    end do
+
+    ! Check if some variables has already converged. This happens when the
+    ! variable is identically zero.
+
+    do i = 1, 16 
+      if ( cont(i) ) then
+        if ( infnormresid(i) < smoothing_eps ) then
+          write ( conv_message, '(a23,i8,a20,a5)' ) 'CG method converged in ', &
+                loop_counter, ' steps for variable ', var_names(i)
+          call CCTK_INFO ( conv_message )
+          cont(i) = .false.
+        end if 
+      end if
+    end do
+
+    if ( .not. any ( cont ) ) then
+      loop_control = 0
+    end if
+
+  end if
+
+end subroutine NoExcision_CGInit_2
+    
+
+subroutine NoExcision_CG_1 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  nx = cctk_lsh(1); ny=cctk_lsh(2); nz = cctk_lsh(3)
+
+  loop_counter = loop_counter + 1
+
+  ! Since d is zero outside of the active region we have:
+  ! q = A d = A' d.
+
+  call residual_all ( dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+                      dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+                      d, dx, dy, dz, &
+                      qgxx, qgxy, qgxz, qgyy, qgyz, qgzz, &
+                      qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, &
+                      q, qx, qy, qz, &
+                      nes_mask, one, smoothing_order )
+
+  ! red = d*q = d*(A d)
+
+  call multiply ( dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+                  dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+                  d, dx, dy, dz, &
+                  qgxx, qgxy, qgxz, qgyy, qgyz, qgzz, &
+                  qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, &
+                  q, qx, qy, qz, &
+                  redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
+                  redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
+                  red, redx, redy, redz, nes_mask )
+
+  sym_selector = 2
+
+end subroutine NoExcision_CG_1
+
+
+subroutine NoExcision_CG_2 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT res_var, ierr
+
+  integer :: i
+
+  if ( loop_control == 1 ) then
+    do i = 1, 16
+      if ( cont(i) ) then
+        call CCTK_VarIndex ( res_var, red_names(i) )
+        if ( res_var < 0 ) then
+          call CCTK_WARN ( 0, 'Could not get index to grid function red' )
+        end if
+
+        ! alpha = delta_new / ( d^T A d ).
+        call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
+                           CCTK_VARIABLE_REAL, alpha(i), 1, res_var )
+      
+        alpha(i) = delta_new(i) / alpha(i)
+
+      end if
+    end do
+
+  end if
+
+end subroutine NoExcision_CG_2
+
+
+subroutine NoExcision_CG_3 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  nx = cctk_lsh(1); ny=cctk_lsh(2); nz = cctk_lsh(3)
+
+  ! x = x + alpha * d.
+
+  call multiply_sum ( gxx, gxy, gxz, gyy, gyz, gzz, &
+                      kxx, kxy, kxz, kyy, kyz, kzz, &
+                      alp, betax, betay, betaz, &
+                      dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+                      dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+                      d, dx, dy, dz, &
+                      (/ (one, i=1, 16) /), alpha, nes_mask )
+                      
+  if ( mod ( loop_counter, 50 ) == 0 ) then
+
+    ! Restart:
+    ! r = b - A x = - A' x
+
+    call residual_all ( gxx, gxy, gxz, gyy, gyz, gzz, &
+                        kxx, kxy, kxz, kyy, kyz, kzz, &
+                        alp, betax, betay, betaz, &
+                        resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                        reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                        res, resx, resy, resz, &
+                        nes_mask, -one, smoothing_order )
+  else
+
+    ! r = r - alpha q.
+
+    call multiply_sum ( resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                        reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                        res, resx, resy, resz, &
+                        qgxx, qgxy, qgxz, qgyy, qgyz, qgzz, &
+                        qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, &
+                        q, qx, qy, qz, &
+                        (/ (one, i=1, 16) /), -alpha, nes_mask )
+                      
+  end if
+
+  ! delta_old = delta_new.
+
+  where ( cont ) delta_old = delta_new
+
+  ! red = r*r.
+  call multiply ( resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                  reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                  res, resx, resy, resz, &
+                  resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                  reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                  res, resx, resy, resz, &
+                  redgxx, redgxy, redgxz, redgyy, redgyz, redgzz, &
+                  redkxx, redkxy, redkxz, redkyy, redkyz, redkzz, &
+                  red, redx, redy, redz, nes_mask )
+
+  sym_selector = 3
+
+end subroutine NoExcision_CG_3
+
+
+subroutine NoExcision_CG_4 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT res_var, ierr
+
+  integer :: i
+
+  if ( loop_control == 1 ) then
+    do i = 1, 16
+      if ( cont(i) ) then
+        call CCTK_VarIndex ( res_var, red_names(i) )
+        if ( res_var < 0 ) then
+          call CCTK_WARN ( 0, 'Could not get index to grid function red' )
+        end if
+
+        ! delta_new = r^T r.
+        call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, &
+                           CCTK_VARIABLE_REAL, delta_new(i), 1, res_var )
+
+        call CCTK_Reduce ( ierr, cctkGH, -1, infnorm_handle, 1, &
+                           CCTK_VARIABLE_REAL, infnormresid(i), 1, res_var )
+        infnormresid(i) = sqrt(infnormresid(i))
+      end if
+    end do
+
+  end if
+
+end subroutine NoExcision_CG_4
+
+
+subroutine NoExcision_CG_5 (CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  integer :: i
+  character(len=192) :: res_message, var_message
+  character(len=56) :: conv_message
+  character(len=18) :: iter_message
+
+  nx = cctk_lsh(1); ny=cctk_lsh(2); nz = cctk_lsh(3)
+
+  ! beta = delta_new / delta_old
+
+  where ( cont ) beta = delta_new / delta_old
+
+  ! d = r + beta d
+
+  call multiply_sum ( dgxx, dgxy, dgxz, dgyy, dgyz, dgzz, &
+                      dkxx, dkxy, dkxz, dkyy, dkyz, dkzz, &
+                      d, dx, dy, dz, &
+                      resgxx, resgxy, resgxz, resgyy, resgyz, resgzz, &
+                      reskxx, reskxy, reskxz, reskyy, reskyz, reskzz, &
+                      res, resx, resy, resz, &
+                      beta, (/ (one, i=1, 16) /), nes_mask )
+  
+  if ( verbose > 0 .and. mod(loop_counter,100) == 0 ) then
+    write (var_message, '(16(a3,a9))' ) (' | ', var_names(i), i=1, 16)
+    call CCTK_INFO ( var_message )
+  end if
+  if ( verbose > 0 .and. mod(loop_counter,10) == 0 ) then
+    write (iter_message, '(a10,i8)' ), 'Iteration ', loop_counter
+    write (res_message, '(16(a3,es9.3))' ) (' | ', infnormresid(i), i=1, 16 )
+    call CCTK_INFO ( iter_message )
+    call CCTK_INFO ( res_message )
+  end if
+  
+  ! Check if any variables have converged during this iteration.
+
+  do i = 1, 16
+    if ( cont(i) ) then
+      if ( infnormresid(i) < smoothing_eps ) then
+        write ( conv_message, '(a23,i8,a20,a5)' ) 'CG method converged in ', &
+                loop_counter, ' steps for variable ', var_names(i)
+        call CCTK_INFO ( conv_message )
+        cont(i) = .false.
+      end if
+    end if
+  end do
+
+  ! If all variables have converged we exit.
+
+  if ( .not. any ( cont ) ) then
+    loop_control = 0
+  end if
+
+  sym_selector = 4
+
+end subroutine NoExcision_CG_5
+
+
+! I leave this in here in case somebody uses CartGrid3D symmetries...
+! Not tested, though...
+
+subroutine NoExcision_SetSym(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+
+  CCTK_INT :: ierr
+  CCTK_INT, dimension(3) :: sym
+
+  sym = 1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgxx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgyy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgzz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskxx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskyy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskzz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgxx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgyy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgzz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkxx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkyy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkzz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgxx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgyy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgzz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkxx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkyy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkzz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::res' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::d' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::q' )
+
+  call SetCartSymGN ( ierr, cctkGH, sym, 'noexcision::cg_red_all' )
+
+  sym(1) = -1; sym(2) = -1; sym(3) = 1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgxy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskxy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgxy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkxy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgxy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkxy' )
+
+  sym(1) = -1; sym(2) = 1; sym(3) = -1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgxz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskxz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgxz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkxz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgxz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkxz' )
+
+  sym(1) = 1; sym(2) = -1; sym(3) = -1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resgyz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::reskyz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dgyz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dkyz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qgyz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qkyz' )
+
+  sym(1) = -1; sym(2) = 1; sym(3) = 1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dx' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qx' )
+
+  sym(1) = 1; sym(2) = -1; sym(3) = 1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dy' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qy' )
+
+  sym(1) = 1; sym(2) = 1; sym(3) = -1
+
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::resz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::dz' )
+  call SetCartSymVN ( ierr, cctkGH, sym, 'noexcision::qz' )
+
+end subroutine NoExcision_SetSym
+
+
+subroutine NoExcision_CGApplySym(CCTK_ARGUMENTS)
+
+  use NoExcision_mod
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: ierr
+
+  if ( loop_control > 0 ) then
+
+    if ( sym_selector <1 .or. sym_selector > 4 ) then
+      call CCTK_WARN ( 0, 'Internal error. Inconsistent symmetry selector' )
+    end if
+  
+    select case ( sym_selector )
+    case (1, 4)
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_d_lapse', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_d_lapse for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_d_shift', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_d_shift for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_d_curv', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_d_curv for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_d_metric', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_d_metric for boundary condition' )
+      end if
+    end select
+  
+    select case ( sym_selector )
+    case (1, 3)
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_res_lapse', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_res_lapse for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_res_shift', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_res_shift for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_res_curv', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_res_curv for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_res_metric', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_res_metric for boundary condition' )
+      end if
+    end select
+  
+    select case ( sym_selector )
+    case (1, 2, 3)
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_red_all', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_red_all for boundary condition' )
+      end if
+    end select
+  
+    select case ( sym_selector )
+    case (2)
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_q_lapse', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_q_lapse for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_q_shift', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_q_shift for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_q_curv', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_q_curv for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                             'NoExcision::cg_q_metric', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select cg_q_metric for boundary condition' )
+      end if
+    end select
+  
+    select case ( sym_selector )
+    case (3)
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                                'ADMBase::lapse', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select lapse for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                                'ADMBase::shift', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select shift for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                                'ADMBase::curv', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select curv for boundary condition' )
+      end if
+  
+      ierr = Boundary_SelectGroupForBC ( cctkGH, CCTK_ALL_FACES, 1, -1, &
+                                                'ADMBase::metric', 'None' )
+      if ( ierr /= 0 ) then
+        call CCTK_WARN ( 0, 'Could not select metric for boundary condition' )
+      end if
+    end select
+
+  end if
+
+end subroutine NoExcision_CGApplySym
diff --git a/src/getlevelinfo.cc b/src/getlevelinfo.cc
new file mode 100644
index 0000000..4b68abb
--- /dev/null
+++ b/src/getlevelinfo.cc
@@ -0,0 +1,45 @@
+/* $Header$ */
+
+#include <cstdio>
+#include <cmath>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_DefineThorn.h"
+
+#include "carpet.hh"
+
+namespace NoExcision {
+
+#ifdef HAVE_CARPET
+  using namespace Carpet;
+#endif
+
+  extern "C"
+    void NoExcision_levelinfo(const cGH * cctkGH,
+                              CCTK_INT * my_level,
+                              CCTK_INT * n_levels)
+    {
+
+      *my_level = 0;
+      *n_levels = 1;
+
+#ifdef HAVE_CARPET
+
+      if ( CCTK_IsThornActive ( "Carpet" ) ) {
+        *my_level = reflevel;
+        *n_levels = reflevels;
+      } 
+     
+#endif
+
+    }
+
+  extern "C"
+    CCTK_FCALL void CCTK_FNAME (NoExcision_levelinfo) (CCTK_POINTER_TO_CONST * cctkGH, CCTK_INT * my_level, CCTK_INT * n_levels)
+    {
+      NoExcision_levelinfo((const cGH *)* cctkGH,  my_level, n_levels);
+    }
+
+} // namespace NoExcision
diff --git a/src/make.code.defn b/src/make.code.defn
index 7dd625c..a174c46 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS = overwrite.F90
+SRCS = NoExcision_mod.F90 overwrite.F90 smooth.F90 reduce.F90 reduce.c getlevelinfo.cc cg.F90
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100644
index 0000000..e7f8124
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1 @@
+cg.F90.o:		NoExcision_mod.F90.o
diff --git a/src/overwrite.F90 b/src/overwrite.F90
index 6c6089f..3df9ffb 100644
--- a/src/overwrite.F90
+++ b/src/overwrite.F90
@@ -10,85 +10,114 @@ subroutine NoExcision_Overwrite (CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
   DECLARE_CCTK_PARAMETERS
-  
+
   CCTK_REAL, parameter :: zero=0
   CCTK_REAL :: dist2 (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
-  CCTK_REAL :: cx, cy, cz, rad, width
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz, width
   integer   :: n
+
+  integer, parameter :: sm_linear = 1
+  integer, parameter :: sm_spline = 2
+  integer, parameter :: sm_cosine = 3
+  integer :: sm_type
+
+  if (CCTK_EQUALS(method,"old")) then
+    do n = 1, num_regions
+
+       cx = centre_x(n)
+       cy = centre_y(n)
+       cz = centre_z(n)
+       if (CCTK_EQUALS(region_shape(n), "sphere")) then
+          radx = radius(n)
+          rady = radius(n)
+          radz = radius(n)
+       else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+          radx = radius_x(n)
+          rady = radius_y(n)
+          radz = radius_z(n)
+       else
+          call CCTK_WARN (0, "internal error")
+       end if
+
+       if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+          sm_type = sm_linear
+       else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+          sm_type = sm_spline
+       else if (CCTK_EQUALS (smoothing_function(n), "cosine")) then
+          sm_type = sm_cosine
+       else
+          call CCTK_WARN (0, "internal error")
+       end if
+       width = smoothing_zone_width(n)
+
+       dist2 =   ((x - cx) / radx)**2 + ((y - cy) / rady)**2 &
+            &  + ((z - cz) / radz)**2
+
+       if (overwrite_geometry(n) /= 0) then
+
+          if (conformal_state >= 1) then
+             where (dist2 <= 1)
+                psi = 1
+             end where
+          end if
+          if (conformal_state >= 2) then
+             where (dist2 <= 1)
+                psix = 0
+                psiy = 0
+                psiz = 0
+             end where
+          end if
+          if (conformal_state >= 3) then
+             where (dist2 <= 1)
+                psixx = 0
+                psixy = 0
+                psixz = 0
+                psiyy = 0
+                psiyz = 0
+                psizz = 0
+             end where
+          end if
   
-  do n = 1, num_regions
-     
-     cx = centre_x(n)
-     cy = centre_y(n)
-     cz = centre_z(n)
-     rad = radius(n)
-     width = smoothing_zone_width(n)
-     
-     dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2
-     
-     if (overwrite_geometry(n) /= 0) then
-        
-        if (conformal_state >= 1) then
-           where (dist2 <= rad**2)
-              psi = 1
-           end where
-        end if
-        if (conformal_state >= 2) then
-           where (dist2 <= rad**2)
-              psix = 0
-              psiy = 0
-              psiz = 0
-           end where
-        end if
-        if (conformal_state >= 3) then
-           where (dist2 <= rad**2)
-              psixx = 0
-              psixy = 0
-              psixz = 0
-              psiyy = 0
-              psiyz = 0
-              psizz = 0
-           end where
-        end if
-        
-        where (dist2 <= rad**2)
-           gxx = smooth (gxx,  Minkowski_scale(n), dist2)
-           gxy = smooth (gxy,  zero              , dist2)
-           gxz = smooth (gxz,  zero              , dist2)
-           gyy = smooth (gyy,  Minkowski_scale(n), dist2)
-           gyz = smooth (gyz,  zero              , dist2)
-           gzz = smooth (gzz,  Minkowski_scale(n), dist2)
-           kxx = smooth (kxx,  zero              , dist2)
-           kxy = smooth (kxy,  zero              , dist2)
-           kxz = smooth (kxz,  zero              , dist2)
-           kyy = smooth (kyy,  zero              , dist2)
-           kyz = smooth (kyz,  zero              , dist2)
-           kzz = smooth (kzz,  zero              , dist2)
-        end where
-        
-     end if
-     
-     if (overwrite_lapse(n) /= 0) then
-        
-        where (dist2 <= rad**2)
-           alp = smooth (alp, lapse_scale(n), dist2)
-        end where
-        
-     end if
-        
-     if (overwrite_shift(n) /= 0) then
-        
-        if (shift_state /= 0) then
-           where (dist2 <= rad**2)
-              betax = smooth (betax, zero, dist2)
-              betay = smooth (betay, zero, dist2)
-              betaz = smooth (betaz, zero, dist2)
-           end where
-        end if
-        
-     end if
-     
-  end do
+          where (dist2 <= 1)
+             gxx = smooth (gxx,  Minkowski_scale(n), dist2)
+             gxy = smooth (gxy,  zero              , dist2)
+             gxz = smooth (gxz,  zero              , dist2)
+             gyy = smooth (gyy,  Minkowski_scale(n), dist2)
+             gyz = smooth (gyz,  zero              , dist2)
+             gzz = smooth (gzz,  Minkowski_scale(n), dist2)
+             kxx = smooth (kxx,  zero              , dist2)
+             kxy = smooth (kxy,  zero              , dist2)
+             kxz = smooth (kxz,  zero              , dist2)
+             kyy = smooth (kyy,  zero              , dist2)
+             kyz = smooth (kyz,  zero              , dist2)
+             kzz = smooth (kzz,  zero              , dist2)
+          end where
+
+       end if
+
+       if (overwrite_lapse(n) /= 0) then
+
+          where (dist2 <= 1)
+             alp = smooth (alp, lapse_scale(n), dist2)
+          end where
+
+       end if
+
+       if (overwrite_shift(n) /= 0) then
+
+          if (shift_state /= 0) then
+             where (dist2 <= 1)
+                betax = smooth (betax, zero, dist2)
+                betay = smooth (betay, zero, dist2)
+                betaz = smooth (betaz, zero, dist2)
+             end where
+          end if
+
+       end if
+
+    end do
+    
+  end if
   
 contains
   
@@ -98,13 +127,32 @@ contains
     CCTK_REAL, intent(in) :: goodval
     CCTK_REAL, intent(in) :: dist2
     
-    if (dist2 <= (rad-width)**2) then
+    CCTK_REAL, parameter :: zero=0, two=2, half=1/two
+    integer,   parameter :: rk = kind(zero)
+    CCTK_REAL, parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068_rk
+    CCTK_REAL :: dist
+    CCTK_REAL :: x, f
+    
+    if (dist2 <= (1 - width)**2) then
        newval = goodval
-    else if (width == 0 .or. dist2 >= rad**2) then
+    else if (width == 0 .or. dist2 >= 1) then
        newval = oldval
     else
-       newval =   goodval &
-            &   + (oldval - goodval) * (sqrt(dist2) - (rad-width)) / width
+       dist = sqrt(dist2)
+       x = 1 - (1 - dist) / width
+       select case (sm_type)
+       case (sm_linear)
+          f = x
+       case (sm_spline)
+          if (x < half) then
+             f = 2 * x**2
+          else
+             f = 1 - 2 * (1-x)**2
+          end if
+       case (sm_cosine)
+          f = (1 - cos (pi * x)) / 2
+       end select
+       newval = oldval * f + goodval * (1 - f)
     end if
     
   end function smooth
diff --git a/src/overwriteBSSN.F90 b/src/overwriteBSSN.F90
new file mode 100644
index 0000000..cef380a
--- /dev/null
+++ b/src/overwriteBSSN.F90
@@ -0,0 +1,248 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_OverwriteBSSN (CCTK_ARGUMENTS)
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  
+  CCTK_REAL, parameter :: zero=0
+  CCTK_REAL :: dist2 (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz, width, scale
+  integer   :: timelevels
+  integer   :: n
+  
+  integer, parameter :: sm_linear = 1
+  integer, parameter :: sm_spline = 2
+  integer :: sm_type
+  
+  do n = 1, num_regions
+     
+     if (cctk_iteration == iteration(n)) then
+        
+        cx = centre_x(n)
+        cy = centre_y(n)
+        cz = centre_z(n)
+        if (CCTK_EQUALS(region_shape(n), "sphere")) then
+           radx = radius(n)
+           rady = radius(n)
+           radz = radius(n)
+        else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+           radx = radius_x(n)
+           rady = radius_y(n)
+           radz = radius_z(n)
+        else
+           call CCTK_WARN (0, "internal error")
+        end if
+
+        if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+           sm_type = sm_linear
+        else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+           sm_type = sm_spline
+        else
+           call CCTK_WARN (0, "internal error")
+        end if
+        width = smoothing_zone_width(n)
+        
+        dist2 =   ((x - cx) / radx)**2 + ((y - cy) / rady)**2 &
+             &  + ((z - cz) / radz)**2
+
+        if (overwrite_geometry(n) /= 0) then
+           
+           scale = Minkowski_scale(n)
+           
+           if (conformal_state >= 1) then
+              where (dist2 <= 1)
+                 psi = 1
+              end where
+           end if
+           if (conformal_state >= 2) then
+              where (dist2 <= 1)
+                 psix = 0
+                 psiy = 0
+                 psiz = 0
+              end where
+           end if
+           if (conformal_state >= 3) then
+              where (dist2 <= 1)
+                 psixx = 0
+                 psixy = 0
+                 psixz = 0
+                 psiyy = 0
+                 psiyz = 0
+                 psizz = 0
+              end where
+           end if
+
+           where (dist2 <= 1)
+              ADM_BS_phi = smooth (ADM_BS_phi, zero , dist2)
+              ADM_BS_gxx = smooth (ADM_BS_gxx, scale, dist2)
+              ADM_BS_gxy = smooth (ADM_BS_gxy, zero , dist2)
+              ADM_BS_gxz = smooth (ADM_BS_gxz, zero , dist2)
+              ADM_BS_gyy = smooth (ADM_BS_gyy, scale, dist2)
+              ADM_BS_gyz = smooth (ADM_BS_gyz, zero , dist2)
+              ADM_BS_gzz = smooth (ADM_BS_gzz, scale, dist2)
+              ADM_BS_Gx  = smooth (ADM_BS_Gx , zero , dist2)
+              ADM_BS_Gy  = smooth (ADM_BS_Gy , zero , dist2)
+              ADM_BS_Gz  = smooth (ADM_BS_Gz , zero , dist2)
+              ADM_BS_K   = smooth (ADM_BS_K  , zero , dist2)
+              ADM_BS_Axx = smooth (ADM_BS_Axx, zero , dist2)
+              ADM_BS_Axy = smooth (ADM_BS_Axy, zero , dist2)
+              ADM_BS_Axz = smooth (ADM_BS_Axz, zero , dist2)
+              ADM_BS_Ayy = smooth (ADM_BS_Ayy, zero , dist2)
+              ADM_BS_Ayz = smooth (ADM_BS_Ayz, zero , dist2)
+              ADM_BS_Azz = smooth (ADM_BS_Azz, zero , dist2)
+           end where
+           call CCTK_ActiveTimeLevelsVN &
+                (timelevels, cctkGH, "ADM_BSSN::ADM_BS_phi")
+           if (timelevels > 1) then
+              where (dist2 <= 1)
+                 ADM_BS_phi_p = smooth (ADM_BS_phi_p, zero , dist2)
+                 ADM_BS_gxx_p = smooth (ADM_BS_gxx_p, scale, dist2)
+                 ADM_BS_gxy_p = smooth (ADM_BS_gxy_p, zero , dist2)
+                 ADM_BS_gxz_p = smooth (ADM_BS_gxz_p, zero , dist2)
+                 ADM_BS_gyy_p = smooth (ADM_BS_gyy_p, scale, dist2)
+                 ADM_BS_gyz_p = smooth (ADM_BS_gyz_p, zero , dist2)
+                 ADM_BS_gzz_p = smooth (ADM_BS_gzz_p, scale, dist2)
+                 ADM_BS_Gx_p  = smooth (ADM_BS_Gx_p , zero , dist2)
+                 ADM_BS_Gy_p  = smooth (ADM_BS_Gy_p , zero , dist2)
+                 ADM_BS_Gz_p  = smooth (ADM_BS_Gz_p , zero , dist2)
+                 ADM_BS_K_p   = smooth (ADM_BS_K_p  , zero , dist2)
+                 ADM_BS_Axx_p = smooth (ADM_BS_Axx_p, zero , dist2)
+                 ADM_BS_Axy_p = smooth (ADM_BS_Axy_p, zero , dist2)
+                 ADM_BS_Axz_p = smooth (ADM_BS_Axz_p, zero , dist2)
+                 ADM_BS_Ayy_p = smooth (ADM_BS_Ayy_p, zero , dist2)
+                 ADM_BS_Ayz_p = smooth (ADM_BS_Ayz_p, zero , dist2)
+                 ADM_BS_Azz_p = smooth (ADM_BS_Azz_p, zero , dist2)
+              end where
+           end if
+           if (timelevels > 2) then
+              where (dist2 <= 1)
+                 ADM_BS_phi_p_p = smooth (ADM_BS_phi_p_p, zero , dist2)
+                 ADM_BS_gxx_p_p = smooth (ADM_BS_gxx_p_p, scale, dist2)
+                 ADM_BS_gxy_p_p = smooth (ADM_BS_gxy_p_p, zero , dist2)
+                 ADM_BS_gxz_p_p = smooth (ADM_BS_gxz_p_p, zero , dist2)
+                 ADM_BS_gyy_p_p = smooth (ADM_BS_gyy_p_p, scale, dist2)
+                 ADM_BS_gyz_p_p = smooth (ADM_BS_gyz_p_p, zero , dist2)
+                 ADM_BS_gzz_p_p = smooth (ADM_BS_gzz_p_p, scale, dist2)
+                 ADM_BS_Gx_p_p  = smooth (ADM_BS_Gx_p_p , zero , dist2)
+                 ADM_BS_Gy_p_p  = smooth (ADM_BS_Gy_p_p , zero , dist2)
+                 ADM_BS_Gz_p_p  = smooth (ADM_BS_Gz_p_p , zero , dist2)
+                 ADM_BS_K_p_p   = smooth (ADM_BS_K_p_p  , zero , dist2)
+                 ADM_BS_Axx_p_p = smooth (ADM_BS_Axx_p_p, zero , dist2)
+                 ADM_BS_Axy_p_p = smooth (ADM_BS_Axy_p_p, zero , dist2)
+                 ADM_BS_Axz_p_p = smooth (ADM_BS_Axz_p_p, zero , dist2)
+                 ADM_BS_Ayy_p_p = smooth (ADM_BS_Ayy_p_p, zero , dist2)
+                 ADM_BS_Ayz_p_p = smooth (ADM_BS_Ayz_p_p, zero , dist2)
+                 ADM_BS_Azz_p_p = smooth (ADM_BS_Azz_p_p, zero , dist2)
+              end where
+           end if
+
+        end if
+
+        if (overwrite_lapse(n) /= 0) then
+           
+           scale = lapse_scale(n)
+           
+           where (dist2 <= 1)
+              alp          = smooth (alp         , scale, dist2)
+              ADM_BS_dtalp = smooth (ADM_BS_dtalp, zero , dist2)
+           end where
+           call CCTK_ActiveTimeLevelsVN &
+                (timelevels, cctkGH, "ADMBase::alp")
+           if (timelevels > 1) then
+              where (dist2 <= 1)
+                 alp_p          = smooth (alp_p         , scale, dist2)
+                 ADM_BS_dtalp_p = smooth (ADM_BS_dtalp_p, zero , dist2)
+              end where
+           end if
+           if (timelevels > 2) then
+              where (dist2 <= 1)
+                 alp_p_p          = smooth (alp_p_p         , scale, dist2)
+                 ADM_BS_dtalp_p_p = smooth (ADM_BS_dtalp_p_p, zero , dist2)
+              end where
+           end if
+           
+        end if
+
+        if (overwrite_shift(n) /= 0) then
+
+           if (shift_state /= 0) then
+              where (dist2 <= 1)
+                 betax     = smooth (betax    , zero, dist2)
+                 betay     = smooth (betay    , zero, dist2)
+                 betaz     = smooth (betaz    , zero, dist2)
+                 ADM_BS_Bx = smooth (ADM_BS_Bx, zero, dist2)
+                 ADM_BS_By = smooth (ADM_BS_By, zero, dist2)
+                 ADM_BS_Bz = smooth (ADM_BS_Bz, zero, dist2)
+              end where
+              call CCTK_ActiveTimeLevelsVN &
+                   (timelevels, cctkGH, "ADMBase::betax")
+              if (timelevels > 1) then
+                 where (dist2 <= 1)
+                    betax_p     = smooth (betax_p    , zero, dist2)
+                    betay_p     = smooth (betay_p    , zero, dist2)
+                    betaz_p     = smooth (betaz_p    , zero, dist2)
+                    ADM_BS_Bx_p = smooth (ADM_BS_Bx_p, zero, dist2)
+                    ADM_BS_By_p = smooth (ADM_BS_By_p, zero, dist2)
+                    ADM_BS_Bz_p = smooth (ADM_BS_Bz_p, zero, dist2)
+                 end where
+              end if
+              if (timelevels > 2) then
+                 where (dist2 <= 1)
+                    betax_p_p     = smooth (betax_p_p    , zero, dist2)
+                    betay_p_p     = smooth (betay_p_p    , zero, dist2)
+                    betaz_p_p     = smooth (betaz_p_p    , zero, dist2)
+                    ADM_BS_Bx_p_p = smooth (ADM_BS_Bx_p_p, zero, dist2)
+                    ADM_BS_By_p_p = smooth (ADM_BS_By_p_p, zero, dist2)
+                    ADM_BS_Bz_p_p = smooth (ADM_BS_Bz_p_p, zero, dist2)
+                 end where
+              end if
+           end if
+
+        end if
+
+     end if
+
+  end do
+  
+contains
+  
+  elemental function smooth (oldval, goodval, dist2) result (newval)
+    CCTK_REAL             :: newval
+    CCTK_REAL, intent(in) :: oldval
+    CCTK_REAL, intent(in) :: goodval
+    CCTK_REAL, intent(in) :: dist2
+    
+    CCTK_REAL, parameter :: two=2, half=1/two
+    CCTK_REAL :: dist
+    CCTK_REAL :: x, f
+    
+    if (dist2 <= (1 - width)**2) then
+       newval = goodval
+    else if (width == 0 .or. dist2 >= 1) then
+       newval = oldval
+    else
+       dist = sqrt(dist2)
+       x = 1 - (1 - dist) / width
+       select case (sm_type)
+       case (sm_linear)
+          f = x
+       case (sm_spline)
+          if (x < half) then
+             f = 2 * x**2
+          else
+             f = 1 - 2 * (1-x)**2
+          end if
+       end select
+       newval = oldval * f + goodval * (1 - f)
+    end if
+    
+  end function smooth
+  
+end subroutine NoExcision_OverwriteBSSN
diff --git a/src/reduce.F90 b/src/reduce.F90
new file mode 100644
index 0000000..ab7b70e
--- /dev/null
+++ b/src/reduce.F90
@@ -0,0 +1,97 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_Reduce (cctk_iteration, cctk_lsh, rhs, x, y, z)
+  implicit none
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  
+  integer :: cctk_iteration
+  integer :: cctk_lsh(3)
+  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: rhs, x, y, z
+  
+  CCTK_REAL, parameter :: zero=0
+  CCTK_REAL, dimension(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) :: dist2
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz, width
+  integer   :: n
+  
+  integer, parameter :: sm_linear = 1
+  integer, parameter :: sm_spline = 2
+  integer :: sm_type
+  
+  do n = 1, num_regions
+     if (reduce_rhs(n) /= 0) then
+
+        cx = centre_x(n)
+        cy = centre_y(n)
+        cz = centre_z(n)
+        if (CCTK_EQUALS(region_shape(n), "sphere")) then
+           radx = radius(n)
+           rady = radius(n)
+           radz = radius(n)
+        else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+           radx = radius_x(n)
+           rady = radius_y(n)
+           radz = radius_z(n)
+        else
+           call CCTK_WARN (0, "internal error")
+        end if
+
+        if (CCTK_EQUALS (smoothing_function(n), "linear")) then
+           sm_type = sm_linear
+        else if (CCTK_EQUALS (smoothing_function(n), "spline")) then
+           sm_type = sm_spline
+        else
+           call CCTK_WARN (0, "internal error")
+        end if
+        width = smoothing_zone_width(n)
+
+        dist2 = (x - cx)**2 + (y - cy)**2 + (z - cz)**2
+
+#warning "TODO: reduce not only the RHS, but also drive the variables towards Minkowski"
+        where (dist2 <= 1)
+           rhs = smooth (rhs, rhs * reduction_factor(n), dist2)
+        end where
+        
+     end if
+  end do
+  
+contains
+  
+  elemental function smooth (oldval, goodval, dist2) result (newval)
+    CCTK_REAL             :: newval
+    CCTK_REAL, intent(in) :: oldval
+    CCTK_REAL, intent(in) :: goodval
+    CCTK_REAL, intent(in) :: dist2
+    
+    CCTK_REAL, parameter :: two=2, half=1/two
+    CCTK_REAL :: dist
+    CCTK_REAL :: x, f
+    
+    if (dist2 <= (1 - width)**2) then
+       newval = goodval
+    else if (width == 0 .or. dist2 >= 1) then
+       newval = oldval
+    else
+       dist = sqrt(dist2)
+       x = 1 - (1 - dist) / width
+       select case (sm_type)
+       case (sm_linear)
+          f = x
+       case (sm_spline)
+          if (x < half) then
+             f = 2 * x**2
+          else
+             f = 1 - 2 * (1-x)**2
+          end if
+       end select
+       newval = oldval * f + goodval * (1 - f)
+    end if
+    
+  end function smooth
+  
+end subroutine NoExcision_Reduce
diff --git a/src/reduce.c b/src/reduce.c
new file mode 100644
index 0000000..0b44ca2
--- /dev/null
+++ b/src/reduce.c
@@ -0,0 +1,34 @@
+/* $Header$ */
+
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void CCTK_FCALL
+CCTK_FNAME(NoExcision_Reduce) (int const * cctk_iteration,
+                               int const * cctk_lsh,
+                               CCTK_REAL * rhs,
+                               CCTK_REAL const * x,
+                               CCTK_REAL const * y,
+                               CCTK_REAL const * z);
+
+void
+NoExcision_Reduce (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  int var, rhs;
+  void * ptr;
+  
+  for (var=0; var<CCTK_NumVars(); ++var) {
+    rhs = MoLQueryEvolvedRHS (var);
+    if (rhs >= 0) {
+      ptr = CCTK_VarDataPtrI (cctkGH, 0, rhs);
+      assert (ptr);
+      CCTK_FNAME(NoExcision_Reduce) (& cctk_iteration, cctk_lsh, ptr, x, y, z);
+    }
+  }
+}
diff --git a/src/smooth.F90 b/src/smooth.F90
new file mode 100644
index 0000000..d6df156
--- /dev/null
+++ b/src/smooth.F90
@@ -0,0 +1,150 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine NoExcision_Smooth (CCTK_ARGUMENTS)
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_REAL, parameter :: zero=0
+  CCTK_REAL :: mask (cctk_lsh(1), cctk_lsh(2), cctk_lsh(3))
+  CCTK_REAL :: cx, cy, cz, radx, rady, radz
+  integer   :: n
+  integer   :: iter
+  integer   :: ierr
+  
+  do iter = 1, smoothing_iterations
+     
+     do n = 1, num_regions
+        
+        cx = centre_x(n)
+        cy = centre_y(n)
+        cz = centre_z(n)
+        if (CCTK_EQUALS(region_shape(n), "sphere")) then
+           radx = radius(n)
+           rady = radius(n)
+           radz = radius(n)
+        else if (CCTK_EQUALS(region_shape(n), "ellipsoid")) then
+           radx = radius_x(n)
+           rady = radius_y(n)
+           radz = radius_z(n)
+        else
+           call CCTK_WARN (0, "internal error")
+        end if
+        
+        mask = 1 - sqrt (((x - cx) / radx)**2 + &
+             &           ((y - cy) / rady)**2 + &
+             &           ((z - cz) / radz)**2)
+        where (mask <= 0)
+           mask = 0             ! outside
+        else where (mask >= smoothing_zone_width(n))
+           mask = 1             ! far inside
+        elsewhere
+           mask = mask / smoothing_zone_width(n) ! a bit inside
+        end where
+        
+        if (overwrite_geometry(n) /= 0) then
+           
+           if (conformal_state >= 1) then
+              call smooth (psi)
+           end if
+           if (conformal_state >= 2) then
+              call smooth (psix)
+              call smooth (psiy)
+              call smooth (psiz)
+           end if
+           if (conformal_state >= 3) then
+              call smooth (psixx)
+              call smooth (psixy)
+              call smooth (psixz)
+              call smooth (psiyy)
+              call smooth (psiyz)
+              call smooth (psizz)
+           end if
+           
+           call smooth (gxx)
+           call smooth (gxy)
+           call smooth (gxz)
+           call smooth (gyy)
+           call smooth (gyz)
+           call smooth (gzz)
+           call smooth (kxx)
+           call smooth (kxy)
+           call smooth (kxz)
+           call smooth (kyy)
+           call smooth (kyz)
+           call smooth (kzz)
+           
+        end if
+        
+        if (overwrite_lapse(n) /= 0) then
+           
+           call smooth (alp)
+           
+        end if
+        
+        if (overwrite_shift(n) /= 0) then
+           
+           if (shift_state /= 0) then
+              call smooth (betax)
+              call smooth (betay)
+              call smooth (betaz)
+           end if
+           
+        end if
+        
+     end do
+     
+     if (overwrite_geometry(n) /= 0) then
+        
+        if (conformal_state >= 1) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac")
+        end if
+        if (conformal_state >= 2) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_1derivs")
+        end if
+        if (conformal_state >= 3) then
+           call CCTK_SyncGroup (ierr, cctkGH, "StaticConformal::confac_2derivs")
+        end if
+        
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::metric")
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::curv")
+        
+     end if
+     
+     if (overwrite_lapse(n) /= 0) then
+        
+        call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::lapse")
+        
+     end if
+     
+     if (overwrite_shift(n) /= 0) then
+        
+        if (shift_state /= 0) then
+           call CCTK_SyncGroup (ierr, cctkGH, "ADMBase::shift")
+        end if
+        
+     end if
+     
+  end do
+  
+contains
+  
+  subroutine smooth (val)
+    CCTK_REAL, intent(inout) :: val(:,:,:)
+    
+    where (mask > 0)
+       val =  (1 - mask * smoothing_factor) * val &
+            & + (mask * smoothing_factor / 6) &
+            &   * (+ eoshift(val, +1, dim=1) + eoshift(val, -1, dim=1) &
+            &      + eoshift(val, +1, dim=2) + eoshift(val, -1, dim=2) &
+            &      + eoshift(val, +1, dim=3) + eoshift(val, -1, dim=3))
+    end where
+  end subroutine smooth
+  
+end subroutine NoExcision_Smooth