program main
  use wemva_kerspl_mod
  use base_spline_mod

  type(vel_grid)                      :: vgrid
  type(img_grid)                      :: igrid
  type(wemva_ftag)                    :: ftag_
  type(axis)                          :: axk, ayk, azk

  logical                             :: adj, add
  integer                             :: stat

  real, dimension(:,:,:), allocatable :: bfnc, dslw
  character(len=256)                  :: bfnc_

  bfnc_ = "bfnc" 

  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()

  call from_param("adj", adj, .false.)
  call from_param("add", add, .false.)

  stat = wemva_prestk_init(vgrid, igrid, ftag_)

  call auxinout(bfnc_)

  if (adj) then
     call from_param("oxk", axk%o, 0.)
     call from_param("dxk", axk%d, 1.)
     call from_param("nxk", axk%n, 1 )

     call from_param("oyk", ayk%o, 0.)
     call from_param("dyk", ayk%d, 1.)
     call from_param("nyk", ayk%n, 1 )

     call from_param("ozk", azk%o, 0.)
     call from_param("dzk", azk%d, 1.)
     call from_param("nzk", azk%n, 1 )

     call houaxis(bfnc_, axk, 1)
     call houaxis(bfnc_, ayk, 2)
     call houaxis(bfnc_, azk, 3)

  else

     call hinaxis(bfnc_, axk, 1)
     call hinaxis(bfnc_, ayk, 2)
     call hinaxis(bfnc_, azk, 3)
     call houaxis(ftag_%dslw_, vgrid%ax, 1)
     call houaxis(ftag_%dslw_, vgrid%ay, 2)
     call houaxis(ftag_%dslw_, vgrid%az, 3)
  end if

  allocate(bfnc(axk%n,ayk%n,azk%n))
  allocate(dslw(vgrid%ax%n,vgrid%ay%n,vgrid%az%n))

  stat = base_spln2d_init(vgrid%ax,vgrid%az,axk,azk)

  if (adj) then
     write(0,*)"======== Adjoint ================"
     stat = wemva_prestk_oper(adj, add, ftag_%dslw_, ftag_%dimg_)
     call sseek(ftag_%dslw_, 0, 0)
     call sep_read(dslw, ftag_%dslw_)
     stat = base_spln2d_oper(adj, add, bfnc, dslw)
     call sep_write(bfnc, bfnc_)
  else
     write(0,*)"======== Forward ================"
     call sep_read (bfnc, bfnc_) 
     stat = base_spln2d_oper(.false., add, bfnc, dslw)
     call sseek(ftag_%dslw_, 0, 0)
     call sep_write(dslw, ftag_%dslw_)
     stat = wemva_prestk_oper(adj, add, ftag_%dslw_, ftag_%dimg_) 
  end if     
  stat = base_spln2d_clos()
  stat = wemva_prestk_clos()

end program







