      USE DFLIB
      integer(2) n_of_args
	dimension trace(100000),detrended(100000),windowed(100000)
	dimension deglitched(100000),padded(100000)
	dimension respremoved_tr(100000,3)
	complex fourier(100000),filtered(100000)
      complex p(100),z(100),accresp, velresp, dspresp
	complex respremoved_sp(100000,3)
	character input*150,root*150,current*150,currdir*150
	character reslis*80(1000),tralis*80(1000),resdir*80,tradir*80
	character sta_path*150(1000),sta_nam*20(1000)
	character event*12,stime*12,staname*4,comp*3,c4*4,c26*26,c45*45
	common /par/igmtn,iglit,itrend,ifilt, iresp,taper,norder,  
     *	flcut, fhcut, ismooth, fbox, freq1, freq2, nfreq, damp
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	common /traceinfo/sr,nsamp,trace,event,stime,staname,comp
	common /stainfo/stalat,stalon,stah
	common /detrend/a,b,detrended

	n_of_args=nargs()
	if(n_of_args.gt.1)then
	call getarg(1,input)
	else
c	Read data file (in seed format):
	write(*,*)'Enter input file name: '
	read(*,*)input
	endif

c	Read parameter file:
	call readpar()

c	make input file name based on file start time:
	call runner
     *   ('rdseed -t -f '//input(1:index(input,' ')-1)//' >tmp.txt',0,1)
	open(unit=7,file='tmp.txt',status='old')
	call skipline(7,2)
	read(7,'(45a)')c45
	close(unit=7)
	call runner('del tmp.txt',0,1)
	call serialn(num,1)
	write(c4,'(a1,i3.3)')'.',num
	c26=c45(21:42)//c4
	c26(5:5)='.'
	c26(9:9)='.'
	c26(12:12)='.'
	c26(15:15)='.'
	call runner('copy '//input(1:index(input,' ')-1)//' '//c26,0,1)
	input=c26//'  '


c	Create event directory, copy event file in it, change to event directory,
c	run rdseed, and clean the copied seed data file:
	root=currdir(0,1)
	call mkdir('dir_'//input(1:index(input,'  ')-1),0,1)
	call choose_copy_path(root,root(1:index(root,'  ')-1)//'\'//'dir_'
     *//input(1:index(input,'  ')-1),input)
	call changedir(root(1:index(root,'  ')-1)//'\'//'dir_'//input,0,1)
	call runner
     *		('rdseed -R -d -o 6 -f '//input(1:index(input,' ')-1),0,1)
	call runner('del '//input,0,1)

c	Make list of SAC data and response files, open and read them:
	call runner('dirf *sac_asc',0,1)
	call runner('ren filenr.lis trace.lis',0,1)
	call runner('dirf resp.*',0,1)
	call runner('ren filenr.lis resps.lis',0,1)
	open(unit=1,file='resps.lis',status='old')
	open(unit=2,file='trace.lis',status='old')
	do 10 i=1,1000,1
	read(1,'(a80)',end=11,err=11)reslis(i)
10	continue
11	n_resp=i-2
	close(unit=1)
	do 20 i=1,1000,1
	read(2,'(a80)',end=21,err=21)tralis(i)
20	continue
21	n_tra=i-2
	close(unit=2)

c	Check if every trace has a response file (SIGNAL IF NOT):
	if(n_resp.ne.n_tra)write(*,*)'Error, event and trace mismatch!!!'

c	Generate station/component directories, and copy corresponding event and
c	response files in that directory.
	current=currdir(0,1)
	do 50 i=1,n_resp,1
	read(reslis(i)(8:),'(a)')resdir
	call mkdir(resdir(6:index(resdir,' ')-1),0,1)
c	call choose_copy_path(current,current(1:index(current,' ')-1)//'\'
c     *//resdir(6:index(resdir,' ')-1),'*'//resdir(6:index(resdir,' ')-1)
c     *//'*')
	call runner('copy '//current(1:index(current,' ')-1)//'\*'//resdir
     *(6:index(resdir,' ')-1)//'* '//current(1:index(current,' ')-1)//'\
     *'//resdir(6:index(resdir,' ')-1)//'\*.*',0,1)

c	USEFUL AND IMPORTANT LINES			These two lines provide directory names and paths
c	VVVVVVVVVVVVVVVVVVVVVVVVVV			to subfolders of recorded traces and spectra
	sta_nam(i)=resdir(6:index(resdir,' ')-1)
	sta_path(i)=current(1:index(current,' ')-1)//'\'//resdir(6:index(r
     *esdir,' ')-1)
c	^^^^^^^^^^^^^^^^^^^^^^^^^^
c	USEFUL AND IMPORTANT LINES
	call runner('del '//'*'//resdir(6:index(resdir,' ')-1)//'*',0,1)
50	continue

c	Start processing data of each station/component:
	do 100 i=1,n_resp,1
	call changedir(sta_path(i),0,1)
	call readresp('resp.'//sta_nam(i))
	call repair_trace(tralis(i)(8:))
	call readtrace(tralis(i)(8:))
	deglitched=trace												!KA20080508 Order of steps changed after discussion with Hadi
	call deglitcher(deglitched,nsamp,iglit)							!KA20080508 Order of steps changed after discussion with Hadi
	call detrender(deglitched,nsamp,itrend)							!KA20080508 Order of steps changed after discussion with Hadi
	call windower(detrended,nsamp,taper,windowed)					!KA20080508 Order of steps changed after discussion with Hadi
	call padder(windowed,nsamp,padded,npoints)						!KA20080508 Order of steps changed after discussion with Hadi
c	call detrender(trace,nsamp)										!KA20080508 Order of steps changed after discussion with Hadi
c	call windower(detrended,nsamp,taper,windowed)					!KA20080508 Order of steps changed after discussion with Hadi
c	deglitched=windowed												!KA20080508 Order of steps changed after discussion with Hadi
c	call deglitcher(deglitched,nsamp,iglit)							!KA20080508 Order of steps changed after discussion with Hadi
c	call padder(deglitched,nsamp,padded,npoints)					!KA20080508 Order of steps changed after discussion with Hadi

	fourier=cmplx(padded,0.)
	call fork(npoints,fourier,-1.)
	call filter(fourier,sr,npoints,ifilt,norder, flcut,fhcut,filtered)
	call respremover(filtered,sr,npoints,iresp,respremoved_sp)
	call tracemaker(respremoved_sp,npoints,respremoved_tr)
	call make_GSE_resp()
	call make_PITSA_PAZ()
	call make_SAC_PAZ()
	call make_FAP_resp()
	call make_PSA(sr,npoints,respremoved_tr(:,3),staname,comp)
	call writer(npoints,deglitched,windowed,fourier,filtered,
     *			respremoved_sp,respremoved_tr)
c	call writer(npoints,trace,deglitched,detrended,windowed,
c     *			fourier,filtered,respremoved_sp,respremoved_tr)
100	continue
	end




	subroutine writer(npoints,deglitched,windowed,fourier,filtered,
     *			respremoved_sp,respremoved_tr)
c	subroutine writer(nsamp,trace,deglitched,detrended,windowed,
c     *			fourier,filtered,respremoved_sp,respremoved_tr)
	dimension trace(100000),detrended(100000),windowed(100000)
	dimension deglitched(100000),respremoved_tr(100000,3)
	complex fourier(100000),filtered(100000)
	complex respremoved_sp(100000,3)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	complex p(100),z(100)
	character respname*30,event*12,stime*12,staname*4,comp*3,tsr*9
	common /par/igmtn,iglit,itrend,ifilt, iresp,taper,norder,  
     *	flcut, fhcut, ismooth, fbox, freq1, freq2, nfreq, damp
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	common /traceinfo/sr,nsamp,trace,event,stime,staname,comp
	common /stainfo/stalat,stalon,stah
	common /detrend/a,b,detrended
	common /writeinfo/iyr,idoy,mon,iday,respname
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	open(unit=3,file='result.txt',status='unknown')
	write(3,'(a19)')'RECORD INFORMATION:'
	write(3,20)iyr,mon,iday,stime
20	format('Beginning of record:  ',i4,2('/',i2.2),' ',a12)
	write(3,21)sr
21	format('Sampling rate:        ',f9.4)	
	write(3,22)nsamp,npoints
22	format('# of samples:         ',i6,/,'# after padding:      ',i6)
	write(3,'(a20)')'STATION INFORMATION:'
	write(3,23)staname,comp
23	format('Station name:         ',a4,/,'Component:            ',a3)
	write(3,231)stalat,stalon,stah
231	format('Station latitude:     ',f8.4,/,'Station longitude:    ',f9
     *.4,/,'Station height:       ',f8.2)
		if(iglit.eq.0)then
		write(3,24)
24		format('GLITCHES:',/,'Glitches removed:     No')
		else
		write(3,25)
25		format('GLITCHES:',/,'Glitches removed:     Yes')
		endif
	if(itrend.eq.0)then
	write(3,26)
26	format('TREND AND OFFSET:',/,'Trend/offset removed: No')
	else
	write(3,27)
27	format('TREND AND OFFSET:',/,'Trend/offset removed: Yes')
	write(3,'(a)')'A and B values in equation: y=Ax+B'
	write(3,28)a,b
28	format('A and B values:       ',2(g13.5,2x))
	endif
	write(3,29)taper
29	format('TAPERING:',/,'Tapering window type: Cosine',/,'Tapering wi
     *ndow ratio: ',f6.5)
		if(ifilt.eq.0)then
		write(3,30)
30		format('FILTER:',/,'Filter applied:       No')
		else
		write(3,31)
31		format('FILTER:',/,'Filter applied:       Yes')
		write(3,32)
32		format('Filter type:          Butterworth')
		write(3,33)norder
33		format('Filter order:         ',i2)
		write(3,34)flcut, fhcut
34		format('Frequency range:      ',f7.3,' to ',f7.3,' Hertz')
		endif
	write(3,'(a32)')'INSTRUMENT RESPONSE INFORMATION:'
	if(iresp.eq.0)then
	write(3,'(a)')'Instrument correction: Not applied'
	else
	write(3,'(a)')'Instrument correction: Applied'
	write(3,35)A0
35	format('Normalization factor: ',g13.5)
	write(3,36)f0
36	format('Normalization freq.:  ',g13.5)
	write(3,37)nzero
37	format('Number of zeros:      ',i2,' in:')
	do i=1,nzero
	write(3,'(2(2x,g13.5))')real(z(i)),imag(z(i))
	enddo	
	write(3,38)npole
38	format('Number of poles:      ',i2,' in:')
	do i=1,npole
	write(3,'(2(2x,g13.5))')real(p(i)),imag(p(i))
	enddo	
	write(3,39)sensitivity
39	format('Overall sensitivity:  ',g13.5)
	if(A0ok.eq.1)then
	write(3,'(a)')'Normalization factor: Matches'
	else
	write(3,'(a)')'Normalization factor: Does not match'
	endif
	endif
	write(3,40)
40	format('Calibration files in SAC, SEED, GSE2, and PITSA formats cr
     *eated in channel folder')
		write(3,'(a30)')'RESPONSE SPECTRUM INFORMATION:'
		write(3,'(a)')'PSA information file: PSA.out'
		if(ismooth.eq.0)then
		write(3,'(a)')'Response spectrum:    Not smoothed'
		else
		write(3,'(a)')'Response spectrum:    Smoothed'
		write(3,41)fbox
41		format('Smoothing box width:  ',f8.5,'Hz in each side')
		endif
	Write(3,'(a)')'END_HEADER'
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	write(3,61)
61	format('   Sample #   |   Original   |  Deglitched  |  Detrended   
     * |   Windowed   |<--   Fourier transform   -->| Fourier(Amp) |   F
     *iltered   | Instr. corr. Disp., Vel., Acc. Amp. Spect. | Instr. co
     *rr. Disp., Vel., Acc. time history|',/,'                  (Count) 
     *       (Count)        (Count)        (Count)        Re(Count*s)   
     *Im(Count*s)     (Count*s)      (Count*s)      (Meter*s)       (Met
     *er)       (Meter/s)       (Meter)       (Meter/s)      (CM/s^2)')
	do 10 i=1,npoints,1
	write(3,'(15(g13.6,2x))')i,trace(i),deglitched(i),detrended(i),
     *windowed(i),fourier(i),cabs(fourier(i)),cabs(filtered(i)),
     *(cabs(respremoved_sp(i,j)),j=1,3),(respremoved_tr(i,j),j=1,3)
10	continue
	close(unit=3)
	return
	end




	subroutine make_PSA(sr,npoints,a,staname,comp)
	dimension a(100000),b(100000),freq(1000),psa(1000),smooth(1000)
	character staname*4,comp*3
	common /par/igmtn,iglit,itrend,ifilt, iresp,taper,norder,  
     *	flcut, fhcut, ismooth, fbox, freq1, freq2, nfreq, damp
	open(unit=15,file='PSA.out',status='unknown')
	do 5 i=1,npoints,1
	b(i)=a(i)*100
	a(i)=b(i)
5	continue
      call rscomp(sr,npoints,b,freq,psa)
	call PSA_smoother(nfreq,freq,psa,ismooth,fbox,smooth)
	write(15,*)'Calculated PSA for the frequency range: '
	write(15,*)'Station:   ',staname
	write(15,*)'Component: ',comp
	write(15,*)freq1,' Hz and ',freq2,' Hz'
	write(15,*)nfreq,'frequencies distributed logarithmically'
	write(15,*)'Oscillator damping is',damp,'% of critical'
	write(15,'(a64)')'  Freq.(Hz)       Period(s)      PSA(cm/s^2)  Sm
     *ooth PSA(cm/s^2)'
	do 10 i=nfreq,1,-1
	write(15,*)freq(i),1/freq(i),psa(i),smooth(i)
10	continue
	close(unit=15)
	return
	end




	subroutine PSA_smoother(nfreq,freq,psa,ismooth,fbox,smooth)
	dimension freq(1000),psa(1000),smooth(1000)
	if(ismooth.eq.0)then
	smooth=psa
	return
	endif
	do 10 i=1,nfreq,1
		n=0
		sum=0.
		fl=freq(i)-fbox
		fh=freq(i)+fbox
		do 20 j=1,nfreq,1
			if(freq(j).ge.fl.and.freq(j).le.fh)then
			sum=sum+psa(j)
			n=n+1
			endif
20		continue
	smooth(i)=sum/float(n)
10	continue
	return
	end




      subroutine rscomp(sr,npoints,a,freq,psa)
c     calculates response spectrum from acceleration time history
      dimension freq(1000),Omega(1000),psa(1000),a(*)
c     damp is % critical damping. freq2 is maximum freq. to consider.
	common /par/igmtn,iglit,itrend,ifilt, iresp,taper,norder,  
     *	flcut, fhcut, ismooth, fbox, freq1, freq2, nfreq, damp
      pi=3.141593
      if (iabs(nfreq).gt.1000) then
		write (15,*)
		write (15,*) ' Nfreq cannot exceed 1000'
		return
      endif
c     Convert % damping to fraction
      beta=damp*0.01
c     Generate some logarithmically-spaced frequencies
      freq(1)=freq1
      omega(1)=2.*pi*freq(1)
      if (nfreq.ge.2) then
		frinc=alog(freq2/freq1)/(nfreq-1)
		do 120 k=2,nfreq
			freq(k)=freq1*exp((k-1)*frinc)
			omega(k)=2.*pi*freq(k)
120		continue
      endif
c     Call new response spectrum routine for each desired freq
      do 200 k=1,nfreq
		call sdcomp(a,npoints,omega(k),beta,1/sr,sd)
		psa(k)=sd*omega(k)*omega(k)
200   continue
      return
      end




      subroutine sdcomp(accg,na,omn,beta,dt,sd)
c This is a modified version of "Quake.For", written by
c Stavros A. Anagnostopoulos, Oct. 1986.  The formulation is that of
c Nigam and Jennings (BSSA, v. 59, 909-922, 1969).  This modification 
c eliminates the computation of the relative velocity and absolute 
c acceleration; it returns only the relative displacement.  
c Dates: 05/06/95 - Modified by David M. Boore
c This subroutine was implemented in agram developed by Gail Atkinson
      dimension accg(*)
      omt=omn*dt
      d2=1-beta*beta
      d2=sqrt(d2)
      bom=beta*omn
      d3=2.*bom
      omd=omn*d2
      om2=omn*omn
      omdt=omd*dt
      c1=1./om2
      c2=2.*beta/(om2*omt)
      c3=c1+c2
      c4=1./(omn*omt)
      ss=sin(omdt)
      cc=cos(omdt)
      bomt=beta*omt
      ee=exp(-bomt)
      ss=ss*ee
      cc=cc*ee
      s1=ss/omd
      s2=s1*bom
      s3=s2+cc
      a11=s3
      a12=s1
      a21=-om2*s1
      a22=cc-s2
      s4=c4*(1.-s3)
      s5=s1*c4+c2
      b11=s3*c3-s5
      b12=-c2*s3+s5-c1
      b21=-s1+s4
      b22=-s4
      sd=0.
      n1=na-1
      y=0.
      ydot=0.
      do 1 i=1,n1
      y1=a11*y+a12*ydot+b11*accg(i)+b12*accg(i+1)
      ydot=a21*y+a22*ydot+b21*accg(i)+b22*accg(i+1)
c next two lines have been added for hardware portability
      if (y1.lt.1.e-30.and.y1.gt.0.) y1=0.
      if (ydot.lt.1.e-30.and.ydot.gt.0.) ydot=0.
      y=y1
      z=abs(y)
      z1=abs(ydot)
      if (z.gt.sd) sd=z
1     continue
      return
      end




	subroutine make_FAP_resp()
      complex p(100),z(100),accresp,velresp,dspresp
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	open(unit=12,file='Response.FAP',status='unknown')
	write(12,20)
20	format(' Freq.(Hz)     Disp.Ampl.    Disp.Phase     Vel.Ampl.     
     *Vel.Phase     Acc.Ampl.    Acc.Phase')	
	do 10 i=0,200,1
	f=10.**(-2.+i*0.02)
	call response(A0,f,nzero,npole,z,p,sensitivity,accresp, 
     *			 velresp,dspresp)
	write(12,21)f,cabs(dspresp),phase(dspresp),cabs(velresp),
     *phase(velresp),cabs(accresp),phase(accresp)
21	format(7(1x,g12.5,1x))
10	continue
	close(unit=12)
	return
	end




	subroutine make_SAC_PAZ()
      complex p(100),z(100)
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	open(unit=14,file='SAC.PAZ',status='unknown')
		if(nzero.lt.10)then
		write(14,141)nzero
141		format('ZEROS ',i1)
		else
		write(14,142)nzero
142		format('ZEROS ',i2)
		endif
	if(nzero.ge.1)then
	do 10 i=1,nzero,1
	write(14,'(2(f11.4,2x))')real(z(i)),imag(z(i))
10	continue
	endif
		if(npole.lt.10)then
		write(14,143)npole
143		format('POLES ',i1)
		else
		write(14,144)npole
144		format('POLES ',i2)
		endif
	if(npole.ge.1)then
	do 20 i=1,npole,1
	write(14,'(2(f11.4,2x))')real(p(i)),imag(p(i))
20	continue
	endif
		write(14,145)sensitivity*A0
145		format('CONSTANT ',g13.5)
		close(unit=14)
	return
	end




	subroutine make_PITSA_PAZ()
      complex p(100),z(100)
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	open(unit=13,file='PITSA.PAZ',status='unknown')
	write(13,'(a34)')'CAL1                           PAZ'
		if(npole.lt.10)then
		write(13,'(i1)')npole
		else
		write(13,'(i2)')npole
		endif
	if(npole.ge.1)then
	do 10 i=1,npole,1
	write(13,'(2(g8.3E1))')real(p(i)),imag(p(i))
10	continue
	endif
		if(nzero.lt.10)then
		write(13,'(i1)')nzero
		else
		write(13,'(i2)')nzero
		endif
	if(nzero.ge.1)then
	do 20 i=1,nzero,1
	write(13,'(2(g8.3E1))')real(z(i)),imag(z(i))
20	continue
	endif
		write(13,'(g8.3)')sensitivity*A0
		close(unit=13)
	return
	end




	function phase(x)
	complex x
	if(real(x).gt.0.)then
	phase=atan(imag(x)/real(x))
	elseif(real(x).eq.0.)then
		if(imag(x).gt.0.)phase= 2*atan(1.)
		if(imag(x).lt.0.)phase=-2*atan(1.)
	elseif(real(x).lt.0.)then
		if(imag(x).gt.0.)phase=atan(imag(x)/real(x))+4*atan(1.)
		if(imag(x).eq.0.)phase=4*atan(1.)
		if(imag(x).lt.0.)phase=atan(imag(x)/real(x))-4*atan(1.)
	endif
	phase=phase*45./atan(1.)
	return
	end




	subroutine make_GSE_resp()
	dimension trace(100000)
      complex p(100),z(100)
	character event*12,stime*12,staname*4,comp*3,respname*30
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	common /traceinfo/sr,nsamp,trace,event,stime,staname,comp
	common /writeinfo/iyr,idoy,mon,iday,respname
	PI=4*atan(1.)
	read(event(1:4),'(i4)')iyr
	read(event(6:8),'(i3)')idoy
	call doy2md(iyr,idoy,mon,iday)
	respname(1:4)=staname
	do 10 i=1,4,1
	if(respname(i:i).eq.' ')respname(i:i)='_'
10	continue
	respname(5:5)='_'
	respname(6:7)=comp(1:2)
	respname(8:8)='_'
	respname(9:9)=comp(3:3)
	respname(10:10)='.'
	respname(11:14)=event(1:4)
	write(respname(15:29),15)mon,iday
15	format('-',i2.2,'-',i2.2,'-0000_GSE')
	sensGSE=10.**9/(sensitivity*2*PI*f0)
	sfactor=A0/10.**9
	open(unit=11,file=respname,status='unknown')
	write(11,21)staname,comp,sensGSE,1/f0,iyr,mon,iday
21	format('CAL2 ',a5,1x,a3,13x,e10.2,1x,f7.3,18x,i4,'/',i2.2,'/',i2.2
     *,' 00:00')
c21	format('CAL2 ',a5,1x,a3,13x,e15.8,1x,f7.3,13x,i4,'/',i2.2,'/',i2.2
c     *,' 00:00')
	write(11,22)sfactor,npole,nzero+1
22	format('PAZ2  1 V ',e15.8,15x,i3,1x,i3,' Laplace transform')
	do 100 i=1,npole,1
	write(11,23)real(p(i)),imag(p(i))
100	continue
	do 110 i=1,nzero,1
	write(11,23)real(z(i)),imag(z(i))
110	continue
	write(11,23)0.,0.
23	format(2(1x,e15.8))
	write(11,24)sensitivity,sr
24	format('DIG2  2 ',e15.8,1x,f11.5)
	close(unit=11)
	return
	end




	subroutine doy2md(iyr,idoy,mon,iday)
	dimension imonth(12)
	data (imonth(i),i=1,12)/31,28,31,30,31,30,31,31,30,31,30,31/
	if(mod(iyr,4).ne.0) goto 10
	imonth(2)=29
	if(mod(iyr,100).eq.0)then
	imonth(2)=28		
		if(mod(iyr,400).eq.0)then
		imonth(2)=29
		endif
	endif
10	iday=idoy
	do 20 i=1,12,1
		if(iday-imonth(i).le.0)then
		mon=i
		goto 30
		endif
	iday=iday-imonth(i)
20	continue
30	return
	end




	subroutine tracemaker(respremoved_sp,npoints,respremoved_tr)
	dimension respremoved_tr(100000,3)
	complex respremoved_sp(100000,3),temp(100000)
	do 10 i=1,3,1
	temp=respremoved_sp(:,i)
	call fork(npoints,temp,+1.)
	respremoved_tr(:,i)=real(temp)
10	continue
	return
	end




	subroutine respremover(filtered,sr,npoints,iresp,respremoved_sp)
	complex filtered(100000),respremoved_sp(100000,3),S
      complex p(100),z(100),accresp, velresp, dspresp
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	PI=4.*atan(1.)
	ncalc=npoints/2+1
	W0=2.*PI*sr/npoints
	if(iresp.eq.0)then
	do 10 i=1,ncalc,1
	W=(i-1)*W0
	S=cmplx(0.,W)
		if(i.ne.1)then
			respremoved_sp(i,1)=filtered(i)/S
			respremoved_sp(npoints+2-i,1)=respremoved_sp(i,1)
		else
			respremoved_sp(i,1)=filtered(i)/cmplx(0.,W0/100)
		endif
			respremoved_sp(i,2)=filtered(i)
			respremoved_sp(i,3)=filtered(i)*S
		if(i.gt.1)then
			respremoved_sp(npoints+2-i,1)=conjg(respremoved_sp(i,1))
			respremoved_sp(npoints+2-i,2)=conjg(respremoved_sp(i,2))
			respremoved_sp(npoints+2-i,3)=conjg(respremoved_sp(i,3))
		endif
10	continue
	else
	do 20 i=1,ncalc,1
	f=(i-1)*sr/npoints
	call response(A0,f,nzero,npole,z,p,sensitivity,accresp,velresp,
     *			  dspresp)
		respremoved_sp(i,1)=filtered(i)/dspresp
		respremoved_sp(i,2)=filtered(i)/velresp
		respremoved_sp(i,3)=filtered(i)/accresp
		if(i.gt.1)then
		respremoved_sp(npoints+2-i,1)=conjg(respremoved_sp(i,1))
		respremoved_sp(npoints+2-i,2)=conjg(respremoved_sp(i,2))
		respremoved_sp(npoints+2-i,3)=conjg(respremoved_sp(i,3))
		endif
20	continue
	endif
	return
	end




	subroutine filter(fourier,sr,npoints,ifilt,
     *			norder, flcut,fhcut,filtered)
	complex fourier(100000),filtered(100000)
	if(ifilt.eq.0)then
	filtered=fourier
	return
	endif
	ncalc=npoints/2+1
	do 10 i=1,ncalc,1
	f=(i-1)*sr/npoints
        	filtered(i) = cmplx(buttrlcf(f, flcut, norder),0)*fourier(i)
        	filtered(i) = cmplx(buttrhcf(f, fhcut, norder),0)*filtered(i)
		if(i.eq.1)goto 10
        	filtered(npoints+2-i) = conjg(filtered(i))
10	continue
	return
	end

      function buttrlcf(f, fcut, norder)
c**     norder Butterworth low-cut filter response from AGRAM
c**     taken from Boore, 1996
      real buttrlcf
      buttrlcf = 1.
      if (fcut .eq. 0.) return
      buttrlcf = 0.
      if (f .eq. 0.) return
      buttrlcf = 1./ (1. + (fcut/f)**(2*norder))
      return
      end

      function buttrhcf(f, fcut, norder)
c**     norder Butterworth high-cut filter response from AGRAM
c**     taken from Boore, 1996
      real buttrhcf
      buttrhcf = 1.
      if (fcut .eq. 0.) return
      buttrhcf = 0.
      if (f .eq. 0.) return
      buttrhcf = 1./ (1. + (f/fcut)**(2*norder))
      return
      end




      subroutine fork(lx,cx,signi)
c     fast fourier transform routine from Dave Boore.
c     result of sequence from time to freq to time requires no scaling.
      complex cx,carg,cexp,cw,ctemp
      dimension cx(lx)
      j=1
      sc=sqrt(1./lx)
c      write(*,*)'***** FAST FOURIER TRANSFORM CALLED *****'
      do 5 i=1,lx
          if(i .gt. j) go to 2
          ctemp = cx(j)*sc
          cx(j) = cx(i)*sc
          cx(i) = ctemp
2         m = lx/2
3         if(j .le. m) go to 5
            j = j - m
            m = m / 2
            if(m .ge. 1) go to 3
5     j = j + m
      l = 1
c      write(*,*)'First step completed'
6     istep = 2 * l
      temp = 3.14159265 * signi/l
      do 8 m = 1, l
          carg=(0., 1.) * temp * (m-1)
          cw=cexp(carg)
          do 8 i = m, lx, istep
               ctemp = cw *cx(i+l)
               cx(i+l) = cx(i) - ctemp
8     cx(i) = cx(i) + ctemp
      l = istep
      if(l .lt. lx) go to 6
c      write(*,*)'***** FAST FOURIER TRANSFORM DONE *****'
9     return
      end




	subroutine padder(windowed,nsamp,padded,npoints)
	dimension windowed(100000),padded(100000)
	padded=0.
	do 10 i=1,nsamp,1
	padded(i)=windowed(i)
10	continue
	do 20 i=1,16,1
	if(2**i.ge.nsamp)then
	npoints=2**i
	goto 50
	endif
20	continue
50	return
	end




      subroutine deglitcher(x,npts,iglit)
c
c     Removes glitches from array x and makes corresponding
c     correction to valmax.
c     Anything >10* the (log) avg of a running 20-pt. amplitude 
c     average is assumed to be a glitch.  Replaced with avg of 
c     neighbour values.
c     Constant amplitude steps (more than 20 identical values
c     in a row) are assumed to be glitches.  Replaced with 0
c
c     G.M. Atkinson, Oct. 1990
c     Revised Aug. 1991
c
      dimension x(100000)
c      valmax = 0.
c
c     First remove constant amplitude steps.
c     If next 20 values are exactly the same as this one,
c     we assume its an erroneous step.
c
	if(iglit.eq.0)goto 105

      jstop = 0
      do 500 j=1,npts-20
        iflag = 1
        do 400 jj = 1,20
          if(x(j+jj) .ne. x(j)) iflag = 0
400     continue
			if (iflag .eq. 1) then
c			This is a step.
				jstop = j + 20
			endif
        if (jstop .ge. j) x(j) = x(j+20)
500   continue
c      write(*,*)'Finished looking for a Long Glitch'
c
        ntotal = 0
49      continue
c        write(*,*)'Starting glitch removal'
        runavg = 0.
        nspike = 0
c        write(*,*)'Calculating Run Average'
        do 50 jj = 1,21
          if(x(jj) .eq. 0.) x(jj)=1.
          runavg = runavg + alog10(abs(x(jj)))
50      continue
        runavg=runavg/21.
c      write(*,*)'Looking for Spikes'
      do 100 j=11,npts-11
           if(x(j-10) .eq. 0.) x(j-10)=1.
           if(x(j+11) .eq. 0.) x(j+11)=1.
          runavg =runavg - alog10(abs(x(j-10)))/21.
     *            + alog10(abs(x(j+11)))/21.
          spike = 10.* 10.**runavg
          if(abs(x(j)) .gt. spike) then
             nspike = nspike + 1
             x(j) = (x(j-1) + x(j+1))/2.
             if(abs(x(j)) .gt. spike)then
               x(j) = (x(j-2) + x(j+2))/2.
               if(abs(x(j)) .gt. spike) x(j) = 1.
             endif
           endif
c        if(abs(x(j)) .gt. valmax) valmax = abs(x(j))
100   continue
      ntotal = ntotal + nspike
      if (nspike .gt. 1) go to 49
105   return
      end




	subroutine windower(trace,nsamp,taper,windowed)
	dimension trace(100000),windowed(100000)
	do 10 i=1,nsamp,1
	call win(i,1,nsamp,int4(nsamp*taper),wind)
	windowed(i)=trace(i)*wind
10	continue
	return
	end



      subroutine win(i, nstart, nstop, ntaper, wind)
c     applies cosine tapered window.
c      unit amplitude assumed
c      written by D. M. Boore
c      latest revision: 9/26/95
      wind = 0.0
      if ( i .lt. nstart .or. i. gt. nstop) return
      wind = 1.0
      if ( i .ge. nstart+ntaper .and. i .le. nstop-ntaper ) return
      pi = 4.0 * atan(1.0)
      dum1 = (nstop+nstart)/2.0
      dum2 = (nstop-nstart-ntaper)/2.0
      wind = 0.5 * (1.0 - sin( pi*
     *  ( abs(float(i)-dum1) - dum2 ) /float(ntaper) ) )
      return
      end




	subroutine detrender(deglitched,nsamp,itrend)
	dimension deglitched(100000),detrended(100000)
	doubleprecision x11,x12,x21,x22,y11,y21
	common /detrend/a,b,detrended
	if(itrend.eq.0)then
	detrended=deglitched
	return
	endif
	x11=0.
	x12=real(nsamp)
	x21=0.
	x22=0.
	y11=0.
	y21=0.
	do 10 i=1,nsamp,1
	x11=x11+i
	x21=x21+i**2
	y11=y11+deglitched(i)
	y21=y21+i*deglitched(i)
10	continue
	x22=x11
	a=(x22*y11-x12*y21)/(x11*x22-x12*x21)
	b=(x11*y21-x21*y11)/(x11*x22-x12*x21)
	do 20 i=1,nsamp,1
	detrended(i)=deglitched(i)-(i*a+b)
20	continue
	return
	end




	subroutine response(A0,f,nzero,npole,z,p,sensitivity,
     *					accresp, velresp, dspresp)
	complex z(100), p(100),S
	complex accresp, velresp, dspresp, Xnum, Xdenom
	PI=4.*atan(1.)
	if(f.eq.0.)f=0.01
      W=2.*PI*f
      S=CMPLX(0.0,W)
      Xnum = CMPLX(A0*sensitivity,0.)
	do iz=1, nzero
	  Xnum = Xnum*(S-z(iz))
	enddo 
	Xdenom = CMPLX(1.,0.)
	do ip = 1, npole
	  Xdenom = Xdenom*(S-p(ip))
	enddo
      velresp= Xnum/Xdenom
      dspresp=VELRESP*S
      accresp=VELRESP/S
	return
	end




	subroutine readresp(fname)
      complex p(100),z(100),accresp, velresp, dspresp
	character fname*(*),respread*80(3000)
	common /respinfo/A0,f0,nzero,npole,p,z,sensitivity,A0ok
	open(unit=4,file=fname,status='unknown')
	do 10 i=1,3000,1
	read(4,'(a80)',end=11,err=11)respread(i)
10	continue
11	nline=i-1
	close(unit=4)
	istage1=0
	istageL=0
	iZeroCount=0
	iPoleCount=0
	A0ok=1.
	do 20 i=1,nline,1
	if(respread(i)(1:7).eq.'B060F03')then
	read(respread(i)(52:),'(i6)')nstage
	endif
	if(respread(i)(1:7).eq.'B060F04')then
	read(respread(i)(52:),'(i6)')n_curr_stage
	endif
	if(n_curr_stage.eq.1)istage1=1
	if(n_curr_stage.ne.1)istage1=0
	if(n_curr_stage.eq.0)istageL=1
	if(istage1.eq.1.and.respread(i)(1:7).eq.'B043F08')
     *read(respread(i)(52:),*)A0
	if(istage1.eq.1.and.respread(i)(1:7).eq.'B043F09')
     *read(respread(i)(52:),*)f0
	if(istage1.eq.1.and.respread(i)(1:7).eq.'B043F10')
     *read(respread(i)(52:),*)nzero
	if(istage1.eq.1.and.respread(i)(1:7).eq.'B043F15')
     *read(respread(i)(52:),*)npole
	if(istage1.eq.1.and.respread(i)(1:7).eq.'B043F11')then
	iZeroCount=iZeroCount+1
	read(respread(i)(17:),*)z1,z2
	z(iZeroCount)=cmplx(z1,z2)
	endif
	if(istage1.eq.1.and.respread(i)(1:7).eq.'B043F16')then
	iPoleCount=iPoleCount+1
	read(respread(i)(17:),*)p1,p2
	p(iPoleCount)=cmplx(p1,p2)
	endif
	if(istageL.eq.1.and.respread(i)(1:7).eq.'B048F05')
     *read(respread(i)(52:),*)sensitivity
20	continue
	call response(A0,f0,nzero,npole,z,p,sensitivity,
     *              accresp, velresp, dspresp)
	if(abs(cabs(velresp/sensitivity)-1).gt.0.0001)A0ok=0.
	return
	end




	subroutine readtrace(fname)
	dimension trace(100000)
	character fname*(*),event*12,stime*12,staname*4,comp*3
	common /traceinfo/sr,nsamp,trace,event,stime,staname,comp
	common /stainfo/stalat,stalon,stah
	open(unit=5,file=fname,status='old')
	read(5,*)sr1
      sr=1/sr1
      do 9 ii=1,5,1
      read(5,*)
9     continue
	read(5,*)x,stalat,stalon,stah
      do 10 ii=1,7,1
      read(5,*)
10    continue
      read(5,*)yr,doy,hour,min,sec
 	write(event,50)int(yr),int(doy)
50	format(i4.4,'/',i3.3,'    ')
      write(stime,51)int(hour),int(min),int(sec)
51	format(2(i2.2,':'),i2.2,'.000')
      read(5,'(44x,i6)')nsamp
      do 11 ii1=1,6,1
      read(5,*)
11    continue
      read(5,'(a4)')staname
      do 12 ii2=1,5,1
      read(5,*)
12    continue
      read(5,'(16x,a3)')comp
      read(5,*)
      read(5,*) (trace(i),i=1,nsamp)
	close(unit=5)
	return
	end




	subroutine repair_trace(fname)
	character fname*(*),line*100(100000)
	call runner('copy '//fname//' original.sac',0,1)
	open(unit=21,file=fname,status='unknown')
	do 10 i=1,100000,1
	read(21,'(a100)',end=11,err=11)line(i)
10	continue
11	nline=i-1
	do 20 i=31,nline,1
		do 25 j=15,75,15
26			if(line(i)(j:j+1).eq.'00')then
				do 30 k=j,99,1
				line(i)(k:k)=line(i)(k+1:k+1)
30				continue
			goto 26
			elseif(line(i)(j:j).eq.'0')then
			line(i)(j:j)=' '
			endif
25		continue
20	continue
	rewind(unit=21)
	do 50 i=1,nline,1
	write(21,'(a)')line(i)(1:len(line(i)))
50	continue
	close(unit=21)
	return
	end




	subroutine readpar()
	common /par/igmtn,iglit,itrend,ifilt, iresp,taper,norder,  
     *	flcut, fhcut, ismooth, fbox, freq1, freq2, nfreq, damp
	open(unit=6,file='parameter.dat',status='old')
	open(unit=61,file='parameter.out',status='unknown')
      read(6,*)iglit,itrend,ifilt,iresp
      read(6,*)taper
      read(6,*)norder, flcut, fhcut
      read(6,*) freq1, freq2, nfreq, damp
      read(6,*)ismooth, fbox
      write(61,*)iglit,itrend,ifilt,iresp
      write(61,*)taper
      write(61,*)norder, flcut, fhcut
      write(61,*) freq1, freq2, nfreq, damp
      write(61,*)ismooth, fbox
      if (taper .gt. .50)write(*,*)' ERROR. Use taper<.5'
	close(unit=6)
	close(unit=61)
      return
      end




	function rounder(x,n)
      y=x*10.**n
      k=nint(y)
      rounder=real(k/10.**n)
      return
      end




	subroutine skipline(iunit,iline)
	do 10 i=1,iline,1
	read(iunit,*)
10	continue
	return
	end




	subroutine serialn(num,i)
	open(unit=71,file='ren_eve.b',err=5,form='binary',status='old')
	goto 10
5	open(unit=71,file='c:\windows\ren_eve.b',err=7,form='binary',
     *	status='old')
	goto 10
7	write(*,*)'Error finding the file (ren_eve.b), default serial=1'
	num=1
	open(unit=71,file='ren_eve.b',form='binary',status='new')
	write(71)num
	rewind(71)
10	read(71)num
	write(*,*)num
	new_num=num+i
	rewind(71)
	write(71)new_num
	close(71)
20	return
	end