c.....................................................................
c
c       Extract mag data from nav.afx.ascii files
c
c.....................................................................
c

        character*72 file_nav,file_mag
        character*8  cf*80,arg*80,klat*12,klon*12,kar*20,blan*10
        real*4  rlat,rlon
        CHARACTER NSTA*10
c
        dimension ab(6,200),nsta(200),dist(200),az(200),azin(200)
        double precision az,dist
        common/dados/ab,nd,nsta,dist,az,azin,vp,vs
        character cdata*10,ctime*8
c
        namelist /par/ file_nav,file_mag,n_aver
c.......................................................................
        vp=6.0
        vs=3.46
        n_aver=1
        blan='          '
c
        open(unit=10, file='nav03m.in', status='old')
        read(10,par)
        close(10)
c..................................................navigation file *.nav
        open(9,file=file_nav,status='old')
c.........................................................mag file *.mag
        open(10,file=file_mag,status='unknown')
c
        ipon=0
        do 1 i=1,1000000
          read(9,'(a)',err=1,end=2) cf
          do 777 kk=1,80
            if (ichar(cf(kk:kk)).eq.9) cf(kk:kk)=' '
777       continue
          numarg=iargc2(cf)
          if (numarg.le.2) goto 1
          ipon=ipon+1
c
          call getar2(1,arg,cf)
          kar=arg
          ll=lenchr(kar)
          llb=12-ll
          klon=blan(1:llb)//kar(1:ll)
          read (klon,'(g14.7)') rlon
c
          call getar2(2,arg,cf)
          kar=arg
          ll=lenchr(kar) 
          llb=11-ll
          klat=blan(1:llb)//kar(1:ll)
          read (klat,'(g14.7)') rlat
c
          if (ipon.eq.1) then
            rlat0=rlat
            rlon0=rlon
            ab(1,1)=rlat0
            ab(2,1)=0.0
            ab(3,1)=0.0
            ab(4,1)=rlon0
            ab(5,1)=0.0
            ab(6,1)=0.0
          endif
          ab(1,2)=rlat
          ab(2,2)=0.0
          ab(3,2)=0.0
          ab(4,2)=rlon
          ab(5,2)=0.0
          ab(6,2)=0.0
          nd=2
          call subdis()
          dddd=sngl(dist(2))
c
          call getar2(5,arg,cf)
          cdata=arg(1:10)
          call getar2(6,arg,cf)
          ctime=arg(1:8)
          call getar2(3,arg,cf)
          read(arg,'(g14.7)',err=1) depth
          call getar2(4,arg,cf)
          read(arg,'(g14.7)',err=1) rmag
          if (rmag.gt.10.) write(10,2000) cdata,ctime,
     +    klon,klat,rmag,dddd
2000      format(a10,1x,a8,1x,a12,1x,a11,1x,f7.1,1x,f7.3)
1       continue
2       npon=ipon
        write(*,*) npon,'points found'
        close(9)
c
999     stop
        end
c
c***********************************************************************
c
      subroutine getar2(iarg,arg,cf)
c-----------------------------------------------------------------------
c
c   Emulates the UNIX-FORTRAN subroutine "getarg".
c
c   Input via call:
c
c   iarg - index of desired command line argument
c
c
c   Output via call:
c
c   arg - a character*80 variable with the command line argument
c
c   Note: This version reads and parses the command line each time
c         the routine is called.  This is not very satisfactory,
c         however I am at a loss as to how to do this and still
c         preserve the UNIX arguments.
c
c   Requires you link with /lib/st.lib
c
c-----------------------------------------------------------------------
      implicit integer*2 (z)
      character*1 comlin(127)
      character*80 comlst(50), arg
      character*127 comand
      character*80 cf
c
      if (iarg .le. 0) stop 'GETARG Error: Illegal iarg'
      if (iarg .eq. 0) stop 'GETARG Error: no command line parameters'
      comand = ' '
      call parse2(cf,comlst,numarg)
      if (numarg .lt. iarg) stop 'GETARG Error: iarg too large'
      arg = comlst(iarg)
c
      return
      end
      function iargc2(cf)
c----------------------------------------------------------------------
c
c   Emulates the UNIX-FORTRAN fuction "iargc"
c
c   Returns the number of arguments on the command line.
c
c
c----------------------------------------------------------------------
      implicit integer*2 (z)
      character comand*127, comlst*80, comlin*1
      dimension comlst(50), comlin(127)
      character*80 cf
c
        call parse2(cf,comlst,num)
        iargc2 = num
c
      return
      end
      subroutine parse2(comand,comlst,numarg)
c----------------------------------------------------------------------
c
c   Parses command line.
c
c   Input via call:
c   comand -=- command string (must be character*127)
c
c   Output via call
c   comlst -=- array containing each blank delimited element
c              of "comand" (must be dimensioned as below)
c   numarg -=- number of elements in comlst
c
c----------------------------------------------------------------------
      implicit integer*2 (z)
      character comand*80, comlst*80
      dimension comlst(50)
c
c   Remove leading blanks (return if all blank), count trailing blanks
c
      if (comand .eq. ' ') then
        numarg = 0
        return
      endif
      zblank = 0
      do 10 zi = 1, 80
        if (comand(zi:zi) .ne. ' ') goto 20
        zblank = zblank + 1
10    continue
20    if (zblank .ne. 0) then
        comand(1:80-zblank)=comand(zblank+1:80)
        comand(80-zblank+1:80)=' '
      endif
      zblank = 0
      zindex = 81
      do 30 zi = 1, 80
        zindex = zindex - 1
        if (comand(zindex:zindex) .ne. ' ') goto 40
        zblank = zblank + 1
30    continue
c
c   Parse comand into individual pieces...
c
40              zend = 80 - zblank
                numarg = 0
                zsave = 1
c
c          First, look for end of first word
c
                do 70 zi = 1, zend + 1
                if (comand(zi:zi) .eq. ' ' .and. zi .gt. zsave) then
c
c          Have reached end of word, so save it
c
                  numarg = numarg + 1
                  comlst(numarg) = comand(zsave:zi-1)
c
c          Now, look for beginning of next word
c
                  do 50 zsave = zi+1, zend
                    if (comand(zsave:zsave) .ne. ' ') goto 60
50                continue
60                continue
                  endif
c
c          Go back up and look for end of this word
c
70              continue
      return
c
      end
c
c************************************************************
c
      subroutine subdis()
c
      INTEGER RR,SS,PD,RRR,S,TD,ALFA,CON
      CHARACTER NSTA*10,IFILEN*9,NEPI*9,NSTA1*10,idtiro*4
      CHARACTER*1 ICAR1,fout*20
      dimension ab(6,200),nsta(200),dist(200),az(200),azin(200)
      DIMENSION BL(200),BP(200),ALTI(200),ICW(200)
      DIMENSION TH(200),PHI(200),XDEG(200)
      DIMENSION DELTA(200),DELTE(200),T4(200),T5(200)
      DOUBLE PRECISION THK,PHK,SINTK,COSTK,TANTK,TH2,V1,THG,D,E,F,A,
     %B,C,THC,PHC,TANTC,TH2I,V2,AL,DL,A12,AZ,COSA12,DM,EO,EO2,EO3,
     %EOP,EOQ,U1,BO,CO,C2,C4,Z1,Z2,X2,Y2,U2,DIST,D1,E1,F1,C1,A1,B1,
     %SC,SD,XDEG,TEST,DELTA,T4,T5,T2,DELTE,RADI,RAD,EC,ELL,EL,EP,
     %BL,EEL,EEP,BP,TH,PHI,PAXIS,EAXIS
      common/dados/ab,nd,nsta,dist,az,azin,vp,vs
c
      DO 113 I=1,ND
        BP(I)=AB(1,I)+AB(2,I)/60.+AB(3,I)/3600.
        ICW(I)=1
        if(AB(6,I).LT.0.) ICW(I)=-1
        if(AB(6,I).LT.0.) AB(6,I)=-AB(6,I)
        if(AB(5,I).LT.0.) ICW(I)=-1
        IF(AB(5,I).LT.0.) AB(5,I)=-AB(5,I)
        IF(AB(4,I).LT.0) ICW(I)=-1
        IF(ICW(I).LT.0) GO TO 232
        BL(I)=AB(4,I)+AB(5,I)/60.+AB(6,I)/3600.
        GO TO 113
 232    BL(I)=AB(4,I)-AB(5,I)/60.-AB(6,I)/3600.
 113  CONTINUE
C
 107  DO 108 I=1,ND
        TH(I)=BP(I)
        PHI(I)=BL(I)
 108  CONTINUE
c
      EEL=BL(1)
      EEP=BP(1)
      NEPI=NSTA(1)
c
      RADI=57.2957795
      RAD=1.0/RADI
      ELL=.99327733
      EC=.672267002E-02
C     ELL=.993277
C     EC=.672267E-02
      EL=EC/ELL
      EP=1.0+EL
      PAXIS=6356.912
      EAXIS=6378.388
      THK=RAD*EEP
      PHK=RAD*EEL
      SINTK= DSIN(THK)
      COSTK= DCOS(THK)
      TANTK=SINTK/COSTK
      TH2=THK+THK
      V1=DEXP(.230258509E+01*(.380544268E+01-(.7323684E-3)*DCOS(TH2)
     1+(.6175 E-6)*DCOS(TH2+TH2)-(.7E-9)*DCOS(TH2+TH2+TH2)))
C     V1=EXP(.230259E+01*(.380544E+01-(.732368E-3)*COS(TH2)
C    1+(.6175 E-6)*COS(TH2+TH2)-(.7E-9)*COS(TH2+TH2+TH2)))
      THG=DATAN(ELL*TANTK)
      D=DSIN(PHK)
      E=-DCOS(PHK)
      F=-DCOS(THG)
      A=F*E
      B=-D*F
      C=DSIN(THG)
      DO 38 I=2,ND
        THC=RAD*TH(I)
        PHC=RAD*PHI(I)
        TANTC=DSIN(THC)/DCOS(THC)
        TH2I= THC+THC
       V2=DEXP(.230258509E+01*(.380544268E+01-(.7323684E-3)*DCOS(TH2I)
     1+(.6175 E-6)*DCOS(TH2I+TH2I)-(.7E-9)*DCOS(TH2I+TH2I+TH2I)))
C     V2=EXP(.230259E+01*(.380544E+01-(.732368E-3)*COS(TH2I)
C    1+(.6175 E-6)*COS(TH2I+TH2I)-(.7E-9)*COS(TH2I+TH2I+TH2I)))
        AL= TANTC/(EP*TANTK)+EC*DSQRT((EP+TANTC*TANTC)/(EP+TANTK*TANTK))
        DL= PHC-PHK
        A12=DATAN2(DSIN(DL),((AL-DCOS(DL))*SINTK))
        AZ(I)=RADI*A12
        IF(AZ(I).LT.0.0) AZ(I)=AZ(I)+360.
        COSA12=DCOS(A12)
        DM=COSTK*COSA12
        DM=DM*DM+SINTK*SINTK
        EO= EL*DM
        EO2=EO*EO
        EO3=EO*EO2
        EOP=1.+EO
        EOQ= DSQRT(EOP)
        U1=DATAN2(SINTK,(COSTK*EOQ*COSA12))
        BO=V1*DSQRT(1.+EL*COSTK*COSTK*COSA12*COSA12)/EOP
        CO=1.+0.25*EO-(3.*EO2-1.25*EO3)/64.
        C2=-0.125*EO+EO2/32.-15.0*EO3/1024.
        C4=(-EO2+0.75*EO3)/256.
        Z1=V1*(1.-EC)*SINTK
        Z2=V2*(1.-EC)*DSIN(THC)
        X2=V2*DCOS(THC)*DCOS(DL)
        Y2=V2*DCOS(THC)*DSIN(DL)
        U2=DATAN((V1*SINTK+EOP*(Z2-Z1))/(EOQ*(X2*COSA12-Y2*SINTK*
     &     DSIN(A12))))
        DIST(I)=BO*(CO*(U2-U1)+C2*(DSIN(U2+U2)-DSIN(U1+U1))
     &          +C4*(DSIN(U2+U2+U2+U2)-DSIN(U1+U1+U1+U1)))
        DIST(I)=DABS(DIST(I))
        THG= DATAN(ELL*TANTC)
        D1=DSIN(PHC)
        D1=DSIN(PHC)
        E1=-DCOS(PHC)
        F1=-DCOS(THG)
        C1=DSIN(THG)
        A1=F1*E1
        B1=-D1*F1
        SC=A*A1+B*B1+C*C1
        SD=DSQRT(0.25*((A-A1)*(A-A1)+(B-B1)*(B-B1)+(C-C1)*(C-C1))*
     &     ((A+A1)*(A+A1)+(B+B1)*(B+B1)+(C+C1)*(C+C1)))
        XDEG(I)=RADI*DATAN2(SD,SC)
        TEST=DIST(I)-111.0*XDEG(I)
        IF(DABS(TEST)-100.)24,23,23
  23    U2=U2+0.314159265E+01
        DIST(I)=BO*(CO*(U2-U1)+C2*(DSIN(U2+U2)-DSIN(U1+U1))
     &          +C4*(DSIN(U2+U2+U2+U2)-DSIN(U1+U1+U1+U1)))
  24    DIST(I)=DABS(DIST(I))
        DELTA(I)=DIST(I)/vp
        T4(I)=T3+DELTA(I)
        T5(I)=T2
 202    IF(T4(I).GT.60.)T5(I)=T5(I)+1.
        IF(T4(I).GT.60.)T4(I)=T4(I)-60.
        IF(T4(I).GT.60.)GO TO 202
        DELTE(I)=DIST(I)/vs
  25  AL=TANTK/(EP*TANTC)+EC*DSQRT((EP+TANTK*TANTK)/(EP+TANTC*TANTC))
        AZIN(I)=RADI*DATAN2(DSIN(-DL),((AL-DCOS(DL))*DSIN(THC)))
        IF(AZIN(I).LT.0.0)AZIN(I)=AZIN(I)+360.0
  38  CONTINUE
c
c
         return
         end
c 
c******************************************************
c
        function lenchr(fn)
        character *(*) fn
        do 200 lenchr=len(fn),1,-1
          if (fn(lenchr:lenchr).ne.' ') return
 200    continue
        lenchr=0
        return
        end
c