Program contour
c                       CONTOUR  
c               A contour plotting package
c
c  Robert L. Parker, David W. Caress and Peter Salameh
c  Institute of Geophysics and Planetary Physics
c  Scripps Institution of Oceanography
c  La Jolla, CA 92093
c
c  This program was adapted from the original version of contour
c  created by Robert L. Parker (April 1984) using a version of
c  the main routine flevel altered  by Peter Salameh (March 1989).
c
c  Easily operated contouring program based upon Anderson's
c  contouring code for smooth contours of a rectangular array
c  of real values.  Program accepts commands from a limited
c  vocabulary (see menu below) that cause various actions -
c    read   - read in the array
c    levels - specify contour values
c    status - given current status of parameters
c  etc, etc.
c---------------------------------------------------------------------
c                       This is Unit A 
c  Unit A contains routines for language interpretation and data input
c  Unit B contains page arrangement, axes and other decoration
c  Unit C contains contour generation and labeling routines
c---------------------------------------------------------------------
c
c$$$$ calls decode fancy halt help init inmat major scribe split status
c
c     common variables controlling plotting
c     of axes, labels, and all options which are
c     not directly used in plotting the contours,
c     xydata, or vector plots
c
c inp     --  unit to read from screen
c iout    --  unit to write to screen
c idat    --  unit for data input from a file
c idev    --  output device number
c inf     --  data input unit - can be equal to inp or idat
c
c ms      --  number of x values of grid to read (before rotation)
c             if ms < 0 then grid reflected about vertical plane
c ns      --  number of y values of grid to read (before rotation)
c             if nx < 0 then grid reflected about horizontal plane
c m       --  number of x values of grid read
c n       --  number of y values of grid read
c klock   --  grid rotation parameter - number of clockwise quarter
c             rotations
c nskip   --  number of lines in next data file to skip before
c             starting read
c nauto   --  contour values determined automatically if nauto=1
c kindax  --  axis style flag
c             kindax=0    :  linear x and y axes
c             kindax=1    :  log x axis and linear y axis
c             kindax=2    :  linear x axis and log y axis
c             kindax=3    :  log x and y axes
c ntxt    --  number of lines of text in variable text for caption
c
c height  --  height in inches of plot
c width   --  width in inches of plot
c ax(4)   --  axis bound values
c             ax(1) : xmin
c             ax(2) : xmax
c             ax(3) : ymin
c             ax(4) : ymax
c hlet(5) --  height of labels in inches
c             hlet(1) : title height
c             hlet(2) : axis numerals height
c             hlet(3) : axis labels height
c             hlet(4) : caption height
c             hlet(5) : unspecified
c invoke  --  flag to turn on fancy option - 0=off, 1=on
c vcrit   --  if invoke=1, vcrit is drawn with heavy line
c             contours < vcrit are drawn dashed and contours
c             > vcrit are drawn solid
c
c note(20)--  contains up to 20 80 character notes to be plotted
c nnotes  --  number of notes
c leng(20)--  length in characters of each note - >0 if in axes coords
c             and <0 if in inches
c xyh(7,20)-  contains 7 note parameters for each note
c             1 : x position of arrowhead
c             2 : y position of arrowhead
c             3 : x position of note beginning
c             4 : y position of note beginning
c             5 : ht of note lettering in inches
c             6 : angle of note
c             7 : colorcode of note
c
c xyorg(2)--  origin of plot
c xyold(2)--  old origin of plot
c
c
      parameter (nulb=500)
      common /xylabs/ xlb(nulb),ylb(nulb),nlb
      common /back/ backy,box
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      character*116 format,file,title,output,xlab,ylab
      common /charco/  format,file,title,output,xlab,ylab,text,note
c
c     common variables controlling contouring
c
c delx    --  inch distance between x-grid values
c dely    --  inch distance between y-grid values
c nx      --  number of x-grid points
c ny      --  number of y-grid points
c
c qnllhi --  logical flag turns on null for values >= anulhi
c anulhi  --  if qnllhi true then array values >= anulhi are nulled
c qnlllo --  logical flag turns on null for values <= anullo
c anullo  --  if qnllhi true then array values <= anullo are nulled
c
c maxpts  --  maximum dimension of one side of grid
c maxcnt  --  maximum number of contour intervals
c maxcmp  --  maximum number of grid points in each contour line
c maxflg  --  maximum number of segments/level
c f       --  array to be contoured
c value   --  real array of contour values
c kode    --  integer array of contour line types
c kolor   --  integer array of contour line colors
c nlines  --  number of contour intervals
c iline   --  index of current contour interval
c
c scalev  --  scaling factor for contour line labels
c hlab    --  label character height in inches
c smooth  --  line smoothing factor (no. segments/radian arc)
c interp  --  interpolation mode flag - 0=linear, 1=parabolic
c nsgmnt  --  minimum number of points allowed for contour
c xsave   --  buffered points in contour
c ysave   --  buffered points in contour
c nsave   --  number of buffered points in contour
c
c flag    --  list of segments which current contour value
c	      crosses, packed 100000*j + 10*i + k
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      common /flags/ nflag, kntflg(maxflg)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
c     common variables controlling plotting of
c     xydata pairs
c
c nxyp             --  number of xydata reads made
c newp             --  newest xydata read
c nstart(maxset+1) --  starting index of each xydata series in xp and yp
c nkind(maxset)    --  plotting style of xydata series
c                      if nkind(ixy) > 1000 then
c                        data symbol id=nkind(ixy)/1000 - 1
c                        data symbol size=(nkind(ixy)-1000*(id+1))/100
c                                        in inches
c                      if nkind(ixy) < 1000 then
c                        data line style id=nkind(ixy)
c xp(maxp)         --  x data
c yp(maxp)         --  y data
c
c
      parameter(maxp=20000,maxset=5000)
      common /xyco/ newp,nxyp,nstart(maxset+1),nkind(maxset),
     $      ncolor(maxset), xp(maxp),yp(maxp)
c
c eps     --  small compared to min/max z value
c epps    --  small compared to min/max x-y values
      common /little/ eps,epps
c
c Stuff for evenly spaced contours (Command 'inte')
      common /evenly/ neven,zz1,zz2,dzz
c  Postscript control variables/
      common /post/ boxx,boxy,land,infill
c  Affine scaling of input array
      common /affine/ acoef, bcoef
c
c  Flags that level lines must be saved to disk
      common /disk/ isave,dxy(2)
c
      dimension  varble(maxcnt)
      character*4 task,  line*120, tail*116, try*12
c
      parameter (items=38)
      character*4 menu(items)
c
      data menu/'file','form','capt','titl','outp',
     $          'xlab','ylab','widt','heig','scal',
     $          'devi','skip','INTE','smoo','segm',
     $          'Colo','bord','land','    ','inte',
     $          'leve','code','colo','read','lett',
     $          'axes','xyda','orie','null','note',
     $          'stat','help','fanc','plot','stop',
     $          'affi','save','quit'/
      data try/' - try again'/, lipsis/0/
c
c
ce
      nstart(1)=1
      write(*, '(a)') ' Enter contour commands'
c
c  Read a command line.  Comment if char 1 is blank or %
c
   10 read (*, '(a4, a116)', end=3500) task, line
      if (task(1:1).eq.' ' .or. task(1:1).eq.'%') goto 10
      call split(line(1:116), tail, ntail)
c
c  Check the menu of commands
c
   11 do 12 job=1, items
        if (task .eq. menu(job)) goto 15
   12 continue
      write(*, '(/4a)') '>>> Unrecognized command <',task,'>',try
      goto 10
c
c  The task is on the menu.  Perform it
   15 goto ( 100,  200,  300,  400,  500,
     $       600,  700,  750,  750,  750,
     $       750,  750,  750,  750,  750,
     $       750, 1700, 1800,  750,  750,
     $      2050, 2050, 2050, 2400, 2500,
     $      2600, 2700, 2800, 2900, 3000,
     $      3100, 3200, 3300, 3400, 3500,
     $      3600, 3700, 3500), job                              


c
c  job=1-7.  Save various character strings
c
c  job=1: 'file'
  100 file=tail
      goto 10
c  job=2: 'format'
  200 if (ntail .eq. 0) then
         format(1:1)=' '
      else
         if (index(tail, 'b') .ne. 0) format(1:1)='b'
         if (index(tail, '3') .ne. 0) format(2:2)='3'
         if (index(tail, '2') .ne. 0) format(2:2)='2'
      endif
      goto 10
c  job=3: 'caption'
  300 text(1)=tail
      do 305 ntxt=1, 9
        if (tail(ntail:ntail) .eq. '&') then
          text(ntxt)=tail(1:ntail-1)
          read (*, '(a80)') tail
          do 302 ntail=80, 1 ,-1
            if (tail(ntail:ntail) .ne. ' ') goto 303
  302     continue
          ntail=1
  303     text(ntxt+1)=tail
         else
         goto 10
        endif
  305 continue
      ntxt=9
      goto 10
c  job=4: 'title'
  400 title=tail
      goto 10
c  job=5: 'output'
  500 output=tail
      goto 10
c  job=6: 'xlabel'
  600 xlab=tail
      goto 10
c  job=7: 'ylabel'
  700 ylab=tail
      goto 10
c
c  job=8-20.  Save various simple constants
c
  750 call decode(tail(1:ntail), varble, 3, nfound)
      if (nfound .le. 0) then
        write(*,'(4a)') '>>> ',task,
     $  ' must be followed by a number',try
        goto 10
      endif
      goto (800,900,1000,1100,1200,1300,1400,1500,
     $      1600,1700,1800,1900,2000), (job-7)
c  job=8: 'width'
  800 width=varble(1)
      goto 10
c  job=9: 'height'
  900 height=varble(1)
      goto 10
c  job=10: 'scale'
 1000 scalev=varble(1)
      goto 10
c  job=11: device
 1100 idev=varble(1)
      goto 10
c  job=12: 'skip'
 1200 nskip=varble(1)
      goto 10
c  job=13: INTErpolate
 1300 interp=varble(1)
      goto 10
c  job=14: 'smooth'
 1400 mooth=varble(1)
      goto 10
c  job=15: 'segment'
 1500 if (nint(varble(1)).ge.2) nsgmnt=nint(varble(1))
      goto 10
c  job=16: 'Color'
 1600 ikolor=varble(1)
      goto 10
c  job=17: 'border'
 1700 continue
      box=-box
      goto 10
c  job=18: 'landscape' mode
 1800 continue
      land=1
      goto 10
c  19 not used...
 1900 continue
      goto 10
c  job=20: 'interval'
c  Evenly spaced contour intervals requested
 2000 neven=nfound
      if (neven .eq. 3) then
        zz1=varble(1)
        zz2=varble(2)
        if (zz1 .ge. zz2) neven=1
        dzz=varble(3)
        nauto=1
      elseif (neven .eq. 1) then
        dzz=varble(1)
        nauto=1
      else
        neven=0
      endif
      goto 10
c
c  job=21, 22, 23: 'level', 'code', 'color'
c  Save contour levels or their code integers
c
c  Decode ellipsis ... in levels and code cmds
 2050 if (job.eq. 21 .or. job.eq. 22) then
        lipsis=index(tail,'...')
        if (lipsis.gt. 0) tail(lipsis:lipsis+2) = '   '
        if (lipsis.gt. 0 .and. job.eq. 22) tail(lipsis:lipsis+2)='999'
      endif
      call decode(tail(1:ntail), varble, maxcnt, nfound)
      nvalus=abs(nfound)
      if (nfound .eq. 0 .and. job .eq. 21) then
         nlines=0
         nauto=1
         goto 10
      endif
      if (nfound .eq. 0 .and. job .eq. 22) then
         call fancy(-1)
         goto 10
      endif
      if (nfound .eq. 0 .and. job .eq. 23) then
         do 2070 i=1, maxcnt
 2070      kolor(i)=0
         goto 10
      endif
      if (nfound .le. 0) then
         write(*,'(a)')'>>> Unreadable list.  Repeat the command'
         goto 10
      endif
c  No errors in list - continue reading from the  terminal
 2080 read (*, '(a120)') line
      call decode(line, varble(nvalus+1), maxcnt-nvalus, nfound)
      nvalus=nvalus + abs(nfound)
      if (nfound .gt. 0) goto 2080
c  An error in reading is detected.  The list is over
c  Copy list to proper variables
      if (job .eq. 22) call fancy(-1)
      do 2090 i=1, nvalus
        if (job .eq. 21) value(i)=varble(i)
        if (job .eq. 22) kode(i) =varble(i)
        if (job .eq. 23) kolor(i)=varble(i)
 2090 continue
      if (nfound .lt. 0) write(*, '(a/i4,a)')
     $'>>> Unable to read the complete list',
     $ nvalus,' values have been read'
      if (job.eq. 21) then
c  If levels contains ... continue with last increment
         if (lipsis.gt. 0 .and. nvalus.gt.1) then
           vmax=value(nvalus)
           if (ntail.gt. lipsis+2) nvalus=nvalus-1
           idel = max(2,nvalus)
           delta = value(idel)-value(idel-1)
           if (ntail.eq. lipsis+2) vmax=vmax+10*delta
           if (nvalus.eq. 1 .and. ntail.gt. lipsis+2) delta=delta/10
           if (delta.le. 0.0) write(*,'(a)')
     $    '>>> Positive increment required in levels with ...'
           do 2091 i=nvalus+1, maxcnt
             value(i) = value(i-1) + delta
             if (value(i).ge.vmax .or. delta.le. 0.0) goto 2092
 2091      continue
 2092      nvalus=i
         endif
         nlines=nvalus
         nauto =0
         neven =0
c  Fill kode with repeating pattern
      elseif (job.eq. 22 .and. kode(nvalus).eq. 999) then
         k1=nvalus - 1
         do 2095 k=k1, maxcnt-k1, k1
           do 2094 j=1, k1
             kode(k+j)=kode(j)
 2094      continue
 2095    continue
      endif
      if (nfound.lt. 0 .or. line.eq. ' ') goto 10
      call split(line(5:120), tail, ntail)
      task=line(1:4)
      goto 11
c
c  job=24: 'read'
c  Read the matrix.  This is a big job that is done by  inpmat
c
 2400 call decode(tail(1:ntail), varble, 3, nfound)
      if (nfound .eq. 0) then
        write(*,'(2a)')
     $ '>>> read must be followed by 1,2 or 3 numbers',try
        goto 10
      endif
c
      ms=varble(1)
      ns=0
      if (nfound .ge. 2) ns=varble(2)
      if (nfound .eq. 2) klock=0
      if (nfound .eq. 3) klock=varble(3)
      call inpmat
      goto 10
c
c  job=25: 'letter'
c  Change height of (1) contour labels (2) title or caption
c
 2500 call decode(tail(1:ntail), varble, 5, nfound)
      if (nfound .lt. 1) then
         write(*,'(2a)')
     $   '>>> letter must be followed by 1-5 numbers',try
         goto 10
      endif
      hlab=varble(1)
      do 2550 l=2, nfound
        hlet(l-1)=varble(l)
 2550 continue
      goto 10
c
c  job=26: 'axes'
 2600 call decode(tail(1:ntail), varble, 5, nfound)
      if (nfound .eq. 0) then
        do 2610 l=1, 4
          ax(l)=0.0
 2610   continue
        goto 10
      endif
      do 2620 l=1, min(4, nfound)
        ax(l)=varble(l)
 2620 continue
      if (nfound .eq. 5) kindax=varble(5)
      goto 10
c
c  job=27: 'xydata'
c
 2700 call decode(tail(1:ntail), varble, 2, nfound)
      nxtk=min(maxset, nxyp+1)
      if (nfound .eq. 1) nkind(nxtk)=abs(varble(1))
      if (nfound .eq. 2) nkind(nxtk)=100*varble(2)+1000*(varble(1)+1.0)
      ncolor(nxtk)=ikolor
      goto 10
c
c  job=28: transform
c  Re-orient array after reading (plants numbers in /orient/
c
 2800 call group(tail(1:ntail))
      goto 10
c
c  job=29: 'null'
c
 2900 call decode(tail(1:ntail), varble, 2, nfound)
      if (nfound.eq.0) then
        qnlllo=.FALSE.
        qnllhi=.FALSE.
      else if (nfound.eq.2.and.varble(2).gt.varble(1)) then
        qnlllo=.TRUE.
        qnllhi=.TRUE.
        anullo=varble(1)
        anulhi=varble(2)
      else if (nfound.eq.2.and.varble(2).le.varble(1)) then
        write(*,'(2a)')'>>> 2nd null value must exceed 1st',try
      else
       write(*,'(2a)')
     $  '>>> null must be followed by 0 or 2 values',try
      endif
      goto 10
c
c  job=30: 'note'
c
 3000 if (tail(1:1).eq. ' ') then
        nnotes=0
        goto 10
      endif
      n1 =min(20, nnotes+1)
      if (n1.eq.1) xyh(5,1)=hlab
      xyh(5,n1)=xyh(5, max(n1-1,1))
      xyh(7,n1)=ikolor
      call scribe(ntail, tail(1:ntail), xyh(1,n1), leng(n1), note(n1))
      if (leng(n1) .eq. -1) then
        write(*,'(a)')'>>> Error reading note - note ignored'
        goto 10
      endif
      if (xyh(5,n1).gt.10) xyh(5,n1)=hlet(5)
      if (xyh(5,n1).lt.-10) xyh(5,n1)=-hlet(5)
      nnotes=n1
      goto 10
c
c  job=31: 'status'
c
 3100 call status
      goto 10
c
c  job=32  Help
c
c3200 call help(job,menu,items,tail,ntail,iout)
 3200 write(*,'(a)')' Help not provided - sorry'
      goto 10
c
c  job=33: 'fancy'
c
 3300 call decode(tail(1:ntail), varble, 1, nfound)
      if (nfound .ge. 1) vcrit=varble(1)
      invoke=1
      if (nfound .le. 0) call fancy(-1)
      goto 10
c
c  job=34: 'plot'
c  Contour the array.  This task is left to a major routine
c
 3400 call decode(tail(1:ntail), varble, 2, nfound)
      xyorg(1)=0.0
      xyorg(2)=1000.0
      do 3410 j=1, nfound
        xyorg(j)=varble(j)
 3410 continue
c
      call major
c
      write(*,'(/a)')' Contouring completed.  ',
     $' Enter further commands'
      nlb=0
      goto 10
c
c  job=35: 'stop' or 'quit'
c
 3500 call init(idev, ' ')
      call halt('Normal termination')
      stop
c
c  job=36; 'affine'
 3600 call decode(tail(1:ntail), varble, 2, nfound)
      if (nfound .eq. 0) then
        acoef=1.0
        bcoef=0.0
      elseif (nfound .eq. 2) then
        acoef=varble(1)
        bcoef=varble(2)
      else
        write(*,'(a)')'>>> Error reading affine - values unchanged'
      endif
      goto 10
c
c  job=37; 'save'
 3700 isave=1
      if (tail(1:1).eq. '-') then
        isave=0
      elseif (tail(1:1).ne. ' ') then
        open(unit=12, file=tail)
      endif
      goto 10
c
      end
c______________________________________________________________
      blockdata dfault
c  Sets various parameters into common for the user's convenience
c
      logical qnllhi, qnlllo
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
      common /orient/ irient,maphr(3)
      common /affine/ acoef, bcoef
c
      data scalev /1.0/
      data hlab /0.08/
      data mooth /15/
      data interp /1/
      data nsgmnt /2/
      data kode /maxcnt*0/
      data kolor /maxcnt*1/
      data nlines /0/
c
      data qnllhi /.FALSE./
      data qnlllo /.FALSE./
      data anulhi,anullo /2*0.0/
c
      data irient/0/
c
      data acoef, bcoef/1,0/
      end
c______________________________________________________________
      blockdata cfault
c$$$$ calls nothing
c     Loads basic common blocks
c
      common /back/ backy,box
c
      common /evenly/ neven,zz1,zz2,dzz
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      character*116 format,file,title,output,xlab,ylab
      common /charco/  format,file,title,output,xlab,ylab,text,note
c
      parameter(maxp=20000,maxset=5000)
*     common /xyco/ newp,nxyp,nstart(maxset+1),nkind(maxset),
*    $      ncolor(maxset), xp(maxp),yp(maxp)
      common /disk/ isave,dxy(2)
      parameter (nulb=500)
      common /xylabs/ xlb(nulb),ylb(nulb),nlb
c
      data nlb/0/
c /back/
      data box /1.0/
c /evenly/
      data neven/0/
c /plotco/
      data idev,ikolor /-1,1/
      data ms,ns,m,n/4*0/
      data klock/0/
      data nskip,nauto /0,1/
      data kindax,ntxt,invoke /0,0,-1/
      data height,width /0.0,0.0/
      data ax /4*0.0/
      data hlet /5*0.08/
      data vcrit/0.0/
      data xyorg /0.0,1000.0/
      data xyold /0.0,0.0/
      data nnotes /0/
c /charco/
      data format,file,title,output,xlab,ylab  /6*'  '/
      data text(1) /'  '/
c  /xyco/
*     data nstart(1)/1/
c /disk/
      data isave/0/
c
      end
c______________________________________________________________
      subroutine halt(messag)
c$$$$ calls nothing
c  Program prints the message  messag  then stops
      character*(*) messag
c
ce
      write(*, '(/1x, a)') messag
      stop
      end
c_____________________________________________________________________
      subroutine parse(l1, ltxt, text, lwd, word)
c$$$$ calls nothing
c  From the string  text,  beginning at character l1, finds the
c  next word delimited by spaces and inserts into  word.  The
c  length of  text  is  ltxt  and the length of  word  is  lwd.
c  l1  is set past the current word.
c  If text(l1:ltxt) is all blank, returns lwd=0.
      character *(*) text,word
c
c  Skip leading blanks
ce
      do 1100 i=l1, ltxt
        if (text(i:i) .ne. ' ') goto 1150
 1100 continue
c  Only blanks in the field
      lwd=0
      l1=0
      return
 1150 continue
c
c  Locate end of the word
      j=index(text(i:ltxt), ' ') + i - 1
c  No blanks found - whole field is new word
      if (j .eq. 0) j=ltxt+1
      word=text(i:j-1)
      lwd=j - i
      l1=j + 1
      return
      end
c______________________________________________________________
      subroutine split(image, tail, ntail)
c$$$$ calls nothing
c  Discards leading characters until blank.  Then discards blanks until
c  first nonblank.  Then returns the rest in tail.   ntail  is the
c  length of the string  tail  when trailing blanks are trimmed.
      character*116 image, tail
c
ce
      ibegin=0
      do 1100 l=1, 116
      if (image(l:l) .eq. ' ') ibegin=1
      if (image(l:l).ne.' ' .and. ibegin.eq.1) then
        do 1050 k=116, l, -1
         if (image(k:k) .ne. ' ') then
           ntail=k - l + 1
           tail=image(l:k)
           return
         endif
 1050   continue
      endif
 1100 continue
      ntail=1
      tail=' '
      return
      end
c______________________________________________________________
      subroutine decode(char, values, nwant, nfound)
c$$$$ calls nothing
c  Evaluates numbers in the string  char.  nwant  numbers are expected,
c  nfound  are actually found in the character string.  If an error is
c  discovered  nfound=-n,  where  n  is the number of numbers properly
c  decoded.
      dimension values(*)
      character*(*) char, local*40
c
ce
      kn=len(char)
      k1=1
c  Up to  nwant  numbers are sought
      do 1800 nval=1, nwant
        do 1100 k=k1, kn
          if (char(k:k) .ne. ' ') goto 1200
 1100   continue
        nfound=nval-1
        return
 1200   do 1300 l=k, kn
          if (char(l:l).eq. ',' .or. char(l:l) .eq. ' ') goto 1500
 1300   continue
 1500   k2=l+1
        local=char(k:l-1)
        read (local, '(f40.0)', err=1900) values(nval)
 1800 k1=k2
      nval=1 - nwant
 1900 nfound=1 - nval
      return
      end
c______________________________________________________________
      subroutine status
c$$$$ calls assess fancy
c  Reports the status of the various quantities in memory
c
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      character*116 format,file,title,output,xlab,ylab
      common /charco/  format,file,title,output,xlab,ylab,text,note
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      parameter(maxp=20000,maxset=5000)
      common /xyco/ newp,nxyp,nstart(maxset+1),nkind(maxset),
     $      ncolor(maxset), xp(maxp),yp(maxp)
c
c  Report on current array limitations
ce
      write(*,'(/(a,i6))')
     $'             Contour revision  November ',2005,
     $'    Largest array dimension (MAXPTS)',maxpts,
     $'       Max number of levels (MAXCNT)',maxcnt,
     $' Max number of x-y data pairs (MAXP)',maxp,
     $'   Max number of x-y series (MAXSET)',maxset,
     $' '
c
c  Report on the various strings
      if (format(1:2).ne.' ') write(*, '(a)')'  format:', format
      if (file(1:1)  .ne.' ') write(*, '(a)')'    file:', file
      if (title(1:1) .ne.' ') write(*, '(a)')'   title:', title
      if (ntxt .gt. 0)        write(*, '(a)')' caption:', text(1)
      if (output(1:1).ne.' ') write(*, '(a)')'  output:', output
      if (xlab(1:1)  .ne.' ') write(*, '(a)')'  xlabel:', xlab
      if (ylab(1:1)  .ne.' ') write(*, '(a)')'  ylabel:', ylab
      if (nnotes.gt.0) write(*,'(a,i3)')' Number of notes ',nnotes
c
c  Report on numerical axis annotation
      if (ax(1) .ne. ax(2))
     $  write(*, '(a,1p,2g12.4)') ' x-axis interval', ax(1),ax(2)
      if (ax(3) .ne. ax(4))
     $  write(*, '(a,1p,2g12.4)') ' y-axis interval', ax(3),ax(4)
c
c  Report on the x-y data series
      if (nxyp .gt. 0) write(*,'(a,i3,a)')
     $  ' There are',nxyp, ' sets of x-y data series'
c
c  Report on the array itself
      if (m .eq. 0) then
         write(*, '(a)') ' No array has been entered yet'
      else
         write(*, '(a,t40,2i5)') ' Array is stored as ', m,n,
     $                         ' Number of quarter turns', klock
         fmin=f(1,1)
         fmax=fmin
         do 1100 j=1, n
           do 1100 i=1, m
             fmin=min(f(i,j), fmin)
             fmax=max(f(i,j), fmax)
 1100    continue
         write(*, '(a,t40,1p,2g13.5)') ' Min, max values',fmin,fmax
      endif
c
c  Report on contour level status
      if (nauto .eq. 1) write(*, '(a)')
     $    ' No contours have been supplied by user'
      if (invoke .eq. 1) write(*, '(a,1p,g13.5)') 
     $    ' Fancy option enabled with ',vcrit
c
      if (nauto .eq. 1 .and. m .gt. 0) call assess
      if (nlines .gt. 0) then
        call fancy(invoke)
        k1=min(4, (4+nlines)/5)
        k1=min(nlines, 5*((4+(nlines+k1-1)/k1 )/5))
        write(*, '(/a)')' Codes, colors, levels of contours:'
        do 1102 k=1, k1
          write(*,'(1p,4(2i3,g13.5))')
     $                 (kode(j),kolor(j),value(j),j=k,nlines,k1)
 1102   continue
        write(*,*)' '
      endif
      if (nauto .eq. 1) nlines=0
      if (scalev .ne. 1.0) write(*, '(a,1p,g15.5)')
     $   ' Contour labels are scaled by a factor',scalev
      write(*,'(a,i5)') ' Contour smoothing factor',mooth
      if (interp.eq.0) then
        write(*,'(a)')' Contour interpolation:  linear'
      else
        write(*,'(a)')' Contour interpolation:  parabolic'
      endif
      write(*,'(a,i5)')' Minimum segment size of contours',nsgmnt
      if (qnlllo) write(*,'(a,1p,g15.5)')
     $    ' Null enabled for values <= ',anullo
      if (qnllhi) write(*,'(a,1p,g15.5)')
     $    ' Null enabled for values >= ',anulhi
c
c  Report on the sizes of various things
      write(*,'(/a,t35,1p,2g12.4)')
     $    ' Height and width parameters',height,width
      if (m .gt. 1 .and. n .gt. 1) then
         h=height
         w=width
         if (h .le. 0 .and. w .eq. 0) then
           if (m .ge. n) w=6
           if (n .gt. m) h=6
         endif
         if (h .le. 0) h=w*(n-1)/(m-1.0)
         if (w .le. 0) w=h*(m-1)/(n-1.0)
         write(*, '(a,t35,2g12.4)') ' Inferred height and width',h,w
      endif
      return
      end
c______________________________________________________________
      subroutine inpmat
c$$$$ calls inplin matrab matraf
c  Routine for reading matrices and re-orienting them on the page.
c  anyone who has used a slide projector knows there are 8 ways of
c  viewing a slide, 7 of them wrong.  Similarly, there are 8 different
c  ways to view a rectangular array of data.  This routine together
c  with  matran  allows the user to get an array into the desired
c  orientation.
c   It is assumed that there is a file containing an array written
c  it to your specifications.  (Actually the routine just does one
c  clever read.)
c
c  ms, ns  dimensions of array as it was written
c  klock   number of clockwise quarter-turns
c  m, n    array dimensions of matrix as actually used in program
c
      parameter(inp=5,iout=6,idat=1)
      parameter (maxpts=350)
c
      character*116 format,file,title,output,xlab,ylab
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      common /charco/  format,file,title,output,xlab,ylab,text,note
      common /orient/ irient,maphr(3)
c
      character*116 name
      save name
c
      data name/'   '/
c
c  If both  name and file begin with a blank there's no filename!
c
ce
      if (name(1:1) .eq. ' ' .and. file(1:1) .eq. ' ') then
        write(*, '(a)')'>>> There is no filename for the array'
        return
      endif
c
c  Check if this is a new filename.  If it is open the file
c  name is the locally saved filename
      if (file .ne. name) then
         if (file(1:1) .eq. ' ') file=name
         if (name(1:1) .ne. ' ' .and. inf.ne.inp) then
           rewind (unit=inf)
           close  (unit=inf)
         endif
         if (file(1:1) .eq. '*') then
           inf=inp
           name='*'
         else
           inf=idat
           if (format(1:1) .eq. ' ') then
             open (unit=inf, file=file,status='OLD', err=1000)
             do 100 iskip=1,nskip
               read(inf,*)
  100        continue
           else if (format(1:1) .eq. 'b') then
             open (unit=inf, file=file,
     $       status='OLD', err=1000, form='UNFORMATTED')
           endif
           name=file
         endif
      endif
c
c  x-y data read here  ns=0, ms >= 0
c
      if (ns .eq. 0 .and. ms .ge. 0) then
        call inplin
        return
      endif
c
c  If orient command has been used, make equivalent [ms,ns,klock]
c  from values in /orient/
c
      if (irient .eq. 1) then
         ms=sign(ms, maphr(1))
         ns=abs(ns)
         klock=maphr(2)
      endif
c
c  Impossible array dimensions trapped
c
      if (abs(ms) .le. 1  .or. abs(ns) .le. 1 .or.
     $   abs(ms).gt.maxpts .or. abs(ns) .gt.maxpts) then
         write(*, '(a,2i5/a)')'>>> Improper array dimensions', ms,ns,
     $   '>>> Array cannot be contoured.  Please enter another'
         if(abs(ms).gt.maxpts.or.abs(ns).gt.maxpts) write(*,'(a,i5)')
     $   '>>> Array dimension exceeds largest permitted, ie',maxpts
         m=0
         n=0
         return
      endif
c
c
c  Determine the array dimensions as they really are
c
      if (mod(klock,2) .eq. 0) then
        m=abs(ms)
        n=abs(ns)
      else
        m=abs(ns)
        n=abs(ms)
      endif
c
c  Code the sign configuration into a single integer
c
      jsign=(5 - sign(2,ns) - sign(1,ms))/2
c
c  Now read the thing in,  either formatted or binary
c
      if (format(1:1) .eq. ' ') call matraf(jsign)
      if (format(1:1) .eq. 'b') call matrab(jsign)
c
      return
c
c  Can't open the file for some reason
 1000 write(*,'(a)')'>>> Nonexistent file or read permission denied'
      file='      '
      return
      end
c______________________________________________________________
      subroutine matraf(jsign)
c$$$$ calls nothing
c  Routine actually to read formatted diskfile with * format
c  it is assumed the file was written by rows with a statement like
c     write(iunit, format) ((a(i,j), j=1, n), i=1, m)
c  See comments in getmat.
c
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      common /affine/ acoef, bcoef
c
      parameter (inp=5)
      dimension         lookup(4,4)
c
      data lookup /1,2,3,4,5,6,7,8,
     $             4,3,2,1,8,7,6,5/
c
c  Transform the number of clockwise rotations to one of  0,1,2,3
ce
      klick=mod(mod(klock,4) + 4, 4)
c
c  All possible orientations reduce to one of the 8 below.  lookup
c  has this coded into it
c
      igo=lookup(jsign, klick+1)
      goto (1010, 1020, 1030, 1040, 1050, 1060, 1070, 1080), igo
c
c  No transformations.  Read  a  as written
 1010 read (inf,*, err=2010,end=2000) ((f(i,j), j=1,n), i=1, m)
      goto 1500
c
c  v - reflect in vertical plane
 1020 read (inf,*, err=2010,end=2000) ((f(i,j), j=1,n), i=m,1,-1)
      goto 1500
c
c  h - reflect in horizontal plane
 1030 read (inf,*, err=2010,end=2000) ((f(i,j), j=n,1,-1), i=1,m)
      goto 1500
c
c  vh=r**2  - reflect in v and h= 2 clockwise 90 deg turns
 1040 read (inf,*, err=2010,end=2000) ((f(i,j), j=n,1,-1), i=m,1,-1)
      goto 1500
c
c  r - 1 clockwise 90 deg turn
 1050 read (inf,*, err=2010,end=2000) ((f(i,j), i=1,m), j=n,1,-1)
      goto 1500
c
c  rv - v reflection followed by r
 1060 read (inf,*, err=2010,end=2000) ((f(i,j), i=m,1,-1), j=n,1,-1)
      goto 1500
c
c  rh - h reflection followed by r=transpose of a
 1070 read (inf,*, err=2010,end=2000) ((f(i,j), i=1,m), j=1,n)
      goto 1500
c
c  r**3=r**-1 - 3 clockwise 90 deg turns=one anticlockwise
 1080 read (inf,*, err=2010,end=2000) ((f(i,j), i=m,1,-1), j=1,n)
c
c  Transform the array values according to the affine transformation
 1500 if (acoef.eq.1 .and. bcoef.eq.0) return
      do 1550 i=1, m
        do 1510 j=1, n
          f(i,j)=acoef * f(i,j) + bcoef
 1510   continue
 1550 continue
      return
c
c  Error section.  Improper dimensions have been supplied
c
 2000 write(*,200)
 200  format(/'>>> Unexpected endfile encountered.  Improper'
     $,' array size has been given'/' file rewound')
      goto 2015
c
 2010 write(*,'(a)') '>>> Error reading array data ********'
 2015 if (inf .ne. inp) rewind (unit=inf)
      m=0
      return
      end
c______________________________________________________________
      subroutine matrab(jsign)
c$$$$ calls nothing
c  Routine actually to read binary diskfile.
c  It is assumed the file was written by rows with a statement like
c        do 1000 i=1, m
c   1000 write(iunit) (a(i,j), j=1, n)
c  in an ftn77  compiled program
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      common /affine/ acoef, bcoef
c
      dimension         lookup(4,4)
c
      data lookup /1,2,3,4,5,6,7,8,
     $             4,3,2,1,8,7,6,5/
c
c  Transform the number of clockwise rotations to one of  0,1,2,3
ce
      klick=mod(mod(klock,4) + 4, 4)
c
c  All possible orientations reduce to one of the 8 below.  lookup
c  has this coded into it
c
      igo=lookup(jsign, klick+1)
      goto (1010, 1020, 1030, 1040, 1050, 1060, 1070, 1080), igo
c
c  No transformations.  Read  a  as written
 1010 do 1015 i=1, m
 1015 read (inf, end=2000) (f(i,j), j=1,n)
      goto 1500
c
c  v - reflect in vertical plane
 1020 do 1025 i=m,1,-1
 1025 read (inf, end=2000) (f(i,j), j=1,n)
      goto 1500
c
c  h - reflect in horizontal plane
 1030 do 1035 i=1,m
 1035 read (inf, end=2000) (f(i,j), j=n,1,-1)
      goto 1500
c
c  vh=r**2  - reflect in v and h= 2 clockwise 90 deg turns
 1040 do 1045 i=m,1,-1
 1045 read (inf, end=2000) (f(i,j), j=n,1,-1)
      goto 1500
c
c  r - 1 clockwise 90 deg turn
 1050 do 1055 j=n,1,-1
 1055 read (inf, end=2000) (f(i,j), i=1,m)
      goto 1500
c
c  rv - v reflection followed by r
 1060 do 1065 j=n,1,-1
 1065 read (inf, end=2000) (f(i,j), i=m,1,-1)
      goto 1500
c
c  rh - h reflection followed by r=transpose of a
 1070 do 1075 j=1,n
 1075 read (inf, end=2000) (f(i,j), i=1,m)
      goto 1500
c  r**3=r**-1 - 3 clockwise 90 deg turns=one anticlockwise
 1080 do 1085 j=1,n
 1085 read (inf, end=2000) (f(i,j), i=m,1,-1)
c
c  Transform the array values according to the affine transformation
 1500 if (acoef.eq.1 .and. bcoef.eq.0) return
      do 1550 i=1, m
        do 1510 j=1, n
          f(i,j)=acoef * f(i,j) + bcoef
 1510   continue
 1550 continue
      return
c
c  Error section.  Improper dimensions have been supplied
c
 2000 write(*,200)
 200  format(/'>>> Unexpected endfile encountered.  Improper'
     $,' array size has been given'/' File rewound')
      rewind (unit=inf)
      m=0
      return
      end
c______________________________________________________________
      subroutine inplin
c$$$$ calls nothing
c  Reads in x-y pair data, scaling and clipping according to the current
c  axis window and using kode(1) for type of line.
c
      parameter(inp=5,iout=6,idat=1)
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      character*116 format,file,title,output,xlab,ylab
      common /charco/  format,file,title,output,xlab,ylab,text,note
c
      parameter(maxp=20000,maxset=5000)
      common /xyco/ newp,nxyp,nstart(maxset+1),nkind(maxset),
     $      ncolor(maxset), xp(maxp),yp(maxp)
c
c  Skip if axes are not specified
ce
      if ((ax(1).eq.ax(2) .or. ax(3).eq.ax(4)) .and. ms.gt.0) then
        write(*,'(a)')
     $  '>>> You must invoke  axes  before reading x-y data'
        return
      endif
c
c  Purge x-y data series if ms=0
      if (ms .eq. 0) then
        if (nxyp .ne. 0) write(*,'(a)') ' Old x-y data purged'
        nkind(1)=nkind(nxyp+1)
        nxyp=0
        newp=0
        return
      endif
c
c     write(*,'(a,a)') ' Reading x-y data pairs from diskfile ',file
      nxyp=newp
      n1=nstart(nxyp+1)
      n2=min(maxp, n1 + ms -1)
      ncolor(nxyp+1)=ikolor
      if (n2 .lt. n1+ms-1) write(*,'(a,i5)')
     $' Input truncated. Total number of points must be less than ',maxp
c  Reads data PAIRS from input
      if (format(2:2) .ne. '3') then
        if (inf .ne. inp) then
          if (format(1:1).eq.'b') read (inf,end=2000)
     $    (xp(i),yp(i),i=n1,n2)
          if (format(1:1).ne.'b') read (inf,*,end=2000)
     $    (xp(i),yp(i),i=n1,n2)
        else
          read (inf,*,err=2010)(xp(i),yp(i),i=n1,n2)
        endif
c  Reads data x,y and a dummy from input.  Discards the dummy
      else
        if (inf .ne. inp) then
          if (format(1:1).eq.'b') read (inf,end=2000)  
     $    (xp(i),yp(i),skip,i=n1,n2)
          if (format(1:1).ne.'b') read (inf,*,end=2000)
     $    (xp(i),yp(i),skip,i=n1,n2)
        else
          read (inf,*,err=2010)(xp(i),yp(i), skip, i=n1,n2)
        endif
      endif
      nxyp=min(maxset, nxyp+1)
      newp=nxyp
      nstart(nxyp+1)=n2 + 1
      a=1.0/(ax(2) - ax(1))
      if (kindax.eq.1 .or. kindax.eq.3) a=1.0/log(ax(2)/ax(1))
      b=1.0/(ax(4) - ax(3))
      if (kindax.eq.2 .or. kindax.eq.3) b=1.0/log(ax(4)/ax(3))
      xmx=max(ax(2), ax(1))
      xmn=min(ax(2), ax(1))
      ymx=max(ax(4), ax(3))
      ymn=min(ax(4), ax(3))
c  Scale coordinates.  Points out of range set xp=999
      do 1100 i=n1, n2
        if  (xp(i) .lt. xmn .or. xp(i) .gt. xmx
     $  .or. yp(i) .lt. ymn .or. yp(i) .gt. ymx) then
          xp(i)=999.0
          yp(i)=999.0
        else
c  Scale logarithmically if log axis
          if (kindax.eq.0 .or. kindax.eq.2) then
            xp(i)=(xp(i)-ax(1))*a
          else
            xp(i)=log(xp(i)/ax(1))*a
          endif
          if (kindax.eq.0 .or. kindax.eq.1) then
            yp(i)=(yp(i)-ax(3))*b
          else
            yp(i)=log(yp(i)/ax(3))*b
          endif
        endif
 1100 continue
      return
c
c  Error in read
 2000 write(*,'(a)') '>>> Premature EOF while reading x-y line data',
     $' Data ignored and file rewound'
      goto 2015
 2010 write(*,'(a)') '>>> Error in x-y data list.  Data ignored'
 2015 if (inf .ne. inp) rewind (unit=inf)
      return
      end
c______________________________________________________________
      subroutine group(string)
c$$$$ calls nothing
c  Given a string representing repeated application on an array of 
c  the actions:
c    R rotate 90 deg clockwise; V reflect in a verical plane
c    H reflect in a horizontal plane; T transpose,
c  finds an integer representing one of the states of the array
c    1=unaltered; 2=V; 3=H; 4=RR; 5=R; 6=T; 7=RH; 8=RRR.
c  This in turn is decoded into maphr in /orient/ for use by inplin
c
c  The order of operations is as one reads - thus RTV means first
c  rotate array once, then transpose the result, then reflect the
c  result in a vertical plane.
      common /orient/ irient,maphr(3)
      dimension multab(8,8), map(8,2)
      character *(*) string
c  This is the group multiplication table
      data ((multab(i,j),j=1,8),i=1,8)/
     $1,2,3,4,5,6,7,8,  2,1,4,3,7,8,5,6,  3,4,1,2,6,5,8,7,
     $4,3,2,1,8,7,6,5,  5,6,7,8,4,3,2,1,  6,5,8,7,2,1,4,3,
     $7,8,5,6,3,4,1,2,  8,7,6,5,1,2,3,4/
      data ((map(i,j),j=1,2),i=1,8)/
     $1,1, -1,1, -1,-1, 1,-1, 1,2, -1,2, -1,0, 1,0/
c
c  Orientation given as one of 8 choices: n v h u c t k a
      igroup=index('nvhuctka',string(1:1))
      if (igroup.gt. 0) goto 1600 
c  Orientation found by combination of 4 operations: VHRT 
      igroup=1
      l=len(string)
      do 1500 i=1, l
        k=index('.VH.RT',string(i:i))
        if(k .ge. 2) igroup=multab(igroup,k)
 1500 continue
c
c  Map orientation into [m, n, r] instructions
 1600 maphr(1)=map(igroup,1)
      maphr(2)=map(igroup,2)
      irient=1
      return
      end
c_________________________________________________________________
c---------------------------------------------------------------------
c                       This is Unit B 
c  Unit A contains routines for language interpretation and data input.
c  Unit B contains page arrangement, axes and other decoration.
c  Unit C contains contour generation and labeling routines.
c---------------------------------------------------------------------
      subroutine major
c$$$$ calls assess axes captio fancy flevel glyph init justy 
c$$$$ letter, line, newpn, plot, pnotes, remark
c  The major routine for contouring the array  f.
c  This routine positions the plot on the paper, opens the plotfile,
c  invokes the contouring routine  flevel,  plots titles, boxes etc.
c
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      character*116 format,file,title,output,xlab,ylab
      common /charco/  format,file,title,output,xlab,ylab,text,note
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      parameter(maxp=20000,maxset=5000)
      common /xyco/ newp,nxyp,nstart(maxset+1),nkind(maxset),
     $      ncolor(maxset), xp(maxp),yp(maxp)
c
      common /back/ backy,box
      common /disk/ isave,dxy(2)
      save oy,kase
      dimension s(4)
c
      data oy/0.0/, kase/0/
      data iup,idn /3,2/
c
ce
      kase=1 + kase
      if (kase .eq. 1) then
c  PostScript default plotfile name
        if (output(1:1) .eq. ' ') call init(idev, 'mypost')
        if (output(1:1) .ne. ' ') call init(idev, output)
      endif
c
c  No array has been provided to contour. Abandon ship
c
      if (m .eq. 0 .and. nxyp .eq. 0) then
         write(*, '(/a)')'>>> Nothing to plot - read an array'
         return
      endif
c
c  Automatic provision of levels if none is supplied
c
      if (nlines .eq. 0) nauto=1
      if (n .gt. 0) then
        if (nauto .eq. 1) call assess
        call fancy(invoke)
        k1=min(4, (4+nlines)/5)
        k1=min(nlines, 5*((4+(nlines+k1-1)/k1 )/5))
        write(*, '(/a)')' Codes, colors, levels of contours:'
        do 1020 k=1, k1
          write(*,'(1p,4(2i3,g13.5))')
     $                 (kode(j),kolor(j),value(j),j=k,nlines,k1)
 1020   continue
        write(*,*)' '
      endif
c
c
c  Decide on the size of the graph
c
      h=height
      w=width
      if (h .le. 0 .and. w .le. 0) then
         if (m .ge. n) w=6.0
         if (m .lt. n) h=6.0
      endif
      if (h .le. 0) h=w*(n-1)/(m-1.0)
      if (w .le. 0) w=h*(m-1)/(n-1.0)
      delx=w/(m-1.)
      dely=h/(n-1.)
      dxy(1)=w
      dxy(2)=h
      nx=m
      ny=n
c
c  Set drawing pen or color
      call remark('Begin a new contour map')
      call newpn(ikolor)
c
c  Set the origin
c
c  If the origin of this plot is explicitly set, move
c  it there and cancel the origin shift after previous plot.
      if (xyorg(2) .ne. 1000.0) then
        call plot(xyorg(1)-xyold(1), xyorg(2)-xyold(2), -3)
      else
c  No explicit origin requested.  Set standard position
        ox=0.5
c       oy=7.0 - h -- old position at page bottom
        oy=10.0-h
c  If there is a ylabel or y axis move origin to the right.
c  If there is a title move origin down.
        if (ylab(1:1) .ne. ' ')  ox=ox + 2.5*hlet(3)
        if (ax(3) .ne. ax(4))    ox=ox + 7.0*hlet(2)
        if (title(1:1) .ne. ' ') oy=oy - 2.0*hlet(1)
        call plot (ox, oy, -3)
      endif
c
c  Plot axes, labels captions etc
c
c  Justy forces all fonts to initial caption's unless explicitly
c  changed
      call justy(s, hlet(4), text(1))
c  If there is to be an axis with numbers plot it. The justy call forces
c  same font in axis as the last thing in the x label if any.
      dy=0.0
      if (ax(1) .ne. ax(2)) then
        dy=-2.0*hlet(2)
        call justy(s, hlet(3), xlab)
        call axes(kindax, hlet(2), w, ax, 1)
      else
        if (box .gt. 0) then
          call plot(w,   0.0, iup)
          call plot(0.0, 0.0, idn)
        endif
      endif
c  Plot x axis label if there is one
      if (xlab(1:1) .ne. ' ') then
        dy=dy - 2.0*hlet(3)
        call justy(s, hlet(3), xlab)
        call letter(0.5*w-s(2), dy, hlet(3), 0.0, xlab)
      endif
c  Plot y axis or vertical frame line
      backy=0.0
      call justy(s, hlet(3), ylab)
      if (ax(3) .ne. ax(4)) then
        call axes(kindax, hlet(2), h, ax(3), 2)
      else
        if (box .gt. 0) then
          call plot(0.0,   h, iup)
          call plot(0.0, 0.0, idn)
        endif
      endif
c  Finish plotting contour box
      if (box .gt. 0) then
        call plot(0., h, iup)
        call plot(w, h, idn)
        call plot(w, 0., idn)
      endif
c
c  Plot justified caption if there is one
      if (ntxt .gt. 0 ) then
        dy=dy - 3.0*max(hlet(3), hlet(4))
        call captio(dy, w, hlet(4), ntxt, text)
        if (ntxt.gt.1) ntxt=0
      endif
c  If there is a ylabel plot it
      if (ylab(1:1) .ne. ' ') then
        call justy(s, hlet(3), ylab)
        call letter(-backy-hlet(3), h/2-s(2), hlet(3), 90.0, ylab)
      endif
c  If there is a title plot it
      if (title(1:1) .ne. ' ' ) then
           call justy(s, hlet(1), title)
           call letter(w/2-s(2), h+hlet(1), hlet(1), 0.0, title)
      endif
      call justy(s, hlet(4), text(1))
c
c  Draw x-y line data: pattern lies "under" the contours
c
      do 1200 k=1, nxyp
        call newpn(ncolor(k))
        nk=nstart(k)
        ne=nstart(k+1) - 1
        do 1100 i=nk, ne
          if (xp(i) .lt. 900.0) xp(i)=w*xp(i)
          if (yp(i) .lt. 900.0) yp(i)=h*yp(i)
 1100   continue
c  Ordinary lines when nkind .le. 1000
        i1=nk
        call remark('x-y line data segment begins')
        do 1175 i=nk, ne+1
c  Check for pen lift or end of point series
          if (xp(i) .gt. 900.0 .or. i .eq. ne+1) then
            if (nkind(k).le. 1000 .and. i.gt.i1) then
              call line(nkind(k), i-i1, xp(i1), yp(i1))
            else
c  Special symbols kind=(nkind-1)/1000, height=mod(nkind,1000)
c  in .01 inches
              ksymb=(nkind(k)-1000)/1000
              glht=0.01*mod(nkind(k), 1000)
              do 1150 j=i1, i-1
                call glyph(xp(j), yp(j), glht, ksymb)
 1150         continue
            endif
            i1=i+1
          endif
 1175   continue
 1200 continue
c
c  Undo scaling in case there is another contour map
      do 1300 i=1, nstart(nxyp+1)-1
        if (xp(i) .ne. 10.0) xp(i)=xp(i)/w
        yp(i)=yp(i)/h
 1300 continue
      newp=0
c
c  Invoke the contouring routine
c
      if (m.gt.1 .and. n.gt.1) then
        call remark('Preliminaries over - contours begin')
        write(*,'(/a,i5,a)')
     $  ' Contouring now underway.  There are ',nlines,' levels.'
        call flevel
        call remark('Contouring over')
      endif
c
c  Plot any notes
c
      if (nnotes.gt.0) call pnotes
c
c  Restore origin and other things to proper state
c
      xyold(1)=w + 1.0
      xyold(2)=-oy
      call plot(xyold(1), xyold(2), -3)
      if (nauto .eq. 1) nlines=0
      ms=0
      return
      end
c______________________________________________________________
      subroutine assess
c$$$$ calls cover
c  Finds nice contour levels for the array  f. Returns
c  nlines  and contour levels in  value.
c  Avoids clustering of contours around strong peaks.
c
c
      common /evenly/ neven,zz1,zz2,dzz
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      dimension divs(7)
c
      data mlines/12/
      data divs /.05,.1,.2,.25,.5,1.0,2.0/, ten/10.0/
c
c  No array - return dummy
ce
      nlines=0
      if (n .eq. 0) return
c
c  If evenly spaced lines requested, look into that
      if (neven.eq. 3 .and.  zz2-zz1.lt.(maxcnt-1)*dzz) then
        do 1010 nlines=1, maxcnt
          zlev=zz1 + (nlines-1)*dzz
          value(nlines)=zlev
          if (abs(zlev) .lt. 1e-6*(zz2-zz1)) value(nlines)=0.0
          if (zlev .gt. zz2-dzz/2) return
 1010   continue
      endif
c
c  Generate default contour levels
      fhi=1e25
      flo=-fhi
      if (qnllhi) fhi=anulhi
      if (qnlllo) flo=anullo
      call cover(m,n,f, flo,fhi,fbar,fsig,tot)
c  Trim extreme values
      do 1101 ksp=1, 3
        flo=max(flo, fbar - 3.0*fsig)
        fhi=min(fhi, fbar + 3.0*fsig)
        call cover(m,n,f, flo,fhi,fbar,fsig,tot1)
        if (tot1 .lt. tot*0.99) goto 1150
 1101 continue
c  Make a nice covering interval
 1150 if (flo .eq. fhi) fhi=1.0 + flo + abs(0.1*flo)
      plus=1000.0 + log10(fhi - flo)
      da=10.0**mod(plus, 1.0)
      do 2100 i=1, 7
        ii=i
        if (mlines .ge. da/divs(ii)) goto 2110
 2100 continue
 2110 units=divs(ii)*ten**(int(plus)-1000)
      if (neven.eq. 1 .and.  fhi-flo.lt.(maxcnt-1)*dzz) units=dzz
      uns=aint(flo/units - 0.01)
c  Generates evenly space levels
      nlines=0
      do 2200 i=1, maxcnt
         v=units*(uns + i - 2)
         if (v.ge.flo) then
           if (v.gt.fhi .or. nlines.gt.maxcnt) goto 2300
           nlines=nlines + 1
           value(nlines)=v
         endif
 2200 continue
c  Print them out
 2300 write(*,'(a,i4,a)') 
     $' The program has selected',nlines,' contour levels'
      return
      end
c______________________________________________________________
      subroutine cover(m,n,f, flo,fhi,fbar,fsig,tot)
c$$$$ calls nothing
c  Finds extremes, mean and standard dev for f restricted by
c  flo < f < fhi.  Overwrites flo, fhi
      parameter (maxpts=350)
      dimension   f(maxpts, maxpts)
c
      fmin=f(1,1)
      fmax=fmin
      fsum=0.0
      fssq=0.0
      tot=0.001
      do 1200 j=1, n
        do 1100 i=1, m
          fij=f(i, j)
          if (fij.ge.flo.and.fij.le.fhi)  then
            fmin=min(fij, fmin)
            fmax=max(fij, fmax)
            fsum=fij + fsum
            fssq=fij**2 + fssq
            tot=tot + 1.0
          endif
 1100   continue
 1200 continue
      fbar=fsum/tot
      fsig=sqrt(fssq/tot - fbar**2)
      fhi=fmax
      flo=fmin
      return
      end
c_______________________________________________________________________
      subroutine todisk(nxy, xx, yy, zz)
c$$$$ calls nothing
c  Saves level lines to a diskfile
c  Writes coordinates xx, yy and contour height zz on each line
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $ vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      common /disk/ isave,dxy(2)
      dimension xx(*),yy(*)
c
      if (isave.le. 0) return
c
      axes=(ax(1)-ax(2))*(ax(3)-ax(4))
c  Save contours with origin at bottom left and x,y in inches
      if (axes.eq. 0.0) then
        write(12, '(/(3g12.4))') (xx(j),yy(j),zz, j=1, nxy)
c  Save contours referred to axes provided in axes command
      else
        xscal=(ax(2)-ax(1))/dxy(1)
        yscal=(ax(4)-ax(3))/dxy(2)
        write(12, '(/(3g12.4))')
     $  (ax(1)+xscal*xx(j), ax(3)+yscal*yy(j),zz, j=1,nxy)
      endif
      return
      end
c_______________________________________________________________________
c______________________________________________________________
      subroutine init(idev, plotfl)
c$$$$ calls plopen plot
c******** this routine is UNIX-specific *****************************
c  Initializes the plotfile with name  plotfl   or if this is
c  blank,  terminates plotfile in an orderly manner.
c  This routine must be replaced by appropriate routines on another
c  system.
c**********************************************************************
      character*(*) plotfl
      save isnt
c
      data isnt/0/
c
c  Terminate plotfile in unix standard plotfile format
ce
      if (plotfl .eq. ' ') then
         if (isnt .ne. 0) call plopen(-1, plotfl)
         return
       endif
c  Initiate plotfile in unix standard plotfile format
      call plopen(1, plotfl)
      call plot(0.2, 0.2, -3)
      isnt=1
      return
      end
c______________________________________________________________
      subroutine captio(y0, xlen, ht, ntxt, text)
c$$$$ calls justy letter parse
c  A series of  ntxt  80-character lines  held in text  are plotted
c  justified on a page of width  xlen.  Character height is  ht.
c  The first line of text starts at (0, y0)
      character*80 text(ntxt), word, words*128
      dimension s(4),w(80),lw(81)
c
ce
      delta=0.8*ht
      width=-delta
      y=y0
      k=0
      lw(1)=1
c
c   Examine each of the lines of text in turn
      do 1300 in=1, ntxt
        l1=1
c  Parse strips out individual words from the text. These are stuffed
c  into the long character string  words  with pointers  lw.
 1100   k=k + 1
        call parse(l1, 80, text(in), lword, word)
        lw(k+1)=lw(k) + lword
        words(lw(k):lw(k+1)-1)=word(1:lword)
c When l1>0 a new word has been successfully read from text
        if (l1 .gt. 0) then
          call justy(s, ht, word(1:lword))
          w(k)=s(3)
c Current line is shorter than the page width
          if (width+delta+w(k) .lt. xlen) then
            width=width + delta + w(k)
c Line too long for width -  plot it, justified
          else
            k1=max(1, k-1)
            space=delta + (xlen - width)/max(1, k1-1)
            x=0.0
            do 1200 ik=1, k1
              call letter(x, y, ht, 0.0, words(lw(ik):lw(ik+1)-1))
              x=x + space + w(ik)
 1200       continue
            y=y - 2.0*ht
            lw(2)=lword + 1
            words(1:lword)=word(1:lword)
            w(1)=w(k)
            k=k - k1
            width=k*w(1)
          endif
c  Go back for more text
          goto 1100
        endif
c  Current input text buffer finished - go round for the next one
        k=k - 1
 1300 continue
c
c  Finish the last line unjustified
      x=0.0
      do 1500 ik=1, k
        call letter(x, y, ht, 0.0, words(lw(ik):lw(ik+1)-1))
        x=x + delta + w(ik)
 1500 continue
      return
      end
c______________________________________________________________________
      subroutine fancy(mode)
c$$$$ calls nothing
c  Routine sets code values so that levels below vcrit are dashed, above
c  vcrit are regular, and at vcrit are drawn heavy.  A level vcrit is
c  added if necessary.  State of fancy saved in invoke
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
ce
c  If mode=0 take no action
      if (mode .eq. 0) return
c  If mode=-1 turn off the fancy option and return all kodes to 0
      if (mode .eq. -1) then
        invoke=0
        do 1400 j=1, nlines
           kode(j)=0
 1400   continue
c
c  If mode is nonnegative apply the fancy option
      else
        invoke=1
        if (nlines .eq. 0) return
        dcrit=abs(value(nlines) - value(1))/100.0
        ncrit=0
        do 1300 j=1, nlines
            if (abs(value(j)-vcrit) .le. dcrit) then
               kode(j)=2
               ncrit=1
             else
               if (value(j) .gt. vcrit) kode(j)=1
               if (value(j) .lt. vcrit) kode(j)=6
             endif
 1300   continue
c  Add vcrit to the contour list if it's not there
        if (ncrit .eq. 0) then
           nlines=nlines + 1
           kode(nlines)=2
           value(nlines)=vcrit
         endif
c
      endif
      return
      end
c______________________________________________________________
      subroutine axes (kindax, ht, xlen, xx, ixy)
c$$$$ calls linax logax
c  Decides whether to draw linear or log axes and does it
      dimension xx(*)
c
ce
      icode=10*min(3, max(0, kindax)) + ixy
      if (icode .le.  2) call linax(ht, xlen, xx, ixy)
      if (icode .eq. 11) call logax(ht, xlen, xx, ixy)
      if (icode .eq. 12) call linax(ht, xlen, xx, ixy)
      if (icode .eq. 21) call linax(ht, xlen, xx, ixy)
      if (icode .eq. 22) call logax(ht, xlen, xx, ixy)
      if (icode .ge. 31) call logax(ht, xlen, xx, ixy)
c
      return
      end
c_____________________________________________________________________
      subroutine linax(ht, xlen, xx, ixy)
c$$$$ calls justy letter nicer plot
c  Plots either an x- or a y-axis of length  xlen  inches with tick
c  marks at reasonable places between the limits  xx(1), xx(2).
c  If  xx(1) .gt. xx(2) the numerals decrease along the axis.
c  lettering height is  ht.
c  ixy =1 means  x axis,  ixy=2  means y axis.
c  The axis is drawn starting from the plot origin (x=0, y=0), which is
c  normally at the lower left corner.  This position may be moved by
c  the user in the calling routine.
c
      character*40 label,ifmt*20
      dimension xx(2),xa(2),s(4)
      common /back/ backy,box
c
      data ticfac/1.6/, iup,idn /3,2/
c
c  In-line function
      iten(x)=int(500.001 + log10(x)) - 500
c
ce
      if (xlen .eq. 0.0) return
      revers=sign(1.0, xx(2) - xx(1))
      xx(1)=revers*xx(1)
      xx(2)=revers*xx(2)
      htix=0.4*ht
      xa(1)=xx(1)
      xa(2)=ticfac*min(1.0,10.0/xlen)*(xx(2)-xx(1)) + xx(1)
      call nicer(xa, xa, xt)
      if (xx(2)-xx(1).eq.180 .or. xx(2)-xx(1).eq.360) xt=30
c  Set up format for numerical annotation on the axis
      nsiz=iten(max(abs(xx(2)), abs(xx(1))))
      ntix=iten(xt)
      nfld=3 + max(nsiz, -ntix, nsiz-ntix)
      if (nfld.gt.7) goto 1000
c  Width of field less than 8 characters - use an f-format
      ndp=max(0, -ntix)
      write(ifmt, 100) '(f', nfld, '.', ndp, ')'
 100  format(a, i1, a, i1, a)
      if (ndp .eq. 0) nfld=nfld - 1
      goto 1100
c  Width of field more than 7 characters - use a g format
 1000 ndp=nsiz - ntix + 1
      nfld=7 + ndp
      write(ifmt, 110) nfld, ndp
 110  format('(1pg', i2, '.', i2, ')' )
c
 1100 call plot(0.0, 0.0, iup)
      xs= xlen/(xx(2) - xx(1))
      n1=xa(1)/xt
      n2=xx(2)/xt + 1.0
      goto (2100, 3100), ixy
c
c  Draw the x axis
 2100 kskip=1.3 + ht*nfld/(xs*xt)
      do 2500 n=n1, n2
        x=n*xt
        xinch=xs*(x - xx(1))
c  Plot next section of axis and the tick on the right
        if (xinch.lt.-0.01 .or. xinch.gt.xlen+0.01) goto 2500
        call plot(xinch, 0.0,  idn)
        call plot(xinch, htix, idn)
c  Write numerical annotation
        if (mod(n, kskip) .ne. 0) goto 2400
        write(label, ifmt) x*revers
        call justy(s, ht, label(1:nfld))
        call letter(xinch-s(2), -2.0*ht, ht, 0.0, label(1:nfld))
c  Move back onto axis with pen up
 2400   call plot(xinch, 0.0, iup)
 2500 continue
c  Plot the last little piece of the axis
      call plot(xlen,  0.0, idn)
      xx(1)=revers*xx(1)
      xx(2)=revers*xx(2)
      return
c
c  Draw the y axis
 3100 half=-0.5*ht
      backy=half
      kskip=1.5 + ht/(xs*xt)
      do 3500 n=n1, n2
        y=n*xt
        yinch=xs*(y - xx(1))
c  Plot next section of axis and the tick on the right
        if (yinch.lt.-0.01 .or. yinch.gt.xlen+0.01) goto 3500
        call plot(0.0,  yinch, idn)
        call plot(htix, yinch, idn)
c  Write numerical annotation
        if (mod(n, kskip) .ne. 0) goto 3400
        write(label, ifmt) y*revers
        call justy(s, ht, label(1:nfld))
        backy=max(backy, s(3)-s(1)-half)
        call letter(half-s(3), yinch+half, ht, 0.0, label(1:nfld))
 3400   call plot(0.0, yinch, iup)
 3500 continue
c  Plot last little piece of axis
      call plot(0.0,  xlen, idn)
      xx(1)=revers*xx(1)
      xx(2)=revers*xx(2)
      return
      end
c______________________________________________________________
      subroutine logax(hgt, xlen, xx, ixy)
c$$$$ calls justy letter linax plot
c  Plots either an x- or a y-axis of length  xlen  inches with log
c  spaced tick marks and annotation between the limits  xx(1), xx(2).
c  ixy=1  means draw an x axis,  ixy=2  means draw a y axis.
c  The axis is drawn starting from the plot origin (x=0, y=0), which is
c  Normally at the lower left corner.  This position may be moved by
c  the user in the calling routine.
c
      character*4 label, no*2
      dimension xx(2),s(4), no(10)
      common /back/ backy,box
c
      data iup/3/,idn/2/
      data no/'1 ','2 ','3 ','4 ','5 ','6 ','7 ','8 ','9 ','10'/
c
c  In-line funcyion
      iten(x)=int(500.01 + log10(x)) - 500
c
ce
      if (xlen .eq. 0.0) return
      if (xx(1) .ge. xx(2)  .or. xx(1) .le. 0.0) goto 4000
      log1=iten(xx(1))
      log2=iten(xx(2))
      xs=xlen/log10(xx(2)/xx(1))
      if (log2.eq.log1) goto 4000
      hr=0.75*hgt
      htix=0.4*hgt
      call justy(s, hgt, no(10))
      wid10=s(3)
      call plot(0.0, 0.0, iup)
c  Select x- or y- axis
      goto (2100, 3100), ixy
c
c  Plot x axis
 2100 aback=0.5*wid10
      kskip=1.1 + 6.0*hgt/xs
      do 2500 kdec=log1, log2
        do 2400 numb=1, 9
          xinch=xs*(log10(numb/xx(1)) + kdec)
          if (xinch.lt.-0.01 .or. xinch.gt.0.01+xlen) goto 2400
c  Draw next section of axis and its tick mark
          call plot(xinch, 0.0, idn)
          if (numb.eq.1 .or. xs.gt.6.7*hgt) then
            call plot(xinch, htix, idn)
          endif
          if (numb .ne. 1) goto 2200
          if (mod(kdec-log1, kskip) .ne. 0) goto 2450
c  Write out the notation 10**kdec below the appropriate tick
          write(label, '(i4)') kdec
          call justy(s, hr, label)
          call letter(xinch+aback-s(1), -1.5*hgt, hr, 0.0, label)
          call letter(xinch-aback, -2.0*hgt, hgt, 0.0, no(10))
          goto 2450
c  If number of decades .lt.2 write the integers next to their ticks
 2200     if (log2-log1 .le. 2  .and. .02*xs .gt. hgt)
     $    call letter(xinch-0.5*hr, -2.0*hr, hr, 0.0, no(numb))
c  Move back onto the axis with the pen raised
 2450     call plot(xinch, 0.0, iup)
 2400   continue
 2500 continue
c  Plot last piece of the axis
      call plot(xlen,  0.0, idn)
      return
c
c  Plot y axis
 3100 backy=2.0*wid10
      kskip=1.1 + 2.5*hgt/xs
      do 3500 kdec=log1, log2
        do 3400 numb=1, 9
          xinch=xs*(log10(numb/xx(1)) + kdec)
          if (xinch.lt.-0.01 .or. xinch.gt.0.01+xlen) goto 3400
c  Draw next section of axis and its tick mark
          call plot(0.0, xinch,  idn)
          if (numb.eq.1 .or. xs.gt.10.*hgt) then
              call plot(htix, xinch, idn)
          endif
          if (numb .ne. 1) goto 3200
          if (mod(kdec-log1, kskip) .ne. 0) goto 3450
c  Write out the notation 10**kdec below the appropriate tick
          write(label, '(i4)') kdec
          call justy(s, hr, label)
          aback=1.5*wid10 + s(3) - s(1)
          backy=max(aback, backy)
          call letter(wid10-aback-s(1), xinch, hr, 0.0, label)
          call letter(-aback, xinch-0.5*hgt, hgt, 0.0, no(10))
          goto 3450
c  If number of decades .lt.2 write the integers next to their ticks
 3200     if (log2-log1 .le. 2  .and. .02*xs .gt. hgt)
     $    call letter(-1.5*hr, xinch-0.5*hr, hr, 0.0, no(numb))
c  Move back onto the axis with the pen raised
 3450     call plot(0.0, xinch, iup)
 3400   continue
 3500 continue
c  Plot last piece of the axis
      call plot(0.0,  xlen, idn)
      return
c
c  Interval too short to contain an exact power of ten, or
c  nonpositive limit.  Use linear axis instead of log.
 4000 call linax(hgt, xlen, xx, ixy)
      return
      end
c______________________________________________________________
      subroutine nicer(xin, xout, xtick)
c$$$$ calls nothing
c  Routine for scaling intervals and providing tick marks for axes.
c  between 7 and 15 ticks are made, which is suitable for 10in axes.
c    Input parameters
c  xin(1),xin(2)  extremes of variable  x  in its own units
c    output parameters
c  xout(1),xout(2)  adjusted extremes, made to be round numbers.  Note
c    the new interval always covers old one.
c  xtick  distance between tick marks in  x  units (not inches).  This
c    number is always a round number.
      dimension  divs(4),xin(2),xout(2)
c
      data e/1.0e-7/,divs/.1,.2,.5,1.0/
c
ce
      xout(1)=xin(1)
      xout(2)=xin(2)
      if (xout(2).eq.xout(1)) xout(2)=1.0 + 1.1*xout(2)
      plus=1000.0 + log10(xout(2)-xout(1))
      index=1.4969 + 2.857*mod(plus, 1.0)
      units=divs(index)*10.0**(int(plus)-1000)
      bias=(1.+e)*units*aint(1.+max(abs(xout(1)), abs(xout(2)))/units)
      xout(1)=xout(1) - mod(xout(1)+bias, units)
      xout(2)=xout(2) - mod(xout(2)-bias, units)
      if (abs(xout(1)/units) .le. .01) xout(1)=0.0
      if (abs(xout(2)/units) .le. .01) xout(2)=0.0
      xtick=units
      return
      end
c______________________________________________________________________
      subroutine scribe(nchar, line, xyh, leng, note)
c$$$$ calls decode
c  The string in  line  is expected to be of the form '(x,y) text'  or
c  '(x,y in) text' or '(x,y,u,v in)'.  Extracts the numbers
c  and puts text into  note.  'i' means inches, flagged by -ve xyh(5).
      character*80 note, line*(*)
      real xyh(*)
c
      kpi(i)=i + 1000 -sign(1000, i-1)
ce
      leng=-1
c  Seeks ')' in image.  Also records presence of 'i' meaning inches
c  by reversing sign of  xyh(5)
      i1=index(line, '(')
      i2=index(line, ')')
      if (i1.eq.0 .or. i1.gt.i2) return
      ii=index(line(i1:i2), 'i')
      ir=index(line(i1:i2), 'r')
      ic=index(line(i1:i2), 'c')
      iii=min (i2, kpi(ii), kpi(ir), kpi(ic))
c  For left, center, right justification xyh(8)=0, 0.5, 1
      xyh(8)=0.5*min(2, min(ic,1)+min(ir,2))
c
c  Translate up to 7 numbers found in parentheses in  line
      call decode(line(i1+1:iii-1), xyh, 7, nfound)
c  2 numbers: x,y         - no arrow, same height, angle 0, old color
c  3 numbers: x,y, height - no arrow, angle 0, old color
c  4 numbers: p,q,x,y     - same height, angle zero
      if (nfound.le.1) return
      if (nfound.eq.3) xyh(5)=xyh(3)
      if (nfound.eq.2 .or. nfound.eq. 3) then
        xyh(3)=xyh(1)
        xyh(4)=xyh(2)
      endif
      if (ii .ne. 0) xyh(5)=-abs(xyh(5))
      if (nfound.lt.6) xyh(6)=0.0
c  Copy remaining portion of line into the character variable  note
      note=line(i2+1:nchar)
      leng=min(80, nchar - i2)
      return
      end
c______________________________________________________________
      subroutine pnotes
c$$$$ calls justy letter newpn plot
c
c     Plot the notes file and the arrows
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      character*80 note(20), text(10)
      character*116 format,file,title,output,xlab,ylab
      common /charco/  format,file,title,output,xlab,ylab,text,note
c
      real ss(4),p(4)
      character*80 tnote
c
      data iup/3/,idn/2/
      data cosa,tana,tip/0.866,0.577,0.6/, zone/1.2/
c
ce
      h=height
      w=width
      if (h.le.0 .and. w.le.0) then
         if (m.ge.n) w=6
         if (m.lt.n) h=6
      endif
      if (h.le.0) h=w*(n-1)/(m-1.0)
      if (w.le.0) w=h*(m-1)/(n-1.0)
c
c  Run through the notes in turn
      do 3550 inote=1, nnotes
        tnote=note(inote)
        long=leng(inote)
c
c  Set drawing color
        call newpn(nint(xyh(7, inote)))
c
        do 3100 i=1, 4
          p(i)=xyh(i,inote)
 3100   continue
        hn=abs(xyh(5,inote))
        if (long.eq.1 .and. tnote(1:1).eq.' ') hn=0.0
        if (xyh(5,inote) .ge.0.0) then
c  Rescale coordinates accounting for log or linear axes
c  Suppress if outside the frame
        do 3200 i=1, 3, 2
          if ((p(i  )-ax(1))*(p(i  )-ax(2)) .gt.0.0) goto 3550
          if ((p(i+1)-ax(3))*(p(i+1)-ax(4)) .gt.0.0) goto 3550
          if (kindax.eq.0 .or. kindax.eq.2) then
            p(i  )= (p(i  ) - ax(1))*w/(ax(2) - ax(1))
          else
            p(i  )=log(p(i  )/ax(1))*w/log(ax(2)/ax(1))
          endif
          if (kindax.eq.0 .or. kindax.eq.1) then
            p(i+1)= (p(i+1) - ax(3))*h/(ax(4) - ax(3))
          else
            p(i+1)=log(p(i+1)/ax(3))*h/log(ax(4)/ax(3))
          endif
 3200   continue
        endif
        shaft=abs(p(1)-p(3))+abs(p(2)-p(4))
c  Recenter text according to justification in xyh(8,:)
        call justy(ss, hn, tnote(1:long))
        cosr=cos(0.01745329*xyh(6,inote))
        sinr=sin(0.01745329*xyh(6,inote))
        p(3)=p(3) - xyh(8,inote)*ss(3)*cosr
        p(4)=p(4) - xyh(8,inote)*ss(3)*sinr
c  There is an arrow if the stalk is of positive length
        if (shaft .gt.0) then
c  Plot the arrow if there is one
          p5=ss(2)  - (p(1)-p(3))*cosr - (p(2)-p(4))*sinr
          p6=0.3*hn + (p(1)-p(3))*sinr - (p(2)-p(4))*cosr
          el=ss(3) - ss(2) + 0.5*hn
          te=p5 - p5*(el - hn)/(1.0e-6 + el + abs(p5))
          p8=p6 - sign(zone*hn, p6)
          p7=te*p8/(p6 + sign(1.0e-6, p6))
          if (abs(p7 - p5) .gt. el) then
            p7=p5 - sign(el, p5)
            p8=p6*p7/te
          endif
          r=sqrt(p7**2 + p8**2)
          c=tip*cosa*abs(xyh(5,inote))/r
          s=tana*c
          q7=p7*cosr - p8*sinr
          q8=p7*sinr + p8*cosr
          call plot(p(1) + c*q7 + s*q8, p(2) - s*q7 + c*q8, iup)
          call plot(p(1),               p(2),               idn)
          call plot(p(1) + c*q7 - s*q8, p(2) + s*q7 + c*q8, idn)
          call plot(p(1),               p(2),               iup)
          call plot(p(1)+q7,            p(2)+q8,            idn)
        endif
c  Write the text of the note
        call letter(p(3), p(4), hn, xyh(6,inote), tnote(1:long))
 3550 continue
c
c     Reset drawing color
      call newpn(ikolor)
c
      return
      end
c______________________________________________________________________
c---------------------------------------------------------------------
c                       This is Unit C 
c  Unit A contains routines for language interpretation and data input.
c  Unit B contains page arrangement, axes and other decoration.
c  Unit C contains contour generation and labeling routines.
c---------------------------------------------------------------------
      subroutine flevel
c$$$$ calls circ dotted draw1 draw2 dumpit flag halt remark 
c$$$$ sflag, uflag
c
c
c     flevel--  draws contour lines on gridded array
c
c  Principal routine for contouring rectangular arrays of real data.
c  Algorithm by  Dr David P. Anderson
c    modified by R. L. Parker (1984, 1989, 1990)
c    modified by Peter Salameh (1989)
c    modified by David W. Caress (1989)
c
c
c
c  Conventions
c  The standard plotting convention is that  x  increases with i,  the
c  row index of  f(i,j)  and  y  increases with   j.  The point (0, 0)
c  with respect to current plot origin corresponds to the matrix
c  element  f(1,1).
c
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /flags/ nflag, kntflg(maxflg)
c
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      common /little/ eps,epps
c
      common /affine/ acoef, bcoef
c
      real bx(3),by(3),p1(2),p2(2),c1(2),c2(2)
      real ang1(2),ang2(2),ang3(2),ang4(2),ang5(2),ang6(2)
      real small,big,range,val,a,b,axx,ayy
      real dist1,dist3,xs,ys
c
      integer ixof(3,2,2),iyof(3,2,2),koff(3,2,2),doff(3,2,2)
c
      logical ifcase(3),closed
      logical flag,qq
      character*30 messag
c
      c(a, b)=(a - val)/(a - b)
c
      data ixof/0,0,1,0,1,0, 0, 0, 1,-1,-1,-1/
      data iyof/0,1,0,1,0,0,-1,-1,-1, 0, 0, 1/
      data koff/2,1,2,1,2,1, 2, 1, 2, 1, 2, 1/
      data doff/2,1,1,1,1,2, 2, 2, 1, 2, 2, 1/
c
c  Bound used later in dumpit
c  epps used in circ
c
ce
      bound(1)=min(0.0, (nx-1)*delx)
      bound(2)=max(0.0, (nx-1)*delx)
      bound(3)=min(0.0, (ny-1)*dely)
      bound(4)=max(0.0, (ny-1)*dely)
      epps=1e-5*max(abs((nx-1)*delx), abs((ny-1)*dely))**2
c
c  Find largest, smallest function values to aid user
c
      big=f(1,1)
      small=big
      qq=.FALSE.
      do 1000 i=1, nx
        do 1001 j=1, ny
          if ((qnllhi.and.f(i,j).ge.anulhi).or.
     $      (qnlllo.and.f(i,j).le.anullo)) goto 1001
          if (qq) then
            big=max(big, f(i,j))
            small=min(small, f(i,j))
          else
            big=f(i,j)
            small=f(i,j)
            qq=.TRUE.
          endif
 1001   continue
 1000 continue
      if (.not.qq) then
        write(*,*)'>>> Entire array nulled out - no contours plotted'
        return
      endif
      range=big - small
c
c
c  Make sure that no grid points have exact contour values
      eps=range*1e-6
      tol=1e-7
      do 1004 i=1,nx
        do 1003 j=1,ny
          do 1002 k=1,nlines
            if (abs(f(i,j)-value(k)).le.abs(value(k))*tol) then
              if (value(k).ne. 0.0) f(i,j)=value(k)*(1.0+tol)
              if (value(k).eq. 0.0) f(i,j)=range*tol
            endif
 1002     continue
 1003   continue
 1004 continue
c
c  Loop for contour values
c  (check for values outside of function range)
c
      if (acoef .ne. 1.0 .or. bcoef .ne. 0.0)
     $write(*,'(2a,2g13.5/a)') 
     $' An affine mapping applied to input array, ',
     $'parameters:', acoef,bcoef,' After the mapping -'
      write(*,'(a,g13.5,a,g13.5)')
     $' Grid min value:',small,'  max value: ',big
c
      do 1008 iline=1, nlines
        val=value(iline)
        if (val.lt.small.or.val.gt.big) goto 1008
        write(messag,'(a,g12.4)')'Begin level ',val
        call remark(messag)
c
c  Find the points where contour intersects grid lines
c  (reset number of segment found to zero)
c
        nflag=0
        do 1009 i=1, nx
          do 1010 j=1, ny
            a=f(i,j)
            if ((qnllhi.and.a.ge.anulhi).or.
     $        (qnlllo.and.a.le.anullo)) goto 1010
c
c  Look for horizontal element k=1
c
            if (i.lt.nx) then
              b=f(i + 1, j)
              if ((.not.qnllhi.or.b.lt.anulhi).and.
     $            (.not.qnlllo.or.b.gt.anullo)) then
                if (a.le.val.and.b.gt.val.or.b.lt.val.and.a.ge.val)
     $            call sflag(1, i, j)
              endif
            endif
c
c  Look for vertical element k=2
c
            if (j.lt.ny) then
              b=f(i, j + 1)
              if ((.not.qnllhi.or.b.lt.anulhi).and.
     $           (.not.qnlllo.or.b.gt.anullo)) then
                 if (a.le.val.and.b.gt.val.or.b.lt.val.and.a.ge.val)
     $             call sflag(2, i, j)
              endif
            endif
 1010     continue
 1009   continue
c
c  All points found -- now look for components
c
 1012   continue
c
c  The variables in this part are
c  ix, iy           indices of start of current edge
c  ik               direction of edge, 1=x, 2=y
c  idir              direction out of edge, 1=up/right, 2=left/down
c  closed           true if component is a loop
c  ix0,iy0          if closed component, starting point
c
c  Search for a point on border to start from
c
c ix=1 edge (left)
c
        closed=.false.
          do 1013 iy=1, ny-1
            if(flag(2,1,iy)) then
              ix=1
              ik=2
              idir=1
              goto 100
            endif
 1013     continue
c
c iy=1 edge top
c
          do 1014 ix=1, nx-1
            if(flag(1,ix,1)) then
              iy=1
              ik=1
              idir=1
              goto 100
            endif
 1014     continue
c
c iy=ny edge (bottom)
c
          do 1015 ix=1, nx-1
            if (flag(1,ix,ny)) then
              iy=ny
              ik=1
              idir=2
              goto 100
            endif
 1015     continue
c
c ix=nx edge (right)
c
          do 1016 iy=1, ny-1
            if(flag(2,nx,iy)) then
              ix=nx
              ik=2
              idir=2
              goto 100
            endif
 1016     continue
c
c  If none, search for point in middle
c
        do 1017 i=1, nflag
          if (kntflg(i).gt.0) then
            kflag=kntflg(i)
            ix=kflag/100000
            kflag=kflag - 100000*ix
            iy=kflag/10
            ik=kflag - 10*iy
            if (ix.lt.1.or.ix.gt.nx) goto 1017
            if (iy.lt.1.or.iy.gt.ny) goto 1017
            idir=1
            closed=.true.
            ix0=ix
            iy0=iy
            ik0=ik
            goto 100
          endif
 1017   continue
c
c  If can't find a starting point, we're done with this contour value
c
        goto 1018
c
c  Here when found starting point of a component
c  Compute its xy coordinates
c
  100   continue
        axx=(ix-1)*delx
        ayy=(iy-1)*dely
        if(ik.eq.1) then
c
          if (interp.eq.0) then
c           do linear interpolation
            axx=axx + c(f(ix,iy), f(ix+1,iy))*delx
c
          else
c         Do parabolic interpolation
c
            if (ix.eq.1) then
c           Single parabola approximation
              axx=axx + (1.0
     $            - parab(f(ix+2,iy),f(ix+1,iy),f(ix,iy),val))*delx
            else if (ix.ge.(nx-1)) then
c           Single parabola approximation
              axx=axx
     $            + parab(f(ix-1,iy),f(ix,iy),f(ix+1,iy),val)*delx
            else
c           Double parabola approximation
              a1=axx
     $            + parab(f(ix-1,iy),f(ix,iy),f(ix+1,iy),val)*delx
              a2=axx + (1.0
     $            - parab(f(ix+2,iy),f(ix+1,iy),f(ix,iy),val))*delx
              axx=(a1 + a2)/2
            endif
          endif
c
        else
c
          if (interp.eq.0) then
c         Do linear interpolation
            ayy=ayy + c(f(ix,iy), f(ix,iy+1))*dely
c
          else
c         Do parabolic interpolation
c
            if (iy.eq.1) then
c           Single parabola approximation
              ayy=ayy + (1.0
     $            - parab(f(ix,iy+2),f(ix,iy+1),f(ix,iy),val))*dely
            else if (iy.ge.(ny-1)) then
c           Single parabola approximation
              ayy=ayy
     $            + parab(f(ix,iy-1),f(ix,iy),f(ix,iy+1),val)*dely
            else
c           Double parabola approximation
              a1=ayy
     $            + parab(f(ix,iy-1),f(ix,iy),f(ix,iy+1),val)*dely
              a2=ayy + (1.0
     $            - parab(f(ix,iy+2),f(ix,iy+1),f(ix,iy),val))*dely
              ayy=(a1 + a2)/2
            endif
          endif
        endif
c
c  Start a list (in x, y) of the points on the component
c
        xsave(1)=axx
        ysave(1)=ayy
        nsave=1
        if(.not.closed)call uflag(ik, ix, iy)
c
c  Loop for points in this component
c
 1019   continue
c
c  From the current edge, there are three possible edges to go to.
c  See which of these edges are on the grid and contain a point.
c
          do 1021 icase=1, 3
            ifcase(icase)=.false.
            jx=ix + ixof(icase,ik,idir)
            jy=iy + iyof(icase,ik,idir)
            jk=koff(icase,ik,idir)
            if(jx.gt.0.and.jy.gt.0) then
c
c  Note flag is false for nonexistent edges
c
              if(flag(jk,jx,jy)) then
                ifcase(icase)=.true.
                bx(icase)=(jx-1)*delx
                by(icase)=(jy-1)*dely
                if(jk.eq.1) then
                  if (interp.eq.0) then
c                 Do linear interpolation
                    bx(icase)=bx(icase)
     $                         + c(f(jx,jy),f(jx+1,jy))*delx
                  else
c                 Do parabolic interpolation
                    if (jx.eq.1) then
c                   Single parabola approximation
                      bx(icase)=bx(icase) + (1.0 - parab(f(jx+2,jy),
     $                              f(jx+1,jy),f(jx,jy),val))*delx
                    else if (jx.ge.(nx-1)) then
c                   Single parabola approximation
                      bx(icase)=bx(icase) + parab(f(jx-1,jy),
     $                              f(jx,jy),f(jx+1,jy),val)*delx
                    else
c                   Double parabola approximation
                      b1=bx(icase) + (1.0 - parab(f(jx+2,jy),
     $                      f(jx+1,jy),f(jx,jy),val))*delx
                      b2=bx(icase) + parab(f(jx-1,jy),
     $                       f(jx,jy),f(jx+1,jy),val)*delx
                      bx(icase)=(b1 + b2)/2
                    endif
                  endif
c
                else
                  if (interp.eq.0) then
c                 Do linear interpolation
                    by(icase)=by(icase)
     $                          + c(f(jx,jy),f(jx,jy+1))*dely
                  else
c                 Do parabolic interpolation
                    if (jy.eq.1) then
c                   Single parabola approximation
                      by(icase)=by(icase) + (1.0 - parab(f(jx,jy+2),
     $                            f(jx,jy+1),f(jx,jy),val))*dely
                    else if (jy.ge.(ny-1)) then
c                   Single parabola approximation
                      by(icase)=by(icase) + parab(f(jx,jy-1),
     $                            f(jx,jy),f(jx,jy+1),val)*dely
                    else
c                   Double parabola approximation
                      b1=by(icase) + (1.0 - parab(f(jx,jy+2),
     $                     f(jx,jy+1),f(jx,jy),val))*dely
                      b2=by(icase) + parab(f(jx,jy-1),
     $                     f(jx,jy),f(jx,jy+1),val)*dely
                      by(icase)=(b1 + b2)/2
                    endif
                  endif
                endif
c
              endif
            endif
 1021     continue
c
c  See if there wasn't a point where one was expected
c  This is possible because of exclusion values anullo and anulhi.
c  If .not. closed then end of curve.
c  If closed, then reverse order of points found, make .not. closed,
c  and continue from starting point in oposite direction.
c
          if(.not.(ifcase(1).or.ifcase(2).or.ifcase(3))) then
            if (.not.closed) then
              goto 1020
            else
              do 10 i=1, nsave/2
                xs=xsave(i)
                ys=ysave(i)
                xsave(i)=xsave(nsave-i+1)
                ysave(i)=ysave(nsave-i+1)
                xsave(nsave-i+1)=xs
                ysave(nsave-i+1)=ys
 10           continue
              call uflag(ik0, ix0, iy0)
              closed=.false.
              ix=ix0
              iy=iy0
              ik=ik0
              idir=2
              goto 1019
            endif
          endif
c
c  If an choice, go across the square
c  (looks like across is LAST choice!!)
c  (1=left side, 3=right side, 2=straight across)
c
          if(ifcase(1).and.ifcase(3)) then
            dist1=(axx-bx(1))**2 + (ayy-by(1))**2
            dist3=(axx-bx(3))**2 + (ayy-by(3))**2
            if(dist1.lt.dist3) then
              ibest=1
            else
              ibest=3
            endif
          else if(ifcase(1)) then
            ibest=1
          else if(ifcase(3)) then
            ibest=3
          else
            ibest=2
          endif
c
c  Update the coords, etc.
c
          ix=ix + ixof(ibest,ik,idir)
          iy=iy + iyof(ibest,ik,idir)
          idir=doff(ibest,ik,idir)
          ik=koff(ibest,ik,idir)
c
c  Add point to the list
c
          if(nsave.gt.maxcmp)
     $      call halt('>>> Out of room - raise parameter  MAXCMP')
c
          nsave=nsave + 1
          axx=bx(ibest)
          ayy=by(ibest)
          xsave(nsave)=axx
          ysave(nsave)=ayy
c
c  Remove the point from grid so you don't use it again
c
          call uflag(ik, ix, iy)
c
c  Check for termination of component
c
          if(closed) then
            if(ix.eq.ix0.and.iy.eq.iy0.and.ik.eq.ik0) then
              goto 1020
            endif
          else
            if(ix.eq.1.and.ik.eq.2.or.
     +      ix.eq.nx.and.ik.eq.2.or.
     +      iy.eq.1.and.ik.eq.1.or.
     +      iy.eq.ny.and.ik.eq.1) then
              goto 1020
            endif
          endif
          goto 1019
c
 1020   continue
c
c  Here when finished finding points in component.
c  Points are in arrays xsave, ysave.  "closed" tells component type.
c
c
c Dump contour line to plotter here
c
c        if (closed) then
c          write(*, 5200)nsave
c 5200     format(' Drawing closed contour: ',i4,' segments')
c        else
c          write(*, 5300)nsave
c 5300     format(' Drawing open contour:   ',i4,' segments')
c        endif
c
c       If nsave < minimum number of segments, don't plot it!
        if (nsave.lt.nsgmnt) then
          continue
c
c       Else if component has only two points just draw line segment
        else if (nsave.eq.2) then
          p1(1)=xsave(1)
          p1(2)=ysave(1)
          p2(1)=xsave(2)
          p2(2)=ysave(2)
          call dotted(iline, p1, p2)
c
c       Else if not to be smoothed, just dumpit
        else if (mooth.le.1) then
          call dumpit(kode(iline),val,kolor(iline),nsave,xsave,ysave)
c
c       Else draw nice curves
        else
c
          if (closed) then
c         Draw closed curve
c
c           First remove last point (which duplicates first point)
            nsave=nsave - 1
c
            call circ(nsave, 1, 2, c1, rad1, ang1, ang2, ang3)
            do 1031 i=1,nsave
              j=mod(i,nsave)+1
              k=mod(i+1,nsave)+1
              call circ(i, j, k, c2, rad2, ang4, ang5, ang6)
              call draw2(c1, c2, rad1, rad2, ang2, ang3, ang4, ang5)
              c1(1)=c2(1)
              c1(2)=c2(2)
              rad1=rad2
              ang2(1)=ang5(1)
              ang2(2)=ang5(2)
              ang3(1)=ang6(1)
              ang3(2)=ang6(2)
 1031       continue
c
          else
c         Draw open curve
            call circ(1, 2, 3, c1, rad1, ang1, ang2, ang3)
            call draw1(c1, rad1, ang1, ang2)
            do 1030 i=2,nsave-2
              call circ(i, i+1, i+2, c2, rad2, ang4, ang5, ang6)
              call draw2(c1, c2, rad1, rad2, ang2, ang3, ang4, ang5)
              c1(1)=c2(1)
              c1(2)=c2(2)
              rad1=rad2
              ang2(1)=ang5(1)
              ang2(2)=ang5(2)
              ang3(1)=ang6(1)
              ang3(2)=ang6(2)
 1030       continue
            call draw1(c1, rad1, ang2, ang3)
          endif
        endif
c
        goto 1012
c
 1018   continue
c       Flush last buffer in dotted
c
        p1(1)=-99999.9
        p1(2)=-99999.9
        call dotted(-1, p1, p2)
c
 1008 continue
c
c End flevel
c
      return
      end
c______________________________________________________________
      subroutine dumpit(kind, val, icol, nxy, xx, yy)
c$$$$ calls justy letter line mabel newpn sort todisk
c  Arranges for the contour defined by the vectors  xx, yy  to be
c  plotted. Major activity here is the positioning of the label if
c  any.  kind controls line weight and labeling by kind >=0 a label is
c  requested; kind < 0 no label.
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      common /plotco/ idev,ikolor,inf,ms,ns,m,n,klock,nskip,nauto,
     $         kindax,ntxt,invoke,height,width,ax(4),hlet(5),
     $                 vcrit,xyorg(2),xyold(2),nnotes,xyh(9,20),leng(20)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      common /hatch/ klosed,lftrt
      parameter (nulb=500)
      common /xylabs/ xlb(nulb),ylb(nulb),nlb
c
      character label*20
      dimension xx(nxy),yy(nxy),s(4),ls(20),tag(20),key(20),adsite(5)
      save gap,centre,oldval,label
c
      data oldval/1.0e25/, adsite/0.5,0.34,0.67,0.25,0.75/
c
c  In-line function
      slope(x1, x2, y1, y2)=abs(y2-y1)/(abs(x2-x1) + 1.0e-8)
c
c  Send level lines to a file without labels, etc
       call todisk(nxy,xx,yy, val) 
c
c  Set plotting color
ce
      call newpn(icol)
c  Open or closed contour? Sent through /hatch/ to routine line
      klosed=1.0-min(1.0, (abs(xx(1)-xx(nxy))+abs(yy(1)-yy(nxy)))*1e3)
c  Decide on direction for downhill
      i1=xx(1)/delx + 1.1
      j1=yy(1)/dely + 1.1
      lftrt=sign(1.0, ((yy(2)-yy(1))*((i1-1)*delx - xx(1)) -
     $(xx(2)-xx(1))*((j1-1)*dely - yy(1)))*(val- f(i1,j1)))
c
c  If this is a new contour build the string to label it.  Note
c  that the label is scaled by  scalev.
c
      if (oldval .ne. val) then
         call mabel(scalev*val, len, label)
         call justy(s, hlab, label(1:len))
         gap=s(4) - s(1) + hlab
         centre=2.0*s(2)
         oldval=val
      endif
c
c  No label requested.  Just dump the contour as it is
c
      if (kind .lt. 0) then
         call line(-kind, nxy, xx, yy)
         return
      endif
c
c  Make a list of potential sites for a label.
c  Identify places where slope changes sign; save y gradient in tag
c
      k=0
      do 1100 i=2, nxy-1
        if ((yy(i)-yy(i-1))*(yy(i+1)-yy(i)) .le. 0.0) then
          k=min(1 + k, 15)
          ls(k)=i
          i1=xx(i)/delx + 1.49
          j1=max(2.0, yy(i)/dely + 1.49)
          tag(k)=abs(f(i1,j1) - f(i1,j1-1))
        endif
 1100 continue
      call sort(k, tag, key)
c
c  Add 5 more sites; save local slope in tag
c
      nsites=k+5
      do 1300 i=k+1, nsites
        ls(i)=1 + min(int(nxy*adsite(i-k)), nxy-2)
        tag(i)=slope(xx(ls(i)), xx(ls(i)+1), yy(ls(i)), yy(ls(i)+1))
 1300 continue
      call sort(nsites-k, tag(k+1), key(k+1))
c
c  Scan potential sites, testing 3 criteria as identified in the code
c
      do 2500 k1=1, nsites
c
c  Run through list in order of label site quality
        i1=ls(key(k1))
c
c  Create a gap wide enough to hold the label
        do 2100 j=1, nxy
          i1=max(1, i1-1)
          i2=min(nxy, i1 + 2*j)
          dsq=(xx(i1)-xx(i2))**2+(yy(i1)-yy(i2))**2
          if (dsq .gt. gap**2) goto 2200
 2100   continue
c
c  Dropped through loop: contour too small - omit the label and plot
c  in one piece.
         call line(kind, nxy, xx, yy)
         return
c
c   h1 is deviation of contour arc from chord over the gap;
c   d1 is the distance to the nearest wall.
 2200   mid=(i1 + i2)/2
        h1=(xx(i1)*(yy(mid)-yy(i2))+xx(i2)*(yy(i1)-yy(mid))+
     $       xx(mid)*(yy(i2)-yy(i1)))/sign(sqrt(dsq), xx(i1)-xx(i2))
c
         d1=min(xx(i1)-bound(1), bound(2)-xx(i1), yy(i1)-bound(3),
     $   bound(4)-yy(i1), xx(i2)-bound(1),bound(2)-xx(i2),
     $   yy(i2)-bound(3), bound(4)-yy(i2))
c  Check for label collisions
        tolsq=(1.5*hlab)**2
        do 2300 l=1, min(nulb, nlb)
          dist=sqrt((xx(mid)-xlb(l))**2+(yy(mid)-ylb(l))**2)
          if ((xx(mid)-xlb(l))**2+(yy(mid)-ylb(l))**2
     $          .lt. tolsq) goto 2500
 2300   continue
c
c  If curvature is mild enough and label does not collide with the wall
c  the label passes the test and will be drawn in the gap.
        if (abs(h1) .lt. 3*hlab .and. d1 .gt. 0.7*hlab) goto 3000
c
 2500 continue
c
c  Dropped through loop: unable to meet requirements.  Omit label.
c
      call line(kind, nxy, xx, yy)
      return
c
c
c  Plot a contour with a label.
c  First dump the contour in two pieces.
c
 3000 call line(kind, i1, xx, yy)
      call line(kind, nxy-i2+1, xx(i2), yy(i2))
c
c  Draw the label tastefully rotated and positioned in the gap
c
      hprime=hlab - h1
      dx=xx(i1) - xx(i2)
      angl=atan((yy(i1)-yy(i2))/(dx + sign(.001, dx)))
      ox=0.5*(xx(i1) + xx(i2) - centre*cos(angl) + hprime*sin(angl))
      oy=0.5*(yy(i1) + yy(i2) - centre*sin(angl) - hprime*cos(angl))
      call letter(ox, oy, hlab, 57.29*angl, label)
c
c  Record label center coords
      nlb=nlb+1
      k=mod(nlb-1, nulb)+1
      xlb(k)=xx(mid)
      ylb(k)=yy(mid)

c     Reset plotting color
      call newpn(ikolor)
c
      return
      end
c______________________________________________________________
      subroutine sort(nx, x, key)
c$$$$ calls nothing
c     In-place rapid sorter with key
c  x()   orders x  to be increasing, overwriting original.
c  key() index with key(1)=index in unsorted array of min(x),
c        key(2)=index of next largest, ... key(nx)=index of max(x).
c
      dimension x(nx), key(nx)
c
      do 1 i=1,nx
 1    key(i)=i
      mo=nx
 2    if (mo-15) 21,21,23
 21   if (mo.le.1) return
      mo=2*(mo/4) + 1
      go to 24
 23   mo=2*(mo/8) + 1
 24   ko=nx - mo
      jo=1
 25   i=jo
 26   if (x(i) - x(i+mo)) 28,28,27
 27   temp=x(i)
      x(i)=x(i+mo)
      x(i+mo)=temp
      kemp=key(i)
      key(i)=key(i+mo)
      key(i+mo)=kemp
      i=i - mo
      if (i-1) 28,26,26
 28   jo=1 + jo
      if (jo-ko) 25,25,2
      end
c______________________________________________________________
      subroutine sflag(k, i, j)
c$$$$ calls halt
c  Set flag
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /flags/ nflag, kntflg(maxflg)
c
ce
      kflag=i*100000 + j*10 + k
      nflag=nflag + 1
      if (nflag .gt. maxflg)
     $  call halt('>>> Out of room - raise parameter MAXFLG')
      kntflg(nflag)=kflag
 999  continue
      return
      end
c______________________________________________________________
      subroutine uflag(k, i, j)
c$$$$ calls iflag
c  Unset flag
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /flags/ nflag, kntflg(maxflg)
c
ce
      nn=iflag(k, i, j)
      if (nn.gt.0) then
        kntflg(nn)=-kntflg(nn)
      else
        write(*,*)' Request to clear unset flag!!'
      endif
      return
      end
c______________________________________________________________
      subroutine dotted(line, p1,  p2)
c$$$$ calls dumpit
c  Collects a buffer of points to be plotted in the arrays  xx, yy.
c  If a discontinuity is encountered or the buffer is filled, flushes
c  the buffer with routine  dumpit.
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      logical qnllhi, qnlllo
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      real p1(2),p2(2),xx(maxcmp),yy(maxcmp)
      save xx,yy,nxy,lastln
c
      data nxy/0/, lastln/-1/
c
c  Advance the buffer count
ce
      nxy=1 + nxy
c
c  If flagged to flush but nothing to flush return
      if (nxy.eq.1.and.line.le.0) return
c
c  Initialize a new buffer
 1000 if (nxy.eq.1) then
       lastln=line
       xx(nxy)=p1(1)
       yy(nxy)=p1(2)
       nxy=2
       xx(nxy)=p2(1)
       yy(nxy)=p2(2)
       return
      endif
c
c  Add another point to the buffer.  Only done if this line piece
c  is contiguous with the previous one and if the buffer isn't full
      if ((xx(nxy-1)-p1(1))**2+(yy(nxy-1)-p1(2))**2.lt.1e-6
     $   .and. nxy.le.maxcmp)  then
        xx(nxy)=p2(1)
        yy(nxy)=p2(2)
        return
      endif
c
c  Flush the buffer.  Noncontiguous line piece has been encountered.
      nxy=nxy - 1
      call dumpit(kode(iline),value(iline),kolor(iline),nxy,xx,yy)
c  Re-initialize the buffer
      nxy=1
      if (line .gt. 0) goto 1000
c
c  If  nline is not positive reinitialize buffer
c
      nxy=0
      lastln=-1
      return
      end
c______________________________________________________________
      subroutine circ(i, j, k, center, radius, ang1, ang2, ang3)
c$$$$ calls nothing
c
c  Routine to find the center, radius and angular limits of the
c  circular arc determined by the points with indices i, j, k in the
c  arrays x and y.  If colinear return radius=-1 and put the
c  three points in ang1, ang2 and ang3.
c
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /levlco/  xsave(maxcmp),ysave(maxcmp),value(maxcnt),
     $     f(maxpts, maxpts)
c
      common /little/ eps,epps
      real center(2),ang1(2),ang2(2),ang3(2)
c
      data pi2/6.283185/
c
c  The formulas here are hard to derive, at least they were for me
c
ce
      a=(xsave(i)+xsave(j))/2.0
      b=(ysave(i)+ysave(j))/2.
      c=ysave(j)-ysave(i)
      d=xsave(i)-xsave(j)
      e=(xsave(j)+xsave(k))/2.
      ff=(ysave(j)+ysave(k))/2.
      g=ysave(k)-ysave(j)
      h=xsave(j)-xsave(k)
      disc=d*g-c*h
c
c  Check for colinearity
c
      if(abs(disc).lt.epps) then
        radius=-1.
        ang1(1)=xsave(i)
        ang1(2)=ysave(i)
        ang2(1)=xsave(j)
        ang2(2)=ysave(j)
        ang3(1)=xsave(k)
        ang3(2)=ysave(k)
      else
        t=(d*a+ff*c-b*c-d*e)/disc
        center(1)=e+g*t
        center(2)=ff+h*t
        radius=sqrt((center(1)-xsave(i))**2+(center(2)-ysave(i))**2)
        a=xsave(i)-center(1)
        b=ysave(i)-center(2)
        c=xsave(j)-center(1)
        d=ysave(j)-center(2)
        e=xsave(k)-center(1)
        ff=ysave(k)-center(2)
        ang1(1)=atan2(b,a)
        ang2(1)=atan2(d,c)
        ang3(1)=atan2(ff,e)
c
c  Make sure angles are in the right order by adding pi*2 to
c  some of them, if needed.
c
        if(ang1(1).lt.ang2(1)) then
          if(ang1(1).lt.ang3(1).and.ang3(1).lt.ang2(1)) then
            ang1(1)=ang1(1)+pi2
          else if(ang3(1).lt.ang1(1)) then
            ang3(1)=ang3(1)+pi2
          endif
        else
          if(ang2(1).lt.ang3(1).and.ang3(1).lt.ang1(1)) then
            ang1(1)=ang1(1)-pi2
          else if(ang3(1).gt.ang1(1)) then
            ang3(1)=ang3(1)-pi2
          endif
        endif
      endif
      return
      end
c______________________________________________________________
      subroutine draw1(center, radius, ang1, ang2)
c$$$$ calls dotted vmov2
c  Routine to draw a circular arc
c  radius=0 means draw segments from point ang1 to ang2
c
      logical qnllhi, qnlllo
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      real center(2),radius,ang1(2),ang2(2),p1(2),p2(2)
c
ce
      if(radius.ge.0) then
        dang=ang2(1)-ang1(1)
        arclng=abs(dang)*radius
        nsteps=arclng*mooth+1
        da=dang/nsteps
        do 1000 i=0, nsteps
          angle=ang1(1)+i*da
          p2(1)=radius*cos(angle)+center(1)
          p2(2)=radius*sin(angle)+center(2)
          if(i.ne.0) call dotted(iline, p1, p2)
          call vmov2(p2, p1)
1000    continue
      else
        call dotted(iline, ang1, ang2)
      endif
      return
      end
c______________________________________________________________
      subroutine draw2(cent1, cent2, rad1, rad2, ang1, ang2, ang3, ang4)
c$$$$ calls dotted vmov2
c  Routine to draw the average of two circular arcs
c
c  cent1 and cent2 are the centers
c  rad1 and rad2 are the radii
c  If either radius is zero it means the arc is really a line
c  segment whose endpoints are in (ang1,ang2) or (ang3,ang4)
c
      logical qnllhi, qnlllo
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /lineco/ kode(maxcnt),kolor(maxcnt),nsave,delx,dely,
     $     nx,ny,qnllhi,qnlllo,anulhi,anullo,nlines,iline,
     $     scalev,hlab,mooth,interp,nsgmnt,bound(4)
c
      real cent1(2),cent2(2),ang1(2),ang2(2),ang3(2),ang4(2)
      real p1(2),p2(2)
c
c  Find angular ranges (or vector)
c
ce
      dang11=ang2(1)-ang1(1)
      dang12=ang2(2)-ang1(2)
      dang21=ang4(1)-ang3(1)
      dang22=ang4(2)-ang3(2)
c
c  Find approximate length of arc or segment
c
      if(rad1.ge.0) then
        arclng=abs(dang11)*rad1
      else
        arclng=sqrt(dang11**2+dang12**2)
      endif
c
c  Compute number of steps based on that
c
      nsteps=arclng*mooth+1
c
c  Loop for points along curve
c
      do 1000 i=0, nsteps
        t=real(i)/nsteps
c
c  Compute point on first arc
c
        if(rad1.ge.0) then
          angle=ang1(1)+t*(dang11)
          r1=cent1(1)+rad1*cos(angle)
          r2=cent1(2)+rad1*sin(angle)
        else
          r1=ang1(1)+t*dang11
          r2=ang1(2)+t*dang12
        endif
c
c  Compute point on second arc
c
        if(rad2.ge.0) then
          angle=ang3(1)+t*dang21
          q1=cent2(1)+rad2*cos(angle)
          q2=cent2(2)+rad2*sin(angle)
        else
          q1=ang3(1)+t*dang21
          q2=ang3(2)+t*dang22
        endif
c
c  Take weighted average.  Note weighting coefficients
c
        p2(1)=r1*(1.-t)+q1*t
        p2(2)=r2*(1.-t)+q2*t
        if(i.ne.0) then
          call dotted(iline, p1, p2)
        endif
        call vmov2(p2, p1)
1000  continue
      return
      end
c______________________________________________________________
      logical function flag(k, i, j)
c$$$$ calls iflag
c  Routines below to set, unset or test flag
c  (note flag set in j, i, k order,
c   so in order of increasing j, then by i, then k)
c
ce
      nn=iflag(k, i, j)
      if (nn.gt.0) then
        flag=.true.
      else
        flag=.false.
      endif
      return
      end
c______________________________________________________________
      subroutine vmov2(a, b)
c$$$$ calls nothing
      real a(2),b(2)
c
ce
      b(1)=a(1)
      b(2)=a(2)
      return
      end
c______________________________________________________________
      function iflag(k, i, j)
c$$$$ calls nothing
c Return index of flag by binary search
c (note flag set in i, j, k order,
c  so in order of increasing i, then by j, then k)
      parameter (maxpts=350, maxcnt=200, maxcmp=3000, maxflg=20000)
      common /flags/ nflag, kntflg(maxflg)
c
c (note: unset flag set to -flag, so binary search still possible)
ce
      kflag=i*100000 + j*10 + k
      iflag=-1
      low=1
      khigh=nflag
      do 10 nn=1, nflag
        if (low.gt.khigh) goto 999
        mid=(low + khigh)/2
        if (kflag .lt. abs(kntflg(mid))) then
          khigh=mid - 1
        else if (kflag .gt. abs(kntflg(mid))) then
          low=mid + 1
        else
          goto 20
        endif
 10   continue
      goto 999
 20   continue
      if (kntflg(mid).gt.0) iflag=mid
 999  continue
      return
      end
c______________________________________________________________
      real function parab(r1, s1, t1, val) 
c$$$$ calls nothing
c  Suppose a quadratic function f has f(-1)=r1, f(0)=s1, and f(1)=t1,
c  with val between s and t.  Where between 0 and 1 does the parabola
c  assume the value "val"?  Rewritten by Bob Parker, 2013 to
c  eliminate dependance on external scale factor.
c
      r=r1-val
      c=s1-val
      t=t1-val
      a=(r+t)/2.0-c
      if (a.eq. 0.0) then 
        parab=c/(c-t)
        return
      endif
      b=(t-r)/2.0
      disc=sqrt(b**2 - 4.0*a*c)
      parab = -(b + sign(disc, b))/(2.0*a)
      if (parab.ge. 0.0 .and. parab.le. 1.0) return
      parab= c/(a*parab)
      return
      end
c______________________________________________________________
      subroutine line(kind, nxy, x, y)
c$$$$ calls plot
c  Given the 2-vectors of pen positions  (x(i), y(i)), i=1, nxy,
c  Draws a straight line interpolation along the sequence.
c  kind=0,1 the line is a single continuous stroke
c      =2   the line is drawn heavily by going over it twice
c      =3-9 a dashed line is drawn
c      =10  a hatched line is drawn, prongs downslope
c      =11  a hatched line, prongs uphill
c
      dimension x(nxy),y(nxy),dash(2)
      common /hatch/ klosed,lftrt
c
      data iup/3/, idn/2/, epsiln/.002/
c
c
c  Move to the start of the course
ce
       call plot(x(1), y(1), iup)
c
c  Select kind of line
      iroute=min(3, max(1,kind)) + kind/10
      goto (1000, 2000, 3000, 4000), iroute
c
c  Regular line
c
 1000  do 1100 i=1, nxy
         call plot(x(i), y(i), idn)
 1100  continue
       return
c
c  Draws a thickened line.  Runs over the vector twice.
c  First in the increasing index direction.
c
 2000 do 2100 i=2, nxy
        tad=epsiln/sqrt((x(i)-x(i-1))**2 + (y(i)-y(i-1))**2 + 1e-7)
        epx=tad*(y(i) - y(i-1))
        epy=tad*(x(i-1) - x(i))
        call plot(x(i-1)+epx, y(i-1)+epy, idn)
        call plot(x(i) + epx, y(i) + epy, idn)
 2100 continue
c
c  Now reverse the course, but not precisely
      do 2200 i=nxy, 2, -1
        tad=-epsiln/sqrt((x(i)-x(i-1))**2 + (y(i)-y(i-1))**2 + 1e-7)
        epx=tad*(y(i) - y(i-1))
        epy=tad*(x(i-1) - x(i))
        call plot(x(i) + epx, y(i) + epy, idn)
        call plot(x(i-1)+epx, y(i-1)+epy, idn)
 2200 continue
c
      return
c
c  Draw a dashed line with visible segments  dash(1) inches long and
c  gaps  dash2  long. These are computed from  kind >= 3.
c
 3000 ink=3
      r=0.0
      dash(1)=0.2/(kind-2)**2
      dash(2)=min(0.1, 2.0*dash(1))
      do 3500 i=2, nxy
        s=sqrt((x(i)-x(i-1))**2 + (y(i)-y(i-1))**2) + 1e-7
        nseg=2.0 + 2.0*s/(dash(1) + dash(2))
        do 3100 l=1, nseg
          q=min(1.0, r/s)
          call plot(x(i-1) + q*(x(i)-x(i-1)), y(i-1) + q*(y(i)-y(i-1)),
     $    ink)
          if (s .lt. r) goto 3200
          ink=5 - ink
          r=dash(ink-1) + r
 3100   continue
 3200   r=r - s
 3500 continue
c
      return
c
c  Hatched contour, prongs downhill (kind=10), or uphill (kind=11)
 4000 r=0.0
      stroke=0.1
      prong=0.2*stroke*lftrt*sign(1.0, 10.5-kind)
      ssq=(1-klosed)*stroke**2
      ink=iup
      do 4500 i=2, nxy
        s=sqrt((x(i)-x(i-1))**2 + (y(i)-y(i-1))**2) + 1e-7
        nseg=2.0 + s/stroke
        dx= (y(i)-y(i-1))*prong/s
        dy=-(x(i)-x(i-1))*prong/s
        do 4100 l=1, nseg
          q=min(1.0, r/s)
          xq=x(i-1) + q*(x(i)-x(i-1))
          yq=y(i-1) + q*(y(i)-y(i-1))
          call plot(xq, yq, ink)
          ink=idn
          if (s .lt. r) goto 4200
          if ((xq-x(1))**2+(yq-y(1))**2.gt.ssq .and.
     $        (xq-x(nxy))**2+(yq-y(nxy))**2.gt.ssq) then
            call plot(xq+dx, yq+dy, idn)
            call plot(xq, yq, idn)
          endif
          r=stroke + r
 4100   continue
 4200   r=r - s
 4500 continue
c
      end
c_____________________________________________________________
      subroutine mabel(val, len, label)
c$$$$ calls nothing
c  For a real number  val  assumed to be a decimal fraction rounded to
c  6 significant figures (or fewer), makes an equivalent character
c  string  label  of minimal length  len.
      character*8 form,local*20,label*(*)
      double precision ten, v
c
      data ten/10.0d00/
c
c  Find the number of significant figures in val
ce
      nexp=0.0
      if (val .ne. 0.0) nexp=int(500.01 + log10(abs(val))) - 500
      v=int(abs(val)/ten**(nexp-5) + 0.5)
      nres=nexp - 5
      do 1100 ipow=0, 4
        if (mod(int(v/ten**ipow+ 0.5), 10) .ne. 0) goto 1110
        nres=nres + 1
 1100 continue
 1110 nfld=3 + max(nexp, -nres, nexp-nres)
c
c  Build a format for writing the string
c
      ndpt=min(9,max(0, -nres))
      write(form, '( 2h(f, i2, 1h., i1, 1h))') nfld,ndpt
      if (ndpt .eq. 0) nfld=nfld - 1
c
c  Construct a character string of minimal length for the label
c
      write(local, form) val
c  Remove leading blanks
      is=1
      do 1200 i=1, nfld
        if (local(i:i) .eq. ' ') is=1 + is
 1200 continue
      label=local(is:nfld)
      len=nfld - is + 1
      return
      end
c________________________________________________________________
      subroutine letter(x, y, height, theta, text)
c$$$$ calls chrcod plot glyph
c  Generates text and symbol in a variety of fonts to be plotted
c  on a pen-plotter with the standard calcomp subroutine 'plot'.
c    1) four Hershey letter fonts--simplex,complex,italic, and duplex--
c       are provided in upper and lower case Roman
c    2) two Hershey letter fonts--simplex and complex--are provided in
c       upper and lower case Greek
c    3) 47 special mathematical symbols, e.g. integral sign,del, are
c       provided
c    4) super- and sub-scripting is possible within a character string
c       without separate calls to symbol
c
c    Change of font is made by enclosing the name of the font in upper
c  case in backslashes, e.g \simplex\.  Three letters suffice to
c  specify the font.  Simplex is the default font on the initial call
c  to symbol.  A font remains in effect until explicitly changed.
c  Math font \$$\ forces italic letters, retains previous font for other
c  characters.
c    Super- or sub-scripting is accomplished by enclosing the expression
c  to be super- or sub-scripted in curly brackets and preceding it by
c  sup or sub.  The closing curly bracket terminates the super- or
c  sub-scripting and returns to normal character plotting.  Note that
c  super- and sub-script letters are plotted with a different character
c  size.
c    Greek letters are drawn by enclosing the English name of the
c  letter in backslashes, e.g. \alpha\.  The case of the first letter
c  determines the case of the Greek letter.  The closing backslash must
c  be included.
c    The special graphical symbols (box, triangle, etc) can be included
c  text by enclosing the appropriate symbol integer+2000 in backslashes.
c  See subroutine glyph for the special symbol integers.
c    Any regular symbol may be drawn by enclosing the symbol number+1000
c  in backslashes.  This is the only way to call some symbols, notably
c  special mathematical symbols.
c  These regular symbols have the following integer codes
c   1-26   upper case Roman simplex
c  27-52   lower case Roman simplex
c  53-72   simplex numbers and symbols
c  73-96   upper case Greek simplex
c  97-120  lower case Greek simplex
c  121-146 upper case Roman complex
c  147-172 lower case Roman complex
c  173-192 complex numbers and symbols
c  193-216 upper case Greek complex
c  217-240 lower case Greek complex
c  241-266 upper case Roman italic
c  267-292 lower case Roman italic
c  293-312 italic numbers and symbols
c  313-338 upper case Roman duplex
c  339-364 lower case Roman duplex
c  365-384 duplex numbers and symbols
c  385-432 special mathematical symbols
c
c  Symbol parameters taken from N.M.Wolcott, Fortran IV enhanced
c  Character Graphics, NBS
c
c  A. Chave IGPP/UCSD Aug 1981, improved by R. L. Parker
c
c  x, y   are coordinates in inches from current origin to 
c         lower left corner of 1st character to be plotted.  If either
c         is set to 999.0 the saved next character position is used.
c  height is character height in inches
c  text   is an character string containing text to be plotted
c  theta  is positive angle w.r.t. the x-axis in degrees
c
c  Programmed in FORTRAN-77
c
      character*(*) text
      integer istart(432),isstar(22),symbcd(4711),ssymbc(128)
      real width(432),supsub(2),raise(20)
      common /ofset/ ioff,just1,just2
      common /ajust/ nchr,ichr(132)
      common /ialph/ width,symbcd,istart,ssymbc,isstar
      common /crown/ yx(432)
      save xo,yo
      data rad/.017453292/,is/1/
      data factor/0.75/,  supsub/0.50,-0.50/,  iup/3/
c  ichr(j) contains the symbol number of the jth symbol or a
c  code to indicate space (1000),begin super-scripting (1001),
c  begin sub-scripting (1002), or end super/sub-scripting (1003),
c  or back-space (1004).
c  1005: put hat on previous character; or 1007: tilde; or 1008: bar.
c  istart(ichr(j)) contains the address in symbol of the jth
c  character.  symbcd contains pen instructions stored in a
c  special format.  isstar and ssymbc contain addresses and pen
c  instructions for special centered symbols.  width contains
c  widths of the characters.
c
c  ixtrct gets nbits from iword starting at the nstart bit from the
c  right in an array of 32-bit integers.
      ixtrct(nstart,nbits,iword)=mod(iword/(2**(nstart-nbits)),
     $                           2**nbits)+((1-sign(1,iword))/2)*
     $                           (2**nbits-min(1,mod(-iword,
     $                           2**(nstart-nbits))))
c
ce
      ntext=len(text)
      n=ntext
      yoff=0.0
      si=sin(rad*theta)
      co=cos(rad*theta)
      high=height
      scale=high/21.0
      if (scale.eq.0.) return
      if (x.ge.999.0) then
        xi=xo
      else
        xi=x
      endif
      if (y.ge.999.0) then
        yi=yo
      else
        yi=y
      endif
c  Plot a character string.
c  First find pointer array  ichr  containing the starts of characters-
c  but only if  just1 and just2  are not 1, when  ichr is assumed
c  correctly transmitted through common /ajust/.
      if (just1.ne.1 .or. just2.ne.1) 
     $ call chrcod(text, ntext, ichr, nchr)
      just2=2
      oldwid=0.0
      ik=1
      l=1
      rscale=scale
c  Plot each character
      do 100 i=1,nchr
         ic=ichr(i)
c  ic < 0:  Change the character height
         if (ic.le.0) then
            rscale=-0.01*rscale*ic/high
            high  =-0.01*ic
c  ic = 1000: Plot a space
         elseif (ic.eq.1000)  then
           xi=xi + 20.*rscale*co
           yi=yi + 20.*rscale*si
           xo=xi
           yo=yi
           call plot(xi, yi, iup)
c  2000 <= ic <= 2021: Special graphics symbol code 2000-2021
          elseif (ic.ge.2000 .and. ic.le.2021) then
            hgt=20.0*rscale
            xo=xi + hgt*co
            yo=yi + hgt*si
            call glyph(0.5*(xi+xo-hgt*si), 0.5*(yi+yo+hgt*co),
     $      0.9*hgt, ic -2000)
            xi=xo
            yi=yo
c  ic = 1001 or 1002: Begin super-scripting or sub-scripting
         elseif (ic.eq.1001 .or. ic.eq.1002) then
            raise(l)=supsub(ic-1000)*high*rscale/scale
            rscale=factor*rscale
            yoff=raise(l)+yoff
            l=l+1
c  ic =1003: End super/sub-scripting
         elseif (ic.eq.1003) then
           rscale=rscale/factor
           l=l-1
           yoff=yoff-raise(l)
c  ic = 1004:  Backspace - use width of previous letter in oldwid
         elseif (ic.eq.1004) then
           xi=xi - co*oldwid
           yi=yi - si*oldwid
           xo=xi
           yo=yi
c  ic = 1005: Decorate previous character with a hat
c  ic = 1007: Decorate previous character with a tilde
c  ic = 1008: Decorate previous character with a bar
         elseif (ic.eq.1005 .or. ic.eq.1007 .or. ic.eq.1008) then
           if (ic .eq. 1005) is=istart(408)
           if (ic .eq. 1007) is=istart(424)
           if (ic .eq. 1008) is=istart( 64)
           yyoff=yoff + rscale*(int(yx(ik)) - 2.5)
           dxhat=oldwid - rscale*(mod(100*yx(ik),100.0)- 9.0)
           if (ic .eq. 1008) dxhat=1.2*dxhat - 0.2*oldwid
           ib=30
   60      ipen=ixtrct(ib,3,symbcd(is))
           if (ipen.eq.0) goto 100
           ix=ixtrct(ib-3,6,symbcd(is))
           iy=ixtrct(ib-9,6,symbcd(is))
           xx=(ix-10)*rscale - dxhat
           yy=(iy-11)*rscale + yyoff
           call plot(xi+co*xx-si*yy, yi+co*yy+si*xx, ipen)
           ib=45-ib
           if (ib.eq.30) is=is+1
           goto 60
c  ic = 1006: Space forward about 1/3 regular amount
         elseif (ic.eq.1006)  then
           xi=xi + 6.66*rscale*co
           yi=yi + 6.66*rscale*si
           xo=xi
           yo=yi
         else
c
c  ic = everything else: Plot a single symbol
           is=istart(ic)
           ik=ic
           ib=30
   70      ipen=ixtrct(ib,3,symbcd(is))
           if (ipen.eq.0) then
             xi=xi+co*rscale*width(ic)
             yi=yi+si*rscale*width(ic)
             xo=xi
             yo=yi
             oldwid=width(ic)*rscale
             goto 100
           endif
           ix=ixtrct(ib-3,6,symbcd(is))
           iy=ixtrct(ib-9,6,symbcd(is))
           xx=(ix-10)*rscale
           yy=(iy-11)*rscale+yoff
           call plot(xi+co*xx-si*yy, yi+co*yy+si*xx, ipen)
           ib=45-ib
           if (ib.eq.30) is=is+1
           goto 70
         endif
  100   continue
      return
      end
c_______________________________________________________________________
      subroutine chrcod(text, ntext, ichr, nchr)
c  Given text string in text, ntext characters
c  returns ichr containing nchr symbol numbers or codes for
c  space (1000), begin superscripting (1001), begin
c  subscripting (1002), or end super/sub-scripting (1003)
c  change of font commands are decoded and executed internally
c
      common /ofset/ ioff,just1,just2
      character*1 bkslsh
      character*(*) text
      integer ichr(132),irlu(95),iilu(95),iglu(26)
c      data ioff/0/
c  irlu is a look-up table for Roman characters arranged by
c  integer value for the ascii character set with an
c  offset to remove the 31 nonprinting control characters.
c  irlu returns with the symbol number or, if no symbol
c  exists, the code for space = 1000
      data irlu/1000,416, 428,411,72,418,419, 432,67,68,69,63,70,
     $          64,71,65,53,54,55,56,57,58,59,60,61,62,414,415,
     $          385,66,386,417,407,1,2,3,4,5,6,7,8,9,10,11,12,13,
     $          14,15,16,17,18,19,20,21,22,23,24,25,26,409,1000,
     $          410, 408,1000,1000,27,28,29,30,31,32,33,34,35,36,
     $          37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,
     $          405,427,406,424/
c  iilu is a look-up table for italic characters only.  It is
c  identical to irlu with four italic special symbols substituted
c  for regular ones.
      data iilu/1000,422, 428,411,72,418,419, 432,67,68,69,63,70,
     $          64,71,65,53,54,55,56,57,58,59,60,61,62,420,421,
     $          385,66,386,423,407,1,2,3,4,5,6,7,8,9,10,11,12,13,
     $          14,15,16,17,18,19,20,21,22,23,24,25,26,409,1000,
     $          410,1000,1000,1000,27,28,29,30,31,32,33,34,35,36,
     $          37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,
     $          405,427,406,424/
c  iglu is a look-up table for Greek characters arranged by the
c  integer value of their Roman expression with a=1, b=2, etc.
c  Ambiguous cases give 25 for epsilon or eta, 26 for omega or
c  omicron, 27 for phi,pi,or psi, and 28 for tau or theta.  Additional
c  letters must be checked for in these cases. A value of 50 is returned
c  for those Roman letters which have no corresponding Greek letter.
      data iglu/1,2,22,4,25,50,3,50,9,50,10,11,12,13,26,27,50,17,18,
     $             28,20,50,50,14,50,6/
      data math/-1/
c
c  Due to special nature of backslash in Unix, generate it from ASCII.
      bkslsh=char(92)
c
c Finds length of string with blanks trimmed from right end
      do 10 n=ntext,1,-1
        if (text(n:n) .ne. ' ') goto 15
 10   continue
      nchr=0
      return
 15   nt=n
c  Scan text character by character
      k=1
      j=1
c  k is current address of character in text
c  j is index of next symbol code in ichr
   20 if (k.gt.n) then
        nchr=j-1
        return
      endif
      if (text(k:k).ne. bkslsh) then
c  Roman character or keyboard symbol
        if (text(k:k).eq.'}') then
c  Check for closing curly bracket-if found, return 1003
           ichr(j)=1003
           j=j+1
           k=k+1
           goto 20
         endif
c 
c  The function ichar returns integer ASCII value of character 
c  (this is an intrinsic function in FORTRAN-77) offset by nonprinting
c   characters to get entry in look-up table
        ic=ichar(text(k:k))-ichar(' ')+1
        joff=0
c  Force italic font on alphabetic chars when math >= 0.  Save old
c  offest in math
        if (math .ge.0 .and. (34 .le.ic .and. ic.le. 59 .or.
     $      66.le.ic .and. ic.le.91)) then
          math=ioff
          joff=1
          ioff=240
        endif
c  Nonprinting control character-error return
        if (ic.le.0) then
          ichr(j)=1000
c  Not italic font
        elseif (ioff.ne.240) then
           ichr(j)=irlu(ic)
        else
c  Italic font
          ichr(j)=iilu(ic)
        endif
c  Add offset for font if not a special symbol
        if (ichr(j).lt.385) ichr(j)=ichr(j)+ioff
        j=j+1
        k=k+1
        if (joff .eq. 1) ioff=math 
        goto 20
      else
c
c  Backslash found
         k=k+1
c
c  Check for sundry mathematical symbols, substitute 4-digit codes
c  Regular theta: 'rgth' 
         if (text(k:k+3) .eq. 'rgth') then
           text(k:k+3)='1404'
c  Regular phi:  'rgph'
         elseif (text(k:k+3) .eq. 'rgph') then
           text(k:k+3)='1431'
c  Infinity: 'infi'
         elseif (text(k:k+3) .eq. 'infi') then
           text(k:k+3)='1395'
c  Integral: 'inte'
         elseif (text(k:k+3) .eq. 'inte') then
           text(k:k+3)='1393'
c  Partial: 'part'
         elseif (text(k:k+3) .eq. 'part') then
           text(k:k+3)='1387'
c  Grad:  'grad'
         elseif (text(k:k+3) .eq. 'grad') then
           text(k:k+3)='1388'
c  Times: 'time' 
         elseif (text(k:k+3) .eq. 'time') then
           text(k:k+3)='1398'
c  Degrees: 'degr' 
         elseif (text(k:k+3) .eq. 'degr') then
           text(k:k+3)='1429'
c  Square-root: 'sqrt'
         elseif (text(k:k+3) .eq. 'sqrt') then
           text(k:k+3)='1402'
         endif
c
c
c  Check next four characters for 4-digit number
         read (text(k:k+3), '(i4)', err=50) number
c  Number found - valid numbers are
c         0<n<500  or  999<n<1433  or 1999<n<2022
c
c  3-digit number changes letter height
        if (number.gt.0 .and. number.lt.500) then
          ichr(j)=-number
c  Special graphics symbol for glyph
        elseif (number.gt.1999 .and. number.le.2021) then
           ichr(j)=number
c  Valid symbol code
         elseif (number.gt.999 .and. number.lt.1433) then
           ichr(j)=number - 1000
         else
c  Not recognized - plot a blank
           ichr(j)=1000
         endif
         j=j+1
c  Move beyond closing backslash-ignore extra characters
c  Function index returns offset of second substring in first
c  Returns 0 if substring not found
        l=index(text(k:nt), bkslsh)
        if (l.eq.0) then
           k=nt+1
         else
            k=k+l
         endif
         goto 20
c
c  Not a number in the backslashes
   50  continue
c
c  Check for font change command
c  Simplex font
       if (text(k:k+2).eq.'SIM'.or.text(k:k+2).eq.'sim') then
         ioff=0
         math=-1
c  Complex font
       elseif (text(k:k+1).eq.'CO'.or.text(k:k+1).eq.'co') then
         ioff=120
         math=-1
c  Italic font
       elseif (text(k:k+1).eq.'IT'.or.text(k:k+1).eq.'it') then
         ioff=240
         math=-1
c  Duplex font
       elseif (text(k:k+1).eq.'DU'.or.text(k:k+1).eq.'du') then
         ioff=312
         math=-1
c  Math font (=> italics for letters, previous font otherwise)
       elseif (text(k:k).eq.'$') then
         math=ioff
c  Found the back-space code
       elseif (text(k:k+1).eq.'BS'.or.text(k:k+1).eq.'bs') then
         ichr(j)=1004
         j=j+1
         k=k+3
         go to 20
c  Check for super/sub-script command
       elseif (text(k:k+3).eq.'SUP{'.or. text(k:k+3).eq.'sup{') then
c  Begin superscripting
         ichr(j)=1001
         j=j+1
         k=k+4
         goto 20
       elseif (text(k:k+3).eq.'SUB{'.or. text(k:k+3).eq.'sub{') then
c  Begin subscripting
         ichr(j)=1002
         j=j+1
         k=k+4
         goto 20
c  Put a hat over the previous character
       elseif(text(k:k) .eq. '^') then
         ichr(j)=1005
         j=j+1
         k=k+2
         goto 20
c  Put a tilde over the previous character
       elseif(text(k:k) .eq. '~') then
         ichr(j)=1007
         j=j+1
         k=k+2
         goto 20
c  Put a bar over the previous character
       elseif(text(k:k) .eq. '-') then
         ichr(j)=1008
         j=j+1
         k=k+2
         goto 20
c  Space forward 1/3 of regular space 
       elseif(text(k:k) .eq. ' ') then
         ichr(j)=1006
         j=j+1
         k=k+2
         goto 20
       else
c  Greek character or invalid character
      ic=ichar(text(k:k))
      igoff=min(ioff, 120)
      if (ic.ge.ichar('A') .and. ic.le.ichar('Z')) then
c  Upper case
c  Greek character or invalid character
        igr=72
        ico=ichar('A')-1
      elseif (ic.ge.ichar('a') .and. ic.le.ichar('z')) then
c  Lower case
        igr=96
        ico=ichar('a')-1
      else
c  Not a letter-error return
        ichr(j)=1000
        j=j+1
        l=index(text(k:nt), bkslsh)
        if (l.eq.0) then
          k=nt+1
        else
          k=k+l
        endif
        goto 20
      endif
c  Look up the character
      ig=iglu(ic-ico)
      if (ig.lt.25) then
c  Unambiguous Greek letter
        ichr(j)=ig+igr+igoff
      elseif (ig.eq.25) then
c  Epsilon or eta
        ib=ichar(text(k+1:k+1))-ico
        if (ib.eq.16) then
c  Epsilon
          ichr(j)=5+igr+igoff
        elseif (ib.eq.20) then
c  Eta
          ichr(j)=7+igr+igoff
        else
c  Not a Greek character--error return
          ichr(j)=1000
        endif
      elseif (ig.eq.26) then
c  Omega or omicron
        ib=ichar(text(k+1:k+1))-ico
        if (ib.ne.13) then
c  Not a Greek character-error return
          ichr(j)=1000
        else
          ic=ichar(text(k+2:k+2))-ico
          if (ic.eq.5) then
c  Omega
            ichr(j)=24+igr+igoff
          elseif (ic.eq.9) then
c  Omicron
            ichr(j)=15+igr+igoff
          else
c  Not a Greek character-error return
            ichr(j)=1000
          endif
        endif
        elseif (ig.eq.27) then
c  Phi,pi, or psi
          ib=ichar(text(k+1:k+1))-ico
          if (ib.eq.8) then
c  Phi
            ichr(j)=21+igr+igoff
            elseif (ib.eq.9) then
c  Pi
              ichr(j)=16+igr+igoff
            elseif (ib.eq.19) then
c  Psi
              ichr(j)=23+igr+igoff
            else
c  Not a Greek character-error return
              ichr(j)=1000
            endif
          elseif (ig.eq.28) then
c  Tau or theta
            ib=ichar(text(k+1:k+1))-ico
          if (ib.eq.1) then
c  Tau
            ichr(j)=19+igr+igoff
          elseif (ib.eq.8) then
c  Theta
            ichr(j)=8+igr+igoff
          else
c  Not a Greek character-error return
            ichr(j)=1000
          endif
        else
c  Not a Greek character-error return
          ichr(j)=1000
        endif
          j=j+1
      endif
      l=index(text(k:nt), bkslsh)
      if (l.eq.0) then
        k=nt+1
      else
        k=k+l
      endif
      goto 20
      endif
      return
      end
c_______________________________________________________________________
      subroutine justy(s, height, text)
c$$$$ calls chrcod
c  Given the text string  text  with  ntext  characters,
c  height  high in inches,  this routine
c  gives 4 distances in inches, all from the left end of the string -
c  s(1)  to the left edge of the 1st nonblank character
c  s(2)  to the center of the the string, blanks removed from the ends
c  s(3)  to the right edge of the last nonblank character
c  s(4)  to the right edge of the last character of the string.
      character*(*) text
      dimension s(4),ipower(3)
      common /ialph/ width(432),sym1(4711),is1(432),sym2(128),is2(22)
      common /ofset/ ioff,just1,just2
      common /ajust/ nchr,ichr(132)
      data ipower/1,1,-1/,  factor/0.75/
c
      if (text.eq. ' ') then
        do 1010 j=1, 4
          s(j)=0.0
 1010   continue
        return
      endif
c
      ntext=len(text)
      scale=height/21.0
      call chrcod(text, ntext, ichr, nchr)
c
c  Count leading blanks
      do 1100 lead=1, nchr
        if (ichr(lead) .ne. 1000) go to 1110
 1100 continue
      lead=nchr
 1110 s(1)=20.0*scale*(lead-1)
      s(3)=s(1)
c
c  Sum the widths of the remaining text, recalling that trailing blanks
c  were lopped off by  chrcod
      oldwid=0.0
      do 1200 i=lead, nchr
        l=ichr(i)
c  Change character height
        if (l .lt. 0) scale=-0.01*l/21.0
c  Regular character
        if (l.gt.0 .and. l.lt.1000) then
          oldwid=width(l)*scale
          s(3)=s(3) + oldwid
        endif
c  Glyph symbol or space
        if (l.eq.1000.or.l.ge.2000) s(3)=s(3) + 20.0*scale
c  Sub/super script size change
        if (l.ge.1001.and.l.le.1003) scale=scale*factor**ipower(l-1000)
c  Backspace
        if (l.eq.1004) s(3)=s(3) - oldwid
 1200 continue
c
c  Add on width of surplus trailing blanks
      s(4)=s(3) + 20.0*scale*(ntext-nchr)
c
c  Find center of nonblank text
      s(2)=(s(1) + s(3))/2.0
      just2=1
      return
      end
c_______________________________________________________________________
      subroutine glyph(x, y, height, nglyph)
c$$$$ calls plot
c  Draws centered symbols.
c  At the coordinates  (x, y)  in inches measured from the current
c  origin, draws a single centered symbol, with height  height
c  and identity specified by  nglyph  as follows -
c  0 square,    1 triangle, 2 octagon,        3  diamond, 4  plus,
c  5 asterisk,  6 cross,    7 slashed square, 8 up-arrow, 9 hourglass,
c  10 campstool,11 hexagon, 12 Y,         13 |,       14 star of David
c  15 dot,      16 sm circ, 17 circle,    18 lg circ, 19 filled sm circ
c  20 filled sm square,     21 filled sm triangle.
c  Data for these glyphs come through common /ialph/.
      integer istart(432),isstar(22),symbcd(4711),ssymbc(128),map(22)
      common /ialph/ width(432),symbcd,istart,ssymbc,isstar
      data map/0,2,1,5,3,11,4,10,6,12,7,8,9,13,14,15,
     $16,17,18,19,20,21/
c  ixtrct gets nbits from iword starting at the nstart bit from the
c  right.  An alternative version using only arithmetic operations
c  in subroutine letter.
      ixtrct(nstart,nbits,iword)=mod(iword/(2**(nstart-nbits)),
     $                           2**nbits)+((1-sign(1,iword))/2)*
     $                           (2**nbits-min(1,mod(-iword,
     $                           2**(nstart-nbits))))
c
c  Plot a single special centered symbol
      scale=height/21.0
       ia=nglyph + 1
c  Re-orders special symbols to Prime convention.
      ia=map(ia) + 1
      is=isstar(ia)
       ib=30
c  Unpack the pen position from ssymbc.  ipen=0  means quit.
   20   ipen=ixtrct(ib, 3, ssymbc(is))
       if (ipen.eq.0) return
c  Unpack and scale coordinates.
       ix=ixtrct(ib-3, 6, ssymbc(is))
       iy=ixtrct(ib-9, 6, ssymbc(is))
       call plot(x+scale*(ix-32), y+scale*(iy-32), ipen)
       ib=45-ib
       if (ib.eq.30) is=is+1
       goto 20
      end
c_______________________________________________________________________
      blockdata hatdat
c  yx  contains character heights and centers of topmost portion
      common /crown/ yx(432)
      data (yx(j),j=1,162)/
     $21.09,21.07,21.11,21.06,21.08,21.08,21.11,21.11,21.04,
     $21.12,21.11,21.04,21.12,21.09,21.10,21.07,21.10,21.07,
     $21.10,21.08,21.11,21.09,21.12,21.10,21.09,21.12,14.11,
     $21.04,14.10,21.15,14.10,21.09,14.11,21.04,22.04,22.06,
     $21.04,21.04,14.14,14.08,14.09,14.08,14.11,14.08,14.08,
     $21.05,14.10,14.08,14.11,14.08,14.08,14.10,21.10,21.11,
     $21.10,21.10,21.13,21.10,21.11,21.12,21.09,21.10,18.13,
     $9.13,25.2,12.13,25.11,25.03,15.08,2.05,2.05,25.10,
     $21.09,21.07,21.08,21.09,21.08,21.12,21.11,21.10,21.04,
     $21.11,21.09,21.12,21.09,21.09,21.10,21.11,21.07,21.07,
     $21.08,21.09,21.10,21.10,21.11,21.10,14.12,21.13,14.09,
     $21.10,14.08,21.10,14.09,21.13,14.06,14.12,21.02,14.12,
     $14.08,21.10,14.09,14.12,14.10,14.12,14.12,14.09,14.16,
     $14.08,21.16,14.13,21.10,21.08,21.12,21.07,21.10,21.10,
     $21.12,21.12,21.06,21.10,21.11,21.06,21.12,21.10,21.11,
     $21.08,21.11,21.08,21.12,21.10,21.11,21.09,21.12,21.09,
     $21.10,21.11,14.09,21.05,14.10,21.15,14.10,21.09,14.10,
     $21.05,21.05,21.06,21.05,21.05,14.13,14.08,14.10,14.08/
      data (yx(j+162),j=1,144)/
     $14.12,14.07,14.10,21.06,14.10,14.08,14.12,14.09,14.09,
     $14.10,21.10,21.11,21.11,21.11,21.13,21.10,21.11,21.10,
     $21.10,21.10,18.13,9.13,25.2,12.13,25.11,25.03,21.08,
     $2.05,2.05,25.10,21.10,21.08,21.09,21.10,21.10,21.11,
     $21.12,21.11,21.06,21.11,21.10,21.12,21.10,22.12,21.11,
     $21.12,21.08,21.08,21.10,21.10,21.10,21.09,21.12,21.11,
     $14.12,21.14,14.12,22.11,14.08,21.11,14.10,21.15,14.06,
     $14.12,21.04,14.13,14.10,21.11,14.10,14.14,14.11,14.12,
     $14.12,14.08,14.16,14.07,21.16,14.13,21.13,21.12,21.15,
     $21.11,21.13,21.13,21.15,21.16,21.10,21.14,21.15,21.10,
     $21.16,21.14,21.13,21.12,21.13,21.12,21.16,21.14,21.14,
     $21.12,21.16,21.13,21.13,21.15,14.12,21.08,14.10,21.18,
     $14.10,21.14,14.12,21.08,21.09,21.10,21.08,21.08,14.15,
     $14.10,14.10,14.09,14.12,14.09,14.10,21.10,14.11,14.10,
     $14.15,14.12,14.11,14.11,21.13,21.14,21.13,21.13,21.17,
     $21.14,21.13,21.12,21.12,21.12,18.13,9.13,25.24,12.13/
      data (yx(j+306),j=1,126)/
     $25.15,25.09,21.10,2.03,2.03,25.14,21.10,21.07,21.11,
     $21.06,21.10,21.10,21.11,21.11,21.04,21.12,21.11,21.04,
     $21.12,21.12,21.10,21.07,21.10,21.07,21.10,21.06,21.11,
     $21.10,21.13,21.10,21.10,21.11,14.13,21.04,14.10,21.15,
     $14.10,21.10,14.13,21.04,21.04,21.04,21.04,21.04,14.12,
     $14.07,14.09,14.07,14.13,14.08,14.08,21.05,14.10,14.08,
     $14.12,14.09,14.08,14.09,21.10,21.11,21.10,21.10,21.14,
     $21.08,21.11,21.10,21.09,21.10,18.12,10.10,25.2,14.10,
     $25.10,25.03,21.08,3.06,3.05,25.09,18.2,18.04,21.11,
     $21.09,14.10,21.2,21.04,14.2,25.2,25.2,13.12,17.12,
     $17.12,16.11,18.13,25.17,10.05,25.22,25.11,21.12,25.09,
     $25.05,21.14,11.09,25.06,25.08,21.14,21.08,21.08,14.05,
     $14.05,21.05,21.09,21.13,21.10,14.06,14.06,21.08,21.14,
     $12.12,25.10,25.04,25.04,21.10,21.07,14.16,21.14,21.10/
      end
c_______________________________________________________________________
c  If your computer does not use ASCII representation for characters
c  you will need to translate them for chrcod.  Only 95 characters 
c  have significance; here are their ASCII integers:
c          32    33 !  34 "  35 #  36 $  37 %  38 &  39 '  40 (
c    41 )  42 *  43 +  44 ,  45 -  46 .  47 /  48 0  49 1  50 2
c    51 3  52 4  53 5  54 6  55 7  56 8  57 9  58 :  59 ;  60 <
c    61 =  62 >  63 ?  64 @  65 A  66 B  67 C  68 D  69 E  70 F
c    71 G  72 H  73 I  74 J  75 K  76 L  77 M  78 N  79 O  80 P
c    81 Q  82 R  83 S  84 T  85 U  86 V  87 W  88 X  89 Y  90 Z
c    91 [  92 \  93 ]  94 ^  95 _  96 `  97 a  98 b  99 c 100 d
c   101 e 102 f 103 g 104 h 105 i 106 j 107 k 108 l 109 m 110 n
c   111 o 112 p 113 q 114 r 115 s 116 t 117 u 118 v 119 w 120 x
c   121 y 122 z 123 { 124 | 125 } 126 ~ 
c  You must provide a substitute integer-valued function for  ichar
c  (say  jchar) so that, for example,  jchar('W')=87. The simplest way 
c  to do this is to set up an integer array kode mapping your
c  character set into the above scheme; thus  kode(ichar('W'))=87,
c  where  ichar  is your local instrinsic FORTRAN-77 function.
c  Then  jchar  would be defined as follows:
c        function jchar(kar)
c        character*1 kar
c        dimension kode(128)
c        data kode/ ........../
c        j=ichar(kar)
c        jchar=kode(j)
c        return
c        end
c  Finally, you must substitute your integer code for backslash in place
c  of 92 in the line defining the character variable   bkslsh   in
c  subroutine chrcd.
c_______________________________________________________________________
      blockdata blocka
c  Wolcott's blockdata statement reordered for subroutine letter.
c  The new ordering is as follows
c   1-26   upper case roman simplex
c  27-52   lower case roman simplex
c  53-72   simplex numbers and symbols
c  73-96   upper case greek simplex
c  97-120  lower case greek simplex
c  121-146 upper case roman complex
c  147-172 lower case roman complex
c  173-192 complex numbers and symbols
c  193-216 upper case greek complex
c  217-240 lower case greek complex
c  241-266 upper case roman italic
c  267-292 lower case roman italic
c  293-312 italic numbers and symbols
c  313-338 upper case roman duplex
c  339-364 lower case roman duplex
c  365-384 duplex numbers and symbols
c  385-432 special mathematical symbols
c
      integer symbcd,ssymbc
      common /ialph/ width(432),symbcd(4711),istart(432),ssymbc(128),
     $isstar(22)
      common /ofset/ ioff,just1,just2
c     data ioff,just1,just2/0,0,0/
c  Complex font standard: ioff=120
      data ioff,just1,just2/120,0,0/
      data (symbcd(j), j=1, 114)/
     $443556555,443557579,432612882,        0,433070987,433071584,
     $323987166,328083226,325854871,317404054,317400725,325723922,
     $327657165,323364299,298156032,462268125,321889760,309339231,
     $300852123,296493907,298329038,304489675,317040204,325527312,
     $        0,433070987,433071456,319792797,325953304,327788240,
     $323429900,312845195,        0,433070987,433071840,432743830,
     $432383691,        0,433070987,433071840,432743830,        0,
     $462268125,321889760,309339231,300852123,296493907,298329038,
     $304489675,317040204,325527312,327792083,327778304,433070987,
     $462432011,432744214,        0,433070987,        0,449848720,
     $312911116,306553867,298197837,294134546,        0,433070987,
     $462431122,443262731,        0,433070987,432383627,        0,
     $433070987,433071499,466625931,466626443,        0,433070987,
     $433071883,462432011,        0,443556959,300852123,296493907,
     $298329038,304489675,317040204,325527312,329885528,328050397,
     $321889760,309329920,433070987,433071584,323987166,328083225,
     $325822102,317367189,        0,443556959,300852123,296493907,
     $298329038,304489675,317040204,325527312,329885528,328050397,
     $321889760,309343631,327450624,433070987,433071584,323987166/
      data (symbcd(j), j =    115,  228)/
     $328083226,325854871,317399958,447424267,        0,460236383,
     $315630752,300917597,296592281,300688471,317367892,323593937,
     $325527116,314942603,300294990,        0,441459851,426780256,
     $        0,433070993,300360780,310748555,321267406,327722784,
     $        0,426779851,460334283,        0,428876875,449848395,
     $449849035,470820555,        0,430974667,460333899,        0,
     $426779862,308655840,309002240,460333899,430974688,430286539,
     $        0,443556555,443557579,432612882,        0,433070987,
     $433071584,323987166,328083226,325854871,317404054,317400725,
     $325723922,327657165,323364299,298156032,433070987,433071776,
     $        0,443556555,443557579,426092235,        0,433070987,
     $433071840,432743830,432383691,        0,460333899,430974688,
     $430286539,        0,433070987,462432011,432744214,        0,
     $443556959,300852123,296493907,298329038,304489675,317040204,
     $325527312,329885528,328050397,321889760,309343382,319488000,
     $433070987,        0,433070987,462431122,443262731,        0,
     $443556555,443557579,        0,433070987,433071499,466625931,
     $466626443,        0,433070987,433071883,462432011,        0,
     $428877472,436938134,428189323,        0,443556959,300852123/
      data (symbcd(j), j =    229,  342)/
     $296493907,298329038,304489675,317040204,325527312,329885528,
     $328050397,321889760,309329920,433070987,462432011,433071904,
     $        0,433070987,433071584,323987166,328083225,325822102,
     $317367189,        0,428877014,293974816,324023051,323321856,
     $441459851,426780256,        0,428712733,296723360,303047775,
     $307143897,308655771,323921503,319825312,313500957,309100544,
     $445654283,441295834,298623831,296362898,300459152,315106897,
     $323561172,325822105,321725851,307068928,430974667,430286560,
     $        0,447751499,428680026,298623957,302621778,310945169,
     $321463955,325756697,330114970,        0,430285899,298394454,
     $296559517,303015136,313533983,323921626,325789330,317040331,
     $        0,455910987,455812568,313304217,302785430,296330065,
     $298263564,306554187,317072974,        0,433070987,432743448,
     $307012953,317466198,323593873,321332684,312845451,302392206,
     $        0,455812568,313304217,302785430,296330065,298263564,
     $306554187,317072974,        0,456140363,455812568,313304217,
     $302785430,296330065,298263564,306554187,317072974,        0,
     $430548563,321562135,317465945,307012632,298525523,296264590,
     $302392459,312845772,321323008,445654176,303014876,300266265/
      data (symbcd(j), j =    343,  456)/
     $309100544,455910985,318973381,312616068,302167638,317465945,
     $307012632,298525523,296264590,302392459,312845772,321323008,
     $433070987,432710744,309110169,319563349,321224704,430973855,
     $300950433,296760217,298156032,435168287,305144865,300954649,
     $302261189,295838404,        0,433070987,453813135,441034315,
     $        0,433070987,        0,432841611,432710744,309110169,
     $319563349,321238613,327952281,338471128,344631563,        0,
     $432841611,432710744,309110169,319563349,321224704,441230360,
     $298525523,296264590,302392459,312845772,321332881,323593814,
     $317465945,307003392,432841604,432743448,307012953,317466198,
     $323593873,321332684,312845451,302392206,        0,455910980,
     $455812568,313304217,302785430,296330065,298263564,306554187,
     $317072974,        0,432841611,432645078,304882905,315392000,
     $453715416,311207001,298591062,298460179,313075153,319268366,
     $317072651,304456588,296157184,435168207,302392459,310752025,
     $309100544,432841615,300295243,310748556,321369689,321224704,
     $428647563,453813387,        0,430744651,447521867,447522379,
     $464299595,        0,430745099,453813067,        0,428647563,
     $453813387,302228357,293741252,        0,453813067,430745113/
      data (symbcd(j), j =    457,  570)/
     $430286347,        0,443327576,300622740,296264526,298198027,
     $306554124,317171282,325789465,443327833,315368918,321332876,
     $325429003,        0,449848607,307143705,300622738,296100612,
     $449848864,323954331,321693208,315335895,443262294,317335058,
     $319268301,314975499,306553868,300327824,        0,426451800,
     $300721177,306980055,311043344,308655833,323692116,308651079,
     $302120960,447521945,302785430,296330064,298230732,304456907,
     $312878542,319333908,317433177,309175453,307209440,313533919,
     $321814528,451650968,311207001,300688342,302654675,443130834,
     $296231758,298198027,308651340,317128704,445654175,305079389,
     $307111259,319665691,311206999,298459985,296199053,302359753,
     $310617349,308421700,302186496,426418967,298624025,304882774,
     $302588811,436806806,311174553,319596183,323626575,314703872,
     $426418967,298624025,304882774,302556174,304489611,310748556,
     $319268433,323626713,325985951,319825312,313468252,315401750,
     $323626834,        0,437035922,296166220,298165259,306619599,
     $        0,437035787,457975385,319595928,306848787,300528595,
     $304686225,310781259,314942924,        0,426779488,300917790,
     $319141017,293961728,439132868,436904912,300328011,308651340/
      data (symbcd(j), j =    571,  684)/
     $317138514,460105298,319235596,321234635,329688975,        0,
     $430744601,300524430,296072857,321594900,315139278,302392139,
     $        0,445654175,305079389,307111259,319665499,307045401,
     $300655573,304719122,315176210,302556048,296166220,300229832,
     $310617349,306324484,        0,441230360,298525523,296231821,
     $300295243,308651340,317138449,319432151,315368729,307003392,
     $443327435,453813843,323430091,428549016,304916377,        0,
     $432645008,300327948,306554123,314975758,321431124,319530456,
     $313304281,304882646,298427012,        0,462202009,302785430,
     $296330064,298230732,304456907,312878542,319333908,317433240,
     $311197696,447521931,428549016,304916249,        0,426418967,
     $298624025,304882774,300426189,304456907,314975758,323561174,
     $325877760,441197591,298492754,296199053,300295243,310748620,
     $323430161,329918295,325887577,317433171,308749316,        0,
     $428647321,302753158,318908036,460105367,319431561,293806788,
     $        0,458237060,426418967,298624025,304882774,302556174,
     $304489675,312845836,323430161,332081113,        0,441230360,
     $298492754,296199052,300262475,308684111,449422671,314975691,
     $321234636,329754514,332048216,327974912,445653835,445654731/
      data (symbcd(j), j =    685,  798)/
     $445556363,434677265,426091595,451258187,        0,435168203,
     $437265419,428877344,326084382,330180442,327952087,319501856,
     $323987166,328083226,325854871,319501334,319497941,327821138,
     $329754381,325461515,293975574,323659476,327755535,325494412,
     $319127552,460236570,328214237,321889696,311436383,300852123,
     $296493907,298329038,304489739,314943052,325527312,445654175,
     $302949339,298591123,300426254,306586891,        0,435168203,
     $437265419,428877216,321890013,328050520,329885456,325527116,
     $314942219,449848863,323921627,327952147,325592718,319169931,
     $        0,435168203,437265419,449652114,428877600,328017632,
     $436938134,428189451,327722699,        0,435168203,437265419,
     $449652114,428877600,328017632,436938134,428188875,        0,
     $460236570,328214237,321889696,311436383,300852123,296493907,
     $298329038,304489739,314943052,325530912,307209245,300786584,
     $298427344,302457996,310752979,325433107,327530003,334069760,
     $435168203,437265419,462432011,464529227,428877024,456140832,
     $436938518,428188875,455452683,        0,435168203,437265419,
     $428877024,428188875,        0,445654287,308683851,300262220,
     $294069008,296264592,296203488,308782220,304460832,317718528/
      data (symbcd(j), j =    799,  912)/
     $435168203,437265419,464528403,447457099,445359883,428877024,
     $456140768,428188875,455452619,        0,435168203,437265419,
     $428877024,428189387,325625483,        0,435168203,437265806,
     $435168651,464528779,464529227,466626443,428876832,464529504,
     $428188811,457549899,        0,435168203,437266189,437200651,
     $462432011,428876832,456140768,428188811,        0,445654111,
     $300852123,296461140,298329038,304489739,314943052,325527312,
     $329918295,328050397,321889696,311440672,307209245,300786583,
     $298460112,302457996,310752651,319170190,325592852,327919323,
     $323921439,315621376,435168203,437265419,428877344,326084382,
     $330180441,327919318,319464469,454043295,326051612,327984855,
     $323692053,428188875,        0,445654111,300852123,296461140,
     $298329038,304489739,314943052,325527312,329918295,328050397,
     $321889696,311440672,307209245,300786583,298460112,302457996,
     $310752651,319170190,325592852,327919323,323921439,315634765,
     $304555152,310945105,317203982,321103494,327362376,329561614,
     $321201800,325297927,329515008,435168203,437265419,428877344,
     $326084382,330180442,327952087,319497238,454043295,326051612,
     $328017624,323724822,428188875,447423957,319432397,327558988/
      data (symbcd(j), j =    913,  1026)/
     $331789781,319399564,325429067,331786126,        0,458139360,
     $325920413,319792480,307241951,296657755,298623960,304850389,
     $321529554,430810073,304883158,321562260,325658318,321267083,
     $308651020,298263377,296067982,        0,443557067,445654283,
     $430973722,294659808,325920416,436577739,        0,435168209,
     $302457996,312845771,323364622,329820000,437265425,304555212,
     $312849184,309343904,336592896,430974219,433071374,460334347,
     $426779744,451946336,        0,433071243,435168400,449848459,
     $449848971,451946128,466626187,426779808,460335200,        0,
     $430974603,433071819,460333899,426779744,451946336,426091595,
     $451258187,        0,430974229,310752160,313173323,462431573,
     $426779744,454043552,438674955,        0,458236747,460333963,
     $433070938,296756960,430286539,325625483,        0,445653835,
     $445654731,445556363,434677265,426091595,451258187,        0,
     $435168203,437265419,428877344,326084382,330180442,327952087,
     $319501856,323987166,328083226,325854871,319501334,319497941,
     $327821138,329754381,325461515,293975574,323659476,327755535,
     $325494412,319127552,435168203,437265419,428877536,325920416,
     $428188875,        0,445653771,445654795,445556427,430319308/
      data (symbcd(j), j =    1027,  1140)/
     $428189451,        0,435168203,437265419,449652114,428877600,
     $328017632,436938134,428189451,327722699,        0,458236747,
     $460333963,433070938,296756960,430286539,325625483,        0,
     $435168203,437265419,462432011,464529227,428877024,456140832,
     $436938518,428188875,455452683,        0,445654111,300852123,
     $296461140,298329038,304489739,314943052,325527312,329918295,
     $328050397,321889696,311440672,307209245,300786583,298460112,
     $302457996,310752651,319170190,325592852,327919323,323921439,
     $315634841,306787865,319370390,319501461,319455232,435168203,
     $437265419,428877024,428188875,        0,435168203,437265419,
     $464528403,447457099,445359883,428877024,456140768,428188875,
     $455452619,        0,445653835,445654731,445556363,426091595,
     $451258187,        0,435168203,437265806,435168651,464528779,
     $464529227,466626443,428876832,464529504,428188811,457549899,
     $        0,435168203,437266189,437200651,462432011,428876832,
     $456140768,428188811,        0,433103708,464561948,441197651,
     $455878163,432513866,463972106,433039135,433006366,441132566,
     $441099797,432449293,432416524,        0,445654111,300852123,
     $296461140,298329038,304489739,314943052,325527312,329918295/
      data (symbcd(j), j =    1141,  1254)/
     $328050397,321889696,311440672,307209245,300786583,298460112,
     $302457996,310752651,319170190,325592852,327919323,323921439,
     $315621376,435168203,437265419,462432011,464529227,428877856,
     $428188875,455452683,        0,435168203,437265419,428877344,
     $326084382,330180441,327919318,319464469,454043295,326051612,
     $327984855,323692053,428188875,        0,430974230,293974816,
     $309015328,326117146,324023116,323367691,325429009,323321856,
     $443557067,445654283,430973722,294659808,325920416,436577739,
     $        0,428712733,296723360,303047775,307143897,308654877,
     $298820639,307148507,326018719,321922528,315598173,311207179,
     $460236383,317695325,436577739,        0,445654283,447751499,
     $441295834,298623831,296362898,300459152,317204113,325658388,
     $327919321,323823067,307082395,302851033,298558356,300491793,
     $306722256,321431186,325723863,323790426,317568096,319829067,
     $319127552,430974603,433071819,460333899,426779744,451946336,
     $426091595,451258187,        0,447751499,449848715,428647258,
     $300721173,304718994,310948698,298623957,302621778,310945233,
     $323561171,327853913,332215761,321463955,325756697,332212185,
     $441460320,440772171,        0,430384011,306553871,298427222/
      data (symbcd(j), j =    1255,  1368)/
     $296559517,303015136,317728415,328116058,329983763,323462667,
     $327526222,436708306,298525594,300852319,309343712,321890013,
     $328017686,325658255,432415820,455485196,        0,434873302,
     $298525591,300688473,313304536,319530581,321332876,325432855,
     $319235660,325429003,453682644,304718738,296231758,298198091,
     $310748556,319239251,300491664,298263500,304447488,435168203,
     $437265419,436937880,311207321,321660630,327788305,325527116,
     $314942731,306586638,449619480,323692243,325625486,319169931,
     $428876832,        0,455812629,321529493,323692056,315401433,
     $302785430,296330065,298263564,308651339,319170190,443327576,
     $300622739,298361806,304489675,        0,456140363,458237579,
     $455812568,313304281,302785430,296330065,298263564,308651339,
     $317072974,443327576,300622739,298361806,304489675,449848992,
     $455452491,        0,432645779,323659351,319563161,309109784,
     $298525523,296264590,302392523,312845836,323434067,321594904,
     $443327576,300622739,298361806,304489675,        0,445621470,
     $311338334,313500960,307242015,300852171,441459807,302949387,
     $428647705,428188875,        0,441230360,300655509,298427345,
     $302523535,310879632,317236755,319464919,315368729,307016728/
      data (symbcd(j), j =    1369,  1482)/
     $300622802,302527888,317269462,315373015,319563417,323757592,
     $434676624,296166221,298165322,314910281,323236685,298198091,
     $314943050,323233415,321037700,302129989,293839624,296035339,
     $        0,435168203,437265419,436937880,313304537,323757782,
     $325432793,321660566,323334944,303051531,308655563,331710464,
     $435168159,300885023,300954585,300266521,302363417,302822155,
     $308641792,437265375,302982239,303051865,304325637,297935620,
     $291676870,293839686,293778457,302228421,297939801,304906240,
     $435168203,437265419,458007567,447325899,445228683,428876832,
     $451716953,428188875,451258187,        0,435168203,437265419,
     $428876832,428188875,        0,434938827,437036043,436937880,
     $313304537,323757782,325432793,321660566,323335894,330049561,
     $340568408,348858763,474786072,346761547,428647449,428188875,
     $451258251,474327627,        0,434938827,437036043,436937880,
     $313304537,323757782,325432793,321660566,323334937,302822155,
     $308655563,331710464,443327512,298525523,296264590,302392523,
     $312845836,323430097,325691030,319563097,309114073,304882646,
     $298427281,300360780,308655435,317072974,323528339,321594840,
     $313294848,434938820,437036036,436937880,311207321,321660630/
      data (symbcd(j), j =    1483,  1596)/
     $327788305,325527116,314942731,306586638,449619480,323692243,
     $325625486,319169931,428647449,427959492,        0,455910980,
     $458008196,455812568,313304281,302785430,296330065,298263564,
     $308651339,317072974,443327576,300622739,298361806,304489675,
     $448931652,        0,434938827,437036043,436839510,309077337,
     $319596120,321627670,317433368,428647449,428188875,        0,
     $451651097,319464919,315368729,302818200,296461141,298460179,
     $313042384,319271766,298492948,313075153,319301133,317072715,
     $304456652,298230607,296067981,        0,435168207,302392459,
     $310748556,317142048,302490700,306557721,311197696,434938830,
     $302392523,312845836,323433497,302457932,308655769,323335897,
     $325432089,302822873,325891723,331710464,430744779,432841933,
     $455910603,426550361,447522521,        0,432841867,434939022,
     $449619083,449619595,451716750,466396811,426550425,460105817,
     $        0,432842315,434939531,458007435,428647577,449619737,
     $428188811,449160971,        0,432841995,434939149,458007819,
     $306422789,297935684,293774150,297972505,307017113,327974912,
     $453813067,455910283,432841557,296527449,430286411,321365515,
     $        0,445424728,300622740,296264526,298198091,308651340/
      data (symbcd(j), j =    1597,  1710)/
     $319268498,327886681,445424792,302719956,298361742,300295243,
     $445425049,319563350,325527308,329627033,317466134,323430092,
     $329623435,        0,451945759,307143705,300622738,296100612,
     $451945823,309240921,302719954,298197828,451946080,326084382,
     $328050393,323757527,309048928,326051547,323790424,317437143,
     $317400660,323561103,321299980,312845515,304489485,300430551,
     $315303444,321463887,319202764,312836096,426451800,300721241,
     $309077271,313140560,310780996,428581784,306980119,462202582,
     $323626317,306455556,460105366,321529165,        0,451683673,
     $309109784,298492754,296199053,300295243,308651404,319268434,
     $321562135,311305438,309339425,315663904,323957977,304882645,
     $298394510,300299467,312878543,319366678,317465947,311338271,
     $313533920,323944448,455812568,313304153,300688342,304751891,
     $439133208,302720148,311014675,300491600,296166284,304456971,
     $314975758,445228050,298328974,300295243,        0,447751391,
     $307176605,309208475,325953244,319661337,304849812,296264527,
     $298230859,310682951,312648964,306324549,449651863,300557201,
     $298296269,304447488,426418967,298624089,306979990,304686027,
     $437036120,304817170,298169426,309011800,317498969,325854999/
      data (symbcd(j), j =    1711,  1824)/
     $327821007,318912089,325822164,323462596,        0,426418967,
     $298624089,306979990,304653390,306586827,437036120,304817169,
     $302457932,308651339,317072974,325625620,330082141,328181408,
     $319825310,315499993,321595092,331953612,321365649,325723929,
     $328115935,324009984,437035922,296166220,298165323,308716815,
     $439133138,298263436,300253184,437035787,439133003,458008280,
     $327952089,321693144,308946003,300528723,308880716,314946643,
     $306783500,312845771,321267407,        0,430973920,305112222,
     $309208654,323364555,435168350,307111438,321267403,327529753,
     $293975321,296058880,439132868,441230084,439034896,302425227,
     $310748556,319235729,462202446,321267339,329623501,336050009,
     $323430028,325419008,437035915,439133203,300360587,460105365,
     $319338265,325789332,319333775,308716620,298169177,304906240,
     $447751391,307176605,309208475,321762715,307045401,300655573,
     $304719122,317273499,309142617,302752789,306816274,445195281,
     $298328910,296100810,310650183,312648900,304231698,304653264,
     $298263436,302327048,        0,443327512,298492754,296199053,
     $300295243,308651404,319268434,321562135,317465945,309114073,
     $304882645,298394510,300299467,312878543,319366678,317456384/
      data (symbcd(j), j =    1825,  1938)/
     $443294667,443294731,455878219,455878283,428549016,304916377,
     $428549015,304883608,        0,432546765,302392459,310748620,
     $321365650,323659351,319563161,311207000,300589970,289551627,
     $314975759,321463894,319567129,306979861,300491460,        0,
     $464299225,302785429,296297295,298230732,304456907,314975759,
     $321463893,319530456,313308377,304882645,298394510,300299467,
     $312878543,319366678,317470168,330039296,447489163,447489227,
     $428549016,304916249,428549015,304883480,        0,426418967,
     $298624089,306979990,302523405,306557977,304882774,300426189,
     $302392459,308651404,319235729,325723863,323790424,323725012,
     $457746135,        0,441197591,298492754,296199053,300295243,
     $310748620,323430161,329918295,325887577,317433171,308749316,
     $430416845,304489740,317105807,327726935,325854808,317400403,
     $308716612,        0,428647321,302785622,314811845,318911385,
     $300688406,312714629,318908036,460105367,319431561,293806788,
     $        0,456139972,458237060,426418967,298624089,306979990,
     $304653390,308684172,319203024,329888793,304882774,302556174,
     $304489675,314942988,323430161,329885657,        0,432710679,
     $309077145,302785429,296297295,298197963,304456908,312976786/
      data (symbcd(j), j =    1939,  2052)/
     $430416781,300295244,308716879,447292751,314975691,321234636,
     $329754514,332048216,327984856,330016661,447194509,317072972,
     $325494607,        0,451945099,451945995,449783243,432580049,
     $419799947,444966539,        0,443556683,445653899,437266144,
     $332376029,334342040,330016406,460334943,332310427,330049303,
     $323695702,323692309,329885521,327624332,314942091,457909973,
     $327788305,325527116,314933248,462366558,332408666,330180382,
     $326084192,315630815,305046490,298558291,296231821,300295307,
     $312845772,321332880,449848607,307143706,300655507,298329037,
     $302392459,        0,443556683,445653899,437266016,328181598,
     $332244887,329885391,321299916,308650635,456140511,328148827,
     $330016531,323462669,314975435,        0,443556683,445653899,
     $453846418,437266400,332212128,439035350,423994955,325592587,
     $        0,443556683,445653899,453846418,437266400,332212128,
     $439035350,423994443,        0,462366558,332408666,330180382,
     $326084192,315630815,305046490,298558291,296231821,300295307,
     $310748620,321332946,449848607,307143706,300655507,298329037,
     $302392459,444966284,319235730,451487634,        0,443556683,
     $445653899,470820491,472917707,437265888,464529696,439035734/
      data (symbcd(j), j =    2053,  2166)/
     $423994443,451258251,        0,443556683,445653899,437265888,
     $423994443,        0,456140047,308716684,302359435,294003406,
     $292037393,296231695,454042831,306619403,447751968,        0,
     $443556683,445653899,472917011,451651275,449554059,437265888,
     $464529632,423994443,451258187,        0,443556683,445653899,
     $437265888,423994955,325625355,        0,443556683,443557131,
     $445654349,472917259,472917707,475014923,437265696,472918368,
     $423994379,453355467,        0,443556683,443557518,443459211,
     $470820491,437265632,464529632,423994379,        0,449848543,
     $305046490,298558291,296231821,300295243,310748620,321332945,
     $327821144,330147614,326084192,315635104,311403677,302851031,
     $298427280,300328011,444966284,319235729,325723928,328050398,
     $321912832,443556683,445653899,437266208,334473245,336439256,
     $329983573,304789280,332376029,334342040,327886421,423994443,
     $        0,449848543,305046490,298558291,296231821,300295243,
     $310748620,321332945,327821144,330147614,326084192,315635104,
     $311403677,302851031,298427280,300328011,444966284,319235729,
     $325723928,328050398,321926093,300360720,306750673,313009550,
     $314811846,321070728,323270030,316941831,321103496,        0/
      data (symbcd(j), j =    2167,  2280)/
     $443556683,445653899,437266144,332376029,334342040,330016406,
     $304821984,330278813,332244824,327919254,449521173,321529484,
     $325429067,331786126,455747277,327558988,331788939,304447488,
     $464463774,334505882,332277598,328181344,313533599,302949403,
     $304915608,321529554,437101721,321562260,325658319,323397196,
     $314942603,300295053,296198993,293970765,298221568,451945547,
     $454042763,439362458,303048672,332212128,432383307,        0,
     $441459669,298361742,300295307,314943052,325527313,336606432,
     $302687185,300360716,306557920,315635552,342884352,437265483,
     $439362701,466625611,433071392,458237984,        0,441459723,
     $443556941,458236939,458237451,460334669,475014667,435168672,
     $468724064,        0,439363083,441460299,468722379,435168608,
     $460335200,421897163,447063755,        0,437265686,304460896,
     $313205899,468723030,433071392,460335200,432383307,        0,
     $466625227,468722443,441459674,305145824,426092107,325625355,
     $        0,466527124,331710464,432973716,298156032,455747095,
     $317465945,309109784,298492754,296199053,300295243,308651404,
     $319235665,323692187,321857055,315630816,305112094,302949469,
     $305083609,304882645,298394510,300299467,312878542,319333974/
      data (symbcd(j), j =    2281,  2394)/
     $321758750,315621376,428877067,430974221,462431499,428877600,
     $430941919,        0,453780889,309109784,298525523,296231821,
     $300295307,312845772,443327576,300622739,298329037,302392459,
     $432612754,        0,466625433,331953040,331887499,331710464,
     $433072025,298398608,331887499,331710464,468166479,325592658,
     $315303255,309077080,300655509,298427345,304620752,313042322,
     $321595096,330082265,        0,468821922,334538786,336701412,
     $330442467,321955359,317597080,310781128,306394786,321922588,
     $315106636,310682823,304260036,295838469,293806919,298001221,
     $        0,468821922,334538786,336701412,330442467,321955359,
     $317597080,310781128,306394786,321922588,315106636,310682823,
     $304260036,295838469,293806919,298001221,447587482,302785493,
     $300524560,306652493,317105806,327690067,329951000,323823067,
     $313360384,470394833,329787088,321431058,313206039,306979864,
     $298558293,296330129,302523536,310945106,319497815,325855064,
     $334211093,336166912,449717643,432678804,432383883,        0,
     $449717643,432940956,432678804,        0,432908045,462267277,
     $        0,451847580,317564444,317633428,336213453,314975691,
     $319169997,        0,439493700,441590916,479340804,481438020/
      data (symbcd(j), j =    2395,  2508)/
     $431106660,430056836,469903940,        0,434807700,300524564,
     $300580864,430744665,317109273,317044772,317030400,435299926,
     $297939876,319501156,319468388,345123229,343028677,344109956,
     $344074635,341966848,447751327,302916570,298558290,296166284,
     $302359691,312878543,319333972,323790493,321889760,313537888,
     $309306460,302851031,298394510,300295179,440771852,315074001,
     $319432281,321824287,317731798,319488000,443688035,303113184,
     $300885020,304981145,306947093,439460897,303015005,307111130,
     $309077142,298460306,308815054,306586699,302294023,304264211,
     $306750607,304522252,300229576,302195781,308412416,435299427,
     $307307744,309273756,304981017,302752917,439461025,307209309,
     $302916570,300688406,311043090,300426190,302392395,306488455,
     $304264339,302556175,304522380,308618440,306390085,300023808,
     $462169818,321758619,311239897,306914451,308847952,319301265,
     $325694875,311207126,308913425,313014043,325691089,329787344,
     $338241685,340502618,336471966,328181344,315630815,305079260,
     $298656599,296362897,300393549,308684171,321234700,331786190,
     $464365331,327722832,        0,426321109,325661394,309012178,
     $        0,298394766,308651209,306390020,300032901,295936842/
      data (symbcd(j), j =    2509,  2622)/
     $298263570,306881880,317498969,327952214,329852686,323364363,
     $317040012,315041231,319235533,455911128,327886610,325527180,
     $        0,458008082,317138380,319137483,329688975,460105298,
     $319235596,321238546,319464920,313304281,302785429,296297295,
     $298230732,304456907,312878543,319370457,304882645,298394510,
     $300285952,441459603,298329037,302396640,300528595,302720152,
     $311207321,319563351,323659410,321365452,310748299,302392271,
     $300529176,321594962,319268236,310752224,309329920,453715477,
     $321562198,319563161,309109784,298492754,296199053,300295243,
     $308651404,319272153,304882645,298394510,300285952,462431762,
     $317138380,319137483,329688975,464528978,319235596,321238546,
     $319464920,313304281,302785429,296297295,298230732,304456907,
     $312878543,319370457,304882645,298394510,300299872,330301440,
     $432546961,313075220,321594904,315401433,302785429,296297295,
     $298230732,304456907,314975758,443327576,300589970,298263500,
     $        0,456107550,321824414,323987040,317728095,311370972,
     $307012555,298033989,451945822,311305432,304587787,300163974,
     $295871172,287449605,285418055,289612357,432842265,        0,
     $460105163,314844421,304227204,293774022,291742472,295936774/
      data (symbcd(j), j =    2623,  2736)/
     $458007947,312747205,304231954,319464920,313304281,302785429,
     $296297295,298230732,304456907,312878543,319370457,304882645,
     $298394510,300285952,441459467,443556683,434709590,309077337,
     $317498968,323724949,319268364,321238489,321627733,317171148,
     $319137483,329688975,435168480,        0,443557023,309273887,
     $309342933,294364057,304915608,306881551,302392395,437036120,
     $304784335,300295179,308651341,315064320,445654239,311371103,
     $311440149,296461273,307012824,308978699,300163974,295871172,
     $287449605,285418055,289612357,439133336,306881483,298066758,
     $291635200,441459467,443556683,457975383,323692247,325854873,
     $321693144,308946003,300528723,308880716,314946643,306783500,
     $312845771,321267407,435168480,        0,441459602,296166220,
     $298165323,308716815,443556818,298263436,300266464,309329920,
     $426418967,298624089,306979990,304686027,437036120,304817170,
     $298169426,309011800,317498969,325854999,327853643,455911127,
     $325756427,459876182,334243929,342665560,348891541,344434956,
     $346405081,346794325,342337740,344304075,354855567,        0,
     $426418967,298624089,306979990,304686027,437036120,304817170,
     $298169426,309011800,317498969,325854999,327853711,323364555/
      data (symbcd(j), j =    2737,  2850)/
     $455911127,325756495,321267339,329623501,336035840,443327512,
     $298492754,296199053,300295243,308651404,319268434,321562135,
     $317465945,309114073,304882645,298394510,300299467,312878543,
     $319366678,317456384,426418967,298624089,306979990,304685892,
     $437036120,304817170,293745746,306881816,315401753,323757783,
     $327853842,325559884,314942731,306586703,304690840,325789394,
     $323462668,314946116,302120960,458007812,460105028,453584405,
     $317465945,309109784,298492754,296199053,300295243,308651340,
     $317171218,443327576,300589970,298263500,438445572,        0,
     $426418967,298624089,306979990,304686027,437036120,304817170,
     $298169426,309011800,317498969,323757719,321594903,321650688,
     $453748246,321594967,319563097,307012568,298558357,300557712,
     $317174678,300590481,317203917,314975435,302359372,294036238,
     $296166221,        0,443556818,298263436,300262539,310814031,
     $445654034,300360652,302363481,315392000,426418967,298624089,
     $306979989,302490637,306557977,304882773,300393421,302392459,
     $310748556,319235730,462202514,321332812,323331915,333883407,
     $464299730,323430028,325419008,426418967,298624089,306979989,
     $302490637,306557977,304882773,300393421,302392459,308651404/
      data (symbcd(j), j =    2851,  2964)/
     $319235729,325756633,323790551,        0,426418967,298624089,
     $306979989,302490637,306557977,304882773,300393421,302392459,
     $310748556,319235664,460105296,321300108,327526283,335947918,
     $342370580,344762585,344700697,323495565,327516160,430613464,
     $304915737,313238868,443327767,311043280,306652172,298165067,
     $294003469,296166285,296105168,308716811,317040204,325564120,
     $323725014,327919384,325887641,319563158,313140496,310814027,
     $        0,426418967,298624089,306979989,302490637,306557977,
     $304882773,300393421,302392459,310748556,319235730,464299595,
     $319038853,308421636,297968454,295936904,300131206,462202379,
     $316941637,308412416,460105367,319464463,298230603,432710615,
     $304915737,319534039,304882968,319530647,432448525,310781388,
     $321303565,310748619,321300111,        0,433202052,435299268,
     $433202532,432153924,        0,443688132,445785348,431105316,
     $430056708,        0,447751044,460334340,432711445,430417615,
     $        0,447653148,313370012,315532639,309339232,300917661,
     $298689497,304850324,434939158,315237842,317203854,310785048,
     $298525524,296297360,302458187,432547021,312845705,314811717,
     $308421700,300065671,298066889,302261191,        0,441459806/
      data (symbcd(j), j =    2965,  3078)/
     $307111134,307246240,306328725,304686212,308880533,428647320,
     $302818202,294433561,319599897,315368985,315434265,        0,
     $434938776,300655640,300725197,298197963,302392269,        0,
     $434938776,300655640,300725195,298197965,302392330,300163975,
     $        0,435168158,300491806,300954590,300692429,298197963,
     $302392269,        0,432939995,298656603,296625054,300917856,
     $311436767,319759964,321725976,317433045,308884768,315598302,
     $319694362,317465942,442934412,308651276,308707328,468722507,
     $441459998,311305434,304915417,296592221,298820640,307242271,
     $317662878,330278880,459875921,319268365,323331851,331753422,
     $333981522,325648384,468461463,334178327,336340953,332179288,
     $327886481,319235468,310748235,298197838,296264595,311141785,
     $317564381,315598112,307209309,304981144,311076430,325461899,
     $333817868,335983691,300295054,298361811,304788571,307013262,
     $327559051,        0,437035992,302752856,302822221,294003531,
     $298188800,437035992,302752856,302822219,294003533,298197899,
     $296002247,        0,441459807,300528799,300528800,309306323,
     $430351116,296067980,296124416,439231643,304948251,302916702,
     $307209568,321922847,330213211,327984856,313205973,308913426/
      data (symbcd(j), j =    3079,  3192)/
     $315176544,326084381,328050393,323757591,440837196,306554060,
     $306610176,430482259,298525719,306947350,319399570,327755667,
     $334148435,298492950,306914581,319366801,327722898,334145495,
     $        0,445784916,310509568,433202516,297926656,433202052,
     $        0,435168153,437265305,451945881,454043033,        0,
     $323397323,441131922,296231758,298197835,430449612,432612240,
     $300360652,296072531,323761693,319628888,325854938,321758749,
     $453944922,325844992,437265311,296657755,298624024,306980121,
     $313369949,311403680,303038464,464201748,329856665,334112399,
     $432678868,        0,454042756,456139844,445424664,298525523,
     $296231822,302392523,314943116,327624529,329918230,323757529,
     $311211289,304882646,298427280,300360780,308655499,321267406,
     $327722772,325789272,317489152,443557017,445654169,        0,
     $306787478,304751824,306652240,308946070,441001092,440673350,
     $306324678,306459417,298591257,298656537,428647961,445425048,
     $319595930,311210763,298132491,298197771,428189195,444966282,
     $319137164,310738944,443556895,298722135,296362895,302392523,
     $312845836,323462868,325822108,319792480,309329920,437134493,
     $313533771,        0,432907164,300885023,307242400,319792734/
      data (symbcd(j), j =    3193,  3306)/
     $323888794,321660373,296068811,        0,435168928,311174616,
     $321627798,325691089,323429900,312845451,300295053,296189952,
     $451945298,327759328,317030400,456139744,298558424,307012953,
     $319563414,325691089,323429900,312845451,300295053,296189952,
     $458139231,315630880,305112028,298558354,300360780,310748491,
     $319170190,325625554,323659287,313271576,304849877,298385408,
     $460334155,430974688,        0,441459679,298754971,300721240,
     $313239062,323626706,325559949,321267083,306553804,298230607,
     $296297364,302720215,317466201,323856029,321889696,307232768,
     $458008150,317334803,308913172,298525529,296559517,303015136,
     $311436767,321824409,323626575,317072651,306553804,298254336,
     $451847627,432678932,        0,432678932,        0,466756356,
     $        0,432777239,432580625,        0,447882466,305112027,
     $298525586,300328009,308487492,        0,431104994,305112283,
     $311108882,308716617,300098372,        0,441263246,430679505,
     $451650385,        0,436609995,298197965,302392330,300163975,
     $        0,434545548,300262412,300318720,441590919,449979783,
     $460236383,315630752,300917597,296592281,300688471,317367892,
     $323593937,325527116,314942603,300294990,        0,443556895/
      data (symbcd(j), j =    3307,  3420)/
     $298722135,296362895,302392523,312845836,323462868,325822108,
     $319792480,309343456,305112094,300819351,298460111,302425164,
     $308655435,317072909,321365652,323724892,319759839,313524224,
     $437134493,313533771,445621515,436577867,        0,432939995,
     $298656603,296625054,300917920,315631199,323954396,325920408,
     $317400212,302621585,296166219,449848863,321857180,323823192,
     $315303060,430351246,302458188,319170189,325530638,312845899,
     $323364558,325582848,432939995,298656603,296625054,300917920,
     $315631199,323921562,321660311,309048736,319792733,321725976,
     $315340183,319497876,325658319,323397196,314942603,300295053,
     $296198992,298361808,298301013,323561103,321299980,314933248,
     $449783179,451945931,451945233,327726283,323321856,435168086,
     $430646232,307012953,319563414,325691089,323429900,312845451,
     $300295053,296198992,298361808,298300761,317466198,323593873,
     $321332684,312849376,321926111,311404128,        0,456042012,
     $321758876,323921503,317728032,305112029,298689367,296264590,
     $302392523,312845836,323430097,325658261,319530328,311174231,
     $300589970,445654175,302949339,298558353,300360780,308655435,
     $317072974,323528338,321562071,313262080,430973786,430842782/
      data (symbcd(j), j =    3421,  3534)/
     $303047840,317630045,323954400,433005599,307209693,460334813,
     $323822997,313107728,310752922,313173267,308815051,        0,
     $441459679,298754970,300688535,315336280,323823261,321889696,
     $307246240,303014877,300753944,306951575,319563354,321824287,
     $315634839,300622741,296330063,298230732,306554251,321267341,
     $325560019,323659350,315339927,302719957,298427279,300327948,
     $306558347,319170125,323462803,321562134,315326464,458008150,
     $317334803,308913172,298525529,296559517,303015136,313533983,
     $323921626,325723792,321332684,310748235,300295054,298296272,
     $302490574,443130964,300622745,298656733,305112288,447751647,
     $321824410,323626576,319235468,310738944,451847627,432678932,
     $        0,432678932,        0,466756356,        0,432777239,
     $432580625,        0,447882466,305112027,298525586,300328009,
     $308487492,443622494,302883798,300491789,304424134,        0,
     $431104994,305112283,311108882,308716617,300098372,435233886,
     $307078358,308880525,304423878,        0,441459860,430876119,
     $451846999,        0,434480012,300327948,302326728,298024960,
     $434545548,300262412,300318720,441590919,449979783,458139228,
     $323856092,326018655,315630752,300917597,296592281,300688471/
      data (symbcd(j), j =    3535,  3648)/
     $317367892,325661531,300721240,317400661,323626706,325527116,
     $314942603,300294990,296199056,300393358,        0,449848543,
     $305046490,298558291,296231821,300295243,308651404,319235729,
     $325723928,328050398,323986976,315635104,311403677,302851031,
     $298427280,300328011,442869068,317138513,323626712,325953182,
     $319815680,449717323,454042763,454042973,307078170,451847387,
     $302841856,439231643,304948251,302916702,307209568,319825631,
     $328115995,325887575,315270291,300458831,291878432,323987165,
     $325953177,319530131,428254030,300360972,317072973,323466190,
     $310748619,321267343,        0,439231643,304948251,302916702,
     $307209568,319825631,328115995,325887511,313210400,323987165,
     $325953177,319534294,313206293,321529490,323462733,319169867,
     $304456588,296133391,294134609,298328911,447423957,319432274,
     $321365517,317072715,        0,458204427,460334411,460333841,
     $327712768,443556758,443557728,443524639,330314646,300655768,
     $313271831,321595028,323528270,317072651,304456588,296133391,
     $294134609,298328911,447489495,319497812,321431054,314975499,
     $        0,460236444,325953308,328115935,321922464,309306461,
     $300753815,296330063,298230732,304456971,317072974,323495571/
      data (symbcd(j), j =    3649,  3762)/
     $321562134,315335895,304817108,298399136,311403677,302851031,
     $298427278,300299531,314975758,321398356,319488000,437265306,
     $464529181,323822932,308847759,304461466,311043217,304587787,
     $435070112,311436893,437200031,311404125,326018846,330301440,
     $447751327,305079324,302818391,309011862,323725016,328017693,
     $326084128,313537888,309306526,305013849,306947286,449521239,
     $323757786,326018719,319829206,300589907,294167310,296100875,
     $310748684,321300111,323561044,319464854,443229205,298427217,
     $296166284,302363915,317072909,321365587,319455232,460105367,
     $319464852,308946005,302719960,300786717,307209568,319825567,
     $326051612,327952084,323528206,314975435,302359436,296166223,
     $298329039,298267733,302752795,305046751,313538207,326018776,
     $323626577,317138252,308641792,451847627,432678932,        0,
     $432678932,        0,475144708,        0,432777239,432580625,
     $        0,456271201,307176475,298558290,296166281,300098564,
     $447784093,302818262,298361740,300131332,        0,443688226,
     $313501082,315303249,308716618,298033796,443688225,313402711,
     $310977743,304456583,        0,445654292,435070551,456041431,
     $        0,430285580,296133516,298165065,291733504,430351116/
      data (symbcd(j), j =    3763,  3876)/
     $296067980,296124416,449979271,460465351,462300891,328017755,
     $330180382,326084128,311436383,300852187,302818392,319432338,
     $435004505,319465044,323561103,321299980,312845387,298197837,
     $294101776,296264592,296189952,443556895,298722135,296362895,
     $302392523,312845836,323462868,325822108,319792480,309343327,
     $300819351,298460111,304493581,308684108,319206860,321365652,
     $323724892,317699614,313500895,302972928,437134493,313533771,
     $437134363,307111198,310748491,        0,432907164,300885023,
     $307242400,319792734,323888794,321660373,298169243,300786652,
     $302982303,315598366,321791578,319563157,296072076,325461707,
     $430286539,        0,435168928,309048288,300918367,456139927,
     $443295064,319530645,325658321,323429900,312845451,300295053,
     $296199055,441165143,319497875,449554005,323561105,321332620,
     $457713165,312878220,300327823,438707086,        0,451847627,
     $319141408,319141408,296232720,451847056,432580369,327680000,
     $435168151,437232600,435168864,321893407,321893336,307012953,
     $319563414,325691089,323429900,312845451,300295053,296199055,
     $432776151,304883032,319530644,449586774,323593873,321332620,
     $457713165,312878220,300327823,438707086,        0,454010461/
      data (symbcd(j), j =    3877,  3990)/
     $323921503,315630880,305112028,298558354,300360780,310748491,
     $319170190,325625554,323659287,313271576,304849877,456074655,
     $311403614,441426972,300655570,302458060,434644045,310781260,
     $319202960,449193550,323528338,321562007,457811478,313238807,
     $304817107,443261973,300482560,430974688,304460640,296724127,
     $458236939,304447488,441459679,298754971,300721176,306947478,
     $319465044,323561103,321299852,306586573,298296210,300557333,
     $306914711,319563353,323856029,321889696,307246111,300852187,
     $302818456,315336214,323626706,325559949,321267083,306553804,
     $298230607,296297364,302720151,315368985,321758813,319796830,
     $315597983,300888974,304494028,323420160,455812564,311010515,
     $302654358,296526682,298755103,309339424,317695581,323790484,
     $321365452,310748299,300295054,300360716,455910934,313144920,
     $317367572,308945941,298595476,300622745,298656733,307213211,
     $302982367,311403998,321762655,319727193,321529359,314979789,
     $310781068,300318720,449750412,317076893,317629900,432711637,
     $334115733,298461140,        0,432711637,334115733,298461140,
     $        0,466756356,295843748,334635844,        0,432842713,
     $334246809,298592216,432580561,333984657,298330064,        0/
      data (symbcd(j), j =    3991,  4104)/
     $445785250,303014811,296428370,298230793,306390276,312620324,
     $313664738,305112027,298525586,300328009,308487492,        0,
     $431104994,305112283,311108882,308716617,300098372,297939812,
     $298984482,307209499,313206098,310813833,302195588,        0,
     $441459807,308978836,441459860,441459935,304784532,430875549,
     $315336151,430876119,430875484,317466071,451847581,298558295,
     $451846999,451847644,296493911,        0,438707211,300262284,
     $298230734,302457933,304423944,298038221,300295180,302425037,
     $436577354,438707208,        0,434578317,298197963,302359628,
     $304522254,300364749,300295180,302425037,        0,443688135,
     $310621412,311567623,453944989,319792480,307241951,296657755,
     $298623960,317335059,321431119,319202636,306586637,300365341,
     $317662559,307209182,298754971,300721621,321496721,323462733,
     $319169867,306553804,296166350,455550348,        0,445653771,
     $445555531,293975325,325429003,445654795,434677329,432547472,
     $        0,433070987,435135436,433071520,321889950,325986009,
     $323724886,315274207,315598430,323888793,321627542,434840982,
     $321562260,325658319,323397196,314942347,434808213,321529490,
     $323462733,314975180,        0,462268125,321889760,309339231/
      data (symbcd(j), j =    4105,  4218)/
     $300852123,296493907,298329038,304489675,317040204,325527312,
     $462268123,323921502,317695199,305079259,298591123,300426317,
     $308684236,321300110,325592848,        0,433070987,435135436,
     $433071456,319792797,325953304,327788240,323429900,312845195,
     $435135839,319759965,323856088,325691024,321332749,312878028,
     $        0,433070987,435135436,433071776,435136159,324023254,
     $313206101,434808149,434513548,323335051,323321856,433070987,
     $435135435,298169248,324023263,323987104,434840918,313177045,
     $313163776,462268125,321889760,309339231,300852123,296493907,
     $298329038,304489675,317040204,325527312,327820756,462268123,
     $323921502,317695199,305079325,300786584,298427344,302457933,
     $308684236,321300110,325592787,317302228,        0,433070987,
     $433071072,300262283,462431968,325429003,462432011,434841302,
     $434808533,        0,433070987,300266400,300950475,        0,
     $449848720,312911052,304489421,298328912,449848800,317203853,
     $312878283,304456652,298230608,        0,433070987,300266400,
     $300950475,462431968,300562208,300528791,325429003,443262731,
     $        0,433070987,433071072,300299212,323364491,432383627,
     $        0,433070987,435004363,298169307,314946464,315045792/
      data (symbcd(j), j =    4219,  4332)/
     $315045723,314947419,329623435,466626443,        0,433070987,
     $435069899,298169309,327529376,325531360,325531360,328214283,
     $        0,443556959,300852123,296493907,298329038,304489675,
     $317040204,325527312,329885528,328050397,321889760,309343519,
     $305079259,298591123,300426317,310781324,321300176,327788312,
     $325953118,315598111,        0,433070987,435135435,298169248,
     $317728351,323954396,325887639,321594837,300594143,317695582,
     $323888793,321627606,300613632,443556959,300852123,296493907,
     $298329038,304489675,317040204,325527312,329885528,328050397,
     $321889760,309343519,305079259,298591123,300426317,310781324,
     $321300176,327788312,325953118,315598111,449259209,327464334,
     $317138697,        0,433070987,435135435,298169248,315631199,
     $323954396,325887639,321594773,300594143,315598430,323888793,
     $321627542,300627221,323331787,447391435,        0,460236383,
     $315630752,300917597,296592281,300688471,315270676,321496721,
     $323429965,314975372,302425038,296171229,321824286,315597983,
     $300884893,298689497,304883094,319465107,325625550,321267083,
     $306553804,296157184,441427083,443524299,306557728,321922655,
     $428876575,321880064,433070993,300360780,310748555,321267406/
      data (symbcd(j), j =    4333,  4446)/
     $327722784,433071072,300459022,304522508,314975821,323430097,
     $326117152,        0,428877067,428876640,310851360,326116622,
     $462431499,        0,428876939,428876640,306656736,306656733,
     $306558429,327529952,327628960,338700046,475014923,        0,
     $430974603,325432160,298854091,460334752,296072928,298165067,
     $        0,428877014,308651275,428876640,311113440,324019414,
     $460334358,310738944,458236747,460333963,430974688,430973791,
     $323990412,325461707,430286539,        0,455910987,323335769,
     $323790475,455812568,313304217,302785430,296330065,298263564,
     $306554187,317072974,455812440,306979863,300622739,298361806,
     $302425228,312878670,        0,433070987,300266400,300950475,
     $434840664,309110169,319563414,325691089,323429900,314942667,
     $304489422,434840792,315368983,321595027,323528270,319202700,
     $308683726,        0,455812568,313304217,302785430,296330065,
     $298263564,306554187,317072974,455812629,317433176,306979863,
     $300622739,298361806,302425228,312878541,319268430,        0,
     $456140363,323335776,324019851,455812568,313304217,302785430,
     $296330065,298263564,306554187,317072974,455812440,306979863,
     $300622739,298361806,302425228,312878670,        0,432612946/
      data (symbcd(j), j =    4447,  4560)/
     $321562135,317465945,307012632,298525523,296264590,302392459,
     $312845772,321336211,319399445,317433176,306979863,300622739,
     $298361806,302425228,312878541,319268430,        0,447751392,
     $305112092,302359627,447751519,309306462,441427036,304460633,
     $311207192,430744408,311164928,458008153,321201671,316876101,
     $308454470,302228359,458008202,321103301,312616068,302162823,
     $455812568,313304217,302785430,296330065,298263564,306554187,
     $317072974,455812440,306979863,300622739,298361806,302425228,
     $312878670,        0,433070987,300266400,300950475,434807960,
     $311207385,321660565,323335125,306947352,315368983,321562187,
     $323321856,433070943,296690589,300852254,303014880,298857375,
     $298787806,300917663,432841611,300266393,300721099,        0,
     $433070943,296690589,300852254,303014880,298857375,298787806,
     $300917663,432841604,300037017,300721092,        0,433070987,
     $300266400,300950475,458008153,300398233,300364946,319137419,
     $443131531,        0,433070987,300266400,300950475,        0,
     $432841611,300266393,300721099,434807960,311207385,321660565,
     $323335125,306947352,315368983,321562187,323335829,330049497,
     $340568344,346728779,457877335,334243928,342599957,344303947/
      data (symbcd(j), j =    4561,  4674)/
     $        0,432841611,300266393,300721099,434807960,311207385,
     $321660565,323335125,306947352,315368983,321562187,323321856,
     $441230360,298525523,296264590,302392459,312845772,321332881,
     $323593814,317465945,307016856,302752726,298427281,300360717,
     $306586956,317105678,321431123,319497687,313271448,        0,
     $432841604,300037017,300721092,434840664,309110169,319563414,
     $325691089,323429900,314942667,304489422,434840792,315368983,
     $321595027,323528270,319202700,308683726,        0,455910980,
     $323106393,323790468,455812568,313304217,302785430,296330065,
     $298263564,306554187,317072974,455812440,306979863,300622739,
     $298361806,302425228,312878670,        0,432841611,300266393,
     $300721099,434742294,306980121,317502419,302687383,311174616,
     $317489152,453715416,311207001,298591062,298460179,313042384,
     $449357263,317138316,451323148,304489357,434512782,296171030,
     $317400472,451650840,304882583,434906006,300561301,302654802,
     $317236751,319235532,310748235,298197838,        0,435168203,
     $302363616,303047691,428647641,309080857,294397144,        0,
     $432841615,300295243,310748556,321368985,300721103,302425228,
     $310781325,321369689,321234571,455911065,323321856,428647563/
      data (symbcd(j), j =    4675,  4711)/
     $428647257,306624025,317498509,453813387,        0,430744715,
     $430744473,306656665,306656662,306558358,323335577,323434457,
     $332179086,468493963,        0,430745099,321237849,298624587,
     $455910937,296072793,298165067,        0,428647563,428647257,
     $306624025,317498509,297940505,306553796,297926656,451683147,
     $455910348,430745177,430744408,317469644,321267275,430286411,
     $        0/
c
      data (istart(j), j=1,229)/
     $   1,   5,  16,  26,  34,  39,  43,  54,  58,  60,  66,  70,  73,
     $78,  82,  93, 100, 112, 120, 131, 134, 140, 143, 148, 151, 154,
     $ 296, 305, 314, 322, 331, 340, 344, 355, 360, 364, 370, 374, 376,
     $385, 390, 399, 408, 417, 421, 430, 434, 439, 442, 447, 450, 455,
     $3177,3186,3189,3197,3205,3208,3217,3229,3232,3247,3259,3262,3264,
     $3266,3269,3275,3281,3285,3290,3293,
     $ 158, 162, 173, 176, 180, 185, 189, 193, 205, 207, 211, 214, 219,
     $223, 227, 238, 242, 249, 253, 256, 265, 275, 278, 287,
     $ 459, 471, 486, 494, 506, 515, 526, 535, 549, 554, 563, 567, 577,
     $584, 598, 607, 613, 623, 632, 636, 644, 655, 662, 672,
     $ 683, 690, 710, 726, 740, 749, 757, 775, 785, 790, 799, 809, 815,
     $826, 834, 855, 868, 898, 918, 935, 942, 952, 958, 967, 975, 983,
     $1272,1290,1305,1319,1335,1350,1360,1388,1399,1406,1417,1427,1432,
     $1450,1461,1478,1494,1509,1519,1535,1542,1553,1559,1568,1576,1585,
     $3306,3325,3330,3351,3373,3378,3396,3419,3433,3462,3485,3488,3490,
     $3492,3495,3505,3515,3519,3523,3526,
     $ 990, 997,1017,1023,1029,1038,1045,1055,1080,1085,1095,1101,1112,
     $1120,1133,1154,1162,1175,1183,1190,1205,1226,1234,1252,
     $1592,1611,1637,1650,1671,1686,1701,1716,1737,1744,1757,1767,1779/
      data (istart(j), j=230, 432)/
     $1789,1810,1825,1834,1849,1865,1872,1887,1905,1916,1932,
     $1953,1960,1978,1995,2009,2018,2026,2046,2056,2061,2071,2081,2087,
     $2098,2106,2126,2138,2167,2185,2202,2209,2220,2226,2235,2243,2251,
     $2522,2540,2556,2568,2587,2600,2617,2637,2651,2663,2678,2693,2701,
     $2725,2742,2757,2776,2791,2803,2817,2825,2842,2855,2874,2894,2913,
     $3546,3566,3572,3592,3616,3620,3638,3660,3673,3702,3724,3727,3729,
     $3731,3734,3744,3754,3758,3762,3765,
     $4074,4082,4102,4121,4136,4146,4154,4176,4185,4189,4199,4208,4214,
     $4224,4232,4252,4264,4287,4302,4323,4329,4341,4347,4357,4364,4371,
     $4379,4396,4413,4429,4446,4464,4474,4497,4508,4519,4530,4539,4543,
     $4562,4573,4591,4608,4625,4634,4656,4663,4674,4680,4690,4697,4704,
     $3784,3803,3809,3825,3846,3853,3876,3904,3909,3941,3969,3976,3980,
     $3984,3991,4003,4015,4031,4042,4050,
     $2258,2260,2262,2283,2289,2301,2305,2309,2320,2336,2360,2373,2377,
     $2381,2384,2391,2399,2402,2406,2415,2435,2454,2473,2500,
     $2927,2932,2937,2942,2964,2977,2983,2990,2997,3012,3027,3051,3056,
     $3063,3070,3086,3098,3100,3102,3104,3123,3130,3135,3154/
      data (width(j), j=1,216)/
     $18.,21.,21.,21.,19.,18.,21.,22., 8.,16.,21.,17.,24.,22.,22.,21.,
     $22.,21.,20.,16.,22.,18.,24.,20.,18.,20.,
     $19.,19.,18.,19.,18.,12.,19.,19., 8.,10.,17., 8.,30.,19.,19.,19.,
     $19.,13.,17.,12.,19.,16.,22.,17.,16.,17.,
     $20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,26.,26.,22.,26.,14.,14.,
     $16.,10.,10.,20.,
     $18.,21.,17.,18.,19.,20.,22.,22., 8.,21.,18.,24.,22.,18.,22.,22.,
     $21.,18.,16.,18.,20.,20.,22.,20.,
     $21.,19.,19.,18.,16.,15.,20.,21.,11.,18.,16.,21.,18.,16.,17.,22.,
     $18.,20.,20.,20.,22.,18.,23.,23.,
     $20.,22.,21.,22.,21.,20.,23.,24.,11.,15.,22.,18.,25.,23.,22.,22.,
     $22.,22.,20.,19.,24.,20.,24.,20.,21.,20.,
     $20.,21.,19.,21.,19.,13.,19.,22.,11.,11.,21.,11.,33.,22.,20.,21.,
     $20.,17.,17.,15.,22.,18.,24.,20.,19.,18.,
     $20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,26.,26.,22.,26.,14.,14.,
     $16.,10.,10.,20.,
     $20.,22.,18.,20.,21.,20.,24.,22.,11.,22.,20.,25.,23.,22.,22.,24.,
     $22.,21.,19.,19.,21.,20.,23.,22./
      data (width(j), j= 217, 432)/
     $23.,21.,20.,19.,18.,18.,22.,23.,12.,20.,20.,23.,20.,17.,18.,22.,
     $19.,21.,20.,20.,22.,18.,23.,23.,
     $20.,24.,21.,23.,23.,22.,22.,26.,13.,18.,23.,20.,27.,25.,22.,23.,
     $22.,24.,23.,21.,25.,20.,26.,22.,21.,22.,
     $21.,19.,18.,21.,18.,15.,20.,21.,13.,13.,20.,12.,33.,23.,18.,21.,
     $20.,17.,17.,14.,23.,20.,29.,20.,21.,20.,
     $21.,21.,21.,21.,21.,21.,21.,21.,21.,21.,26.,26.,22.,26.,15.,15.,
     $17.,11.,11.,21.,
     $20.,20.,21.,21.,19.,18.,21.,22., 9.,17.,21.,17.,24.,22.,22.,20.,
     $22.,20.,20.,17.,22.,20.,26.,20.,19.,20.,
     $20.,20.,18.,20.,18.,14.,20.,20., 9., 9.,19., 9.,31.,20.,19.,20.,
     $20.,14.,17.,11.,20.,16.,24.,18.,16.,18.,
     $20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,25.,25.,23.,25.,14.,14.,
     $16.,11.,11.,19.,
     $24.,24.,19.,20.,17.,24.,24.,25.,24.,24.,25.,24.,24.,22.,26.,34.,
     $10.,22.,31.,19.,14.,14.,27.,22.,
     $14.,14.,21.,16.,16.,10.,10.,10.,18.,24.,25.,11.,11.,11.,21.,24.,
     $14.,14., 8.,16.,14.,26.,22., 8./
      data (ssymbc(j), j=1,120)/
     $             471149226,357246358,315959338, 68157440,470825002,
     $ 345320100,357443862,327886236,315762474,336920576,470825002,
     $ 355313115,336920576,470493226,449850016,0,455911911,456370649,0,
     $ 471149216,336274848,336922656,0,470493226,357574048,336920576,
     $ 449522346,315959958,0,470825002,355641947,336274907,317892650,0,
     $ 456370208,336279584,351502336,481470811,325953253,347256234,
     $ 326284694,325958294,346929184,357892096,449850016,470493226,
     $ 455911911,485271143,0,450177706,315304598,315949056,470493226,0,
     $ 470825002,355313115,336935525,336274917,355631104,470853600,
     $             336570464,336625664,468592477,328181537,330409956,
     $ 338831587,345024799,342796380,334364672,466265814,319563163,
     $ 313468258,315794984,326444971,341158250,353643173,359738078,
     $ 357411352,346761365,332038144,465905227,312910991,300491605,
     $ 292332190,290530023,297116654,307799411,322611126,341518837,
     $ 360295345,372714731,380874146,382676313,376089682,365406925,
     $ 350595210,331677696,468592477,328181537,330409956,338831587,
     $ 345024799,342796380,334378847,330344289,466560930,468625379,
     $ 470722595,472819811,474949794,477079777,0,462300964,345123100,
     $ 328087389,330413981,332511197,334608413,336705629,338802845/
      data (ssymbc(j), j=121,128)/
     $             340900061,342982656,470623971,347187226,464594973,
     $                            342964256,334571552,338755584/
      data isstar /1,5,11,14,17,20,24,27,30,35,38,45,50,53,55,60,63,70,
     $ 81,98,113,123/
c
      end
c______________________________________________________________________
c***************************************************************
c POSTSCRIPT DRIVERS FOR THE STANDARD PLOT SUBROUTINES
c***************************************************************
      subroutine plot(x, y, ipen)
c$$$$ calls nothing
      common /post/ boxx,boxy,land,infill
      common /outsid/ box(4),xorig,yorig
      save moves
      data moves/0/
c  PostScript output generator for the old CalComp line-drawing routine
c  Writes to unit 11; abbreviations  k, l, m, n, f  must be defined
ce
      ix=1000.0*x
      iy=1000.0*y
c
c  Extend current polygon; if 1000 reached, stroke it and start over
      if (ipen .eq. 2) then
c  Keep record of maximum excursion for Encapsulation
        box(1)=min(box(1), x+xorig)
        box(2)=min(box(2), y+yorig)
        box(3)=max(box(3), x+xorig)
        box(4)=max(box(4), y+yorig)
        if (moves .le. 1000) then
          write(11,*) ix,iy,' l'
          moves=moves + 1
        else
          moves=0
          write(11,*) ix,iy,' n'
        endif
c  Stroke and terminate current ploygon; start a new one at (x,y)
      elseif (ipen .eq. 3) then
        if (moves .gt. 0) then
c  Fill the polygon if infill set to 1; otherwise just stroke
          if (infill .eq. 1) write(11,*) ix,iy,' f'
          if (infill .eq. 0) write(11,*) ix,iy,' k'
          moves=0
        else
          write(11,*) ix,iy,' m'
        endif
      elseif (ipen .lt.  0) then
        write(11,*) ix,iy,' translate',
     $  0 ,0,  ' moveto'
        xorig=xorig + 0.001*ix
        yorig=yorig + 0.001*iy
      elseif (ipen .eq. 0) then
c  If ipen = 0, eject a page;  insert signal for psview too
        write(11,'(2i8,a/(a))') ix,iy,' moveto',
     $  'currentpoint translate',
     $  'showpage',
     $  '%%Page: fresh page started ---------------------------'
c  If ipen = 4 write the fill command
      elseif (ipen .eq. 4) then 
         write (11,*) 'currentpoint f'
      else
        write(*,*)'Unintelligible pen command in plot ignored'
      endif
      return
      end
c_______________________________________________________________
      subroutine plopen(k, name)
c$$$$ calls nothing
c  When k > 0 opens unit 11 for PostScript and writes the short prolog.
c  Prolog includes landscape rotation if land > 0 in /caps/
c  When k < 0 closes the file
c
      common /post/ boxx,boxy,land,infill
      common /outsid/ box(4),xorig,yorig
      character*(*) name
c
c  When k < 0 closes the file: flush polygon buffer and display
c  info for encapsulating written after showpage
      if (k .lt. 0) then
         write(11,'(a)') 'currentpoint n','showpage'
         write(11,'(a/a,4i8)') '%!PS-Adobe-2.0 EPSF-1.2',
     $   '%%BoundingBox:',(nint(72*box(j)),j=1,4)
         close(unit=11)
         do 1100 j=1, 4
           box(j)=(2.5-j)*1.0e6
 1100   continue
        xorig=0.0
        yorig=0.0
c
c  When k > 0 opens unit 11 for PostScript; writes prolog
      else
         if (name .eq. 'myplot') name='mypost'
         open(unit=11, file=name)
         write(*,'(/a/2a/a)')
     $   ' --------------------------------------',
     $   ' PostScript file written to: ',name,
     $   ' ========== ---------------------------'
c
c  Write PostScript prolog containing PS procedures
         write(11,'(a)')'%!PS PostScript file',
     $   '/k {currentpoint stroke moveto newpath moveto} def',
     $   '/l {lineto} def',
     $   '/m {newpath moveto} def',
     $   '/n {lineto currentpoint stroke moveto} def',
     $   '/f {currentpoint closepath fill moveto newpath moveto} def',
     $   '0.072 0.072 scale % Coords are 1/1000 th inch'
         if (land .gt. 0) write(11,'(a)')
     $  '8500 0 translate 90 rotate    % anticlockwise',
     $  '%0 10750 translate -90 rotate % clockwise'
c  Set up a background color
         write(11,'(a)')
     $   '2 setlinejoin 6 setlinewidth 1 setlinecap',
     $   '% Uncomment next 2 lines for light blue background',
     $   '% 0.67 0.2 1.0 sethsbcolor % Light blue',
     $   '% 0 0 m 9000 0 l 9000 11000 l 0 11000 l closepath fill'
      endif
      return
      end
c_______________________________________________________________
      subroutine newpn(ipen)
c$$$$ calls nothing
c  Uses Postscript command to reset the color and line width.
c  Color table held in array below as hsb values
      parameter (ncol=20)
      common /hues/ table(3,ncol)
      save lastpn
      data lastpn/-9/
c
c  0 or 1: black; 2: red;   3: blue; 4: green; 5: brown;
c  6: orange;  7: yellow;  8: purple;  9:gray;    10: white.
c
c  Do not bother to change color if it doesn't change
      if (ipen .eq. lastpn) return
c  Flush current polygon buffer with previous parameters
      write(11,*) 'currentpoint k'
c  Ipen/1000 is interpreted as line width in 1/000 inches
      iwide=ipen/1000
      ip=min(max(mod(ipen, 1000), 1), ncol)
      if (iwide.gt.0) write(11,'(i5,a)') iwide,' setlinewidth'
      if (iwide.eq.0) write(11,'(i5,a)')     6,' setlinewidth'
c
c  By convention ipen=0, 1 => black
      write(11,'(3f8.3,a)') (table(i,ip),i=1,3),' sethsbcolor'
      lastpn=ipen
      return
      end
c_______________________________________________________________
      subroutine remark(text)
c$$$$ calls nothing
c  Inserts the remark text into the PostScript file
      character*(*) text
c
      write(11,'(2a)')'%',text
      return
      end
c_______________________________________________________________
      blockdata encaps
c  Sets a initial page size for encapsulated PS variables
      common /outsid/ box(4),xorig,yorig
      common /post/ boxx,boxy,land,infill
      data boxx,boxy/ 570.0,756.0/, land/0/, infill/0/
      data  box/1e6,1e6,-1e6,-1e6/
      data xorig,yorig/0.0,0.0/
      end
c_______________________________________________________________
      blockdata shades
c  Sets the colors in common /hues/
c 0 or 1: black; 2: red;   3: blue; 4: green; 5: brown;
c    6: orange;  7: yellow;  8: purple;  9:gray; 10: white.
c    11: light gray;    12-20: black
      parameter (ncol=20)
      common /hues/ table(3,ncol)
      data ((table(i,j),i=1,3), j=1, ncol)/
     $0,0,0,  0,1,1,  0.667,1,1,  0.400,1,.7, 0.1,1,.7,
     $0.1,0.8,1,  0.15,1,1,  0.833,1,1, 0,0,0.5,  0,0,1,
     $0,0,0.8,        27*0/
      end
c_______________________________________________________________