anchor.diff:  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
   
  
positive.code: 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----
        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)
  C----
  C Calculate net displacements at each station
  C----
  
negative.code: 
  C======================================================================C
  
        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)
        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 = f - 1
        close(21)
  
        RETURN
        END
  