c      interface to subroutine setcrt [c] (ial)
c      integer ial [reference]
c      end
        interface to subroutine setmod(a)
        integer*2 a
        end
	subroutine initc
        common /fnames/gfile,dfile,rfile,bfile,sfile,fmtlab,isymb
	common /exist/  idata,iseqex
	common /dfile/  igrd,   iran,   iran2,  ingrd,  indat,
     1                  irf(14),ibf(14),isf(14),igf(14),idf(14)
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncwrk(10),nrwrk(10)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,dum
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,iflip,jflip
	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 /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c
      character*1 isymb
        character*12 fmtlab
	character*56 gfile,dfile,rfile, bfile, sfile
	data dv/1.7e38/
c exist common
	idata = 0
	iseqex = 0
c dfile common
	igrd  = 10
	iran  = 11
	iran2 = 12
	ingrd = 13
	indat = 14
	rfile = 'grafrow.tmp'
	bfile = 'grafblk.tmp'
	sfile = 'grafseq.tmp'
	isymb='6'
c data common
	nran=400
	nval=4
c header common
	dval=dv
c contou common
	ncontd=20
c edit common
	fract=0.05
	lswt1=1
	itrack=0
c stalab common
	sizsta=0.03
	sizlab=0.08
	angle=0.0
	labfmt(1)=0
c boxall common
      icsw=1
	return
	end
c******************************************************************************
	subroutine inits
	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)
	equivalence (dval,idval), (x,ix)
	dimension ix(4)
c
	is(1) = igrd
	is(2) = iran
	is(3) = ibsiz
	is(4) = nsec
	is(5) = ntier
	is(6) = nran
	is(7) = nval
	is(8) = idval
	is(9) = nc
	is(10) = nr
	is(11) = ix(1)
	is(12) = ix(2)
	is(13) = ix(3)
	is(14) = ix(4)
	is(15) = iflip
	is(16) = jflip
	is(17) = ingrd
	is(18) = indat
	is(19) = maxrec
	is(20) = nseq
	return
	end
c******************************************************************************
	subroutine lablpt(labswt,xyz)
        common /fnames/gfile,dfile,rfile,bfile,sfile,fmtlab,isymb
	common /stalab/ sizsta,sizlab,angle,labfmt(3)
      character*1 isymb
        character*12 fmtlab
	character*56 gfile,dfile,rfile, bfile, sfile
	dimension xyz(3)
	character plab*16
c
	call symbol(xyz(1),xyz(2),1,isymb,sizsta,angle)
	if(labswt.eq.1) then
	  write(plab,fmtlab,err=99) xyz(3)
	  nch=leftj(plab,16)
	    do 10  ich=nch,2,-1
	    if(plab(ich:ich).ne.'0') then
	      call du2pu(xyz(1),xyz(2),xin,yin)
	      xin=xin+sizlab
	      call text(xin,yin,plab(1:ich),ich,sizlab,angle)
	      return
	    endif
   10       continue
	endif
99	return
	end	
c******************************************************************************
	function lenrec(idev)
c  find the number of words in a data file record
	integer*4 iz(24)
	do 15 j = 24,1,-1
	  read(idev,err=10) (iz(i),i=1,j)
	  go to 20
10	  rewind idev
15	continue
20	rewind idev
      lenrec=j
	return
	end
c******************************************************************************
	subroutine makran(za,nzrec,iz)
c  data are stored in nsec*ntier+1 pidgeon holes.
c  the last hole receives the data that fall outside the grid.
        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 /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /datai/  loca(200),ioff(200),extra(1000)
	common /array/  wrk(20000)
	dimension idsta(2),za(nzrec)
      character*1 isymb
        character*12 fmtlab
	character*56 gfile,dfile,rfile, bfile, sfile
c        character*79 prompt
	logical xyzrec
	data nwrk/20000/
c
	xyzrec = .true.
	if (nzrec.gt.1) xyzrec = .false.
	open (iran,   form='unformatted', status='scratch', 
     1        access='direct', recl=(nran+2)*4)
c
	del   = dx
	ibsiz = 10
10	nsec  = 1 + nc/ibsiz
	ntier = 1 + nr/ibsiz
	nblk  = nsec*ntier
	nblk1 = 1 + nblk
	if (nblk1.le.200) go to 20
	ibsiz = 1 + ibsiz
	go to 10
c
20	npb = int (float(nwrk)/float(nval*nblk1))
	if (npb.gt.nran/nval) npb = nran/nval
	nvb  = nval*npb
	nchk = nvb-nval
	rdb  = 1.0/(del*float(ibsiz))
	x2   = 1.0 - xo*rdb
	y2   = 1.0 - yo*rdb
c
c   'loca' contains the address where a block of data will be written.
c   a linked list is formed by 'next'.
	do 30  i=1,nblk1
	  ioff(i)=1
	  loca(i)=i
   30   continue
c   A data value resides in 4 contiguous words: x,y,z,and a sequence number
	do 40  i=1,nvb*nblk1
	  wrk(i)=0.0
   40   continue
c
c  read data, find pigeon hole
	next=nblk1+1
	ic=0
100	if (xyzrec) then
	  read(indat,end=200) xin,yin,z
	  else
	  read(indat,end=200) idsta,xin,yin,za
	  z = za(iz)
	endif
	if(iflip.ne.0) xin = -xin
	if(jflip.ne.0) yin = -yin
	ic = ic + 1
	if(mod(ic,1000).eq.0) then
        write(*,110) ic
  110   format(' count = ',i5)
        endif
	ibx  = int(rdb*xin+x2)
	iby  = int(rdb*yin+y2)
	mblk = (iby-1)*nsec+ibx
	if (ibx.lt.1 .or. ibx.gt.nsec)  mblk=nblk1
	if (iby.lt.1 .or. iby.gt.ntier) mblk=nblk1
c
c  put data in pigeon hole, output when full
	ixs=(mblk-1)*nvb
	ip=ixs+ioff(mblk)
	wrk(ip) = xin
	wrk(ip+1)=yin
	wrk(ip+2)=z
	wrk(ip+3)=float(ic)
	ioff(mblk)=ioff(mblk)+nval
	if(ioff(mblk).lt.nvb) go to 100
	ndp=(ioff(mblk)-1)/nval
	call wrblk2(iran,loca(mblk),next,ndp,nran,wrk(ixs+1))
	loca(mblk)=next
	next=next+1
	ioff(mblk)=1
	go to 100
c
c   output all pigeon holes
200	maxrec = next - 1
	ip=1
	next=0
	lstrec = 0
	do 50  i=1,nblk1
	  ndp=(ioff(i)-1)/nval
	  call wrblk2(iran,loca(i),next,ndp,nran,wrk(ip))
	  ip=ip+nvb
	  if (loca(i).gt.lstrec) lstrec = loca(i)
   50   continue
	if (maxrec.ne.lstrec) then 
        write(*,210) maxrec,lstrec
  210   format(' makran: next-1.ne.last record',2i5)
        endif
c
        write(*,220) ic
  220   format(i5,' data points')
	nseq = ic
c  the last outside (nblk1) record collects added points
	lstout = loca(nblk1)
	return
	end
c******************************************************************************
	subroutine makseq(ierr)
        common /fnames/gfile,dfile,rfile,bfile,sfile,fmtlab,isymb
	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 /datai/  xyz(400),extra(1000)
	dimension abc(3),i2(2)
      character*1 isymb
        character*12 fmtlab
	character*56 gfile,dfile,rfile, bfile, sfile
        character*79 prompt
        ierr=0
	if (iexist.eq.1) go to 1
c  create sequence file
	iexist = 1
	open (iran2,  status='scratch', recl=20,
     1        form='unformatted', access='direct')
	do 5  i=1,3
	  abc(i)=0.0
	  if(i.ne.3) i2(i)=0
    5   continue
	do 7  jrec = nseq+1, nseq+100
	 write (iran2,rec=jrec) abc,i2
    7   continue
c
c  update sequence file from blocked file
c  a deleted data point is flagged by a negative position pointer
1	iwrt = 0
	do 40 iblk = 1, nsec*ntier+1
	irec = iblk
10	read (iran,rec=irec,err=40) next,ndp,xyz
	ix = 1
	do 30 ipt = 1,ndp
	  ix3 = ix + 3
	  ipt2 = isign (ipt, int(xyz(ix3)) )
	  jrec = int (abs (xyz(ix3)) + .01)
	  write (iran2,rec=jrec,err=20) xyz(ix),xyz(ix+1),xyz(ix+2),
     1                                  irec,ipt2
	  ix = ix + nval
	  go to 30
20	  iwrt = iwrt + 1
30	continue
	irec = next
	if (irec.gt.0 .and. irec.le.maxrec) go to 10
40	continue
c
        if(iwrt.ne.0) then
        write(prompt,45) iwrt
45      format('makseg : ',i5,' write errors')
        call tpromp(0,1,0,prompt)
        endif
	return
	end
c******************************************************************************
	subroutine menu(nbox,label)
	common /plot/   dumm(12),dxpw(2),dypw(2),xpw(4),ypw(4),
     1                  xpix,ypix,dum(2)
	common /menu1/  mbox,menux,menudy
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	character*12 label(nbox)
c      m1=2
c      m2=0
c      m3=0
c      m4=0
c      call tgin2(m1,m2,m3,m4)
	mbox = nbox
	xsca = (dxpw(2)-dxpw(1)) / xpw(1)
	ysca = (dypw(2)-dypw(1)) / ypw(1)
	xleft = dxpw(2) - 1.0*xsca
	call line(xleft,  dypw(2),1,0,icolor)
	call line(dxpw(2),dypw(2),1,1,icolor)
	call line(dxpw(2),dypw(1),1,1,icolor)
	call line(xleft,  dypw(1),1,1,icolor)
	call line(xleft,  dypw(2),1,1,icolor)
	dy   = (dypw(2)-dypw(1)) / float(nbox)
	ybox = dypw(1) + dy
	do 10  i = 1, nbox-1
	  call line(xleft,  ybox,1,0,icolor)
	  call line(dxpw(2),ybox,1,1,icolor)
	  ybox = ybox + dy
   10   continue
	x = xleft + .1*xsca
	y = dypw(1) + dy/2.
	do 20  i = 1, nbox
	  call vchar(x,y,label(i),12,icsw*2,.1,0.,0.,0.)
	  y = y + dy
   20   continue
	menux  = (xleft-dxpw(1)) / xpix
	menudy = nvert / nbox
c        m1=1
c        call tgin2(m1,m2,m3,m4)
	return
	end
c******************************************************************************
        subroutine menudvr(x,y,ifunc,ixo,iyo)
	common /plot/ dumm(12),dxpw(2),dypw(2),dumm2(8),xpix,ypix,
     1                dum(2)
	common /menu1/ nbox,menux,menudy
	ifunc=0
	x=0.
	y=0.
        call tgin(ix,iy,ixo,iyo)
	if(ix.gt.menux) then
	  ifunc=1+iy/menudy
	  if(ifunc.gt.nbox) ifunc=nbox
	else
	  x=float(ix)*xpix+dxpw(1)
	  y=float(iy)*ypix+dypw(1)
	endif
	return
	end
c******************************************************************************
        subroutine movbox(p1,p2,ixo,iyo)
	dimension p1(2),p2(2)
        character*79 prompt
c  get new center
        prompt='new box center : '
	call tpromp(0,0,0,prompt)
        call cursdu(xcen,ycen,ixo,iyo)
c  undraw old box
	call tdraw(1)
	call box(p1,p2)
	call tdraw(0)
c  draw new box
	hx=0.5*(p2(1)-p1(1))
	hy=0.5*(p2(2)-p1(2))
	p1(1)=xcen-hx
	p2(1)=xcen+hx
	p1(2)=ycen-hy
	p2(2)=ycen+hy
	call box(p1,p2)
	return
	end
c******************************************************************************
	function noyes(idum)
	character inp*20
        character*79 prompt
        idum=0
	ic=0
	noyes=-1
1	continue
	read(*,2,err=5) inp
2	format(a20)
	if(ic.ge.3) go to 9
4	if(inp(1:1).eq.'y' .or. inp(1:1).eq.'Y') noyes=1
	if(inp(1:1).eq.'n' .or. inp(1:1).eq.'N') noyes=0
	if(noyes.gt.-1) return
5	write(*,6)
6	format(' y or n: ',\)
	ic=ic+1
	go to 1
9       write(prompt,10)
10      format('count exceeded, answering no')
        call tpromp(0,1,0,prompt)
	noyes=0
	return
	end
c******************************************************************************
	subroutine opendf (nz,iz)
c  open data file and find type
        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)
      character*1 isymb
        character*12 fmtlab
	character*56 gfile,dfile,rfile, bfile, sfile
	nz = 1
	iz = 1
	open (indat, file=dfile, form='unformatted', status='old')
	nword = lenrec(indat)
c  default record type  -- x,y,z --
	if (nword.gt.3) then
c  setup for record type -- id(2),x,y,z(nz) --
	  nz = nword - 4
	  if (nz.lt.1) stop ' 4 word data record not presently defined'
   10     write(*,11) nz
   11     format(' found ',i3,' data channels',
     &    ';which do you want ? : ',\)
	  read(*,*,end=10,err=10) iz
	  if (iz.gt.nz .or. iz.lt.1) go to 10
	endif
	return
	end
c******************************************************************************
	subroutine output(isave,za,nzrec,iz)
	common /dfile/  igrd,iran,iran2,ingrd,indat,dum(70)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /header/ ifill(23),dy,iflip,jflip
	common /datai/  xyz(400),extra(1000)
	dimension abc(3),za(nzrec)
        character prefix*56,ofile*56,ofile2*56
	character  idsta*8,suff*4,suff2*4
        character*79 prompt
	data iout,idel/15,16/, suff,suff2/'.edt', '.del'/
c
	prefix=' '
	ofile=' '
	ofile2=' '
	isave=0
	iwarn=0
	tol = 0.001*dy
	ic=0
	ic2=0
        prompt='do you want to output data : '
        call tpromp(0,0,0,prompt)
	if(noyes(ldum).eq.0) return
	write(prompt,5)
5	format('enter output filename prefix : ')
        call tpromp(0,0,0,prompt)
	read(*,37) prefix
37	format(a56)
c
c  skip over directory section, check for file prefix, strip extensions
	npre = leftj(prefix,56)
c      nbrac = 1 + index (prefix,']')
c      ndot = index (prefix(nbrac:npre),'.')
c      if(ndot.gt.0) then
c        do 10  i = ndot, npre
c        prefix(i:i)=' '
   10     continue
c        npre = ndot-1
c      endif
c
	ofile(1:npre) = prefix(1:npre)
        ofile(npre+1:npre+4) = suff
	open(iout,file=ofile,status='unknown',form='unformatted',err=88)
	ofile2=' '
	ofile2(1:npre) = prefix(1:npre)
        ofile2(npre+1:npre+4) = suff2
	open(idel,file=ofile2,status='unknown',err=99,iostat=ios)
c
c  update sequence file
	call makseq(ierr)
c
	if (nzrec.eq.1) then
c  output binary xyz file
	  do 30  iseq = 1, nseq
	    read(iran2,rec=iseq) abc,irec,ipt
	    if (iflip.ne.0) abc(1) = -abc(1)
	    if (jflip.ne.0) abc(2) = -abc(2)
	    if(ipt.gt.0) then
	      write(iout) abc
	      ic = ic + 1
	    else
	      write(idel,20) iseq, abc
20	      format(' #',i6,1x,1p3g13.5)
	      ic2 = ic2 + 1
	    endif
   30     continue
	  else
c  output a posting type record
	  rewind (indat)
	  do 60  iseq = 1, nseq
	    read(indat,end=50,iostat=istat) idsta,x,y,za
	    go to 57
50	    idsta=' '
	    do 55  i=1,nzrec
	      za(i)=0.0
   55       continue
57	    read(iran2,rec=iseq) abc,irec,ipt
	    if (iflip.ne.0) abc(1) = -abc(1)
	    if (jflip.ne.0) abc(2) = -abc(2)
	    za(iz) = abc(3)
	    xt = abs(x-abc(1))
	    yt = abs(y-abc(2))
	    if (istat.eq.0 .and. (xt.gt.tol .or. yt.gt.tol) ) then
	      iwarn = iwarn + 1
	      if(iwarn.eq.1) then
                write(prompt,51) x,y,abc(1),abc(2)
51              format('external x,y ',2f9.2,' internal x,y ',2f9.2)
                call tpromp(0,0,0,prompt)
                write(prompt,52)
52              format('lost registration of output file',
     &          '; continuing output')
                call tpromp(0,1,0,prompt)
	      endif
	    endif
	    if(ipt.gt.0) then
	      write(iout) idsta,abc(1),abc(2),za
	      ic = ic + 1
	    else
	      write(idel,53) iseq, idsta, abc(1),abc(2),za
53	      format(' #',i6,1x,a8,1p2g13.5,20g13.5)
	      ic2 = ic2 + 1
	    endif
   60     continue
	endif
c
100	close(iout)
	close(idel)
        write(prompt,110) ic,ic2
110     format('number saved, deleted = ',2i5)
        call tpromp(0,0,0,prompt)
        write(prompt,120) ofile(1:npre+4),ofile2(1:npre+4)
120     format('edit file is : ',a,' delete file is : ',a)
        call tpromp(0,1,0,prompt)
	isave=1
        if(iwarn.gt.0) then
        write(prompt,140) iwarn
140     format(i5,' possibly misregistered stations')
        call tpromp(0,2,0,prompt)
        endif
	return
88      write(prompt,150) ofile
150     format('open error on edit file : ',a)
        call tpromp(0,1,0,prompt)
	return
99      write(prompt,160) ios,ofile2
160     format('open error ',i5,' on delete file : ',a)
        call tpromp(0,1,0,prompt)
	close(iout)
	return
	end
c******************************************************************************
	function leftj(a,nch)
c  left justifies a string and returns the position
c  of the last nonblank character
	character a*(*)
	n = nch
	ich = ichar(a(1:1))
	if (ich.ne.0 .and. ich.ne.32) go to 15
	do 1 m = 2, n
	ich = ichar(a(m:m))
1	if (ich.ne.0 .and. ich.ne.32) go to 5
	leftj = 0
	return
c
5	i2 = 0
	do 10 i = m, n
	i2 = i2 + 1
	a(i2:i2) = a(i:i)
10	continue
	do 11 i3 = i2+1, n
11	a(i3:i3) = ' '
	n = n-m+1
c
15	do 20 leftj = n, 1, -1
	ich = ichar(a(leftj:leftj))
20	if (ich.ne.0 .and. ich.ne.32) go to 25
25	return
	end
c******************************************************************************
	subroutine packxy(nval,n,c1,c2,xyz,n2)
c  nval - number of values associated with each data point
c  n,n2 - are the number of input and output data points
	dimension c1(2),c2(2),xyz(400)
        character*79 prompt
	data jval/3/
	ival = nval
	nin  = n*ival
	jx   = 0
	do 10 iflag = ival, nin, ival
	if(xyz(iflag).gt.0.0) then
	  ix=iflag-3
	  if(xyz(ix).lt.c1(1) .or. xyz(ix).gt.c2(1)) go to 10
	  iy=ix+1
	  if(xyz(iy).lt.c1(2) .or. xyz(iy).gt.c2(2)) go to 10
	  do 5  k = ix, ix+2
	    jx = jx+1
	    xyz(jx) = xyz(k)
    5     continue
	endif
10	continue
	n2 = jx/jval
	itest = jx - n2*jval
        if(itest.ne.0) then
        write(prompt,20) jx,n2
20      format('packxy: index error jx,n2 : ',2i5)
        call tpromp(0,1,0,prompt)
        endif
	return
	end
c******************************************************************************
	subroutine pltpt(labswt)
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,dum1
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),if2(50) 
	common /datai/  xyz(400),extra(1000)
	dimension ibx(2),iby(2)
	character prompt*79
c
	lab1=labswt
	xmn=c1(1,iwrkg)
	xmx=c2(1,iwrkg)
	ymn=c1(2,iwrkg)
	ymx=c2(2,iwrkg)
	ic=0
	call getblk(c1(1,iwrkg),c2(1,iwrkg),ibx,iby,iout)
	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
10	read(iran,rec=irec,err=999) next,np,xyz
	do 20 i=4,nval*np,nval
	if(xyz(i).gt.0.0) then
	  ix=i-3
	  if(xyz(ix).lt.xmn .or. xyz(ix).gt.xmx) go to 20
	  iy=ix+1
	  if(xyz(iy).lt.ymn .or. xyz(iy).gt.ymx) go to 20
	  call lablpt(lab1,xyz(ix))
	  ic=ic+1
c        if(ic.ge.100) lab1=0
	endif
20	continue
	if(next.eq.0) go to 50
	irec=next
	go to 10
50	continue
c
200	write(prompt,210) ic
210	format(1x,i7,' data points')
	call tpromp(0,2,0,prompt)
	return
999	write(prompt,888) irec
888	format(' pltpt: attempt to read record',i15)
	call tpromp(0,1,0,prompt)
	return
	end
c******************************************************************************
	subroutine plugm(n,z,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 regrd(ierr)
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,dum1
	common /wrkgrd/ iw,icount,c1(2,10),c2(2,10),nskip(10),
     1                  icol(10),irow(10),ncw(10),nrw(10)
	common /header/ dval,id(14),ipgm(2),nc,nr,nz,
     1                  xo,dx,yo,dy,dum2(2)
	common /array/  zg(2500),iqd(2500),b(15000)
	common /datai/  row(1400)
	dimension zr(2)
	character prompt*79
	data ngmax/2500/
c
	ierr = 0
	nc2 = ncw(iw)
	nr2 = nrw(iw)
	nn = nc2*nr2
	if(nn.gt.ngmax .or. nn.le.25) then
	  write (prompt,1) nc2,nr2,ngmax
1	  format(' cannot regrid ',i3,' *',i3,' > ',i5,' grid points')
	  call tpromp(0,1,0,prompt)
	  ierr=1
	  return
	endif
	if(dx.ne.dy) then
          prompt=' cannot regrid unless dx=dy'
	  call tpromp(0,1,0,prompt)
	  ierr=1
	  return
	endif
c
c  grid the working area
	xo2 = xo + dx*float(icol(iw)-1)
	yo2 = yo + dy*float(irow(iw)-1)
	call grdctl(nc2,nr2,xo2,yo2,dx,ierr)
	if(ierr.ne.0) then
          write(prompt,3) ierr
3         format('regrd: grdctl error = ',i5) 
          call tpromp(0,1,0,prompt)
	  return
	endif
	call dvmm (nn,zg,zr(1),zr(2),1.e20,idval)
	if (idval.ne.0) ierr=1
	if (zr(1).lt.-1.e20) ierr=2
	if (ierr.ne.0) then
          write(prompt,5)
5         format('regrd: range or dval error')
          call tpromp(0,1,0,prompt)
	  return
	endif
c
c  update grid file
	ird = 0
	iwt = 0
	irec = irow(iw)
	do 30  j = 1,nr2
	  do 10  i = 1,nc
	    row(i)=dval
   10     continue
	  call rowda(nc,row,-1,irec,igrd,ierr)
	  if(ierr.ne.0) ird = ird + 1
	  ig = 1 + nc2*(j-1)
	  do 20  i = icol(iw), icol(iw)+nc2-1
	    row(i) = zg(ig)
	    ig = ig+1
   20     continue
	  call rowda(nc,row,0,irec,igrd,ierr)
	  if(ierr.ne.0) iwt = iwt + 1
	  irec = irec+1
   30   continue
c
        if(ird.ne.0) then
        write(prompt,40) ird
40      format('regrd : ',i5,' read errors on grid file')
        call tpromp(0,1,0,prompt)
        endif
        if(iwt.ne.0) then
        write(prompt,50) iwt
50      format('regrd : ',i5,' write errors on grid file')
        call tpromp(0,1,0,prompt)
        endif
	return
	end
c******************************************************************************
	subroutine rowda(n,z[huge],iop,irec,idev,ierr)
	dimension z(n)
	ierr=0
	if (iop.eq.0) then
	  write (idev,rec=irec,iostat=ierr) z
	  else
	  read  (idev,rec=irec,iostat=ierr) z
	endif
	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 setgrf(iscal,dxp,dyp,xp,yp,ihalfy)
c  input min,max data area coordinates, performs graph scaling call,
c  labels axis, and returns board info for pc screens.
c  iscal = 1 requires x,y scaling be the same
c  ihalfy = 0 use whole screen, ihalfy = -1,1 use lower or upper half
      common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
        dimension dxp(2),dyp(2),xp(4),yp(4)
	character fmtx*20,fmty*20
	dx=0.
	dy=0.
	call pltset(iplotr,xp(4),yp(4),1)
        yp(4)=yp(4)-.8
	call setax(dxp,dx,20,nx,fmtx)
	call setax(dyp,dy,20,ny,fmty)
	halfy = 0.5*yp(4) 
	ymarg = 0.0
	if (ihalfy.eq.1) ymarg = halfy
	xp(3) = .08*( float( max(nx,ny) ) + 1.0 )
	yp(3) = xp(3) + ymarg
	ywid  = yp(4)
	if (ihalfy.eq.-1) ywid = halfy
	xp(1) = xp(4) - xp(3) - 1.5
	yp(1) = ywid  - yp(3) - .05
	if (iscal.eq.1) then
	  xd = dxp(2) - dxp(1)
	  yd = dyp(2) - dyp(1)
	  xsc  = xd / xp(1)
	  ysc  = yd / yp(1)
	  scal = amax1(xsc,ysc)
	  if (abs(scal).lt.1.e-10) scal = 1.0e-10
	  xp(1) = xd / scal
	  yp(1) = yd / scal
	endif
	call scale(dxp,dyp,xp,yp,4,ier)
	if(iscal.eq.1) then
	  if(dy.lt.dx) then
	   dx   = dy
	   fmtx = fmty
	   nx   = ny
	   else
	   dy   = dx
	   fmty = fmtx
	   ny   = nx
	  endif
	endif
	call xaxis(dxp,dyp,xp,dx,5,.08,fmtx,nx)
	call yaxis(dyp,dxp,yp,dy,5,.08,fmty,ny)
        return
	end
c******************************************************************************
	subroutine setax(x,dx,maxint,nch,fmt)
c  adjust interval and labeling format
c  input x-min max range, optional dx-spacing and maxint-intervals
c  returns dx, nch-number of significant figures, fmt-labeling format
	dimension x(2)
	character fmt*9
	ixpn(r) = int( alog10(abs(r)) + 100. ) - 100
	fmt = '(1pe13.6)'
	nch = 13
	if (maxint.le.0)  maxint = 20
	dt = abs( x(2) - x(1) ) / float(maxint)
        if(dx.eq.0.0) then
          dx=1.0
          if(dt .lt. 1.0e-30) return
          p10 =sign(1.0,dt) * 10.**ixpn(dt)
	  t1  = dt / p10
	  if (t1.le.1.0) dx = p10
	  if (t1.gt.1.0) dx = 2.*p10
	  if (t1.gt.2.0) dx = 5.*p10
	  if (t1.gt.5.0) dx = 10.*p10
        endif
	m10 = ixpn( dx )
	m20 = ixpn( amax1(abs(x(1)),abs(x(2))) )
	idecm = 0
	if (m10.lt.0) idecm = iabs(m10)
	iw = 4 + idecm
	if (m20.ge.1) iw = iw + m20
	if (iw.gt.9 .or. idecm.gt.9) go to 99
	write (fmt,10) iw,idecm
10	format('(f',i1,'.',i1,')' )
	nch = iw
	if (idecm.eq.0) nch = nch - 1
99	return
	end
c******************************************************************************
	subroutine solv3(a,b,c,ierr)
	dimension a(9),b(3),c(3)
	ierr=1
	tol=1.e-20
	if(abs(a(1)).lt.tol) return
	if(abs(a(4)).lt.tol) return
	if(abs(a(3)).gt.tol) then
	  piv = a(1)/a(3)
	  aa = piv*a(6) - a(4)
	  bb = piv*a(9) - a(7)
	  cc = piv*b(3) - b(1)
	else
	  aa = a(6)
	  bb = a(9)
	  cc = b(3)
	endif
	if(abs(a(2)).gt.tol) then
	  piv = a(1)/a(2)
	  dd = piv*a(5) - a(4)
	  ee = piv*a(8) - a(7)
	  ff = piv*b(2) - b(1)
	else
	  dd = a(5)
	  ee = a(8)
	  ff = b(2)
	endif
	if(abs(aa).gt.tol) then
	  piv = dd/aa
	  gg = piv*bb - ee
	  hh = piv*cc - ff
	else
	  gg = bb
	  hh = cc
	endif
	if (abs(gg).lt.tol .or. abs(dd).lt.tol) return
	c(3) = hh/gg
	c(2) = (ff - ee*c(3)) / dd
	c(1) = (b(1) - a(7)*c(3) - a(4)*c(2)) / a(1)
	ierr=0
	return
	end
c******************************************************************************
	subroutine text(x,y,string,nstr,size,angle)
c  text driver where x,y are in plot units.
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c        dimension is(25)
	character s*100
	character string*(*)
	data icode/3/, xoff,yoff/0.0,0.0/
	lenstr=nstr
	do 5  nch=lenstr,1,-1
	  ich=ichar(string(nch:nch))
	  if(ich.gt.32 .and. ich.lt.127) go to 10
    5   continue
        go to 20
10	s(1:nch)=string(1:nch)
	radian=1.7453292e-2*angle
	call vchar(x,y,s,nch,icsw*icode,size,radian,xoff,yoff)
   20   return
	end
c******************************************************************************
	subroutine symbol(x,y,nxy,ich,size,angle)
c  symbol driver where x,y are in data units (clipping occurs outside
c  the data window).
c  'ich' is an integer akin to fortran 77 coding ich=038=ichar('&').
c  'vchar'  uses 1-13, 24-31 for symbols such as diamonds and squares, 
c  and 32-126 as the normal ascii character sequence when icode=1.
      character ich*(*)
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	dimension x(nxy),y(nxy)
	data icode/1/, xoff,yoff/0.0,0.0/
	radian=1.7453292e-2*angle
	call vchar(x,y,ich,nxy,icsw*icode,size,radian,xoff,yoff)
	return
	end
c******************************************************************************
	subroutine trackz(ichang,iseq)
c  Track sequential data points with zold and update to znew.
c  note xyz is different from the active block xyz and this routine updates
c  point iseq.  ipyram(1) is a pyramid pointing down, (2) points up.
	common /dfile/  igrd,iran,iran2,dum(72)
	common /data/   ibsiz,nsec,ntier,nran,maxrec,nval,nseq,lstout
	common /wrkgrd/ iwk,icount,c1(2,10),c2(2,10),dum1(50)
	common /edit/   fract,labswt,iold,zold,znew,dummy(8)
	common /stalab/ sizsta,sizlab,angle,ifmt(3)
	common /actblk/ iblock,irecin
	common /array/  xyz(400),dum2(19600)
	dimension nt(2)
	dimension xt(20),yt(20),zt(20),irt(20),ipt(20)
      character*2 ipyram(2)
	character prompt*79
	data maxtrk/11/, ipyram/'12','11'/
c
	size2  = sizlab
	xmin = c1(1,iwk)
	xmax = c2(1,iwk)
	ymin = c1(2,iwk)
	ymax = c2(2,iwk)
	if (iseq.le.0 .or. iseq.gt.nseq) return
	it=1
	read (iran2,rec=iseq,err=999) xt(1),yt(1),zt(1),irt(1),ipt(1)
	call symbol(xt(it),yt(it),1,ipyram(1),size2,angle)
c
c  track negative then positive direction
	do 100 idir = 1,2
	ione = -1
	maxt = 1 + maxtrk/2
	if (idir.eq.2) then
	  ione = 1
	  maxt = maxtrk
	endif
c  step through the sequential file
	jrec = iseq + ione
10	if (jrec.le.0 .or. jrec.gt.nseq) go to 100
	it = it + 1
	read (iran2,rec=jrec,err=50) 
     1        xt(it),yt(it),zt(it),irt(it),ipt(it)
	if (zt(it).eq.zold) then
	  if (xt(it).lt.xmin .or. xt(it).gt.xmax) go to 50
	  if (yt(it).lt.ymin .or. yt(it).gt.ymax) go to 50
	  call symbol(xt(it),yt(it),1,ipyram(idir),size2,angle)
 	  jrec = jrec + ione
	  if (it.lt.maxt) go to 10
	endif
	go to 100
50	it = it - 1
100	nt(idir) = it
	nt(2) = nt(2) - nt(1)
c
        write(prompt,105) zold,nt(1),nt(2)
  105   format('tracking',f9.3,' found',i2,' on down side',
     1  ' and',i2,' on up side')
        call tpromp(0,0,0,prompt)
  110   write(prompt,106)
  106   format(' enter number to change on down side : ')
        call tpromp(0,1,0,prompt)
	read(*,*,err=110) nt1
	if (nt1.eq.0) nt1 = 1
        write(prompt,107)
  107   format(' enter number to change on up side : ')
        call tpromp(0,1,0,prompt)
	read(*,*,err=110) nt2
        prompt=' '
        call tpromp(0,1,0,prompt)
	if (nt1.lt.0 .or. nt1.gt.nt(1)) nt1 = nt(1)
	if (nt2.lt.0 .or. nt2.gt.nt(2)) nt2 = nt(2)
	if(nt1.lt.nt(1) .or. nt2.lt.nt(2)) then
	  call tdraw(1)
	  do 120  i = nt(1), nt1+1, -1
	    call symbol(xt(i),yt(i),1,ipyram(1),size2,angle)
  120     continue
	  do 130  i = nt(1)+nt(2), nt(1)+nt2+1, -1
	    call symbol(xt(i),yt(i),1,ipyram(2),size2,angle)
  130     continue
	  call tdraw(0)
	endif
c
c  update block and sequence files
c  "deleted" points are updated also
	ierrt = 0
	do 150  idir = 1, 2
	is = 1
	ie = nt1
	if (idir.eq.2) then
	  is = nt(1)+1
	  ie = nt(1)+nt2
	endif
c
	do 140  i = is, ie
	  read(iran,rec=irt(i),err=200,iostat=ierr) next,ndp,xyz
	  ix = 1 + nval*(ipt(i)-1)
	  if (xyz(ix).eq.xt(i)) then
	    call tdraw(1)
	    call lablpt(1,xyz(ix))
	    xyz(ix+2) = znew
	    call tdraw(0)
	    call lablpt(1,xyz(ix))
	    jrec = int (abs (xyz(ix+3)) + .01)
	    write (iran, rec=irt(i),err=200,iostat=ierr) next,ndp,xyz
	    if (i.eq.1) ichang=0
	    write (iran2,rec=jrec,  err=200,iostat=ierr) 
     1      xt(i),yt(i),znew,irt(i),ipt(i)
	  else
            write(prompt,135)
  135       format('trackz : xyz address error')
            call tpromp(0,1,0,prompt)
	  endif
200	  if (ierr.ne.0) ierrt = ierrt + 1
  140   continue
  150   continue
c
	if (ierrt.ne.0) then
	  write(prompt,400) ierrt
400	  format(' trackz:',i5,' random access read/write errors')
	  call tpromp(0,1,0,prompt)
	endif
999	return
	end
c******************************************************************************
	subroutine tdraw(iop)
c  drawing mode:  iop= 0 pixel on, 1 off, 2 toggled
        common/boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
        iop1=iop+1
        go to (10,20,30),iop1
   10   icolor=0
        icsw=1
        go to 40
   20   icolor=-7
        icsw=-1
        go to 40
   30   if(icolor.eq.0) then
        icolor=-7
        icsw=-1
        else
        icolor=0
        icsw=1
        endif
   40   return
	end
c******************************************************************************
	subroutine tclear(iop1)
c  clear screen
        integer*2 jop1
        jop1=iop1
      call setmod(jop1)
	return
	end
c******************************************************************************
        subroutine tgin(ix,iy,ixo,iyo)
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c  gets cursor position
        m1=0
        call mouse(m1,m2,m3,m4)
        m2=0
        m1=4
        call mouse(m1,m2,ixo,iyo)
        m1=1
        call mouse(m1,m2,m3,m4)
        m2=0
        m1=3
c  loop on n added to fix Windows95 bug
   11   n=0
   10   call mouse(m1,m2,ix,iy)
        n=n+1
        if(n.lt.10000) go to 10
        if(m2.eq.0) go to 10
        if(n.eq.10000) go to 11
        iyo=iy
        ixo=ix
        iy=nvert+tspace-iy
        return
      end
c******************************************************************************
      subroutine tgin2(m1,m2,m3,m4)
      call mouse(m1,m2,m3,m4)
      return
      end
c******************************************************************************
	subroutine tpromp(iwrt,line,icol,prompt)
        character*16 ifmt
        character bpromp*70
	   character prompt*79
        data bpromp/' '/
        iwrt=0
        call setcur(line,icol)
        write(*,*) bpromp
        call setcur(line-1,icol)
        if(icol.eq.0) icol=1
        num=lench(prompt)
c        write(ifmt,10) icol,inum
        write(ifmt,10) icol,num
   10   format('(',i2,'x,a',i2,'\)')
        write(*,ifmt) prompt(1:num)
	return
	end
c***********************************************************************
        function lench(data)
        character data*79
        do 10 i=79,1,-1
        if(data(i:i).ne.' ') go to 20
   10   continue
   20   lench=i
        if(lench.lt.1) lench=1
        return
        end
c******************************************************************************
	subroutine whscal
c  rescale entire board in coordinate system of dxp,dyp
c  xscale, yscale are data units/inch
	common /plot/ dxp(2),dyp(2),xp(4),yp(4),dxpw(2),dypw(2),
     1                xpw(4),ypw(4),xpix,ypix,xscale,yscale
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	dimension a(16)
        character*79 prompt
	equivalence (a,dxpw(1))
	do 1 i=1,16
1	a(i)=0.
	xpw(1)=xp(4)
	ypw(1)=yp(4)
	xpw(4)=xp(4)
	ypw(4)=yp(4)
	xscale=(dxp(2)-dxp(1))/xp(1)
	yscale=(dyp(2)-dyp(1))/yp(1)
	xpix=xscale*xp(4)/nhor
	ypix=yscale*yp(4)/nvert
	dxpw(1)=dxp(1)-xp(3)*xscale
	dypw(1)=dyp(1)-yp(3)*yscale
	dxpw(2)=xscale*xp(4)+dxpw(1)
	dypw(2)=yscale*yp(4)+dypw(1)
	call scale(dxpw,dypw,xpw,ypw,3,ier)
        if(ier.ne.0) then
        write(prompt,10)
10      format('whscal : scaling error')
        call tpromp(0,1,0,prompt)
        endif
	return
	end
c******************************************************************************
	subroutine wrblk2(iunit,iorec,next,ndp,nxyz,xyz)
	dimension xyz(nxyz)
        character*79 prompt
	write(iunit,rec=iorec,err=99,iostat=ierr) next,ndp,xyz
	return
99      write(prompt,100) iorec,ierr
100     format('wrblk2 : write error; record = ',i5,
     &  ' io status = ',i5)
        call tpromp(0,1,0,prompt)
	return
	end
c******************************************************************************
	subroutine curvmn(zg,iqd,b[huge],nc,nr,epsmx,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.
c  Mike Webring, USGS Open-File report 81-1224.
	dimension zg(2500),iqd(2500),b(15000)
	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
	


d

