! ! 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