	  subroutine anacon(xarr,nx,ny,xlev,numlev,itvord,
     &		 xl,yl,xsc,ysc,cut)
c
c   xarr  is data array to contour (real*4) 
c   nx,ny  is size of xarr
c   xlev   is a table of contour levels real*4
c   numlev are the number of levels to contour
c   itvord  is the desired flip of the array - see inccon below for details
c   xsc,ysc are the coords. of the upper rhs
c   xl,yl are the coordinates of the lower left corner for the contour
c	     
	  real*4 xarr(nx*ny),aa(4),bb(4)
	  real*4 xlev(numlev),xa(4),ya(4)
	  integer*4 ltvcol(8)
	  real*4 cell(5),xmax,xmin

c	type *,'numlev, nx,ny =',numlev, nx,ny
c	do i=1,numlev
c	type *,i,xlev(i)
c	enddo
c--	compute internal values for spacing
c	type *,'xl, yl, xsc, ysc', xl, yl, xsc, ysc
	xsc=xsc-xl
	ysc=ysc-yl
	if(mod(itvord,2).eq.0)then !match tektronix x,y with data x,y
	ysc=ysc/ny
	xsc=xsc/nx
	else
	ysc=ysc/nx
	xsc=xsc/ny
	endif
c must reset ix,iy to center the contour on the raster
	  xz=xl+xsc/2.
	  yz=yl+ysc/2.
c	type *,'xz, yz, xsc, ysc', xz, yz, xsc, ysc
c
c set up increments and limits for array processing according to itvord
c
	  call inccon(itvord,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4)
	  minc=ny-1
	  if(ninc.eq.minc)minc=nx-1
	  do irow=1,minc	!row counter for position on grinnell
	     do icol=1,ninc   !col counter for position of grinnell
		  cell(1)=xarr(ibase)
		  cell(2)=xarr(ibase+n2)
		  cell(3)=xarr(ibase+n3)
		  cell(4)=xarr(ibase+n4)
		  cell(5)=cell(1)
c set ibase for next cell
		  ibase=ibase+ibindx
		  xmax=cell(1)
		  xmin=cell(1)
		  imax=1
		  do j=2,4				    !find max,min of cell
		     xmax=amax1(xmax,cell(j))
		     xmin=amin1(xmin,cell(j))
		     if(cell(j).eq.xmax)imax=j	  !keep track of max
		  enddo
c
c  now check each cell for contours
c
		  do k=1,numlev
		     if(xlev(k).le.xmax.and.xlev(k).ge.xmin)then
c  level goes through cell 
			  ndsh=1
			  if(xlev(k).lt.cut)ndsh=2
			  ns=0			     !number of sides 
			  do l=1,4
			     if(xlev(k).ge.cell(l).and.xlev(k).lt.
     &			  cell(l+1).or.xlev(k).lt.cell(l).and.
     &			  xlev(k).ge.cell(l+1))then
c   
c  calculate endpoints of line segment on the side of the cell
c
			 ns=ns+1 
			if(l.eq.1)then	 !top
				xa(ns)=float(icol-1)+(xlev(k)-cell(l))
     &				  /(cell(l+1)-cell(l))
				 ya(ns)=float(irow-1)
			 else if(l.eq.2)then  !right side
				xa(ns)=float(icol)
				ya(ns)=float(irow-1)+(xlev(k)-cell(l))
     &				 /(cell(l+1)-cell(l))
			else if (l.eq.3)then !bottom side
				xa(ns)=float(icol-1)+(xlev(k)-
     &				  cell(l+1))/(cell(l)-cell(l+1))
				ya(ns)=float(irow)
			 else		     !left side 
				xa(ns)=float(icol-1)
				ya(ns)=float(irow-1)+(xlev(k)-
     &				   cell(l+1))/(cell(l)-cell(l+1))
				endif
			     endif
			  enddo
c
			  if(ns.gt.2)then
			     ntim=2
			     xmean=cellavg(cell,4)
c  if the present level is greater than the mean, then the contour
c  is drawn between the max point and the center. (imax gives max
c  point)
c  first calculate all screen coordinates
			     in=imax
			     do ii=1,4
				  aa(ii)=xz+xa(in)*xsc
				  bb(ii)=yz+ya(in)*ysc
				  in=in-1
				  if(in.lt.1)in=4
			     enddo
			     if(xlev(k).le.xmean)then
c switch 2 and 4 so that point 1 connects to point 2 and point 3 to 4.
				  tmp=aa(4)
				  aa(4)=aa(2)
				  aa(2)=tmp
				  tmp=bb(4)
				  bb(4)=bb(2)
				  bb(2)=tmp
			     endif
			  else if(ns.eq.2)then
			     ntim=1
			     do ii=1,2
				  aa(ii)=xz+xa(ii)*xsc
				  bb(ii)=yz+ya(ii)*ysc
			     enddo
			  endif				!ns.gt.2 
c now draw the contours
				if (ns.ge.2) then
				  call tkdash(aa,bb,ndsh,ntim)
				endif
		     endif				   !xlev(k).gt.... 
		  enddo					!k=1,numlev 
	     enddo					   !icol=1,nx-1
	     ibase=ibase+incr     !reset for next row or col
	  enddo						!irow=1,ny-1
999     return
	  end
c
	  function cellavg(cell,nav)
	  real*4 cell(*)
	  cellavg=0.0
	  do ii=1,nav
	     cellavg=cellavg+cell(ii)
	  enddo
	  cellavg=cellavg/float(nav)
999     return
	  end
c
	  subroutine inccon(itv,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4)
c get limits and increments for processing the raw array
c inputs:
c itv is the order to process :
c	    0  left to right, bottom to top
c	    1  bottom to top, left to right
c	    2  left to right, top to bottom
c	    3  top to bottom, left to right
c	    4  right to left, top to bottom
c	    5  top to bottom, right to left
c	    6  right to left, bottom to top
c	    7  bottom to top, right to left
c nx,ny  size of raw array (columns,rows)
c 
c outputs:
c ibase base index to begin processing
c ibindx index used to increment base address during processing
c incr used at end of a row or column to reset base address
c n2,n3,n4 offsets in base address for 4 point cell
c
c note that even though itv may indicate array is to be processed
c vertically in a certain direction (either top to bottom or bottom to
c top), the outputs are such that the array is processed in the
c opposite direction.  this is because the contours are actually drawn
c from the bottom of the screen to the top.
c
c--	7/7/86        orders 1 and 5 seem to be reversed, patch to fix
c	the comments and logic in parts of this seem wrong, should be cleaned
c	up someday, beware !
	if(mod(itv,2).eq.0)then    !process horizontal before vertical
	   ibindx=1
	   if(itv.gt.3)ibindx=-1   !right to left
	   incr=1
	   if(itv.eq.2.or.itv.eq.6)incr=2*nx-1
	   if(itv.eq.2.or.itv.eq.4)incr=-incr
	   ninc=nx-1		   !process to end of row before reset
	else			     !process vertical before horizontal
	   ibindx=nx    !processing top to bottom
	   if(itv.eq.5.or.itv.eq.7)ibindx=-nx    !nope, bottom to top
	   incr=(ny-1)*nx+1  !reset a col to begin processing next col
	   if(itv.eq.1.or.itv.eq.5)incr=incr-2
	   if(itv.eq.3.or.itv.eq.1)incr=-incr
	   ninc=ny-1   !process to end of column before reset
	endif
c now get initial base address
	if(itv.eq.2.or.itv.eq.7)then
	   ibase=(ny-1)*nx+1	 !start at lower left
	else if(itv.eq.4.or.itv.eq.5)then
	   ibase=nx*ny		 !lower right
	else if(itv.eq.0.or.itv.eq.1)then
	   ibase=1		     !upper left
	else
	   ibase=nx		    !upper right
	endif
c now need to fill the cell properly - figure offsets from base address
	if(itv.eq.0)then
	   n2=1
	   n3=nx+1
	   n4=nx
	else if(itv.eq.3)then
	   n2=nx
	   n3=nx-1
	   n4=-1
	else if(itv.eq.2)then
	   n2=1
	   n3=-nx+1
	   n4=-nx
	else if(itv.eq.5)then
	   n2=-nx
	   n3=-nx-1
	   n4=-1
	else if(itv.eq.4)then
	   n2=-1
	   n3=-nx-1
	   n4=-nx
	else if(itv.eq.7)then
	   n2=-nx
	   n3=-nx+1
	   n4=1
	else if(itv.eq.6)then
	   n2=-1
	   n3=nx-1
	   n4=nx
	else if(itv.eq.1)then
	   n2=nx
	   n3=nx+1
	   n4=1
	endif
c	type *,'variables in inccon:'
c	type 1001,itv,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4
c1001  format(1x,'itv,nx,ny,ibase,ibindx,incr,ninc,n2,n3,n4',/,1x,10i6)
	return
	end
