c/* %W%      latest revision %G% %U% */
      subroutine lagrng(x,arg,y,val,ndim,nfs,nptx,maxarg,maxfs)
************************************************************************

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


c >>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      npts=nptx
      ni=1
      nf=ndim
   10 if((x .le. arg(ni)) .or.  (x .ge. arg(nf)))go to 30
      if((nf-ni+1) .eq. 2) go to 40
      nmid=(nf+ni)/2
      if( x .gt. arg(nmid)) go to 20
      nf=nmid
      go to 10
   20 ni=nmid
      go to 10
c------ x is one of the tabltd values
   30 if( x .le. arg(ni)) go to 35
      nn=nf
   31 nused=0
      do 32 n=1,nfs
   32 y(n)=val(n,nn)
      return
   35 nn=ni
      go to 31
c------- 2 pts left, chpose smaller one
   40 n0=ni
   50 nn=npts-1
      go to (200,300,400,500,600),nn
  600 continue
      if(((n0+3) .gt. ndim) .or. ((n0-2) .lt. 1)) go to 500
      nused=6
      go to 1000
  500 continue
      if((n0+2) .gt. ndim) go to 300
      if((n0-2) .lt. 1) go to 400
      nused=5
      go to 1000
  400 continue
      if(((n0+2).gt. ndim) .or. ((n0-1) .lt. 1)) go to 300
      nused=4
      go to 1000
  300 if((n0+1) .lt. ndim) go to 305
c    n0=ndim,special case
  302 nn=ndim
      go to 31
  305 nused=3
      if((n0-1) .lt. 1) nused=2
      go to 1000
  200 if((n0+1) .gt. ndim) go to 302
      nused=2
 1000 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) go to 2000
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) go to 3000
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) go to 4000
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) go to 5000
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
      go to 6000
 2000 continue
      do 2100 n=1,nfs
 2100  y(n)=c0*val(n,n0)+c1*val(n,n0+1)
      go to 7000
 3000 continue
      do 3100 n=1,nfs
 3100 y(n)= cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)
      go to 7000
 4000 continue
      do 4100 n=1,nfs
 4100 y(n)=cm1*val(n,n0-1)+c0*val(n,n0)+c1*val(n,n0+1)+c2*val(n,n0+2)
      go to 7000
 5000 continue
      do 5100 n=1,nfs
 5100 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)
      go to 7000
 6000 continue
      do 6100 n=1,nfs
 6100 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)
 7000 return
      end