FORTRAN语言潮流计算程序

$DEBUG$LARGE c for 2002 studensimplicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1 /BNDA/ZB(80,5),ZN(50,7),PVA(30,2) /GBDA/GBII(50,4),1 GBIJ(80,3) /FACT/JF(2,100),DJ(100),UJ(2,2000) OPEN(1,FILE=lf02.res)CALL INPUT(SB,U0,REFA)CALL YMCALL DEVR(U0)DO 18 I=1,NPVK=PVA(I,1)18 V(1,K)=PVA(I,2)VS=V(1,NN)V(1,NN)=VS*DCOS(REFA/180*3.1415926)V(2,NN)=VS*DSIN(REFA/180*3.1415926)KT=027 CALL PQ(1,KT) IF(KT.GT.20) STOPIF(DABS(AXM).LE.ES) GOTO 58CALL JACO(1)CALL SBW(1,1)KT=KT+1GOTO 2758 CONTINUEWRITE(1,*)ITERATING NUMBER: ,KTWRITE(1,7)MAXMUM ERROR,F(1,KKK),F(2,KKK), AT NODE,KKK7 FORMAT(1X,A,2F18.11,A,I5) CALL PQ(2,1) close(1)STOPEND*SUBROUTINE INPUT(SB,U0,REFA)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1 /BNDA/ZB(80,5),ZN(50,7),PVA(30,2) /NNOO/NEW(150),NOLD(50,2) COMMON /GBDA/GBII(50,4),GBIJ(80,3)OPEN(7,FILE=lf02.DAT,STATUS=OLD)REWIND(7)READ(7,*)NN,NB,NPV,NSWRITE(1,*)NN,NB,NPV,NS:WRITE(1,1)NN,NB,NPV,NSWRITE(1,*)ES,U0,SB,REFA:C READ(7,*)ES,U0,SB,REFAU0=1.0SB=1.0REFA=0.0READ(7,*)ESWRITE(1,2)ES,U0,SB,REFADO 8 I=1,NBREAD(7,*)(ZB(I,J),J=1,5)8 CONTINUEDO 18 I=1,NN18 READ(7,*)(ZN(I,J),J=1,5)DO 38 I=1,NPV38 READ(7,*)(PVA(I,J),J=1,2)DO 48 I=1,3*NN48 NEW(I)=0DO 58 I=1,NBK=DABS(ZB(I,1)L=DABS(ZB(I,2)NEW(K)=-158 NEW(L)=-1DO 68 I=1,NPVK=PVA(I,1)68 NEW(K)=-2NEW(NS)=-3KL=1DO 338 I=1,3*NNIF(NEW(I).NE.-1) GOTO 338NEW(I)=KLKL=KL+1338 CONTINUEDO 348 I=1,3*NNIF(NEW(I).NE.-2) GOTO 348NEW(I)=KLKL=KL+1348 CONTINUENEW(NS)=KLDO 358 I=1,3*NNIF(NEW(I).EQ.0) GOTO 358K=NEW(I)NOLD(K,1)=I358 CONTINUE DO 78 I=1,NBK=ABS(ZB(I,1) L=ABS(ZB(I,2)ZB(I,1)=ZB(I,1)/K*NEW(K)78 ZB(I,2)=ZB(I,2)/L*NEW(L) DO 98 I=1,NNK=DABS(ZN(I,1)98 ZN(I,1)=NEW(K) DO 128 I=1,NPVK=PVA(I,1)128 PVA(I,1)=NEW(K) DO 148 I=1,NB K=DABS(ZB(I,1) L=DABS(ZB(I,2) IF(K.LE.L) GOTO 148 T=ZB(I,1) ZB(I,1)=ZB(I,2) ZB(I,2)=T148 CONTINUE CALL REGDA(ZB,80,5,2,NB) CALL REGDA(ZN,50,7,1,NN) CALL REGDA(PVA,30,2,1,NPV)C RETURN WRITE(1,*) (ZB(I,J),J=1,5):DO 28 I=1,NB WRITE(1,2)(ZB(I,J),J=1,5)28 CONTINUE WRITE(1,*) (ZN(I,J),J=1,5): DO 168 I=1,NN168 WRITE(1,2)(ZN(I,J),J=1,5) WRITE(1,*) (PVA(I,J),J=1,2): DO 188 I=1,NPV188 WRITE(1,2)(PVA(I,J),J=1,2)1 FORMAT(10I5)2 FORMAT(6F12.6) END*SUBROUTINE REGDA(ARY,KR,KC,K2,KNB) implicit real*8 (a-h,p-z) real*8 ARY(KR,KC) DO 38 K=1,KNB-1 I=dabs(ARY(K,1)I1=dabs(ARY(K,K2)II=I II1=I1 DO 48 K1=K+1,KNB J=dabs(ARY(K1,1) J1=dabs(ARY(K1,K2) IF(I.GT.J.OR.(I.EQ.J.AND.I1.GT.J1) THEN I=JI1=J1KK=K1 ENDIF48 CONTINUE IF(I.EQ.II.AND.I1.EQ.II1) GOTO 38 DO 58 K1=1,KC AI=ARY(KK,K1)ARY(KK,K1)=ARY(K,K1)58 ARY(K,K1)=AI38 CONTINUEEND*SUBROUTINE DEVR(U0) implicit real*8 (a-h,p-z) REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F 1 /BNDA/ZB(80,5),ZN(50,7),PVA(30,2)DO 118 I=1,NN J=ZN(I,1) V(1,J)=U0 V(2,J)=0.D0 P(1,J)=ZN(I,2)-ZN(I,4) P(2,J)=ZN(I,3)-ZN(I,5)C WRITE(1,(F5.0,6F12.6) (ZN(I,K),K=1,5),P(1,J),P(2,J) 118 CONTINUEEND*SUBROUTINE YM implicit real*8 (a-h,p-z) COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK COMMON /GBDA/GBII(50,4),GBIJ(80,3)1 /BNDA/ZB(80,5),ZN(50,7),PVA(30,2) DO 8 I=1,NNDO 8 J=1,48 GBII(I,J)=0.D0DO 58 I=1,NB58 GBIJ(I,3)=0.D0NC=0N=1IOLD=0JOLD=0DO 18 K=1,NBI=ZB(K,1)J=ZB(K,2)AK=ZB(K,5)TI=1.0TJ=1.0IF(I.LT.0) TI=AKIF(J.LT.0) TJ=AKIF(I.LT.0.OR.J.LT.0) AK=0I=dabs(ZB(K,1)J=dabs(ZB(K,2)IF(I.EQ.J) THENNC=NC+1GBII(I,2)=GBII(I,2)+ZB(K,5)GOTO 18ELSEGIJ=ZB(K,3)BIJ=ZB(K,4)RR=GIJ*GIJ+BIJ*BIJGIJ=GIJ/RRBIJ=-BIJ/RRGBII(I,1)=GBII(I,1)+GIJ/TI/TIGBII(I,2)=GBII(I,2)+BIJ/TI/TI+AKGBII(J,1)=GBII(J,1)+GIJ/TJ/TJGBII(J,2)=GBII(J,2)+BIJ/TJ/TJ+AKGBII(I,3)=GBII(I,3)+1IF(I.EQ.IOLD.AND.J.EQ.JOLD) THENN=N-1GBII(I,3)=GBII(I,3)-1ENDIFGBIJ(N,1)=GBIJ(N,1)-GIJ/TI/TJGBIJ(N,2)=GBIJ(N,2)-BIJ/TI/TJGBIJ(N,3)=JENDIFIF(I.NE.J) N=N+1IOLD=IJOLD=J18 CONTINUEGBII(1,4)=1.DO 38 I=1,NN38 GBII(I+1,4)=GBII(I,3)+GBII(I,4)C RETURN15 FORMAT(3X,6F12.6)WRITE(1,*) (GBII(I,J),J=1,4):DO 108 I=1,NNWRITE(1,15)(GBII(I,J),J=1,4)108 CONTINUEKI=GBII(NN+1,4)-1 WRITE(1,*) (GBIJ(I,J),J=1,3):DO 118 I=1,KI118 WRITE(1,15)(GBIJ(I,J),J=1,3)END*SUBROUTINE PQ(III,KT) implicit real*8 (a-h,p-z) REAL*8 V(2,50),P(2,50),F(2,50) COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F1/GBDA/GBII(50,4),GBIJ(80,3) /NNOO/NEW(150),NOLD(50,2)1/BNDA/ZB(80,5),ZN(50,7),PVA(30,2)IF(III.EQ.2) GOTO 30DO 18 I=1,NNF(1,I)=0.D018 F(2,I)=0.D0DO 28 I=1,NN ei=v(1,i) fi=v(2,i) VI=EI*EI+FI*FIF(1,I)=F(1,I)+VI*GBII(I,1)F(2,I)=F(2,I)-VI*GBII(I,2)IF(I.EQ.NN) GOTO 28M=GBII(I,4) N=GBII(I+1,4)-1 DO 38 L=M,N J=GBIJ(L,3)EJ=V(1,J)FJ=V(2,J)GIJ=GBIJ(L,1)BIJ=GBIJ(L,2)F(1,I)=F(1,I)+EI*(GIJ*EJ-BIJ*FJ)+FI*(GIJ*FJ+BIJ*EJ)F(2,I)=F(2,I)+FI*(GIJ*EJ-BIJ*FJ)-EI*(GIJ*FJ+BIJ*EJ)F(1,J)=F(1,J)+EJ*(GIJ*EI-BIJ*FI)+FJ*(GIJ*FI+BIJ*EI)F(2,J)=F(2,J)+FJ*(GIJ*EI-BIJ*FI)-EJ*(GIJ*FI+BIJ*EI)38 CONTINUE28 CONTINUEL=1AXM=0.D0KKK=0DO 48 I=1,NN-1F(1,I)=P(1,I)-F(1,I)F(2,I)=P(2,I)-F(2,I)IPV=PVA(L,1)IF(I.EQ.IPV) THENP(2,I)=P(2,I)-F(2,I)EI=V(1,I)FI=V(2,I)F(2,I)=PVA(L,2)*PVA(L,2)-EI*EI-FI*FIL=L+1ENDIFIF(DABS(F(1,I).GT.DABS(AXM) THENAXM=F(1,I)KKK=IENDIFIF(DABS(F(2,I).GT.DABS(AXM) THENAXM=F(2,I)KKK=IENDIF48 CONTINUEP(1,NN)=F(1,NN)P(2,NN)=F(2,NN)F(1,NN)=0.D0F(2,NN)=0.D0WRITE(1,5)NO.=,KT,: ERR=,F(1,KKK)1 ,F(2,KKK),AT BUS:,NOLD(KKK,1) return5 FORMAT(1X,A,I2,A,2F18.11,5X,A,I4) WRITE(1,*)F(I)-DP(I) goto 5030 CONTINUEDO 58 I=1,NNAI=V(1,I)BI=V(2,I)F(1,I)=DSQRT(AI*AI+BI*BI)BI=BI/AIF(2,I)=dATAN(BI)58 CONTINUEWRITE(1,*)P(I)=Pi+j Qi:WRITE(1,15)(P(1,NEW(I),P(2,NEW(I),I=1,NN) WRITE(1,*)Vi=Ei+jFiWRITE(1,17)(F(1,NEW(KL),F(2,NEW(KL),KL=1,NN) DO 68 I=1,NNF(2,I)=F(2,I)*180./3.141592668 CONTINUEWRITE(1,*)V(I)=Vi,angle:50 WRITE(1,16)(F(1,NEW(KL),F(2,NEW(KL),KL=1,NN) IF(III.EQ.1) WRITE(1,*)V(I)=Ei+j Fi:IF(III.EQ.1) WRITE(1,15)(V(1,NEW(I),V(2,NEW(I),I=1,NN)15 FORMAT(1X,8F9.4)16 FORMAT(1X,8F9.4)17 FORMAT(1X,4(F9.4,F9.6)END*SUBROUTINE JACO(KK) implicit real*8 (a-h,p-z) real*8 H(100),E(100) REAL*8 V(2,50),P(2,50),F(2,50) COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,F COMMON /GBDA/GBII(50,4),GBIJ(80,3)1/BNDA/ ZB(80,5),ZN(50,7),PVA(30,2)1/FACT/ JF(2,100),DJ(100),UJ(2,2000)N2=2*NNII1=3II2=3JF(1,1)=2JF(2,1)=2DO 408 K=1,NN-1RXI=0.D0RXJ=0.D0DO 18 I=1,N2H(I)=0.D018 E(I)=0.D0L=0DO 28 I=1,NPVKPV=PVA(I,1)IF(KPV.NE.K) GOTO 28L=1GOTO 3028 CONTINUE30 R=0.D0X=0.D0DO 38 I=1,K-1M=GBII(I,4)N=GBII(I+1,4)-1DO 48 J=M,NKJ=GBIJ(J,3)IF(KJ.NE.K) GOTO 48CALL JAA(H,E,I,J,K,N2,L,KK) RXI=RXI+GBIJ(J,1)*V(1,I)-GBIJ(J,2)*V(2,I) RXJ=RXJ+GBIJ(J,1)*V(2,I)+GBIJ(J,2)*V(1,I) GOTO 3848 CONTINUE38 CONTINUEM=GBII(K,4)N=GBII(K+1,4)-1DO 58 J=M,NI=GBIJ(J,3)CALL JAA(H,E,I,J,K,N2,L,KK)RXI=RXI+GBIJ(J,1)*V(1,I)-GBIJ(J,2)*V(2,I) RXJ=RXJ+GBIJ(J,1)*V(2,I)+GBIJ(J,2)*V(1,I)58 CONTINUERXI=RXI+GBII(K,1)*V(1,K)-GBII(K,2)*V(2,K)RXJ=RXJ+GBII(K,1)*V(2,K)+GBII(K,2)*V(1,K)GII=GBII(K,1)BII=GBII(K,2)EI=V(1,K)FI=V(2,K)J=2*KH(J-1)=BII*EI-GII*FI-RXJH(J)=-GII*EI-BII*FI-RXIIF(L.NE.1.OR.KK.NE.1) THENE(J-1)=GII*EI+BII*FI-RXIE(J)=BII*EI-GII*FI+RXJELSEE(J-1)=-2.*FIE(J)=-2.*EIENDIFR=F(1,K)X=F(2,K)IF(KK.EQ.2) THENR=0.D0X=0.D0ENDIFH(N2-1)=RE(N2-1)=XI1=2*K-2CALL JBB(H,N2,II1,II2,I1,KK)I1=2*K-1CALL JBB(E,N2,II1,II2,I1,KK)408 CONTINUERETURN5 FORMAT(3X,8I9)6 FORMAT(3X,8F9.4)WRITE(1,*)JF,UJWRITE(1,5)(JF(1,KL),KL=1,10) WRITE(1,6)(UJ(1,KL),KL=1,50) WRITE(1,*) Factor table pointer is,ii1END*SUBROUTINE JAA(H,E,I,J,K,N2,L,KK) implicit real*8 (a-h,p-z) real*8 H(N2),E(N2)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,FCOMMON /GBDA/GBII(50,4),GBIJ(80,3)G=GBIJ(J,1)B=GBIJ(J,2)MM=2*IH(MM-1)=B*V(1,K)-G*V(2,K) H(MM)=-G*V(1,K)-B*V(2,K)IF (L.EQ.1.AND.KK.EQ.1) GOTO 10 E(MM-1)=-H(MM)E(MM)=H(MM-1)10 CONTINUERETURNEND*SUBROUTINE JBB(H,N2,II1,II2,K,KK) implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50),H(N2)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,FCOMMON /FACT/ JF(2,100),DJ(100),UJ(2,2000)K1=0IF(II1.EQ.3) GOTO 40DO 38 I=1,KIF (dabs(H(I).LT.1.0E-15) GOTO 38 IF(KK.EQ.2) THEN UJ(2,II2)=H(I)UJ(2,II2+1)=III2=II2+2K1=K1+1IF(II2.GE.2000-10) WRITE(1,*) Factor table pointer is,ii2IF(II2.GE.2000) WRITE(1,*) Add UJ(2,2000) IF(II2.GE.2000) stopENDIFM=JF(1,I)N=JF(1,I+1)-1DO 18 J=M,NLL=UJ(1,2*J)H(LL)=H(LL)-H(I)*UJ(1,2*J-1)18 CONTINUE38 CONTINUE40 K=K+1S=H(K)IF(KK.EQ.2) THENJF(2,K+1)=JF(2,K)+K1DJ(K)=SWRITE(1,*) k,sENDIFK1=0DO 48 I=K+1,N2-1IF (dabs(H(I).LT.1.0E-15) GOTO 48UJ(1,II1)=H(I)/SUJ(1,II1+1)=IK1=K1+1II1=II1+248 CONTINUEIF(II1.GE.2000-10) WRITE(1,*) Factor table pointer is,ii1IF(II1.GE.2000) WRITE(1,*) Add UJ(1,2000) IF(II1.GE.2000) stopJF(1,K+1)=JF(1,K)+K1RETURNEND*SUBROUTINE SBW(KK,KK1)implicit real*8 (a-h,p-z)REAL*8 V(2,50),P(2,50),F(2,50)COMMON NN,NB,NG,NL,NPV,NGL,NT,ES,AXM,KKK /VPF/V,P,FCOMMON /FACT/ JF(2,100),DJ(100),UJ(2,2000)1/BNDA/ZB(80,5),ZN(50,7),PVA(30,2) /BH/ H(100)IF(KK.EQ.2) THENDO 68 I=KK1,2*NN-2K=JF(2,I)DO 78 J=K,JF(2,I+1)-1LL=UJ(2,2*J)78H(I)=H(I)-H(LL)*UJ(2,2*J-1)H(I)=H(I)/DJ(I)68 CONTINUEELSEDO 18 I=1,2*NN18 H(I)=0.D0ENDIFK1=2*NN-2DO 38 I=K1,1,-1K=JF(1,I)K2=JF(1,I+1)-1LL=UJ(1,2*K2)IF(LL.NE.2*NN-1.or.k.gt.k2) THENK2=K2+1GOTO 30ELSEH(I)=UJ(1,2*K2-1)ENDIF30 DO 28 J=K,K2-1L=UJ(1,2*J)28 H(I)=H(I)-UJ(1,2*J-1)*H(L)38 CONTINUEIF(KK.EQ.2) GOTO 60DO 48 I=1,NN-1K=2*IV(1,I)=V(1,I)-H(K)V(2,I)=V(2,I)-H(K-1)48 CONTINUE60 CONTINUEEND*RESNN,NB,NPV,NS:5711ES,U0,SB,REFA:0.0002001.0000001.0000000.000000(ZB(I,J),J=1,5):1.0000002.0000000.0600000.1800000.0000001.0000003.0000000.0600000.1800000.0000001.0000004.0000000.0400000.1200000.0000001.0000005.0000000.0200000.0600000.0000002.0000003.0000000.0100000.0300000.0000002.0000005.0000000.0800000.2400000.0000003.0000004.0000000.0800000.2400000.000000(ZN(I,J),J=1,5):1.0000000.2000000.2000000.0000000.0000002.0000000.0000000.0000000.4500000.1500003.0000000.0000000.0000000.4000000.0500004.0000000.0000000.0000000.6000000.1000005.0000000.0000000.0000000.0000000.000000(PVA(I,J),J=1,2):5.000000 1.060000(GBII(I,J),J=1,4):10.833333-32.5000004.0000001.00000012.0 -38.7500002.0000005.00000012.0 -38.7500001.0000007.0000003.0 -11.2500000.0000008.0000006.250000-18.7500000.0000008.000000(GBIJ(I,J),J=1,3):-1.6666675.0000002.000000-1.6666675.0000003.000000-2.5000007.5000004.000000-5.00000015.0000005.000000-10.00000030.0000003.000000-1.2500003.7500005.000000-1.2500003.7500004.000000NO.= 0:ERR=0.500000000001.10000000000AT BUS:2NO.= 1:ERR=-0.07704425659-0.02203816492AT BUS:2NO.= 2:ERR=-0.00000341934-0.00084234440AT BUS:5NO.= 3:ERR=-0.00000002401-0.00000014864AT BUS:5ITERATING NUMBER: 3MAXMUM ERROR -0.00000002401 -0.00000014864 AT NODE 4 P(I)=Pi+j Qi:1.2982 0.2445 0.2000 0.2000 -0.4500 -0.1500-0.6000 -0.1000Vi=Ei+jFi1.0600 0.000000 1.0365-0.046069 1.0088-0.0839071.0016-0.104414V(I)=Vi,angle:1.0600 0.0000 1.0365 -2.6396 1.0088 -4.80751.0016 -5.9825-0.4000 -0.05001.0073-0.0896091.0073 -5.1342DATE5,7,1,10.0002,1.00,1,0.01,2,0.02,0.06,0.001,3,0.08,0.24,0.00 2,3,0.06,0.18,0.00 2,4,0.06,0.18,0.00 2,5,0.04,0.12,0.00 3,4,0.01,0.03,0.00 4,5,0.08,0.24,0.00 1,0.0,0.0,0.0,0.0 2,0.2,0.2,0.0,0.0 3,0.0,0.0,0.45,0.15 4,0.0,0.0,0.4,0.05 5,0.0,0.0,0.6,0.1 1,1.06610.83333312.91666712.9166673.7500006.2500000.000000-32.500000-38.750000-38.750000-11.250000-18.7500000.000000-1.666667-1.666667-2.500000-5.000000-10.000000-1.250000-1.2500005.0000005.0000007.50000015.00000030.0000003.7500003.7500001.0000005.0000007.0000008.0000008.0000008.0000002.0000003.0000004.0000005.0000003.0000005.0000004.000000NN,NB,NPV,NS:5711ES,U0,SB,REFA:.0002001.0000001.000000 .000000(GBII(I,J),J=1,3):10.833333-32.5000001.00000012.916667-38.7500005.00000012.916667-38.7500007.0000003.750000-11.2500008.0000006.250000-18.7500008.000000.000000.0000008.000000(GBIJ(I,J),J=1,3):-1.6666675.0000002.000000-1.6666675.0000003.000000-2.5000007.5000004.000000-5.00000015.0000005.000000-10.00000030.0000003.000000-1.2500003.7500005.000000-1.2500003.7500004.000000(ZB(I,J),J=1,5):1.0000002.000000.060000.180000.0000001.0000003.000000.060000.180000.0000001.0000004.000000.040000.120000.0000001.0000005.000000.020000.060000.0000002.0000003.000000.010000.030000.0000002.0000005.000000.080000.240000.0000003.0000004.000000.080000.240000.000000(ZN(I,J),J=1,5):1.000000.200000.200000.000000.0000002.000000.000000.000000.450000.1500003.000000.000000.000000.400000.0500004.000000.000000.000000.600000.1000005.000000.000000.000000.000000.000000(PVA(I,J),J=1,2):5.000000 1.060000NO.= 0:ERR=.500001000001.10000000000AT BUS:2NO.= 1:ERR=-.07704425188-.02203816770AT BUS:2NO.= 2:ERR=-.00000341921-.00084234361AT BUS:5NO.= 3:ERR=-.00000002401-.00000014864AT BUS:5ITERATING NUMBER: 3MAXMUM ERROR-.00000002401-.00000014864 AT NODE 4P(I)=Pi+j Qi: