Subroutine lagrng (x,arg,y,val,ndim,nfs,nptx,maxarg,maxfs)
c *** see RH Landau, Program lpott, 1981
c     lagrange interpolation,unequally spaced points
c     npts=2,3,4,5,6, nfs functions(y*s) simultaneously interpltd
c     x= value of argument, arg is the tabulated x*s
c     y= a vector of interpolad functions, from tabulted val
c     ndim= Dimension of table
c     nfs= = of functions simult interpolated
c     maxarg and maxfs are maximum values of subscripts,ie Dimensions
      Implicit  Real*8   (a-h,o-z)
      Dimension arg(maxarg), val(maxfs,maxarg), y(maxfs)
c     -----find x0 the Closest point to x
      npts = nptx
c                                                                   2=st
      ni = 1
      nf = ndim
 10   If ((x.le.arg(ni)).or.(x.ge.arg(nf))) GoTo 30
      If ((nf-ni+1).eq.2) GoTo 70
      nmid = (nf+ni)/2
      If (x.gt.arg(nmid)) GoTo 20
      nf = nmid
      GoTo 10
 20   ni = nmid
      GoTo 10
c     ------ x is one of the tabltd values
 30   If (x.le.arg(ni)) GoTo 60
      nn = nf
 40   nused = 0
      Do 50 n=1,nfs
         y(n) = val(n,nn)
 50   Continue
      Return
 60   nn = ni
      GoTo 40
c     ------- 2 pts left, chpose smaller one
 70   n0 = ni
      nn = npts-1
      GoTo (140,110,100,90,80), nn
 80   Continue
      If (((n0+3).gt.ndim).or.((n0-2).lt.1)) GoTo 90
      nused = 6
      GoTo 150
 90   Continue
      If ((n0+2).gt.ndim) GoTo 110
      If ((n0-2).lt.1) GoTo 100
      nused = 5
      GoTo 150
 100  Continue
      If (((n0+2).gt.ndim).or.((n0-1).lt.1)) GoTo 110
      nused = 4
      GoTo 150
 110  If ((n0+1).lt.ndim) GoTo 130
c     n0=ndim,special case
 120  nn = ndim
      GoTo 40
 130  nused = 3
      If ((n0-1).lt.1) nused = 2
      GoTo 150
 140  If ((n0+1).gt.ndim) GoTo 120
      nused = 2
 150  Continue
c     at least 2 pts left
      y0 = x-arg(n0)
      y1 = x-arg(n0+1)
      y01 = y1-y0
      c0 = y1/y01
      c1 = -y0/y01
      If (nused.eq.2) GoTo 160
c     at least 3 pts
      ym1 = x-arg(n0-1)
      y0m1 = ym1-y0
      ym11 = y1-ym1
      cm1 = -y0*y1/y0m1/ym11
      c0 = c0*ym1/y0m1
      c1 = -c1*ym1/ym11
      If (nused.eq.3) GoTo 180
c     ------at least 4 pts
      y2 = x-arg(n0+2)
      ym12 = y2-ym1
      y02 = y2-y0
      y12 = y2-y1
      cm1 = cm1*y2/ym12
      c0 = c0*y2/y02
      c1 = c1*y2/y12
      c2 = -ym1*y0*y1/ym12/y02/y12
      If (nused.eq.4) GoTo 200
c     at least 5 pts
      ym2 = x-arg(n0-2)
      ym2m1 = ym1-ym2
      ym20 = y0-ym2
      ym21 = y1-ym2
      ym22 = y2-ym2
      cm2 = ym1*y0*y1*y2/ym2m1/ym20/ym21/ym22
      cm1 = -cm1*ym2/ym2m1
      c0 = -c0*ym2/ym20
      c1 = -c1*ym2/ym21
      c2 = -c2*ym2/ym22
      If (nused.eq.5) GoTo 220
c     at least 6 pts
      y3 = x-arg(n0+3)
      ym23 = y3-ym2
      ym13 = y3-ym1
      y03 = y3-y0
      y13 = y3-y1
      y23 = y3-y2
      cm2 = cm2*y3/ym23
      cm1 = cm1*y3/ym13
      c0 = c0*y3/y03
      c1 = c1*y3/y13
      c2 = c2*y3/y23
      c3 = ym2*ym1*y0*y1*y2/ym23/ym13/y03/y13/y23
      GoTo 240
 160  Continue
      Do 170 n=1,nfs
         y(n) = c0*val(n,n0)+c1*val(n,n0+1)
 170  Continue
      GoTo 260
 180  Continue
      Do 190 n=1,nfs
         y(n) = cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)
 190  Continue
      GoTo 260
 200  Continue
      Do 210 n=1,nfs
         y(n) = cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)+c2*val(n,n0+
     1   2)
 210  Continue
      GoTo 260
 220  Continue
      Do 230 n=1,nfs
         y(n) = cm2*val(n,n0-2)+cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0
     1   +1)+c2*val(n,n0+2)
 230  Continue
      GoTo 260
 240  Continue
      Do 250 n=1,nfs
         y(n) = cm2*val(n,n0-2)+cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0
     1   +1)+c2*val(n,n0+2)+c3*val(n,n0+3)
 250  Continue
 260  Return
      End