e4.doc

边坡稳定上限解程序EMU

	OPEN(5,FILE='     ',STATUS='UNKNOWN')
	OPEN(6,FILE='     ',STATUS='UNKNOWN')
	F1=1.1
	CALL READ_DATA
	CALL  FINDFOS(FOS,F1)
	STOP
	END
	SUBROUTINE READ_DATA
INCLUDE 'IN4.FOR'
	O.PI=3.1415927
	 !水容重，坡外水位，水平地震力系数，拉力缝控制符，安全系数控制符
	 READ(5,*)UWOW,WLOS,COEQ,TENSION,OFOS
	WRITE(6,*)'水容重,坡外水位,水平地震力系数,拉力缝,安全系数控制符',
	& UWOW,WLOS,COEQ,TENSION,OFOS
	 READ(5,*)S.N  !土条界面总数
	 WRITE(6,*)'土条界面总数',S.N
	 WRITE(6,*)'界面有关数据'
	! NU=界面与编号为NU的外边坡线相交，XU,YU,X,Y, 界面顶部和底部坐标，L=界面长度
WRITE(6,*)'I,NU,XU,YU,X,Y,L'
DO I=1,S.N
	 READ (5,*)I,S.NU(I),S.XU(I),S.YU(I),S.X2(I),S.Y2(I),S.L(I)
WRITE(6,27)I,S.NU(I),S.XU(I),S.YU(I),S.X2(I),S.Y2(I),S.L(I)
ENDDO
	 !G=界面上的平均容重,FW,CW=界面强度指标,P=孔压,DELTA=界面倾角
WRITE(6,*)'B:I,G,FW,CW,P,DELTA'
DO I=1,S.N
READ(5,*)I,B.G(I),B.FW(I),B.CW(I),B.P(I),S.D(I)
	 B.FW(I)=DTR(B.FW(I))
	 S.D(I)=DTR(S.D(I))
WRITE(6,25)I,B.G(I),RTD(B.FW(I)),B.CW(I),B.P(I),RTD(S.D(I))
ENDDO
	 ! 土条底面中点数据，X,Y=坐标，ALPHA=倾角,C,F=强度指标
WRITE(6,*)'DATA FOR THE MID-POINT OF SLICES'
WRITE(6,*)'I,Q:X,Y,ALPHA,C,F'
KG=S.N-1
DO I=1,KG
READ(5,*)I,Q.X(I),Q.Y(I),Q.A(I),Q.C(I),Q.F(I)
	 Q.A(I)=DTR(Q.A(I))
	 Q.F(I)=DTR(Q.F(I))
WRITE(6,25)I,Q.X(I),Q.Y(I),RTD(Q.A(I)),Q.C(I),RTD(Q.F(I))
ENDDO
	! 条块数据，W=重量，P=条底孔压，QX，QY=水平，垂直外荷，
WRITE(6,*)'I,Q: WEIGHT,PORE PRESSURE,QX,QY'
DO I=1,KG
READ(5,*)I,Q.W(I),Q.P(I),SL.QX(I),SL.QY(I)
WRITE(6,25)I,Q.W(I),Q.P(I),SL.QX(I),SL.QY(I)
ENDDO
	!PHX，PHY=条块界面浸润线坐标
WRITE(6,*)'I,ID,PHX,PHY'
DO I=1,S.N
	 READ(5,*)I,S.PHX(I),S.PHY(I)
WRITE(6,25)I,S.PHX(I),S.PHY(I)
ENDDO
CONTINUE
25  FORMAT(1X,I5,7F10.3)
27  FORMAT(1X,2I5,6F10.3)
	RETURN
	END
SUBROUTINE FINDFOS(FOS,F1)
	AER=0.0001
	WRITE(6,*)'ITERATION PROCESS'
	WRITE(6,*)'KM,ETA,FOS'
WRITE(6,*)KM,ETA,F1
CALL CALWE(F1,EWT,WSM,WSL,E1,1)
KM=0
F2=F1*1.1
CALL CALWE(F2,EWT,WSM,WSL,E2,1)
AA=(E2-E1)/E1
F0=F1-(F2-F1)/AA
F1=F0
20	CONTINUE
CALL CALWE(F1,EWT,WSM,WSL,ETA,0)
WRITE(6,*)KM,ETA,F1
CALL DERIV(F1,DRF)
DF=-ETA/DRF
IF(ABS(DF).LT.AER)GOTO 10
IF(KM.GT.10)THEN
WRITE(6,*)'ITERATION TROUBLE IN FINDING FOS'
ID=1
RETURN
ENDIF
F1=F1+DF
KM=KM+1
GOTO 20
10  CONTINUE
FOS=F1
	 WRITE(6,*)'FOS=',FOS
RETURN
END
C
SUBROUTINE DERIV(F,DRF)
INCLUDE 'IN4.FOR'
DF=0.05*F
F1=F+DF
CALL CALWE(F1,EWT,WSM,WSL,ETA1,1)
F2=F-DF
CALL CALWE(F2,EWT,WSM,WSL,ETA2,1)
DRF=(ETA2-ETA1)/(F2-F1)
RETURN
END
c
SUBROUTINE CALWE(FOS,EWT,WSM,WSL,ETA,IFD)
DIMENSION FE(60),VS(60),FWE(60),SKS(60),GTEM(60),QQX(60),QQY(60)
#,FWF(60),ICS(60),CE(60),CWE(60),BAX(60),BAY(60),BX(60),BY(60)
COMMON/SUMXY/SIGX,SIGY,LWW,SIGXP,SIGYP,ICS
COMMON/VEL/FE,VS
COMMON/ITEN/ITENSION,KTEN,HTENSION
COMMON/OP6/IP1,IP2,IP3,IP4,IP5,IP6
COMMON/FILES/PFILE,FILE5,FILE6,FILE7,FILE8
CHARACTER*70 PFILE,FILE5,FILE6,FILE7,FILE8
DIMENSION WW(60)
	INCLUDE 'IN4.FOR'
real*8 ax,dea, dei,dei1,vs
IF(LWW.NE.0)IFD=1
IF(IFD.EQ.0)OPEN(10,FILE='DEB.EMU',STATUS='UNKNOWN')
FC.G(1)=0.
WG=0.
WSL=0.
WSF=0.
WSM=0.
FM=0.
EANC=0.
ECB=0.
NT=S.N-1
VS(1)=1.0d0
FE(1)=ATAN(TAN(Q.F(1))/FOS)
ce(1)=q.c(1)/fos
DO I=1,NT
CE(I+1)=Q.C(I+1)/FOS
FE(I+1)=ATAN(TAN(Q.F(I+1))/FOS)
IF(I.EQ.NT)FE(I+1)=FE(I)
AX=Q.A(I+1)-Q.A(I)-FE(I+1)+FE(I)
IF(AX.GE.-0.01.OR.IP4.EQ.1)THEN  !1.7.1993
ICS(I+1)=0
IF(I.EQ.1)FWE(1)=ATAN(TAN(B.FW(1))/FOS)
FWE(I+1)=ATAN(TAN(B.FW(I+1))/FOS)
c      ii=i+1
c      write(6,*)'i,fos ',i,fos,b.fw(i+1),fwe(i+1)
FWF(I+1)=FWE(I+1)
IF(I.EQ.1)FWF(1)=FWE(1)
ELSE
ICS(I+1)=1
IF(I.EQ.1)FWE(1)=O.PI-ATAN(TAN(B.FW(1))/FOS)
FWE(I+1)=O.PI-ATAN(TAN(B.FW(I+1))/FOS)
FWF(I+1)=-ATAN(TAN(B.FW(I+1))/FOS)
IF(I.EQ.1)FWF(1)=-ATAN(TAN(B.FW(1))/FOS)
ENDIF
IF(I.EQ.NT)ICS(I+1)=ICS(I)
SKS(I)=Q.A(I)+S.D(I)-FWF(I+1)-FE(I)
IF(I.EQ.NT)THEN
II=I+1
SKS(II)=Q.A(II)+S.D(II)-FWF(II+1)-FE(II)
ENDIF
IF(ABS(Q.F(I+1)-Q.F(I)).LT.0.01.AND.
% ABS(Q.A(I+1)-Q.A(I)).LT.0.05)THEN
DEA=Q.A(I+1)-Q.A(I)
AVA=(Q.A(I+1)+Q.A(I))/2.
DEI=AVA+O.PI/2.+S.D(I+1)-FWE(I+1)-FE(I+1)
DEI1=Q.A(I)+O.PI/2.+S.D(I+1)-FWE(I+1)-FE(I)
IF(AX.GE.-0.01.AND.DEI1.LT.0.05)ID=1  !1.7.1993
IF(AX.LT.-0.01.AND.DEI1.GT.-0.05)ID=2 !1.7.1993
IF(IP4.EQ.1)ID=0
IF(ID.NE.0)GOTO 110
VS(I+1)=VS(I)*dexp(-DEA/dTAN(DEI))
IF(VS(I+1).LT.-0.01.AND.IP4.NE.1)ID=3
ELSE
EXPL=Q.A(I)+O.PI/2.+S.D(I+1)-FWE(I+1)-FE(I)
EXPR=Q.A(I+1)+O.PI/2.+S.D(I+1)-FWE(I+1)-FE(I+1)
VS(I+1)=VS(I)*SIN(EXPL)/SIN(EXPR)
ENDIF
IF(VS(I+1).LT.-0.01.AND.IP4.NE.1)ID=4
IF(ID.EQ.4)THEN
WRITE(6,*)'EXP=Q.A+PI/2+D-FWE-FE'
WRITE(6,*)'A,A+1,D,FW,F',Q.A(I),Q.A(I+1),S.D(I+1),FWE(I+1)
$  ,FE(I+1)
WRITE(6,*)'EXPL,EXPR=', EXPL,EXPR
ENDIF
10  DX=S.X2(I+1)-S.X2(I)
IF(VS(I+1).LT.-0.01.AND.IP4.NE.1)GOTO 110
GOTO 113
110  CONTINUE
IF(I.EQ.1.OR.I.EQ.S.N)GOTO 113
WRITE(6,*)'KINEMATIC CONDITION VIOLATED AT NO. ',I
WRITE(6,*)
%  'DEI1=Q.A(I)+O.PI/2.+S.D(I+1)-FWE(I+1)-FE(I),AX=DALFA-DPHI'
WRITE(6,*)'AX,DEI1=',AX,DEI1
WRITE(6,*)'Q.A(I),Q.A(I+1),S.D(I+1),FWE(I+1),FE(I+1)=, IN RADIUS'
$   ,Q.A(I),Q.A(I+1),S.D(I+1),FWE(I+1),FE(I+1)
	   ID=1
	 RETURN
113  CONTINUE
	BAX(I)=0.
	BAY(I)=0.
	BX(I)=0.
	BY(I)=0.
	WW(I)=UWOW*(S.PHY(I)-S.Y2(I)+S.PHY(I+1)-S.Y2(I+1))/2.
	IF(ANC.N.EQ.0)GOTO 815
	DO IH=1,ANC.N
	IF(ABS(ANC.NID(IH)).EQ.1)CYCLE
237    IF(ANC.NL(IH).EQ.I)THEN
	TAA=COS(Q.A(I)-FE(I)+ANC.A(IH))
	IF(TAA.GT.0.0)EANC=EANC+VS(I)*TAA*ANC.TA(IH)
	ENDIF
	IF(ANC.NF(IH).EQ.I)THEN
	BAY(I)=BAY(I)+ANC.TA(IH)*SIN(ANC.A(IH))*SIN(Q.A(I)-FE(I))
	BAX(I)=BAX(I)-ANC.TA(IH)*COS(ANC.A(IH))*COS(Q.A(I)-FE(I))
	BY(I)=BY(I)+ANC.TA(IH)*SIN(ANC.A(IH))
	BX(I)=BX(I)-ANC.TA(IH)*COS(ANC.A(IH))
	ENDIF
	ENDDO
815  CONTINUE
	QQX(I)=SL.QX(I)
	QQY(I)=SL.QY(I)
DECB=(CE(I)*COS(FE(I))-Q.P(I)*SIN(FE(I)))
%*DX/COS(Q.A(I))
ECB=ECB+DECB*VS(I)
DWG=(Q.W(I)+QQY(I))*SIN(Q.A(I)-FE(I))
DWG1=DWG+BAY(I)
WG=WG+DWG*VS(I)
C     WSL- SURFACE LOAD      ; WSM- BODY FORCE
IF(OFOS.EQ.1.OR.OFOS.EQ.-1)WSL=WSL+QQY(I)*SIN(Q.A(I)
% -FE(I))*VS(I)
C      WSF=WSF+(COEQ*Q.W(I)+QQX(I))*COS(Q.A(I)-FE(I))*VS(I)
DWSF=(COEQ*(Q.W(I)+WW(I))+QQX(I))*COS(Q.A(I)-FE(I))
DWSF1=DWSF+BAX(I)
WSF=WSF+DWSF*VS(I)
IF(OFOS.EQ.1.OR.OFOS.EQ.-1)
% WSL=WSL+QQX(I)*COS(Q.A(I)-FE(I))*VS(I)
IF(OFOS.EQ.2.OR.OFOS.EQ.-2.OR.OFOS.EQ.0)
% WSM=WSM+Q.W(I)*COS(Q.A(I)-FE(I))*VS(I)
GTEM(I)=DECB-DWG1-DWSF1
ENDDO
FC.G(S.N)=0.
B.PW(S.N)=0.
PTEN=0.
FTEN=0.
IF(ITENSION.NE.0)THEN
HHH=S.L((S.N))
IF(KTEN.EQ.1)HHH=HTENSION
PTEN=0.5*HHH**2*UWOW*COS(Q.A(NT)-FE(NT))*VS(NT)
FTEN=0.5*HHH**2*UWOW
ENDIF
EJ=0.
FC.G(1)=0.
FC.GN(1)=0.
FC.GT(1)=0.
DC1=0.
DK1=0.
IF(IFD.EQ.0)SIGX=PTEN
IF(IFD.EQ.0)SIGY=0.
TW=0
DO I=2,S.N
DX=S.X2(I)-S.X2(I-1)
I1=I-1
TW=TW+Q.W(I1)+QQY(I1)+BY(I1)
c      write(6,*)'i1,tw,qqy,by',i1,tw,qqy(i1),by(i1)
CWE(I)=B.CW(I)/FOS
IF(I.EQ.2)CWE(1)=CWE(2)
DA=Q.A(I)-Q.A(I-1)-FE(I)+FE(I-1)
IF(DA.LT.0.AND.DA.GT.-0.05)DA=0.
AG=Q.A(I)+O.PI/2.-FWE(I)+S.D(I)-FE(I)
AX=SIN(DA)/SIN(AG)
IF(I.EQ.S.N)AX=1.0
IF(AX.GE.-0.01)THEN   !1.7.1993
EJ=EJ+S.L(I)*(CWE(I)*ABS(COS(FWE(I)))-B.PW(I)*SIN(FWE(I)))*AX
$*VS(I-1)
IF(ANC.N.NE.0)THEN
DO IK=1,ANC.N
	IF(ANC.NID(IK).LT.0)CYCLE
	IF(I1.GE.ANC.NF(IK).AND.I1.LE.(ANC.NL(IK)-1))THEN
	TQ=FWE(I)-O.PI/2.
	IF(ICS(I).EQ.1)TQ=-TQ
	TTA=COS(TQ-S.D(I)+ANC.A(IK))
	IF(TTA.GT.0)EANC=EANC+ANC.TA(IK)*TTA*AX*VS(I-1)
	ENDIF
ENDDO
ENDIF
ELSE
ID=1
ENDIF
ENDDO
EWT=ECB-WG-WSF+EJ-PTEN+EANC
IF(OFOS.EQ.1.OR.OFOS.EQ.-1)ETA=EWT/WSL
IF(OFOS.EQ.2.OR.OFOS.EQ.0.OR.OFOS.EQ.-2)ETA=EWT/WSM
193   format(1x,5f8.2)
RETURN
END
FUNCTION RTD(A)
RTD=A*180./3.14159
END
FUNCTION DTR(A)
DTR=A*3.14159/180.
END

	STRUCTURE /SLIPW/   !S
	INTEGER*4 N,NU(60)
	REAL*4 X2(60),Y2(60),XU(60),YU(60),D(60),PHX(60),PHY(60)
	REAL*4 L(60)
	END STRUCTURE
	STRUCTURE/SLIPB/!B
	REAL*4 G(60),PWA(60),PW(60),CW(60),FW(60),CCW(60),P(60)
	END STRUCTURE
C	N=TOTAOL NUMBER OF POINTS IDENTIFYING THE BORDER
C	G=AVERAGE UNIT WEIGHT, P=PORE PRESSURE AT THE BASE
C	PW=TOTAL PORE PRESSURE ON THE WALL, PH=Y VALUE OF PHREATIC LINE
C	PWA=PW, IN TERMS OF PORE PRESSURE ANGLE
C	C(W),F(W),P(W)= COHESION,FRICTION, PORE PRESURE AT BASE OR WALL
	STRUCTURE /NODE/!N
	INTEGER*4 N,NN(20),NALF,NST,DF(20),NNALF(20),NNST(20),NU(20)
	INTEGER*4 ID,MOP,ITR,NAB(20),INTF(20)
	REAL*4 X(20),Y(20),D(20),M(20),W(20),XU(20),YU(20),DX(20)
C	X,Y= COORDINATES. DF=0, UNMOVED POINTS; DF=1, MOVE ALONG THE SPECIFIED
C	DIRECTION; DF=-?, MOVE ALONG THE LINE NUMBERED ?.
C	NN=SLICE NUMBER AT A NODE. NABR=NUMBER OF DISCONTINUITY IN ALPHA
C	NST=NOMBER OF STRAIGHT LINE SEGMENTS, NU=LINE NO. OF INTERCEPTION PT.
C	D=DELTA MEASURED WITH RESPECT TO POSITIVE X AXIS
C	M=DIRECTION OF MOVEMENT WITH RESPECT TO POSITIVE X AXIS
C	W=BANDWITH IN RANDOM SEARCH, XU,YU=COORDINATES OF INTERCEPTION PT.
	END STRUCTURE
	STRUCTURE /OPTION/    !O
	INTEGER*4 P1,P2,P3,P4,P5,P6,W1,W2,W3,W4,W5,W6
	REAL*8 PI,E
	END STRUCTURE
C	GERROR=1, TROUBLE IN SUB.
	STRUCTURE /SLIPM/     !Q
	INTEGER*4 N,SOIL(60)
	REAL*4 X(60),A(60),Y(60),C(60)
	REAL*4 F(60),W(60),P(60),CP(60),CPA(60),WX(60),WY(90)
	END STRUCTURE
	STRUCTURE /ANCHOR/    !ANC
	INTEGER*4 N,NF(20),NL(20),NID(20)
	REAL*4 L(20),TA(20),XF(20),YF(20),XL(20),YL(20),A(20)
	END STRUCTURE
	STRUCTURE /FORCE/   !FC
	REAL*4 G(60),FN(60),FT(60),GN(60),GT(60)
	END STRUCTURE
	STRUCTURE /SFLD/     !SL
	INTEGER*4 ID(60)
	REAL*4  QX(60),QY(60)
	END STRUCTURE
	STRUCTURE /PILE/    !PL
	INTEGER*4 N,NA,NF(20),NL(20),NID(20)
	REAL*4 L(20),TA(20),XF(20),YF(20),XL(20),YL(20),A(20)
	REAL*4 H(20),PA(20),C(20),F(20)
	END STRUCTURE
	RECORD/SFLD/SL
	COMMON/A10/SL
	RECORD/SLIPW/S
	COMMON/A5/S
	RECORD/NODE/N
	COMMON/A6/N
	RECORD/OPTION/O
	COMMON/A7/O
	RECORD/SLIPB/B
	COMMON/A8/B
	INTEGER*4 TENSION,OFOS
	COMMON/CONPAR/UWOW,WLOS,COEQ,TENSION,OFOS
	RECORD/SLIPM/Q
	COMMON/A9/Q
	INTEGER SLOAD,SLINE(15)
	COMMON/SFLOAD1/SLOAD,SLINE,VALUEL(15),VALUER(15),QANGLE(15)
	RECORD/FORCE/FC
	COMMON/CFC/FC
	RECORD/ANCHOR/ANC
	COMMON/ACR/ANC
	RECORD/PILE/PL
	COMMON/PLR/PL