c
c *** Seismic Waves are computed at arbitrary points on the receiver plane
c     using the bicubic spline interpolation
c
      parameter (mv=11,mh=101,mt=2050,mob=1010)
      dimension xv(mv),xh(mh)
      dimension x(mt,mob),y(mt,mob),z(mt,mob),xxv(mob),xxh(mob)
      dimension ggx(mh,mv,mt),ggy(mh,mv,mt),ggz(mh,mv,mt)
      dimension gx(mh,mv),gy(mh,mv),gz(mh,mv)
     *         ,gx2(mh,mv),gy2(mh,mv),gz2(mh,mv)
      dimension gxtmp(mh),gytmp(mh),gztmp(mh)
      dimension gx2tmp(mh),gy2tmp(mh),gz2tmp(mh)

C
C * Read data from 'grflt8.in' *
C
      PI = ACOS(-1.0)
      write(*,*)'Data from "grflt8.in"'

      OPEN (7,FILE='grflt8.in',STATUS='OLD')
      read(7,1100)dt,nt
      read(7,1110)period
 1100 format(/f8.0,i8)
 1110 format(f8.0)
      duration=dt*nt
      if(nt.gt.mt)then
        write(*,*)'nt > mt !'
        stop
      end if
      write(*,*)
      write(*,*)'Delta Time (s) and Number of Times =',dt,nt
      write(*,*)'Minimum Period to be computed (s)  =',period
C
      READ(7,1200)NL
c      write(*,*)'nl=',nl
      DO 100 IL=1,NL
        READ(7,*)
  100 CONTINUE
 1200 FORMAT(/I10)
C
      READ(7,1300)NF
c      write(*,*)'nf=',nf
 1300 FORMAT(/2I10)
      DO 250 I=1,NF
        READ(7,*)
        READ(7,*)
        READ(7,*)
        READ(7,*)
        READ(7,*)
  250 CONTINUE
C
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
C
      READ(7,1700)ISIGN
 1700 FORMAT(/I10)
c
      read(7,1500)
 1500 format(//)
      read(7,1510)nobh,olnh
      read(7,1510)nobv
 1510 format(i10,f10.0)
      if(nobh.gt.mh)then
        write(*,*)'nobh > mh !'
        stop
      end if
      if(nobv.gt.mv)then
        write(*,*)'nobv > mv !'
        stop
      end if
      write(*,*)
c      write(*,*)'See Locations of Vertical Nodes in fort.11'
c      write(11,*)'Location of Vertical Nodes (m)'
      xvmax=0.0
      DO 260 I=1,nobv
        read(7,1550)xv(I)
        if(xv(i).gt.xvmax)xvmax=xv(i)
c        write(11,*)'i, xv(i)=',i,xv(i)
 260  CONTINUE
 1550 format(f10.0) 
      write(*,*)
      write(*,*)'Number of Receivers along V and H axes =',NOBV,NOBH
      write(*,*)'Horizontal Length of the Plane (m)=',olnh
      write(*,*)'Vertical Length of the Plane (m)  =',xvmax

      close(7)
c
c   * Configure the horizontal axis *
c      write(*,*)
c      write(*,*)'See Locations of Horizontal Nodes in fort.11'
c      write(11,*)'Location of Horizontal Nodes (m)'
      if(nobh.eq.1)then
        xh(1)=0.0
      else
        dh=olnh/(nobh-1)
        do 270 i=1,nobh
          xh(i)=dh*(i-1)
c          write(11,*)'i, xh(i)=',i,xh(i)
 270    continue
      end if


c
c *** Read Seismic Waves on Grids on the Plane from "wave_dat" ***
c
      OPEN(12,FILE='wave_dat',STATUS='unknown')
c
      write(*,*)
c      write(*,*)'See Values at All Nodes in fort.11'
c      write(11,*)'ih,iv,it,and Gx,Gy,GZ:'
      do 300 iv=1,nobv
      do 300 ih=1,nobh
      read(12,*)
      do 310 i=1,nt
      read(12,4000)ggx(ih,iv,i),ggy(ih,iv,i),ggz(ih,iv,i)
c      write(*,4100)ih,iv,i,ggx(ih,iv,i),ggy(ih,iv,i),ggz(ih,iv,i)
 310  continue
 300  continue
 4000 format(3e12.4)
 4100 format(3i4,3e12.4)
c
      close(12)
c
c *** Read number of arbitrary observation points ***
c
      open(7,file='splin8.in',status='old')
      read(7,*)nob
      write(*,*)
      write(*,*)' Number of Receivers =',nob
c      write(*,*)
c      write(*,*)'See Location of Receivers in fort.11'
c      write(11,*)'Location of Receivers;'
c      write(11,*)'io,xxv,xxh'
      if(nob.gt.mob)then
        write(*,*)'nob > mob !'
        stop
      end if
      do io=1,nob
      read(7,*)xxv(io),xxh(io)
      if((nobv.ne.1).and.(xxv(io).gt.xvmax))then
        write(*,*)'xxv > xvmax !!',xxv(io)
        stop
      end if
      if((nobh.ne.1).and.(xxh(io).gt.olnh))then
        write(*,*)'xxh > olnh !!',xxh(io)
        stop
      end if
      if((nobv.eq.1).and.(xxv(io).ne.0.0))then
        write(*,*)'nobv=1, but xxv.ne.0 !!'
c        stop
      end if
      if((nobh.eq.1).and.(xxh(io).ne.0.0))then
        write(*,*)'nobh=1, but xxh.ne.0 !!'
c        stop
      end if
c      write(11,*)'io,xxv,xxh',io,xxv(io),xxh(io)
      end do
      close(7)
c
c *** Spline Interpolation ***
c
      do 400 it=1,nt
c
        do 410 iv=1,nobv
        do 410 ih=1,nobh
          gx(ih,iv)=ggx(ih,iv,it)
          gy(ih,iv)=ggy(ih,iv,it)
          gz(ih,iv)=ggz(ih,iv,it)
 410    continue
c
        if(nobv.eq.1)then

          if(nobh.eq.1)then
            x(it,1)=gx(1,1)
            y(it,1)=gy(1,1)
            z(it,1)=gz(1,1)
          else
            do ih=1,nobh
            gxtmp(ih)=gx(ih,1)
            gytmp(ih)=gy(ih,1)
            gztmp(ih)=gz(ih,1)
            call spline(xh,gxtmp,gytmp,gztmp,nobh,1.e30,1.e30
     *                 ,gx2tmp,gy2tmp,gz2tmp)
            do io=1,nob
              call splint(xh,gxtmp,gytmp,gztmp,gx2tmp,gy2tmp
     *        ,gz2tmp,nobh,xxh(io),x(it,io),y(it,io),z(it,io))
            end do
            end do
          end if

        else

          if(nobh.eq.1)then

            do iv=1,nobv
            gxtmp(iv)=gx(1,iv)
            gytmp(iv)=gy(1,iv)
            gztmp(iv)=gz(1,iv)
            call spline(xv,gxtmp,gytmp,gztmp,nobv,1.e30,1.e30
     *                 ,gx2tmp,gy2tmp,gz2tmp)
            do io=1,nob
              call splint(xv,gxtmp,gytmp,gztmp,gx2tmp,gy2tmp
     *        ,gz2tmp,nobv,xxv(io),x(it,io),y(it,io),z(it,io))
            end do
            end do

          else
c
c         * Calculate 2nd derivatives "y2" using a natural spline *
            call splie2(xh,xv,gx,gy,gz,nobh,nobv,mh,mv
     *                 ,gx2,gy2,gz2)
c         * Calculate seismic wave at arbitrary observation points *
            do io=1,nob
              call splin2(xh,xv,gx,gy,gz,gx2,gy2,gz2,nobh,nobv
     *         ,mh,mv,xxh(io),xxv(io),x(it,io),y(it,io),z(it,io))
            end do

          end if

        end if

 400  continue
c
c *** Output ***
c
      open(8,file='out_dat',status='unknown')
      do 500 io=1,nob
      write(8,5000)io
      do 500 it=1,nt
      write(8,5010)x(it,io),y(it,io),z(it,io)
 500  continue
 5000 format(i5)
 5010 format(3e13.4)
      close(8)

      write(*,*)
      write(*,*)'Output for Matlab? (Yes=1)'
      read(*,*)mat
      if(mat.eq.1)then
        open(8,file='mat0_dat',status='unknown')
        write(8,*)nob,nt,dt
        close(8)
        open(7,file='matx_dat',status='unknown')
        open(8,file='maty_dat',status='unknown')
        open(9,file='matz_dat',status='unknown')
        do 600 io=1,nob
        do 600 it=1,nt
          write(7,5010)x(it,io)
          write(8,5010)y(it,io)
          write(9,5010)z(it,io)
 600    continue
        close(7)
        close(8)
        close(9)
      end if
c
      stop
      end
c
c
c *** Bicubic Spline Interpolaion ***
c
      subroutine splie2(xh,xv,gx,gy,gz,m,n,mh,mv,gx2,gy2,gz2)
      parameter(nn=22)
      dimension xh(mh),xv(mv),gx(mh,mv),gy(mh,mv),gz(mh,mv)
     *         ,gx2(mh,mv),gy2(mh,mv),gz2(mh,mv)
     *         ,gxtmp(nn),gytmp(nn),gztmp(nn)
     *         ,gx2tmp(nn),gy2tmp(nn),gz2tmp(nn)

      if(n.gt.nn)then
        write(*,*)'From splie2: n > nn !'
        stop
      end if

      if(m.gt.mh)then
        write(*,*)'From splie2: m > mh !'
        stop
      end if

      if(n.gt.mv)then
        write(*,*)'From splie2: n > mv !'
        stop
      end if

      do 13 j=1,m
        do 11 k=1,n
          gxtmp(k)=gx(j,k)
          gytmp(k)=gy(j,k)
          gztmp(k)=gz(j,k)
 11     continue
        call spline(xv,gxtmp,gytmp,gztmp
     *             ,n,1.e30,1.e30,gx2tmp,gy2tmp,gz2tmp)
        do 12 k=1,n
          gx2(j,k)=gx2tmp(k)
          gy2(j,k)=gy2tmp(k)
          gz2(j,k)=gz2tmp(k)
 12     continue
 13   continue

      return
      end
c
c
      subroutine splin2(xh,xv,gx,gy,gz,gx2,gy2,gz2
     *                 ,m,n,mh,mv,x1,x2,x,y,z)
      parameter (nn=22)
      dimension xh(mh),xv(mv),gx(mh,mv),gy(mh,mv),gz(mh,mv)
     *         ,gx2(mh,mv),gy2(mh,mv),gz2(mh,mv)
     *         ,gxtmp(nn),gytmp(nn),gztmp(nn)
     *         ,gx2tmp(nn),gy2tmp(nn),gz2tmp(nn)
     *         ,xtmp(nn),ytmp(nn),ztmp(nn)

      if(n.gt.nn)then
        write(*,*)'From splin2: n > nn !'
        stop
      end if

      if(m.gt.mh)then
        write(*,*)'From splin2: m > mh !'
        stop
      end if

      if(n.gt.mv)then
        write(*,*)'From splin2: n > mv !'
        stop
      end if

      do 12 j=1,m
        do 11 k=1,n
          gxtmp(k)=gx(j,k)
          gytmp(k)=gy(j,k)
          gztmp(k)=gz(j,k)
          gx2tmp(k)=gx2(j,k)
          gy2tmp(k)=gy2(j,k)
          gz2tmp(k)=gz2(j,k)
 11     continue
        call splint(xv,gxtmp,gytmp,gztmp,gx2tmp,gy2tmp,gz2tmp
     *             ,n,x2,xtmp(j),ytmp(j),ztmp(j))
 12   continue
      call spline(xh,xtmp,ytmp,ztmp
     *           ,m,1.e30,1.e30,gx2tmp,gy2tmp,gz2tmp)
      call splint(xh,xtmp,ytmp,ztmp
     *           ,gx2tmp,gy2tmp,gz2tmp,m,x1,x,y,z)

      return
      end
c
c 
c *** Cubic Spline Interpolation ***
c
      subroutine spline(x,gx,gy,gz,n,yp1,ypn,gx2,gy2,gz2)
      parameter (nmax=1000)
      dimension x(n),gx(n),gy(n),gz(n),gx2(n),gy2(n),gz2(n)
     *         ,ua(nmax),ub(nmax),uc(nmax)

      if(n.gt.nmax)then
        write(*,*)'From spline: n > nmax !'
        stop
      end if

      if(yp1.gt..99e30)then
        gx2(1)=0.
        gy2(1)=0.
        gz2(1)=0.
        ua(1)=0.
        ub(1)=0.
        uc(1)=0.
      else
        gx2(1)=-0.5
        gy2(1)=-0.5
        gz2(1)=-0.5
        ua(1)=(3./(x(2)-x(1)))*((gx(2)-gx(1))/(x(2)-x(1))-yp1)
        ub(1)=(3./(x(2)-x(1)))*((gy(2)-gy(1))/(x(2)-x(1))-yp1)
        uc(1)=(3./(x(2)-x(1)))*((gz(2)-gz(1))/(x(2)-x(1))-yp1)
      end if

      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        pa=sig*gx2(i-1)+2.
        pb=sig*gy2(i-1)+2.
        pc=sig*gz2(i-1)+2.
        gx2(i)=(sig-1.)/pa
        gy2(i)=(sig-1.)/pb
        gz2(i)=(sig-1.)/pc
        ua(i)=(6.*((gx(i+1)-gx(i))/(x(i+1)-x(i))-(gx(i)-gx(i-1))
     *       /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*ua(i-1))/pa
        ub(i)=(6.*((gy(i+1)-gy(i))/(x(i+1)-x(i))-(gy(i)-gy(i-1))
     *       /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*ub(i-1))/pb
        uc(i)=(6.*((gz(i+1)-gz(i))/(x(i+1)-x(i))-(gz(i)-gz(i-1))
     *       /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*uc(i-1))/pc
 11   continue

      if(ypn.gt..99e30)then
        qna=0.
        qnb=0.
        qnc=0.
        una=0.
        unb=0.
        unc=0.
      else
        qna=0.5
        qnb=0.5
        qnc=0.5
        una=(3./(x(n)-x(n-1)))*(ypn-(gx(n)-gx(n-1))/(x(n)-x(n-1)))
        unb=(3./(x(n)-x(n-1)))*(ypn-(gy(n)-gy(n-1))/(x(n)-x(n-1)))
        unc=(3./(x(n)-x(n-1)))*(ypn-(gz(n)-gz(n-1))/(x(n)-x(n-1)))
      end if

      gx2(n)=(una-qna*ua(n-1))/(qna*gx2(n-1)+1.)
      gy2(n)=(unb-qnb*ub(n-1))/(qnb*gy2(n-1)+1.)
      gz2(n)=(unc-qnc*uc(n-1))/(qnc*gz2(n-1)+1.)
      do 12 k=n-1,1,-1
        gx2(k)=gx2(k)*gx2(k+1)+ua(k)
        gy2(k)=gy2(k)*gy2(k+1)+ub(k)
        gz2(k)=gz2(k)*gz2(k+1)+uc(k)
 12   continue

      return
      end
c
c
c
      subroutine splint(xa,gx,gy,gz,gx2,gy2,gz2,n,x1,x,y,z)
      dimension xa(n),gx(n),gy(n),gz(n),gx2(n),gy2(n),gz2(n)

      klo=1
      khi=n
      
 1    if(khi-klo.gt.1)then
        k=(khi+klo)/2
        if(xa(k).gt.x1)then
          khi=k
        else
          klo=k
        end if
        goto 1
      end if

      h=xa(khi)-xa(klo)
      if(h.eq.0)then
        x=gx(x1)
        y=gy(x1)
        z=gz(x1)
      else
        a=(xa(khi)-x1)/h
        b=(x1-xa(klo))/h
        x=a*gx(klo)+b*gx(khi)+
     *    (a*(a*a-1.0)*gx2(klo)+b*(b*b-1.0)*gx2(khi))*h*h/6.
        y=a*gy(klo)+b*gy(khi)+
     *    (a*(a*a-1.0)*gy2(klo)+b*(b*b-1.0)*gy2(khi))*h*h/6.
        z=a*gz(klo)+b*gz(khi)+
     *    (a*(a*a-1.0)*gz2(klo)+b*(b*b-1.0)*gz2(khi))*h*h/6.
      end if

      return
      end



