!-----------------------------------------------------------------------
Subroutine MK_Sen(Y,N,Miss,SenSlope,prob,N1)
Integer :: N,i,j,S, Ng, T(N), N1
!Real, Intent (In) :: X(N)
Real :: X(N), Y(N)
Real :: VarS, z, Miss, prob, SenSlope, Q(N*(N-1)/2)
Logical :: Empty
X(:) = Y(:)   ! Dam bao mang dau vao khong bi thay doi
Empty = .true.; N1=0
Do i=1, N
   if (X(i) /= Miss) then
      N1 = N1+1
   endif
Enddo
if (N1 >= NumYear) Empty = .false.
if (Empty) then
   SenSlope = Miss; prob = Miss ; Return
else  ! Calculate Sen's slope first
  Ng = 0
  Do i=1, N-1
     Do j=i+1, N
       if (X(j) /= Miss .And. X(i) /= Miss) then
        Ng = Ng+1
        Q(Ng) = (X(j)-X(i))/Real(j-i)
       endif
     Enddo
  Enddo
  Call My_Sort (Q, Ng)
  if (mod(Ng,2) /= 0) then
     SenSlope = Q(Ng/2+1)
  else
     SenSlope = (Q(Ng/2) + Q(Ng/2+1))/2.
  endif
  S = 0
  Do i=1,N-1
     Do j=i+1,N
       if (X(j) /= Miss .And. X(i) /= Miss) then
        if (X(j) > X(i)) then
           S = S+1
        else if (X(j) < X(i)) then
           S = S-1
        endif
       endif
     Enddo
  Enddo
  Call My_Sort(X, N)
  Ng = 0; T(:) = 0; i=0
  Do
     if (i>=N) exit
     i = i+1
     if (i == 1) then
        Ng = Ng+1; T(Ng) = T(Ng)+1
     else
        if (X(i)==X(i-1)) then
           T(Ng) = T(Ng)+1
        else
           Ng = Ng+1; T(Ng) = T(Ng)+1
        endif
     endif
  Enddo
  VarS = 0
  if (X(1) == Miss) then ! Loai nhom dau tien neu co Miss
     Do i=2, Ng
        VarS = VarS + T(i)*(T(i)-1)*(2*T(i)+5)
     Enddo
  else
     Do i=1, Ng
        VarS = VarS + T(i)*(T(i)-1)*(2*T(i)+5)
     Enddo
  endif
  VarS = (N*(N-1)*(2*N+5)-VarS)/18.
  if (S > 0) then
     z = (S-1)/SQRT(VarS)
  else if (S < 0) then
     z = (S+1)/SQRT(VarS)
  else
     z = 0
  endif
  prob = erfcc(abs(z)/1.4142136)
  if (prob>=1.) then !print "(2(A,F10.3))","S=",SenSlope," Pro=",prob
     SenSlope = Miss; prob = Miss
  endif
endif
End Subroutine MK_Sen
Subroutine My_Sort(X,N)
! Sap xep X(n) thanh chuoi tang dan
Integer :: N, i,j
Real :: X(N), tmp

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
End Subroutine My_Sort
Real FUNCTION erfcc (x)
REAL :: x, tmp
REAL :: t,z
z = abs(x)
t = 1./(1.+0.5*z)
 tmp  = t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+       &
        t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ &
        t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
if (x.lt.0.) tmp=2.-tmp
erfcc=tmp
END FUNCTION erfcc
!-----------------------------------------------------------------------