<p><b>xylar@lanl.gov</b> 2010-04-26 16:58:02 -0600 (Mon, 26 Apr 2010)</p><p>Branch Commit<br>
<br>
module_interp_tests: Adding test functions to make sure the RBF interpolation code works properly.  So far, so good...<br>
<br>
removed references to dmpar_abort in anticipation of a dmpar_global_abort function to come<br>
<br>
Changes to the Makefiles are copied from the trunk<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/Makefile
===================================================================
--- branches/ocean_projects/IBinterp/mpas/Makefile        2010-04-26 18:52:09 UTC (rev 212)
+++ branches/ocean_projects/IBinterp/mpas/Makefile        2010-04-26 22:58:02 UTC (rev 213)
@@ -1,7 +1,7 @@
 #MODEL_FORMULATION = -DNCAR_FORMULATION
 MODEL_FORMULATION = -DLANL_FORMULATION
 
-EXPAND_LEVELS = -DEXPAND_LEVELS=26
+#EXPAND_LEVELS = -DEXPAND_LEVELS=26
 #FILE_OFFSET = -DOFFSET64BIT
 
 #########################

Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F        2010-04-26 18:52:09 UTC (rev 212)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F        2010-04-26 22:58:02 UTC (rev 213)
@@ -18,6 +18,7 @@
    use time_integration
    use IB_interpolation
    use vector_reconstruction
+   use interp_tests
 
 
    implicit none
@@ -29,6 +30,9 @@
    call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
 
    call ibInterp_initialize(mesh)
+
+   call doTests(mesh)
+
    call init_reconstruct(mesh)
 
 ! mrp 100316 In order for this to work, we need to pass domain % dminfo as an 

Modified: branches/ocean_projects/IBinterp/mpas/src/operators/Makefile
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/Makefile        2010-04-26 18:52:09 UTC (rev 212)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/Makefile        2010-04-26 22:58:02 UTC (rev 213)
@@ -1,6 +1,6 @@
 .SUFFIXES: .F .o
 
-OBJS = module_IB_interpolation.o module_vector_reconstruction.o
+OBJS = module_IB_interpolation.o module_interp_tests.o module_vector_reconstruction.o
 
 all: operators
 

Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F        2010-04-26 18:52:09 UTC (rev 212)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F        2010-04-26 22:58:02 UTC (rev 213)
@@ -9,10 +9,6 @@
    save
 
 
-   integer, public, parameter :: &amp;
-     MPAS_i4 = selected_int_kind(6), &amp;
-     MPAS_r8 = selected_real_kind(13)
-
    public :: ibInterp_initialize, ibInterp_get2DRBFCoefficients, &amp;
      ibInterp_evaluate2DRBFWithDerivs, &amp;
      ibInterp_computeScalarRBFAndConstantCoefficients, &amp;
@@ -23,7 +19,7 @@
      ibInterp_computeVectorFreeSlipRBFCoefficients, &amp;
      ibInterp_computeVectorFreeSlipPlanarRBFCoefficients
 
-   integer(MPAS_i4), public, parameter :: &amp;
+   integer, public, parameter :: &amp;
      kMPAS_IBRBFPolynomialConstant = 0, &amp;
      kMPAS_IBRBFPolynomialLinear = 1
 
@@ -50,15 +46,15 @@
 
     type (grid_meta), intent(inout) :: grid 
 
-    integer(MPAS_i4) :: nCells, nEdges
-    integer(MPAS_i4), dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
-    integer(MPAS_i4) :: iEdge, iCell, cell1, cell2
-    real(MPAS_r8), dimension(:), pointer :: xCell, yCell, zCell
-    real(MPAS_r8), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &amp;
+    integer :: nCells, nEdges
+    integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+    integer :: iEdge, iCell, cell1, cell2
+    real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+    real(kind=RKIND), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &amp;
       edgeLocations
-    real(MPAS_r8), dimension(:,:,:), pointer :: cellTangentPlane
-    real(MPAS_r8), dimension(3) :: xHatPlane, yHatPlane, rHat
-    real(MPAS_r8) :: normalDotRHat
+    real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
+    real(kind=RKIND), dimension(3) :: xHatPlane, yHatPlane, rHat
+    real(kind=RKIND) :: normalDotRHat
 
     xCell       =&gt; grid % xCell % array
     yCell       =&gt; grid % yCell % array
@@ -115,25 +111,24 @@
   end subroutine ibInterp_initialize
 
   subroutine ibInterp_get2DRBFCoefficients(pointCount, coeffCount, dimensions, points, &amp;
-    fieldValues, polynomialFlag, alpha, coefficients, dminfo)
+    fieldValues, polynomialFlag, alpha, coefficients)
  
-    integer(MPAS_i4), intent(in) :: pointCount, coeffCount, dimensions
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: points
-    real(MPAS_r8), dimension(pointCount), intent(in) :: fieldValues
-    integer(MPAS_i4), intent(in) :: polynomialFlag
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(coeffCount), intent(out) :: coefficients
-    type (dm_info), intent(in) :: dminfo
+    integer, intent(in) :: pointCount, coeffCount, dimensions
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: points
+    real(kind=RKIND), dimension(pointCount), intent(in) :: fieldValues
+    integer, intent(in) :: polynomialFlag
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(coeffCount), intent(out) :: coefficients
 
-    integer(MPAS_i4) :: i, j, matrixSize
-    real(MPAS_r8), dimension(pointCount+3,pointCount+3) :: matrix
-    real(MPAS_r8), dimension(pointCount+3) :: rhs
-    integer(MPAS_i4), dimension(pointCount+3) :: pivotIndices
-    real(MPAS_r8) :: rSquared
+    integer :: i, j, matrixSize
+    real(kind=RKIND), dimension(pointCount+3,pointCount+3) :: matrix
+    real(kind=RKIND), dimension(pointCount+3) :: rhs
+    integer, dimension(pointCount+3) :: pivotIndices
+    real(kind=RKIND) :: rSquared
 
-    coefficients = 0.0_MPAS_r8
-    matrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
+    coefficients = 0.0
+    matrix = 0.0
+    rhs = 0.0
 
     rhs(1:pointCount) = fieldValues
 
@@ -147,19 +142,19 @@
     if(polynomialFlag == kMPAS_IBRBFPolynomialConstant) then
       matrixSize = pointCount+1
       do j=1,pointCount
-        matrix(pointCount+1,j) = 1.0_MPAS_r8
-        matrix(j,pointCount+1) = 1.0_MPAS_r8
+        matrix(pointCount+1,j) = 1.0
+        matrix(j,pointCount+1) = 1.0
       end do
     else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
       matrixSize = pointCount+3
       do j=1,pointCount
-        matrix(pointCount+1,j) = 1.0_MPAS_r8
+        matrix(pointCount+1,j) = 1.0
         matrix(pointCount+2,j) = points(j,1)
         matrix(pointCount+3,j) = points(j,2)
         matrix(j,pointCount+1:pointCount+3) = matrix(pointCount+1:pointCount+3, j)
       end do
     else
-      call exit_MPAS(&quot;ibInterp_get2DRBFCoefficients: not a valid polynomialFlag&quot;, dminfo)
+      call exit_MPAS(&quot;ibInterp_get2DRBFCoefficients: not a valid polynomialFlag&quot;)
     end if
 
     call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &amp;
@@ -168,21 +163,20 @@
   end subroutine ibInterp_get2DRBFCoefficients
 
   subroutine ibInterp_evaluate2DRBFWithDerivs(fieldCount, coeffCount, pointCount, &amp;
-    coefficients, evaluationPoint, points, polynomialFlag, alpha, derivs, dminfo)
-    integer(MPAS_i4), intent(in) :: fieldCount, coeffCount, pointCount
-    real(MPAS_r8), dimension(coeffCount, fieldCount), intent(in) :: coefficients
-    real(MPAS_r8), dimension(2), intent(in) :: evaluationPoint
-    real(MPAS_r8), dimension(pointCount,2), intent(in) :: points
-    integer(MPAS_i4), intent(in) :: polynomialFlag
-    real(MPAS_r8), intent(in) :: alpha
-    type (dm_info), intent(in) :: dminfo
+    coefficients, evaluationPoint, points, polynomialFlag, alpha, derivs)
+    integer, intent(in) :: fieldCount, coeffCount, pointCount
+    real(kind=RKIND), dimension(coeffCount, fieldCount), intent(in) :: coefficients
+    real(kind=RKIND), dimension(2), intent(in) :: evaluationPoint
+    real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
+    integer, intent(in) :: polynomialFlag
+    real(kind=RKIND), intent(in) :: alpha
 
-    real(MPAS_r8), dimension(6,fieldCount), intent(out) :: derivs
+    real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
 
-    integer(MPAS_i4) :: pointIndex
-    real(MPAS_r8) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
+    integer :: pointIndex
+    real(kind=RKIND) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
 
-    derivs = 0.0_MPAS_r8
+    derivs = 0.0
     do pointIndex = 1, pointCount
       x = (evaluationPoint(1) - points(pointIndex,1))
       y = (evaluationPoint(2) - points(pointIndex,2))
@@ -217,30 +211,29 @@
       derivs(2,:) = derivs(2,:) + coefficients(pointCount+2,:)
       derivs(3,:) = derivs(3,:) + coefficients(pointCount+3,:)
     else
-      call exit_MPAS(&quot;ibInterp_get2DRBFCoefficients: not a valid polynomialFlag&quot;, dminfo)
+      call exit_MPAS(&quot;ibInterp_get2DRBFCoefficients: not a valid polynomialFlag&quot;)
     end if
   end subroutine ibInterp_evaluate2DRBFWithDerivs
 
   subroutine ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &amp;
     sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
-    alpha, dirichletCoefficients, neumannCoefficients, dminfo)
+    alpha, dirichletCoefficients, neumannCoefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isInterface
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount), intent(out) :: &amp;
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount), intent(out) :: &amp;
       dirichletCoefficients, neumannCoefficients
-    type (dm_info), intent(in) :: dminfo
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
-    real(MPAS_r8), dimension(:), pointer :: rhs, rhsCopy, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+    real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+1 !! 1 extra space for constant 
 
@@ -251,11 +244,11 @@
     allocate(coeffs(matrixSize))  
     allocate(pivotIndices(matrixSize))  
 
-    dirichletMatrix = 0.0_MPAS_r8
-    neumannMatrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    rhsCopy = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    dirichletMatrix = 0.0
+    neumannMatrix = 0.0
+    rhs = 0.0
+    rhsCopy = 0.0
+    coeffs = 0.0
 
     call setUpScalarRBFMatrixAndRHS(pointCount, &amp;
       sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
@@ -263,9 +256,9 @@
       neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
 
     do i = 1, pointCount
-      dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+      dirichletMatrix(i,pointCount+1) = 1.0
       if(isInterface(i)) then
-        neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+        neumannMatrix(i,pointCount+1) = 0.0
       else
         neumannMatrix(i,pointCount+1) = dirichletMatrix(i,pointCount+1)
       end if
@@ -273,7 +266,7 @@
       neumannMatrix(pointCount+1,i) = neumannMatrix(i,pointCount+1)
     end do
 
-    rhs(pointCount+1) = 1.0_MPAS_r8
+    rhs(pointCount+1) = 1.0
 
     ! solve each linear system
     rhsCopy = rhs
@@ -294,25 +287,24 @@
 
   subroutine ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &amp;
     sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
-    alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+    alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isInterface
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(2,3) :: planeBasisVectors
-    real(MPAS_r8), dimension(pointCount), intent(out) :: &amp;
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+    real(kind=RKIND), dimension(pointCount), intent(out) :: &amp;
       dirichletCoefficients, neumannCoefficients
-    type (dm_info), intent(in) :: dminfo
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
-    real(MPAS_r8), dimension(:), pointer :: rhs, rhsCopy, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+    real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+3 !! 3 extra space for constant and 2 planar dimensions
 
@@ -323,11 +315,11 @@
     allocate(coeffs(matrixSize))  
     allocate(pivotIndices(matrixSize))  
 
-    dirichletMatrix = 0.0_MPAS_r8
-    neumannMatrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    rhsCopy = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    dirichletMatrix = 0.0
+    neumannMatrix = 0.0
+    rhs = 0.0
+    rhsCopy = 0.0
+    coeffs = 0.0
 
     call setUpScalarRBFMatrixAndRHS(pointCount, &amp;
       sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
@@ -335,11 +327,11 @@
       neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
 
     do i = 1, pointCount
-      dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+      dirichletMatrix(i,pointCount+1) = 1.0
       dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVectors(1,:))
       dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVectors(2,:))
       if(isInterface(i)) then
-        neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+        neumannMatrix(i,pointCount+1) = 0.0
         neumannMatrix(i,pointCount+2) = sum(interfaceNormals(i,1:3)*planeBasisVectors(1,:))
         neumannMatrix(i,pointCount+3) = sum(interfaceNormals(i,1:3)*planeBasisVectors(2,:))
       else
@@ -352,7 +344,7 @@
         = neumannMatrix(i,pointCount+1:pointCount+3)
     end do
 
-    rhs(pointCount+1) = 1.0_MPAS_r8
+    rhs(pointCount+1) = 1.0
     rhs(pointCount+2) = sum(destinationPoint(1:3)*planeBasisVectors(1,:))
     rhs(pointCount+3) = sum(destinationPoint(1:3)*planeBasisVectors(2,:))
 
@@ -375,24 +367,23 @@
 
   subroutine ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &amp;
     sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
-    alpha, dirichletCoefficients, neumannCoefficients, dminfo)
+    alpha, dirichletCoefficients, neumannCoefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isInterface
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount), intent(out) :: &amp;
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount), intent(out) :: &amp;
       dirichletCoefficients, neumannCoefficients
-    type (dm_info), intent(in) :: dminfo
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
-    real(MPAS_r8), dimension(:), pointer :: rhs, rhsCopy, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+    real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
 
@@ -403,11 +394,11 @@
     allocate(coeffs(matrixSize))  
     allocate(pivotIndices(matrixSize))  
 
-    dirichletMatrix = 0.0_MPAS_r8
-    neumannMatrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    rhsCopy = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    dirichletMatrix = 0.0
+    neumannMatrix = 0.0
+    rhs = 0.0
+    rhsCopy = 0.0
+    coeffs = 0.0
 
     call setUpScalarRBFMatrixAndRHS(pointCount, &amp;
       sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
@@ -415,10 +406,10 @@
       neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
 
     do i = 1, pointCount
-      dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+      dirichletMatrix(i,pointCount+1) = 1.0
       dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
       if(isInterface(i)) then
-        neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+        neumannMatrix(i,pointCount+1) = 0.0
         neumannMatrix(i,pointCount+2:pointCount+4) = interfaceNormals(i,1:3)
       else
         neumannMatrix(i,pointCount+1:pointCount+4) &amp;
@@ -430,7 +421,7 @@
         = neumannMatrix(i,pointCount+1:pointCount+4)
     end do
 
-    rhs(pointCount+1) = 1.0_MPAS_r8
+    rhs(pointCount+1) = 1.0
     rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
 
     ! solve each linear system
@@ -452,22 +443,21 @@
 
   subroutine ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &amp;
     sourcePoints, unitVectors, destinationPoint, &amp;
-    alpha, coefficients, dminfo)
+    alpha, coefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
-    type (dm_info), intent(in) :: dminfo
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8), dimension(:,:), pointer :: matrix
-    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: matrix
+    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+3 ! extra space for constant vector 
 
@@ -476,9 +466,9 @@
     allocate(coeffs(matrixSize,3))  
     allocate(pivotIndices(matrixSize))  
 
-    matrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    matrix = 0.0
+    rhs = 0.0
+    coeffs = 0.0
 
     call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 3, &amp;
       sourcePoints, unitVectors, destinationPoint, &amp;
@@ -490,7 +480,7 @@
         = matrix(i,pointCount+1:pointCount+3)
     end do
     do i = 1, 3
-      rhs(pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+      rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
     end do
 
     ! solve each linear system
@@ -508,29 +498,28 @@
 
   subroutine ibInterp_computeVectorNoSlipPlanarRBFCoefficients(pointCount, &amp;
     sourcePoints, unitVectors, destinationPoint, &amp;
-    alpha, planeBasisVectors, coefficients, dminfo)
+    alpha, planeBasisVectors, coefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(2,3) :: planeBasisVectors
-    real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
-    type (dm_info), intent(in) :: dminfo
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+    real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
 
-    real(MPAS_r8), dimension(pointCount,2) :: planarSourcePoints
-    real(MPAS_r8), dimension(pointCount,2) :: planarUnitVectors
-    real(MPAS_r8), dimension(2) :: planarDestinationPoint
+    real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
+    real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
+    real(kind=RKIND), dimension(2) :: planarDestinationPoint
 
-    real(MPAS_r8), dimension(:,:), pointer :: matrix
-    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: matrix
+    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+2 ! space for constant vector in plane
 
@@ -539,9 +528,9 @@
     allocate(coeffs(matrixSize,2))  
     allocate(pivotIndices(matrixSize)) 
 
-    matrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    matrix = 0.0
+    rhs = 0.0
+    coeffs = 0.0
 
     do i = 1, pointCount
       planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:)) 
@@ -583,24 +572,23 @@
 
   subroutine ibInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &amp;
     sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
-    alpha, coefficients, dminfo)
+    alpha, coefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isTangentToInterface
-    integer(MPAS_i4), dimension(pointCount), intent(in) :: normalVectorIndex
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
-    type (dm_info), intent(in) :: dminfo
+    integer, dimension(pointCount), intent(in) :: normalVectorIndex
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8), dimension(:,:), pointer :: matrix
-    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: matrix
+    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+3 ! extra space for constant vector 
 
@@ -609,9 +597,9 @@
     allocate(coeffs(matrixSize,3))  
     allocate(pivotIndices(matrixSize))  
 
-    matrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    matrix = 0.0
+    rhs = 0.0
+    coeffs = 0.0
 
     call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &amp;
       sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
@@ -625,7 +613,7 @@
         = matrix(i,pointCount+1:pointCount+3)
     end do
     do i = 1, 3
-      rhs(pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+      rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
     end do
 
     ! solve each linear system
@@ -643,29 +631,28 @@
 
   subroutine ibInterp_computeVectorFreeSlipPlanarRBFCoefficients(pointCount, &amp;
     sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
-    alpha, planeBasisVectors, coefficients, dminfo)
+    alpha, planeBasisVectors, coefficients)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isTangentToInterface
-    integer(MPAS_i4), dimension(pointCount), intent(in) :: normalVectorIndex
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: unitVectors
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(2,3) :: planeBasisVectors
-    real(MPAS_r8), dimension(pointCount, 3), intent(out) :: coefficients
-    type (dm_info), intent(in) :: dminfo
+    integer, dimension(pointCount), intent(in) :: normalVectorIndex
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
 
-    integer(MPAS_i4) :: i, j
-    integer(MPAS_i4) :: matrixSize
+    integer :: i, j
+    integer :: matrixSize
 
-    real(MPAS_r8), dimension(pointCount,2) :: planarSourcePoints
-    real(MPAS_r8), dimension(pointCount,2) :: planarUnitVectors
-    real(MPAS_r8), dimension(2) :: planarDestinationPoint
+    real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
+    real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
+    real(kind=RKIND), dimension(2) :: planarDestinationPoint
 
-    real(MPAS_r8), dimension(:,:), pointer :: matrix
-    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
-    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+    real(kind=RKIND), dimension(:,:), pointer :: matrix
+    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+    integer, dimension(:), pointer :: pivotIndices
 
     matrixSize = pointCount+2 ! space for constant vector in plane
 
@@ -674,9 +661,9 @@
     allocate(coeffs(matrixSize,2))  
     allocate(pivotIndices(matrixSize)) 
 
-    matrix = 0.0_MPAS_r8
-    rhs = 0.0_MPAS_r8
-    coeffs = 0.0_MPAS_r8
+    matrix = 0.0
+    rhs = 0.0
+    coeffs = 0.0
 
     do i = 1, pointCount
       planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:)) 
@@ -723,17 +710,16 @@
 ! private subroutines
 !!!!!!!!!!!!!!!!!!!!!
 
-  subroutine exit_MPAS(string, dminfo)
+  subroutine exit_MPAS(string)
     character(len=*), intent(in) :: string
-    type (dm_info), intent(in) :: dminfo
 
     write(0,*) string
-    call dmpar_abort(dminfo)
+    stop
   end subroutine exit_MPAS
 
   function evaluateRBF(rSquared) result(rbfValue)
-    real(MPAS_r8), intent(in) :: rSquared
-    real(MPAS_r8) :: rbfValue
+    real(kind=RKIND), intent(in) :: rSquared
+    real(kind=RKIND) :: rbfValue
 
     ! inverse multiquadratic
     rbfValue = 1/sqrt(1 + rSquared)
@@ -741,8 +727,8 @@
   end function evaluateRBF
 
   subroutine evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
-    real(MPAS_r8), intent(in) :: rSquared
-    real(MPAS_r8), intent(out) :: rbfValue, rbfDerivOverR
+    real(kind=RKIND), intent(in) :: rSquared
+    real(kind=RKIND), intent(out) :: rbfValue, rbfDerivOverR
 
     ! inverse multiquadratic
     rbfValue = 1/sqrt(1 + rSquared)
@@ -751,8 +737,8 @@
   end subroutine evaluateRBFAndDeriv
 
   subroutine evaluateRBFAndDerivs(rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv)
-    real(MPAS_r8), intent(in) :: rSquared
-    real(MPAS_r8), intent(out) :: rbfValue, rbfDerivOverR, rbfSecondDeriv
+    real(kind=RKIND), intent(in) :: rSquared
+    real(kind=RKIND), intent(out) :: rbfValue, rbfDerivOverR, rbfSecondDeriv
 
     ! inverse multiquadratic
     rbfValue = 1/sqrt(1 + rSquared)
@@ -765,19 +751,19 @@
     sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
     alpha, dirichletMatrix, neumannMatrix, rhs)
 
-    integer(MPAS_i4), intent(in) :: pointCount
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isInterface
-    real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
-    real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: &amp;
+    real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: &amp;
       dirichletMatrix, neumannMatrix
-    real(MPAS_r8), dimension(pointCount), intent(out) :: rhs
+    real(kind=RKIND), dimension(pointCount), intent(out) :: rhs
 
-    integer(MPAS_i4) :: i, j
+    integer :: i, j
 
-    real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+    real(kind=RKIND) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
 
     do j = 1, pointCount
       if(isInterface(j)) then
@@ -811,17 +797,17 @@
     sourcePoints, unitVectors, destinationPoint, &amp;
     alpha, matrix, rhs)
 
-    integer(MPAS_i4), intent(in) :: pointCount, dimensions
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
-    real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: matrix
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(out) :: rhs
+    integer, intent(in) :: pointCount, dimensions
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: unitVectors
+    real(kind=RKIND), dimension(dimensions), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: matrix
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(out) :: rhs
 
-    integer(MPAS_i4) :: i, j
+    integer :: i, j
 
-    real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+    real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
 
     do j = 1, pointCount
       do i = j, pointCount
@@ -844,19 +830,19 @@
     sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
     alpha, matrix, rhs)
 
-    integer(MPAS_i4), intent(in) :: pointCount, dimensions
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+    integer, intent(in) :: pointCount, dimensions
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: sourcePoints
     logical, dimension(pointCount), intent(in) :: isTangentToInterface
-    integer(MPAS_i4), dimension(pointCount), intent(in) :: normalVectorIndex
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
-    real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
-    real(MPAS_r8), intent(in) :: alpha
-    real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: matrix
-    real(MPAS_r8), dimension(pointCount,dimensions), intent(out) :: rhs
+    integer, dimension(pointCount), intent(in) :: normalVectorIndex
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: unitVectors
+    real(kind=RKIND), dimension(dimensions), intent(in) :: destinationPoint
+    real(kind=RKIND), intent(in) :: alpha
+    real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: matrix
+    real(kind=RKIND), dimension(pointCount,dimensions), intent(out) :: rhs
 
-    integer(MPAS_i4) :: i, j
+    integer :: i, j
 
-    real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalVector(dimensions), &amp;
+    real(kind=RKIND) :: rSquared, rbfValue, rbfDerivOverR, normalVector(dimensions), &amp;
       normalDotX, unitVectorDotProduct
 
     do j = 1, pointCount
@@ -951,12 +937,12 @@
 ! Copyright (c) Tao Pang 2001.
 !
   IMPLICIT NONE
-  integer(MPAS_i4), INTENT (IN) :: N
-  integer(MPAS_i4) :: I,J
-  integer(MPAS_i4), INTENT (OUT), DIMENSION (N) :: INDX
-  real(MPAS_r8), INTENT (INOUT), DIMENSION (N,N) :: A
-  real(MPAS_r8), INTENT (INOUT), DIMENSION (N) :: B
-  real(MPAS_r8), INTENT (OUT), DIMENSION (N) :: X
+  integer, INTENT (IN) :: N
+  integer :: I,J
+  integer, INTENT (OUT), DIMENSION (N) :: INDX
+  real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+  real(kind=RKIND), INTENT (INOUT), DIMENSION (N) :: B
+  real(kind=RKIND), INTENT (OUT), DIMENSION (N) :: X
 !
   CALL ELGS (A,N,INDX)
 !
@@ -1002,12 +988,12 @@
 ! in X(N,N) in the output.  Copyright (c) Tao Pang 2001.
 !
   IMPLICIT NONE
-  integer(MPAS_i4), INTENT (IN) :: N
-  integer(MPAS_i4) :: I,J,K
-  integer(MPAS_i4), INTENT (OUT), DIMENSION (N) :: INDX
-  real(MPAS_r8), INTENT (INOUT), DIMENSION (N,N):: A
-  real(MPAS_r8), INTENT (OUT), DIMENSION (N,N):: X
-  real(MPAS_r8), DIMENSION (N,N) :: B
+  integer, INTENT (IN) :: N
+  integer :: I,J,K
+  integer, INTENT (OUT), DIMENSION (N) :: INDX
+  real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
+  real(kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
+  real(kind=RKIND), DIMENSION (N,N) :: B
 !
   DO I = 1, N
     DO J = 1, N
@@ -1049,12 +1035,12 @@
 ! INDX(N) records the pivoting order.  Copyright (c) Tao Pang 2001.
 !
   IMPLICIT NONE
-  integer(MPAS_i4), INTENT (IN) :: N
-  integer(MPAS_i4) :: I,J,K,ITMP
-  integer(MPAS_i4), INTENT (OUT), DIMENSION (N) :: INDX
-  real(MPAS_r8) :: C1,PI,PI1,PJ
-  real(MPAS_r8), INTENT (INOUT), DIMENSION (N,N) :: A
-  real(MPAS_r8), DIMENSION (N) :: C
+  integer, INTENT (IN) :: N
+  integer :: I,J,K,ITMP
+  integer, INTENT (OUT), DIMENSION (N) :: INDX
+  real(kind=RKIND) :: C1,PI,PI1,PJ
+  real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+  real(kind=RKIND), DIMENSION (N) :: C
 !
 ! Initialize the index
 !

Added: branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F                                (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F        2010-04-26 22:58:02 UTC (rev 213)
@@ -0,0 +1,657 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+module interp_tests
+   use dmpar
+   use grid_types
+   use IB_interpolation
+
+   implicit none
+   private
+   save
+
+
+   public :: doTests
+
+ contains
+
+  subroutine doTests(grid)
+
+    implicit none
+
+    type (grid_meta), intent(inout) :: grid 
+
+print *, &quot;initialization test&quot;
+    call test_ibInterp_initialize(grid) 
+
+print *, &quot;2D RBF test&quot;
+    call test_ibInterp_get2DRBF
+
+print *, &quot;ScalarRBFAndConstant test&quot;
+    call test_ibInterp_computeScalarRBFAndConstant
+
+    print *, &quot;tests finished, stopping...&quot;
+    stop
+
+  end subroutine doTests
+
+  subroutine test_ibInterp_initialize(grid)
+
+    implicit none
+
+    type (grid_meta), intent(inout) :: grid 
+
+    integer :: iCell, iEdge
+
+    real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+    real(kind=RKIND), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &amp;
+      edgeLocations
+    real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
+
+    radialUnitVectors =&gt; grid % radialUnitVectors % array
+    edgeNormalVectors =&gt; grid % edgeNormalVectors % array
+    edgeLocations =&gt; grid % edgeLocations % array
+    cellTangentPlane =&gt; grid % cellTangentPlane % array
+
+    xCell       =&gt; grid % xCell % array
+    yCell       =&gt; grid % yCell % array
+    zCell       =&gt; grid % zCell % array
+
+    iCell = 1
+    iEdge = grid%edgesOnCell%array(1,iCell) 
+
+    print *, &quot;cell: &quot;, iCell, &quot;radial unit vector: &quot;, radialUnitVectors(:,iCell)
+    print *, &quot;x, y, z: &quot;, xCell(iCell), yCell(iCell), zCell(iCell)
+    print *, &quot;tangentPlane: &quot;, cellTangentPlane(:,:,iCell)
+
+    print *, &quot;edge: &quot;, iEdge, &quot;normal: &quot;, edgeNormalVectors(:,iEdge)
+    print *, &quot;location: &quot;, edgeLocations(:,iEdge)
+
+  end subroutine test_ibInterp_initialize
+
+  subroutine test_ibInterp_get2DRBF
+ 
+    integer, parameter :: pointCount = 16, coeffCount = 19, dimensions = 2, fieldCount = 1
+    integer, parameter :: polynomialFlag = kMPAS_IBRBFPolynomialLinear
+    real(kind=RKIND), parameter :: alpha = 1.35, epsilon = 1e-6
+    real(kind=RKIND), dimension(pointCount,dimensions) :: points
+    real(kind=RKIND), dimension(pointCount) :: fieldValues
+    real(kind=RKIND), dimension(coeffCount) :: coefficients
+    real(kind=RKIND), dimension(9,2) :: evaluationPoint
+    real(kind=RKIND), dimension(9,6) :: derivs
+    real(kind=RKIND) :: x0, y0
+
+    integer :: xIndex, yIndex, pointIndex
+
+    pointIndex = 1
+    do yIndex=-1,2
+      do xIndex = -1,2
+        points(pointIndex, 1) = xIndex
+        points(pointIndex, 2) = yIndex
+        fieldValues(pointIndex) = xIndex**2*yIndex
+        pointIndex = pointIndex + 1
+      end do
+    end do
+
+    print *, &quot;points:&quot;, points
+    print *, &quot;fieldValues:&quot;, fieldValues
+    call ibInterp_get2DRBFCoefficients(pointCount, coeffCount, dimensions, points, &amp;
+      fieldValues, polynomialFlag, alpha, coefficients)
+
+    print *, &quot;coefficients:&quot;, coefficients
+
+    x0 = 0.3
+    y0 = 0.6
+
+    pointIndex = 1
+    do yIndex=-1,1
+      do xIndex = -1,1
+        evaluationPoint(pointIndex,1) = epsilon*xIndex + x0
+        evaluationPoint(pointIndex,2) = epsilon*yIndex + y0
+        call ibInterp_evaluate2DRBFWithDerivs(fieldCount, coeffCount, pointCount, &amp;
+          coefficients, evaluationPoint(pointIndex,:), points, polynomialFlag, alpha, &amp;
+          derivs(pointIndex,:))
+        pointIndex = pointIndex + 1
+      end do
+    end do
+
+    print *, &quot;evaluation points:&quot;, evaluationPoint
+
+    print *, &quot;exact derivs:&quot;, x0**2*y0, 2*x0*y0, x0**2, 2*y0, 2*x0, 0.0
+    print *, &quot;rbf derivs:&quot;, derivs(5,:)
+    print *, &quot;finite diff derivs:&quot;, derivs(5,1), (derivs(6,1)-derivs(4,1))/(2*epsilon), &amp;
+      (derivs(8,1)-derivs(2,1))/(2*epsilon), (derivs(6,1)-2*derivs(5,1)+derivs(4,1))/epsilon**2, &amp;
+      (derivs(9,1)-derivs(7,1)-derivs(3,1)+derivs(1,1))/(4*epsilon**2), &amp;
+      (derivs(8,1)-2*derivs(5,1)+derivs(2,1))/epsilon**2
+
+  end subroutine test_ibInterp_get2DRBF
+
+  subroutine test_ibInterp_computeScalarRBFAndConstant
+
+    integer, parameter :: pointCount = 32
+    real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
+    logical, dimension(pointCount) :: isInterface
+    real(kind=RKIND), dimension(pointCount,3) :: interfaceNormals
+    real(kind=RKIND), dimension(3) :: destinationPoint
+    real(kind=RKIND), dimension(pointCount) :: &amp;
+      dirichletCoefficients, neumannCoefficients
+
+    real(kind=RKIND), dimension(pointCount) :: testFunction
+    real(kind=RKIND), dimension(2) :: functionAtDestination
+
+    real(kind=RKIND) :: epsilon, alpha
+
+    epsilon = 1e-6
+    alpha = 1.35
+
+    ! set up a fake set of fluid points, interface points and interfce normals
+    sourcePoints(1,:) = (/0, 0, -1/)
+    sourcePoints(2,:) = (/1, 0, -1/)
+    sourcePoints(3,:) = (/0, 1, -1/)
+    sourcePoints(4,:) = (/1, 1, -1/)
+    sourcePoints(5,:) = (/0, -1, 0/)
+    sourcePoints(6,:) = (/1, -1, 0/)
+    sourcePoints(7,:) = (/-1, 0, 0/)
+    sourcePoints(8,:) = (/0, 0, 0/)
+    sourcePoints(9,:) = (/1, 0, 0/)
+    sourcePoints(10,:) = (/2, 0, 0/)
+    sourcePoints(11,:) = (/-1, 1, 0/)
+    sourcePoints(12,:) = (/0, 1, 0/)
+    sourcePoints(13,:) = (/1, 1, 0/)
+    sourcePoints(14,:) = (/2, 1, 0/)
+    sourcePoints(15,:) = (/0, 2, 0/)
+    sourcePoints(16,:) = (/1, 2, 0/)
+    sourcePoints(17,:) = (/0, -1, 1/)
+    sourcePoints(18,:) = (/1, -1, 1/)
+    sourcePoints(19,:) = (/-1, 0, 1/)
+    sourcePoints(20,:) = (/0, 0, 1/)
+    sourcePoints(21,:) = (/1, 0, 1/)
+    sourcePoints(22,:) = (/2, 0, 1/)
+    sourcePoints(23,:) = (/-1, 1, 1/)
+    sourcePoints(24,:) = (/0, 1, 1/)
+    sourcePoints(25,:) = (/1, 1, 1/)
+    sourcePoints(26,:) = (/0, 0, 2/)
+    sourcePoints(27,:) = (/1.92857142857143, 0.857142857142857, 0.785714285714286/)
+    sourcePoints(28,:) = (/-0.0714285714285715, 1.85714285714286, 0.785714285714286/)
+    sourcePoints(29,:) = (/0.857142857142857, 1.71428571428571, 0.571428571428571/)
+    sourcePoints(30,:) = (/0.928571428571429, -0.142857142857143, 1.78571428571429/)
+    sourcePoints(31,:) = (/-0.142857142857143, 0.714285714285714, 1.57142857142857/)
+    sourcePoints(32,:) = (/0.785714285714286, 0.571428571428571, 1.35714285714286/)
+
+    interfaceNormals(:,1) = 0.267261241912424
+    interfaceNormals(:,2) = 0.534522483824849
+    interfaceNormals(:,3) = 0.801783725737273
+
+    isInterface(1:26) = .false.
+    isInterface(27:32) = .true.
+
+    testFunction(1:26) = sourcePoints(1:26,1)**2*sourcePoints(1:26,2)
+    testFunction(27:32) = interfaceNormals(27:32,1)*2.0*sourcePoints(27:32,1) &amp;
+      *sourcePoints(27:32,2) + interfaceNormals(27:32,2)*sourcePoints(27:32,1)**2
+print *, &quot;testFunction:&quot;, testFunction
+
+    destinationPoint = sourcePoints(13,:)
+    call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &amp;
+      sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+      alpha, dirichletCoefficients, neumannCoefficients)
+
+    functionAtDestination(1) = sum(testFunction*neumannCoefficients)
+
+print *, &quot;function value:&quot;, functionAtDestination(1)
+print *, &quot;function should be:&quot;, testFunction(13)
+
+    destinationPoint = sourcePoints(32,:) + epsilon*interfaceNormals(1,:)
+
+    call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &amp;
+      sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+      alpha, dirichletCoefficients, neumannCoefficients)
+
+!print *, &quot;dirichlet:&quot;, dirichletCoefficients
+!print *, &quot;neumann:&quot;, neumannCoefficients
+
+    ! test to make sure it reconstructs a function properly
+
+    functionAtDestination(1) = sum(testFunction*neumannCoefficients)
+
+    destinationPoint = sourcePoints(32,:) - epsilon*interfaceNormals(1,:)
+
+    call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &amp;
+      sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+      alpha, dirichletCoefficients, neumannCoefficients)
+
+    functionAtDestination(2) = sum(testFunction*neumannCoefficients)
+print *, &quot;function value:&quot;, functionAtDestination
+print *, &quot;function normal deriv:&quot;, (functionAtDestination(1)-functionAtDestination(2))/(2.0*epsilon)
+print *, &quot;deriv should be:&quot;, testFunction(32)
+
+  end subroutine test_ibInterp_computeScalarRBFAndConstant
+
+!  subroutine test_ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &amp;
+!    sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+!    alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+!
+!    integer, intent(in) :: pointCount
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+!    logical, dimension(pointCount), intent(in) :: isInterface
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+!    real(kind=RKIND), intent(in) :: alpha
+!    real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+!    real(kind=RKIND), dimension(pointCount), intent(out) :: &amp;
+!      dirichletCoefficients, neumannCoefficients
+!    type (dm_info), intent(in) :: dminfo
+!
+!    integer :: i, j
+!    integer :: matrixSize
+!
+!    real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+!    real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+!    integer, dimension(:), pointer :: pivotIndices
+!
+!    matrixSize = pointCount+3 !! 3 extra space for constant and 2 planar dimensions
+!
+!    allocate(dirichletMatrix(matrixSize,matrixSize))  
+!    allocate(neumannMatrix(matrixSize,matrixSize))  
+!    allocate(rhs(matrixSize))  
+!    allocate(rhsCopy(matrixSize))  
+!    allocate(coeffs(matrixSize))  
+!    allocate(pivotIndices(matrixSize))  
+!
+!    dirichletMatrix = 0.0
+!    neumannMatrix = 0.0
+!    rhs = 0.0
+!    rhsCopy = 0.0
+!    coeffs = 0.0
+!
+!    call setUpScalarRBFMatrixAndRHS(pointCount, &amp;
+!      sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+!      alpha, dirichletMatrix(1:pointCount,1:pointCount), &amp;
+!      neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+!
+!    do i = 1, pointCount
+!      dirichletMatrix(i,pointCount+1) = 1.0
+!      dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVectors(1,:))
+!      dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVectors(2,:))
+!      if(isInterface(i)) then
+!        neumannMatrix(i,pointCount+1) = 0.0
+!        neumannMatrix(i,pointCount+2) = sum(interfaceNormals(i,1:3)*planeBasisVectors(1,:))
+!        neumannMatrix(i,pointCount+3) = sum(interfaceNormals(i,1:3)*planeBasisVectors(2,:))
+!      else
+!        neumannMatrix(i,pointCount+1:pointCount+3) &amp;
+!          = dirichletMatrix(i,pointCount+1:pointCount+3)
+!      end if
+!      dirichletMatrix(pointCount+1:pointCount+3,i) &amp;
+!        = dirichletMatrix(i,pointCount+1:pointCount+3)
+!      neumannMatrix(pointCount+1:pointCount+3,i) &amp;
+!        = neumannMatrix(i,pointCount+1:pointCount+3)
+!    end do
+!
+!    rhs(pointCount+1) = 1.0
+!    rhs(pointCount+2) = sum(destinationPoint(1:3)*planeBasisVectors(1,:))
+!    rhs(pointCount+3) = sum(destinationPoint(1:3)*planeBasisVectors(2,:))
+!
+!    ! solve each linear system
+!    rhsCopy = rhs
+!    call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+!    dirichletCoefficients = coeffs(1:pointCount)
+!
+!    call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+!    neumannCoefficients = coeffs(1:pointCount)
+!
+!    deallocate(dirichletMatrix)  
+!    deallocate(neumannMatrix)  
+!    deallocate(rhs)  
+!    deallocate(rhsCopy)  
+!    deallocate(coeffs)  
+!    deallocate(pivotIndices) 
+!
+!  end subroutine test_ibInterp_computeScalarRBFAndPlaneCoefficients
+!
+!  subroutine test_ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &amp;
+!    sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+!    alpha, dirichletCoefficients, neumannCoefficients, dminfo)
+!
+!    integer, intent(in) :: pointCount
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+!    logical, dimension(pointCount), intent(in) :: isInterface
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+!    real(kind=RKIND), intent(in) :: alpha
+!    real(kind=RKIND), dimension(pointCount), intent(out) :: &amp;
+!      dirichletCoefficients, neumannCoefficients
+!    type (dm_info), intent(in) :: dminfo
+!
+!    integer :: i, j
+!    integer :: matrixSize
+!
+!    real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+!    real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+!    integer, dimension(:), pointer :: pivotIndices
+!
+!    matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
+!
+!    allocate(dirichletMatrix(matrixSize,matrixSize))  
+!    allocate(neumannMatrix(matrixSize,matrixSize))  
+!    allocate(rhs(matrixSize))  
+!    allocate(rhsCopy(matrixSize))  
+!    allocate(coeffs(matrixSize))  
+!    allocate(pivotIndices(matrixSize))  
+!
+!    dirichletMatrix = 0.0
+!    neumannMatrix = 0.0
+!    rhs = 0.0
+!    rhsCopy = 0.0
+!    coeffs = 0.0
+!
+!    call setUpScalarRBFMatrixAndRHS(pointCount, &amp;
+!      sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+!      alpha, dirichletMatrix(1:pointCount,1:pointCount), &amp;
+!      neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+!
+!    do i = 1, pointCount
+!      dirichletMatrix(i,pointCount+1) = 1.0
+!      dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
+!      if(isInterface(i)) then
+!        neumannMatrix(i,pointCount+1) = 0.0
+!        neumannMatrix(i,pointCount+2:pointCount+4) = interfaceNormals(i,1:3)
+!      else
+!        neumannMatrix(i,pointCount+1:pointCount+4) &amp;
+!          = dirichletMatrix(i,pointCount+1:pointCount+4)
+!      end if
+!      dirichletMatrix(pointCount+1:pointCount+4,i) &amp;
+!        = dirichletMatrix(i,pointCount+1:pointCount+4)
+!      neumannMatrix(pointCount+1:pointCount+4,i) &amp;
+!        = neumannMatrix(i,pointCount+1:pointCount+4)
+!    end do
+!
+!    rhs(pointCount+1) = 1.0
+!    rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+!
+!    ! solve each linear system
+!    rhsCopy = rhs
+!    call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+!    dirichletCoefficients = coeffs(1:pointCount)
+!
+!    call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+!    neumannCoefficients = coeffs(1:pointCount)
+!
+!    deallocate(dirichletMatrix)  
+!    deallocate(neumannMatrix)  
+!    deallocate(rhs)  
+!    deallocate(rhsCopy)  
+!    deallocate(coeffs)  
+!    deallocate(pivotIndices) 
+!
+!  end subroutine test_ibInterp_computeScalarRBFAndLinearCoefficients
+!
+!  subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &amp;
+!    sourcePoints, unitVectors, destinationPoint, &amp;
+!    alpha, coefficients, dminfo)
+!
+!    integer, intent(in) :: pointCount
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+!    real(kind=RKIND), intent(in) :: alpha
+!    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+!    type (dm_info), intent(in) :: dminfo
+!
+!    integer :: i, j
+!    integer :: matrixSize
+!
+!    real(kind=RKIND), dimension(:,:), pointer :: matrix
+!    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+!    integer, dimension(:), pointer :: pivotIndices
+!
+!    matrixSize = pointCount+3 ! extra space for constant vector 
+!
+!    allocate(matrix(matrixSize,matrixSize))  
+!    allocate(rhs(matrixSize,3))  
+!    allocate(coeffs(matrixSize,3))  
+!    allocate(pivotIndices(matrixSize))  
+!
+!    matrix = 0.0
+!    rhs = 0.0
+!    coeffs = 0.0
+!
+!    call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 3, &amp;
+!      sourcePoints, unitVectors, destinationPoint, &amp;
+!      alpha, matrix, rhs)
+!
+!    do i = 1, pointCount
+!      matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
+!      matrix(pointCount+1:pointCount+3,i) &amp;
+!        = matrix(i,pointCount+1:pointCount+3)
+!    end do
+!    do i = 1, 3
+!      rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
+!    end do
+!
+!    ! solve each linear system
+!    do i = 1, 3
+!      call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+!    end do
+!    coefficients = coeffs(1:pointCount,:)
+!
+!    deallocate(matrix)  
+!    deallocate(rhs)  
+!    deallocate(coeffs)  
+!    deallocate(pivotIndices) 
+!
+!  end subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients 
+!
+!  subroutine test_ibInterp_computeVectorNoSlipPlanarRBFCoefficients(pointCount, &amp;
+!    sourcePoints, unitVectors, destinationPoint, &amp;
+!    alpha, planeBasisVectors, coefficients, dminfo)
+!
+!    integer, intent(in) :: pointCount
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+!    real(kind=RKIND), intent(in) :: alpha
+!    real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+!    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+!    type (dm_info), intent(in) :: dminfo
+!
+!    integer :: i, j
+!    integer :: matrixSize
+!
+!    real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
+!
+!    real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
+!    real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
+!    real(kind=RKIND), dimension(2) :: planarDestinationPoint
+!
+!    real(kind=RKIND), dimension(:,:), pointer :: matrix
+!    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+!    integer, dimension(:), pointer :: pivotIndices
+!
+!    matrixSize = pointCount+2 ! space for constant vector in plane
+!
+!    allocate(matrix(matrixSize,matrixSize))  
+!    allocate(rhs(matrixSize,2))  
+!    allocate(coeffs(matrixSize,2))  
+!    allocate(pivotIndices(matrixSize)) 
+!
+!    matrix = 0.0
+!    rhs = 0.0
+!    coeffs = 0.0
+!
+!    do i = 1, pointCount
+!      planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:)) 
+!      planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:)) 
+!      planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:)) 
+!      planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:)) 
+!    end do
+!    planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:)) 
+!    planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:)) 
+!    call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 2, &amp;
+!      planarSourcePoints, planarUnitVectors, planarDestinationPoint, &amp;
+!      alpha, matrix, rhs)
+!
+!    do i = 1, pointCount
+!      matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+!      matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+!      matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
+!    end do
+!    rhs(pointCount+1,:) = planeBasisVectors(1,:)
+!    rhs(pointCount+2,:) = planeBasisVectors(2,:)
+!
+!    ! solve each linear system
+!    call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+!    call LEGS(matrix, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+!
+!    coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &amp;
+!      + planeBasisVectors(2,1)*coeffs(1:pointCount,2) 
+!    coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &amp;
+!      + planeBasisVectors(2,2)*coeffs(1:pointCount,2) 
+!    coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &amp;
+!      + planeBasisVectors(2,3)*coeffs(1:pointCount,2) 
+!
+!    deallocate(matrix)  
+!    deallocate(rhs)  
+!    deallocate(coeffs)  
+!    deallocate(pivotIndices)
+!
+!  end subroutine test_ibInterp_computeVectorNoSlipPlanarRBFCoefficients 
+!
+!  subroutine test_ibInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &amp;
+!    sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
+!    alpha, coefficients, dminfo)
+!
+!    integer, intent(in) :: pointCount
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+!    logical, dimension(pointCount), intent(in) :: isTangentToInterface
+!    integer, dimension(pointCount), intent(in) :: normalVectorIndex
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+!    real(kind=RKIND), intent(in) :: alpha
+!    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+!    type (dm_info), intent(in) :: dminfo
+!
+!    integer :: i, j
+!    integer :: matrixSize
+!
+!    real(kind=RKIND), dimension(:,:), pointer :: matrix
+!    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+!    integer, dimension(:), pointer :: pivotIndices
+!
+!    matrixSize = pointCount+3 ! extra space for constant vector 
+!
+!    allocate(matrix(matrixSize,matrixSize))  
+!    allocate(rhs(matrixSize,3))  
+!    allocate(coeffs(matrixSize,3))  
+!    allocate(pivotIndices(matrixSize))  
+!
+!    matrix = 0.0
+!    rhs = 0.0
+!    coeffs = 0.0
+!
+!    call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &amp;
+!      sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
+!      alpha, matrix, rhs)
+!
+!    do i = 1, pointCount
+!      if(.not. isTangentToInterface(i)) then
+!        matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
+!      end if
+!      matrix(pointCount+1:pointCount+3,i) &amp;
+!        = matrix(i,pointCount+1:pointCount+3)
+!    end do
+!    do i = 1, 3
+!      rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
+!    end do
+!
+!    ! solve each linear system
+!    do i = 1, 3
+!      call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+!    end do
+!    coefficients = coeffs(1:pointCount,:)
+!
+!    deallocate(matrix)  
+!    deallocate(rhs)  
+!    deallocate(coeffs)  
+!    deallocate(pivotIndices)
+!
+!   end subroutine test_ibInterp_computeVectorFreeSlipRBFCoefficients 
+!
+!  subroutine test_ibInterp_computeVectorFreeSlipPlanarRBFCoefficients(pointCount, &amp;
+!    sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &amp;
+!    alpha, planeBasisVectors, coefficients, dminfo)
+!
+!    integer, intent(in) :: pointCount
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+!    logical, dimension(pointCount), intent(in) :: isTangentToInterface
+!    integer, dimension(pointCount), intent(in) :: normalVectorIndex
+!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+!    real(kind=RKIND), intent(in) :: alpha
+!    real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+!    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+!    type (dm_info), intent(in) :: dminfo
+!
+!    integer :: i, j
+!    integer :: matrixSize
+!
+!    real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
+!    real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
+!    real(kind=RKIND), dimension(2) :: planarDestinationPoint
+!
+!    real(kind=RKIND), dimension(:,:), pointer :: matrix
+!    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+!    integer, dimension(:), pointer :: pivotIndices
+!
+!    matrixSize = pointCount+2 ! space for constant vector in plane
+!
+!    allocate(matrix(matrixSize,matrixSize))  
+!    allocate(rhs(matrixSize,2))  
+!    allocate(coeffs(matrixSize,2))  
+!    allocate(pivotIndices(matrixSize)) 
+!
+!    matrix = 0.0
+!    rhs = 0.0
+!    coeffs = 0.0
+!
+!    do i = 1, pointCount
+!      planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:)) 
+!      planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:)) 
+!      planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:)) 
+!      planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:)) 
+!    end do
+!    planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:)) 
+!    planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:)) 
+!    call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &amp;
+!      planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &amp;
+!      planarDestinationPoint, alpha, matrix, rhs)
+!
+!    do i = 1, pointCount
+!      if(.not. isTangentToInterface(i)) then
+!        matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+!        matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+!      end if
+!      matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
+!    end do
+!    rhs(pointCount+1,:) = planeBasisVectors(1,:)
+!    rhs(pointCount+2,:) = planeBasisVectors(2,:)
+!
+!    ! solve each linear system
+!    call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+!    call LEGS(matrix, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+!
+!    coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &amp;
+!      + planeBasisVectors(2,1)*coeffs(1:pointCount,2) 
+!    coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &amp;
+!      + planeBasisVectors(2,2)*coeffs(1:pointCount,2) 
+!    coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &amp;
+!      + planeBasisVectors(2,3)*coeffs(1:pointCount,2) 
+!    coefficients = coeffs(1:pointCount,:)
+!
+!    deallocate(matrix)  
+!    deallocate(rhs)  
+!    deallocate(coeffs)  
+!    deallocate(pivotIndices)
+!
+!   end subroutine test_ibInterp_computeVectorFreeSlipPlanarRBFCoefficients 
+!
+
+end module interp_tests
+

</font>
</pre>