program PseSacGenMH
c Generate minimum-header SAC data files for PSE data
c ** 13-Mar-94 **
c f77 PseSacGenMH.f Ntimex.f lpedit.f spedit.f -lsac
c setenv SACAUX /usr/local/src/iris/sac10.6f/aux
c
	real data(65536)
	integer in(4864),lp(360,3),sp(2882),it(4)
	integer*2 exhdr(8)
	logical newfmt
	character cmp*3,ext*8,fname*16
	equivalence (exhdr(1),in(1))
c
c open input PSE data file in Exabyte format
	write(*,"('Open psedata.'$)")
	read(*,'(a)') ext
	fname='psedata.'//ext
	open(1,file=fname,access='direct',form='unformatted',recl=19456)
c name output file
	write(*,"('Output file psac'$)")
	read(*,'(a)') fname
10	write(*,"('lpx,lpy,lpz or spz? '$)")
	read(*,'(a)') cmp
	if (cmp.eq.'lpx') then
	 kp=1
	else if (cmp.eq.'lpy') then
	 kp=2
	else if (cmp.eq.'lpz') then
	 kp=3
	else if (cmp.eq.'spz') then
	 kp=0
	else
	 go to 10
	 end if
	i=2
	do while (fname(i:i).ne.' ')
	 i=i+1
	 end do
	fname='psac'//fname(1:i-1)//'.'//cmp
c specify data length and set sample interval
20	write(*,"('Record range: '$)")
	read(*,*) n1,n2
	nb=n2-n1+1
	if (kp.eq.0) then
	 npt=8640*nb
	 delt=.01886719
	else
	 delt=.1509375
	 end if
c process data
	j=1
	do nr=n1,n2
	 read(1,rec=nr) in
	 if (nr.eq.n1) then
	  nlb=(exhdr(6)+1)*3
	  if (kp.gt.0) npt=360*nlb*nb
	  if (npt.gt.65536) then
	   write(*,"('Data too long')")
	   go to 20
	   end if
	  lbs=4860/nlb
	  newfmt=exhdr(6).eq.1
	  ist=exhdr(2)
	  it(1)=exhdr(5)
	  call time03(in(5),it(2))
	  sec=it(4)/1000.
	  write(*,"('Data start at',3i5,f7.3)") (it(i),i=1,3),sec
	  end if
	 do lr=1,nlb
	  ip=(lr-1)*lbs+5
	  if (kp.eq.0) then
	   call spedit(in(ip),sp,ist)
	   if (nr.eq.n1.and.lr.eq.1) sp(3)=sp(4)
	   do i=3,2882
	    data(j)=sp(i)
	    j=j+1
	    end do
	  else
	   call lpedit(in(ip),lp,newfmt)
	   do i=1,360
	    data(j)=lp(i,kp)
	    j=j+1
	    end do
	   end if
	  end do
	 end do
c write a SAC data file
	call wsac1(fname,data,npt,0.,delt,nerr)
	end