p9.for

关于simple算法

*=======================================================================
      SUBROUTINE USER(K)
*-----------------------------------------------------------------------
$INCLUDE:'SIMPLE.INC'
      DIMENSION TH(22),THU(22),THDIF(22),THCV(22),THCVS(22)
      EQUIVALENCE(X,TH),(XU,THU),(XDIF,THDIF),(XCV,THCV),
     + (XCVS,THCVS),(XL,THL)
*-----------------------------------------------------------------------
      COMMON/CNTL/LSTOP
      COMMON/SORC/SMAX,SSUM
      COMMON/COEF/FLOW,DIFF,ACOF
*-----------------------------------------------------------------------
*--------------------------   PROBLEM  9   -----------------------------
*-----------------------------------------------------------------------
	DIMENSION FLOX(22,22),FUEL(22,22),OX(22,22),H(22,22),
     + T(22,22),RSW(22,22)
	EQUIVALENCE (F(1,1,4),FLOX(1,1)),(F(1,1,5),FUEL(1,1)),
     + (F(1,1,6),OX(1,1)),(F(1,1,7),H(1,1)),(F(1,1,8),T(1,1)),
     + (F(1,1,9),RSW(1,1))
*-----------------------------------------------------------------------
	GOTO (1,2,3,4,5,6),K
*-----------------------------------------------------------------------
*	ENTRY GRID
1	XL=1.
	YL=.2
	L1=7
	M1=10
	R(1)=0.
	MODE=2
	CALL SUPPLY(1)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY START
2	CONTINUE
	DO 14 I=1,3
	LSOLVE(I)=.TRUE.
14	LPRINT(I)=.TRUE.
	LSOLVE(4)=.TRUE.
	LSOLVE(7)=.TRUE.
	LSOLVE(9)=.TRUE.
	LPRINT(5)=.TRUE.
	LPRINT(6)=.TRUE.
	LPRINT(8)=.TRUE.
	LPRINT(9)=.TRUE.
	LPRINT(11)=.TRUE.
	LPRINT(12)=.TRUE.
	TITLE(1)=' VEL U '
	TITLE(2)=' VEL V '
	TITLE(3)=' STR FN'
	TITLE(8)=' TEMP  '
	TITLE(11)='PRESSURE'
	TITLE(5)='  FUEL   '
	TITLE(6)=' OXYGEN '
	TITLE(9)='  R.VTH '
	TITLE(12)='  RHO   '
	RELAX(1)=0.5
	RELAX(2)=0.5
	RELAX(11)=0.8
	RELAX(12)=0.5
	LAST=30
	AI=4.
	CP=1000.
	TIN=300.
	RHOCON=10.
	WFU=16.
	WOX=32.
	WPR=28.
	HFU=3.E+7
	ARFU=20.
	DO 100 J=1,M1
	DO 100 I=1,L1
	U(I,J)=20.
	U(2,J)=0.
	U(I,M1)=0.
	V(I,J)=0.
	RSW(I,J)=0.
	P(I,J)=0.
	FLOX(I,J)=-.232/AI
	H(I,J)=CP*TIN
100	CONTINUE
	PR=.7
	AMU=.1
	AMUP=AMU/PR
	U(2,M2)=20.
	U(2,3)=50.
	RATIO=1./(WFU*(.232/WOX+.768/WPR))
	U(2,2)=(U(2,3)*YCVR(3)+U(2,M2)*YCVR(M2))*RATIO/(YCVR(2)*ARFU)
	H(1,2)=CP*TIN+HFU
	FLOX(1,2)=1.
	RSW(1,3)=R(3)*U(2,3)*.2
	RETURN
*-----------------------------------------------------------------------
*	ENTRY DENSE
3	CONTINUE
	DO 200 J=1,M1
	DO 200 I=1,L1
	FUEL(I,J)=AMAX1(0.,FLOX(I,J))
	OX(I,J)=AMAX1(0.,-AI*FLOX(I,J))
	T(I,J)=(H(I,J)-FUEL(I,J)*HFU)/CP
	PROD=1.-FUEL(I,J)-OX(I,J)
	WAV=1./(FUEL(I,J)/WFU+OX(I,J)/WOX+PROD/WPR)
	DEN=RHOCON*WAV/T(I,J)
	IF(ITER.EQ.0) RHO(I,J)=DEN
	RHO(I,J)=DEN*RELAX(NRHO)+RHO(I,J)*(1.-RELAX(NRHO))
200	CONTINUE
	RETURN
*-----------------------------------------------------------------------
*	ENTRY BOUND
4	CONTINUE
	IF(ITER.NE.0)GOTO 300
	FLOWIN=0.
	DO 301 J=2,M2
301	FLOWIN=FLOWIN+RHO(1,J)*U(2,J)*YCVR(J)
300	FL=0.
	DO 302 J=2,M2
	FL=FL+RHO(L1,J)*U(L2,J)*YCVR(J)
	H(L1,J)=H(L2,J)
302	FLOX(L1,J)=FLOX(L2,J)
	FACTOR=FLOWIN/FL
	DO 303 J=2,M2
303	U(L1,J)=U(L2,J)*FACTOR
	RETURN
*-----------------------------------------------------------------------
*	ENTRY OUTPUT
5	CONTINUE
	IF(ITER.NE.0) GO TO 400
	WRITE(*,401)
401	FORMAT('   ITER',6X,'SMAX',9X,'SSUM',8X,'U(4,2)',7X,
     + 'T(4,2)',7X,'OX(4,6)')
400	WRITE(*,403) ITER,SMAX,SSUM,U(4,2),T(4,2),OX(4,6)
403	FORMAT(I6,1P5E13.4)
	IF(ITER.EQ.LAST)CALL SUPPLY(2)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY GAMSOR
6	CONTINUE
	DO 500 J=1,M1
	DO 500 I=1,L1
	GAM(I,J)=AMU
	IF(NF.NE.4.AND.NF.NE.7) GO TO 501
	GAM(I,J)=AMUP
	GAM(I,M1)=0.
	GAM(1,J)=0.
501	GAM(L1,J)=0.
500	CONTINUE
	IF(NF.NE.2)GOTO 510
	DO 502 J=3,M2
	DO 502 I=2,L2
	RSWM=FY(J)*RSW(I,J)+FYM(J)*RSW(I,J-1)
	RHOM=FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1)
	CON(I,J)=RHOM*RSWM**2/RMN(J)**3
	AP(I,J)=-AMU/RMN(J)**2
502	CONTINUE
510	CONTINUE
	IF(NF.NE.9)RETURN
	DO 512 J=2,M2
	DO 512 I=2,L2
	AR=2.*AMU/YCVR(J)
	CON(I,J)=AR*RSW(I,J-1)
512	AP(I,J)=-AR
	RETURN
	END