<p><b>xylar@lanl.gov</b> 2010-04-21 16:12:14 -0600 (Wed, 21 Apr 2010)</p><p>Checking in the new interpolation module I'm working on just for safe keeping...<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F                                (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F        2010-04-21 22:12:14 UTC (rev 197)
@@ -0,0 +1,952 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+module IB_interpoation
+   use dmpar
+
+   implicit none
+   private
+   save
+
+   integer, parameter, public :: &amp;
+     MPAS_i4 = selected_int_kind(6), &amp;
+     MPAS_r8 = selected_real_kind(13)
+
+   public :: ibInterp_initialize, ibInterp_get2DRBFCoefficients, &amp;
+     ibInterp_evaluate2DRBFWithDerivs, &amp;
+     ibInterp_computeScalarRBFCoefficients, &amp;
+     ibInterp_computeVectorNoSlipCoefficientsRBF, &amp;
+     ibInterp_computeVectorFreeSlipCoefficientsRBF
+
+   integer(MPAS_i4), public, parameter :: &amp;
+     kMPAS_IBRBFPolynomialConstant = 0, &amp;
+     kMPAS_IBRBFPolynomialLinear = 1
+
+
+ contains
+
+   subroutine ibInterp_initialize(grid)
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+   ! Purpose: 
+   !
+   ! Input: 
+   !
+   ! Output:
+   !       
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+
+      implicit none
+
+      type (grid_meta), intent(inout) :: grid 
+
+      integer(MPAS_i4) :: nCells, nEdges
+      integer(MPAS_i4), dimension(:,:), pointer :: cellsOnEdge
+      integer(MPAS_i4) :: iEdge, iCell, cell1, cell2
+      real(MPAS_r8), dimension(:), pointer :: xCell, yCell, zCell
+      real(MPAS_r8), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &amp;
+        edgeLocations
+      real(MPAS_r8), dimension(3) :: xHatPlane, yHatPlane
+      real(MPAS_r8) :: normalDotRHat
+
+      xCell       =&gt; grid % xCell % array
+      yCell       =&gt; grid % yCell % array
+      zCell       =&gt; grid % zCell % array
+      cellsOnEdge =&gt; grid % cellsOnEdge % array
+      nCells      = grid % nCells
+      nEdges      = grid % nEdges
+
+      radialUnitVectors =&gt; grid % radialUnitVectors % array
+      edgeNormalVectors =&gt; grid % edgeNormalVectors % array
+      edgeLocations =&gt; grid % edgeLocations % array
+      cellTangentPlane =&gt; grid % cellTangentPlane % array
+
+      ! init arrays
+      edgeNormalVectors = 0
+      radialUnitVectors = 0
+
+      ! loop over all cells to be solved on this block
+      do iCell=1,nCells
+         radialUnitVectors(iCell,1) = xCell(iCell)
+         radialUnitVectors(iCell,2) = yCell(iCell)
+         radialUnitVectors(iCell,3) = zCell(iCell)
+         call unit_vector_in_R3(radialUnitVectors(iCell,:))
+      end do
+
+      do iEdge = 1,nEdges
+        cell1 = cellsOnEdge(1,iEdge)
+        cell2 = cellsOnEdge(2,iEdge)
+        edgeNormalVectors(iEdge,1) = xCell(cell2) - xCell(cell1)
+        edgeNormalVectors(iEdge,2) = yCell(cell2) - yCell(cell1)
+        edgeNormalVectors(iEdge,3) = zCell(cell2) - zCell(cell1)
+        edgeLocations(iEdge,1) = 0.5*(xCell(cell2) + xCell(cell1))
+        edgeLocations(iEdge,2) = 0.5*(yCell(cell2) + yCell(cell1))
+        edgeLocations(iEdge,3) = 0.5*(zCell(cell2) + zCell(cell1))
+        call unit_vector_in_R3(edgeNormalVectors(iEdge,:))
+      end do
+
+      do iCell=1,nCells
+        iEdge = edgesOnCell(1,iCell)
+        ! xHat and yHat are a local basis in the plane of the horizontal cell
+        ! we arbitrarily choose xHat to point toward the first edge
+        normalDotRHat = sum(edgeNormalVectors(iEdge,:)*radialUnitVectors(iCell,:))
+        xHatPlane = edgeNormalVectors(iEdge,:) - normalDotRHat*radialUnitVectors(iCell,1)
+         call unit_vector_in_R3(xHatPlane)
+
+         call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
+         call unit_vector_in_R3(yHatPlane) ! just to be sure...
+         cellTangentPlane(iCell, 1, :) = xHatPlane
+         cellTangentPlane(iCell, 2, :) = yHatPlane
+      end do
+
+   end subroutine ibInterp_initialize
+
+  subroutine ibInterp_getRBFCoefficients(pointCount, coeffCount, dimensions, points, fieldValues, &amp;
+    polynomialFlag, alpha, coefficients, dminfo)
+ 
+    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(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
+
+    coefficients = 0.0_MPAS_r8
+    matrix = 0.0_MPAS_r8
+    rhs = 0.0_MPAS_r8
+
+    rhs(1:pointCount) = fieldValues
+
+    do j=1,pointCount
+      do i=j,pointCount
+        rSquared = sum((points(j,:) - points(i,:))**2)/alpha**2
+        matrix(i,j) = evaluateRBF(rSquared)
+        matrix(j,i) = matrix(i,j)
+      end do
+    end do
+    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
+      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+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)
+    end if
+
+    call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &amp;
+      coefficients(1:matrixSize), pivotIndices(1:matrixSize))
+
+  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
+
+    real(MPAS_r8), dimension(6,fieldCount), intent(out) :: derivs
+
+    integer(MPAS_i4) :: pointIndex
+    real(MPAS_r8) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+    derivs = 0.0_MPAS_r8
+    do pointIndex = 1, pointCount
+      x = (evaluationPoint(1) - points(pointIndex,1))
+      y = (evaluationPoint(2) - points(pointIndex,2))
+      rSquared = x**2 + y**2
+      call evaluateRBFAndDerivs(rSquared/alpha**2,  rbfValue, rbfDerivOverR, rbfSecondDeriv)
+      rbfDerivOverR = rbfDerivOverR/alpha**2
+      rbfSecondDeriv = rbfSecondDeriv/alpha**2
+      if(rSquared/alpha**2 &lt; 1e-7) then
+        ! terms 2,3 and 5 are zero by radial symmetry of the RBF
+        derivs(1,:) = derivs(1,:) + coefficients(pointIndex,:)*rbfValue
+        derivs(4,:) = derivs(4,:) + coefficients(pointIndex,:)*rbfSecondDeriv
+        derivs(6,:) = derivs(6,:) + coefficients(pointIndex,:)*rbfSecondDeriv
+      else
+        call evaluateRBFAndDerivs(rSquared,  rbfValue, rbfDerivOverR, rbfSecondDeriv)
+        derivs(1,:) = derivs(1,:) + coefficients(pointIndex,:)*rbfValue
+        derivs(2,:) = derivs(2,:) + coefficients(pointIndex,:)*rbfDerivOverR*x
+        derivs(3,:) = derivs(3,:) + coefficients(pointIndex,:)*rbfDerivOverR*y
+        derivs(4,:) = derivs(4,:) + coefficients(pointIndex,:) &amp;
+          * (rbfSecondDeriv*x**2 + rbfDerivOverR*y**2)/rSquared
+        derivs(5,:) = derivs(5,:) + coefficients(pointIndex,:) &amp;
+          * (rbfSecondDeriv - rbfDerivOverR)*x*y/rSquared
+        derivs(6,:) = derivs(6,:) + coefficients(pointIndex,:) &amp;
+          * (rbfSecondDeriv*y**2 + rbfDerivOverR*x**2)/rSquared
+      end if
+    end do
+    if(polynomialFlag == kMPAS_IBRBFPolynomialConstant) then
+      derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:)
+    else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
+      derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:) &amp;
+        + coefficients(pointCount+2,:)*evaluationPoint(1) &amp;
+        + coefficients(pointCount+3,:)*evaluationPoint(2)
+      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)
+    end if
+  end subroutine ibInterp_evaluate2DRBFWithDerivs
+
+  subroutine ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &amp;
+    sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+    polynomialFlag, alpha, dirichletCoefficients, neumannCoefficients, dminfo)
+
+    integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+    real(MPAS_r8), 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;
+      dirichletCoefficients, neumannCoefficients
+    type (dm_info), intent(in) :: dminfo
+
+    integer(MPAS_i4) :: i, j
+    integer(MPAS_i4), parameter :: matrixSize = pointCount+1 !! 1 extra space for constant 
+
+
+    real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
+    real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
+    integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+    rhsCopy = 0.0_MPAS_r8
+    coeffs = 0.0_MPAS_r8
+
+    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_MPAS_r8
+      if(isInterface(i)) then
+        neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+      else
+        neumannMatrix(i,pointCount+1) = dirichletMatrix(i,pointCount+1)
+      end if
+      dirichletMatrix(pointCount+1,i) = dirichletMatrix(i,pointCount+1)
+      neumannMatrix(pointCount+1,i) = neumannMatrix(i,pointCount+1)
+    end do
+
+    rhs(pointCount+1) = 1.0_MPAS_r8
+
+    ! 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)
+
+  end subroutine ibInterp_computeScalarRBFAndConstantCoefficients
+
+  subroutine ibInterp_computeScalarRBFAndPlane(pointCount, &amp;
+    sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+    alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+
+    integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+    real(MPAS_r8), 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;
+      dirichletCoefficients, neumannCoefficients
+    type (dm_info), intent(in) :: dminfo
+
+    integer(MPAS_i4) :: i, j
+    integer(MPAS_i4), parameter :: matrixSize = pointCount+3 ! extra space for constant and plane
+
+    real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
+    real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
+    integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+    rhsCopy = 0.0_MPAS_r8
+    coeffs = 0.0_MPAS_r8
+
+    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_MPAS_r8
+      dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVector(1,:))
+      dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVector(2,:))
+      if(isInterface(i)) then
+        neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+        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_MPAS_r8
+    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)
+
+  end subroutine ibInterp_computeScalarRBFCoefficients
+
+  subroutine ibInterp_computeScalarRBFAndLinear(pointCount, &amp;
+    sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+    alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+
+    integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+    real(MPAS_r8), 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(3,2) :: planeBasisVectors
+    real(MPAS_r8), dimension(pointCount), intent(out) :: &amp;
+      dirichletCoefficients, neumannCoefficients
+    type (dm_info), intent(in) :: dminfo
+
+    integer(MPAS_i4) :: i, j
+    integer(MPAS_i4), parameter :: matrixSize = pointCount+4 ! extra space for constant and plane
+
+
+    real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
+    real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
+    integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+    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
+
+    do j = 1, pointCount
+      if(isInterface(j)) cycle
+      do i = 1, pointCount
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        rbfValue = evaluateRBF(rSquared)
+        dirichletMatrix(i,j) = rbfValue
+        neumannMatrix(i,j) = rbfValue
+      end do
+    end do
+
+    do j = 1, pointCount
+      if(.not. isInterface(j)) cycle
+      do i = 1, pointCount 
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        normalDotX = sum(interfaceNormals(j,:) &amp;
+          * (sourcePoints(j,:)-sourcePoints(i,:)))
+        call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+        rbfDerivOverR = rbfDerivOverR/alpha**2
+        dirichletMatrix(i,j) = rbfValue
+        neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+      end do
+    end do
+
+    do j = 1, pointCount
+      rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+      rhs(j) = evaluateRBF(rSquared)
+    end do
+
+      do i = 1, pointCount
+        dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+        if(isInterface(i)) then
+          neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+        else
+          neumannMatrix(i,pointCount+1) = dirichletMatrix(i,pointCount+1)
+        end if
+        dirichletMatrix(pointCount+1,i) = dirichletMatrix(i,pointCount+1)
+        neumannMatrix(pointCount+1,i) = neumannMatrix(i,pointCount+1)
+      end do
+
+      rhs(pointCount+1) = 1.0_MPAS_r8
+    else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
+      do i = 1, pointCount
+        dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+        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+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_MPAS_r8
+      rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+    else
+      call exit_MPAS(&quot;ibInterp_computeScalarRBFCoefficients: not a valid polynomialFlag&quot;, dminfo)
+    end if
+
+
+    ! 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)
+
+  end subroutine ibInterp_computeScalarRBFCoefficients
+
+  subroutine ibInterp_computeVectorNoSlipCoefficientsRBF(pointCount, &amp;
+    dimensions, sourcePoints, unitVectors, destinationPoint, &amp;
+    alpha, coefficients, dminfo)
+
+    integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
+    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, dimensions), intent(out) :: coefficients
+    type (dm_info), intent(in) :: dminfo
+
+    integer(MPAS_i4) :: matrixSize, i, j
+
+    real(MPAS_r8), dimension(:,:), pointer :: matrix
+    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+    matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+
+    allocate(matrix(matrixSize,matrixSize))
+    allocate(rhs(matrixSize,dimensions))
+    allocate(coeffs(matrixSize,dimensions))
+    allocate(pivotIndices(matrixSize))
+
+    matrix = 0.0_MPAS_r8
+    rhs = 0.0_MPAS_r8
+    coeffs = 0.0_MPAS_r8
+
+    do j = 1, pointCount
+      do i = 1, pointCount
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+        rbfValue = evaluateRBF(rSquared)
+        matrix(j,i) = rbfValue*unitVectorDotProduct
+      end do
+      rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+      rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+    end do
+
+    do i = 1, pointCount
+      matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
+      matrix(pointCount+1:pointCount+dimensions,i) &amp;
+        = matrix(i,pointCount+1:pointCount+dimensions)
+    end do
+    do i = 1, dimensions
+      rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+    end do
+
+    ! solve each linear system
+    do i = 1, dimensions
+      call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+    end do
+    coefficients = coeffs(1:pointCount,:)
+
+  end subroutine ibInterp_computeVectorNoSlipCoefficientsRBF 
+
+  subroutine ibInterp_computeVectorNoSlipCoefficientsRBF(pointCount, &amp;
+    dimensions, sourcePoints, unitVectors, destinationPoint, &amp;
+    alpha, coefficients, dminfo)
+
+    integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
+    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, dimensions), intent(out) :: coefficients
+    type (dm_info), intent(in) :: dminfo
+
+    integer(MPAS_i4) :: matrixSize, i, j
+
+    real(MPAS_r8), dimension(:,:), pointer :: matrix
+    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+    matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+
+    allocate(matrix(matrixSize,matrixSize))
+    allocate(rhs(matrixSize,dimensions))
+    allocate(coeffs(matrixSize,dimensions))
+    allocate(pivotIndices(matrixSize))
+
+    matrix = 0.0_MPAS_r8
+    rhs = 0.0_MPAS_r8
+    coeffs = 0.0_MPAS_r8
+
+    do j = 1, pointCount
+      do i = 1, pointCount
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+        rbfValue = evaluateRBF(rSquared)
+        matrix(j,i) = rbfValue*unitVectorDotProduct
+      end do
+      rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+      rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+    end do
+
+    do i = 1, pointCount
+      matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
+      matrix(pointCount+1:pointCount+dimensions,i) &amp;
+        = matrix(i,pointCount+1:pointCount+dimensions)
+    end do
+    do i = 1, dimensions
+      rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+    end do
+
+    ! solve each linear system
+    do i = 1, dimensions
+      call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+    end do
+    coefficients = coeffs(1:pointCount,:)
+
+  end subroutine ibInterp_computeVectorNoSlipCoefficientsRBF 
+
+  !!!!!! Not implemented yet, identical to No Slip !!!!!!!!!
+  subroutine ibInterp_computeVectorFreeSlipCoefficientsRBF(pointCount, &amp;
+    dimensions, sourcePoints, unitVectors, destinationPoint, &amp;
+    alpha, coefficients, dminfo)
+
+    integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
+    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, dimensions), intent(out) :: coefficients
+    type (dm_info), intent(in) :: dminfo
+
+    integer(MPAS_i4) :: matrixSize, i, j
+
+    real(MPAS_r8), dimension(:,:), pointer :: matrix
+    real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+    integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+    matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+
+    allocate(matrix(matrixSize,matrixSize))
+    allocate(rhs(matrixSize,dimensions))
+    allocate(coeffs(matrixSize,dimensions))
+    allocate(pivotIndices(matrixSize))
+
+    matrix = 0.0_MPAS_r8
+    rhs = 0.0_MPAS_r8
+    coeffs = 0.0_MPAS_r8
+
+    do j = 1, pointCount
+      do i = 1, pointCount
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+        rbfValue = evaluateRBF(rSquared)
+        matrix(j,i) = rbfValue*unitVectorDotProduct
+      end do
+      rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+      rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+    end do
+
+    do i = 1, pointCount
+      matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
+      matrix(pointCount+1:pointCount+dimensions,i) &amp;
+        = matrix(i,pointCount+1:pointCount+dimensions)
+    end do
+    do i = 1, dimensions
+      rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+    end do
+
+    ! solve each linear system
+    do i = 1, dimensions
+      call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+    end do
+    coefficients = coeffs(1:pointCount,:)
+
+  end subroutine ibInterp_computeVectorFreeSlipCoefficientsRBF 
+
+
+
+!!!!!!!!!!!!!!!!!!!!!
+! private subroutines
+!!!!!!!!!!!!!!!!!!!!!
+
+  subroutine exit_MPAS(string, dminfo)
+    char(len=*), intent(in) :: string
+    type (dm_info), intent(in) :: dminfo
+
+    write(0,*) string
+    call dmpar_abort(dminfo)
+  end subroutine exit_MPAS
+
+  function evaluateRBF(rSquared) result(rbfValue)
+    real(MPAS_r8), intent(in) :: rSquared
+    real(MPAS_r8) :: rbfValue
+
+    ! inverse multiquadratic
+    rbfValue = 1/sqrt(1 + rSquared)
+
+  end function evaluateRBF
+
+  subroutine evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+    real(MPAS_r8), intent(in) :: rSquared
+    real(MPAS_r8), intent(out) :: rbfValue, rbfDerivOverR
+
+    ! inverse multiquadratic
+    rbfValue = 1/sqrt(1 + rSquared)
+    rbfDerivOverR = -rbfValue**3
+
+  end subroutine evaluateRBFAndDeriv
+
+  subroutine evaluateRBFAndDerivs(rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv)
+    real(MPAS_r8), intent(in) :: rSquared
+    real(MPAS_r8), intent(out) :: rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+    ! inverse multiquadratic
+    rbfValue = 1/sqrt(1 + rSquared)
+    rbfDerivOverR = -rbfValue**3
+    rbfSecondDeriv = (2*rSquared-1)*rbfValue**5
+
+  end subroutine evaluateRBFAndDerivs
+
+  subroutine setUpScalarRBFMatrixAndRHS(pointCount, &amp;
+    sourcePoints, isInterface, interfaceNormals, destinationPoint, &amp;
+    alpha, dirichletMatrix, neumannMatrix, rhs)
+
+    integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+    real(MPAS_r8), 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;
+      dirichletCoefficients, neumannCoefficients
+    real(MPAS_r8), dimension(pointCount), intent(out) :: rhs
+
+    integer(MPAS_i4) :: i, j
+
+
+    integer(MPAS_i4), dimension(matrixSize) :: pivotIndices
+
+    real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+    dirichletMatrix = 0.0_MPAS_r8
+    neumannMatrix = 0.0_MPAS_r8
+    rhs = 0.0_MPAS_r8
+
+    do j = 1, pointCount
+      if(isInterface(j)) cycle
+      do i = 1, pointCount
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        rbfValue = evaluateRBF(rSquared)
+        dirichletMatrix(i,j) = rbfValue
+        neumannMatrix(i,j) = rbfValue
+      end do
+    end do
+
+    do j = 1, pointCount
+      if(.not. isInterface(j)) cycle
+      do i = 1, pointCount 
+        rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+        normalDotX = sum(interfaceNormals(j,:) &amp;
+          * (sourcePoints(j,:)-sourcePoints(i,:)))
+        call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+        rbfDerivOverR = rbfDerivOverR/alpha**2
+        dirichletMatrix(i,j) = rbfValue
+        neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+      end do
+    end do
+
+    do j = 1, pointCount
+      rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+      rhs(j) = evaluateRBF(rSquared)
+    end do
+
+  end subroutine setUpScalarRBFMatrixAndRHS
+
+   subroutine unit_vector_in_R3(xin)
+     implicit none
+     real (kind=RKIND), intent(inout) :: xin(3)
+     real (kind=RKIND) :: mag
+     mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
+     xin(:) = xin(:) / mag
+   end subroutine unit_vector_in_R3
+
+!======================================================================
+! BEGINNING OF CROSS_PRODUCT_IN_R3
+!======================================================================
+        subroutine cross_product_in_R3(p_1,p_2,p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE: compute p_1 cross p_2 and place in p_out
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+        real (kind=RKIND), intent(in) ::                            &amp;
+                        p_1 (3),                                      &amp;
+                        p_2 (3)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+        real (kind=RKIND), intent(out) ::                           &amp;
+                        p_out (3)
+
+        p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+        p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+        p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+
+        end subroutine cross_product_in_R3
+!======================================================================
+! END OF CROSS_PRODUCT_IN_R3
+!======================================================================
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 4.3   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!                                                                       !
+! Please Note:                                                          !
+!                                                                       !
+! (1) This computer program is written by Tao Pang in conjunction with  !
+!     his book, &quot;An Introduction to Computational Physics,&quot; published   !
+!     by Cambridge University Press in 1997.                            !
+!                                                                       !
+! (2) No warranties, express or implied, are made for this program.     !
+!                                                                       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!PROGRAM EX43
+!
+!
+! An example of solving linear equation set A(N,N)*X(N) = B(N)
+! with the partial-pivoting Gaussian elimination scheme.  The
+! numerical values are for the Wheatstone bridge example discussed
+! in Section 4.1 in the book with all resistances being 100 ohms
+! and the voltage 200 volts.
+!
+!  IMPLICIT NONE
+!  INTEGER, PARAMETER :: N=3
+!  INTEGER :: I,J
+!  INTEGER, DIMENSION (N) :: INDX
+!  REAL, DIMENSION (N) :: X,B
+!  REAL, DIMENSION (N,N) :: A
+!  DATA B /200.0,0.0,0.0/, &amp;
+!       ((A(I,J), J=1,N),I=1,N) /100.0,100.0,100.0,-100.0, &amp;
+!                         300.0,-100.0,-100.0,-100.0, 300.0/
+!
+!  CALL LEGS (A,N,B,X,INDX)
+!
+!  WRITE (6, &quot;(F16.8)&quot;) (X(I), I=1,N)
+!END PROGRAM EX43
+
+
+SUBROUTINE LEGS (A,N,B,X,INDX)
+!
+! Subroutine to solve the equation A(N,N)*X(N) = B(N) with the
+! partial-pivoting Gaussian elimination scheme.
+! 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
+!
+  CALL ELGS (A,N,INDX)
+!
+  DO I = 1, N-1
+    DO J = I+1, N
+      B(INDX(J)) = B(INDX(J))-A(INDX(J),I)*B(INDX(I))
+    END DO
+  END DO
+!
+  X(N) = B(INDX(N))/A(INDX(N),N)
+  DO I = N-1, 1, -1
+    X(I) = B(INDX(I))
+    DO J = I+1, N
+      X(I) = X(I)-A(INDX(I),J)*X(J)
+    END DO
+    X(I) =  X(I)/A(INDX(I),I)
+  END DO
+!
+END SUBROUTINE LEGS
+!
+
+
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 4.4   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!                                                                       !
+! Please Note:                                                          !
+!                                                                       !
+! (1) This computer program is written by Tao Pang in conjunction with  !
+!     his book, &quot;An Introduction to Computational Physics,&quot; published   !
+!     by Cambridge University Press in 1997.                            !
+!                                                                       !
+! (2) No warranties, express or implied, are made for this program.     !
+!                                                                       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+SUBROUTINE MIGS (A,N,X,INDX)
+!
+! Subroutine to invert matrix A(N,N) with the inverse stored
+! 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
+!
+  DO I = 1, N
+    DO J = 1, N
+      B(I,J) = 0.0
+    END DO
+  END DO
+  DO I = 1, N
+    B(I,I) = 1.0
+  END DO
+!
+  CALL ELGS (A,N,INDX)
+!
+  DO I = 1, N-1
+    DO J = I+1, N
+      DO K = 1, N
+        B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
+      END DO
+    END DO
+  END DO
+!
+  DO I = 1, N
+    X(N,I) = B(INDX(N),I)/A(INDX(N),N)
+    DO J = N-1, 1, -1
+      X(J,I) = B(INDX(J),I)
+      DO K = J+1, N
+        X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
+      END DO
+      X(J,I) =  X(J,I)/A(INDX(J),J)
+    END DO
+  END DO
+END SUBROUTINE MIGS
+
+
+SUBROUTINE ELGS (A,N,INDX)
+!
+! Subroutine to perform the partial-pivoting Gaussian elimination.
+! A(N,N) is the original matrix in the input and transformed matrix
+! plus the pivoting element ratios below the diagonal in the output.
+! 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
+!
+! Initialize the index
+!
+  DO I = 1, N
+    INDX(I) = I
+  END DO
+!
+! Find the rescaling factors, one from each row
+!
+  DO I = 1, N
+    C1= 0.0
+    DO J = 1, N
+      !C1 = AMAX1(C1,ABS(A(I,J)))
+      C1 = MAX(C1,ABS(A(I,J)))
+    END DO
+    C(I) = C1
+  END DO
+!
+! Search the pivoting (largest) element from each column
+!
+  DO J = 1, N-1
+    PI1 = 0.0
+    DO I = J, N
+      PI = ABS(A(INDX(I),J))/C(INDX(I))
+      IF (PI.GT.PI1) THEN
+        PI1 = PI
+        K   = I
+      ENDIF
+    END DO
+!
+! Interchange the rows via INDX(N) to record pivoting order
+!
+    ITMP    = INDX(J)
+    INDX(J) = INDX(K)
+    INDX(K) = ITMP
+    DO I = J+1, N
+      PJ  = A(INDX(I),J)/A(INDX(J),J)
+!
+! Record pivoting ratios below the diagonal
+!
+      A(INDX(I),J) = PJ
+!
+! Modify other elements accordingly
+!
+      DO K = J+1, N
+        A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
+      END DO
+    END DO
+  END DO
+!
+END SUBROUTINE ELGS
+
+end module IB_interpoation
+

</font>
</pre>