anchor.diff:      ldb = max(m,n)   ! Leading dimension of b;  ldb >= max(1,m,n)
   
       ! Compute optimal workspace for least squares problem
  +    if (.not.allocated(work)) then
  +        allocate(work(1))
  +    endif
       lwork = -1
  -    call dgels(trans,m,n,nrhs,A,lda,b,ldb,work,lwork,info)
  +    call dgels(trans,m,n,nrhs,Alocal,lda,blocal,ldb,work,lwork,info)
       lwork = int(work(1))
  -
  -    ! Copy displacement vector, b to btmp because it is replaced in dgels
  -    do i = 1,nrows
  -        btmp(i,1) = b(i,1)
  -    enddo
  +    if (verbosity.ge.2) then
  +        write(stderr,'(A,I10)') 'lsqr_solve_dgels says: lwork=',lwork
  +    endif
  +    deallocate(work)
  +    allocate(work(lwork))
   
       ! Solve least squares problem for x
  -    call dgels(trans,m,n,nrhs,a,lda,btmp,ldb,work,lwork,info)
  -    ! do i = 1,nrows
  -    !     write(0,*) (A(i,j),j=1,ncols)
  -    ! enddo
  +    call dgels(trans,m,n,nrhs,alocal,lda,btmp,ldb,work,lwork,info)
  +    deallocate(work)
   
       ! Put solution into x
       do i = 1,ncols
  
positive.code:     implicit none
      ! Local variables
      integer :: i, j, n, dn, ndsp, nflt, nsmo
  
      if (verbosity.ge.2) then
          write(stderr,'(A)') "load_slip_constraints says: starting"
      endif
  
      nflt = fault%nrecords
      ndsp = displacement%nrecords
      nsmo = smoothing%nrecords
  
      ! Move contribution from constrained fault segments to observation side of equation (do not
      ! renumber yet.....will do that later in subroutine)
      do i = 1,nflt
          ! Is the strike-slip (or constant rake) component of this fault constrained?
          if (dabs(slip_constraint%array(i,1)).lt.99998.0d0) then
  
negative.code:             endif
  
              ! Indicate this column has been adjusted by making first entry ridiculous
              cols2delete(i+nflt) = 1
          endif
      enddo
  
      ! Renumber columns of arrays
      if (verbosity.ge.2) then
          write(stderr,'(A)') 'load_slip_constraints says: renumbering matrix columns'
          write(stderr,'(A,I10)') '    ncols=',ncols
      endif
      ! Start counters for existing columns and number of shifted columns
      j = 0
      n = 0
      do i = 1,ncols
          ! If this slip value has not been fixed, then add it to the modified array
          if (cols2delete(i).eq.0) then
              j = j + 1
              A(:,j) = A(:,i)
          else
              n = n + 1
          endif
      enddo
      ncols = ncols - n
      if (verbosity.ge.2) then
          write(stderr,'(A)') 'load_slip_constraints says: new number of columns:'
          write(stderr,'(A,I10)') '    ncols=',ncols
      endif
  
      ! Renumber rows of arrays
      if (verbosity.ge.2) then
          write(stderr,'(A)') 'load_slip_constraints says: renumbering matrix rows'
          write(stderr,'(A,I10)') '    nrows=',nrows
      endif
      ! Start counters for existing rows and number of shifted rows
      j = 0
      n = 0
      do i = 1,nrows
          ! If this row is retained, add it to the modified array
          if (rows2delete(i).eq.0) then
              j = j + 1
              A(j,:) = A(i,:)
              b(j,:) = b(i,:)
          else
              n = n + 1
          endif
      enddo
      nrows = nrows - n
      if (verbosity.ge.2) then
          write(stderr,'(A)') 'load_slip_constraints says: new number of rows:'
          write(stderr,'(A,I10)') '    nrows=',nrows
      endif
  
      if (verbosity.ge.2) then
          write(stderr,'(A)') "load_slip_constraints says: finished"
  