# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module kxmig_mod
  !! fk-fx migration
  !! receive wfld at (izs,iws)
  !! return wfld downward/upward continued
  !! needs Fop, FXop

  use util_mod
  use external_mod
  use down_mod
  use S_mod
  use fft_mod
  use taper_mod

  !---------------
  implicit none

  !! temporary interpolation arrays
  complex, dimension(:,:,:), allocatable, private :: tmp1,tmp2

contains
!----------------------------------------------------------------
  subroutine kxmig(izs,iws,ith,wfld_r,wfld_s,WCop,FKop,FXop)
    complex, dimension(:,:,:), allocatable :: wfld_s,wfld_r
    integer,external :: FKop,FXop
    integer,intent(in) :: izs, iws!,ith
    integer :: st,iw,ith



    interface !!-------------------------
       integer function WCop(iws,izs,ith,wfld,FKop,FXop) result(st)
         complex, dimension(:,:,:), allocatable ::wfld
         integer, intent(in) :: iws,izs,ith
         integer,external :: FKop,FXop
       end function WCop
    end interface !!---------------------



    ith=1
    if(down%sourcefn.or.down%passive) then !!D.C. "/homes/sep/claudio/research.ln/Others/ShotMig/kxmig.F90" WAVEFIELD
       call sellect_slow_comp(.true.) !!selects slowness from source field

!$OMP PARALLEL PRIVATE(ith)
       ith=1+omp_get_thread_num()
!$OMP DO PRIVATE(iw,st)

       do iw=1,down%nws
          st = WCop(izs,iw,ith,wfld_s,FKop,FXop)
       end do

!$OMP END DO
!$OMP END PARALLEL

    end if
    if(down%gathers) then !!D.C. RECEIVER WAVEFIELD
       call sellect_slow_comp(.false.)!!selects slowness from recvr field

!$OMP PARALLEL PRIVATE(ith)
       ith=1+omp_get_thread_num()
!$OMP DO PRIVATE(iw,st)

       do iw=1,down%nws
          st = WCop(izs,iw,ith,wfld_r,Fkop,Fxop)
       end do

!$OMP END DO
!$OMP END PARALLEL

    end if



  end subroutine kxmig


!----------------------------------------------------------------
!! N reference slownesses
  integer function kxmigN_init() result(st)
    integer :: er




    call init_wave_temps()
    if(down%antion) then
       call cos_taper_cmplx_allocate()
    end if



    st=OK
  end function kxmigN_init
!----------------------------------------------------------------
!! N reference slownesses SSR (WCop)
  integer function kxmigN_SSR(izs,iws,ith,wfld,FKop,FXop) result(st)
    complex, dimension(:,:,:) , allocatable :: wfld
    integer, intent(in) :: izs,iws,ith
    integer :: ierr,ivel
    real:: w



    interface !!----------------------------------------------
       integer function FKop(izs,iws,ith,ivel,wfld) result(st)
         integer, intent(in) :: izs,iws,ith,ivel
         complex, dimension(:,:,:), allocatable:: wfld
       end function FKop
       integer function FXop(izs,iws,ith,ivel,wfld) result(st)
         integer, intent(in) :: izs,iws,ith,ivel
         complex, dimension(:,:,:), allocatable:: wfld
       end function FXop
    end interface !!------------------------------------------



    !----------------------------------------------------------------
       call store_wavefield(wfld,izs,iws,ith) !STORE THE CURRENT WAVEFIELD in tmp2
       refvels:do ivel=1,slo%nrs(izs) !ref. vels
          call grab_part_wvfld(wfld ,ivel,izs,iws,ith) !GRAB PORTION OF THE WFLD FOR REF.V
          st = FXop(izs,iws,ith,ivel, wfld) !! f-x operator
          call fft(iws,ith,wfld)
          st = FKop(izs,iws,ith,ivel,wfld) !! f-k operator
          call add_wfld_extraps(wfld ,izs,iws,ith) !SUM OUR EFFORTS TOGETHER into tmp1
       end do refvels
       call grab_final_wfld(wfld,izs,iws,ith) !copy summed up field from tmp1

       if(down%antion.and.izs.eq.1) then
          w = down%bw%o + (iws-1)*down%aw%d
          call cos_taper_cmplx(wfld,iws,ith,w)
          if (iws==down%nws) write(run%io,*) '        *    Finished wavefield antialias taper.  '
       end if
       call ift(iws,ith,wfld)




    st=OK
  end function kxmigN_SSR
!----------------------------------------------------------------
!! multiple fk ops
!!
!! D(z)
!! |
!! v
!! T1
!!
!! D1 D2 D3 ...
!!
!! S1 S2 S3 ...
!! \ | /
!! T2
!! |
!! v
!! D(z+dz)
!---------------------------------------------
! Initialize storage for wavefields
subroutine init_wave_temps()
 integer ::er
 real::a
 !! allocate tmp arrays for interpolation of multiple fk ops
 if(down%nref .gt.1) then
    allocate(tmp1(down%amx%n,down%amy%n,run%o_nth),stat=er)
    if(er/=0) call seperr('cant allocate tmp1 in kxmig.F90')
    allocate(tmp2(down%amx%n,down%amy%n,run%o_nth),stat=er)
    if(er/=0) call seperr('cant allocate tmp2 in kxmig.F90')

    a=8*(size(tmp1)+size(tmp2))
    down%mem=down%mem+a



 end if
 end subroutine init_wave_temps
!----------------------------------------------------------------
!! tmp2 = original wfld before any propagation
 subroutine store_wavefield(wfld,izs,iws,ith)
   complex, dimension(:,:,:), allocatable:: wfld
   integer,intent(in)::izs,iws,ith
   integer :: i,j

   if(slo%nrs(izs) .gt.1) then

      do j=1,size(tmp2,2)
         do i=1,size(tmp2,1)
            tmp2(i,j,ith) = wfld(i,j,iws)
         end do
      end do

      do j=1,size(tmp1,2)
         do i=1,size(tmp1,1)
            tmp1(i,j,ith)=cmplx(0.,0.)
         end do
      end do
   end if

 end subroutine store_wavefield
!----------------------------------------------------------------
!! wfld seeded only at locations needing current ref. vel.
!! no need for shuffling if only 1 ref.vel. for this depth step
 subroutine grab_part_wvfld(wfld,ifk,izs,iws,ith)
   complex, dimension(:,:,:), allocatable :: wfld
   integer,intent(in)::ifk,izs,iws,ith
   integer :: i,j
   if(slo%nrs(izs) .gt.1) then
      do j=1,size(tmp2,2)
         do i=1,size(tmp2,1)
            wfld(i,j,iws)=cmplx(0.,0.)
         end do
      end do
      where(slo%si(:,:,izs) == ifk)
         wfld(:,:,iws) = tmp2(:,:,ith)
      end where
   end if
 end subroutine grab_part_wvfld
!------------------------------------------------------
!! current ref. vel. extrapolations added to tmp1
!! no need for shuffling if only 1 ref.vel. for this depth step
 subroutine add_wfld_extraps(wfld,izs,iws,ith)
   complex, dimension(:,:,:), allocatable :: wfld
   integer :: i,j
   integer,intent(in)::izs,iws,ith
   if(slo%nrs(izs) .gt.1) then
      do j=1,size(tmp2,2)
         do i=1,size(tmp2,1)
            tmp1(i,j,ith) = tmp1(i,j,ith) + wfld(i,j,iws)
         end do
      end do
   end if
 end subroutine add_wfld_extraps
!----------------------------------------------------------------
!! wfld = sum of extrapolations with each ref. vel.
!! no need for shuffling if only 1 ref.vel. for this depth step
 subroutine grab_final_wfld(wfld,izs,iws,ith)
   complex, dimension(:,:,:), allocatable:: wfld
   integer ::i,j
   integer,intent(in):: iws,ith,izs
   if(slo%nrs(izs) .gt.1) then
      do j=1,size(tmp2,2)
         do i=1,size(tmp2,1)
            wfld(i,j,iws) = tmp1(i,j,ith)
         end do
      end do
   end if
 end subroutine grab_final_wfld
!----------------------------------------------------------------
! subroutine kxmig_model(izs,iws,ith,wfld_in,WCop,FKop,FXop)
! complex, dimension(:,:,:), allocatable:: wfld_in
! integer,external :: FKop,FXop
! integer,intent(in) :: izs,iws!,ith
! integer :: st,iw,ith
!#ifdef DEBUG
! character(len=128), parameter :: unit='kxmig'
!#endif
! interface !!-------------------------
! integer function WCop(izs,iws,ith,wfld,FKop,FXop) result(st)
! complex, dimension(:,:,:), allocatable ::wfld
! integer, intent(in) :: izs,iws,ith
! integer,external :: FKop,FXop
! end function WCop
! end interface !!---------------------
!#ifdef DEBUG
! if(run%depar) call in(unit)
!#endif
! ith=1
! if(down%datum) then
! !!D.C.RECEIVER WAVEFIELD
! call sellect_slow_comp(.false.) !!selects slowness from source field
! else
! !!D.C. "/homes/sep/claudio/research.ln/Others/ShotMig/kxmig.F90" WAVEFIELD
! call sellect_slow_comp(.true.) !!selects slowness from source field
! end if
!
!#ifdef 1
!!$OMP PARALLEL PRIVATE(ith)
! ith=1+omp_get_thread_num()
!!$OMP DO PRIVATE(iw)
!#endif
! do iw=1,down%nws
! st = WCop(izs,iw,ith,wfld_in,FKop,FXop)
! end do
!#ifdef 1
!!$OMP END DO
!!$OMP END PARALLEL
!#endif
!
!#ifdef DEBUG
! if(run%depar) call out(unit)
!#endif
! end subroutine kxmig_model
!----------------------------------------------------------------
end module kxmig_mod
