c
c*******************************************************************************
c*******************************************************************************
c
c
c
subroutine contour_rec (z, nx, ny, parms, link, con_disp)
c
c
real*4 z(*)
integer*4 nx, ny
real*4 parms(*)
byte link(2,*)
external con_disp
c
c
c*******************************************************************************
c*******************************************************************************
c
c Contour the data field Z using the rectangular subdivisions
c indicated by the 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 2 * 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_rec.inc'
byte ijdx(2,2), ijdy(2,2)
logical do_border, dead_end
data ijdx / 0, -1, 0, 1 /
data ijdy / -1, 0, 1, 0 /
if (nx .lt. 2) return
if (ny .lt. 2) return
c
c*******************************************************************************
c
c Initialize things
c
c*******************************************************************************
c
npt = nx * ny
do ipt = 1, npt
link(1,ipt) = 0
link(2,ipt) = 0
end do
do ipt = npt-nx+1, npt
link(2,ipt) = -1
end do
do ipt = nx, npt, nx
link(1,ipt) = -1
end do
do ileg = 1, 2
do idir = 1, 2
ijdpt(ileg,idir) = ijdx(ileg,idir) +
+ (ijdy(ileg,idir) * nx)
end do
end do
idp(1) = 1
idp(2) = nx
nxdim = nx
ip = 0
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. If these two values are
c equal, the data minimum and maximum will be
c used instead. This range will be divided into
c |NLEV| evenly spaced contour 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
nlev = -nlev
ip = ip + 1
zmin = parms(ip)
ip = ip + 1
zmax = parms(ip)
else if (nlev .eq. 0) then
ip = ip + 1
zint = parms(ip)
end if
if (zmax .eq. zmin) then
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
end if
if (zint .eq. 0.0) then
if ((nlev .eq. 1) .or. (zmax .le. zmin)) then
zint = 1.0
zmax = zmin
else
zint = (zmax - zmin) / float (nlev - 1)
end if
else
zmin = float (nint (zmin / zint)) * zint
zmax = float (nint (zmax / zint)) * zint
end if
if (zint .eq. 0.0) return
end if
c
c*******************************************************************************
c
c Signal the start of contouring.
c
c*******************************************************************************
c
x = nlev
y = 0.0
call con_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
end do
c
c-------------------------------------------------------------------------------
c
c Signal the start of this contour level
c
c-------------------------------------------------------------------------------
c
x = ilev
y = zlev
call con_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 = 0
11 if (ileg .ge. 2) go to 10
ileg = ileg + 1
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 con_find (ipt, jpt, z, zlev, xhit, yhit, nhit)
xstart = xhit
ystart = yhit
call con_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 con_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)
mhit = nhit
call con_find (ipt, jpt, z, zlev, xhit, yhit, nhit)
if (nhit .gt. mhit) call con_disp (0, xhit, yhit)
xprev = xhit
yprev = 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
itmplnk = link(ilstrt,ipstrt)
link(ilstrt,ipstrt) = 0
ipt = kpt
ileg = kleg
idir = kdir
call con_next (ipt, ileg, idir, z, zlev, link)
link(ilstrt,ipstrt) = itmplnk
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
link(ileg,ipt) = 1
call con_next (ipt, ileg, idir, z, zlev, link)
if (ipt .gt. 0) then
call con_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 con_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 con_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 con_disp (99, x, y)
c
c
c
return
end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
subroutine con_next (ipt, ileg, idir, z, zlev, link)
c
c
integer*4 ipt, ileg, idir
real*4 z(*), zlev
byte link(2,*)
c
c
c*******************************************************************************
c*******************************************************************************
c
c Test for an intersection between a contour line and a grid
c
c*******************************************************************************
c*******************************************************************************
c
c
include 'contour_rec.inc'
jpt = ipt + idp(ileg)
if (ipt .gt. 0) then
else
end if
if (z(ipt) .lt. z(jpt)) then
ilo = ipt
ihi = jpt
else
ilo = jpt
ihi = ipt
end if
kleg = 3 - ileg
jpt = 0
jlo = ilo + ijdpt(ileg,idir)
klo = ilo
if (jlo .lt. ilo) klo = jlo
if (klo .gt. 0) then
if (link (kleg,klo) .eq. 0) then
if (z(jlo) .ge. zlev) then
jpt = klo
go to 10
else
link(kleg,klo) = -1
end if
end if
end if
jhi = ihi + ijdpt(ileg,idir)
khi = ihi
if (jhi .lt. ihi) khi = jhi
if (khi .gt. 0) then
if (link(kleg,khi) .eq. 0) then
if (z(jhi) .lt. zlev) then
jpt = khi
go to 10
else
link(kleg,khi) = 1
end if
end if
end if
10 if (jpt .gt. 0) then
jleg = kleg
jdir = 1
if (jpt .ne. ipt) then
kpt = ipt + ijdpt(ileg,idir)
if (kpt .ne. jpt) jdir = 2
end if
ipt = jpt
ileg = jleg
idir = jdir
link(ileg,ipt) = 1
go to 99
end if
ipt = 0
if ((jlo .gt. 0) .and. (jhi .gt. 0)) then
if (link(kleg,khi) .ne. link(kleg,klo)) then
if ((z(jlo) .lt. zlev) .and. (z(jhi) .ge. zlev)) then
kpt = jhi
if (jlo .lt. jhi) kpt = jlo
if (link(ileg,kpt) .eq. 0) then
ipt = kpt
link(ileg,ipt) = 1
end if
end if
end if
end if
99 continue
if (ipt .gt. 0) then
else
end if
c
c
c
return
end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
subroutine con_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_rec.inc'
call con_i2xy (ipt, ix, iy)
call con_i2xy (jpt, jx, jy)
zz = zlev
dz = z(jpt) - z(ipt)
if (zz .eq. z(ipt)) zz = zz + (dz * 0.1)
if (zz .eq. z(jpt)) zz = zz - (dz * 0.1)
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
if (nhit .gt. 0) then
if ((xhit .eq. xprev) .and. (yhit .eq. yprev)) return
end if
xprev = xhit
yprev = yhit
nhit = nhit + 1
c
c
c
return
end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
subroutine con_i2xy (ipt, ix, iy)
c
c
integer*4 ipt, ix, iy
c
c
c*******************************************************************************
c*******************************************************************************
c
c Test for an intersection between a contour line and a grid
c
c*******************************************************************************
c*******************************************************************************
c
c
include 'contour_rec.inc'
iy = (ipt - 1) / nxdim
ix = ipt - (iy * nxdim)
iy = iy + 1
c
c
c
return
end
c
c*******************************************************************************
c*******************************************************************************
c
c
c
subroutine con_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