subroutine SortR8(N1,N2,M1,M2,RA)
            integer             N1
            integer             N2
            integer             M1
            integer             M2
            double precision    RA(M1:M2)

        double precision        RRA

        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 = max(N1,IR-1)   ! 'max' needed to sort one element (N1=N2)
            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) then
                if (RA(J) .lt. RA(J+1)) J=J+1
            end if
            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