        common /fnames/gfile,rfile
        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 /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,ncsam,nrsam,clevel(50),
     1               ltype(50),xydist(400)
	common /mparms/ efld,einc,edec,tfdir(3)
	common /dparms/ ndike,dstore(7,20),drange(2,20)
        common /pparms/ iprf,nprf,nseg(10),sx(10,10),sy(10,10)
c
	common /plot2/  dypg(2),ypg(4),g2msca,g2myo
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
	common /tty1/   iterm,mode,line2
	common /tty/    ipen,ich
	common /menu1/  nbox,menux,menudy
	common /array/  wrk(40000)
c
	dimension    iflag(20000)
	dimension    xlmt(2),ylmt(2)
	dimension t(500),h(500),rx(500),ry(500),rz(500)
        character*12 label(7)
        character    prog*8
        character*56 gfile,rfile
	equivalence  (wrk(20001),iflag)
c        equivalence  (igf,gfile), (plot1,dxp(1)),(menus,nbox)
        equivalence (wrk,t), (wrk(501),h), (wrk(1001),rx),
     1              (wrk(1501),ry), (wrk(2001),rz)
        equivalence  (plot1,dxp(1)),(menus,nbox)
c
	data ity/5/, nwrk/40000/,
     1 nmenu/7/, label/'exit',   'contour', 'profile' , 'delete pt',
     1                 'expand', 'contract', 'save grid'/
        data iwrk/3001/
        data r2d/57.29578/
c
	prog = 'profilex'
cc        call user(prog)
        call initpx
	ngmax = nwrk / 2
cc        call tmode(1)
cc        call tclear(0)
c
c       print 11
c11	format(' data: grid = 0, synthetic = 1 : '$)
    5 write(*,10)
   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
c
c       read(ity,*) isyn
c        isyn = 0
c        if(isyn.eq.0) then
          write(*,20)
20        format(' grid filename: ',\)
	  read(ity,37) gfile
37	  format(a)
	  call grd2da
c        endif
c
c  draw from igrd
c
cc50      call tclear(0)
cc        call tmode(3)
      ixo=600
      iyo=300
50      icoper=0
	icplot=0
	call cdvr(icplot,ncontd,ngmax,wrk,iflag)
	call boxall
	call pltseg
	call whscal
	call menu(nmenu,label)
c
c  function switching
c
100     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  profile extraction
c300     if (isyn.eq.0) then
300       xlmt(1) = c1(1,1) + dx
	  ylmt(1) = c1(2,1) + dy
	  xlmt(2) = c2(1,1) - dx
	  ylmt(2) = c2(2,1) - dy
          call getseg(ipen,xlmt,ylmt,ifunc2,ixo,iyo)
c        endif
        if(ifunc2.ne.1) go to 50
        dr = 0.0
        call getprf(igrd,nc,nr,xo,dx,yo,dy, nseg(iprf),sx(1,iprf),
     1              sy(1,iprf), dr,nt,rx,ry,t, wrk(iwrk) )
        call outprf (nseg(iprf),sx(1,iprf),sy(1,iprf),
     1               dr,nt,rx,ry,rz,t, wrk(iwrk) )
cc        call prfctl(isyn,ifunc2)
	go to 50
c
c  delete point on the profile handled in getseg
400	go to 100
c
c  expand view
500     call expand(ipen,iret,ixo,iyo)
	go to (50,100) iret
c
c  contract view
600	call contrk(iret)
	go to (50,100) iret
c
c  output active grid
700	call outgrd
	go to 100
c
cc999     call tclear(0)
cc        call tmode(1)
999     close(igrd)
        call tclear(ial)
        m1=2
        call tgin2(m1,m2,m3,m4)
        print*
        if(nprf.gt.0) then
          do 998 iprf=1,nprf
          az=r2d * atan2(sx(nseg(iprf),iprf)-sx(1,iprf),
     1       sy(nseg(iprf),iprf)-sy(1,iprf))
          print*,'azimuth of profile ',iprf,' = ',az
 998      continue
        endif
	stop
	end
c****************************************
        subroutine initpx
	common /dfile/  igrd,   iter,   iobs,   ingrd,  indat,
     1               ihf(14),itf(14),iof(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 /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,ncsam,nrsam,clevel(50),
     1               ltype(50),xydist(400)
	common /mparms/ efld,einc,edec,tfdir(3)
	common /dparms/ ndike,dstore(7,20),drange(2,20)
	common /plot2/  dypg(2),ypg(4),g2msca,g2myo
	common /tty1/   iterm,mode,line2
	common /tty/    ipen,ich
	common /menu1/  nbox,menux,menudy
        common /boxcol/icolor,icsw,iplotr,nvert,nhor,tspace
c
	character hfile*56
	equivalence  (ihf,hfile)
c        data dv/'ffff7fff'x/
        data dv/1.e38/
c
c  dfile common
	igrd  = 10
	iter  = 11
	iobs  = 12
	ingrd = 13
	indat = 14
	hfile = 'pxmag.tmp'
c
c  mparms common 
	efld = 50000.0
	einc = 90.0
	edec = 0.0
c
	dval = dv
	ipen = 1
	ncontd = 20
c boxcol common
        icsw=1
	return
	end
c****************************************
        subroutine grd2da
c  convert input grid to random access, initialize wrkgrd common
        common /fnames/gfile,rfile
	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(40000)
        character*56 gfile,rfile
        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'
c        open (igrd, file=rfile, form='unformatted', status='unknown',
        open (igrd, 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 '
cc          call tpromp(0,1,0,prompt)
          call tpromp(0,0,0,prompt)
	endif
	if(dy.lt.0.0) then
	  dy=-dy
	  yo=-yo
	  jflip=1
          prompt='flipping y direction '
c          call tpromp(0,1,0,prompt)
          call tpromp(0,0,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,' ')
c          call tpromp(0,1,0,prompt)
          call tpromp(0,0,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,' ')
c          call tpromp(0,2,0,prompt)
          call tpromp(0,0,0,prompt)
          endif
   10   continue
	close(ingrd)
	return
        end

c****************************************
        subroutine rowio(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
c****************************************
        subroutine rowda(n,z,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 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 /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(4),labfmt(3)
	dimension zrange(2)
	dimension z(ngmax),iflag(ngmax)
	character id*56,p*8,fmtz*12,fmttit*72,title*72
	equivalence (labfmt,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
c          call tpromp(1,1,1,' cdvr: getz returned nc*nr=0')
          call tpromp(0,0,0,' cdvr: getz returned nc*nr=0')
	  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 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
	enddo
c
c  plot title
	write(fmttit,100) fmtz,fmtz
100	format(' (''min prime contour ='',', a,
     1      ',''  interval ='',', a, ')')
	write (title,fmttit) c(iprime), zint
	call text(.5,0.08,title,.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 getz(ngmax,nc2,nr2,xy,z[huge])
c  prepare contourable z array
        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)
	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 iz = ngrd+1, ngrd+nc2
	  z(iz) = dval	  
	enddo
	iz = 1
	ioff = ngrd+1
	do j=1,nr2
	  jrec = irow(iwrkg) + nskip(iwrkg)*(j-1)
	  call rowda (nc,z(ioff),-1,jrec,igrd,ierr)
          if (ierr.ne.0) print *,' getz: ierr=', ierr
          ix = icol(iwrkg)
	  do i=1,nc2
	    z(iz) = z(ngrd+ix)
	    ix = ix + nskip(iwrkg)
	    iz = iz+1
	  enddo		
	enddo
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 i=1,nc2
	  xy(i)= xo2 + float(i-1)*dx2
	enddo
	ixy = nc2+1
	do i=1,nr2
	  xy(ixy)= yo2 + float(i-1)*dy2
	  ixy = ixy+1
	enddo
	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
c        print'(a)',prompt
        num=lench(prompt)+1
c        write(ifmt,10) icol,inum
        write(ifmt,10) icol,num
   10   format('(',i2,'x,a',i2,'\)')
c        print*,icol,num
        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
        return
        end
c****************************************
        subroutine dvmm(n,a,amn,amx,dval,idval)
        dimension a(n)
        amn=dval
        amx=-dval
1       do 4 i=1,n
        if(a(i).ge.dval) go to 3
        if(a(i).ge.amn) go to 2
        amn=a(i)
2       if(a(i).le.amx) go to 3
        amx=a(i)
        go to 4
3       idval=1
4       continue
        return
        end
c****************************************
c        subroutine dvmm(n,a[huge],amn,amx,dval,idval)
cc  find min&max grid values and whether dval's are present
c        dimension a(n)
c        idval=0
c        amn=1.0e30
c        amx=-amn
c        do 10 i=1,n
c          if(a(i).ge.dval) then
c            idval=1
c          else
c            if(a(i).lt.amn) amn=a(i)
c            if(a(i).gt.amx) amx=a(i)
c          endif
c   10   continue
c        if (idval.eq.1) then
c          aset = amn - (amx-amn)/20.
c          do 20 i=1,n
c            if (a(i).ge.dval) a(i)=aset
c   20   continue
c        endif
c        return
c        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*(*)
	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 text(x,y,string,size,angle)
c  text driver where x,y are in plot units.
	dimension is(25)
	character s*100
	character string*(*)
	equivalence (is,s)
	data icode/3/, xoff,yoff/0.0,0.0/
	lenstr=len(string)
	do nch=lenstr,1,-1
	  ich=ichar(string(nch:nch))
	  if(ich.gt.32 .and. ich.lt.127) go to 10
	enddo
	return
10	s(1:nch)=string(1:nch)
	radian=1.7453292e-2*angle
	call vchar(x,y,is,nch,icode,size,radian,xoff,yoff)
	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 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
        icolor=300
	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 boxall
	common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),if2(50)
        common /plot/   dxp(2),dyp(2),xp(4),yp(4),dxpw(2),dypw(2),
     1               xpw(4),ypw(4),xpix,ypix,xscale,yscale
	character label*2
	if(icount.le.1) return
	size=.1
	do 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,size,0.0)
	enddo
	return
	end
c****************************************
        subroutine box(p1,p2)
	dimension p1(2),p2(2)
	do i = 1,2
	if(p1(i).gt.p2(i)) then
	  tmp=p1(i)
	  p1(i)=p2(i)
	  p2(i)=tmp
	endif
	enddo
	call line(p1(1),p1(2),1,0,0)
	call line(p1(1),p2(2),1,1,0)
	call line(p2(1),p2(2),1,1,0)
	call line(p2(1),p1(2),1,1,0)
	call line(p1(1),p1(2),1,1,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 pltseg
c  plot segments defining profile locations
	common /pparms/ iprf,nprf,nseg(10),sx(10,10),sy(10,10)
	do j = 1, nprf
	  if (nseg(j).gt.1) 
     1     call line (sx(1,j),sy(1,j), nseg(j), 0,0)
	enddo	    
	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 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)
      m1=2
      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
	return
        end
c****************************************
      subroutine tgin2(m1,m2,m3,m4)
c        print*,m1,m2,m3,m4,'tgin2 in'
      call mouse(m1,m2,m3,m4)
c        print*,m1,m2,m3,m4,'tgin2 out'
      return
      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 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)
        call tgin2(2,0,0,0)
	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
c        print*,ix,iy,ifunc,menux
	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)
        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 getseg(ipen,xlmt,ylmt,ifunc,ixo,iyo)
c get profile location 
c active main menu functions: 1=exit, 4=delete pt, 2,3,5,6,7=cancel
	common /pparms/ iprf,nprf,nseg(10),sx(10,10),sy(10,10)
	dimension xlmt(2),ylmt(2)
        character*79 prompt
	data ity,nsmax /5,10/
	ns = 0
cc        call tmode(1)
cc        call tclear(1)
cc        call tpromp(1,1,1,' ')
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
1       prompt=' enter 1 to generate profile '
        call tpromp(0,0,0,prompt)
        prompt=' enter 2 to read profile coordinates from a file '
        call tpromp(0,1,0,prompt)
        if (nprf.gt.0) then
          prompt=' enter 3 to activate previous profile '
          call tpromp(0,2,0,prompt)
        endif
	read(ity,*,err=1,end=1) ians
	if ( ians.lt.1 .or. ians.gt.3 .or. 
     1   (nprf.eq.0 .and. ians.eq.3) ) go to 1
c
	if (ians.le.2) then
	  nprf = nprf + 1 
	  iprf = nprf
	endif
	if (ians.eq.2) then
	  call rdseg(ns,sx(1,iprf),sy(1,iprf),nsmax,ierr)
	  if (ierr.eq.0 .and. ns.gt.1) 
     1     call line(sx(1,iprf),sy(1,iprf),ns,0,0)
	endif
	if (ians.eq.3) then
          prompt=' '
          call tpromp(0,2,0,prompt)
          write(prompt,15) nprf
15        format(' number of defined profiles =',i5,' ')
          call tpromp(0,0,0,prompt)
          prompt=' enter profile number to activate: '
          call tpromp(0,1,0,prompt)

	  read(ity,*,end=1,err=1) ip
	  if (ip.ge.1 .and. ip.le.nprf) then
	    iprf = ip
	    ns = nseg(iprf)
	    else
            print *,' not defined'
	    go to 1
	  endif
	endif
c
cc        call tmode(3)
cc        call tpromp(1,1,1,' Begin profile, hit "exit" to process')
cc        call tpromp(1,2,1,' or "delete pt" to erase. Any other')
cc        call tpromp(1,3,1,' menu function will cancel profile.')
        prompt=' Begin profile, hit "exit" to process or "delete pt" to
     1erase.  '
        call tpromp(0,0,0,prompt)
        prompt=' Any other menu function will cancel profile. '
        call tpromp(0,1,0,prompt)
c
10      call menudvr(tsx,tsy,ifunc,ixo,iyo)
cc        if (ns.eq.0) call tclear(1)
	go to (100,50,50,40,50,50,50) ifunc
c
c  store and draw segment
20	if(ns.eq.nsmax) go to 999
	if (tsx.lt.xlmt(1)) tsx = xlmt(1)
	if (tsx.gt.xlmt(2)) tsx = xlmt(2)
	if (tsy.lt.ylmt(1)) tsy = ylmt(1)
	if (tsy.gt.ylmt(2)) tsy = ylmt(2)
	ns = ns + 1
	sx(ns,iprf) = tsx
	sy(ns,iprf) = tsy
	if (ns.gt.1) then
	  call line(sx(ns-1,iprf),sy(ns-1,iprf),1,0,0)
	  call line(sx(ns,iprf),  sy(ns,iprf),  1,1,0)
	endif
	go to 10
c
c  erase last segment
40	if(ns.gt.1) then
	  call tdraw(1)
          call line(sx(ns-1,iprf),sy(ns-1,iprf),1,0,-7)
          call line(sx(ns,iprf),  sy(ns,iprf),  1,1,-7)
	  call tdraw(0)
	endif
	ns = ns - 1
	go to 10
c
c  erase profile, return
50	if (ns.gt.1) then
          call tdraw(1)
          call line( sx(1,iprf),sy(1,iprf), 1,    0,-7)
          call line( sx(2,iprf),sy(2,iprf), ns-1, 1,-7)
          call tdraw(0)
	endif
	ns = 0
	if (ians.lt.3) then
	  nprf = nprf - 1
	  iprf = nprf
	endif
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
cc        call tpromp(1,5,1,' function canceled')
        prompt=' function canceled '
        call tpromp(0,2,0,prompt)
cc        call waiter(1.0)
c
100	nseg(iprf) = 0
	if (ns.gt.1) nseg(iprf) = ns
cc999     call tclear(1)
999     return
	end
c****************************************
        subroutine rdseg(ns,sx,sy,nmax,ierr)
	dimension sx(nmax),sy(nmax)
	character ifile*56
        character*79 prompt
        data ity,iseg/5,20/
	ierr = 1
cc        call tmode(1)
cc        call tclear(1)
c        call tpromp(1,1,1,' ')
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
        call tpromp(0,2,0,prompt)
c
10      prompt=' enter input profile coordinate filename: '
        call tpromp(0,0,0,prompt)
	read(ity,37,end=10,err=10) ifile
37	format(a)
	open (iseg, file=ifile, status='old', form='formatted',
     1     err=777)
c
	do ns = 1, nmax
	  read (iseg,*,end=20,err=666) sx(ns),sy(ns)
	enddo
        write(prompt,39) nmax
39      format(' limit of',i6,' profile points ')
        call tpromp(0,0,0,prompt)
	go to 555
20	ns = ns - 1
c
555     write(prompt,21) ns
21      format(' profile has', i6, 'points ')
        call tpromp(0,1,0,prompt)
	close(iseg)
	ierr = 0
	go to 999
666     prompt=' error in data file '
        call tpromp(0,0,0,prompt)
        prompt=' file should consist of x,y coordinate pair records '
        call tpromp(0,1,0,prompt)
	close(iseg)
	go to 999
777     prompt=' file not found or error in filename '
        call tpromp(0,0,0,prompt)
        prompt=' try again ? '
        call tpromp(0,1,0,prompt)
	if (noyes(ldum).eq.1) go to 10
        prompt=' function canceled '
        call tpromp(0,2,0,prompt)
c
cc999     call waiter(2.0)
cc        call tclear(1)
cc        call tmode(3)
999     return
	end
c****************************************
        function noyes(idum)
	character inp*8
	icnt = 0
	noyes = -1
1	inp = ' '
	read (5,2,end=5,err=5) inp
2	format(a8)
	do n = 1,8
	  if (inp(n:n).ne.' ') then
	    if (inp(n:n).eq.'y' .or. inp(n:n).eq.'Y') noyes = 1
	    if (inp(n:n).eq.'n' .or. inp(n:n).eq.'N') noyes = 0
	    if (noyes.gt.-1) return
	    go to 5
	  endif
	enddo
5	if (icnt.lt.3) then
	  icnt = icnt + 1
          print  6
6	  format(' y or n:'$)
	  go to 1
	endif
        print *,' count exceeded, answering no'
	noyes = 0
	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 getprf(igrd,nc,nr,xo,dx,yo,dy, ns,sx,sy,
     1                 dp,nprf,px,py,pz,wrk)
c  interpolate profile from grid
	dimension sx(ns),sy(ns)
	dimension px(1),py(1),pz(1),wrk(1)
        character*79 prompt
	data ity/5/, npmax/500/
	if (ns.le.1) go to 999
cc        call tmode(1)
cc        call tpromp(1,1,1,' ')
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
        call tpromp(0,2,0,prompt)
c
	if(dp.eq.0.0) then
10        write(prompt,15) dx
15        format('grid interval = ',f10.5,' ')
          call tpromp(0,0,0,prompt)
19        prompt=' enter interpolating interval: '
20        call tpromp(0,1,0,prompt)
	  read(ity,*,err=19) dp
cc          type *,' interpolating interval =', dp
	  dp = abs(dp)
	  r = 0.0
	  do i = 1, ns-1
	    x = sx(i+1) - sx(i)
	    y = sy(i+1) - sy(i)
	    r = r + sqrt( x*x + y*y )
	  enddo
	  rdp = r/dp
	  np = 99999
	  if (rdp.lt.99999)  np = 2 + int( r/dp )
	  if (np.gt.npmax) then
            write(prompt,21) npmax
21          format(' maximum number of interpolated values is',i5,' ')
            call tpromp(0,0,0,prompt)
            write(prompt,22) r
22          format(' the profile length is',f10.5,' ')
            call tpromp(0,1,0,prompt)
            prompt=' do you want to change the interpolating interval? '
            call tpromp(0,2,0,prompt)
	    if (noyes(ldum).eq.1) go to 10
c            type *,' ok, will do what can'
	  endif
	endif
c
	ncm    = nc - 1
	nrm    = nr - 1
	lastr  = 0
	nprf   = 0
	dstart = 0.0
c
	do 200 iseg = 1, ns-1
	distx = sx(iseg+1) - sx(iseg)
	disty = sy(iseg+1) - sy(iseg)
	dist  = sqrt(distx**2 + disty**2)
	dpx   = distx / dist
	dpy   = disty / dist
	ipseg = 0
c
50	iok   = 0
	ipseg = ipseg + 1
	r     = dp*(ipseg-1) + dstart
	prfx  = r*dpx + sx(iseg)
	prfy  = r*dpy + sy(iseg)
	gx    = (prfx-xo) / dx
	gy    = (prfy-yo) / dy
	igx   = int(gx) + 1
	igy   = int(gy) + 1
	if (igx.le.1   .or. igy.le.1)   go to 100
	if (igx.ge.ncm .or. igy.ge.nrm) go to 100
	delx  = gx - float(igx-1)
	dely  = gy - float(igy-1)
	if (delx.gt.1.0 .or. delx.lt.0.0 .or.
     1   dely.gt.1.0 .or. dely.lt.0.0) then
          print *, 'getprf error:', delx,dely, prfx,prfy, igx,igy
	  go to 100
	endif
c
	if (igy.ne.lastr) then
	  iw = 1
	  loca = igy-1
	  do irow = loca, loca+3
	    call rowda(nc,wrk(iw),-1,irow,igrd,ierr)
	    if (ierr.ne.0) go to 100
	    iw = iw + nc
	  enddo
	  lastr = igy
	endif
	iok = 1
c
100	nprf = nprf + 1
	px(nprf) = prfx
	py(nprf) = prfy
	pz(nprf) = 1.0e-30
	if (iok.eq.1) call poly3(wrk,nc,4, igx,2, delx,dely,pz(nprf))
c
	if (nprf.eq.npmax) go to 999
	rtest = r + dp
	if (rtest.le.dist) go to 50
	dstart = rtest - dist
200	continue
c
cc999     call tmode(3)
999     return
	end

c****************************************
      subroutine poly3(y,nc,nr,igx,igy,dx,dy,yint)
      dimension t(5),y(1)
      dv=1.e30
      if(igx.le.0 .or. igy.le.0) go to 99
      if(igx.gt.nc .or. igy.gt.nr) go to 99
      if(igx.lt.2 .or. igy.lt.2) go to 88
      if(igx.gt.nc-1 .or. igy.gt.nr-1) go to 88
      ix=(igy-2)*nc+(igx-1)
      x=1.0+dx
      j=0
1     if(y(ix).ge.dv) go to 99
      if(y(ix+1).ge.dv) go to 99
      if(y(ix+2).ge.dv) go to 99
      if(y(ix+3).ge.dv) go to 99
      y1=y(ix)
      y2=y(ix+1)-y1
      y3=y(ix+2)-y1
      y4=y(ix+3)-y1
2     y2h=0.5*y2
      y3h=0.5*y3
      y2t=y2+y2+y2h
      a= y2h -y3h +.1666667*y4
      b= (y3+y3) -(y2t +0.5*y4)
      c= (y2t+y2h) -(y3+y3h) +.3333333*y4
      j=j+1
      t(j)= y1 +x*(c +x*(b +x*a))
      ix=ix+nc
      if(j-4) 1,3,4
3     y1=t(1)
      y2=t(2)-t(1)
      y3=t(3)-t(1)
      y4=t(4)-t(1)
      x=1.+dy
      go to 2
4     yint=t(5)
      return
88    ix=(igy-1)*nc+igx
      z1=(y(ix+1)-y(ix))*dx+y(ix)
      z2=(y(ix+nc+1)-y(ix+nc))*dx+y(ix+nc)
      yint=(z2-z1)*dy+z1
      return
99    yint=1.e30
      return
      end
c****************************************
        subroutine outprf(ns,sx,sy, dr,np,px,py,pz,pf, wrk)
	dimension sx(ns),sy(ns)
	dimension px(1),py(1),pz(1),pf(1),wrk(1)
	character ifile*56,ofile*56, id2*56,ipgm2*8
        character*79 prompt
        data ity, idat, idat2, nsmax /5, 20, 21, 10/
        data dv/1.e30/
cc        call tmode(1)
cc        call tclear(1)
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
        call tpromp(0,2,0,prompt)
37	format(a)
	ifield = 1
	np2 = np
c
	if(ns.le.1)then
	  call rdseg(ns,sx,sy,nsmax,ierr)
	  if (ierr.eq.1) go to 999
	  ifield = 0
	endif
30 	if (dr.le.0.0) then
          prompt=' enter interpolation interval: '
          call tpromp(0,0,0,prompt)
	  read(ity,*,err=30,end=30) dr
	  ifield = 0
	endif
c
100	if (ifield.eq.0) then
          prompt=' enter field grid filename: '
          call tpromp(0,1,0,prompt)
	  read(ity,37) ifile	
	  open (idat, file=ifile, status='old', form='unformatted',
     1     err=777)
	  read(idat) id2,ipgm2, nc2,nr2,nz2, xo2,dx2,yo2,dy2
          open (idat2, status='scratch', err=777,
     1   form='unformatted', access='direct', recl=nc2*4)
	  do j = 1, nr2
	    call rowio(nc2,wrk,-1,idat,idat,iend)
	    if (iend.ne.0) go to 555
	    call rowda(nc2,wrk,0,j,idat2,ierr)
	  enddo
	  call getprf(idat2,nc2,nr2,xo2,dx2,yo2,dy2, ns,sx,sy,
     1           dr,np2,px,py,pf, wrk)
	  close(idat)
	  close(idat2)
          ifield = 1
	endif
c
        prompt=' do you have an observation surface grid ? '
        call tpromp(0,0,0,prompt)
	if(noyes(ldum).eq.0) then
          prompt=' enter a constant for the observation elevation: '
          call tpromp(0,1,0,prompt)
          read(ity,*,err=100,end=100) zc
	  do i = 1, np2
            pz(i) = -zc
	  enddo
          const=0.0
          conv=1.0
        else
          prompt=' enter filename: '
          call tpromp(0,1,0,prompt)
          read(ity,37) ifile
	  open (idat, file=ifile, status='old', form='unformatted',
     1     err=777)
	  read(idat) id2,ipgm2, nc2,nr2,nz2, xo2,dx2,yo2,dy2
          open (idat2, status='scratch', err=777,
     1     form='unformatted', access='direct', recl=nc2*4)
	  do j = 1, nr2
	    call rowio(nc2,wrk,-1,idat,idat,iend)
	    call rowda(nc2,wrk,0,j,idat2,ierr)
	  enddo
	  call getprf(idat2,nc2,nr2,xo2,dx2,yo2,dy2, ns,sx,sy,
     1           dr,np3,px,py,pz, wrk)
	  close(idat)
          close(idat2)
          prompt=' enter grid units/elevation units: '
          call tpromp(0,0,0,prompt)
          read(ity,*) conv
          conv=-conv
          prompt=' enter additive constant in elevation units: '
          call tpromp(0,1,0,prompt)
          read(ity,*) const
        endif
c
        write(prompt,101) np2,np3
101     format(' number of field values and elevations',i6,i6,' ')
        call tpromp(0,0,0,prompt)
        prompt=' enter an output filename: '
        call tpromp(0,1,0,prompt)
	read(ity,37) ofile
	open (idat, file=ofile, status='new', form='formatted',
     1     err=888)
        prompt=' do you want to output dvals? '
        call tpromp(0,2,0,prompt)
        iodv=noyes(ldum)
cc110     type *,' coordinate output type, enter either '
cc        type *,' 1 for distance along profile, elevation, field value '
cc        type *,' 2 for x,y pairs, elevation, field'
cc        type *,' 3 for distance, elev, field, x,y pairs'
cc        read(ity,*,err=110,end=110) ifun
cc        if(ifun.lt.1 .or. ifun.gt.3) go to 110
120	format(1p5g13.5)
c
cc        if(ifun.eq.1) then
cc          do i = 1, np2
cc            r = dr*float(i-1)
cc            write(idat,120) r,pz(i),pf(i)
cc          enddo
cc          else if(ifun.eq.2) then
cc          do i = 1, np2
cc            write(idat,120) px(i),py(i),pz(i),pf(i)
cc          enddo
cc          else
	  do i = 1, np2
            if(iodv.eq.1.or.(pz(i).lt.dv.and.pf(i).lt.dv)) then
              r = dr*float(i-1)
              zp=(pz(i)+const)*conv
              write(idat,120) r,zp,pf(i),px(i),py(i)
            endif
          enddo
cc        endif
	close(idat)
	ifield = 0
c
        prompt=' '
        call tpromp(0,2,0,prompt)
        prompt=' do you want the same profile location '
        call tpromp(0,0,0,prompt)
        prompt=' with different field data? '
        call tpromp(0,1,0,prompt)
        if (noyes(ldum).eq.1) go to 100
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
        prompt=' do you want to save profile coordinates? '
        call tpromp(0,0,0,prompt)
        if (noyes(ldum).eq.1) call wtseg(ns,sx,sy)
	go to 1000
c
555     prompt=' error in grid file '
        call tpromp(0,2,0,prompt)
	close(idat)
	close(idat2)
	go to 999
777     prompt=' could not open input file '
        call tpromp(0,2,0,prompt)
	go to 999
888     prompt=' could not open output file '
        call tpromp(0,2,0,prompt)
c
cc999     call waiter(2.0)
cc        call tclear(1)
cc        call tmode(3)
999     continue
1000	return
	end
c****************************************
        subroutine wtseg(ns,sx,sy)
	dimension sx(1),sy(1)
	character ofile*56
        character*79 prompt
	if (ns.le.1) return
        prompt=' enter output profile coordinate filename: '
        call tpromp(0,1,0,prompt)
	read(5,1) ofile
1	format(a)
	open (22, file=ofile, status='new', form='formatted', err=888)
	do i = 1, ns
	  write(22,*,err=777) sx(i),sy(i)
	enddo
	close(22)
	go to 999
777     print *,' error writing output file'
	close(22)
	go to 999
888     print *,' error opening output file'
999	return
	end

c****************************************
        subroutine expand(ipen,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,iflip,jflip
	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
5	call tdraw(1)
	call box(p1,p2)
	call tdraw(0)
cc10      call tpromp(1,1,1,' lower-left corner (''exit'' to cancel)')
        prompt=' lower-left corner (''exit'' to cancel) '
10      call tpromp(0,0,0,prompt)
        call menudvr(p1(1),p1(2),ifunc,ixo,iyo)
	if(ifunc.gt.0) then
cc          call tpromp(1,1,1,' function canceled')
          prompt=' function canceled '
          call tpromp(0,0,0,prompt)
	  go to 999
	endif
cc15      call tpromp(1,1,1,' upper-right corner')
        prompt=' upper-right corner '
15      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)
cc        call tpromp(1,1,1,' ''exit'' if ok')
        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 i=1,2
	  p1(i)=0.0
	  p2(i)=0.0
	enddo
	return
	end
c****************************************
        subroutine contrk(iret)
        common /wrkgrd/ iwrkg,icount,c1(2,10),c2(2,10),nskip(10),
     1               icol(10),irow(10),ncw(10),nrw(10)
	character prompt*79
cc        call tmode(1)
cc        call tclear(1)
cc1       call tpromp(1,1,1,' ')
1       call tpromp(0,0,0,' ')
        write(prompt,20) icount,iwrkg
20      format(' number of active grids is',i5,' present grid number',
     1         i5,' ')
        call tpromp(0,0,0,prompt)
        prompt=' enter desired grid number (Gn < 0, will reset active gr
     1ids to abs(Gn)): '
        call tpromp(0,1,0,prompt)
	read (5,*,end=99,err=99) iw2
	if(iabs(iw2).gt.10 .or. iw2.gt.icount) go to 1
        prompt=' '
        call tpromp(0,1,0,prompt)
	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,20) icount
          call tpromp(0,0,0,prompt)
	endif
	iwrkg=iw2
	iret=1
	go to 99
c
88      prompt=' no change '
        call tpromp(0,0,0,prompt)
	iret=2
cc        call waiter(1.0)
cc99      call tclear(1)
cc        call tmode(3)
99      return
	end
c****************************************
        subroutine outgrd
	common /dfile/  igrd,   iout,   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/  wrk(40000)
	character ofile*56
	character prompt*79
	equivalence (ibf,ofile)
	data ity/5/
c
cc        call tmode(1)
cc        call tpromp(1,1,1,' > ')
c        call tpromp(0,0,0,' > ')
c
	if (icount.eq.1) then
          prompt=' no subgrids defined '
          call tpromp(0,0,0,prompt)
	  go to 999
	endif
c
9       write(prompt,10) icount
10      format(' subgrids available are 2 thru',i3,' ')
        call tpromp(0,0,0,prompt)
        prompt=' which do you want written out? '
        call tpromp(0,1,0,prompt)
	read(ity,*,err=9,end=9) isel
	if (isel.lt.2 .or. isel.gt.icount) go to 9
c
        prompt=' enter filename: '
        call tpromp(0,2,0,prompt)
	read(ity,37) ofile
37	format(a)
        prompt=' '
        call tpromp(0,0,0,prompt)
        call tpromp(0,1,0,prompt)
        call tpromp(0,2,0,prompt)
	open (iout, file=ofile, status='new', form='unformatted',
     1     err=777)
c
	xo2 = c1(1,isel)
	dx2 = dx
	if (iflip.eq.1) then
	  xo2 = -xo2
	  dx2 = -dx2
	endif
	yo2 = c1(2,isel)
	dy2 = dy
	if (jflip.eq.1) then
	  yo2 = -yo2
	  dy2 = -dy2
	endif
	write(iout) id,ipgm, ncw(isel),nrw(isel),nz, xo2,dx2,yo2,dy2
c
	j2 = irow(isel)
	i2 = icol(isel)
	do j = 1, nrw(isel)
	  call rowda(nc,wrk,-1,j2,igrd,ierr)
	  if (ierr.ne.0) go to 888
	  call rowio(ncw(isel),wrk(i2),0,iout,iout,iend)
	  j2 = j2 + 1
	enddo
c
	close(iout)
	go to 999
c
777     prompt=' could not open output file '
        call tpromp(0,0,0,prompt)
	go to 999
888     write(prompt,88) igrd
88      format(' error reading grid file, switch number ',i5,' ')
        call tpromp(0,0,0,prompt)
        write(prompt,89) j-1
89      format(' exiting after writing', i5,' rows ')
        call tpromp(0,1,0,prompt)
c
cc999     call tclear(1)
cc        call tmode(3)
999     return
	end
c****************************************
        subroutine tclear(iop1)
c  clear screen
        integer*2 jop1
        jop1=iop1
      call setmod(jop1)
	return
        end

c****************************************
c****************************************
