This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/Boundary/trunk@3 6a38eb6e-646e-4a02-a296-d141613ad6c4

9 files changed, 816 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..09c425e
--- /dev/null
+++ b/README
@@ -0,0 +1,56 @@
+Cactus Code Thorn Boundary
+Authors    : Gabrielle Allen, Gerd Lanfermann
+Managed by : ... <cactusmaint@cactuscode.org>
+Version    : 1  
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+This thorn provides some standard outer boundary conditions. 
+
+Currently:
+
+  - flat boundary conditions (a copy of the point just inside the
+    boundary)
+ 
+  - radiation boundary conditions
+
+2. Dependencies of the thorn
+
+This thorn additionally requires:
+
+implementations:
+	driver 
+thorns: 
+ 	CartGrid3D (needs to get at a GH extension)
+
+3. Thorn distribution
+
+This thorn is publically available
+
+4. Additional information
+
+This thorn only currently works with a 3D Cartesian Grid
+The boundary conditions can be called from Fortran or C, 
+and can be passed either groups of grid functions, or single
+grid variables.
+The flat boundary conditions work with any size stencil width,
+the radiation boundary conditions only currently work with 
+a stencil width of one.
+
+To apply a flat boundary condition from Fortran
+
+	call ApplyFlatBC(cctkGH,sw,"imp::group_or_var")
+
+to apply a radiative boundary condition from Fortran
+
+	call ApplyRadiativeBC(cctkGH,num,sw,"imp::group_or_var","imp::group_or_var2")
+
+Where:
+ stencil width is a dimension 3 integer array containing the stencil 
+ width in each direction
+ num is a coefficient in the implemented radiative boundary condition
+ that is usually one. 
+
+See the documentation for more details.
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..3c6a779
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,4 @@
+# Interface definition for thorn Boundary
+# $Header$
+
+implements: boundary
\ No newline at end of file
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..9645d46
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,2 @@
+# Parameter definitions for thorn Boundary
+# $Header$
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..1af3345
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,2 @@
+# Schedule definitions for thorn Boundary
+# $Header$
diff --git a/src/FlatBoundary.F b/src/FlatBoundary.F
new file mode 100644
index 0000000..779220d
--- /dev/null
+++ b/src/FlatBoundary.F
@@ -0,0 +1,101 @@
+/*@@
+   @file      FlatBoundary.F
+   @date      Mon Mar 15 15:50:14 1999
+   @author    Gerd Lanfermann
+   @desc 
+      File contains the Fortran routine, which perform the 
+      actual GF field assignments for the boundary routines.
+   @enddesc 
+@@*/
+
+/*#define DEBUG_BOUND*/
+
+#include "cctk.h"
+
+/*@@
+   @routine    FortranFlatBC
+   @date       Mon Mar 15 15:51:57 1999
+   @author     Gerd Lanfermann
+   @desc 
+      Routine performs the flat boundary operations.
+      Apply BC, if the BC flag is set to 1. 
+      this is done in FlatBoundaryWrappers.c if
+       * the grid chunk has a physical boundary (bbox)
+       * the size in that direction is gt 1
+       * the lower BC is supposed to be applied (no symmetry)
+      passing: size of GFarray, GF, stencil_width, doBC flags
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+@@*/
+   
+      subroutine FortranFlatBC(nxyz,var,stencil_size,doBC);
+
+      implicit none
+
+      INTEGER nxyz(3),stencil_size(3)
+      CCTK_REAL var(nxyz(1),nxyz(2),nxyz(3))
+      INTEGER doBC(6)
+      INTEGER nx,ny,nz,sw
+
+      nx = nxyz(1)
+      ny = nxyz(2)
+      nz = nxyz(3)
+
+c     LOWER BC
+      if (doBC(1) == 1) then
+
+         do sw = 1, stencil_size(1)
+            var(sw,:,:) = var(stencil_size(1)+1,:,:)
+         end do
+
+      end if
+
+      if (doBC(3) == 1) then
+
+         do sw = 1, stencil_size(2)
+            var(:,sw,:) = var(:,stencil_size(2)+1,:)
+         end do
+
+      endif
+ 
+      if (doBC(5) == 1) then 
+
+         do sw = 1, stencil_size(3)
+            var(:,:,sw) = var(:,:,stencil_size(3)+1)
+         end do
+
+      end if
+
+c     UPPER BC
+      if (doBC(2) == 1 ) then
+         do sw = 0, stencil_size(1)-1
+            var(nx-sw,:,:) = var(nx-stencil_size(1),:,:)
+         end do
+      end if
+
+      if (doBC(4).eq.1) then
+
+         do sw = 0, stencil_size(2)-1
+            var(:,ny-sw,:) = var(:,ny-stencil_size(2),:)
+         end do
+
+      end if
+
+      if (doBC(6).eq.1) then 
+
+         do sw=0, stencil_size(3)-1
+            var(:,:,nz-sw) = var(:,:,nz-stencil_size(3))
+         end do
+
+      end if
+      
+      end subroutine
+
+
+
+
+
diff --git a/src/FlatBoundaryWrappers.c b/src/FlatBoundaryWrappers.c
new file mode 100644
index 0000000..35bc424
--- /dev/null
+++ b/src/FlatBoundaryWrappers.c
@@ -0,0 +1,145 @@
+ /*@@
+   @file      FlatBoundaryWrappers.c
+   @date      Mon Mar 15 15:09:00 1999
+   @author    Gerd Lanfermann
+   @desc 
+     This file contains the C part of the BC routines. They are 
+     called by C directly or by F via the wrappers below.
+   @enddesc 
+ @@*/
+
+/*#define DEBUG_BOUND*/
+
+#include <stdio.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <stdarg.h>
+/*#include <math.h>*/
+#include <string.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "cctk_parameters.h"
+#include "Comm.h"
+#include "Groups.h"
+#include "GHExtensions.h"
+#include "CactusBase/CartGrid3D/src/Symmetry.h"
+#include "Misc.h"
+#include "WarnLevel.h"
+#include "FortranString.h"
+
+ /*@@
+   @routine    ApplyFlatBC
+   @date       Mon Mar 15 15:19:54 1999
+   @author     Gerd Lanfermann
+   @desc 
+     Routine applies the flat (copying) BCs when called from C or 
+     Fortran (via wrapper) with the groupname. The actual GF assignment is
+     carried out in Fortran by FortranFlatBC.
+   @enddesc 
+   @calls      FortranFlatBC
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+/* Call with implementation::group notation, eg: ApplyFlatBC(admgerd::metric) */
+
+void ApplyFlatBC(cGH *GH, int *stencil_size, char *name) {
+
+  DECLARE_CCTK_PARAMETERS
+
+  SymmetryGHex *sGHex; /* the Boundary GHextension, needed for ref*/
+  int gf,first,last,index;
+  int type;            /* type 0 for a group, type 1 for a variable */
+  int num;             /* index number of the group (not the gf)*/
+  int doBC[6]; /* flags if lower/upper BCs are applied (1) or not (0)      */
+               /* indexing as in bbox: 0 xlow, 1 xup, 2 ylow, ...(C!)      */
+               /* this has 3 dimension hardwired! will break for GH->dim>3 */
+
+  /* Get the pointer to the SymmetryGHextension */
+  sGHex  = (SymmetryGHex*)GH->extensions[CCTK_GHExtensionHandle("Symmetry")];
+
+  /* Decide if we have a group or a variable, and get the index */
+  /* type = 1 (group), type = 0 (var), type = -1 (neither) */
+  num  = CCTK_GroupIndex(name);
+
+  if (num < 0) 
+  {
+    num = CCTK_VarIndex(name);
+    if (num > 0)
+    {
+      type = 1; /* Variable */
+    }
+    else
+    {
+      type = -1;
+      CCTK_WARN(0,"Name in ApplyRadiativeBC is neither a group nor a variable");
+    }
+  }
+  else
+  {
+    type = 0; /* Group */
+  }
+
+  /*get the index of the first GF in the group and how many Vars there are*/
+  if (type == 0)
+  {
+    first   = CCTK_FirstVarIndexI(num);
+    last    = first+CCTK_NumVarsInGroupI(num)-1;
+    if (first == -1 || last == -1) 
+      CCTK_WARN(0,"Invalid group number used");
+  }
+  else
+  {
+    first = num;
+    last = num;
+  }
+
+  /* loop over the vars in the group */
+  for (index=first; index<=last; index++) {
+
+    /* and check that we actually have a grid function (and not a scalar)*/
+
+    if (CCTK_GroupTypeFromVarI(index)==GROUP_GF) {
+
+      /* The rule is: IF we have no symmetries for the lower grid boundary,
+	 (GFSym==ESYM_NOSYM or ==ESYM_UNSET) AND we have a lower(upper) bound
+	 AND we have enough gridpoints THEN apply BC  */
+
+      doBC[0]=(((sGHex->GFSym[index][0]==GFSYM_NOSYM)||
+	       (sGHex->GFSym[index][0]==GFSYM_UNSET)) &&
+	      GH->cctk_lsh[1]>1 && GH->cctk_bbox[0]);
+      doBC[2]=(((sGHex->GFSym[index][1]==GFSYM_NOSYM)||
+	       (sGHex->GFSym[index][1]==GFSYM_UNSET)) &&
+	      GH->cctk_lsh[1]>1 && GH->cctk_bbox[2]);
+      doBC[4]=(((sGHex->GFSym[index][2]==GFSYM_NOSYM)||
+	       (sGHex->GFSym[index][2]==GFSYM_UNSET)) &&
+	      GH->cctk_lsh[2]>1 && GH->cctk_bbox[4]);
+         
+      doBC[1] = GH->cctk_lsh[0]>1 && GH->cctk_bbox[1];
+      doBC[3] = GH->cctk_lsh[1]>1 && GH->cctk_bbox[3];
+      doBC[5] = GH->cctk_lsh[2]>1 && GH->cctk_bbox[5];
+
+      FORTRAN_NAME(FortranFlatBC)(
+	  (GH->cctk_lsh),        /* xyz-size of PE local grid */
+	  (GH->data[index][0]),  /* pointer to start of data array GF[index]*/
+	  stencil_size,          /* stencil_width */
+	  (doBC)                 /* 6 dim. array to flag whether to apply low/up BCs */
+      );
+   
+    }
+  }  
+}  
+
+
+void FMODIFIER FORTRAN_NAME(ApplyFlatBC)(cGH *GH, int *stencil_size, ONE_FORTSTRING_ARG) {
+
+  ONE_FORTSTRING_CREATE(name)
+  ApplyFlatBC(GH,stencil_size,name);
+  free(name);
+
+} 
+
diff --git a/src/RadiationBoundary.F b/src/RadiationBoundary.F
new file mode 100644
index 0000000..b317f2b
--- /dev/null
+++ b/src/RadiationBoundary.F
@@ -0,0 +1,277 @@
+/*@@
+   @file      RadiationBoundary.F
+   @date      Mon Mar 15 15:50:14 1999
+   @author    Gerd Lanfermann
+   @desc 
+      File contains the Fortran routine, which perform the 
+      actual GF field assignments for the radiation boundary 
+      routines.
+   @enddesc 
+@@*/
+
+#include "cctk.h"
+
+#define DEBUG_BOUND
+
+/*@@
+   @routine    FortranRadiativeBC
+   @date       Mon Mar 15 15:51:57 1999
+   @author     Miguel Alcubierre
+   @desc Routine performs the radiative boundary operations.
+      
+   @enddesc 
+   @calls     
+   @calledby   
+   @history Cactus 4.0 by Gerd Lanfermann
+ 
+   @endhistory 
+@@*/ 
+      subroutine FortranRadiativeBC(nxyz,dxyz,dt,var_n,var_p,x,y,z,r,
+     $     sw,doBC,var0)
+      implicit none
+      INTEGER nxyz(3)
+      CCTK_REAL    dxyz(3)
+      CCTK_REAL    dt
+      CCTK_REAL    var_n(nxyz(1),nxyz(2),nxyz(3))
+      CCTK_REAL    var_p(nxyz(1),nxyz(2),nxyz(3))
+      CCTK_REAL    x(nxyz(1),nxyz(2),nxyz(3))
+      CCTK_REAL    y(nxyz(1),nxyz(2),nxyz(3))
+      CCTK_REAL    z(nxyz(1),nxyz(2),nxyz(3))
+      CCTK_REAL    r(nxyz(1),nxyz(2),nxyz(3))
+      INTEGER sw(3)
+      INTEGER doBC(6)
+      CCTK_REAL    var0
+      
+      CCTK_REAL    half,one,two,idx,idy,idz,rhox,rhoy,rhoz,dth,dx,dy,dz
+      INTEGER nx,ny,nz
+      
+c     we keep close to the variables used in Cactus3.2
+      nx=nxyz(1)
+      ny=nxyz(2)
+      nz=nxyz(3)
+      
+      dx=dxyz(1)
+      dy=dxyz(2)
+      dz=dxyz(3)
+
+c     Numbers.
+
+      half = 0.5D0
+      one  = 1.0D0
+      two  = 2.0D0
+
+      idx = one/dx
+      idy = one/dy
+      idz = one/dz
+
+c     Find Courant parameters.
+
+      rhox = dt*idx
+      rhoy = dt*idy
+      rhoz = dt*idz
+
+      dth = half*dt
+
+      if (sw(1).ne.1 .or. sw(2).ne.1 .or. sw(3).ne.1) then 
+         call CCTK_WARN(0,"Radiative BC only with stencil_size=1")
+      endif
+      
+
+c     **** LOWER X-bound
+      if (doBC(1).eq.1) then
+
+c     Second order radial wave boundaries.
+         var_n(1,:,:) = (dt*var0*(x(1,:,:)/r(1,:,:)**2
+     .        + x(2,:,:)/r(2,:,:)**2)
+     .        - var_n(2,:,:)*(rhox + x(2,:,:)/r(2,:,:)
+     .        *(one + dth/r(2,:,:)))
+     .        + var_p(1,:,:)*(rhox + x(1,:,:)/r(1,:,:)
+     .        *(one - dth/r(1,:,:)))
+     .        - var_p(2,:,:)*(rhox - x(2,:,:)/r(2,:,:)
+     .        *(one - dth/r(2,:,:))))
+     .        /(-rhox + x(1,:,:)/r(1,:,:)*(one + dth/r(1,:,:)))
+
+c     Second order  normal wave  boundaries.
+c         var_n(1,:,:) = (-dt*var0*(x(1,:,:)/r(1,:,:)**2
+c     .       + x(2,:,:)/r(2,:,:)**2)
+c     .       + var_n(2,:,:)*(rhox - one
+c     .       + dth*x(2,:,:)/r(2,:,:)**2)
+c     .       + var_p(1,:,:)*(one - rhox
+c     .       + dth*x(1,:,:)/r(1,:,:)**2)
+c     .       + var_p(2,:,:)*(one + rhox
+c     .       + dth*x(2,:,:)/r(2,:,:)**2))
+c     .       / (one + rhox - dth*x(1,:,:)/r(1,:,:)**2)
+c     First order  normal wave  boundaries.
+c     
+c     var_n(1,:,:) = (x(1,:,:)*dt*var0 - r(1,:,:)**2
+c     .                         *(var_p(1,:,:) + rhox*var_n(2,:,:)))
+c     .                         /(x(1,:,:)*dt - r(1,:,:)**2*(one+rhox))
+
+      end if
+
+
+c     **** LOWER X-bound
+      if (doBC(3).eq.1) then
+c     Second order  radial wave  boundaries.
+
+         var_n(:,1,:) = (dt*var0*(y(:,1,:)/r(:,1,:)**2
+     .        + y(:,2,:)/r(:,2,:)**2)
+     .        - var_n(:,2,:)*(rhoy + y(:,2,:)/r(:,2,:)
+     .        *(one + dth/r(:,2,:)))
+     .        + var_p(:,1,:)*(rhoy + y(:,1,:)/r(:,1,:)
+     .        *(one - dth/r(:,1,:)))
+     .        - var_p(:,2,:)*(rhoy - y(:,2,:)/r(:,2,:)
+     .        *(one - dth/r(:,2,:))))
+     .        /(-rhoy + y(:,1,:)/r(:,1,:)*(one + dth/r(:,1,:)))
+         
+c     Second order  normal wave  boundaries.
+c                 var_n(:,1,:) = (-dt*var0*(y(:,1,:)/r(:,1,:)**2
+c     .                         + y(:,2,:)/r(:,2,:)**2)
+c     .                         + var_n(:,2,:)*(rhoy - one
+c     .                         + dth*y(:,2,:)/r(:,2,:)**2)
+c     .                         + var_p(:,1,:)*(one - rhoy
+c     .                         + dth*y(:,1,:)/r(:,1,:)**2)
+c     .                         + var_p(:,2,:)*(one + rhoy
+c     .                         + dth*y(:,2,:)/r(:,2,:)**2))
+c     .                         / (one + rhoy - dth*y(:,1,:)/r(:,1,:)**2)
+c     First order  normal wave  boundaries.
+c
+c                 var_n(:,1,:) = (y(:,1,:)*dt*var0 - r(:,1,:)**2
+c     .                         *(var_p(:,1,:) + rhoy*var_n(:,2,:)))
+c     .                         /(y(:,1,:)*dt - r(:,1,:)**2*(one+rhoy))
+      end if
+
+
+
+
+c     **** LOWER Z-bound
+      if (doBC(5).eq.1) then
+c     Second order  radial wave  boundaries.
+
+         var_n(:,:,1) = (dt*var0*(z(:,:,1)/r(:,:,1)**2
+     .        + z(:,:,2)/r(:,:,2)**2)
+     .        - var_n(:,:,2)*(rhoz + z(:,:,2)/r(:,:,2)
+     .        *(one + dth/r(:,:,2)))
+     .        + var_p(:,:,1)*(rhoz + z(:,:,1)/r(:,:,1)
+     .        *(one - dth/r(:,:,1)))
+     .        - var_p(:,:,2)*(rhoz - z(:,:,2)/r(:,:,2)
+     .        *(one - dth/r(:,:,2))))
+     .        /(-rhoz + z(:,:,1)/r(:,:,1)*(one + dth/r(:,:,1)))
+
+c     Second order  normal wave  boundaries.
+c                 var_n(:,:,1) = (-dt*var0*(z(:,:,1)/r(:,:,1)**2
+c     .                         + z(:,:,2)/r(:,:,2)**2)
+c     .                         + var_n(:,:,2)*(rhoz - one
+c     .                         + dth*z(:,:,2)/r(:,:,2)**2)
+c     .                         + var_p(:,:,1)*(one - rhoz
+c     .                         + dth*z(:,:,1)/r(:,:,1)**2)
+c     .                         + var_p(:,:,2)*(one + rhoz
+c     .                         + dth*z(:,:,2)/r(:,:,2)**2))
+c     .                         / (one + rhoz - dth*z(:,:,1)/r(:,:,1)**2)
+c     First order  normal wave  boundaries.
+c                 var_n(:,:,1) = (z(:,:,1)*dt*var0 - r(:,:,1)**2
+c     .                         *(var_p(:,:,1) + rhoz*var_n(:,:,2)))
+c     .                         /(z(:,:,1)*dt - r(:,:,1)**2*(one+rhoz))
+
+      end if
+
+
+c     **** UPPER X-bound
+      if (doBC(2).eq.1) then
+
+c     Second order  radial wave  boundaries.
+         var_n(nx,:,:) = (dt*var0*(x(nx,:,:)/r(nx,:,:)**2
+     .        + x(nx-1,:,:)/r(nx-1,:,:)**2)
+     .        + var_n(nx-1,:,:)*(rhox - x(nx-1,:,:)/r(nx-1,:,:)
+     .        *(one + dth/r(nx-1,:,:)))
+     .        + var_p(nx,:,:)*(-rhox + x(nx,:,:)/r(nx,:,:)
+     .        *(one - dth/r(nx,:,:)))
+     .        + var_p(nx-1,:,:)*(rhox + x(nx-1,:,:)/r(nx-1,:,:)
+     .        *(one - dth/r(nx-1,:,:))))
+     .        /(rhox + x(nx,:,:)/r(nx,:,:)*(one + dth/r(nx,:,:)))
+         
+c     Second order  normal wave  boundaries.
+c              var_n(nx,:,:) = (dt*var0*(x(nx,:,:)/r(nx,:,:)**2
+c    .                       + x(nx-1,:,:)/r(nx-1,:,:)**2)
+c    .                       + var_n(nx-1,:,:)*(rhox - one
+c    .                       - dth*x(nx-1,:,:)/r(nx-1,:,:)**2)
+c    .                       + var_p(nx,:,:)*(one - rhox
+c    .                       - dth*x(nx,:,:)/r(nx,:,:)**2)
+c    .                       + var_p(nx-1,:,:)*(one + rhox
+c    .                       - dth*x(nx-1,:,:)/r(nx-1,:,:)**2))
+c    .                       / (one + rhox + dth*x(nx,:,:)/r(nx,:,:)**2)
+c     First order  normal wave  boundaries.
+c
+c              var_n(nx,:,:) = (x(nx,:,:)*dt*var0 + r(nx,:,:)**2
+c     .                       *(var_p(nx,:,:) + rhox*var_n(nx-1,:,:)))
+c     .                       /(x(nx,:,:)*dt + r(nx,:,:)**2*(one+rhox))
+         
+      end if
+
+
+c     **** UPPER Y-bound
+      if(doBC(4).eq.1) then
+
+c     Second order  radial wave  boundaries.
+         var_n(:,ny,:) = (dt*var0*(y(:,ny,:)/r(:,ny,:)**2
+     .        + y(:,ny-1,:)/r(:,ny-1,:)**2)
+     .        + var_n(:,ny-1,:)*(rhoy - y(:,ny-1,:)/r(:,ny-1,:)
+     .        *(one + dth/r(:,ny-1,:)))
+     .        + var_p(:,ny,:)*(-rhoy + y(:,ny,:)/r(:,ny,:)
+     .        *(one - dth/r(:,ny,:)))
+     .        + var_p(:,ny-1,:)*(rhoy + y(:,ny-1,:)/r(:,ny-1,:)
+     .        *(one - dth/r(:,ny-1,:))))
+     .        /(rhoy + y(:,ny,:)/r(:,ny,:)*(one + dth/r(:,ny,:)))
+         
+c     Second order  normal wave  boundaries.
+c             var_n(:,ny,:) = (dt*var0*(y(:,ny,:)/r(:,ny,:)**2
+c     .                       + y(:,ny-1,:)/r(:,ny-1,:)**2)
+c     .                       + var_n(:,ny-1,:)*(rhoy - one
+c     .                       - dth*y(:,ny-1,:)/r(:,ny-1,:)**2)
+c     .                       + var_p(:,ny,:)*(one - rhoy
+c     .                       - dth*y(:,ny,:)/r(:,ny,:)**2)
+c     .                       + var_p(:,ny-1,:)*(one + rhoy
+c     .                       - dth*y(:,ny-1,:)/r(:,ny-1,:)**2))
+c     .                       / (one + rhoy + dth*y(:,ny,:)/r(:,ny,:)**2)
+c     First order  normal wave  boundaries.
+c
+c              var_n(:,ny,:) = (y(:,ny,:)*dt*var0 + r(:,ny,:)**2
+c     .                       *(var_p(:,ny,:) + rhoy*var_n(:,ny-1,:)))
+c     .                       /(y(:,ny,:)*dt + r(:,ny,:)**2*(one+rhoy))
+
+      end if
+
+c     **** UPPER Z-bound
+      if (doBC(6).eq.1) then  
+
+c     Second order  radial wave  boundaries.
+         var_n(:,:,nz) = (dt*var0*(z(:,:,nz)/r(:,:,nz)**2
+     .        + z(:,:,nz-1)/r(:,:,nz-1)**2)
+     .        + var_n(:,:,nz-1)*(rhoz - z(:,:,nz-1)/r(:,:,nz-1)
+     .        *(one + dth/r(:,:,nz-1)))
+     .        + var_p(:,:,nz)*(-rhoz + z(:,:,nz)/r(:,:,nz)
+     .        *(one - dth/r(:,:,nz)))
+     .        + var_p(:,:,nz-1)*(rhoz + z(:,:,nz-1)/r(:,:,nz-1)
+     .        *(one - dth/r(:,:,nz-1))))
+     .        /(rhoz + z(:,:,nz)/r(:,:,nz)*(one + dth/r(:,:,nz)))
+         
+c     Second order  normal wave  boundaries.
+c             var_n(:,:,nz) = (dt*var0*(z(:,:,nz)/r(:,:,nz)**2
+c     .                       + z(:,:,nz-1)/r(:,:,nz-1)**2)
+c     .                       + var_n(:,:,nz-1)*(rhoz - one
+c     .                       - dth*z(:,:,nz-1)/r(:,:,nz-1)**2)
+c     .                       + var_p(:,:,nz)*(one - rhoz
+c     .                       - dth*z(:,:,nz)/r(:,:,nz)**2)
+c     .                       + var_p(:,:,nz-1)*(one + rhoz
+c     .                       - dth*z(:,:,nz-1)/r(:,:,nz-1)**2))
+c     .                       / (one + rhoz + dth*z(:,:,nz)/r(:,:,nz)**2)
+         
+c     First order  normal wave  boundaries.
+c
+c              var_n(:,:,nz) = (z(:,:,nz)*dt*var0 + r(:,:,nz)**2
+c     .                       *(var_p(:,:,nz) + rhoz*var_n(:,:,nz-1)))
+c     .                       /(z(:,:,nz)*dt + r(:,:,nz)**2*(one+rhoz))
+      end if
+
+      return
+      end subroutine
diff --git a/src/RadiationBoundaryWrappers.c b/src/RadiationBoundaryWrappers.c
new file mode 100644
index 0000000..e5c7bf5
--- /dev/null
+++ b/src/RadiationBoundaryWrappers.c
@@ -0,0 +1,217 @@
+ /*@@
+   @file      RadiationBoundaryWrappers.c
+   @date      Mon Mar 15 15:09:00 1999
+   @author    Gerd Lanfermann
+   @desc 
+     This file contains the C part of the BC routines. They are 
+     called by C directly or by F via the wrappers below.
+   @enddesc 
+ @@*/
+
+/*#define DEBUG_BOUND*/
+
+#include <stdio.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <stdarg.h>
+/*#include <math.h>*/
+#include <string.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "cctk_parameters.h"
+#include "Comm.h"
+#include "Coord.h"
+#include "Groups.h"
+#include "GHExtensions.h"
+#include "CactusBase/CartGrid3D/src/Symmetry.h"
+#include "Misc.h"
+#include "WarnLevel.h"
+#include "FortranString.h"
+
+ /*@@
+   @routine    ApplyRadiativeBC
+   @date       Mon Mar 15 15:21:45 1999
+   @author     Gerd Lanfermann
+   @desc 
+     Radiative BC require GFs at past timesteps:
+     Call with implementation::group and implementation_group_previous
+     eg: ApplyFlatBC(admgerd::metric,admgerd::metric_previous) 
+     It is not checked, that these two groups make sense!
+     Actual assigment is carried out in Fortran by  FortranRadiativeBC
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void ApplyRadiativeBC(cGH *GH, CCTK_REAL var0, int *sw, char *imp_group, char *imp_group_p) {
+  DECLARE_CCTK_PARAMETERS
+
+  SymmetryGHex *sGHex;  /* the symmetry boundary GHextension, needed for ref*/
+  int xi,yi,zi,ri;           /* index of the cartesian coordinate fields */
+  int first,last,index,first_p,last_p;
+  int num,num_p;   /* index number of the group (not the gf)*/
+  int doBC[6]; /* flags if lower/upper BCs are applied (1) or not (0)      */
+               /* indexing as in bbox: 0 xlow, 1 xup, 2 ylow, ...(C!)      */
+               /* this has 3 dimension hardwired! will break for GH->dim>3 */
+  int type;    /* 0 groups; 1 variables; -1 neither */
+
+  /* Get the pointer to the SymmetryGHextension */
+  sGHex  = (SymmetryGHex*)malloc(sizeof(SymmetryGHex));
+  sGHex  = (SymmetryGHex*)GH->extensions[CCTK_GHExtensionHandle("Symmetry")];
+
+  /* Decide if we have a group or a variable, and get the index */
+  /* type = 1 (group), type = 0 (var), type = -1 (neither) */
+  num  = CCTK_GroupIndex(imp_group);
+
+  if (num < 0) 
+  {
+    num = CCTK_VarIndex(imp_group);
+    if (num > 0)
+    {
+      type = 1; /* Variable */
+    }
+    else
+    {
+      type = -1;
+      CCTK_WARN(0,"Name in ApplyRadiativeBC is neither a group nor a variable");
+    }
+  }
+  else
+  {
+    type = 0; /* Group */
+  }
+
+  num_p= CCTK_GroupIndex(imp_group_p);
+  if (num_p < 0) 
+  {
+    num_p = CCTK_VarIndex(imp_group_p);
+    if (num_p > 0)
+    {
+      if (type != 1) /* Variable */
+      {
+	CCTK_WARN(0,"Group and Variable given in ApplyRadiativeBC");
+      }
+    }
+    else
+    {
+      CCTK_WARN(0,"Name in ApplyRadiativeBC is neither a group nor a variable");
+    }
+  }
+  else
+  {
+    if (type != 0) /* Group */
+    {
+      CCTK_WARN(0,"Group and Variable given in ApplyRadiativeBC");
+    }
+  }
+
+  /*get the index (offset) of the first GF in the group and how many Vars there are*/
+  if (type == 0)
+  {
+    first   = CCTK_FirstVarIndexI(num);
+    last    = first+CCTK_NumVarsInGroupI(num)-1;
+    first_p = CCTK_FirstVarIndexI(num_p);
+    last_p  = first_p+CCTK_NumVarsInGroupI(num_p)-1;
+    if (first == -1 || last == -1 || first_p == -1) 
+      CCTK_WARN(0,"Invalid group number used");
+    if (last-first != last_p-first_p)
+      CCTK_WARN(0,"Groups have different number of variables in ApplyRadiativeBC");
+  }
+  else
+  {
+    first = num;
+    last = num;
+    first_p = num_p;
+    last_p = num_p;
+  }
+
+#ifdef DEBUG_BOUND
+  CCTK_PRINTSEPARATOR
+  printf("In RadiationBoundaryWrappers\n----------------------------\n");
+  if (type == 0)
+  {
+    printf("Group: %s, has indices %d-%d\n",CCTK_GroupName(num),first,last);
+    printf("Group: %s, has indices %d-%d\n",CCTK_GroupName(num_p),first_p,last_p);
+  }
+  else
+  {
+    printf("Variable: %s\n",CCTK_VarName(num));
+    printf("Variable: %s\n",CCTK_VarName(num_p));
+  }
+
+#endif
+
+  /* radiative needs reference to the XYZR datafield defined by Cartesian grid 
+     and we get the index for those fields here: 
+     (I wonder how this will work with multiple grids). */
+  xi = CCTK_CoordIndex("x");
+  yi = CCTK_CoordIndex("y");
+  zi = CCTK_CoordIndex("z");
+  ri = CCTK_CoordIndex("r");
+
+  /* loop over the vars in the group (index + offset) */
+  for (index=0;index<last-first+1;index++) {
+    /* and check that we actually have a grid function (and not a scalar)*/
+    if (CCTK_GroupTypeFromVarI(first+index) == GROUP_GF) {
+
+      /* The rule is: IF we have no symmetries for the lower grid boundary,
+	 (GFSym==GFSYM_NOSYM or ==GFSYM_UNSET) AND we have a lower(upper) bound
+	 AND we have enough gridpoints THEN apply BC  */
+      doBC[0]=(((sGHex->GFSym[first+index][0]==GFSYM_NOSYM)||
+	       (sGHex->GFSym[first+index][0]==GFSYM_UNSET)) &&
+	      GH->cctk_lsh[1]>1 && GH->cctk_bbox[0]);
+      doBC[2]=(((sGHex->GFSym[first+index][1]==GFSYM_NOSYM)||
+	       (sGHex->GFSym[first+index][1]==GFSYM_UNSET)) &&
+	      GH->cctk_lsh[1]>1 && GH->cctk_bbox[2]);
+      doBC[4]=(((sGHex->GFSym[first+index][2]==GFSYM_NOSYM)||
+	       (sGHex->GFSym[first+index][2]==GFSYM_UNSET)) &&
+	      GH->cctk_lsh[2]>1 && GH->cctk_bbox[4]);
+         
+      doBC[1] = GH->cctk_lsh[0]>1 && GH->cctk_bbox[1];
+      doBC[3] = GH->cctk_lsh[1]>1 && GH->cctk_bbox[3];
+      doBC[5] = GH->cctk_lsh[2]>1 && GH->cctk_bbox[5];
+
+#ifdef DEBUG_BOUND
+      printf("Applying radiative boundary condition to %s (%d)\n",CCTK_FullName(first+index),first+index);
+#endif
+
+      /* Now make the call to F-boundary routine, passing all relevant information,
+         and let F do the physics */
+      FORTRAN_NAME(FortranRadiativeBC)(
+               (GH->cctk_lsh),                     /* yzx-size of PE local grid */
+	       (GH->cctk_delta_space),             /* grid spacing */
+	      &(GH->cctk_delta_time),              /* dt */
+	       (GH->data[index+first]  [0]),       /* pointer to start of data array GF[] */
+	       (GH->data[index+first_p][0]),       /* pointer to start of prev data array */
+	        GH->data[xi][0],GH->data[yi][0],   /* pointer to the XYZR fields */
+	        GH->data[zi][0],GH->data[ri][0],
+	        sw,                                /* stencil_width */
+	       (doBC),                             /* 6 dim. array to flag applying low/up BCs */
+	      &(var0));                            /* radiative BC limit ? */
+
+    }
+  }
+
+#ifdef DEBUG_BOUND
+  CCTK_PRINTSEPARATOR
+#endif
+
+
+}  
+    
+void FMODIFIER FORTRAN_NAME(ApplyRadiativeBC)(cGH *GH,CCTK_REAL *var0,int *sw, TWO_FORTSTRINGS_ARGS) {
+
+  TWO_FORTSTRINGS_CREATE(imp_group,imp_group_p)
+
+  ApplyRadiativeBC(GH,*var0,sw,imp_group,imp_group_p);
+  free(imp_group);
+  free(imp_group_p);
+
+}
+
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..4ee7a71
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,12 @@
+# Main make.code.defn file for thorn Boundary
+# $Header$
+
+# Source files in this directory
+SRCS = FlatBoundary.F\
+       FlatBoundaryWrappers.c\
+       RadiationBoundary.F\
+       RadiationBoundaryWrappers.c 
+
+# Subdirectories containing source files
+SUBDIRS = 
+