c  2-d surface fitting by orthogonal polynomials
c  coded by ron wahl,us geological survey
c  modified for incomplete grids by mike webring, usgs
c
	common /surfa/ fx(20,20),fy(20,20),zm(300),vx(300),vy(1000),
     1 alk(20,20),zsq(20,20),resd(40),zc(300),zres(500)
	common /surfb/ ave,ihdr(23),if1(14)
	common /iflag/ iflag,xi(300),yi(1000)
	dimension iopt(10),ys(301)
        double precision fx,fy,zm,vx,vy,alk,sumz,sumzsq,zsq,zc,zres,
     1  resd  
        character pgm*8,tmp*10
	character*56 ifile
	integer dimt,cm,cn,cnm
	equivalence (if1(1),ifile)
	equivalence (iopt(10),ires),(ihdr(15),pgm),(ihdr(17),cm),
     1 (ihdr(18),cn),(ihdr(21),delx),(ihdr(23),dely),(zc(1),ys(1))
        data tmp/'surfit.res'/,
     1 is/1/,il/1/,inp/5/,isw/12/,jsw/13/,ksw/14/,lsw/15/
c
c	call user(pgm1)
	nwork=8000
	iflag=1
	do 1 i=1,10
1	iopt(i)=0
	iopt(7)=1
	write (*,2)
2	format(' enter grid filename :'$)
	read(inp,37) ifile
37	format(a56)
	open(unit=isw,file=ifile,status='old',form='unformatted')
	write(*,38)
38    format(' start, end, increment of orders, [e.g., 2,5,1]')
      read(inp,*) il,im,is
      print*,'Want grids of residuals? [0=no, 1=yes]'
      read*,ires
      if(ires.ne.1)ires=0
	if(il.le.0) il=1
	if(im.gt.19) im=19
	if(is.le.0) is=1
	im=((im-il)/is)*is+il
	write(*,103)il,im,is
103	format(' start, end, increment :',3i5)
	open(unit=8,file='surfit.lis',status='unknown',form='formatted')
	open(unit=jsw,file=tmp,status='unknown',form='unformatted')
        call slevel(nwork,isw,jsw,1,ave)
	rewind jsw
	read (jsw) ihdr
	if (cn.le.1000.and.cm.le.300)  go to 5
	stop ' nc>300 or nr>1000'
c5	pgm=pgm1
5       continue
c
c     set up x and y coordinates and compute the orthogonal
c     polynomial coefficients
	dimt=im+1
	if (delx.ne.0)  go to 10
	call read(xi,cm,jsw)
	iflag=0
	go to 30
10	xi(1)=-1.
	step=2.0/float(cm-1)
	do 20 i=2,cm
	xi(i)=xi(i-1)+step
20	continue
30	if (dely.ne.0)  go to 50
	ncol1=cm+1
	do 40 j=1,cn
	call read(ys,ncol1,jsw)
	yi(j)=ys(1)
40	continue
	go to 70
50	yi(1)=-1.
	step=2./float(cn-1)
	do 60 j=2,cn
	yi(j)=yi(j-1)+step
60	continue
70	do 80 i=1,cm
80	vx(i)=dble(xi(i))
	do 90 j=1,cn
	vy(j)=dble(yi(j))
90	continue
	call opoly(vx,cm,dimt,fx)
	call opoly(vy,cn,dimt,fy)
	call coef(jsw,fx,fy,vx,vy,zm,alk,zsq,dimt,cn,cm,iopt,sumz,sumzsq)
	call rsdl(sumzsq,zsq,resd,il,is,im,iout,cn,iopt,alk)
	cnm=cn*cm
	write (*,102)
102	format(' least_square coefficients completed')
	do 100 i=il,iout,is
	idim=i+1
	call scoef(alk,fy,fx,zsq,idim,iopt,cn,cnm)
	call zval(isw,sumz,vx,vy,zsq,cn,cm,cnm,zm,zc,zres,resd,i,iopt)
	write (*,101)i
101	format(' order',i3,' finished')
	close(ksw)
	if(ires.ne.0) close(lsw)
100	continue
	close(8)
	close(isw)
	close(jsw,status='delete')
110	stop
	end
	subroutine slevel(nw,is,js,iave,ave)
        common /surfa/w(8080)
	dimension ihdr(23)
	equivalence (ihdr(17),nc),(ihdr(18),nr)
	rewind is
	rewind js
	iout=js
	if(iave.eq.0) go to 11
	iout=22
	open(unit=iout,file='grid1.tmp',status='unknown',
     1 form='unformatted')
11	read(is) ihdr
	write(iout) ihdr
	nt=1
	nrb=nr
1	nn=max(nc,nrb)+nrb+nc*nrb
	if(nn.lt.nw) go to 2
	nrb=nrb-1
	if(nrb.lt.1) stop 11
	go to 1
2	nn=(nt-1)*(nrb-1)+nrb
	if(nn.ge.nr) go to 3
	nt=nt+1
	if(nt.gt.20) stop 12
	go to 2
3	nrb=nrb-(nn-nr)/nt
	nout=(nt-1)*(nrb-1)+nrb
	if(nout.ge.nr) go to 4
	write (*,5)
5	format(' slevel parm warning')
	nrb=nrb+1
4	iy=nc*nrb+1
	iaux=iy+nrb
	call slevl2(w,w(iaux),w(iy),nc,nr,nrb,is,iout,iave,ave)
	rewind is
	rewind iout
	if(iave.eq.0) return
	read(iout) ihdr
	write(js) ihdr
	do 7 j=1,nr 
	call rw(nc,w,y,-1,iout,js,ie)
	do 6 i=1,nc
6	w(i)=w(i)-ave
	call rw(nc,w,y,0,iout,js,ie)
7	continue
	close(iout,status='delete')
	rewind js
	return
	end
	subroutine slevl2(z,w,y,nc,nr,nrb,is,js,iave,ave)
	dimension z(1),w(1),y(1)
	nn=nc*nr
	rwt=1.0/float(nn)
	irt=0
	ir=0
	ave=0.0
1	if(ir.eq.nrb .or. irt.eq.nr) go to 2
	m=ir*nc+1                    
	call rw(nc,z(m),y(ir+1),-1,is,js,ier)
	ir=ir+1
	irt=irt+1
	go to 1
2	call gridr(nc,nrb,z,w,ier)
	if(iave.eq.0) go to 4
	do 3 i=1,ir*nc
3	ave=ave+z(i)*rwt
4	nro=ir-1
	if(irt.eq.nr) nro=ir
	m=1
	do 5 i=1,nro     
	call rw(nc,z(m),y(i),0,is,js,ier)
5	m=m+nc
	if(irt.eq.nr) return
	m=(nrb-1)*nc+1
	do 6 i=1,nc
	z(i)=z(m)
6	m=m+1
	ir=1
	go to 1
	end
	subroutine rw(n,z,y,iop,idev,jdev,iend)
	dimension z(n)
	iend=0      
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
	subroutine write(a,mx,lu)
	dimension a(mx)
	write (lu)  a
	return
	end
      subroutine read(a,mx,insw)
      dimension a(mx)
      read (insw) a
      return
      end
	subroutine inputz(zm,ny,mx,insw)
	double precision zm
	common /iflag/ iflag,xi(300),yi(1000)
	dimension zm(1),a(301)
	mxx=mx+1
	if (ny.gt.1)  go to 10
	if (ny.lt.1) stop 1111
	rewind insw
	read (insw)
	if (iflag.eq.0)  read (insw)	
10	call read(a,mxx,insw)
	do 20 l=1,mx	
	zm(l)=dble(a(l+1))
20	continue
	return
	end
	subroutine zcrout(zc,zres,ns,mx,ideg,ir)
	common /surfb/ ave,ihdr(23),if1(14)
	common /iflag/ iflag,xi(300),yi(1000)
	dimension a(301)
	double precision zc(1),zres(1)
	character*56 ifile,outfil
	character*4 extent(19),rxtent(19)
	equivalence (if1,ifile)
	data extent/'.s01','.s02','.s03','.s04','.s05','.s06',
     1  '.s07','.s08','.s09','.s10','.s11','.s12','.s13','.s14',
     1  '.s15','.s16','.s17','.s18','.s19'/
	data rxtent/'.r01','.r02','.r03','.r04','.r05','.r06',
     1  '.r07','.r08','.r09','.r10','.r11','.r12','.r13','.r14',
     1  '.r15','.r16','.r17','.r18','.r19'/
	ksw=14
	lsw=15
	mxx=mx+1
	if (ns.gt.1)  go to 10
	if (ns.lt.1)  stop 2222
c  setup output files
c  concatenate ifile and extent(ideg)
	n=index(ifile,'.')-1
	if(n.le.0) stop ' short input filename'
	outfil=ifile(1:n)//extent(ideg)
	open(unit=ksw,file=outfil,status='unknown',
     1 form='unformatted')
	write (ksw) ihdr
	if (iflag.eq.0)  call write(xi,mx,ksw)
	if (ir.eq.0)  go to 10
	outfil=ifile(1:n)//rxtent(ideg)
	open(unit=lsw,file=outfil,status='unknown',
     1 form='unformatted')
	write (lsw) ihdr
	if (iflag.eq.0)  call write(xi,mx,lsw)
c  output data row
10	a(1)=yi(ns)
	do 20 i=1,mx
	a(i+1)=sngl(zc(i))
20	continue
	call write(a,mxx,ksw)
	if (ir.eq.0)  go to 900
	do 30 i=1,mx
	a(i+1)=sngl(zres(i))
30	continue  
	call write(a,mxx,lsw)
900	return
	end
	subroutine zval(insw,sumz,vx,vy,b,ny,mx,nymx,zm,zc,
     1 zres,resd,ideg,iopt)
c
c     this routine using the 'b' coefficients produced by scoef,
c     computes the value of the zc(x,y)=f(x,y) surface at each data
c     point and the residual values zres(x,y)=zm(x,y)-zc(x,y). the
c     computed surface values and residual values are output by a user
c     supplied routine. statistics produced concerning the measured and
c     computed 'z' values are output in this routine.
c
c     declarations         
	common /surfb/ ave,ihdr(23),if1(14)
	character*56 ifile
	double precision sumz,vx,vy,b,zm,zc,zres,resd,xy,poly,
     1 sigsqr,sumsqr,zrmax,zrmin,zmax,zmin,zresl,vyns,
     1 vxm,ypower,sumzc,sigsrt,stdr,sig
	dimension b(20,20),vx(1),vy(1),zc(1),zm(1),zres(1),ypower(20),
     1   iopt(10),resd(1)
	equivalence (if1,ifile)
c     ******************************************************************
c     *  sumz--sum of measured 'z' (zm) values.                        *
c     *  vx--an array containing the 'x' coordinates.                  *
c     *  vy--an array containing the 'y' coordinates.                  *
c     *  b--an array containing the coefficients, b, of the f(x,y)     *
c     *      surface.                                                  *
c     *  ny--the number of data rows.                                  *
c     *  mx--the number of data columns.                               *
c     *  nymx--the number of data points.                              *
c     *  zm--an array containing a row of measured 'z' values.         *
c     *  zc--an array containing a row of computed 'z' values.         *
c     *  zres--an array containing the actual residual values.         *
c     *  resd--an array containing the sum sq res(est).                *
c     *  ideg--the degree of the surface to be generated.              *
c     *  iopt--an array containing the program options.                 *
c     *  ------------------------------------------------------------  *
c     *  called suprograms:                                            *
c     *      inputz(zm,ns,mx)--a user supplied routine to input        *
c     *      measured 'z' by rows.                                     *
c     *      zcrout(zc,zres,ns,mx)--a user supplied routine to output  *
c            computed 'z' and residuals by rows.                       *
c     *      poly*8(a,w,k,l)--computes the value of a polynomial.      *
c     ******************************************************************
c
c     initialize necessary vaiables.
c     data dv/'ffff7fff'x/ 
        data dv /1.e+38/
	lout=8
	rewind insw
	xy=dble(nymx)
	sigsqr=0.0d0
	sumsqr=0.0d0
	idim=ideg+1
	zmean=sumz/xy
	zrmax=1.0d-35
	zrmin=1.0d35
	assign 10 to nysw
	if (ny.eq.1)  assign 30 to nysw
	assign 60 to ksw
	if (iopt(5).ne.0)  assign 50 to ksw
c
c     evaluate surface a row at a time and compute the residuals
c     zres(m)=zm(m)-zc(m)
	do 110 ns=1,ny
	vyns=vy(ns)
	call inputz(zm,ns,mx,insw)
	ypower(1)=1.0d0
	go to nysw,(10,30)
10	do 20 n=2,idim
20	ypower(n)=ypower(n-1)*vyns
	idimt=idim
	go to 40
30	idimt=1
40	do 100 m=1,mx
	vxm=vx(m)
	sumzc=0.0d0
	do 80 l=1,idimt
	go to ksw,(50,60)
50	k=idim
	go to 70
60	k=idim-l+1
70	sumzc=poly(b,vxm,k,l)*ypower(l)+sumzc
80	continue
	zc(m)=sumzc+ave
	if(zm(m).lt.1.e35) go to 81
	zres(m)=dv
	go to 100
81	zresl=zm(m)-zc(m)
	zres(m)=zresl
	sigsrt=zm(m)-zmean
	sigsqr=sigsrt*sigsrt+sigsqr
	sumsqr=zresl*zresl+sumsqr
	if (zrmax.ge.zresl)  go to 90
	zrmax=zresl
	zmax=zm(m)
	go to 100
90	if (zrmin.le.zresl)  go to 100
	zrmin=zresl
	zmin=zm(m)
100	continue
105	call zcrout(zc,zres,ns,mx,ideg,iopt(10))
110	continue
	sig=dsqrt(sigsqr/(xy-1.0d0))
	stdr=dsqrt(dabs(sumsqr/(xy-0.5d0*dfloat(ideg*idim))))
	write (lout,120) resd(ideg),zmean,stdr,sumsqr,sigsqr,
     1 sig,zmax,zrmax,zmin,zrmin
120	format(////' estimated sum sqrs res  ',e12.5,9x,
     1 'mean of observed values  ',e12.5,8x,
     1 'std err of est  ',e12.5/'    actual sum sqrs res  ',
     1 e12.5,11x,'sum sqr dev from mean  ',e12.5,8x,
     1 'std dev of obs  ',e12.5//1x,10x,'ranges'/1x,5x,
     1 'max',2(3x,e12.5)/1x,5x,'min',2(3x,e12.5))
	return
	end
	subroutine scoef(alk,fy,fx,b,mdeg,iopt,ny,nymx)
c
c     this routine computes the coefficients, b(l,k) in the equation
c     f(x,y)=sum(sum(b(l,k)*(x**k)*(y**l))) preparatory to evaluation
c     of the surface, f(x,y) at each point xi,y4.
c     ******************************************************************
c     *  alk--an array containing the least-square coefficients.       *
c     *  fy--an array containing the 'y' orthogonal polynomial         *
c     *      coefficients.                                             *
c     *  fx--an array containing the 'x' orthogonal polynomial         *
c     *      coefficients.                                             *
c     *  b--(formerly, if desired zsq) an array containing the 'b' s.  *
c     *  mdeg--the degree of the surface to be generated.              *
c     *  iopt--an array containing the program options.              *
c     *  ny--the number of data rows.                                  *
c     *  nymx--the number of data points.                              *
c     *  ------------------------------------------------------------  *
c     *  called subprograms :  none by name.                           *
c     ******************************************************************
c
c     declarations
	double precision alk,fy,fx,b,fyjl,sumb
	dimension alk(20,20),fy(20,20),fx(20,20),b(20,20),iopt(10)
c
c     initialize necessary variables
	lout=8
	ldeg=mdeg
	if (ny.eq.1)  ldeg=1
	assign 20 to lsw
	if (iopt(5).ne.0)  assign 10 to lsw
c
c     gen the 'b's' directly from sum(sum(alk(l,k)*psil(y)*phik(x)))
c     start the four loops to gen the 'b's'.
c
	do 50 l=1,ldeg
        goto (10,20),lsw
c	go to lsw,(10,20)
10	kdeg=mdeg
	go to 30
20	kdeg=mdeg-l+1
30	do 50 k=1,kdeg
	sumb=0.0d0
	jdeg=mdeg-k+1
	do 40 j=l,jdeg
	fyjl=fy(j,l)
	ideg=mdeg-j+1
	do 40 i=k,ideg
40	sumb=alk(j,i)*fx(i,k)*fyjl+sumb
50	b(l,k)=sumb
c
c     output the 'b's'.
	ldegm=mdeg-1
	write (lout,60)  ldegm,nymx
60	format ('1',10x,'degree =',i3,2x,'no. of obs =',
     1 i5/' least squares solution...coefficients of ascending',
     1 ' powers of x across and y down'//)
	i1=mdeg
	do 80 j=1,ldeg
	write (lout,70)  (b(j,i),i=1,i1)
70	format (10(1x,e12.5))
80	i1=mdeg-j
	return
	end
	subroutine rsdl(sumzsq,zsq,resd,dmin,dint,dmax,
     1 dout,ny,iopt,alk)
c     this subroutine computes the sum of the squares of the residuals,
c     (estimates) for the desired polynomial surfaces, or if iopt(6)
c     and (7).eq.0 the computation is carried out until the above sum <
c     7 orders of magnitude of the sum of the zsq matrix or the above
c     sum of the current surface > the above sum for the prev surface.
c     ******************************************************************
c     *  sumzsq--sum of the squares of the measured 'z' values.        *
c     *  zsq--an array containing the 'z-square' matrix.               *
c     *  resd--an array containing the sum sq resdl for each surface   *
c     *      to be generated                                           *
c     *  dmin-- the minimum degree surface to be generated             *
c     *  dint-- the change in degree between each surface generated.   *
c     *  dmax-- the maximum degree surface to be generated.            *
c     *  dout-- if dmax not input, ie if iopt(6) and (7)=0, this is *
c     *      the maximum degree surface to be generated.               *
c     *  ny--the number of rows of 'z'.                                *
c     *  iopt-- an array containing the program options.             *
c     *  alk--an array containing the least-square coefficients.       *
c     *  ------------------------------------------------------------  *
c     *  subprograms called:--dlog10                                   *
c     ******************************************************************
c
c     declarations
	double precision sumzsq,zsq,resd,alk,test,sumd,rlastd
	integer dmin,dint,dmax,dmint,dminp,dminpt,dimp,
     1 dmi,dout,dmaxt,cr
	logical temp
	dimension alk(20,20),zsq(20,20),temp(20,20),resd(40),iopt(10)
	cr=5
	lout=8
	write(lout,1)
1	format(////)
c	10**40 by indirection
	sumd=1.0d38
	dmaxt=dmax+1
	dmint=dmin-dint
	assign 150 to msw
	if (iopt(5).eq.0)  go to 10
	assign 100 to ksw
	itest=2
	go to 20
10	assign 110 to ksw
	itest=1
20	if (iopt(6).eq.0.and.iopt(7).eq.0)  go to 30
	assign 90 to msw
	dout=dmax
30	if (iopt(8).eq.0)  go to 90
	do 40 l=1,dmax
	do 40 k=1,dmax
40	temp(l,k)=.false.
50	write (*,51)
51	format(' enter irow,icol :'$)
	read (cr,60)  irow,icol
60	format (2i3)
	if (irow.eq.999)  go to 70
	temp(irow,icol)=.true.
	go to 50
70	do 80 l=1,dmax
	do 80 k=1,dmax
	if (temp(l,k))  go to 80
	alk(l,k)=0.0d0
	zsq(l,k)=0.0d0
80	continue
90	rlastd=sumd
	sumd=sumzsq
	dmint=dmint+dint
	if (dmint.gt.dmax)  return
	dmi=dmint*itest
	dminp=dmint+1
	dminpt=dminp
	if (ny.eq.1)  dminpt=1
	do 130 l=1,dminpt
	go to ksw,(100,110)
100	dimp=dminp
	go to 120
110	dimp=dminp-l+1
120	do 130 k=1,dimp
130	sumd=sumd-zsq(l,k)
	resd(dmi)=sumd
	write (lout,140)  dmi,resd(dmi)
140	format (' for degree= ',i2,' the sum sq res (est) = ',e13.5,/)
	go to msw,(90,150)
150	if(sumd.gt.rlastd)  go to 180
	if (sumd.gt.1.0d0)  go to 160
	sumd=1.0d0
160	test=dlog10(sumzsq)-dlog10(sumd)
	if (test.gt.7.0d0)  go to 90
	write (lout,170)
170	format (' higher order surface not computed since sum',
     1 ' sq res(est) < 7 orders of mag than sum of z-sq matrix.')
	go to 200
180	write (lout,190)
190	format (' higher order surface not computed since sum',
     1 ' sq res(est) (current) > last value.')
200	dout=dmint
	return
	end
	subroutine coef(insw,fx,fy,vx,vy,zm,alk,zsq,dimt,ny,mx,
     1 iopt,sumz,sumzsq)
c     this subroutine generates the least-square coefficients and the
c     z-square matrix for fitting a polynomial surface to the 'z' data
c     distributed in 'x' and 'y' by orthogonal polynomials.
c     ******************************************************************
c     *  fx--an array containing, in double precicion, the orthogonal  *
c     *      polynomial coefficients for the 'x' direction.            *
c     *  fy--an array containing, in double precicion, the orthogonal  *
c     *      polynomial coefficients for the 'y' direction.            *
c     *  vx--an array containing the the 'x' coordinates.              *
c     *  vy--an array containing the the 'y' coordinates.              *
c     *  alk--an array containg the least-square coefficients.         *
c     *  zsq--an array containing the computed 'z-square' matrix.      *
c     *  dimt--the maximum  degree of the surface to be generated + 1. *
c     *  ny--the number of 'y' values.                                 *
c     *  mx--the number of 'x' values.                                 *
c     *  iopt   --an integer array used to convey program options.      *
c     *  sumz and sumzsq--sum of measured 'z' values and the sum of    *
c     *      the squares of 'z' respectively.                          *
c     *  ------------------------------------------------------------  *
c     *  called subprograms:                                           *
c     *    inputz(zm,i,j)--a user supplied routine to input the        *
c     *      measured 'z' values, row by row.                          *
c     *    real poly*8--a function to evaluate a polynomial with known *
c     *      coefficients.                                             *
c     ******************************************************************
c
c     declarations
	double precision fx,fy,vx,vy,alk,zsq,sumz,sumzsq,terml,
     1 fycom,fxycom,xnumer,xdenom,poly,zm,zmis,vyns,vxm
	integer dimt,dimp
	dimension fx(400),alk(20,20),fy(20,20),zsq(20,20),iopt(10),
     1 zm(1),vx(1),vy(1)
c     main pgm to generate the least-square coefficients-alk(l,k).
c
c     initialize necessary variables and arrays.
	do 10 l=1,dimt
	do 10 k=1,dimt
	alk(l,k)=0.0d0
10	zsq(l,k)=0.0d0
	sumz=0.0d0
	sumzsq=0.0d0
	lout=8
	if (iopt(5).eq.0)  go to 13
	assign 30 to ksw
	assign 83 to ksw2
	go to 17
13	assign 40 to ksw
	assign 85 to ksw2
c
c     start main loop--computations are by rows of 'z'.
c     user must input zm as real*8.
17	do 80 ns=1,ny
c     read in a row of data.
	call inputz(zm,ns,mx,insw)
	vyns=vy(ns)
c     sum z and z**2 for the row of data.
	do 20 is=1,mx
	zmis=zm(is)
	sumz=sumz+zmis
	sumzsq=zmis*zmis+sumzsq
20	continue
	do 70 l=1,dimt
        goto (30,40),ksw
c	go to ksw,(30,40)
30	dimp=dimt
	go to 50
40	dimp=dimt-l+1
c     compute  fl(y) for degree l-1 and row ns.
50	fycom=poly(fy,vyns,l,0)
	do 70 k=1,dimp
	xnumer=0.0d0
	xdenom=0.0d0
	do 60 m=1,mx
	vxm=vx(m)
c     compute fl(y)*fk(x)--fk(x) degree=k-1 for column m.
c	fxycom=poly(fx,vxm,k,0)*fycom
c  equivalant code
	nfx=(k-1)*20+k
	fxycom=0.d0
	do 55 i=1,k
	fxycom=fxycom*vxm+fx(nfx)
55	nfx=nfx-20
	fxycom=fxycom*fycom
	xnumer=zm(m)*fxycom+xnumer
	xdenom=fxycom*fxycom+xdenom
60	continue
c     temp use of alk and zsq for accum numerators and denomenators.
	alk(l,k)=xdenom+alk(l,k)
70	zsq(l,k)=xnumer+zsq(l,k)
80	continue
c     now compute values of alk and zsq
	do 89 l=1,dimt
	go to ksw2,(83,85)
83	dimp=dimt
	go to 87
85	dimp=dimt-l+1
87	do 89 k=1,dimp
	terml= zsq(l,k)/alk(l,k)
	zsq(l,k)=terml*zsq(l,k)
89	alk(l,k)=terml
c     output results
	write (lout,90)
90	format ('1',10x,'least square matrix'/' coefficients of ',
     1 'orthogonal polynomials in ascending degrees of f(x)',
     1 ' across and g(y) down'//)
	do 100 i=1,dimt
100	write (lout,110)  (alk(i,j),j=1,dimt)
110	format (10(1x,e12.5))
	write (lout,120)
120	format (///,10x,'z-square matrix'/' components of orthogonal',
     1 ' polynomials in ascending degrees of f(x)',
     1 ' across and g(y) down'//)
	do 130 i=1,dimt
130	write (lout,110)  (zsq(i,j),j=1,dimt)
	return
	end
	double precision function poly(a,w,k,l)
c     this function evaluates polynomials with coefficients in a 2-d
c     array.
	double precision a,w,sump
	dimension a(20,20)
	m=l
	sump=0.0d0
	if (l.eq.0)  m=k
	do 10 i=1,k
	is=k-i+1
10	sump=sump*w+a(m,is)
	poly=sump
	return
	end
	subroutine opoly(w,num,ip,fw)
c     this subroutine generates the coefficients for orthogonal
c     polynomials up to the 19th degree.
c     coded by--ronald r. wahl, special projects br, usgs, denver co
c     ******************************************************************
c     *  w-input array                                                 *
c     *  num-the length of the input array                             *
c     *  ip-the desired maximum degree of the polynomials+1--max=20    *
c     *  fw-the coefficient array, in double precision. in fw(i,j) the *
c     *     ith position is the location for the coefficients of the   *
c     *     polynomial of the ith-1 degree.                            *
c     *  called subprograms----real poly*8                             *
c     ******************************************************************
	double precision fw,a1,a11,a11sq,a2,alfa,b1,beta,poly,sumw,w,wl
	dimension w(1),fw(20,20)
	np=ip
c     initialize array for coefficients of polynomials--other variables
	if (np.gt.20)  np=20
	id=np-1
	sumw=0.0d0
	do 20 j=1,np
	do 10 i=1,np
10	fw(i,j)=0.0d0
20	fw(j,j)=1.0d0
	b1=dfloat(num)
	do 30 i=1,num
30	sumw=sumw+w(i)
	fw(2,1)=-sumw/b1
c     start main loop to gen coefficients
	do 60 k=2,id
	a1=0.0
	a2=0.0
	kp=k+1
	km=k-1
	do 40 l=1,num
	wl=w(l)
	a11=poly(fw,wl,k,0)
	a11sq=a11*a11
	a1=a11sq*wl+a1
40	a2=a11sq+a2
	alfa=-a1/a2
	beta=-a2/b1
	b1=a2
	fw(kp,1)= alfa*fw(k,1)+beta*fw(km,1)
	do 50 l=2,k
	lm=l-1
50	fw(kp,l)= alfa*fw(k,l)+beta*fw(km,l)+fw(k,lm)
60	continue
	return
	end
	subroutine gridr(nc,nr,zg,wz,ier)
c  initialize grid with reasonable anomalies
c  wz at least max(nc,nr) 
	dimension iw1(3)
	dimension zg(1),wz(1)
        data dval/1.e+38/,iw1/3,7,9/
c	data dval/'ffff7fff'x/,iw1/3,7,9/
	fltbig=dval*1.e-2
	ier=0
	nn=nc*nr
	nsep=(iw1(3)-1)/2
	ns1=nsep+1
	mxpass=max(nc,nr)/(2*nsep)
	i=1
201	if(i.gt.nn) return
	if(zg(i).gt.fltbig) go to 200
	i=i+1
	go to 201
c
c  insert control point net
c     distances for ring averages specified by iw1
200	if(nc.lt.nsep+ns1 .or. nr.lt.nsep+ns1) go to 140
	ipass=0
	iflag=1
130	iw=iw1(iflag)
	ihw=(iw-1)/2
	iset=ihw*nc+ihw
131	nass=0
	do 110 jj=ns1,nr-ihw,nsep
	ip=(jj-1)*nc+ns1
	do 100 ii=ns1,nc-ihw,nsep
	if(zg(ip).lt.fltbig) go to 100
	ips=ip-iset
	it=0
	t=0.0
	do 121 j=1,iw
	ip2=ips
	do 120 i=1,iw
	if(zg(ip2).gt.fltbig) go to 120
	t=zg(ip2)+t
	it=it+1
120	ip2=ip2+1
121	ips=ips+nc
	if(it-3.lt.0) go to 100
	zg(ip)=t/float(it)
	nass=nass+1
100	ip=ip+nsep
110	continue
	ipass=ipass+1
	if(nass.eq.0 .and. iflag.eq.3) go to 140
	if(ipass.gt.mxpass) go to 140
	if(iflag-2)133,132,131
133	iflag=2
	go to 130
132	iflag=3
	go to 130
140	continue
c	call gprint(1,1,nc,nr,zg,' ')
c
c  fill holes
	inc=nc*nsep
	j=inc+1
	if(nr.lt.ns1) go to 21
	do 20 irow=ns1,nr-nsep,nsep
	call plugm3(nc,zg(j),dval)
20	j=j+inc
21	do 28 icol=1,nc
	j=icol
	do 22 k=1,nr
	wz(k)=zg(j)
22	j=j+nc
	call plugm3(nr,wz,dval)
	j=icol
	do 24 k=1,nr
	zg(j)=wz(k)
24	j=j+nc
28	continue
c  final check
	do 40 i=1,nn
        if(zg(i).ge.dval) go to 41
40	continue
	return
41	t=0.0
	it=0
	do 42 i=1,nn
        if(zg(i).ge.dval) go to 42
	t=t+zg(i)
	it=it+1
42	continue
	if(it.eq.0) stop ' cannot init grid'
	t=t/float(it)
	write (*,43)t
43	format(' gridr init with',1pe15.5)
	do 44 i=1,nn
        if(zg(i).ge.dval) zg(i)=t
44	continue
	return
	end
	subroutine plugm3(n,z,dv)
c  plug holes using linear interpolation
	dimension z(n)
	fltbig=dv*.01
	do 1 is=1,n
	if(z(is) .lt. fltbig) go to 2
1	continue
	return
2	ix=is
3	ix=ix-1
	if(ix.lt.1) go to 4
	z(ix)=z(is)
	go to 3
4	do 5 idv=is,n
	if(z(idv) .gt. fltbig) go to 6
5	continue
	return
6	is=idv-1
	do 7 ie=idv,n
	if(z(ie) .lt. fltbig) go to 10
7	continue
	ix=is
9	ix=ix+1
	if(ix.gt.n) return
	z(ix)=z(is)
	go to 9
10	dz=(z(ie)-z(is))/float(ie-is)
	do 11 i=is+1,ie-1
11	z(i)=z(i-1)+dz
	is=ie
	go to 4
	end
	subroutine rowio(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end


