Hi, I am in trouble with using the routine hpsort in fortran 77.
What I want to do is to rearrange an array of integer in order.
hpsort works when the size of the array is small, but when I
increase the size (about 400), it does not work well.
When I call the routine hpsort several times, it still goes
successfully, but I believe the proper sorting routine rearrange
an array at one time.
I show the test code here. An integer array is generated
in random order at first. The output is the index n and the
value of the array lb(n), with some errors at n=382.
Would anyone find any hints of this problem?

! -------------------------------------------------------------------
program main

implicit double precision (a-h,o-z)
parameter (ndat=400)
dimension la(ndat),lb(ndat),lc(ndat)

idum=103
do n=1,ndat
la(n)=int(ran0(idum)*ndat)+1
enddo

do n=1,ndat
lb(n)=la(n)
enddo

call hpsort (ndat,lb)

write(6,*) ' n lb(n)'
do n=1,ndat
write(6,*) n,lb(n)
enddo

ld=lb(1)
do n=2,ndat
if(lb(n).lt.ld) then
write(6,*) 'Error at',n-1
endif
ld=lb(n)
enddo

stop
END

! -------------------------------------------------------------------
SUBROUTINE hpsort (n,ra)
! -------------------------------------------------------------------
integer n
integer ra(n)
! Sorts an array ra(1:n) into ascending numerical order using
! the Heapsort algorithm. n is input; ra is replaced on output
! by its sorted rearrangement.
integer i,ir,j,l
integer rra

if (n.lt.2) return
! The index l will be decremented from its initial value down
! to l during the "hiring" (heap creation) phase. Once it
! reaches l, the index ir will be decremented from its initial
! value down to l during the "retirement-and-promotion" (heap
! selection) phase.
l=n/2+1
ir=n
10 continue
if(l.gt.1) then
l=l-1
rra=ra(l)
else
rra=ra(ir)
ra(ir)=ra(1)
ir=ir-1
if(ir.eq.1) then
ra(1)=rra
return
endif
endif
i=l
j=l+1
20 if(j.le.ir) then
if(j.lt.ir) then
if(ra(j).lt.ra(j+1)) j=j+1
endif
if(rra.lt.ra(j)) then
ra(i)=ra(j)
i=j
j=j+j
else
j=ir+1
endif
goto 20
endif
ra(i)=rra
goto 10
END

! -------------------------------------------------------------------
function ran0 (idum)
! -------------------------------------------------------------------
implicit double precision (a-h,o-z)
parameter (ia=16807,im=2147483647,iq=127773,ir=2836)
parameter (mask=123459876)

am=1.d0/float(im)
idum=ieor(idum,mask)
k=idum/iq
idum=ia*(idum-k*iq)-ir*k
if(idum.lt.0) idum=idum+im
ran0=am*idum
idum=ieor(idum,mask)

return
END

There was a mistake in my code.
(It was "my" mistake, j=l+1 should be j=l+l in hpsort.)
I fixed it up and confirmed the routine works perfectly.
Sorry.