rotate_sac.f -- Preprocessing of SAC data for receiver function analysis

c program rotate_sac_puteventinheader  (sac_rotate assumes the event was already inserted)
c 7/8/99 JJP -- adapted for JWEED output
c f77 -o /Users/jjpark/bin/rotate_sac rotate_sac.f /Users/jjpark/Ritz/jlib.a
c to read infiles1=filenames of data files (E,N,Z) for a single event
c cross-ref with file *events* of source locs and times
c writes a SAC macro to read three components of data, synchronize times, cut to equal length,
c put the event location in the header, and rotate the horizontals to radial and transverse.
c
c there must be three lines in infiles1 for every line in events, and the events must match the data
c take out the dots and replace with spaces before compiling
c
......character*80 name,nami,namo*15,dumf(3)*2,bigbanner*85,chh*51
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,nami(3),namo,name,dumf
c we read in data in scratch space
......common/sing2/a(240000)
......common/header/ahead(158)
......dimension iah(158),chead(158),mmonth(12)
......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
......print *,'open(7,file=infiles1) ! sacfilenames'
......print *,'open(18,file=events) ! source locations and times'
......open(7,file='infiles1',form='formatted') ! sacfilenames
......open(18,file='events',form='formatted') ! source locations and times
......open(9,file='outfile',form='formatted')
......write(9,101) 'qdp 2000'
10 print *,'starting again'
c......if(jjc.eq.1) go to 112
......read(18,101,end=110) name
......print *,name
c......go to 10
c 112 continue
......read(name,102) jy,jm,jd,jh,jmi,ssec
......print *, jy,jm,jd,jh,jmi,ssec
c......nami(1:38)=name(43:80)
c......read(name,105) slat,slon,sdep
......read(name(44:80),*) slat,slon,sdep
......print *, slat,slon,sdep
101 format(a)
103 format(7a)
102 format(19x,i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f6.3)
c here reads to column 45 -- if evdp can be read OK here, doesnt matter
c if the spacing of columns is perfect
105 format(22x,f8.2,f9.2,f4.0)
c 105 format(25x,3f8.2)
104 format(a,f10.3,a,f8.3,a,f10.2)
......leap=1-(mod(jy,4)+3)/4
......jd_julian=mmonth(jm)+jd
......if(jm.gt.2) jd_julian=jd_julian+leap
......write(9,103) 'cuterr fillz'
......sdep=sdep*1000.
......print *,slat,slon,sdep,' is depth in meters?'
......do k=1,3
...... read(7,101) nami(k)
......end do
......write(9,103) '/bin/rm ae an az'
......do k=1,3
...... name=nami(k)
...... ii=80
...... do while (name(ii:ii).eq.' ')
...... ii=ii-1
...... end do
c...... ii=ii-4 ! earlier SEED-output name convention
...... ii=ii-6 ! newer SEED-output name convention
...... print *,ii
...... if(name(ii:ii).eq.'Z') then
...... j=3
...... elseif(name(ii:ii).eq.'E'.or.name(ii:ii).eq.'2') then
...... j=1
...... elseif(name(ii:ii).eq.'N'.or.name(ii:ii).eq.'1') then
...... j=2
...... else
...... print *,'bad character: ',name(ii:ii)
...... endif
...... write(9,103) 'sc ln ',name(1:ii+7),dumf(j)
......end do
c LOVROK is the flag for allowing the user to overwrite a data file
......chh='synch;rmean;rtrend;p1;chnhdr LOVROK TRUE;write over'
......write(9,103)'read ae an az;',chh
......write(9,103) 'cut 0 600'
......write(9,103)'read;write over;cut off'
......do k=1,3
...... name=nami(k)
c read in the file to obtain its reference time
...... call sacread(name,a,ierr)
...... if(ierr.ne.0) go to 111 ! in case of mistyped filename
...... iy=iah(71)
...... id=iah(72)
...... ih=iah(73)
...... im=iah(74)
...... iss=iah(75)
...... ims=iah(76)
c compute origin time - year, month are assumed equal
...... torg=ssec-iss-float(ims)/1000.+
x......60.*(jmi-im+60*(jh-ih+24*(jd_julian-id)))
...... write(9,103) 'readhdr ',nami(k)
...... write(9,104) 'chnhdr evla ',slat,' evlo ',slon,' evdp ',sdep
...... write(9,104) 'chnhdr o ',torg
...... write(6,104) 'chnhdr o ',torg
...... write(9,103) 'writehdr'
......end do
......write(9,103) 'read ae an az'
......write(9,103) 'plot1'
......write(9,103) 'pause period 0.25'
......write(9,103) ' '
......write(9,103) 'read ae an'
......write(9,103) 'rotate to gcp reversed'
......write(9,103) 'read more az'
......write(9,103) 'plot1'
......write(9,103) 'pause period 0.25'
......write(9,103) ' '
......namo=name(1:15)
......write(9,103) 'write ',namo,'bhr ',namo,'bht ',namo,'bhz'
......go to 10
110 close(7)
......close(18)
......close(9)
......stop
111 print *,'error! ',name
......stop
......end