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_int_mod
  use wei_cig_mod
  use wei_type_mod
  use wei_report_mod
  implicit none

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

  !!.. input history
  call wei_input_history()
  call wei_image_history()
  call wei_io_imag_init()
  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")

  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) cycle
     if(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, ogrid, pgrid)

     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_crec(wblock%indx(iwb), aws%n, ish)
        call wei_io_csou(wblock%indx(iwb), aws%n)

        call wei_wfld_init(aws%n)
 
        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)

           do izs=1, azs%n
              modl_slice => model(:,:,:,:,izs)
              vmod_slice => vmod(:,:,izs)
              vref_slice => vref(:,izs)
              stat = wei_vel_oper(vmod_slice, vref_slice)
              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_multi_ref_kernel(.true., +1, ws, vref_slice, vmod_slice, wfld_slice_d, FXop, WIop)
                 !stat = wei_multi_ref_kernel(.true., -1, ws, vref_slice, vmod_slice, wfld_slice_u, FXop, WIop)
                 stat = wei_ssf_oper(.true., +1, ws, vref_slice(1), vmod_slice, wfld_slice_d)
                 stat = wei_ssf_oper(.true., -1, ws, vref_slice(1), vmod_slice, wfld_slice_u)
                 stat = wei_cig_oper(wfld_slice_d, wfld_slice_u, modl_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

     end do

     stat = wei_vel_clos()
     stat = wei_ssf_clos()
     stat = wei_cig_clos()

     call wei_modl_dealloc()
     call wei_wfld_dealloc()

  end do


  do ish=1, nsh
     call wei_global_calc(ish)

     do iws=1, nws
        do izs=1, nzs
           stat = wei_ssf_oper(wfld_slice_d)
           call wei_wfld_io("rite")
        end do
     end do

     do iws=1, nws
        call wfld_copying()
        do izs=nzs, 1, -1
           stat = wei_ssf_oper(wfld_slice_w)
           call wei_wfld_io("reed")
           wfld_slice_u = wfld_slice_u + wfld_slice_d*wfld_slice_w
        end do
     end do 

     call wei_data_io()

  end do













end program
