subroutine tarjan( n, nz, num_block, col_ind, row_ptr, row_len,
     1	    stack, block_start, block_end, stack_pos, prev_node,
     2	    lowest, jused, visited )
	implicit none
c***********************************************************************
c	subroutine tarjan
c
c	Apply Tarjan's algorithm
c***********************************************************************
c	Parameters
c
c	n - integer: order of the matrix
c	nz - integer: number of nonzero elements in the matrix
c	num_block - integer: upon return is the number of strongly
c	    connected components
c	col_ind - integer(nz): column index for each nonzero element
c	row_ptr - integer(n): index of the first element of each row
c	row_len - integer(n): length of each row
c	stack - integer(n): used as the stack for Tarjan's algorithm.
c	block_start - integer(n): upon return it will contain
c	    the index to the first row of each block in the stack
c	block_end - integer(n): work space used for the top element
c	    in each block, (is numbered backwards)
c	stack_pos - integer(n): work space used for the position of
c	    each row in the stack. Upon return it will contain the 
c           rows in each block
c	prev_node - integer(n): work space for the node previous
c	    to each node
c	lowest - integer(n): work space for the lowest node reachable
c	    for each node
c	jused - integer(n): work space for the number of edges already
c	    traced from the node
c	visited - integer(n): work space indicating if a node has been
c	    visited:
c		 0 - node has not been visited
c		 1 - node has been visited
c***********************************************************************
	integer n, nz, num_block
	integer col_ind(nz), row_ptr(n), row_len(n)
	integer stack(n), block_start(n), block_end(n)
	integer stack_pos(n), prev_node(n), lowest(n)
	integer jused(n), visited(n)
c***********************************************************************
c	Local variables
c***********************************************************************
	integer n1, i, j, jstart, jstop, jrow
	integer next_save, next_stack, cur_node, last_node
D	character ifmt*10
D	ifmt = '(16i5)'
c***********************************************************************
c	Initialize the variables
c***********************************************************************
	n1 = n + 1
	num_block = 0
	next_save = n 
	next_stack = 1
	do i = 1, n
	    visited(i) = 0
	    jused(i) = 0
	    lowest(i) = n1
	end do

c	Loop through all the nodes

	do i = 1, n
D	    print *,'loop i', i

c	    If a node has already been visited, then go to the next node

	    if ( visited(i) .ne. 0 ) goto 1000

D	    print *,'adding i',i
	    cur_node = i
	    prev_node(cur_node) = -1

100	    continue

c	    Place the current node on the stack

D	    print *,'placing i at st', cur_node, next_stack
	    visited(cur_node) = 1
	    lowest(cur_node) = next_stack
	    stack(next_stack) = cur_node
	    stack_pos(cur_node) = next_stack
	    next_stack = next_stack + 1

200	    continue

c	    Loop through the edges for the current node until a new
c	    node is found or all the edges from the current node have
c	    been checked

	    jstart = row_ptr(cur_node) + jused(cur_node)
	    jstop = row_ptr(cur_node) + row_len(cur_node) - 1
D	    print *,'looping for ',cur_node,jstart,jstop
	    do j = jstart, jstop
		jrow = col_ind(j)
D		print *,'index is',jrow,visited(jrow)

c		If visited <> 0 then the node has already been visited,
c		so just update the pointer to the lowest node

		if ( visited(jrow) .ne. 0 ) then
D		    print*,'lowest cur, jrow',lowest(cur_node),lowest(jrow)
		    if (lowest(jrow).lt.lowest(cur_node))then
			lowest(cur_node) = lowest(jrow)
		    end if
		else

c		    This a a new node to be added to the stack

		    prev_node(jrow) = cur_node
		    jused(cur_node) = jused(cur_node) + j-jstart+1
D		    print *,'adding jrow',jrow,jused(cur_node)
		    cur_node = jrow
		    goto 100
		end if
	    end do

c	    All the edges for the current node have been checked

	    jused(cur_node) = row_len(cur_node)

D	    print *,'through edges',cur_node,stack_pos(cur_node),
D     1	        lowest(cur_node)

c	    If the lowest node that can be reached is the current
c	    node, then this forms a block.  The block is all nodes
c	    above the current node on the stack (including the
c	    current node).

	    if ( lowest(cur_node) .eq. stack_pos(cur_node) ) then
		num_block = num_block + 1
D		print *,'new block',num_block,next_save,next_stack
		block_end(num_block) = next_save
		do j = next_stack-1, stack_pos(cur_node), -1
		    stack(next_save) = stack(j)
		    lowest(stack(j)) = n+1
		    next_save = next_save - 1
		end do
		next_stack = stack_pos(cur_node)
	    end if

c	    Backtrack to the previous node

	    last_node = cur_node
	    cur_node = prev_node(last_node)
	    if ( cur_node .gt. 0 ) then
		if (lowest(last_node).lt.lowest(cur_node)) then
		    lowest(cur_node) = lowest(last_node)
		end if
		goto 200
	    end if

1000	    continue
	end do

c	Now set up the the pointers to the beginning of the blocks

c	block_start(1) = 1
c	do i = 2, num_block
c	    block_start(i) = block_end(num_block-i+2) + 1
c	end do
	block_start(1) = 1
	do i = 2, num_block
	    block_start(i) = n - block_end(i) + 1
	end do
	block_start(num_block+1) = n+1
c	reverse the stack.
	do i = 1,n
	    stack_pos(i) = stack(n-i+1)
	enddo
	return
	end