c  GHIST.FOR
c
c  This program plots a crude histogram on the terminal to aid in picking
c  contour intervals.  Also gives minimum, maximum and standard deviation.
c  nz must = 1, max ncols = 5000
c
      include 'fgraph.fi'
      include 'fgraph.fd'
      integer*2 dum2,fgd,fgd1,row
      integer*4 dum1,bgd,bgd1
      character*80 string

      character infile*50,id*56,pgm*8,answer*1,iexs*50,fmt*40
      character*1 iex(50)
      dimension dell(2),range(2)
      common/block/a(5000),iclass(200),iclasf(200)
      equivalence (iexs,iex(1))
      record/rccoord/cp
      character*8 buff,passwrd
      data amax/-1.0e37/,amin/1.0e37/,dcval/1.0e38/
      data bgd/2/,fgd/15/,bgd1/1/,fgd1/14/
      data passwrd/'sB8yr0hF'/
      call getarg(1,buff,length)
      if(buff.ne.passwrd) then
        print *,'program can only be run from program GEOCON'
        stop
      endif


      row=1
      call settextposition(row,1,cp)
6     print 801
801   format(1x,'Enter input file: ',$)
      read(5,'(a)')infile
      call opngrd(10,infile,ierr,*6,*300)
      read(10) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      if(nz.ne.1.or.ncol.gt.5000) go to 299
c
c  read file, find min, max, and calculate mean 
c
      icount=0
      sum=0
      do 10 j=1,nrow
      call row1io(a,ncol,10,1)
      do 8 i=1,ncol
      if(a(i).ge.dcval) go to 8
      if(a(i).gt.amax) amax=a(i)
      if(a(i).lt.amin) amin=a(i)
      sum=sum+a(i)
      icount=icount+1
8     continue
10    continue
c  calculate and print parameters
      if(icount.eq.0) go to 110
      en=1./float(icount)
      amean=sum*en
      percent=float(icount)/float(ncol*nrow)
      percent=100.e0*percent
      if((amax-amin).eq.0.0) go to 100
	diff=amax-amin
      row=row+2
      call settextposition(row,1,cp)
      print 802,amin,amax,amean,diff,icount,percent
802   format(/,1x,'minimum=',g14.7,10x,'maximum=',g14.7,/,1x,'mean=',
     &g14.7,10x,'max - min =',g14.7,/1x,i12,' valid points,',3x,f6.2,
     &'% of grid')
      row=row+4
      call settextposition(row,1,cp)
      print 806
806   format(/,1x,'Want to see histogram, standard dev., both, or ',
     1 'none [h,s,b,n]? ',\)
      read(5,901) answer
	iflag=3
      if(answer.eq.'h'.or.answer.eq.'H') iflag=1
      if(answer.eq.'b'.or.answer.eq.'B'.or.answer.eq.'y'.or.answer.
     1eq.'Y') iflag=1
      if(answer.eq.'s'.or.answer.eq.'S') iflag=2
	go to (699,9,300),iflag
c
c find default interval
c
699   lcount=0
      div=.1e0
c round off amin & amax & rename as range(1) & range(2)
      range(1)=amin
      range(2)=amax
      do 20 i=1,2
      if(range(i).eq.0.0) go to 20
      ntest=2-ifix(alog10(abs(range(i))))
      test=10.e0**float(ntest)
      rtest=float(nint(range(i)*test))
      range(i)=rtest/test
20    continue
550   lcount=lcount+1
      dell(1)=range(2)-amean
      dell(2)=amean-range(1)
      do 15 i=1,2
      dell(i)=div*dell(i)
      ntest=2-ifix(alog10(dell(i)))
      test=10.e0**float(ntest)
      dtest=float(nint(dell(i)*test))
      dell(i)=dtest/test
15    continue
      atest=range(2) - range(1)
      if(atest/dell(1).le.20.and.dell(1).lt.dell(2)) go to 555
      if(atest/dell(2).le.20.and.dell(2).lt.dell(1)) go to 556
      if(lcount.gt.20) then
	del=dell(1)
        if(dell(1).gt.dell(2)) del=dell(2)
	go to 552
      endif
      div=div+.1e0
      go to 550
555   del=dell(1)
      go to 552
556   del=dell(2)
552   icheck=nint((range(2)-range(1))/del) + 1
      row=row+2
      call settextposition(row,1,cp)
	print 999,del,icheck
999	format(/' I found histogram interval of ',g14.7,', giving '
     1 ,i3,' histogram classes')
      row=row+2
      call settextposition(row,1,cp)
	write(6,700)
700	format(1x,'Want to redefine interval ? [n] '\)
	read(5,901)answer
	if(answer.ne.'y'.and.answer.ne.'Y')go to 705
      row=row+2
      call settextposition(row,1,cp)
	write(6,701)
701	format(1x,'Enter new interval: '\)
	read(5,*)del
      icheck=nint((range(2)-range(1))/del) + 1
      row=row+1
      call settextposition(row,1,cp)
      if(icheck.gt.200.) then
         del=(range(2)-range(1))/200.
         print*,'No. of classes will exceed 200.  Changing del to ',del
      else
	 print 998,del,icheck
998	format(' Interval of ',g14.7,' gives ',i3,' histogram classes')
      endif
901   format(a1)
705   row=row+2
      call settextposition(row,1,cp)
	write(6,703)
703	format(1x,'Want to see percent or count of values [p or c]? ',
     1 \)
	read(5,901)answer
	if(answer.ne.'p'.and.answer.ne.'c'.and.answer.ne.'P'.and.
     1answer.ne.'C')go to 705
	newopt=1
	if(answer.eq.'c'.or.answer.eq.'C')newopt=2
c
c  count number of points in each interval, calc variance
c
9     k=0
      do 5 i=1,50
      iclass(i)=0
      iex(i)='X'
5     continue
      sum2=0.
      kmax=0
      ihistmx=0
      rewind 10
      read(10) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      do 31 j=1,nrow
      call row1io(a,ncol,10,1)
      line=k
      do 30 i=1,ncol
      if(a(i).ge.dcval) go to 30
      aa=a(i)-amean
      sum2=sum2+aa*aa
      top=range(1)-(del*.5)
      k=0
25    top=top+del
      k=k+1
      if(a(i).ge.top.and.iflag.eq.1) go to 25
      iclass(k)=iclass(k)+1
      if(k.gt.kmax) kmax=k
      if(k.eq.line.or.i.eq.1.or.iflag.eq.2) go to 30
      line=k
30    continue
31    continue
      s=sqrt(sum2/float(icount-1))      
	if(iflag.eq.2) go to 441
      do 35 i=1,kmax
      if(iclass(i).gt.ihistmx) ihistmx=iclass(i)
35    continue
      fac=50./float(ihistmx)
      if(ihistmx.le.50) fac=1.
      do 40 i=1,kmax
      iclasf(i)=ifix(float(iclass(i))*fac)
40    continue
c
c  print std dev
c
441   row=row+2
      call settextposition(row,1,cp)
      print 993,s
993	format(/' Standard Deviation = ',g14.7)
	if(iflag.eq.2) then
		print*
		go to 300
	endif
c
c  plot classes
c
41    row=row+2
      call settextposition(row,1,cp)
      print 803,del
803   format(/,' midpoint',1x,'(interval=',g14.4,')',3x,'HISTOGRAM',/)
      dum1= setbkcolor(bgd)
      dum2=settextcolor(fgd)
      pt=range(1)-del
      row=row+2
      do 60 i=1,kmax
      row=row+1
      if(row.ge.26) then
        dum1= setbkcolor(bgd1)
        call scrolltextwindow(1)
        dum1= setbkcolor(bgd)
      endif
       call settextposition(row,1,cp)
	label=iclass(i)
	if(newopt.eq.1)rlabel=(float(label)*100.0)/float(icount)
      pt=pt+del
      if(iclasf(i).eq.0) go to 59
	if(newopt.eq.2)then
         write(fmt,804)iclasf(i)
         write(string,fmt)pt,iexs,label
         len=index(string,'#') - 1
         call outtext(string(1:len))
       endif
cprint 804,pt,iexs,label
	if(newopt.eq.1)then
         write(fmt,714)iclasf(i)
         write(string,fmt)pt,iexs,rlabel
         len=index(string,'%')
         call outtext(string(1:len))
        endif
cprint 714,pt,iexs,rlabel
714     format('(1x,g10.4,1x,''|'',a',i6,',1x,f6.1,''%'')')
c714	format(1x,g10.4,1x,'|',a<iclassf(i)>,1x,f6.1,'%')
804     format('(1x,g10.4,1x,''|'',a',i6,',1x,i6,''#'')')
c804   format(1x,g10.4,1x,'|',a<iclassf(i)>,1x,i6)
      go to 60
59    if(iclass(i).eq.0) go to 599
c	if(newopt.eq.2)print 805,pt,label
c	if(newopt.eq.1)print 715,pt,rlabel
       if(newopt.eq.2) then
         write(string,805) pt,label
         call outtext(string(1:20))
       endif
       if(newopt.eq.1) then
         write(string,715) pt,rlabel
         call outtext(string(1:19))
       endif
715	format(1x,g10.4,1x,'|',1x,f6.1,'%')
805   format(1x,g10.4,1x,'|',1x,i6)
      go to 60
c599   print 805,pt
599    write(string,805) pt
       call outtext(string(1:11))
60    continue
      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)
      print 812
812   format(' REPRINT info, CHANGE histogram, or STOP ',
     1'[r, c, or s]? '\)
	read(5,901) answer
	if(answer.eq.'r'.or.answer.eq.'R') go to 441
	if(answer.eq.'c'.or.answer.eq.'C') go to 699
      go to 400
100   print 807,icount,amax,percent
807   format(1x,'All ',i12,' valid data values are = ',g14.7,', ',
     &f6.2,'% of grid')
      go to 300
110   print*,'No valid data points in grid'
      go to 300
c
c  ERROR message
c
299   print*, 'nz must =1, max ncol = 5000 '
      print*, 'something wrong...'
      print*
      close(10)
      go to 300
99    print*,'Input file not found '
      print*,'TRY AGAIN [y] '
      read(*,'(a)') answer
      if(answer.ne.'n'.and.answer.ne.'N') go to 6
c
c END OF PROGRAM
c
300   pause 'press ENTER to return to menu'
400   close(10)
      return
      end
c
      subroutine row1io(z,ncol,iunit,key)
c  for reading and writing rows when nz=1
      dimension z(ncol)
      go to (1,2),key
1     read(iunit) dum,z
      go to 90
2     write(iunit) dum,z
90    return
      end
c***************************************
      subroutine opngrd(igrid,ifile,ierr,*,*)
      character*(*) ifile
      character*64 title
      logical opn
        inquire(unit=igrid,opened=opn)
          if(opn) close(igrid)
        open(unit=igrid,file=ifile,
     &  status='old',form='unformatted',err=442,iostat=icheck)
        read(igrid,err=442,end=442,iostat=icheck) title,ncol
        rewind igrid
        return
442     call fcheck(igrid,ifile,icheck,ierr)
        if(ierr.eq.0) return 1
        if(ierr.eq.1) return 2
        if(ierr.eq.2) then
          ierr = 0
          return
        endif
c           if(ncol.gt.10000) then
c             open(unit=igrid,file=ifile,access='sequential',
c     &       status='old',form='unformatted')
c             access='sequen'
c           else
c             nbytes=(ncol+1)*4
c             open(unit=igrid,file=ifile,access='direct',
c     &       recl=nbytes,status='old',form='unformatted')
c             access='direct'
       return
       end
c***************************************
      subroutine fcheck(igrid,fname,ierror,ierr)
      character*(*) fname
      character*1 resp
      logical exists
      inquire(file=fname,exist=exists)
        if(.not.exists) then
          print *,'can''t find file ',fname
          fname=' '
          write(*,'(a\)') ' TRY AGAIN ? [y]: '
          read(*,'(a)') resp
          if(resp.eq.'n'.or.resp.eq.'N') then
            ierr=1
            return
          else
            ierr=0
            return
          endif
        endif
      call iof(igrid,fname,ierror,ierr)
      if(ierr.eq.1) go to 10
      return
10    print *,'error no. ',ierror,' with file: ',fname
      fname=' '
      return
      end
c***************************************************************
      subroutine iof(igrid,fname,ierror,ierr)
      parameter (numz=4500)
      character*(*) fname
      character*150 string
      dimension z(4500)
      close(igrid)
      jgrid=igrid+1
      kgrid=jgrid+1
      open(jgrid,file=fname,status='old',err=90,iostat=ierror,
     & form='formatted')
      read(jgrid,'(a)',err=90,iostat=ierror) string
c*****FORMATTED INPUT DATA FILE
c*****CREATE A SCRATCH UNFORMATTED FILE FOR INPUT TO PROGRAM
      open(igrid,status='scratch',form='unformatted')
      open(kgrid,status='scratch',form='formatted')
      lens = len_trim(string)
      write(kgrid,'(4a)') '''',string(1:64),'''',string(65:lens)
      rewind kgrid
      read(kgrid,*,err=90,iostat=ierror) string(1:64),nc,nr,nz,
     & xo,dx,yo,dy
      close(kgrid)
      if(nc.eq.0.or.nr.eq.0.or.nz.ne.1) then
        print *,'grid file has wrong format'
        go to 90
      endif
      if(nc.gt.numz) then
        print *,'number of input columns ',nc,' > ',numz
        go to 90
      endif
      write(igrid) string(1:64),nc,nr,nz,xo,dx,yo,dy
      iend=0
10    call iorw(jgrid,igrid,z,nc,iend,ierror,ierr)
      if(iend.ne.0) go to 20
      if(ierr.eq.1) go to 20
      go to 10
20    close(jgrid)
      rewind igrid
      return
90    ierr = 1
      close(jgrid)
      close(kgrid)
      return
      end
c***************************************
      subroutine iorw(jgrid,igrid,z,nc,iend,ierror,ierr)
      dimension z(nc)
      read(jgrid,*,end=10,err=20,iostat=ierror) y,z
      write(igrid) y,z
      return
10    iend=1
      ierr = 2
      return
20    ierr = 1
      return
      end
