c
c*******************************************************************************
c*******************************************************************************
c
c
c
      subroutine contour_tri (z, nx, ny, parms, link, cont_disp)
c
c
      real*4	z(*)
      integer*4	nx, ny
      real*4	parms(*)
      byte	link(3,*)
      external	cont_disp
c
c
c*******************************************************************************
c*******************************************************************************
c
c	Contour the data field Z using triangular subdivisions of the
c	rectangular grid.
c
c	For each of the NLEV contour levels, scan the field for an intersection
c	between the contour level ZLEV and one of the grid lines.  Such
c	an intersection, or "hit", will occur if one point at the end of
c	a grid line segment is greater than or equal to ZLEV and the other
c	is less than ZLEV, or one point is less than or equal to ZLEV and the
c	other is greater than ZLEV.
c
c	Sometimes, data fields contain very steep or step gradients that
c	would normally result in several contours occupying the same or very
c	nearly the same space.  For each user-specified contouring level,
c	the user-specified ZHIGH value sets the limit on the maximum Z value
c	for the areas through which contour lines for that level may pass.
c	This capability may be over-ridden by specifying the ZLIMIT value
c	as 1.0e+30, thereby permitting very steep gradients to be contoured
c	in the usual manner.
c
c
c	CONTOUR does no plotting or displaying itself.  Whenever a critical
c	step in the contouring process is encountered, CONTOUR calls the
c	user-supplied subroutine whose name was passed as CON_DISP.  All
c	calls have the form:
c
c		CALL CON_DISP (KODE, X, Y)
c
c	where the three arguments have the following definitions:
c
c	    KODE    X    Y	Remarks
c
c	    0	    XPT  YPT	The point (XPT,YPT) is a point that
c				continues the current contour line.
c
c	    90	    NLEV 0.0	Signals the start of the contouring process.
c				NLEV is the number of levels that will be
c				contoured.
c
c	    91	    ILEV ZLEV	Signals the start of contouring for the
c				level ILEV whose value is ZLEV.
c
c	    92      XPT  YPT	The point (XPT,YPT) is the start of a
c				contour line at the current level.
c
c	    95      XSTRT YSTRT	The current contour line has hit a dead end.
c				Contouring will resume at the start of the
c				current contour line.  This signal is provided
c				so that CON_DISP can reverse the order of the
c				points on the current contour line, thereby
c				making a single line out of what would
c				otherwise be two fragments.
c
c	    96      XSTRT YSTRT	Signals that the current contour line has
c				closed on the starting point.
c
c	    97      0.0  0.0	Signals that the current contour line has
c				ended without closing on the starting point.
c				The previous values for XPT and YPT represent
c				the end point of this line.
c
c	    98	    ILEV ZLEV	Signals the end of contouring for the level
c				ILEV whose value is ZLEV.
c
c	    99	    NLEV 0.0	Signals the end of the contouring process.
c				NLEV is the number of levels that were
c				contoured.
c
c	Values for KODE other than the ones specified above may be used
c	by CON_DISP for user-defined purposes (for instance, initializing
c	the display, changing colors or line types, etc.).
c
c
c	Inputs:
c
c	    Z		real*4		Data field to be contoured.
c					Dimensioned NX * NY
c
c	    NX		integer*4	Dimensions of the data field Z.
c	    NY
c
c	    PARMS	real*4		An array used to pass contouring
c					parameters.
c
c	    LINK	byte		An array used as work space.  Must 
c					be dimensioned at least:
c					    3 * NX * NY
c
c	    CON_DISP	external	The name of the subroutine that
c					will be called to dispose of the
c					contour information.
c
c
c
c	Outputs:
c
c	    None
c
c
c
c	History:
c
c	    01 May 88	Phil McDonald, NOAA/ERL/FSL	Operational version
c
c	    11 Jun 91	Phil McDonald, NOAA/ERL/FSL	Substantially revamped.
c
c
c*******************************************************************************
c*******************************************************************************
c
c
c

      include		'contour_tri.inc'



      byte		lnktmp(3)

      logical		do_border, dead_end



      if (nx .lt. 2) return
      if (ny .lt. 2) return
c
c*******************************************************************************
c
c	Initialize things
c
c*******************************************************************************
c
      nxdim = nx

      npt   = nx * ny

      ip    = 0
c-------------------------------------------------------------------------------
c
c	Unpack the triangulation flag passed through PARMS.
c
c	If the triangulation flag is positive, the grid squares will
c	be divided into triangles along lines running in the positive
c	direction (i.e. from point (1,1) to point (2,1)).  If the
c	triangulation flag is negative, the grid squares will be
c	divided into triangles along lines running in the negative
c	direction (i.e. from point (1,2) to point (2,1)).
c
c-------------------------------------------------------------------------------
c
      ip = ip + 1
      if (parms(ip) .gt. 0.0) then
          idp(1) = 1
          idp(2) = nx + 1
          idp(3) = nx
      else
          idp(1) = -1
          idp(2) = nx - 1
          idp(3) = nx
      end if


      do ipt = 1, npt
          link(1,ipt) = 0
          link(2,ipt) = 0
          link(3,ipt) = 0
      end do
      do ipt = npt-nx+1, npt
          link(2,ipt) = -2
          link(3,ipt) = -2
      end do
      if (idp(1) .gt. 0) then
          do ipt = nx, npt, nx
              link(1,ipt) = -2
              link(2,ipt) = -2
          end do
      else
          do ipt = 1, npt, nx
              link(1,ipt) = -2
              link(2,ipt) = -2
          end do
      end if
c
c-------------------------------------------------------------------------------
c
c	Unpack the ZLIMIT information passed through PARMS.
c
c		ZLIMIT		The upper limit for all contours.  A 
c				contour line will not enter an area
c				with Z values greater than this value.
c
c				If ZLIMIT = 1.0e+30 and NLEV > 0 then
c				the user must supply a value for ZHIGH
c				for each of the NLEV contouring levels.
c				Otherwise, the value of ZLIMIT will be
c				used as for the value of ZHIGH for all
c				contouring levels.
c
c
c-------------------------------------------------------------------------------
c
      ip     = ip + 1
      zlimit = parms(ip)
c
c-------------------------------------------------------------------------------
c
c	Unpack the contouring information passed through PARMS.
c
c		NLEV > 0	There will be NLEV contouring levels, one
c				at each of the NLEV values of ZLEV specified
c				by the user.  Starting with the PARMS
c				element after the one containing NLEV, there
c				must be NLEV ZLEV values.  If ZLIMIT is not
c				equal to 1.0e+30, then these values must be
c				followed by NLEV ZHIGH values.
c
c		NLEV = 0	The next element of PARMS contains ZINT, the
c				contouring interval.  The entire range of
c				Z values will be contoured, starting with
c				the lowest possible integer multiple of ZINT,
c				and ending with the highest possible integer
c				multiple of ZINT.
c
c		NLEV < 0	The next element of PARMS contains the lowest
c				contour level and the element after that one
c				contains the highest.  This range will be
c				divided into |NLEV| evenly spaced contour
c				intervals.
c
c-------------------------------------------------------------------------------
c
      ip     = ip + 1
      nlev   = parms(ip)

      zmin = 1.0e+30
      zmax = zmin
      zint = 0.0
      if (nlev .gt. 0) then
          zmin = -zmax
      else if (nlev .lt. 0) then
          ip   = ip + 1
          zmin = parms(ip)
          ip   = ip + 1
          zmax = parms(ip)
          if (zmax .gt. zmin) then
              zint = (zmax - zmin) / float ((-nlev) - 1)
          else
              zint = 1.0
          end if
      else if (nlev .eq. 0) then
          ip   = ip + 1
          zint = parms(ip)
          if (zint .eq. 0.0) return
          zmin = z(1)
          zmax = z(1)
          do ipt = 2, npt
              if (z(ipt) .lt. zmin) then
                  zmin = z(ipt)
              else if (z(ipt) .gt. zmax) then
                  zmax = z(ipt)
              end if
          end do
          zmin = float (nint (zmin / zint)) * zint
          zmax = float (nint (zmax / zint)) * zint
      end if
c
c*******************************************************************************
c
c	Signal the start of contouring.
c
c*******************************************************************************
c
      x = nlev
      y = 0.0
      call cont_disp (90, x, y)
c
c*******************************************************************************
c
c	Process each contour level
c
c*******************************************************************************
c
      ilev = 0
      zlev = zmin - zint
      do while (zlev .lt. zmax)

          ilev  = ilev + 1
c
c-------------------------------------------------------------------------------
c
c	For user-specified contour levels, unpack the parameters for this
c	level from PARMS:
c
c		ZLEV		The value for this contour
c
c
c		ZHIGH		The upper limit for this contour.  A 
c				contour line will not enter an area
c				with Z values greater than this value.
c
c-------------------------------------------------------------------------------
c
          zhigh = zlimit
          if (zint .eq. 0.0) then
              ip   = ip + 1
              zlev = parms(ip)
              if (ilev .eq. nlev) zmax = zlev
              if (zlimit .ne. 1.0e+30) zhigh = parms(ip + nlev)
          else
              zlev = zmin + (zint * float (ilev - 1))
          end if
c
c-------------------------------------------------------------------------------
c
c	Clear the flags used to prevent grid points from being re-used
c
c-------------------------------------------------------------------------------
c
          do ipt = 1, npt
              if (link(1,ipt) .gt. 0) link(1,ipt) = 0
              if (link(2,ipt) .gt. 0) link(2,ipt) = 0
              if (link(3,ipt) .gt. 0) link(3,ipt) = 0
          end do
c
c-------------------------------------------------------------------------------
c
c	Signal the start of this contour level
c
c-------------------------------------------------------------------------------
c
          x = ilev
          y = zlev
          call cont_disp (91, x, y)
c
c*******************************************************************************
c
c	Find the contour lines for this contour level
c
c*******************************************************************************
c
c	Scan the data field, one grid interval at a time, looking for
c	a "hit" -- an intersection between the contour line and a grid
c	line.  Such and intersection exists if a Z value is greater than
c	or equal to ZLEV and an adjacent one is is less, or a Z value is
c	less than or equal to ZLEV and an adjacent one is greater.
c
c	When the whole field has been scanned, go to the next contour level.
c
c-------------------------------------------------------------------------------
c
          ipt = 0
   10     if (ipt .ge. npt) go to 100
          ipt = ipt + 1

          ileg = -1
   11     if (ileg .ge. 3) go to 10
          ileg = ileg + 2

          if (link(ileg,ipt) .ne. 0) go to 11

          jpt = ipt + idp(ileg)
c
c...............................................................................
c
c	Check for a hit
c
c...............................................................................
c
          if (z(jpt) .eq. zlev) go to 11
          if (z(ipt) .lt. zlev) then
              if (z(jpt) .lt. zlev) then
                  link(ileg,ipt) = -1
                  go to 11
              end if
          else
              if (z(jpt) .ge. zlev) then
                  link(ileg,ipt) = 1
                  go to 11
              end if
          end if
          link(ileg,ipt) = 1
c
c-------------------------------------------------------------------------------
c
c	A hit has been found between points IPT and JPT.
c	Save the current scan position and these starting positions.
c
c-------------------------------------------------------------------------------
c
          ipstrt = ipt
          ilstrt = ileg
          ipscan = ipt
          ilscan = ileg
          idir   = 2
c
c-------------------------------------------------------------------------------
c
c	Start a contour line.  Signal the start of the contour line
c	and initialize pointers to the high and low points.
c
c-------------------------------------------------------------------------------
c
          dead_end = .FALSE.
          nhit     = 0

          call cont_find (ipt, jpt, z, zlev, xhit, yhit, nhit)

          xstart = xhit
          ystart = yhit

          call cont_disp (92, xstart, ystart)
c
c-------------------------------------------------------------------------------
c
c	Trace the contour through the field.
c
c	Find the available unused points linked to both point IPT and to
c	point JPT.  Point JPT is in the direction IDIR from point IPT.
c	The number of available points may be 0, 1, or 2.
c
c-------------------------------------------------------------------------------
c
   20     kpt  = ipt
          kleg = ileg
          kdir = idir

          call cont_next (ipt, ileg, idir, z, zlev, link)
c
c-------------------------------------------------------------------------------
c
c	If a hit was found, pass the new point on the contour line to 
c	CON_DISP, and then continue tracing the contour.
c
c-------------------------------------------------------------------------------
c
   25    if (ipt .gt. 0) then

              jpt = ipt + idp(ileg)

              call cont_find (ipt, jpt, z, zlev, xhit, yhit, nhit)

              call cont_disp (0, xhit, yhit)

              go to 20

          end if
c
c-------------------------------------------------------------------------------
c
c	If a hit was not found, then the contour line has hit a dead end.
c	If this is the second dead end, then the line has hit dead ends
c	in both directions and cannot be closed.  Presumably these dead
c	ends are at the border.  Terminate the current contour line and
c	indicate that the line has failed to close.
c
c-------------------------------------------------------------------------------
c
          kode = 97
          xhit = 0.0
          yhit = 0.0

          if (dead_end) go to 90
c
c...............................................................................
c
c	Since this is the first dead end, temporarily reopen the link at
c	the starting point and see if a connection can be made.  If it
c	can, then the line has closed on the starting point.  Terminate
c	the current contour line and indicate that the line has closed.
c
c...............................................................................
c
          dead_end = .TRUE.

          if (nhit .gt. 1) then

              ipt            = kpt
              ileg           = kleg
              idir           = kdir

              call cont_copy (link(1,ipstrt), lnktmp)

              call cont_next (ipt, ileg, idir, z, zlev, link)

              call cont_copy (lnktmp, link(1,ipstrt))

              if (ipt .gt. 0) then
                  kode = 96
                  xhit = xstart
                  yhit = ystart
                  go to 90
              end if

          end if
c
c...............................................................................
c
c	The current contour line cannot be extended in the current direction.
c	Presumably, this is because the current end point is on the border.
c	Return to the starting point and see if the line can be extended
c	in the opposite direction.  If this is not possible, presumably
c	the starting point is on the border.  At any rate, the line has
c	failed to close.  However, if the line can be extended in to 
c	opposite direction, do so and continue to trace the contour.
c
c...............................................................................
c
          ipt            = ipstrt
          ileg           = ilstrt
          idir           = 1

          call cont_copy (link(1,ipstrt), lnktmp)

          call cont_next (ipt, ileg, idir, z, zlev, link)

          call cont_copy (lnktmp, link(1,ipstrt))

          if (ipt .gt. 0) then
              link(ileg,ipt) = 1
              call cont_disp (95, xstart, ystart)
              go to 25
          end if
c
c-------------------------------------------------------------------------------
c
c	Signal the end of the contour line and then resume scanning the
c	data field.
c
c-------------------------------------------------------------------------------
c
   90     call cont_disp (kode, xstart, ystart)

          ipt  = ipscan
          ileg = ilscan

          go to 11
c
c*******************************************************************************
c
c	The entire data field has been scanned.  Signal the end of this 
c	contour level
c
c*******************************************************************************
c
  100     x = ilev
          y = zlev
          call cont_disp (98, x, y)
c
      end do
c
c*******************************************************************************
c
c	All contour levels have been processed.  Signal the end of contouring.
c
c*******************************************************************************
c
      x = nlev
      y = 0.0
      call cont_disp (99, x, y)
c
c
c
      return
      end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
      subroutine cont_next (ipt, ileg, idir, z, zlev, link)
c
c
      integer*4	ipt, ileg, idir
      real*4	z(*), zlev
      byte	link(3,*)
c
c
c*******************************************************************************
c*******************************************************************************
c
c	Test for an intersection between a contour line and a grid
c
c*******************************************************************************
c*******************************************************************************
c
c



      include		'contour_tri.inc'



      logical	same_side



      jpt  = ipt
      jleg = ileg
      jdir = idir

      ipt  = 0
c
c*******************************************************************************
c
c	Find the next point to test.  If this point is outside the grid,
c	ignore it.
c
c*******************************************************************************
c
      kpt  = jpt
      kleg = jleg
      if (jdir .eq. 1) then
          kleg = kleg - 1
          if (kleg .lt. 1) then
              kleg = 3
              kpt  = kpt - idp(3)
          end if
      else
          kleg = kleg + 1
          if (kleg .gt. 3) then
              kleg = 1
              kpt  = kpt - idp(1)
          end if
      end if

      if (kpt .le. 0) return
      if (link(kleg,kpt) .eq. -2) return

      if (kpt .eq. jpt) kpt = kpt + idp(kleg)
c
c*******************************************************************************
c
c	Now find the leg to follow.
c
c*******************************************************************************
c
      if (z(jpt) .lt. zlev) then
          if (z(kpt) .ge. zlev) jdir = 3 - jdir
      else
          if (z(kpt) .lt. zlev) jdir = 3 - jdir
      end if
      same_side = jdir .eq. idir

      jleg = jleg + jdir
      if (jleg .gt. 3) then
          jleg = jleg - 3
          jdir = 3 - jdir
      end if
      if (idir .eq. 2) jdir = 3 - jdir

      if (same_side) then
          if (jdir .eq. idir) then
              jpt = kpt
          else
              if (idir .eq. 1) then
                  jpt = jpt + idp(3)
              else
                  jpt = jpt + idp(1)
              end if
          end if
      else
          if (jdir .ne. idir) jpt = kpt
      end if


      if (jpt .gt. 0) then
          if (link(jleg,jpt) .eq. 0) then
              ipt            = jpt
              ileg           = jleg
              idir           = jdir
              link(ileg,ipt) = 1
          end if
      end if
c
c
c
      return
      end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
      subroutine cont_copy (link1, link2)
c
c
      byte	link1(3), link2(3)
c
c
c*******************************************************************************
c*******************************************************************************
c
c	Save the values of LINK1 in LINK2 and then reopen the links
c	of LINK1.
c
c*******************************************************************************
c*******************************************************************************
c
c



      link2(1) = link1(1)
      link2(2) = link1(2)
      link2(3) = link1(3)

      if (link1(1) .ne. -2) link1(1) = 0
      if (link1(2) .ne. -2) link1(2) = 0
      if (link1(3) .ne. -2) link1(3) = 0
c
c
c
      return
      end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
      subroutine cont_find (ipt, jpt, z, zlev, xhit, yhit, nhit)
c
c
      integer*4	ipt, jpt
      real*4	z(*), zlev, xhit, yhit
      integer*4	nhit
c
c
c*******************************************************************************
c*******************************************************************************
c
c	Test for an intersection between a contour line and a grid
c
c*******************************************************************************
c*******************************************************************************
c
c



      include		'contour_tri.inc'



      call cont_i2xy (ipt, ix, iy)
      call cont_i2xy (jpt, jx, jy)

      zz   = zlev
      dz   = z(jpt) - z(ipt)
      if (zz .eq. z(ipt)) then
          zz = zz + (dz * 0.1)
      else if (zz .eq. z(jpt)) then
          zz = zz - (dz * 0.1)
      end if
      zrat = (zz - z(ipt)) / (z(jpt) - z(ipt))

      dx   = zrat * float (jx - ix)
      dy   = zrat * float (jy - iy)
      xhit = float (ix) + dx
      yhit = float (iy) + dy

      nhit = nhit + 1
c
c
c
      return
      end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
      subroutine cont_i2xy (ipt, ix, iy)
c
c
      integer*4	ipt, ix, iy
c
c
c*******************************************************************************
c*******************************************************************************
c
c	Convert a point index IPT into grid IX and IY.
c
c*******************************************************************************
c*******************************************************************************
c
c



      include		'contour_tri.inc'



      iy = (ipt - 1) / nxdim
      ix = ipt - (iy * nxdim)
      iy = iy + 1
c
c
c
      return
      end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
      subroutine cont_disp (kode, x, y)
c
c
      integer*4	kode
      real*4	x, y
c
c
c*******************************************************************************
c*******************************************************************************
c
c	A sample subroutine to dispose of contour line information.
c
c	The contouring subroutine CONTOUR does no plotting or displaying
c	itself.  Whenever a critical step in the contouring process is
c	encountered, CONTOUR calls the user-supplied subrutine whose name
c	was passed as CON_DISP.  All calls have the form:
c
c		CALL CON_DISP (KODE, X, Y)
c
c	where the three arguments have the following definitions:
c
c	    KODE    X    Y	Remarks
c
c	    0	    XPT  YPT	The point (XPT,YPT) is a point that
c				continues the current contour line.
c
c	    90	    NLEV 0.0	Signals the start of the contouring process.
c				NLEV is the number of levels that will be
c				contoured.
c
c	    91	    ILEV ZLEV	Signals the start of contouring for the
c				level ILEV whose value is ZLEV.
c
c	    92      XPT  YPT	The point (XPT,YPT) is the start of a
c				contour line at the current level.
c
c	    95      XSTRT YSTRT	The current contour line has hit a dead end.
c				Contouring will resume at the start of the
c				current contour line.  This signal is provided
c				so that CON_DISP can reverse the order of the
c				points on the current contour line, thereby
c				making a single line out of what would
c				otherwise be two fragments.
c
c	    96      XSTRT YSTRT	Signals that the current contour line has
c				closed on the starting point.
c
c	    97      0.0  0.0	Signals that the current contour line has
c				ended without closing on the starting point.
c				The previous values for XPT and YPT represent
c				the end point of this line.
c
c	    98	    ILEV ZLEV	Signals the end of contouring for the level
c				ILEV whose value is ZLEV.
c
c	    99	    NLEV 0.0	Signals the end of the contouring process.
c				NLEV is the number of levels that were
c				contoured.
c
c	Values for KODE other than the ones specified above may be used
c	by CON_DISP for user-defined purposes (for instance, initializing
c	the display, changing colors or line types, etc.).
c
c
c
c	Inputs:
c
c	    KODE	integer*4	Defined above.
c
c	    X		real*4		Defined above.
c	    Y
c
c
c
c	Outputs:
c
c	    None
c
c
c
c	History:
c
c	    01 May 88	Phil McDonald, NOAA/ERL/FSL	Operational version
c
c	    11 Jun 91	Phil McDonald, NOAA/ERL/FSL	Substantially revamped.
c
c
c*******************************************************************************
c*******************************************************************************
c
c
      logical	closed, reverse
c
c
c
      if (kode .eq. 0) then
          npt = npt + 1
c         call vector (x, y)
      else
          if ((kode .lt. 90) .or. (kode .gt. 99)) then
              user1   = x
              user2   = y
          else if (kode .eq. 90) then
              nlev    = x
          else if (kode .eq. 91) then
              ilev    = x
              zlev    = y
          else if (kode .eq. 92) then
              npt     = 1
              xstart  = x
              ystart  = y
              closed  = .false.
              reverse = .false.
c             call frstpt (x, y)
          else if (kode .eq. 95) then
              reverse = .true.
          else if (kode .eq. 96) then
              closed  = .true.
          else if (kode .eq. 97) then
              closed  = .false.
          else if (kode .eq. 98) then
              ilev    = x
              zlev    = y
          else if (kode .eq. 99) then
              nlev    = x
          end if
      end if
c
c
c
      return
      end