
! From settei: kf=kmaxp1, k1,k2=1,kmax
!     CALL TRSOLV(S1,S2,S3,S4,S7,F(1,NTEK),LEN1,3,LEN1-2,KMAXP1,1,KMAX
!    1,1,j,0)

      SUBROUTINE TRSOLV(A,B,C,F,W,X,IF,I1,I2,KF,K1,K2,NF,j,iprint)
      implicit none
C     ****     SOLVES TRIDIAGONAL SYSTEM
C     ****     A(I,K)*X(I,K-1)+B(I,K)*X(I,K)+C(I,K)*X(I,K+1)=F(I,K)
C     ****     WHERE I=I1,I2,1 AND K=K1,K2,1 AND A(I,K1)=C(I,K2)=0.
C     ****     ARRAYS A,B,C,F,W,X ARE DIMENSIONED (IF,KF)
C     ****     WHERE 1.LE.I1.LT.I2.LE.IF AND 1.LE.K1.LT.K2.LE.KF
C     ****     W IS WORK SPACE
C     ****     THE ROUTINE USES VECTOR OPERATIONS ON ARRAYS OF LENGTH
C     ****     IMAX=I2-I1+1
! Args:
      integer,intent(in) :: if,i1,i2,kf,k1,k2,nf,j,iprint
      real,intent(in) :: A(IF,KF),B(IF,KF),C(IF,KF),F(IF,KF)
      real,intent(out) :: W(IF,KF),X(IF,KF)
!
! Local:
      integer :: imax,k1p,i,k,kk
      integer :: zimxp,zkmxp,zkmx ! for addfsech
!
      zimxp=(i2+2)-(i1-2)+1

!     if (iprint > 0) 
!    |  write(6,"('trsolv: lat=',i3,' k1,k2=',2i3,' kf=',i3,' zimxp=',
!    |    i3,' zkmxp=',i3,' zkmx=',i3)") j,k1,k2,kf,zimxp,zkmxp,zkmx
!     if (iprint > 0) then
!       write(6,"('trsolv: j=',i3,' kf=',i3)") j,kf
!       write(6,"('a(:,kf)=',/,(6e12.4))") a(:,kf)
!       write(6,"('b(:,kf)=',/,(6e12.4))") b(:,kf)
!       write(6,"('c(:,kf)=',/,(6e12.4))") c(:,kf)
!       write(6,"('f(:,kf)=',/,(6e12.4))") f(:,kf)
!     endif 

      IMAX=I2-I1+1
      K1P=K1+1
      GO TO(1,2),NF
    1 CONTINUE
C     ****          LOWER BOUNDARY
C     ****          W(K1)=B(K1)
      DO 3 I=I1,I2
      W(I,K1)=B(I,K1)
    3 CONTINUE
C     ****          SWEEP THROUGH K
      DO 4 K=K1P,K2
C     ****          W(KF+K-1)=C(K-1)/W(K-1)
        do i=1,imax
          w(i1+i-1,kf+k-1) = c(i1+i-1,k-1) / w(i1+i-1,k-1)
        enddo
C     ****          W(K)=A(K)*W(KF+K-1)
        do i=1,imax
          w(i1+i-1,k) = a(i1+i-1,k) * w(i1+i-1,kf+k-1)
        enddo
C     ****          W(K)=B(K)-W(K)
        do i=1,imax
          w(i1+i-1,k) = b(i1+i-1,k) - w(i1+i-1,k)
        enddo
    4 CONTINUE

!     if (iprint > 0) then
!       call addfsech('W1',' ',' ',w        ,zimxp,kf,kf-1,j)
!       call addfsech('W2',' ',' ',w(1,kf+1),zimxp,kf,kf-1,j)
!     endif

    2 CONTINUE
C     ****          LOWER BOUNDARY
C     ****          W(2*KF+K1)=F(K1)/W(K1)
      do i=1,imax
        w(i1+i-1,2*kf+k1) = f(i1+i-1,k1) / w(i1+i-1,k1)
      enddo
      DO 5 K=K1P,K2
C     ****          W(2*KF+K)=A(K)*W(2*KF+K-1)
      do i=1,imax
        w(i1+i-1,2*kf+k) = a(i1+i-1,k) * w(i1+i-1,2*kf+k-1)
      enddo
      
C     ****          W(2*KF+K)=F(K)-W(2*KF+K)
      do i=1,imax
        w(i1+i-1,2*kf+k) = f(i1+i-1,k) - w(i1+i-1,2*kf+k)
      enddo
C     ****          W(2*KF+K)=W(2*KF+K)/W(K)
      do i=1,imax
        w(i1+i-1,2*kf+k) = w(i1+i-1,2*kf+k) / w(i1+i-1,k)
      enddo
    5 CONTINUE

!     if (iprint > 0) 
!    |  call addfsech('W3',' ',' ',w(1,2*kf+1),zimxp,kf,kf-1,j)

C     ****          BACK SUBSTITUTION
C     ****          X(K2)=W(2*KF+K2)
      DO 6 I=I1,I2
      X(I,K2)=W(I,2*KF+K2)
    6 CONTINUE
      DO 7 KK=K1P,K2
      K=K1+K2-KK
C     ****          X(K)=W(KF+K)*X(K+1)
      do i=1,imax
        x(i1+i-1,k) = w(i1+i-1,kf+k) * x(i1+i-1,k+1)
      enddo
C     ****          X(K)=W(2*KF+K)-X(K)
      do i=1,imax
        x(i1+i-1,k) = w(i1+i-1,2*kf+k) - x(i1+i-1,k)
      enddo
    7 CONTINUE

!     if (iprint > 0) 
!    |  call addfsech('X'  ,' ',' ',x,zimxp,kf,kf-1,j)

      RETURN
      END
C
