program sortcub
      implicit none
      integer*4 n0
      parameter (n0=100000)
      integer*4 ijk(n0,4),i,n,j,istack(n0),i1,i2,i3
      real*8 c(n0)
      character*80 fn
c
      write(6,*)'sorting anharmonic cubic constants'
      write(6,*)'filename:'
      read(5,*)fn
      open(8,file=fn)
      read(8,*,end=99,err=99)( ( (ijk(i,j),j=1,3) ,c(i) ) ,i=1,n0)
      write(6,*)'buffer overflew'
99    close(8)
      n=i-1
      write(6,*)n,' constants'
      do 1 i=1,n
      ijk(i,4)=1
      if(c(i).lt.0.0d0)ijk(i,4)=-1
1     c(i)=dabs(c(i))
      open(8,file='udia.lst')
      do 5 i=1,n
      i1=ijk(i,1)
      i2=ijk(i,2)
      i3=ijk(i,3)
5     if(i1.eq.i2.and.i2.eq.i3)
     1write(8,6000)(ijk(i,j),j=1,3),c(i)*dble(ijk(i,4))
      close(8)
      write(6,*)'unsorted diagonals written into udia.lst'
      call sort2(n0,n,c,ijk,istack)
      write(6,*)'sorted'
      open(8,file='cdia.lst')
      do 2 i=1,n
      i1=ijk(i,1)
      i2=ijk(i,2)
      i3=ijk(i,3)
2     if(i1.eq.i2.and.i2.eq.i3)
     1write(8,6000)(ijk(i,j),j=1,3),c(i)*dble(ijk(i,4))
6000  format(3i3,f12.4)
      close(8)
      write(6,*)'diagonals written into cdia.lst'
      open(8,file='c2.lst')
      do 3 i=1,n
      i1=ijk(i,1)
      i2=ijk(i,2)
      i3=ijk(i,3)
3     if((i1.eq.i2.and.i1.ne.i3).or.
     1   (i1.eq.i3.and.i1.ne.i2).or.
     2   (i2.eq.i3.and.i1.ne.i2))
     3write(8,6000)(ijk(i,j),j=1,3),c(i)*dble(ijk(i,4))
      close(8)
      write(6,*)'2-center written into c2.lst'
      open(8,file='c3.lst')
      do 4 i=1,n
      i1=ijk(i,1)
      i2=ijk(i,2)
      i3=ijk(i,3)
4     if(i1.ne.i2.and.i1.ne.i3.and.i2.ne.i3)
     1write(8,6000)(ijk(i,j),j=1,3),c(i)*dble(ijk(i,4))
      close(8)
      write(6,*)'3-center written into c3.lst'
      stop
      end
c     =========================================================
      subroutine sort2(n0,n,arr,brr,istack)
c     Sorts an array arr(1:n) into ascending order while
c     making rearangement of the array brr(1:n)
c     M ... when buble sort is switched on
c
      implicit none
      integer*4 n,M,n0
      parameter (M=10)
      real*8 arr(*),a,temp
      integer*4 brr(n0,3),b(3),itemp
      integer*4 i,ir,j,jstack,k,l,istack(n0),i3
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
c      insertion sort for small Ms:
       do j=l+1,ir
        a=arr(j)
        do i3=1,3
         b(i3)=brr(j,i3)
        enddo
        do i=j-1,l,-1
         if(arr(i).le.a)goto 2
         arr(i+1)=arr(i)
         do i3=1,3
          brr(i+1,i3)=brr(i,i3)
         enddo
        enddo
        i=l-1
2       arr(i+1)=a
        do i3=1,3
         brr(i+1,i3)=b(i3)
        enddo
       enddo
       if(jstack.eq.0)return
c      pop stack and begind new round of partitioning
       ir=istack(jstack)
       l=istack(jstack-1)
       jstack=jstack-2
      else
c      choose median of left, center and right elements as partitioning
c      element a. Also rearrange so that a(l)<=a(l+1)<=a(ir)
       k=(l+ir)/2
       temp=arr(k)
       arr(k)=arr(l+1)
       arr(l+1)=temp
       do i3=1,3
        itemp=brr(k,i3)
        brr(k,i3)=brr(l+1,i3)
        brr(l+1,i3)=itemp
       enddo
       if(arr(l).gt.arr(ir))then
        temp=arr(l)
        arr(l)=arr(ir)
        arr(ir)=temp
        do i3=1,3
         itemp=brr(l,i3)
         brr(l,i3)=brr(ir,i3)
         brr(ir,i3)=itemp
        enddo
       endif
       if(arr(l+1).gt.arr(ir))then
        temp=arr(l+1)
        arr(l+1)=arr(ir)
        arr(ir)=temp
        do i3=1,3
         itemp=brr(l+1,i3)
         brr(l+1,i3)=brr(ir,i3)
         brr(ir,i3)=itemp
        enddo
       endif
       if(arr(l).gt.arr(l+1))then
        temp=arr(l)
        arr(l)=arr(l+1)
        arr(l+1)=temp
        do i3=1,3
         itemp=brr(l,i3)
         brr(l,i3)=brr(l+1,i3)
         brr(l+1,i3)=itemp
        enddo
       endif
c      initialize pointers for partitioning
       i=l+1
       j=ir
c      partitioning element
       a=arr(l+1)
       do i3=1,3
        b(i3)=brr(l+1,i3)
       enddo
c      beginninf of innermost loop
3      continue
c      scan up to find element >a
       i=i+1
       if(arr(i).lt.a)goto 3
4      continue
c      scan down to find element <a
       j=j-1
       if(arr(j).gt.a)goto 4
c      pointers crossed; exit whe partitioning complete
       if(j.lt.i)goto 5
c      exchange elements of both arrays
       temp=arr(i)
       arr(i)=arr(j)
       arr(j)=temp
       do i3=1,3
        itemp=brr(i,i3)
        brr(i,i3)=brr(j,i3)
        brr(j,i3)=itemp
       enddo
       goto 3
c      end of innermost loop
c      insert partitioning elements in both arrays
5      arr(l+1)=arr(j)
       arr(j)=a
       do i3=1,3
        brr(l+1,i3)=brr(j,i3)
        brr(j,i3)=b(i3)
       enddo
       jstack=jstack+2
c      push pointers to larger subarray on stack,
c      process smaller subarray immediately
       if(jstack.gt.n0)pause 'nstack too small'
       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