c
c * 27-node brick (triquadratic interpolation)
c * http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/th04.htm
c
      program grintbx1
      parameter(mt=257,mob=250)
c
      dimension G(27),Fx(27),Fy(27),Fz(27),NODE(27)
      dimension wx(mt,mob),wy(mt,mob),wz(mt,mob)
      dimension wavex(mt,27),wavey(mt,27),wavez(mt,27)
      dimension wintx(mt),winty(mt),wintz(mt)
c
c *** Wave Input Data at Box Nodes ***
c
      open(1,file='wave_d',status='unknown')
       read(1,*)xoo,yoo,zoo
       read(1,*)BOXX,BOXY,BOXZ
       read(1,*)NBOXX,NBOXY,NBOXZ
       read(1,*)nob,nt,dt
       write(*,*)'nob,nt,dt=',nob,nt,dt
       if(nob.gt.mob)then
         write(*,*)'nob>mob !!',nob,mob
         stop
       end if
       if(nt.gt.mt)then
         write(*,*)'nt>mt !!',nt,mt
         stop
       end if
       do j=1,nob
         do i=1,nt
           read(1,*)wx(i,j),wy(i,j),wz(i,j)
         end do
       end do       
      close(1)
c
      x1=xoo
      x2=x1+BOXX*NBOXX
      y1=yoo
      y2=y2+BOXY*NBOXY
      z1=zoo
      z2=z1+BOXZ*NBOXZ
c
      write(*,*)'x1,x2',x1,x2
      write(*,*)'y1,y2',y1,y2
      write(*,*)'z1,z2',z1,z2
c
c *** Input of Output Location (x,y,z), and find the target Box Num. ***
c
      write(*,*)'Enter Location (x,y,z):'
      read(*,*)x,y,z
      if((x.lt.x1).or.(x.gt.x2))then
        write(*,*)'x is not within the domain !!'
        write(*,*)'x1,x2,x =',x1,x2,x
        stop
      end if
      if((y.lt.y1).or.(y.gt.y2))then
        write(*,*)'y is not within the domain !!'
        write(*,*)'y1,y2,y =',y1,y2,y
        stop
      end if
      if((z.lt.z1).or.(z.gt.z2))then
        write(*,*)'z is not within the domain !!'
        write(*,*)'z1,z2,z =',z1,z2,z
        stop
      end if
c
      do i=1,NBOXX
        xx1=xoo+BOXX*(i-1)
        xx2=xoo+BOXX*i
        if((x.ge.xx1).and.(x.le.xx2))iBOXx=i
      end do
      do i=1,NBOXY
        yy1=yoo+BOXY*(i-1)
        yy2=yoo+BOXY*i
        if((y.ge.yy1).and.(y.le.yy2))iBOXy=i
      end do
      do i=1,NBOXZ
        zz1=zoo+BOXZ*(i-1)
        zz2=zoo+BOXZ*i
        if((z.ge.zz1).and.(z.le.zz2))iBOXz=i
      end do
      write(*,*)'iBOXx,y,z',iBOXx,iBOXy,iBOXz
      NBOX=(iBOXz-1)*NBOXX*NBOXY+(iBOXy-1)*NBOXX+iBOXx
      write(*,*)'NBOX(Box Num)=',NBOX
c
c *** Finding the 27 Node numbers of the target Box ***
c
      i=NBOX/(NBOXX*NBOXY)
      j=mod(NBOX,NBOXX*NBOXY)
      write(*,*)'i,j',i,j
      ii=i+1
      jj=j      
      if(j.eq.0)then
        ii=i
        jj=NBOXX*NBOXY
      end if
      write(*,*)'ii(row),jj(box)',ii,jj
c
      NODEx=NBOXX*2+1
      NODEy=NBOXY*2+1
      NODExy=NODEx*NODEy
      iNODE=(ii-1)*NODExy*2
      k=jj/NBOXX
      l=MOD(jj,NBOXX)
      kk=k
      ll=l
      if(l.eq.0)then
        kk=k-1
        ll=NBOXX
      end if
      write(*,*)'kk(xBOX),ll(+BOX)',iBOX,kk,ll
      jNODE=NODEx*2*kk+(ll*2-1)
      write(*,*)'iNODE(row),jNODE(+NODE)',iNODE,jNODE
c
      NODE(1)=iNODE+jNODE
      NODE(2)=NODE(1)+1
      NODE(3)=NODE(2)+1
      NODE(4)=NODE(1)+NODEx
      NODE(5)=NODE(2)+NODEx
      NODE(6)=NODE(3)+NODEx
      NODE(7)=NODE(4)+NODEx
      NODE(8)=NODE(5)+NODEx
      NODE(9)=NODE(6)+NODEx
      NODE(10)=NODE(1)+NODExy
      NODE(11)=NODE(2)+NODExy
      NODE(12)=NODE(3)+NODExy
      NODE(13)=NODE(4)+NODExy
      NODE(14)=NODE(5)+NODExy
      NODE(15)=NODE(6)+NODExy
      NODE(16)=NODE(7)+NODExy
      NODE(17)=NODE(8)+NODExy
      NODE(18)=NODE(9)+NODExy
      NODE(19)=NODE(10)+NODExy
      NODE(20)=NODE(11)+NODExy
      NODE(21)=NODE(12)+NODExy
      NODE(22)=NODE(13)+NODExy
      NODE(23)=NODE(14)+NODExy
      NODE(24)=NODE(15)+NODExy
      NODE(25)=NODE(16)+NODExy
      NODE(26)=NODE(17)+NODExy
      NODE(27)=NODE(18)+NODExy
      write(*,*)'NODE 1, 2, 3',NODE(1),NODE(2),NODE(3)
      write(*,*)'NODE 4, 5, 6',NODE(4),NODE(5),NODE(6)
      write(*,*)'NODE 7, 8, 9',NODE(7),NODE(8),NODE(9)
      write(*,*)'NODE10,11,12',NODE(10),NODE(11),NODE(12)
      write(*,*)'NODE13,14,15',NODE(13),NODE(14),NODE(15)
      write(*,*)'NODE16,17,18',NODE(16),NODE(17),NODE(18)
      write(*,*)'NODE19,20,21',NODE(19),NODE(20),NODE(21)
      write(*,*)'NODE22,23,24',NODE(22),NODE(23),NODE(24)
      write(*,*)'NODE25,26,27',NODE(25),NODE(26),NODE(27)
c
c *** Extracting Seismic Waves at the 27 Nodes of the target Box ***
c
       do 100 j=1,nob
         do inode=1,27
           if(j.eq.NODE(inode))then
             do i=1,nt
               wavex(i,inode)=wx(i,j)
               wavey(i,inode)=wy(i,j)
               wavez(i,inode)=wz(i,j)
             end do
           end if
         end do       
 100   continue
c
c *** Calcuration of Wave using triquadratic interpolation ***
c
      x0=xoo+BOXX*(ll-1)+BOXX/2.0
      y0=yoo+BOXY*kk+BOXY/2.0
      z0=zoo+BOXZ*(ii-1)+BOXZ/2.0
      r=2.0*(x-x0)/BOXX
      s=2.0*(y-y0)/BOXY
      t=2.0*(z-z0)/BOXZ
      write(*,*)'r,s,t =',r,s,t
c
      G(1) =-r*s*t*(1.0-r)*(1.0-s)*(1.0-t)/8.0
      G(2) = s*t*(1.0-r*r)*(1.0-s)*(1.0-t)/4.0
      G(3) = r*s*t*(1.0+r)*(1.0-s)*(1.0-t)/8.0
      G(4) = r*t*(1.0-r)*(1.0-s*s)*(1.0-t)/4.0
      G(5) =-t*(1.0-r*r)*(1.0-s*s)*(1.0-t)/2.0
      G(6) =-r*t*(1.0+r)*(1.0-s*s)*(1.0-t)/4.0
      G(7) = r*s*t*(1.0-r)*(1.0+s)*(1.0-t)/8.0
      G(8) =-s*t*(1.0-r*r)*(1.0+s)*(1.0-t)/4.0
      G(9) =-r*s*t*(1.0+r)*(1.0+s)*(1.0-t)/8.0
      G(10)= r*s*(1.0-r)*(1.0-s)*(1.0-t*t)/4.0
      G(11)=-s*(1.0-r*r)*(1.0-s)*(1.0-t*t)/2.0
      G(12)=-r*s*(1.0+r)*(1.0-s)*(1.0-t*t)/4.0
      G(13)=-r*(1.0-r)*(1.0-s*s)*(1.0-t*t)/2.0
      G(14)= (1.0-r*r)*(1.0-s*s)*(1.0-t*t)
      G(15)= r*(1.0+r)*(1.0-s*s)*(1.0-t*t)/2.0
      G(16)=-r*s*(1.0-r)*(1.0+s)*(1.0-t*t)/4.0
      G(17)= s*(1.0-r*r)*(1.0+s)*(1.0-t*t)/2.0
      G(18)= r*s*(1.0+r)*(1.0+s)*(1.0-t*t)/4.0
      G(19)= r*s*t*(1.0-r)*(1.0-s)*(1.0+t)/8.0
      G(20)=-s*t*(1.0-r*r)*(1.0-s)*(1.0+t)/4.0
      G(21)=-r*s*t*(1.0+r)*(1.0-s)*(1.0+t)/8.0
      G(22)=-r*t*(1.0-r)*(1.0-s*s)*(1.0+t)/4.0
      G(23)= t*(1.0-r*r)*(1.0-s*s)*(1.0+t)/2.0
      G(24)= r*t*(1.0+r)*(1.0-s*s)*(1.0+t)/4.0
      G(25)=-r*s*t*(1.0-r)*(1.0+s)*(1.0+t)/8.0
      G(26)= s*t*(1.0-r*r)*(1.0+s)*(1.0+t)/4.0
      G(27)= r*s*t*(1.0+r)*(1.0+s)*(1.0+t)/8.0
c
      do 200 i=1,nt
c
        do j=1,27
          Fx(j)=wavex(i,j)
          Fy(j)=wavey(i,j)
          Fz(j)=wavez(i,j)
        end do
c
        Fx1=0.0
        Fy1=0.0
        Fz1=0.0
        do j=1,27
          Fx1=Fx1+G(j)*Fx(j)
          Fy1=Fy1+G(j)*Fy(j)
          Fz1=Fz1+G(j)*Fz(j)
        end do
        wintx(i)=Fx1
        winty(i)=Fy1
        wintz(i)=Fz1
c
 200  continue
c
      open(1,file='wave1_dat',status='unknown')
       write(1,2000)nob,nt,dt
       do i=1,nt
         write(1,2010)wintx(i),winty(i),wintz(i)
       end do
      close(1)
 2000 format(2i12,f12.5)
 2010 format(3e12.4)
c
      stop
      end

