      subroutine st_sln(imagn,imevt,mcode,me_stn_adj_file)
c
c ...... write hypocentral solution files
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
c
      character*80 fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1             epi_dt,me_stn_adj_file
      character*9 epi_id
      common /fnamr/ fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1               epi_id,epi_dt
      common/files/ncr,nsp,ns1,ns2,ns3,ns4
      integer*4 imagn,imevt
      character*1 mcode
c
c ...... write relp solution file
c
      open(nsp,file=fnamsl(1:lnblnk(fnamsl)),status='unknown')
      call eqloc(imagn,imevt,mcode,me_stn_adj_file)
      close(nsp)
c
c ...... write relp event parameter file
c
      call eqevl(imagn,mcode)
c
      return
      end

      subroutine eqevl(imagn,mcode)
c
c ...... write event parameter file
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
c
      common /mags/ amla,seamla,amlb,seamlb,nmla,nmlb
      character*80 fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1             epi_dt
      character*9 epi_id
      common /fnamr/ fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1               epi_id,epi_dt
      character*80 geo_info
      common /ginfo/geo_info
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/dim/m,n,nmax,nit
      common/obs/pt(1000),st(1000),nsreq
      integer*4 imagn
      character*1 mcode
      common /eqelp/ del_ns,hor_se,iaz_gap,np_loc,ns_loc
	integer*4 evt_flag
	common /etyp3/ evt_flag
c
      open(26,file=fnamel(1:lnblnk(fnamel)),status='unknown')
      if(imagn.lt.1) then
        write(26,'(a2,".",a3,2x,a2,":",a2,3x,"ML = ",f3.1,3x,a)') 
     1       epi_id(1:2),epi_id(3:5),epi_id(6:7),epi_id(8:9),amla,
     2       geo_info(1:lnblnk(geo_info))
      else
        write(26,'(a2,".",a3,2x,a2,":",a2,14x,a)') 
     1         epi_id(1:2),epi_id(3:5),epi_id(6:7),epi_id(8:9),
     2         geo_info(1:lnblnk(geo_info))
      endif
      sqfact=0.
      if(cov(1,1).gt.0.) then
        seex=sqrt(cov(1,1))
        sqfact=sqfact+cov(1,1)
      else
        seex=0.
      endif
      if(cov(2,2).gt.0.) then
        seey=sqrt(cov(2,2))
        sqfact=sqfact+cov(2,2)
      else
        seey=0.
      endif
      if(cov(3,3).gt.0.) then
        seet=sqrt(cov(3,3))
        sqfact=sqfact+25.d0*cov(3,3)
      else
        seet=0.
      endif
      if(n.eq.3) then
        seez=0.
      else
        if(cov(4,4).gt.0.) then
          seez=sqrt(cov(4,4))
        sqfact=sqfact+0.25d0*cov(4,4)
        else
          seez=sqrt(cov(4,4))
        endif
      endif
c
c ...... cov_sma
c
	cov_sma=cov(1,1)+cov(2,2)+cov(3,3)+cov(4,4)
	csmaxis=dsqrt(cov_sma)
	if(csmaxis.le.0.00001d0) then
	  fl_csma=-5.0
	else
	  fl_csma=dlog10(csmaxis)
	endif
c
      ieth=et/3600.d0
      ietm=(et-3600.d0*dble(ieth))/60.d0
      fets=et-3600.d0*dble(ieth)-60.d0*dble(ietm)
      ieth=mod(ieth,24)
c      write(*,*) 'eqevl - et, ieth, ietm & fets = ',et,ieth,ietm,fets
      call icon5(fets,iets1,iets2,iets3,iets4,iets5)
      if(seet.gt.9.999d0) seet=9.999d0
      if(seey.gt.99.99d0) seey=99.99d0
      if(seex.gt.99.99d0) seex=99.99d0
      if((n.eq.4).and.(seez.gt.999.99d0)) seez=999.99d0
      if(n.eq.4) then
        write(26,'(8x,i2.2,":",i2.2,":",2i1,".",3i1,"  ",f5.3,
     1        3x,f8.4," ",f5.2,3x,f9.4," ",f5.2,
     2        3x,f6.2," ",f6.2)') ieth,ietm,iets1,iets2,iets3,iets4,
     3        iets5,seet,elatd,seey,elongd,seex,edepth,seez
      else
        write(26,'(8x,i2.2,":",i2.2,":",2i1,".",3i1,"  ",f5.3,
     1        3x,f8.4," ",f5.2,3x,f9.4," ",f5.2,
     2        3x,f6.2," (R)")') ieth,ietm,iets1,iets2,iets3,iets4,
     3        iets5,seet,elatd,seey,elongd,seex,edepth
      endif
      if(imagn.eq.1) go to 1
      if((nmla.eq.0).and.(nmlb.eq.0)) then
        go to 1
      elseif(nmlb.eq.0) then
        write(26,'(8x,"ML(bdsn) = ",f4.2," ",f5.3," ",
     1        i2,29x,"Model = ",a1)') 
     2        amla,seamla,nmla,mcode
      else
        write(26,'(8x,"ML(bdsn) = ",f4.2," ",f5.3," ",
     1        i2,3x,"ML(brk) = ",f4.2," ",f5.3," ",
     2        i2,3x,"Model = ",a1)') 
     3        amla,seamla,nmla,amlb,seamlb,nmlb,mcode
      endif
    1 continue
c
c ...... add sqfact - solution quality factor
c        add following parameters for REDI
c
c            del_ns  = distance in km to nearest station
c            hor_se  = location horizontal standard error in km
c            iaz_gap = azimuthal gap in degrees
c            np_loc  = number of phases used in location (P+S)
c
      if(sqfact.gt.0.) then
        nnn=5
        sqfact=dsqrtf(nnn,sqfact)
      endif
      if(sqfact.le.2.) then
        isqfact=0
      elseif(sqfact.le.5.) then
        isqfact=1
      elseif(sqfact.le.10.) then
        isqfact=2
      elseif(sqfact.le.20.) then
        isqfact=3
      else
        isqfact=4
      endif
c      if(np_loc.ge.ns_loc+2) ns_loc=ns_loc+1
      if(se.gt.99.999d0) se=99.999d0
      if(del_ns.gt.999.99d0) del_ns=999.99d0
      write(26,'(8x,"se = ",f6.3,"   nstn = ",i3,"   qual = ",i1,
     1      "   del_ns = ",f6.2)') 
     2      se,ns_loc,isqfact,del_ns
      if(hor_se.gt.99.999d0) hor_se=99.999d0
      write(26,'(8x,"hor_se = ",f6.3,"   az_gap = ",i3.3,
     1      "   npha = ",i4,"   evt_flag = ",i1,"   log_sma = ",
     2      f4.1)') hor_se,iaz_gap,np_loc,evt_flag,fl_csma
      close(26)
c
      return
      end

      subroutine eqloc(imagn,imevt,mcode,me_stn_adj_file)
c
c ...... output event location parameters
c
      implicit real*8 (a-h,o-z)
      character*1 polar,gfc1,gc1
      character*72 id
      character*4 stn,gfc4,gc4
      common/files/ncr,nsp,ns1,ns2,ns3,ns4
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/obs/pt(1000),st(1000),nsreq
      common/obs1/stn(1000),polar(1000)
      common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns
      common/model/ vpa,vpb,vpc,vsa,vsb,vsc,vd
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/condn/cn
      common/dim/m,n,nmax,nit
      common/nm/nitmax
      common/gfile/gf(1000,20)
      common/gfile1/gfc4(1000),gfc1(1000)
      common/psiv/psi(1000),psiscl
      common/prec/precis,cndn
      common/gapaz/gap1,gap,gap2
      common/info/sinfo(1000),sclinf
      common/time/eqhour,isyear,isdoy,ishour,ismin,fssec
      character*80 geo_info,me_stn_adj_file
      common /ginfo/geo_info
      common /mags/ amla,seamla,amlb,seamlb,nmla,nmlb
      character*80 fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1             epi_dt
      character*9 epi_id
      common /fnamr/ fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam,
     1               epi_id,epi_dt
      common /eqelp/ del_ns,hor_se,iaz_gap,np_loc,ns_loc
      real*8 g(1000)
      real*8 cor(4,4)
      integer*4 imagn,imevt
      character*1 mcode
c
      ext=ex
      eyt=ey
      j=0
      do 3 i=1,nsreq
      sy=syf(i)
      sx=sxf(i)
      psa=psaf(i)
      ssa=ssaf(i)
      call dazd(-ext,eyt,-sx,sy,delta,azim)
      do 1 k=4,20
      g(k)=0.d0
    1 continue
      gc4=stn(i)
      g(2)=delta
      g(3)=azim
      g(6)=-1.d0
      if(pt(i).gt.0.d0) then
c        pt(i)=pt(i)+60.d0*dreal(ismin)+fssec
        j=j+1
c        if(pt(i).lt.(et-600.d0)) pt(i)=pt(i)+3600.d0
        pti=pt(i)-dint(pt(i)/3600.d0)*3600.d0
        g(4)=idint(pti)/60
        g(5)=pti-60.d0*g(4)
        g(6)=w(j)
        g(7)=b(j)
        g(16)=psi(j)
        g(18)=sinfo(j)
        gc1=polar(j)
        g(7)=bmax(j)
        z0=0.d0
        err=1.d-3
        call ptcalc(vpa,vpb,vpc,vd,z0,ez,err,g(2),d1,d2,g(12),g(13),
     1              d3)
      endif
      g(10)=-1.d0
      if(st(i).gt.0.d0) then
c        st(i)=st(i)+60.d0*dreal(ismin)+fssec
        j=j+1
c        if(st(i).lt.(et-600.d0)) st(i)=st(i)+3600.d0
        sti=st(i)-dint(st(i)/3600.d0)*3600.d0
        g(8)=idint(sti)/60
        g(9)=sti-60.d0*g(8)
        g(10)=w(j)
        g(11)=b(j)
        g(17)=psi(j)
        g(19)=sinfo(j)
        g(11)=bmax(j)
        z0=0.
        err=1.d-3
        call ptcalc(vsa,vsb,vsc,vd,z0,ez,err,g(2),d1,d2,g(14),g(15),
     1              d3)
      endif
      do 2 k=2,19
      gf(i,k)=g(k)
    2 continue
      gfc1(i)=gc1
      gfc4(i)=gc4
    3 continue
c
      if(imagn.lt.1) then
c        if(amla.gt.0.) then
          write(nsp,'(a8," (",a2,".",a3,") ",a5," - ML ",f3.1," - ",a)') 
     1        epi_dt(1:8),epi_id(1:2),epi_id(3:5),epi_dt(10:14),amla,
     2        geo_info(1:lnblnk(geo_info))
c        else
c          write(nsp,'(a8," (",a2,".",a3,") ",a5," - ML ",f3.1," - ",a)') 
c     1        epi_dt(1:8),epi_id(1:2),epi_id(3:5),epi_dt(10:14),amla,
c     2        geo_info(1:lnblnk(geo_info))
c        endif
      else
        write(nsp,'(a8," (",a2,".",a3,") ",a5," - ",a)') 
     1        epi_dt(1:8),epi_id(1:2),epi_id(3:5),epi_dt(10:14),
     2        geo_info(1:lnblnk(geo_info))
      endif
c
c      write(nsp,'(a)') id(1:lnblnk(id))
c      write(5,'(a)') id(1:lnblnk(id))
      if(nit.gt.1) then
        if(nit.ge.nitmax) then
          write(nsp,6) nitmax
        else
          write(nsp,7) nit
        endif
      endif
    6 format(/,1x,'*dnc* did not converge within ',i2,
     1           ' iterations *****')
    7 format(/,1x,'*c* converged in ',i2,' iterations')
      write(nsp,8) se
    8 format(/,1x,'standard error = ',f7.3,
     1           '   number of observations')
      smaxis=0.d0
      do 9 i=1,n
      smaxis=smaxis+qi(i)*qi(i)
    9 continue
      if(smaxis.gt.0.d0) then
        nnn=6
        smaxis=se*dsqrtf(nnn,smaxis)
      else
        smaxis=0.d0
      endif
      call necomp(smaxis,eff,inp,ins,int)
      np_loc=inp+ins
      write(nsp,10) smaxis,inp,ins,int
   10 format(1x,'semimajor axis = ',f6.2,
     1'    p=',i3,', s=',i3,', & total=',i3)
      write(nsp,11) eff,precis
   11 format(1x,'efficiency     = ',f5.1,' %',
     1'   precision = ',e7.1)
      call azgap(nsreq)
      iaz_gap=gap
      igap1=gap1+0.5d0
      igap=gap
      igap2=gap2+0.5d0
      write(nsp,12) igap2,igap,igap1,cndn
   12 format(1x,'az gap = ',i3,' (',i3,') ',i3,
     1'     condition number = ',f8.2)
      write(nsp,13) (q(iqq),iqq=1,n)
   13 format(1x,'singular values = ',4f9.5)
      if(ez.le.vd) then
        zvp=vpa+vpb*ez
        zvs=vsa+vsb*ez
      else
        zvp=vpc
        zvs=vsc
      endif
      write(nsp,14) zvp,zvs,mcode
   14 format(1x,'vp = ',f5.2,' & vs = ',f5.2,
     1       ' at focal depth   (',a1,' model)')
      write(nsp,15) psiscl
   15 format(1x,'data covariance scale factor  = ',f6.3)
      write(nsp,16) sclinf
   16 format(1x,'data information scale factor = ',f6.3)
      arg = datan(1.d0)/45.d0
c      eth=eqhour
c      et=et+60.d0*dreal(ismin)+fssec
c      if(et.lt.0.d0) then
c        et=et+3600.d0
c        eth=eth-1.d0
c        if(eth.lt.0.d0) then
c          eth=23.d0
c        endif
c      elseif(et.gt.3600.d0) then
c        et=et-3600.d0
c        eth=eth+1.0d0
c        if(eth.gt.23.d0) then
c          eth=0.d0
c        endif
c      endif
c      etm=idint(et/60.d0)
c      ets=et-60.d0*etm
      edepth = ez
      ey=datand(dtand(ey)/0.993277d0)
      eyt=ey
      ext = ex
      if(cov(1,1).gt.0.d0) then
        nnn=7
        seex=dsqrtf(nnn,cov(1,1))
      else
        seex=0.d0
      endif
      if(cov(2,2).gt.0.d0) then
        nnn=8
        seey=dsqrtf(nnn,cov(2,2))
      else
        seey=0.d0
      endif
      hor_se=cov(1,1)+cov(2,2)+2.*dabs(cov(1,2))
      if(hor_se.gt.0.d0) then
        nnn=9
        hor_se=dsqrtf(nnn,hor_se)
      endif
      if(cov(3,3).gt.0.d0) then
        nnn=10
        seet=dsqrtf (nnn,cov(3,3))
      else
        seet=0.d0
      endif
      if(n.eq.3) then
        seez=0.d0
        cov(4,4)=0.d0
      else
        if(cov(4,4).gt.0.d0) then
          nnn=11
          seez=dsqrtf(nnn,cov(4,4))
        else
          seez=0.d0
        endif
      endif
      ieth=et/3600.d0
      ietm=(et-3600.d0*dble(ieth))/60.d0
      fets=et-3600.d0*dble(ieth)-60.d0*dble(ietm)
      ieth=mod(ieth,24)
c      write(*,*) 'eqloc - et, ieth, ietm & fets = ',et,ieth,ietm,fets
      call corr(cov,n,cor,nit)
      call icon(fets,iets1,iets2,iets3,iets4)
c   17 format(6i1,'.',2i1,2x,f7.3,2x,f8.3,2x,
c     1       f4.1,2x,i2,1x,i2,1x,f5.3,1x,f5.2,3(1x,f4.1))
      write(nsp,18)
   18 format(/,9x,'parameter',11x,'std err',8x,'correlation matrix')
      write(nsp,19) ieth,ietm,iets1,iets2,iets3,iets4,seet,
     1              (cor(1,iq),iq=1,n)
   19 format(1x,
     1'origin time = ',i2.2,'-',i2.2,'-',2i1,'.',2i1,3x,f6.3,
     2' s   ',4f6.2)
      write(nsp,20) ey,seey,(cor(2,iq),iq=1,n)
   20 format(1x,
     1'latitude    = ',f11.4,3x,f5.2,' km   ',4f6.2)
      write(nsp,21) ex,seex,(cor(3,iq),iq=1,n)
   21 format(1x,
     1'longitude   = ',f11.4,3x,f5.2,' km   ',4f6.2)
      if(n.eq.3) then
        write(nsp,22) edepth
      else
        write(nsp,23) edepth,seez,(cor(4,iq),iq=1,n)
      endif
   22 format(1x,
     1       'depth       = ',f8.2,' km (restrained)')
   23 format(1x,
     1       'depth       = ',f8.2, ' km   ',f5.2,' km   ',4f6.2)
      if(nit.ge.1) then
        write(nsp,24)
      else
        write(nsp,25)
      endif
   24 format(/,1x,'stn   delta  az p   time   h ai f  res   wt',
     1           '   cov info')
   25 format(/,1x,'stn   delta  az  p-time   p-res  s-time   s-res')
      ex=ext
      ey=eyt
      call sdprnt(imevt,me_stn_adj_file)
      write(2,26)
   26 format('******************************',
     1         '******************************',/)
c
      return
      end

      subroutine icon(x,ix1,ix2,ix3,ix4)
c
c ...... find individual digits in a four digit number
c
      implicit real*8 (a-h,o-z)
      data i1,i2,i3/10,100,1000/
      ix=idint(100.d0*x)
      ix1=ix/i3
      ix2=(ix-i3*ix1)/i2
      ix3=(ix-i3*ix1-i2*ix2)/i1
      ix4=(ix-i3*ix1-i2*ix2-i1*ix3)
c
      return
      end

      subroutine icon5(x,ix1,ix2,ix3,ix4,ix5)
c
c ...... find individual digits in a five digit number
c
      implicit real*8 (a-h,o-z)
      data i1,i2,i3,i4/10,100,1000,10000/
      ix=idint(1000.d0*x)
      ix1=ix/i4
      ix2=(ix-i4*ix1)/i3
      ix3=(ix-i4*ix1-i3*ix2)/i2
      ix4=(ix-i4*ix1-i3*ix2-i2*ix3)/i1
      ix5=(ix-i4*ix1-i3*ix2-i2*ix3-i1*ix4)
c
      return
      end

      subroutine azgap(n)
c
c ...... calculate azimuthal gap
c
c        find largest difference in azimuth among adjacent stations
c
      implicit real*8 (a-h,o-z)
      character*1 gfc1
      character*4 gfc4
      common/gfile/gf(1000,20)
      common/gfile1/gfc4(1000),gfc1(1000)
      common/gapaz/gap1,gap,gap2
      common/files/ ncr,nsp,ns1,ns2,ns3,ns4
      real*8 t(1000)
      do 1 i=1,n
      t(i)=gf(i,3)
c      if(t(i).lt.0.d0) then
c        t(i)=t(i)+360.d0
c      elseif(t(i).gt.360.d0) then
c        t(i)=t(i)-360.d0
c      endif
    1 continue
      nm1=n-1
      do 4 k=1,nm1
      k1=k+1
      do 3 i=k1,n
      if(t(i).gt.t(k)) go to 2
      temp=t(i)
      t(i)=t(k)
      t(k)=temp
    2 continue
    3 continue
    4 continue
      np1=n+1
      t(np1)=t(1)+360.d0
      gap1=t(1)
      gap=t(2)-t(1)
      gap2=t(2)
      do 5 j=2,np1
      k=j-1
      arg=t(j)-t(k)
      if(arg.ge.gap) then
        gap1=t(k)
        gap=arg
        gap2=t(j)
      endif
    5 continue
      if(gap2.gt.360.d0) gap2=gap2-360.d0
      gap1=360.d0-gap1
      gap2=360.d0-gap2
c
      return
      end

      subroutine necomp(sma,pcteff,inp,ins,int)
c
c ...... tally number of observations and compute efficiency
c
c        Efficiency is a measure of how well the given stations perform
c        compared to the theoretically optimal network geometry with the
c        same number of observations.  See above Uhrhammer reference.
c
      implicit real*8 (a-h,o-z)
      character*1 polar
      character*72 id
      character*4 stn
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/obs/pt(1000),st(1000),nsreq
      common/obs1/stn(1000),polar(1000)
      common/model/ vpa,vpb,vpc,vsa,vsb,vsc,vd
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/dim/m,n,nmax,nit
c
      znp=0.d0
      zns=0.d0
      do 1 i=1,nsreq
      if(pt(i).gt.0.d0) znp=znp+1.d0
      if(st(i).gt.0.d0) zns=zns+1.d0
    1 continue
      znt=znp+zns
      znp=idint(znp+0.1d0)
      zns=idint(zns+0.1d0)
      znt=idint(znt+0.1d0)
      inp=znp
      ins=zns
      int=znt
      zvp=vpa+vpb*edepth
      zvs=vsa+vsb*edepth
      znum=3.d0*dreal(n)-3.d0
      arg=0.5d0*znum*zvp*zvp*zvs*zvs/(znp*zvs*zvs+zns*zvp*zvp)
     1   +1.d0/(znt)
c ...... pcteff = efficiency in percent
      if(arg.gt.0.d0) then
        if(sma.ne.0.d0) then
          nnn=12
          pcteff=100.d0*dsqrtf(nnn,arg)*se/sma
          if(pcteff.gt.100.d0) pcteff=100.d0
        else
          pcteff=100.d0
        endif
      else
        pcteff=0.d0
      endif
c
      return
      end

      subroutine sdprnt(imevt,me_stn_adj_file)
c
c ...... output data for each station/observation
c
      implicit real*8 (a-h,o-z)
      integer*4 imevt
      character*80 me_stn_adj_file
      character*1 tc2,tc3,tc4,hl,hu,hp,hs,gfc1,rc1
      character*4 tc1,tc5,stnm,gfc4,rc4
      common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4),
     1            v(4,4),cov(4,4),x(4),bmax(1000),se
      common/obs/pt(1000),st(1000),nsreq
      common/files/ncr,nsp,ns1,ns2,ns3,ns4
      common/gfile/gf(1000,20)
      common/gfile1/gfc4(1000),gfc1(1000)
      common/miss/nmstn
      common/miss1/stnm(100)
      common/dim/m,n,nmax,nit
      common/time/eqhour,isyear,isdoy,ishour,ismin,fssec
      real*8 g(1000),fng(1000),r(1000)
      common /eqelp/ del_ns,hor_se,iaz_gap,np_loc,ns_loc
      data hu,hl/'U','L'/
c
      do 1 i=1,m
      g(i)=b(i)
      fng(i)=w(i)
      r(i)=bmax(i)
    1 continue
      do 2 i=1,nsreq
      r(2)=gf(i,2)
      g(i)=r(2)
      fng(i)=dreal(i)
    2 continue
      nm1=nsreq-1
      do 4 i=1,nm1
      i1=i+1
      do 3 j=i1,nsreq
      if(g(j).le.g(i)) then
        gt=g(i)
        fngt=fng(i)
        g(i)=g(j)
        fng(i)=fng(j)
        g(j)=gt
        fng(j)=fngt
      endif
    3 continue
    4 continue
c
      if(imevt.gt.0) then
c
c ...... open master event station adjustment file
c
        open(27,file=me_stn_adj_file(1:lnblnk(me_stn_adj_file)),
     1       status='unknown')
      endif
c
      do 10 i=1,nsreq
      ji=i
      nr=idint(fng(i)+0.5d0)
      do 6 j=1,nr
      do 5 jkl=2,19
      r(jkl)=gf(j,jkl)
    5 continue
      rc4=gfc4(j)
      rc1=gfc1(j)
    6 continue
      iaz=r(3)+0.5d0
      if(iaz.gt.360) iaz=0
      iaz=360-iaz
      iazh=iaz/100
      iazt=(iaz-iazh*100)/10
      iazo=iaz-iazh*100-iazt*10
c      if(r(4).gt.59.d0) then
c        r(4)=r(4)-60.d0
c      endif
      ipmt=r(4)/10.d0
      ipmo=r(4)-dreal(ipmt)*10.d0
      ip=r(5)*100.d0+0.5d0
      ips1=ip/1000
      ips2=(ip-ips1*1000)/100
      ips3=(ip-ips1*1000-ips2*100)/10
      ips4=ip-ips1*1000-ips2*100-ips3*10
      if(r(8).gt.59.d0) then
        r(8)=r(8)-60.d0
      endif
      ismt=r(8)/10.d0
      ismo=r(8)-dreal(ismt)*10.d0
      is=r(9)*100.d0+0.5d0
      iss1=is/1000
      iss2=(is-iss1*1000)/100
      iss3=(is-iss1*1000-iss2*100)/10
      iss4=is-iss1*1000-iss2*100-iss3*10
      hp=hu
      hs=hu
      if(r(12).gt.0.d0) hp=hl
      if(r(14).gt.0.d0) hs=hl
      iap=r(13)
      ias=r(15)
      tc1 = rc4
      tc2 = hp
      tc3 = rc1
      tc4 = hs
      smp=0.
      if(ichar(tc3).lt.43) then
        tc3=char(32)
      endif
      if((r(6).ge.0.d0).and.(r(10).ge.0.d0)) then
        smp=60.d0*(r(8)-r(4))+(r(9)-r(5))
      endif
      if(smp.lt.0.d0) then
        smp=smp+3600.d0
      endif
      if(ji.eq.1) then
        del_ns=r(2)
      endif
      if((r(6).ge.0.d0).and.(r(10).ge.0.d0)) then
        write(nsp,7)
     1  tc1,r(2),iazh,iazt,iazo,ipmt,ipmo,ips1,ips2,ips3,ips4,tc2,iap,
     2  tc3,r(7),r(6),r(16),r(18),smp,ismt,ismo,iss1,iss2,iss3,iss4,
     3  tc4,ias,r(11),r(10),r(17),r(19)
        ncparam=1
        call msa_wrt(ncparam,tc1,r(7),r(11))
      endif
    7 format(1x,a4,f7.2,1x,3i1,' p ',2i1,'-',2i1,'.',2i1,1x,a1,1x,
     1       i2,1x,a1,1x,f5.2,3(1x,f4.2),/,'     s-p =',f6.2,
     2       ' s ',2i1,'-',2i1,'.',2i1,1x,a1,1x,i2,3x,f5.2,3(1x,f4.2))
      if((r(6).ge.0.d0).and.(r(10).lt.0.d0)) then
        write(nsp,8)
     1  tc1,r(2),iazh,iazt,iazo,ipmt,ipmo,ips1,ips2,ips3,ips4,tc2,iap,
     2  tc3,r(7),r(6),r(16),r(18)
        r11tmp=0.
        ncparam=1
        call msa_wrt(ncparam,tc1,r(7),r11tmp)
      endif
    8 format(1x,a4,f7.2,1x,3i1,' p ',2i1,'-',2i1,'.',2i1,1x,a1,1x,
     1       i2,1x,a1,1x,f5.2,3(1x,f4.2))
      if((r(6).lt.0.d0).and.(r(10).ge.0.d0)) then
        write(nsp,9)
     1  tc1,r(2),iazh,iazt,iazo,ismt,ismo,iss1,iss2,iss3,iss4,tc4,ias,
     2  r(11),r(10),r(17),r(19)
        r7tmp=0.
        ncparam=1
        call msa_wrt(ncparam,tc1,r7tmp,r(11))
      endif
    9 format(1x,a4,f7.2,1x,3i1,' s ',2i1,'-',2i1,'.',2i1,1x,a1,1x,
     1       i2,3x,f5.2,3(1x,f4.2))
   10 continue
      if(nmstn.ge.1) then
        do 12 i=1,nmstn
        tc5 = stnm(i)
        write(nsp,11) tc5
   11   format(' no location on file for - ',a4)
   12   continue
        nmstn=0
      endif
c
      ncparam=2
      call msa_wrt(ncparam,tc1,r(7),r(11))
      if(imevt.gt.0) then
        ncparam=3
        call msa_wrt(ncparam,tc1,r(7),r(11))
        close(27)
      endif
c
      return
      end

      subroutine msa_wrt(ncparam,astn,presid,sresid)
c
c ...... write master event station adjustments to file
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*4 astn
      real*8 presid,sresid
      character*4 stn
      character*1 polar
      integer*4 nstns
      common/obs1/stn(1000),polar(1000)
      common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns
      character*4 tstn(1000),qstn(1000)
      real*8 tpres(1000),tsres(1000),qpres(1000),qsres(1000)
      common /eqelp/ del_ns,hor_se,iaz_gap,np_loc,ns_loc
      data cnall /0/
      if(ncall.eq.0) then
        ijk=0
        do 4 i=1,1000
        qstn(i)='    '
        qpres(i)=0.d0
        qsres(i)=0.d0
    4   continue
        ncall=1
      endif
c
      if(ncparam.eq.1) then
        if(astn(1:1).ne.' ') then
          ijk=ijk+1
          tstn(ijk)=astn
          tpres(ijk)=presid
          tsres(ijk)=sresid
        endif
        return
      elseif(ncparam.ge.2) then
        if(ncparam.eq.2) then
          ns_loc=0
        endif
c        write(*,'(a)') 'st_sln parameters -'
c        do 12 i=1,ijk
c        write(*,'(a,i3,3x,a4,3x,2f10.3)') 'i, stn, padj & sadj =',
c     1        i,tstn(i),tpres(i),tsres(i)
c   12   continue
        iii=1
        qstn(1)=tstn(1)
        qpres(1)=tpres(1)
        qsres(1)=tsres(1)
        do 1 i=2,ijk
        itstn=16777216*ichar(tstn(i)(1:1))+
     1           65536*ichar(tstn(i)(2:2))+
     2             256*ichar(tstn(i)(3:3))+
     3                 ichar(tstn(i)(4:4))
        iqstn=16777216*ichar(qstn(iii)(1:1))+
     1           65536*ichar(qstn(iii)(2:2))+
     2             256*ichar(qstn(iii)(3:3))+
     3                 ichar(qstn(iii)(4:4))
        if(itstn.eq.iqstn) then
          if(tpres(i).ne.0.) qpres(iii)=tpres(i)
          if(tsres(i).ne.0.) qsres(iii)=tsres(i)
        else
          iii=iii+1
          qstn(iii)=tstn(i)
          if(tpres(i).ne.0.) qpres(iii)=tpres(i)
          if(tsres(i).ne.0.) qsres(iii)=tsres(i)
        endif
    1   continue
      endif
c
      do 3 j=1,iii
      iqstn=16777216*ichar(qstn(j)(1:1))+
     1         65536*ichar(qstn(j)(2:2))+
     2           256*ichar(qstn(j)(3:3))+
     3               ichar(qstn(j)(4:4))
      do 2 i=1,nstns
      istn=16777216*ichar(stn(i)(1:1))+65536*ichar(stn(i)(2:2))+
     1          256*ichar(stn(i)(3:3))+      ichar(stn(i)(4:4))
c      write(*,'(1h ,a,a4,3x,i10)') 'stn & istn = ',stn,istn
      if(iqstn.eq.istn) then
        psadj=psaf(i)+qpres(j)
        ssadj=ssaf(i)+qsres(j)
        if(ncparam.eq.2) then
          ns_loc=ns_loc+1
        endif
        if(ncparam.eq.3) then
          write(27,'(1x,a4,2f6.2)') qstn(j),psadj,ssadj
        endif
        go to 25
      endif
    2 continue
   25 continue
    3 continue
c
      return
      end

      subroutine nsfig(q,n)
c
c ...... estimate number of significant figures remaining in solution
c
      implicit real*8 (a-h,o-z)
      common/prec/precis,cndn
      real*8 q(4)
c
c ...... sfibm = precision of double precision 
c                mantissa for a 64 bit IBM word
c
      sfibm=1.0d16
      spmp=dlog10(sfibm)
c      write(*,*) 'spmp = ',spmp
      qmax=q(1)
      qmin=q(1)
      do 1 i=1,n
      if(q(i).gt.qmax) qmax=q(i)
      if((q(i).gt.0.d0).and.(q(i).lt.qmin)) qmin=q(i)
    1 continue
      if(qmin.gt.0.d0) then
        cn=qmax/qmin
        if(cn.gt.99999.99) cn=99999.99
      else
        cn=99999.99d0
      endif
      sigfig=spmp-dlog10(cn)
      precis=1.d0/10**sigfig
      cndn=cn
c
      return
      end

      subroutine daz(x1,y1,x2,y2,d,az)
c
c ...... distance and azimuth calculation
c
c        input -
c           x1,x2 = longitude in radians (+ east)
c           y1,y2 = geocentric latitude in radians (+ north)
c
c        returns -
c           d  = distance in kilometers
c           az = azimuth in radians (CW from north)
c
c        Note - The distance and azimuth are calculated using a
c               tangent plane approximation.  This method is
c               accurate to within 0.1 km at 400 km and 0.3 km
c               at 700 km.
c
      implicit real*8 (a-h,o-z)
      data con/6371.0d0/
      pi=4.d0*datan(1.d0)
      yy=con*(y2-y1)
      xx=-con*(x2-x1)*dcos((y1+y2)/2.)
      arg=xx*xx+yy*yy
      if(arg.gt.0.d0) then
        nnn=13
        d=dsqrtf(nnn,arg)
      else
        d=0.
      endif
      if((xx.eq.0.d0).and.(yy.eq.0.d0)) then
        az=0.d0
      else
        az=datan2(xx,yy)
      endif
      if(az.lt.0.d0) az=az+2.d0*pi
c
      return
      end

      subroutine dazd(x1,y1,x2,y2,d,az)
c
c ...... distance and azimuth calculation
c
c        input -
c           x1,x2 = longitude in degrees (+ east)
c           y1,y2 = geocentric latitude in degrees (+ north)
c
c        returns -
c           d  = distance in kilometers
c           az = azimuth in degrees (CW from north)
c
c        Note - The distance and azimuth are calculated using a
c               tangent plane approximation.  This method is
c               accurate to within 0.1 km at 400 km and 0.3 km
c               at 700 km.
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 kpd
      data kpd /111.195d0/
      yy=kpd*(y2-y1)
      xx=kpd*(x2-x1)*dcosd((y1+y2)/2.)
      arg=xx*xx+yy*yy
      if(arg.gt.0.d0) then
        nnn=14
        d=dsqrtf(nnn,arg)
      else
        d=0.
      endif
      if((xx.eq.0.d0).and.(yy.eq.0.d0)) then
        az=0.d0
      else
        az=datan2d(xx,yy)
      endif
      if(az.lt.0.d0) az=az+360.d0
c
      return
      end

