344土质边坡稳定分析一原理·方法·程序 AI(14),滑裂面上控制点处四个与样条插值有关的数,K,K2,K3,K4; A(),滑裂面上控制点处的一阶导数值m B(D),一维数组,工作单元 2.源程序 SUBROUTINE SPI(LS, LSA, X2, Y2) DIMENSION KO2(50),X1(50,Y1(50),AI(50,4)LSA(60),A(50),B(50) DIMENSION X2(80),Y2(80) COMMON /A23/NS1, KQ2, X1, Y1/A13/DD, IWRl, IWR2, IWR3 COMMON/WALL/TWALL, G WALL, HMW,EWALL, ETA INTEGER4 OPTION(6),OP1(6) COMMON/OPP/OPTION. OPI OPEN(S, FILE STATUS=UNKNOWN OPEN(6, FILE STATUS=UNKNOWN READ(S NSI DO=LNSI READO, )KQ2), Xl,YI) ENDDO 如果滑裂面上控制点总数NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的xy值 LSAFO ENDDO IF(LS.NEO)READ(, (LSA I=l, LS) sA(S+IFNS C 滑裂面上第NS个控制点对应的土条侧面边界线的编号 Xl(①),Yl(),滑裂面上第I个控制点的x,y值 N1=0 DO 100IL=LLS+ LI=NI+1 N2=NI-1 IF(LSEQ0GOTO83!滑裂面上没有直线段 IF(NIGE NS1)GOTO 40 AI(Nl,1=0.!计算直线段左端点处K1,K2,K,K4 AI(Nl,3)=Y1(N1+1}Y1(N1) AI(N1, 4FYI(ND) 40IF(N1L1EQ0GOTO100!滑裂面上没有直线段 83B(Ll=1.5°(Y1(L+1Y1(L1)/(X1(L1+1)-Xl(L1) A(L1)=0.5!插值区间左端点处的一阶导数值 CN=3.0*(Y1(Nl)Y1(N2)(X1(Nl}X1(N2) IF(N1-L1EQ.1)GOTO55!插值区间只具有左右两个端点 L2=L1+1 DO 50 J=L2 N2 HI=XIO)(-1) H=XIO+1)-XI) AF=H1/(H+H1)!计算α值 BT=3.*(1.-AF)°Y1(J-Y1(J-1)MHI+AF(Y1(J+1}Y1)/H) 计算β值 (JAF(2.+(1.-AF)*A(J-)!计算插值区间内插值节点处的一阶导数值m 50B(J=(BT-(1.-AF)*B(J-1)/(2+(1.-AF)A(J-1) 5A(NI=(CN-B(N2)/(2.+A(N2) A(N2=A(N2)*A(N2+1)+B(N2)!计算插值区间右端点处的一阶导数值mn344 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 AI(I,4) 滑裂面上控制点处四个与样条插值有关的数 κ1, κ2, κ3, κ4 A(I) 滑裂面上控制点处的一阶导数值m j B(I) 一维数组 工作单元 2. 源程序 C SUBROUTINE SP1(LS,LSA,X2,Y2) DIMENSION KQ2(50),X1(50),Y1(50),AI(50,4),LSA(60),A(50),B(50) DIMENSION X2(80),Y2(80) COMMON /A23/NS1,KQ2,X1,Y1/A13/IDD,IWR1,IWR2,IWR3 COMMON/WALL/IWALL,GWALL,HMW,EWALL,ETA INTEGER*4 OPTION(6),OP1(6) COMMON/OPP/OPTION,OP1 OPEN(5,FILE=' ',STATUS='UNKNOWN') OPEN(6,FILE=' ',STATUS='UNKNOWN') READ(5,*)NS1 DO I=1,NS1 READ(5,*)KQ2(I),X1(I),Y1(I) ENDDO C 如果滑裂面上控制点总数 NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的 x,y 值; READ(5,*)LS DO I=1,NS1-1 LSA(I)=0 ENDDO IF(LS.NE.0)READ(5,*)(LSA(I),I=1,LS) LSA(LS+1)=NS1 C 滑裂面上第 NS1 个控制点对应的土条侧面边界线的编号 C X1(I) Y1(I) 滑裂面上第 I 个控制点的 x y 值 N1=0 DO 100 IL=1,LS+1 L1=N1+1 N1=LSA(IL) N2=N1-1 IF(LS.EQ.0)GOTO 83 ! 滑裂面上没有直线段 IF(N1.GE.NS1) GOTO 40 AI(N1,1)=0. ! 计算直线段左端点处κ1,κ2,κ3,κ4 AI(N1,2)=0. AI(N1,3)=Y1(N1+1)-Y1(N1) AI(N1,4)=Y1(N1) 40 IF(N1-L1.EQ.0)GOTO 100 ! 滑裂面上没有直线段 83 B(L1)=1.5*(Y1(L1+1)-Y1(L1))/(X1(L1+1)-X1(L1)) A(L1)=-0.5 ! 插值区间左端点处的一阶导数值 CN=3.0*(Y1(N1)-Y1(N2))/(X1(N1)-X1(N2)) IF(N1-L1.EQ.1)GOTO 55 ! 插值区间只具有左右两个端点 L2=L1+1 DO 50 J=L2,N2 H1=X1(J)-X1(J-1) H=X1(J+1)-X1(J) AF=H1/(H+H1) ! 计算α值 BT=3.*((1.-AF)*(Y1(J)-Y1(J-1))/H1+AF*(Y1(J+1)-Y1(J))/H) !计算β值 A(J)=-AF/(2.+(1.-AF)*A(J-1)) ! 计算插值区间内插值节点处的一阶导数值 mj 50 B(J)=(BT-(1.-AF)*B(J-1))/(2.+(1.-AF)*A(J-1)) 55 A(N1)=(CN-B(N2))/(2.+A(N2)) 60 A(N2)=A(N2)*A(N2+1)+B(N2) ! 计算插值区间右端点处的一阶导数值 mn