program eqsel3
c
c ...... 25 August 1992
c
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,nverts,selcrit
       integer*4 eqcent
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs,open
       character*1 comment(2442),qualdat,felt,ans
       character*12 infile
       character*3 user
c
c ...... define common blocks
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /eqchar/ comment,qualdat,felt
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
c
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
c
       common /userinit/ user
c
       open (20,status='scratch')
c
c ...... set default values for all variables
c
       call default
c
c ...... write title screen, get user name
c
       write (*,*) ' '
       write (*,*) 'University of California Seismographic Stations '
       write (*,*) 'Earthquake Catalog Search Program       v. 3.02 '
       write (*,*) ' '
       write (*,1000)
       read (*,'(a3)') user
       write (*,*) ' '
c
c ...... get next operation
c    
  1    write (*,1001)
       read (*,'(a1)') ans
       write (*,*) ' '
c
c ...... do operation
c
       if ((ans.eq.'i').or.(ans.eq.'I')) then
      		call input(infile)
		open (10,file=infile,status='old',err=2)
		write (*,*) ' '
       endif
       if ((ans.eq.'o').or.(ans.eq.'O')) then
		call output
       		write (*,*) ' '
       endif
       if ((ans.eq.'p').or.(ans.eq.'P')) then
       		call crit
		write (*,*) ' '
       endif
       if ((ans.eq.'s').or.(ans.eq.'S')) then
		inquire(file=infile,opened=open)
		if (open) then
			call search
			rewind 10
			rewind 20
		else
			write (*,*) char(7)
			write (*,*) 'No input files are open!'
			goto 1
		endif
		write (*,*) ' '
       endif
       if ((ans.eq.'e').or.(ans.eq.'E')) goto 999
       goto 1
 2     write (*,*) char(7)
       write (*,*) 'File not found!'
       goto 1
 1000  format('Please enter your initials --> ',$)
 1001  format('Next Operation: (Input,Params,Search,Output,Exit) --> ',
     + $)
 999   stop
       end
c
c ...... input subroutine, get catalog type and filename
c
       subroutine input(infile)
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs
       character*1 ans
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       character*12 infile
       taca = .false.
       noca = .false.
       gs = .false.
   2   write (*,1000)
       read (*,'(a1)',err=2) ans
       if ((ans.eq.'b').or.(ans.eq.'B')) then
       		noca = .true.
		goto 4
       endif
       if ((ans.eq.'t').or.(ans.eq.'T')) then
       		taca = .true.
		goto 4
       endif
       if ((ans.eq.'u').or.(ans.eq.'U')) then
       		gs = .true.
		goto 4
       endif
       if ((ans.eq.'o').or.(ans.eq.'O')) then
       		noca = .true.
		goto 4
       endif
       if ((ans.eq.'q').or.(ans.eq.'Q')) goto 5
       goto 2
   4   write (*,1001)
       read (*,'(a12)',err=4) infile
 1000  format('Catalog Format (Townley-Allen,BRK,USGS,Other,Quit) --> '
     +   ,$)
 1001  format('Filename --> ',$)
 5     return
       end
c
c ....... selection criteria
c
       subroutine crit
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,nverts,selcrit
       logical*4 radius,latlong,polygon,mag,time,depth,fexist
       character*12 cfile
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
c
       radius = .false.
       polygon = .false.
       latlong = .false.
       mag = .false.
       time = .false.
       depth = .false.
   1   write (*,1000)
       read (*,'(a1)',err=1) ans
       if ((ans.eq.'q').or.(ans.eq.'Q')) goto 24
       write (*,*) ' '
c
c ....... next action
c
       if ((ans.eq.'r').or.(ans.eq.'R')) then
       		radius = .true.
  2		write (*,1001)
		read (*,*,err=2) maxrad
		if (maxrad.le.0.0) then
			radius = .false.
			write (*,*) ' '
			goto 1
		endif
  3		write (*,1002)
   		read (*,*,err=3) cenlat
		if (abs(cenlat).gt.90.0) then
			write (*,*) char(7)
			goto 3
		endif
  4  		write (*,1003)
  		read (*,*,err=4) cenlong
		if (abs(cenlong).gt.180.0) then
			write (*,*) char(7)
			goto 4
		endif
		write (*,*) ' '
       endif
       if ((ans.eq.'p').or.(ans.eq.'P')) then
  		polygon = .true.
  		write (*,1004)
  5		read (*,*,err=5) nverts
		if (nverts.eq.0) then
			polygon = .false.
			write (*,*) ' '
			goto 1
		endif
  6		write (*,1005)
   		read (*,'(a12)',err=6) cfile
		inquire(file=cfile,exist=fexist)
		if (fexist) then
  			open (40,file=cfile)
		else
			write (*,*) char(7)
			write (*,*) ' File not found!'
			goto 6
		endif
  		do 10 i = 1,nverts
  			read (40,*,err=6) plat(i),plong(i)
  10		continue
		write (*,*) ' '
       endif
       if ((ans.eq.'l').or.(ans.eq.'L')) then
  		latlong = .true.
   7		write (*,1006)
  		read (*,*,err=7) nlat
		if (abs(nlat).gt.90.0) then
			latlong = .false.
			write (*,*) ' '
			goto 1
		endif
   8		write (*,1007)
   		read (*,*,err=8) slat
		if (abs(slat).gt.90.0) then
			write (*,*) char(7)
			goto 8
		endif
   9		write (*,1008)
  		read (*,*,err=9) wlong
		if (abs(wlong).gt.180.0) then
			write (*,*) char(7)
			goto 9
		endif
  11		write (*,1009)
  		read (*,*,err=11) elong
		if (abs(elong).gt.180.0) then
			write (*,*) char(7)
			goto 11
		endif
  		if (nlat.le.slat) goto 7
  		if (elong.gt.wlong) goto 7
		write (*,*) ' '
       endif
       if ((ans.eq.'m').or.(ans.eq.'M')) then
		mag = .true.
 12		write (*,1010)
     		read (*,*,err=12) minmag
		if (minmag.lt.0.0) then
			mag = .false.
			write (*,*) ' '
			goto 1
		endif
 13  		write (*,1011)
     		read (*,*,err=13) maxmag
 14   		write (*,1012)
    		read (*,*,err=14) selcrit
    		if (minmag.gt.maxmag) goto 12
		write (*,*) ' '
       endif
       if ((ans.eq.'t').or.(ans.eq.'T')) then
    		time = .true.
 15		write (*,1013)
     		read (*,*,err=15) begyr
		if (begyr.le.0) then
			time = .false.
			write (*,*) ' '
			goto 1
		endif
 16   		write (*,1014)
    		read (*,*,err=16) begmo
		if ((begmo.lt.1).or.(begmo.gt.12)) then
			write (*,*) char(7)
			goto 16
		endif
 17   		write (*,1015)
     		read (*,*,err=17) begdy
		if ((begdy.lt.1).or.(begdy.gt.31)) then
			write (*,*) char(7)
			goto 17
		endif
 18   		write (*,1016)
     		read (*,*,err=18) endyr
 19   		write (*,1017)
    		read (*,*) endmo
		if ((endmo.lt.1).or.(endmo.gt.12)) then
			write (*,*) char(7)
			goto 19
		endif
 20   		write (*,1018)
    		read (*,*,err=20) enddy
		if ((enddy.lt.1).or.(enddy.gt.31)) then
			write (*,*) char(7)
			goto 20
		endif
		if (endmo.eq.2) then
			if (enddy.gt.29) goto 15
		endif
		if (begmo.eq.2) then
			if (begdy.gt.29) goto 15
		endif
		if (endyr.lt.begyr) goto 15
		if (enddy.lt.begdy) then
	           if ((endmo.eq.begmo).and.(endyr.eq.begyr)) goto 15
		endif
		if ((endmo.lt.begmo).and.(endyr.eq.begyr)) goto 15
		write (*,*) ' '
       endif
       if ((ans.eq.'d').or.(ans.eq.'D')) then
		depth = .true.
 21		write (*,1019)
     		read (*,*,err=21) mindep
		if (mindep.lt.0.0) then
			depth = .false.
			write (*,*) ' '
			goto 1
		endif
 22   		write (*,1020)
 
    		read (*,*,err=22) maxdep
    		if (mindep.gt.maxdep) goto 21
		write (*,*) ' '
       endif
       if ((ans.eq.'s').or.(ans.eq.'S')) then
    		if (radius) then
    			write (*,*) 'Circular Area Search -->'
			write (*,*) 'Maximum radius:   ',maxrad,' km'
			write (*,*) 'Center Latitude:  ',cenlat,' degrees N'
			write (*,*) 'Center Longitude: ',cenlong,' degrees W'
		endif
		if (polygon) then
			write (*,*) 'Polygon Search -->'
   			write (*,*) 'Vertex #     Latitude  Longitude'
   			do 23 i = 1,nverts
   				write (*,1021) i,plat(i),plong(i)
  23			continue
   		endif
   		if (latlong) then
   			write (*,*) 'Latitude/Longitude Box Search -->'
   			write (*,*) 'North Latitude: ',nlat,' degrees N'
  			write (*,*) 'South Latitude: ',slat,' degrees N'
  			write (*,*) 'West Longitude: ',wlong,' degrees W'
  			write (*,*) 'East Longitude: ',elong,' degrees W'
  		endif
  		if (mag) then
  			write (*,*) 'Magnitude Search -->'
  			write (*,*) 'Minimum Magnitude: ',minmag
  			write (*,*) 'Maximum Magnitude: ',maxmag
  			if (selcrit.eq.0) then
			write (*,*) 'constant magnitude threshold'
			else
			write (*,*) 'variable magnitude threshold'
			endif
		endif
		if (time) then
			write (*,*) 'Date Range Search -->'
			write (*,1022) begmo,begdy,begyr
			write (*,1023) endmo,enddy,endyr
		endif
		if (depth) then
			write (*,*) 'Depth Range Search -->'
			write (*,*) 'Minimum Depth: ',mindep,' kilometers'
			write (*,*) 'Maximum Depth: ',maxdep,' kilometers'
		endif
		write (*,*) ' '
       endif
       goto 1
 1000  format('Search Parameters: (Rad,Poly,LatLong,Mag,Time,Depth',
     +         ',Summary,Quit) --> ',$)
c
c ... radius search formats
c
 1001  format('Maximum radius (km) [enter 0.0 to exit] --> ',$)
 1002  format('Center Latitude (decimal degrees N) --> ',$)
 1003  format('Center Longitude (decimal degrees W) --> ',$)
c 
c ... polygon search formats
c
 1004  format('Number of vertices (enter 0 to exit) --> ',$)
 1005  format('Enter polygon coordinate filename --> ',$)
c
c ... lat/long box search formats
c
 1006  format('North Latitude (decimal degrees N) [enter 99.0 to exit]
     + --> ',$)
 1007  format('South Latitude (decimal degrees N) --> ',$)
 1008  format('West Longitude (decimal degrees W) --> ',$)
 1009  format('East Longitude (decimal degrees W) --> ',$)
c
c ... magnitude search formats
c
 1010  format('Minimum magnitude [enter -9.9 to exit] --> ',$)
 1011  format('Maximum magnitude --> ',$)
 1012  format('Constant or variable ML threshold (0 or 1) --> ',$)
c
c ... date search formats
c
 1013  format('Start year (yyyy) [enter 0 to exit] --> ',$)
 1014  format('Start month (mm) --> ',$)
 1015  format('Start day (dd) --> ',$)
 1016  format('End year (yyyy) --> ',$)
 1017  format('End month (mm) --> ',$)
 1018  format('End day (dd) --> ',$)
c
c ... depth search formats
c
 1019  format('Minimum depth (km) [enter -9.9 to exit] --> ',$)
 1020  format('Maximum depth (km) --> ',$)
c
c ... summary section formats
c
 1021  format(3x,i2,10x,f6.3,4x,f7.3)
 1022  format('Start Date: ',3x,i2,'/',i2,'/',i4)
 1023  format('End Date:   ',3x,i2,'/',i2,'/',i4)
 24    return
       end
c
c
c ...... reading,searching,writing subroutine
c
c
       subroutine search
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts,eqcent
       integer*4 selcrit,eqcen1
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs
       logical*4 year84
       character*1 comment(2442),qualdat,felt
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
       common /eqchar/ comment,qualdat,felt
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       common /yearflag/ year84
       year84 = .false.
       eqcent = 1900
       nhypo = 0
       nselect = 0
 1     exit = 0
c
c ...... read event data
c
       if (taca) call eqrd83(exit)
       if (noca) then
		if ((begyr.le.1983).and.(.not.year84)) then
                        call eqrd83(exit)
			if (exit.eq.1) goto 1
		else
			call eqrd84(exit)
			if (exit.eq.1) goto 1
		endif
       endif
       if (gs) call gsread(exit)
       if (exit.eq.2) goto 9999

c
c ...... run sorting routines
c
       if (time) call datesort(exit)
       if (exit.eq.1) goto 1
       if (exit.eq.2) goto 9999
       if (depth) then
		if (eqdep.lt.mindep) goto 1
		if (eqdep.gt.maxdep) goto 1
       endif
       if (mag) then
		if (eqmag.lt.minmag) goto 1
		if (eqmag.gt.maxmag) goto 1
       endif
       if (latlong) then
		if (eqlat.gt.nlat) goto 1
		if (eqlat.lt.slat) goto 1
		if (eqlong.gt.wlong) goto 1
		if (eqlong.lt.elong) goto 1
       endif
       if (radius) then
		call boxrad(exit)
		if (exit.eq.1) goto 1
		call cirsort(exit)
		if (exit.eq.1) goto 1
       endif
       if (polygon) then
		call boxpoly(exit)
		if (exit.eq.1) goto 1
		call polysort(exit)
		if (exit.eq.1) goto 1
       endif
c
c ...... write intermediate file, increment number of selected eqs
c
       nselect = nselect + 1
       if (mod(nselect,500).eq.0) then
       		write (*,1111) nselect,nhypo,eqhr,eqmn,eqmo,eqdy,eqyr
       endif
       eqyr = eqyr - eqcent
       if (gs) then
		write (20,1001) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,eqlat
     +		,eqlong,eqmag,eqdep
       endif
       if ((noca).or.(taca)) then
		if (taca) then
			if (eqcent.ne.eqcen1) then
				eqcen1 = eqcent
				write (20,'(2x,i4)') eqcent
			endif
		endif
       		write (20,1006) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,eqlat
     +		,eqlong,qualdat,eqmag,numsta,felt,(comment(ii),ii=1,42)
     +		,nc
     		if (nc .gt. 1) then
         		do 100 i=2,nc
          			ij=80*(i-2)+43
          			ik=80*(i-1)+42
          			write (20,1007) (comment(ii),ii=ij,ik)
 100	 		continue
       		endif
       endif
       goto 1
 1001  format(5i2.2,f5.2,f6.3,f7.3,1x,f3.1,9x,f5.2)
 1006  format(5i2.2,f5.2,f6.3,f7.3,a1,f3.1,1x,i2,a1,42a1,i2)
 1007  format(80a1)
 1111  format ('Events selected/read',i6,'/',i6,4x,'O-Time: ',i2,'h'
     + ,i2,'m',2x,i2,'/',i2,'/',i4)
9999  if ((eqyr-eqcent).le.0) then
 		eqyr = eqyr+eqcent
       endif
       write (*,1111) nselect,nhypo,eqhr,eqmn,eqmo,eqdy,eqyr
       return
       end

c
c ......... read in earthquake data
c
       subroutine eqrd83(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat
       real*4 elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts,eqcent
       integer*4 selcrit
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs,year84
       character*1 comment(2442),qualdat,felt
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
       common /eqchar/ comment,qualdat,felt
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       common /yearflag/ year84
       data blank/' '/
       exit = 0
 1     read (10,1005,err=11,end=10) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec
     + ,eqlat,eqlong,qualdat,eqmag,numsta,felt,(comment(ii),ii=1,42)
     + ,nc
       if (eqmo.gt.12) then
    	  	eqcent=100*eqmo
      		goto 1
       endif
       eqyr = eqyr + eqcent
       if ((eqyr.eq.1983).and.(eqmo.eq.12).and.(eqdy.eq.31).and.
     +	(eqhr.eq.22).and.(eqmn.eq.39)) then
		year84 = .true.
       endif
       if (eqyr .gt. endyr) exit = 2
       if (nc .gt. 1) then
       	 do 2 i=2,nc
        	ij=80*(i-2)+43
         	ik=80*(i-1)+42
          	read (10,1007) (comment(ii),ii=ij,ik)
 2      continue
       endif
       nhypo = nhypo + 1
    9  return
   11  write (*,*) eqyr,eqmo,eqdy
       write (*,*) eqhr,eqmn
   10  exit=2
 1005  format (5i2,f4.1,1x,f5.2,1x,f6.2,a1,1x,f3.1,1x,i2,a1,42a1,i2)
 1006  format (1x,5i2,f4.1,1x,f5.2,1x,f6.3,a1)
c,1x,f3.1,1x,i2,a1,42a1,i2)
 1007  format(80a1)
       return
       end
c
c ...... read gs format data
c
       subroutine gsread(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts,eqcent
       integer*4 selcrit
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       read (10,1000,err=10,end=10) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec
     +		,eqlat,eqlong,eqmag,eqdep
       nhypo = nhypo + 1
       eqyr = eqyr + eqcent
       if (eqyr .gt. endyr) exit = 2
       return
 10    exit=2
 1000  format(5i2,f5.2,f6.3,f7.3,1x,f3.1,9x,f5.2)
       return
       end
c
c ........ read nc format post 1984
c
       subroutine eqrd84(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts,eqcent
       integer*4 selcrit
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs,year84
       character*1 comment(2442),qualdat,felt
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
       common /eqchar/ comment,qualdat,felt
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       common /yearflag/ year84
       data blank/' '/
       read (10,1006,err=10,end=10) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,
     + eqlat,eqlong,qualdat,eqmag,numsta,felt,(comment(ii),ii=1,42),nc
 1006  format (5i2,f5.2,f6.3,f7.3,a1,f3.1,1x,i2,a1,42a1,i2)
       eqyr = eqyr + eqcent
       nhypo = nhypo + 1
       if (eqyr .gt. endyr) exit = 2
       if (nc .gt. 1) then
         do 600 i=2,nc
          	ij=80*(i-2)+43
          	ik=80*(i-1)+42
          	read (10,1007) (comment(ii),ii=ij,ik)
 1007	  format(80a1)
 600	 continue
       endif
       return
 10    exit=2
       return
       end
c
c ...... sort by date
c
       subroutine datesort(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 eqmd,selcrit,eqcent
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,emd,bmd,nverts,exit
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
c
c
c
       emd = endmo*100+enddy
       bmd = begmo*100+begdy
       eqmd = eqmo*100+eqdy
       if (eqyr.lt.begyr) exit = 1
       if ((eqmd .lt. bmd) .and. (eqyr .eq. begyr)) exit = 1
       if ((eqmd .gt. emd) .and. (eqyr .ge. endyr)) exit = 2
       return
       end

c
c ...... make lat/long box surrounding radius search area
c
       subroutine boxrad(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag,tlatmin,tlatmax
       real*4 tlonmin,tlonmax,dlat
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts,selcrit
       integer*4 eqcent
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
c
c ...... find min/max latitude
c
       dlat=maxrad/111.195
       tlatmin=cenlat-dlat
       tlatmax=cenlat+dlat
c
c ...... find min/max longitude
c
       dlong=maxrad/(111.195*cos(atan(1.0)*cenlat/45.0))
       tlonmin=cenlong-dlong
       tlonmax=cenlong+dlong
c
c ...... apply box search to event data
c
       if (eqlat.gt.tlatmax) exit = 1
       if (eqlat.lt.tlatmin) exit = 1
       if (eqlong.gt.tlonmax) exit = 1
       if (eqlong.lt.tlonmin) exit = 1
       return
       end
c
c ...... apply circular sort routine
c
       subroutine cirsort(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag,at,arg,po2,pi,xeq,yeq
       real*4 eqgc,xcen,ycen,cengc,eqrad
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts
       integer*4 selcrit,eqcent
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       at = atan(1.)
       arg = at/45.0
       po2 = 2.0*at
       pi = 4.0*at
c
c ...... calculate event geocentric colatitude
c
       xeq = eqlong*arg*(-1.0)
       yeq = eqlat*arg
       eqgc = po2-(atan(0.993277*tan(yeq)))
c
c ...... calculate search center geocentric colatitude
c
       xcen = cenlong*arg*(-1.0)
       ycen = cenlat*arg
       cengc = po2-(atan(0.993277*tan(ycen)))
       call rad(eqgc,xeq,cengc,xcen,eqrad)
c
c ...... work with magnitude threshold
c
       if (eqrad .gt. maxrad) then
       		exit = 1
       elseif (mag) then
       		if (selcrit.eq.0) then
			exit = 0
		else
 			call mlthresh(eqrad,minmag,threshml)
 			if (eqmag.lt.threshml) exit = 1
		endif
       endif
       return
       end
c
c ...... compute distance from center to earthquake
c
       subroutine rad(ec,el,sc,sl,d)
       real*4 at,po2,pi,tpi,c1,c2
       real*4 sc,sl,ec,el,d
       real*8 darg
       at = atan(1.)
       po2 = 2.0*at
       pi = 4.0*at
       tpi = 8.0*at
       c2 = 45.0/at
       c1 = 111.195*c2
       darg = dcos(dble(sc))*dcos(dble(ec))+
     1        dsin(dble(sc))*dsin(dble(ec))*dcos(dble(sl-el))
       d = c1*sngl(dacos(darg))
       return
       end
c
c .......... compute ML threshold with increasing distance
c
       subroutine mlthresh(radius,ml0,mlt)
       real*4 radius,ml0,mlt,mao(34),dist(34)
       data mao/1.4,1.4,1.5,1.6,1.7,1.9,2.1,2.3,2.4,2.5,2.6,2.7,2.8,2.9,
     + 3.0,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.0,4.1,4.2,4.3,4.4,4.5,
     + 4.6,4.7,4.8,4.9/
       data dist/0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,70.,85.,
     +100.,120.,140.,160.,180.,200.,210.,240.,260.,280.,300.,320.,340.,
     +370.,390.,420.,460.,500.,550.,600./
       icount = 0
       do 10 i = 1,34
		icount = icount +1
 		if (radius .lt. dist(i)) goto 20
 10    continue
       goto 30
 20    icount = icount - 1
 30    mlt = ml0 + mao(icount) - 1.4
c       write(*,'(1h ,a,3f10.2)') 'mlthresh: rad, ml0 & mlt = ',
c     1                            radius,ml0,mlt
       return
       end
c
c ...... generate box surrounding polygon search area
c
       subroutine boxpoly(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag,minlat,maxlat,minlong
       real*4 maxlong
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts,selcrit
       integer*4 eqcent
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
c
c
c
       maxlat = plat(1)
       minlat = plat(1)
       maxlong = plong(1)
       minlong = plong(1)
       do 10 i = 2,nverts
		if (plat(i).gt.maxlat) maxlat = plat(i)
		if (plat(i).lt.minlat) minlat = plat(i)
		if (plong(i).gt.maxlong) maxlong = plong(i)
		if (plong(i).lt.minlong) minlong = plong(i)
 10    continue
c
c ...... apply box search to event data
c
       if (eqlat.gt.maxlat) exit = 1
       if (eqlong.gt.maxlong) exit = 1
       if (eqlat.lt.minlat) exit = 1
       if (eqlong.lt.minlong) exit = 1
       return
       end
c
c ........... polygon sort
c
       subroutine polysort(exit)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,exit,nverts
       integer*4 selcrit,eqcent
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       real*4 xi,xj,yi,yj,a,yq
       integer*4 nl,nr,exit,icount
       nl = 0
       do 5 i = 1,nverts
       		icount = i + 1
      		xi = plat(i)
      		xj = plat(icount)
      		yi = plong(i)
      		yj = plong(icount)
      	   	if (((eqlong.ge.yi).and.(eqlong.ge.yj)).or.((eqlong.lt.yi)
     +	   	.and.(eqlong.lt.yj))) goto 5
      		if (yi.eq.yj) goto 5
      		if (xi.ne.xj) goto 1
      		if (eqlat.ge.xi) nl = nl +1
      		if (eqlat.lt.xi) nr = nr +1
      		goto 5
  1   		a = (yj - yi) / (xj - xi)
      		yq = eqlong - a * (eqlat - xi)
      		if (xj.lt.xi) goto 3
      		if (yj.lt.yi) goto 2
c
c ........ (x .gt. 0, y .gt. 0) quadrant test
c
      		if (yq .le. yi) nl = nl + 1
      		goto 5
c
c ........ (x .gt. 0, y .lt. 0) quadrant test
c
 2		if (yq .ge. yi) nl = nl + 1
      		goto 5
 3    		if (yj .lt. yi) goto 4
c
c ........ (x .lt. 0, y .gt. 0) quadrant text
c
      		if (yq.ge.yi) nl = nl + 1
      		goto 5
c
c ........ (x .lt. 0, y .lt. 0) quadrant test
c
 4    		if (yq .le. yi) nl = nl + 1
 5     continue
       nl = nl - 2*(nl/2)
       if (nl .eq. 1) exit = 1
       return
       end
c
c ....... output subroutines
c
       subroutine output
       character*1 ans
c
c ...... get next action
c
   1   write (*,1000)
       read (*,'(a1)',err=1) ans
       write (*,*) ' '
c
c ...... do actions
c
       if ((ans.eq.'F').or.(ans.eq.'f')) then
		call fileout
       endif
       if ((ans.eq.'L').or.(ans.eq.'l')) then
		call listout
       endif
       if ((ans.eq.'Q').or.(ans.eq.'q')) goto 9999
       goto 1
 9999  return
 1000  format('Output:  (File,List,Quit) --> ',$)
       end
c
c ...... subroutine for file output
c
       subroutine fileout
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 eqcent,mnc,nnc
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs
       character*1 comment(2442),qualdat,felt,ans,ans1
       character*12 outfile
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
       common /eqchar/ comment,qualdat,felt
c
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
c
c ...... select output file type, name
c
  1    write (*,1000)
       read (*,'(a12)') outfile
       open (30,file=outfile)
  2    write (*,1001)
       read (*,'(a1)') ans
       if ((ans.eq.'B').or.(ans.eq.'b')) then
       		write (*,1002)
       		read (*,'(a1)') ans1
       endif
       write (*,*) ' '
c
c ...... do appropriate action
c
       if ((ans.eq.'Q').or.(ans.eq.'q')) goto 9999
       if ((ans.eq.'B').or.(ans.eq.'b')) then
		if (gs) then
			write (*,*) char(7)
			write (*,*) 'You selected the BRK format with a '//
     +			'USGS-format data file.'
			write (*,*) 'Please select the USGS format.'
			write (*,*) ' '
			goto 2
		endif
		rewind 20
  3		read (20,1003,end=9999,err=999) eqyr,eqmo,eqdy,eqhr
     +          ,eqmn,eqsec,eqlat,eqlong,qualdat,eqmag,numsta,felt
     +          ,(comment(ii),ii=1,42),nc
     		if (nc .gt. 1) then
        		do 100 i=2,nc
          			ij=80*(i-2)+43
          			ik=80*(i-1)+42
          			read (20,1004) (comment(ii),ii=ij,ik)
 100	 		continue
		endif
		if (taca) then
			if (eqmo.gt.12) then
				write (30,*) eqmo*100
 				goto 3
			endif
		endif
		if ((ans1.eq.'F').or.(ans1.eq.'f')) then
       			write (30,1003) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,eqlat
     +			,eqlong,qualdat,eqmag,numsta,felt,(comment(ii),ii=1,42)
     +			,nc
     			if (nc .gt. 1) then
        			do 101 i=2,nc
          				ij=80*(i-2)+43
          				ik=80*(i-1)+42
          				write (30,1004) (comment(ii),ii=ij,ik)
 101	 			continue
       			endif
		elseif ((ans1.eq.'M').or.(ans1.eq.'m')) then
			mnc = 1
			write (30,1003) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,eqlat
     +			,eqlong,qualdat,eqmag,numsta,felt,(comment(ii),ii=1,42)
     +			,mnc
		else
			nnc = 1
			write (30,1100) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,eqlat 
     +			,eqlong,qualdat,eqmag,numsta,felt,nnc
		endif
		goto 3
       endif
       if ((ans.eq.'U').or.(ans.eq.'u')) then
		rewind 20
  4		if (gs) then
  			read (20,1005,err=999,end=9999) eqyr,eqmo,eqdy,eqhr
     +			,eqmn,eqsec,eqlat,eqlong,eqmag,eqdep
		endif
		if ((noca).or.(taca)) then
  5			read (20,1003,end=9999,err=999) eqyr,eqmo,eqdy,eqhr
     +          	,eqmn,eqsec,eqlat,eqlong,qualdat,eqmag,numsta,felt
     +          	,(comment(ii),ii=1,42),nc
     			if (nc .gt. 1) then
        		do 110 i=2,nc
          			ij=80*(i-2)+43
          			ik=80*(i-1)+42
          			read (20,1004) (comment(ii),ii=ij,ik)
 110	 		continue
			endif
			if (taca) then
			if (eqmo.gt.12) then
				write (30,*) eqmo*100
 				goto 5
			endif
			endif
		endif
		write (30,1005) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,eqlat
     +		,eqlong,eqmag,eqdep
		goto 4
       endif
 999   write (*,*) eqyr,eqmo,eqdy
       write (*,*) eqhr,eqmn,eqsec
 9999  return
 1000  format('Enter output filename --> ',$)
 1001  format('Enter output file format: (BRK,USGS,Quit) --> ',$)
 1002  format('Type of comments: (None,Full,Map) --> ',$)
 1003  format(5i2.2,f5.2,f6.3,f7.3,a1,f3.1,1x,i2,a1,42a1,i2)
 1004  format(80a1)
 1005  format(5i2.2,f5.2,f6.3,f7.3,1x,f3.1,9x,f5.2)
 1100  format(5i2.2,f5.2,f6.3,f7.3,a1,f3.1,1x,i2,a1,42x,i2)
       end
c
c ...... listing output routine
c
       subroutine listout
       implicit integer*4 (i)
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,nverts,selcrit
       integer*4 eqcent
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs
       character*1 comment(2442),qualdat,felt,ans1,blank,ans2
       character*3 user,amo(12)
       character*30 pfile
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
       common /eqchar/ comment,qualdat,felt
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,
     + endyr,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       common /userinit/ user
c
       data amo/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug',
     + 		'Sep','Oct','Nov','Dec'/
       data blank/' '/
       
c
c ...... get listing type,printer type
c
 1     write (*,1000)
       read (*,'(a1)') ans1
       if ((ans1.eq.'c').or.(ans1.eq.'C')) then
       	if (gs) then
       		write (*,*) char(7)
		write (*,*) 'You selected the "Comments" format with a '//
     +		'USGS-format data file.'
		write (*,*) 'Please select the "NoComments" format.'
		write (*,*) ' '
		goto 1
	endif
	ans1 = 'c'
       endif
       if ((ans1.eq.'n').or.(ans1.eq.'N')) ans1 = 'n'
       write (*,1001)
       read (*,'(a30)') pfile
       open (11,file=pfile)
       write (*,1100)
       read (*,'(a1)') ans2
       write (*,*) ' '
c
c ...... write printer file to be printed later
c
       if (noca) call nchead
       if (taca) call tahead
       write (11,'(a2)') char(13)//char(12)
       if (ans2.eq.'y') call headpage
       write (11,'(a2)') char(13)//char(12)
       rewind 20
  2    if (gs) then
      	read (20,1002,err=999,end=9999) eqyr,eqmo,eqdy,eqhr,eqmn,eqsec,
     + 	eqlat,eqlong,eqmag,eqdep
       endif
       if ((noca).or.(taca)) then
  		read (20,1003,end=9999,err=999) eqyr,eqmo,eqdy,eqhr
     +          ,eqmn,eqsec,eqlat,eqlong,qualdat,eqmag,numsta,felt
     +          ,(comment(ii),ii=1,42),nc
     		if (nc .gt. 1) then
        		do 100 i=2,nc
          			ij=80*(i-2)+43
          			ik=80*(i-1)+42
          			read (20,1004) (comment(ii),ii=ij,ik)
 100	 		continue
		endif
		if ((taca).and.(eqmo.gt.12)) goto 1
       endif
       imot=eqmo/10
       imoo=eqmo-10*imot
       idyt=eqdy/10
       idyo=eqdy-10*idyt
       ihrt=eqhr/10
       ihro=eqhr-10*ihrt
       imnt=eqmn/10
       imno=eqmn-10*imnt
       itsec=10.*(eqsec+0.05)
       isech=itsec/100
       isect=(itsec-100*isech)/10
       iseco=itsec-100*isech-10*isect
       if (ans1.eq.'n') then
      		write (11,1005) idyt,idyo,amo(eqmo),eqyr
     +		,ihrt,ihro,imnt,imno,isech,isect,iseco,eqlat,eqlong,eqmag
     +		,eqdep
       endif
       if (ans1.eq.'c') then
 	if(numsta.eq.0) then
 	   if(eqmag.ne.0.) then
                write (11,1006) idyt,idyo,amo(eqmo),eqyr
     +		,ihrt,ihro,imnt,imno,isech,isect,iseco,eqlat,eqlong,eqmag
     +		,qualdat,felt,(comment(ii),ii=1,42)
             else
                write (11,1007) idyt,idyo,amo(eqmo),eqyr
     +		,ihrt,ihro,imnt,imno,isech,isect,iseco,eqlat,eqlong
     +		,qualdat,felt,(comment(ii),ii=1,42)
           endif
          else
           if(eqmag.ne.0) then
        	write (11,1008) idyt,idyo,amo(eqmo),eqyr
     +		,ihrt,ihro,imnt,imno,isech,isect,iseco,eqlat,eqlong,eqmag
     +		,numsta,felt,(comment(ii),ii=1,42)
             else
        	write (11,1009) idyt,idyo,amo(eqmo),eqyr
     +		,ihrt,ihro,imnt,imno,isech,isect,iseco,eqlat,eqlong,numsta
     +		,felt,(comment(ii),ii=1,42)
            endif
         endif
            if(nc.eq.1) goto 300
               do 200 i=2,nc
               ij=43+80*(i-2)
               is=42+80*(i-1)
               write (11,1010) (comment(ii),ii=ij,is)
 200           continue
 300        continue
       endif
       goto 2
 1000  format('Listing Types: (Comments,NoComments,Quit) --> ',$)
 1001  format('Listing filename --> ',$)
 1002  format(5i2.2,f5.2,f6.3,f7.3,1x,f3.1,9x,f5.2)
 1003  format(5i2.2,f4.1,1x,f6.3,f7.3,a1,f3.1,1x,i2,a1,42a1,i2)
 1004  format(80a1)
 1005  format(1h ,2i1,a3,i2,1x,6i1,'.',i1,1x,f6.3,1x,f7.3,1x,f3.1,f5.2)
 1006  format(1h ,2i1,a3,i2,1x,6i1,'.',i1,1x,f6.3,1x,f7.3,1x,f3.1,3x,a1,
     + 1x,a1,1x,42a1)
 1007  format(1h ,2i1,a3,i2,1x,6i1,'.',i1,1x,f6.3,1x,f7.3,8x,a1,1x,a1,1x,
     + 42a1)
 1008  format(1h ,2i1,a3,i2,1x,6i1,'.',i1,1x,f6.3,1x,f7.3,1x,f3.1,1x,i2,
     + 1x,a1,1x,42a1)
 1009  format(1h ,2i1,a3,i2,1x,6i1,'.',i1,1x,f6.3,1x,f7.3,7x,i2,1x,a1,1x,
     + 42a1)
 1010  format(1h ,80a1)
 1100  format('Print header page (y/n) ',$)
 999   write (*,*) eqyr,eqmo,eqdy
       write (*,*) eqhr,eqmn,eqsec
       write (*,*) eqlat,eqlong,eqmag,eqdep
 9999  return
       end
 
c
c ........ northern california header
c
      subroutine nchead
      character*61 head(61)
      data (head(i),i=1,10)/
     1 '   EARTHQUAKES IN NORTHERN CALIFORNIA AND ADJACENT AREAS     ',
     2 '    SELECTED FROM THE MASTER LIST (Jan 1910 to Current)      ',
     3 '               OF THE SEISMOGRAPHIC STATIONS                 ',
     4 '            UNIVERSITY OF CALIFORNIA AT BERKELEY             ',
     5 '                                                             ',
     6 '                                                             ',
     7 '                                                             ',
     8 '                        - SOURCES -                          ',
     9 '                                                             ',
     1 '                                                             '/
      data (head(i),i=11,20)/
     1 '"Catalogue of Earthquakes in Northern California and Adjoin-',
     2 '     ing Areas: 1 January 1910 - 31 December 1972", B. A.   ',
     3 '     Bolt and R. D. Miller, 567pp., Seismographic Stations, ',
     4 '     University of California, Berkeley, California, 1975.  ',
     5 '                                                            ',
     6 '"Bulletin of the Seismographic Stations", January 1, 1973 to',
     7 '     December 31, 1991, University of California, Berkeley, ',
     8 '     California.                                            ',
     9 '                                                            ',
     1 'Data on earthquakes occurring after December 31, 1991 are   '/
      data (head(i),i=21,30)/
     1 '     from the master list at the Seismographic Station which',
     2 '     is subject to revision pending publication of the Bull-',
     3 '     etin of the Seismographic Stations.                    ',
     4 '                                                            ',
     5 '                                                            ',
     6 '                         - NOTES -                          ',
     7 '                                                            ',
     8 '                                                            ',
     9 'On December 20, 1932 the Intensity Scale changed from the   ',
     1 '     Rossi-Forel to the Modified Mercalli.  Observed inten- '/
      data (head(i),i=31,40)/
     1 '     sities are given by Roman numerals in the comments.    ',
     2 '                                                            ',
     3 'HEADER -                                                    ',
     4 '     DATE - day, month & year                               ',
     5 '     TIME - origin time in UTC (subtract 8 hours for PST)   ',
     6 '     LAT  - latitude in decimal degrees north               ',
     7 '     LONG - longitude in decimal degrees west               ',
     8 '     ML   - Richter Magnitude (local magnitude scale)       ',
     9 '     QU   - quality of solution / number of stations used   ',
     1 '     F    - F in column indicates that earthquake was felt  '/
      data head(41)/
     1 '                                                            '/
      data (head(i),i=42,51)/
     1 '************************************************************',
     2 '*  The listed earthquake data were assembled from the      *',
     3 '*  above sources and compiled for in-house use. Although   *',
     4 '*  considerable care was taken in the transfer procedures, *',
     5 '*  some errors were introduced from this source. There-    *',
     6 '*  fore, for precise and detailed studies using these      *',
     7 '*  lists, the original sources should be consulted. Some   *',
     8 '*  amendments and corrections have been successively       *',
     9 '*  incorporated in the lists. Any remaining inconsist-     *',
     1 '*  tencies or suspicious values should be drawn to the     *'/
      data (head(i),i=52,61)/
     1 '*  attention of the director.  Checks can then be carried  *',
     2 '*  out against the original seismograms kept in store at   *',
     3 '*  Berkeley.                                               *',
     4 '*                                                          *',
     5 '*                     - CONTACT -                          *',
     6 '*            Director of the Seismographic Stations        *',
     7 '*            475 Earth Sciences Building                   *',
     8 '*            University of California                      *',
     9 '*            Berkeley, California  94720                   *',
     1 '************************************************************'/
      write (11,*) ' '
      write (11,*) ' '
      write (11,*) ' '
      write (11,*) ' '
      write(11,2) (head(i),i=1,61)
    2 format(1h ,7x,a61)
      return
      end
c
c........Townley-Allen Catalog Header
c
      subroutine tahead
      write (11,*) ' '
      write (11,1000)
1000  format (1h ,/,46x,'EARTHQUAKES IN CALIFORNIA (1769-1927)'//38x,
     1 'SEISMICITY  LIST  PREPARED FOR THE NATIONAL OCEANIC AND',/,38x,
     2 'ATMOSPHERIC ADMINISTRATION BY DONALD A. SEEBURGER UNDER',/,38x,
     3 'THE SUPERVISION OF BRUCE A. BOLT, DIRECTOR OF THE SEIS-',/,38x,
     4 'MOGRAPHIC STATIONS, UNIVERSITY OF CALIFORNIA, BERKELEY.',/,38x,
     5 'PROGRAMMING AID BY ROBERT UHRHAMMER. STUDY COMPLETED IN',/,38x,
     6 'JUNE, 1976.',//)
      write (11,1001)
1001  format (1h ,59x,'- SOURCES -',//,17x,'B. A. Bolt and R. D. Miller,
     1 "Catalogue of Earthquakes in Northern California and Adjoining Ar
     2eas,',/,17x,'1 January 1910 - 31 December 1972", Seismographic Sta
     3tions, University of California at Berkeley, 1975.',//,
     4 17x,'"Earthquake History of the United States",  Revised  1970,
     5U. S. Department  of  Commerce,  U. S.',/,17x,'Printing Office, 19
     673.',//,
     7 17x,'B. Gutenberg  and  C. F. Richter,  "Seismicity of the Earth"
     8,  1954,  Princeton University  Press.',//,
     9 17x,'E. S. Holden,  "Catalogue of Earthquakes on the Pacific Coas
     At,  1769 - 1897",  Reprinted from  the',/,17x,'Smithsonian Miscell
     Baneous Collections, No. 1087.',//,
     C 17x,'S. D. Townley  and  M. W. Allen,  "Descriptive Catalog of Ea
     Drthquakes of the Pacific Coast of  the',/,17x,'United States, 1869
     E to 1928", B.S.S.A., V. 29, No. 1, Jan 1939.',//,
     F 17x,'H. O. Wood and N. H. Heck, "Earthquake History of the United
     G States, Part 22",  U.S.C. and G.S.',/,17x,'Serial No. 609, Rev
     Hised 1951.')
      write (11,*) ' '
      write (11,1002)
1002  format (1h ,//,58x,'- INTRODUCTION -',//,22x,'This  list  is  esse
     1ntially  an  edited  abstract of the Townley-Allen Catalog of  Cal
     2ifornia',/,17x,'Earthquakes,  1769-1927.  The Townley-Allen  Catal
     3og has been revised where necessary after making',/,17x,
     4 'seismological  inferences  based on  the modern understanding of
     5 eartjquakes,  tectonic structure,',/,17x,
     6 'published  field observations  and  damage reports available,  a
     7nd after making assessments of the',/,17x,'historical accuracy  of
     8  the catalog entries by means of  intercomparison and internal co
     9nsistency',/,17x,'checks.  Account was also taken of some recent s
     Apecial studies of individual earthquakes.',//,
     B 22x,'The list entries give the Date,  Time (GMT),  Location, and
     Ca summary of the more significant',/,17x,'observations associated
     D with  the  earthquake.  Included in the description is an estimat
     Ee of the',/,17x,
     F 'maximum  Modified Mercalli Intensity experienced during the eart
     Ghquake.  The Townley-Allen Catalog',/,17x,'uses the Rossi-Forel In
     Htensity Scale.  A listing of the Modified Mercalli Scale and an ap
     2proximate',/,17x,'conversion scheme from  Rossi-Forel  to Modified
     J Mercalli are to be found in Richter,  "Elementary')
      write (11,1003)
1003  format (1h ,16x,'Seismology",  pp. 135-39,  650-52.  Modified Merc
     1alli  Intensities  were assigned  based upon  the'/,17x,'following
     2 criteria:',//,22x,'1)  If a description of earthquake effects exi
     3sts, a Modified Mercalli Intensity was assigned',/,22x,
     4 'after considering both the given Rossi-Forel Intensity and the d
     5escription.',//,22x,
     6 '2)  If no description  was  available,  the conversion was made
     7 based upon relations such as',/,22x,'those presented  in  Richter
     8,  referred to above.  If the conversion relations give an inter-'
     9,/,22x,'mediate value, then, lacking further evidence, the given M
     Aodified Mercalli Intensity  is  the'/,22x,
     B 'lower of the two values.',//,17x,'All magnitudes listed are take
     Cn from Gutenberg and Richter, "Seismicity of the Earth".',//)
      write (11,1004)
1004  format (1h ,61x,'- NOTES -',//,17x,'1)  Time listed is Greenwich M
     1ean Time.',//,17x,'2)  Earthquake locations are placed at the infe
     2rred center of the zone of highest intensity.',//,17x,
     3 '3)  Entry Reference:  F - The  Townley-Allen entry  is geographi
     4cal, the list entry gives the Time',/,43x,'and Location',//,39x,
     5 'G - The Townley-Allen entry gives a description of the  earthqua
     6ke  effects,',/,43x,'the list entry is an abstract of the Townley-
     7Allen account.',//,39x,
     8 'H - The list entry is a  major revision  of the  Townley-Allen a
     9ccount,  the',/,43x,'source of the revision is noted.'//)
      return
      end
c
c ...... print header page for search results
c
       subroutine headpage
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,nverts,selcrit
       integer*4 eqcent
       logical*4 radius,latlong,polygon,mag,time,depth
       logical*4 noca,taca,gs
       character*3 user
c
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
c
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,
     + endyr,endmo,enddy,nverts,selcrit
       common /flag/ radius,latlong,polygon,mag,time,depth,noca
     + ,taca,gs
       common /userinit/ user
c
	write (11,*) ' '
	write (11,*) ' '
	write (11,*) ' '
c
	if (gs) write (11,1000)
	if (noca) write (11,1001)
	if (taca) write (11,1002)
	write (11,*) ' '
	write (11,*) 'Search Parameters: '
	write (11,*) ' '
	if (radius) then
		write (11,*) 'Circular Area Search '
		write (11,*) ' '
		write (11,1003) maxrad
		write (11,1004) cenlat
		write (11,1005) cenlong
		write (11,*) ' '
	endif
	if (polygon) then
		write (11,*) 'Polygonal Area Search '
		write (11,*) ' '
		if (nverts.le.20) then
			write (11,*) ' Vert #  Latitude  Longitude '
			do 10 i = 1,20
				write (11,1006) i,plat(i),plong(i)
  10			continue
		elseif (nverts.le.40) then
		write (11,1100)  			
			do 20 i = 1,20
			write (11,1007) i,plat(i),plong(i),i+20,
     +			plat(i+20),plong(i+20)
  20			continue
		else
			write (11,1101) 			
			do 30 i = 1,20
			write (11,1008) i,plat(i),plong(i),i+20,
     +			plat(i+20),plong(i+20),i+40,plat(i+40),
     +			plong(i+40)
  30			continue
		endif
		write (11,*) ' '
	endif
	if (latlong) then
		write (11,*) 'Latitude/Longitude Box Area Search '
		write (11,*) ' '
		write (11,1009) nlat
		write (11,1010) slat
		write (11,1011) wlong
		write (11,1012) elong
		write (11,*) ' '
	endif
	if (mag) then	
		write (11,*) 'Magnitude Range Search '
		write (11,*) ' '
		write (11,1013) minmag
		write (11,1014) maxmag
		if (selcrit.eq.0) then
			write (11,*) 'constant magnitude threshold'
		else
			write (11,*) 'variable magnitude threshold'
		endif
		write (11,*) ' '
	endif
	if (time) then
		write (11,*) 'Date Range Search '
		write (11,*) ' '
		write (11,1015) begmo,begdy,begyr
		write (11,1016) endmo,enddy,endyr
		write (11,*) ' '
	endif
	if (depth) then
		write (11,*) 'Depth Range Search '
		write (11,*) ' '
		write (11,1017) mindep
		write (11,1018) maxdep
	endif
	write (11,*) ' '
	write (11,*) 'Events Selected: ',nselect
	write (11,*) ' '
	write (11,1020) user
	write (11,*) ' '
c
c ...... formats
c
 1000	format('Catalog Searched:     usgs.cat')
 1001   format('Catalog Searched:     brk1090.cat')
 1002   format('Catalog Searched:     townalln.cat')
 1003   format('Maximum Radius:    ',f6.2,' kilometers')
 1004   format('Center Latitude:   ',f6.3,' degrees N')
 1005   format('Center Longitude: ',f7.3,' degrees W')
 1006   format(3x,i2,6x,f6.3,4x,f7.3)
 1007   format(3x,i2,6x,f6.3,4x,f7.3,10x,i2,6x,f6.3,4x,f7.3)
 1008   format(3x,i2,6x,f6.3,4x,f7.3,10x,i2,6x,f6.3,4x,f7.3,
     +  10x,i2,6x,f6.3,4x,f7.3)
 1009   format('North Latitude:  ',f6.3,' degrees N')
 1010   format('South Latitude:  ',f6.3,' degrees N')
 1011   format('West Longitude: ',f7.3,' degrees W')
 1012   format('East Longitude: ',f7.3,' degrees W')
 1013   format('Minimum Magnitude: ',f3.1)
 1014   format('Maximum Magnitude: ',f3.1)
 1015   format('Start Date: ',3x,i2,'/',i2,'/',i4)
 1016   format('End Date:   ',3x,i2,'/',i2,'/',i4)
 1017   format('Minimum Depth:  ',f4.1,' kilometers')
 1018   format('Maximum Depth:  ',f4.1,' kilometers') 
 1020   format('UC Berkeley Seismographic Stations',
     +  '     Printed:  ',a3)
 1100	format(' Vert #  Latitude  Longitude ','     Vert #',
     +  '  Latitude  Longitude ')
 1101   format(' Vert #  Latitude  Longitude ','     Vert #',
     +  '  Latitude  Longitude ','     Vert #  Latitude  Longitude ')
	return
	end


c ...... set default values
c
       subroutine default
       real*4 eqlat,eqlong,eqmag,eqdep,eqsec
       real*4 cenlat,cenlong,maxrad,plat(51),plong(51),nlat,slat,elong
       real*4 wlong,mindep,maxdep,minmag,maxmag
       integer*4 eqyr,eqmo,eqdy,eqhr,eqmn,numsta,nc,nch,nhypo,nselect
       integer*4 begyr,begmo,begdy,endyr,endmo,enddy,nverts,selcrit
       integer*4 eqcent
       common /eqdata/ eqlat,eqlong,eqmag,eqdep,eqyr,eqmo,eqdy,eqhr,eqmn
     + ,eqsec,nc,nch,numsta,nhypo,nselect,eqcent
       common /srchparm/ cenlat,cenlong,maxrad,plat,plong,nlat,slat
     + ,elong,wlong,mindep,maxdep,minmag,maxmag,begyr,begmo,begdy,endyr
     + ,endmo,enddy,nverts,selcrit
       eqlat = 0.0
       eqlong = 0.0
       eqmag = 0.0
       eqdep = 0.0
       eqsec = 0.0
       cenlat = 0.0
       cenlong = 0.0
       maxrad = 0.0
       nlat = 0.0
       slat = 0.0
       elong = 0.0
       wlong = 0.0
       mindep = 0.0
       maxdep = 0.0
       minmag = 0.0
       maxmag = 0.0
       selcrit = 0
       nverts = 0
       do 10 i = 1,51
	plat(i) = 0.0
	plong(i) = 0.0
  10   continue
       eqyr = 0
       eqmo = 1
       eqdy = 1
       eqhr = 0
       eqmn = 0
       numsta = 0
       nc = 0
       nch = 0
       begyr = 1
       begmo = 1
       begdy = 1
       endyr = 2000
       endmo = 12
       enddy = 31
       return
       end