e4.for

边坡稳定上限解程序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