      subroutine st_pfd(t1p_time)
c
c ...... decode seistool pickfile format
c
c
c ...... 06 March 1995 - modified glitch detection algorithm
c
c Modified 21 nov 94 to acquire file names from config.h
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
c
      integer*4 ncr,nsp,ns1,ns2,ns3,ns4
      real*8 t1p_time
      common/files/ ncr,nsp,ns1,ns2,ns3,ns4
      character*80 fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1             epi_dt
      character*9 epi_id
      common /fnamr/ fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1               epi_id,epi_dt
c
      ncr=45
      open(ncr,file=fnampf(1:lnblnk(fnampf)),status='old')
c
c ...... read phase data
c
      call obtime(t1p_time)
      close(ncr)
c
      return
      end

      subroutine obtime(t1p_time)
c
c ...... input observed phases and onset times
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      include "config.h"
      character*4 stn
      character*3 strm,stream(1000)
      character*72 id
      character*1 iwtp,iwts,iwt,polar
      real*8 pt,st,syf,sxf,psaf,ssaf,t1p_time
      integer*4 nsreq
      common/files/ ncr,nsp,ns1,ns2,ns3,ns4
      common/obs/pt(1000),st(1000),nsreq
      common/obs1/stn(1000),polar(1000)
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns
      common/time/eqhour,isyear,isdoy,ishour,ismin,fssec
      common/wght/iwtp(1000),iwts(1000),iwt(1000)
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      real*8 etsmp(1000)
      character*80 aline1,aline2
      character*4 astn,aistn
      character*1 apoii,aptii,apfii
      character*1 apo(1000),aptype(1000),apfm(1000)
      real*8 fpotime
      data con /0.993277d0/
      do 1 i=1,1000
      etsmp(i)=0.d0
      apo(i)=' '
      iwtp(i)=' '
      iwts(i)=' '
      aptype(i)=' '
      apfm(i)=' '
      pt(i)=0.d0
      st(i)=0.d0
    1 continue
      isp=0
      nmstn=0
      ietc=0
c
      inpot=0
      insot=0
   11 continue
      aline1(1:1)=' '
      read(ncr,'(a)') aline1
      if(aline1(1:1).eq.'P') go to 11
      ilnb=lnblnk(aline1)
      if(ilnb.gt.72) ilnb=72
      id=aline1(1:ilnb)
c      write(*,'(2a)') 'id = ',id(1:ilnb)
      call l1dec(aline1,isyear,isdoy,ishour,ismin,fssec)
      ii = 0
    2 ii = ii + 1
      read(ncr,'(a)',end=5) aline2
c      write(*,'(a,i5,2a)') 'aline2 = ',ii,'   ',aline2
      ier=0
      call l2dec(aline2,astn,strm,apoii,aptii,apfii,fpotime,ier)
c      write(*,'(a,1x,a4,1x,a3,1x,a1,1x,a1,1x,a1,1x,f10.4)')
c     1      'l2dec = ',astn,strm,apoii,aptii,apfii,fpotime
      if(ier.gt.0) then
        ii=ii-1
        go to 2
      endif
      apo(ii)=apoii
      aptype(ii)=aptii
      if((aptype(ii).ne.'P').and.(aptype(ii).ne.'S')) then
        ii=ii-1
        go to 2
      endif
c
      if(ii.eq.1) then
        tbors=3600.d0*dble(ishour)+60.d0*dble(ismin)+fssec
        fpotime=tbors+fpotime
c        call ts_calc(fpotime,t_shift)
        t1p_time=fpotime
      else
        fpotime=tbors+fpotime
      endif
c
      apfm(ii)=apfii
      open(38,file=instr_coord(1:lnblnk(instr_coord)),status='old')
    3 continue
      read(38,'(a4,1x,2f9.4)',end=4) aistn,flat,flong
      if(astn.eq.aistn) then
        stream(ii)=strm
        stn(ii)=astn
        sy=flat
        syf(ii)=datand(con*dtand(sy))
        sxf(ii)=flong
c        write(*,'(a,i5,3x,a4,3x,2f15.7)') 
c     1        'station ',ii,stn(ii),syf(ii),sxf(ii)
        psaf(ii)=0.
        ssaf(ii)=0.
        if(aptype(ii).eq.'P') then
          pt(ii)=fpotime
c          write(*,'(a,i2,a,f15.4)') 'pt(',ii,') = ',pt(ii)
          iwtp(ii)='1'
        elseif(aptype(ii).eq.'S') then
          st(ii)=fpotime
c          write(*,'(a,i2,a,f15.4)') 'st(',ii,') = ',st(ii)
          iwts(ii)='1'
        endif
        polar(ii)=apfm(ii)
        close(38)
c
c ...... check whether station is in Northern and Central California
c	 in region bounded by: (41,-124.5),(41.,-121.5),(38.,-119.5),
c	 (35.5,-119.5),(35.5,-124.5)
c
c	write(*,*) 'stn, flat, flong = ',astn," ",flat,flong
c
	if(flat.gt.41.d0) then
	  ii=ii-1
	elseif(flat.lt.35.5d0) then
	  ii=ii-1
	elseif(flong.gt.-119.5d0) then
	  ii=ii-1
	elseif(flat.ge.38.d0) then
	  if(flong.gt.(-94.1667d0-2.d0*flat/3.d0)) ii=ii-1
	endif
c
        go to 2
      else
        go to 3
      endif
      go to 5
    4 continue
      close(38)
      write(*,'(3a)') 'station ',astn,
     1      ' not in station coordinate file'
      ii=ii-1
      go to 2
    5 continue
      close(38)
      nsreq=ii-1
      nstns=nsreq
c      write(*,*) 'nsreq = ',nsreq
c      do 32 k=1,nsreq
c      write(*,'(a,a4,1x,a3,1x,2f15.7,2f5.1,2f15.4,5(1x,a1))')
c     1 'call phsort - ',
c     2 stn(k),stream(k),syf(k),sxf(k),psaf(k),ssaf(k),pt(k),st(k),
c     3 iwtp(k),iwts(k),apo(k),aptype(k),apfm(k)
c   32 continue
c      write(*,'(a)') 'call phsort'
      call phsort(nsreq,stn,stream,syf,sxf,psaf,ssaf,pt,st,iwtp,iwts,
     1            apo,aptype,apfm)
      if(nsreq.lt.1000) then
        do 40 i=nsreq+1,1000
        etsmp(i)=0.d0
        apo(i)=' '
        iwtp(i)=' '
        iwts(i)=' '
        aptype(i)=' '
        apfm(i)=' '
        pt(i)=0.d0
        st(i)=0.d0
        stn(i)='    '
        sxf(i)=0.d0
        syf(i)=0.d0
   40   continue
      endif
      call stnfil(nsreq,stn,syf,sxf)
c
      return
      end

      subroutine ts_calc(fpotime,t_shift) 
c
c ...... subtract modulo 86400 from time
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
c
      i_time=idint(fpotime)
      n_days=i_time/86400
      t_shift=86400.d0*dble(n_days)
c
      return
      end

      subroutine stnfil(nsreq,stn,syf,sxf)
c
c ...... write station coordinates to file for plotting
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      integer*4 getpid,ipid
      real*8 syf(1000),sxf(1000)
      character*4 stn(1000)
      character*80 tscfile
c
c ...... open temporary station coordinate file
c        used in plotting solution
c
      ipid=getpid()
      write(tscfile,'(a15,i5.5)') '/tmp/stn_coord.',ipid
c      write(*,*) 'tscfile = ',tscfile(1:lnblnk(tscfile))
      lfa=lnblnk(tscfile)
      open(39,file=tscfile(1:lfa),status='unknown')
c      write(*,*) 'nsreq = ',nsreq
      do 3 i=1,nsreq
      slatd=datand(dtand(syf(i))/0.993277d0)
      write(39,'(a4,1x,2f9.4)') stn(i),slatd,sxf(i)
c      write(*,'(i3,1x,a4,1x,2f9.4)') i,stn(i),slatd,sxf(i)
    3 continue
      close(39)
c
      return
      end

      subroutine l1dec(aline1,isyear,isdoy,ishour,ismin,fssec)
c
c ...... decode *.pf input line 2
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 aline1
      integer*4 icl(2),iblb(3),ible(3)
c
      inbb=1
      inbe=0
      iblb(inbb)=2
      do 1 i=4,lnblnk(aline1)
      if(aline1(i:i).eq.' ') then
        if(aline1(i-1:i-1).ne.' ') then
          inbe=inbe+1
          ible(inbe)=i
        endif
      else
        if(aline1(i-1:i-1).eq.' ') then
          inbb=inbb+1
          iblb(inbb)=i-1
          if(inbb.eq.3) go to 2
        endif
      endif
    1 continue
    2 continue
      inc=0
      do 3 i=iblb(3)+1,lnblnk(aline1)
      if(aline1(i:i).eq.':') then
        inc=inc+1
        icl(inc)=i
        if(inc.eq.2) go to 4
      endif
    3 continue
    4 continue
      read((aline1(iblb(1)+1:ible(1)-1)),'(i)') isyear
      read((aline1(iblb(2)+1:ible(2)-1)),'(i)') isdoy
      read((aline1(iblb(3)+1:icl(1) -1)),'(i)') ishour
      read((aline1(icl(1) +1:icl(2) -1)),'(i)') ismin
      il1a=lnblnk(aline1)
      read((aline1(icl(2)+1:il1a)),'(f)') fssec
c      write(*,'(a,4i10,f10.4)') 
c     1      'aline1 = ',isyear,isdoy,ishour,ismin,fssec
c
      return
      end

      subroutine l2dec(aline2,stn,stream,apo,aptype,apfm,
     1                 fpotime,ier)
c
c ...... decode *.pf input lines 3-n
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 aline2
      character*4 stn
      character*3 stream
      character*1 apo,aptype,apfm,apbt
      real*8 fpotime
      ier=0
      read(aline2(2:5),'(a)') stn
      if(stn(4:4).eq.'_') stn(4:4)=' '
      read(aline2(7:9),'(a)') stream
      inb=0
c
      jp=15
    1 jp=jp+1
      read(aline2(jp:jp),'(a)') apbt
      if(apbt.ne.' ') go to 1
c
      jp=jp+1
      read(aline2(jp:jp),'(a)') apo
      jp=jp+2
      read(aline2(jp:jp),'(a)') aptype
c
c ...... assume that a phase id of "p" or "s" are typos
c        and that the phase is really "P or "S"
c
      if(aptype.eq.'p') aptype='P'
      if(aptype.eq.'s') aptype='S'
c
      if((aptype.ne.'P').and.(aptype.ne.'S')) then
        ier=1
        return
      endif
    3 jp=jp+1
      if(aline2(jp:jp).ne.' ') go to 3
    4 jp=jp+1
      if(aline2(jp:jp).eq.' ') go to 4
      jp1=jp
    5 jp=jp+1
      if(aline2(jp:jp).ne.' ') go to 5
      jp2=jp-1
      read(aline2(jp1:jp2),'(f)') fpotime
      read(aline2(jp2+2:jp2+2),'(a)') apfm
c      write(*,'(a,3x,a4,3x,a3,3x,a1,3x,a1,3x,a1,3x,f10.4)')
c     1      'aline2 = ',stn,stream,apo,aptype,apfm,fpotime
c
      return
      end

      subroutine phsort(nsreq,stn,stream,syf,sxf,psaf,ssaf,pt,st,
     1                  wtp,wts,apo,aptype,apfm)
c
c ...... sort out phase data by station and time order
c
c        this routine converts the phase file output by SeisTool
c        into a format that RELP can use
c

      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      common/time/eqhour,isyear,isdoy,ishour,ismin,fssec
      integer*4 nsreq
      character*1 apo(1000),aptype(1000),apfm(1000),
     1            wtp(1000),wts(1000)
      character*1 apot,aptt,apft
      character*3 stream(1000)
      character*4 stn(1000),astn
      real*8 syf(1000),sxf(1000),psaf(1000),
     1       ssaf(1000),pt(1000),st(1000)
      real*8 syft(1000),sxft(1000),psat(1000),ssat(1000),pott(1000),
     1       sott(1000)
      character*1 apot(1000),aptt(1000),apft(1000),wtpt(1000),
     1            wtst(1000)
      character*4 stnt(1000)
      character*3 strt(1000)
	real*8 etyp_x,etyp_y,etyp_t
	integer*4 etyp_n,ol_ns
	common /etyp1/ etyp_x(1000),etyp_y(1000),etyp_t(1000),
     1                 etyp_n,ol_ns
	character*4 etyp_s,ol_stn
	common /etyp2/ etyp_s(1000),ol_stn(1000)
c
c      tbor=3600.d0*dble(ishour)+60.d0*dble(ismin)+fssec
	kk=0
      inumb=nsreq
      do 88 i=1,nsreq
      stnt(i)='    '
      strt(i)='   '
      syft(i)=0.d0
      sxft(i)=0.d0
      psat(i)=0.d0
      pott(i)=0.d0
      wtpt(i)='0'
      apot(i)=' '
      aptt(i)=' '
      apft(i)=' '
      ssat(i)=0.d0
      sott(i)=0.d0
      wtst(i)='0'
   88 continue
      j=1
      astn=stn(1)
      do 10 i=1,inumb
      jastn=1000000*ichar(astn(1:1))+10000*ichar(astn(2:2))+
     1      100*ichar(astn(3:3))+ichar(astn(4:4))
      istn=1000000*ichar(stn(i)(1:1))+10000*ichar(stn(i)(2:2))+
     1     100*ichar(stn(i)(3:3))+ichar(stn(i)(4:4))
      if(istn.eq.jastn) then
        if(aptype(i).eq.'P') then
	  kk=kk+1
	  etyp_s(kk)=stn(i)
	  etyp_x(kk)=sxf(i)
	  etyp_y(kk)=syf(i)
	  etyp_t(kk)=pt(i)
          stnt(j)=stn(i)
          strt(j)=stream(i)
          syft(j)=syf(i)
          sxft(j)=sxf(i)
          psat(j)=psaf(i)
          pott(j)=pt(i)
          wtpt(j)='1'
          apot(j)=apo(i)
          aptt(j)=aptype(i)
          apft(j)=apfm(i)
        elseif(aptype(i).eq.'S') then
	  kk=kk+1
	  etyp_s(kk)=stn(i)
	  etyp_x(kk)=sxf(i)
	  etyp_y(kk)=syf(i)
	  etyp_t(kk)=st(i)
          stnt(j)=stn(i)
          strt(j)=stream(i)
          syft(j)=syf(i)
          sxft(j)=sxf(i)
          ssat(j)=ssaf(i)
          sott(j)=st(i)
          wtst(j)='1'
          if(aptt(j).eq.' ') aptt(j)='S'
        endif
      else
        j=j+1
        astn=stn(i)
        if(aptype(i).eq.'P') then
	  kk=kk+1
	  etyp_s(kk)=stn(i)
	  etyp_x(kk)=sxf(i)
	  etyp_y(kk)=syf(i)
	  etyp_t(kk)=pt(i)
          stnt(j)=stn(i)
          strt(j)=stream(i)
          syft(j)=syf(i)
          sxft(j)=sxf(i)
          psat(j)=psaf(i)
          pott(j)=pt(i)
          wtpt(j)='1'
          apot(j)=apo(i)
          aptt(j)=aptype(i)
          apft(j)=apfm(i)
        elseif(aptype(i).eq.'S') then
	  kk=kk+1
	  etyp_s(kk)=stn(i)
	  etyp_x(kk)=sxf(i)
	  etyp_y(kk)=syf(i)
	  etyp_t(kk)=st(i)
          stnt(j)=stn(i)
          strt(j)=stream(i)
          syft(j)=syf(i)
          sxft(j)=sxf(i)
          ssat(j)=ssaf(i)
          sott(j)=st(i)
          wtst(j)='1'
          if(aptt(j).eq.' ') aptt(j)='S'
        endif
      endif
   10 continue
c
c ...... test pickfile data for event type and validity
c
	etyp_n=kk
	call chkevt
c
      nsreq=j
c
      do 13 i=1,nsreq
      stn(i)=stnt(i)
      stream(i)=strt(i)
      syf(i)=syft(i)
      sxf(i)=sxft(i)
      psaf(i)=psat(i)
      if(pott(i).gt.0.d0) then
        pt(i)=pott(i)
      else
        pt(i)=0.d0
      endif
      wtp(i)=wtpt(i)
      apo(i)=apot(i)
      aptype(i)=aptt(i)
      apfm(i)=apft(i)
      ssaf(i)=ssat(i)
      if(sott(i).gt.0.d0) then
        st(i)=sott(i)
      else
        st(i)=0.d0
      endif
      wts(i)=wtst(i)
   13 continue
c
c      do 12 i=1,nsreq
c      write(*,'(a,1x,a4,1x,a3,2f15.7,2f5.1,2f15.4,1x,5(a1,1x))')
c     1 'phsort - ',
c     2 stn(i),stream(i),syf(i),sxf(i),psaf(i),ssaf(i),pt(i),st(i),
c     3 wtp(i),wts(i),apo(i),aptype(i),apfm(i)
c   12 continue
c
      return
      end

	subroutine chkevt
c
c ...... flag events which appear not to be local or regional events
c
c	evt_flag
c		= 0 -> local or regional earthquake
c		= 1 -> noise glitch
c		= 2 -> non-clustered stations
c		= 4 -> teleseism
c		= 6 -> 2 + 4 above
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 x(1000),y(1000),t(1000)
	real*8 etyp_x,etyp_y,etyp_t
	integer*4 etyp_n,ol_ns
	common /etyp1/ etyp_x(1000),etyp_y(1000),etyp_t(1000),
     1                 etyp_n,ol_ns
	character*4 etyp_s,ol_stn
	common /etyp2/ etyp_s(1000),ol_stn(1000)
	integer*4 evt_flag
	common /etyp3/ evt_flag
	real*8 slow,se_slow,teq_del,teq_az,teq_lat,teq_lon
	integer*4 iteq_hr,iteq_mn,iteq_sc
	common /etyp4/ slow,se_slow,teq_del,teq_az,teq_lat,teq_lon,
     1                 iteq_hr,iteq_mn,iteq_sc
	real*8 ep_dt,ep_tt,ep_st,ep_ct
	common /eqtyp5/ ep_dt,ep_tt,ep_st,ep_ct
	real*8 krv,sekrv
	common /etyp6/ krv,sekrv
c
c	write(*,'(a,i5,3f15.5)') 
c     1       ('etyp - i,x,y,t = ',i,etyp_x(i),etyp_y(i),etyp_t(i),
c     2        i=1,etyp_n)
c
c ...... convert to cartesian coordinates
c
	call cordcon(x0,y0,t0,x,y,t,n)
c
c ...... test for time clustering
c
	evt_flag=0
	call cluster(t,n,ep_tt,qti)
	if(qti.gt.0.d0) then
c
c ...... probable telemetry noise glitch
c
	  evt_flag=evt_flag+1
	  return
	endif
c
c ...... test for spatial anti-clustering
c
	call spaceout(etyp_s,x,y,n,ol_ns,ol_stn)
	if(ol_ns.ge.1) then
	  evt_flag=evt_flag+2
	endif
c
c ...... test for teleseism
c
	if(evt_flag.ge.2) then
	  call remols(etyp_s,x,y,t,n,ol_ns,ol_stn)
	endif
	call pase(n,x,y,t)
c	write(*,*) 'slow & se_slow = ',slow,se_slow
	arg1=slow+2.d0*se_slow
	arg2=dabs(krv)+2.d0*sekrv
	if((arg1.lt.ep_st).and.(arg2.le.ep_ct)) then
	  evt_flag=evt_flag+4
	endif
c
	return
	end

	  subroutine remols(etyp_s,x,y,t,n,ol_ns,ol_stn)
c
c ...... remove outlying stations from list
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 x(1000),y(1000),t(1000)
	real*8 xx(1000),yy(1000),tt(1000)
	character*4 ss(1000)
	character*4 etyp_s(1000),ol_stn(1000)
	integer*4 n,ol_ns
c
	i=0
	k=0
    1   i=i+1
	if(i.gt.n) go to 3
	do 2 j=1,ol_ns
	if(etyp_s(i).eq.ol_stn(j)) then
	  go to 1
	endif
    2   continue
	k=k+1
	ss(k)=etyp_s(i)
	xx(k)=x(i)
	yy(k)=y(i)
	tt(k)=t(i)
	go to 1
    3   continue
	n=k
	do 4 i=1,n
	etyp_s(i)=ss(i)
	x(i)=xx(i)
	y(i)=yy(i)
	t(i)=tt(i)
    4   continue	
c
	return
	end

	subroutine cluster(t,n,ep_tt,qti)
c
c ...... determine median inter-phase time differences
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 t(1000)
	integer*4 n
	real*8 tt(1000)
c
c ...... check for ascending time order
c
	do 1 i=1,n
	tt(i)=t(i)
    1   continue
	do 3 i=1,n-1
	do 2 j=i+1,n
	if(tt(j).lt.tt(i)) then
	  tmp=tt(j)
	  tt(j)=tt(i)
	  tt(i)=tmp
	endif
    2   continue
    3   continue
c
c ...... look for telemetry glitch which is identified by several
c	 stations in a row (5 or more, say) which have the same (within
c	 ep_tt sec) onset times
c
	qti=0.d0
	ns=1
	i=0
   11   i=i+1
	if(ns.ge.5) then
	  qti=1.d0
	  go to 13
	endif
	ns=1
	if(i.ge.n) return
	j=i
   12   j=j+1
	if(j.gt.n) then
	  if(ns.ge.5) then
	    qti=1.d0
	    go to 13
	  endif
	  go to 11
	endif
	if((tt(j)-tt(i)).le.ep_tt) then
	  ns=ns+1
	  go to 12
	else
	  go to 11
	endif
c
   13   continue
c
	return
	end

	subroutine spaceout(etyp_s,x,y,n,ol_ns,ol_stn)
c
c ...... test for non-clustered stations
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 x(1000),y(1000),xt(1000),yt(1000)
	character*4 etyp_s(1000),ol_stn(1000),ol_stn_x(1000),
     1              ol_stn_y(1000),ol_stn_t
	integer*4 n,ol_ns
	real*8 ep_dt,ep_tt,ep_st,ep_ct
	common /eqtyp5/ ep_dt,ep_tt,ep_st,ep_ct
c
	do 1 i=1,n
	xt(i)=x(i)
	yt(i)=y(i)
	ol_stn_x(i)=etyp_s(i)
	ol_stn_y(i)=etyp_s(i)
    1   continue
c
c ...... put in ascending order
c
	do 3 i=1,n-1
	do 2 j=i+1,n
	if(xt(j).lt.xt(i)) then
	  xtt=xt(i)
	  xt(i)=xt(j)
	  xt(j)=xtt
	  ol_stn_t=ol_stn_x(i)
	  ol_stn_x(i)=ol_stn_x(j)
	  ol_stn_x(j)=ol_stn_t
	endif
	if(yt(j).lt.yt(i)) then
	  ytt=yt(i)
	  yt(i)=yt(j)
	  yt(j)=ytt
	  ol_stn_t=ol_stn_y(i)
	  ol_stn_y(i)=ol_stn_y(j)
	  ol_stn_y(j)=ol_stn_t
	endif
    2   continue
    3   continue
c	do 23 ij=1,n
c	write(*,'(a,2f15.4)') 
c     1    'spaceout - x & y = ',xt(ij),yt(ij)
c   23   continue
c
c ...... quartiles
c
	if(n.ge.4) then
	  q1x=xt(n/4)
	  q2x=xt(n/2)
	  q3x=xt((3*n)/4)
	else
	  q1x=xt(1)
	  q2x=(xt(1)+xt(n))/2.d0
	  q3x=xt(n)
	endif
	if(n.ge.4) then
	  q1y=yt(n/4)
	  q2y=yt(n/2)
	  q3y=yt((3*n)/4)
	else
	  q1y=yt(1)
	  q2y=(yt(1)+yt(n))/2.d0
	  q3y=yt(n)
	endif
c
c ....... stations more than ep_dt times the inter-quartile distance
c	  are assumed to be co-incident detections and not associated
c	  with the event
c
	qiqx=ep_dt*dabs(q3x-q1x)
	if(qiqx.le.50.d0) qiqx=50.d0
	qiqy=ep_dt*dabs(q3y-q1y)
	if(qiqy.le.50.d0) qiqy=50.d0
c	write(*,'(a,4f10.4)') 'q1x,q2x,q3x,qiqx = ',q1x,q2x,q3x,qiqx
c	write(*,'(a,4f10.4)') 'q1y,q2y,q3y,qiqy = ',q1y,q2y,q3y,qiqy
	ol_ns=0
	do 4 i=1,n
	if(xt(i).lt.(q2x-qiqx)) then
	  ol_ns=ol_ns+1
	  ol_stn(ol_ns)=ol_stn_x(i)
	elseif(xt(i).gt.(q2x+qiqx)) then
	  ol_ns=ol_ns+1
	  ol_stn(ol_ns)=ol_stn_x(i)
	elseif(yt(i).lt.(q2y-qiqy)) then
	  ol_ns=ol_ns+1
	  ol_stn(ol_ns)=ol_stn_y(i)
	elseif(yt(i).gt.(q2y+qiqy)) then
	  ol_ns=ol_ns+1
	  ol_stn(ol_ns)=ol_stn_y(i)
	endif
    4   continue
c
	return
	end

	subroutine cordcon(x0,y0,t0,x,y,t,n)
c
c ...... convert to a local cartesian coordinate frame
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 x(1000),y(1000),t(1000)
	real*8 etyp_x,etyp_y,etyp_t
	integer*4 etyp_n,ol_ns
	common /etyp1/ etyp_x(1000),etyp_y(1000),etyp_t(1000),
     1                 etyp_n,ol_ns
c
	n=etyp_n
	sumy=0.d0
	sumt=0.d0
	do 1 i=1,n
	sumy=sumy+etyp_y(i)
	sumt=sumt+etyp_t(i)
    1   continue
	y0d=sumy/dble(n)
	t0=sumt/dble(n)
c
	sumx=0.d0
	sumy=0.d0
	do 2 i=1,n
	sumx=sumx+111.195d0*etyp_x(i)*dcosd((etyp_y(i)+y0d)/2.d0)
	sumy=sumy+111.195d0*etyp_y(i)
    2   continue
	x0=sumx/dble(n)
	x0deg=x0/(111.195d0*dcosd(y0d))
	y0=sumy/dble(n)
	y0deg=y0/111.195d0
c
c	write(*,'(a,3f15.5)') 'x0, y0 & t0 = ',x0,y0,t0
c
	do 3 i=1,n
	x(i)=111.195d0*dcosd((etyp_y(i)+y0d)/2.d0)*etyp_x(i)-x0
	y(i)=111.195d0*etyp_y(i)-y0
	t(i)=etyp_t(i)-t0
c
c	write(*,'(a,i5,3f15.5)') 
c     1        'i, x(i), y(i), t(i) = ',i,x(i),y(i),t(i)
    3   continue
c
	return
	end

	subroutine pase(n,x,y,t)
c
c ...... determine apparent slowness and azimuth of assumed teleseism
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 x(1000),y(1000),t(1000)
	real*8 slow,se_slow,teq_del,teq_az,teq_lat,teq_lon
	integer*4 iteq_hr,iteq_mn,iteq_sc
	common /etyp4/ slow,se_slow,teq_del,teq_az,teq_lat,teq_lon,
     1                 iteq_hr,iteq_mn,iteq_sc
	real*8 krv,sekrv
	common /etyp6/ krv,sekrv
c
c ...... remove means
c
	x0=0.d0
	y0=0.d0
	t0=0.d0
	do 1 i=1,n
	x0=x0+x(i)
	y0=y0+y(i)
	t0=t0+t(i)
    1   continue
	x0=x0/dble(n)
	y0=y0/dble(n)
	t0=t0/dble(n)
	do 2 i=1,n
	x(i)=x(i)-x0
	y(i)=y(i)-y0
	t(i)=t(i)-t0
    2   continue
c
	atax11=0.d0
	atax12=0.d0
	atax21=0.d0
	atax22=0.d0
	atbx1=0.d0
	atbx2=0.d0
	do 4 i=1,n
	atax11=atax11+x(i)*x(i)
	atax12=atax12+x(i)*x(i)*x(i)
	atax21=atax21+x(i)*x(i)*x(i)
	atax22=atax22+x(i)*x(i)*x(i)*x(i)
	atbx1=atbx1+x(i)*t(i)
	atbx2=atbx2+x(i)*x(i)*t(i)
    4   continue
	denom=atax11*atax22-atax12*atax21
	atai11= atax22/denom
	atai12=-atax12/denom
	ata121=-atax21/denom
	atai22= atax11/denom
	ax=atai11*atbx1+atai12*atbx2
	bx=atai21*atbx1+atai22*atbx2
c
	sum=0.d0
	do 5 i=1,n
	sum=sum+(t(i)-ax*x(i)-bx*x(i)*x(i))*
     1          (t(i)-ax*x(i)-bx*x(i)*x(i))
    5   continue
	varx=sum/dble(n-1)
	varax=atai11*varx
	varbx=atai22*varx
	seax=dsqrt(varax)
	sebx=dsqrt(varbx)
c
	atay11=0.d0
	atay12=0.d0
	atay21=0.d0
	atay22=0.d0
	atby1=0.d0
	atby2=0.d0
	do 6 i=1,n
	atay11=atay11+y(i)*y(i)
	atay12=atay12+y(i)*y(i)*y(i)
	atay21=atay21+y(i)*y(i)*y(i)
	atay22=atay22+y(i)*y(i)*y(i)*y(i)
	atby1=atby1+y(i)*t(i)
	atby2=atby2+y(i)*y(i)*t(i)
    6   continue
	denom=atay11*atay22-atay12*atay21
	atai11= atay22/denom
	atai12=-atay12/denom
	ata121=-atay21/denom
	atai22= atay11/denom
	ay=atai11*atby1+atai12*atby2
	by=atai21*atby1+atai22*atby2
c
	sum=0.d0
	do 7 i=1,n
	sum=sum+(t(i)-ay*y(i)-by*y(i)*y(i))*
     1          (t(i)-ay*y(i)-by*y(i)*y(i))
    7   continue
	vary=sum/dble(n-1)
	varay=atai11*vary
	varby=atai22*vary
	seay=dsqrt(varay)
	seby=dsqrt(varby)
c
c ...... determine curvature of wavefront
c
c	write(*,'(a,4f15.5)') 
c     1    '  ax,   bx,   ay &   by = ',ax,bx,ay,by
c	write(*,'(a,4f15.5)') 
c     1    'seax, sebx, seay & seby = ',seax,sebx,seay,seby
c
	call kurve(ax,bx,ay,by,seax,sebx,seay,seby,krv,sekrv)
	varkrv=sekrv*sekrv
	if(varkrv.gt.0.d0) then
	  sekrv=dsqrt(varkrv/dble(n))
	else
	  sekrv=0.d0
	endif
c
	p=ax*ax+ay*ay
	varp=varax+varay
	if(p.gt.0.d0) p=dsqrt(p)
	if(varp.gt.0.d0) then
	  semp=dsqrt(varp/dble(n))
	else
	  semp=0.d0
	endif
	azim=datan2d(ax,ay)
	if(azim.lt.0.d0) azim=azim+360.d0
	u=ax/ay
	con=1.d0/(1.d0+u*u)
	semaz=con*con*(varax/(ay*ay)+varay*ax*ax/(ay*ay*ay*ay))
	if(semaz.gt.0.d0) semaz=dsqrt(semaz)
c	write(*,'(a,6f10.4)') 
c     1    'p, semp, az, semaz, k & semk = ',
c     2    p,semp,azim,semaz,krv,sekrv
c
	slow=p
	se_slow=semp
c
	if((p-2.d0*semp).lt.(1.d0/8.2d0)) then
	  p=111.195d0*p
	  semp=111.195d0*semp
	  semaz=semaz*45.d0/datan(1.d0)
	  call teloc(x0deg,y0deg,t0,p,azim,semp,semaz,
     1               eqdel,eqaz,eqlat,eqlon,eqtime,
     2               seqdel,seqlat,seqlon,seqtime)
c	  write(*,*) ' '
c	  write(*,*) '---------------------------------------------'
c	  write(*,*) 'Possible Teleseism'
c	  write(*,'(a,f8.5,a,f8.5,a)') 
c     1                               'p      = ',p,    ' +/- ',semp,
c     2          ' sec/deg'
c	  write(*,'(a,f8.2,a,f8.2)') 'az     = ',eqaz, ' +/- ',semaz
c	  write(*,'(a,f8.2,a,f8.2)') 'eqdel  = ',eqdel,' +/- ',seqdel
c	  write(*,'(a,f8.2)')        'eqlat  = ',eqlat
c	  write(*,'(a,f8.2)')        'eqlon  = ',eqlon
	  if(eqtime.lt.0.d0) eqtime=eqtime+86400.d0
	  ieqhr=eqtime/3600.d0
	  ieqmn=(eqtime-3600.d0*dble(ieqhr))/60.d0
	  ieqsc= eqtime-3600.d0*dble(ieqhr)-60.d0*dble(ieqmn)
c
	  teq_del=eqdel
	  teq_az=eq_az
	  teq_lat=eq_lat
	  teq_lon=eqlon
	  iteq_hr=ieqkr
	  iteq_mn=ieqmn
	  iteq_sc=ieqsc
c
c	  write(*,'(a,i2.2,a,i2.2,a,i2.2,a,f8.2)')
c     1          'eqtime = ',ieqhr,':',ieqmn,':',ieqsc,' +/- ',seqtime
c	  write(*,*) '---------------------------------------------'
c	  write(*,*) ' '
	endif
c
	return
	end

	subroutine kurve(ax,bx,ay,by,seax,sebx,seay,seby,k,sek)
c
c ...... calculate wavefront curvature & standard deviation
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 ax,bx,ay,by,seax,sebx,seay,seby,k,sek
c
	varax=seax*seax
	varbx=sebx*sebx
	xp=-1.d0/ax
	dxpdax=-xp/ax
	varxp=dxpdax*dxpdax*varax
	xpp=2.d0*bx/(ax*ax*ax)
	dxppdax=-3.d0*xpp/ax
	dxppdbx=xpp/bx
	varxpp=dxppdax*dxppdax*varax+dxppdbx*dxppdbx*varbx
c
	varay=seay*seay
	varby=seby*seby
	yp=-1.d0/ay
	dypday=-yp/ay
	varyp=dypday*dypday*varay
	xpp=2.d0*by/(ay*ay*ay)
	dyppday=-3.d0*xpp/ay
	dyppdby=xpp/by
	varypp=dyppday*dyppday*varay+dyppdby*dyppdby*varby
c
	u=xp*ypp-xpp*yp
	v=dsqrt((xp*xp+yp*yp)**3)
	k=u/v
	dvdxp=3.d0*xp*dsqrt(xp*xp+yp*yp)
	dvdyp=3.d0*yp*dsqrt(xp*xp+yp*yp)
	dkdxp=( v*ypp-u*dvdxp)/(v*v)
	dkdxpp=-yp/v
	dkdyp=(-v*xpp-u*dvdyp)/(v*v)
	dkdypp=xp/v
	vark=dkdxp*dkdxp*varxp+
     1       dkdxpp*dkdxpp*varxpp+
     2       dkdyp*dkdyp*varyp+
     3       dkdypp*dkdypp*varypp
	sek=dsqrt(vark)
c
	return
	end

	subroutine teloc(x0deg,y0deg,t0,p,azim,semp,semaz,
     1             eqdel,eqaz,eqlat,eqlon,eqtime,
     2             seqdel,seqlat,seqlon,seqtime)
c
c ...... locate probable teleseism
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 x0deg,y0deg,t0,p,azim,semp,semaz,
     1             eqdel,eqaz,eqlat,eqlon,eqtime,
     2             seqlat,seqlon,seqtime
c
	call eqtdel(p,semp,eqdel,seqdel,ptt,septt)
	eqaz=azim+180.d0
	if(eqaz.ge.360.d0) eqaz=eqaz-360.d0
	call strig(eqdel,seqdel,eqaz,semaz,x0deg,y0deg,
     1             eqlat,eqlon,seqlat,seqlon)
	eqtime=t0-ptt
	seqtime=septt
c
	return
	end

	subroutine eqtdel(p,semp,eqdel,seqdel,ptt,septt)
c
c ...... calculate delta & P travel time (assuming surface focus)
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 p,semp,eqdel,seqdel,ptt,septt
c
	vp=semp*semp
c
	call dtcomp(p,vp,d,vd,t,vt)
c
	eqdel=d
	seqdel=vd
	if(seqdel.gt.0.d0) seqdel=dsqrt(seqdel)
	ptt=t
	septt=vt
	if(vt.gt.0.d0) septt=dsqrt(vt)
c
	return
	end

	subroutine strig(eqdel,seqdel,eqaz,semaz,x0deg,y0deg,
     1             eqlat,eqlon,seqlat,seqlon)
c
c ...... spherical trig
c
	implicit real*8 (a-h,o-z)
	implicit integer*4 (i-n)
	real*8 eqdel,seqdel,eqaz,semaz,x0deg,y0deg,
     1         eqlat,eqlon,seqlat,seqlon
c
	arg=dcosd(90.d0-y0deg)*dcosd(eqdel)+
     1      dsind(90.d0-y0deg)*dsind(eqdel)*dcosd(eqaz)
	eqlat=90.d0-dacosd(arg)
	sarg=dsind(eqdel)*dsind(eqaz)/dsind(90.d0-y0deg)
	carg=(dcosd(eqdel)-dcosd(90.d0-y0deg)*dcosd(90.d0-eqlat))/
     1       (dsind(90.d0-y0deg)*dsind(90.d0-eqlat))
	eqlon=x0deg+datan2d(sarg,carg)
	if(eqlon.lt.-180.d0) eqlon=eqlon+360.d0
	if(eqlon.gt. 180.d0) eqlon=eqlon-360.d0
c
	
c
	return
	end

      subroutine dtcomp(p,vp,d,vd,t,vt)
c
c ...... compute distance and origin time using JB P Tables
c
      implicit real*8 (a-h,o-z)
      common/pdt/pjb(220),djb(220),tjb(220)
c
      i=1
    1 i=i+1
      if(pjb(i).gt.p) go to 1
      dddp=(djb(i)-djb(i-1))/(pjb(i)-pjb(i-1))
      d=djb(i-1)+dddp*(pjb(i-1)-p)
      vd=dddp*dddp*vp
      dtdp=(tjb(i)-tjb(i-1))/(pjb(i-1)-pjb(i))
      t=tjb(i-1)+dtdp*(pjb(i-1)-p)
      vt=dtdp*dtdp*vp
c
c	write(*,*) 'dtcomp - '
c	write(*,*) 'p,vp,d,vd,t,vt = ',p,vp,d,vd,t,vt
c
      return
      end

      block data
c
c ...... JB P-wave travel time, distance and slowness tables
c
      implicit real*8 (a-h,o-z)
      common/pdt/pjb(220),djb(220),tjb(220)
c
      data (pjb(i),i=1,81)/
     *19.856,17.060,14.263,14.232,14.191,14.141,14.082,14.012,13.930,
     *13.835,13.726,13.602,13.465,13.315,13.153,12.970,12.747,12.455,
     *12.075,11.612,11.109,10.633,10.236, 9.933, 9.705, 9.520, 9.360,
     * 9.220, 9.099, 8.997, 8.909, 8.831, 8.760, 8.691, 8.623, 8.555,
     * 8.490, 8.431, 8.376, 8.320, 8.262, 8.203, 8.150, 8.106, 8.068,
     * 8.023, 7.964, 7.889, 7.809, 7.731, 7.657, 7.584, 7.510, 7.433,
     * 7.355, 7.274, 7.191, 7.106, 7.023, 6.941, 6.861, 6.781, 6.700,
     * 6.621, 6.544, 6.470, 6.397, 6.324, 6.250, 6.175, 6.100, 6.025,
     * 5.950, 5.875, 5.800, 5.725, 5.650, 5.575, 5.500, 5.425, 5.350/
      data (pjb(i),i=82,162)/
     * 5.276, 5.204, 5.135, 5.067, 5.000, 4.933, 4.867, 4.802, 4.742,
     * 4.689, 4.649, 4.621, 4.602, 4.587, 4.571, 4.554, 4.540, 4.529,
     * 4.519, 4.505, 4.485, 4.459, 4.455, 4.451, 4.447, 4.447, 4.445,
     * 4.443, 4.441, 4.439, 4.437, 4.435, 4.433, 4.431, 4.427, 4.419,
     * 4.406, 4.390, 4.374, 4.362, 4.355, 4.350, 4.346, 4.339, 4.330,
     * 4.319, 4.309, 4.298, 4.286, 4.270, 4.250, 4.229, 4.208, 4.187,
     * 4.158, 4.116, 4.058, 3.983, 3.893, 3.786, 3.656, 3.501, 3.180,
     * 3.180, 2.837, 2.650, 2.463, 2.224, 1.978, 1.970, 1.966, 1.962,
     * 1.957, 1.950, 1.944, 1.941, 1.940, 1.941, 1.940, 1.938, 1.935/
      data (pjb(i),i=163,220)/
     * 1.932, 1.928, 1.926, 1.925, 1.924, 1.920, 1.914, 1.905, 1.893,
     * 1.880, 1.867, 1.854, 1.844, 1.836, 1.829, 1.818, 1.801, 1.780,
     * 1.759, 1.741, 1.729, 1.719, 1.705, 1.683, 1.653, 1.618, 1.580,
     * 1.539, 1.496, 1.452, 1.410, 1.375, 1.344, 1.313, 1.273, 1.221,
     * 1.160, 1.101, 1.052, 1.014, 0.978, 0.934, 0.873, 0.800, 0.727,
     * 0.666, 0.620, 0.582, 0.540, 0.486, 0.422, 0.354, 0.287, 0.224,
     * 0.168, 0.123, 0.087, 0.000/
      data (djb(i),i=1,120)/
     *  0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11.,
     * 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23.,
     * 24., 25., 26., 27., 28., 29., 30., 31., 32., 33., 34., 35.,
     * 36., 37., 38., 39., 40., 41., 42., 43., 44., 45., 46., 47.,
     * 48., 49., 50., 51., 52., 53., 54., 55., 56., 57., 58., 59.,
     * 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., 70., 71.,
     * 72., 73., 74., 75., 76., 77., 78., 79., 80., 81., 82., 83.,
     * 84., 85., 86., 87., 88., 89., 90., 91., 92., 93., 94., 95.,
     * 96., 97., 98., 99.,100.,101.,102.,103.,104.,105.,180.,179.,
     *178.,177.,176.,175.,174.,173.,172.,171.,170.,169.,168.,167./
      data (djb(i),i=121,220)/
     *166.,165.,164.,163.,162.,161.,160.,159.,158.,157.,156.,155.,
     *154.,153.,152.,151.,150.,149.,148.,147.,146.,145.,144.,143.,
     *143.,144.,145.,146.,147.,110.,111.,112.,113.,114.,115.,116.,
     *117.,118.,119.,120.,121.,122.,123.,124.,125.,126.,127.,128.,
     *129.,130.,131.,132.,133.,134.,135.,136.,137.,138.,139.,140.,
     *141.,142.,143.,144.,145.,146.,147.,148.,149.,150.,151.,152.,
     *153.,154.,155.,156.,157.,158.,159.,160.,161.,162.,163.,164.,
     *165.,166.,167.,168.,169.,170.,171.,172.,173.,174.,175.,176.,
     *177.,178.,179.,180./
      data (tjb(i),i=1,81)/
     *   6.8,  21.1,  35.4,  49.7,  63.9,  78.1,  92.2, 106.3, 120.3,
     * 134.2, 148.0, 161.7, 175.3, 188.7, 201.9, 215.0, 228.0, 240.7,
     * 253.2, 265.5, 277.0, 287.4, 297.5, 307.4, 317.1, 326.8, 336.2,
     * 345.4, 354.5, 363.5, 372.5, 381.3, 390.1, 398.8, 407.5, 416.1,
     * 424.6, 433.0, 441.4, 449.8, 458.1, 466.3, 474.5, 482.7, 490.3,
     * 498.9, 506.8, 514.7, 522.6, 530.3, 538.0, 545.6, 553.2, 560.7,
     * 568.0, 575.4, 582.6, 589.8, 596.8, 603.8, 610.7, 617.5, 624.3,
     * 630.9, 637.5, 644.0, 650.4, 656.8, 663.1, 669.3, 675.4, 681.5,
     * 687.5, 693.4, 699.2, 705.0, 710.7, 716.3, 721.8, 727.3, 732.7/
      data (tjb(i),i=82,162)/
     * 738.0, 743.2, 748.4, 753.5, 758.5, 763.5, 768.4, 773.2, 778.0,
     * 782.7, 787.3, 791.9, 796.5, 801.1, 805.7, 810.3, 814.8, 819.3,
     * 823.8, 828.4, 832.9, 837.4, 841.8, 846.2, 850.6,1330.6,1326.2,
     *1321.8,1317.4,1312.9,1308.5,1304.1,1299.6,1295.2,1290.8,1286.3,
     *1281.9,1277.5,1273.1,1268.8,1264.4,1260.1,1255.7,1251.4,1247.0,
     *1242.7,1238.4,1234.1,1229.8,1225.5,1221.2,1217.0,1212.8,1208.6,
     *1204.4,1200.2,1196.1,1192.0,1188.1,1184.2,1180.4,1176.8,1173.5,
     *1173.5,1176.5,1179.3,1181.8,1184.1,1113.2,1115.2,1117.1,1119.1,
     *1121.1,1123.0,1125.0,1126.9,1128.8,1130.8,1132.7,1134.7,1136.6/
      data (tjb(i),i=163,220)/
     *1138.5,1140.5,1142.4,1144.3,1146.2,1148.2,1150.1,1152.0,1153.9,
     *1155.8,1157.7,1159.5,1161.4,1163.2,1165.0,1166.9,1168.7,1170.5,
     *1172.3,1174.0,1175.7,1177.4,1179.2,1180.9,1182.6,1184.2,1185.8,
     *1187.4,1188.9,1190.4,1191.8,1193.2,1194.5,1195.8,1197.2,1198.5,
     *1199.7,1200.8,1201.8,1202.8,1203.8,1204.8,1205.8,1206.7,1207.4,
     *1208.0,1208.6,1209.2,1209.8,1210.4,1210.9,1211.2,1211.5,1211.8,
     *1212.0,1212.1,1212.1,1212.2/
c
      end

