C NCLFORTSTART
      SUBROUTINE my_sort(n,arr,arr_i,arr_j)
      INTEGER n
      REAL arr(n)
      INTEGER arr_i(n),arr_j(n)
C NCLEND
      INTEGER M,NSTACK
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp,b,c,a_i,a_j,temp_i,temp_j
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          b=arr_i(j)
          c=arr_j(j)
          do 11 i=j-1,l,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
            arr_i(i+1)=arr_i(i)
            arr_j(i+1)=arr_j(i)
11        continue
          i=l-1
2         arr(i+1)=a
          arr_i(i+1)=b
          arr_j(i+1)=c
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        temp_i=arr_i(k)
        temp_j=arr_j(k)
        arr(k)=arr(l+1)
        arr_i(k)=arr_i(l+1)
        arr_j(k)=arr_j(l+1)
        arr(l+1)=temp
        arr_i(l+1)=temp_i
        arr_j(l+1)=temp_j
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          temp_i=arr_i(l)
          temp_j=arr_j(l)
          arr(l)=arr(ir)
          arr_i(l)=arr_i(ir)
          arr_j(l)=arr_j(ir)
          arr(ir)=temp
          arr_i(ir)=temp_i
          arr_j(ir)=temp_j
        endif
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          temp_i=arr_i(l+1)
          temp_j=arr_j(l+1)
          arr(l+1)=arr(ir)
          arr_i(l+1)=arr_i(ir)
          arr_j(l+1)=arr_j(ir)
          arr(ir)=temp
          arr_i(ir)=temp_i
          arr_j(ir)=temp_j
        endif
        if(arr(l).gt.arr(l+1))then
          temp=arr(l)
          temp_i=arr_i(l)
          temp_j=arr_j(l)
          arr(l)=arr(l+1)
          arr_i(l)=arr_i(l+1)
          arr_j(l)=arr_j(l+1)
          arr(l+1)=temp
          arr_i(l+1)=temp_i
          arr_j(l+1)=temp_j
        endif
        i=l+1
        j=ir
        a=arr(l+1)
        a_i=arr_i(l+1)
        a_j=arr_j(l+1)
        a=arr(l+1)
        a_i=arr_i(l+1)
        a_j=arr_j(l+1)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        temp_i=arr_i(i)
        temp_j=arr_j(i)
        arr(i)=arr(j)
        arr_i(i)=arr_i(j)
        arr_j(i)=arr_j(j)
        arr(j)=temp
        arr_i(j)=temp_i
        arr_j(j)=temp_j
        goto 3
5       arr(l+1)=arr(j)
        arr_i(l+1)=arr_i(j)
        arr_j(l+1)=arr_j(j)
        arr(j)=a
        arr_i(j)=a_i
        arr_j(j)=a_j
        jstack=jstack+2
C       if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END