anchor.diff:                    gf_euler, &
                     displacement, &
                     disp_components, &
  +                  input_disp_unit, &
                     los, &
                     prestress, &
                     cov_matrix, &
  
positive.code: integer :: nflt, nslip, ndsp, ndsp_dof, nlos, nsts, nnbr, nfixed, nfree
  integer :: i, j, ierr, icomp, iflt, inbr
  double precision, allocatable :: A(:,:), b(:), x(:), atmp(:,:), btmp(:)
  double precision :: damping_squared, smoothing_squared
  logical, allocatable :: isSlipFixed(:)
  
  
  ! Initialize solution
  
negative.code:     enddo
  endif
  
  ! Load solution into euler pole vector, scaling by correct units
  factor = 0.0d0
  if (input_disp_unit.eq.'m/s') then
      factor = 1.0d0
  elseif (input_disp_unit.eq.'m/yr') then
      factor = 1.0d0*spy
  elseif (input_disp_unit.eq.'mm/s') then
      factor = 1.0d3
  elseif (input_disp_unit.eq.'mm/yr') then
      factor = 1.0d3*spy
  else
      call usage('invert_lsqr: unit '//trim(input_disp_unit)//' not compatible')
  endif
  do i = 1,npoles
      euler_pole(i,1) = x(j+i-1)/factor
      euler_pole(i,2) = x(j+i  )/factor
      euler_pole(i,3) = x(j+i+1)/factor
  enddo
  
  return
  