anchor.diff:  
         SUBROUTINE calcdisp(uNnet,uEnet,uZnet,stlo,stla,stdp,nflt,
        1                    evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  -     2                    vp,vs,dens,flttyp)
  +     2                    vp,vs,dens,flttyp,xy)
   C----
   C Compute north, east, and vertical (positive up) displacements at
   C (stlo,stla,stdp) from shear dislocations.
  
positive.code:           write(STO,*)
        endif
        write(STO,*)
       1 '-fn          (default)  Treat subfaults as finite sources'
        write(STO,*)
       1 '-pt                     Treat subfaults as point sources'
        if (str.eq.'long') then
            write(STO,*)
        endif
  
negative.code: 
  C----------------------------------------------------------------------C
  
        SUBROUTINE getautostnlims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
       1                          evlo,evla,evdp,str,dip,rak,dx,dy,slip,
       2                          nflt,xy,depth)
  C----
  C Get limits of receiver grid by taking steps along longitude and
  C latitudes until strains smaller than minimum threshold value.
  C Variable 'dir' determines direction: 1=N, 2=S, 3=E, 4=W
  C----
        IMPLICIT NONE
        INTEGER i,NTHR
        PARAMETER (NTHR=3)
        REAL*8 dl(NTHR),thr(NTHR)
        REAL*8 hylo,hyla,WMAX,EMAX,SMAX,NMAX,W,E,S,N,lo,la
        REAL*8 stn(3,3),stdp,vp,vs,dens,stnmax,stnmin,depth
        INTEGER nflt,FMAX
        PARAMETER (FMAX=1500)
        REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
       1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
        INTEGER dir,flttyp,xy
        if (xy.eq.0) then
            dl(1) = 1.0d0 ! Amounts to move test location (deg)
            dl(2) = 0.1d0
            dl(3) = 0.05d0
        else
            dl(1) = 100.0d0 ! Amounts to move test location (km)
            dl(2) = 10.0d0
            dl(3) = 5.0d0
        endif
        thr(1) = 1.0d-7 ! Strain thresholds (m)
        thr(2) = 5.0d-8
        thr(3) = 1.0d-8
        flttyp = 1
        call autohalfspace(vp,vs,dens)
        stdp = depth
        do 913 dir = 1,4
            lo = hylo
            la = hyla
  C         Take first step (by increment of dl(1))
            if (dir.eq.1) la = dnint((la+dl(1))*1.0d2)*1.0d-2
            if (dir.eq.2) la = dnint((la-dl(1))*1.0d2)*1.0d-2
            if (dir.eq.3) lo = dnint((lo+dl(1))*1.0d2)*1.0d-2
            if (dir.eq.4) lo = dnint((lo-dl(1))*1.0d2)*1.0d-2
  C         Take steps of decreasing size and threshold value
            i = 1
    911     if (i.gt.NTHR) goto 912
  C             Check if lon, lat exceed maximum allowed values
                if (dir.eq.1.and.la.gt.NMAX) then
                    la = NMAX
                    goto 912
                elseif (dir.eq.2.and.la.lt.SMAX) then
                    la = SMAX
                    goto 912
                elseif (dir.eq.3.and.lo.gt.EMAX) then
                    lo = EMAX
                    goto 912
                elseif (dir.eq.4.and.lo.lt.WMAX) then
                    lo = WMAX
                    goto 912
                endif
  C             Calculate strain and compare to current threshold
                call calcstn(stn,lo,la,stdp,nflt,evlo,evla,evdp,
       1                     str,dip,rak,dx,dy,slip,vp,vs,dens,
       2                     flttyp,xy)
                stnmax = dmax1(stn(1,1),stn(2,2),stn(3,3),stn(1,2),
       1                       stn(1,3),stn(2,3))
                stnmin = dmin1(stn(1,1),stn(2,2),stn(3,3),stn(1,2),
       1                       stn(1,3),stn(2,3))
                stnmax = dmax1(dabs(stnmax),dabs(stnmin))
                if (stnmax.gt.thr(i)) then
                    if (dir.eq.1) la = la + dl(i)
                    if (dir.eq.2) la = la - dl(i)
                    if (dir.eq.3) lo = lo + dl(i)
                    if (dir.eq.4) lo = lo - dl(i)
                    goto 911
                else
                    i = i + 1
                    goto 911
                endif
    912     continue
            if (dir.eq.1) N = la
            if (dir.eq.2) S = la
            if (dir.eq.3) E = lo
            if (dir.eq.4) W = lo
    913 continue
        RETURN
        END
  
  C----------------------------------------------------------------------C
  
        SUBROUTINE writestafile(stafile,W,E,S,N,incr,depth)
  C----
  C Write receiver grid to STAFILE.
  C----
        IMPLICIT NONE
        CHARACTER*40 stafile
        REAL*8 W,E,S,N,incr,lat,lon,stdp,depth
        INTEGER i,j
        stdp = depth
        open(unit=95,file=stafile,status='unknown')
        i = 0
        lon = W
  