Generalized boundary condition routines for applying to

arbitrary CCTK data types (except CCTK_COMPLEX).

Added/completed grdoc.


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

1 files changed, 1352 insertions, 1110 deletions

diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index 8ebc6d4..dfe0e46 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -1,1318 +1,1560 @@
  /*@@
-   @file      RadiationBoundaryWrappers.c
+   @file      RadiationBoundary.c
    @date      Mon Mar 15 15:09:00 1999
-   @author    Gabrielle Allen, Gerd Lanfermann
-   @desc 
-   Radiation boundary conditions for 3D only callable from Fortran or C
-   @enddesc 
-   @version $Header$
+   @author    Miguel Alcubierre, Gabrielle Allen, Gerd Lanfermann
+   @desc
+              <PRE>
+              Routines for applying radiation boundary conditions
+
+              The radiative boundary condition that is implemented is:
+
+                f  =  f0  +  u(r - v*t) / r  +  h(r + v*t) / r
+
+              That is, I assume outgoing radial waves with a 1/r
+              fall off, and the correct asymptotic value f0, plus
+              I include the possibility of incoming waves as well
+              (these incoming waves should be modeled somehow).
+
+              The condition above leads to the differential equation:
+
+                (x / r) d f  +  v d f  + v x (f - f0) / r^2  =  v x H / r^2
+                  i      t         i        i                      i
+
+              where x_i is the normal direction to the given boundaries,
+              and H = 2 dh(s)/ds.
+
+              So at a given boundary I only worry about derivatives in
+              the normal direction.  Notice that u(r-v*t) has dissapeared,
+              but we still do not know the value of H.
+
+              To get H what I do is the following:  I evaluate the
+              expression one point in from the boundary and solve for H
+              there.  We now need a way of extrapolation H to the boundary.
+              For this I assume that H falls off as a power law:
+
+                H = k/r**n   =>  d H  =  - n H/r
+                                  i
+
+              The value of n is is defined by the parameter "radpower".
+              If this parameter is negative, H is forced to be zero (this
+              corresponds to pure outgoing waves and is the default).
+              <P>
+              The behaviour I have observed is the following:  Using H=0
+              is very stable, but has a very bad initial transient. Taking
+              n to be 0 or positive improves the initial behaviour considerably,
+              but introduces a drift that can kill the evolution at very late
+              times.  Empirically, the best value I have found is n=2, for
+              which the initial behaviour is very nice, and the late time drift
+              is quite small.
+
+              Another problem with this condition is that it does not
+              use the physical characteristic speed, but rather it assumes
+              a wave speed of v, so the boundaries should be out in
+              the region where the characteristic speed is constant.
+              Notice that this speed does not have to be 1.  For gauge
+              quantities {alpha,phi,trK} we can have a different asymptotic
+              speed, which is why the value of v is passed as a parameter.
+              </PRE>
+   @enddesc
+   @history
+   @hdate     unknown
+   @hauthor   Gerd Lanfermann
+   @hdesc     Ported to Cactus 4.0
+   @hdate     Fri 6 Apr 2001
+   @hauthor   Thomas Radke
+   @hdesc     BC routines generalized for applying to arbitrary CCTK data types
+   @endhistory
+   @version   $Id$
  @@*/
 
-/*#define DEBUG_BOUNDARY*/
- 
-#include <stdio.h>
-#include <assert.h>
 #include <stdlib.h>
-#include <ctype.h>
-#include <stdarg.h>
 #include <string.h>
 
 #include "cctk.h"
 #include "cctk_FortranString.h"
 #include "cctk_Parameters.h"
-#include "cctk_Stagger.h"
 
-#include "BoundarySymmetries.h"
+#include "Symmetry.h"
 #include "Boundary.h"
 
-#define SQR(a) ((a)*(a))
-
+/* the rcs ID and its dummy function to use it */
 static char *rcsid = "$Header$";
-
-CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundaryWrappers_c)
-
-/* Internal routine prototypes */
-
-static int BndApplyRadiative3Di(cGH *GH, 
-				int doBC, 
-				int *lssh,
-				int *stencil_size,  
-				CCTK_REAL *dxyz, 
-				CCTK_REAL dt, 
-				CCTK_REAL *var_n, 
-				CCTK_REAL *var_p,
-				CCTK_REAL *x, 
-				CCTK_REAL *y, 
-				CCTK_REAL *z, 
-				CCTK_REAL *r,
-				CCTK_REAL var0, 
-				CCTK_REAL v0) ;
-
-static int ApplyBndRadiative(cGH *GH, 
-			     int stencil,
-			     int *stencil_array,
-			     CCTK_REAL var0, 
-			     CCTK_REAL v0, 
-			     int var_current, 
-			     int var_previous,
-			     int num_vars,
-			     int dir);
-
-/* Local variables */
+CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundary_c)
 
 
 /********************************************************************
  ********************    External Routines   ************************
  ********************************************************************/
+/* prototypes for external C routines are declared in header Boundary.h
+   here only follow the fortran wrapper prototypes */
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirGI)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *gi_to,
+                            const int *gi_from);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeGI)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *gi_to,
+                            const int *gi_from);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirGN)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeGN)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirVI)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *vi_to,
+                            const int *vi_from);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeVI)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *vi_to,
+                            const int *vi_from);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirVN)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS);
+void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS);
+
+
+/********************************************************************
+ ********************    Internal Routines   ************************
+ ********************************************************************/
+static int ApplyBndRadiative (cGH *GH,
+                              int stencil_dir,
+                              int stencil_alldirs[],
+                              int dir,
+                              CCTK_REAL var0,
+                              CCTK_REAL speed,
+                              int first_var_to,
+                              int first_var_from,
+                              int num_vars);
+
 
 /*@@
    @routine    BndRadiativeDirGI
    @date       Sat Jan 20 2001
    @author     Gabrielle Allen
-   @desc 
-      Aply radiative BC's by group index in given direction
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Aply radiative BC's by group index in given direction
+   @enddesc
+   @calls      ApplyBndRadiative
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil_size
+   @vdesc      stencil size in this direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to apply BC
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        gi_to
+   @vdesc      index of group to apply BC to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        gi_from
+   @vdesc      index of group to apply BC from
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine ApplyBndRadiative <BR>
+               -1 if invalid group indices are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeDirGI(cGH *GH, 
-		   int stencil_size,
-		   int dir,
-		   CCTK_REAL var0, 
-		   CCTK_REAL speed, 
-		   int gi, 
-		   int gi_p) 
+int BndRadiativeDirGI (cGH *GH,
+                       int stencil_size,
+                       int dir,
+                       CCTK_REAL var0,
+                       CCTK_REAL speed,
+                       int gi_to,
+                       int gi_from)
 {
-  int numvars, first_vi, first_vi_p; 
+  int first_vi_to, first_vi_from, retval;
 
-  first_vi   = CCTK_FirstVarIndexI(gi);
-  first_vi_p = CCTK_FirstVarIndexI(gi_p);
-  numvars    = CCTK_NumVarsInGroupI(gi);
 
-  return ApplyBndRadiative(GH,stencil_size,NULL,var0,
-			   speed,first_vi,first_vi_p,numvars,dir);
+  first_vi_to   = CCTK_FirstVarIndexI (gi_to);
+  first_vi_from = CCTK_FirstVarIndexI (gi_from);
+  if (first_vi_to >= 0 && first_vi_from >= 0)
+  {
+    retval = ApplyBndRadiative (GH, stencil_size, NULL, dir, var0, speed,
+                                first_vi_to, first_vi_from,
+                                CCTK_NumVarsInGroupI (gi_to));
+  }
+  else
+  {
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid group indices %d and/or %d in BndRadiativeDirGI",
+                gi_to, gi_from);
+    retval = -1;
+  }
 
+  return (retval);
 }
 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeDirGI)
-     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
-      CCTK_REAL *speed, int *gi, int *gi_p)
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirGI)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *gi_to,
+                            const int *gi_from)
 {
-  *ierr = BndRadiativeDirGI(GH, *stencil_size, *dir, *var0, *speed, *gi, *gi_p);
+  *ierr = BndRadiativeDirGI (GH, *stencil_size, *dir, *var0, *speed,
+                             *gi_to, *gi_from);
 }
 
+
 /*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
-   @desc 
-      Aply radiative BC's by group index
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Aply radiative BC's by group index
+   @enddesc
+   @calls      ApplyBndRadiative
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil
+   @vdesc      stencil width
+   @vtype      int [ dimension of group ]
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        gi_to
+   @vdesc      index of group to apply BC to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        gi_from
+   @vdesc      index of group to apply BC from
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine ApplyBndRadiative <BR>
+               -1 if invalid group indices are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeGI(cGH *GH, 
-		   int *stencil_array,
-		   CCTK_REAL var0, 
-		   CCTK_REAL speed, 
-		   int gi, 
-		   int gi_p) 
+int BndRadiativeGI (cGH *GH,
+                    int stencil[],
+                    CCTK_REAL var0,
+                    CCTK_REAL speed,
+                    int gi_to,
+                    int gi_from)
 {
-  int numvars, first_vi, first_vi_p; 
+  int first_vi_to, first_vi_from, retval;
 
-  first_vi   = CCTK_FirstVarIndexI(gi);
-  first_vi_p = CCTK_FirstVarIndexI(gi_p);
-  numvars    = CCTK_NumVarsInGroupI(gi);
 
-  return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,first_vi,first_vi_p,numvars,0);
+  first_vi_to   = CCTK_FirstVarIndexI (gi_to);
+  first_vi_from = CCTK_FirstVarIndexI (gi_from);
+  if (first_vi_to >= 0 && first_vi_from >= 0)
+  {
+    retval = ApplyBndRadiative (GH, -1, stencil, 0, var0, speed,
+                                first_vi_to, first_vi_from,
+                                CCTK_NumVarsInGroupI (gi_to));
+  }
+  else
+  {
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid group indices %d and/or %d in BndRadiativeGI",
+                gi_to, gi_from);
+    retval = -1;
+  }
 
+  return (retval);
 }
 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeGI)
-     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
-      CCTK_REAL *speed, int *gi, int *gi_p)
+void CCTK_FCALL CCTK_FNAME (BndRadiativeGI)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *gi_to,
+                            const int *gi_from)
 {
-  *ierr = BndRadiativeGI(GH, stencil_array, *var0, *speed, *gi, *gi_p);
+  *ierr = BndRadiativeGI (GH, stencil, *var0, *speed, *gi_to, *gi_from);
 }
 
 
 /* ===================================================================== */
 
-
 /*@@
    @routine    BndRadiativeDirGN
    @date       Sat Jan 20 2001
    @author     Gabrielle Allen
-   @desc 
-      Apply radiative BC's by group name in given direction
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Apply radiative BC's by group name in given direction
+   @enddesc
+   @calls      BndRadiativeDirGI
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil_size
+   @vdesc      stencil size in this direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to apply BC
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        gn_to
+   @vdesc      name of group to apply BC to
+   @vtype      char []
+   @vio        in
+   @endvar
+   @var        gn_from
+   @vdesc      name of group to apply BC from
+   @vtype      char []
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine BndRadiativeDirGI <BR>
+               -1 if invalid group names are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeDirGN(cGH *GH, 
-		      int stencil_size,
-		      int dir,
-		      CCTK_REAL var0, 
-		      CCTK_REAL speed, 
-		      const char *gn, 
-		      const char *gn_p)
+int BndRadiativeDirGN (cGH *GH,
+                       int stencil_size,
+                       int dir,
+                       CCTK_REAL var0,
+                       CCTK_REAL speed,
+                       char gn_to[],
+                       char gn_from[])
 {
-  int retval;
-  int gi   = CCTK_GroupIndex(gn);
-  int gi_p = CCTK_GroupIndex(gn_p);
-  
-  if ((gi>-1)&&(gi_p>-1))
-    retval = BndRadiativeDirGI(GH, stencil_size, dir, var0, speed, gi, gi_p);
-  else 
+  int gi_to, gi_from, retval;
+
+
+  gi_to   = CCTK_GroupIndex (gn_to);
+  gi_from = CCTK_GroupIndex (gn_from);
+  if (gi_to >= 0 && gi_from >= 0)
   {
-    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-	      "BndRadiativeDirGN: Grid variable groups %s or %s not found",
-	      gn, gn_p);
+    retval = BndRadiativeDirGI (GH, stencil_size, dir, var0, speed,
+                                gi_to, gi_from);
+  }
+  else
+  {
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid group names '%s' and/or '%s' in BndRadiativeDirGN",
+                gn_to, gn_from);
     retval = -1;
   }
 
-  return retval;
+  return (retval);
 }
- 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeDirGN)
-     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
-      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
+
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirGN)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS)
 {
-  TWO_FORTSTRINGS_CREATE(gn, gn_p)
-  *ierr = BndRadiativeDirGN(GH, *stencil_size, *dir, *var0, *speed, gn, gn_p);
-  free(gn);
-  free(gn_p);
-  return;
+  TWO_FORTSTRINGS_CREATE (gn_to, gn_from)
+  *ierr = BndRadiativeDirGN (GH, *stencil_size, *dir, *var0, *speed,
+                             gn_to, gn_from);
+  free (gn_to);
+  free (gn_from);
 }
 
 
 /*@@
-   @routine    BndRadiativeGI
+   @routine    BndRadiativeGN
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
-   @desc 
-      Aply radiative BC's by group name   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Aply radiative BC's by group name
+   @enddesc
+   @calls      BndRadiativeGI
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil
+   @vdesc      stencil width
+   @vtype      int [ dimension of group ]
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        gn_to
+   @vdesc      name of group to apply BC to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        gn_from
+   @vdesc      name of group to apply BC from
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine BndRadiativeGI <BR>
+               -1 if invalid group names are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeGN(cGH *GH, 
-		   int *stencil_array, 
-		   CCTK_REAL var0, 
-		   CCTK_REAL speed, 
-		   const char *gn, 
-		   const char *gn_p)
+int BndRadiativeGN (cGH *GH,
+                    int stencil[],
+                    CCTK_REAL var0,
+                    CCTK_REAL speed,
+                    char gn_to[],
+                    char gn_from[])
 {
-  int retval;
-  int gi   = CCTK_GroupIndex(gn);
-  int gi_p = CCTK_GroupIndex(gn_p);
-  
-  if ((gi>-1)&&(gi_p>-1))
-    retval = BndRadiativeGI(GH, stencil_array, var0, speed, gi, gi_p);
-  else 
+  int gi_to, gi_from, retval;
+
+
+  gi_to   = CCTK_GroupIndex (gn_to);
+  gi_from = CCTK_GroupIndex (gn_from);
+  if (gi_to >= 0 && gi_from >= 0)
   {
-    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-	      "BndRadiativeGN: Grid variable groups %s or %s not found",
-	      gn, gn_p);
+    retval = BndRadiativeGI (GH, stencil, var0, speed, gi_to, gi_from);
+  }
+  else
+  {
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid group names '%s' and/or '%s' in BndRadiativeGN",
+                gn_to, gn_from);
     retval = -1;
   }
 
-  return retval;
+  return (retval);
 }
- 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeGN)
-     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
-      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
+
+void CCTK_FCALL CCTK_FNAME (BndRadiativeGN)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS)
 {
-  TWO_FORTSTRINGS_CREATE(gn, gn_p)
-  *ierr = BndRadiativeGN(GH, stencil_array, *var0, *speed, gn, gn_p);
-  free(gn);
-  free(gn_p);
-  return;
+  TWO_FORTSTRINGS_CREATE (gn_to, gn_from)
+  *ierr = BndRadiativeGN (GH, stencil, *var0, *speed, gn_to, gn_from);
+  free (gn_to);
+  free (gn_from);
 }
 
 
 /* ===================================================================== */
 
-
 /*@@
    @routine    BndRadiativeDirVI
    @date       Sat Jan 20 2001
    @author     Gabrielle Allen
-   @desc 
-      Apply radiative BC's by variable index in given direction
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Apply radiative BC's by variable index in given direction
+   @enddesc
+   @calls      ApplyBndRadiative
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil_size
+   @vdesc      stencil size in this direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to apply BC
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        vi_to
+   @vdesc      index of variable to apply BC to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        vi_from
+   @vdesc      index of variable to apply BC from
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine ApplyBndRadiative <BR>
+               -1 if invalid variable indices are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeDirVI(cGH *GH, 
-		      int stencil_size, 
-		      int dir,
-		      CCTK_REAL var0, 
-		      CCTK_REAL speed, 
-		      int vi, 
-		      int vi_p) 
+int BndRadiativeDirVI (cGH *GH,
+                       int stencil_size,
+                       int dir,
+                       CCTK_REAL var0,
+                       CCTK_REAL speed,
+                       int vi_to,
+                       int vi_from)
 {
-  return ApplyBndRadiative(GH,stencil_size,NULL,var0,speed,vi,vi_p,1,dir);
+  int retval, num_vars;
+
+
+  num_vars = CCTK_NumVars ();
+  if (vi_to >= 0 && vi_to < num_vars && vi_from >= 0 && vi_from < num_vars)
+  {
+    retval = ApplyBndRadiative (GH, stencil_size, NULL, dir, var0, speed,
+                                vi_to, vi_from, 1);
+  }
+  else
+  {
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid variable indices %d and/or %d in BndRadiativeDirVI",
+                vi_to, vi_from);
+    retval = -1;
+  }
+
+  return (retval);
 }
 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeDirVI)
-     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
-      CCTK_REAL *speed, int *vi, int *vi_p)
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirVI)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *vi_to,
+                            const int *vi_from)
 {
-  *ierr = BndRadiativeDirVI(GH, *stencil_size, *dir, *var0, 
-			    *speed, *vi, *vi_p);
-  return;
+  *ierr = BndRadiativeDirVI (GH, *stencil_size, *dir, *var0, *speed,
+                             *vi_to, *vi_from);
 }
 
 
-
 /*@@
    @routine    BndRadiativeVI
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
-   @desc 
-      Apply radiative BC's by variable index 
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Apply radiative BC's by variable index
+   @enddesc
+   @calls      ApplyBndRadiative
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil
+   @vdesc      stencil width
+   @vtype      int [ dimension of variable ]
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        vi_to
+   @vdesc      index of variable to apply BC to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        vi_from
+   @vdesc      index of variable to apply BC from
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine ApplyBndRadiative <BR>
+               -1 if invalid variable indices are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeVI(cGH *GH, 
-                    int *stencil_array, 
-                    CCTK_REAL var0, 
-                    CCTK_REAL speed, 
-                    int vi, 
-                    int vi_p) 
+int BndRadiativeVI (cGH *GH,
+                    int stencil[],
+                    CCTK_REAL var0,
+                    CCTK_REAL speed,
+                    int vi_to,
+                    int vi_from)
 {
-  return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,vi,vi_p,1,0);
-}
-
-void CCTK_FCALL CCTK_FNAME(BndRadiativeVI)
-     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
-      CCTK_REAL *speed, int *vi, int *vi_p)
-{
-  *ierr = BndRadiativeVI(GH, stencil_array, *var0, *speed, *vi, *vi_p);
-  return;
-}
+  int retval, num_vars;
 
 
-/* ======================================================================= */
-
-
-/*@@
-   @routine    BndRadiativeDirGI
-   @date       Sat Jan 20 2001
-   @author     Gabrielle Allen
-   @desc 
-      Apply radiative BC's by variable name in given direction
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
-@@*/
-
-int BndRadiativeDirVN(cGH *GH, 
-                   int stencil_size,
-		   int dir,
-                   CCTK_REAL var0, 
-                   CCTK_REAL speed, 
-                   const char *vn, 
-                   const char *vn_p) 
-{
-  int retval;
-  int vi   = CCTK_VarIndex(vn);
-  int vi_p = CCTK_VarIndex(vn_p);
-  
-  if (vi>=0 && vi_p>=0)
+  num_vars = CCTK_NumVars ();
+  if (vi_to >= 0 && vi_to < num_vars && vi_from >= 0 && vi_from < num_vars)
   {
-    retval = BndRadiativeDirVI(GH, stencil_size, dir, var0, speed, vi, vi_p);
+    retval = ApplyBndRadiative (GH, -1, stencil, 0, var0, speed,
+                                vi_to, vi_from, 1);
   }
   else
   {
-    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-	      "BndRadiativeDirVN: Grid variable %s or %s not found",
-	      vn, vn_p);
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid variable indices %d and/or %d in BndRadiativeVI",
+                vi_to, vi_from);
     retval = -1;
   }
 
-  return retval;
-
+  return (retval);
 }
 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeDirVN)
-     (int *ierr, cGH *GH, int *stencil_size, int *dir, CCTK_REAL *var0, 
-      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
+void CCTK_FCALL CCTK_FNAME (BndRadiativeVI)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            const int *vi_to,
+                            const int *vi_from)
 {
-  TWO_FORTSTRINGS_CREATE(vn,vn_p)
-  *ierr = BndRadiativeDirVN(GH, *stencil_size, *dir, *var0, *speed, vn, vn_p);
-  free(vn);
-  free(vn_p);
-  return;
-}  
+  *ierr = BndRadiativeVI (GH, stencil, *var0, *speed, *vi_to, *vi_from);
+}
 
 
-/*@@
-   @routine    BndRadiativeGI
-   @date       Tue Jul 18 18:04:07 2000
-   @author     Gerd Lanfermann
-   @desc 
-      Apply radiative BC's by variable name 
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
+/* ======================================================================= */
 
+/*@@
+   @routine    BndRadiativeDirVN
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc
+               Apply radiative BC's by variable name in given direction
+   @enddesc
+   @calls      BndRadiativeDirVI
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil_size
+   @vdesc      stencil size in this direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to apply BC
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        vn_to
+   @vdesc      name of variable to apply BC to
+   @vtype      char []
+   @vio        in
+   @endvar
+   @var        vn_from
+   @vdesc      name of variable to apply BC from
+   @vtype      char []
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine BndRadiativeDirVI <BR>
+               -1 if invalid variable names are given
+   @endreturndesc
 @@*/
-
-int BndRadiativeVN(cGH *GH, 
-                   int *stencil_array,
-                   CCTK_REAL var0, 
-                   CCTK_REAL speed, 
-                   const char *vn, 
-                   const char *vn_p) 
+int BndRadiativeDirVN (cGH *GH,
+                       int stencil_size,
+                       int dir,
+                       CCTK_REAL var0,
+                       CCTK_REAL speed,
+                       char vn_to[],
+                       char vn_from[])
 {
-  int retval;
-  int vi   = CCTK_VarIndex(vn);
-  int vi_p = CCTK_VarIndex(vn_p);
-  
-  if (vi>=0 && vi_p>=0)
+  int vi_to, vi_from, num_vars, retval;
+
+
+  vi_to    = CCTK_VarIndex (vn_to);
+  vi_from  = CCTK_VarIndex (vn_from);
+  num_vars = CCTK_NumVars ();
+
+  if (vi_to >= 0 && vi_to < num_vars && vi_from >= 0 && vi_from < num_vars)
   {
-    retval = BndRadiativeVI(GH, stencil_array, var0, speed, vi, vi_p);
+    retval = BndRadiativeDirVI (GH, stencil_size, dir, var0, speed,
+                                vi_to, vi_from);
   }
   else
   {
-    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-	      "BndRadiativeVN: Grid variable %s or %s not found",
-	      vn, vn_p);
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid variable names '%s' and/or '%s' in BndRadiativeDirVN",
+                vn_to, vn_from);
     retval = -1;
   }
 
-  return retval;
-
+  return (retval);
 }
 
-void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
-     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
-      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
+void CCTK_FCALL CCTK_FNAME (BndRadiativeDirVN)
+                           (int *ierr,
+                            cGH *GH,
+                            const int *stencil_size,
+                            const int *dir,
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS)
 {
-  TWO_FORTSTRINGS_CREATE(vn,vn_p)
-  *ierr = BndRadiativeVN(GH, stencil_array, *var0, *speed, vn, vn_p);
-  free(vn);
-  free(vn_p);
-  return;
-}  
-
+  TWO_FORTSTRINGS_CREATE (vn_to, vn_from)
+  *ierr = BndRadiativeDirVN (GH, *stencil_size, *dir, *var0, *speed,
+                             vn_to, vn_from);
+  free (vn_to);
+  free (vn_from);
+}
 
 
-/********************************************************************
- *********************     Local Routines   *************************
- ********************************************************************/
 /*@@
-   @routine    ApplyBndRadiative
+   @routine    BndRadiativeVN
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
-   @desc 
-      Apply radiative BC's - internal routine
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
+   @desc
+               Apply radiative BC's by variable name
+   @enddesc
+   @calls      BndRadiativeVI
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil
+   @vdesc      stencil width
+   @vtype      int [ dimension of variable ]
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        vn_to
+   @vdesc      name of variable to apply BC to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        vn_from
+   @vdesc      name of variable to apply BC from
+   @vtype      int
+   @vio        in
+   @endvar
+
+   @returntype int
+   @returndesc
+               return code of @seeroutine BndRadiativeVI <BR>
+               -1 if invalid variable names are given
+   @endreturndesc
 @@*/
-
-static int ApplyBndRadiative(cGH *GH, 
-			     int stencil,
-			     int *stencil_array,
-			     CCTK_REAL var0, 
-			     CCTK_REAL speed, 
-			     int var_current, 
-			     int var_previous,
-			     int num_vars,
-			     int dir) 
+int BndRadiativeVN (cGH *GH,
+                    int stencil[],
+                    CCTK_REAL var0,
+                    CCTK_REAL speed,
+                    char vn_to[],
+                    char vn_from[])
 {
-  int xi=-1;
-  int yi=-1;
-  int zi=-1;
-  int ri=-1;    
-  int count;
-  int gi;
-  int vi1, vi2;
-  int symmetry_handle;
-  int idim, dim;
-  int time, time_p;
-  int berr,ierr;
-  int dirloop1,dirloop2;
-  int doBC;
-  int retval;
-  int *sw;
-  int *dstag;
-  int *lssh;
-  CCTK_REAL dx[3], dt;
-  SymmetryGHex *sGHex;
+  int vi_to, vi_from, num_vars, retval;
 
-#ifdef DEBUG_BOUNDARY
-  printf("Input arguments for ApplyBndRadiative:\n");
-  printf("GH = %x\n",GH);
-  printf("stencil = %d\n",stencil);
-  printf("stencil_array = %x\n",stencil_array);
-  printf("var0 = %f\n",var0);
-  printf("speed = %f\n",speed);
-  printf("var_current = %d\n",var_current);
-  printf("var_previous = %d\n",var_previous);
-  printf("num_vars = %d\n",num_vars); 
-  printf("dir = %d\n",dir);
-  printf("-----------------------------------\n");
-#endif
-  
-  retval = 0;
 
-  /* Get group dimension */
-  dim = CCTK_GroupDimFromVarI(var_current);
+  vi_to    = CCTK_VarIndex (vn_to);
+  vi_from  = CCTK_VarIndex (vn_from);
+  num_vars = CCTK_NumVars ();
 
-  if (dim<3 || dim>3)
+  if (vi_to >= 0 && vi_to < num_vars && vi_from >= 0 && vi_from < num_vars)
   {
-    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
-	       "ApplyBndRadiative: Invalid grid function dimension %d",
-	       dim);
-    retval=-1;
+    retval = BndRadiativeVI (GH, stencil, var0, speed, vi_to, vi_from);
   }
   else
   {
-    /* See if we have a symmetry array */
-    symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
-    if (symmetry_handle < 0)
-    {
-      sGHex = NULL;
-    }
-    else
-    {
-      sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
-    }
-
-    /* Set up stencil width array if needed */
-    if (stencil_array && dir == 0)
-    {
-      sw = stencil_array;
-    }
-    else if (dir==0)
-    {
-      CCTK_WARN(0,"ApplyBndRadiative: Direction 0 invalid for directional BC");
-    }
-    else if (dir>dim||-dir>dim)
-    {
-      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
-		 "ApplyBndRadiative: Direction %d greater than dimension %d",
-		 dir,dim);
-    }
-    else
-    {
-
-      int i;
-      /* Should be single direction */
-      sw  = (int *)malloc(dim*sizeof(int));
-      if (sw)
-      {
-	for (i=0;i<dim;i++)
-	{
-	  sw[i] = 0;
-	}
-	if (dir > 0)
-	{
-	  sw[dir-1] = stencil;
-	}
-	else
-	{
-	  sw[-dir-1] = stencil;
-	}
-      }
-      else
-      {
-	CCTK_WARN(0,"ApplyBndRadiative: Malloc sw failed");
-      }
-    }
-
-    /* Radiative boundaries need the underlying Cartesian coordinates */
-    switch (dim) 
-    {
-    case 1: 
-      xi = CCTK_CoordIndex(-1,"x","cart1d");
-      break;
-    case 2:
-      xi = CCTK_CoordIndex(-1,"x","cart2d");
-      yi = CCTK_CoordIndex(-1,"y","cart2d");
-      break;
-    case 3:
-      xi = CCTK_CoordIndex(-1,"x","cart3d");
-      yi = CCTK_CoordIndex(-1,"y","cart3d");
-      zi = CCTK_CoordIndex(-1,"z","cart3d");
-      ri = CCTK_CoordIndex(-1,"r","spher3d");
-      break;
-    default:
-      CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
-		 "ApplyBndRadiative: Radiative boundary only for dim=3: grid variable '%s'",
-		 CCTK_VarName(var_current));
-      break;
-    }
-
-    /* Allocate temporary arrays */
-    dstag = (int *)malloc(dim*sizeof(int));
-    lssh  = (int *)malloc(dim*sizeof(int));
-
-    /* get the directional staggering of the group */
-    gi   = CCTK_GroupIndexFromVarI(var_current);
-    ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
-
-    
-    /* Use next time level, if available */
-    time = CCTK_NumTimeLevelsFromVarI(var_current) - 1;
-    if (time < 0)
-    {
-      time = 0;
-    }
-    /* Use current time level, if available */
-    time_p = CCTK_NumTimeLevelsFromVarI(var_previous) - 2;
-    if (time_p < 0)
-    {
-      time_p = 0;
-    }
-
-    /* Treat all boundaries or just one? */
-    if (dir>0)
-    {
-      dirloop1 = 2*(dir-1)+1;
-      dirloop2 = dirloop1+1;
-    }
-    else if (dir<0)
-    {
-      dirloop1 = 2*(-dir-1);
-      dirloop2 = dirloop1+1;
-    }
-    else
-    {
-      dirloop1 = 0;
-      dirloop2 = 2*dim;
-    }
-  
-    for (vi1=var_current; vi1<var_current+num_vars; vi1++) 
-    {
-      
-      vi2 = (vi1-var_current+var_previous);
-
-      /* Apply condition if:
-	 + boundary is not a symmetry boundary 
-	 (no symmetry or unset(=unsed))
-	 + boundary is a physical boundary
-	 + have enough grid points            */
-
-      for (idim=0;idim<dim;idim++)
-      {
-	/* FIXME: THIS LINE FAILS WAVE TESTSUITES */
-        /*lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];*/
-	lssh[idim] = GH->cctk_lsh[idim];
-      }
-
-      for (count=dirloop1;count<dirloop2;count++)
-      {
-	if (sGHex)
-	{
-	  doBC = (
-		  ( (sGHex->GFSym[vi1][count]==GFSYM_NOSYM) ||
-		    (sGHex->GFSym[vi1][count]==GFSYM_UNSET) ) &&
-		  GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
-	    count : -1;
-	}
-	else
-	{
-	  doBC  = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
-	    count : -1 ;
-	}
-	
-	if (doBC >=0)
-	{
-	  switch (dim)
-	  {
-	  case 1: 
-	    berr = -1; 
-	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
-		       "ApplyBndRadiative: Radiative boundary conditions" 
-		       "not implemented for dimension = 1, grid variable %s",
-		       CCTK_VarName(vi1));
-	    break;
-	  case 2: 
-	    berr = -1; 
-	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
-		       "ApplyBndRadiative: Radiative boundary conditions"
-		       "not implemented for dimension = 2, grid variable %s",
-		       CCTK_VarName(vi1));
-	    break;	  
-	  case 3:   
-	    /* According to the Cactus spec, the true dx and dt values for a
-	       grid are calculated as follows: */
-	    dx[0] = GH->cctk_delta_space[0] / GH->cctk_levfac[0];
-	    dx[1] = GH->cctk_delta_space[1] / GH->cctk_levfac[1];
-	    dx[2] = GH->cctk_delta_space[2] / GH->cctk_levfac[2];
-	    dt    = GH->cctk_delta_time;
-    
-	    berr = BndApplyRadiative3Di
-	      (GH,
-	       doBC,		      /* apply BC on which face */
-	       lssh, 	              /* local shape, respects staggering */
-	       sw,	              /* stencil width array */
-	       dx,	              /* dx,dy,dz */
-	       dt,		      /* dt */
-	       GH->data[vi1][time  ], /* start of data array GF[] */
-	       GH->data[vi2][time_p], /* start of prev data array */
-	       GH->data[xi][0],       /* pointer to the X field */
-	       GH->data[yi][0],       /* pointer to the Y field */
-	       GH->data[zi][0],       /* pointer to the Z field */
-	       GH->data[ri][0],       /* pointer to the R field */
-	       var0,		      /* asymptotic limit */
-	       speed);		      /* wave speed */
-	    break;
-	  default : 
-	    berr=-1; 
-	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
-		       "ApplyBndRadiative: No Radiative BC for dim>3, grid variable %s",
-		       CCTK_VarName(vi1));
-	  }
-	}
-      }
-      retval=(berr<0)?-1:0;
-    }
-    free(dstag);
-    free(lssh);
-
-    if (!stencil_array)
-    {
-      free(sw);
-    }
+    CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "Invalid variable names '%s' and/or '%s' in BndRadiativeVN",
+                vn_to, vn_from);
+    retval = -1;
   }
 
-  return retval;
+  return (retval);
 }
 
-/*@@
-   @routine BndApplyRadiative3Di
-   @date    May 2000
-   @author  Miguel Alcubierre
-   @desc
-            Implements radiative boundary condition.
-
-            The radiative boundary condition that is implemented is:
-
-
-               f  =  f0  +  u(r - v*t) / r  +  h(r + v*t) / r
-
-
-            That is, I assume outgoing radial waves with a 1/r
-            fall off, and the correct asymptotic value f0, plus
-            I include the possibility of incoming waves as well
-            (these incoming waves should be modeled somehow).
-
-            The condition above leads to the differential equation:
+void CCTK_FCALL CCTK_FNAME (BndRadiativeVN)
+                           (int *ierr,
+                            cGH *GH,
+                            int stencil[],
+                            const CCTK_REAL *var0,
+                            const CCTK_REAL *speed,
+                            TWO_FORTSTRINGS_ARGS)
+{
+  TWO_FORTSTRINGS_CREATE (vn_to, vn_from)
+  *ierr = BndRadiativeVN (GH, stencil, *var0, *speed, vn_to, vn_from);
+  free (vn_to);
+  free (vn_from);
+}
 
-                                                                  
-            (x / r) d f  +  v d f  + v x (f - f0) / r^2  =  v x H / r^2
-              i      t         i        i                      i
 
-            where x_i is the normal direction to the given boundaries,
-            and H = 2 dh(s)/ds.
+/********************************************************************
+ *********************     Local Routines   *************************
+ ********************************************************************/
 
-            So at a given boundary I only worry about derivatives in
-            the normal direction.  Notice that u(r-v*t) has dissapeared,
-            but we still do not know the value of H.
+/* shortcut for multiplying a variable with itself */
+#define SQR(a) ((a) * (a))
 
-            To get H what I do is the following:  I evaluate the
-            expression one point in from the boundary and solve for H
-            there.  We now need a way of extrapolation H to the boundary.
-            For this I assume that H falls off as a power law:
+/* the maximum dimension we can deal with */
+#define MAXDIM 3
 
 
-            H = k/r**n   =>  d H  =  - n H/r
-                              i
+/*@@
+   @routine    LOWER_RADIATIVE_BOUNDARY_3D
+   @date       Mon 9 Apr 2001
+   @author     Thomas Radke
+   @desc
+               Macro to apply radiative BC to a lower bound of a 3D variable
+   @enddesc
 
-            The value of n is is defined by the parameter "radpower".
-            If this parameter is negative, H is forced to be zero (this
-            corresponds to pure outgoing waves and is the default).
+   @var        istart, jstart, kstart
+   @vdesc      start index for the x,y,z direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        xyz
+   @vdesc      coordinate arrays for x, y, and z
+   @vtype      CCTK_REAL [ MAXDIM ]
+   @vio        in
+   @endvar
+   @var        radius
+   @vdesc      radius coordinate array
+   @vtype      CCTK_REAL []
+   @vio        in
+   @endvar
+   @var        offset
+   @vdesc      offset to the next element in this direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var_to
+   @vdesc      target variable array
+   @vtype      cctk_type []
+   @vio        in
+   @endvar
+   @var        var_from
+   @vdesc      source variable array
+   @vtype      cctk_type []
+   @vio        in
+   @endvar
+   @var        cctk_type
+   @vdesc      CCTK datatypes of the source and target variable
+   @vtype      <cctk_type>
+   @vio        in
+   @endvar
+@@*/
+#define LOWER_RADIATIVE_BOUNDARY_3D(istart, jstart, kstart,                   \
+                                    xyz,                                      \
+                                    radius,                                   \
+                                    offset,                                   \
+                                    var_to,                                   \
+                                    var_from,                                 \
+                                    cctk_type)                                \
+{                                                                             \
+  int _i, _j, _k;                                                             \
+                                                                              \
+                                                                              \
+  for (_k = kstart-1; _k >= 0; _k--)                                          \
+  {                                                                           \
+    for (_j = jstart-1; _j >= 0; _j--)                                        \
+    {                                                                         \
+      int _0 = 0*offset,                                                      \
+          _1 = 1*offset,                                                      \
+          _2 = 2*offset;                                                      \
+      int _index = CCTK_GFINDEX3D (GH, istart-1, _j, _k);                     \
+      CCTK_REAL *_radius = radius   + _index,                                 \
+                *_xyz    = xyz      + _index;                                 \
+      cctk_type *_to     = (cctk_type *) (var_to)   + _index,                 \
+                *_from   = (cctk_type *) (var_from) + _index;                 \
+                                                                              \
+                                                                              \
+      for (_i = istart-1; _i >= 0; _i--)                                      \
+      {                                                                       \
+        CCTK_REAL _radius0_inv = 1/_radius[_0],                               \
+                  _radius1_inv = 1/_radius[_1];                               \
+                                                                              \
+                                                                              \
+        if (radpower > 0)                                                     \
+        {                                                                     \
+          CCTK_REAL H;                                                        \
+                                                                              \
+                                                                              \
+          H  = 0.25 * radpower * dxyz[0]                                      \
+               * (_xyz[_0]*SQR (_radius0_inv) + _xyz[_1]*SQR (_radius1_inv)); \
+          H  = (1 + H) / (1 - H);                                             \
+          H *= dtv * (0.25*(_to[_1] + _to[_2] + _from[_1] + _from[_2]) - var0)\
+               + 0.5*(  _radius[_1] * (_to[_1] - _from[_1])                   \
+                      + _radius[_2] * (_to[_2] - _from[_2]))                  \
+               + 0.25*(_to[_2] - _to[_1] + _from[_2] - _from[_1]) * rho[0]    \
+                     *(  SQR (_radius[_1]) / _xyz[_1]                         \
+                       + SQR (_radius[_2]) / _xyz[_2]);                       \
+          dtvvar0H = dtvvar0 + H;                                             \
+        }                                                                     \
+                                                                              \
+        _to[_0] = (cctk_type) (                                               \
+                   (dtvvar0H * (  _xyz[_0] * SQR (_radius0_inv)               \
+                                + _xyz[_1] * SQR (_radius1_inv))              \
+                    - _to[_1]   * ( rho[0] + _xyz[_1]*_radius1_inv            \
+                                                     *(1 + dtvh*_radius1_inv))\
+                    + _from[_0] * ( rho[0] + _xyz[_0]*_radius0_inv            \
+                                                     *(1 - dtvh*_radius0_inv))\
+                    - _from[_1] * ( rho[0] - _xyz[_1]*_radius1_inv            \
+                                                     *(1 - dtvh*_radius1_inv))\
+                   )                                                          \
+                   /              (-rho[0] + _xyz[_0]*_radius0_inv            \
+                                                     *(1 + dtvh*_radius0_inv))\
+                  );                                                          \
+        _radius--; _xyz--; _to--; _from--;                                    \
+      }                                                                       \
+    }                                                                         \
+  }                                                                           \
+}
 
-            The behaviour I have observed is the following:  Using H=0
-            is very stable, but has a very bad initial transient. Taking
-            n to be 0 or positive improves the initial behaviour considerably,
-            but introduces a drift that can kill the evolution at very late
-            times.  Empirically, the best value I have found is n=2, for
-            which the initial behaviour is very nice, and the late time drift
-            is quite small.
 
-            Another problem with this condition is that it does not
-            use the physical characteristic speed, but rather it assumes
-            a wave speed of v, so the boundaries should be out in
-            the region where the characteristic speed is constant.
-            Notice that this speed does not have to be 1.  For gauge
-            quantities {alpha,phi,trK} we can have a different asymptotic
-            speed, which is why the value of v is passed as a parameter.
+/*@@
+   @routine    UPPER_RADIATIVE_BOUNDARY_3D
+   @date       Mon 9 Apr 2001
+   @author     Thomas Radke
+   @desc
+               Macro to apply radiative BC to an upper bound of a 3D variable
+   @enddesc
 
+   @var        istart, jstart, kstart
+   @vdesc      start index for the x,y,z direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        xyz
+   @vdesc      coordinate arrays for x, y, and z
+   @vtype      CCTK_REAL [ MAXDIM ]
+   @vio        in
+   @endvar
+   @var        radius
+   @vdesc      radius coordinate array
+   @vtype      CCTK_REAL []
+   @vio        in
+   @endvar
+   @var        offset
+   @vdesc      offset to the next element in this direction
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var_to
+   @vdesc      target variable array
+   @vtype      cctk_type []
+   @vio        in
+   @endvar
+   @var        var_from
+   @vdesc      source variable array
+   @vtype      cctk_type []
+   @vio        in
+   @endvar
+   @var        cctk_type
+   @vdesc      CCTK datatypes of the source and target variable
+   @vtype      <cctk_type>
+   @vio        in
+   @endvar
+@@*/
+#define UPPER_RADIATIVE_BOUNDARY_3D(istart, jstart, kstart,                   \
+                                    xyz,                                      \
+                                    radius,                                   \
+                                    offset,                                   \
+                                    var_to,                                   \
+                                    var_from,                                 \
+                                    cctk_type)                                \
+{                                                                             \
+  int _i, _j, _k;                                                             \
+                                                                              \
+                                                                              \
+  for (_k = kstart; _k < GH->cctk_lsh[2]; _k++)                               \
+  {                                                                           \
+    for (_j = jstart; _j < GH->cctk_lsh[1]; _j++)                             \
+    {                                                                         \
+      int _0 = -0*offset,                                                      \
+          _1 = -1*offset,                                                      \
+          _2 = -2*offset;                                                      \
+      int _index = CCTK_GFINDEX3D (GH, istart, _j, _k);                       \
+      CCTK_REAL *_radius = radius   + _index,                                 \
+                *_xyz    = xyz      + _index;                                 \
+      cctk_type *_to     = (cctk_type *) (var_to)   + _index,                 \
+                *_from   = (cctk_type *) (var_from) + _index;                 \
+                                                                              \
+                                                                              \
+      for (_i = istart; _i < GH->cctk_lsh[0]; _i++)                           \
+      {                                                                       \
+        CCTK_REAL _radius0_inv = 1/_radius[_0],                               \
+                  _radius1_inv = 1/_radius[_1];                               \
+                                                                              \
+                                                                              \
+        if (radpower > 0)                                                     \
+        {                                                                     \
+          CCTK_REAL H;                                                        \
+                                                                              \
+                                                                              \
+          H  = 0.25 * radpower * dxyz[0]                                      \
+               * (_xyz[_0]*SQR (_radius0_inv) + _xyz[_1]*SQR (_radius1_inv)); \
+          H  = (1 - H) / (1 + H);                                             \
+          H *= dtv * (0.25*(_to[_1] + _to[_2] + _from[_1] + _from[_2]) - var0)\
+               + 0.5 *(  _radius[_1] * (_to[_1] - _from[_1])                  \
+                       + _radius[_2] * (_to[_2] - _from[_2]))                 \
+               + 0.25*(_to[_1] - _to[_2] + _from[_1] - _from[_2]) * rho[0]    \
+                     * (  SQR (_radius[_1]) / _xyz[_1]                        \
+                        + SQR (_radius[_2]) / _xyz[_2]);                      \
+          dtvvar0H = dtvvar0 + H;                                             \
+        }                                                                     \
+                                                                              \
+        _to[_0] = (cctk_type) (                                               \
+                   (dtvvar0H * (  _xyz[_0] * (SQR (_radius0_inv))             \
+                                + _xyz[_1] * (SQR (_radius1_inv)))            \
+                    + _to[_1]   * ( rho[0] - _xyz[_1]*_radius1_inv            \
+                                                     *(1 + dtvh*_radius1_inv))\
+                    + _from[_0] * (-rho[0] + _xyz[_0]*_radius0_inv            \
+                                                     *(1 - dtvh*_radius0_inv))\
+                    + _from[_1] * ( rho[0] + _xyz[_1]*_radius1_inv            \
+                                                     *(1 - dtvh*_radius1_inv))\
+                   )                                                          \
+                   /              ( rho[0] + _xyz[_0]*_radius0_inv            \
+                                                     *(1 + dtvh*_radius0_inv))\
+                  );                                                          \
+        _radius++; _xyz++; _to++; _from++;                                    \
+      }                                                                       \
+    }                                                                         \
+  }                                                                           \
+}
 
-	    Some remarks on the C version for Fortran programmers:
 
-	    This is 1:1 translation of the F code, which used the very
-            nice looking F array assignments x(2,:,:), etc. In C this
-            needs to be written as loops, but to my surprise it doesn't
-            look too bad actually. F statments like x(1,:,:) become x(xgp0)
-            (the integer index "xgp0" meaning the zeroth gridpoint from
-            the boundary), x(2,:,:) --> x(xgp1), where:
+/*@@
+   @routine    RADIATIVE_BOUNDARY
+   @date       Mon 9 Apr 2001
+   @author     Thomas Radke
+   @desc
+               Macro to apply radiative BC to a variable
+               Currently it is limited to 3D variables only.
+   @enddesc
+   @calls      LOWER_RADIATIVE_BOUNDARY_3D
+               UPPER_RADIATIVE_BOUNDARY_3D
+
+   @var        lsh
+   @vdesc      local shape of the variable
+   @vtype      int [ dim ]
+   @vio        in
+   @endvar
+   @var        doBC
+   @vdesc      flags telling whether to apply BC in a given direction or not
+   @vtype      int [ 2*dim ]
+   @vio        in
+   @endvar
+   @var        stencil
+   @vdesc      stencils in every direction
+   @vtype      int [ 2*dim ]
+   @vio        in
+   @endvar
+   @var        coords
+   @vdesc      coordinate arrays for x, y, z, and the radius
+   @vtype      CCTK_REAL [ MAXDIM+1 ]
+   @vio        in
+   @endvar
+   @var        offset
+   @vdesc      offsets to the next element in each dimension
+   @vtype      int [ dim ]
+   @vio        in
+   @endvar
+   @var        var_to
+   @vdesc      target variable array
+   @vtype      cctk_type []
+   @vio        in
+   @endvar
+   @var        var_from
+   @vdesc      source variable array
+   @vtype      cctk_type []
+   @vio        in
+   @endvar
+   @var        dim
+   @vdesc      dimension of the source and target variable
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        cctk_type
+   @vdesc      CCTK datatypes of the source and target variable
+   @vtype      <cctk_type>
+   @vio        in
+   @endvar
+@@*/
+#define RADIATIVE_BOUNDARY(lsh,                                               \
+                           doBC,                                              \
+                           stencil,                                           \
+                           coords,                                            \
+                           offset,                                            \
+                           var_to,                                            \
+                           var_from,                                          \
+                           dim,                                               \
+                           cctk_type)                                         \
+{                                                                             \
+  /* check the dimensionality */                                              \
+  if (dim != 3)                                                               \
+  {                                                                           \
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,                      \
+                "ApplyBndRadiative: variable dimension of %d not supported",  \
+                dim);                                                         \
+    return (-5);                                                              \
+  }                                                                           \
+                                                                              \
+  /* Lower x-bound */                                                         \
+  if (doBC[0])                                                                \
+  {                                                                           \
+    LOWER_RADIATIVE_BOUNDARY_3D (stencil[0], lsh[1], lsh[2],                  \
+                                 coords[0], coords[MAXDIM], offset[0],        \
+                                 var_to, var_from, cctk_type);                \
+  }                                                                           \
+                                                                              \
+  /* Upper x-bound */                                                         \
+  if (doBC[1])                                                                \
+  {                                                                           \
+    UPPER_RADIATIVE_BOUNDARY_3D (lsh[0] - stencil[0], 0, 0,                   \
+                                 coords[0], coords[MAXDIM], offset[0],        \
+                                 var_to, var_from, cctk_type);                \
+  }                                                                           \
+                                                                              \
+  /* Lower y-bound */                                                         \
+  if (doBC[2])                                                                \
+  {                                                                           \
+    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], stencil[1], lsh[2],                  \
+                                 coords[1], coords[MAXDIM], offset[1],        \
+                                 var_to, var_from, cctk_type);                \
+  }                                                                           \
+                                                                              \
+  /* Upper y-bound */                                                         \
+  if (doBC[3])                                                                \
+  {                                                                           \
+    UPPER_RADIATIVE_BOUNDARY_3D (0, lsh[1] - stencil[1], 0,                   \
+                                 coords[1], coords[MAXDIM], offset[1],        \
+                                 var_to, var_from, cctk_type);                \
+  }                                                                           \
+                                                                              \
+  /* Lower z-bound */                                                         \
+  if (doBC[4])                                                                \
+  {                                                                           \
+    LOWER_RADIATIVE_BOUNDARY_3D (lsh[0], lsh[1], stencil[2],                  \
+                                 coords[2], coords[MAXDIM], offset[2],        \
+                                 var_to, var_from, cctk_type);                \
+  }                                                                           \
+                                                                              \
+  /* Upper z-bound */                                                         \
+  if (doBC[5])                                                                \
+  {                                                                           \
+    UPPER_RADIATIVE_BOUNDARY_3D (0, 0, lsh[2] - stencil[2],                   \
+                                 coords[2], coords[MAXDIM], offset[2],        \
+                                 var_to, var_from, cctk_type);                \
+  }                                                                           \
+}
 
-	    xgp0 = CCTK_GFINDEX3D(GH,0,j,k)
-	    xgp1 = CCTK_GFINDEX3D(GH,1,j,k)
 
-            and i,j are looping of the grid sizes.
+/*@@
+   @routine    ApplyBndRadiative
+   @date       Tue Jul 18 18:04:07 2000
+   @author     Gerd Lanfermann
+   @desc
+               Apply radiation boundary conditions to a group of grid functions
+               given by their indices
+               This routine is called by the various BndRadiativeXXX wrappers.
 
+               Although it is currently limited to handle 3D variables only
+               it can easily be extended for other dimensions
+               by adapting the appropriate macros.
    @enddesc
-   @calls 
-   @calledby
+
+   @var        GH
+   @vdesc      Pointer to CCTK grid hierarchy
+   @vtype      cGH *
+   @vio        in
+   @endvar
+   @var        stencil_dir
+   @vdesc      stencil width in direction dir
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        stencil_alldirs
+   @vdesc      stencil widths for all directions
+   @vtype      int [ dimension of variable(s) ]
+   @vio        in
+   @endvar
+   @var        dir
+   @vdesc      direction to copy boundaries (0 for copying all directions)
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        var0
+   @vdesc      asymptotic value of function at infinity
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        speed
+   @vdesc      wave speed
+   @vtype      CCTK_REAL
+   @vio        in
+   @endvar
+   @var        first_var_to
+   @vdesc      index of first variable to copy boundaries to
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        first_var_from
+   @vdesc      index of first variable to copy boundaries from
+   @vtype      int
+   @vio        in
+   @endvar
+   @var        num_vars
+   @vdesc      number of variables
+   @vtype      int
+   @vio        in
+   @endvar
+               CCTK_VarTypeI
+               CCTK_GroupDimFromVarI
+               CCTK_NumTimeLevelsFromVarI
+               RADIATIVE_BOUNDARY
    @history
-      Cactus 4.0 by Gerd Lanfermann
+   @hdate      Mon 9 Apr 2001
+   @hauthor    Thomas Radke
+   @hdesc      Merged separate routines for 1D, 2D, and 3D
+               into a single generic routine
    @endhistory
-@@*/
 
-static int BndApplyRadiative3Di(cGH *GH,
-				int doBC, 
-				int *lssh,
-				int *sw,  
-				CCTK_REAL *dxyz, 
-				CCTK_REAL dt, 
-				CCTK_REAL *var_n, 
-				CCTK_REAL *var_p,
-				CCTK_REAL *x, 
-				CCTK_REAL *y, 
-				CCTK_REAL *z, 
-				CCTK_REAL *r,
-				CCTK_REAL var0, 
-				CCTK_REAL speed) 
+   @returntype int
+   @returndesc
+                0 for success
+               -1 if abs(direction) is greater than variables' dimension
+               -2 if variable dimension is not supported
+               -3 if NULL pointer passed as stencil array
+               -4 if variable type is not supported
+               -5 if variable dimension is other than 3D
+   @endreturndesc
+@@*/
+static int ApplyBndRadiative (cGH *GH,
+                              int stencil_dir,
+                              int stencil_alldirs[],
+                              int dir,
+                              CCTK_REAL var0,
+                              CCTK_REAL speed,
+                              int first_var_to,
+                              int first_var_from,
+                              int num_vars)
 {
-
   DECLARE_CCTK_PARAMETERS
+  int i, gdim;
+  int var_to, var_from;
+  int timelvl_to, timelvl_from;
+  SymmetryGHex *sGHex;
+  char coord_system_name[10];
+  CCTK_REAL dxyz[MAXDIM], rho[MAXDIM], *coords[MAXDIM+1];
+  int doBC[2*MAXDIM], stencil[MAXDIM], offset[MAXDIM];
+  CCTK_REAL dtv, dtvh, dtvvar0, dtvvar0H;
 
-  int i,j,k;
-
-  CCTK_REAL dtv,dtvh,dtvvar0;
-  CCTK_REAL dx,dy,dz;
-  CCTK_REAL rhox,rhoy,rhoz;
-  CCTK_REAL H,fac,aux;
-  CCTK_REAL half,one;
-
-  /* Linear grid point index, used to access the 0th/1st gridpoint away 
-     from the BC in the  calculation. 
-     In Fortran, 
-       x(1,:,:) --> x[xgp0],
-       x(2,:,:) --> x[xgp1], 
-     where xgp0/1 are calculated with CCTK_GFINDEX3D */
-
-  int xgp0,xgp1,xgp2,ygp0,ygp1,ygp2,zgp0,zgp1,zgp2;
-
-  CCTK_REAL u0,u1,u2;
-  CCTK_REAL ui1,ui2;
-  CCTK_REAL r0,r1,r2;
-  CCTK_REAL ri0,ri1;
-
-  /* Grid sizes */
-
-  int lssh0, lssh1, lssh2;
-
-  /* stencil widths */
-
-  int sw0, sw1, sw2;
-
-  /* Grid parameters. */
-
-  dx = dxyz[0];
-  dy = dxyz[1];
-  dz = dxyz[2];
-
-  /* Numbers */
-
-  half = 0.5;
-  one  = 1.0;
-
-  /* Store the grid shape in local variables to help the optimiser */
-
-  lssh0 = lssh[0];
-  lssh1 = lssh[1];
-  lssh2 = lssh[2];
-
-  /* Ditto with the stencil widths */
-
-  sw0 = sw[0];
-  sw1 = sw[1];
-  sw2 = sw[2];
-
-  /* Find Courant parameters. */
-
-  dtv  = speed*dt;
-  dtvh = half*dtv;
 
-  dtvvar0 = dtv*var0;
+  /* check the direction parameter */
+  if (abs (dir) > MAXDIM)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "ApplyBndRadiative: direction %d is greater than maximum "
+                "dimension %d", dir, MAXDIM);
+    return (-1);
+  }
 
-  rhox = dtv/dx;
-  rhoy = dtv/dy;
-  rhoz = dtv/dz;
+  /* get the dimensionality */
+  gdim = CCTK_GroupDimFromVarI (first_var_to);
 
-  /* Extrapolation power */
+  /* check the dimensionality */
+  if (gdim > MAXDIM)
+  {
+    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "ApplyBndRadiative: variable dimension of %d not supported",
+                gdim);
+    return (-2);
+  }
 
-  if (radpower<0) 
+  /* set up stencil width array */
+  if (dir)
   {
-     fac = 0;
-  } 
-  else 
+    stencil[abs (dir) - 1] = stencil_dir;
+  }
+  else if (stencil_alldirs)
   {
-     fac = radpower;
+    memcpy (stencil, stencil_alldirs, gdim * sizeof (int));
   }
-
-  /* Lower x-bound: x(2,:,:) --> x[xgp2] */
-
-  if (doBC == 0) 
+  else
   {
-#ifdef DEBUG_BOUNDARY
-    printf("Applying radiative boundary (lower x)\n");
-#endif
-    for (k=0;k<lssh2;k++) 
-    {
-      for (j=0;j<lssh1;j++) 
-      {
-	for (i=sw0-1;i>=0;i--) 
-	{
-	  xgp0 = CCTK_GFINDEX3D(GH,i  ,j,k);
-	  xgp1 = CCTK_GFINDEX3D(GH,i+1,j,k);
-	  xgp2 = CCTK_GFINDEX3D(GH,i+2,j,k);
-
-          u0 = x[xgp0];
-          u1 = x[xgp1];
-          u2 = x[xgp2];
-
-          r0 = r[xgp0];
-          r1 = r[xgp1];
-          r2 = r[xgp2];
-
-          ui1 = 1.0/u1;
-          ui2 = 1.0/u2;
-
-          ri0 = 1.0/r0;
-          ri1 = 1.0/r1;
-
-          if (fac==0) 
-	  {
-
-             H = 0;
-
-          } 
-	  else 
-	  {
-
-             H = 0;
-
-             H = dtv*(-var0 + 0.25
-               *(var_n[xgp1] + var_n[xgp2]
-               + var_p[xgp1] + var_p[xgp2]))
-               +(r1*(var_n[xgp1] - var_p[xgp1])
-               + r2*(var_n[xgp2] - var_p[xgp2]))*half
-               + 0.25*rhox*(SQR(r1)*ui1 + SQR(r2)*ui2)
-               *(var_n[xgp2] - var_n[xgp1]
-               + var_p[xgp2] - var_p[xgp1]);
-
-             aux = 0.25*radpower*dx*(u0*SQR(ri0) + u1*SQR(ri1));
-
-             H = H*(one + aux)/(one - aux);
-
-          }
-
-	  var_n[xgp0] =
-	    ((dtvvar0 + H)*(u0*SQR(ri0) + u1*SQR(ri1))
-	     - var_n[xgp1]*( rhox + u1*ri1*(one + dtvh*ri1))
-	     + var_p[xgp0]*( rhox + u0*ri0*(one - dtvh*ri0))
-	     - var_p[xgp1]*( rhox - u1*ri1*(one - dtvh*ri1)))
-	    /              (-rhox + u0*ri0*(one + dtvh*ri0));
-
-	}
-      }
-    }
+    CCTK_WARN (1, "ApplyBndRadiative: NULL pointer passed "
+                  "for stencil width array");
+    return (-3);
   }
 
-  /* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */
-
-  if (doBC == 1)
+  /* Use next time level, if available */
+  timelvl_to = CCTK_NumTimeLevelsFromVarI (first_var_to) - 1;
+  if (timelvl_to < 0)
   {
-#ifdef DEBUG_BOUNDARY
-    printf("Applying radiative boundary (upper x)\n");
-#endif
-    for (k=0;k<lssh2;k++) 
-    {
-      for (j=0;j<lssh1;j++) 
-      {
-	for (i=lssh0-sw0;i<lssh0;i++) 
-	{
-
-	  xgp0 = CCTK_GFINDEX3D(GH,i  ,j,k);
-	  xgp1 = CCTK_GFINDEX3D(GH,i-1,j,k);
-	  xgp2 = CCTK_GFINDEX3D(GH,i-2,j,k);
-
-          u0 = x[xgp0];
-          u1 = x[xgp1];
-          u2 = x[xgp2];
-
-          r0 = r[xgp0];
-          r1 = r[xgp1];
-          r2 = r[xgp2];
-
-          ui1 = 1.0/u1;
-          ui2 = 1.0/u2;
-
-          ri0 = 1.0/r0;
-          ri1 = 1.0/r1;
-
-          if (fac==0) 
-	  {
-
-             H = 0;
-
-          } 
-	  else 
-	  {
-
-             H = 0;
-
-             H = dtv*(-var0 + 0.25
-               *(var_n[xgp1] + var_n[xgp2]
-               + var_p[xgp1] + var_p[xgp2]))
-               +(r1*(var_n[xgp1] - var_p[xgp1])
-               + r2*(var_n[xgp2] - var_p[xgp2]))*half
-               + 0.25*rhox*(SQR(r1)*ui1 + SQR(r2)*ui2)
-               *(var_n[xgp1] - var_n[xgp2]
-               + var_p[xgp1] - var_p[xgp2]);
-
-             aux = 0.25*radpower*dx*(u0*SQR(ri0) + u1*SQR(ri1));
-
-             H = H*(one - aux)/(one + aux);
-
-          }
-
-	  var_n[xgp0] =
-	    ((dtvvar0 + H)*(u0*(SQR(ri0)) + u1*(SQR(ri1)))
-	     + var_n[xgp1]*( rhox - u1*ri1*(one + dtvh*ri1))
-	     + var_p[xgp0]*(-rhox + u0*ri0*(one - dtvh*ri0))
-	     + var_p[xgp1]*( rhox + u1*ri1*(one - dtvh*ri1)))
-	    /              ( rhox + u0*ri0*(one + dtvh*ri0));
-
-	}
-      }
-    }
+    timelvl_to = 0;
   }
-
-  /* Lower y-bound */
-
-  if (doBC == 2) 
+  /* Use current time level, if available */
+  timelvl_from = CCTK_NumTimeLevelsFromVarI (first_var_from) - 2;
+  if (timelvl_from < 0)
   {
-#ifdef DEBUG_BOUNDARY
-    printf("Applying radiative boundary (lower y)\n");
-#endif
-    for (k=0;k<lssh2;k++) 
-    {
-      for (i=0;i<lssh0;i++) 
-      {
-	for (j=sw1-1;j>=0;j--) 
-        {
-	  ygp0 = CCTK_GFINDEX3D(GH,i,j  ,k);
-	  ygp1 = CCTK_GFINDEX3D(GH,i,j+1,k);
-	  ygp2 = CCTK_GFINDEX3D(GH,i,j+2,k);
-
-          u0 = y[ygp0];
-          u1 = y[ygp1];
-          u2 = y[ygp2];
-
-          r0 = r[ygp0];
-          r1 = r[ygp1];
-          r2 = r[ygp2];
-
-          ui1 = 1.0/u1;
-          ui2 = 1.0/u2;
-
-          ri0 = 1.0/r0;
-          ri1 = 1.0/r1;
-
-          if (fac==0) 
-	  {
-
-             H = 0;
-
-          } 
-	  else 
-	  {
-
-             H = 0;
-
-             H = dtv*(-var0 + 0.25
-               *(var_n[ygp1] + var_n[ygp2]
-               + var_p[ygp1] + var_p[ygp2]))
-               +(r1*(var_n[ygp1] - var_p[ygp1])
-               + r2*(var_n[ygp2] - var_p[ygp2]))*half
-               + 0.25*rhoy*(SQR(r1)*ui1 + SQR(r2)*ui2)
-               *(var_n[ygp2] - var_n[ygp1]
-               + var_p[ygp2] - var_p[ygp1]);
+    timelvl_from = 0;
+  }
 
-             aux = 0.25*radpower*dy*(u0*SQR(ri0) + u1*SQR(ri1));
+  /* Find Courant parameters. */
+  dtv      = speed * GH->cctk_delta_time;
+  dtvh     = 0.5 * dtv;
+  dtvvar0  = dtv * var0;
+  dtvvar0H = dtvvar0;
 
-             H = H*(one + aux)/(one - aux);
+  sprintf (coord_system_name, "cart%dd", gdim);
+  for (i = 0; i < gdim; i++)
+  {
+    /* Radiative boundaries need the underlying Cartesian coordinates */
+    coords[i] = GH->data[CCTK_CoordIndex (i + 1, NULL, coord_system_name)][0];
 
-          }
+    /* According to the Cactus spec, the true delta_space values for a
+       grid are calculated as follows: */
+    dxyz[i] = GH->cctk_delta_space[i] / GH->cctk_levfac[i];
 
-	  var_n[ygp0] =
-	    ((dtvvar0 + H)*(u0*SQR(ri0) + u1*SQR(ri1))
-	     - var_n[ygp1]*( rhoy + u1*ri1*(one + dtvh*ri1))
-	     + var_p[ygp0]*( rhoy + u0*ri0*(one - dtvh*ri0))
-	     - var_p[ygp1]*( rhoy - u1*ri1*(one - dtvh*ri1)))
-	    /              (-rhoy + u0*ri0*(one + dtvh*ri0));
+    rho[i] = dtv / dxyz[i];
 
-	}
-      }
-    }
+    offset[i] = i == 0 ? 1 : offset[i-1] * GH->cctk_lsh[i-1];
   }
+  sprintf (coord_system_name, "spher%dd", gdim);
+  coords[MAXDIM] = GH->data[CCTK_CoordIndex (-1, "r", coord_system_name)][0];
 
-  /* Upper y bound */
+  /* see if we have a symmetry array */
+  sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry");
 
-  if (doBC == 3) 
+  /* now loop over all variables */
+  for (var_to = first_var_to, var_from = first_var_from;
+       var_to < first_var_to + num_vars;
+       var_to++, var_from++)
   {
-#ifdef DEBUG_BOUNDARY
-    printf("Applying radiative boundary (upper y)\n");
-#endif
-    for (k=0;k<lssh2;k++) 
+    /* Apply condition if:
+       + boundary is not a symmetry boundary
+         (no symmetry or unset(=unsed))
+       + boundary is a physical boundary
+       + have enough grid points
+    */
+    memset (doBC, 1, sizeof (doBC));
+    if (sGHex)
     {
-      for (i=0;i<lssh0;i++) 
+      for (i = 0; i < 2 * MAXDIM; i++)
       {
-	for (j=lssh1-sw1;j<lssh1;j++) 
-        {
-
-	  ygp0 = CCTK_GFINDEX3D(GH,i,j  ,k);
-	  ygp1 = CCTK_GFINDEX3D(GH,i,j-1,k);
-	  ygp2 = CCTK_GFINDEX3D(GH,i,j-2,k);
-
-          u0 = y[ygp0];
-          u1 = y[ygp1];
-          u2 = y[ygp2];
-
-          r0 = r[ygp0];
-          r1 = r[ygp1];
-          r2 = r[ygp2];
-
-          ui1 = 1.0/u1;
-          ui2 = 1.0/u2;
-
-          ri0 = 1.0/r0;
-          ri1 = 1.0/r1;
-
-          if (fac==0) 
-          {
-
-             H = 0;
-
-          } 
-	  else 
-	  {
-
-             H = 0;
-
-             H = dtv*(-var0 + 0.25
-               *(var_n[ygp1] + var_n[ygp2]
-               + var_p[ygp1] + var_p[ygp2]))
-               +(r1*(var_n[ygp1] - var_p[ygp1])
-               + r2*(var_n[ygp2] - var_p[ygp2]))*half
-               + 0.25*rhoy*(SQR(r1)*ui1 + SQR(r2)*ui2)
-               *(var_n[ygp1] - var_n[ygp2]
-               + var_p[ygp1] - var_p[ygp2]);
-
-             aux = 0.25*radpower*dy*(u0*SQR(ri0) + u1*SQR(ri1));
-
-             H = H*(one - aux)/(one + aux);
-
-          }
-
-	  var_n[ygp0] =
-	    ((dtvvar0 + H)*(u0*(SQR(ri0)) + u1*(SQR(ri1)))
-	     + var_n[ygp1]*( rhoy - u1*ri1*(one + dtvh*ri1))
-	     + var_p[ygp0]*(-rhoy + u0*ri0*(one - dtvh*ri0))
-	     + var_p[ygp1]*( rhoy + u1*ri1*(one - dtvh*ri1)))
-	    /              ( rhoy + u0*ri0*(one + dtvh*ri0));
-
-	}
+        doBC[i] = sGHex->GFSym[var_to][i] == GFSYM_NOSYM ||
+                  sGHex->GFSym[var_to][i] == GFSYM_UNSET;
       }
     }
-  }
-
-  /* Lower z-bound */
-
-  if (doBC==4) 
-  {
-#ifdef DEBUG_BOUNDARY
-    printf("Applying radiative boundary (lower z)\n");
-#endif
-    for (j=0;j<lssh1;j++) 
+    for (i = 0; i < MAXDIM; i++)
     {
-      for (i=0;i<lssh0;i++) 
+      doBC[i*2]   &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2];
+      doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1];
+      if (dir != 0)
       {
-	for (k=sw2-1;k>=0;k--) 
-	{
-
-	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
-	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k+1);
-	  zgp2 = CCTK_GFINDEX3D(GH,i,j,k+2);
-
-          u0 = z[zgp0];
-          u1 = z[zgp1];
-          u2 = z[zgp2];
-
-          r0 = r[zgp0];
-          r1 = r[zgp1];
-          r2 = r[zgp2];
-
-          ui1 = 1.0/u1;
-          ui2 = 1.0/u2;
-
-          ri0 = 1.0/r0;
-          ri1 = 1.0/r1;
-
-          if (fac==0) 
-	  {
-
-             H = 0;
-
-          } 
-	  else 
-	  {
-
-             H = 0;
-
-             H = dtv*(-var0 + 0.25
-               *(var_n[zgp1] + var_n[zgp2]
-               + var_p[zgp1] + var_p[zgp2]))
-               +(r1*(var_n[zgp1] - var_p[zgp1])
-               + r2*(var_n[zgp2] - var_p[zgp2]))*half
-               + 0.25*rhoz*(SQR(r1)*ui1 + SQR(r2)*ui2)
-               *(var_n[zgp2] - var_n[zgp1]
-               + var_p[zgp2] - var_p[zgp1]);
-
-             aux = 0.25*radpower*dz*(u0*SQR(ri0) + u1*SQR(ri1));
-
-             H = H*(one + aux)/(one - aux);
-
-          }
-
-	  var_n[zgp0] =
-	    ((dtvvar0 + H)*(u0*SQR(ri0) + u1*SQR(ri1))
-	     - var_n[zgp1]*( rhoz + u1*ri1*(one + dtvh*ri1))
-	     + var_p[zgp0]*( rhoz + u0*ri0*(one - dtvh*ri0))
-	     - var_p[zgp1]*( rhoz - u1*ri1*(one - dtvh*ri1)))
-	    /              (-rhoz + u0*ri0*(one + dtvh*ri0));
-
-	}
+        doBC[i*2]   &= (dir < 0 && (i + 1 == abs (dir)));
+        doBC[i*2+1] &= (dir > 0 && (i + 1 == abs (dir)));
       }
     }
-  }
 
-  /* Upper z-bound */
-
-  if (doBC ==  5) 
-  {
-#ifdef DEBUG_BOUNDARY
-    printf("Applying radiative boundary (upper z)\n");
-#endif
-    for (j=0;j<lssh1;j++) 
+    switch (CCTK_VarTypeI (var_to))
     {
-      for (i=0;i<lssh0;i++) 
-      {
-	for (k=lssh2-sw2;k<lssh2;k++) 
-	{
-
-	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
-	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k-1);
-	  zgp2 = CCTK_GFINDEX3D(GH,i,j,k-2);
-
-          u0 = z[zgp0];
-          u1 = z[zgp1];
-          u2 = z[zgp2];
-
-          r0 = r[zgp0];
-          r1 = r[zgp1];
-          r2 = r[zgp2];
-
-          ui1 = 1.0/u1;
-          ui2 = 1.0/u2;
-
-          ri0 = 1.0/r0;
-          ri1 = 1.0/r1;
-
-          if (fac==0) 
-	  {
-
-             H = 0;
-
-          } 
-	  else 
-	  {
-
-             H = 0;
+      case CCTK_VARIABLE_CHAR:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_CHAR);
+        break;
+
+      case CCTK_VARIABLE_INT:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_INT);
+        break;
+
+      case CCTK_VARIABLE_REAL:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_REAL);
+        break;
+
+#ifdef CCTK_VARIABLE_INT2
+      case CCTK_VARIABLE_INT2:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_INT2);
+        break;
+#endif
 
-             H = dtv*(-var0 + 0.25
-               *(var_n[zgp1] + var_n[zgp2]
-               + var_p[zgp1] + var_p[zgp2]))
-               +(r1*(var_n[zgp1] - var_p[zgp1])
-               + r2*(var_n[zgp2] - var_p[zgp2]))*half
-               + 0.25*rhoz*(SQR(r1)*ui1 + SQR(r2)*ui2)
-               *(var_n[zgp1] - var_n[zgp2]
-               + var_p[zgp1] - var_p[zgp2]);
+#ifdef CCTK_VARIABLE_INT4
+      case CCTK_VARIABLE_INT4:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_INT4);
+        break;
+#endif
 
-             aux = 0.25*radpower*dz*(u0*SQR(ri0) + u1*SQR(ri1));
+#ifdef CCTK_VARIABLE_INT8
+      case CCTK_VARIABLE_INT8:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_INT8);
+        break;
+#endif
 
-             H = H*(one - aux)/(one + aux);
+#ifdef CCTK_REAL4
+      case CCTK_VARIABLE_REAL4:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_REAL4);
+        break;
+#endif
 
-          }
+#ifdef CCTK_REAL8
+      case CCTK_VARIABLE_REAL8:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim, CCTK_REAL8);
+        break;
+#endif
 
-	  var_n[zgp0] =
-	    ((dtvvar0 + H)*(u0*(SQR(ri0)) + u1*(SQR(ri1)))
-	     + var_n[zgp1]*( rhoz - u1*ri1*(one + dtvh*ri1))
-	     + var_p[zgp0]*(-rhoz + u0*ri0*(one - dtvh*ri0))
-	     + var_p[zgp1]*( rhoz + u1*ri1*(one - dtvh*ri1)))
-	    /              ( rhoz + u0*ri0*(one + dtvh*ri0));
+#ifdef CCTK_REAL16
+      case CCTK_VARIABLE_REAL16:
+        RADIATIVE_BOUNDARY (GH->cctk_lsh, doBC, stencil, coords, offset,
+                            GH->data[var_to][timelvl_to],
+                            GH->data[var_from][timelvl_from], gdim,CCTK_REAL16);
+        break;
+#endif
 
-	}
-      }
+      default:
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Unsupported variable type %d for variable '%s'",
+                    CCTK_VarTypeI (var_to), CCTK_VarName (var_to));
+        return (-4);
     }
   }
-  return 0;
-}
-
-
 
+  return (0);
+}