plot.f90

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

!C############################################################
PROGRAM PLOTP
!C############################################################
!C     This code produces plots of grid, velocity vectors, and
!C     profiles, contours lines and color fills for any quantity.
!C     The output are postscript files for each page. The code
!C     can easily be adapted for interactive use and screen window
!C     output, but as this is to some extent hardware-dependent,
!C     only postscript output is provided here. The same code
!C     is used for both Cartesian and non-orthogonal grids, and
!C     for both single-grid and multi-grid solutions. Array
!C     dimensions NX and NY should be equal to or greater than
!C     the maximum number of nodes in respective directions on
!C     the finest grid. Up to 128 contours can be ploted. The
!C     plots are saved on files which carry the name+grid number,
!C     eg. VECT1.PS for the plot of velocity vectors from grid 1.
!C
!C                     M. Peric, IfS, Hamburg, 1996
!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 /VEL/ U(NXY),V(NXY)
	DIMENSION P(NXY),T(NXY),PSI(NXY),F2(NXY)
	LOGICAL GRID,VEL,PRES,SF,TEMP,LPROF
	CHARACTER FILRES*10,FILIN*10,FILREF*10
!C
!C==========================================================
!C.....OPEN FILES
!C==========================================================
!C     Input file contains control parameters (what to plot);
!C     results file contains data written by the flow solver
!C     (grid, variable values); file with reference data
!C     contains profiles of some variables (from experiment or
!C     other computations) which shall be plotted along with
!C     profiles from result file. A dummy name is used if no
!C     such data is used. NGRD defines how many data sets the
!C     results file contains (the number of grids from which
!C     data has been recorded). If one wants to plot the results
!C     from say only the forth grid, one line of input is needed
!C     for grids to be skiped, see below.
!C
	PRINT *, ' ENTER NAME OF INPUT FILE:  '
	FILIN = 'plot.inp'
!	READ(*,1) FILIN
	PRINT *, ' ENTER NAME OF RESULTS FILE:  '
	FILRES = 'pstag.plot'
!	READ(*,1) FILRES
	PRINT *, ' ENTER NAME OF FILE WITH REF. DATA:  '
	FILREF = 'test.dat'
!	READ(*,1) FILREF
	PRINT *,' ENTER NGRID:  '
	NGRD = 1
!	READ(*,*) NGRD
!C
	1 FORMAT(A10)
	OPEN (UNIT=5,FILE=FILIN)
	OPEN (UNIT=3,FILE=FILRES,FORM='UNFORMATTED')
	OPEN (UNIT=2,FILE=FILREF)
	REWIND 3
	REWIND 2
	REWIND 5
!C
!C.....LOOP OVER ALL GRIDS: READ GRID DATA AND RESULTS FROM UNIT 3
!C
	DO LK=1,NGRD
!C
		READ(3) ITIM,TIME,NI,NJ,NIM,NJM,NIJ,			&
		&        (X(IJ),IJ=1,NIJ),(Y(IJ),IJ=1,NIJ),	&
		&        (XC(IJ),IJ=1,NIJ),(YC(IJ),IJ=1,NIJ),	&
		&        (PSI(IJ),IJ=1,NIJ),(F2(IJ),IJ=1,NIJ),	&
		&        (U(IJ),IJ=1,NIJ),(V(IJ),IJ=1,NIJ),	&
		&        (P(IJ),IJ=1,NIJ),(T(IJ),IJ=1,NIJ)
!C
		DO I=1,NI
			LI(I)=(I-1)*NJ
		END DO
!C
!C==========================================================
!C.....SOLUTION DOMAIN DIMENSIONS, SCALING FACTOR
!C==========================================================
!C     The solution domain is maped on an area 7 x 7 inches;
!C     Postscript unit is a point (72 points in an inch), and
!C     in order to enable accurate drawing on printers with
!C     a resolution of 1200 x 1200 dot per inch, the coordinates
!C     are scaled appropriately (7 x 1200 = 8400, thus 7 inches
!C     are represented by 8400 dots; a scaling factor of 0.06
!C     is used to convert dots of maximum printer resolution
!C     to postscript units). XPAGE and YPAGE define the ploting
!C     area of 7 x 7 inches.
!C
		XMAX=-1.E20
		YMAX=-1.E20
		XMIN=1.E20
		YMIN=1.E20
		XPAGE=8400.
		YPAGE=8400.
!C
		DO IJ=1,NIJ
			XMAX=MAX(XMAX,X(IJ))
			YMAX=MAX(YMAX,Y(IJ))
			XMIN=MIN(XMIN,X(IJ))
			YMIN=MIN(YMIN,Y(IJ))
		END DO
!C
		XTOT=XMAX-XMIN
		YTOT=YMAX-YMIN
		SCFX=8400.0/XTOT
		SCFY=8400.0/YTOT
		SCFG=MIN(SCFX,SCFY)
		XTOT=XTOT*SCFG
		YTOT=YTOT*SCFG
!C
!C.....GRID COORDINATES AND CENTRAL POINTS COORINATES ON PLOTING PAGE
!C
		DO IJ=1,NIJ
			XC(IJ)=(XC(IJ)-XMIN)*SCFG
			YC(IJ)=(YC(IJ)-YMIN)*SCFG
			X(IJ)=(X(IJ)-XMIN)*SCFG
			Y(IJ)=(Y(IJ)-YMIN)*SCFG
		END DO
!C
		XMIN=0.
		YMIN=0.
		XMAX=XTOT*0.06
		YMAX=YTOT*0.06
!C
!C==========================================================
!C.....CONTROL PARAMETERS AND FLUID PROPERTIES
!C==========================================================
!C     Logical variables GRID, VEL, PRES, TEMP and SF define
!C     whether the grid, velocity vectors or profiles, pressure
!C     contours or colour fill, temperature contours or
!C     colour fill, and streamlines or colour fill will be 
!C     ploted. If one grid level is to be skiped, set all
!C     these variables to FALSE.
!C
		READ(5,*) GRID,VEL,PRES,TEMP,SF
!C
		IF (GRID.OR.VEL.OR.PRES.OR.TEMP.OR.SF) THEN 
!C
!C     UIN and DENSIT are the reference velocity and density
!C     used to scale pressure.
!C
			READ(5,*) UIN,DENSIT
!C
!C     Every IARth vector in I-direction and JARth vector in 
!C     J-direction will be ploted if VEL=TRUE; NPRX profiles
!C     at constant X and NPRY profiles at constant Y will be
!C     ploted; maximum velocity vector of profile width will
!C     be AROMAX inches.
!C
			READ(5,*) IAR,JAR,NPRX,NPRY,AROMAX
			AROMAX=AROMAX*1200.
!C
!C     XBOX and YBOX define the position of the bottom left corner
!C     of the box which contains information about contours
!C     (colour scale interpretation). The values are given in inches
!C     relative to bottom left corner of the ploting area (e.g. if
!C     the geometry is long in X and shorter in Y direction, there
!C     is room within the 7 x 7 inches area for the info-box).
!C
			READ(5,*) XBOX,YBOX
			XBOX=XBOX*1200.
			YBOX=YBOX*1200.
!C
!C     XPRO(I) are the coordinates X at which profiles of variables
!C     are to be ploted. IEPRX(I) defines for which profile there is
!C     comparison data (1 -> comparison data exists, 0 -> no data).
!C     XPRO(I) is to be given in meters.
!C
			IF (NPRX.NE.0) THEN
				READ(5,*) (XPRO(I),I=1,NPRX)
				READ(5,*) (IEPRX(I),I=1,NPRX)
				DO I=1,NPRX
					XPRO(I)=XPRO(I)*SCFG
				END DO
			end if
!C
!C     YPRO(I) are the coordinates Y at which profiles of variables
!C     are to be ploted. IEPRY(I) defines for which profile there is
!C     comparison data (1 -> comparison data exists, 0 -> no data).
!C     YPRO(I) is to be given in meters.
!C
			IF (NPRY.NE.0) THEN
				READ(5,*) (YPRO(J),J=1,NPRY)
				READ(5,*) (IEPRY(J),J=1,NPRY)
				DO J=1,NPRY
					YPRO(J)=YPRO(J)*SCFG
				END DO
			END IF
!C
			LPROF=(NPRX.GT.0.OR.NPRY.GT.0)
!C
!C======================================================
!C.....plot the grid
!C======================================================
!C
			if (grid) then
				call bplot('grid',lk)
				call gridpl
				WRITE(7,*) 'p'
				CLOSE(UNIT=7)
			end if
!c
!C======================================================
!c.....plot velocity vectors 
!C======================================================
!C     NCOL is the number of colours used for colouring
!C     velocity vector according to their magnitude; if
!C     NCOL=1, black vectors are ploted.
!C
			if (vel) then
				read(5,*) NCOL
				call setcol(NCOL)
				call bplot('vect',lk)
				call velpl(NCOL)
				WRITE(7,*) 'p'
				CLOSE(UNIT=7)
!C
!C======================================================
!C......PLOT VELOCITY PROFILES	
!C======================================================
!C
				if (lprof) then
					scfv=1.
					do ij=1,nij
						fi(ij)=u(ij)
					END DO
					call profil('U-velocity     ','upro',lk)
!C
					do ij=1,nij
						fi(ij)=v(ij)
					END DO
					call profil('V-velocity     ','vpro',lk)
				end if
!C
			end if
!c
!C======================================================
!C.....plot isobars
!C======================================================
!C     ICON defines how the contour values are defined:
!C     if ICON=0, the values are read from input file,
!C     otherwise they are calculated between maximum and
!C     minimum value, see the routine CONT. NCON is the 
!C     number of contours to be ploted. ICOL defines 
!C     whether the contours are to be ploted in colour
!C     or black (ICOL=0 -> black, otherwise colour).
!C     IVEC and IGRID define whether the vector or/and
!C     grid plot is to be overlayed over the colour fill
!C     plot (0 -> no overlay, 1 -> yes).
!C
			if (pres) then
				read(5,*) icon,ncon,icol,ivec,igrid
				if (icon.eq.0) then
					read(5,*) (cval(i),i=1,ncon)
				end if
				scfv=1./(densit*uin**2)
				do ij=1,nij
					fi(ij)=p(ij)*scfv
				END DO
!C
				call bplot('isob',lk)
				call cont(icon,ncon,icol,'ISOBARS        ')
				WRITE(7,*) 'p'
				CLOSE(UNIT=7)
!C
!C======================================================
!C.....PLOT PRESSURE COLOR FILL (IF ICOL.GT.0)
!C======================================================
!C
				if (icol.gt.0) then
					call setcol(ncon)
					call bplot('pres',lk)
					call paint(ncon,'ISOBARS        ')
					if(ivec.eq.1) call velpl(0)
					if(igrid.eq.1) call gridpl
					WRITE(7,*) 'p'
					CLOSE(UNIT=7)
				end if
			end if
!c
!C======================================================
!c.....plot temperature contours
!C======================================================
!C
			if (temp) then
				read(5,*) icon,ncon,icol,ivec,igrid
				if (icon.eq.0) then
					read(5,*) (cval(i),i=1,ncon)
				end if
				do ij=1,nij
					fi(ij)=t(ij)
				END DO
!C
				call bplot('isot',lk)
				call cont(icon,ncon,icol,'ISOTHERMS      ')
				WRITE(7,*) 'p'
				CLOSE(UNIT=7)
!C
!C======================================================
!C.....PLOT TEMPERATURE COLOR FILL
!C======================================================
!C
				if (icol.gt.0) then
					call setcol(ncon)
					call bplot('temp',lk)
					call paint(ncon,'ISOTHERMS      ')
					if(ivec.eq.1) call velpl(0)
					if(igrid.eq.1) call gridpl
					WRITE(7,*) 'p'
					CLOSE(UNIT=7)
				end if
!C
!C======================================================
!C.....PLOT TEMPERATURE PROFILES
!C======================================================
!C
				if (lprof) then
					call profil('Temperature    ','tpro',lk)
				end if
			end if
!c
!C======================================================
!c.....calculate stream function values
!C======================================================
!C     Streamfunction values are calculated at CV corners
!C     using mass fluxes throufg cell faces (since the
!C     mass flow rate between two streamlines is constant).
!C     The value at the south-west corner is zero. Since
!C     the contour and colour fill routines operate on
!C     values at XC and YC, and PSI is calculated at 
!C     locations defined by X and Y, the true XC and YC
!C     are copied into P and T arrays and then overwritten
!C     by X and Y.    
!C
			if (sf) then
				PSI(1)=0.
				DO I=1,NIM
					II=LI(I)
					IF(I.NE.1) PSI(II+1)=PSI(II-NJ+1)-F2(II+1)
					DO J=2,NJM
						IJ=II+J
						PSI(IJ)=PSI(IJ-1)+PSI(IJ)
					END DO
				END DO
!C
!C.....SET XC=X AND YC=Y FOR PSI-PLOTS
!C
				do ij=1,nij
					fi(ij)=psi(ij)
					p(ij)=xc(ij)
					t(ij)=yc(ij)
					xc(ij)=x(ij)
					yc(ij)=y(ij)
				END DO
!c
!C=================================================================
!c.....READ CONTROL PARAMETERS AND CONTOUR LEVELS
!C=================================================================
!C     If ICON=0, contour levels are read from input file; otherwise
!C     they are calculated below considering the minimum and 
!C     maximum value at the interior and at boundaries.
!C
				read(5,*) icon,ncon,icol,ivec,igrid
				if (icon.eq.0) then
					read(5,*) (cval(i),i=1,ncon)
				else
!C
!c.....FIND MIN AND MAX VALUES OF PSI IN INTERIOR
!c
					psimax=-1.e20
					psimin=1.e20
					DO I=1,NIM
						DO J=1,NJM