2D code use simple algorithm

boling2000

I found this code when I was visiting other website.
The only thing is that there is no document for this code.
I hope it is useful for us:
___________________________________________________________________________


LOGICAL LSTOP
COMMON/CNTL/LSTOP *-----------------------------------------------------------------------
CALL USER(1)
CALL SETUP(1)
CALL USER(2) 100 CALL USER(3)
CALL USER(4)
CALL USER(5)
IF(LSTOP) STOP
CALL SETUP(2)
GOTO 100
END *=====================================================================
SUBROUTINE DIFLOW *-----------------------------------------------------------------------
COMMON/COEF/FLOW,DIFF,ACOF *-----------------------------------------------------------------------
ACOF=DIFF
IF(FLOW.EQ.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 *-----------------------------------------------------------------------
COMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22,22),
+ AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22,22),
+ X(22),XU(22),XDIF(22),XCV(22),XCVS(22),
+ Y(22),YV(22),YDIF(22),YCV(22),YCVS(22),
+ YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22),
+ R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22),
+ DU(22,22),DV(22,22),FV(22),FVP(22),
+ FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22)
*----------- Arrays U, V and PC may be absent in "SOLVE" -------------
DIMENSION U(22,22),V(22,22),PC(22,22)
EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) *-----------------------------------------------------------------------
CHARACTER*10 TITLE
COMMON/A/TITLE(13)
LOGICAL LSOLVE,LPRINT,LBLK ! &,LSTOP
COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,
+ IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL,
+ IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON
*-----------------------------------------------------------------------
ISTF=IST-1
JSTF=JST-1
IT1=L2+IST
IT2=L3+IST
JT1=M2+JST
JT2=M3+JST *-----------------------------------------------------------------------
DO 100 NT=1,NTIMES(NF)
DO 100 N=NF,NF *-----------------------------------------------------------------------
IF(.NOT.LBLK(NF)) GOTO 110
PT(ISTF)=0.
QT(ISTF)=0.
DO 120 I=IST,L2
BL=0.
BLP=0.
BLM=0.
BLC=0.
DO 130 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) 130 CONTINUE
DENOM=BL-PT(I-1)*BLM
IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E35
PT(I)=BLP/DENOM
QT(I)=(BLC+BLM*QT(I-1))/DENOM 120 CONTINUE
BL=0.
DO 140 II=IST,L2
I=IT1-II
BL=BL*PT(I)+QT(I)
DO 140 J=JST,M2 140 F(I,J,N)=F(I,J,N)+BL *-----------------------------------------------------------------------
PT(JSTF)=0.
QT(JSTF)=0.
DO 150 J=JST,M2
BL=0.
BLP=0.
BLM=0.
BLC=0.
DO 160 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) 160 CONTINUE
DENOM=BL-PT(J-1)*BLM
IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E+35
PT(J)=BLP/DENOM
QT(J)=(BLC+BLM*QT(J-1))/DENOM 150 CONTINUE
BL=0.
DO 170 JJ=JST,M2
J=JT1-JJ
BL=BL*PT(J)+QT(J)
DO 170 I=IST,L2 170 F(I,J,N)=F(I,J,N)+BL 110 CONTINUE *-----------------------------------------------------------------------
DO 180 J=JST,M2
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J,N)
DO 190 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 190 CONTINUE
DO 200 II=IST,L2
I=IT1-II 200 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) 180 CONTINUE *-----------------------------------------------------------------------
DO 210 JJ=JST,M3
J=JT2-JJ
PT(ISTF)=0.
QT(ISTF)=F(ISTF,J,N)
DO 220 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 220 CONTINUE
DO 230 II=IST,L2
I=IT1-II 230 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) 210 CONTINUE *-----------------------------------------------------------------------
DO 240 I=IST,L2
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF,N)
DO 250 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 250 CONTINUE
DO 260 JJ=JST,M2
J=JT1-JJ 260 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) 240 CONTINUE *-----------------------------------------------------------------------
DO 270 II=IST,L3
I=IT2-II
PT(JSTF)=0.
QT(JSTF)=F(I,JSTF,N)
DO 280 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 280 CONTINUE
DO 290 JJ=JST,M2
J=JT1-JJ 290 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) 270 CONTINUE *----------------------------------------------------------------------- 100 CONTINUE
DO 300 J=2,M2
DO 300 I=2,L2
CON(I,J)=0.
AP(I,J)=0. 300 CONTINUE
RETURN
END *=====================================================================
SUBROUTINE SETUP(K) *-----------------------------------------------------------------------
COMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22,22),
+ AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22,22),
+ X(22),XU(22),XDIF(22),XCV(22),XCVS(22),
+ Y(22),YV(22),YDIF(22),YCV(22),YCVS(22),
+ YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22),
+ R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22),
+ DU(22,22),DV(22,22),FV(22),FVP(22),
+ FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22)
*----------- Arrays U, V and PC may be absent in "SOLVE" -------------
DIMENSION U(22,22),V(22,22),PC(22,22)
EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) *-----------------------------------------------------------------------
CHARACTER*10 TITLE
COMMON/A/TITLE(13)
LOGICAL LSOLVE,LPRINT,LBLK,LSTOP
COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,
+ IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL,
+ IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON
COMMON/CNTL/LSTOP
COMMON/SORC/SMAX,SSUM
COMMON/COEF/FLOW,DIFF,ACOF *----------------------------------------------------------------------- 10 FORMAT(/15X,' COMPUTATION IN CARTESIAN COORDINATES') 20 FORMAT(/15X,'COMPUTATION FOR AXISYMMETRIC SITUATION') 30 FORMAT(/15X,' COMPUTATION IN POLAR COORDINATES') 40 FORMAT(14X,40(1H*),/) *-----------------------------------------------------------------------
GOTO (1000,2000) K *----------------------------------------------------------------------- * ENTRY SETUP1 1000 NFMAX=10
NP=11
NRHO=12
NGAM=13
LSTOP=.FALSE.
DO 1010 I=1,10
LSOLVE(I)=.FALSE.
NTIMES(I)=1 1010 LBLK(I)=.TRUE.
DO 1020 I=1,13
LPRINT(I)=.FALSE. 1020 RELAX(I)=1.
LAST=5
TIME=0.
ITER=0
DT=1.E+10
IPREF=1
JPREF=1
RHOCON=1. *-----------------------------------------------------------------------
L2=L1-1
L3=L2-1
M2=M1-1
M3=M2-1
X(1)=XU(2)
DO 1030 I=2,L2 1030 X(I)=0.5*(XU(I)+XU(I+1))
X(L1)=XU(L1)
Y(1)=YV(2)
DO 1040 J=2,M2 1040 Y(J)=0.5*(YV(J+1)+YV(J))
Y(M1)=YV(M1)
DO 1050 I=2,L1 1050 XDIF(I)=X(I)-X(I-1)
DO 1060 I=2,L2 1060 XCV(I)=XU(I+1)-XU(I)
DO 1070 I=3,L2 1070 XCVS(I)=XDIF(I)
XCVS(L2)=XCVS(L2)+XDIF(L1)
XCVS(3)=XCVS(3)+XDIF(2)
DO 1080 I=3,L3
XCVI(I)=0.5*XCV(I) 1080 XCVIP(I)=XCVI(I)
XCVIP(2)=XCV(2)
XCVI(L2)=XCV(L2)
DO 1090 J=2,M1 1090 YDIF(J)=Y(J)-Y(J-1)
DO 1100 J=2,M2 1100 YCV(J)=YV(J+1)-YV(J)
DO 1110 J=3,M2 1110 YCVS(J)=YDIF(J)
YCVS(3)=YCVS(3)+YDIF(2)
YCVS(M2)=YCVS(M2)+YDIF(M1)
IF(MODE.NE.1) GOTO 1120
DO 1130 J=1,M1
RMN(J)=1.0 1130 R(J)=1.0
GOTO 1140 1120 DO 1150 J=2,M1 1150 R(J)=R(J-1)+YDIF(J)
RMN(2)=R(1)
DO 1160 J=3,M2 1160 RMN(J)=RMN(J-1)+YCV(J-1)
RMN(M1)=R(M1) 1140 CONTINUE
DO 1170 J=1,M1
SX(J)=1.0
SXMN(J)=1.0
IF(MODE.NE.3) GOTO 1170
SX(J)=R(J)
IF(J.NE.1) SXMN(J)=RMN(J) 1170 CONTINUE
DO 1180 J=2,M2
YCVR(J)=R(J)*YCV(J)
ARX(J)=YCVR(J)
IF(MODE.NE.3) GOTO 1180
ARX(J)=YCV(J) 1180 CONTINUE
DO 1190 J=4,M3 1190 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J)
YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3)
YCVRS(M2)=.5*(R(M1)+R(M3))*YCVS(M2)
IF(MODE.NE.2) GOTO 1200
DO 1210 J=3,M3
ARXJ(J)=.25*(1.+RMN(J)/R(J))*ARX(J) 1210 ARXJP(J)=ARX(J)-ARXJ(J)
GOTO 1220 1200 DO 1230 J=3,M3
ARXJ(J)=.5*ARX(J) 1230 ARXJP(J)=ARXJ(J) 1220 ARXJP(2)=ARX(2)
ARXJ(M2)=ARX(M2)
DO 1240 J=3,M3
FV(J)=ARXJP(J)/ARX(J) 1240 FVP(J)=1.0-FV(J)
DO 1250 I=3,L2
FX(I)=.5*XCV(I-1)/XDIF(I) 1250 FXM(I)=1.-FX(I)
FX(2)=0.0
FXM(2)=1.
FX(L1)=1.
FXM(L1)=0.
DO 1260 J=3,M2
FY(J)=.5*YCV(J-1)/YDIF(J) 1260 FYM(J)=1.-FY(J)
FY(2)=0.
FYM(2)=1.0
FY(M1)=1.0
FYM(M1)=0. *---------- CON,AP,U,V,RHO,PC,P ARRAYS ARE INITIALIZED HERE ----------
DO 1270 J=1,M1
DO 1270 I=1,L1
PC(I,J)=0.
U(I,J)=0.
V(I,J)=0.
CON(I,J)=0.
AP(I,J)=0.
RHO(I,J)=RHOCON
P(I,J)=0. 1270 CONTINUE
IF(MODE.EQ.1) WRITE(10,10)
IF(MODE.EQ.2) WRITE(10,20)
IF(MODE.EQ.3) WRITE(10,30)
WRITE(10,40)
RETURN *----------------------------------------------------------------------- * ENTRY SETUP2 *----------------- COEFFICIENTS FOR THE U EQUATION ----------------- 2000 NF=1
IF(.NOT.LSOLVE(NF)) GOTO 2110
IST=3
JST=2
CALL USER(6)
REL=1.-RELAX(NF)
DO 2120 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 2120 AJM(I,2)=ACOF+AMAX1(0.,FLOW)
DO 2130 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 2130 I=3,L2
IF(I.EQ.L2) GOTO 2140
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))
GOTO 2150 2140 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J)) 2150 CALL DIFLOW
AIM(I+1,J)=ACOF+AMAX1(0.,FLOW)
AIP(I,J)=AIM(I+1,J)-FLOW
IF(J.EQ.M2) GOTO 2160
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)
GOTO 2170 2160 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) 2170 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) AP0/DX/DY
AP(I,J)=AP(I,J)-APT AP->Sp
CON(I,J)=CON(I,J)+APT*U(I,J) CON->Sc must be given firstly
AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJM(I,J)+AJP(I,J))/
+ RELAX(NF)
CON(I,J)=CON(I,J)*VOL+AP(I,J)*U(I,J)*REL
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) 2130 CONTINUE
CALL SOLVE 2110 CONTINUE *----------------- COEFFICIENTS FOR THE V EQUATION -----------------
NF=2
IF(.NOT.LSOLVE(NF)) GOTO 2210
IST=2
JST=3
CALL USER(6)
REL=1.-RELAX(NF)
DO 2220 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 2220 AJM(I,3)=ACOF+AMAX1(0.,FLOW)
DO 2230 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 2230 I=2,L2
IF(I.EQ.L2) GOTO 2240
DIFF=2.*(GM+GMM)/SXMN(J)
GOTO 2250 2240 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)) 2250 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) GOTO 2260
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)
GOTO 2270 2260 AREA=R(M1)*XCV(I)
FLOW=AREA*V(I,M1)*RHO(I,M1)
DIFF=AREA*GAM(I,M1)/YCV(M2) 2270 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+AIM(I,J)+AIP(I,J)+AJM(I,J)+AJP(I,J))/
+ RELAX(NF)
CON(I,J)=CON(I,J)*VOL+AP(I,J)*V(I,J)*REL
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) 2230 CONTINUE
CALL SOLVE 2210 CONTINUE *--------- COEFFICIENT FOR THE PRESSURE CORRECTION EQUATION ----------
NF=3
IF(.NOT.LSOLVE(NF)) GOTO 2310
IST=2
JST=2
CALL USER(6)
SMAX=0.
SSUM=0.
DO 2320 J=2,M2
DO 2320 I=2,L2
VOL=YCVR(J)*XCV(I) 2320 CON(I,J)=CON(I,J)*VOL
DO 2330 I=2,L2
ARHO=R(1)*XCV(I)*RHO(I,1)
CON(I,2)=CON(I,2)+ARHO*V(I,2) 2330 AJM(I,2)=0.
DO 2340 J=2,M2
ARHO=ARX(J)*RHO(1,J)
CON(2,J)=CON(2,J)+ARHO*U(2,J)
AIM(2,J)=0.
DO 2340 I=2,L2
IF(I.EQ.L2) GOTO 2350
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)
AIM(I+1,J)=AIP(I,J)
GOTO 2360 2350 ARHO=ARX(J)*RHO(L1,J)
CON(I,J)=CON(I,J)-ARHO*U(L1,J)
AIP(I,J)=0. 2360 IF(J.EQ.M2) GOTO 2370
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
AJP(I,J)=ARHO*DV(I,J+1)
AJM(I,J+1)=AJP(I,J)
GOTO 2380 2370 ARHO=RMN(M1)*XCV(I)*RHO(I,M1)
CON(I,J)=CON(I,J)-ARHO*V(I,M1)
AJP(I,J)=0. 2380 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) 2340 CONTINUE
CALL SOLVE *--------- COME HERE TO CORRECT THE PRESSURE AND VELOCITIES ----------
DO 2390 J=2,M2
DO 2390 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)) 2390 CONTINUE 2310 CONTINUE *----------------- COEFFICIENTS FOR OTHER EQUATIONS ------------------
IST=2
JST=2
DO 2410 N=4,NFMAX
NF=N
IF(.NOT.LSOLVE(NF)) GOTO 2410
CALL USER(6)
REL=1.-RELAX(NF)
DO 2420 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 2420 AJM(I,2)=ACOF+AMAX1(0.,FLOW)
DO 2430 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 2430 I=2,L2
IF(I.EQ.L2) GOTO 2440
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.E-30)*SX(J))
GOTO 2450 2440 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)
DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J)) 2450 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) GOTO 2460
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+1)*GAM(I,J)/(YCV(J)*GAM(I,J+1)+
+ YCV(J+1)*GAM(I,J)+1.E-30)
GOTO 2470 2460 FLOW=AREA*V(I,M1)*RHO(I,M1)
DIFF=AREA*GAM(I,M1)/YDIF(M1) 2470 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+AP(I,J)*F(I,J,NF)*REL 2430 CONTINUE
CALL SOLVE 2410 CONTINUE
TIME=TIME+DT
ITER=ITER+1
IF(ITER.EQ.LAST) LSTOP=.TRUE.
RETURN
END *=====================================================================
SUBROUTINE SUPPLY(K) *-----------------------------------------------------------------------
COMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22,22),
+ AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22,22),
+ X(22),XU(22),XDIF(22),XCV(22),XCVS(22),
+ Y(22),YV(22),YDIF(22),YCV(22),YCVS(22),
+ YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22),
+ R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22),
+ DU(22,22),DV(22,22),FV(22),FVP(22),
+ FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22)
*----------- Arrays U, V and PC may be absent in "SOLVE" -------------
DIMENSION U(22,22),V(22,22),PC(22,22)
EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) *-----------------------------------------------------------------------
CHARACTER*10 TITLE
COMMON/A/TITLE(13)
LOGICAL LSOLVE,LPRINT,LBLK ! &,LSTOP
COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,
+ IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL,
+ IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON
*----------------------------------------------------------------------- 110 FORMAT(1X,26(1H*),3X,A10,3X,26(1H*)) 120 FORMAT(1X,'I =',I8,6I10) 130 FORMAT(1X,'J') 140 FORMAT(1X,I2,2X,1P7E10.2) 150 FORMAT(1X,' ') 160 FORMAT(1X,'I =',I8,6I10) 170 FORMAT(1X,'X = ',1P7E10.2) 180 FORMAT(1X,'XU =',1P7E10.2) 190 FORMAT(1X,'TH =',1P7E10.2) 200 FORMAT(1X,'J =',I8,6I10) 210 FORMAT(1X,'Y = ',1P7E10.2) 220 FORMAT(1X,'YV =',1P7E10.2) *-----------------------------------------------------------------------
GOTO (1000,2000) K *----------------------------------------------------------------------- * ENTRY UGRID 1000 XU(2)=0.
DX=XL/FLOAT(L1-2)
DO 1010 I=3,L1 1010 XU(I)=XU(I-1)+DX
YV(2)=0.
DY=YL/FLOAT(M1-2)
DO 1020 J=3,M1 1020 YV(J)=YV(J-1)+DY
RETURN *----------------------------------------------------------------------- * ENTRY PRINT 2000 IF(.NOT.LPRINT(3)) GOTO 2010 *------------------ CALCULATE THE STREAM FUNCTION ------------------
F(2,2,3)=0.
DO 2020 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 2020 J=3,M1
RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1) 2020 F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1) 2010 CONTINUE
IF(.NOT.LPRINT(NP)) GOTO 2030 *----------- CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION -----------
DO 2040 J=2,M2
P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3) 2040 P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2)
DO 2050 I=2,L2
P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3) 2050 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 2060 J=1,M1
DO 2060 I=1,L1 2060 P(I,J)=P(I,J)-PREF 2030 CONTINUE *-----------------------------------------------------------------------
WRITE(10,150)
IEND=1 2080 IF(IEND.EQ.L1) GOTO 2070
IBEG=IEND+1
IEND=IEND+7
IEND=MIN0(IEND,L1)
WRITE(10,160) (I,I=IBEG,IEND)
WRITE(10,180) (XU(I),I=IBEG,IEND)
GOTO 2080 2070 WRITE(10,150)
JEND=1 2100 IF(JEND.EQ.M1) GOTO 2090
JBEG=JEND+1
JEND=JEND+7
JEND=MIN0(JEND,M1)
WRITE(10,200) (J,J=JBEG,JEND)
WRITE(10,220) (YV(J),J=JBEG,JEND)
GOTO 2100 2090 CONTINUE *-----------------------------------------------------------------------
WRITE(10,150)
IEND=0 2140 IF(IEND.EQ.L1) GOTO 2110
IBEG=IEND+1
IEND=IEND+7
IEND=MIN0(IEND,L1)
WRITE(10,160) (I,I=IBEG,IEND)
IF(MODE.EQ.3) GOTO 2120
WRITE(10,170) (X(I),I=IBEG,IEND)
GOTO 2130 2120 WRITE(10,190) (X(I),I=IBEG,IEND) 2130 GOTO 2140 2110 WRITE(10,150)
JEND=0 2160 IF(JEND.EQ.M1) GOTO 2150
JBEG=JEND+1
JEND=JEND+7
JEND=MIN0(JEND,M1)
WRITE(10,200) (J,J=JBEG,JEND)
WRITE(10,210) (Y(J),J=JBEG,JEND)
GOTO 2160 2150 CONTINUE *-----------------------------------------------------------------------
DO 2170 N=1,NGAM
NF=N
IF(.NOT.LPRINT(NF)) GOTO 2170
WRITE(10,150)
WRITE(10,110) 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 2190 CONTINUE
IBEG=IBEG+7
IEND=IBEG+6
IEND=MIN0(IEND,L1)
WRITE(10,150)
WRITE(10,120) (I,I=IBEG,IEND)
WRITE(10,130)
JFL=JFST+M1
DO 2180 JJ=JFST,M1
J=JFL-JJ
WRITE(10,140) J,(F(I,J,NF),I=IBEG,IEND) 2180 CONTINUE
IF(IEND.LT.L1) GOTO 2190 2170 CONTINUE
RETURN
END

雨烟子

这个程序好像成了一个程序，我用得是四个子程序组成，main1,main1,main2,user组成，同时有几个例子，如果是固体和流体耦合，广义扩散系数不知道如何取，它给出得例子对于动量老是取1，如果取实际（10得负几次方）老是提示被零除，不知道为什么？