anchor.diff:        RETURN
         END
   
  -C----------------------------------------------------------------------C
  -
  -      SUBROUTINE gcmdln(fltfile,stafile,haffile,dspfile,auto)
  +C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C
  +
  +      SUBROUTINE autogrid(stafile,nflt,evlo,evla,evdp,
  +     1                    str,dip,rak,dx,dy,area,slip,vp,vs,dens,flttyp)
  +C----
  +C Automatically determine grid of stations to calculate displacements
  +C----
  +      IMPLICIT none
  +      REAL*8 pi,r2d,d2r
  +      PARAMETER (pi=4.0d0*atan(1.0d0),r2d=180.0d0/pi,d2r=pi/180.0d0)
  +      INTEGER STER
  +      PARAMETER (STER=0)
  +      CHARACTER*40 stafile
  +      INTEGER nflt,FMAX,i,j,NTHR,dir
  +      PARAMETER (FMAX=1500,NTHR=3)
  +      INTEGER flttyp
  +      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  +     1       rak(FMAX),slip(FMAX),area(FMAX),dx(FMAX),dy(FMAX),
  +     2       dl(NTHR),thr(NTHR)
  +      REAL*8 stdp,lon,lat,avglon,avglat,az,Nmax,Smax,Emax,Wmax,
  +     1       vp,vs,dens,dN,dE,dZ,N,S,E,W,incr,distmax
  +
  +      write(STER,'(A)') 'Determining receiver grid limits...'
  +      do 110 i = 1,nflt
  +          avglon = avglon + evlo(i)
  +          avglat = avglat + evla(i)
  +  110 continue
  +      avglon = avglon/dble(nflt)
  +      avglat = avglat/dble(nflt)
  +C----
  +C Limit maximum extent of receiver grid
  +C NB: Nmax and Smax are returned in RADIANS!
  +C----
  +      distmax = 750.0d0 ! km
  +      az = 0.0d0
  +      call dlola(lon,Nmax,avglon,avglat,distmax,az)
  +      az = 180.0d0
  +      call dlola(lon,Smax,avglon,avglat,distmax,az)
  +      Emax = avglon+distmax*3.6d2/(2.0d0*pi*6.371d3*dcos(avglat*d2r))
  +      Wmax = avglon-distmax*3.6d2/(2.0d0*pi*6.371d3*dcos(avglat*d2r))
  +C----
  +C Set up distance steps and corresponding displacement thresholds
  +C----
  +      dl(1) = 1.0d0 ! Amounts to move test location
  +      dl(2) = 0.1d0
  +      dl(3) = 0.05d0
  +      thr(1) = 0.5d0 ! Displacement thresholds
  +      thr(2) = 0.01d0
  +      thr(3) = 0.001d0
  +C----
  +C Find limits to north, south, east, and west
  +C Variable 'dir' controls direction: 1=N, 2=S, 3=E, 4=W
  +C----
  +      stdp = 0.0d0
  +      do 103 dir = 1,4
  +          lon = avglon
  +          lat = avglat
  +C         Take first step (by increment of dl(1))
  +          if (dir.eq.1) lat = dnint((lat+dl(1))*1.d2)*1.d-2
  +          if (dir.eq.2) lat = dnint((lat-dl(1))*1.d2)*1.d-2
  +          if (dir.eq.3) lon = dnint((lon+dl(1))*1.d2)*1.d-2
  +          if (dir.eq.4) lon = dnint((lon-dl(1))*1.d2)*1.d-2
  +C         Take steps of decreasing size and threshold value
  +          i = 1
  +  101     if (i.gt.NTHR) goto 102
  +C             Check if lon, lat exceed maximum allowed values
  +              if (dir.eq.1.and.lat.gt.Nmax*r2d) then
  +                  lat = Nmax*r2d
  +                  goto 102
  +              elseif (dir.eq.2.and.lat.lt.Smax*r2d) then
  +                  lat = Smax*r2d
  +                  goto 102
  +              elseif (dir.eq.3.and.lon.gt.Emax) then
  +                  lat = Emax
  +                  goto 102
  +              elseif (dir.eq.4.and.lon.lt.Wmax) then
  +                  lon = Wmax
  +                  goto 102
  +              endif
  +C             Calculate displacement and compare to current threshold
  +              call calcdisp(dN,dE,dZ,lon,lat,stdp,nflt,evlo,evla,evdp,
  +     1                      str,dip,rak,dx,dy,area,slip,vp,vs,dens,
  +     2                      flttyp)
  +              dN = dabs(dN)
  +              dE = dabs(dE)
  +              dZ = dabs(dZ)
  +              if (dN.gt.thr(i).or.dE.gt.thr(i).or.dZ.gt.thr(i)) then
  +                  if (dir.eq.1) lat = lat + dl(i)
  +                  if (dir.eq.2) lat = lat - dl(i)
  +                  if (dir.eq.3) lon = lon + dl(i)
  +                  if (dir.eq.4) lon = lon - dl(i)
  +                  goto 101
  +              else
  +                  i = i + 1
  +                  goto 101
  +              endif
  +  102     continue
  +          if (dir.eq.1) N = lat
  +          if (dir.eq.2) S = lat
  +          if (dir.eq.3) E = lon
  +          if (dir.eq.4) W = lon
  +  103 continue
  +C----
  +C Create gridded stafile with separation between points = 'incr'
  +C----
  +      incr = 0.05d0
  +      open(unit=51,file=stafile,status='unknown')
  +      i = 0
  +      lon = W
  +  104 if (lon.gt.E) goto 107
  +          i = i + 1
  +          lat = S
  +          j = 0
  +  105     if (lat.gt.N) goto 106
  +              write(51,8001) lon,lat,stdp
  +              lat = lat + incr
  +              j = j + 1
  +              goto 105
  +  106     continue
  +          lon = lon + incr
  +          goto 104
  +  107 continue
  +      open(unit=52,file='autogrid.dat',status='unknown')
  +      write(STER,9001) W,E,S,N
  +      write(STER,9002) incr
  +      write(STER,9003) i*j
  +      write(52,9001) W,E,S,N
  +      write(52,9002) incr
  +      write(52,9003) i*j
  +      write(52,9004) i
  +      write(52,9005) j
  + 8001 format(3F10.2)
  + 9001 format('Grid limits (WESN):',4F8.2)
  + 9002 format('Grid spacing: ',F10.2)
  + 9003 format('Number of points: ',I10)
  + 9004 format('nlon: ',I10)
  + 9005 format('nlat: ',I10)
  +      close(51)
  +
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE mt2flt(fltfile,mtfile)
  +      IMPLICIT none
  +      CHARACTER*40 fltfile,mtfile
  +      REAL*8 lon,lat,dep,str,dip,rak,mag,wid,len,vp,vs,dens,mu,slip
  +      INTEGER typ
  +      call sethafdef(vp,vs,dens)
  +      mu = vs*vs*dens
  +      open(unit=61,file=mtfile,status='old')
  +      open(unit=62,file=fltfile,status='unknown')
  +  101 read(61,*,end=102) lon,lat,dep,str,dip,rak,mag
  +          if (45.0d0.lt.rak.and.rak.lt.135.0d0) then
  +              typ = 2
  +          elseif (-135.0d0.lt.rak.and.rak.lt.-45.0d0) then
  +              typ = 3
  +          else
  +              typ = 1
  +          endif
  +          call wellscoppersmith(wid,len,mag,typ)
  +          slip = ((1.d1**(1.5d0*(mag+10.7d0)))*1.d-7)/(mu*wid*len*1.d6)
  +          write(62,*)lon,lat,dep,str,dip,rak,slip,wid,len
  +          goto 101
  +  102 continue
  +      close(61)
  +      close(62)
  +      RETURN
  +      END
  +
  +C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C
  +
  +      SUBROUTINE wellscoppersmith(wid,len,mag,typ)
  +      IMPLICIT none
  +      INTEGER typ
  +      REAL*8 wid,len,mag
  +
  +      if (typ.eq.1) then
  +          len = 10.0d0**(-2.57d0+0.62d0*mag)
  +          wid = 10.0d0**(-0.76d0+0.27d0*mag)
  +      elseif (typ.eq.2) then
  +          len = 10.0d0**(-2.42d0+0.58d0*mag)
  +          wid = 10.0d0**(-1.61d0+0.41d0*mag)
  +      elseif (typ.eq.3) then
  +          len = 10.0d0**(-1.88d0+0.50d0*mag)
  +          wid = 10.0d0**(-1.14d0+0.35d0*mag)
  +      else
  +          len = 10.0d0**(-2.44d0+0.59d0*mag)
  +          wid = 10.0d0**(-1.01d0+0.32d0*mag)
  +      endif
  +
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE makegmtfile(gmtfile,nflt,evlo,evla,str,dip,dx,dy,slip)
  +      IMPLICIT none
  +      REAL*8 pi,d2r
  +      PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/180.0d0)
  +      CHARACTER*40 gmtfile
  +      INTEGER i,nflt,FMAX
  +      PARAMETER (FMAX=1500)
  +      REAL*8 evlo(FMAX),evla(FMAX),str(FMAX),dip(FMAX),
  +     1       slip(FMAX),dx(FMAX),dy(FMAX),cosd
  +      open(unit=71,file=gmtfile,status='unknown')
  +      do 301 i = 1,nflt
  +          cosd = dcos(dip(i)*d2r)
  +          write(71,1001) evlo(i),evla(i),slip(i),str(i),
  +     1                   dx(i)*1.0d-3,dy(i)*1.0d-3*cosd
  +  301 continue
  +      close(71)
  + 1001 format(2F10.3,F8.2,F8.1,2F10.3)
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE progbar(prog,prog100,progtag)
  +      IMPLICIT none
  +      INTEGER prog,prog100,progtag
  +      CHARACTER*1 CR
  +      CR = char(13) ! Carriage return
  +      if (100*prog/prog100.ge.progtag) then
  +          write (*,1000,advance='no') 100*prog/prog100,CR
  +          progtag = progtag + 1
  +      endif
  +      if (100*prog/prog100.ge.100) write (*,1000) 100*prog/prog100
  + 1000 format ('flt2displacement: [',I3,'% Complete]',A)
  +      RETURN
  +      END
  +
  +C----------------------------------------------------------------------C
  +
  +      SUBROUTINE gcmdln(fltfile,mtfile,stafile,haffile,dspfile,flttyp,
  +     1                  auto,mt,gmtfile,gmt)
         IMPLICIT none
         CHARACTER*40 tag
         INTEGER narg,i
  -      CHARACTER*40 fltfile,stafile,haffile,dspfile
  -      INTEGER auto
  +      CHARACTER*40 fltfile,mtfile,stafile,haffile,dspfile,gmtfile
  +      INTEGER auto,mt,flttyp,gmt
   
         fltfile = 'faults.txt'
         stafile = 'stations.txt'
         haffile = 'structure.txt'
         dspfile = 'disp.out'
  +      mtfile = 'none'
  +      gmtfile = 'none'
  +      mt = 0
         auto = 0
  +      gmt = 0
  +      flttyp = 1
   
         narg = iargc()
         i = 0
  
positive.code:       REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),
       1       dip(FMAX),rak(FMAX),slip(FMAX),dy(FMAX),dx(FMAX),
       2       area(FMAX)
        INTEGER flttyp(FMAX)
  
        open (unit=21,file=fltfile,status='old')
        f = 0
    201 f = f + 1
        read (21,*,END=202) evlo(f),evla(f),evdp(f),str(f),
       1                    dip(f),rak(f),slip(f),dy(f),dx(f),
       2                    flttyp(f)
            evdp(f) = 1.0d3*evdp(f)
            dx(f) = 1.0d3*dx(f)
            dy(f) = 1.0d3*dy(f)
            area(f) = dx(f)*dy(f)
            goto 201
    202 continue
  
negative.code:           dist = dist*6.371d6
            x = dist*( dcos(az-d2r*str(f)))
            y = dist*(-dsin(az-d2r*str(f)))
            if (flttyp.eq.0) then
                call o92pt(ux,uy,uz,x,y,stdp,evdp(f),dip(f),rak(f),
       1                   area(f),slip(f),vp,vs,dens)
            else
  