c
c * 27-node brick (triquadratic interpolation)
c * http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/th04.htm
c
      program grintbx5
      parameter(mt=257,mob=3400,mbox=290,mloc=20)
c
      dimension G(27),Fx(27),Fy(27),Fz(27),NODE(27)
      dimension wx(mt,mob,mbox),wy(mt,mob,mbox),wz(mt,mob,mbox)
      dimension wavex(mt,27),wavey(mt,27),wavez(mt,27)
      dimension wintx(mt),winty(mt),wintz(mt)
      dimension xo(mbox),yo(mbox),zo(mbox)
      dimension x1(mbox),y1(mbox),z1(mbox)
      dimension x2(mbox),y2(mbox),z2(mbox)
      dimension BXX(mbox),BXY(mbox),BXZ(mbox)
      dimension NBXX(mbox),NBXY(mbox),NBXZ(mbox)
      dimension nob(mbox),nt(mbox),dt(mbox)
      dimension xloc(mloc),yloc(mloc),zloc(mloc)
      character*20 fname(mbox)
c
c *** Input of Output Location (x,y,z) ***
c
      nloc=5
      xloc(1)=611.0
      yloc(1)=293.0
      zloc(1)=12.0
      xloc(2)=903.0
      yloc(2)=493.0
      zloc(2)=0.0
      xloc(3)=711.0
      yloc(3)=907.0
      zloc(3)=112.0
      xloc(4)=290.0
      yloc(4)=555.0
      zloc(4)=149.9
      xloc(5)=611.0
      yloc(5)=293.0
      zloc(5)=152.0
c      write(*,*)'i,x,y,z=',i,xloc(i),yloc(i),zloc(i)
c
c *** Input of Wave Data at 1st Box Nodes (from wave_d1) ***
c
      nbox=5
      fname(1)='wave_d1'
      fname(2)='wave_d2'
      fname(3)='wave_d3'
      fname(4)='wave_d4'
      fname(5)='wave_d5'
c
      do ibox=1,nbox
        open(1,file=fname(ibox),status='unknown')
         read(1,*)xo(ibox),yo(ibox),zo(ibox)
         read(1,*)BXX(ibox),BXY(ibox),BXZ(ibox)
         read(1,*)NBXX(ibox),NBXY(ibox),NBXZ(ibox)
         read(1,*)nob(ibox),nt(ibox),dt(ibox)
         write(*,*)'1st nob,nt,dt=',nob(ibox),nt(ibox),dt(ibox)
         if(nob(ibox).gt.mob)then
           write(*,*)'nob>mob !!',nob(ibox),mob
           stop
         end if
         if(nt(ibox).gt.mt)then
           write(*,*)'nt>mt !!',nt(ibox),mt
           stop
         end if
         do j=1,nob(ibox)
           do i=1,nt(ibox)
             read(1,*)wx(i,j,ibox),wy(i,j,ibox),wz(i,j,ibox)
           end do
         end do       
        close(1)
      end do
c
c *** Boundaries of Boxes ***
c
      do i=1,nbox
        x1(i)=xo(i)
        x2(i)=x1(i)+BXX(i)*NBXX(i)
        y1(i)=yo(i)
        y2(i)=y1(i)+BXY(i)*NBXY(i)
        z1(i)=zo(i)
        z2(i)=z1(i)+BXZ(i)*NBXZ(i)
c        write(*,*)'x1,x2',i,x1(i),x2(i)
c        write(*,*)'y1,y2',i,y1(i),y2(i)
c        write(*,*)'z1,z2',i,z1(i),z2(i)
      end do
c
c *** Big-Loop of Observation Points and Opne Output File ***
c
      open(1,file='wave2_dat',status='unknown')
      write(1,2000)nloc,nt(1),dt(1)
c
      do 100 iloc=1,nloc
        x=xloc(iloc)
        y=yloc(iloc)
        z=zloc(iloc)
c        write(*,*)
        write(*,*)'iloc,x,y,z',iloc,x,y,z
c
c *** Finding the target Box Num. ***
c
        do i=1,nbox
          if((x.ge.x1(i)).and.(x.le.x2(i)))then
            if((y.ge.y1(i)).and.(y.le.y2(i)))then
              if((z.ge.z1(i)).and.(z.le.z2(i)))then
                iiBOX=i
c                write(*,*)'iiBOX',iiBOX
                goto 10
              end if
            end if
          end if
        end do
c
        write(*,*)'(x,y,z) is outside of the Boxes !!'
        write(*,*)'iloc,x,y,z',iloc,x,y,z
        write(*,*)'x1,x2',x1(i),x2(i)
        write(*,*)'y1,y2',y1(i),y2(i)
        write(*,*)'z1,z2',z1(i),z2(i)
        stop
c
 10     NBOXX=NBXX(iiBOX)
        NBOXY=NBXY(iiBOX)
        NBOXZ=NBXZ(iiBOX)
        BOXX=BXX(iiBOX)
        BOXY=BXY(iiBOX)
        BOXZ=BXZ(iiBOX)
        xoo=xo(iiBOX)
        yoo=yo(iiBOX)
        zoo=zo(iiBOX)
        nobs=nob(iiBOX)
        ntime=nt(iiBOX)
        dtime=dt(iiBOX)
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
c        write(*,*)'iBOXx,y,z',iBOXx,iBOXy,iBOXz
        NBOX1=(iBOXz-1)*NBOXX*NBOXY+(iBOXy-1)*NBOXX+iBOXx
c        write(*,*)'NBOX(Box Num)=',NBOX1
c
c *** Finding the 27 Node numbers of the target Box ***
c
        i=NBOX1/(NBOXX*NBOXY)
        j=mod(NBOX1,NBOXX*NBOXY)
c        write(*,*)'i,j',i,j
        ii=i+1
        jj=j      
        if(j.eq.0)then
          ii=i
          jj=NBOXX*NBOXY
        end if
c        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
c        write(*,*)'kk(xBOX),ll(+BOX)',iBOX,kk,ll
        jNODE=NODEx*2*kk+(ll*2-1)
c        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
c        write(*,*)'NODE 1, 2, 3',NODE(1),NODE(2),NODE(3)
c        write(*,*)'NODE 4, 5, 6',NODE(4),NODE(5),NODE(6)
c        write(*,*)'NODE 7, 8, 9',NODE(7),NODE(8),NODE(9)
c        write(*,*)'NODE10,11,12',NODE(10),NODE(11),NODE(12)
c        write(*,*)'NODE13,14,15',NODE(13),NODE(14),NODE(15)
c        write(*,*)'NODE16,17,18',NODE(16),NODE(17),NODE(18)
c        write(*,*)'NODE19,20,21',NODE(19),NODE(20),NODE(21)
c        write(*,*)'NODE22,23,24',NODE(22),NODE(23),NODE(24)
c        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 200 j=1,nobs
           do inode=1,27
             if(j.eq.NODE(inode))then
                 do i=1,ntime
                   wavex(i,inode)=wx(i,j,iiBOX)
                   wavey(i,inode)=wy(i,j,iiBOX)
                   wavez(i,inode)=wz(i,j,iiBOX)
                 end do
             end if
           end do         
 200    continue
c
c *** Calcuration of Wave using triquadratic interpolation and output ***
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
c        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 300 i=1,ntime
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
 300    continue
c
         do i=1,ntime
           write(1,2010)wintx(i),winty(i),wintz(i)
         end do
 2000   format(2i12,f12.5)
 2010   format(3e12.4)
c
 100  continue
      close(1)
c
      stop
      end

