program main
  use sep
  use grid_mod
  use wei_ssf_mod
  use image_gather_mod
  use filter_mod
  use myfft_mod
  use wei_global_mod
  use wei_param_mod

  implicit none

  real, allocatable, dimension(:,:,:)    :: data, wfld_d_data
  complex, allocatable, dimension(:,:,:) :: cdata

  real, allocatable, dimension(:)        :: sfun
  complex, allocatable, dimension(:)     :: csfun, ctrace

  complex, allocatable, dimension(:,:,:)   ,target :: wfld_d, wfld_u
  real   , allocatable, dimension(:,:,:)   ,target :: model
  complex, dimension(:,:,:)  , pointer :: pwfld_d, pwfld_u
  real   , dimension(:,:), pointer :: pmodel

  integer    :: nz, iz, stat, izs
  real       :: vel 
  type(axis) :: agx, agy, aw, asx, asy
  type(axis) :: az, ahx, ahy, akx, aky
  
  integer, parameter  :: opsign = 1
  real   , parameter  :: PI = 3.1415926
  logical, parameter  :: T = .true., F = .false.

  integer             :: iw, igx, igy, nt, it, nw1, nw2, nw3, nw4, nw
  real                :: oz, ow, ot, dz, dt, dw, f1, f2, f3, f4, w1, w2, w3, w4

  !!.. nbw, the number of w blocks
  !!.. nbz, the number of z blocks

  nwblock = aw%n / dwblock
  nzblock = az%n / dzblock

  do iwblock=1, nwblock
     allocate(cdata(agx%n,agy%n,dwblock))
     allocate(csfun(dwblock))
     write(0,*)"get grid"
     call get_computation_grid(asx, asy, agx, agy, az)
     allocate(wfld_u(grid%x%n,grid%y%n,dwblock))
     allocate(wfld_d(grid%x%n,grid%y%n,dwblock))
     allocate(model(grid%x%n,grid%y%n,grid%z%n)) 
     wfld_d(grid%iosx, grid%iosy, :) = csfun
     wfld_u(grid%iogx:grid%iogx+size(wfld_u, 1)-1, grid%iogy:grid%iogy+size(wfld_u, 2)-1, :) = cdata
     deallocate(cdata, csfun)
     pwfld_d => wfld_d
     pwfld_u => wfld_u
     !pmodel  => model

     akx%n=grid%x%n
     akx%d=2.0*PI/(akx%n*grid%x%d)
     akx%o= 0.

     aky%n=grid%y%n
     aky%d=2.0*PI/(aky%n*grid%y%d)
     aky%o= 0.

     awblock%o = aw%o + (iwblock-1) * dwblock
     awblock%d = aw%d
     awblock%n = dwblock
     write(0,*)"ssf init"
     call wei_ssf_init(akx, aky, az%d, awblock)

     do izblock=1, azblock%n  !!.. izblock should be in the outer loop
        do iz=1, azblock%d
           call sep_read(vmod, "vmod")
           call wei_get_vref_slice( F )
           stat = wei_ssf_down( F, pwfld_d )
           stat = wei_ssf_down( T, pwfld_u )

           if(iwblock.gt.1) call sep_read(model, "image")
           do iw=1, nw
              model(:,:) = model(:,:) + real(conjg(pwfld_d(:,:,iw))*pwfld_u(:,:,iw))
           end do
           call sep_write(model, "image")
        end do
     end do
  end do

end program









