SUBROUTINE ZSPLDR(X,NPTS,XA,YA,N,Y,DERIV,STATUS)
*
*  Module number:
*
*  Module name: ZSPLDR
*
*  Keyphrase:
*  ----------
*       cubic spline interpolation
*  Description:
*  ------------
*       This routine computes the cubic spline interpolation using
*       Numerical recipe routines spline and splint. It also numerically
*       computes the partial derivatives of y with respect to ya
*
*  FORTRAN name: zsplin.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*       none
*  Subroutines Called:
*  -------------------
*  CDBS:
*       zsplin
*  SDAS:
*       umsput
*  Others:
*
*
*  History:
*  --------
*  Version      Date        Author          Description
*    1          Dec 87  D. Lindler      Designed and coded
*-------------------------------------------------------------------------------
C
C INPUTS
C       X - VECTOR OF X-POINTS TO INTERPOLATE FOR (REAL*8)
C       NPTS - NUMBER OF POINTS IN X (INTEGER)
C       XA - VECTOR OF X-NODE POSITIONS OF THE SPLINE (REAL*8)
C       YA - VECTOR OF Y-NODE POSITIONS OF THE SPLINE (REAL*8)
C       N - NUMBER OF POINTS IN X AND Y
C
C OUTPUTS
C
C       Y - VECTOR OF INTERPOLATED VALUES OF LENGTH N (REAL*8)
C       DERIV - ARRAY OF PARTIAL DERIVIATIVES REAL*8 OF SIZE (30,2200)
C               DERIV(I,J) HAS THE DERIVATIVE OF Y(J) WITH RESPECT TO YA(I)
C       STATUS - ERROR STATUS
C
C     INCLUDE FILE FOR THE IRAF77 FORTRAN INTERFACE TO THE IRAF VOS
C
C
C     FILE I/O ACCESS MODES
C
      INTEGER   RDONLY
      PARAMETER (RDONLY = 1)
      INTEGER   RDWRIT
      PARAMETER (RDWRIT = 2)
      INTEGER   WRONLY
      PARAMETER (WRONLY = 3)
      INTEGER   APPEND
      PARAMETER (APPEND = 4)
C
C     CODES FOR DATA TYPES
C
      INTEGER   TYBOOL
      PARAMETER (TYBOOL = 1)
      INTEGER   TYCHAR
      PARAMETER (TYCHAR = 2)
      INTEGER   TYINT
      PARAMETER (TYINT = 4)
      INTEGER   TYREAL
      PARAMETER (TYREAL = 6)
      INTEGER   TYDOUB
      PARAMETER (TYDOUB = 7)
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C     UHDAS HEADER PARM TYPES -- CB, DAO, 5-SEP-87
C
      INTEGER GENHDR
      PARAMETER (GENHDR = 0)
      INTEGER IMSPEC
      PARAMETER (IMSPEC = 1)
C
C     THIS SECTION IS FOR PARAMETERS RELEVANT TO TABLE I/O.
C
C     THESE MAY BE SET BY UTPPTI AND/OR READ BY UTPGTI:
C
C                                       LENGTH OF ROW (UNIT = SIZE OF REAL)
      INTEGER   TBRLEN
      PARAMETER (TBRLEN = 1)
C                                       INCREASE ROW LENGTH
      INTEGER   TBIRLN
      PARAMETER (TBIRLN = 2)
C                                       NUMBER OF ROWS TO ALLOCATE
      INTEGER   TBALLR
      PARAMETER (TBALLR = 3)
C                                       INCREASE ALLOC NUM OF ROWS
      INTEGER   TBIALR
      PARAMETER (TBIALR = 4)
C                                       WHICH TYPE OF TABLE? (ROW OR COLUMN)
      INTEGER   TBWTYP
      PARAMETER (TBWTYP = 5)
C                                       MAXIMUM NUMBER OF USER PARAMETERS
      INTEGER   TBMXPR
      PARAMETER (TBMXPR = 6)
C                                       MAXIMUM NUMBER OF COLUMNS
      INTEGER   TBMXCL
      PARAMETER (TBMXCL = 7)
C                                       TYPE = ROW-ORDERED TABLE
      INTEGER   TBTYPR
      PARAMETER (TBTYPR = 11)
C                                       TYPE = COLUMN-ORDERED TABLE
      INTEGER   TBTYPC
      PARAMETER (TBTYPC = 12)
C
C     THESE MAY BE READ BY UTPGTI BUT MAY NOT BE SET:
C
C                                       NUMBER OF ROWS WRITTEN TO
      INTEGER   TBNROW
      PARAMETER (TBNROW = 21)
C
C     END IRAF77.INC
*
        DOUBLE PRECISION X(1),Y(1),XA(1),YA(1),DERIV(30,2200)
        INTEGER N,NPTS
        INTEGER STATUS,ISTAT
*
* LOCAL VARIABLES
*
        INTEGER I,J
        CHARACTER*130 CONTXT
        DOUBLE PRECISION YAPLUS(30),DIF
        DOUBLE PRECISION Y2(2500)
*--------------------------------------------------------------------------------
C
C CHECK FOR VALID SIZES
C
        IF(NPTS.GT.2200)THEN
                STATUS=1
                CONTXT='NUMBER OF DATA POINTS CAN NOT EXCEED 2200'
                GO TO 999
        ENDIF
C
        IF(N.GT.30)THEN
                CONTXT='NUMBER OF NODES IN SPLINE CAN NOT EXCEED 30'
                GO TO 999
        ENDIF
C
C COMPUTE FUNCTION VALUE AT XA,YA
C
        CALL ZSPLIN(X,NPTS,XA,YA,N,Y,STATUS)
        IF(STATUS.NE.0)THEN
                CONTXT='UNABLE TO COMPUTE SPLINE FUNCTION'
                STATUS=1
                GO TO 999
        ENDIF
C
C LOOP ON NUMBER OF NODES AND COMPUTE NUMERICAL PARTIAL DERIVATIVES
C       FOR CHANGES TO EACH YA(J)
C
        DO 100 J=1,N
C
C COPY PRESENT YA VALUES TO YAPLUS
C
                DO 10 I=1,N
10                      YAPLUS(I)=YA(I)
C
C INCREASE YAPLUS(J) BY 0.1%
C
                DIF=ABS(YAPLUS(J))*0.0001
                IF (DIF.LT.1.0E-20) DIF=1.0E-20
                YAPLUS(J)=YAPLUS(J)+DIF
C
C COMPUTE FUNCTION VALUE AT YAPLUS
C
                CALL ZSPLIN(X,NPTS,XA,YAPLUS,N,Y2,STATUS)
                IF(STATUS.NE.0)THEN
                  CONTXT='ERROR COMPUTING SPLINE INTERPOLATION VALUES'
                  STATUS=1
                  GO TO 999
                ENDIF
C
C COMPUTE NUMERICAL DERIVATIVES
C
                DO 20 I=1,NPTS
                        DERIV(J,I)=(Y2(I)-Y(I))/DIF
20              CONTINUE
100     CONTINUE
        STATUS=0
        GO TO 1000
999     CALL UMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
1000    RETURN
        END