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