       subroutine polint ( xa, ya, n, x, y, dy )
c
c    dac:zeus3d.polint <---------------------- interpolates gridded data
c    from whp:numerical recipes                           december, 1992
c
c    written by: David Clarke
c    modified 1:
c
c  PURPOSE:  Given arrays "xa" and "ya", each of length "n", and a given
c  value "x", this routine returns a value "y" and an error estimate
c  "dy" (see Numerical Recipes, 1st edition for FORTRAN, page 82).
c
c  INPUT VARIABLES:
c    xa       input independent variable array
c    ya       input dependent   variable array
c    n        length of "xa" and "ya"
c    x        value of independent variable at which dependent variable
c             is to be interpolated
c
c  OUTPUT VARIABLES:
c    y        interpolated value of dependent variable
c    dy       estimate of error of interpolation.
c
c  EXTERNALS: [NONE]
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c------------------------ IMPLICIT STATEMENTS --------------------------
c-----------------------------------------------------------------------
c
       implicit      none
c
c-----------------------------------------------------------------------
c----------------------------- PARAMETERS ------------------------------
c-----------------------------------------------------------------------
c
c      Primary parameters:
c
       integer       in      , jn      , kn      , ijkn    , nxpx
     1             , nypx    , nxrd    , nyrd    , niov    , nios
     2             , ntsl    , mn      , nmd     , isig
       real*8        pi  , tiny  , huge  , one  , mone  , haf  , zro
c
       parameter     ( in = 5005, jn = 1, kn = 10, ijkn = 5005)
       parameter     ( nxpx = 1, nypx = 1, nxrd =   1, nyrd =   1 )
       parameter     ( niov =  20, nios =  20, ntsl =1000 )
       parameter     ( mn   =   1, nmd  =   1 )
       parameter     ( pi   = 3.1415926535898 )
       parameter     ( isig = 9, tiny = 1.0d-99, huge = 1.0d+99 )
       parameter     ( one  = 1.0d+00, mone = -1.0d+00
     1                               , haf  =  0.5d+00, zro = 0.0d+00)
c
c      Secondary parameters:
c
       integer        idim    , jdim
     1              , nt01    , nt02    , nt03    , nt04    , nt05
     2              , nt06    , nt07    , nt08    , nt09    , nt10
     3              , nt11    , nt12    , nt13    , nt14    , nt15
     4              , nt16    , nt17    , nt18    , nt19    , nt20
       parameter     ( idim =   kn, jdim =   in )
       parameter     ( nt01 = 600, nt02 = 600, nt03 =   7, nt04 =   5 )
       parameter     ( nt05 =  18, nt06 =   1, nt07 =  10, nt08 =  12 )
       parameter     ( nt09 =   5, nt10 =  54, nt11 =   8 )
       parameter     ( nt12 =        16, nt17 =         4 )
       parameter     ( nt13 = nt12 +  4, nt18 = nt17 +  1 )
       parameter     ( nt14 = nt13     , nt19 = nt18      )
       parameter     ( nt15 = nt14 + 30, nt20 = nt19 +  3 )
       parameter     ( nt16 = nt14 + 42 )
c
       integer       nmax
       parameter     ( nmax=20 )
c
       integer       n       , ns      , i       , m
       real*8        x       , y       , dy      , dif     , dift
     1             , ho      , hp      , w       , den
c
       real*8        xa      (   n), ya      (   n)
     1             , c       (nmax), d       (nmax)
c
c-----------------------------------------------------------------------
c
       ns  = 1
       dif = abs ( x - xa(1) )
       do 10 i=1,n
         dift = abs ( x - xa(i) )
         if (dift .lt. dif) then
           ns  = i
           dif = dift
         endif
         c(i) = ya(i)
         d(i) = ya(i)
10     continue
       y  = ya(ns)
       ns = ns - 1
       do 30 m=1,n-1
         do 20 i=1,n-m
           ho   = xa(i  ) - x
           hp   = xa(i+m) - x
           w    = c (i+1) - d(i)
           den  = w / ( ( ho - hp ) + tiny )
           c(i) = ho * den
           d(i) = hp * den
20       continue
         if (2*ns .lt. n-m) then
           dy = c(ns+1)
         else
           dy = d(ns  )
           ns = ns - 1
         endif
         y = y + dy
30     continue
c
       return
       end
