SUBROUTINE POISTGG(NPEROD,N,MPEROD,M,vecA,vecB,vecC,IDIMY,matY, &
                         IERROR,W,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(OUT) :: IERROR
      DOUBLE PRECISION, DIMENSION(M), INTENT(IN)         :: vecA,vecB,vecC
      DOUBLE PRECISION, DIMENSION(IDIMY,N),INTENT(INOUT) :: matY
      DOUBLE PRECISION, DIMENSION(IW)                    :: W
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IWBA, IWBB, IWBC, IWB2, IWB3, IWW1, IWW2, IWW3, &
                 IWD,IWTCOS, IWP, I, K, J, NP, MP, IPSTOR, &
                 IREV, MH, MHM1, MODD, NBY2, MSKIP
      DOUBLE PRECISION :: A1
!-----------------------------------------------
      IERROR = 0

      IWBA = M + 1
      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
      DO I = 1, M
         K = IWBA + I - 1
         ! write(*,*) '0: i,k:',i,k
         ! write(*,*) '0: vecA(I):',vecA(I)
         W(K) = -vecA(I)
         ! write(*,*) 'A: i,k,w(k):',i,k,w(k)
         K = IWBC + I - 1
         W(K) = -vecC(I)
         ! write(*,*) 'B: i,k,w(k):',i,k,w(k)
         K = IWBB + I - 1
         W(K) = TWO - vecB(I)
         ! write(*,*) 'C: i,k,w(k):',i,k,w(k)
         matY(I,:N) = -matY(I,:N)
         ! write(*,*) 'D: Y(i,:N):',Y(I,:N)
      END DO
      NP = NPEROD
      MP = MPEROD + 1
      GO TO (110,107) MP
  107 CONTINUE
      GO TO (108,108,108,119) NPEROD
  108 CONTINUE
      CALL POSTG2 (NP, N, M, W(IWBA:IWBA+M-1), W(IWBB:IWBB+M-1), &
                   W(IWBC:IWBC+M-1), IDIMY, matY, W(1:M), W(IWB2:IWB2+M-1),  &
                   W(IWB3:IWB3+M-1), W(IWW1:IWW1+M-1), W(IWW2:IWW2+M-1), &
                   W(IWW3:IWW3+M-1), W(IWD:IWD+M-1), W(IWTCOS:IWTCOS+4*N-1), &
                   W(IWP:),IW-IWP)
      IPSTOR = W(IWW1)
      IREV = 2
      IF (NPEROD == 4) GO TO 120
  109 CONTINUE
      GO TO (123,129) MP
  110 CONTINUE
      MH = (M + 1)/2
      MHM1 = MH - 1
      MODD = 1
      IF (MH*2 == M) MODD = 2
      DO J = 1, N
         DO I = 1, MHM1
            W(I) = matY(MH-I,J) - matY(I+MH,J)
            W(I+MH) = matY(MH-I,J) + matY(I+MH,J)
         END DO
         W(MH) = TWO*matY(MH,J)
         GO TO (113,112) MODD
  112    CONTINUE
         W(M) = TWO*matY(M,J)
  113    CONTINUE
         matY(:M,J) = W(:M)
      END DO
      K = IWBC + MHM1 - 1
      I = IWBA + MHM1
      W(K) = 0.
      W(I) = 0.
      W(K+1) = TWO*W(K+1)
      SELECT CASE (MODD) 
      CASE DEFAULT
         K = IWBB + MHM1 - 1
         W(K) = W(K) - W(I-1)
         W(IWBC-1) = W(IWBC-1) + W(IWBB-1)
      CASE (2) 
         W(IWBB-1) = W(K+1)
      END SELECT
      GO TO 107
  119 CONTINUE
      IREV = 1
      NBY2 = N/2
      NP = 2
  120 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 (108,109) IREV
  123 CONTINUE
      DO J = 1, N
         W(MH-1:MH-MHM1:(-1)) = HALF*(matY(MH+1:MHM1+MH,J)+matY(:MHM1,J))
         W(MH+1:MHM1+MH) = HALF*(matY(MH+1:MHM1+MH,J)-matY(:MHM1,J))
         W(MH) = HALF*matY(MH,J)
         GO TO (126,125) MODD
  125    CONTINUE
         W(M) = HALF*matY(M,J)
  126    CONTINUE
         matY(:M,J) = W(:M)
      END DO
  129 CONTINUE
      W(1) = IPSTOR + IWP - 1
!
      END SUBROUTINE POISTGG

      SUBROUTINE POSTG2(NPEROD, N, M, vecA, vecBB, vecC, IDIMQ, matQ, &
                        vecB, vecB2, vecB3, vecW,vecW2, vecW3, vecD,  &
                        TCOS, vecP,IvecP)
      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,M,IDIMQ,IvecP
      DOUBLE PRECISION, DIMENSION(M),INTENT(IN)        :: vecA,vecBB,vecC
      DOUBLE PRECISION, DIMENSION(M),INTENT(OUT)       :: vecB
      DOUBLE PRECISION, DIMENSION(IDIMQ,N),INTENT(INOUT) :: matQ
      DOUBLE PRECISION, DIMENSION(M),INTENT(INOUT)     :: vecB2,vecB3,vecD, &
                                                          vecW2,vecW3, vecW
      DOUBLE PRECISION, DIMENSION(IvecP),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, NP, MR, IP, IPSTOR, I2R, JR, NR, NLAST, &
                 KR, LR, NROD, JSTART, JSTOP, I2RBY2, &
                 J, IJUMP, JP1, JP2, JP3, JM1, JM2, JM3, I, NRODPR, II,  &
                 NLASTP, JSTEP
      DOUBLE PRECISION :: FNUM, FNUM2, FI, T
!-----------------------------------------------
!
!     SUBROUTINE TO SOLVE POISSON'S EQUATION ON A STAGGERED GRID.
!
      EQUIVALENCE (K(1), K1), (K(2), K2), (K(3), K3), (K(4), K4)
      NP = NPEROD
      FNUM = HALF*DBLE(NP/3)
      FNUM2 = HALF*DBLE(NP/2)
      MR = M
      IP = -MR
      IPSTOR = 0
      I2R = 1
      JR = 2
      NR = N
      NLAST = N
      KR = 1
      LR = 0
      IF (NR > 3) THEN
  101    CONTINUE
         JR = 2*I2R
         NROD = 1
         IF ((NR/2)*2 == NR) NROD = 0
         JSTART = 1
         JSTOP = NLAST - JR
         IF (NROD == 0) JSTOP = JSTOP - I2R
         I2RBY2 = I2R/2
         IF (JSTOP < JSTART) THEN
            J = JR
         ELSE
            IJUMP = 1
            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
                  CALL COSGEN (I2R, 1, FNUM, HALF, TCOS, size(TCOS) )
                  IF (I2R == 1) THEN
                     vecB(:MR)   = matQ(:MR,1)
                     matQ(:MR,1) = matQ(:MR,2)
                     GO TO 112
                  ENDIF
                  vecB(:MR)   = matQ(:MR,1) + &
                                HALF*(matQ(:MR,JP2)-matQ(:MR,JP1)-matQ(:MR,JP3))
                  matQ(:MR,1) = matQ(:MR,JP2) + matQ(:MR,1) - matQ(:MR,JP1)
                  GO TO 112
               ENDIF
               GO TO (107,108) IJUMP
  107          CONTINUE
               IJUMP = 2
               CALL COSGEN (I2R, 1, HALF, ZERO, TCOS, size(TCOS) )
  108          CONTINUE
               IF (I2R == 1) THEN
                  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
  112          CONTINUE
               CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS, &
                          size(TCOS), 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)   = 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) CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(KR+1), &
                                         size(tcos)-KR)
            ENDIF
            CALL COSGEN (KR, 1, FNUM2, HALF, TCOS, size(tcos))
            CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                       size(TCOS), 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)
               TCOS(1) = 0.
               CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
               IP = 0
               IPSTOR = MR
               vecP(:MR) = vecB(:MR)
               vecB(:MR) = vecB(:MR) + matQ(:MR,N)
               TCOS(1) = -ONE + TWO*DBLE(NP/2)
               TCOS(2) = 0.
               CALL TRIX (1, 1, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
               matQ(:MR,J) = matQ(:MR,JM2) + vecP(:MR) + vecB(:MR)
            ELSE
               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 COSGEN (I2R, 1, HALF, ZERO, TCOS, size(tcos))
               CALL TRIX (I2R, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), 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, FNUM2, HALF, TCOS(I2R+1), size(tcos)-I2R)
                  CALL MERGE (TCOS, 0, I2R, I2R, LR, KR, size(TCOS))
               ELSE
                  DO I = 1, I2R
                     II = KR + I
                     TCOS(II) = TCOS(I)
                  END DO
               ENDIF
               CALL COSGEN (KR, 1, FNUM2, HALF, TCOS, size(TCOS))
               CALL TRIX (KR, KR, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
               matQ(:MR,J) = matQ(:MR,JM2) + vecP(IP+1:MR+IP) + vecB(:MR)
            ENDIF
            LR = KR
            KR = KR + JR
         ENDIF
         NR = (NLAST - 1)/JR + 1
         IF (NR <= 3) GO TO 142
         I2R = JR
         NRODPR = NROD
         GO TO 101
      ENDIF ! IF (NR > 3)
  142 CONTINUE
      J = 1 + JR
      JM1 = J - I2R
      JP1 = J + I2R
      JM2 = NLAST - I2R
      IF (NR /= 2) THEN
         IF (LR == 0) THEN
            IF (N == 3) THEN
!
!     CASE N = 3.
!
               GO TO (143,148,143) NP
  143          CONTINUE
               vecB(:MR) = matQ(:MR,2)
               vecB2(:MR) = matQ(:MR,1) + matQ(:MR,3)
               vecB3(:MR) = 0.
               SELECT CASE (NP) 
               CASE DEFAULT
                  TCOS(1) = -ONE
                  TCOS(2) = ONE
                  K1 = 1
               CASE (1:2) 
                  TCOS(1) = -TWO
                  TCOS(2) = ONE
                  TCOS(3) = -ONE
                  K1 = 2
               END SELECT
               K2 = 1
               K3 = 0
               K4 = 0
               GO TO 150
  148          CONTINUE
               vecB(:MR) = matQ(:MR,2)
               vecB2(:MR) = matQ(:MR,3)
               vecB3(:MR) = matQ(:MR,1)
               CALL COSGEN (3, 1, HALF, ZERO, TCOS, size(TCOS))
               TCOS(4) = -ONE
               TCOS(5) = ONE
               TCOS(6) = -ONE
               TCOS(7) = ONE
               K1 = 3
               K2 = 2
               K3 = 1
               K4 = 1
  150          CONTINUE
               CALL TRI3(MR,vecA,vecBB,vecC,K,vecB,vecB2,vecB3,TCOS, &
                         size(TCOS),vecD,vecW,vecW2,vecW3)
               vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR)
               GO TO (153,153,152) NP
  152          CONTINUE
               TCOS(1) = TWO
               CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
  153          CONTINUE
               matQ(:MR,2) = vecB(:MR)
               vecB(:MR) = matQ(:MR,1) + vecB(:MR)
               TCOS(1) = -ONE + FOUR*FNUM
               CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
               matQ(:MR,1) = vecB(:MR)
               JR = 1
               I2R = 0
               GO TO 188
            ENDIF
!
!     CASE N = 2**P+1
!
            vecB(:MR)=matQ(:MR,J)+matQ(:MR,1)-matQ(:MR,JM1)+ &
                      matQ(:MR,NLAST)-matQ(:MR,JM2)
            GO TO (158,160,158) NP
  158       CONTINUE
            vecB2(:MR) = matQ(:MR,1) + matQ(:MR,NLAST) + &
                         matQ(:MR,J) - matQ(:MR,JM1) -matQ(:MR,JP1)
            vecB3(:MR) = 0.
            K1 = NLAST - 1
            K2 = NLAST + JR - 1
            CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(NLAST), size(TCOS)-NLAST+1)
            TCOS(K2) = TWO*DBLE(NP - 2)
            CALL COSGEN (JR, 1, HALF - FNUM, HALF, TCOS(K2+1),size(TCOS)-K2)
            K3 = (3 - NP)/2
            CALL MERGE (TCOS, K1, JR - K3, K2 - K3, JR + K3, 0, size(TCOS))
            K1 = K1 - 1 + K3
            CALL COSGEN (JR, 1, FNUM, HALF, TCOS(K1+1),size(TCOS)-K1)
            K2 = JR
            K3 = 0
            K4 = 0
            GO TO 162
  160       CONTINUE
            DO I = 1, MR
               FI = HALF*(matQ(I,J)-matQ(I,JM1)-matQ(I,JP1))
               vecB2(I) = matQ(I,1) + FI
               vecB3(I) = matQ(I,NLAST) + FI
            END DO
            K1 = NLAST + JR - 1
            K2 = K1 + JR - 1
            CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K1+1), size(TCOS)-K1)
            CALL COSGEN (NLAST, 1, HALF, ZERO, TCOS(K2+1), size(TCOS)-K2)
            CALL MERGE (TCOS, K1, JR - 1, K2, NLAST, 0, size(TCOS))
            K3 = K1 + NLAST - 1
            K4 = K3 + JR
            CALL COSGEN (JR, 1, HALF, HALF, TCOS(K3+1),SIZE(TCOS)-K3)
            CALL COSGEN (JR, 1, ZERO, HALF, TCOS(K4+1),SIZE(TCOS)-K4)
            CALL MERGE (TCOS, K3, JR, K4, JR, K1, size(TCOS) )
            K2 = NLAST - 1
            K3 = JR
            K4 = JR
  162       CONTINUE
            CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3,  &
                       TCOS, size(TCOS), vecD, vecW, vecW2, vecW3)
            vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR)
            IF (NP == 3) THEN
               TCOS(1) = TWO
               CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
            ENDIF
            matQ(:MR,J) = vecB(:MR) +  &
                          HALF*(matQ(:MR,J)-matQ(:MR,JM1)-matQ(:MR,JP1))
            vecB(:MR) = matQ(:MR,J) + matQ(:MR,1)
            CALL COSGEN (JR, 1, FNUM, HALF, TCOS,size(TCOS))
            CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
            matQ(:MR,1) = matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR)
            GO TO 188
         ENDIF
!
!     CASE OF GENERAL N WITH NR = 3 .
!
         vecB(:MR) = 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) + T
         END DO
         K1 = KR + 2*JR
         CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K1+1),size(TCOS)-K1)
         K2 = K1 + JR
         TCOS(K2) = TWO*DBLE(NP - 2)
         K4 = (NP - 1)*(3 - NP)
         K3 = K2 + 1 - K4
         CALL COSGEN(KR+JR+K4,1,HALF*DBLE(K4),ONE-DBLE(K4),TCOS(K3), &
                     size(TCOS)-K3+1 )
         K4 = 1 - NP/3
         CALL MERGE (TCOS, K1, JR - K4, K2 - K4, KR + JR + K4, 0, size(TCOS) )
         IF (NP == 3) K1 = K1 - 1
         K2 = KR + JR
         K4 = K1 + K2
         CALL COSGEN (KR, 1, FNUM2, HALF, TCOS(K4+1),size(TCOS)-K4)
         K3 = K4 + KR
         CALL COSGEN (JR, 1, FNUM, HALF, TCOS(K3+1),size(TCOS)-K3)
         CALL MERGE (TCOS, K4, KR, K3, JR, K1, size(TCOS) )
         K4 = K3 + JR
         CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(K4+1),size(TCOS)-K4)
         CALL MERGE (TCOS, K3, JR, K4, LR, K1 + K2, size(TCOS) )
         CALL COSGEN (KR, 1, FNUM2, HALF, TCOS(K3+1),size(TCOS)-K3)
         K3 = KR
         K4 = KR
         CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, &
                    TCOS, SIZE(TCOS), vecD, vecW, vecW2, vecW3)
         vecB(:MR) = vecB(:MR) + vecB2(:MR) + vecB3(:MR)
         IF (NP == 3) THEN
            TCOS(1) = TWO
            CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
         ENDIF
         matQ(:MR,J) = matQ(:MR,J) + vecB(:MR)
         vecB(:MR) = matQ(:MR,1) + matQ(:MR,J)
         CALL COSGEN (JR, 1, FNUM, HALF, TCOS,size(TCOS))
         CALL TRIX (JR, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                          size(TCOS), vecD, vecW)
         IF (JR == 1) THEN
            matQ(:MR,1) = vecB(:MR)
            GO TO 188
         ENDIF
         matQ(:MR,1) = matQ(:MR,1) - matQ(:MR,JM1) + vecB(:MR)
         GO TO 188
      ENDIF
      vecB3(:MR) = 0.
      vecB (:MR) = matQ(:MR,1) + vecP(IP+1:MR+IP)
      matQ (:MR,1) = matQ(:MR,1) - matQ(:MR,JM1)
      vecB2(:MR) = matQ(:MR,1) + matQ(:MR,NLAST)
      K1 = KR + JR
      K2 = K1 + JR
      CALL COSGEN (JR - 1, 1, ZERO, ONE, TCOS(K1+1),SIZE(TCOS)-K1)
      GO TO (182,183,182) NP
  182 CONTINUE
      TCOS(K2) = TWO*DBLE(NP - 2)
      CALL COSGEN (KR, 1, ZERO, ONE, TCOS(K2+1),SIZE(TCOS)-K2)
      GO TO 184
  183 CONTINUE
      CALL COSGEN (KR + 1, 1, HALF, ZERO, TCOS(K2),SIZE(TCOS)-K2)
  184 CONTINUE
      K4 = 1 - NP/3
      CALL MERGE (TCOS, K1, JR - K4, K2 - K4, KR + K4, 0, size(TCOS) )
      IF (NP == 3) K1 = K1 - 1
      K2 = KR
      CALL COSGEN (KR, 1, FNUM2, HALF, TCOS(K1+1),SIZE(TCOS)-K1)
      K4 = K1 + KR
      CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(K4+1),SIZE(TCOS)-K4)
      K3 = LR
      K4 = 0
      CALL TRI3 (MR, vecA, vecBB, vecC, K, vecB, vecB2, vecB3, &
                 TCOS, SIZE(TCOS), vecD, vecW, vecW2, vecW3)
      vecB(:MR) = vecB(:MR) + vecB2(:MR)
      IF (NP == 3) THEN
         TCOS(1) = TWO
         CALL TRIX (1, 0, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                    size(TCOS), vecD, vecW)
      ENDIF
      matQ(:MR,1) = matQ(:MR,1) + vecB(:MR)
  188 CONTINUE
      J = NLAST - JR
      vecB(:MR) = matQ(:MR,NLAST) + matQ(:MR,J)
      JM2 = NLAST - I2R
      IF (JR == 1) THEN
         matQ(:MR,NLAST) = 0.
      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, FNUM2, HALF, TCOS,SIZE(TCOS))
      CALL COSGEN (LR, 1, FNUM2, HALF, TCOS(KR+1),SIZE(TCOS)-KR)
      CALL TRIX (KR, LR, MR, vecA, vecBB, vecC, vecB, TCOS,  &
                 size(TCOS), vecD, vecW)
      matQ(:MR,NLAST) = matQ(:MR,NLAST) + vecB(:MR)
      NLASTP = NLAST
  197 CONTINUE
      JSTEP = JR
      JR = I2R
      I2R = I2R/2
      IF (JR == 0) GO TO 210
      JSTART = 1 + JR
      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,SIZE(TCOS))
      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) = 0.
         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,  &
                    size(TCOS), 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 188
      GO TO 197
  210 CONTINUE
      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 POSTG2