c  grid_stat.for
c
c  To calculate mean, RMS, and standard deviation on input
c  grid f, or on difference grid d=f-g of two registered
c  grids, if the second grid g is given. Dvals are skipped.
c
c  Lin Cordell 18 June 88.
      character*50 ifile1,ifile2
	character*1 ans
      dimension f(2000),g(2000),id(14),pgm(2)	
c      data dval/0.170412e+39/
	dval2=1.e20 
      write(6,1)
1     format(1x,'Enter input grid f name:'/' *'$)
	read(5,100)ifile1
	open(10,file=ifile1,form='unformatted',status='old')
	read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
	if(ncol.gt.2000.or.nz.ne.1.or.dx.eq.0.e0) go to 900
	do 7 i=1,ncol
7	g(i)=0.0
2	write(6,3)
3	format(' Want to enter 2nd grid g, for f-g ?:'/' *'$)
	read(5,4)ans
4	format(a1)
	if(ans.eq.'n')go to 6
	if(ans.ne.'y')go to 2
      write(6,5)
5     format(' Enter input grid g name:'/' *'$)
	read(5,100)ifile2
100	format(a50)
	open(11,file=ifile2,form='unformatted',status='old')
	read(11)id,pgm,ncol2,nrow2,nz,x02,dx2,y02,dy2
	if(ncol.gt.2000.or.nz.ne.1.or.dx.eq.0.e0) go to 900
	if(ncol2.ne.ncol.or.nrow2.ne.nrow)go to 900
	if(dx2.ne.dx.or.dy2.ne.dy)go to 900
	if(x02.ne.x0.or.y02.ne.y0)go to 900
6	continue
	n=0
	sum=0.0
	sumsq=0.0
	do 20 j=1,nrow
	call rowio(f,y,1,ncol,10,1)
	if(ans.eq.'y')call rowio(g,y,1,ncol,11,1)
	do 10 i=1,ncol
	if(f(i).gt.dval2)go to 10
	if(g(i).gt.dval2)go to 10
	t=f(i)-g(i)
	tsq=t*t
	sum=sum+t
	sumsq=sumsq+tsq
	n=n+1
10	continue
20	continue
	rn=n
	if(n.gt.1)go to 30
	write(6,21)
21	format(' Empty grid(s).')
	go to 999
30	continue
	rmean=sum/rn
	rms=sqrt(sumsq/rn)
	sigma=(sumsq-(sum*sum)/rn)/(rn-1)
	sigma=sqrt(abs(sigma))
	rmax=ncol*nrow
	pcent=(rn*100)/rmax
	write(6,300)rmean,rms,sigma,pcent
300     format(1x,' Mean value =',e12.5,'; RMS value =',e12.5,/,
     1 '  Standard deviation = ',e12.5,'; Percent data =',f7.2)
      go to 999
900   write(6,901)
901   format(' Can''t handle it.')
999	close(10)
	if(ans.eq.'y')close(11)
      stop
      end
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
