anchor.diff:  logical :: iWantTraction
   character(len=16) :: coord_type         ! cartesian-m, cartesian-km, geographic
   integer :: nfaults
  +integer :: nstations
   double precision, allocatable :: faults(:,:)
  -integer :: nstations
   double precision, allocatable :: stations(:,:)
  +double precision, allocatable :: targets(:,:)
   double precision :: poisson, lame, mu
  +integer :: verbosity
   
   end module
   
  
positive.code:     elseif (trim(tag).eq.'-stress') then
          i = i + 1
          call get_command_argument(i,stress_file)
          isOutputFileDefined = .true.
          iWantStrain = .true.
          iWantStress = .true.
  
negative.code:     tri_coord(1,3) = faults(iFlt,7)
      tri_coord(2,3) = faults(iFlt,8)
      tri_coord(3,3) = faults(iFlt,9)
      slip(1) = -faults(iFlt,10)        ! triutil: ll positive;  tri_disloc_disp(): rl positive
      slip(2) = -faults(iFlt,11)        ! triutil: thr positive; tri_disloc_disp(): nor positive
      slip(3) = faults(iFlt,12)
  
      if (coord_type.eq.'geographic') then
          ! Convert lon lat dep(km) to x(m) y(m) z(m) from triangle center
          call tri_center(center,tri_coord(:,1),tri_coord(:,2),tri_coord(:,3))
          call ddistaz(dist,az,center(1),center(2),x,y)
          sta_coord(1) = dist*6.371d6*dsin(az)
          sta_coord(2) = dist*6.371d6*dcos(az)
          do iTri = 1,3
              call ddistaz(dist,az,center(1),center(2),tri_coord(1,iTri),tri_coord(2,iTri))
              tri_coord_new(1,iTri) = dist*6.371d6*dsin(az)
              tri_coord_new(2,iTri) = dist*6.371d6*dcos(az)
              tri_coord_new(3,iTri) = tri_coord(3,iTri)*1.0d3
          enddo
      else
          ! Convert triangle x(km) y(km) z(km) to x(m) y(m) z(m)
          tri_coord_new = tri_coord*1.0d3
      endif
  
      call tri_disloc_strain(straintmp, sta_coord, tri_coord_new, poisson, slip)
      strain = strain + straintmp
  enddo
  return
  end subroutine calc_strain
  
  !--------------------------------------------------------------------------------------------------!
  
  subroutine strain2stress(stress,strain,lame_param,shear_modulus)
  implicit none
  ! Arguments
  double precision :: stress(3,3), strain(3,3), lame_param, shear_modulus
  ! Local variables
  double precision :: diag
  diag = strain(1,1) + strain(2,2) + strain(3,3)
  stress(1,1) = lame_param*diag + 2.0d0*shear_modulus*strain(1,1)
  stress(2,2) = lame_param*diag + 2.0d0*shear_modulus*strain(2,2)
  stress(3,3) = lame_param*diag + 2.0d0*shear_modulus*strain(3,3)
  stress(1,2) = 2.0d0*shear_modulus*strain(1,2)
  stress(1,3) = 2.0d0*shear_modulus*strain(1,3)
  stress(2,3) = 2.0d0*shear_modulus*strain(2,3)
  stress(2,1) = stress(1,2)
  stress(3,1) = stress(1,3)
  stress(3,2) = stress(2,3)
  return
  end subroutine strain2stress
  
  !--------------------------------------------------------------------------------------------------!
  
  subroutine calc_estress(estress,stress)
  implicit none
  ! Arguments
  double precision :: estress, stress(3,3)
  ! Local variables
  double precision :: s11_s22, s11_s33, s22_s33, s12s12, s13s13, s23s23
  s11_s22 = stress(1,1)-stress(2,2)
  s11_s33 = stress(1,1)-stress(3,3)
  s22_s33 = stress(2,2)-stress(3,3)
  s12s12 = stress(1,2)*stress(1,2)
  s13s13 = stress(1,3)*stress(1,3)
  s23s23 = stress(2,3)*stress(2,3)
  estress = (1.0d0/6.0d0)*(s11_s22*s11_s22+s11_s33*s11_s33+s22_s33*s22_s33) + s12s12+s13s13+s23s23
  estress = dsqrt(estress)
  return
  end subroutine calc_estress
  
  !--------------------------------------------------------------------------------------------------!
  
  subroutine calc_tractions(trac_vector,shear,shearmax,normal,coulomb,stress,trg)
  !----
  ! Calculate various tractions resolved onto a fault geometry from a stress tensor
  !----
  
  implicit none
  ! Arguments
  double precision :: trac_vector(3), shear, shearmax, normal, coulomb, stress(3,3), trg(4)
  ! Local variables
  integer :: i
  double precision :: str, dip, rak, fric
  double precision :: n(3), s(3), r(3)
  double precision, parameter :: pi=datan(1.0d0)*4.0d0, d2r=pi/180.0d0
  
  str = trg(1)*d2r
  dip = trg(2)*d2r
  rak = trg(3)*d2r
  fric = trg(4)
  
  ! Unit normal vector to target plane
  n(1) = dsin(dip)*dsin(str+pi/2.0d0)
  n(2) = dsin(dip)*dcos(str+pi/2.0d0)
  n(3) = dcos(dip)
  
  ! Traction vector is stress matrix times normal vector: t = S*n
  do i = 1,3
      trac_vector(i) = stress(i,1)*n(1) + stress(i,2)*n(2) + stress(i,3)*n(3)
  enddo
  
  ! Normal component of traction is parallel to unit normal vector; take dot product
  normal = 0.0d0
  do i = 1,3
      normal = normal + trac_vector(i)*n(i)
  enddo
  
  ! Shear component of traction is difference between total traction vector and normal traction
  s(1) = trac_vector(1) - normal*n(1)
  s(2) = trac_vector(2) - normal*n(2)
  s(3) = trac_vector(3) - normal*n(3)
  shearmax = dsqrt(s(1)*s(1) + s(2)*s(2) + s(3)*s(3))
  
  ! Compute unit slip vector (parallel to rake)
  r(1) =  dcos(str-pi/2.0d0)*dcos(rak) + dsin(str-pi/2.0d0)*dsin(rak)*dcos(dip)
  r(2) = -dsin(str-pi/2.0d0)*dcos(rak) + dcos(str-pi/2.0d0)*dsin(rak)*dcos(dip)
  r(3) =                                                    dsin(rak)*dsin(dip)
  
  ! Project shear component on slip vector
  shear = s(1)*r(1) + s(2)*r(2) + s(3)*r(3)
  
  ! Coulomb stress (recall sign convention: pos = dilation)
  coulomb = shear + fric*normal
  
  return
  end subroutine calc_tractions
  
  !--------------------------------------------------------------------------------------------------!
  
  