SUBROUTINE RATVA(A,X,Y,NXP,NYP,NX,NY,PP,XX,YY,F)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     Computation of the interpolated functional value at (XX,YY)
C     using birational splines.
C
C  INPUT
C     A        Coefficient matrix (NX,NY,4,4).
C     X,Y      Coordinate vectors.
C     NXP,NYP  Physical length of the above.
C     NX,NY    Length of the above to be processed.
C     PP       Parameter for birational spline interpolation as
C              described in ARAT2D.
C     XX,YY    Coordinates of the desired point.
C
C  OUTPUT
C     F        Interpolated functional value at the desired point.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL A(NXP,NYP,4,4),X(NXP),Y(NYP),H(4,4),HH(4,4)
C
C...LOCALIZE THE X INDEX OF THE POINT
C
      I=1
      J=1
      IF (XX.LT.X(I)) GOTO 10
      IF (XX.LE.X(I+1)) GOTO 40
      L=NX
      GOTO 30
   10 L=I
      I=1
   20 K=(I+L)/2
      IF (XX.LT.X(K)) THEN
         L=K
      ELSE
         I=K
      ENDIF
   30 IF (L.GT.I+1) GOTO 20
C
C...LOCALIZE THE Y INDEX OF THE POINT
C
   40 IF (YY.LT.Y(J)) GOTO 50
      IF (YY.LE.Y(J+1)) GOTO 80
      L=NY
      GOTO 70
   50 L=J
      J=1
   60 K=(J+L)/2
      IF (YY.LT.Y(K)) THEN
         L=K
      ELSE
         J=K
      ENDIF
   70 IF (L.GT.J+1) GOTO 60
C
C...COMPUTE THE COEFFICIENTS OF THE INTERPOLATION
C
   80 XXX=(XX-X(I))/(X(I+1)-X(I))
      YYY=(YY-Y(J))/(Y(J+1)-Y(J))
      H1=1.-XXX
      H2=1.-YYY
      DO 90 L=1,4
         H(1,L)=H1
         H(2,L)=XXX
         H(3,L)=(H1*H1*H1)/(PP*XXX+1.)
         H(4,L)=(XXX*XXX*XXX)/(PP*H1+1.)
         HH(L,1)=H2
         HH(L,2)=YYY
         HH(L,3)=(H2*H2*H2)/(PP*YYY+1.)
         HH(L,4)=(YYY*YYY*YYY)/(PP*H2+1.)
   90 CONTINUE
C
C...INTERPOLATE AT XX AND YY
C
      F=0.
      DO 100 K=1,4
         DO 110 L=1,4
            F=F+A(I,J,K,L)*H(K,L)*HH(K,L)
  110    CONTINUE
  100 CONTINUE
      RETURN
      END