module wei_image
  !! WEI imaging condition

  use wei_util
  use wei_process
  !---------------
  implicit none

  complex, dimension(:,:,:,:,:), pointer, private :: imwf
contains
  !----------------------------------------------------------------  
  subroutine weiimage_init()
    integer :: er
#ifdef DEBUG
    character(len=128), parameter :: unit='weiimage_init'
    if(pro%debug) call in(unit)
#endif

    !! imaging wavefield
    allocate(imwf(pro%amx%n,pro%amy%n,pro%ahx%n,pro%ahy%n,pro%nth),stat=er)
    if(er/=0) call seperr('cant allocate imwf')

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weiimage_init
  !----------------------------------------------------------------  
  subroutine weiimgW2R(izs,iws,ith,image,wfld,IGop)
    integer, intent(in) :: izs,iws,ith
    complex, dimension(:,:,:,:,:), pointer :: wfld, image
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weiimgW2R'
#endif
    interface !!-------------------------
       integer function IGop(wfld,iws,ith) result(st)
         complex, dimension(:,:,:,:,:), pointer :: wfld
         integer, intent(in) :: iws,ith
       end function IGop
    end interface !!---------------------
    integer :: i1,i2,i3,i4
#ifdef DEBUG
    if(pro%depar) call in(unit)
#endif

    imwf(:,:,:,:,ith) = wfld(:,:,:,:,iws)
    st = IGop(imwf,iws,ith) !! CIG h->h or h->ph

    do i4=1,size(image,4)
    do i3=1,size(image,3)
    do i2=1,size(image,2)
    do i1=1,size(image,1)
    !$OMP ATOMIC
    image(i1,i2,i3,i4,izs) = &
    image(i1,i2,i3,i4,izs) + &
    imwf (i1,i2,i3,i4,ith)/pro%aw%n
   end do
   end do
   end do
   end do


#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
  end subroutine weiimgW2R
  !----------------------------------------------------------------
  subroutine weiimgR2W(izs,iws,ith,image,wfld,IGop)
    integer, intent(in) :: izs,iws,ith
    complex, dimension(:,:,:,:,:), pointer :: wfld, image
    integer :: st
#ifdef DEBUG
    character(len=128), parameter :: unit='weiimgR2W'
#endif
    interface !!-------------------------
       integer function IGop(wfld,iws,ith) result(st)
         complex, dimension(:,:,:,:,:), pointer :: wfld
         integer, intent(in) :: iws,ith
       end function IGop
    end interface !!---------------------
#ifdef DEBUG
    if(pro%depar) call in(unit)
#endif

    imwf(:,:,:,:,ith) = image(:,:,:,:,izs)
    st = IGop(imwf,iws,ith) !! CIG: h->h or ph->h


    wfld(:,:,:,:,iws) = &
    wfld(:,:,:,:,iws) + &
    imwf(:,:,:,:,ith)/pro%aw%n

#ifdef DEBUG
    if(pro%depar) call out(unit)
#endif
  end subroutine weiimgR2W
  !----------------------------------------------------------------
end module wei_image





