plot.f90

This program incorporates the FV method for solving the Navier-Stokes equations using 2D, Cartesian

!C     be provided on a file in the following format:
!C     The first line is the number of reference points, KPR. It
!C     can be zero, so one does not have to have data for each
!C     variable at all cross-sections.
!C     There follow KPR lines with the Y-coordinate of the reference
!C     point and the variable value at that point, FIPR. A diamond
!C     symbol is ploted at each point.
!C
			IF (IEPRX(K).NE.0) THEN
				WRITE(7,*) '16 w'
				READ(2,*) KPR
				IF (KPR.NE.0) THEN
					DO KK=1,KPR
						READ(2,*) YPR(KK),FIPR(KK)
					END DO
!C
					DO KK=1,KPR
						YPR(KK)=YPR(KK)*SCFG
						FIPR(KK)=XINT+FIPR(KK)*SCF*SCFV
						WRITE(7,*) int(FIPR(KK)-100.),int(YPR(KK)),' m'
						WRITE(7,*) int(FIPR(KK)),int(YPR(KK)-100.),' l'
						WRITE(7,*) int(FIPR(KK)+100.),int(YPR(KK)),' l'
						WRITE(7,*) int(FIPR(KK)),int(YPR(KK)+100.),' l'
						WRITE(7,*) int(FIPR(KK)-100.),int(YPR(KK)),' l s'
					END DO
				end if
			end if
		END DO
!C
!C------------------------------------------------------------
!C.....WRITE THE TITLE AND PROFILE SCALE
!C------------------------------------------------------------
!C
		TSIZE=300.
		x1=0.35*xtot
		y1=ytot+2.*TSIZE
		x2=x1+aromax
		y2=y1
		write(7,*) x1,y1,' m ',x2,y2,' l s'
		x2=x2+tsize
		y2=y2-0.5*tsize
		write(7,*) '/Helvetica findfont 330.00 scalefont setfont'
		write(7,*) int(x2),int(y2),'  m'
		write(tt,'(5h  =  ,1pe8.2,4h m/s)') fimax
		write(7,*) '(',tt,') show'
		y1=y1+2.*tsize
		write(7,*) int(x1),int(y1),'  m'
		write(7,*) '(',title,') show'
		WRITE(7,*) 's p'
		CLOSE(UNIT=7)
!C
	end if
!C
!C===============================================================
!C.....PROFILES AT SPECIFIED Y-CROSS-SECTION
!C===============================================================
!C
	IF (NPRY.NE.0) THEN
		CALL BPLOT(FILN,LK)
!C
		DO K=1,NPRY
			YINT=YPRO(K)
			DO I=1,NIM
				DO J=1,NJM
					IJ=LI(I)+J
					CALL PROFY(IJ,IP,YINT)
!C
					IF (IP.EQ.2) THEN
						WRITE(7,*) '20 w ',INT(XP(1)),INT(YINT+YP(1)),' m ',	&
						&                    INT(XP(2)),INT(YINT+YP(2)),' l s'
						WRITE(7,*) '10 w ',INT(XP(1)),INT(YINT),' m ',			&
						&                    INT(XP(2)),INT(YINT),' l s'
					end if
				END DO
			END DO
!C
!C.....COMPARISON WITH REFERENCE DATA
!C
			IF (IEPRY(K).NE.0) THEN
				WRITE(7,*) '16 w'
				READ(2,*) KPR
				IF (KPR.NE.0) THEN
					READ(2,*) (XPR(KK),KK=1,KPR)
					READ(2,*) (FIPR(KK),KK=1,KPR)
					DO KK=1,KPR
						XPR(KK)=XPR(KK)*SCFG
						FIPR(KK)=XINT+FIPR(KK)*SCF*SCFV
						WRITE(7,*) int(XPR(KK)),int(FIPR(KK)-100.),' m'
						WRITE(7,*) int(XPR(KK)-100.),int(FIPR(KK)),' l'
						WRITE(7,*) int(XPR(KK)),int(FIPR(KK)+100.),' l'
						WRITE(7,*) int(XPR(KK)+100.),int(FIPR(KK)),' l'
						WRITE(7,*) int(XPR(KK)),int(FIPR(KK)-100.),' l s'
					END DO
				end if
			end if
		END DO
!C
!C.....WRITE THE TITLE AND PROFILE SCALE
!C
		x1=0.35*xtot
		y1=ytot+2.*tsize
		x2=x1+aromax
		y2=y1
		write(7,*) x1,y1,' m ',x2,y2,' l s'
		x2=x2+tsize
		y2=y2-0.5*tsize
		write(7,*) '/Helvetica findfont 330.00 scalefont setfont'
		write(7,*) int(x2),int(y2),'  m'
		write(tt,'(5h  =  ,1pe8.2,4h m/s)') fimax
		write(7,*) '(',tt,') show'
		y1=y1+2.*tsize
		write(7,*) int(x1),int(y1),'  m'
		write(7,*) '(',title,') show'
		WRITE(7,*) 's p'
		CLOSE(UNIT=7)
!C
	end if
!C
	
END subroutine
!C
!C
!C##########################################################
SUBROUTINE PROFX(IJ,IP,XINT)
!C##########################################################
!C     This routine checks if the line X=XINT crosses an
!C     element, defined by the south-west corner index IJ.
!C     For each side of the element, it calls routine PROPX
!C     to find the crossing point. The coordinates of the
!C     crossing point are stored as XP and YP. IP is the
!C     number of crossing points found.
!C==========================================================
	PARAMETER (nx=101,ny=101,NXY=NX*NY)
	COMMON /IPAR/ NI,NJ,NIM,NJM,NIJ,LI(NX),NPRX,NPRY,IAR,	& 
	&       JAR,IEPRX(10),IEPRY(10)
	COMMON /RPAR/ XMIN,YMIN,XMAX,YMAX,XTOT,YTOT,XBOX,YBOX,	& 
	&       XPAGE,YPAGE,AROMAX,SCFG,SCFV,					& 
	&       XPRO(10),YPRO(10),CVAL(128),FI1,FI2
	COMMON /GEO/ X(NXY),Y(NXY),XC(NXY),YC(NXY),FI(NXY),		& 
	&       CVX(20),CVY(20)
	COMMON /PRO/ XP(2),YP(2)
!C
!C.....FIND IF XINT CROSSES THE ELEMENT WHOSE SW-CORNER IS IJ
!C
	IP=0
	XMN=MIN(XC(IJ),XC(IJ+NJ),XC(IJ+NJ+1),XC(IJ+1))
	XMX=MAX(XC(IJ),XC(IJ+NJ),XC(IJ+NJ+1),XC(IJ+1))
	IF (XINT.GE.XMN.AND.XINT.LE.XMX) THEN
		CALL PROPX(IJ,IJ+NJ,IP,XINT)
		CALL PROPX(IJ+NJ,IJ+NJ+1,IP,XINT)
		IF(IP.NE.2) CALL PROPX(IJ+NJ+1,IJ+1,IP,XINT)
		IF(IP.NE.2) CALL PROPX(IJ+1,IJ,IP,XINT)
	end if
	
END subroutine
!C
!C#############################################################
SUBROUTINE PROPX(IJ1,IJ2,IP,XINT)
!C#############################################################
!C     This routine finds the crossing point of the line X=XINT
!C     with an element side defined by node indices at its ends,
!C     IJ1 and IJ2. IP is incremented by one if the crossing
!C     point is found, and the coordinates are stored in XP(IP)
!C     and YP(IP).
!C==============================================================
	PARAMETER (nx=101,ny=101,NXY=NX*NY)
	COMMON /IPAR/ NI,NJ,NIM,NJM,NIJ,LI(NX),NPRX,NPRY,IAR,	& 
	&       JAR,IEPRX(10),IEPRY(10)
	COMMON /RPAR/ XMIN,YMIN,XMAX,YMAX,XTOT,YTOT,XBOX,YBOX,	& 
	&       XPAGE,YPAGE,AROMAX,SCFG,SCFV,					& 
	&       XPRO(10),YPRO(10),CVAL(128),FI1,FI2
	COMMON /GEO/ X(NXY),Y(NXY),XC(NXY),YC(NXY),FI(NXY),		& 
	&       CVX(20),CVY(20)
	COMMON /PRO/ XP(2),YP(2)
!C
!C.....SEARCH FOR CROSSING POINT BETWEEN IJ1 AND IJ2 NODES
!C
	IF (XINT.GE.MIN(XC(IJ1),XC(IJ2)).AND.		&
		&   XINT.LE.MAX(XC(IJ1),XC(IJ2))) THEN
		IP=IP+1
		FAC=(XINT-XC(IJ1))/(XC(IJ2)-XC(IJ1)+1.E-20)
		YP(IP)=YC(IJ1)+FAC*(YC(IJ2)-YC(IJ1))
		XP(IP)=FI(IJ1)+FAC*(FI(IJ2)-FI(IJ1))
	end if
	
END subroutine
!C         
!C
!C##########################################################
SUBROUTINE PROFY(IJ,IP,YINT)
!C##########################################################
!C     This routine is analogous to PROFX, only for Y=YINT.
!C==========================================================
	PARAMETER (nx=101,ny=101,NXY=NX*NY)
	COMMON /IPAR/ NI,NJ,NIM,NJM,NIJ,LI(NX),NPRX,NPRY,IAR,	&
	&       JAR,IEPRX(10),IEPRY(10)
	COMMON /RPAR/ XMIN,YMIN,XMAX,YMAX,XTOT,YTOT,XBOX,YBOX,	&
	&       XPAGE,YPAGE,AROMAX,SCFG,SCFV,					&
	&       XPRO(10),YPRO(10),CVAL(128),FI1,FI2
	COMMON /GEO/ X(NXY),Y(NXY),XC(NXY),YC(NXY),FI(NXY),		&
	&       CVX(20),CVY(20)
	COMMON /PRO/ XP(2),YP(2)
!C
!C.....FIND IF YINT CROSSES THE CELL WHOSE SW-CORNER IS IJ
!C
	IP=0
	YMN=MIN(YC(IJ),YC(IJ+NJ),YC(IJ+NJ+1),YC(IJ+1))
	YMX=MAX(YC(IJ),YC(IJ+NJ),YC(IJ+NJ+1),YC(IJ+1))
	IF (YINT.GE.YMN.AND.YINT.LE.YMX) THEN
		CALL PROPY(IJ,IJ+NJ,IP,YINT)
		CALL PROPY(IJ+NJ,IJ+NJ+1,IP,YINT)
		IF(IP.NE.2) CALL PROPY(IJ+NJ+1,IJ+1,IP,YINT)
		IF(IP.NE.2) CALL PROPY(IJ+1,IJ,IP,YINT)
	end if
	
END subroutine
!C
!C##########################################################
SUBROUTINE PROPY(IJ1,IJ2,IP,YINT)
!C##########################################################
!C     This routine is analogous to PROPX, only for Y=YINT.
!C==========================================================
	PARAMETER (nx=101,ny=101,NXY=NX*NY)
	COMMON /IPAR/ NI,NJ,NIM,NJM,NIJ,LI(NX),NPRX,NPRY,IAR,	&
	&       JAR,IEPRX(10),IEPRY(10)
	COMMON /RPAR/ XMIN,YMIN,XMAX,YMAX,XTOT,YTOT,XBOX,YBOX,	&
	&       XPAGE,YPAGE,AROMAX,SCFG,SCFV,					&
	&       XPRO(10),YPRO(10),CVAL(128),FI1,FI2
	COMMON /GEO/ X(NXY),Y(NXY),XC(NXY),YC(NXY),FI(NXY),		&
	&       CVX(20),CVY(20)
	COMMON /PRO/ XP(2),YP(2)
!C
!C.....SEARCH FOR CROSSING POINT BETWEEN IJ1 AND IJ2 NODES
!C
	IF (YINT.GE.MIN(YC(IJ1),YC(IJ2)).AND.	&
		&   YINT.LE.MAX(YC(IJ1),YC(IJ2))) THEN
		IP=IP+1
		FAC=(YINT-YC(IJ1))/(YC(IJ2)-YC(IJ1)+1.E-20)
		XP(IP)=XC(IJ1)+FAC*(XC(IJ2)-XC(IJ1))
		YP(IP)=FI(IJ1)+FAC*(FI(IJ2)-FI(IJ1))
	end if
	
END subroutine
!C
!C
!C###############################################################
SUBROUTINE ARROW(XC1,YC1,X,Y)
!C###############################################################
!C     This routine draws the arrow representing velocity vector.
!C     XC1 and YC1 are page coordinates of vector start, X and Y
!C     are coordinates of vector end. The arrow is long 20% of the
!C     vector length and wide 5% of the vector length.
!C===============================================================
!C
	VL=SQRT((X-XC1)**2+(Y-YC1)**2)
	IF (VL.GT.1.E-8) THEN
		WRITE(7,*) int(XC1),int(YC1),' m',int(X),int(Y),' l'
		DX=X-XC1
		DY=Y-YC1
		X1=X-0.2*DX
		Y1=Y-0.2*DY
		DA=0.025*VL
		SAL=DY/VL
		CAL=DX/VL
		DX=DA*SAL
		DY=DA*CAL
		X1=X1-DX
		Y1=Y1+DY
		X2=X1+2.*DX
		Y2=Y1-2.*DY
		WRITE(7,*) int(X1),int(Y1),' l',int(X2),int(Y2),' l',	&
		&             int(X),int(Y),' l s'
	end if
	
END subroutine
!C
!C
!C###########################################################
SUBROUTINE GRIDPL
!C###########################################################
!C     This routine plots the grid lines other than boundary
!C     lines.
!C===========================================================
	PARAMETER (nx=101,ny=101,NXY=NX*NY)
	COMMON /IPAR/ NI,NJ,NIM,NJM,NIJ,LI(NX),NPRX,NPRY,IAR,	&
	&       JAR,IEPRX(10),IEPRY(10)
	COMMON /RPAR/ XMIN,YMIN,XMAX,YMAX,XTOT,YTOT,XBOX,YBOX,	&
	&       XPAGE,YPAGE,AROMAX,SCFG,SCFV,					&
	&       XPRO(10),YPRO(10),CVAL(128),FI1,FI2
	COMMON /GEO/ X(NXY),Y(NXY),XC(NXY),YC(NXY),FI(NXY),		&
	&       CVX(20),CVY(20)
!C
	WRITE(7,*) '10 w'
!C
!C.....PLOT I-LINES (J=CONST)
!C
	DO J=2,NJM-1
		IJ=LI(1)+J
		WRITE(7,'(2I5,A3)') INT(X(IJ)),INT(Y(IJ)),' m '
!C
		DO I=2,NIM
			IJ=LI(I)+J
			WRITE(7,'(2I5,A3)') INT(X(IJ)),INT(Y(IJ)),' l '
		END DO
		WRITE(7,*) 's'
	END DO
!C
!C.....PLOT  J - LINES (I=CONST)
!C
	DO I=2,NIM-1
		IJ=LI(I)+1
		WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' m '
!C
		DO J=2,NJM
			IJ=LI(I)+J
			WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' l '
		END DO
		WRITE(7,*) 's'
	END DO
!C
	
END subroutine
!C
!C
!C##########################################################
SUBROUTINE BPLOT(HEAD,L)
!C##########################################################
!C     This routine opens the postscript file with the name
!C     specified by HEAD and L (e.g. isob1.ps if HEAD=isob,
!C     and L=1), writes the postscript header to it and plots
!C     the boundary grid lines with a thick line. This is
!C     the first routine that has to be called for any plot.
!C==========================================================
	PARAMETER (nx=101,ny=101,NXY=NX*NY)
	COMMON /IPAR/ NI,NJ,NIM,NJM,NIJ,LI(NX),NPRX,NPRY,IAR,	&
	&       JAR,IEPRX(10),IEPRY(10)
	COMMON /RPAR/ XMIN,YMIN,XMAX,YMAX,XTOT,YTOT,XBOX,YBOX,	&
	&       XPAGE,YPAGE,AROMAX,SCFG,SCFV,					&
	&       XPRO(10),YPRO(10),CVAL(128),FI1,FI2
	COMMON /GEO/ X(NXY),Y(NXY),XC(NXY),YC(NXY),FI(NXY),		&
	&       CVX(20),CVY(20)
	COMMON /RGB/ R(255),G(255),B(255)
	CHARACTER FILOUT*8,HEAD*4
!C
!C.....DEFINE PLOT AREA AND COORDINATE ORIGIN
!C
	WRITE(FILOUT,'(a4,I1,3H.ps)') HEAD,L
	OPEN(UNIT=7,FILE=FILOUT)
	REWIND 7
	CALL PSHEAD(XMIN,XMAX,YMIN,YMAX)
	WRITE(7,*) '0. 0. 0. col'
!C
!C.....Move to south-west corner
!C
	WRITE(7,*) '40 w'
	IJ=LI(1)+1
	WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' m '
!C
!C.....Draw west boundary
!C
	DO J=2,NJM
		IJ=LI(1)+J
		WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' l '
	END DO
!C
!C.....Draw north boundary
!C
	DO I=2,NIM
		IJ=LI(I)+NJM
		WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' l '
	END DO
!C
!C.....Draw east boundary
!C
	DO J=NJM-1,1,-1
		IJ=LI(NIM)+J
		WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' l '
	END DO
!C
!C.....Draw south boundary
!C