!----------------------------------------------------------
SUBROUTINE LU(a, n, order)
       ! Factorise a square matrix a = LU
       ! Elements of (a) are replaced by elments of factors
       ! Pivoting is effected by row swaps. The row order
       ! is an output array                J.R.Dormand 1/96
!---------------------------------------------------------- 
INTEGER :: n, order(n), i, j, pivot, piv(1)
REAL (KIND =8) :: a(n, n), row(SIZE(order))
DO i = 1, n;  order(i) = i; END DO
DO i = 1, n-1
   piv = MAXLOC(ABS(a(i:n, i)))
   pivot = piv(1) + i - 1
   IF(pivot /= i) THEN
     row = a(i, :); a(i, :) = a(pivot, :); a(pivot, :) = row
     j = order(i); order(i) = order(pivot); order(pivot) = j
   END IF
   IF(ABS(a(i, i)) < 1.0D-10) THEN
      PRINT*, 'ZERO PIVOT'
      RETURN
   END IF
   DO j = i+1, n
      a(j, i) = a(j, i)/a(i, i)
      a(j, i+1:) = a(j, i+1:) - a(j, i)*a(i, i+1:)
   END DO
END DO
END SUBROUTINE LU
!----------------------------------------------------------
SUBROUTINE FBSubst(a, b, n, order)
 ! Carry out forward and backward substitution for LUx = b
 ! where (a) contains the elements of L and U. The pivoting
 ! order is supplied. The RHS vector (b) is replaced by the 
 ! solution vector.                        J.R.Dormand 1/96
!----------------------------------------------------------
INTEGER :: n, order(n), i
REAL (KIND = 8) :: a(n, n), b(n), x(SIZE(order))
DO i = 1, n
   x(i) = b(order(i))
END DO
DO i = 2, n
   x(i:n) = x(i:n) - a(i:n, i-1)*x(i-1)
END DO
DO i = n, 1, -1
   b(i) = x(i)/a(i, i)
   x(1:i-1) = x(1:i-1) - a(1:i-1, i)*b(i)
END DO
END SUBROUTINE FBSubst
!---------------------------------------------------------- 

