
c  program megaplug 

	parameter ( nwork=100000 )
	dimension iwork(nwork), work(nwork)
        equivalence(iwork,work)
	character ifile*56, ofile*56, title*56, pgm*8
	data idv, jdv, kdv, ldv/10,11,12,13/, dval/0.1701412e39/

      nwork2 = nwork / 2

1	write(0,'(a,$)') ' enter grid filename:'
	read(0,37) ifile
37	format( a )
	open( idv, file=ifile, status='old', form='unformatted',
     1     share='denywr')
c       replace following with an err=
c	if ( ios .ne. 0 ) then
c	  type *,' file not found'
c	  go to 1
c	endif
	call rdhdr( idv, title, pgm, nc, nr, nz, xo, dx, yo, dy, 
     1           iproj, cm, baslat )

	nn = nc * nr
	if ( nn .gt. nwork2 ) then
        write(0,*) ' ncol*nrow >',nwork2
	  stop
	endif

	write(0,'(a,$)') ' enter output filename : '
	read(0,37) ofile
	length = nc + 1
	if ( length .lt. 30 ) length = 30

	write(0,'(a,$)')
     1 ' number of minimum curvature iterations to use :'
	read(0,*) nim

	write(0,*) ' enter 1 to completely fill'
	write(0,*) '    or 2 generate hulls around contiguous data'
	read(0,*) ifunc
	
	open( jdv, file=ofile, status='unknown', form='unformatted',
     1     recl=length*4)     

	write( jdv ) title, pgm, nc, nr, nz, xo, dx, yo, dy,
     1                    iproj, cm, baslat 
	n2     = 1 + nwork2
	write(0,*) dval, idv, jdv
	write(0,*) nc, nr, nim, nwork2, n2
	call megafill( dval, idv, jdv, nc, nr, nim,
     1              work, iwork(n2), nwork2 )

	if ( ifunc .ne. 2 ) stop

	close( jdv ) 
c       copy output to direct access file
	open( jdv, file=ofile, status='old', form='unformatted',
     1     recl=length*4 )  
        open (ldv,file='mega.tmp',status='unknown',
     1  form='unformatted',access='direct',recl=length*4)
        read(jdv)title,pgm,nc,nr,nz,xo,dx,yp,dy,iproj,cm baslat
        write(ldv,rec=1)title,pgm,nc,nr,nz,xo,dx,yp,dy,iproj,cm,baslat
        do 2000 i=1,nr
        call rowio(nc,work,-1,jdv,0,iend)
        call rowda(nc,iwork,0,i,ldv,ierr)
        if (ierr .ne. 0) then
                write(0,'(1x,a,i4)')
     1  'Error writing temp file, mega.tmp. Error is ',ierr
                stop
        endif
2000    continue
c        close (jdv)

	open( kdv ,file='masking.tmp',status='unknown',form='unformatted',
     1     access='direct', recl=length*4 )
	write( kdv, rec=1 ) title, pgm, nc, nr, nz, xo, dx, yo, dy,
     1                    iproj, cm, baslat 
	call genmsk( idv, kdv, nc, nr, work, iwork(nc+1) )

c  set outside polygon to 2.
	ihwind = 0
	nwork3 = 50000
	call outside( kdv, ihwind, iwork, nwork3 )

	istart = 1
	iflag  = 2
	call templat( istart, istart, iflag, dval, 
     1             kdv,     iwork, nc, nr,
     1             ldv, work(n2), nc, nr )

        close (kdv,status='delete')
    
        rewind (jdv)
        read(ldv,rec=1)title,pgm,nc,nr,nz,xo,dx,yo,dy,iproj,cm,baslat
        write(jdv)title,pgm,nc,nr,nz,xo,dx,yo,dy,iproj,cm,baslat
        do 1000 i=1,nr
        call rowda(nc,iwork,1,i,ldv,ierr) 
        if(ierr.ne.0) then
        write(0,'(1x,a,i4)')
     1  'Error reading temp file, mega.tmp. Error is ',ierr
        stop
        endif
        call rowio(nc,work,0,0,jdv,iend)       
1000    continue
        close(ldv,status='delete')      
        close(jdv)
	stop
	end
c*******************************************************************************
	subroutine genmsk( iswt, kswt, nc, nr, work[huge], mask[huge] )
	dimension work(nc), mask(nc)
	rewind( iswt )
	read( iswt )
	do 10 jrec = 2, nr + 1
	  read( iswt ) y, work
	  do 20 i = 1, nc
	    mask(i) = 0
	    if ( work(i) .lt. 1.0e30 ) mask(i) = 1
20      continue
	  write( kswt, rec=jrec ) jrec, mask
10      continue
	return
	end
c******************************************************************************
	subroutine megafill( dval, idv, jdv, nc, nr, nim,
     1                    zg[huge], iqd[huge], nwork )
c  plug holes using minimum curvature interpolation
c  gridr and curvmn are from the minc program, usgs open file 81-1224

	dimension wz(2000), zg(nwork), iqd(nwork)

	nn = nc * nr
	do 10 i = 1, nn
	  iqd(i) = -1
10      continue

	ndx = 1
	do 20 j = 1, nr
	call rowio( nc, zg(ndx), -1, idv, jdv, iend )
	if ( iend .eq. 1 ) go to 99
	i2 = ndx
	do 30 i = 1, nc
	  if ( zg(i2) .gt. 1.e29 ) then
	    zg(i2)  = dval
	    iqd(i2) = 0
	  endif
	 i2 = i2 + 1
30      continue
	ndx = ndx + nc
20	continue

	call gridr( nc, nr, zg, wz, dval, ier )
	if ( ier .eq. 1 ) then
	  write(0,*) ' gridr: initialization error'
	  stop
	endif
	if ( nim .gt. 0 ) then
	  eps = 0.
	  call curvmn( zg, iqd, wz, nc, nr, eps, nim, a, b, ni )
	endif

	ndx = 1
	do 40 i = 1, nr
	  call rowio( nc, zg(ndx), 0, idv, jdv, ie )
	  ndx = ndx + nc
40      continue
	return

99	write(0,*) ' %%megautp: premature EOF'
	return
	end
c******************************************************************************
	subroutine rowio(n,z[huge],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
c******************************************************************************
	subroutine gridr(nc,nr,zg[huge],wz,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(0,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
c******************************************************************************
	subroutine curvmn(zg[huge],iqd[huge],b,nc,nr,epsmx,
     1 nim,eps1,dn1,ni)
c   applies minimum curvature equations to the first
c  nc*nr elements of array zg.
c   array iqd contains nc*nr elements which indicate
c  for each mesh location the quadrant where a data
c  value is located. an iqd value of zero indicates
c  no data and -1 locks the present mesh value.
c   array b should contain 6*nc*nr elements used for
c  weighting when iqd is 1 to 4, in the case where
c  iqd is only 0 or -1, b can be of length one.
c   the over-relaxation parameter w increases
c  as the system converges until 1.7 is reached.
	dimension zg(1),iqd(1),b(1)
	data nimn/5/,lmtc/1/
	if(nc.lt.5 .or. nr.lt.5) return
	ni=0
	dn=1.e20
	w=1.3
	eps=0.
	eps1=0.
	epsm=abs(epsmx)
111	continue
	if(ni.ge.nim) go to 72
	eps=0.
c first row
	if(iqd(1))2,1,1
1	zg(1)=(( (2.*(zg(2)+zg(nc+1))-zg(nc+nc+1)-zg(3))*.5 )-
     1 zg(1))*w+zg(1)
2	j1=nc+2
	j2=j1+nc
	if(iqd(2))4,3,3
3	zg(2)=(( (4.*(zg(3)+zg(j1))+2.*zg(1)-zg(4)-zg(j1-1)-
     1 zg(j1+1)-zg(j2))*.16666667 )-zg(2))*w+zg(2)
4	do 6 i=3,nc-2
	j1=i+nc
	j2=j1+nc
	if(iqd(i))6,5,5
5	zg(i)=(( (4.*(zg(i-1)+zg(j1)+zg(i+1))-zg(j1+1)-zg(j1-1)-
     1 zg(j2)-zg(i+2)-zg(i-2))*.14285714 )-zg(i))*w+zg(i)
6	continue
	if(iqd(nc-1))8,7,7
7	i=nc-1
	j1=i+nc
	zg(i)=(( (4.*(zg(i-1)+zg(j1))+2.*zg(i+1)-zg(i-2)-
     1 zg(j1+1)-zg(j1-1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
8	if(iqd(nc))10,9,9
9	j1=nc+nc
	zg(nc)=(( (2.*(zg(j1)+zg(nc-1))-zg(nc-2)-zg(j1+nc))*.5 )-
     1 zg(nc))*w+zg(nc)
c second row
10	if(iqd(nc+1))12,11,11
11	i=nc+1
	j1=i+nc
	zg(i)=(( (4.*(zg(j1)+zg(i+1))+2.*zg(1)-zg(2)-
     1 zg(i+2)-zg(j1+1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
12	if(iqd(nc+2))14,13,13
13	i=nc+2
	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i+1))+4.*(zg(jm)+zg(i-1))-
     1 2.*zg(j1+1)-zg(jm+1)-zg(j1-1)-zg(i+2)-zg(j1+nc))*
     1 5.5555556e-2 )-zg(i))*w+zg(i)
14	do 16 i=nc+3,nc+nc-2
	j1=i+nc
	jm=i-nc
	if(iqd(i))16,15,15
15	zg(i)=(( (8.*(zg(i-1)+zg(j1)+zg(i+1))+4.*(zg(jm))-
     1 2.*(zg(j1-1)+zg(j1+1))-zg(jm-1)-zg(jm+1)-
     1 zg(j1+nc)-zg(i+2)-zg(i-2))*5.263158e-2 )-zg(i))*w+zg(i)
16	continue
	i=nc+nc-1
	if(iqd(i))18,17,17
17	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i-1))+4.*(zg(jm)+zg(i+1))-2.*zg(j1-1)-
     1 zg(jm-1)-zg(j1+1)-zg(i-2)-zg(j1+nc))*5.5555556e-2 )-
     1 zg(i))*w+zg(i)
18	i=nc+nc
	if(iqd(i))20,19,19
19	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(j1)+zg(i-1))+2.*zg(jm)-zg(jm-1)-
     1 zg(i-2)-zg(j1-1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
c rows 3 to nr-2
20	do 39 j=3,nr-2
	i=(j-1)*nc+1
	if(iqd(i))22,21,21
21	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(i+1)+zg(j1)+zg(jm))-zg(j1+nc)-zg(j1+1)-zg(i+2)-
     1 zg(jm+1)-zg(jm-nc))*.14285714 )-zg(i))*w+zg(i)
22	i=i+1
	if(iqd(i))24,23,23
23	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i+1)+zg(jm))+4.*zg(i-1)
     1 -2.*(zg(j1+1)+zg(jm+1))-zg(j1-1)-zg(j1+nc)-zg(i+2)-
     1 zg(jm-nc)-zg(jm-1))*5.2631578e-2 )-zg(i))*w+zg(i)
24	do 35 j2=3,nc-2
	i=i+1
	if(iqd(i))35,25,25
25	j1=i+nc
	jm=i-nc
	d=zg(i)
	if(iqd(i))26,26,27
26	d=(( (8.*(zg(i+1)+zg(i-1)+zg(jm)+zg(j1))-2.*(zg(j1+1)+zg(jm+1)+
     1 zg(jm-1)+zg(j1-1))-zg(j1+nc)-zg(jm-nc)-zg(i-2)-zg(i+2))*
     1 .05 )-d)*w+d
	go to 33
27	ndx=(i-1)*6+1
	b1=b(ndx)
	b2=b(ndx+1)
	b3=b(ndx+2)
	b4=b(ndx+3)
	b5=b(ndx+4)
	b6=b(ndx+5)
	go to (28,29,30,31)iqd(i)
28	bu=b1*zg(jm+1)+b2*zg(jm)+b3*zg(i-1)+b4*zg(j1-1)
	go to 32
29	bu=b1*zg(jm-1)+b2*zg(jm)+b3*zg(i+1)+b4*zg(j1+1)
	go to 32
30	bu=b1*zg(j1-1)+b2*zg(j1)+b3*zg(i+1)+b4*zg(jm+1)
	go to 32
31	bu=b1*zg(j1+1)+b2*zg(j1)+b3*zg(i-1)+b4*zg(jm-1)
32	t=.25*(zg(j1+nc)+zg(i-2)+zg(jm-nc)+zg(i+2))
     1 +.5*(zg(j1-1)+zg(jm-1)+zg(jm+1)+zg(j1+1))-
     1 (zg(j1)+zg(i-1)+zg(jm)+zg(i+1))
	d=(( (bu+b5-t)*b6 )-d)*w+d
33	epsln=d-zg(i)
	if(abs(epsln).lt.abs(eps)) go to 34
	eps=epsln
	ieps=i
34	zg(i)=d
35	continue
	i=i+1
	if(iqd(i))37,36,36
36	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i-1)+zg(jm))+4.*zg(i+1)-2.*(zg(j1-1)+
     1 zg(jm-1))-zg(jm+1)-zg(jm-nc)-zg(i-2)-
     1 zg(j1+nc)-zg(j1+1))*5.2631578e-2 )-zg(i))*w+zg(i)
37	i=i+1
	if(iqd(i))39,38,38
38	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(j1)+zg(i-1)+zg(jm))-zg(jm-nc)-zg(jm-1)-zg(i-2)-
     1 zg(j1-1)-zg(j1+nc))*.14285714 )-zg(i))*w+zg(i)
39	continue
c row nr-1
40	i=(nr-2)*nc+1
	if(iqd(i))42,41,41
41	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(jm)+zg(i+1))+2.*zg(j1)-zg(jm-nc)-zg(jm+1)-
     1 zg(i+2)-zg(j1+1))*.16666667 )-zg(i))*w+zg(i)
42	i=i+1
	if(iqd(i))44,43,43
43	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(i+1)+zg(jm))+4.*(zg(i-1)+zg(j1))-
     1 2.*zg(jm+1)-zg(jm-1)-zg(jm-nc)-zg(i+2)-
     1 zg(j1+1))*5.5555556e-2 )-zg(i))*w+zg(i)
44	do 46 j=3,nc-2
	i=i+1
	if(iqd(i))46,45,45
45	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(i-1)+zg(jm)+zg(i+1))+4.*zg(j1)-
     1 2.*(zg(jm-1)+zg(jm+1))-zg(j1-1)-zg(i-2)-
     1 zg(jm-nc)-zg(i+2)-zg(j1+1))*5.2631578e-2 )-zg(i))*w+zg(i)
46	continue
	i=(nr-1)*nc-1
	if(iqd(i))48,47,47
47	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(i-1)+zg(jm))+4.*(zg(j1)+zg(i+1))-2.*zg(jm-1)-
     1 zg(j1-1)-zg(i-2)-zg(jm-nc)-zg(jm+1))*5.5555556e-2 )-
     1 zg(i))*w+zg(i)
48	i=i+1
	if(iqd(i))50,49,49
49	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(i-1)+zg(jm))+2.*zg(j1)-zg(jm-nc)-zg(jm-1)-
     1 zg(i-2)-zg(j1-1))*.16666667 )-zg(i))*w+zg(i)
c last row
50	i=i+1
	if(iqd(i))52,51,51
51	jm=i-nc
	zg(i)=(( (2.*(zg(i+1)+zg(jm))-zg(i+2)-zg(jm-nc))*.5 )-
     1 zg(i))*w+zg(i)
52	i=i+1
	if(iqd(i))54,53,53
53	jm=i-nc
	zg(i)=(( (4.*(zg(i+1)+zg(jm))+2.*zg(i-1)-zg(i+2)-zg(jm+1)-
     1 zg(jm-nc)-zg(jm-1))*.16666667 )-zg(i))*w+zg(i)
54	do 56 j=3,nc-2
	i=i+1
	if(iqd(i))56,55,55
55	jm=i-nc
	zg(i)=(( (4.*(zg(i-1)+zg(i+1)+zg(jm))-zg(i-2)-zg(jm-1)-
     1 zg(jm-nc)-zg(jm+1)-zg(i+2))*.14285714 )-zg(i))*w+zg(i)
56	continue
	i=i+1
	if(iqd(i))58,57,57
57	jm=i-nc
	zg(i)=(( (4.*(zg(i-1)+zg(jm))+2.*zg(i+1)-zg(i-2)-
     1 zg(jm-1)-zg(jm-nc)-zg(jm+1))*.16666667 )-zg(i))*w+zg(i)
58	i=i+1
	if(iqd(i))60,59,59
59	jm=i-nc
	zg(i)=(( (2.*(zg(i-1)+zg(jm))-zg(i-2)-zg(jm-nc))*.5 )-
     1 zg(i))*w+zg(i)
60	if(ni)70,70,71
70	eps1=abs(eps/w)
71	ni=ni+1
	if(eps.eq.0) go to 72
	dn1=abs(eps/w)
	if(dn1.le.epsm .and. ni.ge.nimn) go to 72
	dlam=dn1/dn
	dn=dn1
	if(dlam.gt.1.) go to 74
	if(dlam.lt..8) go to 75
	if(w.ge.1.6) go to 75
	w=w+.1
	go to 75
74	if(iconv.eq.lmtc) go to 76
	iconv=iconv+1
	go to 75
76	w=w-.1*aint(dlam*10.-9.11)
	iconv=0
	if(w.lt.1.)w=1.
75	continue
	go to 111
72	return
	end
c*******************************************************************************
	subroutine outside( imask, ihwind, iwork[huge], nwork )

c  fill the polygons of no-data that intersect the grid boundary.

	dimension iwork(nwork)
	character title*56, pgm*8

	read( imask, rec=1 ) title, pgm, ncm, nrm
	call setpoly( imask, ihwind, iwork, ncm, nrm )

	nrblk = nwork / ncm
	call fill( imask, iwork, ncm, nrm, nrblk )

	return
	end
c*******************************************************************************
	subroutine setpoly( iswt, ihwind, m[huge], nc, nr )
	dimension m(nc)
	character title*56, pgm*8

c  set no-data sections on boundary to seed of 2.
c  flag half window section with 3.

c  col&row locations for the boundary of the output grid.
	icolm = ihwind + 1
	irowm = ihwind + 1
	icolx = nc - ihwind
	irowx = nr - ihwind

c  bottom half window
	do 300 j = 1, irowm - 1
	  irec = j + 1
	  read( iswt, rec=irec ) iy, m
	  do 310 i = 1, nc
	    m(i) = 3
310    continue
	  write( iswt, rec=irec ) iy, m
300     continue

c  middle section 

	do 320 j = irowm, irowx
	  irec = j + 1
	  read( iswt, rec=irec ) iy, m
c  left&right half window
	  i2 = nc
	  do 330 i = 1, ihwind
	    m(i)  = 3
	    m(i2) = 3
	    i2    = i2 - 1
330     continue
c  enter seed values on output grid boundary.
	  if ( m(icolm) .eq. 0 ) m(icolm) = 2
	  if ( m(icolx) .eq. 0 ) m(icolx) = 2
	  if ( j .eq. irowm  .or.  j .eq. irowx ) then
	    do 340 i = icolm, icolx
	      if ( m(i) .eq. 0 ) m(i) = 2
340     continue
	  endif
	  write( iswt, rec=irec ) iy, m
320     continue

c  top half window
	do 20 j = irowx + 1, nr
	  irec = j + 1
	  read( iswt, rec=irec, err=200 ) iy, m
	  do 10 i = 1, nc
	    m(i) = 3
10      continue
	  write( iswt, rec=irec ) iy, m
20      continue

	return
200	read( iswt, rec=1 ) title, pgm, ncg, nrg
	write(0,*) ' %%setpoly: EOF at row', j, ' of', nrg
	return
	end
c*******************************************************************************
	subroutine fill( iswt, m[huge], nc, nrmax, nrbuf )
c  seed fill driver. 
	dimension js(20), je(20), m(nc,nrbuf)

c  setup buffer dimensions
	jbuf  = 1
	jdir  = 1
	nbuf  = 1
	js(1) = 1
	je(1) = nrmax
	if ( nrmax .le. nrbuf ) go to 20

c  each buffer will share its ending row with the start of the next.
	nbuf = nrmax / nrbuf + 1
10	nr   = ( nrmax + nbuf ) / nbuf + 1
	je(1) = nr
	do 30 i = 2, nbuf
	  js(i) = je(i-1)
	  je(i) = js(i) + nr - 1
	  if ( je(i) .ge. nrmax ) then
	    je(i) = nrmax
	    if ( i .ne. nbuf ) nbuf = i
	    go to 20
	  endif
30      continue
c  increase number of buffers
	nbuf = nbuf + 1
	go to 10  

c  cycle buffers

20	nr   = je(jbuf) - js(jbuf) + 1
	jrow = js(jbuf)
	do 40 j = 1, nr
	  call rowda( nc, m(1,j), -1, jrow, iswt, ierr )
	  jrow = jrow + 1
40      continue

	call fillblk( m, nc, nr, nfill ) 
	ndir = ndir + nfill

	if ( nfill .gt. 0 ) then
	  jrow = js(jbuf)
	  do 50 j = 1, nr
	    call rowda( nc, m(1,j), 0, jrow, iswt, ierr )
	    jrow = jrow + 1
50      continue
	endif

	if ( nbuf .eq. 1 ) go to 999
	if ( ( jdir .gt. 0  .and.  jbuf .eq. nbuf )  .or. 
     1    ( jdir .lt. 0  .and.  jbuf .eq. 1 ) )   then
c  pass completed 
	  if ( ndir .eq. 0 ) then
	    if ( noass .eq. 1 ) go to 999
	    noass = 1
	    else
	    noass = 0
	  endif
	  jdir = -jdir
	  ndir = 0
	endif
	
	jbuf = jbuf + jdir
	go to 20

999	return
	end
c*******************************************************************************
	subroutine fillblk( m[huge], nc, nr, nfill )
c  seed fill controller for one block.
c  at least one element of m should equal 2.
	common /filcom/ jcurr, jlast
	dimension       m(nc,nr)
	logical         nearby

	nfill = 0
	noass = 0
	npass = 0

	jmin  = 1
	jmax  = nr
	jmina =  9999
	jmaxa = -9999
	nearby = .false.

10	npass = npass + 1
	ndir  = 0

c  scan up
	jcurr = jmin
c	do while ( nearby  .or.  jcurr .le. jmax + 1 )
15      if(.not.nearby .and. jcurr.gt.jmax+1) goto 20
	  if ( jcurr .gt. nr ) go to 20
	  if ( jcurr .lt. 2  ) jcurr = 2
	  nass  = 0
	  jlast = jcurr - 1
	  call fillrow( nc, m(1,jcurr), m(1,jlast), nass )
	  if ( nass .gt. 0 ) then
	    nearby = .true. 
	    ndir   = ndir + nass
	    if ( jcurr .gt. jmaxa ) jmaxa = jcurr
	    else
	    nearby = .false.
	  endif
	  jcurr = jcurr + 1
c	enddo
        goto 15

c  scan down
20	jcurr = jmax
c	do while ( nearby  .or.  jcurr .ge. jmin - 1 )
25     if(.not.nearby .and. jcurr.lt.jmin-1) goto 30
	  if ( jcurr .lt. 1 ) go to 30
	  if ( jcurr .gt. nr - 1 ) jcurr = nr - 1
	  nass  = 0
	  jlast = jcurr + 1
	  call fillrow( nc, m(1,jcurr), m(1,jlast), nass )
	  if ( nass .gt. 0 ) then
	    nearby = .true. 
	    ndir   = ndir + nass
	    if ( jcurr .lt. jmina ) jmina = jcurr
	    else
	    nearby = .false.
	  endif
	  jcurr = jcurr - 1
c	enddo
        goto 25

30	nfill = nfill + ndir

	if ( ndir .eq. 0 ) then
	  if ( noass .eq. 1 ) go to 999
	  noass = 1
	  jmin  = 1
	  jmax  = nr
	  else
	  noass = 0
	  jmin  = jmina
	  jmax  = jmaxa
	endif
	go to 10

999	return
	end
c*******************************************************************************
	subroutine fillrow( nc, mcurr[huge], mlast[huge], nass ) 
c  scan row for seed values.
	common /filcom/ jcurr, jlast
	dimension mcurr(nc), mlast(nc)

	do 100 icurr = 2, nc - 1
	  if ( mcurr(icurr) .eq. 0 ) then

	    if ( mcurr(icurr-1) .eq. 2 ) then
	      mcurr(icurr) = 2
	      nass         = nass + 1
	      go to 100
	    endif

	    if ( mlast(icurr) .eq. 2 ) then
	      mcurr(icurr) = 2
	      nass         = nass + 1
	      if ( mcurr(icurr-1) .eq. 0 ) 
     1          call backfill( mcurr, icurr-1, nass )
	      go to 100
	    endif

	    if ( mcurr(icurr+1) .eq. 2 ) 
     1        call backfill( mcurr, icurr, nass )
	  endif
100	continue

	if ( mcurr(1) .eq. 0  .and. 
     1  ( mcurr(2) .eq. 2  .or.  mlast(1) .eq. 2 ) ) then
	  mcurr(1) = 2
	  nass     = nass + 1
	endif

	if ( mcurr(nc)   .eq. 0  .and.
     1  ( mcurr(nc-1) .eq. 2  .or.  mlast(nc) .eq. 2 ) ) then
	  mcurr(nc) = 2
	  nass      = nass + 1
	endif

	return
	end
c*******************************************************************************
        subroutine backfill( mcurr[huge], icurr, nass )
	dimension mcurr(icurr)
	do 10 i = icurr, 1, -1
	  if( mcurr(i) .ne. 0 ) return
	  mcurr(i) = 2
	  nass     = nass + 1
10      continue
	return
	end
c*******************************************************************************
	subroutine templat( istart, jstart, iflag, dval,
     1                   iswt, m[huge], ncm, nrm,
     1                   jswt, z[huge], nc,  nr   )

c iswt contains the template where dval's are inserted into jswt.
c jswt area must be a subset of iswt.

	dimension m(ncm), z(nc)
	if ( istart .le. 0 ) go to 200
	if ( jstart .le. 0 ) go to 200
	ioff = istart - 1
	joff = jstart - 1
	iend = ioff + nc
	jend = joff + nr
	if ( iend .gt. ncm ) go to 200
	if ( jend .gt. nrm ) go to 200

	do 10 j = 1, nr
c  grid header is record 1, row 1 is record 2.
	  irec = ( j + 1 ) + joff
	  jrec = ( j + 1 )
	  read( iswt, rec=irec ) iy, m
	  read( jswt, rec=jrec ) y,  z

	  im = 1 + ioff
	  do 20 i = 1, nc
	    if ( m(im) .eq. iflag ) z(i) = dval
	    im = im + 1
20      continue

	  write( jswt, rec=jrec ) y, z
10      continue

	return
200	write(0,*) ' %%templat: grid is not a subset of the template.'
	return
	end
c*******************************************************************************
	subroutine rdhdr( iunit, id, pgm, nc,nr,nz, xo,dx,yo,dy,
     1                        iproj, cmerid, baslat )
	character id*56, pgm*8
	read( iunit, err=9 )  id, pgm, nc,nr,nz, xo,dx,yo,dy,
     1                     iproj, cmerid, baslat
	return
9	rewind iunit
	read( iunit, err=99 ) id, pgm, nc,nr,nz, xo,dx,yo,dy
	iproj  = 0
	cmerid = 999.
	baslat = 999.
	return
99    write(0,*) ' cannot read grid header'
	stop
	end
c******************************************************************************
	subroutine rowda( n, iz[huge], iop, irow, idev, ierr )
c  header is record 1, so irec = irow + 1
	dimension iz(n)
	ierr = 0
	irec = irow + 1
	if ( iop .eq. 0 ) then
	  write( idev, rec=irec, iostat=ierr ) ydum, iz
	  else
	  read ( idev, rec=irec, iostat=ierr ) ydum, iz
	endif
	return
	end
