|
plotclaw.f.html |
 |
|
Source file: plotclaw.f
|
|
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/amrclaw/2d/lib
|
|
Converted: Mon Jan 21 2013 at 20:15:12
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
program main
c
c This is the graphics driver program for AMR, the Adaptive
c Mesh Refinement code. It reads commands and data from
c the file connected to unit 5, which should have been
c produced by AMR.
c Currently this program produces (1) an NCAR metafile
c which can be viewed on a tex4014 terminal or sent to the
c a high quality film processor, and (2) a log file indicating
c what has been produced on each frame of the graphics file.
c
c This program has the capability to plot:
c (1) The physical orientation of the grids at a given timestep
c (2) "Cluster" plots produced during error estimation
c (3) Contour plots of solution values and derived quantities
c at timesteps in which these values are available.
c A composite contour plot is produced for each selected
c variable using the best data available at each point
c in the physical domain. Other plots are available by
c interactively selecting them from a menu.
c
c At present, only the third capibility (contour plots) is
c on by default. The other two may be selected interactively
c by the "options" menu.
c Note: even when option 1 is turned off, a frame containing
c the physical orientation of the grids within the
c problem domain is produced before each series of
c contour plots.
c
c When the program is run, several menus are presented to the
c user. The first asks which variables and derived values (such
c as Density or Log of Entropy) should be plotted. The second menu
c lists some plotting options. These include the two options
c listed above plus (3) an option to produce a sequence of
c contour plots, one for each indicidual grid, (4) an option
c to show the orientation of the individual grids superimposed
c on the composite plot, (5) an option to produce a composite
c plot for each level, using only the data available on that
c level, and (6) an option to produce a composite plot of
c all levels greater than 1.
c
c Finally, the user is prompted for a threshold time, before
c which no plots are produced.
c
c The menu of variables to plot can be extended and modified
c by changing the the definitions of "numvars", "names",
c and "namelen" defined in the common blocks "userdef" and
c "userdef2" below.
c The new variables can be computed from the given variables
c in subroutine "derive". See also subroutine "getinput"
c for how the menu is produced, although this should not have
c to be changed.
include "call.i"
c
character*30 timestmp
character*10 command
character*50 cline, line, bline
character*3 cmd, cluster, drawtree, plotgrid
character*8 str
integer clen, framecnt, numptc(maxcl)
integer oldmode
logical gprint2, plotclus, plottree, echo
logical allgrids, levplots, fineplot, plotbox
logical nocomp, lineplot
c
c
c :::::::: User defined variable definitions
c :::::::: "numvars" is the number of variables that
c :::::::: can be viewed.
c :::::::: "name" is a character array containing the
c :::::::: names of these variables
c :::::::: "namelen" is an array of the lengths of the above
c :::::::: variable names
c :::::::: "noplot" if a logical array set in "getinput" that
c :::::::: determines whether the corresponding
c :::::::: variable has been selected for plotting
c
integer numvars, namelen(15), blen, blinelen
character name(15)*15, bname*15
logical noplot(15)
common /userdef/ numvars, namelen, noplot
common /userdef2/ name
data numvars/10/
data name/"DENSITY ","X VELOCITY ","Y VELOCITY ",
1 "TOTAL VELOCITY ","TOTAL ENERGY ","PRESSURE ",
1 "ENTROPY ","MACH NUMBER ","ENTHALPY ",
1 "LOG DENSITY "," "," ",
1 " "," "," "/
data namelen/7,10,10,14,12,8,7,11,8,11,0,0,0,0,0/
data bname/"BODYPLOT "/
data blen/8/
data gprint2/.true./
data echo/.false./
c
c *****************************************************************
c graphics driver - reads in file generated by amr output routines:
c and generates the requested graphics output
c the commands are:
c basic time - what follows is lfine,levellist, and grid
c info. from levels lst to lend
c vals mgrid, nvar - soln values follow
c if (mgrid <> 0) only grid mgrid values else entire
c level from lst to lend
c nvar = number of variables to be read for each grid point
c clusters nclust, lpar - cluster info. follows
c nclust = num. of clusters, lpar = level being examined
c in each option, more info. if read in according to that option.
c *****************************************************************
c
c
oldmode = ieee_handler("set","common",SIGFPE_ABORT)
if (oldmode .ne. 0) then
write(6,*) ' did not set floatingpt. traps right '
stop
endif
open(3,file='plotfile',status='old',form='formatted')
open(16,file='ploterr',status='unknown',form='formatted')
open(7,file='plotlog',status='unknown',form='formatted')
call opngks
call ctable
call initspac
c ::::::::: Get variable list and time threshold
call getinput( thresh, allgrids, levplots, fineplot, plotbox,
1 plotclus, plottree, nocomp ,lineplot)
c ::::::::: Initialize command patterns
cluster = "CLU"
drawtree = "DTR"
plotgrid = "VAL"
framecnt = 0
c ::::::::: Read time stamp
10 str = " "
timestmp = "TIME = "
lentime = 7
read(3,101,end=99) str, time
101 format( a8, f10.5 )
if (index(str,'*TIME') .le. 0) then
write(6,*) "Illegal time stamp: ",str
call clsgks
stop
endif
str = " "
write(str,111) time
111 format( f8.3 )
c ::::::: skip leading blanks in time string
do 20 K = 1, len(str)
if (str(K:K) .ne. ' ') then
lentime = lentime + 1
timestmp(lentime:lentime) = str(K:K)
endif
20 continue
c ::::::::: Now read in the graph structure
c ::::::::: NOTE: The graph structure is dumped by AMR
c ::::::::: after every timestamp is produced.
call basic(lst,lend)
c ::::::::: Now read the command
command = " "
cmd = " "
read(3,102,end=99) command, ivar1, ivar2, ivar3
102 format( a10, 3i10 )
if (echo) then
write(16,103) str, command, ivar1, ivar2, ivar3, lst, lend
103 format(' Time Stamp: ',a15,/,' Input Command: ',a10,3i10,
* /,' lst, lend: ',2i10 )
endif
I = index(command,'*')
if (I .le. 0) then
write(6,*) "Illegal Command Line: ",command
call clsgks
stop
endif
cmd = command(I+1:I+3)
c ::::::::: Now process the command
if (cmd .eq. plotgrid) then
c read the grid values and plot contours
nvar = ivar1
call invals(lst,lend,nvar,thresh)
if (time .ge. thresh) then
if (plottree) then
c :::::::: Draw the grid overlap structure
c :::::::: before drawing the variables
c :::::::: since plottree is turned off
call dtree(1,lfine,timestmp,lentime,gprint2)
framecnt = framecnt + 1
write(7,299) framecnt, "DRAWTREE", timestmp
call frame
endif
c :::::::: Produce a series of plots for each
c :::::::: variable selected.
do 50 ivar = 1, numvars
if (noplot(ivar)) goto 50
line = name(ivar)(1:namelen(ivar))// ", " //timestmp
bline = bname(1:blen)// ", " //timestmp
linelen = lentime + namelen(ivar) + 2
blinelen = lentime + blen + 2
c :::::::: cycle through grids computing quantities
c :::::::: to be plotted.
level = lst
150 if (level .le. lend) then
mptr = lstart(level)
160 if (mptr .ne. 0) then
c nrows = node(maxnumrow,mptr) + 1
c ncols = node(maxnumcol,mptr) + 1
nrows = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncols = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
locvar = node(tempptr,mptr)
loc = igetsp(nrows*ncols)
node(store1,mptr) = loc
call derive(alloc(locvar),alloc(loc),
1 nrows,ncols,nvar,ivar)
mptr = node(levelptr,mptr)
goto 160
endif
level = level + 1
goto 150
endif
c ::::::: plot all grids for variable ivar
if (lineplot) then
call lplotx(framecnt,line,linelen)
else
call dvals(lst,lend,line,linelen,nvar,
1 ivar,framecnt,allgrids,levplots,fineplot,
2 plotbox,bline,blinelen,nocomp)
endif
50 continue
endif
c ::::::::: Reclaim space used by each grid
level = lst
60 if (level .le. lend) then
mptr = lstart(level)
70 if (mptr .ne. 0) then
c nrows = node(maxnumrow,mptr) + 1
c ncols = node(maxnumcol,mptr) + 1
nrows = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncols = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
locvar = node(tempptr,mptr)
call reclam( locvar, nrows*ncols*nvar )
node(tempptr,mptr) = 0
mptr = node(levelptr,mptr)
go to 70
endif
level = level + 1
goto 60
endif
c :::::::: read the next command
go to 10
endif
c
if (cmd .eq. cluster) then
c ::::::: plot the clusters
nclust = ivar1
lpar = ivar2
nxypts = ivar3
if (nclust .gt. maxcl) then
write(16,104) maxcl, nclust
104 format(" Cannot read more than ",i3,
1 " clusters: nclust = ",i3)
go to 99
endif
locbad = igetsp( 2*nxypts )
call inclus(nclust,numptc,nxypts,alloc(locbad),mkid,echo)
if ( plotclus .and. ( time .ge. thresh)) then
if (nclust .gt. maxcl) then
write(16,109) nclust
109 format( " Cannot plot so many clusters:",
1 " nclust = ",i4,/,
2 " ... Skipping to next command")
write(7,119) nclust
119 format( " ... Skipping cluster plot: nclust = ",
1 i4," > maxcl")
else
cline = timestmp//" "
clen = lentime + 10
cline(lentime+1:clen) = ", (BEFORE)"
call dclust(nclust,numptc,nxypts,alloc(locbad),
1 lpar,lstart(lpar+1),cline,clen,.false.)
framecnt = framecnt + 1
write(7,199) framecnt, "CLUSTER ", cline
call frame
cline(lentime+1:clen) = ", (AFTER) "
call dclust(nclust,numptc,nxypts,alloc(locbad),
1 lpar,mkid,cline,clen,.true.)
framecnt = framecnt + 1
write(7,199) framecnt, "CLUSTER ", cline
call frame
199 format( "Frame# ",i3,2x,a8,2x,a50)
endif
endif
call reclam( locbad, 2*nxypts )
c read next command
go to 10
endif
c
if (cmd .eq. drawtree) then
c draw the overlaping rectangles to show grid arrangement
if (plottree .and. (time .ge. thresh)) then
call dtree(1,lfine,timestmp,lentime,gprint2)
framecnt = framecnt + 1
write(7,299) framecnt, "DRAWTREE", timestmp
299 format( "Frame# ",i3,2x,a8,2x,a30 )
call frame
endif
go to 10
endif
c
c invalid command read, issue error message and try again
write(16,106) command
106 format( "Unknown Command: ",a10," Try again..." )
go to 10
c
99 call clsgks
stop
end
c
c ------------------------------------------------------------------
c
subroutine getinput(thresh, allgrids, levplots, fineplot,
1 plotbox, plotclus, plottree,
2 nocomp, lineplot )
character varlist*30, temp*20, str*30, str2*10
character ch, lowch, upch, ych
logical allgrids, levplots, fineplot, plotbox
logical plotclus, plottree, nocomp, lineplot
real thresh
integer numvars, namelen(15)
character name(15)*15
logical noplot(15)
common /userdef/ numvars, namelen, noplot
common /userdef2/ name
data ych/'y'/
do 5 I = 1, 15
noplot(I) = .true.
5 continue
varlist = " "
do 10 I = 1, numvars
ch = char( ichar( 'A' ) + I - 1 )
write(6,900) ch, name(I)
10 continue
900 format( "Variable Symbol ",a1," = ",a15 )
write(6,*) "Enter symbols of variables to be plotted"
read(5,910) varlist
910 format( a30 )
c ::::::: now parse varlist to find variable indices
c ::::::: both lower and upper case inputs are accepted
ilen = 0
temp = " "
do 20 I = 1, numvars
lowch = char( ichar( 'a' ) + I - 1 )
upch = char( ichar( 'A' ) + I - 1 )
ixch = index(varlist,lowch) + index(varlist,upch)
if (ixch .gt. 0) then
noplot(I) = .false.
ilen = ilen + 1
temp(ilen:ilen) = upch
endif
20 continue
write(6,950)
950 format(
1 / "A composite contour plot will be generated by default",
2 /,"for each of these variables.",
3 //,"In addition, here is a list of plotting options:",
4 /," S = do not draw a composite contour plot",
6 /," U = draw grid tree on EACH available timestep",
7 /," V = show CLUSTER plots on EACH available timestep",
8 /," W = produce a contour plot of each grid separately",
9 /," X = draw boxes around individual grids in composite",
1 /," Y = produce a contour plot at each level",
2 /," Z = produce a contour plot of (2) finest levels",
3 /," L = produce a line plot through domain in x",
4 /,"Enter symbols of selected options (if any):" )
varlist = " "
read(5,910) varlist
ix = index(varlist,'S') + index(varlist,'s')
nocomp = (ix .gt. 0)
ix = index(varlist,'U') + index(varlist,'u')
plottree = (ix .gt. 0)
ix = index(varlist,'V') + index(varlist,'v')
plotclus = (ix .gt. 0)
ix = index(varlist,'W') + index(varlist,'w')
allgrids = (ix .gt. 0)
ix = index(varlist,'X') + index(varlist,'x')
plotbox = (ix .gt. 0)
ix = index(varlist,'Y') + index(varlist,'y')
levplots = (ix .gt. 0)
ix = index(varlist,'Z') + index(varlist,'z')
fineplot = (ix .gt. 0)
ix = index(varlist,'L') + index(varlist,'l')
lineplot = (ix .gt. 0)
c
30 str = " "
str2 = " "
write(6,920)
920 format(//,' Enter minimum timestep to be plotted')
read(5,110) str(10:30)
110 format( a20 )
ix = index(str,'.')
if (ix .le. 0) then
write(6,*) "ERROR: Real number must have a decimal point"
goto 30
endif
str2(1:10) = str(ix-4:ix+5)
read(str2,120,err=999) thresh
120 format( f10.5 )
c write(6,909) temp(1:ilen), thresh
909 format("Selected Variables = ",a,", Minumum Timestep = ",f10.5)
return
999 write(6,*) "Error in reading TIME threshold"
goto 30
end
c
c ----------------------------------------------------------------
c
c ::::::: This subroutine is used to derive the values
c ::::::: of the variables to be plotted from the values
c ::::::: read in by subroutine "invars"
c
subroutine derive(val,derval,nrows,ncols,nvar,ivar)
dimension val(nrows,ncols,nvar), derval(nrows,ncols)
integer nrows, ncols, nvar, ivar
data gamma/1.4/
data spval2/-1.0e10/, spval/1.0e10/
gamm1 = gamma - 1.0
if (ivar .eq. 1 .or. ivar .eq. 10) then
c ::::::: Return Density (or Log density)
if (ivar .eq. 1) then
do 100 J = 1, ncols
do 110 I = 1, nrows
derval(I,J) = val(I,J,1)
110 continue
100 continue
else
do 101 J = 1, ncols
do 111 I = 1, nrows
if (val(I,J,1) .eq. spval) then
derval(I,J) = spval2
else
derval(I,J) = alog(val(I,J,1))
endif
111 continue
101 continue
endif
elseif (ivar .le. 3) then
c ::::::: return Velocity
do 200 J = 1, ncols
do 210 I = 1, nrows
if (abs(val(I,J,ivar)) .eq. spval) then
derval(I,J) = spval2
else
derval(I,J) = val(I,J,ivar)/val(I,J,1)
endif
210 continue
200 continue
elseif ( ivar .eq. 4) then
c ::::::: Return Total Velocity
do 300 J = 1, ncols
do 310 I = 1, nrows
if (abs(val(I,J,1)) .eq. spval) then
derval(I,J) = spval2
else
derval(I,J) = (val(I,J,2)**2 + val(I,J,3)**2)
1 /(val(I,J,1)**2)
endif
310 continue
300 continue
elseif ( ivar .eq. 5 ) then
c ::::::: Return Energy per unit volume
do 400 J = 1, ncols
do 410 I = 1, nrows
if (abs(val(I,J,1)) .eq. spval) then
derval(I,J) = spval2
else
derval(I,J) = val(I,J,4)/val(I,J,1)
endif
410 continue
400 continue
elseif ( ivar .le. 7 .or. ivar .eq. 9) then
c ::::::: Return Pressure or Entropy or Enthalpy...
c ::::::: First compute Pressure
do 500 J = 1, ncols
do 510 I = 1, nrows
if (abs(val(I,J,1)) .eq. spval) then
derval(I,J) = spval2
else
derval(I,J) = gamm1*(val(I,J,4) - (val(I,J,2)**2 +
1 val(I,J,3)**2)/(2.0*val(I,J,1)))
endif
510 continue
500 continue
if ( ivar .eq. 7 ) then
c ::::::: Compute Entropy using computed pressure
do 600 J = 1, ncols
do 610 I = 1, nrows
if (abs(val(I,J,1)) .ne. spval) then
derval(I,J) = alog(derval(I,J)) -
1 gamma*alog(val(I,J,1))
endif
610 continue
600 continue
endif
if ( ivar .eq. 9 ) then
c ::::::: Compute Enthalpy
do 603 J = 1, ncols
do 613 I = 1, nrows
if (abs(val(I,J,1)) .ne. spval) then
derval(I,J) = (derval(I,J)+val(I,J,4))/(val(I,J,1))
endif
613 continue
603 continue
endif
elseif (ivar .eq. 8) then
c ::::::: Return Mach Number
do 700 J = 1, ncols
do 710 I = 1, nrows
if (abs(val(I,J,1)) .eq. spval) then
derval(I,J) = spval2
else
vel2 = (val(I,J,2)**2 + val(I,J,3)**2)
1 /(val(I,J,1)**2)
p = gamm1*(val(I,J,4)-(val(I,J,2)**2 +
1 val(I,J,3)**2)/(2.0*val(I,J,1)))
c2 = gamma*p / val(I,J,1)
derval(I,J) = sqrt(vel2/c2)
endif
710 continue
700 continue
c
endif
return
end
subroutine initspac
include "call.i"
c
c ::::::::: Initialize the data space "alloc" and its free list "lfree"
c ::::::::: to be one contiguous block from which chunks will be
c ::::::::: allocated on a first fit basis by routine "igetsp".
c
c
c ::::::::: clear data space
do 10 i = 1, memsize
alloc(i) = 0.0
10 continue
c
c ::::::::: initialize the free list so that:
c ::::::::: Location 1 is a dummy header node
c ::::::::: Location 2 starts at position 1 and contains memsize elements
c ::::::::: Location 3 starts at position 645000 and contains 0 elements
c ::::::::: All remaining loctions are marked as empty
lenf = 3
idimf = lfdim
do 20 i = 1, idimf
lfree( i, 1 ) = 0
lfree( i, 2 ) = 0
20 continue
lfree( 2, 1 ) = 1
lfree( 2, 2 ) = memsize
lfree( 3, 1 ) = memsize + 2
c
return
end
c
c ------------------------------------------------------------------
c
function igetsp (nwords)
c
include "call.i"
logical sprint2
data sprint2/.false./
c
c allocate contiguous space of length nword in main storage array
c alloc. user code (or pointer to the owner of this storage)
c is mptr. lenf = current length of lfree list.
c
c
c find first fit from free space list
c
itake = 0
do 20 i = 1, lenf
if (lfree(i,2) .ge. nwords) then
itake = i
goto 25
endif
20 continue
c ::::::::: if we got here, there is not enough memory left :::::
c ::::::::: to satisfy the request. :::::
write(16,901) nwords
901 format(' require ',i10,' words - either none left or not big',
1 ' enough space')
call clsgks
stop
c
c ::::::::: request is satisfied, now allocate memory :::::
25 left = lfree(itake,2) - nwords
igetsp = lfree(itake,1)
c
if (left .gt. 0) then
c ::::::::: retain the unused part of this memory :::::
lfree(itake,2) = left
lfree(itake,1) = lfree(itake,1) + nwords
else
c ::::::::: item is totally removed, move others in the list :::::
c ::::::::: up one space :::::
lenf = lenf - 1
do 40 i = itake, lenf
lfree(i,1) = lfree(i+1,1)
lfree(i,2) = lfree(i+1,2)
40 continue
endif
c
if (sprint2) write(16,100) nwords, igetsp
100 format(' allocating ',i8,' words in location ',i8)
c
return
end
subroutine reclam (index, nwords)
c
c return of space. add to free list.
c iplace points to next item on free list with larger index than
c the item reclaiming, unless said item is greater then
c everything on the list.
c
include "call.i"
logical sprint2
data sprint2/.false./
c
c
c find position in free list
c
do 20 i = 1, lenf
iplace = i
if (lfree(i,1) .gt. index) goto 30
20 continue
c
c no position found, error
c (?? what is the last piece of memory is returned ??)
c
write(16,902)
902 format(' no insertion pointer into freelist. error stop')
call clsgks
stop
c
c try to merge this memory back into the free list
c
30 iprev = iplace - 1
if ((iprev .ne. 0) .and.
1 ((lfree(iprev,1) + lfree(iprev,2)) .eq. index) ) then
c
c merge with previous segment
c
lfree(iprev,2) = lfree(iprev,2) + nwords
c
c try to merge with following segment as well
c
if (lfree(iplace,1) .eq. (index + nwords)) then
lfree(iprev,2) = lfree(iprev,2) + lfree(iplace,2)
ipp1 = iplace + 1
do 60 i = ipp1, lenf
lfree(i-1,1) = lfree(i,1)
lfree(i-1,2) = lfree(i,2)
60 continue
lenf = lenf - 1
endif
c
c no previous merge, try to merge with next segment
c
elseif (lfree(iplace,1) .eq. (index + nwords)) then
lfree(iplace,1) = index
lfree(iplace,2) = lfree(iplace,2) + nwords
c
c no merges at all, try to insert new segment in the
c free list by bumping subsequent segments
c
elseif (lenf .ne. idimf) then
do 80 ii = iplace, lenf
i = lenf + 1 - ii + iplace
lfree(i,1) = lfree(i-1,1)
lfree(i,2) = lfree(i-1,2)
80 continue
lenf = lenf + 1
lfree(iplace,1) = index
lfree(iplace,2) = nwords
c
c no room in free list to add segment. error.
c
else
write(16,901) idimf
901 format(' free list full with ',i5,' items')
call clsgks
stop
endif
if (sprint2) write(16,100) nwords, index
100 format(' reclaiming ',i8,' words at loc. ',i8)
return
end
subroutine basic (lst,lend)
c
include "call.i"
common /userdt/ pinf,rhoinf,vinf,chord,alpha,foil,xprob,yprob
c
c read in basic info. generated for every graphics call
c
read(3,102) lfine, (lstart(i),i=1,lfine), lwidth
102 format(10i6)
read(3,105) xupper,yupper,xlower,ylower
105 format(4e15.8)
read(3,102) lst, lend
read(3,106)(hxposs(i),i=1,lfine)
read(3,106)(hyposs(i),i=1,lfine)
106 format(6e13.8)
c
xprob = xupper
yprob = yupper
level = lst
10 if (level .le. lend) then
mptr = lstart(level)
20 if (mptr .ne. 0) then
read(3,103) mptr,(node(i,mptr),i=1,nsize)
103 format(10i7)
read(3,104) (rnode(i,mptr),i=1,rsize)
104 format(5e14.8)
c # here as a reminder
c node(maxnumrow,mptr) = node(ndihi,mptr)-node(ndilo,mptr)+3
c node(maxnumcol,mptr) = node(ndjhi,mptr)-node(ndjlo,mptr)+3
mptr = node(levelptr,mptr)
go to 20
endif
30 level = level + 1
go to 10
endif
c
return
end
subroutine inclus(nclust,numptc,nxypts,badpts,mkid,echo)
c
include "call.i"
dimension numptc(nclust), badpts(2,nxypts)
logical echo
c
c read in info. to plot grids, flagged pts., clusters, new grids.
c
read(3,105) iregsz,jregsz
105 format(6i10)
read(3,100) (numptc(i),i=1,nclust)
100 format(10i6)
if (echo) then
write(16,102) nclust, (numptc(i),i=1,nclust)
102 format(' no. of clusters ',i5, ' with no. of points ',
1 4(/,10x,15i5))
endif
read(3,101)((badpts(i,j),i=1,2),j=1,nxypts)
101 format(4e12.4)
c
c there is always at least one subgrid, since the output routine
c is only called if there is at least one flagged pt.
mkid = 0
20 continue
read(3,103) minput, (node(i,minput),i = 1,nsize)
103 format(10i6)
if (mkid .eq. 0) mkid = minput
read(3,104) (rnode(i,minput), i = 1,rsize)
104 format(5e12.4)
minput = node(levelptr, minput)
if (minput .ne. 0) go to 20
c
return
end
subroutine dclust(nclust,numptc,nxypts,badpts,lpar,mkid,
1 cline,clen,second)
c
c draw existing grids level lpar, flagged points in their clusters,
c and fine subgrids starting with mkid.
c mkid may be start of old or new (i.e. not yet in place in lstart)
c level of subgrids.
c
include "call.i"
character*1 isign(50), usesign
character*50 cline
integer clen
logical gprint2, second
c
dimension numptc(nclust), badpts(2,nxypts)
data isign/'A','B','C','D','E','F','G','H',
1 'I','J','K','L','M','N','O','P',
2 'Q','R','S','T','U','V','W','X',
3 'Y','Z','1','2','3','4','5','6',
4 '7','8','9','0','+','-','*', '/',
5 '<','>','?','!','"','#','$','%',
6 '&',':'/
data gprint2/.false./
c
call dtree(lpar,lpar,cline,clen,gprint2)
c
ist = 0
do 10 icl = 1, nclust
npts = numptc(icl)
do 20 i = 1, npts
xpt = xlower + badpts(1,ist+i)*hxposs(lpar)
ypt = ylower + badpts(2,ist+i)*hyposs(lpar)
if (second) then
usesign = '+'
else
usesign = isign(icl)
endif
call pwrit(xpt,ypt,usesign,1,0,0,0)
20 continue
ist = ist + numptc(icl)
10 continue
c
c output new fine subgrids
c
mnew = mkid
30 if (mnew .ne. 0) then
call drect(mnew,gprint2)
mnew = node(levelptr,mnew)
go to 30
endif
c
return
end
subroutine invals (lst,lend,nvar,thresh)
c
include "call.i"
integer charst(10)
data lstar/1h*/,lexs/1hS/
c
spval = 1.0e10
c
c read in soln values for val option of graphics package
c
level = lst
10 if (level .le. lend) then
mptr = lstart(level)
20 if (mptr .ne. 0) then
read(3,100) charst,minput
100 format(10a1,i10)
if (charst(1).ne.lstar .or. charst(2).ne.lexs .or.
* minput .ne. mptr) then
write(16,901) charst,mptr,minput,lst,lend
901 format(1h ,10a1,/,1x,4i10)
call clsgks
stop
endif
c ::::::::: Data dumped in an array (2..maxrows,2..maxcols)
c ::::::::: need to allocate an array (1..maxrows+1,1..maxcols+1)
c ::::::::: to add possible boundary values
c nrows = node( maxnumrow,mptr ) + 1
c ncols = node( maxnumcol,mptr ) + 1
nrows = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncols = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
num = nrows*ncols
ist = igetsp(num*nvar)
node(tempptr,mptr) = ist
call readvals(alloc(ist),nrows,ncols,nvar)
mptr = node(levelptr,mptr)
go to 20
endif
level = level + 1
go to 10
endif
c
c add boundary values for each grid, all variables, for smooth
c contour plotting at grid interfaces
c
c don't bother if below threshold time for plotting
c
time = rnode(timemult,lstart(lst))
if (time .lt. thresh) go to 99
c
c find corners
c
xlo = 1.e23
ylo = 1.e23
xhi = -1.e23
yhi = -1.e23
mptr = lstart(1)
40 if (mptr .ne. 0) then
xhi = amax1(xhi,rnode(cornxhi,mptr))
yhi = amax1(yhi,rnode(cornyhi,mptr))
xlo = amin1(xlo,rnode(cornxlo,mptr))
ylo = amin1(ylo,rnode(cornylo,mptr))
mptr = node(levelptr,mptr)
go to 40
endif
c
level = lst
110 if (level .le. lend) then
hx = hxposs(level)
hy = hyposs(level)
mptr = lstart(level)
120 if (mptr .ne. 0) then
c nrow = node(maxnumrow,mptr) + 1
c ncol = node(maxnumcol,mptr) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
call addbndry(level,mptr,alloc(node(tempptr,mptr)),
. nrow,ncol,nvar,xlo,ylo,xhi,yhi,spval,lst,
. hx,hy)
mptr = node(levelptr,mptr)
go to 120
endif
level = level + 1
go to 110
endif
c
c set special value to signal exterior cells within domain
c
99 return
end
c
c -------------------------------------------------------
c
subroutine readvals(a,nrows,ncols,nvar)
dimension a(nrows,ncols,nvar)
do 10 ivar = 1, nvar
read(3,100) ((a(i,j,ivar),i=2,nrows-1),j=2,ncols-1)
10 continue
100 format(5e12.6)
99 return
end
c
c --------------------------------------------------------
c
subroutine dvals(lst,lend,line,linelen,nvar,ivar,
1 framecnt,allgrids,levplots,fineplot,
2 plotbox,bline,blinelen,nocomp)
c
include "call.i"
character line*50, header*50, strtmp*5, bline*50
integer headlen, hlen, framecnt, blinelen
logical nocomp
logical allgrids, levplots, fineplot, plotbox
common /conre4/ sizel, sizem, sizep, nrep, ncrt, ilab,
1 nulbll, ioffd, ext, ioffm, isolid, nla,
1 nlm, xlt, ybt, side
c
c ******************************************************************
c dvals - generate a series of surface, contour, or body plots of every
c grid in a subtree,
c ******************************************************************
c
c
c ::::::::: set NCAR parameter to turn off printing of contour labels :::::
ilab = 0
c
if (.not. allgrids) goto 500
c
c ::::::::: Set up header for contour plots
header = line
header(linelen+1:linelen+7) = ", GRID "
headlen = linelen + 7
c ::::::::: walk through levels from coursest to finest generating a :::::
c ::::::::: separate contour plot for each grid at each level. :::::
level = lst
10 if (level .le. lend) then
mptr = lstart(level)
20 if (mptr .ne. 0) then
c nrow = node(maxnumrow,mptr) + 1
c ncol = node(maxnumcol,mptr) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
loc = node(store1,mptr)
header(headlen:headlen+5) = " "
hlen = headlen
write(strtmp,111) mptr
111 format( i5 )
do 25 K = 1, 5
if (strtmp(K:K) .ne. ' ') then
hlen = hlen + 1
header(hlen:hlen) = strtmp(K:K)
endif
25 continue
c # work with copy of array, don't mess up original.
c # due to setting spval, then too much restored later.
num = nrow*ncol
lcopy = igetsp(num)
do 26 i = 0, num-1
alloc(lcopy+i) = alloc(loc+i)
26 continue
call contor(alloc(lcopy),nrow,ncol,header,hlen,mptr)
call reclam(lcopy,num)
framecnt = framecnt + 1
write(7,199) framecnt, "GRIDPLOT", header
199 format( "Frame# ",i3,2x,a8,2x,a50 )
mptr = node(levelptr,mptr)
goto 20
endif
level = level + 1
go to 10
endif
c
c ::::::::: Produce a composite contour plot in which the grids are :::::
c ::::::::: plotted from finest to coarsest and in such a way that :::::
c ::::::::: coarse grids are "blacked out" in regions where fine :::::
c ::::::::: grids exist. :::::
500 call composit(lst,lend,nvar,ivar,line,linelen,framecnt,
1 levplots, fineplot, plotbox,
2 bline,blinelen,nocomp)
c
return
end
c
c ------------------------------------------------------------
c
subroutine composit(lst, lend, nvar, ivar, line, linelen,
1 framecnt, levplots, fineplot, plotbox,
2 bline, blinelen, nocomp)
c
include "call.i"
character line*50, cline*50, lline*50, fline*50, bline*50
real xmin(maxlv),xmax(maxlv),ymin(maxlv),ymax(maxlv)
integer clen, framecnt, flen, llen, blinelen, blen
logical levplots, fineplot, plotbox, nocomp
common /conre1/ ioffp, spval
common /conre4/ sizel, sizem, sizep, nrep, ncrt, ilab,
1 nulbll, ioffd, ext, ioffm, isolid, nla,
2 nlm, xlt, ybt, side
data ncontor /30/
c ::::::::: set NCAR parameter to turn off printing of contour labels :::::
ilab = 0
nlm = 40
call cpseti('LLP - line label position',0)
call cpsetc('HLT - high/low label text',' ')
call cpsetc('ILT - informational label text',' ')
call cpseti('NCL - number of contour levels',ncontor)
call cpseti('SET - do set call flag',0)
call cpsetr('SPV - set special value ',spval)
c ::::::::: set NCAR parameter to blackout grid points set to spval :::::
ioffp = 1
spval = 1.0e10
spval2 = -1.0e10
c ::::::::: set up header for composite plots
cline(1:50) = line(1:50)
lline(1:50) = line(1:50)
fline(1:50) = line(1:50)
clen = linelen + 11
llen = linelen + 9
flen = linelen + 12
blen = blinelen + 11
cline(linelen+1:clen) = ", COMPOSITE"
lline(linelen+1:llen) = ", LEVEL "
fline(linelen+1:flen) = ", FINE GRIDS"
bline(blinelen+1:blen) = ", COMPOSITE"
c :::::::: compute min and max values, level by level :::::
totmin = 1.0e+10
totmax = -1.0e+10
level = lst
100 if (level .le. lend) then
mptr = lstart(level)
xmin(level) = 1.0e+10
ymin(level) = 1.0e+10
xmax(level) = -1.e+10
ymax(level) = -1.e+10
105 if (mptr .ne. 0) then
c nrows = node(maxnumrow,mptr) + 1
c ncols = node(maxnumcol,mptr) + 1
nrows = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncols = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
loc = node(store1,mptr)
call mnmax(alloc(loc),nrows,ncols,vmin,vmax,spval2)
totmin = amin1( totmin, vmin )
totmax = amax1( totmax, vmax )
xmin(level) = amin1(xmin(level),rnode(cornxlo,mptr))
ymin(level) = amin1(ymin(level),rnode(cornylo,mptr))
xmax(level) = amax1(xmax(level),rnode(cornxhi,mptr))
ymax(level) = amax1(ymax(level),rnode(cornyhi,mptr))
mptr = node(levelptr,mptr)
goto 105
endif
level = level + 1
goto 100
endif
c ::::::: statistics on the lowest level grids
xmin1 = xmin(lst)
ymin1 = ymin(lst)
xmax1 = xmax(lst)
ymax1 = ymax(lst)
lev1 = lstart(lst)
dx = hxposs(lev1)/100.0
dy = hyposs(lev1)/100.0
xlo = xmin1 - dx
ylo = ymin1 - dy
xhi = xmax1 + dx
yhi = ymax1 + dy
c ::::::: set up parameters for contour increments
cinc = (totmax - totmin) / float(ncontor)
cmin = totmin + .5*cinc
cmax = totmax - .5*cinc
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c ::::::::::::::: produce composite plots of each level :::::::::::::::
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
if (levplots) then
level = lst
300 if (level .le. lend) then
mptr = lstart(level)
delx = xmax(level) - xmin(level)
dely = ymax(level) - ymin(level)
xlen = .8*amin1( 1.0, delx/dely )
ylen = .8*amin1( 1.0, dely/delx )
x1set = .1
x2set = .1+xlen
y1set = .1
y2set = .1+ylen
call set(.1,.1+xlen,.1,.1+ylen,xmin(level),xmax(level),
1 ymin(level), ymax(level), 1)
dx = hxposs(level)/2.0
dy = hyposs(level)/2.0
320 if (mptr .ne. 0) then
loc = node(store1,mptr)
c nrow = node(maxnumrow,mptr) + 1
c ncol = node(maxnumcol,mptr) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
x1 = rnode(cornxlo,mptr)
y1 = rnode(cornylo,mptr)
x3 = rnode(cornxhi,mptr)
y3 = rnode(cornyhi,mptr)
call drawbox(x1,y1,x3,y3,level)
x1 = x1 - dx
y1 = y1 - dy
x3 = x3 + dx
y3 = y3 + dy
x1ratio = (x1-xmin(level))/(xmax(level)-xmin(level))
x3ratio = (x3-xmin(level))/(xmax(level)-xmin(level))
y1ratio = (y1-ymin(level))/(ymax(level)-ymin(level))
y3ratio = (y3-ymin(level))/(ymax(level)-ymin(level))
call set(x1set+x1ratio*(x2set-x1set),
. x1set+x3ratio*(x2set-x1set),
. y1set+y1ratio*(y2set-y1set),
. y1set+y3ratio*(y2set-y1set),
. 1.5,float(nrow)-.5,1.5,float(ncol)-.5,1)
c ::::::::: plot the contour lines :::::
call setspv(alloc(loc),nrow,ncol,spval)
call conrec( alloc(loc), nrow, nrow, ncol,
1 cmin, cmax, cinc, -1, -1, 0 )
c call cpcnrc( alloc(loc), nrow, nrow, ncol,
c 1 cmin, cmax, cinc, 1, -1, 0 )
c ::::::::: reset the window to that of the coarsest grid :::::
call set(.1,.1+xlen,.1,.1+ylen,xmin(level),xmax(level),
1 ymin(level), ymax(level), 1)
mptr = node(levelptr,mptr)
goto 320
endif
lline(llen:llen) = char( ichar('0') + level )
call set(0.,1.,0.,1.,xmin(level),xmax(level),
1 ymin(level),ymin(level)+1.2*(ymax(level)-ymin(level)),1)
xplace = xmin(level) + (xmax(level)-xmin(level))/4.
yplace = ymin(level)+ (ymax(level)-ymin(level))*1.1
call pwrit(xplace,yplace,lline,llen,2,0,-1)
call gstxci(1)
framecnt = framecnt + 1
write(7,199) framecnt, "GRIDPLOT", lline
call frame
level = level + 1
goto 300
endif
endif
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c :::::::::::::::::: BLACKOUT OVERLAP BY FINE GRIDS :::::::::::::::::::
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
level = lst
130 if (level .le. lend) then
mptr = lstart(level)
140 if (mptr .ne. 0) then
c nrow = node(maxnumrow,mptr) + 1
c ncol = node(maxnumcol,mptr) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
num = nrow*ncol
lorg = node(store1,mptr)
lblk = igetsp(num)
do 125 i = 0, num-1
alloc(lblk+i) = alloc(lorg+i)
125 continue
if (level .ne. lend) then
call blackout(mptr,level,alloc(lblk),nrow,ncol,spval)
c write(6,502) mptr, nrow, ncol
502 format("Blackout: mptr, nrow, ncol = ",3(i4,3x) )
endif
c :::::::: create space for the final grid :::::
lrst = igetsp(num)
node(store2,mptr) = lrst
c :::::::: restore boarders of blackened areas :::::
call restore(alloc(lorg),alloc(lblk),alloc(lrst),nrow,
1 ncol,mptr,spval,xlo,ylo,xhi,yhi)
c :::::::: reclaim space used for blackened array :::::
call reclam(lblk,num)
mptr = node(levelptr,mptr)
goto 140
endif
level = level + 1
goto 130
endif
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c ::::::::::: PRODUCE A COMPOSITE PLOT OF ALL LEVELS ::::::::::::::::::
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
if (nocomp) go to 95
delx = xmax1 - xmin1
dely = ymax1 - ymin1
xlen = .8*amin1( 1.0, delx/dely )
ylen = .8*amin1( 1.0, dely/delx )
x1set = .1
x2set = .1+xlen
y1set = .1
y2set = .1+ylen
call set( .1,.1+xlen,.1,.1+ylen,xmin1,xmax1,ymin1,ymax1,1)
call drawbox(xmin1,ymin1,xmax1,ymax1,level)
level = lend
80 if (level .ge. lst) then
dx = hxposs(level) / 2.0
dy = hyposs(level) / 2.0
mptr = lstart(level)
90 if (mptr .ne. 0) then
loc = node(store2,mptr)
c nrow = node(maxnumrow,mptr) + 1
c ncol = node(maxnumcol,mptr) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
x1 = rnode(cornxlo,mptr)
y1 = rnode(cornylo,mptr)
x3 = rnode(cornxhi,mptr)
y3 = rnode(cornyhi,mptr)
if (plotbox) call drawbox(x1,y1,x3,y3,level)
xrat1 = (x1-xmin1)/(xmax1-xmin1)
xrat2 = (x3-xmin1)/(xmax1-xmin1)
yrat1 = (y1-ymin1)/(ymax1-ymin1)
yrat2 = (y3-ymin1)/(ymax1-ymin1)
call set(x1set+xrat1*(x2set-x1set),
1 x1set+xrat2*(x2set-x1set),
2 y1set+yrat1*(y2set-y1set),
3 y1set+yrat2*(y2set-y1set),
4 x1,x3,y1,y3,1)
c ::::::::: window into the rectangle defined by this grid :::::
x1 = x1 - dx
y1 = y1 - dy
x3 = x3 + dx
y3 = y3 + dy
x1ratio = (x1-xmin1)/(xmax1-xmin1)
x3ratio = (x3-xmin1)/(xmax1-xmin1)
y1ratio = (y1-ymin1)/(ymax1-ymin1)
y3ratio = (y3-ymin1)/(ymax1-ymin1)
call set(x1set+x1ratio*(x2set-x1set),
. x1set+x3ratio*(x2set-x1set),
. y1set+y1ratio*(y2set-y1set),
. y1set+y3ratio*(y2set-y1set),
. 1.5,float(nrow)-.5,1.5,float(ncol)-.5,1)
c ::::::::: plot the contour lines :::::
call setspv(alloc(loc),nrow,ncol,spval)
c call cpcnrc( alloc(loc), nrow, nrow, ncol,
c 1 cmin, cmax, cinc, 1, -1, 0 )
call conrec( alloc(loc), nrow, nrow, ncol,
1 cmin, cmax, cinc, -1, -1, 0 )
c ::::::::: reset the window to that of the coarsest grid :::::
call set( .1,.1+xlen,.1,.1+ylen,xmin1,xmax1,ymin1,ymax1,1)
mptr = node(levelptr,mptr)
goto 90
endif
level = level - 1
goto 80
endif
call set(0.,1.,0.,1.,xmin1,xmax1,ymin1,ymin1+1.2*(ymax1-ymin1),1)
xplace =xmin1 + (xmax1-xmin1)/4.
yplace = ymin1 + 1.1*(ymax1-ymin1)
call gstxci(1)
call pwrit(xplace,yplace,cline,clen,2,0,-1)
framecnt = framecnt + 1
write(7,199) framecnt, "GRIDPLOT", cline
call frame
c
95 continue
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c :::::::::::::: PRODUCE A PLOT OF ONLY THE FINEST LEVELS :::::::::::::
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
if (.not. fineplot .or. lfine .eq. 1) goto 999
xmin2 = 1.0e+10
ymin2 = 1.0e+10
xmax2 = -1.0e+10
ymax2 = -1.0e+10
level = max0(lfine-1,2)
205 if (level .le. lend) then
xmin2 = amin1( xmin2, xmin(level) )
ymin2 = amin1( ymin2, ymin(level) )
xmax2 = amax1( xmax2, xmax(level) )
ymax2 = amax1( ymax2, ymax(level) )
level = level + 1
goto 205
endif
delx = xmax2 - xmin2
dely = ymax2 - ymin2
xlen = .8*amin1( 1.0, delx/dely )
ylen = .8*amin1( 1.0, dely/delx )
x1set = .1
x2set = .1+xlen
y1set = .1
y2set = .1+ylen
call set(.1,.1+xlen,.1,.1+ylen,xmin2,xmax2,ymin2,ymax2,1)
c ::::::::: Now plot the finest grids with boxes :::::
level = max0(lfine-1,2)
225 if (level .le. lend) then
dx = hxposs(level) / 2.0
dy = hyposs(level) / 2.0
mptr = lstart(level)
220 if (mptr .ne. 0) then
loc = node(store2,mptr)
c nrow = node( maxnumrow, mptr ) + 1
c ncol = node( maxnumcol, mptr ) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
x1 = rnode(cornxlo,mptr)
y1 = rnode(cornylo,mptr)
x3 = rnode(cornxhi,mptr)
y3 = rnode(cornyhi,mptr)
call drawbox(x1,y1,x3,y3,level)
c ::::::::: window into the rectangle defined by this grid :::::
x1 = x1 - dx
y1 = y1 - dy
x3 = x3 + dx
y3 = y3 + dy
x1ratio = (x1-xmin2)/(xmax2-xmin2)
x3ratio = (x3-xmin2)/(xmax2-xmin2)
y1ratio = (y1-ymin2)/(ymax2-ymin2)
y3ratio = (y3-ymin2)/(ymax2-ymin2)
call set(x1set+x1ratio*(x2set-x1set),
. x1set+x3ratio*(x2set-x1set),
. y1set+y1ratio*(y2set-y1set),
. y1set+y3ratio*(y2set-y1set),
. 1.5,float(nrow)-.5,1.5,float(ncol)-.5,1)
c ::::::::: plot the contour lines :::::
call setspv(alloc(loc),nrow,ncol,spval)
call conrec( alloc(loc), nrow, nrow, ncol,
1 cmin, cmax, cinc, -1, -1, 0 )
c call cpcnrc( alloc(loc), nrow, nrow, ncol,
c 1 cmin, cmax, cinc, 1, -1, 0 )
c call set(minx,maxx,miny,maxy,xmin2,xmax2,ymin2,ymax2,1)
call set(.1,.1+xlen,.1,.1+ylen,xmin2,xmax2,ymin2,ymax2,1)
mptr = node(levelptr,mptr)
goto 220
endif
level = level + 1
goto 225
endif
call set(0.,1.,0.,1.,xmin2,xmax2,ymin2,ymin2+1.2*(ymax2-ymin2),1)
xplace = xmin2 + (xmax2-xmin2)/4.
yplace = ymin2 + 1.1*(ymax2-ymin2)
call gstxci(1)
call pwrit(xplace,yplace,fline,flen,2,0,-1)
framecnt = framecnt + 1
write(7,199) framecnt, "GRIDPLOT", fline
call frame
c ::::::::: reclaim the space used by temporary grids (store1&2) :::::
999 level = lst
230 if (level .le. lend) then
mptr = lstart(level)
240 if ( mptr .ne. 0 ) then
c nrow = node(maxnumrow,mptr) + 1
c ncol = node(maxnumcol,mptr) + 1
nrow = node(ndihi,mptr) - node(ndilo,mptr) + 3
ncol = node(ndjhi,mptr) - node(ndjlo,mptr) + 3
num = nrow*ncol
loc = node(store2,mptr)
call reclam(loc,num)
node(store2,mptr) = 0
loc = node(store1,mptr)
call reclam(loc,num)
node(store1,mptr) = 0
mptr = node(levelptr,mptr)
goto 240
endif
level = level + 1
goto 230
endif
c
199 format( "Frame# ",i3,2x,a8,2x,a50 )
return
end
c
c ---------------------------------------------------------c
c
subroutine setspv(a,m,n,spval)
dimension a(m,n)
data spval2/-1.e10/
c
c exterior points inside the rectangle have been marked by spval2.
c reset to spval.
c
do 25 i = 1,m
do 25 j = 1,n
if (a(i,j) .eq. spval2) a(i,j) = spval
20 continue
25 continue
return
end
c
c ------------------------------------------------------------
c
subroutine mnmax(a,nrow,ncol,vmin,vmax,spval2)
dimension a(nrow,ncol)
vmin = 1.0e+10
vmax = -1.0e+10
do 10 i = 2, nrow-1
do 20 j = 2, ncol-1
val = a(i,j)
c spval2 is special "blackout" value - don't count it.
if (val .eq. spval2) go to 20
vmin = amin1(val, vmin)
vmax = amax1(val, vmax)
20 continue
10 continue
return
end
c
c ------------------------------------------------------------
c
subroutine blackout(mptr,level,val,nrow,ncol,spval)
c
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c ::::::::: Array "val" contains the values of grid "mptr" at "level".
c ::::::::: For each cell of mptr, determine if this cell overlaps
c ::::::::: a grid at level+1. If so, black out the cell's value.
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
dimension val(nrow,ncol)
c
include "call.i"
c
dx = hxposs(level)
dy = hyposs(level)
x = rnode(cornxlo,mptr)
y = rnode(cornylo,mptr)
nptr = lstart(level+1)
10 if (nptr .ne. 0) then
xst = rnode(cornxlo,nptr)
yst = rnode(cornylo,nptr)
xend = rnode(cornxhi,nptr)
yend = rnode(cornyhi,nptr)
ist = max(ifix((xst-x)/dx+2.1), 1)
jst = max(ifix((yst-y)/dy+2.1), 1)
iend = min(ifix((xend-x)/dx+1.1), nrow)
jend = min(ifix((yend-y)/dy+1.1), ncol)
if ((ist .gt. iend) .or. (jst .gt. jend)) go to 40
do 20 i = ist, iend
do 20 j = jst, jend
val(i,j) = spval
20 continue
40 nptr = node(levelptr,nptr)
go to 10
endif
return
end
c
c ------------------------------------------------------------
c
subroutine addbndry( level, mptr, a, nrow, ncol, nvar,
1 xlo, ylo, xhi, yhi, spval, lst,
2 deltax, deltay)
c
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c ::::::::: Array "a" contains the values of the grid at "level"
c ::::::::: pointed to by "mptr".
c ::::::::: This routine fills in the boundry values array "a".
c ::::::::: The boundry value of a cell is determined by:
c ::::::::: (1) it is "spval" if the cell is outside the physical
c ::::::::: bounds of the problem defined by xlow, ylow, xhi, yhi.
c ::::::::: (2) if the cell overlaps another grid at "level" it is
c ::::::::: given the value of the overlapping cell.
c ::::::::: (3) if it does not overlap a grid at the same level
c ::::::::: it must overlap a cell at a coarser level. The
c ::::::::: four nearest neighbor cells of this coarser grid
c ::::::::: are interpolated to give a value to the cell.
c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
dimension a(nrow,ncol,nvar)
real xlo, ylo, xhi, yhi, interp
real q(4),v11(4),v21(4),v12(4),v22(4),regval(4)
integer level, mptr
logical outside, xinside, yinside
c
include "call.i"
c
interp(xx1,vv1,xx2,vv2,xx) = vv1 + (vv2-vv1)*(xx-xx1)/(xx2-xx1)
outside(x,y) = ( (x .lt. xlo) .or. (x .gt. xhi) .or.
1 (y .lt. ylo) .or. (y .gt. yhi) )
xbot = rnode(cornxlo,mptr)
ybot = rnode(cornylo,mptr)
do 70 i = 1, nrow
xinside = ( (i .gt. 1) .and. (i .lt. nrow) )
xcen = xbot + deltax*(float(i) - 1.5)
do 60 j = 1, ncol
yinside = ( (j .gt. 1) .and. (j .lt. ncol) )
if ( xinside .and. yinside ) goto 60
ycen = ybot + deltay*(float(j) - 1.5)
if (outside(xcen,ycen)) then
do 10 ivar = 1, nvar
c if (i .eq. 1) then
c if (j .eq. 1 .or. j .eq. ncol) then
c a(i,j,ivar) = spval
c else
c a(1,j,ivar) = a(2,j,ivar)
c endif
c else if (i .eq. nrow) then
c a(nrow,j,ivar) = a(nrow-1,j,ivar)
c else if (j .eq. ncol) then
c a(i,ncol,ivar) = a(i,ncol-1,ivar)
c else if (j .eq. 1) then
c a(i,1,ivar) = a(i,2,ivar)
c endif
10 a(i,j,ivar) = spval
c10 continue
goto 60
endif
call srchlev(xcen,ycen,level,ngrd,i1,j1,x1,y1)
if (ngrd .ne. 0) then
c :::::::: found at same level, just use that value :::::
c nr = node(maxnumrow,ngrd) + 1
c nc = node(maxnumcol,ngrd) + 1
nr = node(ndihi,ngrd) - node(ndilo,ngrd) + 3
nc = node(ndjhi,ngrd) - node(ndjlo,ngrd) + 3
loc = node(tempptr,ngrd)
call getval(alloc(loc),nr,nc,i1,j1,q,nvar)
do 15 ivar = 1, nvar
15 a(i,j,ivar) = q(ivar)
go to 60
endif
c :::::::: not found at same level, search lower level
c :::::::: and interplolate from four nearest points
lev = level - 1
call srchlev(xcen,ycen,lev,ngrd,i1,j1,x1,y1)
if (ngrd .ne. 0) then
dx = hxposs(lev)
dy = hyposs(lev)
xsgn = sign( 1.0, xcen - x1 )
ysgn = sign( 1.0, ycen - y1 )
x2 = x1 + xsgn*dx
y2 = y1 + ysgn*dy
i2 = i1 + nint(xsgn)
j2 = j1 + nint(ysgn)
c nr = node(maxnumrow,ngrd) + 1
c nc = node(maxnumcol,ngrd) + 1
nr = node(ndihi,ngrd) - node(ndilo,ngrd) + 3
nc = node(ndjhi,ngrd) - node(ndjlo,ngrd) + 3
loc = node(tempptr,ngrd)
call getval(alloc(loc),nr,nc,i1,j1,v11,nvar)
if (outside(x1,y2)) then
c :::: (x2,y1) must be inside. Use 2 point interp.
call srchlev(x2,y1,lev,ng,ii,jj,xx,yy)
c nr = node(maxnumrow,ng) + 1
c nc = node(maxnumcol,ng) + 1
nr = node(ndihi,ng) - node(ndilo,ng) + 3
nc = node(ndjhi,ng) - node(ndjlo,ng) + 3
loc = node(tempptr,ng)
call getval(alloc(loc),nr,nc,ii,jj,v21,nvar)
if ((abs(v11(1)).eq.spval).and.
. (abs(v21(1)).ne.spval)) then
call fix(v11,v21)
elseif ((abs(v11(1)).ne.spval).and.
. (abs(v21(1)).eq.spval)) then
call fix(v21,v11)
endif
do 35 ivar = 1, nvar
a(i,j,ivar)=interp(x1,v11(ivar),x2,v21(ivar),xcen)
35 continue
elseif (outside(x2,y1)) then
c :::: (x1,y2) must be inside. Use 2 point interp.
call srchlev(x1,y2,lev,ng,ii,jj,xx,yy)
c nr = node(maxnumrow,ng) + 1
c nc = node(maxnumcol,ng) + 1
nr = node(ndihi,ng) - node(ndilo,ng) + 3
nc = node(ndjhi,ng) - node(ndjlo,ng) + 3
loc = node(tempptr,ng)
call getval(alloc(loc),nr,nc,ii,jj,v12,nvar)
if ((abs(v11(1)).eq.spval).and.
. (abs(v12(1)).ne.spval)) then
call fix(v11,v12)
elseif ((abs(v11(1)).ne.spval).and.
. (abs(v12(1)).eq.spval)) then
call fix(v12,v11)
endif
do 40 ivar = 1, nvar
a(i,j,ivar)=interp(y1,v11(ivar),y2,v12(ivar),ycen)
40 continue
else
c :::: All points inside. Use 4 point interp.
call srchlev(x1,y2,lev,ng,ii,jj,xx,yy)
c nr = node(maxnumrow,ng) + 1
c nc = node(maxnumcol,ng) + 1
nr = node(ndihi,ng) - node(ndilo,ng) + 3
nc = node(ndjhi,ng) - node(ndjlo,ng) + 3
loc = node(tempptr,ng)
call getval(alloc(loc),nr,nc,ii,jj,v12,nvar)
call srchlev(x2,y1,lev,ng,ii,jj,xx,yy)
c nr = node(maxnumrow,ng) + 1
c nc = node(maxnumcol,ng) + 1
nr = node(ndihi,ng) - node(ndilo,ng) + 3
nc = node(ndjhi,ng) - node(ndjlo,ng) + 3
loc = node(tempptr,ng)
call getval(alloc(loc),nr,nc,ii,jj,v21,nvar)
call srchlev(x2,y2,lev,ng,ii,jj,xx,yy)
c nr = node(maxnumrow,ng) + 1
c nc = node(maxnumcol,ng) + 1
nr = node(ndihi,ng) - node(ndilo,ng) + 3
nc = node(ndjhi,ng) - node(ndjlo,ng) + 3
loc = node(tempptr,ng)
call getval(alloc(loc),nr,nc,ii,jj,v22,nvar)
c ::: are all coarse values in the regular domain? :::
if ( ((abs(v11(1)).ne.spval)
. .and.(abs(v12(1)).ne.spval)
. .and.(abs(v21(1)).ne.spval)
. .and.(abs(v22(1)).ne.spval)) .or.
. ((abs(v11(1)).eq.spval)
. .and.(abs(v12(1)).eq.spval)
. .and.(abs(v21(1)).eq.spval)
. .and.(abs(v22(1)).eq.spval)) ) then
do 45 ivar = 1, nvar
v1 = interp(y1,v11(ivar),y2,v12(ivar),ycen)
v2 = interp(y1,v21(ivar),y2,v22(ivar),ycen)
c :::: interpolate in x direction :::::
a(i,j,ivar) = interp(x1,v1,x2,v2,xcen)
45 continue
else
c ::: just find one non-special value to use
c ::: piecewise constant interp.
if (abs(v11(1)).ne.spval) then
call fix(regval,v11)
elseif (abs(v12(1)).ne.spval) then
call fix(regval,v12)
elseif (abs(v21(1)).ne.spval) then
call fix(regval,v21)
elseif (abs(v22(1)).ne.spval) then
call fix(regval,v22)
endif
do 234 ivar = 1, nvar
a(i,j,ivar) = regval(ivar)
234 continue
endif
go to 60
endif
c
c :::::::: should never fall through this loop :::::
else
write(16,100) xcen, ycen
100 format( "no coarse grid contains ",2(e13.6,2x))
do 55 ivar = 1, nvar
if (i .eq. 1) then
if (j .eq. 1 .or. j .eq. ncol) then
a(i,j,ivar) = spval
else
a(1,j,ivar) = a(2,j,ivar)
endif
else if (i .eq. nrow) then
a(nrow,j,ivar) = a(nrow-1,j,ivar)
else if (j .eq. ncol) then
a(i,ncol,ivar) = a(i,ncol-1,ivar)
else if (j .eq. 1) then
a(i,1,ivar) = a(i,2,ivar)
endif
55 continue
call clsgks
c stop
endif
c
60 continue
70 continue
return
end
c
c ---------------------------------------------------
c
subroutine fix(v1,v2)
dimension v1(4), v2(4)
c
c set v1 equal to v2 so that interpolation in addbndry
c routine either sets boundary value to spval or to a "real"
c (but obviously not very accurate) value
c
do 10 i = 1, 4
v1(i) = v2(i)
10 continue
return
end
c
c ------------------------------------------------------------
c
subroutine getval(a,nrows,ncols,i,j,q,nvar)
dimension a(nrows,ncols,nvar), q(nvar)
c
do 10 ivar = 1, nvar
q(ivar) = a(i,j,ivar)
10 continue
c
return
end
c
c ------------------------------------------------------------
c
subroutine srchlev(xcen,ycen,level,ngrd,i,j,x,y)
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c ::::::::: Search all the grids at "level" until one is found
c ::::::::: that contains the point (xcen,ycen).
c ::::::::: If no such grid is found, ngrd will be returned as zero.
c ::::::::: If one is found, ngrd will point to that grid,
c ::::::::: (i,j) will be the indicies of the cell it is found in,
c ::::::::: and (x,y) will be the coordinates of the center of the cell.
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
include "call.i"
logical inside
c
ngrd = lstart(level)
10 if (ngrd .ne. 0) then
xlow = rnode(cornxlo,ngrd)
ylow = rnode(cornylo,ngrd)
xhigh = rnode(cornxhi,ngrd)
yhigh = rnode(cornyhi,ngrd)
inside = ( (xcen .ge. xlow) .and. (xcen .le. xhigh) .and.
1 (ycen .ge. ylow) .and. (ycen .le. yhigh) )
if (inside) then
dx = hxposs(level)
dy = hyposs(level)
i = int( (xcen - xlow) / dx ) + 2
j = int( (ycen - ylow) / dy ) + 2
x = xlow + dx*(float(i) - 1.5)
y = ylow + dy*(float(j) - 1.5)
go to 99
endif
ngrd = node(levelptr,ngrd)
goto 10
endif
99 return
end
c
c ------------------------------------------------------------
c
subroutine restore( orig, black, fixed, nrow, ncol,
1 mptr, spval, xlo, ylo, xhi, yhi )
c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c ::::::::: Array "black" is identical to array "orig" except
c ::::::::: that some cells may be blackened out because they
c ::::::::: overlap cells in finer grids.
c ::::::::: This routine creates the array "fixed" to be identical
c ::::::::: to "black" except that if restores the origional
c ::::::::: value to a blackened cell if that cell satisfies
c ::::::::: the following conditions:
c ::::::::: (1) the cell is within the physical problem domain
c ::::::::: (2) the cell is blackened out.
c ::::::::: (3) the cell is adjacent to a cell that is not
c ::::::::: blackened out.
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
parameter (eps=1.0)
real orig(nrow, ncol), fixed(nrow, ncol), black(nrow, ncol)
logical notblack, outside
c
include "call.i"
c
notblack(xx) = ((xx .lt. spval-eps) .or. (xx .gt. spval+eps))
c ::::::::: get the physical coordinates of the grid we are restoring
level = node(nestlevel,mptr)
dx = hxposs(level)
dy = hyposs(level)
x = rnode(cornxlo,mptr)
y = rnode(cornylo,mptr)
do 10 i = 1, nrow
klow = -1
if (i .eq. 1) klow = 0
khigh = 1
if (i .eq. nrow) khigh = 0
do 20 j = 1, ncol
fixed(i,j) = black(i,j)
c ::::::: check only blackened cells :::::
if (notblack( black(i,j) )) goto 20
c ::::::: check if cell is within bounds of physical domain :::::
xcen = x + dx*(float(i) - 1.5)
ycen = y + dy*(float(j) - 1.5)
outside = ( (xcen .le. xlo) .or. (xcen .ge. xhi) .or.
1 (ycen .le. ylo) .or. (ycen .ge. yhi) )
if (outside) goto 20
llow = -1
if (j .eq. 1) llow = 0
lhigh = 1
if (j .eq. ncol) lhigh = 0
c ::::::: check nearest neighbors for non-black cells :::::
do 30 k = klow,khigh
do 40 l = llow, lhigh
if (notblack(black(i+k,j+l))) then
c ::::::: restore to origional value :::::
c write(16,200) mptr, i, j
200 format("Restoring value: mptr, I, J ",3(i4,3x) )
fixed(i,j) = orig(i,j)
goto 20
endif
40 continue
30 continue
20 continue
10 continue
c
return
end
c
c ------------------------------------------------------------------
c
subroutine drawbox(xlo,ylo,xhi,yhi,level)
dimension icl(5)
data icl/10,40,100,110,120/
call gsplci(icl(level))
call gstxci(icl(level))
call line(xlo,ylo,xlo,yhi)
call line(xlo,yhi,xhi,yhi)
call line(xhi,yhi,xhi,ylo)
call line(xhi,ylo,xlo,ylo)
return
end
c
c ------------------------------------------------------------------
c
subroutine dtree (lst,lend,line,linelen,gprint2)
c
include "call.i"
character*50 line
integer linelen
logical gprint2
c
c *******************************************************************
c dtree - draw the grids in the subtree starting at level lst to lend.
c (just outline the grids, not their values)
c ******************************************************************
c
c
c set scaling using (x,y) min/max of coarsest level rectangles.
xmin = 10.e32
xmax = - 10.e32
ymin = 10.e32
ymax = - 10.e32
c
mptr = lstart(lst)
10 if (mptr .ne. 0) then
x1 = rnode(cornxlo,mptr)
x2 = rnode(cornxhi,mptr)
xmin = amin1(x1,x2,xmin)
xmax = amax1(x1,x2,xmax)
y1 = rnode(cornylo,mptr)
y2 = rnode(cornyhi,mptr)
ymin = amin1(y1,y2,ymin)
ymax = amax1(y1,y2,ymax)
mptr = node(levelptr,mptr)
go to 10
endif
c
c add border
ratxy =.8*amin1(1.,(xmax-xmin)/(ymax-ymin))
ratyx =.8*amin1(1.,(ymax-ymin)/(xmax-xmin))
call set(.1,.1+ratxy,.1,.1+ratyx,xmin,xmax,ymin,ymax,1)
c
c actually generate call to draw each rectangle now.
level = lst
30 if (level .le. lend) then
mptr = lstart(level)
40 if (mptr .ne. 0) then
call drect(mptr,gprint2)
mptr = node(levelptr,mptr)
go to 40
endif
level = level + 1
go to 30
endif
c
call set(0.,1.,0.,1.,xmin,xmax,ymin,ymin+1.2*(ymax-ymin),1)
xthird = xmin + (xmax-xmin)/3.
call pwrit(xthird,ymin+1.1*(ymax-ymin),line,linelen,2,0,-1)
c
c reset plotting coords. in case call came from dclust - more to draw.
c
call set(.1,.1+ratxy,.1,.1+ratyx,xmin,xmax,ymin,ymax,1)
c
return
end
subroutine drect(mptr,gprint2)
c
include "call.i"
dimension x(4), y(4)
logical gprint2
character*3 numg
c
c ******************************************************************
c drect - draw rectangle of grid mptr.
c if gprint2 true, print grid number in grid.
c
c this routine assumes dtree already called, appropriate graphics
c calls already made.
c ******************************************************************
c
x(1) = rnode(cornxlo,mptr)
y(1) = rnode(cornylo,mptr)
x(2) = rnode(cornxlo,mptr)
y(2) = rnode(cornyhi,mptr)
x(3) = rnode(cornxhi,mptr)
y(3) = rnode(cornyhi,mptr)
x(4) = rnode(cornxhi,mptr)
y(4) = rnode(cornylo,mptr)
level = node(nestlevel, mptr)
hx = hxposs(level)
hy = hyposs(level)
do 10 i = 1,4
call line(x(i),y(i),x((mod(i,4))+1),y((mod(i,4))+1))
10 continue
c
c print grid number?
c
if (gprint2) then
xplace = x(1) + 2.*hx
yplace = y(2) - 3.*hy
write(numg(1:3),'(i3)') mptr
call pwrit(xplace,yplace,numg,3,0,0,0)
endif
return
end
c
c -------------------------------------------------------
c
subroutine contor(a,m,n,header,hlen,mptr)
dimension a(m,n)
character header*50
include "call.i"
common /conre1/ ioffp, spval
common /conre4/ sizel, sizem, sizep, nrep, ncrt, ilab,
1 nulbll, ioffd, ext, ioffm, isolid, nla,
2 nlm, xlt, ybt, side
data ncons/40/
c
ioffp = 1
nlm = 40
spval = 1.0e10
spval2 = -1.e10
call cpseti('LLP - line label position',0)
c call cpsetc('HLT - high/low label text',' ')
c call cpsetc('ILT - informational label text',' ')
call cpseti('NCL - number of contour levels',ncons)
call cpseti('SET - do set call flag',0)
call cpsetr('SPV - set special value ',spval)
c
c ::::::::: only work with data in the subarray (2..m-1,2..n-1) :::::
flm = float(m-2)
fln = float(n-2)
xlen = .8*amin1(1.,flm/fln)
ylen = .8*amin1(1.,fln/flm)
xmax = .1+xlen
ymax = .1+ylen
fncons = float(ncons)
totmax = -1.e+10
totmin = 1.e+10
do 100 i=2,m-1
do 100 j=2,n-1
if (a(i,j) .eq. spval2) go to 100
totmax = amax1(a(i,j),totmax)
totmin = amin1(a(i,j),totmin)
100 continue
write(7,*)" grid ",mptr," ranges from ",totmin," to ", totmax
c
if (ioffp .eq. 0) go to 30
do 25 i = 2,m-1
do 25 j = 2,n-1
if (a(i,j) .eq. spval2) a(i,j) = spval
25 continue
30 continue
c
33 dela = totmax - totmin
finc = dela/fncons
fhi = totmax -.5*finc
flo = totmin +.5*finc
call set(.1,xmax,.1,ymax,.5,flm+.5,.5,fln+.5,1)
call perim(10,1,10,1)
call conrec(a(2,2),m,m-2,n-2,flo,fhi,finc,-1,-1,0)
c call cpcnrc(a(2,2),m,m-2,n-2,flo,fhi,finc,1,-1,0)
c
c prepare titles
c
xthird = flm/3.
xplace = flm/10.
yplace = 1.1*fln
call set(0.,1.,0.,1.,0.,flm,0.,1.2*fln,1)
call pwrit(xplace,yplace,header,hlen,2,0,-1)
c
call frame
return
end
SUBROUTINE CONREC (Z,L,M,N,FLO,HI,FINC,NSET,NHI,NDOT)
C
C +-----------------------------------------------------------------+
C | |
C | Copyright (C) 1987 by UCAR |
C | University Corporation for Atmospheric Research |
C | All Rights Reserved |
C | |
C | NCARGRAPHICS Version 2.00 |
C | |
C +-----------------------------------------------------------------+
C
C
C
C DIMENSION OF Z(L,N)
C ARGUMENTS
C
C LATEST REVISION June, 1987
C
C PURPOSE CONREC draws a contour map from data stored
C in a rectangular array, labeling the lines.
C
C USAGE If the following assumptions are met, use
C
C CALL EZCNTR (Z,M,N)
C
C ASSUMPTIONS:
C --All of the array is to be contoured.
C --Contour levels are picked
C internally.
C --Contouring routine picks scale
C factors.
C --Highs and lows are marked.
C --Negative lines are drawn with a
C dashed line pattern.
C --EZCNTR calls FRAME after drawing the
C contour map.
C
C If these assumptions are not met, use
C
C CALL CONREC (Z,L,M,N,FLO,HI,FINC,NSET,
C NHI,NDOT)
C
C ARGUMENTS
C
C ON INPUT Z
C FOR EZCNTR M by N array to be contoured.
C
C M
C First dimension of Z.
C
C N
C Second dimension of Z.
C
C ON OUTPUT All arguments are unchanged.
C FOR EZCNTR
C
C ON INPUT Z
C FOR CONREC The (origin of the) array to be
C contoured. Z is dimensioned L by N.
C
C L
C The first dimension of Z in the calling
C program.
C
C M
C The number of data values to be contoured
C in the X-direction (the first subscript
C direction). When plotting an entire
C array, L = M.
C
C N
C The number of data values to be contoured
C in the Y-direction (the second subscript
C direction).
C
C FLO
C The value of the lowest contour level.
C If FLO = HI = 0., a value rounded up from
C the minimum Z is generated by CONREC.
C
C HI
C The value of the highest contour level.
C If HI = FLO = 0., a value rounded down
C from the maximum Z is generated by
C CONREC.
C
C FINC
C > 0 Increment between contour levels.
C = 0 A value, which produces between 10
C and 30 contour levels at nice values,
C is generated by CONREC.
C < 0 The number of levels generated by
C CONREC is ABS(FINC).
C
C NSET
C Flag to control scaling.
C = 0 CONREC automatically sets the
C window and viewport to properly
C scale the frame to the standard
C configuration.
C The GRIDAL entry PERIM is
C called and tick marks are placed
C corresponding to the data points.
C > 0 CONREC assumes that the user
C has set the window and viewport
C in such a way as to properly
C scale the plotting
C instructions generated by CONREC.
C PERIM is not called.
C < 0 CONREC generates coordinates so as
C to place the (untransformed) contour
C plot within the limits of the
C user's current window and
C viewport. PERIM is not called.
C
C NHI
C Flag to control extra information on the
C contour plot.
C = 0 Highs and lows are marked with an H
C or L as appropriate, and the value
C of the high or low is plotted under
C the symbol.
C > 0 The data values are plotted at
C each Z point, with the center of
C the string indicating the data
C point location.
C < 0 Neither of the above.
C
C NDOT
C A 10-bit constant designating the desired
C dashed line pattern.
C If ABS(NDOT) = 0, 1, or 1023, solid lines
C are drawn.
C > 0 NDOT pattern is used for all lines.
C < 0 ABS(NDOT) pattern is used for nega-
C tive-valued contour lines, and solid is
C used for positive-valued contours.
C CONREC converts NDOT
C to a 16-bit pattern and DASHDB is used.
C See DASHDB comments in the DASHLINE
C documentation for details.
C
C
C
C ON OUTPUT All arguments are unchanged.
C FOR CONREC
C
C
C ENTRY POINTS CONREC, CLGEN, REORD, STLINE, DRLINE,
C MINMAX, PNTVAL, CALCNT, EZCNTR, CONBD
C
C COMMON BLOCKS INTPR, RECINT, CONRE1, CONRE2, CONRE3,
C CONRE4,CONRE5
C
C REQUIRED LIBRARY Standard version: DASHCHAR, which at
C ROUTINES NCAR is loaded by default.
C Smooth version: DASHSMTH which must be
C requested at NCAR.
C Both versions require GRIDAL, the
C ERPRT77 package, and the SPPS.
C
C REQUIRED GKS LEVEL 0A
C
C I/O Plots contour map.
C
C PRECISION Single
C
C LANGUAGE FORTRAN 77
C
C HISTORY Replaces old contouring package called
C CALCNT at NCAR.
C
C ALGORITHM Each line is followed to completion. Points
C along a line are found on boundaries of the
C (rectangular) cells. These points are
C connected by line segments using the
C software dashed line package, DASHCHAR.
C DASHCHAR is also used to label the
C lines.
C
C NOTE To draw non-uniform contour levels, see
C the comments in CLGEN. To make special
C modifications for specific needs see the
C explanation of the internal parameters
C below.
C
C TIMING Varies widely with size and smoothness of
C Z.
C
C INTERNAL PARAMETERS NAME DEFAULT FUNCTION
C ---- ------- --------
C
C ISIZEL 1 Size of line labels,
C as per the size definitions
C given in the SPPS
C documentation for WTSTR.
C
C ISIZEM 2 size of labels for minimums
C and maximums,
C as per the size definitions
C given in the SPPS
C documentation for WTSTR.
C
C ISIZEP 0 Size of labels for data
C point values as per the size
C definitions given in the SPPS
C documentation for WTSTR.
C
C NLA 16 Approximate number of
C contour levels when
C internally generated.
C
C NLM 40 Maximum number of contour
C levels. If this is to be
C increased, the dimensions
C of CL and RWORK in CONREC
C must be increased by the
C same amount.
C
C XLT .05 Left hand edge of the plot
C (0.0 is the left edge of
C the frame and 1.0 is the
C right edge of the frame.)
C
C YBT .05 Bottom edge of the plot
C (0.0 is the bottom of the
C frame and 1.0 is the top
C of the frame.)
C
C SIDE 0.9 Length of longer edge of
C plot (see also EXT).
C
C NREP 6 Number of repetitions of
C the dash pattern between
C line labels.
C
C NCRT 4 Number of CRT units per
C element (bit) in the dash
C pattern.
C
C ILAB 1 Flag to control the drawing
C of line labels.
C . ILAB non-zero means label
C the lines.
C . ILAB = 0 means do not
C label the lines.
C
C NULBLL 3 Number of unlabeled lines
C between labeled lines. For
C example, when NULBLL = 3,
C every fourth level is
C labeled.
C
C IOFFD 0 Flag to control
C normalization of label
C numbers.
C . IOFFD = 0 means include
C decimal point when
C possible (do not
C normalize unless
C required).
C . IOFFD non-zero means
C normalize all label
C numbers and output a
C scale factor in the
C message below the graph.
C
C EXT .25 Lengths of the sides of the
C plot are proportional to M
C and N (when CONREC sets
C the window and viewport).
C In extreme cases, when
C MIN(M,N)/MAX(M,N) is less
C than EXT, CONREC
C produces a square plot.
C
C IOFFP 0 Flag to control special
C value feature.
C . IOFFP = 0 means special
C value feature not in use.
C . IOFFP non-zero means
C special value feature in
C use. (SPVAL is set to the
C special value.) Contour
C lines will then be
C omitted from any cell
C with any corner equal to
C the special value.
C
C SPVAL 0. Contains the special value
C when IOFFP is non-zero.
C
C IOFFM 0 Flag to control the message
C below the plot.
C . IOFFM = 0 if the message
C is to be plotted.
C . IOFFM non-zero if the
C message is to be omitted.
C
C ISOLID 1023 Dash pattern for
C non-negative contour lines.
C
C
EXTERNAL CONBD
C
SAVE
CHARACTER*1 IGAP ,ISOL ,RCHAR
CHARACTER ENCSCR*22 ,IWORK*126
DIMENSION LNGTHS(5) ,HOLD(5) ,WNDW(4) ,VWPRT(4)
DIMENSION Z(L,N) ,CL(50) ,RWORK(50) ,LASF(13)
COMMON /INTPR/ PAD1, FPART, PAD(8)
COMMON /CONRE1/ IOFFP ,SPVAL
COMMON /CONRE3/ IXBITS ,IYBITS
COMMON /CONRE4/ ISIZEL ,ISIZEM ,ISIZEP ,NREP ,
1 NCRT ,ILAB ,NULBLL ,IOFFD ,
2 EXT ,IOFFM ,ISOLID ,NLA ,
3 NLM ,XLT ,YBT ,SIDE
COMMON /CONRE5/ SCLY
COMMON /RECINT/ IRECMJ ,IRECMN ,IRECTX
DATA LNGTHS(1),LNGTHS(2),LNGTHS(3),LNGTHS(4),LNGTHS(5)
1 / 12, 3, 20, 9, 17 /
DATA ISOL, IGAP /'$', ''''/
C
C ISOL AND IGAP (DOLLAR-SIGN AND APOSTROPHE) ARE USED TO CONSTRUCT PAT-
C TERNS PASSED TO ROUTINE DASHDC IN THE SOFTWARE DASHED-LINE PACKAGE.
C
C
C
C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
C
CALL Q8QST4 ('GRAPHX','CONREC','CONREC','VERSION 01')
C
C NONSMOOTHING VERSION
C
C
C
C CALL RESET FOR COMPATIBILITY WITH ALL DASH ROUTINES(EXCEPT DASHLINE)
C
CALL RESET
C
C GET NUMBER OF BITS IN INTEGER ARITHMETIC
C
IARTH = I1MACH(8)
IXBITS = 0
DO 101 I=1,IARTH
IF (M .LE. (2**I-1)) GO TO 102
IXBITS = I+1
101 CONTINUE
102 IYBITS = 0
DO 103 I=1,IARTH
IF (N .LE. (2**I-1)) GO TO 104
IYBITS = I+1
103 CONTINUE
104 IF ((IXBITS*IYBITS).GT.0 .AND. (IXBITS+IYBITS).LE.24) GO TO 105
C
C REPORT ERROR NUMBER ONE
C
IWORK = 'CONREC - DIMENSION ERROR - M*N .GT. (2**IARTH) M =
+ N = '
WRITE (IWORK(56:62),'(I6)') M
WRITE (IWORK(73:79),'(I6)') N
CALL SETER( IWORK, 1, 1 )
RETURN
105 CONTINUE
C
C INQUIRE CURRENT TEXT AND LINE COLOR INDEX
C
CALL GQTXCI ( IERR, ITXCI )
CALL GQPLCI ( IERR, IPLCI )
C
C SET LINE AND TEXT ASF TO INDIVIDUAL
C
CALL GQASF ( IERR, LASF )
LSV3 = LASF(3)
LSV10 = LASF(10)
LASF(3) = 1
LASF(10) = 1
CALL GSASF ( LASF )
C
GL = FLO
HA = HI
GP = FINC
MX = L
NX = M
NY = N
IDASH = NDOT
NEGPOS = ISIGN(1,IDASH)
IDASH = IABS(IDASH)
IF (IDASH.EQ.0 .OR. IDASH.EQ.1) IDASH = ISOLID
C
C SET CONTOUR LEVELS.
C
CALL CLGEN (Z,MX,NX,NY,GL,HA,GP,NLA,NLM,CL,NCL,ICNST)
C
C FIND MAJOR AND MINOR LINES
C
IF (ILAB .NE. 0) CALL REORD (CL,NCL,RWORK,NML,NULBLL+1)
IF (ILAB .EQ. 0) NML = 0
C
C SAVE CURRENT NORMALIZATION TRANS NUMBER NTORIG AND LOG SCALING FLAG
C
CALL GQCNTN ( IERR, NTORIG )
CALL GETUSV ('LS',IOLLS)
C
C SET UP SCALING
C
CALL GETUSV ( 'YF' , IYVAL )
SCLY = 1.0 / ISHIFT ( 1, 15 - IYVAL )
C
IF (NSET) 106,107,111
106 CALL GQNT ( NTORIG,IERR,WNDW,VWPRT )
X1 = VWPRT(1)
X2 = VWPRT(2)
Y1 = VWPRT(3)
Y2 = VWPRT(4)
C
C SAVE NORMALIZATION TRANS 1
C
CALL GQNT (1,IERR,WNDW,VWPRT)
C
C DEFINE NORMALIZATION TRANS AND LOG SCALING
C
CALL SET(X1, X2, Y1, Y2, 1.0, FLOAT(NX), 1.0, FLOAT(NY), 1)
GO TO 111
107 CONTINUE
X1 = XLT
X2 = XLT+SIDE
Y1 = YBT
Y2 = YBT+SIDE
X3 = NX
Y3 = NY
IF (AMIN1(X3,Y3)/AMAX1(X3,Y3) .LT. EXT) GO TO 110
IF (NX-NY) 108,110,109
108 X2 = SIDE*X3/Y3+XLT
GO TO 110
109 Y2 = SIDE*Y3/X3+YBT
C
C SAVE NORMALIZATION TRANS 1
C
110 CALL GQNT ( 1, IERR, WNDW, VWPRT )
C
C DEFINE NORMALIZATION TRANS 1 AND LOG SCALING
C
CALL SET(X1,X2,Y1,Y2,1.0,X3,1.0,Y3,1)
C
C DRAW PERIMETER
C
CALL PERIM (NX-1,1,NY-1,1)
111 IF (ICNST .NE. 0) GO TO 124
C
C SET UP LABEL SCALING
C
IOFFDT = IOFFD
IF (GL.NE.0.0 .AND. (ABS(GL).LT.0.1 .OR. ABS(GL).GE.1.E5))
1 IOFFDT = 1
IF (HA.NE.0.0 .AND. (ABS(HA).LT.0.1 .OR. ABS(HA).GE.1.E5))
1 IOFFDT = 1
ASH = 10.**(3-IFIX(ALOG10(AMAX1(ABS(GL),ABS(HA),ABS(GP)))-5000.)-
1 5000)
IF (IOFFDT .EQ. 0) ASH = 1.
IF (IOFFM .NE. 0) GO TO 115
IWORK ='CONTOUR FROM TO CONTOUR INTERVAL
1 OF PT(3,3)= LABELS SCALED BY'
HOLD(1) = GL
HOLD(2) = HA
HOLD(3) = GP
HOLD(4) = Z(3,3)
HOLD(5) = ASH
NCHAR = 0
DO 114 I=1,5
WRITE ( ENCSCR, '(G13.5)' ) HOLD(I)
NCHAR = NCHAR+LNGTHS(I)
DO 113 J=1,13
NCHAR = NCHAR+1
IWORK(NCHAR:NCHAR) = ENCSCR(J:J)
113 CONTINUE
114 CONTINUE
IF (ASH .EQ. 1.) NCHAR = NCHAR-13-LNGTHS(5)
C
C SET TEXT INTENSITY TO LOW, AND WRITE TITLE USING NORMALIZATION
C TRANS NUMBER 0
C
CALL GSTXCI (IRECTX)
CALL GETUSV('LS',LSO)
CALL SETUSV('LS',1)
CALL GSELNT (0)
CALL WTSTR ( 0.5, 0.015625, IWORK(1:NCHAR), 0, 0, 0 )
CALL SETUSV('LS',LSO)
CALL GSELNT (1)
C
C
C
C * * * * * * * * * *
C * * * * * * * * * *
C
C
C PROCESS EACH LEVEL
C
115 FPART = .5
C
DO 123 I=1,NCL
CONTR = CL(I)
NDASH = IDASH
IF (NEGPOS.LT.0 .AND. CONTR.GE.0.) NDASH = ISOLID
C
C CHANGE 10 BIT PATTERN TO 10 CHARACTER PATTERN.
C
DO 116 J=1,10
IBIT = IAND(ISHIFT(NDASH,(J-10)),1)
RCHAR = IGAP
IF (IBIT .NE. 0) RCHAR = ISOL
IWORK(J:J) = RCHAR
116 CONTINUE
IF (I .GT. NML) GO TO 121
C
C SET UP MAJOR LINE (LABELED)
C
C
C NREP REPITITIONS OF PATTERN PER LABEL.
C
NCHAR = 10
IF (NREP .LT. 2) GO TO 119
DO 118 J=1,10
NCHAR = J
RCHAR = IWORK(J:J)
DO 117 K=2,NREP
NCHAR = NCHAR+10
IWORK(NCHAR:NCHAR) = RCHAR
117 CONTINUE
118 CONTINUE
119 CONTINUE
C
C PUT IN LABEL.
C
CALL ENCD (CONTR,ASH,ENCSCR,NCUSED,IOFFDT)
DO 120 J=1,NCUSED
NCHAR = NCHAR+1
IWORK(NCHAR:NCHAR) = ENCSCR(J:J)
120 CONTINUE
GO TO 122
C
C SET UP MINOR LINE (UNLABELED).
C
121 CONTINUE
C
C SET LINE INTENSITY TO LOW
C
CALL GSPLCI ( IRECMN )
NCHAR = 10
122 CALL DASHDC ( IWORK(1:NCHAR),NCRT, ISIZEL )
C
C
C DRAW ALL LINES AT THIS LEVEL.
C
CALL STLINE (Z,MX,NX,NY,CONTR)
C
C SET LINE INTENSITY TO HIGH
C
CALL GSPLCI ( IRECMJ )
C
123 CONTINUE
C
C FIND RELATIVE MINIMUMS AND MAXIMUMS IF WANTED, AND MARK VALUES IF
C WANTED.
C
IF (NHI .EQ. 0) CALL MINMAX (Z,MX,NX,NY,ISIZEM,ASH,IOFFDT)
IF (NHI .GT. 0) CALL MINMAX (Z,MX,NX,NY,ISIZEP,-ASH,IOFFDT)
FPART = 1.
GO TO 127
124 CONTINUE
IWORK = 'CONSTANT FIELD'
WRITE( ENCSCR, '(G22.14)' ) GL
DO 126 I=1,22
IWORK(I+14:I+14) = ENCSCR(I:I)
126 CONTINUE
C
C WRITE TITLE USING NORMALIZATION TRNS 0
C
CALL GETUSV('LS',LSO)
CALL SETUSV('LS',1)
CALL GSELNT (0)
CALL WTSTR ( 0.09765, 0.48825, IWORK(1:36), 3, 0, -1 )
C
C RESTORE NORMALIZATION TRANS 1, LINE AND TEXT INTENSITY TO ORIGINAL
C
127 IF (NSET.LE.0) THEN
CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
- WNDW(1),WNDW(2),WNDW(3),WNDW(4),IOLLS)
END IF
CALL GSPLCI ( IPLCI )
CALL GSTXCI ( ITXCI )
C
C SELECT ORIGINAL NORMALIZATION TRANS NUMBER NTORIG, AND RESTORE ASF
C
CALL GSELNT ( NTORIG )
LASF(3) = LSV3
LASF(10) = LSV10
CALL GSASF ( LASF )
C
RETURN
C
C
END
c
c -------------------------------------------------------
c
subroutine ctable
c
c set up the color table
c
ncolors = 125
neach = ncolors/3 - 1
c
c
c # interpolate between blue and blue-green
ilower = 2
iupper = ilower + neach
irange = iupper - ilower
range = float(irange)
blue = 1.0
red = 0.0
do 20 i = ilower, iupper
green = float(i-ilower)/range
call gscr(1,i,red,green,blue)
20 continue
c
c # interpolate between blue-green and pure green
ilower = iupper + 1
iupper = ilower + neach/2
irange = iupper - ilower
range = float(irange)
green = 1.0
red = 0.0
do 30 i = ilower, iupper
blue = float(iupper-i)/range
call gscr(1,i,red,green,blue)
30 continue
c
c # interpolate between pure green and green-red
ilower = iupper + 1
iupper = ilower + neach/2
irange = iupper - ilower
range = float(irange)
blue = 0.0
green = 1.0
do 40 i = ilower, iupper
red = float(i-ilower)/range
call gscr(1,i,red,green,blue)
40 continue
c
c # interpolate between green-red and pure red
ilower = iupper + 1
iupper = ncolors + 2
irange = iupper - ilower
range = float(irange)
blue = 0.0
red = 1.0
do 50 i = ilower, iupper
green = float(iupper-i)/range
call gscr(1,i,red,green,blue)
50 continue
c
return
end
c
c ----------------------------------------------------------------
c
subroutine lplotx(framecnt,line,linelen)
c
include "call.i"
c
integer framecnt, llen
character lline*50, line*50
c FIX DIMENSIONING
dimension xcord(2933),val(2933)
data lline/15HLINEPLOT LEVEL /,llen/16/
c
iadd(i,j) = loc + i - 1 + maxrows*(j-1)
c
c plot chosen variable in a straight line through domain.
c (here it's a horizontal line, use mid value for y)
c
level = 1
npts = 0
lline = "LINEPLOT "//line(1:linelen)//" LEVEL "
lentot = linelen + 9
eps = hxposs(lfine) /10.
ymid = (ylower + yupper)/2.0
c
xloc = xlower + eps
yloc = ymid
c
2 level = lfine
3 call srchlev(xloc,yloc,level,ngrd,i,j,x,y)
if (ngrd .eq. 0) then
level = level - 1
if (level .lt. 1) then
write(6,*)" couldn't find any grid in lplotx"
call clsgks
stop
endif
go to 3
endif
c
c ngrd is starting grid
c
loc = node(store1,ngrd)
maxrows = node(ndihi,ngrd) - node(ndilo,ngrd) + 3
maxcols = node(ndjhi,ngrd) - node(ndjlo,ngrd) + 3
15 npts = npts + 1
xloc = x
xcord(npts) = xloc
val(npts) = alloc(iadd(i,j))
c
c find next grid to continue with
c
if (level .eq. lfine) then
x = x + hxposs(level)
i = i + 1
if (i .lt. maxrows) go to 15
endif
c
c in general, increment just enough to possibly start next finer level
c
xloc = xloc + hxposs(level)/2.0 + eps
if (xloc .lt. xupper) go to 2
c
c the line is set. now plot
c
40 continue
valmax = -1.e10
valmin = 1.e10
do 115 i = 1, npts
valmin = amin1(valmin,val(i))
115 valmax = amax1(valmax,val(i))
if (valmax .eq. valmin) then
valdif = .05*valmin
else
valdif = .05*(valmax - valmin)
endif
valmin = valmin - valdif
valmax = valmax + valdif
c
valmax = float(int(valmax+.01) + 1)
valmin = float(int(valmin))
call set(.2,.8,.2,.8,xlower,xupper,valmin,valmax,1)
call labmod('(f4.1)','(f6.2)',4,6,1,1,0,0,0)
call gridal(4,4,2,5,1,1,5,dummy,dummy)
c call curve(xcord,val,npts)
ich = ichar('.')
do i = 1, npts
call point(xcord(i),val(i))
end do
c
c setup and print title
c
valdif = (valmax-valmin)/3.
call set(.2,.9,.2,1.,xlower,xupper,valmin,valmax+valdif,1)
xplace = xlower
yplace = valmax +valdif/2.
call pwrit(xplace,yplace,lline,lentot,2,0,-1)
framecnt = framecnt + 1
write(7,299) framecnt, "LINEPLOT",lline
299 format( "Frame# ",i3,2x,a8,2x,a50 )
call frame
c
return
end