diff --git a/src/util/util_solve_linear_system.f90 b/src/util/util_solve_linear_system.f90
index 1f12c0711..b52f08796 100644
--- a/src/util/util_solve_linear_system.f90
+++ b/src/util/util_solve_linear_system.f90
@@ -20,6 +20,7 @@ function util_solve_linear_system(A,b,n,lerr) result(x)
    real(QP), dimension(:), allocatable :: qx
 
    qx = solve_wbs(ge_wpp(real(A, kind=QP), real(b, kind=QP)),lerr)
+   where(abs(qx(:)) < tiny(1._DP)) qx(:) = 0._QP
    x = real(qx, kind=DP)
    return
 
@@ -37,11 +38,15 @@ function solve_wbs(u, lerr) result(x) ! solve with backward substitution
          integer(I4B)             :: i,n
          real(DP), parameter :: epsilon = 10 * tiny(1._DP) 
 
+         if (allocated(x)) deallocate(x)
+         allocate(x(n))
+
          lerr = any(abs(u(:,:)) < epsilon)
-         if (lerr) return
+         if (lerr) then
+            x(:) = 0._DP
+            return
+         end if
 
-         n = size(u,1)
-         allocate(x(n))
          do i = n,1,-1 
             x(i) = (u(i, n + 1) - sum(u(i, i + 1:n) * x(i + 1:n))) / u(i, i)
          end do