anchor.diff:  C Parse command line
   C----
         call gcmdln(ffmf,fltf,magf,haff,staf,trgf,dspf,stnf,stsf,norf,
  -     1            shrf,coulf,gmtf,flttyp,auto,incr,depth,xy,prog,long)
  +     1            shrf,coulf,gmtf,flttyp,auto,incr,val,xy,prog,long)
         call dbggcmdln(ffmf,fltf,magf,haff,staf,trgf,dspf,stnf,stsf,norf,
  -     2               shrf,coulf,gmtf,flttyp,auto,incr,depth,xy,prog,
  +     2               shrf,coulf,gmtf,flttyp,auto,incr,val,xy,prog,
        3               long)
   
   C----
  
positive.code:           write(*,*)
        endif
        write(*,*)
       1 '-long                   Increase displacement output precision'
        if (str.eq.'long') then
            write(*,*)
  
negative.code:       print *,'FLTTYP',flttyp
        print *,'AUTO',auto
        print *,'INCR',incr
        print *,'VAL',val
        print *,'XY?',xy
        print *,'PROG',prog
        print *,'LONG',long
        RETURN
        END
  
  C||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||C
  CVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVC
  C THESE ARE ALL THE OLD SUBROUTINES FOR AUTOMATIC OPTIONS
  C
  C      SUBROUTINE autodispgrid(staf,evlo,evla,evdp,str,dip,rak,dx,dy,
  C     1                        slip,hylo,hyla,nflt,incr,val,dspf,xy)
  CC----
  CC Automatically determine receiver grid for displacement computation,
  CC write grid to STAFILE.
  CC----
  C      IMPLICIT NONE
  C      CHARACTER*40 staf,dspf
  C      INTEGER nflt,FMAX,xy
  C      PARAMETER (FMAX=1500)
  C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  C      REAL*8 W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,incr,val
  C      write(*,'(A)') ' Creating receiver grid in file: '//staf
  C      dmax = 750.0d0 ! km
  C      call getmaxlims(WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,xy)
  C      call getautodisplims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  C     1                     evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  C     2                     nflt,xy,val)
  C      call writestaf(staf,W,E,S,N,incr,val)
  C      write(*,'(A)') ' Displacements written to file: '//dspf
  C      RETURN
  C      END
  C
  C
  C      SUBROUTINE autocoulgrid(staf,evlo,evla,evdp,str,dip,rak,dx,dy,
  C     1                        slip,hylo,hyla,nflt,incr,val,coulf,xy)
  CC----
  CC Automatically determine receiver grid for Coulomb stress computation,
  CC write grid to COUFILE.
  CC----
  C      IMPLICIT NONE
  C      CHARACTER*40 staf,coulf
  C      INTEGER nflt,FMAX,xy
  C      PARAMETER (FMAX=1500)
  C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  C      REAL*8 W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,incr,val
  C      write(*,'(A)') ' Creating receiver grid in file: '//staf
  C      dmax = 750.0d0 ! km
  C      call getmaxlims(WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,xy)
  C      call getautostnlims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  C     1                    evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  C     2                    nflt,xy,val)
  C      call writestaf(staf,W,E,S,N,incr,val)
  C      write(*,'(A)') ' Coulomb stress written to file: '//coulf
  C      RETURN
  C      END
  C
  C
  C      SUBROUTINE getmaxlims(WMAX,EMAX,SMAX,NMAX,hylo,hyla,dmax,xy)
  C      IMPLICIT NONE
  C      REAL*8 pi,r2d,d2r
  C      PARAMETER (pi=4.0d0*atan(1.0d0),r2d=180.0d0/pi,d2r=pi/180.0d0)
  C      REAL*8 az,dmax,WMAX,EMAX,SMAX,NMAX,hylo,hyla,junk
  C      INTEGER xy
  CC----
  CC Determine maximum geographic limits of receiver grid to W,E,S,N.
  CC----
  C      if (xy.eq.0) then
  C          WMAX = hylo-dmax*3.6d2/(2.0d0*pi*6.371d3*dcos(hyla*d2r))
  C          EMAX = hylo+dmax*3.6d2/(2.0d0*pi*6.371d3*dcos(hyla*d2r))
  C          az = 180.0d0
  C          call dlola(junk,SMAX,hylo,hyla,dmax,az)
  C          SMAX = SMAX*r2d
  C          az = 0.0d0
  C          call dlola(junk,NMAX,hylo,hyla,dmax,az)
  C          NMAX = NMAX*r2d
  C      else
  C          WMAX = hylo - dmax
  C          EMAX = hylo + dmax
  C          SMAX = hyla - dmax
  C          NMAX = hyla + dmax
  C      endif
  C      RETURN
  C      END
  C
  C
  C      SUBROUTINE getautodisplims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  C     1                           evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  C     2                           nflt,xy,val)
  CC----
  CC Get limits of receiver grid by taking steps along longitude and
  CC latitudes until displacements smaller than minimum threshold value.
  CC Variable 'dir' determines direction: 1=N, 2=S, 3=E, 4=W
  CC----
  C      IMPLICIT NONE
  C      INTEGER i,NTHR
  C      PARAMETER (NTHR=3)
  C      REAL*8 dl(NTHR),thr(NTHR)
  C      REAL*8 hylo,hyla,WMAX,EMAX,SMAX,NMAX,W,E,S,N,lo,la
  C      REAL*8 uE,uN,uZ,stdp,vp,vs,dens,val
  C      INTEGER nflt,FMAX,xy
  C      PARAMETER (FMAX=1500)
  C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  C      INTEGER dir,flttyp
  C      if (xy.eq.0) then
  C          dl(1) = 1.0d0 ! Amounts to move test location (deg)
  C          dl(2) = 0.1d0
  C          dl(3) = 0.05d0
  C      else
  C          dl(1) = 100.0d0 ! Amounts to move test location (km)
  C          dl(2) = 10.0d0
  C          dl(3) = 5.0d0
  C      endif
  C      thr(1) = 0.5d0 ! Displacement thresholds (m)
  C      thr(2) = 0.01d0
  C      thr(3) = 0.001d0
  C      flttyp = 1
  C      call autohaf(vp,vs,dens)
  C      stdp = val
  C      do 903 dir = 1,4
  C          lo = hylo
  C          la = hyla
  C         Take first step (by increment of dl(1))
  C          if (dir.eq.1) la = dnint((la+dl(1))*1.0d2)*1.0d-2
  C          if (dir.eq.2) la = dnint((la-dl(1))*1.0d2)*1.0d-2
  C          if (dir.eq.3) lo = dnint((lo+dl(1))*1.0d2)*1.0d-2
  C          if (dir.eq.4) lo = dnint((lo-dl(1))*1.0d2)*1.0d-2
  C         Take steps of decreasing size and threshold value
  C          i = 1
  C  901     if (i.gt.NTHR) goto 902
  C             Check if lon, lat exceed maximum allowed values
  C              if (dir.eq.1.and.la.gt.NMAX) then
  C                  la = NMAX
  C                  goto 902
  C              elseif (dir.eq.2.and.la.lt.SMAX) then
  C                  la = SMAX
  C                  goto 902
  C              elseif (dir.eq.3.and.lo.gt.EMAX) then
  C                  lo = EMAX
  C                  goto 902
  C              elseif (dir.eq.4.and.lo.lt.WMAX) then
  C                  lo = WMAX
  C                  goto 902
  C              endif
  C             Calculate displacement and compare to current threshold
  C              call calcdisp(uN,uE,uZ,lo,la,stdp,nflt,evlo,evla,evdp,
  C     1                      str,dip,rak,dx,dy,slip,vp,vs,dens,
  C     2                      flttyp,xy)
  C              uN = dabs(uN)
  C              uE = dabs(uE)
  C              uZ = dabs(uZ)
  C              if (uN.gt.thr(i).or.uE.gt.thr(i).or.uZ.gt.thr(i)) then
  C                  if (dir.eq.1) la = la + dl(i)
  C                  if (dir.eq.2) la = la - dl(i)
  C                  if (dir.eq.3) lo = lo + dl(i)
  C                  if (dir.eq.4) lo = lo - dl(i)
  C                  goto 901
  C              else
  C                  i = i + 1
  C                  goto 901
  C              endif
  C  902     continue
  C          if (dir.eq.1) N = la
  C          if (dir.eq.2) S = la
  C          if (dir.eq.3) E = lo
  C          if (dir.eq.4) W = lo
  C  903 continue
  C      RETURN
  C      END
  C
  C
  C      SUBROUTINE getautostnlims(W,E,S,N,WMAX,EMAX,SMAX,NMAX,hylo,hyla,
  C     1                          evlo,evla,evdp,str,dip,rak,dx,dy,slip,
  C     2                          nflt,xy,val)
  CC----
  CC Get limits of receiver grid by taking steps along longitude and
  CC latitudes until strains smaller than minimum threshold value.
  CC Variable 'dir' determines direction: 1=N, 2=S, 3=E, 4=W
  CC----
  C      IMPLICIT NONE
  C      INTEGER i,NTHR
  C      PARAMETER (NTHR=3)
  C      REAL*8 dl(NTHR),thr(NTHR)
  C      REAL*8 hylo,hyla,WMAX,EMAX,SMAX,NMAX,W,E,S,N,lo,la
  C      REAL*8 stn(3,3),stdp,vp,vs,dens,stnmax,stnmin,val
  C      INTEGER nflt,FMAX
  C      PARAMETER (FMAX=1500)
  C      REAL*8 evlo(FMAX),evla(FMAX),evdp(FMAX),str(FMAX),dip(FMAX),
  C     1       rak(FMAX),dx(FMAX),dy(FMAX),slip(FMAX)
  C      INTEGER dir,flttyp,xy
  C      if (xy.eq.0) then
  C          dl(1) = 1.0d0 ! Amounts to move test location (deg)
  C          dl(2) = 0.1d0
  C          dl(3) = 0.05d0
  C      else
  C          dl(1) = 100.0d0 ! Amounts to move test location (km)
  C          dl(2) = 10.0d0
  C          dl(3) = 5.0d0
  C      endif
  C      thr(1) = 1.0d-7 ! Strain thresholds (m)
  C      thr(2) = 5.0d-8
  C      thr(3) = 1.0d-8
  C      flttyp = 1
  C      call autohaf(vp,vs,dens)
  C      stdp = val
  C      do 913 dir = 1,4
  C          lo = hylo
  C          la = hyla
  C         Take first step (by increment of dl(1))
  C          if (dir.eq.1) la = dnint((la+dl(1))*1.0d2)*1.0d-2
  C          if (dir.eq.2) la = dnint((la-dl(1))*1.0d2)*1.0d-2
  C          if (dir.eq.3) lo = dnint((lo+dl(1))*1.0d2)*1.0d-2
  C          if (dir.eq.4) lo = dnint((lo-dl(1))*1.0d2)*1.0d-2
  C         Take steps of decreasing size and threshold value
  C          i = 1
  C  911     if (i.gt.NTHR) goto 912
  C             Check if lon, lat exceed maximum allowed values
  C              if (dir.eq.1.and.la.gt.NMAX) then
  C                  la = NMAX
  C                  goto 912
  C              elseif (dir.eq.2.and.la.lt.SMAX) then
  C                  la = SMAX
  C                  goto 912
  C              elseif (dir.eq.3.and.lo.gt.EMAX) then
  C                  lo = EMAX
  C                  goto 912
  C              elseif (dir.eq.4.and.lo.lt.WMAX) then
  C                  lo = WMAX
  C                  goto 912
  C              endif
  C             Calculate strain and compare to current threshold
  C              call calcstn(stn,lo,la,stdp,nflt,evlo,evla,evdp,
  C     1                     str,dip,rak,dx,dy,slip,vp,vs,dens,
  C     2                     flttyp,xy)
  C              stnmax = dmax1(stn(1,1),stn(2,2),stn(3,3),stn(1,2),
  C     1                       stn(1,3),stn(2,3))
  C              stnmin = dmin1(stn(1,1),stn(2,2),stn(3,3),stn(1,2),
  C     1                       stn(1,3),stn(2,3))
  C              stnmax = dmax1(dabs(stnmax),dabs(stnmin))
  C              if (stnmax.gt.thr(i)) then
  C                  if (dir.eq.1) la = la + dl(i)
  C                  if (dir.eq.2) la = la - dl(i)
  C                  if (dir.eq.3) lo = lo + dl(i)
  C                  if (dir.eq.4) lo = lo - dl(i)
  C                  goto 911
  C              else
  C                  i = i + 1
  C                  goto 911
  C              endif
  C  912     continue
  C          if (dir.eq.1) N = la
  C          if (dir.eq.2) S = la
  C          if (dir.eq.3) E = lo
  C          if (dir.eq.4) W = lo
  C  913 continue
  C      RETURN
  C     END
  C
  C
  C      SUBROUTINE writestaf(staf,W,E,S,N,incr,val)
  CC----
  CC Write receiver grid to STAFILE.
  CC----
  C      IMPLICIT NONE
  C      CHARACTER*40 staf
  C      REAL*8 W,E,S,N,incr,lat,lon,stdp,val
  C      INTEGER i,j
  C      stdp = val
  C      open(unit=95,file=staf,status='unknown')
  C      i = 0
  C      lon = W
  C  951 if (lon.gt.E) goto 954
  C          i = i + 1
  C          lat = S
  C          j = 0
  C  952     if (lat.gt.N) goto 953
  C              write(95,8001) lon,lat,stdp
  C              lat = lat + incr
  C              j = j + 1
  C              goto 952
  C  953     continue
  C          lon = lon + incr
  C          goto 951
  C  954 continue
  C      open(unit=96,file='autogrid.dat',status='unknown')
  C      write(*,9001) W,E,S,N
  C      write(*,9002) incr
  C      write(*,9003) i*j
  C      write(*,9004) i
  C      write(*,9005) j
  C      write(96,9001) W,E,S,N
  C      write(96,9002) incr
  C      write(96,9003) i*j
  C      write(96,9004) i
  C      write(96,9005) j
  C      write(*,*) 'This information is saved in autogrid.dat'
  C 8001 format(3F10.2)
  C 9001 format(' Grid limits (WESN):',4F8.2)
  C 9002 format(' Grid spacing: ',F10.2)
  C 9003 format(' Number of points: ',I10)
  C 9004 format(' nlon: ',I10)
  C 9005 format(' nlat: ',I10)
  C      close(95)
  C      close(96)
  C      RETURN
  C      END
  
  