* major changes to driver routines to find multiple horizons in parallel

  across multiple processors -- see src/driver/README.parallel for details
* drop convergence checks on ||Delta_h|| in param.ccl because they don't
  fit well with parallelization changes

==> With this changes, AHFinderDirect is now (I think)
    multiprocessor-ready!!


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@946 f88db872-0e4f-0410-b76b-b9085cfa78c5

36 files changed, 2463 insertions, 1135 deletions

diff --git a/README b/README
index 4d4de14..1d3c93d 100644
--- a/README
+++ b/README
@@ -29,7 +29,7 @@ See below in this file for notes on compiling this thorn.
 Copyright
 =========
 
-This thorn is copyright (C) 2001-2002
+This thorn is copyright (C) 2001-2003
 by Jonathan Thornburg <jthorn@aei.mpg.de>.
 
 This thorn is free software; you can redistribute it and/or modify
@@ -70,7 +70,7 @@ Most of this thorn's relativity code is machine-generated using Maple
 relativity code.
 
 By default, this thorn doesn't use any external libraries.  However,
-if HAVE_DENSE_JACOBIAN__LAPACK is defined in src/include/configure.h,
+if HAVE_DENSE_JACOBIAN__LAPACK is defined in src/include/config.h,
 then this thorn uses the LAPACK library (which in turn uses the BLAS
 library), so you will need to configure your Cactus to use LAPACK.
 Instructions on doing this are in the file README.library in this
diff --git a/README.library b/README.library
index 2b7fad0..93f2d3c 100644
--- a/README.library
+++ b/README.library
@@ -5,6 +5,9 @@ $Header$
 By default this thorn doesn't use any external libraries, and you
 can ignore the instructions in this file.
 
+
+Using the LAPACK and BLAS Libraries
+===================================
 However, if this thorn is configured to use the LAPACK and BLAS
 libraries (see README for details), then you need to configure Cactus
 to use these libraries; this file describes how to do this.
@@ -54,3 +57,17 @@ or as assignments in your ~/.cactus/config file,
    LAPACK=yes
    LAPACK_EXTRA_LIBS=g2c
    LAPACK_EXTRA_LIB_DIRS=/usr/lib/gcc-lib/i386-redhat-linux/2.96
+
+
+Compiler Version Compatability
+==============================
+All of Cactus -- including any external libraries -- need not be
+compiled with the same compilers, but the compilers must be link-compatible.
+In practice this isn't usually a problem.  However, as an example of what
+*not* to do, consider using libraries compiled with a gcc 2.* version
+(e.g. the system default libraries on many GNU/Linux distributions as
+of late 2002), combined with thorns compiled with a gcc 3.* version.
+gcc 2.* and 3.* are *not* link-compatible, so if you're lucky this
+combination will give all sorts of wierd errors in linking.  If you're
+unlucky it will link ok but then crash (or even worse, just give wrong
+results) at run-time.
diff --git a/interface.ccl b/interface.ccl
index 8e2dfe6..fd71fbc 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -12,28 +12,28 @@ inherits: Grid ADMBase StaticConformal
 #
 protected:
 
-#
-# all the remaining diagnostics are arrays subscripted by the
-# (1-origin) horizon number  hn ; subscript 0 is unused
-#
-int BH_diagnostics__int TYPE=array DIM=1 SIZE=N_horizons+1 \
-                        DISTRIB=constant TIMELEVELS=1
-{
-# this is a Boolean flag: 0=false, 1=true
-AH_found	# was this AH found the last time we searched for it?
-		# (if this thorn is active, we search at each time step)
-} "CCTK_INT diagnostics calculated for each apparent horizon found"
-
-#
-# these diagnostics are only defined for those apparent horizons
-# where  AH_found[hn]  is true
-#
-real BH_diagnostics__real TYPE=array DIM=1 SIZE=N_horizons+1 \
-                          DISTRIB=constant TIMELEVELS=1
-{
-# FIXME: it would be nice to gather these together in some way
-centroid_x, centroid_y, centroid_z	# centroid position
-
-area		# area of apparent horizon
-m_irreducible	# irreducible mass = sqrt(area/(16*pi))
-} "CCTK_REAL diagnostics calculated for each apparent horizon found"
+###
+### all the remaining diagnostics are arrays subscripted by the
+### (1-origin) horizon number  hn ; subscript 0 is unused
+###
+##int BH_diagnostics__int TYPE=array DIM=1 SIZE=N_horizons+1 \
+##                        DISTRIB=constant TIMELEVELS=1
+##{
+### this is a Boolean flag: 0=false, 1=true
+##AH_found	# was this AH found the last time we searched for it?
+##		# (if this thorn is active, we search at each time step)
+##} "CCTK_INT diagnostics calculated for each apparent horizon found"
+##
+###
+### these diagnostics are only defined for those apparent horizons
+### where  AH_found[hn]  is true
+###
+##real BH_diagnostics__real TYPE=array DIM=1 SIZE=N_horizons+1 \
+##                          DISTRIB=constant TIMELEVELS=1
+##{
+### FIXME: it would be nice to gather these together in some way
+##centroid_x, centroid_y, centroid_z	# centroid position
+##
+##area		# area of apparent horizon
+##m_irreducible	# irreducible mass = sqrt(area/(16*pi))
+##} "CCTK_REAL diagnostics calculated for each apparent horizon found"
diff --git a/par/Kerr-tiny.par b/par/Kerr-tiny.par
index 90b4a08..3c91a90 100644
--- a/par/Kerr-tiny.par
+++ b/par/Kerr-tiny.par
@@ -35,7 +35,7 @@ Exact::Kerr_KerrSchild__spin = 0.6
 
 ########################################
 
-ActiveThorns = "LocalInterp PUGHInterp AHFinderDirect"
+ActiveThorns = "LocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 AHFinderDirect::find_AHs_at_poststep = "false"	# no time evolution
 AHFinderDirect::h_base_file_name     = "Kerr-tiny.h"
 
diff --git a/par/Kerr.par b/par/Kerr.par
index abcc41c..c9e1077 100644
--- a/par/Kerr.par
+++ b/par/Kerr.par
@@ -32,7 +32,7 @@ Exact::Kerr_KerrSchild__spin = 0.6
 
 ########################################
 
-ActiveThorns = "LocalInterp PUGHInterp AHFinderDirect"
+ActiveThorns = "LocalInterp PUGHInterp PUGHReduce AHFinderDirect"
 
 AHFinderDirect::find_AHs_at_postinitial = "true"
 AHFinderDirect::find_AHs_at_poststep    = "false"
diff --git a/param.ccl b/param.ccl
index f5ce1e5..e54df30 100644
--- a/param.ccl
+++ b/param.ccl
@@ -43,14 +43,21 @@ boolean find_AHs_at_poststep \
 keyword method "what should this thorn do for each apparent horizon?"
 {
 # these options are mostly for testing/debugging
-"evaluate expansion"	  :: "evaluate the LHS function Theta(h)"
-"test expansion Jacobian" :: "compute/print the J[Theta(h)] Jacobian matrix \
-			      (possibly in several ways, depending on \
-			      the  test_all_Jacobian_methods  parameter"
+# ... in a multiprocessor Cactus run, the horizons are done sequentually
+#     on processor #0; the other processors do dummy computations
+"evaluate expansions"	  :: "evaluate the LHS function Theta(h)"
+
+# ... in a multiprocessor Cactus run, the Jacobian is computed on
+#     processor #0; the other processors do dummy computations
+"test expansion Jacobians" :: \
+  "compute/print horizon 1's J[Theta(h)] Jacobian matrix (possibly in \
+   multiple ways, depending on the  test_all_Jacobian_methods  parameter)"
 
 # this is for normal apparent horizon finding
-"find horizon"		  :: "find the apparent horizon"
-} "find horizon"
+# ... in a multiprocessor Cactus run, the horizons are done in parallel
+#     across processors; see src/driver/README.parallel for details
+"find horizons"		  :: "find the apparent horizon"
+} "find horizons"
 
 #
 # We support searching for up to N_horizons distinct apparent horizons
@@ -72,10 +79,8 @@ int N_horizons "number of apparent horizons to search for"
 } 1
 
 #
-# This parameter controls which (how many) informational (non-error)
-# messages are printed describing the operation of the apparent horizon
-# finder.  This is analogous to the -W Cactus command-line option, except
-# that this parameter controls (only) *non*-error messages.
+# This parameter controls how verbose this thorn is in printing
+# informational (non-error) messages describing what it's doing.
 #
 keyword verbose_level \
   "controls which (how many) messages to print describing AH finding"
@@ -200,35 +205,20 @@ int max_Newton_iterations__subsequent \
 # move in any single Newton iteration (i.e. the infinity-norm of Delta_h)
 # to <= this fraction of the mean horizon radius
 #
-real max_Delta_h_over_h \
+real max_allowable_Delta_h_over_h \
   "don't let horizon move > this fraction of mean radius in a Newton iteration"
 {
 (0.0:* :: "any positive real number"
 } 0.1
 
 #
-# we declare convergence if *either* of the following two criteria are met
+# convergence tolerance for the Newton iteration
 #
-real Theta_norm_for_convergence "declare convergence if ||Theta||_inf <= this"
+real Theta_norm_for_convergence \
+  "we declare the horizon to be found if ||Theta||_infinity <= this"
 {
 (0.0:* :: "any positive real number"
 } 1.0e-8
-real Delta_h_norm_for_convergence \
-  "declare convergence after any Newton step with ||Delta_h||_inf <= this"
-{
-(0.0:* :: "any positive real number"
-} 1.0e-8
-
-#
-# This only needs to be set to true for very careful convergence studies
-# etc.  On the other hand, setting it to true probably only slows down
-# the apparent horizon finder by a few percent.
-#
-boolean final_Theta_update_if_Delta_h_converged \
-  "should we do a final Theta(h) update if we terminate the \
-   Newton iteration by the small-||Delta_h|| convergence criterion?"
-{
-} "false"
 
 ################################################################################
 
diff --git a/run/test-ahfinderdirect/Kerr/Kerr-eval.par b/run/test-ahfinderdirect/Kerr/Kerr-eval.par
new file mode 100644
index 0000000..013d148
--- /dev/null
+++ b/run/test-ahfinderdirect/Kerr/Kerr-eval.par
@@ -0,0 +1,56 @@
+# This parameter file sets up Kerr/Kerr-Schild initial data, then
+# evalutes Theta for the "exact" initial guess.  The local coordinate
+# origin is deliberately de-centered with respect to the black hole,
+# to make this a non-trivial test for the apparent horizon finder.
+
+# flesh
+cactus::cctk_itlast = 0
+
+ActiveThorns = "PUGH"
+driver::ghost_size = 3
+driver::global_nx = 31
+driver::global_ny = 31
+driver::global_nz = 19
+
+ActiveThorns = "CartGrid3D"
+grid::domain = "bitant"
+grid::avoid_origin = "false"
+grid::type = "byspacing"
+grid::dxyz = 0.2
+
+ActiveThorns = "ADMBase ADMCoupling StaticConformal CoordGauge Exact"
+ADMBase::initial_lapse = "exact"
+ADMBase::initial_shift = "exact"
+ADMBase::initial_data = "exact"
+ADMBase::lapse_evolution_method = "static"
+ADMBase::shift_evolution_method = "static"
+ADMBase::metric_type = "physical"
+Exact::exact_model = "Kerr/Kerr-Schild"
+Exact::Kerr_KerrSchild__mass = 1.0
+Exact::Kerr_KerrSchild__spin = 0.6
+
+########################################
+
+ActiveThorns = "LocalInterp PUGHInterp PUGHReduce AHFinderDirect"
+
+AHFinderDirect::find_AHs_at_postinitial = "true"
+AHFinderDirect::find_AHs_at_poststep    = "false"
+AHFinderDIrect::method = "evaluate expansions"
+AHFinderDirect::print_timing_stats = "true"
+AHFinderDirect::verbose_level = "algorithm details"
+
+AHFinderDirect::how_often_to_output_Theta = 1
+AHFinderDirect::h_base_file_name     = "Kerr-eval.h"
+AHFinderDirect::Theta_base_file_name = "Kerr-eval.Theta"
+
+AHFinderDirect::N_horizons = 1
+AHFinderDirect::origin_x[1] = 0.5
+AHFinderDirect::origin_y[1] = 0.7
+AHFinderDirect::origin_z[1] = 0.0
+
+AHFinderDirect::initial_guess_method[1] = "Kerr/Kerr-Schild"
+AHFinderDirect::initial_guess__Kerr_Kerr__mass[1] = 1.0
+AHFinderDirect::initial_guess__Kerr_Kerr__spin[1] = 0.6
+AHFinderDirect::initial_guess__Kerr_Kerr__x_posn[1] = 0.0
+AHFinderDirect::initial_guess__Kerr_Kerr__y_posn[1] = 0.0
+AHFinderDirect::initial_guess__Kerr_Kerr__z_posn[1] = 0.0
diff --git a/run/test-ahfinderdirect/Kerr/Kerr-find.par b/run/test-ahfinderdirect/Kerr/Kerr-find.par
new file mode 100644
index 0000000..50b9720
--- /dev/null
+++ b/run/test-ahfinderdirect/Kerr/Kerr-find.par
@@ -0,0 +1,55 @@
+# This parameter file sets up Kerr/Kerr-Schild initial data, then
+# finds the apparent horizon in it.  The local coordinate system origin
+# and the initial guess are both deliberately de-centered with respect
+# to the black hole, to make this a non-trivial test for the apparent
+# horizon finder.
+
+# flesh
+cactus::cctk_itlast = 0
+
+ActiveThorns = "PUGH"
+driver::ghost_size = 3
+driver::global_nx = 31
+driver::global_ny = 31
+driver::global_nz = 19
+
+ActiveThorns = "CartGrid3D"
+grid::domain = "bitant"
+grid::avoid_origin = "false"
+grid::type = "byspacing"
+grid::dxyz = 0.2
+
+ActiveThorns = "ADMBase ADMCoupling StaticConformal CoordGauge Exact"
+ADMBase::initial_lapse = "exact"
+ADMBase::initial_shift = "exact"
+ADMBase::initial_data = "exact"
+ADMBase::lapse_evolution_method = "static"
+ADMBase::shift_evolution_method = "static"
+ADMBase::metric_type = "physical"
+Exact::exact_model = "Kerr/Kerr-Schild"
+Exact::Kerr_KerrSchild__mass = 1.0
+Exact::Kerr_KerrSchild__spin = 0.6
+
+########################################
+
+ActiveThorns = "LocalInterp PUGHInterp PUGHReduce AHFinderDirect"
+
+AHFinderDirect::find_AHs_at_postinitial = "true"
+AHFinderDirect::find_AHs_at_poststep    = "false"
+AHFinderDirect::print_timing_stats = "true"
+AHFinderDirect::verbose_level = "algorithm details"
+
+AHFinderDirect::how_often_to_output_Theta = 1
+AHFinderDirect::h_base_file_name     = "Kerr-find.h"
+AHFinderDirect::Theta_base_file_name = "Kerr-find.Theta"
+
+AHFinderDirect::N_horizons = 1
+AHFinderDirect::origin_x[1] = 0.5
+AHFinderDirect::origin_y[1] = 0.7
+AHFinderDirect::origin_z[1] = 0.0
+
+AHFinderDirect::initial_guess_method[1] = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2
+AHFinderDirect::initial_guess__coord_sphere__y_center[1] =  0.3
+AHFinderDirect::initial_guess__coord_sphere__z_center[1] =  0.0
+AHFinderDirect::initial_guess__coord_sphere__radius[1] = 2.0
diff --git a/run/test-ahfinderdirect/Kerr/Kerr.par b/run/test-ahfinderdirect/Kerr/Kerr.par
new file mode 100644
index 0000000..1d9a09a
--- /dev/null
+++ b/run/test-ahfinderdirect/Kerr/Kerr.par
@@ -0,0 +1,55 @@
+# This parameter file sets up Kerr/Kerr-Schild initial data, then
+# finds the apparent horizon in it.  The local coordinate system origin
+# and the initial guess are both deliberately de-centered with respect
+# to the black hole, to make this a non-trivial test for the apparent
+# horizon finder.
+
+# flesh
+cactus::cctk_itlast = 0
+
+ActiveThorns = "PUGH"
+driver::ghost_size = 2
+driver::global_nx = 31
+driver::global_ny = 31
+driver::global_nz = 17
+
+ActiveThorns = "CartGrid3D"
+grid::domain = "bitant"
+grid::avoid_origin = "false"
+grid::type = "byspacing"
+grid::dxyz = 0.2
+
+ActiveThorns = "ADMBase ADMCoupling StaticConformal CoordGauge Exact"
+ADMBase::initial_lapse = "exact"
+ADMBase::initial_shift = "exact"
+ADMBase::initial_data = "exact"
+ADMBase::lapse_evolution_method = "static"
+ADMBase::shift_evolution_method = "static"
+ADMBase::metric_type = "physical"
+Exact::exact_model = "Kerr/Kerr-Schild"
+Exact::Kerr_KerrSchild__mass = 1.0
+Exact::Kerr_KerrSchild__spin = 0.6
+
+########################################
+
+ActiveThorns = "LocalInterp PUGHInterp PUGHReduce AHFinderDirect"
+
+AHFinderDirect::find_AHs_at_postinitial = "true"
+AHFinderDirect::find_AHs_at_poststep    = "false"
+AHFinderDirect::print_timing_stats = "true"
+AHFinderDirect::verbose_level = "algorithm details"
+
+AHFinderDirect::how_often_to_output_Theta = 1
+AHFinderDirect::h_base_file_name     = "Kerr.h"
+AHFinderDirect::Theta_base_file_name = "Kerr.Theta"
+
+AHFinderDirect::N_horizons = 1
+AHFinderDirect::origin_x[1] = 0.5
+AHFinderDirect::origin_y[1] = 0.7
+AHFinderDirect::origin_z[1] = 0.0
+
+AHFinderDirect::initial_guess_method[1] = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2
+AHFinderDirect::initial_guess__coord_sphere__y_center[1] =  0.3
+AHFinderDirect::initial_guess__coord_sphere__z_center[1] =  0.0
+AHFinderDirect::initial_guess__coord_sphere__radius[1] = 2.0
diff --git a/run/test-ahfinderdirect/Kerr/Kerr3-find.par b/run/test-ahfinderdirect/Kerr/Kerr3-find.par
new file mode 100644
index 0000000..3885c21
--- /dev/null
+++ b/run/test-ahfinderdirect/Kerr/Kerr3-find.par
@@ -0,0 +1,74 @@
+# This parameter file sets up Kerr/Kerr-Schild initial data, then
+# finds the apparent horizon in it.  The local coordinate system origin
+# and the initial guess are both deliberately de-centered with respect
+# to the black hole, to make this a non-trivial test for the apparent
+# horizon finder.
+
+# flesh
+cactus::cctk_itlast = 0
+
+ActiveThorns = "PUGH"
+driver::ghost_size = 3
+driver::global_nx = 31
+driver::global_ny = 31
+driver::global_nz = 19
+
+ActiveThorns = "CartGrid3D"
+grid::domain = "bitant"
+grid::avoid_origin = "false"
+grid::type = "byspacing"
+grid::dxyz = 0.2
+
+ActiveThorns = "ADMBase ADMCoupling StaticConformal CoordGauge Exact"
+ADMBase::initial_lapse = "exact"
+ADMBase::initial_shift = "exact"
+ADMBase::initial_data = "exact"
+ADMBase::lapse_evolution_method = "static"
+ADMBase::shift_evolution_method = "static"
+ADMBase::metric_type = "physical"
+Exact::exact_model = "Kerr/Kerr-Schild"
+Exact::Kerr_KerrSchild__mass = 1.0
+Exact::Kerr_KerrSchild__spin = 0.6
+
+########################################
+
+ActiveThorns = "LocalInterp PUGHInterp PUGHReduce AHFinderDirect"
+
+AHFinderDirect::find_AHs_at_postinitial = "true"
+AHFinderDirect::find_AHs_at_poststep    = "false"
+AHFinderDirect::print_timing_stats = "true"
+AHFinderDirect::verbose_level = "algorithm details"
+
+AHFinderDirect::how_often_to_output_Theta = 1
+AHFinderDirect::h_base_file_name     = "Kerr3-find.h"
+AHFinderDirect::Theta_base_file_name = "Kerr3-find.Theta"
+
+AHFinderDirect::N_horizons = 3
+
+AHFinderDirect::origin_x[1] = 0.5
+AHFinderDirect::origin_y[1] = 0.7
+AHFinderDirect::origin_z[1] = 0.0
+AHFinderDirect::initial_guess_method[1] = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[1] = -0.2
+AHFinderDirect::initial_guess__coord_sphere__y_center[1] =  0.3
+AHFinderDirect::initial_guess__coord_sphere__z_center[1] =  0.0
+AHFinderDirect::initial_guess__coord_sphere__radius[1] = 2.0
+
+AHFinderDirect::origin_x[2] = 0.9
+AHFinderDirect::origin_y[2] = 1.1
+AHFinderDirect::origin_z[2] = 0.0
+AHFinderDirect::initial_guess_method[1] = "coordinate sphere"
+AHFinderDirect::initial_guess__coord_sphere__x_center[2] = 0.1
+AHFinderDirect::initial_guess__coord_sphere__y_center[2] = 0.2
+AHFinderDirect::initial_guess__coord_sphere__z_center[2] = 0.0
+AHFinderDirect::initial_guess__coord_sphere__radius[2] = 1.7
+
+AHFinderDirect::origin_x[3] = 0.9
+AHFinderDirect::origin_y[3] = 1.1
+AHFinderDirect::origin_z[3] = 0.0
+AHFinderDirect::initial_guess_method[1] = "Kerr/Kerr-Schild"
+AHFinderDirect::initial_guess__Kerr_KerrSchild__x_posn[3] = 0.0
+AHFinderDirect::initial_guess__Kerr_KerrSchild__y_posn[3] = 0.0
+AHFinderDirect::initial_guess__Kerr_KerrSchild__z_posn[3] = 0.0
+AHFinderDirect::initial_guess__Kerr_KerrSchild__mass[3] = 1.0
+AHFinderDirect::initial_guess__Kerr_KerrSchild__spin[3] = 0.6
diff --git a/schedule.ccl b/schedule.ccl
index 8fc6720..a95879d 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,9 +1,9 @@
 # Schedule definitions for thorn AHFinderDirect
 # $Header$
 
-# these are small, so there's no harm in having storage for them always on
-STORAGE: BH_diagnostics__int
-STORAGE: BH_diagnostics__real
+### these are small, so there's no harm in having storage for them always on
+##STORAGE: BH_diagnostics__int
+##STORAGE: BH_diagnostics__real
 
 schedule AHFinderDirect_setup at CCTK_BASEGRID after SpatialCoordinates
 {
diff --git a/src/README b/src/README
index 0ee1180..b6cafd0 100644
--- a/src/README
+++ b/src/README
@@ -2,12 +2,12 @@ This is the top-level source directory for the AHFinderDirect thorn.
 See ../doc/ for further documentation.
 
 The main subdirectories under this directory are:
-driver/		high-level driver routines to solve the H(h) = 0
+driver/		high-level driver routines to solve the Theta(h) = 0
 		equations and interface to the rest of Cactus
 elliptic/	code to solve elliptic equations on the multipatch $S^2$
 gr/		relativity code; all knowledge of the actual apparent
 		horizon equation lives in the code in this directory
-gr.cg/		Maple-generated code to compute the H(h) function
+gr.cg/		Maple-generated code to compute the Theta(h) function
 		and its Jacobian coefficients
 patch/		basic multipatch infrastructure for storing and
 		finite differencing gridfns in angular coordinates
diff --git a/src/driver/BH_diagnostics.cc b/src/driver/BH_diagnostics.cc
new file mode 100644
index 0000000..ca411f8
--- /dev/null
+++ b/src/driver/BH_diagnostics.cc
@@ -0,0 +1,313 @@
+// BH_diagnostics.cc -- compute/print BH diagnostics
+// $Header$
+//
+// BH_diagnostics::BH_diagnostics - initialize a  struct BH_diagnostics
+//
+// compute_BH_diagnostics - compute BH diagnostics
+/// surface_integral - compute surface integral over a patch system
+//
+// print_BH_diagnostics - print a summary of BH diagnostics
+//
+// setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file
+// output_BH_diagnostics_fn - append BH diagnostics to file
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "config.h"
+#include "stdc.h"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "horizon_sequence.hh"
+#include "driver.hh"
+
+//******************************************************************************
+
+//
+// prototypes for functions local to this file
+//
+
+namespace {
+fp surface_integral(const patch_system& ps, int src_gfn,
+		    enum patch::integration_method method);
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function initializes a  struct BH_diagnostics  to all zeros.
+//
+BH_diagnostics::BH_diagnostics()
+	: centroid_x(0.0), centroid_y(0.0), centroid_z(0.0),
+	  mean_radius(0.0),
+	  circumference_xy(0.0), circumference_xz(0.0), circumference_yz(0.0),
+	  area(0.0),
+	  m_irreducible(0.0)
+{ }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// Given that an apparent horizon has been found, this function computes
+// various black hole diagnostics.
+//
+// Inputs (gridfns)
+// h		# ghosted
+// one		# nominal
+// global_[xyz]	# nominal
+//
+void compute_BH_diagnostics
+	(const patch_system& ps,
+	 const struct BH_diagnostics_info& BH_diagnostics_info,
+	 const struct verbose_info& verbose_info,
+	 struct BH_diagnostics& BH_diagnostics)
+{
+//
+// compute raw surface integrals
+//
+fp integral_one = surface_integral(ps, gfns::gfn__one,
+				   BH_diagnostics_info.integral_method);
+fp integral_h = surface_integral(ps, gfns::gfn__h,
+				 BH_diagnostics_info.integral_method);
+fp integral_x = surface_integral(ps, gfns::gfn__global_x,
+				 BH_diagnostics_info.integral_method);
+fp integral_y = surface_integral(ps, gfns::gfn__global_y,
+				 BH_diagnostics_info.integral_method);
+fp integral_z = surface_integral(ps, gfns::gfn__global_z,
+				 BH_diagnostics_info.integral_method);
+
+//
+// adjust integrals to take into account patch system not covering full sphere
+//
+switch	(ps.type())
+	{
+case patch_system::patch_system__full_sphere:
+	break;
+case patch_system::patch_system__plus_z_hemisphere:
+	integral_one *= 2.0;
+	integral_h *= 2.0;
+	integral_x *= 2.0;
+	integral_y *= 2.0;
+	integral_z = integral_one * ps.origin_z();
+	break;
+case patch_system::patch_system__plus_xy_quadrant_mirrored:
+case patch_system::patch_system__plus_xy_quadrant_rotating:
+	integral_one *= 4.0;
+	integral_h *= 4.0;
+	integral_x = integral_one * ps.origin_x();
+	integral_y = integral_one * ps.origin_y();
+	integral_z *= 4.0;
+	break;
+case patch_system::patch_system__plus_xz_quadrant_rotating:
+	integral_one *= 4.0;
+	integral_h *= 4.0;
+	integral_x = integral_one * ps.origin_x();
+	integral_y *= 4.0;
+	integral_z = integral_one * ps.origin_z();
+	break;
+case patch_system::patch_system__plus_xyz_octant_mirrored:
+case patch_system::patch_system__plus_xyz_octant_rotating:
+	integral_one *= 8.0;
+	integral_h *= 8.0;
+	integral_x = integral_one * ps.origin_x();
+	integral_y = integral_one * ps.origin_y();
+	integral_z = integral_one * ps.origin_z();
+	break;
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n"
+"                             (this should never happen!)"
+,
+		   int(ps.type()));				/*NOTREACHED*/
+	}
+
+BH_diagnostics.centroid_x = integral_x / integral_one;
+BH_diagnostics.centroid_y = integral_y / integral_one;
+BH_diagnostics.centroid_z = integral_z / integral_one;
+
+BH_diagnostics.area = integral_one;
+BH_diagnostics.circumference_xy = ps.xy_circumference
+    (gfns::gfn__h,
+     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+					     gfns::gfn__g_dd_33,
+     BH_diagnostics_info.integral_method);
+BH_diagnostics.circumference_xz = ps.xz_circumference
+    (gfns::gfn__h,
+     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+					     gfns::gfn__g_dd_33,
+     BH_diagnostics_info.integral_method);
+BH_diagnostics.circumference_yz = ps.yz_circumference
+    (gfns::gfn__h,
+     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+					     gfns::gfn__g_dd_33,
+     BH_diagnostics_info.integral_method);
+BH_diagnostics.mean_radius = integral_h / integral_one;
+BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
+}
+
+//******************************************************************************
+
+//
+// This function computes the surface integral of a gridfn over the
+// horizon.
+//
+namespace {
+fp surface_integral(const patch_system& ps, int src_gfn,
+		    enum patch::integration_method method)
+{
+return ps.integrate_gridfn
+	   (src_gfn,
+	    gfns::gfn__h,
+	    gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
+				gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
+						    gfns::gfn__g_dd_33,
+	    method);
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function prints a summary of the BH diagnostics, using CCTK_VInfo().
+//
+void print_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
+			  int N_horizons, int hn,
+			  const struct verbose_info& verbose_info)
+{
+if (verbose_info.print_physics_details)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "AH %d/%d: r=%g at (%f,%f,%f)",
+		   hn, N_horizons,
+		   double(BH_diagnostics.mean_radius),
+		   double(BH_diagnostics.centroid_x),
+		   double(BH_diagnostics.centroid_y),
+		   double(BH_diagnostics.centroid_z));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "AH %d/%d: area=%.10g m_irreducible=%.10g",
+		   hn, N_horizons,
+		   double(BH_diagnostics.area),
+		   double(BH_diagnostics.m_irreducible));
+	}
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function creates a BH-diagnostics output file, writes a suitable
+// header comment, and returns a stdio file pointer which can be used to
+// append data to the file.
+//
+FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
+				       int N_horizons, int hn)
+{
+const int N_file_name_buffer = 200;
+char file_name_buffer[N_file_name_buffer];
+
+snprintf(file_name_buffer, N_file_name_buffer,
+	 "%s.ah%d.%s",
+	 IO_info.BH_diagnostics_base_file_name,
+	 hn, IO_info.BH_diagnostics_file_name_extension);
+FILE *fileptr = fopen(file_name_buffer, "w");
+if (fileptr == NULL)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   setup_BH_diagnostics_output_file():\n"
+"        can't open BH-diagnostics output file\n"
+"        \"%s\"!"
+		   ,
+		   file_name_buffer);				/*NOTREACHED*/
+
+fprintf(fileptr, "# apparent horizon %d/%d\n", hn, N_horizons);
+fprintf(fileptr, "#\n");
+fprintf(fileptr, "# column 1 = cctk_iteration\n");
+fprintf(fileptr, "# column 2 = cctk_time\n");
+fprintf(fileptr, "# column 3 = centroid_x\n");
+fprintf(fileptr, "# column 4 = centroid_y\n");
+fprintf(fileptr, "# column 5 = centroid_z\n");
+fprintf(fileptr, "# column 6 = mean radius\n");
+fprintf(fileptr, "# column 7 = xy-plane circumference\n");
+fprintf(fileptr, "# column 8 = xz-plane circumference\n");
+fprintf(fileptr, "# column 9 = yz-plane circumference\n");
+fprintf(fileptr, "# column 10 = area\n");
+fprintf(fileptr, "# column 11 = m_irreducible\n");
+fflush(fileptr);
+
+return fileptr;
+}
+
+//******************************************************************************
+
+//
+// This function appends a BH-diagnostics line to a data file.  It
+// attempts to ensure that the newly-written line is flushed to disk,
+// so the output file can be examined while the Cactus run is still in
+// progress.
+//
+// (The "_fn" in the function name is to distinguish it from the Cactus
+// parameter  output_BH_diagnostics .)
+//
+// Arguments:
+// BH_diagnostics = The BH diagnostics to be written
+// fileptr = The stdio file pointer to append to
+//
+void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
+			      const struct IO_info& IO_info,
+			      int hn, FILE* fileptr)
+{
+assert(fileptr != NULL);
+
+fprintf(fileptr,
+//  cctk_iteration        mean radius                     area    m_irreducible
+//  ==  cctk_time         ======  {xy,xz,yz}-plane        ======  ======
+//  ==  ====  centroid_[xyz]      circumferences          ======  ======
+//  ==  ====  ==========  ======  ======================  ======  ======
+   "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
+	IO_info.time_iteration, double(IO_info.time),
+	BH_diagnostics.centroid_x, BH_diagnostics.centroid_y,
+				   BH_diagnostics.centroid_z,
+	BH_diagnostics.mean_radius,
+	BH_diagnostics.circumference_xy, BH_diagnostics.circumference_xz,
+					 BH_diagnostics.circumference_yz,
+	BH_diagnostics.area,
+	BH_diagnostics.m_irreducible);
+
+fflush(fileptr);
+}
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index f75e407..1af34a0 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -1,11 +1,15 @@
 // Newton.cc -- solve Theta(h) = 0 via Newton's method
 // $Header$
 //
-// Newton_solve - driver to solve Theta(h) = 0 via Newton's method
+// Newton_solve - driver to solve Theta(h) = 0 for all our horizons
+/// reduce_and_broadcast - reduce/broadcast status to other processors
+/// print_Theta_norms - print the Theta norms
+/// Newton_step - take the Newton step, scaling it down if it's too large
 //
 
 #include <stdio.h>
 #include <assert.h>
+#include <limits.h>
 #include <math.h>
 
 #include "util_Table.h"
@@ -34,117 +38,328 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
+#define DEBUG
+
 //******************************************************************************
 
 //
-// This function solves Theta(h) = 0 for h via Newton's method.
+// prototypes for functions local to this file
 //
-// Arguments:
-// initial_find_flag = true ==> This is the first attempt to find this
-//				horizon, or this is a subsequent attempt
-//				and the immediately previous attempt failed.
-//		       false ==> This isn't the first attempt to find this
-//				 horizon, and we found it successfully on
-//				 the immediately previous attempt.
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We find the horizon as usual, but don't write
-//			   any output files.
-// hn = the horizon number (used only in naming output files)
+
+namespace {
+bool reduce_and_broadcast(cGH* GH,
+			  int N_procs, int N_active_procs,
+			  int my_proc, bool my_active_flag,
+			  bool we_need_more_iterations,
+			  int hn, fp rms_norm, fp infinity_norm,
+			  int hn_buffer[],
+			  fp rms_norm_buffer[],
+			  fp infinity_norm_buffer[]);
+void print_Theta_norms(int N_active_procs,
+		       const int hn_buffer[],
+		       const fp rms_norm_buffer[],
+		       const fp infinity_norm_buffer[]);
+void Newton_step(patch_system& ps,
+		 fp max_allowable_Delta_h_over_h,
+		 const struct verbose_info& verbose_info);
+	  }
+
+//******************************************************************************
+
 //
-// Results:
-// This function returns a success flag: this is true if (and only if)
-// it found an h satisfying Theta(h) = 0 to within the error tolerances,
-// otherwise it's false.
-//
-bool Newton_solve(patch_system& ps,
-		  Jacobian& Jac,
-		  const struct cactus_grid_info& cgi,
-		  const struct geometry_info& gi,
-		  const struct Jacobian_info& Jacobian_info,
-		  const struct solver_info& solver_info,
-		  bool initial_find_flag,
-		  const struct IO_info& IO_info,
-		  bool active_flag,
-		  int hn, const struct verbose_info& verbose_info)
+// This function tries to finds each horizon assigned to this processor,
+// by solving Theta(h) = 0 via Newton's method.  For each horizon found,
+// it computes the BH diagnostics, optionally prints them (via CCTK_VInfo()),
+// and optionally writes them to a BH diagnostics output file.
+//
+// This function must be called synchronously across all processors,
+// and it operates synchronously, returning only when every processor
+// is done with all its horizons.  See ./README.parallel for a discussion
+// of the parallel/multiprocessor issues and algorithm.
+//
+// If this is a dummy processor (i.e. hs has no genuine horizons), then
+// my_AH_info will be a NULL pointer.
+//
+// FIXME: Since the BH diagnostics printing is done with CCTK_VInfo(),
+//        in a multiprocessor Cactus run the output will typically be
+//        lost for processors other than processor #0.
+//
+void Newton(cGH* GH,
+	    int N_procs, int N_active_procs, int my_proc,
+	    horizon_sequence& hs, struct AH_info* const my_AH_info[],
+	    const struct cactus_grid_info& cgi,
+	    const struct geometry_info& gi,
+	    const struct Jacobian_info& Jacobian_info,
+	    const struct solver_info& solver_info,
+	    const struct IO_info& IO_info,
+	    const struct BH_diagnostics_info& BH_diagnostics_info,
+	    const struct verbose_info& verbose_info)
 {
-const int max_Newton_iterations
-    = initial_find_flag ? solver_info.max_Newton_iterations__initial
-			: solver_info.max_Newton_iterations__subsequent;
+const bool my_active_flag = hs.has_genuine_horizons();
+const int N_horizons = hs.N_horizons();
+
+// buffers for  reduce_and_broadcast()
+int*            hn_buffer = NULL;
+fp *      rms_norm_buffer = NULL;
+fp * infinity_norm_buffer = NULL;
+if (hn_buffer == NULL)
+   then {
+	// allocate on first call
+	           hn_buffer = new int[N_active_procs];
+	     rms_norm_buffer = new fp [N_active_procs];
+	infinity_norm_buffer = new fp [N_active_procs];
+	}
+
+    //
+    // each pass through this loop finds a single horizon,
+    // or does corresponding dummy-horizon calls
+    //
+    // note there is no explicit exit test, rather we exit from the middle
+    // of the loop (only) when all processors are done with all their genuine
+    // horizons
+    //
+    for (int hn = hs.init_hn() ; ; hn = hs.next_hn())
+    {
+    if (verbose_info.print_algorithm_details)
+       then CCTK_VInfo(CCTK_THORNSTRING,
+		       "Newton_solve(): processor %d working on horizon %d",
+		       my_proc, hn);
 
-	for (int iteration = 1 ;
-	     iteration <= max_Newton_iterations ;
-	     ++iteration)
+    const bool horizon_is_genuine               = hs.is_genuine();
+    const bool there_is_another_genuine_horizon = hs.is_next_genuine();
+    if (verbose_info.print_algorithm_details)
+       then {
+	    CCTK_VInfo(CCTK_THORNSTRING,
+		       "                horizon_is_genuine=%d",
+		       int(horizon_is_genuine));
+	    CCTK_VInfo(CCTK_THORNSTRING,
+		       "                there_is_another_genuine_horizon=%d",
+		       int(there_is_another_genuine_horizon));
+	    }
+
+    struct AH_info* AH_info_ptr
+	                  = horizon_is_genuine ? my_AH_info[hn]      : NULL;
+    patch_system*  ps_ptr = horizon_is_genuine ? AH_info_ptr->ps_ptr : NULL;
+    Jacobian*     Jac_ptr = horizon_is_genuine ? AH_info_ptr->Jac_ptr: NULL;
+
+    const int max_iterations
+       = horizon_is_genuine
+	 ? (AH_info_ptr->initial_find_flag
+	    ? solver_info.max_Newton_iterations__initial
+	    : solver_info.max_Newton_iterations__subsequent)
+	 : INT_MAX;
+    
+	//
+	// each pass through this loop does a single Newton iteration
+	// on the current horizon (which might be either genuine or dummy)
+	//
+	for (int iteration = 1 ; ; ++iteration)
 	{
-	if (verbose_info.print_algorithm_details)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "Newton iteration %d/%d",
-			   iteration, max_Newton_iterations);
-	if (active_flag
+      #ifdef DEBUG
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "Newton_solve(): iteration %d",
+		   iteration);
+      #endif
+
+
+	//
+	// evaluate the expansion Theta(h) and its norms
+	// *** this is a synchronous operation across all processors ***
+	//
+	if (horizon_is_genuine
 	    && solver_info.debugging_output_at_each_Newton_iteration)
-	   then output_gridfn(ps, gfns::gfn__h,
+	   then output_gridfn(*ps_ptr, gfns::gfn__h,
 			      IO_info, IO_info.h_base_file_name,
 			      hn, verbose_info.print_algorithm_details,
 			      iteration);
 
-	//
-	// evaluate the expansion Theta(h)
-	//
+      #ifdef DEBUG
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   computing Theta: horizon_is_genuine=%d",
+		   int(horizon_is_genuine));
+      #endif
 	jtutil::norm<fp> Theta_norms;
-	if (! expansion(ps,
-			cgi, gi,
-			true,		// yes, we want the Jacobian coeffs
-			verbose_info.print_algorithm_details,
-			&Theta_norms))
-	   then return false;				// *** ERROR RETURN ***
-	if (verbose_info.print_algorithm_highlights)
+	const bool Theta_is_ok
+	   = expansion(ps_ptr,
+		       cgi, gi,
+		       true,	// yes, we want the Jacobian coeffs
+		       verbose_info.print_algorithm_details,
+		       &Theta_norms);
+      #ifdef DEBUG
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "Newton_solve(): Theta_is_ok=%d",
+		   Theta_is_ok);
+	if (horizon_is_genuine && Theta_is_ok)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
-			  "   iteration %d: Theta ||rms||=%.1e, ||infty||=%.1e",
-			   iteration,
-			   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
-	if (active_flag
+			   "   Theta rms-norm %.1e, infinity-norm %.1e",
+			   double(Theta_norms.rms_norm()),
+			   double(Theta_norms.infinity_norm()));
+      #endif
+	if (horizon_is_genuine && Theta_is_ok
 	    && solver_info.debugging_output_at_each_Newton_iteration)
-	   then output_gridfn(ps, gfns::gfn__Theta,
+	   then output_gridfn(*ps_ptr, gfns::gfn__Theta,
 			      IO_info, IO_info.Theta_base_file_name,
 			      hn, verbose_info.print_algorithm_details,
 			      iteration);
 
+
+	//
+	// compute and reduce-across-processors the are-we-done flags
+	// also broadcast and print error norms
+	//
+
+	// have we found *this* horizon?
+	const bool found_this_horizon
+	   = horizon_is_genuine && Theta_is_ok
+	     && (Theta_norms.infinity_norm() <= solver_info
+						.Theta_norm_for_convergence);
+	if (horizon_is_genuine)
+	   then AH_info_ptr->found_flag = found_this_horizon;
+	if (found_this_horizon)
+	   then {
+		compute_BH_diagnostics(*ps_ptr, BH_diagnostics_info,
+				       verbose_info,
+				       AH_info_ptr->BH_diagnostics);
+		// FIXME: see header comment
+		print_BH_diagnostics(AH_info_ptr->BH_diagnostics,
+				     N_horizons, hn,
+				     verbose_info);
+		if (IO_info.output_BH_diagnostics)
+		   then {
+			if (AH_info_ptr->BH_diagnostics_fileptr == NULL)
+			   then AH_info_ptr->BH_diagnostics_fileptr
+				  = setup_BH_diagnostics_output_file(IO_info,
+								     N_horizons,
+								     hn);
+			output_BH_diagnostics_fn(AH_info_ptr->BH_diagnostics,
+						 IO_info,
+						 hn, AH_info_ptr
+						     ->BH_diagnostics_fileptr);
+			}
+		}
+
+	// does *this* horizon need more iterations?
+	// i.e. has this horizon's Newton iteration not yet converged?
+	const bool this_horizon_needs_more_iterations
+	   = horizon_is_genuine && Theta_is_ok
+	     && !found_this_horizon
+	     && (iteration <= max_iterations);
+
+	// do we (this processor) need to do more iterations
+	// on this or a following horizon?
+	const bool we_need_more_iterations
+	   = this_horizon_needs_more_iterations
+	     || there_is_another_genuine_horizon;
+
+	if (verbose_info.print_algorithm_details)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "flags: found_this_horizon=%d",
+			   int(found_this_horizon));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "       this_horizon_needs_more_iterations=%d",
+			   int(this_horizon_needs_more_iterations));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "       we_need_more_iterations=%d",
+			   int(we_need_more_iterations));
+		}
+
+	// does any processor need more iterations
+	// on its current or a following horizon?
+	const bool any_proc_needs_more_iterations
+	   = reduce_and_broadcast
+		(GH, N_procs, N_active_procs,
+		 my_proc, my_active_flag,
+		 we_need_more_iterations,
+		 hn, horizon_is_genuine ? Theta_norms.rms_norm() : 0.0,
+		     horizon_is_genuine ? Theta_norms.infinity_norm() : 0.0,
+		 hn_buffer, rms_norm_buffer, infinity_norm_buffer);
+
+	if (verbose_info.print_algorithm_details)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "       ==> any_proc_needs_more_iterations=%d",
+			   int(any_proc_needs_more_iterations));
+		}
+	if (my_active_flag && verbose_info.print_algorithm_highlights)
+	   then print_Theta_norms(N_active_procs,
+				  hn_buffer, rms_norm_buffer,
+					     infinity_norm_buffer);
+
+
 	//
-	// convergence test on ||Theta||
+	// are all processors done with all their genuine horizons?
+	// or if this is a single-processor run, are we done with this horizon?
 	//
-	if (Theta_norms.infinity_norm() <= solver_info
-					   .Theta_norm_for_convergence)
+	if (! any_proc_needs_more_iterations)
 	   then {
 		if (verbose_info.print_algorithm_details)
 		   then CCTK_VInfo(CCTK_THORNSTRING,
-				   "==> finished (||Theta|| < %g)",
-				   double(solver_info
-					  .Theta_norm_for_convergence));
-		return true;				// *** NORMAL RETURN ***
+				   "==> all processors are finished!");
+		return;					// *** NORMAL RETURN ***
+		}
+	if ((N_procs == 1) && !this_horizon_needs_more_iterations)
+	   then {
+	      #ifdef DEBUG
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "   single processor: skipping to next horizon");
+	      #endif
+		break;					// *** LOOP EXIT ***
+		}
+
+
+	//
+	// compute the Jacobian matrix
+	// *** this is a synchronous operation across all processors ***
+	//
+      #ifdef DEBUG
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   computing Jacobian: genuine/dummy flag %d",
+		   int(this_horizon_needs_more_iterations));
+      #endif
+	const bool Jacobian_is_ok
+	   = expansion_Jacobian(this_horizon_needs_more_iterations ? ps_ptr
+								   : NULL,
+				this_horizon_needs_more_iterations ? Jac_ptr
+								   : NULL,
+				cgi, gi, Jacobian_info,
+				verbose_info.print_algorithm_details);
+
+
+	//
+	// skip to the next horizon unless
+	// this is a genuine Jacobian computation, and it went ok
+	//
+	if (! (this_horizon_needs_more_iterations && Jacobian_is_ok))
+	   then {
+	      #ifdef DEBUG
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "   skipping to next horizon");
+	      #endif
+		break;					// *** LOOP EXIT ***
 		}
 
+
 	//
 	// compute the Newton step
 	//
-	if (! expansion_Jacobian(ps, Jac,
-				 cgi, gi, Jacobian_info,
-				 verbose_info.print_algorithm_details))
-	   then return false;				// *** ERROR RETURN ***
 	if (verbose_info.print_algorithm_details)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "solving linear system for Delta_h (%d unknowns)",
-			   Jac.N_rows());
-	const fp rcond = Jac.solve_linear_system(gfns::gfn__Theta,
-						 gfns::gfn__Delta_h,
-						 verbose_info
-						 .print_algorithm_details);
+			   Jac_ptr->N_rows());
+	const fp rcond
+	   = Jac_ptr->solve_linear_system(gfns::gfn__Theta,
+					  gfns::gfn__Delta_h,
+					  verbose_info.print_algorithm_details);
 	if	(rcond == 0.0)
 	   then {
 		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "Newton_solve: Jacobian matrix is numerically singular!");
-		return false;				// *** ERROR RETURN ***
+		// give up on this horizon
+		break;				// *** LOOP CONTROL ***
 		}
 	if (verbose_info.print_algorithm_details)
 	   then {
@@ -156,96 +371,242 @@ const int max_Newton_iterations
 				   "done solving linear system");
 		}
 
-	if (active_flag
-	    && solver_info.debugging_output_at_each_Newton_iteration)
-	   then output_gridfn(ps, gfns::gfn__Delta_h,
+	if (solver_info.debugging_output_at_each_Newton_iteration)
+	   then output_gridfn(*ps_ptr, gfns::gfn__Delta_h,
 			      IO_info, IO_info.Delta_h_base_file_name,
 			      hn, verbose_info.print_algorithm_details,
 			      iteration);
 
-	//
-	// if the Newton step is too large, scale it down
-	//
-	jtutil::norm<fp> h_norms;
-	ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
-	const fp max_allowable_Delta_h
-		= solver_info.max_Delta_h_over_h * h_norms.mean();
-	jtutil::norm<fp> Delta_h_norms;
-	ps.gridfn_norms(gfns::gfn__Delta_h, Delta_h_norms);
-	const fp max_Delta_h = Delta_h_norms.infinity_norm();
-	const fp scale = (max_Delta_h > max_allowable_Delta_h)
-			 ? max_allowable_Delta_h / max_Delta_h
-			 : 1.0;
 
 	//
 	// take the Newton step (scaled if need be)
 	//
-		for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
-		{
-		patch& p = ps.ith_patch(pn);
-
-			for (int irho = p.min_irho() ;
-			     irho <= p.max_irho() ;
-			     ++irho)
-			{
-			for (int isigma = p.min_isigma() ;
-			     isigma <= p.max_isigma() ;
-			     ++isigma)
-			{
-			p.ghosted_gridfn(gfns::gfn__h, irho,isigma)
-			   -= scale * p.gridfn(gfns::gfn__Delta_h, irho,isigma);
-			}
-			}
-		}
-	if (verbose_info.print_algorithm_details)
-	   then {
-		if (scale == 1.0)
-		   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "h += Delta_h (rms-norm=%.1e, infinity-norm=%.1e)",
-				   Delta_h_norms.rms_norm(),
-				   Delta_h_norms.infinity_norm());
-		   else CCTK_VInfo(CCTK_THORNSTRING,
-			   "h += %g * Delta_h (infinity-norm clamped to %.2g)",
-				   scale,
-				   scale * Delta_h_norms.infinity_norm());
-		}
+	Newton_step(*ps_ptr,
+		    solver_info.max_allowable_Delta_h_over_h,
+		    verbose_info);
 
-	//
-	// convergence test on ||Delta_h||
-	//
-	if ( scale * Delta_h_norms.infinity_norm()
-	     <= solver_info.Delta_h_norm_for_convergence )
-	   then {
-		if (verbose_info.print_algorithm_details)
-		   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "==> finished (||Delta_h|| < %g)",
-				   double(solver_info
-					  .Delta_h_norm_for_convergence));
-		if (solver_info.final_Theta_update_if_Delta_h_converged)
-		   then {
-			if (verbose_info.print_algorithm_details)
-			   then CCTK_VInfo(CCTK_THORNSTRING,
-					 "    doing final Theta(h) evaluation");
-			if (! expansion(ps,
-					cgi, gi,
-					false,		// no Jacobian coeffs
-					verbose_info.print_algorithm_details))
-			   then return false;		// *** ERROR RETURN ***
-			}
-		return true;				// *** NORMAL RETURN ***
-		}
+	// end of this Newton iteration
 	}
 
-if (solver_info.final_Theta_update_if_Delta_h_converged)
+    // end of this horizon
+    }
+
+// we should never get to here
+assert( false );
+}
+
+//******************************************************************************
+
+//
+// This function does an inclusive-or--reduction of the "we need more
+// iterations (either on this or on another genuine horizon" flags across
+// all active processors.  It also broadcasts which horizon we're working
+// on, and the current error norms, to processor #0.
+//
+// The reduction and broadcast are combined in this function to allow
+// doing them all with a single reduction operation (and thus only a
+// single set of MPI overheads).
+//
+// Arguments:
+// GH --> The Cactus grid hierarchy.
+// N_procs = The total number of processors.
+// N_active_procs = The number of active processors.
+// my_active_flag = Is this processor an active processor?
+// we_need_more_iterations = "we need more iterations (either on this
+//			      or on another genuine horizon" flag to be
+//			     inclusive-or--reduced
+// hn,{rms,infinity}_norm = horizon number and norms: if this is an active
+//			    processor these are broadcast to processor #0,
+//			    otherwise these are ignored
+// {hn,{rms,infinity}_norm}_buffer[N_active_procs]
+//	= buffers for horizon numbers and norms: if this is processor #0
+//	  these are set to the horizon numbers and norms for all active
+//	  processors, otherwise these are set to undefined (garbage) values
+//
+// Results:
+// This function returns the inclusive-or--reduction of the
+//  we_need_more_iterations  flags across all active processors.
+//
+namespace {
+bool reduce_and_broadcast(cGH* GH,
+			  int N_procs, int N_active_procs,
+			  int my_proc, bool my_active_flag,
+			  bool we_need_more_iterations,
+			  int hn, fp rms_norm, fp infinity_norm,
+			  int hn_buffer[],
+			  fp rms_norm_buffer[],
+			  fp infinity_norm_buffer[])
+{
+assert( my_proc >= 0 );
+assert( my_proc < N_procs );
+//
+// We do the combined reduce/broadcast in a single reduction operation
+// by setting the buffers to all 0s, except that on each processor the
+// [my_proc] entries have the values we want to reduce/broadcast, then
+// doing a sum-reduction of the buffers.
+//
+// Alas, to do everything in a single operation, all the values have to
+// be of a single datatype in a single array, so we convert  hn  to
+// CCTK_REAL, and encode the "we need more iterations" in its sign.
+//
+const int N_reduce_buffer = 3*N_active_procs;
+static CCTK_REAL* reduce_buffer = NULL;
+if (reduce_buffer == NULL)
    then {
-	if (verbose_info.print_algorithm_details)
+	// allocate on first call
+	reduce_buffer = new CCTK_REAL[N_reduce_buffer];
+	}
+
+//
+// pack the values into the reduction buffer
+//
+	for (int i = 0 ; i < N_reduce_buffer ; ++i)
+	{
+	reduce_buffer[i] = 0.0;
+	}
+if (my_active_flag)
+   then {
+	assert( my_proc < N_active_procs );
+	reduce_buffer[3*my_proc + 0] = we_need_more_iterations ? hn : -hn;
+	reduce_buffer[3*my_proc + 1] = rms_norm;
+	reduce_buffer[3*my_proc + 2] = infinity_norm;
+	}
+
+//
+// do the reduction
+//
+
+// this name is appropriate for PUGHReduce, caveat user for other drivers :)
+const int reduction_handle = CCTK_ReductionArrayHandle("sum");
+if (reduction_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "reduce_and_broadcast(): can't get sum-reduction handle!");
+								/*NOTREACHED*/
+
+const int status
+   = CCTK_ReduceLocArrayToArray1D(GH,
+				  -1,		// all processors get result
+				  reduction_handle,
+				  (/*const*/ void*) reduce_buffer,
+				  (          void*) reduce_buffer,
+				  N_reduce_buffer,
+				  CCTK_VARIABLE_REAL);
+if (status < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "reduce_and_broadcast(): error status %d from reduction!",
+		   status);					/*NOTREACHED*/
+
+//
+// unpack the reduction buffer back to our results
+//
+bool any_proc_needs_more_iterations = false;
+	for (int proc = 0 ; proc < N_active_procs ; ++proc)
+	{
+	const int i_proc = 3*proc;
+	const bool proc_needs_more_iterations = reduce_buffer[i_proc + 0] > 0.0;
+	any_proc_needs_more_iterations |= proc_needs_more_iterations;
+	           hn_buffer[proc] = (int) fabs(reduce_buffer[i_proc + 0]);
+	     rms_norm_buffer[proc] =            reduce_buffer[i_proc + 1];
+	infinity_norm_buffer[proc] =            reduce_buffer[i_proc + 2];
+	}
+
+return any_proc_needs_more_iterations;
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function is called on processor #0 to print the Theta norms.
+//
+// Arguments:
+// N_active_procs = The number of active processors.
+// hn_buffer[N_active_procs] = Gives the horizon each active processor
+//			       is currently working on.
+// {rms,infinity}_norm_buffer[N_active_procs] = Gives the Theta norms for
+//						each active processor.
+//
+namespace {
+void print_Theta_norms(int N_active_procs,
+		       const int hn_buffer[],
+		       const fp rms_norm_buffer[],
+		       const fp infinity_norm_buffer[])
+{
+	for (int proc = 0 ; proc < N_active_procs ; ++proc)
+	{
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   proc %d/horizon %d: ||Theta|| rms=%.1e, infty=%.1e",
+		   proc, int(hn_buffer[proc]),
+		   double(rms_norm_buffer[proc]),
+		   double(infinity_norm_buffer[proc]));
+	}
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function takes the Newton step, scaling it down if it's too large.
+//
+// Arguments:
+// ps = The patch system containing the gridfns h and Delta_h.
+// max_allowable_Delta_h_over_h = The maximum allowable
+//				     ||Delta_h||_infinity / ||h||_mean
+//				  Any step over this is internally clamped
+//				  (scaled down) to this size.
+//
+namespace {
+void Newton_step(patch_system& ps,
+		 fp max_allowable_Delta_h_over_h,
+		 const struct verbose_info& verbose_info)
+{
+//
+// compute scale factor (1 for small steps, <1 for large steps)
+//
+
+jtutil::norm<fp> h_norms;
+ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
+const fp max_allowable_Delta_h = max_allowable_Delta_h_over_h * h_norms.mean();
+
+jtutil::norm<fp> Delta_h_norms;
+ps.gridfn_norms(gfns::gfn__Delta_h, Delta_h_norms);
+const fp max_Delta_h = Delta_h_norms.infinity_norm();
+
+const fp scale = (max_Delta_h <= max_allowable_Delta_h)
+		 ? 1.0
+		 : max_allowable_Delta_h / max_Delta_h;
+
+if (verbose_info.print_algorithm_details)
+   then {
+	if (scale == 1.0)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "    doing final Theta(h) evaluation");
-	if (! expansion(ps,
-			cgi, gi,
-			false,		// no Jacobian coeffs
-			verbose_info.print_algorithm_details))
-	   then return false;				// *** ERROR RETURN ***
+			   "h += Delta_h (rms-norm=%.1e, infinity-norm=%.1e)",
+			   Delta_h_norms.rms_norm(),
+			   Delta_h_norms.infinity_norm());
+	   else CCTK_VInfo(CCTK_THORNSTRING,
+			   "h += %g * Delta_h (infinity-norm clamped to %.2g)",
+			   scale,
+			   scale * Delta_h_norms.infinity_norm());
+	}
+
+
+//
+// take the Newton step (scaled if necessary)
+//
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
+		{
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
+		     ++isigma)
+		{
+		p.ghosted_gridfn(gfns::gfn__h, irho,isigma)
+			-= scale * p.gridfn(gfns::gfn__Delta_h, irho,isigma);
+		}
+		}
 	}
-return false;						// *** ERROR RETURN ***
 }
+	  }
diff --git a/src/driver/README b/src/driver/README
index 053da7a..7fac708 100644
--- a/src/driver/README
+++ b/src/driver/README
@@ -3,24 +3,28 @@ and interface with the rest of Cactus.
 
 The main source files in this directory are as follows:
 
+README.parallel
+	design notes for the parallel/multiprocessor algorithms
+	and data structures
+
 driver.hh
 	overall header file for routines in this directory
 
 state.cc
 	state information (data structures) kept across Cactus scheduler calls
 
-setup.cc
+setup.cc		# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
 	this is called from the scheduler to setup our data structures
 
-find_horizons.cc
+find_horizons.cc	# sees CCTK_ARGUMENTS, CCTK_PARAMETERS
 	this is called from the scheduler to find the apparent horizon(s)
 	in a slice
 
 initial_guess.cc
 	this sets up the initial guess(es) for the horizon position(s)
 
-Newton.cc
-	this solves the H(h) = 0 equations by Newton's method
+Newton.cc		# sees cGH (for interprocessor broadcasts)
+	this solves the Theta(h) = 0 equations by Newton's method
 
 io.cc
 	I/O routines
diff --git a/src/driver/README.parallel b/src/driver/README.parallel
new file mode 100644
index 0000000..3d553e0
--- /dev/null
+++ b/src/driver/README.parallel
@@ -0,0 +1,145 @@
+Design notes for the parallel/multiprocessor algorithms and data structures
+===========================================================================
+
+$Header$
+
+The Problem
+-----------
+Doing horizon finding on multiple processors is difficult, because when
+we interpolate the geometry to the trial horizon surface, we must use a
+global interpolator (CCTK_InterpGridArrays()), which is (necessarily) a
+a collective operation, i.e. the interpolator must be called synchronously
+on every processor, with identical arguments except for the number of
+interpolation points, the interpolation coordinates, and the output
+array pointers.  It's quite tricky to arrange this synchronicity when
+different processors are locating different horizons.
+
+
+Newton's Method
+---------------
+We find the apparent horizon by solving   Theta(h) = 0   by a global
+Newton's-method iteration.  ("global" meaning we update all the horizon
+surface points simultaneously at each iteration.)  That is, the basic
+uniprocessor solution algorithm is
+
+	while (true)
+	{
+	compute Theta(h)
+	if (||Theta(h)|| small or too many iterations)
+		exit
+	compute Jacobian[Theta(h)]
+	solve linear system J.dh = -Theta
+	update h += dh
+	}
+
+As mentioned above, the computation of Theta(h) must be done synchronously
+(across all processors).  It turns out that the evaluation of the Jacobian
+matrix must also be done synchronously -- it too involves a geometry
+interpolation.
+
+
+High-Level Outline of the Solution
+----------------------------------
+We use the following design to do this:
+
+Horizons are identified by an integer "horizon number".  This runs from
+1 to N_horizons for the "genuine" horizons we want to find, or 0 for a
+"dummy" horizon used if there are more processors than horizons to find.
+	[We actually implement the dummy horizon as
+	passing NULL horizon and/or Jacobian pointers
+	to expansion() and expansion_Jacobian(), but
+	this is an implementation detail which doesn't
+	matter here.]
+
+Each processor is assigned the dummy horizon, and zero or more genuine
+horizons; each genuine horizon is assigned to exactly one procesor.
+A processor with one or more genuine horizons is called an "active"
+processor; a processor with only the dummy horizon is called a "dummy"
+processor.  The assignment of horizons to processors is static for the
+whole Cactus run.
+	[A dynamic assignment might be more efficient,
+	but then we'd have to move our state for a horizon
+	from one processor to processor.]
+
+All the processors do their Newton iteration synchronously, working on
+genuine horizons if they have any, or dummy horizons if not.
+
+For example, suppose we have 3 horizons, which are found (the Newton
+iteration converges) after respectively 3, 5, and 4 iterations.  Then
+with 2 processors the iterations would look like this (where h1/2/3
+means horizon 1/2/3, h- means the dummy horizon, and a * after Theta
+means convergence):
+
+	processor #0	processor #1
+	------------	------------
+1	h1 Theta	h2 Theta
+2	h1 Jacobian	h2 Jacobian
+3	h1 Theta	h2 Theta
+4	h1 Jacobian	h2 Jacobian
+5	h1 Theta*	h2 Theta
+6	h- Jacobian	h2 Jacobian
+7	h3 Theta	h2 Theta
+8	h3 Jacobian	h2 Jacobian
+9	h3 Theta	h2 Theta*
+10	h3 Jacobian	h- Jacobian
+11	h3 Theta	h- Theta
+12	h3 Jacobian	h- Jacobian
+13	h3 Theta*	h- Theta
+
+(Notice that at line 6, processor #0 does a dummy-horizon Jacobian
+computation before starting its next genuine horizon.  This is to keep
+the Theta and Jacobian computations synchronized across all processors.
+In a single-processor run we'd skip this and go directly to the next
+genuine horizon.)
+
+With 3 processors this same example would look like this:
+
+	processor #0	processor #1	processor #2
+	------------	------------	------------
+1	h1 Theta	h2 Theta	h3 Theta
+2	h1 Jacobian	h2 Jacobian	h3 Jacobian
+3	h1 Theta	h2 Theta	h3 Theta
+4	h1 Jacobian	h2 Jacobian	h3 Jacobian
+5	h1 Theta*	h2 Theta	h3 Theta
+6	h- Jacobian	h2 Jacobian	h3 Jacobian
+7	h- Theta	h2 Theta	h3 Theta*
+8	h- Jacobian	h2 Jacobian	h- Jacobian
+9	h- Theta	h2 Theta*	h- Theta
+
+With 4 processors it would look like this:
+
+	processor #0	processor #1	processor #2	processor #3
+	------------	------------	------------	------------
+1	h1 Theta	h2 Theta	h3 Theta	h- Theta
+2	h1 Jacobian	h2 Jacobian	h3 Jacobian	h- Jacobian
+3	h1 Theta	h2 Theta	h3 Theta	h- Theta
+4	h1 Jacobian	h2 Jacobian	h3 Jacobian	h- Jacobian
+5	h1 Theta*	h2 Theta	h3 Theta	h- Theta
+6	h- Jacobian	h2 Jacobian	h3 Jacobian	h- Jacobian
+7	h- Theta	h2 Theta	h3 Theta*	h- Theta
+8	h- Jacobian	h2 Jacobian	h- Jacobian	h- Jacobian
+9	h- Theta	h2 Theta*	h- Theta	h- Theta
+
+Any additional processors would similarly do all Theta and Jacobian
+computations on the dummy horizon.
+
+
+How and What to Synchronize
+---------------------------
+To implement this algorithm, after each Theta evaluation each active
+processor computes a "we need more iterations (either on this or another
+genuine horizon)" flag, and all the processors do an inclusive-or-reduction
+of these flags, with the results broadcast to all processors.
+
+(If any error occurs when computing Theta or the Jacobian, or in solving
+the linear equations, we treat this as failing to find this horizon.)
+
+Each processor then uses the inclusive-or-reduction result to determine
+whether or not to continue iterating.
+
+After each Theta evaluation each active processor also broadcasts
+* which horizon it's working on
+* the iteration number
+* its error norms
+to processor 0, which uses this information to print a status report
+of how the iterations are converging.
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index bbee5aa..e5ad9aa 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -14,9 +14,9 @@
 //
 enum	method
 	{
-	method__evaluate_expansion,
-	method__test_expansion_Jacobian,
-	method__find_horizon // no comma
+	method__evaluate_expansions,
+	method__test_expansion_Jacobians,
+	method__find_horizons // no comma
 	};
 
 //
@@ -88,10 +88,8 @@ struct	solver_info
 	bool debugging_output_at_each_Newton_iteration;
 	int max_Newton_iterations__initial,
 	    max_Newton_iterations__subsequent;
-	fp max_Delta_h_over_h;
+	fp max_allowable_Delta_h_over_h;
 	fp Theta_norm_for_convergence;
-	fp Delta_h_norm_for_convergence;
-	bool final_Theta_update_if_Delta_h_converged;
 	};
 
 //
@@ -172,6 +170,9 @@ struct	BH_diagnostics
 	fp circumference_xy, circumference_xz, circumference_yz;
 	fp area;
 	fp m_irreducible;
+
+	// constructor initializes all diagnostics to 0.0
+	BH_diagnostics();
 	};
 
 //
@@ -183,38 +184,59 @@ struct	AH_info
 	patch_system* ps_ptr;
 	Jacobian* Jac_ptr;
 
+	// are we finding this horizon "from scratch"?
+	// ... true if this is the first time we've tried to find it,
+	//          or if we've tried before but failed to find it the
+	//          last time
+	//     false if we've tried before and succeeded the last time
+	//           (so we have that position as a very good initial guess)
+	bool initial_find_flag;
+
 	struct initial_guess_info initial_guess_info;
 
-	bool AH_found;
+	bool found_flag;	// did we find this horizon (successfully)
 	struct BH_diagnostics BH_diagnostics;
 	FILE *BH_diagnostics_fileptr;
 	};
 
 //
-// (A single copy of) this struct holds all of our state that's
-// persistent across Cactus scheduler calls.  This copy lives in "state.cc".
+// (A single copy of) this struct holds all of our state on this processor
+// that's persistent across Cactus scheduler calls.  This copy lives in
+// "state.cc".
 //
 struct	state
 	{
 	enum method method;
 	struct verbose_info verbose_info;
 	int timer_handle;
-	int my_proc;			// our processor number (0 to N_procs-1)
 	int N_procs;			// total number of processors
+	int my_proc;			// processor number of this processor
+					// (0 to N_procs-1)
+
+	int N_horizons;			// total number of genuine horizons
+					// being searched for
+	int N_active_procs;		// total number of active processors
+					// (the active processors are processor
+					//  numbers 0 to N_active_procs-1)
 
-	struct IO_info IO_info;
-	struct Jacobian_info Jac_info;
-	struct solver_info solver_info;
 	struct cactus_grid_info cgi;
 	struct geometry_info gi;
+	struct Jacobian_info Jac_info;
+	struct solver_info solver_info;
+	struct IO_info IO_info;
 
 	struct BH_diagnostics_info BH_diagnostics_info;
 
-	int N_horizons;
+	horizon_sequence *my_hs;	// --> new-allocated object describing
+					//     the sequence of horizons
+					//     assigned to this processor
 
-	// --> array of size N_horizons+1,
-	// indexed with "horizon number" hn (should be in range [1,N_horizons]
-	AH_info* AH_info_array;
+	// horizon numbers ("hn") run from 1 to N_horizons inclusive
+	struct AH_info** my_AH_info;	// --> new[]-allocated array of size
+					//     N_horizons+1, subscripted by hn,
+					//     of --> info or NULL for horizons
+					//            not assigned to this
+					//            processor
 	};
 
 //******************************************************************************
@@ -244,16 +266,31 @@ enum initial_guess_method
 
 // Newton.cc
 // returns true for success, false for failure to converge
-bool Newton_solve(patch_system& ps,
-		  Jacobian& Jac,
-		  const struct cactus_grid_info& cgi,
-		  const struct geometry_info& gi,
-		  const struct Jacobian_info& Jacobian_info,
-		  const struct solver_info& solver_info,
-		  bool initial_find_flag,
-		  const struct IO_info& IO_info,
-		  bool active_flag,
-		  int hn, const struct verbose_info& verbose_info);
+void Newton(cGH* GH,
+	    int N_procs, int N_active_procs, int my_proc,
+	    horizon_sequence& hs, struct AH_info* const my_AH_info[],
+	    const struct cactus_grid_info& cgi,
+	    const struct geometry_info& gi,
+	    const struct Jacobian_info& Jacobian_info,
+	    const struct solver_info& solver_info,
+	    const struct IO_info& IO_info,
+	    const struct BH_diagnostics_info& BH_diagnostics_info,
+	    const struct verbose_info& verbose_info);
+
+// BH_diagnostics.cc
+void compute_BH_diagnostics
+	(const patch_system& ps,
+	 const struct BH_diagnostics_info& BH_diagnostics_info,
+	 const struct verbose_info& verbose_info,
+	 struct BH_diagnostics& BH_diagnostics);
+void print_BH_diagnostics(const struct BH_diagnostics& BH_diagnostics,
+			  int N_horizons, int hn,
+			  const struct verbose_info& verbose_info);
+FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
+				       int N_horizons, int hn);
+void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
+			      const struct IO_info& IO_info,
+			      int hn, FILE* fileptr);
 
 // io.cc
 void input_gridfn(patch_system& ps, int unknown_gfn,
@@ -267,10 +304,3 @@ void output_Jacobians(const patch_system& ps,
 		      const Jacobian* Jac_SD_FDdr_ptr,
 		      const struct IO_info& IO_info, const char base_file_name[],
 		      int hn, bool print_msg_flag, int AHF_iteration = 0);
-FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
-				       int hn, int N_horizons);
-void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
-			      const struct IO_info& IO_info,
-			      int hn, FILE* fileptr);
-enum horizon_file_format
-  decode_horizon_file_format(const char horizon_file_format_string[]);
diff --git a/src/driver/ellipsoid.log b/src/driver/ellipsoid.log
index 71cad20..e69de29 100644
--- a/src/driver/ellipsoid.log
+++ b/src/driver/ellipsoid.log
@@ -1,38 +0,0 @@
-    |\^/|     Maple 7 (IBM INTEL LINUX)
-._|\|   |/|_. Copyright (c) 2001 by Waterloo Maple Inc.
- \  MAPLE  /  All rights reserved. Maple is a registered trademark of
- <____ ____>  Waterloo Maple Inc.
-      |       Type ? for help.
-# ellipsoid.maple -- compute equations for offset ellipsoid setup
-# $Header$
-> 
-#
-# ellipsoid has center (A,B,C), radius (a,b,c)
-# angular coordinate system has center (U,V,W)
-#
-# direction cosines wrt angular coordinate center are (alpha,beta,gamma)
-# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
-# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
-#
-# then the equation of the ellipsoid is
-#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
-#	-----------------  +  ----------------  +  -----------------  =  1
-#	        a^2                  b^2                   c^2
-#
-# to solve this, we introduce intermediate variables
-#	AU = A - U
-#	BV = B - V
-#	CW = C - W
-#
-> eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
-                              2                2                2
-                 (xcos r - AU)    (ycos r - BV)    (zcos r - CW)
-          eqn := -------------- + -------------- + -------------- = 1
-                        2                2                2
-                       a                b                c
-
-> 
-> read "../maple/util.mm";
-Error, unable to read `../maple/util.mm`
-> quit
-bytes used=129844, alloc=196572, time=0.03
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 01e28e2..fd35afd 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -11,11 +11,7 @@
 /// find_horizon - find a horizon
 /// do_evaluate_expansion
 /// do_test_expansion_Jacobian
-/// do_find_horizon
 ///
-/// compute_BH_diagnostics - compute BH diagnostics for a horizon
-/// surface_integral - compute surface integral of a gridfn over the horizon
-//
 
 #include <stdio.h>
 #include <assert.h>
@@ -48,6 +44,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
 //******************************************************************************
@@ -68,28 +65,25 @@ const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
 					const char gridfn_name[],
 					bool check_for_NULL = true);
 
-bool do_evaluate_expansion
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps,
-	 bool active_flag, int hn, int N_horizons);
-bool do_test_expansion_Jacobian
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info,
-	 bool test_all_Jacobian_compute_methods,
-	 struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons);
-bool do_find_horizon
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 const struct solver_info& solver_info, bool initial_find_flag,
-	 const struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons);
+void do_evaluate_expansions(int my_proc, int N_horizons,
+			    horizon_sequence& hs,
+			       struct AH_info* const my_AH_info[],
+			    const struct cactus_grid_info& cgi,
+			    const struct geometry_info& gi,
+			    const struct IO_info& IO_info,
+			    bool output_h, bool output_Theta,
+			    const struct verbose_info& verbose_info,
+			    int timer_handle);
+void do_test_expansion_Jacobians(int my_proc, int N_horizons,
+				 struct AH_info* const my_AH_info[],
+				 const struct cactus_grid_info& cgi,
+				 const struct geometry_info& gi,
+				       struct Jacobian_info& Jac_info,
+				 bool test_all_Jacobian_compute_methods,
+				 const struct IO_info& IO_info,
+				 const struct verbose_info& verbose_info,
+				 int timer_handle);
+
 
 void compute_BH_diagnostics
 	(const patch_system& ps,
@@ -112,27 +106,30 @@ extern "C"
 DECLARE_CCTK_ARGUMENTS
 DECLARE_CCTK_PARAMETERS
 
-const struct verbose_info& verbose_info = state.verbose_info;
-      struct IO_info&      IO_info      = state.IO_info;
-const struct solver_info&  solver_info  = state.solver_info;
-
 if (state.timer_handle >= 0)
    then CCTK_TimerResetI(state.timer_handle);
 
+const int my_proc = state.my_proc;
+horizon_sequence& hs = *state.my_hs;
+const bool active_flag = hs.has_genuine_horizons();
+
+      struct cactus_grid_info&          cgi =  state.cgi;
+const struct    geometry_info&           gi =  state.gi;
+      struct    Jacobian_info&     Jac_info =  state.Jac_info;
+      struct          IO_info&      IO_info =  state.IO_info;
+const struct     verbose_info& verbose_info =  state.verbose_info;
+
 // what are the semantics of the Cactus gxx variables? (these may
 // change from one call to another, so we have to re-check each time)
 if      (CCTK_Equals(metric_type, "physical"))
-   then state.cgi.Cactus_conformal_metric = false;
+   then cgi.Cactus_conformal_metric = false;
 else if (CCTK_Equals(metric_type, "static conformal"))
-   then state.cgi.Cactus_conformal_metric = (conformal_state > 0);
+   then cgi.Cactus_conformal_metric = (conformal_state > 0);
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "AHFinderDirect_find_horizons(): unknown metric_type=\"%s\"!",
 		   metric_type);				/*NOTREACHED*/
 
-// only processor #0 is "active"
-const bool active_flag = (state.my_proc == 0);
-
-// get the Cactus time step and decide if we want to output [hH] now
+// get the Cactus time step and decide if we want to output h and/or Theta now
 IO_info.time_iteration = cctk_iteration;
 IO_info.time           = cctk_time;
 const bool output_h
@@ -142,32 +139,22 @@ const bool output_Theta
    = (IO_info.how_often_to_output_Theta > 0)
      && ((IO_info.time_iteration % IO_info.how_often_to_output_Theta) == 0);
 
-//
 // if we're using them,
 //    we need to re-fetch the Cactus data pointers at least each time step,
 //    because they change each time Cactus rotates the time levels
-//
-if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid)
-   then setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
+if (gi.geometry_method == geometry__local_interp_from_Cactus_grid)
+   then setup_Cactus_gridfn_data_ptrs(cctkGH, cgi);
 
-	for (int hn = 1 ; hn <= state.N_horizons ; ++hn)
+// set initial guess for any (genuine) horizons that need it,
+// i.e. for any (genuine) horizons where we didn't find the horizon previously
+	for (int hn = hs.init_hn() ; hs.is_genuine() ; hn = hs.next_hn())
 	{
-	struct AH_info& AH_info = state.AH_info_array[hn];
-	patch_system& ps = *AH_info.ps_ptr;
-
-	//
-	// If this is our first attempt to find this horizon, or
-	//    if we've tried to find it before but we failed on our
-	//    immediately previous attempt, then we need to (re)set
-	//    the initial guess.
-	// Otherwise (i.e. if we've tried to find this horizon before,
-	//    and we succeeded on our immediately previous attempt),
-	//    then we can just reuse the previous horizon position as
-	//    our initial guess for this time around.
-	//
-	const bool initial_find_flag = ! AH_info.AH_found;
-	if (initial_find_flag)
-	   then {
+	assert( state.my_AH_info[hn] != NULL );
+	struct AH_info& AH_info = *state.my_AH_info[hn];
+	if (AH_info.found_flag)
+	   then AH_info.initial_find_flag = false;
+	   else {
+		patch_system& ps = *AH_info.ps_ptr;
 		setup_initial_guess(ps,
 				    AH_info.initial_guess_info,
 				    IO_info,
@@ -177,129 +164,51 @@ if (state.gi.geometry_method == geometry__local_interp_from_Cactus_grid)
 				      IO_info, IO_info.h_base_file_name,
 				      hn, verbose_info
 					  .print_algorithm_highlights);
+		AH_info.initial_find_flag = true;
 		}
+	}
 
-	//
-	// create the Jacobian matrix if we're going to need it
-	// and it doesn't already exist
-	//
-	switch	(state.method)
-		{
-	case method__evaluate_expansion:
-		// no Jacobian needed ==> no-op here
-		break;
-	case method__test_expansion_Jacobian:
-	case method__find_horizon:
-		if (AH_info.Jac_ptr == NULL)
-		   then AH_info.Jac_ptr
-			   = new_Jacobian(state.Jac_info
-					       .Jacobian_store_solve_method,
-					  ps,
-					  verbose_info
-					  .print_algorithm_highlights);
-		break;
-	default:
-		CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "\n"
-			   "   find_horizons(): unknown method=(int)%d!\n"
-			   "                    in create-Jacobian switch\n"
-			   "                    (this should never happen!)"
-			   ,
-			   int(state.method));			/*NOTREACHED*/
-		}
+//
+// now the main horizon finding (or other computation)
+//
+switch	(state.method)
+	{
+case method__evaluate_expansions:
+	do_evaluate_expansions(my_proc, state.N_horizons,
+			       *state.my_hs, state.my_AH_info,
+			       cgi, gi,
+			       IO_info, output_h, output_Theta,
+			       state.verbose_info, state.timer_handle);
+	break;
 
-	AH_info.AH_found = false;
-	switch	(state.method)
-		{
-	case method__evaluate_expansion:
-		do_evaluate_expansion(verbose_info, state.timer_handle,
-				      IO_info, output_h, output_Theta,
-				      state.cgi, state.gi,
-				      ps,
-				      active_flag, hn, state.N_horizons);
-		break;
-
-	case method__test_expansion_Jacobian:
-		do_test_expansion_Jacobian(verbose_info, state.timer_handle,
-					   IO_info,
-				       (test_all_Jacobian_compute_methods != 0),
-					   state.Jac_info, state.cgi, state.gi,
-					   ps, *AH_info.Jac_ptr,
-					   active_flag, hn, state.N_horizons);
-		break;
-
-	case method__find_horizon:
-		AH_info.AH_found
-		   = do_find_horizon(verbose_info, state.timer_handle,
-				     IO_info, output_h, output_Theta,
-				     solver_info, initial_find_flag,
-				     state.Jac_info, state.cgi, state.gi,
-				     ps, *AH_info.Jac_ptr,
-				     active_flag, hn, state.N_horizons);
-
-		// store found flag in Cactus array variable
-		AH_found[hn] = AH_info.AH_found;
-
-		if (AH_info.AH_found)
-		   then {
-// compute BH diagnostics
-compute_BH_diagnostics(ps,
-		       state.BH_diagnostics_info,
-		       verbose_info, AH_info.BH_diagnostics);
-const struct BH_diagnostics& BH_diagnostics = AH_info.BH_diagnostics;
-
-			// print BH diagnostics?
-if (verbose_info.print_physics_details)
-   then {
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "AH %d/%d: r=%g at (%f,%f,%f)",
-		   hn, state.N_horizons,
-		   double(BH_diagnostics.mean_radius),
-		   double(BH_diagnostics.centroid_x),
-		   double(BH_diagnostics.centroid_y),
-		   double(BH_diagnostics.centroid_z));
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "AH %d/%d: area=%.10g m_irreducible=%.10g",
-		   hn, state.N_horizons,
-		   double(BH_diagnostics.area),
-		   double(BH_diagnostics.m_irreducible));
-	}
+case method__test_expansion_Jacobians:
+	do_test_expansion_Jacobians(my_proc, state.N_horizons,
+				    state.my_AH_info,
+				    cgi, gi, Jac_info, 
+				    (test_all_Jacobian_compute_methods != 0),
+				    IO_info, verbose_info, state.timer_handle);
+	break;
 
-// store BH diagnostics in Cactus array variables
-centroid_x[hn]    = BH_diagnostics.centroid_x;
-centroid_y[hn]    = BH_diagnostics.centroid_y;
-centroid_z[hn]    = BH_diagnostics.centroid_z;
-area[hn]          = BH_diagnostics.area;
-m_irreducible[hn] = BH_diagnostics.m_irreducible;
+case method__find_horizons:
+	if (state.timer_handle >= 0)
+	   then CCTK_TimerStartI(state.timer_handle);
+	Newton(cctkGH,
+	       state.N_procs, state.N_active_procs, my_proc,
+	       *state.my_hs, state.my_AH_info,
+	       cgi, gi, Jac_info, state.solver_info,
+	       IO_info, state.BH_diagnostics_info,
+	       verbose_info);
+	if (state.timer_handle >= 0)
+	   then CCTK_TimerStopI(state.timer_handle);
+	break;
 
-// output BH diagnostics?
-if (active_flag && IO_info.output_BH_diagnostics)
-   then {
-	if (AH_info.BH_diagnostics_fileptr == NULL)
-	   then AH_info.BH_diagnostics_fileptr
-		   = setup_BH_diagnostics_output_file(IO_info,
-						      hn, state.N_horizons);
-	output_BH_diagnostics_fn(BH_diagnostics,
-				 IO_info, hn,
-				 AH_info.BH_diagnostics_fileptr);
-	}
-			}
-		   else {
-			if (verbose_info.print_physics_details)
-			   then CCTK_VInfo(CCTK_THORNSTRING,
-					   "no apparent horizon found");
-			}
-		break;
-
-	default:
-		CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "\n"
-			   "   find_horizons(): unknown method=(int)%d!\n"
-			   "                    in main switch\n"
-			   "                    (this should never happen!)"
-			   ,
-			   int(state.method));			/*NOTREACHED*/
-		}
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   find_horizons(): unknown method=(int)%d!\n"
+"                    (this should never happen!)"
+		   ,
+		   int(state.method));				/*NOTREACHED*/
 	}
 
 if (state.timer_handle >= 0)
@@ -387,8 +296,10 @@ return data_ptr;
 //******************************************************************************
 
 //
-// This function implements  AHFinderDirect::method == "horizon function",
-// that is, it evaluates the Theta(h) function for a single apparent horizon.
+// This function implements  AHFinderDirect::method == "horizon function":
+// On processor #0 it evaluates the Theta(h) function for each apparent
+// horizon (and does any I/O desired); on other processors it does N_horizons
+// dummy evaluations on horizon #0.
 //
 // Note that if we decide to output h, we output it *after* any Theta(h)
 // evaluation or horizon finding has been done, to ensure that all the
@@ -398,63 +309,85 @@ return data_ptr;
 // timer_handle = a valid Cactus timer handle if we want to time the
 //		  apparent horizon process, or -ve to skip this
 //		  (we only time the computation, not the file I/O)
-// output_[hH] = flags to control whether or not we should output [hH]
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We evaluate Theta(h) as usual, but don't write
-//			   any output files.
-// hn = the horizon number (used only in naming output files)
-// N_horizons = the total number of horizon(s) being searched for number
-//		(used only in formatting info messages)
-//
-// Results:
-// This function returns true if the evaluation was successful, false
-// otherwise.
 //
 namespace {
-bool do_evaluate_expansion
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps,
-	 bool active_flag, int hn, int N_horizons)
+void do_evaluate_expansions(int my_proc, int N_horizons,
+			    horizon_sequence& hs,
+			       struct AH_info* const my_AH_info[],
+			    const struct cactus_grid_info& cgi,
+			    const struct geometry_info& gi,
+			    const struct IO_info& IO_info,
+			    bool output_h, bool output_Theta,
+			    const struct verbose_info& verbose_info,
+			    int timer_handle)
 {
-jtutil::norm<fp> Theta_norms;
-
-if (timer_handle >= 0)
-   then CCTK_TimerStartI(timer_handle);
-const bool status = expansion(ps, cgi, gi, false, true, &Theta_norms);
-if (timer_handle >= 0)
-   then CCTK_TimerStopI(timer_handle);
-if (!status)
-   then return false;					// *** ERROR RETURN ***
-
-if (Theta_norms.is_nonempty())	// might be empty if Theta(h) eval failed
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   Theta(h) rms-norm %.2e, infinity-norm %.2e",
-		   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
-
-if (active_flag && output_h)
-   then output_gridfn(ps, gfns::gfn__h,
-		      IO_info, IO_info.h_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-if (active_flag && output_Theta)
-   then output_gridfn(ps, gfns::gfn__Theta,
-		      IO_info, IO_info.Theta_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-
-return true;						// *** NORMAL RETURN ***
+const bool active_flag = (my_proc == 0);
+
+if (active_flag)
+   then {
+	assert( hs.N_horizons() == N_horizons );
+	assert( hs.my_N_horizons() == N_horizons );
+
+		for (int hn = hs.init_hn() ;
+		     hs.is_genuine() ;
+		     hn = hs.next_hn())
+		{
+		assert( my_AH_info[hn] != NULL );
+		struct AH_info& AH_info = *my_AH_info[hn];
+		patch_system& ps = *AH_info.ps_ptr;
+
+		if (timer_handle >= 0)
+		   then CCTK_TimerStartI(timer_handle);
+		jtutil::norm<fp> Theta_norms;
+		const bool Theta_ok = expansion(&ps,
+						cgi, gi,
+						false,	// no Jacobian coeffs
+						true,	// yes, print msgs
+						&Theta_norms);
+		if (timer_handle >= 0)
+		   then CCTK_TimerStopI(timer_handle);
+
+		if (output_h)
+		   then output_gridfn(ps, gfns::gfn__h,
+				      IO_info, IO_info.h_base_file_name,
+				      hn, verbose_info.print_algorithm_details);
+
+		if (Theta_ok)
+		   then {
+			CCTK_VInfo(CCTK_THORNSTRING,
+			   "   Theta(h) rms-norm %.2e, infinity-norm %.2e",
+			   Theta_norms.rms_norm(), Theta_norms.infinity_norm());
+			if (output_Theta)
+			   then output_gridfn(ps, gfns::gfn__Theta,
+					      IO_info, IO_info
+						       .Theta_base_file_name,
+					      hn, verbose_info
+						  .print_algorithm_details);
+			}
+		}
+	}
+   else {
+		for (int i = 0 ; i < N_horizons ; ++i)
+		{
+		expansion(NULL,	// dummy computation
+			  cgi, gi);
+		}
+	}
 }
 	  }
 
 //******************************************************************************
 
 //
-// This function implements  AHFinderDirect::method == "test Jacobian",
-// that is, it computes and prints the Jacobian matrix J[Theta(h)] for a
-// single apparent horizon.  The computation may optionally be done
-// in several different ways, in which case all the resulting Jacobian
-// matrices are printed, as are their differences.  Alternatively, only
+// This function implements
+//  AHFinderDirect::method == "test expansion Jacobians":
+// On processor #0 it computes and prints the Jacobian matrix J[Theta(h)]
+// function for horizon #1; on other processors it does dummy Jacobian
+// computations.
+//
+// The Jacobian computation may optionally be done in several different
+// ways, in which case all the resulting Jacobian matrices are printed,
+// as are their differences.  Alternatively, only
 // the numerical perturbation computation may be done/printed.
 //
 // Arguments:
@@ -468,259 +401,60 @@ return true;						// *** NORMAL RETURN ***
 //	  false ==> Just test/print the numerical perturbation calculation.
 //		    (This may be useful if one or more of the other methods
 //		    is broken.) 
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We compute the Jacobian(s) as usual, but don't
-//			   write any output files.
-// hn = the horizon number (used only in naming output files)
-// N_horizons = the total number of horizon(s) being searched for number
-//		(used only in formatting info messages)
-//
-// Results:
-// This function returns true if the Jacobian computation was successful,
-// false otherwise.
 //
 namespace {
-bool do_test_expansion_Jacobian
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info,
-	 bool test_all_Jacobian_compute_methods,
-	 struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons)
+void do_test_expansion_Jacobians(int my_proc, int N_horizons,
+				 struct AH_info* const my_AH_info[],
+				 const struct cactus_grid_info& cgi,
+				 const struct geometry_info& gi,
+				       struct Jacobian_info& Jac_info,
+				 bool test_all_Jacobian_compute_methods,
+				 const struct IO_info& IO_info,
+				 const struct verbose_info& verbose_info,
+				 int timer_handle)
 {
-// numerical perturbation
-Jacobian* Jac_NP_ptr = & Jac;
-if (!  expansion(ps, cgi, gi, true))
-   then return false;					// *** ERROR RETURN ***
+const bool active_flag = (my_proc == 0);
+assert(N_horizons >= 1);
+
+const bool print_msg_flag = true;
+const int hn = 1;
+
+struct AH_info* AH_info_ptr = active_flag ? my_AH_info[hn]      : NULL;
+patch_system*        ps_ptr = active_flag ? AH_info_ptr->ps_ptr : NULL;
+
+//
+// numerical-perturbation Jacobian
+//
+Jacobian* Jac_NP_ptr = active_flag ? AH_info_ptr->Jac_ptr : NULL;
+expansion(ps_ptr,
+	  cgi, gi);
 Jac_info.Jacobian_compute_method = Jacobian__numerical_perturbation;
-if (! expansion_Jacobian(ps, *Jac_NP_ptr,
-			 cgi, gi, Jac_info,
-			 true))	// print msgs
-   then return false;					// *** ERROR RETURN ***
+expansion_Jacobian(ps_ptr, Jac_NP_ptr,
+		   cgi, gi, Jac_info,
+		   print_msg_flag);
 
 Jacobian* Jac_SD_FDdr_ptr = NULL;
 if (test_all_Jacobian_compute_methods)
    then {
 	// symbolic differentiation with finite diff d/dr
-	Jac_SD_FDdr_ptr = new_Jacobian(Jac_info.Jacobian_store_solve_method,
-				       ps,
-				       verbose_info.print_algorithm_details);
-	if (! expansion(ps, cgi, gi, true))
-	   then return false;				// *** ERROR RETURN ***
+	Jac_SD_FDdr_ptr = active_flag
+			  ? new_Jacobian(Jac_info.Jacobian_store_solve_method,
+					 *ps_ptr,
+					 verbose_info.print_algorithm_details)
+			  : NULL;
+	expansion(ps_ptr,
+		  cgi, gi,
+		  true);		// compute SD Jacobian coeffs
 	Jac_info.Jacobian_compute_method = Jacobian__symbolic_diff_with_FD_dr;
-	if (! expansion_Jacobian(ps, *Jac_SD_FDdr_ptr,
-				 cgi, gi, Jac_info,
-				 true))	// print msgs
-	   then return false;				// *** ERROR RETURN ***
+	expansion_Jacobian(ps_ptr, Jac_SD_FDdr_ptr,
+			   cgi, gi, Jac_info,
+			   print_msg_flag);
 	}
 
 if (active_flag)
-   then output_Jacobians(ps,
+   then output_Jacobians(*ps_ptr,
 			 Jac_NP_ptr, Jac_SD_FDdr_ptr,
 			 IO_info, IO_info.Jacobian_base_file_name,
-			 hn, true);
-
-return true;						// *** NORMAL RETURN ***
-}
-	  }
-
-//******************************************************************************
-
-//
-// This function implements  AHFinderDirect::method == "find horizon",
-// that is, it finds the (an) apparent horizon.
-//
-// Arguments:
-// timer_handle = a valid Cactus timer handle if we want to time the
-//		  apparent horizon process, or -ve to skip this
-//		  (we only time the computation, not the file I/O)
-// initial_find_flag = true ==> This is the first time we've tried to find
-//				this horizon, or we've tried before but
-//				failed on our most recent previous attempt.
-//				Thus we should use "initial" parameters
-//				for the Newton iteration.
-//		       flase ==> We've tried to find this horizon before,
-//				 and we succeeded on our most recent
-//				 previous attempt.  Thus we should use
-//				 "subsequent" parameters for the Newton
-//				 iteration.
-// active_flag = Controls whether this processor is "active":
-//		 true ==> Normal operation.
-//		 false ==> We find the horizon as usual, but don't write
-//			   any output files.
-// hn = the horizon number (used only in naming output files)
-// N_horizons = the total number of horizon(s) being searched for number
-//		(used only in formatting info messages)
-//
-// Results:
-// This function returns true if the apparent horizon was found
-// successfully, false otherwise.
-//
-namespace {
-bool do_find_horizon
-	(const struct verbose_info& verbose_info, int timer_handle,
-	 struct IO_info& IO_info, bool output_h, bool output_Theta,
-	 const struct solver_info& solver_info, bool initial_find_flag,
-	 const struct Jacobian_info& Jac_info,
-	 struct cactus_grid_info& cgi, struct geometry_info& gi,
-	 patch_system& ps, Jacobian& Jac,
-	 bool active_flag, int hn, int N_horizons)
-{
-if (verbose_info.print_algorithm_highlights)
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "searching for horizon %d/%d",
-		   hn, N_horizons);
-
-if (timer_handle >= 0)
-   then CCTK_TimerStartI(timer_handle);
-const bool status = Newton_solve(ps, Jac,
-				 cgi, gi, Jac_info,
-				 solver_info, initial_find_flag, IO_info,
-				 active_flag, hn, verbose_info);
-if (timer_handle >= 0)
-   then CCTK_TimerStopI(timer_handle);
-if (! status)
-   then return false;					// *** ERROR RETURN ***
-
-if (active_flag && output_h)
-   then output_gridfn(ps, gfns::gfn__h,
-		      IO_info, IO_info.h_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-if (active_flag && output_Theta)
-   then output_gridfn(ps, gfns::gfn__Theta,
-		      IO_info, IO_info.Theta_base_file_name,
-		      hn, verbose_info.print_algorithm_details);
-
-return true;						// *** NORMAL RETURN ***
-}
-	  }
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// Given that an apparent horizon has been found, this function computes
-// various black hole diagnostics.
-//
-// Inputs (gridfns)
-// h		# ghosted
-// one		# nominal
-// global_[xyz]	# nominal
-//
-namespace {
-void compute_BH_diagnostics
-	(const patch_system& ps,
-	 const struct BH_diagnostics_info& BH_diagnostics_info,
-	 const struct verbose_info& verbose_info,
-	 struct BH_diagnostics& BH_diagnostics)
-{
-//
-// compute raw surface integrals
-//
-fp integral_one = surface_integral(ps, gfns::gfn__one,
-				   BH_diagnostics_info.integral_method);
-fp integral_h = surface_integral(ps, gfns::gfn__h,
-				 BH_diagnostics_info.integral_method);
-fp integral_x = surface_integral(ps, gfns::gfn__global_x,
-				 BH_diagnostics_info.integral_method);
-fp integral_y = surface_integral(ps, gfns::gfn__global_y,
-				 BH_diagnostics_info.integral_method);
-fp integral_z = surface_integral(ps, gfns::gfn__global_z,
-				 BH_diagnostics_info.integral_method);
-
-//
-// adjust integrals to take into account patch system not covering full sphere
-//
-switch	(ps.type())
-	{
-case patch_system::patch_system__full_sphere:
-	break;
-case patch_system::patch_system__plus_z_hemisphere:
-	integral_one *= 2.0;
-	integral_h *= 2.0;
-	integral_x *= 2.0;
-	integral_y *= 2.0;
-	integral_z = integral_one * ps.origin_z();
-	break;
-case patch_system::patch_system__plus_xy_quadrant_mirrored:
-case patch_system::patch_system__plus_xy_quadrant_rotating:
-	integral_one *= 4.0;
-	integral_h *= 4.0;
-	integral_x = integral_one * ps.origin_x();
-	integral_y = integral_one * ps.origin_y();
-	integral_z *= 4.0;
-	break;
-case patch_system::patch_system__plus_xz_quadrant_rotating:
-	integral_one *= 4.0;
-	integral_h *= 4.0;
-	integral_x = integral_one * ps.origin_x();
-	integral_y *= 4.0;
-	integral_z = integral_one * ps.origin_z();
-	break;
-case patch_system::patch_system__plus_xyz_octant_mirrored:
-case patch_system::patch_system__plus_xyz_octant_rotating:
-	integral_one *= 8.0;
-	integral_h *= 8.0;
-	integral_x = integral_one * ps.origin_x();
-	integral_y = integral_one * ps.origin_y();
-	integral_z = integral_one * ps.origin_z();
-	break;
-default:
-	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   compute_BH_diagnostics(): unknown ps.type()=(int)%d!\n"
-"                             (this should never happen!)"
-,
-		   int(ps.type()));				/*NOTREACHED*/
-	}
-
-BH_diagnostics.centroid_x = integral_x / integral_one;
-BH_diagnostics.centroid_y = integral_y / integral_one;
-BH_diagnostics.centroid_z = integral_z / integral_one;
-
-BH_diagnostics.area = integral_one;
-BH_diagnostics.circumference_xy = ps.xy_circumference
-    (gfns::gfn__h,
-     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-					     gfns::gfn__g_dd_33,
-     BH_diagnostics_info.integral_method);
-BH_diagnostics.circumference_xz = ps.xz_circumference
-    (gfns::gfn__h,
-     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-					     gfns::gfn__g_dd_33,
-     BH_diagnostics_info.integral_method);
-BH_diagnostics.circumference_yz = ps.yz_circumference
-    (gfns::gfn__h,
-     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-			 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-					     gfns::gfn__g_dd_33,
-     BH_diagnostics_info.integral_method);
-BH_diagnostics.mean_radius = integral_h / integral_one;
-BH_diagnostics.m_irreducible = sqrt(BH_diagnostics.area / (16.0*PI));
-}
-	  }
-
-//******************************************************************************
-
-//
-// This function computes the surface integral of a gridfn over the
-// horizon.
-//
-namespace {
-fp surface_integral(const patch_system& ps, int src_gfn,
-		    enum patch::integration_method method)
-{
-return ps.integrate_gridfn
-	   (src_gfn,
-	    gfns::gfn__h,
-	    gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
-				gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
-						    gfns::gfn__g_dd_33,
-	    method);
+			 hn, print_msg_flag);
 }
 	  }
diff --git a/src/driver/horizon_sequence.cc b/src/driver/horizon_sequence.cc
new file mode 100644
index 0000000..ca47489
--- /dev/null
+++ b/src/driver/horizon_sequence.cc
@@ -0,0 +1,71 @@
+// horizon_sequence.cc -- describes sequence of horizons a processor works on
+// $Header$
+
+//
+// horizon_sequence::horizon_sequence
+// horizon_sequence::~horizon_sequence
+// horizon_sequence::append_hn
+// horizon_sequence::next_posn
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include "horizon_sequence.hh"
+
+//******************************************************************************
+
+//
+// This function constructs an horizon sequence, which can hold a maximum
+// of N_horizons horizon numbers.
+//
+horizon_sequence::horizon_sequence(int N_horizons_in)
+	: N_horizons_(N_horizons_in),
+	  my_N_horizons_(0),		// sequence starts out empty
+	  posn_(-1),
+	  my_hn_(new int[N_horizons_in])
+{ }
+
+//******************************************************************************
+
+//
+// This function destroys a horizon sequence.
+//
+horizon_sequence::~horizon_sequence()
+{
+delete[] my_hn_;
+}
+
+//******************************************************************************
+
+//
+// This function appends  hn  to the sequence.  It returns the new value
+// of my_N_horizons().
+//
+int horizon_sequence::append_hn(int hn)
+{
+assert( hn > 0 );			// can only append genuine horizons
+assert( my_N_horizons_ < N_horizons_ );	// make sure there's space for it
+my_hn_[my_N_horizons_++] = hn;
+posn_ = 0;
+return my_N_horizons_;
+}
+
+//******************************************************************************
+
+//
+// This function computes the internal position immediately following
+// a given internal position in the sequence.
+//
+// Arguments:
+// p = (in) The current internal position, with posn_ semantics
+//
+// Results:
+// This function returns the next internal position after p.
+//
+int horizon_sequence::next_posn(int pos)
+	const
+{
+return   (pos < 0) ? pos-1
+       : (pos+1 < my_N_horizons_) ? pos+1
+       : -1;
+}
diff --git a/src/driver/horizon_sequence.hh b/src/driver/horizon_sequence.hh
new file mode 100644
index 0000000..398f687
--- /dev/null
+++ b/src/driver/horizon_sequence.hh
@@ -0,0 +1,117 @@
+// horizon_sequence.hh -- describes sequence of horizons a processor works on
+// $Header$
+
+//
+// A  horizon_sequence  object describes the sequence of horizons a
+// (the current) procesor works on.  This is some sequence of genuine
+// horizons, followed by a dummy horizon, the latter repeating indefinitely.
+//
+// A horizon is specified by its (globally unique) "horizon number",
+// which is 0 for a dummy horizon, or [1,N_horizons] for a genuine horizon.
+//
+// Thus a typical sequence of horizon numbers might be
+//	1, 3, 4, 0, 0, 0, 0, ...
+//
+class	horizon_sequence
+	{
+public:
+	//
+	// ***** query functions *****
+	//
+
+	// how many (genuine) horizons are there in total?
+	int N_horizons() const { return N_horizons_; }
+
+	// how many genuine horizons are in the sequence (for this processor)?
+	int my_N_horizons() const { return my_N_horizons_; }
+
+	// are there any genuine horizons in the sequence (for this processor)?
+	bool has_genuine_horizons() const { return my_N_horizons_ > 0; }
+
+	// is the current horizon in the sequence dummy/genuine?
+	bool is_dummy  () const { return posn_is_dummy  (posn_); }
+	bool is_genuine() const { return posn_is_genuine(posn_); }
+
+	// what will  is_genuine()  return after  next_hn()  is called?
+	// i.e. is this *not* the final genuine  hn  in the sequence?
+	bool is_next_genuine() const
+		{ return posn_is_genuine( next_posn(posn_) ); }
+
+	// return 0 if current hn is genuine,
+	// or 1-origin ordinal number of dummy if it's dummy,
+	//    i.e. 1 for 1st dummy, 2 for 2nd dummy, 3 for 3rd dummy, ...
+	int dummy_number() const { return is_genuine() ? 0 : -posn_; }
+
+	// get current hn in sequence
+	int get_hn() const
+		{ return posn_is_genuine(posn_) ? my_hn_[posn_] : 0; }
+
+	//
+	// ***** traverse the sequence *****
+	//
+	// the idiom to traverse the genuine horizons is
+	//    for (int hn = hs.init_hn() ; hs.is_genuine() ; hn = hs.next_hn())
+	//    {
+	//    }
+	//
+	// the idiom to traverse the genuine horizons, followed by
+	// infinite repetition of the dummy horizon, is
+	//    for (int hn = hs.init_hn() ; ; hn = hs.next_hn())
+	//    {
+	//    }
+	//
+
+	// reset sequence, return starting genuine hn
+	int init_hn()
+		{
+		posn_ = (my_N_horizons_ == 0) ? -1 : 0;
+		return get_hn();
+		}
+
+	// get next hn in sequence
+	int next_hn() { posn_ = next_posn(posn_); return get_hn(); }
+
+
+	//
+	// ***** set up the sequence *****
+	//
+
+	// construct an empty horizon sequence
+	// which can hold <= N_horizons horizon numbers
+	horizon_sequence(int N_horizons);
+	~horizon_sequence();
+
+	// append hn to the sequence
+	// ... returns new my_N_horizons()
+	int append_hn(int hn);
+
+private:
+	// is a specified posn dummy/genuine
+	bool posn_is_genuine(int pos) const
+		{ return (pos >= 0) && (pos < my_N_horizons_); }
+	bool posn_is_dummy(int pos) const
+		{ return !posn_is_genuine(pos); }
+
+	// what is the next posn in the sequence
+	int next_posn(int pos) const;
+
+private:
+	// we forbid copying and passing by value
+	// by declaring the copy constructor and assignment operator
+	// private, but never defining them
+	horizon_sequence(const horizon_sequence& rhs);
+	horizon_sequence& operator=(const horizon_sequence& rhs);
+
+private:
+	const int N_horizons_;
+	int my_N_horizons_;
+
+	// "internal position" in sequence
+	// this is a [0,my_N_horizons_) index into my_hn_[]
+	//           for genuine horizons,
+	//         or < 0 for the dummy horizon, with the absolute value
+	//            counting the number of times we've returned hn=0
+	int posn_;
+
+	int* my_hn_;	// --> new[]-allocated array of genuine horizon numbers
+	};
diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc
index a58617b..bc9ab89 100644
--- a/src/driver/initial_guess.cc
+++ b/src/driver/initial_guess.cc
@@ -41,6 +41,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
 //******************************************************************************
diff --git a/src/driver/io.cc b/src/driver/io.cc
index defbe0f..2bddaa9 100644
--- a/src/driver/io.cc
+++ b/src/driver/io.cc
@@ -1,14 +1,9 @@
 // io.cc -- I/O routines for this thorn
 // $Header$
 //
-// <<<prototypes for functions local to this file>>>
 // input_gridfn - read an angular grid function from an input file
 // output_gridfn - write an angular grid function to an output file
 // output_Jacobians - write a Jacobian matrix or matrices to an output file
-// decode_horizon_file_format - decode the horizon_file_format parameter
-//
-// setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file
-// output_BH_diagnostics_fn - append BH diagnostics to file
 //
 /// io_file_name
 //
@@ -45,12 +40,13 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
 //******************************************************************************
 
 //
-// ***** prototypes for functions local to this file ***** 
+// prototypes for functions local to this file
 //
 
 namespace {
@@ -301,110 +297,6 @@ fclose(fileptr);
 }
 
 //******************************************************************************
-
-//
-// This function decodes the  horizon_file_format  parameter (string) into an
-// internal enum for future use.
-//
-enum horizon_file_format
-  decode_horizon_file_format(const char horizon_file_format_string[])
-{
-if	(STRING_EQUAL(horizon_file_format_string, "ASCII (gnuplot)"))
-   then return horizon_file_format__ASCII_gnuplot;
-else if (STRING_EQUAL(horizon_file_format_string, "HDF5"))
-   then return horizon_file_format__HDF5;
-else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-"decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!",
-		   horizon_file_format_string);			/*NOTREACHED*/
-}
-
-//******************************************************************************
-//******************************************************************************
-//******************************************************************************
-
-//
-// This function creates a BH-diagnostics output file, writes a suitable
-// header comment, and returns a stdio file pointer which can be used to
-// append data to the file.
-//
-FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
-				       int hn, int N_horizons)
-{
-const int N_file_name_buffer = 200;
-char file_name_buffer[N_file_name_buffer];
-
-snprintf(file_name_buffer, N_file_name_buffer,
-	 "%s.ah%d.%s",
-	 IO_info.BH_diagnostics_base_file_name,
-	 hn, IO_info.BH_diagnostics_file_name_extension);
-FILE *fileptr = fopen(file_name_buffer, "w");
-if (fileptr == NULL)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   setup_BH_diagnostics_output_file():\n"
-"        can't open BH-diagnostics output file\n"
-"        \"%s\"!"
-		   ,
-		   file_name_buffer);				/*NOTREACHED*/
-
-fprintf(fileptr, "# apparent horizon %d/%d\n", hn, N_horizons);
-fprintf(fileptr, "#\n");
-fprintf(fileptr, "# column 1 = cctk_iteration\n");
-fprintf(fileptr, "# column 2 = cctk_time\n");
-fprintf(fileptr, "# column 3 = centroid_x\n");
-fprintf(fileptr, "# column 4 = centroid_y\n");
-fprintf(fileptr, "# column 5 = centroid_z\n");
-fprintf(fileptr, "# column 6 = mean radius\n");
-fprintf(fileptr, "# column 7 = xy-plane circumference\n");
-fprintf(fileptr, "# column 8 = xz-plane circumference\n");
-fprintf(fileptr, "# column 9 = yz-plane circumference\n");
-fprintf(fileptr, "# column 10 = area\n");
-fprintf(fileptr, "# column 11 = m_irreducible\n");
-fflush(fileptr);
-
-return fileptr;
-}
-
-//******************************************************************************
-
-//
-// This function appends a BH-diagnostics line to a data file.  It
-// attempts to ensure that the newly-written line is flushed to disk,
-// so the output file can be examined while the Cactus run is still in
-// progress.
-//
-// (The "_fn" in the function name is to distinguish it from the Cactus
-// parameter  output_BH_diagnostics .)
-//
-// Arguments:
-// BH_diagnostics = The BH diagnostics to be written
-// fileptr = The stdio file pointer to append to
-//
-void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
-			      const struct IO_info& IO_info,
-			      int hn, FILE* fileptr)
-{
-assert(fileptr != NULL);
-
-fprintf(fileptr,
-//  cctk_iteration        mean radius                     area    m_irreducible
-//  ==  cctk_time         ======  {xy,xz,yz}-plane        ======  ======
-//  ==  ====  centroid_[xyz]      circumferences          ======  ======
-//  ==  ====  ==========  ======  ======================  ======  ======
-   "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
-	IO_info.time_iteration, double(IO_info.time),
-	BH_diagnostics.centroid_x, BH_diagnostics.centroid_y,
-				   BH_diagnostics.centroid_z,
-	BH_diagnostics.mean_radius,
-	BH_diagnostics.circumference_xy, BH_diagnostics.circumference_xz,
-					 BH_diagnostics.circumference_yz,
-	BH_diagnostics.area,
-	BH_diagnostics.m_irreducible);
-
-fflush(fileptr);
-}
-
-//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
diff --git a/src/driver/make.code.defn b/src/driver/make.code.defn
index bb7cb46..879219e 100644
--- a/src/driver/make.code.defn
+++ b/src/driver/make.code.defn
@@ -2,9 +2,10 @@
 # $Header$
 
 # Source files in this directory
-SRCS = setup.cc find_horizons.cc \
-       state.cc \
-       initial_guess.cc Newton.cc io.cc
+SRCS = setup.cc find_horizons.cc state.cc \
+       initial_guess.cc Newton.cc \
+       BH_diagnostics.cc io.cc \
+       horizon_sequence.cc
 
 # Subdirectories containing source files
 SUBDIRS =
diff --git a/src/driver/makefile b/src/driver/makefile
index e95ddc3..bbc3443 100644
--- a/src/driver/makefile
+++ b/src/driver/makefile
@@ -1,22 +1,28 @@
-# Makefile for GR code
-# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/src/driver/makefile,v 1.2 2002-09-16 14:17:16 jthorn Exp $
+# Makefile for driver code
+# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/src/driver/makefile,v 1.3 2003-02-15 18:43:08 jthorn Exp $
 
 #
-# Environment Variables:
-#	MAPLE_VERSION	used via @ifdef for version control in Maple code;
-#			typically set to something like  MAPLE_V_RELEASE_4
-#			(not presently used, but may be needed in the future)
-#
 # Targets:
 #	ellipsoid	==> compute ellipsoid equations
+#	test		==> build standalone test drivers
+#	clean		==> remove standalone test drivers
 #
 
-ifneq ($(MAPLE_VERSION),)
-MPP_FLAGS := -D$(MAPLE_VERSION)
-endif
+CXXFLAGS := -g $(STD_GXX_FLAGS)
 
 ###############################################################################
 
 .PHONY		: ellipsoid
 ellipsoid	:
 	maple <ellipsoid.maple >ellipsoid.log
+
+.PHONY		: test
+test		: test_horizon_sequence
+
+.PHONY			: test_horizon_sequence
+test_horizon_sequence	:
+	$(CXX) $(CXXFLAGS) -o $@ test_horizon_sequence.cc horizon_sequence.cc
+
+.PHONY		: clean
+clean		:
+	-rm *core test_horizon_sequence
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 444e2b1..df6e76e 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -4,11 +4,13 @@
 // <<<prototypes for functions local to this file>>>
 // <<<access to persistent data>>>
 //
-// AHFinderDirect_driver - top-level driver to setup persistent data structures
-// Cactus_gridfn_data_ptr - get Cactus data pointer for CCTK_REAL gridfn
+// AHFinderDirect_setup - top-level driver to setup persistent data structures
 ///
-/// decode_method - decode the method parameter
-/// decode_verbose_level - decode the verbose_level parameter
+/// decode_method - decode the  method  parameter
+/// decode_verbose_level - decode the  verbose_level  parameter
+/// decode_horizon_file_format - decode the  horizon_file_format  parameter
+///
+/// allocate_horizons_to_processor - choose which horizons this proc will find
 ///
 /// choose_patch_system_type - choose patch system type
 ///
@@ -45,6 +47,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
 //******************************************************************************
@@ -58,6 +61,13 @@ enum method
   decode_method(const char method_string[]);
 enum verbose_level
   decode_verbose_level(const char verbose_level_string[]);
+enum horizon_file_format
+  decode_horizon_file_format(const char horizon_file_format_string[]);
+
+int allocate_horizons_to_processor(int N_procs, int my_proc,
+				   int N_horizons, bool multiproc_flag,
+				   horizon_sequence& my_hs,
+				   const struct verbose_info& verbose_info);
 
 enum patch_system::patch_system_type
   choose_patch_system_type(const char grid_domain[],
@@ -89,16 +99,10 @@ DECLARE_CCTK_PARAMETERS
 CCTK_VInfo(CCTK_THORNSTRING,
 	   "setting up AHFinderDirect data structures");
 
-state.N_horizons = N_horizons;
-CCTK_VInfo(CCTK_THORNSTRING,
-           "           to search for %d horizon%s",
-	   state.N_horizons, ((state.N_horizons == 1) ? "" : "s"));
-
 
 //
-// decode/copy parameters into structures
+// basic setup
 //
-
 state.method = decode_method(method);
 
 struct verbose_info& verbose_info = state.verbose_info;
@@ -113,50 +117,20 @@ verbose_info.print_algorithm_details
    = (state.verbose_info.verbose_level >= verbose_level__algorithm_details);
 
 state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1;
-state.my_proc = CCTK_MyProc(cctkGH);
-state.N_procs = CCTK_nProcs(cctkGH);
-
-struct IO_info& IO_info = state.IO_info;
-IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
-IO_info.output_initial_guess = (output_initial_guess != 0);
-IO_info.how_often_to_output_h     = how_often_to_output_h;
-IO_info.how_often_to_output_Theta = how_often_to_output_Theta;
-IO_info.output_ghost_zones_for_h  = (output_ghost_zones_for_h != 0);
-IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_file_name_extension;
-IO_info.HDF5_file_name_extension          = HDF5_file_name_extension;
-IO_info.h_base_file_name         = h_base_file_name;
-IO_info.Theta_base_file_name     = Theta_base_file_name;
-IO_info.Delta_h_base_file_name   = Delta_h_base_file_name;
-IO_info.Jacobian_base_file_name  = Jacobian_base_file_name;
-IO_info.output_BH_diagnostics              = (output_BH_diagnostics != 0);
-IO_info.BH_diagnostics_base_file_name      = BH_diagnostics_base_file_name;
-IO_info.BH_diagnostics_file_name_extension = BH_diagnostics_file_name_extension;
-IO_info.time_iteration = 0;
-IO_info.time           = 0.0;
 
-struct Jacobian_info& Jac_info = state.Jac_info;
-Jac_info.Jacobian_compute_method
-	= decode_Jacobian_compute_method(Jacobian_compute_method);
-Jac_info.Jacobian_store_solve_method
-	= decode_Jacobian_store_solve_method(Jacobian_store_solve_method);
-Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude;
+state.N_procs = CCTK_nProcs(cctkGH);
+state.my_proc = CCTK_MyProc(cctkGH);
 
-struct solver_info& solver_info = state.solver_info;
-solver_info.debugging_output_at_each_Newton_iteration
-			= (debugging_output_at_each_Newton_iteration != 0);
-solver_info.max_Newton_iterations__initial
-			= max_Newton_iterations__initial;
-solver_info.max_Newton_iterations__subsequent
-			= max_Newton_iterations__subsequent;
-solver_info.max_Delta_h_over_h           = max_Delta_h_over_h;
-solver_info.Theta_norm_for_convergence       = Theta_norm_for_convergence;
-solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
-solver_info.final_Theta_update_if_Delta_h_converged
-	= (final_Theta_update_if_Delta_h_converged != 0);
+state.N_horizons = N_horizons;
+state.N_active_procs = 0;	// dummy value, will be set properly later
+CCTK_VInfo(CCTK_THORNSTRING,
+           "           to search for %d horizon%s on %d processor%s",
+	   state.N_horizons, ((state.N_horizons == 1) ? "" : "s"),
+	   state.N_procs, ((state.N_procs == 1) ? "" : "s"));
 
 
 //
-// set up the Cactus grid info
+// Cactus grid info
 //
 if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
@@ -211,7 +185,7 @@ cgi.psi_data     = NULL;		// dummy value, will be reset later
 
 
 //
-// set up the geometry parameters and the geometry interpolator
+// geometry info
 //
 if (verbose_info.print_algorithm_highlights)
    then CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
@@ -243,6 +217,53 @@ gi.check_that_geometry_is_finite = (check_that_geometry_is_finite != 0);
 
 
 //
+// Jacobian info
+//
+struct Jacobian_info& Jac_info = state.Jac_info;
+Jac_info.Jacobian_compute_method
+	= decode_Jacobian_compute_method(Jacobian_compute_method);
+Jac_info.Jacobian_store_solve_method
+	= decode_Jacobian_store_solve_method(Jacobian_store_solve_method);
+Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude;
+
+
+//
+// elliptic-solver info
+//
+struct solver_info& solver_info = state.solver_info;
+solver_info.debugging_output_at_each_Newton_iteration
+			= (debugging_output_at_each_Newton_iteration != 0);
+solver_info.max_Newton_iterations__initial
+			= max_Newton_iterations__initial;
+solver_info.max_Newton_iterations__subsequent
+			= max_Newton_iterations__subsequent;
+solver_info.max_allowable_Delta_h_over_h = max_allowable_Delta_h_over_h;
+solver_info.Theta_norm_for_convergence   = Theta_norm_for_convergence;
+
+
+//
+// I/O info
+//
+struct IO_info& IO_info = state.IO_info;
+IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
+IO_info.output_initial_guess = (output_initial_guess != 0);
+IO_info.how_often_to_output_h     = how_often_to_output_h;
+IO_info.how_often_to_output_Theta = how_often_to_output_Theta;
+IO_info.output_ghost_zones_for_h  = (output_ghost_zones_for_h != 0);
+IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_file_name_extension;
+IO_info.HDF5_file_name_extension          = HDF5_file_name_extension;
+IO_info.h_base_file_name         = h_base_file_name;
+IO_info.Theta_base_file_name     = Theta_base_file_name;
+IO_info.Delta_h_base_file_name   = Delta_h_base_file_name;
+IO_info.Jacobian_base_file_name  = Jacobian_base_file_name;
+IO_info.output_BH_diagnostics              = (output_BH_diagnostics != 0);
+IO_info.BH_diagnostics_base_file_name      = BH_diagnostics_base_file_name;
+IO_info.BH_diagnostics_file_name_extension = BH_diagnostics_file_name_extension;
+IO_info.time_iteration = 0;
+IO_info.time           = 0.0;
+
+
+//
 // other misc setup
 //
 state.BH_diagnostics_info.integral_method
@@ -250,7 +271,36 @@ state.BH_diagnostics_info.integral_method
 
 
 //
-// create AH-specific info for each AH
+// (genuine) horizon sequence for this processor
+//
+state.my_hs = new horizon_sequence(state.N_horizons);
+horizon_sequence& hs = *state.my_hs;
+
+//
+// if we're going to actually find horizons
+//    we spread the horizons over multiple processors for maximum efficiency,
+//    otherwise (we're just doing testing/debugging computations, so) we
+//    allocate all the horizons to processor #0 for simplicity
+//
+const bool multiproc_flag = (state.method == method__find_horizons);
+state.N_active_procs
+   = allocate_horizons_to_processor(state.N_procs, state.my_proc,
+				    state.N_horizons, multiproc_flag,
+				    hs,
+				    verbose_info);
+
+// horizon numbers run from 1 to N_horizons inclusive
+state.my_AH_info = (N_horizons == 0) ? NULL : new AH_info*[N_horizons+1];
+	  {
+	for (int hn = 1 ; hn <= N_horizons ; ++hn)
+	{
+	state.my_AH_info[hn] = NULL;
+	}
+	  }
+
+
+//
+// horizon-specific info for each horizon
 //
 
 // set up the interpatch interpolator
@@ -262,23 +312,23 @@ if (interp_handle < 0)
 "AHFinderDirect_setup(): couldn't find interpolator \"%s\"!",
 		   interpatch_interpolator_name);		/*NOTREACHED*/
 const int interp_param_table_handle
-	= Util_TableCreateFromString(interpatch_interpolator_pars);
+   = Util_TableCreateFromString(interpatch_interpolator_pars);
 if (interp_param_table_handle < 0)
    then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "AHFinderDirect_setup(): bad interpatch-interpolator parameter(s) \"%s\"!",
 		   interpatch_interpolator_pars);		/*NOTREACHED*/
 
-// horizon numbers run from 1 to N_horizons; we don't use horizon number hn=0
-state.AH_info_array = new AH_info[N_horizons+1];
-
-	for (int hn = 1 ; hn <= N_horizons ; ++hn)
+// setup all the genuine horizons on this processor
+	  {
+	for (int hn = hs.init_hn() ; hs.is_genuine() ; hn = hs.next_hn())
 	{
-	struct AH_info& AH_info = state.AH_info_array[hn];
+	state.my_AH_info[hn] = new AH_info;
+	struct AH_info& AH_info = *state.my_AH_info[hn];
 
 	if (verbose_info.print_algorithm_highlights)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "   setting up data structures for horizon %d/%d",
-			   hn, int(N_horizons));
+			   "   setting up data structures for horizon %d",
+			   hn);
 
 	// decide what type of patch system this one should be
 	const enum patch_system::patch_system_type ps_type
@@ -298,27 +348,28 @@ state.AH_info_array = new AH_info[N_horizons+1];
 	     : patch_system::type_of_name(patch_system_type[hn]);
 
 	// create the patch system
-	patch_system& ps
-	   = *new
-	      patch_system(origin_x[hn], origin_y[hn], origin_z[hn],
-			   ps_type,
-			   ghost_zone_width, patch_overlap_width,
-			   N_zones_per_right_angle[hn],
-			   gfns::nominal_min_gfn, gfns::nominal_max_gfn,
-			   gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
-			   interp_handle, interp_param_table_handle,
-			   true, verbose_info.print_algorithm_details);
-	AH_info.ps_ptr = &ps;
+	AH_info.ps_ptr
+	   = new patch_system(origin_x[hn], origin_y[hn], origin_z[hn],
+			      ps_type,
+			      ghost_zone_width, patch_overlap_width,
+			      N_zones_per_right_angle[hn],
+			      gfns::nominal_min_gfn, gfns::nominal_max_gfn,
+			      gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
+			      interp_handle, interp_param_table_handle,
+			      true, verbose_info.print_algorithm_details);
+	patch_system& ps = *AH_info.ps_ptr;
 	ps.set_gridfn_to_constant(1.0, gfns::gfn__one);
 
 	// Jacobian matrix will be created later
-	AH_info.Jac_ptr = NULL;
+	AH_info.Jac_ptr = new_Jacobian(Jac_info.Jacobian_store_solve_method,
+				       ps,
+				       verbose_info.print_algorithm_details);
+
+	AH_info.initial_find_flag = true;
 
 	if (verbose_info.print_algorithm_details)
 	   then CCTK_VInfo(CCTK_THORNSTRING,
 			   "      setting initial guess parameters etc");
-
-	// store the initial guess parameters
 	AH_info.initial_guess_info.method
 		= decode_initial_guess_method(initial_guess_method[hn]);
 	// ... Kerr/Kerr
@@ -366,17 +417,10 @@ state.AH_info_array = new AH_info[N_horizons+1];
 	AH_info.initial_guess_info.coord_ellipsoid_info.z_radius
 		= initial_guess__coord_ellipsoid__z_radius[hn];
 
-	// initialize the BH diagnostics
-	AH_info.AH_found = false;
-	AH_info.BH_diagnostics.centroid_x = 0.0;
-	AH_info.BH_diagnostics.centroid_y = 0.0;
-	AH_info.BH_diagnostics.centroid_z = 0.0;
-	AH_info.BH_diagnostics.area = 0.0;
-	AH_info.BH_diagnostics.m_irreducible = 0.0;
-
-	// other misc stuff
+	AH_info.found_flag = false;
 	AH_info.BH_diagnostics_fileptr = NULL;
 	}
+	  }
 }
 
 //******************************************************************************
@@ -391,12 +435,12 @@ namespace {
 enum method
   decode_method(const char method_string[])
 {
-if	(STRING_EQUAL(method_string, "evaluate expansion"))
-   then return method__evaluate_expansion;
-else if (STRING_EQUAL(method_string, "test expansion Jacobian"))
-   then return method__test_expansion_Jacobian;
-else if (STRING_EQUAL(method_string, "find horizon"))
-   then return method__find_horizon;
+if	(STRING_EQUAL(method_string, "evaluate expansions"))
+   then return method__evaluate_expansions;
+else if (STRING_EQUAL(method_string, "test expansion Jacobians"))
+   then return method__test_expansion_Jacobians;
+else if (STRING_EQUAL(method_string, "find horizons"))
+   then return method__find_horizons;
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "decode_method(): unknown method_string=\"%s\"!",
 		   method_string);				/*NOTREACHED*/
@@ -428,6 +472,81 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	  }
 
 //******************************************************************************
+
+//
+// This function decodes the  horizon_file_format  parameter (string) into an
+// internal enum for future use.
+//
+namespace {
+enum horizon_file_format
+  decode_horizon_file_format(const char horizon_file_format_string[])
+{
+if	(STRING_EQUAL(horizon_file_format_string, "ASCII (gnuplot)"))
+   then return horizon_file_format__ASCII_gnuplot;
+else if (STRING_EQUAL(horizon_file_format_string, "HDF5"))
+   then return horizon_file_format__HDF5;
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!",
+		   horizon_file_format_string);			/*NOTREACHED*/
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function allocates horizons to a processor.  That is, it chooses
+// which horizons a given processor (in practice, the current processor)
+// will attempt to find.  It records this choice in a  horizon_sequence
+// object.
+//
+// Arguments:
+// N_procs = The total number of processors.
+// my_proc = My processor number (0-origin).
+// N_horizons = The total number of (genuine) horizons to be found.
+// multiproc_flag = true to spread the horizons over multiple processors,
+//		    false to allocate all horizons to processor #0.
+// my_hs = (out) This will be set to the sequence of (genuine) horizons
+//		 allocated to the specified processor.
+// verbose_info = Specifies what information to print about the allocation.
+//
+// Results:
+// This function returns the total number of active processors, i.e. the
+// total number of processors which are allocated one or more genuine
+// horizons allocated.
+//
+namespace {
+int allocate_horizons_to_processor(int N_procs, int my_proc,
+				   int N_horizons, bool multiproc_flag,
+				   horizon_sequence& my_hs,
+				   const struct verbose_info& verbose_info)
+{
+const int N_active_procs = multiproc_flag ? jtutil::min(N_procs, N_horizons)
+					  : 1;
+
+//
+// Implementation note:
+// We allocate the horizons to active processors in round-robin order.
+//
+int proc = 0;
+	for (int hn = 1 ; hn <= N_horizons ; ++hn)
+	{
+	if (verbose_info.print_algorithm_highlights)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "   allocating horizon %d to processor #%d",
+			   hn, proc);
+	if (proc == my_proc)
+	   then my_hs.append_hn(hn);
+	if (++proc >= N_active_procs)
+	   then proc = 0;
+	}
+
+return N_active_procs;
+}
+	  }
+
+//******************************************************************************
 //******************************************************************************
 //******************************************************************************
 
diff --git a/src/driver/state.cc b/src/driver/state.cc
index 77dc1f2..a6092c7 100644
--- a/src/driver/state.cc
+++ b/src/driver/state.cc
@@ -32,6 +32,7 @@ using jtutil::error_exit;
 #include "../gr/gfns.hh"
 #include "../gr/gr.hh"
 
+#include "horizon_sequence.hh"
 #include "driver.hh"
 
 //******************************************************************************
diff --git a/src/driver/test_horizon_sequence.cc b/src/driver/test_horizon_sequence.cc
new file mode 100644
index 0000000..04e0978
--- /dev/null
+++ b/src/driver/test_horizon_sequence.cc
@@ -0,0 +1,223 @@
+// test_horizon_sequence.cc -- test driver for horizon_sequence class
+// $Header$
+
+#include <stdio.h>
+#include <assert.h>
+#include "horizon_sequence.hh"
+
+// prototypes for functions local to this file
+namespace {
+void test0();
+void test1();
+void test2();
+void test345();
+	  }
+
+int main(void)
+{
+test0();
+test1();
+test2();
+test345();
+
+printf("all ok!\n");
+return 0;
+}
+
+//******************************************************************************
+
+namespace {
+void test0()
+{
+printf("capacity 0 ...\n");
+horizon_sequence hs0(0);
+assert( hs0.N_horizons() == 0 );
+assert( hs0.my_N_horizons() == 0 );
+assert( ! hs0.has_genuine_horizons() );
+assert( hs0.is_dummy() );
+assert( ! hs0.is_genuine() );
+assert( ! hs0.is_next_genuine() );
+assert( hs0.dummy_number() == 1 );
+assert( hs0.init_hn() == 0 );
+	assert( hs0.get_hn() == 0 );
+	assert( hs0.dummy_number() == 1 );
+assert( hs0.next_hn() == 0 );
+	assert( hs0.get_hn() == 0 );
+	assert( hs0.dummy_number() == 2 );
+assert( hs0.next_hn() == 0 );
+	assert( hs0.get_hn() == 0 );
+	assert( hs0.dummy_number() == 3 );
+assert( hs0.next_hn() == 0 );
+	assert( hs0.get_hn() == 0 );
+	assert( hs0.dummy_number() == 4 );
+assert( hs0.N_horizons() == 0 );
+assert( hs0.my_N_horizons() == 0 );
+}
+	  }
+
+//******************************************************************************
+
+namespace {
+void test1()
+{
+printf("capacity 1 ...\n");
+horizon_sequence hs1(1);
+assert( hs1.N_horizons() == 1 );
+assert( hs1.my_N_horizons() == 0 );
+assert( ! hs1.has_genuine_horizons() );
+assert( hs1.append_hn(42) == 1 );
+assert( hs1.N_horizons() == 1 );
+assert( hs1.my_N_horizons() == 1 );
+assert( hs1.has_genuine_horizons() );
+assert( hs1.is_genuine() );
+assert( ! hs1.is_next_genuine() );
+assert( hs1.dummy_number() == 0 );
+assert( hs1.init_hn() == 42 );
+	assert( hs1.get_hn() == 42 );
+	assert( hs1.is_genuine() );
+	assert( ! hs1.is_next_genuine() );
+	assert( hs1.dummy_number() == 0 );
+assert( hs1.next_hn() == 0 );
+	assert( hs1.get_hn() == 0 );
+	assert( ! hs1.is_genuine() );
+	assert( ! hs1.is_next_genuine() );
+	assert( hs1.dummy_number() == 1 );
+assert( hs1.next_hn() == 0 );
+	assert( hs1.get_hn() == 0 );
+	assert( ! hs1.is_genuine() );
+	assert( ! hs1.is_next_genuine() );
+	assert( hs1.dummy_number() == 2 );
+assert( hs1.next_hn() == 0 );
+	assert( hs1.get_hn() == 0 );
+	assert( ! hs1.is_genuine() );
+	assert( ! hs1.is_next_genuine() );
+	assert( hs1.dummy_number() == 3 );
+assert( hs1.next_hn() == 0 );
+	assert( hs1.get_hn() == 0 );
+	assert( ! hs1.is_genuine() );
+	assert( ! hs1.is_next_genuine() );
+	assert( hs1.dummy_number() == 4 );
+assert( hs1.N_horizons() == 1 );
+assert( hs1.my_N_horizons() == 1 );
+}
+	  }
+
+//******************************************************************************
+
+namespace {
+void test2()
+{
+printf("capacity 2 ...\n");
+horizon_sequence hs2(2);
+assert( hs2.N_horizons() == 2 );
+assert( hs2.my_N_horizons() == 0 );
+assert( ! hs2.has_genuine_horizons() );
+assert( hs2.append_hn(42) == 1 );
+assert( hs2.append_hn(69) == 2 );
+assert( hs2.N_horizons() == 2 );
+assert( hs2.my_N_horizons() == 2 );
+assert( hs2.has_genuine_horizons() );
+assert( hs2.is_genuine() );
+assert( hs2.is_next_genuine() );
+assert( hs2.dummy_number() == 0 );
+assert( hs2.init_hn() == 42 );
+	assert( hs2.get_hn() == 42 );
+	assert( hs2.is_genuine() );
+	assert( hs2.is_next_genuine() );
+	assert( hs2.dummy_number() == 0 );
+assert( hs2.next_hn() == 69 );
+	assert( hs2.get_hn() == 69 );
+	assert( hs2.is_genuine() );
+	assert( ! hs2.is_next_genuine() );
+	assert( hs2.dummy_number() == 0 );
+assert( hs2.next_hn() == 0 );
+	assert( hs2.get_hn() == 0 );
+	assert( ! hs2.is_genuine() );
+	assert( ! hs2.is_next_genuine() );
+	assert( hs2.dummy_number() == 1 );
+assert( hs2.next_hn() == 0 );
+	assert( hs2.get_hn() == 0 );
+	assert( ! hs2.is_genuine() );
+	assert( ! hs2.is_next_genuine() );
+	assert( hs2.dummy_number() == 2 );
+assert( hs2.next_hn() == 0 );
+	assert( hs2.get_hn() == 0 );
+	assert( ! hs2.is_genuine() );
+	assert( ! hs2.is_next_genuine() );
+	assert( hs2.dummy_number() == 3 );
+assert( hs2.N_horizons() == 2 );
+assert( hs2.my_N_horizons() == 2 );
+}
+	  }
+
+//******************************************************************************
+
+namespace {
+void test345()
+{
+	for (int capacity = 3 ; capacity <= 5 ; ++capacity)
+	{
+	printf("capacity %d ...\n", capacity);
+	horizon_sequence hs(capacity);
+	assert( hs.N_horizons() == capacity );
+	assert( hs.my_N_horizons() == 0 );
+	assert( ! hs.has_genuine_horizons() );
+
+	assert( hs.init_hn() == 0 );
+	assert( hs.dummy_number() == 1 );
+
+	assert( hs.next_hn() == 0 );
+	assert( hs.dummy_number() == 2 );
+
+	assert( hs.next_hn() == 0 );
+	assert( hs.dummy_number() == 3 );
+
+	assert( hs.append_hn(42) == 1 );
+	assert( hs.append_hn(69) == 2 );
+	assert( hs.append_hn(105) == 3 );
+	assert( hs.N_horizons() == capacity );
+	assert( hs.my_N_horizons() == 3 );
+	assert( hs.has_genuine_horizons() );
+	assert( hs.is_genuine() );
+	assert( hs.is_next_genuine() );
+	assert( hs.dummy_number() == 0 );
+
+		for (int i = 1 ; i <= 4 ; ++i)
+		{
+		printf("   try %d...\n", i);
+		assert( hs.N_horizons() == capacity );
+		assert( hs.my_N_horizons() == 3 );
+		assert( hs.init_hn() == 42 );
+			assert( hs.get_hn() == 42 );
+			assert( hs.is_genuine() );
+			assert( hs.is_next_genuine() );
+			assert( hs.dummy_number() == 0 );
+		assert( hs.next_hn() == 69 );
+			assert( hs.get_hn() == 69 );
+			assert( hs.is_genuine() );
+			assert( hs.is_next_genuine() );
+			assert( hs.dummy_number() == 0 );
+		assert( hs.next_hn() == 105 );
+			assert( hs.get_hn() == 105 );
+			assert( hs.is_genuine() );
+			assert( ! hs.is_next_genuine() );
+			assert( hs.dummy_number() == 0 );
+		assert( hs.next_hn() == 0 );
+			assert( hs.get_hn() == 0 );
+			assert( ! hs.is_genuine() );
+			assert( ! hs.is_next_genuine() );
+			assert( hs.dummy_number() == 1 );
+		assert( hs.next_hn() == 0 );
+			assert( hs.get_hn() == 0 );
+			assert( ! hs.is_genuine() );
+			assert( ! hs.is_next_genuine() );
+			assert( hs.dummy_number() == 2 );
+		assert( hs.next_hn() == 0 );
+			assert( hs.get_hn() == 0 );
+			assert( ! hs.is_genuine() );
+			assert( ! hs.is_next_genuine() );
+			assert( hs.dummy_number() == 3 );
+		}
+	}
+}
+	  }
diff --git a/src/elliptic/Jacobian.hh b/src/elliptic/Jacobian.hh
index 7d70f5b..0cb43a2 100644
--- a/src/elliptic/Jacobian.hh
+++ b/src/elliptic/Jacobian.hh
@@ -140,9 +140,7 @@ public:
 	// constructor, destructor
 	//
 protected:
-	// ... derived class constructors have a  print_msg_flag
-	//     argument controlling messages about the construction
-	//     and element-setting process
+	// the constructor only uses  ps  to get the size of the matrix
 	Jacobian(patch_system& ps)
 		: ps_(ps),
 		  N_rows_(ps.N_grid_points())
diff --git a/src/elliptic/dense_Jacobian.hh b/src/elliptic/dense_Jacobian.hh
index b56d226..e6227a6 100644
--- a/src/elliptic/dense_Jacobian.hh
+++ b/src/elliptic/dense_Jacobian.hh
@@ -59,6 +59,7 @@ public:
 	// constructor, destructor
 	//
 protected:
+	// the constructor only uses ps to get the size of the matrix
 	// ... print_msg_flag controls messages about data structure setup
 	//     during construction and element-setting process
 	dense_Jacobian(patch_system& ps,
@@ -92,6 +93,7 @@ public:
 
 	// constructor, destructor
 public:
+	// the constructor only uses ps to get the size of the matrix
 	dense_Jacobian__LAPACK(patch_system& ps,
 			       bool print_msg_flag = false);
 	~dense_Jacobian__LAPACK();
diff --git a/src/elliptic/row_sparse_Jacobian.hh b/src/elliptic/row_sparse_Jacobian.hh
index e66a962..5bbc8e7 100644
--- a/src/elliptic/row_sparse_Jacobian.hh
+++ b/src/elliptic/row_sparse_Jacobian.hh
@@ -100,6 +100,7 @@ public:
 	// constructor, destructor
 	//
 protected:
+	// the constructor only uses ps to get the size of the matrix
 	// ... print_msg_flag controls messages about data structure setup
 	//     during construction and element-setting process
 	row_sparse_Jacobian(patch_system& ps,
@@ -190,6 +191,7 @@ public:
 
 	// constructor, destructor
 public:
+	// the constructor only uses ps to get the size of the matrix
 	row_sparse_Jacobian__ILUCG(patch_system& ps,
 				   bool print_msg_flag = false);
 	~row_sparse_Jacobian__ILUCG();
diff --git a/src/gr/expansion.cc b/src/gr/expansion.cc
index 0ca256e..490ff99 100644
--- a/src/gr/expansion.cc
+++ b/src/gr/expansion.cc
@@ -56,7 +56,7 @@ using jtutil::pow4;
 
 namespace {
 void setup_xyz_posns(patch_system& ps, bool print_msg_flag);
-bool interpolate_geometry(patch_system& ps,
+bool interpolate_geometry(patch_system* ps_ptr,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
 			  bool print_msg_flag);
@@ -81,9 +81,12 @@ bool compute_Theta(patch_system& ps,
 //******************************************************************************
 
 //
-// This function computes the LHS function Theta(h), and optionally also
-// its Jacobian coefficients (from which the Jacobian matrix may be
-// computed later).
+// If ps_ptr != NULL, this function computes the LHS function Theta(h),
+// and optionally also its Jacobian coefficients (from which the Jacobian
+// matrix may be computed later).
+//
+// If ps_ptr == NULL, this function does a dummy computation, described
+// below.
 //
 // Inputs (angular gridfns, on ghosted grid):
 // ... defined on ghosted grid
@@ -114,6 +117,11 @@ bool compute_Theta(patch_system& ps,
 //	Theta						# $\Theta = \Theta(h)$
 //
 // Arguments:
+// ps_ptr --> The patch system, or == NULL to do (only) a dummy computation
+//	      in which only the parameter-table setup and a dummy geometry
+//	      interpolator call are done, the latter with the number of
+//	      interpolation points is set to 0 and all the output array
+//	      pointers set to NULL.
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
 // print_msg_flag = true to print status messages,
@@ -124,64 +132,88 @@ bool compute_Theta(patch_system& ps,
 //			   to compute various (gridwise) norms of Theta.
 //
 // Results:
-// This function returns true for a successful computation, or false
-// if the computation failed for any of the following reasons:
-// - the geometry interpolation would need data outside the Cactus grid,
-//   or from an excised region
-//   FIXME: excision isn't implemented yet :(
-// - NaNs are found in h or in the interpolate geometry variables
-// - Theta_D <= 0 (this means the interpolated g_ij aren't positive definite)
+// This function returns
+// * true for a successful computation
+// * true for a dummy computation (ps_ptr == NULL)
+// * false if the computation failed for any of the following reasons:
+//   - the geometry interpolation would need data outside the Cactus grid,
+//     or from an excised region
+//     FIXME: excision isn't implemented yet :(
+//   - NaNs are found in h or in the interpolate geometry variables
+//   - Theta_D <= 0 (this means the interpolated g_ij aren't positive definite)
 //
-bool expansion(patch_system& ps,
+bool expansion(patch_system* ps_ptr,
 	       const struct cactus_grid_info& cgi,
 	       const struct geometry_info& gi,
 	       bool Jacobian_flag /* = false */,
 	       bool print_msg_flag /* = false */,
 	       jtutil::norm<fp>* Theta_norms_ptr /* = NULL */)
 {
+const bool active_flag = (ps_ptr != NULL);
 if (print_msg_flag)
-   then CCTK_VInfo(CCTK_THORNSTRING, "   expansion");
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   %sexpansion",
+		   active_flag ? "" : "dummy ");
 
-// fill in values of all ghosted gridfns in ghost zones
-ps.synchronize();
+if (active_flag)
+   then {
+	//
+	// normal computation
+	//
+
+	// fill in values of all ghosted gridfns in ghost zones
+	ps_ptr->synchronize();
 
-if (gi.check_that_h_is_finite && !h_is_finite(ps, print_msg_flag))
-   then return false;					// *** ERROR RETURN ***
+	if (gi.check_that_h_is_finite
+	    && !h_is_finite(*ps_ptr, print_msg_flag))
+	   then return false;				// *** ERROR RETURN ***
 
-// set up xyz positions of grid points
-setup_xyz_posns(ps, print_msg_flag);
+	// set up xyz positions of grid points
+	setup_xyz_posns(*ps_ptr, print_msg_flag);
+	}
 
 // compute the "geometry" g_ij, K_ij, and partial_k g_ij
 switch	(gi.geometry_method)
 	{
 case geometry__global_interp_from_Cactus_grid:
 case geometry__local_interp_from_Cactus_grid:
-	if (!interpolate_geometry(ps, cgi, gi, print_msg_flag))
+	// this is the only function we call unconditionally; it looks at
+	// ps_ptr (non-NULL vs NULL) to choose a normal vs dummy computation
+	if (!interpolate_geometry(ps_ptr,
+				  cgi, gi,
+				  print_msg_flag))
 	   then return false;				// *** ERROR RETURN ***
-	if (cgi.Cactus_conformal_metric)
-	   then convert_conformal_to_physical(ps, print_msg_flag);
+	if (active_flag && cgi.Cactus_conformal_metric)
+	   then convert_conformal_to_physical(*ps_ptr, print_msg_flag);
 	break;
 
 case geometry__Schwarzschild_EF:
-	Schwarzschild_EF_geometry(ps, gi, print_msg_flag);
+	if (active_flag)
+	   then Schwarzschild_EF_geometry(*ps_ptr, gi, print_msg_flag);
 	break;
+
 default:
 	error_exit(PANIC_EXIT,
 "***** expansion(): unknown gi.geometry_method=(int)%d!\n"
-"                          (this should never happen!)\n"
-,
+"                   (this should never happen!)\n"
+		   ,
 		   int(gi.geometry_method));			/*NOTREACHED*/
 	}
 
-if (gi.check_that_geometry_is_finite
-    && !geometry_is_finite(ps, print_msg_flag))
-   then return false;					// *** ERROR RETURN ***
+if (active_flag)
+   then {
+	if (gi.check_that_geometry_is_finite
+	    && !geometry_is_finite(*ps_ptr, print_msg_flag))
+	   then return false;				// *** ERROR RETURN ***
 
-// compute remaining gridfns --> $\Theta$
-// and optionally also the Jacobian coefficients
-// by algebraic ops and angular finite differencing
-if (!compute_Theta(ps, Jacobian_flag, Theta_norms_ptr, print_msg_flag))
-   then return false;					// *** ERROR RETURN ***
+	// compute remaining gridfns --> $\Theta$
+	// and optionally also the Jacobian coefficients
+	// by algebraic ops and angular finite differencing
+	if (!compute_Theta(*ps_ptr,
+			   Jacobian_flag, Theta_norms_ptr,
+			   print_msg_flag))
+	   then return false;				// *** ERROR RETURN ***
+	}
 
 return true;						// *** NORMAL RETURN ***
 }
@@ -234,7 +266,7 @@ if (print_msg_flag)
 //******************************************************************************
 
 //
-// This function interpolates the Cactus gridfns
+// If ps_ptr != NULL, this function interpolates the Cactus gridfns
 //	gxx...gzz
 //	kxx...kzz
 //	psi			# optional
@@ -245,13 +277,18 @@ if (print_msg_flag)
 //	psi			# optional
 //	partial_d_psi_k		# optional
 // at the nominal-grid trial horizon surface positions given by the
-// global_(x,y,z) angular gridfns.  The psi interpolation is done if and
-// only if the cgi.Cactus_conformal_metric flag is set.  Note that this
-// function ignores the physical-vs-conformal semantics of the gridfns;
-// it just interpolates and takes derivatives of the stored gridfn values.
+// global_(x,y,z) angular gridfns in the patch system *ps_ptr.  The psi
+// interpolation is only done if the cgi.Cactus_conformal_metric flag
+// is set.  Note that this function ignores the physical-vs-conformal
+// semantics of the gridfns; it just interpolates and takes derivatives
+// of the stored gridfn values.
+//
+// If ps_ptr == NULL, this function does (only) the parameter-table
+// setup and a a dummy interpolator call, as described in the comments
+// to  expansion()  above.
 //
-// The interpolation may be either via  CCTK_InterpGridArrays() or via
-// CCTK_InterpLocalUniform() , according to gi.geometry_method.
+// The interpolation may be either via  CCTK_InterpGridArrays()  or via
+//  CCTK_InterpLocalUniform() , according to gi.geometry_method.
 //
 // Inputs (angular gridfns, all on the nominal grid):
 //	global_[xyz]			# xyz positions of grid points
@@ -272,18 +309,23 @@ if (print_msg_flag)
 // This function may also modify the interpolator parameter table.
 //
 // Results:
-// This function returns true for a successful computation, or false
-// if the computation failed because the geometry interpolation would
-// need data outside the Cactus grid, or data from an excised region.
-//
-// FIXME: excision isn't implemented yet :(
+// This function returns
+// * true for a successful computation
+// * true for a dummy computation (ps_ptr == NULL)
+// * false if the computation failed for any of the following reasons:
+//   - the geometry interpolation would need data outside the Cactus grid,
+//     or from an excised region
+//     FIXME: excision isn't implemented yet :(
+// Any other interpolator error return codes are treated as a fatal error.
 //
 namespace {
-bool interpolate_geometry(patch_system& ps,
+bool interpolate_geometry(patch_system* ps_ptr,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_info& gi,
 			  bool print_msg_flag)
 {
+const bool active_flag = (ps_ptr != NULL);
+
 //
 // Implementation Notes:
 //
@@ -293,8 +335,8 @@ bool interpolate_geometry(patch_system& ps,
 // to either include or exclude the psi entries as appropriate.
 //
 // We remember whether or not psi was interpolated on the previous call,
-// and only modify the interpolator parameter table if this changes (and
-// on our first call).
+// and only modify the interpolator parameter table if this changes (or
+// if this is our first call).
 //
 
 const char* const global_local_str
@@ -311,17 +353,20 @@ if (print_msg_flag)
 
 int status;
 
+#define CAST_PTR_OR_NULL(type_,ptr_)	\
+	(ps_ptr == NULL) ? NULL : static_cast<type_>(ptr_)
+
 
 //
 // ***** interpolation points *****
 //
-const int N_interp_points = ps.N_grid_points();
+const int N_interp_points = (ps_ptr == NULL) ? 0 : ps_ptr->N_grid_points();
 const int interp_coords_type_code = CCTK_VARIABLE_REAL;
 const void* const interp_coords[N_GRID_DIMS]
   = {
-    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_x)),
-    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_y)),
-    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_z)),
+    CAST_PTR_OR_NULL(const void*, ps_ptr->gridfn_data(gfns::gfn__global_x)),
+    CAST_PTR_OR_NULL(const void*, ps_ptr->gridfn_data(gfns::gfn__global_y)),
+    CAST_PTR_OR_NULL(const void*, ps_ptr->gridfn_data(gfns::gfn__global_z)),
     };
 
 
@@ -424,45 +469,45 @@ const CCTK_INT operation_codes[]
 
 void* const output_arrays[]
   = {
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_11)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_211)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_311)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_12)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_112)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_212)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_312)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_13)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_113)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_213)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_313)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_22)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_122)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_222)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_322)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_23)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_123)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_223)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_323)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_33)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_133)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_233)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_333)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_11)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_12)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_13)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_22)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_23)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_33)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__psi)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_1)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_2)),
-    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_psi_3)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_11)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_111)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_211)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_311)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_12)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_112)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_212)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_312)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_13)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_113)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_213)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_313)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_22)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_122)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_222)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_322)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_23)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_123)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_223)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_323)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__g_dd_33)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_133)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_233)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_g_dd_333)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_11)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_12)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_13)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_22)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_23)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__K_dd_33)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__psi)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_psi_1)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_psi_2)),
+    CAST_PTR_OR_NULL(void*, ps_ptr->gridfn_data(gfns::gfn__partial_d_psi_3)),
     };
 
 const int N_output_arrays_for_psi = 4;
-const int N_output_arrays_dim = sizeof(output_arrays)
-				/ sizeof(output_arrays[0]);
+const int N_output_arrays_dim
+	= sizeof(output_arrays) / sizeof(output_arrays[0]);
 const int N_output_arrays_use
 	= psi_flag ? N_output_arrays_dim
 		   : N_output_arrays_dim - N_output_arrays_for_psi;
@@ -526,15 +571,12 @@ if (par_table_setup && (psi_flag == par_table_psi_flag))
 //
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling %s interpolator (%d points)",
-		   global_local_str, N_interp_points);
+		   "         calling %s interpolator (%s%d points)",
+		   global_local_str,
+		   (active_flag ? "" : "dummy: "), N_interp_points);
 switch	(gi.geometry_method)
 	{
 case geometry__global_interp_from_Cactus_grid:
-	if (print_msg_flag)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling CCTK_InterpGridArrays() (%d points)",
-			   N_interp_points);
 	status = CCTK_InterpGridArrays(cgi.GH, N_GRID_DIMS,
 				       gi.operator_handle,
 				       gi.param_table_handle,
@@ -549,10 +591,6 @@ case geometry__global_interp_from_Cactus_grid:
 					  output_arrays);
 	break;
 case geometry__local_interp_from_Cactus_grid:
-	if (print_msg_flag)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "         calling CCTK_InterpLocalUniform() (%d points)",
-			   N_interp_points);
 	status = CCTK_InterpLocalUniform(N_GRID_DIMS,
 					 gi.operator_handle,
 					 gi.param_table_handle,
@@ -581,39 +619,39 @@ if (status == CCTK_ERROR_INTERP_POINT_OUTSIDE)
    then {
 	// look in interpolator output table entries
 	// to see *which* point is out-of-range
-	CCTK_INT out_of_range_pt, out_of_range_axis, out_of_range_end;
+	CCTK_INT error_pt, error_axis, error_direction;
 	if (    (Util_TableGetInt(gi.param_table_handle,
-				  &out_of_range_pt,
-				  "out_of_range_pt") < 0)
+				  &error_pt,
+				  "error_pt") < 0)
 	     || (Util_TableGetInt(gi.param_table_handle,
-				  &out_of_range_axis,
-				  "out_of_range_axis") < 0)
+				  &error_axis,
+				  "error_axis") < 0)
 	     || (Util_TableGetInt(gi.param_table_handle,
-				  &out_of_range_end,
-				  "out_of_range_end") < 0)    )
+				  &error_direction,
+				  "error_direction") < 0)    )
 	   then error_exit(ERROR_EXIT,
 "***** interpolate_geometry():\n"
 "        point out of range when interpolating geometry info from 3-D grid!\n"
 "        (one or more points on the trial horizon surface\n"
 "         are outside the grid or too close to the boundary),\n"
-"        but we can't get info about the out-of-range point!\n"
+"        but we can't get info about the bad point(s)!\n"
 "        ==> maybe an interpolator problem?\n");		/*NOTREACHED*/
 
-	assert(out_of_range_pt >= 0);
-	assert(out_of_range_pt < ps.N_grid_points());
+	assert(error_pt >= 0);
+	assert(error_pt < ps_ptr->N_grid_points());
 	const double global_x
-	   = ps.gridfn_data(gfns::gfn__global_x)[out_of_range_pt];
+	   = ps_ptr->gridfn_data(gfns::gfn__global_x)[error_pt];
 	const double global_y
-	   = ps.gridfn_data(gfns::gfn__global_y)[out_of_range_pt];
+	   = ps_ptr->gridfn_data(gfns::gfn__global_y)[error_pt];
 	const double global_z
-	   = ps.gridfn_data(gfns::gfn__global_z)[out_of_range_pt];
+	   = ps_ptr->gridfn_data(gfns::gfn__global_z)[error_pt];
 
-	assert(out_of_range_axis >= 0);
-	assert(out_of_range_axis < N_GRID_DIMS);
-	const char axis = "xyz"[out_of_range_axis];
+	assert(error_axis >= 0);
+	assert(error_axis < N_GRID_DIMS);
+	const char axis = "xyz"[error_axis];
 
-	assert((out_of_range_end == -1) || (out_of_range_end == +1));
-	const char end = (out_of_range_end == -1) ? '-' : '+';
+	assert((error_direction == -1) || (error_direction == +1));
+	const char end = (error_direction == -1) ? '-' : '+';
 
 	if (print_msg_flag)
 	   then {
@@ -636,7 +674,6 @@ if (status < 0)
 ,
 		   status);					/*NOTREACHED*/
 
-
 return true;						// *** NORMAL RETURN ***
 }
 	  }
diff --git a/src/gr/expansion_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index 082387b..805ed97 100644
--- a/src/gr/expansion_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -67,8 +67,8 @@ void add_ghost_zone_Jacobian(const patch_system& ps,
 			     const patch& xp, const ghost_zone& xmgz,
 			     int x_II,
 			     int xm_irho, int xm_isigma);
-bool expansion_Jacobian_dr_FD(patch_system& ps,
-			      Jacobian& Jac,
+bool expansion_Jacobian_dr_FD(patch_system* ps_ptr,
+			      Jacobian* Jac_ptr,
 			      const struct cactus_grid_info& cgi,
 			      const struct geometry_info& gi,
 			      const struct Jacobian_info& Jacobian_info,
@@ -78,32 +78,57 @@ bool expansion_Jacobian_dr_FD(patch_system& ps,
 //******************************************************************************
 
 //
-// This function is a top-level driver to compute the Jacobian matrix
-// J[Theta(h)] of the expansion Theta(h).  It just decodes the
-// Jacobian_compute_method parameter and calls the appropriate subfunction.
+// If ps_ptr != NULL and Jac_ptr != NULL, this function computes the
+// Jacobian matrix J[Theta(h)] of the expansion Theta(h).  We assume
+// that Theta(h) has already been computed.
 //
-// We assume that Theta(h) has already been computed.
+// If ps_ptr == NULL and Jac_ptr == NULL, this function does a dummy
+// computation, in which only any expansion() (and hence geometry
+// interpolator) calls are done, these with the number of interpolation
+// points set to 0 and all the output array pointers set to NULL.
+//
+// It's illegal for one but not both of ps_ptr and Jac_ptr to be NULL.
+//
+// Only some values of  Jacobian_info.Jacobian_compute_method  support
+// the dummy computation.
+//
+// Arguments:
+// ps_ptr --> The patch system, or == NULL to do (only) a dummy computation.
+// Jac_ptr --> The Jacobian, or == NULL to do (only) a dummy computation.
 //
 // Results:
-// This function returns true for a successful computation, or false
-// for failure.  See  expansion()  (in "expansion.cc")
-// for details of the various ways the computation may fail.
+// This function returns
+// * true for a successful computation
+// * true for a dummy computation (ps_ptr == NULL)
+// * false if the computation failed for any of the following reasons
+//   documented for  expansion()  (in "expansion.cc").
 //
-bool expansion_Jacobian(patch_system& ps,
-			Jacobian& Jac,
+bool expansion_Jacobian(patch_system* ps_ptr,
+			Jacobian* Jac_ptr,
 			const struct cactus_grid_info& cgi,
 			const struct geometry_info& gi,
 			const struct Jacobian_info& Jacobian_info,
 			bool print_msg_flag /* = false */)
 {
+const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL);
+
 switch	(Jacobian_info.Jacobian_compute_method)
 	{
 case Jacobian__numerical_perturbation:
-	if (!  expansion_Jacobian_NP(ps, Jac,
-				   cgi, gi, Jacobian_info,
-				   print_msg_flag))
-	   then return false;				// *** ERROR RETURN ***
-	break;
+	if (active_flag)
+	   then {
+		if (! expansion_Jacobian_NP(*ps_ptr, *Jac_ptr,
+					    cgi, gi, Jacobian_info,
+					    print_msg_flag))
+		   then return false;			// *** ERROR RETURN ***
+		break;
+		}
+	   else CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   expansion_Jacobian():\n"
+"        dummy computation isn't supported for\n"
+"        Jacobian_compute_method = \"numerical perturbation\"!");/*NOTREACHED*/
+
 case Jacobian__symbolic_diff:
 	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
@@ -112,15 +137,20 @@ case Jacobian__symbolic_diff:
 "        isn't implemented (yet); reverting to\n"
 "        \"symbolic differentiation with finite diff d/dr\"");
 	// fall through
+
 case Jacobian__symbolic_diff_with_FD_dr:
-	expansion_Jacobian_partial_SD(ps, Jac,
-				      cgi, gi, Jacobian_info,
-				      print_msg_flag);
-	if (! expansion_Jacobian_dr_FD(ps, Jac,
+	if (active_flag)
+	   then expansion_Jacobian_partial_SD(*ps_ptr, *Jac_ptr,
+					      cgi, gi, Jacobian_info,
+					      print_msg_flag);
+	// this function looks at ps_ptr and Jac_ptr (non-NULL vs NULL)
+	// to choose a normal vs dummy computation
+	if (! expansion_Jacobian_dr_FD(ps_ptr, Jac_ptr,
 				       cgi, gi, Jacobian_info,
 				       print_msg_flag))
 	   then return false;				// *** ERROR RETURN ***
 	break;
+
 default:
 	error_exit(PANIC_EXIT,
 "***** expansion_Jacobian():\n"
@@ -129,6 +159,7 @@ default:
 ,
 		   int(Jacobian_info.Jacobian_compute_method));	/*NOTREACHED*/
 	}
+
 return true;						// *** NORMAL RETURN ***
 }
 
@@ -203,8 +234,9 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 
 	const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma);
 	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon;
-	if (! expansion(ps, cgi, ii))
+	if (! expansion(&ps, cgi, ii))
 	   then return false;				// *** ERROR RETURN ***
+
 		for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
 		{
 		patch& xp = ps.ith_patch(xpn);
@@ -226,6 +258,7 @@ ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
 		}
 		}
 		}
+
 	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y;
 	}
 	}
@@ -459,8 +492,16 @@ patch& yp = ye.my_patch();
 //******************************************************************************
 
 //
-// This function sums the d/dr terms into the Jacobian matrix of the
-// expansion Theta(h), computing those terms by finite differencing.
+// If ps_ptr != NULL and Jac_ptr != NULL, this function sums the d/dr
+// terms into the Jacobian matrix of the expansion Theta(h), computing
+// those terms by finite differencing.
+//
+// If ps_ptr == NULL and Jac_ptr == NULL, this function does a dummy
+// computation, in which only any expansion() (and hence geometry
+// interpolator) calls are done, these with the number of interpolation
+// points set to 0 and all the output array pointers set to NULL.
+//
+// It's illegal for one but not both of ps_ptr and Jac_ptr to be NULL.
 //
 // The basic algorithm is that
 //	Jac += diag[ (Theta(h+epsilon) - Theta(h)) / epsilon ]
@@ -473,53 +514,65 @@ patch& yp = ye.my_patch();
 //	Jac += d/dr terms
 //
 // Results:
-// This function returns true for a successful computation, or false
-// if the computation failed because the geometry interpolation would
-// need data outside the Cactus grid, or data from an excised region.
-// FIXME: excision isn't implemented yet :(
+// This function returns
+// * true for a successful computation
+// * true for a dummy computation (ps_ptr == NULL)
+// * false if the computation failed for any of the following reasons:
+//   - the geometry interpolation would need data outside the Cactus grid,
+//     or from an excised region
+//     FIXME: excision isn't implemented yet :(
 //
 namespace {
-bool expansion_Jacobian_dr_FD(patch_system& ps,
-			      Jacobian& Jac,
+bool expansion_Jacobian_dr_FD(patch_system* ps_ptr,
+			      Jacobian* Jac_ptr,
 			      const struct cactus_grid_info& cgi,
 			      const struct geometry_info& gi,
 			      const struct Jacobian_info& Jacobian_info,
 			      bool print_msg_flag)
 {
+const bool active_flag = (ps_ptr != NULL) && (Jac_ptr != NULL);
 if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   horizon Jacobian: d/dr terms (finite diff)");
+		   "   horizon Jacobian: %sd/dr terms (finite diff)",
+		   active_flag ? "" : "dummy ");
 
 const fp epsilon = Jacobian_info.perturbation_amplitude;
 
 // compute Theta(h+epsilon)
-ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
-ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
-if (! expansion(ps, cgi, gi))
+if (active_flag)
+   then {
+	ps_ptr->gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
+	ps_ptr->add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
+	}
+if (! expansion(ps_ptr, cgi, gi))
    then return false;					// *** ERROR RETURN ***
 
-	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
-	{
-	patch& p = ps.ith_patch(pn);
+if (active_flag)
+   then {
+	    for (int pn = 0 ; pn < ps_ptr->N_patches() ; ++pn)
+	    {
+	    patch& p = ps_ptr->ith_patch(pn);
 		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
 		{
 		for (int isigma = p.min_isigma() ;
 		     isigma <= p.max_isigma() ;
 		     ++isigma)
 		{
-		const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma);
+		const int II = ps_ptr->gpn_of_patch_irho_isigma(p, irho,isigma);
 		const fp old_Theta = p.gridfn(gfns::gfn__save_Theta,
 					      irho,isigma);
 		const fp new_Theta = p.gridfn(gfns::gfn__Theta,
 					      irho,isigma);
-		Jac.sum_into_element(II,II, (new_Theta - old_Theta) / epsilon);
+		const fp d_dr_term = (new_Theta - old_Theta) / epsilon;
+		Jac_ptr->sum_into_element(II,II, d_dr_term);
 		}
 		}
-	}
+	    }
 
-// restore h and Theta
-ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
-ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
+	// restore h and Theta
+	ps_ptr->add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+	ps_ptr->gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
+	}
 
 return true;						// *** NORMAL RETURN ***
 }
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index c20a720..3014cd7 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -144,6 +144,7 @@ struct	Jacobian_info
 
 // expansion.cc
 // ... returns true for successful computation,
+// ... returns true for dummy computation (ps_ptr == NULL)
 // ... returns false for failure (eg horizon outside grid);
 //     if (print_msg_flag) then also reports CCTK_VWarn() at the
 //     appropriate warning level
@@ -156,7 +157,7 @@ namespace expansion_warning_levels
 					// ==> skipping nonfinite check
 	       };
 	  }
-bool expansion(patch_system& ps,
+bool expansion(patch_system* ps_ptr,
 	       const struct cactus_grid_info& cgi,
 	       const struct geometry_info& gi,
 	       bool Jacobian_flag = false,
@@ -164,9 +165,11 @@ bool expansion(patch_system& ps,
 	       jtutil::norm<fp>* H_norms_ptr = NULL);
 
 // expansion_Jacobian.cc
-// ... returns true for successful computation, false for failure
-bool expansion_Jacobian(patch_system& ps,
-			Jacobian& Jac,
+// ... returns true for successful computation
+// ... returns true for dummy computation (ps_ptr == NULL && Jac_ptr == NULL)
+// ... returns false for failure
+bool expansion_Jacobian(patch_system* ps_ptr,
+			Jacobian* Jac_ptr,
 			const struct cactus_grid_info& cgi,
 			const struct geometry_info& gi,
 			const struct Jacobian_info& Jacobian_info,
@@ -180,6 +183,8 @@ void Schwarzschild_EF_geometry(patch_system& ps,
 // misc.cc
 enum geometry_method
   decode_geometry_method(const char geometry_method_string[]);
+#ifdef NOT_USED
 bool geometry_method_is_interp(enum geometry_method geometry_method);
+#endif
 enum Jacobian_compute_method
   decode_Jacobian_compute_method(const char Jacobian_compute_method_string[]);
diff --git a/src/gr/misc.cc b/src/gr/misc.cc
index 7cdc018..8a0ced4 100644
--- a/src/gr/misc.cc
+++ b/src/gr/misc.cc
@@ -2,7 +2,9 @@
 // $Header$
 //
 // decode_geometry_method - decode the geometry_method parameter
+#ifdef NOT_USED
 // geomery_method_is_interp - does  enum geometry_method  specify interpolation?
+#endif
 // decode_geometry_method - decode the geometry_method parameter
 //
 
@@ -61,6 +63,7 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 
 //******************************************************************************
 
+#ifdef NOT_USED
 //
 // This function determines whether or not an  enum geometry_method
 // value specifies an interpolation.
@@ -80,6 +83,7 @@ default:
 		   geometry_method);				/*NOTREACHED*/
 	}
 }
+#endif
 
 //******************************************************************************
 
diff --git a/src/include/config.h b/src/include/config.h
index 8a4f84d..a5334e9 100644
--- a/src/include/config.h
+++ b/src/include/config.h
@@ -79,9 +79,9 @@ typedef CCTK_INT integer;
 /**************************************/
 
 /*
- * The following #ifdefs and #defines shouldn't need to be changed;
- * they decode the previous set of Jacobian storage format and linear
- * solver #defines to determine which storage formats we are using.
+ * The following #ifdefs and #defines shouldn't need changing; they
+ * decode the previous set of Jacobian storage format and linear solver
+ * #defines to determine which storage formats we are using.
  */
 #ifdef HAVE_DENSE_JACOBIAN__LAPACK
 #define HAVE_DENSE_JACOBIAN