anchor.diff:  !--------------------------------------------------------------------------------------------------!
   
   subroutine strain2stress(stress,strain,lame_param,shear_modulus)
  -implicit none
  +!----
  +! Calculate stress tensor from strain tensor, assuming linear isotropic medium
  +!----
  +
  +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)
  +
  +! Diagonal terms
   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)
  +
  +! Off-diagonal terms
   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
  +!----
  +! Calculate maximum shear 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
  +
  +! Make sure it has units of Pa
   estress = dsqrt(estress)
  +
   return
   end subroutine calc_estress
   
  
positive.code: endif
  
  ! Count number of target geometries and allocate memory to target array (str dip rak fric)
  call line_count(ntargets,target_file)
  if (ntargets.ne.1.and.ntargets.ne.nstations) then
      call usage('read_targets: number of target geometries must be 1 or nstations')
  endif
  
negative.code: 
  !--------------------------------------------------------------------------------------------------!
  
  subroutine gcmdln()
  
  use io, only: stderr, verbosity
  use triutil, only: fault_file, &
                     station_file, &
                     target_file, &
                     halfspace_file, &
                     displacement_file, &
                     strain_file, &
                     stress_file, &
                     estress_file, &
                     normal_file, &
                     shear_file, &
                     coulomb_file, &
                     coord_type, &
                     isFaultFileDefined, &
                     isStationFileDefined, &
                     isTargetFileDefined, &
                     isOutputFileDefined, &
                     iWantDisp, &
                     iWantStrain, &
                     iWantStress, &
                     iWantTraction, &
                     poisson, &
                     lame, &
                     shearmod
  implicit none
  
  ! Local variables
  integer :: i, narg, ios
  character(len=512) :: tag
  
  ! Initialize control variables
  fault_file = ''
  station_file = ''
  target_file = ''
  