subroutine plot3d(f,x,y,nx,ny,xw,yw,xview,yview,zview,line,scale)
c
	dimension f(nx,ny),x(nx,ny),y(nx,ny)
c
	common /hwindo/ xscr,yscr,ixorig,iyorig,iascr,ibscr,icscr,idscr
	common/pscale/xv,yv,dx,dy,ca,sa,ly,lx,xscale,xn,yscale,yn
	common/color/isubch
c
	data small,pt5/1.e-06,.5/
c
	if(line.le.2)goto 10
	line=line-3
	goto 305
c
c  normalize constants
c
10	scrx = xscr
	scry = yscr
	xb=1.
	yb=yw/xw
	xv=xview/xw
	yv=yview/xw
	h=zview/xw
	ir=1
	nxx=nx+1
	nyy=ny+1
c
c  determine projection plane
c
	a=atan2(.5*xb-xv,.5*yb-yv)
	sa=sin(a)
	ca=cos(a)
c
c  y component of the vertices of the projected grid
c
	y2=xb*sa
	y3=yb*ca
	y4=y2+y3
	yc=xv*sa+yv*ca
	yn=amin1(0.,y2,y3,y4)
	if(yc.ge.yn)yc=yn-.01*abs(yb)
	yc=-yc
	xc=xv*ca-yv*sa
	dx=xb/(nx-1)
	dxx=1./dx
	vx=dxx*xv+1.00001
	lx=vx
	dy=yb/(ny-1)
	dyy=1./dy
	vy=dyy*yv+1.00001
	ly=vy
	yn=0.
	yx=yn
	d=.5*y4
	yd=yc+d
	yp=0.
c
c  determine minimum and maximum of the surface
c
	tn=f(1,1)
	tx=tn
	do 50 j=1,ny
	do 50 i=1,nx
	tn=amin1(tn,f(i,j))
50	tx=amax1(tx,f(i,j))
c
c  compute scaling factor
c
	s=tx-tn
	if(s.eq.0.)s=tx
	s=amax1(xb,yb)/s
	if(scale.ne.0.)s=scale
c
c  compute projected x and y coordinates of surface
c
	do 120 j=1,ny
	xp=0.
	c=yp*sa+xc
	do 100 i=1,nx
	t=f(i,j)*s
	hl=xp*sa+yp*ca
	x(i,j)=(xp*ca-c)*yd/(hl+yc)+small
	y(i,j)=((hl-d)*h+t*yd)/(hl+yc)+small
	yn=amin1(yn,y(i,j))
	yx=amax1(yx,y(i,j))
100	xp=xp+dx
120	yp=yp+dy
c
c  x component of the vertices of the grid
c
	x1=x(1,1)
	x2=x(nx,1)
	x3=x(1,ny)
	x4=x(nx,ny)
	tn=tn*yd*s*1.1
c
c  y component of the vertices of the projected grid
c
	y1=(tn-d*h)/yc
	y2=((y2-d)*h+tn)/(y2+yc)
	y3=((y3-d)*h+tn)/(y3+yc)
	y4=((y4-d)*h+tn)/(y4+yc)
c
	xn=amin1(x1,x2,x3,x4)
	xx=amax1(x1,x2,x3,x4)
	yn=amin1(yn,y1,y2,y3,y4)
	xscale=scrx*xw/(xx-xn)
	yscale=scry*yw/(yx-yn)
	xn=xscale*xn
	yn=yscale*yn
c
c  convert to screen cordinates
c
	do 200 j=1,ny
	do 200 i=1,nx
	x(i,j)=aint(xscale*x(i,j)-xn+pt5)
200	y(i,j)=aint(yscale*y(i,j)-yn+pt5)
c
	x1=aint(xscale*x1-xn+pt5)
	y1=aint(yscale*y1-yn+pt5)
	x2=aint(xscale*x2-xn+pt5)
	y2=aint(yscale*y2-yn+pt5)
	x3=aint(xscale*x3-xn+pt5)
	y3=aint(yscale*y3-yn+pt5)
	x4=aint(xscale*x4-xn+pt5)
	y4=aint(yscale*y4-yn+pt5)
	xscale=xscale*yd
c
c  begin hidden line removal
c
	do 300 j=1,ny
	do 300 i=1,nx
	ix=isign(1,lx-i)
	iy=isign(1,ly-j)
	ih=i
	iplot=1
	if(0.eq.lx-i)goto 250
215	ih=ih+ix
	if(ih.lt.1.or.ih.gt.nx)goto 250
	xp=(ih-1)*dx
	xm=dx/dy*(dy*(j-1)-yv)/(dx*(i-1)-xv)
	yp=xm*(ih-i)+j
	jj=ifix(yp)
	jh=jj+1
	if(jj.lt.1.or.jh.gt.ny)goto 250
	yp=(yp-1)*dy
	hl=yc+xp*sa+yp*ca
	xi=aint(xscale*(xp*ca-yp*sa-xc)/hl-xn+pt5)
	yi=(y(ih,jj)-y(ih,jh))*(xi-aint(x(ih,jh)))/
	1(aint(x(ih,jj))-aint(x(ih,jh)))+y(ih,jh)
	if(aint(yi+pt5).lt.y(i,j))goto 215
	iplot=0
	goto 290
250	if(ly-j.eq.0)goto 290
	jh=j
260	jh=jh+iy
	if(jh.lt.1.or.jh.gt.ny)goto 290
	yp=(jh-1)*dy
	ym=dy/dx*(dx*(i-1)-xv)/(dy*(j-1)-yv)
	xp=ym*(jh-j)+i
	ii=ifix(xp)
	ih=ii+1
	if(ii.lt.1.or.ih.gt.nx)goto 290
	xp=(xp-1)*dx
	hl=yc+xp*sa+yp*ca
	xi=aint(xscale*(xp*ca-yp*sa-xc)/hl-xn+pt5)
	yi=(y(ii,jh)-y(ih,jh))*(xi-aint(x(ih,jh)))/
	1(aint(x(ii,jh))-aint(x(ih,jh)))+y(ih,jh)
	if(aint(yi+pt5).lt.y(i,j))goto 260
	iplot=0
290	x(i,j)=x(i,j)+.101*iplot
300	continue
c
c  plot x grid lines
c
305	call window( 0., xw, 0., yw )
	call vuport( 0., scrx * xw, 0., scry * yw )
	if(line.eq.2)goto 350
	isubch=1
	do 310 j=1,ny
	iplot0=ifix(10*(x(1,j)-ifix(x(1,j))))
	if(iplot0.eq.1)call moveto(aint(x(1,j)),aint(y(1,j)),0,0)
	do 310 i=2,nx
	iplot=ifix(10*(x(i,j)-aint(x(i,j))))
	call insect(i,j,iplot,iplot0,x,y,nx,ny,0)
310	iplot0=iplot
c
c  plot visible verticle bars
c
	if(lx.ge.1.and.lx.le.nx)goto 350
	if(lx.lt.1)goto 320
	ii=nx
	aa=(y2-y4)/(x2-x4)
	bb=y4-x4*aa
	ix=x2
	iy=y2
	goto 330
320	ii=1
	aa=(y3-y1)/(x3-x1)
	bb=y1-x1*aa
	ix=x1
	iy=y1
330	do 340 j=1,ny
	call moveto(aint(x(ii,j)),aint(y(ii,j)),0,0)
340	call moveto(aint(x(ii,j)),aint(aa*x(ii,j)+bb),1,0)
	call moveto(float(ix),float(iy),1,0)
c
c  plot y grid lines
c
350	if(line.eq.1)return
	isubch=2
	do 400 i=1,nx
	iplot0=ifix(10*(x(i,1)-ifix(x(i,1))))
	if(iplot0.eq.1)call moveto(aint(x(i,1)),aint(y(i,1)),0,0)
	do 400 j=2,ny
	iplot=ifix(10*(x(i,j)-aint(x(i,j))))
	call insect(i,j,iplot,iplot0,x,y,nx,ny,1)
400	iplot0=iplot
c
c  plot visible verticle bars
c
	if(ly.ge.1.and.ly.le.ny)return
	if(ly.lt.1)goto 420
	jj=ny
	aa=(y3-y4)/(x3-x4)
	bb=y4-x4*aa
	ix=x3
	iy=y3
	goto 430
420	jj=1
	aa=(y2-y1)/(x2-x1)
	bb=y1-x1*aa
	ix=x1
	iy=y1
430	do 440 i=1,nx
	call moveto(aint(x(i,jj)),aint(y(i,jj)),0,0)
440	call moveto(aint(x(i,jj)),aint(aa*x(i,jj)+bb),1,0)
	call moveto(float(ix),float(iy),1,0)
	return
	end

	subroutine insect(i1,j1,iplot,iplot0,x,y,nx,ny,l)
c
	dimension x(nx,ny),y(nx,ny)
	byte eq1,eq2
	common /hwindo/ scrx,scry,ixorig,iyorig,iascr,ibscr,icscr,idscr
	common/pscale/xv,yv,dx,dy,ca,sa,ly,lx,xscale,xn,yscale,yn
	data pt5/.5/
	if(iplot.eq.iplot0)goto 290
	i=i1
	j=j1
	i0=i1
	j0=j1
	ix=isign(1,lx-i)
	iy=isign(1,ly-j)
	insct=0
	if(l.eq.0)goto 10
	j=j1-iplot0
	j0=j1-iplot
	xf=aint(x(i0,j0))
	yf=y(i0,j0)
	xl=aint(x(i,j))
	yl=y(i,j)
	imin=i1
	jmin=iy*min0(iy*j,iy*j0)
	if(jmin.eq.j0.and.(y(i+ix,j)-yl)/(x(i+ix,j)-xl).lt.
	1(yf-yl)/(xf-xl))goto 270
	goto 20
10	i=i1-iplot0
	i0=i1-iplot
	xf=aint(x(i0,j0))
	yf=y(i0,j0)
	xl=aint(x(i,j))
	yl=y(i,j)
	imin=ix*min0(ix*i,ix*i0)
	jmin=j1
	if(imin.eq.i0.and.(y(i,j+iy)-yl)/(x(i,j+iy)-xl).gt.
	1(yf-yl)/(xf-xl))goto 270
20	xint=xf-xl
	yint=yf-yl
	ih=imin
	if(lx-i.eq.0)goto 250
	xm1=dx/dy*(dy*(j-1)-yv)/(dx*(i-1)-xv)
	xm2=dx/dy*(dy*(j0-1)-yv)/(dx*(i0-1)-xv)
215	ih=ih+ix
	if(ih.lt.1.or.ih.gt.nx)goto 250
	yp1=xm1*(ih-i)+j
	if(i.eq.ih)yp1=yp1+iy
	yp2=xm2*(ih-i0)+j0
	if(i0.eq.ih)yp2=yp2+iy
	jn=ifix(amin1(yp1,yp2))
	jx=ifix(amax1(yp1,yp2))
	do 220 jj=jn,jx
	jh=jj+1
	if(jj.lt.1.or.jh.gt.ny)goto 220
	if((x(ih,jh)-aint(x(ih,jh)))+(x(ih,jj)-aint(x(ih,jj))).eq.0)
	1goto 220
	s1=(yf-yl)/(xf-xl)
	s2=(y(ih,jj)-y(ih,jh))/(aint(x(ih,jj))-aint(x(ih,jh)))
	if(s1.eq.s2)goto 220
	xi=(y(ih,jh)-yl+xl*s1-aint(x(ih,jh))*s2)/(s1-s2)
	if(abs(aint(xi+pt5)-.5*(xl+xf)).gt..5*abs(xf-xl))goto 220
	xmid=.5*aint(x(ih,jj)+x(ih,jh))
	xdif=.5*abs(aint(x(ih,jj))-aint(x(ih,jh)))
	if(abs(aint(xi+pt5)-xmid).gt.xdif)goto 220
	jh=(jj+jh+iy*(jj-jh))/2
	yi=s1*(xi-xl)+yl
	if(xint*(y(ih,jh+iy)-y(ih,jh))*ix*iy.lt.
	1yint*(aint(x(ih,jh+iy))-aint(x(ih,jh)))*ix*iy)goto 220
	xf=aint(xi+pt5)
	yf=aint(yi+pt5)
	insct=1
	if(xf-xl.eq.0)goto 270
220	continue
	goto 215
250	if(ly-j.eq.0)goto 270
	ym1=dy/dx*(dx*(i-1)-xv)/(dy*(j-1)-yv)
	ym2=dy/dx*(dx*(i0-1)-xv)/(dy*(j0-1)-yv)
	jh=jmin
260	jh=jh+iy
	if(jh.lt.1.or.jh.gt.ny)goto 270
	xp1=ym1*(jh-j)+i
	if(j.eq.jh)xp1=xp1+ix
	xp2=ym2*(jh-j0)+i0
	if(j0.eq.jh)xp2=xp2+ix
	jn=ifix(amin1(xp1,xp2))
	jx=ifix(amax1(xp1,xp2))
	do 265 ii=jn,jx
	ih=ii+1
	if(ii.lt.1.or.ih.gt.nx)goto 265
	if((x(ih,jh)-aint(x(ih,jh)))+(x(ii,jh)-aint(x(ii,jh))).eq.0)
	1goto 265
	eq1=i.eq.ih.or.i.eq.ii
	if(eq1.and.eq2)goto 265
	s1=(yf-yl)/(xf-xl)
	s2=(y(ii,jh)-y(ih,jh))/(aint(x(ii,jh))-aint(x(ih,jh)))
	if(s1.eq.s2)goto 265
	xi=(y(ih,jh)-yl+xl*s1-aint(x(ih,jh))*s2)/(s1-s2)
	if(abs(aint(xi+pt5)-.5*(xl+xf)).gt..5*abs(xf-xl))goto 265
	xmid=.5*aint(x(ii,jh)+x(ih,jh))
	xdif=.5*abs(aint(x(ii,jh))-aint(x(ih,jh)))
	if(abs(aint(xi+pt5)-xmid).gt.xdif)goto 265
	ih=(ih+ii+ix*(ii-ih))/2
	yi=s1*(xi-xl)+yl
	if(xint*(y(ih+ix,jh)-y(ih,jh))*ix*iy.gt.
	1yint*(aint(x(ih+ix,jh))-aint(x(ih,jh)))*ix*iy)goto 265
	xf=aint(xi+pt5)
	yf=aint(yi+pt5)
	insct=1
	if(xf.eq.xl)goto 270
265	continue
	goto 260
270	if(insct.eq.1)goto 280
	if(l.eq.0.and.i.eq.1.or.i.eq.nx)return
	if(l.eq.1.and.j.eq.1.or.j.eq.ny)return
	call moveto(aint(xl),aint(yl),0,0)
	return
280	ix=ifix(xf)
	iy=ifix(yf)
	call moveto(float(ix),float(iy),iplot0,0)
	if(xf.eq.xl)return
290	if(iplot.eq.0)return
	ix=ifix(x(i1,j1))
	iy=ifix(y(i1,j1))
	call moveto(float(ix),float(iy),1,0)
	return
	end