examine_sac.f -- write a SAC macro to visually assess your body-wave data before receiver function estimation

c program examine_sac
c 7/8/99 JJP - adapted for WEED output, after running event_prune & rotate_sac
c 4/1/02 altered to order the input filenames by epicentral distance
c...... to facilitate quality control (P and S wave pattern varies smoothly)
c
c xf77 -o /Users/jjpark/bin/examine_sac examine_sac.f /Users/jjpark/Ritz/jlib.a /Users/jjpark/Plotxy/plotlib.a
c to read filenames of data files (bhr,t,z) for a single event
c
......character*80 name,nami,namo*15,dumf(3)*2,bigbanner*85
c slat, slon are in radians, sdep is in km, slat=colatitude
c SAC headers are degrees and meters, respectively, and lat is lat
c ---> so must convert
......common/sour/slat,slon,sdep,sol(6),am0,strike,dip,slip
......common/ptitle/bigbanner,name,namo,dumf
c we read in data in scratch space
......common/sing2/a(240000)
......common/allfiles/nami(3,6000),gcarc(6000),gcord(6000),ip(6000)
......common/header/ahead(158)
......dimension iah(158),chead(158),mmonth(12),xx(3)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data dumf/'ae','an','az'/
......con=180./3.14159265358979
......open(7,file='infiles2',form='formatted') ! sacfiles 19*.bh[rtz] 20*.bh[rtz]
......open(18,file='outfile1',form='formatted') ! to be in_recfunk
......open(9,file='outfile',form='formatted') ! to be sac_examine
......open(13,file='outfile3',form='formatted') ! to be infiles3 (for pruning SKS data sets)
......write(18,103) '0.4'
......j=1
10 print *,'starting again'
......do k=1,3
...... read(7,101,end=998) nami(k,j)
......end do
......name=nami(1,j)
......print *,name
......call sacread(name,a,ierr)
......gcarc(j)=ahead(54)
......j=j+1
......go to 10
998 nevents=j-1
......close(7)
......xx(1)=1.
......xx(2)=1.
c order the events by epicentral distance
......call plotit(xx,gcarc,dum,nevents,'Events','event number',
x 'epicentral distance',1,0,0.05,2,21)
......do l=1,nevents
...... blob=200.
...... do m=1,nevents
...... if(blob.gt.gcarc(m)) then
............blob=gcarc(m)
............mm=m
...... endif
...... end do
...... ip(l)=mm
...... gcord(l)=gcarc(mm)
...... gcarc(mm)=200.
......end do
......call plotit(xx,gcord,dum,nevents,'Events','event number',
x 'epicentral distance',1,0,0.05,2,22)
......do j=1,nevents
...... name=nami(1,ip(j))
...... write(18,103) name(1:17),'?'
...... write(18,103)'117 75 75'
......end do
......write(18,103) 'stop'
......close(18)
......write(9,101) 'qdp 2000'
......do j=1,nevents
...... name=nami(1,ip(j))
........write(13,103) nami(1,ip(j))
........write(13,103) nami(2,ip(j))
........write(13,103) nami(2,ip(j))
...... write(9,103) 'read ',name(1:17),'?; chnhdr a 115 '
...... write(9,103)'rtrend; lp c 0.15 np 8; write lowr lowt lowz'
...... write(9,103) 'read ',name(1:17),'?'
...... write(9,103)'rtrend; hp c .5 np 8; write hir hit hiz'
...... write(9,103) 'readhdr lowz;listhdr gcarc baz evla evlo evdp ',
x 'kzdate'
c to pick arrival times for pulses, change hi?;plot1;pause
c to hi?;plotpk markall on; write over
...... write(9,103) 'read ',name(1:17),'? low? hi?;plot1;pause'
......end do
999 close(9)
......close(13)
101 format(a)
103 format(7a)
102 format(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f4.1)
105 format(22x,f8.2,1x,f9.2,1x,f4.0)
104 format(a,f10.3,a,f8.3,a,f10.2)
......print *,'to pick arrival times, change '
......print *,'to '
......stop
......end