program main
  use sep
  use base_struct_mod
  use base_blkcal_mod
  use base_ioutil_mod
  use wem_velref_mod
  use wem_wfdexp_mod
  use wem_tapers_mod
  use wem_fftwv3_mod
  implicit none

  type(axis)                                 :: ax, ay, az, aw, aws
  type(axis)                                 :: arx, ary, avx, avy
  type(block)                                :: wblock
  complex  , dimension(:,:,:)  , allocatable :: wfdr_slice
  complex  , dimension(:,:,:)  , allocatable :: crec
  real     , dimension(:,:,:)  , allocatable :: modl
  real     , dimension(:,:,:)  , allocatable :: vmod
  real     , dimension(:,:)    , allocatable :: vref

  real     , dimension(:)      , allocatable :: kx2, ky2
  real     , dimension(:)      , allocatable :: tapx, tapy
  real     , dimension(:)      , allocatable :: hamw

  integer  , dimension(:,:,:)  , allocatable :: iwrk2d
  real     , dimension(:,:)    , allocatable :: rwrk1d
  complex  , dimension(:,:)    , allocatable :: cwrk1d
  complex  , dimension(:,:,:,:), allocatable :: cwrk3d

  integer*8, dimension(:)      , allocatable :: plan_2drF, plan_2drI

  character(len=5)                           :: oper
  integer                                    :: nxpad_beg, nxpad_end
  integer                                    :: nypad_beg, nypad_end
  integer                                    :: pdip, nvrf, ntpx, ntpy
  real                                       :: trick, ws, PI
  real                                       :: rx_min, rx_max
  real                                       :: ry_min, ry_max
  real                                       :: px_min, px_max
  real                                       :: py_min, py_max
  real                                       :: f1, f2, f3, f4
  real                                       :: w1, w2, w3, w4

  character(len=256)                         :: vmod_, crec_, imag_

  logical                                    :: adj, radian
  integer                                    :: ix, iy, iz, iwb, iws, indw
  integer                                    :: nw1, nw2, nw3, nw4, iw
  integer                                    :: ith, nth, ompchunk, stat 
  integer                         , external :: omp_get_num_procs
  integer                         , external :: omp_get_num_threads
  integer                         , external :: omp_set_num_threads
  integer                         , external :: omp_get_thread_num

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

  PI = 3.14159265359

  vmod_ = "vmod"
  crec_ = "crec"
  imag_ = "imag"

  ompchunk = 1
  nth      = 1
  ith      = 1

!$OMP PARALLEL
    if (omp_get_thread_num() == 0) then
       nth = omp_get_num_threads()
       write(0,*)"nthread =", nth, omp_get_num_procs()
    end if
!$OMP END PARALLEL


  call from_param("pdip"     , pdip     ,      45)
  call from_param("nvrf"     , nvrf     ,       1)
  call from_param("trick"    , trick    ,    0.14)
  call from_param("ntpx"     , ntpx     ,       1)
  call from_param("ntpy"     , ntpy     ,       1)
  call from_param("oper"     , oper     , "ssfpi")
  call from_param("nws"      , wblock%ns,       1)
  call from_param("adj"      , adj      ,  .true.)
  call from_param("radian"   , radian   ,  .true.)

  call from_param("nxpad_beg", nxpad_beg,       0)
  call from_param("nxpad_end", nxpad_end,       0)
  call from_param("nypad_beg", nypad_beg,       0)
  call from_param("nypad_end", nypad_end,       0)

  call from_param("f1"       , f1       ,      5.)
  call from_param("f2"       , f2       ,     10.)
  call from_param("f3"       , f3       ,     35.)
  call from_param("f4"       , f4       ,     40.)


  call auxin(vmod_)
  call hinaxis(vmod_, avx, 1); call hinaxis(vmod_, avy, 2); call hinaxis(vmod_, az , 3)


  if (adj) then
     write(0,*)"postack migration"
     call auxin (crec_); call auxout(imag_)
     call hinaxis(crec_, arx, 1); call hinaxis(crec_, ary, 2); call hinaxis(crec_, aw, 3)
  else
     write(0,*)"postack exploding reflector modeling"
     call auxin (imag_); call auxout(crec_)
     call hinaxis(imag_, arx, 1); call hinaxis(imag_, ary, 2)
     call from_param("nw", aw%n, 1); call from_param("dw", aw%d, 1.); call from_param("ow", aw%o, 0.)
  end if 

  if (.not.radian) then
     aw%d = aw%d*2*PI; aw%o = aw%o*2*PI 
  end if

  if (arx%n/=avx%n .or. ary%n/=avy%n) call seperr("ERROR")

  rx_min = arx%o; rx_max = arx%o + (arx%n-1)*arx%d
  ry_min = ary%o; ry_max = ary%o + (ary%n-1)*ary%d

  px_min = rx_min - nxpad_beg*arx%d; px_max = rx_max + nxpad_end*arx%d
  py_min = ry_min - nypad_beg*ary%d; py_max = ry_max + nypad_end*ary%d

  ax%o = px_min; ax%d = arx%d; ax%n = 1.5 + (px_max-px_min)/ax%d
  ay%o = py_min; ay%d = ary%d; ay%n = 1.5 + (py_max-py_min)/ay%d

  if (adj) then
     call houaxis(imag_, arx, 1); call houaxis(imag_, ary, 2); call houaxis(imag_, az, 3)
  else
     call houaxis(crec_, arx, 1); call houaxis(crec_, ary, 2); call houaxis(crec_, aw, 3)
     call auxputch("esize", "i", 8, crec_)
  end if

  call base_blkcal(wblock, aw)

  allocate(modl(ax%n,ay%n,az%n))
  allocate(vmod(ax%n,ay%n,az%n))
  allocate(vref(nvrf,az%n)) 
  allocate(crec(ax%n,ay%n,wblock%ns))
  allocate(wfdr_slice(ax%n,ay%n,nth))
  allocate(kx2(ax%n),ky2(ay%n))
  allocate(tapx(ntpx),tapy(ntpy))
  allocate(hamw(aw%n))

  allocate(plan_2drF(nth),plan_2drI(nth))

  allocate(cwrk1d(ax%n*5+ay%n*5,nth)) !!.. ifd or ffd
  allocate(cwrk3d(ax%n,ay%n,2,nth))   !!.. ssf-epi
  allocate(iwrk2d(ax%n,ay%n,nth))     !!.. ssf-epi
  allocate(rwrk1d(ax%n+ay%n,nth))     !!.. ifd or ffd

  call sep_read(vmod(nxpad_beg+1:nxpad_beg+arx%n,nypad_beg+1:nypad_beg+ary%n,:), vmod_)

  !! expand domain in x direction
  do ix=1, nxpad_beg
     vmod(       ix,:,:) = vmod(   1+nxpad_beg,:,:)
  end do
  do ix=1, nxpad_end
     vmod(ax%n-ix+1,:,:) = vmod(ax%n-nxpad_end,:,:)
  end do

  !! expand domain in y direction
  do iy=1, nypad_beg
     vmod(:,       iy,:) = vmod(:,   1+nypad_beg,:)
  end do
  do iy=1, nypad_end
     vmod(:,ay%n-iy+1,:) = vmod(:,ay%n-nypad_end,:)
  end do

  vmod = vmod * 0.5

  stat = wem_velref_init(nvrf)
  do iz=1, az%n
     stat = wem_velref_oper(vmod(:,:,iz), vref(:,iz))
  end do

  stat = wem_wfdexp_init(kx2, ky2, ax%n, ay%n, ax%d, ay%d, az%d, pdip, nvrf, oper, trick)
  stat = wem_tapers_init(tapx, tapy, ntpx, ntpy)

  w1 = 2*PI*f1; nw1 = 1.5+(w1-aw%o)/aw%d
  w2 = 2*PI*f2; nw2 = 1.5+(w2-aw%o)/aw%d 
  w3 = 2*PI*f3; nw3 = 1.5+(w3-aw%o)/aw%d
  w4 = 2*PI*f4; nw4 = 1.5+(w4-aw%o)/aw%d

  do iw=1, aw%n
     if      (iw>=nw1 .and. iw<=nw2) then
        hamw(iw) = 0.54+0.46*cos(PI*(iw-nw1)/(nw2-nw1)-PI)
     else if (iw>=nw3 .and. iw<=nw4) then
        hamw(iw) = 0.54+0.46*cos(PI*(nw3-iw)/(nw4-nw3)   )
     else if (iw <nw1 .or.  iw> nw4) then
        hamw(iw) = 0.
     else 
        hamw(iw) = 1.
     end if
  end do

  write(0,*)"hamw=", hamw

  call k2_map(kx2, ax%n, ax%d)
  call k2_map(ky2, ay%n, ay%d)
  call wem_fftw_init(plan_2drF, plan_2drI, wfdr_slice, ax%n, ay%n, nth)


  if (adj) then  !!.. postack migration
     modl = 0.
     do iwb=1, wblock%nb
        aws%o = wblock%os(iwb); aws%d = wblock%ds; aws%n = wblock%blen(iwb)
        
        call sep_read(crec(nxpad_beg+1:nxpad_beg+arx%n,nypad_beg+1:nypad_beg+ary%n,1:aws%n), crec_)

        wfdr_slice = 0.; iwrk2d = 0; rwrk1d = 0.; cwrk1d = 0.; cwrk3d = 0.

!$OMP PARALLEL DO SCHEDULE(DYNAMIC,OMPCHUNK) DEFAULT(SHARED) PRIVATE(iws,iz,ith,ws)
        do iws=1, aws%n
           ith  = omp_get_thread_num() + 1
           ws   = aws%o + (iws-1)*aws%d
           indw = 1.5 + (ws-aw%o)/aw%d

           if (iws/10*10==iws) write(0,*)"iws=", iws, ws, ith

           wfdr_slice(:,:,ith) = crec(:,:,iws)
           do iz=1, az%n
              stat = wem_wfdexp_oper(wfdr_slice(:,:,ith), vmod(:,:,iz), vref(:,iz), ws, -1,                          &
                     iwrk2d(:,:,ith), rwrk1d(:,ith), cwrk1d(:,ith), cwrk3d(:,:,:,ith), plan_2drF(ith), plan_2drI(ith))
              stat = wem_tapers_oper(wfdr_slice(:,:,ith), ax%n, ay%n)
              modl(:,:,iz) = modl(:,:,iz) + wfdr_slice(:,:,ith)*hamw(indw)
           end do
        end do
!$OMP END PARALLEL DO 
     end do
     call sep_write(modl(nxpad_beg+1:nxpad_beg+arx%n,nypad_beg+1:nypad_beg+ary%n,:), imag_)

  else           !!.. exploding reflector modeing

     call sep_read(modl(nxpad_beg+1:nxpad_beg+arx%n,nypad_beg+1:nypad_beg+ary%n,:), imag_)
     do iwb=1, wblock%nb
        aws%o = wblock%os(iwb); aws%d = wblock%ds; aws%n = wblock%blen(iwb)

        crec = 0.; iwrk2d = 0; rwrk1d = 0.; cwrk1d = 0.; cwrk3d = 0.
!$OMP PARALLEL DO SCHEDULE(DYNAMIC,OMPCHUNK) DEFAULT(SHARED) PRIVATE(iws,iz,ith,ws)
        do iws=1, aws%n
           ith  = omp_get_thread_num() + 1
           ws   = aws%o + (iws-1)*aws%d
           indw = 1.5 + (ws-aw%o)/aw%d 

           if (iws/10*10==iws) write(0,*)"iws=", iws, ws, ith

           wfdr_slice(:,:,ith) = 0.
           do iz=az%n, 1, -1
              wfdr_slice(:,:,ith) = wfdr_slice(:,:,ith) + modl(:,:,iz)
              stat = wem_tapers_oper(wfdr_slice(:,:,ith), ax%n, ay%n)  
              stat = wem_wfdexp_oper(wfdr_slice(:,:,ith), vmod(:,:,iz), vref(:,iz), ws, +1,                          &
                     iwrk2d(:,:,ith), rwrk1d(:,ith), cwrk1d(:,ith), cwrk3d(:,:,:,ith), plan_2drF(ith), plan_2drI(ith))
           end do
           crec(:,:,iws) = wfdr_slice(:,:,ith)*hamw(indw)
        end do
!$OMP END PARALLEL DO
        call sep_write(crec(nxpad_beg+1:nxpad_beg+arx%n,nypad_beg+1:nypad_beg+ary%n,1:aws%n), crec_)
     end do
  end if

  stat = wem_velref_clos()
  stat = wem_wfdexp_clos()
  stat = wem_tapers_clos()

  call wem_fftw_clos(plan_2drF, plan_2drI, nth)

  deallocate(crec)
  deallocate(vmod, vref)
  deallocate(wfdr_slice)
  deallocate(tapx, tapy, hamw)
  deallocate(kx2, ky2)
  deallocate(plan_2drF, plan_2drI)


end program
