📄 3dmain.for
字号:
$DEBUG
*----------------------------MAIN PROGRAM-----------------------------------
****************************************************************************
LOGICAL LSTOP
COMMON/CNTL/LSTOP
****************************************************************************
* CALL SETUP0
OPEN(08,FILE='RESULT')
CALL GRID
CALL SETUP1
CALL START
10 CALL DENSE
CALL BOUND
CALL OUTPUT
IF(.NOT.LSTOP) GO TO 15
C* CLOSE(08)
STOP
15 CALL SETUP2
GO TO 10
END
*---------------------------------------------------------------------------
SUBROUTINE DIFLOW
****************************************************************************
COMMON/COEF/FLOW,DIFF,ACOF
****************************************************************************
ACOF=DIFF
IF(FLOW .LT. 1.E-35)RETURN
TEMP=DIFF-ABS(FLOW)*0.1
ACOF=0.
IF(TEMP .LE. 0. ) RETURN
TEMP=TEMP/DIFF
ACOF=DIFF*TEMP**5
RETURN
END
*--------------------------------------------------------------------------
SUBROUTINE SOLVE
****************************************************************************
PARAMETER (ID=25,JD=55,KD=25)
DOUBLE PRECISION TITLE
LOGICAL LSOLVE,LPRINT,LBLK,LSTOP
COMMON F(ID,JD,KD,10),P(ID,JD,KD),RHO(ID,JD,KD),GAM(ID,JD,KD),
&CON(ID,JD,KD),AP(ID,JD,KD),AKP(ID,JD,KD),AKM(ID,JD,KD),
&AIP(ID,JD,KD),AIM(ID,JD,KD),AJP(ID,JD,KD),AJM(ID,JD,KD),
&X(ID),XU(ID),XDIF(ID),XCV(ID),XCVS(ID),XCVI(ID),XCVIP(ID),
&Y(JD),YV(JD),YDIF(JD),YCV(JD),YCVS(JD),YCVJ(JD),YCVJP(JD),
&Z(KD),ZW(KD),ZDIF(KD),ZCV(KD),ZCVS(KD),ZCVK(KD),ZCVKP(KD),
&AX(JD,KD),AY(ID,KD),AZ(ID,JD)
COMMON DU(ID,JD,KD),DV(ID,JD,KD),DW(ID,JD,KD),
& FX(ID),FXM(ID),FY(JD),FYM(JD),FZ(KD),FZM(KD),PT(JD),QT(JD)
COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,N1,N2,N3,
& IST,JST,KST,ITER,LAST,TITLE(13),RELAX(13),TIME,DT,XL,YL,ZL,
& IPREF,JPREF,KPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),
& RHOCON
*****************************************************************************
ISTF=IST-1
JSTF=JST-1
KSTF=KST-1
IT1=L2+IST
JT1=M2+JST
KT1=N2+KST
****************************************************************************
DO 999 NT=1,NTIMES(NF)
DO 999 N=NF,NF
*-----------------------------------------------------------------------
DO 90 K=KST,N2
DO 90 J=JST,M2
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J,K,N)
DO 70 I=IST,L2
DENOM=AP(I,J,K)-PT(I-1)*AIM(I,J,K)
PT(I)=AIP(I,J,K)/DENOM
TEMP=CON(I,J,K)+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
&+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
QT(I)=(TEMP+AIM(I,J,K)*QT(I-1))/DENOM
70 CONTINUE
DO 80 II=IST,L2
I=IT1-II
80 F(I,J,K,N)=F(I+1,J,K,N)*PT(I)+QT(I)
90 CONTINUE
*-----------------------------------------------------------------------
DO 190 KK=KST,N2
K=KT1-KK
DO 190 JJ=JST,M2
J=JT1-JJ
IF(J.EQ.M2.AND.K.EQ.N2) GOTO 190
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J,K,N)
DO 170 I=IST,L2
DENOM=AP(I,J,K)-PT(I-1)*AIM(I,J,K)
PT(I)=AIP(I,J,K)/DENOM
TEMP=CON(I,J,K)+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
&+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
QT(I)=(TEMP+AIM(I,J,K)*QT(I-1))/DENOM
170 CONTINUE
DO 180 II=IST,L2
I=IT1-II
180 F(I,J,K,N)=F(I+1,J,K,N)*PT(I)+QT(I)
190 CONTINUE
*-----------------------------------------------------------------------
DO 290 K=KST,N2
DO 290 I=IST,L2
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF,K,N)
DO 270 J=JST,M2
DENOM=AP(I,J,K)-PT(J-1)*AJM(I,J,K)
PT(J)=AJP(I,J,K)/DENOM
TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
&+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
QT(J)=(TEMP+AJM(I,J,K)*QT(J-1))/DENOM
270 CONTINUE
DO 280 JJ=JST,M2
J=JT1-JJ
280 F(I,J,K,N)=F(I,J+1,K,N)*PT(J)+QT(J)
290 CONTINUE
*-----------------------------------------------------------------------
DO 390 KK=KST,N2
K=KT1-KK
DO 390 II=IST,L2
I=IT1-II
IF(I.EQ.L2.AND.K.EQ.N2) GOTO 390
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF,K,N)
DO 370 J=JST,M2
DENOM=AP(I,J,K)-PT(J-1)*AJM(I,J,K)
PT(J)=AJP(I,J,K)/DENOM
TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
&+AKP(I,J,K)*F(I,J,K+1,N)+AKM(I,J,K)*F(I,J,K-1,N)
QT(J)=(TEMP+AJM(I,J,K)*QT(J-1))/DENOM
370 CONTINUE
DO 380 JJ=JST,M2
J=JT1-JJ
380 F(I,J,K,N)=F(I,J+1,K,N)*PT(J)+QT(J)
390 CONTINUE
*-----------------------------------------------------------------------
DO 490 J=JST,M2
DO 490 I=IST,L2
PT(KSTF)=0.
QT(KSTF)=F(I,J,KSTF,N)
DO 470 K=KST,N2
DENOM=AP(I,J,K)-PT(K-1)*AKM(I,J,K)
IF(DENOM.EQ.0.) THEN
WRITE(*,998) DENOM,PT(K-1),AKM(I,J,K),AP(I,J,K),NF,I,J,K
WRITE(2,998) DENOM,PT(K-1),AKM(I,J,K),AP(I,J,K),NF,I,J,K
998 FORMAT(1X,'DENOM=',F15.5,4X,'PT(K-1)=',
& F15.5,4X,'AKM(I,J,K)=',F15.5,4X,'AP(I,J,K)=',F15.5,4X,'NF=',I3,
& 'I=',I3,2X,'J=',I3,2X,'K=',I3)
ENDIF
PT(K)=AKP(I,J,K)/DENOM
IF(DENOM.EQ.0.) WRITE(*,*) 'OK'
TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
&+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
QT(K)=(TEMP+AKM(I,J,K)*QT(K-1))/DENOM
470 CONTINUE
DO 480 KK=KST,N2
K=KT1-KK
480 F(I,J,K,N)=F(I,J,K+1,N)*PT(K)+QT(K)
490 CONTINUE
*-----------------------------------------------------------------------
DO 590 JJ=JST,M2
J=JT1-JJ
DO 590 II=IST,L2
I=IT1-II
IF(I.EQ.L2.AND.J.EQ.M2) GOTO 590
PT(KSTF)=0.
QT(KSTF)=F(I,J,KSTF,N)
DO 570 K=KST,N2
DENOM=AP(I,J,K)-PT(K-1)*AKM(I,J,K)
PT(K)=AKP(I,J,K)/DENOM
TEMP=CON(I,J,K)+AIP(I,J,K)*F(I+1,J,K,N)+AIM(I,J,K)*F(I-1,J,K,N)
&+AJP(I,J,K)*F(I,J+1,K,N)+AJM(I,J,K)*F(I,J-1,K,N)
QT(K)=(TEMP+AKM(I,J,K)*QT(K-1))/DENOM
570 CONTINUE
DO 580 KK=KST,N2
K=KT1-KK
580 F(I,J,K,N)=F(I,J,K+1,N)*PT(K)+QT(K)
590 CONTINUE
************************************************************************
999 CONTINUE
DO 600 K=2,N2
DO 600 J=2,M2
DO 600 I=2,L2
CON(I,J,K)=0.
AP(I,J,K)=0.
600 CONTINUE
RETURN
END
SUBROUTINE SETUP
************************************************************************
DOUBLE PRECISION TITLE
PARAMETER (ID=25,JD=55,KD=25)
LOGICAL LSOLVE,LPRINT,LBLK,LSTOP
COMMON F(ID,JD,KD,10),P(ID,JD,KD),RHO(ID,JD,KD),GAM(ID,JD,KD),
&CON(ID,JD,KD),AP(ID,JD,KD),AKP(ID,JD,KD),AKM(ID,JD,KD),
&AIP(ID,JD,KD),AIM(ID,JD,KD),AJP(ID,JD,KD),AJM(ID,JD,KD),
&X(ID),XU(ID),XDIF(ID),XCV(ID),XCVS(ID),XCVI(ID),XCVIP(ID),
&Y(JD),YV(JD),YDIF(JD),YCV(JD),YCVS(JD),YCVJ(JD),YCVJP(JD),
&Z(KD),ZW(KD),ZDIF(KD),ZCV(KD),ZCVS(KD),ZCVK(KD),ZCVKP(KD),
&AX(JD,KD),AY(ID,KD),AZ(ID,JD)
COMMON DU(ID,JD,KD),DV(ID,JD,KD),DW(ID,JD,KD),
& FX(ID),FXM(ID),FY(JD),FYM(JD),FZ(KD),FZM(KD),PT(ID),QT(ID)
COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,N1,N2,N3,
& IST,JST,KST,ITER,LAST,TITLE(13),RELAX(13),TIME,DT,XL,YL,ZL,
& IPREF,JPREF,KPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),
& RHOCON
COMMON/CNTL/LSTOP
COMMON/SORC/SMAX,SSUM
COMMON/COEF/FLOW,DIFF,ACOF
DIMENSION U(ID,JD,KD),V(ID,JD,KD),W(ID,JD,KD),PC(ID,JD,KD),
&T(ID,JD,KD)
DIMENSION XPLUS(ID,JD,KD),AKE(ID,JD,KD)
& ,DIS(ID,JD,KD),AMUT(ID,JD,KD),GEN(ID,JD,KD)
EQUIVALENCE (F(1,1,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)),
&(F(1,1,1,3),W(1,1,1)),(F(1,1,1,4),PC(1,1,1)),(F(1,1,1,5),T(1,1,1))
&,(F(1,1,1,6),AKE(1,1,1)),(F(1,1,1,7),DIS(1,1,1))
& ,(F(1,1,1,8),AMUT(1,1,1)),(F(1,1,1,9),GEN(1,1,1))
************************************************************************
1 FORMAT(//15X,'COMPUTATION IN CARTISIAN COORDINATES(SS)')
2 FORMAT(//15X,'COMPUTATION FOR AXISYMMETRICAL SITUATION')
3 FORMAT(//15X,' COMPUTATION IN POLAR COORDINATES ')
4 FORMAT(1X,14X,40(1H*),//)
*
DATA NFMAX,NP,NRHO,NGAM/10,11,12,13/
DATA LSTOP,LSOLVE,LPRINT/24*.FALSE./
DATA LBLK/10*.FALSE./
DATA MODE,LAST,TIME,ITER/1,5,0.,0/
DATA RELAX,NTIMES/13*1.,10*1/
DATA DT,IPREF,JPREF,KPREF,RHOCON/1.E+10,1,1,1,1./
*-----------------------------------------------------------------------
ENTRY SETUP1
L2=L1-1
L3=L2-1
M2=M1-1
M3=M2-1
N2=N1-1
N3=N2-1
X(1)=XU(2)
DO 5 I=2,L2
5 X(I)=0.5*(XU(I+1)+XU(I))
X(L1)=XU(L1)
Y(1)=YV(2)
DO 10 J=2,M2
10 Y(J)=0.5*(YV(J+1)+YV(J))
Y(M1)=YV(M1)
Z(1)=ZW(2)
DO 12 K=2,N2
12 Z(K)=0.5*(ZW(K+1)+ZW(K))
Z(N1)=ZW(N1)
DO 15 I=2,L1
15 XDIF(I)=X(I)-X(I-1)
DO 18 I=2,L2
18 XCV(I)=XU(I+1)-XU(I)
DO 20 I=3,L2
20 XCVS(I)=XDIF(I)
XCVS(3)=XCVS(3)+XDIF(2)
XCVS(L2)=XCVS(L2)+XDIF(L1)
DO 22 I=3,L3
XCVI(I)=0.5*XCV(I)
22 XCVIP(I)=XCVI(I)
XCVIP(2)=XCV(2)
XCVI(L2)=XCV(L2)
DO 35 J=2,M1
35 YDIF(J)=Y(J)-Y(J-1)
DO 40 J=2,M2
40 YCV(J)=YV(J+1)-YV(J)
DO 45 J=3,M2
45 YCVS(J)=YDIF(J)
YCVS(3)=YCVS(3)+YDIF(2)
YCVS(M2)=YCVS(M2)+YDIF(M1)
DO 48 J=3,M3
YCVJ(J)=0.5*YCV(J)
48 YCVJP(J)=YCVJ(J)
YCVJP(2)=YCV(2)
YCVJ(M2)=YCV(M2)
DO 50 K=2,N1
50 ZDIF(K)=Z(K)-Z(K-1)
DO 52 K=2,N2
52 ZCV(K)=ZW(K+1)-ZW(K)
DO 55 K=3,N2
55 ZCVS(K)=ZDIF(K)
ZCVS(3)=ZCVS(3)+ZDIF(2)
ZCVS(N2)=ZCVS(N2)+ZDIF(N1)
DO 58 K=3,N3
ZCVK(K)=0.5*ZCV(K)
58 ZCVKP(K)=ZCVK(K)
ZCVKP(2)=ZCV(2)
ZCVK(N2)=ZCV(N2)
DO 60 K=2,N2
DO 62 J=2,M2
62 AX(J,K)=YCV(J)*ZCV(K)
DO 64 I=2,L2
64 AY(I,K)=XCV(I)*ZCV(K)
60 CONTINUE
DO 66 J=2,M2
DO 66 I=2,L2
66 AZ(I,J)=XCV(I)*YCV(J)
DO 85 I=3,L2
FX(I)=0.5*XCV(I-1)/XDIF(I)
85 FXM(I)=1.-FX(I)
FX(2)=0.
FXM(2)=1.
FX(L1)=1.
FXM(L1)=0.
DO 90 J=3,M2
FY(J)=0.5*YCV(J-1)/YDIF(J)
90 FYM(J)=1.-FY(J)
FY(2)=0.
FYM(2)=1.
FY(M1)=1.
FYM(M1)=0.
DO 92 K=3,N2
FZ(K)=0.5*ZCV(K-1)/ZDIF(K)
92 FZM(K)=1.-FZ(K)
FZ(2)=0.
FZM(2)=1.
FZ(N1)=1.
FZM(N1)=0.
*---CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE----
DO 95 K=1,N1
DO 95 J=1,M1
DO 95 I=1,L1
PC(I,J,K)=0.
U(I,J,K)=0.
V(I,J,K)=0.
W(I,J,K)=0.
CON(I,J,K)=0.
AP(I,J,K)=0.
RHO(I,J,K)=RHOCON
P(I,J,K)=0.
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -