计算传热学程序.doc

计算传热学课程报告一、问题概述：有限单元法是上个世纪五、六十年代首先在力学中发展起来的数值计算方法，由于它是基于变分原理，理论基础统一，对于复杂边界的适应程度比较好，所以很快的在其它领域得到运用，其中就包括了在传热学中的运用。本次计算传热学的课程就是对有限单元法在传热学中运用的一个学习与练习。有限单元法处理问题的步骤，首先是建立有限元模型也即是将问题离散化，它的主要步骤之一就是将要计算的物体进行有限元的划分；第二步，进行单元分析也就是将变分原理运用到问题的方程与单元中，形成单元刚度矩阵；第三步，进行整体刚度矩阵的组集；最后就是引入边界条件进行求解的过程。在计算传热学的课程中，主要完成了两个任务：第一，是将一个比较复杂的活塞进行了网格划分，并编译成一个通用性比较好的程序。第二，在前一个程序的基础上，加入计算过程，运用焓法，对一个比较简单的平面相变问题进行了计算。二、划分单元网格：划分单元网格是将问题进行有限元法分析的基础，但是如果在图纸上进行手工的单元划分，不但繁琐、容易出错，而且也不利于进一步计算程序的利用。因此有必要编辑一个程序，以自动完成划分网格的目的。网格的自动划分必须遵循以下的几条规则：(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,