subroutine simpin1_init(ixi,tlout,alocalout,blocalout,iord,lbig,x1grid,n1grid)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    simpin1_init
!     prgmmr:    parrish     org: np22                date: 2000-01-01
!
! abstract:  This routine computes all Taylor matrices ever needed
!            for 1-d interpolation on grid x1grid.
!
! program history log:
!   2000-01-01  parrish
!   2004-06-22  treadon - update documentation
!
!   input argument list:
!     iord   - order of interpolation (1=linear, 2=quadratic, etc)
!     lbig   - number of interpolating points (=iord+1)
!     x1grid - coordinates of interpolator grid
!                (can be monotonic increasing or decreasing)
!     n1grid - dimension of interpolator grid
!
!   output argument list:
!     ixi    - local coordinate order counter
!     tl     - output taylor matrices
!     alocal,blocal - used to transform to local coordinates around each grid-point
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use constants, only: zero, half, one, two, one_tenth
  implicit none

! Declare local parameters
  real(r_kind),parameter:: r0_001=0.001_r_kind

  integer(i_kind),intent(in   ) :: iord,lbig,n1grid
  integer(i_kind),intent(  out) :: ixi(0:iord)
  real(r_kind)   ,intent(  out) :: tlout(lbig,lbig,2*n1grid),&
       alocalout(2*n1grid),blocalout(2*n1grid)
  real(r_kind)   ,intent(in   ) :: x1grid(n1grid)

  real(r_kind) x1in(2*n1grid)
  real(r_kind) dx1gridi(-3:n1grid+3)
  real(r_kind) tl(max(64,2*n1grid/lbig),lbig,lbig)
  integer(i_kind) i1ref(2*n1grid)
  integer(i_kind) ix1sign(2*n1grid)
  real(r_kind) x1temp(2*n1grid),x1p(2*n1grid)
  integer(i_kind) iflag(2*n1grid)
  real(r_kind) alocal(max(64,2*n1grid/lbig)),blocal(max(64,2*n1grid/lbig))
  real(r_kind) xmaxlocal(max(64,2*n1grid/lbig)),xminlocal(max(64,2*n1grid/lbig))
  integer(i_kind), allocatable, dimension(:) :: ix1grid
  
  real(r_kind) delx1,dxa,dxb,dxc,dxd,dxe,dxf,dxg,dxfinei,dxmax,rhigh,rlow,x1ref,xlocal
  real(r_kind) xboundright,xboundleft,dxfine,dxminmin,dxmin,dxthis
  integer(i_kind) i1fine,i1ref0,ih,ii,ip,iximx,j,l,lp,n1fine,nend,np,nstart,nthis,ntl
  integer(i_kind) nmaxright,nn,nin,n,iximax,nminleft,iximn,nm,i
  
  nin=2*n1grid
  do n=1,n1grid
     nm=max(1,n-1)
     np=min(n1grid,n+1)
     delx1=x1grid(np)-x1grid(np-1)
     x1in(n)=x1grid(n)+one_tenth*delx1
     delx1=x1grid(nm+1)-x1grid(nm)
     x1in(n1grid+n)=x1grid(n)-one_tenth*delx1
  end do
  
  ntl=max(64,nin/lbig)
  ! set ixi, the coordinate order counter
 
  do i=0,iord
     ih=(i+1)/2
     if(2*ih==i) then
        ixi(i)=-ih
     else
        ixi(i)=ih
     end if
!     write(6,100)i,ixi(i)
!     100 format(' xi(',i2,')=',f10.0)
  end do

  !  find and mark all points that are outside interval of interpolation

  iximx=maxval(ixi) ; iximn=minval(ixi)
  iximax=max(abs(iximx),abs(iximn))
  nminleft=abs(iximn)+1 ; nmaxright=n1grid-iximax
  if(iord==1) then
     xboundleft=x1grid(1)
     xboundright=x1grid(n1grid)
  else
     xboundleft=x1grid(1+iximax)-.49999_r_kind*(x1grid(1+iximax)-x1grid(iximax))
     xboundright=x1grid(n1grid-iximax)+ &
          .49999_r_kind*(x1grid(n1grid-iximax+1)-x1grid(n1grid-iximax))
  end if
  if(x1grid(n1grid)>x1grid(1)) then
     iflag=1
     do n=1,nin
        if(x1in(n)<=xboundleft.or.x1in(n)>=xboundright) iflag(n)=0
     end do
  end if
  if(x1grid(n1grid)<x1grid(1)) then
     iflag=1
     do n=1,nin
        if(x1in(n)>=xboundleft.or.x1in(n)<=xboundright) iflag(n)=0
     end do
  end if

  !  set up uniform fine grid to use in finding interpolation coordinates
 
  dxmax=-huge(dxmax) ; dxmin=huge(dxmin)
  do i=1,n1grid-1
     dxthis=x1grid(i+1)-x1grid(i)
     dx1gridi(i)=one/dxthis
     dxmax=max(dxthis,dxmax) ; dxmin=min(dxthis,dxmin)
  end do
  dx1gridi(-3:0)=dx1gridi(1)
  dx1gridi(n1grid:n1grid+3)=dx1gridi(n1grid-1)
  if(dxmax*dxmin<=zero) then
     write(6,*)' INTERPOLATION GRID NOT MONOTONIC IN SIMPIN1'
     do i=1,n1grid-1
        write(6,*)' i,x1grid,dx=',i,x1grid(i),x1grid(i+1)-x1grid(i)
     end do
     stop
  end if
  dxminmin=min(abs(dxmax),abs(dxmin))
  dxfine=sign(dxminmin,dxmax)
  dxfinei=one/dxfine
  n1fine=1+ceiling((x1grid(n1grid)-x1grid(1))/dxfine)
!  write(6,*)' in simpin1, n1grid,n1fine=',n1grid,n1fine
  allocate (ix1grid(n1fine))
  ii=1
  do i=1,n1fine
     x1ref=x1grid(1)+(i-1)*dxfine
     if(dxfinei*(x1ref-x1grid(ii+1))>=-r0_001) ii=min(ii+1,n1grid-1)
     ix1grid(i)=ii
!     write(6,)' x1fine,x1grid=',x1ref,x1grid(ix1grid(i))
  end do
 
  ! now get i1ref, index of interval containing point
 
  rlow=-epsilon(rlow)
  rhigh=one+epsilon(rhigh)
  do n=1,nin
     i1fine=1+nint((x1in(n)-x1grid(1))*dxfinei)
     i1fine=max(1,min(i1fine,n1fine))
     i1ref0=max(4,min(ix1grid(i1fine),n1grid-4))
     dxa=(x1in(n)-x1grid(i1ref0-3))*dx1gridi(i1ref0-3)
     dxb=(x1in(n)-x1grid(i1ref0-2))*dx1gridi(i1ref0-2)
     dxc=(x1in(n)-x1grid(i1ref0-1))*dx1gridi(i1ref0-1)
     dxd=(x1in(n)-x1grid(i1ref0  ))*dx1gridi(i1ref0  )
     dxe=(x1in(n)-x1grid(i1ref0+1))*dx1gridi(i1ref0+1)
     dxf=(x1in(n)-x1grid(i1ref0+2))*dx1gridi(i1ref0+2)
     dxg=(x1in(n)-x1grid(i1ref0+3))*dx1gridi(i1ref0+3)
     if(dxa<=rhigh) i1ref(n)=i1ref0-3
     if(dxb>=rlow.and.dxb<=rhigh) i1ref(n)=i1ref0-2
     if(dxc>=rlow.and.dxc<=rhigh) i1ref(n)=i1ref0-1
     if(dxd>=rlow.and.dxd<=rhigh) i1ref(n)=i1ref0
     if(dxe>=rlow.and.dxe<=rhigh) i1ref(n)=i1ref0+1
     if(dxf>=rlow.and.dxf<=rhigh) i1ref(n)=i1ref0+2
     if(dxg>=rlow) i1ref(n)=i1ref0+3
     i1ref(n)=max(nminleft,min(i1ref(n),nmaxright))
  end do
  deallocate (ix1grid)

  !             i1ref now has index of interval containing point
  !               adjust to be closest one to interpolating point
  ix1sign=1
  do n=1,nin
     dxa=(x1in(n)-x1grid(i1ref(n)))*dx1gridi(i1ref(n))
     if(dxa>half) then
        i1ref(n)=min(i1ref(n)+1,nmaxright)
        ix1sign(n)=-1
     end if
  end do

  ! get taylor matrices and invert
  
  nstart=1
  nend=min(nin,ntl)
  do while (nstart<=nend)
     nthis=nend-nstart+1
     
!  compute limits of local x coordinate

     xmaxlocal=-huge(xmaxlocal) ; xminlocal=huge(xminlocal)
     do i=0,iord
        do n=nstart,nend
           if(iflag(n)>0) then
              nn=n-nstart+1
              xlocal=x1grid(i1ref(n)+ix1sign(n)*ixi(i))-x1grid(i1ref(n))
              xmaxlocal(nn)=max(xmaxlocal(nn),xlocal)
              xminlocal(nn)=min(xminlocal(nn),xlocal)
           end if
        end do
     end do
     alocal=zero
     blocal=zero
     do n=nstart,nend
        if(iflag(n)>0) then 
           nn=n-nstart+1
           alocal(nn)=two/(xmaxlocal(nn)-xminlocal(nn))
           blocal(nn)=-one-two*xminlocal(nn)/(xmaxlocal(nn)-xminlocal(nn))
           alocalout(n)=alocal(nn)
           blocalout(n)=blocal(nn)
        end if
     end do
     
     tl=zero
     do i=1,lbig
        tl(1:nthis,i,i)=one
     end do
     l=0
     do i=0,iord
        l=l+1
        x1temp(nstart:nend)=zero
        do n=nstart,nend
           nn=n-nstart+1
           x1temp(n)=alocal(nn)*(x1grid(i1ref(n)+ix1sign(n)*ixi(i))-x1grid(i1ref(n))) &
                +blocal(nn)
        end do
        lp=0
        x1p(nstart:nend)=one
        do ip=0,iord
           lp=lp+1
           do n=nstart,nend
              if(iflag(n)>0) tl(n-nstart+1,l,lp)=x1p(n)
           end do
           x1p(nstart:nend)=x1p(nstart:nend)*x1temp(nstart:nend)
        end do
     end do
     
     call vinvmm(tl,tl,lbig,lbig,lbig,nthis,ntl)
     
     do n=nstart,nend
        nn=n-nstart+1
        do j=1,lbig
           do i=1,lbig
              tlout(i,j,n)=tl(nn,i,j)
           end do
        end do
     end do
     
     nstart=nstart+ntl
     nend=min(nend+ntl,nin)
     
  end do
  
  return
end subroutine simpin1_init