plot.f90

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

							IJ=LI(I)+J
							psimax=max(psimax,fi(ij))
							psimin=min(psimin,fi(ij))
						END DO
					END DO
					write(*,*) 'psimin = ',psimin,' , psimax = ',psimax
!C
!C.....FIND MIN AND MAX VALUES OF PSI AT BOUNDARY
!C
					PSBMIN=1.e20
					PSBMAX=-1.E20
					DO I=1,NIM
						IJ=LI(I)+1
						PSBMIN=MIN(PSBMIN,FI(IJ))
						PSBMAX=MAX(PSBMIN,FI(IJ))
						IJ=LI(I)+NJM
						PSBMIN=MIN(PSBMIN,FI(IJ))
						PSBMAX=MAX(PSBMIN,FI(IJ))
					END DO
!C
					DO J=1,NJM
						IJ=LI(1)+J
						PSBMIN=MIN(PSBMIN,FI(IJ))
						PSBMAX=MAX(PSBMIN,FI(IJ))
						IJ=LI(NIM)+J
						PSBMIN=MIN(PSBMIN,FI(IJ))
						PSBMAX=MAX(PSBMIN,FI(IJ))
					END DO
					write(*,*) 'psbmin = ',psbmin,' , psbmax = ',psbmax
!C
!C=================================================================
!C.....RANGES OF PSI VALUES, SET NO. OF CONTOURS IN RECIRC. REGIONS
!C=================================================================
!C     NCMIN contours are defined between the minimum value of PSI
!C     at boundary and in the interior; it is chosen as a fraction
!C     of NCON, depending on how big is the range DPSMIN.
!C   
					DELPSB=PSBMAX-PSBMIN
					DPSMIN=PSBMIN-PSIMIN
					DPSMAX=PSIMAX-PSBMAX
					DM1=MAX(DELPSB,DPSMAX)
!C
					IF (DPSMIN.LT.0.01*DM1) THEN
						NCMIN=0
					elseif(DPSMIN.GE.0.01*DM1.AND.DPSMIN.LT.0.05*DM1) THEN
						NCMIN=NCON/4
					elseif(DPSMIN.GE.0.05*DM1.AND.DPSMIN.LT.0.1*DM1) THEN
						NCMIN=NCON/3
					elseif(DPSMIN.GE.0.1*DM1.AND.DPSMIN.LT.0.5*DM1) THEN
						NCMIN=NCON/2
					ELSE
						NCMIN=NCON
					END IF
!C
!C     NCMAX contours are defined between the maximum value of PSI
!C     at boundary and in the interior; it is chosen as a fraction
!C     of NCON, depending on how big is the range DPSMAX.
!C   
					DM1=MAX(DELPSB,DPSMIN)
					IF (DPSMAX.LT.0.01*DM1) THEN
						NCMAX=0
					elseif(DPSMAX.GE.0.01*DM1.AND.DPSMAX.LT.0.05*DM1) THEN
						NCMAX=NCON/4
					elseif(DPSMAX.GE.0.05*DM1.AND.DPSMAX.LT.0.1*DM1) THEN
						NCMAX=NCON/3
					elseif(DPSMAX.GE.0.1*DM1.AND.DPSMAX.LT.0.5*DM1) THEN
						NCMAX=NCON/2
					ELSE
						NCMAX=NCON
					END IF
!C
!c.....SET STREAMFUNCTION CONTOUR LEVELS
!c
					LC=0
					DPSI=DPSMIN/REAL(NCMIN+1)    
					DO L=1,NCMIN
						lc=lc+1
						cval(lc)=psimin+DPSI*real(l)
					END DO 
					lc=lc+1
					cval(lc)=PSBMIN-1.E-20
!C
!C     Boundary values are defined as contour lines. If
!C     the difference between maximum and minimum values at 
!C     boundary is small compared to the full range, no
!C     contour values are set between these two values;
!C     otherwise, NCON values are defined.
!C
					IF (DELPSB.GT.0.1*(DPSMIN+DPSMAX)) THEN
						DPSI=DELPSB/REAL(NCON+1)
						DO L=1,NCON
							lc=lc+1
							cval(lc)=PSBMIN+DPSI*real(l)
						END DO
					END IF
!C
					DPSI=DPSMAX/REAL(NCMAX+1)
					DO L=1,NCMAX
						LC=LC+1
						CVAL(LC)=PSIMAX-DPSI*REAL(L)
					END DO
!C
					LC=LC+1
					CVAL(LC)=PSBMAX+1.E-20
					icon=0
					ncon=lc
				end if
!C
!C======================================================
!C.....PLOT STREAMLINES
!C======================================================
!C
				call bplot('strl',lk)
				call cont(icon,ncon,icol,'STREAMLINES   ')
				WRITE(7,*) 'p'
				CLOSE(UNIT=7)
!C
!C======================================================
!C.....PLOT STREAMLINE COLOR FILL
!C======================================================
!C
				if (icol.gt.0) then
					call setcol(ncon)
					call bplot('strc',lk)
					call paint(ncon,'STREAMLINES    ')
					if (ivec.eq.1) then
						icol=1
						do ij=1,nij
							xc(ij)=p(ij)
							yc(ij)=t(ij)
						end do
						call velpl(0)
					end if
					if(igrid.eq.1) call gridpl
					WRITE(7,*) 'p'
					CLOSE(UNIT=7)
				end if
			end if
!C
		END IF
	END DO
!C
	
END Program
!C
!C
!C##################################################################
SUBROUTINE VELPL(NCOL)
!C##################################################################
!C     This routine plots velocity vectors (every IARth and JARth
!C     value in I and J direction, respectively). Colour is defined
!C     by RGB values, which range between 0. and 1. for each of 
!C     the three base colours. The appropriate values for R, G and B
!C     are calculated in the routine SETRGB.
!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)
	COMMON /RGB/ R(255),G(255),B(255)
	CHARACTER TT*17
!c
!c.....calculate scaling factor for velocities (max. value = aromax)
!c
	vmax=-1.0e10
	vmin=1.e10
	do ij=1,nij
		vmean=sqrt(u(ij)**2+v(ij)**2)
		vmin=min(vmin,vmean)
		vmax=max(vmax,vmean)
	end do
	scf=aromax/(vmax+1.e-15)
!c
!c.....Set colour levels (for NCOL=1, black is used)
!c
	delv=(vmax-vmin)/max(real(ncol),1.)
	do n=1,ncol+1
		cval(n)=vmin+real(n-1)*delv
	end do
	WRITE(7,*) '10 w'
	IF(NCOL.EQ.1) WRITE(7,*) '0. 0. 0. col'
!c
!C=======================================================
!c.....Plot vectors at each  IAR-th and JAR-th position
!C=======================================================
!c
	do i=2,nim,iar
		do j=2,njm,jar
			ij=li(i)+j
			x1=u(ij)*scf+xc(ij)
			y1=v(ij)*scf+yc(ij)
!c
!C.....Assign colour according to velocity vector magnitude
!C
			vm=sqrt(u(ij)**2+v(ij)**2)
			if (ncol.gt.1) then
				do n=2,ncol+1
					if (vm.le.cval(n).and.vm.ge.cval(n-1)) THEN
						WRITE(7,*) R(N-1),G(N-1),B(N-1),' col'
					end if
				end do
			end if
!C
!C....Draw the arrow
!C 
			call arrow(xc(ij),yc(ij),x1,y1)
		end do
	end do
!c
!C=======================================================
!c.....write title and vector scale to postscript file
!C=======================================================
!c
	tsize=300.
	x1=0.35*xtot
	y1=ytot+3.*tsize
	x2=x1+aromax
	y2=y1
	call arrow(x1,y1,x2,y2)
	x2=x2+tsize
	y2=y2-0.5*tsize
	write(7,*) int(x2),int(y2),'  m'
	write(7,*) '/Times-Roman findfont 250.00 scalefont setfont'
	write(tt,'(4h =  ,1pe8.2,5h m/s )') vmax
	write(7,*) '(',tt,') show'
	y1=y1+2.*tsize
	write(7,*) int(x1),int(y1),'  m'
	write(7,*) '/Helvetica findfont 300.00 scalefont setfont'
	write(tt,'(17hVelocity  Vectors)')
	write(7,*) '(',tt,') show s'
!C
	
END subroutine
!C
!C##################################################
SUBROUTINE SETCOL(NCOL)
!C##################################################
!C     This routine sets the values of R, G, and B
!C     composition of the colour table. NCOL values
!C     are set so that the range of values from MAX
!C     to MIN corresponds to the colour scale from
!C     pink to deep blue.
!C==================================================
	COMMON /RGB/ R(255),G(255),B(255)
!C
	IF(NCOL.LT.5) RETURN
!C
!C----------------------------------
!C.....FIVE RANGES IN SPECTRUM
!C----------------------------------
	NDC=NCOL/5
	DC=1./REAL(NDC)
!C----------------------------------
!C.....RANGE 1: PINK TO RED
!C----------------------------------
	DO L=1,NDC
		R(L)=1.
		G(L)=0.
		B(L)=1.-(L-1)*DC
	END DO
	NL=NDC
!C----------------------------------
!C.....RANGE 2: RED TO YELLOW
!C----------------------------------
	DO L=1,NDC
		R(NL+L)=1.
		G(NL+L)=(L-1)*DC
		B(NL+L)=0.
	END DO
	NL=NL+NDC
!C----------------------------------
!C.....RANGE 3: YELLOW TO GREEN
!C----------------------------------
	DO L=1,NDC
		R(NL+L)=1.-(L-1)*DC
		G(NL+L)=1.
		B(NL+L)=0.
	END DO
	NL=NL+NDC
!C----------------------------------
!C.....RANGE 4: GREEN TO LIGHT BLUE
!C----------------------------------
	DO L=1,NDC
		R(NL+L)=0.
		G(NL+L)=1.
		B(NL+L)=(L-1)*DC
	END DO
	NL=NL+NDC
!C
	NDC=NCOL-4*NDC
	DC=1./REAL(NDC)
!C
!C-------------------------------------
!C.....RANGE 5: LIGHT BLUE TO DARK BLUE
!C-------------------------------------
	DO L=1,NDC
		R(NL+L)=0.
		G(NL+L)=1.-(L-1)*DC
		B(NL+L)=1.
	END DO
!C
END subroutine
!C
!C#####################################################################
SUBROUTINE PROFIL(TITLE,FILN,LK)
!C#####################################################################
!C     This routine draws profiles of the variable FI at given X and
!C     Y cross-sections. The cross-sections are defined in input data
!C     as XPRO(I) and YPRO(I); ten profiles can be ploted. Up to 1000
!C     points can be in a profile (unlikely to be exceeded, unless
!C     very fine grid is used).
!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)
	COMMON /RGB/ R(255),G(255),B(255)
	DIMENSION YPR(1000),FIPR(1000),XPR(1000)
	CHARACTER TITLE*15,TT*17,FILN*4
!C
!C.....CALCULATE SCALING FACTOR, SCALE VARIABLE (accorning to AROMAX)
!C
	FIMAX=0.0
	DO IJ=1,NIJ
		PC=ABS(FI(IJ))
		FIMAX=MAX(FIMAX,PC)
	END DO
	SCF=AROMAX/(FIMAX+1.E-20)
	DO IJ=1,NIJ
		FI(IJ)=SCF*FI(IJ)
	END DO
!C
!C===============================================================
!C.....PROFILES AT GIVEN  X - COORDINATE
!C===============================================================
!C     We consider quadrilateral elements, made by CV centers;
!C     if the line X=XINT crosses two sides of such an element,
!C     a segment of the profile between the two crossing points
!C     is obtained (since nodes are connected by straight lines
!C     and linear interpolation is used, the same crossing points
!C     will be found in neighbor elements, so the profile will be
!C     continous). The variable values at the crossing points
!C     are obtained by interpolation from two corner nodes. The
!C     profile line is thick, the base line is thin.
!C
	IF (NPRX.NE.0) THEN
		CALL BPLOT(FILN,LK)
!C
		DO K=1,NPRX
			XINT=XPRO(K)
			DO I=1,NIM
				DO J=1,NJM
					IJ=LI(I)+J
					CALL PROFX(IJ,IP,XINT)
!C
!C     IP is the number of crossing points; it has to be two if the
!C     line crosses the element. 
!C
					IF (IP.EQ.2) THEN
						WRITE(7,*) '20 w ',INT(XINT+XP(1)),INT(YP(1)),' m ',	&
						&                    INT(XINT+XP(2)),INT(YP(2)),' l s'
						WRITE(7,*) '10 w ',INT(XINT),INT(YP(1)),' m ',			&
						&                    INT(XINT),INT(YP(2)),' l s'
					end if
				END DO
			END DO
!C
!C------------------------------------------------------------
!C.....COMPARISON WITH REFERENCE DATA
!C------------------------------------------------------------
!C     Reference data for each profile for which IEPRX=1 should