c********************PROFGRD DRIVER********************
c
c  This is the driver for subroutine SPROFGRD.
c  All files are opened, namelist &parms read in, and upon completion,
c  all files are closed.
c
c                         Written by Janet Bernard, USGS,
c                          Denver, Colorado, 1/80
c
c******************************************************************
c  The user has the option to input information by:
c     1)  Set up a command file to read in the namelist parms
c
c     2)  Input these parms directly from the terminal by typing
c         'term' when the command file name is requested
c
c  The following are read in by NAMELIST &parms
c     iseg           Number of connected line segments to be profiled
c                    iseg <= 10  (DEFAULT is 1)
c
c     title          For identification, first line of heading on graph.
c                    title <= 40 characters  (DEFAULT is gtitle) 
c
c     ddxy           Interpolation interval along any line segment.
c                    (DEFAULT is 0.5)
c
c     xscale         Units of x per inch:  If xscale=0.0, maximum
c                    board size used.  (DEFAULT is 0.0)
c
c     yscale         Units of y per inch:  If yscale=0.0, maximum
c                    board size used.  (DEFAULT is 0.0)
c
c     idrive         -1  for no plot
c                     5  HP7475
c                     8  CGA monitor output
c                     9  EGA monitor output (default)
c                    10  VGA monitor output
c
c     flag           flag=-1  coefficients generated and saved on
c                    coeffil, infil must be data grid
c
c                    flag=0  coefficients generated but not saved,
c                    infil must be data grid, coeffil=blank
c                    (DEFAULT is 0)
c
c                    flag=1  coefficients not generated in program
c                    infil must be file of coefficients, coeffil=blank
c
c                    NOTE:  for subsequent runs during the same
c                    execution, the user should be cautioned against
c                    using the same command file with the flag=-1.
c                    When the command segment is requested, type 'term',
c                    and input new namelist &parms.
c
c     coeffil        File name where computed coefficients of subgrid saved
c                    coeffil <= 50 characters (DEFAULT is 'profgrd.cof')
c
c     infil          Input file of coefficients or grid data, dependent
c                    on value of flag.
c                    infil <= 50 characters
c
c     binfil         Binary ouput file of interpolated data.
c                    binfil <= 50 characters  (DEFAULT is '    ',       
c                    no file created)
c
c     ascifil        ASCII output file of actual first and last
c                    x,y coordinate of each line segment plotted, and
c                    a list of all x,y coordinates used in the interpolation
c                    ascifil <= 50 characters (DEFAULT is '    ',
c                    no file created)
c
c     seg            array of iseg+1 x,y-coordinate pairs (DEFAULT
c                    is terminal input)
c
c  This program is terminated by typing 'ex' when the command
c  segment is requested or using a carraige return.
c
c  If the array 'seg' is not included in the namelist file,
c  coordinates of the line segments are read directly from the
c  terminal in 'v' format, one x-y pair per line.
c
c     NOTE:  These coordinates are GRID coordinates, with
c            xo <= x <= xo + number of columns * x spacing
c            yo <= y <= yo + number of rows * y spacing
c            (xo and yo are minimum x and y value at MESH point 1)
c
c***********************************************************************
c
c  BEGIN PROFGRD: 
c
c
c
      dimension xi(22),yi(22)
c
      integer flag
c
      character coeffil*50
      character blank*4, exit*4, infil*50, binfil*50, cmdseg*50, 
     -           gtitle*56, pgm*8, title*40, dumfil*50, ascifil*50
c       the following declaration added may, 1988 for terminal input
        character aiseg*10,adrive*5,aflag*5,addxy*15,axscale*15,
     1  ayscale*15
c
      common/prof/ infil,binfil,title,nrow,dx,dy,ddxy,xscale,
     -  yscale,idrive,inn,pgm,idum,iout,xo,yo,iasci,
     -  flag,gtitle,ascifil
c
c       changed namelist to common for use with namemc
c       will also need common area, prof, in namemc
      common/parms/ iseg, coeffil, seg(22)
c
c      data in/10/, inn/12/, ir/5/, iw/6/, blank/'   '/, exit/'ex'/,
c     -    iout/14/, idum/18/, iasci/16/, dval/1.0e+38/
c
        write (*,'(//1x,A)')'Entering USGS program, PROFGRD'
c  Set DEFAULT VALUES.
c
        in=10
        inn=12
        ir=5
        iw=6
        blank='   '
        exit='exit'
        iout=14
        idum=18
        iasci=16
        dval=1.0e+38
c
      xscale=0.0
        yscale=0.0
          idrive=9
            ddxy=0.5
              iseg=1
                title=blank
c
        binfil=blank        
          infil=blank
            ascifil=blank          
              coeffil=blank
                seg(1)=dval
                  flag=0
c
      dumfil='dumfil.dum'
      open(idum,access='direct',status='unknown',form='unformatted'
     1,file=dumfil,recl=64)
c
c  Read command segment.
c
10    write(iw,12) 'Enter command file name',
     1 '("term" for terminal input,',' "exit" or return to quit): '
12    format(/1x,A23/1x,a27,a28,$)
        read(ir,15) cmdseg
15        format(a50)
c
      if(cmdseg.eq.blank.or.cmdseg.eq.exit) go to 99
        if(cmdseg.eq.'term') go to 30
      open(in,status='old',form='formatted',file=cmdseg)
c
c  Read namelist parameter from terminal.
c
c      read(in, parms, end=97)
        call namemc(in)
	if (coeffil.eq.blank.and.flag.eq.-1)coeffil='profgrd.cof'
      close(in)
          go to 40
c
c       terminal input , added may 1988
c30    write(iw,32)
c32    format(' enter parms')
30      write (*,'(//1x,a)')
     1  'Entering input parameters.  Return selects default.'
        iseg=1
        write(*,'(/1x,a58,a15,$)')
     1  'Number of connected line segments to be profiled (<=10).  ',
     2  '(Default = 1): '   
        read (*,'(a10)') aiseg
        if (aiseg .ne. ' ') then
        call shiftr(aiseg,10)
        read(aiseg,'(i10)') iseg
        endif
        write(*,'(/1x,a)')
     1  'Graph title (<=40 chars, default is grid file identification):'
        read (*,'(a40)') title
        write (*,'(/1x,a,$)')
     1  'Interpolation interval along any line segment (default = 
     2  0.5): '
        read(*,'(a15)') addxy
        if (addxy .ne .' ') then
        call shiftr(addxy,15)
        read (addxy,'(g15.7)') ddxy
        endif
        write (*,'(/1x,a,$)')
     1  'Units x per inch. If 0.0, max board size used (default = 
     2  0.0): '
        read (*,'(a15)') axscale
        if  (axscale .ne. ' ') then
        call shiftr(axscale,15) 
        read (axscale,'(g15.7)') xscale
        endif
        write (*,'(1x,a,$)')
     1  'Units y per inch. If 0.0, max board size used (default = 
     2  0.0): '
        read(*,'(a15)') ayscale
        if (ayscale .ne. ' ') then
        call shiftr(ayscale,15)
        read (ayscale,'(g15.7)') yscale
        endif
        write (*,'(/1x,a,5(/5x,a),$)') 'Plotter:',' -1   No plot',
     1  '  5   HP7475','  8   CGA monitor output',
     2  '  9   EGA monitor output (default)',
     3  ' 10   VGA monitor output '
        write (*,'(//1x,a,$)') 'Enter plotter selection: '
        read (*,'(a5)') adrive
        if (adrive .ne. ' ') then
        call shiftr(adrive,5)
        read(adrive, '(i5)') idrive
        endif
        write (*,'(/1x,a,/(5x,a))') 'Output file options:',
     1  ' -1   coefficients generated and saved on coefficient file,',
     2  '      input file must be a data grid',
     3  '  0   coefficients generated but not saved',
     4  '      input file must be data grid, coefficient file name is bl
     5 ank', '  1   coefficients not generated in program, (default)',
     6  '      input file must be file of coefficients, coefficient file
     7  name blank'        
        flag=1
        write (*,'(/1x,a,$)') 'Enter file option (default = 1): '
        read (*,'(a5)') aflag 
        if (aflag .ne. ' ') then
        read (aflag,'(i5)') flag
        endif 
        if (flag .eq.-1)then
	write (*,'(/(1x,a),$)')
     1  'Name of file where computed coefficients are saved',
     2  '   <=50 chars, default = "profgrd.cof": '
        read (*,'(a50)') coeffil
        if (coeffil .eq. ' ') coeffil='profgrd.cof'
	endif
        write (*,'(2(/1x,a),$)')
     1  'Name of input file of coefficients or grid data',
     2  '   depending on file option (above), <= 50 chars: '
        read (*,'(a50)') infil
        write (*,'(/(1x,a),$)')
     1  'Name of binary output file of interpolated data,',       
     2  '   <= 50 chars (default is blank name, no file created): '
        read (*,'(a50)') binfil
        write (*,'(4(/1x,a),$)')
     1  'Name of ASCII output file of actual first and last x,y coordina
     2tes' ,'   of each segment plotted, and a list of all x,y coordinat
     3es' ,'   used in the interpolation (<= 50 chars, default is blank 
     4 name,','   no file created): '
        read (*,'(a50)') ascifil
c       finished terminal input of parameters
c        read(ir, parms)
c
40    if(iseg.gt.10) go to 9999
c
c
45    open(inn,status='old',form='unformatted',file=infil)
        read(inn) gtitle, pgm, ncol, nrow, nz, xo, dx, yo, dy
        write (*,'(//1x,a17,a56/1x,a12,a8,a20,i4,a13,i4,a6)')
     1  'Selected grid is ',gtitle,'created by ',pgm,
     2  ' and consisting of ',ncol,' columns and ',nrow,' rows.'
        write (*,'(/1x,a44,e12.5,a1,e12.5,a1)')
     1  'Coordinates of lower left corner of grid = (',xo,',',
     2  yo,')'
        write (*,'(1x,a5,e12.7,3x,a5,e12.7)')'dx = ',dx,' dy = ',dy
c
      if(nz.ne.16.and.flag.eq.1) go to 9777
      if(nz.ne.1.and.flag.ne.1) go to 9777
      if(dy.eq.0.0.and.flag.ne.1) go to 9666
      if(dx.eq.0.0.and.flag.ne.1) go to 9555
      if(title.eq.blank) title=gtitle
      if(gtitle.eq.blank)title='profile'
        pgm='profgrid'
          npoint=iseg+1
c
      if(flag.eq.-1) open(15,status='unknown',form='unformatted',
     1file=coeffil)
c
65    continue
      if(seg(1).ne.dval) go to 59
      write(iw, 50) npoint
50      format(//' Enter', i3,' nodes (x,y) defining line on which the',
     1 ' profile will be constructed'/
     -       '   Enter one (x,y) pair per line , e.g., 20.4,100 ')
      do 60 i=1, npoint
        read(ir,*) xi(i), yi(i)
60    continue
      go to 57
59    continue
      do 58 i=1,2*npoint,2
      j=i/2+1
      xi(j)=seg(i)
      yi(j)=seg(i+1)
58    continue
57    continue
c
      if(ascifil.ne.blank) open(iasci,status='unknown',form='formatted',
     1file=ascifil)
c
      if(binfil.ne.blank) open(iout,status='unknown',form='unformatted',
     1file=binfil)
        call sprofgrd (iseg,npoint,xi,yi,ncol)
          if(binfil.ne.blank)close(iout)
      close(inn)
              if(ascifil.ne.blank)close(iasci)
                if(flag.eq.-1)close(15)
                    go to 10
c
97    write(iw,98)
98      format('**odd eof in command segment')
      close(in)
            go to 10
c
9999  iseg1=iseg*2+2
      write(iw,9000) iseg1,iseg1,iseg,iseg,iseg,iseg
9000    format(' number of line segments exceeds limit:'/
     -' redimension the following arrays--'/
     -'  in profgrd:  xi(',i4,'), yi(',i4,')'/
     -'  in sprofgrd: ax(',i4,'), ay(',i4,'), bx(',i4,'), by(',i4,')'/ 
     -' and rerun program')
      go to 99
c
9777  write(iw,7000) nz
7000    format(' error in input file header record, nz=',i3)
          go to 10
c
9666  write(iw,6000) dy
6000    format(' error in input file header record, dy=',e16.7)
          go to 10
c
9555  write(iw,5000) dx
5000    format(' error in input file header record, dx=',e16.7)
          go to 10
 99   continue
      close(idum,status='delete')
       stop
       end
***************************SUBROUTINE SPROFGRID**************************
c
c  SUBROUTINE SPROFGRD  generates selected profiles from a contour map
c  via spline coefficients generated within the program from a subgrid
c  of the original grid, or generates selected profiles from
c  an input file of spline coefficients, depending on value of flag.
c
c*******************************************************************************
c
c          NOTE:     This program is a version of 'PROFGRID',
c                    written for the same purpose on the IBM 360 by
c                    W. A. Anderson, USGS, Denver, Colorado.
c
c                      This version revised by Janet Bernard, USGS,
c                      Denver, Colorado, 1/80.
c
c
c******************************************************************************
c
c
c     Selected lines are defined by coordinates at the end of a line segment
c     in terms of GRID UNITS.
c            ORIGIN (x=0,y=0) at MESH POINT (ix=1, iy=1)
c
c     The coordinates of A and B are BOTH expressed in terms of GRID UNITS;
c     if the ORIGIN is not (0, 0), the offset is added.  Thus, the
c     lower left grid point is (xo, yo), and the upper right corner
c     is (xo+(ncol-1)*dx, yo+(nrow-1)*dy).
c
c***********************LOCAL VARIABLES********************************
c
c
c     POINTS IN THE PARAMETERS ax, ay, bx, by, MUST SATISFY THE FOLLOWING
c     RELATIONSHIPS FOR index=1..iseg:
c
c          LET
c                    axi=ax(index)
c                    ayi=ay(index)
c                    bxi=bx(index)
c                    byi=by(index)
c
c          THEN
c                    axi.le.bxi
c                    dx.le.bxi.le.ncol*dx
c                    dx.le.axi.le.ncol*dx
c                    dy.le.byi.le.nrow*dy
c                    dy.le.ayi.le.nrow*dy
c                    axi.ne.bxi.and.ayi.ne.byi
c
c     This means point A MUST be to the LEFT or directly ABOVE or BELOW
c     point B.
c
c
c     nrow           Number of rows in original data grid.
c
c     ncol           Number of columns in original data grid.
c
c     dx             Original GRID x-spacing.  (DEFAULT is 1.0)
c
c     dy             Original GRID y-spacing.  (DEFAULT IS 1.0)
c
c     dval           Grid rejection tolerance, set to 'ffff7fff'x
c
c     nc1            Original no. columns in data grid less one
c
c     nr1            Original no. rows in data grid less one
c
c     xmax           Maximum x-coordinate value on GRID.
c
c     ymax           Maximum y-coordinate value on GRID.
c
c     coef           ARRAY with SPECIFIC spline coefficients for selected
c                    cell.
c
c     npts           Counter of points interpolated on EACH line
c                    segment.  (maximum per line segment is 750,
c                    TOTAL MAXIMUM is 7500)
c
c     nptlast        Last count of points interpolated.  Initial value is 0.
c                    Used as offset to index of X and Y ARRAYS.
c
c     maxpts         Total number of points interpolated from all line
c                    segments.
c
c     ifirst         Used as pointer when copying xp and yp values into
c                    ARRAYS xplot and yplot. 
c
c     inumber        Counter of ACTUAL line segments plotted.
c
c     d              Delta x.  (point B minus point A)
c
c     slope          Slope of profile line segment.
c
c     theta          Angle between profile line segment and X-axis
c
c     ddx            Increment on X-axis for next interpolation point.
c
c     ddy            Increment on Y-axis for next interpolation point.
c
c     xp             ARRAY of intermediate, actual length of line segment.
c                    Initial value equals zero.  Values to be
c                    X-coordinates passed to plot graph of profile.
c
c     xplast         Last value of length of line segment plotted.
c                    Initial value equals zero.  Used as an offset for
c                    X values to plot a continuous graph for iseg>1.
c
c     yp             ARRAY of 'z' values interpolated with bicubic spline
c                    function.  Values to be Y-coordinates passed to plot
c                    graph of profile.
c
c     xplot          ARRAY of all xp values, passed to plot graph.
c
c     yplot          ARRAY of al yp values, passed to plot graph.
c
c     x              X in 'x-x(i-1)'
c
c     y              Y in 'y-y(i-1)'
c
c     ix             X-coordinate of MESH point for selected cell of
c                    coefficients.
c
c     iy             Y-coordinate of MESH point for selected cell of
c                    coefficients.
c
c     irec           Key to read correct set of coefficients from indexed file.
c
c     lrec           Key to last set of coefficients used.
c
c     xx             Value of 'x-x(i-1)' in bicubic spline function.
c
c     yy             Value of 'y-y(i-1)' in bicubic spline function.
c
c     beta           ARRAY of intermediate results in interpolation of 'z'
c                    values with bicubic spline function.
c
c     nseg           Index for main DO loop of subprogram PROFILE.
c
c     i,j,ic,jc,k    Indices for DO and IMPLIED DO loops.
c
c     nptplot        ARRAY of number of points interpolated for each line
c                    segment.
c
c     startx         ARRAY of x-coordinates of actual starting point.
c
c     starty         ARRAY of y-coordinates of actual starting point.
c
c     endx           ARRAY of x-coordinates of actual ending point.
c
c     endy           ARRAY of y-coordinates of actual ending point.
c
c     delx,dely      Delx and dely are used to determine the start
c                    of the next line segment when two line segments are
c                    connected so as to preserve the 'ddxy' interval
c
c  The next two arrays are dimensioned in SPROFGRD  and used for
c  computed values in subroutine PRPLOT:
c
c     px and py
c
c******************************************************************************
c
c     BEGIN SUBROUTINE SPROFGRD: 
c
c
c
      subroutine sprofgrd (iseg,npoint,xi,yi,ncol)
c
       character infil*50, title*40, binfil*50, pltfil*50,  pgm*8,
     -  blank*4, ascifil*50
       character titl*56, pgmid*8, gtitle*56
c
       real minx,maxx,miny,maxy
       integer subcol1, subrow1, flag
       integer rowmin,rowmax,colmin,colmax        
c
        dimension yyo(150),splcoef(16,150)
c      dimension yyo(250),splcoef(16,250)
c       reduce dims of 20000 to 7500,100 to 500,2000 to 750
      dimension px(7500),py(7500),coef1(16,500),extra(500)
      dimension coef(16), xp(750), yp(750),
     -   xplot(7500), yplot(7500), beta(4), nptplot(10),
     -   ax(10), ay(10), bx(10), by(10),
     -   startx(10), starty(10), endx(10), endy(10),
     -   xi(npoint), yi(npoint)            
c
      common/prof/infil,binfil,title,nrow,dx,dy,ddxy,xscale,
     -  yscale,idrive,inn,pgm,idum,iout,xo,yo,iasci,
     -  flag, gtitle, ascifil
c
c       reduce subgrid to 150 x 100  mrm 5/88
        common/matrix/u(150,100),p(150,100),q(150,100),s(150,100)
c      common/matrix/u(250,200),p(250,200),q(250,200),s(250,200)
c
c       change storage locations to match new dimensions
      equivalence (coef1(1),u(1,1)), (extra(1),u(86,1)), (px(1),p(1,1)),
     -  (py(1),p(30,1)), (xplot(1),q(1,1)), (yplot(1),q(30,1)),
     -  (splcoef(1),s(1,1)), (yyo(1),s(30,1))
c      equivalence (coef1(1),u(1)), (extra(1),u(13000)), (px(1),p(1)),
c     -  (py(1),p(7501)), (xplot(1),q(1)), (yplot(1),q(7501)),
c     -  (splcoef(1),s(1)), (yyo(1),s(7501))
c
c         data dval/'ffff7fff'x/
        data dval/.1701412e+39/
         data lrec/-1/, blank/'    '/
c
c  Check if arrays extra and coef1 in range
c  If an error message is received, redimension the appropiate array
c  and change the following statements which apply:
c
      if(ncol.gt.500.and.flag.eq.1) go to 7770
      if(nrow.gt.500.and.flag.eq.1) go to 7771
c
c  Put data into ax and ay arrays and check if x-coordinates in
c  ascending order:
c
c
      do 10 i=1, iseg
        ax(i)=xi(i)
          ay(i)=yi(i)
            bx(i)=xi(i+1)
              by(i)=yi(i+1)
                if(ax(i).gt.bx(i)) go to 9006
   10 continue
      pltfil=binfil
        if(binfil.eq.blank) pltfil=infil
c
c  Check for PARAMETER errors.
c
       if(nrow.le.1.or.ncol.le.1.or.dx.le.0.0.or.dy.le.0.0)          
     -   go to 9001
c
c  Compute maximum values.
c
       nc1=ncol-1
         nr1=nrow-1
            xmax=xo+nc1*dx
             ymax=yo+nr1*dy
c
c  Initialize counters and offsets to zero
c
      nptlast=0
        xplast=0.0
          maxpts=0
            delx=0.0
              dely=0.0
c
c  Check that line segments lie within grid
c
      do 50 nseg=1,iseg
        if(ax(nseg).lt.xo.or.ax(nseg).gt.xmax) go to 9002
        if(bx(nseg).lt.xo.or.bx(nseg).gt.xmax) go to 9002
        if(ay(nseg).lt.yo.or.ay(nseg).gt.ymax) go to 9002
        if(by(nseg).lt.yo.or.by(nseg).gt.ymax) go to 9002
        if(ax(nseg).eq.bx(nseg).and.ay(nseg).eq.by(nseg)) go to 9002
50    continue
c
c  Read in spline coefficients on keyed sequential file if flag=1
c
      if(flag.ne.1) go to 700
      n=0
      do 1 i=1,nr1
        read(inn,end=99)extra(i),((coef1(k,ic),k=1,16),ic=1,nc1)
        do 1 jc=1,nc1
          n=n+1
          write(idum,rec= n) (coef1(k,jc),k=1,16)
   1  continue
      minx=xo
        miny=yo
          maxx=xmax
            maxy=ymax
              go to 75
c
700   continue
c
c**********compute grid coordinates of subgrid**********
c
c  minx is subgrid xo
c  maxx is subgrid xmax
c  miny is subgrid yo
c  maxy is subgrid ymax
c
c    The following parameters are computed in SUB_GRD and used in
c    BICUBIC_SPL to compute coefficients
c
c  rowmin is first row to be read from input 'z' values
c  rowmax is last row to be read from input 'z' values
c  colmin is first column to be read from input 'z' values
c  colmax is last column to be read from input 'z' values
c
c*******************************************************
c
      maxx=-dval
      minx=dval
      maxy=-dval
      miny=dval
      do 800 nseg=1,iseg
      if(ax(nseg).gt.maxx)maxx=ax(nseg)
      if(ay(nseg).gt.maxy)maxy=ay(nseg)
      if(ax(nseg).lt.minx)minx=ax(nseg)
800   if(ay(nseg).lt.miny)miny=ay(nseg)
      if(bx(iseg).gt.maxx)maxx=bx(iseg)
      if(bx(iseg).lt.minx)minx=bx(iseg)
      if(by(iseg).gt.maxy)maxy=by(iseg)
      if(by(iseg).lt.miny)miny=by(iseg)
      maxx=maxx+10.*dx
      if(maxx.gt.xmax)maxx=xmax
      minx=minx-10.*dx
      if(minx.lt.xo)minx=xo
      maxy=maxy+10.*dy
      if(maxy.gt.ymax)maxy=ymax
      miny=miny-10.*dy
      if(miny.lt.yo)miny=yo
      rowmin=ifix((miny-yo)/dy+1.0000001)
      rowmax=ifix((maxy-yo)/dy+1.0000001)
      colmin=ifix((minx-xo)/dx+1.0000001)
      colmax=ifix((maxx-xo)/dx+1.0000001)
c
c**********compute spline coefficients**********
c
c  compute number of columns and rows in subgrid
      nc=colmax-colmin+1
      nr=rowmax-rowmin+1
c
c  Check if dimensions of subgrid is within array bounds.
c
c     NOTE:  The maximum dimensions of subgrid are 250 by 200.
c            This program cannot be used for any thing larger.
c
c       changed size of subgrid mrm 5/88
        if (nc .gt. 150 .or. nr.gt.100) goto 9004
c      if(nc.gt.250.or.nr.gt.200) go to 9004
c
c  Compute spline coefficients
c
      call bicubspl(nc,nr,dx,dy,dval,rowmin,colmin)
c
*********************************************
c
c  subcol1 is number of columns in subgrid less one
c  subrow1 is number of rows in subgrid less one
      subcol1=nc-1
      subrow1=nr-1
c
c  Save coefficients if flag equals -1:
c
      if(flag.ne.-1) go to 75
        nza=16
          pgmid='splncoef'
            index=0
              yyo(1)=yo
                write(15) gtitle,pgmid,nc,nr,nza,minx,dx,miny,dy
                do 500 i=2,nr
                  yyo(i)=yyo(i-1)+dy
500             continue
                do 550 i=1,subrow1
                  do 576 j=1,subcol1
                    index=index+1
                    read(idum,rec= index) (splcoef(k,j),k=1,16)
                    if(j.ne.subcol1) go to 576
                    do 575 jj=1,16
                      splcoef(jj,nc)=dval
575               continue                    
576               continue
                  write(15) yyo(i),((splcoef(kk,ii),kk=1,16),ii=1,nc)
550             continue
      do 600  i=1,16
        do 600 j=1,nc
          splcoef(i,j)=dval
600   continue
      write(15) yyo(nr),((splcoef(kk,ii),kk=1,16),ii=1,nc)
c
75    continue
c
c  BEGIN LOOP TO INTERPOLATE VALUES.
c
c
       do 100 nseg=1, iseg
        npts=0
         axi=ax(nseg)
           ayi=ay(nseg)
             bxi=bx(nseg)
               byi=by(nseg)
c
c  Compute DELTA X.  If DELTA X is zero, line segment is parallel to Y-axis.
c  Jump to set related X and Y increments.
c
       d=bxi-axi
         if(d.eq.0.0) go to 200
c
c  Compute slope and increments of non-vertical line segments.
c
      slope=(byi-ayi)/d
        theta=atan(slope)
          ddx=ddxy*cos(theta)
            ddy=ddxy*sin(theta)
              if(nseg.eq.1) go to 202
                dist=sqrt((x-axi)**2+(y-ayi)**2)
                  delx=dist*cos(theta)
                    dely=dist*sin(theta)
                      go to 202
c
c  Set increments for vertical line segment.
c
  200  continue
      slope=1.0e37
        ddx=0.0
          ddy=ddxy
            if(ayi.gt.byi) ddy=-ddy
              if(nseg.eq.1) go to 202
                delx=0.0
                  dely=sqrt((x-axi)**2+(y-ayi)**2)
                    if(ayi.gt.byi) dely=-dely
c
c  Initialize X and Y coordinates.
c
  202  continue
      axi=axi+delx
        ayi=ayi+dely
          x=axi
            y=ayi
c
c  TOP OF INNER LOOP FOR EACH LINE SEGMENT.
c
  3   continue
        npts=npts+1
c
c  Check if maximum number of points interpolated.
c  If an error message is received, redimension arrays
c  xp,yp,and px,py,xplot,yplot if needed.
c
c  NOTE:  Arrays px,py,xplot, and yplot would not have
c         to be redimensioned if the maximum number of line segments
c         ie, 10, is not used.
c
c
      if(npts.gt.750) go to 9003
           if(y.eq.maxy) y=y-1.0e-5
             if(x.eq.maxx) x=x-1.0e-5
                 xp(npts)=sqrt((x-axi)**2+(y-ayi)**2)+xplast
c
c  Compute MESH point coordinates.
       iy=int((y-miny)/dy)+1
         ix=int((x-minx)/dx)+1
c
c  Compute KEY for record.
c
       irec=(iy-1)*subcol1+ix
         if(irec.eq.lrec) go to 8
           read(idum,rec= irec) (coef(j), j=1,16)
             lrec=irec
c
c  Compute 'x-x(i-1)' and 'y-y(i-1)'.
c

    8 continue
         xx=(x-minx)-float(ix-1)*dx
           yy=(y-miny)-float(iy-1)*dy

c  Compute intermediate values of bicubic spline function,
c  and check for missing data.  The first coefficient has been
c  set to dval if missing data in cell.
c
      if(coef(1).ge.1.0e+30) go to 11
      do 9 j=1,4
        beta(j)=((coef(j+12)*yy+coef(j+8))*yy+coef(j+4))*yy+coef(j)
    9 continue
        yp(npts)=((xx*beta(4)+beta(3))*xx+beta(2))*xx+beta(1)
        go to 12
11    yp(npts)=dval
c
c  Increment X and Y in 'x-x(i-1)' and 'y-y(i-1)'.
c
12    x=x+ddx
        y=y+ddy
c
c  Check for end of inner loop.  Jump to different check if horizontal
c  line segment.
c
      if(slope.eq.0.0) go to 31
        if(y.le.byi.and.byi.gt.ayi) go to 3
          if(y.ge.byi.and.byi.lt.ayi) go to 3
              go to 33
c
   31  continue
        if(x.le.bxi) go to 3
      if(ascifil.eq.blank) go to 275
c
c  DONE WITH LINE SEGMENT(nseg).
c
c  Save points xp and yp in ARRAYS xplot and yplot for graph.
c  Initialize index k to zero.
c
   33 continue
      xplast=xp(npts) + ddxy
        ifirst=nptlast+1   
          maxpts=maxpts+npts
            k=0
      do 300 j=ifirst, maxpts
        k=k+1
          xplot(j)=xp(k)
            yplot(j)=yp(k)
  300 continue
c
c  Save number of points interpolated for line segment(nseg).
c
      nptplot(nseg)=npts
c
c  Save grid coordinates of actual start and end of line segment profiled.
c
      startx(nseg)=axi
        starty(nseg)=ayi
          endx(nseg)=x-ddx
            endy(nseg)=y-ddy
c
c  Increment nptlast.  nptlast becomes ifirst when
c  indexing ARRAYS xplot and yplot.
c
      nptlast=nptlast+npts
c
c  JUMP TO TOP OF LOOP FOR SPECIFICATIONS OF NEXT LINE SEGMENT.
c
  100 continue
c
c  DONE WITH ALL LINE SEGMENTS.
c  Save ALL information on file.
      if(ascifil.eq.blank) go to 275
      write(iasci,3333) infil,title,binfil,ddxy,dx,dy,minx,miny
3333  format(5x,'infil:',5x,a50/5x,'title:',5x,a56/4x,'binfil:',5x,a50/
     -6x,'ddxy:',  e16.7/8x,'dx:',  e16.7/8x,'dy:',  e16.7/
     -8x,'xo:',  e16.7/8x,'yo:',  e16.7)  
c
      do 411 i=1,iseg
        write(iasci,4444)i,startx(i),starty(i),endx(i),endy(i),          
     -nptplot(i)
4444    format(//2x,'line segment ',i2,':'/6x,'start x&y=',2(e16.7)/
     -8x,'end x&y=',2(e16.7)/2x,'number points=',i4/)
      if(i.ne.1) go to 68
      ibegin=1
      iend=7 
        number=nptplot(1)
          go to 4000
68    ibegin=number+1
      iend=number+7
      number=number+nptplot(i)
4000  if(iend.gt.number) iend=number
      write(iasci,5555)(xplot(j),j=ibegin,iend)
5555    format(5x,'xplot:',5x,7(e16.7))
      write(iasci,5566)(yplot(k),k=ibegin,iend)
5566    format(6x,'data:',5x,7(e16.7)/)
      if(iend.eq.number) go to 411
      ibegin=iend+1
        iend=iend+7 
          go to 4000
411   continue
275   continue
      if(binfil.eq.blank) go to 250
        nrprf=1
          dum=0.0
            nz=1
c              dxprf=0.0
c                xoprf=ax(1)  
c                  yoprf=ay(1)  
      xoprf=0.0
      yoprf=0.0
      titl=title
      write(iout) titl,  pgm, maxpts, nrprf, nz, xoprf, ddxy,
     -     yoprf, ddxy
c        write(iout) (xplot(i), i=1,maxpts)
          write(iout) dum, (yplot(i), i=1,maxpts)
c
c  Call plot routine.
c
250   if (idrive .ge. 0) call prplot(xplot, yplot, maxpts, 
     1 xscale, yscale, idrive,title,pltfil,ddxy,px,py)
c
c  RETURN TO MAIN PROGRAM.
c
      return
c
c  ERROR MESSAGES.
c
7770  write(6,7000) ncol
7000    format(' redimension array coef1 to coef1(16,',i5,')',/
     -' and rerun program')
        go to 999
c
7771  write(6,7001) nrow
7001    format(' redimension array extra to extra(',i6,')',/
     -' and rerun program')
        go to 999
c
 9006 write(6,6666)
 6666   format(' x coordinates must be in ascending order')
          go to 999
c
 9004 write(6,9994)nc,nr
 9994    format(' dimension of subgrid(',i6,',',i6,'), exceeds limit'/
     -' this program cannot be used')
c
 9002 write(6,9999)
 9999   format(' line segment lies outside grid area.')
          go to 999
c
c
 9001 write(6,9998)
9998   format('  paramater error in namelist')
          go to 999
c
c
c
9003  write(6,8000) npts
8000  format(' maximum points per segment is 750, redimension arrays'/
     -' xp and yp greater than ',i5,' and rerun program'/
     -' NOTE:  if number of segments plotted equals 10, arrays'/
     -' px, py, xplot, and yplot must be redimensioned also')
        go to 999
c
   99 write(6,109)
  109   format(' eof reached before completion of data record.')
c
c  Return to calling program.
c
999   return
      end








********************SUBROUTINE BICUBIC_SPL********************
c
c  Subroutine BICUBIC_SPL computes spline coefficients from a selected
c  subgrid, to be used for selected profiles.
c
c      Based on splin2 by W. L. Anderson, USGS,
c      revised by Janet Bernard, USGS, Denver, Colorado, 1/80.
c
c******************************************************************
c
c  Variables used:
c
c     dat            Array of data points read from input file
c
c     u              Array into which selected data for subgrid is read
c
c     p,q,s          Work arrays into which partial derivatives are put
c
c     rowmin         First row of data read from input grid
c
c     rowmax         Last row of data read from input grid
c
c     colmin         First column of data read from input grid
c
c     colmax         Last column of data read from input grid
c
c     nc             Number of columns in subgrid used
c
c     nr             Number of rows in subgrid used
c
c     dx             x interval of grid
c
c     dy             y interval of grid
c
c     z              Array of 'z' values passed to CUBIC
c
c     d              First or second derivative array
c
c     t              Flag for first or second derivatives.
c                    t=0 for second derivative boundary conditions for
c                    first 4 one dimensional splines.
c                    t=1 for first derivative boundary conditions for
c                    last one dimensional spline.
c
c     ax,ay,wk       Intermediate work arrays to compute coefficients
c
c     dd             Array used for 16 spline coefficients
c
c     idum           file18, keyed sequential file onto which each set
c                    of 16 coefficients read
c
c     dx1,dx2,dx3    Based on x interval, used to set ax array
c
c     dy1,dy2,dy3    Based on y interval, used to set ay array
c
c
c  Subroutine BICUBIC_SPL calls subroutines DTEST, CUBIC, and MATMULT.
c     DTEST          Sets dvals
c
c     CUBIC          Computes one dimensional splines on given column or row
c
c     MATMULT        Multiplies 2 x 2 matrices
c
c  The following arrays are dimensioned in BICUBIC_SPL but used in
c  CUBIC:
c
c     a,b,c,w1       Results are returned in array c
c
c
c*********************************************************
c
c  BEGIN BICUBIC_SPL:
c
      subroutine bicubspl(nc,nr,dx,dy,dval,rowmin,colmin)
        integer t,rowmin,colmin,rowmi1,colmi1
      dimension ddum(500),
     -    a(250),b(250),c(250),z(250),w1(250),
     -    d(2), dd(16), wk(16), ax(16), ay(16)
c
c
c       change subgrid dimension mrm 5/88
        common/matrix/u(150,100),p(150,100),q(150,100),s(150,100)
c      common/matrix/u(250,200),p(250,200),q(250,200),s(250,200)
        real kk(16)
      data ax/1.0,4*0.0,1.0,10*0.0/
      data ay/1.0,4*0.0,1.0,10*0.0/
c
      colmi1=colmin-1
      rowmi1=rowmin-1
c
c****Check dimensions of array ddum.  ddum is a 'dummy' array
c    used to read in preceding data on a given row, until the data
c    which is required is reached.  If an error message is received,
c    redimension ddum greater than or equal to colmi1.
c
c    NOTE:  Change the following test statement to match the
c           new dimension:
c
      if(colmi1.gt.500) go to 9999
c
c****input data****
c
      if(rowmi1.eq.0.and.colmi1.eq.0) go to 10
      if(rowmi1.eq.0) go to 30
      if(colmi1.eq.0) go to 20
c**rowmi1 and colmi1 are not equal to zero
      do 500 i=1,rowmi1
        read(12) dum
500   continue
c**rowmi1 equals zero, colmi1 not equal zero
 30   continue
      do 600 i=1,nr
        read(12) dum,(ddum(j),j=1,colmi1), (u(j,i),j=1,nc)
600   continue
        go to 25
c**colmi1 equals zero, rowmi1 not equal zero
 20   continue
      do 700 i=1,rowmi1
        read(12) dum
700   continue
      do 800 i=1,nr
          read(12) dum, (u(j,i),j=1,nc)
800   continue
        go to 25
c**romi1 and colmi1 both equal zero
 10   continue
      do 1000 i=1,nr
        read(12) dum,(u(j,i),j=1,nc)
1000  continue
c****************
 25   continue
      nc1=nc-1
      nr1=nr-1
c**second derivative boundaries assumed 0.0
      d(1)=0.0
      d(2)=0.0
      t=0
      n=0
      idum=18
      do 5 j=1,nr
        do 4 i=1,nc
          z(i)=u(i,j)
 4      continue
        call dtest(z,nc)
        call cubic(nc,dx,z,a,b,c,t,d,w1)
c**check for error return**
        if(nc.lt.0) return
        do 41 l=1,nc
          p(l,j)=a(l)
41      continue
 5    continue
      do 7 i=1,nc
        do 6 j=1,nr
          z(j)=u(i,j)
 6      continue
        call dtest(z,nr)
        call cubic(nr,dy,z,a,b,c,t,d,w1)   
c**check for error return**
        if(nr.lt.0) return
        do 61 l=1,nr
          q(i,l)=a(l)
61      continue
 7    continue
      do 8 i=1,nc
        z(i)=q(i,1)
 8    continue
      call cubic(nc,dx,z,a,b,c,t,d,w1)   
c**check for error return**
      if(nc.lt.0) return
      do 81 l=1,nc
        s(l,1)=a(l)
81    continue
      do 9 i=1,nc
        z(i)=q(i,nr)
 9    continue
      call cubic(nc,dx,z,a,b,c,t,d,w1)   
c**check for error return**
      if(nc.lt.0) return
      do 91 l=1,nc
        s(l,nr)=a(l)
91    continue
c**flag reset for first derivative boundaries**
      t=1
      do 100 i=1,nc
        do 95 j=1,nr
          z(j)=p(i,j)
95      continue
        d(1)=s(i,1)
        d(2)=s(i,nr)
        call cubic(nr,dy,z,a,b,c,t,d,w1)   
c**check for error return**
        if(nr.lt.0) return
        do 99 l=1,nr
          s(i,l)=a(l)
99      continue
100   continue
c**set the ax array**
      dx1=1.0/dx
      dx2=dx1*dx1
      dx3=dx2*dx1
      ax(3)=-3.0*dx2
      ax(4)=2.0*dx3
      ax(7)=-2.0*dx1
      ax(8)=dx2
      ax(11)=-ax(3)
      ax(12)=-ax(4)
      ax(15)=-dx1
      ax(16)=dx2
c
c**set the ay array**
      dy1=1.0/dy
      dy2=dy1*dy1
      dy3=dy2*dy1
      ay(9)=-3.0*dy2
      ay(10)=-2.0*dy1
      ay(11)=-ay(9)
      ay(12)=-dy1
      ay(13)=2.0*dy3
      ay(14)=dy2
      ay(15)=-ay(13)
      ay(16)=dy2
c**if a dval is encountered in the data grid points, the first
c  coefficient is set to dval as a flag**
      do 200 j=1,nr1
        j1=j+1
        do 200 i=1,nc1
          i1=i+1
          if(u(i,j).ge.1.0e+30.or.u(i1,j).ge.1.0e+30.or.u(i,j1)
     -.ge.1.0e+30.or.u(i1,j1).ge.1.0e+30) go to 198
          kk(1)=u(i,j)
          kk(2)=p(i,j)
          kk(3)=u(i1,j)
          kk(4)=p(i1,j)
          kk(5)=q(i,j)
          kk(6)=s(i,j)
          kk(7)=q(i1,j)
          kk(8)=s(i1,j)
          kk(9)=u(i,j1)
          kk(10)=p(i,j1)
          kk(11)=u(i1,j1)
          kk(12)=p(i1,j1)
          kk(13)=q(i,j1)
          kk(14)=s(i,j1)
          kk(15)=q(i1,j1)
          kk(16)=s(i1,j1)
          call matmult(ax,kk,wk,4,4,4)
          call matmult(wk,ay,dd,4,4,4)
          go to 199
198         dd(1)=dval
199   continue
      n=n+1
c****write coefficients on keyed sequential file****
      write(idum,rec= n)(dd(m),m=1,16)
300       continue
200   continue
      return
c
c****Write error message and set flag for error return****
c
9999  write(6,9009) colmi1
9009    format(' Array ddum in subroutine BICUBIC_SPL must'/
     -' be redimensioned greater than or equal to ', i3)
      nr=-1
      return
      end
c********************SUBROUTINE CUBIC**********************
c
c  One dimensional cubic spline interpolation, called by
c  subroutine BICUBIC_SPL.
c
c          Based on splin1, by W. L. Anderson, USGS,
c          revised by Janet Bernard, USGS, Denver, Colorado, 1/80.
c
c*************************************************************
c
c  Parameters used:
c
c     m              Number of data points, m > 2
c
c     h              Length of EQUALLY spaced interval between data points
c
c                    NOTE:  If grid is unequally spaced, program SPLN2D
c                           must be used and then, program PROFGRD,
c                           not PROFGRD2.
c
c     a,b,c          Coefficient arrays.  Results are returned in first
c                    (m-1) elements of a, b, c.
c
c     p              Work array
c
c     t              Type of boundary condition supplied in d array.
c                    t=1, first derivatives given at end points
c                    t=0, second derivatives given at end points
c
c     d              Derivative array
c
c     y              Data of selected row or column
c
c**************************************************************
c
c  BEGIN CUBIC:
c
       subroutine cubic(m,h,y,a,b,c,t,d,p)  
         dimension y(1), a(1), b(1), c(1), p(1)
         dimension d(2)
         real mul
         integer t
c**check for parameter error**
           if(t.lt.0.or.m.lt.3.or.h.eq.0.0) go to 600
             hh=3.0/h
             n=m-1
             if(t.eq.0) go to 500
c*****1st derivative boundaries given*****
             ne=n-1
             do 100 j=1,ne
               b(j)=4.0
               c(j)=1.0
               a(j)=1.0
               p(j)=hh*(y(j+3)-y(j))
100          continue
             p(1)=p(1)-d(1)
             p(ne)=p(ne)-d(2)
             go to 700
c*****2nd derivative boundaries given*****
500      ne=n+1
         do 200 j=2,n
           b(j)=4.0
           c(j)=1.0
           a(j)=1.0
           p(j)=hh*(y(j+1)-y(j-1))
200       continue
           b(1)=2.0
           b(ne)=2.0
           c(1)=1.0
           c(ne)=1.0
           a(ne)=1.0
           p(1)=hh*(y(2)-y(1))-.5*h*d(1)
           p(ne)=hh*(y(m)-y(n))+.5*h*d(2)
c*****solution of tridiagonal matrix*****
700      fa=1.0/b(1)
         c(1)=c(1)*fa
         p(1)=p(1)*fa
           do 300 j=2,ne
             mul=1.0/(b(j)-a(j)*c(j-1))
             c(j)=mul*c(j)
             p(j)=mul*(p(j)-a(j)*p(j-1))
300        continue
c*****obtain spline coefficients*****
         a(ne+t)=p(ne)
         j=ne-1
800             continue
                a(j+t)=p(j)-c(j)*a(j+t+1)
                j=j-1
                if(j.ge.1) go to 800
         if(t.eq.0) go to 900
           a(1)=d(1)
           a(m)=d(2)
900      continue
         hh=1.0/h
         do 1000 j=1,n
           mul=hh*(y(j+1)-y(j))
           b(j)=hh*(3.0*mul-(a(j+1)+2.0*a(j)))
           c(j)=hh*hh*(-2.0*mul+a(j+1)+a(j))
1000     continue
         return
c*****error return*****
600      m=-abs(m)
        write (6,999)
999       format(' paramater error in subroutine cubic')
         return
         end
c********************SUBROUTINE DTEST************************
c
c    by W. L. Anderson, USGS
c
c   Subroutine DTEST is called by BICUBIC_SPL and presets row
c   or column vector containing any dvals.
c
c**************************************************************
c
      subroutine dtest(y,lim)
      dimension y(1)
      if(y(1).ge.1.0e+30) go to 30
   10 do 20 i=2,lim
      if(y(i).lt.1.0e+30) go to 20
      y(i)=y(i-1)
      if(i.eq.lim) go to 20
      if(y(i+1).lt.1.0e+30) y(i)=y(i+1)
   20 continue
      return
   30 do 40 i=2,lim
      if(y(i).ge.1.0e+30) go to 40
      y(1)=y(i)
      go to 10
   40 continue
      return
      end
c********************SUBROUTINE MATMULT************************
c
c  multiplication of an (n x m) matix by a (m x l) matrix
c
c  by Janet Bernard, 10/79, USGS
c
c**parameters**
c
c     a              First input matrix
c
c     b              Second input matrix
c
c     c              Output matrix
c
c     n              Number of rows in matrix a
c
c     m              Number of columns in matrix a and rows in matrix b
c
c     l              Number of columns in matrix b
c
c***************************************************************
c
c  BEGIN MATRIX_MULT:
c
        subroutine matmult(a,b,c,n,m,l)    
        dimension a(n,m), b(m,l), c(n,l)
          do 100 i=1,n
            do 100 j=1,l
              c(i,j)=0.0
              do 100 k=1,m
                c(i,j)=a(i,k)*b(k,j)+c(i,j)
100       continue
        return
        end














c*********************SUBROUTINE PRPLOT************************
c
c  Subroutine PRPLOT is a version of AMPLOT.  It is called from
c  SPROFGRD and SPROFGRD2, and is the driver for CRTPLT.
c  This routine plots profiles of field data from selected line segments
c
c***********************************************************
c
c  BEGIN PRPLOT:
c
      subroutine prplot(x, y, max, xscale, yscale, idrive,
     -    title,pltfil,ddxy,xp,yp)
      dimension xp(max),yp(max),titles(5)
      dimension x(max),y(max),xx(2),yy(2)
      character titles*40,title*40,pltfil*50
      titles(3)=title
      iw=6
c
c     determine min & max values of amplitutes and distances.
c
      bval=9999999.0
      ymin=1.0e30
      ymax=-1.0e30
      do 10 i=1,max
      if(y(i).gt.bval) go to 10
      if(y(i).lt.ymin)ymin=y(i)
      if(y(i).gt.ymax)ymax=y(i)
10    continue
      dif=(ymax-ymin)*0.1
      ymax=ymax+dif
      ymin=ymin-dif
      print*,'1449 ',ymin
      dati=ymin-1.0000001
      print*,'1451 ',dati
      idatm=ifix(dati)
      ymin=float(idatm)
      dati=ymax+1.0000001
      idatm=ifix(dati)
      ymax=float(idatm)
      if(ymax.ne.ymin) go to 12
      ymax=dati
12    xmax=x(max)
c
c   set plot parameters
c
      dum=0.0
      x3=0.
      y3=1.
      m=-2
      iplot=1
      isym=20
      xx(1)=0.
      xx(2)=xmax
c       commented out tektronix jun 88
c      if(idrive.eq.1.or.idrive.eq.4) go to 40
c      dati=xx(2)/xscale+0.5
c      nxinch=ifix(dati)
c      dati=(ymax-ymin)/yscale+0.5
c      nyinch=ifix(dati)
c      nxl=nxinch*4+1
c      nyl=nyinch*4+1
c      go to 50
40    nxl=24
      nyl=17
50    yy(1)=ymin
      yy(2)=ymax
      titles(1)='distance                                '
      titles(2)='amplitude                               '
      write(titles(4),60)pltfil
60    format(' segment name= ', a25)
      write(titles(5),70)max,ddxy
70    format(' no. values= ',i5,2x,'interval=',f8.3,4x) 
      nt=5
      ndec=2
      if(xmax.ge.100000.or.xmax.lt.1.)ndec=-2
      iclose=0
      if(xscale.eq.0.) go to 55
      ys=(ymax-ymin)/yscale+1.
      xs=xx(2)/xscale+1.
      go to 65
55    xs=0.
      ys=0.
c
c    plot axes and titles
c
65    call crtplt(m,xx,yy,dum,dum,isym,iplot,xs,ys,nxl,nyl,
     -       titles,nt,x3,y3,ndec,idrive,iclose)
      nt=2
      iplot=-1
      isym=1
      do 66 i=1,max
66      if(y(i).lt.bval) go to 67
67    ib=i
c
c   plot values
c
75    m=0
      do 80 i=ib,max
        m=m+1
        yp(m)=y(i)
        xp(m)=x(i)
80      if(yp(m).ge.1.0e+30) go to 90
85    iclose=1
      call crtplt(m,xp,yp,dum,dum,isym,iplot,xs,ys,nxl,nyl,
     -       titles,nt,x3,y3,ndec,idrive,iclose)
      go to 99
90    m=m-1
      do 100 ii=i+1,max
100     if(y(ii).lt.1.0e+30) go to 110
      go to 85
110   ib=ii
      iclose=0
      call crtplt(m,xp,yp,dum,dum,isym,iplot,xs,ys,nxl,nyl,
     -       titles,nt,x3,y3,ndec,idrive,iclose)
      if(ib.eq.max) go to 99
      go to 75
c
99    return
      end
      subroutine numchr(var,nchr)
      character*40 var
      do 10 i=40,1,-1
      if(var(i:i).ne.' ')go to 50
10    continue
50    nchr=i
      return
      end
      subroutine maxmin(a,n,amin,amax)
      dimension a(1)
      amin=a(1)
      amax=amin
      do 1 i=1,n
      if(a(i).lt.amin)amin=a(i)
      if(a(i).gt.amax)amax=a(i)
    1 continue
      return
      end
      subroutine namemc(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     The program does not handle the repeat indicator *
c     nor does it handle subscripted array variables or
c     substrings of character variables.
c     It also does not handle complex or logical variables.
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv,vset
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     check for continuation of array values on succeding lines
c
      item=ichar(var(i))
      if(item.eq.43.or.item.eq.45.or.item
     & .eq.46.or.(item.gt.47.and.item.lt.58).or.var(i).eq.apos) then
      k=i
      go to 95
      endif
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      vset=.false.
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      vset=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
   95 if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if(var(l).eq.blank.and.nn.eq.0) go to 100
      if(var(l-1).eq.blank.and.var(l).eq.comma.and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
c
c     check for consecutive commas
c
      if(var(l).eq.comma.and.vset) then
      vset=.false.
      inum=inum+1
      go to 90
      else if(var(l).eq.comma) then
      inum=inum + 1
      go to 100
      endif
c
      if(var(l).eq.blank
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57.and.nn.eq.0) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
      vset=.true.
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
c
c     check for consecutive commas
c
      if(var(j).eq.comma.and.vset) then
      vset=.false.
      go to 105
      else if(var(j).eq.comma) then
      inum=inum + 1
      go to 105
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      vset=.true.
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      print*,'1756 ',l
      print'(a2)',var(l)
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 130
  120 m=l
  130 if(m.lt.81) then
      go to 15
      else
      go to 10
      endif
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
        subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c       assigns values to proper variable
c       variables are passed to program profgrd thru common blocks
c       numr=position in the array var where real variables start
c       numa=position in the array var where arrays start
c       nnvar=number of variables in program profgrd
c
        parameter (nnvar=12,numr=4)
        character*6 pvar,var(nnvar)
        character*56 tvar,kvar,cfmt
        logical chv
        character infil*50,binfil*50,pgm*8,gtitle*56,ascifil*50,
     1  coeffil*50,title*40
        integer flag
        common/prof/infil,binfil,title,nrow,dx,dy,ddxy,xscale,
     1  yscale,idrive,inn,pgm,idum,iout,xo,yo,iasci,flag,gtitle,
     2  ascifil
        common /parms/iseg,coeffil,seg(22)
c
        data var/'iseg','idrive','flag','ddxy','xscale','yscale',
     1  'title','infil','binfil','coeffi','ascifi','seg'/
        numa = 12
c
c     inum is set to one only to prevent a compilier
c     warning about a variable not being used.
c     it must be removed if any arrays would be used
c     in the future.
c
c      inum=1
        do 190 i=1,nnvar
        if (pvar .ne. var(i)) goto 190
        if (.not. chv) then
c
c       non character value
        m = 57 - nn
        im = m - 1
        kvar(m:56)=tvar(1:nn)
        if (i .lt. numr) then
c
c       integer value
c
        write (cfmt,50) im,nn
50      format ('(',i2,'x,i',i2,')')
        read (kvar,cfmt) jvar
        else
c
c       real value
c
        write (cfmt,60) im,nn
60      format('(',i2,'x,g',i2,'.0)')
        read (kvar,cfmt) xvar
        endif
        endif
        goto(101,102,103,104,105,106,107,108,109,110,111,112),i
101     iseg = jvar
        goto 200
102     idrive = jvar
        goto 200
103     flag = jvar
        goto 200
104     ddxy = xvar
        goto 200
105     xscale=xvar
        goto 200
106     yscale = xvar
        goto 200
107     title = tvar(1:nn)
        goto 200
108     infil = tvar(1:nn)
        goto 200
109     binfil = tvar(1:nn)
        goto 200
110     coeffil = tvar(1:nn)
        goto 200
111     ascifil = tvar(1:nn)
        goto 200
112     seg(inum) = xvar
        go to 200
190     continue
        write (*,*) ' error in namelist - ',pvar,
     1  ' variable not included'
        stop
200     nvar=i
        return
        end


        subroutine shiftr(a,ll)
c       added mar,1988 for pc
c       shift list directed real input to right of field for internal read
        character*20 a,c
        character*1 b(20)
        equivalence (c,b)
        c=a
        do 10 i=ll,1,-1
        if (b(i) .ne. ' ') goto 11
10      continue
        i=1
11      ishft=ll-i
        a=' '
        a(ishft+1:ll)=c(1:i)
        return
        end





      subroutine crtplt(m,x,y,xoff,yoff,isym,iplot,xs,ys,
     $ nxl,nyl,titles,nt,x3,y3,ndec,idrive,iclose)
c
c--tektronix-4010 crt general curve plots & labeling-- 11/2/77.
c  (also other driver options--see parm 'idrive')
c  by w.l.anderson, u.s.geological survey, denver, colorado.
c--routine uses a honeywell multics plot package (evenden,1977).
c
c--operational notes:
c  (1) after call crtplt(...), the user should program a 'wait'
c  or pause to obtain crt hardcopy (for idrive=1 and iclose=1)
c  (2) for idrive=0 (calcomp option), the user should make
c  sure the plot file is closed after last call (see parm iclose)
c
c
c--parameters:
c
c     m=        no.points given in x,y arrays..may be splined
c               before call
c               if m>0 points connected vectorally;
c               if m<0 iabs(m) points char. plotted via isym..
c     x=        abscissa array(dim.m). must be supplied and should
c               be in ascending order.
c               zero or neg.values ignored in logx options only.
c     y=        ordinate array(dim.m) corresponding to x array.
c               zero or neg.values ignored in logy options only.
c     xoff=     x-translation option. use 0.0 if x array
c               not to be translated.
c     yoff=     y-translation option. use 0.0 if y array
c               not to be translated.
c                       note: xoff,yoff are 'added' at
c                       plot time only and does not
c                       destroy the x,y arrays--xoff,yoff
c                       are in linear or powers-of-10 units
c                       depending on iplot parameter..
c                       when iplot>0,xoff=yoff=0. should be used.
c     isym=     software pen or char symbol selection.
c               (1<=isym<=8) e.g. isym=1 selects
c               solid line if m>0,..isym=2 selects
c               dashed line, etc.  when m<0 box,diamond
c               char. used...note: isym>8 ok if m<0...
c               see plot package symbol table, where
c               isym=table value + 1.
c     iplot=1   to initilize plot area & plot (x,y) curve.
c          =-1 to plot (x,y) curve over previous plot.
c                i.e., after initial iplot=1 was used.
c          =2 to initilize plot area & plot (logx,y).
c          =-2 to plot (logx,y) over previous plot.
c                i.e., after initial iplot=2 was used.
c          =3 to initilize plot area & plot (logx,logy).
c          =-3 to plot (logx,logy) over prev. plot.
c               i.e., after initial iplot=3 was used.
c          =4 to initilize plot area & plot (x,logy).
c          =-4 to plot (x,logy) over prev. plot.
c               i.e., after initial iplot=4 was used.
c               note: when iplot>0 the (x,y) data & other parms
c               are used to initilize the image area first. when
c               iplot<0 the (x,y) data is simply plotted
c               over last image & scale established--
c               thus points off-scale will be ignored.
c     xs=       xsize (in.): if xs=0., max.board size used.
c               if 0<xs<max.board, xs-1.' is x-neatline size.
c               xs>0.0 must be specified when idrive=0 or 6,
c               else a very large board size is assumed.
c               when iplot=3, xsize may be reduced automatically
c               to force equal log cycle size on both axes.
c     ys=       ysize (in.): if ys=0., max.board size used.
c               if 0<ys<max.board, ys-1.' is y-neatline size.
c               ys>0.0 must be specified when idrive=0 or 6,
c               else a very large board size is assumed.
c               when iplot=3, ysize may be reduced automatically
c               to force equal log cycle size on both axes.
c     nxl=      no.vert.line intervals..for x or logx axis
c               labels. nxl is fixed for iplot=1,4 (1<=nxl).
c               nxl is max.no.log cycle intervals desired
c               when iplot=2,3 (actual no.cycles should be .le. nxl)
c     nyl=      no.horiz.line intervals. for y or logy axis
c               labels.  nyl is fixed for iplot=1,2 cases.
c               nyl is max.no.log cycle intervals desired
c               when iplot=3,4 (actual no.cycles should be .le. nyl).
c                       note: in log axis option, no. cycles
c                       should be .le.10 (for crt)..the subr
c                       will adjust the no.cycles to actual range.
c                       the easiest thing to do for any log option
c                       is to use nxl or nyl = 10 and let routine
c                       determine the actual no. of log cycles....
c     titles(nt)=  optional title character arrays (if nt>0)
c                       (use dummy name if nt=0)
c               titles(1)= max.40 char.x-axis title (if nt>=1).
c               titles(2)= max.40 char.y-axis title (if nt>=2).
c               titles(3)= max.40 char. data line#1 title (if nt>=3),
c                       starting at relative position (x3,y3) in plot.
c                       (see parms x3,y3 below)
c               titles(4)= max.40 char. data line#2 title (if nt>=4),
c                       plotted under line#1 at (x3,y3-ydel)--ydel is
c                       determined by subroutine crtplt automatically..
c               titles(5)= max.40 char. data line#3 title (if nt=5),
c                       plotted under line#2 at (x3,y3-2*ydel)..
c               notes on titles:
c               (1). all titles arrays should be left-justified, and
c                    padded by blank char's to fill all 40-char's.
c               (2). for x and y-axis titles (nt>=2), the routine will
c                   automatically center titles axes.
c               (3). all blank char's could be used to suppress an
c                   intermediate title line..
c     nt=0      if no titles. (titles array name still required)
c               =1 if only x-axis title.
c               =2 if both x and y-axis titles.
c               =3 x & y axis titles + 1st data titles(3) at (x3,y3)
c               =4 same as =3 + 2nd data titles(4) at (x3,y3-ydel)
c               =5 same as =4 + 3rd data titles(5) at (x3,y3-2*ydel)
c               (ydel is determined by 'crtplt')
c     x3=       relative x-position of 1st char. of titles(3), where
c               0<=x3<=1.0 means relative position in (xmin,xmax);e.g.
c               x3=.5 is mid-point,x3=0 at xmin+xdel,x3=1.0 at xmax-xdel
c               where xdel depends on xrange & size.
c     y3=       relative y-position of 1st char. of titles(3),where
c               0<=y3<=1.0 means relative position in (ymin,ymax);e.g.
c               y3=.5 is mid-point,y3=0 at ymin+ydel,y3=1.0 at ymax-ydel
c               where ydel depends on yrange & size.
c     ndec=     no.dec.places in linear axis label values (usually 1 is
c               sufficient for most plots). 0<=ndec<8 may be used for
c               format (f08.ndec).  ndec<0 uses (1pe8.1).
c               ndec used only when iplot.ne.3 or -3.
c     idrive=   5  HP 7475
c               8  CGA monitor output
c               9  EGA monitor output
c              10  VGA monitor output
c
c     iclose=   0 to not close plot file via call endpt(ie).
c                 (use for overplot mode -- iplot<0 cases).
c               1 (or .ne.0) to close plot file via call endpt(ie).
c                 (use to terminate a series of overplots or
c                 to terminate a single iplot>0).
c                note:  if iplot>0 is used after 1st call, then an
c                 automatic close is performed with each new iplot>0.
c
c--subprograms called:
c         scale,line,char,pltset,endpt,neatl,vchar (in
c              Godson's plot.lib)
c--an error message is typed & return without action in cases where
c  some parameters or data are in error....
c
      dimension x(1),y(1),dxp(2),dyp(2),xp(4),yp(4)
      dimension p(200),tik(4),tik2(4)
      integer is(2),ie(2)
      character*40 titles(5),buf
      character*4 fmt(2),elab,fmt0(2),fmte(2)
c       fix isym type for vchar
        character*4 aisym
        equivalence(aisym,iisym)
      data tik/.06,.12,.06,.12/,tik2/.03,.06,.03,.06/
      data new/1/,nwth/8/,is/1,0/,ie/2,0/
      data xp/2*0.0,0.96,0.0/,yp/2*0.0,.96,0.0/
      data fmt0/'(f08','.1) '/,fmte/'(1pe','8.1)'/
      iplt=iabs(iplot)
      if(iplt.gt.4.or.iplt.eq.0.or.nt.gt.5
     1 .or. idrive.lt.5.or.idrive.eq.6.or.idrive.eq.7
     2 .or. idrive .gt. 10) goto 9999
c       jun 88 changed error check
c     &.or.((idrive.eq.0.or.idrive.eq.6).and.(xs.eq.0.0.or.ys.eq.0.0))
c     $.or.idrive.lt.0.or.idrive.gt.6) go to 9999
      indec=iabs(ndec)
      if(indec.gt.8) go to 9999
      nfmt=nwth
      isym1=isym-1
      if(isym1.lt.0) go to 9999
      im=iabs(m)
      if(iplot.lt.0) go to 1111
c--set up new plot area (iplot.gt.0) 
      if(new.eq.1) go to 1
      ie(2)=2
      call endpt(ie)
1     call pltset(idrive,xb,yb,is)   
34    new=0
      xp(4)=xb
      yp(4)=yb
      if(xs.gt.0.0.and.xs.lt.xb) xp(4)=xs
      if(ys.gt.0.0.and.ys.lt.yb) yp(4)=ys
      xp(1)=xp(4)-xp(3)-.04
      yp(1)=yp(4)-yp(3)-.04
      call maxmin(x,im,dxp(1),dxp(2))
      xmin=dxp(1)
      xmax=dxp(2)
      call maxmin(y,im,dyp(1),dyp(2))
      ymin=dyp(1)
      ymax=dyp(2)
      go to (402,402,10,401),iplot
402   continue
      go to (401,20,9999,9999),iplot
401   if(iplot.eq.4.and.ymax.le.0.0) go to 9999
      if(iplot.eq.4) go to 10
c--linear(x),linear(y)
      xp(2)=0.0
      yp(2)=0.0 
      call scale(dxp,dyp,xp,yp,4,ierr) 
      if(ierr.ne.0) go to 9999
      nx=2
      ixp=nx*nxl
      ndel=(dxp(2)-dxp(1))/ixp
       if (ndel.ne.0) go to 14
       xdel=(dxp(2)-dxp(1))/ixp
       n=0
8      xdel=xdel*10
       n=n+1
       if (xdel.lt.1.) go to 8
       ndel=xdel
       xdel=ndel/(10.**n)
       go to 5
14     xdel=ndel
5     ny=2
      iyp=ny*nyl
      ndel=(dyp(2)-dyp(1))/iyp
       if (ndel.ne.0) go to 16
       ydel=(dyp(2)-dyp(1))/iyp
       n=0
2      ydel=ydel*10
       n=n+1
       if (ydel.lt.1.) go to 2
       ndel=ydel
       ydel=ndel/(10.**n)
       go to 51
16     ydel=ndel
51    if(ndec.lt.0) go to 511
      fmt(1)=fmt0(1)
      if(nwth.gt.0) write(fmt(1),60)nwth
60    format('(f',i2)
      fmt(2)=fmt0(2)
      if(indec.ne.1) write(fmt(2),61)indec
61    format('.',i1,')')
      go to 512
511   fmt(1)=fmte(1)
      fmt(2)=fmte(2)
c--label linear(x) and/or linear(y) axes
512   if(iplt.eq.1.or.iplt.eq.4) 
     & call xaxis(dxp,dyp,xp,xdel,nx,.08,fmt,nfmt)
      if(iplt.eq.1.or.iplt.eq.2)
     & call yaxis(dyp,dxp,yp,ydel,ny,.08,fmt,nfmt)
      call neatl
100   if(nt.le.0) go to 1111
c--label plot as selected by parms nt & titles:
      call numchr(titles(1),nchr)
36    go to (41,42,42,41),iplot
41    xx=.5*(dxp(1)+dxp(2))
      go to 44
42    xx=10.**(.5*(ix+jx))
44    buf=titles(1)
       call vchar(xx,dyp(1),buf,
     $ nchr,2,.1,0.,-.05*nchr,-.75)
      if(nt.le.1) go to 1111
      call numchr(titles(2),nchr)
38    go to (50,50,52,52),iplot
50    yy=.5*(dyp(1)+dyp(2))
      go to 55
52    yy=10.**(.5*(iy+jy))
55    buf=titles(2)
       call vchar(dxp(1),yy,buf,
     $ nchr,2,.1,1.5707963,-.05*nchr,.75)
      if(nt.le.2) go to 1111
      go to (71,72,72,71),iplot
71    xrange=dxp(2)-dxp(1)
      xdel=.24*xrange/xp(1)
      if(x3.lt.1.0.and.x3.gt.0.0) xx=dxp(1)+x3*xrange
      if(x3.le.0.0) xx=dxp(1)+xdel
      if(x3.ge.1.0) xx=dxp(2)-xdel
      go to (75,9999,9999,76),iplot
72    xrange=jx-ix
      xdel=10.**(.24*xrange/xp(1))
      if(x3.lt.1.0.and.x3.gt.0.0) xx=10.**(ix+x3*xrange)
      if(x3.le.0.0) xx=xdel*10.**ix
      if(x3.ge.1.0) xx=(10.**jx)/xdel
      go to (9999,75,76,9999),iplot
75    yrange=dyp(2)-dyp(1)
      ydel=.24*yrange/yp(1)
      if(y3.lt.1.0.and.y3.gt.0.0) yy=dyp(1)+y3*yrange
      if(y3.le.0.0) yy=dyp(1)+ydel
      if(y3.ge.1.0) yy=dyp(2)-ydel
      yy1=yy-ydel
      yy2=yy1-ydel
      go to (80,80,9999,9999),iplot
76    yrange=jy-iy
      ydel=10.**(.24*yrange/yp(1))
      if(y3.lt.1.0.and.y3.gt.0.0) yy=10.**(iy+y3*yrange)
      if(y3.le.0.0) yy=ydel*10.**iy
      if(y3.ge.1.0) yy=(10.**jy)/ydel
      yy1=yy/ydel
      yy2=yy1/ydel
80    call numchr(titles(3),nchr)
      buf=titles(3)
       call vchar(xx,yy,buf,nchr,2,.12,0.0,0.0,0.0)
      if(nt.le.3) go to 1111
      call numchr(titles(4),nchr)
      buf=titles(4)
       call vchar(xx,yy1,buf,nchr,2,.12,0.0,0.0,0.0)
      if(nt.le.4) go to 1111
      call numchr(titles(5),nchr)
      buf=titles(5)                    
       call vchar(xx,yy2,buf,nchr,2,.12,0.0,0.0,0.0)
      go to 1111
10    if(ymin.gt.0.0) go to 20
      if(ymax.le.0.0) go to 9999
      y1=y(1)
      do 11 i=2,im
      y2=y(i)
      if(sign(1.0,y1).ne.sign(1.0,y2)) go to 12
11    y1=y2
      go to 9999
12    if(sign(1.0,y1).eq.1.0) ymin=y1
      if(sign(1.0,y2).eq.1.0) ymin=y2
20    if(iplot.eq.4) go to 302
      if(xmin.gt.0.0) go to 30
      if(xmax.le.0.0) go to 9999
      do 21 i=2,im
      xmin=x(i)
      if(xmin.gt.0.0) go to 30
21    continue
c--set up log(x)
30    ix=alog10(xmin)
      if(xmin.lt.10.**ix) ix=ix-1
      dxp(1)=10.**ix
      jx=alog10(xmax)
      if(xmax.gt.10.**jx) jx=jx+1
      if(jx-ix.ge.nxl) jx=ix+nxl
      dxp(2)=10.**jx
      xp(2)=8.0
      yp(2)=0.0
      if(iplot.eq.2) go to 40
c--set up  log(y)
302   iy=alog10(ymin)
      if(ymin.lt.10.**iy) iy=iy-1
      dyp(1)=10.**iy
      jy=alog10(ymax)
      if(ymax.gt.10.**jy) jy=jy+1
      if(jy-iy.ge.nyl) jy=iy+nyl
      dyp(2)=10.**jy
      xp(2)=0.0
      if(iplot.eq.3) xp(2)=8.0
      yp(2)=8.0
c--force equal log cycle size when iplot=3
      nxc=jx-ix
      nyc=jy-iy
      xsc=xp(1)/nxc
      ysc=yp(1)/nyc
      if(xsc.eq.ysc) go to 40
      if(xsc.lt.ysc) go to 500
      xp(1)=nxc*ysc
      xp(4)=xp(1)+1.
      go to 40
500   yp(1)=nyc*xsc
      yp(4)=yp(1)+1.
40    call scale(dxp,dyp,xp,yp,4,ierr)
      if(ierr.ne.0) go to 9999
      if(iplot.eq.4) go to 1101
c--tick horiz. log(x) lines
      xl=dxp(1)
      ixc=ix+1
      xc=10.**ixc
      xd=xl
      np=0
      k=1
105   xx=xl+k*xd
1051  if(k.eq.9) go to 106
      np=np+1
      if(np.gt.200) go to 9999
      p(np)=xx
      k=k+1
      go to 105
106   xl=xc
      ixc=ixc+1
      xc=10.**ixc
      if(xc.gt.dxp(2))  go to 107
      xd=xl
      k=0
      xx=xl
      go to 1051
107   yc=dyp(2)
      do 108 k=1,np
      n=1
      if(mod(k,9).eq.0) n=2
108   call vchar(p(k),yc,'l',1,1,tik(n),0.0,0.0,-tik2(n))
      yc=dyp(1)
      do 109 k=1,np
      n=1
      if(mod(k,9).eq.0) n=2
109   call vchar(p(k),yc,'l',1,1,tik(n),0.0,0.0,tik2(n))
c--label log(x) axis
      yc=dyp(1)
      do 110 nx=ix,jx
      xc=10.**nx
      call vchar(xc,yc,'10',2,2,.08,0.0,-.24,-.2)
      write(elab,4)nx
4     format(' ',i3)
110   call vchar(xc,yc,elab,4,2,.08,0.0,-.24,-.1)
      if(iplot.eq.2)  go to 5
c--tick vert. log(y) lines
1101  yl=dyp(1)
      iyc=iy+1
      yc=10.**iyc
      yd=yl
      np=0
      k=1
115   yy=yl+k*yd
1151  if(k.eq.9) go to 116
      np=np+1
      if(np.gt.200) go to 9999
      p(np)=yy
      k=k+1
      go to 115
116   yl=yc
      iyc=iyc+1
      yc=10.**iyc
      if(yc.gt.dyp(2)) go to 117
      yd=yl
      k=0
      yy=yl
      go to 1151
117   xc=dxp(1)
      do 118 k=1,np
      n=3
      if(mod(k,9).eq.0) n=4
118   call vchar(xc,p(k),'-',1,1,tik(n),0.0,tik2(n),0.0)
      xc=dxp(2)
      do 119 k=1,np
      n=3
      if(mod(k,9).eq.0) n=4
119   call vchar(xc,p(k),'-',1,1,tik(n),0.0,-tik2(n),0.0)
c--label log(y) axis
      xc=dxp(1)
      do 121 ny=iy,jy
      yc=10.**ny
      call vchar(xc,yc,'10',2,2,.08,0.0,-.32,-.1)
      write(elab,4)ny
121   call vchar(xc,yc,elab,4,2,.08,0.0,-.32,0.0)
      if(iplot.eq.4) go to 122
      call neatl
      go to 100
122   nx=2
      ixp=nx*nxl
      xdel=(dxp(2)-dxp(1))/ixp
      go to 51
c--plot requested curve
1111  if(m.lt.0) go to 300
      if(xoff.ne.0.0.or.yoff.ne.0.0) go to 111
      if(idrive.eq.0.or.idrive.eq.6) go to 1002
      ip=0
      do 1000 i=1,m
      xx=x(i)
      if(xx.lt.dxp(1)) go to 1001
      yy=y(i)
      if(yy.lt.dyp(1)) go to 1001
      call line(xx,yy,1,ip,isym1)
      ip=1
      go to 1000
1001  ip=0
1000  continue
      go to 99
1002  call line(x,y,m,0,isym1)
99    call vchar(.05,.48,'.',1,3,.08,0.0,0.0,0.0)
      if(iclose.eq.0) return
      ie(2)=2
      call endpt(ie)
      new=1
      if(ie(2).ne.0) write(6,66) ie(2)
66    format(' endpt: ',i10,' plotting errors.')
      return
111   ip=0
      x10=10.**xoff
      y10=10.**yoff
      do 200 i=1,m
      if(iplt.eq.1.or.iplt.eq.4) xx=x(i)+xoff
      if(iplt.eq.2.or.iplt.eq.3) xx=x(i)*x10
      if(xx.lt.dxp(1).or.xx.gt.dxp(2)) go to 201
      if(iplt.le.2) yy=y(i)+yoff
      if(iplt.ge.3) yy=y(i)*y10
      if(yy.lt.dyp(1).or.yy.gt.dyp(2)) go to 201
      call line(xx,yy,1,ip,isym1)
      ip=1
      go to 200
201   ip=0
200   continue
      go to 99
c--plot 'isym' points (m.lt.0)
300   if(xoff.ne.0.0.or.yoff.ne.0.0) go to 333
      do 301 i=1,im
      if(x(i).lt.dxp(1).or.y(i).lt.dyp(1)) go to 301
        iisym=isym
      call vchar(x(i),y(i),aisym,1,1,.08,0.0,0.0,0.0)
301   continue
      go to 99
333   x10=10.**xoff
      y10=10.**yoff
      do 400 i=1,im
      if(iplt.eq.1.or.iplt.eq.4) xx=x(i)+xoff
      if(iplt.eq.2.or.iplt.eq.3) xx=x(i)*x10
      if(xx.lt.dxp(1).or.xx.gt.dxp(2)) go to 400
      if(iplt.le.2) yy=y(i)+yoff
      if(iplt.ge.3) yy=y(i)*y10
      if(yy.lt.dyp(1).or.yy.gt.dyp(2)) go to 400
        iisym=isym
      call vchar(xx,yy,aisym,1,1,.08,0.0,0.0,0.0)
400   continue
      go to 99
9999  write(6,69999)
69999 format('0error--in call crtplt')
      return
      end





