CFD Online Discussion Forums - 2D CFD code using SIMPLE algorithm

I know many people can have this code. Do you have document to explain it? If yes, pls give a copy.

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,NT IMES(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,NT IMES(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)*R HO(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)*RH O(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,NT IMES(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