program main
  use base_struct_mod
  use base_ioutil_mod
  use base_spline_mod
  use cgstep_mod
  use solver_smp_mod 
  implicit none
  type(axis)                      :: axd, ayd, azd
  type(axis)                      :: axk, ayk, azk 
  real, dimension(:), allocatable :: modl, data
  integer                         :: stat, niter
  logical                         :: adj, add
  character(len=256)              :: modl_, data_

  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()

  data_ = "data"
  modl_ = "modl"

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


  if (adj) then
     call hinaxis(data_, axd, 1)
     call hinaxis(data_, ayd, 2)
     call hinaxis(data_, azd, 3)

     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 )

  else
     call hinaxis(modl_, axk, 1) 
     call hinaxis(modl_, ayk, 2)
     call hinaxis(modl_, azk, 3)

     call from_param("oxd", axd%o, 0.)
     call from_param("dxd", axd%d, 1.)
     call from_param("nxd", axd%n, 1 )

     call from_param("oyd", ayd%o, 0.)
     call from_param("dyd", ayd%d, 1.)
     call from_param("nyd", ayd%n, 1 )

     call from_param("ozd", azd%o, 0.)
     call from_param("dzd", azd%d, 1.)
     call from_param("nzd", azd%n, 1 )

  end if

write(0,*)"axk = ", axk
write(0,*)"ayk = ", ayk
write(0,*)"azk = ", azk
write(0,*)"axd = ", axd
write(0,*)"ayd = ", ayd
write(0,*)"azd = ", azd 

  allocate(modl(axk%n*ayk%n*azk%n))
  allocate(data(axd%n*ayd%n*azd%n))

  if (adj) then
     call sep_read(data, data_)
  else
     call sep_read(modl, modl_)
  end if

  stat = base_spln3d_init(axd,ayd,azd,axk,ayk,azk)

  if (adj) then
     call solver_smp(m=modl, d=data, Fop=base_spln3d_oper, stepper=cgstep, &
                     niter=niter, verb=.true.)
  else
     stat = base_spln3d_oper(adj, add, modl, data)
  end if

  stat = base_spln3d_clos()


  if (adj) then
     call houaxis(modl_, axk, 1)
     call houaxis(modl_, ayk, 2)
     call houaxis(modl_, azk, 3)
     call sep_write(modl, modl_)
  else
     call houaxis(data_, axd, 1)
     call houaxis(data_, ayd, 2)
     call houaxis(data_, azd, 3)
     call sep_write(data, data_)
  end if


end program
