diff --git a/fish90/src/gnbnaux.f90 b/fish90/src/gnbnaux.f90 new file mode 100644 index 0000000000000000000000000000000000000000..4640121ecf744e48588a3222b58c153710f7eb37 --- /dev/null +++ b/fish90/src/gnbnaux.f90 @@ -0,0 +1,399 @@ +! +! 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