SUBROUTINE RATPE(U,MX,MY,P,RX,NNX,NNY,IVJ,IPU,P3,
     &                 AX,BX,CX,DX,NX,NY)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Solves the equation system for the subroutine ARAT2D.
C
C  INPUT
C     U          First matrix (MX,MY).
C     MX,MY      Physical length of the above.
C     P          Matrix (NNX,NNY) of last sweep.
C     NNX,NNY    Physical length of P.
C     IVJ,IPU    Counting parameters.
C     P3         Parameter of the birational interpolation.
C     AX,BX,CX   Diagonals of L.
C     DX         Grid vector.
C     NX,NY      Length of the calculation.
C
C  OUTPUT
C     P          Second matrix (NNX,NNY).
C     RX         Solution vector (NNX).
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL U(MX,MY),P(NNX,NNY),AX(NNX),BX(NNX),CX(NNX),RX(NNX),
     &     DX(NNX)
      real R(NY,NX), Rxx(NY,NNX)
      integer, parameter :: BLKSZ=64
C
      NX2=NX-2
      IF (IVJ.NE.0) THEN
!     ipu is normally 1 but sometimes ny-1
!     actually the only time ipu is ny-1, ivj is 0 
        if (IPU.eq.1) then
          do j1=1,NY,BLKSZ
            do i=1,NX-1
              do j=j1,min(NY,j1+BLKSZ)
                r(j,i) = P3*DX(i)*DX(i)*(U(j,i+1)-U(j,i))
              end do
            end do
          end do
        else ! ipu will be ny-1, so only 2 points of j
          do j=1,NY,IPU
            do i=1,NX-1
              r(j,i) = P3*DX(i)*DX(i)*(U(j,i+1)-U(j,i))
            end do
          end do
        end if
! it is independent in j, can i take advantage of that?
! rx gets value at final j (but it is not used by calling routine)
        DO J=1,NY,IPU
          rxx(j,1) = r(j,1) + r(j,2) - DX(1)*P(J,1)
        end do
        do I=3,NX-2
          DO J=1,NY,IPU
            RXx(j,i-1) = r(j,i-1) + r(j,i)
          end do
        end do
        DO J=1,NY,IPU
          RXx(j,nx-2) = r(j,nx-2) + r(j,nx-1) -DX(NX-1)*P(J,NX)
          RXx(j,1)=RXx(j,1)/BX(1)
        end do
        do K=2,NX2
          DO J=1,NY,IPU
            RXx(j,K)=(RXx(j,K)-AX(K-1)*RXx(j,K-1))/BX(K)
          end do
        end do
        do K=NX2-1,1,-1
          DO J=1,NY,IPU
            RXx(j,K)=RXx(j,K)-CX(K)*RXx(j,K+1)
          end do
        end do
        do I=2,NX-1
          DO J=1,NY,IPU
            P(J,I)=RXx(j,I-1)
          end do
        end do

!         DO 10 J=1,NY,IPU
!            R2=P3*DX(1)*DX(1)*(U(J,2)-U(J,1))
!            R1=R2
!            R2=P3*DX(2)*DX(2)*(U(J,3)-U(J,2))
!            RX(1)=R1+R2
!            RX(1)=RX(1)-DX(1)*P(J,1)
!            R1=R2
!            DO 20 I=3,NX-2
!               I1=I-1
!               R2=P3*DX(I)*DX(I)*(U(J,I+1)-U(J,I))
!               RX(I1)=R1+R2
!               R1=R2
!   20       CONTINUE
!            R2=P3*DX(NX-1)*DX(NX-1)*(U(J,NX)-U(J,NX-1))
!            RX(NX2)=R1+R2
!            RX(NX2)=RX(NX2)-DX(NX-1)*P(J,NX)
!            RX(1)=RX(1)/BX(1)
!            DO 30 K=2,NX2
!               KM1=K-1
!               RX(K)=(RX(K)-AX(KM1)*RX(KM1))/BX(K)
!   30       CONTINUE
!            DO 40 K=NX2-1,1,-1
!               RX(K)=RX(K)-CX(K)*RX(K+1)
!   40       CONTINUE
!            DO 50 I=2,NX-1
!               P(J,I)=RX(I-1)
!   50       CONTINUE
!   10    CONTINUE
      ELSE
! ipu is either 1 or ny-1
        if (IPU.eq.1) then
          do j1=1,NY,BLKSZ
            do i=1,NX-1
              do j=j1,min(NY,j1+BLKSZ)
                r(j,i) = P3*DX(i)*DX(i)*(U(i+1,j)-U(i,j))
              end do
            end do
          end do
        else ! ipu will be ny-1, so only 2 points of j
          do j=1,NY,IPU
            do i=1,NX-1
              r(j,i) = P3*DX(i)*DX(i)*(U(i+1,j)-U(i,j))
            end do
          end do
        end if

        DO J=1,NY,IPU
          rxx(j,1) = r(j,1) + r(j,2) - DX(1)*P(1,J)
        end do
        do I=3,NX-2
          DO J=1,NY,IPU
            RXx(j,i-1) = r(j,i-1) + r(j,i)
          end do
        end do
        DO J=1,NY,IPU
          RXx(j,nx-2) = r(j,nx-2) + r(j,nx-1) -DX(NX-1)*P(NX,J)
          RXx(j,1)=RXx(j,1)/BX(1)
        end do
        do K=2,NX2
          DO J=1,NY,IPU
            RXx(j,K)=(RXx(j,K)-AX(K-1)*RXx(j,K-1))/BX(K)
          end do
        end do
        do K=NX2-1,1,-1
          DO J=1,NY,IPU
            RXx(j,K)=RXx(j,K)-CX(K)*RXx(j,K+1)
          end do
        end do
        do i1=2,NX-1,BLKSZ
        DO J=1,NY,IPU
          do I=i1,min(BLKSZ,NX-1)
            P(I,J)=RXx(j,I-1) 
          end do
        end do
        end do


!         DO 60 J=1,NY,IPU
!            R2=P3*DX(1)*DX(1)*(U(2,J)-U(1,J))
!            R1=R2
!            R2=P3*DX(2)*DX(2)*(U(3,J)-U(2,J))
!            RX(1)=R1+R2
!            RX(1)=RX(1)-DX(1)*P(1,J)
!            R1=R2
!            DO 70 I=3,NX-2
!               I1=I-1
!               R2=P3*DX(I)*DX(I)*(U(I+1,J)-U(I,J))
!               RX(I1)=R1+R2
!               R1=R2
!   70       CONTINUE
!            R2=P3*DX(NX-1)*DX(NX-1)*(U(NX,J)-U(NX-1,J))
!            RX(NX2)=R1+R2
!            RX(NX2)=RX(NX2)-DX(NX-1)*P(NX,J)
!            RX(1)=RX(1)/BX(1)
!            DO 80 K=2,NX2
!               KM1=K-1
!               RX(K)=(RX(K)-AX(KM1)*RX(KM1))/BX(K)
!   80       CONTINUE
!            DO 90 K=NX2-1,1,-1
!               RX(K)=RX(K)-CX(K)*RX(K+1)
!   90       CONTINUE
!            DO 100 I=2,NX-1
!               P(I,J)=RX(I-1)
!  100       CONTINUE
!   60    CONTINUE
      ENDIF
      RETURN
      END