c PROFPLOT.FOR
c
c  To plot a profile.
c  Modified from PROFSPEC by Lin Cordell
c
      character*50 infile
      character*60 ititle
      character*8 jfmt,jfmt1
      character answer*1
      dimension isl(5),dxp(2),dyp(2),xp(4),yp(4),ie(10)
      dimension x(20000),f(20000),id(14),pgm(2),p(6)
      integer kfmt(2),kfmt1(2),jtitle(15),isym(20000)
      equivalence (jfmt,kfmt(1)),(jfmt1,kfmt1(1)),(ititle,jtitle(1))
      data jfmt/'(f8.2)'/,jfmt1/'(g14.6)'/
      i=0
102   print*,' Enter format: 1= ascii xyz; 2=binary xyz; 3=ascii xz;'
      print 101,' 4=binary xz; 5=post; 6=one-row grid:'
      read*,mode
      if(mode.lt.1.or.mode.gt.6)then
        print*,' Got to be 1,2,3,4,5 or 6.'
        go to 102
      endif
      print 101,' Enter input file name:'
101   format(a,/' *'$)
      read(5,100)infile
      print 101,' Enter plotsys symbol code (4=x,6=+,etc):'
      read(5,*) jsym
      print 101,' Enter vertical offset:'
      read(5,*) offset
100   format(a)
      go to(161,162,163,164,165,166),mode
161   open(10,file=infile,form='formatted',status='old',err=298)
c      i=0
170   i=i+1
      if(i.gt.20000)go to 181
      read(10,*,end=171)x(i),zilch,f(i)
      f(i)=f(i)+offset
      isym(i)=jsym
      go to 170
171   ncol=i-1
      go to 180
162   open(10,file=infile,form='unformatted',status='old',err=298)
c      i=0
172   i=i+1
      if(i.gt.20000)go to 181
      read(10,end=173)x(i),zilch,f(i)
      f(i)=f(i)+offset
      isym(i)=jsym
      go to 172
173   ncol=i-1
      go to 180
163   open(10,file=infile,form='formatted',status='old',err=298)
c      i=0
174   i=i+1
      if(i.gt.20000)go to 181
      read(10,*,end=175)x(i),f(i)
      f(i)=f(i)+offset
      isym(i)=jsym
      go to 174
175   ncol=i-1
      go to 180
164   open(10,file=infile,form='unformatted',status='old',err=298)
c      i=0
176   i=i+1
      if(i.gt.20000)go to 181
      read(10,end=177)x(i),f(i)
      f(i)=f(i)+offset
      isym(i)=jsym
      go to 176
177   ncol=i-1
      go to 180
165   open(10,file=infile,form='unformatted',status='old',err=298)
      print 101,' Enter channel no (1,2,3,4,5 or 6):'
      read*,ichan
c      i=0
178   i=i+1
      if(i.gt.20000)go to 181
      read(10,end=179)pgm(1),pgm(2),xc,yc,(p(k),k=1,6)
      x(i)=xc
      f(i)=p(ichan)
      f(i)=f(i)+offset
      isym(i)=jsym
      go to 178
179   ncol=i-1
180   print 101,' Plot additional profiles?'
      read(5,100)answer
      if(answer.eq.'y'.or.answer.eq.'Y') go to 102
      xmax=-1.e30
      xmin=1.e30
      do 183 j=1,ncol
      if(x(j).gt.xmax)xmax=x(j)
183   if(x(j).lt.xmin)xmin=x(j)
      x0=xmin
      dx=(xmax-xmin)/(ncol-1)
      go to 184
166   open(10,file=infile,form='unformatted',status='old',err=298)
      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      if(ncol.le.20000.and.dx.ne.0.e0)go to 3
181   write(6,2)
2     format(1x,'Can''t handle it.')
      go to 999
3     continue
      read(10)dum,(f(j),j=1,ncol)
      do 182 i=1,ncol
182   x(j)=x0+(j-1)*dx
      xmin=x0
      xmax=x(ncol)
184   continue
      icode=0
      isl(1)=2
      iplotr=9
      xp(3)=1.5e0
      yp(3)=2.0e0
      xp(2)=0.
      print *,' Want program to find min & max for vertical scale? '
      read 888,answer
888   format(a1)
      if(answer.eq.'y'.or.answer.eq.'Y') goto 555
      print 101,' Enter min & max for vertical scale'
      read*,dyp(1),dyp(2)
      goto 666
555   fmax=-1.e30
      fmin=1.e30
c  find max and mins
      do 150 j=1,ncol
      if(f(j).gt.1.e20)go to 150
      if(f(j).gt.fmax)fmax=f(j)
      if(f(j).lt.fmin)fmin=f(j)
150   continue
      if(fmin.gt.0) then
        dyp(1)=fmin*.9
      else
        dyp(1)=fmin*1.1
      endif
      if(fmax.gt.0) then
        dyp(2)=fmax*1.1
      else
        dyp(2)=fmax*.9
      endif
666   continue
ccc      dxp(2)=1.1*.5/dx
ccc      dxp(1)=0.
      dxp(1)=xmin
      dxp(2)=xmax
      call fixscale(dxp(2),dxp(1),xdel)
10    call fixscale(dyp(2),dyp(1),ydel)
      print 815, dxp(1),dxp(2),dyp(1),dyp(2)
815   format(' x values:  minimum = ',g10.4,' maximum = ',g10.4,
     1 /' y values:  minimum = ',g10.4,' maximum = ',g10.4)
      print*,'Enter title (max 60 chars)'
      read 805, ititle
805   format(a)
      yp(2)=0.0
15    print 811
811   format(' Plot on screen(cga,ega,vga) or hp plotter(hp): '$)
      read 806, answer
806   format(a1)
      if(answer.ne.'h'.and.answer.ne.'H') go to 25
      iplotr=5
      call pltset(iplotr,xboard,yboard,isl)
      print 810,xboard-1.5,yboard-2.
810   format(' Max size available (in cm) is ',g10.4,' X ',g10.4)
      print*,' Enter vertical axis length in cm'
      read*,yp(1)
20    print*,' Enter horiz axis length in cm'
      read*,xp(1)
      go to 30
25    if(answer.eq.'c'.or.answer.eq.'C') then
      iplotr=8
      else if(answer.eq.'e'.or.answer.eq.'E') then
      iplotr=9
      else if(answer.eq.'v'.or.answer.eq.'V') then
      iplotr=10
      else
      go to 15
      endif
      call pltset(iplotr,xboard,yboard,isl)
      xp(1)=xboard-xp(3)-1
      yp(1)=yboard-yp(3)-1
30    nopts=4
      xp(4)=xp(1)+xp(3)+1
      yp(4)=yp(1)+yp(3)+1
      call scale(dxp,dyp,xp,yp,nopts,ier)
      do 40 j=1,ncol
      call vchar(x(j),f(j),isym(j),1,1,.16,0.,0.,0.)
40    continue
50    call xaxis(dxp,dyp,xp,xdel,5,.2,kfmt,8)
      if (yp(2).ne.0.0) then
            call yaxsl(dyp,dxp,yp,.2,kfmt1,14)
      else
            call yaxis(dyp,dxp,yp,ydel,5,.2,kfmt,8)
      endif
60    xd=3.5
      call vchar(xd,.5,jtitle,60,3,.25,0.,0.,0.)
      call neatl
      ie(2)=2
      call endpt(ie)
      go to 300
298   print*,'Input file not found or not binary.'
999   continue
300   close(10)
       stop
       end
      subroutine fixscale(xmax,xmin,del)
c
c  For plot_sys programs to find parameter del for call to xaxis or
c  yaxis knowing end points of axis
c
      power=-1.0
      constant=.5
      idel=alog10(xmax-xmin)
      tmp=float(1-idel)
      ten=10.e0**tmp
      quantity=ten*(xmax-xmin)
      if(quantity.gt.30.0) constant=1.0
      if(quantity.lt.15.0) constant=.2
      del=ten**power
      del=del*constant
      return
      end
