<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 + &
- 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) <= 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) <= 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 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
- if(cell2 <= 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))) / &
- 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) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
- gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
- 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 + &
- 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) <= 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) <= 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 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
- if(cell2 <= 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))) / &
- 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) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
- gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
- 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))) / &
+ 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) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
+ do k = 1,nVertLevels
+ grad(k,iEdge) = (u_cell(k,cellsOnEdge(2,iEdge)) - u_cell(k,cellsOnEdge(1,iEdge))) / &
+ 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 <= nCells) then
+ do k=1,nVertLevels
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ enddo
+ endif
+ if(cell2 <= 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) <= 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) <= 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 + &
+ 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) >= 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) >= 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) :: &
+ p_1 (3), &
+ p_2 (3)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real (kind=RKIND), intent(out) :: &
+ 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, "An Introduction to Computational Physics," 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) :: &
- p_1 (3), &
- p_2 (3)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(out) :: &
- 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, "An Introduction to Computational Physics," 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>