************************************************************
*
*                                               Robert Renka
*                                       Oak Ridge Natl. Lab.
*
*   This subroutine uses an order N*LOG(N) quick sort to sort an integer
* array X into increasing order.  The algorithm is as follows.  IND is
* initialized to the ordered sequence of indices 1,...,N, and all
* interchanges are applied to IND.  X is divided into two portions by
* picking a central element T.  The first and last elements are compared
* with T, and interchanges are applied as necessary so that the three
* values are in ascending order.  Interchanges are then applied to that
* all elements greater than T are in the upper portion of the array and
* all elements less than T are in the lower portion.  The upper and lower
* indices of one of the portions are saved in local arrays, and the
* process is repeated iteratively on the other portion.  When a portion
* is completely sorted, the process begins again by retrieving the 
* indices bounding another unsorted portion.
*
* INPUT PARAMETERS -   N - Length of the array X.
*
*                      X - Vector of length N to be sorted.
*
*                    IND - Vector of length .GE. N.
*
* N and X are not altered by this routine.
*
* OUTPUT PARAMETER - IND - Sequence of indices 1,...,N
*                          permuted in the same fashion as X
*                          would be.  Thus, the ordering on
*                          X is defined by Y(I) = X(IND(I)).
*
* INTRINSIC FUNCTIONS CALLED BY QSORTI - IFIX, FLOAT
*
************************************************************
*
* NOTE -- IU and IL must be dimensioned .GE. LOG(N) where LOG has base 2.
*
************************************************************
      SUBROUTINE QSORTI (X, IND, N)
      INTEGER N, X(N), IND(N)

      INTEGER IU(32), IL(32)
      INTEGER M, I, J, K, L, IJ, IT, ITT, INDX, T
      REAL    R

* LOCAL PARAMETERS -
*
* IU,IL =  temporary storage for the upper and lower
*            indices of portions of the array
* M =      index for IU and IL
* I,J =    lower and upper indices of a portion of X
* K,L =    indices in the range I,...,J
* IJ =     randomly chosen index between I and J
* IT,ITT = temporary storage for interchanges in IND
* INDX =   temporary index for X
* R =      pseudo random number for generating IJ
* T =      central element of X

      IF (N .LE. 0) RETURN

* Initialize IND, M, I, J, and R
      DO 1 I = 1,N
  1     IND(I) = I
      M = 1
      I = 1
      J = N
      R = .375

* top of loop
  2   IF (I .GE. J) GO TO 10
      IF (R > .5898437) GO TO 3
      R = R + .0390625
      GO TO 4
  3   R = R - .21875

* initialize K
  4   K = I

* select a central element of X and save it in T
      IJ = I + IFIX(R*FLOAT(J-I))
      IT = IND(IJ)
      T = X(IT)

* If the first element of the array is greater than T,
*   interchange it with T
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 5
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)

* initialize L
  5   L = J

* if the last element of the array is less than T,
*   interchange it with T
      INDX = IND(J)
      IF (X(INDX) .GE. T) GO TO 7
      IND(IJ) = INDX
      IND(J) = IT
      IT = INDX
      T = X(IT)

* if the first element of the array is greater than T,
*   interchange it with T
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 7
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)
      GO TO 7

* interchange elements K and L
  6   ITT = IND(L)
      IND(L) = IND(K)
      IND(K) = ITT

* Find an element in the upper part of the array which is
*   not larger than T
  7   L = L - 1
      INDX = IND(L)
      IF (X(INDX) > T) GO TO 7

* Find an element in the lower part of the array which is 
*   not smaller than T
  8   K = K + 1
      INDX = IND(K)
      IF (X(INDX) < T) GO TO 8

* if K .LE. L, interchange elements K and L
      IF (K .LE. L) GO TO 6

* save the upper and lower subscripts of the portion of the
*   array yet to be sorted
      IF (L-I .LE. J-K) GO TO 9
      IL(M) = I
      IU(M) = L
      I = K
      M = M + 1
      GO TO 11

  9   IL(M) = K
      IU(M) = J
      J = L
      M = M + 1
      GO TO 11

* begin again on another unsorted portion of the array
 10   M = M - 1
      IF (M .EQ. 0) RETURN
      I = IL(M)
      J = IU(M)

 11   IF (J-I .GE. 11) GO TO 4
      IF (I .EQ. 1) GO TO 2
      I = I - 1

* sort elements I+1,...,J.  Note that 1 .LE. I < J and J-I < 11.
 12   I = I + 1
      IF (I .EQ. J) GO TO 10
      INDX = IND(I+1)
      T = X(INDX)
      IT = INDX
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 12
      K = I

   13 IND(K+1) = IND(K)
      K = K - 1
      INDX = IND(K)
      IF (T < X(INDX)) GO TO 13
      IND(K+1) = IT
      GO TO 12
      END ! of QSORTI


************************************************************************
      SUBROUTINE QSORTR (X, IND, N)
      INTEGER N, IND(N)
      REAL    X(N)

      INTEGER IU(32), IL(32)
      INTEGER M, I, J, K, L, IJ, IT, ITT, INDX
      REAL    R, T

* LOCAL PARAMETERS -
*
* IU,IL =  temporary storage for the upper and lower
*            indices of portions of the array
* M =      index for IU and IL
* I,J =    lower and upper indices of a portion of X
* K,L =    indices in the range I,...,J
* IJ =     randomly chosen index between I and J
* IT,ITT = temporary storage for interchanges in IND
* INDX =   temporary index for X
* R =      pseudo random number for generating IJ
* T =      central element of X
************************************************************************
      IF (N .LE. 0) RETURN

* initialize IND, M, I, J, and R
      DO 1 I = 1, N
  1     IND(I) = I
      M = 1
      I = 1
      J = N
      R = .375

* top of loop
  2   IF (I .GE. J) GO TO 10
      IF (R > .5898437) GO TO 3
      R = R + .0390625
      GO TO 4
  3   R = R - .21875

* initialize K
  4   K = I

* select a central element of X and save it in T
      IJ = I + IFIX(R*FLOAT(J-I))
      IT = IND(IJ)
      T = X(IT)

* if the first element of the array is greater than T, interchange it with T
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 5
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)

* initialize L
  5   L = J

* if the last element of the array is less than T, interchange it with T
      INDX = IND(J)
      IF (X(INDX) .GE. T) GO TO 7
      IND(IJ) = INDX
      IND(J) = IT
      IT = INDX
      T = X(IT)

* if the first element of the array is greater than T, interchange it with T
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 7
      IND(IJ) = INDX
      IND(I) = IT
      IT = INDX
      T = X(IT)
      GO TO 7

* interchange elements K and L
  6   ITT = IND(L)
      IND(L) = IND(K)
      IND(K) = ITT

* find an element in the upper part of the array which is not larger than T
  7   L = L - 1
      INDX = IND(L)
      IF (X(INDX) > T) GO TO 7

* find an element in the lower part of the array which is not smaller than T
  8   K = K + 1
      INDX = IND(K)
      IF (X(INDX) < T) GO TO 8

* if K .LE. L, interchange elements K and L
      IF (K .LE. L) GO TO 6

* save the upper and lower subscripts of the portion of the 
*   array yet to be sorted
      IF (L-I .LE. J-K) GO TO 9
      IL(M) = I
      IU(M) = L
      I = K
      M = M + 1
      GO TO 11

  9   IL(M) = K
      IU(M) = J
      J = L
      M = M + 1
      GO TO 11

* begin again on another unsorted portion of the array
  10  M = M - 1
      IF (M .EQ. 0) RETURN
      I = IL(M)
      J = IU(M)

  11  IF (J-I .GE. 11) GO TO 4
      IF (I .EQ. 1) GO TO 2
      I = I - 1

* sort elements I+1,...,J.  Note that 1 .LE. I < J and J-I < 11.
  12  I = I + 1
      IF (I .EQ. J) GO TO 10
      INDX = IND(I+1)
      T = X(INDX)
      IT = INDX
      INDX = IND(I)
      IF (X(INDX) .LE. T) GO TO 12
      K = I

  13  IND(K+1) = IND(K)
      K = K - 1
      INDX = IND(K)
      IF (T < X(INDX)) GO TO 13
      IND(K+1) = IT
      GO TO 12
      END ! of QSORTR


************************************************************
*
*                                               Robert Renka
*                                       Oak Ridge Natl. Lab.
*
*   This routine applies a set of permutations to a vector.
*
* INPUT PARAMETERS -  N - Length of A and IP.
*
*                    IP - Vector containing the sequence of 
*                        Integers 1,...,N permuted in the 
*                         same fashion that A is to be permuted
*
*                     A - Vector to be permuted.
*
* N and IP are not altered by this routine.
*
* OUTPUT PARAMETER -  A - Reordered vector reflecting the 
*                         permutations defined by IP.
*
************************************************************
      SUBROUTINE IORDER (A, IP, N)
      INTEGER N, A(N), IP(N)

      INTEGER NN, K, J, IPJ, TEMP

* LOCAL PARAMETERS -
*
* NN =   Local copy of N
* K =    Index for IP and for the first element of A[] in a permutation
* J =    Index for IP and A, J .GE. K
* IPJ =  IP[J]
* TEMP = Temporary storage for A[K]
      NN = N
      IF (NN < 2) RETURN
      K = 1

* Loop on permutations
  1   J = K
      TEMP = A(K)

* Apply permutation to A.  IP[J] is marked (made negative)
*   as being included in the permutation.
  2   IPJ = IP(J)
      IP(J) = -IPJ
      IF (IPJ .EQ. K) GO TO 3
      A(J) = A(IPJ)
      J = IPJ
      GO TO 2
  3   A(J) = TEMP

* Search for an unmarked element of IP
  4   K = K + 1
      IF (K > NN) GO TO 5
      IF (IP(K) > 0) GO TO 1
      GO TO 4

* All permutations have been applied.  Unmark IP.
  5   DO 6 K = 1,NN
  6     IP(K) = -IP(K)
      RETURN
      END ! of IORDER


************************************************************
      SUBROUTINE RORDER (A, IP, N)
      INTEGER N, IP(N)
      REAL    A(N)

      INTEGER NN, K, J, IPJ
      REAL    TEMP

* LOCAL PARAMETERS -
* NN =   Local copy of N
* K =    Index for IP and for the first element of A in a permutation
* J =    Index for IP and A, J .GE. K
* IPJ =  IP(J)
* TEMP = Temporary storage for A(K)
************************************************************************
      NN = N
      IF (NN < 2) RETURN
      K = 1

* Loop on permutations
    1 J = K
      TEMP = A(K)

* Apply permutation to A.  IP(J) is marked (made negative)
*   as being included in the permutation.
    2 IPJ = IP(J)
      IP(J) = -IPJ
      IF (IPJ .EQ. K) GO TO 3
      A(J) = A(IPJ)
      J = IPJ
      GO TO 2
    3 A(J) = TEMP

* Search for an unmarked element of IP
    4 K = K + 1
      IF (K > NN) GO TO 5
      IF (IP(K) > 0) GO TO 1
      GO TO 4

* All permutations have been applied.  Unmark IP.
    5 DO 6 K = 1,NN
    6   IP(K) = -IP(K)
      RETURN
      END


************************************************************************
      SUBROUTINE DORDER (A, IP, N)
      INTEGER N, IP(N)
      DOUBLE PRECISION A(N)

      INTEGER NN, K, J, IPJ
      DOUBLE PRECISION TEMP

      NN = N
      IF (NN .LT. 2) RETURN
      K = 1

* Loop on permutations
    1 J = K
      TEMP = A(K)

* Apply permutation to A.  IP(J) is marked (made negative)
*   as being included in the permutation.
    2 IPJ = IP(J)
      IP(J) = -IPJ
      IF (IPJ .EQ. K) GO TO 3
      A(J) = A(IPJ)
      J = IPJ
      GO TO 2
    3 A(J) = TEMP

* Search for an unmarked element of IP
    4 K = K + 1
      IF (K .GT. NN) GO TO 5
      IF (IP(K) .GT. 0) GO TO 1
      GO TO 4

* All permutations have been applied.  Unmark IP.
    5 DO 6 K = 1,NN
    6   IP(K) = -IP(K)
      RETURN
      END


*-----------------------------------------------------------------------
* ISORT - sort integers
*-----------------------------------------------------------------------
*          Sort an array and optionally make the same interchanges in
*            an auxiliary array.  The array may be sorted in increasing
*            or decreasing order.  A slightly modified QUICKSORT
*            algorithm is used.  From the SLATEC library, public domain.
*
*      by   Jones, R. E., (SNLA)
*           Kahaner, D. K., (NBS)
*           Wisniewski, J. A., (SNLA)
*
*   Description of Parameters
*      IX - integer array of values to be sorted
*      N  - number of values in integer array IX to be sorted
*
*  References:  R. C. Singleton, Algorithm 347, An efficient algorithm
*                 for sorting with minimal storage, Communications of
*                 the ACM, 12, 3 (1969), pp. 185-187.
*
*  REVISION HISTORY  (YYMMDD)
*   761118  Date written
*   810801  Modified by David K. Kahaner.
*   890531  Changed all specific intrinsics to generic.  (WRB)
*   890831  Modified array declarations.  (WRB)
*   891009  Removed unreferenced statement labels.  (WRB)
*   891009  REVISION DATE from Version 3.2
*   891214  Prologue converted to Version 4.0 format.  (BAB)
*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
*   901012  Declared all variables; changed X,Y to IX,IY. (M. McClain)
*   920501  Reformatted the REFERENCES section.  (DWL, WRB)
*   920519  Clarified error messages.  (DWL)
*   920801  Declarations section rebuilt and code restructured to use
*           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
*   111005  Removed call to XERMSG, simplified. (AJA)
*-----------------------------------------------------------------------
      SUBROUTINE ISORT (IX, N)

*     .. Scalar Arguments ..
      INTEGER N
*     .. Array Arguments ..
      INTEGER IX(N)
*     .. Local Scalars ..
      REAL R
      INTEGER I, IJ, J, K, L, M, NN, T, TT
*     .. Local Arrays ..
      INTEGER IL(32), IU(32)

*        Begin
      NN = N
      IF (NN < 1) THEN
         PRINT *, 'The number of values to be sorted is not positive.'
         RETURN
      END IF

*     Sort IX only
      M = 1
      I = 1
      J = NN
      R = 0.375E0

   20 IF (I .EQ. J) GO TO 60
      IF (R .LE. 0.5898437E0) THEN
         R = R+3.90625E-2
      ELSE
         R = R-0.21875E0
      ENDIF

   30 K = I

*     Select a central element of the array and save it in location T
      IJ = I + INT(FLOAT(J-I)*R)
      T = IX(IJ)

*     If first element of array is greater than T, interchange with T
      IF (IX(I) > T) THEN
         IX(IJ) = IX(I)
         IX(I) = T
         T = IX(IJ)
      ENDIF
      L = J

*     If last element of array is less than than T, interchange with T
      IF (IX(J) < T) THEN
         IX(IJ) = IX(J)
         IX(J) = T
         T = IX(IJ)

*        If first element of array is greater than T, interchange with T
         IF (IX(I) > T) THEN
            IX(IJ) = IX(I)
            IX(I) = T
            T = IX(IJ)
         ENDIF
      ENDIF

*     Find an element in the second half of the array which is smaller than T
   40 L = L-1
      IF (IX(L) > T) GO TO 40

*     Find an element in the first half of the array which is greater than T
   50 K = K+1
      IF (IX(K) < T) GO TO 50

*     Interchange these elements
      IF (K .LE. L) THEN
         TT = IX(L)
         IX(L) = IX(K)
         IX(K) = TT
         GO TO 40
      ENDIF

*     Save upper and lower subscripts of the array yet to be sorted
      IF (L-I > J-K) THEN
         IL(M) = I
         IU(M) = L
         I = K
         M = M+1
      ELSE
         IL(M) = K
         IU(M) = J
         J = L
         M = M+1
      ENDIF
      GO TO 70

*     Begin again on another portion of the unsorted array
   60 M = M-1
      IF (M .EQ. 0) GO TO 200
      I = IL(M)
      J = IU(M)

   70 IF (J-I .GE. 1) GO TO 30
      IF (I .EQ. 1) GO TO 20
      I = I-1

   80 I = I+1
      IF (I .EQ. J) GO TO 60
      T = IX(I+1)
      IF (IX(I) .LE. T) GO TO 80
      K = I

   90 IX(K+1) = IX(K)
      K = K-1
      IF (T < IX(K)) GO TO 90
      IX(K+1) = T
      GO TO 80

 200  END ! of ISORT


*-----------------------------------------------------------------------
* SSORT - sort single precision numbers
*-----------------------------------------------------------------------
*   Sort an array in increasing order.  A slightly modified QUICKSORT 
* algorithm is used.  From the SLATEC library, public domain. 
*
*      By  R.E. Jones, (SNLA) 
*      and J. A. Wisniewski, (SNLA)
*
*   Description of arguments
*      X[]     array of values to be sorted
*      N       number of values in array X to be sorted
*
*  REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
*                 for sorting with minimal storage, Communications of
*                 the ACM, 12, 3 (1969), pp. 185-187.
*
*  REVISION HISTORY  (YYMMDD)
*   761101  Date written
*   761118  Modified to use the Singleton quicksort algorithm.  (JAW)
*   890531  Changed all specific intrinsics to generic.  (WRB)
*   890831  Modified array declarations.  (WRB)
*   891009  Removed unreferenced statement labels.  (WRB)
*   891024  Changed category.  (WRB)
*   891024  REVISION DATE from Version 3.2
*   891214  Prologue converted to Version 4.0 format.  (BAB)
*   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
*   901012  Declared all variables; changed X,Y to SX,SY. (M. McClain)
*   920501  Reformatted the REFERENCES section.  (DWL, WRB)
*   920519  Clarified error messages.  (DWL)
*   920801  Declarations section rebuilt and code restructured to use
*           IF-THEN-ELSE-ENDIF.  (RWC, WRB)
*   110517  removed call to XERMSG, simplified (AJA)
*-----------------------------------------------------------------------
      SUBROUTINE SSORT (X, N)

*     .. Scalar Arguments ..
      INTEGER  N
*     .. Array Arguments ..
      REAL X(N)
*     .. Local Scalars ..
      REAL R, T, TT
      INTEGER I, IJ, J, K, L, M, NN
*     .. Local Arrays ..
      INTEGER IL(32), IU(32)
*-----------------------------------------------------------------------
      NN = N
      IF (NN < 1) THEN
          PRINT *, 'The number of values to be sorted is not positive.'
         RETURN
      ENDIF

*     Sort X only
      M = 1
      I = 1
      J = NN
      R = 0.375E0

   20 IF (I .EQ. J) GO TO 60
      IF (R .LE. 0.5898437E0) THEN
         R = R+3.90625E-2
      ELSE
         R = R-0.21875E0
      ENDIF

   30 K = I

*     Select a central element of the array and save it in location T
      IJ = I + INT(FLOAT(J-I)*R)
      T = X(IJ)

*     If first element of array is greater than T, interchange with T
      IF (X(I) > T) THEN
         X(IJ) = X(I)
         X(I) = T
         T = X(IJ)
      ENDIF
      L = J

*     If last element of array is less than than T, interchange with T
      IF (X(J) < T) THEN
         X(IJ) = X(J)
         X(J) = T
         T = X(IJ)

*        If first element of array is greater than T, interchange with T
         IF (X(I) > T) THEN
            X(IJ) = X(I)
            X(I) = T
            T = X(IJ)
         ENDIF
      ENDIF

*     Find an element in the second half of the array which is smaller than T
   40 L = L-1
      IF (X(L) > T) GO TO 40

*     Find an element in the first half of the array which is greater than T
   50 K = K+1
      IF (X(K) < T) GO TO 50

*     Interchange these elements
      IF (K .LE. L) THEN
         TT = X(L)
         X(L) = X(K)
         X(K) = TT
         GO TO 40
      ENDIF

*     Save upper and lower subscripts of the array yet to be sorted
      IF (L-I > J-K) THEN
         IL(M) = I
         IU(M) = L
         I = K
         M = M+1
      ELSE
         IL(M) = K
         IU(M) = J
         J = L
         M = M+1
      ENDIF
      GO TO 70

*     Begin again on another portion of the unsorted array
   60 M = M-1
      IF (M .EQ. 0) GO TO 200
      I = IL(M)
      J = IU(M)

   70 IF (J-I .GE. 1) GO TO 30
      IF (I .EQ. 1) GO TO 20
      I = I-1

   80 I = I+1
      IF (I .EQ. J) GO TO 60
      T = X(I+1)
      IF (X(I) .LE. T) GO TO 80
      K = I

   90 X(K+1) = X(K)
      K = K-1
      IF (T < X(K)) GO TO 90
      X(K+1) = T
      GO TO 80

  200 END ! of SSORT


*-----------------------------------------------------------------------
* rsorti - radix sort and carry along an integer array
* by Andy Allinger, 2011-2016, released to the public domain
* This program may be used by any person for any purpose
*
* origin:  Herman Hollerith, 1887
*
* NOTE:  This routine is a sort-and-carry.  To obtain a permutation 
* vector to apply to other arrays, you must first initialize IND to 
* 1, 2, 3...N
*
* NOTE ALSO:
*   if B is positive, then the integers will be interpreted as UNSIGNED
*   note that FORTRAN does not support unsigned integers, so make sure
*   you know what you are doing if you use B=32
*   if B is negative, the sign bit will be inspected as a final step
*   So, for the range -256...255 set B = -8
*
*__Name____Type_____In/Out_________Description__________________________
*  D[N]    Integer  Both           Array, returned sorted
*  IND[N]  Integer  Both           Array to carry along, meaning
*                                         permuted in the same fashion as D
*  N       Integer  In             Length of the array
*  B       Integer  In             Number of significant bits
*  W[2@N]  Integer  Neither        Workspace
*-----------------------------------------------------------------------
      SUBROUTINE RSORTI (D, IND, N, B, W)    ! partition workspace
       IMPLICIT NONE
       INTEGER D, IND, N, B, W
       DIMENSION D(N), IND(N), W(2*N)
       
       IF (N < 2) RETURN
       CALL RSI01 (D, IND, N, B, W(1), W(N+1))
       RETURN
      END  ! of RSORTI


*-----------------------------------------------------------------------
      SUBROUTINE RSI01 (D, D2, N, B, W, W2)
       IMPLICIT NONE
       INTEGER D, D2, N, B, W, W2
       DIMENSION D(N), D2(N), W(N), W2(N)

       INTEGER I, J, K, ILIM, P1OLD, P0OLD, P1, P0, S
       LOGICAL INW                         ! data is in W
       INTEGER BITSZ                       ! external function

*          validate
       IF (B .EQ. 0) STOP 'Bit size must not be zero'
       IF (ABS(B) > BITSZ(N)) STOP 'Bit size to large'

*-----------------------------------------------------------------------
*      sort starts on the least significant bit, position 0
*-----------------------------------------------------------------------
       P1 = N+1
       P0 = 0
       DO J = 1, N

*      test for negative values
         IF (D(J) < 0 .AND. B > 0 .AND. B .NE. BITSZ(N)) 
     &       STOP 'Set argument B negative to sort negative numbers'

         IF ( BTEST(D(J), 0) ) THEN
           P1 = P1 - 1
           W(P1) = D(J)
           W2(P1) = D2(J)
         ELSE 
           P0 = P0 + 1
           W(P0) = D(J)
           W2(P0) = D2(J)
         END IF
       END DO
       INW = .TRUE.

*-----------------------------------------------------------------------
*      now alternate between putting data into D and putting data into W
*-----------------------------------------------------------------------
       ILIM = 0
       IF (B < 0) ILIM = MIN( ABS(B)-1, BITSZ(I)-2 )
       IF (B > 0) ILIM = B - 1
       DO I = 1, ILIM
         P1OLD = P1
         P0OLD = P0         ! save the value from previous bit
         P1 = N+1
         P0 = 0                 ! start a fresh count for next bit

*         odd I
         IF (INW) THEN

*          copy data from the zeros
           DO J = 1, P0OLD, 1
             IF ( BTEST(W(J), I) ) THEN
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             ELSE
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             END IF
           END DO
          
*          copy data from the ones
           DO J = N, P1OLD, -1           
             IF ( BTEST(W(J), I) ) THEN
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             ELSE
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             END IF
           END DO
          
*        even I
         ELSE 

*          copy data from the zeros
           DO J = 1, P0OLD, 1
             IF ( BTEST(D(J), I) ) THEN
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             ELSE
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             END IF
           END DO

*            copy data from the ones
           DO J = N, P1OLD, -1
             IF ( BTEST(D(J), I) ) THEN
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             ELSE
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             END IF
           END DO
         END IF ! even or odd i
        
         INW = .NOT. INW
       END DO   ! next i

*-----------------------------------------------------------------------
*        the sign bit
*-----------------------------------------------------------------------
       IF (B < 0) THEN
         P1OLD = P1
         P0OLD = P0
         P1 = N+1
         P0 = 0 

*          if sign bit is set, send to the zero end of the work array
         IF (INW) THEN
           DO J = 1, P0OLD, 1
             IF ( BTEST(W(J), BITSZ(I)-1) ) THEN 
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             ELSE
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             END IF
           END DO          
           DO J = N, P1OLD, -1
             IF ( BTEST(W(J), BITSZ(I)-1) ) THEN
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             ELSE
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             END IF
           END DO
          
         ELSE
           DO J = 1, P0OLD, 1
             IF ( BTEST(D(J), BITSZ(I)-1) ) THEN 
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             ELSE
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             END IF
           END DO          
           DO J = N, P1OLD, -1
             IF ( BTEST(D(J), BITSZ(I)-1) ) THEN
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             ELSE
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             END IF
           END DO
         END IF

         INW = .NOT. INW  
       END IF
          
*-----------------------------------------------------------------------
*      put the data back into the input array
*-----------------------------------------------------------------------
       K = P1
       IF (INW) THEN     
         DO J = 1, P0, 1
           D(J) = W(J)
           D2(J) = W2(J)
         END DO
         DO J = N, P1, -1        
           D(K) = W(J)
           D2(K) = W2(J)
           K = K + 1
         END DO
       
*      Swap values in the ones portion
       ELSE 
         DO J = N, (N+P1+2)/2, -1
           S = D(J)
           D(J) = D(K)
           D(K) = S
           S = D2(J)
           D2(J) = D2(K)
           D2(K) = S
           K = K + 1
         END DO
       END IF
       RETURN
      END ! of RSI01


***********************************************************************
* msortr - merge sort on real numbers
* by Andy Allinger, 2012-2016, released to the  public domain
* This program may be used by any person for any purpose
*
* origin:  John von Neumann, 1945
*
* NOTE:  This routine is a sort-and-carry.  To obtain a permutation
* vector to apply to other arrays, you must first initialize INDA to
* 1, 2, 3...N
*
*__Name_____Type_________I/O______Description_______________________
*  A[N]     Real         Both     The array, returned sorted
*  INDA[N]  Integer      Both     Integer array to carry along, meaning
*                                   permuted in the same fashion as A
*  N        Integer      In       Length of array
*  RW[3@N]  Real         Neither  Workspace
*  IW[3@N]  Integer      Neither  Workspace
***********************************************************************
      SUBROUTINE MSORTR (A, INDA, N, RW, IW)   ! partition workspace
       IMPLICIT NONE
       INTEGER INDA, N, IW
       REAL A, RW
       DIMENSION A(N), INDA(N), RW(3*N), IW(3*N)

       IF (N < 2) RETURN       
       CALL MSR01 (A, INDA, N, RW(1), RW(N+1), RW(2*N+1),
     &  IW(1), IW(N+1), IW(2*N+1))
       RETURN
      END  ! of MSORTR
       
       
***********************************************************************
      SUBROUTINE MSR01 (A, INDA, N, B, C, D, INDB, INDC, INDD)
       IMPLICIT NONE
       INTEGER INDA, N, INDB, INDC, INDD
       REAL A, B, C, D
       DIMENSION A(N), INDA(N), B(N), C(N), D(N), 
     &             INDB(N), INDC(N), INDD(N)
       
*            local variables
       INTEGER g, h, hlim, iA, iB, iC, iD, ilimA, ilimB, ilimC, 
     &   ilimD, j, jA, jB, jC, jD, kA, kB, kC, kD, LA, LB, LC, LD, 
     &   m, o, o0, indx, indy
       INTEGER BITSZ      ! external function
       REAL x, y

*        validate
       IF (n < 2) RETURN
       
*         first pass through data
       jC = 0 
       jD = 0
       DO h = 2, n, 2
         x = A(h-1)
         indx = indA(h-1)
         y = A(h)
         indy = indA(h)
         IF (BTEST(h, 1)) THEN
           IF (x > y) THEN
             jC = jC + 1
             C(jC) = y
             indC(jC) = indy
             jC = jC + 1
             C(jC) = x
             indC(jC) = indx
           ELSE
             jC = jC + 1
             C(jC) = x
             indC(jC) = indx
             jC = jC + 1
             C(jC) = y
             indC(jC) = indy
           END IF
         ELSE
           IF (x > y) THEN
             jD = jD + 1
             D(jD) = y
             indD(jD) = indy
             jD = jD + 1
             D(jD) = x
             indD(jD) = indx
           ELSE
             jD = jD + 1
             D(jD) = x
             indD(jD) = indx
             jD = jD + 1
             D(jD) = y
             indD(jD) = indy
           END IF
         END IF
       END DO ! next h    
       
*          odd n
       IF (BTEST(n, 0)) THEN
         IF (BTEST(n, 1)) THEN
           jD = jD + 1
           D(jD) = A(n)
           indD(jD) = indA(n)
         ELSE
           jC = jC + 1
           C(jC) = A(n)
           indC(jC) = indA(n)
         END IF
       END IF
       LC = jC
       LD = jD

*             number of passes through data
       M = BITSZ(N)
       DO WHILE(.NOT. BTEST(n-1, m-1))
         m = m - 1
       END DO

***********************************************************************
*              sort
***********************************************************************
       o = 2
       LA = 0
       LB = 0                 ! avoid compiler warnings
       DO g = 2, m                      ! begin pass through data
         o0 = o
         o = o * 2 
         hlim = (n-1) / o +1

         IF (BTEST(g, 0)) THEN
           kA = 0
           kB = 0
           jC = 0
           jD = 0

           DO h = 1, hlim   
             ilimA = MIN(o0, LA - kA)
             ilimB = MIN(o0, LB - kB)
             iA = 0
             iB = 0
                    
***********************************************************************
             IF (BTEST(h, 0)) THEN      ! alternate between output arrays
***********************************************************************   
               IF (iA .GE. ilimA) PRINT *, '!!empty buffer A iA=', iA
               iA = iA + 1
               kA = kA + 1
               x = A(kA)
               indx = indA(kA)
            
               IF (iB .GE. ilimB) THEN
                 jC = jC + 1
                 C(jC) = x
                 indC(jC) = indx
                 GO TO 20
               END IF
              
               iB = iB + 1
               kB = kB + 1
               y = B(kB)
               indy = indB(kB)
            
***********************************************************************
*              main comparison statement
***********************************************************************
  10           IF (x > y) THEN              ! send y first
                 jC = jC + 1
                 C(jC) = y
                 indC(jC) = indy
                 IF (iB < ilimB) THEN      ! take another from B[]
                   iB = iB + 1
                   kB = kB + 1
                   y = B(kB)
                   indy = indB(kB)
                   GO TO 10
                 ELSE 
                   jC = jC + 1
                   C(jC) = x
                   indC(jC) = indx
                 END IF
               ELSE                          ! send x first
                 jC = jC + 1
                 C(jC) = x
                 indC(jC) = indx
                 IF (iA < ilimA) THEN        ! take another from A[]
                   iA = iA + 1
                   kA = kA + 1
                   x = A(kA)
                   indx = indA(kA)
                   GO TO 10
                 ELSE 
                   jC = jC + 1
                   C(jC) = y
                   indC(jC) = indy
                 END IF
               END IF
 
***********************************************************************
*           one side of the input is bound to become empty before the other 
***********************************************************************
  20           DO WHILE (iA < ilimA) 
                 iA = iA + 1
                 kA = kA + 1
                 x = A(kA)
                 indx = indA(kA)
                 jC = jC + 1
                 C(jC) = x
                 indC(jC) = indx
               END DO
             
               DO WHILE (iB < ilimB)
                 iB = iB + 1
                 kB = kB + 1
                 y = B(kB)
                 indy = indB(kB)
                 jC = jC + 1
                 C(jC) = y
                 indC(jC) = indy
               END DO
             
***********************************************************************
             ELSE            ! even h
***********************************************************************
               IF (iA .GE. ilimA) PRINT *, 'empty buffer A iA=', iA
               iA = iA + 1
               kA = kA + 1
               x = A(kA)
               indx = indA(kA)
              
               IF (iB .GE. ilimB) THEN
                 jD = jD + 1
                 D(jD) = x
                 indD(jD) = indx
                 GO TO 40
               END IF
              
               iB = iB + 1
               kB = kB + 1
               y = B(kB)
               indy = indB(kB)
              
***********************************************************************
*               main comparison statement
***********************************************************************
  30           IF (x > y) THEN
                 jD = jD + 1
                 D(jD) = y
                 indD(jD) = indy
                 IF (iB < ilimB) THEN
                   iB = iB + 1
                   kB = kB + 1
                   y = B(kB)
                   indy = indB(kB)
                   GO TO 30
                 ELSE 
                   jD = jD + 1
                   D(jD) = x
                   indD(jD) = indx
                 END IF
               ELSE
                 jD = jD + 1
                 D(jD) = x
                 indD(jD) = indx
                 IF (iA < ilimA) THEN
                   iA = iA + 1
                   kA = kA + 1
                   x = A(kA)
                   indx = indA(kA)
                   GO TO 30
                 ELSE 
                   jD = jD + 1
                   D(jD) = y
                   indD(jD) = indy
                 END IF
               END IF
  
***********************************************************************
*           one side of the input is bound to become empty before the other 
***********************************************************************
  40           DO WHILE (iA < ilimA) 
                 iA = iA + 1
                 kA = kA + 1
                 x = A(kA)
                 indx = indA(kA)
                 jD = jD + 1
                 D(jD) = x
                 indD(jD) = indx
               END DO
            
               DO WHILE (iB < ilimB)
                 iB = iB + 1
                 kB = kB + 1
                 y = B(kB)
                 indy = indB(kB)
                 jD = jD + 1
                 D(jD) = y
                 indD(jD) = indy
               END DO
             END IF ! even or odd h
           END DO ! next h
         
*            remember lengths of arrays
           LC = jC
           LD = jD
       
***********************************************************************
         ELSE                   ! even g
***********************************************************************
           kC = 0
           kD = 0
           jA = 0
           jB = 0

           DO h = 1, hlim   
             ilimC = MIN(o0, LC - kC)
             ilimD = MIN(o0, LD - kD)
             iC = 0
             iD = 0
                    
             IF (BTEST(h, 0)) THEN        ! alternate between output arrays
            
               IF (iC .GE. ilimC) PRINT *, '!!empty buffer C iC=', iC
               iC = iC + 1
               kC = kC + 1
               x = C(kC)
               indx = indC(kC)
            
               IF (iD .GE. ilimD) THEN
                 jA = jA + 1
                 A(jA) = x
                 indA(jA) = indx
                 GO TO 60
               END IF
              
               iD = iD + 1
               kD = kD + 1
               y = D(kD)
               indy = indD(kD)
            
***********************************************************************
*              main comparison statement
***********************************************************************
  50           IF (x > y) THEN
                 jA = jA + 1
                 A(jA) = y
                 indA(jA) = indy
                 IF (iD < ilimD) THEN
                   iD = iD + 1
                   kD = kD + 1
                   y = D(kD)
                   indy = indD(kD)
                   GO TO 50
                 ELSE 
                   jA = jA + 1
                   A(jA) = x
                   indA(jA) = indx
                 END IF
               ELSE
                 jA = jA + 1
                 A(jA) = x
                 indA(jA) = indx
                 IF (iC < ilimC) THEN
                   iC = iC + 1
                   kC = kC + 1
                   x = C(kC)
                   indx = indC(kC)
                   GO TO 50
                 ELSE 
                   jA = jA + 1
                   A(jA) = y
                   indA(jA) = indy
                 END IF
               END IF
 
***********************************************************************
*           one side of the input is bound to become empty before the other 
***********************************************************************
  60           DO WHILE (iC < ilimC) 
                 iC = iC + 1
                 kC = kC + 1
                 x = C(kC)
                 indx = indC(kC)
                 jA = jA + 1
                 A(jA) = x
                 indA(jA) = indx
               END DO
             
               DO WHILE (iD < ilimD)
                 iD = iD + 1
                 kD = kD + 1
                 y = D(kD)
                 indy = indD(kD)
                 jA = jA + 1
                 A(jA) = y
                 indA(jA) = indy
               END DO
             
***********************************************************************
             ELSE            ! even h
***********************************************************************
               IF (iC .GE. ilimC) PRINT *, '!!empty buffer C iC=', iC
               iC = iC + 1
               kC = kC + 1
               x = C(kC)
               indx = indC(kC)
              
               IF (iD .GE. ilimD) THEN
                 jB = jB + 1
                 B(jB) = x
                 indB(jB) = indx
                 GO TO 80
               END IF
              
               iD = iD + 1
               kD = kD + 1
               y = D(kD)
               indy = indD(kD)
              
***********************************************************************
*               main comparison statement
***********************************************************************
  70           IF (x > y) THEN
                 jB = jB + 1
                 B(jB) = y
                 indB(jB) = indy
                 IF (iD < ilimD) THEN
                   iD = iD + 1
                   kD = kD + 1
                   y = D(kD)
                   indy = indD(kD)
                   GO TO 70
                 ELSE 
                   jB = jB + 1
                   B(jB) = x
                   indB(jB) = indx
                 END IF
               ELSE
                 jB = jB + 1
                 B(jB) = x
                 indB(jB) = indx
                 IF (iC < ilimC) THEN
                   iC = iC + 1
                   kC = kC + 1
                   x = C(kC)
                   indx = indC(kC)
                   GO TO 70
                 ELSE 
                   jB = jB + 1
                   B(jB) = y
                   indB(jB) = indy
                 END IF
               END IF
  
***********************************************************************
*           one side of the input is bound to become empty before the other 
***********************************************************************
  80           DO WHILE (iC < ilimC) 
                 iC = iC + 1
                 kC = kC + 1
                 x = C(kC)
                 indx = indC(kC)
                 jB = jB + 1
                 B(jB) = x
                 indB(jB) = indx
               END DO
            
               DO WHILE (iD < ilimD)
                 iD = iD + 1
                 kD = kD + 1
                 y = D(kD)
                 indy = indD(kD)
                 jB = jB + 1
                 B(jB) = y
                 indB(jB) = indy
               END DO

             END IF ! even or odd h
           END DO ! next h
         
*            remember lengths of arrays
           LA = jA
           LB = jB
           
         END IF ! even or odd g
       END DO ! next g
       
*         move output back into A
       IF (BTEST(m, 0)) THEN
         DO j = 1, n
           A(j) = C(j)
           indA(j) = indC(j)
         END DO
       END IF
       RETURN
      END ! of MSR01


*-----------------------------------------------------------------------
* bisection search for the first element greater than or equal to a value
* a fragment of SEVAL, by Rondall Jones, Sandia National Laboratory
*                    for unlimited distribution
*-----------------------------------------------------------------------
*   Name        Type     In/Out     Description
*    X[N]       REAL     in         array (in increasing order)
*    N          INTEGER  in         length of X
*    V          REAL     in         the value sought
*    I          INTEGER  out        index where X .GE. V found
*-----------------------------------------------------------------------
      SUBROUTINE RBSGE(X, N, V, I)      
       IMPLICIT NONE
       INTEGER N, I
       REAL X, V
       DIMENSION X(N)

       INTEGER IL, IR
       
*       check input
       IF (N < 1) RETURN

*        V is greater than X(N) or less than X(1)
       IF (V < X(1)) THEN
         I = 1
         RETURN
       ELSE IF (V > X(N)) THEN
         I = N
         RETURN
       END IF

*        bisection search
       IL = 1
       IR = N

  10   I = (IL + IR) / 2
       IF (I .EQ. IL) THEN
         IF (X(I) .GE. V) THEN
           I = IL
         ELSE 
           I = IR
         END IF
         RETURN
       END IF

       IF (X(I) < V) THEN
         IL = I
       ELSE
         IR = I
       END IF 
       GO TO 10
       
      END ! of RBSGE


*-----------------------------------------------------------------------
* bisection search an array for the last element less than or equal to a value
* a fragment of SEVAL, by Rondall Jones, Sandia National Laboratory
*                    for unlimited distribution
*
*  note that type is changed to INTEGER and .GE. is changed to .LE.
*-----------------------------------------------------------------------
*   Name        Type     In/Out     Description
*    J[N]       INTEGER  in         array (in increasing order)
*    N          INTEGER  in         length of J
*    V          INTEGER  in         the value sought
*    I          INTEGER  out        index where J .LE. V found
*                                    returns 0 when J[1] > V
*-----------------------------------------------------------------------
      SUBROUTINE IBSLE(J, N, V, I)
       IMPLICIT NONE
       INTEGER J, N, V, I
       DIMENSION J(N)

       INTEGER IL, IR
       
*       check input
       IF (N < 1) RETURN

*         V is greater than J[N] or less than J[1]
       IF (V < J(1)) THEN
         I = 0
         RETURN
       ELSE IF (V > J(N)) THEN
         I = N
         RETURN
       END IF

*        bisection search
       IL = 1
       IR = N

  10   I = (IL + IR) / 2
       IF (I .EQ. IL) THEN
         IF (J(IR) .LE. V) THEN
           I = IR
         ELSE 
           I = IL
         END IF
         RETURN
       END IF

       IF (J(I) .LE. V) THEN 
         IL = I
       ELSE 
         IR = I
       END IF
       GO TO 10
       
      END ! of IBSLE


*-----------------------------------------------------------------------
* Bisection search for the first element greater than a value
*
*___Name_______Type_____In/Out_____Description__________________________
*   X(N)       REAL     In         Array (in increasing order)
*   N          INTEGER  In         Length of X
*   V          REAL     In         The value sought
*   I          INTEGER  Out        Index where X > V found
*-----------------------------------------------------------------------
      SUBROUTINE RBSGT (X, N, V, I)
       IMPLICIT NONE
       INTEGER N, I
       REAL X, V
       DIMENSION X(N)

       INTEGER IL, IR

*       check input
       IF (N < 1) RETURN

*        V is greater than X(N) or less than X(1)
       IF (V < X(1)) THEN
         I = 1
         RETURN
       ELSE IF (V > X(N)) THEN
         I = N
         RETURN
       END IF

*        bisection search
       IL = 1
       IR = N

  10   I = (IL + IR) / 2
       IF (I .EQ. IL) THEN
         IF (X(I) > V) THEN
           I = IL
         ELSE
           I = IR
         END IF
         RETURN
       END IF

       IF (X(I) .LE. V) THEN
         IL = I
       ELSE
         IR = I
       END IF
       GO TO 10
      END  ! of RBSGT

