anchor.diff:       1       '=========================')
    9999 format('-------------------------------------------------------',
        1       '-------------------------')
  -c      write(*,*)
  -c     1 '123456789012345678901234567890123456789012345678901234567890',
  -c     2 '12345678901234567890'
         STOP
         END
   
   
  -      SUBROUTINE dbggcmdln(ffmf,fltf,magf,haff,staf,trgf,dspf,stnf,stsf,
  -     1                     norf,shrf,coulf,gmtf,flttyp,auto,incr,val,
  -     2                     xy,prog,long)
  -      IMPLICIT none
  -      CHARACTER*40 ffmf,fltf,magf,haff,staf,trgf,
  -     1             dspf,stnf,stsf,norf,shrf,coulf,
  -     2             gmtf
  -      INTEGER flttyp,auto,xy,prog,long
  -      REAL*8 incr,val
  -      INTEGER kffm,kflt,kmag,ksta,khaf,ktrg
  -      INTEGER kdsp,kstn,ksts,knor,kshr,kcou,kgmt
  -      COMMON /ICHECK/ kffm,kflt,kmag,ksta,khaf,ktrg
  -      COMMON /OCHECK/ kdsp,kstn,ksts,knor,kshr,kcou,kgmt
  -      print *,'DEBUGGING SUBROUTINE GCMDLN'
  -      print *,'---------------------------'
  -      print *,'KFFM',kffm,ffmf
  -      print *,'KFLT',kflt,fltf
  -      print *,'KMAG',kmag,magf
  -      print *,'KSTA',ksta,staf
  -      print *,'KHAF',khaf,haff
  -      print *,'KTRG',ktrg,trgf
  -      print *,'KDSP',kdsp,dspf
  -      print *,'KSTN',kstn,stnf
  -      print *,'KSTS',ksts,stsf
  -      print *,'KNOR',knor,norf
  -      print *,'KSHR',kshr,shrf
  -      print *,'KCOU',kcou,coulf
  -      print *,'KGMT',kgmt,gmtf
  -      print *,'FLTTYP',flttyp
  -      print *,'AUTO',auto
  -      print *,'INCR',incr
  -      print *,'VAL',val
  -      print *,'XY?',xy
  -      print *,'PROG',prog
  -      print *,'LONG',long
  -      RETURN
  -      END
  -
  -C||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||C
  -CVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVC
  -C THESE ARE ALL THE OLD SUBROUTINES FOR AUTOMATIC OPTIONS
  -C
  -C      SUBROUTINE autodispgrid(staf,evlo,evla,evdp,str,dip,rak,dx,dy,
  -C     1                        slip,hylo,hyla,nflt,incr,val,dspf,xy)
  -CC----
  -CC Automatically determine receiver grid for displacement computation,
  -CC write grid to STAFILE.
  -CC----
  -C      IMPLICIT NONE
  -C      CHARACTER*40 staf,dspf
  -C      INTEGER nflt,FMAX,xy
  -C      PARAMETER (FMAX=1500)
  -C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  -C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  -C      REAL*8 W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,incr,val
  -C      write(*,'(A)') ' Creating receiver grid in file: '//staf
  -C      dmax = 750.0d0 ! km
  -C      call getmaxlims(WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,xy)
  -C      call getautodisplims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  -C     1                     evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  -C     2                     nflt,xy,val)
  -C      call writestaf(staf,W,E,S,N,incr,val)
  -C      write(*,'(A)') ' Displacements written to file: '//dspf
  -C      RETURN
  -C      END
  -C
  -C
  -C      SUBROUTINE autocoulgrid(staf,evlo,evla,evdp,str,dip,rak,dx,dy,
  -C     1                        slip,hylo,hyla,nflt,incr,val,coulf,xy)
  -CC----
  -CC Automatically determine receiver grid for Coulomb stress computation,
  -CC write grid to COUFILE.
  -CC----
  -C      IMPLICIT NONE
  -C      CHARACTER*40 staf,coulf
  -C      INTEGER nflt,FMAX,xy
  -C      PARAMETER (FMAX=1500)
  -C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  -C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  -C      REAL*8 W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,incr,val
  -C      write(*,'(A)') ' Creating receiver grid in file: '//staf
  -C      dmax = 750.0d0 ! km
  -C      call getmaxlims(WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,xy)
  -C      call getautostnlims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  -C     1                    evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  -C     2                    nflt,xy,val)
  -C      call writestaf(staf,W,E,S,N,incr,val)
  -C      write(*,'(A)') ' Coulomb stress written to file: '//coulf
  -C      RETURN
  -C      END
  -C
  -C
  -C      SUBROUTINE getmaxlims(WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,xy)
  -C      IMPLICIT NONE
  -C      REAL*8 pi,r2d,d2r
  -C      PARAMETER (pi=4.0d0*atan(1.0d0),r2d=180.0d0/pi,d2r=pi/180.0d0)
  -C      REAL*8 az,dmax,WMAX,EMAX,SMAX,NMAX,hylo,hyla,junk
  -C      INTEGER xy
  -CC----
  -CC Determine maximum geographic limits of receiver grid to W,E,S,N.
  -CC----
  -C      if (xy.eq.0) then
  -C          WMAX = hylo-dmax*3.6d2/(2.0d0*pi*6.371d3*dcos(hyla*d2r))
  -C          EMAX = hylo+dmax*3.6d2/(2.0d0*pi*6.371d3*dcos(hyla*d2r))
  -C          az = 180.0d0
  -C          call dlola(junk,SMAX,hylo,hyla,dmax,az)
  -C          SMAX = SMAX*r2d
  -C          az = 0.0d0
  -C          call dlola(junk,NMAX,hylo,hyla,dmax,az)
  -C          NMAX = NMAX*r2d
  -C      else
  -C          WMAX = hylo - dmax
  -C          EMAX = hylo + dmax
  -C          SMAX = hyla - dmax
  -C          NMAX = hyla + dmax
  -C      endif
  -C      RETURN
  -C      END
  -C
  -C
  -C      SUBROUTINE getautodisplims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  -C     1                           evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  -C     2                           nflt,xy,val)
  -CC----
  -CC Get limits of receiver grid by taking steps along longitude and
  -CC latitudes until displacements smaller than minimum threshold value.
  -CC Variable 'dir' determines direction: 1=N, 2=S, 3=E, 4=W
  -CC----
  -C      IMPLICIT NONE
  -C      INTEGER i,NTHR
  -C      PARAMETER (NTHR=3)
  -C      REAL*8 dl(NTHR),thr(NTHR)
  -C      REAL*8 hylo,hyla,WMAX,EMAX,SMAX,NMAX,W,E,S,N,lo,la
  -C      REAL*8 uE,uN,uZ,stdp,vp,vs,dens,val
  -C      INTEGER nflt,FMAX,xy
  -C      PARAMETER (FMAX=1500)
  -C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  -C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  -C      INTEGER dir,flttyp
  -C      if (xy.eq.0) then
  -C          dl(1) = 1.0d0 ! Amounts to move test location (deg)
  -C          dl(2) = 0.1d0
  -C          dl(3) = 0.05d0
  -C      else
  -C          dl(1) = 100.0d0 ! Amounts to move test location (km)
  -C          dl(2) = 10.0d0
  -C          dl(3) = 5.0d0
  -C      endif
  -C      thr(1) = 0.5d0 ! Displacement thresholds (m)
  -C      thr(2) = 0.01d0
  -C      thr(3) = 0.001d0
  -C      flttyp = 1
  -C      call autohaf(vp,vs,dens)
  -C      stdp = val
  -C      do 903 dir = 1,4
  -C          lo = hylo
  -C          la = hyla
  -C         Take first step (by increment of dl(1))
  -C          if (dir.eq.1) la = dnint((la+dl(1))*1.0d2)*1.0d-2
  -C          if (dir.eq.2) la = dnint((la-dl(1))*1.0d2)*1.0d-2
  -C          if (dir.eq.3) lo = dnint((lo+dl(1))*1.0d2)*1.0d-2
  -C          if (dir.eq.4) lo = dnint((lo-dl(1))*1.0d2)*1.0d-2
  -C         Take steps of decreasing size and threshold value
  -C          i = 1
  -C  901     if (i.gt.NTHR) goto 902
  -C             Check if lon, lat exceed maximum allowed values
  -C              if (dir.eq.1.and.la.gt.NMAX) then
  -C                  la = NMAX
  -C                  goto 902
  -C              elseif (dir.eq.2.and.la.lt.SMAX) then
  -C                  la = SMAX
  -C                  goto 902
  -C              elseif (dir.eq.3.and.lo.gt.EMAX) then
  -C                  lo = EMAX
  -C                  goto 902
  -C              elseif (dir.eq.4.and.lo.lt.WMAX) then
  -C                  lo = WMAX
  -C                  goto 902
  -C              endif
  -C             Calculate displacement and compare to current threshold
  -C              call calcdisp(uN,uE,uZ,lo,la,stdp,nflt,evlo,evla,evdp,
  -C     1                      str,dip,rak,dx,dy,slip,vp,vs,dens,
  -C     2                      flttyp,xy)
  -C              uN = dabs(uN)
  -C              uE = dabs(uE)
  -C              uZ = dabs(uZ)
  -C              if (uN.gt.thr(i).or.uE.gt.thr(i).or.uZ.gt.thr(i)) then
  -C                  if (dir.eq.1) la = la + dl(i)
  -C                  if (dir.eq.2) la = la - dl(i)
  -C                  if (dir.eq.3) lo = lo + dl(i)
  -C                  if (dir.eq.4) lo = lo - dl(i)
  -C                  goto 901
  -C              else
  -C                  i = i + 1
  -C                  goto 901
  -C              endif
  -C  902     continue
  -C          if (dir.eq.1) N = la
  -C          if (dir.eq.2) S = la
  -C          if (dir.eq.3) E = lo
  -C          if (dir.eq.4) W = lo
  -C  903 continue
  -C      RETURN
  -C      END
  -C
  -C
  -C      SUBROUTINE getautostnlims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  -C     1                          evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  -C     2                          nflt,xy,val)
  -CC----
  -CC Get limits of receiver grid by taking steps along longitude and
  -CC latitudes until strains smaller than minimum threshold value.
  -CC Variable 'dir' determines direction: 1=N, 2=S, 3=E, 4=W
  -CC----
  -C      IMPLICIT NONE
  -C      INTEGER i,NTHR
  -C      PARAMETER (NTHR=3)
  -C      REAL*8 dl(NTHR),thr(NTHR)
  -C      REAL*8 hylo,hyla,WMAX,EMAX,SMAX,NMAX,W,E,S,N,lo,la
  -C      REAL*8 stn(3,3),stdp,vp,vs,dens,stnmax,stnmin,val
  -C      INTEGER nflt,FMAX
  -C      PARAMETER (FMAX=1500)
  -C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  -C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  -C      INTEGER dir,flttyp,xy
  -C      if (xy.eq.0) then
  -C          dl(1) = 1.0d0 ! Amounts to move test location (deg)
  -C          dl(2) = 0.1d0
  -C          dl(3) = 0.05d0
  -C      else
  -C          dl(1) = 100.0d0 ! Amounts to move test location (km)
  -C          dl(2) = 10.0d0
  -C          dl(3) = 5.0d0
  -C      endif
  -C      thr(1) = 1.0d-7 ! Strain thresholds (m)
  -C      thr(2) = 5.0d-8
  -C      thr(3) = 1.0d-8
  -C      flttyp = 1
  -C      call autohaf(vp,vs,dens)
  -C      stdp = val
  -C      do 913 dir = 1,4
  -C          lo = hylo
  -C          la = hyla
  -C         Take first step (by increment of dl(1))
  -C          if (dir.eq.1) la = dnint((la+dl(1))*1.0d2)*1.0d-2
  -C          if (dir.eq.2) la = dnint((la-dl(1))*1.0d2)*1.0d-2
  -C          if (dir.eq.3) lo = dnint((lo+dl(1))*1.0d2)*1.0d-2
  -C          if (dir.eq.4) lo = dnint((lo-dl(1))*1.0d2)*1.0d-2
  -C         Take steps of decreasing size and threshold value
  -C          i = 1
  -C  911     if (i.gt.NTHR) goto 912
  -C             Check if lon, lat exceed maximum allowed values
  -C              if (dir.eq.1.and.la.gt.NMAX) then
  -C                  la = NMAX
  -C                  goto 912
  -C              elseif (dir.eq.2.and.la.lt.SMAX) then
  -C                  la = SMAX
  -C                  goto 912
  -C              elseif (dir.eq.3.and.lo.gt.EMAX) then
  -C                  lo = EMAX
  -C                  goto 912
  -C              elseif (dir.eq.4.and.lo.lt.WMAX) then
  -C                  lo = WMAX
  -C                  goto 912
  -C              endif
  -C             Calculate strain and compare to current threshold
  -C              call calcstn(stn,lo,la,stdp,nflt,evlo,evla,evdp,
  -C     1                     str,dip,rak,dx,dy,slip,vp,vs,dens,
  -C     2                     flttyp,xy)
  -C              stnmax = dmax1(stn(1,1),stn(2,2),stn(3,3),stn(1,2),
  -C     1                       stn(1,3),stn(2,3))
  -C              stnmin = dmin1(stn(1,1),stn(2,2),stn(3,3),stn(1,2),
  -C     1                       stn(1,3),stn(2,3))
  -C              stnmax = dmax1(dabs(stnmax),dabs(stnmin))
  -C              if (stnmax.gt.thr(i)) then
  -C                  if (dir.eq.1) la = la + dl(i)
  -C                  if (dir.eq.2) la = la - dl(i)
  -C                  if (dir.eq.3) lo = lo + dl(i)
  -C                  if (dir.eq.4) lo = lo - dl(i)
  -C                  goto 911
  -C              else
  -C                  i = i + 1
  -C                  goto 911
  -C              endif
  -C  912     continue
  -C          if (dir.eq.1) N = la
  -C          if (dir.eq.2) S = la
  -C          if (dir.eq.3) E = lo
  -C          if (dir.eq.4) W = lo
  -C  913 continue
  -C      RETURN
  -C     END
  -C
  -C
  -C      SUBROUTINE writestaf(staf,W,E,S,N,incr,val)
  -CC----
  -CC Write receiver grid to STAFILE.
  -CC----
  -C      IMPLICIT NONE
  -C      CHARACTER*40 staf
  -C      REAL*8 W,E,S,N,incr,lat,lon,stdp,val
  -C      INTEGER i,j
  -C      stdp = val
  -C      open(unit=95,file=staf,status='unknown')
  -C      i = 0
  -C      lon = W
  -C  951 if (lon.gt.E) goto 954
  -C          i = i + 1
  -C          lat = S
  -C          j = 0
  -C  952     if (lat.gt.N) goto 953
  -C              write(95,8001) lon,lat,stdp
  -C              lat = lat + incr
  -C              j = j + 1
  -C              goto 952
  -C  953     continue
  -C          lon = lon + incr
  -C          goto 951
  -C  954 continue
  -C      open(unit=96,file='autogrid.dat',status='unknown')
  -C      write(*,9001) W,E,S,N
  -C      write(*,9002) incr
  -C      write(*,9003) i*j
  -C      write(*,9004) i
  -C      write(*,9005) j
  -C      write(96,9001) W,E,S,N
  -C      write(96,9002) incr
  -C      write(96,9003) i*j
  -C      write(96,9004) i
  -C      write(96,9005) j
  -C      write(*,*) 'This information is saved in autogrid.dat'
  -C 8001 format(3F10.2)
  -C 9001 format(' Grid limits (WESN):',4F8.2)
  -C 9002 format(' Grid spacing: ',F10.2)
  -C 9003 format(' Number of points: ',I10)
  -C 9004 format(' nlon: ',I10)
  -C 9005 format(' nlat: ',I10)
  -C      close(95)
  -C      close(96)
  -C      RETURN
  -C      END
  -
  +
  
positive.code:      1 '    Format: evlo evla evdp(km) str dip rak slip(m)  wid(km)   ',
       2                                                       'len(km)'
            write(*,*)
       1 '                                                   along-dip ',
       2                                                    'along-str'
            write(*,*)
            write(*,9999)
  
negative.code: C----------------------------------------------------------------------C
  
        SUBROUTINE gcmdln(ffmf,fltf,magf,haff,staf,trgf,dspf,stnf,stsf,
       1                  norf,shrf,coulf,gmtf,flttyp,auto,incr,adist,xy,
       2                  prog,long)
        IMPLICIT NONE
        CHARACTER*40 tag
  