subroutine sort(n, ra, nloc)
      implicit double precision (o-z, a-h)
      dimension ra(n), nloc(n)
      if(n .eq. 1) return
      l = (n / 2) + 1
      ir = n
   10 continue
         if (l .gt. 1) then
            l = l - 1
            rra = ra(l)
            na = nloc(l)
         else
            rra = ra(ir)
            na = nloc(ir)
            ra(ir) = ra(1)
            nloc(ir) = nloc(1)
            ir = ir - 1
            if (ir .eq. 1) then
               ra(1) = rra
               nloc(1) = na
               return 
            end if
         end if
         i = l
         j = l + l
   20    if (j .le. ir) then
            if (j .lt. ir) then
               if (ra(j) .gt. ra(j + 1)) j = j + 1
            end if
            if (rra .gt. ra(j)) then
               ra(i) = ra(j)
               nloc(i) = nloc(j)
               i = j
               j = j + j
            else
               j = ir + 1
            end if
            goto 20
         end if
         ra(i) = rra
         nloc(i) = na
      goto 10
      end