anchor.diff:      real(8), dimension(1:nGrains, 1:nIsom+nReac, 1:nIsom)  ::  pa
   
       ! Indices
  -    integer i, n, r
  +    integer i, j, m, n, r
   
       real(8), dimension(1:nIsom,1:nGrains) :: eqDist
       real(8), dimension(1:nIsom) :: bres
  
positive.code:     do i = 1, nIsom
          bres(i) = sum(eqDist(i,1:nRes(i)))
      end do
      
      ! Determine pseudo-steady state populations of active state
      if (nIsom == 1) then
          ! Not worth it to call banded solver when only one well
  
negative.code:     end do
      ! Rows relating to products
      do n = nReac+1, nReac+nProd
          ! Reaction from isomer or reactant j to product n
          do j = 1, nIsom+nReac
              do i = 1, nIsom
                  K(n+nIsom,j) = K(n+nIsom,j) + sum(Gnj(n,i,nRes(i)+1:nGrains) * pa(nRes(i)+1:nGrains,j,i))
              end do
          end do
      end do
  
      do i = 1, nIsom
          K(:,i) = K(:,i) * bres(i)
      end do
  
  end subroutine
  