anchor.diff:  
   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)
  +      SUBROUTINE readfaults(fltfile,nflt,evlo,evla,evdp,str,dip,rak,
  +     1                      slip,dy,dx,area,flttyp)
  +
  +      IMPLICIT none
  +      CHARACTER*40 fltfile
  +      INTEGER f,nflt,FMAX
  +      PARAMETER (FMAX=1500)
         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),
  -     1                    dip(flt),rak(flt),slip(flt),dy(flt),dx(flt),
  -     2                    flttyp(flt)
  -          evdp(flt) = 1.0d3*evdp(flt)
  -          dx(flt) = 1.0d3*dx(flt)
  -          dy(flt) = 1.0d3*dy(flt)
  -          area(flt) = dx(flt)*dy(flt)
  +      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
  -      nflt = flt - 1
  +      nflt = f - 1
  +      close(21)
   
         RETURN
         END
  
positive.code: C an isotropic elastic half-space due to shear dislocations.
  C----
        IMPLICIT NONE
        CHARACTER*40 fltfile,stafile,haffile,dspfile
        LOGICAL ex
        INTEGER auto
  
negative.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
  
  C----------------------------------------------------------------------C
  
        SUBROUTINE sethafdef(vp,vs,dens)
        IMPLICIT none
        REAL*8 vp,vs,dens
        vp = 6800.0d0
        vs = 3926.0d0
        dens = 2800.0d0
        RETURN
        END
  
  