cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      program phist
c  This program plots a crude histogram on the terminal from a post file
c     Uses guts of 'this.fortran' by Tien Grauch.
c
      parameter (maxpst=52, dval=1.e38)
      character id*8
      character*50 afile
      character*1 iex(50)
      dimension iclass(50), post(8)
c      integer istack(1024)
        logical    binary, expres
c      parameter (dval=0.170412E+39)

      amax = -dval
      amin = +dval
      kmax=0
      ihistmx=0
      do 5 i=1,50
      iclass(i)=0
      iex(i)='X'
5     continue

c      call getfile2(10,' Infile ','in','unformatted','.pos')

      print '(a\)',' Input binary post or xyz file name: '
      read(5,100) afile
  100 format(a)
        call openf(10,afile,istack,maxpst,binary,expres,
     1        nword,maxlog,ntotal)
c  set function switches
        ifunc = 0
	if ( binary  .and.  nword .gt. 3 ) then
	  ifunc = 2
	endif
	if ( binary  .and.  nword .eq. 3 ) then
	  ifunc = 3
	endif
        if ( expres ) ifunc = 1
      if(ifunc.eq.2) then
        print '(a\)',' Enter z-field (1 to 6):'
        read(5,*) nopt
      endif

c      call getint (' Give post option (1 to 6)', 1, nopt)
      nopt=nopt+2
c
c  read file, calculate mean
c
      icount=0
      sum=0.
 10   if(ifunc.eq.2) then
        call getpost (10, id, post, *20)
        aval=post(nopt)
      else if(ifunc.eq.3) then
        read(10,end=20) x,y,aval
      endif
      icount=icount+1
      if(aval.gt.amax) amax=aval
      if(aval.lt.amin) amin=aval
      sum=sum+aval
      go to 10
 20   rewind 10
c
c  determine interval
c
      amean=sum/icount
      div=.1e0
550   del1=amax-amean
      del2=amean-amin
      del1=div*del1
      del2=div*del2
      if((amax-amin)/del1.le.50.and.del1.lt.del2) go to 555
      if((amax-amin)/del2.le.50.and.del2.lt.del1) go to 556
      div=div+.1e0
      go to 550
555   del=del1
      go to 9
556   del=del2
c
c  count number of points in each interval, calc variance and line density
c
9     k=0
      lden=0
      s2=0.
 22   if(ifunc.eq.2) then
        call getpost (10, id, post, *32)
        aval=post(nopt)
      else if(ifunc.eq.3) then
        read(10,end=32) x,y,aval
      endif
      line=k
      top=amin-(del*.5)
      k=0
      s2=(aval-amean)**2.+s2
25    top=top+del
      k=k+1
      if(aval.ge.top) go to 25
      iclass(k)=iclass(k)+1
      if(k.gt.kmax) kmax=k
      if(k.eq.line.or.i.eq.1) go to 30
      line=k
      lden=lden+1
30    continue
      go to 22
 32   close(10)
      s2=s2/(icount-1)
      s=sqrt(s2)
      do 35 i=1,kmax
      if(iclass(i).gt.ihistmx) ihistmx=iclass(i)
35    continue
c Addition: to print out a phist.tmp file.
      afile='phist.tmp'
      open(11,file=afile,form='formatted',status='unknown')
      if(ihistmx.le.50) go to 50
      fac=50./ihistmx
      do 40 i=1,kmax
      iclass(i)=iclass(i)*fac
40    continue
c
c  plot classes, print parameters
c
50    print 802,amin,amax,amean,s,icount
      write(11,802)amin,amax,amean,s,icount
802   format(/,' min=    ',f14.7,10x,'max=',f14.7,/,' mean=   ',f14.7,
     1   10x,'std dev=',f14.7,/i12,' valid points,',3x,f6.2)
      print 803,del
      write(11,803)del
803   format(//,' midpoint ',1x,'(interval=',f14.7,')',1x,'HISTOGRAM')
      print*,'<Enter> to continue'
      read(5,100) afile
      pt=amin-del
      do 60 i=1,kmax
      pt=pt+del
      if(iclass(i).eq.0) go to 59
      print 804,pt,(iex(j),j=1,iclass(i))
      write(11,804)pt,(iex(j),j=1,iclass(i))
804   format(1x,g10.4,1x,'|',50(a1))
      go to 60
59    print 804,pt
      write(11,804)pt
60    continue
      read(5,100) afile
      print*,'Output also contained in file PHIST.TMP.'
      go to 300
c
c  ERROR message
c
 300  stop ' '
      end
c*******************************************************************************
      subroutine openf( iswt, ifile, ibuff, maxpst, binary, expres,
     1                 nword, maxlog, ntotal )

c  open data file, derive and checkout file structure.
c
c  input file types:
c  xyz      single record xyz data.
c  post     single record containing an id (8 characters), x, y, and
c           1-n dependant-variable channels. Normally n=6, but this
c           program accepts more channels.
c
c  Note: express file options deleted for PC version 1 series.

c        parameter  ( maxbuf=1024 )
c      integer*2  ibuff(maxbuf)
	character  ifile*(*)
	logical    binary, expres
	nword  = 0
	maxlog = 0
	ntotal = 0
	binary = .false.
	expres = .false.

c  open binary
	open( unit=iswt, file=ifile, status='old', form='unformatted',
     1     mode='read', iostat=ios )
	if ( ios .ne. 0 ) then
        print *,' error opening input data file ( cannot find )'
	  stop
	endif

c  binary or ascii input ?
	read( iswt, err=10, iostat=ios ) itest
10	rewind( iswt )
	if ( ios .eq. 0 ) then
	  binary = .true. 
	  else
	  close( iswt ) 
	  go to 100
	endif

c  file is binary
c  lenrec returns zero if buffer is too short
c      nword = lenrec( iswt, ibuff, maxpst )
      nword = lenrec( iswt, maxpst )
	if ( nword .eq. 1  .or.  nword .eq. 2 ) go to 200

c  simple xyz file.
	if ( nword .eq. 3 ) then
        print *, ' binary xyz data'
	  go to 999
	endif

c  simple post file.
	if ( nword .ge. 4   .and.   nword .le. maxpst ) then
        nchan  = nword - 2
        nnchan=nchan-2
        print *, ' post data,  found', nnchan, ' data channels'
	  go to 999
	endif


c  open for ascii input
100	  open( unit=iswt, file=ifile, status='old', form='formatted', 
     1       mode='read' )
	go to 999

c  record length undefined.
200   print *, ' input record length =', nword
      print *, ' input is neither xyz or post'

999	return
        end
c******************************************************************************
c        function lenrec( iswt, irec, nwork )
        function lenrec( iswt, nwork )
c  brute force recovery of record length
        integer*4 irec(52)

	istep = 10
	if ( nwork .gt. 100 ) istep = 20
	if ( nwork .gt. 500 ) istep = 50
	rewind iswt

	do max = istep, nwork, istep
	  read( iswt, err=10 ) ( irec(i), i=1, max )
	  rewind iswt
	enddo
	lenrec = -1
	go to 99

10	rewind iswt
	do lenrec = max, max - istep, -1
	  read( iswt, err=20 ) ( irec(i), i=1, lenrec )
	  go to 99
20	  rewind iswt
	enddo

99	rewind iswt
	return
        end

