急求fortran 修改，有两个错误怎么也找不出来

我用F文件的格式调试了一下，解决了几个BUG，不知道能不能行，你试试

! F File，Fixed Format，Mar/18/2012

!C 梯形断面明槽水面曲线计算源程序：WTSF.FOR

EXTERNAL FHK,FHO,FE

REAL I,N,M,J1,J2

DIMENSION H(201),V(201),S(201)

COMMON Q,I,N,M,B,csn,srm,alfa,DS,DR,V2,J1,J2,ES1,ES2

OPEN(2,FILE='RESULTS.TXT') !输出数据文件

WRITE(*,1001)

WRITE(2,1001)

1001 FORMAT(4(5X,A//))

WRITE(*,*)'键盘输入Q,I,N,M,B'

READ(*,*) Q,I,N,M,B

WRITE(*,1000) Q,I,N,M,B

WRITE(2,1000) Q,I,N,M,B

1000 FORMAT(5X,'Q=',F6.2,'(m**3/s)',4x,'I=',F8.5,4X,

1 'N=',F6.4,5X,'M=',F4.2,4X,'B=',F7.2,'(m)')

alfa=1.05

csn=(1-I*I)**0.5

srm=(1+M*M)**0.5

!C 计算临界水深

HK=ERFENFA(FHK,0.0,40.0,0.0005)

WRITE(*,1010) HK

WRITE(2,1010) HK

1010 FORMAT(5X,'临界水深 HK=',F9.6,'(m)')

!c 计算正常水深，平坡、逆坡时取HO为一大值

IF(I.LE.0) THEN

HO=100

WRITE(*,*)'底坡i<=0,无正常水深'

WRITE(2,*)'底坡i<=0，无正常水深'

ELSE

HO=ERFENFA(FHO,0.0,40.0,0.0005)

WRITE(*,1020) HO

WRITE(2,1020) HO

ENDIF

1020 FORMAT(5X,'HO=',F7.4,'(m)')

WRITE(*,1021)

1021 FORMAT(//3(5X,A/)/)

WRITE(*,*)'input HD,DS,NS'

READ(*,*) HD,DS,NS

H(1)=HD

A=(B+M*H(1))*H(1)

R=A/(B+2*H(1)*srm)

V(1)=Q/A

JI=(V(1)*N)**2/R**(4.0/3)

ES1=csn*H(1)+alfa*V(1)**2/19.6

!c 判断水面线计算方向：DR=1，控制断面在上游；DR=-1,控制断面在下游

IF ((HD>HK).OR.((HD==HK).AND.(HD<HO))) THEN

DR=-1

ELSE

DR=1

ENDIF

!C 二分法区间端点HB取值

IF (((HK>HO).AND.(HD>HK)).OR.

1 ((HO>HK).AND.(HD<HK))) THEN

HB=HK

ELSE

HB=HO

ENDIF

S(1)=0.0

!C 计算各断面水深H(L)和流速V(L)

DO 10 L=2,NS+1

IF(ABS(H(L-1)-HO)<0.0005) THEN

H(L)=HO

V(L)=V(L-1)

ELSE

H(L)=ERFENFA(FE,H(L-1),HB,0.00001)

V(L)=V2

ES1=ES2

ENDIF

S(L)=(L-1)*DS

10 CONTINUE

!C 输出计算结果

WRITE(2,*)

WRITE(2,1030) HD,DS,NS

1030 FORMAT(5X,'控制水深HD=',F6.3,'(m)',5X,'断面间距DS=',F7.2,'(m)',

1 5X,'间距总数NS=',I3/)

IF (DR<0) THEN

WRITE(2,*)'缓流,控制断面在下游'

ELSE

WRITE(2,*)'急流,控制断面在上游'

ENDIF

WRITE(2,1035)

1035 FORMAT(/7X,'L',10X,'H(L)',10X,'V(L)',9X,'S(L)'/18X,

1 '(m)',10X,'(m/s)',9X,'(m)'/5X,

2 '-------------------------------------')

WRITE(2,1040) (L,H(L),V(L),S(L),L=1,NS+1)

1040 FORMAT(5X,I3,7X,F7.3,7X,F7.3,4X,F10.2)

END

!c 求临界水深函数

FUNCTION FHK(H)

REAL I,N,M

COMMON Q,I,N,M,B,csn,srm,alfa,DS,DR,V2,J1,J2,ES1,ES2

FHK=9.8*csn*((B+M*H)*H)**3-alfa*Q*Q*(B+2*M*H)

END

FUNCTION FHO(H)

REAL I,N,M

COMMON Q,I,N,M,B,csn,srm,alfa,DS,DR,V2,J1,J2,ES1,ES2

FHO=Q*N/I**0.5*(B+2*srm*H)**(2.0/3)-((B+M*H)*H**(5.0/3))

END

FUNCTION FE(H)

REAL I,N,M,J1,J2

COMMON Q,I,N,M,B,csn,srm,alfa,DS,DR,V2,J1,J2,ES1,ES2

A=(B+M*H)*H

V2=Q/A

J2=(N*V2)**2/(A/(B+2*H*srm))**(4.0/3)

ES2=csn*H+alfa*V2*V2/19.6

FE=ES1-ES2+DR*(I-(J1+J2)/2)*DS

END

! 二分法函数程序

FUNCTION ERFENFA (F,X1,X2,EPS)

A=X1

B=X2

10 FA=F(A)

FB=F(B)

IF(FA*FB.GT.0)THEN

WRITE(*,*)'(X1，X2)不是有根区间,请重新输入X1,X2'

READ(*,*) A,B

GOTO 10

ENDIF

DO 50 I=1,30

ERFENFA=(A+B)*0.5

IF (ABS(B-A).LT.EPS) RETURN

FM=F(ERFENFA)

IF (FM*FA.LT.0) THEN

B=ERFENFA

ELSE

A=ERFENFA

ENDIF

50 CONTINUE

END

---