<p><b>duda</b> 2010-05-27 13:41:43 -0600 (Thu, 27 May 2010)</p><p>BRANCH COMMIT <br>
<br>
Use routines from linear algebra and geometry operator modules<br>
in hydrostatic atmosphere core; remove local copies of these<br>
routines.<br>
<br>
<br>
M src/core_hyd_atmos/module_advection.F<br>
M src/core_hyd_atmos/module_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/operator_addition/src/core_hyd_atmos/module_advection.F
===================================================================
--- branches/operator_addition/src/core_hyd_atmos/module_advection.F 2010-05-27 15:19:53 UTC (rev 314)
+++ branches/operator_addition/src/core_hyd_atmos/module_advection.F 2010-05-27 19:41:43 UTC (rev 315)
@@ -3,6 +3,8 @@
use grid_types
use configure
use constants
+ use linear_algebra
+ use geometry
contains
@@ -329,178 +331,6 @@
end subroutine initialize_advection_rk
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! 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 poly_fit_2(a_in,b_out,weights_in,m,n,ne)
implicit none
@@ -558,131 +388,4 @@
end subroutine poly_fit_2
-
-! 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
-
end module advection
Modified: branches/operator_addition/src/core_hyd_atmos/module_test_cases.F
===================================================================
--- branches/operator_addition/src/core_hyd_atmos/module_test_cases.F 2010-05-27 15:19:53 UTC (rev 314)
+++ branches/operator_addition/src/core_hyd_atmos/module_test_cases.F 2010-05-27 19:41:43 UTC (rev 315)
@@ -3,6 +3,7 @@
use grid_types
use configure
use constants
+ use geometry
contains
@@ -466,25 +467,6 @@
end subroutine hyd_test_case_1
- 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
</font>
</pre>