c
c * 27-node brick (triquadratic interpolation)
c * http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/th04.htm
c
      program grintbx1
c
      dimension G(27),Fx(27),Fy(27),Fz(27)
      dimension wavex(257,28),wintx(257)
      dimension wavey(257,28),winty(257)
      dimension wavez(257,28),wintz(257)
c
      open(1,file='xyz_dat',status='unknown')
       read(1,*)nt,nob
       write(*,*)'nt,nob=',nt,nob
      close(1)
c
      x1=0.0
      x2=1200.0
      y1=0.0
      y2=1200.0
      z1=0.0
      z2=1200.0
      dx=x2-x1
      dy=y2-y1
      dz=z2-z1
      x0=(x1+x2)/2.0
      y0=(y1+y2)/2.0
      z0=(z1+z2)/2.0
c
      write(*,*)'Enter Location (x,y,z):'
      read(*,*)x,y,z
      if((x.lt.x1).and.(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).and.(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).and.(z.gt.z2))then
        write(*,*)'z is not within the domain !!'
        write(*,*)'z1,z2,z =',z1,z2,z
        stop
      end if
c
      r=2.0*(x-x0)/dx
      s=2.0*(y-y0)/dy
      t=2.0*(z-z0)/dz
      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
      open(1,file='x_dat',status='unknown')
      open(2,file='y_dat',status='unknown')
      open(3,file='z_dat',status='unknown')
       do j=1,nob
         do i=1,nt
           read(1,*)wavex(i,j)
           read(2,*)wavey(i,j)
           read(3,*)wavez(i,j)
         end do
       end do
      close(1)
      close(2)
      close(3)
c
      do 100 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
 100  continue
c
      open(1,file='wave1_dat',status='unknown')
       do i=1,nt
         write(1,*)wintx(i),winty(i),wintz(i)
       end do
      close(1)
c
      stop
      end

