!!.. propagates the wavefield excited at a single point
module grn_kernel_mod
  implicit none

contains


  integer function grn_kernel_init() result(stat)
  end function


  integer function grn_kernel_oper() result(stat)
  end function


  subroutine grn_kernel_oneway()
    do iwb=1, wblock%nb
       aws%o = wblock%os(iwb); aws%d = wblock%ds; aws%n = wblock%blen(iwb)
        call wem_fio_csou(csou, dgrid, wblock%indx(iwb), aws%n, ish, csou_, "reed")
        call grn_wfs_init(wfds, csou, dgrid, pgrid)

!$OMP PARALLEL DO SCHEDULE(DYNAMIC,OMPCHUNK) DEFAULT(SHARED) PRIVATE(iws,izs,ith,ws)
          do iws=1, aws%n

             ith = omp_get_thread_num() + 1
             if (iws/50*50==iws) write(0,*)"iws=", iws, ith

             ws = aws%o + (iws-1)*aws%d
             wfds_slice(:,:,ith) = wfds(:,:,1,iws)
             wfdr_slice(:,:,ith) = wfdr(:,:,1,iws)
             do izs=1, pgrid%az%n
                stat = wem_wfdexp_oper(wfds_slice(:,:,ith), vmod(:,:,izs), vref(:,izs), ws, +1,                          &
                       iwrk2d(:,:,ith), rwrk1d(:,ith), cwrk1d(:,ith), cwrk3d(:,:,:,ith), plan_2dsF(ith), plan_2dsI(ith))
                stat = wem_wfdexp_oper(wfdr_slice(:,:,ith), vmod(:,:,izs), vref(:,izs), ws, -1,                          &
                       iwrk2d(:,:,ith), rwrk1d(:,ith), cwrk1d(:,ith), cwrk3d(:,:,:,ith), plan_2drF(ith), plan_2drI(ith))
                stat = wem_tapers_oper(wfds_slice(:,:,ith), pgrid%ax%n, pgrid%ay%n)
                stat = wem_tapers_oper(wfdr_slice(:,:,ith), pgrid%ax%n, pgrid%ay%n)
                wfds(:,:,izs,iws) = wfds_slice(:,:,ith)
                wfdr(:,:,izs,iws) = wfdr_slice(:,:,ith)
             end do

!$OMP CRITICAL
             if (igrid%wntimg) stat = wem_imgcon_oper(modl, wfds(:,:,:,iws), wfdr(:,:,:,iws), ws)
             if (hgrid%wnthes) stat = wem_hesapp_oper(hess, wfds(:,:,:,iws), wfdr(:,:,:,iws))
!$OMP END CRITICAL

          end do

  end subroutine

  subroutine grn_twoway_kernel()
  end subroutine

end module
