!---------------- 滑动t检验 ---------------------------!
!要求输入的有:                                         !
!  X(n):需作检验的序列;                                !
!  n1:滑动长度;                                        !
!  Alfa:信度;                                          !
!  Year:原始资料的第一个年份(如果资料从1950--1990,     !
!       则Year=1950);注意Year为整数                    !
!输出的有:                                             !
!  第一列:检验的年份;                                  !
!  第二列:各个滑动的t检验值;                           !
!  第三列:信度Alfa下的t值;                             !
!  第四列:各个t检验值所能达到的信度                    !
!注意:                                                 !
!  本程序引用了fortran 5.0 内部sstats.lib中的统计函数  !
!---------------------------------- 程正泉 2000.5 -----!
subroutine sub_SmoothT_test(n,n1,X,Alfa,Year,filename)
implicit none
real,intent(in)::Alfa
integer,intent(in)::n,n1,Year
real,dimension(n),intent(in)::X
real,allocatable,dimension(:)::X1,X2
real::t,df,tin,tdf
integer::i,j,Iyear
character*50 filename

open(90,file=filename)

allocate(X1(n1),X2(n1)) 
do i=1,n-2*n1+1
   do j=1,n1
      X1(j)=X(j+i-1)
      X2(j)=X(j+i+n1-1)
   enddo
   call t_test(n1,n1,X1,X2,t)
   write(90,'(i5,3f7.3,f10.5)')i+Year+n1-1,t,tin(1-Alfa/2.0,real(n1+n1-2)),  &
         -1.0*tin(1-Alfa/2.0,real(n1+n1-2)),(1-tdf(abs(t),real(n1+n1-2)))*2.0
enddo
close(90)
deallocate(X1,X2)
end

!----------------------------------------
subroutine t_test(n1,n2,X1,X2,t)
implicit none
integer,intent(in)::n1,n2
real,dimension(n1),intent(in)::X1,X2
real,intent(out)::t
real::ave1,ave2,var1,var2,Sp

call avevar(n1,X1,ave1,var1)
call avevar(n2,X2,ave2,var2)
Sp=((N1-1)*VAR1+(N2-1)*VAR2)/(N1+N2-2)
Sp=SQRT(Sp)
t=(ave2-ave1)/(Sp*(sqrt(1.0/n1+1.0/n2)))

end

!--------------------------------
subroutine AveVar(n,X,ave,var)
implicit none
integer,intent(in)::n
real,dimension(n)::X(n)
real,intent(out)::ave,var
integer::i

ave=0.0
do i=1,n
   ave=ave+X(i)/float(n)
enddo
var=0.0
do i=1,n
   var=var+(X(i)-ave)*(X(i)-ave)/float(n)
enddo
end