📄 coneq.f
字号:
COE3(1) = (ODT(2,1)+CONVO(1,1))/DDX(1)
DO 648 I = 2,NBL-1
COE3(I) = (DDY(I)*ODT(I+1,1)-DDY(I-1)*ODT(I,1)+
* CONVO(I,1)*DDY(I)-DDY(I-1)*CONVO(I-1,1))
* /(DDX(I)*DDY(I))
648 CONTINUE
COE3(NBL) = (-DDY(NBL-1)*ODT(NBL,1)+DDY(NBL)*CONVO
* (NBL,1)-DDY(NBL-1)*CONVO(NBL-1,1))/(DDY(NBL)
* *DDX(NBL))
ELSE
COE3(1) = 2.*RRP(1)*(ODT(2,1)+CONVO(1,1))/DDX(1)
DO 649 I = 2,NBL-1
COE3(I) = 2.*(RRP(I)*(ODT(I+1,1)+CONVO(I,1))-
* RRP(I-1)*(ODT(I,1)+CONVO(I-1,1)))/DDX(I)
649 CONTINUE
COE3(NBL) = 2.*(RR(NBL)*CONVO(NBL,1)-RRP(NBL-1)*
* (ODT(NBL,1)+CONVO(NBL-1,1)))/DDX(NBL)
ENDIF
IF (NY.GT.1) THEN
DO 650 I = 1,NX
COE4(I) = (ODT(I+NX,2)+CONVO(I,2))/DDY(I)
650 CONTINUE
DO 651 I = NX+1,NBL-NX
COE4(I) = (ODT(I+NX,2)-ODT(I,2)+CONVO(I,2)-
cbug 9.0
* CONVO(I-NX,2))/DDY(I)
651 CONTINUE
DO 652 I = NBL-NX+1,NBL
COE4(I) = (-ODT(I,2)+CONVO(I,2)-CONVO(I-NX,2))
* /DDY(I)
652 CONTINUE
ENDIF
IF (NZ.GT.1) THEN
DO 653 I = 1,NXNY
COE5(I) = (ODT(I+NXNY,3)+CONVO(I,3))/DDZ(I)
653 CONTINUE
IF (NZ.EQ.2) THEN
DO 654 I = NXNY+1,NBL
COE5(I) = (-ODT(I,3)+CONVO(I,3)-CONVO(I-NXNY
* ,3))/DDZ(I)
654 CONTINUE
ELSE
DO 655 I = NXNY+1,NBLM3
COE5(I) = (ODT(I+NXNY,3)-ODT(I,3)+CONVO(I,3)-
* CONVO(I-NXNY,3))/DDZ(I)
655 CONTINUE
DO 656 I = NBLM3+1,NBL
COE5(I) = (-ODT(I,3)+CONVO(I,3)-CONVO(I-NXNY
* ,3))/DDZ(I)
656 CONTINUE
ENDIF
ENDIF
DO 657 I = 1,NBL
DCO(I) = -DT*(COE3(I)+COE4(I)+COE5(I))/POR(I)
657 CONTINUE
ENDIF
ENDIF
C
C ----FLEXI-GRID
C SCALE DISPERSIVE FLUX BY CELL AREAS RE-CALCULATED IN "SOLMAT"
C
IF (ICOORD.EQ.4) THEN
DO 658 I = 1,NBL-1
COEX(I+1) = COEX(I+1)*DUM5(I)
658 CONTINUE
DO 659 I = 1,NBL-NX
COEYU(I+NX) = COEYU(I+NX)*DUM6(I)
659 CONTINUE
DO 660 I = 1,NBL-NXNY
COEZT(I+NXNY) = COEZT(I+NXNY)*DUM7(I)
660 CONTINUE
C
C SCALE CONVECTIVE FLUX BY CELL AREAS RE-CALCULATED IN "SOLMAT"
C
DO 661 I = 1,NBL
CONVX(I,KC,1) = CONVX(I,KC,1)*DUM5(I)
CONVY(I,KC,1) = CONVY(I,KC,1)*DUM6(I)
CONVZ(I,KC,1) = CONVZ(I,KC,1)*DUM7(I)
661 CONTINUE
C
C CALC' TRUE SURFACE FLUX DIFFERENCES
C WHERE COEX = DISPERSION : CONVX = CONVECTION
C
COE3(1) = (COEX(2)+CONVX(1,KC,1))
DO 662 I = 2,NBL-1
COE3(I) = COEX(I+1)-COEX(I)+
* CONVX(I,KC,1)-CONVX(I-1,KC,1)
662 CONTINUE
COE3(NBL) = -COEX(NBL)+CONVX(NBL,KC,1)
* -CONVX(NBL-1,KC,1)
IF (NY.GT.1) THEN
DO 663 I = 1,NX
COE4(I) = COEYU(I+NX)+CONVY(I,KC,1)
663 CONTINUE
DO 664 I = NX+1,NBL-NX
COE4(I) = COEYU(I+NX)-COEYU(I)+CONVY(I,KC,1)-
* CONVY(I-NX,KC,1)
664 CONTINUE
DO 665 I = NBL-NX+1,NBL
COE4(I) = -COEYU(I)+CONVY(I,KC,1)-CONVY(I-NX,KC,1)
665 CONTINUE
ENDIF
IF (NZ.GT.1) THEN
DO 666 I = 1,NXNY
COE5(I) = COEZT(I+NXNY)+CONVZ(I,KC,1)
666 CONTINUE
IF (NZ.EQ.2) THEN
DO 667 I = NXNY+1,NBL
COE5(I) = -COEZT(I)+CONVZ(I,KC,1)-CONVZ(I-NXNY
* ,KC,1)
667 CONTINUE
ELSE
DO 668 I = NXNY+1,NBLM3
COE5(I) = COEZT(I+NXNY)-COEZT(I)+CONVZ(I,KC,1)-
* CONVZ(I-NXNY,KC,1)
668 CONTINUE
DO 669 I = NBLM3+1,NBL
COE5(I) = -COEZT(I)+CONVZ(I,KC,1)-CONVZ(I-NXNY
* ,KC,1)
669 CONTINUE
ENDIF
ENDIF
ELSE
IF (ICOORD.NE.2) THEN
COE3(1) = (COEX(2)+CONVX(1,KC,1))/DDX(1)
DO 670 I = 2,NBL-1
COE3(I) = (DDY(I)*COEX(I+1)-DDY(I-1)*COEX(I)+
* CONVX(I,KC,1)*DDY(I)-DDY(I-1)*CONVX(I-1,KC,1))
* /(DDX(I)*DDY(I))
670 CONTINUE
C-------------------------------------------------------
IF(NBL==1)THEN
COE3(NBL) = (-0.0*COEX(NBL)+DDY(NBL)*CONVX(NBL
* ,KC,1)-DDY(NBL-1)*0.0)/(DDY(NBL)
* *DDX(NBL))
ELSE
C-------------------------------------------------------
COE3(NBL) = (-DDY(NBL-1)*COEX(NBL)+DDY(NBL)*CONVX(NBL
* ,KC,1)-DDY(NBL-1)*CONVX(NBL-1,KC,1))/(DDY(NBL)
* *DDX(NBL))
C-------------------------------------------------------
END IF
C-------------------------------------------------------
ELSE
COE3(1) = 2.*RRP(1)*(COEX(2)+CONVX(1,KC,1))/DDX(1)
DO 671 I = 2,NBL-1
COE3(I) = 2.*(RRP(I)*(COEX(I+1)+CONVX(I,KC,1))-
* RRP(I-1)*(COEX(I)+CONVX(I-1,KC,1)))/DDX(I)
671 CONTINUE
COE3(NBL) = 2.*(RR(NBL)*CONVX(NBL,KC,1)-
* RRP(NBL-1)*(COEX(NBL)+CONVX(NBL-1,KC,1)))/DDX(NBL)
ENDIF
IF (NY.GT.1) THEN
DO 672 I = 1,NX
COE4(I) = (COEYU(I+NX)+CONVY(I,KC,1))/DDY(I)
672 CONTINUE
DO 680 I = NX+1,NBL-NX
COE4(I) = (COEYU(I+NX)-COEYU(I)+CONVY(I,KC,1)-CONVY
* (I-NX,KC,1))/DDY(I)
680 CONTINUE
DO 690 I = NBL-NX+1,NBL
COE4(I) = (-COEYU(I)+CONVY(I,KC,1)-CONVY(I-NX,KC,1))
* /DDY(I)
690 CONTINUE
ENDIF
IF (NZ.GT.1) THEN
DO 700 I = 1,NXNY
COE5(I) = (COEZT(I+NXNY)+CONVZ(I,KC,1))/DDZ(I)
700 CONTINUE
IF (NZ.EQ.2) THEN
DO 710 I = NXNY+1,NBL
COE5(I) = (-COEZT(I)+CONVZ(I,KC,1)-CONVZ(I-NXNY
* ,KC,1))/DDZ(I)
710 CONTINUE
ELSE
DO 720 I = NXNY+1,NBLM3
COE5(I) = (COEZT(I+NXNY)-COEZT(I)+CONVZ(I,KC,1)-
* CONVZ(I-NXNY,KC,1))/DDZ(I)
720 CONTINUE
DO 730 I = NBLM3+1,NBL
COE5(I) = (-COEZT(I)+CONVZ(I,KC,1)-CONVZ(I-NXNY
* ,KC,1))/DDZ(I)
730 CONTINUE
ENDIF
ENDIF
ENDIF
C
C ADSORPTION OF TRACERS (RETARDATION)
C
IF (KC.GT.MTWM1.AND.KC.LT.MTWM1+1+NTW) THEN
IT = KC-MTWM1
IF (ABS(RET(IT)).GE.PRCSN) THEN
DO 740 I = 1,NBL
IF (SF(I,1).GT.ZERO) THEN
CON = ONE/(1.+RET(IT)/SF(I,1))
COE3(I) = COE3(I)*CON
COE4(I) = COE4(I)*CON
COE5(I) = COE5(I)*CON
ENDIF
740 CONTINUE
ENDIF
ENDIF
IF (KC.GT.MTAM1.AND.KC.LT.MTAM1+1+NTA) THEN
IT = KC-MTAM1
IF (ABS(RET(IT)).GE.PRCSN) THEN
DO 742 I = 1,NBL
IF (S(I,NPHAS).GT.ZERO) THEN
CON = ONE/(1.+RET(IT)/S(I,NPHAS))
COE3(I) = COE3(I)*CON
COE4(I) = COE4(I)*CON
COE5(I) = COE5(I)*CON
ENDIF
742 CONTINUE
ENDIF
ENDIF
C
IF (ICOORD.EQ.4) THEN
C
C ----FLEXI-GRID
C DIVIDE FLUX INTEGRAL BY TRUE VOLUME
C
DO 745 I = 1,NBL
DC(I) = (-DT/EPHI1)*(COE3(I)+COE4(I)+COE5(I))/VB(I)
745 CONTINUE
ELSE
DO 750 I = 1,NBL
DC(I) = (-DT/EPHI1)*(COE3(I)+COE4(I)+COE5(I))/POR(I)
750 CONTINUE
ENDIF
C
C ADD SINK AND SOURCE TERM (ASSUME THE CORRECTION DUE TO
C FLUID COMPRESSIBILITY IS SMALL AT THIS TIME (9/19/89)
C FOR COMPONENTS WITH MEQ/ML UNITS, USE M^3 INSTEAD OF FT3
C (CONV. FROM FT3 TO M3 = 0.028316847)
C
WCONS = 24.
IF (KC.EQ.5.OR.KC.EQ.6.OR.(IREACT.EQ.1.AND.KC.EQ.MGLM1+NG).OR.
* (LGC.AND.KC.GE.MGCM1+1.AND.KC.LE.MGCM1+NGC)) THEN
CONFAC = 0.028316847
ELSE
CONFAC = 1.0
ENDIF
DO 760 M = 1,NWELL
ID = IDW(M)
DO 760 IWB = 1,NWBC(ID)
IJK2 = IJKPOS(IWC(ID,IWB),JWC(ID,IWB),KWC(ID,IWB))
DO 760 L = 1,NPHAS
IF (Q(ID,IWB,L).GT.ZERO) THEN
CW = C(NBL+ID,KC,L)
CUMI(KC) = CUMI(KC)+Q(ID,IWB,L)*CONFAC*CW*DT
ELSE
CW = C(IJK2,KC,L)
C
C CONCENTRATION WELL MODEL ONLY FOR VERTICAL WELL
C
IF (ICWM.EQ.0.OR.IDIR(ID).NE.3) GO TO 755
IF (NY.NE.1) THEN
TERM0 = 28.*CW
TERM1 = CW
IF (JW(ID).NE.1) TERM1 = C(IJK2-NX,KC,L)
TERM2 = CW
IF (IW(ID).NE.1) TERM2 = C(IJK2-1,KC,L)
TERM3 = CW
IF (IW(ID).NE.NX) TERM3 = C(IJK2+1,KC,L)
TERM4 = CW
IF (JW(ID).NE.NY) TERM4 = C(IJK2+NX,KC,L)
CW = (TERM0-TERM1-TERM2-TERM3-TERM4)/WCONS
ENDIF
755 CUMP(KC) = CUMP(KC)+Q(ID,IWB,L)*CONFAC*CW*DT
ENDIF
DC(IJK2) = DC(IJK2)+(DT*CW/EPHI1)*Q(ID,IWB,L)
* /VB(IJK2)
IF (IMASS.EQ.2.AND.(KC.EQ.2.OR.(KC.GT.MMOM1.AND.
* KC.LE.MMOM1+NO)).AND.L.EQ.1) THEN
DCO(IJK2) = DCO(IJK2)+DT*CW*Q(ID,IWB,L)/VB(IJK2)
ENDIF
760 CONTINUE
C
C EXTERNAL BOUNDARY FOR RADIAL COORDINATE
C
IF (ICOORD.EQ.2) THEN
M = NWELL+1
IDW(M) = NWELL+1
ID = IDW(M)
DO 765 K = 1,NZ
IJK2 = K*NX
DO 765 L = 1,NPHAS
IF (Q(ID,K,L).GE.ZERO) THEN
QE = Q(ID,K,L)*CE(IJK2,KC,L)*DT*CONFAC
CUMI(KC) = CUMI(KC)+QE
ELSE
QE = Q(ID,K,L)*C(IJK2,KC,L)*DT*CONFAC
CUMP(KC) = CUMP(KC)+QE
ENDIF
DC(IJK2) = DC(IJK2)+QE/(EPHI1*CONFAC*VB(IJK2))
IF (IMASS.EQ.2.AND.(KC.EQ.2.OR.(KC.GT.MMOM1.AND.
* KC.LE.MMOM1+NO)).AND.L.EQ.1) THEN
DCO(IJK2) = DCO(IJK2)+QE/VB(IJK2)
ENDIF
765 CONTINUE
ENDIF
C
C BOUNDARIES
C
IF (ICOORD.NE.2.AND.IBOUND.EQ.1) THEN
IF (IBL.EQ.1) THEN
DO 766 J = 1,NY
DO 766 K = 1,NZ
IL = IJKPOS(1,J,K)
DO 766 L = 1,NPHAS
IF (QL(IL,L).GE.ZERO) THEN
QE = QL(IL,L)*CE(IL,KC,L)*DT*CONFAC
CUMI(KC) = CUMI(KC)+QE
ELSE
QE = QL(IL,L)*C(IL,KC,L)*DT*CONFAC
CUMP(KC) = CUMP(KC)+QE
C IF (KC.EQ.2.AND.C(IL,KC,L).GT.EPAC) THEN
C
C IF (MOD(ICNT,250).EQ.0)
C * WRITE (8,*) 'OIL IS LEAVING THE LEFT ',
C * 'BOUNDARY AT TIME=',T,' AND BLOCK=',IL
C IFLAGL = 1
C ENDIF
ENDIF
DC(IL) = DC(IL)+QE/(EPHI1*CONFAC*VB(IL))
IF (IMASS.EQ.2.AND.(KC.EQ.2.OR.(KC.GT.MMOM1.AND.
* KC.LE.MMOM1+NO)).AND.L.EQ.1) THEN
DCO(IL) = DCO(IL)+QE/VB(IL)
ENDIF
766 CONTINUE
ENDIF
IF (IBR.EQ.1) THEN
DO 767 J = 1,NY
DO 767 K = 1,NZ
IR = IJKPOS(NX,J,K)
DO 767 L = 1,NPHAS
IF (QL(IR,L).GE.ZERO) THEN
QE = QL(IR,L)*CE(IR,KC,L)*DT*CONFAC
CUMI(KC) = CUMI(KC)+QE
ELSE
QE = QL(IR,L)*C(IR,KC,L)*DT*CONFAC
CUMP(KC) = CUMP(KC)+QE
C IF (KC.EQ.2.AND.C(IR,KC,L).GT.EPAC) THEN
C IF (M⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -