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

13 files changed, 1926 insertions, 88 deletions

diff --git a/interface.ccl b/interface.ccl
index 0aa4916..968f4cf 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -3,4 +3,98 @@
 
 IMPLEMENTS: NoExcision
 
-INHERITS: ADMBase StaticConformal grid
+INHERITS: ADMBase StaticConformal grid SphericalSurface Boundary
+
+USES INCLUDE: Boundary.h
+
+USES INCLUDE: carpet.hh
+
+CCTK_INT FUNCTION MoLQueryEvolvedRHS (CCTK_INT IN EvolvedIndex)
+USES FUNCTION MoLQueryEvolvedRHS
+
+CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \
+                                            CCTK_INT IN faces, \
+                                            CCTK_INT IN boundary_width, \
+                                            CCTK_INT IN table_handle, \
+                                            CCTK_STRING IN var_name, \
+                                            CCTK_STRING IN bc_name)
+
+USES FUNCTION Boundary_SelectGroupForBC
+
+CCTK_INT smask type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+  nes_mask
+} "mask for smoothing"
+
+CCTK_REAL cg_res_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+  resgxx, resgxy, resgxz, resgyy, resgyz, resgzz
+} "Conjugate Gradient residual for the metric"
+
+CCTK_REAL cg_res_curv type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+  reskxx, reskxy, reskxz, reskyy, reskyz, reskzz
+} "Conjugate Gradient residual for the extrinsic curvature"
+
+CCTK_REAL cg_res_shift type=GF timelevels=1 tags='tensortypealias="u" Prolongation="None"'
+{
+  resx, resy, resz
+} "Conjugate Gradient residual for the shift"
+
+CCTK_REAL cg_res_lapse type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+  res
+} "Conjugate Gradient residual for the lapse"
+
+
+CCTK_REAL cg_d_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+  dgxx, dgxy, dgxz, dgyy, dgyz, dgzz
+} "Conjugate Gradient d for the metric"
+
+CCTK_REAL cg_d_curv type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+  dkxx, dkxy, dkxz, dkyy, dkyz, dkzz
+} "Conjugate Gradient d for the extrinsic curvature"
+
+CCTK_REAL cg_d_shift type=GF timelevels=1 tags='tensortypealias="u" Prolongation="None"'
+{
+  dx, dy, dz
+} "Conjugate Gradient d for the shift"
+
+CCTK_REAL cg_d_lapse type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+  d
+} "Conjugate Gradient d for the lapse"
+
+
+CCTK_REAL cg_q_metric type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+  qgxx, qgxy, qgxz, qgyy, qgyz, qgzz
+} "Conjugate Gradient q for the metric"
+
+CCTK_REAL cg_q_curv type=GF timelevels=1 tags='tensortypealias="dd_sym" Prolongation="None"'
+{
+  qkxx, qkxy, qkxz, qkyy, qkyz, qkzz
+} "Conjugate Gradient q for the extrinsic curvature"
+
+CCTK_REAL cg_q_shift type=GF timelevels=1 tags='tensortypealias="u" Prolongation="None"'
+{
+  qx, qy, qz
+} "Conjugate Gradient q for the shift"
+
+CCTK_REAL cg_q_lapse type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+  q
+} "Conjugate Gradient q for the lapse"
+
+
+CCTK_REAL cg_red_all type=GF timelevels=1 tags='tensortypealias="scalar" Prolongation="None"'
+{
+  redgxx, redgxy, redgxz, redgyy, redgyz, redgzz
+  redkxx, redkxy, redkxz, redkyy, redkyz, redkzz
+  red, redx, redy, redz
+} "Conjugate Gradient red for all variables"
+
+# Control variable for smoothing loop
+CCTK_INT loop_control TYPE=SCALAR
diff --git a/param.ccl b/param.ccl
index 1fe7342..6d93bd3 100644
--- a/param.ccl
+++ b/param.ccl
@@ -5,13 +5,18 @@ BOOLEAN verbose "Produce some screen output"
 {
 } "no"
 
+KEYWORD method "Method to use"
+{
+  "old" :: "Use old method"
+  "new" :: "Use new method"
+} "old"
+
+
 INT num_regions "Number of no-excision regions"
 {
   0:10 :: ""
 } 0
 
-
-
 BOOLEAN overwrite_geometry[10] "Set the geometry to Minkowski"
 {
 } "yes"
@@ -40,27 +45,100 @@ BOOLEAN overwrite_shift[10] "Set the shift to zero"
 
 
 
-REAL centre_x[10] "X-coordinate of the centre of the region"
+REAL centre_x[10] "x-coordinate of the centre of the region"
 {
   (*:*) :: ""
 } 0.0
 
-REAL centre_y[10] "Y-coordinate of the centre of the region"
+REAL centre_y[10] "y-coordinate of the centre of the region"
 {
   (*:*) :: ""
 } 0.0
 
-REAL centre_z[10] "Z-coordinate of the centre of the region"
+REAL centre_z[10] "z-coordinate of the centre of the region"
 {
   (*:*) :: ""
 } 0.0
 
+KEYWORD region_shape[10] "Shape of the region"
+{
+  "sphere"    :: "use radius"
+  "ellipsoid" :: "use radius_x, radius_y, and radius_z"
+  "surface"   :: "use a spherical surface shape"
+} "sphere"
+
 REAL radius[10] "Radius of the region"
 {
   0.0:*) :: ""
 } 1.0
 
-REAL smoothing_zone_width[10] "Width of smoothing zone inside the region"
+REAL radius_x[10] "x-radius of the region"
+{
+  0.0:*) :: ""
+} 1.0
+
+REAL radius_y[10] "y-radius of the region"
 {
   0.0:*) :: ""
+} 1.0
+
+REAL radius_z[10] "z-radius of the region"
+{
+  0.0:*) :: ""
+} 1.0
+
+INT surface_index[10] "Spherical surface index"
+{
+  0:* :: "must be an index of a spherical surface"
+} 0
+
+
+
+BOOLEAN reduce_rhs[10] "Reduce RHS"
+{
+} "no"
+
+REAL reduction_factor[10] "Reduction factor for RHS (0=complete, 1=no reduction)"
+{
+  *:* :: ""
+} 0.0
+
+
+
+KEYWORD smoothing_function[10] "Smoothing function"
+{
+  "linear" :: "linear ramp"
+  "spline" :: "cubic spline ramp"
+  "cosine" :: "cosine ramp"
+} "linear"
+
+REAL smoothing_zone_width[10] "Relative width of smoothing zone inside the region"
+{
+  0.0:1.0 :: ""
 } 0.0
+
+INT smoothing_iterations "Smoothing iterations"
+{
+  0:* :: ""
+} 10
+
+REAL smoothing_factor "Initial moothing factor"
+{
+  (0:2) :: ""
+} 1.2
+
+
+
+BOOLEAN smooth_regions "Smooth overwritten regions?"
+{
+} "no"
+
+INT smoothing_order "Order of the derivatives used for CG smoothing"
+{
+  2:6:2 :: ""
+} 6
+
+REAL smoothing_eps "CG smoothing stop criteria"
+{
+  (0.0:* :: ""
+} 1e-6
diff --git a/schedule.ccl b/schedule.ccl
index ac7b7b1..54d402b 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,7 +1,150 @@
 # Schedule definitions for thorn NoExcision
 # $Header$
 
-SCHEDULE NoExcision_Overwrite AT initial AFTER (ADMBase_InitialData ADMBase_InitialGauge)
+SCHEDULE NoExcision_Overwrite IN ADMBase_PostInitial
 {
   LANG: Fortran
 } "Overwrite regions with Minkowski"
+
+if (CCTK_Equals(method,"new")) {
+  schedule NoExcision_SetSym at CCTK_BASEGRID
+  {
+    LANG: Fortran
+  } "Register the symmetries for the conjugate gradient functions."
+}
+
+if (smooth_regions) {
+  if (CCTK_Equals(method,"old")) {
+    SCHEDULE NoExcision_Smooth IN ADMBase_PostInitial AFTER NoExcision_Overwrite
+    {
+      LANG: Fortran
+    } "Smooth regions"
+  } else {
+
+    SCHEDULE GROUP NoExcision_CGSmoothing IN ADMBase_PostInitial
+    {
+      STORAGE: cg_res_metric, cg_res_curv, cg_res_shift, cg_res_lapse
+      STORAGE: cg_d_metric, cg_d_curv, cg_d_shift, cg_d_lapse
+      STORAGE: cg_q_metric, cg_q_curv, cg_q_shift, cg_q_lapse
+      STORAGE: cg_red_all
+      STORAGE: smask
+      STORAGE: loop_control
+
+    } "Conjugate Gradient smoothing"
+
+    SCHEDULE NoExcision_CGInit_1 IN NoExcision_CGSmoothing
+    {
+      LANG: Fortran
+    } "Initialise the conjugate gradient method 1"
+
+    SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p1 IN NoExcision_CGSmoothing AFTER NoExcision_CGInit_1
+    {
+      LANG: Fortran
+      SYNC: cg_d_metric
+      SYNC: cg_d_curv
+      SYNC: cg_d_shift
+      SYNC: cg_d_lapse
+      SYNC: cg_res_metric
+      SYNC: cg_res_curv
+      SYNC: cg_res_shift
+      SYNC: cg_res_lapse
+      SYNC: cg_red_all
+    } "Select variables for boundary conditions 1"
+
+    SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p1 IN NoExcision_CGSmoothing AFTER NoExcision_CGApplySym_p1
+    {
+    } "Apply boundary conditions (Symmetries) 1"
+
+    SCHEDULE NoExcision_CGInit_2 IN NoExcision_CGSmoothing AFTER NoExcision_CGApplySym_p1
+    {
+      LANG: Fortran
+      OPTIONS: level
+    } "Initialise the conjugate gradient method 2"
+
+    SCHEDULE GROUP NoExcision_Smoothing IN NoExcision_CGSmoothing AFTER NoExcision_CGInit_2 while NoExcision::loop_control
+    {
+    } "Smooth regions"
+
+    SCHEDULE NoExcision_CG_1 IN NoExcision_Smoothing
+    {
+      LANG: Fortran
+    } "Conjugate gradients step 1"
+
+    SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p2 IN NoExcision_Smoothing AFTER NoExcision_CG_1
+    {
+      LANG: Fortran
+      SYNC: cg_q_metric
+      SYNC: cg_q_curv
+      SYNC: cg_q_shift
+      SYNC: cg_q_lapse
+      SYNC: cg_red_all
+    } "Select variables for boundary conditions 2"
+
+    SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p2 IN NoExcision_Smoothing AFTER NoExcision_CGApplySym_p2
+    {
+    } "Apply boundary conditions (Symmetries) 2"
+
+    SCHEDULE NoExcision_CG_2 IN NoExcision_Smoothing AFTER NoExcision_CG_1
+    {
+      LANG: Fortran
+      OPTIONS: level
+    } "Conjugate gradients step 2"
+
+    SCHEDULE NoExcision_CG_3 IN NoExcision_Smoothing
+    {
+      LANG: Fortran
+    } "Conjugate gradients step 3"
+
+    SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p3 IN NoExcision_Smoothing AFTER NoExcision_CG_3
+    {
+      LANG: Fortran
+      SYNC: cg_res_metric
+      SYNC: cg_res_curv
+      SYNC: cg_res_shift
+      SYNC: cg_res_lapse
+      SYNC: metric
+      SYNC: curv
+      SYNC: shift
+      SYNC: lapse
+      SYNC: cg_red_all
+    } "Select variables for boundary conditions 3"
+
+    SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p3 IN NoExcision_Smoothing AFTER NoExcision_CGApplySym_p3
+    {
+    } "Apply boundary conditions (Symmetries) 3"
+
+    SCHEDULE NoExcision_CG_4 IN NoExcision_Smoothing AFTER NoExcision_CG_1
+    {
+      LANG: Fortran
+      OPTIONS: level
+    } "Conjugate gradients step 4"
+
+    SCHEDULE NoExcision_CG_5 IN NoExcision_Smoothing
+    {
+      LANG: Fortran
+    } "Conjugate gradients step 5"
+
+    SCHEDULE NoExcision_CGApplySym AS NoExcision_CGApplySym_p4 IN NoExcision_CGSmoothing AFTER NoExcision_CG_5
+    {
+      LANG: Fortran
+      SYNC: cg_d_metric
+      SYNC: cg_d_curv
+      SYNC: cg_d_shift
+      SYNC: cg_d_lapse
+    } "Select variables for boundary conditions 4"
+
+    SCHEDULE GROUP ApplyBCs AS NoExcision_CGApplyBCs_p4 IN NoExcision_CGSmoothing AFTER NoExcision_CGApplySym_p4
+    {
+    } "Apply boundary conditions (Symmetries) 4"
+  }
+}
+
+#SCHEDULE NoExcision_OverwriteBSSN AT poststep
+#{
+#  LANG: Fortran
+#} "Overwrite regions with Minkowski"
+
+SCHEDULE NoExcision_Reduce IN MoL_PostRHS
+{
+  LANG: C
+} "Reduce RHS"
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