!----------------------------------------------
 Real Function percentile(Y,N1,p,Miss)
! Tham khao:
! http://www.childrens-mercy.org/stats/definitions/percentiles.htm
  Integer :: N, N1, i, j
  Real :: Y(n1), X(n1), Miss, tmp, p, per  ! (p value mustbe in [0-1])

! Filter and store Y into X:
  N=0
  Do i=1, N1
     if (Y(i) /= Miss) then
        N=N+1; X(N) = Y(i)
     endif
  Enddo
  !print*,"percentile: ",N1,N; pause
  if (N<5) then
     per=Miss
  else
     ! Work on X:
     !Do i = 1, n-1
     !   Do j = i+1, n
     !      if (x(j) < x(i)) then
     !         tmp = x(j)
     !         x(j) = x(i)
     !         x(i) = tmp
     !      endif
     !   Enddo
     !Enddo
     ! Sort X:
     Call T_sort(n,x)
     tmp = p*(n+1)
     if (tmp < 1) then
        per = x(1)
     elseif (tmp > n) then
        per = x(n)
     else
        i = int(tmp)
        j = i+1
        per = x(i) + (x(j)-x(i))*(tmp-real(i))
     endif
  endif 
  percentile = per
  End Function Percentile
!-----------------------------------------------------------------------
      SUBROUTINE T_sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=0
2         arr(i+1)=a
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END SUBROUTINE T_sort
!----------------------------------------------