SUBROUTINE U_LU_BACKSUB(a,n,ndim,indx,b)
!***********************************************************************
!U_LU_BACKSUB solves the matrix equation a x = b, where a is in LU form as
!  generated by a prior call to U_LU_DECOMP
!References:
!  Flannery, Teukolsky, Vetterling, Numerical Recipes
!  W.A.Houlberg 1/99
!Input:
!  a-matrix in lu decomposed form
!  n-dimension of matrix
!  ndim-first dimension of a array
!  indx-vector showing row permutations due to partial pivoting
!  b-right hand side of equation (overwritten on return)
!Output:
!  b-solution vector
!Comments: 
!  The complete procedure is, given the equation a*x = b
!    CALL U_LU_DECOMP(a,n,ndim,indx,d,iflag)
!    IF(iflag.eq.0) CALL U_LU_BACKSUB(a,n,ndim,indx,b) 
!  and b is now the solution vector
!  To solve with a different right hand side, just reload b as  desired
!  and use the same LU decomposition, CALL U_LU_BACKSUB(a,n,ndim,indx,b)                                 *
!***********************************************************************
      IMPLICIT NONE
!Declaration of input variables
      INTEGER        indx(*),               n,
     #               ndim
      REAL           a(ndim,*),             b(*)
!Declaration of local variables
      INTEGER        i,                     ii,
     #               j,                     k
      REAL           sum
!Find the index of the first nonzero element of b
      ii=0
      DO i=1,n
        k=indx(i)
        sum=b(k)
        b(k)=b(i)
        IF(ii.ne.0) THEN
          DO j=ii,i-1
            sum=sum-a(i,j)*b(j)
          ENDDO   
        ELSEIF (sum.ne.0.0) THEN
          ii=i
        ENDIF
        b(i)=sum
      ENDDO    
!  Begin back substitution
      DO i=n,1,-1
        sum=b(i)
        IF(i.lt.n) THEN
          DO j=i+1,n
            sum=sum-a(i,j)*b(j)
          ENDDO   
        ENDIF
        b(i)=sum/a(i,i)
      ENDDO   
      RETURN
      END