anchor.diff:              call get_command_argument(i,tag)
               char_index = index(tag,'-')
               if (char_index.eq.0) then
  -                read(tag,'(BN,I10)') verbosity
  +                call get_command_argument(i,halfspace%flag)
               else
                   i = i - 1
               endif
           endif
  +    elseif (trim(tag).eq.'-v') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) verbosity
  +
  +    ! Least squares options
  +    elseif (trim(tag).eq.'-lsqr:mode') then
  +        i = i + 1
  +        call get_command_argument(i,lsqr_mode)
  +    elseif (trim(tag).eq.'-lsqr:stress_weight') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) stress_weight
  +
  +    ! Simulated annealing options
  +    elseif (trim(tag).eq.'-anneal:init_mode') then
  +        i = i + 1
  +        call get_command_argument(i,anneal_init_mode)
  +    elseif (trim(tag).eq.'-anneal:it_max') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) max_iteration
  +    elseif (trim(tag).eq.'-anneal:it_reset') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) reset_iteration
  +    elseif (trim(tag).eq.'-anneal:log_file') then
  +        i = i + 1
  +        call get_command_argument(i,anneal_log_file)
  +    elseif (trim(tag).eq.'-anneal:temp_0') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) temp_start
  +    elseif (trim(tag).eq.'-anneal:temp_min') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) temp_minimum
  +    elseif (trim(tag).eq.'-anneal:cool') then
  +        i = i + 1
  +        call get_command_argument(i,tag)
  +        read(tag,*) cooling_factor
  +
  +    ! No option
       else
  -        call usage('!! Error: No option '//trim(tag))
  +        call print_usage('!! Error: No option '//trim(tag))
       endif
       i = i + 1
   enddo
  
positive.code: endif
  if (verbosity.ge.2) then
      write(stderr,'("Parsed command line inputs")')
      write(stderr,'("    output_file:        ",A)') trim(output_file)
      write(stderr,'("    fault_file:         ",A)') trim(fault_file)
      write(stderr,'("    displacement_file:  ",A)') trim(displacement_file)
      write(stderr,'("    prestress_file:     ",A)') trim(prestress_file)
      write(stderr,'("    halfspace_file:     ",A)') trim(halfspace_file)
      write(stderr,'("    smoothing_file:     ",A)') trim(smoothing_file)
      write(stderr,'("    inversion_mode:     ",A)') trim(inversion_mode)
      write(stderr,'("    damping constant:   ",F20.8)') damping_constant
      write(stderr,'("    smoothing constant: ",F20.8)') smoothing_constant
      write(stderr,'("    geographic flag:    ",I20)') geographic
      write(stderr,'("    disp component:     ",A)') disp_comp
      write(stderr,'("    slip constraint:    ",A)') slip_constraint_file
      write(stderr,'("    rake file:          ",A)') trim(rake_file)
      write(stderr,*)
  endif
  if (verbosity.ge.1) then
      write(stderr,*)
  endif
  
  return
  end
  
  !--------------------------------------------------------------------------------------------------!
  
  subroutine usage(string)
  !----
  ! Print program usage statement and exit
  !----
  use io
  implicit none
  character(len=*) :: string
  integer :: string_length
  if (string.ne.'') then
      string_length = len(string)
      write(stderr,*) trim(string)
      write(stderr,*)
  endif
  write(stderr,*) 'Usage: fltinv ...options...'
  write(stderr,*)
  write(stderr,*) '-o|--output OUTPUT_FILE'
  write(stderr,*) '-m|--mode MODE'
  write(stderr,*) '-f|--faults FAULT_FILE'
  write(stderr,*) '-d|--displacements DISP_FILE'
  write(stderr,*) '-s|--prestresses PRESTS_FILE'
  write(stderr,*) '-h|--halfspace HAFSPC_FILE'
  write(stderr,*) '--damping DAMP'
  write(stderr,*) '--smoothing SMOOTH SMOOTH_FILE'
  write(stderr,*) '--geographic'
  write(stderr,*) '--disp_comp [1][2][3]'
  write(stderr,*) '--rake RAKE'
  write(stderr,*) '--slip_constraint CONST_FILE'
  write(stderr,*) '-v|--verbose VRB_LEVEL'
  stop
  end
  
negative.code: program main
  implicit none
  
  call initialize_fltinv_variables()
  call get_command_line()
  
  call read_fltinv_inputs()
  
  call calc_greens_functions()
  
  call run_inversion()
  call write_solution()
  
  end program main
  
  !==================================================================================================!
  
  subroutine initialize_fltinv_variables()
  use io_module, only: verbosity, initialize_program_data
  use variable_module, only: output_file, displacement, disp_components, prestress, &
                             fault, slip_constraint, rake_constraint, &
                             gf_type, gf_disp, gf_stress, &
                             inversion_mode, damping_constant, smoothing_constant, smoothing, &
                             coord_type, halfspace
  use lsqr_module, only: lsqr_mode, stress_weight
  use anneal_module, only: anneal_init_mode, anneal_log_file, max_iteration, reset_iteration, &
                           temp_start, temp_minimum, cooling_factor
  implicit none
  
  ! Initialize program behavior variables
  output_file = 'stdout'
  inversion_mode = 'lsqr'
  gf_type = 'okada_rect'
  coord_type = 'cartesian'
  disp_components = '123'
  
  ! Initialize regularization variables
  damping_constant = -1.0d0
  smoothing_constant = -1.0d0
  
  ! Initialize derived data variable values
  call initialize_program_data(fault)
  call initialize_program_data(displacement)
  call initialize_program_data(prestress)
  call initialize_program_data(gf_disp)
  call initialize_program_data(gf_stress)
  call initialize_program_data(slip_constraint)
  call initialize_program_data(rake_constraint)
  call initialize_program_data(smoothing)
  
  ! Same type of derived data variable, but not always used
  call initialize_program_data(halfspace)
  
  ! Pretty obvious, but initialize program verbosity variable
  verbosity = 0
  
  ! Initialize least-squares variables
  lsqr_mode = 'gels'
  stress_weight = 1.0d-9
  
  ! Initialize annealing variables
  anneal_init_mode = 'mean'
  anneal_log_file = 'none'
  max_iteration = 1000
  reset_iteration = 1000000
  temp_start = 2.0d0
  temp_minimum = 0.00d0
  cooling_factor = 0.98d0
  
  return
  end
  
  !--------------------------------------------------------------------------------------------------!
  
  subroutine get_command_line()
  use io_module, only: stderr, verbosity
  use variable_module, only: output_file, displacement, disp_components, prestress, &
                             fault, slip_constraint, rake_constraint, &
                             gf_type, gf_disp, gf_stress, &
                             inversion_mode, damping_constant, smoothing_constant, smoothing, &
                             coord_type, halfspace
  use lsqr_module, only: lsqr_mode, stress_weight
  use anneal_module, only: anneal_init_mode, anneal_log_file, max_iteration, reset_iteration, &
                           temp_start, temp_minimum, cooling_factor
  implicit none
  ! Local variables
  integer :: i, narg
  character(len=256) :: tag
  integer :: char_index
  
  narg = command_argument_count()
  if (narg.eq.0) call print_usage('')
  
  i = 1
  do while (i.le.narg)
  
      call get_command_argument(i,tag)
  
      ! Output options
      if (trim(tag).eq.'-o') then
          i = i + 1
          call get_command_argument(i,output_file)
  
      ! Input options
      elseif (trim(tag).eq.'-disp') then
          i = i + 1
          call get_command_argument(i,displacement%file)
      elseif (trim(tag).eq.'-disp:components') then
          i = i + 1
          call get_command_argument(i,tag)
          read(tag,*) disp_components
      elseif (trim(tag).eq.'-prests') then
          i = i + 1
          call get_command_argument(i,prestress%file)
      elseif (trim(tag).eq.'-flt') then
          i = i + 1
          call get_command_argument(i,fault%file)
      elseif (trim(tag).eq.'-flt:rake') then
          i = i + 1
          call get_command_argument(i,rake_constraint%file)
      elseif (trim(tag).eq.'-flt:slip') then
          i = i + 1
          call get_command_argument(i,slip_constraint%file)
  
      ! Green's functions options
      elseif (trim(tag).eq.'-gf:model') then
          i = i + 1
          call get_command_argument(i,gf_type)
      elseif (trim(tag).eq.'-gf:disp') then
          i = i + 1
          call get_command_argument(i,gf_disp%file)
      elseif (trim(tag).eq.'-gf:stress') then
          i = i + 1
          call get_command_argument(i,gf_stress%file)
  
      ! Inversion options
      elseif (trim(tag).eq.'-mode') then
          i = i + 1
          call get_command_argument(i,inversion_mode)
      elseif (trim(tag).eq.'-damp') then
          i = i + 1
          call get_command_argument(i,tag)
          read(tag,*) damping_constant
      elseif (trim(tag).eq.'-smooth') then
          i = i + 1
          call get_command_argument(i,tag)
          read(tag,*) smoothing_constant
          i = i + 1
          call get_command_argument(i,smoothing%file)
  
      ! Miscellaneous options
      elseif (trim(tag).eq.'-geo') then
          coord_type = 'geographic'
      elseif (trim(tag).eq.'-haf') then
          i = i + 1
          call get_command_argument(i,halfspace%file)
          i = i + 1
          if (i.gt.narg) then
              return
  