program test_inter3d

  use fvn_interpol

  implicit none

  integer(kind=ip_kind),parameter :: nx=21,ny=42,nz=18
  integer(kind=ip_kind) :: i,j,k
  real(kind=dp_kind) :: f,fdata(nx,ny,nz),dble,pi,q,sin,x,xdata(nx),y,ydata(ny),z,zdata(nz)
  real(kind=dp_kind) :: tv

  intrinsic dble,sin
  f(x,y,z)=sin(x+2.*y+3.*z)
  do i=1,nx
        xdata(i)=2.*(dble(i-1)/dble(nx-1))
  end do
  do i=1,ny
        ydata(i)=2.*(dble(i-1)/dble(ny-1))
  end do
  do i=1,nz
        zdata(i)=2.*(dble(i-1)/dble(nz-1))
  end do
  do i=1,nx
        do j=1,ny
              do k=1,nz
                    fdata(i,j,k)=f(xdata(i),ydata(j),zdata(k))
              end do
        end do
  end do
  call init_random_seed()
  call random_number(x)
  call random_number(y)
  call random_number(z)
  q=fvn_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,fdata)
  tv=f(x,y,z)
  write(*,'("x y z ",3(f8.5))') x,y,z
  write(*,'("Calculated (real) value :",f8.5)')  tv
  write(*,'("fvn interpolation : ",f8.5)') q
  write(*,'("Relative fvn error :",e12.5)') abs((q-tv)/tv)
  end program