program main
  use wei_global_defi_mod
  use wei_global_calc_mod
  use wei_io_mod
  use wei_shotmig_kernel_mod
  use wei_multi_ref_kernel_mod
  use wei_vel_mod
  use wei_ssf_mod
  use wei_int_mod
  use wei_cig_mod
  use wei_type_mod
  use wei_report_mod
  use sep
  implicit none

  integer    :: ish, iwb, izb, stat
  integer    :: nws_bak, nzs_bak
  integer    :: indx_dept, indx_freq
  type(axis) :: aws
  type(axis) :: azs


  real   , dimension(:)      , pointer :: vref_slice
  real   , dimension(:,:)    , pointer :: vmod_slice
  real   , dimension(:,:,:,:), pointer :: model_slice
  complex, dimension(:,:)    , pointer :: wfld_slice_d
  complex, dimension(:,:)    , pointer :: wfld_slice_u
  complex, dimension(:,:,:)  , pointer :: pwfld_d, pwfld_u
  real :: rvref, ws
  integer :: izs, iws
  !!.. nbw, the total number of w blocks
  !!.. nbz, the total number of z blocks

  !!.. input history
  call wei_input_history()
  call wei_image_history()
  call wei_io_imag_init()
  if (dgrid%hdr_rec) call wei_io_head_init("hrec")
  if (dgrid%hdr_sou) call wei_io_head_init("hsou")


  call to_history("n1", dgrid%agx%n, "test")
  call to_history("n2", dgrid%agy%n, "test")
  call to_history("n3", dgrid%aw%n, "test")

  write(0,*)"ntr", dgrid%ntr

  do ish=1, dgrid%nsh

     write(0,*)"ish=",ish, dgrid%nsh

     if (dgrid%hdr_rec) call wei_io_head_lop("hrec")
     if (dgrid%hdr_sou) call wei_io_head_lop("hsou")

     call wei_global_calc(ish)

     call wei_grids_report()

     if(igrid%ax%n>1.and.ogrid%ax%n == 0) cycle
     if(igrid%ay%n>1.and.ogrid%ay%n == 0) cycle 

     nws_bak = wblock%ns
     nzs_bak = zblock%ns

     call wei_wfld_alloc(1)
     call wei_modl_alloc(1)

     call wei_memory_report()

     stat = wei_vel_init(lloyd, nref)
     stat = wei_ssf_init(kx2, ky2, zblock%ds)
     stat = wei_cig_init(ogrid, pgrid)

     !call sep_read(crec, "rec", esize=8)
     !call sep_read(csou, "sou", esize=8)
     !call wei_wfld_init(nws_bak)

     !call wei_io_velo(1, nzs_bak)
     
     do iws=1, wblock%ns
        write(0,*)"iws=", iws
        ws = dgrid%aw%o + (iws-1)*dgrid%aw%d
        call sep_read(crec, "rec", esize=8)
        call sep_read(csou, "sou", esize=8)
        call wei_wfld_init(1)
        wfld_slice_d => wfld_d(:,:,1)
        wfld_slice_u => wfld_u(:,:,1)
        do izs=1, zblock%ns
           call wei_io_velo(izs, 1)
           call wei_io_imag("reed", izs, 1)
           vmod_slice => vmod(:,:,1)
           model_slice => model(:,:,:,:,1)
           rvref = sum(vmod_slice)/size(vmod_slice)
           stat = wei_ssf_lop(.true., +1, ws, rvref, vmod_slice, wfld_slice_d)
           stat = wei_ssf_lop(.true., -1, ws, rvref, vmod_slice, wfld_slice_u)
           stat = wei_cig_lop(wfld_slice_d, wfld_slice_u, model_slice)
           call wei_io_imag("rite", izs, 1)
        end do
     end do

     !call wei_io_imag("rite", 1, nzs_bak)    

     call wei_modl_dealloc()
     call wei_wfld_dealloc()

  end do

end program
