📄 simplec.f
字号:
$DEBUG
****************************************************************************
*----------------------------MAIN PROGRAM-----------------------------------
****************************************************************************
LOGICAL LSTOP
COMMON/CNTL/LSTOP
****************************************************************************
real(4) kk, ta(2)
*********************************************************************
OPEN(8,FILE='RESULT.DAT')
CALL SETUP0
CALL GRID
CALL SETUP1
CALL START
10 CALL DENSE
DO WHILE(.NOT.LSTOP)
CALL DENSE
CALL BOUND
CALL OUTPUT
CALL SETUP2
ENDDO
kk = ETIME(TA)
write(8,*) 'Program has used', KK, 'seconds of CPU time.'
write(8,*) ' This includes', TA(1), 'seconds of user time and',
+ TA(2), 'seconds of system time.'
PRINT*, 'Program has used', KK, 'seconds of CPU time.'
PRINT*,' This includes', TA(1), 'seconds of user time and',
+ TA(2), 'seconds of system time.'
CLOSE(8)
END
*---------------------------------------------------------------------------
SUBROUTINE DIFLOW
IMPLICIT REAL*8(A-H,O-Z)
****************************************************************************
COMMON/COEF/FLOW,DIFF,ACOF
****************************************************************************
ACOF=DIFF
IF(FLOW .EQ.0.0)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
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER(NI=300,NJ=300,NIJ=NI,NFMAX=10,NFX3=NFMAX+3)
****************************************************************************
DOUBLE PRECISION TITLE
LOGICAL LSOLVE,LPRINT,LBLK,LSTOP
COMMON F(NI,NJ,NFMAX),P(NI,NJ),RHO(NI,NJ),GAM(NI,NJ),CON(NI,NJ),
1 AIP(NI,NJ),AIM(NI,NJ),AJP(NI,NJ),AJM(NI,NJ),AP(NI,NJ),
2 X(NI),XU(NI),XDIF(NI),XCV(NI),XCVS(NI),
3 Y(NJ),YV(NJ),YDIF(NJ),YCV(NJ),YCVS(NJ),
4 YCVR(NJ),YCVRS(NJ),ARX(NJ),ARXJ(NJ),ARXJP(NJ),
5 R(NJ),RMN(NJ),SX(NJ),SXMN(NJ),XCVI(NI),XCVIP(NI)
COMMON DU(NI,NJ),DV(NI,NJ), FV(NI),FVP(NI),
1 FX(NI),FXM(NI),FY(NJ),FYM(NJ),PT(NIJ),QT(NIJ)
COMMON/INDX/NF,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,
1IST,JST,ITER,LAST,TITLE(NFX3),RELAX(NFX3),TIME,DT,XL,YL,
2IPREF,JPREF,LSOLVE(NFX3),LPRINT(NFX3),LBLK(NFX3),MODE
3,NTIMES(NFX3),RHOCON
****************************************************************************
ISTF=IST-1
JSTF=JST-1
IT1=L2+IST
IT2=L3+IST
JT1=M2+JST
JT2=M3+JST
****************************************************************************
DO 999 NT=1,NTIMES(NF)
DO 999 N=NF,NF
*---------------------------------------------------------------------------
IF(.NOT. LBLK(NF)) GO TO 10
PT(ISTF)=0.
QT(ISTF)=0.
DO 11 I=IST,L2
BL=0.
BLP=0.
BLM=0.
BLC=0.
DO 12 J=JST,M2
BL=BL+AP(I,J)
IF(J .NE. M2) BL=BL-AJP(I,J)
IF(J .NE. JST) BL=BL-AJM(I,J)
BLP=BLP+AIP(I,J)
BLM=BLM+AIM(I,J)
BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
& +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)
12 CONTINUE
DENOM=BL-PT(I-1)*BLM
DENO=1.E15
IF(ABS(DENOM/BL) .LT. 1.E-10) DENOM=1.E20*DENO
PT(I)=BLP/DENOM
QT(I)=(BLC+BLM*QT(I-1))/DENOM
11 CONTINUE
BL=0.
DO 13 II=IST,L2
I=IT1-II
BL=BL*PT(I)+QT(I)
DO 13 J=JST,M2
13 F(I,J,N)=F(I,J,N)+BL
*---------------------------------------------------------------------------
PT(JSTF)=0.
QT(JSTF)=0.
DO 21 J=JST,M2
BL=0.
BLP=0.
BLM=0.
BLC=0.
DO 22 I=IST,L2
BL=BL+AP(I,J)
IF(I .NE. L2) BL=BL-AIP(I,J)
IF(I .NE. IST) BL=BL-AIM(I,J)
BLP=BLP+AJP(I,J)
BLM=BLM+AJM(I,J)
BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
& +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)
22 CONTINUE
DENOM=BL-PT(J-1)*BLM
IF (ABS(DENOM/BL) .LT. 1E-10) DENOM=1.E20*DENO
PT(J)=BLP/DENOM
QT(J)=(BLC+BLM*QT(J-1))/DENOM
21 CONTINUE
BL=0.
DO 23 JJ=JST,M2
J=JT1-JJ
BL=BL*PT(J)+QT(J)
DO 23 I=IST,L2
23 F(I,J,N)=F(I,J,N)+BL
10 CONTINUE
*-----------------------------------------------------------------------
DO 90 J=JST,M2
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J,N)
DO 70 I=IST,L2
DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
PT(I)=AIP(I,J)/DENOM
TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)
QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
70 CONTINUE
DO 80 II=IST,L2
I=IT1-II
80 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)
90 CONTINUE
*-----------------------------------------------------------------------
DO 190 JJ=JST,M3
J=JT2-JJ
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J,N)
DO 170 I=IST,L2
DENOM=AP(I,J)-PT(I-1)*AIM(I,J)
PT(I)=AIP(I,J)/DENOM
TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)
QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM
170 CONTINUE
DO 180 II=IST,L2
I=IT1-II
180 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)
190 CONTINUE
*-----------------------------------------------------------------------
DO 290 I=IST,L2
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF,N)
DO 270 J=JST,M2
DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
PT(J)=AJP(I,J)/DENOM
TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
270 CONTINUE
DO 280 JJ=JST,M2
J=JT1-JJ
280 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)
290 CONTINUE
*-----------------------------------------------------------------------
DO 390 II=IST,L3
I=IT2-II
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF,N)
DO 370 J=JST,M2
DENOM=AP(I,J)-PT(J-1)*AJM(I,J)
PT(J)=AJP(I,J)/DENOM
TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)
QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM
370 CONTINUE
DO 380 JJ=JST,M2
J=JT1-JJ
380 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)
390 CONTINUE
************************************************************************
999 CONTINUE
DO 400 J=2,M2
DO 400 I=2,L2
CON(I,J)=0.
AP(I,J)=0.
400 CONTINUE
RETURN
END
************************************************************************
SUBROUTINE SETUP
IMPLICIT REAL*8(A-H,O-Z)
************************************************************************
PARAMETER(NI=300,NJ=300,NIJ=NI,NFMAX=10,NFX3=NFMAX+3)
DOUBLE PRECISION TITLE
LOGICAL LSOLVE,LPRINT,LBLK,LSTOP
COMMON F(NI,NJ,NFMAX),P(NI,NJ),RHO(NI,NJ),GAM(NI,NJ),CON(NI,NJ),
1 AIP(NI,NJ),AIM(NI,NJ),AJP(NI,NJ),AJM(NI,NJ),AP(NI,NJ),
2 X(NI),XU(NI),XDIF(NI),XCV(NI),XCVS(NI),
3 Y(NJ),YV(NJ),YDIF(NJ),YCV(NJ),YCVS(NJ),
4 YCVR(NJ),YCVRS(NJ),ARX(NJ),ARXJ(NJ),ARXJP(NJ),
5 R(NJ),RMN(NJ),SX(NJ),SXMN(NJ),XCVI(NI),XCVIP(NI)
COMMON DU(NI,NJ),DV(NI,NJ), FV(NI),FVP(NI),
1 FX(NI),FXM(NI),FY(NJ),FYM(NJ),PT(NIJ),QT(NIJ)
COMMON/INDX/NF,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,
1IST,JST,ITER,LAST,TITLE(NFX3),RELAX(NFX3),TIME,DT,XL,YL,
2IPREF,JPREF,LSOLVE(NFX3),LPRINT(NFX3),LBLK(NFX3),MODE
3,NTIMES(NFX3),RHOCON
COMMON /CNTL/LSTOP
COMMON /SORC/SMAX,SSUM
COMMON /COEF/FLOW,DIFF,ACOF
DIMENSION U(NI,NJ),V(NI,NJ),PC(NI,NJ)
EQUIVALENCE (F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1))
************************************************************************
1 FORMAT(//15X,'COMPUTATION IN CARTISIAN COORDINATES')
2 FORMAT(//15X,'COMPUTATION FOR AXISYMMETRICAL SITUATION')
3 FORMAT(//15X,' COMPUTATION IN POLAR COORDINATES ')
4 FORMAT(1X,14X,40(1H*),//)
*-----------------------------------------------------------------------
ENTRY SETUP0
NP=NFMAX+1
NRHO=NP+1
NGAM=NRHO+1
LSTOP=.FALSE.
DO 779 I=1,10
LSOLVE(I)=.FALSE.
LBLK(I)=.TRUE.
779 NTIMES(I)=1
DO 889 I=1,13
LPRINT(I)=.FALSE.
889 RELAX(I)=1.
MODE=1
LAST=5
TIME=0.
ITER=0
DT=1.0E+10
IPREF=1
JPREF=1
RHOCON=1
RETURN
*-----------------------------------------------------------------------
ENTRY SETUP1
L2=L1-1
L3=L2-1
M2=M1-1
M3=M2-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)
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)
IF (MODE .NE. 1) GO TO 55
DO 52 J=1,M1
RMN(J)=1.
52 R(J)=1.
GO TO 56
55 DO 50 J=2,M1
50 R(J)=R(J-1)+YDIF(J)
RMN(2)=R(1)
DO 60 J=3,M2
60 RMN(J)=RMN(J-1)+YCV(J-1)
RMN(M1)=R(M1)
56 CONTINUE
DO 57 J=1,M1
SX(J)=1.
SXMN(J)=1.
IF(MODE .NE. 3) GO TO 57
SX(J)=R(J)
IF(J .NE. 1) SXMN(J)=RMN(J)
57 CONTINUE
DO 62 J=2,M2
YCVR(J)=R(J)*YCV(J)
ARX(J)=YCVR(J)
IF (MODE .NE. 3) GO TO 62
ARX(J)=YCV(J)
62 CONTINUE
DO 64 J=4,M3
64 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J)
YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3)
YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2)
IF(MODE .NE. 2) GO TO 67
DO 65 J=3,M3
ARXJ(J)=0.25*(1.+RMN(J)/R(J))*ARX(J)
65 ARXJP(J)=ARX(J)-ARXJ(J)
GO TO 68
67 DO 66 J=3,M3
ARXJ(J)=0.5*ARX(J)
66 ARXJP(J)=ARXJ(J)
68 ARXJP(2)=ARX(2)
ARXJ(M2)=ARX(M2)
DO 70 J=3,M3
FV(J)=ARXJP(J)/ARX(J)
70 FVP(J)=1.-FV(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.
*---CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE----
DO 95 J=1,M1
DO 95 I=1,L1
PC(I,J)=0.
U(I,J)=0.
V(I,J)=0.
CON(I,J)=0.
AP(I,J)=0.
C RS(I,J)=0.
RHO(I,J)=RHOCON
P(I,J)=0.
95 CONTINUE
IF(MODE .EQ. 1) WRITE(8,1)
IF(MODE .EQ. 2) WRITE(8,2)
IF(MODE .EQ. 3) WRITE(8,3)
WRITE(8,4)
RETURN
*----------------------------------------------------------------------
ENTRY SETUP2
*---COEFFICIENTS FOR THE U EQUATION----
NF=1
IF(.NOT. LSOLVE(NF)) GO TO 100
IST=3
JST=2
CALL GAMSOR
REL=1.-RELAX(NF)
DO 102 I=3,L2
FL=XCVI(I)*V(I,2)*RHO(I,1)
FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1)
FLOW=R(1)*(FL+FLM)
DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2)
CALL DIFLOW
102 AJM(I,2)=ACOF+AMAX1(0.,FLOW)
DO 103 J=2,M2
FLOW=ARX(J)*U(2,J)*RHO(1,J)
DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J))
CALL DIFLOW
AIM(3,J)=ACOF+AMAX1(0.,FLOW)
DO 103 I=3,L2
IF(I .EQ. L2) GO TO 104
FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J))
FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
FLOW=ARX(J)*0.5*(FL+FLP)
DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J))
GO TO 105
104 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J))
105 CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -