|
new file 100644
|
|
|
!
|
|
|
! file gnbnaux.f
|
|
|
!
|
|
|
!
|
|
|
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
|
! . .
|
|
|
! . copyright (c) 2004 by UCAR .
|
|
|
! . .
|
|
|
! . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH .
|
|
|
! . .
|
|
|
! . all rights reserved .
|
|
|
! . .
|
|
|
! . .
|
|
|
! . FISHPACK version 5.0 .
|
|
|
! . .
|
|
|
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|
|
|
!
|
|
|
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
|
|
|
! * *
|
|
|
! * F I S H P A C K *
|
|
|
! * *
|
|
|
! * *
|
|
|
! * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF *
|
|
|
! * *
|
|
|
! * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS *
|
|
|
! * *
|
|
|
! * (Version 5.0 , JUNE 2004) *
|
|
|
! * *
|
|
|
! * BY *
|
|
|
! * *
|
|
|
! * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET *
|
|
|
! * *
|
|
|
! * OF *
|
|
|
! * *
|
|
|
! * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH *
|
|
|
! * *
|
|
|
! * BOULDER, COLORADO (80307) U.S.A. *
|
|
|
! * *
|
|
|
! * WHICH IS SPONSORED BY *
|
|
|
! * *
|
|
|
! * THE NATIONAL SCIENCE FOUNDATION *
|
|
|
! * *
|
|
|
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
|
|
|
!
|
|
|
!
|
|
|
! PACKAGE GNBNAUX
|
|
|
!
|
|
|
! LATEST REVISION June 2004
|
|
|
!
|
|
|
! PURPOSE TO PROVIDE AUXILIARY ROUTINES FOR FISHPACK
|
|
|
! ENTRIES GENBUN AND POISTG.
|
|
|
!
|
|
|
! USAGE THERE ARE NO USER ENTRIES IN THIS PACKAGE.
|
|
|
! THE ROUTINES IN THIS PACKAGE ARE NOT INTENDED
|
|
|
! TO BE CALLED BY USERS, BUT RATHER BY ROUTINES
|
|
|
! IN PACKAGES GENBUN AND POISTG.
|
|
|
!
|
|
|
! SPECIAL CONDITIONS NONE
|
|
|
!
|
|
|
! I/O NONE
|
|
|
!
|
|
|
! PRECISION SINGLE
|
|
|
!
|
|
|
!
|
|
|
! LANGUAGE FORTRAN 90
|
|
|
!
|
|
|
! HISTORY WRITTEN IN 1979 BY ROLAND SWEET OF NCAR'S
|
|
|
! SCIENTIFIC COMPUTING DIVISION. MADE AVAILABLE
|
|
|
! ON NCAR'S PUBLIC LIBRARIES IN JANUARY, 1980.
|
|
|
! Revised by John Adams in June 2004 incorporating
|
|
|
! Fortran 90 features
|
|
|
!
|
|
|
! PORTABILITY FORTRAN 90
|
|
|
! ********************************************************************
|
|
|
SUBROUTINE COSGEN(N, IJUMP, FNUM, FDEN, vecA, IA)
|
|
|
implicit none
|
|
|
|
|
|
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, &
|
|
|
ONE = 1.0D0, TWO = 2.0D0, &
|
|
|
FOUR = 4.0D0
|
|
|
|
|
|
!-----------------------------------------------
|
|
|
! D u m m y A r g u m e n t s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER, INTENT(IN) :: N, IJUMP, IA
|
|
|
DOUBLE PRECISION, INTENT(IN) :: FNUM, FDEN
|
|
|
DOUBLE PRECISION, DIMENSION(IA),INTENT(OUT) :: vecA
|
|
|
!-----------------------------------------------
|
|
|
! L o c a l V a r i a b l e s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER :: K3, K4, K, K1, K5, I, K2, NP1
|
|
|
DOUBLE PRECISION :: PI, PIBYN, X, Y
|
|
|
!-----------------------------------------------
|
|
|
!
|
|
|
!
|
|
|
! THIS SUBROUTINE COMPUTES REQUIRED COSINE VALUES IN ASCENDING
|
|
|
! ORDER. WHEN IJUMP .GT. 1 THE ROUTINE COMPUTES VALUES
|
|
|
!
|
|
|
! 2*COS(J*PI/L) , J=1,2,...,L AND J .NE. 0(MOD N/IJUMP+1)
|
|
|
!
|
|
|
! WHERE L = IJUMP*(N/IJUMP+1).
|
|
|
!
|
|
|
!
|
|
|
! WHEN IJUMP = 1 IT COMPUTES
|
|
|
!
|
|
|
! 2*COS((J-FNUM)*PI/(N+FDEN)) , J=1, 2, ... ,N
|
|
|
!
|
|
|
! WHERE
|
|
|
! FNUM = 0.5, FDEN = 0.0, FOR REGULAR REDUCTION VALUES
|
|
|
! FNUM = 0.0, FDEN = 1.0, FOR B-R AND C-R WHEN ISTAG = 1
|
|
|
! FNUM = 0.0, FDEN = 0.5, FOR B-R AND C-R WHEN ISTAG = 2
|
|
|
! FNUM = 0.5, FDEN = 0.5, FOR B-R AND C-R WHEN ISTAG = 2
|
|
|
! IN POISN2 ONLY.
|
|
|
!
|
|
|
!
|
|
|
PI = 4.0*ATAN(1.0)
|
|
|
IF (N /= 0) THEN
|
|
|
IF (IJUMP /= 1) THEN
|
|
|
K3 = N/IJUMP + 1
|
|
|
K4 = K3 - 1
|
|
|
PIBYN = PI/FLOAT(N + IJUMP)
|
|
|
DO K = 1, IJUMP
|
|
|
K1 = (K - 1)*K3
|
|
|
K5 = (K - 1)*K4
|
|
|
DO I = 1, K4
|
|
|
X = K1 + I
|
|
|
K2 = K5 + I
|
|
|
vecA(K2) = -2.*COS(X*PIBYN)
|
|
|
END DO
|
|
|
END DO
|
|
|
ELSE
|
|
|
NP1 = N + 1
|
|
|
Y = PI/(FLOAT(N) + FDEN)
|
|
|
DO I = 1, N
|
|
|
X = FLOAT(NP1 - I) - FNUM
|
|
|
vecA(I) = 2.*COS(X*Y)
|
|
|
END DO
|
|
|
ENDIF
|
|
|
ENDIF
|
|
|
!
|
|
|
END SUBROUTINE COSGEN
|
|
|
|
|
|
SUBROUTINE MERGE(TCOS, I1, M1, I2, M2, I3, itcos)
|
|
|
implicit none
|
|
|
|
|
|
|
|
|
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, &
|
|
|
ONE = 1.0D0, TWO = 2.0D0, &
|
|
|
FOUR = 4.0D0
|
|
|
|
|
|
!-----------------------------------------------
|
|
|
! D u m m y A r g u m e n t s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER, INTENT(IN) :: I1, M1, I2, M2, I3, ITCOS
|
|
|
DOUBLE PRECISION, DIMENSION(ITCOS), INTENT(INOUT) :: TCOS
|
|
|
!-----------------------------------------------
|
|
|
! L o c a l V a r i a b l e s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER :: J11, J3, J1, J2, J, L, K, M
|
|
|
DOUBLE PRECISION :: X, Y
|
|
|
!-----------------------------------------------
|
|
|
!
|
|
|
! THIS SUBROUTINE MERGES TWO ASCENDING STRINGS OF NUMBERS IN THE
|
|
|
! ARRAY TCOS. THE FIRST STRING IS OF LENGTH M1 AND STARTS AT
|
|
|
! TCOS(I1+1). THE SECOND STRING IS OF LENGTH M2 AND STARTS AT
|
|
|
! TCOS(I2+1). THE MERGED STRING GOES INTO TCOS(I3+1).
|
|
|
!
|
|
|
!
|
|
|
J1 = 1
|
|
|
J2 = 1
|
|
|
J = I3
|
|
|
IF (M1 == 0) GO TO 107
|
|
|
IF (M2 == 0) GO TO 104
|
|
|
101 CONTINUE
|
|
|
J11 = J1
|
|
|
J3 = MAX(M1,J11)
|
|
|
DO J1 = J11, J3
|
|
|
J = J + 1
|
|
|
L = J1 + I1
|
|
|
X = TCOS(L)
|
|
|
L = J2 + I2
|
|
|
Y = TCOS(L)
|
|
|
IF (X - Y > 0.) GO TO 103
|
|
|
TCOS(J) = X
|
|
|
END DO
|
|
|
GO TO 106
|
|
|
103 CONTINUE
|
|
|
TCOS(J) = Y
|
|
|
J2 = J2 + 1
|
|
|
IF (J2 <= M2) GO TO 101
|
|
|
IF (J1 > M1) GO TO 109
|
|
|
104 CONTINUE
|
|
|
K = J - J1 + 1
|
|
|
DO J = J1, M1
|
|
|
M = K + J
|
|
|
L = J + I1
|
|
|
TCOS(M) = TCOS(L)
|
|
|
END DO
|
|
|
GO TO 109
|
|
|
106 CONTINUE
|
|
|
IF (J2 > M2) GO TO 109
|
|
|
107 CONTINUE
|
|
|
K = J - J2 + 1
|
|
|
DO J = J2, M2
|
|
|
M = K + J
|
|
|
L = J + I2
|
|
|
TCOS(M) = TCOS(L)
|
|
|
END DO
|
|
|
109 CONTINUE
|
|
|
!
|
|
|
END SUBROUTINE MERGE
|
|
|
|
|
|
|
|
|
SUBROUTINE TRIX(IDEGBR, IDEGCR, M, vecA, vecB, vecC, vecY, TCOS, ITCOS, &
|
|
|
vecD, vecW)
|
|
|
implicit none
|
|
|
|
|
|
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, &
|
|
|
ONE = 1.0D0, TWO = 2.0D0, &
|
|
|
FOUR = 4.0D0
|
|
|
|
|
|
!-----------------------------------------------
|
|
|
! D u m m y A r g u m e n t s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER, INTENT(IN) :: IDEGBR, IDEGCR, M, ITCOS
|
|
|
DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA, vecB, vecC
|
|
|
DOUBLE PRECISION, DIMENSION(ITCOS), INTENT(IN) :: TCOS
|
|
|
DOUBLE PRECISION, DIMENSION(M), INTENT(INOUT) :: vecY, vecD, vecW
|
|
|
!-----------------------------------------------
|
|
|
! L o c a l V a r i a b l e s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER :: MM1, IFB, IFC, L, LINT, K, I, IP
|
|
|
DOUBLE PRECISION :: X, XX, Z
|
|
|
!-----------------------------------------------
|
|
|
!
|
|
|
! SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS WHERE THE
|
|
|
! COEFFICIENT MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY
|
|
|
! TRIDIAGONAL ( . . . , vecA(I), vecB(I), vecC(I), . . . ).
|
|
|
!
|
|
|
MM1 = M - 1
|
|
|
IFB = IDEGBR + 1
|
|
|
IFC = IDEGCR + 1
|
|
|
L = IFB/IFC
|
|
|
LINT = 1
|
|
|
DO K = 1, IDEGBR
|
|
|
X = TCOS(K)
|
|
|
IF (K == L) THEN
|
|
|
I = IDEGBR + LINT
|
|
|
XX = X - TCOS(I)
|
|
|
vecW(:M) = vecY(:M)
|
|
|
vecY(:M) = XX*vecY(:M)
|
|
|
ENDIF
|
|
|
Z = 1./(vecB(1)-X)
|
|
|
vecD(1) = vecC(1)*Z
|
|
|
vecY(1) = vecY(1)*Z
|
|
|
DO I = 2, MM1
|
|
|
Z = 1./(vecB(I)-X-vecA(I)*vecD(I-1))
|
|
|
vecD(I) = vecC(I)*Z
|
|
|
vecY(I) = (vecY(I)-vecA(I)*vecY(I-1))*Z
|
|
|
END DO
|
|
|
Z = vecB(M) - X - vecA(M)*vecD(MM1)
|
|
|
IF (Z == 0.) THEN
|
|
|
vecY(M) = 0.
|
|
|
ELSE
|
|
|
vecY(M) = (vecY(M)-vecA(M)*vecY(MM1))/Z
|
|
|
ENDIF
|
|
|
DO IP = 1, MM1
|
|
|
vecY(M-IP) = vecY(M-IP) - vecD(M-IP)*vecY(M+1-IP)
|
|
|
END DO
|
|
|
IF (K /= L) CYCLE
|
|
|
vecY(:M) = vecY(:M) + vecW(:M)
|
|
|
LINT = LINT + 1
|
|
|
L = (LINT*IFB)/IFC
|
|
|
END DO
|
|
|
!
|
|
|
END SUBROUTINE TRIX
|
|
|
|
|
|
|
|
|
SUBROUTINE TRI3(M, vecA, vecB, vecC, ivecK, vecY1, vecY2, vecY3, TCOS, &
|
|
|
ITCOS, vecD, vecW1, vecW2, vecW3)
|
|
|
implicit none
|
|
|
|
|
|
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, &
|
|
|
ONE = 1.0D0, TWO = 2.0D0, &
|
|
|
FOUR = 4.0D0
|
|
|
|
|
|
!-----------------------------------------------
|
|
|
! D u m m y A r g u m e n t s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER, INTENT(IN) :: M,ITCOS
|
|
|
INTEGER, DIMENSION(4),INTENT(IN) :: ivecK
|
|
|
DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA, vecB, vecC
|
|
|
DOUBLE PRECISION, DIMENSION(ITCOS), INTENT(IN) :: TCOS
|
|
|
DOUBLE PRECISION, DIMENSION(M), INTENT(INOUT) :: vecY1, vecY2, vecY3, &
|
|
|
vecD, vecW1, vecW2, vecW3
|
|
|
!-----------------------------------------------
|
|
|
! L o c a l V a r i a b l e s
|
|
|
!-----------------------------------------------
|
|
|
INTEGER :: MM1, K1, K2, K3, K4, IF1, IF2, IF3, IF4, K2K3K4, &
|
|
|
L1, L2, L3, LINT1, LINT2, LINT3, KINT1, KINT2, KINT3, &
|
|
|
N, I, IP
|
|
|
DOUBLE PRECISION :: X, Z, XX
|
|
|
!-----------------------------------------------
|
|
|
!
|
|
|
! SUBROUTINE TO SOLVE THREE LINEAR SYSTEMS WHOSE COMMON COEFFICIENT
|
|
|
! MATRIX IS A RATIONAL FUNCTION IN THE MATRIX GIVEN BY
|
|
|
!
|
|
|
! TRIDIAGONAL (...,vecA(I),vecB(I),vecC(I),...)
|
|
|
!
|
|
|
MM1 = M - 1
|
|
|
K1 = ivecK(1)
|
|
|
K2 = ivecK(2)
|
|
|
K3 = ivecK(3)
|
|
|
K4 = ivecK(4)
|
|
|
IF1 = K1 + 1
|
|
|
IF2 = K2 + 1
|
|
|
IF3 = K3 + 1
|
|
|
IF4 = K4 + 1
|
|
|
K2K3K4 = K2 + K3 + K4
|
|
|
IF (K2K3K4 /= 0) THEN
|
|
|
L1 = IF1/IF2
|
|
|
L2 = IF1/IF3
|
|
|
L3 = IF1/IF4
|
|
|
LINT1 = 1
|
|
|
LINT2 = 1
|
|
|
LINT3 = 1
|
|
|
KINT1 = K1
|
|
|
KINT2 = KINT1 + K2
|
|
|
KINT3 = KINT2 + K3
|
|
|
ELSE
|
|
|
write(*,*) 'warning tri3: l1,l2,l3,kint1,kint2,kint3 uninitialized'
|
|
|
stop 'stop in tri3: l1,l2,l3,kint1,kint2,kint3 uninitialized'
|
|
|
ENDIF
|
|
|
DO N = 1, K1
|
|
|
X = TCOS(N)
|
|
|
IF (K2K3K4 /= 0) THEN
|
|
|
IF (N == L1) THEN
|
|
|
vecW1(:M) = vecY1(:M)
|
|
|
ENDIF
|
|
|
IF (N == L2) THEN
|
|
|
vecW2(:M) = vecY2(:M)
|
|
|
ENDIF
|
|
|
IF (N == L3) THEN
|
|
|
vecW3(:M) = vecY3(:M)
|
|
|
ENDIF
|
|
|
ENDIF
|
|
|
Z = 1./(vecB(1)-X)
|
|
|
vecD(1) = vecC(1)*Z
|
|
|
vecY1(1) = vecY1(1)*Z
|
|
|
vecY2(1) = vecY2(1)*Z
|
|
|
vecY3(1) = vecY3(1)*Z
|
|
|
DO I = 2, M
|
|
|
Z = 1./(vecB(I)-X-vecA(I)*vecD(I-1))
|
|
|
vecD(I) = vecC(I)*Z
|
|
|
vecY1(I) = (vecY1(I)-vecA(I)*vecY1(I-1))*Z
|
|
|
vecY2(I) = (vecY2(I)-vecA(I)*vecY2(I-1))*Z
|
|
|
vecY3(I) = (vecY3(I)-vecA(I)*vecY3(I-1))*Z
|
|
|
END DO
|
|
|
DO IP = 1, MM1
|
|
|
vecY1(M-IP) = vecY1(M-IP) - vecD(M-IP)*vecY1(M+1-IP)
|
|
|
vecY2(M-IP) = vecY2(M-IP) - vecD(M-IP)*vecY2(M+1-IP)
|
|
|
vecY3(M-IP) = vecY3(M-IP) - vecD(M-IP)*vecY3(M+1-IP)
|
|
|
END DO
|
|
|
IF (K2K3K4 == 0) CYCLE
|
|
|
IF (N == L1) THEN
|
|
|
I = LINT1 + KINT1
|
|
|
XX = X - TCOS(I)
|
|
|
vecY1(:M) = XX*vecY1(:M) + vecW1(:M)
|
|
|
LINT1 = LINT1 + 1
|
|
|
L1 = (LINT1*IF1)/IF2
|
|
|
ENDIF
|
|
|
IF (N == L2) THEN
|
|
|
I = LINT2 + KINT2
|
|
|
XX = X - TCOS(I)
|
|
|
vecY2(:M) = XX*vecY2(:M) + vecW2(:M)
|
|
|
LINT2 = LINT2 + 1
|
|
|
L2 = (LINT2*IF1)/IF3
|
|
|
ENDIF
|
|
|
IF (N /= L3) CYCLE
|
|
|
I = LINT3 + KINT3
|
|
|
XX = X - TCOS(I)
|
|
|
vecY3(:M) = XX*vecY3(:M) + vecW3(:M)
|
|
|
LINT3 = LINT3 + 1
|
|
|
L3 = (LINT3*IF1)/IF4
|
|
|
END DO
|
|
|
RETURN
|
|
|
!
|
|
|
! REVISION HISTORY---
|
|
|
!
|
|
|
! SEPTEMBER 1973 VERSION 1
|
|
|
! APRIL 1976 VERSION 2
|
|
|
! JANUARY 1978 VERSION 3
|
|
|
! DECEMBER 1979 VERSION 3.1
|
|
|
! OCTOBER 1980 CHANGED SEVERAL DIVIDES OF FLOATING INTEGERS
|
|
|
! TO INTEGER DIVIDES TO ACCOMODATE CRAY-1 ARITHMETIC.
|
|
|
! FEBRUARY 1985 DOCUMENTATION UPGRADE
|
|
|
! NOVEMBER 1988 VERSION 3.2, FORTRAN 77 CHANGES
|
|
|
!-----------------------------------------------------------------------
|
|
|
END SUBROUTINE TRI3
|