plot.f90

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

	DO I=NIM-1,1,-1
		IJ=LI(I)+1
		WRITE(7,'(2I5,A3)') int(X(IJ)),int(Y(IJ)),' l '
	END DO
!C
	WRITE(7,*) 's'
!C
	
END subroutine
!C
!C############################################################
subroutine cont(icon,ncon,icol,title)
!C############################################################
!C     This routine plots contour lines of a constant value of
!C     the variable stored in the array FI(NXY). It operates on
!C     elements defined by CV centers. Contour lines are
!C     made of segments of straight lines connecting points at
!C     element edges where the interpolated varible value equals
!C     the contour level. Since these segments are small, the
!C     lines apear smooth for grids with more than 20 nodes in
!C     one direction.
!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 title*15
!C
!C.... find extrema
!C
	fimax=-1.e20
	fimin=1.e20
	do ij=1,nij
		fimax=max(fimax,fi(ij))
		fimin=min(fimin,fi(ij))
	end do
!C
!C.....set contour levels (if icon.ne.0; linear distribution)
!C
	if (icon.ne.0) then
		dfi=(fimax-fimin)/real(ncon)
		cval(1)=fimin+0.5*dfi
		do nc=2,ncon
			cval(nc)=cval(nc-1)+dfi
		end do
	end if
!C
!C.....set line thicknes and colors (black if ICOL=0)
!C   
	WRITE(7,*) '20 w '
	if (icol.eq.0) then
		write(7,*) '0. 0. 0. col'
	else
		call setcol(ncon)
	end if
!C
!C=====================================================================
!C.....Search each quadrilateral element (IJ node is south-west corner)
!C=====================================================================
!C     There are NIM elements in I and NJM elements in J direction
!C     (boundary nodes are also used in defining the elements).
!C     IJ is the south-west, IJE the south-east, IJNE the north-east
!C     and IJN the north-west corner of the element. The variable
!C     values are available at corners (these are CV centers). If the
!C     contour level is within the minimum and maximum value at the
!C     four vertices of the element, the four element edges are
!C     checked to find out if the contour level falls between values
!C     at ends (routine NEWC). If so, the location is found by linear
!C     interpolation, the coordinates are stored in CVX and CVY arrays
!C     and the counter IVC is increased. The lines between points are
!C     then drawn. If more than two points are found, the line is
!C     closed from last to first point.
!C
	do k=1,ncon
		if(icol.gt.0) WRITE(7,*) R(K),G(K),B(K),' col'
!C
		do i=1,nim
			do j=1,njm
				ij=li(i)+j
				ije=ij+nj
				ijn=ij+1
				ijne=ije+1
!C
!C.....Search for contour AND PLOT IT
!C
				fi1=cval(k)
				fimn=min(fi(ij),fi(ije),fi(ijn),fi(ijne))
				fimx=max(fi(ij),fi(ije),fi(ijn),fi(ijne))
				if (.not.((fi1.lt.fimn).or.(fi1.gt.fimx))) then
					ivc=0
					call newc(ij,ije,ivc)
					call newc(ije,ijne,ivc)
					call newc(ijne,ijn,ivc)
					call newc(ijn,ij,ivc)
					if (ivc.gt.1) then
						WRITE(7,*) int(cvx(1)),int(cvy(1)),' m'
						do kk=2,ivc
							WRITE(7,*) int(cvx(kk)),int(cvy(kk)),' l'
						end do
						if(ivc.gt.2) write(7,*) int(cvx(1)),int(cvy(1)),' l'
						WRITE(7,*) 's'
					end if
				end if
			end do
		end do
!C
	end do
!C
!C.....PLOT THE BOX with contour levels and colour scale
!C
	write(7,*) '0. 0. 0. col'
	call plbox(ncon,icol,title)
!C
	
end subroutine
!C
!C############################################################
subroutine newc(ij1,ij2,ivc)
!C############################################################
!C     This routine finds the contour point on one element
!C     edge (it checks also if the corner value at IJ1 equals
!C     the contour level value).
!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
!C.....Check if value at IJ1 equals contour value
!C
	if (fi(ij1).eq.fi1) then
		ivc=ivc+1
		cvx(ivc)=xc(ij1)
		cvy(ivc)=yc(ij1)
!C
!C.....Check if contur value FI1 lies between values at IJ1 and IJ2;
!C     if so, interpolate to find location between nodes IJ1 and IJ2
!C
		elseif ((fi1.ge.min(fi(ij1),fi(ij2))).and.	&
		&   (fi1.le.max(fi(ij1),fi(ij2)))) then
		ivc=ivc+1
		fac1=(fi1-fi(ij1))/(fi(ij2)-fi(ij1)+1.e-20)
		cvx(ivc)=xc(ij1)+fac1*(xc(ij2)-xc(ij1))
		cvy(ivc)=yc(ij1)+fac1*(yc(ij2)-yc(ij1))
	end if
!C
	
end subroutine
!C
!c#############################################################
subroutine plbox(ncon,icol,title)
!c#############################################################
!C     This routine plots the box with contour levels and colour 
!C     scale.
!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 /RGB/ R(255),G(255),B(255)
	character title*15,tt*9
!C
	x1=xbox
	y1=ybox+3300.
	write(7,*) '/Helvetica findfont 250.00 scalefont setfont'
	write(7,*) int(x1),int(y1),' m (',title,') show'
	y1=y1-300.
	write(7,*) '/Helvetica findfont 200.00 scalefont setfont'
	write(7,*) int(x1),int(y1),' m (Contour levels:) show'
	write(7,*) '100 w'
	x2=x1+700.
	y1=y1-270.
!c
	istep=(ncon+5)/10
!c
	do n=1,ncon,istep
		if(icol.eq.1) write(7,*) r(n),g(n),b(n),' col'
		write(7,*) int(x1),int(y1),' m ',int(x2),int(y1),' l s' 
		write(tt,'(1pe9.2)') cval(n)
		write(7,*) int(x2+200.),int(y1-80.),' m ','(',tt,') show'
		y1=y1-250.
	end do
!c
	
end subroutine
!c
!c#############################################################
subroutine paint(ncon,title)
!c#############################################################
!C     This routine paints each element by filling one colour 
!C     between two contour levels. It operates on the same
!C     element as the contouring routine. Two contour levels
!C     are considered, and the polygon is defined using contour
!C     lines and element edges.
!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 title*15
!c      
!c.....Find extrema, choose contour colours, plot the box
!c
	fimax=-1.e20
	fimin=1.e20
	do ij=1,nij
		fimax=max(fimax,fi(ij))
		fimin=min(fimin,fi(ij))
	end do
	dfi=(fimax-fimin)/real(ncon)
!c
	cval(1)=fimin+0.5*dfi
	do nc=2,ncon
		cval(nc)=cval(nc-1)+dfi
	end do
!c
	call setcol(ncon)
	icol=1
	call plbox(ncon,icol,title)
!c
!C============================================================
!c.....Paint each element (south-west corner is node IJ)
!C============================================================
!c     One fills first the colour between FIMIN and the first
!C     contour level, then between contours, and finally between
!C     the last contour level and PHIMAX.
!C
	WRITE(7,*) '20 w '
	fi1=fimin
	do k=1,ncon
		WRITE(7,*) R(K),G(K),B(K),' col'
		fi2=fi1+dfi
		if(k.eq.ncon) fi2=fimax
!C
		do i=1,nim
			do j=1,njm
				ij=li(i)+j
				ije=ij+nj
				ijn=ij+1
				ijne=ije+1
!c
!c.....Search for filling range between FI1 and FI2
!c
				fimn=min(fi(ij),fi(ije),fi(ijn),fi(ijne))
				fimx=max(fi(ij),fi(ije),fi(ijn),fi(ijne))
				IF (.not.(((fi1.lt.fimn).and.(fi2.lt.fimn)).or.	&
					&            ((fi1.gt.fimx).and.(fi2.gt.fimx)))) then 
					ivc=0
					call newvc(ij,ije,ivc)
					call newvc(ije,ijne,ivc)
					call newvc(ijne,ijn,ivc)
					call newvc(ijn,ij,ivc)
					if (ivc.gt.2) then
						write(7,*) 'n ',int(cvx(1)),int(cvy(1)),' m'
						do l=2,ivc
							write(7,*) int(cvx(l)),int(cvy(l)),' l'
						end do
						write(7,*) int(cvx(1)),int(cvy(1)),' l cp f s'
					end if
				END IF
			END DO
		END DO
!C
		fi1=fi2
	END DO
!C
	
end subroutine
!c
!c#############################################################
subroutine newvc(ij1,ij2,ivc)
!c#############################################################
!C     This routine defines and fills with colour a polygon
!C     bounded by two contours (FI1 and FI2) and element edges.
!C     If an element corner is within the range, it becomes a
!C     vertex in the polygon. When points corresponding to both
!C     FI1 and FI2 are found on one edge, they are sorted out so
!C     that polygon edges define a closed surface (one searces
!C     along edges of an element counterclockwise; on one edge,
!C     from IJ1 node to IJ2 node).
!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
!c.....Check if node IJ1 is within painted range 
!c
	if (fi(ij1).ge.fi1.and.fi(ij1).le.fi2) then
		ivc=ivc+1
		cvx(ivc)=xc(ij1)
		cvy(ivc)=yc(ij1)
	end if
!c
!c.....Value FI1 (lower bound) between nodes IJ1 and IJ2
!c
	fac1=0.
	if ((fi1.ge.min(fi(ij1),fi(ij2))).and.			&
		&   (fi1.le.max(fi(ij1),fi(ij2)))) then
		ivc=ivc+1
		fac1=(fi1-fi(ij1))/(fi(ij2)-fi(ij1)+1.e-20)
		cvx(ivc)=xc(ij1)+fac1*(xc(ij2)-xc(ij1))
		cvy(ivc)=yc(ij1)+fac1*(yc(ij2)-yc(ij1))
	end if
!c
!c....Value FI2 (upper bound) between nodes IJ1 and IJ2
!c
	fac2=0.
	if ((fi2.ge.min(fi(ij1),fi(ij2))).and.		&
		&   (fi2.le.max(fi(ij1),fi(ij2)))) then
		ivc=ivc+1
		fac2=(fi2-fi(ij1))/(fi(ij2)-fi(ij1)+1.e-20)
		cvx(ivc)=xc(ij1)+fac2*(xc(ij2)-xc(ij1))
		cvy(ivc)=yc(ij1)+fac2*(yc(ij2)-yc(ij1))
!c
!c.....if both points found, sort them out
!c
		if (fac2.lt.fac1) then
			x1=cvx(ivc-1)
			y1=cvy(ivc-1)
			cvx(ivc-1)=cvx(ivc)
			cvy(ivc-1)=cvy(ivc)
			cvx(ivc)=x1
			cvy(ivc)=y1
		end if
	end if
!C
	
end subroutine
!C
!C
!C##########################################################
SUBROUTINE PSHEAD(XMIN,XMAX,YMIN,YMAX)
!C##########################################################
!C     This routine writes postscript header to the plot
!C     file. Bounding box coveres only the solution domain
!C     boundary, not the contour information. 
!C==========================================================
!C
	X1=XMIN+50.
	Y1=YMIN+50.
	X2=XMAX+55.
	Y2=YMAX+55.
!C
	WRITE(7,*) '%!PS-Adobe-2.0'
	WRITE(7,*) '%%Creator: PLOTP'
	WRITE(7,*) '%%BoundingBox: ',X1,Y1,X2,Y2
	WRITE(7,*) '%%EndComments'
	WRITE(7,*) '/c {currentpoint} def /f {fill} def '
	WRITE(7,*) '/gr {grestore} def /gs {gsave} def /l {lineto} def '
	WRITE(7,*) '/m {moveto} def /n {newpath} def /p {showpage} def '
	WRITE(7,*) '/s {stroke} def /sg {setgray} def '
	WRITE(7,*) '/w {setlinewidth} def /cp {closepath} def'
	WRITE(7,*) '/col {setrgbcolor} def'
	WRITE(7,*) '50 50 translate 0.06 0.06 scale '
	WRITE(7,*) '1 setlinecap 1 setlinejoin '
!C
	
END subroutine