c lentru.for    []
      integer function lentru(alph)
c     finds the true length of a character variable
      character alph*(*)
      l = len(alph)
      do 100 i = l, 1, -1
       if(alph(i:i).ne.' ' .and. alph(i:i).ne.'\0') then
         lentru = i
         return
        endif
100   continue
      lentru = 0
      return
      end
c end lentru
c line3.for    []
      subroutine line3(nr, p, s, pwt, swt, test,
     *                         vpvs1, sint1, orig1, se1, ses, oses,
     *                         vpvs2, sint2, orig2, se2, sep, osep,
     *                         vpvs3, sint3, orig3, se3)
c bi-weighted regression:
c assume p and s both have errors.  regress s on p.  (madansky, 1959)
c
      include 'params.inc'
c calling arguments:
      integer           nr
c                                 ! number of arrival time pairs
      real              p(nr)
c                                 ! p arrival times
      real              s(nr)
c                                 ! s arrival times
      real              pwt(nr)
c                                 ! p weights
      real              swt(nr)
c                                 ! s weights
      real              vpvs1
c                                 ! [output] vpvs ratio   (s vs p)
      real              vpvs2
c                                 ! [output] vpvs ratio   (p vs s)
      real              vpvs3
c                                 ! [output] vpvs ratio   (bi-weight)
      real              sint1
c                                 ! [output] s intercept  (s vs p)
      real              sint2
c                                 ! [output] s intercept  (p vs s)
      real              sint3
c                                 ! [output] s intercept  (bi-weight)
      real              orig1
c                                 ! [output] orig time  (s vs p)
      real              orig2
c                                 ! [output] orig time  (p vs s)
      real              orig3
c                                 ! [output] orig time  (bi-weight)
      real              oses
c                                 ! rms difference between s obs & comp
      real              osep
c                                 ! rms difference between p obs & comp
      real              se1
c                                 ! [output]standard error (s vs p)
      real              se2
c                                 ! [output]standard error (p vs s)
      real              se3
c                                 ! [output]standard error (bi-weight)
c
c                                 ! sep & ses are not allowd to go below test(29
      real              sep
c                                 ! [output]s.e. of p used in s.e. of vp/vs (p v
      real              ses
c                                 ! [output]s.e. of s used in s.e. of vp/vs (s v
      real              test(50)
c                                 ! test parameters
c
      double precision  dubl
c                                 ! statement function
      double precision  arg
c                                 ! variance
      double precision  dslop1
c                                 ! vpvs ratio (s vs p)
      double precision  dslop2
c                                 ! vpvs ratio (p vs s)
      double precision  eqwt(npa)
c                                 ! equation weight
      double precision  dorig1
c                                 ! origin (s vs p)
      double precision  dorig2
c                                 ! origin (p vs s)
      double precision  dorig3
c                                 ! origin (bi-weight)
      double precision  dp(npa)
c                                 ! p
      double precision  ds(npa)
c                                 ! s
      double precision  dse1
c                                 ! se (s vs p)
      double precision  dse2
c                                 ! se (p vs s)
      double precision  dse3
c                                 ! se (bi-weight)
      double precision  dsint1
c                                 ! sinter (s vs p)
      double precision  dsint2
c                                 ! sinter (p vs s)
      double precision  dsint3
c                                 ! sinter (bi-weight)
      double precision  dvpvs
c                                 ! increment in vpvs
      double precision  pav
c                                 ! bi-weighted average of p
      double precision  pavp
c                                 ! average value of p weighted with p weights
      double precision  pavs
c                                 ! average value of p weighted with s weights
      double precision  pssum
c                                 ! sum weighted p*s
      double precision  psum
c                                 ! sum p times
      double precision  p2sum
c                                 ! sum weighted p**2
      double precision  sav
c                                 ! bi-weighted average of s
      double precision  savp
c                                 ! average value of s weighted with p weights
      double precision  savs
c                                 ! average value of s weighted with s weights
      double precision  ssum
c                                 ! sum s times
      double precision  s2sum
c                                 ! sum weighted s**2
      double precision  tsum
c                                 ! trial sum to be minimized
      double precision  tsumin
c                                 ! minimum sum
      double precision  vpvsp
c                                 ! current preferred value of vpvs
      double precision  vpvsp1
c                                 ! next preferred value of vpvs
      double precision  vpvst
c                                 ! trial lower limit of vpvs
      double precision  wpsum
c                                 ! bi-weighted sum of p
      double precision  wpsump
c                                 ! sum of p times weighted with p weights
      double precision  wpsums
c                                 ! sum of p times weighted with s weights
      double precision  wssum
c                                 ! bi-weighed sum of s
      double precision  wssump
c                                 ! sum of s times weighted with p weights
      double precision  wssums
c                                 ! sum of s times weighted with s weights
      double precision  wtsum
c                                 ! sum of bi-weights
      double precision  wtsump
c                                 ! sum of p weights
      double precision  wtsums
c                                 ! sum of s weights
c
c initialize some variables
      psum = 0.d0
      ssum = 0.d0
      wpsum = 0.d0
      wpsump = 0.d0
      wpsums = 0.d0
      wssump = 0.d0
      wssums = 0.d0
      wssum = 0.d0
      p2sum = 0.d0
      s2sum = 0.d0
      pssum = 0.d0
      wtsum = 0.d0
      wtsump = 0.d0
      wtsums = 0.d0
c
c main loop to compute vp/vs and origin time
      do 40 i=1,nr
        pi = p(i)
        si = s(i)
        dp(i) = dubl(pi)
        ds(i) = dubl(si)
        wpsump = wpsump + pwt(i)*dp(i)
        wpsums = wpsums + swt(i)*dp(i)
        wtsump = wtsump + pwt(i)
        wssump = wssump + pwt(i)*ds(i)
        wssums = wssums + swt(i)*ds(i)
        wtsums = wtsums + swt(i)
40    continue
c compute weighted and unweighted averages
      pavp = wpsump/wtsump
      pavs = wpsums/wtsums
      savp = wssump/wtsump
      savs = wssums/wtsums
c
c compute slop (vp/vs) based on s dependent on p
      do 50 i=1,nr
        p2sum = p2sum + swt(i)*(dp(i)-pavs)**2
        pssum = pssum + swt(i)*(dp(i)-pavs)*(ds(i)-savs)
        s2sum = s2sum + swt(i)*(ds(i)-savs)**2
50    continue
cd    print *, 's2sum, savs, pssum, pavs, p2sum'
cd    print *,  s2sum, savs, pssum, pavs, p2sum
c
      dslop1 = pssum/p2sum
c compute s.e. of s readings with weight code of zero (ses)
      arg = ( s2sum - dslop1*pssum ) / (nr - 2.d0)
      if (arg .lt. 0.d0) then
        print *, 'probable round off error, arg = ', arg
        print *, ' but arg should not be less than 0.'
        arg = 0.d0
      endif
      oses = dsqrt(arg)
c use the larger estimate of ses
      ses = oses
      if (ses .lt. abs(test(29))) ses = abs(test(29))
      dse1 = ses/dsqrt(p2sum)
cd    print *, 'ses, test(29) ', ses, test(29)
c
      dsint1 = savs - dslop1*pavs
      if(dslop1 .eq. 1.0) dslop1 = 1.01
      dorig1 = (dslop1*pavs-savs)/(dslop1-1.0)
cd     write(7, 51) dslop1, dse1
cd51    format(' for s vs p regression, vpvs = ', d10.4, ' and se = ',
cd   *   d10.4)
c
c compute slope (vp/vs) based on p dependent on s
      p2sum = 0.d0
      pssum = 0.d0
      s2sum = 0.d0
      do 52 i=1,nr
        p2sum = p2sum + pwt(i)*(dp(i)-pavp)**2
        pssum = pssum + pwt(i)*(dp(i)-pavp)*(ds(i)-savp)
        s2sum = s2sum + pwt(i)*(ds(i)-savp)**2
52    continue
cd    print *, 's2sum, savp, pssum, pavp, p2sum'
cd    print *,  s2sum, savp, pssum, pavp, p2sum
c
      dslop2 = pssum/s2sum
c compute s.e. of p readings with weight code of zero (sep)
      arg = ( p2sum - dslop2*pssum ) / (nr - 2.d0)
      if (arg .lt. 0.d0) then
        print *, 'probable round off error, arg = ', arg
        print *, ' but arg should not be less than 0.'
        arg = 0.d0
      endif
      osep = dsqrt(arg)
      sep = osep
      if (sep .lt. abs(test(29))) sep = abs(test(29))
      dse2 = sep/( dsqrt(s2sum)*(dslop2**2.d0) )
      dslop2 = 1.d0/dslop2
      dsint2 = savp - dslop2*pavp
      if(dslop2 .eq. 1.0) dslop2 = 1.01
      dorig2 = (dslop2*pavp-savp)/(dslop2-1.0)
cd     write(7, 53) dslop2, dse2
cd53    format(' for p vs s regression, vpvs = ', d10.4, ' and se = ',
cd   *   d10.4)
c
c find minimum sum
      dvpvs = .6
c     initial central value
      vpvsp = (dslop1 + dslop2)/2.d0
      tsumin = 99999.
c the central value need only be computed the first time
      do 60 j = 1, 8
        do 54 i = 1, 5
          vpvst = vpvsp - dvpvs*(-3 + i)
          tsum = 0.d0
c         compute averages
          do 43 k = 1, nr
            eqwt(k) = 1./( 1./swt(k) + (vpvst**2.d0)/pwt(k) )
            wpsum = wpsum + eqwt(k)*dp(k)
            wssum = wssum + eqwt(k)*ds(k)
            wtsum = wtsum + eqwt(k)
43        continue
          pav = wpsum/wtsum
          sav = wssum/wtsum
          do 44 k = 1, nr
            tsum = tsum + eqwt(k)*
     *                ( ds(k)-sav - vpvst*(dp(k)-pav) )**2
44        continue
          if (tsum .lt. tsumin) then
            tsumin = tsum
            vpvsp1 = vpvst
          endif
cd        print *, vpvst, tsum, tsumin
54      continue
cd     write(injump, 55) vpvsp1, tsum
cd55    format(d12.4, 3x, d12.4)
      vpvsp = vpvsp1
      dvpvs = dvpvs*.4
60    continue
c
      dse3 = dsqrt(dse1**2.d0 + dse2**2.d0)
      dsint3 = sav - vpvsp1*pav
      if(vpvsp1 .eq. 1.0) vpvsp1 = 1.01
      dorig3 = (vpvsp1*pav-sav)/(vpvsp1-1.0)
c convert to real for output
      vpvs1 = dslop1
      sint1 = dsint1
      se1 = dse1
      orig1 = dorig1
c
      vpvs2 = dslop2
      sint2 = dsint2
      se2 = dse2
      orig2 = dorig2
c
      vpvs3 = vpvsp1
      sint3 = dsint3
      se3 = dse3
      orig3 = dorig3
c
      return
      end
c end line3
c linv.for    []
c test driver for subroutine linv
c      data elevmx/3.0/, ielest/2000/, x/10./, zeq/3.0/
c      data v/3.0/, ak/1.0/
c      print *, 'welcome to linv!  this program tests subroutine'
c      print *, 'linv, which computes travel times in a half space'
c      print *, 'with linearly increasing velocity.'
c
c20    elevmx = raskk(
c     *'maximum topographic elevation (reference elevation in km)',
c     *          elevmx)
c      ielest = iaskk(
c     *'station elevation above sea level', ielest)
c      stz = elevmx - ielest*.001
c      v = raskk('velocity at reference elevation (km/s)', v)
c      ak = raskk('velocity gradient (km/s per km)', ak)
c30    x = raskk('epicentral distance (km)', x)
c      zeq = raskk('depth of eq beneath reference elevation (km)', zeq)
c
c      call linv(x, zeq, v, ak, t, ain, dtdd, dtdh, stz)
c      vsta = v + ak*(elevmx - ielest/1000.)
c      zeff = zeq - (elevmx - ielest/1000.)
c      call linvjl(x, zeff, vsta, ak, t, ain, dtdd, dtdh, stz,
c    *  vst, veq)
c
c      print *, 'linvjl results to test linv when eq is deeper',
c     * ' than station:'
c      print *, 'travel time (s) = ', t
c      ain =  asin(c      if(ain .lt. 0.) ain = 180. + ain
c      ain = 180. - ain
c      print *, 'ain in degrees = ', ain
c      print *, 'partial of tt wrt dist and depth = ', dtdd, dtdh
c
c      if(x .gt. 0.) goto 20
c      stop
c      end
      subroutine linv(x, zeq, v, ak, t, ain, dtdd, dtdh, stz,
     *  vst, veq)
c compute travel time and partial derivatives for a linear increasing
c velocity model in a region of significant topography.
      real x
c          x      epicentral distance in km
      real zeq
c          zeq    depth in km of earthquake beneath reference elevation
      real v
c          v      velocity at reference elevation in km/sec
      real ak
c          ak     velocity gradient km/sec per km
      real t
c          t      computed travel time in sec
      real ain
c          ain    angle-of-incidence
      real dtdd
c          dtdd   partial derivative of tt wrt distance
      real dtdh
c          dtdh   partial derivative of tt wrt depth
      real stz
c          stz    depth in km of station beneath reference elevation
      real vst
c          vst    velocity at station
      real veq
c          veq    velocity at earthquake
      real dz
c          dz     depth of eq beneath station
      real dc
c          dc     distance in km from eq or station, whichever is
c                 at higher elevation, to center of raypath circle
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      parameter (deg = 1.0/rad)
c      print *, 'station elevation (m)    = ', ielest
c      print *, 'eq depth (km)            = ', zeq
c      print *, 'epicentral distance      = ', x
c      print *, 'velocity at ref elev and grad = ', v, ak
      vst = v + ak*stz
      veq = v + ak*zeq
      dz = (zeq - stz)
      dist = sqrt(dz*dz + x*x)
      if (dist .lt. .0001) then
c       distance less than 10 cm
        t = 0.0
        dtdd = 1./veq
        dtdh = 1./veq
        anin = 1.
        ain = 90.
        return
      endif
      p = ak*ak*(x*x + dz*dz)/(2.*vst*veq) + 1.
      t = coshi(p)/ak
      if (dz .lt. 0) then
c earthqauke is above station (ray is down-going)
c        print *, 'eq is above station, dz = ', dz
        if (x .eq. 0) then
          anin = 0.0
          ain = 0.0
        else
          dc = (x*x + 2.*abs(dz)*veq/ak + dz*dz)/(2.*x)
          anin = (veq/ak) / sqrt(dc*dc + (veq/ak)**2)
          ain = deg*asin(anin)
c          print *, 'anin = ', anin
        endif
      else
c station is above earthquake (ray may be up- or down-going)
c        print *, 'station is above eq, dz = ', dz
        alam = vst*(cosh(ak*t) - 1.)/ak
        if(ak .le. .0001) alam = 0.0
        dif = dz - alam
        anin = x/sqrt(dif*dif + x*x)
      endif
      dtdd = x*ak/( vst*(vst + ak*dz)*sqrt(p*p - 1.) )
      dtdh = (dz*ak/(vst*(vst+ak*dz)) -
     *         ak*ak*(x*x + dz*dz)/(2.*vst*((vst + ak*dz)**2.)))/
     *         ((p*p - 1.)**0.5)
      if(dtdh .gt. 0.0) then
c       upgoing
        ain = 180. - deg*asin(anin)
      else
c       downgoing
        ain = deg*asin(anin)
      endif
      return
      end
c end linv
c linvol.for    []
      subroutine linvol(delta, depthsv, stzsv, ak, voa, aha, vi, t,
     *  ain, dtdd, dtdh)
c model:  linearly increasing velocity over a halfspace
c bill gawthrop  (modified by j.c. lahr  3/6/91)
c     delta  epicentral distance (km)
c     depthsv eq depth (km)
c     stzsv   station depth (km)
c     ak     velocity gradient
c     voa    v at surface
c     aha    thickness of linear layer (depth to half space)
c     vi     velocity of half space
c     t      travel time
c     ain    angle of incidence
c     dtdd   partial derivative of tt wrt distance
c     dtdh   partial derivative of tt wrt depth
      integer punt
      common /punt/ punt
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      parameter (deg = 1.0/rad)
      logical flip
      dist = sqrt(delta**2 + (depthsv - stzsv)**2)
      if (dist .lt. .0001) then
c       distance less than 10 cm
        t = 0.0
        dtdd = 1./(voa + ak*depthsv)
        dtdh = 1./(voa + ak*depthsv)
        ain = 90.
        return
      endif
      if((stzsv .ge. aha - .0001) .and.
     *  (depthsv .ge. aha - .0001)) then
c both eq and station are within (or within 10 cm) of the
c constant velocity halfspace
c       write(punt, *) 'eq and station both in halfspace'
        t = dist/vi
        dtdd = delta/(dist*vi)
        dtdh = (depthsv - stzsv)/(dist*vi)
        if (stzsv .eq. depthsv) then
          ain = 90.0
        else
          if(stzsv .lt. depthsv) then
c           upward
            ain = 180 + deg*atan(delta/(stzsv - depthsv))
          else
c           downward
            ain = deg*atan(delta/(stzsv - depthsv))
          endif
        endif
        return
      endif
      if(depthsv .ge. stzsv) then
        flip = .false.
        deptha = depthsv
        stza = stzsv
      else
        flip = .true.
        deptha = stzsv
        stza = depthsv
      endif
c if near boundary, then move to boundary
      if(abs(deptha - aha) .le. .0001) deptha = aha
      if(abs(stza - aha) .le. .0001) stza = aha
c set vsource
      if(flip) then
        if(stza .le. aha) then
          vsource = voa + stza*ak
        else
          vsource = vi
        endif
      else
        if(deptha .le. aha) then
	  vsource = voa + deptha*ak
	else
	  vsource = vi
	endif
      endif
c create modified model with station at surface
      ah = aha - stza
      vo = voa + ak*stza
      depth = deptha - stza
      stz = 0.0
      if (depth .le. ah) then
c source: within the layer
c       write(punt, *) ' source: within the layer'
        vz = vo + depth*ak
c       compute time of direct wave
        call tatime(delta,depth,vo,ak,ah,tat,aina)
c       compute time of refracted wave
        call tbtime(delta,depth,vo,ak,ah,vi,tbt,ainb)
        if (tbt .lt. tat) then
c         refracted wave is faster
          t = tbt
          if(flip) then
c           downward ray
            anin = vsource*sin(ainb)/vz
            ain = asin(anin)
          else
c           downward ray
            anin = sin(ainb)
            ain = ainb
          endif
        else if (tat .eq. 900000.) then
c         neither direct nor refracted wave could be computed!
          print *,
     *      'neither direct nor refracted wave could be computed!'
          print *, 'so stop'
	  print *, 'delta,depth,vo,ak,ah,vi'
	  print *, delta,depth,vo,ak,ah,vi
          stop '1: abort from linvol'
        else
c         direct wave is faster
          t = tat
          if(flip) then
c           downward ray
            anin = vsource*sin(aina)/vz
            ain = asin(anin)
          else
c           upward or downward ray
            anin = sin(aina)
            ain = aina
          endif
        endif
      else
c source: within the halfspace
c       write(punt, *) ' source: within the halfspace'
c	write(punt, *) 'delta,depth,vo,ak,ah,vi'
c	write(punt, *)  delta,depth,vo,ak,ah,vi
        call tdtime(delta,depth,vo,ak,ah,vi,td,aina)
        if (td.eq.900000.) then
c         direct wave could not be computed!
          print *,
     *      'direct wave could not be computed, so stop!'
          stop '2: abort from linvol'
        else
c         direct wave from halfspace
          t = td
          if(flip) then
c           downward ray
            anin = vsource*sin(aina)/vi
            ain = asin(anin)
          else
c           upward ray
            anin = sin(aina)
            ain = aina
          endif
        endif
      endif
      dtdd = anin/vsource
      dtdh = -cos(ain)/vsource
      ain = deg*ain
      return
      end
c end linvol
c lisinc.for    []
      subroutine lisinc(qno, sumrms, rms, jav, iq, nrwt, avwt, nr,
     * kdx, wt, nrp, xmag, fmag, blank, avxm, avfm, nxm, nfm, sxm,
     * sxmsq, sfm, sfmsq, sw, ipcod, ww, ksmp, nsmpr, ssmpr, x,
     * ssmpq, ssmpw, nres, sr, srsq, srwt, nresp, srp, srpsq, srpwt,
     * iscod, nsres, ssr, ssrsq, ssrwt, nsresp, ssrp, ssrpsq, ssrpwt)
c increment summary after each earthquake
      save
      include 'params.inc'
      dimension kdx(npa), wt(npa), xmag(npa), fmag(npa), sw(nsn)
      dimension ipcod(nsn), ww(npa), ksmp(npa), x(4,npa), nres(nsn)
      dimension iscod(nsn)
      dimension nxm(nsn), nfm(nsn), sr(nsn), srsq(nsn),
     * srwt(nsn),sxm(nsn),sxmsq(nsn),sfm(nsn),sfmsq(nsn),qno(4)
      dimension nresp(nsn), srp(nsn), srpsq(nsn), srpwt(nsn)
      dimension nsresp(nsn), ssrp(nsn), ssrpsq(nsn), ssrpwt(nsn)
      dimension nsres(nsn),ssr(nsn),ssrsq(nsn),ssrwt(nsn),nsmpr(nsn),
     * ssmpr(nsn),ssmpq(nsn),ssmpw(nsn)
c
      qno(jav) = qno(jav) + 1
      sumrms = sumrms + rms
      if(jav .gt. iq) return
c compute the sqrt(sum of weights squared) for pavlis average residual calc
      srswq = 0.0
      if(nrwt .gt. 4) srswq = avwt*(nrwt - 4)
c accumulate summary of travel time residuals --- include p, s, and s-p
      do 825 i = 1, nr
c ji is index to station list for phase number i
        ji = kdx(i)
        wtu = wt(i)
        if(i .le. nrp) then
          if((xmag(i) .ne. blank) .and. (avxm .ne. blank)) then
            dxmag = xmag(i)-avxm
            nxm(ji) = nxm(ji)+1
            sxm(ji) = sxm(ji)+dxmag
            sxmsq(ji) = sxmsq(ji)+dxmag**2
          endif
          if((fmag(i) .ne. blank) .and. (avfm .ne. blank)) then
            dfmag = fmag(i)-avfm
            nfm(ji) = nfm(ji)+1
            sfm(ji) = sfm(ji)+dfmag
            sfmsq(ji) = sfmsq(ji)+dfmag**2
          endif
          if( (sw(ji) .eq. 0.) .or. ((ipcod(ji) .gt. 3) .and.
     *    (ipcod(ji) .lt. 9)) ) wtu = ww(i)
          if(wtu .eq. 0.) goto 825
          if(ksmp(i) .eq. 0) then
c smp residual calculation
            nsmpr(ji) = nsmpr(ji) + 1
            ssmpr(ji) = ssmpr(ji) + x(4,i)*wtu
            ssmpq(ji) = ssmpq(ji) + x(4,i)**2*wtu
            ssmpw(ji) = ssmpw(ji) + wtu
            goto 825
          else
c p-residual calculation
            nres(ji) = nres(ji) + 1
            sr(ji) = sr(ji) + x(4,i)*wtu
            srsq(ji) = srsq(ji) + x(4,i)**2*wtu
            srwt(ji) = srwt(ji) + wtu
c pavlis p-residual calculation
            if(srswq .gt. 0.0) then
              wtu = wtu*srswq
              nresp(ji) = nresp(ji) + 1
              srp(ji) = srp(ji) + x(4,i)*wtu
              srpsq(ji) = srpsq(ji) + x(4,i)**2*wtu
              srpwt(ji) = srpwt(ji) + wtu
            endif
            goto 825
          endif
        else
c s-residual calculation
          if( (sw(ji) .eq. 0.) .or. ((iscod(ji) .gt. 3) .and.
     *    (iscod(ji) .lt. 9)) ) wtu = ww(i)
          if(wtu .eq. 0.0) goto 825
          nsres(ji) = nsres(ji) + 1
          ssr(ji) = ssr(ji) + x(4,i)*wtu
          ssrsq(ji) = ssrsq(ji) + x(4,i)**2*wtu
          ssrwt(ji) = ssrwt(ji) + wtu
c pavlis s-residual calculation
          if(srswq .gt. 0.0) then
            wtu = wtu*srswq
            nsresp(ji) = nsresp(ji) + 1
            ssrp(ji) = ssrp(ji) + x(4,i)*wtu
            ssrpsq(ji) = ssrpsq(ji) + x(4,i)**2*wtu
            ssrpwt(ji) = ssrpwt(ji) + wtu
          endif
        endif
  825 continue
      return
      end
c end lisinc
c lissum.for    []
      subroutine lissum(lisco)
c output summary of time and magnitude residuals
      save
      include 'params.inc'
      parameter (ndly = 11)
      real lat,lon
      character fixsw*1, fixp*1, fixs*1, jfmt*4, line*80
      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 /dhil/ iq,ilat,kms,icat
      common /dmost/ ipun,ivlr,blank
      common /hl/ nwad, tslope, tsqsl
      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 /ilotv/ elvdly(npa)
      common /ilpx/ calr(5, nsn),fmgc(nsn),xmgc(nsn)
      common /ilv/ c(nsn),e(nsn)
      common /ilt/ vthk(2,nsn),ipdly(nsn),mod(nsn),ipthk(nsn)
      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)
      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 /lc/ nres(nsn)
      common /ohbl/ jav
      common /olnx/ xmag(npa),fmag(npa),avxm,avfm,xmmed,fmmed
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /pnl/ ww(npa)
      common /pnoqtx/ ldx(npa)
      common /qmost/ wt(npa),z
      common /reloc/ irelo, nreloc
      common /rq/ xmean(4),avwt,aveaj
      common /tmost/ x(4,npa)
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      dimension nxm(nsn), nfm(nsn), sr(nsn), srsq(nsn),
     * srwt(nsn),sxm(nsn),sxmsq(nsn),sfm(nsn),sfmsq(nsn),qno(4)
      dimension nresp(nsn), srp(nsn), srpsq(nsn), srpwt(nsn)
      dimension nsresp(nsn), ssrp(nsn), ssrpsq(nsn), ssrpwt(nsn)
      dimension nsres(nsn),ssr(nsn),ssrsq(nsn),ssrwt(nsn),nsmpr(nsn),
     * ssmpr(nsn),ssmpq(nsn),ssmpw(nsn),qnoper(4)
      dimension avres(7),sdres(7)
      data avres/7*0./
      if (iprn .ge. 5) write(punt, '(a, i5)') ' begin lissum ', lisco
      logu = 6
c initialize arrays
      if(lisco .eq. 1) then
        nwad = 0
        tslope = 0.
        tsqsl = 0.
        sumrms = 0.0
        do 48 l=1,ns
          nres(l)=0
          nresp(l)=0
          nres(l)=0
          nxm(l)=0
          nfm(l)=0
          sr(l)=0.
          srp(l)=0.
          srsq(l)=0.
          srpsq(l)=0.
          srwt(l)=0.
          srpwt(l)=0.
          nsres(l) = 0
          nsresp(l) = 0
          ssr(l) = 0.0
          ssrp(l) = 0.0
          ssrsq(l) = 0.0
          ssrpsq(l) = 0.0
          ssrwt(l) = 0.0
          ssrpwt(l) = 0.0
          nsmpr(l) = 0
          ssmpr(l) = 0.0
          ssmpq(l) = 0.0
          ssmpw(l) = 0.0
          sxm(l)=0.
          sxmsq(l)=0.
          sfm(l)=0.
          sfmsq(l)=0.
   48   continue
        do 49 i=1,4
          qno(i)=0.
   49   continue
        return
      endif
c
    4 if(lisco .eq. 0) then
        call lisinc(qno, sumrms, rms, jav, iq, nrwt, avwt, nr,
     *  kdx, wt, nrp, xmag, fmag, blank, avxm, avfm, nxm, nfm, sxm,
     *  sxmsq, sfm, sfmsq, sw, ipcod, ww, ksmp, nsmpr, ssmpr, x,
     *  ssmpq, ssmpw, nres, sr, srsq, srwt, nresp, srp, srpsq, srpwt,
     *  iscod, nsres, ssr, ssrsq, ssrwt, nsresp, ssrp, ssrpsq, ssrpwt)
        return
      endif
c
c write summary of time and mag residuals (statements 0 to 70)
      rewind 2
      qsum=qno(1)+qno(2)+qno(3)+qno(4)
      if(qsum .eq. 0.) then
        write(logu, 3)
        if(punt .ne. logu) write(punt, 3)
    3   format(' there were no events located in this run.')
        return
      endif
      averms = 0.0
      if(qsum .ne. 0.0) averms = sumrms/qsum
      avpvs = 0.
      if(nwad .gt. 0) avpvs = tslope/nwad
      sdpvs = 0.
      if(nwad .gt. 1) then
        arg = tsqsl/(nwad-1) - tslope*tslope/(nwad*(nwad-1))
        if(arg .lt. 0.0) arg = 0.0
        sdpvs = sqrt(arg)
      endif
      write(logu,5) averms, avpvs, nwad, sdpvs, (qno(i),i=1,4), qsum
      if(punt .ne. logu)
     *write(punt,5) averms, avpvs, nwad, sdpvs, (qno(i),i=1,4), qsum
    5 format(' average rms of all events = ', f10.5, /,
     *       ' average vp/vs ratio = ', f10.2, ' for ', i4, ' events.',
     *       '  standard deviation of ratio = ', f10.2,
     * //, ' ***** class/     a     b     c     d total *****',
     * //, '      number/', 5f6.1)
      do 10 i=1,4
   10 qnoper(i)=100.*qno(i)/qsum
      write(logu,15)(qnoper(i),i=1,4)
      if(punt .ne. logu) write(punt,15)(qnoper(i),i=1,4)
   15 format(/2x,'percentage/',4f6.1)
      if(iqcls .eq. '0') return
      write(logu,20) iqcls
      if(punt .ne. logu) write(punt,20) iqcls
   20 format(/,'   include only class ',a1,' and better in the',
     *  ' following statistics.' /)
   71 write(logu,1020)
      if(punt .ne. logu) write(punt,1020)
 1020 format(
     *'          ----------- p residuals --------------',
     *'    ----------- s residuals -------------', /,

     *'           no event wting        event wting',
     *'         no event wting        event wting', /,

     *' station  n  wt  ave   sd     n   wt  ave   sd',
     *'      n  wt   ave  sd     n   wt   ave  sd    station')

      do 70 i=1,ns

        do 30 j=1,7
          avres(j)=0.
          sdres(j)=0.
   30   continue

        avwt1 = 0.
        avwt6 = 0.
        if(nres(i) .gt. 0) then
c compute average p residual
          avres(1)=sr(i)/srwt(i)
c compute average weight
          avwt1 = srwt(i)/nres(i)
          if(nres(i) .gt. 1) then
	    temp = srsq(i)/srwt(i) - avres(1)**2
	    if(temp .gt. 0.0) then
              sdres(1)=sqrt(temp)
	    else
	      sdres(1)=0.0
	    endif
          endif
          if(nresp(i) .gt. 0) then
            if(srpwt(i) .gt. 0.0) then
c compute average residual, weight, and std. deviation a la pavlis
c see pavlis and hokanson (1985) - jgr, v 90, p 12777-12789
              avres(6)=srp(i)/srpwt(i)
              avwt6 = srpwt(i)/nresp(i)
	      temp = srpsq(i)/srpwt(i) - avres(6)**2
	      if(temp .gt. 0.0) then
                sdres(6)=sqrt(temp)
	      else
	        sdres(6)=0.0
	      endif
            endif
          endif
        endif

        avwt2 = 0.
        avwt7 = 0.
        if(nsres(i) .gt. 0) then
          avwt2 = ssrwt(i) / nsres(i)
          avres(2) = ssr(i)/ssrwt(i)

          if(nsres(i) .gt. 1) then
	    temp = ssrsq(i)/ssrwt(i) - avres(2)**2
	    if(temp .gt. 0.0) then
              sdres(2) = sqrt(temp)
	    else
	      sdres(2) = 0.0
	    endif
          endif

          if(nsresp(i) .gt. 0) then
            if(ssrpwt(i) .gt. 0.0) then
c compute average residual, weight, and std. deviation a la pavlis
              avres(7) = ssrp(i)/ssrpwt(i)
              avwt7 = ssrpwt(i) / nsresp(i)
	      temp = ssrpsq(i)/ssrpwt(i) - avres(7)**2
	      if(temp .gt. 0.0) then
                sdres(7) = sqrt(temp)
	      else
	        sdres(7) = 0.0
	      endif
            endif
          endif
        endif

        ipt = nres(i) + nsres(i)
        if(ipt .eq. 0) goto 70

        fixsw = ' '
        if (sw(i) .eq. 0.) fixsw = 'w'
        fixp = ' '
        if ( (ipcod(i) .gt. 3) .and.
     *  (ipcod(i) .lt. 9) ) fixp = 'p'
        fixs = ' '
        if ( (iscod(i) .gt. 3) .and.
     *  (iscod(i) .lt. 9) ) fixs = 's'

        write(logu, 65) fixsw, nsta(i),
     *  nres(i), avwt1, avres(1), sdres(1),
     *  nresp(i), avwt6, avres(6), sdres(6), fixp,
     *  nsres(i), avwt2, avres(2), sdres(2),
     *  nsresp(i), avwt7, avres(7), sdres(7), fixs, nsta(i)

        if(punt .ne. logu) write(punt, 65) fixsw, nsta(i),
     *  nres(i), avwt1, avres(1), sdres(1),
     *  nresp(i), avwt6, avres(6), sdres(6), fixp,
     *  nsres(i), avwt2, avres(2), sdres(2),
     *  nsresp(i), avwt7, avres(7), sdres(7), fixs, nsta(i)

   65   format(1x, a1, 1x, a5,
     *  2(i4, f4.1, 2f6.3, i4, f5.1, 2f6.3, a1), 2x, a4)

      if(irelo .gt. 0) then
        if(nreloc .eq. irelo) then
c write out revised primary station record
          call unfold2(lat(i),lon(i),la,ins,ala,lo,iew,alo)
          write(line, 67) nsta(i), la, ins,  ala, lo, iew,  alo,
     *    ielv(i), mod(i), ipthk(i)
67        format(              a5, i2,  a1, f5.2, i4,  a1, f5.2,
     *         i5,     i2,       i1)
c is a5 correct?  Shouldn't it be a4????????????

          call riorbk(vthk(1, i), ivt, jfmt, 4, 2)
          write(line(31:34), jfmt) ivt
          call riorbk(vthk(2, i), ivt, jfmt, 4, 2)
          write(line(35:38), jfmt) ivt
          write(line(39:39), '(i1)') ipdly(i)
          call riorbk(dly(1, i), idly, jfmt, 4, 2)
          write(line(40:43), jfmt) idly
          call riorbk(sdly(1, i), idly, jfmt, 4, 2)
          write(line(44:47), jfmt) idly
          write(13, '(a)') line
        else
          dly(1, i) = dly(1, i) + avres(6)
          sdly(1, i) = sdly(1, i) + avres(7)
        endif
      endif
   70 continue

      write(logu,1021)
      if(punt .ne. logu) write(punt,1021)
 1021 format(/,
     *'           s-p residuals        x-mag res       f-mag res', /
     *' station  n  wt  ave   sd     n  ave   sd     n  ave   sd')

      do 75 i=1,ns

        do 73 j=1,7
          avres(j)=0.
          sdres(j)=0.
   73   continue

        if(nxm(i) .gt. 0) then
          avres(3)=sxm(i)/nxm(i)
          if(nxm(i) .gt. 1) then
	    temp = sxmsq(i)/nxm(i) - avres(3)**2
	    if(temp .gt. 0.0) then
              sdres(3)=sqrt(temp)
	    else
	      sdres(3) = 0.0
 	    endif
          endif
        endif

        if(nfm(i) .gt. 0) then
          avres(4)=sfm(i)/nfm(i)
          if(nfm(i) .gt. 1) then
	    temp = sfmsq(i)/nfm(i) - avres(4)**2
	    if(temp .gt. 0.0) then
              sdres(4)=sqrt(temp)
	    else
	      sdres(4) = 0.0
	    endif
          endif
        endif

        asmpw = 0.
        if(nsmpr(i) .gt. 0) then
          asmpw = ssmpw(i) / nsmpr(i)
          avres(5) = ssmpr(i)/ssmpw(i)
          if(nsmpr(i) .gt. 1) then
	    temp = ssmpq(i)/ssmpw(i) - avres(5)**2
	    if(temp .gt. 0.0) then
              sdres(5) = sqrt(temp)
	    else
	      sdres(5) = 0.0
	    endif
          endif
        endif

        ipt = nsmpr(i) + nxm(i) + nfm(i)
        if(ipt .eq. 0) goto 75

        fixsw = ' '
        if (sw(i) .eq. 0.) fixsw = 'w'

        write(logu, 74) fixsw, nsta(i),
     *  nsmpr(i), asmpw, avres(5), sdres(5),
     *  nxm(i), avres(3), sdres(3), nfm(i), avres(4), sdres(4)

        if(punt .ne. logu) write(punt, 74) fixsw, nsta(i),
     *  nsmpr(i), asmpw, avres(5), sdres(5),
     *  nxm(i), avres(3), sdres(3), nfm(i), avres(4), sdres(4)

   74   format(1x, a1, 1x, a5,
     *  (i4, f4.1, 2f6.3), 2(i4, 2f6.3), 4x, a4)

   75 continue
      return
      end
c end lissum
c locate.for    []
      subroutine locate(iterm)
      save sockets
      include 'params.inc'
      parameter (ndly = 11)
      parameter (pi = 3.14159265)
      parameter (rad = pi/180.)
      parameter (deg = 1.0/rad)
      real bat2,bon2,lat,lon
      character*120 erout
      character*4 malaz, malla, malor, dnstrg
      logical good, eoff, supout, sockets, needsum
      character*1 kwr, iqdo
c     if instset .ne. ' ', then put this new value of inst on summary record
      character*1 instset
      integer punt
      common /punt/ punt
      common /hnu/ instset
      common /anox/ keyd(npa)
      common /bqz/ avrps,avuse
      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 /dhin/ iglob, zup, zdn
      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 /gg/ altla(3), altlo(3), altz(3), altrms(3), nsolut,
     * fla(20), flo(20), fz(20), frms(20), forg(20), maxf, altorg(3)
      common /ghnq/ iexit
      common /gmost/ az(npa)
      common /hf/ cp(72),sp(72)
      common /hl/ nwad, tslope, tsqsl
      real*8 time1, time2
      common /hop/ time1,time2,nopu,notim
      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 /igl/ nglobalzs, globalzs(20)
      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 /pbnoq/ kwr(npa),iqdo(npa)
      common /ph/ nfirst
      common /pno/ lph,keyph(npa),nsum
      character phcard*117, krm*2, seq*5
      common /pno1/ phcard(npa), krm(npa), seq
      character*4 krms
      common /po/ krms(npa)
      common /pqt/ near
      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 /qgnotx/ delta(npa)
      common /qw/ azmv(5), iuse(npa), duse(npa), ndwt, iboxc
      common /qmost/ wt(npa),z
      common /reloc/ irelo, nreloc
      common /ro/ yse,seorg,phi
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      data sockets/.false./
      if (iterm .gt. 0) then
        sockets = .true.
      endif
      if(iprn .ge. 5) write(punt, '(a)') 'begin subroutine locate'
20    if(itest(38) .eq. 11) itest(38) = 1
c read in data for one earthquake
      call phasin
      needsum = .true.
      instset = ' '
      if((lph .gt. 0) .and. (nr .eq. 0)) then
c case of no phase data
        iexit = 1
        write(logfil, 22) kdate, ihr, kmin
        write(punt, 22) kdate, khrmn/100, khrmn - 100*(khrmn/100)
22      format(' could not locate event of ', i6, ' ', i2, ':', i2)
c---- write out warning messages
        rewind(16)
        read(16, '(a)') erout
2296    read(16, '(a)', end=2298) erout
        if(erout(1:3) .eq. 'end') goto 2298
        write(punt, '(a)') erout
        write(logfil, '(a)') erout
        goto 2296
2298    continue

        if((ipun .eq. 2) .or. (ipun .eq. 3)) call npunch('final')
        if(eoff) goto 180
        goto 20
      endif
c at end of file
      if(eoff) goto 180
c no data here, but try main input file again
      if(nr .eq. 0) goto 20
      if ((nr .eq. 1) .and. (ipro .eq. 'stop')) then
c print summary
        call lissum(2)
c terminate, even if using sockets
        stop
      endif
      if ((nr .eq. 1) .and. (ipro .eq. 'rese')) then
c "reset" record found
        if (ichec .eq. 't s ') then
c "reset s" record found
c print summary
          call lissum(2)
c reset summary
          call lissum(1)
c show elapsed time
          call timit(1)
          time1 = 0.0d0
        endif
c get new parameters
        call input1
        goto 20
      endif
c normal earthquake case
      kkdate = kdate
      khr = ihr
      kmin = lbastm
      sec = 0.0
      dum = 0.0
      call tshift(kkdate,khr,kmin,sec,dum)
      kkyr=kkdate/10000
      kkmo=(kkdate-kkyr*10000)/100
      kkday=(kkdate-kkyr*10000-kkmo*100)
      if(iprn .lt. -1) then
c no print in this case
      else if(iprn .eq. -1) then
c one line per event in output and log file
        write(punt, 24) kkyr,kkmo,kkday,khr,kmin
        write(logfil, 24) kkyr,kkmo,kkday,khr,kmin
24      format(1x, i2.2, '/', i2.2, '/', i2.2 , 4x, i2.2, ':', i2.2,
     *  5x, a)
      else
        write(punt,26)
26      format('                                              --------',
     * '--------------------------------------------------------',
     *      /'                                              ---------',
     * '-------------------------------------------------------',
     *      /'                                              ---------',
     * '-------------------------------------------------------')
        write(punt,24)kkyr,kkmo,kkday,khr,kmin,iahead
      endif
c if desired, call wadati
      wslope = 0.
      if(test(49) .ne. 0.) then
        call wadati(nr, nrp, w, ldx, tp, kdate, khrmn, iprn,
     *  ilis, krmp, krms, msta, test, wslope, worig, se3)
        if(wslope .ge. 1.2) then
          tslope = tslope + wslope
          tsqsl = tsqsl + wslope*wslope
          nwad = nwad + 1
          if(freor) savor = worig
          if(test(49) .lt. 0.) then
            savor = worig
            inst = 8
          endif
        else
          if((wslope .ne. 0.) .and. (iprn .ge. 1)) then
            write(punt, '(a,/,a)')
     *      ' VP/VS slope was less than 1.2, so it will not be used to',
     *      ' set initial origin time nor to compute average VP/VS.'
	  endif
        endif
      endif
c find improved first trial location using half space velocity = test(48)
      vhalf = test(48)
      nmax = 10
c if any initial values = 99999. at this time, then compute them from p times!
      if(iprn .ge. 5) write(punt, *) ' savla, savlo = ', savla, savlo
      if( (savla .eq. 99999.) .or. (savez .eq. 99999.) ) then
        if(iprn .ge. 5) then
          write(punt, *) ' number of p times ', nrp
          write(punt, *) ' p times ', (tp(i), i = 1, nrp)
          write(punt, *) ' weights ', (w(i), i = 1, nrp)
          write(punt, *) ' savla, savlo, savez, savor ',
     *          savla, savlo, savez, savor
          write(punt, *) ' call strtep'
        endif
        npinv = nrp
        if(npinv .gt. 9) npinv = 9
        if((iprn .ge. 1) .and. (ilis .gt. 0)) write(punt, 28) npinv
28      format(' compute starting parameters from the first ',
     *  i3, ' p-arrival times.')
        call strtep(c, e, kdx, lat(1), lon(1), nrp, nmax, vez,
     *                   vla, vlo, savor, tp, vhalf, w, test(54))
        if(iprn .ge. 5) then
          write(punt, *) ' strtep input and output:'
          write(punt, *) ' lat(1), lon(1) = ', lat(1)*deg, lon(1)*deg
          write(punt, *) ' nrp, nmax = ', nrp, nmax
          write(punt, *) ' vhalf = ', vhalf
          write(punt, *) ' savor ', savor
          write(punt, *) ' vla, vlo, vez ', vla*deg, vlo*deg, vez
        endif
c
c define savla and savlo
        if(savla .eq. 99999.) then
c in this case, redefine savla and savlo from strtep results
	  if(iprn .ge. 5)
     *     write(punt, *) 'savla = 99999. so use vla and vlo'
          savla = vla
          savlo = vlo
        endif
c
c set starting depth
        if(savez .eq. 99999.) then
c usedly uses the current values of savla and savlo.
c zuse is set to 99999. by usedly if usedly does not set starting depth.
c if zuse is set to 99999. then test(36) is the maximum starting depth.
cd        write(punt, *) 'call usedly to get starting depth'
          call usedly(0, zuse, modset, savla, savlo)
          if(zuse .ne. 99999.) then
            savez = zuse
          else if(test(5) .ne. -99.) then
            savez = test(5) + test(8)
          else
c use vez determined by strtep
c limit vez to the range 15 to test(36) km and round to the nearest 5 km
            if(vez .gt. test(36) + test(8)) vez = test(36) + test(8)
            if(vez .lt. 15. + test(8)) vez = 15. + test(8)
            nz = vez + .5
            savez = nz
          endif
        endif
      endif
      if(savor .eq. 99999.) savor = 0.
c
c find location for one earthquake
      if(itest(50) .gt. 22) itest(50) = 22
      irep = itest(50)
30    continue
cd    write(punt, *) 'call usedly from statement 91'
      call usedly(0, zuse, modset, savla, savlo)
      if((dnstrg(evtype) .eq. 't') .or.
     *   (dnstrg(evtype) .eq. 'n') .or.
     *   (dnstrg(evtype) .eq. 'r')) then
c       teleseism, nuclear or regional
        call velazm
        iexit = 1
        if((ipun .eq. 2) .or. (ipun .eq. 3)) call npunch('final')
        if(eoff) goto 180
        if(ipro .eq. 'more') goto 96
        goto 20
      endif

      if(nglobalzs .ne. 0) then
	call glob_new
	goto 40
      endif

      if((iglob .eq. 0) .and.
     *   ((inst .eq. 0) .or. (inst .eq. 8))) then
        call global
        if ((nsolut .gt. 1) .and. (iprn .ge. 0) .and. (ilis .gt. 0))
     *  write(punt, 32) nsolut, (altz(i), altrms(i), i = 1, nsolut)
32      format(' the number of solutions = ', i2, /,
     *         '     depth      rms', /,
     *         ( 2f10.2))
        call output(1)
        if(iexit .eq. 0) then
          if((kms.eq.0) .and. (needsum)) call mising
c keep running sum of magnitude residuals for sumout
c but do not include multiple solutions of same event in summary
          if(needsum) call lissum(0)
          if(nfirst .ge. iabs(itest(7)) )  call fmplot
c write summary data to temporary files 14 and 15, call geterr and
c then add the the data to files 4 and 11
          if(ipun .ne. 0) call npunch('temp ')
          call geterr(zup, zdn, rmslim)
c add zup and zdn to summary record and write output to .sum and .arc
cd        write(punt, *) 'call adderr with zup, zdn = ', zup, zdn
          if(ipun .ne. 0) call adderr(zup, zdn)
          if (iprn .ge. 0) write(punt, 34) zup, zdn, yse, rmslim
34        format(' depth may decrease by ', f10.2,
     *    ' or increase by ', f10.2, ' given a reading ', /,
     *    ' standard error of ', f8.3, ' and rms limit of', f8.3)
        else
c iexit = 1, so there is insufficient data for a solution
          if((ipun .eq. 2) .or. (ipun .eq. 3)) then
            call npunch('final')
          endif
        endif
        if((itest(44) .eq. 1) .and. (inst .ne. 9)) then
c since (44) = 1, check if this is a debug event to rerun
          if((nedit .eq. 3) .and. (.not. good)) then
            itest(44) = 100
            inst = 0
            goto 40
          endif
        endif
c check for critical solution based on itest(44) = 2
        if ((itest(44) .eq. 2) .and. (inst .ne. 9)) then
c  continue iteration with only critical stations
          if (iprn .ge. 5) write(punt, '(a)')
     *    ' continue iteration with only critical stations'
          itest(44) = 200
          inst = 0
          goto 40
        endif
        goto 50
      endif
c
c normal solution (not global)
40    call quakes
c if iexit is set = 2 in quakes there are too few phases to
c rerun event with only critical stations.
      if(iexit .eq. 2) then
        if(itest(44) .eq. 100) itest(44) = 1
        if(itest(44) .eq. 200) itest(44) = 2
        goto 50
      endif
      call output(1)
c if iexit is set = 1 in quakes there is insufficient data
      if(iexit .eq. 1) then
        if((ipun .eq. 2) .or. (ipun .eq. 3)) then
          call npunch('final')
        endif
        goto 20
      endif
      if (abs(test(6)) .gt. 0.01) then
        call boxau
      endif
      if(itest(44) .eq. 100) then
        itest(44) = 1
        goto 50
      endif
      if (itest(44) .eq. 200) then
        itest(44) = 2
        goto 50
      endif
      if((kms.eq.0) .and. (needsum) .and.
     * (.not. supout)) call mising
      if( (nfirst .ge. iabs(itest(7))) .and.
     * (.not. supout)) call fmplot
c keep running sum of magnitude residuals for sumout
c do not include multiple solutions of same event in summary
      if((needsum) .and. (.not. supout)) call lissum(0)
      if(ipun .ne. 0) call npunch('final')
      if(itest(50) .ne. 0) then
        if(irep .eq. 0) then
c place a blank record between sets of summary records
          write(4, '(1h1)')
c go on to the next set of phase data
          supout = .false.
          goto 50
        endif
        if(irep .eq. itest(50)) then
c this is the first extra run of this event
c set sepout to true to supress calling output and
c computing magnitude
          supout = .true.
c set inst = 1 to fix z
          inst = 1
c set starting depth to surface of model
          savez = 0.0
          zdif = 1.
        else
c on later extra runs
          savez = savez + zdif
          zdif = 1. + savez/5.
        endif
        irep = irep - 1
        goto 30
      endif
      if((itest(44) .eq. 1) .and. (inst .ne. 9)) then
c since (44) = 1, check if this is a debug event to rerun
        if((nedit .eq. 3) .and. (.not. good)) then
          itest(44) = 100
          goto 30
        endif
      endif
      if((itest(44) .eq. 2) .and. (inst .ne. 9)) then
c since (44) = 2, rerun
        itest(44) = 200
        goto 30
      endif
50    continue
      time1=time2
      if(itest(38) .eq. 1) then
        itest(38) = 11
        knst = 0
        goto 30
      else if(itest(38) .eq. 11) then
        itest(38) = 1
      endif
      if(ipro .ne. 'more') goto 20
c
c read next instruction record
96    read(inpt, '(a)') icard
      phcard(lph) = icard
      if(isa.eq.1) write(1,'(a)') icard
      read(icard, 100)
     *     ipro, ichec,  knst, inst, zres,     ala1,
     *    isnla,  ala2,    alo1, isnlo,  alo2,      org1, org2
100   format(a4,    a4, 9x,i1,   i1, f5.2, 16x,f2.0,
     *       a1,  f5.2, 5x,f3.0,    a1,  f5.2, 11x, f2.0, f5.2)
      instset = ' '
      malaz = icard(21:24)
      malla = icard(43:46)
      malor = icard(74:77)
      if(icard(9:9) .ne. ' ') evstat = icard(9:9)
      if(icard(10:10) .ne. ' ') evtype = icard(10:10)
      write(punt,110) ipro, ichec, evstat, evtype,   knst,   inst,
     *     zres,     ala1, isnla, ala2,     alo1, isnlo, alo2,
     *     org1, org2

110   format(' ipro check evstat evtype knst inst depth   ',
     *' latitude   longitude, origin time', /,
     * 1x,             a4,     a4, 5x,a1, 6x,  a1, 5x, i1, 4x, i1,
     * 2x, f6.2, 2x, f4.0,    a1, f5.2, 1x, f5.0,    a1, f5.2,
     * 3x, f3.0, f5.2, /, 25x, '--- run again ---  ')

      ipro = dnstrg(ipro)
      ichec = dnstrg(ipro)
      if(itest(38) .eq. 2) knst = 1
      if(itest(38) .eq. 3) knst = 0
c calculate first trial hypocenter
      if(malla .ne. '    ') then
        la = ala1 + .0001
        lo = alo1 + .0001
        call fold2(savla,savlo,la,isnla,ala2,lo,isnlo,alo2)
      endif
      if(malaz .ne. '    ') savez = zres + test(8)
c define starting origin time
      savor = 0.
      freor = .true.
      if(malor .ne. '    ') then
        savor = 60.*(org1 - lbastm) + org2
        freor = .false.
      else if(inst .eq. 8) then
        inst = 0
      endif
c convert inst=7 to inst=9 with origin free - this is
c the option to use for quary blasts with unknown origin time.
      if(inst .eq. 7) then
        inst = 9
        freor = .true.
      endif
      if((dnstrg(evtype) .eq. 't') .or.
     *   (dnstrg(evtype) .eq. 'n') .or.
     *   (dnstrg(evtype) .eq. 'r')) then
c       teleseism, nuclear or regional
        call velazm
        iexit = 1
        if((ipun .eq. 2) .or. (ipun .eq. 3)) call npunch('final')
        if(ipro .eq. 'more') goto 96
        goto 20
      endif
      needsum = .false.
      if((ipro .eq. '    ') .or. (ipro .eq. 'more')) goto 30
160   write(punt,170) ipro, ichec
170   format(/' xxxerrorxxx ', 2a4, ' skipped.  must be blank or more,')
      goto 20
c
c after all earthquakes are processed (eoff is true)
180   if(sockets) then
        return
      endif
      call lissum(2)
c     show elapsed time
      call timit(1)
      write(punt, *) 'irelo = ', irelo, ' nreloc = ', nreloc
      if(irelo .gt. 0) then
c relocate events with new station delays irelo times
        nreloc = nreloc + 1
        if(nreloc .gt. irelo) stop
c rewind input(inpt), archive(11), summary(4), output(9)
        write(punt, '(a)') ' rewind files and rerun events'
        if(inpt .eq. 5) then
          write(punt, *)
     *     ' input is standard input, which can not be rewound'
          write(punt, *)
     *     ' therefore, the relocate option may not be used'
          stop 'abort from locate'
        endif
        rewind (inpt)
        rewind (11)
        rewind (4)
        rewind (9)
c       reset summary
        call lissum(1)
        eoff = .false.
        time1 = 0.0d0
        call input1
        goto 20
      endif
      return
      end
c end locate
c median.f    []
c	parameter (n = 8)
c	dimension x(n)

c	x(1) = 1.90560
c	x(2) = 1.27506
c	x(3) = 2.35488
c	x(4) = 2.22930
c	x(5) = 1.90560
c	x(6) = 1.75282
c	x(7) = 2.01101
c	x(8) = 1.84767

c	write(6,'(''Data values:'')')
c	write(6,'(4e20.10)') x
c	call median(x,n,xmed,ierr)
c	print *, 'xmed = ', xmed
c	stop
c	end

	subroutine median(x,n,xmed,ierr)

c  use iterative method to determine the median:

c      N                                          __ N     xj
c     ___     xj - xmed                           >     -----------
c     \      ----------- = 0    ==>    xmed =     -- 1  |xj - xmed|
c     /      |xj - xmed|                         -----------------
c     ---                                         __ N      1
c     j=1                                         >     -----------
c                                                 -- 1  |xj - xmed|

c  modified by C. Stephens from routine mdian2 in Numerical Recipes, by
c  Press, Flannery, Teukolsky, and Vetterling, Cambridge Univ Press, 1986,
c  pp. 460-462.

c  the variable nequal was added to ensure convergence in the case where 1 or
c  more data values are equal to the current guess for the median (= a)

c  when more than one data value is equal to the true median and the current
c  guess (= a) is near but not equal to this value, then convergence can be
c  slowed by successive iterations that oscillate about the true median.
c  under these conditions, either xp is equal to xm from the previous iteration,
c  xmlast are used to test for this case.

c  error return code ierr = 0 for no errors
c			  = maxit (= max(100,100*log(n)) if convergence does not
c				  occur within this number of iterations, in
c				  which case xmed is set to the current guess

	integer*4 n
	real*4 x(n)

	integer*4 np, nm, nequal, j, halfn
	real*4 xx, xp, xm, sumx, sum, eps, dum
	real*4 ap, am, aa, a

	parameter (big=1.0e30, afac=1.5, amp=1.5)

	ierr = 0

	a = 0.5*(x(1)+x(n))
	eps = abs(x(n)-x(1))
	am = -big
	ap = big
	halfn = n/2
	xplast = big
	xmlast = -big
	maxit = 100*nint(log(1.*n))
	if (maxit .lt. 100) maxit = 100
	nit = 0

1	nit = nit + 1
	if (nit .gt. maxit) then
	  ierr = maxit
	  xmed = a
	  return
	endif
	sumx = 0.0
	sum = 0.0
	np = 0
	nm = 0
	nequal = 0
	xm = -big
	xp = big

	j = 0
	do while (j .lt. n)
	  j = j + 1
	  xx = x(j)

	  if (xx .eq. a) then
	    nequal = nequal + 1

	  else

	    if (xx .gt. a) then
	      np = np + 1
	      if (xx .lt. xp) xp = xx

	    else if (xx .lt. a) then
	      nm = nm + 1
	      if (xx .gt. xm) xm = xx
	    endif

	    dum = 1.0/(eps+abs(xx-a))
	    sum = sum + dum
	    sumx = sumx + xx*dum

	  endif

	enddo

c -- check for oscillations about true median

	if (xm .eq. xplast) then
	  a = xm
	  go to 1
	else if (xp .eq. xmlast) then
	  a = xp
	  go to 1
	else
	  xmlast = xm
	  xplast = xp
	endif

c -- adjust xp,np and xm,nm if some data values equal current guess

	if (nequal .gt. 0) then

	  if (np .lt. halfn) then
	    np = np + nequal
	    if (np .gt. halfn) np = halfn
	    xp = a
	  endif

	  if (nm .lt. halfn) then
	    nm = nm + nequal
	    if (nm .gt. halfn) nm = halfn
	    xm = a
	  endif

	endif

c -- check distribution of points about current guess

	if (np-nm .gt. 1) then
c					guess too low
	  am = a
	  aa = xp + max(0., sumx/sum-a)*amp
	  if (aa .gt. ap) aa = 0.5*(a+ap)
	  eps = afac*abs(aa-a)
	  a = aa
	  go to 1

	else if (nm-np .gt. 1) then
c					guess too high
	  ap = a
	  aa = xm + min(0., sumx/sum-a)*amp
	  if (aa .lt. am) aa = 0.5*(a+am)
	  eps = afac*abs(aa-a)
	  a = aa
	  go to 1

c -- success

	else

	  if (mod(n,2) .eq. 0) then
c					n is even
	    if (np .eq. nm) then
	      xmed = 0.5*(xp+xm)
	    else if (np .gt. nm) then
	      xmed = 0.5*(a+xp)
	    else
	      xmed = 0.5*(xm+a)
	    endif

	  else
c					n is odd
	    if (np .eq. nm) then
	      xmed = a
	    else if (np .gt. nm) then
	      xmed = xp
	    else
	      xmed = xm
	    endif

	  endif

	endif

	return

	end
c mising.for    []
      subroutine mising
c check to see if data from stations that were not used in
c------- this solution would significantly improve hypocenter
      include 'params.inc'
      parameter (ndly = 11)
      real lat,lon,latep,lonep,mag
      character*4 iahead*60, msta*5, nsta*5, icard*110
      integer punt
      common /punt/ punt
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /amo/ tempg(npa), ntgap
      common /dmost/ ipun,ivlr,blank
      common /iclmpq/ lat(nsn),lon(nsn)
      common /ilmpu/ ns
      common /ilpu/ sw(nsn),ndate(nsn),nhr(nsn),mno(nsn),ielv(nsn)
      common /ilpx/ calr(5, nsn),fmgc(nsn),xmgc(nsn)
      common /ilt/ vthk(2,nsn),ipdly(nsn),mod(nsn),ipthk(nsn)
      common /ilv/ c(nsn),e(nsn)
      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 /iox/ prr(nsn),iuses
      integer fmwt, xmwt
      common /ioxl/ fmwt(nsn), xmwt(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 /ip1/ rsew(4)
      common /omnfh/ dmin,dmin3,sminp
      common /pfnoqv/ kdate,khrmn,khr
      common /pm/ jdx(nsn)
      common /pox/ fmp(npa),prx(npa),icent1,icent2
      common /qmost1/ lonep,ni,latep
      common /reloc/ irelo, nreloc
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      common /xfmno/ mag
      ihd=0
      nj = ntgap + 1
      tempg(nj)=tempg(1)+360.
      if (mag .eq. blank) then
        tdel=100.
      else
        tdel = 10.**(1.48 + mag/2.)
      endif
      do 130 i=1,ns
    5 if (ndate(i) .eq. 99999998) goto 130
      if (jdx(i) .eq. 1) goto 130
      call delaz(latep,lonep,deli,deldeg,azi,lat(i),lon(i))
      if (deli .gt. tdel) goto 130
      if (azi .le. tempg(1)) azi=azi+360.
      do 10 j=2,nj
      if (azi .lt. tempg(j)) goto 20
   10 continue
      j=nj
   20 exgap=tempg(j)-tempg(j-1)
      rdgap=tempg(j)-azi
      tgap=azi-tempg(j-1)
      if (tgap .lt. rdgap) rdgap=tgap
      if (deli.gt.dmin3 .and. rdgap.lt.30.) goto 130
      if (azi .ge. 360.) azi=azi-360.
      if (ihd .eq. 1) goto 22
      write(punt,21)
   21 format(/10x,45hmissing station  delta   azim  ex-gap  rd-gap)
      ihd=1
c     make sure station record has not expired.
   22 if (ndate(i) .gt. icent2*1000000+kdate) goto 100
      if ((ndate(i) .eq. icent2*1000000+kdate) .and.
     *  (nhr(i) .ge. khrmn/100)) goto 100
      call update(indexs, nsta, ielv,
     *         lat, lon, c, e, sw,
     *         klas, calr, xmgc, fmwt, xmwt, fmgc, ipcod, iscod,
     *         tpdly, revp, ndate, nhr, ns,
     *         exdly, icent2, kdate, khrmn,
     *         infil, iofil)
      goto 5
  100 write(punt,125) nsta(i),deli,azi,exgap,rdgap
  125 format(21x,a5,2f7.1,2f8.1)
  130 continue
      return
      end
c end mising
c npunch.for    [pc]
      subroutine npunch(style)
c----
c print summary record with regress information
c write summary record with regress information
c write summary record without regress information
c write station data
      include 'params.inc'
      parameter (ndly = 11)
      integer sumout, arcout
      character*5 style, upstrg*1
c style can be 'final' or 'temp '.  if 'temp ', then write output
c on 4 to temporaty file 14 and output on 11 to temporary file 15
      logical sumpr
      real bat2,bon2,mag
      character*1 iqpr, iwrk, fmit*6, dnstrg
      character aline*133, pline*117
      character*4 iahead*60, msta*5, nsta*5, icard*110
      common /an/ n14, n15
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /dbiln/ ioldq
      common /dhin/ iglob, zup, zdn
      logical medmag
      common /dinx/ imag,imagsv,medmag
      common /dmost/ ipun,ivlr,blank
      common /gmost/ az(npa)
      common /omnfh/ dmin,dmin3,sminp
      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 /pot/ ts(npa),ihr,model(npa), keyphi(npa)
      common /pox/ fmp(npa),prx(npa),icent1,icent2
      integer punt
      common /punt/ punt
      common /ghnq/ iexit
      common /hpn/ wslope
      common /idno/ ksel,ksort
      common /ilotv/ elvdly(npa)
      common /imost/ test(100)
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /in/ irmo,iryr,odmag
      character*1  magke
      common /in1/ magke
      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)
      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)
c if instset .ne. ' ', then put this new value of inst on summary record
      character*1 instset
      common /hnu/ instset
      common /ofnv/ kmin
      character*1 nq
      common /ofnv1/ nq
      common /ofln/ sec
      common /olnx/ xmag(npa),fmag(npa),avxm,avfm,xmmed,fmmed
      common /on/ ix,iy,iz,ax1,axz
      character iqax*1
      common /on1/ iqax
      common /onx/ kluse(npa), cluse(npa), peruse(npa), ampuse(npa)
      common /pfnoqv/ kdate,khrmn,khr
      common /pgnoqv/ p(npa),jmin(npa)
      common /pgnov/ dt(npa)
      common /pgoq/ s(npa)
      character*4 ipro, ichec, evtype*1, evstat*1, root*256
      common /dbhip/ ipro, ichec, evtype, evstat, root
      common /phoqn/ inst,knst
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /pn/ itpused(npa)
      common /pnl/ ww(npa)
      common /pno/ lph,keyph(npa),nsum
      character phcard*117, krm*2, seq*5
      common /pno1/ phcard(npa), krm(npa), seq
      character*1 kwr, iqdo
      common /pbnoq/ kwr(npa),iqdo(npa)
      common /pnoqtx/ ldx(npa)
      common /qmost/ wt(npa),z
      common /qgnotx/ delta(npa)
      common /qgo/ org
      common /rbno/ idip(3),iaaz(3),a(3,3)
      common /rfgnoq/ se(4)
      common /ro/ yse,seorg,phi
      common /tmost/ x(4,npa)
      common /tonxb/ t(npa),fms(npa)
      common /xfmno/ mag
      character*1 mgndx
      common /xfmno1/ mgndx
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      character*(4) tmpsum_type, tmparc_type, fname*256
      character acent2*2, fooline*133
      integer tmpsum_unit, tmparc_unit
      parameter (tmpsum_type  = '.3sc')
      parameter (tmparc_type  = '.4sc')
      parameter (tmpsum_unit  = 14)
      parameter (tmparc_unit  = 15)

      write(acent2, '(i2.2)') icent2
c* (unix
c      close(tmpsum_unit)
c      length_root = lentru(root)
c      fname = root(1:length_root)//tmpsum_type
c      call openfl(tmpsum_unit, fname, 'scratch', 'null', 'none',
c     *      'noprint', 300)
c
c      close(tmparc_unit)
c      fname = root(1:length_root)//tmparc_type
c      call openfl(tmparc_unit, fname, 'scratch', 'null', 'none',
c     *      'noprint', 0)
c* unix)
c* (pc
      rewind tmpsum_unit
      rewind tmparc_unit
c* pc)
c* (vax
c      rewind tmpsum_unit
c      rewind tmparc_unit
c* vax)

      sumpr = .false.
      rmag = mag
c n14 and n15 count the number of records saved temporarily, until
c geterr has been run to compute zup and zdn
      n14 = 0
      n15 = 0
c set up unit numbers for output
      if(style .eq. 'final') then
        sumout = 4
        arcout = 11
      else
        sumout = tmpsum_unit
        arcout = tmparc_unit
      endif
c	write(punt, '(a, 2i5)') 'sumout, arcout = ', sumout, arcout
      if((magke .eq. ' ') .or. (magke .eq. 'x') .or.
     *   (magke .eq. 'f') .or. (magke .eq. 'a') .or.
     *   (magke .eq. 'k')) goto 3
c     in this case, save old magnitude value
      rmag = odmag
      mgndx = magke
    3 nfile = 4
      if(ipun .lt. 0 .or. ipun .gt. 4) return
c fix up the date and time
      khr = ihr
      kmin = lbastm
      sec = org
      kkdate = kdate
      call tshift(kkdate,khr,kmin,sec,dum)

c take care of the fake solution cases
      if(iexit .eq. 1) then
c       don't write out fake arrival data if there wasn't a solution
        if (ipun .eq. 4) return
	if ((inst .eq. 9) .and. (keyph(1) .eq. -2)) then
c         for inst=9, use original summary record if there was one
          pline = acent2//phcard(1)(1:115)
	  goto 137
        else if((dnstrg(evtype) .eq. 't') .or.
     *    (dnstrg(evtype) .eq. 'n') .or.
     *    (dnstrg(evtype) .eq. 'r')) then
c         for tele, nuclear, or regional, use old summary record
          pline = acent2//phcard(1)(1:115)
	  goto 137
	else
          goto 225
	endif
      endif
c
c key to keyph values
c      -1  comment record
c      -2  summary record
c      -3  instruction record
c      -4  deleted station record
c     -11  decode error
c      .ge. 1, index number of phase used in solution
c
c print and write summary information with regress information
   10 call formf(sec,iisec,4,2)
      call formf(bat2,iilat,4,2)
      call formf(bon2,iilon,4,2)
      call formf(z-test(8),iiz,5,2)
      call formf(dmin,iidmin,3,0)
      call formf(rms,iirms,4,2)
      iqpr = iqax
      if(ioldq .eq. 1) iqpr = nq
      if(iprn .lt. 3) goto 130
      write(punt, '(3a)')
     *  '     date   origin      lat      long    depth    mag no',
     *  ' gap  d3   rms az/dp    se az/dp    se az/dp    se  iqpr',
     *  ' isetno nswt  seq  s-p'
c
c*********************create summary record for print file**************
      aline = ' '

c     write(punt, '(a,2i10, f10.2, 2i10 )')
c    *  'khr, kmin, sec, kdate, kkdate',
c    *  khr, kmin, sec, kdate, kkdate
      write(aline, 126)  kkdate, khr, kmin, sec, lat1, isnla, bat2,
     * lon1, isnlo, bon2, z-test(8)
  126 format(1x, i6.6, 1x ,i2, i2, f6.2, i3, a1, f5.2, i4, a1,
     * f5.2, 1x, f6.2)
c
      if (rmag .ne. blank) then
        call formit(rmag, frmag, fmit, 3, 1)
        write(aline(49:51), fmit) frmag
      endif
c
      write(aline(52:74), 127) nrwt, igap, dmin, rms, iaaz(ix), idip(ix)
  127 format(i3, i4, f5.0, f5.2, i3, '/' ,i2)
c
      if(se(ix) .ne. blank) then
        call formit(se(ix), frsei, fmit, 5, 1)
        write(aline(76:80), fmit) frsei
      endif
c
      write(aline(81:86), 128) iaaz(iy), idip(iy)
  128 format(i3, '/' ,i2)
c
      if (se(iy) .ne. blank) then
        call formit(se(iy), frsej, fmit, 5, 1)
        write(aline(88:92), fmit) frsej
      endif
c
      write(aline(93:98), 128) iaaz(iz), idip(iz)
c
      if (se(iz) .ne. blank) then
        call formit(se(iz), frsek, fmit, 5, 1)
        write(aline(100:104), fmit) frsek
      endif
c
c     write(aline(105:132), 129) iqpr, mgndx, nswt, seq, sminp
c 129 format(5x ,a1, 6x ,a1, 3x ,i2, a5, f5.2)
      write(aline(105:127), 129) iqpr, mgndx, nswt, seq
  129 format(5x ,a1, 6x ,a1, 3x ,i2, a5)
      if(sminp .ne. blank) then
        write(aline(128:132), '(f5.2)') sminp
      endif
c
c     add century
      fooline = aline(2:132)
      aline = acent2//fooline
      write(punt, '(1x, a)') aline
c********************finished with printed summary record****************
c
  130 continue
c*********************create summary record******************************
      pline = ' '
      write(pline, 131) kkdate, khr, kmin, iisec, lat1, isnla, iilat,
     * lon1, isnlo, iilon
  131 format(            i6.6,  i2,   i2,    i4,   i2,    a1,    i4,
     *   i3,    a1,    i4)
      write(pline(111:115), '(i5)') iiz
      if(test(9) .eq. 0.0) then
        write(pline(30:34), '(i5)') iiz
      else
        if(iiz .lt. 0) then
          write(pline(30:34), '(a)') '  -00'
        else
          write(pline(30:34), '(i5)') iiz
        endif
      endif
c
      if (rmag .ne. blank) then
        call riorbk(rmag, jmag, fmit, 2, 1)
        write(pline(35:36), fmit) jmag
      endif
c
      write(pline(37:54), 133) nrwt, igap, iidmin, iirms, iaaz(ix),
     * idip(ix)
  133 format(                    i3,   i3,     i3,    i4,       i3,
     *       i2)
c
      if (se(ix) .ne. blank) then
        call riorbk(se(ix), iisei, fmit, 4, 2)
        write(pline(55:58), fmit) iisei
      endif
c
      write(pline(59:63), '(i3, i2)') iaaz(iy), idip(iy)
c
      if (se(iy) .ne. blank) then
        call riorbk(se(iy), iisej, fmit, 4, 2)
        write(pline(64:67), fmit) iisej
      endif
c
      if(medmag) then
c use median xmag
        if (xmmed .ne. blank) then
          call riorbk(xmmed, iavxm, fmit, 2, 1)
          write(pline(68:69), fmit) iavxm
        endif
      else
c use average xmag
        if (avxm .ne. blank) then
          call riorbk(avxm, iavxm, fmit, 2, 1)
          write(pline(68:69), fmit) iavxm
        endif
      endif
c
      if(medmag) then
c use median fmag
        if (fmmed .ne. blank) then
          call riorbk(fmmed, iavfm, fmit, 2, 1)
          write(pline(70:71), fmit) iavfm
        endif
      else
c use average fmag
        if (avfm .ne. blank) then
          call riorbk(avfm, iavfm, fmit, 2, 1)
          write(pline(70:71), fmit) iavfm
        endif
      endif
c
      write(pline(72:72), '(a1)') evstat
c
      if (se(iz) .ne. blank) then
        call riorbk(se(iz), iisek, fmit, 4, 2)
        write(pline(73:76), fmit) iisek
      endif
c
      write(pline(77:96), 134) iqpr, mgndx, nswt,    ipro,
     * irmo, iryr, evtype, phcard(lph)(19:19), seq
  134 format(                    a1,    a1,   i2, '/', a4,
     *   i2,   i2,     a1,                 a1,  a5)
      if(instset .ne. ' ') pline(91:91) = instset
c
      if(sminp .ne. blank) then
        call riorbk(sminp, isminp, fmit, 4, 2)
        write(pline(97:100), fmit) isminp
      endif
c
c     if(iglob .eq. 0) then
      if((iglob .eq. 0) .and. ((inst .eq. 0) .or.
     *  (inst .eq. 8))) then
        call formal(zup, izup, 2, 0, fmit, azup)
        if (fmit .eq. ' ') then
          write(pline(101:102), '(i2)') izup
        else
          write(pline(101:102), fmit) azup
        endif
        call formal(zdn, izdn, 2, 0, fmit, azdn)
        if (fmit .eq. ' ') then
          write(pline(103:104), '(i2)') izdn
        else
          write(pline(103:104), fmit) azdn
        endif
      endif
c
      if (wslope .ne. 0.) then
        call riorbk(wslope, iiwslp, fmit, 4, 2)
        write(pline(105:108), fmit) iiwslp
      endif
c
c compute the number of phases weighted out due to large residuals
      nwtout = 0
      do 136 i = 1, lph
        if(keyph(i) .gt. 0) then
          k = keyph(i)
c	  print *, 'nwtout, k, kwr(k) ', nwtout, k, kwr(k)
c p phases
          if((kwr(k) .eq. 'b') .or. (kwr(k) .eq. 'm') .or.
     *      (kwr(k) .eq. 'j')) nwtout = nwtout + 1
          if(ldx(k) .ne. 0) then
c s phases
            kk = ldx(k)
            if((kwr(kk) .eq. 'b') .or. (kwr(kk) .eq. 'm') .or.
     *        (kwr(kk) .eq. 'j')) nwtout = nwtout + 1
          endif
        endif
136   continue
      if(nwtout .gt. 99) nwtout = 99
      write(pline(109:110), '(i2)') nwtout
c
c     for now (12/21/91) keep upper case characters
      pline(17:17) = upstrg(pline(17:17))
      pline(25:25) = upstrg(pline(25:25))
      pline(77:77) = upstrg(pline(77:77))
      pline(78:78) = upstrg(pline(78:78))
c
c     add century
      fooline = pline(1:115)
      pline = acent2//fooline
c
c*************write out current summary record parameters********
c	write(punt, '(a, i5)') 'ipun = ', ipun
137   if (ipun .lt. 3) then	
        write(sumout, '(a)') pline
c keep the summary file up to date
c* (unix
c	if(sumout .eq. 4) call flush(sumout)
c* unix)
c* (pc
c* pc)
c* (vax
c* vax)
c	write(punt, '(a, i5)') 'just flushed sumout unit = ', sumout
        if(style .ne. 'final') n14 = n14 + 1
      endif
      if (ipun .eq. 1) return
      goto 260
c
c****************come here if no solution was found**************
  225 continue
c before writing out a fake summary record, check if the current
c summary record (if there is one) is also a fake, based on having
c a blank rms field (cols 46:49)
      if(((phcard(1)(46:49) .ne. '    ') .and.
     *   (keyph(1) .eq. -2)) .or.
     *   (keyph(1) .ne. -2)) then
c there is a valid summary record from a previous run to preserve, or
c there is no summary record from previous run, so
c generate a fake summary record
        ibegin = 1
        pline = ' '
	if(keyph(1) .eq. -2) then
	  pline(1:10) = phcard(1)(1:10)
	else
          write(pline, '(i6.6, i4)') kkdate, khrmn
	endif
c
        pline(17:17) = 'n'
        pline(25:25) = 'w'
        if(rmag .ne. blank) then
          call riorbk(rmag, jmag, fmit, 2, 1)
          write(pline(35:36), fmit) jmag
        endif
c
        if(avxm .ne. blank) then
          call riorbk(avxm, iavxm, fmit, 2, 1)
          write(pline(68:69), fmit) iavxm
        endif
c
        if (avfm .ne. blank) then
          call riorbk(avfm, iavfm, fmit, 2, 1)
          write(pline(70:71), fmit) iavfm
        endif
c
        write(pline(72:72), '(a1)') evstat
c
        pline(78:81) = 'k  /'
c
        write(pline(82:96), 227) ipro, irmo, iryr, evtype,
     *  phcard(lph)(19:19), seq
  227   format(a4, 2i2, 2a1, a5)
        if(instset .ne. ' ') pline(91:91) = instset

c       for now (12/21/91) keep upper case characters
        pline(17:17) = upstrg(pline(17:17))
        pline(25:25) = upstrg(pline(25:25))
        pline(77:77) = upstrg(pline(77:77))
        pline(78:78) = upstrg(pline(78:78))
	fooline = acent2//pline(1:115)
	pline = fooline
      else
c there was already a fake summary record, so
c do not generate another fake summary record
        ibegin = 2
        if(phcard(1)(81:81) .eq. '/') then
	  pline = acent2//phcard(1)
	else
	  pline = phcard(1)
	endif
      endif
c
      if(ipun .lt. 3) then
        write(sumout, '(a)') pline
c keep the summary file up to date
c* (unix
c	if(sumout .eq. 4) call flush(sumout)
c* unix)
c* (pc
c* pc)
c* (vax
c* vax)
        if(style .ne. 'final') n14 = n14 + 1
      endif
c
      if(ipun .eq. 1) return
c for archive option, write out first old summary record if it was
c real, followed by all privious summary and phase records
      write(arcout, '(a)') pline(1:lentru(pline))
      if(style .ne. 'final') n15 = n15 + 1
      do 228 i = ibegin, lph
	if(keyph(i) .eq. -2) then
c* (unix
c          phcard(i)(81:81) = '\\'
c* unix)
c* (vax
c* (pc
         phcard(i)(81:81) = '\'
c* pc)
c* vax)
          write(arcout, 144) acent2//phcard(i)(1:lentru(phcard(i)))
	else
          if(upstrg(phcard(i)(1:2)) .ne. 'C*')
     *      phcard(i)(65:65) = upstrg(phcard(i)(65:65))
          write(arcout, 144) phcard(i)(1:lentru(phcard(i)))
          if(style .ne. 'final') n15 = n15 + 1
	endif
  228 continue
      return
c
c **************output archive-phase data for each station*************
  260 iqpr = iqax
      if(ioldq .eq. 1) iqpr = nq
c archive current summary record parameters
      write(arcout, '(a)') pline(1:lentru(pline))
      if(style .ne. 'final') n15 = n15 + 1
      nott = 0
c
c********************** loop through stations
      do 310 i = 1, lph
      pline = ' '
c skip on down if ipun = 4 and generate perfect data
      if(ipun .eq. 4) then
        if(keyph(i) .lt. 1) goto 310
        goto 302
      endif
c old summary records, skip the first one, be sure col 81 has \
c     for the secondary summary records
      if(keyph(i) .eq. -2) then
        nott = nott + 1
        if(nott .eq. 1) goto 310
c* (unix
c        phcard(i)(81:81) = '\\'
c* unix)
c* (vax
c* (pc
        phcard(i)(81:81) = '\'
c* pc)
c* vax)
        pline = acent2//phcard(i)(1:115)
        write(arcout, 144) pline(1:lentru(pline))
  144   format(a)
        if(style .ne. 'final') n15 = n15 + 1
        goto 310
      endif
      if(keyph(i) .le. 0) then
c   write out comment record, deleted station record or instruction record
        write(arcout, 144) phcard(i)(1:lentru(phcard(i)))
        if(style .ne. 'final') n15 = n15 + 1
        goto 310
      endif
      k = keyph(i)
      kji = kdx(k)
      call formf(x(4,k),ipr,5,2)
      iain = ain(k)  + 0.5001
      iwrk = '    '
      isr = 0
      isse = 0
      if(ldx(k) .ne. 0) then
c s data
        kk = ldx(k)
        iwrk = kwr(kk)
        call formf(x(4,kk),isr,5,2)
        if(wt(kk) .ne. 0.0) then
          setmp = yse/sqrt(wt(kk))
          call formal(setmp, isse, 3, 2, fmit, xsse)
          if(fmit .eq. '      ') then
            write(pline(81:83), '(i3)') isse
          else
            write(pline(81:83), fmit) xsse
          endif
        else
          write(pline(81:83), '(a)') '99.'
        endif
        goto 300
      endif
      if(ksmp(k) .ne. 1) then
c s-p data
        call formf(x(4,k),isr,5,2)
      endif
  300 continue
      call formf(dly(kno,kji),ipdly,3,1)
      call formf(sdly(kno,kji),isdly,3,1)
      call formf(elvdly(k),ieldly,3,1)
      call formf(delta(k),idelt,4,1)
      call formf(az(k),iaz,3,0)
c
      write(pline(1:50), 301) phcard(i)(1:24), idelt, iaz,
     * phcard(i)(32:40), iain, phcard(i)(44:50)
  301 format(a24, i4, i3, a9, i3, a7)
c
      call formal(t(k), iptt, 4, 2, fmit, xiptt)
      if(fmit .eq. '      ') then
        write(pline(51:54), '(i4)') iptt
      else
        write(pline(51:54), fmit) xiptt
      endif
      if(wt(k) .ne. 0.0) then
        setmp = yse/sqrt(wt(k))
        call formal(setmp, ise, 3, 2, fmit, xse)
        if(fmit .eq. '      ') then
          write(pline(55:57), '(i3)') ise
        else
          write(pline(55:57), fmit) xse
        endif
      else
        write(pline(55:57), '(a)') '99.'
      endif
c
      write(pline(58:58), '(a1)') kwr(k)
c save the instrument period and gain characters
      write(pline(59:60), '(a)') phcard(i)(59:60)
c
c replace calr with siemens gain state and a1vco gain range state
      write(pline(61:62), '(a)') phcard(i)(61:62)
c
      write(pline(63:80), 303) phcard(i)(63:75), ipr
  303 format(a13, i5)
c
      write(pline(84:98), 3031) iwrk, isr, ipdly,
     * isdly, ieldly
3031  format(a1, i5, 3i3)
c
      if (xmag(k) .ne. blank) then
        write(pline(99:100), '(i2)') kluse(k)
        call riorbk(xmag(k), ixmag, fmit, 2, 1)
        write(pline(101:102), fmit) ixmag
      endif
c
      if (fmag(k) .ne. blank) then
        call riorbk(fmag(k), ifmag, fmit, 2, 1)
        write(pline(103:104), fmit) ifmag
      endif
c
c satellite hops
c value of nhop is no longer preserved, because the station history can
c accomodate two different telemetry delays
c      if(phcard(i)(110:110) .ne. ' ') then
c preserve any non blank value
c        pline(105:110) = phcard(i)(105:110)
c      else
c note that itpused is based soley on the p-phase
      if (itpused(k) .ne. 0) then
        nhop = -tpdly(itpused(k), kji)/.27 + .1
      else
        nhop = 0
      endif
      write(pline(105:110), '(a, i1)') phcard(i)(105:109), nhop

c     for now (12/21/91) keep upper case
c     if(pline(1:2) .ne. 'C*') then
      if(upstrg(pline(1:2)) .ne. 'C*') then
        pline(58:58) = upstrg(pline(58:58))
        pline(65:65) = upstrg(pline(65:65))
        pline(84:84) = upstrg(pline(84:84))
      endif
c
      write(arcout, '(a)') pline(1:lentru(pline))
      if(style .ne. 'final') n15 = n15 + 1
      goto 310
c
c***********************genterate perfect data*********************
  302 k = keyph(i)
      read(phcard(i), 3021) ldate, lhr, lmin
3021  format(9x, i6.6, 2i2)
      psec = p(k) - x(4,k)
      if(ldx(k) .ne. 0) ssec = s(k) - x(4,ldx(k))
      call tshift(ldate,lhr,lmin,psec,ssec)
      call formf(psec, ipsec, 5, 2)
      if(ldx(k) .ne. 0) then
c       take care of s phase
        call formf(ssec, issec, 5, 2)
        write(arcout, 304) phcard(i)(1:9), ldate, lhr, lmin, ipsec,
     *  issec, phcard(i)(37:40)
        if(style .ne. 'final') n15 = n15 + 1
  304   format(a9, i6.6, 2i2, i5, 7x, i5, a4)
      else
c       only have p to write out
        write(arcout, 305) phcard(i)(1:9), ldate, lhr, lmin, ipsec
  305   format(a9, i6.6, 2i2, i5)
        if(style .ne. 'final') n15 = n15 + 1
      endif
  310 continue
      if(ipun .eq. 4) then
        write(arcout,320)
  320   format('                      ')
        if(style .ne. 'final') n15 = n15 + 1
      endif
      return
      end
c end npunch
c openfl.for    [pc]
      subroutine openfl(iunit, ifile, istat, izero, ishr, iform, irecl)
c system dependent program to open files for hypoelllipse
      integer iunit
c             iunit is unit number for this file
      character*(*) ifile
c                   ifile is name of file to be opened
      character*(*) istat
c                   istat is open status ('new' or 'old' or 'scratch'
c                     or 'unknown')
      character*(*) izero
c                   izero must be 'zero' or 'null'; if 'zero' then
c                     open with "blank='zero'"
      character*(*) ishr
c                   ishr indicates share status:
c                     must be 'readonly' or 'none' ('shared' is not coded)
      character*(*) iform
c                   iform = 'noprint' for some files on masscomp
c                     uacal file needs 'unformatted'
c                     otherwise = 'none'
      integer*4     irecl
c                   irecl is used as record length on vax, unless
c                     it is zero.
c
      if (ishr .ne. 'readonly' .and. ishr .ne.
     &       'none') then
        print *, 'openfl: invalid argument for ishr: ', ishr
	print *, 'iunit = ', iunit
	print *, 'ifile = ', ifile
	print *, 'istat = ', istat
	print *, 'izero = ', izero
	print *, 'ishr  = ', ishr
	print *, 'irecl = ', irecl

        stop
      endif
      if (izero .ne. 'zero' .and. izero .ne. 'null') then
        print *, 'openfl: invalid argument for izero: ', izero
        stop
      endif
      if (istat .ne. 'new' .and. istat .ne. 'old' .and.
     &    istat .ne. 'scratch' .and. istat .ne. 'unknown') then
        print *, 'openfl: invalid argument for istat: ', istat
        stop
      endif
c* (vax
c      if (irecl .le. 0) then
c        if (ishr .eq. 'readonly') then
c          if (izero .eq. 'zero') then
c            open(unit=iunit, file=ifile, status=istat,
c     &      blank='zero', readonly)
c          else
c            open(unit=iunit, file=ifile, status=istat,
c     &      readonly)
c          endif
c        else
c          if (izero .eq. 'zero') then
c            open(unit=iunit, file=ifile, status=istat,
c     &      blank='zero')
c          else
c            open(unit=iunit, file=ifile, status=istat)
c          endif
c        endif
c      else
c        if (ishr .eq. 'readonly') then
c          if (izero .eq. 'zero') then
c            open(unit=iunit, file=ifile, status=istat,
c     &      blank='zero', readonly, recl=irecl)
c          else
c            open(unit=iunit, file=ifile, status=istat,
c     &                    readonly, recl=irecl)
c          endif
c        else
c          if (izero .eq. 'zero') then
c            open(unit=iunit, file=ifile, status=istat,
c     &      blank='zero', recl=irecl)
c          else
c            open(unit=iunit, file=ifile, status=istat,
c     &      recl=irecl)
c          endif
c        endif
c      endif
c* vax)
c* (unix
cc#if masscomp
cc      if (iform .eq. 'none') then
cc#else
c      if (iform .eq. 'none' .or. iform .eq. 'noprint') then
cc#endif
c          if (izero .eq. 'zero') then
c            open(unit=iunit, file=ifile, status=istat,
c     &      blank='zero')
c          else
c            open(unit=iunit, file=ifile, status=istat)
c          endif
c      else
c          if (izero .eq. 'zero') then
c            open(unit=iunit, file=ifile, status=istat,
c     &      blank='zero', form=iform)
c          else
c            open(unit=iunit, file=ifile, status=istat,
c     &      form=iform)
c          endif
c      endif
cc     write (0,'(" trying to open ",a," with status = ",a)')
cc    &  ifile, istat
c* unix)
c* (pc
      if (istat .eq. 'scratch') istat = 'new'
      if (irecl .le. 0) then
          if (izero .eq. 'zero') then
            open(unit=iunit, file=ifile, status=istat,
     &      blank='zero')
          else
            open(unit=iunit, file=ifile, status=istat)
          endif
      else
          if (izero .eq. 'zero') then
            open(unit=iunit, file=ifile, status=istat,
     &      blank='zero', recl=irecl)
          else
            open(unit=iunit, file=ifile, status=istat,
     &      recl=irecl)
          endif
      endif
c* pc)
      return
      end
c end openfl
c opfls.for    [pc]
c ******************************************************************
c ***          subroutine opfls of program hypoellipse          ***
c ******************************************************************
c ***   a rewritten version of the opfls subroutine from usgs.   ***
c ***  this version works satisfactorily when run interactively  ***
c ***    or with its input "redirected" from a file or shell     ***
c ***               script.  ghcs at uaf, 5/17/88.               ***
c ******************************************************************
      subroutine opfls (iplt, inmain, root, nttab)
c* (pc
$notruncate
c* pc)
c
c     for the masscomp compiler, open statements must have form='noprint'
c     for all files to be written, except the print and log files.
c --- the iplt and inmain arguments will be set to return the
c --- filecode from which the main input file should be read.
c --- that is, either from a separate file (unit 8), or the
c --- standard input stream (unit 5) which may be a redirected
c --- c-shell "here document" containing shell variables.
c --- the ilis argument will be set to return the initial, or
c --- default value of ilis, a new option variable controlling
c --- the verboseness of the "printer" output file (unit 9).
c --- it provides a mechanism for generating an output file
c --- which contains only the information vital to evaluating
c --- the "correctness" of a phaselist and location.
c --- 7/29/89 note, jcl:
c      ilis is now set by the constants print option (see 2.2.3.20 in
c      the hypoellipse manual.
      integer iplt, inmain
c --- lengths of filename and type strings and the input
c --- and output unit numbers for this routine.
      integer max_length, max_type, in, out
      parameter (max_length = 255)
      parameter (max_type = 4)
      parameter (in = 5)
      parameter (out = 6)
c --- logical unit numbers for the files to be opened.
      integer stdin_unit, main_unit, stderr_unit, punt
      integer print_unit, summary_unit, arch_unit, tb11_unit
      integer tb12_unit, tb13_unit, stacorr_unit, log_unit
      integer tmp1_unit, tmp2_unit, tmpsum_unit, tmparc_unit
      integer tmpmess_unit
c* (vax
c* (pc
      parameter (stderr_unit  = 6)
c* pc)
c* vax)
c* (unix
c      parameter (stderr_unit  = 0)
c* unix)
      parameter (stdin_unit   = 5)
      parameter (main_unit    = 8)
      parameter (print_unit   = 9)
      parameter (summary_unit = 4)
      parameter (arch_unit    = 11)
      parameter (tb11_unit    = 21)
      parameter (tb12_unit    = 22)
      parameter (tb13_unit    = 23)
      parameter (stacorr_unit = 13)
      parameter (log_unit     = 6)
      parameter (tmp1_unit    = 2)
      parameter (tmp2_unit    = 3)
      parameter (tmpsum_unit  = 14)
      parameter (tmparc_unit  = 15)
      parameter (tmpmess_unit = 16)
c --- default input file and root/base strings.
      character*(max_length) in_def, root_def
      parameter (in_def = 'stdin')
      parameter (root_def = 'hypoe')
c --- the default file types (extensions) for output files.
      character*(max_type) print_type, summary_type, arch_type
      character*(max_type) stacorr_type,log_type, tmp1_type
      character*(max_type) tmp2_type, tmpsum_type, tmparc_type
      character*(max_type) tmpmess_type
      parameter (print_type   = '.out')
      parameter (summary_type = '.sum')
      parameter (arch_type    = '.arc')
      parameter (stacorr_type = '.nst')
      parameter (log_type     = '.log')
      parameter (tmp1_type    = '.1st')
      parameter (tmp2_type    = '.2st')
      parameter (tmpsum_type  = '.3sc')
      parameter (tmparc_type  = '.4sc')
      parameter (tmpmess_type = '.5sc')
c --- the lentru function returns the number of
c --- non-blank characters in its string argument.
c --- it returns zero if the string is completely blank.
      integer lentru
c --- local variable declarations.
      character*(*) root
      character*(max_length) fname, default, log_file
      character*(max_length) dnstrg
      character answer*1, openstat*7
      integer status, length
      integer length_def, length_root, length_log
      common /punt/ punt
      default = in_def
      length_def = lentru(default)
      write(out, '(a)') ' begin hypoellipse'
      write (out,
     *  '(/'' hypoellipse input filename ('',a,'')? '',$)')
     *  default(1:length_def)
      read (in, '(a)', iostat=status) fname
      length = lentru(fname)
      if ((status.ne.0) .or. (length.eq.0)) then
        fname = default
        length = length_def
      endif
      write (out, '(1x,a)') fname(1:length)
      if (dnstrg(fname(1:length)).eq.'stdin') then
        inmain = stdin_unit
        iplt = stdin_unit
      else
        inmain = main_unit
        iplt = main_unit
c       call openfl(iunit, ifile, istat,  izero,   ishr,
c    *          iform, irecl)
        call openfl( iplt, fname, 'old', 'zero', 'readonly',
     *          'none', 0)
c       open (unit=iplt, file=fname(1:length), status='old',
c    *  blank='zero')
c ---   previous version used blank='zero' which is very troublesome
c ---   (e.g. integers must be right justified in their fields
c ---   which is not obvious when looking at an input file).
c ---   however, some summary records have imbedded blanks with the
c ---   yrmodyhrmn field, and will not be read correctly unless
c ---   blank = 'zero' is specified!
c ---   also, iostat=status gives less information than the operating
c ---   system when the open fails.
c       if (status.ne.0) then
c         write (out,
c    *      '('' hypoellipse failed to open input file "'',a,''"'')')
c    *      fname(1:length)
c         stop 'abort from module opfls.'
c       endif
      endif
c --- should the output files be opened with a status of
c --- "new" to prohibit overwriting existing ones of the
c --- same name, or "unknown" to allow overwriting?
      openstat = 'new'
c* (unix
c      write (out,
c     * '('' allow overwriting of existing output files (no)? '',$)')
c      read (in, '(a)', iostat=status) answer
c      answer = dnstrg( answer )
c      if ((status.ne.0) .or. (answer.eq.' ')) answer = 'n'
c      if (answer.eq.'y') then
c        answer = 'y'
c        openstat = 'unknown'
c      else
c        answer = 'n'
c        openstat = 'new'
c      endif
c      write (out, '(1x,a)') answer
c* unix)
c --- output filename's root or basename used in the default
c --- names for the remaining files.
      default = root_def
      length_def = lentru(default)
      write (out, '('' root for output filenames ('',a,'')? '',$)')
     *  default(1:length_def)
      read (in, '(a)', iostat=status) root
      length_root = lentru(root)
      if ((status.ne.0) .or. (length_root.eq.0)) then
        root = root_def
      endif
      length_root = lentru(root)
c --- a kludge to fix the file names not passed as arguments.
c     if (root(length_root:length_root) .eq. 'p') then
c       root = root(1:length_root-1)
c       length_root = length_root - 1
c     endif
      write (out, '(1x,a)') root(1:length_root)
c --- summary of warning messages and final summary/log filename
      default = root(1:length_root)//log_type
      length_def = lentru(default)
c* (vax
c* (pc
      write (out, '('' log filename or "screen" ('',a,'')? '',$)')
c* pc)
c* vax)
c* (unix
c      write (out, '('' log filename or "stdout" ('',a,'')? '',$)')
c* unix)
     *  default(1:length_def)
      read (in, '(a)', iostat=status) log_file
      length_log = lentru(log_file)
      if ((status.ne.0) .or. (length_log.eq.0)) then
        log_file = default
        length_log = length_def
      endif
      write (out, '(1x,a)') log_file(1:length_log)
c --- we can't actually open the log file now, since we're still
c --- in the process of using unit 6 interactively...
c
c --- printer output file --------------------------------------
      punt = print_unit
      default = root(1:length_root)//print_type
      length_def = lentru(default)
c* (unix
c      write (out,
c     * '('' filename for printer output or "screen" ('',
c     * a,'')? '',$)')
c     *  default(1:length_def)
c      read (in, '(a)', iostat=status) fname
c* unix)
c* (vax
c* (pc
      write (out, '(a)') ' filename for printer output or "log"'
      write (out,
     * '('' answer "log" to combine with log file ('',
     * a,'')? '',$)')
     *  default(1:length_def)
      read (in, '(a)', iostat=status) fname
      fname = dnstrg( fname )
      if (fname .eq. 'log') fname = 'screen'
c* pc)
c* vax)
      length = lentru(fname)
      if ((status.ne.0) .or. (length.eq.0)) then
        fname = default
        length = length_def
      endif
      write (out, '(1x,a)') fname(1:length)
      if (dnstrg(fname(1:length)) .eq. 'screen') then
        punt = stderr_unit
      else
c       call openfl(     iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
        call openfl(punt, fname, openstat, 'null', 'none',
     *         'none', 0)
      endif
c     open (print_unit, file=fname(1:length), status=openstat)
c
c --- summary record filename --------------------------------------
      write (out,
     * '('' will this job generate a summary file (no)? '',$)')
      read (in, '(a)', iostat=status) answer
      if ((status.ne.0) .or. (answer.eq.' ')) answer = 'n'
      write (out, '(1x,a)') answer
      answer = dnstrg( answer )
      if (answer .eq. 'y') then
        default = root(1:length_root)//summary_type
        length_def = lentru(default)
        write (out,
     *    '('' summary output filename ('',a,'')? '',$)')
     *    default(1:length_def)
        read (in, '(a)', iostat=status) fname
        length = lentru(fname)
        if ((status.ne.0) .or. (length.eq.0)) then
          fname = default
          length = length_def
        endif
        write (out, '(1x,a)') fname(1:length)
c       call openfl(       iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
        call openfl(summary_unit, fname, openstat, 'null', 'none',
     *      'noprint', 0)
c       open (summary_unit,file=fname(1:length),status=openstat)
c    *        form = 'noprint')
      endif
c
c --- archive phase filename --------------------------------------
      write (out,
     *  '('' will this job generate an archive file (no)? '',$)')
      read (in, '(a)', iostat=status) answer
      if ((status.ne.0) .or. (answer.eq.' ')) answer = 'n'
      write (out, '(1x,a)') answer
      answer = dnstrg( answer )
      if (answer .eq. 'y') then
        default = root(1:length_root)//arch_type
        length_def = lentru(default)
        write (out, '('' archive output filename ('',a'')? '',$)')
     *    default(1:length_def)
        read (in, '(a)', iostat=status) fname
        length = lentru(fname)
        if ((status.ne.0) .or. (length.eq.0)) then
          fname = default
          length = length_def
        endif
        write (out, '(1x,a)') fname(1:length)
c       call openfl(    iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
        call openfl(arch_unit, fname, openstat, 'null', 'none',
     *      'noprint', 0)
c       open (arch_unit, file=fname, status=openstat)
c     *        form = 'noprint')
      endif
c
c --- travel time table filenames --------------------------------------
      write (out, '('' use travel time tables (no)? '',$)')
      read (in, '(a)', iostat=status) answer
      if ((status.ne.0) .or. (answer.eq.' ')) answer = 'n'
      write (out, '(1x,a)') answer
      answer = dnstrg( answer )
      if (answer .eq. 'y') then
        default = 'none'
        length_def = lentru(default)
        write (out, '('' filename for first table ('',a,'')? '',$)')
     *    default(1:length_def)
        read (in, '(a)', iostat=status) fname
        length = lentru(fname)
        if ((status.ne.0) .or. (length.eq.0)) then
          fname = default
          length = length_def
        endif
        write (out, '(1x,a)') fname(1:length)
c       call openfl(    iunit, ifile, istat,  izero,   ishr,
c    *          iform, irecl)
        call openfl(tb11_unit, fname, 'old', 'null', 'readonly',
     *         'none', 0)
c       open (tb11_unit, file=fname(1:length), status='old')
	nttab = 1
        default = 'none'
        length_def = lentru(default)
        write (out, '('' filename for second table ('',a,'')? '',$)')
     *    default(1:length_def)
        read (in, '(a)', iostat=status) fname
        length = lentru(fname)
        if ((status.ne.0) .or. (length.eq.0)) then
          fname = default
          length = length_def
        endif
        write (out, '(1x,a)') fname(1:length)
        if (dnstrg(fname(1:4)) .ne. 'none') then
	  nttab = 2
          call openfl(tb12_unit, fname, 'old', 'null', 'readonly',
     *      'none', 0)
c         call openfl(    iunit, ifile, istat,  izero,   ishr,
c    *      iform, irecl)
	endif
        default = 'none'
        length_def = lentru(default)
        write (out, '('' filename for third table ('',a,'')? '',$)')
     *    default(1:length_def)
        read (in, '(a)', iostat=status) fname
        length = lentru(fname)
        if ((status.ne.0) .or. (length.eq.0)) then
          fname = default
          length = length_def
        endif
        write (out, '(1x,a)') fname(1:length)
        if (dnstrg(fname(1:4)) .ne. 'none') then
          call openfl(tb13_unit, fname, 'old', 'null', 'readonly',
     *      'none', 0)
c         call openfl(    iunit, ifile, istat,  izero,   ishr,
c    *      iform, irecl)
	  if(nttab .ne. 2) then
	    write (out, '(a)') 'You may not define tables 11 and 13'
	    write (out, '(a)') 'without defining travel time table 12!'
	    stop
	  endif
	  nttab = 3
	endif
      endif
c
c --- optional relocation of events with revised delays --------------
      write (out, '('' does this job relocate events with revised '',
     *  ''station delays (no)? '',$)')
      read (in, '(a)', iostat=status) answer
      if ((status.ne.0) .or. (answer.eq.' ')) answer = 'n'
      write (out, '(1x,a)') answer
      answer = dnstrg( answer )
      if (answer .eq. 'y') then
        default = root(1:length_root)//stacorr_type
        length_def = lentru(root)
        write (out,
     *    '('' station delay output filename ('',a,'')? '',$)')
     *    default(1:length_def)
        read (in, '(a)', iostat=status) fname
        length = lentru(fname)
        if ((status.ne.0) .or. (length.eq.0)) then
          fname = default
          length = length_def
        endif
        write (out, '(1x,a)') fname(1:length)
c       open (stacorr_unit, file=fname(1:length), status=openstat)
c    *        form = 'noprint')
c       call openfl(       iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
        call openfl(stacorr_unit, fname, openstat, 'null', 'none',
     *      'noprint', 0)
      endif
c --- write out a blank record --------------------------------------
      write (out, *)
c --- we're done with interactive i/o using stdout (unit 6) so
c --- we can finally open the "log" file, if there is to be one.
c     fname = dnstrg( log_file(1:length_log) )
      fname = log_file(1:length_log)
      if((fname .ne. 'stdout') .and.
     *   (fname .ne. 'nolog')  .and.
     *   (fname .ne. 'screen'))
     *  call openfl(log_unit, log_file(1:length_log), openstat, 'null',
     *   'none', 'none', 0)
c       call openfl(   iunit,                  ifile,    istat,  izero,
c         ishr,   iform, irecl)
c    *  open (log_unit,file=log_file(1:length_log),status=openstat)
c --- use of status='scratch' is not standard and causes the files
c --- to be opened with "unique" names, which are uninformative.
c --- furthermore, it's sometimes useful to keep the temporary files
c --- when the program aborts.  therefore, you may wich to handle the
c --- "automatic" deletion another way (external to hypoellipse).
      write (out, *) 'root name = ', root(1:length_root)
c --- station scratch file number 1
      fname = root(1:length_root)//tmp1_type
c     call openfl(    iunit, ifile,     istat,  izero,   ishr,
c    *          iform, irecl)
      call openfl(tmp1_unit, fname, 'scratch', 'null', 'none',
     *      'noprint', 0)
c     open (tmp1_unit, file=root(1:length_root)//tmp1_type,
c    *      status='scratch')
c    *      form = 'noprint')
c --- station scratch file number 2
      fname = root(1:length_root)//tmp2_type
c     call openfl(    iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
      call openfl(tmp2_unit, fname, 'scratch', 'null', 'none',
     *      'noprint', 0)
c     open (tmp2_unit, file=root(1:length_root)//tmp2_type,
c    *      status='scratch')
c    *      form = 'noprint')
c --- temporary summary storage and uamag storage.  length must
c --- be 300 for uamag records.
      fname = root(1:length_root)//tmpsum_type
c     call openfl(       iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
      call openfl(tmpsum_unit, fname, 'scratch', 'null', 'none',
     *      'noprint', 300)
c     open (tmpsum_unit, file=root(1:length_root)//tmpsum_type,
c    *      status='scratch', recl=300)
c    *      form = 'noprint')
c --- temporary archive storage
      fname = root(1:length_root)//tmparc_type
c     call openfl(      iunit, ifile,     istat,  izero,   ishr,
c    *          iform, irecl)
      call openfl(tmparc_unit, fname, 'scratch', 'null', 'none',
     *      'noprint', 0)
c     open (tmparc_unit, file=root(1:length_root)//tmparc_type,
c    *      status='scratch')
c    *      form = 'noprint')
c --- temporary message storage
      fname = root(1:length_root)//tmpmess_type
c     call openfl(       iunit, ifile,    istat,  izero,   ishr,
c    *          iform, irecl)
      call openfl(tmpmess_unit, fname, 'scratch', 'null', 'none',
     *      'noprint', 0)
c     open (tmpmess_unit, file=root(1:length_root)//tmpmess_type,
c    *      status='scratch')
c    *      form = 'noprint')
      return
      end
c end opfls
c output.for    []
      subroutine output(kp)
c---- output hypocenter information
      include 'params.inc'
      parameter (ndly = 11)
      real bat2,bon2,latep,lonep,mag
      logical good, supout
c                 comp - z, n, or e component
      character*1 comp, dnstrg
c                 pha  - 3 letter phase description
      character*3 pha
      character*1 krmo, krm4, krm5, iprmk, iqs, iqd
      character erout*120, isorp*3
      character*133 fline, aline, bline, fmit*6, fmit1*6, fmit2*6, fmit3*24
      character*4 iahead*60, msta*5, nsta*5, icard*110
      integer punt
      common /dhip/ inpt,isa,ilis,inmain,injump
      common /punt/ punt
      common /char/ iahead, msta(npa), nsta(nsn), icard
      common /amo/ tempg(npa), ntgap
      common /anox/ keyd(npa)
      character*4 ipro, ichec, evtype*1, evstat*1, root*256
      common /dbhip/ ipro, ichec, evtype, evstat, root
      common /dbiln/ ioldq
      character*1 iclass
      common /dbio/ iclass(0:4)
      common /dhio/ nedit,good
      common /dio/ rmsmx,presmx,sresmx,noutmx,nimx,iprun,semx
      common /dmost/ ipun,ivlr,blank
      common /gmost/ az(npa)
      real*8 time1, time2
      common /hop/ time1,time2,nopu,notim
      logical freor
      common /hopq/ savla,savlo,savez,savor,freor
      common /hpn/ wslope
      common /ilpx/ calr(5, nsn),fmgc(nsn),xmgc(nsn)
      common /imost/ test(100)
      common /imost1/ dly(ndly,nsn),iprn,kno,klas(5, nsn)
      common /idno/ ksel,ksort
      common /ihfgpq/ itest(100)
      common /ilotv/ elvdly(npa)
      common /in/ irmo,iryr,odmag
      character*1  magke
      common /in1/ magke
      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 /iot/ flt(2,nsn),thks(nsn)
      common /iox/ prr(nsn),iuses
      integer fmwt, xmwt
      common /ioxl/ fmwt(nsn), xmwt(nsn)
      common /logfil/ logfil
      common /it/ vpvs(lmax2),vpvsm(mmax+3),nttab
      common /omnfh/ dmin,dmin3,sminp
      common /onx/ kluse(npa), cluse(npa), peruse(npa), ampuse(npa)
      common /ox/ sysmag(npa), gndmot(npa)
      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 /ofnv/ kmin
      character*1 nq
      common /ofnv1/ nq
      common /ohq/ gap, supout
      common /ohbl/ jav
      common /ofln/ sec
      common /olnx/ xmag(npa),fmag(npa),avxm,avfm,xmmed,fmmed
      common /on/ ix,iy,iz,ax1,axz
      character iqax*1
      common /on1/ iqax
      common /oqr/ y(4),kz,at(3),tl(3),pdrms,b(4)
      common /orz/ onf
      common /pfnoqv/ kdate,khrmn,khr
      character msym*1
      common /pfo/ msym(npa)
      common /pgoq/ s(npa)
      common /pgnoqv/ p(npa),jmin(npa)
      character*1 kwr, iqdo
      common /pbnoq/ kwr(npa),iqdo(npa)
      common /pgnov/ dt(npa)
      common /phoqn/ inst,knst
      common /pno/ lph,keyph(npa),nsum
      character phcard*117, krm*2, seq*5
      common /pno1/ phcard(npa), krm(npa), seq
      common /pmost/ nr,nrp,lbastm,tp(npa),ksmp(npa),kdx(npa)
      common /pnoqtx/ ldx(npa)
      character*4 krmp
      common /pnoqv/ krmp(npa)
      character*4 krms
      common /po/ krms(npa)
      common /pot/ ts(npa),ihr,model(npa), keyphi(npa)
      common /povx/ amx(npa)
      common /pox/ fmp(npa),prx(npa),icent1,icent2
      common /pt/ nlay(npa)
      common /qmost/ wt(npa),z
      common /qmost1/ lonep,ni,latep
      common /qgo/ org
      common /qgnotx/ delta(npa)
      common /qo/ dmax,ndec,adjsq,iph
      common /rbno/ idip(3),iaaz(3),a(3,3)
      common /rfgnoq/ se(4)
      common /ro/ yse,seorg,phi
      common /rioq/ been,damp,dmpinc,igo
      common /rob/ v(4,4),noaxp
      common /rq/ xmean(4),avwt,aveaj
      common /tmost/ x(4,npa)
      common /tonxb/ t(npa),fms(npa)
      common /xo/ nm,sdxm,nf,sdfm
      character*1 kodsym
      common /xo1/ kodsym(npa)
      common /xfmno/ mag
      common /zmost/ nrwt,rms,nswt,avr,aar,nsmp
      dimension demp(npa)
      dimension pruse(npa),sruse(npa),keyi(npa)
      dimension vv(2,2), ael(3), vec(2,2), eval(2), key(npa)
      data dum/1.0/
      if(iprn .ge. 5) write(punt, '(a)') 'begin subroutine output'
      iqax = ' '
      noaxp = 0
      rpd = 1.74533e-2
c---- initialize some remarks and calculate magnitudes if necessary
      krm1 = ' '
      krmo = ' '
c---- onf is sum of p and s weights. if onf = 0. put * next to origin time
      if(onf .eq. 0.) krmo = '*'
      krm2 = ' '
c---- kz is the number of fixed component in regress
      if(kz .eq. 3) krm2 = '*'
c---- define ix, iy, and iz in order of dip
      ix = 1
      iy = 2
      iz = 3
      if(idip(1) .le. idip(2)) goto 600
      ix = 2
      iy = 1
  600 if(idip(iy) .le. idip(3)) goto 610
      iz = iy
      iy = 3
  610 if(idip(ix) .le. idip(iy)) goto 620
      k = ix
      ix = iy
      iy = k
  620 continue
c---- convert lat and long to degrees and minutes
c---- convert lat and long to degrees and minutes
      if(iprn .ge. 5) then
        write(punt, '(a)') ' convert lat and long to deg & min'
        write(punt, *) ' latep, lonep = ', latep, lonep
      endif
      call unfold2(latep,lonep,lat1,isnla,bat2,lon1,isnlo,bon2)
      if(iprn .ge. 5) write(punt, *) lat1, isnla, bat2, lon1,
     *                               isnlo, bon2
c---- calculate time
      khr = ihr
      kmin = lbastm
      sec = org
      kkdate = kdate
c     write(punt, '(a, 2i10, f10.2, 2i10)')
c    *  'ihr, lbastm, org, kdate, kkdate',
c    *  ihr, lbastm, org, kdate, kkdate
      call tshift(kkdate,khr,kmin,sec,dum)
      adj = sqrt(adjsq)
      if((se(ix) .eq. blank) .or. (se(iy) .eq. blank)) then
        eqh = 10.0
      else
        eqh=sqrt(se(ix)**2+se(iy)**2)
      endif
c---- calculate gap
  431 j=0
      do 10 i = 1,nrp
        if (wt(i) .gt. 0.) then
c         continue on because there is p data
        else
c         check for s data
          if (ldx(i) .eq. 0) then
c           no s data
            goto 10
          else
c           if the s data has zero weight, skip it
            if (wt(ldx(i)) .eq. 0.) goto 10
          endif
        endif
        j=j+1
        tempg(j)=az(i)
   10 continue
      ntgap = j
      gap = 360.
      if(j .lt. 2) goto 22
      call sort(tempg,key,j)
      gap=tempg(1)+360.-tempg(j)
      do 20 i=2,j
      dtemp=tempg(i)-tempg(i-1)
      if(dtemp .gt. gap) gap=dtemp
   20 continue
   22 igap = gap + 0.5001
c---- calculate minimum distance
c---- sruse and pruse are arrays of s and p residuals used in solution
c---- nwtout is number of readings weighted out by the program
      iminp = 1
      dminp = 999.
      imins = 1
      dmins = 999.
      sruse(1) = 0.0
      pruse(1) = 0.0
      nwtout = 0
      iwt = 0
      ipwt = 0
      iswt = 0
      do 26 i = 1,nrp
        ianywt = 0
        kk = ldx(i)
        if(wt(i) .eq. 0.0) goto 23
        ianywt = 1
        ipwt = ipwt + 1
        if(delta(i) .lt. dminp) then
          dminp = delta(i)
          iminp = i
        endif
        pruse(ipwt) = abs(x(4,i))
   23   if (kwr(i) .ne. ' ' .and. dnstrg(kwr(i)) .ne. 'g')
     *   nwtout = nwtout + 1
        if(kk .eq. 0) goto 24
        if (kwr(kk) .ne. ' ' .and. dnstrg(kwr(i)) .ne. 'g')
     *   nwtout = nwtout + 1
        if(wt(kk) .eq. 0.0)  goto 24
        ianywt = 1
        iswt = iswt + 1
        if(delta(i) .lt. dmins) then
          dmins = delta(i)
          imins = i
        endif
        sruse(iswt) = abs(x(4,kk))
   24   if(ianywt .eq. 0) goto 26
        iwt = iwt + 1
        demp(iwt) = delta(i)
   26 continue
      dmin = dminp
      if(dmins .lt. dmin) dmin = dmins
      dmin3=dmin
      if(iwt .gt. 1) call sort(demp,key,iwt)
      if(iwt .ge. 2) dmin3 = demp(2)
      if(iwt .ge. 3) dmin3 = demp(3)
c---- compute s minus p time for closest station used in solution
      sminp = blank
      if((ipwt .eq. 0) .or. (iswt .eq. 0)) goto 950
c---- dminp/dmins is distance to closest station with p/s used
c---- iminp/imins is index of the closest stations with p/s
      if(msta(iminp) .ne. msta(imins)) goto 950
      sminp = s(imins) - p(iminp)
      if((sminp .gt. 99.9) .or. (sminp .lt. 0.)) sminp = 99.9
  950 if(ipwt .eq. 0) ipwt = 1
      if(iswt .eq. 0) iswt = 1
c---- calculate magnitude
      mag = blank
      if((iprn .ge. 2) .or. (kp .eq. 1)) call xfmags
      rmag = mag
c---- calculate q, iqs, and iqd
      ofd=z
      tfd=2.*z
      if(ofd .lt. 5.) ofd=5.
      if(tfd .lt. 10.) tfd=10.
      js=4
      if((rms.lt.0.50).and.(eqh.le.5.0)) js=3
      if(se(iz) .eq. blank) goto 28
      if((rms.lt.0.30).and.(eqh.le.2.5).and.(se(iz).le.5.0)) js=2
      if((rms.lt.0.15).and.(eqh.le.1.0).and.(se(iz).le.2.0)) js=1
   28 jd = 4
      if(nrwt .lt. 6) goto 30
      if((gap.le.180.).and.(dmin.le.50.)) jd=3
      if((gap.le.135.).and.(dmin.le.tfd)) jd=2
      if((gap.le. 90.).and.(dmin.le.ofd)) jd=1
   30 jav=(js+jd+1)/2
      nq=iclass(jav)
      iqs=iclass(js)
      iqd=iclass(jd)
      idmin=dmin + .5
c---- write heading and check if eq is out of order
      time2=1.d+0*sec+1.d+02*kmin+1.d+04*khr+1.d+06*kkdate
      if(iprn .lt. 0) goto 720
      if((iprn .eq. 0) .or. (inst .eq. 9)) goto 52
      if(ni .ne. 1) goto 60
      if(ndec .ge. 1) goto 60
   52 continue
      if((icent2 .eq. icent1 .and. time2 .lt. time1 - 20.d0) .or.
     *  (icent2 .lt. icent1)) then
        write(punt,54) icent2, kkdate, khr, kmin, sec
        write(logfil,54) icent2, kkdate, khr, kmin, sec
   54   format(' xxxxx this event is out of order xxxxx', /,
     *  i2, i6.6, 1x, i2, ':', i2, 1x, f6.2)
      endif
c---- write step output
   60 if(iprn .eq. 0) goto 720
      if(iph .eq. 1) goto 62
c---- step line  output
      write(punt,61)
   61 format(/,84x,
     *'( adjustments )(adjst.  taken)', /,
     *'  i  lat long depth    rms   prms    damp     no ',
     *'  phi  --------eigenvalues---------',
     *' dlat dlon   dz dlat dlon   dz')
      if(iprn .eq. 1) iph=1
c
   62 continue
      aline = ' '
c      write(aline, 162) ni,bat2,bon2,z-test(8),krm2,rms,nrwt
c  162 format(1x ,i2, f5.1, f5.1, f6.1, a1, f7.4, i3)
      write(aline, 162) ni,bat2,bon2,z-test(8),krm2
  162 format(1x ,i2, f5.1, f5.1, f6.1, a1)
c
      call formit(rms, rrms, fmit, 6, 1)
      write(aline(22:27), fmit) rrms

      if (pdrms .ne. blank) then
        call formit(pdrms, rpdrms, fmit, 6, 1)
        write(aline(29:34), fmit) rpdrms
      endif
c
      write(aline(36:42), '(e7.2)') damp
      write(aline(46:84), '(i3, 1x, f5.1, 3(1x, e9.4))')
     *  nrwt, phi, a(1,1), a(2,2), a(3,3)
c
      if (b(2) .ne. blank) then
        call formit(b(2), rb, fmit, 4, 1)
        write(aline(86:89), fmit) rb
      endif
c
      if (b(1) .ne. blank) then
        call formit(b(1), rb, fmit, 4, 1)
        write(aline(91:94), fmit) rb
      endif
c
      if (b(3) .ne. blank) then
        call formit(b(3), rb, fmit, 4, 1)
        write(aline(96:99), fmit) rb
      endif
c
      call formit(y(2), rb, fmit, 4, 1)
      write(aline(101:104), fmit) rb
      call formit(y(1), rb, fmit, 4, 1)
      write(aline(106:109), fmit) rb
      call formit(y(3), rb, fmit, 4, 1)
      write(aline(111:114), fmit) rb

      if( (ni .eq. 1) .or. (ni .eq. itest(37)) )
     *  write(punt, 64) lat1, lon1
   64 format(3x, 2i5)
c
c     step output
      write(punt, '(a132)') aline
c
      if((test(41) .eq. 1.0) .and. (iabs(ipun) .eq. 1))
     *call npunch('final')
      if((kp .eq. 0) .and. (iprn .lt. 3)) goto 100
c---- calculate single variable standard deviations
c---- first calculate shadow on surface
c-------- vv(1,1)*x**2 + vv(1,2)*x*y + vv(2,2)*y**2 = 1
  720 if(v(3,3) .lt. 0.000009) v(3,3) = 0.000009
      vv(1, 1) = (v(1,1)-v(1,3)**2/v(3,3))
      vv(1, 2) =    (v(1,2)-v(1,3)*v(2,3)/v(3,3))
      vv(2, 1) = vv(1, 2)
      vv(2, 2) = (v(2,2)-v(2,3)**2/v(3,3))
c---- compute principal axis
c---- find eigenvalues [eval(2)] and eigenvectors [vec(2,2)] for the
c     upper left nxn portion of vv(2, 2)
      call eigen1(ael, 2, 2, 2, vec, eval, vv, 0.0)
c
      if(eval(2) .lt. 0.000009) then
        ax1 = 99.
      else
        ax1 = yse/sqrt(eval(2))
        if(ax1 .gt. 99.) ax1 = 99.
      endif
c
      if(eval(1) .lt. 0.000009) then
        ax2 = 99.
      else
        ax2 = yse/sqrt(eval(1))
        if(ax2 .gt. 99.) ax2 = 99.
      endif
c
      aze2 = atan2(-vec(1,1),vec(2,1))/rpd
      aze1 = atan2(-vec(1,2),vec(2,2))/rpd
c
c---- compute maximum in z (standard error of z)
      if((v(1,1) .lt. 0.000009) .or. (nrwt .le. 3)) then
        axz = 99.
      else
        dnom = v(2,2) - v(1,2)*v(1,2)/v(1,1)
        if(dnom .lt. 0.000009) then
          axz = 99.
        else
          prob =  v(3,3) - v(1,3)**2/v(1,1)
     *    -( v(2,3) - v(1,3)*v(1,2)/v(1,1) )**2/dnom
          if(prob .lt. 0.000009) then
            axz = 99.
          else
            axz = yse/sqrt(prob)
            if(axz .gt. 99.) axz = 99.
          endif
        endif
      endif
c---- calculate quality based upon ax1 and ax2
      jq = 4
      bigst = axz
      if(ax1 .gt. axz) bigst = ax1
      if(bigst .le. 5.35) jq = 3
      if(bigst .le. 2.67) jq = 2
      if(bigst .le. 1.34) jq = 1
      iqax = iclass(jq)
      if(ioldq .eq. 0) jav = jq
      if(iprn .lt. 0) return
      if((ksel.ne.1) .and. ((nedit.eq.1).or.(nedit.eq.2))) goto 70

c write out az, dip, step and se for the principal directions
      write(punt, '(9x, a)')
     * '-az/dp--step---se =az/dp==step===se -az/dp--step---se'
      aline = ' '
      write(aline(10:15), '(i3, 1h/ ,i2)') iaaz(ix),idip(ix)
c
      if (at(ix) .ne. blank) then
        call formit(at(ix), rat, fmit, 5, 1)
        write(aline(17:21), fmit) rat
        call formit(se(ix), rse, fmit, 4, 1)
        write(aline(23:26), fmit) rse
      endif
c
      write(aline(28:33), '(i3, 1h/ ,i2)') iaaz(iy),idip(iy)
c
      if (at(iy) .ne. blank) then
        call formit(at(iy), rat, fmit, 5, 1)
        write(aline(35:39), fmit) rat
        call formit(se(iy), rse, fmit, 4, 1)
        write(aline(41:44), fmit) rse
      endif
c
      write(aline(46:51), '(i3, 1h/ ,i2)') iaaz(iz),idip(iz)
c
      if (at(iz) .ne. blank) then
        call formit(at(iz), rat, fmit, 5, 1)
        write(aline(53:57), fmit) rat
        call formit(se(iz), rse, fmit, 4, 1)
        write(aline(59:62), fmit) rse
      endif
      write(punt, '(a)') aline

c write out comment records
      do 530 i = 1,lph
        if(keyph(i) .eq. -1) write(punt,520) phcard(i)(1:115)
  520   format(1x,a)
  530 continue

      if(iprn .ge. 0)
     * write(punt,65) ax2,ax1,axz,iqax,aze2,aze1
   65 format(/1x,'horizontal and vertical single variable standard',
     *   ' deviations (68% - one degree of freedom; max 99 km)'/,
     *  7x,6hseh = ,f6.2,13x,6hseh = ,f6.2,13x,6hsez = ,f6.2,
     * 13h   quality = ,a1/7x,6haz  = ,f6.0,13x,6haz  = ,f6.0/)
      write(punt,78) seorg,ni,dmax,seq
   78 format(' se of orig = ',f6.2,'; # of iterations = ',
     *       i3,'; dmax = ',f10.2,'; sequence number = ',a5)
      if(sminp .eq. blank) then
        write(punt, 975) evtype, evstat
  975   format(' event type = "', a1, '"',
     *           '; processing status = "', a1, '"', /,
     *           ' closest station did not use both p and s')
      else
        write(punt, 980) evtype, evstat, sminp
  980   format(' event type = "', a1, '"',
     *           '; processing status = "', a1, '"', /
     *    ' s minus p interval for closest station = ', f10.2)
      endif
c
c---- check for debug event
   70 continue
      if(nedit .ne. 0) then
        nopu = 1
        good = .true.
        call sort(pruse,keyi,ipwt)
        call sort(sruse,keyi,iswt)
        if(bigst .gt. semx) goto 107
        if(ni .gt. nimx) goto 107
        if(rms .gt. rmsmx) goto 107
        if(sruse(iswt) .gt. sresmx) goto 107
        if(pruse(ipwt) .gt. presmx) goto 107
        if(nwtout .gt. noutmx) goto 107
c----   event is good
        if(iabs(nedit) .eq. 1) then
          write(punt,74) iqax
   74     format('  quality based on standard errors is ',a1/)
          goto 73
        endif
        if(iabs(nedit) .eq. 3) goto 73
        noaxp = 1
        return
  107   good = .false.
        write(punt,76) ni,rms,pruse(ipwt),sruse(iswt),nwtout,bigst
   76   format(1h0,' d e b u g  e v e n t     ni   rms presmx sresmx',
     *  ' nwtout  bigst   ',/,
     *  24x,                              i5, f6.1,  f7.1,  f7.1,
     *  i7,  f7.1)
        if(nedit .lt. 0) nopu=0
        if(iabs(nedit) .eq. 3) goto 73
      endif
   73 continue
c
c---- check for preferred magnitude
      if((magke .ne. ' ') .and.
     * (dnstrg(magke) .ne. 'x') .and.
     * (dnstrg(magke) .ne. 'f') .and.
     * (dnstrg(magke) .ne. 'a')) then
        write(punt,77) odmag,magke
   77   format('  preferred magnitude on summary record will be: ',
     *  f5.2,1x,a1)
      endif
c
c---- set up final output line
      write(punt,75)
   75 format(/'    date    origin      lat      long    depth    mag',
     * ' no d1 gap d  rms    avwt   se')
      fline = ' '
      write(fline, 71) icent2, kkdate, krmo, khr, kmin, sec,
     *  lat1, isnla, bat2, lon1, isnlo, bon2, krm1, z-test(8), krm2
      jnst = knst*10 + inst
   71 format(1x, i2, i6.6, a1, 2i2, f6.2, i3, a1, f5.2, i4, a1,
     *       f5.2, a1, f6.2, a1)
c---- add magnitude
      if (rmag .ne. blank) then
        call formit(mag, rmag, fmit, 5, 1)
        write(fline(49:53), fmit) rmag
      endif
c----- add up through rms
      write(fline(54:84), 82) nrwt, idmin, igap, kno, rms, avwt, yse
   82 format(2i3, i4, i2, f7.4, f6.2, f6.2)
      write(punt, '(a)') fline(1:84)
      write(punt, '(t22, f7.4, 2x, f8.4)') lat1+bat2/60.,lon1+bon2/60.
      fline = ' '
c
      write(punt, 821)
821   format('    seh  sez q sqd  adj in nr   avr  aar nm',
     * ' avxm mdxm sdxm nf avfm mdfm sdfm   vpvs')
      write(fline(1:29), 822) ax1, axz, nq, iqs, iqd, adj, jnst, nr
822   format( 1x, f6.1, f5.1, 3(1x,a1), f5.2, 2i3)
c----- add avr, aar, and nm
      call formit(avr, avrlm, fmit1, 5, 1)
      call formit(aar, aarlm, fmit2, 4, 1)
      fmit3 = fmit1(1:5)//',1x,'//fmit2(2:5)//',i3)'
      write(fline(31:44), fmit3) avrlm, aarlm, nm
c----- add average xmag
      if (avxm .ne. blank) then
        call formit(avxm, ravxm, fmit, 4, 1)
        write(fline(45:48), fmit) ravxm
      endif
c----- add median xmag
      if (xmmed .ne. blank) then
        call formit(xmmed, ravxm, fmit, 4, 1)
        write(fline(50:53), fmit) ravxm
      endif
c----- add two more
      write(fline(54:62), 83) sdxm,nf
   83 format(f5.1, i3)
c----- add average fmag
      if (avfm .ne. blank) then
        call formit(avfm, ravfm, fmit, 4, 1)
        write(fline(63:66), fmit)ravfm
      endif
c----- add median fmag
      if (fmmed .ne. blank) then
        call formit(fmmed, ravfm, fmit, 4, 1)
        write(fline(68:71), fmit)ravfm
      endif
c----- add rest
      write(fline(72:83), 84) sdfm, wslope
   84 format(f5.1, f7.3)
c
      write(punt, '(a)') fline(1:83)
c
c
c---- write out warning messages
      rewind(16)
      read(16, 97) erout
   96 read(16, 97, end=98) erout
        if(erout(1:3) .eq. 'end') goto 98
        write(punt, 97) erout
        write(logfil, 97) erout
   97   format(a)
      goto 96
   98 continue
c
c for now, remove short-circuit of detailed printout
c     if(nedit .eq. 0) goto 100
c     if( .not. good) goto 108
c     if(iprn .ge. 1) goto 108
c     return
c
  100 if(kp .eq. 1) goto 108
      if(iprn .le. 1) return
  108 isorp = 'fmp'
      if(iuses .eq. 1) isorp = 'fms'
      write(punt,110)
  110 format (/'                      -- travel times and delays --', /,
     *'  stn c pha remk p p-sec s-sec resid  std-er   dist  azm ain',
     *'    tc c vthk  ttob-ttcal-dlay-edly=resid rmk stn pha sources')
czzz
c***     1         2         3         4         5         6         7         8
c***5678901234567890123456789012345678901234567890123456789012345678901234567890
c stn c pha remk p p-sec s-sec resid  std-er   dist  azm ain    tc c vthk  ttob-
cbrse     1 ipu3 d199.25       -0.99 * 1.011  199.3  241 180  -.26 1 4.00 49.32
cbrse   s 1  s 2         49.25 -0.99 r 1.011             180  -.26 1 4.00 49.32
c spu     2 ip+9 -  3.22                      551.1  123  82  -.26 1 2.00 49.32
c spu   s 2  s 2         13.22                                -.26 1 2.00 49.32
c spu   smp                     1.55   1.222                              14.00
c
      if(nrp .eq. 0) return
      do 500 i=1,nrp
      aline = ' '
      k=i
      if(ksort .eq. 0) k = keyd(i)
      kk = ldx(k)
c
c set up component
c     if(dnstrg(krmp(k)(2:2)) .eq. 'p') then
c       comp = 'z'
c       comp = ' '
c     else
c       comp = krmp(k)(2:2)
c     endif
      comp = msta(k)(5:5)
c set up p-phase identifier
c     pha = 'p  '
      pha = '   '
      if(nlay(k) .ne. 0) then
        write(pha(3:3), '(i1)') nlay(k)
      endif
c
c      write(aline(1:24), 91) msta(k),  comp,   pha, krmp(k), msym(k),
c     1 p(k)
      write(aline(1:24), 91) msta(k)(1:4),  comp,   pha, krmp(k),
     *  msym(k)
91    format(                  1x,a4, 1x,a1, 1x,a3,   1x,a4,   1x,a1)
c add p arrival time
      aline(20:24) = phcard(keyphi(k))(20:24)
c
      pres = x(4,k)
      apres = abs(pres)
      iprmk = kwr(k)
      if(iprun .ne. 1) then
        if(wt(k) .gt. 0.) then
c fix up p residuals for routine processing
          if(iprmk .eq. ' ') then
            if((apres.gt.0.6).and.(i.le.5)) iprmk = '*'
            if((apres.gt.0.9).and.(delta(k).lt.150.)) iprmk = '*'
            if((apres.gt.1.5).and.(delta(k).lt.350.)) iprmk = '*'
          endif
        endif
        if(apres .lt. 2.25) then
          l = (apres + 0.25)/0.5
          pres = l*0.5
        endif
      endif
      if(pres .gt. 999.99) pres = 999.99
      if(pres .lt. -99.99) pres = -99.99
c add p-residual
      if(wt(k) .eq. 0.0) then
        write(aline(31:44), '(f6.2, 1x, a1, a)') pres, iprmk, ' -----'
      else
        setmp = yse/sqrt(wt(k))
        if(setmp .gt. 99.99) setmp = 99.99
        write(aline(31:44), '(f6.2, 1x, a1, f6.2)') pres, iprmk, setmp
      endif
c
      write(aline(45:52), 93) delta(k), iqdo(k)
93    format(                     f7.1,      a1)
c
      call formf(az(k), iaz, 3, 0)
      write(aline(53:56), '(i4)') iaz
c
      call formf(ain(k), iain, 4, 0)
      write(aline(57:60), '(i4)') iain
c add time correction
      if (dt(k) .ne. 0.) then
        dtk = dt(k)
        call formit(dtk, dtk, fmit, 5, 1)
        write(aline(62:66), fmit) dtk
      endif
c add model
      write(aline(67:68), '(i2)') model(k)
c leave one blank
c add variable layer thickness
      if (thks(k) .ne. blank) then
        call formit(thks(k), rthks, fmit, 4, 1)
        write(aline(70:73), fmit) rthks
      endif
c add p-travel time observed
      tpk = tp(k) - org
      if(tpk .lt. -3500.) tpk = tpk + 3600.
      write(aline(74:85), 94) tpk, t(k)
94    format(2f6.2, i2)
c leave one blank
c add p-delay
      kji=kdx(k)
      if (dly(kno, kji) .ne. 0.) then
        call formit(dly(kno,kji), dlyk, fmit, 4, 1)
        write(aline(87:90), fmit) dlyk
      endif
c leave one blank
c add elevation delay
      if (elvdly(k) .ne. 0.) then
        call formit(elvdly(k), eldly, fmit, 4, 1)
        write(aline(92:95), fmit) eldly
      endif
c repeate p-residual and add remark
      write(aline(96:104), '(f6.2, 1x, a2)') pres, krm(k)
c leave one blank
c repeate station name
      aline(106:109) = msta(k)(1:4)
c leave one blank
c repeate phase type
      aline(111:113) = aline(9:11)
c leave two blanks
c give sources and number of hops from phase record
      aline(116:121) = phcard(keyphi(k))(105:110)
c
c aline is now set up correctly for a normal p phase
      if(ldx(k) .eq. 0) then
        if(ksmp(k) .eq. 1) then
c no s data - write p record and go on to next phase
          write(punt, '(a)') aline(1:121)
          goto 500
        else
c s minus p data -- remove residual and std-er
c set up bline as a revised p record and write out
98764     bline = aline
          bline(25:45) = ' '
c compute the p-residual
          pres = tp(k) - t(k) - org - dly(kno, kji) - elvdly(k)
          if(pres .gt. 999.99) pres = 999.99
          if(pres .lt. -99.99) pres = -99.99
          write(bline(96:101), '(f6.2)') pres
c write out the p line for an s-p
          write(punt, '(a)') bline(1:121)
c
c next set up bline as an s record
          bline(9:9) = 's'
          write(bline(13:16), '(a4)') krms(k)
c remove p seconds
          bline(20:24) = ' '
c add s seconds
          bline(26:30) = phcard(keyphi(k))(32:36)
c remove dist, azm, and ain
          bline(47:60) = ' '
c add observed and computed travel times
          tsk = ts(k) - org
          if(tsk .lt. -3500.) tsk = tsk + 3600.
          write(bline(74:86), '(2f6.2)') tsk, vpvsm(model(k))*t(k)
c add s-delay
          if (sdly(kno, kji) .ne. 0.) then
            call formit(sdly(kno,kji), sdlyk, fmit, 4, 1)
            write(bline(87:90), fmit) sdlyk
          else
            bline(87:90) = ' '
          endif
c add elevation delay
          if (elvdly(k) .ne. 0.) then
            call formit(vpvsm(model(k))*elvdly(k), eldly, fmit, 4, 1)
            write(bline(92:95), fmit) eldly
          else
            bline(92:95) = '    '
          endif
c compute the s-residual
98765     sres = ts(k) - org - vpvsm(model(k))*(t(k) + elvdly(k))
     *    - sdly(kno, kji)
          if(sres .gt. 999.99) sres = 999.99
          if(sres .lt. -99.99) sres = -99.99
          write(bline(96:101), '(f6.2)') sres
          bline(111:111) = 's'
          write(punt, '(a)') bline(1:113)
c
c finally set up aline as the s minus p record
          aline(9:30) = 'smp'
          aline(45:73) = ' '
c add observed and computed s-p
          csmp = (vpvsm(model(k))-1.)*(t(k) + elvdly(k))
     *    - dly(kno, kji) + sdly(kno, kji)
          write(aline(74:85), '(2f6.2)') ts(k)-tp(k), csmp
          aline(86:104) = ' '
c repeate phase type
          aline(111:113) = aline(9:11)
          write(punt, '(a)') aline(1:113)
        endif
      else
c
c this is the case of a normal p and s
c first write p record
          write(punt, '(a)') aline(1:121)
c
c next convert aline into an s record
        aline(9:9) = 's'
        aline(13:16) = krms(k)
c remove corrected first motion and p-res
        aline(18:24) = ' '
c add s seconds
        aline(25:30) = phcard(keyphi(k))(32:36)
        sres = x(4,kk)
        krm5 = kwr(kk)
        if(iprun .ne. 1) then
c fix up s residuals for routine processing
          asres = abs(sres)
          if(wt(kk) .gt. 0.0) then
            if(krm5 .eq. ' ') then
              if((asres.gt.0.9).and.(i.le.5)) krm5 = '*'
              if((asres.gt.1.5).and.(delta(k).lt.150.)) krm5 = '*'
              if((asres.gt.2.0).and.(delta(k).lt.350.)) krm5 = '*'
            endif
          endif
          if(asres .lt. 2.25) then
            l = (asres + 0.25)/0.5
            sres = l*0.5
          endif
        endif
        if(sres .gt. 999.99) sres = 999.99
        if(sres .lt. -99.99) sres = -99.99
        if(wt(kk) .eq. 0.0) then
          write(aline(31:44), '(f6.2, 1x, a1, a)') sres, krm5, ' -----'
        else
          setmp = yse/sqrt(wt(kk))
          if(setmp .gt. 99.99) setmp = 99.99
          write(aline(31:44), '(f6.2, 1x, a1, f6.2)') sres, krm5, setmp
        endif
c remove dist and azm
        aline(47:56) = ' '
        call formf(ain(kk), iain, 4, 0)
        write(aline(57:60), '(i4)') iain
        write(aline(67:68), '(i2)') model(kk)
c add observed and computed travel times
        tsk = ts(k) - org
        if(tsk .lt. -3500.) tsk = tsk + 3600.
        write(aline(74:86), '(2f6.2)') tsk, t(kk)
c add s-delay
        if (sdly(kno, kji) .ne. 0.) then
          call formit(sdly(kno,kji), sdlyk, fmit, 4, 1)
          write(aline(87:90), fmit) sdlyk
        else
          aline(87:90) = ' '
        endif
c add elevation delay
        if (elvdly(kk) .ne. 0.) then
          call formit(elvdly(kk), eldly, fmit, 4, 1)
          write(aline(92:95), fmit) eldly
        else
          aline(92:95) = '    '
        endif
c repeate s-residual and leave remark
        write(aline(96:101), '(f6.2)') sres
c repeate phase type
        aline(111:113) = aline(9:11)
        write(punt, '(a)') aline(1:113)
      endif
500   continue
c
c loop through again for magnitude data
c
      write(punt, 550)
550   format(/,
     * '                                     -- magnitude data --',/,
     * '  stn c source sys    c10   amx gr ink     amf    per ',
     * '  unit/mm  gnd mot u xmgc xmag  fmp fmag')
      do 900 i = 1, nrp
        aline = ' '
        k=i
        if(ksort .eq. 0) k = keyd(i)
        kji=kdx(k)
        if (amx(k) .ne. 0.) then
          write(aline(16:17), '(i2)') kluse(k)
          aline(28:31) = phcard(keyphi(k))(44:47)
c add corrected amplitude
          write(aline(39:47), '(f9.0)') ampuse(k)
        endif
        if (cluse(k) .gt. 0.) then
          write(aline(18:25), '(f8.2)') cluse(k)
	endif
	if ((sysmag(k) .ne. blank) .and. (gndmot(k) .ne. blank)) then
          write(aline(55:74), '(1pe10.4, 1x, 1pe9.3)')
     *    sysmag(k), gndmot(k)
        endif
c a1vco gain range state 0, 1, or 2
        aline(34:34) = phcard(keyphi(k))(62:62)
c semens playback state 0 or 1
        aline(37:37) = phcard(keyphi(k))(61:61)
c add xmag
        if (xmag(k) .ne. blank) then
          write(aline(75:80), '(f5.2, 1x)') xmgc(kji)
          call formit(xmag(k), rxmag, fmit, 4, 1)
          write(aline(81:84), fmit) rxmag
        endif
c flag deviant xmags
        if(xmag(k) .ne. blank) then
	  if(xmwt(kji) .eq. 0) then
c flag excluded xmags
c 	    aline(65:65) = 'e'
	    aline(85:85) = 'e'
	  else if(abs(xmag(k) - avxm) .ge. 0.5) then
c flag deviant xmags
c           aline(65:65) = '*'
            aline(85:85) = '*'
	  endif
        endif
c add f-p
        rfmp = fmp(k)
        if(iuses .eq. 1) rfmp = fms(k)
        if (rfmp .ne. blank) then
          call riorbk(rfmp, ifmp, fmit, 4, 0)
c         write(aline(66:69), fmit) ifmp
          write(aline(86:89), fmit) ifmp
        endif
c leave 2 blanks
c add fmag
        rfmag = fmag(k)
        if (fmag(k) .ne. blank) then
          call formit(fmag(k), rfmag, fmit, 4, 1)
c         write(aline(71:74), fmit) rfmag
          write(aline(91:94), fmit) rfmag
        endif
c
        if(dnstrg(kodsym(k)) .eq. 's') then
c fmp too short with respect to s-p, so fmag not computed
          krm4 = kodsym(k)
        else
          krm4 = ' '
          if(fmag(k) .ne. blank) then
            if(fmwt(kji) .eq. 0) then
c flag excluded fmags
              krm4 = 'e'
            else if(abs(fmag(k) - avfm) .ge. 0.5) then
c flag deviant fmags with krm4
              krm4 = '*'
            endif
          endif
        endif
c       aline(75:75) = krm4
        aline(95:95) = krm4
	if(aline(1:33) .ne. ' ')
     *    write(aline(49:53), '(f5.2)') peruse(k)
        if((aline(1:33) .ne. ' ') .or.
     *    (aline(38:133) .ne. ' ')) then
c set up component
c         if(dnstrg(krmp(k)(2:2)) .eq. 'p') then
c           comp = 'z'
c           comp = ' '
c         else
c           comp = krmp(k)(2:2)
c         endif
	  comp = msta(k)(5:5)
          aline(2:5) = msta(k)(1:4)
          aline(7:7) = comp
c add amplitude source
          aline(11:11) = phcard(keyphi(k))(108:108)
c         write(punt, '(a)') aline(1:75)
          write(punt, '(a)') aline(1:95)
        endif
900   continue
      return
      end
c end output
c phagt.for    [pc]
      subroutine phagt(phcard, keyph, lph, inpt, injump, eoff, icent2,
     *  itest)
c read the data associated with one earthquake
c
c     phcard(npa)     array of records read
c     keyph(npa)      key to record type
c                       0 --- phase record
c                      -1 --- comment record
c                      -2 --- summary record
c                      -3 --- instruction record
c                             (begins with more, dist, or 4 blanks)
c                      -4 --- bad time
c                      -5 --- reset record
c                      -6 --- save record - no longer used
c                      -7 --- rerun record - no longer used
c                      -8 --- scatter record
c                      -9 --- not on station list
c                     -10 --- jump record
c                     -11 --- decode error
c                     -12 --- sleep record
c                     -13 --- stop record
c                     -14 --- nongaus record
c                     -15 --- master scatter record
c     lph             number of records returned
c     inpt            unit number for reading data
c     eoff            .false. --- end of file not reached
c                     .true. --- found end of file
c     npa            maximum number of phases allowed
c
      save onesav
      include 'params.inc'
      character*10 begtim, dnstrg*4, pname*4, foocard*115, acent2*2
      integer punt, icentsave, icent2, itest(100)
      common /punt/ punt
      common /logfil/ logfil
      common /dmost/ ipun,ivlr,blank
      character*117 phcard(npa), fname*50
      dimension keyph(npa)
      logical onesav, eoff, fsteq
      data onesav/.false./,fsteq/.true./
      eoff = .false.
      icentsave = 0
      npurep = 0
c
      lph = 1
c     write(punt, *) 'onesav = ', onesav
      if (onesav) then
        phcard(1) = phcard(isv)
        keyph(1) = keyph(isv)
        onesav = .false.
        if ((keyph(1) .ne. -2) .and. (keyph(1) .ne. 0)) then
          return
        endif
        goto 25
      endif
c
c loop from lph = 1 to npa
c     print *,'phaget about to read unit ',inpt
15    read(inpt, 20, end = 70) phcard(lph)
20    format(a)
c
c classify record
25    continue
      pname = dnstrg(phcard(lph)(1:4))
c
c comment record
      if (pname(1:2) .eq. 'c*') then
        if( (ipun .eq. 2 .or. ipun .eq. 3) .and.
     *    (lph .eq. 1) .and. (fsteq) ) then
c comment records prior to the first summary or phase record should
c be written directly to the archive phase file
          write(11, '(a)') phcard(lph)
          goto 15
        endif
        keyph(lph) = -1
        goto 50
c
c summary record in 115 column format
      else if ((phcard(lph)(81:81) .eq. '/') .or.
c* (vax
c* (pc
     *  (phcard(lph)(81:81) .eq. '\')) then
c* pc)
c* vax)
c* (unix
c     *  (phcard(lph)(81:81) .eq. '\\')) then
c* unix)
        keyph(lph) = -2
	icent2 = itest(55)
	write(phcard(lph)(116:117), '(i2)') icent2
	if(icentsave .ne. 0) then
c two summary records with different centuries will bomb out the program
	  if(icentsave .ne. icent2) then
	    write(punt, '(a,/)')
     *        'xxxerrorxxx in phagt.  Two summary records for the ',
     *        'save event with different centuries. ', phcard(lph)
	    write(logfil, '(a,/)')
     *        'xxxerrorxxx in phagt.  Two summary records for the ',
     *        'save event with different centuries. ', phcard(lph)
	    stop
	  endif
	endif
        icentsave = icent2

        if (npurep .gt. 0) goto 90
        goto 50
c
c summary record in 117 column format
      else if ((phcard(lph)(83:83) .eq. '/') .or.
c* (vax
c* (pc
     *  (phcard(lph)(83:83) .eq. '\')) then
c* pc)
c* vax)
c* (unix
c     *  (phcard(lph)(83:83) .eq. '\\')) then
c* unix)
        keyph(lph) = -2
c must be in new 117-column format, so switch back for internal use
        read(phcard(lph)(1:2), '(i2, t1, a2)') icent2, acent2
        foocard = phcard(lph)(3:117)
        phcard(lph) = foocard//acent2
        if(icentsave .ne. 0) then
c two summary records with different centuries will bomb out the program
          if(icentsave .ne. icent2) then
            write(punt, '(a,/)')
     *        'xxxerrorxxx in phagt.  Two summary records for the ',
     *        'save event with different centuries. ', phcard(lph)
            write(logfil, '(a,/)')
     *        'xxxerrorxxx in phagt.  Two summary records for the ',
     *        'save event with different centuries. ', phcard(lph)
	    stop
          endif
        endif
        icentsave = icent2

        if (npurep .gt. 0) goto 90
        goto 50
c
c instruction record
      else if ((pname .eq. '    ') .or.
     *        (pname .eq. 'more')) then
c    *   .or. (pname .eq. 'dist')) then
        keyph(lph) = -3
        goto 80
c
c reset record
      else if (pname .eq. 'rese') then
        keyph(lph) = -5
        if (lph .gt. 1) goto 90
        return
c
c save record
      else if (pname .eq. 'save') then
        keyph(lph) = -6
        if (lph .gt. 1) goto 90
        return
c
c rerun record
      else if (pname .eq. 'reru') then
        keyph(lph) = -7
        if (lph .gt. 1) goto 90
        return
c
c scatter record
      else if (pname .eq. 'scat') then
        keyph(lph) = -8
        if (lph .gt. 1) goto 90
        return
c
c nongaus record
      else if (pname .eq. 'nong') then
        keyph(lph) = -14
        if (lph .gt. 1) goto 90
        return
c
c master scatter record
      else if (pname .eq. 'mast') then
        keyph(lph) = -15
        if (lph .gt. 1) goto 90
        return
c
c sleep record
      else if (pname .eq. 'slee') then
        keyph(lph) = -12
        if (lph .gt. 1) goto 90
        return
c
c stop record
      else if (pname .eq. 'stop') then
        keyph(lph) = -13
        if (lph .gt. 1) goto 90
        return
c
c jump record
      else if (pname .eq. 'jump') then
        keyph(lph) = -10
        fname = phcard(lph)(6:55)
        begtim = phcard(lph)(56:65)
        if (inpt .eq. injump) then
          close (unit = injump, iostat = iostat)
          write(punt, 40)
          write(logfil, 40)
40        format(' xxxerrorxxx can not jump from a jump file.')
          stop
        else
          inpt = injump
          write(punt, 41) fname
          write(logfil, 41) fname
41        format(' jump to ', a)
          close (unit = injump, iostat = iostat)
c         call openfl( iunit, ifile, istat,  izero, ishr,
c    *    iform, irecl)
          call openfl(injump, fname, 'old', 'zero', 'readonly',
     *    'none', 0)
c          open(unit=injump, file=fname, blank='zero', status='old',
c     *    iostat = iostat)
c          if (iostat .ne. 0) goto 498
          if(begtim .ne. ' ') then
            write(punt, 412) begtim
            write(logfil, 412) begtim
412         format(' begin processing with event:  ', a)
            if(lph .ne. 1) then
              write(punt, 413)
              write(logfil, 413)
413           format(' *** jump record must be between events ***', /,
     *        ' *** instruction record missing, so stop ***')
              stop
            endif
            do 417 i = 10, 1, -1
              if(begtim(i:i) .ne. ' ') goto 418
417         continue
418         lenbt = i
420         read(injump, '(a)', end = 430) phcard(lph)
            if(phcard(lph)(1:lenbt) .eq. begtim(1:lenbt)) goto 25
            goto 420
430         write(punt, 435) begtim
            write(logfil, 435) begtim
435         format(' *** could not find ', a, ' so stop ***')
            stop
          endif
        endif
        if (lph .gt. 1) goto 90
        goto 15
c498     write(punt, 499) fname
c        write(logfil, 499) fname
c499     format(' *** jump file ', a, /, ' could not be opened',
c     *  ' *** so stop ***')
c        stop
      endif
c
c must be regular phase record
      keyph(lph) = 0
      npurep = npurep + 1
50    fsteq = .false.
      lph = lph + 1
      if (lph .le. npa) goto 15
c
55    write(logfil, 60) npa
      write(punt, 60) npa
60    format(///,' xxxxx exceeded maximum of ',i4,
     *           ' records per event, so stop!  xxxxx')
      stop
c
c come here on eof detected
70    eoff = .true.
c     write(punt, *) 'just set eoff to true in phagt.  eoff = ', eoff
      onesav = .false.
      if (lph .gt. 1) then
        phcard(lph) = ' '
        keyph(lph) = -3
        write(logfil, 72)
        write(punt, 72)
72      format(//, ' detected end of file prior to final instruction',
     *  ' record.', /, ' check input data for completeness.')
        return
      else
        lph = lph - 1
        return
      endif
c
c come here after an instruction record or other special record
80    continue
      return
c
c come here if missing an instruction record
c     ('add' blank instruction record)
90    isv = lph + 1
      phcard(isv) = phcard(lph)
      keyph(isv) = keyph(lph)
      phcard(lph) = ' '
      keyph(lph) = -3
      onesav = .true.
      write(logfil, 92) phcard(isv)
      write(punt, 92) phcard(isv)
92    format(/, ' xxxxx warning - missing instruction record - xxxxx',
     *       /, ' prior to this record:', /, ' ', a)
      return
      end
c end phagt
c phaind.for    []
      integer function phaind(msta,nsta,ns,itarg,ierr)
c     returns itarg, the index of nsta that equals target.
      character*5 msta, nsta(*)
      ierr = 0
      i = 1
      j = ns
   10 mid = (i + j)/2
      if (msta .eq. nsta(mid) ) then
        itarg = mid
        goto 50
      else
        if (msta .lt. nsta(mid) ) then
          j = mid - 1
          if (j .lt. i) goto 20
          goto 10
        else
          i = mid + 1
          if (i .gt. j) goto 20
          goto 10
        end if
      end if
   20 ierr = 1
   50 phaind = ierr
      return
      end
c
cnch      integer function phaind(msta,nsta,ns,i,ierr)
c     returns i, the index of nsta that equals msta.
c     modified so that one does not need to assume the stations are in order
c     9 June 1994 jas/vtso small corrections
c     6 September 2003: corrected character*4 to character*5
cnch      dimension nsta(*)
cnch      character*5 msta, nsta
chch      ierr = 0
cnch      do 40 i=1,ns
cnch      if (msta .eq. nsta(i)) goto 50
cnch   40 continue
cnch      ierr = 1
cnch   50 phaind = ierr
cnch      return
cnch      end
cnch
c end phaind