如何用fortran编写一个线性插值程序

这是一元全区间等距插值子程序，X1和H为实数，分别为等距节点中的第一个节点值和等距节点的步长；N为整数，等距节点的个数；Y(N)，存放N个等距节点上的函数值；T是指定插值点的值；Z返回指定插值点T处的函数近似值。
SUBROUTINE EELGR(X1,H,N,Y,T,Z)
DIMENSION Y(N)
DOUBLE PRECISION X1,H,Y,T,Z,S,XX,P,Q
Z=0.0
IF (N.LE.0) RETURN
IF (N.EQ.1) THEN
Z=Y(1)
RETURN
END IF
IF (N.EQ.2) THEN
Z=(Y(2)*(T-X1)-Y(1)*(T-X1-H))/H
RETURN
END IF
IF (T.GT.X1) THEN
P=(T-X1)/H
I=P
Q=1.0*I
I=I+1
IF (P.GT.Q) I=I+1
ELSE
I=1
END IF
K=I-4
IF (K.LT.1) K=1
M=I+3
IF (M.GT.N) M=N
Z=0.0
DO 30 I=K,M
S=1.0
XX=X1+(I-1)*H
DO 20 J=K,M
IF (J.NE.I) THEN
XJ=X1+(J-1)*H
S=S*(T-XJ)/(XX-XJ)
END IF
20 CONTINUE
Z=Z+S*Y(I)
30 CONTINUE
RETURN
END
这是测试用的主程序代码。函数为F(x) = e^(-x)，求x = 0.63处的函数近似值：
PROGRAM EELGRO
DIMENSION Y(10)
DOUBLE PRECISION Y,H,X1,T,Z
DATA Y/0.904837,0.818731,0.740818,0.670320,0.606531,
* 0.548812,0.496585,0.449329,0.406570,0.367879/
H=0.1
N=10
X1=0.1
T=0.63
CALL EELGR(X1,H,N,Y,T,Z)
WRITE(*,20) T,Z
20 FORMAT(1X,'T=',F6.3,10X,'Z=',D15.6)
END