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  Tien Grauch
c
      character infile*50,id*56,answer*1,iexs*50,fmt*40,prompt*80
      character*1 iex(50)
      dimension dell(2),range(2)
      common/block/a(5000),iclass(200),iclasf(200)
      equivalence (iexs,iex(1))
      data amax/-1.0e37/,amin/1.0e37/,dcval/1.0e38/
      kmax=0
      ihistmx=0
      call askin
      do 5 i=1,50
      iclass(i)=0
      iex(i)='X'
5     continue
      infile=' '
6     call askc('Enter input file',infile,ierr)
      if(ierr.eq.-2) stop
      call gopen(10,infile,'old','read',ierr)
      if(ierr.ne.0) go to 99
900   format(a50)
      call gheader('r',10,id,ncol,nrow,x0,dx,y0,dy,ierr)
      if(ierr.ne.0) stop 'error reading grid header'
      if(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 grow('r',10,idum,a,ncol,ierr)
      if(ierr.ne.0) stop 'error reading row'
      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
      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')
806   answer='n'
      prompt=
     1'Want to see histogram, standard dev., both, or none (h,s,b,n)?'
      call askc(prompt,answer,ierr)
      if(ierr.eq.-2) then
        close(10)
        go to 6
      endif
      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
550   lcount=lcount+1
      dell(1)=amax-amean
      dell(2)=amean-amin
      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=amax-amin
      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=(amax-amin)/del
	print 999,del,icheck
999	format(/' I found histogram interval of ',g14.7,', giving '
     1 ,i3,' histogram classes')
700     answer='n'
        call askc('Want to redefine interval?',answer,ierr)
        if(ierr.eq.-2) go to 806
	if(answer.ne.'y'.and.answer.ne.'Y')go to 705
701   call askf4('Enter new interval',del,ierr)
      if(ierr.eq.-2) go to 700
      icheck=(amax-amin)/del
      if(icheck.gt.200.) then
         del=(amax-amin)/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     answer='p'
        prompt='Want to see percent or count of values (p or c)?'
        call askc(prompt,answer,ierr)
        if(ierr.eq.-2) go to 806
	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
      sum2=0.
c round off amin & amax & rename as range(1) & range(2)
      range(1)=amin
      range(2)=amax
      do 11 i=1,200
11    iclass(i)=0
      kmax=0
      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
      call gheader('r',10,id,ncol,nrow,x0,dx,y0,dy,ierr)
      do 31 j=1,nrow
      call grow('r',10,idum,a,ncol,ierr)
      if(ierr.ne.0) stop 'error reading row'
      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   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      print 803,del
803   format(/,' midpoint',1x,'(interval=',g14.4,')',3x,'HISTOGRAM',/)
      pt=range(1)-del
      do 60 i=1,kmax
	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(*,fmt)pt,iexs,label
        endif
	if(newopt.eq.1)then
        write(fmt,714)iclasf(i)
        write(*,fmt)pt,iexs,rlabel
        endif
714     format('(1x,g10.4,1x,''|'',a',i6,',1x,f6.1,''%'')')
804     format('(1x,g10.4,1x,''|'',a',i6,',1x,i6)')
      go to 60
59    if(iclass(i).eq.0) go to 599
	if(newopt.eq.2)print 805,pt,label
	if(newopt.eq.1)print 715,pt,rlabel
715	format(1x,g10.4,1x,'|',1x,f6.1,'%')
805   format(1x,g10.4,1x,'|',1x,i6)
      go to 60
599   print 805,pt
60    continue
        answer='s'
        prompt='REPRINT info, CHANGE histogram, or STOP (r, c, or s)?'
        call askc(prompt,answer,ierr)
        if(ierr.eq.-2) go to 699
	if(answer.eq.'r'.or.answer.eq.'R') go to 441
	if(answer.eq.'c'.or.answer.eq.'C') go to 699
      go to 300
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...'
      call gclose(10,'keep')
      go to 300
99    print*,'Input file not found or busy.  Try again'
      go to 6
c
c END OF PROGRAM
c
300   call gclose(10,'keep')
c      stop
      end
