c*****PROGRAM PLOT3D
c*****THE NORMAL CONVENTION OF X VALUES REFERING TO COLUMNS
c*****AND Y VALUES TO ROWS IS REVERSED IN THIS PROGRAM
c*****INPUT OF ROWS AND COLUMNS,HOWEVER IS CONVENTIONAL
      parameter (maxd=90000)
      character*1 resp,resp2
      character*8 buff,passwrd
      character*56 ifile
      character*64 title
	dimension work(maxd),idim(1)
       logical opaque,block,array
       integer startx,skipx,dimx,starty,skipy,dimy
	common/pl3/amax,amin,startx,skipx,dimx,starty,skipy,dimy,
     & theta,phi,opaque,block,array
      data amin/1.0e38/,amax/-1.0e38/,dval/1.0e37/,
     & idim(1)/0/,sdval/-1.0e30/
       data passwrd/'sB8yr0hF'/

       call getarg(1,buff,length)
       if(buff.ne.passwrd) then
        print *,'program can only be run from program GEOCON'
        stop
       endif
c*****OPEN INPUT FILE AND READ HEADER OF DATA
1     write(*,'(a\)') ' enter input grid file name: '
      read(*,'(a)') ifile
      call opngrd(10,ifile,ierr,*1,*100)
      read(10) title,nc,nr,nz,xo,dx,yo,dy
	nvalue = nc *  nr
      if(nvalue.gt.maxd) then
        print *,'no. of columns * no. of rows > ',maxd
        pause 'press ENTER to return to menu'
        stop
      endif

      write(*,'(a)') ' enter viewing angle from the horizon (0-90)'
      write(*,'(a\)') ' 60 is a good start:'
      read *,thedeg
      if(thedeg.lt.0..or.thedeg.gt.90.) thedeg=60.
      thedeg = 90. - thedeg
      write(*,'(a)') ' enter horizontal viewing azimuth (0-360)'
      write(*,'(2a\)') ' 0 is north (top of the grid) and values',
     & ' are clockwise:'
      read *,phideg
      if(phideg.lt.0..or.phideg.gt.360.) phideg=135.
      phideg =  180. - phideg
      mode=3
      write(*,'(a\)') ' want to plot the entire data set ? [y]'
      read(*,'(a)') resp
      if(resp.eq.'n'.or.resp.eq.'N') then
        write(*,'(2a\)') ' enter row subset indexes (start,end',
     &  ' & increment):'
        read *,startx,dimx,skipx
        write(*,'(2a\)') ' enter column subset indexes (start,end',
     &  ' & increment):'
        read *,starty,dimy,skipy
        if(startx.le.0.or.startx.gt.nr) startx=1
        if(dimx.le.0.or.dimx.gt.nr) dimx=nr
        if(skipx.le.0) skipx = 1
	 if(mod((dimx-startx),skipx).ne.0) 
     &   dimx = (dimx-startx)/skipx * skipx + startx
c*****REVERSE THE ROW INDEXES
        stx = nr - dimx
        dimx = nr - startx
        startx = stx
        if(starty.le.0.or.starty.gt.nc) starty=1
        if(dimy.le.0.or.dimy.gt.nc) dimy=nc
        if(skipy.le.0) skipy = 1
	 if(mod((dimy-starty),skipy).ne.0) 
     &   dimy = (dimy-starty)/skipy * skipy + starty
      else
        startx=1
        dimx=nr
        skipx=1
        starty=1
        dimy=nc
        skipy=1
      endif
c*****READ DATA - REVERSING THE ROW ORDER
c*****FIND MIN AND MAX OF DATA AND SET DVAL'S TO -DVAL'S
	 inx = nvalue - nc + 1
        do 20 i = nr,1,-1
	 call r3d(work(inx),nc,10)
        do 10 j = inx,inx+nc-1
        if(work(j).gt.dval) then
          work(j) = -dval
          go to 10 
        endif
        if(work(j).lt.amin) amin=work(j)
        if(work(j).gt.amax) amax=work(j)
10      continue
        inx = inx - nc
20      continue
c*******NORMALIZE IF NECESSARY
	 if(amin.lt.0.) then
          do 30 i = 1,nvalue
          if(work(i).lt.sdval) go to 30
          work(i) = work(i) + abs(amin) 
30        continue
          amax = amax + abs(amin)
          amin = 0.
        endif

       block = mode.eq.2.or.mode.eq.4
       opaque = mode.eq.3.or.mode.eq.4
       array=nc.ne.1
40     theta=(3.1416/180.)*thedeg
	phi =(3.1416/180.)*phideg
	call ptt3d(work,nc,nr,resp)
	call endpt(idim)
      write(*,'(a\)') ' want to change viewing angles ? [y]'
      read(*,'(a)') resp2
      if(resp2.eq.'n'.or.resp2.eq.'N') go to 90
      write(*,'(2a\)') ' enter new viewing angle from the',
     & ' horizon (0-90):'
      read *,thedeg
      if(thedeg.lt.0..or.thedeg.gt.90.) thedeg=60.
      thedeg = 90. - thedeg
      write(*,'(a)') ' enter new horizontal viewing azimuth (0-360)'
      write(*,'(2a\)') ' 0 is north (top of the grid) and values',
     & ' are clockwise:'
      read *,phideg
      if(phideg.lt.0..or.phideg.gt.360.) phideg=135.
      phideg =  180. - phideg
      go to 40
90    stop
100   pause 'press ENTER to return to menu'
      stop
      end
c***************************************
      subroutine opngrd(igrid,ifile,ierr,*,*)
      character*(*) ifile
      character*64 title
      logical opn
        inquire(unit=igrid,opened=opn)
          if(opn) close(igrid)
        open(unit=igrid,file=ifile,
     &  status='old',form='unformatted',err=442,iostat=icheck)
        read(igrid,err=442,end=442,iostat=icheck) title,ncol
        rewind igrid
        return
442     call fcheck(igrid,ifile,icheck,ierr)
        if(ierr.eq.0) return 1
        if(ierr.eq.2) then
          ierr = 0
          return
        endif
        if(ierr.ne.0) return 2
c           if(ncol.gt.10000) then
c             open(unit=igrid,file=ifile,access='sequential',
c     &       status='old',form='unformatted')
c             access='sequen'
c           else
c             nbytes=(ncol+1)*4
c             open(unit=igrid,file=ifile,access='direct',
c     &       recl=nbytes,status='old',form='unformatted')
c             access='direct'
       return
       end
c***************************************
      subroutine fcheck(igrid,fname,ierror,ierr)
      character*(*) fname
      character*1 resp
      logical exists
      inquire(file=fname,exist=exists)
        if(.not.exists) then
          print *,'can''t find file ',fname
          fname=' '
          write(*,'(a\)') ' TRY AGAIN ? [y]: '
          read(*,'(a)') resp
          if(resp.eq.'n'.or.resp.eq.'N') then
            ierr=1
            return
          else
            ierr=0
            return
          endif
        endif
      call iof(igrid,fname,ierror,ierr)
      if(ierr.eq.1) go to 10
      return
10    print *,'error no. ',ierror,' with file: ',fname
      fname=' '
      return
      end
c***************************************************************
      subroutine iof(igrid,fname,ierror,ierr)
      parameter (numz=4500)
      character*(*) fname
      character*150 string
      dimension z(4500)
      close(igrid)
      jgrid=igrid+1
      kgrid=jgrid+1
      open(jgrid,file=fname,status='old',err=90,iostat=ierror,
     & form='formatted')
      read(jgrid,'(a)',err=90,iostat=ierror) string
c*****FORMATTED INPUT DATA FILE
c*****CREATE A SCRATCH UNFORMATTED FILE FOR INPUT TO PROGRAM
      open(igrid,status='scratch',form='unformatted')
      open(kgrid,status='scratch',form='formatted')
      lens = len_trim(string)
      write(kgrid,'(4a)') '''',string(1:64),'''',string(65:lens)
      rewind kgrid
      read(kgrid,*,err=90,iostat=ierror) string(1:64),nc,nr,nz,
     & xo,dx,yo,dy
      close(kgrid)
      if(nc.eq.0.or.nr.eq.0.or.nz.ne.1) then
        print *,'grid file has wrong format'
        go to 90
      endif
      if(nc.gt.numz) then
        print *,'number of input columns ',nc,' > ',numz
        go to 90
      endif
      write(igrid) string(1:64),nc,nr,nz,xo,dx,yo,dy
      iend=0
10    call iorw(jgrid,igrid,z,nc,iend,ierror,ierr)
      if(iend.ne.0) go to 20
      if(ierr.eq.1) go to 20
      go to 10
20    close(jgrid)
      rewind igrid
      return
90    ierr = 1
      close(jgrid)
      close(kgrid)
      return
      end
c***************************************
      subroutine iorw(jgrid,igrid,z,nc,iend,ierror,ierr)
      dimension z(nc)
      read(jgrid,*,end=10,err=20,iostat=ierror) y,z
      write(igrid) y,z
      return
10    iend=1
      ierr=2
      return
20    ierr = 1
      return
      end
c***********************************************************************
	subroutine r3d(z,nc,in)
	dimension z(nc)
	read(in) dum,z
       return
	end
c***********************************************************************
	subroutine ptt3d(a[huge],nc,nr,resp)
	character*1 resp
	common/pl3/amax,amin,startx,skipx,dimx,starty,skipy,dimy,
     & theta,phi,opaque,block,array
	dimension dxp(2),dyp(2),xp(4),yp(4)
       dimension a(nc,nr),sv(500),sh(500)
	integer icolor(500)
       integer startx,skipx,dimx,starty,skipy,dimy
	logical see,saw,opaque,block,array,scanx
       character*1 ctype
       common/captur/icapt,iwhite,ctype
	equivalence (m1,j1),(m12,j12),(j13,m13),(m2,j2),(mc,jc),(mp,jp)
	data sdval/-1.0e30/
c*****SET UP PLOTTER
      call pltset(9,xbrd,ybrd,isl)
      xp(1) = xbrd - 1.2
      xp(3) = .1
      xp(4) = xp(1) + xp(3) + 1.1
      yp(1) = ybrd - 1.2
      yp(3) = .1
      yp(4) = yp(1) + yp(3) + 1.1
       size = 10.
       scalez = size
	scalex=dimx-startx
	scaley=dimy-starty
	if(scalex.ge.scaley) go to 1 
	scalex=size*scalex/scaley
	scaley = size
	go to 2
1      scaley=size*scaley/scalex
	scalex = size
2 	tdimx=dimx-startx
	tdimy=dimy-starty
	tstrtx=startx-1
	tstrty=starty-1
	xh=scalex*sin(phi)
	yh=scaley*cos(phi)
	xv=scalex*cos(phi)*cos(theta)
	yv=scaley*sin(phi)*cos(theta)
	if(xh.eq.0.)xh=.000001
	if(yh.eq.0.)yh=.000001
	if(xv.eq.0.)xv=.000001
	if(yv.eq.0.)yv=.000001
	dh=tdimy/yh
 	if(theta.le..000001) theta=.000001
	dv=amax/(scalez*sin(theta))
	if(abs(yh).lt..01) dh=500.*amax
	if(theta.lt..01) dv=500.*amax
	dely=abs(xh*tdimy/yh/tdimx)
	delx=abs(yh*tdimx/xh/tdimy)
	scanx=tdimx/abs(xv).lt.tdimy/abs(yv)
	sto3=dv*xv/tdimx
	sto4=dv*yv/tdimy
	sto7=dh*xh/tdimx
	sto9=dh*yh/tdimy
c*****FIND VERTICAL SCALE DATA 
c*****THE MINIMUM VALUE FOUND IS OFTEN HIDDEN SO SCALING
c*****IN THE VERTICAL DIRECTION IS NOT ALWAYS CORRECT,ESPESCIALLY
c*****WITH LARGE THETA
      bmin=1.0e38
      bmax=-1.0e38
      if(resp.ne.'n'.and.resp.ne.'N') then
        do 4 i = 1,nr
	 st8 = xv*(i-1)/tdimx
        do 3 j = 1,nc
        if(a(j,i).lt.sdval) go to 3
	 dum1 = a(j,i) - (dv*(yv*(j-1)/tdimy + st8))
        if(dum1.lt.bmin) bmin=dum1
        if(dum1.gt.bmax) bmax=dum1
3       continue
4       continue
      else
        do 6 i=startx,dimx,skipx
	 st8 = xv*(i-startx)/tdimx
        do 5 j=starty,dimy,skipy
        if(a(j,i).lt.sdval) go to 5
	 dum1 = a(j,i) - (dv*(yv*(j-starty)/tdimy + st8))
        if(dum1.lt.bmin) bmin=dum1
        if(dum1.gt.bmax) bmax=dum1
5       continue
6       continue
      endif
c*****FIND HORIZONTAL SCALE FACTOR
	xhdh = xh * dh
      dxp(1) = amin1(0.,-xhdh)
      dxp(2) = amax1(tdimy,tdimy-xhdh)
      dyp(1) = bmin
      dyp(2) = bmax
	ctype = 'c'
	call pcolor(14)
        dcval = (amax-amin)/13.
      call scale(dxp,dyp,xp,yp,4,ier)
      if(ier.ne.0) then
        print *,'error in subroutine scale'
        pause 'press ENTER to return to menu'
        stop
      endif
      call neatl
	st81 = yv/tdimy
	xvtdim = xv/tdimx
	tstx15 = tstrtx + 1.5
	tsty15 = tstrty + 1.5
       if(scanx) then
         if(-yh.gt.0.) then
           dm=-1.
           n1=startx
         else
           dm=1.
           n2=dimx
         endif
         if(yv.ge.0.) then
           dj=dely
         else
           dj = -dely
         endif
       else
         if(-xh.gt.0.) then
           dj=-1.
           n1=starty
         else
           dj=1.
           n2=dimy
         endif
         if(xv.ge.0.) then
           dm=delx
         else
           dm = -delx
         endif
       endif
       dmsto3 = dm*sto3
	djsto4 = dj*sto4

	do 220 m=startx,dimx,skipx
	tm=m-startx
	sto2=tm+tstx15
	sto5=tm/tdimx*xhdh
	st8 = xvtdim*tm
	see=.true.
	saw=.false.
	j1=0
	do 210 jc=starty,dimy
	j=jc
c      point=0.
c	if(.not.array) go to 20
c	point=a(m,j)
	point=a(j,m)
c	go to 30
c20     call extsch(a,m,j,point)
c30     continue
	tj=j-starty
	sto1=tj+tsty15
c	sto8=dv*(xv*tm/tdimx+yv*tj/tdimy)
	sto8=dv*(st8+st81*tj)
	if(.not.opaque) go to 190
       if(scanx) then
         if(-yh.gt.0.) then
           n2=m
         else
           n1=m
         endif
       else
         if(-xh.gt.0.) then
           n2=j
         else
           n1=j
         endif
       endif
90     ti=0.
	do 120 i=n1,n2
	ti=ti+1.
	jip=dj*ti+sto1
	mip=dm*ti+sto2
	if(jip.gt.dimy.or.jip.lt.starty.or.mip.gt.dimx.or.mip.lt.startx)
     & go to 120
c	pointp=0.
c	if(.not.array) go to 100
c	pointp=a(mip,jip)
	pointp=a(jip,mip)
c	go to 110
c100    call extsch(a,mip,jip,pointp)
c110    continue
	see=pointp-abs(dmsto3*ti)-abs(djsto4*ti).le.point
	if(.not.see) go to 130
120    continue
130    if(see) go to 190
	if(.not.saw) go to 210
140    j=j-1
150    if(j1.lt.2) j1=0
	if(j1.eq.0) go to 210
	j2=j-j1+1
	if(block) j2=j-j1/2+1
	l=1
	do 160 jp=j2,j
	tjp=jp-starty
	sh(l)=tjp-sto5
	if(.not.block) go to 160
	l=l+1
	sh(l)=sh(l-1)+1.
160    l=l+1
       do 165 k = 1,j1-1
       if(sv(k).lt.sdval.or.sv(k+1).lt.sdval) go to 165
	call line(sh(k),sv(k),2,0,icolor(k))
165    continue
c 	call line(sh,sv,j1,0,0)
167    if(.not.block) go to 180
	do 170 i=1,j1
	sv(i)=sv(i)-sto3
170    sh(i)=sh(i)-sto7
	do 175 k=1,j1-1
       if(sv(k).lt.sdval.or.sv(k+1).lt.sdval) go to 175
	call line(sh(k),sv(k),2,0,icolor(k))
175    continue
c 	call line(sh,sv,j1,0,0)
180    j1=0
	if(jc.eq.dimy) go to 220
	saw=.false.
	go to 210
c190    if(point.lt.amin.or.point.gt.amax) go to 140
190    j1=j1+1
	saw=see
	sv(j1)=point-sto8
c       print *,'j1,jc',j1,jc,sv(j1),point,sto8
	if(point.gt.sdval) 
     & icolor(j1) = (ifix((point-amin)/dcval) + 1) * 100 
	if(.not.block) go to 200
	j1=j1+1
	sv(j1)=sv(j1-1)-sto4
	if(point.gt.sdval) 
     & icolor(j1) =  icolor(j1-1)
200    if(jc.eq.dimy) go to 150
210    continue
220    continue
	st81 = xv/tdimx
	yvtdim = yv/tdimy
	do 430 j=starty,dimy,skipy
	tj=j-starty
	sto1=tj+tsty15
	st8 = yvtdim*tj
	see=.true.
	saw=.false.
	m1=0
	do 420 mc=startx,dimx
	m=mc
c	point=0.
c	if(.not.array) go to 230
c	point=a(m,j)
	point=a(j,m)
c	go to 240
c230    call extsch(a,m,j,point)
c240    continue
	tm=m-startx
	sto2=tm+tstx15
c	sto8=dv*(xv*tm/tdimx+yv*tj/tdimy)
	sto8=dv*(st81*tm+st8)
	if(.not.opaque) go to 400
       if(scanx) then
         if(-yh.gt.0.) then
           n2=m
         else
           n1=m
         endif
       else
         if(-xh.gt.0.) then
           n2=j
         else
           n1=j
         endif
       endif
300    ti=0.
	do 330 i=n1,n2
	ti=ti+1.
	jip=dj*ti+sto1
	mip=dm*ti+sto2
	if(jip.gt.dimy.or.jip.lt.starty.or.mip.gt.dimx.or.mip.lt.startx)
     & go to 330
c      pointp=0.
c	if(.not.array) go to 310
c	pointp=a(mip,jip)
	pointp=a(jip,mip)
c	go to 320
c310    call extsch(a,mip,jip,pointp)
c320    continue
	see = pointp-abs(dmsto3*ti)-abs(djsto4*ti).le.point
	if(.not.see) go to 340
330    continue
340    if(see) go to 400
	if(.not.saw) go to 420
350    m=m-1
360    if(m1.lt.2) m1=0
	if(m1.eq.0) go to 420
	m2=m-m1+1
	if(block) m2=m-m1/2+1
	l=1
	do 370 mp=m2,m
	tmp=mp-startx
	sh(l)=tj-tmp*sto7
	if(.not.block) go to 370
	l=l+1
	sh(l)=sh(l-1)-sto7
370    l=l+1
       do 375 k = 1,m1-1
       if(sv(k).lt.sdval.or.sv(k+1).lt.sdval) go to 375
	call line(sh(k),sv(k),2,0,icolor(k))
375    continue
c 	call line(sh,sv,m1,0,0)
377    if(.not.block) go to 390
	do 380 i=1,m1
	sv(i)=sv(i)-sto4
380    sh(i)=sh(i)+sto9
       do 385 k = 1,m1-1
       if(sv(k).lt.sdval.or.sv(k+1).lt.sdval) go to 385
	call line(sh(k),sv(k),2,0,icolor(k))
385    continue
c 	call line(sh,sv,m1,0,0)
390    m1=0
	if(mc.eq.dimx) go to 430
	saw=.false.
	go to 420
c400    if(point.lt.amin.or.point.gt.amax) go to 350
400    m1=m1+1
	sv(m1)=point-sto8
	if(point.gt.sdval) 
     & icolor(m1) = (ifix((point-amin)/dcval) + 1) * 100 
	saw=see
	if(.not.block) go to 410
	m1=m1+1
	sv(m1)=sv(m1-1)-sto3
	if(point.gt.sdval) 
     & icolor(m1) = icolor(m1-1) 
410    if(mc.eq.dimx) go to 360
420    continue
430    continue
	return
	end
c	subroutine extsch(data,x,y,z)
c       integer x,y
c	This is a dummy routine to be written if data are not
c	in the form of a two dimensional array
c	data = list of x y z points
c	x = x coordinate ranging from startx to dimx
c	y = y coordinate ranging from starty to dimy
c	z = value at x and y
c       return
c       end
