diff --git a/fish90/src/genbunal.f90 b/fish90/src/genbunal.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c4d025d8c880227ef1b1b6bbc5f3c28ac9d6d4ac
--- /dev/null
+++ b/fish90/src/genbunal.f90
@@ -0,0 +1,1325 @@
+!
+!     file genbun.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                  *
+!     *                                                               *
+!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+!
+!     SUBROUTINE GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR)
+!
+!
+! DIMENSION OF           A(M),B(M),C(M),Y(IDIMY,N)
+! ARGUMENTS
+!
+! LATEST REVISION        JUNE 2004
+!
+! PURPOSE                THE NAME OF THIS PACKAGE IS A MNEMONIC FOR THE
+!                        GENERALIZED BUNEMAN ALGORITHM.
+!
+!                        IT SOLVES THE REAL LINEAR SYSTEM OF EQUATIONS
+!
+!                        A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J)
+!                        + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J)
+!
+!                        FOR I = 1,2,...,M  AND  J = 1,2,...,N.
+!
+!                        INDICES I+1 AND I-1 ARE EVALUATED MODULO M,
+!                        I.E., X(0,J) = X(M,J) AND X(M+1,J) = X(1,J),
+!                        AND X(I,0) MAY EQUAL 0, X(I,2), OR X(I,N),
+!                        AND X(I,N+1) MAY EQUAL 0, X(I,N-1), OR X(I,1)
+!                        DEPENDING ON AN INPUT PARAMETER.
+!
+! USAGE                  CALL GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,
+!                                     IERROR)
+!
+! ARGUMENTS
+!
+! ON INPUT               NPEROD
+!
+!                          INDICATES THE VALUES THAT X(I,0) AND
+!                          X(I,N+1) ARE ASSUMED TO HAVE.
+!
+!                          = 0  IF X(I,0) = X(I,N) AND X(I,N+1) =
+!                               X(I,1).
+!                          = 1  IF X(I,0) = X(I,N+1) = 0  .
+!                          = 2  IF X(I,0) = 0 AND X(I,N+1) = X(I,N-1).
+!                          = 3  IF X(I,0) = X(I,2) AND X(I,N+1) =
+!                               X(I,N-1).
+!                          = 4  IF X(I,0) = X(I,2) AND X(I,N+1) = 0.
+!
+!                        N
+!                          THE NUMBER OF UNKNOWNS IN THE J-DIRECTION.
+!                          N MUST BE GREATER THAN 2.
+!
+!                        MPEROD
+!                          = 0 IF A(1) AND C(M) ARE NOT ZERO
+!                          = 1 IF A(1) = C(M) = 0
+!
+!                        M
+!                          THE NUMBER OF UNKNOWNS IN THE I-DIRECTION.
+!                          N MUST BE GREATER THAN 2.
+!
+!                        A,B,C
+!                          ONE-DIMENSIONAL ARRAYS OF LENGTH M THAT
+!                          SPECIFY THE COEFFICIENTS IN THE LINEAR
+!                          EQUATIONS GIVEN ABOVE.  IF MPEROD = 0
+!                          THE ARRAY ELEMENTS MUST NOT DEPEND UPON
+!                          THE INDEX I, BUT MUST BE CONSTANT.
+!                          SPECIFICALLY, THE SUBROUTINE CHECKS THE
+!                          FOLLOWING CONDITION .
+!
+!                            A(I) = C(1)
+!                            C(I) = C(1)
+!                            B(I) = B(1)
+!
+!                          FOR I=1,2,...,M.
+!
+!                        IDIMY
+!                          THE ROW (OR FIRST) DIMENSION OF THE
+!                          TWO-DIMENSIONAL ARRAY Y AS IT APPEARS
+!                          IN THE PROGRAM CALLING GENBUN.
+!                          THIS PARAMETER IS USED TO SPECIFY THE
+!                          VARIABLE DIMENSION OF Y.
+!                          IDIMY MUST BE AT LEAST M.
+!
+!                        Y
+!                          A TWO-DIMENSIONAL COMPLEX ARRAY THAT
+!                          SPECIFIES THE VALUES OF THE RIGHT SIDE
+!                          OF THE LINEAR SYSTEM OF EQUATIONS GIVEN
+!                          ABOVE.
+!                          Y MUST BE DIMENSIONED AT LEAST M*N.
+!
+!
+!  ON OUTPUT             Y
+!
+!                          CONTAINS THE SOLUTION X.
+!
+!                        IERROR
+!                          AN ERROR FLAG WHICH INDICATES INVALID
+!                          INPUT PARAMETERS  EXCEPT FOR NUMBER
+!                          ZERO, A SOLUTION IS NOT ATTEMPTED.
+!
+!                          = 0  NO ERROR.
+!                          = 1  M .LE. 2  .
+!                          = 2  N .LE. 2
+!                          = 3  IDIMY .LT. M
+!                          = 4  NPEROD .LT. 0 OR NPEROD .GT. 4
+!                          = 5  MPEROD .LT. 0 OR MPEROD .GT. 1
+!                          = 6  A(I) .NE. C(1) OR C(I) .NE. C(1) OR
+!                               B(I) .NE. B(1) FOR
+!                               SOME I=1,2,...,M.
+!                          = 7  A(1) .NE. 0 OR C(M) .NE. 0 AND
+!                                 MPEROD = 1
+!                          = 20 If the dynamic allocation of real and
+!                               complex work space required for solution
+!                               fails (for example if N,M are too large
+!                               for your computer)
+!
+!
+! SPECIAL CONDITONS      NONE
+!
+! I/O                    NONE
+!
+! PRECISION              SINGLE
+!
+! REQUIRED FILES         comf.f,gnbnaux.f,fish.f
+! FILES
+!
+! 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 in June 2004 by John Adams using
+!                        Fortran 90 dynamically allocated work space.
+!
+! ALGORITHM              THE LINEAR SYSTEM IS SOLVED BY A CYCLIC
+!                        REDUCTION ALGORITHM DESCRIBED IN THE
+!                        REFERENCE.
+!
+! PORTABILITY            FORTRAN 90 --
+!                        THE MACHINE DEPENDENT CONSTANT PI IS
+!                        DEFINED IN FUNCTION PIMACH.
+!
+! REFERENCES             SWEET, R., "A CYCLIC REDUCTION ALGORITHM FOR
+!                        SOLVING BLOCK TRIDIAGONAL SYSTEMS OF ARBITRARY
+!                        DIMENSIONS," SIAM J. ON NUMER. ANAL., 14 (1977)
+!                        PP. 706-720.
+!
+! ACCURACY               THIS TEST WAS PERFORMED ON a platform with
+!                        64 bit floating point arithmetic.
+!                        A UNIFORM RANDOM NUMBER GENERATOR WAS USED
+!                        TO CREATE A SOLUTION ARRAY X FOR THE SYSTEM
+!                        GIVEN IN THE 'PURPOSE' DESCRIPTION ABOVE
+!                        WITH
+!                          A(I) = C(I) = -0.5*B(I) = 1, I=1,2,...,M
+!
+!                        AND, WHEN MPEROD = 1
+!
+!                          A(1) = C(M) = 0
+!                          A(M) = C(1) = 2.
+!
+!                        THE SOLUTION X WAS SUBSTITUTED INTO THE
+!                        GIVEN SYSTEM  AND, USING DOUBLE PRECISION
+!                        A RIGHT SIDE Y WAS COMPUTED.
+!                        USING THIS ARRAY Y, SUBROUTINE GENBUN
+!                        WAS CALLED TO PRODUCE APPROXIMATE
+!                        SOLUTION Z.  THEN RELATIVE ERROR
+!                          E = MAX(ABS(Z(I,J)-X(I,J)))/
+!                              MAX(ABS(X(I,J)))
+!                        WAS COMPUTED, WHERE THE TWO MAXIMA ARE TAKEN
+!                        OVER I=1,2,...,M AND J=1,...,N.
+!
+!                        THE VALUE OF E IS GIVEN IN THE TABLE
+!                        BELOW FOR SOME TYPICAL VALUES OF M AND N.
+!
+!                   M (=N)    MPEROD    NPEROD        E
+!                   ------    ------    ------      ------
+!
+!                     31        0         0         6.E-14
+!                     31        1         1         4.E-13
+!                     31        1         3         3.E-13
+!                     32        0         0         9.E-14
+!                     32        1         1         3.E-13
+!                     32        1         3         1.E-13
+!                     33        0         0         9.E-14
+!                     33        1         1         4.E-13
+!                     33        1         3         1.E-13
+!                     63        0         0         1.E-13
+!                     63        1         1         1.E-12
+!                     63        1         3         2.E-13
+!                     64        0         0         1.E-13
+!                     64        1         1         1.E-12
+!                     64        1         3         6.E-13
+!                     65        0         0         2.E-13
+!                     65        1         1         1.E-12
+!                     65        1         3         4.E-13
+! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+      SUBROUTINE GENBUN(NPEROD, N, MPEROD, M, A, B, C, IDIMY, Y, IERROR)
+      USE fish
+      implicit none
+      
+      DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, &
+                                     ONE  = 1.0D0, TWO  = 2.0D0, &
+                                     FOUR = 4.0D0
+
+      TYPE(fishworkspace) :: w
+!-----------------------------------------------
+!   D u m m y   A r g u m e n t s
+!-----------------------------------------------
+      INTEGER, INTENT(IN)  :: NPEROD, N, MPEROD, M, IDIMY
+      INTEGER, INTENT(OUT) :: IERROR
+      DOUBLE PRECISION, DIMENSION(:)  :: A, B, C
+      DOUBLE PRECISION, INTENT(INOUT) :: Y(IDIMY,*)
+!-----------------------------------------------
+!   L o c a l   V a r i a b l e s
+!-----------------------------------------------
+      INTEGER :: IRWK, ICWK
+!-----------------------------------------------
+
+      write(*,*) 'fish90: genbun'
+      IERROR = 0
+!     check input arguments
+      IF (M <= 2) then
+         ierror = 1
+         return
+      end if
+      IF (N <= 2) then
+         ierror = 2
+         return
+      end if
+      IF (IDIMY < M) then
+         ierror = 3
+         return
+      end if
+      IF (NPEROD<0 .OR. NPEROD>4) then
+         ierror = 4
+         return
+      end if
+      IF (MPEROD<0 .OR. MPEROD>1) then
+         ierror = 5
+         return
+      end if
+!     compute and allocate real work space for genbun
+      CALL GEN_SPACE (N, M, IRWK)
+      ICWK = 0
+      CALL ALLOCATFISH (IRWK, ICWK, W, IERROR)
+!     return if allocation failed (e.g., if n,m are too large)
+      IF (IERROR == 20)  THEN
+         write(*,*) 'error call ALLOCATFISH'
+         RETURN 
+      END IF
+      call genbunn(NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR,w%rew,IRWK)
+!     release allocated work space
+      CALL FISHFIN (W,IERROR)
+      IF (IERROR == 20)  THEN
+         write(*,*) 'error call FISHFIN'
+         RETURN 
+      END IF
+!
+      END SUBROUTINE GENBUN
+
+ 
+      SUBROUTINE GENBUNN(NPEROD,N,MPEROD,M,vecA,vecB,vecC,IDIMY,matY, &
+                         IERROR,vecW,IW)
+      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)    :: NPEROD, N, MPEROD, M, IDIMY,IW
+      INTEGER, INTENT(INOUT) :: IERROR
+      DOUBLE PRECISION, DIMENSION(M), INTENT(IN)   :: vecA, vecB, vecC
+      DOUBLE PRECISION, DIMENSION(IW), INTENT(OUT) :: vecW
+      DOUBLE PRECISION, DIMENSION(IDIMY,N),INTENT(INOUT) :: matY
+!-----------------------------------------------
+!   L o c a l   V a r i a b l e s
+!-----------------------------------------------
+      INTEGER :: I, MP1, IWBA, IWBB, IWBC, IWB2, IWB3, IWW1, IWW2, IWW3, &
+                 IWD, IWTCOS, IWP, IW2,  K, J, MP, NP, &
+                 IPSTOR, IREV, MH, MHM1, MODD, NBY2, MSKIP
+      DOUBLE PRECISION :: A1
+!-----------------------------------------------
+      IF (MPEROD /= 1) THEN
+         DO I = 2, M
+            IF (vecA(I) /= vecC(1)) GO TO 103
+            IF (vecC(I) /= vecC(1)) GO TO 103
+            IF (vecB(I) /= vecB(1)) GO TO 103
+         END DO
+         GO TO 104
+      ENDIF
+      IF (vecA(1)/=ZERO .OR. vecC(M)/=ZERO) IERROR = 7
+      GO TO 104
+  103 CONTINUE
+      IERROR = 6
+  104 CONTINUE
+      IF (IERROR /= 0) RETURN 
+      MP1 = M + 1
+      IWBA = MP1
+      IWBB = IWBA + M
+      IWBC = IWBB + M
+      IWB2 = IWBC + M
+      IWB3 = IWB2 + M
+      IWW1 = IWB3 + M
+      IWW2 = IWW1 + M
+      IWW3 = IWW2 + M
+      IWD = IWW3 + M
+      IWTCOS = IWD + M
+      IWP = IWTCOS + 4*N
+      vecW(IWBA:M-1+IWBA) = -vecA(:M)
+      vecW(IWBC:M-1+IWBC) = -vecC(:M)
+      vecW(IWBB:M-1+IWBB) = TWO - vecB(:M)
+      matY(:M,:N) = -matY(:M,:N)
+      MP = MPEROD + 1
+      NP = NPEROD + 1
+      GO TO (114,107) MP
+  107 CONTINUE
+      GO TO (108,109,110,111,123) NP
+  108 CONTINUE
+      IW2 = IW - IWP + 1
+      CALL POISP2 (M, N, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), &
+                   vecW(IWBC:IWBC+M-1), matY, IDIMY, &
+                   vecW(1:M), vecW(IWB2:IWB2+M-1), vecW(IWB3:IWB3+M-1), &
+                   vecW(IWW1:IWW1+M-1), vecW(IWW2:IWW2+M-1), &
+                   vecW(IWW3:IWW3+M-1), vecW(IWD:IWD+M-1), &
+                   vecW(IWTCOS:IWTCOS+4*N-1), vecW(IWP:),IW2)
+      GO TO 112
+  109 CONTINUE
+      IW2 = IW - IWP + 1
+      CALL POISD2 (M, N, 1, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1),  &
+                   vecW(IWBC::IWBC+M-1), matY, IDIMY, &
+                   vecW(1:M), vecW(IWW1:IWW1+M-1), vecW(IWD:IWD+M-1), &
+                   vecW(IWTCOS:IWTCOS+4*N-1), vecW(IWP:), IW2)
+      GO TO 112
+  110 CONTINUE
+      IW2 = IW - IWP + 1
+      CALL POISN2 (M, N, 1, 2, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), &
+                   vecW(IWBC:IWBC+M-1), matY, IDIMY, vecW(1:M),  &
+                   vecW(IWB2:IWB2+M-1), vecW(IWB3:IWB3+M-1), &
+                   vecW(IWW1:IWW1+M-1), vecW(IWW2:IWW2+M-1), &
+                   vecW(IWW3::IWW3+M-1), vecW(IWD:IWD+M-1), &
+                   vecW(IWTCOS:IWTCOS+4*N-1),vecW(IWP:),IW2)
+      GO TO 112
+  111 CONTINUE
+      IW2 = IW - IWP + 1
+      CALL POISN2 (M, N, 1, 1, vecW(IWBA:IWBA+M-1), vecW(IWBB:IWBB+M-1), &
+                   vecW(IWBC::IWBC+M-1), matY, IDIMY, vecW(1:M), &
+                   vecW(IWB2:IWW1+M-1), vecW(IWB3:IWB3+M-1), &
+                   vecW(IWW1:IWW1+M-1), vecW(IWW2:IWW2+M-1), &
+                   vecW(IWW3::IWW3+M-1), vecW(IWD:IWD+M-1), &
+                   vecW(IWTCOS:IWTCOS+4*N-1),vecW(IWP:),IW2)
+  112 CONTINUE
+      IPSTOR = vecW(IWW1)
+      IREV = 2
+      IF (NPEROD == 4) GO TO 124
+  113 CONTINUE
+      GO TO (127,133) MP
+  114 CONTINUE
+      MH = (M + 1)/2
+      MHM1 = MH - 1
+      MODD = 1
+      IF (MH*2 == M) MODD = 2
+      DO J = 1, N
+         vecW(:MHM1) = matY(MH-1:MH-MHM1:(-1),J) - matY(MH+1:MHM1+MH,J)
+         vecW(MH+1:MHM1+MH) = matY(MH-1:MH-MHM1:(-1),J) + matY(MH+1:MHM1+MH,J)
+         vecW(MH) = TWO*matY(MH,J)
+         GO TO (117,116) MODD
+  116    CONTINUE
+         vecW(M) = TWO*matY(M,J)
+  117    CONTINUE
+         matY(:M,J) = vecW(:M)
+      END DO
+      K = IWBC + MHM1 - 1
+      I = IWBA + MHM1
+      vecW(K) = ZERO
+      vecW(I) = ZERO
+      vecW(K+1) = TWO*vecW(K+1)
+      SELECT CASE (MODD) 
+      CASE DEFAULT
+         K = IWBB + MHM1 - 1
+         vecW(K) = vecW(K) - vecW(I-1)
+         vecW(IWBC-1) = vecW(IWBC-1) + vecW(IWBB-1)
+      CASE (2) 
+         vecW(IWBB-1) = vecW(K+1)
+      END SELECT
+      GO TO 107
+!
+!     REVERSE COLUMNS WHEN NPEROD = 4.
+!
+  123 CONTINUE
+      IREV = 1
+      NBY2 = N/2
+  124 CONTINUE
+      DO J = 1, NBY2
+         MSKIP = N + 1 - J
+         DO I = 1, M
+            A1 = matY(I,J)
+            matY(I,J) = matY(I,MSKIP)
+            matY(I,MSKIP) = A1
+         END DO
+      END DO
+      GO TO (110,113) IREV
+  127 CONTINUE
+      DO J = 1, N
+         vecW(MH-1:MH-MHM1:(-1)) = HALF*(matY(MH+1:MHM1+MH,J)+matY(:MHM1,J))
+         vecW(MH+1:MHM1+MH) = HALF*(matY(MH+1:MHM1+MH,J)-matY(:MHM1,J))
+         vecW(MH) = HALF*matY(MH,J)
+         GO TO (130,129) MODD
+  129    CONTINUE
+         vecW(M) = HALF*matY(M,J)
+  130    CONTINUE
+         matY(:M,J) = vecW(:M)
+      END DO
+  133 CONTINUE
+      vecW(1) = IPSTOR + IWP - 1
+      RETURN 
+      END SUBROUTINE GENBUNN
+
+
+      SUBROUTINE POISD2(MR,NR,ISTAG,vecBA,vecBB,vecBC,matY,IDIMY, &
+                        vecB,vecW,vecD,TCOS,vecP,IDIMP)
+      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) :: MR, NR, ISTAG, IDIMY,IDIMP
+      DOUBLE PRECISION, DIMENSION(MR), INTENT(IN)         :: vecBA, vecBB, &
+                                                             vecBC
+      DOUBLE PRECISION, DIMENSION(IDIMY,NR),INTENT(INOUT) :: matY
+      DOUBLE PRECISION, DIMENSION(MR), INTENT(INOUT)      :: vecB, vecD, vecW
+      DOUBLE PRECISION, DIMENSION(4*NR), INTENT(INOUT)    :: TCOS
+      DOUBLE PRECISION, DIMENSION(IDIMP), INTENT(INOUT)   :: vecP
+!-----------------------------------------------
+!   L o c a l   V a r i a b l e s
+!-----------------------------------------------
+      INTEGER :: M, N, JSH, IP, IPSTOR, KR, IRREG, JSTSAV, I, LR, NUN, &
+                 JST, JSP, L, NODD, J, JM1, JP1, JM2, JP2, JM3, JP3, &
+                 NODDPR, KRPI, IDEG, JDEG
+      DOUBLE PRECISION :: FI, T
+!-----------------------------------------------
+!
+!     SUBROUTINE TO SOLVE POISSON'S EQUATION FOR DIRICHLET BOUNDARY
+!     CONDITIONS.
+!
+!     ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A.
+!     ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A+I.
+!
+      M = MR
+      N = NR
+      JSH = 0
+      FI = ONE/FLOAT(ISTAG)
+      IP = -M
+      IPSTOR = 0
+      SELECT CASE (ISTAG) 
+      CASE DEFAULT
+         KR = 0
+         IRREG = 1
+         IF (N > 1) GO TO 106
+         TCOS(1) = ZERO
+      CASE (2) 
+         KR = 1
+         JSTSAV = 1
+         IRREG = 2
+         IF (N > 1) GO TO 106
+         TCOS(1) = -ONE
+      END SELECT
+      vecB(:M) = matY(:M,1)
+      CALL TRIX (1, 0, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+      matY(:M,1) = vecB(:M)
+      GO TO 183
+  106 CONTINUE
+      LR = 0
+      vecP(:M) = ZERO
+      NUN = N
+      JST = 1
+      JSP = N
+!
+!     IRREG = 1 WHEN NO IRREGULARITIES HAVE OCCURRED, OTHERWISE IT IS 2.
+!
+  108 CONTINUE
+      L = 2*JST
+      NODD = 2 - 2*((NUN + 1)/2) + NUN
+!
+!     NODD = 1 WHEN NUN IS ODD, OTHERWISE IT IS 2.
+!
+      SELECT CASE (NODD) 
+      CASE DEFAULT
+         JSP = JSP - L
+      CASE (1) 
+         JSP = JSP - JST
+         IF (IRREG /= 1) JSP = JSP - L
+      END SELECT
+      CALL COSGEN (JST, 1, HALF, ZERO, TCOS, 4*NR)
+      IF (L <= JSP) THEN
+         DO J = L, JSP, L
+            JM1 = J - JSH
+            JP1 = J + JSH
+            JM2 = J - JST
+            JP2 = J + JST
+            JM3 = JM2 - JSH
+            JP3 = JP2 + JSH
+            IF (JST == 1) THEN
+               vecB(:M) = TWO*matY(:M,J)
+               matY(:M,J) = matY(:M,JM2) + matY(:M,JP2)
+            ELSE
+               DO I = 1, M
+                  T = matY(I,J) - matY(I,JM1) - matY(I,JP1) + matY(I,JM2) + &
+                      matY(I,JP2)
+                  vecB(I) = T + matY(I,J) - matY(I,JM3) - matY(I,JP3)
+                  matY(I,J) = T
+               END DO
+            ENDIF
+            CALL TRIX (JST, 0, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+            matY(:M,J) = matY(:M,J) + vecB(:M)
+         END DO
+      ENDIF
+!
+!     REDUCTION FOR LAST UNKNOWN
+!
+      SELECT CASE (NODD) 
+      CASE DEFAULT
+         GO TO (152,120) IRREG
+!
+!     ODD NUMBER OF UNKNOWNS
+!
+  120    CONTINUE
+         JSP = JSP + L
+         J = JSP
+         JM1 = J - JSH
+         JP1 = J + JSH
+         JM2 = J - JST
+         JP2 = J + JST
+         JM3 = JM2 - JSH
+         GO TO (123,121) ISTAG
+  121    CONTINUE
+         IF (JST /= 1) GO TO 123
+         vecB(:M) = matY(:M,J)
+         matY(:M,J) = ZERO
+         GO TO 130
+  123    CONTINUE
+         SELECT CASE (NODDPR) 
+         CASE DEFAULT
+            vecB(:M) = HALF*(matY(:M,JM2)-matY(:M,JM1)-matY(:M,JM3)) + &
+                       vecP(IP+1:M+IP)+ matY(:M,J)
+         CASE (2) 
+            vecB(:M) = HALF*(matY(:M,JM2)-matY(:M,JM1)-matY(:M,JM3)) +  &
+                       matY(:M,JP2) - matY(:M,JP1) + matY(:M,J)
+         END SELECT
+         matY(:M,J) = HALF*(matY(:M,J)-matY(:M,JM1)-matY(:M,JP1))
+  130    CONTINUE
+         CALL TRIX (JST, 0, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+         IP = IP + M
+         IPSTOR = MAX0(IPSTOR,IP + M)
+         vecP(IP+1:M+IP) = matY(:M,J) + vecB(:M)
+         vecB(:M) = matY(:M,JP2) + vecP(IP+1:M+IP)
+         IF (LR == 0) THEN
+            DO I = 1, JST
+               KRPI = KR + I
+               TCOS(KRPI) = TCOS(I)
+            END DO
+         ELSE
+            CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(JST+1), 4*NR-JST)
+            CALL MERGE (TCOS, 0, JST, JST, LR, KR, 4*NR)
+         ENDIF
+         CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR)
+         CALL TRIX (KR, KR, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+         matY(:M,J) = matY(:M,JM2) + vecB(:M) + vecP(IP+1:M+IP)
+         LR = KR
+         KR = KR + L
+!
+!     EVEN NUMBER OF UNKNOWNS
+!
+      CASE (2) 
+         JSP = JSP + L
+         J = JSP
+         JM1 = J - JSH
+         JP1 = J + JSH
+         JM2 = J - JST
+         JP2 = J + JST
+         JM3 = JM2 - JSH
+         SELECT CASE (IRREG) 
+         CASE DEFAULT
+            JSTSAV = JST
+            IDEG = JST
+            KR = L
+         CASE (2) 
+            CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR)
+            CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(KR+1), 4*NR-KR)
+            IDEG = KR
+            KR = KR + JST
+         END SELECT
+         IF (JST == 1) THEN
+            IRREG = 2
+            vecB(:M) = matY(:M,J)
+            matY(:M,J) = matY(:M,JM2)
+         ELSE
+            vecB(:M) = matY(:M,J) + HALF*(matY(:M,JM2)-matY(:M,JM1)- &
+                       matY(:M,JM3))
+            SELECT CASE (IRREG) 
+            CASE DEFAULT
+               matY(:M,J) = matY(:M,JM2) + HALF*(matY(:M,J)- &
+                            matY(:M,JM1)-matY(:M,JP1))
+               IRREG = 2
+            CASE (2) 
+               SELECT CASE (NODDPR) 
+               CASE DEFAULT
+                  matY(:M,J) = matY(:M,JM2) + vecP(IP+1:M+IP)
+                  IP = IP - M
+               CASE (2) 
+                  matY(:M,J) = matY(:M,JM2) + matY(:M,J) - matY(:M,JM1)
+               END SELECT
+            END SELECT
+         ENDIF
+         CALL TRIX (IDEG, LR, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+         matY(:M,J) = matY(:M,J) + vecB(:M)
+      END SELECT
+  152 CONTINUE
+      NUN = NUN/2
+      NODDPR = NODD
+      JSH = JST
+      JST = 2*JST
+      IF (NUN >= 2) GO TO 108
+!
+!     START SOLUTION.
+!
+      J = JSP
+      vecB(:M) = matY(:M,J)
+      SELECT CASE (IRREG) 
+      CASE DEFAULT
+         CALL COSGEN (JST, 1, HALF, ZERO, TCOS, 4*NR)
+         IDEG = JST
+      CASE (2) 
+         KR = LR + JST
+         CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR)
+         CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(KR+1), 4*NR-KR)
+         IDEG = KR
+      END SELECT
+      CALL TRIX (IDEG, LR, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+      JM1 = J - JSH
+      JP1 = J + JSH
+      SELECT CASE (IRREG) 
+      CASE DEFAULT
+         matY(:M,J) = HALF*(matY(:M,J)-matY(:M,JM1)-matY(:M,JP1)) + vecB(:M)
+      CASE (2) 
+         SELECT CASE (NODDPR) 
+         CASE DEFAULT
+            matY(:M,J) = vecP(IP+1:M+IP) + vecB(:M)
+            IP = IP - M
+         CASE (2) 
+            matY(:M,J) = matY(:M,J) - matY(:M,JM1) + vecB(:M)
+         END SELECT
+      END SELECT
+  164 CONTINUE
+      JST = JST/2
+      JSH = JST/2
+      NUN = 2*NUN
+      IF (NUN > N) GO TO 183
+      DO J = JST, N, L
+         JM1 = J - JSH
+         JP1 = J + JSH
+         JM2 = J - JST
+         JP2 = J + JST
+         IF (J <= JST) THEN
+            vecB(:M) = matY(:M,J) + matY(:M,JP2)
+         ELSE
+            IF (JP2 <= N) GO TO 168
+            vecB(:M) = matY(:M,J) + matY(:M,JM2)
+            IF (JST < JSTSAV) IRREG = 1
+            GO TO (170,171) IRREG
+  168       CONTINUE
+            vecB(:M) = matY(:M,J) + matY(:M,JM2) + matY(:M,JP2)
+         ENDIF
+  170    CONTINUE
+         CALL COSGEN (JST, 1, HALF, ZERO, TCOS, 4*NR)
+         IDEG = JST
+         JDEG = 0
+         GO TO 172
+  171    CONTINUE
+         IF (J + L > N) LR = LR - JST
+         KR = JST + LR
+         CALL COSGEN (KR, JSTSAV, ZERO, FI, TCOS, 4*NR)
+         CALL COSGEN (LR, JSTSAV, ZERO, FI, TCOS(KR+1), 4*NR-KR)
+         IDEG = KR
+         JDEG = LR
+  172    CONTINUE
+         CALL TRIX (IDEG, JDEG, M, vecBA, vecBB, vecBC, vecB, TCOS, 4*NR, vecD, vecW)
+         IF (JST <= 1) THEN
+            matY(:M,J) = vecB(:M)
+         ELSE
+            IF (JP2 > N) GO TO 177
+  175       CONTINUE
+            matY(:M,J) = HALF*(matY(:M,J)-matY(:M,JM1)-matY(:M,JP1)) + vecB(:M)
+            CYCLE 
+  177       CONTINUE
+            GO TO (175,178) IRREG
+  178       CONTINUE
+            IF (J + JSH <= N) THEN
+               matY(:M,J) = vecB(:M) + vecP(IP+1:M+IP)
+               IP = IP - M
+            ELSE
+               matY(:M,J) = vecB(:M) + matY(:M,J) - matY(:M,JM1)
+            ENDIF
+         ENDIF
+      END DO
+      L = L/2
+      GO TO 164
+  183 CONTINUE
+      vecW(1) = IPSTOR
+      RETURN 
+      END SUBROUTINE POISD2
+
+
+      SUBROUTINE POISN2(M, N, ISTAG, MIXBND, vecA, vecBB, vecC, &
+                        matQ, IDIMQ, vecB, vecB2,vecB3, &
+                        vecW, vecW2, vecW3, vecD, TCOS, vecP, IDIMP)
+      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, N, ISTAG, MIXBND, IDIMQ,IDIMP
+      DOUBLE PRECISION, DIMENSION(M), INTENT(IN)   :: vecA, vecBB, vecC
+      DOUBLE PRECISION, DIMENSION(IDIMQ,N), INTENT(INOUT) :: matQ
+      DOUBLE PRECISION, DIMENSION(M),INTENT(INOUT) :: vecB, vecB2, vecB3,  &
+                                                      vecD, vecW, vecW2, vecW3
+      DOUBLE PRECISION, DIMENSION(IDIMP),INTENT(INOUT) :: vecP
+      DOUBLE PRECISION, DIMENSION(4*N),INTENT(INOUT)   :: TCOS
+!-----------------------------------------------
+!   L o c a l   V a r i a b l e s
+!-----------------------------------------------
+      INTEGER, DIMENSION(4) :: K
+      INTEGER :: K1, K2, K3, K4, MR, IP, IPSTOR, I2R, JR, NR, NLAST, &
+                 KR, LR, I, NROD, JSTART, JSTOP, I2RBY2, &
+                 J, JP1, JP2, JP3, JM1,JM2, JM3, NRODPR, II, I1, I2, &
+                 JR2, NLASTP, JSTEP
+      DOUBLE PRECISION :: FISTAG, FNUM, FDEN, FI, T
+!-----------------------------------------------
+!
+!     SUBROUTINE TO SOLVE POISSON'S EQUATION WITH NEUMANN BOUNDARY
+!     CONDITIONS.
+!
+!     ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS A.
+!     ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS A-I.
+!     MIXBND = 1 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTH BOUNDARIES.
+!     MIXBND = 2 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTTOM AND
+!     DIRICHLET CONDITION AT TOP.  (FOR THIS CASE, MUST HAVE ISTAG = 1.)
+!
+      EQUIVALENCE (K(1), K1), (K(2), K2), (K(3), K3), (K(4), K4)
+      FISTAG = 3 - ISTAG
+      FNUM = ONE/FLOAT(ISTAG)
+      FDEN = HALF*FLOAT(ISTAG - 1)
+      MR = M
+      IP = -MR
+      IPSTOR = 0
+      I2R = 1
+      JR = 2
+      NR = N
+      NLAST = N
+      KR = 1
+      LR = 0
+      GO TO (101,103) ISTAG
+  101 CONTINUE
+      matQ(:MR,N) = HALF*matQ(:MR,N)
+      GO TO (103,104) MIXBND
+  103 CONTINUE
+      IF (N <= 3) GO TO 155
+  104 CONTINUE
+      JR = 2*I2R
+      NROD = 1
+      IF ((NR/2)*2 == NR) NROD = 0
+      SELECT CASE (MIXBND) 
+      CASE DEFAULT
+         JSTART = 1
+      CASE (2) 
+         JSTART = JR
+         NROD = 1 - NROD
+      END SELECT
+      JSTOP = NLAST - JR
+      IF (NROD == 0) JSTOP = JSTOP - I2R
+      CALL COSGEN (I2R, 1, HALF, ZERO, TCOS, 4*NR)
+      I2RBY2 = I2R/2
+      IF (JSTOP < JSTART) THEN
+         J = JR
+      ELSE
+         DO J = JSTART, JSTOP, JR
+            JP1 = J + I2RBY2
+            JP2 = J + I2R
+            JP3 = JP2 + I2RBY2
+            JM1 = J - I2RBY2
+            JM2 = J - I2R
+            JM3 = JM2 - I2RBY2
+            IF (J == 1) THEN
+               JM1 = JP1
+               JM2 = JP2
+               JM3 = JP3
+            ENDIF
+            IF (I2R == 1) THEN
+               IF (J == 1) JM2 = JP2
+               vecB(:MR) = TWO*matQ(:MR,J)
+               matQ(:MR,J) = matQ(:MR,JM2) + matQ(:MR,JP2)
+            ELSE
+               DO I = 1, MR
+                  FI = matQ(I,J)
+                  matQ(I,J)=matQ(I,J)-matQ(I,JM1)-matQ(I,JP1)+ &
+                            matQ(I,JM2)+matQ(I,JP2)
+                  vecB(I) = FI + matQ(I,J) - matQ(I,JM3) - matQ(I,JP3)
+               END DO
+            ENDIF
+            CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+            matQ(:MR,J) = matQ(:MR,J) + vecB(:MR)
+!
+!     END OF REDUCTION FOR REGULAR UNKNOWNS.
+!
+         END DO
+!
+!     BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN.
+!
+         J = JSTOP + JR
+      ENDIF
+      NLAST = J
+      JM1 = J - I2RBY2
+      JM2 = J - I2R
+      JM3 = JM2 - I2RBY2
+      IF (NROD /= 0) THEN
+!
+!     ODD NUMBER OF UNKNOWNS
+!
+         IF (I2R == 1) THEN
+            vecB(:MR) = FISTAG*matQ(:MR,J)
+            matQ(:MR,J) = matQ(:MR,JM2)
+         ELSE
+            vecB(:MR) = matQ(:MR,J) + HALF*(matQ(:MR,JM2)- &
+                        matQ(:MR,JM1)-matQ(:MR,JM3))
+            IF (NRODPR == 0) THEN
+               matQ(:MR,J) = matQ(:MR,JM2) + vecP(IP+1:MR+IP)
+               IP = IP - MR
+            ELSE
+               matQ(:MR,J) = matQ(:MR,J) - matQ(:MR,JM1) + matQ(:MR,JM2)
+            ENDIF
+            IF (LR /= 0) THEN
+               CALL COSGEN (LR, 1, HALF, FDEN, TCOS(KR+1), 4*NR-KR)
+            ELSE
+               vecB(:MR) = FISTAG*vecB(:MR)
+            ENDIF
+         ENDIF
+         CALL COSGEN (KR, 1, HALF, FDEN, TCOS, 4*NR)
+         CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,J) = matQ(:MR,J) + vecB(:MR)
+         KR = KR + I2R
+      ELSE
+         JP1 = J + I2RBY2
+         JP2 = J + I2R
+         IF (I2R == 1) THEN
+            vecB(:MR) = matQ(:MR,J)
+            CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+            IP = 0
+            IPSTOR = MR
+            SELECT CASE (ISTAG) 
+            CASE DEFAULT
+               vecP(:MR) = vecB(:MR)
+               vecB(:MR) = vecB(:MR) + matQ(:MR,N)
+               TCOS(1) = ONE
+               TCOS(2) = ZERO
+               CALL TRIX (1, 1, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+               matQ(:MR,J) = matQ(:MR,JM2) + vecP(:MR) + vecB(:MR)
+               GO TO 150
+            CASE (1) 
+               vecP(:MR) = vecB(:MR)
+               matQ(:MR,J) = matQ(:MR,JM2) + TWO*matQ(:MR,JP2) + 3.*vecB(:MR)
+               GO TO 150
+            END SELECT
+         ENDIF
+         vecB(:MR) = matQ(:MR,J) + HALF*(matQ(:MR,JM2)- &
+                     matQ(:MR,JM1)-matQ(:MR,JM3))
+         IF (NRODPR == 0) THEN
+            vecB(:MR) = vecB(:MR) + vecP(IP+1:MR+IP)
+         ELSE
+            vecB(:MR) = vecB(:MR) + matQ(:MR,JP2) - matQ(:MR,JP1)
+         ENDIF
+         CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         IP = IP + MR
+         IPSTOR = MAX0(IPSTOR,IP + MR)
+         vecP(IP+1:MR+IP) = vecB(:MR) + HALF*(matQ(:MR,J)- &
+                            matQ(:MR,JM1)-matQ(:MR,JP1))
+         vecB(:MR) = vecP(IP+1:MR+IP) + matQ(:MR,JP2)
+         IF (LR /= 0) THEN
+            CALL COSGEN (LR, 1, HALF, FDEN, TCOS(I2R+1), 4*NR-I2R)
+            CALL MERGE (TCOS, 0, I2R, I2R, LR, KR, 4*NR)
+         ELSE
+            DO I = 1, I2R
+               II = KR + I
+               TCOS(II) = TCOS(I)
+            END DO
+         ENDIF
+         CALL COSGEN (KR, 1, HALF, FDEN, TCOS, 4*NR)
+         IF (LR == 0) THEN
+            GO TO (146,145) ISTAG
+         ENDIF
+  145    CONTINUE
+         CALL TRIX (KR, KR, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         GO TO 148
+  146    CONTINUE
+         vecB(:MR) = FISTAG*vecB(:MR)
+  148    CONTINUE
+         matQ(:MR,J) = matQ(:MR,JM2) + vecP(IP+1:MR+IP) + vecB(:MR)
+  150    CONTINUE
+         LR = KR
+         KR = KR + JR
+      ENDIF
+      SELECT CASE (MIXBND) 
+      CASE DEFAULT
+         NR = (NLAST - 1)/JR + 1
+         IF (NR <= 3) GO TO 155
+      CASE (2) 
+         NR = NLAST/JR
+         IF (NR <= 1) GO TO 192
+      END SELECT
+      I2R = JR
+      NRODPR = NROD
+      GO TO 104
+  155 CONTINUE
+      J = 1 + JR
+      JM1 = J - I2R
+      JP1 = J + I2R
+      JM2 = NLAST - I2R
+      IF (NR /= 2) THEN
+         IF (LR /= 0) GO TO 170
+         IF (N == 3) THEN
+!
+!     CASE N = 3.
+!
+            GO TO (156,168) ISTAG
+  156       CONTINUE
+            vecB(:MR) = matQ(:MR,2)
+            TCOS(1) = ZERO
+            CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+            matQ(:MR,2) = vecB(:MR)
+            vecB(:MR) = 4.*vecB(:MR) + matQ(:MR,1) + TWO*matQ(:MR,3)
+            TCOS(1) = -TWO
+            TCOS(2) = TWO
+            I1 = 2
+            I2 = 0
+            CALL TRIX (I1, I2, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+            matQ(:MR,2) = matQ(:MR,2) + vecB(:MR)
+            vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,2)
+            TCOS(1) = ZERO
+            CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+            matQ(:MR,1) = vecB(:MR)
+            JR = 1
+            I2R = 0
+            GO TO 194
+         ENDIF
+!
+!     CASE N = 2**P+1
+!
+         GO TO (162,170) ISTAG
+  162    CONTINUE
+         vecB(:MR) = matQ(:MR,J) + HALF*matQ(:MR,1) - &
+                     matQ(:MR,JM1) + matQ(:MR,NLAST) - matQ(:MR,JM2)
+         CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR)
+         CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,J) = HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1)) + vecB(:MR)
+         vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,NLAST) + 4.*matQ(:MR,J)
+         JR2 = 2*JR
+         CALL COSGEN (JR, 1, ZERO, ZERO, TCOS, 4*NR)
+         TCOS(JR+1:JR*2) = -TCOS(JR:1:(-1))
+         CALL TRIX (JR2, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,J) = matQ(:MR,J) + vecB(:MR)
+         vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,J)
+         CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR)
+         CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,1) = HALF*matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR)
+         GO TO 194
+!
+!     CASE OF GENERAL N WITH NR = 3 .
+!
+  168    CONTINUE
+         vecB(:MR) = matQ(:MR,2)
+         matQ(:MR,2) = ZERO
+         vecB2(:MR) = matQ(:MR,3)
+         vecB3(:MR) = matQ(:MR,1)
+         JR = 1
+         I2R = 0
+         J = 2
+         GO TO 177
+  170    CONTINUE
+         vecB(:MR) = HALF*matQ(:MR,1) - matQ(:MR,JM1) + matQ(:MR,J)
+         IF (NROD == 0) THEN
+            vecB(:MR) = vecB(:MR) + vecP(IP+1:MR+IP)
+         ELSE
+            vecB(:MR) = vecB(:MR) + matQ(:MR,NLAST) - matQ(:MR,JM2)
+         ENDIF
+         DO I = 1, MR
+            T = HALF*(matQ(I,J)-matQ(I,JM1)-matQ(I,JP1))
+            matQ(I,J) = T
+            vecB2(I) = matQ(I,NLAST) + T
+            vecB3(I) = matQ(I,1) + TWO*T
+         END DO
+  177    CONTINUE
+         K1 = KR + 2*JR - 1
+         K2 = KR + JR
+         TCOS(K1+1) = -TWO
+         K4 = K1 + 3 - ISTAG
+         CALL COSGEN (K2 + ISTAG - 2, 1, ZERO, FNUM, TCOS(K4), 4*NR-K4+1)
+         K4 = K1 + K2 + 1
+         CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K4), 4*NR-K4+1)
+         CALL MERGE (TCOS, K1, K2, K1 + K2, JR - 1, 0, 4*NR)
+         K3 = K1 + K2 + LR
+         CALL COSGEN (JR, 1, HALF, ZERO, TCOS(K3+1), 4*NR-K3)
+         K4 = K3 + JR + 1
+         CALL COSGEN (KR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1)
+         CALL MERGE (TCOS, K3, JR, K3 + JR, KR, K1, 4*NR)
+         IF (LR /= 0) THEN
+            CALL COSGEN (LR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1)
+            CALL MERGE (TCOS, K3, JR, K3 + JR, LR, K3 - LR, 4*NR)
+            CALL COSGEN (KR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1)
+         ENDIF
+         K3 = KR
+         K4 = KR
+         CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, &
+                    TCOS, 4*NR, vecD, vecW, vecW2, vecW3)
+         vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR)
+         TCOS(1) = TWO
+         CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,J) = matQ(:MR,J) + vecB(:MR)
+         vecB(:MR) = matQ(:MR,1) + TWO*matQ(:MR,J)
+         CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR)
+         CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         IF (JR == 1) THEN
+            matQ(:MR,1) = vecB(:MR)
+            GO TO 194
+         ENDIF
+         matQ(:MR,1) = HALF*matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR)
+         GO TO 194
+      ENDIF
+      IF (N == 2) THEN
+!
+!     CASE  N = 2
+!
+         vecB(:MR) = matQ(:MR,1)
+         TCOS(1) = ZERO
+         CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,1) = vecB(:MR)
+         vecB(:MR) = TWO*(matQ(:MR,2)+vecB(:MR))*FISTAG
+         TCOS(1) = -FISTAG
+         TCOS(2) = TWO
+         CALL TRIX (2, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,1) = matQ(:MR,1) + vecB(:MR)
+         JR = 1
+         I2R = 0
+         GO TO 194
+      ENDIF
+      vecB3(:MR) = ZERO
+      vecB(:MR) = matQ(:MR,1) + TWO*vecP(IP+1:MR+IP)
+      matQ(:MR,1) = HALF*matQ(:MR,1) - matQ(:MR,JM1)
+      vecB2(:MR) = TWO*(matQ(:MR,1)+matQ(:MR,NLAST))
+      K1 = KR + JR - 1
+      TCOS(K1+1) = -TWO
+      K4 = K1 + 3 - ISTAG
+      CALL COSGEN (KR + ISTAG - 2, 1, ZERO, FNUM, TCOS(K4), 4*NR-K4+1)
+      K4 = K1 + KR + 1
+      CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K4), 4*NR-K4+1)
+      CALL MERGE (TCOS, K1, KR, K1 + KR, JR - 1, 0, 4*NR)
+      CALL COSGEN (KR, 1, HALF, FDEN, TCOS(K1+1), 4*NR-K1)
+      K2 = KR
+      K4 = K1 + K2 + 1
+      CALL COSGEN (LR, 1, HALF, FDEN, TCOS(K4), 4*NR-K4+1)
+      K3 = LR
+      K4 = 0
+      CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, TCOS, &
+                 4*NR, vecD, vecW, vecW2, vecW3)
+      vecB(:MR) = vecB(:MR) + vecB2(:MR)
+      TCOS(1) = TWO
+      CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+      matQ(:MR,1) = matQ(:MR,1) + vecB(:MR)
+      GO TO 194
+  192 CONTINUE
+      vecB(:MR) = matQ(:MR,NLAST)
+      GO TO 196
+  194 CONTINUE
+      J = NLAST - JR
+      vecB(:MR) = matQ(:MR,NLAST) + matQ(:MR,J)
+  196 CONTINUE
+      JM2 = NLAST - I2R
+      IF (JR == 1) THEN
+         matQ(:MR,NLAST) = ZERO
+      ELSE
+         IF (NROD == 0) THEN
+            matQ(:MR,NLAST) = vecP(IP+1:MR+IP)
+            IP = IP - MR
+         ELSE
+            matQ(:MR,NLAST) = matQ(:MR,NLAST) - matQ(:MR,JM2)
+         ENDIF
+      ENDIF
+      CALL COSGEN (KR, 1, HALF, FDEN, TCOS, 4*NR)
+      CALL COSGEN (LR, 1, HALF, FDEN, TCOS(KR+1), 4*NR-KR)
+      IF (LR == 0) THEN
+         vecB(:MR) = FISTAG*vecB(:MR)
+      ENDIF
+      CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+      matQ(:MR,NLAST) = matQ(:MR,NLAST) + vecB(:MR)
+      NLASTP = NLAST
+  206 CONTINUE
+      JSTEP = JR
+      JR = I2R
+      I2R = I2R/2
+      IF (JR == 0) GO TO 222
+      SELECT CASE (MIXBND) 
+      CASE DEFAULT
+         JSTART = 1 + JR
+      CASE (2) 
+         JSTART = JR
+      END SELECT
+      KR = KR - JR
+      IF (NLAST + JR <= N) THEN
+         KR = KR - JR
+         NLAST = NLAST + JR
+         JSTOP = NLAST - JSTEP
+      ELSE
+         JSTOP = NLAST - JR
+      ENDIF
+      LR = KR - JR
+      CALL COSGEN (JR, 1, HALF, ZERO, TCOS, 4*NR)
+      DO J = JSTART, JSTOP, JSTEP
+         JM2 = J - JR
+         JP2 = J + JR
+         IF (J == JR) THEN
+            vecB(:MR) = matQ(:MR,J) + matQ(:MR,JP2)
+         ELSE
+            vecB(:MR) = matQ(:MR,J) + matQ(:MR,JM2) + matQ(:MR,JP2)
+         ENDIF
+         IF (JR == 1) THEN
+            matQ(:MR,J) = ZERO
+         ELSE
+            JM1 = J - I2R
+            JP1 = J + I2R
+            matQ(:MR,J) = HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1))
+         ENDIF
+         CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS, 4*NR, vecD, vecW)
+         matQ(:MR,J) = matQ(:MR,J) + vecB(:MR)
+      END DO
+      NROD = 1
+      IF (NLAST + I2R <= N) NROD = 0
+      IF (NLASTP /= NLAST) GO TO 194
+      GO TO 206
+  222 CONTINUE
+      vecW(1) = IPSTOR
+      RETURN 
+      END SUBROUTINE POISN2
+
+
+      SUBROUTINE POISP2(M, N, vecA, vecBB, vecC, matQ, IDIMQ, &
+                        vecB, vecB2, vecB3, vecW, vecW2, vecW3, &
+                        vecD, TCOS, vecP, IP)
+      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, N, IDIMQ, IP
+      DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecA, vecBB, vecC
+      DOUBLE PRECISION, DIMENSION(4*N), INTENT(INOUT) :: TCOS
+      DOUBLE PRECISION, DIMENSION(IDIMQ,N), INTENT(INOUT) :: matQ
+      DOUBLE PRECISION, DIMENSION(IP),INTENT(INOUT) :: vecP
+      DOUBLE PRECISION, DIMENSION(M),INTENT(INOUT) :: vecB, vecB2, vecB3, &
+                                                      vecD, vecW, vecW2, vecW3
+!-----------------------------------------------
+!   L o c a l   V a r i a b l e s
+!-----------------------------------------------
+      INTEGER          :: MR, NR, NRM1, J, NRMJ, NRPJ, I, IPSTOR, LH
+      DOUBLE PRECISION :: S, T
+!-----------------------------------------------
+!
+!     SUBROUTINE TO SOLVE POISSON EQUATION WITH PERIODIC BOUNDARY
+!     CONDITIONS.
+!
+      MR = M
+      NR = (N + 1)/2
+      NRM1 = NR - 1
+      IF (2*NR == N) THEN
+!
+!     EVEN NUMBER OF UNKNOWNS
+!
+         DO J = 1, NRM1
+            NRMJ = NR - J
+            NRPJ = NR + J
+            DO I = 1, MR
+               S = matQ(I,NRMJ) - matQ(I,NRPJ)
+               T = matQ(I,NRMJ) + matQ(I,NRPJ)
+               matQ(I,NRMJ) = S
+               matQ(I,NRPJ) = T
+            END DO
+         END DO
+         matQ(:MR,NR) = TWO*matQ(:MR,NR)
+         matQ(:MR,N)  = TWO*matQ(:MR,N)
+         CALL POISD2 (MR, NRM1, 1, vecA, vecBB, vecC, matQ, IDIMQ, &
+                      vecB, vecW, vecD, TCOS, vecP, IP)
+         IPSTOR = vecW(1)
+         CALL POISN2 (MR, NR + 1, 1, 1, vecA, vecBB, vecC, matQ(:,NR:),  &
+                      IDIMQ, vecB, vecB2, vecB3, vecW, vecW2, vecW3, vecD, &
+                      TCOS, vecP, IP)
+         IPSTOR = MAX0(IPSTOR,INT(vecW(1)))
+         DO J = 1, NRM1
+            NRMJ = NR - J
+            NRPJ = NR + J
+            DO I = 1, MR
+               S = HALF*(matQ(I,NRPJ)+matQ(I,NRMJ))
+               T = HALF*(matQ(I,NRPJ)-matQ(I,NRMJ))
+               matQ(I,NRMJ) = S
+               matQ(I,NRPJ) = T
+            END DO
+         END DO
+         matQ(:MR,NR) = HALF*matQ(:MR,NR)
+         matQ(:MR,N)  = HALF*matQ(:MR,N)
+      ELSE
+         DO J = 1, NRM1
+            NRPJ = N + 1 - J
+            DO I = 1, MR
+               S = matQ(I,J) - matQ(I,NRPJ)
+               T = matQ(I,J) + matQ(I,NRPJ)
+               matQ(I,J) = S
+               matQ(I,NRPJ) = T
+            END DO
+         END DO
+         matQ(:MR,NR) = TWO*matQ(:MR,NR)
+         LH = NRM1/2
+         DO J = 1, LH
+            NRMJ = NR - J
+            DO I = 1, MR
+               S = matQ(I,J)
+               matQ(I,J) = matQ(I,NRMJ)
+               matQ(I,NRMJ) = S
+            END DO
+         END DO
+         CALL POISD2 (MR, NRM1, 2, vecA, vecBB, vecC, matQ, IDIMQ,  &
+                      vecB, vecW, vecD, TCOS, vecP, IP)
+         IPSTOR = vecW(1)
+         CALL POISN2 (MR, NR, 2, 1, vecA, vecBB, vecC, matQ(:,NR:), IDIMQ, &
+                      vecB, vecB2, vecB3, vecW, vecW2, vecW3, vecD, TCOS, &
+                      vecP, IP)
+         IPSTOR = MAX0(IPSTOR,INT(vecW(1)))
+         DO J = 1, NRM1
+            NRPJ = NR + J
+            DO I = 1, MR
+               S = HALF*(matQ(I,NRPJ)+matQ(I,J))
+               T = HALF*(matQ(I,NRPJ)-matQ(I,J))
+               matQ(I,NRPJ) = T
+               matQ(I,J) = S
+            END DO
+         END DO
+         matQ(:MR,NR) = HALF*matQ(:MR,NR)
+         DO J = 1, LH
+            NRMJ = NR - J
+            DO I = 1, MR
+               S = matQ(I,J)
+               matQ(I,J)    = matQ(I,NRMJ)
+               matQ(I,NRMJ) = S
+            END DO
+         END DO
+      ENDIF
+      vecW(1) = IPSTOR
+!
+! REVISION HISTORY---
+!
+! SEPTEMBER 1973    VERSION 1
+! APRIL     1976    VERSION 2
+! JANUARY   1978    VERSION 3
+! DECEMBER  1979    VERSION 3.1
+! FEBRUARY  1985    DOCUMENTATION UPGRADE
+! NOVEMBER  1988    VERSION 3.2, FORTRAN 77 CHANGES
+! June      2004    Version 5.0, Fortran 90 changes
+!-----------------------------------------------------------------------
+      END SUBROUTINE POISP2