* modified rdsetupsimulf from DWF. c rdsetupsimulf.yang.f do not do any corrections c except correcting c then make each station amplitude correction and the attenuation correction c This program is the part of simannerr that reads in the data and does the c Fourier analysis. It is separated so that simannerr# can run on katmai or c other machines that do not have sac. Writes out the data ready for use in c simannerr# c Pipe in data rdsetupsimul < eqlistper16 parameter (maxnfreq=1, maxnsta=1000, 1 maxpts = 100000, maxnodes = 700, maxevnts = 1000, 2 nusedsta = 1000, nusedtype = 2, ifreqbp = 11) real*4 staph(maxevnts,maxnsta,maxnfreq) real*4 staamp(maxevnts,maxnsta,maxnfreq) real*4 freq(maxnfreq) real*4 stadist(maxevnts,maxnsta), staazi(maxevnts,maxnsta) real*4 stacor(maxnsta), tempcor(maxnsta),geomsprd(maxnsta) real*4 bazi(maxevnts,maxnsta), stadelt(maxevnts,maxnsta) real*4 tdata(maxpts),beg(maxevnts,maxnsta),delt(maxnsta) real*4 stalat(maxevnts,maxnsta),stalon(maxevnts,maxnsta) real*4 attnfac(maxnfreq),fattnfac(10,2) real*4 staampcor(nusedsta,ifreqbp) real*4 adiattfactor1(ifreqbp) , adiattfactor2(ifreqbp) real*4 typeampcor(nusedtype,ifreqbp),typephcor(nusedtype,ifreqbp) real*4 avevel(ifreqbp),idfreq(ifreqbp) real*4 istatype(maxnsta) integer*4 nsta(maxevnts), dot, idot integer*4 nfreq,nstapts(maxnsta),idnum(maxevnts) integer*4 istanum(maxevnts,maxnsta),istacor(maxnsta), nevents integer*4 nevntsta(maxnsta) integer*4 findid integer*4 findstatype character*130 foutput, fn(maxevnts,maxnsta) , fsummary character dirpath1 *2, dirpath2 *37, evtdatatime *14,fnname *80 character*70 finvrsnodes, fftoutput, fvariance,fmaxavamp character*70 ftemp,fvelarea c character*14 staampcor character*44 dummy character*80 dummy1 character*3 idstaname pi = 3.1415928 convdeg = 3.1415928/180. circ = 6371.*3.1415928/180. twopi = 3.1415928*2. dirpath1 = "./" c dirpath2 = "earthquake_data/surfacewave/Rayleigh/" c ************************************************** c WARNING: The following two statements need to be switched depending on which c compiler is used c ************************************************ 210 format(a67, a3) c200 format(1H ,a44, i2) c200 format(a44, i2) c read in frequency and its corresponding attenuation coefficient ( data from c Brian Mitchell, 1995) c data(fattnfac(i,1), fattnfac(i,2), i=1,6)/0.05, 0.15e-3, 0.025, c * 0.15e-3, 0.01667, 0.15e-3, 0.01429, 0.1e-3, 0.0125, 0.05e-3, c * 0.0111, 0.05e-3/ c cattnfac=0.25e-3 c read list of files and frequencies to be analyzed and file to output results c Usually will pipe in data from some file like fasearrayinp read(*,*) nevents write(*,*) nevents nobs = 0 do iev = 1, nevents read(*,*) nsta(iev),idnum(iev) nobs = nobs+ 2*nsta(iev) do i = 1, nsta(iev) read(*, '(a)') fnname C evtdatatime = fnname(1:84) fn(iev,i) = dirpath1//fnname write(*,'(a)') fn(iev,i) enddo c read(*,'(a)') (fn(iev,i), i=1,nsta(iev)) c write(*,'(a)') (fn(iev,i), i=1,nsta(iev)) write(*,*) iev enddo c a lot of the following input is unneeded, just here to parallel simannerr# c input, so can use same fasearrayinp# or eqlistper# file for both read(*,*) nfreq write(*,*) nfreq read(*,*) (freq(j), j= 1, nfreq) read(*,*) foutput read(*,*) fsummary read(*,*) finvrsnodes read(*,*) fftoutput read(*,*) fvariance read(*,*) fmaxavamp read(*,*) ftemp read(*,*) fvelarea read(*,*) iterlimit, wlambda, dampvel,dampaniso write(*,*) iterlimit, wlambda, dampvel,dampaniso open(11, file = fftoutput) open(13, file = ftemp) open(34, file = 'bad_headdat.dat') idfreq(1) = 0.007 idfreq(2) = 0.008 idfreq(3) = 0.009 idfreq(4) = 0.010 idfreq(5) = 0.012 idfreq(6) = 0.015 idfreq(7) = 0.020 idfreq(8) = 0.025 idfreq(9) = 0.030 idfreq(10) = 0.035 idfreq(11) = 0.040 do i = 1, ifreqbp if ( freq(1) .eq. idfreq(i) ) then idnumfreq = i goto 102 endif enddo 102 write(*,*) idnumfreq c start input loop over events do iev = 1, nevents do ista = 1, nsta(iev) call rsac1(fn(iev,ista),tdata,nstapts(ista),beg(iev,ista), 1 delt(ista),maxpts,nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('b', beg(iev,ista),nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('dist',stadist(iev,ista),nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('az', staazi(iev,ista),nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('baz', bazi(iev,ista),nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('gcarc', stadelt(iev,ista), nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('stla', stalat(iev,ista),nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif call getfhv('stlo', stalon(iev,ista),nerr) if( nerr .ne. 0 ) then write(*,*) 'cannot read ', fn(iev,ista) stop endif geomsprd(ista) = sqrt(sin(stadelt(iev,ista)*convdeg)) do ifreq = 1, nfreq call frt(tdata, freq(ifreq), staamp(iev,ista,ifreq), 1 staph(iev,ista,ifreq), 1 nstapts(ista),delt(ista)) staamp(iev,ista,ifreq)= staamp(iev,ista,ifreq) 1 *geomsprd(ista) enddo enddo enddo c end of input loop over events c output loop over events do iev = 1, nevents write (11,*) iev, idnum(iev) write(*,*) iev,idnum(iev) do ista = 1, nsta(iev) write(11,*) beg(iev,ista) write (11,*) stadist(iev,ista),staazi(iev,ista), bazi(iev,ista), 1 stadelt(iev,ista), stalat(iev,ista), stalon(iev,ista) write(11,*) (staamp(iev,ista,ifreq),staph(iev,ista,ifreq), 1 ifreq=1,nfreq) enddo enddo c close(unit = 13) close(unit = 11) c close(unit = 12) end c---positive fourier transform c SUBROUTINE FRT(UP,FR,ARZ,PRZ,NALL,DELT) dimension UP(40000) DIMENSION W(40000) THETA=6.283185*FR*DELT C=COS(THETA)*2.0 NR1=1 NR2=NALL NDR1=NR2-1 W(1)=UP(NR2) W(2)=C*W(1)+UP(NDR1) NDATA=NR2-NR1+1 NTR1=NDATA-1 NTR2=NDATA-2 DO 97 I=3,NDATA I1=I-1 I2=I-2 NDRI=NR2-I+1 W(I)=C*W(I1)-W(I2)+UP(NDRI) 97 CONTINUE ZRE=(W(NDATA)-W(NTR2)+UP(NR1))*DELT/2.0 ZIM=W(NTR1)*SIN(THETA)*DELT CALL PHASE(ZRE,ZIM,PRZ) ARZ=SQRT(ZRE*ZRE+ZIM*ZIM) RETURN END SUBROUTINE PHASE(X,Y,PHI) IF(X) 21,20,22 20 IF(Y) 23,24,25 23 PHI=1.5*3.141592 GO TO 28 24 PHI=0.0 GO TO 28 25 PHI=0.5*3.141592 GO TO 28 21 PHI=ATAN(Y/X) +3.141592 GO TO 28 22 IF(Y) 26,27,27 27 PHI=ATAN(Y/X) GO TO 28 26 PHI=ATAN(Y/X)+2.0*3.141592 GO TO 28 28 CONTINUE PHI=PHI/6.283184 PHI=PHI-AINT(PHI) RETURN END integer function dot(file) character file*70 do 50 i=1,70 if(file(i:i).ne.'.') goto 50 dot=i-1 return 50 continue write(1,100) file 100 format(' no dots found in ',a70) dot = 0 return end integer function findid(inpstaname) parameter ( nsta = 41) character *3 inpstaname,staname cc 123 format( a3,i2) open(125,file = 'stationid.dat',status = 'old') do ista = 1, nsta read(125,*) staname, idsta, idummy cc write(*,*) staname,idsta if ( inpstaname .eq. staname) then findid = idsta goto 111 endif enddo 111 continue close(125) return end integer function findstatype(stanumber) parameter ( nsta =41) integer statype,stanumber character*3 dummy open(125,file = 'stationid.dat',status = 'old') do ista = 1, nsta read(125,*) dummy,idsta,statype if ( stanumber .eq.idsta ) then findstatype = statype goto 1111 endif enddo findstatype = 0 1111 continue close(125) return end