      program Splintab2
c-----------------------------------------------
c
c     Splintab < in > out 
c
c-----------------------------------------------------------------
c     Purpose: Make values of 1st derivatives, and 2nd cross derivatives
c              at nodes in a 2-D gridded model.  The values are input
c              to bcucof.f, described in "Numerical Recipies in FORTRAN,
c              which computes coefficient needed in bicubic interpolation
c--------------------------------------------------------------------------
c     Input file: physical parameters in nodes. fastest dimension, tau-direction
c
c     Input file: physical parameters and derivatives in nodes.
c                 Fastest dimension, tau-direction
c                 Sekond dimension, x-direction
c                 Third direction:  n3=1 parameters as in input 
c                                   n3=2 1st derivatives in x-directionm
c                                   n3=3 1st derivatives in tau-directionm
c                                   n3=4 cross derivatives. (dx dtau)
c                                   n3=5 sigma (as in input)
c                                   n3=6 dsigma/ dx
c                                   n3=7 dsigma/ dtau
c                                   n3=8 sigma cross derivatives. (dx dtau)
c------------------------------------------------------------------------
c%     
      implicit NONE
      integer maxnod,maxn6

      parameter (maxnod = 6000000, maxn6 = 6*maxnod)

      real*4 velpar(maxnod),veldpar(maxn6),eps/1.0e-15/
      real*4 sigpar(maxnod)
      complex input(maxnod)
      
      integer n1,n2,n3,esize,esizeout,numn
      real*4  d1,d2,o1,o2
      integer fetch,sreed
      integer numd1,numd2,numd3,numd4,numd5,numd6,numd7
      integer n3n,inum,ixdir,itaudir,ixn1,ixn2
      real*4 rdelx,rdeltau,r2dx,r2dtau,rdxdtau

      call initpar

      call doc("Splintab.f")

c     fetch dimensions and esize.

      if (fetch('n1','d',n1).eq.0) call erexit('n1 must be supplied')
      if (fetch('n2','d',n2).eq.0) call erexit('n2 must be supplied')
      if (fetch('n3','d',n3).gt.1) n3=1
      if (fetch('d1','f',d1).lt.eps) call erexit('d1 must be supplied')
      if (fetch('d2','f',d2).lt.eps) call erexit('d2 must be supplied')
      if (fetch('o1','f',o1).lt.eps) o1 = 0.0
      if (fetch('o2','f',o2).lt.eps) o2 = 0.0

      if (fetch('esize','d',esize).eq.0) esize=8
      esizeout=4

      if (n1*n2 .gt. maxnod) call seperr('n1*n2 > maxnod')

c     write updated trace length to new header

      n3n = 8
      call putch('n1','d',n1)
      call putch('n2','d',n2)
      call putch('n3','d',n3n)
      call putch('esize','d',esizeout)
      call putch('d1','f',d1)
      call putch('d2','f',d2)
      call putch('o1','f',o1)
      call putch('o2','f',o2)

c     close history file

      call hclose()

c     read input file

      if( sreed('in',input,n1*n2*esize) .ne. n1*n2*esize) then
           call seperr("Error in reading input")
      end if
      do inum=1,n1*n2
            velpar(inum)=real(input(inum))
            sigpar(inum)=imag(input(inum))
      end do

      if (maxnod.lt.(n1*n2)) 
     *     call erexit('max input trace length exceeded ')

      if (maxn6.lt.(n1*n2*5)) 
     *     call erexit('max output trace length exceeded ')
      
      rdelx = 0.5/d2
      r2dx = 2.0*rdelx
      rdeltau = 0.5/d1
      r2dtau = 2.0*rdeltau
      rdxdtau = rdelx*rdeltau
      numd1 = n1*n2
      numd2 = numd1*2
      numd3 = numd1*3
      numd4 = numd1*4
      numd5 = numd1*5
      numd6 = numd1*6
      numd7 = numd1*7
      numn = numd1-n1+1

c     evaluation for upper left corner.
      
      veldpar(1) = (velpar(n1+1)-velpar(1))*r2dx
      veldpar(numd1+1) = (velpar(2)-velpar(1))*r2dtau
      veldpar(numd2+1) = 
     *       (velpar(n1+2)-velpar(2)-velpar(n1+1)+velpar(1))*4.0*rdxdtau

      veldpar(numd3+1) = (sigpar(n1+1)-sigpar(1))*r2dx
      veldpar(numd4+1) = (sigpar(2)-sigpar(1))*r2dtau
      veldpar(numd5+1) = 
     *       (sigpar(n1+2)-sigpar(2)-sigpar(n1+1)+sigpar(1))*4.0*rdxdtau

c     evaluation for left column except upper and lower corners.

      do inum=2,n1-1
         veldpar(inum) = (velpar(n1+inum)-velpar(inum))*r2dx
         veldpar(numd1+inum) = (velpar(inum+1) - velpar(inum-1))*rdeltau
         veldpar(numd2+inum) = (velpar(n1+inum+1)-velpar(inum+1)
     *                    -velpar(n1+inum-1)+velpar(inum-1))*2.0*rdxdtau

         veldpar(numd3+inum) = (sigpar(n1+inum)-sigpar(inum))*r2dx
         veldpar(numd4+inum) = (sigpar(inum+1) - sigpar(inum-1))*rdeltau
         veldpar(numd5+inum) = (sigpar(n1+inum+1)-sigpar(inum+1)
     *                    -sigpar(n1+inum-1)+sigpar(inum-1))*2.0*rdxdtau

      enddo

c     evaluation for lower left corner.

      veldpar(n1) = (velpar(2*n1)-velpar(n1))*r2dx
      veldpar(numd1+n1) = (velpar(n1) - velpar(n1-1))*r2dtau
      veldpar(numd2+n1) = (velpar(2*n1)-velpar(n1)
     *                     -velpar(2*n1-1)+velpar(n1-1))*4.0*rdxdtau

      veldpar(numd3+n1) = (sigpar(2*n1)-sigpar(n1))*r2dx
      veldpar(numd4+n1) = (sigpar(n1) - sigpar(n1-1))*r2dtau
      veldpar(numd5+n1) = (sigpar(2*n1)-sigpar(n1)
     *                     -sigpar(2*n1-1)+sigpar(n1-1))*4.0*rdxdtau

c     loop over columns except for left and right columns.

      do ixdir =2,n2-1
         ixn1 = (ixdir-1)*n1
         
c     evaluation of upper node on column no. ixdir
         
         veldpar(ixn1+1) = (velpar(ixn1+n1+1)
     *                                -velpar(ixn1-n1+1))*rdelx
         veldpar(numd1+ixn1+1) = (velpar(ixn1+2) 
     *                                -velpar(ixn1+1))*r2dtau
         veldpar(numd2+ixn1+1) = (velpar(ixn1+n1+2)
     *                                -velpar(ixn1-n1+2)
     *                                -velpar(ixn1+n1+1)
     *                                +velpar(ixn1-n1+1))*2.0*rdxdtau

         veldpar(numd3+ixn1+1) = (sigpar(ixn1+n1+1)
     *                                -sigpar(ixn1-n1+1))*rdelx
         veldpar(numd4+ixn1+1) = (sigpar(ixn1+2) 
     *                                -sigpar(ixn1+1))*r2dtau
         veldpar(numd5+ixn1+1) = (sigpar(ixn1+n1+2)
     *                                -sigpar(ixn1-n1+2)
     *                                -sigpar(ixn1+n1+1)
     *                                +sigpar(ixn1-n1+1))*2.0*rdxdtau
         
c     evaluation of all nodes  except for the nodes on the edges.

         do itaudir =2,n1-1
            ixn2 = ixn1+itaudir

            veldpar(ixn2) = (velpar(ixn2+n1)
     *                             -velpar(ixn2-n1))*rdelx
            veldpar(numd1+ixn2) = (velpar(ixn2+1) 
     *                             - velpar(ixn2-1))*rdeltau
            veldpar(numd2+ixn2) = (velpar(ixn2+n1+1)
     *                            -velpar(ixn2-n1+1)
     *                            -velpar(ixn2+n1-1)
     *                            +velpar(ixn2-n1-1))*rdxdtau

            veldpar(numd3+ixn2) = (sigpar(ixn2+n1)
     *                             -sigpar(ixn2-n1))*rdelx
            veldpar(numd4+ixn2) = (sigpar(ixn2+1) 
     *                             - sigpar(ixn2-1))*rdeltau
            veldpar(numd5+ixn2) = (sigpar(ixn2+n1+1)
     *                            -sigpar(ixn2-n1+1)
     *                            -sigpar(ixn2+n1-1)
     *                            +sigpar(ixn2-n1-1))*rdxdtau

         enddo

c     evaluation of lower node on column no. ixdir
         
         veldpar(ixn1+n1) = (velpar(ixn1+2*n1)
     *                             -velpar(ixn1))*rdelx
         veldpar(numd1+ixn1+n1) = (velpar(ixn1+n1) 
     *                               -velpar(ixn1+n1-1))*r2dtau
         veldpar(numd2+ixn1+n1) = (velpar(ixn1+2*n1)
     *                               -velpar(ixn1)
     *                               -velpar(ixn1+2*n1-1)
     *                               +velpar(ixn1-1))*2.0*rdxdtau

         veldpar(numd3+ixn1+n1) = (sigpar(ixn1+2*n1)
     *                             -sigpar(ixn1))*rdelx
         veldpar(numd4+ixn1+n1) = (sigpar(ixn1+n1) 
     *                               -sigpar(ixn1+n1-1))*r2dtau
         veldpar(numd5+ixn1+n1) = (sigpar(ixn1+2*n1)
     *                               -sigpar(ixn1)
     *                               -sigpar(ixn1+2*n1-1)
     *                               +sigpar(ixn1-1))*2.0*rdxdtau

      enddo


c     evaluation for upper right corner.
      
      veldpar(numn) = (velpar(numn)
     *                       -velpar(numn-n1))*r2dx
      veldpar(numd1+numn) = (velpar(numn+1)
     *                       -velpar(numn))*r2dtau
      veldpar(numd2+numn) = (velpar(numn+1)
     *                             -velpar(numn-n1+1)
     *                             -velpar(numn)
     *                             +velpar(numn-n1))*4.0*rdxdtau

      veldpar(numd3+numn) = (sigpar(numn)
     *                       -sigpar(numn-n1))*r2dx
      veldpar(numd4+numn) = (sigpar(numn+1)
     *                       -sigpar(numn))*r2dtau
      veldpar(numd5+numn) = (sigpar(numn+1)
     *                             -sigpar(numn-n1+1)
     *                             -sigpar(numn)
     *                             +sigpar(numn-n1))*4.0*rdxdtau

c     evaluation for right column except upper and lower corners.
      
      do inum=2,n1-1
         veldpar(numn+inum) = (velpar(numn+inum)
     *                              -velpar(numn+inum-1))*r2dx
         veldpar(numd1+numn+inum) = (velpar(numn+inum+1) 
     *                              -velpar(numn+inum-1))*rdeltau
         veldpar(numd2+numn+inum) = (velpar(numn+inum+1)
     *                              -velpar(numn+inum-n1+1)
     *                              -velpar(numn+inum-1)
     *                              +velpar(numn+inum-n1-1))*2.0*rdxdtau

         veldpar(numd3+numn+inum) = (sigpar(numn+inum)
     *                              -sigpar(numn+inum-1))*r2dx
         veldpar(numd4+numn+inum) = (sigpar(numn+inum+1) 
     *                              -sigpar(numn+inum-1))*rdeltau
         veldpar(numd5+numn+inum) = (sigpar(numn+inum+1)
     *                              -sigpar(numn+inum-n1+1)
     *                              -sigpar(numn+inum-1)
     *                              +sigpar(numn+inum-n1-1))*2.0*rdxdtau
      enddo

c     evaluation for lower right corner.

      veldpar(numd1) = (velpar(numd1)-velpar(numd1-n1))*r2dx
      veldpar(numd1+numd1) = (velpar(numd1) - velpar(numd1-1))*r2dtau
      veldpar(numd2+numd1) = (velpar(numd1)
     *                     -velpar(numd1-n1)
     *                     -velpar(numd1-1)
     *                     +velpar(numd1-n1-1))*4.0*rdxdtau

      veldpar(numd3+numd1) = (sigpar(numd1)-sigpar(numd1-n1))*r2dx
      veldpar(numd4+numd1) = (sigpar(numd1) - sigpar(numd1-1))*r2dtau
      veldpar(numd5+numd1) = (sigpar(numd1)
     *                     -sigpar(numd1-n1)
     *                     -sigpar(numd1-1)
     *                     +sigpar(numd1-n1-1))*4.0*rdxdtau

c     write to file
      
      call srite('out',velpar,numd1*esizeout)
      call srite('out',veldpar(1),numd1*3*esizeout)
      call srite('out',sigpar,numd1*esizeout)
      call srite('out',veldpar(numd3+1),numd1*3*esizeout)

      end

