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 (dgrid%hdr_rec) call wei_io_head_init("hrec")
  if (dgrid%hdr_sou) call wei_io_head_init("hsou")

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

  !!.. shot-profile migration
  do ish=1, dgrid%nsh
     if (dgrid%hdr_rec) call wei_io_head_oper("hrec")
     if (dgrid%hdr_sou) call wei_io_head_oper("hsou")
     call wei_counter_init()
     call wei_global_calc(ish)
     call wei_grids_report(ish, dgrid%nsh)
     if(igrid%ax%n>1.and.ogrid%ax%n == 0 .or. &
        igrid%ay%n>1.and.ogrid%ay%n == 0) cycle
     nws_bak = wblock%ns
     nzs_bak = zblock%ns
     call wei_wfld_alloc(wblock%ns)
     call wei_modl_alloc(zblock%ns)
     call wei_memory_report()
     stat = wei_vel_init(lloyd, nref)
     stat = wei_ssf_init(kx2, ky2, zblock%ds, maxang)
     stat = wei_cig_init(decon, eps, hessian_type, ogrid, pgrid)
     if (taper) then
        stat = wei_taper_2d_init(taper_type, pgrid%nxpad_beg, pgrid%nxpad_end, pgrid%nypad_beg, pgrid%nypad_end)
        write(0,*)"======== taper on ============"
     end if
     do iwb=1, wblock%nb
        aws%o = wblock%os(iwb); aws%d = wblock%ds; aws%n = wblock%blen(iwb)
        if (aws%n /= nws_bak) then
           call wei_wfld_dealloc()
           call wei_wfld_alloc(aws%n)
           nws_bak = aws%n
        end if
        call wei_io_csou(wblock%indx(iwb), aws%n)
        call wei_wfld_d_init(aws%n)
        call wei_wfsd_history(aws)
        if (taper) stat = wei_taper_2d_oper(wfld_d, taper_pct)
        !!.. downward propagate the source wavefiled for one frequency block
        !!.. all the way to each depth step and write them to the disk
        do izb=1, zblock%nb
           azs%o = zblock%os(izb); azs%d = zblock%ds; azs%n = zblock%blen(izb)
           if (azs%n /= nzs_bak) then
              call wei_modl_dealloc()
              call wei_modl_alloc(azs%n)
              nzs_bak = azs%n
           end if
           call wei_io_velo(zblock%indx(izb), azs%n)
           do izs=1, azs%n
              vmod_slice => vmod(:,:,izs)
              vref_slice => vref(:,izs)
              stat = wei_vel_oper(vmod_slice, vref_slice)
              !if (taper) stat = wei_taper_2d_oper(wfld_d, taper_pct)
              do iws=1, aws%n
                 ws = aws%o + (iws-1)*aws%d
                 wfld_slice_d => wfld_d(:,:,iws)
                 stat = wei_ssf_oper(T, +1, ws, vref_slice(1), vmod_slice, wfld_slice_d)
              end do
              indz_wfsd = zblock%indx(izb) + izs - 1
              call wei_io_wfsd("rite", wfld_d, aws%n, indz_wfsd)
           end do
        end do
        !!...........................................................................
        !!.. shot profile migration
        if (adj) then  
           call wei_io_crec("reed", wblock%indx(iwb), aws%n, ish)
           call wei_wfld_u_init(adj, F, aws%n)
           if (taper) stat = wei_taper_2d_oper(wfld_u, taper_pct)
           do izb=1, zblock%nb
              azs%o = zblock%os(izb); azs%d = zblock%ds; azs%n = zblock%blen(izb)
              if (azs%n /= nzs_bak) then
                 call wei_modl_dealloc()
                 call wei_modl_alloc(azs%n)
                 nzs_bak = azs%n
              end if
              call wei_io_velo(zblock%indx(izb), azs%n)
              call wei_io_imag("reed", zblock%indx(izb), azs%n)
              if (hessian_type > 0) call wei_io_illu("reed", zblock%indx(izb), azs%n)
              do izs=1, azs%n
                 zs = azs%o + (izs-1)*azs%d
                 model_slice => model(:,:,:,:,izs)
                 vmod_slice  => vmod(:,:,izs)
                 vref_slice  => vref(:,izs)
                 stat = wei_vel_oper(vmod_slice, vref_slice)
                 indz_wfsd = zblock%indx(izb) + izs - 1
                 call wei_io_wfsd("reed", wfld_d, aws%n, indz_wfsd)
                 !if (taper) stat = wei_taper_2d_oper(wfld_u, taper_pct)
                 do iws=1, aws%n
                    ws = aws%o + (iws-1)*aws%d
                    wfld_slice_d => wfld_d(:,:,iws)
                    wfld_slice_u => wfld_u(:,:,iws)
                    stat = wei_ssf_oper(adj, -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)
                    istep = istep + 1
                    if (istep == astep) then
                       percent = nint(real(astep)/real(nstep)*100.)
                       write(0,*)"Finished:",percent,"%"
                       astep = astep + rstep
                    end if
                 end do
              end do
              call wei_io_imag("rite", zblock%indx(izb), azs%n)
           end do
        !!.. shot profile demigration
        else 
           wfld_u = 0.
           crec   = 0.
           do izb=zblock%nb, 1, -1
              azs%o = zblock%os(izb); azs%d = zblock%ds; azs%n = zblock%blen(izb)
              if (azs%n /= nzs_bak) then
                 call wei_modl_dealloc()
                 call wei_modl_alloc(azs%n)
                 nzs_bak = azs%n
              end if
              call wei_io_velo(zblock%indx(izb), azs%n)
              call wei_io_imag("reed", zblock%indx(izb), azs%n)
              do izs=azs%n, 1, -1
                 zs = azs%o + (izs-1)*azs%d
                 vmod_slice => vmod(:,:,izs)
                 vref_slice => vref(:,izs)
                 model_slice => model(:,:,:,:,izs)
                 stat = wei_vel_oper(vmod_slice, vref_slice)
                 indz_wfsd = zblock%indx(izb) + izs - 1
                 call wei_io_wfsd("reed", wfld_d, aws%n, indz_wfsd)
                 do iws=1, aws%n
                    ws = aws%o + (iws-1)*aws%d
                    wfld_slice_d => wfld_d(:,:,iws)
                    wfld_slice_u => wfld_u(:,:,iws)
                    stat = wei_cig_oper(adj, T, ws, zs, wfld_slice_d, wfld_slice_u, model_slice)
                    stat = wei_ssf_oper(adj, +1, ws, vref_slice(1), vmod_slice, wfld_slice_u)
                    istep = istep + 1
                    if (istep == astep) then
                       percent = nint(real(astep)/real(nstep)*100.)
                       write(0,*)"Finished:",percent,"%"
                       astep = astep + rstep
                    end if
                 end do
                 !if (taper) stat = wei_taper_2d_oper(wfld_u, taper_pct)
              end do
           end do
           if (taper) stat = wei_taper_2d_oper(wfld_u, taper_pct)
           call wei_wfld_u_init(adj, F, aws%n)
           call wei_io_crec("rite", wblock%indx(iwb), aws%n, ish) 
        end if
     end do
     stat = wei_vel_clos()
     stat = wei_ssf_clos()
     stat = wei_cig_clos()
     stat = wei_taper_2d_clos()
     call wei_modl_dealloc()
     call wei_wfld_dealloc()
  end do

  call sep_end_prog()
  call exit()

end program








