Module m_indmed
Integer, Parameter :: kdp = selected_real_kind(15)
public :: indmed
private :: kdp
private :: R_indmed, I_indmed, D_indmed
private :: r_med, i_med, d_med
Integer, Allocatable, Dimension(:), Private, Save :: IDONT
interface indmed
  module procedure d_indmed, r_indmed, i_indmed
end interface indmed
contains

Subroutine D_indmed (XDONT, INDM)
!  Returns index of median value of XDONT.
! __________________________________________________________
      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
      Integer, Intent (Out) :: INDM
! __________________________________________________________
      Integer :: IDON
!
      Allocate (IDONT (SIZE(XDONT)))
      Do IDON = 1, SIZE(XDONT)
         IDONT (IDON) = IDON
      End Do
!
      Call d_med (XDONT, IDONT, INDM)
!
      Deallocate (IDONT)
End Subroutine D_indmed
   Recursive Subroutine d_med (XDATT, IDATT, ires_med)
!  Finds the index of the median of XDONT using the recursive procedure
!  described in Knuth, The Art of Computer Programming,
!  vol. 3, 5.3.3 - This procedure is linear in time, and
!  does not require to be able to interpolate in the
!  set as the one used in INDNTH. It also has better worst
!  case behavior than INDNTH, but is about 30% slower in
!  average for random uniformly distributed values.
! __________________________________________________________
      Real (kind=kdp), Dimension (:), Intent (In) :: XDATT
      Integer, Dimension (:), Intent (In) :: IDATT
      Integer, Intent (Out):: ires_med
! __________________________________________________________
!
      Real (kind=kdp), Parameter :: XHUGE = HUGE (XDATT)
      Real (kind=kdp) :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
      Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
      Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
      Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
      Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
      Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
      NDAT = Size (IDATT)
      NMED = (NDAT+1) / 2
      IWRKT = IDATT
!
!  If the number of values is small, then use insertion sort
!
     If (NDAT < 35) Then
!
!  Bring minimum to first location to save test in decreasing loop
!
         IDCR = NDAT
         If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
            IWRK = IWRKT (1)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT (1)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, NDAT - 2
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                XWRK = XWRK1
                IWRK = IWRK1
            Endif
         End Do
         IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
         Do ICRS = 3, NMED
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            IDCR = ICRS - 1
            Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
!  Insert any value less than the current median in the first half
!
         XWRK1 = XDATT (IWRKT (NMED))
         Do ICRS = NMED+1, NDAT
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            If (XWRK < XWRK1) Then
               IDCR = NMED - 1
               Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
               End Do
               IWRKT (IDCR+1) = IWRK
               XWRK1 = XDATT (IWRKT (NMED))
            End If
         End Do
         ires_med = IWRKT (NMED)
         Return
      End If
!
!  Make sorted subsets of 7 elements
!  This is done by a variant of insertion sort where a first
!  pass is used to bring the smallest element to the first position
!  decreasing disorder at the same time, so that we may remove
!  remove the loop test in the insertion loop.
!
      IMAX = 1
      IMIN = 1
      XMAX = XDATT (IWRKT(IMAX))
      XMIN = XDATT (IWRKT(IMIN))
      DO IDEB = 1, NDAT-6, 7
         IDCR = IDEB + 6
         If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
            IWRK = IWRKT(IDEB)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT(IDEB)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, 5
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                IWRK = IWRK1
                XWRK = XWRK1
            Endif
         End Do
         IWRKT (IDEB) = IWRK
         If (XWRK < XMIN) Then
             IMIN = IWRK
             XMIN = XWRK
         End If
         Do ICRS = IDEB+1, IDEB+5
            IWRK = IWRKT (ICRS+1)
            XWRK = XDATT (IWRK)
            IDON = IWRKT(ICRS)
            If (XWRK < XDATT(IDON)) Then
               IWRKT (ICRS+1) = IDON
               IDCR = ICRS
               IWRK1 = IWRKT (IDCR-1)
               XWRK1 = XDATT (IWRK1)
               Do
                  If (XWRK >= XWRK1) Exit
                  IWRKT (IDCR) = IWRK1
                  IDCR = IDCR - 1
                  IWRK1 = IWRKT (IDCR-1)
                  XWRK1 = XDATT (IWRK1)
               End Do
               IWRKT (IDCR) = IWRK
            EndIf
         End Do
         If (XWRK > XMAX) Then
             IMAX = IWRK
             XMAX = XWRK
         End If
      End Do
!
!  Add-up alternatively MAX and MIN values to make the number of data
!  an exact multiple of 7.
!
      IDEB = 7 * (NDAT/7)
      NTRI = NDAT
      If (IDEB < NDAT) Then
!
         Do ICRS = IDEB+1, NDAT
            XWRK1 = XDATT (IWRKT (ICRS))
            IF (XWRK1 > XMAX) Then
               IMAX = IWRKT (ICRS)
               XMAX = XWRK1
            End If
            IF (XWRK1 < XMIN) Then
               IMIN = IWRKT (ICRS)
               XMIN = XWRK1
            End If
         End Do
         IWRK1 = IMAX
         Do ICRS = NDAT+1, IDEB+7
               IWRKT (ICRS) = IWRK1
               If (IWRK1 == IMAX) Then
                  IWRK1 = IMIN
               Else
                  NMED = NMED + 1
                  IWRK1 = IMAX
               End If
         End Do
!
         Do ICRS = IDEB+2, IDEB+7
            IWRK = IWRKT (ICRS)
            XWRK = XDATT (IWRK)
            Do IDCR = ICRS - 1, IDEB+1, - 1
               If (XWRK >= XDATT (IWRKT(IDCR))) Exit
               IWRKT (IDCR+1) = IWRKT (IDCR)
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
         NTRI = IDEB+7
      End If
!
!  Make the set of the indices of median values of each sorted subset
!
         IDON1 = 0
         Do IDON = 1, NTRI, 7
            IDON1 = IDON1 + 1
            IMEDT (IDON1) = IWRKT (IDON + 3)
         End Do
!
!  Find XMED7, the median of the medians
!
         Call d_med (XDATT, IMEDT(1:IDON1), IMED7)
         XMED7 = XDATT (IMED7)
!
!  Count how many values are not higher than (and how many equal to) XMED7
!  This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
!  subsets where the median is lower than the median of medians. For similar
!  reasons, we also have at least 2N/7 values not lower than XMED7. At the
!  same time, we find in each subset the index of the last value < XMED7,
!  and that of the first > XMED7. These indices will be used to restrict the
!  search for the median as the Kth element in the subset (> or <) where
!  we know it to be.
!
         IDON1 = 1
         NLEQ = 0
         NEQU = 0
         Do IDON = 1, NTRI, 7
            IMED = IDON+3
            If (XDATT (IWRKT (IMED)) > XMED7) Then
                  IMED = IMED - 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                  IMED = IMED + 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Endif
            If (XDATT (IWRKT (IMED)) > XMED7) Then
               NLEQ = NLEQ + IMED - IDON
               IENDT (IDON1) = IMED - 1
               ISTRT (IDON1) = IMED
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
               NLEQ = NLEQ + IMED - IDON + 1
               IENDT (IDON1) = IMED
               ISTRT (IDON1) = IMED + 1
            Else                    !       If (XDATT (IWRKT (IMED)) == XMED7)
               NLEQ = NLEQ + IMED - IDON + 1
               NEQU = NEQU + 1
               IENDT (IDON1) = IMED - 1
               Do IMED1 = IMED - 1, IDON, -1
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     IENDT (IDON1) = IMED1 - 1
                  Else
                     Exit
                  End If
               End Do
               ISTRT (IDON1) = IMED + 1
               Do IMED1 = IMED + 1, IDON + 6
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     NLEQ = NLEQ + 1
                     ISTRT (IDON1) = IMED1 + 1
                  Else
                     Exit
                  End If
               End Do
            Endif
            IDON1 = IDON1 + 1
         End Do
!
!  Carry out a partial insertion sort to find the Kth smallest of the
!  large values, or the Kth largest of the small values, according to
!  what is needed.
!
!
         If (NLEQ - NEQU + 1 <= NMED) Then
            If (NLEQ < NMED) Then   !      Not enough low values
                IWRK1 = IMAX
                XWRK1 = XDATT (IWRK1)
                NORD = NMED - NLEQ
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                IDCR = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) < XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT (ICRS1)
                            XWRK1 = XDATT (IWRK1)
                         Else
                           If (ICRS2 < NORD) Then
                              IWRKT (ICRS1) = IWRKT (ICRS)
                              IWRK1 = IWRKT (ICRS1)
                              XWRK1 = XDATT (IWRK1)
                           Endif
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT (ICRS1)
                         XWRK1 = XDATT (IWRK1)
                      End Do
                   End If
                End Do
                ires_med = IWRK1
                Return
            Else
                ires_med = IMED7
                Return
            End If
         Else                       !      If (NLEQ > NMED)
!                                          Not enough high values
                XWRK1 = -XHUGE
                NORD = NLEQ - NEQU - NMED + 1
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
!
                      Do ICRS = IDON, IENDT (IDON1)
                         If (XDATT(IWRKT (ICRS)) > XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            IDCR = ICRS1 - 1
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT(ICRS1)
                            XWRK1 = XDATT(IWRK1)
                         Else
                            If (ICRS2 < NORD) Then
                               IWRKT (ICRS1) = IWRKT (ICRS)
                               IWRK1 = IWRKT(ICRS1)
                               XWRK1 = XDATT(IWRK1)
                            End If
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = IENDT (IDON1), IDON, -1
                         If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         IDCR = ICRS1 - 1
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT(ICRS1)
                         XWRK1 = XDATT(IWRK1)
                      End Do
                   Endif
                End Do
!
                ires_med = IWRK1
                Return
         End If
!
   END Subroutine d_med
!
Subroutine R_indmed (XDONT, INDM)
!  Returns index of median value of XDONT.
! __________________________________________________________
      Real, Dimension (:), Intent (In) :: XDONT
      Integer, Intent (Out) :: INDM
! __________________________________________________________
      Integer :: IDON
!
      Allocate (IDONT (SIZE(XDONT)))
      Do IDON = 1, SIZE(XDONT)
         IDONT (IDON) = IDON
      End Do
!
      Call r_med (XDONT, IDONT, INDM)
!
      Deallocate (IDONT)
End Subroutine R_indmed
   Recursive Subroutine r_med (XDATT, IDATT, ires_med)
!  Finds the index of the median of XDONT using the recursive procedure
!  described in Knuth, The Art of Computer Programming,
!  vol. 3, 5.3.3 - This procedure is linear in time, and
!  does not require to be able to interpolate in the
!  set as the one used in INDNTH. It also has better worst
!  case behavior than INDNTH, but is about 30% slower in
!  average for random uniformly distributed values.
! __________________________________________________________
      Real, Dimension (:), Intent (In) :: XDATT
      Integer, Dimension (:), Intent (In) :: IDATT
      Integer, Intent (Out) :: ires_med
! __________________________________________________________
!
      Real, Parameter :: XHUGE = HUGE (XDATT)
      Real :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
      Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
      Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
      Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
      Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
      Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
      NDAT = Size (IDATT)
      NMED = (NDAT+1) / 2
      IWRKT = IDATT
!
!  If the number of values is small, then use insertion sort
!
     If (NDAT < 35) Then
!
!  Bring minimum to first location to save test in decreasing loop
!
         IDCR = NDAT
         If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
            IWRK = IWRKT (1)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT (1)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, NDAT - 2
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                XWRK = XWRK1
                IWRK = IWRK1
            Endif
         End Do
         IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
         Do ICRS = 3, NMED
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            IDCR = ICRS - 1
            Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
!  Insert any value less than the current median in the first half
!
         XWRK1 = XDATT (IWRKT (NMED))
         Do ICRS = NMED+1, NDAT
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            If (XWRK < XWRK1) Then
               IDCR = NMED - 1
               Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
               End Do
               IWRKT (IDCR+1) = IWRK
               XWRK1 = XDATT (IWRKT (NMED))
            End If
         End Do
         ires_med = IWRKT (NMED)
         Return
      End If
!
!  Make sorted subsets of 7 elements
!  This is done by a variant of insertion sort where a first
!  pass is used to bring the smallest element to the first position
!  decreasing disorder at the same time, so that we may remove
!  remove the loop test in the insertion loop.
!
      IMAX = 1
      IMIN = 1
      XMAX = XDATT (IWRKT(IMAX))
      XMIN = XDATT (IWRKT(IMIN))
      DO IDEB = 1, NDAT-6, 7
         IDCR = IDEB + 6
         If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
            IWRK = IWRKT(IDEB)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT(IDEB)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, 5
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                IWRK = IWRK1
                XWRK = XWRK1
            Endif
         End Do
         IWRKT (IDEB) = IWRK
         If (XWRK < XMIN) Then
             IMIN = IWRK
             XMIN = XWRK
         End If
         Do ICRS = IDEB+1, IDEB+5
            IWRK = IWRKT (ICRS+1)
            XWRK = XDATT (IWRK)
            IDON = IWRKT(ICRS)
            If (XWRK < XDATT(IDON)) Then
               IWRKT (ICRS+1) = IDON
               IDCR = ICRS
               IWRK1 = IWRKT (IDCR-1)
               XWRK1 = XDATT (IWRK1)
               Do
                  If (XWRK >= XWRK1) Exit
                  IWRKT (IDCR) = IWRK1
                  IDCR = IDCR - 1
                  IWRK1 = IWRKT (IDCR-1)
                  XWRK1 = XDATT (IWRK1)
               End Do
               IWRKT (IDCR) = IWRK
            EndIf
         End Do
         If (XWRK > XMAX) Then
             IMAX = IWRK
             XMAX = XWRK
         End If
      End Do
!
!  Add-up alternatively MAX and MIN values to make the number of data
!  an exact multiple of 7.
!
      IDEB = 7 * (NDAT/7)
      NTRI = NDAT
      If (IDEB < NDAT) Then
!
         Do ICRS = IDEB+1, NDAT
            XWRK1 = XDATT (IWRKT (ICRS))
            IF (XWRK1 > XMAX) Then
               IMAX = IWRKT (ICRS)
               XMAX = XWRK1
            End If
            IF (XWRK1 < XMIN) Then
               IMIN = IWRKT (ICRS)
               XMIN = XWRK1
            End If
         End Do
         IWRK1 = IMAX
         Do ICRS = NDAT+1, IDEB+7
               IWRKT (ICRS) = IWRK1
               If (IWRK1 == IMAX) Then
                  IWRK1 = IMIN
               Else
                  NMED = NMED + 1
                  IWRK1 = IMAX
               End If
         End Do
!
         Do ICRS = IDEB+2, IDEB+7
            IWRK = IWRKT (ICRS)
            XWRK = XDATT (IWRK)
            Do IDCR = ICRS - 1, IDEB+1, - 1
               If (XWRK >= XDATT (IWRKT(IDCR))) Exit
               IWRKT (IDCR+1) = IWRKT (IDCR)
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
         NTRI = IDEB+7
      End If
!
!  Make the set of the indices of median values of each sorted subset
!
         IDON1 = 0
         Do IDON = 1, NTRI, 7
            IDON1 = IDON1 + 1
            IMEDT (IDON1) = IWRKT (IDON + 3)
         End Do
!
!  Find XMED7, the median of the medians
!
         Call r_med (XDATT, IMEDT(1:IDON1), IMED7)
         XMED7 = XDATT (IMED7)
!
!  Count how many values are not higher than (and how many equal to) XMED7
!  This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
!  subsets where the median is lower than the median of medians. For similar
!  reasons, we also have at least 2N/7 values not lower than XMED7. At the
!  same time, we find in each subset the index of the last value < XMED7,
!  and that of the first > XMED7. These indices will be used to restrict the
!  search for the median as the Kth element in the subset (> or <) where
!  we know it to be.
!
         IDON1 = 1
         NLEQ = 0
         NEQU = 0
         Do IDON = 1, NTRI, 7
            IMED = IDON+3
            If (XDATT (IWRKT (IMED)) > XMED7) Then
                  IMED = IMED - 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                  IMED = IMED + 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Endif
            If (XDATT (IWRKT (IMED)) > XMED7) Then
               NLEQ = NLEQ + IMED - IDON
               IENDT (IDON1) = IMED - 1
               ISTRT (IDON1) = IMED
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
               NLEQ = NLEQ + IMED - IDON + 1
               IENDT (IDON1) = IMED
               ISTRT (IDON1) = IMED + 1
            Else                    !       If (XDATT (IWRKT (IMED)) == XMED7)
               NLEQ = NLEQ + IMED - IDON + 1
               NEQU = NEQU + 1
               IENDT (IDON1) = IMED - 1
               Do IMED1 = IMED - 1, IDON, -1
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     IENDT (IDON1) = IMED1 - 1
                  Else
                     Exit
                  End If
               End Do
               ISTRT (IDON1) = IMED + 1
               Do IMED1 = IMED + 1, IDON + 6
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     NLEQ = NLEQ + 1
                     ISTRT (IDON1) = IMED1 + 1
                  Else
                     Exit
                  End If
               End Do
            Endif
            IDON1 = IDON1 + 1
         End Do
!
!  Carry out a partial insertion sort to find the Kth smallest of the
!  large values, or the Kth largest of the small values, according to
!  what is needed.
!
!
         If (NLEQ - NEQU + 1 <= NMED) Then
            If (NLEQ < NMED) Then   !      Not enough low values
                IWRK1 = IMAX
                XWRK1 = XDATT (IWRK1)
                NORD = NMED - NLEQ
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                IDCR = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) < XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT (ICRS1)
                            XWRK1 = XDATT (IWRK1)
                         Else
                           If (ICRS2 < NORD) Then
                              IWRKT (ICRS1) = IWRKT (ICRS)
                              IWRK1 = IWRKT (ICRS1)
                              XWRK1 = XDATT (IWRK1)
                           Endif
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT (ICRS1)
                         XWRK1 = XDATT (IWRK1)
                      End Do
                   End If
                End Do
                ires_med = IWRK1
                Return
            Else
                ires_med = IMED7
                Return
            End If
         Else                       !      If (NLEQ > NMED)
!                                          Not enough high values
                XWRK1 = -XHUGE
                NORD = NLEQ - NEQU - NMED + 1
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
!
                      Do ICRS = IDON, IENDT (IDON1)
                         If (XDATT(IWRKT (ICRS)) > XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            IDCR = ICRS1 - 1
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT(ICRS1)
                            XWRK1 = XDATT(IWRK1)
                         Else
                            If (ICRS2 < NORD) Then
                               IWRKT (ICRS1) = IWRKT (ICRS)
                               IWRK1 = IWRKT(ICRS1)
                               XWRK1 = XDATT(IWRK1)
                            End If
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = IENDT (IDON1), IDON, -1
                         If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         IDCR = ICRS1 - 1
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT(ICRS1)
                         XWRK1 = XDATT(IWRK1)
                      End Do
                   Endif
                End Do
!
                ires_med = IWRK1
                Return
         End If
!
   END Subroutine r_med
Subroutine I_indmed (XDONT, INDM)
!  Returns index of median value of XDONT.
! __________________________________________________________
      Integer, Dimension (:), Intent (In) :: XDONT
      Integer, Intent (Out) :: INDM
! __________________________________________________________
      Integer :: IDON
!
      Allocate (IDONT (SIZE(XDONT)))
      Do IDON = 1, SIZE(XDONT)
         IDONT (IDON) = IDON
      End Do
!
      Call i_med (XDONT, IDONT, INDM)
!
      Deallocate (IDONT)
End Subroutine I_indmed
   Recursive Subroutine i_med (XDATT, IDATT, ires_med)
!  Finds the index of the median of XDONT using the recursive procedure
!  described in Knuth, The Art of Computer Programming,
!  vol. 3, 5.3.3 - This procedure is linear in time, and
!  does not require to be able to interpolate in the
!  set as the one used in INDNTH. It also has better worst
!  case behavior than INDNTH, but is about 30% slower in
!  average for random uniformly distributed values.
! __________________________________________________________
      Integer, Dimension (:), Intent (In) :: XDATT
      Integer, Dimension (:), Intent (In) :: IDATT
      Integer, Intent (Out) :: ires_med
! __________________________________________________________
!
      Integer, Parameter :: XHUGE = HUGE (XDATT)
      Integer :: XWRK, XWRK1, XMED7, XMAX, XMIN
!
      Integer, Dimension (7*(((Size (IDATT)+6)/7+6)/7)) :: ISTRT, IENDT, IMEDT
      Integer, Dimension (7*((Size(IDATT)+6)/7)) :: IWRKT
      Integer :: NTRI, NMED, NORD, NEQU, NLEQ, IMED, IDON, IDON1
      Integer :: IDEB, ITMP, IDCR, ICRS, ICRS1, ICRS2, IMAX, IMIN
      Integer :: IWRK, IWRK1, IMED1, IMED7, NDAT
!
      NDAT = Size (IDATT)
      NMED = (NDAT+1) / 2
      IWRKT = IDATT
!
!  If the number of values is small, then use insertion sort
!
     If (NDAT < 35) Then
!
!  Bring minimum to first location to save test in decreasing loop
!
         IDCR = NDAT
         If (XDATT (IWRKT (1)) < XDATT (IWRKT (IDCR))) Then
            IWRK = IWRKT (1)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT (1)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, NDAT - 2
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                XWRK = XWRK1
                IWRK = IWRK1
            Endif
         End Do
         IWRKT (1) = IWRK
!
! Sort the first half, until we have NMED sorted values
!
         Do ICRS = 3, NMED
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            IDCR = ICRS - 1
            Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
!  Insert any value less than the current median in the first half
!
         XWRK1 = XDATT (IWRKT (NMED))
         Do ICRS = NMED+1, NDAT
            XWRK = XDATT (IWRKT (ICRS))
            IWRK = IWRKT (ICRS)
            If (XWRK < XWRK1) Then
               IDCR = NMED - 1
               Do
                  If (XWRK >= XDATT (IWRKT(IDCR))) Exit
                  IWRKT (IDCR+1) = IWRKT (IDCR)
                  IDCR = IDCR - 1
               End Do
               IWRKT (IDCR+1) = IWRK
               XWRK1 = XDATT (IWRKT (NMED))
            End If
         End Do
         ires_med = IWRKT (NMED)
         Return
      End If
!
!  Make sorted subsets of 7 elements
!  This is done by a variant of insertion sort where a first
!  pass is used to bring the smallest element to the first position
!  decreasing disorder at the same time, so that we may remove
!  remove the loop test in the insertion loop.
!
      IMAX = 1
      IMIN = 1
      XMAX = XDATT (IWRKT(IMAX))
      XMIN = XDATT (IWRKT(IMIN))
      DO IDEB = 1, NDAT-6, 7
         IDCR = IDEB + 6
         If (XDATT (IWRKT(IDEB)) < XDATT (IWRKT(IDCR))) Then
            IWRK = IWRKT(IDEB)
         Else
            IWRK = IWRKT (IDCR)
            IWRKT (IDCR) = IWRKT(IDEB)
         Endif
         XWRK = XDATT (IWRK)
         Do ITMP = 1, 5
            IDCR = IDCR - 1
            IWRK1 = IWRKT (IDCR)
            XWRK1 = XDATT (IWRK1)
            If (XWRK1 < XWRK) Then
                IWRKT (IDCR) = IWRK
                IWRK = IWRK1
                XWRK = XWRK1
            Endif
         End Do
         IWRKT (IDEB) = IWRK
         If (XWRK < XMIN) Then
             IMIN = IWRK
             XMIN = XWRK
         End If
         Do ICRS = IDEB+1, IDEB+5
            IWRK = IWRKT (ICRS+1)
            XWRK = XDATT (IWRK)
            IDON = IWRKT(ICRS)
            If (XWRK < XDATT(IDON)) Then
               IWRKT (ICRS+1) = IDON
               IDCR = ICRS
               IWRK1 = IWRKT (IDCR-1)
               XWRK1 = XDATT (IWRK1)
               Do
                  If (XWRK >= XWRK1) Exit
                  IWRKT (IDCR) = IWRK1
                  IDCR = IDCR - 1
                  IWRK1 = IWRKT (IDCR-1)
                  XWRK1 = XDATT (IWRK1)
               End Do
               IWRKT (IDCR) = IWRK
            EndIf
         End Do
         If (XWRK > XMAX) Then
             IMAX = IWRK
             XMAX = XWRK
         End If
      End Do
!
!  Add-up alternatively MAX and MIN values to make the number of data
!  an exact multiple of 7.
!
      IDEB = 7 * (NDAT/7)
      NTRI = NDAT
      If (IDEB < NDAT) Then
!
         Do ICRS = IDEB+1, NDAT
            XWRK1 = XDATT (IWRKT (ICRS))
            IF (XWRK1 > XMAX) Then
               IMAX = IWRKT (ICRS)
               XMAX = XWRK1
            End If
            IF (XWRK1 < XMIN) Then
               IMIN = IWRKT (ICRS)
               XMIN = XWRK1
            End If
         End Do
         IWRK1 = IMAX
         Do ICRS = NDAT+1, IDEB+7
               IWRKT (ICRS) = IWRK1
               If (IWRK1 == IMAX) Then
                  IWRK1 = IMIN
               Else
                  NMED = NMED + 1
                  IWRK1 = IMAX
               End If
         End Do
!
         Do ICRS = IDEB+2, IDEB+7
            IWRK = IWRKT (ICRS)
            XWRK = XDATT (IWRK)
            Do IDCR = ICRS - 1, IDEB+1, - 1
               If (XWRK >= XDATT (IWRKT(IDCR))) Exit
               IWRKT (IDCR+1) = IWRKT (IDCR)
            End Do
            IWRKT (IDCR+1) = IWRK
         End Do
!
         NTRI = IDEB+7
      End If
!
!  Make the set of the indices of median values of each sorted subset
!
         IDON1 = 0
         Do IDON = 1, NTRI, 7
            IDON1 = IDON1 + 1
            IMEDT (IDON1) = IWRKT (IDON + 3)
         End Do
!
!  Find XMED7, the median of the medians
!
         Call i_med (XDATT, IMEDT(1:IDON1), IMED7)
         XMED7 = XDATT (IMED7)
!
!  Count how many values are not higher than (and how many equal to) XMED7
!  This number is at least 4 * 1/2 * (N/7) : 4 values in each of the
!  subsets where the median is lower than the median of medians. For similar
!  reasons, we also have at least 2N/7 values not lower than XMED7. At the
!  same time, we find in each subset the index of the last value < XMED7,
!  and that of the first > XMED7. These indices will be used to restrict the
!  search for the median as the Kth element in the subset (> or <) where
!  we know it to be.
!
         IDON1 = 1
         NLEQ = 0
         NEQU = 0
         Do IDON = 1, NTRI, 7
            IMED = IDON+3
            If (XDATT (IWRKT (IMED)) > XMED7) Then
                  IMED = IMED - 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                  IMED = IMED + 2
                  If (XDATT (IWRKT (IMED)) > XMED7) Then
                     IMED = IMED - 1
                  Else If (XDATT (IWRKT (IMED)) < XMED7) Then
                     IMED = IMED + 1
                  Endif
            Endif
            If (XDATT (IWRKT (IMED)) > XMED7) Then
               NLEQ = NLEQ + IMED - IDON
               IENDT (IDON1) = IMED - 1
               ISTRT (IDON1) = IMED
            Else If (XDATT (IWRKT (IMED)) < XMED7) Then
               NLEQ = NLEQ + IMED - IDON + 1
               IENDT (IDON1) = IMED
               ISTRT (IDON1) = IMED + 1
            Else                    !       If (XDATT (IWRKT (IMED)) == XMED7)
               NLEQ = NLEQ + IMED - IDON + 1
               NEQU = NEQU + 1
               IENDT (IDON1) = IMED - 1
               Do IMED1 = IMED - 1, IDON, -1
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     IENDT (IDON1) = IMED1 - 1
                  Else
                     Exit
                  End If
               End Do
               ISTRT (IDON1) = IMED + 1
               Do IMED1 = IMED + 1, IDON + 6
                  If (XDATT (IWRKT (IMED1)) == XMED7) Then
                     NEQU = NEQU + 1
                     NLEQ = NLEQ + 1
                     ISTRT (IDON1) = IMED1 + 1
                  Else
                     Exit
                  End If
               End Do
            Endif
            IDON1 = IDON1 + 1
         End Do
!
!  Carry out a partial insertion sort to find the Kth smallest of the
!  large values, or the Kth largest of the small values, according to
!  what is needed.
!
!
         If (NLEQ - NEQU + 1 <= NMED) Then
            If (NLEQ < NMED) Then   !      Not enough low values
                IWRK1 = IMAX
                XWRK1 = XDATT (IWRK1)
                NORD = NMED - NLEQ
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                IDCR = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) < XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT (ICRS1)
                            XWRK1 = XDATT (IWRK1)
                         Else
                           If (ICRS2 < NORD) Then
                              IWRKT (ICRS1) = IWRKT (ICRS)
                              IWRK1 = IWRKT (ICRS1)
                              XWRK1 = XDATT (IWRK1)
                           Endif
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = ISTRT (IDON1), IDON + 6
                         If (XDATT (IWRKT (ICRS)) >= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK >= XDATT (IWRKT (IDCR))) Exit
                               IWRKT  (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT (ICRS1)
                         XWRK1 = XDATT (IWRK1)
                      End Do
                   End If
                End Do
                ires_med = IWRK1
                Return
            Else
                ires_med = IMED7
                Return
            End If
         Else                       !      If (NLEQ > NMED)
!                                          Not enough high values
                XWRK1 = -XHUGE
                NORD = NLEQ - NEQU - NMED + 1
                IDON1 = 0
                ICRS1 = 1
                ICRS2 = 0
                Do IDON = 1, NTRI, 7
                   IDON1 = IDON1 + 1
                   If (ICRS2 < NORD) Then
!
                      Do ICRS = IDON, IENDT (IDON1)
                         If (XDATT(IWRKT (ICRS)) > XWRK1) Then
                            IWRK = IWRKT (ICRS)
                            XWRK = XDATT (IWRK)
                            IDCR = ICRS1 - 1
                            Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                            End Do
                            IWRKT (IDCR+1) = IWRK
                            IWRK1 = IWRKT(ICRS1)
                            XWRK1 = XDATT(IWRK1)
                         Else
                            If (ICRS2 < NORD) Then
                               IWRKT (ICRS1) = IWRKT (ICRS)
                               IWRK1 = IWRKT(ICRS1)
                               XWRK1 = XDATT(IWRK1)
                            End If
                         End If
                         ICRS1 = MIN (NORD, ICRS1 + 1)
                         ICRS2 = MIN (NORD, ICRS2 + 1)
                      End Do
                   Else
                      Do ICRS = IENDT (IDON1), IDON, -1
                         If (XDATT(IWRKT (ICRS)) <= XWRK1) Exit
                         IWRK = IWRKT (ICRS)
                         XWRK = XDATT (IWRK)
                         IDCR = ICRS1 - 1
                         Do IDCR = ICRS1 - 1, 1, - 1
                               If (XWRK <= XDATT(IWRKT (IDCR))) Exit
                               IWRKT (IDCR+1) = IWRKT (IDCR)
                         End Do
                         IWRKT (IDCR+1) = IWRK
                         IWRK1 = IWRKT(ICRS1)
                         XWRK1 = XDATT(IWRK1)
                      End Do
                   Endif
                End Do
!
                ires_med = IWRK1
                Return
         End If
!
   END Subroutine i_med
end module m_indmed
