	
	subroutine psv_output

	include 'common.h'

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

	character*80 file


	if(isnap<>0.and.it==1)then
	file=seisfile(1:lnblnk(seisfile))//'_t_snp'
	open(unit=98,file=file,status='unknown')
	file=seisfile(1:lnblnk(seisfile))//'_r_snp'
	open(unit=99,file=file,status='unknown')
	endif

	write(*,*)izrec
	if(mod(it,isamp)==0)then
	do i=1,nr
	iout=(i-1)*nt+it
	forall(ix=1:nx)tmp1(ix)=theta1(ix,izrec)
	forall(ix=1:nx)tmp2(ix)=u21(ix,izrec)
c	call rec_int(tmp1,nx,rec_loc(i)/180.*pi,tmp2,seis(1,iout))
	forall(ix=1:nx)tmp1(ix)=theta2(ix,izrec)
	forall(ix=1:nx)tmp2(ix)=u22(ix,izrec)
c	call rec_int(tmp1,nx,rec_loc(i)/180.*pi,tmp2,seis(2,iout))
	enddo
	endif

	if(mod(it,icheckpoint)==0)then
	write(*,*)' Flushing out seismograms ... '
	file=seisfile(1:lnblnk(seisfile))//'_t'
	open(unit=10,file=file)
	file=seisfile(1:lnblnk(seisfile))//'_r'
	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


c snap output

	if(isnap<>0.and.mod(it,isnap).eq.0)then
	write(*,*)' Outputting snap : ',it/isnap
	do ix=1,nx,2
	do iz=1,nz,2
	write(98,*)u21(ix,iz)
	write(99,*)u22(ix,iz)
	enddo	
	enddo
	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


