1      subroutine DILUP (X, Y, NTAB, XT, YT, NDEG, LUP, IOPT, EOPT)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>> 2010-08-26 DILUP  Krogh  Dimension of IDAT increased to 4.
6c>> 2010-03-05 DILUP  Krogh  Added more checks for NDEG too big.
7c>> 2010-03-03 DILUP  Krogh  Fixed bug in flagging extrpolation.
8c>> 2006-05-11 DILUP  Krogh  No flag of extrapolation on equality.
9c>> 2006-04-01 DILUP  Krogh  Fixed minor bug just introduced.
10c>> 2006-03-31 DILUP  Krogh  Don't flag extrap. when on end point.
11c>> 2003-04-17 DILUP  Krogh  Fixed bugs in derivatives with vectors.
12c>> 2000-12-03 DILUP  Krogh  Fixed so no reference to YT when LEXIT=0.
13c>> 2000-12-01 DILUP  Krogh  Removed unused parameter MENTXT.
14c>> 1998-01-26 DILUP  Krogh  Extrap. due to too many bad pts. flagged.
15c>> 1997-09-30 DILUP  Krogh  Decremented NDEGI when LUP=4 & GETERR.
16c>> 1996-04-08 DILUP  Krogh  With LUP=4 & GETERR, NDEGI was 1 too small.
17c>> 1996-03-30 DILUP  Krogh  Added external statement.
18C>> 1995-12-01 DILUP  Krogh  Fixed bugs connected with option 6.
19C>> 1995-11-10 DILUP  Krogh  Fixed so char. data at col. 72 is not ' '.
20C>> 1995-03-03 DILUP  Krogh  Bug in last change made DILUPM fail.
21C>> 1994-12-22 DILUP  Krogh  Added Hermite interpolation.
22C>> 1994-11-11 DILUP  Krogh  Declared all vars.
23c>> 1994-10-20 DILUP  Krogh  Changes to use M77CON
24c>> 1994-09-12 DILUP  Krogh  Added CHGTYP code.
25c>> 1993-04-28 DILUP  Krogh  Additions for Conversion to C.
26c>> 1992-05-27 DILUP  Krogh  Fixed bug in error estimate.
27c>> 1992-04-08 DILUP  Krogh  Removed unused labels 510 and 2080.
28c>> 1991-10-17 DILUP  Krogh  Initial Code.
29c
30c--D replaces "?": ?ILUP, ?ILUPM, C?ILUP, ?MESS
31c
32c Polynomial Interpolation with look up.
33c Design/Code by  Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA
34c
35c In addition to doing standard 1 dimensional interpolation, this
36c subroutine supports efficient interpolation of several functions
37c defined at the same values of the independent variable and supports
38c DILUPM, a subroutine for doing multidimensional interpolation.  Error
39c estimates, Hermite interpolation, and derivatives of the interpolant
40c can be obtained via options.
41c Algorithms used are described in "Efficient Algorithms for Polynomial
42c Interpolation and Numerical Differentiation", by Fred T. Krogh, Math.
43c of Comp. Vol. 24, #109 (Jan. 1970), pp. 185-190.
44c
45c     *************     Formal Arguments     ***************************
46c
47c X      Independent variable where value of interpolant is desired.
48c Y      Value of interpolant, computed by this subroutine.  The
49c        interpolant is always a piecewise polynomial.
50c NTAB   Number of points in the table.
51c XT     Array of independent variable values.  Must be monotone
52c        increasing or monotone decreasing.  If XT(I) = XT(I+1) for some
53c        I, then the XT(J)'s that are used in the interpolation will
54c        either have all J's .le. I, or all J's .ge. I+1.  If the XT's
55c        are equally spaced an option allows one to provide only XT(1),
56c        and the increment between the XT's.
57c YT     Array of dependent variable values.  Y(XT(I)) = YT(I).
58c NDEG   If NDEG < 2 or odd, then it gives the degree of the polynomial
59c        used.  Else the polynomial used is a linear combination of two
60c        polynomials of degree NDEG, such that the resulting polynomial
61c        is of degree NDEG+1 and has a continuous first derivative.
62c        Typically accuracy will improve as NDEG increases up to some
63c        point where it will start to get worse because of either
64c        rounding errors or the inherant instability of high degree
65c        polynoimial interpolation. If more than MAXDEG-th degree is
66c        desired, parameter MAXDEG must be changed below.  It is
67c        currently 15.  If X is so close to the end of the table that
68c        the same number of points can not be selected on both sides of
69c        x the degree of the interpolant is NDEG exactly.  When
70c        extrapolating, the degree used is max(2, 2*(NDEG/2)), where the
71c        divide is truncated to the nearest integer.
72c LUP    Defines the type of look up method.  (Changed if LUP < 1.)
73c    < 0   Use a sequential search starting with an index of -LUP, and
74c          set LUP to -k on exit where k minimizes abs(X-XT(k)).
75c    = 0   As for < 0, except start with a binary search.
76c    = 1   Use a binary search.
77c    = 2   Start a sequential search with an index = [1.5 +
78c          (NTAB-1) * (X-XT(1)) / (XT(NTAB)-XT(1))].
79c    = 3   YT(k) corresponds to XT(1) + (k-1)*XT(2), no search is needed
80c          to do the look up and XT can have dimension 2.
81c    = 4   Internal information connected with X-XT(k) values used in
82c          in the last interpolation is reused.  (Only use if there are
83c          no intervening calls to DILUP.)  No options should be
84c          specified and only YT should be different.  Intended for
85c          interpolating components after the first of a vector valued
86c          function.
87c IOPT   IOPT(1) is used to return a status as follows.
88c    -10   An option index is out of range.
89c    -9    NTAB is outside of allowed limits.
90c    -8    NDEG is outside of allowed limits.
91c    -7    LUP > 3 (use old parameters), used when not ready for it.
92c    -6    Option 3 (compute derivatives), has requested more than
93c          MAXDEG derivatives.
94c    -5    LUP = 3 and XT(2) = 0.
95c    -4    XT(1) = XT(NTAB), and NTAB is not 1.
96c    -3    < 2 points available because of bad points.
97c    -2    Only one table entry available, req. err. est. not computed.
98c    -1    The accuracy requested was not obtained.
99c     0    Nothing special to flag.
100c     1    X was outside the domain of the table, extrapolation used.
101c     2    NTAB is so small, it restricted the degree of the polynomial.
102c
103c        Starting with IOPT(2) options are specified by integers in the
104c        range 0-6, followed in some cases by integers providing
105c        argument(s).  Each option, together with its options its
106c        arguments if any is followed in IOPT by the next option (or 0).
107c     0    No more options; this must be last in the option list.
108c     1    An error estimate is to be returned in EOPT(1).
109c     2    (Argument: K2)  K2 gives the polynomial degree to use when
110c          extrapolating.
111c     3    (Argument: K3, L3) Save (k-th derivative of interpolating
112c          polynomial) / (k!) in EOPT(K3+K-1) for k = 1, 2, ..., L3.
113c          These values are the coefficients of the polynomial in the
114c          monomial basis expanded about X.  One must have 0<L3<NDEG+1.
115c     4    (Argument K4) The absolute and relative errors expected in
116c          YT entries are specified in EOPT(K4) and EOPT(K4+1)
117c          respectively.  The values provided here are used in
118c          estimating the error in the interpolation which is stored
119c          in EOPT(1).
120c     5    (Argument K5, L5) Do the interpolation to the accuracy
121c          requested by the absolute error tolerance specified in
122c          EOPT(K5) and the relative error tolerance in EOPT(K5+1)
123c          respectively.  An attempt is made to keep the final error
124c          < EOPT(K5) + EOPT(K5+1) * (abs(YT(i1) +abs(YT(i2)), where
125c          i1 and i2 are indices for table values close to X.
126c          Interpolation is done as for the default case, except that
127c          in this case NDEG gives the maximal degree polynomial to use
128c          in the interpolation, and polynomial interpolation of even
129c          degree can happen.  (The form that gives the C1 continuity
130c          is not available in this case, and in fact continuity of the
131c          interpolant itself is not to be expected.  The actual degree
132c          used in doing the interpolation is stored in the space for
133c          the argument L5.  An error estimate is returned in EOPT(1).
134c     6    (Argument K6) Do not use point k in the interpolation if
135c          YT(k) = EOPT(K6).
136c     7    (Argument K7) YT(K7+i) gives the first derivative
137c          corresponding to the function value in YT(i).  These
138c          derivatives are to be used in doing the interpolation.  One
139c          gets a continuous interpolant only for NDEG = 3, 7, 11, and
140c          15.  The interpolating polynomial satisfies p(XT(i)) = YT(i),
141c          p'(XT(i)) = YT(K7+i) for values of i that give values of XT
142c          close to X.  (Subject to keeping within one of the same
143c          number of XT's on either side of X, the maximum of |XT(i)-X|
144c          is minimized.)  If NDEG is even a value of YT(i) is used
145c          without using the corresponding value of YT(K7+i).
146c  >253    Used for calls from DILUPM, a number of variables in the
147c          common block have been set.  If this is used, it must be the
148c          first option, and common variables must be set.  The rest
149c          of IOPT is not examined for options.
150c EOPT   Array used to return an error estimate and also used for
151c        options.
152c     EOPT(1)  if an error estimate is returned, this contains a crude
153c              estimate of the error in the interpolation.
154c
155c      ************     External Procedures      ***********************
156c
157c D1MACH Returns system parameters.
158c DMESS  Prints error messages.
159c
160c      ************     Variables referenced     ***********************
161c
162c ADJSAV Saved difference of XT values used to adjust error estimate on
163c  variable order when derivatives are being computed.
164c BADPT  In common CDILUP.  YT(K) is not to be used if YT(K) = BADPT,
165c        and LBADPT = .true.
166c CY     Array giving the YT values used to construct the interpolant.
167c  Also used to store divided differences.
168c D1MACH Function to get parameters of floating point arithmetic.
169c DMESS  Program calling MESS to print error messages.
170c DX     Array giving the differences X - XT(I), where I runs through
171c  the points that are used in the interpolation.
172c DXS    Saved value of DX when computing derivatives
173c E1     Error estimate from one iteration ago (for variable order)
174c E2     Error estimate from this iteration (for variable order)
175c E2L    Used in computing error estimate for variable order.
176c EBND   If estimated error is < EBND then order is suff. high.
177c EBNDI  Internal value for bounding error.
178c EBNDR  Used in getting part of EBND due to relative accuracy request.
179c EF     Factor used in estimating errors for variable order.
180c EN     Used in computing error estimate for variable order.
181c EOPT   Formal argument, see above.
182c EPSR   Relative error level (D1MACH(4)).
183c ERRDAT Holds floating point data for error messages.
184c ERREST Estimate of the error in the interpolation.
185c GETERR Logical variable that is .true. if we are getting an error est.
186c H      In indexed lookup contains the difference between XT values.
187c I      Temporary index.
188c IO     Index used in processing options.
189c IDX    In common CDILUP.  Gives indices of points selected for the
190c        interpolation in order selected.
191c IIFLG  Internal value for IOPT(1).
192c ILI    Lower (sometimes upper) bound on the XT indices that will be
193c  used in the interpolation.
194c INDXED Logical variable that is .true. if we have an indexed XT.
195c IOPT   Formal argument, see above.
196c IUI    As for ILI, except an upper (sometimes lower) bound.
197c K      Index usually into YT and XT, into CY for derivatives.
198c KAOS   In common CDILUP.  Keeps track of state on variable order.
199c    = 0   No variable order -- this value only set in DILUPM.
200c    = 1   First time
201c    = 2   Second time
202c    = 3   Just had an increase or a very weak decrease in the error.
203c    = 4   Just had a strong decrease in the error.
204c   On an error with IIFLG .lt. 0, this is set to -10 * stop level -
205c   the print level.
206c KDER   0 if not doing Hermite interpolation, else is < 0 if have next
207c  higher order difference computed (in TP3), > 0 if it isn't.
208c KEXTRP In common CDILUP.  Gives degree to use when extrapolating.
209c KGO    Indicates type of interpolation as follows:
210c    1     Standard polynomial interpolation.
211c    2     Get error estimate after 1.
212c    3     Compute an interpolant with some desired accuracy.
213c    4     Compute an interpolant with a continuous derivative.
214c  >10     Same as if KGO were 10 smaller, except XT values have already
215c          been selected.
216c KGOC   In common CDILUP.  Value saved for -KGO.  If YT contains values
217c        of Y computed elsewhere in the order specified in IDX, then
218c        KGOC should be set to abs(KGOC) before calling DILUP.  If -2,
219c        an extra value was computed for an error estimate.
220c KK     Used in selecting data points to use.
221c     1  decrease ILI only
222c     0  increase IUI
223c    -1  decrease ILI
224c   <-1, increase IUI only
225c L      Temporary index used in searching the XT array.  Also one less
226c  than the number of points used in the current interpolant.
227c LBADPT Logical variable set = .true. if checking for bad points.
228c LDER   Offset of y' values from the y values.
229c LDERIV In common CDILUP. Loc. where first deriv. is stored in EOPT().
230c LEXERR In common CDILUP. Loc. where absolute and relative error
231c        information on Y is stored in EOPT.
232c LEXIT  In common CDILUP.  Defines action after finding location in XT.
233c    0     Don't compute anything for Y, just return. (Used for DILUPM.)
234c    1     Usual case, just interpolate.
235c   >1     Compute LEXIT-1 derivatives of the interpolant.
236c LINC   Logical variable used in the sequential search.  Set .true. if
237c  the values in XT are increasing with I, and set .false. otherwise.
238c LNDEG  Location in IOPT to save degree when variable order is used.
239c LOPT   Index of the last option.
240c LUP    Formal argument, see above.
241c MACT   Array used for error message actions, see DMESS for details.
242c MACTV  Array used for part of error message for IOPT.
243c MAXDEG Parameter giving the maximum degree polynomial interpolation
244c  supported.  MAXDEG must be odd and > 2.
245c MESS   Message program called from DMESS to print error messages.
246c MEEMES Parameter giving value for printing an error message in MESS.
247c MEIVEC Parameter giving value for printing an integer vector in MESS.
248c MEMDA1 Parameter giving value for indicating value in MACT for MESS.
249c MEMDA2 Parameter giving value for indicating value in MACT for MESS.
250c MEMDAT Parameter giving value for specifying MACT index for MESS.
251c MERET  Parameter giving value to indicate no more actions for MESS.
252c METEXT Parameter giving value to tell MESS to print from MTXTAA.
253c LTXTxx Parameter names of this form were generated by PMESS in making
254c        up text for error messages and are used to locate various parts
255c        of the message.
256c MLOC   Array giving starting loctions for error message text.
257c MTXTAA Character array holding error message text for MESS.
258c N      Used in logic for deciding bounds.  If N = 0, ILI can not be
259c  reduced further.  If N = 3*NTABI, IUI can not be increased further.
260c NDEG   Formal argument, see above.
261c NDEGE  Degree of polynomial to be used.  Starts out = NDEG.
262c NDEGEC In common CDILUP.  Degree actually used, saved in common.
263c NDEGI  Degree of polynomial up to which y & differences are computed.
264c NDEGQ  Used when KGO > 10.  In this case If L .ge. NDEGQ special
265c        action is needed.  (Usually means quitting.)
266c NTAB   Formal argument, see above.
267c NTABI  Internal value of NTAB, = number of points in XT and YT.
268c PI     Array use to store interpolation coefficients
269c PID    Temporary storage used when computing derivatives.
270c SAVED  Logical variable set = .true. if a DX value needs to be
271c  replaced and restored when computing derivatives.
272c TP1    Used for temporary accumulation of values and temp. storage.
273c TP2    Used for temporary storage.
274c TP3    Used for divided difference when doing Hermite interpolation.
275c X      Formal argument, see above.
276c XI     Internal value for X, the place where interpolation is desired.
277c XL     Value of "left hand" X when selecting indexed points.
278c XT     Formal argument, see above.
279c XU     Value of "right hand" X when selecting indexed points.
280c Y      Formal argument, see above.
281c YL     Last value of Y when selecting the order.
282c YT     Formal argument, see above.
283c YTORD  Logical variable set .true. when YT values are ordered on entry
284c
285c     *************     Formal Variable Declarations     ***************
286c
287      external         D1MACH
288      integer          IOPT(*), LUP, NDEG, NTAB
289      double precision X, XT(*), Y, YT(*), EOPT(*)
290c
291c     *************     Common Block and Parameter     *****************
292c
293c                               MAXDEG must be odd and > 2.
294      integer MAXDEG
295      parameter (MAXDEG = 15)
296c
297      logical          GETERR
298      integer          KAOS, KEXTRP, KGOC,LDERIV,LEXERR,LEXIT,NDEGEC,
299     1                 IDX(0:MAXDEG+1)
300      double precision BADPT, DX(0:MAXDEG+1), EBND, EBNDR
301      common / CDILUP / BADPT, DX, EBND, EBNDR, KAOS, KEXTRP, KGOC,
302     1                  LDERIV, LEXERR,LEXIT,NDEGEC,IDX,GETERR
303c
304c     *************     Local Variables     ****************************
305c
306      logical          INDXED, LBADPT, LINC, YTORD, SAVED
307      integer          I, IIFLG, ILI, IUI, IO, K, KDER, KGO, KK, L,
308     1                LDER, LNDEG, LOPT, N, NDEGE, NDEGI, NDEGQ, NTABI
309      double precision ADJSAV, CY(0:MAXDEG+1), DXS, E1, E2, E2L, EBNDI,
310     1                 EF, EM, EPSR, ERREST, H, PI(0:MAXDEG+1),
311     2                 PID(0:MAXDEG), TP1, TP2, TP3, XI, XL, XU, YL
312      double precision D1MACH
313      save             EPSR, /CDILUP/
314c
315c ************************ Error Message Stuff and Data ****************
316c
317c Parameter defined below are all defined in the error message program
318c DMESS.
319c
320      integer MECONT,MEMDA1,MERET,MEEMES,METEXT,MEIVEC
321      parameter (MEMDA1 =27)
322      parameter (MECONT = 50)
323      parameter (MERET =51)
324      parameter (MEEMES =52)
325      parameter (METEXT =53)
326      parameter (MEIVEC =57)
327c
328      double precision ERRDAT(2)
329      integer MLOC(9), IDAT(4), MACT(5), MACTV(6)
330
331c ********* Error message text ***************
332c[Last 2 letters of Param. name]  [Text generating message.]
333cAA DILUP$B
334cAB Estimated error = $F, Requested error = $F.$E
335cAC NTAB = 1, no error estimate computed.$E
336cAD Too many bad points, only $I point(s) available.$E
337cAE XT(1) = XT(NTAB=$I) = $F ??$E
338cAF LUP = 3, and XT(2) = 0.$E
339cAG Order of requested derivative, $I, is > bound of $I.$E
340cAH LUP > 3, with unprepared common block.$E
341cAI NDEG = $I must be .le. $I or  NTAB = $I must be at $C
342c   least $I.  You must either decrease NDEG, increase NTAB, $C
343c   or if requesting an error estimate, turn off that request.$E
344cAJ NDEG = $I is not in the allowed interval of [0, $I].$E
345cAK NTAB = $I is not in the allowed interval of [1, $I].$E
346cAL IOPT($I) = $I is not a valid option.$E
347c   $
348cAM IOPT(1:$M): $E
349      integer LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,LTXTAG,LTXTAH,
350     * LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM
351      parameter (LTXTAA=  1,LTXTAB=  8,LTXTAC= 53,LTXTAD= 92,LTXTAE=142,
352     * LTXTAF=171,LTXTAG=196,LTXTAH=250,LTXTAI=290,LTXTAJ=458,
353     * LTXTAK=512,LTXTAL=566,LTXTAM=  1)
354      character MTXTAA(3) * (201)
355      character MTXTAB(1) * (14)
356      data MTXTAA/'DILUP$BEstimated error = $F, Requested error = $F.$EN
357     *TAB = 1, no error estimate computed.$EToo many bad points, only $I
358     * point(s) available.$EXT(1) = XT(NTAB=$I) = $F ??$ELUP = 3, and XT
359     *(2) = 0.$EOrder ','of requested derivative, $I, is > bound of $I.$
360     *ELUP > 3, with unprepared common block.$ENDEG = $I must be .le. $I
361     * or  NTAB = $I must be at least $I.  You must either decrease NDEG
362     *, increase NTAB, or if',' requesting an error estimate, turn off t
363     *hat request.$ENDEG = $I is not in the allowed interval of [0, $I].
364     *$ENTAB = $I is not in the allowed interval of [1, $I].$EIOPT($I) =
365     * $I is not a valid option.$E'/
366      data MTXTAB/'IOPT(1:$M): $E'/
367c **** End of automatically generated text
368c
369c                      1 2 3 4      5
370      data MACT / MEEMES,0,0,0, MECONT /
371c                        2        3      4   5     6
372      data MACTV /MEMDA1, 0, METEXT, MEIVEC, 0, MERET /
373
374c
375      data MLOC / LTXTAB, LTXTAC, LTXTAD, LTXTAE, LTXTAF, LTXTAG,
376     1            LTXTAH, LTXTAI, LTXTAJ /
377c
378      data EPSR / 0.D0 /
379c
380c
381c      ************     Start of executable code     *******************
382c
383      IIFLG = 0
384      KDER = 0
385      PI(0) = 1.D0
386      LBADPT = .false.
387      ILI = -LUP
388      if (ILI .le. -4) go to 1400
389      LEXERR = 0
390      NTABI = NTAB
391      if ((NTABI .le. 0) .or. (NTABI .gt. 9999999)) then
392         IIFLG = -9
393         IDAT(1) = NTABI
394         IDAT(2) = 9999999
395         go to 2120
396      end if
397      NDEGI = NDEG
398      NDEGE = NDEGI
399      XI = X
400      KGO = 1
401      IO = 1
402      if (IOPT(2) .ge. 254) then
403         LBADPT = IOPT(2) .eq. 255
404         if (KAOS .le. 0) go to 50
405         LNDEG = 0
406         KGO = 3
407         KAOS = 1
408         NDEGI = min(MAXDEG, IOPT(3)+1)
409         NDEGE = -NDEGI
410         go to 60
411      end if
412      GETERR = .false.
413      KEXTRP = -1
414      LEXIT = 1
415c                                      Loop to take care of options
416   10 IO = IO + 1
417      LOPT = IOPT(IO)
418      if (LOPT .ne. 0) then
419        go to (1930, 1500, 1600, 1900, 1950, 1970, 1980), LOPT
420        IIFLG = -10
421        IDAT(1) = IO
422        IDAT(2) = LOPT
423        go to 2120
424      end if
425   50 continue
426      K = NDEGI + 1
427      if (KDER .ne. 0) then
428        K = K / 2
429      else if (mod(NDEGI,2) .eq. 0) then
430        K = K + 1
431      end if
432      if (GETERR) K = K + 1
433      if ((NDEGI.lt.0) .or. (NDEGI.gt.MAXDEG) .or. (K.gt.NTAB)) then
434         IIFLG = -8
435         IDAT(1) = NDEGI
436         IDAT(2) = MAXDEG
437         IDAT(3) = NTAB
438         IDAT(4) = K
439         go to 2120
440      end if
441      if (KGO .eq. 1) then
442         if (NDEGI .ge. 2) then
443            if (mod(NDEGI, 2) .eq. 0) then
444               if (KDER .eq. 0) then
445                  NDEGE = NDEGE + 1
446                  KGO = 4
447               end if
448            end if
449         end if
450      end if
451   60 if (KEXTRP .lt. 0) then
452         KEXTRP = max(NDEGE-1, min(NDEGE, 2))
453         if (KEXTRP .lt. 0) KEXTRP = NDEGI
454      end if
455      if (ILI .gt. 0) go to 200
456      if (ILI + 2) 1300, 1200, 100
457c
458c                                      Binary search, then sequential
459  100    continue
460c                         In binary search XT(ILI) .le. XI .le.  XT(IUI)
461c                         or we are extrapolating
462            ILI = 1
463            IUI = NTABI
464            if (XT(NTABI) - XT(1)) 110, 1220, 130
465  110       ILI = NTABI
466            L = 1
467  120       IUI = L
468  130       L = (IUI - ILI) / 2
469            if (L .eq. 0) go to 210
470            L = ILI + L
471            if (XT(L) .gt. XI) go to 120
472            ILI = L
473            go to 130
474c
475c                                      Sequential search, then exit
476  200    continue
477            if (ILI .gt. NTABI) go to 100
478            if (XT(NTABI) .eq. XT(1)) go to 1220
479  210       INDXED = .false.
480            LINC = XT(NTABI) .gt. XT(1)
481            if ((XT(ILI) .gt. XI) .eqv. LINC) go to 230
482  220       if (ILI .eq. NTABI) go to 1230
483            ILI = ILI + 1
484            if ((XT(ILI) .lt. XI) .eqv. LINC) go to 220
485            N = 2*ILI
486            if (abs(XT(ILI-1)-XI) .lt. abs(XT(ILI) - XI)) then
487               ILI = ILI - 1
488               N = N - 1
489            end if
490            go to 240
491  230       if (ILI .eq. 1) go to 1230
492            ILI = ILI - 1
493            if ((XT(ILI) .gt. XI) .eqv. LINC) go to 230
494            N = 2*ILI + 1
495            if (abs(XT(ILI+1)-XI) .lt. abs(XT(ILI) - XI)) then
496               ILI = ILI + 1
497               N = N + 1
498            end if
499  240       if (LUP .le. 0) LUP = -ILI
500  250       DX(0) = XI - XT(ILI)
501c                                  Get bounding indices and interpolate
502  260       IUI = ILI
503            K = ILI
504            L = -1
505            if (LEXIT .eq. 0) then
506               NDEGI = 0
507               if (GETERR) then
508                  if (KGO .ne. 3) then
509                     if (KGO .eq. 1) NDEGE = NDEGE + 1
510                     KEXTRP = min(KEXTRP+1, NDEGE)
511                  else
512                     NDEGE = -NDEGE
513                  end if
514               end if
515            end if
516c                                  Just got index for next XT
517  270       if (LEXIT .eq. 0) then
518               L = L + 1
519               IDX(L) = K
520               go to 350
521            end if
522            TP2 = YT(K)
523            if (LBADPT) then
524c                                  Check if point should be discarded
525               if (TP2 .eq. BADPT) then
526                  if (IUI .ne. ILI) then
527                     if (abs(KK) .eq. 1) then
528                        N = N - 1
529                     else
530                        N = N + 1
531                     end if
532                  end if
533                  go to 1020
534               end if
535            end if
536  280       L = L + 1
537            IDX(L) = K
538c
539  290       if (KDER .ne. 0) then
540c                           Hermite interpolation
541               KDER = -KDER
542               if (KDER .lt. 0) then
543                  TP3 = YT(LDER + K)
544                  if (L .eq. 0) then
545                     TP1 = TP2
546                  else
547                     do 300 I = 1, L
548                        TP2 = (TP2 - CY(I-1)) / (DX(I-1) - DX(L))
549                        TP3 = (TP3 - TP2) / (DX(I-1) - DX(L))
550  300                continue
551                     PI(L) = PI(L-1) * DX(L-1)
552                     TP1 = TP1 + PI(L) * TP2
553                  end if
554               else
555                  DX(L) = DX(L-1)
556                  PI(L) = PI(L-1) * DX(L)
557                  TP2 = TP3
558                  TP1 = TP1 + PI(L) * TP2
559               end if
560            else if (L .eq. 0) then
561                  TP1 = TP2
562            else
563               PI(L) = PI(L-1) * DX(L-1)
564               if (L .le. NDEGI) then
565c                                 Get divided differences & interpolate.
566                  do 320 I = 1, L
567                     TP2 = (TP2 - CY(I-1)) / (DX(I-1) - DX(L))
568  320             continue
569                  TP1 = TP1 + PI(L) * TP2
570               end if
571            end if
572            CY(L) = TP2
573  350       if (L .lt. NDEGE) go to 1000
574  360       if (LEXIT .eq. 0) go to 2110
575            go to (500, 600, 900, 800), KGO
576            if (L .ge. NDEGQ) go to (500,600,900,800), KGO-10
577c                        Already got the points selected.
578  400       L = L + 1
579            if (YTORD) then
580               TP2 = YT(L+1)
581            else
582               K = IDX(L)
583               TP2 = YT(K)
584            end if
585            go to 290
586c Got simple interpolated value.
587  500       Y = TP1
588            if (.not. GETERR) go to  2000
589            NDEGI = NDEGI + 1
590            KGO = KGO + 1
591            if (NDEGE .ge. 0) go to 1000
592            go to 400
593c Got info. for error estimate.
594  600       if (L .eq. 0) then
595               IIFLG = -2
596            else
597               ERREST=1.5D0*(abs(TP1-Y)+.03125D0*abs(PI(L-1)*CY(L-1)))
598               L = L - 1
599            end if
600            go to  2000
601c C1 interpolant.
602  800       do 810 I = 1, L-1
603               TP2 = (TP2 - CY(I-1)) / (DX(I-1) - DX(L))
604  810       continue
605            TP2 = (TP2 - CY(L-1)) / (DX(0) - DX(1))
606            CY(L) = TP2
607            PI(L) = PI(L-1) * DX(0)
608            Y = TP1 + PI(L) * TP2
609            if (GETERR) ERREST = 1.5D0 * abs(PI(L-1)) * abs(((DX(L-1) *
610     1         (DX(0) - DX(1)) / (DX(L-1) - DX(L)) - DX(0)) * TP2) +
611     2         abs(.03125D0*CY(L-1)))
612            go to 2000
613c    Variable order, check estimated error and convergence.
614  900       continue
615            if (KAOS .ge. 3) go to 930
616            if (KAOS .eq. 2) go to 920
617c                                        First time
618            if ((KDER.ne.0) .and. (LEXIT.gt.1) .and. (L.eq.0)) go to 970
619            KAOS = 2
620            E2L = abs(TP1)
621            E2 = EBND + 1.D30
622            go to 970
623c                                        Second time
624  920       KAOS = 4
625            EBNDI = .66666D0*(EBND + EBNDR*(abs(CY(0))+abs(YT(IDX(1)))))
626            EF = DX(0)
627            if (LEXIT .gt. 1) then
628               EF = DX(L) - DX(0)
629               ADJSAV = EF
630            end if
631            E2 = abs(EF * TP2)
632            EM = .75D0
633            go to 950
634c                                        Usual case
635  930       E2L = E2
636            E1 = E2L * (5.D0 * EM / dble(L))
637            EF = EF * DX(L-1)
638            E2 = abs(EF*TP2)
639            EM = 0.5D0 * EM + E2 / (E2L + E2 + 1.D-20)
640            if (E2 .ge. E1) then
641               if (KAOS .eq. 3) then
642c                          Apparently diverging, so quit.
643                  L = L - 1
644                  go to 960
645               else
646                  KAOS = 3
647               end if
648            else
649               KAOS = 4
650            end if
651  950       YL = TP1
652            if (E2L + E2 .gt. EBNDI) then
653               if ((L+NDEGE .lt. 0) .and. (IIFLG .ne. 2)) go to 970
654               if (KGO .eq. 13) then
655c                                 May need early exit to get next point.
656                  IOPT(2) = 0
657                  if (L .lt. NDEG) go to 2110
658               end if
659            end if
660  960       TP2 = 1.5D0
661            if (LEXIT .gt. 1) then
662               if (KDER .eq. 0) TP2 = 1.5D0 * abs(DX(0) / ADJSAV)
663            end if
664            ERREST = TP2 * (E2  + .0625D0 * E2L)
665            if (LNDEG .ne. 0) IOPT(LNDEG) = L
666            Y = YL
667            go to 2000
668  970       if (KGO .gt. 10) go to 400
669c
670 1000       if (KDER .lt. 0) go to 280
671 1020       KK = min(N - IUI - ILI, 2) - 1
672c In this section of code, KK=: 1, decrease ILI only; 0, increase IUI;
673c -1, decrease ILI; and <-1, increase IUI only
674            if (abs(KK) .eq. 1) then
675               ILI = ILI - 1
676               K = ILI
677               if (ILI .ne. 0) then
678                  if (INDXED) then
679                     XL = XL + H
680                     DX(L+1) = XL
681                     go to 270
682                  end if
683                  DX(L+1) = XI - XT(K)
684                  if (XT(ILI+1) .ne. XT(ILI)) go to 270
685               end if
686               if (KK .eq. 1) go to 1100
687               N = 0
688c                If all points will be on one side, flag extrapolation.
689c               if (L .le. 0) go to 1250
690            else
691               IUI = IUI + 1
692               K = IUI
693               if (IUI .le. NTABI) then
694                  if (INDXED) then
695                     XU = XU - H
696                     DX(L+1) = XU
697                     go to 270
698                  end if
699                  DX(L+1) = XI - XT(K)
700                  if (XT(IUI-1) .ne. XT(IUI)) go to 270
701               end if
702               if (KK .ne. 0) go to 1100
703               N = 3*NTABI
704c                If all points will be on one side, flag extarpolation.
705c               if (L .le. 0) go to 1250
706            end if
707            if (KGO .lt. 4) go to 1020
708c If we can't get the same number of points on either side, the
709c continuous derivative case becomes standard interpolation.
710            KGO = 1
711            NDEGE = NDEGE - 1
712            if (L .lt. NDEGE) go to 1020
713            go to 360
714c                                     No more data accessible.
715 1100       if (L .le. 0) then
716c                             Too many bad points, couldn't find points.
717                Y = TP1
718                IDAT(1) = L + 1
719                IIFLG = -3
720                go to 2110
721            end if
722            NDEGI = min(NDEGI, L)
723            NDEGE = 0
724            IIFLG = 2
725            go to 360
726c
727c                                     Secant start, then use sequential
728 1200    continue
729            if (XT(1) .eq. XT(NTABI)) go to 1220
730            ILI = max(1, min(NTABI, int(1.5D0+dble(NTABI-1)*(XI-XT(1))/
731     1         (XT(NTABI) - XT(1)))))
732         go to 210
733c
734c                         Special cases
735c                                  1 entry in XT
736 1220    ILI = 1
737         IUI = 1
738         KGO = 1
739         K = 1
740         if (NDEGE .ne. 0) IIFLG = 2
741         NDEGE = 0
742         if (NTABI .eq. 1) go to 250
743c                         Error -- XT(1) .eq. XT(NTAB), and NTAB .ne. 1
744         IIFLG = -4
745         IDAT(1) = NTABI
746         ERRDAT(1) = XT(1)
747         KGO = 0
748         go to 2120
749c                  Check if really extrapolating
750 1230    N = 6 * (ILI - 1)
751         if ((XI - XT(1)) * (XT(NTABI) - XI) .ge. 0.D0) go to 240
752c                                  Extrapolating
753 1250    IIFLG = 1
754         if (KGO .ge. 4) KGO = 1
755         NDEGI = KEXTRP
756         NDEGE = sign(NDEGI, NDEGE)
757         if (INDXED) go to 260
758         go to 240
759c
760c                                      Index search, then exit
761 1300    continue
762         INDXED = .true.
763         H = XT(2)
764         if (H .eq. 0) then
765            IIFLG = -5
766            go to 2120
767         end if
768         TP1 = 1.D0 + (XI - XT(1)) / H
769         ILI = min(max(1, nint(TP1)), NTABI)
770         XU = (TP1 - dble(ILI)) * H
771         XL = XU
772         DX(0) = XU
773         N = ILI + ILI
774         if (TP1 .gt. dble(ILI))  N = N + 1
775         if ((XI - XT(1)) * ((XT(1) + H*dble(NTABI-1) - XI)) .ge. 0.E0)
776     1     go to 260
777         N = 6 * (ILI - 1)
778         go to 1250
779c
780c                                 Already set up
781 1400    KGO = KGOC
782         if (KGO .eq. 0) then
783            IIFLG = -7
784            go to 2120
785         end if
786         YTORD = KGO .gt. 0
787         KGO = abs(KGO)
788         if (KGO .lt. 10) KGO = KGO + 10
789         if (KGO .eq. 12) KGO = 11
790         NDEGI = NDEGEC
791         NDEGQ = NDEGI
792         if (KGO .eq. 13) NDEGQ = 0
793         NDEGE = -NDEGI
794         L = -1
795         if (KGO .eq. 14) NDEGI = NDEGI - 1
796         go to 400
797c
798c                                 Set number of points for extrapolation
799 1500    continue
800            I = I + 1
801            KEXTRP = min(IOPT(I), NDEGE)
802            go to 10
803c
804c                                 Get derivatives of interpolant
805 1600    continue
806            I = I + 2
807            LDERIV = IOPT(I-1)
808            LEXIT = IOPT(I) + 1
809            go to 10
810c
811c                                     Get expected errors in YT.
8121900     continue
813            I = I + 1
814            LEXERR = IOPT(I)
815c
816c                                     Set to get error estimate.
817 1930    continue
818            GETERR = .true.
819            go to 10
820c
821c                                     Set for automatic order selection.
822 1950    continue
823            IO = IO + 2
824            LNDEG = IO
825            EBND = max(0.D0, EOPT(IOPT(IO-1)))
826            EBNDR = max(0.D0, EOPT(IOPT(IO-1)+1))
827            KGO = 3
828            KAOS = 1
829            NDEGI = min(MAXDEG, NDEGE+1)
830            NDEGE = -NDEGI
831         go to 1930
832c
833c                                    Set to ignore special points
834 1970    continue
835            LBADPT = .true.
836            IO = IO + 1
837            BADPT = EOPT(IOPT(IO))
838            go to 10
839c
840c                                    Set up for Hermite interpolation.
841 1980    continue
842         IO = IO + 1
843         LDER = IOPT(IO)
844         KDER = abs(LDER)
845         go to 10
846c
847c                                    Put any new options just above here
848c                    End of the loop
849c
850 2000 if (LEXIT .lt. 2) go to 2100
851c                       Compute derivatives of Y
852      SAVED = mod(KGO, 10) .eq. 4
853      if (SAVED) then
854         DXS = DX(L-1)
855         DX(L-1) = DX(0)
856      end if
857      N = LEXIT - 1
858      if (N .gt. L) then
859        if (N .gt. MAXDEG) then
860           IDAT(1) = N
861           IDAT(2) = MAXDEG
862           N = MAXDEG
863           IIFLG = -6
864        end if
865        do 2020 I = L+1, N
866           EOPT(I+LDERIV-1) = 0.D0
867 2020   continue
868        N = L
869      end if
870      TP1 = CY(1)
871      PID(0) = 1.D0
872      do 2030 I = 1, L-1
873         PID(I) = PI(I) + DX(I) * PID(I-1)
874         TP1 = TP1 + PID(I) * CY(I+1)
875 2030 continue
876      EOPT(LDERIV) = TP1
877      do 2070 K = 2, N
878         TP1 = CY(K)
879         do 2060 I = 1, L-K
880            PID(I) = PID(I) + DX(I+K-1) * PID(I-1)
881            TP1 = TP1 + PID(I) * CY(I+K)
882 2060    continue
883         EOPT(K+LDERIV-1) = TP1
884 2070 continue
885      if (SAVED) DX(L-1) = DXS
886c                                    Save info. and return
887 2100 continue
888      if (GETERR) then
889         if (EPSR .eq. 0.D0) then
890            EPSR = D1MACH(4)
891         end if
892         TP1 = EPSR
893         TP2 = 0.D0
894         if (LEXERR .ne. 0) then
895            TP2 = max(0.D0, EOPT(LEXERR))
896            TP1 = max(TP1, EOPT(LEXERR+1))
897         end if
898         if (IIFLG .eq. 2) ERREST = 32.D0 * ERREST
899         EOPT(1) = ERREST + TP2 + TP1*(abs(CY(0)) + abs(DX(0)*CY(1)))
900         if ((KGO .eq. 3) .or. (KGO .eq. 13)) then
901            if (EBNDI .ne. 0.D0) then
902               if (EOPT(1) .gt. EBND) then
903                  ERRDAT(1) = EOPT(1)
904                  ERRDAT(2) = EBND
905                  IIFLG = -1
906               end if
907            end if
908         end if
909      end if
910 2110 KGOC = -KGO
911      NDEGEC = L
912 2120 IOPT(1) = IIFLG
913      if (IIFLG .ge. 0) return
914c Take care of error messages, IIFLG < -2 should stop.
915      MACT(2) = 88
916      if (IIFLG .lt. -1) KAOS = -MACT(2)
917      if (IIFLG .ge. -3) MACT(2) = min(26, 24 - IIFLG)
918      if (IOPT(2) .ge. 254)  then
919         if (IIFLG .eq. -1) return
920         MACT(2) = 28
921      end if
922      MACT(3) = -IIFLG
923      MACT(4) = MLOC(-IIFLG)
924      call DMESS(MACT, MTXTAA, IDAT, ERRDAT)
925      MACTV(2) = IO
926      MACTV(5) = IO
927      K = 1
928      if (IO .eq. 1) K = 6
929      call MESS(MACTV(K), MTXTAB, IOPT)
930      return
931c                       End of Interpolation routine -- DILUP
932      end
933