anchor.diff:  C Parse the command line, check for FFM file, and read FFM
   C----
         call gcmdln(fltfile,stafile,haffile,dspfile,auto)
  +      inquire(file=fltfile,exist=ex)
  +      if (.not.ex) call usage('!! Error: no file named '//fltfile)
  +      call readfaults(fltfile,nflt,evlo,evla,evdp,str,dip,rak,slip,dy,
  +     1                dx,area,flttyp)
  +
  +C----
  +C If auto=0: use existing input files
  +C If auto=1: generate half-space and station inputs
  +C----
  +      if (auto.eq.0) then
  +          inquire(file=haffile,exist=ex)
  +          if (.not.ex) then
  +              write(STER,*) '!! Warning: no file named '//haffile
  +              write(STER,*) '!! Using default half-space parameters'
  +              call sethafdef(vp,vs,dens)
  +          else
  +              open (unit=23,file=haffile,status='old')
  +              read (23,*) vp,vs,dens
  +          endif
  +          inquire(file=stafile,exist=ex)
  +          if (.not.ex) call usage('!! Error: no file named '//stafile)
  +          open(unit=22,file=stafile,status='old')
  +      elseif (auto.eq.1) then
  +          call sethafdef(vp,vs,dens)
  +C          call autogrid(stafile,nflt,hypolon,hypolat,evlo,evla,evdp,str,
  +C     1                  dip,rak,dx,dy,area,slip,vp,vs,dens,flttyp)
  +          open (unit=22,file=stafile,status='old')
  +      endif
   
   C----
   C Output file
   C----
  -      open (unit=11,file='disp.out',status='unknown')
  -C----
  -C Parse fault file
  -C List of:
  -C  lon lat dep(km) str dip rak slip(m) dy(km) dx(km) flttyp
  -C Last column indicates fault type:
  -C   0 for point source
  -C   1 for finite source
  -C----
  -      call readfaults(nflt,evlo,evla,evdp,str,dip,rak,slip,dy,dx,area,
  -     1                flttyp,FMAX)
  +      open (unit=11,file=dspfile,status='unknown')
  +
   C----
   C Calculate net displacements at each station
   C----
  
positive.code: C----------------------------------------------------------------------C
  
        SUBROUTINE xy2NE(uN,uE,ux,uy,str)
  
        IMPLICIT none
        REAL*8 uN,uE,ux,uy,str,theta,uhor
        REAL*8 pi,d2r
        PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
  
        theta = datan2(uy,ux)
        uhor = dsqrt(ux*ux+uy*uy)
        theta = d2r*str - theta
        uN = uhor*dcos(theta)
        uE = uhor*dsin(theta)
  
        RETURN
        END
  
  
negative.code: C an isotropic elastic half-space due to shear dislocations.
  C----
        IMPLICIT NONE
        INTEGER STER
        PARAMETER (STER=0)
        CHARACTER*40 fltfile,stafile,haffile,dspfile
        LOGICAL ex
        INTEGER auto
  