c     QCKSRT (a version that sorts a column of integers and does the
c             same index shift for a second column of real numbers)
c
c     Reference: Press, et. al.
c
c     Ellen McGrattan, 4-5-95
c
      subroutine qcksrt3(n,arr,x)
      implicit real*8 (a-h,o-z)
      parameter(m=7,nstack=100,fm=7875.,fa=211.,fc=1663.,fmi=1./fm)
      double precision x(n)
      integer arr(n),istack(nstack)
      jstack = 0
      l      = 1
      ir     = n
      fx     = 0.
 10   if (ir-l.lt.m) then
        do 13 j=l+1,ir
          a  = arr(j)
          p  = x(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a) go to 12
            arr(i+1) = arr(i)
            x(i+1)   = x(i)
 11       continue
          i        = 0
 12       arr(i+1) = a
          x(i+1)   = p 
 13     continue
        if (jstack.eq.0) then
          do 55 k=1,n-1
            if(arr(k).gt.arr(k+1)) write(*,*) 'Qcksrt failed'
 55       continue
          return
        endif
        ir     = istack(jstack)
        l      = istack(jstack-1)
        jstack = jstack - 2
      else
        i      = l
        j      = ir
        fx     = mod(fx*fa+fc,fm)
        iq     = l+(ir-l+1)*(fx*fmi)
        a      = arr(iq)
        p      = x(iq)
        arr(iq)= arr(l)
        x(iq)  = x(l)
 20     continue
 21     if (j.gt.0) then
          if (a.lt.arr(j)) then
            j = j-1
            go to 21
          endif
        endif
        if (j.le.i) then
          arr(i) = a
          x(i)   = p 
          go to 30
        endif
        arr(i) = arr(j)
        x(i)   = x(j)
        i      = i+1
 22     if (i.le.n) then
          if (a.gt.arr(i)) then
            i = i+1
            go to 22
          endif
        endif
        if (j.le.i) then
          arr(j) = a
          x(j)   = p 
          i      = j
          go to 30
        endif
        arr(j) = arr(i)
        x(j)   = x(i)
        j      = j-1
        go to 20
 30     jstack = jstack +2
        if (jstack.gt.nstack) write(*,*) 'Increase size of NSTACK'
        if (ir-i.ge.i-l) then
          istack(jstack)   = ir
          istack(jstack-1) = i+1
          ir               = i-1
        else
          istack(jstack)   = i-1
          istack(jstack-1) = l
          l                = i+1
        endif
      endif
      go to 10
      end