recfunk_pick_through.f -- read list of data files and retain only those which have a P-wave window marked for RF estimation

click title to see Google Drive directory JParkCodes

c program recfunk_pick_through
c version of recfunk_pick that reads seismic record start times from
c the sac header, and writes to in_recpick the filenames that have been picked
c 10/13/05 JJP
c input file: in_recpick_raw
c output file: outfile (to be renamed in_recpick)
c
c xf77 -o /park/backup/bin/recfunk_pick_through recfunk_pick_through.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/recfunk_pick_through recfunk_pick_through.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
......implicit real*4(a-h,o-z)
......real*8 ar,ai,el(12),aaa
......complex*8 zc,zero
......complex*8 afft,rf,rfs
......character*80 name,subtext,string
......character*28 title(2)
......character*10 cmps(3)
......character comp(3)
......common/nnaamm/iwflag
......common/npiprol/anpi
......common/stap2/tap(16400,16)
......common/taperz/ta(16400,8,3),tai(16400,8,3)
......common/stap0/ar(16400),ai(16400),bb(16400),time(99000)
......common/specs/afft(4100,8,3),rf(4100,2),sss(4100,6),crf(4100,2)
......common/saveit/rft(4100,799,2),crft(4100,2),rfs(2050,799,2),
x drfs(2050,799,2),bazs(799),epic(799),s2n(799),ip(799)
......common/saveit2/crfts(4100,799,2)
......common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
......common/header/ahead(158)
......common/distrib/xi(72,36),yi(72,36),sumij(72,36)
......common/chisq/chi(160000),enchi(160000),bbaz(160000),
x chim(200,2),enchim(200,2)
......dimension iah(158),chead(158),mmonth(12),fr(3),tt(3),ttt(3)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data comp/'r','t','z'/
......data cmps/'Radial ','Transverse','Vertical '/
......data title/'Radial Receiver Function ',
x 'Transverse Receiver Function'/
......con=180./3.14159265358979
......pi=3.14159265358979
......pih=pi/2.
......zero=cmplx(0.,0.)
......fmax=5.
......print *,'This code reads start of interval from sac header'
......print *,'it writes filenames of those files that have been picked'
c some default values (often superceded)
......fmax0=fmax
......fmin=0
......anpi=2.5
......nwin=3
......npad=16384
......nnf=npad/2
......do i=1,1000
...... crft(i,1)=0.
...... crft(i,2)=0.
......end do
......irec=0
......open(10,file='in_recpick_raw',form='formatted')
......open(11,file='in_recpick',form='formatted')
10 print *,'enter filename e.g. 1998.361.bhz'
......read(10,101) name
101 format(a)
......if(name(1:4).eq.'stop') go to 111
......subtext=name
c we assume that the last character is either 'r','t',or 'z'
......do i=80,1,-1
...... if(name(i:i).ne.' ') then
...... namlen=i
...... go to 98
...... endif
......end do
......print *,name
......print *,'no filename?'
......stop
98 continue
......do kcmp=3,3
...... name(namlen:namlen)=comp(kcmp)
...... print *,name
...... call sacread(name,a(1,kcmp),ierr)
...... if(ierr.eq.1) stop
...... dt=ahead(1)
...... tt(1)=0.
...... tt(2)=dt
...... tt(3)=dt
...... nscan=iah(80)
...... if(nscan.gt.99000) then
...... print *,'Careful! data length greater than 99000:', nscan
...... pause
...... endif
......end do
c naturally, you are in deep poop if dt is not the same for all records
......npre=10./dt
......npost=30./dt
......ntot=npre+npost
......ttt(1)=-npre*dt
......ttt(2)=dt
......ttt(3)=dt
......do i=1,ntot
...... time(i)=-ttt(1)-(i-1)*ttt(2)
......end do
......irec=irec+1
......if(irec.gt.799) then
...... print *,'max records exceeded'
...... stop
......endif
......baz=ahead(53)
......bazs(irec)=ahead(53)
......epic(irec)=ahead(54)
......print *,baz,'= back azimuth'
c ddf is cycles/sec
......ddf=1./(dt*npad)
......nf1=1
......nf2=fmax/ddf+1
......nf=nf2-nf1+1
......if(nf.le.0) stop
......fr(1)=0.
......fr(2)=ddf
......fr(3)=ddf
c FOR READING TIME INTERVAL FROM SAC HEADER
c ahead(9) is SAC's A parameter, for arrival time. Set by pressing A in plotpk
c......print *,'P window start time, duration (sec)'
c......read(10,*) pstart,postwin
......pstart=ahead(9)-ahead(6)-3.
......print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart
......if(ahead(9).gt.0.1) write(11,101) subtext
......go to 10
111 write(11,101) name
......close(10)
......close(11)
......stop
......end