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: 
      ! Allocate matrices and vectors used in this subroutine
      allocate( A(1:nIsom, 1:nIsom) )
      allocate( b(1:nIsom) )
      allocate( iPiv(1:nIsom) )
      allocate( pa(1:nGrains, 1:nIsom) )
  
      ! Zero the phenomenological rate coefficient matrix
      do i = 1, nIsom + nReac + nProd
  
negative.code:     ! Set msg to successful; will be changed later if an error occurs
      msg = ""
  
      ! Determine the starting grain for the calculation based on the
      ! active-state cutoff energy
      start = ceiling((minval(Eres) - minval(E)) / (E(2) - E(1))) + 1
      if (start < 1 .or. start > nGrains) then
          msg = 'Unable to determine starting grain; 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
              do n = 1, nIsom+nReac
                  pa(r,i,n) = 0.0
              end do
          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
              do n = 1, nIsom+nReac
                  b(i,n) = 0.0
              end do
          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
  
          ! Populate RHS vectors, one per isomer and reactant
          do n = 1, nIsom + nReac
              if (n <= nIsom) then
                  ! Thermal activation via collisions
                  b(n,n) = collFreq(n) * densStates(n, r) * exp(-E(r) / 8.314472 / T)
              else
                  ! Chemical activation via association reaction
                  do j = 1, nIsom
                      b(j,n) = Fim(j,n-nIsom,r)
                  end do
              end if
          end do
  
          ! Solve for steady-state population
          call DGESV( nIsom, nIsom+nReac, 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
              do n = 1, nIsom+nReac
                  if (pa(r,i,n) < 0.0) then
                      msg = "A negative steady-state concentration was encountered."
                      return
                  end if
              end do
          end do
  
      end do
  
  
      do src = 1, nIsom + nReac
  
          ! 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,src))
                  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,src))
                      K(n+nIsom,src) = K(n+nIsom,src) + val
                      K(src,src) = K(src,src) - val
                  end if
  