c hymain.for    []
      program hypoel
c hypoellipse -- john c. lahr
c non-unix version of main routine
c
c************************* notes for programmers **************************
c
c     a binary search of the station list is used in the function phaind.
c         if the search does not work on your computer, use the version
c         of phaind that is commented out, which does a complete search.
c
c     subroutines init, opfls, trvcon, trvdrv, and uamag must have
c     $notruncate statements added for pc systems so that variable names
c     longer than 6 characters may be used.
c
c     all non-unix systems use hymain.for and init.for.  the equivalent
c         routines for unix systems are hypoe.c, setup_server.c,
c         listen_serv.c, fdgetstr.c, cleanup.f, initial.f, and
c         getbin.f.
c
c     subroutines openfl and opfls, which open files, have differences between
c         various systems.
c
c     subroutine openfl also has differences between the unix-masscomp and
c         the unix-sun systems.
c
c     subroutines phagt and npunch use a back slash character, which
c         must be doubled on unix systems.
c
c     subroutines dubl, erset, jdate, and timit use non-standard fortran
c         and must be modified for use with vax, pc, or unix computers.
c
c           dubl     - sets a double precision number equal to a
c                        single precision number
c           erset    - resets error limits on vax/vms computers
c           jdate    - gets current date and time from the operating system
c           timit    - times the execution on vax/vms computers
c
c     alternat versions of the code is enclosed by 'c* unix', 'c* pc',
c         or 'c* vax' comment statements in each of the above subroutines.
c

c get filenames, open files and write greeting
      intype = 1
      call init

c read station list, crustal model, and control records
      call input1

c initialize summary of residuals
      call lissum(1)

c read arrival times and locate earthquakes
      call locate(-1)
      stop
      end
c end hymain
c adddly.for  []
      subroutine adddly(nmodel)
c read in delays for model number greater than 5
      include 'params.inc'
      parameter (ndly = 11)
      character*4 iahead*60, msta*5, nsta*5, icard*110
      integer punt, phaind
      common /punt/ punt
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /ilmpu/ ns
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /int/ thk(lmax+1),lbeg(mmax),lend(mmax),vi(lmax),vsq(lmax),
     *  vsqd(lmmax,lmax),f(lmax,lmmax),kl,ivway,sthk(mmax),sthk1(mmax),
     *  sdly(ndly,nsn)
      common /logfil/ logfil
      character*5 stard, dnstrg, rshft
      if(ns .eq. 0) then
        write(punt, 18)
        write(logfil, 62)
18      format('xxxerrorxxx delay models may not preceed the ',
     *  'station list')
        stop 'abort from adddly'
      endif
20    read(inpt, '(a)', end=60) icard
      stard = dnstrg(rshft(icard(1:5)))
      if(stard .eq. '  end') return
c get station number (i) for station stard
      if(phaind(stard, nsta, ns, i, ierr) .ne. 0) then
        write(punt, 34) stard, icard
        write(logfil, 34) stard, icard
34      format(' ***>', a5, ' is not on station list, so these delays ',
     *  ' will not be used:', /, a)
        goto 20
      endif
c because standard fortran 77 to can not read an internal file with
c free format, file 14 must be used in the following code!
      rewind 14
      write(14, '(a)') icard(6:110)
      rewind 14
c     read(icard(5:110), *) dly(nmodel, i), sdly(nmodel, i)
      read(14, *) dly(nmodel, i), sdly(nmodel, i)
      goto 20
60    write(punt, 62)
      write(logfil, 62)
62    format('xxxerrorxxx end of file found while reading additional',
     * ' delays, so stop')
      stop 'abort from adddly'
      end
c end adddly
c adderr.for    [pc]
      subroutine adderr(zup, zdn)
c adds zup and zdn to the summary record and adds the event
c to file 4, and file 11
      character*117 sumrec
      character*117 arcrec
      character*6 fmit
      common /an/ n14, n15
c* (unix
c      call flush (14)
c      call flush (15)
c* unix)
c* (pc
c* pc)
c* (vax
c* vax)
      rewind 14
      rewind 15
      if(n14 .gt. 0) then
        if(n14 .gt. 1) then
          print *, ' n14 = ', n14
          print *, ' logic error: n14 should be 0 or 1'
        endif
        do 20 i = 1, n14
          read(14, '(a)') sumrec
          call formal(zup, izup, 2, 0, fmit, azup)
          if (fmit .eq. ' ') then
c           write(sumrec(101:102), '(i2)') izup
            write(sumrec(103:104), '(i2)') izup
          else
c           write(sumrec(101:102), fmit) azup
            write(sumrec(103:104), fmit) azup
          endif
          call formal(zdn, izdn, 2, 0, fmit, azdn)
          if (fmit .eq. ' ') then
c           write(sumrec(103:104), '(i2)') izdn
            write(sumrec(105:106), '(i2)') izdn
          else
c           write(sumrec(103:104), fmit) azdn
            write(sumrec(105:106), fmit) azdn
          endif
          write(4, '(a)') sumrec(1:lentru(sumrec))
20      continue
      endif
      if(n15 .gt. 1) then
        do 30 i = 1, n15
          read(15, '(a)') arcrec
c         if(arcrec(81:81) .ne. '/') then
          if(arcrec(83:83) .ne. '/') then
            write(11, '(a)') arcrec(1:lentru(arcrec))
          else
c summary record
            sumrec = arcrec
            call formal(zup, izup, 2, 0, fmit, azup)
            if (fmit .eq. ' ') then
c             write(sumrec(101:102), '(i2)') izup
              write(sumrec(103:105), '(i2)') izup
            else
c             write(sumrec(101:102), fmit) azup
              write(sumrec(103:105), fmit) azup
            endif
            call formal(zdn, izdn, 2, 0, fmit, azdn)
            if (fmit .eq. ' ') then
c             write(sumrec(103:104), '(i2)') izdn
              write(sumrec(105:106), '(i2)') izdn
            else
c             write(sumrec(103:104), fmit) azdn
              write(sumrec(105:106), fmit) azdn
            endif
c	    print *, 'in adderr, about to write summary rec to unit 11'
c	    print *, 'sumrec = ', sumrec
c	    print *, 'lentru = ', lentru(sumrec)
            write(11, '(a)') sumrec(1:lentru(sumrec))
          endif
30      continue
      endif
c* (unix
c      call flush(4)
c      call flush(11)
c* unix)
c* (pc
c* pc)
c* (vax
c* vax)
      return
      end
c end adderr
c azwtos.for   []
      subroutine azwtos
c azimuthal weighting of stations by quadrants
      include 'params.inc'
      character*4 iahead*60, msta*5, nsta*5, icard*110
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /amo/ tempg(npa), ntgap
      common /gmost/ az(npa)
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /qmost/ wt(npa),z
      dimension tx(4),txn(4),ktx(4),kemp(npa),key(npa)
c count and sort by azimuth stations with weight .ne. 0
      j=0
      do 10 i=1,nr
      if (wt(i) .eq. 0.) goto 10
      j=j+1
      tempg(j)=az(i)
   10 continue
c divide into 4 zones with one axis bisecting largest gap between stations
      call sort(tempg,key,j)
      gap=tempg(1)+360.-tempg(j)
      ig=1
      do 20 i=2,j
      dtemp=tempg(i)-tempg(i-1)
      if (dtemp .le. gap) goto 20
      gap=dtemp
      ig=i
   20 continue
      tx(1)=tempg(ig)-0.5*gap
      tx(2)=tx(1)+90.
      tx(3)=tx(1)+180.
      tx(4)=tx(1)+270.
      do  30 i=1,4
      txn(i)=0.
      if (tx(i) .lt. 0.) tx(i)=tx(i)+360.
      if (tx(i).gt.360.) tx(i)=tx(i)-360.
   30 continue
      call sort(tx,ktx,4)
      do  80 i=1,nr
      if (wt(i) .eq. 0.) goto  80
      if (az(i) .gt. tx(1)) goto  50
   40 txn(1)=txn(1)+1.
      kemp(i)=1
      goto  80
   50 if (az(i) .gt. tx(2)) goto  60
      txn(2)=txn(2)+1.
      kemp(i)=2
      goto  80
   60 if (az(i) .gt. tx(3)) goto  70
      txn(3)=txn(3)+1.
      kemp(i)=3
      goto  80
   70 if (az(i) .gt. tx(4)) goto  40
      txn(4)=txn(4)+1.
      kemp(i) = 4
   80 continue
      xn=4
      if (txn(1).eq.0.) xn=xn-1
      if (txn(2).eq.0.) xn=xn-1
      if (txn(3).eq.0.) xn=xn-1
      if (txn(4).eq.0.) xn=xn-1
      fj=j/xn
      do  90 i=1,nr
      if (wt(i) .eq. 0.) goto  90
      ki=kemp(i)
c new weight equals old weight times total number of stations
c divided by number of zones contatning stations and divided
c------- again by number of zones containing this station
      wt(i)=wt(i)*fj/txn(ki)
   90 continue
      return
      end
c end azwtos
c back.for   []
      subroutine back (delat, delon, newlat, newlon, slat, slon)
c
c-------- back - calculate geocentric coordinates of secondary point from
c            step in latitude (km) and longitude(km)
c
c input:  delat     change in earthquake latitude in km (northward positive)
c         delon     change in earthquake longitude in km (westward positive)
c output: newlat    new earthquake geocentric latitude in radians
c         newlon    new earthquake longitude in radians
c
      real newlat,newlon
      parameter (pi = 3.14159265)
      parameter (twopi = 2.0*pi)
      parameter (halfpi = 0.5*pi)
      parameter (rad = pi/180.)
      parameter (deg = 1.0/rad)
      parameter (equrad = 6378.2064)
      parameter (polrad = 6356.5838)
      parameter (flat = (equrad - polrad)/equrad)
c
      st0 = cos(slat)
      ct0 = sin(slat)
      call cvrtop(delat,delon,delta,az1)
      if (az1 .lt. 0.0) az1 = az1 + twopi
c use approximation of local radius for derivative of surface
c     distance with geocentric latitude.
c  more accurate formulation would be:
c      drdth = -r**3 * cos(alat)*sin(alat)*( (1.-flat)**(-2) - 1. )/a**2
c      dsd = sqrt(drdth**2 + r**2)
      radius = (cos(slat)**2/equrad**2 + sin(slat)**2/polrad**2)**(-.5)
      sdelt = sin(delta/radius)
      cdelt = cos(delta/radius)
      cz0 = cos(az1)
      ct1 = st0*sdelt*cz0+ct0*cdelt
      call cvrtop(st0*cdelt-ct0*sdelt*cz0, sdelt*sin(az1), st1, dlon)
      newlat = atan2(ct1, st1)
      newlon = slon + dlon
      if (abs(newlon) .gt. pi) newlon = newlon - sign(twopi, newlon)
      return
      end
c end back
c block.for   []
      block data
c initialize constants in common statements
      include 'params.inc'
      parameter (ndly = 11)
      logical good, eoff, supout
      character*4 iahead*60, msta*5, nsta*5, icard*110, uacal*50
      common /char/ iahead, msta(npa), nsta(nsn), icard
      character*1 iclass
      common /dbio/ iclass(0:4)
      common /dbiln/ ioldq
      common /dhio/ nedit,good
      common /dio/ rmsmx,presmx,sresmx,noutmx,nimx,iprun,semx
      common /dph/ noswt, eoff
      common /dhin/ iglob, zup, zdn
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /dhil/ iq,ilat,kms,icat
      character*1 bksrc
      common /dipu/ ipkdly, igsum, bksrc
      common /dit/ tid(lmax,lmmax),did(lmax,lmmax),lowv,modin(mmax+3)
      logical medmag
      common /dinx/ imag,imagsv,medmag
      common /dix/ iexcal, uacal
      common /dmost/ ipun,ivlr,blank
      common /imost/ test(100)
      character*1 iqcls
      common /il1/ iqcls
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /logfil/ logfil
      common /idno/ ksel,ksort
      common /idt/ v(lmax2)
      common /iox/ prr(nsn),iuses
      common /it/ vpvs(lmax2),vpvsm(mmax+3),nttab
      common /ix/ ir,qspa(9,40)
      common /ohq/ gap, supout
      common /pot/ ts(npa),ihr,model(npa), keyphi(npa)
      common /reloc/ irelo, nreloc
      common /rioq/ been,damp,dmpinc,igo
      character*1 mgndx
      common /xfmno1/ mgndx
      common /ip1/ rsew(4)
      common /ilmpu/ ns
      data uacal/'pub1:[alaska.data]uofacal.dat'/, iexcal/0/
c     data rsew/1., 3., 7.5, 15./, inpt/8/, inmain/8/, injump/12/
c revised default relative weights 4/22/89
      data rsew/1., 5., 10., 20./, inpt/8/, inmain/8/, injump/12/
      data logfil/6/,isa/0/,iqcls/'b'/,damp/0.0010/,bksrc/' '/
      data iprun/1/, supout/.false./, ns/0/
      data iclass/'0','a','b','c','d'/, nttab/0/
      data iahead/'    '/,mgndx/' '/,good/.true./
      data blank/1.e20/
      data test/100*1.23456/, v/lmax2*0.001/, model/npa*1/
      data ivlr,kms,iprn,ipun,imag,iq,ksort,ksel,icat/0,1,1,0,0,2,0,1,1/
      data ioldq/0/,ir/0/,lowv/1/,imagsv/0/,iuses/0/,ipkdly/1/
cds      data tid,did/narray*0.0/, modin/13*0/, iglob/1/
      data tid,did/narray*0.0/, modin/mmaxp3*0/, iglob/1/
      data nedit/0/, eoff/.false./, igsum/1/, irelo/0/, nreloc/0/
      data ilis/1/,medmag/.false./
      end
c end block
c boxau.for    []
      subroutine boxau
c compute rms at auxiliary points and estimate quality
      include 'params.inc'
      parameter (ndly = 11)
      real latep,lonep,latsv,lonsv
      save latsv, lonsv, rmssv, zsv
      save dez, t6, drsv, kdrsv, diag, kdiag, sum, aalz, aala
      save aalo, idirec, idrc
      character*1 iqa, ins, iew
      character*4 iahead*60, msta*5, nsta*5, icard*110
      integer punt
      common /punt/ punt
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /bz/ na
      common /bqz/ avrps,avuse
      common /dmost/ ipun,ivlr,blank
      character*1 iclass
      common /dbio/ iclass(0:4)
      common /gmost/ az(npa)
      logical freor
      common /hopq/ savla,savlo,savez,savor,freor
      common /imost/ test(100)
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /obcfn/ ain(npa)
      character*1 kwr, iqdo
      common /pbnoq/ kwr(npa),iqdo(npa)
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /phoqn/ inst,knst
      common /qmost/ wt(npa),z
      common /qmost1/ lonep,ni,latep
      common /qgnotx/ delta(npa)
      common /qgo/ org
      common /rob/ v(4,4),noaxp
      common /rbno/ idip(3),iaaz(3),a(3,3)
      common /rfgnoq/ se(4)
      common /tmost/ x(4,npa)
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      common /zqr/ fno
      dimension drsv(14),kdrsv(14),diag(7),kdiag(7)
      dimension sum(5)
      dimension aalz(14),aala(14),aalo(14)
      dimension vec(3)
      character*2 idirec(7), idrc(7)
      data idirec/' n','se','sw','nw','ne',' z',' e'/
      data aala/1.732,1.,1.,1.,1.,0.,0.,0.,0.,-1.,-1.,-1.,-1.,-1.732/,
     *     aalo/0.,1.,-1.,1.,-1.,0.,1.732,-1.732,0.,1.,-1.,1.,-1.,0./,
     *     aalz/0.,-1.,-1.,1.,1.,-1.732,0.,0.,1.732,-1.,-1.,1.,1.,0./
      kgoon = 0
c
c write heading
   17 if( (iprn .gt. 2)  .or. (noaxp .eq. 1) ) write(punt,18)
   18 format(/50h       lat       lon         z     avrps        no ,
     *  7x,3hrms,25x,4hdrms/)
c
c save current values
      rmssv = rms
      latsv = latep
      lonsv = lonep
      zsv = z
      orgsv = org
      instsv = inst
c set inst = -9 so quakes will compute fixed location rms but will
c not call regres.
      inst = -9
      freor = .true.
      t6 = abs(test(6))/1.732
      dez = t6
      do 70 na = 1, 14
        call back(t6*aala(na),t6*aalo(na),savla,savlo,latsv,lonsv)
        savez = zsv +aalz(na)*dez
c       if(z .lt. 0.) z = 0.0
        if(savez .lt. 0.) savez = 0.0
        call quakes
c output rms error af auxiliary points
        rmsx = rms
        drms = rmsx - rmssv
        if(iprn .eq. 4) then
c calcute predicted rms at aux points
          tot = 0.0
          vec(1) = aalo(na)/1.732
          vec(2) = aala(na)/1.732
          vec(3) = aalz(na)/1.732
          do 48 i = 1,3
            do 46 j = 1,3
              tot = tot + vec(i)*vec(j)*v(i,j)
   46       continue
   48     continue
          pdrms = sqrt(rmssv**2 + tot*test(6)*test(6)/fno)
          drms = rmsx - pdrms
        endif
        call unfold2(savla,savlo,la,ins,ala,lo,iew,alo)
        drsv(na) = drms
        if( (iprn .ge. 2) .and. (noaxp .ne. 1) ) then
          goto (1,2,3,4,5,6,7,8,9,10,11,4,13,1), na
    2     write(punt,801) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  801     format(4f10.2,i10,f10.2,10x,f5.2)
          goto  50
    3     write(punt,802) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  802     format(4f10.2,i10,f10.2,24x,1h-,8x,f5.2/96x,1h.)
          goto  50
    1     write(punt,803) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  803     format(/4f10.2,i10,f10.2,23x,f5.2/)
          goto  50
    4     write(punt,804) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  804     format(4f10.2,i10,f10.2,13x,f5.2,19x,1h.)
          goto  50
    5     write(punt,805) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  805     format(4f10.2,i10,f10.2,12x,'|',23x,f5.2)
          goto  50
    6     write(punt,806) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  806     format(4f10.2,i10,f10.2,20x,f5.2,10x,'|')
          goto  50
    7     write(punt,807) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  807     format(4f10.2,i10,f10.2,3x,f5.2/)
          write(punt,815) rmssv
  815     format(50x,f10.2,12x,'|',10x,5h 0.00,10x,'|'/95x,'|')
          goto  50
    8     write(punt,808) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  808     format(4f10.2,i10,f10.2,43x,f5.2)
          goto  50
    9     write(punt,809) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  809     format(4f10.2,i10,f10.2,26x,f5.2)
          goto 50
   10     write(punt,810) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  810     format(4f10.2,i10,f10.2,10x,f5.2,23x,'|')
          goto 50
   11     write(punt,811) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  811     format(4f10.2,i10,f10.2,13x,1h.,10x,1h-,8x,f5.2/
     *    74x,1h.,21x,1h.)
          goto 50
   13     write(punt,813) ala,alo,z-test(8),avrps,nrwt,rmsx,drms
  813     format(4f10.2,i10,f10.2,27x,1h-,8x,f5.2)
        endif
   50   continue
   70 continue
c ouality output
      i = 0
      ii = 15
      do  80 iii = 1,7
      i = i + 1
      ii = ii - 1
   80 diag(iii) = (drsv(i) + drsv(ii))/2.0
      call sort(diag, kdiag, 7)
      call sort(drsv, kdrsv, 14)
      do  90 i = 1,7
        j = kdiag(i)
        idrc(i) = idirec(j)
   90 continue
      if((noaxp .eq. 0) .and. (iprn .ge. 0)) write(punt,100) idrc, diag
  100 format(/25x,18hquality evaluation//,
     * 34h diagonals in order of strength   ,7(4x,a2)/,
     *  19h ave. of end points,15x,7f6.2//)
      sum(1) = 0.0
      do 110 i=1,14
c do not use upper point in average if limited by surface
        if((z.lt.abs(test(6))*1.732).and.(kdrsv(i).eq. 6)) goto 110
        sum(1) = sum(1) + drsv(i)
  110 continue
      avdrms = sum(1)/14.
      if(z .lt. abs(test(6))*1.732) avdrms = sum(1)/13.
      drmin = drsv(1)
      icmin = kdrsv(1)
      if((z .ge. abs(test(6))*1.732).or.(kdrsv(1) .ne. 6)) goto 115
      drmin = drsv(2)
      icmin = kdrsv(2)
  115 jab = 4
      if((nrwt.ge. 4.).and.(rmssv.le. 0.4).and.(avdrms.ge. 0.5)) jab = 3
      if((nrwt.ge.5.).and.(rmssv.le. 0.4).and.(drmin  .ge. 0.15)) jab=2
      if((nrwt.ge.6.).and.(rmssv.le. 0.2).and.(drmin  .ge. 0.30)) jab=1
      iqa = iclass(jab)
      if((noaxp.eq.0) .and. (iprn .ge. 0))
     * write(punt,120)nrwt,rmssv,drmin,avdrms,iqa
  120 format(10x,50h    number       rms  min drms  ave drms   quality/
     * 10x,i10,3f10.2,9x,a1/)
c
c restore variable values
      latep = latsv
      lonep = lonsv
      z = zsv
      org = orgsv
      inst = instsv
      if (drmin .gt. -.03 .or. test(6) .gt. 0.) then
c restore all values, such as azimuth, angle of incidence, and
c standard error to those of final solution for npunch.
130     iprnsv = iprn
        iprn = -2
        inst = 9
        freor = .false.
        savla = latsv
        savlo = lonsv
        savez = zsv
        savor = orgsv
        call quakes
c       call output(0)
c       call output(1)
        iprn = iprnsv
        inst = instsv
        return
      else
        kgoon = kgoon + 1
c the first time an event reaches here, kgoon will equal 1.
c the 2nd time kgoon will equal 2, so the event will not be rurun again.
        if((inst .ne. 9) .and. (inst .ne. 8)) return
        if(kgoon .ge. 2) then
c do the same as above, starting with statement 130
c         goto 130
          iprnsv = iprn
          iprn = -2
          inst = 9
          freor = .false.
          savla = latsv
          savlo = lonsv
          savez = zsv
          savor = orgsv
          call quakes
c         call output(0)
c         call output(1)
          iprn = iprnsv
          inst = instsv
          return
        endif
        call back(t6*aala(icmin),t6*aalo(icmin),savla,savlo,latsv,lonsv)
        savez = zsv + aalz(icmin)*dez
        if(z .lt. 0.0) z = 0.0
        write(punt,1000)
 1000   format(///' *** run again starting at best nearby point ***',/)
        call quakes
        goto 17
      endif
      end
c end boxau
c cosh.for    []
      function cosh(x)
      cosh = (exp(x) + exp(-x))/2.
      return
      end
c end cosh
c coshi.for    []
      function coshi(x)
      p = x + sqrt(x*x - 1.)
      coshi = alog(p)
      return
      end
c end coshi
c critic.for    []
      subroutine critic(delta,az,w,iuse,nrp,nr,ldx)
c select critical stations
c   1) closest 4 with p-phase readings.
c   2) additional with p- or s-phase readings if they reduce the a gap > 72 deg
c      by 5 deg or more.
c   3) s is used at critical stations.  if there are no s phases selected,
c      then s is used from the closest non-critical station.
      include 'params.inc'
      logical swtchp, swtchs
      dimension dtemp(npa),key(npa)
      dimension delta(npa),az(npa),w(npa),iuse(npa),ldx(npa)
      integer punt
      common /punt/ punt
      write(punt,20)
20    format(' determining which stations are critical.')
      do 40 j = 1, nr
        iuse(j) = 0
        if (w(j) .eq. 0.) iuse(j) = 1
40    continue
      do 50 j = 1, nrp
        dtemp(j) = delta(j)
50    continue
      call sort(dtemp,key,nrp)
      ns = 0
      nwt = 0
c always use first 4 with p-wt .gt. 0.0
      do 100 i = 1, nrp
        j = key(i)
c check for p
        if (w(j) .gt. 0.0) then
          nwt = nwt + 1
          iuse(j) = 1
c add s only if p is selected
          if (ldx(j) .ne. 0) then
            k = ldx(j)
            if (w(k) .gt. 0.0) then
              iuse(k)  = 1
              ns = ns + 1
            endif
          endif
        endif
        if (nwt .ge. 5) goto 125
100   continue
c fall through if there are 4 or fewer p phases that have wt <> 0
125   if (nwt .eq. 0) then
        sge72 = 400.
        nge72 = 0.
      else
        call sumgap(iuse, az, nr, sge72, nge72, w)
      endif
c try one at a time for stations that reduce gap.
      do 200 i = 1, nrp
        j = key(i)
        k = ldx(j)
        if (iuse(j) .eq. 1) then
c p already in use
          if (k .gt. 0) then
            if (iuse(k) .eq. 1) then
c   s also in use, so skip on to next phase
              goto 200
            endif
          else
c   no s
            goto 200
          endif
        endif
c either p and/or s are not now being used - check weights
        if (w(j) .eq. 0.) then
          if (k .gt. 0) then
c there is an s phase
            if (w(k) .eq. 0.) then
c p and s weights are 0, so skip to next phase
              goto 200
            endif
          endif
        endif
c either p and/or s has non zero weight and is not yet being used
c try adding the p or s to test gap reduction
        swtchp = .false.
        swtchs = .false.
        if ((iuse(j) .eq. 0) .and. (w(j) .gt. 0.)) then
          swtchp = .true.
          iuse(j) = 1
        else
          swtchs = .true.
          iuse(k) = 1
        endif
        call sumgap(iuse, az, nr, sgtr, ngtr, w)
c check if the number of gaps larger than 72 has been increased
        if (ngtr .gt. nge72) then
          nge72 = ngtr
c check if the reduction in sum of gaps > 72 has been reduced by 5 or more deg
        else if (sgtr .le. (sge72 - 5.)) then
          nge72 = ngtr
          sge72 = sgtr
        else
c there was no or not enough improvement
          if(swtchp) iuse(j) = 0
          if(swtchs) iuse(k) = 0
          goto 200
        endif
c this phase reduced the gap, so use it!
        if (swtchp) then
          nwt = nwt + 1
c add s if available
          if (k .ne. 0) then
            if (w(k) .gt. 0.0) then
              iuse(k)  = 1
              ns = ns + 1
            endif
          endif
        endif
        if (swtchs) then
          ns = ns + 1
        endif
200   continue
400   if (ns .gt. 1) goto 500
c use closest s if none found at any critical station
      do 350 i = 1,nrp
        j = key(i)
        if (ldx(j) .eq. 0) goto 350
        k = ldx(j)
        if (w(k) .le. 0.) goto 350
        iuse(k) = 1
        iuse(i) = 1
        goto 500
350   continue
500   return
      end
c end critic
c cvrtop.for    []
      subroutine cvrtop(x, y, r, theta)
c
c-------- bruce julian
c
c-------- cvrtop - convert from rectangular to polar coordinates
c
c (output - may overlay x, y)
c
c-------- standard fortran funct. required:  atan2
c-------- funct. required:  hypot
c
      r = hypot(x, y)
      theta = 0.
      if ((y .ne. 0.) .or. (x .ne. 0.)) theta = atan2(y, x)
      return
      end
c
c end cvrtop
c cyldly.for    []
      subroutine cyldly
     *  (kno, alat, alon, lat1, lon1, x0, y0, z, dly, sdly, iprn,
     *   modset, test8)
      save cyldy, cylrd, cylrd1, cylup, cylup1, cyldn, cyldn1,
     *  xc, yc, ncyl, setcyl
      real lat1, lon1
      character*1 ins, iew
      include 'params.inc'
      parameter (ndly = 11)
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      parameter (mxcyl = 50)
      logical setcyl
      real     cyld
c              cyld		horizontal distance to inner cylinder
      real     cyldin(mxcyl)
c              cyldin(i)    	distance to cylinder center for eq within
c                               inner radius
      integer  cyldy(mxcyl)
c              cyldy(i)		delay assigned to cylinder i
      integer  cylmd(mxcyl)
c              cylmd(i)		velocity model assigned to cylinder i
      real     cylrd(mxcyl)
c              cylrd(i)         inner radius of cylinder number i
      real     cylrd1(mxcyl)
c              cylrd1(i)        outer radius of cylinder number i
      integer  cyldm(mxcyl)
c              cyldm(i)		delay model for transition table entry i
      integer  cyldmin(mxcyl)
c              cyldmin(i)	delay model for inner table entry i
      real     cylup(mxcyl)
c              cylup(i)         upper limit of inner cylinder i
      real     cylup1(mxcyl)
c              cylup1(i)        upper limit of outer cylinder i
      real     cyldn(mxcyl)
c              cyldn(i)         lower limit of inner cylinder i
      real     cyldn1(mxcyl)
c              cyldn1(i)        lower limit of outer cylinder i
      real     cylwt(mxcyl)
c              cylwt(i)         weight derived for transition zone i
      real     dly(ndly,nsn)
c              dly(i, j)         p-delay for model i, station j
      integer  ntrans
c              ntrans		number of entries in transition table
      real     sdly(ndly,nsn)
c              sdly(i, j)        s-delay for model i, station j
      real     xc(mxcyl)
c              xc(i)            x coordinate of center of cylinder i
      real     yc(mxcyl)
c              yc(i)            y coordinate of center of cylinder i
      real     z
c              z		current trial eq depth
      real     zmarg
c              zmarg		vertical distance from event to endcap
c                               (negative for lower endcap)
      data setcyl /.false./
      modset = 0
      if(iprn .ge. 5) then
        print *, ' '
        print *, 'begin sub. cyldly to select delay and velocity'
	print *, 'model(s) based on cylindrical domains.'
      endif
      if (.not. setcyl) then
cd       print *, 'read in cylinder definitions with cylget'
        call cylget
     *  ( mxcyl, cyldy, cylmd, cylrd, cylrd1, cylup, cylup1, cyldn,
     *    cyldn1, xc, yc, ncyl, lat1, lon1, test8)
        setcyl = .true.
      endif
c compute the x,y coordinates of the current trial epicenter
      call delaz(lat1, lon1, delt, deldeg, azz, alat, alon)
      x0 = delt*sin(azz*rad)
      y0 = delt*cos(azz*rad)
      if(iprn .ge. 5) then
	call unfold2(alat, alon, la, ins, ala, lo, iew, alo)
	print *, 'current epicenter = ', la, ins, ala, lo, iew, alo
        print *, 'location of epicenter wrt reference station is'
        print *, 'azimuth (deg) = ', azz
        print *, 'distance (km) = ', delt
        print *, 'x, y = ', x0, y0
        print *, 'loop through ', ncyl, ' regions. regions must be in '
        print *, 'order of preference in cases of overlap.'
      endif
      ntrans = 0
      ninner = 0
      do 20 i = 1, ncyl
        cyld = sqrt((x0 - xc(i))**2 + (y0 - yc(i))**2)
         if(iprn .ge. 5) then
	   print *, 'x, y of cylinder ', i, ' is ', xc(i), yc(i)
           print *, 'dist to this cylinder ',
     *     'which uses delay model ', cyldy(i), ' is ', cyld
           print *, 'z = ', z, ' cylup,dn(i) = ', cylup(i), cyldn(i)
        endif
	if ((cyld .le. cylrd(i)) .and. (z .ge. cylup(i)) .and.
     *    (z .le. cyldn(i))) then
          if(iprn .ge. 5)
     *      print *, 'location is within inner cylinder'
	  ninner = ninner + 1
	  cyldin(ninner) = cyld
	  cyldmin(ninner) = i
        else if ((cyld .le. cylrd1(i)) .and. (z .ge. cylup1(1)) .and.
     *    (z .le. cyldn1(i))) then
c
cd         print *, 'location is within transition zone, so compute '
cd         print *, 'weight and add to list'
          ntrans = ntrans + 1
          cyldm(ntrans) = cyldy(i)
cd         print *, 'zmarg is distance into upper or lower cap'
          zmarg = 0.0
          if ((z .lt. cylup(i)) .and. (z .gt. cylup1(i)))
     *    zmarg = cylup(i) - z
          if ((z .gt. cyldn(i)) .and. (z .lt. cyldn1(i)))
     *    zmarg = cyldn(i) - z
cd         print *, 'note that zmarg will be negative for lower ',
cd    *             'cap region'
          if (zmarg .eq. 0.0) then
c
cd           print *, 'adjacent to cylinder '
cd           print *, 'cyld - cylrd(i) ', cyld - cylrd(i)
cd           print *, 'pi ', pi
cd           print *, 'cylrd1(i) - cylrd(i) ', cylrd1(i) - cylrd(i)
            cylwt(ntrans) = 0.5 + 0.5*cos( pi*(cyld - cylrd(i)) /
     *      (cylrd1(i) - cylrd(i)) )
          else if (cyld .le. cylrd(i)) then
c
cd           print *, 'within end caps'
            if (zmarg .gt. 0) then
cd             print *, 'within upper cap'
              cylwt(ntrans) = 0.5 + 0.5*cos( pi*zmarg /
     *        (cylup(i) - cylup1(i)) )
            else
c
cd             print *, 'within lower cap'
              cylwt(ntrans) = 0.5 + 0.5*cos( pi*zmarg /
     *        (cyldn(i) - cyldn1(i)) )
            endif
          else
c
cd           print *, 'within corner zone'
            if (zmarg .gt. 0) then
cd             print *, 'within upper corner zone'
              cylwt(ntrans) = 0.5 + 0.5*cos(pi*sqrt
     *        ( ((cyld - cylrd(i)) /
     *          (cylrd1(i) - cylrd(i)))**2 +
     *          (zmarg /
     *          (cylup(i) - cylup1(i)))**2 ) )
            else
cd             print *, 'within lower corner zone'
              cylwt(ntrans) = 0.5 + 0.5*cos(pi*sqrt
     *        ( ((cyld - cylrd(i)) /
     *          (cylrd1(i) - cylrd(i)))**2 +
     *          (zmarg /
     *          (cyldn(i) - cyldn1(i)))**2 ) )
            endif
          endif
        endif
20    continue

c use the parameters for the cylinder for which the eq is
c closest to the center
      if (ninner .ne. 0) then
	if(ninner .eq. 1) then
          kno = cyldy(cyldmin(ninner))
 	  modset = cylmd(cyldmin(ninner))
          return
        else
	  ntouse = 1
	  do 22 i = 2, ninner
	    if(cyldin(i) .lt. cyldin(ntouse)) then
	      ntouse = i
	    endif
22        continue
          kno = cyldy(cyldmin(ntouse))
 	  modset = cylmd(cyldmin(ntouse))
          return
	endif
      endif
	
	
      if (ntrans .eq. 0) then
cd       print *, 'not within any cylinders, so use default model (1)'
        kno = 1
        return
      else
cd       print *, 'weights: ', (cylwt(i), i = 1, ntrans)
cd       print *, 'models:  ', (cyldm(i), i = 1, ntrans)
c reduce list if there are 2 or more entries
        if (ntrans .gt. 1) call cylred(cyldm, cylwt, ntrans, sumwt)
        kno = ndly
c
        if (ntrans .eq. 1) then
          if (cylwt(1) .ge. 1.0) then
            kno = cyldm(1)
            return
          else
c           fill in with default model
            ntrans = 2
            cylwt(2) = 1.0 - cylwt(1)
            cyldm(2) = 1
            if (iprn .ge. 5) then
              print *, 'weights: ', (cylwt(i), i = 1, ntrans)
              print *, 'models:  ', (cyldm(i), i = 1, ntrans)
            endif
          endif
c
        else if (ntrans .eq. 3) then
          cylwt(1) = cylwt(1)/sumwt
          cylwt(2) = cylwt(2)/sumwt
          cylwt(3) = cylwt(3)/sumwt
          sumwt = 1.0
cd         print *, 'weights: ', (cylwt(i), i = 1, ntrans)
cd         print *, 'models:  ', (cyldm(i), i = 1, ntrans)
        else
c
c         ntrans = 2
          if (sumwt .gt. 1.0) then
            cylwt(1) = cylwt(1)/sumwt
            cylwt(2) = cylwt(2)/sumwt
          else
c           fill in with default model
            ntrans = 3
            cylwt(3) = 1. - sumwt
            sumwt = 1.0
            cyldm(3) = 1
cd           print *, 'weights: ', (cylwt(i), i = 1, ntrans)
cd           print *, 'models:  ', (cyldm(i), i = 1, ntrans)
          endif
        endif
      endif
c compute combined delays
      do 30 i = 1, nsn
        dly(ndly,i) = 0.0
        sdly(ndly,i) = 0.0
        do 28 j = 1, ntrans
          dly(ndly,i) =   dly(ndly,i) +  dly(cyldm(j), i)*cylwt(j)
          sdly(ndly,i) = sdly(ndly,i) + sdly(cyldm(j), i)*cylwt(j)
28      continue
30    continue
      return
      end
c end cyldly
c cylget.for    []
      subroutine cylget
     * ( mxcyl, cyldy, cylmd, cylrd, cylrd1, cylup, cylup1, cyldn,
     *   cyldn1, xc, yc, ncyl, lat1, lon1, test8 )
      real lat1, lon1
      parameter (ndly = 11)
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      character*1 ins, iew
      character*110 record
c read in the cylinders defining weight model regions
      integer  cyldy(mxcyl)
c              cyldy(i)		delay assigned to cylinder i
      integer  cylmd(mxcyl)
c              cylmd(i)         velocity model assigned to cylinder i
      real     cylrd(mxcyl)
c              cylrd(i)         inner radius of cylinder number i
      real     cylrd1(mxcyl)
c              cylrd1(i)        outer radius of cylinder number i
      real     cylup(mxcyl)
c              cylup(i)         upper limit of inner cylinder i
      real     cylup1(mxcyl)
c              cylup1(i)        upper limit of outer cylinder i
      real     cyldn(mxcyl)
c              cyldn(i)         lower limit of inner cylinder i
      real     cyldn1(mxcyl)
c              cyldn1(i)        lower limit of outer cylinder i
      real     xc(mxcyl)
c              xc(i)            x coordinate of center of cylinder i
      real     yc(mxcyl)
c              yc(i)            y coordinate of center of cylinder i
      ncyl = 0
cd     print *, 'read cylinder parameters'
      do 20 i = 1, mxcyl
18      read(17, '(a)', end = 90) record
        if(record(1:2) .eq. 'c*') goto 18
c because standard fortran 77 to can not read an internal file with
c free format, file 14 must be used in the following code!
        rewind 14
        write(14, '(a)') record
        rewind 14
c       read(record, *)
        read(14, *)
     *  cyldy(i),
     *  cylmd(i), blat, blon, cylrd(i), cylrd1(i), cylup(i),
     *  cylup1(i), cyldn(i), cyldn1(i)

c       adjust for the depth of the top of the model
        cylup(i) = cylup(i) + test8
        cylup1(i) = cylup1(i) + test8
        cyldn(i) = cyldn(i) + test8
        cyldn1(i) = cyldn1(i) + test8

c       compute the x,y coordinates of the cylinder center
        nlat = abs(blat)
        rlat = (abs(blat) - nlat)*60.
        ins = 'n'
        if (blat .lt. 0) ins = 's'
        nlon = abs(blon)
        rlon = (abs(blon) - nlon)*60.
        iew = 'w'
        if (blon .lt. 0) iew = 'e'
c convert to geocentric coordinates> alat, alon
cd       print *, 'blat, blon = ', blat, blon
cd       print *, 'nlat, ins, rlat, nlon, iew, rlon'
cd       print *, nlat, ins, rlat, nlon, iew, rlon
        call fold2(alat, alon, nlat, ins, rlat, nlon, iew, rlon)
        call delaz(lat1, lon1, delt, deldeg, azz, alat, alon)
cd       print *, 'location of cylinder wrt reference station is'
cd       print *, 'cyl lat lon = ', alat, alon
cd 	 print *, 'ref sta lat lon = ', lat1, lon1
cd       print *, 'azimuth (deg) = ', azz
cd       print *, 'distance (km) = ', delt
        xc(i) = delt*sin(azz*rad)
        yc(i) = delt*cos(azz*rad)
cd       print *, 'converted to x, y = ', xc(i), yc(i)
cd       print *, 'i, cyldy(i), cylmd(i), cylrd(i), cylrd1(i),',
cd    *  'cylup(i), cylup1(i), cyldn(i), cyldn1(i)'
cd       print *,  i, cyldy(i), cylmd(i), cylrd(i), cylrd1(i),
cd    *  cylup(i), cylup1(i), cyldn(i), cyldn1(i)
        if (cyldy(i) .gt. ndly) then
          print *, 'sub. cylget attempting to assign delay model',
     *    cyldy(i)
          print *, 'but only ', ndly, ' are allowed.'
          stop
        endif
20    continue
      ncyl = mxcyl
      return
90    ncyl = i - 1
      close (17)
      return
      end
c end cylget
c cylred.for    []
      subroutine cylred(cyldm, cylwt, ntrans, sumwt)
c reduce delay model list, eliminating duplicates and sorting by weight
      integer cyldm(ntrans)
c             cyldm(i)            delay model for # i
      real    cylwt(ntrans)
c             cylwt(i)            weight for # i
      integer key(50)
c             key                 key to original order of weights
      real    sumwt
c             sumwt               sum of weights of up to first 3
      integer tmpdm(50)
c             tmpdm               temp array of delay models
      real    tmpwt(50)
c             tmpwt               temp array of weights
      ntr = 0
      do 20 i = 1, ntrans
        if(cylwt(i) .gt. 0.0) then
          ntr = ntr + 1
          if(i .lt. ntrans) then
            do 18 j = i+1, ntrans
              if((cyldm(j) .eq. cyldm(i)) .and. (cylwt(j) .ne. 0.)) then
                cylwt(i) = cylwt(i) + cylwt(j)
                cylwt(j) = 0.0
              endif
18          continue
          endif
        endif
20    continue
      ntrans = ntr
c sort the weights
      call sort(cylwt, key, ntrans)
      do 30 i = 1, ntrans
        tmpdm(i) = cyldm(i)
30    continue
      do 40 i = 1, ntrans
        cyldm(ntrans + 1 - i) = tmpdm(key(i))
40    continue
c save weights in temporary arrays
      do 50 i = 1, ntrans
        tmpwt(i) = cylwt(i)
50    continue
c change to ascending order
      if (ntrans .gt. 3) ntrans = 3
      sumwt = 0.0
      do 60 i = 1, ntrans
        cylwt(i) = tmpwt(ntrans + 1 - i)
        sumwt = sumwt + cylwt(i)
60    continue
      return
      end
c end cylred
c delaz.for    []
      subroutine delaz(eqlat, eqlon, dekm, dedeg, az0, slat, slon)
c
c-------- delaz - calculate the distance in km (approx equal to geocentric
c      distance times local radius), and azimuths in radians
c
c input:  slat     station geocentric latitude in radians
c         slon     station longitude in radians
c         eqlat    earthquake geocentric latitude in radians
c         eqlon    earthquake longitude in radians
c output: dekm     distance from earthquake to station in kilometers
c         dedeg    distance from earthqauke to station in degrees
c         az0      azimuth from earthquake to station measured clockwise
c                     from north in degrees
c
      parameter (pi = 3.14159265)
      parameter (twopi = 2.0*pi)
      parameter (halfpi = 0.5*pi)
      parameter (rad = pi/180.)
      parameter (deg = 1.0/rad)
      parameter (equrad = 6378.2064)
      parameter (polrad = 6356.5838)
      parameter (flat = (equrad - polrad)/equrad)
c
      st0 = cos(eqlat)
      ct0 = sin(eqlat)
c use approximation of local radius for derivative of surface
c     distance with geocentric latitude.
c  more accurate formulation would be:
c      drdth = -r**3 * cos(slat)*sin(slat)*( (1.-flat)**(-2) - 1. )/a**2
c      dsd = sqrt(drdth**2 + r**2)
      radius = (cos(eqlat)**2/equrad**2 +
     *sin(eqlat)**2/polrad**2)**(-.5)
      ct1 = sin(slat)
      st1 = cos(slat)
      sdlon = sin(eqlon-slon)
      cdlon = cos(eqlon-slon)
      cdelt = st0*st1*cdlon+ct0*ct1
      call cvrtop(st0*ct1-st1*ct0*cdlon, st1*sdlon, sdelt, az0)
      dedeg = atan2(sdelt, cdelt)*deg
      dekm = radius*atan2(sdelt, cdelt)
      if (az0 .lt. 0.0) az0 = az0 + twopi
      if (az0 .ge. twopi) az0 = az0 - twopi
      az0 = az0*deg
c calculation of back azimuth if needed
c      call cvrtop(st1*ct0 - st0*ct1*cdlon, (-sdlon)*st0, sdelt, az1)
c      if (az1 .lt. 0.0) az1 = az1 + twopi
      return
      end
c end delaz
c diftim.for    []
      subroutine diftim(icent1, iymd1,ihr1,icymd2,ihr2,iporm,difhr)
c calculate first time minus second time
c iporm  -1  first cnyrmody is smaller than second, set difhr = -1.0
c             0  first cnyrmody is equal to second, calculate difhr
c            +1  first cnyrmody is larger, set difhr = +1.0
c difhr       this is the time difference in hours, if dates are =
c
      if (icent1*1000000+iymd1 .lt. icymd2) goto 100
      if (icent1*1000000+iymd1 .gt. icymd2) goto 200
c equal dates
      iporm = 0
      difhr = ihr1 - ihr2
      return
c smaller first date
  100 iporm = -1
      difhr = -1.
      return
c larger first date
  200 iporm = +1
      difhr = +1.
      return
      end
c end diftim
c dnstrg.for    []
      character*(*) function dnstrg(array)
c
c  program to change string to uppercase
c
c  array - a character variable
      character*(*) array
      integer offset
      data offset/32/
c
c  get length of array
c
      dnstrg = ' '
      lenstr = len(array)
      if (lenstr .eq. 0) return
      do 10 i = 1, lenstr
        ic = ichar(array(i:i))
        if ((ic .ge. 65) .and. (ic .le. 90)) then
          dnstrg(i:i) = char(ic + offset)
        else
          dnstrg(i:i) = array(i:i)
        endif
   10 continue
      return
      end
c end dnstrg
c dubl.for    [pc]
      double precision function dubl(x)
c     the wadati sub uses p and s arrival times that are entered
c     to the nearest .01 sec in single precision.  the purpose of this
c     sub is to convert these numbers to the nearest double
c     precision number, setting less significant decmial digits to zero.
c     this version is 16 times faster than dubl1, which uses internal
c     write and read statements.
c     lahr & stephens   feb 1986
c
c* (vax
c      double precision dfact
c      if (x .eq. 0.) then
c        dubl = 0.d0
c        return
c      endif
c      al = alog10(abs(x))
c      if (al .lt. 0.) al = al - .9999995
c      i = al
c      dfact = 10.d0**(7-i)
c      ix = x*dfact + sign(.5, x)
c      dubl = dflotj(ix)/dfact
c* vax)
c* (pc
c* (unix
      dubl = x
c* unix)
c* pc)
      return
      end
c end dubl
c dwnwrd.for    []
      subroutine dwnwrd(alrms, altla, altlo, altz,
     * rmslim, zdn, axz)
      logical freor
      common /hopq/ savla,savlo,savez,savor,freor
      common /phoqn/ inst,knst
      real lonep, latep
      common /qmost1/ lonep,ni,latep
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      dimension zinc(8)
      data zinc/20, 20, 20, 20, 20, 20, 20, 20/
      inst = 1
c set initial steps according to erz estimate
      zbas = axz
      if(zbas .lt. .2) zbas = .2
      if(zbas .gt. 20.) zbas = 20.
      zinc(1) = zbas
      zinc(2) = zbas*2.
c find downward shift in z with respect to altz that results in rms = rmslim
      savla = altla
      savlo = altlo
      savez = altz + zinc(1)
      alrmsl = alrms
      altzl = altz
      n = 1
38    call quakes
      if(rms .ge. rmslim) then
        zdn = (savez - altzl)*(rmslim - alrmsl) /
     *          (rms - alrmsl) + altzl - altz
      else
        n = n + 1
        altzl = savez
        savez = savez + zinc(n)
        if(savez - altz .gt. 110.) then
          zdn = 99.
          return
        endif
        alrmsl = rms
        savla = latep
        savlo = lonep
        goto 38
      endif
      return
      end
c end dwnwrd
c eigen1.for    []
      subroutine eigen1(a, mev, mv, n, ev, s, v, damp)
c used to find eigenvalues and eigenvectors for the upper left nxn
c     portion of v(mv, mv)
c
c modified from an ibm sub. by j. c. lahr
c
c     dimension of r must be greater than n**2,
c     so n may not exceed 10 in this version
c              mev              (input) actual dimensions of ev
      integer  mev
c              mv               (input) actual dimensions of v
      integer  mv
c              n                (input) tensor dimension
      integer  n
c              a(n + (n*n-n)/2) input v transformed to storage mode 1
      real     a(n + (n*n-n)/2)
c              damp             (input) term to be added to diagonal elements
      real     damp
c              ev(mev, mev)     (output) eigenvectors
      real     ev(mev, mev)
c              r(100)           eigenvectors as linear array
      real     r(100)
c              s(n)             (output) eignevalues
      real     s(n)
c              v(mv, mv)        (input) symetric tensor to be diagonalized
      real     v(mv, mv)
c
      do 10 j = 1, n
      do 10 i = 1, j
        ni = i + (j*j-j)/2
        a(ni) = v(i, j)
        if (i .eq. j) a(ni) = a(ni) + damp
10    continue
      iq = -n
      do 20 j = 1, n
        iq = iq+n
        do 20 i = 1, n
          ij = iq+i
          r(ij) = 0.0
          if (i-j) 20, 15, 20
15        r(ij) = 1.0
20    continue
      anorm = 0.0
      do 35 i = 1, n
      do 35 j = i, n
        if (i-j) 30, 35, 30
30      ia = i+(j*j-j)/2
        anorm = anorm+a(ia)*a(ia)
35    continue
      if (anorm) 165, 165, 40
40    anorm = 1.414*sqrt(anorm)
      anrmx = anorm*1.0e-6/float(n)
      ind = 0
      thr = anorm
45    thr = thr/float(n)
50    l = 1
55    m = l+1
60    mq = (m*m-m)/2
      lq = (l*l-l)/2
      lm = l+mq
      if ( abs(a(lm))-thr) 130, 65, 65
65    ind = 1
      ll = l+lq
      mm = m+mq
      x = 0.5*(a(ll)-a(mm))
      y = -a(lm)/ sqrt(a(lm)*a(lm)+x*x)
      if (x) 70, 75, 75
70    y = -y
75    sinx = y/sqrt(2.0*(1.0+(sqrt((1.-y)*(1.+y)))))
      sinx2 = sinx*sinx
      cosx = sqrt((1.-sinx)*(1.+sinx))
      cosx2 = cosx*cosx
      sincs  = sinx*cosx
      ilq = n*(l-1)
      imq = n*(m-1)
      do 125 i = 1, n
        iq = (i*i-i)/2
        if (i-l) 80, 120, 80
80      if (i-m) 85, 120, 90
85      im = i+mq
        goto 95
90      im = m+iq
95      if (i-l) 100, 105, 105
100     il = i+lq
        goto 110
105     il = l+iq
110     x = a(il)*cosx-a(im)*sinx
        a(im) = a(il)*sinx+a(im)*cosx
        a(il) = x
120     ilr = ilq+i
        imr = imq+i
        x = r(ilr)*cosx-r(imr)*sinx
        r(imr) = r(ilr)*sinx+r(imr)*cosx
        r(ilr) = x
125   continue
      x = 2.0*a(lm)*sincs
      y = a(ll)*cosx2+a(mm)*sinx2-x
      x = a(ll)*sinx2+a(mm)*cosx2+x
      a(lm) = (a(ll)-a(mm))*sincs+a(lm)*(cosx-sinx)*(cosx+sinx)
      a(ll) = y
      a(mm) = x
130   if (m-n) 135, 140, 135
135   m = m+1
      goto 60
140   if (l-(n-1)) 145, 150, 145
145   l = l+1
      goto 55
150   if (ind-1) 160, 155, 160
155   ind = 0
      goto 50
160   if (thr-anrmx) 165, 165, 45
165   iq = -n
      do 185 i = 1, n
        iq = iq+n
        ll = i+(i*i-i)/2
        jq = n*(i-2)
        do 185 j = i, n
          jq = jq+n
          mm = j+(j*j-j)/2
          if (a(ll)-a(mm)) 170, 185, 185
170       x = a(ll)
          a(ll) = a(mm)
          a(mm) = x
          do 180 k = 1, n
            ilr = iq+k
            imr = jq+k
            x = r(ilr)
            r(ilr) = r(imr)
180       r(imr) = x
185   continue
      do 190 j = 1, n
        ni = j + (j*j-j)/2
        s(j) = a(ni)
190   continue
      ni = 0
      do 200 i = 1, n
        do 200 j = 1, n
          ni = ni + 1
          ev(j, i) = r(ni)
200   continue
      return
      end
c end eigen1
c erset.for    [pc]
      subroutine erset (p, a, b, c, d)
      integer p
      logical a, b, c, d
c* (vax
c      call errset (p, a, b, c, d)
c* vax)
c* (pc
c* pc)
c* (unix
c* unix)
      return
      end
c end erset
c fmplot.for    []
      subroutine fmplot
c plot first-motion directions on the lower focal hemisphere
c------- in equal area projection
      include 'params.inc'
      logical repeat
      real bat2,bon2,mag
      character*4 iahead*60, msta*5, nsta*5, icard*110
      integer punt
      common /punt/ punt
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /dmost/ ipun,ivlr,blank
      common /gmost/ az(npa)
      common /hf/ cp(72),sp(72)
      common /imost/ test(100)
      common /obcfn/ ain(npa)
      common /ocfn/ adj,seh,lat1,bat2,lon1,bon2,igap
      character*1 isnla, isnlo, krm1, krm2
      common /ocfn1/ isnla, isnlo, krm1, krm2
      common /ofln/ sec
      common /ofnv/ kmin
      common /omnfh/ dmin,dmin3,sminp
      character msym*1
      common /pfo/ msym(npa)
      common /pfnoqv/ kdate,khrmn,khr
      common /pno/ lph,keyph(npa),nsum
      character phcard*117, krm*2, seq*5
      common /pno1/ phcard(npa), krm(npa), seq
      common /qmost/ wt(npa),z
      common /xfmno/ mag
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
c
c     for multics printer use xscale = 0.101064
c     for slac printer use xscale = 0.1379
      xscale = test(45)
      repeat = .false.
      if (test(7) .lt. 0.) repeat = .true.
      write(punt,2)
    2 format(52h1  date    origin    lat n    long w    depth    mag,
     * 17h no gap dmin  rms)
      rmag = mag
      if (mag .eq. blank) rmag = 0.0
      write(punt,5) kdate,khr,kmin,sec,lat1,bat2,lon1,bon2,krm1,
     * z-test(8),krm2,rmag,nrwt,igap,dmin,rms
    5 format(2x,i6.6,1x,2i2,f6.2,i3,1h-,f5.2,i4,1h-,f5.2,a1,f6.2,a1
     *,f6.2,i3,i4,f5.1,f5.2)
      call fmplt(msta, az, ain, msym, lph, xscale, repeat, punt,
     *                 cp, sp)
      return
      end
c end fmplot
c fmplt.for    []
      subroutine fmplt(msta, az, ain, fm, lph, xscale, repeat, punt,
     *                 cp, sp)
c plot first-motion directions on the lower focal hemisphere
c------- in equal area projection
      include 'params.inc'
      integer punt
      character*1 fm(npa), jtem, igrap(95, 95)
      character*5 msta(npa)
      dimension az(npa), ain(npa)
      dimension cp(72), sp(72)
      logical repeat
      data nox, noy, iy, noy1, noy2/95, 59, 24, 57, 30/
      data rmax, yscale, add/3.937008, 0.169643, 4.75/
      ntin = 0
    6 ntin = ntin + 1
      nfmr = 0
c zero graph
      do 10 i = 1, nox
      do 10 j = 1, noy
   10 igrap(i, j) = ' '
c make circle of *'s
      do 20 i = 1,72
        jx = (rmax*cp(i)+add)/xscale + 1.5
        jy = (rmax*sp(i)+add)/yscale + 0.5
        jy = noy - jy - 1
        ii = i - (i/2)*2
        igrap(jx, jy) = '*'
        if (ii .eq. 0) igrap(jx, jy) = ' '
   20 continue
      nox2 = add/xscale + 1.5
      it = (-rmax + add)/xscale + 1.5
      igrap(it, noy2) = '-'
      it = (rmax + add)/xscale + 2.5
      igrap(it, noy2) = '-'
      it = noy2 - iy - 1
      igrap(nox2, it) = 'i'
      it = noy2 + iy + 1
      igrap(nox2, it) = 'i'
c     center is nox2, noy2
      igrap(nox2, noy2) = '*'
      do 50 i = 1, lph
        if (fm(i) .eq. ' ') goto 50
        if (ain(i) .gt. 90.) goto 30
        ann = ain(i)
        azz = az(i)*.0174533
        goto 32
   30   ann = 180. - ain(i)
        azz = (180. + az(i) )*.0174533
   32   r = rmax*1.414214*sin(ann*.0087266)
        x = r*sin(azz) + add
        y = r*cos(azz) + add
        jx = x/xscale + 1.5
        jy = y/yscale + .5
        jy = noy - jy - 1
        if (ntin .eq. 2) goto 52
        jtem = igrap(jx, jy)
c if previous symbol is weak, overwrite it.
        if ((jtem.eq.' ').or.(jtem.eq.'*').or.(jtem.eq.'+')
     *  .or.(jtem.eq.'-').or.(jtem.eq.' ').or.(jtem.eq.'?')) goto 47
c if new symbol is weak, do not overwrite strong symbol.
        if ((fm(i).eq.'+').or.(fm(i).eq.'-').or.(fm(i).eq.'?'))
     *  goto 50
        if (fm(i) .eq. 'c' ) goto 40
c plot d on top of previous strong symbol
        if (igrap(jx, jy) .ne. 'd') goto 35
        igrap(jx, jy) = 'e'
        goto 50
   35   if (igrap(jx, jy) .ne. 'e') goto 37
        igrap(jx, jy) = 'f'
        goto 50
   37   if (igrap(jx, jy) .eq. 'f') goto 50
        igrap(jx, jy) = 'x'
        goto 50
c plot c on top of privious strong symbol
   40   if (igrap(jx, jy) .ne. 'c') goto 43
        igrap(jx, jy) = 'b'
        goto 50
   43   if (igrap(jx, jy) .ne. 'b') goto 45
        igrap(jx, jy) = 'a'
        goto 50
   45   if (igrap(jx,jy) .eq. 'a') goto 50
        igrap(jx, jy) = 'x'
        goto 50
   47   igrap(jx, jy) = fm(i)
        goto 50
c write station name on focal sphere
   52   njx = jx - 1
        do 56 j = 1, 4
          if (msta(i)(j:j) .ne. ' ') igrap(njx, jy) = msta(i)(j:j)
          njx = njx + 1
   56   continue
   50 continue
      nm1 = noy1 - 1
      do 80 i = 4, nm1
        write(punt, 55) (igrap(j, i), j = 1, nox)
   55   format(21x, 95a1)
   80 continue
      if ((repeat) .and. (ntin .eq. 1)) then
        write(punt, 82)
   82   format(1h1, /)
        goto 6
      endif
      return
      end
c end fmplt
c fold2.for    []
      subroutine fold2(alat,alon,la,ins,ala,lo,iew,alo)
c
c-------- given geographic coordinates compute geocentric lat and lon
c
c input:  la     degree portion of latitude in degrees
c         ins    n for north, s for south
c         ala    minutes portion of latitude
c         lo     degree portion of longitude
c         iew    e for east, w for west
c         alo    minutes portion of longitude
c output: alat   geocentric latitude in radians
c         alon   longitude in radians
      parameter (pi = 3.14159265)
      parameter (twopi = 2.0*pi)
      parameter (halfpi = 0.5*pi)
      parameter (rad = pi/180.)
      parameter (deg = 1.0/rad)
      parameter (equrad = 6378.2064)
      parameter (polrad = 6356.5838)
      parameter (flat = (equrad - polrad)/equrad)
      parameter (c1 = (1.0 - flat)**2)
      parameter (c2 = halfpi*(1.0/c1 - 1.0))
c
      character*1 ins, iew, dnstrg
c
      alat = (la + ala*1.6666667e-2)*rad
c ggtogc - convert from geographic to geocentric latitude
      if (halfpi-abs(alat) .ge. 0.02) goto 201
         alat = alat/c1-sign(c2,alat)
         goto 202
  201    alat = atan(c1*tan(alat))
  202    continue
      if (dnstrg(ins) .eq. 's') alat = -alat
      alon = (lo + alo*1.6666667e-2)*rad
      if (dnstrg(iew) .eq. 'e') alon = -alon
      return
      end
c end fold2
c formal.for    []
      subroutine formal(x, ix, n, nsig, fmit, xout)
c test program for subroutine formal
c      character*6 fmit
c      print *, 'test formal'
c30    n = iaskk ('give number of columns for number', n)
c      print *, 'format f3.1 uses 1 decimal digit.'
c      nsig = iaskk ('give number of decimal digits in read statement',
c     1nsig)
c      x = raskk ('give number to be printed out', x)
c      call formal(x, ix, n, nsig, fmit, xout)
c      if (fmit .eq. ' ') then
c        print *, 'write ', ix, ' with format i',n
c      else
c        print *, 'write ', xout, ' with format ', fmit
c      endif
c      goto 30
c      end
c---- converts real number x into the integer ix to be
c---- written with format(in) and read with format(fn.nsig)
c---- corrections by willy aspinall, principia testing  july 1983
c
c---- if the number is too large to write as an integer and
c---- n - nsig is greater than 1, then find xout and fmit
c---- so that the number may be written as xout with fmit.
c
      character*6 fmit
      ix = x*(10.**nsig) + sign(0.50001,x)
      imax = 10**n
      fmit = '      '
      if(x .gt. 0.) then
c       positive number
        if(ix .ge. imax) then
          if(nsig .lt. 2) then
            ix = imax - 1
            return
          else
c           in this case pass back a format to be used in writing
            call formit(x, xout, fmit, n, 1)
            return
          endif
        endif
        if(ix .le. imax/100) then
          if((n - nsig) .gt. 1) then
c           in this case pass back a format to be used in writing
            call formit(x, xout, fmit, n, 1)
            return
          endif
        endif
      else
c       negative number
        if(ix .le. -imax/10) then
          if(nsig .lt. 2) then
            ix = -imax/10 + 1
            return
          else
c           in this case pass back a format to be used in writing
            call formit(x, xout, fmit, n, 1)
            return
          endif
        endif
        if(ix .ge. -imax/1000) then
          if((n - nsig) .gt. 2) then
c           in this case pass back a format to be used in writing
            call formit(x, xout, fmit, n, 1)
            return
          endif
        endif
      endif
      end
c end formal
c formf.for    []
      subroutine formf(x, ix, n, nsig)
c converts real number x into the integer ix to be
c written with format(in) and read with format(fn.nsig)
c corrections by willy aspinall, principia testing  july 1983
      ix = x*(10.**nsig) + sign(0.50001,x)
      imax = 10**n
      if(ix .ge. imax) ix = imax - 1
      if(ix .le. -imax/10) ix = -imax/10 + 1
      return
      end
c end formf
c formit.for    []
      subroutine formit(a, aout, fmit, n, ifix)
c        squeez a positive real number (a) into a field of
c        length n where 2<n<9, and maintain the maximum number of
c        significant digits.
c     input: n (integer) and a (real).
c            ifix = 0, never change input number
c            ifix .ne. 0, if number is out of range, set equal to
c            maximum or minimum number that can be printed.
c     output: best format (fmit) of form (fn.x) (a6)
c             aout is the number to be printed
c
      character*6 fmit
      character*1 int(10), ilem(4)
      data int/'0','1','2','3','4','5','6','7','8','9'/
c
      ilem(1) = 'f'
      ilem(2) = int(n+1)
      ilem(3) = '.'
      nm1 = n - 1
      nm2 = n - 2
      aout = a

      if(ifix .eq. 0) go to 20

      top = 10.**nm1 - 0.5
      bot = -(10.**nm2) + 0.5
      if(aout .ge. top) aout = top - .1
      if(aout .le. bot) aout = bot + .1

20    if(aout .eq. -0.) then
        i = n - 1
      else if(aout .eq. 0.) then
	i = n
      else if(aout .lt. 0.) then
        do 30 i = 1,nm2
          nm = nm2 - i
          ip1 = i + 1
          if(aout .le. -(10.**nm - 5./10.**ip1)) go to 50
30      continue
        i = nm1
      else
        do 40 i = 1,nm1
          nm = nm1 - i
          ip1 = i + 1
          if(aout .ge. (10.**nm - 5./10.**ip1)) go to 50
40      continue
        i = n
      endif
50    ilem(4) = int(i)

      write(fmit, 70) ilem
70    format('(', 4a1, ')')
      return
      end
c end formit
c geterr.for    []
      subroutine geterr(zup, zdn, rmslim)
c this subroutine will find the error limits for depth.
c j. c. lahr 4/22/87
      include 'params.inc'
      parameter (ndly = 11)
      real bat2,bon2,lat,lon
      logical good, eoff, supout
      common /gg/ altla(3), altlo(3), altz(3), altrms(3), nsolut,
     * fla(20), flo(20), fz(20), frms(20), forg(20), maxf, altorg(3)
      character*4 iahead*60, msta*5, nsta*5, icard*110
      common /char/ iahead, msta(npa), nsta(nsn), icard
      character*4 ipro, ichec, evtype*1, evstat*1, root*256
      common /dbhip/ ipro, ichec, evtype, evstat, root
      common /dmost/ ipun,ivlr,blank
      common /dhil/ iq,ilat,kms,icat
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /dhio/ nedit,good
      common /dio/ rmsmx,presmx,sresmx,noutmx,nimx,iprun,semx
      common /dph/ noswt, eoff
      common /ghnq/ iexit
      common /gmost/ az(npa)
      common /hf/ cp(72),sp(72)
      common /hl/ nwad, tslope, tsqsl
      common /hpn/ wslope
      logical freor
      common /hopq/ savla,savlo,savez,savor,freor
      common /iclmpq/ lat(nsn),lon(nsn)
      common /logfil/ logfil
      common /idno/ ksel,ksort
      common /ilv/ c(nsn), e(nsn)
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /lm/ mapend
      common /ocfn/ adj,seh,lat1,bat2,lon1,bon2,igap
      character*1 isnla, isnlo, krm1, krm2
      common /ocfn1/ isnla, isnlo, krm1, krm2
      common /ohq/ gap, supout
      common /omnfh/ dmin,dmin3,sminp
      common /on/ ix,iy,iz,ax1,axz
      common /ph/ nfirst
      common /pno/ lph,keyph(npa),nsum
      character*4 krms
      common /po/ krms(npa)
      common /pfnoqv/ kdate,khrmn,khr
      common /pgqv/ w(npa)
      common /phoqn/ inst,knst
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /pnoqtx/ ldx(npa)
      character*4 krmp
      common /pnoqv/ krmp(npa)
      common /pot/ ts(npa),ihr,model(npa), keyphi(npa)
      common /qmost/ wt(npa),z
      real latep,lonep
      common /qmost1/ lonep,ni,latep
      common /ro/ yse,seorg,phi
      common /reloc/ irelo, nreloc
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      common /zqr/ fno
c
c given a local minimum, find the depth limits within which the
c rms will not increase by more than expected, given yse, the
c reading standard error.
c
c set iprn = -2 to supress printing
      isvprn = iprn
      iprn = -2
      rmslim = sqrt(altrms(1)**2 + (yse**2)/fno)
c
c first consider the case of only one solution
      if(nsolut .eq. 1) then
c
c
        call upward(altrms(1), altla(1), altlo(1), altz(1),
     *  frms(1), rmslim, zup, axz)
c
        call dwnwrd(altrms(1), altla(1), altlo(1), altz(1),
     *  rmslim, zdn, axz)
c
      else
c next consider the case of two solutions
        if(altz(1) .lt. altz(2)) then
c in this case, the primary solution is more shallow.
c get the upper limit above the more shallow solution.
          call upward(altrms(1), altla(1), altlo(1), altz(1),
     *    frms(1), rmslim, zup, axz)
c
c get the lower limit below the deeper solution.
          call dwnwrd(altrms(2), altla(2), altlo(2), altz(2),
     *    rmslim, zdn, axz)
          zdn = zdn + (altz(2) - altz(1))
        else
c in this case, the secondary solution is more shallow.
c get the upper limit above the more shallow solution.
          call upward(altrms(2), altla(2), altlo(2), altz(2),
     *    frms(1), rmslim, zup, axz)
          zup = zup + (altz(1) - altz(2))
c
c get the lower limit below the deeper solution.
          call dwnwrd(altrms(1), altla(1), altlo(1), altz(1),
     *    rmslim, zdn, axz)
        endif
c
      endif
      iprn = isvprn
      return
      end
c end geterr
c getsta.for    []
      subroutine getsta(   indexs, nsta, ielv, mod,
     *         ipthk, vthk, ipdly, dly, sdly, lat, lon, c, e, sw,
     *         klas, calr, xmgc, fmwt, xmwt, fmgc, ipcod, iscod, tpdly,
     *         revp, ndate, nhr, ns,
     *         prr, inpt, nreloc, inmain, injump, exdly, ibate,
     *         test53)
c cole sonafrank notes of august, 1988:
c restructured the code a bit while trying to follow the logic...
c first reading of station list for parameters in effect on ibate
      include 'params.inc'
      parameter (ndly = 11)
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      parameter (irecsz = 132)
      integer punt
      common /punt/ punt
      real lat,lon
      dimension ielv(nsn),nsta(nsn),dly(ndly,nsn),mod(nsn),ipthk(nsn)
      dimension vthk(2,nsn),ipdly(nsn),sdly(ndly,nsn),lat(nsn),lon(nsn)
      dimension c(nsn),e(nsn),sw(nsn),ndate(nsn),nhr(nsn),klas(5, nsn)
      dimension calr(5, nsn),fmgc(nsn),xmgc(nsn),ipcod(nsn),iscod(nsn)
      integer fmwt, xmwt
      dimension fmwt(nsn),xmwt(nsn),tpdly(2,nsn),prr(nsn)
      character*1 exdly(4, nsn), revp(6, nsn)
      character icard*(irecsz), icardup*(irecsz), dnstrg*(irecsz)
      character*4 nsta4, nsta*5, name, iast*1, rshft, nsta5
      character*1 ins, iew, izne
      integer iostat
      l = 1
      logfil = punt
c     write(logfil, '(a)') 'begin getsta.for - unit logfil'
  312 read(inpt, '(a)', end=999) icardup
      icard = dnstrg(icardup)
      name = rshft(icard(1:4))
      if(name(1:2) .eq. 'c*') goto 312
      if(name .eq. 'jump') then
        if(inpt .eq. inmain) then
          inpt = injump
c	  write (logfil, '(2a)') 'jump to station list: ', icardup(6:55)
          call openfl( injump, icardup(6:55), 'old', 'zero', 'readonly',
     *    'none', 0)
          goto 312
        else
          write(punt, 319) icard(1:55)
  319     format(' xxxerrorxxx can not nest jump statements ', /, 1x, a)
          stop
        endif
      endif
c
c ---   read primary station data values.
  320 continue
      if(nreloc .eq. 0) then
c        write (logfil,
c    *   '('' reading primary station values.  nreloc='',i2)')
c    *    nreloc
c        write (logfil, '(1x, a)') icard
        read(icard, 325)
     *    nsta4,llat,ins,alat,llon,iew,alon,ielv(l),
     *    mod(l),ipthk(l),vthk(1,l),vthk(2,l),ipdly(l),(dly(i,l),
     *    sdly(i,l),i=1,5), izne
  325   format(a4,i2,a1,f5.3,1x,i3,a1,f5.3,i5,i2,i1,2f4.2,i1,
     *    10f4.2,a1)
c ---   shift name to right and add component
	nsta(l)(1:4) = rshft(nsta4)

c for 4-letter codes, if the last letter is 'e' or 'n' then
c assume this is a horizontal component station!


        if(test53 .eq. 1.0) then
          if(nsta4(1:1) .ne. ' ') then
            if((nsta4(4:4) .eq. 'n') .or. (nsta4(4:4) .eq. 'e')) then
              izne = nsta4(4:4)
            endif
          endif
        endif

	if((izne .ne. 'e') .and. (izne .ne. 'n')) then
	  nsta(l)(5:5) = 'z'
	else
	  nsta(l)(5:5) = izne
	endif

c zero the extra delay models
        do 3255 i = 6, ndly-1
          dly(i,l) = 0.0
          sdly(i,l) = 0.0
3255    continue
      else
c        write (logfil,
c    *  '('' reading primary values without delays.  nreloc='', i2)')
c    *  nreloc
c ---   do not read delays if job is now recycling.
        read(icard, 326)
     *  nsta4,llat,ins,alat,llon,iew,alon,ielv(l),
     *  mod(l),ipthk(l),vthk(1,l),vthk(2,l),ipdly(l),izne
  326   format(a4,i2,a1,f5.3,1x,i3,a1,f5.3,i5,i2,i1,2f4.2,i1,40x,a1)
c ---   shift name to right and add component
	nsta(l)(1:4) = rshft(nsta4)

c for 4-letter codes, if the last letter is 'e' or 'n' then
c assume this is a horizontal component station!

        if(test53 .eq. 1.0) then
          if(nsta4(1:1) .ne. ' ') then
            if((nsta4(4:4) .eq. 'n') .or. (nsta4(4:4) .eq. 'e')) then
              izne = nsta4(4:4)
            endif
          endif
        endif

	if((izne .ne. 'e') .and. (izne .ne. 'n')) then
	  nsta(l)(5:5) = 'z'
	else
	  nsta(l)(5:5) = izne
	endif
      endif
c --- just read primary record
      if (nsta(l)(1:4) .eq. ' end') goto 991
      if (mod(l) .lt. 1) mod(l) = 1
      prr(l) = 0.1
      call fold2(lat(l),lon(l),llat,ins,alat,llon,iew,alon)
      if (l .eq. 1) then
        c(1) = 0.0
        e(1) = 0.0
      else
        call delaz(lat(1),lon(1),delt,deldeg,azz,lat(l),lon(l))
        c(l) = delt*sin(azz*rad)
        e(l) = delt*cos(azz*rad)
      endif
      if (ipthk(l) .ne. 2) ipthk(l) = 1
      if ((ipdly(l) .lt. 1) .or. (ipdly(l) .gt. 10)) ipdly(l) = 1
c
c --- read time dependent station data.
330   read(inpt, '(a)', iostat=iostat) icard
      icard = dnstrg(icard)
      if(icard(1:2) .eq. 'c*') goto 330
      if (iostat.eq.-1) then
        if(inpt .eq. injump) then
          inpt = inmain
          icard(1:4) = ' end'
          goto 390
        else
          goto 999
        endif
      endif
c decode a time dependent station record (icard) for station number l
      if(icard(5:5) .ne. '*') goto 395
      call rdtmdp(l, icard, nsta, ielv,
     *         lat, lon, c, e, sw,
     *         klas, calr, xmgc, fmwt, xmwt, fmgc, ipcod, iscod, tpdly,
     *         revp, ndate, nhr,
     *         exdly)
      if (ndate(l) .lt. ibate) goto 350
c
c --- station not expired, so read next record
  340 read(inpt, '(a)', iostat=iostat) icard
      icard = dnstrg(icard)
      if(icard(1:2) .eq. 'c*') goto 340
c --- shift name to right
      nsta4 = rshft(icard(1:4))
      if (iostat.eq.-1) then
        if(inpt .eq. injump) then
          inpt = inmain
          icard(1:4) = ' end'
          iast = ' '
          goto 990
        else
          goto 999
        endif
      endif
      iast = icard(5:5)
      izne = icard(80:80)
c     write(logfil, '(a, 1x, a, 1x, a, 1x, i5)')
c    *  nsta4, iast, nsta(l), l
      if (icard(1:4) .eq. ' end') goto 990
      if ((iast .ne. '*') .or. (nsta4 .ne. nsta(l)(1:4))) then
c --- this is NOT another time-dependent record
	if (iast .ne. '*') then
c         this SHOULD be a new station

c for 4-letter codes, if the last letter is 'e' or 'n' then
c assume this is a horizontal component station!

          if(test53 .eq. 1.0) then
            if(nsta4(1:1) .ne. ' ') then
              if((nsta4(4:4) .eq. 'n') .or. (nsta4(4:4) .eq. 'e')) then
                izne = nsta4(4:4)
              endif
            endif
          endif

	  if((izne .ne. 'e') .and. (izne .ne. 'n')) then
	    nsta5 = nsta4//'z'
	  else
	    nsta5 = nsta4//izne
	  endif
          if (nsta5 .ne. nsta(l))  then
c           this IS a new station
            goto 360
          else
c	    the same station AGAIN!  A clear error.
	    goto 392
	  endif
	else
c	  it is an error to have an * with a different station name
	  goto 392
	endif
      endif
c --- write extra time dependent record
      write(2,'(2a)') nsta(l), icard(6:lentru(icard))
      goto 340

c --- station expired so far, look for unexpired one.
  350 continue
c     write (logfil, '('' station expired so far'')')
      read(inpt, '(a)', iostat=iostat) icard
      icard = dnstrg(icard)
c     write (logfil,'(x,a)') icard
      if(icard(1:2) .eq. 'c*') goto 350
      nsta4 = rshft(icard(1:4))
      if (iostat.eq.-1) then
        if (inpt .eq. injump) then
          inpt = inmain
          goto 991
        else
          goto 999
        endif
      endif
      iast = icard(5:5)
      izne = icard(80:80)
      if (nsta4 .eq. ' end') goto 991
      if ((iast .ne. '*') .or. (nsta4 .ne. nsta(l)(1:4))) then
c --- this is NOT another time-dependent record
	if (iast .ne. '*') then
c         this SHOULD be a new station

c for 4-letter codes, if the last letter is 'e' or 'n' then
c assume this is a horizontal component station!

          if(test53 .eq. 1.0) then
            if(nsta4(1:1) .ne. ' ') then
              if((nsta4(4:4) .eq. 'n') .or. (nsta4(4:4) .eq. 'e')) then
                izne = nsta4(4:4)
              endif
            endif
          endif

	  if((izne .ne. 'e') .and. (izne .ne. 'n')) then
	    nsta5 = nsta4//'z'
	  else
	    nsta5 = nsta4//izne
	  endif
          if (nsta5 .ne. nsta(l))  then
c           this IS a new station
            goto 320
          else
c	    the same station AGAIN!  A clear error.
	    goto 393
	  endif
	else
c	  it is an error to have an * with a different station name
	  goto 392
	endif
      endif
      if(icard(5:5) .ne. '*') goto 395
      call rdtmdp(l, icard, nsta, ielv,
     *         lat, lon, c, e, sw,
     *         klas, calr, xmgc, fmwt, xmwt, fmgc, ipcod, iscod, tpdly,
     *         revp, ndate, nhr,
     *         exdly)
      if (ndate(l) .lt. ibate) goto 350

c --- found non-expired data for this one, save extra data.
      goto 340

  360 if (indexs .gt. 0) call iprst(l,llat,ins,alat,llon,iew,alon)
      l = l + 1
c     write (logfil, '('' 360: l='',i3)') l
      if (l .le. nsn) goto 320
      write(punt,370) nsn
  370 format(///,' xxxerrorxxx  station list exceeds max of ',
     *  i5, ' stations', /, ' so stop.')
      stop
  390 write(punt,380) l,nsta(l),name,iast
  380 format(///,'  xxxerrorxxx 390 in station list format, so stop',/,
     *  1x,i5,'th station ',2('"',a5,'" '),a1)
      stop
  392 write(punt,382) l, nsta(l), name, iast
  382 format(///,'  xxxerrorxxx 392 in station list format, so stop',/,
     *  1x,i5,'th station ',2('"',a5,'" '),a1)
      stop
  393 write(punt,383) l,nsta(l),name,iast
  383 format(///,'  xxxerrorxxx 393 in station list format, so stop',/,
     *  1x,i5,'th station ',2('"',a5,'" '),a1)
      stop
  394 write(punt,384) l,nsta(l),name,iast
  384 format(///,'  xxxerrorxxx 394 in station list format, so stop',/,
     *  1x,i5,'th station ',2('"',a5,'" '),a1)
      stop
  395 write(punt, 385) l, nsta(l), name, iast
  385 format(///,'  xxxerrorxxx in station list format, so stop',/,
     *  1x,i5,'th station ',2('"',a5,'" '),a1,
     *  ' missing time-dependent record')
 990  if (indexs .gt. 0) call iprst(l,llat,ins,alat,llon,iew,alon)
      l = l + 1
      if (iast .eq. '*') goto 394
 991  ns = l - 1
      write(2, 992)
 992  format(' end')
      return
 999  write(punt, 389)
 389  format (' xxxerrorxxx unexpected end of file on station list')
      stop 'error - getsta'
      end
c end getsta
c glob_new.for    []
      subroutine glob_new
c this subroutine attempt to find the global minimum when there is
c more than one local minimum at different depths.  j. c. lahr 4/22/87
c
      	include 'params.inc'
      	parameter (ndly = 11)
      	integer punt
      	common /punt/ punt
      	common /gg/ altla(3), altlo(3), altz(3), altrms(3), nsolut,
     * fla(20), flo(20), fz(20), frms(20), forg(20), maxf, altorg(3)
      	logical ingulf
      	common /gu/ ingulf
      	logical freor
      	common /hopq/ savla,savlo,savez,savor,freor
      	common /igl/ nglobalzs, globalzs(20)
      	common /imost/ test(100)
      	common /ihfgpq/ itest(100)
      	common /phoqn/ inst,knst
      	common /qmost/ wt(npa),z
      	real latep,lonep
      	common /qmost1/ lonep,ni,latep
      	common /ro/ yse,seorg,phi
      	common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      	common /zqr/ fno

c
c deep depth should be deep enough to catch deepest local mimima.
c we use 75. km in alaska.  add test(8) so that deepz is in km wrt
c the top of the model, not wrt sea level.

      	i = nglobalzs
	
20	savez = globalzs(i) + test(8)
	write(punt, *) 'savez =', savez
      	instsv = inst
      	if(instsv .eq. 8) freor = .false.
c set inst = 1 to fix depth
      	inst = 1
      	call quakes

c free up the depth
	inst = 0
      	savla = latep
      	savlo = lonep
	savez = z
	write(punt, *) 'savez for free run is ', savez
	call quakes

	write(punt, *) 'z =', z

c save this depth and rms
      	fz(i) = z
      	fla(i) = latep
      	flo(i) = lonep
        frms(i) = rms
c     	rmslm1 = frms(1)
c     	if(fno .gt. 0.) rmslm1 = sqrt(frms(1)**2 + (yse**2)/fno)

c for now, just force a start at every depth specified
	if (i .gt. 1) then
	  i = i - 1
	  goto 20
	endif

c this version would allow depths to be skipped
c	if (i .gt. 1) then
c	  find next depth to test
c	  do j = i - 1, 1, -1
c	    if (fz(i) .gt. globalzs(j)) then
c	      i = j
c	      goto 20
c	    else
c	      skip this starting depth
c	      frms(j) = 999999.
c	      fz(j) = -99.
c	    endif
c	  enddo
c	endif

c now decide which was the best solution
	write(punt, *) frms
	write(punt, *) fz
	bestrms = frms(1)
	ibest = 1
	do i = 2, nglobalzs
	  if (frms(i) .lt. bestrms) then
	    bestrms = frms(i)
	    ibest = i
	  endif
	enddo

c get ready to make final run again
	savela = fla(ibest)
	savelo = flo(ibest)
	savez = fz(ibest)
	inst = instsv
	return
	end

c global.for    []
      subroutine global
c this subroutine attempt to find the global minimum when there are
c more than one local minimum at different depths.  j. c. lahr 4/22/87
c
      include 'params.inc'
      parameter (ndly = 11)
      real bat2,bon2,lat,lon
      logical good, eoff, supout
      integer punt
      common /punt/ punt
      common /gg/ altla(3), altlo(3), altz(3), altrms(3), nsolut,
     * fla(20), flo(20), fz(20), frms(20), forg(20), maxf, altorg(3)
      character*4 iahead*60, msta*5, nsta*5, icard*110
      common /char/ iahead, msta(npa), nsta(nsn), icard
      character*4 ipro, ichec, evtype*1, evstat*1, root*256
      common /dbhip/ ipro, ichec, evtype, evstat, root
      common /dmost/ ipun,ivlr,blank
      common /dhil/ iq,ilat,kms,icat
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /dhio/ nedit,good
      common /dio/ rmsmx,presmx,sresmx,noutmx,nimx,iprun,semx
      common /dph/ noswt, eoff
      common /ghnq/ iexit
      common /gmost/ az(npa)
      logical ingulf
      common /gu/ ingulf
      common /hf/ cp(72),sp(72)
      common /hl/ nwad, tslope, tsqsl
      common /hpn/ wslope
      logical freor
      common /hopq/ savla,savlo,savez,savor,freor
      common /iclmpq/ lat(nsn),lon(nsn)
      common /logfil/ logfil
      common /idno/ ksel,ksort
      common /ilv/ c(nsn), e(nsn)
      common /imost/ test(100)
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /ihfgpq/ itest(100)
      common /lm/ mapend
      common /ocfn/ adj,seh,lat1,bat2,lon1,bon2,igap
      character*1 isnla, isnlo, krm1, krm2
      common /ocfn1/ isnla, isnlo, krm1, krm2
      common /ohq/ gap, supout
      common /omnfh/ dmin,dmin3,sminp
      common /ph/ nfirst
      common /pno/ lph,keyph(npa),nsum
      character*4 krms
      common /po/ krms(npa)
      common /pfnoqv/ kdate,khrmn,khr
      common /pgqv/ w(npa)
      common /phoqn/ inst,knst
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /pnoqtx/ ldx(npa)
      character*4 krmp
      common /pnoqv/ krmp(npa)
      common /pot/ ts(npa),ihr,model(npa), keyphi(npa)
      common /qgo/ org
      common /qmost/ wt(npa),z
      real latep,lonep
      common /qmost1/ lonep,ni,latep
      common /ro/ yse,seorg,phi
      common /reloc/ irelo, nreloc
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      common /zqr/ fno
c
      ingulf = .false.
c
c deep depth should be deep enough to catch deepest local mimima.
c we use 75. km in alaska.  add test(8) so that deepz is in km wrt
c the top of the model, not wrt sea level.
c     test(42) = 75. is the default.
      deepz = test(42) + test(8)
c
c a solution is first computed for z = 0 and for z free, starting at deepz.
c
c if the deep solution converges to a depth less than some cutoff,
c cutz, then the solution with smaller rms is taken as the best.
c we use 20 in alaska.  add test(8) so that deepz is in km wrt
c the top of the model, not wrt sea level.
c     test(27) = 20. is the default value
      cutz = test(27) + test(8)
c
c if the deep solution converges to a depth greater than cutz, then
c another solution is found starting at shalz = cutz/2.0.
c     shalz = 9. (old value)   new default value is 10.
      shalz = cutz/2.0
c
c emperical tests indicate that for alaska data this algorithm with
c the default values of test(27) and test(42) usually
c finds the global minimum.  for other networks, different values of
c deepz, cutz, and shalz may need to be used.
c
      if(iprn .ge. 5) then
        write(punt,91)
   91   format(' call global')
c show elapsed time
        call timit(1)
      endif
      instsv = inst
      if(instsv .eq. 8) freor = .false.
c set inst = 1 to fix depth
      inst = 1
c set savez to surface
      savez = 0.
      call quakes
c save this location and rms
      fla(1) = latep
      flo(1) = lonep
      fz(1) = z
      frms(1) = rms
      forg(1) = org
      rmslm1 = frms(1)
      if(fno .gt. 0.) rmslm1 = sqrt(frms(1)**2 + (yse**2)/fno)
c set savez to deep depth
      savez = deepz
      call quakes
c save this location and rms
      fla(2) = latep
      flo(2) = lonep
      fz(2) = z
      frms(2) = rms
      forg(2) = org
c now free up the depth and begin at deep depth again
      inst = 0
      savla = latep
      savlo = lonep
      call quakes
c if solution is in the gulf of alaska, then save it and return
      if(ingulf) then
	if(iprn .ge. 5) then
          write(punt, '(a)')
     *    'Global deep-starting free solution ends up in gulf.'
        endif
        altla(1) = latep
        altlo(1) = lonep
        altz(1) = z
        altrms(1) = rms
        altorg(1) = org
        nsolut = 1
        return
      endif
c save this location and rms
      fla(3) = latep
      flo(3) = lonep
      fz(3) = z
      frms(3) = rms
      forg(3) = org
      rmslm3 = frms(3)
      if(fno .gt. 0.) rmslm3 = sqrt(frms(3)**2 + (yse**2)/fno)
      if(fz(3) .lt. cutz) then
c the solution must be shallow.
        maxf = 3
        if(fz(3) .le. 0.1) then
c in this case there is only one minimum, and it is at the surface
          altla(1) = fla(3)
          altlo(1) = flo(3)
          altz(1) = fz(3)
          altrms(1) = frms(3)
          altorg(1) = forg(3)
          nsolut = 1
          return
        else if(rmslm1 .lt. frms(3)) then
c in this case near the surface is significantly better,
c so run free location to derive information for output
          savez = fz(1)
          savla = fla(1)
          savlo = flo(1)
          inst = 0
          call quakes
c save this location and rms
          altla(1) = latep
          altlo(1) = lonep
          altz(1) = z
          altrms(1) = rms
          altorg(1) = org
          nsolut = 1
          return
        else if(frms(1) .gt. rmslm3) then
c the deeper solution is significantly better, so use it
          altla(1) = latep
          altlo(1) = lonep
          altz(1) = z
          altrms(1) = rms
          altorg(1) = org
          nsolut = 1
          return
        else
c in this case the difference in rms is less than limit
c these will be joined if there is not a boundary in between them.
          nsolut = 2
c the prefered solution will have the lesser rms
          if(frms(1) .le. rms) then
c the surface is at least slightly better or equal to deeper z.
c save current solution as secondary minimum.
            altla(2) = latep
            altlo(2) = lonep
            altz(2) = z
            altrms(2) = rms
            altorg(2) = org
c recompute information for output.
            savez = fz(1)
            savla = fla(1)
            savlo = flo(1)
            if (instsv .ne. 8) freor = .true.
            inst = 9
            call quakes
c save this location and rms
            altla(1) = latep
            altlo(1) = lonep
            altz(1) = z
            altrms(1) = rms
            altorg(1) = org
            return
          else
c the deeper solution is at least slightly better.
            altla(1) = latep
            altlo(1) = lonep
            altz(1) = z
            altrms(1) = rms
            altorg(1) = org
            altla(2) = fla(1)
            altlo(2) = flo(1)
            altz(2) = fz(1)
            altrms(2) = frms(1)
            altorg(2) = forg(1)
            return
          endif
        endif
      else
c the final depth is greater than cutz km when starting at deep depth
        altla(1) = latep
        altlo(1) = lonep
        altz(1) = z
        altrms(1) = rms
        altorg(1) = org
c compute a fixed depth solution at shallow depth
        savez = shalz
        savla = fla(1)
        savlo = flo(1)
        inst = 1
        call quakes
c save this location and rms
        fla(4) = latep
        flo(4) = lonep
        fz(4) = z
        frms(4) = rms
        forg(4) = org
        maxf = 4
c compare rms at surface and shallow depth and start free solution at better.
        if(frms(4) .le. frms(1)) then
c shallow depth is better than surface so start free solution there.
          inst = 0
          savla = latep
          savlo = lonep
          call quakes
c save this location and rms
          fla(5) = latep
          flo(5) = lonep
          fz(5) = z
          frms(5) = rms
          forg(5) = org
          rmslm5 = frms(5)
          if(fno .gt. 0.) rmslm5 = sqrt(frms(5)**2 + (yse**2)/fno)
          maxf = 5
        else
c surface is better than shallow depth so start free soltuion there.
          savez = fz(1)
          savla = fla(1)
          savlo = flo(1)
          inst = 0
          call quakes
c save this location and rms
          fla(5) = latep
          flo(5) = lonep
          fz(5) = z
          frms(5) = rms
          forg(5) = org
          rmslm5 = frms(5)
          if(fno .gt. 0.) rmslm5 = sqrt(frms(5)**2 + (yse**2)/fno)
          maxf = 5
        endif
c       if(abs(fz(5) - fz(3)) .lt. 5.) then
        if(abs(fz(5) - fz(3)) .lt. deepz/10.) then
c the depths are so close together that these solutions are equivalent.
          nsolut = 1
          if(frms(5) .le. frms(3)) then
            altla(1) = fla(5)
            altlo(1) = flo(5)
            altz(1) = fz(5)
            altrms(1) = frms(5)
            altorg(1) = forg(5)
            return
          else
c compute fixed depth solution for output values
            savla = fla(3)
            savlo = flo(3)
            savez = fz(3)
            if (instsv .ne. 8) freor = .true.
            inst = 9
            call quakes
c save this location and rms
            altla(1) = latep
            altlo(1) = lonep
            altz(1) = z
            altrms(1) = rms
            altorg(1) = org
            return
          endif
        else
c the depths are far enough appart that there may be two solutions.
          if(frms(5) .gt. frms(3)) then
c the deep solution is at least slightly better
            altla(1) = fla(3)
            altlo(1) = flo(3)
            altz(1) = fz(3)
            altrms(1) = frms(3)
            altorg(1) = forg(3)
            altla(2) = fla(5)
            altlo(2) = flo(5)
            altz(2) = fz(5)
            altrms(2) = frms(5)
            altorg(2) = forg(5)
            nsolut = 1
            savla = fla(3)
            savlo = flo(3)
            savez = fz(3)
            if (instsv .ne. 8) freor = .true.
            inst = 9
            call quakes
          else
c the shallow solution is at least slightly better
            altla(2) = fla(3)
            altlo(2) = flo(3)
            altz(2) = fz(3)
            altrms(2) = frms(3)
            altorg(2) = forg(3)
            altla(1) = fla(5)
            altlo(1) = flo(5)
            altz(1) = fz(5)
            altrms(1) = frms(5)
            altorg(1) = forg(5)
          endif
          if(frms(5) .gt. rmslm3) then
c the deep solution is significantly better
            nsolut = 1
            return
          else if(frms(3) .gt. rmslm5) then
c the shallow solution is significantly better
            nsolut = 1
            return
          else
c report two solutions
            nsolut = 2
            return
          endif
        endif
      endif
      end
c end global
c halfsp.for    []
      subroutine halfsp(m, t, tz, vh, x, xz, y, yz, zz)
      common /logfil/ logfil
c solve for best halfspace solution, given m p arrival times and the
c origin time.
c              mx               maximum number of equations (m .eq. mx)
      parameter (mx = 10)
c              a(6)             storage mode 1 version of sum of squares
      real     a(6)
c              m                (input) number of p arrival times
      integer  m
c              c(4, mx)         normal equation coeficients and constants
      real     c(4, mx)
c              ev(3, 3)         eigenvector matrix
      real     ev(3, 3)
c              r(3)             normal equation constants in rotated system
      real     r(3)
c              s(3)             eignevalues
      real     s(3)
c              t(mx)            (input) arrival time at station m
      real     t(mx)
c              tz               origin time
      real     tz
c              vh               halfspace velocity
      real     vh
c              v(4, 4)          symetric tensor  - upper 3x3 to be diagonalized
      real     v(4, 4)
c              x(mx)            (input) x coordinate of station m
      real     x(mx)
c              xyt(3)           eq coordinates in original system
      real     xyt(3)
c              xytr(3)          eq coordinates in rotated system
      real     xytr(3)
c              xz               (output) x coordinate of earthquake
      real     xz
c              y(mx)            (input) y coordinate of station m
      real     y(mx)
c              yz               (output) y coordinate of earthquake
      real     yz
c              zz               (output) z coordinate of earthquake
      real     zz
c
      if (m .gt. mx) then
        print *, 'halfsp can not have more than ', mx, ' equations,'
        print *, 'so ', m, ' is too many.'
        stop
      endif
      n = 3
      vsq = vh*vh
      nfix = 1
      if (tz .eq. 99999.) nfix = 0
c
c set up cooeficients of equations, c
c
cd    print *, 'original equations for ', m, ' arrival times:'
      do 20 i = 1, m-1
        c(1, i) = x(m) - x(i)
        c(2, i) = y(m) - y(i)
        if (tz .eq. 99999.) then
c         in this case, tz is unknown
          c(3, i) = -vsq*(t(m) - t(i))
          c(4, i) = ( x(m)*x(m) + y(m)*y(m) - t(m)*t(m)*vsq
     *             -x(i)*x(i) - y(i)*y(i) + t(i)*t(i)*vsq )/2.
        else
c         in this case, tz is given
          c(3, i) = ( x(m)*x(m) + y(m)*y(m) - t(m)*t(m)*vsq
     *             -x(i)*x(i) - y(i)*y(i) + t(i)*t(i)*vsq )/2.
     *             +vsq*( t(m) - t(i) )*tz
          c(4, i) = 0.0
        endif
cd      print *, i, (c(j, i), j = 1, 4)
20    continue
c
c compute sum of squares matrix from coeficients
c
      do 25 j = 1, 4-nfix
        do 25 k = j, 4-nfix
          v(j, k) = 0.0
          do 23 i = 1, m-1
23          v(j, k) = v(j, k) + c(j, i)*c(k, i)
25    v(k,j) = v(j,k)
cd    print *, 'normal equations:'
cd    do 26 j = 1, 4-nfix
cd      print *, (v(j, k), k = 1, 4-nfix)
cd26   continue
c
c compute determinant
c
cd    det = deter(v, 4, n-nfix)
cd    print *, 'determinant = ', det
c
c compute eigenvalues and eigenvectors of nxn tensor v
c ev is the nxn eigenvector tensor
      call eigen1(a, 3, 4, n-nfix, ev, s, v, 0.0)
c
cd    write (logfil, *) ' eigenvectors'
cd    do 45 i = 1, 3
cd      write (logfil, *) (ev(i, j), j = 1, 3)
cd45   continue
cd    write (logfil, *) ' eigenvalues'
cd    write (logfil, *) s
c
c
c compute response vector in rotated coordinates
c
      do 30 j= 1, n-nfix
        r(j) = 0.0
        do 30 k = 1, n-nfix
          r(j) = r(j) + ev(k, j)*v(k, 4-nfix)
30    continue
cd    print *, 'response vector in rotated coordinates'
cd    print *, r
c
c compute the eq location in rotated system
c
      do 50 i = 1, n-nfix
        if(s(i) .gt. .000001) then
          xytr(i) = r(i)/s(i)
        else
          xytr(i) = 0.0
        endif
50    continue
c
c rotate location back into original coordinates
c
      do 55 j = 1, n-nfix
        xyt(j) = 0.0
        do 55 k = 1, n-nfix
55        xyt(j) = xyt(j) + ev(j,k)*xytr(k)
c
c define output variables
c
      xz = xyt(1)
      yz = xyt(2)
      if (tz .eq. 99999.) tz = xyt(3)
      arg =  vsq*(t(m) - tz)**2 - (x(m) - xz)**2 - (y(m) - yz)**2
      if (arg .ge. 0.) then
        zz = sqrt (arg)
      else
cd      print *, 'arg = ', arg, ' so z set to zero'
        zz = 0.
      endif
cd    print *, 'xz        yz        zz          tz'
cd    print *, xz, yz, zz, tz
c
      return
      end
c end halfsp
c hycrt.for    []
      subroutine hycrt (mod, itable, vpvs)
c    written by f. klein for hypoinverse
c    modified slightly for hypoellipse - j.c. lahr
c    further improvements - j.a. snoke - january 1991 - spherical earth
c--reads a travel time table for a linear gradient crust model
c--called by hypoinv
c--data used only for linear gradient travel time tables
c--ln is number of models
      parameter (ln=3)
c  Changed nlyr from 12 to 20 on 2/9/2000 jcl
      parameter (nlyr=20)
      integer punt
      common /punt/ punt
      character modnam*20
      common /mc/ modnam(ln)
c                            model name or label
      common /m/ lay(ln)
c                            number of layers or v-d points
      common /m/ d(nlyr,ln)
c                            depth to layer top or velocity point
      common /m/ vel(nlyr,ln)
c                            layer or point velocity
      common /m/ thk(nlyr,ln)
c                            thickness of homogeneous layer
      common /m/ vsq(nlyr,ln)
c                            squared velocity of homogeneous layer
      common /m/ modtyp(ln)
c                            mdl type (-1=undef, 0=grad, 1=homo layer)
c                            the number of linear-grad models allowed
      integer kt
      logical gd1,gd2,gz1,gz2
      common /m/ redv(ln)
c                            one over the reducing velocity
      common /m/ nz(ln)
c                            number of depth grid points
      common /m/ nz1(ln),dz1(ln),nz2(ln),dz2(ln)
c                            depth grid params
      common /m/ nd(ln)
c                            number of distance grid points
      common /m/ nd1(ln),dd1(ln),nd2(ln),dd2(ln)
c                            dist grid params
      common /m/ gd1(ln),gd2(ln),gz1(ln),gz2(ln)
c                            grid flags
      common /m/ kdhr(ln,28)
c                            for each model (ln) and for each depth,
c                            the distance at which horizontal ray emerges
c                            (in units of 0.1 km)
      common /m/ kt(ln,28,42)
c                            the travel times (up to 28 z's & 42 dist's)
c--identify the model mod as a linear gradient table
      modtyp(mod)=0
c
c--read header info
      rewind(itable)
      read (itable,1000) modnam(mod),lay(mod),redv(mod),vpvs,radius
1000  format (a20, i5, 2f10.5, f10.3)
      write(punt, '(1x, a20, 5x, a, f10.5)')
     *  modnam(mod), 'vp/vs =', vpvs
c
c--read depths & velocities of model
      read (itable,1001) (d(i,mod),i=1,lay(mod))
      read (itable,1001) (vel(i,mod),i=1,lay(mod))
1001  format (3x, 15f7.2)
c
c--read distance & depth grid info
      read (itable,1003) dd1(mod),nd1(mod),dd2(mod),nd2(mod)
      read (itable,1003) dz1(mod),nz1(mod),dz2(mod),nz2(mod)
1003  format (3x, 2(f10.4, i5))
      write(punt, '('' z '', 15f7.2)') (d(i, mod), i = 1, lay(mod))
      write(punt, '('' v '', 15f7.2)') (vel(i, mod), i = 1, lay(mod))
      write(punt, '('' dd '', 2(f10.4, i5))')
     *  dd1(mod), nd1(mod), dd2(mod), nd2(mod)
      write(punt, '('' dz '', 2(f10.4, i5))')
     *  dz1(mod), nz1(mod), dz2(mod), nz2(mod)
      nd(mod)=nd1(mod)+nd2(mod)+1
      nz(mod)=nz1(mod)+nz2(mod)+1
c
c--read reduced travel times, grouped by depth
      do 20 j=1,nz(mod)
      read (itable,1004) kdhr(mod,j)
1004  format (20x,i10)
20    read (itable,1005) (kt(mod,j,i),i=1,nd(mod))
1005  format (15i8)
      gd1(mod)=nd1(mod).ne.0 .and. dd1(mod).ne.0
      gd2(mod)=nd2(mod).ne.0 .and. dd2(mod).ne.0
      gz1(mod)=nz1(mod).ne.0 .and. dz1(mod).ne.0
      gz2(mod)=nz2(mod).ne.0 .and. dz2(mod).ne.0
      return
      end
c end hycrt
c hypot.for    []
      real function hypot(a, b)
c
c-------- bruce julian
c
c
c-------- hypot - calculates euclidian distance, accurately and
c            avoids overflow
c
      real a, b
      real abs, l, s, t, sqrt
      l = abs(a)
      s = abs(b)
      if (s .le. l) goto 1
         t = s
         s = l
         l = t
   1  if (l .ne. 0.0) goto 2
         hypot = 0.0
         return
   2  s = s/l
      hypot = l*sqrt(s*s+1.0)
      return
      end
c end hypot
c hyset.for    []
      subroutine hyset(nttab, z1)
c    written by f. klein for hypoinverse
c    modified slightly for hypoellipse - j.c. lahr
c    further improvements - j.a. snoke - january 1991 - spherical earth
c--initialize values for a given depth, z1.
c--uses a condensed and reduced travel time table generated
c--by the program ttgen.
      parameter (ln=3)
c                            the number of linear-grad models allowed
c  Changed nlyr from 12 to 20 on 2/9/2000 jcl
      parameter (nlyr=20)
      character modnam*20
      common /mc/ modnam(ln)
c                            model name or label
      common /m/ lay(ln)
c                            number of layers or v-d points
      common /m/ d(nlyr,ln)
c                            depth to layer top or velocity point
      common /m/ vel(nlyr,ln)
c                            point velocity
      common /m/ thk(nlyr,ln)
c                            thickness of homogeneous layer
      common /m/ vsq(nlyr,ln)
c                            squared velocity of homogeneous layer
      common /m/ modtyp(ln)
c                            mdl type (-1=undef, 0=grad, 1=homo layer)
c--data used only for linear gradient travel time tables
      integer kt
      logical gd1,gd2,gz1,gz2
      common /m/ redv(ln)
c                            one over the reducing velocity
      common /m/ nz(ln)
c                            number of depth grid points
      common /m/ nz1(ln),dz1(ln),nz2(ln),dz2(ln)
c                            depth grid params
      common /m/ nd(ln)
c                            number of distance grid points
      common /m/ nd1(ln),dd1(ln),nd2(ln),dd2(ln)
c                            dist grid params
      common /m/ gd1(ln),gd2(ln),gz1(ln),gz2(ln)
c                            grid flags
      common /m/ kdhr(ln,28)
c                            for each model (ln) and for each depth,
c                            the distance at which horizontal ray emerges
c                            (in units of 0.1 km)
      common /m/ kt(ln,28,42)
c                            the travel times (up to 28 z's & 42 dist's)
      common /tlook/ vh(ln)
c                            velocity at hypocenter
      common /tlook/ dhrz(ln)
c                            dist at which horiz ray emerges
      common /tlook/ tz(ln,42)
c                            tt at each distance point
      common /tlook/ dtz(ln,42)
c                            tt depth deriv at each distance point
      data scfac/2000./
c--loop over the ln possible models, skipping loop if model not used.
      do 35 md=1,nttab
c--perform some preliminary calcs and interpolations which depend only on
c  depth. do this for each linear grad model once.
      nzm=nz(md)
      nz1m=nz1(md)
      nz2m=nz2(md)
      dz1m=dz1(md)
      dz2m=dz2(md)
      laym=lay(md)
      nd1m=nd1(md)
c--perform depth interpolation first
c--find depth index = index of table entry nearest hypocenter
c--also depth spacing h and fraction of interval x
      temp=nz1m*dz1m
      if ((.not.gz1(md) .or. z1.gt.temp) .and. gd2(md)) then
c--hypo in lower part of table
        i=nz1m+(z1-temp)/dz2m+1.5001
        if (i.ge.nzm) i=nzm-1
        if (i.lt.nz1m+2) i=nz1m+2
        h=dz2m
        x=(z1-temp)/dz2m-(i-nz1m-1)
        i1 = nz1m + (z1-temp)/dz2m + 1.0001
        if (i1.ge.nzm) i1 = nzm - 1
        y = (z1 - temp)/dz2m - (i1 - nz1m - 1)
      else
c--hypo in upper part of table
        i=z1/dz1m+1.5001
        if (i.gt.nz1m) i=nz1m
        if (i.lt.2) i=2
        h=dz1m
        x=z1/dz1m-(i-1)
        i1 = z1/dz1m + 1.0001
        y = z1/dz1m - (i1 - 1)
      end if
c--find exact velocity at hypocenter
      l = 1
      do 20 k=1,laym
      if (d(k,md).lt.z1) l=k
20    continue
      vh(md)=vel(laym,md)
      if (l.lt.laym) vh(md)=vel(l,md)+(vel(l+1,md)-vel(l,md))
     *   *(z1-d(l,md))/(d(l+1,md)-d(l,md))
c--interpolate dist at which a horizontal ray emerges
      dhrz(md)=(kdhr(md,i1)+y*(kdhr(md,i1+1)-kdhr(md,i1)))*.1
c--depth interpolation
      if (z1.le.temp+dz2m*nz2m) then
c--use 3 point interpolation
        ca=x*.5*(x-1.)
        cb=1.-x**2
        cc=x*.5*(x+1.)
        da=(x-.5)/h
        db=-2.*x/h
        dc=(x+.5)/h
      else
c--use linear extrapolation
        ca=0.
        cb=1.-x
        cc=x
        da=0.
        dc=1./h
        db=-dc
        dhrz(md)=1000.
      end if
c--interpolate tt and its depth derivative for all distance grid points
      temp=nd1m*dd1(md)
      do 30 j=1,nd(md)
      dx=(j-1)*dd1(md)
      if (j.gt.nd1m+1) dx=temp+(j-nd1m-1)*dd2(md)
      tz(md,j)=(ca*kt(md,i-1,j)+cb*kt(md,i,j)
     * +cc*kt(md,i+1,j)+32000.)/scfac+dx*redv(md)
30    dtz(md,j)=(da*kt(md,i-1,j)+db*kt(md,i,j)+dc*kt(md,i+1,j))/scfac
35    continue
      return
      end
c end hyset
c hytab.for    []
      subroutine hytab(md, dx, travti, ain, dtdd, dtdh, vt, vs)
c    written by f. klein for hypoinverse
c    modified slightly for hypoellipse - j.c. lahr
c    further improvements - j.a. snoke - january 1991 - spherical earth
c--md is model for this station
c--dx is station distance
      parameter (ln=3)
c                            the number of linear-grad models allowed
c
c  Changed nlyr from 12 to 20 on 2/9/2000 jcl
      parameter (nlyr=20)
      character modnam*20
      common /mc/ modnam(ln)
c                            model name or label
      common /m/ lay(ln)
c                            number of layers or v-d points
      common /m/ d(nlyr,ln)
c                            depth to layer top or velocity point
      common /m/ vel(nlyr,ln)
c                            point velocity
      common /m/ thk(nlyr,ln)
c                            thickness of homogeneous layer
      common /m/ vsq(nlyr,ln)
c                            squared velocity of homogeneous layer
      common /m/ modtyp(ln)
c                            mdl type (-1=undef, 0=grad, 1=homo layer)
c--data used only for linear gradient travel time tables
      integer kt
      logical gd1,gd2,gz1,gz2
      common /m/ redv(ln)
c                            one over the reducing velocity
      common /m/ nz(ln)
c                            number of depth grid points
      common /m/ nz1(ln),dz1(ln),nz2(ln),dz2(ln)
c                            depth grid params
      common /m/ nd(ln)
c                            number of distance grid points
      common /m/ nd1(ln),dd1(ln),nd2(ln),dd2(ln)
c                            dist grid params
      common /m/ gd1(ln),gd2(ln),gz1(ln),gz2(ln)
c                            grid flags
      common /m/ kdhr(ln,28)
c                            for each model (ln) and for each depth,
c                            the distance at which horizontal ray emerges
c                            (in units of 0.1 km)
      common /m/ kt(ln,28,42)
c                            the travel times (up to 28 z's & 42 dist's)
      common /tlook/ vh(ln)
c                            velocity at hypocenter
      common /tlook/ dhrz(ln)
c                            dist at which horiz ray emerges
      common /tlook/ tz(ln,42)
c                            tt at each distance point
      common /tlook/ dtz(ln,42)
c                             tt depth deriv at each distance point
c--now find travel time, dtdr, dtdz, and angle of emergence for
c--start distance interpolation
      ndm=nd(md)
      nd1m=nd1(md)
      nd2m=nd2(md)
      dd1m=dd1(md)
      dd2m=dd2(md)
      temp=nd1m*dd1m
c--find distance index = index of table entry nearest the station dist
c--also distance spacing h and fraction of interval x
      if ((.not.gd1(md) .or. dx.gt.temp) .and. gd2(md)) then
c--hypo in far part of table
        j=nd1m+(dx-temp)/dd2m+1.5001
        if (j.ge.ndm) j=ndm-1
        if (j.lt.nd1m+2) j=nd1m+2
        h=dd2m
        x=(dx-temp)/dd2m-(j-nd1m-1)
      else
c--hypo in near part of table
        j=dx/dd1m+1.5
        if (j.gt.nd1m) j=nd1m
        if (j.lt.2) j=2
        h=dd1m
        x=dx/dd1m-(j-1)
      end if
c--distance interpolation
      if (dx.le.temp+dd2m*nd2m) then
c--use 3 point interpolation
        anin=x*.5*(x-1.)
        cb=1.-x**2
        cc=x*.5*(x+1.)
        da=x-.5
        db=-2.*x
        dc=x+.5
      else
c--use linear extrapolation
        anin=0.
        cb=1.-x
        cc=x
        da=0.
        dc=1.
        db=-1.
      end if
c--interpolate tt and its 2 derivatives
      travti=anin*tz(md,j-1)+cb*tz(md,j)+cc*tz(md,j+1)
      dtdh  =anin*dtz(md,j-1)+cb*dtz(md,j)+cc*dtz(md,j+1)
      dtdd=(da*tz(md,j-1)+db*tz(md,j)+dc*tz(md,j+1))/h
      if (dtdd.lt.0.) dtdd=0.
c--calc emergence angle
c-- for spherical models, the following equation would be more accurate:
c--   anin=vh(md)*dtdd*radius/(radius-depth)
c   however in most situations this change will be quite small.
      anin=vh(md)*dtdd
      if (anin.gt..99999) anin=.99999
      ain=57.29578*asin(anin)
      if (dx.lt.dhrz(md)) ain=180.-ain
      if (dx.ge.dhrz(md)) anin = -anin
      vt = vel(1, md)
      vs = vh(md)
      return
      end
c end hytab
c init.for    [pc]
      subroutine init
c* (pc
$notruncate
c* pc)
c get filenames, open files and write greeting for hypoellipse
      include 'params.inc'
      character*256 root
      common /dhip/ inpt,isa,ilis,inmain,injump
      character*4 ipro, ichec, evtype*1, evstat*1
      common /dbhip/ ipro, ichec, evtype, evstat, root
      common /hf/ cp(72),sp(72)
      real*8 time1, time2
      common /hop/ time1,time2,nopu,notim
      common /it/ vpvs(lmax2),vpvsm(mmax+3),nttab
      common /lm/ mapend
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      integer punt
      common /punt/ punt
      integer print_unit, summary_unit, arch_unit
      parameter (print_unit   = 9)
      parameter (summary_unit = 4)
      parameter (arch_unit    = 11)
      character*1 print_type, summary_type, arch_type
      parameter (print_type   = 'o')
      parameter (summary_type = 's')
      parameter (arch_type    = 'a')
      character openstat*7
      openstat = 'unknown'
c
c open the files for hypoellipse (returns root name)
        call opfls(inpt, inmain, root, nttab)
        call timit(0)
c get current date
        call jdate(irmo, irdy, iryr, ihr, imn, isec)
c do not print output conversion errors
        call erset (63,.true.,.false.,.false.,.false.)
        write(punt, 20) nsn, npa, iryr, irmo, irdy, ihr, imn
20      format (
c* (pc
     *' *** Hypoellipse: PC/Non-Xpick/Y2K version 3.9 11/1/2001  ***',
c* pc)
c* (unix
c     *' *** Hypoellipse: Unix/Non-Xpick/Y2K version 3.9 9/8/2003  ***',
c* unix)
     */, '  Configured for up to ', i4, ' stations in station list',
     */, '  and up to ', i4, ' records per earthquake.',
     */, '  Run on ', 2(i2.2, '/'), i2.2, ' at ', i2.2, ':', i2.2)
c calculate tables to be used in focal mechanism plots
        do 30 i=1,72
c         pi=i*.0349066*2.5
          pi = i*0.0872665
          cp(i)=cos(pi)
          sp(i)=sin(pi)
30      continue
        ichec = 't s '
        mapend=0
        nr = 0
        time1 = 0.0d0
      return
      end
c end initial
c input1.for    []
      subroutine input1
c input station list, crustal model, test and control values
c     to change max number of stations from nsn to some other number,
c       151 for example:
c                1)  change all occurrences of nsn) to 151)
c                2)  change nmax from nsn to 151 in next data statement.
c                3)  note:  be sure the new dimension is a unique number
c                      which has not been previously used so that
c                      subsequent changes in dimension will not be
c                      difficult to make.
c --- added inmain in place of the constant 8.  its now in the /dhip/
c --- common block and also passed to getsta().
      include 'params.inc'
      parameter (ndly = 11)
      real lat,lon,latr,lonr
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      character*4 iscan
      character*4 iahead*60, msta*5, nsta*5, icard*110, dnstrg*110
      character*50 uacal, uaini, dlyfil, icardup*110
      integer punt
      common /punt/ punt
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /dmost/ ipun,ivlr,blank
      common /dhil/ iq,ilat,kms,icat
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /dbiln/ ioldq
      character*1 iclass
      common /dbio/ iclass(0:4)
      common /dhin/ iglob, zup, zdn
      common /dhio/ nedit,good
      common /dio/ rmsmx,presmx,sresmx,noutmx,nimx,iprun,semx
      character*1 bksrc
      common /dipu/ ipkdly, igsum, bksrc
      logical medmag
      common /dinx/ imag,imagsv,medmag
      common /dit/ tid(lmax,lmmax),did(lmax,lmmax),lowv,modin(mmax+3)
      common /dix/ iexcal, uacal
      common /igl/ nglobalzs, globalzs(20)
      common /imost/ test(100)
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /idt/ v(lmax2)
      integer fmwt, xmwt
      common /ioxl/ fmwt(nsn),xmwt(nsn)
      common /logfil/ logfil
      common /ihfgpq/ itest(100)
      character*1 iqcls
      common /il1/ iqcls
      common /ilpu/ sw(nsn),ndate(nsn),nhr(nsn),mno(nsn),ielv(nsn)
      common /ilmpu/ ns
      common /iclmpq/ lat(nsn),lon(nsn)
      common /ilpx/ calr(5, nsn),fmgc(nsn),xmgc(nsn)
      common /in/ irmo,iryr,odmag
      character*1  magke
      common /in1/ magke
      common /idno/ ksel,ksort
      common /iot/ flt(2,nsn),thks(nsn)
      common /iox/ prr(nsn),iuses
      common /iiq/ wf(51)
      logical fsteq
      character*1 revp, exdly
      common /ip/ latr,lonr,tpdly(2,nsn),infil,iofil,indexs,
     *  iscod(nsn),ipcod(nsn), fsteq, exdly(4, nsn), revp(6, nsn)
      common /ip1/ rsew(4)
      common /int/ thk(lmax+1),lbeg(mmax),lend(mmax),vi(lmax),vsq(lmax),
     *  vsqd(lmmax,lmax),f(lmax,lmmax),kl,ivway,sthk(mmax),sthk1(mmax),
     *  sdly(ndly,nsn)
      common /ilt/ vthk(2,nsn),ipdly(nsn),mod(nsn),ipthk(nsn)
      common /ilv/ c(nsn),e(nsn)
      common /it/ vpvs(lmax2),vpvsm(mmax+3),nttab
      character*4 mtyp*2
      common /iu/ mtyp(mmax)
      common /ix/      ir,qspa(9,40)
      common /lvz/ jref(lmax)
      common /reloc/ irelo, nreloc
      common /rioq/ been,damp,dmpinc,igo
      character*4 inote(100)
      dimension atest(100),d(lmax2)
c initialize test variables used throughout the program
      character*1 ins, iew
      integer iostat
      data d/lmax2*0./
      data atest/1.78, 5.0, 0.0, 0.0, -99.0, 0., 10., 0., 0., 50.,
     * 50.0,100.0 ,50.0,50.0, 10.0,50.0, 2.0, 50.0, 0.05,0.05,
     * 9.0, 35., 0.7, 35., 40., 0.0025, 20.0, 0., -.1, 0.0,
     * -1.15, 2.00, 0.0, 0.0, 0.001, 100.0, 3.0, 0.0, 1.0, 0.007,
     * 0.0, 75.0, 2*0., 0.1379, 0.0, 0.0, 6.5, 2*0., 1000.,
     * 2800., 1.0, 200., 19., 45*0.0/
      data inote/100*'    '/
      data uaini/' '/, kl/0/
      data nglobalzs/0/, globalzs/20*0.0/
      fsteq = .true.
200   continue
      nkrst = 0
      been = 0.0
205   read(inpt, '(a)', iostat = iostat) icardup
      icard = dnstrg(icardup)
      if (iostat.eq.-1) then
        if (ilis.gt.0) then
          write(punt, 220) inpt
          write(logfil, 220) inpt
220       format(' subroutine input1 found end of file on unit ', i5)
        endif
        if(inpt .eq. inmain) then
          write(logfil, 225)
          write(punt, 225)
225       format(' xxxerrorxxx abnormal termination while ',
     *    'reading control file.')
          stop 'abort from input1'
        else
c         return to reading from unit inmain
          inpt = inmain
          goto 205
        endif
      endif
      iscan = icard(1:4)
      if (icard .eq. ' ') goto 205
      if (iscan(1:4) .eq. 'stop') then
        write(punt, '(a)') ' stop processing with stop instruction'
        stop 'normal termination'
      endif
      if (iscan(1:2) .eq. 'c*') then
        write(punt, 235) icard
235     format(1x, a)
        goto 205
      endif
      if (iscan .eq. 'head') then
        iahead = icard(19:lentru(icard))
        if (ilis.gt.0) write(punt, '(//35x,a)') iahead
        goto 205
      endif
c
c remove everything on record beginning with an ! mark
      do 2515 i = 1, 110
        if(icard(i:i) .eq. '!') then
          if(i .eq. 1) goto 205
          icard = icard(1:i-1)
          goto 2517
        endif
2515  continue
c
2517  if (iscan .eq. 'arc ') then
        if(ilis .gt. 0)
     *  write(punt, '(a, a)')
     *    ' open archive file named:  ', icardup(19:78)
        close(11)
        call openfl(   11, icardup(19:78), 'unknown', 'null',
     *  'none', 'noprint', 0)
        goto 205
      endif
c
      if (iscan .eq. 'uofa') then
        uacal = icardup(19:78)
        goto 205
      endif
      if (iscan .eq. 'jump') then
        if (inpt .eq. inmain) then
          close (unit = injump,iostat = iostat)
          if (ilis.gt.0) write(punt, '(1x, a)') icardup
          call openfl(injump, icardup(6:55), 'old', 'zero', 'readonly',
     *    'none', 0)
          inpt = injump
        else
          write(logfil, 260) inmain
          write(punt, 260) inmain
260       format(' xxxwarningxxx can not nest jumps from control file,'
     *    ,/, ' returning to main input file (unit ',i2,')')
          inpt = inmain
        endif
        goto 205
      endif
c because standard fortran 77 to can not read an internal file with
c free format, file 14 must be used in the following code!
      rewind 14
      write(14, '(a)') icard(19:80)
      rewind 14
      if (iscan .eq. 'rese') then
c       read(icard(19:80), *, err = 5243) j, ax
        read(14, *, iostat = iostat) j, ax
        if (iostat.ne.0) then
          ierror = 1
          goto 5243
        endif
        if (( j .gt. 100) .or. ( j .lt. 1 )) then
          ierror = 1
          goto 530
        endif
        inote(j) = '****'
        test(j) = ax
        goto 205
      endif
c
      if ((iscan .eq. 'crus') .or. (iscan .eq. 'velo')) then
c       read(icard(19:80), *, err = 5243)  ax, bx, cx
        read(14, *, iostat = iostat)  ax, bx, cx
        if (iostat.ne.0) then
          ierror = 2
          goto 5243
        endif
        if (nkrst .eq. 0) then
          lh = lmax + 2
          do  290 k = 1,lh
            d(k) = 0.0
            v(k) = -1.
290       continue
          do 295 k = 1, mmax
            mtyp(k) = 'ps'
295       continue
        endif
        nkrst = nkrst + 1
        if (nkrst .gt. lmax) then
          write(logfil, 310) icard
          write(punt, 310) icard(1:30)
310       format(' too many velocity layers. ',
     *    ' layer ', a, ' was not used.')
        else
          v(nkrst) = ax
          d(nkrst) = bx
          vpvs(nkrst) = cx
        endif
        goto 205
      endif
c
      if (iscan .eq. 'weig') then
c       read(icard(19:80), *, err = 5243) ax, bx, cx
        read(14, *, iostat = iostat) ax, bx, cx
        if (iostat.ne.0) then
          ierror = 3
          goto 5243
        endif
        if ((ax .lt. 1.).or.(bx .lt. 1.).or.(cx .lt. 1.)) then
          if (ilis.gt.0) write(punt, 320) icard
320       format(' weight option parameters represent the ratio of', /,
     *    ' the standard error of weight code 1, 2, and 3 ', /,
     *    ' readings to the standard error of 0 weight code', /,
     *    ' readings.  none can be less than 1. so stop.', /,
     *    1x, a)
          stop 'abort from input1.'
        endif
        rsew(2) = ax
        rsew(3) = bx
        rsew(4) = cx
        goto 205
      endif
c
      if (iscan .eq. 'vari') then
c       read(icard(19:80), *, err = 5243) ivlr, ivway, lowv
        read(14, *, iostat = iostat) ivlr, ivway, lowv
        if (iostat.ne.0) then
          ierror = 4
          goto 5243
        endif
        goto 205
      endif
c
      if (iscan .eq. 'miss') then
c       read(icard(19:80), *, err = 5243) kms
        read(14, *, iostat = iostat) kms
        if (iostat.ne.0) then
          ierror = 5
          goto 5243
        endif
        if(ilis .gt. 0) then
          if (kms .eq. 0) write(punt,335)
335       format(' scan for missing stations.  code = 0')
          if (kms .ne. 0) write(punt,340) kms
340       format(' do not scan for missing stations.  code = ', i3)
        endif
        goto 205
      endif
c
      if (iscan .eq. 'prin') then
c       read(icard(19:80), *, err = 5243) iprn
        read(14, *, iostat = iostat) iprn
        if (iostat.ne.0) then
          ierror = 6
          goto 5243
        endif
        if ((iprn .lt. -2) .or. (iprn .gt. 5)) then
          ierror = 6
          goto 530
        endif
        goto 205
      endif
c
      if ((iscan .eq. 'punc') .or. (iscan .eq. 'summ')) then
c       read(icard(19:80), *, err = 5243) ipun
        read(14, *, iostat = iostat) ipun
        if (iostat.ne.0) then
          ierror = 7
          goto 5243
        endif
        if ((ipun .lt. -4) .or. (ipun .gt. 4)) then
          ierror = 7
          goto 530
        endif
        goto 205
      endif
c
      if (iscan .eq. 'magn') then
c       read(icard(19:80), *, err = 5243) j
        read(14, *, iostat = iostat) j
        if (iostat.ne.0) then
          ierror = 8
          goto 5243
        endif
        iabj = iabs(j)
        if (((iabj .gt. 4) .and. (iabj .lt. 10)) .or.
     *       (iabj .gt. 14)) then
          ierror = 8
          goto 530
        endif
        imagsv = j
        imag = iabj
        medmag = .false.
        if(imag .gt. 9) then
          medmag = .true.
	  imag = imag - 10
        endif
        iuses = 0
        if (j .lt. 0) iuses = 1
        goto 205
      endif
c
      if ((iscan .eq. 'qual') .or. (iscan .eq. 'tabu')) then
c       read(icard(19:80), *, err = 5243) j
        read(14, *, iostat = iostat) j
        if (iostat.ne.0) then
          ierror = 9
          goto 5243
        endif
        if ((j.lt. -4 ).or.(j.gt. 4)) then
          ierror = 9
          goto 530
        endif
        ioldq = 0
        if (j .lt. 0) ioldq = 1
        iq = iabs(j)
        iqcls = iclass(iq)
        goto 205
      endif
c
365   if (iscan .eq. 'sort') then
c       read(icard(19:80), *, err = 5243) ksort
        read(14, *, iostat = iostat) ksort
        if (iostat.ne.0) then
          ierror = 10
          goto 5243
        endif
        if(ilis .gt. 0) then
          if (ksort .eq. 0) write(punt,370)
370       format(' sort stations by distance.  code = 0' )
          if (ksort .ne. 0) write(punt,375) ksort
375       format(' do not sort stations by distance.  code = ', i3)
        endif
        goto 205
      endif
c
      if (iscan .eq. 'stan') then
        write(punt,395)
395     format(46h change all program options to standard values )
        do 400 i = 1, 100
          inote(i) = '    '
          test(i) = 1.23456
400     continue
        do 405 i = mmax+1, mmax+3
          modin(i) = 0
405     continue
        bksrc = ' '
        iexcal = 0
        ns = 0
        rsew(2) = 5.
        rsew(3) = 10.
        rsew(4) = 20.
        igsum = 1
        ivlr = 0
        lowv = 1
        kms = 1
        ir = 0
        iprn = 1
        ilis = 1
        ipun = 0
        iuses = 0
        imagsv = 0
	medmag = .false.
        imag = 0
        iq = 2
        iqcls = iclass(2)
        ioldq = 0
        ksort = 0
        inpt = inmain
        icat = 1
        nedit = 0
        iprun = 1
        ipkdly = 1
        irelo = 0
        iglob = 1
        uacal = uaini
        nglobalzs = 0
	do i = 1, 20
	  globalzs(i) = 0.0
	enddo
        goto 205
      endif
      if (iscan .eq. 'blan') then
c source code to use if phase record has blank field
        do 412 i = 19, 80
          if(icard(i:i) .ne. ' ') goto 413
412     continue
413     bksrc = dnstrg(icard(i:i))
        if (ilis .gt. 0) write(punt, 414) bksrc
414     format(' arrival-time record blank source fields will be',
     *  ' assumed to be source "', a, '"')
        goto 205
      endif
c
      if (iscan .eq. 'comp') then
c       read(icard(19:80), *, err = 5243) ksel
        read(14, *, iostat = iostat) ksel
        if (iostat.ne.0) then
          ierror = 12
          goto 5243
        endif
        if (ilis .gt. 0) then
          if (ksel .ne. 0) write(punt,415)
415       format(42h do not compress printed output.  code = 1)
          if (ksel .eq. 0) write(punt,420) ksel
420       format(33h compress printed output.  code =,i3)
        endif
        goto 205
      endif
c
425   if (iscan .eq. 'igno') then
c       read(icard(19:80), *, err = 5243) igsum
        read(14, *, iostat = iostat) igsum
        if (iostat.ne.0) then
          ierror = 13
          goto 5243
        endif
        if (igsum .eq. 0) then
          if (ilis.gt.0) write(punt, 430) igsum
430       format(/,' ignore sum code = ', i5, /,
     *    ' (ignore starting locations on summary records)')
        else
          if (ilis.gt.0) write(punt, 435) igsum
435       format(/,' ignore sum code = ', i5, /,
     *    ' (use starting locations from summary records)')
        endif
        goto 205
      endif
c
      if (iscan .eq. 'cali') then
        if (ilis.gt.0) write(punt,445)
445     format(/' input extra calibration curves.')
c       read(icard(19:80), *, err = 5243) ir
        read(14, *, iostat = iostat) iexcal
        if (iostat.ne.0) then
          ierror = 14
          goto 5243
        endif
        if (iexcal .lt. 1 .or. iexcal .gt. 9) then
          write(punt, '(a, a, /, a)') ' error in: ', icard,
     *    ' the number of extra calibrations must be between 1 and 9'
          write(logfil, '(a, a, /, a)') ' error in: ', icard,
     *    ' the number of extra calibrations must be between 1 and 9'
          stop 'abort from input1'
        else
          do 460 i = 1,iexcal
            read(inpt,450) (qspa(i,k),k = 1,40)
450         format(20f4.2)
            ip8 = i + 8
            if (ilis.gt.0) write(punt,455) ip8,(qspa(i,k),k = 1,40)
455         format(/,' qspa(', i2, ') ', 20f5.2, /, 10x, 20f5.2)
460       continue
        endif
        goto 205
      endif
c
465   if (iscan .eq. 'begi') goto 540
      if (iscan .eq. 'arri') then
        if(inpt .eq. inmain) then
          goto 580
        else
          write(logfil, 470) inmain
          write(punt, 470) inmain
470       format(' xxxerrorxxx input data structure error,' /,
     *    ' the arrival times next record must be in main', /,
     *    ' input stream (unit ',i2,').')
          stop 'abort from input1'
        endif
      endif
c
      if (iscan .eq. 'debu') then
c       read(icard(19:80), *, err = 5243) nedit
        read(14, *, iostat = iostat) nedit
        if (iostat.ne.0) then
          ierror = 15
          goto 5243
        endif
        if (nedit .eq. 0) then
          if (ilis .gt. 0) write(punt,480) nedit
480       format(33h do not use debug option.  code =,i3)
          goto 205
        endif
        read(inpt,485) icard
485     format(a)
        if (ilis .gt. 0) write(punt,490) icard
490     format(' debug option record is',/1x,a)
        read(icard, 495) rmsmx,presmx,sresmx,noutmx,nimx,semx
495     format(5x,3(5x,f5.2),2(5x,i5),5x,f5.0)
        goto 205
      endif
      if (iscan .eq. 'cal ') then
        if (ilis.gt.0) write(punt, 496)
496     format(' xxxx warning xxxx  calibration option no longer valid')
        goto 205
      endif
      if (iscan .eq. 'resi') then
c       read(icard(19:80), *, err = 5243) iprun
        read(14, *, iostat = iostat) iprun
        if (iostat.ne.0) then
          ierror = 16
          goto 5243
        endif
        goto 205
      endif
      if (iscan .eq. 'sele') then
c       read(icard(19:80), *, err = 5243) ipkdly
        read(14, *, iostat = iostat) ipkdly
        if (iostat .ne. 0) then
          ierror = 17
          goto 5243
        endif
        if(ipkdly .le. 0) then
          if (ilis.gt.0) write(punt,515)
515       format(' select delays based on sub. usedly regions.')
          if(ipkdly .lt. 0) then
            read(inpt, '(a)') dlyfil
c           call openfl(iunit,  ifile, istat,  izero, ishr,
c    *      iform, irecl)
            call openfl(   17, dlyfil, 'old', 'null', 'none',
     *      'none', 0)
            if (ilis.gt.0) write(punt, 5155) dlyfil
5155        format(' delay cylinders defined in the file:', /, 1x, a)
          endif
        else
          write(punt,516)
516       format(' select delays based on closest station.')
        endif
        goto 205
      endif
      if (iscan .eq. 'relo') then
c       read(icard(19:80), *, err = 5243) irelo
        read(14, *, iostat = iostat) irelo
        if (iostat.ne.0) then
          ierror = 18
          goto 5243
        endif
        if (ilis.gt.0) write(punt,522) irelo
522     format(' relocate events ', i2, ' times, each time', /,
     *  ' updating the station delays for delay model 1.')
        goto 205
      endif
      if (iscan .eq. 'newg') then
	read(14, *, iostat = iostat) (globalzs(i), i = 1, 20)
	do i = 20, 1, -1
	  if(globalzs(i) .ne. 0.0) then
	    nglobalzs = i
	    goto 5225
	  endif
	enddo
	write(punt, *)
     *    ' new global option record included with no non-zero depths'
	write(logfil, *)
     *    ' new global option record included with no non-zero depths'
	goto 205
5225    write(punt, 5226) nglobalzs, (globalzs(i), i = 1, nglobalzs)
        write(logfil, 5226) nglobalzs, (globalzs(i), i = 1, nglobalzs)
5226	format(' running new global option with ', i3,
     *    ' starting depths of ', /, 2x, 10f10.2)
	goto 205
      endif
      if (iscan .eq. 'glob') then
c       read(icard(19:80), *, err = 5243) iglob
        read(14, *, iostat = iostat) iglob
        if (iostat.ne.0) then
          ierror = 19
          goto 5243
        endif
        if (ilis .gt. 0) then
          if (iglob .eq. 0) then
            write(punt, 5241)
            write(logfil, 5241)
5241        format(' global minimum search turned on')
          else
            write(punt, 5242)
            write(logfil, 5242)
5242        format(' global minimum search turned off')
          endif
        endif
        goto 205
      endif
      if (iscan .eq. 'cons') then
        read(14, *, iostat = iostat) j
        if (iostat .ne. 0) then
          ierror = 21
          goto 5243
        endif
        ilis = j
        write(logfil, '(a, i5)') ' constants print = ', ilis
        goto 205
      endif
      if (iscan .eq. 'dela') then
c       read an additional delay model
        read(14, *, iostat = iostat) nmodel
        if (iostat .ne. 0) then
          ierror = 22
          goto 5243
        endif
        call adddly(nmodel)
        goto 205
      endif
5243  write(punt,525)  icard, ierror
      write(logfil,525)  icard, ierror
525   format(' xxxerrorxxx control record in improper format.',
     * /,1x, a,/' so stop (error #', i3, ')')
      stop 'abort from input1'
530   write(logfil,535) ierror, icard
      write(punt,535) ierror, icard
535   format(
     * ' this record skipped because j is improper value (error #',
     * i3,'). ', /, 1x, a)
      goto 205
c input station list
540   if (ilis.gt.0) write(punt,550) iahead
550   format(' list of stations available for these solutions',
     *   /,35x,a)
560   continue
c     read(icard(19:80), *, err = 5243) indexs, ibate
      read(14, *, iostat = iostat) indexs, ibate
      if(ibate .lt. 19000000) then
	write(logfil, *) 'xxxerrorxxx century not included on ',
     *    'begin station list record.  date = ', ibate
	write(punt, *) 'xxxerrorxxx century not included on ',
     *    'begin station list record.  date = ', ibate
	stop
      endif
      if (iostat.ne.0) then
        ierror = 20
        goto 5243
      endif
      write(punt, '(1x, a)') icard(1:80)
      ihrmn = 0
c indexs = 0 or 1 -> stations list (index option has been removed)
c negative indexs -> no print
      if (indexs .eq. 0) indexs = 1
      if (ilis.gt.0) then
        if (iabs(indexs) .eq. 1) write(punt,565) indexs
565     format(' station list code =',i5)
        if (iabs(indexs) .ne. 1) then
          write(punt,'(a, i5)')
     *    ' xxxerrorxxx station list code must be 1, not ', indexs
          write(logfil,'(a)')
     *    ' xxxerrorxxx station list code must be 1, not ', indexs
          stop 'abort from input1'
        endif
        write(punt,575) ibate
575     format(' set up for events starting on ',i8,/,
     * ' name latitude  longitude  elev p thickness p p  pdy1 sdy1',
     * '  pdy2 sdy2  pdy3 sdy3  pdy4 sdy4  pdy5 sdy5     calr xmgc',
     * ' mgwt fmgc wt ',/,
     * ' * continuation  record *      thk 1   2  mod dly          ',
     * '                                            sys          ',
     * '           ps',/,
     * '  polarity stawt teldy code altdy cnyrmody hr')
      endif
      rewind 2
      l = 1
      if (test(2) .eq. 1.23456) test(2) = atest(2)
      if (test(53) .eq. 1.23456) test(53) = atest(53)
c     read station list
c     write(logfil, '(a)') ' just about to call getsta - unit logfil'
      call getsta(   indexs, nsta, ielv, mod,
     *         ipthk, vthk, ipdly, dly, sdly, lat, lon, c, e, sw,
     *         klas, calr, xmgc, fmwt, xmwt, fmgc, ipcod, iscod, tpdly,
     *         revp, ndate, nhr, ns,
     *         prr, inpt, nreloc, inmain, injump, exdly, ibate,
     *         test(53))
      infil = 3
      iofil = 2
      goto 205
c----
580   continue
c write test variables
      if (ilis.gt.0) write(punt,595) iahead
595   format(1h ,35x,a)
      call jdate(irmo, irdy, iryr, ihr, imn, isec)
      iseed = irmo*irdy*iryr + imn + isec
      dum = ran3(-iseed)
      do 600 i = 1, 55
        if (test(i).eq.1.23456) test(i) = atest(i)
        itest(i) = test(i) + sign(0.0005,test(i))
600   continue
      if(test(34) .ne. 0.0) test(34) = 1.0
      itest(38) = abs(test(38)) +0.0005
      if ((iglob .eq. 0) .and. (test(47) .ne. 0.0)) then
        write(punt, 590)
        write(logfil, 590)
590     format(' xxxerrorxxx the global option is not compatible with',/
     *         '   the option to fix the hypocenter on a plane. ',/
     *         '   either set test(47) to 0.0 or turn off the',
     *         ' global option')
        stop 'abort from input1'
      endif
c define jeffreys weighting funct.
      if (iprn .ge. 3) write(punt, 605) test(20)
605   format(/,47h jeffreys weighting funct         i      weight,
     *  7h   mu =,f10.3)
      do 615 i = 1,51
        wf(i) = (1.0 + test(20))/(1.0+test(20)*exp(((i-1)*0.1)**2/2.0))
        if (iprn .ge. 3) then
          write(punt,610) i,wf(i)
610       format(30x,i5,f12.3)
        endif
615   continue
      if (iprn .ge. 3) write(punt,595)
      if (ilis.gt.0) then
        write(punt, 620)  ( i,atest(i), test(i), inote(i), i =  1, 10)
        write(punt, 625)  ( i,atest(i), test(i), inote(i), i = 11, 20)
        write(punt, 630)  ( i,atest(i), test(i), inote(i), i = 21, 30)
        write(punt, 635)  ( i,atest(i), test(i), inote(i), i = 31, 40)
        write(punt, 640)  ( i,atest(i), test(i), inote(i), i = 41, 50)
        write(punt, 642)  ( i,atest(i), test(i), inote(i), i = 51, 55)
      endif
620   format(/7x, 'test variables', 14x, 'description', /, 5x,
     *'standard  reset to',
     */i3,2f10.4,a4,'ratio of p-wave velocity to s-wave velocity.',
     */i3,2f10.4,a4,'lt 0 no elev cor/ =0 use 1st vel/ gt 0 use this.',
     */i3,2f10.4,a4,'first trial latitude in degrees.',
     */i3,2f10.4,a4,'first trial longitude in degrees.',
     */i3,2f10.4,a4,'first trial depth in kilometers, unless = -99.',
     */i3,2f10.4,a4,'sphere rad for aux rms values. if neg cont iterat',
     *'ion at most neg point.',
     */i3,2f10.4,a4,'minimum number of first motions required to plot.',
     */i3,2f10.4,a4,'elevation of top of layered models (km).',
     */i3,2f10.4,a4,'if 0 allow neg depths in summary and archive ',
     * 'files.',
     */i3,2f10.4,a4,'apply distance weighting on this iteration.')
625   format(
     *i3,2f10.4,a4,'xnear = greatest distance with weight of 1.0',
     */i3,2f10.4,a4,'xfar = least distnace with weight of 0.0',
     */i3,2f10.4,a4,'apply azimuthal weighting on this iteration.',
     */i3,2f10.4,a4,'weight out large residuals on this iteration.',
     */i3,2f10.4,a4,'give zero weight to residuals gt this.',
     */i3,2f10.4,a4,'apply boxcar weighting on this iteration.',
     */i3,2f10.4,a4,'give zero weight to residuals gt this*stand. dev.',
     */i3,2f10.4,a4,'begin jeffreys weighting on this iteration.',
     */i3,2f10.4,a4,'use jeffreys weighting only if rms gt this.',
     */i3,2f10.4,a4,'mu of jeffreys weighting funct.')
630   format(
     * i3,2f10.4,a4,'maximum number of iterations.',
     */i3,2f10.4,a4,'limit change in focal depth to this amount (km).',
     */i3,2f10.4,a4,'if delz would make z neg, set delz = -this*z ',
     *'(km).',
     */i3,2f10.4,a4,'limit change in epicenter to this. (km).',
     */i3,2f10.4,a4,'fix depth if epicentral change gt this. (km).',
     */i3,2f10.4,a4,'stop iterating if square of adjustment lt this.',
     */i3,2f10.4,a4,'global opt: if deep solution z > this, continue',
     *' with z 1/2 way to surface.',
     */i3,2f10.4,a4,'for fixed hypo on plane, set = plunge azimuth.',
     *'  if neg. continue as free sol.',
     */i3,2f10.4,a4,'set std err of res=+this if degrees of freedom =',
     *'0 or =-this if this lt 0.'
     */i3,2f10.4,a4,'dip of plunge vector for epi. fixed on plane. ',
     *' see test(28) & (47) also.')
635   format(
     * i3,2f10.4,a4,'duration magnitude c1, constant.',
     */i3,2f10.4,a4,'duration magnitude c2, *log((f - p)*fmgc).',
     */i3,2f10.4,a4,'duration magnitude c3, *delta.',
     */i3,2f10.4,a4,'if not 0, scale the normal equations.',
     */i3,2f10.4,a4,'minimum damping of normal equations.  ',
     */i3,2f10.4,a4,'maximum first trial depth if computed from p-',
     *'arrival times.',
     */i3,2f10.4,a4,'if termination occurs before this iteration, set ',
     *'iteration number to this and continue.',
     */i3,2f10.4,a4,'if this =1, run all with and then without s/ =2,',
     *'run with s/ =3, run without s/ =4, fix hypo', /,
     * 27x, '   / neg, use s to fix origin.',
     */i3,2f10.4,a4,'multiply the s and s-p weights by this factor.',
     */i3,2f10.4,a4,'duration magnitude c4, *depth.')
640   format(
     * i3,2f10.4,a4,'if this =1, print opt. ge 1, & summary ',
     *'opt. =+ or -1, then write sum. record each itteration.',
     */i3,2f10.4,a4,'global opt: deep starting z wrt top of model.',
     */i3,2f10.4,a4,'duration magnitude c5, *(log((f - p)*fmgc)**2).',
     */i3,2f10.4,a4,'if =1 rerun debug eqs with critical sta/ =2 ',
     *'continue iter with crit sta.',
     */i3,2f10.4,a4,'x scale factor for focal mechanism plot.',
     */i3,2f10.4,a4,'xfar set ge dist of test(46)th station + 10.  if ',
     *'lt 0 then fill gap.',
     */i3,2f10.4,a4,'weight for fix on plane.  see test(28) and (30).',
     */i3,2f10.4,a4,'half-space velocity for first trial location.',
     */i3,2f10.4,a4,'if .ne. 0 calculate vp/vs ratio; if abs val >1 ma',
     *'ke wadati plot; if neg, use wadati origin in solution.',
     */i3,2f10.4,a4,'for exploring rms space, compute this number of',
     *' fixed depth solutions (up to 22).')
642   format(
     */i3,2f10.4,a4,'for epicentral distance beyond this, use first',
     *' travel-time table.',
     */i3,2f10.4,a4,'Wood Anderson static magnification assumed for',
     *' local magnitude determination.',
     */i3,2f10.4,a4,'if .eq. 1 stations with 4-letter codes ending',
     *' e or n treated as horizontals.',
     */i3,2f10.4,a4,'if 1st computed trial location > this (km) from',
     *' closest station, start at closest station.',
     */i3,2f10.4,a4,'assumed century for events without summary',
     *' record.')

      j = iq
      if (ioldq .eq. 1) j = -iq
      if (ilis.gt.0) then
        write(punt, 645) rsew
645     format(/,
     *  ' weight option - relative standard errors for code:   0'
     *, '      1      2      3', /,
     *         '                                                   ',
     *    4f7.3, /)
        write(punt,650) iprn,ipun,imagsv,j
650     format(' printer option ', i4, 2x, ' summary option  ', i4, 5x,
     * ' magnitude option ', i4, 5x, ' tabulation option ',i2,//
     * ' no event output  -2', /,
     * ' one line/eq      -1', /,
     * ' final solution    0   no sum records     0      use xmag',
     * 12x, '0      no summary        0', /
     * ' one line per iter 1   summary records    1      use fmag',
     * 12x, '1      a                 1', /
     * ' sta res each iter 2   sum + archive file 2      use (xmag+fm'
     *,'ag)/2   2      a + b             2', /
     * 70h regres each iter  3   archive file       3      prefer fmag /
     *xmag   3,
     * 6x, 19ha,b + c           3,/
     * '                       "corrected" input  4  ',
     * '    prefer xmag /fmag   4      a,b,c + d         4',
     * /, 23x, '                          if neg use fms not fmp',
     *  5x, 'positive/q from std errors', /
     * 45x, 31x, 'negative/q from sol+sta')
        if(uacal .ne. ' ') write(punt, 654) uacal
654     format(' u of a cal data file:  ', a)
      endif
      if (nedit .ne. 0 .and. ilis .gt. 0)
     * write(punt,655) nedit,rmsmx,presmx,sresmx,
     *                              noutmx,nimx,semx
655   format(' debug code = ',i1,'.    debug events have:',
     * /,'                                           rms .gt. ',f5.2,
     * /,'                                 or a p res is .gt. ',f5.2,
     * /,'                                or an s res is .gt. ',f5.2,
     * /,'        or the no. of readings weighted out is .gt. ',i5,
     * /,'                or the number of iterations is .gt. ',i5,
     * /,'  or the max single variable standard error is .gt. ',f5.0,
     * /,'                     or the change in depth is .lt. ',f5.2)
c----
      if (nkrst .eq. 0) then
        if (kl .eq. 0) then
          write(punt, '(a,/,a)')
     *    ' xxx error xxx:  ',
     *    ' at least one layered crustal model must be specified.'
          stop 'abort from input1'
        endif
        goto 790
      endif
c scan and print velocity structure data read in above
      kl = 1
      lbeg(1) = 1
      if (test(1) .le. 0.0) then
        write(logfil, 660) test(1)
        write(punt, 660) test(1)
660     format(' xxx error xxx', /,
     *  ' test(1) is the vp/vs ratio and may not equal ', f10.2)
        stop 'abort from input1'
      endif
      if (ivlr.ne.0.and.ilis.gt.0) write(punt,665) ivlr,ivway
665   format(54h variable layer option is used.  layer to be varied is
     *,7h layer ,i3,9h    vmod=,i3)
      if (lowv .eq. 0.and.ilis.gt.0) write(punt,670)
670   format(47h do not make compensating change in layer below,
     * 16h variable layer.)
      if (lowv .eq. 1.and.ilis.gt.0) write(punt,671)
671   format(' make compensating change in layer below',
     * ' variable layer.')
      if (iprun .ne. 1.and.ilis.gt.0) write(punt,680)
680   format(' residual option is in effect.  for residuals between + '
     *  ,49hand - 2.25 only the absoulte value rounded to the,
     * /,41h      nearest 0.5 second will be printed.)
      if (ilis.gt.0) write(punt,685) kl
685   format(//7x,14hvelocity model ,i3
     */29h   layer  velocity     depth ,
     * 21h  thickness    vpvs  /11x,25hkm/sec       km        km/)
      lp1 = lmax + 1
      do 690 l = 1, nkrst
690   if (vpvs(l) .eq. 0.0) vpvs(l) = test(1)
      do 785 l = 1,lp1
        thk(l) = d(l+1)-d(l)
        if (thk(l).lt.0.0) thk(l) = 1000.000
        if ( (d(l) .ne. 0.) .or. (l .eq. 1) ) goto 770
c found depth of zero, so begin new model
        lend(kl) = l-1
        if (lend(kl) .ne. lbeg(kl)) goto 700
        write(punt,695)
        write(logfil,695)
695     format(//' xxxerrorxxx half space model not allowed.'/
     *         ' so stop.')
        stop 'abort from input1'
c check for variation in vp/vs
700     vpvsm(kl) = vpvs(lbeg(kl))
c        print *, 'vpvs for model ', kl, ' is ', vpvsm(kl)
c        print *, 'lbeg = ', lbeg
c        print *, 'lend = ', lend
c        print *, 'vpvs = ', vpvs
        do 705 i = lbeg(kl) + 1,lend(kl)
          if (vpvsm(kl) .ne. vpvs(i)) goto 710
705     continue
        goto 745
c define s model
c need to define s model
710     vpvsm(kl) = 0.0
        ishft = lend(kl) + 1 - lbeg(kl)
        nkrstn = nkrst + ishft
        if (nkrstn .le. lmax) goto 720
        write(punt,715) nkrstn,lmax
        write(logfil,715) nkrstn,lmax
715     format(//' xxxerrorxxx ',i5,' exceeds',i5,' the max. number'
     *  ,' of layers.',/,' so stop.')
        stop 'abort from input1'
720     ntomv = nkrst - lend(kl)
        nkrst = nkrstn
        if (ntomv .eq. 0) goto 730
c shift down remaining models
        do 725 i = 1,ntomv
          ii = lend(kl) + ntomv + 1 - i
          iii = nkrstn + 1 - i
          v(iii) = v(ii)
          d(iii) = d(ii)
          vpvs(iii) = vpvs(ii)
725     continue
730     do 735 i = 1,ishft
          ii = lbeg(kl) - 1 + i
          iii = lend(kl) + i
          v(iii) = v(ii)/vpvs(ii)
          d(iii) = d(ii)
          vpvs(iii) = 0.0
735     continue
        if (ilis.gt.0) write(punt,740)
740     format(//' the next model is for s only:')
        mtyp(kl) = 'p '
        mtyp(kl+1) = 's '
745     if (v(l) .eq. -1.) goto 790
        if (v(l).lt.0.01) then
          write(logfil, 750) v(l), kl
          write(punt, 750) v(l), kl
750       format(' xxxerrorxxx velocity may not be less than .01',
     *    ' but was ', f10.4, ' in model ', i4, ', so stop!')
          stop 'abort from input1'
        endif
        kl = kl+1
        if (kl.le.mmax) goto 760
        write(punt,755) mmax
        write(logfil,755) mmax
755     format(//' xxxerrorxxx too many velocity models. only ',i3,
     *  ' models allowed. so stop.')
        stop 'abort from input1'
760     lbeg(kl) = l
        if (ilis.gt.0) write(punt,685) kl
        thk(l) = d(l+1) - d(l)
        if (d(l) .eq. 0.0) goto 770
        write(punt,765)
        write(logfil,765)
765     format(//' xxxerrorxxx each structure must begin at a',
     *         ' depth of 0.0 km.  so stop.')
        stop 'abort from input1'
770     vi(l) = 1/v(l)
        vsq(l) = v(l)**2
        if (ilis.gt.0) write(punt,775) l,v(l),d(l),thk(l),vpvs(l)
775     format(1x,i5,2x,4f10.3)
        if ((l-lbeg(kl)+1).le.lmmax) goto 785
        write(punt,780) kl, lmmax
        write(logfil,780) kl, lmmax
780     format(//' xxxerrorxxx velocity model',i3,' has too many ',
     *         'layers.  only',i3,' layers are allowed.',/,
     *         ' so stop.')
        stop 'abort from input1'
785   continue
790   continue
c     take care of travel time tables.
      if( nttab .ne. 0) then
	do 792 i = 1, nttab
          call hycrt(i, i+20, vpvsm(i+mmax) )
	  if ( (vpvsm(i+mmax) .lt. 0.) .and. (nttab .lt. i+1) ) then
            write(logfil, '(3a)')
     *        ' negative vp/vs in third travel-time table requires',
     *        ' that a forth table be include for s travel times.',
     *        '  However, program is currently limited to 3 tables.'
            stop 'abort from input1'
	  endif
          modin(i+mmax) = i
          if (vpvsm(i+mmax) .eq. 0.0) vpvsm(i+mmax) = test(1)
	  if (test(1) .lt. 0.) then
	    write(logfil, '(a)')
     *        ' test(1) (vp/vs ratio) may not be negative'
            stop 'abort from input1'
	  endif
	  if ((i .gt. 1) .and.
     *      (vpvsm(i+mmax-1) .lt. 0.)) modin(i+mmax) = 0
792	continue
      endif

      do 810 l = 1, ns
        if (mod(l) .le. kl) then
c         this is fine, unless it is an s model
          if ( mtyp(mod(l)) .eq. 's ') then
            write(punt, 795) nsta(l), mod(l)
            write(logfil, 795) nsta(l), mod(l)
795         format(' xxxerrorxxx velocity model specified for station ',
     *      a, ' was ', i4, ', an s-model, so stop!')
            stop 'abort from input1'
          endif
        else if ( ((mod(l) .gt. kl) .and. (mod(l) .lt. 11)) .or.
     *            (mod(l) .gt. mmax+3) ) then
c         model out of range, so switch to highest model
          msav = mod(l)
          mod(l) = kl
c         if the last model is an s model, use the previous one.
          if (mtyp(kl) .eq. 's ') mod(l) = kl - 1
          write(punt, 800) nsta(l), msav, mod(l)
          write(logfil, 800) nsta(l), msav, mod(l)
800       format(' xxx warning xxx velocity model specified for station'
     *    ,a, ' was too large (equaled ',i3, ').  reset to ', i3/)
        else if ( modin(mod(l)) .eq. 0 ) then
c         mod(l) is mmax+1, +2, or +3, and it doesn't exist or is not a p model
          write(punt, 801) nsta(l), mod(l)
          write(logfil, 801) nsta(l), mod(l)
801	  format(' xxx error xxx velocity model specified for station',
     *      a, '(', i3,
     *      ') was not defined or was assigned to an s table.')
          stop 'abort from input1'
	endif
810   continue
      if (nkrst .eq. 0) goto 871
c set up arrays to be used in trvdrv
      do 815 l = 1, lmax
815   jref(l) = 0
      if (iprn .le. 3) goto 825
      if (ilis.gt.0) write(punt, 820)
820   format(//1x, '   leq    m      tid(leq, m)      did(leq, m) ',
     * 38x, 'input1 formats 700 and 760.')
825   do 870 imod = 1, kl
c       loop through models, imod
        nlayers = lend(imod) - lbeg(imod)+1
        do 840 leq = lbeg(imod), lend(imod)
c         loop through eq layers, leq
          do 839 mrefr = 1, nlayers
            mref = lbeg(imod) - 1 + mrefr
            if ( (mref .gt. leq) .and. (v(mref) .gt. v(leq)) ) then
              vsqd(mrefr, leq) =
     *          sqrt((v(mref)/v(leq)-1.)*(v(mref)/v(leq)+1.))
            else
              vsqd(mrefr, leq) = 0.
            endif
            if (leq .ge. mref) then
              f(leq, mrefr) = 2.0
            else
              f(leq, mrefr) = 1.0
            endif
839       continue
840     continue
        ivl = lbeg(imod) + ivlr - 1
        do 870 leq = lbeg(imod), lend(imod)
c         loop through eq layers, leq
          leqr = leq - lbeg(imod) + 1
          if (leq .eq. ivl) then
c           preserve original thicknesses of variable layers
            sthk(imod) = thk(leq)
            sthk1(imod) = thk(leq + 1)
          endif
          do 870 mrefr = leqr, nlayers
c           loop through refracter layers, mrefr
            if (mrefr .eq. 1) goto 870
            sumt = 0.0
            sumd = 0.0
            do 850 l = lbeg(imod), lbeg(imod) + mrefr - 2
c             true if mref .le. l .or. v(mref) .le. v(l)
c             skip if refracter velocity .le. any overlaying velocity
c             if (vsqd(mrefr, l) .eq. 0.) goto 860
              if (v(lbeg(imod) + mrefr - 1) .le. v(l)) goto 860
850         continue
            jref(mrefr + lbeg(imod) - 1) = 1
            do 855 l = lbeg(imod), lbeg(imod) + mrefr - 2
c             loop from top layer to layer above refractor
              if (((l.eq.ivl) .or. (l.eq.ivl+1))
     *        .and. (ivlr .gt. 0))goto 855
              sumt = sumt + f(l, leqr)*thk(l)*vsqd(mrefr, l)
              sumd = sumd + f(l, leqr)*thk(l)/vsqd(mrefr, l)
855         continue
860         tid(leq, mrefr) = sumt*vi(mrefr + lbeg(imod) - 1)
            did(leq, mrefr) = sumd
            if (iprn .le. 3) goto 870
            if (ilis.gt.0) write(punt, 865)
     *        leq, mrefr, tid(leq, mrefr), did(leq, mrefr)
865         format(1x, 2i5, 2f15.3)
870   continue
871   do 875 l = 1, ns
875   thks(l) = blank
c set up trial hypocenter read in as test variables
      latr = 0.0
      lonr = 0.0
      if ((abs(test(3))+abs(test(4))).le.0.00001) goto 880
      la = abs(test(3))
      lo = abs(test(4))
      ala = (abs(test(3)) - la)*60.
      alo = (abs(test(4)) - lo)*60.
      ins = 'n'
      iew = 'w'
      if (test(3) .lt. 0.) ins = 's'
      if (test(4) .lt. 0.) iew = 'e'
      call fold2(latr, lonr, la, ins, ala, lo, iew, alo)
880   continue
885   continue
      return
      end
c end input1
c inside.for    []
      integer function inside(x0, y0, px, py, n)
c     check if point x0, y0 is inside polygon px(i), p(y), i = 1 to n
c     (n is the number of vertaces - eg. for a rectangle, n = 4)
c     based on godkin & pulli, bssa 74, p. 1845-1848, 1984.
c     converted to fortran from b. julian's c implementation
c     by j. c. lahr  -  18 may 1986
c     returns 0 if point outside polygon
c             1 if point at vertex
c           +-4 if point on an edge
c           +-2 if point inside
c
      logical vertex
      dimension px(n), py(n)
      inside = 0
      vertex = .false.
      x1 = px(n) - x0
      y1 = py(n) - y0
      do 20 i = 1, n
c       print *, '******************>> line segment ', i, vertex
        x2 = px(i) - x0
        y2 = py(i) - y0
        y1y2 = y1*y2
        if(y1y2 .le. 0.) then
c         print *, ' line crosses or at least touches x axis'
          ksi =  ksicr(x1, y1, x2, y2, y1y2, vertex)
c         print *, 'ksi = ', ksi
          if(iabs(ksi) .eq. 4) then
            inside = ksi
            return
          endif
          if(vertex) then
            inside = 1
            return
          endif
          inside = inside + ksi
        endif
        x1 = x2
        y1 = y2
20    continue
      return
      end
c end inside
c iprst.for    []
      subroutine iprst(l,llat,ins,alat,llon,iew,alon)
c     print out station data
      include 'params.inc'
      parameter (ndly = 11)
      character*1 ins, iew
      character*4 iahead*60, msta*5, nsta*5, icard*110
      integer punt
      common /punt/ punt
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      integer fmwt, xmwt
      common /ioxl/ fmwt(nsn), xmwt(nsn)
      common /ilpu/ sw(nsn),ndate(nsn),nhr(nsn),mno(nsn),ielv(nsn)
      common /ilpx/ calr(5, nsn),fmgc(nsn),xmgc(nsn)
      logical fsteq
      character*1 revp, exdly
      common /ip/ latr,lonr,tpdly(2,nsn),infil,iofil,indexs,
     *  iscod(nsn),ipcod(nsn), fsteq, exdly(4, nsn), revp(6, nsn)
      common /int/ thk(lmax+1),lbeg(mmax),lend(mmax),vi(lmax),vsq(lmax),
     *  vsqd(lmmax,lmax),f(lmax,lmmax),kl,ivway,sthk(mmax),sthk1(mmax),
     *  sdly(ndly,nsn)
      common /ilt/ vthk(2,nsn),ipdly(nsn),mod(nsn),ipthk(nsn)
c
      if (ilis .gt. 0)
     * write(punt,100) nsta(l), llat, ins,   alat, llon, iew,    alon,
     *  ielv(l), ipthk(l), vthk(1,l), vthk(2,l), mod(l), ipdly(l),
     * (dly(i,l), sdly(i,l),i=1,5), klas(1, l), calr(1, l), xmgc(l),
     * xmwt(l), fmwt(l), fmgc(l), ipcod(l), iscod(l), (revp(i,l),i=1,6),
     * sw(l), tpdly(1,l), (exdly(ii, l), ii = 1, 4),
     * tpdly(2,l), ndate(l), nhr(l)
  100 format(        1x,a5, 1x,i2,  a1,1x,f5.2,1x,i3,   a1, 1x,f5.2,
     *      i5,        i2,             1x,2f4.2,     i3,       i2,
     * 1x,             5(2f5.2,1x),        i2,    1x,f5.2,  1x,f4.2,
     *      i2,      i2, 1x,f4.2,             1x,2i2,      /,' * ', 6a1,
     *  3x,f4.1, 1x, f5.2,                      1x,4a1,
     *   1x,f5.2,   1x, i8,  1x, i2)
      return
      end
c end iprst
c jdate.for    [pc]
      subroutine jdate(irmo, irdy, iryr, ihr, imn, isec)
c* (vax
c      character*8 atime
c      call idate(irmo,irdy,iryr)
c      call time(atime)
c      read(atime, 10)  ihr, imn, isec
c10    format(i2, 1x, i2, 1x, i2)
c* vax)
c* (unix
c      integer time
c      character*24 ctime, nowtime
c      character*3 month(12), dnstrg, mo
c      data month/'jan', 'feb', 'mar', 'apr', 'may', 'jun',
c     *           'jul', 'aug', 'sep', 'oct', 'nov', 'dec'/
c* unix)
c* (pc
      call getdat(iryr, irmo, irdy)
      iryr = iryr - (iryr/100)*100
      call gettim(ihr, imn, isec, i100th)
c* pc)
c* (unix
c      itime = time()
c      nowtime = ctime(itime)
c      print *, nowtime
c      read(nowtime, 20) mo, irdy, ihr, imn, isec, iryr
c20    format(4x, a3, 1x, i2, 1x, i2, 1x, i2, 1x, i2, 3x, i2)
c      do 30 i = 1, 12
c        if (dnstrg(mo) .eq. month(i)) irmo = i
c30    continue
c* unix)
      return
      end
c end jdate
c ksicr.for    []
      integer function ksicr(x1, y1, x2, y2, y1y2, vertex)
c     signed crossing number - converted to fortran by j. c. lahr from
c       b. julian's c code.    18 may 1986
c
c       0                          :  does not cross -x axis
c       0 with vertex set to true  :  one end at point
c       +/-4                       :  goes through point
c       +/-2                       :  crosses -x axis
c       +/-1                       :  half crosses -x axis
c     [test for y1*y2 .le. 0. is performed in calling routine.]
c     [otherwise y1*y2 .gt. 0. would be placed here and return zero.]
      logical vertex
      t = x1*y2 - x2*y1
c     print *, 't = ', t
      if( (t .eq. 0.) .and. (x1*x2 .le. 0.) ) then
c       print *, 'line passes through or to point'
        if(y1y2 .ne. 0.) then
c         print *, 'neither end of line terminates at point'
          if(y2 .gt. 0.) then
c           print *, 'line passes up through point'
            ksicr = 4
          else
c           print *, 'line passes down through point'
            ksicr = -4
          endif
          return
        else if(x1*x2 .eq. 0.) then
c         print *, 'y1*y2 = 0. and x1*x2 = 0.'
          vertex = .true.
c         print *, 'vertex = ', vertex
        else
          if(x1 .lt. x2) then
c           print *, 'line passes to right through point'
            ksicr = 4
          else
c           print *, 'line passes to left through point'
            ksicr = -4
          endif
          return
        endif
      else if(y1y2 .lt. 0.) then
c       print *, 'complete crossing of x axis'
        if(t*y2 .lt. 0) then
c         print *, 'complete crossng of -x asix'
          if(y2 .gt. 0.) then
            ksicr = 2
          else
            ksicr = -2
          endif
          return
        endif
      else if(y1 .eq. 0.) then
c       print *, 'half crossing, y1 equals 0'
        if((x1 .lt. 0) .and. (y2 .ne. 0.) ) then
          if(y2 .gt. 0.) then
            ksicr = 1
          else
            ksicr = -1
          endif
          return
        endif
      else
c       print *, 'half crossing, y2 must equal 0'
        if(x2 .lt. 0) then
          if(y1 .lt. 0) then
            ksicr = 1
          else
            ksicr = -1
          endif
          return
        endif
      endif
      ksicr = 0
      return
      end
c end ksicr