      dimension f1(1000),f2(1000),dif(1000)
      character fil1*50,fil2*50,title*56,pgm*8,exit*2,ofil*50
c      common /dummy/a(28500,3)
c      common /dummy2/u(28500),y(28500)
      dimension a(10000,3),u(10000),y(10000)
      double precision a,u,y,out(3)
      data exit/'ex'/,ir/5/
      bval=99999999.
  1   write(6,*) ' Enter name of 1st input file (one to be modified) or 
     1"ex" '
      read(ir,10)fil1
10    format(a50)
      if(fil1.eq.exit) go to 99
      write(6,*) ' Enter name of 2nd input file '
      read(ir,10)fil2
      open(10,file=fil1,status='old',form='unformatted')
c     call fatt('file10',fil1,'-append')
      read(10)title,pgm,nc1,nr1,nz,x1,dx1,y1,dy1
      open(12,file=fil2,status='old',form='unformatted')
c     call fatt('file12',fil2,'-append')
      read(12)title,pgm,nc2,nr2,nz,x2,dx2,y2,dy2
      write(6,*) 'Enter name of output file '
      read(ir,10)ofil
      open(11,file=ofil,status='new',form='unformatted')
      if(dx1.ne.dx2.or.dy1.ne.dy2) go to 98
c     comout='compgrd.out'
c     call fatt('file20',comout,'-append')
c     open(20,file=comout,status='new')
c     write(20,20)fil1,nr1,nc1,y1,x1,dy1,dx1,fil2,nr2,nc2,y2,x2,dy2,dx2
20    format(///,' comparison of two grids'/' file 1= ',a50
     &/' no. of rows & cols=',
     &2i4,/,' first row & col coord=',2e12.5/' spacing along rows & 
     &cols=',2e12.5//' file 2= ',a50/' no. of rows & cols=',2i4/
     &' first row & col coord.=',2e12.5/' spacing along rows & cols=',
     &2e12.5,///)
      ye1=y1+(nr1-1)*dy1
      ye2=y2+(nr2-1)*dy1
      ye=ye1
      if(ye1.gt.ye2)ye=ye2
      if(y1.lt.y2) go to 30
      mf=12
      nrb=(y1-y2)/dy1+1.000001
      nre=(ye-y1)/dy1+1.000001
      go to 40
30      mf=10
      nrb=(y2-y1)/dy1+1.000001
      nre=(ye-y2)/dy2+1.000001
40    xe1=x1+(nc1-1)*dx1
      xe2=x2+(nc2-1)*dx1
      xe=xe1
      if(xe1.gt.xe2)xe=xe2
      if(x1.lt.x2) go to 50
      nce=(xe-x1)/dx1+1.000001
      ncb1=1
      ncb2=(x1-x2)/dx1+1.000001
      nce1=(xe-x1)/dx1+1.000001
      nce2=(xe-x2)/dx1+1.000001
      go to 60
50      ncb1=(x2-x1)/dx1+1.000001
      nce=(xe-x2)/dx1+1.000001
      ncb2=1
      nce1=(xe-x1)/dx1+1.000001
      nce2=(xe-x2)/dx1+1.000001
60    if(nrb.eq.1) go to 72
      do 70 i=1,nrb-1
70    read(mf)dum
72    sumt=0.
      nt=0
      nsum=0
      sum=0.
      kk=0
      do 80 i=1,nre
      read(10)dum,(f1(j),j=1,nce1)
      read(12)dum,(f2(j),j=1,nce2)
      j1=ncb1-1
      j2=ncb2-1
      sum=0.
      nsum=0
      do 90 j=1,nce
      dif(j)=0.1701411e+39
      j1=j1+1
      j2=j2+1
      if(f1(j1).gt.bval.or.f2(j2).gt.bval) go to 90
      dif(j)=f1(j1)-f2(j2)
      sum=sum+dif(j)
      nsum=nsum+1
      kk=kk+1
      a(kk,1)=j1-1
      a(kk,2)=i+nrb-2
      a(kk,3)=1.
      y(kk)=dif(j)
90    continue
      sumt=sumt+sum
      nt=nt+nsum
      if(nsum.eq.0) go to 80
      asum=sum/nsum
      write(6,100)i,asum
c     write(20,100)i,asum,(dif(j),j=1,nce)
80    continue
100   format(' row =',i4,'    average difference=',f12.2,5x,1000f10.2)
      if (nt.eq.0) go to 95
      asumt=sumt/nt
      write(6,110)asumt,nt
c     write(20,110)asumt,nt
110   format(//,' total average difference=',f10.2,4x,i5,'  values 
     & used')
      call hdecom(10000,kk,3,a,u)
      call hsolve(10000,kk,3,a,u,y)
      do 83 i=1,3
83    out(i)=y(i)
      rewind(10)
      write(6,103) out
103   format(3d24.9)
      read(10) title
      write(11) title,pgm,nc1,nr1,nz,x1,dx1,y1,dy1
      do 97 j=1,nr1
      read(10) dum,(f1(i),i=1,nc1)
      do 96 i=1,nc1
      if(f1(i).ge.bval) go to 96
      f1(i)=f1(i)-out(1)*(i-1)-out(2)*(j-1)-out(3)
96    continue
      write(11) dum,(f1(i),i=1,nc1)
97    continue
c     call fdet('file20')
c     close(20)
c98   call fdet('file10')
  98  close(10)
c     call fdet('file12')
      close(12)
      close(11)
      go to 1
95    print *,'no overlapping points'
c99   call closet
  99  stop
      end
      subroutine hdecom(mdim,m,n,a[HUGE],u[HUGE])
      integer mdim,m,n
      real*8 a(mdim,n),u(m)
c
c  householder reduction of rectangular matrix to upper triangular form
c  use with hsolve for least squares solution of overdetermined systems
c  parameters..
c    mdim = declared row dimension of a
c    m = number of rows of a
c    n = number of columns of a
c    a = m-by-n matrix, m.ge.n
c        input, matrix to be reduced
c        output, reduced matrix and information about reduction
c    u = m-vector
c        input, ignored and destroyed
c        output, 'diagonal' of reduced matrix
c
      real*8 alpha,beta,gamma,dsqrt
c
      do 6 k=1,n
c
c        find reflection which zeroes a(i,k),i=k+1,...,m
c
      alpha=0.
      do 1 i=k,m
      u(i)=a(i,k)
      alpha=alpha+u(i)**2
c      print*,'alpha ',alpha
    1 continue
      alpha=dsqrt(alpha)
      if(u(k).lt.0.)alpha=-alpha
      u(k)=u(k)+alpha
      beta=alpha*u(k)
      if(k.eq.n)go to 5
c
c        apply reflection to remaining columns of a
c
      kp1=k+1
      do 4 j=kp1,n
      gamma=0.
      do 2 i=k,m
c      print*,'gamma ',gamma
      gamma=gamma+a(i,j)*u(i)
    2 continue
      gamma=gamma/beta
      do 3 i=k,m
      a(i,j)=a(i,j)-gamma*u(i)
    3 continue
    4 continue
    5 a(k,k)=u(k)
      u(k)=-alpha
    6 continue
      return
c
c  triangular result stored in u(i) and a(i,j), i.lt.j
c  vectors defining reflections stored in a(i,k), i.ge.k
c
      end
      subroutine hsolve(mdim,m,n,a[HUGE],u[HUGE],y[HUGE])
      integer mdim,m,n
      real*8 a(mdim,n),u(m),y(m)
c
c  least squares solution of overdetermined system
c  find x that minimizes norm(a*x-y)
c  parameters..
c    mdim,m,n,a,u..results from hdecom
c    y = m-vector
c        input, right hand side
c        output, forst n components = the solution, x
c            last m-n components = transformed residual
c
      real*8 beta,gamma
c
c        apply reflections to y
c
      do 9 k=1,n
      beta=-u(k)*a(k,k)
      gamma=0.
      do 7 i=k,m
      gamma=gamma+y(i)*a(i,k)
    7 continue
      gamma=gamma/beta
      do 8 i=k,m
      y(i)=y(i)-gamma*a(i,k)
    8 continue
    9 continue
c
c        back substitution
c
      y(n)=y(n)/u(n)
      nm1=n-1
      do 11 ib=1,nm1
      i=n-ib
      ip1=i+1
      do 10 j=ip1,n
      y(i)=y(i)-a(i,j)*y(j)
   10 continue
      y(i)=y(i)/u(i)
   11 continue
      return
      end
