      SUBROUTINE matrix(A,B,C,RHS)

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INTEGER jmaxm1,jmaxm2

      REAL    A(imax,jmax),B(imax,jmax),C(imax,jmax),RHS(imax,jmax),
     +        D(imax,jmax),tmpvar(imax)

c-----------------------------------------------------------------------

c *** Define two work variables for array indices:

      jmaxm1 = jmax - 1
      jmaxm2 = jmax - 2

c*********************************************************************
c
c     Subroutine matrix solves the set of linear equations:
c
c       Matrix (X) = (Right-hand side vektor)
c
c     on input:
c
c         A : matrix elements of the first diagonal of the matrix
c         B : matrix elements of the second diagonal of the matrix
c         C : matrix elements of the third diagonal of the matrix
c       RHS : right-hand side vektor
c
c     on output:
c
c       RHS : solution vektor (X)
c
c     Note: using the i-index as an inner-loop makes vectorization
c           possible without problems.
c
c*********************************************************************

c *** Fill the 2 work array's tmpvar and D:

c     -------------------------------------------start i-loop
      do 100 i=1,imax

c *** In the lower-left corner:

        tmpvar(i) = C(i,jmax)

c *** In the upper right corner:

        D(i,1)    = A(i,1)

100   continue
c     -------------------------------------------start i-loop

c *** Now reduce the matrix to an upper triangular matrix.
c *** Exceptions are made for the indices j=jmax-1 and j=jmax.
c *** They are treated within another loop.

c     -------------------------------------------start j-loop
      do 350 j=1,jmax-3
c     -------------------------------------------start i-loop
         do 450 i=1,imax

c ***    Elimination of element A(i,j+1) ----> A(i,j+1) = 0.0

         B(i,j+1)   =  B(i,j+1)   - (A(i,j+1)/B(i,j))  * C(i,j)

         RHS(i,j+1) =  RHS(i,j+1) - (A(i,j+1)/B(i,j))  * RHS(i,j)

         D(i,j+1)   =             - (A(i,j+1)/B(i,j))  * D(i,j)

c ***    Elimination of element tmpvar(i) ----> tmpvar(i) = 0.0
c ***    (This action moves tmpvar(i) from position j=j to j=j+1)

         RHS(i,jmax) =  RHS(i,jmax)  - (tmpvar(i)/B(i,j)) * RHS(i,j)

         B(i,jmax)   =  B(i,jmax)    - (tmpvar(i)/B(i,j)) * D(i,j)

         tmpvar(i)   =               - (tmpvar(i)/B(i,j)) * C(i,j)

450      continue
c     -------------------------------------------end i-loop
350   continue
c     -------------------------------------------end j-loop

c *** Now we arrive at the lower-right corner of the matrix which
c *** gets a special treatment.

c     -------------------------------------------start i-loop
      do 550 i=1,imax

c *** Elimination of element A(i,jmaxm1) ----> A(i,jmaxm1) = 0.0

      B(i,jmaxm1)   = B(i,jmaxm1)   - ( A(i,jmaxm1)/B(i,jmaxm2) ) *
     +                                  C(i,jmaxm2)

      C(i,jmaxm1)   = C(i,jmaxm1)   - ( A(i,jmaxm1)/B(i,jmaxm2) ) *
     +                                  D(i,jmaxm2)

      RHS(i,jmaxm1) = RHS(i,jmaxm1) - ( A(i,jmaxm1)/B(i,jmaxm2) ) *
     +                                  RHS(i,jmaxm2)

c *** Now the last row of the matrix is updated because of the
c *** periodicity of the problem.

c *** Elimination of element tmpvar(i) ----> tmpvar(i) = 0.0
c *** (This action removes tmpvar(i) from position j=jmaxm2)

      A(i,jmax)   = A(i,jmax)   - (tmpvar(i)/B(i,jmaxm2))*C(i,jmaxm2)

      B(i,jmax)   = B(i,jmax)   - (tmpvar(i)/B(i,jmaxm2))*D(i,jmaxm2)

      RHS(i,jmax) = RHS(i,jmax) - (tmpvar(i)/B(i,jmaxm2))*RHS(i,jmaxm2)

c *** As a last step in the procedure of reducing the matrix to an
c *** upper triangular matrix, the last row is now treated in the
c *** "standard" way.

c *** Final elimination of element A(i,jmax) ----> A(i,jmax) = 0.0

      B(i,jmax)   = B(i,jmax)   - (A(i,jmax)/B(i,jmaxm1))*C(i,jmaxm1)

      RHS(i,jmax) = RHS(i,jmax) - (A(i,jmax)/B(i,jmaxm1))*RHS(i,jmaxm1)

c *** The matrix is now completely reduced to an upper matrix and the
c *** solution can be obtained by backward substitution starting with
c *** the element with the highest index (j=jmax):

      RHS(i,jmax) = RHS(i,jmax) / B(i,jmax)

c *** Calculate the solution member with the one-but-highest index:

      RHS(i,jmaxm1) = ( RHS(i,jmaxm1) - C(i,jmaxm1) * RHS(i,jmax) ) /
     +                  B(i,jmaxm1)

550   continue
c     -------------------------------------------end i-loop

c *** All other solution members (j=jmax-2 until j=1) can be obtained
c *** within one loop:

      DO j=jmax-2,1,-1
         DO i=1,imax

            RHS(i,j) = ( RHS(i,j) - C(i,j)*RHS(i,j+1) - D(i,j)*RHS(i
     $           ,jmax) )/ B(i,j)

         END DO
      END DO

c *** Now, the array RHS(i,j) contains the solution.

      RETURN
      END
