c              program aver2d
c  fast 2-d averaging lowpass filter.
c  input and output are the same size.
c  cutoff wavelength approx   1.653 * dx * nsamples_in_window.
c  first sidelobe down 15 dB.
c  m webring, usgs
c
	common /header/  id2,p2
	common /runave/  z(30000)
        character        p*8, p2*8, id*56, id2*56, f1*56, f2*56
	data             nwrk/30000/, dv/0.1701412e+39/, 
     1                isw,jsw,ksw,lsw/13,14,15,16/
        call pfinit('aver2d')
	fltbig = dv*.01


	write(*,'(a,$)') ' input grid :'
	read(5,37) f1
37	format(a56)
c        open( unit=isw, file=f1, status='old', form='unformatted',
c     1     share='denywr' )
        call gopen(isw,f1,'old','read',ierr)

	write(*,'(a,$)') ' smoothed output grid :'
	read(5,37) f2
c        open( unit=jsw, file=f2, status='unknown', form='unformatted' )
        call gopen(jsw,f2,'new','readwrite',ierr)

c        open( unit=ksw, file='filled.tmp', status='unknown',
c     1     form='unformatted' )
        call gopen(ksw,' ','scratch','readwrite',ierr)
c        open( unit=lsw, file='expand.tmp', status='unknown',
c     1     form='unformatted' )
        call gopen(lsw,' ','scratch','readwrite',ierr)
	write(*,'(a,$)')' enter title :'
	read(5,37) id2
c
      write(*,'(a,$)')' enter averaging radius in data units (km e.g.):'
	read(5,*) rad
c        read(isw) id, p, nc, nr, nz, xo, dx, yo, dy
        call gheader('r',isw,id,nc,nr,xo,dx,yo,dy,ierr)
	iw = 1 + 2 * int( abs( rad / dx ) + .5 )
	write(p2,40) iw
40	format( 'aver:', i3 )
	write(*,45)iw
45	format('    window width =',i4,' values')
	if ( nc*(iw+3) .gt. nwrk ) stop ' nc*(width+3) > 30k'
c
50	write(*,60)
60	format(/,' 0 = completely filled',/,
     1        ' 1 = dvals matching input   : '$)
	read(5,*) ifunc
	if ( ifunc .lt. 0 .or. ifunc .gt. 1 ) go to 50
	ifunc = - ( ifunc + 1 )
c
	ihw  = iw / 2
	ihw1 = ihw + 1
	nc2  = nc + 2 * ihw
	nr2  = nr + 2 * ihw
	xo2  = xo - ihw * dx
	yo2  = yo - ihw * dy
c        write(lsw) id2,p2,nc2,nr2,nz,xo2,dx,yo2,dy
        call gheader('w',lsw,id2,nc2,nr2,xo2,dx,yo2,dy,ierr)
	do 150 i = 1, 2 * nc2
	  z(i) = dv	
150      continue
	do 160 j = 1, ihw
c          call rowio( nc2, z(nc2+1), 0, lsw, lsw, ie )
        call grow('w',lsw,j,z(nc2+1),nc2,ie)
160     continue
	do 170 j = 1, nr
c          call rowio( nc, z(ihw1), -1, isw, isw, ie )
        call grow('r',isw,idum,z(ihw1),nc,ie)
c          call rowio(nc2,       z,  0, lsw, lsw, ie )
        call grow('w',lsw,idum,z,nc2,ie)
170     continue
	do 180 j = 1, ihw
c          call rowio( nc2, z(nc2+1), 0, lsw, lsw, ie )
        call grow('w',lsw,idum,z(nc2+1),nc2,ie)
180     continue
c
	call fill2( z, nwrk, dv, lsw, ksw )
c
	n2 = iw * nc2 + 1
	n3 = n2 + nc2
	n4 = n3 + nc2
c        print*,n2,n3,n4
	call rave2( ifunc, iw, nc2, nr2, z, z(n2), z(n3), z(n4),
     1           ksw, jsw, lsw )
c
c        close(isw)
        call gclose(isw,'keep')
c        close(jsw)
        call gclose(jsw,'keep')
c        close(ksw,status='delete')
        call gclose(ksw,'delete')
c        close(lsw,status='delete')
        call gclose(lsw,'delete')
c        stop
	end
c*******************************************************************************
	subroutine fill2( z[huge], nwrk, dv, isw, jsw )
	dimension xx(4), z(nwrk)
	character id*56, p*8
	fltbig = 1.0e29
c
c        rewind isw
c        rewind jsw
c        read(isw) id,p,nc,nr,nz,xx
        call gheader('r',isw,id,nc,nr,xx(1),xx(2),xx(3),xx(4),ierr)
c
	itest = nwrk - max( (nc+nc/10), (nr+nr/10) )
	do 300 jump = 1, min( nc, nr )
	  mc = ( nc + jump ) / jump + 1
	  mr = ( nr + jump ) / jump + 1
	  mg = mc * mr
	  if ( mg .lt. itest ) go to 20
300     continue
	stop ' fill2: work area too small'
c
20	i2 = 1
	do 310 j = 1, mr
	  do 320 i = 1, mc
	    z(i2) = dv
	    i2    = i2 + 1
320     continue
310     continue
c
	jrow = 1
	do 330 j = 1, nr, jump
c        call rowio( nc, z(mg+1), -1, isw, jsw, ie )
        call grow('r',isw,idum,z(mg+1),nc,ie)
	if ( ie .ne. 0 ) go to 100
	i2 = 1 + mc * ( jrow - 1 )
	do 340 i = mg+1, mg+nc, jump
	  z(i2)=z(i)
	  i2=i2+1
340     continue
	if ( jump .gt. 1 ) then
	  do 350 k = 1, jump - 1
c            read( isw, end=100 )
        call grow('r',isw,idum,xdum,1,ierr)
350     continue
        endif
	jrow = jrow + 1
330     continue
c
100	call gridr( mc, mr, z, z(mg+1), fltbig, ier )
	if ( ier .ne. 0 ) stop ' fill2: error in subroutine gridr'
c
	rjump = 1.0 / float( jump )
c        rewind isw
c        read(isw)
         call gheader('r',isw,id,nc,nr,xx(1),xx(2),xx(3),xx(4),ierr)
c        write(jsw) id, p, nc, nr, nz, xx
        call gheader('w',jsw,id,nc,nr,xx(1),xx(2),xx(3),xx(4),ierr)
c
	do 150 j = 1, nr
c        call rowio( nc, z(mg+1), -1, isw, jsw, ie )
        call grow('r',isw,idum,z(mg+1),nc,ie)
	if ( ie .ne. 0 ) go to 200
	iy  = 1 + ( j-1 ) / jump
	idy = mod( j-1, jump )
	if ( iy + 1 .gt. mr ) stop ' fill2: y ran off control grid'
	istart = mc * ( iy - 1 )
	i2     = mg + 1
c
	do 130 i = 1, nc
	if ( z(i2) .lt. fltbig ) go to 130
	ix  = 1 + (i-1) / jump
	idx = mod( i-1, jump )
	if ( ix + 1 .gt. mc ) stop ' fill2: x ran off control grid'
	ig = istart + ix
	if ( idx - idy .gt. 0 ) then
	  ig1 = ig + 1
	  dzx = z(ig1)    - z(ig)
	  dzy = z(ig1+mc) - z(ig1)
	  else
	  ig2 = ig + mc
	  dzx = z(ig2+1) - z(ig2)
	  dzy = z(ig2)   - z(ig)
	endif
	z(i2) = z(ig) + (dzx*float(idx) + dzy*float(idy)) * rjump
130	i2 = i2 + 1
c
c        call rowio( nc, z(mg+1), 0, isw, jsw, ie )
        call grow('w',jsw,idum,z(mg+1),nc,ie)
150	continue
c
	return
200	write(*,210) j
210	format(' fill2: eof on input grid')
	return
	end
c*******************************************************************************
	subroutine rave2( ifunc, iw, nc, nr, z[huge], csum[huge],
     1  arow[huge], comp[huge], isw, jsw, ksw )
c  in-place 2d averaging
c  when ifunc < 0, trim grid one half window width. 
	common /header/ id2,p2
	dimension       z(nc,iw), csum(nc), arow(nc), comp(nc)
	character       id*56, id2*56, p*8, p2*8
	logical         trim
c
	trim = .false.
	if ( ifunc .lt. 0 ) trim = .true.
	ifunc = iabs( ifunc ) - 1
	if ( ifunc .lt. 0 .or. ifunc .gt. 2 ) 
     1    stop ' rave2: 0 <= ifunc <= 2'
	if ( nc.le.iw .or. nr.le.iw) 
     1    stop ' rave2: window larger than grid dimension'
	ihw  = ( iw - 1 ) / 2
	ihw1 = ihw + 1
	if ( ihw + ihw1 .ne. iw ) stop ' rave2: iw must be odd'
c
c        rewind isw
c        rewind jsw
c        rewind ksw
c        read(isw) id, p, nc, nr, nz, xo, dx, yo, dy
        call gheader('r',isw,id,nc,nr,xo,dx,yo,dy,ierr)
	if ( trim ) then
	  xo2 = xo + ihw * dx
	  yo2 = yo + ihw * dy
	  nc2 = nc - 2 * ihw
	  nr2 = nr - 2 * ihw
c          write(jsw) id2, p2, nc2, nr2, nz, xo2, dx, yo2, dy
          call gheader('w',jsw,id2,nc2,nr2,xo2,dx,yo2,dy,ierr)
	  else
c          write(jsw) id2, p2, nc, nr, nz, xo, dx, yo, dy
          call gheader('w',jsw,id2,nc,nr,xo,dx,yo,dy,ierr)
	endif
c        read(ksw)
        call gheader('r',ksw,id,nc,nr,xo,dx,yo,dy,ierr)
c
c  initialize column sums
	do 301 i = 1, nc
	  csum(i) = 0.
301     continue
	do 305 j = 1, ihw
c          call rowio( nc, z(1,j), -1, isw, jsw, iend )
          call grow('r',isw,idum,z(1,j),nc,iend)
	  do 320 i = 1, nc
	    csum(i) = csum(i) + z(i,j)
320     continue
305     continue
c
c  lower half window
	irow=0
	do 50 j=1,ihw1
	jrow=j+ihw
c        call rowio(nc,z(1,jrow),-1,isw,jsw,iend)
        call grow('r',isw,idum,z(1,jrow),nc,iend)
	do 51 i = 1, nc
	  csum(i) = csum(i) + z(i,jrow)
51      continue
c        read(ksw) ydum, comp
        call grow('r',ksw,idum,comp,nc,ie)
	call raopr( nc, jrow, iw, csum, arow, comp, ifunc )
	if ( .not. trim ) then
c          write(jsw) ydum,arow
          call grow('w',jsw,idum,arow,nc,ie)
	  irow=irow+1
	endif
50	continue
	if ( trim ) then
c          call rowio( nc2, arow(ihw1), 0, jsw, jsw, ie )
          call grow('w',jsw,idum,arow(ihw1),nc2,ie)
	  irow = irow + 1
	endif
c
c  full window
c  iptr is the trailing row address
	iptr=1
	jrow=iw
	do 100 j=ihw1+1,nr-ihw
	do 340 i = 1, nc
	  csum(i) = csum(i) - z(i,iptr)
340     continue
c        call rowio(nc,z(1,iptr),-1,isw,jsw,iend)
        call grow('r',isw,idum,z(1,iptr),nc,iend)
	do 350 i = 1, nc
	  csum(i) = csum(i) + z(i,iptr)
350     continue
	iptr = iptr + 1
	if ( iptr .gt. iw ) iptr = 1
c        read(ksw) ydum,comp
        call grow('r',ksw,idum,comp,nc,ie)
	call raopr(nc,jrow,iw,csum,arow,comp,ifunc)
	if(trim) then
c         call rowio(nc2,arow(ihw1),0,jsw,jsw,ie)
         call grow('w',jsw,idum,arow(ihw1),nc2,ie)
	 irow=irow+1
	 else
c         write(jsw) ydum,arow
         call grow('w',jsw,idum,arow,nc,ie)
	 irow=irow+1
	endif
100	continue
	if(trim) go to 300
c
c  upper half window
	do 200 j=nr-ihw+1,nr
	jrow=jrow-1
	do 360 i = 1, nc
	  csum(i) = csum(i) - z(i,iptr)
360     continue
	iptr=iptr+1
	if(iptr.gt.iw) iptr=1
c        read(ksw) ydum,comp
        call grow('r',ksw,idum,comp,nc,ie)
	call raopr(nc,jrow,iw,csum,arow,comp,ifunc)
c        write(jsw) ydum,arow
        call grow('w',jsw,idum,arow,nc,ie)
	irow=irow+1
200	continue
	if(irow.ne.nr) write(*,310)irow,nr
	return
300	if(irow.ne.nr2) write(*,310)irow,nr2
310	format(' rave2: rows output vs rows desired',2i5)
	return
	end
c*******************************************************************************
	subroutine raopr(nc,jrow,iw,csum[huge],arow[huge],
     1 comp[huge],ifunc)
c  ifunc=0 smooth all, =1 smooth & mask, =2 smooth extrapolated
	dimension csum(nc),arow(nc),comp(nc)
	data dv/0.1701412e+39/
	ihw=(iw-1)/2
	ihw1=ihw+1
	t=0.0
	do 10 i=1,ihw
10	t=t+csum(i)
	do 20 i=1,ihw1
	t=t+csum(ihw+i)
20	arow(i)=t/float((ihw+i)*jrow)
c
	rwt=1.0/float(iw*jrow)
	ldr=iw+1
	itrl=1
	do 50 i=ihw1+1,nc-ihw
	t=(t+csum(ldr))-csum(itrl)
	arow(i)=t*rwt
	ldr=ldr+1
50	itrl=itrl+1
c
	icol=iw-1
	do 100 i=nc-ihw+1,nc
	t=t-csum(itrl)
	arow(i)=t/float(icol*jrow)
	itrl=itrl+1
100	icol=icol-1
c
	if(ifunc.eq.0) return
	do 200 i=1,nc
	go to (110,120) ifunc
110	if(comp(i).gt.1.e37) arow(i)=dv
	go to 200
120	if(comp(i).lt.1.e37) arow(i)=comp(i)
200	continue
	return
	end
c******************************************************************************
c        subroutine rowio(n,z[huge],iop,idev,jdev,iend)
cc  read iop<0, write iop=0, r&w iop=1
c        dimension z(n)
c        iend=0
c        if(iop)1,2,1
c1       read(idev,end=10) y,z
c        if(iop)9,9,2
c2       write(jdev) y,z
c9       return
c10      iend=1
c        return
c        end
c******************************************************************************
	subroutine gridr(nc,nr,zg[huge],wz[huge],dval,ier)
c  quickly replace dvals with reasonable values.
c  A control net spaced at 'nsep' is generated from ring averages
c  with radii 1, 3, nsep; linear interpolation completes the process.
c  Array 'wz' is at least max(nc,nr).
	dimension iw1(3)
	dimension zg(1),wz(1)
	data iw1/3,7,9/
	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).ge.dval) 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.dval) 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).ge.dval) 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
c  fill holes
	inc=nc*nsep
	j=inc+1
	if(nr.lt.ns1) go to 21
	do 20 irow=ns1,nr-nsep,nsep
	call plugm(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 plugm(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
c******************************************************************************
	subroutine plugm(n,z[huge],dv)
c  plug holes using linear interpolation
	dimension z(n)
	do 1 is=1,n
	if(z(is) .lt. dv) 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) .ge. dv) go to 6
5	continue
	return
6	is=idv-1
	do 7 ie=idv,n
	if(z(ie) .lt. dv) 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
