program main
  use wei_global_defi_mod
  use wei_global_calc_mod
  use wei_io_mod
  use wei_vel_mod
  use wei_ssf_mod
  use wei_cig_mod
  use wei_type_mod
  use wei_multi_ref_kernel_mod
  use wei_source_boundary_mod
  use wei_report_mod
  implicit none

  integer    :: ish, iwb, izb, iws, izs, stat
  integer    :: nws_bak, nzs_bak, percent
  integer    :: indz_wfsd 
  real       :: ws, zs
  type(axis) :: aws
  type(axis) :: azs
  complex, dimension(:,:)    , pointer :: wfld_slice_d
  complex, dimension(:,:)    , pointer :: wfld_slice_u
  real   , dimension(:,:,:,:), pointer :: model_slice
  real   , dimension(:,:,:,:), pointer :: normd_slice
  real   , dimension(:,:)    , pointer :: vmod_slice
  real   , dimension(:)      , pointer :: vref_slice

  complex, dimension(:,:,:)  , pointer :: wfld_tmp

  !logical, parameter                   :: T=.true., F=.false.

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

  !!.. input history

  call wei_param_init()
  call wei_data_history(adj)
  call wei_imag_history(adj)
  call wei_io_file_init(adj)

  call wei_radians_calc()
  call wei_block_calc()

  if (adj) then
     write(0,*)"************************************"
     write(0,*)"       Shot profile migration       "
     write(0,*)"************************************"
  else
     write(0,*)"************************************"
     write(0,*)"     Shot profile demigration       "
     write(0,*)"************************************"
  end if


  do ish=1, dgrid%nsh
     call wei_counter_init()
     call wei_global_calc(ish)
     call wei_wfld_alloc(1)
     call wei_modl_alloc(1)
     stat = wei_vel_init(lloyd, nref)
     stat = wei_ssf_init(kx2, ky2, igrid%az%d, maxang)
     stat = wei_cig_init(decon, eps, ogrid, pgrid)
     aws%o = dgrid%aw%o
     aws%d = dgrid%aw%d
     aws%n = dgrid%aw%n 
     do iws=1, dgrid%aw%n
        ws = dgrid%aw%o + (iws-1)*dgrid%aw%d
        call wei_io_csou(iws, 1)
        call wei_wfld_d_init(1)
        call wei_wfsd_history(aws)
        wfld_slice_d => wfld_d(:,:,1)
        do izs=1, igrid%az%n
           zs = igrid%az%o + (izs-1)*igrid%az%d
           call wei_io_velo(izs, 1)
           vmod_slice => vmod(:,:,1)
           vref_slice => vref(:,1)
           stat = wei_vel_oper(vmod_slice, vref_slice)
           stat = wei_ssf_oper(T, +1, ws, vref_slice(1), vmod_slice, wfld_slice_d)
           call wei_io_wfsd("rite", wfld_d, 1, izs)
        end do 
        if (adj) then
           model = 0.; wfld_u=0; wfld_d=0; vmod=0; vref=0
           call wei_io_crec("reed", iws, 1, ish)
           call wei_wfld_u_init(adj, F, 1)
           if (associated(wfld_slice_d)) nullify(wfld_slice_d)
           if (associated(wfld_slice_u)) nullify(wfld_slice_u)
           wfld_slice_u => wfld_u(:,:,1)
           wfld_slice_d => wfld_d(:,:,1) 
           do izs=1, igrid%az%n
              zs = igrid%az%o + (izs-1)*igrid%az%d
              call wei_io_velo(izs, 1)
              call wei_io_imag("reed", izs, 1)
              model_slice => model(:,:,:,:,1)
              vmod_slice  => vmod(:,:,1)
              vref_slice  => vref(:,1)
              stat = wei_vel_oper(vmod_slice, vref_slice)
              call wei_io_wfsd("reed", wfld_d, 1, izs)
              stat = wei_ssf_oper(T, -1, ws, vref_slice(1), vmod_slice, wfld_slice_u)
              stat = wei_cig_oper(adj, T, ws, zs, wfld_slice_d, wfld_slice_u, model_slice) 

              !call wei_io_wfsd("reed", wfld_d, 1, izs)
              !call wei_io_imag("reed", izs, 1)
              !model(:,:,1,1,1) = model(:,:,1,1,1) + conjg(wfld_slice_d) * wfld_slice_u

              call wei_io_imag("rite", izs, 1)
           end do
        else
           wfld_u = 0. ; crec = 0.
           wfld_d = 0  ; model = 0.; vmod= 0.; vref= 0.
           if (associated(wfld_slice_d)) nullify(wfld_slice_d)
           if (associated(wfld_slice_u)) nullify(wfld_slice_u)
           wfld_slice_u => wfld_u(:,:,1)
           wfld_slice_d => wfld_d(:,:,1) 
           do izs=igrid%az%n, 1, -1
              zs = igrid%az%o + (izs-1)*igrid%az%d
              call wei_io_velo(izs, 1)
              call wei_io_imag("reed", izs, 1)
              model_slice => model(:,:,:,:,1)
              vmod_slice  => vmod(:,:,1)
              vref_slice  => vref(:,1)
              stat = wei_vel_oper(vmod_slice, vref_slice)
              call wei_io_wfsd("reed", wfld_d, 1, izs)
              stat = wei_cig_oper(adj, T, ws, zs, wfld_slice_d, wfld_slice_u, model_slice)
              stat = wei_ssf_oper(T, +1, ws, vref_slice(1), vmod_slice, wfld_slice_u)

              !call wei_io_wfsd("reed", wfld_d, 1, izs)
              !call wei_io_imag("reed", izs, 1)
              !wfld_slice_u = wfld_slice_u + wfld_slice_d * model(:,:,1,1,1)
           end do
           call wei_wfld_u_init(adj, T, 1)
           call wei_io_crec("rite", iws, 1, ish)
        end if
        if(mod(iws, dgrid%aw%n/5) == 0) then
              write(0,*)"Finished:", int(real(iws)/real(dgrid%aw%n)*100), "%"
        end if
  
     end do
     stat = wei_vel_clos()
     stat = wei_ssf_clos()
     stat = wei_cig_clos() 
  end do

  call sep_end_prog()
  call exit()

end program








