p10.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  10   -----------------------------
*-----------------------------------------------------------------------
	DIMENSION T(22,22),XPLUS(22),DIS(22,22),AKE(22,22),
     + AMUT(22,22),GEN(22,22)
	EQUIVALENCE (F(1,1,4),T(1,1)),(F(1,1,5),AKE(1,1)),
     + (F(1,1,6),DIS(1,1)),(F(1,1,7),AMUT(1,1)),(F(1,1,8),GEN(1,1))
*-----------------------------------------------------------------------
	GOTO (1,2,3,4,5,6),K
*-----------------------------------------------------------------------
*	ENTRY GRID
1	XL=1.
	YL=4.
	L1=7
	M1=9
	MODE=1
	CALL SUPPLY(1)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY START
2	CONTINUE
	DO 14 I=1,6
	LSOLVE(I)=.TRUE.
14	LPRINT(I)=.TRUE.
	LPRINT(11)=.TRUE.
	LPRINT(7)=.TRUE.
	TITLE(1)='  VEL  U  '
	TITLE(2)='  VEL  V  '
	TITLE(3)='  STR FN  '
	TITLE(4)='   TEMP   '
	TITLE(11)=' PRESSURE '
	TITLE(5)=' K.ENERGY '
	TITLE(6)=' DISIPAT. '
	TITLE(7)='  T.VISC  '
	TITLE(12)='  DENSITY '
	RELAX(1)=0.5
	RELAX(2)=0.5
	RELAX(11)=0.8
	RELAX(5)=0.4
	RELAX(6)=0.4
	RELAX(13)=0.5
	LAST=35
	DO 100 J=1,M1
	DO 100 I=1,L1
	U(I,J)=0.
	V(I,J)=10.
	V(1,J)=0.
	V(I,2)=10.
	IF(I.GT.4) V(I,2)=100.
	T(I,J)=100.
	T(1,J)=0.
	IF(I.GT.4) T(I,1)=400.
	AKE(I,J)=.005*V(I,2)**2
	DIS(I,J)=.1*AKE(I,J)**2
100	CONTINUE
	AMU=1.E-6
	CMU=.09
	C1=1.44
	C2=1.92
	PRT=.9
	PRK=1.
	PRD=1.3
	PR=.7
	PRPRT=PR/PRT
	PFN=9.*(PRPRT-1.)/PRPRT**.25
	CMU4=CMU**.25
	RETURN
*-----------------------------------------------------------------------
*	ENTRY DENSE
3	CONTINUE
	RETURN
*-----------------------------------------------------------------------
*	ENTRY BOUND
4	CONTINUE
	IF(ITER.NE.0)GOTO 300
	FLOWIN=0.
	DO 310 I=2,L2
310	FLOWIN=FLOWIN+RHO(I,1)*V(I,2)*XCV(I)
300	FL=0.
	AFL=0.
	VMIN=0.
	DO 301 I=2,L2
	IF(V(I,M2).LT.0.) VMIN=AMAX1(VMIN,-V(I,M2))
	AFL=AFL+RHO(I,M1)*XCV(I)
301	FL=FL+RHO(I,M1)*V(I,M2)*XCV(I)
	FACTOR=FLOWIN/(FL+AFL*VMIN)
	DO 302 I=2,L2
302	V(I,M1)=(V(I,M2)+VMIN)*FACTOR
	DO 303 J=2,M2
	AKE(L1,J)=AKE(L2,J)
303	DIS(L1,J)=DIS(L2,J)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY OUTPUT
5	CONTINUE
	IF(ITER.NE.0) GO TO 400
	WRITE(*,401)
401	FORMAT('   ITER',6X,'SMAX',6X,'SSUM',5X,'V(6,6)',4X,
     + 'T(5,6)',4X,'KE(5,6)')
400	WRITE(*,403) ITER,SMAX,SSUM,V(6,6),T(5,6),AKE(5,6)
403	FORMAT(I6,1P5E11.3)
	IF(ITER.EQ.LAST)CALL SUPPLY(2)
	RETURN
*-----------------------------------------------------------------------
*	ENTRY GAMSOR
6	CONTINUE
	IF(NF.EQ.3)RETURN
	IF(NF.NE.1)GOTO 510
	REL=1.-RELAX(NGAM)
	DO 500 J=1,M1
	DO 500 I=1,L1
	AMT=CMU*RHO(I,J)*AKE(I,J)**2/(DIS(I,J)+1.E-30)
	IF(ITER.EQ.0)AMUT(I,J)=AMT
500	AMUT(I,J)=RELAX(NGAM)*AMT+REL*AMUT(I,J)
510	FACTOR=1.
	IF(NF.EQ.4)FACTOR=1./PRT
	IF(NF.EQ.5)FACTOR=1./PRK
	IF(NF.EQ.6)FACTOR=1./PRD
	DO 520 J=1,M1
	DO 520 I=1,L1
	GAM(I,J)=AMUT(I,J)*FACTOR
	IF(NF.NE.1)GAM(L1,J)=0.
520	GAM(I,M1)=0.
	DO 530 J=2,M2
	GOTO (531,532,531,533,531,531),NF
531	GAM(1,J)=0.
	GOTO 530
532	XPLUS(J)=RHO(2,J)*SQRT(AKE(2,J))*CMU4*XDIF(2)/AMU
	GAM(1,J)=AMU
	IF(XPLUS(J).GT.11.5)GAM(1,J)=AMU*XPLUS(2)/(ALOG(9.*XPLUS(J))*2.5)
	GOTO 530
533	GAM(1,J)=AMU/PR
	IF(XPLUS(J).GT.11.5)GAM(1,J)=AMU/PRT*XPLUS(J)/(2.5*
     + ALOG(9.*XPLUS(J))+PFN)
530	CONTINUE
	IF(NF.NE.1)GOTO 509
	DO 590 J=2,M2
	DO 590 I=3,L2
	CON(I,J)=(GAM(I,J)*(U(I+1,J)-U(I,J))/XCV(I)-GAM(I-1,J)*
     + (U(I,J)-U(I-1,J))/XCV(I-1))/XDIF(I)
	GAMP=GAM(I,J+1)*GAM(I-1,J+1)/(GAM(I,J+1)+GAM(I-1,J+1)+1.E-30)
	GAMP=GAMP+GAM(I,J)*GAM(I-1,J)/(GAM(I,J)+GAM(I-1,J)+1.E-30)
	GAMM=GAM(I,J-1)*GAM(I-1,J-1)/(GAM(I,J-1)+GAM(I-1,J-1)+1.E-30)
	GAMM=GAMM+GAM(I,J)*GAM(I-1,J)/(GAM(I,J)+GAM(I-1,J)+1.E-30)
	CON(I,J)=CON(I,J)+(GAMP*(V(I,J+1)-V(I-1,J+1))-GAMM*(V(I,J)-
     + V(I-1,J)))/(YCV(J)*XDIF(I))
	AP(I,J)=0.
590	CONTINUE
	RETURN
509	CONTINUE
	IF(NF.NE.2)GOTO 519
	DO 595 J=3,M2
	DO 595 I=2,L2
	CON(I,J)=(GAM(I,J)*(V(I,J+1)-V(I,J))/YCV(J)-GAM(I,J-1)*
     + (V(I,J)-V(I,J-1))/YCV(J-1))/YDIF(J)
	GAMP=GAM(I+1,J)*GAM(I+1,J-1)/(GAM(I+1,J)+GAM(I+1,J-1)+1.E-30)
	GAMP=GAMP+GAM(I,J)*GAM(I,J-1)/(GAM(I,J)+GAM(I,J-1)+1.E-30)
	GAMM=GAM(I-1,J)*GAM(I-1,J-1)/(GAM(I-1,J)+GAM(I-1,J-1)+1.E-30)
	GAMM=GAMM+GAM(I,J)*GAM(I,J-1)/(GAM(I,J)+GAM(I,J-1)+1.E-30)
	CON(I,J)=CON(I,J)+(GAMP*(U(I+1,J)-U(I+1,J-1))-GAMM*(U(I,J)-
     + U(I,J-1)))/(XCV(I)*YDIF(J))
	AP(I,J)=0.
595	CONTINUE
	RETURN
519	CONTINUE
	IF(NF.NE.4)GOTO 529
	DO 596 J=2,M2
	DO 596 I=2,L2
	CON(I,J)=0.
596	AP(I,J)=0.
	RETURN
529	CONTINUE
	IF(NF.NE.5)GOTO 539
	DO 598 J=2,M2
	DO 598 I=2,L2
	DUDX=(U(I+1,J)-U(I,J))/XCV(I)
	DVDY=(V(I,J+1)-V(I,J))/YCV(J)
	IF(J.EQ.2)GOTO 581
	IF(J.EQ.M2)GOTO 582
	DUDY=(0.5*(U(I,J+1)-U(I,J-1))+0.5*(U(I+1,J+1)-U(I+1,J-1)))/
     + (YDIF(J)+YDIF(J+1))
	GOTO 583
581	DUDY=(0.5*(U(I,J+1)-U(I,J))+0.5*(U(I+1,J+1)-U(I+1,J)))/(YDIF(J+1))
	GOTO 583
582	DUDY=(0.5*(U(I,J)-U(I,J-1))+0.5*(U(I+1,J)-U(I+1,J-1)))/YDIF(J)
583	CONTINUE
	IF(I.EQ.2)GOTO 584
	IF(I.EQ.L2)GOTO 585
	DVDX=(0.5*(V(I+1,J)-V(I-1,J))+0.5*(V(I+1,J+1)-V(I-1,J+1)))/
     + (XDIF(I)+XDIF(I+1))
	GOTO 586
584	DVDX=(0.5*(V(I+1,J)-V(I,J))+0.5*(V(I+1,J+1)-V(I,J+1)))/XDIF(I+1)
	GOTO 586
585	DVDX=(0.5*(V(I,J)-V(I-1,J))+0.5*(V(I,J+1)-V(I-1,J+1)))/XDIF(I)
586	GEN(I,J)=2.*(DUDX**2+DVDY**2)+(DUDY+DVDX)**2
	CON(I,J)=GEN(I,J)*AMUT(I,J)
	AP(I,J)=-RHO(I,J)*DIS(I,J)/(AKE(I,J)+1.E-30)
598	CONTINUE
	RETURN
539	CONTINUE
	DO 599 J=2,M2
	DO 599 I=2,L2
	CON(I,J)=C1*GEN(I,J)*CMU*RHO(I,J)*AKE(I,J)
	AP(I,J)=-C2*RHO(I,J)*DIS(I,J)/(AKE(I,J)+1.E-30)
599	CONTINUE
	DO 600 J=2,M2
	DISS=CMU*AKE(2,J)**1.5/(0.4*CMU4*XDIF(2))
	CON(2,J)=1.E+30*DISS
	AP(2,J)=0.-1.E+30
600	CONTINUE
	RETURN
	END