fortran程序错误请教

a2883524

c
c     TRAN SOLVES THE LINEAR TRANSPORT EQUATION USING
C     VARIOUS EXPLICIT AND IMPLICIT SCHEMES
C
      DIMENSION R(65),T(65),TD(65),TEX(65),A(5,65)
     1,D(4),X(65)
      OPEN(1,FILE=';TRAN.DAT';)
OPEN(6,FILE=';TRAN.OUT';)
READ(1,1)ME,JMAX,NTIM,MEX,JPR,LN
READ(1,2)C,U,S,Q,EM
    1 FORMAT(6I5)
    2 FORMAT(2F5.2,3E10.3)
C
      IF(JPR.EQ.1)WRITE(6,3)
IF(JPR.EQ.2)WRITE(6,4)
    3 FORMAT('ROPAGATING SINE-WAVE';)
    4 FORMAT('ROPAGATING TEMP-FRONT';)
      JMAF=JMAX-2
JMAP=JMAX-1
AMP=JMAP
DX=1.0/AMP
IF(JPR.EQ.2)DX=4./AMP
DT=C*DX/U
EL=LN
ALPH=S*DX*DX/DT
IF(ALPH.LT.1.0E-10)ALPH=1.0E-10
RCEL=U*DX/ALPH
QQ=Q*C/3
IF(ALPH.LT.3)QQ=0.
MQ=0
IF(ABS(QQ).GT.0.0001)MQ=1
ATIM=NTIM
TIM=0.
TIMAX=DT*ATIM
PI=3.141592654
C
      IF(ME.EQ.1)WRITE(6,5)ME
IF(ME.EQ.2)WRITE(6,6)ME
IF(ME.EQ.3)WRITE(6,7)ME
IF(ME.EQ.4)WRITE(6,8)ME
    5 FORMAT(';ME=';,I2,';   FTCS DIFFERENCING';)
    6 FORMAT(';ME=';,I2,';   LAX-WENDROFF';)
    7 FORMAT(';ME=';,I2,';   EXPLICIT 4PT UPWIND';)
    8 FORMAT(';ME=';,I2,';   GENERAL CRANK-NICOLSON';)
      WRITE(6,9)JMAX,NTIM,C,U,DX,DT
    9 FORMAT(';JMAX=';,I3,'; NTIM=';,I3,';C=';,F5.2,';U=';,F5.2,
     1'; DX=';,F5.3,'; DT=';,F5.3)
WRITE(6,10)S,ALPH,RCEL,Q,QQ
WRITE(6,11)MEX,LN,EM
   10 FORMAT('; S=';,F5.2,'; ALPH=';,E10.3,';RCEL=';,F6.3,
     1'; Q=';,F5.2,';QQ=';,F6.3)
   11 FORMAT(';MEX=';,I5,';EL=';,I5,';EM=';,E10.3)
C
      IF(ME.GT.3)GOTO 12
SS=S
IF(ME.EQ.2)SS=S+0.5*C*C
AA=(0.5*C+SS)+3.*QQ
BB=1.-2.*SS-3.*QQ
CC=-0.5*C+SS+QQ
GOTO 13
   12 AA=EM-0.25*C-0.5*S-1.5*QQ
      BB=1.0-2.0*EM+S+1.5*QQ
CC=EM+0.25*C-0.5*S-0.5*QQ
AE=EM+0.25*C+1.5*QQ+0.5*S
BE=1.-2.*EM-1.5*QQ-S
CE=EM-0.25*C+0.5*QQ+0.5*S
   13 WRITE(6,14)AA,BB,CC,AE,BE,CE
   14 FORMAT(';AA=';,F8.5,';BB=';,F8.5,';CC=';,F8.5,';AE=';,F8.5,';BE=';,
     1F8.5,';CE=';,F8.5,/)
C
C     INITIALISE T AND EVaLUATE TEX
C
      CALL EXSOL(JPR,JMAX,X,T,TEX,NEX,DX,U,ALPH,TIMAX,EL)
C
      DO 16 J=1,JMAX
DO 15 K=1,5
   15 A(K,J)=0.
   16 CONTINUE
      WRITE(6,17)TIM
   17 FORMAT(';INITIAL SOLUTION,TIM=';,F5.3)
      WRITE(6,18)(X(J),J=1,JMAX)
WRITE(6,19)(T(J),J=1,JMAX)
   18 FORMAT(';X=';,12F6.3)
   19 FORMAT(';T=';,12F6.3)
C
C     MARCH SOLUTION IN TIME
C
      DO 26 N=1,NTIM
IF(ME.GT.3)GOTO 21
C
C     EXPLICIT SCHEMES
C
      D(1)=T(1)
D(2)=T(1)
D(3)=T(2)
DO 20 J=2,JMAP
IF(ME.EQ.3)D(4)=D(1)
D(1)=D(2)
D(2)=D(3)
D(3)=T(J+1)
T(J)=AA*D(1)+BB*D(2)+CC*D(3)
IF(ME.EQ.3)T(J)=T(J)-QQ*D(4)
   20 CONTINUE
      GOTO 26
C
C     TRIDIAGONAL SYSTEM FOR IMPLICIT SCHEMES
C
   21 IF(ME.EQ.1)DIM=1.
      DO 22 J=2,JMAP
JM=J-1
A(1,JM)=0.5*QQ
A(2,JM)=AA
A(3,JM)=BB
A(4,JM)=CC
R(JM)=AE*T(JM)+BE*T(J)+CE*T(J+1)
IF(ME.EQ.1.AND.J.GT.2)DIM=T(JM-1)
IF(ME.EQ.1)R(JM)=R(JM)-0.5*QQ*DIM
   22 CONTINUE
      R(1)=R(1)-A(2,1)*T(1)
IF(ME.EQ.0)GOTO 24
R(1)=R(1)-A(1,1)*T(1)
R(2)=R(2)-A(1,2)*T(1)
C
C     REDUCE A TO TRIDIAGONAL FORM
C
      DO 23 JM=3,JMAP
JMM=JM-1
DUM=A(J,JM)/A(2,JMM)
A(2,JM)=A(2,JM)-A(3,JM)*DUM
A(3,JM)=A(3,JM)-A(4,JM)*DUM
A(1,JM)=0.
R(JM)=R(JM)-R(JMM)*DUM
   23 CONTINUE
      A(1,1)=0.
      A(1,2)=0.
   24 A(2,1)=0.
      A(4,JMAF)=0.
C
      CALL BANFAC(A,JMAP,1)
CALL BANSOL(P,TD,A,JMAP,1)
C
      DO 25 J=2,JMAP
   25 T(J)=TD(J-1)
   26 CONTINUE
C
      WRITE(6,27)TIMAX
   27 FORMAT(';FINAL SOLUTION,TIM=';,F5.3)
      WRITE(6,18)(X(J),J=1,JMAX)
WRITE(6,19)(T(J),J=1,JMAX)
WRITE(6,28)(TEX(J),J=1,JMAX)
   28 FORMAT(';TEX=';,12F6.3)
      SUM=0.
DO 29 J=2,JMAP
   29 SUM=SUM+(T(J)-TEX(J))**2
      RMS=SQRT(SUM/(AMP-1.))
WRITE(6,30)RMS
   30 FORMAT(';RMS ERR=';,E12.5)
      STOP
END
c
      SUBROUTINE EXSOL(JPR,JMAX,X,T,TEX,NEX,DX,U,ALPH,TIMAX,EL)
C
C     SETS THE INITIAL T SOLUTION AND FINAL EXACT (TEX) SOLUTION
C
      DIMENSION T(65),TEX(65),X(65)
JMAP=JMAX-1
PI=3.141592654
IF(JPR.EQ.1)XST=0.
IF(JPR.EQ.2)XST=-2.0
DO 1 J=1,JMAX
AJ=J-1
X(J)=XST+AJ*DX
      T(J)=0.
    1 TEX(J)=0.
      IF(JPR.EQ.2)GOTO 3
C
C     EXACT SOLUTION FOR PROPAGATING SINE-WAVE
C
      JM=0.1001/DX+1.0
INC=U*TIMAX/DX+0.001
DO 2 J=1,JM
T(J)=SIN(10.*PI*X(J))
JP=J+INC
    2 TEX(JP)=T(J)
      RETURN
C
C     EXACT SOLUTION FOR PROPAGATING TEMP-FRONT
C
    3 T(1)=1.0
      TEX(1)=1.0
DO 5 J=2,JMAP
IF(X(J).LT.0.)T(J)=1.0
IF(ABS(X(J)).LT.1.0E-04)T(J)=0.5
DEM=0.
DO 4 K=1,NEX
AK=2*K-1
DUM=AK*PI/EL
SNE=SIN(DUM*(X(J)-U*TIMAX))
DIM=-ALPH*DUM*DUM*TIMAX
IF(DIM.LT.-20.)GOTO 5
    4 DEN=DEN+(SNE/AK)*EXP(DIM)
    5 TEX(J)=0.5-2.*DEN/PI
      RETURN
END
这个程序是我在书上看到的一个程序，用来求解一维热传导方程的。我看了好久不大懂！请教各位高手有没有人见过这个程序！大体的错误好像是出在
           CALL BANFAC(A,JMAP,1)
CALL BANSOL(P,TD,A,JMAP,1)
这两句上面，因为我认为是banfac ,bansol子函数都没有定义！请教怎么解决，或者指导一下该怎么解决，确实第一次做，很多不懂的！各位给点建议！

a2883524

补充说明一些：
JPR  ＝1propagating sine_wave
     =2propagating temperature front
c      courant number
A      elements in tridiagonal matrix
banfac   factorises A into upper triangular
bansol   solves A T=R
此程序出自computational techniques for fluid dynamics 1   page 307