环状网平差程序 那位高手能帮忙用c改写一下 谢谢


PROGRAM GS
DIMENSION I(200),J(200),IJ(200),D(200),Q(200),L(200),S(200),
* HF(200),HFZ(50),DQ(50),SQ(200),SQZ(50),A(200),B(200),NO(10)
REAL L
c I(K)-K管段属于小环的环号, J(K)--K管段属于大环的环号,IJ(K)--K管段内流量在小环的方向
C D(K)，Q(K), L(K)--K管段的管径、流量、管长
C IG-- 管段的个数， NK---水源个数，NO(I)--水源编号，
C 水源可写成通式：H=A+B*Q^2 对于水泵 H=Hb-B*Q^2,Hb为水泵的扬程，对于水塔 H=Ht 为水塔高度

OPEN(1,FILE='GSQ.DAT')
READ(1,*)IH,IG,AS,(I(K),J(K),IJ(K),Q(K),D(K),L(K),K=1,IG)
READ(1,*)NK
IF(NK.NE.1) THEN
READ(1,*)(NO(II),A(ABS(NO(II))),B(ABS(NO(II))),II=1,NK)
END IF
CLOSE(1)

c 管段循环
DO K=1,IG
IJ(K)=1
J(K)=ABS(J(K))
D(K)=D(K)/1000
Q(K)=Q(K)/1000*SIGN(1,I(K))
END DO

IV=1
100 DO K=1,IG
IF ((D(K).NE.0.).AND.(Q(K).NE.0.)) THEN
V=(4/3.14159)*Q(K)/D(K)/D(K)
IF(ABS(V).GE.1.2) THEN
S(K)=0.001735*L(K)/D(K)**5.3
ELSE
S(K)=0.00148/D(K)**5.3*(1+0.681*D(K)*D(K)/ABS(Q(K)))
* **0.3*L(K)
END IF
END IF
c 实管段，计算S*Q 和管段水头损失Hf
IF(D(K).GT.0.00001) THEN
SQ(K)=ABS(S(K)*Q(K))
HF(K)=SQ(K)*Q(K)
ELSE
c 虚管段
SQ(K)=ABS(B(K)*Q(K))
DO II=1,NK
IF(ABS(NO(II)).EQ.K) THEN
HF(K)=(A(K)-B(K)*Q(K)*Q(K))*SIGN(1,NO(II))*SIGN(1,I(K))
GOTO 111
END IF
END DO
111 END IF
END DO

c 环循环，计算环的总S*Q之和SQA，水头损失之和HFA
DO N=1,IH
SQA=0.
HFA=0.
DO K=1,IG
IF(ABS(I(K)).EQ.N) THEN
C 所属小环
HFA=HFA+HF(K)
SQA=SQA+2*SQ(K)
END IF
IF(J(K).EQ.N) THEN
C 所属大环
HFA=HFA-HF(K)
SQA=SQA+2*SQ(K)
END IF
END DO
HFZ(N)=HFA
SQZ(N)=SQA
C 计算修正流量
DQ(N)=-1*HFZ(N)/SQZ(N)
END DO

C 对各个管段进行流量修正
DO K=1,IG
IF(J(K).EQ.0) THEN
Q(K)=Q(K)+DQ(ABS(I(K)))/IJ(K)
ELSE
Q(K)=Q(K)+(DQ(ABS(I(K)))-DQ(ABS(J(K))))/IJ(K)
END IF
END DO

C 对各环闭合差进行判断
DO N=1,IH
IF(ABS(HFZ(N)).GT.AS) THEN
IV=IV+1
GOTO 100
END IF
END DO

C 输出结果
OPEN(2,FILE='result.DAT')
WRITE(2,*)'--------------------------------------------------',
* '-------------'
WRITE(2,*)' NO D L S HF(m) ',
* ' Q(l/s)'
WRITE(2,*)'--------------------------------------------------',
* '-------------'
WRITE(2,1)(K,D(K)*1000,L(K),S(K),HF(K),Q(K)*1000,K=1,IG)
1 FORMAT(1X,I3,4F10.1,F12.3)
WRITE(2,*)'--------------------------------------------------',
* '-------------'
WRITE(2,*)'DIE DAI :','IV=',IV
WRITE(2,*)'==================================================',
* '============================='
WRITE(2,2)(N,N=1,IH)
2 FORMAT(1X,' NO !',10I8)
WRITE(2,*)'--------------------------------------------------',
* '-----------------------------'
WRITE(2,3)(HFZ(N),N=1,IH)
3 FORMAT(1X,' HFZ !',10F8.3)
WRITE(2,*)'--------------------------------------------------',
* '-----------------------------'
WRITE(2,4)(DQ(N),N=1,IH)
4 FORMAT(1X,' DQ !',10F8.3)
WRITE(2,*)'--------------------------------------------------',
* '-----------------------------'
CLOSE(2)
END