c  interactive graphic editting of xyz data sets using local
c              gridding and contouring.
c  
c  mike webring, u.s. geological survey, 2/84, 6/84, 1/85
c
        common /fnames/gfile,dfile,rfile,bfile,sfile,fmtlab,isymb
	common /save/   icheck(50)
	common /exist/  idata,iexist
	common /dfile/  igrd,   iran,   iran2,  ingrd,  indat,
     1                  irf(14),ibf(14),isf(14),igf(14),idf(14)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,iflip,jflip
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncwrk(10),nrwrk(10)
	common /contou/ icoper,ncontd,nlevel,ncsam,nrsam,clevel(50),
     1                  ltype(50),xydist(400)
	common /edit/   fract,lswt1,iold,zold,znew,itrack,dummy(7)
	common /stalab/ sizsta,sizlab,angle,labfmt(3)
	common /boxcrn/ p1(2),p2(2)
c
	common /plot/   dxp(2),dyp(2),xp(4),yp(4),dxpw(2),dypw(2),
     1                  xpw(4),ypw(4),xpix,ypix,xscale,yscale
	common /plot2/  dypg(2),ypg(4),g2msca,g2myo
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	common /menu1/  nbox,menux,menudy
c
	common /datai/  xyz(1400)
	common /array/  wrk(20000)
c
	dimension za(20),iflag(10000)
        dimension plot1(26),menus(3)
        character*56 gfile,dfile,rfile,bfile,sfile
      character*1 isymb
	character*12 label(7),fmtlab
        character prompt*79
c      integer lc,ur,rc,lr
c
	equivalence (plot1,dxp(1)),(menus,nbox),
     1    (wrk(10001),iflag)
        data nwrk/20000/,
     1       nmenu/7/,label/'exit', 'contour', 'plt dat', 'expand',
     1       'contract', 'edit', 'setup'/
c initialize common variables
	call initc
	ngmax = nwrk / 2
c initialize graphics
c tspace = number of y pixels for 3 lines of text
c        = 40(.8 x 50.5(pixels/inch)) for Mitsubishi Diamond Scan
c        = 43(.8 x 53.84(pixels/inch)) for IBM Enhanced Color Monitor
c        = the above numbers are for EGA mode
    5 write(*,10)
c   10 format(' enter 0(Tandy 2000),1(EGA) or 2(VGA) : ',\)
   10 format(' enter 1(EGA) or 2(VGA) : ',\)
      read(*,*) imode
      if(imode.eq.0) then
        iplotr=7
        ial=8
        nhor=640
        tspace=46
        nvert=400-tspace
        icolor=500
      else if(imode.eq.1) then
        iplotr=9
        ial=2
        nhor=640
        tspace=40
        nvert=350-tspace
        icolor=300
      else if(imode.eq.2) then
        iplotr=10
        ial=2
        nhor=640
        tspace=55
        nvert=480-tspace
        icolor=300
      else
        print *,'must enter either a 0,1 or 2'
        go to 5
      endif
      ixo=600
      iyo=300
c
	write(*,20)
20	format(' grid filename : ',\)
	read(*,25) gfile
25	format(a56)
	call grd2da
	write(*,30)
30	format(' data filename : ',\)
	dfile=' '
	read(*,25) dfile
	if(dfile(1:5).ne.'     ') then
	  call opendf(nzrec,iz)
	  call makran(za,nzrec,iz)
	  idata=1
      write(*,35)
   35 format(' track sequential points',
     & ' in the edit function ? : ',\)
	  if(noyes(ldum).eq.1) then
	    call makseq(ierrsq)
	    itrack=1
	  endif
	endif
	call inits
	call chkcom(ierr)
	if (ierr.ne.0) go to 999
c
c  draw from igrd
   50 icoper=0
	icplot=0
	call cdvr(icplot,ncontd,ngmax,wrk,iflag)
	call boxall
	call whscal
	call menu(nmenu,label)
c
100	call chkcom(ierr)
	if (ierr.ne.0) then
          prompt=' static variable error'
	  call tpromp(0,0,0,prompt)
	  go to 999
	endif
        call fundvr(ifunc,ixo,iyo)
	go to (999,200,300,400,500,600,700)ifunc
c
c  contour data area
200	icplot=1
	call cdvr(icplot,ncontd,ngmax,wrk,iflag)
	go to 100
c
c  plot data 
300	if(idata.eq.0) then
          prompt=' no data to plot'
	  call tpromp(0,0,0,prompt)
	  go to 100
	endif
	call pltpt(0)
	go to 100
c
c  expand
400     call expand(iret,ixo,iyo)
	go to (50,100) iret
c
c  contract
500	call contrk(iret)
	go to (50,100) iret
c
c  edit function
600	if(idata.eq.0) then
          prompt=' no data to edit'
	  call tpromp(0,0,0,prompt)
	  go to 100
	endif
	if (icoper.ne.2) 
     1     call cdvr(1,ncontd,ngmax,wrk,iflag)
       call tdraw(1)
	call menu(nmenu,label)
       call tdraw(0)
        call edtctl(ial,ixo,iyo)
	call menu(nmenu,label)
	go to 100
c
c  setup function
700	call geset
	go to 100
c
999     if (idata.eq.1) then
	  call output(isave,za,nzrec,iz)
        if(isave.eq.0) then
        prompt='no output file has been created'
        call tpromp(0,0,0,prompt)
        endif
	endif
        prompt='do you want to exit the program : '
        call tpromp(0,2,0,prompt)
	if(noyes(ldum).eq.0) then
	  go to 100
	endif
	close(igrd)
	if (idata.eq.1)  close (indat)
	if (idata.eq.1)  close (iran)
        if (itrack.eq.1) close (iran2)
        call tclear(ial)
        m1=2
c        call tgin2(m1,m2,m3,m4)
	stop
	end
c******************************************************************************
        subroutine addpt(ixo,iyo)
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /edit/   dumm(5),itrack,dummy(7)
	common /stalab/ sizsta,sizlab,angle,labfmt(3)
	common /datai/  xyz(400),extra(1000)
	dimension dnew(3)
        character*79 prompt
        data lab1/1/
c
	npmax=nran/nval
10      prompt='move cursor to new data location : '
	call tpromp(0,0,0,prompt)
        call cursdu(dnew(1),dnew(2),ixo,iyo)
        prompt='enter z value : '
        call tpromp(0,0,0,prompt)
	read(*,*,end=10,err=10) dnew(3)
       write(prompt,15) dnew(3)
15     format('new z=',f11.4,' is this ok ? : ',\)
        call tpromp(0,0,0,prompt)
	if(noyes(ldum).eq.0) go to 10
	call lablpt(lab1,dnew)
c
c  add new point to outside area linked list
	nseq = nseq + 1
	irec = lstout
20	read(iran,rec=irec,err=999) next,ndp,xyz
	if (ndp.lt.npmax) then
c  add point to partially filled record
	  ix = ndp*nval + 1
	  xyz(ix)   = dnew(1)
	  xyz(ix+1) = dnew(2)
	  xyz(ix+2) = dnew(3)
	  xyz(ix+3) = float(nseq)
	  ndp = ndp + 1
	  call wrblk2(iran,irec,next,ndp,nran,xyz)
	else 
c  full record at end of list
c  change linking parameter of old maxrec
	  maxrec = maxrec + 1
	  next   = maxrec
	  call wrblk2(iran,irec,next,ndp,nran,xyz)
	  do 30 i = 1,nran
	    xyz(i) = 0.0
   30     continue
	  xyz(1)=dnew(1)
	  xyz(2)=dnew(2)
	  xyz(3)=dnew(3)
	  xyz(4) = float(nseq)
c  add new maxrec
	  irec = maxrec
	  next = 0
	  ndp  = 1
	  call wrblk2(iran,irec,next,ndp,nran,xyz)
	  lstout = maxrec
	endif
c
c  update sequence file
	if (itrack.eq.1) then
	  write(iran2,rec=nseq,err=80) dnew,irec,ndp
	endif
	go to 90
80	prompt=' error updating sequence file'
        call tpromp(0,0,0,prompt)
        prompt='>> another ? : '
        call tpromp(0,0,32,prompt)
        go to 95
c
90      prompt='>> another ?'
        call tpromp(0,0,0,prompt)
95      if(noyes(ldum).eq.1) go to 10
100	return
999	prompt=' error reading linked list'
        call tpromp(0,0,0,prompt)
      write(prompt,110) irec,maxrec
110   format('irec, maxrec',2i6)
        call tpromp(0,1,0,prompt)
c      call waiter(5.0)
	go to 100
	end
c******************************************************************************
	subroutine ammi(n,a,amn,amx,mn,mx,isave)
	dimension a(n)
	if(isave.gt.1) go to 1
	mn=1
	mx=1
	amn=a(1)
	amx=a(1)
1	do 3 i=1,n
	if(a(i).ge.amn) go to 2
	amn=a(i)
	mn=i
2	if(a(i).le.amx) go to 3
	amx=a(i)
	mx=i
3	continue
	return
	end
c******************************************************************************
	subroutine assign(slope,xogu,yogu,nc,nr,b[huge],np3,xyz)
c   associate np data values with grid locations of one block.
c   'slope' > 0 indicates that all data values within +-one half grid unit
c   will be combined by distance weighting.
c   'slope' = 0 indicates that the closest data value is the only one used.
	dimension xyz(np3),b(1)
	logical lslope
c
	if(np3.lt.3) return
	lslope=.false.
	if(slope.gt.0.0) lslope=.true.
c
	do 100 ixy=1,np3-2,3
	x=xyz(ixy)-xogu
	igx=int(x+.5)
	y=xyz(ixy+1)-yogu
	igy=int(y+.5)
c
	if(igx.gt.nc .or. igx.lt.1) go to 10
	if(igy.gt.nr .or. igy.lt.1) go to 20
	go to 30
10	write(*,11) igx
11	format(i5,'x')
	go to 100
20	write(*,21) igy
21	format(i5,'y')
	go to 100
c
30	z=xyz(ixy+2)
	ig=(igy-1)*nc+igx
	ib1=(ig-1)*6+1
	dx=x-float(igx)
	dy=y-float(igy)
	rsq=dx*dx+dy*dy
c
	if(lslope) go to 50
	if(rsq.ge.b(ib1)) go to 100
	b(ib1)=rsq
	b(ib1+1)=dx
	b(ib1+2)=dy
	b(ib1+5)=z
	go to 100
c
50	dwt=1.0/(rsq+slope)
	ib2=ib1+1
	ib3=ib1+2
	ib6=ib1+5
	b(ib1)=b(ib1)+dwt
	b(ib2)=b(ib2)+dx*dwt
	b(ib3)=b(ib3)+dy*dwt
	b(ib6)=b(ib6)+z*dwt
100	continue
	return
	end
c******************************************************************************
	subroutine bwts(izset,slope,binit,nc,nr,zg,iqd,b[huge])
c   initialize iqd and zg arrays.
c   calculate weights for minimum curvature equations.
c   'b' is dimensioned (6,nc*nr)
c   b(1,ig) contains either r squared or 1./weight
	dimension itabl(4)
	dimension zg(1),iqd(1),b(1)
	logical notslp
	data itabl/3,4,2,1/
c
	notslp=.true.
	if(slope.gt.0.0) notslp=.false.
	ib1=1
	do 100 ig=1,nc*nr
	if(b(ib1).eq.binit) go to 100
	ib2=ib1+1
	ib3=ib1+2
	ib6=ib1+5
	if(notslp) go to 50
	rwt=1./(b(ib1))
	b(ib2)=b(ib2)*rwt
	b(ib3)=b(ib3)*rwt
	b(ib6)=b(ib6)*rwt
	b(ib1)=b(ib2)*b(ib2)+b(ib3)*b(ib3)
c
50	if(b(ib1).gt.0.0025) go to 80
	iqd(ig)=-1
	zg(ig)=b(ib6)
	go to 100
80	ix=1
	if(b(ib2).ge.0.0) ix=2
	iy=0
	if(b(ib3).ge.0.0) iy=2
	iqd(ig)=itabl(ix+iy)
	if(izset.ne.0) zg(ig)=b(ib6)
100	ib1=ib1+6
c
c  calculate weighting
200	is=2*nc+3
	ie=3*nc-2
	do 230 j=3,nr-2
	ib1=(is-1)*6+1
	do 220 i=is,ie
	if(iqd(i))220,220,210
210	ib2=ib1+1
	ib3=ib2+1
	ib4=ib3+1
	ib5=ib4+1
	ib6=ib5+1
	dx=abs(b(ib2))
	dy=abs(b(ib3))
	dy2=dy*dy
	f1=dx*(dx+dy+dy+1.)
	b5=4./(f1+dy2+dy)
	b4=(b5*f1*.5)-1.
	b5dx=b5*dx
	b4b4=b4+b4
	b3=b5dx*(dy+1.)-b4b4
	b(ib2)=2.+b5dx-(b5*dy2+b4b4+b3)
	b(ib1)=b3+b4-b5dx
	b(ib3)=b3
	b(ib4)=b4
	b(ib5)=b5*b(ib6)
	b(ib6)=1./(1.+b(ib1)+b(ib2)+b3+b4+b5)
220	ib1=ib1+6
	is=is+nc
230	ie=ie+nc
	return
	end
c******************************************************************************
	subroutine border(nc,nr,zg,iqd,row)
	common /dfile/  igrd,dum(74)
	common /header/ dval, id(14), ipgm(2), nc1,dum1(8)
	common /wrkgrd/ iw,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncw(10),nrw(10)
	dimension isie(4)
	dimension zg(nc,nr),iqd(nc,nr),row(1)
        character*79 prompt
	data isie/1,2,0,0/
        if(nc.ne.ncw(iw)) then
        write(prompt,3) nc,ncw(iw)
    3   format('border : ',i5,' .ne. ',i5)
        call tpromp(0,1,0,prompt)
        endif
        if(nr.ne.nrw(iw)) then
        write(prompt,5) nr,nrw(iw)
    5   format('border : ',i5,' .ne. ',i5)
        call tpromp(0,2,0,prompt)
        endif
	isie(3) = nc-1
	isie(4) = nc
c  first 2 rows
	irec = irow(iw)
	do 20 j = 1, 2
	  call rowda (nc1,row,-1,irec,igrd,ierr)
	  ic = icol(iw)
	  do 10 i = 1, nc
	    zg(i,j) = row(ic)
	    if (zg(i,j).lt.dval) iqd(i,j) = -1
	    ic = ic+1
   10    continue
	irec = irec + 1
   20   continue
c  2 columns up the sides
	do 40 j = 3, nr-2
	  call rowda (nc1,row,-1,irec,igrd,ierr)
	  do 30 k = 1, 4
	    i = isie(k)
	    ic = icol(iw) + isie(k) - 1
	    zg(i,j) = row(ic)
	    if (zg(i,j).lt.dval) iqd(i,j) = -1
   30     continue
	  irec = irec+1
   40   continue
c  and the last 2 rows
	do 60 j = nr-1, nr
	  call rowda (nc1,row,-1,irec,igrd,ierr)
	  ic = icol(iw)
	  do 50 i = 1, nc
	    zg(i,j) = row(ic)
	    if (zg(i,j).lt.dval) iqd(i,j) = -1
	    ic = ic+1
   50     continue
	irec = irec + 1
   60   continue
	return
	end
c******************************************************************************
	subroutine box(p1,p2)
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	dimension p1(2),p2(2)
        m1=2
c        m2=0
c        m3=0
c        m4=0
        call tgin2(m1,m2,m3,m4)
	do 10 i = 1,2
	if(p1(i).gt.p2(i)) then
	  tmp=p1(i)
	  p1(i)=p2(i)
	  p2(i)=tmp
	endif
   10   continue
	call line(p1(1),p1(2),1,0,icolor)
	call line(p1(1),p2(2),1,1,icolor)
	call line(p2(1),p2(2),1,1,icolor)
	call line(p2(1),p1(2),1,1,icolor)
	call line(p1(1),p1(2),1,1,icolor)
c        m1=1
c        call tgin2(m1,m2,m3,m4)
	return
	end
c******************************************************************************
	subroutine boxall
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),if2(50)
	common /plot/   dxp(2),dyp(2),xp(4),yp(4),dum(16)
	character label*2
	if(icount.le.1) return
	size=.1
	do 20 ibox = 2,icount
	  call box(c1(1,ibox),c2(1,ibox))
	  call du2pu(c1(1,ibox),c1(2,ibox),xpu,ypu)
	  write(label,10) ibox
10	  format(i2)
	  call text(xpu+size,ypu+size,label,2,size,0.0)
   20   continue
	return
	end
c******************************************************************************
	subroutine cdvr(icplot,leveld,ngmax,z[huge],iflag[huge])
c  contour driver input: 
c  icplot, 0=init only, 1=contour only, 2=init&cont,
c  leveld, number of levels default, same as 'ncontd'
c  2 'ngmax' working arrays.
        common /fnames/gfile,dfile,rfile,bfile,sfile,fmtlab,isymb
	common /plot/   dxp(2),dyp(2),xp(4),yp(4),dxpw(2),dypw(2),
     1                  xpw(4),ypw(4),xpix,ypix,xscale,yscale
	common /contou/ icoper,ncontd,nlevel,nc,nr, c(50),ltype(50),
     1                  xy(400)
	common /stalab/   ifill2(3),labfmt(3)
	dimension zrange(2)
	dimension z(ngmax)
      character*1 isymb
        character*12 fmtlab
        character*56 gfile,dfile,rfile,bfile,sfile
        character fmtz*12,title*72
        character*79 prompt
        equivalence(fmtlab,fmtz)
	data maxlev/50/, dval/1.0e30/
	idval=0
	if(icplot.eq.1) go to 200
	icoper=0
c
c  assign levels default
	if (leveld.le.0 .and. ncontd.le.0) ncontd=20
	if (leveld.gt.0 .and. ncontd.le.0) ncontd=leveld
	levels=leveld
	if (levels.le.0) levels=ncontd
c
	call getz(ngmax,nc,nr,xy,z)
	nn = nc*nr
	if (nn.le.0) then
          prompt=' cdvr: getz returned nc*nr=0'
	  call tpromp(0,0,0,prompt)
	  return
	endif
c
c  get grid range and replace dval's
	call dvmm(nc*nr,z,zrange(1),zrange(2),dval,idval)
c
c  initialize plot, scale and axis tic
	dxp(1) = xy(1)
	dxp(2) = xy(nc)
	dyp(1) = xy(nc+1)
	dyp(2) = xy(nc+nr)
	call setgrf(1,dxp,dyp,xp,yp,0)
	xscale = (dxp(2)-dxp(1)) / xp(1)
	yscale = (dyp(2)-dyp(1)) / yp(1)
	call neatl
c
c  contour interval
	zint = 0.0
	max  = levels
	if (max.gt.maxlev) max = maxlev
	call setax(zrange,zint,max,nz,fmtz)
	nlevel = 1
	if (zint.le.0.0) zint = 1.0
	nlevel = 1 + int( (zrange(2)-zrange(1)) / zint )
	if (nlevel.lt.1)      nlevel = 1
	if (nlevel.gt.maxlev) nlevel = maxlev
c
c  assign contour levels and dashing patterns
c  line 1 is solid, 2-7 are dashed.
	idash  = 2
	iprime = 0
	pint = 5.0*zint
	cs = zint*aint(zrange(1)/zint + 0.5)
	do 10 i = 1, nlevel
	  c(i) = cs + float(i-1)*zint
	  t = abs( amod(c(i),pint) )
	  if (t.lt.1.0e-2*zint) then
	    if (iprime.eq.0) iprime = i
	    ltype(i) = idash
	    idash = idash + 1
	    if (idash.gt.7) idash = 2
	  else
	    ltype(i) = 1
	  endif
   10   continue
c
c  plot title
c	write(fmttit,100) fmtz,fmtz
c100	format(' (''min prime contour ='',', a,
c  1         ',''  interval ='',', a, ')')
c	write (title,fmttit) c(iprime), zint
       write(title,150) c(iprime),zint
150    format('min prime contour =',f11.4,
     1   '  interval =',f11.4)
	call text(.5,0.08,title,72,.08,0.0)
	icoper=1
	if (icplot.eq.0) return
c
c  contour z array
200	iy = nc+1
	call contur(nc,nr, xy,xy(iy), z,iflag, nlevel,c,ltype)
	icoper=2
	return
	end
c******************************************************************************
	subroutine changz(ichang,xyzs)
	common /stalab/ sizsta,sizlab,angle,labfmt(3)
	common /edit/   ifill(2),iold,zold,znew,itrack,dummy(7)
	dimension xyzs(4)
	character zchar*12
      character*1 icross
        character*79 prompt
	icross='4'
	call tdraw(0)
	call symbol(xyzs(1),xyzs(2),1,icross,2*sizsta,angle)
       write(prompt,2) xyzs(3)
2      format('> old z =',f11.4)
        call tpromp(0,0,0,prompt)
5	write(*,10)
10	format('> new z = ',\)
	zchar=' '
	znew=1.e-30
	read(*,20,end=5,err=5) zchar
20	format(a12)
	if(zchar.ne.'            ') read(zchar,30,end=40,err=5) znew
30	format(bn,f12.0)
40      call tdraw(1)
	call symbol(xyzs(1),xyzs(2),1,icross,2*sizsta,angle)
	if(znew.ne.1.e-30 .and. znew.ne.xyzs(3)) then
	  ichang=1
	  call lablpt(1,xyzs)
	  zold   = xyzs(3)
	  xyzs(3) = znew
	endif
	call tdraw(0)
	call lablpt(1,xyzs)
	if (ichang.eq.1 .and. itrack.eq.1) then
	  iseq = int (xyzs(4) + 0.01)
	  if (iseq.le.0) return
	  call trackz(ichang,iseq)
	endif
	return
	end		  
c******************************************************************************
	subroutine chkcom(ierr)
	common /save/   is(50)
	common /dfile/  igrd,iran,iran2,ingrd,indat,dum(70)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,dum1
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  x(4),iflip,jflip
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncwrk(10),nrwrk(10)
        character*79 prompt
	equivalence (dval,idval), (x,ix)
	dimension ix(4)
c
	ierr=0
	if (is(1) .ne. igrd)  go to 999
	if (is(2) .ne. iran)  go to 999
	if (is(3) .ne. ibsiz) go to 999
	if (is(4) .ne. nsec)  go to 999
	if (is(5) .ne. ntier) go to 999
	if (is(6) .ne. nran)  go to 999
	if (is(7) .ne. nval)  go to 999
	if (is(8) .ne. idval) go to 999
	if (is(9) .ne. nc)    go to 999
	if (is(10) .ne. nr)   go to 999
	if (is(11) .ne. ix(1)) go to 999
	if (is(12) .ne. ix(2)) go to 999
	if (is(13) .ne. ix(3)) go to 999
	if (is(14) .ne. ix(4)) go to 999
	if (is(15) .ne. iflip) go to 999
	if (is(16) .ne. jflip) go to 999
	if (is(17) .ne. ingrd) go to 999
	if (is(18) .ne. indat)  go to 999
c	is(19) is initial maxrec
c
	if(maxrec.gt.is(19)+2) then
          write(prompt,10) maxrec,is(19)
10        format('chkcom : maxrec, chkrec ',2i5)
          call tpromp(0,1,0,prompt)
	  ierr=1
	endif
	maxseq = maxrec*(nran/nval)
	if(nseq.lt.0 .or. nseq.gt.maxseq) then
          write(prompt,20) nseq,maxrec,nran,nval
20        format('chkcom : nseq,maxrec,nran,nval ',4i5)
          call tpromp(0,2,0,prompt)
	  go to 999
	endif
	return
999	ierr=1
	return
	end
c******************************************************************************
	subroutine contrk(iret)
	common /wrkgrd/ iwrkg,icount,dum(90)
	character prompt*79
1       write(prompt,10) icount,iwrkg
10      format('number of active grids is ',i2,
     &  ' ;present grid number = ',i2)
        call tpromp(0,0,0,prompt)
        write(prompt,20)
20      format('enter desired grid number(Gn);',
     &  '(Gn)<0 resets active grids to abs(Gn) : ',\)
        call tpromp(0,1,0,prompt)
        read(*,*,end=99,err=99) iw2
	if(iabs(iw2).gt.10 .or. iw2.gt.icount) go to 1
	if(iw2.eq.0 .or. iw2.eq.iwrkg) go to 88
	if(iw2.lt.0) then
	  iw2 = iabs(iw2)
	  if(iw2.gt.icount) go to 1
	  icount=iw2
          write(prompt,30) icount
30        format('number of active grids is : ',i2)
          call tpromp(0,1,0,prompt)
	endif
	iwrkg=iw2
	iret=1
	go to 99
c
88      write(prompt,40)
40      format('no change')
        call tpromp(0,1,0,prompt)
	iret=2
c      call waiter(1.0)
99	return
	end
c******************************************************************************
	subroutine contur(nx,ny, x,y, z[huge],un[huge], nc,cont,ltype)
c  specified dashing patterns for each contour require ltype(nc) storage
c  and ltype(1) > 0. 
	common /cfoll/  open,first,last,h,line
	dimension x(nx),y(ny),z(nx,ny),cont(nc),ltype(1)
	integer un(nx,ny)
	logical open,first,last
	jm=ny-1
	im=nx-1
	do 7 k=1,nc
	h=cont(k)
	if(ltype(1).le.0) then
	  line=1
	  else
	  line=ltype(k)
	endif
c  Evenden plotsys line range is 0:6
	line=line-1
	do 1 j=2,jm
	do 1 i=2,im
	un(i,j)=0
	if(z(i-1,j).lt.h .and. z(i,j).ge.h) un(i,j)=1
1	continue
	open=.true.
	do 2 i=2,nx
	if( z(i-1,1).lt.h .and. z(i,1).ge.h ) call follow(i,1,-1,0,
     1    x,y,z,nx,ny,un)
2	continue
	do 3 j=2,ny
	if( z(nx,j-1).lt.h .and. z(nx,j).ge.h ) call follow(nx,j,0,-1,
     1    x,y,z,nx,ny,un)
3	continue
	do 4 l=1,im
	i=nx-l
	if( z(i+1,ny).lt.h .and. z(i,ny).ge.h ) call follow(i,ny,1,0,
     1    x,y,z,nx,ny,un)
4	continue
	do 5 l=1,jm
	j=ny-l
	if( z(1,j+1).lt.h .and. z(1,j).ge.h ) call follow(1,j,0,1,
     1    x,y,z,nx,ny,un)
5	continue
	open=.false.
	do 6 l=2,jm
	j=ny-l+1
	do 6 m=2,im
	i=nx-m+1
	if( un(i,j).eq.1 ) call follow(i,j,-1,0,x,y,z,nx,ny,un)
6	continue
7	continue
	return
	end
c*****************************************************************************
        subroutine cursdu(xdu,ydu,ixo,iyo)
c  scales GIN location in data units
c  dxpw, dypw board area in data units; xpix, ypix data units/pixel
	common /plot/ dumm(12),dxpw(2),dypw(2),dumm2(8),xpix,ypix,
     1                dum(2)
      character*79 prompt
      ix=0
      iy=0
      prompt=' '
      call tpromp(0,1,0,prompt)
       call tgin(ix,iy,ixo,iyo)
	xdu=float(ix)*xpix+dxpw(1)
	ydu=float(iy)*ypix+dypw(1)
	return
	end
c******************************************************************************
	subroutine defbox(fract,p1,p2)
c  default box coordinates
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),if2(50)
	dimension p1(2),p2(2)
	xbox  = fract * (c2(1,iwrkg)-c1(1,iwrkg))
	ybox  = fract * (c2(2,iwrkg)-c1(2,iwrkg))
	p1(1) = 0.5 * (c2(1,iwrkg)+c1(1,iwrkg)) - 0.5*xbox
	p1(2) = 0.5 * (c2(2,iwrkg)+c1(2,iwrkg)) - 0.5*ybox
	p2(1) = xbox + p1(1)
	p2(2) = ybox + p1(2)
	return
	end
c******************************************************************************
	subroutine delpt(iswt,p1,p2,labswt)
c  delete points inside of box p1,p2 by negating the sequence number
c  iswt=-1 delete, iswt=0 change value, iswt=1 restore
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,dum1(2)
	common /actblk/ iblock,irec,next,np
	common /edit/   if2(2),ix,zold,znew,dum2(8)
	common /datai/  xyz(400),extra(1000)
	dimension p1(2),p2(2),ibx(2),iby(2)
	character prompt*79
c
	swt=float(iswt)
	if(iswt.lt.0) call tdraw(1)
	if(iswt.gt.0) call tdraw(0)
	call getblk(p1,p2,ibx,iby,iout)
	last=0
c
	do 50 jblk = iby(1), iby(2)
	if (jblk.eq.iby(2)) last=1
	do 50 iblk = ibx(1), ibx(2)+last
	  iblock = iblk + nsec*(jblk-1)
	  if (iblk .eq. ibx(2)+1) iblock = nsec*ntier+1
	  irec = iblock
10	read(iran,rec=irec,err=777) next,np,xyz
	ichang=0
	do 20 i=4,nval*np,nval
	  ix=i-3
	  if(xyz(ix).lt.p1(1) .or. xyz(ix).gt.p2(1)) go to 20
	  iy=ix+1
	  if(xyz(iy).lt.p1(2) .or. xyz(iy).gt.p2(2)) go to 20
	  if(iswt.eq.0) then
	    if (xyz(i).lt.0.0) go to 20
	    call changz(ichang,xyz(ix))
	  else
	    call lablpt(labswt,xyz(ix))
	    xyz(i)=sign(xyz(i),swt)
	    ichang=1
	  endif
20	continue
	if(ichang.eq.1) write(iran,rec=irec) next,np,xyz
	if(next.eq.0) go to 50
	irec=next
	go to 10
50	continue
	go to 999
c
777	write(prompt,888) iblk,jblk,irec
888	format(' delpt: block=',2i4,' attempt to read record',i15)
	call tpromp(0,0,0,prompt)
999	call tdraw(0)
	return
	end
c******************************************************************************
	subroutine du2pu(xdu,ydu,xpu,ypu)
	common /plot/  dxp(2),dyp(2),xp(4),yp(4),dxpw(2),dypw(2),
     1                 xpw(4),ypw(4),xpix,ypix,xscale,yscale
	xpu = (xdu-dxp(1))/xscale + xp(3)
	ypu = (ydu-dyp(1))/yscale + yp(3)
	return
	end
c******************************************************************************
	subroutine dvmm(n,a[huge],amn,amx,dval,idval)
c  find min&max grid values and whether dval's are present
	dimension a(n)
	idval=0
	amn=1.0e30
	amx=-amn
        itst=0
	do 10 i=1,n
	  if(a(i).ge.dval) then
	    idval=1
	  else
	    if(a(i).lt.amn) amn=a(i)
	    if(a(i).gt.amx) amx=a(i)
            itst=1
	  endif
   10   continue
	if (idval.eq.1) then
	  aset = amn - (amx-amn)/20.
	  do 20 i=1,n
	    if (a(i).ge.dval) a(i)=aset
   20     continue
          if(itst.eq.0) then
            amn=aset
            amx=aset
          endif
	endif
	return
	end
c******************************************************************************
        subroutine edtctl(ial,ixo,iyo)
c  edit level controller
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),if2(50)
	common /edit/   fract,labswt,iold,zold,znew,itrack,dummy(7)
	common /datai/  xyz(400),extra(1000)
	common /array/  wrk(20000)
	dimension p1(2),p2(2), iflag(10000)
	equivalence (wrk(10001),iflag)
	character*12 label(8)
        character*79 prompt
        data  ngmax/20000/
	data nmenu/8/, label/'return', 'size box', 'move  box', 'restore',
     1     'delete pt', 'add data', 'change z', 'grid&cont'/
c
        prompt=' '
        call tpromp(0,1,0,prompt)
        call tpromp(0,2,0,prompt)
        prompt='> label existing points ? :  '
        call tpromp(0,0,0,prompt)
	lab1=noyes(ldum)
	call pltpt(lab1)
	p1(1)=0.0
	p2(1)=0.0
	call defbox(fract,p1,p2)
	call box(p1,p2)
	call menu(nmenu,label)
c
10	call chkcom(ierr)
	if (ierr.ne.0) return
        call fundvr(ifunc,ixo,iyo)
	go to (999,200,300,400,500,600,700,800) ifunc
c
c  adjust box size
200	call tdraw(1)
	call box(p1,p2)
	call tdraw(0)
        call getbox(p1,p2,ixo,iyo)
	go to 10
c
300     call movbox(p1,p2,ixo,iyo)
	go to 10
c
c  restore points within box
400	call delpt(1,p1,p2,labswt)
	go to 10
c
c  delete points within box
500	call delpt(-1,p1,p2,labswt)
	go to 10
c
c  add point at cross-hair
600     call addpt(ixo,iyo)
	go to 10
c
c  change z value
700	call delpt(0,p1,p2,labswt)
	go to 10
c
800	continue
	call regrd(ierr)
	if(ierr.ne.0) go to 10
	call tclear(ial)
	call cdvr(2,0,ngmax,wrk,iflag)
	call pltpt(lab1)
	call boxall
	call box(p1,p2)
	call whscal
	call menu(nmenu,label)
	go to 10
c
999	call tdraw(1)
	call box(p1,p2)
	call menu(nmenu,label)
	call tdraw(0)
        prompt=' '
        call tpromp(0,1,0,prompt)
        call tpromp(0,2,0,prompt)
	return
	end
c******************************************************************************
        subroutine expand(iret,ixo,iyo)
c  wrkgrd common: igrd-  active unit number
c                 iwrkg- active grid index (1-10)
c                 icount-number of grids active (default=1)
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,dum(2)
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncw(10),nrw(10)
	common /boxcrn/ p1(2),p2(2)
        character*79 prompt
c
	iret=1
c  get expand area
      prompt= ' '
      call tpromp(0,0,0,prompt)
5	call tdraw(1)
        call box(p1,p2)
	call tdraw(0)
   10 prompt=' lower-left corner (''exit'' to cancel) : '
  	call tpromp(0,0,0,prompt)
        call menudvr(p1(1),p1(2),ifunc,ixo,iyo)
	if(ifunc.gt.0) then
       prompt=' function canceled'
	  call tpromp(0,0,0,prompt)
	  go to 999
	endif
      prompt= ' '
      call tpromp(0,0,0,prompt)
      prompt=' upper-right corner : '
  	call tpromp(0,0,0,prompt)
        call menudvr(p2(1),p2(2),ifunc,ixo,iyo)
	if(ifunc.gt.0) go to 10
        call box(p1,p2)
      prompt= ' '
      call tpromp(0,0,0,prompt)
      prompt=' ''exit'' if ok : '
	call tpromp(0,0,0,prompt)
        call menudvr(tx,ty,ifunc,ixo,iyo)
	if(ifunc.ne.1) go to 5
c
	if(icount.lt.1) icount=1
	icount=icount+1
	if(icount.gt.10) icount=10
	new = icount
c
	is = 1 + int((p1(1)-xo)/dx)
	js = 1 + int((p1(2)-yo)/dy)
	if (is.lt.1) is=1
	if (js.lt.1) js=1
c  extend across sampling interval 
	if (nskip(iwrkg).le.0) nskip(iwrkg)=1
	ie = nskip(iwrkg) + int((p2(1)-xo)/dx + 0.5)
	je = nskip(iwrkg) + int((p2(2)-yo)/dy + 0.5)
	if (ie.gt.nc) ie=nc
	if (je.gt.nr) je=nr
c  store new area
	c1(1,new) = xo + dx*float(is-1)
	c1(2,new) = yo + dy*float(js-1)
	c2(1,new) = xo + dx*float(ie-1)
	c2(2,new) = yo + dy*float(je-1)
	ncw(new)  = ie-is+1
	nrw(new)  = je-js+1
	icol(new) = is
	irow(new) = js
	iwrkg = new
	call tdraw(1)
	call box(p1,p2)
	call tdraw(0)
	return
c
c  undefined box
999	iret=2
	do 20 i=1,2
	  p1(i)=0.0
	  p2(i)=0.0
   20   continue
	return
	end
c******************************************************************************
	subroutine fit3c (dval,ihw,ip,jp,nc,nr,zg,iset)
c  Generate local 4 coefficient surface, width 2*ihw+1, and 
c  calculate value at zg(ip,jp).  Defined grid values must exist
c  in at least 2 diagonal quadrants.
	dimension a(9),b(3),c(3),iqt(4)
	dimension zg(nc,nr)
	tol=1.e-20
	iset=0
	idata=-1
	do 10 i=1,9
	  a(i)=0.0
	  if (i.le.3) b(i)=0.0
	  if (i.le.4) iqt(i)=0
   10   continue
	jqd = 0
	do 30 j = -ihw, ihw
	  jg = jp + j
	  if (j.gt.0) jqd = 2
	  iqd = 1
	  do 20 i = -ihw, ihw
	    if (i.gt.0) iqd = 2
	    ig = ip + i
	    if(zg(ig,jg).lt.dval) then
	      iqt(iqd+jqd) = 1
	      idata = idata + 1
	      if(idata.eq.0) then
	        zo = zg(ig,jg)
	        io = i
	        jo = j
	      else
	        z = zg(ig,jg) - zo
	        x = float(i-io)
	        y = float(j-jo)
	        x2 = x*x
	        y2 = y*y
	        a(1) = a(1) + x2
	        a(2) = a(2) + x*y
	        a(3) = a(3) + x2*y
	        a(5) = a(5) + y2
	        a(6) = a(6) + x*y2
	        a(9) = a(9) + x2*y2
	        b(1) = b(1) + x*z
	        b(2) = b(2) + y*z
	        b(3) = b(3) + x*y*z
	      endif
	    endif
   20     continue
   30   continue
	if (idata.lt.3) return
	x = -float(io)
	y = -float(jo)
	i1 = iqt(1) + iqt(4)
	i2 = iqt(2) + iqt(3)
c
	if (i1.ne.2 .and. i2.ne.2) return
	a(4) = a(2)
	a(7) = a(3)
	a(8) = a(6)
	a(6) = a(1)*a(6) - a(3)*a(4)
	a(9) = a(1)*a(9) - a(3)*a(7)
	b(3) = a(1)*b(3) - a(3)*b(1)
	a(5) = a(1)*a(5) - a(2)*a(4)
	a(8) = a(1)*a(8) - a(2)*a(7)
	b(2) = a(1)*b(2) - a(2)*b(1)
	a(9) = a(5)*a(9) - a(6)*a(8)
	b(3) = a(5)*b(3) - a(6)*b(2)
	if (abs(a(9)*a(5)).lt.tol) return
	c(3) =  b(3) / a(9)
	c(2) = (b(2) - a(8)*c(3)) / a(5)
	c(1) = (b(1) - a(7)*c(3) - a(4)*c(2)) / a(1)
	zg(ip,jp) = zo + x*c(1)+ y*c(2) + x*y*c(3)
	iset = 1
	return
	end
c******************************************************************************
	subroutine follow(if,jf,iaf,jaf,x,y,z[huge],nx,ny,un[huge])
c  follow one contour to either closure or grid boundary
	dimension x1(20),y1(20)
	dimension x(nx),y(ny),z(nx,ny)
	common /cfoll/ open,first,last,h,ltyp
      common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	integer un(nx,ny)
	logical open,last
	integer first,contin
c first=0 is true, =1 is false
	ixy=0
	first=0
	contin=0
	last=.false.
	i=if
	j=jf
	ia=iaf
	ja=jaf
	z1=z(i,j)
	za=z(i+ia,j+ja)
1	t=0.
	if(z1.ne.za) t=(z1-h)/(z1-za)
c
	ixy=ixy+1
	x1(ixy)=x(i)-t*(x(i)-x(i+ia))
	y1(ixy)=y(j)-t*(y(j)-y(j+ja))
	if(ixy.eq.20) then
	  call line(x1,y1,ixy,contin,ltyp+icolor)
	  contin=1
	  ixy=0
	endif
c
	if(open) go to 2
	if(ia.eq.-1 .and. un(i,j).eq.0) last=.true.
	go to 3
2	if(first.eq.0) go to 4
	if(ja.eq.0 .and. (j.eq.1 .or. j.eq.ny)) last=.true.
	if(ja.ne.0 .and. (i.eq.1 .or. i.eq.nx)) last=.true.
3	if(last) go to 4
	if(ia.eq.-1) un(i,j)=0
c
4	if(last .and. ixy.gt.0) call line(x1,y1,ixy,contin,ltyp+icolor)
c
	if(last) return
	zb=z(i+ja,j-ia)
	if(zb.ge.h) go to 5
	za=zb
	temp=ia
	ia=ja
	ja=-temp
	go to 7
5	zc=z(i+ia+ja,j-ia+ja)
	if(zc.ge.h) go to 6
	z1=zb
	za=zc
	i=i+ja
	j=j-ia
	go to 7
6	z1=zc
	i=i+ia+ja
	j=j-ia+ja
	temp=ja
	ja=ia
	ia=-temp
7	first=1
	go to 1
	end
c******************************************************************************
        subroutine fundvr(ifunc,ixo,iyo)
        character*79 prompt
      prompt=' '
      call tpromp(0,0,0,prompt)
        prompt=' enter menu function : '
 	call tpromp(0,0,0,prompt)
1       call menudvr(x,y,ifunc,ixo,iyo)
	  if(ifunc.eq.0) then
            prompt=' renter function : '
	    call tpromp(0,0,0,prompt)
c          call waiter(1.0)
	    go to 1
	  endif
	return
	end
c******************************************************************************
	subroutine geset
	common /exist/  idata,iexist
	common /contou/ icoper,ncontd,dum(503)
	common /edit/fract,labswt,iold,zold,znew,itrack,dummy(7)
	common /stalab/ sizsta,sizlab,angle,labfmt(3)
        character*79 prompt
c        modein=mode
c
	if(idata.eq.1) then
          write(prompt,1)
1         format('enable tracking in edit function ? : ',\)
          call tpromp(0,0,0,prompt)
	  itrack = noyes(ldum)
	  if (iexist.eq.0) call makseq(ierr)
	endif
c
        write(prompt,3)
3       format('label stations at edit level ? : ',\)
        call tpromp(0,0,0,prompt)
	labswt=noyes(ldum)
c
10      write(prompt,15)
15      format('enter box size (.01 to 1.0 of active grid) : ',\)
        call tpromp(0,0,0,prompt)
	read(*,*,end=10,err=10) fract
	if(fract.le.0.01) go to 10
	if(fract.gt.1.0)  go to 10
c
20      write(prompt,25)
25      format('enter number of contour intervals (1 to 50) : ',\)
        call tpromp(0,0,0,prompt)
	read(*,*,end=20,err=20) ncontd
c
30      write(prompt,35)
35      format('enter angle of data labels : ',\)
        call tpromp(0,0,0,prompt)
	read(*,*,end=30,err=30) angle
	return
	end
c******************************************************************************
	subroutine getblk(p1,p2,ibx,iby,iout)
c  get beginning and end block addresses
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,dum(2)
        character*79 prompt
	dimension p1(2),p2(2),ibx(2),iby(2)
c
	do 10  i=1,2
	  ibx(i)=0
	  iby(i)=0
   10   continue
	if(ibsiz.le.0 .or. maxrec.le.0) then
          write(prompt,20)
   20     format('getblk : error in data file addresses')
          call tpromp(0,1,0,prompt)
	  return
	endif
	rdb=1.0/(dx*float(ibsiz))
	x2 = 1.0 - rdb*xo
	y2 = 1.0 - rdb*yo
	ibx(1)=int(rdb*p1(1)+x2)
	ibx(2)=int(rdb*p2(1)+x2)
	iby(1)=int(rdb*p1(2)+y2)
	iby(2)=int(rdb*p2(2)+y2)
	iout=0
	if(ibx(1).lt.1 .or. ibx(2).gt.nsec)  iout=1
	if(iby(1).lt.1 .or. iby(2).gt.ntier) iout=1
	if(ibx(1).lt.1) ibx(1)=1
	if(ibx(2).gt.nsec)  ibx(2)=nsec
	if(iby(1).lt.1) iby(1)=1
	if(iby(2).gt.ntier) iby(2)=ntier
	return
	end
c******************************************************************************
        subroutine getbox(p1,p2,ixo,iyo)
	dimension p1(2),p2(2)
        character*79 prompt
	if(p1(1).eq.p2(1) .or. p1(2).eq.p2(2)) 
     1                                call defbox(0.10, p1,p2)
	call box(p1,p2)	
        prompt='lower-left box corner is reference'
        prompt(35:62)='; move cursor to adjust size'
	call tpromp(0,0,0,prompt)
	ipen = 0
        call cursdu(xdu,ydu,ixo,iyo)
	call tdraw(1)
	call box(p1,p2)
	call tdraw(0)
	p2(1) = xdu
	p2(2) = ydu
	call box(p1,p2)
	ipen = 1
	return
	end
c******************************************************************************
	subroutine getz(ngmax,nc2,nr2,xy,z[huge])
	common /dfile/   igrd,dum(74)
	common /header/  dval,id(14),ipgm(2),nc,nr,nz,
     1                   xo,dx,yo,dy,dum1(2)
	common /wrkgrd/  iwrkg,icount, c1(2,10),c2(2,10),nskip(10),
     1                   icol(10),irow(10),ncw(10),nrw(10)
        character*79 prompt
	dimension xy(1),z(1)
c
	ngrd = ncw(iwrkg)*nrw(iwrkg)
	nskip(iwrkg) = 1 + int( sqrt( float(ngrd)/float(ngmax) ))
	nc2 = ncw(iwrkg)/nskip(iwrkg)
	nr2 = nrw(iwrkg)/nskip(iwrkg)
c
	ngrd = nc2*nr2
	do 10 iz = ngrd+1, ngrd+nc2
	  z(iz) = dval	  
   10   continue
	iz = 1
	ioff = ngrd+1
	do 30  j=1,nr2
	  jrec = irow(iwrkg) + nskip(iwrkg)*(j-1)
	  call rowda (nc,z(ioff),-1,jrec,igrd,ierr)
          if(ierr.ne.0) then
          write(prompt,15) ierr
   15     format('getz : ierr = ',i5)
          call tpromp(0,1,0,prompt)
          endif
	  ix = icol(iwrkg)
	  do 20  i=1,nc2
	    z(iz) = z(ngrd+ix)
	    ix = ix + nskip(iwrkg)
	    iz = iz+1
   20   continue
   30     continue
c
	dx2 = dx*float(nskip(iwrkg))
	dy2 = dy*float(nskip(iwrkg))
	xo2 = xo + dx*float(icol(iwrkg) - 1)
	yo2 = yo + dy*float(irow(iwrkg) - 1)
	do 40  i=1,nc2
	  xy(i)= xo2 + float(i-1)*dx2
   40   continue
	ixy = nc2+1
	do 50  i=1,nr2
	  xy(ixy)= yo2 + float(i-1)*dy2
	  ixy = ixy+1
   50   continue
	return
	end
c******************************************************************************
	subroutine grd2da
c  convert input grid to random access, initialize wrkgrd common
        common /fnames/gfile,dfile,rfile,bfile,sfile,fmtlab,isymb
	common /dfile/  igrd,   iran,   iran2,  ingrd,  indat,
     1                  irf(14),ibf(14),isf(14),igf(14),idf(14)
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,iflip,jflip
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncw(10),nrw(10)
	common /array/  z(20000)
      character*1 isymb
        character*12 fmtlab
	character*56 gfile,dfile,rfile,bfile,sfile
        character*79 prompt
c
	open (ingrd, file=gfile, status='old', form='unformatted')
	read(ingrd) id,ipgm,nc,nr,nz,xo,dx,yo,dy
	if (dx.eq.0.0 .or. dy.eq.0.0) stop ' input dx or dy = zero'
	open (igrd, file=rfile, form='unformatted', status='unknown',
     1       access='direct', recl=nc*4)
c
c  shift origin to lower-left to avoid truncation problems
	iflip=0
	jflip=0
	if(dx.lt.0.0) then
	  dx=-dx
	  xo=-xo
	  iflip=1
          prompt='flipping x direction'
          call tpromp(0,1,0,prompt)
	endif
	if(dy.lt.0.0) then
	  dy=-dy
	  yo=-yo
	  jflip=1
          prompt='flipping y direction'
          call tpromp(0,1,0,prompt)
	endif
c
	iwrkg=1
	icount=1
	c1(1,1) = xo
	c1(2,1) = yo
	c2(1,1) = xo + dx*float(nc-1)
	c2(2,1) = yo + dy*float(nr-1)
	icol(1)=1
	irow(1)=1
	ncw(1)=nc
	nrw(1)=nr
c
	do 10  j=1,nr
	  call rowio(nc,z,-1,ingrd,ingrd,iend)
          if(iend.ne.0) then
          write(prompt,3) iend
    3     format('grd2da : iend = ',i5)
          call tpromp(0,1,0,prompt)
          endif
	  jrec=j
	  call rowda(nc,z,0,jrec,igrd,ierr)
          if(ierr.ne.0) then
          write(prompt,5) ierr
    5     format('grd2da : ierr = ',i5)
          call tpromp(0,2,0,prompt)
          endif
   10   continue
	close(ingrd)
	return
	end
c******************************************************************************
	subroutine grdctl(nc,nr,xo,yo,del,ierr)
c  minimum curvature gridding control for up to 2500 grid points.
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /header/ dval,dum1(25)
	common /datai/  xyz(400), extra(1000)
	common /array/  zg(2500),iqd(2500),b(15000)
	dimension c1(2),c2(2),ibx(2),iby(2)
	character prompt*79
	data binit/1.0/, slope/0.0/, bxo,byo/0.0,0.0/
	data izset/1/, nim/50/, epsm/0.0/
c
        prompt=' regridding...'
	call tpromp(0,0,0,prompt)
	ierr=0
	ndata=0
	nn=nc*nr
	if(nn.gt.2500 .or. nn.le.25) go to 91
	do 10  i=1,nn
	  iqd(i)=0
	  zg(i)=dval
   10   continue
	do 20  i=1,nn*6
	  b(i)=binit
   20   continue
c
	c1(1) = xo
	c1(2) = yo
	c2(1) = xo + del*float(nc-1)	
	c2(2) = yo + del*float(nr-1)
	call getblk(c1,c2,ibx,iby,iout)
	rdel = 1.0/del
	x2 = 1.0 - xo*rdel
	y2 = 1.0 - yo*rdel
c
c  read data blocks near the working grid
	last = 0
	do 50 jblk = iby(1), iby(2)
	  if (jblk.eq.iby(2)) last = 1
	  do 50 iblk = ibx(1), ibx(2)+last
	    irec = iblk + nsec*(jblk-1)
	    if(last.eq.1 .and. iblk.eq.ibx(2)+1) irec = nsec*ntier+1
c  
30	read (iran,rec=irec,err=99,iostat=ierr) next,np,xyz
	npout=0
c  eliminate data outside of working grid
	if (np.gt.0) call packxy(nval,np,c1,c2,xyz,npout)
	np3 = npout*3
	if (np3.ge.3) then
	  ndata = ndata + npout
c  convert to grid units
	  do 40  ix = 1, np3-2, 3
	    xyz(ix) = xyz(ix)*rdel + x2
	    iy = ix+1
	    xyz(iy) = xyz(iy)*rdel + y2
   40     continue
	  call assign(slope,bxo,byo,nc,nr,b,np3,xyz)
	endif
	if(next.eq.0) go to 50
	irec=next
	go to 30
50	continue
c
	if(ndata.eq.0) go to 99
	call bwts(izset,slope,binit,nc,nr,zg,iqd,b)
	call border(nc,nr,zg,iqd,xyz)
	call grdsrf(nc,nr,zg,xyz,dval,ierr)
	if(ierr.ne.0) go to 94
	call curvmn(zg,iqd,b,nc,nr,epsm,nim,st,end,ni)
	write (prompt,60) ndata,st,end
60	format(' ndata,start,end err',i6,2g15.4)
	call tpromp(0,0,0,prompt)
c      call waiter(2.0)
	return
91	ierr=1
	return
94	ierr=4
99	return
	end
c******************************************************************************
	subroutine grdsrf(nc,nr,zg,wz,dval,ier)
c  quickly replace dvals with reasonable values.
c  A control net spaced at 'nsep' is generated from local surfaces or
c  ring averages with radii 3, nsep; linear interpolation completes 
c  the process.  Array 'wz' is at least max(nc,nr).
        character*79 prompt
	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=2
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
	  call fit3c(dval,ihw,ii,jj,nc,nr,zg,ifit)
	  if(ifit.eq.1 .or. iflag.eq.2) 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(prompt,43) t
43	format(' gridr init with',1pe15.5)
        call tpromp(0,1,0,prompt)
	do 44 i=1,nn
	if(zg(i).ge.dval) zg(i)=t
44	continue
	return
	end
