計(jì)算傳熱學(xué)程序.doc

計(jì)算傳熱學(xué)課程報(bào)告一、問題概述：有限單元法是上個(gè)世紀(jì)五、六十年代首先在力學(xué)中發(fā)展起來的數(shù)值計(jì)算方法，由于它是基于變分原理，理論基礎(chǔ)統(tǒng)一，對(duì)于復(fù)雜邊界的適應(yīng)程度比較好，所以很快的在其它領(lǐng)域得到運(yùn)用，其中就包括了在傳熱學(xué)中的運(yùn)用。本次計(jì)算傳熱學(xué)的課程就是對(duì)有限單元法在傳熱學(xué)中運(yùn)用的一個(gè)學(xué)習(xí)與練習(xí)。有限單元法處理問題的步驟，首先是建立有限元模型也即是將問題離散化，它的主要步驟之一就是將要計(jì)算的物體進(jìn)行有限元的劃分；第二步，進(jìn)行單元分析也就是將變分原理運(yùn)用到問題的方程與單元中，形成單元?jiǎng)偠染仃?；第三步，進(jìn)行整體剛度矩陣的組集；最后就是引入邊界條件進(jìn)行求解的過程。在計(jì)算傳熱學(xué)的課程中，主要完成了兩個(gè)任務(wù)：第一，是將一個(gè)比較復(fù)雜的活塞進(jìn)行了網(wǎng)格劃分，并編譯成一個(gè)通用性比較好的程序。第二，在前一個(gè)程序的基礎(chǔ)上，加入計(jì)算過程，運(yùn)用焓法，對(duì)一個(gè)比較簡(jiǎn)單的平面相變問題進(jìn)行了計(jì)算。二、劃分單元網(wǎng)格：劃分單元網(wǎng)格是將問題進(jìn)行有限元法分析的基礎(chǔ)，但是如果在圖紙上進(jìn)行手工的單元?jiǎng)澐?，不但繁瑣、容易出錯(cuò)，而且也不利于進(jìn)一步計(jì)算程序的利用。因此有必要編輯一個(gè)程序，以自動(dòng)完成劃分網(wǎng)格的目的。網(wǎng)格的自動(dòng)劃分必須遵循以下的幾條規(guī)則：(1).要嚴(yán)格區(qū)分邊界單元與內(nèi)部單元，并且嚴(yán)格區(qū)分邊界單元不同的組；(2).單元標(biāo)號(hào)必須先標(biāo)志內(nèi)部單元，然后依次標(biāo)志第一類邊界條件，第二類邊界條件，第三類邊界條件，如果同一類邊界條件中有不同的組，那么也必須嚴(yán)格先劃分第一組，然后第二組，第三組；(3). 對(duì)于邊界單元，每一個(gè)邊界單元必須只有一條邊在邊界上，而且為了程序的簡(jiǎn)單，一般是j,m邊作為邊界；(4).內(nèi)部單元節(jié)點(diǎn)標(biāo)號(hào)必須遵循逆時(shí)針方向的規(guī)則；(5). 一個(gè)單元中只能有一種材料組成。遵循以上的規(guī)則，用FORTRAN 90編制了一個(gè)對(duì)形狀比較復(fù)雜的活塞的網(wǎng)格劃分，由于在編制過程中考慮了多種情況，所以這個(gè)程序有比較好的通用性，只需要輸入不同的數(shù)據(jù)，程序也可以對(duì)許多其它情況進(jìn)行劃分。需要指出的是，由于FORTRAN 90程序?qū)τ谥茍D功能比較弱，所以下面的圖是用VB 6.0的程序做出的，由于該網(wǎng)格劃分程序集成了后續(xù)對(duì)第一類邊界條件和第三類邊界條件的焓法計(jì)算程序，故該程序源代碼將在最后統(tǒng)一給出。網(wǎng)格劃分的結(jié)果如圖(1)。需要輸入的初始數(shù)據(jù)主要有：邊界單元分組總數(shù)、邊界單元分組中前一組的最后一個(gè)單元號(hào)、各組邊界單元節(jié)點(diǎn)數(shù)、各邊界單元邊界節(jié)點(diǎn)號(hào)、每一條層線的左右端點(diǎn)、坐標(biāo)，每一層單元?jiǎng)澐炙鶎俚念愋汀D（）三、焓法有限單元法原理：帶有相變的傳熱問題，又被稱為斯蒂芬問題，在冶金、鑄造、建筑、冷凍、航天和醫(yī)療等領(lǐng)域有著廣泛的運(yùn)用。由于是這類問題存在的相變過程，使得求解區(qū)域中存在著一個(gè)隨著時(shí)間移動(dòng)的固液或者固氣界面，這一界面使得這種問題的求解非常困難。一般說來，處理這種問題有兩種不同的思路，一種思路是，首先著眼于相變界面的求解，確定相變界面以后再分別處理固相或者液相的溫度分布；另一種思路是將該問題看作是“單相”的非線性導(dǎo)熱問題，首先確定整個(gè)求解區(qū)域上的溫度分布或者焓的分布，然后把達(dá)到相變溫度的位置定為相變界面。在實(shí)際運(yùn)用中，后一種思路比較簡(jiǎn)單和實(shí)用而得到了廣泛的使用；目前后一種處理方式中主要有焓法和顯熱容法，這里主要討論焓法。焓法有限元法的主要思路是，不把溫度作為求解的變量，而是把焓作為求解變量，因此可以在固相和液相的整個(gè)區(qū)域建立統(tǒng)一的方程進(jìn)行求解。然后根據(jù)焓與溫度的關(guān)系確定整個(gè)區(qū)域的溫度與相變的界面。焓法的優(yōu)點(diǎn)就是不需要跟蹤相變的界面而可以對(duì)整個(gè)區(qū)域進(jìn)行統(tǒng)一的求解。焓法的主要推導(dǎo)如下：在一般過程中，傳熱微分方程由(1)式表示：（）其中為導(dǎo)熱系數(shù)，為物體的密度， 為物體的定壓熱容，為物體的溫度，而表示時(shí)間。由于焓與溫度存在如下的關(guān)系：所以（）式可以化為（）式：（）該式子即是在整個(gè)區(qū)域上都適用的焓法的基本方程，對(duì)于這個(gè)方程，在空間上用有限單元法離散，在時(shí)間上用向后差分格式離散，就可以得到如下的(3)式：（）其中：整體剛度矩陣 變焓矩陣 時(shí)間步長(zhǎng)t時(shí)刻熱焓時(shí)刻熱焓t時(shí)刻右端列向量值得注意的是，首先使用變分法的使用要直接對(duì)焓而不是溫度變分；其次整體剛度矩陣和變焓矩陣是與時(shí)間無關(guān)的的矩陣，而右端向量則與邊界條件有關(guān)。最后焓與溫度的關(guān)系由（）式給出：（）四、焓法實(shí)例：為了簡(jiǎn)單起見，運(yùn)用上述的焓法對(duì)如圖(2)的簡(jiǎn)單區(qū)域的相變問題進(jìn)行求解，其條件如下：邊長(zhǎng)為80cm60cm，周圍用溫度為500K，鐵水溫度為1833K，比熱為711.62J/kg.k，潛熱為271100J/kg，密度為7800kg/m3，導(dǎo)熱系數(shù)為33.5W/m. 相溫度為1790K，使用第三類邊界條件，換熱系數(shù)1950W/Km2。圖 （2）圖中1-4，1-15等邊都是第三類邊界，其余單元?jiǎng)t作為邊界單元處理。程序源代碼如下： 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,