      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, z-direction
c
c     Input file: physical parameters and derivatives in nodes.
c                 Fastest dimension, z-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 z-directionm
c                                   n3=4 cross derivatives. (dx dz)
c------------------------------------------------------------------------
c%     
      implicit NONE
      integer maxnod,maxn3

      parameter (maxnod = 2000000, maxn3 = 3*maxnod)

      real*4 velpar(maxnod),veldpar(maxn3),eps/1.0e-15/
      
      integer n1,n2,n3,esize,numn
      real*4  d1,d2,o1,o2
      integer fetch
      integer numd1,numd2,n3n,inum,ixdir,izdir,ixn1,ixn2
      real*4 rdelx,rdelz,r2dx,r2dz,rdxdz

      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=4

c     write updated trace length to new header

      n3n = 4
      call putch('n1','d',n1)
      call putch('n2','d',n2)
      call putch('n3','d',n3n)
      call putch('esize','d',esize)
      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

      call sreed('in',velpar,n1*n2*esize)

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

      if (maxn3.lt.(n1*n2*3)) 
     *     call erexit('max output trace length exceeded ')
      
      rdelx = 0.5/d2
      r2dx = 2.0*rdelx
      rdelz = 0.5/d1
      r2dz = 2.0*rdelz
      rdxdz = rdelx*rdelz
      numd1 = n1*n2
      numd2 = numd1*2
      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))*r2dz
      veldpar(numd2+1) = 
     *         (velpar(n1+2)-velpar(2)-velpar(n1+1)+velpar(1))*4.0*rdxdz

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))*rdelz
         veldpar(numd2+inum) = (velpar(n1+inum+1)-velpar(inum+1)
     *                     -velpar(n1+inum-1)+velpar(inum-1))*2.0*rdxdz
      enddo

c     evaluation for lower left corner.

      veldpar(n1) = (velpar(2*n1)-velpar(n1))*r2dx
      veldpar(numd1+n1) = (velpar(n1) - velpar(n1-1))*r2dz
      veldpar(numd2+n1) = (velpar(2*n1)-velpar(n1)
     *                     -velpar(2*n1-1)+velpar(n1-1))*4.0*rdxdz

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))*r2dz
         veldpar(numd2+ixn1+1) = (velpar(ixn1+n1+2)
     *                                -velpar(ixn1-n1+2)
     *                                -velpar(ixn1+n1+1)
     *                                +velpar(ixn1-n1+1))*2.0*rdxdz
         
c     evaluation of all nodes  except for the nodes on the edges.

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

            veldpar(ixn2) = (velpar(ixn2+n1)
     *                             -velpar(ixn2-n1))*rdelx
            veldpar(numd1+ixn2) = (velpar(ixn2+1) 
     *                             - velpar(ixn2-1))*rdelz
            veldpar(numd2+ixn2) = (velpar(ixn2+n1+1)
     *                            -velpar(ixn2-n1+1)
     *                            -velpar(ixn2+n1-1)
     *                            +velpar(ixn2-n1-1))*rdxdz
         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))*r2dz
         veldpar(numd2+ixn1+n1) = (velpar(ixn1+2*n1)
     *                               -velpar(ixn1)
     *                               -velpar(ixn1+2*n1-1)
     *                               +velpar(ixn1-1))*2.0*rdxdz

      enddo


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

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))*rdelz
         veldpar(numd2+numn+inum) = (velpar(numn+inum+1)
     *                               -velpar(numn+inum-n1+1)
     *                               -velpar(numn+inum-1)
     *                               +velpar(numn+inum-n1-1))*2.0*rdxdz
      enddo

c     evaluation for lower right corner.

      veldpar(numd1) = (velpar(numd1)-velpar(numd1-n1))*r2dx
      veldpar(numd1+n1) = (velpar(numd1) - velpar(numd1-1))*r2dz
      veldpar(numd2+n1) = (velpar(numd1)
     *                     -velpar(numd1-n1)
     *                     -velpar(numd1-1)
     *                     +velpar(numd1-n1-1))*4.0*rdxdz

c     write to file
      
      call srite('out',velpar,numd1*esize)
      call srite('out',veldpar,numd1*3*esize)

      end











