diff --git a/src/algebra.F90 b/src/algebra.F90
index 148141d..d3c21e1 100644
--- a/src/algebra.F90
+++ b/src/algebra.F90
@@ -1103,7 +1103,7 @@ module algebra
 !
 !===============================================================================
 !
-  subroutine tridiag(n, l, d, u, r, x)
+  subroutine tridiag(n, l, d, u, r, x, iret)
 
 ! import external procedures
 !
@@ -1118,6 +1118,7 @@ module algebra
     integer                   , intent(in)  :: n
     real(kind=8), dimension(n), intent(in)  :: l, d, u, r
     real(kind=8), dimension(n), intent(out) :: x
+    integer                   , intent(out) :: iret
 
 ! local variables
 !
@@ -1142,7 +1143,8 @@ module algebra
       if (t == 0.0d+00) then
         write(error_unit,"('[',a,']: ',a)") trim(loc)                          &
                         , "Could not find solution!"
-        stop
+        iret = 1
+        return
       end if
       x(i) = (r(i) - l(i) * x(j)) / t
     end do
@@ -1154,6 +1156,10 @@ module algebra
       x(i) = x(i) - g(j) * x(j)
     end do
 
+! set return value to success
+!
+    iret = 0
+
 !-------------------------------------------------------------------------------
 !
   end subroutine tridiag
diff --git a/src/interpolations.F90 b/src/interpolations.F90
index 88a883d..99c6202 100644
--- a/src/interpolations.F90
+++ b/src/interpolations.F90
@@ -2465,6 +2465,7 @@ module interpolations
 ! local variables
 !
     integer      :: i, im1, ip1, im2, ip2
+    integer      :: iret
     real(kind=8) :: bl, bc, br, tt
     real(kind=8) :: wl, wc, wr, ww
     real(kind=8) :: ql, qc, qr
@@ -2777,7 +2778,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations for the left-side interpolation
 !
-    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
+    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret)
 
 ! substitute the left-side values
 !
@@ -2785,7 +2786,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations for the left-side interpolation
 !
-    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
+    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret)
 
 ! substitute the right-side values
 !
@@ -2857,6 +2858,7 @@ module interpolations
 ! local variables
 !
     integer      :: i, im1, ip1, im2, ip2
+    integer      :: iret
     real(kind=8) :: bl, bc, br, tt
     real(kind=8) :: wl, wc, wr, ww
     real(kind=8) :: ql, qc, qr
@@ -3172,7 +3174,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations for the left-side interpolation
 !
-    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
+    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret)
 
 ! substitute the left-side values
 !
@@ -3180,7 +3182,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations for the left-side interpolation
 !
-    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
+    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret)
 
 ! substitute the right-side values
 !
@@ -3252,6 +3254,7 @@ module interpolations
 ! local variables
 !
     integer      :: i, im1, ip1, im2, ip2
+    integer      :: iret
     real(kind=8) :: bl, bc, br, tt
     real(kind=8) :: wl, wc, wr, ww
     real(kind=8) :: df, lq, l3, zt
@@ -3587,7 +3590,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations for the left-side interpolation
 !
-    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n))
+    call tridiag(n, a(1:n,1), b(1:n,1), c(1:n,1), r(1:n,1), u(1:n), iret)
 
 ! substitute the left-side values
 !
@@ -3595,7 +3598,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations for the left-side interpolation
 !
-    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n))
+    call tridiag(n, a(1:n,2), b(1:n,2), c(1:n,2), r(1:n,2), u(1:n), iret)
 
 ! substitute the right-side values
 !
@@ -3914,7 +3917,7 @@ module interpolations
 
 ! local variables
 !
-    integer :: i
+    integer :: i, iret
 
 ! local arrays for derivatives
 !
@@ -3979,7 +3982,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations
 !
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
 
 ! apply the monotonicity preserving limiting
 !
@@ -4016,7 +4019,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations
 !
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
 
 ! apply the monotonicity preserving limiting
 !
@@ -4091,7 +4094,7 @@ module interpolations
 
 ! local variables
 !
-    integer :: i
+    integer :: i, iret
 
 ! local arrays for derivatives
 !
@@ -4157,7 +4160,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations
 !
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
 
 ! apply the monotonicity preserving limiting
 !
@@ -4194,7 +4197,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations
 !
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
 
 ! apply the monotonicity preserving limiting
 !
@@ -4262,7 +4265,7 @@ module interpolations
 
 ! local variables
 !
-    integer :: i
+    integer :: i, iret
 
 ! local arrays for derivatives
 !
@@ -4333,7 +4336,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations
 !
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
 
 ! apply the monotonicity preserving limiting
 !
@@ -4372,7 +4375,7 @@ module interpolations
 
 ! solve the tridiagonal system of equations
 !
-    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n))
+    call tridiag(n, a(1:n), b(1:n), c(1:n), r(1:n), u(1:n), iret)
 
 ! apply the monotonicity preserving limiting
 !