FORTRAN程序中MK检验的一点小问题

suifeng1106

这几天在家园里从http://bbs.06climate.com/forum.php?mod=viewthread&tid=633下了MK检验的程序，程序如下：

                               
登录/注册后可看大图

C        THIS IS A PROGRAM FOR DETECTING ABRUPT CLIMATIC CHANGE
C        BY USING MANN-KENDALL TECHNIQUE
        PROGRAM MK
        DIMENSION Y(2000),YY(2000),U(2000),UF(2000),UB(2000),
     &        M(2000),MD(2000)
        WRITE(*,10)
  10        FORMAT(2X,'N=1944,NYEAR=1')
        READ(*,*)N,NYEAR
C        ***************************************************
C        * N:     SAMPLE SIZE                              *
C        * NYEAR: FIRST YEAR OF THE TIME SERIES            *
C        * Y(N):  ORIGINAL TIME SERIES                     *
C        * UF(N): ORIGINAL SERIES OF U(LN)                 *
C        * UB(N): COUNTER SERIES OF U(LN)                  *
C        * A,B:   CRITICAL VALUE 1.96 AND -1.96            *
C        ***************************************************          
        OPEN(2,FILE='pc2.txt')
        READ(2,*)(Y(I),I=1,N)
        CALL SMK(Y,M,MD,UF,N)
        DO 20 I=1,N
  20        YY(I)=Y(N+1-I)
        CALL SMK(YY,M,MD,U,N)
        DO 30 I=1,N
  30        UB(I)=-U(N+1-I)
        OPEN(3,FILE='test1.txt',STATUS='NEW')
        A=1.96
        B=-1.96
        DO 40 I=1,N
        WRITE(3,50)NYEAR+I-1,UF(I),UB(I),A,B
  50        FORMAT(1X,I4,4F8.2)
  40        CONTINUE
        CLOSE(3)
        STOP
        END
C***********************************************************
        SUBROUTINE SMK(Y,M,MD,U,N)
        DIMENSION Y(N),M(N),MD(N),U(N)
        M(1)=0
        DO 10 I=2,N
        M(I)=0
        MD(I)=0
        DO 20 J=1,I-1
        IF(Y(I).LT.Y(J))GOTO 20
        M(I)=M(I)+1
  20        CONTINUE
        MD(I)=MD(I-1)+M(I)
  10        CONTINUE
        U(1)=0.0
        DO 30 I=2,N
        E=I*(I-1)/4.00
        VAR=I*(I-1)*(2*I+5)/72.00
        U(I)=(MD(I)-E)/SQRT(VAR)
!        U(I)=(MD(I)-E)/SQRT(abs(VAR))
  30        CONTINUE
        RETURN
        END


                               
登录/注册后可看大图

显示SQRT出错，可能是VAR是一个负值，然后我把 U(I)=(MD(I)-E)/SQRT(VAR) 改成 U(I)=(MD(I)-E)/SQRT(abs(VAR))，程序可以成功运行，出来的数据也看不出来有什么奇怪，我想问的问题是加ABS（）指令生成的结果会出现错误吗？

suifeng1106

请教一下大家