📄 simplec.f
字号:
IF (J .EQ. M2) GOTO 106
FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)*
& RHO(I-1,J))
GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+
& 1.0E-30)*XCVI(I)
GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)*
& GAM(I-1,J)+1.E-30)*XCVIP(I-1)
DIFF=RMN(J+1)*2.*(GM+GMM)
GO TO 107
106 FL=XCVI(I)*V(I,M1)*RHO(I,M1)
FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1)
DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1)
107 FLOW=RMN(J+1)*(FL+FLM)
CALL DIFLOW
AJM(I,J+1)=ACOF+AMAX1(0.,FLOW)
AJP(I,J)=AJM(I,J+1)-FLOW
VOL=YCVR(J)*XCVS(I)
APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1))
& /(XCVS(I)*DT)
AP(I,J)=AP(I,J)-APT
CON(I,J)=CON(I,J)+APT*U(I,J)
AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))
& /RELAX(NF)
CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J)
DU(I,J)=VOL/(XDIF(I)*SX(J))
CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J))
DU(I,J)=DU(I,J)/
&(AP(I,J)-AIM(I,J)-AIP(I,J)-AJM(I,J)-AJP(I,J)+1.0E-30)
103 CONTINUE
c******u方程余量数 ****************************************
c smaxu=0.
c do i=3,l2
c do j=2,m2
c smaxt=ap(i,j)*u(i,j)-aip(i,j)*u(i+1,j)
c + -aim(i,j)*u(i-1,j)-ajp(i,j)*u(i,j+1)-ajm(i,j)
c + *u(i,j-1)-con(i,j)
c smaxu=amax1(smaxu,ABS(smaxt))
c smaxu=smaxu+smaxt*smaxt
c enddo
c enddo
c smaxu=sqrt(smaxu)
c**********************************************************
CALL SOLVE
100 CONTINUE
*---COEFFICIENTS FOR THE V EQUATION----
NF=2
IF(.NOT. LSOLVE(NF)) GO TO 200
IST=2
JST=3
CALL GAMSOR
REL=1.-RELAX(NF)
DO 202 I=2,L2
AREA=R(1)*XCV(I)
FLOW=AREA*V(I,2)*RHO(I,1)
DIFF=AREA*GAM(I,1)/YCV(2)
CALL DIFLOW
202 AJM(I,3)=ACOF+AMAX1(0.,FLOW)
DO 203 J=3,M2
FL=ARXJ(J)*U(2,J)*RHO(1,J)
FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1)
FLOW=FL+FLM
DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J))
CALL DIFLOW
AIM(2,J)=ACOF+AMAX1(0.,FLOW)
DO 203 I=2,L2
IF(I .EQ. L2)GO TO 204
FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)*
& RHO(I,J-1))
GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+
& 1.E-30)*ARXJ(J)
GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)*
& GAM(I,J-1)+1.0E-30)*ARXJP(J-1)
DIFF=2.*(GM+GMM)/SXMN(J)
GO TO 205
204 FL=ARXJ(J)*U(L1,J)*RHO(L1,J)
FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1)
DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))
205 FLOW=FL+FLM
CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
IF(J .EQ. M2) GO TO 206
AREA=R(J)*XCV(I)
FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J)
FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))*RMN(J+1)
FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I)
DIFF=AREA*GAM(I,J)/YCV(J)
GO TO 207
206 AREA=R(M1)*XCV(I)
FLOW=AREA*V(I,M1)*RHO(I,M1)
DIFF=AREA*GAM(I,M1)/YCV(M2)
207 CALL DIFLOW
AJM(I,J+1)=ACOF+AMAX1(0.,FLOW)
AJP(I,J)=AJM(I,J+1)-FLOW
VOL=YCVRS(J)*XCV(I)
SXT=SX(J)
IF(J .EQ. M2) SXT=SX(M1)
SXB=SX(J-1)
IF(J .EQ. 3) SXB=SX(1)
APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)*
& 0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT)
AP(I,J)=AP(I,J)-APT
CON(I,J)=CON(I,J)+APT*V(I,J)
AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))
& /RELAX(NF)
CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J)
DV(I,J)=VOL/YDIF(J)
CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J))
DV(I,J)=DV(I,J)
&/(AP(I,J)-AIM(I,J)-AIP(I,J)-AJM(I,J)-AJP(I,J)+1.0e-30)
203 CONTINUE
c******V方程余量的范数*************************************
c smaxv=0.
c do i=2,l2
c do j=3,m2
c smaxt=ap(i,j)*v(i,j)-aip(i,j)*v(i+1,j)
c + -aim(i,j)*v(i-1,j)-ajp(i,j)*v(i,j+1)-ajm(i,j)
c + *v(i,j-1)-con(i,j)
c smaxv=amax1(smaxv,ABS(smaxt))
c smaxv=smaxv+smaxt*smaxt
c enddo
c enddo
c smaxv=sqrt(smaxv)
c**********************************************************
CALL SOLVE
200 CONTINUE
*---COEFIICIENTS FOR THE PRESSURE CORRECTION EQUATION----
NF=3
IF(.NOT. LSOLVE(NF)) GO TO 500
IST=2
JST=2
CALL GAMSOR
SMAX=0.
SSUM=0.
DO 410 J=2,M2
DO 410 I=2,L2
VOL=YCVR(J)*XCV(I)
410 CON(I,J)=CON(I,J)*VOL
DO 402 I=2,L2
ARHO=R(1)*XCV(I)*RHO(I,1)
CON(I,2)=CON(I,2)+ARHO*V(I,2)
402 AJM(I,2)=0.
DO 403 J=2,M2
ARHO=ARX(J)*RHO(1,J)
CON(2,J)=CON(2,J)+ARHO*U(2,J)
AIM(2,J)=0.
DO 403 I=2,L2
IF(I .EQ. L2) GO TO 404
ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
FLOW=ARHO*U(I+1,J)
CON(I,J)=CON(I,J)-FLOW
CON(I+1,J)=CON(I+1,J)+FLOW
AIP(I,J)=ARHO*DU(I+1,J)
*CCC
AIM(I+1,J)=AIP(I,J)
GO TO 405
404 ARHO=ARX(J)*RHO(L1,J)
CON(I,J)=CON(I,J)-ARHO*U(L1,J)
AIP(I,J)=0.
405 IF(J .EQ. M2) GO TO 406
ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
FLOW=ARHO*V(I,J+1)
CON(I,J)=CON(I,J)-FLOW
CON(I,J+1)=CON(I,J+1)+FLOW
AJP(I,J)=ARHO*DV(I,J+1)
AJM(I,J+1)=AJP(I,J)
GO TO 407
406 ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
CON(I,J)=CON(I,J)-ARHO*V(I,M1)
AJP(I,J)=0.
407 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)
PC(I,J)=0.
SMAX=AMAX1(SMAX,ABS(CON(I,J)))
SSUM=SSUM+CON(I,J)
403 CONTINUE
CALL SOLVE
*---COMEE HERE TO CORRECT THE PRESSURE AND VELOCITIES
DO 501 J=2,M2
DO 501 I=2,L2
P(I,J)=P(I,J)+PC(I,J)*RELAX(NP)
IF(I .NE. 2) U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J))
IF(J .NE. 2) V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J))
501 CONTINUE
500 CONTINUE
*---COEFFICIENTS FOR OTHER EQUATIONS----
IST=2
JST=2
DO 600 N=4,NFMAX
NF=N
IF(.NOT. LSOLVE(NF)) GO TO 600
CALL GAMSOR
REL=1.-RELAX(NF)
DO 602 I=2,L2
AREA=R(1)*XCV(I)
FLOW=AREA*V(I,2)*RHO(I,1)
DIFF=AREA*GAM(I,1)/YDIF(2)
CALL DIFLOW
602 AJM(I,2)=ACOF+AMAX1(0.,FLOW)
DO 603 J=2,M2
FLOW=ARX(J)*U(2,J)*RHO(1,J)
DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J))
CALL DIFLOW
AIM(2,J)=ACOF+AMAX1(0.,FLOW)
DO 603 I=2,L2
IF(I .EQ. L2) GO TO 604
FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))
DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+
& XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J))
GO TO 605
604 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J))
605 CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
AREA=RMN(J+1)*XCV(I)
IF(J .EQ. M2) GO TO 606
FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))
DIFF=AREA*2.*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+
& YCV(J+1)*GAM(I,J)+1.0E-30)
GO TO 607
606 FLOW=AREA*V(I,M1)*RHO(I,M1)
DIFF=AREA*GAM(I,M1)/YDIF(M1)
607 CALL DIFLOW
AJM(I,J+1)=ACOF+AMAX1(0.,FLOW)
AJP(I,J)=AJM(I,J+1)-FLOW
VOL=YCVR(J)*XCV(I)
APT=RHO(I,J)/DT
AP(I,J)=AP(I,J)-APT
CON(I,J)=CON(I,J)+APT*F(I,J,NF)
AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))
& /RELAX(NF)
CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*F(I,J,NF)
603 CONTINUE
CALL SOLVE
600 CONTINUE
*---------------
TIME=TIME+DT
ITER=ITER+1
IF(ITER.GE.LAST) LSTOP=.TRUE.
RETURN
END
*-----------------------------------------------------------------------
SUBROUTINE SUPPLY
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
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))
************************************************************************
10 FORMAT(1X,20(1H*),3X,A10,3X,20(1H*))
20 FORMAT(1X,4H I =,I6,6I9)
30 FORMAT(1X,1HJ)
40 FORMAT(1X,I2,3X,1P7E9.2)
50 FORMAT(1X,1H )
51 FORMAT(1X,' I =',2X,7(I4,5X))
52 FORMAT(1X,' X =',1P7E9.2)
53 FORMAT(1X,'TH =',1P7E9.2)
54 FORMAT(1X,'J =',2X,7(I4,5X))
55 FORMAT(1X,'Y =',1P7E9.2)
************************************************************************
ENTRY UGRID
XU(2)=0.
DX=XL/FLOAT(L1-2)
DO 1 I=3,L1
1 XU(I)=XU(I-1)+DX
YV(2)=0.
DY=YL/FLOAT(M1-2)
DO 2 J=3,M1
2 YV(J)=YV(J-1)+DY
RETURN
************************************************************************
ENTRY PRINT
IF(.NOT. LPRINT(3)) GO TO 80
*---CALCULATE THE STREAM FUNTION----------------------------------------
F(2,2,3)=0.
DO 82 I=2,L1
IF(I .NE. 2) F(I,2,3)=F(I-1,2,3)-RHO(I-1,1)*V(I-1,2)
& *R(1)*XCV(I-1)
DO 82 J=3,M1
RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1)
82 F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1)
80 CONTINUE
*
IF( .NOT. LPRINT(NP)) GO TO 90
*
*---CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION
DO 91 J=2,M2
P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3)
91 P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2)
DO 92 I=2,L2
P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3)
92 P(I,M1)=(P(I,M2)*YCVS(M2)-P(I,M3)*YDIF(M1))/YDIF(M2)
P(1,1)=P(2,1)+P(1,2)-P(2,2)
P(L1,1)=P(L2,1)+P(L1,2)-P(L2,2)
P(1,M1)=P(2,M1)+P(1,M2)-P(2,M2)
P(L1,M1)=P(L2,M1)+P(L1,M2)-P(L2,M2)
PREF=P(IPREF,JPREF)
DO 93 J=1,M1
DO 93 I=1,L1
93 P(I,J)=P(I,J)-PREF
90 CONTINUE
return
*
IF(TIME.GT.0.5*DT) GOTO 320
WRITE (8,50)
IEND=0
301 IF(IEND .EQ. L1) GO TO 310
IBEG=IEND+1
IEND=IEND+7
IEND=MIN0(IEND,L1)
WRITE (8,50)
WRITE(8,51) (I,I=IBEG,IEND)
IF(MODE .EQ. 3) GO TO 302
WRITE(8,52) (X(I),I=IBEG,IEND)
GO TO 303
302 WRITE (8,53) (X(I),I=IBEG,IEND)
303 GO TO 301
310 JEND=0
WRITE(8,50)
311 IF(JEND .EQ. M1) GO TO 320
JBEG=JEND+1
JEND=JEND+7
JEND=MIN0(JEND,M1)
WRITE(8,50)
WRITE(8,54) (J,J=JBEG,JEND)
WRITE(8,55) (Y(J),J=JBEG,JEND)
GO TO 311
320 CONTINUE
*-------------------------
JT1=M2+JST
DO 332 J=2,M2
JJ=M2-J+2
DO 332 I=2,L2
332 WRITE(3,321) X(I),Y(JJ),U(I+1,JJ),V(I,JJ+1)
321 FORMAT(2X,F7.3,2X,F7.3,2X,E15.3,2X,E15.3)
*-------------------------
DO 999 N=1,NGAM
NF=N
IF(.NOT. LPRINT(NF)) GO TO 999
WRITE(8,50)
WRITE(8,10) TITLE(NF)
IFST=1
JFST=1
IF(NF .EQ. 1 .OR. NF .EQ. 3) IFST=2
IF(NF .EQ. 2 .OR. NF .EQ. 3) JFST=2
IBEG=IFST-7
110 CONTINUE
IBEG=IBEG+7
IEND=IBEG+6
IEND=MIN0(IEND,L1)
WRITE(8,50)
WRITE(8,20) (I,I=IBEG,IEND)
WRITE(8,30)
JFL=JFST+M1
DO 115 JJ=JFST,M1
J=JFL-JJ
WRITE(8,40) J,(F(I,J,NF),I=IBEG,IEND)
115 CONTINUE
IF(IEND .LT. L1) GO TO 110
999 CONTINUE
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -