Главная Журналы RT(3)=AJP(I,J)*F(I,J+1,N) RT (4 ) =AJM(I, J) Td, J-1,N) RT(5)=-AP(I,J)*F(I,J,N) RT(6)=CON(I,J) RES = 0 TERM=1.0E-8 DO 30 IRT=1,6 RES = RES + RT CRT) 30 TERM=MAX(TERM,ABS (RTdRT) ) ) IF(ABS(RES/TERM).GT.CRIT(N))ICON=0 BLC=BLC+RES ENDIF 20 CONTINUE DEN0M=BL-PTX(I-1)*BLM IF(ABS(DENOM/BL).LT.SMALLl) DENOM=BIG PTX(I)=BLP/DENOM QTX(I)=(BLC+BLM*QTX(I-1))/DENOM 10 CONTINUE IF(NTT.NE.1.AND.ICON.EQ.1) GO TO 990 IF(KBLOC(NF).EQ.O) GO TO 80 BL=0. DO 40 I=L2,2,-1 BL=BL*PTX(I)+QTX(I) DO 4 0 J=2,M2 IF(AP(I,J).LT.BIGl) F(I,J,N)=F(I,J,N)+BL 40 CONTINUE COME HERE TO PERFORM THE J-DIRECTION BLOCK CORRECTION С--- PTY(1)=0. QTY(1)=0. DO 50 J=2,M2 BL=SMALL BLP=0. BLM=0. BLC=0. DO 60 1=2,L2 IF(AP(I,J).LT.BIGl) THEN BL=BL+AP(I,J) IF(AP(I+1,J).LT.BIGl) BL=BL-AIP(I,J) IF(AP(I-1,J).LT.BIGl) BL=BL-AIM(I,J) IF(AP(I,J+I).LT.BIGl) BLP=BLP+AJP(I,J) IF(AP(I,J-1).LT.BIGl) BLM=BLM+AJM(I,J) BLC=BLC+CON(I,J)+AIP(I,J) *F(I + l,J,N)+AIM(I, J)*F(I-1,J,N) 1 +AJP(I,J)*F(I,J+I,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N) ENDIF 60 CONTINUE DENOM=BL-PTY(J-1)*BLM IF(ABS(DENOM/BL).LT.SMALLl) DENOM=BIG РТУ ; J) =BLP/DEK0[-1 -JTY { J) = (BLC + BLMQTY t J-1 ) ) /DF,NOM 5 0 CONTINUE BLO . DO 7 0 J=M2,2,-1 BL = BL*PTY ( J) + ОТГ ( J) DO 70 1=2,L2 IF (AP ( I, J) . LT .BIGl ) F (I, J, N) = F(I, J, Ii) +BL 7 0 CONTINUE 8 0 CONTINUE CARRY OUT THE I-DIRECTION TDMA DO 9 0 JJ2,MM J=MIN(JJ,MM2-JJ) PTX a)= 0. OTX(1)-0 DO 100 1=2,L2 DENOM=AP(I,J)-PTX(I-l)*AIM(I,J) PTX {I ) -AIP (I, Ji /DENOM TEMP = CON {1, J) tAJP ( I, J ) -*F (I, J+1, !!) +AJM (I, J) * F (I, J-1, tl) QTXI I) = (TEMP + AIM(I,J)"OTX{I-l) ) /DEUOM 10 0 CONTINUE DO 110 I=L2,2,-1 110 F{ I, J, N) = F{I + 1, J, N) "PTX (1) +OTX (I ) 90 CONTINUE CARRY OUT THE J-DIRECTION TDMA r--------- DO 120 11=2,LL I=MIN(II,LL2-1I) PTY(1)=0. QTY(1)=0 DO 130 J=2,M2 DENOM=AP(I,J)-PTY(J-1)*AJM{I, J) PTY(J)=AJP(I,J)/DENOM TEMP=CON (I,J)+AIP(I,J)*F{I + l,J,tl)+AIM(I,J)*F{l-l,J,N) QTY{J)=(TEMP+AJM(I,J)*QTY(J-1))/DENOM 130 CONTINUE DO 14 0 J=M2,2,-i 110 F(I,J,N)=F(I,J+1,N)*PTY{J)+OTY(J) 120 CONTINUE С------- 999 CONTINUE NTC(N)=NTT GO TO 991 990 NTC{N)=NTT-1 991 CONTINUE CALCULATE THE UNKNOWN BOUNDARY VALUES AliD FLUXES С-------- DO 160 1=2,L2 TEMP=AJM(I,1)*(F(I,3,N)-F(I,2,N)) IF(KBCJ1(I).EQ.2) 1 F(I, 1,N) = (AJPd, 1) *F(I, 2,N) -TEMP+CONd, 1) ) /АР(1, 1) FLUXJl (I, N) =AJP (I, l)*(Fd,l,N)-F(I,2,N)) +TEMP TEMP=AJP(I,M1)*(F(I,M3,N)-F(I,M2,N)) IF(KBCMl(I).EQ.2) 1 F(I,M1,N)=(AJM(I,M1)*F(I,M2,N)-TEMP+CON(I,M1))/АР(I,M1) 160 FLUXMl(I,N)=AJM(I,M1)*(F(I,M1,N)-F(I,M2,N))+TEMP DO 17 0 J=2,M2 TEMP=AIM(1,J)*(F(3, J,N)-F(2, J,N) ) IF(KBCI1(J).EQ.2) 1 F(1,J,N)=(AIP(1,J)*F(2,J,N)-TEMP+CON(1,J))/АР(1,J) FLUXIl(J,N)=AIP(1,J)*(F(1,J,N)-F(2,J,N))+TEMP TEMP=AIP(L1,J)*(F(L3,J,N)-F(L2,J,N)) IF(KBCL1(J).EQ.2) 1 F(L1,J,N)=(AIM(L1,J)*F(L2,J,N)-TEMP+CON(L1,J))/AP(L1,J) 17 0 FLUXLl(J,N)=AIM(LI,J)*(F(Ll,J,N)-F(L2,J,N))+TEMP С COME HERE TO RESET CON,AP,KBC,FLXC, AND FLXP С-- DO 18 0 J=2,M2 KBCIl(J)=l KBCLl(J)=l FLXCIl(J)=0. FLXCLl(J)=0. FLXPll(J)=0. FLXPLl(J)=0. DO 180 1=2,L2 CON(I,J)=0. AP (I, J) =0. 180 CONTINUE DO 190 1=2,L2 KBCJl(I)=l KBCMl(I)=l FLXCJl(I)=0. FLXCMl(I)=0. FLXPJl(I)=0. FLXPMl(I)=0. 190 CONTINUE RETURN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TOOLS $INCLUDE:COMMON 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 [ 94 ] 95 96 97 98 99 |