Question

develop a mathematical model using fortran 77 on pulse tube cryocooler using r 407 as a refrigerant please write the code in fortran 77 language only and verify the code generated with the help of photon 95628 from both the virgin the answer is accurate validate the reason with the available theoretical answers

Answer 1

C *** TOTAL COEFFICIENTS

AW(I,J)=DW*(AAW+AMAX1(PECLW,0.0))

AE(I,J)=DE*(AAE+AMAX1(-PECLE,0.0))

AS(I,J)=DS*(AAS+AMAX1(PECLS,0.0))

AN(I,J)=DN*(AAN+AMAX1(-PECLN,0.0))

1 CONTINUE

GO TO 2000

C COEFFICIENTS OF PRESSURE CORRECTION EQUATION

1000 DO 2 J=2,JNM

DO 2 I=2,INM

PP(I,J)=0.0

SP(I,J)=0.0

SU(I,J)=0.0

LW=NTAGW(I,J)/10

LE=NTAGE(I,J)/20

LS=NTAGS(I,J)/30

LN=NTAGN(I,J)/40

LB=1-NTAG(I,J)P1: ICD/GKJ P2: IWV

DXW=X(IN)-X(INM)

DXE=X(2)-X(1)

DYS=Y(JN)-Y(JNM)

DYN=Y(2)-Y(1)

C WEST

SUMW=FINTW(APU,I,J)

AW(I,J)=FINTW(RHO,I,J)*(R(J)*DYP(J))**2/SUMW*(1-LW)*LB

IF(NTAGW(I,J).EQ.17)THEN

JJ=J

IF(IPERIOD.EQ.1)JJ=JN-J+1

RHOW=(DXW*RHO(2,J)+DXE*RHO(INM,JJ))/(DXE+DXW)

SUMW=(DXW*APU(2,J)+DXE*APU(INM,JJ))/(DXE+DXW)

AW(I,J)=RHOW*(R(J)*DYP(J))**2/SUMW*2*LB

ENDIF

C EAST

SUME=FINTE(APU,I,J)

AE(I,J)=FINTE(RHO,I,J)*(R(J)*DYP(J))**2/SUME*(1-LE)*LB

IF(NTAGE(I,J).EQ.27)THEN

JJ=J

IF(IPERIOD.EQ.1)JJ=JN-J+1

RHOE=(DXW*RHO(2,JJ)+DXE*RHO(INM,J))/(DXE+DXW)

SUME=(DXW*APU(2,JJ)+DXE*APU(INM,J))/(DXE+DXW)

AE(I,J)=RHOE*(R(J)*DYP(J))**2/SUME*2.0*LB

ENDIF

C SOUTH

SUMS=FINTS(APV,I,J)

AS(I,J)=FINTS(RHO,I,J)*(RC(J)*DXP(I))**2/SUMS*(1-LS)*LB

IF(NTAGS(I,J).EQ.37)THEN

II=I

IF(JPERIOD.EQ.1)II=IN-I+1

RHOS=(DYS*RHO(I,2)+DYN*RHO(II,JNM))/(DYN+DYS)

SUMS=(DYS*APV(I,2)+DYN*APV(II,JNM))/(DYN+DYS)

AS(I,J)=RHOS*(RC(J)*DXP(I))**2/SUMS*2.0*LB

ENDIF

C NORTH

SUMN=FINTN(APV,I,J)

AN(I,J)=FINTN(RHO,I,J)*(RC(J+1)*DXP(I))**2/SUMN*(1-LN)*LB

IF(NTAGN(I,J).EQ.47)THEN

II=I

IF(JPERIOD.EQ.1)II=IN-I+1

RHON=(DYS*RHO(II,2)+DYN*RHO(I,JNM))/(DYN+DYS)

SUMN=(DYS*APV(II,2)+DYN*APV(I,JNM))/(DYN+DYS)P1: ICD/GKJ P2: IWV

AN(I,J)=RHON*(RC(J+1)*DXP(I))**2/SUMN*2.0*LB

ENDIF

2 CONTINUE

2000 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE SORCE(NNV)

INCLUDE ’COM2D.FOR’

C *******************************************

DIMENSION PROD(IT,JT)

N=NNV

GO TO (10,20,30,40,50,60,70),N

C *** FOR PRESSURE CORRECTION

10 DO 11 J=2,JNM

DO 11 I=2,INM

CW=FINTW(RHO,I,J)*FINTW(U,I,J)*R(J)*DYP(J)

CE=FINTE(RHO,I,J)*FINTE(U,I,J)*R(J)*DYP(J)

CS=FINTS(RHO,I,J)*FINTS(V,I,J)*RC(J)*DXP(I)

CN=FINTN(RHO,I,J)*FINTN(V,I,J)*RC(J+1)*DXP(I)

SM=CE-CW+CN-CS

IF(UNSTDY)SM=SM+(RHO(I,J)-RHOO(I,J))/DELT*VOL(I,J)

SU(I,J)=SU(I,J)-SM*(1-NTAG(I,J))

11 CONTINUE

GO TO 1000

C *** FOR U-VELOCITY

20 DO 21 J=2,JNM

DO 21 I=2,INM

DPDX=(FINTE(P,I,J)-FINTW(P,I,J))/DXP(I)

SU(I,J)=SU(I,J)-DPDX*VOL(I,J)*(1-NTAG(I,J))

21 CONTINUE

GO TO 1000

C *** FOR V-VELOCITY

30 DO 31 J=2,JNM

DO 31 I=2,INM

DPDY=(FINTN(P,I,J)-FINTS(P,I,J))/DYP(J)

SU(I,J)=SU(I,J)-DPDY*VOL(I,J)*(1-NTAG(I,J))

VISP=VIS(I,J)+VIST(I,J)

IF(AXISYMM)SP(I,J)=SP(I,J)+VISP/R(J)**2*VOL(I,J)

31 CONTINUE

GO TO 1000

C *** FOR W-VELOCITYP1: ICD/GKJ P2: IWV

40 DO 41 J=2,JNM

DO 41 I=2,INM

SU(I,J)=SU(I,J)+0.0

41 CONTINUE

GO TO 1000

C *** FOR KINETIC ENERGY

50 DO 51 J=2,JNM

DO 51 I=2,INM

IF(NTAG(I,J).EQ.1)GO TO 51

LW=NTAGW(I,J)

LE=NTAGE(I,J)

LS=NTAGS(I,J)

LN=NTAGN(I,J)

C EXCLUDE NEAR-WALL NODE

IF(LW.EQ.14.OR.LW.EQ.16)GO TO 51

IF(LE.EQ.24.OR.LE.EQ.26)GO TO 51

IF(LS.EQ.34.OR.LS.EQ.36)GO TO 51

IF(LN.EQ.44.OR.LN.EQ.46)GO TO 51

C PRODUCTION TERMS

DUDX=(FINTE(U,I,J)-FINTW(U,I,J))/DXP(I)

DVDX=(FINTE(V,I,J)-FINTW(V,I,J))/DXP(I)

DUDY=(FINTN(U,I,J)-FINTS(U,I,J))/DYP(J)

DVDY=(FINTN(V,I,J)-FINTS(V,I,J))/DYP(J)

TERM=2.0*(DUDX**2+DVDY**2)+(DUDY+DVDX)**2

IF(AXISYMM)TERM=TERM+2*(V(I,J)/R(J))**2

PROD(I,J)=VIST(I,J)*TERM

ENP=AMAX1(E(I,J),0.0)

SU(I,J)=PROD(I,J)*VOL(I,J) +SU(I,J)

SP(I,J)=RHO(I,J)**2*CMU*ENP/(VIST(I,J)+SMALL)*VOL(I,J)+SP(I,J)

51 CONTINUE

GO TO 1000

C *** FOR DISSIPATION

60 DO 61 J=2,JNM

DO 61 I=2,INM

IF(NTAG(I,J).EQ.1)GO TO 61

RHOP=RHO(I,J)

VOLP=VOL(I,J)

DPEP=ABS(D(I,J)/(E(I,J)+SMALL))

SU(I,J)=CD1*RHOP*DPEP*PROD(I,J)*VOLP +SU(I,J)

SP(I,J)=CD2*RHOP*DPEP*VOLP+SP(I,J)

61 CONTINUE

GO TO 1000P1: ICD/GKJ P2: IWV

C *** FOR TEMPERATURE

70 DO 71 J=2,JNM

DO 71 I=2,INM

SU(I,J)=SU(I,J)+0.0

71 CONTINUE

1000 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE APCOF(NN)

INCLUDE ’COM2D.FOR’

C *******************************************

N=NN

RPINV=1./RP(N)

DO 1 J=2,JNM

DO 1 I=2,INM

SUM=AW(I,J)+AE(I,J)+AS(I,J)+AN(I,J)

IF(N.EQ.1)AP1(I,J)=(SUM+SP(I,J))*RPINV

IF(N.GT.1)AP(I,J)=(SUM+SP(I,J))*RPINV

IF(N.EQ.2)APU(I,J)=AP(I,J)

IF(N.EQ.3)APV(I,J)=AP(I,J)

1 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE PROPS

INCLUDE ’COM2D.FOR’

C *******************************************

DO 1 J=1,JN

DO 1 I=1,IN

SPH(I,J)=SPHEAT

RHO(I,J)=DENSIT

VIS(I,J)=VISCOS

IF(TURBUL)THEN

VISO=VIST(I,J)

VIST(I,J)=(CMU*RHO(I,J)*E(I,J)**2/(D(I,J)+SMALL)*RP(8)

1 +(1.-RP(8))*VISO)

IF(VIST(I,J).LE.0.0)VIST(I,J)=SMALL

ENDIF

1 CONTINUE

IF(BSOR(8))CALL ADSORB(8)

RETURNP1: ICD/GKJ P2: IWV

END

C *******************************************

SUBROUTINE UNST(NN)

INCLUDE ’COM2D.FOR’

C *******************************************

N=NN

DO 1 J=2,JNM

DO 1 I=2,INM

SUU=SU(I,J)

TERM=RHO(I,J)*VOL(I,J)/DELT*(1-NTAG(I,J))

IF(UNSTDY)THEN

TERM=TERM*RHOO(I,J)/RHO(I,J)

IF(N.EQ.2)SU(I,J)=TERM*UO(I,J)+SUU

IF(N.EQ.3)SU(I,J)=TERM*VO(I,J)+SUU

IF(N.EQ.4)SU(I,J)=TERM*WO(I,J)+SUU

IF(N.EQ.5)SU(I,J)=TERM*EO(I,J)+SUU

IF(N.EQ.6)SU(I,J)=TERM*DO(I,J)+SUU

IF(N.EQ.7)SU(I,J)=TERM*TO(I,J)+SUU

ELSE IF(FTRAN)THEN

IF(N.EQ.2)SU(I,J)=TERM*U(I,J)+SUU

IF(N.EQ.3)SU(I,J)=TERM*V(I,J)+SUU

IF(N.EQ.4)SU(I,J)=TERM*W(I,J)+SUU

IF(N.EQ.5)SU(I,J)=TERM*E(I,J)+SUU

IF(N.EQ.6)SU(I,J)=TERM*D(I,J)+SUU

IF(N.EQ.7)SU(I,J)=TERM*T(I,J)+SUU

ENDIF

SP(I,J)=TERM +SP(I,J)

1 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE UPDATE

INCLUDE ’COM2D.FOR’

C *******************************************

DO 1 J=1,JN

DO 1 I=1,IN

RHOO(I,J)=RHO(I,J)

PO(I,J)=P(I,J)

UO(I,J)=U(I,J)

VO(I,J)=V(I,J)

WO(I,J)=W(I,J)

EO(I,J)=E(I,J)P1: ICD/GKJ P2: IWV

DO(I,J)=D(I,J)

TO(I,J)=T(I,J)

1 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE INFLUX

INCLUDE ’COM2D.FOR’

C *******************************************

DO 1 N=1,7

1 RNORM(N)=0.0

DO 2 J=2,JNM

DO 2 I=2,INM

IF(NTAGW(I,J).EQ.11) THEN

CW=ABS(RHO(I-1,J)*U(I-1,J)*DYP(J)*R(J))

VT=SQRT(U(I-1,J)**2+V(I-1,J)**2+W(I-1,J)**2)

RNORM(1)=RNORM(1)+CW

RNORM(2)=RNORM(2)+CW*VT

RNORM(3)=RNORM(2)

RNORM(4)=RNORM(2)

RNORM(5)=RNORM(5)+CW*ABS(E(I-1,J))

RNORM(6)=RNORM(6)+CW*ABS(D(I-1,J))

RNORM(7)=RNORM(7)+CW*ABS(T(I-1,J))

ELSE IF(NTAGE(I,J).EQ.21) THEN

CE=ABS(RHO(I+1,J)*U(I+1,J)*DYP(J)*R(J))

VT=SQRT(U(I+1,J)**2+V(I+1,J)**2+W(I+1,J)**2)

RNORM(1)=RNORM(1)+CE

RNORM(2)=RNORM(2)+CE*VT

RNORM(3)=RNORM(2)

RNORM(4)=RNORM(2)

RNORM(5)=RNORM(5)+CE*ABS(E(I+1,J))

RNORM(6)=RNORM(6)+CE*ABS(D(I+1,J))

RNORM(7)=RNORM(7)+CE*ABS(T(I+1,J))

ELSE IF(NTAGS(I,J).EQ.31) THEN

CS=ABS(RHO(I,J-1)*V(I,J-1)*DXP(I)*RC(J))

VT=SQRT(U(I,J-1)**2+V(I,J-1)**2+W(I,J-1)**2)

RNORM(1)=RNORM(1)+CS

RNORM(2)=RNORM(2)+CS*VT

RNORM(3)=RNORM(2)

RNORM(4)=RNORM(2)

RNORM(5)=RNORM(5)+CS*ABS(E(I,J-1))

RNORM(6)=RNORM(6)+CS*ABS(D(I,J-1))P1: ICD/GKJ P2: IWV

RNORM(7)=RNORM(7)+CS*ABS(T(I,J-1))

ELSE IF(NTAGN(I,J).EQ.41) THEN

CN=ABS(RHO(I,J+1)*V(I,J+1)*DXP(I)*RC(J+1))

VT=SQRT(U(I,J+1)**2+V(I,J+1)**2+W(I,J+1)**2)

RNORM(1)=RNORM(1)+CN

RNORM(2)=RNORM(2)+CN*VT

RNORM(3)=RNORM(2)

RNORM(4)=RNORM(2)

RNORM(5)=RNORM(5)+CN*ABS(E(I,J+1))

RNORM(6)=RNORM(6)+CN*ABS(D(I,J+1))

RNORM(7)=RNORM(7)+CN*ABS(T(I,J+1))

ENDIF

2 CONTINUE

DO 3 N=1,7

TERM=ABS(RNORM(N))

3 IF(TERM.LT.10.*SMALL)RNORM(N)=1.0

WRITE(6,*)’ RNORM VALUES’

WRITE(6,*)(RNORM(N),N=1,7)

RETURN

END

C *******************************************

SUBROUTINE MASBAL

INCLUDE ’COM2D.FOR’

C *******************************************

SUMFW=0.0

SUMFE=0.0

SUMFS=0.0

SUMFN=0.0

DO 2 J=2,JNM

DO 2 I=2,INM

IF(NTAGW(I,J).EQ.13.OR.NTAGW(I,J).EQ.15) THEN

CW=RHO(I-1,J)*U(I-1,J)*DYP(J)*R(J)

SUMFW=SUMFW+CW

ELSE IF(NTAGE(I,J).EQ.23.OR.NTAGE(I,J).EQ.25) THEN

CE=RHO(I+1,J)*U(I+1,J)*DYP(J)*R(J)

SUMFE=SUMFE+CE

ELSE IF(NTAGS(I,J).EQ.33.OR.NTAGS(I,J).EQ.35) THEN

CS=RHO(I,J-1)*V(I,J-1)*DXP(I)*RC(J)

SUMFS=SUMFS+CS

ELSE IF(NTAGN(I,J).EQ.43.OR.NTAGN(I,J).EQ.45) THEN

CN=RHO(I,J+1)*V(I,J+1)*DXP(I)*RC(J+1)

SUMFN=SUMFN+CNP1: ICD/GKJ P2: IWV

ENDIF

2 CONTINUE

SUMF=ABS(SUMFW)+ABS(SUMFE)+ABS(SUMFS)+ABS(SUMFN)

FACTOR=RNORM(1)/(SUMF+SMALL)

WRITE(6,8787)FACTOR

WRITE(*,8787)FACTOR

C APPLY MASS CONSERVATION AT EXIT

IF(CONMAS)THEN

DO 3 J=2,JNM

DO 3 I=2,INM

IF(NTAGW(I,J).EQ.13.OR.NTAGW(I,J).EQ.15) THEN

U(I-1,J)=U(I-1,J)*FACTOR

V(I-1,J)=V(I-1,J)*FACTOR

W(I-1,J)=W(I-1,J)*FACTOR

ELSE IF(NTAGE(I,J).EQ.23.OR.NTAGE(I,J).EQ.25) THEN

U(I+1,J)=U(I+1,J)*FACTOR

V(I+1,J)=V(I+1,J)*FACTOR

W(I+1,J)=W(I+1,J)*FACTOR

ELSE IF(NTAGS(I,J).EQ.33.OR.NTAGS(I,J).EQ.35) THEN

V(I,J-1)=V(I,J-1)*FACTOR

U(I,J-1)=U(I,J-1)*FACTOR

W(I,J-1)=W(I,J-1)*FACTOR

ELSE IF(NTAGN(I,J).EQ.43.OR.NTAGN(I,J).EQ.45) THEN

V(I,J+1)=V(I,J+1)*FACTOR

U(I,J+1)=U(I,J+1)*FACTOR

W(I,J+1)=W(I,J+1)*FACTOR

ENDIF

3 CONTINUE

ENDIF

8787 FORMAT(50X,F10.4,F10.4)

RETURN

END

C *******************************************

SUBROUTINE PVCOR

INCLUDE ’COM2D.FOR’

C *******************************************

C **** APPLY SMOOTHING PRESSURE CORRECTION

DO 4 J=2,JNM

DO 4 I=2,INM

PMX=(DXMI(I)*P(I+1,J)+DXMI(I+1)*P(I-1,J))/(DXMI(I)+DXMI(I+1))

PMY=(DYMI(J)*P(I,J+1)+DYMI(J+1)*P(I,J-1))/(DYMI(J)+DYMI(J+1))

PSM(I,J)=(P(I,J)-(PMX+PMY)/2.0)*GAMMAP1: ICD/GKJ P2: IWV

PP(I,J)=(PP(I,J)-PSM(I,J))*(1-NTAG(I,J))

4 CONTINUE

C *** APPLY MASS-CONSERVING PRESSURE CORRECTION

PREF=0

RSP=0.0

DO 6 J=2,JNM

DO 6 I=2,INM

P(I,J)=P(I,J)+(PP(I,J)-PREF)*RP(9)*(1-NTAG(I,J))

IF(ABS(PP(I,J)).GT.RSP)RSP=ABS(PP(I,J))

6 CONTINUE

FDIF(1)=RSP

CALL BOUNDP

C *** CORRECT VELOCITIES

RSU=0.0

RSV=0.0

DO 1 J=2,JNM

DO 1 I=2,INM

PSMW=FINTW(PP,I,J)

PSME=FINTE(PP,I,J)

PSMS=FINTS(PP,I,J)

PSMN=FINTN(PP,I,J)

C CORRECT U-VELOCITY

IF(SLVE(2))THEN

DPDX=(PSME-PSMW)/DXP(I)

UDASH=-DPDX*VOL(I,J)/APU(I,J)*(1-NTAG(I,J))

IF(ABS(UDASH).GT.RSU)RSU=ABS(UDASH)

U(I,J)=U(I,J)+UDASH

ENDIF

C CORRECT V-VELOCITY

IF(SLVE(3))THEN

DPDY=(PSMN-PSMS)/DYP(J)

VDASH=-DPDY*VOL(I,J)/APV(I,J)*(1-NTAG(I,J))

IF(ABS(VDASH).GT.RSV)RSV=ABS(VDASH)

V(I,J)=V(I,J)+VDASH

ENDIF

1 CONTINUE

FDIF(2)=RSU

FDIF(3)=RSV

C CHECK MASS RESIDUAL

SUM=0.0

DO 9 J=2,JNM

DO 9 I=2,INMP1: ICD/GKJ P2: IWV

TERM=AE(I,J)*PP(I+1,J)+AW(I,J)*PP(I-1,J)

1 +AN(I,J)*PP(I,J+1)+AS(I,J)*PP(I,J-1)-AP1(I,J)*PP(I,J)

IF(TERM.GT.GREAT*0.01)TERM=0.0

SUM=SUM+TERM**2*(1-NTAG(I,J))

9 CONTINUE

RSDU(1)=SQRT(SUM)/RNORM(1)

RETURN

END

C *******************************************

SUBROUTINE BOUNDP

INCLUDE ’COM2D.FOR’

C *******************************************

DO 2 J=2,JNM

DO 2 I=2,INM

IF (NTAG(I,J).EQ.1) GO TO 2

LW=NTAGW(I,J)/10

LE=NTAGE(I,J)/20

LS=NTAGS(I,J)/30

LN=NTAGN(I,J)/40

DXW=X(IN)-X(INM)

DXE=X(2)-X(1)

DYS=Y(JN)-Y(JNM)

DYN=Y(2)-Y(1)

C EAST-WEST PERIODICITY

IF(NTAGW(I,J).EQ.17.OR.NTAGE(I,J).EQ.27)THEN

JJ=J

IF(IPERIOD.EQ.1)JJ=JN-J+1

ENDIF

IF(LW.EQ.1) THEN

RATIO=(X(I)-XC(I))/DXMI(I+1)

P(I-1,J)=P(I,J)-RATIO*(P(I+1,J)-P(I,J))

PP(I-1,J)=PP(I,J)-RATIO*(PP(I+1,J)-PP(I,J))

IF(NTAGW(I,J).EQ.17)THEN

PMEAN=(DXW*P(2,J)+DXE*P(INM,JJ))/(DXE+DXW)

P(1,J)=PMEAN+DP1/2.0

PPMEAN=(DXW*PP(2,J)+DXE*PP(INM,JJ))/(DXE+DXW)

PP(1,J)=PPMEAN

ENDIF

ENDIF

IF(LE.EQ.1) THEN

RATIO=(XC(I+1)-X(I))/DXMI(I)

P(I+1,J)=P(I,J)+RATIO*(P(I,J)-P(I-1,J))P1: ICD/GKJ P2: IWV

PP(I+1,J)=PP(I,J)+RATIO*(PP(I,J)-PP(I-1,J))

IF(NTAGE(I,J).EQ.27)THEN

PMEAN=(DXW*P(2,JJ)+DXE*P(INM,J))/(DXE+DXW)

P(IN,J)=PMEAN-DP1/2.0

PPMEAN=(DXW*PP(2,JJ)+DXE*PP(INM,J))/(DXE+DXW)

ENDIF

ENDIF

C NORTH-SOUTH PERIODICITY

IF(NTAGS(I,J).EQ.37.OR.NTAGN(I,J).EQ.47)THEN

II=I

IF(JPERIOD.EQ.1)II=IN-I+1

ENDIF

IF(LS.EQ.1) THEN

RATIO=(Y(J)-YC(J))/DYMI(J+1)

P(I,J-1)=P(I,J)-RATIO*(P(I,J+1)-P(I,J))

PP(I,J-1)=PP(I,J)-RATIO*(PP(I,J+1)-PP(I,J))

IF(NTAGS(I,J).EQ.37)THEN

P(I,1)=(DYS*P(I,2)+DYN*P(II,JNM))/(DYN+DYS)

P(I,1)=P(I,1)+DP2/2.0

PP(I,1)=(DYS*PP(I,2)+DYN*PP(II,JNM))/(DYN+DYS)

ENDIF

ENDIF

IF(LN.EQ.1) THEN

RATIO=(YC(J+1)-Y(J))/DYMI(J)

P(I,J+1)=P(I,J)+RATIO*(P(I,J)-P(I,J-1))

PP(I,J+1)=PP(I,J)+RATIO*(PP(I,J)-PP(I,J-1))

IF(NTAGN(I,J).EQ.47)THEN

P(I,JN)=(DYS*P(II,2)+DYN*P(I,JNM))/(DYN+DYS)

P(I,JN)=P(I,1)-DP2/2.0

PP(I,JN)=(DYS*PP(II,2)+DYN*PP(I,JNM))/(DYN+DYS)

ENDIF

ENDIF

2 CONTINUE

IF(BSOR(9))CALL ADSORB(9)

RETURN

END

C *******************************************

SUBROUTINE INDATA

INCLUDE ’COM2D.FOR’

C *******************************************

WRITE(6,*)’*** THIS IS COLLOCATED GRID PROGRAM ***’

WRITE(6,*)’*************************************’P1: ICD/GKJ P2: IWV

WRITE(6,*)’GRID INFORMATION’

WRITE(6,*)’ IN = ’,IN,’ JN = ’,JN

WRITE(6,*)’ X - COORDINATES ’

CALL PR1D(X,1,IN)

WRITE(6,*)’ Y - COORDINATES ’

CALL PR1D(Y,1,JN)

WRITE(6,*)’ XC - COORDINATES ’

CALL PR1D(XC,1,IN)

WRITE(6,*)’ YC - COORDINATES ’

CALL PR1D(YC,1,JN)

WRITE(6,*)’ DXMI’

CALL PR1D(DXMI,1,IN)

WRITE(6,*)’ DYMI’

CALL PR1D(DYMI,1,JN)

WRITE(6,*)’ PRESSURE REFERENCE POINT IPREF = ’,IPREF,

1 ’ JPREF = ’,JPREF

WRITE(6,*)’RELAXATION PARAMETERS ARE’

WRITE(6,*)’ RP(1) = ’,RP(1),’ RP(2) = ’,RP(2),’ RP(3) = ’,RP(3),

1 ’ RP(4) = ’,RP(4),’ RP(5) = ’,RP(5), ’ RP(6) = ’,RP(6),

1 ’ RP(7) = ’,RP(7),’ RP(8) = ’,RP(8), ’ RP(9) = ’,RP(9)

WRITE(6,*)’FLUID VISCOSITY = ’,VISCOS

WRITE(6,*)’FLUID DENSITY = ’,DENSIT

WRITE(6,*)’FLUID PRANDTL NUMBERS ARE’

WRITE(6,*)’ PR(1) = ’,PR(1),’ PR(2) = ’,PR(2),’ PR(3) = ’,PR(3),

1 ’ PR(4) = ’,PR(4),’ PR(5) = ’,PR(5),’ PR(6) = ’,PR(6),

1 ’ PR(7) = ’,PR(7)

IF(TURBUL) THEN

WRITE(6,*)’TURBULENT PRANDTL NUMBERS ARE’

WRITE(6,*)

1 ’ PRT(1) = ’,PRT(1),’ PRT(2) = ’,PRT(2),’ PRT(3) = ’,PRT(3),

1 ’ PRT(4) = ’,PRT(4),’ PRT(5) = ’,PRT(5), ’ PRT(6) = ’,PRT(6),

1 ’ PRT(7) = ’,PRT(7)

ENDIF

IF(STEADY)WRITE(6,*)’ STEADY FLOW CALCULATIONS’

IF(UNSTDY)WRITE(6,*)’ UNSTEADY FLOW CALCULATIONS’

IF(FTRAN)WRITE(6,*)’ FALSE TRANSIENT DELT = ’,DELT

IF(CONMAS)WRITE(6,*)’ MASS BALANCE IS IMPOSED’

IF(UPWIND)WRITE(6,*)’ CONVECTION SCHEME = UPWIND ’

IF(HYBRID)WRITE(6,*)’ CONVECTION SCHEME = HYBRID ’

IF(POWER)WRITE(6,*)’ CONVECTION SCHEME = POWER LAW’

WRITE(6,*)’ THE FOLLOWING EQUATIONS ARE SOLVED’

IF(SLVE(1))WRITE(6,*)’ PRESSURE CORRECTION EQUN.’P1: ICD/GKJ P2: IWV

IF(SLVE(2))WRITE(6,*)’ U-VELOCITY ’

IF(SLVE(3))WRITE(6,*)’ V-VELOCITY ’

IF(SLVE(4))WRITE(6,*)’ W-VELOCITY ’

IF(SLVE(5))WRITE(6,*)’ T. KINETIC ENERGY’

IF(SLVE(6))WRITE(6,*)’ DISSIPATION ’

IF(SLVE(7))WRITE(6,*)’ TEMPERATURE’

WRITE(6,*)’******************************************’

RETURN

END

C *******************************************

SUBROUTINE SOLVE(F,RPP,RSUM)

INCLUDE ’COM2D.FOR’

C *******************************************

DIMENSION F(IT,JT)

DIMENSION SA(MXGR),SB(MXGR),SS(MXGR),PSI(MXGR)

C *** CALCULATION OF RESIDUALS

RS=0.0

DO 10 J=2,JNM

DO 10 I=2,INM

TERM=AW(I,J)*F(I-1,J)+AE(I,J)*F(I+1,J)

1 +AS(I,J)*F(I,J-1)+AN(I,J)*F(I,J+1)

TERM=TERM+SU(I,J)-F(I,J)*AP(I,J)*RPP

FACTOR=1.0

IF(SP(I,J).GT.GREAT*1.0E-10)FACTOR=0.0

TERM=TERM*FACTOR

10 RS=RS+TERM*TERM

RSUM=SQRT(RS)

C*** J-DIRECTION SWEEP

DO 51 J=2,JNM

DO 52 I=2,INM

SOR=SU(I,J)

DEN=1.0/(AP(I,J)+SMALL)

SOR=SOR+(1.-RPP)/(DEN+SMALL)*F(I,J)

SA(I)=AE(I,J)*DEN

SB(I)=AW(I,J)*DEN

SS(I)=(AS(I,J)*F(I,J-1)+AN(I,J)*F(I,J+1)+SOR)*DEN

52 CONTINUE

PSI1=F(1,J)

PSIN=F(IN,J)

CALL TDMA(2,INM,PSI1,PSIN,SA,SB,SS,PSI)

DO 53 I=2,INM

LP=NTAG(I,J)P1: ICD/GKJ P2: IWV

53 F(I,J)=PSI(I)*(1-LP)+LP*F(I,J)

51 CONTINUE

C*** I-DIRECTION SWEEP

DO 54 I=2,INM

DO 55 J=2,JNM

SOR=SU(I,J)

DEN=1.0/(AP(I,J)+SMALL)

SOR=SOR+(1.-RPP)/DEN*F(I,J)

SA(J)=AN(I,J)*DEN

SB(J)=AS(I,J)*DEN

SS(J)=(AW(I,J)*F(I-1,J)+AE(I,J)*F(I+1,J)+SOR)*DEN

55 CONTINUE

PSI1=F(I,1)

PSIN=F(I,JN)

CALL TDMA(2,JNM,PSI1,PSIN,SA,SB,SS,PSI)

DO 56 J=2,JNM

LP=NTAG(I,J)

56 F(I,J)=PSI(J)*(1-LP)+LP*F(I,J)

54 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE SOLP

INCLUDE ’COM2D.FOR’

C *******************************************

DIMENSION SA(MXGR),SB(MXGR),SS(MXGR),PSI(MXGR)

DO 100 L=1,NSWEEP(1)

C*** J-DIRECTION SWEEP

DO 51 J=2,JNM

DO 52 I=2,INM

SOR=SU(I,J)

DEN=1.0/(AP1(I,J)+SMALL)

SOR=SOR+(1.-RP(1))/DEN*PP(I,J)

SA(I)=AE(I,J)*DEN

SB(I)=AW(I,J)*DEN

SS(I)=(AS(I,J)*PP(I,J-1)+AN(I,J)*PP(I,J+1)+SOR)*DEN

52 CONTINUE

PSI1=PP(1,J)

PSIN=PP(IN,J)

CALL TDMA(2,INM,PSI1,PSIN,SA,SB,SS,PSI)

DO 53 I=2,INM

LP=NTAG(I,J)P1: ICD/GKJ P2: IWV

53 PP(I,J)=PSI(I)*(1-LP)+LP*PP(I,J)

51 CONTINUE

C*** I-DIRECTION SWEEP

DO 54 I=2,INM

DO 55 J=2,JNM

SOR=SU(I,J)

DEN=1.0/(AP1(I,J)+SMALL)

SOR=SOR+(1.-RP(1))/DEN*PP(I,J)

SA(J)=AN(I,J)*DEN

SB(J)=AS(I,J)*DEN

SS(J)=(AW(I,J)*PP(I-1,J)+AE(I,J)*PP(I+1,J)+SOR)*DEN

55 CONTINUE

PSI1=PP(I,1)

PSIN=PP(I,JN)

CALL TDMA(2,JNM,PSI1,PSIN,SA,SB,SS,PSI)

DO 56 J=2,JNM

LP=NTAG(I,J)

56 PP(I,J)=PSI(J)*(1-LP)+LP*PP(I,J)

54 CONTINUE

IF(NPERIOD.EQ.1)CALL BOUNDP

100 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE EQN

INCLUDE ’COM2D.FOR’

C *******************************************

MWRITE=NITER+MFREQ

IF(NITER.EQ.0)NITER=1

NADD=MXIT

5555 NBEGIN=NITER

MXIT=NITER+NADD

DO 2000 NTIME=1,MXSTEP

TTIME=STIME+NTIME*DELT

DO 1000 NITER=NBEGIN,MXIT

C **** U-VELOCITY

IF(SLVE(2))THEN

CALL COEF(2,PR(2),PRT(2))

CALL SORCE(2)

IF(UNSTDY.OR.FTRAN)CALL UNST(2)

CALL BOUND(2)

IF(BSOR(2))CALL ADSORB(2)P1: ICD/GKJ P2: IWV

CALL APCOF(2)

CALL SOLVE(U,RP(2),RSU)

RSDU(2)=RSU/(RNORM(2)+SMALL)

CALL BOUND(2)

ENDIF

C **** V-VELOCITY

IF(SLVE(3))THEN

CALL COEF(3,PR(3),PRT(3))

CALL SORCE(3)

IF(UNSTDY.OR.FTRAN)CALL UNST(3)

CALL BOUND(3)

IF(BSOR(3))CALL ADSORB(3)

CALL APCOF(3)

CALL SOLVE(V,RP(3),RSU)

RSDU(3)=RSU/(RNORM(3)+SMALL)

CALL BOUND(3)

ENDIF

C **** PRESSURE CORRECION

IF(SLVE(1))THEN

CALL MASBAL

CALL COEF(1,PR(1),PRT(1))

CALL SORCE(1)

IF(BSOR(1))CALL ADSORB(1)

CALL APCOF(1)

CALL SOLP

CALL PVCOR

ENDIF

C **** W-VELOCITY

IF(SLVE(4))THEN

CALL COEF(4,PR(4),PRT(4))

CALL SORCE(4)

IF(UNSTDY.OR.FTRAN)CALL UNST(4)

CALL BOUND(4)

IF(BSOR(4))CALL ADSORB(4)

CALL APCOF(4)

CALL SOLVE(W,RP(4),RSU)

RSDU(4)=RSU/(RNORM(4)+SMALL)

CALL BOUND(4)

ENDIF

C **** KINETIC ENERGY

IF(TURBUL)THEN

IF(SLVE(5))THENP1: ICD/GKJ P2: IWV

ALL COEF(5,PR(5),PRT(5))

CALL SORCE(5)

IF(UNSTDY.OR.FTRAN)CALL UNST(5)

CALL BOUND(5)

IF(BSOR(5))CALL ADSORB(5)

CALL APCOF(5)

CALL SOLVE(E,RP(5),RSU)

RSDU(5)=RSU/(RNORM(5)+SMALL)

CALL BOUND(5)

C **** DISSIPATION

CALL COEF(6,PR(6),PRT(6))

CALL SORCE(6)

IF(UNSTDY.OR.FTRAN)CALL UNST(6)

CALL BOUND(6)

IF(BSOR(6))CALL ADSORB(6)

CALL APCOF(6)

CALL SOLVE(D,RP(6),RSU)

RSDU(6)=RSU/(RNORM(6)+SMALL)

CALL BOUND(6)

ENDIF

ENDIF

C **** TEMPERATURE

IF(SLVE(7))THEN

CALL COEF(7,PR(7),PRT(7))

CALL SORCE(7)

IF(UNSTDY.OR.FTRAN)CALL UNST(7)

CALL BOUND(7)

IF(BSOR(7))CALL ADSORB(7)

CALL APCOF(7)

CALL SOLVE(T,RP(7),RSU)

RSDU(7)=RSU/(RNORM(7)+SMALL)

CALL BOUND(7)

ENDIF

C **** SPECIES AND ENTHALPY

CALL OMEGA

C **** PROPERTIES

CALL PROPS

C **** CHECK MAX RESIDUALS

RSTOP=AMAX1(RSDU(1),RSDU(2),RSDU(3),RSDU(4),RSDU(5)

1 ,RSDU(6),RSDU(7))

C STORE RESIDUALS FOR PLOTTING

C RESIU(NITER)=RSDU(2)P1: ICD/GKJ P2: IWV

RESIV(NITER)=RSDU(3)

C RESIM(NITER)=RSDU(1)

IF(STEADY)WRITE(6,1919)NITER,(RSDU(N),N=1,7)

IF(STEADY)WRITE(6,1919)NITER,(FDIF(N),N=1,4)

IF(STEADY)WRITE(*,1919)NITER,(FDIF(N),N=1,4)

IF(STEADY)WRITE(*,1919)NITER,(RSDU(N),N=1,7)

1919 FORMAT(1X,I6,7(E10.3))

IF(RSTOP.LT.CC) GO TO 1100

C INTERMEDIATE WRITE-OUT

IF(MWRITE.EQ.NITER)THEN

MWRITE=NITER+MFREQ

CALL OPT

WRITE(*,*)’OUTPUT IS WRITTEN AT NITER = ’,NITER

WRITE(6,*)’OUTPUT IS WRITTEN AT NITER = ’,NITER

ENDIF

1000 CONTINUE

1100 IF(STEADY)RETURN

CALL UPDATE

WRITE(6,*)’NTIME = ’,NTIME,’ TTIME = ’,TTIME

WRITE(6,*)(RSDU(N),N=1,7)

WRITE(6,*)(FDIF(N),N=1,7)

2000 CONTINUE

RETURN

END

C *******************************************

SUBROUTINE TDMA(IB,IL,Y1,YN,BA,BB,BS,YY)

INCLUDE ’COM2D.FOR’

C *******************************************