anchor.diff:                  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+nReac
  -        write (*,*) i, K(i,i), - sum(K(1:i-1,i)) - sum(K(i+1:nIsom+nReac+nProd,i))
       end do
   
   end subroutine
  
positive.code:     end if
      deallocate( iPiv )
  
      ! Convert solution to pseudo-steady state populations
      do r = minval(nRes)+1, nGrains
          do n = 1, nIsom+nReac
  
negative.code:     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
          call activeStateFull(T, P, E, Mcoll, densStates, Kij, Fim, Gnj, &
  