Replaced all the CCTK_InterpGV() calls by calls to the new interpolation API

CCTK_InterpGridArrays(). Note that you now need to also compile in and activate
a thorn which provides local interpolation operators (eg. LocalInterp from
the CactusBase arrangement).

Also evaluate the return code from CCTK_InterpGridArrays() and print a warning
in case of an error.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@324 89daf98e-ef62-4674-b946-b8ff9de2216c

3 files changed, 197 insertions, 164 deletions

diff --git a/src/AHFinder_dis.F b/src/AHFinder_dis.F
index 374289e..47bcb82 100644
--- a/src/AHFinder_dis.F
+++ b/src/AHFinder_dis.F
@@ -2,9 +2,9 @@
   @file      AHFinder_dis.F
   @date      September 2001
   @author    Miguel Alcubierre
-  @desc 
+  @desc
              Routine to calculate proper distances between horizons.
-  @enddesc 
+  @enddesc
 @@*/
 
 #include "cctk.h"
@@ -22,11 +22,14 @@
       DECLARE_CCTK_FUNCTIONS
 
       integer i,l,m
-      integer np,npoints
-      integer interp_handle,coord_system_handle
-      integer num_arrays,ierror
+      integer np,npoints,ierror
+      integer param_table_handle,interp_handle,coord_system_handle
+      character(30) options_string
 
-      integer, dimension(6) :: in_array_indices
+      CCTK_POINTER, dimension(3) :: interp_coords
+      CCTK_INT, dimension(6) :: in_array_indices
+      CCTK_POINTER, dimension(6) :: out_arrays
+      CCTK_INT, dimension(6) :: out_array_type_codes
 
       CCTK_REAL LEGEN
       CCTK_REAL zero,half,one,two,pi
@@ -40,10 +43,9 @@
       CCTK_REAL, allocatable, dimension(:) :: xa,ya,za
       CCTK_REAL, allocatable, dimension(:) :: txx,tyy,tzz,txy,txz,tyz
 
-      character(len=200) :: disf
+      character(200) disf
       logical disf_exists
 
-
 !     **************************
 !     ***   DEFINE NUMBERS   ***
 !     **************************
@@ -482,51 +484,63 @@
 !     ***   INTERPOLATE METRIC   ***
 !     ******************************
 
-!     Interpolator and coordinate system handle.
-
+!     parameter, local interpolator, and coordinate system handle.
+      param_table_handle = -1
       interp_handle = -1
       coord_system_handle = -1
 
-      if (interpolation_order.eq.1) then
-        call CCTK_InterpHandle(interp_handle,"first-order uniform cartesian")
-      else if (interpolation_order.eq.2) then
-        call CCTK_InterpHandle(interp_handle,"second-order uniform cartesian")
-      else if (interpolation_order.eq.3) then
-        call CCTK_InterpHandle(interp_handle,"third-order uniform cartesian")
+      options_string = "order = " // char(ichar('0') + interpolation_order)
+      call Util_TableCreateFromString (param_table_handle, options_string)
+      if (param_table_handle .lt. 0) then
+        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
       endif
 
+      call CCTK_InterpHandle (interp_handle,"Lagrange polynomial interpolation")
       if (interp_handle.lt.0) then
-        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system")
+        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
       endif
 
       call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
-
-      if (coord_system_handle.lt.0) then
-        call CCTK_WARN(0,"Cannot get handle for interpolation")
+      if (coord_system_handle .lt. 0) then
+        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
       endif
 
+
+!     fill in the input/output arrays for the interpolator
+      interp_coords(1) = CCTK_PointerTo(xa)
+      interp_coords(2) = CCTK_PointerTo(ya)
+      interp_coords(3) = CCTK_PointerTo(za)
+
+      call CCTK_VarIndex (in_array_indices(1), "admbase::gxx")
+      call CCTK_VarIndex (in_array_indices(2), "admbase::gyy")
+      call CCTK_VarIndex (in_array_indices(3), "admbase::gzz")
+      call CCTK_VarIndex (in_array_indices(4), "admbase::gxy")
+      call CCTK_VarIndex (in_array_indices(5), "admbase::gxz")
+      call CCTK_VarIndex (in_array_indices(6), "admbase::gyz")
+
+      out_arrays(1) = CCTK_PointerTo(txx)
+      out_arrays(2) = CCTK_PointerTo(tyy)
+      out_arrays(3) = CCTK_PointerTo(tzz)
+      out_arrays(4) = CCTK_PointerTo(txy)
+      out_arrays(5) = CCTK_PointerTo(txz)
+      out_arrays(6) = CCTK_PointerTo(tyz)
+
+      out_array_type_codes = CCTK_VARIABLE_REAL
+
+
 !     Interpolation.
 
-      call CCTK_VarIndex(in_array_indices(1), "admbase::gxx")
-      call CCTK_VarIndex(in_array_indices(2), "admbase::gyy")
-      call CCTK_VarIndex(in_array_indices(3), "admbase::gzz")
-      call CCTK_VarIndex(in_array_indices(4), "admbase::gxy")
-      call CCTK_VarIndex(in_array_indices(5), "admbase::gxz")
-      call CCTK_VarIndex(in_array_indices(6), "admbase::gyz")
-
-      num_arrays = 6
-      call CCTK_InterpGV(ierror, cctkGH, interp_handle, coord_system_handle,
-     .                   npoints, num_arrays, num_arrays,
-     .                   xa, ya, za,
-     .                   CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                   CCTK_VARIABLE_REAL,
-     .                   in_array_indices(1), in_array_indices(2),
-     .                   in_array_indices(3), in_array_indices(4),
-     .                   in_array_indices(5), in_array_indices(6),
-     .                   txx, tyy, tzz, txy, txz, tyz,
-     .                   CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                   CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                   CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL)
+      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+     .                            param_table_handle, coord_system_handle,
+     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
+     .                            6, in_array_indices,
+     .                            6, out_array_type_codes, out_arrays)
+      if (ierror < 0) then
+        call CCTK_WARN (1, "interpolator call returned an error code");
+      endif
+
+!     release parameter table
+      call Util_TableDestroy (ierror, param_table_handle)
 
 
 !     ******************************
diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index d66c214..8005843 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -2,7 +2,7 @@
   @file      AHFinder_gau.F
   @date      October 1998
   @author    Miguel Alcubierre
-  @desc 
+  @desc
              Find gaussian curvature of surface.  The gaussian
              curvature is defined as R/2, where R is the Ricci
              scalar of the induced 2-geometry of the surface.
@@ -10,7 +10,7 @@
              As an extra "goody", this routine also integrates
              the equatorial circumference of the horizon and
              two polar circumferences at phi=0 and phi=pi/2.
-  @enddesc 
+  @enddesc
 @@*/
 
 #include "cctk.h"
@@ -32,11 +32,13 @@
 
       integer i,j,k,l,m,n,p
       integer npoints
-      integer interp_handle,coord_system_handle,sum_handle
+      integer param_table_handle,interp_handle,coord_system_handle,sum_handle
       integer ierror
 
       CCTK_INT rerror,error1,error2
 
+      character(30) options_string
+
       CCTK_REAL LEGEN
       CCTK_REAL theta,phi,xp,yp,zp,rp
       CCTK_REAL cost,sint,cosp,sinp
@@ -53,8 +55,10 @@
       CCTK_REAL, dimension(2,2,2)   :: d1a,gammad,gamma
       CCTK_REAL, dimension(2,2,2,2) :: d2a,gamma2
 
-      integer    num_arrays
-      integer,   dimension(7)   :: in_array_indices
+      CCTK_POINTER, dimension(3) :: interp_coords
+      CCTK_INT, dimension(7) :: in_array_indices
+      CCTK_POINTER, dimension(7) :: out_arrays
+      CCTK_INT, dimension(7) :: out_array_type_codes
 
       CCTK_REAL, allocatable, dimension(:,:) :: rr,xa,ya,za
       CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz
@@ -205,7 +209,7 @@
 
       call CCTK_ReductionArrayHandle(sum_handle,"sum")
 
-      if (sum_handle.lt.0) then 
+      if (sum_handle.lt.0) then
         call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??")
       end if
 
@@ -271,7 +275,7 @@
       g12 = zero
 
       gaussian = zero
-      
+
 
 !     *********************************
 !     ***   INITIALIZE PARAMETERS   ***
@@ -630,16 +634,16 @@
 !     Reduce the errors across processors (all processors must
 !     know about this since all will participate on the interpolation
 !     below).
-      
-      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle, 
-     .   error1,rerror,CCTK_VARIABLE_INT) 
+
+      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
+     .   error1,rerror,CCTK_VARIABLE_INT)
 
       if (ierror.ne.0) then
          call CCTK_WARN(1,"Reduction failed!")
       end if
       error1 = rerror
 
-      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle, 
+      call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
      .   error2,rerror,CCTK_VARIABLE_INT)
 
       if (ierror.ne.0) then
@@ -670,30 +674,32 @@
 !     should be done in the routine AHFinder_int.F).  I need to fix
 !     this sometime soon! (October 2000).
 
-!     Interpolator and coordinate system handle.
-
+!     parameter, local interpolator, and coordinate system handle.
+      param_table_handle = -1
       interp_handle = -1
       coord_system_handle = -1
 
-      if (interpolation_order .eq. 1) then
-        call CCTK_InterpHandle (interp_handle, "first-order uniform cartesian")
-      else if (interpolation_order .eq. 2) then
-        call CCTK_InterpHandle (interp_handle, "second-order uniform cartesian")
-      else if (interpolation_order .eq. 3) then
-        call CCTK_InterpHandle (interp_handle, "third-order uniform cartesian")
+      options_string = "order = " // char(ichar('0') + interpolation_order)
+      call Util_TableCreateFromString (param_table_handle, options_string)
+      if (param_table_handle .lt. 0) then
+        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
       endif
 
+      call CCTK_InterpHandle (interp_handle,"Lagrange polynomial interpolation")
       if (interp_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
+        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
       endif
 
       call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
-
       if (coord_system_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
+        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
       endif
 
-!     Interpolation.
+
+!     fill in the input/output arrays for the interpolator
+      interp_coords(1) = CCTK_PointerTo(xa)
+      interp_coords(2) = CCTK_PointerTo(ya)
+      interp_coords(3) = CCTK_PointerTo(za)
 
       call CCTK_VarIndex (in_array_indices(1), "admbase::gxx")
       call CCTK_VarIndex (in_array_indices(2), "admbase::gyy")
@@ -703,21 +709,29 @@
       call CCTK_VarIndex (in_array_indices(6), "admbase::gyz")
       call CCTK_VarIndex (in_array_indices(7), "ahfinder::ahfgauss")
 
-      num_arrays = 7
-      call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
-     .                    npoints, num_arrays, num_arrays,
-     .                    xa, ya, za,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL,
-     .                    in_array_indices(1), in_array_indices(2),
-     .                    in_array_indices(3), in_array_indices(4),
-     .                    in_array_indices(5), in_array_indices(6),
-     .                    in_array_indices(7),
-     .                    txx, tyy, tzz, txy, txz, tyz, gaussian,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL)
+      out_arrays(1) = CCTK_PointerTo(txx)
+      out_arrays(2) = CCTK_PointerTo(tyy)
+      out_arrays(3) = CCTK_PointerTo(tzz)
+      out_arrays(4) = CCTK_PointerTo(txy)
+      out_arrays(5) = CCTK_PointerTo(txz)
+      out_arrays(6) = CCTK_PointerTo(tyz)
+      out_arrays(7) = CCTK_PointerTo(gaussian)
+
+      out_array_type_codes = CCTK_VARIABLE_REAL
+
+
+!     Interpolation.
+      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+     .                            param_table_handle, coord_system_handle,
+     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
+     .                            7, in_array_indices,
+     .                            7, out_array_type_codes, out_arrays)
+      if (ierror < 0) then
+        call CCTK_WARN (1, "interpolator call returned an error code");
+      endif
+
+!     release parameter table
+      call Util_TableDestroy (ierror, param_table_handle)
 
 
 !     *********************************
diff --git a/src/AHFinder_int.F b/src/AHFinder_int.F
index b3feabe..060834a 100644
--- a/src/AHFinder_int.F
+++ b/src/AHFinder_int.F
@@ -3,7 +3,7 @@
   @file      AHFinder_int.F
   @date      April 1998
   @author    Miguel Alcubierre
-  @desc 
+  @desc
              Find surface integrals. The integrals are done
              in parallel, but the number of processors is
              assumed to be either the square of an integer
@@ -34,10 +34,16 @@
       integer auxi
       integer ierror
 
-      integer interp_handle,coord_system_handle,reduction_handle
+      integer interp_handle,coord_system_handle
+      integer reduction_handle,param_table_handle
+      character(30) options_string
       CCTK_INT red_tmp
 
-      integer, dimension(9) :: in_array_indices
+      CCTK_POINTER, dimension(3) :: interp_coords
+      CCTK_INT, dimension(9) :: in_array_indices
+      CCTK_POINTER, dimension(9) :: out_arrays
+      CCTK_INT, dimension(9) :: out_array_type_codes
+
 
       CCTK_REAL LEGEN
       CCTK_REAL xp,yp,zp,rp
@@ -50,7 +56,6 @@
       CCTK_REAL sina,cosa,intw,grad,sigma,h0
       CCTK_REAL theta,phi
       CCTK_REAL aux,aux1,aux2
-      INTEGER   num_arrays
 
       CCTK_REAL tempv(0:lmax),tempm(lmax,lmax)
 
@@ -166,7 +171,7 @@
 
       call CCTK_ReductionArrayHandle(reduction_handle,"sum")
 
-      if (reduction_handle.lt.0) then 
+      if (reduction_handle.lt.0) then
         call CCTK_WARN(1,"Cannot get handle for sum reduction ! Forgot to activate an implementation providing reduction operators ??")
       end if
 
@@ -577,7 +582,7 @@
       interror1 = red_tmp
 
       call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .   interror2,red_tmp,CCTK_VARIABLE_INT) 
+     .   interror2,red_tmp,CCTK_VARIABLE_INT)
       if (ierror.ne.0) then
          call CCTK_WARN(1,"Reduction of norm failed!")
       end if
@@ -608,29 +613,32 @@
 
       npoints = (l_ntheta+1)*(l_nphi+1)
 
-!     Interpolator and coordinate system handle.
-
+!     parameter, local interpolator, and coordinate system handle.
+      param_table_handle = -1
       interp_handle = -1
       coord_system_handle = -1
 
-      if (interpolation_order .eq. 1) then
-        call CCTK_InterpHandle (interp_handle, "first-order uniform cartesian")
-      else if (interpolation_order .eq. 2) then
-        call CCTK_InterpHandle (interp_handle, "second-order uniform cartesian")
-      else if (interpolation_order .eq. 3) then
-        call CCTK_InterpHandle (interp_handle, "third-order uniform cartesian")
+      options_string = "order = " // char(ichar('0') + interpolation_order)
+      call Util_TableCreateFromString (param_table_handle, options_string)
+      if (param_table_handle .lt. 0) then
+        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
       endif
 
+      call CCTK_InterpHandle (interp_handle,"Lagrange polynomial interpolation")
       if (interp_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
+        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
       endif
 
       call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
-
       if (coord_system_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
+        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
       endif
 
+!     fill in the input/output arrays for the interpolator
+      interp_coords(1) = CCTK_PointerTo(xa)
+      interp_coords(2) = CCTK_PointerTo(ya)
+      interp_coords(3) = CCTK_PointerTo(za)
+
       call CCTK_VarIndex (in_array_indices(1), "admbase::gxx")
       call CCTK_VarIndex (in_array_indices(2), "admbase::gyy")
       call CCTK_VarIndex (in_array_indices(3), "admbase::gzz")
@@ -641,69 +649,66 @@
       call CCTK_VarIndex (in_array_indices(8), "ahfinder::ahmask")
       call CCTK_VarIndex (in_array_indices(9), "ahfinder::ahfgradn")
 
+      out_arrays(1) = CCTK_PointerTo(txx)
+      out_arrays(2) = CCTK_PointerTo(tyy)
+      out_arrays(3) = CCTK_PointerTo(tzz)
+      out_arrays(4) = CCTK_PointerTo(txy)
+      out_arrays(5) = CCTK_PointerTo(txz)
+      out_arrays(6) = CCTK_PointerTo(tyz)
+      out_arrays(7) = CCTK_PointerTo(exp)
+      out_arrays(8) = CCTK_PointerTo(intmask)
+      out_arrays(9) = CCTK_PointerTo(gradn)
+
+      out_array_type_codes = CCTK_VARIABLE_REAL
+
+
 !     For minimization we need the interpolated metric,
 !     the expansion and the mask.
 
       if (.not.flow) then
 
-         num_arrays = 8
-         call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
-     .                       npoints, num_arrays, num_arrays,
-     .                       xa, ya, za,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL,
-     .                       in_array_indices(1), in_array_indices(2),
-     .                       in_array_indices(3), in_array_indices(4),
-     .                       in_array_indices(5), in_array_indices(6),
-     .                       in_array_indices(7), in_array_indices(8),
-     .                       txx, tyy, tzz, txy, txz, tyz, exp, intmask,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL)
+         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+     .                               param_table_handle, coord_system_handle,
+     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
+     .                               8, in_array_indices,
+     .                               8, out_array_type_codes, out_arrays)
+      if (ierror < 0) then
+        call CCTK_WARN (1, "interpolator call returned an error code");
+      endif
+
 
 !     For N flow, we need the interpolated metric, the expansion,
 !     the norm of the gradient of the horizon function, and the mask.
 
       else if (nw.ne.zero) then
 
-         num_arrays = 9
-         call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
-     .                       npoints, num_arrays, num_arrays,
-     .                       xa, ya, za,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL,
-     .                       in_array_indices(1), in_array_indices(2),
-     .                       in_array_indices(3), in_array_indices(4),
-     .                       in_array_indices(5), in_array_indices(6),
-     .                       in_array_indices(7), in_array_indices(8),
-     .                       in_array_indices(9),
-     .                       txx, tyy, tzz, txy, txz, tyz, exp, intmask, gradn,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL)
+         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+     .                               param_table_handle, coord_system_handle,
+     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
+     .                               9, in_array_indices,
+     .                               9, out_array_type_codes, out_arrays)
+         if (ierror < 0) then
+           call CCTK_WARN (1, "interpolator call returned an error code");
+         endif
 
 !     For H or C flows, we need the interpolated expansion, the
 !     norm of the gradient of the horizon function, and the mask.
 
       else
 
-         num_arrays = 3
-         call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
-     .                       npoints, num_arrays, num_arrays,
-     .                       xa, ya, za,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL,
-     .                       in_array_indices(7), in_array_indices(8),
-     .                       in_array_indices(9),
-     .                       exp, intmask, gradn,
-     .                       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     .                       CCTK_VARIABLE_REAL)
+         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
+     .                               param_table_handle, coord_system_handle,
+     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
+     .                               3, in_array_indices(7),
+     .                               3, out_array_type_codes(7), out_arrays(7))
+         if (ierror < 0) then
+           call CCTK_WARN (1, "interpolator call returned an error code");
+         endif
 
       end if
 
+!     release parameter table
+      call Util_TableDestroy (ierror, param_table_handle)
 
 !     ***************************************
 !     ***   CHECK IF WE ARE INSIDE MASK   ***
@@ -937,27 +942,27 @@
 !        Reduce integrals.
 
          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .      intarea,aux,CCTK_VARIABLE_REAL) 
+     .      intarea,aux,CCTK_VARIABLE_REAL)
          intarea=aux
 
          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .      intexp,aux,CCTK_VARIABLE_REAL) 
+     .      intexp,aux,CCTK_VARIABLE_REAL)
          intexp=aux
 
          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .      intexp2,aux,CCTK_VARIABLE_REAL) 
+     .      intexp2,aux,CCTK_VARIABLE_REAL)
          intexp2=aux
 
          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .      intexpdel2,aux,CCTK_VARIABLE_REAL) 
+     .      intexpdel2,aux,CCTK_VARIABLE_REAL)
          intexpdel2=aux
 
          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .      inside_min_neg_count,red_tmp,CCTK_VARIABLE_INT) 
+     .      inside_min_neg_count,red_tmp,CCTK_VARIABLE_INT)
          inside_min_neg_count=red_tmp
 
          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,reduction_handle,
-     .      inside_min_count,red_tmp,CCTK_VARIABLE_INT) 
+     .      inside_min_count,red_tmp,CCTK_VARIABLE_INT)
          inside_min_count=red_tmp
 
       end if
@@ -1033,9 +1038,9 @@
                      zw = za(i,j) - zc
 
 !                    Calculate metric on surface.
-             
-                     det = txx(i,j)*tyy(i,j)*tzz(i,j) 
-     .                   + two*txy(i,j)*txz(i,j)*tyz(i,j) 
+
+                     det = txx(i,j)*tyy(i,j)*tzz(i,j)
+     .                   + two*txy(i,j)*txz(i,j)*tyz(i,j)
      .                   - txx(i,j)*tyz(i,j)**2 - tyy(i,j)*txz(i,j)**2
      .                   - tzz(i,j)*txy(i,j)**2
 
@@ -1044,7 +1049,7 @@
                      gupij(1,1) = idet*(tyy(i,j)*tzz(i,j)-tyz(i,j)**2)
                      gupij(2,2) = idet*(txx(i,j)*tzz(i,j)-txz(i,j)**2)
                      gupij(3,3) = idet*(txx(i,j)*tyy(i,j)-txy(i,j)**2)
-                  
+
                      gupij(1,2) = idet*(txz(i,j)*tyz(i,j)-tzz(i,j)*txy(i,j))
                      gupij(2,1) = gupij(1,2)
                      gupij(1,3) = idet*(txy(i,j)*tyz(i,j)-tyy(i,j)*txz(i,j))
@@ -1122,21 +1127,21 @@
 
          auxi = lmax+1
 
-	 if (hw.ne.zero) then
+         if (hw.ne.zero) then
             call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
-     .         hflow0,tempv,auxi,CCTK_VARIABLE_REAL) 
+     .         hflow0,tempv,auxi,CCTK_VARIABLE_REAL)
             hflow0 = tempv
          end if
 
-         if (cw.ne.zero) then      
+         if (cw.ne.zero) then
             call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
-     .         cflow0,tempv,auxi,CCTK_VARIABLE_REAL) 
+     .         cflow0,tempv,auxi,CCTK_VARIABLE_REAL)
             cflow0 = tempv
          end if
 
          if (nw.ne.zero) then
             call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,reduction_handle,
-     .         nflow0,tempv,auxi,CCTK_VARIABLE_REAL) 
+     .         nflow0,tempv,auxi,CCTK_VARIABLE_REAL)
             nflow0 = tempv
          end if
 
@@ -1146,41 +1151,41 @@
 
             auxi = lmax*lmax
 
-	    if (hw.ne.zero) then
+            if (hw.ne.zero) then
                call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
-     .            reduction_handle,hflowc,tempm,auxi,CCTK_VARIABLE_REAL) 
+     .            reduction_handle,hflowc,tempm,auxi,CCTK_VARIABLE_REAL)
                hflowc = tempm
-	    end if
+            end if
 
             if (cw.ne.zero) then
                call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
-     .            reduction_handle,cflowc,tempm,auxi,CCTK_VARIABLE_REAL) 
+     .            reduction_handle,cflowc,tempm,auxi,CCTK_VARIABLE_REAL)
                cflowc = tempm
             end if
 
             if (nw.ne.zero) then
                call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
-     .            reduction_handle,nflowc,tempm,auxi,CCTK_VARIABLE_REAL) 
+     .            reduction_handle,nflowc,tempm,auxi,CCTK_VARIABLE_REAL)
                nflowc = tempm
             end if
 
             if (.not.refy) then
 
-	       if (hw.ne.zero) then
+               if (hw.ne.zero) then
                   call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
-     .               reduction_handle,hflows,tempm,auxi,CCTK_VARIABLE_REAL) 
+     .               reduction_handle,hflows,tempm,auxi,CCTK_VARIABLE_REAL)
                   hflows = tempm
-	       end if
+               end if
 
                if (cw.ne.zero) then
                   call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
-     .               reduction_handle,cflows,tempm,auxi,CCTK_VARIABLE_REAL) 
+     .               reduction_handle,cflows,tempm,auxi,CCTK_VARIABLE_REAL)
                   cflows = tempm
                end if
 
                if (nw.ne.zero) then
                   call CCTK_ReduceLocalArray1D(ierror,cctkGH,-1,
-     .               reduction_handle,nflows,tempm,auxi,CCTK_VARIABLE_REAL) 
+     .               reduction_handle,nflows,tempm,auxi,CCTK_VARIABLE_REAL)
                   nflows = tempm
                end if
 
@@ -1236,7 +1241,7 @@
          if (flow) then
 
 !           Axisymmetric.
-	
+
             hflow0 = intw*hflow0
             cflow0 = intw*cflow0
             nflow0 = intw*nflow0