      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
