anchor.diff:          end do
       end do
   
  -    K = 0.0 * K
  +    ! Zero output rate constant matrix
  +    do i = 1, nIsom+nProd
  +        do j = 1, nIsom+nProd
  +            K(i,j) = 0.0
  +        end do
  +    end do
   
       ! Calculate the eigenvalues and eigenvectors
  -    allocate( V0(1:nRows), work(1:6*nRows) )
  -    call DSYEV('V', 'U', nRows, M, nRows, V0, work, 6*nRows, info)
  +    allocate( V(1:nRows), work(1:6*nRows) )
  +    call DSYEV('V', 'U', nRows, M, nRows, V, work, 6*nRows, info)
       if (info > 0) then
           msg = 'DSYEV eigenvalue algorithm failed to converge.'
           return
  
positive.code:         return
      end if
  
      ! Extract the chemically-significant eigenvalues and eigenvectors
      ! These are the last nIsom+nProd columns of M
      allocate( X(1:nIsom+nProd, 1:nIsom+nProd), V(1:nIsom+nProd), Xinv(1:nIsom+nProd, 1:nIsom+nProd))
      do j = 1, nIsom + nProd
          count = nRows - (nIsom + nProd) + j
          V(j) = V0(count)
          do r = 1, nGrains
              do i = 1, nIsom
                  index = indices(r,i)
                  if (index > 0) then
                      X(i,j) = X(i,j) + Sinv(index) * M(index,count) * S(count)
                  end if
              end do
          end do
          do i = 1, nProd
              index = nRows - nProd + i
              X(i+nIsom,j) = X(i+nIsom,j) + Sinv(index) * M(index,count) * S(count)
          end do
      end do
      
      ! Invert the chemical eigenvector matrix
      Xinv = X
      allocate(iPiv(1:nIsom+nProd))
      call DGETRF(nIsom+nProd, nIsom+nProd, Xinv, nIsom+nProd, ipiv, info)
      if (info > 0) then
          msg = 'Chemical eigenvector matrix is singular.'
          return
      elseif (info < 0) then
          msg = 'Illegal argument passed to DGETRF.'
          return
      end if
      call DGETRI(nIsom+nProd, Xinv, nIsom+nProd, ipiv, work, 6*nRows, info)
      if (info > 0) then
          msg = 'Chemical eigenvector matrix is singular.'
          return
      elseif (info < 0) then
          msg = 'Illegal argument passed to DGETRI.'
          return
      end if
  
      ! Construct the rate constant matrix
      do i = 1, nIsom+nProd
          do j = 1, nIsom+nProd
              K(i,j) = 0.0
              do r = 1, nIsom+nProd
                  K(i,j) = K(i,j) + X(i,r) * V(r) * Xinv(r,j)
              end do
          end do
      end do
  
      ! Clean up
      deallocate( Mcoll, Mrxn, M, S, Sinv, V0, work, iPiv, X, V, Xinv )
  
  end subroutine
  
  
negative.code:     integer, dimension(1:nGrains,1:nIsom) :: indices
      
      ! Intermediate matrices
      real(8), dimension(:,:), allocatable :: Mcoll, Mrxn, M, p0, dXij
      real(8), dimension(:), allocatable :: S, Sinv, V, Cj
      
      ! BLAS and LAPACK temporary variables
      real(8), dimension(:), allocatable :: work
      integer info
  
      ! Indices
      integer i, j, r, index, start
  
      ! Construct accounting matrix
  	! Row is grain number, column is well number, value is index into matrix
  