      subroutine d_backsubplus(matrix, invector, outvector, n, ier)
      implicit none
      integer n, ier
      double precision matrix(n,n), invector(n), outvector(n)
      integer i, k

      do i=1, n
         outvector(i)=invector(i)
         do k=1, i-1
            outvector(i) = outvector(i) - matrix(i,k) * outvector(k)
         end do
         if (matrix(i,i) .eq. 0.0) then
            write (0,*) "d_backsubplus:  matrix passed not invertible"
            ier=1
            return
         end if
         outvector(i) = outvector(i) / matrix(i,i)
      end do
      ier=0
      end ! subroutine d_backsubplus

c     ******************************************************************

      subroutine d_backsubminus(matrix, invector, outvector, n, ier)
      implicit none
      integer n, ier
      double precision matrix(n,n), invector(n), outvector(n)
      integer i, k

      do i=n, 1, -1
         outvector(i)=invector(i)
         do k=i+1, n
            outvector(i) = outvector(i) - matrix(i,k) * outvector(k)
         end do
         if (matrix(i,i) .eq. 0.0) then
            write (0,*) "d_backsubplus:  matrix passed not invertible"
            ier=1
            return
         end if
         outvector(i) = outvector(i) / matrix(i,i)
      end do
      ier=0
      end ! subroutine d_backsubminus
c     ******************************************************************

