anchor.diff:            call usage('!! Use -smooth CONST FILE to specify input file')
         else
             inquire(file=smoof,EXIST=ex)
  -          if (.not.ex) call usage('!! Error: no input file: '//smoof)
  +          if (.not.ex) call usage('!! Error: no input file: '//
  +     1                                                      trim(smoof))
         endif
         call linect(nsmoo,smoof)
         if (nsmoo.gt.nflt) then
  
positive.code: 
  C----------------------------------------------------------------------C
  
        SUBROUTINE geo2xy(flt,obs,nflt,nobs)
  C----
  C Convert longitude and latitude to x and y
  C----
        IMPLICIT none
        INTEGER OBSMAX,FLTMAX
        PARAMETER (OBSMAX=2500,FLTMAX=2000)
        REAL*8 obs(OBSMAX,6),flt(FLTMAX,7)
        INTEGER i,nobs,nflt
        REAL*8 dist,az,lon,lat,lon0,lat0
  
negative.code: 
  C----------------------------------------------------------------------C
  
        SUBROUTINE green(gf,obs,nobs,OBSMAX,flt,nflt,FLTMAX,haf)
  C----
  C Compute and store static Green's functions
  C----
        IMPLICIT none
        INTEGER OBSMAX,FLTMAX
        REAL*8 obs(OBSMAX,6),flt(FLTMAX,7)
        REAL*8 gf(OBSMAX,FLTMAX,6)
        INTEGER nobs,nflt
        REAL*8 haf(3)
        REAL*8 stlo,stla,stdp
        REAL*8 evlo,evla,evdp,str,dip,rak,slip,wid,len
        INTEGER i,j
        REAL*8 vp,vs,dens
  C Compute vp, vs, dens from Lame parameters
        vp   = haf(1)
        vs   = haf(2)
        dens = haf(3)
  C Unit slip for GF computation
        slip = 1.0d0
  C GF for each flt-sta pair
        do 302 i = 1,nobs
            stlo = obs(i,1)
            stla = obs(i,2)
            stdp = obs(i,3)
            do 301 j = 1,nflt
                evlo = flt(j,1)
                evla = flt(j,2)
                evdp = flt(j,3)
                str  = flt(j,4)
                dip  = flt(j,5)
                wid  = flt(j,6)
                len  = flt(j,7)
                rak  = 0.0d0 ! strike-slip GF
                call calcdisp(gf(i,j,1),gf(i,j,2),gf(i,j,3),stlo,stla,
       1                      stdp,evlo,evla,evdp,str,dip,rak,slip,wid,
       2                      len,vp,vs,dens)
                rak  = 90.0d0 ! dip-slip GF
                call calcdisp(gf(i,j,4),gf(i,j,5),gf(i,j,6),stlo,stla,
       1                      stdp,evlo,evla,evdp,str,dip,rak,slip,wid,
       2                      len,vp,vs,dens)
    301     continue
    302 continue
        RETURN
        END
  
  C----------------------------------------------------------------------C
  
        SUBROUTINE calcdisp(ux,uy,uz,stlo,stla,stdp,evlo,evla,evdp,str,
       1                    dip,rak,slip,wid,len,vp,vs,dens)
  C----
  C Compute static displacement vector given a source-station pair
  C All units are SI; angles inputs are in degrees
  C----
        IMPLICIT none
        REAL*8 pi,d2r
        PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
        REAL*8 ux,uy,uz
        REAL*8 stlo,stla,stdp
        REAL*8 evlo,evla,evdp
        REAL*8 str,dip,rak
        REAL*8 slip,wid,len
        REAL*8 vp,vs,dens
        REAL*8 delx,dely,dist,az,x,y
        REAL*8 uxp,uyp,theta,uhor
  C Distance and azimuth from source to station
        delx = stlo - evlo
        dely = stla - evla
        dist = dsqrt(delx*delx+dely*dely)
        az   = datan2(delx,dely) ! clockwise from north
  C Rotate to x=along-strike, y=horizontal up-dip
        x = dist*( dcos(az-d2r*str))
        y = dist*(-dsin(az-d2r*str))
  C Compute displacement
        call o92rect(uxp,uyp,uz,x,y,stdp,evdp,dip,rak,wid,len,slip,vp,vs,
       1             dens)
  C Rotate back to original x-y coordinates
        theta = datan2(uyp,uxp)
        uhor = dsqrt(uxp*uxp+uyp*uyp)
        theta = d2r*str - theta
        ux = uhor*dsin(theta)
        uy = uhor*dcos(theta)
        RETURN
        END
  
  C----------------------------------------------------------------------C
  
        SUBROUTINE lstsqr(soln,gf,obs,nobs,OBSMAX,nflt,FLTMAX,damp,smooth,
       1                  smoof)
  C----
  C Compute a least squares inversion for slip and rake
  C----
        IMPLICIT none
        CHARACTER*80 smoof
        REAL*8 damp,smooth
        INTEGER OBSMAX,FLTMAX,dimx,dimy
        REAL*8 obs(OBSMAX,6)
        REAL*8 gf(OBSMAX,FLTMAX,6)
        REAL*8 A(3*OBSMAX+2*FLTMAX,2*FLTMAX)
        REAL*8 b(3*OBSMAX+2*FLTMAX,1)
        REAL*8 x(3*OBSMAX+2*FLTMAX,1)
        REAL*8 soln(FLTMAX,2)
        INTEGER i,j,nobs,nflt
        INTEGER n,m ! n=nrows, m=ncol
  C----
  C Load model matrix, A, and observation vector, b
  C----
        do 402 i = 1,nobs
            do 401 j = 1,nflt
                A(i,j)             = gf(i,j,1) ! ssx
                A(i+nobs,j)        = gf(i,j,2) ! ssy
                A(i+2*nobs,j)      = gf(i,j,3) ! ssz
                A(i,j+nflt)        = gf(i,j,4) ! dsx
                A(i+nobs,j+nflt)   = gf(i,j,5) ! dsy
                A(i+2*nobs,j+nflt) = gf(i,j,6) ! dsz
    401     continue
    402 continue
        do 403 i = 1,nobs
            b(i,1)        = obs(i,4)
            b(i+nobs,1)   = obs(i,5)
            b(i+2*nobs,1) = obs(i,6)
    403 continue
  C Change number of rows, columns
        n = 3*nobs
        m = 2*nflt
  C----
  C Add damping
  C----
        call damping(A,b,n,m,damp)
  C----
  C Add smoothing
  C----
        call smoothing(A,b,n,m,smooth,smoof)
  C----
  C Solve generalized least squares problem
  C----
        dimx = 3*OBSMAX+2*FLTMAX
        dimy = 2*FLTMAX
        call lsqsolve(A,b,x,dimx,dimy,n,m)
        do 404 i = 1,nflt
            soln(i,1) = x(i,1)
            soln(i,2) = x(i+nflt,1)
    404 continue
        RETURN
        END
  
  