c  To calculate scaling factor s and datum c such that
c  sum[(sF+c-G)**2]=min, where F and G are grids, s and c
c  are scalar constants.  Grid F can be rescaled, as an 
c  option
c
c   Grids same size and registered, dvals ok.
c
c  Written by Lin Cordell
c
c
      character*50 ifile1,ifile2
      character prompt*80,id*56
c      dimension f(2000),g(2000),id(14),pgm(2),ans(1)
      dimension f(2000),g(2000),pgm(2),ans(1)
      data dval/0.170412e+39/
      dval2=1.e20
c      write(6,1)
c1     format(1x,'To find s,c | sum[(sF+c-G)**2]=min.'/
c     1  1x,'Enter grid F name:'/' *'$)
c        read(5,100)ifile1
c100     format(a50)
        call pfinit('scale')
        call askin
        print*,'To find s,c | sum[(sF+c-G)**2]=min.'
        ifile1=' '
400     prompt='Enter grid F name'
        call askc(prompt,ifile1,ierr)
c        write(6,201)
        if(ierr.lt.0) stop
c201     format(1x,'Enter grid G name:'/' *'$)
c        read(5,100)ifile2
c        open(10,file=ifile1,form='unformatted',status='old')
c      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
        call gopen(10,ifile1,'old','read',ierr)
        if(ierr.ne.0) stop 'error opening file'
        call gheader('r',10,id,ncol,nrow,x0,dx,y0,dy,ierr)
        if(ierr.ne.0) stop 'error reading header'
        ifile2=' '
401     prompt='Enter grid G name'
        call askc(prompt,ifile2,ierr)
        if(ierr.lt.0) then
          close(10)
          go to 400
        endif
c      if(ncol.gt.2000.or.nz.ne.1.or.dx.eq.0.e0) go to 900
      if(ncol.gt.2000.or.dx.eq.0.e0) go to 900
c        open(11,file=ifile2,form='unformatted',status='old')
c        read(11)id,pgm,ncol2,nrow2,nz2,x02,dx2,y02,dy2
        call gopen(11,ifile2,'old','read',ierr)
        if(ierr.ne.0) stop 'error opening file'
        call gheader('r',11,id,ncol2,nrow2,x02,dx2,y02,dy2,ierr)
        if(ierr.ne.0) stop 'error reading header'
	if(ncol2.ne.ncol.or.nrow2.ne.nrow)then
	  print*,'Grids are not the same size.'
	  go to 900
      endif
c      if(nz.ne.1.or.nz2.ne.1)then
c        print*,'nz not equal to 1.'
c        go to 900
c      endif
      idif=0
c      if(nz2.ne.1.or.dx2.ne.dx.or.dy2.ne.dy)idif=1
      if(dx2.ne.dx.or.dy2.ne.dy)idif=1
      if(x02.ne.x0.or.y02.ne.y0)idif=1
      if(idif.eq.1)then
	  print*,'Grids are same size but origins and/or'
c          print*,'intervals differ.  Want to continue ?'
402          prompt='intervals differ.  Want to continue ?'
          call aski4l(prompt,ians,ierr)
          if(ierr.eq.-2) then
            close(10)
            close(11)
            go to 400
          endif
c          read(5,102)ans(1)
c          if(ans(1).eq.'y'.or.ans(1).eq.'Y')go to 5
          if(ians.eq.1) go to 5
	  go to 999
      endif
5     continue
	sf=0.0
	sg=0.0
	sgf=0.0
	sff=0.0
	do 20 j=1,nrow
c        call rowio(f,y,1,ncol,10,1)
c        call rowio(g,y,1,ncol,11,1)
        call grow('r',10,j,f,ncol,ierr)
        call grow('r',11,j,g,ncol,ierr)
	do 10 i=1,ncol
	if(f(i).gt.dval2)go to 10
	if(g(i).gt.dval2)go to 10
	sf=sf+f(i)
	sg=sg+g(i)
	sgf=sgf+g(i)*f(i)
	sff=sff+f(i)*f(i)
	n=n+1
10      continue
20      continue
	rn=n
	rmax=ncol*nrow
	pcent=(rn*100)/rmax
	if(n.eq.0)go to 900
	if(sff.eq.0.0)go to 900
	denom=rn*sff-(sf*sf)
	if(denom.eq.0.0)go to 900
	s=(rn*sgf-sf*sg)/denom
	c=(sg*sff-sf*sgf)/denom
	write(6,300) s,c,pcent
300     format(1x,'s = ',e12.5,3x,'c = ',e12.5,'  % =',f5.1)
      close(11)
c      print*,'Want to rescale the first grid ?'
c      print 101
c101   format(' *'$)
c      read(5,102)ans(1)
c102   format(a1)
c      if(ans(1).ne.'y'.and.ans(1).ne.'Y')go to 999
403    prompt='Want to rescale the first grid?'
       call aski4l(prompt,ians,ierr)
       if(ierr.eq.-2) then
         close(10)
         close(11)
         go to 400
       endif
       if(ians.ne.1) go to 999
c30    print*,'Enter output grid name:'
c      print 101
c      read(5,100)ifile1
c      open(12,file=ifile1,form='unformatted',status='new')
30     ifile1=' '
       prompt='Enter output grid name'
404    call askc(prompt,ifile1,ierr)
       if(ierr.eq.-2) go to 403
       call gopen(12,ifile1,'new','write',ierr)
       if(ierr.ne.0) stop 'error opening output file'
c      print*,'Enter id for new grid:'
c      print 101
c      read(5,103)(id(i),i=1,14)
c103   format(14a4)
c      write(12)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
405    prompt='Enter id for new grid'
       call askc(prompt,id,ierr)
       if(ierr.eq.-2) then
         close(12)
         go to 404
       endif
       call gheader('w',12,id,ncol,nrow,x0,dx,y0,dy,ierr)
      rewind (unit=10)
      read(10)id
      do 16 j=1,nrow
c      call rowio(f,y,nz,ncol,10,1)
      call grow('r',10,j,f,ncol,ierr)
      do 15 i=1,ncol
      if(f(i).gt.dval2)then
	g(i)=dval
	go to 15
      endif
      g(i)=f(i)*s+c
15    continue
c      call rowio(g,y,nz,ncol,12,2)
      call grow('w',12,j,g,ncol,ierr)
16    continue
c      print*,'Want to rescale another grid using these values?'
c      print 101
c      read(5,102)ans(1)
c      if(ans(1).ne.'y'.and.ans(1).ne.'Y')go to 999
      ians=0
406   prompt='Want to rescale another grid using this multiplier?'
      call aski4l(prompt,ians,ierr)
      if(ierr.eq.-2) then
        rewind(12)
        go to 405
      endif
      if(ians.ne.1) go to 999
      close(10)
      close(12)
c      print*,'Enter grid F name:'
c      print 101
c      read(5,100)ifile1
c      open(10,file=ifile1,form='unformatted',status='old')
      ifile1=' '
407   prompt='Enter grid F name'
      call askc(prompt,ifile1,ierr)
      if(ierr.eq.-2) go to 406
       call gopen(10,ifile1,'old','read',ierr)
      c=0.0
      go to 30
900   write(6,901)
901   format(' Can''t handle it.')
999     close(10)
      stop
      end
c      subroutine rowio(a,y,nz,ncol,ld,key)
c      dimension a(nz,ncol)
c      go to (1,2),key
c1     read(ld) y,a
c      go to 90
c2     write(ld)y,a
c90    return
c      end
