Adding directional boundary conditions for everything apart from Robin.

See thorn documentation for details

Next to come, boundary conditions for grid functions which aren't CCTK_REALs.


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

6 files changed, 1517 insertions, 441 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index d501260..c1c2176 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -177,7 +177,7 @@ integer \> group\_index\\
 \end{tabbing}
 }
 
-\subsection*{Arguments}
+\subsubsection*{Arguments}
 \begin{Lentry}
 \item[{\tt ierr}] Return value, negative value indicates the
 boundary condition was not successfully applied
@@ -201,6 +201,8 @@ boundary value of phi {\tt phi(nx,j,k)}, on the positive x-boundary will
 be copied from {\tt phi(nx-1,j,k)}. 
 
 \subsubsection*{Calling from C:}
+
+{\bf All Coordinate Directions:}
 \begin{verbatim}
 int ierr = BndFlatVN(cGH *cctkGH, int *stencil_size, char *variable_name)
 int ierr = BndFlatGN(cGH *cctkGH, int *stencil_size, char *group_name)
@@ -208,30 +210,62 @@ int ierr = BndFlatVI(cGH *cctkGH, int *stencil_size, int variable_index)
 int ierr = BndFlatGI(cGH *cctkGH, int *stencil_size, int group_index)
 \end{verbatim}
 
+{\bf Individual Coordinate Directions:}
+\begin{verbatim}
+int ierr = BndFlatVN(cGH *cctkGH, int stencil, int dir, char *variable_name)
+int ierr = BndFlatGN(cGH *cctkGH, int stencil, int dir, char *group_name)
+int ierr = BndFlatVI(cGH *cctkGH, int stencil, int dir, int variable_index)
+int ierr = BndFlatGI(cGH *cctkGH, int stencil, int dir, int group_index)
+\end{verbatim}
+
 \subsubsection*{Calling from Fortran:}
 
+{\bf All Coordinate Directions:}
+\begin{verbatim}
+call BndFlatVN(ierr, cctkGH, stencil_array, variable_name)
+call BndFlatGN(ierr, cctkGH, stencil_array, group_name)
+call BndFlatVI(ierr, cctkGH, stencil_array, variable_index)
+call BndFlatGI(ierr, cctkGH, stencil_array, group_index)
+\end{verbatim}
+
+{\bf Individual Coordinate Directions:}
 \begin{verbatim}
-call BndFlatVN(ierr, cctkGH, stencil_size, variable_name)
-call BndFlatGN(ierr, cctkGH, stencil_size, group_name)
-call BndFlatVI(ierr, cctkGH, stencil_size, variable_index)
-call BndFlatGI(ierr, cctkGH, stencil_size, group_index)
+call BndFlatVN(ierr, cctkGH, stencil, dir, variable_name)
+call BndFlatGN(ierr, cctkGH, stencil, dir, group_name)
+call BndFlatVI(ierr, cctkGH, stencil, dir, variable_index)
+call BndFlatGI(ierr, cctkGH, stencil, dir, group_index)
 \end{verbatim}
 where
+{\tt
+\begin{tabbing}
+character*(*) \= variable\_name\=\kill
+integer \> ierr \\
+CCTK\_POINTER \> cctkGH\\
+integer \> dir\\
+integer \> stencil\\
+integer \> stencil\_array(dim)\\
+character*(*) \> variable\_name\\
+character*(*) \> group\_name\\
+integer \> variable\_index\\
+integer \> group\_index\\
+\end{tabbing}
+}
+
+\subsubsection*{Arguments}
 \begin{Lentry}
-\item[{\tt integer ierr}] return value, operation failed when return
-value {\em negative}
-\item[{\tt CCTK\_POINTER cctkGH}] grid hierarchy pointer
-\item[{\tt CCTK\_REAL var0}] scalar value to apply
-\item[{\tt integer stencil\_size(dim)}] array of size {\tt dim}
-(dimension of the gridfunction). To how many points from the outer
-boundary to apply the boundary condition. 
-\item[{\tt character*(*) variable\_name}] Name of the variable.
-\item[{\tt character*(*) group\_name}] Name of the group.
-\item[{\tt integer variable\_index}] Variable index.
-\item[{\tt integer group\_index}] Group index.
+\item[{\tt ierr}] Return value, negative value indicates the
+boundary condition was not successfully applied
+\item[{\tt cctkGH}] Grid hierarchy pointer
+\item[{\tt dir}] Coordinate direction in which to apply boundary condition
+\item[{\tt stencil\_size}] Array with dimension of the grid function, containing the stencil width to apply the boundary at
+\item[{\tt variable\_name}] Name of the variable
+\item[{\tt group\_name}] Name of the group
+\item[{\tt variable\_index}] Variable index
+\item[{\tt group\_index}] Group index
 \end{Lentry}
 
 
+
 \subsection{Radiation Boundary Condition}
 
 This is a two level scheme. Grid functions are given for the current time 
@@ -293,22 +327,42 @@ the region where the characteristic speed is constant.
 Notice that this speed does not have to be 1.  
 
 \subsubsection*{Calling from C:}
+
+{\bf All Coordinate Directions:}
 \begin{verbatim}
 int ierr = BndRadiativeVN(cGH *cctkGH, int *stencil_size, 
-                          CCTK_REAL speed, CCTK_REAL limit, 
+                          CCTK_REAL limit, CCTK_REAL speed, 
                           char *variable_name, char *variable_name_past)
 int ierr = BndRadiativeGN(cGH *cctkGH, int *stencil_size, 
-                          CCTK_REAL speed, CCTK_REAL limit,  
+                          CCTK_REAL limit, CCTK_REAL speed,  
                           char *group_name, char *group_name_past)
 int ierr = BndRadiativeVI(cGH *cctkGH, int *stencil_size, 
-                          CCTK_REAL speed, CCTK_REAL limit,  
+                          CCTK_REAL limit, CCTK_REAL speed,  
                           int variable_index, int variable_index_past)
 int ierr = BndRadiativeGI(cGH *cctkGH, int *stencil_size, 
-                          CCTK_REAL speed, CCTK_REAL limit,  
+                          CCTK_REAL limit, CCTK_REAL speed,  
+                          int group_index, int group_index_past)
+\end{verbatim}
+
+{\bf Individual Coordinate Directions:}
+\begin{verbatim}
+int ierr = BndRadiativeVN(cGH *cctkGH, int stencil, int dir,
+                          CCTK_REAL limit, CCTK_REAL speed, 
+                          char *variable_name, char *variable_name_past)
+int ierr = BndRadiativeGN(cGH *cctkGH, int *stencil, int dir, 
+                          CCTK_REAL limit, CCTK_REAL speed,  
+                          char *group_name, char *group_name_past)
+int ierr = BndRadiativeVI(cGH *cctkGH, int *stencil, int dir, 
+                          CCTK_REAL limit, CCTK_REAL speed,  
+                          int variable_index, int variable_index_past)
+int ierr = BndRadiativeGI(cGH *cctkGH, int *stencil, int dir, 
+                          CCTK_REAL limit, CCTK_REAL speed,  
                           int group_index, int group_index_past)
 \end{verbatim}
 
 \subsubsection*{Calling from Fortran:}
+
+{\bf All Coordinate Directions:}
 \begin{verbatim}
 call BndRadiativeVN(ierr, cctkGH, stencil_size, speed, limit, 
                     variable_name, variable_name_past)
@@ -319,6 +373,39 @@ call BndRadiativeVN(ierr, cctkGH, stencil_size, speed, limit,
 call BndRadiativeVN(ierr, cctkGH, stencil_size, speed, limit, 
                     group_index, group_index_past)
 \end{verbatim}
+
+{\bf Individual Coordinate Directions:}
+\begin{verbatim}
+call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, 
+                    variable_name, variable_name_past)
+call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, 
+                    group_name, group_name_past)
+call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, 
+                    variable_index, variable_index_past)
+call BndRadiativeVN(ierr, cctkGH, stencil, dir, speed, limit, 
+                    group_index, group_index_past)
+\end{verbatim}
+where
+{\tt
+\begin{tabbing}
+character*(*) \= variable\_name\=\kill
+integer \> ierr \\
+CCTK\_POINTER \> cctkGH\\
+integer \> dir\\
+integer \> stencil\\
+integer \> stencil\_array(dim)\\
+character*(*) \> variable\_name\\
+character*(*) \> group\_name\\
+integer \> variable\_index\\
+integer \> group\_index\\
+CCTK_REAL\>speed\\
+CCTK_REAL\>limit\\
+\end{tabbing}
+}
+
+
+
+
 where
 \begin{Lentry}
 \item[{\tt integer ierr}] return value, operation failed when return
@@ -329,7 +416,7 @@ value {\em negative}
 boundary to apply the boundary condition. 
 \item[{\tt CCTK\_REAL speed}] wave speed used for boundary condition ($v$).
 	
-\item[{\tt CCTK\_REAL limit}] radation constant
+\item[{\tt CCTK\_REAL limit}] asymptotic value of function at infinity
 
 \item[{\tt character*(*) variable\_name}] the name of the grid function 
 	to which the boundary condition will be applied
diff --git a/src/Boundary.h b/src/Boundary.h
index a07ebe9..acddeb5 100644
--- a/src/Boundary.h
+++ b/src/Boundary.h
@@ -57,6 +57,27 @@ int BndScalarVN(cGH *GH,
 
 
 /* Copying boundaries  */
+int BndCopyDirGI(cGH *GH, 
+		 int stencil_size, 
+		 int dir,
+		 int gi_to, 
+		 int gi_from);
+int BndCopyDirGN(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 const char *gn_to, 
+		 const char *gn_from);
+int BndCopyDirVI(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 int vi_to, 
+		 int vi_from);
+int BndCopyDirVN(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 const char *vn_to, 
+		 const char *vn_from);
+
 int BndCopyGI(cGH *GH, 
 	      int *stencil, 
 	      int gi_to, 
@@ -76,6 +97,39 @@ int BndCopyVN(cGH *GH,
 
 
 /* Radiative boundaries */
+int BndRadiativeDirGI(cGH *GH, 
+		      int stencil_size,
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL v0, 
+		      int gi, 
+		      int gi_p);
+
+int BndRadiativeDirGN(cGH *GH, 
+		      int stencil_size,
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL v0, 
+		      const char *gn, 
+		      const char *gn_p);
+
+int BndRadiativeDirVI(cGH *GH, 
+		      int stencil_size,
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL v0, 
+		      int vi, 
+		      int vi_p);
+
+int BndRadiativeDirVN(cGH *GH, 
+		      int stencil_size,
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL v0, 
+		      const char *vn, 
+		      const char *vn_p);
+
+
 int BndRadiativeGI(cGH *GH, 
 		   int *sw,
 		   CCTK_REAL var0, 
@@ -132,6 +186,23 @@ int BndRobinVN(cGH *GH,
 
 
 /* Flat boundaries      */
+int BndFlatDirGI(cGH *GH, 
+		 int stencil,
+		 int dir,
+		 int gi);
+int BndFlatDirGN(cGH *GH, 
+		 int stencil, 
+		 int dir,
+		 const char *gn);
+int BndFlatDirVI(cGH *GH, 
+		 int stencil,
+		 int dir,
+		 int vi);
+int BndFlatDirVN(cGH *GH, 
+		 int stencil,
+		 int dir,
+		 const char *vn);
+
 int BndFlatGI(cGH *GH, 
 	      int *stencil, 
 	      int gi);
diff --git a/src/CopyBoundary.c b/src/CopyBoundary.c
index 93ebcb3..117b742 100644
--- a/src/CopyBoundary.c
+++ b/src/CopyBoundary.c
@@ -25,34 +25,33 @@
 /* Internal routine prototypes */
 
 static int BndApplyCopy3Di(cGH *GH,
-			   int  gdim,
-			   int *doBC,
+			   int doBC,
 			   int *lssh,
 			   int *stencil_size,
 			   CCTK_REAL *var,
 			   CCTK_REAL *varp);
 
 static int BndApplyCopy2Di(cGH *GH,
-			   int  gdim,
-			   int *doBC,
+			   int doBC,
 			   int *lssh,
 			   int *stencil_size,
 			   CCTK_REAL *var,
 			   CCTK_REAL *varp);
 
 static int BndApplyCopy1Di(cGH *GH,
-			   int  gdim,
-			   int *doBC,
+			   int doBC,
 			   int *lssh,
 			   int *stencil_size,
 			   CCTK_REAL *var,
 			   CCTK_REAL *varp);
 
 static int ApplyBndCopy(cGH *GH, 
-			int *stencil, 
+			int stencil,
+			int *stencil_array, 
 			int first_var,
 			int second_var,
-			int num_vars);
+			int num_vars,
+			int dir);
 
 /* Local variables */
 
@@ -63,6 +62,43 @@ static int ApplyBndCopy(cGH *GH,
  ********************************************************************/
 
 /*@@
+   @routine    BndCopyDirVI
+   @date       Sat Mar 20
+   @author     Gabrielle Allen
+   @desc 
+   Apply copy boundary routines by var index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndCopyDirVI(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 int vi_to, 
+		 int vi_from) 
+{
+  return ApplyBndCopy(GH,stencil_size,NULL,vi_to,vi_from,1,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndCopyDirVI)
+     (int *ierr, 
+      cGH *GH, 
+      int *stencil_size, 
+      int *dir,
+      int *vi_to, 
+      int *vi_from)
+{
+  *ierr = BndCopyDirVI(GH, *stencil_size, *dir, *vi_to, *vi_from);
+}
+
+
+
+/*@@
    @routine    BndCopyVI
    @date       Thu Mar  2 11:02:10 2000
    @author     Gerd Lanfermann
@@ -82,21 +118,54 @@ int BndCopyVI(cGH *GH,
 	      int vi_to, 
 	      int vi_from) 
 {
-  return ApplyBndCopy(GH,stencil,vi_to,vi_from,1);
+  return ApplyBndCopy(GH,-1,stencil,vi_to,vi_from,1,0);
 }
 
-void CCTK_FCALL CCTK_FNAME(BndCopyVI)
+/* ====================================================== */
+
+ /*@@
+   @routine    BndCopyDirGI
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply copy boundaries by group index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndCopyDirGI(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 int gi_to, 
+		 int gi_from) 
+{
+  int first_vi_to,first_vi_from;
+  int numvars;
+
+  first_vi_to  = CCTK_FirstVarIndexI(gi_to);
+  first_vi_from= CCTK_FirstVarIndexI(gi_from);
+  numvars      = CCTK_NumVarsInGroupI(gi_to);
+
+  return ApplyBndCopy(GH,stencil_size,NULL,first_vi_to,first_vi_from,numvars,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndCopyDirGI)
      (int *ierr, 
       cGH *GH, 
-      int *stencil, 
-      int *vi_to, 
-      int *vi_from)
+      int *stencil_size,
+      int *dir,
+      int *gi_to, 
+      int *gi_from)
 {
-  *ierr = BndCopyVI(GH, stencil, *vi_to, *vi_from);
+  *ierr = BndCopyDirGI(GH, *stencil_size, *dir, *gi_to, *gi_from);
+  return;
 }
 
-
-
  /*@@
    @routine    BndCopyGI
    @date       Thu Mar  2 11:07:11 2000
@@ -124,7 +193,7 @@ int BndCopyGI(cGH *GH,
   first_vi_from= CCTK_FirstVarIndexI(gi_from);
   numvars      = CCTK_NumVarsInGroupI(gi_to);
 
-  return ApplyBndCopy(GH,stencil,first_vi_to,first_vi_from,numvars);
+  return ApplyBndCopy(GH,-1,stencil,first_vi_to,first_vi_from,numvars,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndCopyGI)
@@ -137,7 +206,62 @@ void CCTK_FCALL CCTK_FNAME(BndCopyGI)
   *ierr = BndCopyGI(GH, stencil, *gi_to, *gi_from);
 }
 
+/* ======================================================= */
 
+/*@@
+   @routine    BndCopyDirGN
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply copy boundary routines by group name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndCopyDirGN(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 const char *gn_to, 
+		 const char *gn_from)
+{
+  int retval;
+  int gi_to, gi_from;
+
+  gi_to   = CCTK_GroupIndex(gn_to);
+  gi_from = CCTK_GroupIndex(gn_from);
+  
+  if ((gi_to>-1)&&(gi_from>-1))
+    retval = BndCopyDirGI(GH, stencil_size, dir, gi_to, gi_from);
+  else 
+  {
+    CCTK_VWarn(2,__LINE__,__FILE__,"Boundary",
+	       "BndCopyDirGN: Cannot find group indices for %s, %s ",
+	       gn_to,gn_from);
+    retval = -1;
+  }
+
+  return retval;
+
+}
+
+void CCTK_FCALL CCTK_FNAME(BndCopyDirGN)
+     (int *ierr, 
+      cGH *GH, 
+      int *stencil_size,
+      int *dir,
+      TWO_FORTSTRINGS_ARGS)
+{
+  TWO_FORTSTRINGS_CREATE(gn_to, gn_from)
+  *ierr = BndCopyDirGN(GH, *stencil_size, *dir, gn_to, gn_from);
+  free(gn_to);
+  free(gn_from);
+  return;
+}
 
 /*@@
    @routine    BndCopyGN
@@ -189,9 +313,64 @@ void CCTK_FCALL CCTK_FNAME(BndCopyGN)
   *ierr = BndCopyGN(GH, stencil, gn_to, gn_from);
   free(gn_to);
   free(gn_from);
+  return;
 }
 
+/* ================================================================ */
+
+/*@@
+   @routine    BndCopyDirVN
+   @date       Sat Jan 20 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply copy boundary routines by variable name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
 
+int BndCopyDirVN(cGH *GH, 
+		 int stencil_size,
+		 int dir,
+		 const char *vn_to, 
+		 const char *vn_from)
+{
+  int vi_to, vi_from;
+  int retval;
+
+  vi_to   = CCTK_VarIndex(vn_to);
+  vi_from = CCTK_VarIndex(vn_from);
+  
+  if ((vi_to>-1)&&(vi_from>-1))
+    retval = BndCopyDirVI(GH, stencil_size, dir, vi_to, vi_from);
+  else 
+  {
+    CCTK_VWarn(2,__LINE__,__FILE__,"Boundary",
+	       "BndCopyDirVN: Cannot find variable indices for %s, %s",
+	       vn_to,vn_from);
+    retval = -1;
+  }
+
+  return retval;
+}
+
+void CCTK_FCALL CCTK_FNAME(BndCopyDirVN)
+     (int *ierr, 
+      cGH *GH, 
+      int *stencil_size, 
+      int *dir,
+      TWO_FORTSTRINGS_ARGS)
+{
+  TWO_FORTSTRINGS_CREATE(vn_to, vn_from)
+  *ierr = BndCopyDirVN(GH, *stencil_size, *dir, vn_to, vn_from);
+  free(vn_to);
+  free(vn_from);
+  return;
+}
 
 /*@@
    @routine    BndCopyVN
@@ -267,117 +446,214 @@ void CCTK_FCALL CCTK_FNAME(BndCopyVN)
 @@*/
 
 static int ApplyBndCopy(cGH *GH, 
-			int *stencil, 
+			int stencil,
+			int *stencil_array,
 			int first_var,
 			int second_var,
-			int num_vars)
+			int num_vars,
+			int dir)
 {
   int symmetry_handle;  /* handle for the optional symmetry structure */
   int vi, vi2, gi, dim;
   int idim;
+  int retval=0;
+  int dirloop1,dirloop2;
+  int doBC,count;
   int berr,ierr;
-  int *doBC, *dstag, *lssh;
-  SymmetryGHex *sGHex;
+  int *dstag;
+  int *lssh;
+  int *sw;
   int timelevel;        /* timelevel that condition applied on */
+  SymmetryGHex *sGHex;
+
+#ifdef DEBUG_BOUNDARY
+  printf("Input arguments for ApplyBndCopy:\n");
+  printf("GH = %x\n",GH);
+  printf("stencil = %d\n",stencil);
+  printf("stencil_array = %x\n",stencil_array);
+  printf("first_var = %d\n",first_var);
+  printf("second_var = %d\n",second_var);
+  printf("num_vars = %d\n",num_vars); 
+  printf("dir = %d\n",dir);
+  printf("-----------------------------------\n");
+#endif
+
+  /* Get group dimension */
+  dim = CCTK_GroupDimFromVarI(first_var);
 
-  /* See if we have a symmetry array */
-  symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
-  if (symmetry_handle < 0)
+  if (dim<1 || dim>3)
   {
-    sGHex = NULL;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "ApplyBndCopy: Invalid grid function dimension %d",
+	       dim);
+    retval=-1;
   }
   else
   {
-    sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
-  }
 
-  /* Get group dimension */
-  dim = CCTK_GroupDimFromVarI(first_var);
+    /* 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,"ApplyBndCopy: Direction 0 invalid for directional BC");
+    }
+    else if (dir>dim||-dir>dim)
+    {
+      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+		 "ApplyBndCopy: 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,"ApplyBndCopy: Malloc sw failed");
+      }
+    }
 
-  /* allocate arrays */
-  doBC   = (int *)malloc((2*dim)*sizeof(int)); 
-  dstag  = (int *)malloc(dim*sizeof(int));
-  lssh   = (int *)malloc(dim*sizeof(int));
+    /* allocate arrays */
+    dstag  = (int *)malloc(dim*sizeof(int));
+    lssh   = (int *)malloc(dim*sizeof(int));
 
-  /* get the directional staggering of the group */
-  gi   = CCTK_GroupIndexFromVarI(first_var);
-  ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
+    /* get the directional staggering of the group */
+    gi   = CCTK_GroupIndexFromVarI(first_var);
+    ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
 
-  /* get the current timelevel */
-  timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
-  /*  if (timelevel > 0)
-      timelevel--;*/
+    /* get the current timelevel */
+    timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
 
-  for (vi=first_var; vi<first_var+num_vars; vi++) 
-  {
-    vi2 = (vi-first_var+second_var);
-
-    /* Apply condition if:
-       + boundary is not a symmetry boundary 
-       (no symmetry or unset(=unsed))
-       + boundary is a physical boundary
-       + have enough grid points
-    */
-    if (sGHex)
+    /* Treat all boundaries or just one? */
+    if (dir>0)
     {
-      for (idim=0;idim<dim;idim++) 
-      {
-	doBC[idim*2]   = (((sGHex->GFSym[vi][idim*2]==GFSYM_NOSYM)||
-			   (sGHex->GFSym[vi][idim*2]==GFSYM_UNSET)) &&
-			  GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-	doBC[idim*2+1] = (((sGHex->GFSym[vi][idim*2+1]==GFSYM_NOSYM)||
-			   (sGHex->GFSym[vi][idim*2+1]==GFSYM_UNSET)) &&
-			  GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
-	lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
-      }
+      dirloop1 = 2*(dir-1)+1;
+      dirloop2 = dirloop1+1;
+    }
+    else if (dir<0)
+    {
+      dirloop1 = 2*(-dir-1);
+      dirloop2 = dirloop1+1;
     }
     else
     {
-      for (idim=0;idim<dim;idim++) 
+      dirloop1 = 0;
+      dirloop2 = 2*dim;
+    }
+  
+    for (vi=first_var; vi<first_var+num_vars; vi++) 
+    {
+      vi2 = (vi-first_var+second_var);
+
+      /* 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++)
       {
-	doBC[idim*2]   = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-	doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
-	lssh[idim]     = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
+	lssh[idim] = (int)GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
       }
-    }
+
+      for (count=dirloop1;count<dirloop2;count++)
+      {
+	if (sGHex)
+	{
+	  doBC = (
+		  ( (sGHex->GFSym[vi][count]==GFSYM_NOSYM) ||
+		    (sGHex->GFSym[vi][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 ;
+	}
 
 
-    switch (dim) 
-    {
-      case 1: berr = BndApplyCopy1Di(GH, 
-				     dim,
+	if (doBC >=0)
+	{
+	  switch (dim)
+	  {
+	  case 1: 
+	    berr = BndApplyCopy1Di(GH, 
 				     doBC, 
 				     lssh,
-				     stencil,
+				     sw,
 				     GH->data[vi][timelevel],
-				     GH->data[vi2][timelevel]); break;
-      case 2: berr = BndApplyCopy2Di(GH, 
-				     dim,
+				     GH->data[vi2][timelevel]); 
+	    break;
+	  case 2: 
+	    berr = BndApplyCopy2Di(GH, 
 				     doBC, 
 				     lssh,
-				     stencil,
+				     sw,
 				     GH->data[vi][timelevel],
-				     GH->data[vi2][timelevel]); break;
-      case 3: berr = BndApplyCopy3Di(GH, 
-				     dim,
+				     GH->data[vi2][timelevel]); 
+	    break;
+	  case 3: 
+	    berr = BndApplyCopy3Di(GH, 
 				     doBC,
 				     lssh,
-				     stencil,
+				     sw,
 				     GH->data[vi][timelevel],
-				     GH->data[vi2][timelevel]); break;     
-      default : berr = -1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", 
-		   "ApplyBndCopy: No BC for dim>3: grid variables '%s'/'%s'",
-		    CCTK_VarName(vi),CCTK_VarName(vi2));
+				     GH->data[vi2][timelevel]); 
+	    break;     
+	  default : 
+	    berr = -1; 
+	    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", 
+		       "No copy boundary for dim>3: grid variable '%s'",
+		       CCTK_VarName(vi));
+	  }
+	}
+      }
+
+      /*      retval = (berr>-1) ? 0 : -1;*/
+    }
+    
+    free(dstag);
+    free(lssh);
+
+    if (!stencil_array)
+    {
+      free(sw);
     }
-    berr = (berr>-1) ? 0 : -1;
   }
-  
-  free(dstag);
-  free(doBC);
-  free(lssh);
-  
-  return(ierr);
+
+  return retval;
 }
 
 
@@ -391,7 +667,7 @@ static int ApplyBndCopy(cGH *GH,
       for the 3d case
    @enddesc 
    @calls     
-   @calledby BndCopyGI,  BndCopyVN,  BndCopyGN,  BndCopyVI,  
+   @calledby
    @history 
  
    @endhistory 
@@ -399,8 +675,7 @@ static int ApplyBndCopy(cGH *GH,
 @@*/
 
 static int BndApplyCopy3Di(cGH *GH,
-			   int  gdim,
-			   int *doBC,
+			   int doBC,
 			   int *lssh,
 			   int *stencil_size,
 			   CCTK_REAL *var,
@@ -409,8 +684,11 @@ static int BndApplyCopy3Di(cGH *GH,
   int i,j,k,sw,lin;  
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (lower x)\n");
+#endif
     for (k=0;k<lssh[2];k++)
     {
       for (j=0;j<lssh[1];j++)
@@ -425,8 +703,11 @@ static int BndApplyCopy3Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (upper x)\n");
+#endif
     for (k=0;k<lssh[2];k++)
     {
       for (j=0;j<lssh[1];j++)
@@ -441,8 +722,11 @@ static int BndApplyCopy3Di(cGH *GH,
   }
 
   /* lower y */
-  if (doBC[2] == 1)
+  if (doBC == 2)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (lower y)\n");
+#endif
     for (k=0;k<lssh[2];k++)
     {
       for (i=0;i<lssh[0];i++)
@@ -457,8 +741,11 @@ static int BndApplyCopy3Di(cGH *GH,
   }
 
   /* upper y */
-  if (doBC[3] == 1)
+  if (doBC == 3)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (upper y)\n");
+#endif
     for (k=0;k<lssh[2];k++)
     {
       for (i=0;i<lssh[0];i++)
@@ -473,8 +760,11 @@ static int BndApplyCopy3Di(cGH *GH,
   }
 
   /* lower z */
-  if (doBC[4] == 1)
+  if (doBC == 4)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (lower z)\n");
+#endif
     for (j=0;j<lssh[1];j++)
     {
       for (i=0;i<lssh[0];i++)
@@ -489,8 +779,11 @@ static int BndApplyCopy3Di(cGH *GH,
   }
 
   /* upper z */  
-  if (doBC[5] == 1)
+  if (doBC == 5)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (upper z)\n");
+#endif
     for (j=0;j<lssh[1];j++)
     {
       for (i=0;i<lssh[0];i++)
@@ -526,8 +819,7 @@ static int BndApplyCopy3Di(cGH *GH,
 @@*/
 
 static int BndApplyCopy2Di(cGH *GH,
-			   int  gdim,
-			   int *doBC,
+			   int doBC,
 			   int *lssh,
 			   int *stencil_size,
 			   CCTK_REAL *var,
@@ -537,8 +829,11 @@ static int BndApplyCopy2Di(cGH *GH,
   int i,j,sw,lin;  
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (lower x)\n");
+#endif
     for (j=0;j<lssh[1];j++)
     {
       for (sw=0;sw<stencil_size[0];sw++)
@@ -550,8 +845,11 @@ static int BndApplyCopy2Di(cGH *GH,
   }
    
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (upper x)\n");
+#endif
     for (j=0;j<lssh[1];j++)
     {
       for (sw=0;sw<stencil_size[0];sw++)
@@ -563,8 +861,11 @@ static int BndApplyCopy2Di(cGH *GH,
   }
 
   /* lower y */
-  if (doBC[2] == 1)
+  if (doBC == 2)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (lower y)\n");
+#endif
     for (i=0;i<lssh[0];i++)
     {
       for (sw=0;sw<stencil_size[1];sw++)
@@ -576,8 +877,11 @@ static int BndApplyCopy2Di(cGH *GH,
   }
 
   /* upper y */
-  if (doBC[3] == 1)
+  if (doBC == 3)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (upper y)\n");
+#endif
     for (i=0;i<lssh[0];i++)
     {
       for (sw=0;sw<stencil_size[1];sw++)
@@ -610,8 +914,7 @@ static int BndApplyCopy2Di(cGH *GH,
 @@*/
 
 static int BndApplyCopy1Di(cGH *GH,
-			   int  gdim,
-			   int *doBC,
+			   int doBC,
 			   int *lssh,
 			   int *stencil_size,
 			   CCTK_REAL *var,
@@ -622,8 +925,11 @@ static int BndApplyCopy1Di(cGH *GH,
 
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (lower x)\n");
+#endif
     for (sw=0;sw<stencil_size[0];sw++)
     {
       lin=CCTK_GFINDEX1D(GH,sw);
@@ -632,8 +938,11 @@ static int BndApplyCopy1Di(cGH *GH,
   }
    
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying copy boundary (upper x)\n");
+#endif
     for (sw=0;sw<stencil_size[0];sw++)
     {
       lin=CCTK_GFINDEX1D(GH,lssh[0]-sw-1);
diff --git a/src/FlatBoundary.c b/src/FlatBoundary.c
index 2cc8780..616f2da 100644
--- a/src/FlatBoundary.c
+++ b/src/FlatBoundary.c
@@ -8,7 +8,7 @@
    @version $Header$
  @@*/
 
-/*#define DEBUG_BOUND*/
+/*#define DEBUG_BOUNDARY*/
 
 #include <stdio.h>
 #include <assert.h>
@@ -28,30 +28,29 @@
 /* Internal routine prototypes */
 
 static int BndApplyFlat3Di(cGH *GH, 
-			   int  gdim, 
-			   int *doBC, 
+			   int doBC, 
 			   int *lssh,
-			   int *stencil, 
+			   int *stencilarray, 
 			   CCTK_REAL *var);
 
 static int BndApplyFlat2Di(cGH *GH, 
-			   int  gdim, 
-			   int *doBC, 
+			   int doBC, 
 			   int *lssh,
-			   int *stencil, 
+			   int *stencil_array, 
 			   CCTK_REAL *var);
 
 static int BndApplyFlat1Di(cGH *GH, 
-			   int  gdim, 
-			   int *doBC, 
+			   int doBC, 
 			   int *lssh,
-			   int *stencil, 
+			   int *stencil_array, 
 			   CCTK_REAL *var);
 
 static int ApplyBndFlat(cGH *GH, 
-			int *stencil, 
+			int stencil,
+			int *stencil_array, 
 			int first_var,
-			int num_vars);
+			int num_vars,
+			int dir);
 
 /* Local variables */
 
@@ -62,6 +61,39 @@ static int ApplyBndFlat(cGH *GH,
  ********************************************************************/
 
 /*@@
+   @routine    BndFlatDirGI
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by group index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirGI(cGH *GH, int stencil, int dir, int gi) 
+{
+  int numvars, first_vi;
+
+  first_vi = CCTK_FirstVarIndexI(gi);
+  numvars  = CCTK_NumVarsInGroupI(gi);
+
+  return ApplyBndFlat(GH,stencil,NULL,first_vi,numvars,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndFlatDirGI)
+     (int *ierr, cGH *GH, int *stencil, int *dir, int *gi)
+{
+  *ierr = BndFlatDirGI(GH, *stencil, *dir, *gi);
+  return;
+}
+
+
+/*@@
    @routine    BndFlatGI
    @date       Thu Mar  2 11:11:40 2000
    @author     Gerd Lanfermann
@@ -83,7 +115,7 @@ int BndFlatGI(cGH *GH, int *stencil, int gi)
   first_vi = CCTK_FirstVarIndexI(gi);
   numvars  = CCTK_NumVarsInGroupI(gi);
 
-  return ApplyBndFlat(GH,stencil,first_vi,numvars);
+  return ApplyBndFlat(GH,-1,stencil,first_vi,numvars,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndFlatGI)
@@ -94,6 +126,51 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGI)
 }
 
 
+ /* ********************************************************************/
+
+
+/*@@
+   @routine    BndFlatDirGN
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by group name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirGN(cGH *GH, int stencil, int dir, const char *gn)
+{
+  int retval;
+  int gi;
+  gi   = CCTK_GroupIndex(gn);
+  
+  if (gi>-1)
+  {
+    retval=BndFlatDirGI(GH, stencil, dir, gi);
+  }
+  else 
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "BndFlatDirGN: Grid variable %s not found",gn);
+    retval=-1;
+  }
+  return retval;
+}
+ 
+void CCTK_FCALL CCTK_FNAME(BndFlatDirGN)
+     (int *ierr, cGH *GH, int *stencil, int *dir, ONE_FORTSTRING_ARG)
+{
+  ONE_FORTSTRING_CREATE(gn)
+  *ierr = BndFlatDirGN(GH, *stencil, *dir, gn);
+  free(gn);
+  return;
+}
 
 
 /*@@
@@ -140,6 +217,37 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGN)
 }
 
 
+ /* ********************************************************************/
+
+
+/*@@
+   @routine    BndFlatDirVI
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by variable index in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirVI(cGH *GH, int stencil, int dir, int vi) 
+{
+  return ApplyBndFlat(GH,stencil,NULL,vi,1,dir);
+}
+
+void CCTK_FCALL CCTK_FNAME(BndFlatDirVI)
+     (int *ierr, cGH *GH, int *stencil, int *dir, int *vi)
+{
+  *ierr = BndFlatDirVI(GH, *stencil, *dir, *vi);
+  return;
+}
+
+
 
 /*@@
    @routine    BndFlatVI
@@ -158,7 +266,7 @@ void CCTK_FCALL CCTK_FNAME(BndFlatGN)
 
 int BndFlatVI(cGH *GH, int *stencil, int vi) 
 {
-  return ApplyBndFlat(GH,stencil,vi,1);
+  return ApplyBndFlat(GH,-1,stencil,vi,1,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndFlatVI)
@@ -169,6 +277,54 @@ void CCTK_FCALL CCTK_FNAME(BndFlatVI)
 }
 
 
+ /* ********************************************************************/
+
+
+/*@@
+   @routine    BndFlatDirVN
+   @date       Sun Jan 21 2001
+   @author     Gabrielle Allen
+   @desc 
+   Apply flat boundary conditions by variable name in given direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+int BndFlatDirVN(cGH *GH, int stencil, int dir, const char *vn)
+{
+  int  vi;
+  int  retval;
+
+  vi   = CCTK_VarIndex(vn);
+  
+  if (vi>-1)
+  {
+    retval = BndFlatDirVI(GH, stencil, dir, vi);
+  }
+  else
+  {
+    retval = -1;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "BndFlatDirVN: Grid variable %s not found",vn);
+  }
+  return retval;
+}
+
+void CCTK_FCALL CCTK_FNAME(BndFlatDirVN)
+     (int *ierr, cGH *GH, int *stencil, int *dir, ONE_FORTSTRING_ARG)
+{  
+  ONE_FORTSTRING_CREATE(vn)
+  *ierr = BndFlatDirVN(GH, *stencil, *dir, vn);
+  free(vn);
+  return;
+}
+  
+
 
 /*@@
    @routine    BndFlatVN
@@ -236,118 +392,208 @@ void CCTK_FCALL CCTK_FNAME(BndFlatVN)
 @@*/
 
 int ApplyBndFlat(cGH *GH, 
-		 int *stencil, 
+		 int stencil,
+		 int *stencil_array, 
 		 int first_var,
-		 int num_vars)
+		 int num_vars,
+		 int dir)
 {
   int symmetry_handle;  /* handle for the optional symmetry structure */
   int vi, gi, dim;
   int idim; 
-  int berr=0,ierr=0; 
-  int *doBC, *dstag, *lssh;
+  int berr=0,ierr=0;
+  int retval=0;
+  int dirloop1,dirloop2;
+  int doBC;
+  int count;
+  int *dstag;
+  int *sw;
+  int *lssh;
   int timelevel;        /* timelevel that condition applied on */
   SymmetryGHex *sGHex;
 
-  /* See if we have a symmetry array */
-  symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
-  if (symmetry_handle < 0)
+#ifdef DEBUG_BOUNDARY
+  printf("Input arguments for ApplyBndFlat:\n");
+  printf("GH = %x\n",GH);
+  printf("stencil = %d\n",stencil);
+  printf("stencil_array = %x\n",stencil_array);
+  printf("first_var = %d\n",first_var);
+  printf("num_vars = %d\n",num_vars); 
+  printf("dir = %d\n",dir);
+  printf("-----------------------------------\n");
+#endif
+
+   /* Get group dimension */
+  dim = CCTK_GroupDimFromVarI(first_var);
+
+  if (dim<1 || dim>3)
   {
-    sGHex = NULL;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "ApplyBndScalar: Invalid grid function dimension %d",
+	       dim);
+    retval=-1;
   }
   else
   {
-    sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
-  }
+    /* 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];
+    }
 
-  /* Get group dimension */
-  dim = CCTK_GroupDimFromVarI(first_var);
+    /* Set up stencil width array if needed */
+    if (stencil_array && dir == 0)
+    {
+      sw = stencil_array;
+    }
+    else if (dir==0)
+    {
+      CCTK_WARN(0,"ApplyBndScalar: Direction 0 invalid for directional BC");
+    }
+    else if (dir>dim||-dir>dim)
+    {
+      CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+		 "ApplyBndScalar: 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,"ApplyBndScalar: Malloc sw failed");
+      }
+    }
 
-  /* allocate arrays */
-  doBC   = (int *)malloc((2*dim)*sizeof(int)); 
-  dstag  = (int *)malloc(dim*sizeof(int));
-  lssh   = (int *)malloc(dim*sizeof(int));
-
-  /* get the directional staggering of the group */
-  gi   = CCTK_GroupIndexFromVarI(first_var);
-  ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
-  if (ierr<0) CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-	"Cannot get stagger information for group '%s' (GF '$s')\n",
-	 CCTK_GroupName(gi), CCTK_VarName(first_var));
-  
+    /* allocate arrays */
+    dstag  = (int *)malloc(dim*sizeof(int));
+    lssh   = (int *)malloc(dim*sizeof(int));
 
-  /* get the current timelevel */
-  timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
-  /*  if (timelevel > 0)
-      timelevel--;*/
+    /* get the directional staggering of the group */
+    gi   = CCTK_GroupIndexFromVarI(first_var);
+    ierr = CCTK_GroupStaggerDirArrayGI(dstag, dim, gi);
 
-  for (vi=first_var; vi<first_var+num_vars; vi++) 
-  {
+    /* get the current timelevel */
+    timelevel = CCTK_NumTimeLevelsFromVarI (first_var) - 1;
 
-    /* Apply condition if:
-       + boundary is not a symmetry boundary 
-       (no symmetry or unset(=unsed))
-       + boundary is a physical boundary
-       + have enough grid points
-    */
-    if (sGHex)
+    /* Treat all boundaries or just one? */
+    if (dir>0)
     {
-      for (idim=0;idim<dim;idim++) 
-      {
-	doBC[idim*2]   = (((sGHex->GFSym[vi][idim*2]==GFSYM_NOSYM)||
-			   (sGHex->GFSym[vi][idim*2]==GFSYM_UNSET)) &&
-			  GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-	doBC[idim*2+1] = (((sGHex->GFSym[vi][idim*2+1]==GFSYM_NOSYM)||
-			   (sGHex->GFSym[vi][idim*2+1]==GFSYM_UNSET)) &&
-			  GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
-	lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
-      }
+      dirloop1 = 2*(dir-1)+1;
+      dirloop2 = dirloop1+1;
+    }
+    else if (dir<0)
+    {
+      dirloop1 = 2*(-dir-1);
+      dirloop2 = dirloop1+1;
     }
     else
     {
-      for (idim=0;idim<dim;idim++) 
+      dirloop1 = 0;
+      dirloop2 = 2*dim;
+    }
+  
+    for (vi=first_var; vi<first_var+num_vars; vi++) 
+    {
+      /* 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++)
       {
-	doBC[idim*2]   = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-	doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
-	lssh[idim]     = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
+	lssh[idim] = (int)GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
       }
+
+      for (count=dirloop1;count<dirloop2;count++)
+      {
+	if (sGHex)
+	{
+	  doBC = (
+		  ( (sGHex->GFSym[vi][count]==GFSYM_NOSYM) ||
+		    (sGHex->GFSym[vi][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 = BndApplyFlat1Di(GH, 
+				   doBC, 
+				   lssh,
+				   sw,
+				   GH->data[vi][timelevel]); 
+	    break;
+	  case 2:
+	    berr = BndApplyFlat2Di(GH, 
+				   doBC, 
+				   lssh,
+				   sw,
+				   GH->data[vi][timelevel]); 
+	    break;
+	  case 3:
+	    berr = BndApplyFlat3Di(GH, 
+				   doBC, 
+				   lssh,
+				   sw,
+				   GH->data[vi][timelevel]); 
+	    break;
+	  default : 
+	    berr = -1; 
+	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, 
+		       "ApplyBndFlat: No BC for dim>3: grid variable %s",
+		       CCTK_VarName(vi));
+	  }
+	}
+      }
+      /*	  retval = (berr>-1) ? 0 : -1;*/
     }
+  
+    free(dstag);
+    free(lssh);
 
-    switch (dim) 
+    if (!stencil_array)
     {
-      case 1: berr = BndApplyFlat1Di(GH, 
-				     dim,
-				     doBC, 
-				     lssh,
-				     stencil,
-				     GH->data[vi][timelevel]); break;
-      case 2: berr = BndApplyFlat2Di(GH, 
-				     dim,
-				     doBC, 
-				     lssh,
-				     stencil,
-				     GH->data[vi][timelevel]); break;
-      case 3: berr = BndApplyFlat3Di(GH, 
-				     dim,
-				     doBC,
-				     lssh,
-				     stencil,
-				     GH->data[vi][timelevel]); break;     
-      default : berr = -1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", 
-		   "ApplyBndFlat: No BC for dim>3: grid variable '%s'",
-		   CCTK_VarName(vi));
+      free(sw);
     }
-    berr = (berr>-1) ? 0 : -1;
   }
-  
-  free(dstag);
-  free(doBC);
-  free(lssh);
-  
-  return(berr);
+
+  return retval;
 }
 
 
+
 /*@@
    @routine    BndApplyFlat3Di
    @date       Thu Mar  2 11:14:08 2000
@@ -364,8 +610,7 @@ int ApplyBndFlat(cGH *GH,
 @@*/
 
 int BndApplyFlat3Di(cGH *GH, 
-		    int  gdim, 
-		    int *doBC, 
+		    int doBC, 
 		    int *lssh,
 		    int *stencil, 
 		    CCTK_REAL *var) 
@@ -373,8 +618,11 @@ int BndApplyFlat3Di(cGH *GH,
   int i,j,k,sw;
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower x)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (j=0;j<GH->cctk_lsh[1];j++)
@@ -389,8 +637,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper x)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (j=0;j<GH->cctk_lsh[1];j++)
@@ -405,8 +656,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
   
   /*lower y */
-  if (doBC[2] == 1)
+  if (doBC == 2)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower y)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -421,8 +675,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
   
   /* upper y */
-  if (doBC[3] == 1)
+  if (doBC == 3)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper y)\n");
+#endif
     for (k=0;k<GH->cctk_lsh[2];k++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -437,8 +694,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
 
   /* lower z */
-  if (doBC[4] == 1)
+  if (doBC == 4)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower z)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -453,8 +713,11 @@ int BndApplyFlat3Di(cGH *GH,
   }
 
   /* upper z */
-  if (doBC[5] == 1)
+  if (doBC == 5)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper z)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (i=0;i<GH->cctk_lsh[0];i++)
@@ -487,8 +750,7 @@ int BndApplyFlat3Di(cGH *GH,
 
 @@*/
 int BndApplyFlat2Di(cGH *GH, 
-		    int  gdim, 
-		    int *doBC, 
+		    int doBC, 
 		    int *lssh,
 		    int *stencil, 
 		    CCTK_REAL *var) 
@@ -496,8 +758,11 @@ int BndApplyFlat2Di(cGH *GH,
   int i,j,sw;
     
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower x)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (sw=0;sw<stencil[0];sw++)
@@ -509,8 +774,11 @@ int BndApplyFlat2Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper x)\n");
+#endif
     for (j=0;j<GH->cctk_lsh[1];j++)
     {
       for (sw=0;sw<stencil[0];sw++)
@@ -522,8 +790,11 @@ int BndApplyFlat2Di(cGH *GH,
   }
 
   /* lower y */
-  if (doBC[2] == 1)
+  if (doBC == 2)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower y)\n");
+#endif
     for (i=0;i<GH->cctk_lsh[0];i++)
     {
       for (sw=0;sw<stencil[1];sw++)
@@ -535,8 +806,11 @@ int BndApplyFlat2Di(cGH *GH,
   }
   
   /* upper y */
-  if (doBC[3] == 1)
+  if (doBC == 3)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper y)\n");
+#endif
     for (i=0;i<GH->cctk_lsh[0];i++)
     {
       for (sw=0;sw<stencil[1];sw++)
@@ -568,16 +842,18 @@ int BndApplyFlat2Di(cGH *GH,
 @@*/
 
 int BndApplyFlat1Di(cGH *GH, 
-		    int  gdim, 
-		    int *doBC, 
+		    int doBC, 
 		    int *lssh,
 		    int *stencil, 
 		    CCTK_REAL *var) 
 { int sw;
 
   /* lower x */
-  if (doBC[0] == 1)
+  if (doBC == 0)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (lower x)\n");
+#endif
     for (sw=0;sw<stencil[0];sw++)
     {
       var[CCTK_GFINDEX1D(GH,sw)]
@@ -586,8 +862,11 @@ int BndApplyFlat1Di(cGH *GH,
   }
   
   /* upper x */
-  if (doBC[1] == 1)
+  if (doBC == 1)
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying flat boundary (upper x)\n");
+#endif
     for (sw=0;sw<stencil[0];sw++)
     {
       var[CCTK_GFINDEX1D(GH,GH->cctk_lsh[0]-sw-1)]
diff --git a/src/RadiationBoundary.c b/src/RadiationBoundary.c
index dcc853f..68ab210 100644
--- a/src/RadiationBoundary.c
+++ b/src/RadiationBoundary.c
@@ -8,7 +8,7 @@
    @version $Header$
  @@*/
 
-/*#define DEBUG_BOUND*/
+/*#define DEBUG_BOUNDARY*/
  
 #include <stdio.h>
 #include <assert.h>
@@ -34,10 +34,9 @@ CCTK_FILEVERSION(CactusBase_Boundary_RadiationBoundaryWrappers_c)
 /* Internal routine prototypes */
 
 static int BndApplyRadiative3Di(cGH *GH, 
-				int  gdim,
-				int *sw,  
-				int *doBC, 
+				int doBC, 
 				int *lssh,
+				int *stencil_size,  
 				CCTK_REAL *dxyz, 
 				CCTK_REAL dt, 
 				CCTK_REAL *var_n, 
@@ -50,12 +49,14 @@ static int BndApplyRadiative3Di(cGH *GH,
 				CCTK_REAL v0) ;
 
 static int ApplyBndRadiative(cGH *GH, 
-			     int *sw,
+			     int stencil,
+			     int *stencil_array,
 			     CCTK_REAL var0, 
 			     CCTK_REAL v0, 
 			     int var_current, 
 			     int var_previous,
-			     int num_vars);
+			     int num_vars,
+			     int dir);
 
 /* Local variables */
 
@@ -65,6 +66,46 @@ static int ApplyBndRadiative(cGH *GH,
  ********************************************************************/
 
 /*@@
+   @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 
+
+@@*/
+
+int BndRadiativeDirGI(cGH *GH, 
+		   int stencil_size,
+		   int dir,
+		   CCTK_REAL var0, 
+		   CCTK_REAL speed, 
+		   int gi, 
+		   int gi_p) 
+{
+  int numvars, first_vi, first_vi_p; 
+
+  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);
+
+}
+
+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)
+{
+  *ierr = BndRadiativeDirGI(GH, *stencil_size, *dir, *var0, *speed, *gi, *gi_p);
+}
+
+/*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
@@ -79,9 +120,9 @@ static int ApplyBndRadiative(cGH *GH,
 @@*/
 
 int BndRadiativeGI(cGH *GH, 
-		   int *sw,
+		   int *stencil_array,
 		   CCTK_REAL var0, 
-		   CCTK_REAL v0, 
+		   CCTK_REAL speed, 
 		   int gi, 
 		   int gi_p) 
 {
@@ -91,17 +132,72 @@ int BndRadiativeGI(cGH *GH,
   first_vi_p = CCTK_FirstVarIndexI(gi_p);
   numvars    = CCTK_NumVarsInGroupI(gi);
 
-  return ApplyBndRadiative(GH,sw,var0,v0,first_vi,first_vi_p,numvars);
+  return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,first_vi,first_vi_p,numvars,0);
 
 }
 
 void CCTK_FCALL CCTK_FNAME(BndRadiativeGI)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, int *gi, int *gi_p)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, int *gi, int *gi_p)
 {
-  *ierr = BndRadiativeGI(GH, sw, *var0, *v0, *gi, *gi_p);
+  *ierr = BndRadiativeGI(GH, stencil_array, *var0, *speed, *gi, *gi_p);
 }
 
+
+/* ===================================================================== */
+
+
+/*@@
+   @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 
+
+@@*/
+
+int BndRadiativeDirGN(cGH *GH, 
+		      int stencil_size,
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL speed, 
+		      const char *gn, 
+		      const char *gn_p)
+{
+  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 
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
+	      "BndRadiativeDirGN: Grid variable groups %s or %s not found",
+	      gn, gn_p);
+    retval = -1;
+  }
+
+  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)
+{
+  TWO_FORTSTRINGS_CREATE(gn, gn_p)
+  *ierr = BndRadiativeDirGN(GH, *stencil_size, *dir, *var0, *speed, gn, gn_p);
+  free(gn);
+  free(gn_p);
+  return;
+}
+
+
 /*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
@@ -116,9 +212,9 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeGI)
 @@*/
 
 int BndRadiativeGN(cGH *GH, 
-		   int *sw, 
+		   int *stencil_array, 
 		   CCTK_REAL var0, 
-		   CCTK_REAL v0, 
+		   CCTK_REAL speed, 
 		   const char *gn, 
 		   const char *gn_p)
 {
@@ -127,7 +223,7 @@ int BndRadiativeGN(cGH *GH,
   int gi_p = CCTK_GroupIndex(gn_p);
   
   if ((gi>-1)&&(gi_p>-1))
-    retval = BndRadiativeGI(GH, sw, var0, v0, gi, gi_p);
+    retval = BndRadiativeGI(GH, stencil_array, var0, speed, gi, gi_p);
   else 
   {
     CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
@@ -140,18 +236,58 @@ int BndRadiativeGN(cGH *GH,
 }
  
 void CCTK_FCALL CCTK_FNAME(BndRadiativeGN)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, TWO_FORTSTRINGS_ARGS)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
 {
   TWO_FORTSTRINGS_CREATE(gn, gn_p)
-  *ierr = BndRadiativeGN(GH, sw, *var0, *v0, gn, gn_p);
+  *ierr = BndRadiativeGN(GH, stencil_array, *var0, *speed, gn, gn_p);
   free(gn);
   free(gn_p);
   return;
 }
 
+
+/* ===================================================================== */
+
+
 /*@@
-   @routine    BndRadiativeGI
+   @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 
+
+@@*/
+
+int BndRadiativeDirVI(cGH *GH, 
+		      int stencil_size, 
+		      int dir,
+		      CCTK_REAL var0, 
+		      CCTK_REAL speed, 
+		      int vi, 
+		      int vi_p) 
+{
+  return ApplyBndRadiative(GH,stencil_size,NULL,var0,speed,vi,vi_p,1,dir);
+}
+
+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)
+{
+  *ierr = BndRadiativeDirVI(GH, *stencil_size, *dir, *var0, 
+			    *speed, *vi, *vi_p);
+  return;
+}
+
+
+
+/*@@
+   @routine    BndRadiativeVI
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
    @desc 
@@ -165,23 +301,81 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeGN)
 @@*/
 
 int BndRadiativeVI(cGH *GH, 
-                    int *sw, 
+                    int *stencil_array, 
                     CCTK_REAL var0, 
-                    CCTK_REAL v0, 
+                    CCTK_REAL speed, 
                     int vi, 
                     int vi_p) 
 {
-  return ApplyBndRadiative(GH,sw,var0,v0,vi,vi_p,1);
+  return ApplyBndRadiative(GH,-1,stencil_array,var0,speed,vi,vi_p,1,0);
 }
 
 void CCTK_FCALL CCTK_FNAME(BndRadiativeVI)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, int *vi, int *vi_p)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, int *vi, int *vi_p)
 {
-  *ierr = BndRadiativeVI(GH, sw, *var0, *v0, *vi, *vi_p);
+  *ierr = BndRadiativeVI(GH, stencil_array, *var0, *speed, *vi, *vi_p);
   return;
 }
 
+
+/* ======================================================================= */
+
+
+/*@@
+   @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)
+  {
+    retval = BndRadiativeDirVI(GH, stencil_size, dir, var0, speed, vi, vi_p);
+  }
+  else
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
+	      "BndRadiativeDirVN: Grid variable %s or %s not found",
+	      vn, vn_p);
+    retval = -1;
+  }
+
+  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)
+{
+  TWO_FORTSTRINGS_CREATE(vn,vn_p)
+  *ierr = BndRadiativeDirVN(GH, *stencil_size, *dir, *var0, *speed, vn, vn_p);
+  free(vn);
+  free(vn_p);
+  return;
+}  
+
+
 /*@@
    @routine    BndRadiativeGI
    @date       Tue Jul 18 18:04:07 2000
@@ -197,9 +391,9 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVI)
 @@*/
 
 int BndRadiativeVN(cGH *GH, 
-                   int *sw,
+                   int *stencil_array,
                    CCTK_REAL var0, 
-                   CCTK_REAL v0, 
+                   CCTK_REAL speed, 
                    const char *vn, 
                    const char *vn_p) 
 {
@@ -209,7 +403,7 @@ int BndRadiativeVN(cGH *GH,
   
   if (vi>=0 && vi_p>=0)
   {
-    retval = BndRadiativeVI(GH, sw, var0, v0, vi, vi_p);
+    retval = BndRadiativeVI(GH, stencil_array, var0, speed, vi, vi_p);
   }
   else
   {
@@ -224,11 +418,11 @@ int BndRadiativeVN(cGH *GH,
 }
 
 void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
-     (int *ierr, cGH *GH, int *sw, CCTK_REAL *var0, 
-      CCTK_REAL *v0, TWO_FORTSTRINGS_ARGS)
+     (int *ierr, cGH *GH, int *stencil_array, CCTK_REAL *var0, 
+      CCTK_REAL *speed, TWO_FORTSTRINGS_ARGS)
 {
   TWO_FORTSTRINGS_CREATE(vn,vn_p)
-  *ierr = BndRadiativeVN(GH, sw, *var0, *v0, vn, vn_p);
+  *ierr = BndRadiativeVN(GH, stencil_array, *var0, *speed, vn, vn_p);
   free(vn);
   free(vn_p);
   return;
@@ -240,7 +434,7 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
  *********************     Local Routines   *************************
  ********************************************************************/
 /*@@
-   @routine    BndRadiativeGI
+   @routine    ApplyBndRadiative
    @date       Tue Jul 18 18:04:07 2000
    @author     Gerd Lanfermann
    @desc 
@@ -254,44 +448,117 @@ void CCTK_FCALL CCTK_FNAME(BndRadiativeVN)
 @@*/
 
 static int ApplyBndRadiative(cGH *GH, 
-			     int *sw,
+			     int stencil,
+			     int *stencil_array,
 			     CCTK_REAL var0, 
-			     CCTK_REAL v0, 
+			     CCTK_REAL speed, 
 			     int var_current, 
 			     int var_previous,
-			     int num_vars) 
+			     int num_vars,
+			     int dir) 
 {
-  SymmetryGHex *sGHex;
   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 *doBC, *dstag, *lssh;
+  int dirloop1,dirloop2;
+  int doBC;
+  int retval;
+  int *sw;
+  int *dstag;
+  int *lssh;
   CCTK_REAL dx[3], dt;
+  SymmetryGHex *sGHex;
+
+#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
+  
+  /* Get group dimension */
+  dim = CCTK_GroupDimFromVarI(var_current);
 
-  /* See if we have a symmetry array */
-  symmetry_handle = CCTK_GHExtensionHandle("Symmetry");
-  if (symmetry_handle < 0)
+  if (dim<3 || dim>3)
   {
-    sGHex = NULL;
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+	       "ApplyBndRadiative: Invalid grid function dimension %d",
+	       dim);
+    retval=-1;
   }
   else
   {
-    sGHex  = (SymmetryGHex*)GH->extensions[symmetry_handle];
-  }
+    /* 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];
+    }
 
-  /* Get group dimension */
-  dim = CCTK_GroupDimFromVarI(var_current);
+    /* 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
+    {
 
-  /* Radiative boundaries need the underlying Cartesian coordinates */
-  switch (dim) 
-  {
+      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;
@@ -306,130 +573,148 @@ static int ApplyBndRadiative(cGH *GH,
       ri = CCTK_CoordIndex(-1,"r","spher3d");
       break;
     default:
-      CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		 "Radiative boundary only for dim=3: grid varible '%s'",
+      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));
 
-  /* Allocate temporary arrays */
-  doBC  = (int *)malloc((2*dim)*sizeof(int)); 
-  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);
 
-  /* 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);
 
-  /* 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;
-  }
+      /* Apply condition if:
+	 + boundary is not a symmetry boundary 
+	 (no symmetry or unset(=unsed))
+	 + boundary is a physical boundary
+	 + have enough grid points            */
 
-  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
-    */
-    if (sGHex)
-    {
-      for (idim=0;idim<dim;idim++) 
+      for (idim=0;idim<dim;idim++)
       {
-        doBC[idim*2]   = (((sGHex->GFSym[vi1][idim*2]==GFSYM_NOSYM)||
-                           (sGHex->GFSym[vi1][idim*2]==GFSYM_UNSET)) &&
-                          GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-        doBC[idim*2+1] = (((sGHex->GFSym[vi1][idim*2+1]==GFSYM_NOSYM)||
-                           (sGHex->GFSym[vi1][idim*2+1]==GFSYM_UNSET)) &&
-                          GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
 	/* FIXME: THIS LINE FAILS WAVE TESTSUITES */
         /*lssh[idim] = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];*/
 	lssh[idim] = GH->cctk_lsh[idim];
       }
-    }
-    else
-    {
-      for (idim=0;idim<dim;idim++) 
+
+      for (count=dirloop1;count<dirloop2;count++)
       {
-        doBC[idim*2]   = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2]);
-        doBC[idim*2+1] = (GH->cctk_lsh[idim]>1 && GH->cctk_bbox[idim*2+1]);
-        lssh[idim]     = GH->cctk_lssh[CCTK_LSSH_IDX(dstag[idim],idim)];
+	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);
 
-    switch (dim) 
+    if (!stencil_array)
     {
-      case 1: 
-	berr = -1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		   "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__,"Boundary",
-		   "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,
-	   dim,
-	   sw,		          /* stencil_width */
-	   doBC,		  /* flag applying low/up BCs */
-	   lssh, 	          /* local shape, respects staggering */
-	   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 */
-	   v0);		          /* wave speed */
-	break;
-      default : 
-	berr=-1; 
-	CCTK_VWarn(1,__LINE__,__FILE__,"Boundary",
-		   "No Radiative BC for dim>3, grid variable '%s'",
-		   CCTK_VarName(vi1));
+      free(sw);
     }
-    berr=(berr<0)?-1:0;
   }
-  
-  free(lssh);
-  free(dstag);
-  free(doBC);
-  
-  return ierr;
-}
-
 
+  return retval;
+}
 
 /*@@
    @routine BndApplyRadiative3Di
@@ -515,10 +800,9 @@ static int ApplyBndRadiative(cGH *GH,
 @@*/
 
 static int BndApplyRadiative3Di(cGH *GH,
-				int  gdim,
-				int *sw,  
-				int *doBC, 
+				int doBC, 
 				int *lssh,
+				int *sw,  
 				CCTK_REAL *dxyz, 
 				CCTK_REAL dt, 
 				CCTK_REAL *var_n, 
@@ -528,7 +812,7 @@ static int BndApplyRadiative3Di(cGH *GH,
 				CCTK_REAL *z, 
 				CCTK_REAL *r,
 				CCTK_REAL var0, 
-				CCTK_REAL v0) 
+				CCTK_REAL speed) 
 {
 
   DECLARE_CCTK_PARAMETERS
@@ -588,7 +872,7 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Find Courant parameters. */
 
-  dtv  = v0*dt;
+  dtv  = speed*dt;
   dtvh = half*dtv;
 
   dtvvar0 = dtv*var0;
@@ -610,11 +894,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Lower x-bound: x(2,:,:) --> x[xgp2] */
 
-  if (doBC[0]==1) {
-    for (k=0;k<lssh2;k++) {
-      for (j=0;j<lssh1;j++) {
-	for (i=sw0-1;i>=0;i--) {
-
+  if (doBC == 0) 
+  {
+#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);
@@ -633,11 +923,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -670,10 +963,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */
 
-  if (doBC[1]==1){
-    for (k=0;k<lssh2;k++) {
-      for (j=0;j<lssh1;j++) {
-	for (i=lssh0-sw0;i<lssh0;i++) {
+  if (doBC == 1)
+  {
+#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);
@@ -693,11 +993,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -730,8 +1033,11 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Lower y-bound */
 
-  if (doBC[2] == 1) 
+  if (doBC == 2) 
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (lower y)\n");
+#endif
     for (k=0;k<lssh2;k++) 
     {
       for (i=0;i<lssh0;i++) 
@@ -756,11 +1062,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -793,8 +1102,11 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Upper y bound */
 
-  if (doBC[3] == 1) 
+  if (doBC == 3) 
   {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (upper y)\n");
+#endif
     for (k=0;k<lssh2;k++) 
     {
       for (i=0;i<lssh0;i++) 
@@ -860,10 +1172,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Lower z-bound */
 
-  if (doBC[4]==1) {
-    for (j=0;j<lssh1;j++) {
-      for (i=0;i<lssh0;i++) {
-	for (k=sw2-1;k>=0;k--) {
+  if (doBC==4) 
+  {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (lower z)\n");
+#endif
+    for (j=0;j<lssh1;j++) 
+    {
+      for (i=0;i<lssh0;i++) 
+      {
+	for (k=sw2-1;k>=0;k--) 
+	{
 
 	  zgp0 = CCTK_GFINDEX3D(GH,i,j,k  );
 	  zgp1 = CCTK_GFINDEX3D(GH,i,j,k+1);
@@ -883,11 +1202,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
@@ -920,10 +1242,17 @@ static int BndApplyRadiative3Di(cGH *GH,
 
   /* Upper z-bound */
 
-  if (doBC[5] == 1) {
-    for (j=0;j<lssh1;j++) {
-      for (i=0;i<lssh0;i++) {
-	for (k=lssh2-sw2;k<lssh2;k++) {
+  if (doBC ==  5) 
+  {
+#ifdef DEBUG_BOUNDARY
+    printf("Applying radiative boundary (upper z)\n");
+#endif
+    for (j=0;j<lssh1;j++) 
+    {
+      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);
@@ -943,11 +1272,14 @@ static int BndApplyRadiative3Di(cGH *GH,
           ri0 = 1.0/r0;
           ri1 = 1.0/r1;
 
-          if (fac==0) {
+          if (fac==0) 
+	  {
 
              H = 0;
 
-          } else {
+          } 
+	  else 
+	  {
 
              H = 0;
 
diff --git a/src/ScalarBoundary.c b/src/ScalarBoundary.c
index 02098d5..d905a51 100644
--- a/src/ScalarBoundary.c
+++ b/src/ScalarBoundary.c
@@ -132,7 +132,7 @@ void CCTK_FCALL CCTK_FNAME(BndScalarVI)
   return;
 }
 
-
+/* ===================================================================== */
 
  /*@@
    @routine    BndScalarDirGI
@@ -150,10 +150,10 @@ void CCTK_FCALL CCTK_FNAME(BndScalarVI)
 @@*/
 
 int BndScalarDirGI(cGH *GH, 
-		int stencil_size, 
-		int dir,
-		CCTK_REAL scalar_val, 
-		int gi)
+		   int stencil_size, 
+		   int dir,
+		   CCTK_REAL scalar_val, 
+		   int gi)
 {
   int numvars, first_vi;
 
@@ -217,6 +217,7 @@ void CCTK_FCALL CCTK_FNAME(BndScalarGI)
 }
 
 
+/* ===================================================================== */
 
  /*@@
    @routine    BndScalarDirGN
@@ -469,7 +470,6 @@ static int ApplyBndScalar(cGH *GH,
 {
   int symmetry_handle;  /* handle for the optional symmetry structure */
   int vi, gi, dim;
-  int coord;
   int idim;
   int berr,ierr,retval=0;
   int dirloop1,dirloop2;
@@ -544,12 +544,10 @@ static int ApplyBndScalar(cGH *GH,
 	if (dir > 0)
 	{
 	  sw[dir-1] = stencil;
-	  coord = dir-1;
 	}
 	else
 	{
 	  sw[-dir-1] = stencil;
-	  coord = -dir-1;
 	}
       }
       else
@@ -606,12 +604,12 @@ static int ApplyBndScalar(cGH *GH,
 	  doBC = (
 		  ( (sGHex->GFSym[vi][count]==GFSYM_NOSYM) ||
 		    (sGHex->GFSym[vi][count]==GFSYM_UNSET) ) &&
-		  GH->cctk_lsh[coord]>1 && GH->cctk_bbox[count]) ?
+		  GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
 	    count : -1;
 	}
 	else
 	{
-	  doBC  = (GH->cctk_lsh[coord]>1 && GH->cctk_bbox[count]) ?
+	  doBC  = (GH->cctk_lsh[count/2]>1 && GH->cctk_bbox[count]) ?
 	    count : -1 ;
 	}
 	
@@ -645,8 +643,8 @@ static int ApplyBndScalar(cGH *GH,
 	    break;     
 	  default : 
 	    berr = -1; 
-	    CCTK_VWarn(1,__LINE__,__FILE__,"Boundary", 
-		       "No scalar boundary for dim>3: grid variable '%s'",
+	    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, 
+		       "ApplyBndScalar: No scalar boundary for dim>3: grid variable '%s'",
 		       CCTK_VarName(vi));
 	  }
 	}