anchor.diff:        RETURN
         END
   
  -C----------------------------------------------------------------------C
  +!----------------------------------------------------------------------C
   
         SUBROUTINE usage(str)
         IMPLICIT none
  
positive.code:       REAL*8 shr(2)
        str = strin*d2r
        dip = dipin*d2r
  C Unit normal vector to plane
        n(1) = dsin(dip)*dsin(str+pi/2.0d0)
        n(2) = dsin(dip)*dcos(str+pi/2.0d0)
        n(3) = dcos(dip)
  C Traction is stress matrix times normal vector: t = S*n
        do 241 i = 1,3
            trac(i) = stress(i,1)*n(1)+stress(i,2)*n(2)+stress(i,3)*n(3)
    241 continue
  C Normal component of traction is parallel to unit normal vector
        norml = 0.0d0
        do 242 i = 1,3
            norml = norml + trac(i)*n(i)
    242 continue
  C Shear component of traction is difference between total and normal
        s(1) = trac(1) - norml*n(1)
        s(2) = trac(2) - norml*n(2)
        s(3) = trac(3) - norml*n(3)
  C Compute unit vector in strike and up-dip directions
        ns(1) = sin(str)
        ns(2) = cos(str)
        ns(3) = 0.0d0
        nd(1) = cos(dip)*sin(str-pi*0.5d0)
        nd(2) = cos(dip)*cos(str-pi*0.5d0)
        nd(3) = sin(dip)
  C Project shear traction onto strike and dip directions
        shr(1) = s(1)*ns(1) + s(2)*ns(2) + s(3)*ns(3)
        shr(2) = s(1)*nd(1) + s(2)*nd(2) + s(3)*nd(3)
        RETURN
        END
  
  C----------------------------------------------------------------------C
  
        SUBROUTINE calcstn(strain,stlo,stla,stdp,evlo,evla,evdp,str,dip,
       1                   rak,slip,wid,len,vp,vs,dens)
  
negative.code:       REAL*8 soln(FLTMAX,2)
        INTEGER i,j,nobs,nflt
        INTEGER n,m ! n=nrows, m=ncol
        INTEGER vrb
        COMMON /VERBOSE/ vrb
        n = 0
        m = 2*nflt
  !----
  ! Load model matrix, A, and observation vector, b
  !----
        if (nobs.gt.0) then
            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
            ! Change number of rows, columns
            n = 3*nobs
        endif
  
  !----
  ! Only use horizontal or vertical if specified
  !----
        if (cmpnt.eq.1) then
            ! Only horizontal
            n = 2*nobs
        elseif (cmpnt.eq.2) then
            ! Only vertical
            do i = 1,nobs
                do j = 1,nflt
                    A(i,j) = A(i+2*nobs,j)
                    b(i,1) = b(i+2*nobs,j)
                enddo
            enddo
            n = nobs
        endif
  
  !----
  ! Add stress constraints
  !----
        call addsts(A,b,n,m,sgf,stscon,nsts,OBSMAX,FLTMAX,prests)
  
  !----
  ! Add damping
  !----
        call damping(A,b,n,m,damp,OBSMAX,FLTMAX)
  
  !----
  ! Add smoothing
  !----
        call smoothing(A,b,n,m,smooth,smoof,OBSMAX,FLTMAX)
  
  !----
  ! Print results if verbose mode is on
  !----
        if (vrb.ge.3) then
            write(0,*) 'A ='
            do i = 1,n
                write(0,*) (A(i,j),j=1,m)
            enddo
            write(0,*)
            write(0,*) 'b ='
            do i = 1,n
                write(0,*) b(i,1)
            enddo
        endif
  
  !----
  ! Solve generalized least squares problem
  !----
        dimx = 3*OBSMAX+4*FLTMAX
        dimy = 2*FLTMAX
        call lsqsolve(A,b,x,dimx,dimy,n,m)
  