subroutine INDEXX(N1,N2,M1,M2,ARRIN,INDX)
C-------------------------------------------------------------------------------
C  M1,M2  physical dimensions of arrays
C  N1,N2  logical dimensions of arrays (only this part of the array is indexed).
C-------------------------------------------------------------------------------
      integer*2 N1,N2,M1,M2,I,J,L,IR,INDXT,INDX(M1:M2)
      dimension ARRIN(M1:M2)
      do J=N1,N2
          INDX(J) = J
      end do
      L  = (N2-N1+1)/2+N1
      IR = N2
10    if (L .gt. N1) then
          L = L-1
          INDXT = INDX(L)
          Q = ARRIN(INDXT)
      else
          INDXT = INDX(IR)
          Q = ARRIN(INDXT)
          INDX(IR) = INDX(N1)
          IR = IR-1
          if (IR .eq. N1) then
              INDX(N1) = INDXT
              return
          end if
      end if
      I = L
      J = L+L-N1+1
      do while (J .le. IR)
          if (J .lt. IR .and. ARRIN(INDX(J)) .lt. ARRIN(INDX(J+1)))J=J+1
          if (Q .lt. ARRIN(INDX(J))) then
              INDX(I) = INDX(J)
              I = J
              J = J+J-N1+1
          else
              J = IR+1
          end if
      end do
      INDX(I) = INDXT
      go to 10
      end
C-------------------------------------------------------------------------------
      subroutine RANK(N1,N2,M1,M2,INDX,IRANK)
      dimension INDX(M1:M2),IRANK(M1:M2)
      do J=N1,N2
          IRANK(INDX(J)) = J
      end do
      return
      end
C-------------------------------------------------------------------------------
      subroutine SORT4(N1,N2,M1,M2,RA,RB,RC,RD)
C-------------------------------------------------------------------------------
C  Rearrange arrays RA,RB into ascending numerical order. Adapted from NUMERICAL
C  RECIPES, W.H. PRESS et al., page 231
C  M1,M2  physical dimensions of arrays
C  N1,N2  logical dimensions of arrays (only this part of the array is sorted).
C
C  RA,RB   double precision arrays!!!!!!
C-------------------------------------------------------------------------------
      integer*2 N1,N2,M1,M2,L,IR,I,J
      real*8    RRA,RRB,RA(M1:M2),RB(M1:M2),RC(M1:M2),RD(M1:M2)
      if (N1 .lt. M1 .or. N2 .gt. M2) stop 'SORT2: error in array dimensions'
      L = (N2-N1+1)/2+N1
      IR = N2
10    if (L .gt. N1) then
          L = L-1
          RRA = RA(L)
          RRB = RB(L)
          RRC = RC(L)
          RRD = RD(L)
      else
          RRA = RA(IR)
          RRB = RB(IR)
          RRC = RC(IR)
          RRD = RD(IR)
          RA(IR) = RA(N1)
          RB(IR) = RB(N1)
          RC(IR) = RC(N1)
          RD(IR) = RD(N1)
          IR = IR-1
          if (IR .eq. N1) then
              RA(N1) = RRA
              RB(N1) = RRB
              RC(N1) = RRC
              RD(N1) = RRD
              return
          end if
      end if
      I = L
      J = L+L-N1+1
      do while (J .le. IR)
          if (J .lt. IR .and. RA(J) .lt. RA(J+1)) J=J+1
          if (RRA .lt. RA(J)) then
              RA(I) = RA(J)
              RB(I) = RB(J)
              RC(I) = RC(J)
              RD(I) = RD(J)
              I = J
              J = J+J-N1+1
          else
              J = IR+1
          end if
      end do
      RA(I) = RRA
      RB(I) = RRB
      RC(I) = RRC
      RD(I) = RRD
      go to 10
      end
C-------------------------------------------------------------------------------
      subroutine SORT1(N1,N2,M1,M2,RA)
C-------------------------------------------------------------------------------
C  Rearrange array RA into ascending numerical order. Adapted from NUMERICAL
C  RECIPES, W.H. PRESS et al., page 231
C  N1,N2  logical dimensions of arrays (only this part of the array is sorted).
C  M1,M2  physical dimensions of arrays
C-------------------------------------------------------------------------------
      integer*2 N1,N2,M1,M2,L,IR,I,J
      real*4    RA(M1:M2)
      if (N1 .lt. M1 .or. N2 .gt. M2) stop 'SORT1: error in array dimensions'
      L = (N2-N1+1)/2+N1
      IR = N2
10    if (L .gt. N1) then
          L = L-1
          RRA = RA(L)
      else
          RRA = RA(IR)
          RA(IR) = RA(N1)
          IR = IR-1
          if (IR .eq. N1) then
              RA(N1) = RRA
              return
          end if
      end if
      I = L
      J = L+L-N1+1
      do while (J .le. IR)
          if (J .lt. IR .and. RA(J) .lt. RA(J+1)) J=J+1
          if (RRA .lt. RA(J)) then
              RA(I) = RA(J)
              I = J
              J = J+J-N1+1
          else
              J = IR+1
          end if
      end do
      RA(I) = RRA
      go to 10
      end