anchor.diff:  C an isotropic elastic half-space due to shear dislocations.
   C----
         IMPLICIT NONE
  -      LOGICAL verbose
  -      CHARACTER*20 srcfile,stafile,haffile
  -      INTEGER nflt,maxflt
  -      PARAMETER (maxflt=1500)
  -      INTEGER flttyp(maxflt)
  -      REAL*8 evlo(maxflt),evla(maxflt),evdp(maxflt),str(maxflt),
  -     1       dip(maxflt),rak(maxflt),slip(maxflt),dx(maxflt),dy(maxflt),
  -     2       area(maxflt)
  -
  +      CHARACTER*40 fltfile,stafile,haffile,dspfile
  +      LOGICAL ex
  +      INTEGER auto
  +      INTEGER nflt,FMAX
  +      PARAMETER (FMAX=1500)
  +      INTEGER flttyp(FMAX)
  +      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  +     1       rak(FMAX),slip(FMAX),dx(FMAX),dy(FMAX),area(FMAX)
         REAL*8 stlo,stla,stdp
         REAL*8 vp,vs,dens
  -
  +      INTEGER prog,prog100,progtag
         REAL*8 uNnet,uEnet,uZnet
   
  -      call gcmdln(verbose)
  -C----
  -C Input files
  -C----
  -      call checkctrlfiles(srcfile,stafile,haffile,verbose)
  -      open (unit=21,file=srcfile,status='old')
  -      open (unit=22,file=stafile,status='old')
  -      open (unit=23,file=haffile,status='old')
  -      read (23,*) vp,vs,dens
  +C----
  +C Parse the command line, check for FFM file, and read FFM
  +C----
  +      call gcmdln(fltfile,stafile,haffile,dspfile,auto)
  +
   C----
   C Output file
   C----
  
positive.code:           evdp(flt) = 1.0d3*evdp(flt)
            dx(flt) = 1.0d3*dx(flt)
            dy(flt) = 1.0d3*dy(flt)
            area(flt) = dx(flt)*dy(flt)     
            goto 201
    202 continue
        nflt = flt - 1
        
        RETURN
        END
        
  C----------------------------------------------------------------------C
  
        SUBROUTINE calcdisp(uNnet,uEnet,uZnet,stlo,stla,stdp,nflt,                                
       1                    evlo,evla,evdp,str,dip,rak,dx,dy,area,slip,                           
       2                    vp,vs,dens,flttyp)
  C----                                                                                           
  C Compute NEZ displacements at (stlo,stla,stdp) from faults
  C----                                                                                           
        IMPLICIT none
        REAL*8 pi,d2r
        PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
  
negative.code:           goto 101
    103 continue
  
   9999 format (2(f10.3),3(f14.4))
  
        END
  
  C======================================================================C
  
        SUBROUTINE readfaults(nflt,evlo,evla,evdp,str,dip,rak,slip,dy,dx,
       1                      area,flttyp,FMAX)
  
        IMPLICIT none
        INTEGER flt,nflt,FMAX,flttyp(FMAX)
        REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),
       1       dip(FMAX),rak(FMAX),slip(FMAX),dy(FMAX),dx(FMAX),
       2       area(FMAX)
  
        flt = 0
    201 flt = flt + 1
        read (21,*,END=202) evlo(flt),evla(flt),evdp(flt),str(flt),
  