anchor.diff:      call usage('read_faults: no fault file named '//trim(fault_file)//' found')
   endif
   
  +! Count the number of lines in the file and allocate memory for the fault array
  +! Fault array: x1 y1 z1 x2 y2 z2 x3 y3 z3 ss ds ts
   call line_count(nfaults,fault_file)
   allocate(faults(nfaults,12))
   
  +! Read fault file
   open(unit=12,file=fault_file,status='old')
   do i = 1,nfaults
       read(12,'(A)') input_line
  
positive.code: 
  subroutine gcmdln()
  use triutil_module, only: fault_file, station_file, target_file, halfspace_file, &
                            displacement_file, strain_file, stress_file, normal_file, shear_file, &
                            coulomb_file, coord_type, &
               isFaultFileDefined, isStationFileDefined, isTargetFileDefined, isOutputFileDefined, &
               iWantDisp, iWantStrain, iWantStress, iWantTraction, poisson, lame, mu
  implicit none
  ! Local variables
  integer :: i, narg
  character(len=512) :: tag
  
  fault_file = ''
  
negative.code: 
  !--------------------------------------------------------------------------------------------------!
  
  subroutine calc_strain(strain,x,y,z)
  use tri_disloc_module, only: tri_disloc_strain, tri_center
  use triutil_module, only: nfaults, faults, poisson, coord_type
  implicit none
  ! Arguments
  double precision :: strain(3,3), x, y, z
  ! Local variables
  double precision :: sta_coord(3), tri_coord(3,4), slip(3), tri_coord_new(3,4)
  double precision :: center(3), dist, az
  double precision :: straintmp(3,3)
  integer :: iFlt, iTri
  
  ! Initialize displacement
  strain = 0.0d0
  
  ! Set coordinates to meters (input units are geographic and/or km)
  if (coord_type.eq.'cartesian') then
      sta_coord(1) = x*1.0d3
      sta_coord(2) = y*1.0d3
  endif
  sta_coord(3) = z*1.0d3
  
  do iFlt = 1,nfaults
      tri_coord(1,1) = faults(iFlt,1)
  