anchor.diff:      character(len=128), intent(out) :: msg
   
       ! Steady-state populations
  -    real(8), dimension(:,:), allocatable        ::  pa
  +    real(8), dimension(:,:,:), allocatable      ::  pa
       
       ! Steady-state matrix and vector
       real(8), dimension(:,:), allocatable        ::  A
  -    real(8), dimension(:), allocatable          ::  b
  +    real(8), dimension(:,:), allocatable        ::  b
       ! Indices i and j represent sums over unimolecular wells
       integer                                     ::  i, j
       ! Indices m and n represent sums over bimolecular sources/sinks
  
positive.code:     ! Set msg to successful; will be changed later if an error occurs
      msg = ""
  
      ! Iterate over all isomers and channels in the network
      ! Each iteration fills in one column of the phenomenological rate  
      ! coefficient matrix with the rates of that isomer/channel to all others
      ! in the network
      do src = 1, nIsom + nReac
          
          ! Determine the starting grain for the calculation based on the
          ! active-state cutoff energy
          start = ceiling((Eres(src) - minval(E)) / (E(2) - E(1))) + 1
          if (start < 1 .or. start > nGrains) then
              msg = 'Unable to determine starting grain for an isomer; check active-state energies.'
              return
          end if
          
          ! Zero the steady-state population vector for grains below the cutoff
          do r = 1, start-1
              do i = 1, nIsom
                  pa(r,i) = 0.0
          end do
          end do
  
          ! Iterate over the grains, calculating the PSSA concentrations
          do r = start, nGrains
  
              ! Zero A matrix and b vector
              do i = 1, nIsom
                  do j = 1, nIsom
                      A(i,j) = 0.0
                  end do
                  b(i) = 0.0
              end do
              
              ! Collisional deactivation
              do i = 1, nIsom
                  A(i,i) = A(i,i) - collFreq(i)
              end do
              
              ! Isomerization reactions
              do i = 1, nIsom
                  do j = 1, i - 1
                      A(i,j) = Kij(i,j,r)
                      A(j,j) = A(j,j) - Kij(i,j,r)
                      A(j,i) = Kij(j,i,r)
                      A(i,i) = A(i,i) - Kij(j,i,r)
                  end do
              end do
              
              ! Dissociation reactions
              do n = 1, nReac+nProd
                  do j = 1, nIsom
                      A(j,j) = A(j,j) - Gnj(n,j,r)
                  end do
              end do
              
              ! Activation
              if (src <= nIsom) then
                  ! Thermal activation via collisions
                  b(src) = collFreq(src) * densStates(src, r) * exp(-E(r) / 8.314472 / T)
              else
                  ! Chemical activation via association reaction
                  do j = 1, nIsom
                      b(j) = b(j) + Fim(j,src-nIsom,r)
                  end do
              end if
  
              ! Solve for steady-state population
              call DGESV( nIsom, 1, A, nIsom, iPiv, b, nIsom, info )
              if (info > 0) then
                  msg = "A singular matrix was encountered."
                  return
              end if
              pa(r,:) = -b
  
              ! Check that our populations are all positive
              do i = 1, nIsom
                  if (pa(r,i) < 0.0) then
                      msg = "A negative steady-state concentration was encountered."
                      return
                  end if
              end do
  
          end do
  
          ! Calculate stabilization rates (i.e.) R + R' --> Ai or M --> Ai
          do i = 1, nIsom
              if (i /= src) then
                  val = collFreq(i) * sum(pa(:,i))
                  K(i,src) = K(i,src) + val
                  K(src,src) = K(src,src) - val
              end if
          end do
          
          ! Calculate dissociation rates (i.e.) R + R' --> Bn + Cn or M --> Bn + Cn
          do n = 1, nReac+nProd
              do j = 1, nIsom
                  if (n+nIsom /= src) then
                      val = sum(Gnj(n,j,:) * pa(:,j))
                      K(n+nIsom,src) = K(n+nIsom,src) + val
                      K(src,src) = K(src,src) - val
                  end if
  
negative.code: 
      ! Allocate matrices and vectors used in this subroutine
      allocate( A(1:nIsom, 1:nIsom) )
      allocate( b(1:nIsom, 1:nIsom+nReac) )
      allocate( iPiv(1:nIsom) )
      allocate( pa(1:nGrains, 1:nIsom, 1:nIsom+nReac) )
  
      ! Zero the phenomenological rate coefficient matrix
      do i = 1, nIsom + nReac + nProd
  