p6.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  6   -----------------------------
*-----------------------------------------------------------------------
	DIMENSION T(22,22)
	EQUIVALENCE (F(1,1,4),T(1,1))
*-----------------------------------------------------------------------
	GOTO (1,2,3,4,5,6),K
*-----------------------------------------------------------------------
*	ENTRY GRID
1	XL=.5
	YL=2.
	L1=7
	M1=12
	MODE=1
	CALL SUPPLY(1)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY START
2	CONTINUE
	DO 14 I=1,4
	LSOLVE(I)=.TRUE.
14	LPRINT(I)=.TRUE.
	TITLE(1)='  VEL U  '
	TITLE(2)='  VEL V  '
	TITLE(3)='  STR FN  '
	TITLE(4)='   TEMP   '
	TITLE(11)=' PRESSURE '
	TITLE(12)='  DENSITY '
	RELAX(1)=0.5
	RELAX(2)=0.5
	RELAX(11)=0.8
	LPRINT(11)=.TRUE.
	LPRINT(12)=.TRUE.
	LAST=30
	TIN=500.
	TW=300.
	VIN=100.
	VOUT=VIN*XCV(L2)/X(L1)*TW/TIN
	DO 100 J=1,M1
	DO 100 I=1,L1
	U(I,J)=0.
	V(I,J)=VOUT
	V(1,J)=0.
	V(I,2)=0.
	T(I,J)=TW
100	CONTINUE
	V(L2,2)=VIN
	T(L2,1)=TIN
	PR=.7
	AMU=1.
	AMUP=AMU/PR
	TREF=300.
	RHOREF=1.
	RHOT=RHOREF*TREF
	RETURN
*-----------------------------------------------------------------------
*	ENTRY DENSE
3	CONTINUE
	DO 200 J=1,M1
	DO 200 I=1,L1
	RHO(I,J)=RHOT/T(I,J)
200	CONTINUE
	RETURN
*-----------------------------------------------------------------------
*	ENTRY BOUND
4	CONTINUE
	IF(ITER.NE.0) GO TO 300
	FLOWIN=RHO(L2,1)*V(L2,2)*XCV(L2)
300	FL=0.
	DO 301 I=2,L2
301	FL=FL+RHO(I,M1)*V(I,M2)*XCV(I)
	FACTOR=FLOWIN/FL
	DO 302 I=2,L2
	V(I,M1)=V(I,M2)*FACTOR
302	T(I,M1)=T(I,M2)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY OUTPUT
5	CONTINUE
	IF(ITER.NE.0) GO TO 400
	WRITE(*,401)
401	FORMAT('   ITER',9X,'SMAX',9X,'SSUM',7X,'V(4,7)',7X,
     + 'T(4,7)',5X,'F(4,7,3)')
400	WRITE(*,403) ITER,SMAX,SSUM,V(4,7),T(4,7),F(4,7,3)
403	FORMAT(I6,1P5E13.3)
	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.EQ.4) GAM(I,J)=AMUP
	IF(NF.NE.1) GAM(L1,J)=0. 
	GAM(I,M1)=0.
500	CONTINUE
	RETURN
	END