!------------------------------------------- !Include "Betai.f90" !----------------------------------------------------------------------- Subroutine Linear_Trend(X,Y,N,NumYear,A1,F_F10,RMiss) ! Tinh he so goc cua duong hoi qui TT Y=A0+A1*X va Gia tri Fisher (Alfa=10%) ! Input: Cac mang X(N), Y(N), N, RMiss (=-99) ! NumYear: So nam du dai ! Output: A1, F_F10 Integer :: N, i, j, N1, NumYear Real :: X(N), Y(N), W1(N), W2(N), TBx, TBy, R, A0, A1, F10, F, U, Q, Yhq, F_F10 Real :: Ws1(N), Ws2(N) Real :: Sx, Sy, RR, RMiss !print*,Y TBx = 0; TBy = 0; RR = 0; N1 = 0; Do i=1,N if (Y(i) /= RMiss) then N1 = N1 + 1 TBx = TBx + X(i); TBy = TBy + Y(i) Sx = Sx + X(i)*X(i); Sy = Sy + Y(i)*Y(i) RR = RR + X(i)*Y(i) W1(N1) = X(i); W2(N1) = Y(i) endif Enddo if ((N1 < NumYear).OR.((TBx == 0).AND.(TBy == 0)).OR.(Sx==0).OR.(Sy==0)) then A0 = -99; A1 = -99; R = -99; F_F10 = -99; else TBx = TBx/Real(N1); TBy = TBy/Real(N1) Sx = SQRT(Sx/Real(N1)-TBx*TBx); Sy = SQRT(Sy/Real(N1)-TBy*TBy) R = (RR/Real(N1)-TBx*TBy)/(Sx*Sy) A1 = R*Sy/Sx A0 = TBy - A1*TBx F10 = FINV(0.10,1.0,Real(N1)) U = 0; Q = 0 Do i=1,N1 Yhq = A0 + A1*W1(i) U = U + (Yhq-TBy)*(Yhq-TBy) Q = Q + (W2(i)-Yhq)*(W2(i)-Yhq) Enddo if (Q == 0) then F = -99 else F = U/(Q/(N1-2)) endif if (F /=-99) then; F_F10 = F-F10; else; F_F10 = -99; endif endif End subroutine Linear_Trend !------------------------------------------- Real FUNCTION FINV(P,N1,N2) Real, PARAMETER :: EPS=1.0E-6 REAL :: P,AP,N1,N2, A,B,C,P0,SS,SS1 IF (P.LT.0.0.OR.P.GT.1.0) THEN WRITE(*,*)' INVALID NUMERIC INPUT IN FINV FUNCTION' STOP ENDIF IF (P.EQ.0.0) THEN WRITE(*,*)' THE P VALUE TOO SMALL IN FINV FUNCTION' STOP ELSE IF (P.EQ.1) THEN FINV=0.0 RETURN ENDIF ! A=0.0 B=9999999.0 C=A AP=P 10 P0=FDIST(B,N1,N2) SS=ABS((P0-AP)/AP) SS1=ABS((B-C)/B) IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN IF (AP.GT.P0) THEN C=B B=(A+B)/2.0 GOTO 10 ELSE IF (AP.LT.P0) THEN B=(C+B)/2.0 GOTO 10 ENDIF ENDIF A=(C+B)/2.0 FINV=A END FUNCTION FINV !------------------------------------------- Real FUNCTION FDIST(X0,N1,N2) REAL :: X0,N1,N2,X X=N2/(N2+N1*X0) FDIST=BETAI(0.5*N2,0.5*N1,X) END FUNCTION FDIST !------------------------------------------- Real FUNCTION betai(a,b,x) REAL :: a,b,x REAL :: bt if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai' if(x.eq.0..or.x.eq.1.)then bt=0. else bt=exp(gammln_1(a+b)-gammln_1(a)-gammln_1(b)+a*log(x)+b*log(1.-x)) endif if(x.lt.(a+1.)/(a+b+2.))then betai=bt*betacf(a,b,x)/a return else betai=1.-bt*betacf(b,a,1.-x)/b return endif END FUNCTION betai !------------------------------------------- Real FUNCTION betacf(a,b,x) INTEGER, PARAMETER :: MAXIT=100 REAL :: a,b,x Real, PARAMETER :: EPS=3.e-7,FPMIN=1.e-30 INTEGER :: m,m2 REAL :: aa,c,d,del,h,qab,qam,qap qab=a+b qap=a+1. qam=a-1. c=1. d=1.-qab*x/qap if(abs(d).lt.FPMIN)d=FPMIN d=1./d h=d do 11 m=1,MAXIT m2=2*m aa=m*(b-m)*x/((qam+m2)*(a+m2)) d=1.+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d h=h*d*c aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) d=1.+aa*d if(abs(d).lt.FPMIN)d=FPMIN c=1.+aa/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue pause 'a or b too big, or MAXIT too small in betacf' 1 betacf=h END FUNCTION betacf !------------------------------------------- Real FUNCTION gammln_1(xx) REAL :: xx INTEGER :: j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & -.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln_1=tmp+log(stp*ser/x) END FUNCTION gammln_1 !----------------------------------------------------------------------- Subroutine MK_Sen(Y,N,NumYear,Miss,SenSlope,prob,N1) ! Tinh he so goc Sen (SenSlope) cua chuoi thoi gian Y(N) va Xac suat (muc y nghia - prob) ! Input: Chuoi thoi gian Y(N), Miss (=-99) ! Output: He so goc Sen (SenSlope), Muc y nghia (prob) va So phan tu (thuc) cua X(N) Integer :: N,i,j,S, Ng, T(N), NumYear, 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 ! print*,i," ",X(i) 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) Call sort(Ng, Q) 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) Call sort(N, X) 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 !----------------------------------------------------------------------- Real function percentile(x, length, per, RMiss) integer length real x(length), per, RMiss real xtos(length),bb,cc,oout integer nn, i nn=0 do i=1, length if(x(i) /= RMiss)then nn=nn+1 xtos(nn)=x(i) endif enddo if(nn.eq.0) then oout=RMiss else call sort(nn,xtos) bb=nn*per+per/3.+1/3. cc=real(int(bb)) if(int(cc).ge.nn) then oout=xtos(nn) else oout=xtos(int(cc))+(bb-cc)*(xtos(int(cc)+1)-xtos(int(cc))) endif endif percentile = oout end function percentile !----------------------------------------------------------------------- SUBROUTINE 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 sort !--------------------------End of this file--------------------