program main
  use sep
  use wemva_struct_mod
  use wemva_kerpix_incore_mod
  use cgstep_mod
  use solver_smp_mod
  use wemva_offdso_mod
  implicit none


  type(img_grid)                              :: igrid
  type(vel_grid)                              :: vgrid
  type(wemva_ftag)                            :: ftag_
  integer                                     :: iter, niter, stat
  real, dimension(:), allocatable             :: bslw, dslw, bimg, dimg
  character(len=256)                          :: temp_
 

  temp_ = "temp"

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

  stat = wemva_kerpix_init(vgrid, igrid, ftag_)
  stat = wemva_offdso_init(igrid%ax%n ,igrid%ay%n ,igrid%az%n ,igrid%ahx%n,igrid%ahy%n,igrid%ahz%n, &
                           igrid%ahx%o,igrid%ahy%o,igrid%ahz%o,igrid%ahx%d,igrid%ahy%d,igrid%ahz%d)

  allocate(bslw(vgrid%ax%n*vgrid%ay%n*vgrid%az%n))
  allocate(dslw(vgrid%ax%n*vgrid%ay%n*vgrid%az%n))
  allocate(bimg(igrid%ax%n*igrid%ay%n*igrid%az%n*igrid%ahx%n*igrid%ahy%n*igrid%ahz%n))
  allocate(dimg(igrid%ax%n*igrid%ay%n*igrid%az%n*igrid%ahx%n*igrid%ahy%n*igrid%ahz%n))

  call sseek(ftag_%bslw_, 0, 0)
  call sep_read(bslw, ftag_%bslw_)

  do iter=1, 1
     write(0,*)"========== iter ==========", iter
     stat = wemva_bimpix_oper(.false., ftag_%bslw_, ftag_%bimg_)
     stat = wemva_offdso_oper(.false., .false., ftag_%bimg_, ftag_%dimg_)
     call solver_smp(m=dslw, d=dimg, Fop=wemva_kerpix_oper, stepper=cgstep, niter=5, verb=.true.)
     bslw = bslw + dslw
     call sseek(ftag_%bslw_, 0, 0)
     call sep_write(bslw, ftag_%bslw_)

  end do
  


end program
