FORTRAN 第一类边界条件三次样条插值程序(徐士良版)问题！！！

[feihong2201]

在应用过程中发现插值得到的函数一阶导数、二阶导数明显不对，而用书中例子则没问题，求解答？？？
代码如下？
主程序：
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER(ND=3969,ND1=422)
DOUBLE PRECISION X,Y,XX,DY,DDY,S,DS,DDS,T,DY1,DYN,H !程序中用到的变量
DIMENSION X(ND1),XS(ND1),YS(ND1),ZS(ND1),XV(ND1),YV(ND1),ZV(ND1),
*XX(ND),DY(ND1),DDY(ND1),H(ND1),DS(ND),DDS(ND),S1(ND)

OPEN(UNIT=1,FILE='time.DAT')
DO 11 I=1,ND
READ(1,*) XX(I)
IF(I.GT.1.AND.XX(I).LT.XX(I-1)) THEN
XX(I)=XX(I)+24.D0
ENDIF
11 CONTINUE
CLOSE(1)
OPEN(UNIT=3,FILE='xyz_trans1_std')
DO 23 I=1,ND1
READ(3,*) X1,X(I),X3,XS(I),YS(I),ZS(I),XV(I),YV(I),ZV(I),X1
*0,X11,X12
X(I)=AINT(X(I)/100.D0)+(X(I)-AINT(X(I)/100.D0)*100.D0)/60.D0
IF(I.GT.1.AND.X(I).LT.X(I-1)) X(I)=X(I)+24.D0
23 CONTINUE
CLOSE(1)

N=ND1
M=ND

C 计算XS插值
DYXS1=XV(1)
DYXSN=XV(ND1)
CALL ESPL1(X,XS,N,DYXS1,DYXSN,XX,M,DY,DDY,S1,DS,DDS,T,H)

OPEN(UNIT=4,FILE='STK.VEC')
WRITE(4,10)
10 FORMAT(6X,'X(I)',10X,'Y(I)',10X,'DY(I)',9X,'DDY(I)')
WRITE(4,20) (X(I),XS(I),DY(I),DDY(I),I=1,N)
20 FORMAT(1X,F14.8,2X,F17.6,2X,F17.6,2X,F18.6)
WRITE(*,*)
WRITE(4,30) T
30 FORMAT(1X,'T=',D15.9)
WRITE(*,*)
WRITE(4,40)
40 FORMAT(6X,'XX(I)',9X,' XS(I)',9X,'DS(I)',9X,'DDS(I)')
WRITE(4,20) (XX(I),S1(I),DS(I),DDS(I),I=1,M)
END
插值程序：

SUBROUTINE ESPL1(X,Y,N,DY1,DYN,XX,M,DY,DDY,S,DS,DDS,T,H)
DIMENSION X(N),Y(N),XX(M),DY(N),DDY(N)
DIMENSION S(M),DS(M),DDS(M),H(N)
DOUBLE PRECISION X,Y,XX,DY,DDY,S,DS,DDS,H,DY1,DYN,
* T,H0,H1,BETA,ALPHA
DY(1)=0.0
H(1)=DY1
H0=X(2)-X(1)
DO 10 J=2,N-1
H1=X(J+1)-X(J)
ALPHA=H0/(H0+H1)
BETA=(1.0-ALPHA)*(Y(J)-Y(J-1))/H0
BETA=3.0*(BETA+ALPHA*(Y(J+1)-Y(J))/H1)
DY(J)=-ALPHA/(2.0+(1.0-ALPHA)*DY(J-1))
H(J)=(BETA-(1.0-ALPHA)*H(J-1))
H(J)=H(J)/(2.0+(1.0-ALPHA)*DY(J-1))
H0=H1
10 CONTINUE
DY(N)=DYN
DO 20 J=N-1,1,-1
20 DY(J)=DY(J)*DY(J+1)+H(J)
DO 30 J=1,N-1
30 H(J)=X(J+1)-X(J)
DO 40 J=1,N-1
H1=H(J)*H(J)
DDY(J)=6.0*(Y(J+1)-Y(J))/H1-
* 2.0*(2.0*DY(J)+DY(J+1))/H(J)
40 CONTINUE
H1=H(N-1)*H(N-1)
DDY(N)=6.0*(Y(N-1)-Y(N))/H1+
* 2.0*(2.0*DY(N)+DY(N-1))/H(N-1)
T=0.0
DO 50 I=1,N-1
H1=0.5*H(I)*(Y(I)+Y(I+1))
H1=H1-H(I)*H(I)*H(I)*(DDY(I)+DDY(I+1))/24.0
T=T+H1
50 CONTINUE
DO 70 J=1,M
IF (XX(J).GE.X(N)) THEN
I=N-1
ELSE
I=1
60 IF (XX(J).GT.X(I+1)) THEN
I=I+1
GOTO 60
END IF
END IF
H1=(X(I+1)-XX(J))/H(I)
S(J)=(3.0*H1*H1-2.0*H1*H1*H1)*Y(I)
S(J)=S(J)+H(I)*(H1*H1-H1*H1*H1)*DY(I)
DS(J)=6.0*(H1*H1-H1)*Y(I)/H(I)
DS(J)=DS(J)+(3.0*H1*H1-2.0*H1)*DY(I)
DDS(J)=(6.0-12.0*H1)*Y(I)/(H(I)*H(I))
DDS(J)=DDS(J)+(2.0-6.0*H1)*DY(I)/H(I)
H1=(XX(J)-X(I))/H(I)
S(J)=S(J)+(3.0*H1*H1-2.0*H1*H1*H1)*Y(I+1)
S(J)=S(J)-H(I)*(H1*H1-H1*H1*H1)*DY(I+1)
DS(J)=DS(J)-6.0*(H1*H1-H1)*Y(I+1)/H(I)
DS(J)=DS(J)+(3.0*H1*H1-2.0*H1)*DY(I+1)
DDS(J)=DDS(J)+(6.0-12.0*H1)*Y(I+1)/(H(I)*H(I))
DDS(J)=DDS(J)-(2.0-6.0*H1)*DY(I+1)/H(I)
70 CONTINUE
RETURN
END
结果：
X(I) Y(I) DY(I) DDY(I)
21.60000000 -118481328.330000 -1970.745800 -1473931617.829629
21.61666667 -118599536.810000 -8990821.426489 395269536.146876
21.63333333 -118717672.850000 -6576757.148245 -105581822.757609
21.65000000 -118835736.550000 -7218103.180531 28620298.883312
21.66666667 -118953728.040000 -7040764.329631 -7339636.775252
21.68333333 -119071647.440000 -7082799.700945 2295392.217546
21.70000000 -119189494.850000 -7066062.666591 -286948.095049
21.71666667 -119307270.400000 -7065082.432697 404576.162314
21.73333333 -119424974.190000 -7059888.802621 218659.446773
21.75000000 -119542606.350000 -7055833.356815 267994.049950
21.76666667 -119660166.990000 -7051481.770121 254196.353337
。。。。。。
。。。。。。
结果从第三行开始就明显不对，一阶导数和二阶导数都很离谱。。。