<p><b>cnewman@lanl.gov</b> 2010-05-24 08:49:07 -0600 (Mon, 24 May 2010)</p><p>added files and consolidated common code for ocean and sw in operator_addition branch<br>
</p><hr noshade><pre><font color="gray">Modified: branches/operator_addition/src/core_ocean/module_test_cases.F
===================================================================
--- branches/operator_addition/src/core_ocean/module_test_cases.F        2010-05-21 17:53:43 UTC (rev 300)
+++ branches/operator_addition/src/core_ocean/module_test_cases.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -3,8 +3,8 @@
    use grid_types
    use configure
    use constants
+   use geometry
 
-
    contains
 
 
@@ -448,25 +448,6 @@
    end subroutine sw_test_case_6
 
 
-   real function sphere_distance(lat1, lon1, lat2, lon2, radius)
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
-   !   sphere with given radius.
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
-      real (kind=RKIND) :: arg1
-
-      arg1 = sqrt( sin(0.5*(lat2-lat1))**2 +  &amp;
-                   cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
-      sphere_distance = 2.*radius*asin(arg1)
-
-   end function sphere_distance
-
-
    real function AA(theta)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! A, used in height field computation for Rossby-Haurwitz wave

Modified: branches/operator_addition/src/core_ocean/module_time_integration.F
===================================================================
--- branches/operator_addition/src/core_ocean/module_time_integration.F        2010-05-21 17:53:43 UTC (rev 300)
+++ branches/operator_addition/src/core_ocean/module_time_integration.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -7,6 +7,7 @@
    ! xsad 10-02-05:
    use vector_reconstruction
    ! xsad 10-02-05 end
+   use differential_operators
 
 
    contains
@@ -603,50 +604,18 @@
       ! Compute circulation and relative vorticity at each vertex
       !
       circulation(:,:) = 0.0
-      do iEdge=1,nEdges
-         if (verticesOnEdge(1,iEdge) &lt;= nVertices) then
-            do k=1,nVertLevels
-               circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
-            end do
-         end if
-         if (verticesOnEdge(2,iEdge) &lt;= nVertices) then
-            do k=1,nVertLevels
-               circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
-            end do
-         end if
-      end do
-      do iVertex=1,nVertices
-         do k=1,nVertLevels
-            vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
-         end do
-      end do
 
+      call compute_circulation(circulation, u, verticesOnEdge, dcEdge, nEdges, nVertices, nVertLevels)
+  
+      call compute_vorticity(vorticity, circulation, areaTriangle, nVertices, nVertLevels)
 
       !
       ! Compute the divergence at each cell center
       !
       divergence(:,:) = 0.0
-      do iEdge=1,nEdges
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-         if (cell1 &lt;= nCells) then
-            do k=1,nVertLevels
-              divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
-            enddo
-         endif
-         if(cell2 &lt;= nCells) then
-            do k=1,nVertLevels
-              divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
-            enddo
-         end if
-      end do
-      do iCell = 1,nCells
-        r = 1.0 / areaCell(iCell)
-        do k = 1,nVertLevels
-           divergence(k,iCell) = divergence(k,iCell) * r
-        enddo
-      enddo
 
+      call compute_divergence(divergence, u, dvEdge, areaCell, cellsOnEdge, nEdges, nCells, nVertLevels)
+
       !
       ! Compute kinetic energy in each cell
       !
@@ -718,17 +687,9 @@
       ! Compute gradient of PV in the tangent direction
       !   ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
       !
-      do iEdge = 1,nEdges
-         do k = 1,nVertLevels
-           gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &amp;
-                              dvEdge(iEdge)
-         enddo
-      enddo
 
-      !
-      ! Compute pv at the edges
-      !   ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
-      !
+      call compute_gradient(gradPVt, pv_vertex, verticesOnEdge, dvEdge, nEdges, nVertLevels)
+
       pv_edge(:,:) = 0.0
       do iVertex = 1,nVertices
         do i=1,grid % vertexDegree
@@ -771,17 +732,10 @@
       ! Compute gradient of PV in normal direction
       !   ( this computes gradPVn for all edges bounding real cells )
       !
+
       gradPVn(:,:) = 0.0
-      do iEdge = 1,nEdges
-        if( cellsOnEdge(1,iEdge) &lt;= nCells .and. cellsOnEdge(2,iEdge) &lt;= nCells) then
-          do k = 1,nVertLevels
-            gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &amp;
-                                 dcEdge(iEdge)
-          enddo
-        endif
-      enddo
+      call compute_gradient_normal(gradPVn,pv_cell,cellsOnEdge,dcEdge,nEdges,nVertLevels,nCells)
 
-
       ! Modify PV edge with upstream bias.
       !
       do iEdge = 1,nEdges

Modified: branches/operator_addition/src/core_sw/module_test_cases.F
===================================================================
--- branches/operator_addition/src/core_sw/module_test_cases.F        2010-05-21 17:53:43 UTC (rev 300)
+++ branches/operator_addition/src/core_sw/module_test_cases.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -3,6 +3,7 @@
    use grid_types
    use configure
    use constants
+   use geometry
 
 
    contains
@@ -475,25 +476,6 @@
    end subroutine sw_test_case_6
 
 
-   real function sphere_distance(lat1, lon1, lat2, lon2, radius)
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
-   !   sphere with given radius.
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
-      real (kind=RKIND) :: arg1
-
-      arg1 = sqrt( sin(0.5*(lat2-lat1))**2 +  &amp;
-                   cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
-      sphere_distance = 2.*radius*asin(arg1)
-
-   end function sphere_distance
-
-
    real function AA(theta)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! A, used in height field computation for Rossby-Haurwitz wave

Modified: branches/operator_addition/src/core_sw/module_time_integration.F
===================================================================
--- branches/operator_addition/src/core_sw/module_time_integration.F        2010-05-21 17:53:43 UTC (rev 300)
+++ branches/operator_addition/src/core_sw/module_time_integration.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -5,8 +5,8 @@
    use configure
    use constants
    use dmpar
+   use differential_operators
 
-
    contains
 
 
@@ -521,53 +521,18 @@
       ! Compute circulation and relative vorticity at each vertex
       !
       circulation(:,:) = 0.0
-      do iEdge=1,nEdges
-         if (verticesOnEdge(1,iEdge) &lt;= nVertices) then
-            do k=1,nVertLevels
-               circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
-            end do
-         end if
-         if (verticesOnEdge(2,iEdge) &lt;= nVertices) then
-            do k=1,nVertLevels
-               circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
-            end do
-         end if
-      end do
-      do iVertex=1,nVertices
-         do k=1,nVertLevels
-            vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
-         end do
-      end do
 
-
+      call compute_circulation(circulation, u, verticesOnEdge, dcEdge, nEdges, nVertices, nVertLevels)
+  
+      call compute_vorticity(vorticity, circulation, areaTriangle, nVertices, nVertLevels)
+ 
       !
       ! Compute the divergence at each cell center
       !
       divergence(:,:) = 0.0
-      do iEdge=1,nEdges
-         cell1 = cellsOnEdge(1,iEdge)
-         cell2 = cellsOnEdge(2,iEdge)
-         if (cell1 &lt;= nCells) then
-            do k=1,nVertLevels
-              divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
-            enddo
-         endif
-         if(cell2 &lt;= nCells) then
-            do k=1,nVertLevels
-              divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
-            enddo
-         end if
-      end do
-      do iCell = 1,nCells
-        r = 1.0 / areaCell(iCell)
-        do k = 1,nVertLevels
-           divergence(k,iCell) = divergence(k,iCell) * r
-        enddo
-      enddo
 
-      !
-      ! Compute kinetic energy in each cell
-      !
+      call compute_divergence(divergence, u, dvEdge, areaCell, cellsOnEdge, nEdges, nCells, nVertLevels)
+  
       ke(:,:) = 0.0
       do iCell=1,nCells
          do i=1,nEdgesOnCell(iCell)
@@ -636,13 +601,9 @@
       ! Compute gradient of PV in the tangent direction
       !   ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
       !
-      do iEdge = 1,nEdges
-         do k = 1,nVertLevels
-           gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &amp;
-                              dvEdge(iEdge)
-         enddo
-      enddo
 
+      call compute_gradient(gradPVt, pv_vertex, verticesOnEdge, dvEdge, nEdges, nVertLevels)
+
       !
       ! Compute pv at the edges
       !   ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
@@ -691,15 +652,9 @@
       !   ( this computes gradPVn for all edges bounding real cells )
       !
       gradPVn(:,:) = 0.0
-      do iEdge = 1,nEdges
-        if( cellsOnEdge(1,iEdge) &lt;= nCells .and. cellsOnEdge(2,iEdge) &lt;= nCells) then
-          do k = 1,nVertLevels
-            gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &amp;
-                                 dcEdge(iEdge)
-          enddo
-        endif
-      enddo
 
+      call compute_gradient_normal(gradPVn,pv_cell,cellsOnEdge,dcEdge,nEdges,nVertLevels,nCells)
+
       ! Modify PV edge with upstream bias.
       !
       do iEdge = 1,nEdges
@@ -778,5 +733,4 @@
 
    end subroutine enforce_boundaryEdge
 
-
 end module time_integration

Modified: branches/operator_addition/src/operators/Makefile
===================================================================
--- branches/operator_addition/src/operators/Makefile        2010-05-21 17:53:43 UTC (rev 300)
+++ branches/operator_addition/src/operators/Makefile        2010-05-24 14:49:07 UTC (rev 301)
@@ -1,6 +1,6 @@
 .SUFFIXES: .F .o
 
-OBJS = module_vector_reconstruction.o
+OBJS = module_linear_algebra.o module_geometry.o module_vector_reconstruction.o module_differential_operators.o 
 
 all: operators
 
@@ -9,6 +9,12 @@
 
 module_vector_reconstruction.o:
 
+module_geometry.o:
+
+module_differential_operators.o:
+
+module_linear_algebra.o:
+
 clean:
         $(RM) *.o *.mod *.f90 libops.a
 

Added: branches/operator_addition/src/operators/module_differential_operators.F
===================================================================
--- branches/operator_addition/src/operators/module_differential_operators.F                                (rev 0)
+++ branches/operator_addition/src/operators/module_differential_operators.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -0,0 +1,113 @@
+module differential_operators
+
+   contains
+
+   subroutine compute_gradient(grad, u_vertex, verticesOnEdge, dvEdge, nEdges, nVertLevels)
+      implicit none
+      real (kind=RKIND), dimension(:,:), intent(inout) :: grad
+      real (kind=RKIND), dimension(:,:), intent(in) :: u_vertex
+      real (kind=RKIND), dimension(:), intent(in) :: dvEdge
+      integer, dimension(:,:), intent(in) :: verticesOnEdge
+      integer, intent(in) :: nEdges, nVertLevels
+      integer :: iEdge, k
+
+      do iEdge = 1,nEdges
+         do k = 1,nVertLevels
+            grad(k,iEdge) = (u_vertex(k,verticesOnEdge(2,iEdge)) - u_vertex(k,verticesOnEdge(1,iEdge))) / &amp;
+            dvEdge(iEdge)
+         enddo
+      enddo
+   end subroutine compute_gradient
+
+   subroutine compute_gradient_normal(grad, u_cell, cellsOnEdge, dcEdge, nEdges, nVertLevels, nCells)
+      implicit none
+      real (kind=RKIND), dimension(:,:), intent(inout) :: grad
+      real (kind=RKIND), dimension(:,:), intent(in) :: u_cell
+      real (kind=RKIND), dimension(:), intent(in) :: dcEdge
+      integer, dimension(:,:), intent(in) :: cellsOnEdge
+      integer, intent(in) :: nEdges, nVertLevels
+      integer, intent(in) :: nCells
+      integer :: iEdge, k
+
+      do iEdge = 1,nEdges
+         if( cellsOnEdge(1,iEdge) &lt;= nCells .and. cellsOnEdge(2,iEdge) &lt;= nCells) then
+            do k = 1,nVertLevels
+               grad(k,iEdge) = (u_cell(k,cellsOnEdge(2,iEdge)) - u_cell(k,cellsOnEdge(1,iEdge))) / &amp;
+               dcEdge(iEdge)
+            enddo
+         endif
+      enddo
+   end subroutine compute_gradient_normal
+
+   subroutine compute_divergence(divergence, u, dvEdge, areaCell, cellsOnEdge, nEdges, nCells, nVertLevels)
+      implicit none
+      real (kind=RKIND), dimension(:,:), intent(inout) :: divergence
+      real (kind=RKIND), dimension(:,:), intent(in) :: u
+      real (kind=RKIND), dimension(:), intent(in) ::  dvEdge, areaCell
+      real (kind=RKIND) :: r
+      integer, dimension(:,:), intent(in):: cellsOnEdge
+      integer, intent(in) :: nEdges, nCells, nVertLevels
+      integer :: iEdge, iCell, k, cell1, cell2   
+
+      do iEdge=1,nEdges
+         cell1 = cellsOnEdge(1,iEdge)
+         cell2 = cellsOnEdge(2,iEdge)
+         if (cell1 &lt;= nCells) then
+            do k=1,nVertLevels
+              divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+            enddo
+         endif
+         if(cell2 &lt;= nCells) then
+            do k=1,nVertLevels
+              divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+            enddo
+         end if
+      end do
+      do iCell = 1,nCells
+        r = 1.0 / areaCell(iCell)
+        do k = 1,nVertLevels
+           divergence(k,iCell) = divergence(k,iCell) * r
+        enddo
+      enddo
+   end subroutine compute_divergence
+
+   subroutine compute_circulation(circulation, u, verticesOnEdge, dcEdge, nEdges, nVertices, nVertLevels)
+      implicit none
+      real (kind=RKIND), dimension(:,:),intent(inout) :: circulation
+      real (kind=RKIND), dimension(:,:), intent(in) :: u
+      real (kind=RKIND), dimension(:), intent(in) :: dcEdge
+      integer, dimension(:,:), intent(in) :: verticesOnEdge
+      integer, intent(in) :: nEdges, nVertices, nVertLevels
+      integer :: iEdge, k
+
+      do iEdge=1,nEdges
+         if (verticesOnEdge(1,iEdge) &lt;= nVertices) then
+            do k=1,nVertLevels
+               circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
+            end do
+         end if
+         if (verticesOnEdge(2,iEdge) &lt;= nVertices) then
+            do k=1,nVertLevels
+               circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
+            end do
+         end if
+      end do
+   end subroutine compute_circulation
+
+   subroutine compute_vorticity(vorticity, circulation, areaTriangle, nVertices, nVertLevels)
+      implicit none
+      real (kind=RKIND), dimension(:,:),intent(inout) :: vorticity
+      real (kind=RKIND), dimension(:,:),intent(in) :: circulation
+      real (kind=RKIND), dimension(:), intent(in) :: areaTriangle
+      integer, intent(in) :: nVertices, nVertLevels
+      integer :: iVertex, k
+
+      do iVertex=1,nVertices
+         do k=1,nVertLevels
+            vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+         end do
+      end do
+
+   end subroutine compute_vorticity
+
+end module differential_operators

Added: branches/operator_addition/src/operators/module_geometry.F
===================================================================
--- branches/operator_addition/src/operators/module_geometry.F                                (rev 0)
+++ branches/operator_addition/src/operators/module_geometry.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -0,0 +1,237 @@
+module geometry
+
+   contains
+
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
+   !   sphere with given radius.
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   real function sphere_distance(lat1, lon1, lat2, lon2, radius)
+
+      implicit none
+
+      real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
+
+      real (kind=RKIND) :: arg1
+
+      arg1 = sqrt( sin(0.5*(lat2-lat1))**2 +  &amp;
+                   cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
+      sphere_distance = 2.*radius*asin(arg1)
+
+   end function sphere_distance
+
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! FUNCTION SPHERE_ANGLE
+   !
+   ! Computes the angle between arcs AB and AC, given points A, B, and C
+   ! Equation numbers w.r.t. http://mathworld.wolfram.com/SphericalTrigonometry.html
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   real function sphere_angle(ax, ay, az, bx, by, bz, cx, cy, cz)
+   
+      implicit none
+   
+      real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz
+   
+      real (kind=RKIND) :: a, b, c          ! Side lengths of spherical triangle ABC
+   
+      real (kind=RKIND) :: ABx, ABy, ABz    ! The components of the vector AB
+      real (kind=RKIND) :: mAB              ! The magnitude of AB
+      real (kind=RKIND) :: ACx, ACy, ACz    ! The components of the vector AC
+      real (kind=RKIND) :: mAC              ! The magnitude of AC
+   
+      real (kind=RKIND) :: Dx               ! The i-components of the cross product AB x AC
+      real (kind=RKIND) :: Dy               ! The j-components of the cross product AB x AC
+      real (kind=RKIND) :: Dz               ! The k-components of the cross product AB x AC
+   
+      real (kind=RKIND) :: s                ! Semiperimeter of the triangle
+      real (kind=RKIND) :: sin_angle
+   
+      a = acos(max(min(bx*cx + by*cy + bz*cz,1.0),-1.0))      ! Eqn. (3)
+      b = acos(max(min(ax*cx + ay*cy + az*cz,1.0),-1.0))      ! Eqn. (2)
+      c = acos(max(min(ax*bx + ay*by + az*bz,1.0),-1.0))      ! Eqn. (1)
+   
+      ABx = bx - ax
+      ABy = by - ay
+      ABz = bz - az
+   
+      ACx = cx - ax
+      ACy = cy - ay
+      ACz = cz - az
+   
+      Dx =   (ABy * ACz) - (ABz * ACy)
+      Dy = -((ABx * ACz) - (ABz * ACx))
+      Dz =   (ABx * ACy) - (ABy * ACx)
+   
+      s = 0.5*(a + b + c)
+!      sin_angle = sqrt((sin(s-b)*sin(s-c))/(sin(b)*sin(c)))   ! Eqn. (28)
+      sin_angle = sqrt(min(1.,max(0.,(sin(s-b)*sin(s-c))/(sin(b)*sin(c)))))   ! Eqn. (28)
+   
+      if ((Dx*ax + Dy*ay + Dz*az) &gt;= 0.0) then
+         sphere_angle =  2.0 * asin(max(min(sin_angle,1.0),-1.0))
+      else
+         sphere_angle = -2.0 * asin(max(min(sin_angle,1.0),-1.0))
+      end if
+   
+   end function sphere_angle                   
+         
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! FUNCTION PLANE_ANGLE
+   !
+   ! Computes the angle between vectors AB and AC, given points A, B, and C, and
+   !   a vector (u,v,w) normal to the plane.
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   real function plane_angle(ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w)
+   
+      implicit none
+   
+      real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz, cx, cy, cz, u, v, w
+   
+      real (kind=RKIND) :: ABx, ABy, ABz    ! The components of the vector AB
+      real (kind=RKIND) :: mAB              ! The magnitude of AB
+      real (kind=RKIND) :: ACx, ACy, ACz    ! The components of the vector AC
+      real (kind=RKIND) :: mAC              ! The magnitude of AC
+   
+      real (kind=RKIND) :: Dx               ! The i-components of the cross product AB x AC
+      real (kind=RKIND) :: Dy               ! The j-components of the cross product AB x AC
+      real (kind=RKIND) :: Dz               ! The k-components of the cross product AB x AC
+   
+      real (kind=RKIND) :: cos_angle
+   
+      ABx = bx - ax
+      ABy = by - ay
+      ABz = bz - az
+      mAB = sqrt(ABx**2.0 + ABy**2.0 + ABz**2.0)
+   
+      ACx = cx - ax
+      ACy = cy - ay
+      ACz = cz - az
+      mAC = sqrt(ACx**2.0 + ACy**2.0 + ACz**2.0)
+   
+   
+      Dx =   (ABy * ACz) - (ABz * ACy)
+      Dy = -((ABx * ACz) - (ABz * ACx))
+      Dz =   (ABx * ACy) - (ABy * ACx)
+   
+      cos_angle = (ABx*ACx + ABy*ACy + ABz*ACz) / (mAB * mAC)
+   
+      if ((Dx*u + Dy*v + Dz*w) &gt;= 0.0) then
+         plane_angle =  acos(max(min(cos_angle,1.0),-1.0))
+      else
+         plane_angle = -acos(max(min(cos_angle,1.0),-1.0))
+      end if
+   
+   end function plane_angle
+
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! FUNCTION ARC_LENGTH
+   !
+   ! Returns the length of the great circle arc from A=(ax, ay, az) to 
+   !    B=(bx, by, bz). It is assumed that both A and B lie on the surface of the
+   !    same sphere centered at the origin.
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   real function arc_length(ax, ay, az, bx, by, bz)
+   
+      implicit none
+   
+      real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
+   
+      real (kind=RKIND) :: r, c
+      real (kind=RKIND) :: cx, cy, cz
+   
+      cx = bx - ax
+      cy = by - ay
+      cz = bz - az
+
+!      r = ax*ax + ay*ay + az*az
+!      c = cx*cx + cy*cy + cz*cz
+!
+!      arc_length = sqrt(r) * acos(1.0 - c/(2.0*r))
+
+      r = sqrt(ax*ax + ay*ay + az*az)
+      c = sqrt(cx*cx + cy*cy + cz*cz)
+!      arc_length = sqrt(r) * 2.0 * asin(c/(2.0*r))
+      arc_length = r * 2.0 * asin(c/(2.0*r))
+
+   end function arc_length
+   
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! SUBROUTINE ARC_BISECT
+   !
+   ! Returns the point C=(cx, cy, cz) that bisects the great circle arc from
+   !   A=(ax, ay, az) to B=(bx, by, bz). It is assumed that A and B lie on the
+   !   surface of a sphere centered at the origin.
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   subroutine arc_bisect(ax, ay, az, bx, by, bz, cx, cy, cz)
+   
+      implicit none
+   
+      real (kind=RKIND), intent(in) :: ax, ay, az, bx, by, bz
+      real (kind=RKIND), intent(out) :: cx, cy, cz
+   
+      real (kind=RKIND) :: r           ! Radius of the sphere
+      real (kind=RKIND) :: d           
+   
+      r = sqrt(ax*ax + ay*ay + az*az)
+   
+      cx = 0.5*(ax + bx)
+      cy = 0.5*(ay + by)
+      cz = 0.5*(az + bz)
+   
+      if (cx == 0. .and. cy == 0. .and. cz == 0.) then
+         write(0,*) 'Error: arc_bisect: A and B are diametrically opposite'
+      else
+         d = sqrt(cx*cx + cy*cy + cz*cz)
+         cx = r * cx / d
+         cy = r * cy / d
+         cz = r * cz / d
+      end if
+   
+   end subroutine arc_bisect
+
+   subroutine get_distance(x1,x2,xout)
+     implicit none
+     real(kind=RKIND), intent(in) :: x1(3), x2(3)
+     real(kind=RKIND), intent(out) :: xout
+     xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
+   end subroutine get_distance
+
+   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
+!======================================================================
+ 
+end module geometry

Added: branches/operator_addition/src/operators/module_linear_algebra.F
===================================================================
--- branches/operator_addition/src/operators/module_linear_algebra.F                                (rev 0)
+++ branches/operator_addition/src/operators/module_linear_algebra.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -0,0 +1,138 @@
+module linear_algebra
+
+  contains
+
+! 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, 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
+      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, 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
+!
+  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)))
+    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
+
+   subroutine get_dotproduct(x1,x2,xout)
+     implicit none
+     real(kind=RKIND), intent(in) :: x1(3), x2(3)
+     real(kind=RKIND), intent(out) :: xout
+     xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
+   end subroutine get_dotproduct
+ 
+ end module linear_algebra

Modified: branches/operator_addition/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/operator_addition/src/operators/module_vector_reconstruction.F        2010-05-21 17:53:43 UTC (rev 300)
+++ branches/operator_addition/src/operators/module_vector_reconstruction.F        2010-05-24 14:49:07 UTC (rev 301)
@@ -3,6 +3,8 @@
    use grid_types
    use configure
    use constants
+   use linear_algebra
+   use geometry
 
    implicit none
 
@@ -254,30 +256,6 @@
 
    end subroutine reconstruct
 
-   subroutine get_distance(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
-   end subroutine get_distance
- 
-   subroutine get_dotproduct(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
-   end subroutine get_dotproduct
- 
- 
-   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
-
-
    subroutine evaluate_rbf(xin,xout)
      real(kind=RKIND), intent(in) ::  xin
      real(kind=RKIND), intent(out) :: xout
@@ -293,162 +271,4 @@
 
    end subroutine evaluate_rbf
 
-!======================================================================
-! 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.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, 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
-      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, 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
-!
-  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 vector_reconstruction

</font>
</pre>