	
	subroutine psv_read

	include 'common.h'

	real, dimension (nx) :: tmp1,tmp2
!hpf$	distribute (*) :: tmp1,tmp2

	character*80 file

	if(it==1)then
        file=seisfile(1:lnblnk(seisfile))//'_t'
        open(unit=10,file=file,form='unformatted')
        file=seisfile(1:lnblnk(seisfile))//'_r'
        open(unit=11,file=file,form='unformatted')
	endif

	if(it==1)then
	write(*,*)' Reading buffer at it = ',it
	read(10)seisx
	write(*,*)minval(seisx),' < seisx < ',maxval(seisx)
	write(*,*)minval(seisx(:,nbuf)),' < seisx < ',maxval(seisx(:,nbuf))
	read(11)seisz
	write(*,*)minval(seisz),' < seisz < ',maxval(seisz)
	write(*,*)minval(seisz(:,nbuf)),' < seisz < ',maxval(seisz(:,nbuf))
	endif

	do i=1,nr
	iout=(i-1)*nt+it
	tmp1=theta1(:,izrec)
	if(mod(it,nbuf)<>0)tmp2=seisx(:,mod(it,nbuf))
	if(mod(it,nbuf)==0)tmp2=seisx(:,nbuf)
c	write(*,*)minval(tmp2),' < tmp2 < ',maxval(tmp2)
	call rec_int(tmp1,nx,rec_loc(i)/180.*pi,tmp2,seis(1,iout))
c	write(*,*)seis(1,iout)
	tmp1=theta2(:,izrec)
	if(mod(it,nbuf)<>0)tmp2=seisz(:,mod(it,nbuf))
	if(mod(it,nbuf)==0)tmp2=seisz(:,nbuf)
c	write(*,*)minval(tmp2),' < tmp2 < ',maxval(tmp2)
	call rec_int(tmp1,nx,rec_loc(i)/180.*pi,tmp2,seis(2,iout))
c	write(*,*)seis(2,iout)
	enddo


        if(mod(it,nbuf)==0.and.it<nt)then
        write(*,*)' Reading buffer at it = ',it
        read(10)seisx
        write(*,*)minval(seisx),' < seisx < ',maxval(seisx)
        read(11)seisz
        write(*,*)minval(seisz),' < seisz < ',maxval(seisz)
        endif




	if(it==nt)then
	write(*,*)' Flushing out seismograms ... '
	file=seisfile(1:lnblnk(seisfile))//'_at'
	open(unit=10,file=file)
	file=seisfile(1:lnblnk(seisfile))//'_ar'
	open(unit=11,file=file)
	nrt=nt/isamp
	do i=1,nr*nrt
	write(10,*)seis(1,i)
	write(11,*)seis(2,i)
	enddo
	close(10)
	close(11)
	endif


	return
	end
	

c ---------------------------------------------

	subroutine  rec_int(x,nx,xr,f,out)

c this routine interpolates linearly a scalar field known at 
c x(0:nx) and y(0:ny) locations to an arbitrary location
c within xmax and ymax.
c The receiver locations are xr,  the scalar input
c field is plane and the function value at xr  is
c out.

	real x(nx)
	real xr,out
	real f(nx)

	do i=1,nx-1
	if(xr .ge. x(i) .and. xr .lt. x(i+1) ) then
	ilow=i
	iup =i+1
	dx=x(i+1)-x(i)
	hx=xr-x(i)
	endif
	enddo

	f1=f(ilow)
	f2=f(iup)

	out = f1 + (f2-f1)/dx*hx

	return
	end


