      SUBROUTINE RANGEMENT(NEQ)
      INCLUDE 'com_faces.f'
C      REAL*8 SEC0,SEC1
C
C      CALL TEMPS(SEC0,IS)
C      CALL MYSORT_PLUS(CENTR,CENTRMI,CENTRMA,IORDRE,NEQ)
      CALL MYSORT(CENTR,IORDRE,NEQ)
C      CALL TEMPS(SEC1,IS)
C      PRINT*,'Sort    =',SEC1-SEC0
      END
C-----------------------------------------------------------------------
      SUBROUTINE MYSORT(ARRAY,IARRAY,NEQ)
CCBEGIN
C-----------------------------------------------------------------------
C
C
C     PURPOSE
C     ------
C
C     CLASSE LES ELEMENTS DU TABLEAU DE REEL "ARRAY"
C     ET CONSERVE LA TRACE DE LA PERMUTATION DANS LE TABLEAU
C     D'ENTIERS "IARRAY".
C
CCDOC
C     INPUT ARGUMENTS
C     --------------
C
C     NEQ     = LENGTH OF IARRAY.
C
C     INPUT-OUTPUT ARGUMENTS
C     ----------------------
C
C     IARRAY  = INTEGER ARRAY.
C
C     REFERENCE
C     ---------
C
C     LORIN H.,SORTING AND SORTING SYSTEMS,
C     ADDISON-WESLEY PUBLISHING COMPANY,1975.A-42.
C     ACM ALGORITHM 347.
C
C     NOTES
C     -----
C
C     INTERNAL WORK AREA IWR(2,K) PERMITS SORTING UP TO
C     2**(K+1)-1 ELEMENTS. HERE K IS SET TO 32.
C
C-----------------------------------------------------------------------
C
CCEND
C
C
      DIMENSION IARRAY(NEQ),IWR(2,32)
      REAL*4 ARRAY(NEQ)
C
C     INITIALISATION
C
      DO 200 I=1,NEQ
        IARRAY(I)=I
  200 CONTINUE
      M = 1
      I = 1
      J = NEQ
    5 IF(I.GE.J)  GO TO 70
C
C--------BEGIN SPLIT OF SEGMENT--------------------------------------
C
C        FIND MEDIAN OF FIRST, MIDDLE AND LAST ELEMENTS OF SEGMENT, IMID
C
   10    K = I
         IJ2 = (J+I)/2
         IMID = IARRAY(IJ2)
         AMID = ARRAY(IJ2)
         IF(ARRAY(I).LE.AMID)  GOTO 20
            IARRAY(IJ2) = IARRAY(I)
            IARRAY(I) = IMID
            IMID = IARRAY(IJ2)
            ARRAY(IJ2) = ARRAY(I)
            ARRAY(I) = AMID
            AMID = ARRAY(IJ2)
   20    L = J
         IF(ARRAY(J).GE.AMID)  GO TO 40
            IARRAY(IJ2) = IARRAY(J)
            IARRAY(J) = IMID
            IMID = IARRAY(IJ2)
            ARRAY(IJ2) = ARRAY(J)
            ARRAY(J) = AMID
            AMID = ARRAY(IJ2)
            IF(ARRAY(I).LE.AMID)  GO TO 40
               IARRAY(IJ2) = IARRAY(I)
               IARRAY(I) = IMID
               IMID = IARRAY(IJ2)
               ARRAY(IJ2) = ARRAY(I)
               ARRAY(I) = AMID
               AMID = ARRAY(IJ2)
            GO TO 40
C
C        SPLIT SEGMENT AROUND IMID, SORTING ELEMENTS INTO THOSE
C        GT. AND LT. IMID.
C
   30    IARRAY(L) = IARRAY(K)
         IARRAY(K) = ITT
         ARRAY(L) = ARRAY(K)
         ARRAY(K) = ATT
   40    L = L - 1
         IF(ARRAY(L).GT.AMID)  GO TO 40
         ITT = IARRAY(L)
         ATT = ARRAY(L)
   50    K = K + 1
         IF(ARRAY(K).LT.AMID)  GO TO 50
         IF(K.LE.L)  GO TO 30
C
C--------END OF SPLIT OF SEGMENT---------------------------------------
C
C     STORE LIMITS OF LARGER SEGMENT IN IWR
C
      IF(L-I.LE.J-K)  GO TO 60
         IWR(2,M) = I
         IWR(1,M) = L
         I = K
         M = M + 1
         GO TO 80
   60 IWR(2,M) = K
      IWR(1,M) = J
      J = L
      M = M + 1
      GO TO 80
C
C     WHEN SEGMENT HAS BEEN SORTED, GET LIMITS OF NEXT SEGMENT.
C     WHEN ALL SEGMENTS HAVE BEEN SORTED, SORT IS COMPLETE.
C
   70 M = M - 1
      IF(M.EQ.0)  RETURN
      I = IWR(2,M)
      J = IWR(1,M)
C
C     SORT SHORT SEQUENCES (LT.11 ELEMENTS) BY INTERCHANGE OF ADJACENT
C     PAIRS (EXCEPT INITIAL SEGMENT)
C
   80 IF(J-I.GE.11)  GO TO 10
      IF(I.EQ.1)  GO TO 5
      I = I - 1
   90 I = I + 1
      IF(I.EQ.J)  GO TO 70
      IMID = IARRAY(I+1)
      AMID = ARRAY(I+1)
      IF(ARRAY(I).LE.AMID)  GO TO 90
      K = I
  100 IARRAY(K+1) = IARRAY(K)
      ARRAY(K+1) = ARRAY(K)
      K = K - 1
      IF(AMID.LT.ARRAY(K))  GO TO 100
      IARRAY(K+1) = IMID
      ARRAY(K+1) = AMID
      GO TO 90
C
      END
C-----------------------------------------------------------------------
      SUBROUTINE MYSORT_PLUS(ARRAY,ARRAYMI,ARRAYMA,IARRAY,NEQ)
CCBEGIN
C-----------------------------------------------------------------------
C
C
C     PURPOSE
C     ------
C
C     CLASSE LES ELEMENTS DU TABLEAU DE REEL "ARRAY"
C     ET CONSERVE LA TRACE DE LA PERMUTATION DANS LE TABLEAU
C     D'ENTIERS "IARRAY".
C
CCDOC
C     INPUT ARGUMENTS
C     --------------
C
C     NEQ     = LENGTH OF IARRAY.
C
C     INPUT-OUTPUT ARGUMENTS
C     ----------------------
C
C     IARRAY  = INTEGER ARRAY.
C
C     REFERENCE
C     ---------
C
C     LORIN H.,SORTING AND SORTING SYSTEMS,
C     ADDISON-WESLEY PUBLISHING COMPANY,1975.A-42.
C     ACM ALGORITHM 347.
C
C     NOTES
C     -----
C
C     INTERNAL WORK AREA IWR(2,K) PERMITS SORTING UP TO
C     2**(K+1)-1 ELEMENTS. HERE K IS SET TO 32.
C
C-----------------------------------------------------------------------
C
CCEND
C
C
      INTEGER IARRAY(*),IWR(2,32)
      DIMENSION ARRAY(*),ARRAYMI(*),ARRAYMA(*)
      logical pluspetit,plusgrand
      include 'com_sort.f'
C
C INITIALISATION
C
      nboper = 0
      nbfacile = 0
      DO I=1,NEQ
        IARRAY(I)=I
      ENDDO
      M = 1
      I = 1
      J = NEQ
      do while (i.eq.1)
        IF (I.lt.J) then
C
C--------BEGIN SPLIT OF SEGMENT--------------------------------------
C
C FIND MEDIAN OF FIRST, MIDDLE AND LAST ELEMENTS OF SEGMENT, IMID
C
 10       K = I
          IJ2 = (J+I)/2
          call affectecour(array,arraymi,arrayma,iarray,ij2)
          if (plusgrand(array,arraymi,arrayma,i))
     &         call echange(array,arraymi,arrayma,iarray,i,ij2)
          L = J
          IF (pluspetit(array,arraymi,arrayma,j)) then
            call echange(array,arraymi,arrayma,iarray,j,ij2)
            if (plusgrand(array,arraymi,arrayma,i)) then
              call echange(array,arraymi,arrayma,iarray,i,ij2)
              GOTO 40
C
C SPLIT SEGMENT AROUND IMID, SORTING ELEMENTS INTO THOSE
C GT. AND LT. IMID.
C
Cfj              call affecte(array,arraymi,arrayma,iarray,l,k)
Cfj              call affectetmpinv(array,arraymi,arrayma,iarray,k)
            endif
          endif
 40       L = L-1
          do while (plusgrand(array,arraymi,arrayma,l))
            L = L-1
          enddo
          call affectetmp(array,arraymi,arrayma,iarray,l)
          K = K+1
          do while (pluspetit(array,arraymi,arrayma,k))
            K = K+1
          enddo
          IF (K.LE.L) then
            call affecte(array,arraymi,arrayma,iarray,l,k)
            call affectetmpinv(array,arraymi,arrayma,iarray,k)
            goto 40
          endif
C
C--------END OF SPLIT OF SEGMENT---------------------------------------
C
C     STORE LIMITS OF LARGER SEGMENT IN IWR
C
          IF (L-I.gt.J-K) then
            IWR(2,M) = I
            IWR(1,M) = L
            I = K
            M = M + 1
            IF (J-I.GE.11) GOTO 10
          endif
          IWR(2,M) = K
          IWR(1,M) = J
          J = L
          M = M + 1
          IF (J-I.GE.11) GOTO 10
        endif
C
C WHEN SEGMENT HAS BEEN SORTED, GET LIMITS OF NEXT SEGMENT.
C WHEN ALL SEGMENTS HAVE BEEN SORTED, SORT IS COMPLETE.
C
 70     M = M-1
        IF (M.EQ.0) then
          print*,'Operations',nboper,nbfacile
     &         ,100.*real(nbfacile)/real(nboper)
          RETURN
        endif
        I = IWR(2,M)
        J = IWR(1,M)
C
C SORT SHORT SEQUENCES (LT.11 ELEMENTS) BY INTERCHANGE OF ADJACENT
C PAIRS (EXCEPT INITIAL SEGMENT)
C
        IF (J-I.GE.11) GOTO 10
      enddo
c
      I = I-1
 90   I = I+1
      IF (I.EQ.J) GOTO 70
      call affectecour(array,arraymi,arrayma,iarray,i+1)
      if ( .not.(pluspetit(array,arraymi,arrayma,i)) ) then
        k = i
        call affecte(array,arraymi,arrayma,iarray,k+1,k)
        k = k-1
        do while ( plusgrand(array,arraymi,arrayma,k) )
          call affecte(array,arraymi,arrayma,iarray,k+1,k)
          k = k-1
        enddo
        call affectecourinv(array,arraymi,arrayma,iarray,k+1)
      endif
      GOTO 90
C
      END
C-----------------------------------------------------------------------
      logical function pluspetit(array,arraymi,arrayma,i)
      dimension array(*),arraymi(*),arrayma(*)
      include 'com_sort.f'
c
      pluspetit = (array(i).lt.acour)
      nboper = nboper+1
      if (array(i).lt.acour.and.arrayma(i).lt.acourmi)
     &     nbfacile=nbfacile+1
      end
C-----------------------------------------------------------------------
      logical function plusgrand(array,arraymi,arrayma,i)
      dimension array(*),arraymi(*),arrayma(*)
      include 'com_sort.f'
c
      plusgrand = (array(i).gt.acour)
      nboper = nboper+1
      if (array(i).gt.acour.and.arraymi(i).gt.acourma)
     &     nbfacile=nbfacile+1
      end
C-----------------------------------------------------------------------
      subroutine echange(array,arraymi,arrayma,iarray,i,j)
      dimension array(*),arraymi(*),arrayma(*)
      integer   iarray(*)
      include 'com_sort.f'
c
      iarray(j) = iarray(i)
      iarray(i) = icour
      icour     = iarray(j)
      array(j)  = array(i)
      array(i)  = acour
      acour     = array(j)
      arraymi(j) = arraymi(i)
      arraymi(i) = acourmi
      acourmi    = arraymi(j)
      arrayma(j) = arrayma(i)
      arrayma(i) = acourma
      acourma    = arrayma(j)
      end
C-----------------------------------------------------------------------
      subroutine affecte(array,arraymi,arrayma,iarray,i,j)
      dimension array(*),arraymi(*),arrayma(*)
      integer   iarray(*)
      include 'com_sort.f'
c
      iarray(i) = iarray(j)
      array(i)  = array(j)
      arraymi(i)  = arraymi(j)
      arrayma(i)  = arrayma(j)
      end
C-----------------------------------------------------------------------
      subroutine affectecour(array,arraymi,arrayma,iarray,i)
      dimension array(*),arraymi(*),arrayma(*)
      integer   iarray(*)
      include 'com_sort.f'
c
      icour = iarray(i)
      acour = array(i)
      acourmi = arraymi(i)
      acourma = arrayma(i)
      end
C-----------------------------------------------------------------------
      subroutine affectecourinv(array,arraymi,arrayma,iarray,i)
      dimension array(*),arraymi(*),arrayma(*)
      integer   iarray(*)
      include 'com_sort.f'
c
      iarray(i) = icour
      array(i)  = acour
      arraymi(i)  = acourmi
      arrayma(i)  = acourma
      end
C-----------------------------------------------------------------------
      subroutine affectetmp(array,arraymi,arrayma,iarray,i)
      dimension array(*),arraymi(*),arrayma(*)
      integer   iarray(*)
      include 'com_sort.f'
c
      itmp = iarray(i)
      atmp = array(i)
      atmpmi = arraymi(i)
      atmpma = arrayma(i)
      end
C-----------------------------------------------------------------------
      subroutine affectetmpinv(array,arraymi,arrayma,iarray,i)
      dimension array(*),arraymi(*),arrayma(*)
      integer   iarray(*)
      include 'com_sort.f'
c
      iarray(i) = itmp
      array(i)  = atmp
      arraymi(i)  = atmpmi
      arrayma(i)  = atmpma
      end
C-----------------------------------------------------------------------
      SUBROUTINE CORRIGE(IORDRE,NBON,NSURF,NEIS,NSENS)
      INTEGER IORDRE(*),NEIS(*),NSENS(*)
ccc      real*8 S1,S2
C
      N = NBON+NSURF
      NPROF = 1000
ccc      CALL TEMPS(S1,IS)
      DO I=2,N
        IF (IORDRE(I).LE.NBON) THEN
          IF (NSENS(IORDRE(I)).LT.0) THEN
            NFIN = MAX(1,I-NPROF)
            JSURF = 0
            DO J=I-1,NFIN,-1
              IF (IORDRE(J).LE.NBON) THEN
                GOTO 1
              ELSE
                IF (NEIS(IORDRE(J)-NBON).EQ.-NSENS(IORDRE(I))) JSURF=J
              ENDIF
            ENDDO
 1          JJ = JSURF
            IF (JJ.NE.0) THEN
cc                print*,IORDRE(I),' recule de',i,' a',jj,nsens(IORDRE(I))
              IBID = IORDRE(I)
              DO J=I,JJ+1,-1
                IORDRE(J) = IORDRE(J-1)
              ENDDO
              IORDRE(JJ) = IBID
            ENDIF
          ENDIF
        ENDIF
      ENDDO
      DO I=N-1,1,-1
        IF (IORDRE(I).LE.NBON) THEN
          IF (NSENS(IORDRE(I)).GT.0) THEN
            NFIN = MIN(N,I+NPROF)
            JSURF = 0
            DO J=I+1,NFIN
              IF (IORDRE(J).LE.NBON) THEN
                GOTO 2
              ELSE
                IF (NEIS(IORDRE(J)-NBON).EQ.NSENS(IORDRE(I))) JSURF=J
              ENDIF
            ENDDO
 2          JJ = JSURF
            IF (JJ.NE.0) THEN
cc                print*,IORDRE(I),' avance de',i,' a',jj,nsens(IORDRE(I))
              IBID = IORDRE(I)
              DO J=I,JJ-1
                IORDRE(J) = IORDRE(J+1)
              ENDDO
              IORDRE(JJ) = IBID
            ENDIF
          ENDIF
        ENDIF
      ENDDO
ccc      CALL TEMPS(S2,IS)
ccc      print*,real(s2-s1)
      END
