計算傳熱學課程報告一、問題概述：有限單元法是上個世紀五、六十年代首先在力學中發展起來的數值計算方法，由于它是基于變分原理，理論基礎統一，對于復雜邊界的適應程度比較好，所以很快的在其它領域得到運用，其中就包括了在傳熱學中的運用。本次計算傳熱學的課程就是對有限單元法在傳熱學中運用的一個學習與練習。有限單元法處理問題的步驟，首先是建立有限元模型也即是將問題離散化，它的主要步驟之一就是將要計算的物體進行有限元的劃分；第二步，進行單元分析也就是將變分原理運用到問題的方程與單元中，形成單元剛度矩陣；第三步，進行整體剛度矩陣的組集；最后就是引入邊界條件進行求解的過程。在計算傳熱學的課程中，主要完成了兩個任務：第一，是將一個比較復雜的活塞進行了網格劃分，并編譯成一個通用性比較好的程序。第二，在前一個程序的基礎上，加入計算過程，運用焓法，對一個比較簡單的平面相變問題進行了計算。二、劃分單元網格：劃分單元網格是將問題進行有限元法分析的基礎，但是如果在圖紙上進行手工的單元劃分，不但繁瑣、容易出錯，而且也不利于進一步計算程序的利用。因此有必要編輯一個程序，以自動完成劃分網格的目的。網格的自動劃分必須遵循以下的幾條規則：(1).要嚴格區分邊界單元與內部單元，并且嚴格區分邊界單元不同的組；(2).單元標號必須先標志內部單元，然后依次標志第一類邊界條件，第二類邊界條件，第三類邊界條件，如果同一類邊界條件中有不同的組，那么也必須嚴格先劃分第一組，然后第二組，第三組；(3). 對于邊界單元，每一個邊界單元必須只有一條邊在邊界上，而且為了程序的簡單，一般是j,m邊作為邊界；(4).內部單元節點標號必須遵循逆時針方向的規則；(5). 一個單元中只能有一種材料組成。遵循以上的規則，用FORTRAN 90編制了一個對形狀比較復雜的活塞的網格劃分，由于在編制過程中考慮了多種情況，所以這個程序有比較好的通用性，只需要輸入不同的數據，程序也可以對許多其它情況進行劃分。需要指出的是，由于FORTRAN 90程序對于制圖功能比較弱，所以下面的圖是用VB 6.0的程序做出的，由于該網格劃分程序集成了后續對第一類邊界條件和第三類邊界條件的焓法計算程序，故該程序源代碼將在最后統一給出。網格劃分的結果如圖(1)。需要輸入的初始數據主要有：邊界單元分組總數、邊界單元分組中前一組的最后一個單元號、各組邊界單元節點數、各邊界單元邊界節點號、每一條層線的左右端點、坐標，每一層單元劃分所屬的類型。圖（）三、焓法有限單元法原理：帶有相變的傳熱問題，又被稱為斯蒂芬問題，在冶金、鑄造、建筑、冷凍、航天和醫療等領域有著廣泛的運用。由于是這類問題存在的相變過程，使得求解區域中存在著一個隨著時間移動的固液或者固氣界面，這一界面使得這種問題的求解非常困難。一般說來，處理這種問題有兩種不同的思路，一種思路是，首先著眼于相變界面的求解，確定相變界面以后再分別處理固相或者液相的溫度分布；另一種思路是將該問題看作是“單相”的非線性導熱問題，首先確定整個求解區域上的溫度分布或者焓的分布，然后把達到相變溫度的位置定為相變界面。在實際運用中，后一種思路比較簡單和實用而得到了廣泛的使用；目前后一種處理方式中主要有焓法和顯熱容法，這里主要討論焓法。焓法有限元法的主要思路是，不把溫度作為求解的變量，而是把焓作為求解變量，因此可以在固相和液相的整個區域建立統一的方程進行求解。然后根據焓與溫度的關系確定整個區域的溫度與相變的界面。焓法的優點就是不需要跟蹤相變的界面而可以對整個區域進行統一的求解。焓法的主要推導如下：在一般過程中，傳熱微分方程由(1)式表示：（）其中為導熱系數，為物體的密度， 為物體的定壓熱容，為物體的溫度，而表示時間。由于焓與溫度存在如下的關系：所以（）式可以化為（）式：（）該式子即是在整個區域上都適用的焓法的基本方程，對于這個方程，在空間上用有限單元法離散，在時間上用向后差分格式離散，就可以得到如下的(3)式：（）其中：整體剛度矩陣 變焓矩陣 時間步長t時刻熱焓時刻熱焓t時刻右端列向量值得注意的是，首先使用變分法的使用要直接對焓而不是溫度變分；其次整體剛度矩陣和變焓矩陣是與時間無關的的矩陣，而右端向量則與邊界條件有關。最后焓與溫度的關系由（）式給出：（）四、焓法實例：為了簡單起見，運用上述的焓法對如圖(2)的簡單區域的相變問題進行求解，其條件如下：邊長為80cm60cm，周圍用溫度為500K，鐵水溫度為1833K，比熱為711.62J/kg.k，潛熱為271100J/kg，密度為7800kg/m3，導熱系數為33.5W/m. 相溫度為1790K，使用第三類邊界條件，換熱系數1950W/Km2。圖 （2）圖中1-4，1-15等邊都是第三類邊界，其余單元則作為邊界單元處理。程序源代碼如下： PROGRAM MAIN INTEGER L0,V0,E3,C0,B0,F1,V7,M2,A0,D9,Z0,PANBIE,B3380integer,dimension(:),allocatable:B(:),F(:),W(:),M(:),J(:),I(:),H10(:),H11(:),H12(:),H13(:)90 real,dimension(:),allocatable:P(:),Q(:),X(:),Z(:),Y(:),HH1(:),HH2(:),HH3(:),HH32(:),HH4(:),TT(:) DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE:H9(:)doubleprecision K1(3,3),N1(3,3)doubleprecision,dimension(:,:),allocatable:K2(:,:),N2(:,:),EH(:,:),EH2(:,:)100 REAL,dimension(:),allocatable:C(:),S(:),R(:) REAL KC,CP,MD,TS1 ,TS2,ARF,TF,TM REAL HH33(3) open (1,file=score.dat)70 read(1,*) L0,V0,E3,C0,B0,F175 V7=V0+1 allocate(B(F1),F(B0),W(B0),H9(E3),M(V7),H10(E3),H11(E3),H12(E3),H13(E3),HH1(L0),HH2(L0),HH3(L0),HH32(L0),HH4(L0),TT(L0) allocate(J(V0),P(V0),Q(V0),Y(V0),Z(V0),I(C0),C(C0),S(C0),X(L0),R(L0),k2(L0,L0),n2(L0,L0),EH(L0,L0)110 DO 125 J0=1,B0120 READ(1,*) W(J0),F(J0)125 CONTINUE READ(1,*)(B(A0),A0=1,F1) 150 READ(1,*) (M(A0),A0=1,V7)160 do 190 A0=1,V0165 READ(1,*) J(A0),P(A0),Q(A0),Y(A0),Z(A0)190 CONTINUE195 DO 215 A0=1,C0200 READ(1,*) I(A0),C(A0),S(A0)215 CONTINUE READ(1,*)KC,CP,MD,TS,PANBIE,ARF,TF,TS1,TM READ(1,*)(HH1(I0),I0=1,L0) B33=W(1) DO I0=1,L0 TT(I0)=0 ENDDO DO I0=1,L0 DO J0=1,L0 K2(I0,J0)=0 N2(I0,J0)=0 ENDDO ENDDO DO I0=1,L0 HH3(I0)=0 HH4(I0)=0 ENDDO TS2=0 M2=1250 DO I0=1,V0252 IF (J(I0)=5) GOTO 266254 IF(J(I0)=6) GOTO 262256 IF(J(I0)=7) GOTO 262258 D9=M(I0+1)-1 GO TO 270262 D9=M(I0+1)-5 GO TO 270266 D9=M(I0+1)-4270 DO J0=M(I0),D9272 X(J0)=Y(I0)+(J0-M(I0)*(Z(I0)-Y(I0)/(D9-M(I0)274 R(J0)=P(I0)+(J0-M(I0)*(Q(I0)-P(I0)/(D9-M(I0)276 ENDDO ENDDO290 Z0=0292 DO I0=1,V0-1294 D9=M(I0+1)-M(I0)296 IF(J(I0)/=8) GOTO 302298 CALL JK1(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B,W,H10,H11,H12,H13) GOTO 316302 IF(J(I0)=2) GOTO 314304 IF(J(I0)=4) GOTO 314306 IF(J(I0)=10) GOTO 314308 IF(J(I0)=12) GOTO 314310 CALL JK7(D9,M,I0,V7,B0,Z0,H9,E3,M2,F,F1,B,W,J,V0,H10,H11,H12,H13) GOTO 316314CALL JK4(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B,W,J,V0,H10,H11,H12,H13)316 ENDDO350 DO A0=1,C0352 X(I(A0)=C(A0) R(I(A0)=S(A0)356 ENDDO close(1) open (2,file=out.dat)380 write(2,(coordinate of knots)385 DO A0=1,L0390 write(2,500)A0,X(A0),A0,R(A0)400 ENDDO405 write(2,(INFORMATION OF ELEMENTS)410 DO I0=1,E3415 write(2,600)I0,H10(I0),H11(I0),H12(I0),H13(I0)420 ENDDO close(2) open (3,file=out2.dat) do a0=1,L0 WRITE(3,*)X(A0) WRITE(3,*)R(A0) ENDDO DO I0=1,E3 WRITE (3,700)H10(I0),H11(I0),H12(I0),H13(I0) ENDDO CLOSE(3)500 format(1x,X,I4,=,F16.4,R,I4,=,F16.4)600 FORMAT(1X,H,I4,=,I4,I4,I4,I4)700 FORMAT(1X,I4,I4,I4,I4) open (4,file=out4.dat) DO J0=1,E3 call DANGANG(K1,X,R,H10,H11,H12,J0,L0,KC,E3,CP,N1,PANBIE,ARF,HH33,TF,W,B0,MD,B33) CALL ZHENGGANG(K1,K2,H10,H11,H12,L0,E3,J0,N1,N2,HH33,HH4) ENDDO CALL ZUIHOU(K2,N2,EH,L0) DO WHILE(TS2=500.AND.(TT(7)=1823) THEN DO J0=1,L0 WRITE(4,800)J0,TT(J0) ENDDO ELSE CONTINUE ENDIF DO J0=1,L0 HH1(J0)=HH2(J0) ENDDO TS2=TS2+TS ENDDO CLOSE(4)800 FORMAT(1X,I=,I4,H=,F16.4)425 END SUBROUTINE WENDU(HH2,L0,CP,TM,TT) INTEGER L0 REAL HH2(L0),TT(L0),CP,LS,C,TM LS=271100 C=CP*TM DO I=1,L0 IF(HH2(I)=C+LS) THEN TT(I)=(HH2(I)-LS)/CP ELSE TT(I)=TM ENDIF ENDDO END SUBROUTINE GUASS(EH2,HH32,HH2,L0) INTEGER L0,L REAL HH32(L0),HH2(L0),SU,C(L0) double precision EH2(L0,L0),D(L0,L0),O L=0 SU=0.0 DO K=1,L0-1 M=K S=EH2(K,K) L=K DO WHILE(M=L0-1) M=M+1 S1=ABS(S) S2=EH2(M,K) S3=ABS(S2) IF(S1=S3)THEN L=M S=EH2(M,K) ELSE CONTINUE ENDIF ENDDO M=0 DO J=K,L0 C(J)=EH2(L,J) EH2(L,J)=EH2(K,J) EH2(K,J)=C(J) ENDDO O=HH32(L) HH32(L)=HH32(K) HH32(K)=O DO I=K+1,L0 D(I,K)=EH2(I,K)/EH2(K,K) ENDDO DO I=K+1,L0 DO J=K+1,L0 EH2(I,J)=EH2(I,J)-D(I,K)*EH2(K,J) ENDDO ENDDO DO I=K+1,L0 HH32(I)=HH32(I)-D(I,K)*HH32(K) ENDDO ENDDO S4=EH2(L0,L0) HH2(L0)=HH32(L0)/S4 DO K=L0-1,1,-1 DO J=K+1,L0 SU=SU+EH2(K,J)*HH2(J) ENDDOHH2(K)=(HH32(K)-SU)/EH2(K,K) SU=0 ENDDO END SUBROUTINE ZHENGLI(L0,EH,W,B0,H10,H11,H12,E3,EH2,HH32,HH3,HH1,B33) INTEGER L0,B0,W(B0),B33,E3,H10(E3),H11(E3),H12(E3),J33,M33 REAL HH32(L0),HH1(L0),HH3(L0) double precision EH(L0,L0),EH2(L0,L0) DO I0=1,L0 DO J0=1,L0 EH2(I0,J0)=EH(I0,J0) ENDDO HH32(I0)=HH3(I0) ENDDO DO I0=1,E3 IF(I0=B33) THEN CONTINUE ELSE J33=H11(I0) M33=H12(I0) EH2(J33,J33)=EH(J33,J33)*(10*8) EH2(M33,M33)=EH(M33,M33)*(10*8) HH32(J33)= EH2(J33,J33)*HH1(J33)HH32(M33)=EH2(M33,M33)*HH1(M33) ENDIF ENDDO END SUBROUTINE YOUBIAN(N2,L0,HH1,HH3,PANBIE,HH4) INTEGER L0,PANBIE double precision X,HH(L0) double precision N2(L0,L0) REAL HH1(L0),HH3(L0),HH4(L0) I=0 X=10 DO I=1,L0 HH3(I)=0 HH(I)=0 ENDDO DO I0=1,L0 DO J0=1,L0 HH(I0)=HH(I0)+N2(I0,J0)*HH1(J0) ENDDO ENDDO DO I0=1,L0 HH3(I0)=HH(I0)/X ENDDO DO I0=1,L0 HH3(I0)=HH3(I0)+HH4(I0) ENDDO END SUBROUTINE ZUIHOU(K2,N2,EH,L0) INTEGER L0 double precision K2(L0,L0),N2(L0,L0),EH(L0,L0),EHX(L0,L0) REAL X0 X0=10 DO I0=1,L0 DO J0=1,L0 EHX(I0,J0)=0 ENDDO ENDDO DO I0=1,L0 DO J0=1,L0 EHX(I0,J0)=N2(I0,J0)/X0 ENDDO ENDDO DO I0=1,L0 DO J0=1,L0 EH(I0,J0)=EHX(I0,J0)+K2(I0,J0) ENDDO ENDDO END SUBROUTINE ZHENGGANG(K1,K2,H10,H11,H12,L0,E3,J0,N1,N2,HH33,HH4) INTEGER J0,I33,J33,M33,L0,E3 INTEGER H10(E3),H11(E3),H12(E3) double precision K1(3,3),K2(L0,L0),N1(3,3),N2(L0,L0) REAL HH4(L0),HH33(3) I33=H10(J0) J33=H11(J0) M33=H12(J0) K2(I33,I33)=K2(I33,I33)+K1(1,1) K2(I33,J33)=K2(I33,J33)+K1(1,2) K2(I33,M33)=K2(I33,M33)+K1(1,3) K2(J33,I33)=K2(J33,I33)+K1(2,1) K2(J33,J33)=K2(J33,J33)+K1(2,2) K2(J33,M33)=K2(J33,M33)+K1(2,3) K2(M33,I33)=K2(M33,I33)+K1(3,1) K2(M33,J33)=K2(M33,J33)+K1(3,2) K2(M33,M33)=K2(M33,M33)+K1(3,3) N2(I33,I33)=N2(I33,I33)+N1(1,1) N2(I33,J33)=N2(I33,J33)+N1(1,2) N2(I33,M33)=N2(I33,M33)+N1(1,3) N2(J33,I33)=N2(J33,I33)+N1(2,1) N2(J33,J33)=N2(J33,J33)+N1(2,2) N2(J33,M33)=N2(J33,M33)+N1(2,3) N2(M33,I33)=N2(M33,I33)+N1(3,1) N2(M33,J33)=N2(M33,J33)+N1(3,2) N2(M33,M33)=N2(M33,M33)+N1(3,3) HH4(J33)=HH33(2)+HH4(J33) HH4(M33)=HH33(3)+HH4(M33) END SUBROUTINE DANGANG(K1,X,R,H10,H11,H12,J0,L0,KC,E3,CP,N1,PANBIE,ARF,HH33,TF,W,B0,MD,B33) INTEGER J0,I33,J33,M33,L0,E3,PANBIE,B0,B33 REAL B1,B2,B3,C1,C2,C3,SM,KC,KCF,CP,SI,ARF,MD INTEGER H10(E3),H11(E3),H12(E3),W(B0) REAL X(L0),R(L0),HH33(3) double precision K1(3,3),N1(3,3) I33=H10(J0) J33=H11(J0) M33=H12(J0) B1=R(J33)-R(M33) B2=R(M33)-R(I33) B3=R(I33)-R(J33) C1=X(M33)-X(J33) C2=X(I33)-X(M33) C3=X(J33)-X(I33) SM=0.5*(B1*C2-B2*C1) KCF=KC/(4*SM) SI=SQRT(B1*2+C1*2) IF(J0=B33).and.(J0/=1).AND.(J0/=3).AND.(J0/=5)THEN K1(1,1)=KCF*(B1*2+C1*2) K1(2,2)=KCF*(B2*2+C2*2) K1(3,3)=KCF*(B3*2+C3*2) K1(1,2)=KCF*(B1*B2+C1*C2) K1(2,1)=K1(1,2) K1(1,3)=KCF*(B1*B3+C1*C3) K1(3,1)=K1(1,3) K1(2,3)=KCF*(B2*B3+C2*C3) K1(3,2)=K1(2,3) N1(1,1)=(MD*CP*SM)/6.0 N1(2,2)=N1(1,1) N1(3,3)=N1(1,1) N1(1,2)=(MD*CP*SM)/12.0 N1(1,3)=N1(1,2) N1(2,1)=N1(1,2) N1(2,3)=N1(1,2) N1(3,1)=N1(1,2) N1(3,2)=N1(1,2) HH33(1)=0 HH33(2)=0 HH33(3)=0 ELSE K1(1,1)=KCF*(B1*2+C1*2) K1(2,2)=KCF*(B2*2+C2*2)+(ARF*SI)/3.0 K1(3,3)=KCF*(B3*2+C3*2)+(ARF*SI)/3.0 K1(1,2)=KCF*(B1*B2+C1*C2) K1(2,1)=K1(1,2) K1(2,3)=KCF*(B2*B3+C2*C3)+(ARF*SI)/6.0 K1(3,2)=K1(2,3) K1(1,3)=KCF*(B1*B3+C1*C3) K1(3,1)=K1(1,3) N1(1,1)=(MD*CP*SM)/6.0 N1(2,2)=N1(1,1) N1(3,3)=N1(1,1) N1(1,2)=(MD*CP*SM)/12.0 N1(1,3)=N1(1,2) N1(2,1)=N1(1,2) N1(2,3)=N1(1,2) N1(3,1)=N1(1,2) N1(3,2)=N1(1,2) HH33(2)=(ARF*SI*TF)/2 HH33(3)=HH33(2) HH33(1)=0 ENDIF END SUBROUTINE JK1(D9,M,I0,V7,Z0,H9,E3,M2,B0,F,F1,B,W,H10,H11,H12,H13) INTEGER D9,I0,Z0,V7,J0,Z9,I1,M1,J1,E3,M2,B0,F1 DOUBLE PRECISION H9(E3) INTEGER M(V7),F(B0),B(F1),W(B0),H10(E3),H11(E3),H12(E3),H13(E3)3802 D9=(D9+1)/23804 DO J0=2,D93806 I1=M(I0)+J0-1 J1=M(I0+1)+J0-23810 M1=M(I0)+J0-23812 IF(J0=2)GOTO 38223814 Z9=Z0 CALL CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+13820 GOTO 38243822 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)3824 I1=M(I0)+J0-1 J1=M(I0+1)+J0-13828 M1=M(I0+1)+J0-23830 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)3832 ENDDO3834 DO J0=D9,2*D9-23836 I1=M(I0)+J0-1 J1=M(I0)+J03840 M1=M(I0+1)+J03842 IF(J0=2*D9-2)GOTO 38523844 Z9=Z0 CALL CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+13850 GOTO 38543852 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)3854 I1=M(I0)+J0-1 J1=M(I0+1)+J03858 M1=M(I0+1)+J0-1 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)3862 ENDDO END subroutine CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) INTEGER E3,Z9,I1,J1,M1,M2 DOUBLE PRECISION H9(E3) INTEGER H10(E3),H11(E3),H12(E3),H13(E3)3805 H9(Z9+1)=I1*1E-3+J1*1E-6+M1*1E-9 H10(Z9+1)=I1 H11(Z9+1)=J1 H12(Z9+1)=M1 H13(Z9+1)=M23810 END SUBROUTINE BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13) INTEGER B0,Z9,Z0,A0,Y8,S2,I1,J1,M1,M2,E3,F1 DOUBLE PRECISION H9(E3) INTEGER F(B0),B(F1),W(B0) INTEGER H10(E3),H11(E3),H12(E3),H13(E3)3702 Y8=0 S2=03704 DO A0=1,B03708 CALL MJ(S2,A0,F,B0,Y8,B,F1,H9,E3,I1,J1,M1,M2,W,H10,H11,H12,H13)3710 IF(S2=2) GOTO 37203712 ENDDO3714 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+13720 END SUBROUTINE MJ(S2,A0,F,B0,Y8,B,F1,H9,E3,I1,J1,M1,M2,W,H10,H11,H12,H13) INTEGER S2,E9,B0,A0,Y8,F1,Z9,E3,I1,J1,M1,M2 INTEGER F(B0),B(F1),W(B0) DOUBLE PRECISION H9(E3) INTEGER H10(E3),H11(E3),H12(E3),H13(E3)3602 S2=03604 DO E9=1,F(A0)3608 Y8=Y8+13610 IF (B(Y8)=J1) THEN S2=S2+1 ENDIF IF(B(Y8)=M1) THEN S2=S2+1 ENDIF ENDDO3620 IF(S2/=2)GOTO 36283622 Z9=W(A0) CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13)3626 W(A0)=W(A0)+13628 END subroutine JK7 (D9,M,I0,V7,B0,Z0,H9,E3,M2,F,F1,B,W,J,V0,H10,H11,H12,H13) INTEGER D9,I0,Z0,V7,J0,Z9,I1,M1,J1,E3,M2,B0,F1,V0 DOUBLE PRECISION H9(E3) INTEGER M(V7),F(B0),B(F1),W(B0),J(V0) INTEGER H10(E3),H11(E3),H12(E3),H13(E3)4502 I1=M(I0+1)+1 J1=M(I0+1)4506 M1=M(I0) CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4510 IF(J(I0)=5) D9=D9-54512 IF(J(I0)=6) D9=D9-54518 IF(J(I0)/=7) GOTO 45224520 D9=D9-44522 DO J0=2,D94524 IF(J(I0)=1) GOTO 45324526 IF(J(I0)=7) GOTO 45324528 IF(J(I0)=9) GOTO 45324530 GOTO 45344532 IF(J0=D9) GOTO 45784534 I1=M(I0+1)+J0-1 J1=M(I0)+J0-24538 M1=M(I0)+J0-14540 IF(J(I0)=9) GOTO 45524542 IF(J(I0)=11) GOTO 45524544 Z9=Z0 CALL CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+14550 GOTO 45544552 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13) 4554 IF(J0D9) GOTO 45584556 GOTO 45784558 I1=M(I0)+J0-1 J1=M(I0+1)+J04562 M1=M(I0+1)+J0-14564 IF(J(I0)=9) GOTO 45764566 IF(J(I0)=11) GOTO 45764568 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+14574 GOTO 45784576 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4578 ENDDO4580 IF(J(I0)=1) GOTO 45864582 IF(J(I0)=9) GOTO 45864584 GOTO 45944586 I1=M(I0)+D9-2 J1=M(I0)+D9-14590 M1=M(I0+1)+D9-1 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4594 IF(J(I0)=3) GOTO 46004596 IF(J(I0)=11) GOTO 46004598 GOTO 46084600 I1=M(I0+1)+D9-1 J1=M(I0)+D9-14604 M1=M(I0+1)+D9 CALL BCE(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4608 IF(J(I0)/=5) GOTO 46144610 CALL EAB(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,M,I0,J,V0,V7,H10,H11,H12,H13)4612 CALL EDEC(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,M,I0,V7,H10,H11,H12,H13)4614 IF(J(I0)=6) GOTO 46204616 IF(J(I0)=7) GOTO 46204618 GOTO 46244620 CALL EAB (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,M,I0,J,V0,V7,H10,H11,H12,H13)4622 CALL EFGHI (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,M,I0,V7,H10,H11,H12,H13)4624 END SUBROUTINE EAB (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,M,I0,J,V0,V7,H10,H11,H12,H13) INTEGER B0,Z9,Z0,I1,J1,M1,M2,E3,F1,I0,V0,V7 DOUBLE PRECISION H9(E3) INTEGER F(B0),B(F1),W(B0),J(V0),M(V7) INTEGER H10(E3),H11(E3),H12(E3),H13(E3)4002 I1=M(I0+1)-64004 IF(J(I0)=5) J1=M(I0+1)-34008 IF(J(I0)/=5) J1=M(I0+1)-44012 IF(J(I0+1)7) GOTO 40204016 M1=M(I0+2)-5 GOTO 40224020 M1=M(I0+2)-24022 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4024 I1=M(I0+1)-6 J1=M(I0+1)-54028 IF(J(I0)=5) M1=M(I0+1)-34032 IF(J(I0)/=5) M1=M(I0+1)-44036 Z9=Z0 CALL CM(H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+14042 END SUBROUTINE EFGHI(B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,M,I0,V7,H10,H11,H12,H13) INTEGER B0,Z9,Z0,I1,J1,M1,M2,E3,F1,I0,V7 DOUBLE PRECISION H9(E3) INTEGER F(B0),B(F1),W(B0),M(V7) INTEGER H10(E3),H11(E3),H12(E3),H13(E3)4102 I1=M(I0+1)-3 J1=M(I0+1)-14106 M1=M(I0+1)-4 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4110 I1=M(I0+1)-3 J1=M(I0+1)-24114 M1=M(I0+1)-1 CALL BCE (B0,Z9,Z0,H9,E3,I1,J1,M1,M2,F,F1,B,W,H10,H11,H12,H13)4118 I1=M(I0+1)-5 J1=M(I0+1)-34122 M1=M(I0+1)-44124 Z9=Z0 CALL CM (H9,E3,Z9,I1,J1,M1,M2,H10,H11,H12,H13) Z0=Z0+14130 I1=M(I0+1)-3 J1=M(I0+1)-54134 M1=M(I0+1)-2 CALL BCE (B0,Z9,Z0,H9,