!***************************************************************
!*   Calculate the median value of an array without sorting    *
!* ----------------------------------------------------------- *
!* REFERENCE:                                                  *
!*      "NUMERICAL RECIPES By W.H. Press, B.P. Flannery,       *
!*       S.A. Teukolsky and W.T. Vetterling, Cambridge         *
!*       University Press, 1986" [BIBLI 08].                   *
!* ----------------------------------------------------------- *
!* SAMPLE RUN:                                                 *
!*                                                             *
!* Given table:                                                *
!*                                                             *
!* 407.8 192.8 851.1 604.4 932.3 799.4 914.5 965.8 453.7 295.1 *
!* 154.5 977.4 410.2 916.2 934.7 504.8 823.2 225.2 456.6  49.0 *
!* 933.5 663.0 335.3 346.6 568.7 956.1 654.7 300.7 379.6 591.9 *
!* 992.9 689.6 644.7 305.4 148.2 257.2 664.6 612.1 713.0  99.7 *
!*  46.5 167.6 984.6 847.2  55.4  82.7 999.0  10.7 877.7 929.4 *
!* 398.1 972.8 874.1 755.1 472.1 122.8 671.4  35.5 128.8  76.8 *
!* 454.2 959.2 510.1 791.3 122.8 176.6 237.9 995.8 548.3 309.8 *
!* 162.6 996.5 750.0 250.6 577.7 761.1 101.9 797.1 539.0 723.5 *
!*                                                             *
!* Median value:    558.5251                                   *
!*                                                             *
!*                            F90 Release By J-P Moreau, Paris.*
!*                                   (www.jpmoreau.fr)         *
!***************************************************************
PROGRAM TMDIAN1
real    A(100)     !Table to be sorted
real    MAX_VALUE  !Maximum value of table
integer Count(1)   !Current value of system clock

  !initialize random number generator
  call System_Clock(Count(1))
  call Random_Seed(Put=Count)

  N=80  !initialize size of table
  MAX_VALUE = 1000

  !generate random table of numbers (from 0 to 1000)
  do i=1, N
    call Random_Number(x)  !returns random x >= 0 and <1  
    A(i)=MAX_VALUE*x
  end do

  print *,'Given table:'
  call TWRIT(N,A)
 
  !call MDIAN1 subroutine
  call MDIAN1(A,N,AMED)

  print *,' '
  print *,'Median value: ', AMED
  print *,' '
  stop

END


!*******************************************************
!* Given an array X of N numbers, returns their median *
!* value XMED. The array X is not modified, and is     *
!* accessed sequentially in each consecutive pass.     *
!*******************************************************
SUBROUTINE MDIAN1(X,N,XMED)
  real X(N)
  real, parameter :: BIG = 1.e30, AFAC=1.5, AMP=1.5
! Here, AMP is an overconvergence factor: on each iteration,
! we move the guess by this factor. AFAC is a factor used to
! optimize the size of the "smoothing constant" EPS at each
! iteration.
  A=0.5*(X(1)+X(N))
  EPS=ABS(X(N)-X(1))
  AP=BIG
  AM=-BIG
1 SUM=0.
  SUMX=0.
  NP=0
  NM=0
  XP=BIG
  XM=-BIG
  do J=1, N
    XX=X(J)
    if(XX.ne.A)then
      if(XX.gt.A)then
        NP=NP+1
	if(XX.lt.XP) XP=XX
      else if(XX.lt.A)then
        NM=NM+1
	if(XX.gt.XM) XM=XX
      endif
      DUM=1./(EPS+ABS(XX-A))
      SUM=SUM+DUM
      SUMX=SUMX+XX*DUM
    endif
  end do
  print *,' NP=',NP,'  NM=',NM
  if(NP-NM.ge.2)then      !guess is too low, make another pass
    AM=A
    AA=XP+MAX(0.,SUMX/SUM-A)*AMP
    if(AA.gt.AP) AA=0.5*(A+AP)
    EPS=AFAC*ABS(AA-A)
    A=AA
    goto 1
  else if(NM-NP.ge.2)then !guess is too high
    AM=A
    AA=XM+MIN(0.,SUMX/SUM-A)*AMP
    if(AA.lt.AM) AA=0.5*(A+AM)
    EPS=AFAC*ABS(AA-A)
    A=AA
    goto 1
  else                    !guess is ok
    if(MOD(N,2).eq.0)then !for even N median is an average
      if(NP.eq.NM)then
        XMED=0.5*(XP+XM)
      else if(NP.gt.NM)then
        XMED=0.5*(A+XP)
      else
        XMED=0.5*(XM+A)
      endif
    else
      if(NP.eq.NM)then
        XMED=A
      else if(NP.gt.NM)then
        XMED=XP
      else
        XMED=XM
      endif
    endif
  endif
  return
END

!write table of size N to standard output
SUBROUTINE TWRIT(N,ARR)
real ARR(N)
  print *,' '
  WRITE(*,10) (ARR(I),I=1,N)
  return
10 FORMAT(10F6.1)
END

!end of file tmdian1.f90