      SUBROUTINE INTERP (RJD,X,Y,T,N,RJDINT,XINT,YINT,TINT)
C
C     THIS SUBROUTINE TAKES A SERIES OF X, Y, AND UT1-UTC VALUES
C     AND INTERPOLATES THEM TO AN EPOCH OF CHOICE.  THIS ROUTINE
C     ASSUMES THAT THE VALUES OF X AND Y ARE IN SECONDS OF
C     ARC AND THAT UT1-UTC IS IN SECONDS OF TIME.  AT LEAST
C     ONE POINT BEFORE AND ONE POINT AFTER THE EPOCH OF THE
C     INTERPOLATION POINT ARE NECESSARY IN ORDER FOR THE
C     INTERPOLATION SCHEME TO WORK.
C
C     PARAMETERS ARE :
C     RJD   - ARRAY OF THE EPOCHS OF DATA (GIVEN IN MJD)
C     X     - ARRAY OF X POLAR MOTION (ARCSEC)
C     Y     - ARRAY OF Y POLAR MOTION (ARCSEC)
C     T     - ARRAY OF UT1-UTC (SEC)
C     N     - NUMBER OF POINTS IN ARRAYS
C     RJDINT- EPOCH FOR THE INTERPOLATED VALUE
C     XINT  - INTERPOLATED VALUE OF X
C     YINT  - INTERPOLATED VALUE OF Y
C     TINT  - INTERPOLATED VALUE OF UT1-UTC
C
      DOUBLE PRECISION RJD(N), X(N), Y(N), T(N),
     . RJDINT, XINT, YINT, TINT, CORX, CORY, CORT
      CALL LAGINT (RJD,X,N,RJDINT,XINT)
      CALL LAGINT (RJD,Y,N,RJDINT,YINT)
      CALL LAGINT (RJD,T,N,RJDINT,TINT)
      CALL RAY (RJDINT,CORX,CORY,CORT)
      XINT = XINT + CORX
      YINT = YINT + CORY
      TINT = TINT + CORT
      RETURN
      END
C
C----------------------------------------------------------------
C
      SUBROUTINE LAGINT (X,Y,N,XINT,YOUT)
C 
C     THIS SUBROUTINE PERFORMS LAGRANGIAN INTERPOLATION
C     WITHIN A SET OF (X,Y) PAIRS TO GIVE THE Y
C     VALUE CORRESPONDING TO XINT.  THIS PROGRAM USES A
C     WINDOW OF 4 DATA POINTS TO PERFORM THE INTERPOLATION.
C     IF THE WINDOW SIZE NEEDS TO BE CHANGED, THIS CAN BE
C     DONE BY CHANGING THE INDICES IN THE DO LOOPS FOR
C     VARIABLES M AND J.
C
C     PARAMETERS ARE :
C     X     - ARRAY OF VALUES OF THE INDEPENDENT VARIABLE
C     Y     - ARRAY OF FUNCTION VALUES CORRESPONDING TO X
C     N     - NUMBER OF POINTS
C     XINT  - THE X-VALUE FOR WHICH ESTIMATE OF Y IS DESIRED
C     YOUT  - THE Y VALUE RETURNED TO CALLER
C
      REAL*8 X(N),Y(N),XINT,YOUT,TERM
      INTEGER N,I,J
C
      YOUT = 0.0D0
      DO 5 I = 1,N-1
        IF ( XINT .GE. X(I) .AND. XINT .LT. X(I+1) ) K = I
    5 CONTINUE
      IF ( K .LT. 2 ) K = 2
      IF ( K .GT. N-2 ) K = N-2
      DO 20 M = K-1,K+2
        TERM = Y(M)
        DO 10 J = K-1,K+2
          IF ( M .NE. J ) THEN
            TERM = TERM * (XINT - X(J))/(X(M) - X(J))
          END IF
   10   CONTINUE
        YOUT = YOUT + TERM
   20 CONTINUE
      RETURN
      END

C
C----------------------------------------------------------------
C
      SUBROUTINE RAY (RJD,CORX,CORY,CORT)
C
C   THIS SUBROUTINE IMPLEMENTS THE RAY MODEL FOR
C   DIURNAL/SUBDIURNAL TIDES.  IT USES THE SIMON ET AL.
C   FUNDAMENTAL ARGUMENTS.  THE CORRECTIONS IN X AND Y ARE IN
C   UNITS OF SEC. OF ARC AND UT1-UTC IN SEC. OF TIME.  THESE
C   CORRECTIONS SHOULD BE ADDED TO "AVERAGE" EOP VALUES TO GET
C   ESTIMATES OF THE INSTANTANEOUS VALUES.
C
C     PARAMETERS ARE :
C     RJD   - EPOCH OF INTEREST GIVEN IN MJD
C     CORX  - TIDAL CORRECTION IN X (SEC. OF ARC)
C     CORY  - TIDAL CORRECTION IN Y (SEC. OF ARC)
C     CORT  - TIDAL CORRECTION IN UT1-UTC (SEC. OF TIME)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION
     .   L,        LPRIME
      HALFPI = 1.5707963267948966d0
      T = (RJD - 51544.5D0)/36525.0D0
      L = -0.00024470d0*T**4 + 0.051635d0*T**3 + 31.8792d0*T**2
     .  + 1717915923.2178d0*T + 485868.249036d0
      L = DMOD(L,1296000d0)
      LPRIME = -0.00001149d0*T**4 - 0.000136d0*T**3
     .  -  0.5532d0*T**2
     .  + 129596581.0481d0*T + 1287104.79305d0
      LPRIME = DMOD(LPRIME,1296000d0)
      CAPF = 0.00000417d0*T**4 - 0.001037d0*T**3 - 12.7512d0*T**2
     .  + 1739527262.8478d0*T + 335779.526232d0
      CAPF = DMOD(CAPF,1296000d0)
      CAPD = -0.00003169d0*T**4 + 0.006593d0*T**3 - 6.3706d0*T**2
     .  + 1602961601.2090d0*T + 1072260.70369d0
      CAPD = DMOD(CAPD,1296000d0)
      OMEGA = -0.00005939d0*T**4 + 0.007702d0*T**3
     .  + 7.4722d0*T**2
     .  - 6962890.2665d0*T + 450160.398036d0
      OMEGA = DMOD(OMEGA,1296000d0)
      THETA = (67310.54841d0 +
     .        (876600d0*3600d0 + 8640184.812866d0)*T +
     .         0.093104d0*T**2 -
     .         6.2d-6*T**3)*15.0d0 + 648000.0d0
      ARG7 = DMOD((-L - 2.0D0*CAPF - 2.0D0*OMEGA + THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0) - HALFPI
      ARG1 = DMOD((-2.0d0*CAPF - 2.0d0*OMEGA + THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0) - HALFPI
      ARG2 = DMOD((-2.0d0*CAPF + 2.0d0*CAPD - 2.0d0*OMEGA
     .       + THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0) - HALFPI
      ARG3 = DMOD(THETA *
     .        3.14159265D0/648000.0D0,6.28318530718D0)
     .     + HALFPI
      ARG4 = DMOD((-L - 2.0d0*CAPF - 2.0D0*OMEGA + 2.0d0*THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0)
      ARG5 = DMOD((-2.0D0*CAPF - 2.0D0*OMEGA + 2.0d0*THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0)
      ARG6 = DMOD((-2.0d0*CAPF + 2.0d0*CAPD - 2.0d0*OMEGA
     .     + 2.0d0*THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0)
      ARG8 = DMOD((2.0d0*THETA)
     .     * 3.14159265D0/648000.0D0,6.28318530718D0)
      CORX = - 0.026D0*DSIN(ARG7) + 0.006D0*DCOS(ARG7)
     .       - 0.133D0*DSIN(ARG1) + 0.049D0*DCOS(ARG1)
     .       - 0.050D0*DSIN(ARG2) + 0.025D0*DCOS(ARG2)
     .       - 0.152D0*DSIN(ARG3) + 0.078D0*DCOS(ARG3)
     .       - 0.057D0*DSIN(ARG4) - 0.013D0*DCOS(ARG4)
     .       - 0.330D0*DSIN(ARG5) - 0.028D0*DCOS(ARG5)
     .       - 0.145D0*DSIN(ARG6) + 0.064D0*DCOS(ARG6)
     .       - 0.036D0*DSIN(ARG8) + 0.017D0*DCOS(ARG8)
      CORY = - 0.006D0*DSIN(ARG7) - 0.026D0*DCOS(ARG7)
     .       - 0.049D0*DSIN(ARG1) - 0.133D0*DCOS(ARG1)
     .       - 0.025D0*DSIN(ARG2) - 0.050D0*DCOS(ARG2)
     .       - 0.078D0*DSIN(ARG3) - 0.152D0*DCOS(ARG3)
     .       + 0.011D0*DSIN(ARG4) + 0.033D0*DCOS(ARG4)
     .       + 0.037D0*DSIN(ARG5) + 0.196D0*DCOS(ARG5)
     .       + 0.059D0*DSIN(ARG6) + 0.087D0*DCOS(ARG6)
     .       + 0.018D0*DSIN(ARG8) + 0.022D0*DCOS(ARG8)
      CORT = + 0.0245D0*DSIN(ARG7) + 0.0503D0*DCOS(ARG7)
     .       + 0.1210D0*DSIN(ARG1) + 0.1605D0*DCOS(ARG1)
     .       + 0.0286D0*DSIN(ARG2) + 0.0516D0*DCOS(ARG2)
     .       + 0.0864D0*DSIN(ARG3) + 0.1771D0*DCOS(ARG3)
     .       - 0.0380D0*DSIN(ARG4) - 0.0154D0*DCOS(ARG4)
     .       - 0.1617D0*DSIN(ARG5) - 0.0720D0*DCOS(ARG5)
     .       - 0.0759D0*DSIN(ARG6) - 0.0004D0*DCOS(ARG6)
     .       - 0.0196D0*DSIN(ARG8) - 0.0038D0*DCOS(ARG8)
      CORX = CORX * 1.0d-3
      CORY = CORY * 1.0d-3
      CORT = CORT * 0.1d-3
      RETURN
      END
