subroutines sacread and sacout -- SAC-binary I/O

These subroutines lack the neat feature of the new SAC release, which can detect and compensate for byte-swapped data. Beware the swapped byte! There are variants of the routines included that swap bytes, dogswork to be sure.

c......subroutine sacread_noswap(name,a,ierr)
......subroutine sacread(name,a,ierr)
c read a SAC-format file, dont swap the bytes
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......ierr=0
......open(8,file=name,access='direct',recl=512,err=999)
c read the 158-word header, can use its info
......read(8,rec=1,err=998)(ahead(i),i=1,128)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......irec=2
......read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......else......
...... print *,'reading greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......endif
......if(nrec.gt.3) then
...... do irec=3,nrec-1
...... read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128)
...... end do
......endif
c......print *,'rec=',nrec,' Im trying to read the last partial record'
......read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last)
......close(8)
......return
999 print *,'open error'
......ierr=1
......return
998 print *,'read error: reading',irec,' out of',nrec
......ierr=1
......print *,irec
......close(8)
......return
......end
......subroutine sacread_swap(name,a,ierr)
c......subroutine sacread(name,a,ierr)
c read a SAC-format file -- swap bytes for Intel mac for all numbers, not characters
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......ierr=0
......open(8,file=name,access='direct',recl=512,err=999)
c read the 158-word header, can use its info
......read(8,rec=1,err=998)(ahead(i),i=1,128)
c swap bytes
......do i=1,120
...... chd=chead(i)
........do j=1,4
........ chhead(j,i)=chd(5-j:5-j)
........end do
......end do
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......irec=2
......read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98)
c swap bytes
......do i=1,98
...... ahd=a(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........a(i)=ahhd
......end do
c......print *,'rec=2'
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......else......
...... print *,'reading greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......endif
......if(nrec.gt.3) then
...... do irec=3,nrec-1
c...... print *,'rec=',irec
...... read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128)
c swap bytes
...... do i=1,128
............ahd=a((irec-3)*128+98+i)
...... .....do j=1,4
..............chhd(j:j)=chd(5-j:5-j)
........ end do
........ a((irec-3)*128+98+i)=ahhd
...... end do
...... end do
......endif
......print *,'byteswapping read'
c this here is a big fat kluge --- I cant read the last partial record
c so I skip the last fragment and adjust the header to discard the points.
c......print *,'NOT READING ',last,' points in last record -- SACBUGFIX'
......read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last)
c swap bytes
......do i=1,last
...... ahd=a((nrec-3)*128+98+i)
...... .do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........a((nrec-3)*128+98+i)=ahhd
......end do
c......nscan=nscan-last
c......iah(80)=nscan
......close(8)
......return
999 print *,'open error'
......ierr=1
......return
998 print *,'read error: reading',irec,' out of',nrec
......ierr=1
......print *,irec
......close(8)
......return
......end
c......subroutine sacout_noswap(name,a)
......subroutine sacout(name,a)
c write a SAC-format file
......real*4 ahead
......character*80 name
......character*4 chead(158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead)
......open(8,file=name,access='direct',recl=512)
c write the 158-word header
......write(8,rec=1)(ahead(i),i=1,128)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
......else......
...... print *,'writing greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
......endif
......write(8,rec=2)(ahead(i+128),i=1,30),(a(j),j=1,98)
......do irec=3,nrec
...... write(8,rec=irec)(a((irec-3)*128+98+i),i=1,128)
......end do
......print *,nrec,' records written'
......close(8)
......return
......end
c......subroutine sacout(name,a)
......subroutine sacout_swap(name,a)
c write a SAC-format file
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158),output(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......open(8,file=name,access='direct',recl=512)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
c write the 158-word header
c swap bytes
......do i=1,120
...... chd=chead(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........output(i)=ahhd
......end do
......do i=121,128
...... output(i)=ahead(i)
......end do
......write(8,rec=1)(output(i),i=1,128)
......nscan=iah(80)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
......else......
...... print *,'writing greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
......endif
c swap bytes
......do i=1,30
...... output(i)=ahead(128+i)
......end do
......do i=1,98
...... ahd=a(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........output(30+i)=ahhd
......end do
......write(8,rec=2)(output(i),i=1,128)
......do irec=3,nrec
...... iii=(irec-3)*128+98
...... do i=1,128
...... ahd=a(iii+i)
........ do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........ end do
........ output(i)=ahhd
...... end do
...... write(8,rec=irec)(output(i),i=1,128)
......end do
......print *,nrec,'byteswapped records written'
......close(8)
......return
......end