不知道有没人用过这个程序算过边坡的安全系数，我输入数值后，就是出不了结果，你哪位高手指点一下，输入数值时有什么需要注意的（以下红色部分为我输入的值），还有圆心的搜索范围是自己取的，大概多大的范围比较合适，拜托各位了，谢谢！
C
C PROGRAM OF SLOP STABILITY(BISHOP METHOD)
C
C FUNCTION---
C AUTOMATIC SEARCH CRITICAL CIRCULAR SURFACE USE
C OPTIMUM METHOD AND CALCULATED SLIDING FACTOR
C OF SAFTY
DIMENSION XX(5),YY(5),RA(5),FA(5)
DIMENSION XU(24),YU(10,24),WL(24),QA(23),
* C(10),T(10),GS(10),X(41),YW(41),Y(10,41),Q(41)
REAL XA,XB,YA,YB,EQH,TOL1,TOL2,DS,DX
c COMMON /NB/NB,NS,N/C/C,T,GS
c COMMON /TOL/TOL1,TOL2/EQ/EQH/QA/QA
c COMMON /XU/XU/YU/YU/WL/WL/DS/DS
c COMMON /XA/XA,XB,YA,YB
C INPUT DATA
C
DATA XU/60.,135.,147.,159.,201.74,245.,267.,
* 17*0.0/
DATA YU/2*162.0,164.0,177.0,190.,5*0.0,156.497,162.0,164.0,
* 177.0,190.,5*0.0,152.807,162.0,164.0,177.0,190.,5*0.0,
* 150.99,162.03,164.01,177.0,190.,5*0.0,145.941,162.0,
* 164.0,177.0,190.,5*0.0,146.0,162.0,164.0,177.0,190.,
* 5*0.0,
* 146.03,162.0,164.20,177.0,190.,5*0.0,170*0.0/
DATA QA/23*0.0/WL/7*162.0,17*0.0/
DATA C/1.960,78.3,7.843,2.941,6*0.0/T/35.0,4.93,25.0,30.0,6*0.0/
DATA GS/20.5,11.515,18.,19.5,6*0.0/
DATA NB,NS,N/8,4,20/EQH/0./
DATA XA,XB,YA,YB/0.,200.,0.,200.0/DS/0.0/
DATA TOL1,TOL2/0.02,10./
OPEN(2,FILE='OUTPUT.DAT',STATUS='NEW')
C OPEN(1,FILE='FDATA1',STATUS='OLD')
C READ(1,*)XU
C READ(1,*)YU
C READ(1,*)WL
C READ(1,*)QA
C READ(1,*)C
C READ(1,*)T
C READ(1,*)GS
C READ(1,*)NB,NS,N
C READ(1,*)XA,XB,YA,YB
C READ(1,*)DS
C READ(1,*)EQH,TOL1,TOL2
C WRITE(*,*)'IF RETURN HOME OR NOT'
C WRITE(2,*)'*****XU(1)******',XU(1)
WRITE(*,*)'Begining Calculation'
DX=(XU(NB)-XU(1))/FLOAT(N)
X(1)=XU(1)
Y(1,1)=YU(1,1)
IF(NS.EQ.1) GOTO 26
DO 20 K=2,NS
Y(K,1)=YU(K,1)
IF(Y(K,1).LT.Y(K-1,1)) Y(K,1)=Y(K-1,1)
20 CONTINUE
26 YW(1)=WL(1)
J=2
N1=N+1
DO 21 I=2,N1
X(I)=XU(1)+DX*FLOAT(I-1)
23 IF(X(I).LE.XU(J)) GO TO 22
J=J+1
GOTO 23
22 A=(X(I)-XU(J-1))/(XU(J)-XU(J-1))
Y(1,I)=YU(1,J-1)+A*(YU(1,J)-YU(1,J-1))
C WRITE(*,701) Y(1,I)
C701 FORMAT(10X,'Y=',F8.2)
IF(NS.EQ.1) GOTO 25
DO 24 K=2,NS
Y(K,I)=YU(K,J-1)+A*(YU(K,J)-YU(K,J-1))
24 CONTINUE
25 YW(I)=WL(J-1)+(WL(J)-WL(J-1))*A
IF(YW(I).LT.Y(1,I)) YW(I)=Y(1,I)
Q(I-1)=QA(J-1)
21 CONTINUE
C WRITE(*,*)'IF RETURN OR NOT'
WRITE(2,501)
501 FORMAT(//20X,'COODINATERS OF CIRCLE CENTRE')
WRITE(2,208)XA,XB,YA,YB
208 FORMAT(18X,F7.2,1X,10H*** XO ***F7.2/18X,F7.2,1X,
* 10H*** YO ***,F7.2)
WRITE(2,502)
502 FORMAT(/,10X,60(1H+))
C CALL CENTER
WRITE(*,*)'Running Over!'
C STOP
C END
C SUBROUTINE CENTER
C
C FUNCTION--
C SEARCH CRITICAL CENTER OF SLIDING RADIUS
C
C DIMENSION XX(5),YY(5),RA(5),FA(5)
C COMMON /NB/NB,NS,N/C/C(10),T(10),GS(10)//X(41),
C * YW(41),Y(10,41),Q(41),DX
C COMMON /XA/XA,XB,YA,YB/TOL/TOL1,TOL2/DS/DS
WRITE(2,90)
90 FORMAT(15X,2HXC,14X,2HYC,14X,2HRC,14X,2HFA)
DXA=(XB-XA)/4.0
DYA=(YB-YA)/4.0
XX(1)=(XA+XB)/2.0
YY(1)=(YA+YB)/2.0
C WRITE(2,300)
C300 FORMAT(1X,'ERROR(1)')
c WRITE(2,*)'%%%%%%%%%%XX(1),YY(1)%%%%%',XX(1),YY(1)
CALL RADIUS (XX(1),YY(1),RA(1),FA(1),XU,YU,NB,NS,
* N,C,T,GS,X,YW,Y,Q,DX,XA,XB,YA,YB,TOL1,TOL2,DS)
WRITE(2,100) XX(1),YY(1),RA(1),FA(1)
100 FORMAT(1X,5HFA(1),1X,E15.7,1X,E15.7,1X,E15.7,1X,E15.7)
10 XX(2)=XX(1)-DXA
XX(3)=XX(2)
XX(4)=XX(1)+DXA
XX(5)=XX(4)
YY(2)=YY(1)-DYA
YY(3)=YY(1)+DYA
YY(4)=YY(2)
YY(5)=YY(3)
J=1
FM=FA(1)
FAMIN=FA(1)
DO 14 I=2,5
CALL RADIUS(XX(I),YY(I),RA(I),FA(I),XU,YU,NB,NS,
* N,C,T,GS,X,YW,Y,Q,DX,XA,XB,YA,YB,TOL1,TOL2,DS)
WRITE(2,110) I,XX(I),YY(I),RA(I),FA(I)
110 FORMAT(1X,3HFA(,I1,1H),1X,E15.7,1X,E15.7,1X,E15.7,1X,E15.7)
FM=FM+FA(I)
IF(FAMIN.LT.FA(I)) GOTO 14
J=I
FAMIN=FA(J)
14 CONTINUE
FM=FM/5.0
AA=ABS((FAMIN-FM)/FAMIN)
IF(AA.LE.TOL1) GOTO 16
FA(1)=FA(J)
XX(1)=XX(J)
YY(1)=YY(J)
DXA=DXA/2.0
DYA=DYA/2.0
GOTO 10
C WRITE(*,*)'&&&&&&&&&&&&&&&&& *THE FIRST&************'
16 WRITE(2,55)
55 FORMAT(/,10X,60(1H+))
WRITE(2,200) XX(J),YY(J),RA(J),FA(J)
200 FORMAT(/,10X,'CENTERS COODINATE OF SLIDING',
* 'RADIUS XC=',F7.3/,46X,'YC=',F7.3/47X,'R=',F7.3/47X,
* 'FACTOR OF SAFTY F=',F7.4/,10X,60(1H+))
C RETURN
END
SUBROUTINE RADIUS(XC,YC,RC,FC,XU,YU,NB,NS,
* N,C,T,GS,X,YW,Y,Q,DX,XA,XB,YA,YB,TOL1,TOL2,DS)
C
C FUCTION--
C SURCH CRITICAL RADIUS
C
REAL RB(3),FB(3),XU(24),YU(10,24),C(10),
* T(10),GS(10),X(41),YW(41),Y(41),Q(41),DX,XA,XB,YA,YB,TOL1,TOL2,DS
C COMMON /XU/XU(24)/YU/YU(10,24)/NB/NB,NS,N
C COMMON /C/C(10),T(10),GS(10)//X(41),YW(41),Y(10,41),Q(41),
C * DX/XA/XA,XB,YA,YB
C COMMON /TOL/TOL1,TOL2/DS/DS
C WRITE(2,*)'ENTER THE FIRST RADIUS'
C WRITE(2,*)'THE XU1=',XU(1)
R1=SQRT((XC-XU(1))**2+(YC-YU(1,1))**2)
C WRITE(2,*)'THE R1=',R1
R2=SQRT((XC-XU(NB))**2+(YC-YU(1,NB))**2)
C WRITE(2,901)R1,R2
C901 FORMAT(10X,2F8.2)
RMAX=AMIN1(R1,R2)
RMIN=RMAX
DO 10 I=2,NB
R1=SQRT((XC-XU(I))**2+(YC-YU(1,I))**2)
IF(R1.LT.RMIN) RMIN=R1
10 CONTINUE
WRITE(2,62) RMIN
62 FORMAT(/,4X,'***TO SUB-R 18 RMIN=',F8.3)
DR=(RMAX-RMIN)/4.0
RB(1)=(RMAX+RMIN)/2.0
C WRITE(2,*)'THE ABOVE FIRST FACTOR'
CALL FACTOR(RB(1),XC,YC,FB(1),X,Y,YW,N,NS,Q,C,T,GS,DX)
WRITE(2,*)'***************THE FIRST FACTOR*********'
C WRITE(2,40)XC,YC,RB(1),FB(1)
C40 FORMAT(1X,5HFB(1),1X,E15.7,1X,E15.7,1X,E15.7,1X,E15.7)
12 RB(2)=RB(1)-DR
RB(3)=RB(1)+DR
CALL FACTOR(RB(2),XC,YC,FB(2),X,Y,YW,N,NS,Q,C,T,GS,DX)
C WRITE(2,42) XC,YC,RB(2),FB(2)
C42 FORMAT(1X,5HFB(2),1X,E15.7,1X,E15.7,1X,E15.7,1X,E15.7)
CALL FACTOR(RB(3),XC,YC,FB(3),X,Y,YW,N,NS,Q,C,T,GS,DX)
C WRITE(2,44) XC,YC,RB(3),FB(3)
C44 FORMAT(1X,'FB(3)',1X,E15.7,1X,E15.7,1X,E15.7,1X,E15.7)
FM=(FB(1)+FB(2)+FB(3))/3.0
J=1
FMIN=FB(1)
DO 14 I=2,3
IF(FMIN.LT.FB(I)) GOTO 14
J=I
FMIN=FB(J)
14 CONTINUE
IF(ABS((FMIN-FM)/FMIN).LT.TOL1) GOTO 16
RB(1)=RB(J)
FB(1)=FMIN
DR=DR/2.0
GOTO 12
16 FC=FMIN
RC=RB(J)
C WRITE(*,50)
C50 FORMAT(5X,'*END SUB-RADIUS*')
RETURN
END
SUBROUTINE FACTOR(R,XC,YC,F,X,Y,YW,N,NS,Q,C,T,GS,DX)
C
C FUCTION--
C CIRCULAR SLIDING FACTOR OF SAFTY
C
REAl X(41),Y(10,41),YW(41),Q(41),C(10),T(10),GS(10)
DIMENSION DL(60),CA(60),SA(60),YR(61)
C COMMON /NB/NB,NS,N/C/C(10),T(10),GS(10)
C COMMON /EQ/EQH//X(41),YW(41),Y(10,41),Q(41),DX
C COMMON /XA/XA,XB,YA,YB/TOL/TOL1,TOL2/DS/DS
WRITE(2,901)R,XC,YC
901 FORMAT(10X,'R,XC,YC=',3F8.2)
N1=N+1
DO 10 I=1,N1
IF(ABS(XC-X(I)).GT.R) GOTO 11
YR(I)=SQRT(R**2-(X(I)-XC)**2)+YC
IF(YR(I).GE.Y(1,I)) GOTO 10
11 YR(I)=Y(1,I)
10 CONTINUE
C WRITE(*,*)'**** OUT OF QUESION********'
DO 15 I=1,N
DL(I)=SQRT((X(I+1)-X(I))**2+(YR(I)-YR(I+1))**2)
CA(I)=(X(I+1)-X(I))/DL(I)
SA(I)=(YR(I)-YR(I+1))/DL(I)
WRITE(2,*)I,DL(I),SA(I),CA(I)
500 FORMAT(10X,'I=',I2,'DL(I)=',F6.3,'SA(I)=',F7.4,'CA(I)=',F7.4)
15 CONTINUE
WRITE(2,*)'N TIME RECLE IS OVER'
F=1.0
16 WS=0.0
Z=0.0
DO 20 I=1,N
YRM=(YR(I)+YR(I+1))/2.0
YM=(Y(1,I)+Y(1,I+1))/2.0
YT=YM
YEQ=(YT+YRM)/2.0
YWM=(YW(I)+YW(I+1))/2.0
IF(YT.GE.YWM) U=YRM-YT
U=YRM-YWM
IF(U.LT.0.0)U=0.0
C IF(YRM.LE.YM)GOTO 20
W=DX*Q(I)
IF(NS.EQ.1) GOTO 25
NS1=NS-1
DO 22 J=1,NS1
YM=(Y(J+1,I+1)+Y(J+1,I))/2.0
IF(YRM.GT.YM) GO TO 22
JA=J
GOTO 24
22 W=W+GS(J)*DX*(YM-(Y(J,I+1)+Y(J,I))/2.0)
25 JA=NS
24 W=W+GS(JA)*DX*(YRM-(Y(JA,I+1)+Y(JA,I))/2.0)
WQH=0.25*W*EQH
WQV=WQH/3.0
WS=WS+(W+WQV)*SA(I)+WQH*YEQ/YRM
SIT=SIN(3.141592*T(JA)/180.0)
COT=COS(3.141592*T(JA)/180.0)
ZU=C(JA)*DL(I)*CA(I)+(W+WQV-U*DL(I)*CA(I))*SIT/COT
ZL=CA(I)+SA(I)*SIT/(COT*F)
Z=Z+ZU/ZL
20 CONTINUE
WRITE(2,*)'TO THIS POSITION'
WRITE(2,*)' WS=',WS
FA=Z/WS
250 FORMAT(1X,'FC=',E15.7)
IF(ABS(FA-F)/FA.LE.TOL1) GOTO 30
F=FA
GOTO 16
30 F=FA
IF(F.LE.0.0) F=100000.0
RETURN
END
